球壳是载人潜器重要的耐压结构形式之一，随着载人潜器的轻量化、低造价的技术发展要求^[1]，对于耐压球壳的设计提出了更高的要求。耐压球壳的设计目标是在保证其强度与稳定性的前提下，通过尺寸的优化达到减小其重量、降低建造成本的目标。然而，载人潜器耐压球壳是典型的多开孔结构，其中，人员进出口、观察窗等较大的开孔不仅会严重削弱耐压壳体的强度和刚度，还提高了其设计的难度。围绕着多开孔耐压球壳的设计，众多学者开展了研究工作，孙善萍等^[2]针对开孔耐压球壳的极限强度进行了非线性有限元分析，计算了开孔的大小及加强围壁参数的变化对壳体极限强度的影响；刘峰^[3]采用第二代非支配排序遗传算法，针对载人舱出入孔加强结构进行了质量和极限强度的多目标优化。余俊等^[4]针对开孔耐压球壳，运用多岛遗传算法和Hook-Jeeves法进行了优化。然而，载人潜器耐压球壳需要设置观察窗，观察窗所采用材料的有机玻璃为粘弹性材料，恒压作用下一定时间内容易产生蠕变变形，变形过大很可能造成窗玻璃挤出，进一步对于载人潜器的安全产生不利影响。杨青松等^[5]对于半球形观察窗接触、蠕变等问题进行了试验研究和数值模拟；Guo等^[6]建立了球形观察窗的广义Maxwell模型，利用ABAQUS软件进行了计算；Du等^[7]采用数值法建立了接触有限元模型，确定了观察窗应力破坏准则。Pranesh等^[8-9]开展了观察窗的优化研究。Wang等^[10]针对全海深载人潜水器观察窗完整窗口模型，进行了准静态加载-卸载循环试验，所得出的研究结论为全海深观察窗设计提供了依据和参考；黄浔等^[11]利用有限元分析方法对于观察窗蠕变特性进行了研究，总结了观察窗厚度-直径比、倾角和接触面摩擦系数等设计参数对于观察窗轴向位移和应力的影响。目前针对耐压球壳的优化设计中，未有将观察窗蠕变变形纳入优化范围，仅对球壳质量或极限强度的优化可能忽略了观察窗的大变形，从而带来一系列安全隐患。

本文在考虑观察窗蠕变特性的基础上，设计了耐压球壳参数化分析流程，利用该分析流程进行了样本点的选取及分析，在建立了基于蠕变效应的耐压球壳近似模型的基础上，进行了耐压球壳的优化的求解。

1 考虑蠕变特性的耐压球壳分析模型 1.1 观察窗结构形式

某载人潜器设有2个侧观察窗，1个主观察窗。由于各个观察窗之间球面距离满足《潜器规范》^[12]要求，独立且互不影响，因此，仅针对主观察窗进行研究。耐压球壳观察窗和窗座结构形式及主要设计参数见图 1。

图 1中，初始设计参数为：D₂=470 mm，D₁=410 mm，d=200 mm，d₁=220 mm，t=105 mm，α=β=45°。窗座与壳体均采用高强度钢材，观察窗为有机玻璃，材料参数见表 1。

1.2 观察窗蠕变位移求解

根据试验拟合出常温有机玻璃回归公式为^[13]：

$ \log \left(\frac{\varepsilon}{\varepsilon_{c}}\right)=0.036\;8\left[\log \left(\frac{t}{t_{c}}\right)\right]^{2}+0.278\;11 \log \left(\frac{t}{t_{c}}\right) $

(1)

式中：ε为应变；ε_c为临界断裂应变；t_c为临界断裂时间，ε_c和t_c与应力水平相关：

$ \varepsilon_{c}=1.4 \exp \left(-\frac{\sigma}{15}\right) $

(2)

$ t_{c}=(0.5 E+14) \exp \left(-\frac{\sigma}{1.68}\right) $

(3)

式中：σ为应力；E为弹性模量。

将式(2)、(3)代入式(1)，得到观察窗有机玻璃本构模型为：

$ \begin{aligned} \log \varepsilon=& 0.036\;8(\log t+0.021\;94 \sigma)^{2}+\\ & 0.278\;1(\log t+0.021\;94 \sigma)-0.028\;95 \end{aligned} $

