拓扑优化中采用增材制造填充构件的结构屈曲荷载提升设计

来源：Engineering

增材制造技术(3D打印)可实现复杂几何结构的制备，填补了传统工艺的制备盲区。拓扑优化方法，众所周知已成功用于航空航天和汽车行业的轻量化构件的创新设计，提供了一个可充分利用增材制造带来的设计自由度进行灵巧设计的平台，实现设计与制造两项技术的完美契合。然而，迄今为止，拓扑优化方法在适应增材制造技术的约束和带来机遇方面的进展非常有限。挤压式增材制造技术，可制备不同于传统实体结构的内含多孔结构的固体壳结构(图1)，这也是其独特的制备特点之一。

近期，本文作者提出了一种名为涂层方法的拓扑优化方法。相比于仅可设计实心固体结构  [图1(a)]的传统拓扑优化方法，该方法可设计内部填充多孔结构固体壳结构[图1(b)]。研究表明，涂层方法虽无法提高结构刚度，但可显著改善屈曲载荷，提高结构稳定性。此外，该研究也表明拓扑优化技术与增材制造工艺相结合具有巨大潜力。 利用传统最小柔顺性方法对部件进行拓扑优化时，一般不考虑屈曲约束。因此，拓扑优化所获得的最优设计往往是拉压构件配置结构，回避抗弯构件。

然而，构件的屈曲载荷与其弯曲刚度密切相关(如杆的欧拉屈曲与其弯曲刚度成正比)，这些结构最终破坏很可能是由于屈曲载荷，而不是材料的断裂强度。当结构体分比较少时，该问题尤为显著。 目前，一些工作试图在最小柔顺性拓扑优化中考虑屈曲约束。所建议的方法要么未能获得可信的设计结果，要么陷入最小特征值(屈曲荷载)附近大量模态聚集的困扰 。模态聚集意味着需要计算大量的特征值，导致不可承受的计算负担。强制性的屈曲分析是作为优化后的后续过程，而不是集成在优化过程中的约束。如果发现构件的稳定性不足，则执行后处理过程进一步设计以提升临界屈曲载荷。在这个过程导致设计为较优结构而非最优结构。 自然界中存在着大量具有较高屈曲荷载与重量比的结构构型。最典型的例子是动物骨骼和植物的茎，它们是由柔软多孔的内部和坚固的外壳组成。夹层结构也应用了相同的概念，其共同特征在于具有相对于质量的较高的弯曲比刚度，因此具有较高的屈曲荷载。涂层方法利用相同的机理设计较高屈曲载荷的结构构型。正如本文所证实的，涂层方法提供了一种充分利用增材制造填充结构，以提高结构屈曲载荷的实现方式。本项研究包括数值和实验两部分。

2. 方法
为证明填充结构相对于相同质量的实体结构具有更好的屈曲性能，本文比较了以下两种基于密度的拓扑优化方法：
①基于投影的传统拓扑优化方法，其将产生一个几乎完全黑白的结构；
②涂层方法，其将产生内含多孔填充结构的固体壳结构。采用拓扑优化中的标准案例MBB梁(等厚简支梁上端中部受集中载荷作用，详细参考2.3节)作为研究对象，比较了两种拓扑优化方法所得结构的柔顺性和屈曲荷载，并将屈曲分析的数值和实验结果进行了对比。为方便阐述，本文针对二维结构开展研究，但在实验研究中考虑了三维效应，因而该研究可有效地扩展到三维结构。
2.1.优化问题
本文优化问题为标准的一定体积约束下的最小柔顺性问题，其有限元列式为

                 

(1) 其中，μ是设计变量向量；c是柔顺度；K是整体刚度阵 (采用传统基于密度法拓扑优化中的定义方法：叠加插值后的单元刚度矩阵的累加)；U和F分别是整体位移向量和载荷向量；g是体积约束；V(μ)是材料体积；V*是允许使用的最大材料体积。本文基于解析的敏度和移动渐近线法(MMA)[8]更新设计变量。敏度的具体表达式在此不在赘述，详细内容可参考文献。

