复合材料层合板的杂交有限元方法

卿光辉，贾瑞升（中国民航大学航空工程学院，天津 300300）

复合材料层合板的杂交有限元方法

卿光辉，贾瑞升
（中国民航大学航空工程学院，天津 300300）

结合复合材料修正后的H-R混合变分原理，直接借助应力-应变关系，推导了新的应力模式，建立了复合材料层合板的杂交等参有限元列式。利用Mathematica语言编程进行数值实例分析，其计算结果与相关文献的精确解以及Abaqus软件建模分析结果对比，实例证明该方法所得到的各个静力学量更接近精确解，并且可用较少的网格划分得到较精确的解。

复合材料；杂交有限元；等参元；Mathematica

近年来复合材料由于其结构本身的特性，已经在航空领域得到了广泛的应用和发展。20世纪60年代初期假设位移法对于平板弯曲的问题遇到了困难，主要是假定的形函数不容易在边界上满足C1连续性条件，即可以保证位移的连续性，但不容易保证转角的连续性。而杂交应力有限元以假定单元的应力试解为基本特色，很容易实现单元间位移协调性的要求，这对于克服板壳中要求挠度一阶导数连续的困难有重要的意义。1980年以后，Chen和T.H.H.Pian[1]运用等参坐标来表示单元应力，使用Hellinger-Reissner原理推导时，可以不必考虑平衡条件[2-4]。本文在前人工作的基础上，直接借助应力-应变关系，推导了新的应力模式，建立了复合材料层合板的杂交有限元和等参元及其数值分析模型，并利用Mathematica语言编程进行了具体数值实例的分析，数值结果与Abaqus有限元软件模型和相关文献的精确解进行了对比，发现本文的方法对网格的划分不是很敏感，即网格划分数目较少就可以得到比较精确的解，并且随着网格的加密计算所得的结果变化不大。

1 变分原理和状态方程

1.1 修正的H-R变分原理

假设复合材料是各向同性或正交各向异性的，其本构关系可由下式给出

其中：σx、σy、σz、σyz、σxz、σxy分别为应力分量；C11、C12、C13、C22、C23、C33、C44、C55、C66是复合材料的刚度系数；Sx、Sy、Sz、Syz、Sxz、Sxy是应变分量。

应变-位移关系为

其中：u、v、w为沿x、y、z方向的位移分量。

在单元内假设独立的应力场

其中：P={σ1σ2… σM}为应力矩阵；σi为假设应力模式；βi为对应的应力参数。

在单元内部假设的位移场

为三维八节点等参单元的形函数，在这里自然坐标系ξ，η，ζ∈[-1，1]，而ξi、ηi、ζi是i节点在自然坐标系中相应节点的坐标值，i代表某一层中三维八节点单元的节点号，j代表层合板单元的某一层。

引入记号

是每层子单元的位移矢量。

当位移已得到满足，赫林格-赖斯纳（Hellinger-Reissner）广义变分原理的泛函为

其中：Ωe为单元区域；S为材料柔度矩阵；b为体力；t为面力；D为微分操作符；Γt为应力边界区域。

把式（3）、式（5）代入式（6）可得

其中：H为单元柔度矩阵；G为杠杆矩阵；fe为等效节点载荷，分别为

故式（9）可简化为

其中：单元刚度矩阵

对式（14）泛函对qe变分，即利用驻值条件0可以得到位移的矩阵状态控制方程从而可以利用一般有限元的分析方法求得未知节点的位移qe，进一步由式（4）求得单元位移场，并且由式（3）和式（13）求得单元的应力场

1.2 假设应力模式

因为形成杂交应力元刚度矩阵时必须对柔度矩阵H求逆，所以杂交应力元存在的充分必要条件是柔度矩阵可逆，即柔度矩阵是非奇异的，也就是要求假设的应力模式线性无关。

下面根据等函数方法给出三维八节点正交异性复合材料单元的一种假设应力模式的推导[5-6]。

对于一般的三维八节点单元，可假设位移场如下

在自然坐标系中，依据位移-应变关系，可得

将式（19）代入式（1）可以求得

