潜水器是人类从事海洋开发的重要装备之一，随着海洋开发步伐的加快，提高航速、安全性和设计效率，降低艇体重量等日益受到关注和重视。然而，潜水器的设计过程涉及的学科众多，且需要进行大量的协调和反复修改。传统潜水器设计基于串行模式进行，缺乏有效的解耦及协调手段，存在计算花费大、设计周期长、最优解难以获得等不足。多学科设计优化在航空领域应用^[1]，其在复杂系统的设计过程中集成了各学科知识，充分考虑了学科间的耦合关系，以有效的设计/搜索策略为基础，同时针对复杂系统的整个设计过程，基于具有分布式特性的计算机网络系统，进行有效的组织、管理，并将学科之间的相互作用所产生的协同效应进行有效利用，最终目标是获得设计对象的整体最优解^[2]。

多学科设计优化适用于复杂系统设计的优化求解，已在众多领域进行了成功应用。在潜水器设计领域将其应用多学科设计优化得到很多应用。文献[3-5]应用到AUV设计; 文献[6-7]应用到鱼雷设计；应用到载人潜水器设计；文献[12]应用到深海空间站设计；文献[13]应用到微小型水下航行器设计等。以上研究对于提高潜水器的总体性能、促进其技术进步等起到了重要的作用。准确的学科分析模型是实现潜水器多学科设计优化的关键工作之一，对于其性能的提升起着至关重要的作用。然而，以上研究的学科建模大多基于公式，尽管这种方法计算花费低、建模难度小，但计算精度低。随着技术的发展，近似模型、参数化建模等方法已引入到学科建模中。然而，近似模型的样本点数量与其拟合精度很难平衡，而直接将参数化分析应用于学科建模将导致整个求解过程变得十分复杂。在这种情况下，利用近似模型和参数化分析的优点，建立基于参数化分析和近似模型的分析模型，可有效提高潜水器性能与设计效率。

本研究针对一型潜水器^[5]的主体结构设计，设计了基于参数化的无人水下航行器外形阻力、耐压壳体结构分析流程，并建立了近似模型，进一步基于上述近似模型建立了无人潜水器主体结构优化模型，进行了优化的求解，选取了其中的方案进行了对比。

1 外形阻力分析模型 1.1 阻力求解方法

利用RANS法求解外形阻力，其在时均化的N-S方程中，以某种模型的方式将瞬态脉动量体现出来，雷诺时均化后的不可压连续性方程和RANS方程：

$ \partial \left( {{{\bar u}_i}} \right)/\partial {x_i} = 0 $

(1)

$ \begin{array}{l} \partial \left( {{\rho _i}{{\bar u}_i}} \right)/\partial t + {\rho _j}{{\bar u}_j}({\partial _i}{{\bar u}_i}/\partial {x_j}) = \\ \rho {{\bar F}_i} - \frac{{\partial \bar p}}{{\partial {x_i}}} + \frac{\partial }{{\partial {x_j}}}\left( {\mu \frac{{\partial {{\bar u}_i}}}{{\partial {x_j}}} - \rho \overline {u_i^\prime } \overline {u_j^\prime } } \right) \end{array} $

(2)

式中：u_i、u_j为速度分量时均值，m/s；p为压强时均值，Pa；R_i^N=(λ₀+Δλ_i)P^N-I^N为流体密度，kg/m³；K^NM{v_i^M；c_i^M}={P^N；R_i^N}为动力粘性系数，m²/s；$\tilde v_i^N$雷诺应力项。

由于$\tilde c_i^N$的存在而导致方程不封闭，需要对其做出某种假定。在此引入湍流模型SST k-w模型。k、w的输运方程为^[14]：

$ \frac{{\partial \left( {\rho k} \right)}}{{\partial t}} + \frac{{\partial \left( {\rho k{u_i}} \right)}}{{\partial {x_i}}} = \frac{\partial }{{\partial {x_j}}}\left( {{\mathit{\Gamma }_k}\frac{{\partial k}}{{\partial {x_j}}}} \right) + {G_k} - {Y_k} $

(3)

$ \frac{{\partial \left( {\rho \omega } \right)}}{{\partial t}} + \frac{{\partial \left( {\rho \omega {u_i}} \right)}}{{\partial {x_i}}} = \frac{\partial }{{\partial {x_j}}}\left( {{\mathit{\Gamma }_\omega }\frac{{\partial \omega }}{{\partial {x_j}}}} \right){\rm{ + }}{G_\omega } - {Y_\omega } + {D_\omega } $

