基于混合元的变刚度层合板静力学分析

杨立洲，卿光辉

(中国民航大学航空工程学院，天津 300300)

在工程结构有限元仿真领域，众多学者致力于探索高效简洁且精确可靠的数值分析模型，例如，FRP(Fiber-reinforced Polymers)层合板的降维处理模型[1]、分层壳单元高效的非线性分析模型[2]和用于变刚度层合板的静力学与动力学分析的先进模型等。

20 世纪90 年代初，HYER 等[3]首先提出变刚度层合板(Variable Stiffness Laminates)的概念。变刚度层合板的纤维铺放角度沿某一方向不断变化，这极大增加了变刚度层合板的可设计性。但是，其应力-应变关系复杂，理论分析与数值仿真变得困难[4]。

部分学者采用解析法[5-9]对变刚度层合板结构的力学问题求解，虽然计算效率高，但需要对模型有较多的简化。从弹性力学的基本方程出发，附加若干假设(如KIRCHHOFF-LOVE 假设、KIRCHHOFF 假设等)的板壳理论对变刚度层合板的分析较为多见。如GUPTA 等[10]基于改进的三阶剪切变形理论，建立了用于变刚度层合板壳线性和非线性静力分析的有限元公式，研究了曲率比、层数和边界条件等因素对薄壳结构非线性静力性能的影响。谭萍等[11]基于板壳理论，采用加权残值法对纤维角度沿径向变化的变刚度圆环板的挠度变化进行相关研究。陈晓东等[12]基于经典层合理论，研究了不同边界条件下纤维角度变化时变刚度层合板的挠度及平面内节点位移受到的影响。板壳理论中的经典层合理论和高阶剪切变形理论，都假定板的挠度与z坐标无关，忽略了挤压变形的影响。由于中厚板剪切刚度小，存在较大的剪切变形，会导致经典层合理论的预测结果精度较低。当然逐层理论可以应用于中厚板分析，但基于该理论的模型变量个数依赖层合板的层数[13]。

YAN 等[14]基于改进的分层Legendre 展开式和微分求积有限元法，采用Carrera Unified Formulation模型，提出了一种新的统一的准三维解，对变刚度层合板进行振动分析。HASIM 等[15]基于精化锯齿理论(Refined Zig-zag Theory, RZT)建立了一种新的等几何公式，并应用于变刚度层合板和夹芯板的静力分析。该方法将RZT 与非均匀有理B 样条公式融合，避免了单元模型出现剪切自锁等问题。

广义混合元法[16-17]同时引入了应力边界条件，使得应力精度明显提高。避免了经典混合元法的整体刚度矩阵包含零子矩阵的问题，提高了计算过程的稳定性。QING 等[18-19]建立了非协调广义混合元，刘艳红等[20]将该方法扩展至热弹性复合材料，建立了一种含参数的非协调辛元。王燮等[21]建立了层合板的非协调广义部分混合元模型，分析了层合板的自由边界效应问题。王聿航等[22]扩展了非协调广义部分混合元，用于压电材料静力学分析，使用较少的单元数便可获得了和精确解一致性很高的结果。

考虑到解析法、板壳理论以及位移有限元法的局限性和变刚度层合板的特点，本文为变刚度层合板的静力学分析提出一种修正的非协调广义部分混合元模型。

1 非协调广义部分混合元列式

1.1 常规的非协调广义部分混合元列式

含参数的广义混合变分原理：

式中：V、S分别为层合板模型的体积和表面积；σ =[σxzσyzσzσxσyσxy]T为应力向量；u=[ux uy uz]T为位移向量；C为材料的刚度矩阵；∇为微分操作算子矩阵；T=[TxTyTz]T为给定的应力边界条件。式(1)中参数α 的取值范围为0～1，在后文的算例分析中α 取0.75[23]。

将位移和应力分别表示为：

式中：Nd和Ns分别为位移和应力的形函数矩阵；Nr为非协调项的形函数矩阵；qe和pe分别为单元节点的位移向量和应力向量；re为单元内部节点的非协调位移向量。

将式(2)和式(3)代入式(1)后进行变分处理有：

式中：

其中：

采用文献[21]中方法消去式(4)的第一个方程中的面内应力向量 [σxσyσxy]T后，可得到只包含面外应力和位移的部分混合列式：

