计算条件限制下船舶绕流场数值预报精度解决措施探讨

黄卫刚，邱辽原，姜治芳，陆 超，李 鹏

中国舰船研究设计中心，湖北武汉 430064

0 引 言

船舶绕流场数值预报是当前船舶计算流体力学研究的一个重点，该技术正发展成一种具备较高精度、可用于相似船舶航行性能比较的技术手段。采用CFD技术预报船舶流场存在的难点是：船舶运动时的流体动力会导致船舶出现升沉以及纵倾等运动；高质量的网格划分以及网格差异对计算的影响；数值预报与模型试验之间的误差问题。

船舶运动时出现的升沉以及纵倾等变化是导致数值预报结果与模型试验结果之间出现误差的主要原因之一，且这种误差随着航速的增加而越发明显，如何处理两者之间的差异成为数值预报中需首先解决的问题。国内外对此已开展了各种研究，提出了诸如采用动网格等技术来解决该问题的措施。但在实际工程应用中，动网格技术对计算条件的要求较高。一般情况下，要采用动网格开展船舶绕流数值模拟，需并行几十个CPU开展计算，整个计算周期较长，计算过程中的网格变形有可能会导致网格质量变差，尤其是当计算模型的网格越复杂时，这种情况出现的可能性就越大，而网格质量变差有可能会影响预报精度。如何在只有少量CPU且不采用动网格等技术的前提下实现较高精度的船舶绕流数值预报是本文要解决的主要问题。

数值模拟时的环境条件：计算机为4核CPU、内存为8 G、数值模拟软件的授权数为4个。本文的研究结论均是以此计算限制条件为基础而开展。首先，以KCS和KVLCC为对象标定计算模型，然后以DTMB 5415为对象开展网格类型及规模对计算的影响分析，提出在计算条件受限的条件下宜采用六面体网格开展船舶绕流数值预报的建议。为进一步提高预报精度，本文提出了一种预先给定升沉的船舶绕流数值预报技术，并将以某条船为对象对该技术进行验证性计算。最后，以两条船为对象开展数值分析，进一步验证CFD技术在船型优化中的效用。

1 国内外研究情况

模型尺度的船舶绕流模拟正在不断进步，随着技术的进步，船舶计算对象考虑的也将越来越真实。

在自由液面模拟方面，自由液面跟踪方法非常有效。Gorski［1］以一条穿浪内倾船为对象，采用自由液面捕捉的方式进行了数值模拟，数值模拟确定的沿船体表面及船尾的兴波与模型试验测定的兴波之间吻合较好。Raven等［2］说明了自由液面捕捉技术仍有提升的空间，通过消除方程中的时间项并严格解决稳态自由液面粘性流的方式，可以通过较少的数值计算来确定兴波，该兴波与系列60以及Dyne水池的模型试验结果之间的吻合程度较高。在2007年的国际数值船舶水动力学会议上，西班牙Polytechnic大学的研究者指出，采用六面体网格计算船舶绕流能取得较高的精度。

当傅汝德数较低时，因船舶运动引起的升沉以及纵倾较小，因而可以直接采用不计及升沉及纵倾的方法来计算船舶阻力，此时，模型试验的阻力与数值预报的阻力之间误差较小。当傅汝德数增加时，船舶运动产生的升沉以及纵倾逐渐明显，此时，若还采用不计及升沉及纵倾的方法计算阻力，则该阻力与模型试验的阻力间的误差就会逐渐明显，在Fn较高时可能会引起约10%甚至更高的差异。因此，在当前阶段，许多船舶粘性流数值模拟仍采用预先设定升沉和纵倾的方式。例如，Azcueta［3］阐述了采用VOF方法可以高精度预报一条单体濒海战斗舰在 Fn=0.25～0.9时的阻力情况；Azcueta［4］证明了一条高速艇在 Fn=0.75～4.0间的相似结论。

2 计算原理

1）应用RANS方程和连续性方程解决湍流问题时，由于增加了雷诺应力项而使基本方程不封闭，因而需要根据湍流的运动规律构建合适的湍流模型以解决方程的封闭性问题。为了更好地处理高应变率及流线弯曲度比较大的流动，本项目湍流模型采用了SST模型。

2）由于自由液面需要考虑空气和水的相互作用，因而属多相流问题。本研究处理自由液面问题时使用的是VOF方法。这种方法是一种以流体占据网格单元体积份额的途径来跟踪自由面演化的方法，可以跟踪发生复杂变形的自由面。因VOF方法的流体质量守恒在方程的构建中就得到了满足，容易扩展到三维，体积分数的演变只需要用到相邻的单元，因此得到广泛应用。

