TOPOLOGY OPTIMIZATION OF MULTI-MATERIAL STRUCTURES FOR DYNAMIC STIFFNESS USING POLYGONAL MULTIRESOLUTION SCHEME
-
摘要:
利用多边形有限单元的高精度求解优势，融合多分辨率拓扑优化方法，实现粗糙位移网格条件下的高分辨率构型设计，由此提出多材料结构动刚度问题的拓扑优化方法。将多边形单元(位移场求解单元)劈分为精细的小单元，构造设计变量与密度变量的重叠网格，形成多分辨率-多边形单元的优化建模策略；以平均动柔度最小化为目标和多材料的体积占比为约束，建立多材料结构的动力学拓扑优化模型，通过HHT-α方法求解结构动响应，采用伴随变量法推导目标函数和约束的灵敏度表达式，利用基于敏度分离技术的ZPR设计变量更新方案构建多区域体积约束问题的优化迭代格式；通过典型数值算例分析优化方法的可行性和动态载荷作用时间对优化结果的影响机制。
Abstract:Polygonal elements with high precision are used in finite element modeling along with a multiresolution scheme for topology optimization to achieve high resolution design with a coarse finite element mesh. A polygonal multiresolution scheme is proposed to perform a topology optimization of multi-material structures for dynamic stiffness. This combined modelling strategy of polygonal finite element with higher resolution density and design discretization is extended to optimal multi-material structural design problem, which is based on the concept that a polygonal element for dynamic analysis is split into fine design variable elements for optimization and overlapped density variable elements for representation of material distribution. To minimize the mean dynamic compliance subjected to the volume fraction constraints of each candidate material, a topology optimization model of multi-material structure is established. The HHT-α method is adopted to solve the structural dynamic problem. Following the sensitivity analysis of objective function and constraints by the adjoint variable method, the ZPR design variable update scheme is employed to solve highly nonlinear and non-convex optimization problem with multiple regional volume constraints. Several benchmark numerical examples are presented to analyze the feasibility of the proposed algorithm, and the influence of the duration time of dynamic load on the optimized results is analyzed.
-
在飞机、汽车制造与结构工程领域，多材料设计兼顾每一相基材料的有点，能够保证结构具有综合的力学性能，因而在飞机、汽车与桥梁结构中的应用日益增多[1-3]。拓扑优化方法能够寻求满足性能约束条件下材料的最优布局，极大地丰富了多材料结构的设计手段。与传统单一材料的拓扑优化相比，多相材料的拓扑优化极大地扩大了设计空间，能够提供综合性能更优的最优解。随着多材料增材制造技术的发展，依据最优拓扑的计算模型，使通过逐点和逐层方式打印多材料设计结构成为可能[4-6]。由此，多材料的拓扑优化问题引起了研究机构的广泛关注。根据设计变量的定义方式，多材料拓扑优化方法通常归纳为2类：基于单元的方法[7-18]和基于边界的方法[19-23]。在基于单元的方法框架内，ZUO等[7]、杜义贤等[8]、闫浩和吴晓明[9]通过构造材料性能的序列插值模型，实现了柔性机构和传热机构的多材料最优布局设计。俞燎宏等[10]和朱本亮等[11]采用交替主动相算法，将多相优化布局问题转变为序列两相材料优化子问题，然后通过并行计算策略逐级分层计算。目前，多体积约束下的多材料布局优化已拓展至稳态热传导[12]、热力耦合[13-14]、频域动力学[15]和空间桁架结构[16]的最优布局问题。最近，YANG和LI[17]、HUANG和LI[18]则提出了单一质量约束下的多材料轻量化设计方法。在基于边界方法的框架内，水平集法[19-21]和移动构件法(MMC)[22-23]因其能使边界光滑清晰、便于提取设计构型等独特优势，也受到了一定的关注与研究。但是，现有的基于边界的方法未能充分考虑新孔洞的创建，优化结果过度依赖于初始设计，难于获得优化问题的全局最优解。
拓扑优化框架要求在消耗尽可能小的计算成本下获得高分辨率的优化结果，而上述目标的实现取决于有限元求解器、自由度数量、优化建模策略、材料插值模型以及后处理等诸多因素。SUKUMAR[24]首次提出了多边形有限元法，数值实验表明它非常适于求解具有复杂几何结构区域的连续介质力学问题。将多边形单元应用于拓扑优化问题，能够显著地减小棋盘格和孤岛效应等数值奇异性问题[25-28]。为了减少计算成本与提高求解精度，NGUYEN-XUAN[29]构建了多边形有限元法与自适应网格重划技术融合策略，CHAU等[30]基于上述研究工作求解了多材料结构的静态拓扑优化问题。另外，NGUYEN等[31]提出了高分辨率拓扑优化方法 (multiresolution topology optimization，MTOP)，它采用多层级网格优化建模策略，即利用粗糙网格完成有限元分析，精细的重叠网格描述设计变量和密度变量空间，从而形成高分辨率的拓扑优化结果。最近，FILIPOV等[32]将多边形单元替代传统单元作为分析单元，精确地估计复杂结构的动力学响应，高效地实现了特征值、受迫振动等结构频域动力学问题的高分辨率拓扑优化。
在结构时域动响应方面，以提高结构动刚度[27, 33]，降低结构振动幅度[34-35]和动应力幅值[28]等拓扑优化方法相继提出。SHOBEIRI[33]在敏度分析中忽略了动能对单元删除的影响，难于获得惯性力影响显著的动力学拓扑优化解。ZHAO和WANG[34]、GIRALDO-LONDOÑO等[27-28]分别采用了先微分-后离散和先离散-后微分的伴随敏度分析方法，但是，后者基于多边形单元计算结构动响应，能够处理任意曲线边界结构的动力学优化问题。针对模型减缩技术在动力学拓扑优化过程中的有效性问题，ZHAO和WANG[35]指出模态加速法在减少响应计算规模和预测精度上的平衡能力优于模态位移法，然而，如果高频模态为主要影响模态时，模型减缩技术将失效。
多相材料时域动力学布局优化问题属于强非自伴随问题，需要考虑多材料结构的动态约束条件，致使敏度分析和优化问题求解具有挑战性，因而目前的多相材料布局优化集中于静态优化问题，而多相材料的时域动力学布局优化问题研究相对较少。因此，本文将多边形网格与MTOP的融合方法拓展至多材料结构的动态拓扑优化问题，建立融合框架内的多相材料高分辨率布局的时域动力学优化模型。为了避免先微分-后离散方法引起的灵敏度计算的相容性误差，采用先离散-后微分的伴随敏度分析方法。通过泰勒公式实现近似描述目标函数和体积约束函数，建立原问题的凸规划子模型，采用基于敏度分离技术的ZPR方法求解子问题。最后，通过数值算例检验提出的方法在多相材料动力学布局优化的可行性。
1 多分辨率-多边形网格建模策略
多边形单元在动力学分析方面其求解精度优于传统单元，因而将其作为动力学优化问题的响应计算单元[25-28]。同时，将多边形单元劈分为精细的小单元形成设计变量与密度变量网格，并使两者具有相同的位置坐标，继而构建多分辨率-多边形单元的优化建模策略，能够实现粗糙位移网格条件下的高分辨率的拓扑优化构型设计(如图1所示)。多分辨率-多边形单元建模策略，利用粗糙的多边形网格精确求解位移场，精细的设计变量与密度变量重叠网格用于构型优化设计，将有效地提升优化计算效率与优化结果的品质[26]。
在MTOP框架内，为了精确计算多边形位移单元的单元刚度矩阵，将形函数及其梯度的积分点设置在密度变量所在位置(也是设计变量网格的中心)。对于多材料问题，单元刚度和质量矩阵分别表示为：
{{\boldsymbol{k}}_\ell } \cong \sum\limits_{i = 1}^{{N_n}} {\sum\limits_{j = 1}^m {{{\tilde m}_E}({y_{\ell ,ij}}){\boldsymbol{B}}_\ell ^{\text{T}}{{\boldsymbol{D}}^0}{{\boldsymbol{B}}_\ell }{|_i}{A_{\ell ,i}}} } (1) {{\boldsymbol{m}}_\ell } \cong \sum\limits_{i = 1}^{{N_n}} {\sum\limits_{j = 1}^m {{{\tilde m}_V}({y_{\ell ,ij}}){\rho _0}{{\boldsymbol{N}}_\ell ^{e{\text{T}}}}{{\boldsymbol{N}}_{\ell }^e}{|_i}{A_{\ell ,i}}} } (2) 式中：{N_n}为多边形位移单元积分点的个数；{y_{\ell ,ij}}为密度变量在积分点处的估计值；{{\boldsymbol{B}}_\ell }为应变矩阵，{{\boldsymbol{D}}^0}为线弹性材料本构矩阵；{A_{\ell ,i}}为设计变量网格或密度变量网格的面积； {\rho _0} 为材料密度；{\tilde m_E}({y_{\ell ,ij}})和 {\tilde m_V}({y_{\ell ,ij}}) 为m种材料的刚度与体积差值函数。
基于门槛投影函数[36]，体积插值函数表示为：
{\tilde m_V}({y_{\ell ,ij}}) = \varepsilon + (1 - \varepsilon ){m_V}({y_{\ell ,ij}}) (3) {m_V}({y_{\ell ,ij}}) = \frac{{\tanh (\bar \beta \bar \eta ) + \tanh (\bar \beta ({y_{\ell ,ij}} - \bar \eta ))}}{{\tanh (\bar \beta \bar \eta ) + \tanh (\bar \beta (1 - \bar \eta ))}} (4) 式中： \bar \eta 为门槛密度； \bar \beta 为投影控制参数； \varepsilon \ll 1 为Ersatz数，用于阻止{y_{\ell ,ij}} \to 0时引起的数值失稳。
刚度插值函数{\tilde m_E}({y_{\ell ,ij}})是由密度惩罚函数{m_W}\left( {{y_{\ell ,ij}}} \right)与多材料插值函数{m_M}\left( {{y_{\ell ,ij}}} \right)合成形成。首先，通过惩罚函数将单元密度变量推向0或1，基于RAMP材料插值模型构造密度惩罚函数为[37]：
{w_{\ell ,ij}} = {m_W}\left( {{y_{\ell ,ij}}} \right) = \frac{{{m_V}({y_{\ell ,ij}})}}{{1 + {p_0}(1 - {m_V}({y_{\ell ,ij}}))}} (5) 式中，{p_0} > 0为惩罚参数。然后，构造多材料刚度插值函数，表示为[38]：
{\tilde m_E}({y_{\ell ,ij}}) = \varepsilon + (1 - \varepsilon )\sum\limits_{j = 1}^m {{w_{\ell ,ij}}\prod\limits_{\scriptstyle j = i\atop \scriptstyle j \ne i} ^m {\left( {1 - {w_{\ell ,ij}}} \right)E_i} } (6) 式中，E_i为第i种材料的弹性模量。
由此，根据上述单元刚度与质量矩阵，多边形位移单元的总体矩阵表示为：
{\boldsymbol{M}} = \sum\limits_{\ell = 1}^N {{{\boldsymbol{m}}_\ell }} \;{\text{, }}\;{\boldsymbol{K}} = \sum\limits_{\ell = 1}^N {{{\boldsymbol{k}}_\ell }} (7) 为了抑制棋盘格与孤岛现象，利用线性滤波方法获得网格无关的优化结果，则有：
{y_{\ell ,ij}} = {{\sum\limits_{n \in {S_i}} {{d_{\ell ,nj}}w({{\boldsymbol{x}}_n} - {{\boldsymbol{x}}_i})} } \mathord{\left/ {\vphantom {{\sum\limits_{n \in {S_i}} {{d_{\ell ,nj}}w({{\boldsymbol{x}}_n} - {{\boldsymbol{x}}_i})} } {\sum\limits_{n \in {S_i}} {w({{\boldsymbol{x}}_n} - {{\boldsymbol{x}}_i})} }}} \right. } {\sum\limits_{n \in {S_i}} {w({{\boldsymbol{x}}_n} - {{\boldsymbol{x}}_i})} }} (8) 式中：{S_i}为相应于密度变量单元i所占的子域；{{\boldsymbol{x}}_n}为与设计变量{d_{\ell ,nj}}对应的设计变量网格的中心坐标。线性权重函数定义为：
w({x_n} - {x_i}) = \left\{ \begin{aligned} & {{{\left( {{r_{\min }} - {r_{ni}}} \right)} / {{r_{\min }}}}},\;\;{{r_{ni}} \leqslant {r_{\min }}} \\ & 0,\qquad\qquad\qquad\;\;{{\text{其他}}} \end{aligned}\right. (9) 式中：{r_{ni}}为单元密度单元i和n的中心距；{r_{\min }}为指定的过滤半径。通过式(8)将设计变量线性加权处理而映射为密度变量，以此代入式(1)和式(2)中即可获得单元刚度与质量矩阵。
此外，根据式(8)，密度变量对于设计变量的灵敏度为：
{{\partial {y_{l,ij}}} \mathord{\left/ {\vphantom {{\partial {y_{l,ij}}} {\partial {d_{\ell ,nj}}}}} \right. } {\partial {d_{\ell ,nj}}}} = {{w({{\boldsymbol{x}}_n} - {{\boldsymbol{x}}_i})} \mathord{\left/ {\vphantom {{w({{\boldsymbol{x}}_n} - {{\boldsymbol{x}}_i})} {\sum\limits_{n \in {S_i}} {w({{\boldsymbol{x}}_n} - {{\boldsymbol{x}}_i})} }}} \right. } {\sum\limits_{n \in {S_i}} {w({{\boldsymbol{x}}_n} - {{\boldsymbol{x}}_i})} }} (10) 2 多材料结构动刚度拓扑优化
2.1 优化模型
以平均动柔度最小化为目标和多材料的体积占比为约束，建立多材料结构的动力学拓扑优化模型。假设设计域中包含m种材料，按弹性模量高低线性排列，在有限材料约束下，其动刚度优化模型为：
\left\{ \begin{aligned} & \mathop {\min }\limits_{{\boldsymbol{d}} \in {{[0,1]}^{N \times {N_n} \times m}}} f({\boldsymbol{d}},{{\boldsymbol{u}}_0},\cdots,{{\boldsymbol{u}}_{{N_{\text{t}}}}}) = {1 / {{N_{\text{t}}}}} \cdot \sum\limits_{i = 0}^{_{{N_{\text{t}}}}} {{{\boldsymbol{f}}_{i}^{\text{T}}}{{\boldsymbol{u}}_i}} \\[-2pt]& {\text{ s}}.{\text{t}}.{\text{ }}{g_j}({\boldsymbol{d}}) = {{\sum\limits_{\ell \in {\varepsilon _j}} {\sum\limits_{i \in {\eta _j}} {\sum\limits_{r \in {\chi _j}} {{A_{\ell ,i}}{{\tilde m}_V}({y_{\ell ,ir}})} } } } \Bigg/ {\sum\limits_{\ell \in {\varepsilon _j}} {\sum\limits_{i \in {\eta _j}} {{A_{\ell ,i}}} } }} - \\[-2pt]&\qquad {{\bar v}_j} \leqslant 0,j = 1,2,\cdots,{N_{\text{c}}} \\[-2pt]& {\boldsymbol{M}}{{{\boldsymbol{\ddot u}}}_i} + {\boldsymbol{C}}{{{\boldsymbol{\dot u}}}_i} + {\boldsymbol{K}}{{\boldsymbol{u}}_i} - {{\boldsymbol{f}}_i} = {\boldsymbol{0}},\;\;i = 1,2,3,\cdots,{N_{\text{t}}} \end{aligned} \right. (11) 式中：{\boldsymbol{d}}为设计变量矢量；{N_{\text{t}}}为动态激励作用时间定分的区间数；{N_{\text{c}}}为多材料体积约束的数量；{{\boldsymbol{f}}_i}为t = {t_i}时的动载荷向量； {u}_{i}、{\dot{u}}_{i}、{\ddot{u}}_{i} 分别为相应的结构位移、速度、加速度响应； {\boldsymbol{C}} = {\alpha _r}{\boldsymbol{M}} + {\beta _r}{\boldsymbol{K}} 为阻尼矩阵( {\alpha _r} 、 {\beta _r} 为瑞利阻尼参数)；{\varepsilon _j}、{\eta _j}、{\chi _j}分别为单元指标集、设计变量指标集和多材料相数指标集， {\bar v_j} 为第j种材料占设计域总体积的体积分数。
2.2 HHT-α方法
HHT-α方法作为广义的Newmark-β方法，它通过数值阻尼参数α调控动力学控制方程的积分格式，使迭代过程无条件稳定，非常适于求解结构动力学问题。根据文献[39]，HHT-α方法将优化模型式(11)中的半离散形式的有限元方程修改为：
\begin{split} & {\boldsymbol{M}}{{{\boldsymbol{\ddot u}}}_i} + (1 - \alpha ){\boldsymbol{C}}{{{\boldsymbol{\dot u}}}_i} + \alpha {\boldsymbol{C}}{{{\boldsymbol{\dot u}}}_{i - 1}} + (1 - \alpha ){\boldsymbol{K}}{{\boldsymbol{u}}_i} + \alpha {\boldsymbol{K}}{{\boldsymbol{u}}_{i - 1}} =\\&\qquad (1 - \alpha ){{\boldsymbol{f}}_i} + \alpha {{\boldsymbol{f}}_{i - 1}},i = 1,2,\cdots,{N_t} \end{split} (12) 通过Newmark-β有限差分关系，位移、速度场的更新格式为：
{{\boldsymbol{u}}_i} = {{\boldsymbol{u}}_{i - 1}} + \Delta t{{\boldsymbol{\dot u}}_{i - 1}} + \Delta {t^2}\left[\left(\frac{1}{2} - \beta \right){{\boldsymbol{\ddot u}}_{i - 1}} + \beta {{\boldsymbol{\ddot u}}_i}\right] (13) {{\boldsymbol{\dot u}}_i} = {{\boldsymbol{\dot u}}_{i - 1}} + \Delta t[(1 - \gamma ){{\boldsymbol{\ddot u}}_{i - 1}} + \gamma {{\boldsymbol{\ddot u}}_i}] (14) 式中，\beta {\text{ = }}{\left( {1{\text{ + }}\alpha } \right)^2}/4，\gamma {\text{ = }}\left( {1{\text{ + }}2\alpha } \right){\text{/2}}为算法参数，合理选择参数\alpha 保证算法具有至少二阶精度和无条件稳定。
将式(13)、式(14)代入式(12)，得到离散形式的控制方程的残差，则有：
\begin{split} & {{\boldsymbol{R}}_i} = {{\boldsymbol{M}}_1}{{{\boldsymbol{\ddot u}}}_i} + {{\boldsymbol{M}}_0}{{{\boldsymbol{\ddot u}}}_{i - 1}} + {{\boldsymbol{C}}_0}{{{\boldsymbol{\dot u}}}_{i - 1}} + {\boldsymbol{K}}{{\boldsymbol{u}}_{i - 1}} -\\&\qquad (1 - \alpha ){{\boldsymbol{f}}_i} - \alpha {{\boldsymbol{f}}_{i - 1}} = {\boldsymbol{0}} \end{split} (15) \left\{ \begin{aligned} & {{\boldsymbol{M}}_1} = {\boldsymbol{M}} + (1 - \alpha )\gamma \Delta t{\boldsymbol{C}} + (1 - \alpha )\beta \Delta {t^2}{\boldsymbol{K}} \\[-2pt]& {{\boldsymbol{M}}_0} = (1 - \alpha )(1 - \gamma )\Delta t{\boldsymbol{C}} + (1 - \alpha )\left(\frac{1}{2} - \beta \right)\Delta {t^2}{\boldsymbol{K}} \\[-2pt]& {{\boldsymbol{C}}_0} = {\boldsymbol{C}} + (1 - \alpha )\Delta t{\boldsymbol{K}} \end{aligned} \right. (16) 对于HHT-α残差方程式(15)，其初始时刻的表达式为：
{{\boldsymbol{R}}_0} = {\boldsymbol{M}}{{\boldsymbol{\ddot u}}_0} + {\boldsymbol{C}}{{\boldsymbol{\dot u}}_0} + {\boldsymbol{K}}{{\boldsymbol{u}}_0} - {{\boldsymbol{f}}_0} = {\boldsymbol{0}} (17) 由此，求解式(31)获得加速度 {{\boldsymbol{\ddot u}}_i} ，然后回代至式(13)和式(14)，更新位移 {{\boldsymbol{u}}_i} 和速度 {{\boldsymbol{\dot u}}_i} 。
2.3 敏度分析
根据优化模型式(11)，在MTOP框架内目标函数对于设计变量的灵敏度为：
\frac{{{\rm d}f}}{{{\rm d}{d_{\ell ,nj}}}} = \frac{{\partial f}}{{\partial {d_{\ell ,nj}}}} + \sum\limits_{i = 1}^{Nt} {\frac{{\partial f}}{{\partial {{\boldsymbol{u}}_i}}}} \frac{{\partial {{\boldsymbol{u}}_i}}}{{\partial {y_{\ell ,nj}}}}\frac{{\partial {y_{\ell ,nj}}}}{{\partial {d_{\ell ,nj}}}} (18) 为了避免状态变量对于设计变量的导数计算，使用伴随变量法完成目标函数的敏度分析。类似于HHT-α残差方程式(15)，将Newmark-β有限差分关系式(13)、式(14)变换为残差形式：
\begin{split} & {{\boldsymbol{P}}_i} = - {{\boldsymbol{u}}_i} + {{\boldsymbol{u}}_{i - 1}} + \Delta t{{{\boldsymbol{\dot u}}}_{i - 1}} + \\&\qquad \Delta {t^2}\left[\left(\frac{1}{2} - \beta \right){{{\boldsymbol{\ddot u}}}_{i - 1}} + \beta {{{\boldsymbol{\ddot u}}}_i}\right] = {\boldsymbol{0}} \end{split} (19) {{\boldsymbol{Q}}_i} = - {{\boldsymbol{\dot u}}_i} + {{\boldsymbol{\dot u}}_{i - 1}} + \Delta t[(1 - \gamma ){{\boldsymbol{\ddot u}}_{i - 1}} + \gamma {{\boldsymbol{\ddot u}}_i}] = {\boldsymbol{0}} (20) 由此，引入伴随变量{{\boldsymbol{\xi }}_i}、{{\boldsymbol{\mu }}_i}、{{\boldsymbol{\nu }}_i}， i = 1, 2,\cdots,{N_{\text{t}}} ，原目标函数修改为拉格朗日函数\bar f({\boldsymbol{d}},{{\boldsymbol{u}}_0},{{\boldsymbol{u}}_1},\cdots,{{\boldsymbol{u}}_{{N_{\text{t}}}}})，其灵敏度表示为：
\begin{split} & \frac{{{\rm d}\bar f}}{{{\rm d}{d_{\ell ,nj}}}} = \frac{{{\rm d}f}}{{{\rm d}{d_{\ell ,nj}}}} + \sum\limits_{i = 0}^{Nt} {{{\boldsymbol{\xi }}_i^{\text{T}}}\frac{{{\rm d}{{\boldsymbol{R}}_i}}}{{{\rm d}{d_{\ell ,nj}}}}} + \\&\qquad \sum\limits_{i = 1}^{{N_{\text{t}}}} {{{\boldsymbol{\mu }}_i^{\text{T}}}\frac{{{\rm d}{{\boldsymbol{P}}_i}}}{{{\rm d}{d_{\ell ,nj}}}}} + \sum\limits_{i = 1}^{{N_{\text{t}}}} {{{\boldsymbol{\upsilon }}_i^{\text{T}}}\frac{{{\rm d}{{\boldsymbol{Q}}_i}}}{{{\rm d}{d_{\ell ,nj}}}}} \end{split} (21) 根据式(19)、式(20)，{{\partial {{\boldsymbol{P}}_i}} / {\partial {d_{\ell ,nj}}}} = 0和{{\partial {{\boldsymbol{Q}}_i}} / {\partial {d_{\ell ,nj}}}} = 0以及初始条件{{\partial {{\boldsymbol{u}}_0}}/ {\partial {d_{\ell ,nj}}}} = 0和{{\partial {{{\boldsymbol{\dot u}}}_0}} / {\partial {d_{\ell ,nj}}}} = 0，同时，满足伴随方程式(23)和式(24)使{{\partial {{{\boldsymbol{\ddot u}}}_0}}/ {\partial {d_{\ell ,nj}}}}、{{\partial {{{\boldsymbol{\ddot u}}}_i}} / {\partial {d_{\ell ,nj}}}}、{{\partial {{{\boldsymbol{\dot u}}}_i}} / {\partial {d_{\ell ,nj}}}}和{{\partial {{\boldsymbol{u}}_i}} / {\partial {d_{\ell ,nj}}}}等项消失，则式(21)化简为：
\frac{{{\rm d}\bar f}}{{{\rm d}{d_{\ell ,nj}}}} = \frac{{\partial f}}{{\partial {d_{\ell ,nj}}}} + \sum\limits_{i = 1}^{Nt} {{\boldsymbol{\xi }}_i^{\text{T}}\frac{{\partial {{\boldsymbol{R}}_i}}}{{\partial {y_{\ell ,nj}}}}\frac{{\partial {y_{\ell ,nj}}}}{{\partial {d_{\ell ,nj}}}}} \;\; (22) {{\boldsymbol{\xi }}_0^{\text{T}}}\frac{{\partial {{\boldsymbol{R}}_0}}}{{\partial {{{\boldsymbol{\ddot u}}}_0}}} + {{\boldsymbol{\xi }}_1^{\text{T}}}\frac{{\partial {{\boldsymbol{R}}_1}}}{{\partial {{{\boldsymbol{\ddot u}}}_0}}} + {{\boldsymbol{\mu }}_1^{\text{T}}}\frac{{\partial {{\boldsymbol{P}}_1}}}{{\partial {{{\boldsymbol{\ddot u}}}_0}}} + {{\boldsymbol{\upsilon }}_1^{\text{T}}}\frac{{\partial {{\boldsymbol{Q}}_1}}}{{\partial {{{\boldsymbol{\ddot u}}}_0}}} = 0 (23) \left\{ \begin{gathered} \sum\limits_{\ell = 1}^{Nt} \left( {{\boldsymbol{\xi }}_\ell ^{\text{T}}}\frac{{\partial {{\boldsymbol{R}}_\ell }}}{{\partial {{\boldsymbol{u}}_i}}} + {{\boldsymbol{\mu }}_\ell ^{\text{T}}}\frac{{\partial {{\boldsymbol{P}}_\ell }}}{{\partial {{\boldsymbol{u}}_i}}} + {{\boldsymbol{\upsilon }}_\ell ^{\text{T}}}\frac{{\partial {{\boldsymbol{Q}}_\ell }}}{{\partial {{\boldsymbol{u}}_i}}} + \frac{{\partial f}}{{\partial {{\boldsymbol{u}}_i}}}\right) = 0 \\[-2pt] \sum\limits_{\ell = 1}^{Nt} \left( {{\boldsymbol{\xi }}_\ell ^{\text{T}}}\frac{{\partial {{\boldsymbol{R}}_\ell }}}{{\partial {{{\boldsymbol{\dot u}}}_i}}} + {{\boldsymbol{\mu }}_\ell ^{\text{T}}}\frac{{\partial {{\boldsymbol{P}}_\ell }}}{{\partial {{{\boldsymbol{\dot u}}}_i}}} + {{\boldsymbol{\upsilon }}_\ell ^{\text{T}}}\frac{{\partial {{\boldsymbol{Q}}_\ell }}}{{\partial {{{\boldsymbol{\dot u}}}_i}}} + \frac{{\partial f}}{{\partial {{{\boldsymbol{\dot u}}}_i}}}\right) = 0 \\[-2pt] \sum\limits_{\ell = 1}^{Nt} \left( {{\boldsymbol{\xi }}_\ell ^{\text{T}}}\frac{{\partial {{\boldsymbol{R}}_\ell }}}{{\partial {{{\boldsymbol{\ddot u}}}_i}}} + {{\boldsymbol{\mu }}_\ell ^{\text{T}}}\frac{{\partial {{\boldsymbol{P}}_\ell }}}{{\partial {{{\boldsymbol{\ddot u}}}_i}}} + {{\boldsymbol{\upsilon }}_\ell ^{\text{T}}}\frac{{\partial {{\boldsymbol{Q}}_\ell }}}{{\partial {{{\boldsymbol{\ddot u}}}_i}}} + \frac{{\partial f}}{{\partial {{{\boldsymbol{\ddot u}}}_i}}}\right) = 0 \\[-2pt] \end{gathered} \right. (24) 将HHT-α残差方程式(15)、初始条件式(17)和Newmark-β有限差分关系式(13)、式(14)代入伴随方程式(23)、式(24)，则有：
{{\boldsymbol{\mu }}_{{N_{\text{t}}}}} = \frac{{\partial f}}{{\partial {{\boldsymbol{u}}_{{N_{\text{t}}}}}}}\;,\;{{\boldsymbol{\upsilon }}_{{N_{\text{t}}}}} = {\boldsymbol{0}}\;,\;{{\boldsymbol{M}}_1}{{\boldsymbol{\xi }}_{{N_{\text{t}}}}} = - \beta \Delta {t^2}{{\boldsymbol{\mu }}_{{N_{\text{t}}}}} - \gamma \Delta t{{\boldsymbol{\upsilon }}_{{N_{\text{t}}}}} (25) {{\boldsymbol{\mu }}_{i - 1}} = \frac{{\partial f}}{{\partial {{\boldsymbol{u}}_{i - 1}}}} + {\boldsymbol{K}}{{\boldsymbol{\xi }}_i} + {{\boldsymbol{\mu }}_i}\;,\;{{\boldsymbol{\upsilon }}_{i - 1}} = {{\boldsymbol{C}}_0}{{\boldsymbol{\xi }}_i} + \Delta t{{\boldsymbol{\mu }}_i} + {{\boldsymbol{\upsilon }}_i} (26) \begin{split} {{\boldsymbol{M}}_1}{{\boldsymbol{\xi }}_{i - 1}} =& {{\boldsymbol{M}}_0}{{\boldsymbol{\xi }}_i} - \Delta {t^2}\left[\beta {{\boldsymbol{\mu }}_{i - 1}} + \left(\frac{1}{2} - \beta \right){{\boldsymbol{\mu }}_i}\right] -\\& \Delta t[\gamma {{\boldsymbol{\upsilon }}_{i - 1}} + (1 - \gamma ){{\boldsymbol{\upsilon }}_i}] \end{split} (27) {\boldsymbol{M}}{{\boldsymbol{\xi }}_0} = {{\boldsymbol{M}}_0}{{\boldsymbol{\xi }}_1} - \left(\frac{1}{2} - \beta \right)\Delta {t^2}{{\boldsymbol{\mu }}_1} - (1 - \gamma )\Delta t{{\boldsymbol{\upsilon }}_1} (28) 由于优化模型式(11)中使用了体积插值函数{\boldsymbol{V}} = {\tilde m_V}\left( {\boldsymbol{y}} \right)与刚度插值函数{\boldsymbol{E}} = {\tilde m_E}\left( {\boldsymbol{y}} \right)，那么通过链式求导法则，拉格朗日函数与约束函数的灵敏度为：
\frac{{{\rm d}\bar f}}{{{\rm d}{d_{\ell ,nj}}}} = \frac{{\partial {y_{\ell ,nj}}}}{{\partial {d_{\ell ,nj}}}}\frac{{\partial {E_{\ell ,nj}}}}{{\partial {y_{\ell ,nj}}}}\frac{{{\rm d}\bar f}}{{{\rm d}{E_{\ell ,nj}}}} + \frac{{\partial {y_{\ell ,nj}}}}{{\partial {d_{\ell ,nj}}}}\frac{{\partial {V_{\ell ,nj}}}}{{\partial {y_{\ell ,nj}}}}\frac{{{\rm d}\bar f}}{{{\rm d}{V_{\ell ,nj}}_i}} (29) \frac{{{\rm d}{g_i}}}{{{\rm d}{d_{\ell ,nj}}}} = \frac{{\partial {y_{\ell ,nj}}}}{{\partial {d_{\ell ,nj}}}}\frac{{\partial {E_{\ell ,nj}}}}{{\partial {y_{\ell ,nj}}}}\frac{{{\rm d}{g_i}}}{{{\rm d}{E_{\ell ,nj}}}} + \frac{{\partial {y_{\ell ,nj}}}}{{\partial {d_{\ell ,nj}}}}\frac{{\partial {V_{\ell ,nj}}}}{{\partial {y_{\ell ,nj}}}}\frac{{{\rm d}{g_i}}}{{{\rm d}{V_{\ell ,nj}}_i}} (30) 利用伴随变量法，拉格朗日函数\bar f({\boldsymbol{d}},{{\boldsymbol{u}}_0} ,{{\boldsymbol{u}}_1}, \cdots , {{\boldsymbol{u}}_{{N_{\text{t}}}}})对于{\boldsymbol{V}} = {\tilde m_V}\left( {\boldsymbol{y}} \right)和{\boldsymbol{E}} = {\tilde m_E}\left( {\boldsymbol{y}} \right)的敏度表示为：
\left\{ \begin{aligned} & {\frac{{{\rm d}\bar f}}{{{\rm d}{E_{\ell ,nj}}}} = \frac{{\partial \bar f}}{{\partial {E_{\ell ,nj}}}}{\text{ + }}\sum\limits_{i = 0}^{{N_t}} {{{\boldsymbol{\xi }}_{i}^{\text{T}}}\frac{{\partial {{\boldsymbol{R}}_i}}}{{\partial {E_{\ell ,nj}}}}} } \\ & {\frac{{{\rm d}\bar f}}{{{\rm d}{V_{\ell ,nj}}}} = \frac{{\partial \bar f}}{{\partial {V_{\ell ,nj}}}}{\text{ + }}\sum\limits_{i = 0}^{{N_t}} {{{\boldsymbol{\xi }}_{i}^{\text{T}}}\frac{{\partial {{\boldsymbol{R}}_i}}}{{\partial {V_{\ell ,nj}}}}} } \end{aligned} \right. (31) 另外，根据优化模型式(11)，目标函数和约束函数对于{\boldsymbol{V}} = {\tilde m_V}\left( {\boldsymbol{y}} \right)和{\boldsymbol{E}} = {\tilde m_E}\left( {\boldsymbol{y}} \right)的敏度表示为：
\frac{{\partial f}}{{\partial {E_{\ell ,nj}}}}{\text{ = }}0\;,\;\frac{{\partial f}}{{\partial {V_{\ell ,nj}}}} = 0\;,\frac{{\partial f}}{{\partial {{\boldsymbol{u}}_i}}} = \frac{1}{{{N_t}}}{{\boldsymbol{f}}_i} (32) \frac{{\partial {g_i}}}{{\partial {E_{\ell ,nj}}}} = 0\;,\frac{{\partial {g_i}}}{{\partial {V_{\ell ,nj}}}} = \frac{{{A_{\ell ,i}}}}{{\displaystyle\sum\limits_{\ell \in {\varepsilon _j}} {\displaystyle\sum\limits_{i \in {\eta _j}} {{A_{\ell ,i}}} } }} (33) 结合式(1)、式(2)和式(7)，HHT-α残差方程对于{\boldsymbol{V}} = {\tilde m_V}\left( {\boldsymbol{y}} \right)和{\boldsymbol{E}} = {\tilde m_E}\left( {\boldsymbol{y}} \right)的偏导数表示为：
\frac{{\partial {{\boldsymbol{R}}_0}}}{{\partial {E_{\ell ,nj}}}} = \frac{{\partial {{\boldsymbol{k}}_\ell }}}{{\partial {E_{\ell ,nj}}}}({{\boldsymbol{u}}_{\ell 0}} + {\beta _r}{{\boldsymbol{\dot u}}_{\ell 0}}) (34) \frac{{\partial {{\boldsymbol{R}}_0}}}{{\partial {V_{\ell ,nj}}_\ell }} = \frac{{\partial {{\boldsymbol{m}}_\ell }}}{{\partial {V_{\ell ,nj}}_\ell }}({{\boldsymbol{\ddot u}}_{\ell 0}} + {\alpha _r}{{\boldsymbol{\dot \upsilon }}_{\ell 0}}) (35) \frac{{\partial {{\boldsymbol{R}}_i}}}{{\partial {E_{\ell ,nj}}}} = \frac{{\partial {{\boldsymbol{k}}_\ell }}}{{\partial {E_{\ell ,nj}}}}[(1 - \alpha )({{\boldsymbol{u}}_{\ell i}} + {\beta _r}{{{\boldsymbol{\dot u}}}_{\ell i}}) + \alpha ({{\boldsymbol{u}}_{\ell ,i - 1}} + {\beta _r}{{{\boldsymbol{\dot u}}}_{\ell ,i - 1}})] (36) \frac{{\partial {{\boldsymbol{R}}_i}}}{{\partial {V_{\ell ,nj}}}} = \frac{{\partial {{\boldsymbol{m}}_\ell }}}{{\partial {V_{\ell ,nj}}}}[{{\boldsymbol{\ddot u}}_{\ell i}} + {\alpha _r}((1 - \alpha ){{\boldsymbol{\dot u}}_{\ell i}} + \alpha {{\boldsymbol{\dot u}}_{\ell ,i - 1}})] (37) 由此，根据式(25)~式(28)确定伴随变量，随后将式(31)~式(37)代入式(29)和式(30)，则可确定拉格朗日函数与约束函数的灵敏度。
2.4 ZPR设计变量更新方案
利用凸近似方法，优化模型式(11)的近似子问题：
\left\{ \begin{aligned} & {\mathop {\min }\limits_{{\boldsymbol{d}} \in [{\underline {{\boldsymbol{d}}}_k},{{{\boldsymbol{\bar d}}}_k}]} \tilde f({\boldsymbol{d}}) = {{\boldsymbol{a}}^{\text{T}}}{{\boldsymbol{d}}^{ - \tilde \alpha }} + {{\boldsymbol{b}}^{\text{T}}}{\boldsymbol{d}}} \\ & {\text{ s}}.{\text{t}}.{\text{ }}{g_i}({\boldsymbol{d}}) = \sum\limits_{\ell ,n,j} {{{\left. {\partial {g_i}/\partial {d_{\ell ,nj}}} \right|}_{{{\boldsymbol{d}}_k}}}( {{d_{\ell ,nj}} - d_{\ell ,nj}^k} )} +\\&\qquad {g_j}({{\boldsymbol{d}}_k}) \leqslant 0\;,\;i = 1,2,\cdots,{N_{\text{c}}} \end{aligned} \right. (38) 式中：\tilde \alpha > 0；\underline {d} _{\ell ,nj}^k = \max ( {d_{\ell ,nj}^{\min },d_{\ell ,nj}^k - {\text{move}}} )和\bar d_{\ell ,nj}^k = \max ( {d_{\ell ,nj}^{\max },d_{\ell ,nj}^k{\text{ + move}}} )是第k次迭代时的上界与下界，move为移动限；{\boldsymbol{a}} \geqslant {\boldsymbol{0}}和{\boldsymbol{b}} \geqslant {\boldsymbol{0}}为系数矩阵。
将目标函数灵敏度分解为正、负两项，结合{\left. {{{\partial \tilde f} / {\partial {d_{\ell ,nj}}}}} \right|_{{{\boldsymbol{d}}^k}}} = {\left. {{{\partial \bar f} / {\partial {d_{\ell ,nj}}}}} \right|_{{{\boldsymbol{d}}^k}}}，则有[40-41]：
{a_{\ell ,nj}}{\text{ = }} - \frac{1}{{\tilde \alpha }}{( {d_{\ell ,nj}^k} )^{1 + \tilde \alpha }}{\left. {\frac{{\partial {f^ - }}}{{\partial {d_{\ell ,nj}}}}} \right|_{{{\boldsymbol{d}}^k}}},{b_{\ell ,nj}}{\text{ = }}{\left. {\frac{{\partial {f^ + }}}{{\partial {d_{\ell ,nj}}}}} \right|_{{{\boldsymbol{d}}^k}}} (39) \left\{ \begin{aligned} & {\frac{{\partial {f^ - }}}{{\partial {d_{\ell ,nj}}}}{\text{ = }}\min \left\{ { - \frac{{\left| {h_{\ell ,nj}^k} \right|d_{\ell ,nj}^k}}{{\tilde \alpha + 1}},\frac{{\partial f}}{{\partial {d_{\ell ,nj}}}}} \right\}} \\ & {\frac{{\partial {f^ + }}}{{\partial {d_{\ell ,nj}}}} = \frac{{\partial f}}{{\partial {d_{\ell ,nj}}}} - \frac{{\partial {f^ - }}}{{\partial {d_{\ell ,nj}}}}} \end{aligned}\right. (40) 式中：h_{\ell ,nj}^k为目标函数的Hessian矩阵，可通过PSB准牛顿法估计[40]。子问题式(38)的拉格朗日函数为：
\begin{split} & L\left( {{\boldsymbol{d}},{\lambda _1},{\lambda _2}, \cdots ,{\lambda _{{N_c}}}} \right) = {{\boldsymbol{a}}^{\text{T}}}{{\boldsymbol{d}}^{ - \tilde \alpha }} + {{\boldsymbol{b}}^{\text{T}}}{\boldsymbol{d}} +\\&\qquad \sum\limits_{i = 1}^{{N_c}} {{\lambda _i}{g_i}\left( {\boldsymbol{d}} \right)} = \sum\limits_{i = 1}^{{N_c}} {{L^i}\left( {{\boldsymbol{d}},{\lambda _i}} \right)} \end{split} (41) \begin{split} & {L^i}\left( {{\boldsymbol{d}},{\lambda _i}} \right){\text{ = }}\sum\limits_{\ell ,n,j} \Bigg[ {a_{\ell ,nj}}d_{\ell ,nj}^{ - \tilde \alpha } + {b_{\ell ,nj}}{d_{\ell ,nj}} +\\&\qquad {\lambda _i}{\left. {\frac{{\partial {g_i}}}{{\partial {d_{\ell ,nj}}}}} \right|_{{{\boldsymbol{d}}^k}}}( {{d_{\ell ,nj}} - d_{\ell ,nj}^k} ) \Bigg] {\text{ + }}{\lambda _i}{g_i}( {{{\boldsymbol{d}}^k}} ) \end{split} (42) 式(42)表明，每一个设计变量仅与一个约束函数相关，因而拉格朗日函数L可采用分离变量的方法求极值。鉴于Li的可分离特性，其一阶最优条件为：
\frac{{\partial {L^i}\left( {{\boldsymbol{d}},{\lambda _i}} \right)}}{{\partial {d_{\ell ,nj}}}}{\text{ = }} - {a_{\ell ,nj}}d_{\ell ,nj}^{ - \tilde \alpha - 1} + {b_{\ell ,nj}} + {\lambda _i}{\left. {\frac{{\partial {g_i}}}{{\partial {d_{\ell ,nj}}}}} \right|_{{{\boldsymbol{d}}^k}}} = 0 (43) d_{\ell ,nj}^ * {\text{ = }}d_{\ell ,nj}^k{\left[ {\frac{{{{\left. {{{ - \partial {f^ - }} /{\partial {d_{\ell ,nj}}}}} \right|}_{{{\boldsymbol{d}}^k}}}}}{{{{\left. {\left[ {{{\partial {f^ + }} /{\partial {d_{\ell ,nj}}}} + {\lambda _i}{{\partial {g_i}} / {\partial {d_{\ell ,nj}}}}} \right]} \right|}_{{{\boldsymbol{d}}^k}}}}}} \right]^{1/\left( {1 + \alpha } \right)}} (44) 式中：d_{\ell ,nj}^ * 为子问题设计变量的最优点。由此，满足盒子约束作用下的设计变量的更新格式为：
{d_{\ell ,nj}} = \left\{ \begin{aligned} & {\underline{d} _{\ell ,nj}^k}\;\;\;\;{{\text{if }}\;\;\;\;d_{\ell ,nj}^ * \leqslant \underline{d} _{\ell ,nj}^k} \\ & {d_{\ell ,nj}^ * }\;\;\;\;{{\text{if }}\;\;\;\;\underline{d} _{\ell ,nj}^k \leqslant d_{\ell ,nj}^ * \leqslant \bar d_{\ell ,nj}^k} \\ & {\bar d_{\ell ,nj}^k}\;\;\;\;{{\text{if }}\;\;\;\;d_{\ell ,nj}^ * \geqslant \bar d_{\ell ,nj}^k} \end{aligned}\right. (45) 将式(45)代入式(42)，获得子问题的对偶目标函数，根据{L^i}\left( {{\boldsymbol{d}}\left( {{\lambda _i}} \right),{\lambda _i}} \right)的相互独立性，寻求拉格朗日乘子{\lambda _i}使其达到极大值，则{L^i}\left( {{\boldsymbol{d}}\left( {{\lambda _i}} \right),{\lambda _i}} \right)的驻点条件表示为：
{{\partial {L^i}} /{\partial {\lambda _i}}}{\text{ = }}\sum\limits_{\ell ,n,j} {[ {g_i}( {{{\boldsymbol{d}}^k}} ) + {\left. {{{\partial {g_i}} / {\partial {d_{\ell ,nj}}}}} \right|_{{{\boldsymbol{d}}^k}}}( {{d_{\ell ,nj}}( {{\lambda _i}} ) - d_{\ell ,nj}^k} ) ]} = 0 (46) 式(46)为关于{\lambda _i}的非线性方程，可通过区间收缩法或迭代法求解。
3 典型算例
本节通过悬臂梁、简支梁、L型结构与Michell型悬臂梁4个典型数值算例，对比4种工况：① 静态载荷；② 载荷持续时间为0.05 s；③ 载荷持续时间为0.03 s；④ 2种~10种不同弹性模量(E_i = {E_0}( m - i + 1 )/m，i = 1,2,\cdots, m，m = 2,3,\cdots, 10，如图2所示)多相材料的优化结果。分析该优化方法的可行性和动态载荷作用时间对优化结果的影响机制。
对于任一种多材料优化工况，每一相材料体积约束的上限{\bar v_j} = 0.45/m\left( {j = 1,2, \cdots ,m} \right)。另外，为了确保HHT-α方法至少具有二阶精度和无条件稳定性，所用算例均使用\alpha {\text{ = }}0.05，\beta {\text{ = }}{\left( {1{\text{ + }}\alpha } \right)^2}/4，\gamma {\text{ = }}\left( {1{\text{ + }}2\alpha } \right){\text{/2}}。
3.1 悬臂梁
如图3所示，悬臂梁右侧自由边中点处承受正弦半波载荷f\left( t \right) = {f_0}\sin \left( {\pi t/{t_{\text{f}}}} \right)，载荷幅值{f_0}{\text{ = }} 1 {\;\text{ kN}}，初始设计域的结构尺寸为L \times H \times h = 8 {\;\text{ m}} \times 4 {\;\text{ m}} \times 0.01 {\;\text{ m}}。多相材料的基准弹性模量{E_0} = 200{\;\text{ GPa}}，任一材料相的泊松比{\nu _i}{\text{ = }}0.3，密度{\rho _i}{\text{ = }} 7800{\;\text{ kg/}}{{\text{m}}^3}，瑞利阻尼参数{\alpha _{\text{r}}}{\text{ = }}10，{\beta _{\text{r}}}{\text{ = }}1 \times {10^{-5}}。
图4为悬臂梁双材料设计的迭代历史(载荷作用时间{t_{\text{f}}} = 0.05{\;\text{ s}})。目标函数渐进收敛于0.0294 N·m，软、硬两相材料均渐进收敛于体积约束的上限：{\bar \nu _1} = 0.225，{\bar \nu _2} = 0.225。
![]()
图 4 悬臂梁双材料设计的迭代历史({t_{\text{f}}}{\text{ = }}0.05{\;\text{ s}})Figure 4. Iteration histories of a bimaterial cantilever beam design ({t_{\text{f}}}{\text{ = }}0.05{\;\text{ s}})图5为不同载荷持续时间工况下双材料悬臂梁问题的优化设计结果。静力学的优化拓扑更接近于动力学{t_{\text{f}}} = 0.05{\;\text{ s}}的结果，但是，动力学{t_{\text{f}}} = 0.03{\;\text{ s}}的优化拓扑显著不同于前两者的结果。与静力学解相比，动力学解分布硬的材料相或三角形结构用于强化载荷作用点处的刚度，以抵制动载荷的快速波动。如图6所示，静力学工况的位移与动柔度响应均高于动载荷工况，从而验证了前面的假设。另外，{t_{\text{f}}} = 0.05{\;\text{ s}}工况下的设计将硬材料相置于载荷作用端，使其产生较小的变形，因而具有更高的动刚度。
![]()
图 5 双材料悬臂梁问题动静态载荷优化结果对比Figure 5. Comparison of dynamic and static load optimization results of bimaterial cantilever beam![]()
图 6 最优悬臂梁载荷作用点处垂向位移和动柔度时间历程Figure 6. Time history of vertical displacement and dynamic flexibility at the loading point of optimal cantilever structure图7为{t_{\text{f}}} = 0.05\;{\text{ s}}时悬臂梁多材料拓扑优化结果。外部拓扑结构基本一致，内部拓扑结构依靠交叉肋结构强化，随着材料相数目的增加，低弹性模量材料向固定端分布，特别是对于4种材料和10种材料工况，固定端增加了强化结构平衡惯性力的影响。
图8为悬臂梁多材料最优拓扑的动柔度时间历程，多相材料设计的动柔度变化基本一致，平均动柔度几乎相同(如图7所示)。由于7相多材料设计的最软相材料位于内部加强肋的交叉点处，使结构的平均动柔度高于其他多相材料设计。
![]()
图 8 悬臂梁多材料最优拓扑的动柔度时间历程Figure 8. Time history of dynamic flexibility for optimized multi-material topology of cantilever beam3.2 简支梁
如图9所示，简支梁上端自由边中点处承受余弦半波载荷f\left( t \right) = {f_0}\cos \left( {\pi t/{t_{\text{f}}}} \right)，载荷幅值{f_0}{\text{ = }}1{\;\text{ kN}}，初始设计域的结构尺寸为：L = 12{\;\text{ m}}，H = 2{\;\text{ m}}，h = 0.01{\;\text{ m}}。多相材料的基准弹性模量{E_0}=200{\;\text{ GPa}}，任一材料相的泊松比{\nu _i}{\text{ = }}0.3，密度{\rho _i}{\text{ = }} 7800{\;\text{ kg/}}{{\text{m}}^3}，瑞利阻尼参数{\alpha _{\text{r}}}{\text{ = }}10，{\beta _{\text{r}}}{\text{ = }}1 \times {10^{-5}}。
图10为不同载荷持续时间工况下双材料简支梁问题的优化设计结果。3种工况的优化拓扑皆是由两侧支撑杆和中间环型结构组成，但是两相材料分布迥然不同。静力学工况的优化拓扑在中间环型结构处集中了高弹性模量材料，因而它的动刚度大于动力学工况优化拓扑的动刚度；但是，动力学工况的优化拓扑在两侧支撑杆集中了高弹性模量材料，能够进一步抑制惯性载荷的影响，致使其动刚度大于静力学工况优化拓扑的动刚度。如图11所示，静力学工况的位移与动柔度响应介于两种动载荷工况优化拓扑的相应动响应之间，验证了前面的结论。另外，随着动力学载荷作用时间的减少，阻尼对于结构的能量耗散作用有所减弱。
![]()
图 10 双材料简支梁问题动静态载荷优化结果对比Figure 10. Comparison of dynamic and static load optimization results of bimaterial simply supported beam![]()
图 11 最优简支梁载荷作用处垂向位移和动柔度时间历程Figure 11. Time history of vertical displacement and dynamic flexibility at the loading point of the optimal simply supported beam structure图12为{t_{\text{f}}} = 0.05{\;\text{ s}}时简支梁的多材料拓扑优化结果，多相材料设计的平均动柔度基本一致。多相材料结构的整体拓扑基本一致，中间环型强化结构有所差异，特别是4种材料、6种材料和10种材料的优化拓扑差异显著。随着材料相数目的增加，中间环型强化结构的下端桁架部分布置更多的低弹性模量材料，两侧支撑杆结构布置更多的高弹性模量材料，使两者能够达到刚度的最优匹配。
![]()
图 12 简支梁的多材料拓扑优化Figure 12. Multi-material topology optimization of simply supported beams3.3 L型结构
如图13所示，L型结构右侧自由边中点处承受正弦半波载荷f\left( t \right) = {f_0}\sin \left( {\pi t/{t_{\text{f}}}} \right)，载荷幅值{f_0}= 1{\;\text{ kN}}，初始设计域的结构尺寸为：{L_1} = {H_2} = 12{\;\text{ m}}，{L_2} = {H_1} = 18{\;\text{ m}}，厚度h = 0.01{\;\text{ m}}。相材料的基准弹性模量{E_0} = 200{\;\text{ GPa}}，任一材料相的泊松比{\nu _i}{\text{ = }}0.3，密度{\rho _i}{\text{ = }}7800{\;\text{ kg/}}{{\text{m}}^3}，瑞利阻尼参数{\alpha _{\text{r}}}{\text{ = }}50，{\beta _{\text{r}}}{\text{ = }} 3 \times {10^{-5}}。
图14为不同载荷持续时间工况下双材料L型结构问题的优化设计结果。动力学工况的优化拓扑在内侧直角应力集中处分布了更多的高弹性模量材料，而静力学工况的优化拓扑在直角拐角处增加了加强筋强化结构。与静力学优化拓扑相比，动力学最优拓扑在载荷作用点附近分布了更多的材料，特别是外侧的三角形构型具有更大的截面尺寸，因而具有更大的动刚度。如图15所示，静力学工况的位移与动柔度响应均高于动载荷工况，其动刚度偏低，从而验证了前面的假设。另外，动力学{t_{\text{f}}} = 0.05{\;\text{ s}}工况下的优化拓扑在内侧直角处堆积了更多的高弹性模量材料，使得其动刚度大于{t_{\text{f}}} = 0.03{\;\text{ s}}工况下优化拓扑的动刚度。
![]()
图 14 双材料L型结构问题动静态载荷优化结果对比Figure 14. Comparison of dynamic and static load optimization results of bimaterial L-shaped structure![]()
图 15 最优L型结构载荷作用点处的垂向位移和动柔度时间历程Figure 15. Time history of vertical displacement and dynamic flexibility at the loading point of optimal L-shaped structure图16为{t_{\text{f}}} = 0.05{\;\text{ s}}L型结构的多材料拓扑优化结果。多相材料结构的整体拓扑构型基本一致，载荷作用点附近的三角形构型的强化方式为：增加内侧棱边的截面宽度和内部添加加强筋。随着材料相数目的增加，高弹性模量材料分布于固定端垂向平行杆附近，低弹性模量材料分布于载荷作用点附近的三角形构型。但是，由于上述两种强化作用，致使载荷作用点附近的三角形构型与固定端垂向平行杆构型能够达到刚度匹配，从而保证整体结构的动刚度最优。对于6相和10相多材料设计，最软相材料位于载荷作用点附近的三角形构型的内部棱边，致使结构的平均动柔度均高于其他多相材料设计。
![]()
图 16 L型结构的多材料拓扑优化Figure 16. Multi-material topology optimization of L-shaped structure3.4 Michell型悬臂梁
如图17所示，Michell型悬臂梁右侧自由边中点处承受正弦半波载荷f\left( t \right) = {f_0}\sin \left( {\pi t/{t_{\text{f}}}} \right)，载荷幅值{f_0}{\text{ = }}1{\;\text{ kN}}，初始设计域的结构尺寸为：L = 5{\;\text{ m}}，H = 4{\;\text{m}}，h = 0.01{\;\text{m}}，R = 1{\;\text{ m}}。多相材料的基准弹性模量{E_0} = 200{\;\text{ GPa}}，任一材料相的泊松比{\nu _i}{\;\text{ = }}0.3，密度{\rho _i}{\text{ = }}7800{\;\text{ kg/}}{{\text{m}}^3}，瑞利阻尼参数{\alpha _{\text{r}}}{\text{ = }}2.5，{\beta _{\text{r}}}{\text{ = }}4.5 \times {10^{ - 4}}。
图18为不同载荷持续时间工况下双材料Michell型悬臂梁问题的优化设计。与静力学工况的优化设计相比，动力学工况的优化拓扑内部加强肋的具有更大的截面尺寸，载荷作用点附近聚集了更多的高弹性模量材料，因而后者具有更高的动刚度。此外，动力学{t_{\text{f}}} = 0.05{\;\text{ s}}工况下的优化拓扑内部加强肋的宽度尺寸较大，其方位更接近于垂向，有利于抑制动载荷作用下的变形，相对于{t_{\text{f}}} = 0.03{\;\text{ s}}工况下优化拓扑具有更高的动刚度。如图19所示，静力学工况的位移与动柔度响应均高于动载荷工况，其动刚度偏低，而动力学{t_{\text{f}}} = 0.05{\;\text{ s}}工况的上述动响应最低，其动刚度最高，与优化拓扑的分析结果一致。
![]()
图 18 双材料Michell悬臂梁问题动静态载荷优化结果对比Figure 18. Comparison of dynamic and static load optimization results of bimaterial Michell cantilever beam![]()
图 19 最优Michell型悬臂梁结构载荷作用点处的垂向位移和动柔度时间历程Figure 19. Time history of vertical displacement and dynamic flexibility at the loading point of optimal Michell cantilever beam structure图20为Michell型悬臂梁的多材料拓扑优化结果。外环拓扑极为相似，内部强化构型差异显著。高比例的硬材料分布于外部拓扑，而软材料在内部联接于内部拓扑。另外，载荷作用点附近也集中了高弹性模量材料，抵消惯性载荷的动态影响。对于9相和10项多材料设计，最软相或次最软相材料位于载荷作用点附近的三角形构型的内侧棱边，致使结构的平均动柔度均高于其他多相材料设计。
![]()
图 20 Michell型悬臂梁的多材料拓扑优化Figure 20. Multi-material topology optimization of Michell cantilever beam4 结论
本文基于多分辨率-多边形单元建模策略，提出了求解多材料结构动响应问题的拓扑优化方法，研究结论如下：
(1) 通过悬臂梁、简支梁、L型结构和Michell型悬臂梁等典型数值算例，验证了所提多材料结构时域动刚度拓扑优化方法的可行性。
(2) 动态载荷的作用时间对于优化结果具有重要的影响。随着载荷作用时间的减少，结构内惯性力的影响更加显著，致使多材料结构的优化拓扑存在明显的差异。
(3) 将本文所提出的多材料结构动力学优化方法延伸至超材料结构的动力学设计问题，可以实现高分辨率胞元-宏观结构构型的动力学多尺度协同优化设计。
基于提出的动力学多材料拓扑优化框架，作者将在后面工作中考虑动力学载荷作用下多材料之间的界面损伤与解耦，提出多材料结构界面失效问题的动力学拓扑优化方法。
-
![]()
图 4 悬臂梁双材料设计的迭代历史({t_{\text{f}}}{\text{ = }}0.05{\;\text{ s}})
Figure 4. Iteration histories of a bimaterial cantilever beam design ({t_{\text{f}}}{\text{ = }}0.05{\;\text{ s}})
![]()
图 5 双材料悬臂梁问题动静态载荷优化结果对比
Figure 5. Comparison of dynamic and static load optimization results of bimaterial cantilever beam
![]()
图 6 最优悬臂梁载荷作用点处垂向位移和动柔度时间历程
Figure 6. Time history of vertical displacement and dynamic flexibility at the loading point of optimal cantilever structure
![]()
图 8 悬臂梁多材料最优拓扑的动柔度时间历程
Figure 8. Time history of dynamic flexibility for optimized multi-material topology of cantilever beam
![]()
图 10 双材料简支梁问题动静态载荷优化结果对比
Figure 10. Comparison of dynamic and static load optimization results of bimaterial simply supported beam
![]()
图 11 最优简支梁载荷作用处垂向位移和动柔度时间历程
Figure 11. Time history of vertical displacement and dynamic flexibility at the loading point of the optimal simply supported beam structure
![]()
图 12 简支梁的多材料拓扑优化
Figure 12. Multi-material topology optimization of simply supported beams
![]()
图 14 双材料L型结构问题动静态载荷优化结果对比
Figure 14. Comparison of dynamic and static load optimization results of bimaterial L-shaped structure
![]()
图 15 最优L型结构载荷作用点处的垂向位移和动柔度时间历程
Figure 15. Time history of vertical displacement and dynamic flexibility at the loading point of optimal L-shaped structure
![]()
图 16 L型结构的多材料拓扑优化
Figure 16. Multi-material topology optimization of L-shaped structure
![]()
图 18 双材料Michell悬臂梁问题动静态载荷优化结果对比
Figure 18. Comparison of dynamic and static load optimization results of bimaterial Michell cantilever beam
![]()
图 19 最优Michell型悬臂梁结构载荷作用点处的垂向位移和动柔度时间历程
Figure 19. Time history of vertical displacement and dynamic flexibility at the loading point of optimal Michell cantilever beam structure
-
[1] LÓPEZ C, BURGGRAEVE S, LIETAERT P, et al. Model-based, multi-material topology optimization taking into account cost and manufacturability [J]. Structural and Multidisciplinary Optimization, 2020, 62(6): 2951 − 2973. doi: 10.1007/s00158-020-02641-0
[2] 张秋月, 安鲁陵, 岳烜德, 等. 基于遗传算法的飞机复合材料结构装配压紧力大小与布局的优化[J]. 复合材料学报, 2019, 36(6): 1546 − 1557. doi: 10.13801/j.cnki.fhclxb.20180816.001 ZHANG Qiuyue, AN Luling, YUE Xuande, et al. Optimization of size and layout of pressing force for composite airframe structure assembly based on genetic algorithm [J]. Acta Materiae Compositae Sinica, 2019, 36(6): 1546 − 1557. (in Chinese) doi: 10.13801/j.cnki.fhclxb.20180816.001
[3] LI Y, LAI Y P, LU G, et al. Innovative design of long-span steel-concrete composite bridge using multi-material topology optimization [J]. Engineering Structures, 2022, 269: 114838. doi: 10.1016/j.engstruct.2022.114838
[4] 闫晓磊, 陈佳文, 张树忠, 等. 基于骨架提取的拓扑优化最小尺寸精确控制[J]. 工程力学, 2021, 38(5): 239 − 246. doi: 10.6052/j.issn.1000-4750.2020.02.0108 YAN Xiaolei, CHEN Jiawen, ZHANG Shuzhong, et al. Precise control of minimum length scale in topology optimization based on skeleton extraction [J]. Engineering Mechanics, 2021, 38(5): 239 − 246. (in Chinese) doi: 10.6052/j.issn.1000-4750.2020.02.0108
[5] BAI J T, ZUO W J. Multi-material topology optimization of coated structures using level set method [J]. Composite Structures, 2022, 300: 116074. doi: 10.1016/j.compstruct.2022.116074
[6] HAN D, LEE H. Recent advances in multi-material additive manufacturing: Methods and applications [J]. Current Opinion in Chemical Engineering, 2020, 28: 158 − 166. doi: 10.1016/j.coche.2020.03.004
[7] ZUO W J, SAITOU K. Multi-material topology optimization using ordered SIMP interpolation [J]. Structural and Multidisciplinary Optimization, 2017, 55(2): 477 − 491. doi: 10.1007/s00158-016-1513-3
[8] 杜义贤, 李涵钊, 谢黄海, 等. 基于序列插值模型和多重网格方法的多材料柔性机构拓扑优化[J]. 机械工程学报, 2018, 54(13): 47 − 56. doi: 10.3901/JME.2018.13.047 DU Yixian, LI Hanzhao, XIE Huanghai, et al. Topology optimization of multiple materials compliant mechanisms based on sequence interpolation model and multigrid method [J]. Journal of Mechanical Engineering, 2018, 54(13): 47 − 56. (in Chinese) doi: 10.3901/JME.2018.13.047
[9] 闫浩, 吴晓明. 基于ordered-EAMP模型的多材料传热结构拓扑优化[J]. 航空动力学报, 2021, 36(5): 1007 − 1021. doi: 10.13224/j.cnki.jasp.2021.05.012 YAN Hao, WU Xiaoming. Multi-material topology optimization for heat transfer structure based on ordered-EAMP model [J]. Journal of Aerospace Power, 2021, 36(5): 1007 − 1021. (in Chinese) doi: 10.13224/j.cnki.jasp.2021.05.012
[10] 俞燎宏, 荣见华, 唐承铁, 等. 基于可行域调整的多相材料结构拓扑优化设计[J]. 航空学报, 2018, 39(9): 117 − 133. doi: 10.7527/S1000-6893.2018.22023 YU Liaohong, RONG Jianhua, TANG Chengtie, et al. Multi-phase material struclural topology optimization design based on feasible domain adjustment [J]. Acta Aeronautica et Astronautica Sinica, 2018, 39(9): 117 − 133. (in Chinese) doi: 10.7527/S1000-6893.2018.22023
[11] 朱本亮, 张宪民, 李海, 等. 基于节点密度插值的多材料柔顺机构拓扑优化[J]. 机械工程学报, 2021, 57(15): 53 − 61. doi: 10.3901/JME.2021.15.053 ZHU Benliang, ZHANG Xianmin, LI Hai, et al. Topology optimization of multi-material compliant mechanisms using node-density interpolation scheme [J]. Journal of Mechanical Engineering, 2021, 57(15): 53 − 61. (in Chinese) doi: 10.3901/JME.2021.15.053
[12] 赵清海, 张洪信, 华青松, 等. 周期性多材料结构稳态热传导拓扑优化设计[J]. 工程力学, 2019, 36(3): 247 − 256. doi: 10.6052/j.issn.1000-4750.2018.01.0002 ZHAO Qinghai, ZHANG Hongxin, HUA Qingsong, et al. Multi-material topology optimization of steady-state heat conduction structure under periodic constraint [J]. Engineering Mechanics, 2019, 36(3): 247 − 256. (in Chinese) doi: 10.6052/j.issn.1000-4750.2018.01.0002
[13] ZHOU M D, GENG D. Multi-scale and multi-material topology optimization of channel-cooling cellular structures for thermomechanical behaviors [J]. Computer Methods in Applied Mechanics and Engineering, 2021, 383: 113896. doi: 10.1016/j.cma.2021.113896
[14] VANTYGHEM G, BOEL V, STEEMAN M, et al. Multi-material topology optimization involving simultaneous structural and thermal analyses [J]. Structural and Multidisciplinary Optimization, 2019, 59(3): 731 − 743. doi: 10.1007/s00158-018-2095-z
[15] 龙凯, 谷先广, 王选. 基于多相材料的连续体结构动态轻量化设计方法[J]. 航空学报, 2017, 38(10): 134 − 143. LONG Kai, GU Xianguang, WANG Xuan. Lightweight design method for continuum structure under vibration using multiphase materials [J]. Acta Aeronautica et Astronautica Sinica, 2017, 38(10): 134 − 143. (in Chinese)
[16] ZHANG X S, PAULINO G H, RAMOS A S. Multi-material topology optimization with multiple volume constraints: A general approach applied to ground structures with material nonlinearity [J]. Structural and Multidisciplinary Optimization, 2018, 57(1): 161 − 182. doi: 10.1007/s00158-017-1768-3
[17] YANG X T, LI M. Discrete multi-material topology optimization under total mass constraint [J]. Computer-Aided Design, 2018, 102: 182 − 192. doi: 10.1016/j.cad.2018.04.023
[18] HUANG X D, LI W B. A new multi-material topology optimization algorithm and selection of candidate materials [J]. Computer Methods in Applied Mechanics and Engineering, 2021, 386: 114114. doi: 10.1016/j.cma.2021.114114
[19] LIU J K, MA Y S. A new multi-material level set topology optimization method with the length scale control capability [J]. Computer Methods in Applied Mechanics and Engineering, 2018, 329: 444 − 463. doi: 10.1016/j.cma.2017.10.011
[20] 王选, 胡平, 龙凯. 考虑嵌入移动孔洞的多相材料布局优化[J]. 力学学报, 2019, 51(3): 852 − 862. doi: 10.6052/0459-1879-18-327 WANG Xuan, HU Ping, LONG Kai. Multiphase material layout optimization considering embedding movable holes [J]. Chinese Journal of Theoretical and Applied Mechanics, 2019, 51(3): 852 − 862. (in Chinese) doi: 10.6052/0459-1879-18-327
[21] GANGL P. A multi-material topology optimization algorithm based on the topological derivative [J]. Computer Methods in Applied Mechanics and Engineering, 2020, 366: 113090. doi: 10.1016/j.cma.2020.113090
[22] ZHANG W S, SONG J F, ZHOU J H, et al. Topology optimization with multiple materials via moving morphable component (MMC) method [J]. International Journal for Numerical Methods in Engineering, 2018, 113(11): 1653 − 1675. doi: 10.1002/nme.5714
[23] WANG X, LONG K, MENG Z, et al. Explicit multi-material topology optimization embedded with variable-size movable holes using moving morphable bars [J]. Engineering Optimization, 2021, 53(7): 1212 − 1229. doi: 10.1080/0305215X.2020.1779710
[24] SUKUMAR N. Construction of polygonal interpolants: A maximum entropy approach [J]. International Journal for Numerical Methods in Engineering, 2004, 61(12): 2159 − 2181. doi: 10.1002/nme.1193
[25] TALISCHI C, PAULINO G H, PEREIRA A, et al. PolyTop: A Matlab implementation of a general topology optimization framework using unstructured polygonal finite element meshes [J]. Structural and Multidisciplinary Optimization, 2012, 45(3): 329 − 357. doi: 10.1007/s00158-011-0696-x
[26] SANDERS E D, PEREIRA A, AGUILÓ M A, et al. PolyMat: An efficient Matlab code for multi-material topology optimization [J]. Structural and Multidisciplinary Optimization, 2018, 58(6): 2727 − 2759. doi: 10.1007/s00158-018-2094-0
[27] GIRALDO-LONDOÑO O, PAULINO G H. PolyDyna: A Matlab implementation for topology optimization of structures subjected to dynamic loads [J]. Structural and Multidisciplinary Optimization, 2021, 64(2): 957 − 990. doi: 10.1007/s00158-021-02859-6
[28] GIRALDO-LONDOÑO O, AGUILÓ M A, PAULINO G H. Local stress constraints in topology optimization of structures subjected to arbitrary dynamic loads: A stress aggregation-free approach [J]. Structural and Multidisciplinary Optimization, 2021, 64(6): 3287 − 3309. doi: 10.1007/s00158-021-02954-8
[29] NGUYEN-XUAN H. A polytree-based adaptive polygonal finite element method for topology optimization [J]. International Journal for Numerical Methods in Engineering, 2017, 110(10): 972 − 1000. doi: 10.1002/nme.5448
[30] CHAU Khai N, CHAU Khanh N, TUAN Ngo, et al. A polytree-based adaptive polygonal finite element method for multi-material topology optimization [J]. Computer Methods in Applied Mechanics and Engineering, 2018, 332: 712 − 739. doi: 10.1016/j.cma.2017.07.035
[31] NGUYEN T H, PAULINO G H, SONG J, et al. A computational paradigm for multiresolution topology optimization (MTOP) [J]. Structural and Multidisciplinary Optimization, 2010, 41(4): 525 − 539. doi: 10.1007/s00158-009-0443-8
[32] FILIPOV E T, CHUN J, PAULINO G H, et al. Polygonal multiresolution topology optimization (PolyMTOP) for structural dynamics [J]. Structural and Multidisciplinary Optimization, 2016, 53(4): 673 − 694. doi: 10.1007/s00158-015-1309-x
[33] SHOBEIRI V. Bidirectional evolutionary structural optimization for nonlinear structures under dynamic loads [J]. International Journal for Numerical Methods in Engineering, 2020, 121(5): 888 − 903. doi: 10.1002/nme.6249
[34] ZHAO J P, WANG C J. Topology optimization for minimizing the maximum dynamic response in the time domain using aggregation functional method [J]. Computers & Structures, 2017, 190: 41 − 60.
[35] ZHAO J P, WANG C J. Dynamic response topology optimization in the time domain using model reduction method [J]. Structural and Multidisciplinary Optimization, 2016, 53(1): 101 − 114. doi: 10.1007/s00158-015-1328-7
[36] WANG F W, LAZAROV B S, SIGMUND O. On projection methods, convergence and robust formulations in topology optimization [J]. Structural and Multidisciplinary Optimization, 2011, 43(6): 767 − 784. doi: 10.1007/s00158-010-0602-y
[37] STOLPE M, SVANBERG K. An alternative interpolation scheme for minimum compliance topology optimization [J]. Structural and Multidisciplinary Optimization, 2001, 22(2): 116 − 124. doi: 10.1007/s001580100129
[38] STEGMANN J, LUND E. Discrete material optimization of general composite shell structures [J]. International Journal for Numerical Methods in Engineering, 2005, 62(14): 2009 − 2027. doi: 10.1002/nme.1259
[39] HILBER H M, HUGHES T J R, TAYLOR R L. Improved numerical dissipation for time integration algorithms in structural dynamics [J]. Earthquake Engineering & Structural Dynamics, 1977, 5(3): 283 − 292.
[40] GIRALDO-LONDOÑO O, PAULINO G H. Fractional topology optimization of periodic multi-material viscoelastic microstructures with tailored energy dissipation [J]. Computer Methods in Applied Mechanics and Engineering, 2020, 372: 113307. doi: 10.1016/j.cma.2020.113307
[41] GIRALDO-LONDOÑO O, MIRABELLA L, DALLORO L, et al. Multi-material thermomechanical topology optimization with applications to additive manufacturing: Design of main composite part and its support structure [J]. Computer Methods in Applied Mechanics and Engineering, 2020, 363: 112812. doi: 10.1016/j.cma.2019.112812
-
期刊类型引用(1)
1. 邹胤康，李少华，邱文科，夏凉. 面内周期性结构多分辨率拓扑优化设计. 固体力学学报. 2024(04): 533-546 . 
百度学术
其他类型引用(1)
下载:
计量
- 文章访问数: 299
- HTML全文浏览量: 93
- PDF下载量: 75
- 被引次数: 2
