周期性多孔结构特征值拓扑优化

付君健， 张 跃， 杜义贤， 高 亮

(1. 水电机械设备设计与维护湖北省重点实验室， 湖北 宜昌 443002；2. 三峡大学 机械与动力学院， 湖北 宜昌 443002； 3. 华中科技大学 数字制造装备与技术国家重点实验室， 武汉 430074)

周期性多孔结构是功构一体化的优良载体[1]，在航空航天、汽车、医学植入等领域具有广泛的应用前景。周期性多孔结构可设计性强，能根据功能特性需求，设计出不同形状、尺寸和孔隙率的多孔构型。对于含有多孔结构的复杂机电系统，多孔结构的固有频率要求远离驱动系统的工作频率。

在给定的材料属性和边界条件下，周期性多孔结构的固有频率与拓扑构型相关，通过改变多孔结构的拓扑构型，可以被动地改变结构固有频率，使其脱离工作环境下的激励频率，从而避免剧烈共振导致的结构失稳和破坏。特征值拓扑优化是改变多孔结构固有频率的有效设计方法，与宏观结构特征值拓扑优化[2]的不同之处在于，周期性多孔结构通过改变局部单胞的几何特征，实现整体结构特征值的优化。

由于存在尺度跨越，特征值有限元求解方法对周期性多孔结构拓扑优化有重要影响。对于具有相同自由度数量的结构系统，求解频率和振型的特征值问题，其计算时间将比静力分析高出一个数量级[3]。如果是局部几何特征较多的多孔结构，计算量将更大。由于尺度小、特性多、计算量大，多孔结构宏观等效属性计算一般需采用模型缩减方法。均匀化法是常用的周期性多孔结构宏观等效属性计算方法[4]，主要用于计算宏观等效弹性矩阵，等效弹性矩阵则用于计算等效刚度矩阵和等效质量矩阵[5]。均匀化法的假设之一为尺度分离，即多孔微结构尺寸远远小于宏观结构尺寸。在实际应用中，均匀化法更适合小尺度、大规模的周期性多孔结构分析与优化。对于有尺度要求的多孔结构设计，需采用尺度关联的模型缩减方法进行宏观等效属性计算，如子结构法[6]、尺度关联均匀化法[7]、多尺度有限元法[8]，或者采用迭代法进行求解[9]。

子结构法作为一种较为成熟的模型缩减方法，近年来逐渐被用于细观尺度周期性多孔结构的拓扑优化[10]。子结构法具备尺度关联的特性，适用于尺度关联、中等规模多孔结构分析与优化。目前，基于子结构法的多孔结构拓扑优化研究主要集中于静力学问题[11-12]，实现周期性和多层级多孔结构刚度拓扑优化。虽然已有研究将子结构法用于多层级结构一阶特征值拓扑优化[13]，并取得了一定效果，但子结构静态凝聚求解高阶特征值会产生较大误差。因为子结构静态凝聚忽略了动力学方程中的惯性项，导致缩减的质量矩阵为近似值，只有在频率较低时才有较高的缩减精度，随着频率的增加，数值误差将逐渐增大。O’Callahan[14]提出了针对动力学问题的改进缩减系统(improved reduced system, IRS)，保证了子结构法的动态缩减精度，但采用子结构法进行多孔结构特征值分析与优化仍存在诸多问题有待探索，如集中质量施加、特征值尺度问题、特征值拓扑优化敏度分析等。

本文基于子结构动态凝聚方法，探讨其动态缩减精度，建立了周期性多孔结构特征值最大化拓扑优化模型，实现了尺度关联的二维、三维周期性多孔结构的拓扑优化，并揭示了细观尺度多孔结构特征值拓扑优化的尺度效应。

1 周期性多孔结构几何与物理描述

1.1 周期性多孔结构几何描述

针对周期性多孔结构几何描述，提出局部水平集函数(local level set function, LLSF)的概念，即在全局设计域中定义多个局部水平集函数描述细观尺度下的周期性多孔结构，将多孔结构几何边界隐式地嵌入高一维的函数。局部水平集函数使得多孔结构的几何描述更加简单，也使得多孔结构拓扑优化具有更多的设计自由[15]。如图1所示，将二维周期性结构设计域离散为M个细观结构，在每个细观结构内部定义一个局部水平集函数，Φi表示其中一个细观结构的水平集函数

(a) 零水平面

(1)

式中：下标i为设计域D内的子结构编号；x为空间坐标向量；Ω为结构域；∂Ω为结构边界。

所有细观结构水平集函数Φi(xi)的集合组成了全局水平集函数Φ(x)

(2)

局部水平集函数Φi在子结构内部动态演化的Hamilton-Jacobi方程如下

(3)

式中：vn, i为局部水平集函数的法向速度场；t为时间变量。

