基于等几何有限元法的功能梯度微板热力耦合屈曲预测

邓 阳， 尹硕辉， 赵子衡

(湘潭大学 机械工程学院，湘潭 411105)

1 引 言

功能梯度材料是一种通过连续改变组分材料分布，使自身整体性能呈现梯度变化的非均质复合材料。该材料使结构具有高耐热性、高强度和高刚度等优点，因此在航空航天、生物医疗和机械等领域得到了广泛应用。

随着现代微机电技术的发展，功能梯度材料已经应用于微观结构[1]和系统[2]中。实验表明，当结构的尺寸处于微米量级(1 μm～100 μm)时，结构将表现出比宏观尺度下更高的刚度，即尺度效应现象[3，4]。传统连续介质理论无法解释这种现象，因此，人们发展了非局部弹性理论[5]和偶应力理论[6]等高阶连续介质理论。其中，Yang等[7]基于传统的偶应力理论，提出了一种适用于各向同性材料且仅有一个材料尺度参数的修正偶应力理论。Tsiatas[8]基于该理论建立了Kirchhoff微板模型，该模型能用于分析任意几何形状和复杂边界条件的微板静力学问题；Park等[9]分析了尺度效应对微悬臂梁弯曲刚度的影响。陈万吉等[10]将该理论推广至各向异性材料，提出了新修正偶应力理论。基于该理论，贺丹等[11,12]建立了各向异性功能梯度微板模型。

目前，国内外学者仅考虑了单一物理场中微结构的力学性能。然而，在一些特定的应用环境中，微结构会受到多种物理场作用。因此，亟需研究多场作用下微结构的力学问题。Tan等[13]在热力耦合环境下，研究了尺度效应对层合悬臂微梁弯曲挠度的影响；Jia等[14]同时考虑了温度、静电力和分子间作用力的影响，建立了功能梯度微梁的非线性屈曲模型。

由于偶应力中存在高阶项，采用常规有限元法构造高阶连续单元较复杂。等几何有限元法是由Hughes等[15]提出的一种新型高阶数值计算方法。该方法是基于有限元法等参单元的思想，直接采用CAD几何建模的非均匀有理B样条(NURBS)基函数作为有限元法中的形函数，从而实现CAD与CAE的无缝连接。与常规有限元法相比，等几何有限元法能够轻松构造高阶连续函数，精确地表示几何形状。该数值方法在求解考虑偶应力的微尺度模型中已经得到了成功应用[16，17]。

本文结合修正偶应力理论及等几何有限元法，建立了Kirchhoff功能梯度微板热力耦合屈曲模型。研究了热力耦合环境下，材料尺度参数、板的边长比及梯度指数对功能梯度微板稳定性的影响。

2 基本方程

2.1 功能梯度材料模型

组分材料的连续改变导致功能梯度材料的物性参数呈梯度变化。为了描述物性参数的这种变化，已经提出了一些数学模型。本文采用的幂函数模型[18]表示为

(1)

式中 下标1，2表示不同材料；V1为材料1的体积分数；h为板的厚度；n为材料的梯度指数；p(z)表示功能梯度材料的杨氏模量E、热膨胀系数α和泊松比ν等物性参数。本文泊松比ν假设为常数。

2.2 Kirchhoff板理论

根据Kirchhoff板理论，位移场表示为

(2)

2.3 修正偶应力理论

根据修正偶应力理论，各向同性材料的线弹性本构关系表示为

(3)

(4)

式中 杨氏模量E(z)沿着板的厚度方向变化；ν是泊松比。

(5)

式中ui是位移矢量,转动矢量θi表示为

(6)

其中ei j k是置换符号。

将式(2，6)代入式(5)，并考虑热载荷ΔT引起的变形，则应变张量和旋转梯度张量表示为

ε=εp+zεb-εt

(7)

(8)

式中εp为中面的伸缩剪切应变，εb为中面的弯扭应变，εt为温度引起的应变，α(z)为沿着板的厚度方向变化的热膨胀系数。

3 等几何有限元法

3.1 节点矢量和B样条基函数

