基于变密度法的多材料与结构一体拓扑优化研究*

董 莉，吴晓明

(厦门大学 机电工程系，福建 厦门 361102)

0 引 言

结构拓扑优化是一种结构概念设计方法，它可以在给定约束条件的情况下，根据优化目标，通过优化算法计算得到结构的材料最优布局。拓扑优化也是一种非常有效的结构减重优化方法，它可以得到比传统的形状优化和尺寸优化更轻的结构，因此成为工程界研究热点。

自从1988年BENDSØE等[1]提出了求解连续体结构拓扑优化的均匀化方法以来，针对拓扑优化的问题得到了广泛的研究。依据描述方式的不同, 目前常见的拓扑优化方法包括均匀化方法、变密度法[2-3]、渐进结构优化方法[4]，以及水平集方法[5]等。其中，变密度法是以单元密度作为设计变量，具有概念简单、设计变量少、收敛速度快、程序实现简单的优点。

较为成熟的研究大多是单一材料结构拓扑优化，然而单一材料结构在工程实际中通常难以满足一些需要特定综合性能的工程结构问题，因此进行多材料的结构设计研究非常必要，梯度材料以及在梯度材料上的结构设计优化受到工程界广泛关注。

多材料结构拓扑优化的研究始于THOMSEN[6]，之后一些学者用多种不同性能的材料和空洞进行了多相材料拓扑优化设计。SIGMUND和TORQUATO[7-8]针对两相实体材料和空洞组成的柔性机构拓扑，以及具有极端热膨胀特性的微结构拓扑优化进行了研究；YIN等[9]用单个变量描述了多相材料条件下的材料属性，提出了峰值函数多相材料插值模型；GAO等[10]采用线性对等混合材料插值模型，进行了多相材料结构拓扑优化设计；ZUO等[11]采用序列幂函数插值模型进行了多相材料连续体结构拓扑优化。

另外有一些学者对功能梯度材料的拓扑优化进行了研究。陆丹和刘毅[12]采用双向渐进拓扑优化(BESO)方法,引入了许用应力,以实际应力与许用应力之比作为单元增删的标准，结合梯度单元法实现了对功能梯度材料的拓扑优化；XIA和WANG[13]基于水平集方法，以体分比和结构边界为设计变量,实现了功能梯度结构的材料性能和拓扑布局的并行优化设计；邱克鹏和张卫红[14]采用了凸规划求解策略以及周长控制方法,实现了功能梯度MBB梁和功能梯度夹层结构中夹芯的拓扑构型设计拓扑优化。

这一类研究是在给定的功能梯度材料之上进行的结构设计，而结构需要依据材料的状况而定。现代3D打印的多材料融合技术可以让两种甚至多种材料按照需求，以任意的配比进行融合，且材料可以依据结构的功能要求而变化，从而打破了结构对材料的局限。因此，研究多材料与结构一体设计优化方法具有重要意义。

基于变密度法，笔者提出一种双材料与结构一体拓扑优化模型。

1 材料与结构一体优化模型

1.1 变密度法拓扑优化模型

变密度法将材料密度值设置为在区间[0,1]之间连续变化的变量，其中，密度值为0时表示单元材料不存在，密度值为1时表示单元材料存在。

假设材料密度与材料弹性模量之间的关系如下：

E(xi)=Emin+xi(E0-Emin)xi∈[0,1]

(1)

式中：xi—第i个单元的密度；E(xi)—第i个单元的材料弹性模量；E0—材料存在时的弹性模量；Emin—材料不存在时的弹性模量，为防止出现刚度阵奇异而无法求解，可令Emin=E0/1 000。

由于在结构中出现的中间密度单元是不存在的，为了解决中间密度值单元，笔者引入惩罚插值函数来解决中间单元问题。SIMP法是一个幂函数作为惩罚函数，使得所有单元的密度值惩罚到0或者1附近。

SIMP法模型的表达式如下：

(2)

式中：p—惩罚因子；n—结构单元总数。

引入SIMP插值模型后的材料弹性模量如下：

(3)