式中的相关符号表达式见参考文献[21]。

1.2 修正的非协调广义部分混合元列式

对传统定向纤维层合板的离散处理，因为其单层板内纤维角度不变，可以直接将单层板的材料性质赋予每个单元，即每个单元具有统一的材料性质。因此对于传统的层合板结构，单元角度变化的偏轴应力-应变关系为[24]：

式中： σ′为单层材料主方向坐标系下的应力分量；T为坐标转换矩阵。

而变刚度层合板与传统定向纤维层合板不同，如图1(a)中的模型以及纤维路径函数所示，变刚度层合板的铺层为连续曲线纤维，即纤维角度θ 连续变化，且θ 是关于x的函数。对变刚度层合板作如图1(b)所示离散处理后，变刚度层合板的偏轴应力-应变关系式如下：

图1 变刚度层合板及有限元模型Fig.1 Variable stiffness laminate and its finite element model

其中：

x为几何模型离散后第n列单元中间处的横坐标，坐标原点为板的中心，x轴方向为长边a的方向，L为总的单元列数。

考虑变刚度层合板中单元材料参数刚度的变化情况，在推导列式(4)或式(5)时，按文献[15，25]中通用且可靠的方法，取单元中间位置的纤维角度作为该单元的纤维角度，如图1(b)所示。将式(7)中的 (T(θ(x)))-1C((T(θ(x)))-1)T作为单元的材料刚度参数代入列式(4)或式(5)中。因此，在本文用于数值分析的Wolfram 语言代码编写中，对常规的非协调广义混合元列式(4)和非协调部分广义混合元列式代码进行了修改，即根据单元的几何位置修改单元的材料刚度参数。

为了区分常规的非协调广义部分混合元，本文称通过以上处理后的有限元列式为修正的非协调广义部分混合元列式。

2 算例分析

例1：如图2 所示，考虑由两层等厚石墨-环氧树脂材料制成的正方形变刚度层合板[26]，尺寸为a=b=1.0 m，板厚h=0.1 m，铺层方式为[＜45°/90°＞/＜45°/0°＞]，四边固支，上表面施加法向载荷Pz=-10 kN/m2。材料属性如下：E1=137.9 GPa，E2=8.96 GPa，E3=8.96 GPa；G12=7.1 GPa，G13=7.1 GPa，G23=6.21 GPa；µ12=0.3，µ13=0.3，µ23=0.49。

图2 变刚度层合板模型及坐标系 /mFig.2 Variable stiffness laminate and coordinate system

以变刚度层合板的中心为原点，建立如图2所示坐标系。纤维路径函数为：

式中：β1和β0分别为层合板边界处(x=-0.5a)和中心处的曲线切向与x轴的夹角。单层板的路径曲线标记为＜β1/β0＞。

为了评估本文非协调广义部分混合元法的适用性以及准确性，对比ABAQUS 的参考解，在ABAQUS 中采用二十节点六面体单元网格计算分析，网格划分为80×80×20，共128 000 个单元。

本文的方法采用八节点六面体单元，网格划分为32×32×20，共20 480 个单元。分析沿板厚方向点Q(-0.25a, -0.25b)(参见图3(a))处的位移和应力。从图4 可以看出，本文方法得到的沿板厚方向的位移uy和uz结果与ABAQUS 结果曲线吻合较好，尽管uz结果存在差异，但是最大误差在1%以内。由图5 可以看出，尽管本文方法的低阶单元数量小于ABAQUS 中20 节点位移单元数量。但是，本文方法得到的应力曲线与ABAQUS得到的应力曲线基本吻合。例如，在界面处的σxy结果与ABAQUS 的相对误差仅为1.9700%，其它 应 力 σy、 σz和 σxz包 括 本 文 未 给 出 应 力 σx和σyz与ABAQUS 的相对误差均小于1.9700%。另一方面，由于本文考虑了应力边界条件，因此在层合板上下表面处的 σxz结果为0(如图5(d)所示)，符合真实情况。

图3 变刚度层合板的铺层方式Fig.3 Lamination method of the variable stiffness laminate

图4 四边固支变刚度层合板位移沿厚度方向的分布Fig.4 Distribution of displacement along the thickness of variable stiffness plate with sides clamped

