0 引　言

船舶快速性研究的首要目的在于降低航行阻力。从阻力的形成机理出发，可将总阻力分成兴波阻力及粘性阻力两大类。兴波阻力是由于船体扰动流场形成波浪所致，而粘性阻力则可进一步分成摩擦阻力及漩涡阻力，其中摩擦阻力是湿表面上切向力的合力，漩涡阻力是流场前后不对称而产生的压差力^[1]。

对于高速船舶而言，兴波阻力往往可占总阻力的50％以上，且对船型变化相当敏感，局部的型线修正即可能引起阻力曲线的明显波动，故针对兴波阻力的船型优化具有显著的工程价值。本文主要关注舰首的导流罩型线，此处要兼顾基阵布置与舰船减阻的双重需求，且二者往往相互制约。多数文献仅考虑阻力性能，一般得到前伸或上翘形的球艏，但无法满足声呐的安装条件，因此实用性较低，本文通过引入几何约束来避免这一缺陷。

各方面对于优化设计的定义均不尽相同，但其核心思想却一致：即优化设计是在给定约束条件的基础上，通过调整设计变量，反复迭代直至获得最优解的过程。同时，优化设计也是一种充分探索、利用系统中各子系统相互作用的协同机制来设计复杂系统工程的方法论^[2]。现阶段的实践表明，型线优化已从传统的专家优选发展至海量方案的精细化评估，其核心在于数值模拟技术的进步与评估自动化水平的提高，相关进展可参考刘祖源等^[2]、李胜忠^[3]、Diez等^[4]的工作，近期国内外开展了VLCC油船^[5]、21万吨级散货船等低速肥大船型的阻力与推进性能优化研究^[6]和水下潜航器多学科综合优化研究^[7]等一系列工作，同时在型线变换方法^[8]、设计空间采样方法^[9]等支撑技术方面也取得了有益进展。本文基于上述思想，使用自编程序完成了DTMB-5415球艏型线优化，具体流程如图1所示。

1 优化体系构建 1.1 优化关键技术

船型优化的核心由船型几何变换技术、水动力性能评估技术和优化算法技术组成^[10]，其中几何变换是构造样本、探索设计空间的基础，水动力计算是评估方案优劣的依据，而优化算法则是对改进方向的把握和指导，优化程序将反复执行上述过程直至收敛。本节按照上述顺序介绍对应的理论方法，简要概况可参考表1。

船型几何变换技术可分为船型参数法和几何造型法两类。前者利用主尺度、船型系数等宏观特征来反映型线，而后者则利用参数化的样条曲面描述船型，几何造型法包括母型融合法、叠加Bezier曲面法、FFD自由变形法、径向基函数法等。其中融合法通过不同的权重系数，使新船型的几何特征介于多艘母型船之间，其优势在于设计变量少、变换前后光顺性强，但优化区间受到母型船的约束，不能实现全域探索。叠加Bezier曲面法通过在船体表面叠加样条曲面、调整控制点坐标来修改型线，其优势在于调整精细化程度高、设计空间约束度低，特别适用于局部型线的修改工作^[11]。FFD自由变形法是当前运用最为广泛的几何变换技术，适用于全船范围内的整体变形或局部型线调整，且变换后的光顺性仍能得到保证。径向基函数插值法通过部分型值点变形前后的坐标建立映射关系，再推广至全域，易于处理几何约束条件。

水动力评估技术主要包括三类：基于统计回归的经验公式、基于势流理论的数值模拟、基于黏性流理论的数值模拟。其中经验公式分辨率较低，不能反映型线局部变换对阻力的影响，因此应用有限。势流方法包括Mitchell兴波积分、Noblesse新细长体理论、非线性Dawson法等，王中等^[12]曾利用Mitchell兴波积分法计算了某穿浪双体船的阻力，利用非线性Dawson法计算带有艉倾的DTMB-5415船模兴波阻力。D.Peri等^[13]利用三维边界元程序完成了某护卫舰型线优化，证明势流理论在计算兴波阻力方面确有较强的实用性。而粘性流理论则更为全面、真实地描述了流动的本质，可更精确地考虑流体粘性和湍流对流场的影响，例如STAR-CCM+、Fluent等商用求解器已在工程领域应用广泛，但因计算消耗较大，在求解优化问题时多与代理模型或降维技术结合使用。

