有限元作业第二次作业

土木工程专业有限元第二次作业姓名：班级：学号：指导教师：二〇一五年6月12日习 题：平面应力问题的八节点等参元，已给定8个节点的坐标。

试查资料并论述：1、单元中位移函数u （ξ，η），v （ξ，η）和单元节点位移{δe }的关系式；2、[ B ]矩阵的计算步骤和计算式；3、单元刚度矩阵[ k e ]的一般计算方法和计算步骤；4、论述相邻单元间公共边界上位移的连续性；5、如果给定母单元中点A ，（ξ，η），怎样求实际单元中与A ，相对应的点A （x ，y ）；反之，如果给定实际单元中的点A （x ，y ），怎样求其在母单元中对应点A ，（ξ，η）？ 6、如果已经求解得到单元8个节点的位移值{δe }怎样求单元中某一点B （x ，y ）的应力？实际单元1267Y1243678η= 1η=﹣1母单元ξ= 1ξ=﹣1解：1、此题分两步进行：➢ 单元位移场的表达：如图1所示,在任意四边形的每边中间设一附加节点，则单元边界就变成二次曲线的了。

如果直接在整体坐标系(),x y 下，像八节点矩形元那样，构造双二次多项式的位移插值函数，则因曲边四边形单元边界是二次曲线，故边界上的位移是()x y 或的五次多项式，它不能由曲边上三个节点的位移分量唯一地决定，从而不能保证相邻两个单元在公共边上位移的协调条件，所以在整体坐标系(),x y 下构造完全协调的位移插值函数是很困难的，利用坐标变换，可将曲边四边形单元变换成基本单元，如图2所示的在自然坐标(),ξη下具有边长为2的八节点正方形单元，自然坐标系(),ξη是外节点坐标值为±1的局部坐标系。

在自然坐标系的单元上构造协调的位移插值函数，其形状函数是较普通的，取位移分量为,ξη的双二次多项式, 即：2222123456782222910111213141516u a a a a a a a a v a a a a a a a a ξηξξηηξηξηξηξξηηξηξη⎧=+++++++⎪⎨=+++++++⎪⎩（1-1） 利用8 个节点的16 个位移分量可唯一确定16 个待定常数1216,,a a a …,，图1：在总坐标系中具有二次曲边的四边形单元图2：在自然坐标系中的曲边四边形的基本单元若代入8个节点的局部坐标值，得：11523264536774881-1-1111-1-110-10010011-11-11-1-11101000011111111101001001-111-111-11-1010000u a u a a u u a a u a u a u u a ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦（1-2）195102116123137141548161-1-1111-1-110-10010011-11-11-1-11101000011111111101001001-111-111-11-1010000v a v a v a v a v a v a a v v a ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦（1-3）将解出的16 个待定常数1216,,a a a …,代入式（1-1）即得:811552266337744881811552266337744881i i i i ii u N u N u N u N u N u N u N u N u N u v N v N v N v N v N v N v N v N v N v==⎧=+++++++=⎪⎪⎨⎪=+++++++=⎪⎩∑∑ （1-4a ） 也即：[]{}{}128e eu N N N v δδ⎧⎫===⎨⎬⎩⎭u III N （1-4b ）其中I 为二阶单元矩阵，{}eδ为等参元节点位移列阵，N 为形状函数矩阵。

➢ 形状函数的建立：按等参元思想，在整体坐标系XY 下, 任何形状歪斜四边形单元都将变换到局部坐标系ξη下的正方形单元。

对8节点等参元, 其移模式为：()81,iii u N u ξη==⋅∑ （1-5）式中, i u 为歪斜单元8节点的位移，(),i N ξη为形状函数。

查阅相关资料，得形函数公式公式为：()()()8181,,,kk i kiik F N F ξηξηξη===∏∏ （1-6）又由形状函数的性质可具体地求出i N 的表达式为：()()()()()()()()()()()()()()()()()()()()123425262728=1114=1114=1114=1114=112=112=112=112N N N N N N N N ξηξηξηξηξηξηξηξηξηηξξηηξ⎧-----⎪+---⎪⎪+++-⎪-+-+-⎪⎪--⎨⎪-+⎪⎪⎪-+⎪⎪--⎩（1-7）2、根据平面问题的几何方程，单元应变可用节点位移表示如下：{}[]{}128=x e ey xy εεδδγ⎧⎫⎪⎪=⎨⎬⎪⎪⎩⎭ε=B B B B （2-1）其中：0=0i i i i i N x N y N N y x ⎡⎤∂⎢⎥∂⎢⎥⎢⎥∂⎢⎥∂⎢⎥⎢⎥∂∂⎢⎥∂∂⎣⎦B （2-2）即要求出矩阵i B 中的元素i N x∂∂，iN y ∂∂(1,2,,8)i =。