为克服传统水平集方法的数值问题，便于与梯度优化算法结合，本文采用紧支径向基函数插值替代离散的水平集函数，形成参数化水平集法[16]。在参数化水平集法中引入局部水平集函数，可有效避免初始扩展系数求解困难的问题[17-18]。

(4)

式中：φ为紧支径向基函数插值矩阵；α为扩展系数向量，其中，k=1,…,N，N为细观结构设计域紧支径向基插值控制点的数量。

根据式(4)，参数化之后的局部水平集函数法向速度场为

(5)

1.2 周期性多孔结构物理描述

本文采用子结构法[19]描述多孔结构物理模型，如图2所示。将多孔结构单胞凝聚为具有多个节点的超单元。与周期性多孔结构的几何描述定义类似，子结构法的基本思想是将计算域Ω分解为若干子区域Ωi，先求解边界上的数值信息，然后对子区域内部“各个击破”。在有限元分析中采用子结构法可以实现多孔结构物理模型和几何模型的匹配，缩小有限元分析的计算规模，降低计算内存的消耗，提高周期性多孔结构特征值分析与拓扑优化的效率。

图2 子结构凝聚

特征值问题的有限元平衡方程为一个无阻尼自由振动系统，为了展示子结构法动态凝聚也适用于考虑载荷的动力学问题，在特征值问题中引入载荷向量F=0，则特征值问题平衡方程的矩阵形式为

(6)

(7)

由式(7)第二行可得主节点位移和从节点位移之间的关系

(8)

对于静态凝聚，忽略式(8)中所有的惯性项

KsmUm+KssUs=Fs

假设Kss为正定矩阵，且结构内部的从节点不受外载荷作用，即Fs=0

(10)

由式(10)得到变化关系

(11)

式中，Tc为静态凝聚的变换矩阵

(12)

根据式(11)和式(12)，缩减后的平衡方程为

(13)

其中

(14)

(15)

由于式(9)忽略了全部惯性项，对于特征值问题，式(13)的缩减精度存在较大误差，需要对此进行改进。由式(13)可得

(16)

对式(10)求二阶导数

(17)

将式(16)代入式(17)

(18)

将式(16)和式(18)代入式(8)

(19)

对于无阻尼自由振动问题，主、从节点上的载荷Fs和Fm均为零，式(19)可进一步简化为

(20)

由此得到子结构动态凝聚的变换矩阵

(21)

其中

经过子结构动态凝聚后的无阻尼自由振动平衡方程为

(23)

其中，缩减后的质量矩阵和刚度矩阵分别为

(24)

(25)

单一结构的凝聚可通过式(23)实现，周期性多孔结构则涉及多个子结构的凝聚。由于周期性排布的假设条件，含有重复结构特征的子结构只需进行一次凝聚，多个重复子结构的凝聚过程可通过对单个凝聚子结构的组装完成。由于载荷向量一般取零，多个重复子结构组装而成的缩减平衡方程为

(26)

(27)

(28)

为验证子结构动态凝聚求解特征值问题的精度，与子结构静态凝聚、全自由度有限元模型进行对比分析。如图3所示为对比分析的模型，尺寸为1 m×1 m二维正方形板，厚度为0.1 m，左下角和右下角固定约束，离散为40×40个双线性四边形单元，单元尺寸为0.025 m×0.025 m。材料弹性模量为200 GPa，泊松比为0.3，质量密度为7 800 kg/m3。

(a) 全自由度模型

表1为有限元分析的特征频率值，与全自由度有限元分析对比可知，子结构静态凝聚求解的特征值误差较大，随着阶数的增大，特征频率误差越来越大。子结构动态凝聚求解的特征值误差较小，且误差不会随着阶数的增加而放大。因此，子结构动态凝聚能保证特征值求解的精度，同时能缩减结构自由度数量，提高特征值问题的求解效率。

表1 特征频率对比

2 周期性多孔结构特征值拓扑优化模型

多孔结构的共振很大程度取决于前几阶频率，剧烈的共振一般发生在外界激励频率和某一阶固有频率接近的时候，所以在多孔结构设计时，必须对结构的前几阶固有频率进行限制或改变。本文以前几阶加权特征值最大化为目标函数，构造如式所示的p-norm形式的函数，并以此缓解模态转换引起的目标函数振荡问题[20]。该目标函数可通过改变特征值权重ωj和p-norm函数的p值，灵活调整对全部或特定阶次特征值的影响。

(29)

(30)

(31)

(32)

式中：ε为应变场；D为实体材料的弹性矩阵；ρ为材料质量密度。

由于特征值λ是针对多孔结构整体的特征值，优化模型式(30)中的目标函数并未写成对每个多孔结构积分或累加的形式，但灵敏度表达式可写成对子结构积分或累加的形式，其具体形式见灵敏度分析。