2.2. 涂层设计方法
传统拓扑优化方法和涂层方法均可采用过滤方法实现宏观结构特征尺寸的控制。该过程包括基于偏微分方程(PDF)的密度过滤实现光滑化，以及基于中间密度场映射[10–12]方法实现光滑问题离散化。其中，光滑的程度取决于过滤半径R(定义见参考文献[11])，而映射过程由阈值η和陡峭度β决定。 除了控制宏观特征尺寸，涂层方法可设计内含多孔填充物的固体外壳结构，表面厚度tref决定了结构的固体外壳表面。本文将填充结构等效为均质结构，其材料属性为周期性填充结构的等效宏观属性，为数值模型中计算精细的微观结构提供基础。均质化后的填充结构属性可通过密度和刚度两个均匀化参数来描述，这两个参数可通过固体材料性质λm和λE的比值表示。为满足材料的物理意义，两个参数的关系必须满足Hashin-Shtrikman (HS)边界。本文选取三角形填充结构，并假定其性能满足HS上界[14]。对于二维问题，填充结构的密度和刚度的关系已在文献[15]中给出，如图1(c)所示：

(2) 需要注意的是，该关系基于固体材料的泊松比为1/3的假设。然而，对于体积分数较低，以拉压性能为主导的三角形蜂窝填充结构，材料的泊松比的影响可忽略不计。

２．３． 实验设计
实验选择标准案例MBB梁：一个长宽比为6:1的等厚简支梁，上端中间受集中荷载作用。其数值优化结构如图2(a)和2(b)所示。设计域尺寸为300 mm×50 mm，厚度为15 mm，体分比为25 %，采用大小为0.5 mm的双线性四边形单元离散。表面厚度tref取1 mm，光滑设计域时过滤半径R=10 mm。优化过程中利用结构对称性选用模型的一半，但是对于屈曲分析，为了找到非对称模式，需要使用完整模型。本文选取相对刚度为20%的多孔材料即λE = 0.2，可根据公式(2)求解得到λm = 0.43，过滤和惩罚参数的取值参考文献[16]。

2.4.部件制造
实验试件[如图2(c)和(d)所示]采取美国Maskers Tool Works公司的Mendel Max2打印机制备，该打印机配有E3D v6的挤出热喷头，采用熔丝制造(FFF)技术。打印材料为丁苯乙烯丁烯橡胶(SEBS)†。该材料是一种具有500%的断裂伸长率的热塑性弹性体(IPE)材料，可提供充足的延展性以便在试件不断裂的情况下观察屈曲模态。由于SEBS是黏弹性材料，供应商未提供其杨氏模量，其邵氏硬度A和邵氏硬度D的值分别表示为92和40。熔丝挤出层高为0.2 mm，宽度大约为0.5 mm。填充结构的单胞尺寸是由挤出宽度(等于胞壁厚度)和宏观密度 λm(见2.7节)决定。作为一种熔丝技术，FFF工艺会产生具有一定各向异性的结构部件，在熔融沉积成型(FDM)过程中可观察到相同的现象。然而，为简单起见，假定材料是各向同性线弹性材料，该假设合理性将在第4节中讨论。材料的有效参数通过实验确定。 由于实际打印过程中填充结构的密度比预设值43%高，制造出的多孔试件比预先设置的体积稍重一点。试件的蒙皮厚度可用数字测径器测量，附加体积主要由于切片软件设置与挤出熔丝宽度之间的非完美匹配引起的。因而，在数值模型计算中，通过给填充结构分配额外的重量考虑其体积误差。优化过程中，选用20 %的填充结构刚度，实际制备的填充结构密度是固体的52 %，即填充结构刚度为27 %。在有限元模型中，采用27 %的刚度值，并且与实验进行了比较，发现该有限的刚度增加对性能比较影响较小。

2.5.实验设置
实验测试由两部分组成：
①确定有效的材料参数；
②确定屈曲荷载。该两部分实验时选用不同的三点弯曲实验测试装置，如图3(a)和(c)
所示。为了确定材料参数，在梁的荷载处放置钢支架，如图3(a) 所示。荷载的增加通过在用绳子与钢支架相连的篮子中增加100 g重量来实现(其中篮子在图中不可见)，载荷位移可直接利用运托架上精度为0.01 mm的千分表测量。由于在该试验中，梁的荷载始终处于线性范围以内，因此零点可任意设置。 实验件的屈曲荷载通过具有10 kN负荷传感器的机电试验机(Instron 6022，改造于5500R)确定，如图3(c)所示。屈曲分析中，控制位移加载速度为1.5 mm·min–1，并采用其他速度用于验证(见章节4)。数据采样的速度为每秒20次。 之所以两个实验采用不同的装置，是因为只有当力高于20 N(误差低于2 %)时，10 kN负荷传感器是精确的。但是确定材料常数的力仅仅需要5 N。

