《结构分析中的有限元法》有限元习题参考答案_第1页
《结构分析中的有限元法》有限元习题参考答案_第2页
《结构分析中的有限元法》有限元习题参考答案_第3页
《结构分析中的有限元法》有限元习题参考答案_第4页
《结构分析中的有限元法》有限元习题参考答案_第5页
已阅读5页,还剩47页未读 继续免费阅读

下载本文档

版权说明:本文档由用户提供并上传,收益归属内容提供方,若内容存在侵权,请进行举报或认领

文档简介

1、本科有限元习题参考答案2015年3月10日作业1、简述力学课程中介绍的各种力学模型的简化条件、基本假设和适用范围(包括有拉压杆模型、弯曲梁模型、平面应力和平面应变模型、轴对称模型、板模型、壳模型等)力学模型简化条件基本假设适用范围拉压杆线弹性,无塑性变形拉压平面假设拉伸或压缩的轴向二力杆弯曲梁弯曲为主要变形的曲杆弯曲平面假设单向受力假设弯曲为主要变形的曲杆平面应力平面应变微体所受所有外力均在同一平面;构件内一点外侧变形均在同一平面理想弹性体只有平面应力分量在仅为x,y函数;只有平面应变存在仅为x,y函数;轴对称具有对称性对称的理想弹性体弹性体几何形状、约束情况、及所受外力都对称于某一轴板板厚远

2、小于其他两个方向的尺寸弹性体板件壳曲面薄板弹性体壳体2、 给出弹性力学问题中平衡方程、几何方程、物理方程的表达式及其意义。(1) 平衡方程:物理意义:应力分量与体力分量之间的关系。(2) 几何方程:物理意义:应变分量与位移分量之间的关系。(3) 物理方程:物理意义:应变分量与应力分量之间的关系。3、 简述最小势能原理的主要内容和主要公式。根据虚功原理得到:,由则其中,即为系统的总势能,它是弹性体变形势能和外力势能之和。上面变分为零式表明:在所有区域内满足几何关系,在边界上满足给定位移条件的可能位移中,真实位移使系统的总势能取驻值(可证明此驻值为最小值)。4、 “有限元法”都有哪些名称(包括中文

3、和英文)?有限元法,也叫有限单元法、有限元素法、有限元分析;fem(finite element method)、fea(finite element analysis)5、 简述有限元法的发展和现状。近几十年,伴随着计算机科学和技术的快速发展,有限元法作为工程分析的有效方法在理论、方法的研究、计算机程序的开发以及应用领域的开拓者方面均取得了根本性的发展。(1)单元的类型和形式为了扩大有限元法的应用领域,新的单元类型和形式不断涌现(等参元,梁板壳,复合材料)(2)有限元法的理论基础和离散格式将hellinger-reissner、huwashizu(多场变量变分原理)应用于有限元分析,发展了混

4、合模型、杂交型的有限元表达格式,应研究了各自的收敛条件;将加权余量法用于建立有限元的表达格式;进一步研究发展有限元解的后验误差估计和应力磨平方法。(3)有限元方程的解法(大型复杂工程结构问题静态, 特征值, 瞬态等)(4)有限元法的计算机软件(专用软件, 通用软件)6、 结构力学分析问题中的三种“非线性”都包含哪些,并解释其含义。非线性可以是由材料性质、变形状态和边界接触条件引起的,分别称为材料、几何、边界非线性。材料非线性就是材料的本构关系不是线性的。几何非线性时结构在载荷作用过程中产生大的位移和转动,如板壳结构的大挠度。边界非线性是指高挠度部件或由多个部件组成的结构组合件,渐进位移将会增大

5、部件自身或是部件之间产生接触的可靠性,以此特征的特定类型几何非线性为边界条件或者接触非线性。7、 简述有限元法的未来。以有限元法为代表的计算力学提出一系列新的课题:(1)为了真实地模拟新材料和新结构的行为,需要发展新的材料本构模型和单元型式。例如对于特种合金、复合材料、陶瓷材料、机敏材料、智能材料、生物材料以及纳米材料等。(2)为了分析和模拟各种类型和形式的结构在复杂载荷工况和环境作用下的全寿命过程的响应,需要发展新的数值分析方案。如:多重非线性(材料、几何、边界等)相耦合的分析方法;多场(结构、流体、热、电、化学)耦合作用的分析方法;跨时间空间多尺度;非确定性(随机模糊)的分析方法;自适应的

