基于等参梯度单元的功能梯度材料结构分析

张 龙，邱荣凯，程 俊，刘秉斌

（中国空气动力研究与发展中心设备设计及测试技术研究所，绵阳 621000）

功能梯度材料可设计用于抵抗高温度梯度、减小热应力以及提高界面连接强度［1⁃2］，在航空航天热结构领域得到了广泛的应用研究，如高超飞行器的梯度密度一体化防隔热材料［3］、冲压发动机梯度隔热层［4］和航空发动机涡轮叶片热障涂层［5］等。近年来，随着增材制造技术、化学气相沉积等先进制造工艺不断发展，使得功能梯度材料的可设计性越来越强、制备产业越来越成熟。与此同时，功能梯度材料特殊的力学行为也给结构分析与设计带来了新的挑战［6］。

伏培林等［7］采用解析法，开展了功能梯度悬臂梁的自由振动分析。针对欧拉⁃伯努利梁、瑞利梁和铁木辛柯梁，采用分离变量法推导了功能梯度悬臂梁自由振动的解析解，进而给出了不同跨深比和材料性能梯度变化的悬臂梁前三阶固有频率。许琦等［8］提出了一种考虑横法向热应变的三参数Reddy 型高阶功能梯度梁理论，可提高热力响应的分析精度。李华东等［9］研究了四边简支具有功能梯度芯材的夹层板在分布载荷作用下的弯曲问题，基于Reissner 假设，根据功能梯度材料的本构方程得出了应力、位移及内力的表达式。高晟耀等［10］基于一阶剪切变形理论推导了中等厚度功能梯度球环结构公式，利用里兹法求解得到功能梯度球环结构固有频率，分析了一般边界条件下中等厚度功能梯度球环结构的自由振动特性。高英山等［11］基于一阶剪切变形假设和哈密顿原理开展了碳纳米管增强功能梯度复合板非线性建模与仿真，讨论了碳纳米管体积分数、分布方式、结构宽厚比和载荷对功能梯度复合板的影响。苏盛开等［12］采用经典欧拉梁理论和高阶三角剪切变形理论，研究了多孔功能梯度梁的热力耦合屈曲行为，讨论了材料非均匀参数、孔隙率和长细比等参数对屈曲临界温度的影响。

以上研究均是基于解析或者半解析方法，研究对象为简单的梁、板、壳等结构，受载形式也较为简单。当结构或载荷复杂度增加时，有限元等数值计算方法是强有力的分析工具。但是由于功能梯度材料宏观上的材料性质不均匀性，采用分层法［13］等常规有限元开展分析，势必需要采用极高的网格密度以表征材料梯度，导致计算量巨大。Santare等［14］和Kim 等［15］提出了等参梯度有限元法，其单元本身可以反映材料力学性能梯度。Burlayenko等［16］和Xiong 等［17］分别在此基础上提出了改进的梯度有限元法，可用于功能梯度材料的力⁃热耦合分析和磁⁃电⁃弹性耦合分析。黄立新等［18］采用功能梯度材料平板，研究指出等参梯度单元相比常规单元具有更高的计算精度，并且随着弹性模量梯度增大，梯度单元的优势越大。韩国凯等［19］基于梯度有限元法思路，利用材料用户子程序UMAT 定义材料梯度，开展了防隔热一体化复合材料的整体性能优化设计。

本文发展了平面四节点与八节点等参梯度单元，从收敛性、计算精度和计算效率等方面验证了梯度单元相对于常规单元的优势。最后，基于ABAQUS 平台开发等参梯度单元 UEL（User⁃defined element）子程序，开展了功能梯度材料开孔结构和悬臂梁的计算对比分析。

1 功能梯度材料模型

功能梯度材料由两种或两种以上材料混合制备而成，并且材料组分在空间上渐进均匀变化，不存在含量百分比阶跃的材料界面，消除了传统复合材料组分界面处力学性能突变这一缺陷，从而能更好地发挥复合材料的优势。如图1 所示，在x＜0区域内为材料A；在x＞L区域内为材料B；在0≤x≤L区域内为功能梯度材料，材料A 和材料B 的含量百分比渐进均匀变化。

图1 功能梯度材料示意图Fig.1 Sketch of functional graded materials

在x=0 边界处，材料性能与材料A 相同；在x=L边界处，材料性能与材料B 相同；在0＜x＜L区域内，根据组分含量百分比的变化程度不同，材料性能可为坐标x的线性函数或指数函数，即

式中：φ（x）、φA、φB分别为梯度材料、材料A、材料B的性能指标，可以是密度、弹性模量和泊松比等。

2 等参梯度单元

