Gyroid结构力学性能及数值收敛性研究

2022-12-25蒋创宇张保强王存福罗华耿胡杰翔曹龙超

中国机械工程 2022年23期

关键词：六面体体素模量

蒋创宇张保强陈云王存福罗华耿胡杰翔曹龙超

1.厦门大学航空航天学院，厦门，361102

2.华中科技大学航空航天学院，武汉，430074

3.武汉纺织大学机械工程与自动化学院，武汉，430200

0 引言

三周期极小曲面(triply periodic minimal surface，TPMS)结构存在于自然界中，如部分昆虫、骨头、珊瑚等[1]。近年来，由于TPMS结构在力学、声学、光学等领域具有优异性能，且与传统栅格结构相比，其曲面上每点的平均曲率为0，有效地降低了局部应力集中，因此在理论方法研究及工程应用中受到越来越多关注[2]。增材制造技术克服了TPMS难以落地的问题，为轻质结构的开发提供了更多选择[3-5]。

TPMS结构设计者面临的一个关键挑战是如何选择合适的设计变量使其力学性能更优，其中TPMS结构的刚度特性是设计的主要关注点[1]。在力学特性方面，MASKERY等[6]讨论了螺旋二十四面体Gyroid结构(Gyroid cellular structure，GCS)性能和单元尺寸间的关系，指出较小的胞元尺寸可避免由于局部失效而导致的低应变结构失效。YAN等[7]对TPMS结构进行了试验和理论研究，证明TPMS的力学性能与其体积分数有很好的相关性，并推导得到了相对模量和抗压强度以及体积分数间的相关性方程。杨磊[8]设计并制备了不同梯度和体积分数的GCS，通过压缩试验与仿真对梯度GCS结构各向异性静力特性进行了系统分析。张明康[9]利用选区激光融化技术(selective laser melting，SLM)制造了正六面体试件、三明治多孔试件，通过压缩试验、弯曲试验分别对GCS结构的压缩性能与弯曲性能进行了研究。陈剑勇[10]通过重复试验，对TPMS结构应力应变曲线和结构变形模式进行分析，研究了胞元尺寸、拓扑形式和密度梯度等参数对TPMS结构准静态压缩力学性能的影响。

随着TPMS结构研究的深入，许多学者使用有限元法进行研究，但TPMS结构的复杂性一定程度上阻碍了有限元方法的应用。MAZUR等[11]对准静态压缩工况下不同拓扑形式的单个、多个胞元结构的力学性能进行了模拟，量化了胞元尺寸和约束情况对单个胞元力学性能的影响，对多胞元结构仿真与试验变形结果进行对比，验证了有限元方法的可行性。PENG等[12]针对不同宽高比的GCS进行分析，提出了一种改变几何参数控制GCS各向异性的策略，用有限元方法研究了GCS压缩状态时的各向异性弹性响应。MONTAZERIAN等[13]对TPMS结构等效弹性模量的数值收敛性进行了研究，由试验结果拟合得到单元数与临界误差的关系式，预测相对模量收敛时，每个GCS单元数约为6.4万；AFSHAR等[14]和厉雪等[15]对梯度TPMS结构进行了研究，以压缩过程中的应力应变曲线为指标，对数值收敛性进行了讨论，指出每胞元2万至3万体素单元时即收敛。AREMU等[16]使用简化模型对多种TPMS结构收敛性进行研究，通过有限元研究了多种TPMS结构的力学性能，表明晶格结构的性能在很大程度上取决于胞元的拓扑结构。MASKERY等[17-18]研究表明，每个胞元大约5万个单元时足以将仿真误差降低到可接受的水平。PENG等[12]研究发现，当单元尺寸小于0.2 mm时，GCS的等效弹性模量的值几乎恒定，结构尺寸3×3×3和4×4×4的GCS弹性模量偏差小于5%，这一结论与文献[17-18]的结论一致。

