通过深度神经网络实现多尺度结构和材料的双尺度拓扑优化，以达到局部和全局的抗屈曲性能

《COMPUTER METHODS IN APPLIED MECHANICS AND ENGINEERING》：Two-scale topology optimization of multiscale structures and materials via deep neural networks for local and global buckling

　　索布汉·霍纳瓦尔（Sobhan Honarvar）、诺兰·布莱克（Nolan Black）和艾哈迈德·R·纳贾菲（Ahmad R. Najafi）
德雷塞尔大学机械工程与力学系，美国宾夕法尼亚州费城栗子街3141号，邮编19104

**摘要**
多尺度拓扑优化能够创造出具有定制微观结构的轻量化结构，但同时实现全局（宏观尺度）和局部（微观尺度）的屈曲限制在计算上仍然具有挑战性，尤其是对于参数化复杂的单元格而言。本文提出了一种可微分的、基于深度神经网络（DNN）的双尺度屈曲拓扑优化框架，该框架将宏观尺度布局优化与微观尺度稳定性保护结合在一个基于梯度的迭代过程中。通过在有限元均质化数据上训练深度神经网络代理模型，可以预测参数化单元格的均质化本构张量及其敏感性（通过反向传播实现）。这种方法消除了对插值表的需求，并避免了在高维微观结构设计空间中难以应用的有限差分梯度法。该代理模型被嵌入到一个双尺度公式中，旨在最大化宏观尺度屈曲载荷因子，同时确保在最坏情况下微观尺度屈曲载荷因子得到满足，以防止在任意应力状态下的局部失稳。

数值研究表明：（1）通过双尺度体积分割和权重分配，在宏观与微观稳定性之间实现了稳健的权衡；（2）对于设计维度较高的微观结构，传统的每次迭代有限元均质化方法无法有效优化；（3）新的微观结构观察结果表明：基于平滑嵌件的（椭圆形）单元格在激活微观尺度屈曲保护时，可以达到与晶格基单元格相当的全局屈曲性能，同时表现出更小的性能下降。

**1. 引言**
拓扑优化是一种通过高效分配材料来产生最优结构的计算方法，通常优于传统的基于直觉的设计方法[1]。近年来，增材制造技术的进步使得填充结构的制造成为可能，进一步推动了拓扑优化在多尺度设计中的应用[2][3]。这一发展意味着，曾经仅限于理论研究的具有定制微观结构的结构现在可以在航空航天[4]、汽车[5]、建筑[6]和医疗保健[7]等多个行业的实际工程中得到应用。晶格结构相比传统的固态或均匀填充结构具有显著优势，例如更高的强度重量比、更好的能量吸收能力以及更强的整体抗屈曲性能，从而扩展了创新和高性能设计的潜力[8][9]。然而，晶格结构的复杂细部特征可能容易发生局部屈曲现象。

为了缓解晶格结构中的局部屈曲问题，人们对双尺度优化产生了越来越多的研究兴趣。这种方法能够在单一计算环境中整合宏观和微观尺度信息，无需完整详细的模型[1]，从而可以在两个尺度上优化结构的性能[1]。为了处理这些尺度，通常采用渐近均质化方法来有效桥接全局和局部结构特征之间的差距[10]。渐近均质化的一个常见应用是数值均质化，它通过分析微观尺度单元格来近似复合材料的有效性能[11]。尽管数值均质化在多尺度设计优化中得到广泛应用，但由于需要严格的尺度分离，并且经常依赖于简化的微观结构表示（而非真实的复杂几何形状）[1][13][14]，因此面临重大计算挑战。

许多研究提出了不同的策略来应对这些挑战。其中一种方法是对微观尺度单元格进行参数化表示，以减少大规模设计空间带来的计算挑战[14][15][16]。在最近的一项研究中，翟等人[17]引入了一种使用高度场的参数化微观结构的可微表示方法，该方法能够实现几何形状和材料属性在一系列微观结构中的平滑过渡，从而提高了拓扑优化的效率和连续性。虽然参数化微观结构使优化和分析变得更加可行，但如果选定的参数未能完全捕捉到微观结构行为的全部范围，也可能导致性能不佳。

在参数化微观结构中，需要从均质化设计中重构物理上详细的几何形状，以确保结构性能和可制造性[3][18]。作为后处理方法，去均质化可以将经过优化的均质化设计（其中材料属性被平均化）转换为详细的、可制造的结构，从而实现多尺度布局的实际应用[14][19]。对于简单周期性微观结构，去均质化最为直接，因为其潜在模式是规则且易于理解的[3][18]。然而，最近的进展已将去均质化技术扩展到更复杂或分级的微观结构中，将空间变化和各向异性纳入逆映射过程[20]。这种扩展增加了计算复杂性，通常需要高级的后处理或数据驱动技术[19]。去均质化过程中可能会出现问题，例如由于微观结构粗糙导致体积约束被违反，可能导致连接不连续的构件，从而易于发生局部屈曲或导致性能劣于均质化设计[3]。为了解决这个问题，可以添加矫正壳体来连接边缘处不完整的支柱，从而加强结构。这一过程计算量较大，需要精细的网格投影和跨多种单元格尺寸及位置的广泛验证。

为了减轻多尺度优化中的计算负担，人们采用了代理模型来近似有效的材料属性。这些模型能够高效探索大型设计空间，应对高参数化和结构复杂性带来的挑战。多项式回归和多维代理模型等方法建立了微观尺度特征与宏观尺度行为之间的定量关联[21][22]。在此基础上，最近的发展转向了架构材料，其性能更多由几何形状而非成分决定。基于超材料的策略利用几何设计来定制所需的微观尺度特性，通常通过数据驱动的建模技术进行增强[23][24]。为了应对详细微观几何形状带来的挑战，机器学习（ML），特别是神经网络（NN），已成为数值均质化的有效替代方案[25][26]。多种神经架构被用来将晶格配置映射到均质化弹性张量[27]，捕获多样的微观结构类型和几何形状[28][29]，并预测复杂结构（如Voronoi基材料）的行为[30]。这些模型促进了高效的多尺度设计，并且当与数值均质化或有限元方法结合使用时可以显著降低计算成本[31]。基于这些进展，最近的研究展示了深度神经网络在并行多尺度优化中的参数化均质化可行性[32]，并开发了结合基于NN的均质化与应力评估代理的压力约束多尺度优化框架[33]。最近，布莱克和纳贾菲[34]证明了使用针对雅可比矩阵和应力精度的相对L2 Sobolev损失进行训练的NN代理模型，可以嵌入到并行多尺度优化中，训练案例数量约为10^4个，加速效果接近两个数量级。

除了之前提到的双尺度优化问题的局限性外，屈曲问题的非线性特征进一步增加了多尺度屈曲优化的计算成本，相比之下，合规性等目标问题的计算成本较低。因此，最近关于多尺度结构拓扑优化的研究主要集中在合规性最小化上，而没有考虑屈曲问题。因此，大多数现有的屈曲拓扑优化研究仅限于单尺度公式。费拉里等人[35]探索了线性屈曲分析作为解决大规模拓扑优化问题的方法。在此基础上，引入了一种结合了屈曲阻力和局部延性失效约束的复合公式[36]。然而，传统的单尺度方法往往忽略局部失稳性，导致设计在宏观尺度上表现出很强的抗屈曲能力，但在微观尺度上仍然容易发生屈曲。

最近的多尺度拓扑优化进展试图通过整合插值和去均质化等方法来克服双尺度优化中的局限性，从而同时考虑局部和全局屈曲现象[3][18]。克里斯滕森等人[3]开发了一种方法，该方法结合了各向同性多孔微观结构的优化（针对屈曲稳定性）和传统的刚度优化。该方法将基于Bloch–Floquet波理论的局部密度依赖性屈曲应力约束与全局屈曲准则相结合，从而更准确地预测不同长度尺度下的结构稳定性。此外，优化过程允许去均质化，从而产生的设计更适用于增材制造等制造工艺。

在这些发展的基础上，休布纳等人[18]提出了一种针对分级晶格结构的双尺度优化方法，考虑了微观和宏观尺度的屈曲行为。他们的方法利用渐近均质化来提升周期性微观结构的机械性能和屈曲载荷因子（BLFs）。通过预计算均质化属性并应用插值模型，该方法实现了高效的优化，避免了昂贵的全尺度模拟需求。引入了最坏情况下的均质化BLF模型以缓解局部屈曲问题，提供了更加稳健且计算效率高的解决方案。然而，许多这些策略仍然严重依赖插值和去均质化，这可能会引入近似误差或将分析限制在单个尺度上。虽然它们在稳定性预测方面有了显著改进，但在完全捕捉宏观和微观尺度屈曲之间的相互作用方面仍存在保守的差距。

在这项工作中，我们提出了一种由DNN驱动的双尺度屈曲拓扑优化框架，旨在实现鲁棒性和可扩展性。我们强调，机器学习并非直接作为内部循环求解器嵌入；相反，我们采用了一种基于梯度的、由DNN加速的方法，其中代理模型在离线状态下根据有限元均质化数据进行训练，以低边际成本预测均质化本构张量及其分析敏感性。这消除了基于插值的本构管道，并避免了有限差分微观尺度梯度，即使在单元格由许多几何参数描述时也能实现高效的基于梯度的优化。然后将代理模型与并行双尺度稳定性公式相结合：（i）针对结构布局优化宏观尺度屈曲载荷因子；（ii）通过应力空间下的最坏情况均质化屈曲载荷因子来确保微观尺度稳定性，防止局部失稳，从而避免使全局“最优”设计失效。除了计算加速外，该流程还通过降低本构评估的成本，使得更丰富的微观结构家族和更高维度的单元格参数化能够被优化和比较，揭示了在传统FE-in-the-loop多尺度屈曲优化中难以观察到的稳定性趋势。

本文的其余部分安排如下：第2节介绍了用于近似均质化本构属性及其敏感性的DNN代理模型。第3节阐述了多尺度屈曲拓扑优化问题，包括宏观和微观尺度的屈曲分析。第4节提出了所提出的双尺度优化框架，并详细描述了问题公式、敏感性分析和过滤策略。第5节提供了数值结果，以评估DNN代理模型的有效性，并通过各种案例研究证明了所提框架的有效性。最后，第6节对研究结果进行了总结，并指出了未来研究的方向。二维多尺度表示：(a) 在载荷和位移边界条件下的宏观尺度域 (ΩM) (tponΓMt 和uponΓMu)；(b) 在X处定义的具有周期性变化的微观结构的微观尺度几何形状 (Ωm)；(c) 一种晶格微观结构的示例。一旦计算出位移场，就可以使用方程(4)对微观结构进行均质化处理，该方程是适用于有限元分析(FEA)的方程(1)的矩阵表示形式：(4)CH=CijH=1|Y|∑e=1m∫Yχe0(i)?χe(i)Tkeχe0(j)?χe(j)dY，其中索引i和j是指均质化本构张量分量的张量索引，χe0包含由单位应变载荷产生的位移场，而χe包含实际微观结构问题的位移解。符号m和ke分别代表元素数量和元素刚度矩阵。括号中的索引表示每个单位应变情况对应的列。