6、分析方法。(3)有限元软件和cadcam等软件系统共同集成完整的虚拟产品开发(vpd)系统。这个系统强烈影响着未来工程系统的设计、制造、和运行,主要体现在: 它能提供对所设计的工程系统从加工制造到运行,直至失效和破坏的全寿命过程的更深入认识,从而能更好地识别它的属性和特征。 它能够鉴定和评估所设计对象的性能和质量,并允许以最低的费用在设计过程中就对所设计的对象进行修改和优化。 它能显著地缩短工程对象设计和投产的周期,降低生产成本,提高市场竞争力。2015年3月17日作业1、简述有限元法的基本思想,并结合简单结构来说明。(1) 有限元法,也叫有限单元法,它的基本思想是将一个结构或连续体的求解域离

7、散为若干个子域(单元),并通过它们边界上的结点相互联结成为组合体。(2) 有限元法用每一个单元内所假设的近似函数来分片地表示全求解域内待求的未知场变量。而每个单元内的近似函数由未知函数或其导数在单元各个结点上的数值和与其对应的插值函数来表示。由于在联结相邻单元的结点上,场函数应具有相同的数值,因而将它们用作数值求解的基本未知量。这样一来,求解原来待求场函数的无穷自由度问题转换为求解场函数结点值的有限自由度问题。(3) 有限元法是通过和原问题数学模型(基本方程、边界条件)等效的变分原理或加权余量法,建立求解基本未知量(场函数的结点值)的代数方程组或微分方程组。此方程组称为有限元求解方程,并表示成

8、规范的矩阵形式。接着用数值方法求解此方程,从而得到问题的解答。2、简述结构离散(或有限元建模)的内容和要求。有限元建模的内容:1)网格划分-即把结构按一定规则分割成有限单元2)边界处理-即把作用于结构边界上约束和载荷处理为结点约束和结点载荷有限元建模的要求:1)离散结构必须与原始结构保形-单元的几何特性2)一个单元内的物理特性必须相同-单元的物理特性3、简述结点力和结点载荷的差别。结点力:单元与单元间通过结点的相互作用力。结点载荷:作用于结点上的外载荷。4、列表给出有限元几类基本单元的图形、结点数、结点自由度数和单元总自由度数(包括杆单元、梁单元、平面三角形单元、平面四边形单元、轴对称问题三角

9、形单元、四边形壳单元、四面体单元)。单元类型单元图形结点数结点自由度杆单元22梁单元23平面单元32平面四边形42轴对称问题32板壳单元46四面体单元435、写出适用于插值下表的基函数,并给出拉格朗日插值多项式。基函数:,根据拉格朗日插值多项式:。将带入:2015年3月24日作业1、为了保证有限单元法解答的收敛性,位移函数应满足哪些条件?完备协调元、非协调元和完备元分别是什么意思?为了保证有限单元法解答的收敛性,位移函数应满足:1) 位移函数必须包括单元的刚性位移(即常量项);2) 位移函数必须包括常量应变(即线性项);3) 位移函数在单元内部必须连续(连续性条件);4) 位移函数应使得相邻单

10、元间的位移协调(协调性条件)注:上述四个条件称为有限元解收敛于真实解的充分条件;前三个条件称为必要条件。满足四个条件的位移函数构成的单元称为完备协调元;满足前三个条件的单元称为非协调元;满足前两个条件的单元称为完备元。2、如下图所示4结点平面应力单元,结点1结点4对应的结点坐标分别为(0,0), (0,1), (2,0),(2,1),结点1结点4对应的结点位移分别为(u1,v1), (u2,v2), (u3,v3), (u4,v4),试基于拉格朗日插值基函数构造如下单元的位移函数。解:根据拉格朗日插值基函数:所以:同求:。3、说明有限元方法解误差的主要来源?答:影响有限元解的误差:1) 离散误