2.6.屈曲数值分析

结构发生屈曲时的荷载称为临界载荷，用Pc表示。线性屈曲问题可转换为特征值问题求解，其有限元分析(FEA)列式为：

几何刚度矩阵Kσ(uref)与参考载荷作用下的线性静力位移(应力)相关，可通过求解Kuref = Fref获得。临界载荷Pc等于参考荷载与最小特征值λc的乘积。相关屈曲模态通过Φ获得。需要注意的是对于线性屈曲问题，只能获得相关模态，无法得到实际振幅(和实际位移)。 初步试验表明，多孔结构将会产生面外屈曲模态。为与实际模型进行比较，在Abaqus中建立三维模型替代优化中使用的二维数值模型，并采用八节点六面体单元(Abaqus C3D8单元)离散。二维模型离散所使用的方形单元在三维模型(x和y方向)中继续使用。样本沿z方向拉伸15 mm形成三维模型，为了减少自由度数目，模型在z方向上单元长度是平面单元长度的2倍。

最后离散结构的单元尺寸为0.5 mm×0.5 mm×1.0 mm。数值模型使用具有均质填充性质的固体材料模拟填充结构，而不是实际打印出的三角形结构。这部分假设将在第4节讨论。 载荷均布加载在一个小面域内，其内部含有80个节点。除了二维问题中定义的边界条件，约束中心载荷节点面外运动(u3 = 0)。该边界条件还对应于力传动装置与试验试件之间的摩擦作用，可防止面外屈曲模态。 首先，求解方程(3)定义的特征值问题(Abaqus：“线性扰动，屈曲”)。使用弧长法追踪结构的几何非线性响 应，对特征值分析细化(Abaqus：“静态，弧长”)。对于每个工况，引进非对称载荷扰动模拟缺陷，其大小是特征值分析中预测得到的屈曲荷载的1%。在这里提到屈曲载荷值都是非线性分析的结果，所有工况非线性分析都比线性分析低2%以下。


2.7.填充结构的屈曲性质
尽管我们的假设是当填充结构的百分比减少时，宏观的屈曲载荷将显著提升，但是其结构将易发生局部屈曲。三角形填充结构的“屈曲强度”可由一般宏观面内应力状态预测。 图4(a)为三角形填充结构，图4(b)为单胞放大图，其中三角形边长为L，单胞壁厚为t，并且三角形的方向与所使用的物理模型一致。三角形填充结构的屈曲强度取决于比值(t/L)3[19]。当已知正三角形的形状和常数参数时，L、t和λm中只有两个参数是独立的。忽视三角形在交界处的结构重叠效应，第三个参数与另外两个参数 线性相关。因此，填充结构的屈曲强度取决于λ3m。 尽管假设填充结构的弹性性质为线性各向同性，然而屈曲强度与主应力方向相关。对于沿着三角形蜂窝墙壁方向的单轴载荷，屈曲强度可通过简化的表达式计算。

假设荷载与x轴平行，则屈曲强度为

下面将讨论如何建立本文设计所需要的填充结构的屈曲强度，其中填充密度范围为
λm  (0,1]。需要注意的是，一些假设被引入本文所研究的设计问题中。 完全实体结构可近似为一个三角形布局的框架模型，甚至是桁架模型，这是由于拓扑优化方法的本质造成的，此时材料分布使任意内部结构的弯曲最小。