3 计算条件受限情况下排水型中、高速船流场预报存在的问题及解决措施

3.1 存在的问题

在计算条件受限的情况下，排水型中、高速船流场预报存在的主要问题有：

1）船舶运动时的流体动力导致船舶出现升沉以及纵倾等运动。

排水型高速船与低速船在航行过程中存在明显差异。排水型高速船在运动过程中的流体动力较明显，导致船舶在运动过程中存在较大的纵倾和升沉，而低速船的升沉和纵倾则相对较小。很显然，在固定吃水以及升沉的前提下，船舶在运动过程中形成的升沉和纵倾会成为影响数值预报结果与模型试验结果间差异的主要因素，且会随着傅汝德数的增加而更加明显。

表1所示为一条缩尺比为1∶22的高速集装箱船船模随着傅汝德数的增加其升沉和纵倾的变化情况。在实船状态，当Fn=0.21时，该船的升沉约为80 mm，在 Fn=0.39时，升沉约为360 mm。可见随着傅汝德数的增加，升沉和纵倾的变化也会比较明显，并将成为影响数值预报结果与模型试验结果间误差的关键因素。

表1 某船模在不同傅汝德数下的升沉和纵倾变化Tab.1 The sinkage and trim of a ship at different Fn

2）高质量的网格划分以及网格差异对计算的影响。

3）数值预报与模型试验之间的误差问题。

3.2 解决措施研究

为了解决在计算条件受限情况下排水型中、高速船流场预报存在的问题，本文提出了下列解决措施：

1）在软、硬件条件受限的情况下，建议直接采用六面体网格，且网格质量应满足Angle＞20°，Aspect ratio＜100，Maximum Mesh Expansion Factor≤20等要求。该措施能保证在计算条件受限时数值预报结果具有较高的精度，减少对计算条件的要求。

2）如果设计船舶的傅汝德数较小，可以直接采用不计及升沉和纵倾的方法来计算裸船体阻力。如果设计船舶的傅汝德数较大，则建议采用预先给定升沉的方法，该措施能在避免使用动网格等技术的前提下解决固定吃水时船模所受浮力和重力不平衡的问题，有效提高计算精度。

为了验证上述措施的可行性，本文以DTMB5415，KCS和KVLCC等船模为对象开展了数值预报，DTMB5415，KCS及KVLCC这三型标模的尺度以及三维模型如表2和图1～图3所示。

表2 标模尺度简表Tab.2 Standard model’s dimension

图1 DTMB 5415模型Fig.1 The model of DTMB5415

图2 KCS模型Fig.2 The model of KCS

图3 KVLCC模型Fig.3 The model of KVLCC

3.2.1 计算模型标定

计算模型满足下列条件：

1）网格为六面体网格，Angle＞20°，Aspect ratio ＜100，Maximum Mesh Expansion Factor≤20。

2）湍流模型为SST模型。

3）流体域范围：为避免远方边界条件对近船体流场的干扰，在文献［5-10］的基础上，计算区域入口取为船体艏部向上游延伸至2倍船长处；出口取为艉部向下游延伸至2倍船长处；区域外边界分别为由对称面（船体纵中剖面）向右舷方向和设计水线面向下延伸约0.8倍船长处；区域上边界取为设计水线面向上约2倍吃水高度处。

4）边界条件：入口为Inlet；出口为Opening；主船体壁面为No Slip Wall；流体域的底面、顶面及侧面为Free Slip Wall；流体域的中线面为SYMMETRY。

在满足上述条件的前提下，对标模KCS和KVLCC开展了数值预报。表3和表4所示为数值预报的裸船体阻力与模型试验的裸船体阻力间的误差对比。

表3 KCS模型试验和数值预报的阻力对比Tab.3 The resistance’s difference between KCS’s model test and CFD

表4 KVLCC模型试验和数值预报的阻力对比Tab.4 The resistance’s difference between KVLCC’s model test and CFD

表3和表4表明，在限制条件下，采用六面体网格预报的裸船体阻力与模型试验阻力之间的误差可以控制在工程精度范围内。这验证了各种文献所报道的结论，即在航速较低时，采用六面体网格预报的流场与模型试验之间的吻合程度较好。

3.2.2 网格类型及规模对计算的影响