2.1 有限元基本理论方程

在有限元理论中，单元内任意一点处的位移可以表示为

式中：Ni为形函数；uei为编号为i的节点处位移；n为该单元的节点数量。

确定了单元位移后，利用几何方程和物理方程，可将单元的应力应变表示为

式中：B为应变矩阵；D为材料刚度矩阵。

根据虚功原理，可推导有限单元刚度矩阵方程如下［20］

式中：fe为单元节点载荷；ke为单元刚度矩阵，可由式（7）得到

式中Ωe表示单元的控制域。将以上单元水平内的推导方程扩展到整个结构，即可得到总体的求解方程。

2.2 平面四边形八节点与四节点等参单元

由于实际结构形状复杂，实际单元可以为任意形状。通过坐标变换，可以将总体坐标系xy下任意形状的四边形单元变换成局部坐标系ξη下的正方形单元。以平面八节点单元为例，如图2所示。

图2 坐标变化示意图Fig.2 Sketch of coordinate change

采用与式（3）中相同的形函数Ni作为坐标变换插值函数，则单元几何形状和单元内场函数采用了相同数目的节点参数及相同的插值函数进行变换，即为等参变换，响应的单元被称为等参单元。因此，坐标变换也可表示为

基于以上二维平面问题的等参变换，采用局部坐标系ξη表示式（7）中的单元刚度矩阵为

式中：t为平面单元厚度；J为坐标变换雅各比矩阵，可由式（11）得到

2.3 功能梯度材料的等参梯度单元

式（12）中的D为材料刚度矩阵，各项同性材料在平面应力状态下的材料刚度矩阵为

式中E和ν分别为材料的弹性模量和泊松比，将E替换为E(1-ν2)、ν替换为ν(1-ν)即可得到平面应变状态下的材料刚度矩阵。在单元内部，功能梯度材料的材料性能指标也由形函数Ni插值得到，即

将式（15～16）代入式（14）中得到

再将式（17）代入式（12），即可得到功能梯度材料等参梯度单元的单元刚度矩阵。

2.4 数值处理方法

本文采用高斯积分方法求解式（12）中的单元刚度矩阵，则

式中：m为每个积分方向上的积分点个数；ξi、ηj为积分点坐标；Wi、Wj为积分点处相应的权系数；f(ξi，ηj)为式（19）函数在积分点坐标ξi、ηj处的取值。

对于二维平面四边形四节点和八节点单元，其内部高斯积分点的分布如图3 所示。四节点单元内部共4 个高斯积分点，即每个积分方向2 个积分点；八节点单元内部共9 个高斯积分点，即每个积分方向3 个积分点。积分点坐标和权系数如表1 所示。

图3 单元积分点Fig.3 Integration points of the element

表1 积分点坐标与权系数Table 1 Integration point coordinates and weight coeffi⁃cient

求得单元刚度矩阵ke后，可利用ABAQUS 进行二次开发功能，将等参梯度单元模型编写为单元用户子程序UEL，应用于功能梯度材料结构的数值模拟。单元用户子程序UEL 的介绍及处理流程见参考文献［18，21］。

2.5 后处理云图显示方法

由于UEL 单元只参与ABAQUS 计算，将其单元刚度矩阵用于集成总体刚度矩阵，而单元形状、应力云图无法直接在ABAQUS 的后处理中显示。因此必须采用其他手段进行后处理查看云图。常用的方法是将UEL 计算结果储存到文本文件中，然后利用Matlab 等工具读取文本数据绘制云图。在将UEL 单元与ABAQUS 自带单元的计算结果进行对比时，Matlab 等工具绘制云图时采用的插值方法必须与ABAQUS 绘制云图的插值方法一致，以避免插值方法不一致引入的误差。而ABAQUS 后处理中不同类型单元所采用的插值方法不同，因此使用Matlab 等工具绘制云图时也必须采用多种插值方法，过程比较复杂。

本文的后处理云图显示方法为：

（1）将UEL 单元计算结果储存到文本文件中。

（2）将计算模型中的UEL 单元替换为与UEL节点、积分点一致的ABAQUS 自带单元，得到新的计算模型。

（3）在新的计算模型中采用用户子程序SIGI⁃NI，读取文本文件中与该积分点坐标相匹配的UEL 单元计算结果。

（4）对新的计算模型施加“零载荷”，进行一次“假计算”，则后处理的显示结果即为SIGINI 读取的UEL 单元计算结果。

其中SIGINI 为ABAQUS 用户子程序［21］，可用于定义计算模型的初始应力场。在本文中用于将UEL 单元的计算结果赋予与其节点、积分点一致的ABAQUS 自带单元，然后进行“零载荷假计算”，从而生成的应力云图即为UEL 单元的应力云图。