3 灵敏度分析与数值实施

3.1 灵敏度分析

采用形状导数来推导目标函数和约束条件关于设计变量(扩展系数)的灵敏度。根据形状导数的定义及其引理[21]，目标函数、能量双线性形式和载荷双线性形式关于时间变量t的导数分别为：

(33)

(34)

(35)

将式(30)中的平衡方程两边对时间t求偏导得到

(36)

(37)

将式(34)、式(35)、式(37)代入式(36)可得

(38)

通过构造伴随方程可得，该优化问题为自伴随问题[22]，即v=u。

(39)

根据质量矩阵的正交归一化条件[23]

b(u,u)=1

(40)

式(38)可进一步简化为

(41)

式(41)仍存在边界积分形式，引入式(42)

dΓ=δ(Φ)|∇Φ|dΩ

(42)

式(41)的边界积分转换为体积积分

|∇Φ|vndΩ

(43)

由于本文采用子结构法，式(43)可写成对每个子结构积分的形式

(44)

γi,j(ui,λj)=εT(ui,j)Dε(ui,j)-λjρiui,jui,j

(45)

式中：下标i为子结构的编号；下标j为特征值的阶次。

将式(5)中的局部水平集函数的速度场vn,i代入式(44)可得

(46)

将特征值关于时间t的导数代入式(33)可得

(47)

通过链式法则对目标函数J(Φ)直接求其关于时间变量t的偏导数可以得到

(48)

对比式(47)和式(48)可得目标函数关于扩展系数α的导数

φi,k(xi)dΩidΩ

(49)

同理，体积约束对扩展系数α的导数为

(50)

3.2 数值实施

结构特征值的求解是在缩减的全局有限元模型上实现，该过程暂不涉及子结构内部的扩展。在求解敏度信息时，需要求解子结构内部的特征向量。子结构内部特征向量求解可顺序求解或者CPU并行求解，CPU并行求解在子结构划分数量较多时能节约大量时间。在求解时为避免重复计算，在子结构动态缩减时，可将式(20)中相关的数据重复使用。此外，在敏度计算时，子结构内部特征向量的求解可重复使用子结构动态缩减中式(21)的值，这样TIRS只需要计算一次。

为防止优化过程中随着结构材料的不断删除，出现无限多个特征值，需要在结构的特定位置加入所谓的非结构集中质量。如果不加入非结构集中质量，在结构的特定位置就不可能获得期望的最优拓扑结构，从而影响到整个结构频率的最大化[24]。在周期性多孔结构设计中，单方向的周期数量可能为奇数或者偶数，非结构集中质量的施加位置视具体情况而定，非结构集中质量可施加在如图4(a)所示的子结构内部自由度上，也可施加在如图4(b)所示的缩减结构的自由度上。加在子结构内部时，需在子结构内部质量矩阵自由度所在位置的对角线上添加一个很大的数值，大小约为允许设计质量的10%左右，非结构集中质量会参与子结构的动态凝聚。加在缩减的结构自由度上时，需在子结构凝聚和组装完成之后，在缩减结构的质量矩阵自由度所在位置的对角线上添加设计质量10%左右的数值。

(a) 子结构内部

4 数值案例与结果分析

4.1 悬臂梁结构

如图5所示的二维悬臂梁结构，梁结构设计域长宽比L∶H= 2 m∶1 m，厚度为0.1 m，将宏观结构设计域划分为20×10个子结构，子结构内部划分为40×40个双线性四边形单元，单元的尺寸为0.025 m×0.025 m，材料弹性模量E=200 GPa，泊松比v= 0.3，质量密度7 800 kg/m3，在宏观结构右边界中点处施加大小为设计质量10%的集中质量，以多孔结构体积分数f= 0.5为约束，进行二维周期性多孔结构设计。

图5 悬臂梁设计域

图6所示为经过优化的子结构单胞、水平集函数及其组装而成的宏观结构，经过162次迭代后，其最优的目标函数值收敛到16.103，体积分数收敛到0.5，前6阶频率分别为0.398 Hz、31.621 Hz、55.398 Hz、60.04 Hz、94.633 Hz、110.121 Hz。图7所示为前6阶特征频率对应的特征向量图，即结构的模态振型。值得注意的是，数值算例中周期性多孔结构一阶固有频率较低，通过数值试验分析可知，引入非结构集中质量是导致该现象的主要原因。因为在整个宏观设计域中仅定义了子结构的局部水平集函数，因此只需要对一个子结构单胞进行求解，即可得到如图6(c)所示的20×10个周期性分布的子结构。从图6(c)可以看出，最优的周期性结构在集中质量施加点处存在材料，这是由于子结构动态缩减的尺度关联性所带来的益处。

(a) 最优子结构单胞

(a)