节点矢量U={ξ1,ξ2,…,ξn + p + 1}是由一维参数空间中一系列单调递增的坐标值组成，其中ξi表示第i个节点的坐标值(i=1，2，…，n+p+1)；p是B样条基函数的阶次；n是B样条基函数的数量(与控制点数量相同)。根据节点矢量U，ξ方向的p阶B样条基函数Ni,p(ξ)表示为

(9)

(10)

3.2 NURBS基函数

(11)

3.3 等几何有限元位移模式

等几何有限元法是在有限元法的框架下，采用CAD中构造几何形状的非均匀有理B样条(NURBS)基函数作为有限元法中的形函数，位移场近似表示为

(12)

式中uI=[uIvIwI]T表示单元内第I个控制节点的自由度，NP=(p+1)(q+1)为单元控制点总数，RI是NURBS基函数。

将式(12)代入式(7，8)得

(13,14，15)

式中Bp为伸缩剪切应变矩阵，Bb为弯扭应变矩阵，Bm为旋转梯度矩阵。

式中RI,x和RI,y分别为NURBS基函数对坐标x和y的一阶偏导数，RI,x x，RI,y y和RI,x y分别为NURBS基函数对坐标x和y的二阶偏导数。

4 热力耦合屈曲等几何有限元方程

基于修正偶应力理论，应变能表示为

(16)

将式(3,7,8)代入式(16)得

(17)

外力做功表示为

(18)

由虚功原理：

δΠ=δU-δW=0

(19)

将式(17,18)代入式(19)，得到功能梯度微板热力耦合屈曲等几何有限元方程为

(K+Km-Nc rKG N-ΔTKG T)d=0

(20)

式中K为传统刚度矩阵，Km为与尺度效应相关的刚度矩阵，Nc r为临界屈曲力，ΔT为临界屈曲热载荷，KG T与KG N为几何刚度矩阵。具体表示为

5 数值算例

考虑如图1所示的微板，长度a,宽度b，受到均布热载荷及双向均布压力的作用。为了方便表示微板的边界条件，用C表示板边固支约束，S表示板边简支约束。

5.1 收敛性及正确性验证

为了验证模型的收敛性，本文模型退化为功能梯度微板力屈曲模型。其中，板的长和宽取为a=b=10h，板的厚度h=17.6 μm，材料尺度参数l=h，梯度指数n=1，弹性模量Ec=14.4 GPa，Em=1.44 GPa，泊松比ν=0.38，临界屈曲力的归一化公式为N*=Nc rb2/Emh3。图2为单元控制点数量和NURBS基函数阶次对模型预测结果误差的影响。其中，Error=(N*-N)/N，N*为目前模型预测的结果，N为解析解。可以看出，当NURBS基函数阶次p=1时，本文模型预测结果的相对误差很大；当p≥ 2时，随着单元控制点数量的增加，本文模型十分迅速地收敛到一个较小的相对误差。通过权衡计算精度与效率，本文后续算例的控制点数将采用16×16，NURBS基函数阶次p=q=3。从两个不同角度验证本文模型正确性。通过继续使用上述参数和归一化公式，将本文模型预测的临界屈曲力与文献[19，20]的结果对比，由表1可知，本文模型的预测结果与文献结果吻合 较好。

图1 受双向均布压力的功能梯度微板

当临界屈曲力Nc r=0时，本文模型退化为功能梯度微板热屈曲模型。其中，板的长和宽取为a=b=100h，板厚度h=17.6 μm，材料尺度参数l=h，梯度指数n=1，弹性模量Ec=380 GPa，Em=70 GPa，泊松比ν=0.3，热膨胀系数αc=7.4×10-6/℃，αm=23×10-6/℃。将本文模型预测的临界屈曲热载荷与文献[20]结果对比，由 表2 可知，本文模型的预测结果与文献结果吻合良好，证明了本文模型的正确性。

图2 控制点数量及基函数阶次对预测结果误差的影响

表1 本文模型预测的临界屈曲力与文献[19,20]结果的对比

5.2 材料尺度参数对功能梯(Al/Al2O3)微板稳定性的影响

