数据驱动的曲面构件形状⁃拓扑协同优化方法

高天贺，田阔，黄蕾，张澍，李增聪

大连理工大学 工业装备结构分析优化与CAE 软件全国重点实验室 工程力学系，大连 116024

曲面构件是组成航空航天结构的关键零部件［1-3］。火箭中的整流罩，飞机机翼、尾喷管以及密封舱壁板、舱门等都是曲面构件。曲面构件在运载火箭中占干重比80%以上，在飞机中占干重比50%以上［4］。随着我国航空航天事业的不断发展，对航空航天结构的寿命及可靠性提出更高的要求［5］，所以研究曲面构件的轻量化设计方法具有重要意义。

近年来，拓扑优化方法凭借其不依赖指定的初始结构的优势，已逐渐成为曲面构件概念设计阶段的重要手段［6］。曲面构件的拓扑优化方法主要包括固体各向同性材料惩罚法（Solid Isotro⁃pic Material Penalized，SIMP）［7］、水 平 集 法（Level Set，LS）［8］、可变形孔洞法（Moving Mor⁃phable Void，MMV）［9］等。其中SIMP 法以单元密度为设计变量，具有收敛性好、灵敏度分析简单等优点，但应用于复杂曲面构件时却容易产生棋盘格或灰度单元等问题，通常需要引入一些过滤技术［10］或B 样条参数化方法［11］。相比于SIMP法，LS 方法用一个高维度的水平集隐式表达拓扑结构边界，产生的孔洞或筋条更加清晰。但是，LS 方法收敛性较差，利用该方法分析非线性问题的灵敏度较为困难。近年来，相关学者将LS 方法与等几何分析相结合［12］，提高了拓扑边界应力的计算精度，解决了应力约束拓扑优化问题应力奇异的挑战。除上述方法外，张卫红团队［13］还提出一种几何背景网格法，用于航天三维曲面薄壁壳的加筋设计。Feng 等［11］针对曲壳结构提出了一种有效的B 样条参数化方法，与基于密度的方法相比，结果更清晰且无棋盘格现象。Zhou 等［14］参考蠕虫的几何形态，提出了基于仿生B-样条偏移的加筋设计方法。Dong 等［15］基于自然界分支系统生长机理，提出了面向薄壁结构抗屈曲的最佳加筋布局自适应生长方法。Hu 等［16］对自适应增长方法进行了改进，针对加筋角度、位置和长度等参数进行了尺寸优化，提高了结构的固有频率。综上，目前学者们已针对曲面构件的拓扑优化方法开展了大量的研究，相比传统依赖人工经验的设计方法，拓扑优化往往更容易产生创新性设计构型。

然而，由于曲面构件的某些性能受曲面形状的影响很大，当曲面构件形状固定时，仅依靠拓扑优化方法难以保证获得的材料分布是最优的［17］。目前，学者们针对曲面构件形状-拓扑协同优化问题开展了一定研究。最早的形状-拓扑协同优化方法，是在拓扑优化后，利用CAD 方法［18］或一些光滑的隐式曲线［19］对孔洞形状进行局部优化。进一步，为将壁面形状与拓扑优化协同考虑，相关学者提出了一种形状与拓扑交替进行，直到优化问题收敛的分步迭代优化方法［20］。该方法在优化过程中会改变曲面形状，为拓扑优化提供了更为灵活的设计域。然而，由于形状和拓扑变量被单独处理，形状与拓扑变量的相互作用无法得到充分的考虑。近年来，学者们开始将SIMP 法［21-22］等拓扑优化方法与形状设计变量的灵敏度分析方法相结合，通过推导灵敏度公式［23］，进行形状与拓扑的耦合求解。然而对于复杂问题，灵敏度的推导较为困难，将耦合求解方法应用于实际工程还存在较大难度。如何在考虑形状-拓扑相互作用的同时，避免复杂的设计变量灵敏度推导，使其能较好的应用于工程实际中的复杂问题，是形状-拓扑协同优化的研究挑战。

