方向性热膨胀系数可设计点阵结构均匀化性能分析

王 帅，牛 飞，牛 斌

(1.大连理工大学机械工程学院，辽宁 大连 116024;2.中国运载火箭技术研究院，北京 100076)

0 引言

随着点阵结构的研究和发展，点阵结构的单胞数量越来越多，构型越来越复杂，对此类结构的分析往往很耗时甚至无法计算分析。通常的做法是采用近似模型，将宏观结构中具有非均质特性的点阵结构单胞进行等效处理，以此来简化宏观模型的分析。常见的等效性能预测方法有：以变分原理为基础的定界法，如Hashin[1]通过解析法得到材料等效参数的上限和下限范围。基于夹杂理论的近似模型：自洽模型、广义自洽模型、Mori-Tanaka(M-T)模型等，以及基于此构建的数值方法，如自洽法[2]、广义自洽法[3]、M-T法[4]。考虑到实际计算模型的复杂性和预测方法的通用性，数值求解方法在预测结构等效性能方面比较流行，主要有代表体元法(RVE)和渐近均匀化方法(AH)。

代表体元法(RVE)[5-6]基于能量等效原理，选取周期性点阵结构(代表体元)作为研究对象，通过施加单位位移边界条件(Dirichlet边界)或者力边界条件(Neumann边界)，使得点阵结构应变能和相应的均质材料单胞应变能等效，从而求出等效性能。渐近均匀化(AH)[7～9]方法是另外一种预测周期性点阵结构等效性质的数值方法。渐近均匀化方法基于摄动理论和虚位移原理，具有严谨的数学理论推导，当单胞尺寸相对于宏观结构尺寸非常小时，能够给出精确的预测结果。近年来，许多学者将渐近均匀化与商业软件结合，保留严格的数学理论逻辑，简化程序开发流程，实现快速简单计算，如程耿东等[10]提出一种新渐近均匀化实现方法，可以使用商业软件轻松实现多种点阵结构等效弹性模量的计算，张永存等[11]提出了基于渐进均匀化方法预测周期性复合材料等效热膨胀系数的新算法。

本文采用渐近均匀化方法获得点阵单胞的等效弹性模量和热膨胀系数，比较点阵结构等效分析与离散分析之间的差异。首先，本文针对一种单方向热膨胀系数可设计的点阵单胞，用渐近均匀化方法，获得点阵单胞的等效模量和热膨胀系数属性。然后，分别建立一个点阵结构离散模型和采用均匀化后等效模型，讨论在不同几何外形、边界条件、载荷情况和点阵数目下均匀化方法预测的准确性。

1 渐近均匀化理论

对于周期性点阵结构，渐近均匀化方法可以预测其等效性能。点阵结构的等效弹性模量可以用应变能的形式表示为[10]。

(1)

其中，Y是点阵单胞的定义域，|Y|是点阵单胞的体积，D为弹性矩阵，ε0(kl)为单位应变，ε*(kl)为特征应变。

式(1)可以变形为

(2)

其中，x0(ij)是初始位移场，x*(ij)是特征位移场，f(kl)是初始位移场对应的节点反力，f*(kl)是特征位移场对应的节点反力。以二维为例，特征位移场x0(ij)可以表示为[10]:

(3)

f(kl)和f*(kl)可以分别表示为

f(kl)=Kx0(kl)

(4)

f*(kl)=Kx*(kl)

(5)

所以，只要获得f(kl)和f*(kl)就可以计算有效弹性模量[10]。建立一个单胞有限元模型，对模型施加一个初始位移场x0(kl)进行一次有限元分析，提取节点的反力，即f(kl)。再将获得的f(kl)施加在模型上，施加周期性边界条件进行一次有限元分析，可以获得特征位移场x*(kl)。最后，在模型上施加特征位移场x*(kl)进行一次有限元分析，获得节点的反力，即f*(kl)。将获得的数据代入到公式(2)中即可计算出有效弹性模量。

等效热膨胀系数在等效弹性模量的基础上进行计算。等效热膨胀系数与等效弹性模量EH和热弹性常数β有关[11]，关系如公式(6)所示，所以确定热膨胀系数还需计算热弹性常数β。

α=(EH)-1β

(6)

热弹性常数β的表达式如式(7)所示：

(7)

其中，Rα为单位负温升下的节点反力，Λ为特征位移场，RΛ为在特征位移场Λ下的节点反力。

计算等效热膨胀系数的步骤如下[11]：首先，约束所有节点位移，施加单位负温升，获得节点反力Rα。然后将节点反力Rα施加在单胞的每一个节点上，施加周期性边界条件，通过一次静力分析，获得特征位移场Λ。再将各节点施加特征位移场Λ，获得节点反力RΛ。最终，将获得的节点反力Rα和RΛ代入到公式(7)中计算出热弹性常数β，再由公式(6)获得单胞的等效热膨胀系数。

