0 引　言

运输船作为国际贸易中至关重要的运输工具，其结构设计的合理性直接关系到船舶的安全性、经济性和运营效率。高双层底实肋板结构作为大型运输船的关键组成部分，承担着支撑船体、传递载荷以及保障船体结构强度和稳定性的重要任务。在当前船舶设计不断追求轻量化、高性能化的背景下，对高双层底实肋板结构进行优化设计具有显著的现实意义。拓扑优化作为一种先进的结构优化方法，能够在给定的设计空间内，寻找材料的最优分布形式，从而在满足结构性能要求的前提下，最大程度地减轻结构重量，提高材料利用率。因此，将拓扑优化技术应用于大型运输船高双层底实肋板结构设计，具有重要的工程应用价值和理论研究意义。

在众多拓扑优化方法中，变密度法、水平集法和双向渐进法是当前研究和应用较广泛的3种技术。张聪^[1]为提升三体船总体性能并保证舱壁结构强度，应用变密度拓扑优化方法计算优化前后舱壁结构应力分布，结果表明优化能实现舱壁轻量化，还能为加筋布置提供指导。吴贝尼^[2]针对船舶中横剖面结构拓扑优化问题，选取并改进BESO方法，编写可视化应用程序，对油船和挖泥船分别在不同设计阶段开展结构拓扑优化研究。郭德松^[3]以典型船体结构板材开孔为研究对象，通过改进BESO算法得到适用于开孔优化的I - BESO算法，开发多工况优化算法并应用于邮轮开孔优化，检验优化后开孔结构可靠性。蒋垣腾等^[4]基于参数化水平集方法提出新的边界搜索方法，实现静水压力载荷下的水下耐压结构拓扑优化设计，以结构柔度最小化为目标、结构体积为约束研究优化问题，结果表明该方法对复杂边界水下耐压结构拓扑优化有工程应用价值。伍勇等^[5]以V型高速船用柴油机油底壳为研究对象，采用变密度法，确定拓扑优化设计方案，结果显示优化后油底壳质量下降，性能提升，验证了轻量化设计方法可行性，可为柴油机其他部件轻量化提供借鉴。

在当前船舶行业的研究中，优化设计通常侧重于单一方法的应用，例如变密度法或水平集法，本文聚焦于3种方法对比，找出最适合船舶实肋板拓扑优化的方法，构建的标准化对比平台可以为后续开展多方法耦合优化提供方法论基础。针对大型运输船高双层底实肋板结构，本文通过对变密度法、水平集法和双向渐进结构优化法进行系统的对比分析，旨在探讨其在相同设计条件下的表现差异。这一创新性的研究思路不仅能够更加全面地评估每种方法的效果和效率，而且有助于识别在实际工程应用中，各种算法的优缺点与适用范围，从而为船舶结构优化提供更为多样化和科学的选择依据。

1 拓扑优化原理解析 1.1 变密度法

变密度法的基本思路是将优化区域内部的每个单元视为不同的材料密度或空隙。将待优化的零件的材料密度分布建模为相对密度μ，其取值范围为0～1。在这一过程中，目标是将相对密度的分布看作一个函数，通过算法处理该函数，使得单元的密度逐渐趋近于0或1，以便最终生成所需的结构。在优化过程中，引入中间密度可以帮助将原本离散的问题转化为连续的问题。然而，实际中并不存在中间密度，因此我们需要减少这些非0和1的中间密度值，使其更接近这2个极值。常用的材料插值模型有SIMP和RAMP。其中，SIMP模型的原始公式^[6]为：

$ {E}_{\left({x}_{i}\right)}={E}_{\min}+{\left({x}_{i}\right)}^{P}\left({E}_{0}-{E}_{\min}\right), {x}_{i}\in \left[\mathrm{0,1}\right]。$

(1)

式中：$ {x}_{i} $为单元的相对密度值；$ {E}_{\min} $为当$ {x}_{i} $为0时的弹性模量；$ {E}_{0} $为当$ {x}_{i} $为1时的弹性模量；$ P $为惩罚因子，通常取值为2～4。近似公式则可以简化为：

