各向异性修正偶应力Mindlin层合板的有限元热稳定性分析

张大千，王云鹏，王玺鉴

(沈阳航空航天大学 辽宁省飞行器复合材料结构分析与仿真重点实验室，沈阳 110136)

多数实验[1-3]表明，当结构尺度进入微米量级时，材料刚度与强度等力学性能较之传统理论明显增大，该现象被称为材料的尺寸效应。由于经典连续体力学无法解释此现象，因此，学者们在传统力学的基础上提出了偶应力理论和应变梯度理论。

1963年，Mindlin[4]提出了经典偶应力模型，该模型中应变张量对称，曲率张量不对称，并引入一个细观材料尺度参数。由于其曲率张量的非对称性，故将经典偶应力理论应用于实际工程中具有一定困难。

2002年，Yang等[5]提出修正偶应力理论，通过重新定义曲率张量，使应变张量与应力张量对称。基于此理论：2006年，Park 和Gao[6]建立了修正偶应力Bernoulli-Euler梁模型；2008年，Ma等[7]建立了Timoshenko梁模型；2009年，Tsiatas[8]建立了Kirchhoff板模型；2011年，Ma等[9]基于修正偶应力建立了Mindlin板模型；2012年，Alireza Nateghi等[10]研究了功能梯度梁的热尺度效应；同年，Asari等[11]分析了功能梯度Timoshenko梁的热后屈曲尺度效应。

由于具有优越的力学性能，复合材料层合板在各领域均得到了广泛的应用。尺寸效应不仅存在于金属材料，在复合材料中也同样存在。由于在上述理论的本构关系中，偶应力弯矩、曲率和应变关系均为各向同性，故这些理论并不能应用于复合材料层合板的研究。2012年，陈万吉[12]依据上述理论，提出了适用于各向异性材料的新修正偶应力模型，并推导出了新的本构方程关系。新理论中的曲率张量虽不再对称，但可等价于各向同性的曲率张量，以此为基础，陈万吉[13]对复合材料层合板的弯曲和自由振动问题进行了分析。2013年，陈万吉等[14]建立了Reddy层合梁模型。2014年，陈万吉等[15]基于修正偶应力理论，通过有限元法首次分析了偶应力层合板的尺寸效应。2015年，M.Azhari等[16]利用修正偶应力理论和样条有限元方法分析了对功能梯度板的稳定性。

综上所述，涉及复合材料层合板微细观条件下热尺度效应问题的研究较少，本文基于各向异性修正偶应力理论，建立微尺度下Mindlin层合板的热稳定模型，并应用有限元法构造三节点三角形单元，建立有限元热稳定模型，用上述两种模型，对微尺度下Mindlin层合板进行了热稳定性分析，并研究了Mindlin层合板的尺寸效应。

1 各向异性修正偶应力Mindlin层合板热稳定性模型

1.1 各向异性修正偶应力Mindlin层合板位移场

在整体坐标系下，将板上任一点沿x、y和z方向的位移分别表示为u、v和w，将绕x、y和z轴的转动位移分别表示为ωx、ωy和ωz，故Mindlin层合板的初始位移模式表达式为

(1)

转动位移表达式为

(2)

1.2 各向异性修正偶应力Mindlin层合板应变和曲率

在考虑温度变量ΔT的情况下，根据位移应变关系，可推导出层合板应变公式如式(3)所示

(3)

式(3)中，αx、αy分别表示沿x和y轴方向的线性热膨胀系数，层合板曲率公式为

(4)

1.3 各向异性修正偶应力Mindlin层合板本构方程

根据应力—应变关系，可构建整体坐标系下的本构方程如式(5)所示

σk=Qkεk

(5)

(6)

Qk=TkTCkTk

(7)

其中

(8)

lkb和lkm分别为绕纤维和基体方向转动的细观参数。

转换矩阵Tk表示为

(9)

(10)

式(10)中，m=cosΦk，n=sinΦk，Φk为每层的铺设角。

1.4 各向异性修正偶应力Mindlin层合板平衡方程

