基于改进杜阿梅尔积分的大规模长期瞬态拓扑优化方法

《Engineering Structures》：Large-scale long-term transient topology optimization method based on improved Duhamel integral

　　摘要
尽管瞬态拓扑优化在动态工程应用中至关重要，但由于传统时间积分方法如Newmark-β法的迭代特性，大规模长期问题的计算成本仍然过高。为提高设计效率与精度，我们提出了一种基于Duhamel积分的高效优化框架，用于解决长时间加载下具有高自由度的大规模问题。通过模态简化将大规模系统转化为解耦的单自由度问题，我们得到了高效的递归瞬态响应表达式，并利用模态加速度法确保在高阶模态截断情况下仍能保持精度。与传统数值积分不同，这种解析方法即使在较大的时间步长下也能保持无条件稳定性，且可通过显式公式直接进行灵敏度分析。包含数百万自由度及长时间跨度的数值算例表明，所提方法显著降低了计算成本，使得大规模瞬态拓扑优化得以实际应用。

引言
拓扑优化是一种主流的结构设计方法，旨在在预定的设计域内寻求材料的最佳空间分布[1]。该方法已广泛应用于航空航天[2]、汽车[3]和海洋工程[4]领域。然而，在实际工程中，结构往往要承受动态和瞬态载荷[5]，因此将时变结构响应纳入拓扑优化框架已成为重要的研究方向。

关于动态载荷作用下的拓扑优化研究，主要分为谐波激励下的频域分析[6]以及涉及瞬态载荷的时域分析[7]。为避免频率耦合引起的共振，基于谐波的优化方法通常侧重于最大化基频[8]或最小化指定频率范围内的位移幅值[9]。虽然模态位移法[9]、模态加速度法[9]和全方法[9]等已被广泛用于处理谐波响应，但近期研究也致力于解决大规模约束问题。例如，Long等人[10]在全方法框架中引入p-范数聚合函数，以处理复杂的局部约束。

然而，当结构承受非周期性、冲击或随机载荷时，时域瞬态优化就显得不可或缺。这类方法通常采用时间积分方案来求解结构动力学问题，同时结合伴随变量法进行灵敏度分析[11]。Kang等人[12]记录了瞬态灵敏度分析的历史进展，而Jensen等人[13]则提出了“离散化后再求导”的伴随变量法范式，以解决传统方法中存在的灵敏度不一致问题。该方法已扩展到具有非线性阻尼器的滞后系统等复杂系统[14]。最近，Giraldo-Londo?o等人[15]通过PolyDyna这一开源MATLAB框架普及了这些成果，该框架将“离散化后再求导”的伴随方法与基于ZPR的设计更新相结合。在动态载荷作用下的频域和时域拓扑优化方面都取得了显著进展，包括谐波响应方法的改进，以及基于伴随灵敏度计算和开源框架支持的稳健瞬态分析范式的建立。不过，在实际工程应用中，面对大规模自由度问题时，瞬态拓扑优化仍面临诸多挑战。

通过时间积分进行大规模结构优化往往会导致极高的计算复杂度[16]。为缓解这一问题，Slone等人[17]将模态简化技术应用于瞬态分析。Zhao等人[18]的进一步研究表明，虽然模态位移法在时域中的精度不足，但模态加速度法在处理大规模问题时表现优异。值得注意的是，在较大时间步长下，基于模态加速度法的优化效果优于直接积分法。在此基础上，Kristiansen等人[19]开发了一种基于模态简化的开源C++/MPI框架，用于大规模瞬态优化。为进一步提升效率，Zhu等人[16]提出了准静态响应增强型模态位移法（QMDM），并结合逐次迭代分析与设计法[20]来降低计算成本。此外，等效静载荷法[21][22][23]也被提出作为一种加速手段，以应对灵敏度分析中的难题。Jang等人[24]通过将动态载荷转换为等效静载荷，发展了用于时域动态拓扑优化的等效静载荷法。Su等人[25]则将缩放边界有限元法与双向进化结构优化算法相结合，基于高阶隐式时间积分求解动态方程，再利用等效静载荷法进行拓扑优化。尽管有这些进展，但对于数百万自由度且在长时间跨度内的运动方程（如Newmark-β方程）的求解依然是一项艰巨任务。