主应力(绝对)值分布[如图4(b)所示]也证明了这一结论，除了个别梁的横截面由于离散产生少许不规则现象，其他大多数梁都是满足的。  假定填充密度对于拓扑优化结构没有影响，与完全实体模型获得的拓扑优化结构相同；假设任意填充密度的拓扑优化结构可采用完全相同的线性桁架结构建模，并且桁架的正应力在给定的外载荷下不发生变化。但是，从完全实体结构到具有较小均质刚度的内部多孔结构，桁架的截面积发生变化。这意味着对于任意一个杆件，填充结构的宏观轴向应力会比完全实体结构的应力低，可通过一个与表面厚度和填充结构密度有关的因子表示，尽管实体壳中的应力会高一些。通过对完全固体结构所对应的桁架模型面内应力乘以“填充应力参数” [参见图4(c)]，可获得任意填充密度的多孔结构杆对应的面内应力状态。

图4(d)中箭头指向的杆(或者镜像于垂直对称线的杆)在所有压缩杆件中主应力绝对值最大。这个杆近似满足x轴方向单轴压缩，因此，其屈曲强度可采用表达式(4)表示。这个表达式与杆的正压力之间的比值，对应于一个单位荷载来提供任何填充结构密度的临界载荷。 此外，固体壳也可能在宏观失稳之前发生局部屈曲，但是，本次实验中并没有建模研究壳的屈曲强度。 3.结果 本节主要介绍数值分析结果以及实验验证，结果误差的分析将在第4节中讨论。根据本文优化问题的定义，所有的力和位移值都只考虑中心加载点(顶部中心节点)的垂直分量。

3.1.等效材料参数
本文通过实体模型确定等效的泊松比ν0和杨氏模量E0值，来校核数值模型的准确性。假设材料为线弹性，那么泊松比和杨氏模量均为常数。根据结构优化理论，多孔材料的等效刚度遵循HS上边界[方程(2)]。 已知SEBS的泊松比接近0.5，与橡胶类似。为了评价结果对泊松比的敏感程度，将只改变泊松比的两个相同Abaqus模型运行两次，其中泊松比分别为ν0= 0.45和 ν0 = 0.49。数值结果表明，单位荷载可导致0.1%的位移偏差，该值可表征材料模量E0的不确定程度；第一阶临界荷载相差1%。这些不确定性相比与其他误差来源来说很小。这一结果证实了本文前面的假设，拓扑 优化结构中的构件主要处于纯拉伸或纯压缩状态，泊松比的影响较小。因此，可以简单地认为泊松比ν0 = 0.49。 通过拟合数值模型的响应和实验观察到的响应可确定等效的杨氏模量E0。

图3(b)所示为实验测得的力–位移曲线，可以发现当力加载至5.3 N时，它是完全线性的。测得六个测量点的一系列值，用它们的平均值拟合曲线，误差线表征了测得的最小值和最大值。曲线斜率为aexp = 5.92 N·mm–1 ，相关系数的平方值为R2 = 0.9999。虽然采用了几何非线性梁模型进行数值模拟，但在数值模型中也发现了实验中观测到的线性关系。因此，杨氏模量可以用单个力–位移关系来获得，如采用单位载荷(1 N)，通过关系 (1 N)/uout = aexp 就可以得到，其中，uout是输出位移。拟合得到弹性模量为E0 = 65.7 MPa，考虑到材料的黏弹性性质，这个数字在制造商提供的邵氏硬度值的预期范围内。  图3(b)中的力–位移图还包含了多孔结构的数据。与上述实验一致，可再次测得到一个近乎完美的线性关系(R2 = 0.9994)。多孔结构的测量刚度比实体结构的小10 %，多孔结构的数值分析预测刚度比测量值小6 %，将在第4节中将具体讨论这种微小偏差的原因。需要注意的是，以上模型的数据参考实际打印结构，即选取27%固体填充刚度，而不是固体结构20%的刚度。使用20 %刚度的数值分析表明，填充结构刚度比实体结构小23 %。