3 单元验证

3.1 模型说明与解析解

采用本文等参梯度单元对图4 所示的功能梯度材料正方形平板开展有限元分析，验证单元的收敛性，并与采用常规单元分层法［13］的计算结果进行对比，验证本文等参梯度单元在计算精度和效率上的优势。

图4 计算模型示意图（单位：mm）Fig.4 Sketch of computed model(Unit:mm)

正方形平板为材料A 和材料B 混合制备的功能梯度材料，边长100 mm，厚度取单位厚度。考虑功能梯度材料的弹性模量如式（1）线性变化或式（2）指数变化的两种情况，泊松比为常量0.3。令材料A 的弹性模量为100 GPa、材料B 的弹性模量为800 GPa。此处的材料性能数据用于本节等参梯度单元的验证，不具有实际物理意义。

正方形平板受到的约束条件为在（0，0）点受x向和y向约束，在y=0 边受y向约束。开展3 种载荷条件下的计算分析。

工况1：位移载荷条件，在y=100 边施加位移Δy=1.0。

工况2：应力载荷条件，在y=100 边施加均布拉力f=100 N/mm。

工况3：弯曲载荷条件，在y=100 边施加线性分布力f（x）=100-2x（N/mm）；f（x）的正负分别表示指向y轴正向或负向。

当功能梯度材料的弹性模量线性变化时，3 种工况下y=0 坐标处y向应力的解析解分别如式（20～22）所示，推导过程见参考文献［15］。

3.2 收敛性

采用本文的等参梯度单元计算工况1 下线性功能梯度材料、指数型功能梯度材料的应力分布分别如图5、6 所示。图中S为应力，S22表示y向应力分量。图5 中平面四节点单元与平面八节点单元的计算结果采用的单元数目依次为1×1、3×3、6×6、10×10。从图5 可以看出，随着单元密度增加，应力峰值和云图分布均保持稳定不变，表示已得到收敛结果。图6 中平面四节点单元与平面八节点单元的计算结果采用的单元数目依次为1×1、3×3、10×10、30×30。从图6 可以看出，随着单元密度增加，应力峰值和云图分布均发生变化，并最后趋于稳定，得到收敛结果。

图5 线性功能梯度材料应力云图（工况1）Fig.5 Stress contour of linearly functional graded material(Case 1)

图6 指数型功能梯度材料应力云图（工况1）Fig.6 Stress contour of exponentially functional graded ma⁃terial(Case 1)

从图5、6 可以看出，对于线性功能梯度材料，无论是采用平面四节点还是八节点单元，均只需要采用1 个单元，即可得到收敛结果，而指数型功能梯度材料性能更复杂。采用平面四节点单元时，将10×10 网格进一步细化后，应力峰值变化量小于0.3%，在工程应用中可忽略，即采用平面四节点单元需要10×10 个单元得到收敛解。此外，采用平面八节点单元时，将3×3 网格进一步细化后，应力峰值变化量小于0.2%，在工程应用中可忽略，即采用平面四节点单元仅需要3×3 个单元得到收敛解。图6（e）中四节点单元10×10 网格的节点数量为11×11=121，图6（d）中八节点单元3×3 网格的节点数量为7×7=49，而图6（d）的计算精度高于图6（e），可见相比于四节点梯度单元，八节点梯度单元可以采用更少数量的单元或节点，而得到更高精度的计算结果。

将本文的等参梯度单元用于工况2 和工况3下功能梯度材料的计算，亦得到收敛解。为简洁起见，工况2 和工况3 的结果不再详细展示，仅给出单元数目为10×10 的计算结果，在3.3 节中可查看。

从以上分析可以看出，采用本文的平面四节点或八节点等参梯度单元均可得到功能梯度材料的收敛解，但八节点单元比四节点单元需要更少的单元数/节点数即可得到收敛结果。

3.3 计算精度和效率

本节给出了等参梯度单元与常规单元的计算结果，并与解析结果进行对比，验证本文方法在计算精度和效率上的优势。等参梯度单元内部的材料呈梯度特性，而常规单元内部材料为均匀材料。采用常规单元模拟功能材料，需要采用分层法［13］，对于不同层的单元赋予不同的材料属性以表征材料梯度。无论采用平面八节点还是四节点梯度单元，计算精度均优于平面八节点常规单元，更优于平面四节点常规单元。为简洁起见，本节只给出10×10 平面八节点梯度单元与常规单元的计算结果对比。本文采用的常规单元均为ABAQUS 自带的全积分单元，八节点常规单元类型为CPS8，四节点常规单元类型为CPS4。

