耐压结构为载人潜水器提供了布置空间，对其安全性和总体性能等均有重要的影响．传统的耐压结构设计多采用确定性方法，对于设计中的不确定因素不予考虑．随着耐压结构安全性要求的提高，众多学者对耐压结构设计中的不确定性开展了研究．FAULKNER[1]对影响耐压结构的不确定性因素进行了研究，分析了随机变量的敏感性，确定了结构破坏模式；张伟等[2]总结了设计参数对耐压柱壳结构失效的影响；操安喜等[3]研究了影响耐压球壳可靠性的不确定参数和可靠性指标的获取方法．研究表明：不确定性可能使设计产生细微偏差或波动，从而造成性能的不稳定，甚至存在破坏的风险．Leungtsui[4]在优化过程中考虑了不确定性因素，通过对设计参数的优化使目标响应对参数变化的敏感度降低，可减小性能波动[5-6]．6σ优化是稳健性的重要方法之一[7-8]，通过对设计参数的合理组合降低产品性能对不确定性因素波动的敏感性，可获得高稳健性的产品．程妍雪等[9]考虑了不确定性影响，将6σ优化应用于耐压柱壳；甄春博等[10]将6σ优化应用于圆碟形耐压结构；李彬等[11]将6σ应用于环肋复合材料耐压结构，以上研究对象均为无开孔耐压结构．载人潜水器耐压结构的多个较大开孔，会导致其极限载荷降低和应力集中，须要进行加强．目前开孔加强的研究主要围绕形式对比[12]、参数对结构整体的极限强度影响[13]等开展．刘峰等[14]研究了单开孔耐压球壳的不确定性优化，研究表明为提升球壳的稳健性，须大幅提升其质量．载人潜水器耐压结构的多个大开孔的存在，还使其设计参数增加，使设计过程中的不确定性因素增加，导致问题研究难度增大．此外，耐压结构设计是一个反复迭代、不断优化的过程，如何在保证计算精度的基础上，降低计算成本也是设计中亟待解决的问题之一．近似模型通过数学模型的方法逼近输入变量与输出变量，用于替代真实模型，可有效解决计算精度与计算成本的矛盾．这里以多开孔耐压柱壳为研究对象，在对设计变量进行了灵敏度分析及降维处理的基础上，建立多开孔耐压柱壳近似模型，以质量最小、极限强度最大为优化目标，进行基于多目标确定性优化、3σ优化和6σ优化的求解，以得到更加稳健的方案．1 耐压柱壳设计参数与分析方法1.1　耐压结构形式耐压柱壳工作水深1 000 m，共有4个较大开孔，结构形式及开孔布置情况见图1．图1中，开孔1、开孔2、开孔3采用围壁加强形式，开孔4采用窗座加强形式．以耐压柱壳质量M最小和极限强度Pcr最大作为目标函数；约束为最大应力σmax、最大周向应力σ1、最大轴向应力σ2、最大肋骨应力σ3；设计变量为耐压柱壳壁厚x1、环肋厚度x2、环肋高度x3、开孔围壁厚x4、开孔1围壁高度x5、开孔周围壳体厚度x6、开孔2围壁高度x7、开孔3围壁高度x8、开孔4加强高度x9、开孔4加强角度x10．10.13245/j.hust.210922.F001图1耐压柱壳结构图1.2　耐压结构分析方法对于计算载荷Pj，按每个标准大气压条件下压力为9.8 kPa，可换算成Pj=0.009  8hj (1)式中hj=Khjx为计算深度，K为安全系数，取1.5；hjx=hg/0.9为极限下潜深度，hg为潜水器工作下潜深度．利用式(1)计算得到计算载荷为16.33 MPa，σmax，σ1，σ2和σ3的临界值分别为785，902.75，667.25和471 MPa．采用Abaqus软件进行分析，有限元分析模型见图2．将柱壳简化为固定壳体的刚体位移，约束首、尾两半球端点X，Y轴，中纵剖面任意节点处Y，Z轴的位移均设置为零．稳定性分析包括线性和非线性两部分，先进行线性分析，然后采用弧长法与牛顿-拉法松法结合求解非线性屈曲．10.13245/j.hust.210922.F002图2耐压柱壳有限元模型2 近似模型的建立2.1　样本点的选取最优拉丁超立方方法(Opt LHD)完善了拉丁超立方样本分布的均匀性，将n维空间中的每一维坐标区间[xkmin,xkmax]   (k∈[1,n])等分为m个区间，每个区间记为[xki-1,xki]   (k∈[1,m])．构成n维空间、样本数为m，记为m×n的Opt LHD，Opt LHD可使构建的近似模型更加精确，拟合效果更好．2.2　变量灵敏度分析采用Opt LHD在设计变量可行域内选取点，进行设计变量对Pcr和M的灵敏度分析，见表1．表1中，x3对Pcr的影响最大，随后依次为x1，x2，x5，x6，x4，x8，x7，x1，x9；x1对M的影响最大，随后依次为x2，x4，x3，x6，x9，x5，x10，x7，x8．综合设计变量灵敏度分析结果，将x1，x2，x3，x4，x6作为优化的变量，其余变量设置为固定值．10.13245/j.hust.210922.T001表1设计变量对Pcr和M的灵敏度分析参数PcrMx10.275 50.593 7x20.247 00.188 3x30.321 40.078 3x40.012 90.085 8x50.072 80.931 0x60.003 50.018 1x70.005 40.006 0x80.010 90.000 8x90.001 90.016 1x100.002 10.008 02.3　近似模型及精度判断响应面模型(response surface methodology，RSM)利用多项式函数建立输入变量与输出变量之间的函数关系[15]，RSM可分为一阶、二阶、三阶、四阶，RSM的表达式为ŷ(x)=a0+∑j=1Qa1jxj+∑j=1Qa2jxj2+∑j=1Qa3jxj3+∑j=1Qa4jxj4+∑j=1,k≠jQajxjxk,式中：ŷ(x)为近似值；a为响应面模型系数；xj为设计变量；Q为设计变量个数．径向基神经网络模型(radial basis function，RBF) 第i个隐藏层单元输出的响应为hi=exp[-xl-ci2/(2σi2)](1≤i≤l)，式中：ci和σi分别为第i个隐藏层的中心、单元实际的宽度；xl为第l个输入样本．得到神经网络输出层第j个输出为f(x)=∑i=1bwijhi，其中wij为隐藏层i节点的第j个输出值所占的权重．Kriging模型为通过一组已知独立变量和对应的系统响应，构建变量与响应之间近似关系的一种半参数化插值模拟技术[16]．设未观测到的须估值点x0周围的观测点为x1，x2，…，xn，对应观测值为 y(x1)，y(x2)，…，y(xn)，则x0估计值为y˜(x0)=∑i=1nλiy(xi)，其中λi为未知的、须待定的加权系数，须符合无偏估计和方差结果，具体分别为∑i=1Nλi=1；D[y˜(x0)-y(x0)]=minD[y˜(x0)-y(x0)]=-∑i=1N∑j=1Nλiλjγ(xi,xj)+2∑i=1Nλiγ(xi,x0),式中：γ(xi，xj)为xi和xj间距离为h的情况下，参数的半方差值大小；γ(xi，x0)为xi和x0间距离为h的情况下，参数的半方差值大小．利用复相关系数(R2)判断近似模型的拟合精度R2=1-∑i=1n(yi-ŷi)2/∑i=1n(yi-y¯i)2，式中：yi，ŷi，y¯i分别为第i个样本状态变量的响应值、近似值和响应均值；n为样本的数量．2.4　耐压柱壳近似模型采用Opt LHD选择35个样本点进行分析，部分样本点见表2，质量M根据公式计算得到．10.13245/j.hust.210922.T002表2耐压柱壳样本点序号x1/mmx2/mmx3/mmx4/mmx6/mmσmax/MPaσ1/MPaσ2/MPaσ3/MPaPcr/MPa124.9432.00100.8847.0640.71613.7579.6620.9304.119.2223.1827.6897.9448.2944.00610.9583.0610.6332.915.8︙922.9425.0098.8249.9440.50590.4564.9616.1347.014.91022.5931.7999.4149.7442.76595.6570.3614.1316.817.3︙2425.5326.85102.9446.2442.35623.9584.0624.8325.018.3︙3422.1227.0696.1846.4440.09649.3620.5645.6344.414.43523.0625.82102.3545.4141.53655.9624.2650.7345.116.3采用RBF、Kriging、二阶RSM、三阶RSM、四阶RSM对样本点进行拟合，拟合精度见表3，耐压柱壳近似模型系数见表4．进一步利用上述近似模型对Pcr，σmax，σ1，σ2，σ3进行求解得到20组预测值与有限元分析值进行对比，见图3，图中：δ为相对误差；N为点数．表3中Pcr，σmax，σ1，σ2，σ3拟合精度最高的依次为RBF、二阶RSM、四阶RSM、三阶RSM、二阶RSM．对于Pcr，二阶RSM与RBF拟合精度相差仅为0.002，由于二阶RSM在模型表达上具有优势，采用二阶RSM建立Pcr的近似模型；采用二阶RSM建立σmax近似模型；采用二阶RSM建立σ3的近似模型；采用四阶RSM建立σ1近似模型；采用三阶RSM建立σ2近似模型．从图3可以看出：只有两个点Pcr的预测值误差在2%左右，其余均小于2%，大部分处于1%左右，进一步验证了近似模型拟合精度较高．10.13245/j.hust.210922.T003表3近似模型拟合精度方法Pcrσmaxσ1σ2σ3RBF0.9810.9960.7180.9460.824Kriging0.9170.9820.9890.9630.259二阶RSM0.9790.9970.9840.9580.968三阶RSM0.9600.9920.9930.9820.959四阶RSM0.9720.9820.9980.9340.83510.13245/j.hust.210922.T004表4近似模型系数参数σmax/103σ1σ2σ3/103Pcr/103常数4.15-2.7×1069.8×10691.003.759x6-187.00-3.0×106-1.6×10532.511.041x2-18.205.5×106-1.4×105-73.605.878x3-0.149.8×107-3.8×107-232.00-9.689x1-19.00-1.1×106-7.8×1052 928.008.698x4-91.801.2×107-1.2×104-149.00-3.537x6216.301.1×1086.1×10627.024.136x22-25.60-2.8×1083.6×106167.00-4.792x32-1.22-1.3×1074.9×107141.906.290x12662.806.6×1071.9×107277.6013.110x42513.60-3.6×1082.3×105117.903.480x6x2-0.65-4 016.43.6×10418.003.920x6x320.47-2 733.7-9.2×104-46.60-1.900x6x122.25-3.2×104-9.1×10485.77-1.750x6x48.803.7×104-6.4×104-25.300.487x2x324.43-6.9×1046.6×10481.80-7.340x2x135.804.0×1049.5×10439.8015.300x2x4-3.20-7 977.9-3.6×105-99.60-1.120x3x1-40.056.0×1044.9×105-82.10-11.900x3x430.66-3.2×104-2 783.6162.904.029x1x4175.403.7×105-2.5×106232.103.475x63—-1.8×109-5.3×107——x23—6.7×109-3.9×107——x33—5.4×106-2.1×107——x13—-1.9×109-2.5×108——x43—5.0×1097.6×106——x64—1.2×1010———x24—-5.9×1010———x14—1.9×1010———x44—-2.6×1010———10.13245/j.hust.210922.F003图3预测值精度3 确定性多目标优化耐压柱壳确定性多目标优化模型为：Opt{min:M;max:Pcr};  s.t.{σmax;σ1;σ2;σ3};DV{x1;x2;x3;x4;x6},式中DV表示设计变量．NSGA-II是非支配排序遗传算法的改进[17-18]，求解步骤为：a．随机产生初始种群Q0，进行秩的赋予，随后进行种群的非劣排序；b．进行选择、交叉与变异等操作完成Q0筛选，得到新的种群Q0；c．构造新种群Rt，排序得到非劣前端F1，F2，⋯；d．对于所有的Fk，按照拥挤比较操作排序，将其中最优秀的N个个体组成新种群Qt+1；e．对Qt+1执行选择、交叉和变异等操作，得到新种群Qt+1；f．若满足优化条件，则终止计算程序，若不满足，则令t=t+1，并返回至c．NSGA-II的设置为：种群数量为100，进化代数为50，交叉可能性为0.75，交叉分布指数为10，变异分布指数为100，得到的确定性优化Pareto前沿见图3．4 多目标稳健性优化多目标稳健性优化在多目标优化基础上，考虑了设计变量偏差，同时相关目标增加了均值和方差，多目标稳健性优化模型为      min{f1(Yμ1(X),Yσ1(X)),f2(Yμ2(X),Yσ2(X)),⋯,        fk(YμK(X),YσK(X))};      s.t.  gμj(x)+ηgσj(x)≤0,xL+ηxσ≤xμ≤xU-ηxσ,式中：Yμ1(X)，Yμ2(X)，…，Yμk(X)和Yσ1(X)，Yσ2(X)，…，YσK(X)分别为第K个目标函数的均值和方差；gμj(x)和gσj(x)分别为第j个约束的均值和方差；xμ和xσ分别为x的均值和方差．利用蒙特卡罗进行稳健性评价时，所有不确定性因素概率分布已知，进一步对均值、方差等的概率分布进行估计．均值μ、方差σ分别为μY=1N∑i=1Ng(xi);σY2=1N-1∑i=1Ng(xi)-μY2.考虑参数几何尺寸的物理不确定性[19]，将制造误差引入作为随机变量，认为设计变量服从正态分布，变异系数制造误差定义为0.01μ，变异系数为0.01，不确定因素为随机几何尺寸参数，x1，x2，x3，x4，x6的设计空间分别为[22，26]，[25，32]，[75，125]，[45，52]和[35，42] mm．耐压柱壳多目标稳健性优化模型为 Opt  {min:M=μM+ωσ2M,-Pcr=-μPcr+ωσ2Pcr};     s.t.{μσmax+nσσmax≤785MPa,        μσ1+nσσ1≤902.75MPa,        μσ2+nσσ2≤667.25MPa,        μσ3+nσσ3≤471MPa};     DV{xLi+nσxi≤xi≤xUi-nσxi       (i=1,2,3,4,6)},式中：n=3或6；ω=0.01；xLi和xUi分别为第i个变量的最小值和最大值．NSGA-II设置与确定性优化相同，得到3σ优化、6σ优化Pareto前沿见图4．图4中解集的分布范围按照确定性优化、3σ稳健性优化、6σ稳健性优化的次序变小，在M相等的情况下，Pcr按照以上次序增加．10.13245/j.hust.210922.F004图4Pareto前沿采用最小距离选择的方法(TMDSM)求解Pareto解集中的最优解，表达式为[20]minD=∑τ=1Kfcτ-min(fτ(x))d1/d，式中：K为目标的数目；fcτ为非劣解集中第c个非劣解中第τ个目标分量；d=2，4，…，D为拐点与理想点的距离．对于确定性优化和稳健性优化的优化结果，利用最小距离法确定相应的优化方案，如表5所示．10.13245/j.hust.210922.T005表5设计方案对比参数确定性3σ6σx1/mm22.3824.8125.13x2/mm32.0031.1031.11x3/mm105.00104.03105.11x4/mm43.9148.0147.75x6/mm41.9639.7539.83σmax/MPa697.84600.33599.62σ1MPa689.17568.32569.08σ2/MPa660.12601.40603.73σ3/MPa329.21302.22303.44Pcr/MPa19.4120.1420.74M/(103 kg)14.592 1115.763 2015.882 19M质量水平1.331.976.20M均值/(103 kg)14.950 2615.761 3215.881 10M标准差/kg302.11256.6199.42Pcr质量水平1.722.036.09Pcr均值/MPa19.9720.1420.74Pcr标准差/MPa3.853.162.07表5中3σ方案与确定性优化方案相比M增加 1.170 45×103 kg，Pcr增加0.732 MPa，两者的质量水平也有所提高，σmax，σ1，σ2，σ3均明显降低；6σ方案与3σ方案相比，σmax，σ1，σ2，σ3均变化不大，M增加119.47 kg、Pcr增加0.594 MPa，两者的质量水平提升明显，设计变量变化程度按照x1，x4，x6，x2，x3的次序降低．为进一步分析优化结果的变化情况，对设计变量对约束的灵敏度分析见表6．10.13245/j.hust.210922.T006表6设计变量对约束灵敏度分析参数x1x2x3x4x6σmax-0.229 2-0.027 7-0.054 6-0.650 0-0.028 5σ1-0.262 3-0.036 0-0.081 7-0.539 9-0.080 1σ2-0.169 8-0.025 70.118 7-0.665 50.020 3σ3-0.228 8-0.276 40.270 5-0.075 6-0.088 8结合表1和表6进行分析：x1 对Pcr的影响仅次于x3，对σmax，σ1，σ2，σ3的影响较大，增加x1可降低σmax，σ1，σ2，σ3；x2对σmax，σ1，σ2影响很小，对Pcr的影响小于x1；x3对Pcr的影响最大，为提高耐压结构的稳健性，须提高x3，而x3增加也将导致σ3降低，不利于结构的安全性；x4对Pcr和σ3的影响不大，对σmax，σ1，σ2影响很大，提高x4可降低σmax，σ1，σ2，σ3，从而有利于结构安全；x6对Pcr，M，σ2，σmax，σ1，σ3的影响均不大；因此，为提高Pcr，x1的变化幅度最大，导致质量增加明显，对表5各方案的目标函数质量进行分析见图5，图中p为概率密度．10.13245/j.hust.210922.F005图5目标函数质量分析图图5中6σ方案与3σ方案相比，6σ方案的图形更为修长，说明 6σ方案的目标函数在均值附近更为集中，其在均值附近的概率也最大；3σ方案离下约束边界更近，说明3σ方案较易随着设计变量的波动出现不满足约束条件的情况，所以6σ方案更加稳健；同样3σ方案的稳健性明显优于确定性优化方案．进一步对确定性优化、3σ稳健性优化、6σ稳健性优化进行系统稳健性分析，见表7．10.13245/j.hust.210922.T007表7系统稳健性分析参数变量确定性3σ方案6σ方案σ水平Pcr1.952.926.617σmax8.08.08.0σ10.6873.8347.778σ28.08.08.0σ33.778.08.0x12.0058.08.0x20.7363.0907.403x30.5293.6036.519x42.3468.08.0x60.7368.08.0系统0.0893.6246.569方案成功率Pcr0.9770.9950.996σmax1.01.01.0σ10.3720.9960.997σ21.01.01.0σ30.9921.01.0x10.9551.01.0x20.50.9980.999x30.4030.9230.968x40.9811.01.0x60.5381.01.0系统0.0720.4580.962每百万缺陷Pcr24 0004 0002 000σmax000σ162 6004 0002 000σ2000σ3000系统9.28×1055.46×1055.10×105表7中3σ优化、6σ优化在所有的设计变量、约束函数、目标函数、系统的质量水平，以及方案的成功率均有所提高的同时，每百万缺陷数进一步下降，6σ稳健性优化的系统稳健性改善最为明显，多数设计变量的质量水平提高到6σ以上．通过研究主要得出如下结论．a．通过耐压结构设计变量维数的降低、高拟合精度近似模型的建立，可在保证计算精度的条件下，降低问题分析难度、提高设计效率．b．为提高结构的稳健性必须增加厚度，导致稳健性的优化质量增加幅度较大，因此在耐压结构的设计、建造中须重视壳体厚度的精度控制．c．相比确定性多目标优化，稳健性多目标优化Pareto解集的分布范围更小，解的数量上也明显减少，说明稳健性优化波动较小，方案更加稳健．d．尽管6σ方案质量最大，但稳健性最好，且方案的可靠性、方案的成功率和总体稳健性均最高，即使考虑随机波动因素也可保证稳健性需求．