近年来，数据驱动方法受到国内外学者的广泛关注，为解决形状-拓扑的协同优化问题提供了新的途径。目前已有许多新颖的数据驱动的形状优化［24］和拓扑优化［25］方法，然而如何应用数据驱动思想，将形状优化与拓扑优化协同考虑，还有待进一步研究。针对数据驱动的建模，相关学者先后提出了网格映射［26］和网格变形［27］，并成功应用于回转曲面的拓扑优化模型重构［28］。该方法为参数化建立有限元网格模型，实现形状-拓扑协同优化提供了可能。数据驱动方法可以利用一定数目的样本，拟合具有较高非线性程度的代理模型，并利用代理模型代替真实计算，能有效提高优化求解的效率［29-31］。Meng 等［32］构 建 了Kriging 自 适 应 代 理模型，并开展了多学科协同优化设计，显著提高了高度非线性系统协同优化的精度和效率。因此，基于数据驱动的网格变形建模及代理模型优化方法，提出一种数据驱动的曲面构件形状-拓扑协同优化方法，旨在拓宽航空航天曲面构件设计空间，为航空航天曲面构件提供更为轻量化的创新构型。

本文主要包括3 部分，第1 部分介绍数据驱动的曲面构件形状-拓扑协同优化方法；第2 部分开展简支梁和航天器舱门算例研究，以验证所提出方法的有效性；第3 部分为结论与展望。

1 数据驱动的形状⁃拓扑协同优化方法

首先介绍形状方程-网格变形驱动的形状优化方法，该方法通过用数学方程中的少量参数描述控制点的位移，将形状设计变量用少量的参数表示。然后，介绍考虑冲压约束的SIMP 拓扑优化方法，针对每一个网格变形的设计域开展拓扑优化设计，获得优化结果。最后，介绍形状-拓扑协同优化的技术路线，包括离线、在线、更新3 个阶段。

1.1 形状方程⁃网格变形驱动的形状优化方法

传统的基于节点的形状优化方法往往面临着设计变量规模较大、网格更新困难和容易收敛到不满足工艺要求的构型设计等问题［33］。因此，提出一种以形状方程描述形状设计变量、以网格变形为驱动的形状优化方法。该方法能够减少形状设计变量个数且避免形状改变后网格重新剖分，便于获得满足真实需求的设计构型。如图1 所示，该方法主要包括以下步骤：

图1 形状方程-网格变形驱动的形状优化方法Fig.1 Shape equation-mesh deformation-driven shape optimization method

步骤1：建立有限元模型，划分形状优化的可设计域和不可设计域。在剖分好的网格模型中均匀选取控制点，以保证网格变形后网格边界的光滑性。

步骤2：建立描述壁面或截面形状的数学方程，根据数学方程求解控制点网格变形后的位置坐标。

步骤3：基于径向基函数神经网络，拟合控制点变形前后的两个坐标数据集合，获得二者之间的映射关系。

步骤4：基于训练好的映射关系，对原网格模型所有节点坐标进行预测，生成新的网格模型。为保证网格质量，对变形后的网格进行拉普拉斯网格光顺化处理，获得边界平滑、网格质量较高的有限元网格模型［27］。

经过上述步骤，即可完成网格变形驱动的参数化建模过程。需要说明的是，上述方法将数学方程描述壁面形状和网格变形改变壁面形状相结合，仅需要提供数学方程的若干参数，便可自动的进行控制点的移动、映射关系的训练，生成可用于有限元分析的网格模型。

1.2 固体各向同性材料惩罚拓扑优化方法

SIMP 法［34］作为一种经典的拓扑优化方法，具有收敛性好、灵敏度分析简单等优点。因此，以SIMP 法为例，开展形状-拓扑协同优化的研究。以结构的应变能最小为优化目标，体积为约束，优化问题为

式中：ρi为单元i的伪密度值；N为模型设计域的单元总数；为避免刚度阵奇异，δ取0.001；F为各单元的等效节点力；U为各单元节点位移；K为整体刚度矩阵；V0为结构初始体积；Vf为约束体分比。