$ {E}_{\left({x}_{i}\right)}={\left({x}_{i}\right)}^{P}{E}_{0}, {x}_{i}\in \left[\mathrm{0,1}\right]。$

(2)

基于变密度理论的SIMP法拓扑优化的数学模型可以表述为最小化结构的柔度函数，公式为：

$ \text{Find}\boldsymbol{x}=\left(x_1,x_2,\cdots,x_n\right)^{\rm{T}}，$

(3)

$ \text{Min}\text{}{C}_{\left(x\right)}={\boldsymbol{U}}^{{T}}\boldsymbol{K}\boldsymbol{U}=\sum _{i=1}^{n}{f}_{i}\left({x}_{i}\right){\boldsymbol{u}}_{i}^{{T}}{\boldsymbol{k}}_{0}{\boldsymbol{u}}_{i}。$

(4)

式中：$ \boldsymbol{x} $为单元相对密度的矢量；$ {C}_{\left(x\right)} $为结构的柔度函数；$ \boldsymbol{U} $为位移矢量；$ \boldsymbol{K} $为结构的刚度矩阵；$ {f}_{i}\left({x}_{i}\right) $为惩罚函数；$ {\boldsymbol{k}}_{0} $和$ {\boldsymbol{u}}_{i} $分别为初始的单元刚度矩阵和单元位移矢量。此外，优化过程中会有约束条件，例如优化后单元的体积不能超过初始体积的一定比例，同时每个单元的相对密度必须在一个合理范围内。这一系列过程的目标是找到最佳的材料配置，以实现轻量化和高性能结构的设计。

1.2 水平集法

水平集法是一种重要的数学方法，主要用于隐式地表示曲线或曲面，其核心思想是利用一个高维的函数来描绘这些形状。基于这一方法，参数化水平集法应运而生。该方法通过引入径向基函数插值，将水平集函数的时间与空间变量进行解耦，表达为函数和扩展系数的乘积形式，即：

$ \varPhi \left(X,t\right)=\boldsymbol{\varPhi }\left(X\right)\boldsymbol{\alpha }\left(t\right)=\sum _{i=1}^{N}{\varPhi }_{i}\left(X\right){\alpha }_{i}\left(t\right)。$

(5)

式中：$ \boldsymbol{\varPhi }\left(X\right) $为径向基函数插值矩阵；$ \boldsymbol{\alpha }\left(t\right) $为扩展系数向量。在这个参数化水平集法的框架下，可以建立一个拓扑优化模型，其目标是最小化结构的柔度，同时将结构的体积分数作为约束条件。具体来说，模型的目标是找到一组扩展系数，使得结构应变能最小化，表达^[7]为：

$ \mathrm{M}\mathrm{i}\mathrm{n}:J\left(\boldsymbol{u},\boldsymbol{\varPhi }\right)={\int }_{\Omega }f\left(\boldsymbol{u},\boldsymbol{u}\right)H\left(\varPhi \right)\mathrm{d}\Omega 。$

(6)

在这一过程中，$ f\left(\boldsymbol{u},\boldsymbol{u}\right) $为应变能的计算。为了满足设计要求，模型还设定了一系列约束条件。在这些约束中，首先是体积约束，即$ G\left(\boldsymbol{\varPhi }\right)={\int }_{\Omega }H\left(\boldsymbol{\varPhi }\right)\mathrm d\Omega \leqslant {V}_{\mathrm{m}\mathrm{a}\mathrm{x}} $，这里的$ {V}_{\mathrm{m}\mathrm{a}\mathrm{x}} $为体积分数的最大允许值。此外，还需要满足扩展系数的上下限，确保这些设计变量在合理的范围内，即$ {\alpha }_{{i},\mathrm{m}\mathrm{i}\mathrm{n}}\leqslant {\alpha }_{i}\leqslant {\alpha }_{{i},\mathrm{m}\mathrm{a}\mathrm{x}} $。通过以上步骤，能够有效进行结构的优化与设计，使其在实际应用中更加高效和准确。

1.3 双向渐进优化法