3.2.屈曲性能
对于屈曲分析，需要注意本文是采用二维模型进行优化设计，因而主要研究其面内性能。在第4节中将讨论数值计算和实验结果之间的偏差。 屈曲分析的结果汇总在图5中。在图5(a)中给出了数值计算的实心结构与多孔结构的非线性响应，并与实验得到的屈曲荷载进行比较。 对于实心结构，数值屈曲载荷为23.2 N，实验的屈曲载荷为29.8 N，其对应的模态形状是非常相似的，如图5(d)和(g)所示。跟预想的一样，初始屈曲最早出现在具有最高正应力的杆件处，即图4(d)中的箭头所指杆。 对于多孔结构，平面外扭曲模态出现在一阶面内屈曲模态之前，屈曲载荷为64.3 N，如图5(f)所示。在三维数值模型中，面外的平移运动仅在单一荷载节点处受到约束，因而扭曲模态也能出现。该模型在施加荷载为65.1 N时，平面外模态与实测结果几乎相同，如图5(a)和(c)所示。本文中只考虑多孔结构的一阶面内屈曲模态。为了防止扭转模态出现，两个U形钢支架安装在多孔试样周围，U形钢架间距15.2 mm，如图5(e)所示。支架腿宽19 mm，对称地放置在距离组件65 mm远的中心平面处。一个10 kg的总重量(对应100 N)被附在传递荷载的钢支架处，支架与试样之间的摩擦并未明显约束试样的位移形变。如图5中所示，多孔结构在这个荷载下仍然没有屈服，但是开始发生局部失稳。因此，我们无法从实验获得面内失稳模态，利用数值模型可以确定面内屈曲载荷为126 N。尽管多孔结构和实体结构的拓扑不同，但是多孔结构的失稳模态应与实心结构类似，如图5(b)所示。面外扭曲模态的实验和数值结果的高度一致，也进一步支持了数值模型给出的面内屈曲荷载的正确性。 由于面外屈曲模态并未做任何改变只是简单地由人为从二维模型拉伸为三维模型，需要比较面内屈曲模态所对应的临界荷载值来衡量结构的稳定性。数值计算表明，多孔结构的屈曲载荷比实心模型高5.4倍。图中也给出了面外扭曲模态的临界荷载值，进一步表明数值方法确定面内屈曲临界荷载的正确性。

3.3.填充密度的影响
图6比较了之前提出的两种结构和另外三种结构的数值性能，这三种结构使用相同的参数和体积约束优化获得，但是填充结构的刚度分别为40%、60%和80%，这些刚度值对应的体积密度分别为67%、82%和92%，如图6(a)。图6(b)中屈曲荷载值采用二维线性屈曲分析确定，具体见第2.6节中所述。二维模型给出的屈曲荷载值略低；但是，由固体结构和多孔结构三维模型所得的面内屈曲荷载的比值非常接近(5.3而不是5.4的比例)。 显然，当使用较低的填充率时，结构刚度降低接近线性(即柔度增加)，屈曲荷载显著提高。上述结果的原因在于填充率较低的结构组件宽度较大，弯曲刚度按照组件中心轴垂直距离的三次方增加；而填充比例较高时，对改进屈曲荷载影响很小，其原因是HS边界曲线有一个陡峭的斜坡，接近ρ= 1，这意味着牺牲刚度只能收益一点点梁的宽度。值的注意得是，与其他四种结构相比，λE = 0.2的结构拓扑不同不是高屈曲载荷的唯一原因：采用不同的初始猜测会导致与实体结构相同的拓扑结构。结果表明这种结构的临界荷载比实心结构高4.5倍，而不是5.3倍。 图6(b)给出了填充结构稳定性极限的估计曲线，它是填充密度的函数，具体方法如2.7节中所述。只有相对屈曲荷载低于这条曲线，宏观结构才有望保持局部稳定性。

4. 讨论
实验结果清楚地验证了我们的假设，与传统最小柔顺性方法的优化结构相比，利用涂层方法的增材制造多孔填充结构，屈曲荷载显著提高，但是依旧存在一些小的偏差。 数值模型的两个主要假设分别是SEBS是线弹性材料，以及填充结构是均匀的和各向同性的，且满足HS上边界。

关于假设打印出的SEBS是线弹性材料，这是一种简单的理想化，至少忽略了三方面的影响：
①黏弹性，在整个试验中可以观测到；
②杨氏模量的应变相关性，在276Author name et al. / Engineering 2(2016) xxx–xxx 给定的简单实验装置中与黏弹性的影响难于区分；
③熔 丝打印技术产生的各向异性。 黏弹性表现在依赖于应力–应变曲线的变形率和蠕变的显著程度。