曲面加筋是一种常见的结构形式，为获得符合工艺约束的拓扑结果，往往需要约束冲压方向每一列材料的伪密度相等，来施加冲压约束［35］。以图2 为例，假设冲压方向共四层网格，则其冲压约束表达式为

图2 冲压约束效果示意图Fig.2 Schematic diagram of effect of stamping constraint

冲压方向网格单元伪密度在优化过程中保持相等。

通过求解上述优化问题，获得优化后的伪密度场，即可获得结构的拓扑优化构型。

1.3 数据驱动的形状⁃拓扑协同优化方法

一种数据驱动的形状-拓扑协同优化方法如图3 所示，所提出方法主要包括以下3 个阶段：

图3 数据驱动的形状-拓扑协同优化方法流程图Fig.3 Flow chart of data-driven shape-topology collab⁃orative optimization method

1） 离线阶段。离线阶段需要利用少量的样本点建立全局高精度的代理模型。首先，基于拉丁超立方采样方法［36］在设计空间进行采样获得形状变量样本（即形状方程中的参数组合），使样本点具有较好的空间覆盖性，避免构建的代理模型局部坍塌［37-40］。然后，将形状变量样本代入形状方程，利用网格变形技术参数化、自动化地获得样本点数据对应形状的网格模型。基于SIMP法对获得的网格模型进行拓扑优化，获得拓扑优化后的响应结果（如应变能、最大位移、频率等）。训练径向基函数（Radial Basis Function，RBF）代理模型［41］，获得形状设计变量（数学方程中的若干参数）与拓扑优化响应结果之间的映射关系。其中，代理模型的精度指标［1］为留一验证的R2和RRMSE，其表达式分别为

式中：k为所有样本点的个数；yi为第i个样本点响应的真实值；ŷi为第i个样本点响应的预测值；yˉ为所有样本响应的平均值。R2越接近1，RRMSE越接近0，代理模型的预测精度越高。

在实际应用中，可根据实际问题的难易程度及计算资源确定合适的训练样本数目。需要指出的是，当边界和载荷条件发生变化时，代理模型需要重新训练。为保证在线阶段优化能够稳定快速收敛，离线阶段建立的代理模型需要保证具有较高的全局精度，如R2大于0.9。如果精度过低，可以通过自适应序列加点策略（如CVVoronoi 方 法［38］、EI 加 点 方 法［39］等）进 行 更 新 和逐步提升精度。

2） 在线阶段。基于全局寻优算法协方差矩阵自适应进化策略（Covariance Matrix Adaptive Evolution Strategy，CMA-ES）［42］和RBF 代理模型进行优化。该阶段基于离线阶段训练好的代理模型开展高效的优化，无需进行拓扑优化设计。由于代理模型根据输入获得输出几乎无耗时，所以可在优化过程中使用充分多的种群数目和迭代次数。CMA-ES 优化结束后，以CMAES 优化得到的最优结果为起点，利用序列二次规划 算 法（Sequential Quadratic Programming，SQP）［43］对RBF 进 行 再 次 局 部 寻 优，确 保 真 正收敛。

3） 更新阶段。为避免程序无限制的运行下去导致RBF 代理模型在局部密集加点，甚至产生矩阵奇异等问题。根据韩忠华等［44-45］以及Fang等［46］的研究成果，利用最大迭代次数限制更新次数。最大迭代次数可根据问题难易程度或代理模型精度设置。所提出方法利用形状方程驱动的网格变形技术降低了设计变量维度和优化难度，并通过在离线阶段构建具有较高精度的代理模型、在线阶段通过CMA-ES 全局搜索结合SQP 局部搜索的模式，保证优化算法的稳定收敛。因此通常不需要设置过大的迭代次数就足够保证优化收敛。在线阶段结束后，判断是否达到最大迭代次数，如果未达到最大迭代次数，则将在线阶段得到的优化解真实计算后加入原数据集，并重新训练代理模型以达到提高代理模型精度的目的。代理模型更新结束后，基于重新训练的代理模型重新进行在线阶段的优化，直至达到最大迭代次数。更新阶段可通过加点更新提升代理模型的局部精度，自适应的保证优化结果的稳定，减少拉丁超立方采样随机性对优化结果的影响。