在迭代初期前6阶频率一直处于下降趋势，这是因为初始化的结构体积分数大于指定的体积分数，迭代初期的约束条件不满足，导致前6阶频率较大，如图8所示。当体积分数约束得到满足后，随着优化过程的进行，前6阶频率逐渐收敛。

图8 前6阶频率迭代图

4.2 两端固定梁结构

为验证本方法的高效性和特征值问题的尺度效应，考虑如图9所示的两端固定的二维梁结构，将本方法与没有模型缩减的全自由度拓扑优化[25]的计算时间进行了对比。梁结构设计域长宽比L∶H= 7 m∶1 m，厚度为0.1 m，将宏观结构设计域分别划分为14×2、28×4、56×8和84×12个子结构，子结构内部划分为40×40个双线性四边形单元，单元的尺寸为0.025 m×0.025 m，材料弹性模量E=200 GPa，泊松比v= 0.3，质量密度7 800 kg/m3。在结构中心施加大小为设计质量10%的集中质量，以子结构体积分数f= 0.5为约束，进行二维周期性多孔结构设计。

图9 两端固定二维梁设计域

当宏观结构离散为14×2、28×4、56×8和84×12个子结构时，其结构自由度数量分别约为9万、36万、144万和323万，采用子结构法分别将其缩减至约为0.57万、2万、7.6万、17万。如表2所示为对应4种子结构数量下的拓扑优化结果。其中最优子结构拓扑构型的特征在于上下两端有横向分布的主承载结构，主承载结构之间有交叉的结构。随着子结构数量的增多，最优子结构的拓扑构型在局部展现出一定的变化，其目标函数也随之下降。图10所示对应4种子结构数量下周期多孔结构特征频率的变化趋势，在固定大小的设计域内，当离散的子结构数量越多时，优化结果表现出一定的尺度效应，其最优多孔结构的前6阶特征频率整体呈下降趋势。因此，在实际工程应用中，可

图10 多孔结构特征频率变化

表2 特征值拓扑优化结果

针对结构具体的力学性能需求来设置多孔结构周期性排布的数量，从而避开结构的激振频率。

在拓扑优化的计算效率方面，如图10所示，对比本文子结构法与全自由度模型，当子结构数量为14×2时，在每一步迭代中，两种方法在有限元分析和拓扑优化的平均计算时间相当。当子结构数量大于14×2后，在每一步迭代中，本文方法在有限元分析和拓扑优化方面的平均计算时间大幅下降。当子结构数量为84×12时，本文方法优化更新的计算时间仅为全自由度模型的19.4%，有限元分析的计算时间仅为全自由模型的31.9%，平均每步的迭代时间约为全自由度模型的30.5%。其原因在于，本文采用了子结构动态凝聚缩减了自由度数量，采用了局部水平集函数描述多孔结构降低了几何描述上的复杂度，从几何模型和物理模型角度共同提升了周期性多孔结构特征值拓扑优化的效率。此外，从图11可知，有限元分析占据了拓扑优化大部分的计算时间，研究高效的有限元分析方法将有助于提升拓扑优化的计算效率。

图11 计算效率对比

4.3 三维梁结构

如图12所示的两端固定的三维梁结构，结构设计域长宽高比例为L∶W∶H= 7 m∶1 m∶1 m，将结构设计域划分为14×2×2个三维子结构，每个子结构离散为10×10×10个八节点六面体单元，单元的尺寸为0.1m× 0.1 m × 0.1 m，材料弹性模量E= 200 GPa，泊松比v= 0.3，质量密度7 800 kg/m3。在结构中心施加大小为设计质量10%的集中质量，以多孔结构体积分数f= 0.4为约束，进行三维周期性多孔结构设计。

三维结构的有限元分析和子结构的凝聚相对于二维问题具有更大的计算量，本算例采用三维子结构法将结构自由度从18.6万左右缩减至约4.2万。如图13所示为经过优化的三维周期性结构及其子结构单胞，从子结构单胞拓扑构型可知，组装后的周期性多孔结构在集中质量点处必然存在材料，这是子结构法尺度关联所带来的益处。从图14可知，拓扑优化模型经过93次迭代后，其最优的目标函数值收敛至36.564 7，体积分数收敛至0.4。图15所示为前6阶频率迭代图，分别收敛至1.558 Hz、36.049 Hz、68.633 Hz、82.604 Hz、82.814 Hz、106.298 Hz。

图13 三维周期性多孔结构优化结果

(a)

图15 前6阶特征频率迭代图

5 结 论

子结构法是高效、高精度的模型缩减方法，本文将参数化水平集方法和子结构动态凝聚相结合，研究了周期性多孔结构的特征值拓扑优化方法。研究表明，本方法能有效实现尺度关联的二维和三维周期性多孔结构的特征值拓扑优化，且所设计的周期性多孔结构具有一定的尺度效应。