(4)

式中：G_k和G_w为由于平均速度梯度引起的湍流动能；Y_k和Y_w为关于k和w的湍流耗散项。

1.2 基于STAR-CCM+的外形阻力分析

利用STAR-CCM+进行阻力计算，通过设计参数的设置，可实现几何模型的参数化。潜水器外形见图 1。

图 1中，主要参数包括总长度x₁、半球形封头半径x₂、圆柱壳体长度x₃和艉部长度x₄。外形为对称模型，以中纵剖面为界计算一半模型在水中运动的阻力^[15]。将原点设置于模型的艉部，即艉部断电为原点，计算区域为长方体，尺度为前后均为20 m，上、下、侧面均为5 m，各层的网格均为加密网格，模型设置为水平，采用SSTρ_i湍流模型，强度0.01，粘度比10。模型设置如图 2。

计算模型水深400 m、速度30 m/s条件下的航行阻力为6 873.02 N，则整体阻力为13 746.04 N。同文献[4]对比两者相比误差为5.2%，说明计算结果具有一定准确性。

1.3 外形阻力参数化分析流程

利用STAR-CCM+软件Java宏录制的功能，对Java宏文件进行二次编译。在STAR-CCM+软件中完成所有设置后开始录制宏文件，记录修改需求参数的命令。为便于修改，需要进行外形的参数、计算的时间步长(隐式定长的除外)、内部迭代的次数、计算所需时间等的记录。概括而言，主要包括模型尺度和计算时长的修改、计算结果的输出等。Java宏文件准备完毕后通过*.bat批处理文件直接运行Java宏文件，系统直接调用cmd.exe并运行批处理文件中的各个命令^[16]。参数化的无人潜水器主体结构外形分析流程见图 3。

通过编写*.bat文件调用STAR-CCM+，将前面所述的初始设置以*.sim格式保留在文件中，为了方便Isight软件对于程序进行读取和编辑，Java宏文件进行外形参数设置的保留。分析过程所设计的控制参数、计算结果等借助于Isight软件进行后处理^[17]，将计算所生成的结果文件存放于相应的工作目录下。得到81组样本点列于下表。

2 耐压结构分析模型 2.1 耐压结构非线性分析方法

采用Abaqus有限元软件进行耐压结构的分析，包括最大应力σ_max和极限载荷P_cr的求解。屈曲分析包括线性和非线性2部分，线性分析首先进行，随后采用弧长法(需要考虑和引入初挠度)，即Riks法进行非线性的分析，在此部分中非线性的影响还需要考虑材料和几何，计算得到的极限载荷在载荷比例因子变化曲线中的位置处于最高处，采用Riks法与Newton-Raphson法结合求解非线性屈曲问题。具体过程为^[18]：

$ \Delta {{\rm{ \mathsf{ λ} }}_i} = \Delta {{\rm{ \mathsf{ λ} }}_0}, \Delta u_i^N = \Delta {{\rm{ \mathsf{ λ} }}_0}v_0^N $

(1)

对于第i增量步(i=1, 2, 3, …)：

1)内部结点应力为：

$ {I^N} = \int_V {{\beta ^N}} :\sigma {\rm{d}}V;{K^{NM}} = \partial {I^N}/\partial {u^M} $

(2)

2)检查平衡方程：

$ R_i^N = \left( {{{\rm{ \mathsf{ λ} }}_0} + \Delta {{\rm{ \mathsf{ λ} }}_i}} \right){P^N} - {I^N} $

(3)

若各项值非常小，说明结果收敛，如果不收敛，则需求解方程

$ {K^{NM}}\{ v_i^M;c_i^M\} = \{ {P^N};R_i^N\} $

(4)

2个方向的分量是P^N、R^N，2个方向的位移分量是v_i^N、c_i^N。

3)将矢量($\tilde v_i^N$; 1)缩放，代入到($\tilde c_i^N$; ρ_i)(其中ρ_i=R_i^NP^N/$\overline {{P^2}} $)，得到：

$ \left\{ {\left( {0; - {\rho _i}} \right) + {\rm{ }}\left( {\tilde c_i^N;{\rho _i}} \right) + \mu \left( {\tilde v_i^N;1} \right)} \right\}:\left( {\tilde v_0^N;1} \right) = 0 $