将式（19）代入式（20），调整系数后可以得到

σξ=β1+β2ξ+β3η+β4ηζ+β5ξζ+β6ξη

即

因此该假设应力模式为

2 实例分析

例题1 有一悬臂梁，梁长150 mm，宽2.5 mm，高5 mm；一端固定；自由端中心承受5 N的载荷；材料的杨氏模量E为70 GPa，泊松比为0[7]。其计算结果如表1所示。

从表1的计算结果可以看出本文的解与文献[7]的解析解的误差很小，明显好于使用Abaqus有限元软件计算的结果，并且该方法使用较少的网格划分（2×12）就可以得到很精确的解。

表1 悬臂梁自由端的最大挠度Tab.1 The maximum deflection of the cantilever

例题2 考虑一三层板，且a=b，在x=0，a处固支，在y=0，b两边简支，板的总厚度为h，h1=h3= 0.01h，h2=0.08h，第一层、第三层的材料参数相同，第二层材料参数为，板的上表面受均布压强q的作用。材料参数如下

计算结果如表2所示，其中I1、I2和I3分别为第一、第二和第三层所分的薄层数，从表1可以看出本文的解与文献[8]的精确解比较吻合，并且显然好于Abaqus采用的C3D8R单元类型有限元计算的结果。

表2 三层层合板的挠度和应力Tab.2 Three laminated plates of deflection and stress

3 结语

通过修正的H-R混合变分原理，直接借助应力-应变关系，推导了新的应力模式。在此基础上，建立了复合材料杂交等参元的有限元列式，计算步骤简单，数据累计误差小，提高了计算精度，同时，我们发现该方法对网格的划分不是很敏感，即可用较少的网格划分就能得到比较精确的解，随着网格的加密所计算的结果变化不大，因此该方法不仅可节约大量的计算时间，而且还能降低硬件的配置，从而节省了资源。

[1]卞学璜.杂交应力有限元法的研究进展[J].力学进展，2001，31（3）：344-349.

[2] SZE K Y，GHALL A.Hybrid hexahedral element for solids，plates，shells and beams by selective scaling[J].International Journal of Numerical Methods in Engineering，1993，36：1519-1540.

[3]CHEUNG Y C，DI SHENGLIN.Analysis of laminated composite plate by hybrid stress isoperimetric element[J].International Journal of Solids Structures，1993，30（20）：2843-2857.

[4] SPILKER R L.Hybrid stress eight-node elements for thin and thick multilayer laminated plates[J].International Journal for Numerical Methods in Engineering，1982，18：801-828.

[5] FENG W，HOA S V，Huang Q.Classification of stress modes in assumed stress fields of hybrid finite elements[J].International Journal for Numerical Methods in Engineering，1997，40：4313-4339.

[6]张东锋.压电层合板的杂交元方法[D].成都：西南交通大学，2002.

[7]庄 茁，张 帆，岑 松，等.ABAQUS非线性有限元分析与实例[M].北京：科学出版社，2005.

[8]范家让.强厚度叠层板壳的精确理论[M].北京：科学出版社，1996：130-351.

（责任编辑：杨媛媛）

Hybrid finite element method for laminated composite plate

QING Guang-hui，JIA Rui-sheng
（College of Aeronautical Engineering，CAUC，Tianjin 300300，China）

In this paper based on modified H-R mixed variational principle for composite materials，with the stressstrain relations directly，derivation of a new mode of stress，the hybrid and isoparametric finite element formulation for the laminated composite plate was established.Then the Mathematica was applied for the programming and calculating of a numerical example.Compared with the modeling analysis result using Abaqus software and the exact solution provided in relevant literatures concerning some mechanical quantities，the result obtained in this way is proved to be closer to their exact solutions and satisfactory precision can be obtained with less mesh.

composite materials；hybrid finite element；isoparametric element；Mathematica

O242.21

A 文章编号：1674－5590（2013）01－0082－03

2012-06-08；

2012-07-11

中国民航大学科研基金项目（2012kye07）

卿光辉（1968—），男，湖南新化人，教授，博士，研究方向为结构力学.