核心创新点为提出了一种数据驱动的形状-拓扑协同优化框架，通过建立形状变量与拓扑优化响应的代理模型实现形状与拓扑的协同优化。其中，离线阶段的代理模型的采样方法以及在线计算的优化算法可以根据需要进行替换，并不影响整个优化框架。

2 算例研究

2.1 算例一：简支梁算例

首先，基于一个简支梁算例来验证所提出方法的有效性。简支梁的尺寸和载荷边界如图4所示，考虑载荷及边界的对称性，对梁进行二分之一简化。取梁的几何中心为坐标原点，x轴为梁的跨度方向。梁上端面承受f=0.1 MPa 均布载荷，采用六面体单元对半跨简支梁进行网格划分，单元数目为1 200 个。简支梁的弹性模量和泊松比分别为E= 67 GPa 和υ= 0.3，材料密度为ρ= 2.700 g/cm3。

图4 均布压强载荷作用下的简支梁示意图Fig.4 Schematic diagram of a simply supported beam under sinusoidal pressure load

梁的形状变化如图5 所示，要求梁沿跨中截面保持左右对称，半跨梁的外轮廓必须为直线。梁跨中截面可以上下移动，记移动量为a，沿y轴方向为正，梁跨中截面高度可以增加或减小，记改变量为2b，高度增加为正。在半跨梁的上下端面均匀选取控制点，控制点坐标可由a和b表示：

图5 简支梁形状优化示意图Fig.5 Schematic diagram of shape optimization for the simply supported beam

式中：坐标原点位于变形前梁跨中截面的几何中心；yU和yL分别为梁上下端面控制点的纵坐标。

形状-拓扑优化以控制形状的两个参数a和b为设计变量，以拓扑优化后应变能最小为目标，约束拓扑优化后模型质量小于1.4 g。拓扑优化以应变能最小为目标，约束体分比小于0.4。

离线阶段中，利用拉丁超立方采样方法在形状设计空间内抽取15 个样本点，并利用网格变形技术和拓扑优化方法进行建模和计算，获得样本点的应变能和质量数据，并利用RBF 代理模型拟合设计变量与应变能、质量的映射关系。其中应变能的预测精度R2为0.956 5，质量的预测精度R2大于0.999 9。使用RBF 的原因是其训练收敛速度快，耗时少，且具有较强的逼近能力［47］。根据测试对比，利用高斯过程代理模型预测应变能的精度R2为0.949 9，深度神经网络的预测精度R2为0.818 4。进行一次留一验证耗时分别为：RBF 耗时0.001 114 s；高斯过程耗时2.007 636 s；深度神经网络耗时2 216.64 s。可以发现RBF 的训练速度和精度均较高。然后进入在线阶段的优化，CMA-ES 最大迭代次数设置为100，CMA-ES 结束后进行SQP 的继续寻优，最大迭代次数为20。将最终优化结果真实计算后进入更新阶段并更新代理模型，代理模型更新次数设为5 次。优化迭代过程及优化目标（Objective，Obj）的迭代曲线如图6 所示，随着代理模型的更新迭代，形状-拓扑协同优化的目标函数逐渐减小并最终收敛，代理模型优化解的形状和拓扑结果也趋于收敛，最终获得优化后的形状与孔洞构型。

图6 简支梁形状-拓扑优化目标函数迭代过程Fig.6 Iterative process of objective function in shapetopology optimization of simply supported beam