2.2. 用于均质化的深度神经网络模型
典型的神经网络（NN）由一系列矩阵运算组成。对于输入向量x=h0，每一层Ti根据特定的权重矩阵Wi和偏置向量bi使用公式(5)计算一个中间潜在表示hi：Ti+1=Wihi+bi，其中Ti+1通常代表第i+1层的线性变换。可选地，输出可以通过非线性激活函数A（如Sigmoid、tanh或ReLU）进行变换，从而得到公式(6)Ai(Ti)=Ai。另一方面，h(i+1)是在应用非线性（激活函数）后的结果：h(i+1)=A(i+1)T(i+1)，其中A(i+1)(?)表示第i+1层的激活函数。当多层堆叠时，网络成为DNN，其输出预测y?由L层中激活函数和仿射变换的嵌套应用组成：DNN(x)=y?=AL°AL?1°?°A1(T1(x))。DNN的可训练参数，即权重W和偏置b，通过最小化目标函数L来优化。在数值均质化中，微观结构设计参数作为输入x，而均质化弹性张量CH的分量构成目标输出向量y?。常用的损失函数是均方误差（MSE）：L(W,b,y)=1/Ny(y?y?)?(y?y?)，其中Ny是输出变量的数量，y表示真实值。DNN通过使用预计算的训练数据集来最小化预测输出和已知输出之间的差异进行训练。对于均质化任务，该数据集通常包括微观尺度几何形状及其对应的均质化本构张量对。为了优化网络参数，使用反向传播算法计算梯度。链式法则使得相对于输入的输出导数能够高效计算：?y??x=?y??AL??AL?TL??TL?AL?1??T1?x。

在当前研究中，输入包括晶格结构的几何参数（具体而言是支撑杆的厚度），而输出是均质化刚度张量的分量（图2）。

下载：下载高分辨率图像（340KB）
下载：下载全尺寸图像
图2. 用于计算均质化本构张量的实现的DNN示意图。DNN的输入是晶格单元的几何参数，表示为输入向量x。隐藏层hi包括优化后的变量W和b，以获得输出y?。

3. 多尺度屈曲
在屈曲背景下，多尺度优化导致了多尺度屈曲优化的发展，其中在设计过程中同时考虑了全局（宏观尺度）和局部（微观尺度）的屈曲现象。如图3所示，宏观尺度被离散化为有限元素（图3(b1)）。然后，宏观尺度中的每个元素（图3(b2)）映射到一个代表性的微观结构（图3(b3)）。通过渐近均质化获得每个微观结构的均质化有效弹性张量。接下来解决微观尺度问题，并使用得到的有效弹性张量迭代更新宏观尺度设计（图3(b1)）。这种双尺度策略能够在不显式解决所有微观结构细节的情况下实现高效优化，同时保持对两种类型屈曲的鲁棒性。弹性、均质化和屈曲分析的潜在理论框架和控制方程将在以下小节中介绍。

下载：下载高分辨率图像（527KB）
下载：下载全尺寸图像
图3. 多尺度屈曲设计框架：(a) 完整模型；(b) 使用双尺度信息更新的多尺度设计；(b1) 具有理化元素的宏观尺度域；(b2) 带有积分点的宏观尺度元素；(b3) 微观尺度级别的屈曲分析。

3.1. 宏观尺度屈曲
假设完美结构在制造过程中没有任何几何缺陷且变形很小，可以使用线性屈曲（LB）分析。在LB分析中，考虑线性弹性设置和线性分叉屈曲条件。在这些条件下，屈曲前的位移、应力和应变是施加载荷和载荷因子的线性函数，这表明了稳定性。为了解决线性屈曲问题，首先应该求解在给定边界条件下的结构位移的线性弹性状态方程。考虑一个体Ω?R2及其有限元离散化ΩM=?e=1MΩe，其中Ωe表示第e个元素的域，M是元素的总数。线性弹性在离散化域ΩM上的控制方程由Zienkiewicz等人[39]给出：Ku=F，其中K是全局刚度矩阵，u是屈曲前的节点位移向量，F是规定的外部载荷向量。方程(11)中的弹性刚度矩阵可以计算为：K=?e=1M∑k=1Nipce,kBe,kTCHBe,k，其中ce,k表示元素e的雅可比行列式和第k个积分点的积分权重。此外，积分点的数量用Nip表示。Be,k是在第k个积分点评估的元素e的应变-位移矩阵，它包含有限元形状函数的导数Ni，对于四节点元素表示为：B=?xN10?xN20?xN30?xN40?yN10?yN20?yN30?yN4?yN1?xN1?yN2?xN2?yN3?xN3?yN4?xN4。

在解决了方程(11)中的屈曲前问题后，可以使用均质化弹性张量和以下方程[40]求解宏观尺度屈曲问题：(K?ΛG(u))?=0，该方程求解特征值-特征向量对Λ?和φ?，其中?≤M。这些特征值称为载荷因子，决定了结构的屈曲响应。最小的正特征值是临界载荷因子，通常称为BLF，是稳定性分析中的关键参数。如果BLF小于或等于1，则预测会发生屈曲。相反，如果BLF大于1，则认为结构在施加的参考载荷下是稳定的。方程(14)中的G是一个初始应力刚度矩阵，也称为几何刚度矩阵，与问题的非线性有关，可以计算为：G(u)=?e=1M∑k=1Nipce,kB?e,kTσe(ue)B?e,k，其中B?e,k和σe分别是在积分点xk评估的形状函数矩阵的导数和元素e内的平均应力。元素e内的平均应力表示为：σe(ue(α))=1/Nip∑k=1NipCHBe,kue，其中ue是元素e内的节点位移，B?e,k是形状函数的梯度矩阵：B?=?xN10?xN20?xN30?xN40?yN10?yN20?yN30?yN4?xN10?xN20?xN30?yN40?yN4。

在求解特征值和特征向量时，方程(14)的解通常按照特征值（即载荷因子）的绝对值升序排列。特征向量被归一化以确保比例一致性，以便直接比较屈曲模式形状并简化灵敏度和稳定性分析。因此，特征向量可以归一化：?TG?=1。