(5)

式中$\mu = - \tilde c_i^N\tilde v_o^N/\left( {\tilde v_i^N\tilde v_o^N + 1} \right)$。则Aⁱ的解为：(u₀^N+Δu_i^N+c^N_i+μv_i^N; λ₀+Δλ_i+μ)。

4)下一次迭代：

$ \left\{ \begin{array}{l} \Delta u_{i + 1}^N = \Delta u_i^N + c_i^N + \mu v_i^N\\ \Delta {{\rm{ \mathsf{ λ} }}_{i + 1}} = \Delta {{\rm{ \mathsf{ λ} }}_i} + \mu ;i = i + 1 \end{array} \right. $

(6)

下一次迭代则返回式(1)时，式(2)中需要做出额外修正，此时Abaqus软件会在这个阶段后，令v_o^N=v_i^N。

2.2 耐压壳结构有限元分析

耐压壳体选用圆柱和球壳组合形式，采用Abaqus中的的Shell壳单元建立模型，宽度、高度、艇长方向分别为X轴、Y轴、Z轴，上述3个坐标轴方向依次设置为向左为正、向上为正、艇长，具体如图 4。

假设潜水器在深水中处于静水压力为均匀的条件下，其所承受的载荷是关于3个平面是对称的，只需对关键节点处的线位移进行约束，而不需使用角位移约束，边界条件施加情况见下表。采用25 mm×25 mm的网格尺寸，网格划分见图 5。

2.3 耐压壳体参数化分析流程

耐压壳体参数化分析流程与图 3所的外形参数分析流程类似，具体不同在于：1)设计变量为耐压壳体主体及环肋的尺寸及相应的控制参数；2)Abaqus的脚本编程基于Python语言，建模、属性设置、载荷与边界条件设置、网格划分、分析和后处理操作等基于Python语言。

为实现耐压壳的参数化建模，利用Abaqus软件进行耐压壳的参数化建模，需要对于模型的相关参数在Python文件中进行修改。这就需要借助Abaqus进程外的批处理*.bat文件来完成这个过程，最终实现Abaqus有限元分析软件能在的后台完成分析。需要首先进行耐压壳单元最大应力和模态特征值的计算，然后计算载荷因子曲线，经过后处理得到耐压壳的极限载荷。计算得到81组样本点值见表 3。

3 耐压壳体分析模型的建立 3.1 响应面模型

响应面模型(response surface method，RSM)基本原理是通过变量和响应之间所存在多项式关系，对于两者通过简单函数关系进行表达，具有透明性高、简单、高效的特点。响应面模型按照阶数可分为一阶~四阶。其中，二阶响应面模型可表达为:

$ \hat y\left( x \right) = {a_0} + \sum\limits_{i = 1}^n {{b_i}{x_i}} + \sum\limits_{i = 1}^n {{c_{ii}}x_i^2} + \sum\limits_{ \le i \ne j \le n}^n {{d_{ij}}{x_i}{x_j}} $

(7)

式中:ŷ(x)为响应面近似值；x_i为设计变量；n为设计变量个数；a₀、b_i、c_ii、d_ij为待定系数。

响应面模型的拟合精度的预测可借助复相关系数R²进行预测，R²的取值范围为[0, 1]，R²越接近1则近似模型拟合精度越高。R²可得到：

$ {R^2} = 1 - \left( {\sum\limits_{i = 1}^n {{{\left( {{y_i} - {{\hat y}_i}} \right)}^2}} /\sum\limits_{i = 1}^n {{{\left( {{y_i} - {{\bar y}_i}} \right)}^2}} } \right) $

(8)

式中：y_i为第i个样本对应的状态变量的响应值；ŷ_i为第i个样本对应的状态变量的近似值；${{\bar y}_i}$为第i个样本对应的状态变量的响应均值；N为样本数量。

3.2 外形近似模型

采用二阶响应面模型对于表 1的样本点进行拟合，得到R²为0.989 71，满足工程需要。抽取其中10个预测值与计算值进行对比，见图 6。

图 6中，10个预测值与计算值均非常接近。最终得到的结果如表 4。

3.3 耐压壳体近似模型