目前TPMS结构主要通过隐式函数设计得到STL文件，该格式用于增材制造，但难以直接用于有限元仿真，而由STL获取可自由编辑的实体CAD模型存在许多障碍[19]。刘伟洛[20]联合使用MATLAB、Rhino、HyperMesh等软件对Primitive结构进行处理，将STL转为CAD模型，并划分了规则六面体单元，但由于拓扑形式的不同，GCS难以实现全规则六面体单元的划分。也有学者[8-9]绕过STL转为CAD的步骤，使用商业软件Deform对TPMS结构进行仿真，但该软件对GCS结构仅支持四面体网格划分，难以降低节点、单元数。目前许多研究者[14，16-17]基于完美体素单元方法，通过在封闭区域添加节点生成单元，用于有限元仿真。虽然这种方法的前处理十分便捷，但基于体素的六面体网格使曲面阶梯式化，失去了表面的连续性和光滑性，少量体素单元计算的准确性难以保证。因此，需要对基于TPMS的体素化所引起的计算偏差进行分析。

已有学者针对正六面体GCS的数值收敛性和力学性能做了大量工作，但仍有部分问题需要解决：一是体素化方法存在少量单元难以精确描述TPMS结构，而大量单元难以对尺寸较大结构计算求解的矛盾；二是沿某方向胞元数变化的非正六面体结构的力学特性的相关规律较为少见。

本文通过拉伸试验获得材料参数，为提高计算效率使用基于扭曲单元的体素化方法对GCS试件进行数值收敛性研究；通过拉伸试验验证方法的可行性。以此为基础，选择常见的拉伸与弯曲工况，进一步讨论量化了变厚度GCS与正六面体GCS的力学性能差异，通过解析解与有限元方法的对比，为GCS结构设计给出参考。

1 设计和试验方法

1.1 GCS的设计与构造

GCS曲面的级数近似有不同的表述方式[21]，通常在笛卡儿坐标系中表达的实函数更易实现：

G(x，y，z)=CxSy+CySz+CzSx-c

(1)

其中，Sx、Sy、Sz，Cx、Cy、Cz为正弦与余弦函数，可通过下式得到：

(2)

ki=2πni

(3)

式中，c用于控制胞元体积分数，c=0.615；i=x，y，z；Li为胞元结构在x，y，z方向上的尺寸，Li=5；ki为胞元的周期，ki=2π；ni为x，y，z方向上的胞元个数。

待所需参数确定后，根据点、线、面、体的建模思想，首先编写GCS结构表达式，然后使用MATLAB生成点云、提取点云笛卡儿坐标；再使用MeshLab对点云进行初步处理生成STL面片文件，通过Geomagic Wrap对STL文件进行曲面拟合，并进行体偏差的量化(平均偏差±0.002 mm)；接下来，在SolidWorks中对曲面进行缝合，根据GCS尺寸建立大小为5 mm×5 mm×5 mm的平面区域，使用曲面剪裁工具将GCS曲面多余部分剪裁掉，剩余部分与平面区域进行缝合，形成封闭空间；最后与实体求交得到单个胞元，在此基础上进一步通过阵列命令即可得到拉伸试件。图1为GCS设计与构造示意图。

(a)MATLAB点云 (b)曲面重建与偏差对比

1.2 相对模量

ASHBY等[22]基于小应变假设的前提，对开孔泡沫材料的立方模型提出了表征多孔固体力学性能与体积分数(相对密度)的幂函数模型，该模型被广泛应用于预测多孔结构力学性能的模型中，针对GCS而言，这一理论同样可用于量化其力学性能。以相对模量E*为指标对GCS的数值收敛性及力学性能进行研究。多孔结构的等效弹性模量与实心材料弹性模量之间的关系为

(4)

式中，E*为相对模量；El为GCS的等效实心结构弹性模量；Es为材料的弹性模量。

根据图2所示的载荷与约束，GCS的等效弹性模量为

图2 载荷、约束示意图

(5)

式中，F为GCS顶面沿载荷方向的合力；L为GCS的原始长度；A为栅格区域的横截面积；u为GCS顶面在加载方向上的位移，为满足小变形和材料线性假设，令u为L的1%。

2 数值收敛规律及力学性能