另根据符合函数求导法则，可知：=i i i i i i N x y N N x x N N N xy y y ξξξηηη∂∂∂⎧⎫∂∂⎡⎤⎧⎫⎧⎫⎪⎪⎢⎥⎪⎪⎪⎪∂∂∂∂∂⎪⎪⎪⎪⎪⎪⎢⎥=⎨⎬⎨⎬⎨⎬∂∂∂∂∂⎢⎥⎪⎪⎪⎪⎪⎪⎢⎥∂∂⎪⎪⎪⎪⎪⎪∂∂∂⎩⎭⎩⎭⎣⎦⎩⎭J （2-3） 其中，J 为二维坐标变化下的Jacobi 矩阵，即：xy x y ξξηη∂∂⎡⎤⎢⎥∂∂⎢⎥=∂∂⎢⎥⎢⎥∂∂⎣⎦J （2-4） 其元素计算式为：81i i i N x x ξξ=∂∂=∂∂∑，81i i i N x x ηη=∂∂=∂∂∑，81i i i N y y ξξ=∂∂=∂∂∑， 81i i i N y y ηη=∂∂=∂∂∑ （2-5） 又根据式（2-3），有1i i i i N N x N N y ξη-∂⎧⎫∂⎧⎫⎪⎪⎪⎪∂∂⎪⎪⎪⎪=⎨⎬⎨⎬∂∂⎪⎪⎪⎪∂⎪⎪⎪⎪∂⎩⎭⎩⎭J （2-6） 根据公式（2-2）即可得出i B 矩阵，其中i N 可由问题1方法求出。

3、单元刚度矩阵按普遍公式计算，公式如下：[]eeeT T k dV hdxdy ΩΩ==⎰⎰B DB B DB （3-1）其中e Ω为单元体积域，[]ek 为16×16的方阵（具体形式见下文），D 为材料的弹性系数矩阵，各向同性材料的弹性系数矩阵为：()()()101-1101121-1-2002-E μμμμμμμμμ⎡⎤⎢⎥⎢⎥-⎢⎥=⎢⎥+-⎢⎥⎢⎥⎢⎥⎣⎦（1）D （3-2） 上述积分应在局部坐标系内进行，因此面积元素dxdy 需表示成d d ξη.如图3所示为子单元内任一点(),a x y 处的微小正方形，它是由局部坐标系中点(),ξη处的微元体d d ξη变换而成的。

以,i j 表示,x y 轴的单位基矢量，,ab ac 分别由,d d ξη变换而成，则：x yd ξξξ⎛⎫∂∂+⎪∂∂⎝⎭ab =i j x y d ηηη⎛⎫∂∂+⎪∂∂⎝⎭ac =i j （3-3） 上述2个矢量的叉积表示它们所构成的平行四边形面积，故：dV Jd d ξη=⨯=ab ac （3-4）其中，J 为矩阵J 的行列式，即x y J x y ξξηη∂∂∂∂=∂∂∂∂ 将上式带入式（3-1），并写成分块形式：[]111218212228818288e k ⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎣⎦k k k k k k k k k （3-5） 其中子矩阵的计算公式为：1111eT T ij ij i j k hdV Jhd d ξη--Ω==⎰⎰⎰B DB B DB （3-6）其中h 是板的厚度。

由于被积函数极为复杂,很难得到明显的解析式,必须利用数值积分。

程序中采用高斯求积法,对于二维问题的等参元,高斯求积公式为:图3：子单元内任一点处的微小正方形()()88111111,,ijiii i f d d H H f ξηξηξη--===∑∑⎰⎰ （3-7）式中，i H ，j H 为一维求积点的积分系数，i ξ，i η为沿一维编号的求积积分点的横坐标。

对于8节点等参元取三个积分点，即n=3已足够精确。

4、证明：局部坐标系下的单元是规则的正方形，单元边界上的三个节点按线性变化的位移形式，单元变形后这三个节点确定了位移的单元直线边界。

所以，局部坐标系下单元是协调的。

又由位移插值函数在局部坐标系下的协调性，即可推知坐标变换的协调性（即两个相邻曲边四边形在公共边界上经坐标变换后仍保持连续，不会出现重叠和破缺现象），这也就保证了位移插值函数在整体坐标系下的协调性。

即在相邻单元公共边界上位移是连续的。

5、这里，,u v 是,ξη的函数，在下面的计算中还需知道,u v 和,x y 的关系，因此必须写出,x y 和,ξη之间的坐标变换式，这个坐标变换并不难，因为,x y 在单元的8个节点上应取值,(1,28)i i x y i =…，，而单元四条边应为二次曲线，这与,u v 的要求完全类同，因此可沿用和位移插值函数完全相同的式子作为坐标变换式，即：111526374815263748880000000000000000x y N N N N N N N N x N N N N N N N N y x y ⎧⎫⎪⎪⎪⎪⎡⎤⎧⎫⎪⎪=⎨⎬⎨⎬⎢⎥⎩⎭⎣⎦⎪⎪⎪⎪⎪⎪⎩⎭… （5-1） 式中112288,,,,,x y x y x y …,为节点坐标，形状函数128,,N N N …,与前面相同。

由上可见，在整体坐标系下的曲边四边形单元和自然坐标系下的正方形单元存在着一一对应的映射关系，只要已知,(1,28)i i x y i =…，后，由（5-1）式，利用自然坐标系下的形状函数，即可完全确定(),x y 。

即：如果给定母单元中点()*,A ξη，通过求出形状函数(),(1,28)i N i ξη=…，，利用式（5-1），可求出实际单元中与*A 相对应的点(),A x y ；同理，如果给定实际单元中的点(),A x y ，利用式（5-1），即可求出()*,Aξη的坐标值ξη和。

6、根据平面问题的本构方程，单元应变可用节点位移表示如下：{}{}==x e e yxy σσδδτ⎧⎫⎪⎪=⎨⎬⎪⎪⎩⎭=D DB S σε （6-1） 由式（3-2）和式（2-2）可分别求出矩阵D 和B ，故由上式，若已知单元8个节点的位移值{}eδ，可求出单元中某一点B （x ，y ）的应力σ。