11、差。边界上以直线代曲线导致离散化模型与实际物体的差异;2) 位移函数误差。一般情况下单元位移函数不可能与实际单元的位移场一致;3) 计算机计算误差。计算机字长的限制、相差悬殊的数值加减运算。4、说明用有限单元法解题的主要步骤。答:研究问题的力学建模;结构离散;单元分析;整体分析与求解;结果分析及后处理。5、推导基于变分原理的总势能泛函极值条件。解:有积分形式确立的标量泛函有其中是未知函数,和是特定的算子,是求解域,是的边界。称为未知函数的泛函,随函数的变化而变化。连续介质问题的解使泛函对于微小的变化取驻值,即泛函的“变分”等于零,此为变分法。将虚功原理用于弹性变形时,总功w要包括外力功(t)和

12、内力功(u)两部分,即:w = t - u;内力功(-u)前面有一负号,是由于弹性体在变形过程中,内力是克服变形而产生的,所有内力的方向总是与变形的方向相反,所以内力功取负值。根据虚功原理,总功等于零得:t- u = 0,即外力虚功t = 内力虚功u弹性力学中的虚功原理可表达为:在外力作用下处于平衡状态的弹性体,如果发生了虚位移,那么所有的外力在虚位移上的虚功(外力功)等于整个弹性体内应力在虚应变上的虚功(内力功)。根据虚功原理得到其中的即为总势能泛函。由上面变分为零式表明:在所有区域内满足几何关系,在边界上满足给定位移条件的可能位移中,真实位移使系统的总势能取驻值(可证明此驻值为最小值)。此

13、即总势能泛函的极值条件。2015年3月31日作业1、给出利用最小势能原理建立单元有限元静力平衡方程的一般推导过程。解:由单元位移函数:式中:为插值函数(或称形函数),得到单元内的应变和应力分别为 其中为弹性矩阵,它完全取决于弹性常数和。将位移、应力和应变代入势能泛函有根据最小势能原理,势能泛函取驻值的必要条件:式中称为单元刚度矩阵,称为体积力等效结点力,称为面力等效结点力,则上式可写成单元平衡方程:此式即单元有限元静力平衡方程。2、利用最小势能原理推导结构(假设包含n个单元)有限元静力平衡方程。解:利用弹性体所有单元的总势能最小,将所有单元进行组装后即可形成整个结构的有限元静力平衡方程。为此先

14、需要建立单元结点位移向量和总体结点位移向量之间的对应关系,即式中:表示第个单元的结点位移向量;是总体结点位移向量;矩阵是第个单元的结点位移提取矩阵。对于杆单元可以表示为 利用提取矩阵对单元矩阵进行合同变换与累加后,可得到总体的刚度矩阵、质量矩阵和载荷向量,即 根据最小势能变分原理即可得到,此式即为总体结构的有限元静力平衡方程。2015年4月7日作业1、推导拉(压)杆单元形状函数,利用最小势能原理推导拉(压)杆单元在局部坐标系中的单元刚度方程。解:(1)如图1表示某一杆单元,现约定附属于该单元的局部坐标系(单元坐标系)为,点为原点,轴沿着杆轴线,其正方向为由指向,其余各轴按右手螺旋规则确定。设,

15、为杆元结点位移分量,杆单元结点力分量,一律规定和坐标轴正向一致时为正。设杆的长度为,弹性模量为,横截面积为。图1 杆单元结点位移结点力分量对于铰接杆单元,在小变形假设的前提下,与杆垂直方向的位移并不使杆产生应变和应力。于是,对每一个结点只需考虑一个结点位移和结点力,即如图2所示图2 二力杆单元二力杆单元的位移函数为。式中,是两个待定常数,可由,两结点的位移唯一确定。当;时,将其代入位移函数有:,从而可以得到将,的值代入位移函数得或写成则有式中,称为点、点的形状函数,成为形函数矩阵。(2)根据位移函数与应变的定义得或写成,其中称为应变矩阵。根据应力与应变的关系,将应变代入有,其中称为应力矩阵,对

16、于拉(压)杆有。根据杆单元在局部坐标系单元刚度矩阵的表达式代入矩阵与得到拉(压)杆单元在局部坐标系中的单元刚度方程为2、给出将局部坐标系中的拉(压)杆单元刚度矩阵通过坐标变换得到总体坐标系中的单元刚度矩阵推导过程。解:设为总体坐标系,为局部坐标系,如图1所示。图1 平面杆单元总体坐标系位移在局部轴方向分量规定由总体坐标系平面轴到局部坐标系轴的夹角逆时针为正。杆单元总体坐标系下的结点位移分量用表示,局部坐标下的位移分量用表示。则平面杆单元结点在总体坐标系和局部坐标系下的位移分量关系有同理对于结点有用矩阵表示为对两结点杆单元,当用总体坐标系位移表示局部坐标系中位移时有转换关系式中矩阵称为坐标变换矩

