复合材料层合板热弹性问题的新半解析法

卢 翔，屈鹤鹏，刘彦波

（中国民航大学航空工程学院，天津 300300）

随着复合材料在航空航天、汽车工业、交通运输、化工纺织等领域得到越来越多的应用，复合材料层合板的热弹性问题逐渐成为材料界和力学界的研究重点。研究复合材料层合板的理论主要有以下几种：经典层合板理论、一阶剪切变形理论、高阶剪切变形理论、高阶层合板理论、分层理论以及三维弹性理论[1-3]。上述层合板理论没有考虑复合材料承载能力远远大于普通材料的特点，在研究复合材料层合板问题时均存在不足，不能兼顾精确性和方便性。杨正林等[4]以Hamilton 正则方程理论为基础，给出了层合板热弹性问题体的半解析法；邹贵平[5]在一阶层合板理论的基础上，采用共轭辛正交的方法给出了层合板热应力分析的精确解法；随后，相关学者通过增维积分法建立了热弹性体的齐次状态方程，并给出了相应的精确解法[6-9]。传统的半解析法在求解层合板位移和应力时有一定的震荡性，且随着网格密度增大，高维指数运算需要耗费大量机时，不适应于大规模有限元问题。

以弹性力学最小势能原理和修正的H-R 变分原理为基础，将热弹性体的位移单元[10]与Hamilton 半解析混合单元[11]相结合，构建了新的半解析混合单元，提出了新的半解析法，结合位移法和Hamilton 正则方程半解析法的优点，大大简化了问题的计算。通过实例计算，研究了梯度温度载荷作用下层合板的热弹性问题，通过与不同方法的计算结果[12]相对比，验证了新方法的可靠性。

1 基本理论推导

1.1 基本理论

弹性体在变温情况下的应变分量包括自由热膨胀引起的应变分量和热膨胀时弹性体内各部分之间相互约束引起的应变分量，其与热应力之间同样服从胡克定律。因此，参照弹性力学的基本方程，可得各向正交异性复合材料热弹性问题的本构方程为

其中：σ = [σxxσyyσzzσyzσxzσxy]T为应力变量向量；T 为温度载荷向量；ε=[εxxεyyεzzεyzεxzεxy]T为应变向量；S = C-1为热弹性体的柔度系数矩阵；α =[α1α2α30 0 0]T为线膨胀系数矩阵；

为刚度系数矩阵。

1.2 Hamilton 正则方程和最小势能原理

1.2.1 Hamilton 正则方程

根据Hellinger-Reissner 广义变分原理，层合板热弹性问题的泛函表达式为

其中：Su为应力边界；Sσ为位移边界；Q=[u v w]T为弹性体的位移向量为边界面上已知的位移向量；p=[pxxpyypzz]T为弹性体受到的表面力；为边界面上已知的应力向量。

根据弹性力学的Hamilton 正则方程理论，可得热弹性体修正后的H-R 变分原理为

其中：p = [ pxzpyzpzz]T为平面外应力向量；H 为Hamilton 函数，其表达式为

其中：G1和G2为偏微分算子矩阵；Φij为本构关系式（1）经转置后应力向量与应变向量之间的关系矩阵[6]。

按照文献[4]中的方法，分别以P 和Q 为变量对Hamilton 函数进行变分，可得热弹性体的Hamilton 正则方程为

其中：A，B，C 和D 分别为Hamilton 函数经变分后P和Q 的系数矩阵，即

1.2.2 最小势能原理

弹性力学中的传统有限元位移法是以最小势能原理[13]为基础的，不考虑体积力，则热弹性体最小势能原理的泛函可表示为

热弹性体的几何方程可表示为

其中，X 为热弹性体应变向量ε 与位移量Q 之间的关系矩阵，即

将式（8）代入式（7），可得用位移变量表示的最小势能原理为

2 半解析混合单元模型

将文献[11]中的Hamilton 等参元法运用到求解层合板的热弹性问题中。对于任意形状的四边形单元，以平面四节点Hamilton 等参元为例，将插值后的变量[11]代入式（9）即可得到经Hamilton 等参元离散后任意单元的最小势能原理的泛函表达式，即

由于δΠP=0，因此以Qe为变量，对式（10）进行变分，可以得到

其中

按照同样方法对热弹性体的Hamilton 正则方程式（6）进行Hamilton 等参元离散，即可得到离散后的Hamilton 正则方程为

其中：N 为插值函数矩阵；J 为雅各比矩阵式[11]，选取式（12）中的第1 个方程

将其与式（11）合并，可建立热弹性体经Hamilton等参元离散后的半解析混合单元模型（不计体积力），即

3 层合板热弹性问题混合方程

层合板的任一子层均可用上述方法进行离散并建立半解析混合单元。经过单元组装，可以得到任一子层的混合状态方程。具体对式（15）求和后可得

其中：m 为层合板第m 子层；Am、Bm和Km分别为经单元组装后的系数矩阵；f1m、Γ1m和Γ2m分别为f1e、Γ1e和Γ2e进行单元求和后的结果，即