为进一步分析网格差异对计算的影响，针对DTMB 5415模型，分别对流体域为简单四面体网格、四面体网格+边界层网格以及六面体网格的单流体域开展了数值预报。受内存所限，六面体网格的网格数量为240万，四面体网格的网格数量约为450万，四面体网格+边界层网格的网格数量约为480万。计算的网格如图4所示。

图4 DTMB 5415网格划分情况Fig.4 The mesh of DTMB 5415

计算速度Vm=2.096 m/s，该状态下船模的模型试验阻力为45.1 N。不同的网格对CFD计算裸船体阻力误差的影响的计算结果如表5所示。

表5 不同网格对数值预报阻力的误差影响Tab.5 Influence of different meshes on the ship resistance calculation results

由表5可知，六面体网格的阻力预报结果与模型试验结果之间的误差最小，四面体网格+边界层网格的次之，四面体网格的最大。

为进一步分析造成这种差异的原因，对尾部螺旋桨位置的x方向伴流进行了对比，如图5所示。

由图5可知，在计算条件受限的前提下，与模型试验测得的沿主船体的流场形态相比，六面体网格计算确定的流场形态与模型试验结果最接近，四面体+边界层网格计算确定的流场形态次之，四面体网格计算确定的流场形态与模型试验测得的流场形态间的误差最大。出现该现象的原因为：

图5 网格差异对尾部螺旋桨位置x方向伴流影响Fig.5 Influence of different meshes on the x-direction flow around ship’s propeller

1）在进行数值绕流的边界层计算时，六面体网格沿船体表面分布的网格厚度可以很小，Y+可以小于200甚至是接近于10，能满足湍流模型的需求，有助于精确捕捉边界层的流场。四面体+边界层网格相对于六面体网格而言，其沿船体表面分布的网格能满足湍流模型的需求，但在边界层与流体域之间的连接区域，其网格变化较大，导致离边界层较远区域的流场形态与模型试验间的差异较大。四面体网格相对于六面体网格而言，其沿船体表面分布的网格厚度最大，离满足湍流模型的需求尚存在一定的差距。

2）在开展自由液面计算时，六面体网格在自由液面附近能设置较多的网格，可为捕捉自由液面提供良好的条件。四面体+边界层网格虽然弥补了四面体网格在捕捉边界层附近流场上的缺点，但由于自由液面附近的网格较少，影响了自由液面计算的精度，从而导致计算的自由液面与模型试验确定的自由液面之间的差异较大。四面体网格的自由液面捕捉情况与四面体+边界层网格的情况类似。

3.2.3 预先给定升沉的船舶绕流数值预报技术

如前文所述，傅汝德数的增加会导致数值预报结果与模型试验结果间的误差逐渐加大，因此，本文接下来将要解决减少这种误差的方法。一般情况下，可以采用动网格等技术来解决，但这种方法存在计算时间较长、计算过程中容易出现发散等问题，在计算条件受限的情况下实施难度很大。针对这种情况，这里将介绍一种预先给定升沉的船舶绕流数值预报方法。

根据计算排水量状态下对应的吃水来计算船舶阻力。在正常情况下，计算1000～2000步时，计算的垂向力（即船舶所受浮力）就会出现如图6所示的情况。这时，可以根据船体所受垂向力与模型重量间的差值来重新确定吃水，然后，根据新的吃水重新开展计算，直至船体垂向力与模型重量间的差异小于0.5%为止。此时，就可以认为数值预报时的排水量状态与模型试验的排水量状态是一致的，一般迭代2～3次就可以得到计入升沉影响的阻力。

现以某型集装箱船为例分析计入升沉与不计入升沉对数值预报结果的影响。该船模型重量为14400 N，当Fn=0.28时，如果采用不计入升沉的船舶绕流数值方法开展计算，此时计算确定的浮力为14120 N，比模型重量小了约280 N。如果采用计入升沉的船舶绕流数值方法开展计算，计入升沉后的浮力为14400 N，与模型的重量一致。当计入升沉的影响后，船舶的吃水增加，所受的浮力也随之增加，从而达到重量与浮力的平衡，因此，此时的升沉状态便与模型试验状态基本一致。

同时，还对比了计入升沉与不计入升沉对阻力的影响。当Fn=0.28时，根据模型试验确定的阻力为120 N。计入升沉与不计入升沉后的阻力计算结果如表6所示。

表6 计入升沉与不计入升沉对阻力计算的影响Tab.6 The resistance calculation results with or without sinkage taken into account