与传统时间积分方案相比，Duhamel积分方法在瞬态载荷分析中具有更高的计算效率和数值精度[26][27]。其显式的解析性质确保了在长期动态激励下的稳健性能，有效避免了逐步积分所带来的累积误差和高成本问题[28]。该框架最初由Jean-Marie Constant Duhamel为热传导问题提出[29]，随后在结构动力学领域得到广泛应用与完善[30]，并成功应用于汽车[31]、航空航天[32]和海洋工程[33]行业，用于解决结构振动[34]、疲劳寿命[35]和应力传播[36]等关键问题。近期的理论进展进一步拓展了其应用范围；例如，Zhang等人[37]基于Duhamel积分开发了一种快速的递归半解析方法，用于处理复杂的边界反射问题，而Tomczyk等人[38]则利用Duhamel-Neumann热弹性理论对圆柱壳的动态响应进行了建模。尽管Duhamel方法在动态分析中已展现出明显优势，但将其整合到瞬态问题的拓扑优化框架中仍是一个未被探索的领域。这一研究空白为提升瞬态结构设计的效率与精度提供了新的机遇。

本研究提出了一种改进的Duhamel积分方法（IDM），以解决瞬态动态分析中的难题。通过在模态加速度框架中将Duhamel积分改写为时间步长递归形式，所提方法即使在相对较大的时间步长下也能保持良好的数值稳定性。此外，位移灵敏度也可通过解析方式获得，从而大幅节省计算资源。与传统时间积分方法需要依次进行正向动态响应计算和反向伴随变量求解不同，IDM仅需一次迭代循环即可同时获得结构响应和位移灵敏度。由于采用了直接灵敏度分析，该框架显著提升了大规模长时间跨度瞬态拓扑优化的计算效率。

本文其余结构如下：第2节介绍了瞬态动力学的拓扑优化模型及模态简化理论；第3节详细阐述了基于改进Duhamel方法的优化框架；第4节通过一系列数值算例验证了该算法的有效性；第5节对算例结果进行了分析讨论；最后为结论部分。

片段摘录
基于传统模态简化方法的瞬态动态拓扑优化
此处我们介绍瞬态拓扑优化模型。以动态柔度作为目标函数的拓扑优化能够在瞬态载荷作用下实现结构的最大刚度设计。为提高计算效率，采用了模态简化方法，相关细节将在后续小节中介绍。

长期与大规模问题的高效近似分析
为延长时间步长并确保数值稳定性，本节提出了一种改进的Duhamel积分方法，该方法能够在更长的时间步长下实现高精度计算，同时基于该积分方法推导出一种直接的灵敏度分析方法，以便进一步进行拓扑优化。

数值算例
在本节中，我们通过多个数值算例展示了改进的Duhamel积分方法在拓扑优化中的应用效率。所有计算程序均采用MATLAB实现。在所有算例中，弹性模量E=2.1×10^5 MPa，泊松比ν=0.3，密度ρ=7.85×10^-9 t/mm3，瑞利阻尼系数分别设为α=0.01和β=0.01。设计域采用六面体单元进行离散化，RAMP插值参数取为qE=5。

讨论
在瞬态拓扑优化中，不同的动态载荷会导致不同的最优配置。对于具有大量自由度且时间历程较长的瞬态拓扑优化问题，主要的计算成本来自两个方面：动态微分方程的求解以及目标函数的灵敏度分析。传统的方程求解和灵敏度计算数值方法会大大增加整体计算量。

结论
本文提出了一种用于瞬态拓扑优化问题的改进Duhamel积分方法。在具有大规模自由度和长期优化问题的情况下，求解动态微分方程会耗费大量时间。为解决这一难题，我们首先提出了一种改进的Duhamel积分方法，该方法通过模态简化避免了大规模自由度带来的问题，还将传统的Duhamel积分改为每个时间步长下的递归形式。