进一步，将所提出方法与初始设计（仅拓扑优化）的优化结果进行对比。为保证公平对比在有限计算资源下不同优化方法的寻优能力，本文将计算时间作为控制变量，通过控制样本点数目一致保证各方法计算时间一致。通过对比相同计算耗时下的优化结果，验证所提出方法在有限计算资源下具有较高的寻优能力。对比结果如表1 所示，优化得到的形状及拓扑结果如图7 及图8 所示，其中所提出方法（形状-拓扑协同优化设计Ⅰ）优化得到的应变能为0.019 9 mJ（1 mJ=10−3J），与初始设计相比应变能降低20.08%，说明设计域形状的改变为结构设计带来更大的优化空间。此外，也与形状-拓扑协同优化设计Ⅱ（19 代理模型样本+1 精细校核）和Ⅲ（不构建代理模型，直接CMA-ES 优化）进行了对比，结果显示所提出方法在相同计算资源下可获得力学性能更优的设计。其中精细校核与优化加点的区别在于，优化加点需要将真实计算后的优化解加入训练集并更新代理模型，而精细校核只需要真实计算优化解，不需要更新代理模型。由于达到相同的寻优效果所提出方法用时最少，也验证了所提出方法的计算效率较高。此外，本节也尝试定量对比无代理模型方法（形状-拓扑协同优化设计Ⅴ）与所提出方法的计算效率。然而，无代理模型方法迭代100次（所提出方法耗时的5倍）仍与形状-拓扑优化Ⅰ的优化结果有9%的差距。说明所提出方法的优化效率是无代理模型的5 倍以上。形状设计变量上下限优化结果如图7 及表1 所示。此外，还测试了设计域增大后相同材料用量约束的单一拓扑优化，其中扩大后设计域跨中宽度为120 mm，足够包络所提出方法设计空间的所有形状。优化结果如图7 （c）所示，优化构型不清晰，优化后应变能为0.040 7 mJ，相比所提出方法优化结果增大了1 倍以上。说明拓扑优化设计域并不是越大越好，直接扩大拓扑优化设计域并不总会得到与形状-拓扑协同优化类似的结果，进一步验证了所提出方法的必要性。

表1 所提出方法与其他方法简支梁优化结果对比Table 1 Comparison between results of simply supported beam from proposed method and other methods

图7 简支梁上下限形状及扩大设计域后拓扑结果对比Fig.7 Comparison of topology optimization results be⁃tween upper and lower bounds shapes and ex⁃tended design domain of simply supported beams

图8 所提出方法与其他方法优化后简支梁形状和拓扑结果Fig.8 Shape and topology results of simply supported beam optimized by the proposed method and other methods

最后，本节也尝试测试不同采样方法对优化结果及空间覆盖性的影响。空间覆盖性利用代理模型在验证集上的精度表示。因此，首先利用全因子（Full Factorial， FF）试验设计［48］在设计空间均匀抽取100 个样本点作为验证集。然后，利用FF法均匀抽取15 个样本点进行代理模型构建，最终优化结果为形状-拓扑优化Ⅳ。经验证，FF 法样本点构建的代理模型在验证集上的R2为0.974 6，LHS 的R2为0.993 9，说明FF 抽样法的空间覆盖度较差。在该问题中，两种采样方法的优化结果应变能误差为1.5%，说明算法较为稳定，改变采样方法，并不会对最终优化结果造成较大的影响。

2.2 算例二：航天器舱门算例

为进一步验证所提出方法的有效性，本节针对航天器舱门开展算例研究。舱门的简化模型及截面尺寸如图9 所示，舱门门体安装在门框上，舱门半径为560 mm，包括蒙皮在内的设计域厚度为30 mm，门体中心半径100 mm 的圆面在优化过程中保持水平。门体外表面承受p= 0.09 MPa 均布压强作用，门框底部设置固支边界。根据对称性取四分之一模型进行简化，对简化后模型截面设置对称边界条件。采用六面体单元对舱门结构进行网格划分，单元数目为57 780 个。舱门结构材料属性如表2 所示。

表2 舱门结构材料属性Table 2 Material properties of cabin door

图9 舱门载荷边界及形状变化示意图Fig.9 Schematic diagram of cabin door load boundary and shape change

舱门门体截面控制点y坐标可用超椭圆方程表示，根据模型的尺寸，可以构造出门体截面的数学方程：式中：c和n为舱门算例的两个形状设计变量；c为超椭圆的短半轴长度，可以表示舱门截面边界的内凹程度，c越大，舱门截面边界内凹程度越大，当c接近0 时，舱门截面边界形状趋于水平直线；n为超椭圆的次数，n为1 时舱门截面边界为直线，n为2 时为标准椭圆。通过改变c和n，可实现舱门截面形状的改变。