常见的优化算法包括基于梯度的搜索算法和随机搜索算法两类。前者包括梯度下降算法、二次规划算法等，而后者主要代指遗传算法及粒子群算法等群智能算法。基于梯度的搜索算法包括梯度下降法、共轭梯度法、牛顿法、序列二次规划法等，其优势在于具备严格的数学理论基础、收敛速度快，但全局搜索能力偏弱，因此在面对高维优化、复杂优化时存在一定的局限性。型线优化的最优解往往没有规律可循，特别是多学科优化设计时，目标函数与设计变量之间没有显式函数关系，梯度信息仅能依靠差分法获取。因此若在复杂优化问题中采用梯度算法，则必须面对计算成本大、容易陷入局部最优解等困扰，而随机搜索算法具有不依赖梯度信息、全局搜索能力强等特点，因此得到了广泛的应用^[2]。

1.2 FFD自由变形算法

自编程序以FFD算法作为几何重构手段，NURBS参数曲面为型线描述方法，可确保在任意条件下得到光顺型线并兼顾图形交互需求。FFD自由变形法的核心概念是映射：将笛卡尔坐标转化为参数坐标，再利用某种函数进行运算，从而达到调整坐标的目的。因此算法的关键是找准初始值和映射方式，前者是参数坐标而后者是变形函数。传统的FFD控制点均匀分布在控制体的顶点、棱边、表面及内部。设该控制体的原点为$ {X}_{0} $，边长矢量分别为$ S $、$ T $、$ U $，并将其划分为$ l\times m\times n $个均匀网格，则控制点坐标为：

$ {P}_{i,j,k}={X}_{0}+\frac{i}{l}S+\frac{j}{m}T+\frac{k}{n}U 。$

(1)

设空间内任意一点的坐标为$ \left(x,y,z\right) $，为满足基函数映射要求，先将待求的空间坐标通过表达式(2)作无因次化处理。