渐进结构优化（ESO）是一种受自然界生物结构启发的材料优化方法，旨在通过去除设计域内的低材来实现优化。然而这种方法存在一个问题：误删单元后无法复原。为了克服这一缺陷，Huang等^[8]提出了双向渐进结构优化（BESO）算法。该算法的创新之处在于能在每次迭代中同时移除和添加材料，从而将未充分利用的材料重新分配到最需要的位置。这一特性使得BESO算法广泛应用于建筑与桥梁设计，BESO的目标是通过最小化结构的应变能来实现材料的高效使用。目标函数的表达式^[9]为：

$ C=\frac{1}{2}{\int }_{\Omega }{\boldsymbol{\in}}^{\mathrm T}\boldsymbol{D} \in \, {\mathrm{d}}\Omega。$

(7)

式中：$ \boldsymbol{\in} $为应变向量；$ \boldsymbol{D} $为与材料特性相关的矩阵。为了满足设计的体积约束，BESO算法的限制条件设定为：

$ {\int }_{\Omega }{\boldsymbol{v}}_{i}\mathrm d\varOmega -{V}^{\mathrm{*}}=0 。$

(8)

式中：$ {\boldsymbol{v}}_{i} $为节点位移向量；$ {V}^{\mathrm{*}} $为预设的优化后结构的总体积限制值。在BESO算法中，单元的弹性模量与单元密度之间有特定关系，表示为：

$ E\left({\rho }_{i}\right)={E}_{0}{\rho }_{i}^{p}。$

(9)

式中：$ {E}_{\left(0\right)} $为单元的初始弹性模量；$ {\rho }_{i} $为单元密度；$ p $为单元密度的惩罚因子。单元的密度被定义为1表示单元存在，而0表示单元不存在，通常被删除单元的密度取一个较小值（如0.001）。结构的整体刚度是由所有单元的刚度组合而成，计算公式为：

$ \boldsymbol{K}=\sum _{i=1}^{n}E\left({\rho }_{i}\right){\boldsymbol{K}}_{0}^{i}。$

(10)

式中：$ {\boldsymbol{K}}_{0}^{i} $为初始单元刚度矩阵。灵敏度分析是BESO算法的重要环节，灵敏度表达式为：

$ \frac{\partial C}{\partial\rho_i}=-\frac{1}{2}\rho_i^pE_0\boldsymbol{u}^{\mathrm{T}}\boldsymbol{K}_0^i\boldsymbol{u}。$

(11)

这一表达式用于量化每个单元对结构应变能的敏感程度。单元灵敏度有不同情况，当密度为1、最低值或其他状态时，其灵敏度序号的表达式也有所不同。在每次迭代中，BESO算法根据灵敏度调整每个单元的密度。具体而言，高应变能密度的单元会将其密度从最低值调整为1，而低应变能密度的单元则会将其密度从1降低到最低值。同时，算法采用过滤策略以防止“棋盘格现象”，确保优化过程的连续性和合理性。通过不断的迭代，最终实现材料分布的优化，从而达到最佳设计效果。

在工程结构优化领域，变密度法、水平集法与BESO（双向进化结构优化）是3种典型的拓扑优化技术，各自具备独特的适用场景。变密度方法在多工况约束下的复杂工程结构（如机械臂或航空部件）优化中具有显著优势，但需重视灰度单元的后处理环节；水平集法在处理多物理场耦合问题（如热力耦合）和复杂几何设计（如增材制造超材料）方面表现突出，其生成的光滑边界特征使其在形态精度要求较高的领域具有独特价值；BESO方法通过单元迭代删除与添加实现结构布局优化，在增材制造轻量化设计中展现出明显优势。在船舶结构拓扑优化领域，变密度法适用于方案初步探索，水平集法擅长精细化改型，而BESO法则是生产设计阶段的首选工具。

2 大型运输船高双层底肋板结构多方法拓扑优化实例 2.1 大型运输船高双层底肋板结构介绍

