第三章 等参单元

1、船舶结构有限元分析船舶结构有限元分析Finite Element Analysis of Ship Structure上海海事大学商船学院江国和第三章 等参单元3.13.1等参单元的引入等参单元的引入采用等参变换的单元称之为等参单元，即单元几何形状的变换和单元内的场函数，采用相同数目的节点参数，及相同的插值函数进行变换。基本思想：现在具有规则形状的单元（区域）上构造位移插值函数，然后把这个具有规则形状的坐标变换映射为物理平面上的一个形状比较复杂的单元，因而，等参单元也被称为映射单元。3.23.2四节点矩形双线性单元四节点矩形双线性单元3.21位移插值函数 对于某些有规则边界的区域，可以把所讨论

2、的区域划分成若干个矩形单元，设一个中心位于坐标原点，变长为2的正方形，如图 取正方形的四个角点作为单元节点，单元的位移模式为xyyxvxyyxu87654321 式中所示的位移插值函数的特点：固定x，它是y的线性函数；固定y，它是x的线性函数。因而，我们称它为双线性插值函数，把采用这种位移插值函数的四节点矩形单元称为矩形双线性单元。 矩形双线性单元有着常应变三角形单元所不具有的反映单元弯曲变形的能力，从而计算精度有望提高。 由于上式满足常应变准则即完备性条件，又满足位移连续性条件即协调性条件。因而，矩形双线性单元是完备的协调单元。3.2.23.2.2形函数的确定形函数的确定 建立以节点位移为基

3、本未知量的有限元方程，把u和v表示成形函数与节点位移乘积之和的形式 形函数确定的条件： （1）每个形函数都仍然是坐标的双线性函数； （2）Ni在结点i处等于1，在其他节点处等于0。4141iiiiiiuNvuNu根据条件2，可得根据矩形单元各个节点的坐标，可将形函数Ni统一表示为 （i=1,2,3,4）对任意双线性函数f(x , y),都有 为f(x , y)在节点i的值，分别取f(x , y)=1， f(x , y)=x和f(x , y)=y，得：yxNyxNyxNyxN11411141114111414321yyxxNiii114141,iiifNyxfif414141,1iiiiiiii