用图3(a)所示的手动设置过程确定E0，为了使蠕变效应最小，施加瞬态荷载。用图3(c)中的实验机测量屈曲荷载，随后获得了采用较高应变速率(5 mm·min–1)的测量值和使用手动设置得到的测量值。前者的屈曲荷载稍高，后者的屈曲荷载稍低，但两者测量结果是在报告值的10%以内。基于这些观察，黏弹性效应被认为是试验中不确定性的主要来源。如果不考虑这些影响，测量位移将会趋于过高估计，也就是说拟合的杨氏模量和数值确定的屈曲载荷将会估计偏低，这或许可以解释测量实心结构的一些偏差。


杨氏模量的应变相关性可在两个方向上影响数值计算结果。正如前面提到的，这部分的影响很难在现有的简单测试装置情况下与黏弹性影响区分开。 熔丝打印技术产生的各向异性对于多孔结构的影响可以忽略，因为打印结构表面和内部填充时，挤出路径都是平行于局部传力路径方向，并且只有两个对于表面和一个用于填充结构的熔丝路径。因此，各向异性主要与实体结构相关；然而，这种影响在框架结构中不是关键的，因为打印所有的杆件时都是相同的挤出模式：界面采用纵向打印方式，而内部采用横向打印方式。在某种意义上，由于各向异性，固体结构可看作是内部刚度稍微降低的全密度的多孔材料，然而，刚度减少的估计超过了本文的范围，并且一定程度上，这种影响可通过拟合等效材料参数进行有效处理。需要注意的是，当拓展到三维设计时，各向异性的影响在垂直于层的方向比所述层的平面内更强。然而，各向异性的程度取决于增 材制造技术的选择，并且在一定程度上可以通过热处理缓解。


此外，假设填充结构是均质、各向同性，且满足HS上边界是相当简化的。在实际打印时，三角形填充结构显然不能从宏观结构尺度上分离出来，这就意味着填充结构应当认为是单独的结构部件，而不是均质材料。但是，例如在参考文献，当使用某些单胞时，认为其均质性质可看成是一个合理的近似。此外，沿着结构边缘填充三角几何形状较差，从而导致材料分布不均匀，进而导致结果不精确，特别是对于比较细的结构。内部填充结构的缺陷性是多孔结构的各向异性的重要来源。 数值模型和实验之间两个最重要的偏差是实体结构的数值屈曲荷载比实验值小约20%，多孔结构的数值预测刚度比实验测量刚度低6%。相比于屈曲荷载的5倍的增加量，这些偏差较小，无法改变该研究的结论。显然，屈曲荷载的增益是涂层方法所固有的，并不限制于MBB梁。图6(b)是填充结构的稳定性曲线，可看出20%的填充刚度结构可能是接近于最佳的局部和全局屈曲强度之间的权衡。然而需要注意的是，这个说法是建立在总体假设的前提下，不仅仅涉及所述的填充结构的稳定性曲线的推导假设，还有该填充结构是均匀、各向同性，且满足HS上限的假设。


5. 结论
本文证实了涂层拓扑优化方法可以用于设计增材制造填充构件，并且相比标准最小柔顺性方法进行优化的结构，屈曲性能大幅度提高。在特定的MBB梁算例中，多孔结构的屈曲荷载比实体结构高5倍以上，而刚度下降20 %~25%。屈曲载荷提高的原因是由于多孔结构的引入增大了结构部件的宽度，进而提升了比弯曲刚度。 增材制造实验试件的实验测试结果支撑了数值结果。然而测量结果与一些不确定性相关，主要是关于实验试件材料属性的粗略假设，但结果的趋势十分明显，与性能差异相比，这种偏差很小。 继本文研究工作之后，探讨一下未来的相关工作。多孔组件能够得到如此高的面内屈曲荷载，这是由于结构发生局部屈曲(由于集中荷载)而不是整体失稳。一般来说，如果填充密度选择过低，在结构发生整体失稳之前，填充结构处会出现局部屈曲。因此，今后我们需要设计完全的三维结构实验，研究多孔材料的面外屈曲性能。 本文研究证明了拓扑优化方法适合增材制造这一特殊制造手段的可能性，也彰显了该方法对于显著改善结构性能的重要作用。