17、阵。对于图1中的杆单元就可以很容易将局部坐标系的刚度矩阵转换为总体坐标系的刚度矩阵为3、计算下图所示三杆平面桁架结构的总体坐标系下各单元的单元刚度矩阵。杆单元材料弹性模量e=2.0×1011(n/m2),截面积a=1.0×10(-4)(m2)。注:物理量单位均采用国际单位。解:总体坐标系下的单元刚度矩阵为对杆有: 对杆有: 对杆有: 4、推导平面梁弯曲的挠曲线控制方程。解:如下图为梁单元,结点为和因弯曲而引起的轴向位移为,其中是所讨论点离中面的距离,应变。又根据平衡关系有,其中是切面惯性矩,。根据胡克定律,有又根据剪力与弯矩的关系有现在讨论的梁,为常数,故有,若梁上无分布剪

18、力,则有,则可以判断是的三次函数。设则梁的转角将,结点位移代入上式有用矩阵表示为从而得单元的形状函数矩阵为其中设,则上式可写成以上四个函数即是二结点梁单元的形状函数由于每结点有两个位移参数,则每结点有两个形状函数。是指结点零阶导数(即点位移)对应的形函数, 指结点一阶导数(即结点转角)对应的形函数,位移函数用内插多项式表示为2015年4月14日作业1、推导平面直梁单元形状函数,利用最小势能原理推导平面直梁单元在局部坐标系中的单元刚度方程。 解:先将轴力与剪力、弯矩分开考虑。对于平面直梁单元的两个节点1,2各有一个挠度和转角,故可设其挠度为:其中,为待定常数,为单元形函数,;利用hermit插值

19、或书中解方程的方法,可求得: ,对于轴力引起的位移,设:进而,位移函数改写为位移函数求得后,可得到应变和应力的表达式:若忽略剪切影响,是拉力引起的应变,是弯曲引起的应变对应的应力为下面求梁单元的总势能。 梁应变能为:梁上的结点力,并有分布力作用。故外力势为:总势能式中是分布载荷的等效结点载荷,其中由最小势能原理,由的任意性,知其中刚度矩阵将b带入上式,得到2、计算下表所示直梁单元承受的4种典型载荷分布形式的等效结点载荷(要求计算过程!)解:等效节点载荷的基本原理是功的互等原理,即对任意虚位移,等效节点力做的功与实际力做的功相等。对于分布载荷,可推得等效节点载荷的计算公式(书中p28)为: (1

20、)其中为单元形函数向量。对于梁单元: ,设等效节点载荷为(1) 当梁发生虚位移为形函数时,由形函数定义知,等效节点载荷中只有做功,由功的互等原理,有:故同理求得:对任意虚位移,等效节点力做的功与实际力做的功相等。而对于梁单元,所有容许的位移都可表示为单元形函数的线性组合。所以,当每个形函数作为虚位移,等效节点力与实际力做功相等时,对任意容许位移,等效节点力与实际力做功也必相等。满足功的互等。(2) 将带入公式(1),求得等效节点载荷为:(3) 将带入公式(1),求得等效节点载荷为:(4) 将带入公式(1),求得等效节点载荷为:3、计算如下所示受一个集中力和一个集中力矩的直梁单元的等效结点载荷。

21、解法一:仿照2中(1)的方法,设等效节点载荷列向量为。当梁发生虚位移为时,由形函数定义知,等效节点载荷中只有做功,由功的互等原理,有:同样地,解法二:利用狄拉克函数的性质,可设,带入第2题公式(1)中,得:进而可得出与解法一相同的结果。2015年4月21日作业1、请给出由单元刚度矩阵组装形成总体刚度矩阵的基本原理和过程。(补充公式推导)基本原理:结构总体刚度矩阵的元素是由单元刚度矩阵的元素组成的,只要确定了单元刚度矩阵各元素在结构总体刚度矩阵中的位置,就可以由单元刚度矩阵直接集成结构总体刚度矩阵。基本思想是:把单元杆端位移分量(局部码)所对应的结构结点位移向量的序号(整体位移码)组成一向量,它