(4)

1.3 有限元分析模型

观察窗与窗座之间接触定义采用接触对算法，设定软质材料观察窗为从面，高强度钢材窗座为主面，接触面摩擦系数为0.1。遵循从面网格尺寸小于主面网格，考虑计算量，观察窗全局种子取10，壳体取20，网格划分见图 2。分析步设为静力、粘性，分析时间为潜器标准作业时间6 h。

为消除耐压球壳刚体位移，采用3点约束球壳6个方向自由度，具体见图 3。

1.4 观察窗蠕变位移与耐压球壳强度求解

目标潜器工作水深1 500 m，根据《潜器规范》^[12]，安全系数取1.5，则计算压力P_j为24.5 MPa。耐压球壳表面施加24.5 MPa载荷进行蠕变分析，结果见图 4、图 5。计算得到的耐压球壳最大Mises应力为606.63 MPa，满足规范强度要求，观察窗最大蠕变位移为2.914 mm。

1.5 耐压球壳稳定性分析

基于上述有限元模型，采用弧长法进行耐压球壳极限载荷P_cr的计算，弧长法能很好对于位移-载荷曲线下降段进行捕捉，在结构的后屈曲分析中已经得到了广泛应用。计算中的初始缺陷通过模型特征值屈曲模态的提取得到，并对于几何以及材料非线性进行考虑。将材料的名义应力应变进行转化、拟合得到真实应力应变曲线。替换模型分析步为屈曲和线性摄动，进行单位载荷施加，分析特征值屈曲。提取得到前六阶模态，部分模态如图 6所示。

图 6中，各阶模态特征值相差不大，可见，模型对初始缺陷比较敏感，说明在添加初始缺陷时提取全部六阶模态的节点位移来使计算结果偏于保守，模态比例因子通常取壳体厚度的1%。再次替换模型分析步为静力、弧长，载荷根据第六阶模态特征值取109.85 MPa，进行后屈曲分析，球壳失稳破坏形式如图 7，弧长-载荷比例因子曲线见图 8。

图 7、图 8计算得到的耐压球壳极限载荷P_cr=33.933 MPa，满足规范稳定性要求。

2 考虑蠕变的耐压球壳参数化分析 2.1 基于Python语言的Abaqus二次开发

命令执行前借助Python解释器生成.rpy格式文件提交计算将生成inp格式文件，通过分析求解生成便于后处理的odb文件。

Abaqus创建的模型包括根(root)对象：进行调整模型视角等的视图模块session；存储模型树中包含的部件、材料、装配体、分析步、载荷等数据模型的数据模块mdb和odb；存储所有模型信息、输出数据等的数据输入输出模块odb。每个根对象下面又有很多命令分支，通过编写脚本命令可以准确地对不同对象进行访问和编辑。

2.2 Abaqus的集成

在iSight软件进行Abaqus软件集成，主要包括以下阶段：

1) 前处理阶段(模型参数化)。对所研究问题的模型进行定义，包括模型的几何尺寸、材料参数、载荷边界等进行参数化处理，得到Abaqus输入文件(.inp, .py)。最为简便直接的是操作Abaqus/CAE建立三维模型，对于模型录制.rpy格式文件并修改后缀为.py的Python脚本文件，或直接编写Python的代码作为Abaqus输入文件；

2) 模拟计算阶段。在iSight中使用.bat格式文件调用Abaqus/Standard或Abaqus/Explict求解器，后台自动计算得到输出文件(job.odb、job.dat等)；

3) 后处理。在计算结果文件得到后，进一步通过Abaqus或其他后处理软件进行结果数据读取，或者通过iSight读取输出文件中输出的数据，进行结果分析或建立近似模型。

2.3 参数化分析的实现

iSight软件提供了与Abaqus软件的接口，在iSight软件进行Abaqus软件集成主要通过文件进行连接，具体步骤及方法为：

1) 蠕变计算。输入.py文件。文件前处理部分在Abaqus软件用户界面上操作生成，包含有限元模型和计算设置信息等。数据提取部分主要提取蠕变计算结果中模型总的最大Mises应力、模型质量以及目标位置的位移等信息；采用.bat批处理文件调用Abaqus进行后台计算；输出文件为包含模型所有节点Mises应力信息和最大、最小Mises应力数据的.rpt报告文件；以及包含由.py文件取到的目标位置中的蠕变位移和模型质量的.txt文本文件。