基于各向异性修正偶应力Mindlin层合板的虚功原理，可推导出平衡方程，虚功原理可表示为

δU-δW=0

(11)

内力虚功δU与外力虚功δW表达式分别为

(12)

(13)

将公式(3)、(4)、(6)、(12)和(13)带入公式(11)，利用变分法，可最终推导出平衡方程

(14)

相对于板的横向位移w，膜向位移u0、v0是小量，即方程中可忽略u0、v0，并假设

fu0=fv0=fcx=fcy=0

根据Navier解法，得到满足全部边界条件(四边简支层合板)的位移函数

(15)

联立方程(14)和(15)，可得受温度载荷和面内轴向载荷共同作用的平衡方程如下

(16)

其中

(17)

(18)

(19)

因此，由公式(16)可求解当温度变化量为ΔT时面内轴向载荷N1和温度载荷N2共同作用下的Mindlin型层合板的失稳临界载荷。

2 Mindlin层合板热稳定有限元模型

根据高阶剪切变形理论中的位移和应变的关系，可知层合板的横向剪切应变为二次函数，即理论中出现了w的一阶和二阶导数，故须构造满足C0连续和C1弱连续的单元横向位移函数，从而保证单元具有一定的收敛性。

2.1 满足C0连续条件的横向位移函数w0

本文选用九参数非协调单元BCIZ[17]构造满足C0连续条件的横向位移函数w0，其公式为

w0=Fq

(20)

F=[F1Fx1Fy1F2Fx2Fy2F3Fx3Fy3]

(21)

(22)

(23)

2.2 满足C1弱连续条件的横向位移函数w*

本文选用修正九参数三角形薄板弯曲单元RT9[18]构造满足C1连续条件的横向位移函数w*，其公式为

w*=w0+p(Bc-B0)q=F*q

(24)

将公式(21)带入(25)，可得

F*=F+p(Bc-B0)

(25)

矩阵Bc可表示为：

Bc=Bca+0.25(Bca-Bcb)

(26)

(27)

(28)

上述公式中

(29)

(30)

xi、yi是节点坐标，Bca2和Bca3可通过下标轮换求得。

修正矩阵Bcb可表示为

(31)

(32)

同理可通过下标轮换求得Bcb2和Bcb3。

2.3 三角形单元的位移函数

根据节点变量和形函数，可将三角形单元位移函数表示为

(33)

2.4 三角形单元的应变、单元刚度矩阵及几何刚度矩阵

根据应变与位移之间的关系，可将应变表示为

ε=[∂]u*=[∂]Nqe=Bqe

(34)

其中

(35)

(36)

(37)

单元刚度矩阵表示如下

(38)

(39)

板单元的几何刚度矩阵表示为

(40)

式(40)中

(41)

(42)

(43)

Nxx、Nyy分别表示来自x和y方向的面内轴向载荷，NxT、NyT分别表示来自x和y方向的温度载荷，Nxy代表的是剪切应力。

求解临界失稳载荷方程为

(44)

其中，求解λ即为求解临界失稳载荷，Г表示屈曲形函数的矢量模式。

3 数值算例

利用各向异性修正偶应力理论，求解Mindlin层合板算例的解析解，并与利用有限元法的三角形单元求出的数值解对比，验证有限元法的准确性。算例的边界条件见表1所示。

材料参数[17]：E2=6.98 GPa，E1=25E2，G12=0.5E2，G22=0.2E2，v12=v22=0.25，v21=E2v12/E1=0.01，下标1和2分别表示纤维和基体的方向。

层合板的尺寸：板厚h=20 μm，板长L=200 μm。

表1 算例的边界条件

3.1 算例1：验证有限元法的准确性

为了验证提出的三角形单元的准确性，将四边简支方板以铺设角[0°/90°/0°]的方式正交铺设，网格划分参见图 1所示，网格划分方式与失稳临界载荷关系见图2。

图1 层合板网格划分

图2 网格划分与失稳临界载荷关系