要建立能精确表示TPMS几何构型的有限元模型要从计算精度和效率等方面考虑。弹性体有限元求解中计算误差主要来自两方面，一是由于单元内假设位移场与物体真实位移场不一致；二是求解基本方程时计算过程的舍入误差。可通过细化网格、增加单元与节点数来减小误差，对于试件级的GCS有限元模型而言，节点与单元数可能呈幂次增加，是需要避免的。而且对于复杂曲面，通常需要大大增加单元数，整体刚度矩阵随之扩大，反而会增大舍入误差，此外，对计算机硬件要求很高，会带来额外的计算开支。

2.1 网格划分

目前，已报道过的TPMS结构网格划分方式有5类：壳模型、实体模型、均匀化、超单元模型、体素化。这5种建模技术各有优缺点[23]。

为提高前处理效率，使用适应性广、建模过程高效的体素化方法进行GCS有限元网格的划分。文献中多以图3所示的直边体素网格建立有限元模型，图中深蓝色为体素网格，红色部分为GCS单胞的几何结构，随着单元尺寸的减小，红色部分与体素网格趋于一致，但仍存在阶梯变化。

图3 体素化单元尺寸变化图

在实际仿真中发现，对于TPMS等复杂曲边、曲面结构，仅采用直边的体素单元仍会产生较大的误差。为更好地对几何结构进行近似，需要进行有限单元划分，但曲边单元的产生无法避免[24]。曲边单元的产生需要引入坐标变换，在有限元法中最普遍采用的变换方法是等参变换，即单元的几何形状和单元内的场函数采用相同数目的节点参数和相同的插值函数进行变换[25]。借助等参元可以对任意几何形状的工程问题进行有限元离散，反映理想单元与实际单元变换关系的参数主要是雅可比值[26]，因此在体素化方法的基础上，通过设置不同的雅可比值与单元尺寸参数作为划分单元的约束，研究了不同约束参数组合下的有限元模型的精度和效率。

采用HyperMesh作为GCS结构网格的前处理工具，在网格工具Shrink Wrap中对单元最小雅可比(min Jacobian，MJ)值进行限制[27]。图4所示为单元尺寸相同、MJ值不同的六面体单元。

(a)MJ值为0.34 (b)MJ值为0.55

图5所示为不同雅可比值约束下得到的单胞网格模型。结合图3可看到，单元尺寸对体素网格近似几何的保真度也有较大影响。为研究雅可比值与单元尺寸对计算精度的影响，使用控制变量法对雅可比值与单元尺寸进行了单独分析。

(a)MJ值为1 (b)MJ值为0.8 (c)MJ值为0.5 (d)MJ值为0.3

2.2 数值收敛性分析

对体积分数为0.3、沿X、Y、Z轴排列形式为1×1×1的GCS数值收敛性进行分析，载荷及约束如图2所示。使用ANSYS静力模块进行仿真，所用单元类型为8节点六面体单元。X方向上的位移载荷(1%的应变)随时间线性变化，作用于模型顶面的节点。底部平面的节点在X方向上是固定的，而在Y和Z方向上的平移不受约束。后处理中读取约束处的支反力，代入式(4)、式(5)可得到GCS的相对模量。

2.2.1网格敏感性分析

对体积分数为0.3，排列形式为1×1×1的GCS数值收敛性进行分析，仿真结果如图 6、表 1所示，可以得出以下结论：

图6 相对模量随体素单元大小及雅可比值变化

由图6所示的相对模量变化曲线和表1中的具体数据可以发现，体素单元尺寸与MJ值对GCS仿真结果均有影响，且MJ值对结果的影响相比仅改变单元尺寸的影响更大。MJ值为0.5和0.3时，结果相差最大仅7%，继续减小MJ值对结果影响不显著，且可能因雅可比值过小导致计算精度下降。

表1 不同雅可比值与单元尺寸组合下的结果及误差