2) 线性屈曲分析。输入文件为.py文件。在得到.rpy文件后，修改文件后缀为.py。文件中的前处理部分在前一步蠕变分析的基础上，对关键词、替换分析步设置等进行修改，输出部分为分析结果文件中一阶模态的节点位移信息，作为模型初始缺陷数据；采用.bat批处理文件调用Abaqus进行后台计算；输出文件为包含了线性屈曲计算结果中所得到的各阶模态特征值的.dat文件，特征值与载荷的乘积即为随后进行的非线性分析中的载荷值。

3) 非线性屈曲分析。对于.py、.rpy等格式的文件操作以及修改、设置等操作与第2步相同。输出部分主要涉及极限载荷的结果，输出模型弧长-载荷因子曲线数据信息，提取最大载荷因子；采用.bat批处理文件调用Abaqus进行后台计算；输出包含模型载荷因子-弧长曲线数据信息的.rpt报告文件。

将上述输入、批处理、输出文件写入iSight软件Simcode模块中，通过iSight对于Abaqus软件进行驱动，实现模型参数化并依次进行强度分析、蠕变分析、线性屈曲分析及非线性分析。图 1所示的耐压球壳观察窗结构几何参数，即设计变量有α、β、t、D₂、D₁、d、d₁，输出为耐压球壳极限载荷压力P_cr、观察窗临界应力σ_w、球壳临界应力σ_b、质量M以及观察窗蠕变位移U。上述设计变量中，D₂、D₁、d、d₁由球壳布置决定，参数α、t会影响β，也会使U、P_cr和M发生变化，因此，确定α、t作为设计参数进行参数化，取值范围分别为：38°≤α≤60°，105 mm≤t≤150 mm。

3 耐压球壳优化求解 3.1 近似模型的建立

神经网络具有强大的处理复杂非线性问题的能力、高容错度和无须数学假设等优势，已成为解决许多领域中的复杂优化问题的有力工具。神经网络由输入层、中间层和输出层3部分构成，通过非线性变换完成输入矢量到新空间之间的映射，即中间层空间。中间层空间到输出层空间则为线性映射，输出层在新的线性空间中进行线性加权组合。径向基(radial basis functions, RBF)网络函数为^[14]：

$ f(x)=\sum\limits^{m} \lambda_{i} \phi\left(r^{0 i}\right)+\sum\limits^{K} C_{i} P_{i}(x) $

(5)

式中：m为样本点数；x为设计变量向量；λ_i(i=1, 2, …，m)为待定加权系数；ϕ(r⁰ⁱ)为基函数；r⁰ⁱ=‖x－xⁱ‖为插值函数中任意一点到第i个插值点的欧几里得距离；P_i(x)为多项式项；K为多项式项个数；c_i(i=1, 2, …，K)为P_i(x)所对应系数。

基函数一般选用高斯函数：

$ \varphi\left(r^{o i}\right)=\exp \left(-\frac{\left(r^{o i}\right)^{2}}{c^{2}}\right) $

(6)

响应面模型(response surface methodology，RSM)是一种对于试验样本点进行拟合的多项式函数。RSM的表达式为：

$ \begin{aligned} \hat{y}(x)=& a_{0}+\sum\limits_{i=1}^{n} b_{i} x_{i}+\sum\limits_{1 \leqslant i \neq j \leqslant n}^{n} c_{i j} x_{i} x_{j}+\\ & \sum\limits_{i=1}^{n} d_{i} x_{i}^{2}+\sum\limits_{i=1}^{n} e_{i} x_{i}^{3}+\sum\limits_{i=1}^{n} g_{i} x_{i}^{4} \end{aligned} $

(7)

式中：$\hat{y}(x)$为近似值；x_i为第i个设计变量；n为设计变量数量；a₀、b_i、…、e_i、g_i为拟合系数。

误差分析采用工程领域常用的R²分析，R²取值范围为[0, 1]，越接近1则表明模型拟合误差越小，计算公式为：