本章算例以拓扑优化设计空间内的单元密度为设计变量，以结构剩余体积为约束条件，以结构整体柔度（刚度倒数）最小为目标函数。以74000 t大型运输船船货舱区双层底肋板结构初为优化对象，其双层底高度约为1.8 m，船底设置有2道中纵桁材，除此之外，各有3道旁纵桁材结构，整体结构对称布置。2条中纵桁材间距约为2.5 m；中纵桁材与第一条旁纵桁材间距为2.6 m，各旁纵桁材间距为3.2 m，实肋板结构在两端支撑在舷侧纵舱壁上。所选船舶建造钢材为普通钢材，所选研究结构示意图如图1所示。

首先需要确定拓扑优化的设计区域，通过创建几何模型来描述设计空间，算例使用二维矩形区域，矩形区域长度为12.2 m，高度为1.8 m，与所选优化结构尺寸保持一致。在Matlab中，使用矩阵来表示离散化后的设计空间。对于二维问题，定义一个nelx*nely的矩阵，其中每个元素对应一个有限元单元。为准确评估3种方法的优劣，选择计算成本最高的水平集法来预估网格大小，作为本算例网格尺寸。用水平集法对网格单元大小分别为0.01、0.05、0.08、0.1、0.15 m的悬臂方板进行计算，在其末端加向下的集中力，体积比选为0.5。结果显示，当选用0.05 m网格计算时，优化的精度、效率以及结果的可靠性均达到最佳，故选用该尺寸网格进行优化计算，即nelx为244，nely为36，设计域内总单元数为8784。设计域单元划分示意图见图2。

在船舶工程领域，进行双层底肋板结构设计时，需全面考虑结构设计与施工工艺的基本要求。由于船舶作业工况存在差异，该双层底肋板结构^[10]主要呈现出2种受力状态，且这2种状态具有上下对称的特性。具体来说，受载工况1的情形为，船底板下方作用均布压力载荷，此载荷源于船舶在特定航行条件下与外部介质的相互作用，如海水的压力等。而受载工况2则是内底板上方作用均布压力载荷，这通常与船舶内部货物的分布、液舱压力等因素相关。本研究主要是要探求3种方法的优劣性，不对多工况对优化结果影响进行研究，故选取单一工况采用3种方法进行优化计算，设计边界条件为内底板上方并作用均布压力载荷，在对应纵桁处添加水平与竖直方向固定约束，边界条件如图3所示。

2.2 变密度法优化结果

在变密度法拓扑优化里，容积率、惩罚因子和滤波半径是关键可变参数。容积率^[11]表示设计区域内材料的占有比例，范围通常为0～1，直接关联结构质量与材料利用率，是结构轻量化约束条件。实际取值依设计需求而定，航空航天零部件为实现轻量化，可能设为0.3～0.5；对强度要求高、重量不敏感的建筑结构，可能接近0.9，船舶结构进行轻量化设计时选取容积率为0.3～0.8。惩罚因子^[12]用于解决中间密度问题，其值大于1，旨在修改目标函数与约束条件，让优化算法倾向选择密度为0或1的方案，避免出现物理上不合理的中间密度材料状态。取值越大，对中间密度惩罚越严重。合理选择惩罚因子的取值，可以消除多孔材料，从而得到理想的拓扑优化结果；经过试算，当惩罚因子小于或等于2时，存在大量多孔材料，计算结果没有可制造性；当惩罚因子大于或等于3.5时，最终拓扑结果没有大的改变，当惩罚因子大于或等于4时结构总体柔度的变化非常缓慢，迭代步数增加，计算时间延长，本研究最终选取惩罚因子为3。滤波半径用于对单元密度或灵敏度滤波，避免如棋盘格现象等数值不稳定情况。取值常与有限元模型单元尺寸相关，一般在0.5～4倍单元尺寸之间选择。过小无法有效消除不稳定，过大则会过度平滑，丢失结构细节，查阅资料^[11]选取滤波半径为1.5倍单元尺寸。计算结果如图4所示，每种容积率计算迭代次数如表1所示。