随着MJ值的减小，几何保真度增大，仿真结果与收敛值之间的相对偏差显著减小。且根据MJ值不同取值的趋势来看，取1和取0.8时单元尺寸由0.1 mm细化至0.05mm，但相对模量E*收敛参考值的相对误差eE*的改变量为19%和9%，未收敛；取0.5时，单元尺寸为0.05 mm与0.1 mm时，eE*在3%以内；取0.3时，单元尺寸从0.5 mm缩减至0.05 mm，eE*逐渐减小，单元尺寸取0.05 mm与0.1 mm时，单元数增加6倍，eE*仅为1%，可认为结果收敛。

体素单元尺寸以相对误差不超过5%，单位GCS单元数不超过10万为限制，兼顾计算精度和效率。表1第一列中，a_b表示生成单元时HyperMesh中设置的参数，a为该单元中的MJ值，b为单元尺寸。结合表1结果，可选择0.3_0.2(eE*为4.85%，单元数为11 457)、0.5_0.1(eE*为4.52%，单元数为46 228)、0.3_0.1(eE*为1.06%，单元数为63 600)。考虑到硬件原因(处理器：Intel(R)_i7-4770_CPU_@_3.40 GHz；硬盘：Western Digital，接口类型 SATA2(3Gb/s)容量 1TB，转速7200 r/min，缓存 64 MB；内存：金士顿，容量16 GB，DDR3 1333 MHz)，选择0.3_0.2这一组合进行后续的仿真计算，在精度相差不大的情况下显著减少了单元数量。

对比不同文献中数值收敛性研究结果，因难以找到完全相同的GCS，相对模量值难以度量，故选择收敛时每胞元所含有限单元个数为标准。表2中“与MJ值为0.3的单元数之比”表示达到收敛时每个TPMS胞元所含有限单元数与MJ值为0.3时所含有限单元数的比值，由此可知用扭曲体素单元可显著减少收敛时单元数。

表2 不同收敛下每胞元包含的体素单元数

2.2.2GCS正六面体排列结构收敛性分析

在讨论了网格敏感性的基础上，对体积分数为0.3、排列形式为2×2×2至6×6×6的GCS数值收敛性进行分析。如图7所示，GCS单个胞元尺寸不变，结构以正六面体形式增大，一至六阶的GCS相对模量也表现出增大的趋势，且随着阶次的增加，相对模量变化量逐渐减小。三阶相对模量与收敛值相比，相对误差为0.76%，但体素单元数相差约8倍；四阶相对模量与收敛值相比，相对误差仅为0.32%，体素单元数相差约3倍。这一结果在文献[12]、文献[17]中也得到了证实，因此使用四阶GCS仿真更优。

图7 相对模量随阶次的变化

随着GCS阶次的增加，GCS的相对模量逐渐逼近收敛值。这一现象可看作是文献[17]中给出的关于整体效应的另一种解释：随着结构中GCS数目的增大，变形中引入了一定程度的均匀性，因为GCS边缘的自由面受到相邻GCS的约束，可以给出更精确的变形描述，内部GCS可认为是均匀多孔固体的一部分。

2.3 拉伸试验与仿真

用于试验的GCS表达式如式(1)所示，取单个GCS的尺寸为5 mm×5 mm×5 mm，体积分数为0.3。基于1.1节中的GCS设计与构造方式，通过控制体积分数、单个GCS尺寸、总体结构尺寸等参数，在SolidWorks中生成拉伸试件的几何文件，几何尺寸如图8所示，图8a所示为用于打印加工的STL模型，图8b为加工实物图。

(a)STL模型 (b)实心试件与GCS试件

2.3.1GCS制造

试验试件来自西安铂力特增材技术股份有限公司，采用选区激光融化(SLM)技术，所用设备为BLT-S310打印机。所用材料为Ti6Al4V，粉末颗粒尺寸和形态见表3和图9，粉末颗粒球形度达到0.92，球形度高，颗粒尺寸在15～53 μm，Dv(50)值为38.110 μm，位于颗粒尺寸范围中带，满足SLM工艺粉末颗粒尺寸和形态要求。

表3 粉末元素质量分数

图9 颗粒形态

2.3.2材料参数及模型