由表6可知，不计入升沉时的阻力相对模型试验确定的阻力120 N小了约3.7%，计入升沉后的阻力相对模型试验的阻力小了约1.0%。计入升沉后的排水量更接近于模型试验状态，此时两者之间的阻力误差也更小。计入升沉的影响能有效减少数值预报结果与模型试验结果之间的误差，尤其是当船舶的傅汝德数较大时，这种效果更加明显。当两型存在较大差异的船进行阻力对比时，建议采用文中所述的方法。

4 验证性计算

研究CFD技术的最终目的是应用到船型优化中。为了进一步检验现有计算经验的可靠性，在限制条件下，针对两个集装箱船的相似方案（方案A和方案B）进行了流场计算。方案A和方案B均为长度约8 m的相似模型，船型系数和主尺度基本一致。计算的网格为六面体网格，采用的湍流模型为SST模型。数值预报的阻力对比如表7和表8所示。

由表7、表8可知：

1）在方案A和方案B中，流场计算确定的阻力与模型试验之间的阻力误差均小于5%。

2）两型船的数值预报趋势与模型试验趋势基本一致。

3）数值预报阻力与模型试验阻力的差异符合下列趋势：随着傅汝德数的增加，模型试验阻力与数值预报阻力的误差逐渐加大。

表7 方案A计算阻力和试验结果比较Tab.7 The project A resistance’s difference between model test and CFD

表8 方案B计算阻力和试验结果比较Tab.8 The project B resistance’s difference between model test and CFD

5 结 论

本文探讨了计算条件限制下船舶绕流数值预报精度的解决措施，认为要提高数值预报精度可以采取以下措施：

1）在开展船舶绕流数值预报时，如果其软、硬件环境与本文类似，直接采用六面体网格比较合适。

2）如果所设计船舶的傅汝德数较小，可以直接采用固定吃水和纵倾的方法来计算裸船体阻力。一般情况下，数值预报结果与模型试验结果之间的误差能控制在工程允许的精度范围内。

3）如果所设计船舶的傅汝德数较大，此时船舶的升沉和纵倾已成为影响数值预报与模型试验间差异的主要因素时，建议采用文中提出的预先给定升沉的方法。这种方法虽然会增加一些计算时间，但在计算条件等限制情况下，相对动网格等技术而言更简单、实用。

通过本文的验证性计算还表明，在相同的网格参数以及边界条件设置下，采用数值预报来判别型线变化对阻力以及流场形态影响的方法是可行的。

［1］GORSKI J J.Present state of numerical ship hydrodynamics and validation experiments［J］.Journal of Offshore Mechanics and Arctic Engineering，2002，124：74-80.

［2］RAVEN H C，PLOEG A P，STARKE A R.Computation of free-surface viscous flows at model and full scale by a steady iterative approach［C］//Proceeding 25th Symposium Naval Hydrodynamics.Canada，2004.

［3］AZCUETA R.Steady and unsteady RANSE simulations for planning crafts［C］//The 7th International Conference on Fast Sea Transportation，FAST2003.Italy，2003.

［4］AZCUETA R.Steady and unsteady RANSE simulations for littoral combat ships［C］//Proceeding 25th Symposium on Naval Hydrodynamics.Newfoundland and Labrador，Canada，2004.

［5］王中，卢晓平.基于细分曲面的兴波阻力计算方法［J］.海军工程大学学报，2006，18（增刊）：20-25.WANG Zhong，LU Xiaoping.A calculation method for ship wave making resistance based on subdivision surface theory［J］.Journal of Naval University of Engineering，2006，18（Supplement）：20-25.

［6］PENG Q J，HOU Y J.A VOF based numerical model for breaking waves in surf zone［J］.Chinese Journal of Oceanology and Limnology，2006，24（1）：57-64.

［7］GAO Q X.Numerical simulation of free surface flow around ship hull［J］.Journal of Ship Mechanics，2002，6（3）：1-13.

［8］ZHANG Zhirong，ZHAO Feng，LI Baiqi.Numerical calculation of viscous free-surface flow about ship hull［J］.Journal of Ship Mechanics，2002，6（6）：10-16.

［9］HINO T.A study of grid dependence in navier-stokes solutions for free surface flows around a ship hull［J］.Journal of the Society of Navy Architects of Japan，1994，176：11-18.

［10］SHIOTANI S，KODAMA Y.Numerical computation of viscous flows with free surface around a series 60 model［J］.Journal of the Society of Naval Architects of Japan，1996，180：39-50.