1 引言
浩瀚的海洋蕴藏着各类丰富的资源，随着开发利用海洋资源的进程加快，对各类无人航行器的要求日益提升。其中，水下滑翔机能够依靠调整自身重浮力的大小，以“之”字形的运动轨迹在不同的剖面收集海洋信息，这种航行器能耗较低，但其运动单一，抗干扰力差，工作效率低，无法胜任形式多样的探测任务[1]。水下多旋翼航行器通过推进器差速可以实现各种复杂的运动，适应力强，但由螺旋桨旋转产生的噪声与振动无法避免，且受限于能源供应而无法长时间、远距离地进行作业[2]。其次，潜航器的水动力性能直接取决于其外形，许多学者为提高AUV的运动性能，对外形的优化设计作了深入研究，如：吴利红等[3]通过CFX软件计算了不同攻角和漂角工况下滑翔机所受的位置力，并对滑翔机展开了线型优化；赵凯凯等[4]利用CFX软件对滑翔蛇形机器人进行了水动力仿真，求解了水动力参数，并采用矩形平板翼对翼型进行了优化；李天博等[5]以最大化滑翔机的升阻比为优化目标，建立了Kriging代理模型，采用EGO算法对联翼式滑翔机外形进行了优化设计；楚灯旺等[6]通过FLUENT软件分析了4种浮标的外形方案，建立了壳体外形阻力模型，采用基于惩罚函数的遗传算法优化了壳体线型参数；张代雨等[7]通过CAESES软件对翼身融合滑翔机进行全参数化建模，以升阻比为优化目标，并基于多种组合优化算法提高了外形的水动力性能；高启升等[8]提出了基于组合加权响应面法的多目标优化方法，以折衷规划法对子目标进行归一化处理，对梯形肋骨耐压结构进行了多目标优化设计。
在之前的研究基础上，提出了一种四螺旋桨驱动与重浮力调整装置相结合的混合驱动无人潜航器整体结构方案，通过FLUENT软件分析了AUV在不同攻角和速度下的水动力系数，利用拉丁超立方采样试验设计确定了样本点，建立了二阶响应面模型，以升阻比最大和质量最小为优化目标，采用多目标遗传算法对响应面模型进行优化求解，得到了AUV的优化设计方案，为AUV水动力外形进快速优化设计提供一种高效的方法。
2 结构方案
AUV采用模块化设计，其结构组成可分为壳体、多传感器阵列单元、浮力调节单元、重心调节单元、吊装背鳍、外挂设备单元、主控与驱动单元、推进器单元、尾翼和电力电子单元，结构组成如图1所示。
图1 AUV结构组成示意图
Fig.1 AUV structure composition
其中浮力调节单元通过步进电机带动活塞移动改变单元舱内进水量实现净浮力的调节功能。重心调节单元通过2个步进电机改变舱内物体的相对位置实现重心的迁移。通过上述2个单元的配合可以实现AUV低功耗的运行模式，按照“之”字轨迹在海洋中缓慢运动。推进器单元由4个独立的推进器组成，通过螺旋桨差速可实现多自由度的运动。多传感器阵列单元包含水声定位模块、声纳成像模块、深度测量传感器、GPS等组成。电力电子单元包含航姿模块，无线电信标模块，电源管理模块等组成，为便于搬运，在重心调节单元舱外壳处焊接有背鳍。AUV结构的基本参数如表1所示。
表1 结构的基本参数
Table 1 Structural parameters
参数总长d1/mm最大直径d2/mm排水体积V/dm3沾湿面积s/m2数值1 71217834.500 40.976 5参数质量m/kg艏部段长度d3/mm艉部段长度d4/mm最大航速v/(m·s-1)数值33.681622404
3 水动力参数计算
AUV在进行作业时易受到诸如海浪、洋流等各种环境因素的影响，其水动力性能则是判定AUV设计好坏的重要依据，也是实现精准控制、精准作业的基础[9]。设计时AUV艏部曲线AB采用双参数椭圆型曲线，艉部曲线CD采用圆锥曲线设计，整体线型简单且加工成本低，如图2所示。其中回转半径D为84 mm，中间段长度L为1 230 mm。
图2 AUV线型示意图
Fig.2 AUV Linetype
为便于计算结果收敛，提升计算速度，在不改变AUV整体外形结构的前提下省略了对水动力性能影响较小但易造成复杂网格的零件，通过建模软件CREO建立了AUV的简化模型，并采用右手惯性坐标系统，以AUV在静态时的浮心作为坐标原点，其中沿原点向艏部为x轴正方向，沿原点向左为y轴正方向，垂直于机身向上为z轴正方向。AUV的水动力模型和坐标系如图3所示。
图3 水动力模型示意图
Fig.3 Hydrodynamic model
3.1 数值计算方法验证
为确保数值计算结果的可靠性，以上海交通大学制作的“海鸥”号[10]为例对文中所用到的计算方法的有效性进行验证，“海鸥”号模型及其网格划分结果如图4所示。求解模型在0°攻角下分别以0.2 m/s、0.3 m/s、0.4 m/s、0.5 m/s的速度直航时的阻力值，并将其与水池拖曳试验结果进行对比，如图5所示。二者数值相差较小，证明了所运用的数值计算方法的有效性。
图4 验证模型与网格划分结果示意图
Fig.4 The verified the model and meshing result
3.2 数值计算与结果分析
AUV在水中的航行速度一般不超过4 m/s，属于低速状态，航行中主要受到水的粘性阻力，且随航行速度的减小其粘性阻力占总阻力的比重将会增大，因此AUV的水动力性能主要体现在水阻力计算与分析方面。借助FLUENT软件对AUV在不同攻角和航速下的阻力系数、升力系数展开了计算。
图5 数值验证对比曲线
Fig.5 Comparison curve of numerical verification
为了更好地模拟低湍流雷诺数流动，避免在k-ε湍流模型中重增加复杂的阻尼函数，采用RANS方法进行数值求解，使用k-ω SST湍流模型。设置流体域入口为速度入口，流体域出口为压力出口，出口压力为0，AUV表面设置为无滑移壁面，四周流域边界设置为自由滑移壁面。将计算域分成两部分，将广阔的海域简化成15L×8L×8L的矩形域，将AUV所在位置设置为球域，浮心即为球域的圆心，直径为4L。针对AUV这种不规则的几何体，主要采用非结构化的四面体网格进行划，为使网格均匀，提高阻力系数的预报精度，将靠近AUV表面流动梯度变化较大的边界层区域进行加密处理。边界层网格按照几何级数膨胀的形式，采用棱柱形网格，其中棱柱层层数为10层，层厚设置为0.2 mm，第一层网格高度设置为0.05 mm, y+值为6。全域网格共有105万，单元质量平均值为0.83，相对质量较高。网格划分结果如图6所示。
图6 网格划分结果示意图
Fig.6 Meshing result
数值计算中攻角α范围为-20°～20°，速度v范围 0.5～4 m/s，CFD求解器采用二阶迎风差分格式，残差设定为10-4。通过105组仿真并将数据汇总得到阻力系数Cd随攻角α变化曲线以及升阻比S随攻角α变化曲线，如图7、图8。
由图7可以看出阻力系数Cd随攻角和航速的增大而增大。当攻角α=0°时，AUV的迎流面积最小，所以此时阻力最小，迎流面积随着攻角的增大也相应增大，AUV的阻力同样随之增大。由于AUV并不是完全对称的结构，故阻力系数Cd的变化趋势在α<0°和α>0°时并不是完全对称的。
图7 阻力系数随攻角变化曲线
Fig.7 Variation curve of drag coefficient with angle of attack
图8 升阻比曲线
Fig.8 Lift drag ratio curve
从图8中可以看出升阻比随攻角的增大而增大，但是升阻比最大不超过0.9，即AUV需要消耗更多的能源来维持运动，因此需要对其外形进行优化。从升阻比曲线变化趋势中可以看出当攻角α<-15°和α>15°时曲线增幅变缓，考虑到图7中阻力系数的变化趋势，选择4 m/s、α=15°为后续优化节点，此时阻力系数Cd=0.281 34，升力系数Cl=0.190 58。当攻角α=±15°时，AUV的压力云图如图9。AUV的外形有一定的对称性，故其表面的压力分布状况相似，由云图可知AUV表面出现了局部的压力中心，这导致了诱导阻力的产生，限制了升阻比的提升。
4 优化分析
直接采用CFD方法进行优化迭代计算量过大，耗时较长，因此考虑构建AUV外形结构参数至水动力参数的响应面[11]，通过目标函数计算、算法寻优获得最优的水动力外形，优化流程如图10所示。优化设计的目的:调整壳体几何尺寸，提高升阻比，改善水动力性能，同时降低AUV的质量。该多目标优化问题的数学描述如下：
目标函数：max Sx(x)，min Mx(x)
约束条件：ai≤xi≤bi
其中：设计变量X=[xi,x2,x3,x4]为壳体主要几何参数 [d3,d4,D,L]；目标函数Sx和Mx分别为升阻比和结构质量；ai和bi为各几何参数的取值上下限，即初始参数的10%。
航行器采用了鱼雷状的线型设计，故艏部段长度d3和艉部段长度d4以及回转半径D的大小直接影响了航行器艏部和艉部的形状，中间段长度L是各个功能舱长度之和，对航行器的长纵比起主要影响。
图9 压力云图
Fig.9 Pressure nephogram
图10 优化流程框图
Fig.10 Optimization flow chart
4.1 优化设计方法
响应面分析法是解决多变量问题的一种统计方法。为明确设计变量与响应值之间的关系，本文中采用拉丁超立方采样方法(LHD)[12]，采用等概率随机正交分布的原则布置试验样本点，以便采用较少的样本点得到较高精度的响应面模型。
响应面拟合是在试验样本点的基础上，通过建立某些函数关系来拟合出设计变量与响应值之间的关系[13]。选择二阶多项式进行拟合，数学表达式为：
(1)
式中：
为响应面的预测值；xi为n维自变量的第i个分量； β0、 βi、 βii、 βij为回归系数。
为评价响应面的预测能力，使用决定系数R2和均方根差σ对响应面的拟合精度进行评估[14]，计算公式：
(2)
(3)
式中: yi为设计点的观测值，当R2趋近于1，σ趋近于0时，响应面的拟合精度则越高。
响应面的拟合优度曲线如图11所示。
图11 拟合优度曲线
Fig.11 Goodness of fit diagram
从图11中可以看出，不同目标变量的响应面预测值与其观测值是线性变化的关系，即响应面的拟合精度较高，拟合精度评估结果如表2所示。目标变量的决定系数均大于0.97，均方根误差均远小于0.1，因此响应面的拟合精度是符合要求的。设计变量d3、d4、D、L对升力系数Cl、阻力系数Cd、质量M和升阻比S的相对重要程度由图12给出，中间段长度L较其他3个设计变量对Cl、Cd、M和S的影响均是最大的；艏部段长度d3对升、阻力系数的影响较艉部段长度d4大，且艉部段长度d4对升、阻力系数的作用是相反的，同时升、阻力系数受壳体回转半径D的影响较艉部段长度d4大，且回转半径D对升、阻力系数的作用也是相反的；回转半径D对升阻比S的影响较艏部段长度d3大，而艏部段长度d3较艉部段长度d4大；对于质量M的影响：艏部段长度d3较艉部段长度d4大，而d4较回转半径D大。
表2 拟合精度评估结果
Table 2 Fitting accuracy evaluation results
目标变量R2σS0.9780.064M0.9980.004
图12 灵敏度直方图
Fig.12 Sensitivity diagram
图13、图14分别给出了升阻比S、质量M随艏部段长度d3、艉部段长度d4、壳体回转半径D、中间段长度L的变化趋势曲面。由图13可以看出，升阻比随着艏部段长度d3的增加而增加，而随着艉部段长度d4的增加而降低，因此要使升阻比尽可能达到最大，d3、d4则应在接近区间端点的位置取值。同时升阻比会随着中间段长度L、回转半径D的变化而产生显著变化，L、D同样应该在接近区间端点的位置处取值，且在区间右端点处取值能够实现升阻比的最大化，当4个参数共同作用时，中间段长度L将对升阻比的变化起着主要作用。同理，由图14可以看出，质量M与4个变量呈线性变化的趋势，且受中间段长度L影响显著，4个参变量均应该在区间左端点处取值。
通过以上分析可以看出4个设计变量之间存在着互斥的关系，并考虑到实际加工制造和内部装配空间上的一些限制，升阻比S和质量M无法同时达到最优值，因此需要进行算法寻优，将所有可能解构成Pareto最优解集[16]，它是一个以牺牲其他目标变量为代价来提高某一个目标变量的解的集合。
图13 输出为S时的响应面
Fig.13 Response surface when output is S
图14 输出为M时的响应面
Fig.14 Response surface when output is M
4.2 多目标遗传算法与结果分析
为了找出目标变量之间的最优平衡点，采用偏移哈默斯利抽样(SHS)技术选取400个样本点，通过多目标遗传算法对优化设计问题进行寻优。MOGA算法的其他参数设置如表3所示。计算升阻比S和质量M是否满足终止条件，若满足则产生终止种群，若不满足则通过适值分配等过程得到下一代种群，继续判断直到满足终止条件[18]。通过多次迭代得到如图15所示的Pareto最优解集分布，圆点的颜色代表结果的优劣，深蓝色为最优，黄色为最劣，其中红色椭圆区域为最优边界，表4为择优选出的3个Pareto最优解。
选取数值计算中升阻比最大，且质量相对较小的的优化点3为最终的设计结果，此时的阻力系数Cd=0.272 53，升力系数Cl=0.214 52，优化设计前后的详细数据如表5所示。优化后的升阻比S和质量M都优于初始值，表明优化后的AUV外形的水动力性能优于原外形。
图15 Pareto最优解集分布图
Fig.15 Distribution of Pareto solution sets
表3 MOGA参数
Table 3 MOGA parameters
最大允许Pareto百分比/%收敛稳定百分比/%变异率交叉率701000.10.9
表4 Pareto最优解
Table 4 Pareto optimal solutions
名称SM/kgd3/mmd4/mmD/mmL/mm优化点10.78524.31170.6230.186.41 211.7优化点20.78624.26172.5228.385.91 209.8优化点30.78824.25173.1227.585.51 205.2
表5 优化前后数据
Table 5 Comparison of data before and after optimization
名称SM/kgd3/mmd4/mmD/mmL/mm优化前0.67826.30162.0240.084.01 230.0优化后0.78824.25173.1227.585.51 205.2
为了验证优化结果的准确性，将优化点中的设计变量的数据传入FLUENT重新进行数值计算[19]，其中，升阻比S的验证结果如表6所示。从表中可以看出优化给出的结果相对于数值计算的结果最大相对误差为1.15%，误差较小，表明优化结果可靠。最终制作出来的实际外形如图16所示。
表6 升阻比的验证结果
Table 6 Verification results of lift drag ratio
名称优化值计算值相对误差/%优化点10.7850.7761.15优化点20.7860.7781.02优化点30.7880.7791.14
图16 实际外形图
Fig.16 Actual outline drawing
5 结论
以混合驱动AUV为研究对象，通过CFD数值计算的方式分析了阻力系数和升阻比的变化情况，基于响应面分析法对AUV外形水动力性能的提升展开优化设计，通过正交试验设计、响应面拟合、灵敏度分析及MOGA算法对外形参数进行优化研究，结论如下：
1) 在各参数取值范围内，同等长度条件下，取较大回转半径有利于提高升阻比，但也易使结构质量增加。AUV的长度越大，阻力系数增幅明显，不利于升阻比的提高和结构质量的减轻，因此应适当减小。在同等半径和长度条件下，增加艏部段长度或减小艉部段长度利于提高升阻比。
2) 基于MOGA算法搜索了全局最优解，得到了AUV外形的优化方案，并通过CFD方法验证了优化结果的正确性。优化后的外形较优化前相比，升阻比由0.678提高到0.788，壳体质量由26.3 kg降低到24.25 kg。该优化设计方法具有良好的寻优能力，可为AUV水动力外形设计提供参考。
[1] 王树新,刘方,邵帅,等.混合驱动水下滑翔机动力学建模与海试研究[J].机械工程学报,2014,50(02):19-27.
Wang S X,Liu F,Shao S,et al.Dynamic modeling of hybrid underwater glider based on the theory of differential geometry and sea trails[J].Chinese Journal of Engineering Design,2014,50(02):19-27.
[2] 李永成,马峥,王小庆.水下滑翔机高效滑翔水动力性能研究[J].中国造船,2020,61(04):52-59.
Li Y C,Ma Z,Wang X Q,Hydrodynamic performance of autonomous underwater glider with high efficiency[J].Shipbuilding of China,2020,61(04):52-59.
[3] 吴利红,俞建成,封锡盛.水下滑翔机器人水动力研究与运动分析[J].船舶工程,2006(01):12-16.
Wu L H,Yu J C,Feng X S,Hydrodynamic research and motion analysis of aug[J].Ship Engineering,2006(01):12-16.
[4] 赵凯凯,李斌,张国伟,等.水下滑翔蛇形机器人滑翔状态水动力参数计算及外形优化[J].高技术通讯,2019,28(08):785-790.
Zhao K K,Li B,Zhang G W,et al.Calculation of hydrodynamic parameters and shell optimization of gliding state of underwater gliding snake-like robot[J].Chinese High Technology Letters,2019,28(08):785-790.
[5] 李天博,王鹏,孙斌,等.一种联翼式水下滑翔机外形优化设计方法[J].哈尔滨工业大学学报,2019,51(04):26-32.
Li T B,Wang P,Sun B,et al.A shape optimization design method of the joined-wing underwater glider[J].Journal of Harbin institute of Technology,2019,51(04):26-32.
[6] 楚灯旺,周悦,龚涛,等.浅海域系泊剖面浮标及其外形优化设计[J].机械设计与研究,2020,36(05):169-173.
Chu D W,Zhou Y,Gong T,et al.Mechanical design and shape optimization design of the mooring profile buoy in shallow waters[J].Machine Design and Research,2020,36(05):169-173.
[7] 张代雨,张贝,王志东,等.翼身融合水下滑翔机水动力及外形优化设计[J].船舶工程,2021,43(03):117-123.
Zhang D Y,Zhang B,Wang Z D,et al.Optimal design of hydrodynamic and shape of blended-wing-body underwater glider[J].Ship Engineering,2021,43(03):117-123.
[8] 高启升,朱兴华,于延凯,等.UUV耐压结构多目标优化设计[J].工程设计学报,2020,27(02):232-238.
Gao Q S,Zhu X H,Yu Y K,et al.Multi-objective optimal design of uuv pressure structure[J].Chinese Journal of Engineering Design,2020,27(02):232-238.
[9] 马骉,张宏伟,王延辉,等.混合驱动水下滑翔机动力学特性与效率分析[J].机械设计,2017,34(08):26-31.
Ma Y,Zhang H W,Wang Y H,et al.Dynamic characteristics and efficiency analysis of hybrid underwater glider[J].Journal of Machine Design,2017,34(08):26-31.
[10] 杨磊.“海鸥号”水下滑翔机设计及其水动力性能研究[D].上海：上海交通大学,2017.
Yang L.Research on design andhydrodynamic performance of“seagull” underwater glider[D].Shanghai：Shanghai Jiao Tong University,2017.
[11] 陈亮,刘荣忠,郭锐,等.扭曲尾翼弹箭气动外形多目标优化[J].兵工学报,2016,37(07):1187-1193.
Chen L,Liu R Z,Guo R,et al.Multi-objective optimization on aerodynamic shape of projectile with twisted empennages[J].Acta Armamentarii,2016,37(07):1187-1193.
[12] Pei J,Wang W J,Yuan S Q,et al.Optimization on the impeller of a low-specific-speed centrifugal pump for hydraulic performance improvement[J].Chinese Journal of Mechanical Engineering,2016,29(05):992-1002.
[13] Zhu H G,Liu L,Long T,et al.Global optimization method using sle and adaptive rbf based on fuzzy clustering[J].Chinese Journal of Mechanical Engineering,2012,25(04):768-775.
[14] 郑昱,孟凡伟,杨占立,等.基于试验设计方法和响应面方法的九索并联机构优化设计[J].机械工程学报,2017,53(17):92-102.
Zheng Y,Meng F W,Yang Z L,et al.Optimization design of a parallel mechanism driven by nine cables based on experimental design methods and response surface methods[J].Journal of Mechanical Engineering,2017,53 (17):92-102.
[15] 张富毅,吴钦,赵晓阳,等.基于响应面方法的喷水推进器进水流道多目标优化[J].兵工学报,2020,41(10):2071-2080.
Zhang F Y,Wu Q,Zhao X Y,et al.Multi-objective optimization of inlet duct of water-jet propulsion based on response surface method[J].ActaArmamentarii,2020,41(10):2071-2080.
[16] Jackson I.Convergence properties of radial basis functions[J].Constructive Approximation,1998,4(01):243-264.
[17] 钱龙,常思江,倪旖.旋转稳定弹扰流片气动外形多目标优化设计[J].兵工学报,2021,42(12):2575-2585.
Qian L,Chang S J,Ni Y.Multi-objective optimization of aerodynamic shape of microspoiler for spin-stabilized projectile[J].ActaArmamentarii,2021,42(12):2575-2585.
[18] 杨卓懿,宋磊.基于MOGA的潜器快速性和操纵性综合优化研究[J].船舶力学,2015,19(05):526-533.
Yang Z Y,Song L.Integrated design of resistance and manoeuvring performances for submersible based on moga[J].Journal of Ship Mechanics,2015,19(05):526-533.
[19] 万云发,孙文磊,王宏伟,等.基于Kriging模型与MOGA算法的风力机主轴轻量化设计[J].太阳能学报,2022,43(03):388-395.
Wan Y F,Sun W L,Wang H W,et al.Lightweight design of wind turbin’s main shaft based on krining model and moga algorithm[J].ActaEnergiae Solaris Sinica,2022,43(03):388-395.