4、yNyxNxyxN 形函数是定义在给定单元上的，所以上式对单元的每一点（x , y）都成立。单元内一点的坐标通过单元节点坐标 的表达式，即坐标插值函数。 位移插值函数和坐标插值函数具有完全相同的构造，都是以节点处的函数值作为参数，并具有相同的形函数。因而矩形双线性单元也是一种等参数单元。),(iiyx3.2.23.2.2矩形双线性单元的完备性和协调性矩形双线性单元的完备性和协调性完备性：所采用的位移插值函数的确能反映任意给定的刚体位移和常应变，简称为满足常应变准则或完备性条件上式是任何等参数单元满足常应变准则的充分必要条件414141,1iiiiiiiiyNyxNxyxN协调性：只要相邻单元在

5、其公共边界上具有共同节点，相邻单元间公共边界上任一点的位移就必定相同，所以位移连续条件（即协调性条件）得到满足。在x=1边界上，N1=0，N4=0，有同理，在x=1的边界有 。可知，矩形双线性单元的边界在手里变形后，仍保持为直线边界。此协调性论证适用于任何等参单元。)1(2)1(21)1(21),1(),1(),1(232323322yuuuuyuyuyNuyNyuyvvvv122323.33.3四节点四边形等参单元四节点四边形等参单元前面介绍的四节点矩形双线性单元只适用于规则区域。本届介绍可用于不规则区域的四节点任意四边形单元。3.3.1四节点任意四边形单元对任意四边形单元，不能再直接采用双

6、线性插值函数。通过适当的坐标变换，把物理平面（xy平面）四边形单元变换成平面上以坐标原点为中心的边长为2的正方形单元，xy平面上的节点1，2，3，4分别对应正方形单元的节点1，2，3，4。如图所示对于a所示的四边形单元，将各边的等分点用之线连接，并规定它与局部坐标系下单元的相应对边、等分点的连线对应，这样，就建立了图a所示四边形单元和图b所示正方形单元个点间的一一对应关系。3.3.2自然坐标系下的位移插值函数自然坐标系下的位移插值函数：式中形函数 （i=1,2,3,4）有可知，插值函数在自然坐标系下满足常应变准则和协调性条件4141,iiiiiiuNvuNuiiiN1141,414141,1i

7、iiiiiiiNNN3.3.3坐标变换整体坐标和自然坐标之间的坐标变换关系式：由形函数在节点i为1，在其他节点为0的性质可知，上式把正方形的四个顶点变换成整体坐标下的节点（Xi，Yi）。其次，考察母单元某一条边，如12边=-14141,iiiiiiyNyxNx2121431211211211210,0yyyxxxNN上式是一条以为参数的直线方程，因而变换式把母单元的12边变换成整体坐标下连接节点的直边。对其他边也有同样结论。位移插值公式和坐标变换公式具有完全相同的构造，因此这种四节点任意四边形单元也是等参单元，简称为四节点等参数单元，在整体坐标下也必然满足常应变准则。位移插值函数在自然坐标系下

8、的协调性自动保证了坐标变换式的合理性及位移插值函数在整体坐标系下的协调性今后，对每个单元，我们把立足点放在自然坐标上，一切计算都转换到自然坐标系下进行位移插值函数可用矩阵记号记为：3.3.4应变矩阵B建立对平面应力问题eNvuvuNNNNNNNNvuf44114321432100000000exyyxBvuxyyx00上式B是一个三行8列矩阵 （i=1，2，3，4）利用复合函数求导法则，得写成矩阵形式，xNyNyNxNBiiiii00yyxxyyxxyxJyxyxyx式中坐标变换矩阵（雅可比矩阵）由矩阵求逆法则得式中41414141iiiiiiiiiiiiyNxNyNxNyxyxJ414141

9、41111iiiiiiiiiiiixNxNyNyNJxxyyJJxyyxJ3.3.5应力应变关系S=DB为应力矩阵3.3.6单元刚度矩阵及等效节点载荷1.面元变换公示自然坐标系下的无限小面积元，经过坐标变换式变换为整体坐标系下无限小平行四边形，如图eezyxSDBD变换后微元面积记为dA，则2.单元刚度矩阵把积分区域变换到自然坐标系平面上的相应区域 ddJddyxyxdAdyydxxdydydxdxxdyydxdA,dtdJDBBKDBtdxdyBKTeATe 1111A单元所占区域t单元厚度3.等效节点载荷dtdJQNPTeQ 1111讨论作用在单元边界上的表面里引起的等效节点载荷11122

10、11122111211112122111212121414141400002,12tdPlPtdPlPtdPlPtdPlPPNPNPNPNQNdldslsQtdsNPyeyxexyeyxexyxyxTTeQ在12边上，=-1，N3=0,N4=03.3.7等参变换的条件为了确保坐标变换式能确定整体坐标与自然坐标间一一对应关系，要求：所以，为保证等参数变换得以进行，所谓的任意四边形所有内角都必须小于180当内角接近180，丨J丨接近与零，因为丨J丨出现在求逆矩阵的分母，将导致计算误差迅速增大，因而计算实践不应使四边形单元过于歪斜。0,J3.4等参单元的收敛性为了保证协调性，相邻单元在这些公共边（或

11、面）上应有完全相同的节点，同时每一单元沿这些边(或面)的坐标和未知函数应采用相同的插值函数加以确定。关于单元的安全性，对于 型单元，要求插值函数中包含完全的线性项，本章讨论的所有单元在自然坐标中是满足此要求的。研究经过等参变换后，笛卡尔坐标下有：将上式代入，就得到单元内的函数表达式可以看出，如果插值函数满足 ，说明单元满足完备性要求0C1,111111niiniiiiiiniiiniiiniiiniiiNdzcybxNadzcybxadzcybxaNzNzyNyxNx某些情况下，坐标插值多项式次数低于位移插值多项式，计算更为便捷，这样的单元称为次参数单元，如图为典型的二维次参数单元。坐标插值函

12、数的阶数高于位移插值函数的阶数的单元，称为超参数单元，可以更准确的描述壳体的复杂几何形状3.5数值积分方法3.5.1插值求积法在求积区间-1，1上取n个分点i（i=1,2n）,n 个求积节点的分布可等距也可不等距，记根据节点上的函数值Fi，构造一个n-1次插值多项式Pn(),使得 ，称为F()的插值多项式或拉格朗日插值多项式。写成如下形式：根据Ni()的性质可知 ),2 , 1( ,niFFii iinFPiniinFNP1 niiiiiiiniiiN11211121记则有插值多项式可写为于是有 若记 就得到插值求积公式 niw21 iiiwwN iniiinFwwP1 nniiiniiiin

13、FHFHFHdFdwwHdFwwdPdF 221111111111111式中的系数Hi(i=1,2,3n)称为求积系数或加权系数。它的值与被积函数无关，只依赖于求积节点的个数和位置。当F()本身是次数不超过n-1的多项式时， 即具有n个求积节点的插值求积公式至少具有n-1次的代数精确度。3.5.2高斯求积法对求积节点的位置予以优选，使插值求积公式对任何次数不超过2n-1次的多项式函数都精确成立。这时求积公式的最高代数精确度为2n-1次，称为高斯求积公式 FPn设F()为次数不超过2n-1次的任意多项式，因此有只需要选取i(i=1,2,，n)，使以下n个等式成立，则插值求积公式即为高斯求积公式。

14、积分点的个数n，称为数值积分的阶 1110nnnwPF dwFdwwdFnkkkiniii 101111111 111, 2 , 1 , 0, 0nkdwk 11111),2,1(,nidwwHFHdFiiiniiiN次勒让德多项式在区间-1，1上的n个根就是高斯求积公式中的n个求积节点。利用此性质可确定节点i和系数Hi。因此高斯求积公式又叫高斯勒让德求积法。 n=1情形，得一阶高斯求积公式式中1=0，H1=2，此为中点矩形插值公式，对任意一次多项式都精确成立 n=2情形，二阶高斯求积公式式中 ,它对任意三次多项式函数都精确成立。nnnnnddnL)1(!212 1111FHdF 221111

15、FHFHdF1,31,312121HH 下表列出n=1至n=4时高斯求积公式中的求积节点坐标和相应的求积系数3.5.3二维和三维高斯求积公式因为故二维高斯求积公式为：同理，三维高斯求积公式为：对有限元分析而言，由于几何各向同性要求，个方向的求积节点个数取为相同，故有上述公式 niiiFHdFddFddF11111111111, 111111,ninijijiFHHddF nininikjikjiFHHHdddF111111111,3.6数值积分阶次的选择选择积分阶次的原则：1.保证积分的精度缩减积分：高斯积分阶数低于被积函数所有项次精确积分所需要阶数的积分方案。采用缩减积分可以取得较完全精确积

16、分更好的精度。理由如下： .精确积分常常是由插值函数中非完全项的最高次方要求，而决定有限元精度的是完全多项式的次方。 .有限元计算模型具有较实际结构偏大的整体刚度，缩减积分方案将使模型刚度有所降低，因此提高计算精度3.23.2坐标变换与平面四结点等参元坐标变换与平面四结点等参元 图3-1(a)为一个任意四边形单元，称为实际单元。在实际单元内以对边的中点连线建立起一个局部坐标系，通过坐标转换把实际单元“映射”为如图3-1(b)所示的一个正方形，此坐标系称为单元的自然坐标系或等参数坐标系，正方形称为基本单元，基本单元内任一点P( ， )与实际单元内的一点P(x，y)唯一对应。图图3-1(a)3-1

17、(a)(a) 直角坐标系与实际单元 (b) 自然坐标系与基本单元图图3-1(b)3-1(b)实际单元与基本单元的对应关系实际单元与基本单元的对应关系可写为可写为(3-1)(3-1)(3-2) 用同样的形状函数来插值单元内任意一点用同样的形状函数来插值单元内任意一点( (x, yx, y) )的位移的位移(3-3) 为此单元的结点位移列向量为此单元的结点位移列向量 为形状函数矩阵。为形状函数矩阵。这里采用了同样的形状函数式这里采用了同样的形状函数式(3-2)(3-2)，用同样的结点插值函数表示出单，用同样的结点插值函数表示出单元的几何坐标元的几何坐标x x、y y与与u u、v v，这种单，这种

18、单元称为等参单元。元称为等参单元。 e N N类似地可以推广到具有更多结点的单元，如类似地可以推广到具有更多结点的单元，如图图3-23-2所示所示(a) 直角坐标系与实际单元 (b)自然坐标系与基本单元 图3-2 八结点等参单元该基本单元的位移函数可取为该基本单元的位移函数可取为(3-4)其中在顶角结点与边中点上的形函数在顶角结点与边中点上的形函数分别为(3-5)(3-6)3.3 3.3 单元刚度矩阵单元刚度矩阵首先给出单元内的应变列向量，对平面问题首先给出单元内的应变列向量，对平面问题应有应有 exyyxxyyxxyyxNvu0000(3-7)(3-7)按坐标变换关系式(3-1)，有 yyN

19、xxNNiiiyyNxxNNiii写成矩阵表达式为： yNxNJyNxNyxyxNNiiiiii(3-8)(3-8)由式由式(3-8)(3-8)可解出可解出 iiiiNNJyNxN1 yxyxJ其中其中(3-9)JJ称为称为坐标变换的雅可比（坐标变换的雅可比（JacabianJacabian）矩阵）矩阵，其中其中iixNxiiyNyiixNxiiyNy合写成矩阵形式有合写成矩阵形式有 44332211432143yxyxyxyxNNNNNNNNJii(3-10)(3-10)将式(3-3)代入式(3-7)中，则有 eezyxBNNNNNNNNxyyx432143210000000000 B为应变

20、转换矩阵应变转换矩阵，按结点分块表示，有 4321BBBBB xNyNyNxNBiiiii00而而 i i=1,2,3,4=1,2,3,4(5-11)(5-11)将式(3-9)代入式(3-11)，即可得出此单元的应变转换矩阵 ,进而求出 。 B同上，同上，单元内的应力单元内的应力可表示为可表示为 ezyxBD etdBDBkeT(3-12)(3-12)上述积分在自然坐标系内进行，得上述积分在自然坐标系内进行，得 ddJddyxyxd刚度矩阵刚度矩阵 1111TddJtBDBke(3-13)(3-13)单元刚度矩阵由虚功原理求得，即单元刚度矩阵由虚功原理求得，即式中，式中， 为对应的坐标位置为对

21、应的坐标位置 、 值，值， 、 为权重系数，为权重系数，L L、M M为沿为沿 、 方向的积分点数目。方向的积分点数目。)(iifiiiwjw 一般参数单元的计算都采用数值积分求式(3-13)的近似值，同时，为了减少计算点的数目和便于编写程序，多采用高斯数值积分方法。二维积分二维积分法的高斯求积公式法的高斯求积公式为jiiiMjLiwwfddfI)()(111111 (3-14)(3-14)3.4 ANSYS3.4 ANSYS平面结构计算示例平面结构计算示例3.4.1 3.4.1 问题描述问题描述 一个长方形面板，如图5-3所示，其高AB=1m，宽BC=1.5m，板厚b=0.04，孔半径R0.

22、2m，长方形面板的弹性模量E=210GPa，泊松比=0.3，约束条件：在长方形底边约束全部自由度。载荷：BC边施加垂直向下均布载荷F10000000N/m。图3-3 长方形板结构 3.4.2 ANSYS3.4.2 ANSYS求解操作过程求解操作过程 打开打开AnsysAnsys软件，在软件，在AnsysAnsys环境下做如下操环境下做如下操作。作。 （1 1）选择单元类型选择单元类型运行运行PreprocessorElPreprocessorElement TypeAdd/Editement TypeAdd/Edit/Delete/Delete，弹出，弹出ElemenElement Typet

23、 Type对话框，如图对话框，如图3 3-4-4所示。单击所示。单击AddAdd，弹，弹出出Library of ElemenLibrary of Element Typet Type窗口，如图窗口，如图3-53-5所示，选择所示，选择PLANE82PLANE82。图3-4 单元类型对话框 图3-5单元类型库对话框 在在Element TypesElement Types对话框中单击对话框中单击OptionsOptions对话框，弹出对话框，弹出如图如图3-63-6所示对话框，设置所示对话框，设置K3K3选项栏为选项栏为Plane strsw/thkPlane strsw/thk，设置设置K5

24、K5选项栏为选项栏为Nodal stressNodal stress，设置，设置K6K6选项栏为选项栏为No extNo extra outputra output。表示单元是应用于平面应力问题，且单元。表示单元是应用于平面应力问题，且单元是有厚度的。是有厚度的。 图图3-6 PLANE82 3-6 PLANE82 单元类型选项对话框单元类型选项对话框 图5-7实常数对话框图5-8选择要设置实常数的单元类型 图5-9 PLANE82实常数对话框 运行运行PreprocessorReal Constants Add/Edit/DPreprocessorReal Constants Add/Edi

25、t/Deleteelete。弹出如图。弹出如图3-73-7所示对话框，点击所示对话框，点击AddAdd，弹出如图，弹出如图3 3-8-8所示对话框，点击所示对话框，点击OKOK，弹出如图，弹出如图3-93-9所示对话框。在所示对话框。在T THKHK选项栏中设置板厚度为选项栏中设置板厚度为0.04m0.04m。设置完毕单击。设置完毕单击OKOK按钮按钮完成设置。完成设置。图图3-103-10选择材料属性对话框选择材料属性对话框 图图3-11 3-11 设置材料属性对话框设置材料属性对话框 （2 2）设置材料属性设置材料属性运行运行PreprocessorMaterial Props Mater

26、ial ModelPreprocessorMaterial Props Material Models s，弹出如图，弹出如图3-103-10所示对话框，双击所示对话框，双击IsotropicIsotropic，弹出，弹出图图3-113-11所示对话框，在所示对话框，在EXEX选项栏中设置数值选项栏中设置数值2.1e112.1e11，在在PRXYPRXY选项栏中设置数值选项栏中设置数值0.30.3。设置完毕单击。设置完毕单击OKOK按钮按钮。 （3 3）建立模型建立模型 选择选择Preprocessor Modeling Create Area RectangleBy 2 CornersPre

27、processor Modeling Create Area RectangleBy 2 Corners；弹出如图弹出如图3-123-12所示对话框，设置参数，所示对话框，设置参数，WP XWP X选项栏中填写选项栏中填写0 0，WP YWP Y选项栏中填写选项栏中填写0 0，WidthWidth选项栏中填写选项栏中填写1.51.5，HeightHeight选项栏中填写选项栏中填写1 1，单击，单击OKOK，设置完毕。继续运行，设置完毕。继续运行P Preprocessor Modeling Create Area CircleSolid Circlereprocessor Modeling

28、Create Area CircleSolid Circle；得到如图；得到如图3-133-13所示对话框，在所示对话框，在WP XWP X选项栏中填写选项栏中填写0.750.75，WP YWP Y选项栏中填写选项栏中填写0.50.5，在，在RadiusRadius选项栏中填写选项栏中填写0.20.2，设置完毕点击，设置完毕点击OKOK按钮。得到如图按钮。得到如图3-143-14所示图形。所示图形。 图图5-12 5-12 依两点建依两点建立矩形对话框立矩形对话框图图5-13 5-13 创建实心圆创建实心圆对话框对话框图图5-14 5-14 长方形板模型长方形板模型 （4 4）划分网格划分网格

29、运行运行MeshingSize Cntrls ManualSizeAreasAll AreMeshingSize Cntrls ManualSizeAreasAll Areasas，弹出如图，弹出如图3-153-15所示对话框，在所示对话框，在SIZESIZE选项栏中填写选项栏中填写0.050.05，点击点击OKOK按钮；运行按钮；运行MeshAreasFreeMeshAreasFree划分网格，网格划分划分网格，网格划分如图如图3-163-16所示。所示。图5-15 设置网格尺寸对话框 图5-16划分网格后板的有限元模型 图图5-17 5-17 对线施加全约束对线施加全约束 （5 5）施加约

30、束施加约束 选择菜单选择菜单SolutionDefine LoadsApplyStructurSolutionDefine LoadsApplyStructureDisplacementOn LineseDisplacementOn Lines，选择长方形底边，弹出，选择长方形底边，弹出图图5-175-17所示对话框，选择所示对话框，选择All DOFAll DOF，单击，单击OKOK。 （6 6）施加载荷施加载荷 选择菜单选择菜单SolutionDefine LoadsApply StructSolutionDefine LoadsApply StructurePressureOn LinesurePressureOn Line