$ \left\{\begin{aligned} &s=\frac{x-{x}_{0}}{\left|S\right|}，\\ &t=\frac{y-{y}_{0}}{\left|T\right|}，\\ &u=\frac{z-{z}_{0}}{\left|U\right|}。\end{aligned}\right. $

(2)

空间中任意一点的坐标可根据其无因次坐标和控制点坐标来表示，故改变控制点坐标，再代入映射函数表达式(3)即可得到变形后的坐标，而式(4)给出了Bernstein基函数$ {B}_{i,l} $的计算方法。

$ {X}_{(s,t,u)}=\sum _{i=0}^{l}\sum _{j=0}^{m}\sum _{k=0}^{n}{P}_{i,j,k}{B}_{i,l}\left(s\right){B}_{j,m}\left(t\right){B}_{k,n}\left(u\right) ,$

(3)

$ B_{i,l}(s)=\frac{l!}{l!(l-i)!}s^i (1-s)^{l-i} 。$

(4)

注意到Bernstein基函数之和等同于$ [s+(1-s)]^l= 1 $，这又称为Bernstein基函数的单位分解性。该式为一元基函数，对应张量积形式的多元基函数也具有该性质，这给后续计算带来了很大的好处。此外该基函数还具有非负性和全域支撑性，即在除端点以外的定义域上函数值均大于0。本文使用NURBS方法描述船型曲面，实际变换的是NURBS曲面控制点坐标，相较于传统的船型参数或型值点表述法，参数化曲面允许以更高的精度、更简洁的形式描述型线，同时也便于和外部软件交互数据。

1.3 线性兴波理论

常见的势流方法已在1.1节进行了介绍，程序中使用的是P.R. Couser等^[14]提出的改进型Mitchell线性兴波理论，计算时假定流场无限大，流体无粘无旋并忽略表面张力，且波幅与波长相比为一小量，因而满足线性无破碎条件，于是兴波阻力可通过远场波形间接计算。相比于直接形式的压力积分法，波形法的数值稳定性更好，可以在很大程度上避免因艏艉波形奇异而引发的负阻力问题。取计算域的宽度为$ B $、深度为$ H $，设$ \sigma $为面元中心点处的源强，$ U $为速度，$ S $为面元面积，$ \widehat{n} $为曲面法向量，$ \zeta $为自由面波高，$ \stackrel{-}{k} $表示波数，$ \theta $为波向角，在满足连续性方程、自由液面和物面条件时可推得兴波阻力计算公式(5)，式中所需参数通过式(6)～式(10)计算。

$ \begin{split} {R}_{w}=&\frac{\rho gB}{4}\Bigg\{{{\zeta }_{0}}^{2}\left[1-\frac{2{k}_{0}H}{\mathrm{sin}h\left(2{k}_{0}H\right)}\right]+\\ &\sum _{m=1}^{\infty }{{\zeta }_{m}}^{2}\left[1-\frac{{ \mathrm{cos}}^{2}{\theta }_{m}}{2}\left(1+\frac{2{k}_{m}H}{\mathrm{sin}h\left(2{k}_{m}H\right)}\right)\right]\Bigg\} \end{split}，$

(5)

$ {{\zeta }_{m}}^{2}={{\xi }_{m}}^{2}+{{\eta }_{m}}^{2} ，$

(6)

$ \left|\begin{array}{c}{\xi }_{m}\\ {\eta }_{m}\end{array}\right|=\frac{16\pi U}{Bg}\frac{\stackrel-{k}+{k}_{m}{\mathrm{cos}}^{2}{\theta }_{m}}{1+{\mathrm{sin}}^{2}{\theta }_{m}-\stackrel-{k}H{\mathit{sec}h}^{2}\left({k}_{m}H\right)}D ，$

(7)

$ \begin{split}D= & \displaystyle\sum_{\sigma}^{ }\sigma_{\sigma}e^{-k_mH} \mathrm{cos}h\left[k_m\left(H+z_{\sigma}\right)\right]\times \\ & \left|\begin{array}{c}\mathit{\mathrm{cos}}\left(k_mx_{\sigma}\mathit{\mathrm{cos}}\theta_m\right) \\ \mathrm{\mathit{\mathrm{sin}}}\left(k_mx_{\sigma}\mathit{\mathrm{cos}}\theta_m\right)\end{array}\right|\left\{\begin{array}{c}\mathit{\mathrm{cos}}\left(m{\text{π}} y_{\sigma}/B\right) \\ \mathit{\mathrm{sin}}\left(m{\text{π}} y_{\sigma}/B\right)，\end{array}\right.\end{split} $

(8)

$ \sigma=-\frac{\hat{n}\cdot U\times S}{2\boldsymbol{\text{π}}}，$

(9)

$ \overline{k}=\frac{g}{{U}^{2}} 。$

(10)

表达式(5)可直接用于计算Wigly、Series-60等船型的兴波阻力，但在求解DTMB-5415船型时必须计入方尾修正，以反映虚长度造成的影响。Couser在文献中介绍了压力修正、附加源汇修正、虚长度修正3种方法。其中压力修正法假设船舶的航速足够快，艉板完全暴露在空气中，因此在艉封板上对静水压作积分即可。实践证明该方法不适用于中低航速工况，一是此时艉封板浸水面积不确定，二是线性理论对动压的求解精度不高，即便在艉板湿表面上作压力积分也难以获得可靠结果。附加源汇修正法是在艉封板中心线或底边界上布置额外的点汇，以弥补式(5)的缺失部分，该方法在中高航速时的适用性较好。虚长度修正法则是将艉板至鸡尾流之间的区域视作船体的一部分，将船体延伸并光顺过渡直至封闭，随后直接使用式(5)进行计算。虚长度参数需要通过船模试验或经验公式确定，但一般情况下，只要虚长度区域建模比较准确，兴波阻力的计算精度就有保证。本文采用虚长度法进行方尾修正，设艉封板宽度为$ h $，虚长度取$ \Delta L=3h $，过渡区域的型线可通过插值得到。势流计算域示意图见图2。

1.4 优化粒子群算法

传统的粒子群算法（PSO）是Kennedy和Eberhart等参考鸟类觅食行为而发展起来的一种基于群体协作的随机搜索算法^[2]。该方法因具有原理简单、通用性较强等特点而得到广泛应用。研究人员在原始方案的基础上引入了各类修正模型，以进一步提高该方法的求解效果。主要的改进思路是从优化参数取值、优化寻优方法两方面着手，例如，初始位置由正交表或拉丁方采样确定、初始速度正交化、引入惯性权重参数，或将粒子群算法与其他群智能算法相结合等，最终形成了一系列优化粒子群算法（Improved Partical Swarm Optimization）^[3]。

本文IPSO算法的改进之处有两点，一是在初始化阶段利用正交表采样代替随机选点，使初始样本在设计空间内均匀分布，二是引入惯性权重，以平衡局部与全局搜索能力。算法在N维优化空间中选取m个样本粒子，则某粒子的位置、速度、粒子历史最优解、种群最优解等信息的计算式：

$ \left\{\begin{aligned} &{X}_{i}=\left[{X}_{i1},{X}_{i2},{X}_{i3},{X}_{i4}\dots \dots {X}_{iN}\right]，\\ &{V}_{i}=\left[{V}_{i1},{V}_{i2},{V}_{i3},{V}_{i4}\dots \dots {V}_{iN}\right]，\\ &{P}_{i}=\left[{P}_{i1},{P}_{i2},{P}_{i3},{P}_{i4}\dots \dots {P}_{iN}\right]，\\ &{P}_{g}=\left[{P}_{g1},{P}_{g2},{P}_{g3},{P}_{g4}\dots \dots {P}_{gN}\right]。\end{aligned}\right. $

(11)

设$ n $为迭代次数，每次迭代后使用式(12)、式(14)更新粒子的速度与坐标，其中动态惯性权重由式(13)计算，该系数通过改变粒子对优化方向的偏好来影响搜索性质，例如出现收敛趋势时权重系数增加，使得粒子保持自身运动方向，利于扩大搜索范围，反之呈现无序状态时则减小权重系数、促使粒子向优解靠拢，以加速收敛。

$ {V}_{i} \left(n + 1\right) = w{V}_{i}\left(n\right) + {c}_{1}{r}_{1}\left({P}_{i} - {X}_{i}\left(n\right)\right) + {c}_{2}{r}_{2}\left({P}_{g} - {X}_{i}\left(n\right)\right) ,$

(12)

$ w=\left\{\begin{aligned} & w\mathrm{_{min}}\frac{\left(w\mathrm{_{max}}-w_{\mathrm{min}}\right)(f-f\mathrm{_{min}})}{\overline{f}-f\mathrm{_{min}}}+n^2\frac{w\mathrm{_{max}}-w\mathrm{_{min}}}{n_{\mathrm{max}}^2} ,\\ & w_{\mathrm{max}},\end{aligned}\right. $

(13)

$ {X}_{i}\left(n+1\right)={X}_{i}\left(n\right)+{V}_{i}\left(n\right) 。$

(14)

式中：惯性权重上下限分别取0.9和0.3；$ {c}_{1} $和$ {c}_{2} $为学习因子，一般均取常数2；$ {r}_{1} $和$ {r}_{2} $为介于0～1之间的随机数。此外要设定最大速度$ V\mathrm{\mathrm{_{max}}\mathrm{\mathrm{\mathrm{\mathrm{\mathrm{\mathrm{ }}}}}}} $以防错过最优解，$ V\mathrm{_{max}} $一般取区间大小的0.1～0.2倍，当粒子速度足够小或达到设定迭代次数后结束优化。

2 艏部型线优化过程 2.1 目标函数与设计变量

以降低单位排水量下的兴波阻力为优化目标，构造目标函数为：

$ G\left(Fr\right)=\frac{Res\left(Fr\right)}{\Delta } ，$

(15)

$ F=\sum _{i=1}^{M}{w}_{i}\frac{G\left({Fr}_{i}\right)}{{G}_{0}\left({Fr}_{i}\right)} 。$

(16)

式中：$ M $为目标工况数量；$ G $为单位排水量下的兴波阻力；$ {G}_{0} $为原型方案阻力性能；$ Fr $为傅汝德数；$ w $为权重系数。计算中实取Fr=0.28、0.37、0.44等工况，对应实船20、27、32 kn航速，以便获得在较广航速范围内均能减阻的有效改进方案。

根据相关文献经验，导流罩的阻力性能主要受最大横剖面面积、最大横剖面纵向坐标和沉深3个因素影响^[15]，因此本文将优化的重点放在导流罩主尺度和横剖面形状上。优化过程中共设定20个设计变量，包括5个FFD变量和15个NURBS控制点变量。其中FFD变量主要从宏观上控制球艏几何特征，除1个用于调节球艏前伸距离、1个调节球艏吃水深度外，余下用于控制球艏沿纵向的宽度变化，15个控制点变量分布于导流罩横剖面上，用于调节横剖面局部型线。

2.2 优化约束条件

优化方案需同时满足船型参数约束和局部几何约束条件，由于本文只涉及水线以下球艏区域的变形，因此水线面形状及主尺度固定不变，湿表面积变化很小，可认为已自动满足船型参数约束，故只需给定局部约束条件。从总布置角度考虑，导流罩首先要能容纳基阵并与主船体光顺过渡，其次应尽量前伸并增加沉深，以远离自由液面，但过分前伸不利于艏锚和边锚布置，同时吃水深度过大对舰艇航行性能也有影响^{[15, 16]}，所以导流罩的主尺度和基本参数调节范围比较有限，反而是横剖面形状具有很高的优化自由度，设水线长为$ L $、型宽为$ B $、吃水为$ D $，具体限制条件参考见表2。

2.3 优化结果分析

优化方案的型线特征与文献[15]的研究结果一致，即球艏前伸、吃水加深，横剖面由椭圆形变为梨形，球艏与主船体的过渡更为缓和，导流罩容积增加而湿面积变化不大。图3～图5所示为改进前后的型线图、兴波阻力曲线和纵切波形，其中$ x $是以艏垂线为基准的纵向坐标，$ y $是横向坐标，$ \zeta $为波幅，且三者均作无因次化处理。计算表明，在排水量增加1.5%，湿面积增加0.9%的情况下，Fr=0.28、0.37、0.44工况下的兴波阻力分别下降15.3%、9.9%、7.1%。对比优化前后的纵切波形，可见远场兴波高度均有所下降，特别是艏波在相位和幅值上均发生了明显变化。

球艏的流动机理相当复杂，在此只以波形分析法作定性探讨。考虑到波浪翻卷破碎与气泡卷吸等非线性流动仅发生在舷侧附近，因此在远场区域可假设全船的兴波是由主船体波系和球艏波系线性叠加而来，则以一艘无球艏的DTMB-5415波形为基准，可将球艏波系分离开来。在此以Fr=0.28工况为例进行对比（见图6），可见球艏前伸使得波系相位提前，而宽度和剖面变化所引起的横向扩张增加了球艏的体积，使其对自由面的扰动更为明显，增大了波系幅值，从而更好地抑制艏部兴波。从计算结果来看，主船体的艏波与球艏兴波基本在同一量级，因此通过优化球艏型线来降低兴波阻力完全可行，倘若放宽约束条件以允许导流罩继续前伸或接近自由液面，则有望收到更好的效果。

3 结　语

本文通过自编程序完成了DTMB-5415球艏型线优化，在满足导流罩实用性的前提下将Fr=0.28、0.37、0.44时的兴波阻力分别降低15.3%、9.9%、7.1%，对应排水量增加1.5%、湿面积增加0.9%，研究表明：

1）线性势流理论的计算结果虽与实验值仍有一定差距，但该方法确能有效反映型线与兴波阻力之间的关系，可以作为初步设计阶段的型线优化工具。

2）基于线性假设的波形分析法可直观地反映球艏波系与主船体波系之间的关系，从而为型线优化提供定性指导。

3）导流罩适度前伸、采用梨形剖面有利于减少兴波阻力，该形式下的球艏波系在相位和幅值上与主船体兴波相匹配，消波效果更为突出。

应该指出的是，目前的研究仅限于静水阻力性能而未考虑流噪声因素，因此改进方案的实用性仍有待提高。如何计算近场区域的非线性流动，并将流噪声纳入优化框架，将是后期的研究重点。