3.3.1 工况1 计算结果对比

工况1 下线性和指数型功能梯度材料的应力分布分别如图7、8 所示，提取坐标y=0 上的应力值，与式（20）、式（23）的解析计算结果对比分别如图9、10 所示。从图中可以看出，梯度单元的计算结果与解析解吻合较好，而常规单元与解析解误差较大。对于线性功能梯度材料，10×10 常规单元（图7（b））的计算精度低于1×1 四节点或八节点梯度单元（图5（a）或（b））的计算精度。而对于指数型功能梯度材料，10×10 常规单元（图8（b））的计算精度低于1×1 八节点梯度单元（图6（b））的计算精度；也低于3×3 四节点梯度单元（图6（c））的计算精度。

图7 线性功能梯度材料等梯度单元与常规单元应力云图（工况1）Fig.7 Stress contour of graded element and conventional el⁃ement of linearly functional graded material (Case 1)

图8 指数型功能梯度材料等梯度单元与常规单元应力云图（工况1）Fig.8 Stress contour of graded element and conventional el⁃ement of exponentially functional graded material(Case 1)

图9 坐标y=0 处线性功能梯度材料应力云图（工况1）Fig.9 Stress value of linearly functional graded material at y=0 (Case 1)

图10 坐标y=0 处指数型功能梯度材料应力云图（工况1）Fig.10 Stress value of exponentially functional graded ma⁃terial at y=0 (Case 1)

3.3.2 工况2 计算结果对比

工况2 下线性和指数型功能梯度材料的应力分布分别如图11～12 所示，提取坐标y=0 上的应力值，与式（21）、式（24）的解析计算结果对比分别如图13～14 所示。从图中可以看出，梯度单元的计算结果与解析解吻合较好，而常规单元与解析解误差较大。

图11 线性功能梯度材料等梯度单元与常规单元应力云图（工况2）Fig.11 Stress contour of graded element and conventional element of linearly functional graded material(Case 2)

图12 指数型功能梯度材料等梯度单元与常规单元应力云图（工况2）Fig.12 Stress contour of graded element and conventional element of exponentially functional graded material(Case 2)

图13 坐标y=0 处线性功能梯度材料应力云图（工况2）Fig.13 Stress value of linearly functional graded material at y=0 (Case 2)

3.3.3 工况3 计算结果对比

图14 坐标y=0 处指数型功能梯度材料应力云图（工况2）Fig.14 Stress value of exponentially functional graded ma⁃terial at y=0 (Case 2)

工况3 下线性和指数型功能梯度材料的应力分布分别如图15～16 所示，提取坐标y=0 上的应力值，与式（22）、式（25）的解析计算结果对比分别如图17～18 所示。从图中可以看出，梯度单元的计算结果与解析解吻合较好，而常规单元与解析解误差较大。

图15 线性功能梯度材料等梯度单元与常规单元应力云图（工况3）Fig.15 Stress contour of graded element and conventional element of linearly functional graded material(Case 3)

图16 指数型功能梯度材料等梯度单元与常规单元应力云图（工况3）Fig.16 Stress contour of graded element and conventional element of exponentially functional graded material(Case 3)

图17 坐标y=0 处线性功能梯度材料应力云图（工况3）Fig.17 Stress vaule of linearly functional graded material at y=0 (Case 3)

图18 坐标y=0 处指数型功能梯度材料应力云图（工况3）Fig.18 Stress vaule of exponentially functional graded ma⁃terial at y=0 (Case 3)

根据工况1～3 的计算结果对比可以看出，常规单元由于采用分层法，每层单元的材料属性不同，边界处材料性能发生阶跃，计算的应力场也在边界处也产生阶跃，与实际情况不符，特别是当载荷条件较复杂时，应力场严重失真，如图11（b）与图12（b）所示。而梯度单元将材料梯度定义在了单元内部，避免了单元边界处的材料性能阶跃，得到的应力场光滑连续。梯度单元可以采用较少数目的单元，即可得到比常规单元更为精确的计算结果，无论是在计算精度和效率上均优于常规单元。

4 结构分析算例

4.1 开孔结构拉伸应力分析算例1

4.1.1 模型说明

本节算例模型如图19 所示［15］，为单位厚度的功能梯度材料带孔平板结构件，端部受均布拉伸载荷q作用。该平板的基本组分材料为TiB 和Ti，TiB 的弹性模量和泊松比分别为375 GPa 和0.14，Ti 的弹性模量和泊松比分别为107 GPa 和0.34。该结构件材料性能分布关于x轴对称、沿y轴方向呈指数型变化，如式（26～27）所示。