2 点阵模型等效属性计算

本文使用的点阵结构单胞如图1(a)所示，在双材料金字塔结构[12]的基础上，将底面改为板结构，提高承载性能，由两个双材料板杆金字塔构型对接形成本文点阵结构单胞。其中图1(a)、图1(b)中板材为材料1，杆件为材料2。

图1 点阵结构单胞图

参考文献[13]中的推导过程，金字塔结构各边与高度h、角度β之间满足关系

(8)

分别对L1、L2求微分可得：

(9)

其中，dL1、dL2、dh、dβ分别是L1、L2、h、β的微分形式。

根据热膨胀系数的定义，图1(c)中浅色杆件(材料1)的热膨胀系数α1=dL1/L1dt，深色杆件(材料1)的热膨胀系数α2=dL2/L2dt，公式(9)变形为：

(10)

其中，dh/hdt为四棱锥沿高度h方向的热膨胀系数αh(后面简写为等效热膨胀系数)，对公式(10)约分，再将α1、α2、αh代入到公式(10)中，得到公式(11)。

(11)

(12)

从公式(12)中可以看出，选定两种不同的材料或设计角度β的大小，都可以改变等效热膨胀系数αh的大小，实现正、零和负热膨胀系数的变化。

当固定角度β，选择不同的材料组合来改变αh的数值。点阵结构单胞中材料1选择为1Cr8Ni9，材料2选择为4J33，材料属性如表1中所示。将表1中的热膨胀系数代入到公式(12)，当角β=45°，点阵结构单胞等效热膨胀系数为0。然后固定角度β为45°不变，选择不同的材料组合，获得不同等效热膨胀系数的点阵结构单胞。采用渐近均匀化方法，计算点阵结构单胞的等效材料属性，结果如表2所示。表2中Ex、Ey、Ex表示点阵结构单胞XYZ三个方向的弹性模量，αx、αy、αz表示点阵结构单胞XYZ三个方向的热膨胀系数。

表1 材料属性[11]

表2 单胞等效材料属性

由表2中第一行数据，点阵结构单胞在XY方向上的热膨胀系数与材料1的热膨胀系数相近，在Z方向上的热膨胀系数比XY方向上的热膨胀系数小两个数量级，符合低热膨胀的设计。对比3个点阵结构单胞的等效弹性模量可以看出，单胞在XY方向上的弹性模量相同，Z方向上的弹性模量较低。对比3个点阵结构单胞的等效热膨胀系数可以看出，点阵结构单胞沿XY两方向的热膨胀系数与板材料(材料1)的热膨胀系数相同，Z方向上的热膨胀系数是可以设计的，可以实现由正到低再到负热膨胀系数的变化。

3 点阵结构单胞均匀化性能验证

点阵结构单胞均匀化之后可以看作是一个各向异性的等效材料。先采用简单几何的点阵结构为例，分别建立均匀化后的等效模型和点阵离散梁/板模型，考虑不同边界条件、载荷情况和单胞数目，验证渐近均匀化方法预测的等效弹性模量和热膨胀系数的准确性。再使用复杂拓扑的点阵结构，同样建立均匀化后的实体模型和点阵离散模型，探究几何外形对渐近均匀化方法预测结果准确性的影响。

3.1 不同单胞数目

采用图1(a)中的点阵结构单胞，填充一个10 mm×10 mm×14.14 mm的长方体区域，采用不同大小的点阵结构单胞进行填充，分别形成2×2×2、5×5×5和10×10×10点阵结构，建立点阵离散梁/板模型，如图2所示。随着单胞数的增多，杆径和板厚相应变小，保证单胞的体分比不变。基于点阵结构均匀化等效性能，建立均匀化后的等效模型。载荷及边界条件为：底面固定，顶面受到一个4 MPa的均布载荷，位移计算数据如表3所示。表3中数据为顶面上的所有节点位移的平均值。

图2 点阵填充结构

表3 平压下节点位移平均值

由表3中数据可以看出，随着点阵结构单胞数目的增多，均匀化获得的等效性能越来越能反映实际情况。当单胞数量为10×10×10时，离散分析和均匀化等效分析在主方向上的误差在1%左右。之后的仿真分析就以10×10×10的点阵结构离散模型和均匀化后等效模型展开。

3.2 不同点阵边界条件

均匀化后的等效结构不论如何截取都是同样材料属性的实体，但点阵结构却不同。采用图3(a)所示的3种方式截取点阵结构作为3种不同的单胞，再以此3种单胞阵列成周期性结构，如图3(b)、(c)、(d)所示，这样形成的结构在中心大部分区域并无差异，但在边界上却不同。通过仿真分析，比较不同截取方式对点阵性能的影响。位移计算数据在表4中列出。表4中数据为顶面上的所有节点位移的平均值。

图3 点阵截取方法及阵列结构