在Matlab中进行拓扑优化时，输出的结果图通常被称为“密度分布图”或“材料分布图”。这些图展示了各单元的材料密度，通常使用灰度或颜色来表示不同区域的材料密度。黑色或深色区域通常表示材料的高密度部分，而白色或浅色区域则表示低密度或无材料的部分。从图4结果可以看出，当外轮廓缩减到某一限度值时（容积率为0.4），开始内部出现细小孔洞，随着容积率继续减小，细小孔洞逐渐拉长，变成长条状，船体板开孔一般为圆孔或者方孔，细小长条状开孔不符合船舶建造规范。在容积率为0.3与0.7时需要更多迭代次数，说明在这2种容积率工况下，优化过程较为缓慢。容积率为0.3时是因为较低的容积率会导致更复杂的结构形状，从而需要更多的迭代来收敛到最优解。由于不同的容积率会导致初始设计的敏感性变化，容积率为0.7时需要更多的迭代。

2.3 水平集法优化结果

在水平集法拓扑优化中，涉及的关键参数的数值选取与变密度法差别很大。水平集法拓扑优化所需实体分数，介于0～1，决定最终设计中实体部分占比是重要设计约束，为与前节变密度法工程上一致，实体分数设为0.3～0.8。在算法的每次迭代中求解演化方程的CFL时间步数也是关键参数，取值与设计域单元行数与列数有关，如果取值过小，则设计变化太慢，并收敛到较差的局部最小值；如果取值过大，设计变化太快，可能会移除重要特征，并收敛到较差的局部最小值；根据前面划分的计算域，计算出时间步数为20。除此之外，还要选择在执行水平集重新初始化之前算法的迭代次数，迭代次数选择过少，重新初始化过于频繁，以至于在设计中无法形成新的孔洞；迭代次数选择过多，水平集函数变得非常陡峭，导致在求解演化方程时精度变差；参考文献[13]，迭代次数选择为2～6的整数，本计算选择迭代次数为3。演化方程中拓扑导数项对函数演化影响很大，取值范围为1～4，值小在设计中无法形成孔洞，因为形状导数主导了设计演化，难形成期望拓扑特征，值大拓扑导数主导演化，从而产生带有锯齿边缘的大孔洞，可能致拓扑结构不理想，如边界锯齿化。试算后，本设计选择权重为3。水平集法拓扑优化计算结果见图5，每种实体分数计算迭代次数见表2。

可以看出，水平集法拓扑优化与变密度法在高实体分数时迭代结果相似，材料都是逐渐向载荷边靠拢，但与变密度法不同的是，在低实体分数时，无明显内部孔洞扩张现象，更符合船体板材设计要求。但水平集法拓扑优化从高实体分数开始，都会出现“材料脱离”现象，即每次都会载载荷边对侧形成点状材料残留，最后结果出现多个边界，丢失材料。当这种现象出现在船体板建造过程中时，纵桁之间的肋板通常设计为完整一块，会使优化结果与成品板材产生误差，应特别注意。当实体分数较小时（如实体分数为0.3），由于结构内部无孔洞，材料都聚集到了载荷边界，导致优化算法在收敛时遇到困难，特别是在接近边界条件时，需要更多的迭代来调整设计。当实体分数为0.4～0.8时，对该结构进行水平集法拓扑优化时，迭代次数相差不大。

2.4 双向渐进结构优化优化结果

在双向进化结构优化（BESO）中，体积分数、进化率和过滤半径是关键参数。体积分数指设计域内材料所占体积比例，介于0～1，直接影响结构的重量、成本与承载能力等性能。在BESO算法里，体积分数作为约束条件，过高设置会使优化结果保守，难以实现显著轻量化；过低则可能使结构无法满足性能要求，致使算法难以收敛。为与前面2种方法实现工程上一致，体积分数设为0.3～0.8。进化率^[14]是控制BESO中材料添加和删除的比例参数，其依据单元灵敏度分析结果，较大的进化率虽能使结构拓扑快速变化，但可能导致优化过程不稳定，甚至错过最优解；较小的进化率会使优化收敛缓慢，增加计算成本。进化率通常为较小正数，其中0.02与0.03使用最多，本研究经过试算后选取进化率为0.02。向进化结构优化方法也需要设置过滤半径，Han等^[15]后选择过滤半径为6。水平集法拓扑优化计算结果如图6所示，每种体积分数计算迭代次数如表3所示。

