有限元习题与答案

习题2.1 解释如下的概念：应力、应变，几何方程、物理方程、虚位移原理。解 应力是某截面上的应力在该处的集度。 应变是指单元体在某一个方向上有一个U的伸长量，其相对变化量就是应变。表示在x轴的方向上的正应变，其包括正应变和剪应变。 几何方程是表示弹性体内节点的应变分量与位移分量之间的关系，其完整表示如下： 物理方程：表示应力和应变关系的方程某一点应力分量与应变分量之间的关系如下：虚位移原理：在弹性有一虚位移情况下，由于作用在每个质点上的力系，在相应的虚位移上虚功总和为零，即为：若弹性体在已知的面力和体力的作用下处于平衡状态，那么使弹性体产生虚位移，所有作用在弹性体上的体力在虚位移上所做的工就等于弹性体所具有的虚位能。2.2说明弹性体力学中的几个基本假设。 连续性假设：就是假定整个物体的体积都被组成该物体的介质所填满，不存在任何间隙。 完全弹性假设：就是假定物体服从虎克定律。 各向同性假设：就是假定整个物体是由同意材料组成的。 小变形和小位移假设：就是指物体各点的位移都远远小于物体原来的尺寸，并且其应变和转角都小于1。2.3简述线应变与剪应变的几何含义。线应变：应变和刚体转动与位移导数的关系，剪应变表示单元体棱边之间夹角的变化。2.4 推到平面应变平衡微分方程。解：对于单元体而言其平衡方程：在平面中有 代入上式的 2.5 如题图2.1所示，被三个表面隔离出来平面应力状态中的一点，求和的值。解：x方向上：联立二式得：2.6相对于xyz坐标系，一点的应力如下某表面的外法线方向余弦值为，求该表面的法相和切向应力。解：该平面的正应力全应力该平面的切应力2.7一点的应力如下MP求主应力和每一个主应力方向的方向余弦；球该店的最大剪应力。解：设主平面方向余弦为，由题知将代入得即 ，。最大剪应力（1）当时代入式（2.21）（2）当时代入式（2.21）且 2.8已知一点P的位移场为，求该点p(1,0,2)的应变分量。解：p点沿坐标方向的位移分量为u,v,w点p(1,0,2)处线应变为，剪应变为，2.9一具有平面应力场的物体，材料参数为E、v。有如下位移场 其中，a、b、c、d是常量。求讨论位移场的相容性解： 因为 所以满足相容性条件有广义胡克定律得又则2.10一具有平面应力场的物体，材料性质是E=210GPa,v=0.3.并且有如下位移场当x=0.050m,y=0.020m时，求物体的应力和应变。位移场是否相容？解：由广义胡克定律，满足相容性条件2.11对于一个没有任何体积力的圆盘，处于平面应力状态。其中 a, b, c, d, e, f, g, h是常量。为了使应力满足平衡方程和相容方程，这些常量的约束条件是什么？解：由题意得：，代入平衡方程根据广义胡克定律： 代入相容方程 （2）代入（1）得其中2.13 根据弹性力学平面问题的几何方程，证明应变分量满足下列方程，并解释该方程的意义。证明：弹性力学平面问题的几何方程为： ， ， ，将方程，分别对y和x求二阶偏导并相加得：等式右端项，该方程为相容方程中的第一式，其意义为弹性体内任一点都有确定的位移，且同一点不可能有连个不同的位移，应变分量应满足相容方程，否则，变形后的微元体之间有可能出现开裂与重叠。2.14 假设Airy应力函数为，其中为常数，求，并求这些变量间的约束关系。解：由，对该应力函数求偏导得;对以上两式的偏导可求得：考虑相容性条件,将上式代入可得各常量间的关系如下：2.15 对给定的应力矩阵，求最大Tresca和Von.Mises应力。将Von Mises应力和Tresca应力 20 10 10进行比较，= 10 20 10 Mpa。10 10 20 z xy xz解：由Tresca准则：= y yz 故有s=20Mpa，max=s/2=10Mpa z1=（x+y）/2=30Mpa 2=10Mpa由Von Mises准则：2s2=6（xy2+yz2+yz2）解得s=30Mpa 30 -15 202.16 一点出的应力状态由应力矩阵给出，即= -15 -25 10 Mpa，若E=70Gpa， 20 10 40 =0.33，求单位体积的应变能。解：单位体积应变能：=1/2Ex2+y2+z2-2u(xy+yz+zz)+2(1+u)(xy2 +xz2+yz2)u=（E-2）/2 =0.33带入可得：=420.75J3.11 如图3.11所示的平面三角形单元，厚度t=1cm，弹性模量E=2.0*105mpa，泊松比=0.3，试求插值函数矩阵N，应变矩阵B，应力矩阵S，单元刚度矩阵Ke。解：此三角形单元可得：2=（10-2）*4=32，故有a1=1/32*（8u1-5u2-16u3）a2=1/32*（4u1-4u2）a3=1/32*（-8u1+8u3）a4=1/32*（56v1-8v2-16v3）a5=1/32*（-4v1+4v2）a6=1/32*（-8v1+8v3）而b1=y2-y3=-4 b1=x2-x3=-8 b1=y3-y1=4 b1=x3-x1=0 b1=y1-y2=0 b1=x1-x2=8 b1 0 b2 0 b3 0 -4 0 4 0 0B=1/2* 0 c1 0 c2 0 c3 =1/32* 0 -8 0 0 8 c1 b1 c2 b2 c3 b3 -8 4 0 8 0 1 0 1 0.3 0D=E/(1-2)* 1 0 =E/0.91* 0.3 1 0 0 0 (1-)/2 0 0 0.35 1 0.3 0 -0.125 0 0.125 0 0S=D*B=E/0.91* 0.3 1 0 * 0 -0.25 0 0 0.25 0 0 0.35 -0.25 0.125 0 0.25 0 1.4 0 -1.4 -0.7 0 0.7 0 4 -0.6 -4 0 0K=BT*D*B*t*=E/36.4* -1.4 -0.6 2.4 1.3 0.6 0.7 -0.7 -4 1.3 -0.6 -1 0.35 0 0 0.6 -1 -0.6 0 0.7 0 0.7 -0.35 0 0 1 0 0 0.6 -1 -0.6 0 0.35 0.7 0 -0.7 -0.35 0 0.7 1.4 0 -1.4 -0.7K=BT*D*B*t*=E/36.4* 0.6 0 0 4 -0.6 -4 1 -0.7 -1.4 -0.6 2.4 1.3 0.6 -0.35 -1.4 -4 1.3 3.53.12 求下图中所示的三角形的单元插值函数矩阵及应变矩阵，u1=2.0mm，v1=1.2mm，u2=2.4mm，v2=1.2mm，u3=2.1mm，v3=1.4mm，求单元内的应变和应力，求出主应力及方向。若在单元jm边作用有线性分布面载荷（x轴），求结点的的载荷分量。解：如图2=64/3，解得以下参数：a1=19 a2=-2 a3=6； b1=-3 b2=4 b3=-1；c1=-1 c2=-3 c3=4；N1=64/3*(19-3x-y) N2=64/3*(-2-3x-3y)N3=64/3*(6-x+4y)故N= Ni 0 Nj 0 Nm 0 0 Ni 0 Nj 0 Nm 1 0 1 0 1 0 = 0 1 0 1 0 1 bi 0 bj 0 bm 0B=1/2* 0 ci 0 cj 0 cm ci bi cj bj cm bm -3 0 4 0 -1 0=64/3* 0 -1 0 -3 0 4 -1 -3 -3 4 4 -1 1 0D=E/(1-2)* 1 0 0 0 (1-)/2 1 0 -3 0 4 0 -1 0单元应力矩阵S=D*B= E/13(1-2)* 1 0 * 0 -1 0 -3 0 4 0 0 (1-)/2 -1 -3 -3 4 4 -1 2 1.1 -3 -u 4 3u -1 4u 2.4单元应力=S*q= E/13(1-2)* -3u -1 4u -3 -u 4 * 1.2 (u-1)/2 (3u-3)/2 (3u-3)/2 2-2u 2-2u (u-1)/2 2.4 1.43.13 解：二维单元在x,y坐标平面内平移到不同位置，单元刚度矩阵相同，在平面矩阵180时变化，单元作上述变化时，应力矩阵不变化。（0,1）（2,1）3.14（2,0）（0,0）yx解：令，而，单元单元： 由和扩充KZ（总刚度阵）而，其中，化简得：则，3.15如图所示有限元网格，单元厚度，弹性模量，泊松比。回答下述问题：（1）结点如何编号才能使结构刚度矩阵带宽最小？（2）如何设置位移边界条件才能约束结构的刚体移动？ （3）形成单元刚度矩阵并集成结构刚度矩阵。（4）如果施加一定载荷，拟定求解步骤。 (1) (2) （3）解：1、节点编号如图(2)所示；2、如图(3)设置位移边界条件才能约束结构的刚体移动；3、如图(2)所示各节点的坐标为(以m为单位)：1(0,0)，2(0.08,0)，3(0,0.04)，4(0.08,0.04 )，5(0,0.08)，6(0.08,0.08)，7(0,0.12)，8(0.08,0.12)解：单元号123456相邻结点134557225466343678对于单元号1：；对于单元号2：；对于单元号3：；对于单元号4：；对于单元号5：；对于单元号6：；平面三角形单元的面积均为 弹性矩阵均为 应变矩阵 应力矩阵 单元刚度矩阵 结构刚度矩阵为： 若施加一定载荷，求解步骤为：1、对单元编号，并列出各单元三个结点的结点号；2、计算外载荷的等效结点力，列出结构结点载荷列阵；3、计算单元刚度矩阵，组集结构整体刚度矩阵4、引入边界条件，即根据约束情况修正结构有限元方程，特别是消除整体刚度矩阵的奇异性，得到考虑约束条件的可解的有限元方程。5、利用线性方程组的数值解法，对结构的有限元方程进行求解，得到所有各结点的位移向量。最后根据需要求解单元应力。3.16一长方形薄板如图所示。其两端受均匀拉伸。板长12cm,宽4cm，厚1cm。材料，泊松比。均匀拉力。使用有限元法求解板的内应力，并和精确解比较（提示：可利用结构对称性，并用2个三角形单元对结构进行离散）。解：解：结点编号 12 34单元号12X坐标 012 012相邻结点13Y坐标 00 442234平面三角形单元的面积均为 应力矩阵为：单元1的应变距阵为：单元1的单元刚度矩阵为：单元2的应变距阵为：单元2的单元刚度矩阵为：总刚度矩阵为：位移分量为：载荷列阵为：因为 可以得单元1的单元应力： 单元2的单元应力： 长方形薄板内应力的精确解为：拉应力，用有限元法求解出的结果与精确解大致相等。3.17 验证三角形单元的位移差值函数满足及。解：平面三角形形函数为：，其中，分别是行列式2A中的第一行，第二行和第三行各元素的代数余子式。行列式中，任一行的元素与其相应的代数余子式的乘积之和等于行列式的值，而任一行的元素与其它行对应元素的代数余子式乘积之和为零，故有：当，同时有，同理也有：，即。3.18 推导如图所示的9节点矩形单元的形函数。解:三维杆单元的形状函数， 在局部坐标系中令节点1,5,2所对应的带入式得到节点1,5,2仅在x方向上的形函数： 同理可得： 由，即节点2,6,3，可得到沿着全局坐标系y轴的形状函数(通过变量轮换)，节点1的形函数即x，y方向的乘积：由此可得：同理可整理得：，3.19 如图所示为一个桁架单元，端点力为U1，U2，端点位移为u1，u2，设内部任一点的轴向位移u是坐标x的线性函数：推导其形函数矩阵N。解：轴向位移u是坐标x的线性函数，写成向量形式为，设两个节点的坐标为，代入向量形式的位移函数解得：则由位移函数可得形函数为：4.1 答：轴对称三角形环单元不是常应变单元，如果弹性体的几何形状、约束条件及载荷都对称于某一轴，则所有的位移应变及应力也是对称于此轴，这样问题称为轴对称。轴对称三角形环单元与平面常应变单元是不同的，轴对称三角形环单元的应变不是常数矩阵，其应变矩阵B=B B B,其中B=，（i,j,m）。应变分量，都是常量，但环向应变不是常量，它与，中的r和z有关。4.2 答：轴对称问题中，刚度自由度：环向位移，径向位移，轴向位移。以三角环单元平均半径、平均高度进行计算的单元刚度矩阵，配合以精确积分所得的等效结点载荷矩阵，计算的结果还是不错的！4.3 轴对称问题的两个单元a和b，设材料的弹性模量为E，泊松比为 = 0.15，试手算这两个单元的刚度矩阵。 解：对于单元，由题可知：单元a的截面面积为单元a的刚度矩阵写成分块矩阵形式为：其中子矩阵可写为：所以的刚度矩阵为对于单元，由题可知单元的截面面积为单元的刚度矩阵写成分块矩阵形式为：其中子矩阵可写为：所以单元的刚度矩阵为5.1 答：杆件受到纵向（平行于杆轴）载荷的作用，这样杆件的拉压问题；杆件受到横向（垂直于杆轴）载荷的作用，这是梁的弯曲问题。杆件受到力相似到薄板就有，薄板受到纵向载荷的作用，这是平面应力问题；薄板受到横向载荷的作用，这是薄板的弯曲问题。薄板的弯曲可以认为是梁弯曲的推广，是双向的弯曲问题，中面法线在变形后保持不伸缩，并且成为弹性曲面的法线，中面在变形后，其线段和面积的投影形状保持不变（小挠度薄板）。已知中面的挠度，而纵向位移、，主要应力分量，。某一点的位移：，。某一点的应力：,弹性曲面微分方程,其中板的抗挠刚度。5.2 答：矩形薄板单元：薄板单元位移函数并不满足连续性或相容性要求，采用这种位移函数的单元是非协调单元，这种四节点矩形弯曲单元变形后，其挠度面在单元间虽然互相连续，但其法向导数并不连续，单元间在变形后是不连续光滑（有棱）的，当单元逐渐取小的时候，还能够收敛于精确解。三角形薄板单元：常使用面积坐标，分析表明，只以挠度 及其一阶导数 作为节点的位移函数用一般的形状函数是不可能构造满足相容性的薄板单元，需再加上二阶导数，就可以实现。在相邻单元之间，挠度是连续的，但法向的斜率是不连续的，这种位移模式是非协调单云，收敛不如矩形单元，单元足够小，节点增多，如六节点三角形，九节点三角形等。5.3谈论在平面应力和弯曲状态组合的情况下，三角形刚度矩阵的特点（1） 平面内的作用力产生的变形不影响弯曲变形，反之亦然（2） 节点把转向 在两种应力状态下都不加入到变形中，相应的节点力也不存在，将平面应力状态和弯曲状态加以组合后，单元的每个节点的位移向量和节点力向量是 要指出的是，在局部坐标系中，节点位移不包括 ，但为了下一步将局部坐标系的单元刚度阵换到总体坐标系下进行集成，由于平面应力状态下的节点力和平面应力状态下的节点位移 互不影响，弯曲应力状态下的节点与平面应力状态下的节点位移互不影响，所以组合应力状态下的平板、薄板单元的单元刚度矩阵如下：，=其中矩阵和分别是平面应力问题和薄板弯曲问题的相应子矩阵，三角形单元的单元刚度矩阵是1818矩阵。6.1 结构的动态特性：结构的固有频率及其相应的模型，以及在随着时间而变形的外加激振力的激励下，机器或结构被激起的位移，应力或称被激起的动力响应，机械产品的动态性能是其重要的性能指标，尤其对现代复杂、高速、重载精密机械系统，动态性能是影响其工作性能及产品指标的关键技术指标，机械结构的动态特性问题早在上个世纪30年代就引起人们的重视，动态特性的发展为机械动态设计提供了坚实的基础。6.2 结构离散后，在运动状态各节点的动力平衡为：其中,分别以惯性力、阻尼力和动力载荷均为矢量，为弹性力，弹性力矢量可用节点位移和刚度矩阵表示为：=式中刚度矩阵的元素为节点j的单位位移在节点i引起的弹性力，根据达朗贝尔原理，可利用质量矩阵和节点加速度表示惯性如下：=式中质量矩阵为节点j的单位加速度在节点i引起的惯性力，设结构阻尼（滞粘），可用阻尼矩阵C和节点速度，表示阻尼如下：=，将各式带入：+=，记=，=。则运动方程：+=6.3单元的质量矩阵：= 质量矩阵是对称阵，各节点的质量互相耦合，即平动惯性和转动惯性之间耦合，如果把单元的一致质量集中的分配在它们的节点上，则此质量矩阵成为集中质量矩阵质量分配原则：按静力学平行力的分配法则，将单元的一致质量矩阵用集中于节点外的质量来代替，形函数计算所得的M称为一致质量矩阵。6.5 结构阻尼（只与结构本身材料性质有关）结构在自由振动过程中，如果没有能量的耗散，振动将永远保持由初始条件决定的振幅持续不停，但实际上，结构自由振动的振幅都会随时间而衰减，经过一定时间后，这是因为系统的能量因某些原因而消耗，这种能量的耗散作用称阻尼，由阻尼使振动衰减的系统称为阻尼系统。在结构内部阻尼是非粘线的，但它近似于线性的，弹性材料，特别是金属材料表示一种结构阻尼的性质，这种阻尼是由于材料受力变形而产生的内摩擦力和变形之间产生了相位滞后。产生能量耗散的原因有结构的内摩擦（或粘性）构件接口处的摩擦、周围介质（如空气、建筑物地基）的阻尼影响等，但有关阻尼的作用机理，目前尚未完全研究清楚。1.推导横截面积为A的一维桁架架构单元刚度矩阵。解：设杆件两端点位i，j，为单元局部坐标，表示单元任一截面的位置，则其发生的位移：u=a0+b1，v=b0+b1+b22+b33，即： u 1 0 0 0 0 = *（a0 b0 b1 a1 b2 b3）T v 1 0 0 2 2 H 记U=u,v=H* ，由i，j两端的位移分量可得：=G* ,1 0 0 0 0 0 0 1 0 0 0 0其中G= 0 0 1 0 0 0 给上式左乘G-1,则有 1 0 0 L 0 0 0 1 L 0 L2 L3 0 0 1 0 2L 3L2u=H* G-1*,令N= H* G-1N1=1-/L 0 0 /L 0 0, N2=0 1-3/L2+2/L3 *（1-/L）2 0 3/L2+2/L3 *（/L-1）*/L， 应用几何物理方程可得：= n = *=B* n 利用虚功原理推得：Ke=E*= EA/L 0 12EIZ/L3 对 0 6EIZ/L2 4EIZ/L 称 -EA/L 0 0 -EA/L 0 -12EIZ/L3 -6EIZ/L2 0 -12EIZ/L3 0 -6EIZ/L2 2EIZ/L 0 -6EIZ/L2 -EA/L2.如图2为一个平面超静定桁架结构，在载荷P的作用下，求各个杆的轴力。此结构可以看成由14,24,34三个杆组成的，每个杆单元的两端为杆单元的结点，各结点的水平，铅直位移分别用u、v表示。解：由题意可得：各杆件在局部坐标系下的单元刚度矩阵： 1 0 -1 0 0 0 0 0ke=EA/L -1 0 1 0 e=(14, 24, 34) 0 0 0 0图2 桁架超静定结构对于14杆转角=/2+，cos=-cos，sin=sin， sin -cos 0 0 cos sin 0 0故T14= 0 0 sin -cos 0 0 cos sin sin2 sin* cos -sin2 -sin* cos对于K14=T14 T*K14 *T14=EA/L* -sin* cos cos2 sin* cos -cos2 -sin2 -sin* cos sin2 sin* cos sin* cos -cos2 -sin* cos cos2对于24杆转角=90，则有： 0 0 0 0 0 1 0 -1K24= EA/L* 0 0 0 0 0 -1 0 1对于34杆转角=/2-，cos=cos，sin=-sin，cos -sin 0 0 -sin cos 0 0故T34= 0 0 -cos sin 0 0 -sin cos 对于K34=T34 T*K34 *T34 -sin2 sin* cos -sin2 -sin* cos 0 0 0 0 sin* cos cos2 -sin* cos cos2 0 0 0 0 sin2 sin* cos 0 0 0 0 -sin2 -sin* cos Ke=K14+K24+K34=EA/L* -sin* cos cos2 0 -1 0 -1 -sin* cos cos2 0 0 0 0 0 0 0 0 0 0 0 1 0 -1 0 0 0 0 -sin2 -sin* cos 0 0 sin2 sin* cos 0 0 -sin* cos -cos2 0 0 sin* cos cos2利用K*0 0 u v 0 0 0 0=0 0 0 P 0 0 0 0，可以解得u，v的值。对于杆单元14时，14=K14*q14可以求得；对于杆单元24时，24=K24*q24可以求得；对于杆单元34时，34=K34*q34可以求得。3.如图3所示的钢架中，两杆为尺寸相等的等截面杆件，横截面积为A=0.5m2，截面惯性矩为I=1/24m4，E=3*107kpa，求解此结构。图3 钢架解：将杆件单元标出单元号码及结点号码（如图所示），钢架的单元参数如下：单元数为2，结点数为3，各杆件子局部坐标系下的单元刚度矩阵： -1 0 -1 0Ke=EA/L* 0 0 0 0 e=1,2 -1 0 1 0 0 0 0 0 -1 0 -1 0对于单元转角=0，故K =EA/L* 0 0 0 0 -1 0 1 0 0 0 0 0 0 0 0 0对于单元转角=90，故K = EA/L* 0 1 0 -1 0 0 0 0 0 -1 0 1 1 0 0 0 -1 0 0 1 0 -1 0 0 0 0 0 0 0 0K=K +K =EA/L* 0 -1 0 1 0 0 -1 0 0 0 1 0 0 0 0 0 0 0由K*u v 0 0 0 0=14 -22 0 0 0 0,可以求得u，v。