通过图10所示的实心试件的单轴拉伸试验来确定该材料的弹性模量和单轴抗拉强度。以3 mm/min的加载速率，在电子万能试验机HF-JL-005上进行单轴拉伸试验，执行标准为GB/T 228.1—2021和GB/T 7314—2017，拉伸试验过程如图10所示。

图10 拉伸试件载荷、约束示意图

对同样热处理方案下的实心样件试验数据进行处理，得到了Ti6Al4V材料的塑性应力-应变曲线[28]。其中，线弹性段拟合得到弹性模量Es为105 860 MPa，屈服强度σs为830 MPa，塑性应力应变数值如表4所示。

表4 实心试件拉伸塑性数据

2.3.3仿真模型设置

使用显式动力学模块Ls-Dyna模拟GCS拉伸过程，材料模型为分段线性(MAT_24)，所用单元类型为8节点六面体，常应力单元(SECTION_SOLID_ELFORM=1)。参照图10所示的载荷及约束，试件下夹头端节点施加固定约束(BOUNDARY_SPC_SET)，上夹头施加位移(BOUNDARY_PRESCRIBED_MOTION_SET)为1.6 mm，该数值取自拉伸试验试件断裂时标距段位移改变量，为提高计算效率，仅取GCS试件标距段部分进行仿真。使用材料卡片中MAT_ADD_EROSION中的最大塑性应变判断准则作为判断单元的失效准则，失效应变为0.08。

2.3.4GCS拉伸断裂分析

图11a、图11b为实心样件分别在12倍和1000倍的放大倍数下的扫描电镜(SEM)图，12倍下失效截面中部呈纤维区，外部为剪切唇区。用1000倍放大观察纤维区内截面失效机理，可以看到典型韧窝形貌。图11c、图11d给出了12倍GCS形貌及250倍微观失效界面形貌。低倍下可以看出样件表面粗糙，有很多颗粒夹杂物。在SEM观察下，截面杯口形状和纤维区等更直观。微观结果进一步印证了样件符合典型塑性材料断裂失效机制。另外，高低倍数下都可观测到未融化的钛合金粉末(图中虚线圆处)，通常是失效先开始的地方，由此可能带来材料性能的劣化。

(a)实心试件12放大倍 (b)实心试件放大1000倍

2.3.5GCS拉伸断裂分析

因试验条件有限，GCS表面无法粘贴应变片，未测得应变数据，故此处仅对力随时间变化历程以及抗拉强度、极限载荷等进行对比。

如图12所示，仿真与试验所得力随时间变化逐渐增大，趋势一致，然后进一步增长即强化阶段，到达抗拉强度后试件断裂，力迅速减小，即卸载与破坏，这与典型金属材料拉伸过程相符。试验曲线在5 kN左右进入小的平台期，发生屈服，仿真中未观察到类似过程。表5所示为仿真得到的极限载荷和抗拉强度与试验结果对比，抗拉强度、极限载荷误差为1.45%。

表5 试验与仿真结果对比

图12 仿真与试验力随时间变化曲线

图13所示为仿真与试验过程中GCS的损伤情况，在拉伸前期，仿真与试验均可观察到局部断裂，位置基本一致；继续拉伸，发生整体断裂，整体断口路径位置与方向也趋于一致。

(a)加载过程中局部断裂

3 变厚度GCS力学性能

已有研究表明，排列方式为5×5×5以上的GCS，其等效弹性模量已达稳定值[18]，且SIMSEK等[23，29]以有限元得到的等效弹性模量为有限元仿真的输入，对夹层结构进行等效并分析了GCS板的模态特性。但对非正六面体GCS是否也适用正六面体GCS的等效弹性模量，还未见报道。且结构设计中，薄板、梁等非正六面体结构的应用范围更为广泛，因此有必要对此进行分析。

以2.2节中正六面体收敛的等效弹性模量得到的计算结果为对照组，对拉伸与弯曲情况下的响应进行研究。

3.1 变厚度GCS拉伸仿真