采用二阶响应面模型对于表 3的样本点进行拟合，σ_max和P_cr的R²分别为0.999 32、0.983 51，满足工程需要。抽取σ_max和P_cr的10个预测值与计算值进行对比，见图 7、图 8。

图 7，图 8中预测值和计算值均非常接近。得到的二次多项式系数如表 5。

4 优化模型的建立 4.1 第2代非支配排序遗传算法

多目标优化的目标是得到Pareto最优解集，即Pareto前沿。非支配排序遗传算法(NSGA)是遗传算法针对多目标优化的变形形式，第2代非支配排序遗传算法(NSGA-Ⅱ)是NSGA的改进形式^[19]，NSGA-Ⅱ通过不仅引入了精英策略，其还将拥挤距离引入其中，这样就保留了进化过程中的优良种群，使Pareto解集的解呈均匀分布和种群多样的分布状态。第2代非支配排序遗传算法的求解步骤为：

1)初始种群的赋予的同时，进行排序种群的排序，此代种群为P₀，其是随机产生的，初始种群进行选择、交叉与变异，进一步产生新种群，即Q₀种群，令t=0；

2)构造新种群，对于R_t进行非劣排序，得到非劣前端F₁，F₂，…；

3)按拥挤比较操作p_n对所有F_i进行排序，将最优秀的N个个体组成种群P_t+1；

4)对种群P_t+1执行选择、交叉和变异，进一步得到种群Q_t+1；

5)若满足优化条件，则计算程序终止，若不满足，则令t=t+1，并返回至2)重新执行。

4.2 优化模型

确定外形阻力R_x、极限载荷P_cr、质量M、排水体积Δ为无人潜水器主体结构优化的目标函数。其中，R_x、P_cr通过将表 4、表 5的近似模型系数代入式(7)得到，质量M、排水体积Δ通过以下得到:

$ \left\{ \begin{array}{l} M = {M_f}({x_2}, {x_3}, {x_4}, {\rho _f}, {t_f}) + {M_N}({y_1}, \cdots {y_5}, {\rho _N})\\ \Delta = {\Delta _f}({x_2}, {x_3}, {x_4}, {t_f}) + {\Delta _N}({y_1}, \cdots {y_5}, {\rho _N}) \end{array} \right. $

(9)

式中:ρ为密度；下标f为非耐压结构；N为耐压结构。以外形阻力R_x最小、排水体积最大、极限载荷P_cr最大、质量M最小为目标。建立主体结构多学科优化模型:

$ \left\{ \begin{array}{l} {\rm{max}}:\{ V;{P_{cr}}\} \\ {\rm{min}}:\{ M;{R_x}\} \\ {\rm{s}}{\rm{.t}}:\{ {\sigma _{{\rm{max}}}} \le \left[ \sigma \right]\\ DV:\begin{array}{*{20}{c}} {{x_2}, {x_3}, {x_4};}\\ {{y_1}, {y_2}, {y_3}, {y_4}, {y_5}} \end{array} \end{array} \right. $

(10)

式中[σ]为耐压壳体材料许用应力，MPa。其通过式(11)求出。

$ \left[ \sigma \right] = 0.85{\sigma _s} $

(11)

式中σ_s为壳体材料屈服极限，MPa。

5 优化求解

采用NSGA-II进行式(11)优化的求解，种群数量100、进化代数50、交叉0.9、交叉分布指数10、变异分布指数20。得到的Pareto解集见图 9。

图 9中2个方案与初始方案进行的对比见表 6。

表 6中，与初始方案相比，优化1、优化2在各个优化目标上均取得了优化的效果，但两者在优化目标的优化程度上有所不同。优化方案1对阻力与重量的优化效果更为明显，R_x、M分别降低5.75%、8.59%；优化方案2的P_cr提高104.87%，Δ提高2.13%。

6 结论

1)通过针对STAR-CCM+和Abaqus进行的二次开发实现了外形阻力、耐压壳体的参数化分析，可避免设计过程中频繁的模型修改，实现了设计效率的提高；

2)采用二阶响应面拟合得到分析模型具有极高的精度，进一步将分析模型模型引入到主体结构优化模型，在满足了工程需要的同时，降低了计算成本、提高了设计效率；

3)基于NSGA-II对于优化模型进行优化求解得到了Pareto解集，选取其中的2个方案与初始方案进行对比，结果表明：2个优化方案均实现了优化，且各自在不同性能的提升方面效果明显。