形状-拓扑优化通过优化形状设计变量，使拓扑优化后的应变能最小。拓扑优化以应变能最小为目标，以优化后的体积小于等于2 022.47 cm3（初始设计形状设计域体积的30%）为约束，保证优化方法最终构型体积相同。为保证拓扑优化结果满足工艺制造要求，对拓扑优化问题施加冲压约束，冲压方向为舱门门体的面外方向。此外，施加尺寸约束，最小尺寸为20 mm，最大尺寸为60 mm，最小间隙设为200 mm。

离线阶段中，初始样本点个数为30 个，RBF代 理 模 型 预 测 精 度R2为0.931 0，RRMSE 为0.674 3。在线优化CMA-ES 设置最大迭代次数为100，SQP 最大迭代次数为20。代理模型更新次数设置为15 次，优化迭代过程及优化目标（Obj）的迭代曲线如图10 所示，最终可获得清晰的设计构型。

图10 舱门形状-拓扑协同优化目标函数迭代过程Fig.10 Iterative process of objective function in shapetopology optimization of cabin door

形状设计变量上下限优化结果如图11 及表3 所示，上下限形状优化结果分别为15 040 mJ和14 038 mJ。

表3 所提出方法与其他方法舱门优化结果对比Table 3 Comparison between results of cabin door from proposed method and other methods

图11 舱门上下限形状及扩大设计域后拓扑结果对比Fig.11 Comparison of topology optimization results be⁃tween upper and lower bounds shapes and ex⁃tended design domain of cabin doors

补充设计域足够包络住形状设计空间的单一拓扑优化结果，如图11（c）所示。最终优化结果如表3 所示，所提出方法优化结果（形状-拓扑协同优化设计Ⅰ）为7 437.38 mJ，对比相同质量下初始设计（仅拓扑优化），应变能结果降低了37.93%，验证了形状-拓扑协同优化的有效性。扩大后的设计域上端面水平，下端面斜率为1，门体中心厚度为120 mm，优化后应变能为29 981.2 mJ。由于只进行拓扑优化，忽略了蒙皮形状对拓扑优化结果的影响，该优化结果相比所提出方法优化结果增大了4倍以上。同时，将形状-拓扑协同优化设计Ⅰ与形状-拓扑协同优化设计Ⅱ（无代理模型更新）、Ⅲ（不构建代理模型，直接优化）进行对比，结果显示所提出方法在相同计算资源下可获得力学性能更优的设计，与算例一结论完全相同。由图12 可知，所提出方法优化结果形状规则且筋条清晰，由于冲压约束在拓扑优化中的使用，使拓扑优化产生的筋条与蒙皮保持垂直，符合工艺要求。

图12 所提出方法与其他方法舱门优化后形状和拓扑结果Fig.12 Shape and topology results of cabin door optimized by the proposed method and other methods

3 结 论

通过网格变形技术和代理模型方法，实现了壁面形状和拓扑的协同优化。最终形成了离线-在线-更新三阶段的优化路线。面向工程中简支梁和航天器舱门结构开展了算例研究，结果表明，相比于固定形状的拓扑优化结果，形状-拓扑协同优化的应变能可分别降低20.08%和37.93%，验证了所提出方法的有效性。此外，还将所提出方法与无代理模型的形状-拓扑优化、无代理模型更新阶段的形状-拓扑优化进行了对比，结果表明，在相同计算资源下，相比无代理模型方法，所提出方法优化的应变能分别减小18.11%和4.87%，验证了代理模型可提高相同计算资源下的寻优能力；相比无代理模型更新阶段的方法，所提出方法的应变能分别减小3.86%和1.09%，验证了代理模型更新阶段可保证最终优化收敛。

在后续研究中，将探索设计域内包含不可变结构功能性孔洞的网格变形方法，实现更为智能的参数化建模。此外，还将研究代理模型预测精度对寻优能力的影响，尝试利用代理模型精度自适应平衡采样数目与代理模型更新次数，并扩展到更为复杂的载荷工况和结构形式，进一步提高方法对复杂问题的适用性。