图19 带孔结构拉伸示意图（单位：mm）Fig.19 Sketch of perforated structure stretching (Unit:mm)

式中W为试验件宽度，其值为14 mm。

4.1.2 结果对比

采用八节点梯度单元和常规单元计算的应力结果分别如图20～21 所示，图中单元特征尺寸为该计算模型中最小单元面积的平方根，用于表征网格密度。图中S为应力，S11表示x向应力分量。单元特征尺寸越大，网格密度越小，反之网格密度越大。

图20（a～c）单元密度依次增大，应力结果收敛较快。当单元特征尺寸由 0.531 mm 变为0.266 mm 时，应力峰值变化量小于0.05%，在工程应用中可忽略，即八节点梯度单元采用单元特征尺寸为0.531 mm 的计算模型（图20（b））即可得到收敛解。

图20 8 节点梯度单元应力计算结果（算列1）Fig.20 Computed stress results of 8⁃node graded element(Example 1)

从图21（a～d）可以看出，由于采用了分层单元法，单元边界处材料性能发生阶跃，常规单元计算的应力场呈现明显的锯齿状，局部应力场失真。当采用常规单元时，需要更密的网格才能得到收敛解，计算精度和效率大大低于梯度单元。此外，当网格密度较小时，计算的最大应力位置位于左下方内圆弧上，如图21（a）所示。而随着网格密度增大，计算的最大应力位置转移到由上方外圆弧上，如图21（c，d）所示。可以看出，当网格密度不够时，计算得到的最大应力位置和数值均与实际不符。

图21 8 节点常规单元应力计算结果（算列1）Fig.21 Computed stress results of 8⁃node conventional ele⁃ment(Example 1)

4.2 悬臂梁弯曲应力分析算例2

4.2.1 模型说明

本节算例模型如图22 所示，为单位厚度的功能梯度材料悬臂梁，端部受集中力作用而发生弯曲变形。该悬臂梁的基本组分材料同算例1，料性能分布关于x轴对称、沿y轴方向呈指数型变化，如式（26～27）所示。本算例中W′为悬臂梁高度，其值为10 mm。

图22 悬臂梁弯曲示意图（单位：mm）Fig.22 Sketch of cantilever bending (Unit:mm)

4.2.2 结果对比

采用八节点梯度单元和常规单元计算的应力结果分别如图23～24 所示。图23 和图24 中，（a～c）单元密度依次增大，采用的单元数目依次为2×6、4×12、10×30。随着单元密度增大，图23 应力结果收敛较快。当单元数目由4×12 变为10×30 时，应力峰值变化量小于1.1%，在工程应用中可忽略，即八节点梯度单元采用单元数目为4×12的计算模型（图23（b））即可得到收敛解。而图24随着单元密度增大应力结果收敛较慢，当单元数目由4×12 变为10×30 时，应力峰值变化量大于12.7%，计算结果尚未收敛。由以上分析可以看出，当采用常规单元时，需要更密的网格才能得到收敛解，计算精度和效率大大低于梯度单元。

图23 8 节点梯度单元应力计算结果（算列2）Fig.23 Computed stress results of 8⁃node graded element(Example 2)

图24 8 节点常规单元应力计算结果（算列2）Fig.24 Computed stress results of 8⁃node conventional ele⁃ment(Example 2)

此外，从图24（a～c）还可以看出，由于分层法单元边界处材料性能阶跃，当网格密度较小时，常规单元计算的应力场呈现锯齿状（单元数目2×6、4×12），直到网格密度增大到10×30，应力场的锯齿状才逐步消失。而图23 采用梯度单元，当单元数目为2×6 时，所计算得到的应力场也比较光滑。将图23 和图24 中相同单元数目的计算模型对比可以看出，梯度单元的计算精度大大高于常规单元。

5 结 论

（1）将本文的平面四节点或八节点等参梯度单元用于功能梯度材料在位移载荷、应力载荷、弯曲载荷3 种工况下的计算，均可得到收敛解，但八节点单元比四节点单元需要更少的单元数/节点数即可得到收敛结果。

（2）将本文的梯度单元与常规单元对比，相同单元数目下，无论平面八节点还是四节点梯度单元，计算精度均优于平面八节点常规单元，更优于平面四节点常规单元。

（3）本文的梯度单元将材料梯度定义在单元内部，避免了单元边界处的材料性能阶跃，得到的应力场光滑连续，可以采用较少数目的单元，即可得到比常规单元更为精确的计算结果。

（4）将本文的梯度单元用于结构计算，需要较少数目的单元即可得到收敛解，计算精度和效率大大优于常规单元。