表4 平压下节点位移平均值

横向对比表4中数据，发现基本满足上一节的结论，即点阵单胞数量越多，均匀化结果预测越准确。再纵向对比表4中数据发现，截取方式1和截取方式3，对点阵性能的影响较小，但截取方式2的离散计算与均匀化预测结果误差大。原因是，截取方式2会产生不完整点阵单胞，极大的削弱点阵结构的刚度。

3.3 不同载荷情况

改变模型的载荷及边界条件，探究不同载荷情况下均匀化方法预测的准确性。将3.1节中计算的载荷情况视为第一种载荷条件。第二种载荷条件：底面固定，顶面上一条边线上受到一个均布载荷。因为边界效应的影响，比较位移的时候不考虑施加载荷处的位移情况，比较下一层点阵结构和等效模型对应位置的位移。因为载荷的不对称性，同一个面上位移响应也是不同的，所以不再对比面上的平均位移，而是比较线上平均位移，具体位置如图4(a)所示。线上平均位移计算结果如表5所示。

图4 载荷及边界条件

对比表5中的数据，均匀化后实体模型与点阵离散模型的位移基本吻合，误差在7%左右，相比于平压情况误差有所增加。在一些位移较小的地方，如表5中的X方向位移，相对误差较大，但绝对误差数量级一致。与表3中数据对比发现，载荷条件越复杂，渐近均匀化预测的误差越大。

表5 侧压载荷下节点位移平均值

第三种载荷条件：施加一个400 ℃均匀温度载荷。为了模拟一个自由膨胀的过程，边界条件设置为约束底面和底面上两边的形式。如图4(b)所示，约束底面Z方向的自由度，约束底面上一条X方向边线Y方向自由度，约束底面上一条Y方向边线X方向自由度。均匀温度载荷下节点位移计算结果如表6。表6中数据为点阵结构顶面上所有节点位移平均值。

表6 均匀温度载荷下节点位移平均值

从表6中数据发现，两者的相对误差很大。但是，因为使用的点阵结构单胞在Z方向上的理论热膨胀系数为0，Z方向位移很小。将点阵结构单胞的两种初始材料修改为AL-5A02和1Cr8Ni9，材料属性如表1所示，获得一个新的等效材料属性，如表7，并进行如图4(b)相同的仿真分析，均匀温度载荷下节点位移计算结果如表8。

表7 单胞等效材料属性

表8 均匀温度载荷下节点位移平均值

由表8中数据可以看出，均匀化后的等效模型与点阵离散模型位移吻合良好。比较表6和表8中的数据可以发现，两次的绝对误差均在10-3左右。可以看出，表6中的结果因为本身热膨胀系数较小，导致的结果相对误差较大。

3.4 不同几何外形

本文使用的点阵结构单胞沿高度方向的热膨胀系数具有可设计性，采用表2中3种单胞，填充一简单拓扑结构，如图5所示。分别创建均匀化后的等效模型和点阵结构离散模型，约束底面Z方向自由度、后面Y方向自由度和左端面X方向自由度，施加温度载荷，比较两者的位移响应。位移计算结果如表9所示。表9中数据为点阵结构顶面上所有节点位移平均值。

图5 点阵填充拓扑结构

表9中第一行数据中，实体模型与离散模型Z向位移相对误差很大，与表6的情况相同，都是因为点阵结构Z方向的热膨胀系数为0。通过对比表9中数据，均匀化后的实体模型和点阵结构离散模型位移吻合良好。通过比较Z向位移，发现不同点阵结构可以实现结构单方向正、负热膨胀系数的设计。

表9 均匀温度载荷下节点位移平均值

4 结论

以双材料四棱锥结构为基础，针对一种单方向热膨胀系数可设计的点阵结构单胞，采用渐近均匀化方法计算点阵单胞的等效弹性模量和热膨胀系数。基于此，针对不同的点阵填充结构分别创建均匀化后等效模型和离散模型，验证在不同几何、不同边界和载荷、不同单胞数目下均匀化方法预测的等效性能的准确性。获得如下结论：

1)通过设计不同材料组合和几何构型，可以实现点阵结构热膨胀系数由正到负的变化，并采用渐近均匀化方法验证等效热膨胀系数的变化；

2)基于渐近均匀化方法计算的性能与点阵离散结构的吻合情况受点阵单胞数量的影响，当点阵单胞个数足够多时，在力载荷和温度载荷的作用下基于均匀化和离散计算的位移吻合结果令人满意。随着边界和载荷条件的复杂，相比于点阵离散分析结果，热弹性渐近均匀化的计算结果误差变大。

3)采用不同截取方式形成的单胞阵列而成的结构，中心大部分并无不同，但结构的边界却有所差异，在边界上存在不完整点阵单胞，降低点阵结构的等效刚度，使均匀化实体模型与离散模型在边界上的位移计算结果相差很大。