22、成为单元的定位向量。利用单元定位向量可以完全确定单元刚度矩阵的每个元素在结构(原始)刚度矩阵中的行码和列码。过程:先求出单元在整体坐标系中的刚度矩阵,然后将单元 的定位向量分别写在单元刚度矩阵的上方和右侧(或左侧)。这样, 的每一行或每一列就与单元定位向量的一个分量相对应,这个分量即为 中相应的行和列在结构刚度矩阵 中的行码或列码。于是,按照由单元定位向量中分量的行码和列码,就能够将单元刚度矩阵的元素正确地叠加到结构刚度矩阵中去。2、计算下图所示三杆平面桁架结构的总体刚度矩阵。杆单元材料弹性模量,截面积。 注:物理量单位均采用国际单位。解:先计算各单元相对总体坐标系的刚度矩阵。各单元局部坐标系

23、的选取与到的角度如下表:单元编号局部坐标系原点cossin20102 013-arctan(3/4)0.8-0.6由书中p20公式(3.31)得到各单元的刚度矩阵为: ()组装得到总刚:3、从半带宽最小的角度出发,如下图所示平面桁架或刚架的结点编号哪一种最好,其“半带宽b”是多少?解:对于桁架计算半带宽:(a) 单元节点号之差最大的为3(如节点1和4),节点自由度数为2半带宽: (b) 单元节点号之差最大的为9(如节点1和10),节点自由度数为2半带宽:(c) 单元节点号之差最大的为5(如节点3和8),节点自由度数为2半带宽:对于桁架计算半带宽:(d) 单元节点号之差最大的为3(如节点1和4)

24、,节点自由度数为3半带宽: (e) 单元节点号之差最大的为9(如节点1和10),节点自由度数为3半带宽:(f) 单元节点号之差最大的为5(如节点3和8),节点自由度数为3半带宽:由上可知,(a)编号最好2015年4月28日作业1、计算下图所示三杆平面桁架结构的结点位移、各杆轴力、轴向正应力及支座约束反力。1结点处受x和y方向的载荷大小为。杆单元材料弹性模量,截面积。解:先计算各单元相对总体坐标系的刚度矩阵。各单元局部坐标系的选取与到的角度如下表:单元编号局部坐标系原点cossin20102 013-arctan(3/4)0.8-0.6由书中p20公式(3.31)得到各单元的刚度矩阵为: ()将

25、单元刚度矩阵组成总刚,各节点x、y方向上对应的位移依次设为,引入边界条件,用删行删列法处理,删去第2至第5行与第2至第5列,总刚度方程变为: 下面求各杆轴力、轴向正应力。对于号杆,局部坐标系中,带入局部单元刚度矩阵,轴向应力轴力为:同理可求解号杆:轴向应力:轴力:对于号杆,应先将相关位移分量转到局部坐标系中即( 还可以比较简单地通过1,3节点位移在杆方向上投影得到)所以轴向应力:轴力:下面求解约束力 总刚为:故各节点力为:所以支反力为 2、求下图所示平面刚架2结点处结点位移和梁内力。在2结点处受一集中弯矩作用。其中梁截面积, 惯性矩,弹性模量,尺寸。解:先写出两个单元的刚度矩阵单元,局部坐标系