引入SIMP法插值模型后的结构的刚度矩阵、柔度函数以及敏度函数如下所示：

(4)

式中：K—整体刚度矩阵；ki—单元刚度矩阵；C(x)—结构柔度值；F—载荷向量；U—位移向量；C′(x)—结构灵敏度。

以体积约束下柔度最小化(即刚度最大化)建立基于SIMP插值模型的变密度法拓扑优化数学模型如下所示：

(5)

式中：Ω—结构初始设计区域；n—离散后的单元总数目；xi—设计变量，表示离散后每个单元的密度值；ui—单元位移向量；vi—每个单元初始体积；V*—结构总体体积约束。

为避免求解时出现总刚度矩阵奇异，取xmin=0.001，xi取[xmin，1]之间的连续值。

1.2 材料与结构一体优化模型

本文提出一种基于变密度法的双材料与结构一体拓扑优化模型。结构的每一个有限元单元可由A，B两种材料以任意比例融合而成，材料A的弹性模量为Ea，材料B的弹性模量为Eb，设单元i中材料B的比例为αi，则单元i中材料A的比例为1-αi，则单元i的材料为Bαi+A(1-αi)。

基于SIMP插值模型的结构与材料一体弹性模量如下：

(6)

刚度矩阵和柔度函数为：

(7)

(8)

结构柔度最小化作为目标函数，以结构的体积和材料B在总体结构中所占比例作为约束，构建基于SIMP法插值模型的结构与材料一体拓扑优化模型如下式所示：

(9)

式中：Ω—结构初始设计区域；xi—离散后每个单元的密度值；αi—单元中材料B的占有比例；n—结构离散后的单元总数；C—结构柔度；F—外力向量；U—位移；Eα—材料A的弹性模量；Eb—材料B的弹性模量；K—结构刚度阵；V*—结构目标体积；μ—拓扑完成结构中材料B总的使用比例。

为了避免计算过程中刚度矩阵奇异，用Emin、xmin代替0；取Emin=0.001，xmin=0.001。

2 敏度分析与过滤

2.1 敏度分析

结构的体积对单元密度的灵敏度为：

(10)

结构的柔度对单元密度的灵敏度为：

(11)

结构的柔度对单元中材料B的比例灵敏度为：

(12)

2.2 过滤技术及优化算法

为避免结果出现棋盘格现象和网格依赖性问题，笔者采用敏度过滤技术，即将目标单元和邻接单元进行加权平均，将得到的敏度值代替原来的值，从而实现过滤。

结构的柔度对单元密度的灵敏度过滤如下：

(13)

(14)

式中：rmin—过滤半径；dist(i,f)—目标单元和邻接单元的中心距。

笔者采用OC优化算法求解式(9)优化问题。OC优化算法推导过程简单、收敛速度快，是常用的优化算法。通过引入拉格朗日乘数因子，构造拉格朗日函数。根据当设计变量取极值时拉格朗日函数应满足K-T必要条件，以及拉格朗日函数关于设计变量和乘数因子的驻值条件，可得到设计变量迭代格式(中间推导过程不再赘述)：

(15)

(16)

式中：m—移动限制；η—阻尼系数，可以在保证收敛稳定性的同时提升迭代效率。

Bi，Di如下式：

(17)

(18)

式中：λ1，λ2—拉格朗日乘子，可由二分法计算得到。

迭代终止条件为两次迭代之间设计变量差值的绝对值小于等于一个值change，即定义终止条件为max(|xk+1-xk|)≤change，其中：k—迭代次数。

基于OC法的多材料与结构一体拓扑优化方法流程如图1所示。

图1 基于OC法的多材料与结构一体拓扑优化方法流程图

3 算例及结果分析

笔者选取两个算例进行数值运算，验证提出设计方法的有效性。两个算例均使用笔者所提出新的模型计算，进行结构与材料一体拓扑优化。

将结构设计区域离散为四节点平面应力单元，每个单元的相对密度值为xmin至1之间的连续值，每个单元由A、B两种材料构成；单元中B材料的比例为αi，A材料的比例则为1-αi，αi—0到1之间的连续变量。以每个单元的相对密度xi和单元中材料B的比例αi为设计变量，以体积和拓扑完成结构中材料B总的使用比例μ为约束，运用OC优化算法进行求解。