可以看出，该方法结果也是材料现汇集到一起，然后轮廓逐渐缩小，最后内部出现开孔。但内部孔洞形状多为类半圆与类椭圆，经过修饰后与船舶板材孔洞相似。虽然也有材料分离现象，但只出现在了体积分数取最小的这一算例结果中，其他体积分数结果仍为一体。基于BESO拓扑优化迭代次数随体积分数减小而逐渐增加，当体积分数取得较小时，较低的体积分数会导致更频繁的材料去除和添加，导致迭代次数增加。

2.5 3种优化方法典型结果汇总

选取体积分数为0.3的工况是一个具有代表性的典型案，通过对3种方法在该工况下的优化结果进行局部放大图细节对比（见图7），深入分析每种方法在材料分布和结构形状的具体特征。

通过放大后细节对比可知，变密度法密度分布图中存在灰色过渡单元，这种灰度分布在实际制造中会产生难以确定的结构。水平集法与BESO法没有灰度单元，直接反映出材料的存在与否。这种清晰的材料分布使得制造过程更加简单和直接，更适合于船体板的设计制造。

选取体积分数为0.5、0.6、0.7、0.8工况3种方法变化趋势与迭代次数的过程对比，结果如图8所示。

通过对图中4种体积分数下变化趋势与迭代次数的对比，BESO法在所有体积分数下表现出最快的收敛速度。可以看到，BESO法的迭代次数较少就能达到目标体积分数，显示出其优化过程的高效性。这种快速收敛的特性使得BESO法在实际应用中能够节省计算时间和资源。BESO法的变化趋势非常稳定，迭代过程中几乎没有明显的波动，表明其在优化过程中能够保持一致的性能。

体积分数为0.4时结构差异性最大，更能代表每种方法的结果特点。基于上述的密度分布图，对高双层底实肋板结构进行重建，材料减少了60%，得到经过拓扑优化后的肋板结构如图9所示。对重建后的结构进行有限元分析，施加边界条件与本节的设计工况一致。

3种结构在相同边界条件下应力结果相似，应力较大区域均位于靠近两侧纵桁的肋板上。由图9结果可以看出，变密度法与双向渐进结构优化（BESO）方法结果最大应力位置均位于外侧纵桁与肋板连接底部。重构后变密度法结果最大应力为151.1 MPa，而BESO法结果最大应力为131.0 MPa，相同边界条件下，应力减小了13.3%，所以BESO方法结果强度更好。水平集法结果最大应力为131.3 MPa，位于最外侧肋板中间弧状切口处，结构强度略差于BESO法结果。

3 结　语

本文基于Matlab编程对实肋板模型进行拓扑优化，解决了传统设计中存在的材料冗余问题。对变密度法、水平集法和双向渐进结构优化（BESO）方法的对比研究，发现了这3种算法在优化大型运输船高双层底实肋板结构拓扑优化设计时表现出明显差异。在此基础上，对3种方法形成的结果重构，进行有限元分析。通过对比3种方法的密度分布图，迭代过程，以及重构后的应力结果，得到以下结论：

1）变密度法拓扑优化得到的密度结果图会出现小孔与长条状缝隙，放大后还有大量灰色单元。水平集法在所有结果中都出现了材料脱离现象，最后结果产生了不与结构主体相连的小块，产生了多段边界，在船体板建造过程中，都会使优化结果与成品板材产生较大误差。BESO法结果边界清晰，结构完整，孔洞为类半圆与类椭圆，更符合船舶设计需求。

2）通过对比4种体积分数下的变化趋势与迭代次数，BESO法在所有体积分数中展现出最快的收敛速度。BESO法以较少的迭代次数即可达到目标体积分数，体现了其优化过程的高效性。此外，BESO法的变化趋势稳定，迭代过程中没有明显波动，减少了设计过程中的不确定性。相比之下，变密度法和水平集法在接近目标体积分数时表现出较大波动，导致收敛过程不稳定。

3）保留相同体积的材料，对3种方法的结果进行重构，在相同边界条件下进行有限元分析。结果显示，变密度法强度最低，另外2种方法结果相差不大，BESO法略优于水平集法。