由图2可知，当温度变化量ΔT分别为0 ℃、20 ℃和40 ℃时，随着网格不断细化，有限元法计算得到的临界载荷收敛速度逐渐变缓，在网格数达到一定值后，计算结果趋于稳定且与解析解数值相近。为保证有限元法的计算精度，在后续的计算过程中均采用20×20的方式进行有限元网格划分。

表2和表3分别针对单向轴压和双向轴压作用，在温度变化量ΔT=40°的情况下，对有限元法和解析法计算的每阶失稳临界载荷进行了比较。由表2、3中的结果可知，在考虑温度载荷作用的情况下，通过有限元法计算得到的失稳临界载荷与解析解数值相近，二者之间误差较小。利用此算例可证明：本文选取的三节点三角形有限元单元较为合理，在计算失稳临界载荷时具有一定的准确性，故建立的几何模型和采用的计算方式在实际工程应用中具有使用价值。

表2 单向轴压作用下四边简支板的失稳临界载荷(L/h=10，l=0)

表3 双向轴压作用下四边简支板的失稳临界载荷(L/h=10，l=0)

3.2 算例2：计算考虑ΔT时的失稳临界载荷

本算例考虑温度变化量ΔT对临界失稳载荷的影响，其中ΔT分别取0°、20°、40°、60°和80°，四边简支板铺设角为[0°/90°/0°]，位移模数取值：a=1，b=1。通过改变材料尺度参数的数值，对各向异性修正偶应力Mindlin层合板稳定性的尺寸效应进行分析，分别计算单向轴压和双向轴压作用下有限元法与解析法的失稳临界载荷。图3和图4中实线表示通过各向异性修正偶应力理论计算得到的解析解，虚线表示本文建立的有限元模型计算得到的有限元解。结果表明：在上述两图中，随着材料尺度参数的增大，失稳临界载荷越大，因此，该算例证明了尺寸效应的存在。且失稳临界载荷随着温度变化量ΔT的增加而增大，即材料的尺寸效应更明显。考虑温度变量时，产生的温度载荷对失稳临界载荷具有一定影响，且不可忽略。

图3 四边简支板单向轴压临界失稳载荷随l/h变化

3.3 算例3：验证有限元法的优势

为验证本文提出三角形单元的实用性，在考虑温度变化量ΔT的前提下，通过改变材料尺度参数的数值，对不同边界条件、铺设角度和铺设层数的层合板单向轴压临界失稳载荷分别进行求解。其中，边界条件由表4给出，算例的计算结果如表5所示。

图4 四边简支板双向轴压临界失稳载荷随l/h变化

表4 算例3的边界条件

表5 不同计算条件下的单向轴压临界失稳载荷(L/h=10)

由表5的数值结果可知：边界条件、层合板铺设角度及铺设层数等参数一定时，层合板单向轴压临界失稳载荷随尺度参数的增大而增大，证明有限元法可捕捉尺寸效应，且温度变化量越大，尺寸效应越显著。相较于计算解析解时使用的Navier解法对边界条件的要求十分苛刻(仅适用于四边简支的边界条件)，本文提出的三角形单元在保证计算效率与准确性的同时，适用性更强，能够计算不同的边界条件、铺设角度与铺设层数的偶应力层合板问题。

4 结论

本文首先建立了考虑一阶剪切变形的温度载荷与机械载荷共同作用下的Mindlin层合板的热稳定模型，然后提出了一种新的三角形单元，该单元满足C0连续和C1弱连续的条件，建立了有限元模型，并计算了不同条件下的算例。结果表明：

(1)有限元法计算得到的失稳临界载荷与解析解相近，证明了本文提出的三角形单元具有较高的准确性。

(2)依据本文提出的三角形单元建立的有限元模型可验证尺寸效应的存在。

(3)当温度变化量ΔT存在时，温度变化量一定，失稳临界载荷随着尺度参数的增加而增强，尺度参数不变，失稳临界载荷随温度变化量的增加而增大，故温度载荷对失稳临界载荷产生的影响不可忽略。

此三角形单元的提出，证明本文建立的有限元模型可用来分析各向异性Mindlin层合板的热稳定问题。