在具体求解时，首先对式（16）第2 个方程进行位移求解，然后将所得位移结果一一代入第1 个方程进行平面内应力的求解。因此，式（16）中第1 个方程就变成更为简洁的一阶微分方程，即

其中：Φm（z）=BmQ（z）m+f（z）1m

式（17）的解为

其中

对于n 层的层合板，根据层间的连续关系，则有

其中：m = 2，3，…，n；hm为层合板第m 子层的厚度；Pm（0）为第m 层下表面的平面外应力向量；Pm（hm）为第m 层上表面的平面外应力向量。

引入传递矩阵并结合式（18）～式（19），可以得到层合板上下表面应力之间的关系表达式。然后根据层合板上下表面的应力位移初始值即可求得层合板任意位置处的应力和位移。

以上求解方法实际上就是先求位移，然后求平面外应力的微分方程。很明显，对比式（6）和式（17）的系数矩阵可以发现，后者的指数维数是前者的1/4。由于传统位移法在求解位移时是准确的，因此从理论上讲，新半解析法在保证计算精度的同时，可以有效提高运算速度，减少运算时间。

4 实例分析

选取1 块等厚度的矩形薄板，如图1 所示（a，b 远远大于厚度h），温度载荷T 为y 的函数。无量纲化后的几何性质：a=100，b=20；材料特性：弹性模量E=210，泊松比v=0.3；应力温度系数：α11=0.000 1，T0=100；板的上下表面和侧面均为应力自由表面。分别考虑两种温度载荷条件下矩形薄板的热弹性问题：①T（y）=T0（1-y2/b2）；②T（y）=T0（1-y3/b3）。

图1 自由边界条件下的矩形薄板Fig.1 Rectangular plate with free edges

由文献[13]可知矩形薄板内热应力的精确解为

离散过程中，由于网格的疏密程度对计算结果影响较大，因此只选用32 × 12 网格密度下新半解析解的计算结果进行误差分析，即

A（0，0）和B（0，-20）两点在不同温度载荷条件下的应力计算结果如表1～表4 所示。

表1 和表2 分别给出了施加温度载荷1 时不同网格密度条件下A（0，0）点和B（0，-20）点热应力σxx的计算结果。图2 和图3 分别是温度载荷条件1 作用下A（0，0）点和B（0，-20）点的热应力σxx随网格密度的变化。由图2 和图3 可以直观地看出，相对于ABAQUS解，新半解析解的计算结果受网格疏密程度的影响更小。通过横向对比可以发现，当网格疏密程度一样时，新半解析解比ABAQUS 解更接近精确解，具有更高的计算精度。即使在网格密度较小的情况下，新半解析解的应力计算结果也和精确解十分接近。

表1 温度载荷条件1 作用下点A（0，0）处热应力σxx的3 种解Tab.1 σxxat Point A（0，0）under Thermal Case One

表2 温度载荷条件1 作用下点B（0，-20）处热应力σxx的3 种解Tab.2 σxxat Point B（0，-20）under Thermal Case One

表3 温度载荷条件2 作用下点A（0，0）处热应力σxx的3 种解Tab.3 σxxat Point A（0，0）under Thermal Case Two

表4 温度载荷条件2 作用下点B（0，-20）处热应力σxx的3 种解Tab.4 σxxat Point B（0，-20）under Thermal Case Two

表3 和表4 分别给出了施加温度载荷2 时不同网格密度条件下A（0，0）点和B（0，-20）点热应力σxx的计算结果。由表3 可以看出，在温度载荷条件2 作用下，矩形薄板中心处的热应力σxx=0，无论是新半解析解还是ABAQUS 有限元解，都很好地符合了这一结果。由表4 可以看出，在温度载荷条件2 作用下，随着网格密度增大，两者的计算结果均越来越接近精确解。不同的是，与ABAQUS 有限元解相比，新半解析解的计算结果误差更小、精确度更高，并在较小的网格密度下可保证计算的精确性。网格密度决定了计算机的工作量，因此新半解析解在保证较高计算精度的同时，可有效缩短计算时间。同时对计算机的运行内存和运行速度要求不高，适用于大规模有限元计算。

图2 A（0，0）点σxx随网格密度的变化Fig.2 σxxvariation at Point A along with mesh density change

图3 B（0，-20）点σxx随网格密度的变化Fig.3 σxxvariation at Point B along with mesh density change

在合适的网格密度下，ABAQUS 解也可得到较为精确的结果。由于研究温度梯度对层合板层间热应力影响的文献较少，因此在没有精确解做参考的情况下，也可以将新半解析解和ABAQUS 解进行对比分析。

5 结语

新半解析法承袭了传统位移法和Hamilton 半解析法的优点，在保证计算精度的同时，可以大大降低计算量，有效减少计算时间，相对于原有的半解析混合法和传统有限元法，更适应于大规模有限元问题的分析。分别采用新半解析法和ABAQUS 有限元仿真计算不同温度载荷条件下矩形薄板的热弹性参数，并与精确解进行对比，验证了该方法的优越性。