26、取2为原点,角度局部坐标中单刚(p26,(3.74):转换矩阵(p27):总体坐标系下的单刚为:单元,局部坐标系取2为原点,角度同理可得: 转换矩阵:总体坐标系下的单刚为:组装总刚,引入边界条件,求解得到: 解得: 下面求梁内力对于单元,局部坐标系下的位移为: 局部坐标系下,梁两端受力为:可由这些梁端力求出梁的剪力,轴力与弯矩分布。即轴力30n(拉),剪力56.625n,弯矩 (设以节点2位原点)对于单元,同理,可求得:即轴力19.5n(压),剪力62.25n,弯矩(设以节点2位原点)(注:1本题中,带撇号()的量表局部坐标系中的分量,不带撇号的表整体坐标系中的量;2 求内力的方法还有一种,通

27、过先求整体坐标系中单元梁端力分量,再左乘转换矩阵得出局部坐标中的力分量,即 )3、求下图所示平面刚架2结点和3结点位移。在2号单元上作用均布力。其中梁截面积,惯性矩,弹性模量,尺寸。解:先写出各个单元的刚度矩阵。单元,局部坐标系取1为原点,角度,整体坐标系下刚度矩阵为:单元,局部坐标系取2为原点,角度,整体坐标系下刚度矩阵为:单元,局部坐标系取3为原点,角度,整体坐标系下刚度矩阵为:下面列写刚度方程:载荷列向量为: 位移列向量为:将分布力等效为节点力得到 引入边界条件,删去总刚第一至第三行和列,以及第十至十二行和列,得到有限元方程: 解得2、3节点的位移如下: 2015年5月5日作业1. 如图

28、1和图2所示三角形平面单元,写出其形状函数,。图1 三角形平面应力单元 图2 三角形平面应变单元解:根据公式:则形函数:图1形函数:图2形函数:2. 证明平面问题的常应变三角形单元的形状函数ni,nj,nk 所满足的两个性质:(1)“本点为1,它点为0”;(2) 。证明:(1)对于常应变三角形单元,由定义,i节点的形函数可表示为:,对本点:(本点为1)对他点:(它点为0) (它点为0)对于j、k节点同理可证。(2)证一: 其中:故证二:由常应变三角形单元的位移函数的原始表达式: 知,此位移函数可以表达刚体位移,而位移函数又可以表达为形函数的线性组合,即:现设单元有刚体位移,则三个节点的位移均为

29、,进而:,即3列出图3所示各常应变三角形单元的等效结点载荷列向量。解:以下节点载荷向量的写法都是按i、j、m的顺序,等效面力计算公式:(a)此题与梁中节点力的简化大致相同,具体过程如下: 为方便积分,先以j为原点,jm为坐标轴建立局部坐标系s。因只需考虑jm边的位移,则在局部坐标系中常应变单元各点的形函数如下: ,(由形函数定义,jm上)设单元发生虚位移(整体坐标系中),则由功的互等得:设单元发生虚位移(整体坐标系中),则由功的互等得:同理,得:,故等效节点载荷(b)载荷作用在边上,以压力为正,另这条边长为,与轴夹角为。测压在轴轴方向的分量为和为:。在局部坐标中,由积分公式知,只需要考虑jm的

30、位移,故形函数可写作:(由形函数定义,jm上)则单元的等效节点载荷为:因此。(c)为方便积分,先以i为原点,ij为坐标轴建立局部坐标系s。因只需考虑ij边的位移,则在局部坐标系中常应变单元各点的形函数如下: ,作用在ij边的载荷为: i节点: ,y方向无载荷分量: j节点:,y方向无载荷分量:m节点:故载荷向量为:(d)为方便积分,先以j为原点,jm为坐标轴建立局部坐标系s。则在局部坐标系中常应变单元各点的形函数如下: ,作用在ij边的载荷为:j节点: ,y方向无载荷分量: m节点:,y方向无载荷分量:i节点:故载荷向量为:4 求下图平面应力问题(仅一个平面三角形单元)的节点处竖向位移,其中,弹性模量 泊松比厚度。解:这时边界上面积力写作局部坐标的函数单元等效节点载荷为代入数据后

温馨提示

  • 1. 本站所有资源如无特殊说明,都需要本地电脑安装OFFICE2007和PDF阅读器。图纸软件为CAD,CAXA,PROE,UG,SolidWorks等.压缩文件请下载最新的WinRAR软件解压。
  • 2. 本站的文档不包含任何第三方提供的附件图纸等,如果需要附件,请联系上传者。文件的所有权益归上传用户所有。
  • 3. 本站RAR压缩包中若带图纸,网页内容里面会有图纸预览,若没有图纸预览就没有图纸。
  • 4. 未经权益所有人同意不得将文件中的内容挪作商业或盈利用途。
  • 5. 人人文库网仅提供信息存储空间,仅对用户上传内容的表现方式做保护处理,对用户上传分享的文档内容本身不做任何修改或编辑,并不能对任何下载内容负责。
  • 6. 下载文件中如有侵权或不适当内容,请与我们联系,我们立即纠正。
  • 7. 本站不保证下载资源的准确性、安全性和完整性, 同时也不承担用户因使用这些下载资源对自己和他人造成任何形式的伤害或损失。

评论

0/150

提交评论