3.2. 微观尺度屈曲
微观尺度上的屈曲平衡方程被构建为一个特征值问题，类似于其宏观对应物：[18] (Kμ?λμGμ)φμ=0。与宏观情况类似，最小的正特征值λμ称为BLF。这里，Kμ和φμ分别是微观尺度的全局刚度矩阵和微观尺度特征向量。此外，Gμ表示微观初始应力刚度矩阵：Gμ=?e=1m∑k=1Nipc?e,kB?e,i,kTσeijB?e,j,k，其中m表示微观尺度上的元素数量。微观应力张量σe定义为：?σe=CI?BeXe??。这里的σe捕获了宏观应变???在单元格Y内的分布。C是基材料的弹性张量，I是单位矩阵。注意，方程(21)已经考虑了应力放大的效应[41]。矩阵Xe=[χe1χe2χe3]包含微观尺度上平衡方程的解：Kmχj=fj，j=1,2,3，以及(fj=?e=1m∑k=1Nipc?e,kBe,kTC??j，其中??j是单位应变场，???ij=δij。宏观应力可以使用宏观本构方程[24] ?σ?=CH??获得。

3.3. Bloch–Floquet理论
在解决宏观尺度屈曲（图4(a)）之后，应将信息从宏观尺度转移到微观尺度以解决微观尺度上的屈曲问题（图4(b)）。然而，在一个原始单元格（图4(d)）中，只能检测到波长适合在该单元格内的屈曲模式。因此，为了捕获更长波长的微观屈曲，微观尺度必须包含多个单元格（图4(c)）。具体来说，微观结构可以被视为Yn周期性单元格，其中n2是2D中的单元格数量，并且注意在n×n单元格中出现的所有模式也出现在任何m′n×m′n单元格中（图4(c)）。因此，增加单元格的数量可以揭示更长波长的不稳定性。为了保守起见，在允许的单元格重复集合上取最小的（即临界）BLF：λμ(σ?)=minn∈Nλμn(σ?)，实际上，这是对于选定的上限Q在n≤Q范围内评估的，其中λμn表示单元格集合中的微观尺度屈曲载荷因子。对于大量的单元格重复n，Bloch–Floquet公式非常高效，因为它解决的是单个单元格并且使用相移边界条件，而不是n×n超胞，同时捕获相同的模式[42]。Bloch–Floquet公式将在下文中讨论。

下载：下载高分辨率图像（367KB）
下载：下载全尺寸图像
图4. (a) 具有边界划分ΓMu和ΓMt的宏观尺度域ΩM；(b) 微观尺度域Ωm；(1) 使用n×n单元格的单元格集合：揭示波长跨越多个单元格的模式；(c) 这种集合的示意图；(2) Bloch–Floquet：仅分析具有相移边界条件的原始单元格，并在布里渊区内扫描k；(d) Bloch问题中使用的原始单元格；(e) 方格晶格的布里渊区；(f) 布里渊区中的优化点。

考虑单元格集合类似于在布里渊区上采样的Bloch–Floquet分析，作为布里渊区的特殊离散化（图4(e)）。通过足够细的布里渊区离散化，Bloch–Floquet自动涵盖了“单元格数量”效应，而无需显式组装大型超胞。二维Bloch-Floquet公式[42]规定：φμ(y+Rj)=φμ(y)eik?Rj,j=1,2，其中虚数单位i，波矢k=(k1,k2)T，晶格向量Rj定义了单元格的几何形状，y是晶格单元格内的位置。方程(26)为在由不同数量的单元格组成的集合上施加边界条件提供了框架。由于这样的单元格集合展示了周期性边界条件，因此有：φμ(njRj)=φμ(0)。这个要求导致了一个可公度性条件：沿Rj方向恰好重复nj个单元格的模式满足集合施加的周期性。因此，控制单元格数量的方程自然由此条件得出：(eik?Rj)nj=1?kTRj=2πmjnj,mj∈Z。在连续术语中，每个波长的单元格数量是：Njcont=2π|kTRj|，其中Njcont是在单元格集合中重复的单元格数量。对于图4(c)中所示的正方形单元格配置，晶格向量表示为：R1=(d,0)T,R2=(0,d)T，其中d是正方形单元格的大小。这个周期性约束确保通过解决跨越布里渊区的波矢的稳定性问题来考虑所有潜在的分叉波长（图4(e)）。产生最小临界载荷的波矢决定了不稳定性的特征波长（即临界屈曲模式）。

解决全域屈曲特征值问题是计算密集型的。尽管数值证据表明，仅采样不可约布里渊区（IBZ；图4(e)）的边界就足以识别临界屈曲模式，但这尚未得到普遍证明。因此，在本研究中，我们评估了整个IBZ上规定的一组波矢，而不是将分析限制在其边界上。在适用的情况下，利用几何对称性和宏观加载对称性在优化过程中减少k点集，从而加速计算。对于等边三角晶格（图4(d2)），图4(f)中的绿色和红色点都被包括在内。对于空心三角晶格（图4(d1)），由于对称性，仅使用绿色点。

在微观尺度平衡和不稳定性方程中实施周期性边界条件时，单元格u的边界节点处的位移如图5所示进行索引。一组单元格的周期性约束是根据Floquet–Bloch关系在方程(26)中强制执行的，从而得到表达式(31)：ur=eikTR1ul, ut=eikTR2ub, utl=eikTR2ubl, utr=eikT(R1+R2)ubl，其中ul, ub, ur, ut, ubl, ubr, utl,utr是图5所示单元格边缘和角落处节点的位移。需要注意的是，在Bloch–Floquet设置中，k=(0,0)T对应于零相位移动，因此平衡问题是具有Y周期性边界条件的。为了明确n×n集合的周期性与方程(31)之间的联系，我们用四个段落Γl,Γr,Γb,Γt（分别代表左、右、下、上）来表示单元格的边界，如图5所示。我们将这些段落上的节点位移自由度（DOFs）收集到向量ul, ur, ub, ut中，并类似地定义了角节点位移向量ubl, ubr, utl, utr（图5）。一组单元格在整数平移njRj（j=1,2）上是周期性的，而Bloch–Floquet通过方程(26)允许在一个单元格平移Rj上发生相位移动来强制单个单元格的周期性。在离散形式中，Γr上的每个节点是Γl上相应节点的平移R1，而Γt上的每个节点是Γb上相应节点的平移R2。将方程(26)应用于这些相应的节点对，可以得到方程(31)的第一和第二个方程中总结的边界关系。同样，通过应用平移向量R2和R1+R2，我们可以将左上和右上节点与左下节点关联起来，如图31的第三和第四个方程所示。因此，对于任何给定的波矢k，Γr和Γt上的DOFs（以及相关的角DOFs）并不是独立的，可以用来表示ul, ub, ubl和自由内部DOFs ui。

为了系统地施加这些约束，使用了一种变换方法。构造了一个变换矩阵T来将完整的位移向量u与一组简化的自由度uk相关联，这些自由度对应于由k定义的简化基础：(32)u=T(k)uk，其中位移向量定义为(33)u=[ul;ub;ui;ur;ut;ubr;utl;utr]，uk=[ul;ub;ubl;ui]，ui表示自由节点的位移（在图5中显示为浅绿色点）。此外，变换矩阵T为(34)T(k)=I0000I00000Ieik?R1I0000eik?R2I0000eik?R1I000eik?R2I000eik?(R1+R2)I0。下载：下载高分辨率图像（174KB）下载：下载全尺寸图像图5. 带有节点索引和晶格向量指示的离散化单元格。

将简化的位移uk代入平衡方程和不稳定方程，并预先乘以TH（厄米转置），我们得到平衡方程和特征值问题的简化形式(35)Kμ,0χ0j=f0j, Kμ,0=T(0)HKμT(0),f0j=T(0)Hfj, (36)(Kμ,k?λμGμ,k)φk=0, Kμ,k=T(k)HKμT(k),Gμ,k=T(k)HGμT(k)，其中下标0表示特殊情况k=0，对应于没有相位移动的周期性边界条件。

3.4. 最坏情况模型
微观尺度下的屈曲敏感性不仅高度依赖于微观结构的体积分数和几何形状，还依赖于施加应力的类型和方向。为了获得一个不受应力方向影响的稳健的微观尺度稳定性度量，我们采用了Hübner等人[18]提出的最坏情况公式。因为线性化的屈曲问题与施加的应力成线性关系，所以BLF随宏观应力的大小成反比变化。因此，均质化的BLF完全由其在一个单元应力球S2上的值决定（如图6所示），该球包含所有可能的标准化应力方向σ?。然后可以通过用1/‖σ?‖进行重新缩放来获得任意应力状态σ?的BLF，即(37)λμk(σ?,α)=1‖σ?‖λ?μkσ?,α，其中λ?μk表示对于单元应力方向σ?∈S2评估的BLF，‖σ?‖是表示为Voigt符号的宏观应力张量的欧几里得范数。向量α代表确定结构体积分数的设计参数。这种标准化允许我们将研究限制在单位应力上，并在球面上表征响应，其中每个点对应于特定的应力类型和方向。

图6中所示的单元应力球S2提供了所有可能标准化应力状态的几何表示。球上的每个点可以用两个角度来描述。第一个是天顶角或纬度角，它衡量一个点位于极点和赤道之间的距离。这个角度决定了应力的类型：北极点对应于纯双轴压缩，南极点对应于纯双轴拉伸，赤道点对应于纯剪切，而中间纬度点对应于单轴或混合应力状态。第二个是方位角或经度角，它指定了应力张量相对于微观结构单元格的方向，捕捉不同载荷方向如何影响屈曲行为。通过这种方式，球形表示编码了应力的模式（即压缩、拉伸、剪切或混合）及其与晶格几何的对齐。

在以往的研究中，应力空间通常通过假设各向同性来简化，仅关注无剪切的双轴载荷并改变应力比σxx/σyy，或者通过检查沿选定方向的单轴载荷[43]。相比之下，在最坏情况模型中，计算了球S2上的整个屈曲屈服面，因此适用于具有各向异性屈曲行为的任意微观结构。

为了实施Floquet–Bloch理论[42]，在最坏情况模型中需要对布里渊区B内的波矢k进行最小化(38)λ?wc(α)=minσ?∈S2,k∈Bλ?μk(σ?,α)。使用方程(38)，我们确定最坏情况下的标准化应力方向σ?wc为在单元应力球S2上使Bloch–Floquet临界屈曲载荷因子最小化的那个。在我们的案例中，每个研究过的体积分数下的最小均质化载荷因子出现在单轴压缩应力状态下：对于空心三角形晶格沿x方向，对于等边三角形晶格沿y方向。因此，可以通过使用最坏情况应力状态σ?wc来获得与宏观应力σ?相关的最坏情况微观BLF(39)λwc(α,σ?)=1‖σ?‖mink∈Bλ?μk(σ?wc,α)。

方程(39)中的最坏情况屈曲模型通过最小化方程(38)中所有可能的宏观应力方向的均质化屈曲载荷因子，提供了一个稳健的、与模式无关的标准。与传统的微观尺度屈曲分析不同，后者需要在当前应力下解决特征值问题(41)，最坏情况公式通过消除对应力状态的依赖性来确保可导性和稳定性。这避免了对应力方向的敏感性，并消除了跟踪多个特征模式的必要性，使其在拓扑优化中特别有利，因为在优化过程中应力状态会发生变化。Hübner等人[18]最初为晶格微观结构提出的最坏情况公式可以应用于任何参数化的微观结构。

4. 并行多尺度优化框架
4.1. 多尺度分析
在由M个元素组成的宏观尺度域的有限元离散化中（图3），我们定义宏观尺度设计场η={η1,η2,…,ηM}∈(0,1]来表示域内的材料分布。每个值ηi决定了元素i内的局部材料存在情况。为了表征这些M个元素中的微观结构，我们引入了微观尺度设计向量x={x11,x21,…,xP1,x12,…,xP2,…,x1M,…,xPM}，其中每个微观架构由P个独立参数描述，从而可以单独定义和优化总共M个微观架构。为了方便起见，设计变量结合成α={η,x}。
固体各向同性材料惩罚（SIMP）方法[44],[45]通过将宏观尺度有限元模型与均质化的微观结构响应联系起来来抑制中间密度。对于元素e，有效刚度在接近空洞的基线Cmin和均质化张量CH(xe)之间进行插值，使用密度变量ηe∈[0,1]和惩罚因子p>1：(40)Ce(φe,xe)=Cmin+ηep(CH(xe)?Cmin)。这里，Cmin是一个小值的刚度张量，用于在ηe→0时避免奇异矩阵，p是强制材料存在或不存在的偏好因子（ηe≈1或0）；在这项研究中我们设置p=3。方程(40)中的均质化刚度张量CH可以通过直接数值均质化或通过预训练的替代模型获得。

使用像SIMP方法这样的替代材料模型在拓扑优化中解决特征值问题(40)是具有挑战性的，因为存在低体积分数区域，其刚度可以忽略不计。在这些区域，经常会出现非物理或人为的不稳定性，必须仔细区分这些不稳定性与真正的物理模式。
Neves等人[46]提出了一种解决这个问题的实际方法，即忽略低密度元素的几何刚度。遵循Gao等人[47]的方法，我们采用了一个密度截止值ηcut=0.3。对于ηe<ηcut的元素，方程(41)中的几何刚度贡献被设置为零，特征问题(41)中只保留线性刚度kμ?λμiφμ=0。这抑制了来自接近空洞区域的虚假模式，并提高了数值稳定性。 4.2. 优化问题 并行多尺度拓扑优化旨在通过同时优化宏观尺度材料分布η和潜在的微观结构设计参数x来提高结构性能。本研究使用基于深度神经网络（dnn）的替代模型来优化微结构的屈曲设计策略x。 为了优化结构的抗屈曲性，必须考虑blf，因为屈曲可能以几种不同的模式发生。仅关注最小的blf可能会忽略其他可能危及结构完整性的关键或接近临界模式。为此，采用了修改后的kreisselmeier–steinhauser（ks）函数将相关的blf聚合成一个单一的、平滑且可微的度量。修改后的ks函数定义为：(42)jksλ?i(α)=Λ1(α)+1ρln∑i=1qeρΛ?i(α)?Λ1(α)，其中ρ∈[1,∞)，q是要聚合的特征值数量，Λ1是逆屈曲载荷因子的上限，使得(43)Λ1(α)=maxi=1,…,qΛ?i(α)。这个函数的导数可以计算为(44)?JKSΛ?i(α)?αe=∑i=1qeρΛ?i(α)?Λ1(α)?Λ?i(α)?αe∑i=1qeρΛ?i(α)?Λ1(α)。 在屈曲优化的背景下，必须考虑两种不同的现象：宏观尺度的屈曲和微观尺度的屈曲。因此，为了实施该框架，我们提出了三个优化问题： •问题1. 在宏观和微观尺度上都受体积约束的情况下，优化宏观尺度屈曲。 •问题2. 在宏观和微观尺度上都受体积约束的情况下，优化宏观和微观尺度的屈曲。 •问题3. 在宏观和微观尺度上都受顺从性和体积约束的情况下，优化两种尺度上的屈曲。 在第一个优化问题中，为了改善宏观尺度的屈曲性能，我们最小化了宏观blf的倒数之和，记为λinv。相应的目标函数和约束条件为(45)minαgw(α)=JKS(Λinv(α)), gv=1M∑e=1Mηe?Vη,1M∑e=1Mηevx?Vx≤0,>< /><αimax，其中gw是目标函数，gv是体积约束。vx, vη, 和 vx 分别表示微观结构的体积、宏观和微观尺度上的体积分数：vη∈[0,1]表示宏观尺度域被材料占据的比例（0=空洞，1 = 固体），vx是微观结构特征的平均体积分数。因此，设计域中的净体积分数为 v=VηVx。最小允许的微观尺度体积分数由所选择的参数化确定，并受到几何限制>< /><αimax。 结构可能在任何全局（即宏观）不稳定发生之前经历局部（即微观）屈曲。因此，在优化过程中考虑两个尺度上的屈曲是必要的。在第二个屈曲优化问题中，我们的目标是通过最大化宏观尺度屈曲抵抗来提高整体结构稳定性，同时抑制微观尺度屈曲。不失一般性，我们假设通过将两个尺度结合到一个等权重的目标函数中来实现平衡。此外，通过考虑更高的权重水平（0.75）并比较 hollow 晶格的结果blf值和优化后的结构来检验对权重参数的敏感性。与第一个优化问题类似，方程(39)中的均质化blf的倒数，记为λinv，以及宏观blf的倒数之和被定义为要最小化的目标函数。相应的问题表述为(46)minαgw(α)=ω∑i=1Mλinvi(α)λ0+(1?ω)JKS(Λinv(α))J0KS,ω=0.5gV=1M∑e=1Mηe?Vη,1M∑e=1Mηevx?Vx≤0,>< /><αimax，其中ω是权重的值。j0ks和λ0是归一化特征值的常数值。
除了提高屈曲抵抗性外，优化后的结构还应该表现出更大的刚度。>除了提高屈曲抵抗性外，优化后的结构还应该表现出更大的刚度。>为了解决这个问题，必须将柔顺性（即刚度的倒数）纳入优化框架中。因此，从第二个问题获得的优化结构的柔顺性值被用作参考。在第三个优化问题中，柔顺性被引入作为一个约束条件，并且我们不失一般性地假设第二个问题的柔顺性值的70%是上限。此外，为了研究柔顺性上限的影响，使用空心晶格微结构对同一个问题分别以60%和70%的两个值进行了求解。这种公式化旨在产生更刚性的结构，同时仍然提高抗屈曲性能。因此，带有柔顺性约束的优化问题表示为（47）minαgw(α)=ω∑i=1Mλinvi(α)λ0+(1?ω)JKS(Λinv(α))J0KS,ω=0.5gV=1M∑e=1Mηe?Vη,1M∑e=1Mηevx?Vx≤0,gc=FTuCe??1≤0,s.t.αimin< /><αimax，其中gc和ce?分别是柔顺性约束函数和上限。为了确保解与网格无关并施加最小长度尺度，对设计变量应用了空间过滤技术[48]。对于一个元素i和过滤半径r，所有质心距离i在r以内的元素j都被包括在过滤范围内。过滤后的设计变量α?i是使用标准化的高斯加权核计算的（48）α?i=∑jωfjωfαj,（49）ωfj=max1?djr,0,（50）ωf=∑jωfj。这里，dj是元素i和j的质心之间的距离，ωfj是相应的核权重。这种过滤方法在两个尺度上平滑了几何形状，减轻了微结构中的突然局部变化，并提高了优化设计的可制造性和鲁棒性。为了提高优化设计的清晰度和可制造性，应该应用heaviside投影滤波器。平滑后的heaviside函数hm(η)由wang等人[49]定义（51）η?=hm(η)=tanh(βκm)+tanhβ(η?κm)tanh(βκm)+tanhβ(1?κm)，其中κm是阈值。在这项工作中，阈值被认为是0.5。投影的陡峭程度由β决定。当β增加时，heaviside函数变得非线性更强。因此，对β使用了一种连续方法，即初始化为β=1，并在优化过程中逐渐增加到最大值10以增强投影效果。

4.3. 体积约束
在所有优化问题中，净体积上限被认为是0.6。在第一种情况下，通过在宏观尺度和微观尺度上设置上限来应用体积约束。在第二种情况下，我们将微观尺度的体积分数上限设置为1，允许使用完全实心的单元格，同时保持宏观尺度的体积约束有效。在第三种情况下，宏观尺度的体积约束上限设置为1，有效地移除了全局体积限制，允许材料在所有地方存在。在最后一种情况下，对净体积施加单一约束，唯一的设计变量是支柱的厚度（情况4）。所有不同的体积约束在表1中总结。

表1. 每种情况下的宏观尺度（vη）、微观尺度（vx）和净体积（v）的体积分数。

4. 敏感性分析
4.4.1. 宏观尺度屈曲
宏观尺度blf的敏感性分析已经由rodrigues等人[50]进行了彻底讨论。因此，宏观尺度特征值倒数对设计变量α的敏感性为（52）?λinv?α=??t?k?α+λinv?g?α?+λinvwadjt?k?αu，其中wadj表示通过求解（53）k(α)wadj=?t?ug(α,u(α))?获得的宏观尺度屈曲的伴随向量，（54）?ug=?e=1m∑k=1nipce,kb?e,kt?σe?ue,ijb?e,k，以及（55）?σe?ue,ij=1nip∑k=1nipchbe,keij。通过（56）?ke?αi=∫ωebt?ce?αibdωe，可以获得设计变量对刚度矩阵的梯度（52），其中元素的有效弹性张量ce可以使用方程（40）获得。弹性张量对αi={ηi,xi}的梯度可以表示为（57）?ce?αe=?ce?ηe,?ce?xe1,?ce?xe2,…,?ce?xep，其中（58）?ce?ηe=cmin+pη?ep?1dηe?dηech(xe)?cmin，（59）dη?edηe=βsech2β(ηe?κm)tanh(βκm)+tanhβ(1?κm），以及（60）?ce?xei=cmin+η?ep?ch?xei?cmin,i=1,2,…,p。微观尺度的敏感性通过评估?ch/?xei进入分析，描述了均质化材料属性如何随微观结构参数变化。

在应用设计过滤后，如方程（48）所述，过滤后的设计变量α?通过权重矩阵gw与原始设计变量α相关联（61）α?=gwα。通过（62）?α??α=gwt恢复了对未过滤设计变量的敏感性。方程（52）中应力刚度矩阵对设计变量的梯度由（63）?g?αi=?e=1m∑k=1nipce,kb?e,k?σe(ue(α))?αib?e,k给出，其中（64）?σe(ue(α))?αi=1nip∑k=1nip?ce?αibe,kue(α)。

4.4.2. 微观尺度屈曲
可以使用基于伴随方法[51]计算特征值λinv的倒数敏感性（65）??dλinvdαe=φμh?gμ?αe+?gμ???d??dαe?λinv?kμ?αeφμ+∑j=13(t(0)vmj)h?kμ?αeχj??fj?αe，其中vmj表示通过解决下面的伴随问题获得的微观尺度屈曲的伴随变量。（66）kμ0vmj=?t(0)hφμh?gμ?χjφμh,j=1,2,3。首先在元素级别评估伴随载荷，然后将其组装到全局系统中。以展开形式表示，该过程为（67）φμh?gμ?χjφμ=∑e=1nφμeh?gμe?χe,1jφμe?φμeh?gμe?χe,nnjφμe，其中nn表示元素e的自由度数。元素矩阵gμe对χej的局部导数为（68）?gμe?χe,ij=∫yeb?e,lt?σlm,e?χe,ijb?e,mdy，其中（69）??σe?χe,ij=?cebeeij??，其中eij是一个在位置（i,j）处为1的8 × 3单元素矩阵。当结构受到规定的应变???时，项???gμ???d??dαe消失，因为???与设计变量无关。

4.4.3. 柔顺性
关于柔顺性敏感性的完整推导和实现细节可以在black和najafi[32]的工作中找到。设x表示宏观尺度有限元网格中的节点坐标向量。在这项研究中，网格是固定的，不依赖于设计变量（即?x/?α=?f/?α=0）。在这些假设下，通过伴随方法简化了柔顺性函数对设计变量αi的导数（70）?gc?αi=λadjt?k?αiu，其中λadj是柔顺性的伴随向量，定义为（71）λadjtk=f。

下载：下载高分辨率图像（220kb）
下载：下载完整尺寸图像

图7. （a）所有支柱厚度相同的空心三角形晶格架构的几何形状，以及（b）mape和（c）rmse在训练和验证期间的收敛图。

4.4.4. 用于设计敏感性的dnn
方程（10）不仅实现了函数逼近，还实现了敏感性调整，这在优化问题中至关重要。这项工作建立在基础工作之上，该工作广泛研究了通过增加网络深度对敏感性表示、对低保真度或噪声训练数据的鲁棒性以及具有高度参数化微结构输入的性能的影响；这些方面已经得到了全面的讨论，特别是在数值均质化的替代建模中[32]。主要目标是通过联合改进宏观尺度和微观尺度的屈曲载荷因子来最大化抗屈曲性能，这两个因子都是均质化本构矩阵ch的函数（在前面的小节中描述）。由于这些载荷因子的敏感性取决于梯度?ch（例如方程（52），（65）），因此可以将在优化的循环中嵌入训练好的dnn来提供这些量并降低计算成本。优化需要从局部设计变量xe到均质化刚度张量ch(xe)的显式映射。为此，我们使用在有限元参考数据上训练的dnn代理来近似c?h。因此，通过最小化方程（9）来校准方程（8）中定义的网络（见图2），以映射xe?c?h。因为宏观和微观尺度的屈曲载荷因子都依赖于ch，所以优化需要敏感性?c?h/?xei，我们通过链式法则（72）?c?h?xei=?c?h?al?al?tl?tl?al?1??t1?xei来获得这些敏感性。

5. dnn训练
为了实现足够的准确性和稳定性，我们采用了具有三个隐藏层和每层32个神经元的全连接前馈神经网络，遵循black和najafi[32]的方法。图7(a)、图8(a)和图9(a)展示了用于训练的三种不同晶格架构的几何形状。图8(a)和图9(a)之间的主要区别在于设计参数的数量。在图8(a)所示的微结构中，整个结构使用一个设计参数。相比之下，图9(a)中的微结构允许支柱厚度具有空间变化，从而产生五个独立的设计变量。图7(b)–(c)、图8(b)–(c)和图9(b)–(c)分别显示了训练和验证阶段的平均绝对百分比误差（mape）和均方根误差（rmse）的相应收敛图。在每种情况下，训练和验证误差在初始时期迅速减小，然后趋于平稳，表明学习了成功，没有过拟合的迹象。值得注意的是，mape和rmse都收敛到较低的值，证明模型能够准确捕捉不同晶格结构的复杂行为。机器学习加速使得可以使用具有更多设计参数的更复杂的微结构，同时保持两尺度优化的计算可行性。除了晶格单元格外，我们还为（i）具有三个几何设计参数（主轴、副轴和旋转角度）的椭圆形微结构和（ii）具有八个几何设计参数的生物桁架微结构训练了替代模型。图10和图11分别报告了椭圆形和生物桁架的训练历史的rmse和mape。在两种情况下，误差在早期阶段急剧下降，然后趋于平稳，表明快速学习了主导映射，随后逐渐细化。在后期阶段观察到rmse和mape的第二次下降（与额外的训练细化如学习率调度一致），这对于高维生物桁架参数化特别有益。这些历史记录证实了即使在微结构设计参数数量增加的情况下，替代训练也是稳定和准确的。

下载：下载高分辨率图像（242kb）
下载：下载完整尺寸图像

图8. （a）所有支柱厚度相同的等边三角形晶格架构的几何形状，以及（b）mape和（c）rmse在训练和验证期间的收敛图。

下载：下载高分辨率图像（240kb）
下载：下载完整尺寸图像

图9. （a）支柱厚度不同的等边三角形晶格架构的几何形状，以及（b）mape和（c）rmse在训练和验证期间的收敛图。

为了比较有限元方法和基于dnn的预测在不同晶格上的计算时间，使用配备intel core i7（ultra 9 185h）和32 gb ddr5 ram的单台pc测量了cpu时间。在保持可比准确性的同时，dnn模型的推理时间比有限元模拟快几个数量级（大约快104倍），突出了它们作为实时或大规模模拟的替代模型的有效性（图12）。值得注意的是，每个微结构模型可以在20分钟内完成训练。

下载：下载高分辨率图像（258kb）
下载：下载完整尺寸图像

图10. 椭圆微结构的训练历史（3个设计参数）：rmse和mape与训练和验证集的训练周期。

下载：下载高分辨率图像（271kb）
下载：下载完整尺寸图像

图11. 生物桁架微结构的训练历史（8个设计参数）：rmse和mape与训练和验证集的训练周期。

下载：下载高分辨率图像（171kb）
下载：下载完整尺寸图像

图12. 有限元方法（fem）和dnn对不同晶格的均质化材料属性预测的计算时间比较。

6. 优化示例
6.1. 问题定义
我们考虑图13所示的悬臂结构，在宏观尺度上使用6 × 36的粗糙网格和10 × 60的细密线性4节点有限元进行离散化。总载荷的大小为1，在梁的右侧边缘均匀分布。域的几何形状由长度l=6和宽度w=1定义。在这项工作中，通过将杨氏模量设置为e=1和泊松比设置为ν=0.3来无量纲化材料属性。在微观尺度上，采用了两种不同的晶格微结构。第4.2节介绍的三个优化问题应用于这种多尺度设置，所得设计在后续章节中进行了讨论。图14显示了在体积约束下使用空心三角形晶格微结构的宏观尺度blf的多尺度拓扑优化的代表结果，说明了在所提出的框架中材料在宏观和微观尺度上的分布。图14(a)显示了使用均质密度场获得的优化宏观尺度布局；图14(b)描绘了在选定位置的代表性微结构实现；图14(c)展示了在整个域内实例化的完全非均质化设计。接下来，我们解决几个优化问题并比较了所得到的优化结构。 4.3. 体积约束 在所有优化问题中，净体积上限被认为是0.6。在第一种情况下，通过在宏观尺度和微观尺度上设置上限来应用体积约束。在第二种情况下，我们将微观尺度的体积分数上限设置为1，允许使用完全实心的单元格，同时保持宏观尺度的体积约束有效。在第三种情况下，宏观尺度的体积约束上限设置为1，有效地移除了全局体积限制，允许材料在所有地方存在。在最后一种情况下，对净体积施加单一约束，唯一的设计变量是支柱的厚度（情况4）。所有不同的体积约束在表1中总结。 表1. 每种情况下的宏观尺度（vη）、微观尺度（vx）和净体积（v）的体积分数。 4. 敏感性分析 4.4.1. 宏观尺度屈曲 宏观尺度blf的敏感性分析已经由rodrigues等人[50]进行了彻底讨论。因此，宏观尺度特征值倒数对设计变量α的敏感性为（52）?λinv?α=??T?K?α+Λinv?G?α?+ΛinvwadjT?K?αu，其中wadj表示通过求解（53）K(α)wadj=?T?uG(α,u(α))?获得的宏观尺度屈曲的伴随向量，（54）?uG=?e=1M∑k=1Nipce,kB?e,kT?σe?ue,ijB?e,k，以及（55）?σe?ue,ij=1Nip∑k=1NipCHBe,keij。通过（56）?ke?αi=∫ΩeBT?Ce?αiBdΩe，可以获得设计变量对刚度矩阵的梯度（52），其中元素的有效弹性张量Ce可以使用方程（40）获得。弹性张量对αi={ηi,xi}的梯度可以表示为（57）?Ce?αe=?Ce?ηe,?Ce?xe1,?Ce?xe2,…,?Ce?xep，其中（58）?Ce?ηe=Cmin+pη?ep?1dηe?dηeCH(xe)?Cmin，（59）dη?edηe=βsech2β(ηe?κm)tanh(βκm)+tanhβ(1?κm），以及（60）?Ce?xei=Cmin+η?ep?CH?xei?Cmin,i=1,2,…,P。微观尺度的敏感性通过评估?CH/?xei进入分析，描述了均质化材料属性如何随微观结构参数变化。 在应用设计过滤后，如方程（48）所述，过滤后的设计变量α?通过权重矩阵gw与原始设计变量α相关联（61）α?=Gwα。通过（62）?α??α=GwT恢复了对未过滤设计变量的敏感性。方程（52）中应力刚度矩阵对设计变量的梯度由（63）?G?αi=?e=1M∑k=1Nipce,kB?e,k?σe(ue(α))?αiB?e,k给出，其中（64）?σe(ue(α))?αi=1Nip∑k=1Nip?Ce?αiBe,kue(α)。 4.4.2. 微观尺度屈曲 可以使用基于伴随方法[51]计算特征值λinv的倒数敏感性（65）??dλinvdαe=φμH?Gμ?αe+?Gμ???d??dαe?λinv?Kμ?αeφμ+∑j=13(T(0)vMj)H?Kμ?αeχj??fj?αe，其中vMj表示通过解决下面的伴随问题获得的微观尺度屈曲的伴随变量。（66）Kμ0vMj=?T(0)HφμH?Gμ?χjφμH,j=1,2,3。首先在元素级别评估伴随载荷，然后将其组装到全局系统中。以展开形式表示，该过程为（67）φμH?Gμ?χjφμ=∑e=1NφμeH?Gμe?χe,1jφμe?φμeH?Gμe?χe,nnjφμe，其中nn表示元素e的自由度数。元素矩阵Gμe对χej的局部导数为（68）?Gμe?χe,ij=∫YeB?e,lT?σlm,e?χe,ijB?e,mdY，其中（69）??σe?χe,ij=?CeBeeij??，其中eij是一个在位置（i,j）处为1的8 × 3单元素矩阵。当结构受到规定的应变???时，项???gμ???d??dαe消失，因为???与设计变量无关。 4.4.3. 柔顺性 关于柔顺性敏感性的完整推导和实现细节可以在black和najafi[32]的工作中找到。设x表示宏观尺度有限元网格中的节点坐标向量。在这项研究中，网格是固定的，不依赖于设计变量（即?x ?α=?F/?α=0）。在这些假设下，通过伴随方法简化了柔顺性函数对设计变量αi的导数（70）?gc?αi=λadjT?K?αiu，其中λadj是柔顺性的伴随向量，定义为（71）λadjTK=f。 下载：下载高分辨率图像（220kb） 下载：下载完整尺寸图像 图7. （a）所有支柱厚度相同的空心三角形晶格架构的几何形状，以及（b）mape和（c）rmse在训练和验证期间的收敛图。 4.4.4. 用于设计敏感性的dnn 方程（10）不仅实现了函数逼近，还实现了敏感性调整，这在优化问题中至关重要。这项工作建立在基础工作之上，该工作广泛研究了通过增加网络深度对敏感性表示、对低保真度或噪声训练数据的鲁棒性以及具有高度参数化微结构输入的性能的影响；这些方面已经得到了全面的讨论，特别是在数值均质化的替代建模中[32]。主要目标是通过联合改进宏观尺度和微观尺度的屈曲载荷因子来最大化抗屈曲性能，这两个因子都是均质化本构矩阵ch的函数（在前面的小节中描述）。由于这些载荷因子的敏感性取决于梯度?ch（例如方程（52），（65）），因此可以将在优化的循环中嵌入训练好的dnn来提供这些量并降低计算成本。优化需要从局部设计变量xe到均质化刚度张量ch(xe)的显式映射。为此，我们使用在有限元参考数据上训练的dnn代理来近似c?h。因此，通过最小化方程（9）来校准方程（8）中定义的网络（见图2），以映射xe?c?h。因为宏观和微观尺度的屈曲载荷因子都依赖于ch，所以优化需要敏感性?c?h ?xei，我们通过链式法则（72）?c?h?xei=?C?H?AL?AL?TL?TL?AL?1??T1?xei来获得这些敏感性。 5. dnn训练 为了实现足够的准确性和稳定性，我们采用了具有三个隐藏层和每层32个神经元的全连接前馈神经网络，遵循black和najafi[32]的方法。图7(a)、图8(a)和图9(a)展示了用于训练的三种不同晶格架构的几何形状。图8(a)和图9(a)之间的主要区别在于设计参数的数量。在图8(a)所示的微结构中，整个结构使用一个设计参数。相比之下，图9(a)中的微结构允许支柱厚度具有空间变化，从而产生五个独立的设计变量。图7(b)–(c)、图8(b)–(c)和图9(b)–(c)分别显示了训练和验证阶段的平均绝对百分比误差（mape）和均方根误差（rmse）的相应收敛图。在每种情况下，训练和验证误差在初始时期迅速减小，然后趋于平稳，表明学习了成功，没有过拟合的迹象。值得注意的是，mape和rmse都收敛到较低的值，证明模型能够准确捕捉不同晶格结构的复杂行为。机器学习加速使得可以使用具有更多设计参数的更复杂的微结构，同时保持两尺度优化的计算可行性。除了晶格单元格外，我们还为（i）具有三个几何设计参数（主轴、副轴和旋转角度）的椭圆形微结构和（ii）具有八个几何设计参数的生物桁架微结构训练了替代模型。图10和图11分别报告了椭圆形和生物桁架的训练历史的rmse和mape。在两种情况下，误差在早期阶段急剧下降，然后趋于平稳，表明快速学习了主导映射，随后逐渐细化。在后期阶段观察到rmse和mape的第二次下降（与额外的训练细化如学习率调度一致），这对于高维生物桁架参数化特别有益。这些历史记录证实了即使在微结构设计参数数量增加的情况下，替代训练也是稳定和准确的。 下载：下载高分辨率图像（242kb） 下载：下载完整尺寸图像 图8. （a）所有支柱厚度相同的等边三角形晶格架构的几何形状，以及（b）mape和（c）rmse在训练和验证期间的收敛图。 下载：下载高分辨率图像（240kb） 下载：下载完整尺寸图像 图9. （a）支柱厚度不同的等边三角形晶格架构的几何形状，以及（b）mape和（c）rmse在训练和验证期间的收敛图。 为了比较有限元方法和基于dnn的预测在不同晶格上的计算时间，使用配备intel core i7（ultra 9 185h）和32 gb ddr5 ram的单台pc测量了cpu时间。在保持可比准确性的同时，dnn模型的推理时间比有限元模拟快几个数量级（大约快104倍），突出了它们作为实时或大规模模拟的替代模型的有效性（图12）。值得注意的是，每个微结构模型可以在20分钟内完成训练。 下载：下载高分辨率图像（258kb） 下载：下载完整尺寸图像 图10. 椭圆微结构的训练历史（3个设计参数）：rmse和mape与训练和验证集的训练周期。 下载：下载高分辨率图像（271kb） 下载：下载完整尺寸图像 图11. 生物桁架微结构的训练历史（8个设计参数）：rmse和mape与训练和验证集的训练周期。 下载：下载高分辨率图像（171kb） 下载：下载完整尺寸图像 图12. 有限元方法（fem）和dnn对不同晶格的均质化材料属性预测的计算时间比较。 6. 优化示例 6.1. 问题定义 我们考虑图13所示的悬臂结构，在宏观尺度上使用6 × 36的粗糙网格和10 × 60的细密线性4节点有限元进行离散化。总载荷的大小为1，在梁的右侧边缘均匀分布。域的几何形状由长度l=>
4.3. 体积约束
在所有优化问题中，净体积上限被认为是0.6。在第一种情况下，通过在宏观尺度和微观尺度上设置上限来应用体积约束。在第二种情况下，我们将微观尺度的体积分数上限设置为1，允许使用完全实心的单元格，同时保持宏观尺度的体积约束有效。在第三种情况下，宏观尺度的体积约束上限设置为1，有效地移除了全局体积限制，允许材料在所有地方存在。在最后一种情况下，对净体积施加单一约束，唯一的设计变量是支柱的厚度（情况4）。所有不同的体积约束在表1中总结。

表1. 每种情况下的宏观尺度（vη）、微观尺度（vx）和净体积（v）的体积分数。

4. 敏感性分析
4.4.1. 宏观尺度屈曲
宏观尺度blf的敏感性分析已经由rodrigues等人[50]进行了彻底讨论。因此，宏观尺度特征值倒数对设计变量α的敏感性为（52）?λinv?α=??t?k?α+λinv?g?α?+λinvwadjt?k?αu，其中wadj表示通过求解（53）k(α)wadj=?t?ug(α,u(α))?获得的宏观尺度屈曲的伴随向量，（54）?ug=?e=1m∑k=1nipce,kb?e,kt?σe?ue,ijb?e,k，以及（55）?σe?ue,ij=1nip∑k=1nipchbe,keij。通过（56）?ke?αi=∫ωebt?ce?αibdωe，可以获得设计变量对刚度矩阵的梯度（52），其中元素的有效弹性张量ce可以使用方程（40）获得。弹性张量对αi={ηi,xi}的梯度可以表示为（57）?ce?αe=?ce?ηe,?ce?xe1,?ce?xe2,…,?ce?xep，其中（58）?ce?ηe=cmin+pη?ep?1dηe?dηech(xe)?cmin，（59）dη?edηe=βsech2β(ηe?κm)tanh(βκm)+tanhβ(1?κm），以及（60）?ce?xei=cmin+η?ep?ch?xei?cmin,i=1,2,…,p。微观尺度的敏感性通过评估?ch/?xei进入分析，描述了均质化材料属性如何随微观结构参数变化。

在应用设计过滤后，如方程（48）所述，过滤后的设计变量α?通过权重矩阵gw与原始设计变量α相关联（61）α?=gwα。通过（62）?α??α=gwt恢复了对未过滤设计变量的敏感性。方程（52）中应力刚度矩阵对设计变量的梯度由（63）?g?αi=?e=1m∑k=1nipce,kb?e,k?σe(ue(α))?αib?e,k给出，其中（64）?σe(ue(α))?αi=1nip∑k=1nip?ce?αibe,kue(α)。

4.4.2. 微观尺度屈曲
可以使用基于伴随方法[51]计算特征值λinv的倒数敏感性（65）??dλinvdαe=φμh?gμ?αe+?gμ???d??dαe?λinv?kμ?αeφμ+∑j=13(t(0)vmj)h?kμ?αeχj??fj?αe，其中vmj表示通过解决下面的伴随问题获得的微观尺度屈曲的伴随变量。（66）kμ0vmj=?t(0)hφμh?gμ?χjφμh,j=1,2,3。首先在元素级别评估伴随载荷，然后将其组装到全局系统中。以展开形式表示，该过程为（67）φμh?gμ?χjφμ=∑e=1nφμeh?gμe?χe,1jφμe?φμeh?gμe?χe,nnjφμe，其中nn表示元素e的自由度数。元素矩阵gμe对χej的局部导数为（68）?gμe?χe,ij=∫yeb?e,lt?σlm,e?χe,ijb?e,mdy，其中（69）??σe?χe,ij=?cebeeij??，其中eij是一个在位置（i,j）处为1的8 × 3单元素矩阵。当结构受到规定的应变???时，项???gμ???d??dαe消失，因为???与设计变量无关。

4.4.3. 柔顺性
关于柔顺性敏感性的完整推导和实现细节可以在black和najafi[32]的工作中找到。设x表示宏观尺度有限元网格中的节点坐标向量。在这项研究中，网格是固定的，不依赖于设计变量（即?x/?α=?f/?α=0）。在这些假设下，通过伴随方法简化了柔顺性函数对设计变量αi的导数（70）?gc?αi=λadjt?k?αiu，其中λadj是柔顺性的伴随向量，定义为（71）λadjtk=f。

下载：下载高分辨率图像（220kb）
下载：下载完整尺寸图像

图7. （a）所有支柱厚度相同的空心三角形晶格架构的几何形状，以及（b）mape和（c）rmse在训练和验证期间的收敛图。

4.4.4. 用于设计敏感性的dnn
方程（10）不仅实现了函数逼近，还实现了敏感性调整，这在优化问题中至关重要。这项工作建立在基础工作之上，该工作广泛研究了通过增加网络深度对敏感性表示、对低保真度或噪声训练数据的鲁棒性以及具有高度参数化微结构输入的性能的影响；这些方面已经得到了全面的讨论，特别是在数值均质化的替代建模中[32]。主要目标是通过联合改进宏观尺度和微观尺度的屈曲载荷因子来最大化抗屈曲性能，这两个因子都是均质化本构矩阵ch的函数（在前面的小节中描述）。由于这些载荷因子的敏感性取决于梯度?ch（例如方程（52），（65）），因此可以将在优化的循环中嵌入训练好的dnn来提供这些量并降低计算成本。优化需要从局部设计变量xe到均质化刚度张量ch(xe)的显式映射。为此，我们使用在有限元参考数据上训练的dnn代理来近似c?h。因此，通过最小化方程（9）来校准方程（8）中定义的网络（见图2），以映射xe?c?h。因为宏观和微观尺度的屈曲载荷因子都依赖于ch，所以优化需要敏感性?c?h/?xei，我们通过链式法则（72）?c?h?xei=?c?h?al?al?tl?tl?al?1??t1?xei来获得这些敏感性。

5. dnn训练
为了实现足够的准确性和稳定性，我们采用了具有三个隐藏层和每层32个神经元的全连接前馈神经网络，遵循black和najafi[32]的方法。图7(a)、图8(a)和图9(a)展示了用于训练的三种不同晶格架构的几何形状。图8(a)和图9(a)之间的主要区别在于设计参数的数量。在图8(a)所示的微结构中，整个结构使用一个设计参数。相比之下，图9(a)中的微结构允许支柱厚度具有空间变化，从而产生五个独立的设计变量。图7(b)–(c)、图8(b)–(c)和图9(b)–(c)分别显示了训练和验证阶段的平均绝对百分比误差（mape）和均方根误差（rmse）的相应收敛图。在每种情况下，训练和验证误差在初始时期迅速减小，然后趋于平稳，表明学习了成功，没有过拟合的迹象。值得注意的是，mape和rmse都收敛到较低的值，证明模型能够准确捕捉不同晶格结构的复杂行为。机器学习加速使得可以使用具有更多设计参数的更复杂的微结构，同时保持两尺度优化的计算可行性。除了晶格单元格外，我们还为（i）具有三个几何设计参数（主轴、副轴和旋转角度）的椭圆形微结构和（ii）具有八个几何设计参数的生物桁架微结构训练了替代模型。图10和图11分别报告了椭圆形和生物桁架的训练历史的rmse和mape。在两种情况下，误差在早期阶段急剧下降，然后趋于平稳，表明快速学习了主导映射，随后逐渐细化。在后期阶段观察到rmse和mape的第二次下降（与额外的训练细化如学习率调度一致），这对于高维生物桁架参数化特别有益。这些历史记录证实了即使在微结构设计参数数量增加的情况下，替代训练也是稳定和准确的。

下载：下载高分辨率图像（242kb）
下载：下载完整尺寸图像

图8. （a）所有支柱厚度相同的等边三角形晶格架构的几何形状，以及（b）mape和（c）rmse在训练和验证期间的收敛图。

下载：下载高分辨率图像（240kb）
下载：下载完整尺寸图像

图9. （a）支柱厚度不同的等边三角形晶格架构的几何形状，以及（b）mape和（c）rmse在训练和验证期间的收敛图。

为了比较有限元方法和基于dnn的预测在不同晶格上的计算时间，使用配备intel core i7（ultra 9 185h）和32 gb ddr5 ram的单台pc测量了cpu时间。在保持可比准确性的同时，dnn模型的推理时间比有限元模拟快几个数量级（大约快104倍），突出了它们作为实时或大规模模拟的替代模型的有效性（图12）。值得注意的是，每个微结构模型可以在20分钟内完成训练。

下载：下载高分辨率图像（258kb）
下载：下载完整尺寸图像

图10. 椭圆微结构的训练历史（3个设计参数）：rmse和mape与训练和验证集的训练周期。

下载：下载高分辨率图像（271kb）
下载：下载完整尺寸图像

图11. 生物桁架微结构的训练历史（8个设计参数）：rmse和mape与训练和验证集的训练周期。

下载：下载高分辨率图像（171kb）
下载：下载完整尺寸图像

图12. 有限元方法（fem）和dnn对不同晶格的均质化材料属性预测的计算时间比较。

6. 优化示例
6.1. 问题定义
我们考虑图13所示的悬臂结构，在宏观尺度上使用6 × 36的粗糙网格和10 × 60的细密线性4节点有限元进行离散化。总载荷的大小为1，在梁的右侧边缘均匀分布。域的几何形状由长度l=6和宽度w=1定义。在这项工作中，通过将杨氏模量设置为e=1和泊松比设置为ν=0.3来无量纲化材料属性。在微观尺度上，采用了两种不同的晶格微结构。第4.2节介绍的三个优化问题应用于这种多尺度设置，所得设计在后续章节中进行了讨论。图14显示了在体积约束下使用空心三角形晶格微结构的宏观尺度blf的多尺度拓扑优化的代表结果，说明了在所提出的框架中材料在宏观和微观尺度上的分布。图14(a)显示了使用均质密度场获得的优化宏观尺度布局；图14(b)描绘了在选定位置的代表性微结构实现；图14(c)展示了在整个域内实例化的完全非均质化设计。接下来，我们解决几个优化问题并比较了所得到的优化结构。>多尺度拓扑优化采用等边晶格微观结构，以最大化宏观尺度屈曲载荷因子（BLF）：(a) 具有均匀密度场的优化宏观布局；(b) 在选定位置处的代表性微观结构实例；以及 (c) 在整个域中应用微观结构的完全非均匀化设计。下载：下载高分辨率图像（420KB）下载：下载全尺寸图像

图15. (a) 在不同体积约束下，通过宏观尺度屈曲优化获得的由空心三角形晶格微观结构组成的优化结构；(b) 不同体积约束下的收敛历史图。下载：下载高分辨率图像（548KB）下载：下载全尺寸图像

图16. (a) 在不同体积约束下，通过宏观尺度屈曲优化获得的由等边三角形晶格微观结构组成的优化结构；(b) 不同体积约束下的收敛历史图。

6.2. 宏观尺度屈曲优化

方程（45）中的宏观尺度优化问题是针对表1中的体积分割进行的。空心三角形微观结构的优化布局及其历史图显示在图15中。由于外骨骼的惯性矩大于内骨骼[18]，所有解决方案都显示出类似外骨骼的排列。净体积固定为V=0.60，但改变分割（Vη, Vx）（其中Vη和Vx分别表示宏观和微观尺度上的体积分数）会重新分配尺度之间的材料，并改变哪种屈曲模式被最有效地抑制。如图15所示，在宏观主导的分割（Vη=1.00, Vx=0.60）中，设计形成了更连续的宏观尺度构件，而局部晶格变得更薄，这延迟了整体屈曲，但可能会允许更多的微观级不稳定性[18]。在微观主导的分割（Vη=0.60, Vx=1.00）中，材料在单元格内集中，主要构件变得更厚，连接部分变窄（图15）。这增强了局部单元格的强度，但可能会降低主要构件的惯性矩[52]。平衡分割（Vη=0.75, Vx=0.80）产生了一个类似外骨骼的框架，单元格基本均匀，材料在开口周围和中间深度处局部加厚，形成了加固区域，从而在全局连续性和局部强化之间取得了平衡（图15）。因此，增加Vη可以增强宏观尺度上的连续性，而增加Vx则可以增强局部结构的密度。

使用表1中的体积分割，图16显示了等边晶格的优化布局，并可以与图15中的空心晶格进行比较。在净体积固定为V=0.60的情况下，布局取决于材料在尺度之间的分配方式。当宏观尺度体积分数的上限较大（Vη=1.00, Vx=0.60）时，设计将材料放置在外边界，使得外部构件连续且相互连接，而内部的三角形单元格变得更薄。如图16所示，当微观尺度体积分数的上限较大（Vη=0.60, Vx=1.00）时，材料在单元格内集中；三角形筋条变得更厚，宏观区域之间的连接变窄，这增强了局部响应，但中断了构件之间的连续性。平衡分割（Vη=0.75, Vx=0.80）产生了一个带状排列，其中筋条的厚度在宽度上逐渐变化（图16）。与图15中的空心晶格相比，图16中显示的等边单元格保持了更直、更重复的筋条行，内部空洞较少，因此当宏观尺度体积可用时，它更容易形成清晰的框架，并在设计偏向微观尺度时产生更明显的过渡。在仅限制净体积的情况下（即情况4），有效设计的变化仅限于调整筋条的厚度。由于无法在宏观和微观尺度之间转移体积，优化器使用可用材料来加厚已经出现的外部构件，当宏观尺度体积上限为一时。结果模式与图16中的宏观主导布局非常吻合，优化后的屈曲载荷因子与图17中报告的宏观主导值相同。

相应的屈曲载荷因子在图17中绘制。对于空心三角形晶格，图中的顺序是：平衡分割（情况1）最低，微观主导（情况2）较高，宏观主导（情况3）最高。图17中的这一趋势与图15中的布局一致：将更多体积分配给宏观尺度可以产生更连续的主要构件，这些构件更频繁地得到周围材料的支撑。因此，它们整体上不太像长而细的部件，不易发生侧向弯曲，从而提高了整体屈曲能力。

图17中的屈曲载荷因子与图16中的等边晶格布局一致。如图17所示，随着分配给宏观尺度的体积增加，BLF也增加，宏观主导分割（Vη=1.00, Vx=0.60）获得了最高的值。相对于空心晶格，当偏好宏观尺度连续性时（Vη=1.00, Vx=0.60），以及当分割平衡（0.75, 0.80）时，等边设计的BLF更高；而在微观主导设置中（Vη=0.60, Vx=1.00）时则略低（图17）。综合图15、图16和图17可以看出，两种微观结构都受益于将体积分配给宏观尺度，其中等边晶格在能够形成连续框架时特别有效，而空心晶格在单元格级密度化占主导时仍具有竞争力。

通过将深度神经网络（DNN）整合到双尺度屈曲框架中，我们提高了计算成本，并允许更丰富的微观结构描述。为了说明这种增加的自由度的机械效益，重新解决了具有五个独立微观结构参数的宏观尺度问题，以便每一族等边单元格筋条可以获得不同的厚度（图18）。图18中的优化布局显示了压力加载区域（即图13中的右侧边缘）的最大调整：压缩区域的厚度变得更加均匀，沿加载右侧边缘的筋条显著变厚。如图18（红线）所示，加厚的边界筋条类似于增加了边缘加强件：外部框架更有效地承受弯曲和扭转，并为相邻构件提供了更强的侧向约束，从而提高了它们的稳定性[53]。与这些机制一致，具有五个独立微观结构参数的设计相对于图18中的单参数情况，宏观尺度屈曲载荷因子提高了22%。

下载：下载高分辨率图像（613KB）下载：下载全尺寸图像

图17. 在不同体积约束下，由两种晶格微观结构类型组成的结构的优化屈曲载荷因子（BLF）。

将图15中的粗网格布局与图19中的细网格布局进行比较，细网格使边缘和连接更加平滑且连续，同时保持相同的整体拓扑结构。更平滑的边界减少了角落处的峰值应力，因此单个筋条不太可能在局部发生屈曲。图19中连续的外部构件为内部提供了更好的侧向约束，从而延迟了整体屈曲。细网格中更均匀的筋条厚度也使载荷分布更均匀，因此变形不会集中在一个区域。综合来看，细网格设计略微更硬，对小缺陷不太敏感，并且预测的屈曲载荷高于粗网格设计。

下载：下载高分辨率图像（234KB）下载：下载全尺寸图像

图18. 在第三体积约束下，通过宏观尺度屈曲优化获得的等边三角形晶格微观结构的优化屈曲载荷因子（BLF）值及其对应结构。如红线所示，较厚的筋条向施加载荷的右侧边缘聚集。（有关此图中颜色参考的解释，请参阅本文的网页版本。）

下载：下载高分辨率图像（326KB）下载：下载全尺寸图像

图19. 在不同体积约束下，通过宏观尺度屈曲优化获得的具有更细网格的空心三角形晶格微观结构的优化结构。

6.3. 双尺度屈曲优化

考虑全局和局部稳定性是至关重要的，因为即使整体构件仍然看似稳定，薄晶格筋条也可能在局部发生失效；因此，应在优化中包括微观尺度屈曲。这里讨论了在不同微观结构下考虑两个尺度屈曲的效果。对于空心三角形晶格，在表1中的三种体积分割下，仅优化宏观尺度屈曲与同时优化两个尺度之间的对比非常明显（参见图15中的仅宏观尺度情况和图20中的双尺度结果）。

对于平衡分割（Vη=0.75, Vx=0.80）（其中Vη和Vx分别表示宏观和微观尺度上的体积分数），仅优化宏观尺度的设计包含薄的内部桥接和异质分级，而双尺度解决方案将材料聚集到连续的顶部和底部边界构件中，在压缩区域产生更平滑的密度过渡（参见图15与图20）。更平滑的分级减少了内部的峰值压缩应力，连续的边界构件为附近的构件提供了更强的侧向约束。这种组合提高了局部稳定性，并缩短了内部构件的有效无支撑长度。然而，部分可用体积从全局关键区域重新分配以增强局部稳定性，因此一旦考虑了局部不稳定性，预测的全局容量就会降低。这种机制与图21中报告的平衡分割（Vη=0.75, Vx=0.80）的BLF降低11.6%一致。

如图15和图20所示，在微观主导设置（Vη=0.60, Vx=1.00）中，仅宏观尺度布局中存在的狭窄中央颈部被较短、更厚的筋条和额外的加固所取代。更厚的内部筋条提高了单元格内的局部稳定性，但对微观尺度的重视降低了构件的连续性。连续性的降低意味着侧向约束减弱，从而更显著地降低了整体屈曲裕度。这解释了图21中微观主导设置（Vη=0.60, Vx=1.00）的BLF降低了18.0%。

对于宏观主导设置（Vη=1.00, Vx=0.60），两种方法（即宏观尺度和双尺度屈曲优化）都形成了类似外骨骼的结构；双尺度优化使得外部环变薄，并使剩余筋条的厚度均匀（参见图15与图20）。宏观尺度的连续性在很大程度上得以保持，因此侧向约束仍然很强，内部均匀性有助于避免局部薄弱点。由于整体框架保持连续性，因此在这里考虑局部不稳定性的影响较小，这与图21中报告的宏观主导设置（Vη=1.00, Vx=0.60）的BLF降低12%相符。

综合这些观察结果，可以清楚地看出为什么对于每种分割，双尺度的BLF都低于仅宏观尺度的BLF，以及减少的程度如何随着（Vη, Vx）的变化而变化（其中Vη和Vx分别表示宏观和微观尺度上的体积分数）。当更多体积被引导到微观尺度时，为了增强局部稳定性而转移的材料减少了维持主要构件连续性所需的材料，从而使得整体屈曲下降更多。当宏观尺度连续性占主导时，即使包括微观尺度屈曲，优化后的BLF也会降低，但由于连续的外部框架保持了较强的侧向约束，减少的程度较小；双尺度解决方案主要使内部厚度均匀化，而不破坏宏观尺度的连续性。因此，图20和图21中看到的趋势与这些结构考虑一致。

双尺度优化也针对等边三角形晶格进行了测试。图22中的布局可以与图16中的仅宏观尺度结果进行比较。

下载：下载高分辨率图像（486KB）下载：下载全尺寸图像

图20. (a) 在不同体积约束下，通过耦合的宏观和微观尺度屈曲优化获得的由空心三角形晶格微观结构组成的优化结构（尽管在红色矩形标记的区域实施了Heaviside投影，中间密度仍然存在）；(b) 不同体积约束下的收敛历史图。（有关此图中颜色参考的解释，请参阅本文的网页版本。）

下载：下载高分辨率图像（918KB）下载：下载全尺寸图像

图21. 在不同优化问题下，由具有粗网格的空心三角形晶格微观结构组成的结构的优化屈曲载荷因子（BLF）。

对于平衡分割（Vη=0.75, Vx=0.80），两种方法保留了类似的外部形状，但双尺度解决方案在跨度中部附近引入了更平滑的内部分级和有针对性的加厚（参见图22与图16）。更平稳的过渡减少了局部应力峰值并提高了局部构件的稳定性；然而，一些体积从全局影响力的区域转移到了内部规范化，因此与仅考虑宏观尺度的设计相比，优化的BLF有所降低，这与图23中19.7%的降低值一致。如图16和图22所示，在微主导设置（Vη=0.60，Vx=1.00）下，差异更为明显。仅考虑宏观尺度的设计依赖于大的固体区域，而双尺度结果则在单元格局内重新分配材料，特别是在负载的右侧，从而创建了一个更加连贯且支撑更好的框架。因此，优化的BLF增加了，与图23中报告的11.3%的增益相匹配。对于宏观主导设置（Vη=1.00，Vx=0.60），两种优化都产生了内部几乎均匀的外骨骼（参见图22与图16）。图22中的双尺度解决方案主要使支柱厚度均匀，同时保持了整体框架的完整性。由于全局连续性已经很强，额外的局部调整只会略微改变响应；BLF几乎没有变化，双尺度值低1.9%（见图23）。

下载：下载高分辨率图像（452KB）
下载：下载全尺寸图像

图22. (a) 通过在不同的体积约束下结合宏观尺度和微观尺度屈曲优化得到的等边三角格子微结构组成的优化结构；(b) 不同体积约束下的收敛 históric 图。

下载：下载高分辨率图像（719KB）
下载：下载全尺寸图像

图23. 在不同优化问题中，通过不同体积约束下，由粗网格组成的等边格子微结构结构的优化屈曲载荷因子（BLF）。

下载：下载高分辨率图像（384KB）
下载：下载全尺寸图像

图24. 通过在不同的体积约束和合规性约束下结合宏观尺度和微观尺度屈曲优化得到的两种不同格子微结构类型的优化结构。外骨骼和内骨架区域分别由红色和绿色矩形表示。（关于此图例中颜色的解释，请参阅本文的网页版本。）

6.4. 带有合规性约束的双尺度屈曲优化
尽管主要目标是最大化屈曲抵抗力，但从纯屈曲公式获得的设计可能会导致结构刚度较低。为了控制这一点，我们在双尺度问题中增加了合规性约束。我们首先解决双尺度屈曲问题，评估该解决方案的合规性，然后重新优化，并将上限设定为该基线合规性的70%。对于宏观主导设置（Vη=1.00，Vx=0.60）（其中Vη和Vx分别表示宏观尺度和微观尺度的体积分数），结果配置显示在图24中（参见图20中相应的无约束屈曲解决方案）。

如图24所示，对于空心三角格子，增加合规性约束会将材料从边缘重新分配到内部，形成一种混合的内/外骨骼布局，这与对合规性受限悬臂梁的观察结果一致[18]。这种转变加强了直接压缩下的载荷路径，但降低了外部框架的连续性和截面积的第二矩；因此，优化的屈曲载荷因子从3.34×10^-3降至1.23×10^-3，降低了63%。对于等边三角格子，优化后的结构在视觉上非常接近无约束的双尺度结果（参见图24与图22），屈曲指标仅从4.79×10^-3变为4.77×10^-3。沿外部边界的连续构件同时满足了这两个目标，因为它们提高了屈曲载荷，并提供了大部分的平面刚度。因此，满足合规性限制不需要将材料从边缘移开或改变整体框架；只需进行微调，即将施加压力附近的支柱局部加厚，即可实现更高的刚度。这在实践中非常有用：等边微结构可以在保持基本相同的屈曲性能的同时满足刚度要求，使设计权衡变得简单。图25给出了更广泛的比较，它绘制了在所有体积分割下宏观尺度和双尺度优化的屈曲载荷因子。等边格子始终获得了更高的BLF值。此外，当施加合规性限制以增加刚度时，等边设计基本上保持了相同的屈曲性能，而空心设计则表现出明显的降低（见图25）。这表明等边格子更适合同时实现刚度和屈曲目标，而空心格子在强制刚度时面临更大的权衡。

下载：下载高分辨率图像（614KB）
下载：下载全尺寸图像

图25. 在不同双尺度优化问题中，通过不同体积约束下，由两种格子微结构类型组成的结构的优化屈曲载荷因子（BLF）。

6.5. 权重和合规性约束的参数研究
为了说明所提出的多尺度框架如何处理不同的设计偏好，本节报告了通过改变（i）结合屈曲目标中的宏观-微观权重参数ω和（ii）合规性上限得到的空心三角格子的额外解决方案。结果优化的拓扑结构及其相应的屈曲载荷因子（BLF）在图26和图27中进行了总结。

6.5.1. 宏观-微观权重ω的影响
图26比较了三种代表性的权重设置：宏观和微观BLF权重相等（ω=0.5）、宏观主导权重（ω=0.75）和微观主导权重（ω=0.25）。观察到一个清晰且物理上一致的趋势。当宏观权重增加（ω=0.75）时，优化后的布局更接近于宏观尺度屈曲优化问题：材料集中在主要的全局载荷传递路径上，导致结构类似于宏观尺度优化结果，宏观尺度BLF相对于宏观导向的设计几乎不变。在这种情况下，优化器主要通过加强控制整体不稳定性的主要通道来加强主导的全局机制。相比之下，增加微观贡献（ω=0.25）会产生更厚的构件，并在整个格子中产生更均匀的厚度分布。这表明优化器不仅将材料分配到主导的全局路径上，还分配到多个单元格特征上，以抑制局部微观屈曲模式。结果拓扑结构倾向于不那么“路径集中”，更加均匀稳定，这在局部支柱不稳定性或制造引起的缺陷可能触发过早微观屈曲时通常是所期望的。此外，更宽的厚度分布表明冗余性得到了改善：不需要将材料从边缘移开或改变整体框架；只需进行微调，即通过在施加压力附近局部加厚支柱，就可以实现更高的刚度。这在实践中很有用：等边微结构可以在保持基本相同的屈曲性能的同时满足刚度要求，使得设计权衡变得简单。

图26. 宏观-微观权重ω对优化后的空心三角格子的影响。上图：ω={0.5,0.75,0.25}时的优化结构。下图：相应的BLF。

6.5.2. 合规性上限的影响
图27显示了在（i）无合规性约束，（ii）70%合规性上限（相对于参考优化结构），以及（iii）更严格的60%上限下获得的优化解决方案。随着合规性上限的收紧，拓扑结构朝着增强刚度的特征发展，反映了限制变形的日益增加的需求。在没有合规性约束的情况下，优化器可以优先考虑稳定性，并允许在不会强烈影响目标屈曲目标的区域允许更高的灵活性。当施加合规性约束并将其降低（70%→60%）时，材料会被重新分配以减少全局变形并加强载荷引入区域。特别是对于更严格的限制，设计会发展出一种类似外骨骼的承载结构，该结构连接到施加压力的区域，施加压力区域的支柱变得更厚。这种加强载荷引入和传递区域的做法降低了合规性，但同时改变了稳定性最优的材料分布。观察到一个一致的权衡：随着合规性要求的增强，实现的BLF降低。从机制上看，强制刚度通过迫使额外的材料放置在易变形的区域来消耗设计自由度，这可能会降低可获得的最大屈曲稳定性。这种行为突出了刚度和稳定性在设计空间中的竞争角色：更严格的合规性约束会产生更实用的刚度结构，但通常会降低峰值BLF。

下载：下载高分辨率图像（701KB）
下载：下载全尺寸图像

图27. 合规性上限对优化后的空心三角格子的影响。上图：无合规性约束、70%上限和60%上限时的优化结构。下图：相应的BLF。

6.6. 微结构家族的比较：格子、椭圆和双桁架
本节比较了使用四种微结构家族（空心格子、等边格子、椭圆和双桁架）得到的优化结果。这四种微结构的宏观尺度屈曲优化结果显示在图28中。对于宏观尺度屈曲优化，等边格子提供了最高的BLF，而双桁架在考虑的基准测试中产生了最低的BLF。椭圆微结构的性能接近于基于格子的设计，并且略高于（甚至高于）空心格子（见图29），表明基于平滑嵌入的单元格可以实现与类似桁架的格子相当的全局稳定性。当在双尺度公式中额外强制微观尺度屈曲时，由于增加了对局部不稳定性的约束，通常预期BLF会有所降低。图30显示，当从仅宏观尺度屈曲优化过渡到双尺度屈曲优化时，椭圆微结构的BLF下降幅度较小，这可以视为一种优势：椭圆单元格不太容易受到微观尺度屈曲的影响，因此在考虑了局部不确定性后保持了更大比例的宏观尺度稳定性。这种行为与单元格的几何特性一致：格子微结构包含尖锐的角落和支柱相遇的连接点，这可能会促进局部应力集中并触发局部不稳定性或过早的微观屈曲。相比之下，基于椭圆的单元格引入了更平滑的边界和载荷传递，这减轻了应力集中，并提高了对微观尺度屈曲的鲁棒性，同时保持了较强的全局BLF。

下载：下载高分辨率图像（542KB）
下载：下载全尺寸图像

图28. 在宏观尺度屈曲优化下，四种微结构家族的优化结构。

下载：下载高分辨率图像（431KB）
下载：下载全尺寸图像

图29. 使用不同微结构家族进行宏观尺度屈曲优化的屈曲载荷因子（BLF）比较。

下载：下载高分辨率图像（560KB）
下载：下载全尺寸图像

图30. 空心格子、等边格子和椭圆微结构的双尺度屈曲优化结果和BLF比较。

7. 结论
这项工作的主要贡献是一个实用的、端到端可微的分步流程，用于双尺度屈曲设计，它将（i）宏观尺度屈曲最大化与（ii）最坏情况下的微观尺度屈曲保护结合起来，同时通过DNN代理在计算上是可行的，该代理提供了均匀化的属性和梯度。除了大幅节省计算成本（通过避免每次迭代时的FE均匀化）外，该方法还解锁了对微结构家族和参数化的系统比较，这些参数化否则过于昂贵而无法进行。特别是，结果表明，基于平滑嵌入的单元格（椭圆家族）可以在保持与基于格子的单元格相当的全局屈曲性能的同时，减少对微观尺度屈曲约束的敏感性，表明这是一种有前途的方法，可以实现既全局高效又局部稳健的分层设计。

限制和展望
虽然代理大大加快了当前的研究速度，但扩展到完全三维的架构材料和非常高维的单元格描述将需要更丰富的训练策略和/或也直接近似稳定性测量的代理。下一步自然是要开发基于学习的代理来评估微观尺度屈曲（及其敏感性），以及对高维单元格家族的自适应采样，以进一步扩大在所提出的双尺度稳定性框架内可处理的微架构和加载场景的范围。

cRediT作者贡献声明
Sobhan Honarvar：撰写 – 审查与编辑，撰写 – 原始草稿，可视化，方法论，调查，形式分析，数据管理。
Nolan Black：撰写 – 审查与编辑，方法论，调查，形式分析。
Ahmad R. Najafi：撰写 – 审查与编辑，监督，资源管理，项目管理，方法论，资金获取，概念化。