$ R^{2}=1-\frac{\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}} $

(8)

式中：$y_{i} 、\hat{y}_{i}、\bar{y}_{i}$分别为第i个样本点所对应状态的响应值、近似值、响应均值；n为样本数。

不同设计变量取值范围存在差异，量级和量纲的不同往往会降低优化效率，甚至影响优化结果。将设计变量归一化，使各设计变量均位于[-1, 1]，建立起以原点为中心的立方形设计空间，便于统一处理，有利于解决量纲差异带来的麻烦。本文所选取的设计变量归一化处理为：

$ x_{1}=\frac{t-127.5}{22.5}, x_{2}=\frac{\alpha-49}{11} $

(9)

构建近似模型需要通过试验设计获取一定数量的样本点，最优拉丁超立方设计(optimal Latin hypercube design, OptLHD)是基于拉丁方设计(Latin hypercube design, LHD)的改良版，其能使所有试验点尽可能地均匀分布于采样空间，良好的填充性和均衡性使得该方法更加精准真实^[15]。选取30个样本点进行计算，分别采用二阶响应面模型、四阶响应面模型和径向基神经网络模型对样本点进行拟合，不同近似模型拟合精度见表 2。

表 2中，径向基神经网络模型具有很高的精度，满足工程实际意义，因此，选取径向基神经网络模型建立耐压球壳的近似模型。

3.2 单目标优化

以观察窗倾角和厚度为设计变量，耐压球壳质量最小、观察窗蠕变位移最小为目标函数，壳体临界应力σ小于0.85倍的材料屈服强度σ_s、极限载荷P_cr大于计算压力P_j为约束条件，建立耐压球壳单目标优化模型为：