3.1 悬臂梁

悬臂梁结构设计域及边界条件如图2所示。

图2 悬臂梁结构设计域及边界条件

图2中，设计域尺寸为40×20，悬臂梁结构左端固定，受力作用点集中在右侧中间的一个单元节点上，初始结构设计域离散成40×20个四节点平面应力单元。材料A的弹性模量Ea=0.5，材料B的弹性模量Eb=1。体积分数为50%，过滤半径为1.5，惩罚因子为3，μ为0.4，移动限制为0.2；阻尼因子为1/2，设计变量初始值x为0.5，α为0.4。

根据笔者提出结构与材料拓扑优化算法，得出的悬臂梁优化结果如图3所示。

图3 悬臂梁优化结果

图3(a)为计算所得的拓扑结构形状，颜色为黑色的单元为保留单元，白色单元为去掉单元，灰色单元为密度值处于中间的单元，该类单元可根据要求保留或删去；

图3(b)为计算所得材料B在拓扑形状基础上每个单元中的比例值αi的示意图，颜色为黑色的单元值αi为1，代表该类单元由B材料构成；颜色为灰色的单元由A、B两种材料按照单元中B材料的比例为αi，A材料的比例为1-αi混合构成。

在图3(a)中保留单元基础上的图3(b)白色单元值为0，代表此类单元由A材料构成。

迭代终止条件为change=0.001时，计算历时20.687 825 s，迭代339步，柔度值为82.817 3。

由优化结果可以看出，笔者提出的优化方法所得拓扑结构形状边界清晰、分布合理，整体优化稳定。在材料受力单元附近与结构上下边界受力变形大的部分，都采用了弹性模量较大的B材料，表明得到的材料比例分布也是十分合理的。

3.2 MBB梁

MBB梁结构设计域及边界条件如图4所示。

图4 MBB梁结构设计域及边界条件

图4中，设计域尺寸为40×20，悬臂梁结构左端固定，受力作用点集中在一个单元节点上，初始结构设计域离散成40×20个四节点平面应力单元。

材料A的弹性模Ea=0.5，材料B的弹性模量Eb=1；体积分数为50%，过滤半径为1.5，惩罚函数为3，μ为0.4，移动限制为0.2，阻尼因子为1/2，设计变量初始值x为0.5，α为0.4。

根据笔者提出的结构与材料拓扑优化算法，得出的MBB梁优化结果如图5所示。

图5 MBB梁优化结果

图5(a)为计算所得的拓扑结构形状，黑色的单元为保留单元，白色单元为去掉单元，灰色单元为密度值处于中间的单元，此类单元可根据要求保留或删去；

图5(b)为计算所得在拓扑结构形状的基础上B材料比例值，黑色单元由B材料构成，灰色单元由A、B两种材料按照单元中B材料的比例为αi，A材料的比例为1-αi混合构成，在图5(a)中保留单元基础上的图5(b)白色单元由A材料构成。

迭代终止条件为change=0.001时，计算历时16.914 929 s，迭代266步，柔度值为93.509 0。

由优化结果可以看出，笔者提出的优化方法所得拓扑结构形状边界清晰、分布合理，整体优化稳定；在材料受力单元附近与结构下边界受力变形大的部分，都采用了弹性模量较大的B材料，表明得到的材料比例分布也是十分合理的。

4 结束语

基于SIMP法插值模型，笔者建立了单元中两种材料不同比例与单元弹性模量之间的关系，在变密度法中，以单元密度作为设计变量的基础上，加入材料在单元构成中所占的比例为设计变量，以结构柔度最小化为目标，结构体积和其中一类材料在结构中所占比值为约束，建立了基于变密度法的结构与材料一体拓扑优化方法模型，并采用OC优化算法进行了求解。

相关算例及研究结果表明，该方法可以计算出结构的拓扑构型，同时可得到两种材料比例数值结果；且拓扑构型边界清晰，可以实现结构与材料一体的协同优化。