图5 四边固支变刚度层合板应力沿厚度方向的分布Fig.5 Distribution of stresses along the thickness of variable stiffness with sides clamped

例2：考虑TORNABENE 和BACCIOCCHI[25]提出的三层变刚度夹芯板的弯曲问题，尺寸为a=b=1.5 m、h=0.1 m，四边固支，上表面处(z=h/2)受均布载荷qz=-10 kPa 作用。三层夹芯板由下到上每层厚度依次为0.02 m、0.06 m、0.02 m。夹芯板第一、第三层由曲线玻璃纤维环氧树脂制成。第二层的芯为泡沫橡胶，视为各向同性材料。相关材料属性见表1。以夹芯板的中心为原点，板长方向为x轴方向建立坐标系。第一、第三层变刚度板的纤维路径函数分别为βlayer1、βlayer2，表达式如下：

表1 玻璃纤维环氧树脂和泡沫橡胶的材料属性Table 1 Material parameters of glass fiber reinforced epoxy resin and foam rubber

式中：β 为x=-0.5a处纤维切向与x轴的夹角，称为初始角度，由表2 给定； θ′的表达式为：

表2 铺层方式Table 2 Laying methods

考虑如表2 所示的四种不同的纤维初始角度(β1)的铺层方式(a、b、c、d)。其中β1、β0分别为x=-0.5a、x=0 处的纤维角度。四个案例的网格划分方式均为28×28×30。选择P(-0.25a,-0.25b)点沿厚度方向节点作为研究对象，本文方法的位移结果与文献[15]的位移结果曲线列于图6，与文献[25]的应力结果曲线列于图7。

图6 变刚度夹芯板沿厚度方向的ux 和uy 的分布Fig.6 ux and uy along the thickness of the variable stiffness sandwich plate

图7 变刚度夹芯板沿厚度方向的应力分布Fig.7 Stresses distribution along the thickness of variable stiffness sandwich plate

从图6 可以清楚地看出，本文方法的位移曲线与文献[15]位移曲线变化趋势吻合较好。同时可以看出，随着纤维铺层角度不断增加，即由案例a～案例d，ux曲线逐渐左移，uy曲线逐渐右移。基于案例d 的纤维铺设可解释位移曲线如此变化的原因：案例d 的纤维铺设导致夹芯板第一层在x方向和第三层在y方向的刚度显著降低。故导致此时的夹芯板的ux在第一层最小，uy第三层最大。

图7 为本文方法沿夹芯板厚度方向的应力解与文献[25]结果的对比。由图7(a)和图7(b)可以清楚地看出，本文方法得到的面内应力结果与文献[25]的结果吻合较好。从图7(c)和图7(d)可以推断，本文方法与文献[25]的横向剪切应力在四种铺层方式下都有较好吻合。然而在案例c 和案例d 的铺层方式下，文献[25]的 σxz(图7(d))应力曲线在界面处出现突变，这是不真实的情况，而本文方法在界面处的应力曲线无异常突变。

文献[25]以板壳理论为基础，对该变刚度层合板进行分析。鉴于分析对象属于厚板类结构，因此，本文从最基础的线弹性理论出发，采用非协调广义部分混合三维实体元仿真该变刚度层合板的静力学行为，可以避免板壳理论中的诸多假设。

3 结论

考虑到解析法、板壳理论以及位移有限元法的局限性和变刚度层合板的特点，本文在常规的非协调广义部分混合元理论的基础上，为变刚度层合板的静力学分析提出一种修正的非协调广义部分混合元模型。主要结论如下：

(1) 在四边固支和均布载荷作用下，对两层的变刚度层合板进行分析，数值结果表明，相对于ABAQUS 高阶单元、高网格密度的模型，本文的修正的非协调广义部分混合元模型在网格比较稀疏的情况下表现出了良好的适用性和数值结果精度高等特点。

(2) 有关四边固支的三层变刚度夹芯板的分析表明，平面内的位移受纤维角度的变化影响明显。主要原因是单层板内纤维角度变化可显著影响该层的刚度大小。例如，单层板内纤维角度越接近90 度，相应的其沿x方向的刚度越小，则x方向的位移ux越大。