研究热力耦合环境中，材料尺度参数对功能梯度(Al/Al2O3)微板稳定性的影响。其中，微板的长和宽为a=b=20h，板厚h=17.6 μm，梯度指数n=1，弹性模量Ec=380 GPa，Em=70 GPa，泊松比ν=0.3，热膨胀系数αc=7.4×10-6/℃，αm=23×10-6/℃。临界屈曲力Nc r和临界屈曲热载荷ΔT的归一化公式为N*=Nc rb/Emh2，T*=αmΔTb/h。图3显示了材料尺度参数对功能梯度微板稳定性的影响。

表2 本文模型预测的临界屈曲热载荷与文献[20]结果的对比

图3 材料尺度参数对功能梯度(Al/Al2O3)微板稳定性的影

可以看出，在相同边界条件下，本文模型预测的临界屈曲载荷(l/h≠0)大于宏观理论(l/h=0)的预测结果，此时微板表现出比宏观状态下更高的刚度；随着板厚度h增加，尺度效应逐渐减弱；当板厚度h远大于材料尺度参数l时，尺度效应消失，此时本文模型将退化为宏观理论模型(l/h=0)；当l/h不变时，CCCC功能梯度板的临界屈曲载荷大约是SSSS功能梯度板的临界屈曲载荷的两至三倍；功能梯度微板的临界屈曲热载荷随着临界屈曲力的增加呈线性递减。以图3(a)中l/h=1为例，根据最小二乘法拟合，这种线性关系可以表示为T*=-0.3735N*+0.5468。

5.3 边长比对功能梯度(Al/Al2O3)微板稳定性的影响

图4为在热力耦合环境中，边长比对功能梯度(Al/Al2O3)微板稳定性的影响。其中l/h=1，板的几何尺寸、材料参数以及归一化处理与5.2节相同。可以看出，在相同的边界条件下，功能梯度微板的临界屈曲载荷随着边长比的增加而逐渐增大；从图4(b)可以看出，当b/a≤4时，边长比的增加对CCCC板的稳定性有显著提升；当b/a≥4时，边长比的增加对板稳定性的影响较小；当边长比不变时，CCCC功能梯度微板的临界屈曲载荷总高于SSSS功能梯度微板的临界屈曲载荷；功能梯度微板的临界屈曲热载荷随着临界屈曲力的增加呈线性递减。以图4(a)中b/a=5为例，根据最小二乘法拟合，这种线性关系可以表示为T*=-0.3736N*+7.0318。

图4 边长比对功能梯度(Al/Al2O3)微板稳定性的影响

5.4 梯度指数对功能梯度(Al/Al2O3)微板稳定性的影响

考虑梯度指数对功能梯度(Al/Al2O3)微板稳定性的影响，其中l/h=1，板的几何尺寸、材料参数以及归一化处理与5.2节相同。可以看出，当n=0时板的稳定性最好，此时微板完全由陶瓷材料(Al2O3)构成，表明陶瓷材料具有更好的耐热性和抗屈曲能力；当梯度指数n=1时，板的稳定性显著下降；随着梯度指数升高，金属材料(Al)所占的比例不断增加，板的稳定性不断降低；当固定梯度指数n时，CCCC功能梯度微板的临界屈曲载荷总高于SSSS功能梯度微板的临界屈曲载荷；功能梯度微板的临界屈曲热载荷随着临界屈曲力的增加而线性减少。以图5(a)中n=1为例，根据最小二乘法拟合，这种线性关系可以表示为T*=-0.3735N*+0.5468。

图5 梯度指数对功能梯度(Al/Al2O3)微板稳定性的影响

6 结 论

基于修正偶应力理论，建立了功能梯度Kirchhoff微板热力耦合屈曲模型，并采用等几何有限元法求解。算例结果表明，本文模型预测的临界屈曲载荷总大于宏观理论预测的结果，说明本文模型捕捉到了尺度效应现象；随着板厚度减小，尺度效应现象更加显著，从而使功能梯度(Al/Al2O3)板表现出更高的刚度，提升了微板的稳定性；随着临界屈曲力的增加，临界屈曲热载荷呈线性递减；功能梯度(Al/Al2O3)微板的临界屈曲载荷随着板边长比的增加而逐渐增加。此外，梯度指数和边长比也对微板的稳定性产生一定影响。