选择与2.2节中一致的约束，载荷类型由位移变为力，对沿X、Y、Z轴排列为2×2×1、2×2×2、3×3×1、3×3×2、3×3×3、4×4×1、4×4×2、4×4×3、4×4×4、5×5×1、5×5×2、5×5×3、5×5×4、5×5×5的GCS结构进行分析，体素单元尺寸为0.2 mm，最小雅可比值为0.3。

使用1.2节中介绍的相对模量计算方法，对上述排列形式的GCS模型相对模量进行计算。相对模量随Z向胞元数目变化结果如图14所示，可以得出以下结果：

图14 相对模量随Z方向层数的变化

(1)GCS沿Z轴方向的胞元数目变化对其力学性能影响是显著的。随着二阶到五阶非正六面体厚度方向的增大，GCS相对模量变化量逐渐减小，其中最大相对误差出现在5×5×1和5×5×5组，可达14.15%，其次是4×4×1和4×4×4组，可达14.10%。

(2)2×2×2～5×5×5系列GCS随着Z轴方向的胞元数目的增加，达到正六面体时相对模量E*达到最大，且Z方向胞元数为1和2时，E*改变量总是最大，约为10%。随着厚度方向层数增加，E*改变量逐渐减小。

为更好地对结构设计进行指导，进一步对沿Z轴方向的胞元数目变化的力学性能进行分析，结合1.2节中式(4)，得到等效实体的弹性模量，结合图14分析结果，对GCS排列方式为4×4×1、4×4×2、4×4×3、4×4×4的结构，采用解析法和有限元方法对单轴拉伸工况时微小变形(线性弹性)阶段的位移结果进行对比。

因结构位移改变量远小于结构沿载荷方向的尺寸，材料满足理想线弹性假设，故根据胡克定律，位移的解析解为

(6)

式中，d为加载后载荷方向的位移改变量。

计算所需参数见表6。表6中最后一行给出正六面体GCS收敛的数值，最后一列为变厚度结构的等效弹性模量与正六面体等效弹性模量的相对误差Δe，该误差定义为变厚度GCS等效弹性模量与收敛值之差再与收敛值之比，其中4×4×1与6×6×6达到14.41%。

表6 拉伸载荷解析解计算所用参数

图15所示为对照组与GCS仿真结果。图15a为采用Ti6Al4V材料、排列方式为4×4×1的GCS计算得到的位移结果云图，图15b所示为对照组计算得到的位移结果。由于图15a中GCS边角部分结构为悬空状态，局部刚度较小，在载荷下的变形远大于其他区域，因此对载荷施加面中心区域均匀选取4个点的位移结果取平均值作为GCS最终的整体位移。

(a)GCS仿真结果

表7列出了解析解与GCS、对照组的计算结果。GCS结果与解析解相对误差小于3%，表明在弹性段的计算结果是可信的。对照组模型的计算结果与解析解误差最大值为14.41%，最小为0.31%，表明排列方式为4×4×1的等效弹性模量与对照组差距最大，随着厚度增大，与解析解误差逐渐减小。

表7 拉伸载荷解析解与仿真结果对比

值得注意的是，由4×4×3、4×4×4 GCS与对照组的误差对比可以发现，对接近正立方体的结构使用对照组参数直接近似误差不超过5%，可满足工程上对精度的要求，但对于宽高比过大的结构(如4×4×1、4×4×2)则会对力学性能产生过高的估计。

对比4×4×1到4×4×4，解析解与仿真结果均表明：①在同一载荷下，层数增加对力学性能的改善会逐渐减小，即通过改变GCS排列方式对结构性能进行设计存在边界效应；②当GCS中某方向胞元数小于其他方向胞元数时，会对整体的力学性能带来负面影响，直接使用正六面体GCS的等效弹性模量进行等效可能带来较差的结果。

3.2 变厚度GCS悬臂梁弯曲仿真

结合图14的分析结果，仍分析四层GCS厚度变化的力学性能。对GCS沿X、Y、Z轴胞元排列为20×4×1、20×4×2、20×4×3、20×4×4的结构进行分析，体素单元尺寸为0.2 mm，最小雅可比值为0.3。载荷与边界条件如图16所示，固定端与刚性墙固连，自由端施加沿厚度方向变化、大小为50 N的力。