$ \left\{\begin{array}{lc} \text { Optimization : } & \{\min : M\} \\ \text { s. t } : \sigma_{\text {cr }} \leqslant 0.85 \sigma_{s} ; & P_{\text {cr }} \geqslant P_{j} \\ \text { DV }: & x_{1}, x_{2} \in[-1,1] \end{array}\right. $

(10)

多岛遗传算法(multi island genetic algorithm，MIGA)是一种全局优化算法，基于传统遗传算法对并行分布遗传算法有所改进。MIGA将种群分成多个子群，又称之为“岛”，每个岛上子群独立进化，而非全部种群采用相同的进化机制，一定时间间隔进行“岛屿”间的“迁移”，完成信息交换^[16]。MIGA能够有效提高运算速度，并且若干独立进化的子群丰富了整个种群的遗传多样性，避免落入局部解。图 9为算法流程图。

采用MIGA，子群规模为10，岛数为20，代数为10，交叉概率1.0，变异概率0.01，计算结果与有限元计算结果对比见表 3。

基于近似模型的优化结果与有限元结果对比，质量相对误差为0.13%，位移相对误差为2.25%，说明近似模型具有较高的精度。优化结果较初始设计质量减小0.25%，而蠕变位移增大6.42%，说明不考虑观察窗蠕变效应时，质量优化效果不明显，而蠕变变形明显增大，因此，需要将观察窗蠕变效应作为目标函数，进一步开展多目标优化设计。

3.3 多目标优化

以观察窗倾角和厚度为设计变量，耐压球壳质量最小、观察窗蠕变位移最小为目标函数，壳体临界应力σ小于0.85倍的材料屈服强度σ_s、极限载荷P_cr大于计算压力P_j为约束条件，建立耐压球壳单目标和多目标优化模型为：

$ \left\{\begin{array}{lc} \text { Optimization: }\{\min : M\}, & \{\min : U\} \\ \text { s. t} : \sigma_{\mathrm{cr}} \leqslant 0.85 \sigma_{s} ; & P_{\mathrm{cr}} \geqslant P_{j} \\ \text { DV }: & x_{1}, x_{2} \in[-1,1] \end{array}\right. $

(11)

第2代非支配排序遗传算法NSGA-II^[15]改进了非支配排序方法，选择接近Pareto前沿的个体，增强了Pateto的前进能力，拥有良好的探索能力。导入了“拥挤距离”和“拥挤距离排序”的方法，在具有同样的Pareto顺序的层内，能够对个体进行排序，称为拥挤距离排序。进化过程中，当前父代群体通过交叉和变异生成子群体，将2个群体合并。在目标空间中遵照Pareto最优关系将群体中个体两两按其目标函数向量进行比较，将群体中所有个体分成多个依次控制的前沿层。在属于不同的Pareto层的情况下，利用评价Pareto优越性来评价个体的优劣。属于相同Pareto层的个体，具有更大拥挤距离的个体更优秀。交叉和变异运算机制使用SBX方法。根据该方法生成子个体的交叉运算为：

$ \left\{\begin{array}{l} x_{i}^{(1, t+1)}=0.5\left[(1+\beta) x_{i}^{(1, t)}+(1-\beta) x_{i}^{(2, t)}\right] \\ x_{i}^{(2, t+1)}=0.5\left[(1-\beta) x_{i}^{(1, t)}+(1+\beta) x_{i}^{(2, t)}\right] \end{array}\right. $

(12)

式中：x_i^{(1, t+1)}、x_i^{(1, t)}、x_i^{(2, t)}对应于个体x_i的维数分量；随机分量β对应相应维数进行计算重新取值：

$ \beta=\left\{\begin{array}{ll} (2 u)^{1 /\left(\eta_{c}+1\right)}, & u \leqslant 0.5 \\ {[2(1-u)]^{1 /\left(\eta_{c}+1\right)},} & u>0.5 \end{array}\right. $

(13)

式中：u在[0, 1]中随机取值；η_c为常数，表征交叉分布系数，大小直接决定了子代个体与父代个体的接近程度。

根据该方法生成子个体的变异运算为：

$ x_{i}^{(1, t+1)}=x_{i}^{(1, t)}+\delta\left(x_{i}^{U B}-x_{i}^{L B}\right) $

(14)

$ \delta_{q}=\left\{\begin{array}{l} {\left[2 u+(1-2 u)\left(1-\delta_{1}\right)^{\eta_{m}+1}\right]^{1 / \eta_{m}+1}-1,} \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ u \leqslant 0.5 \\ 1-\left[2(1-u)+2(u-0.5)\left(1-\delta_{2}\right)^{\eta_{m}+1}\right]^{1 / \eta_{m}+1}, \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ u>0.5 \end{array}\right. $

(15)

$ \delta=\min \frac{x_{i}-x_{i}^{L B}, x_{i}^{U B}-x_{i}}{x_{i}^{U B}-x_{i}^{L B}}, u \in[1,0] $

(16)

式中：$\delta_{1}=\frac{x_{i}-x_{i}^{L B}}{x_{i}^{U B}-x_{i}^{L B}} ; \delta_{2}=\frac{x_{i}^{L B}-x_{i}}{x_{i}^{U B}-x_{i}^{L B}}$；η_m为常数，表征变异分布系数。

采用NSGA-II算法，种群数为100、代数为50、交叉概率0.8、变异概率0.1。计算得出Pareto解集如图 10。图中，U与球壳M基本呈负相关，Pareto解集分布均匀，在耐压球壳观察窗设计过程中，设计者可根据实际情况进行设计方案的选择。选择若干优化方案与初始设计方案对比，列于表 4。

表 4中，各优化方案与初始设计方案相比，球壳质量稍有增大，但增幅不大，最大仅为2.07%，而观察窗蠕变变形有很大的改善，方案2减小最为明显，达到-18.09%。综合来看，在载人潜器耐压球壳优化设计时，考虑观察窗蠕变效应，能大幅减小蠕变变形，减小潜器的安全隐患，使其安全性能显著提高，而对球壳质量增加不大。

4 结论

1) 通过对于Abaqus软件进行二次开发和集成，实现了考虑蠕变特性耐压球壳的参数化建模与自动分析，提高了耐压球壳设计效率。

2) 得到了满足工程需要的近似模型，统一了设计变量的量级和量纲，计算表明合理选择近似模型及处理设计变量能够提高耐压结构设计和优化效率。

3) 考虑蠕变效应与未考虑时相比，虽然耐压球壳质量有所增加，但观察窗蠕变变形减小明显，说明了在载人潜器耐压球壳设计优化过程中，应加以考虑观察窗蠕变效应，有利于改善观察窗变形不协调问题，从而提高载人潜器的安全性。