邓磊，彭弘宇

1中国人民解放军92692部队装备部，广东湛江5240642海军驻上海江南造船（集团）有限责任公司军事代表室，上海201913

小水线面双体船耐波性能CFD不确定度分析

邓磊1，彭弘宇2

1中国人民解放军92692部队装备部，广东湛江524064
2海军驻上海江南造船（集团）有限责任公司军事代表室，上海201913

基于ITTC推荐的CFD不确定度分析规程，开展基于RANS方程的小水线面双体船（SWATH）波浪中纵向运动计算结果的不确定度分析。分别验证计算模型的网格收敛性和时间步长的收敛性，发现相比于网格的疏密，该数值求解模型对于时间步长更为敏感。通过与试验数据进行对比，计算结果基本得到确认，并由此建立了顶浪规则波中SWATH船纵向运动数值计算结果的确认等级。在此基础上，进一步开展了该SWATH船在2个波高多个波长条件下的运动求解，计算结果均与试验结果吻合较好，说明所采用的数值计算模型对于求解SWATH船在波浪中的运动问题具有较好的适用性和可靠性。

小水线面双体船；RANS；不确定度分析；波浪中运动；验证和确认

0　引　言

近年来，随着计算机技术的发展，CFD数值计算在船舶与海洋工程研究领域的应用日益广泛，与此同时，其计算结果的可信度也受到了越来越多的关注。作为CFD数值计算结果可信度的量化指标，CFD不确定度最早于1968年被提出。在船舶CFD领域，Coleman等［1］于1997年结合AIAA的CFD不确定度规程开创性地提出了适合船舶CFD不确定度分析的方法。以该方法为基础，ITTC于1999年颁布了临时规程，该规程在2000年的Gothenburg数值船舶流体力学研讨会上被推荐使用。经过广泛的实践，ITTC又于2002年发布了修订的推荐规程。此后，Simonsen等［2］应用推荐的规程对油轮ESSO Osaka数值模拟的不确定度进行了分析；Weymouth等［3］采用基于RANS方程的数值计算方法对WigleyⅢ模型在波浪中的运动进行了数值模拟，并结合ITTC推荐的规程对耐波性数值计算结果进行了不确定度分析；Guo等［4］采用ITTC推荐的规程对KVLCC2的耐波性数值计算结果开展了不确定分析。

针对CFD不确定度的研究，国内的起步相对较晚。朱德祥等［5］应用ITTC的CFD不确定度分析推荐规程对SUBOFF潜艇模型粘性绕流场的数值计算结果进行了不确定度分析；张楠等［6］采用5套网格对SUBOFF的阻力及流场的CFD不确定度进行了探讨。杨仁友等［7］参考ITTC的CFD不确定度分析推荐规程，对螺旋桨敞水性能数值计算结果进行了验证和确认；杨春蕾等［8］分别采用基于RANS和DES方程的数值计算方法，对船舶在静水中航行进行了数值模拟，并采用ITTC的CFD不确定度分析推荐规程，对船体阻力及兴波波形的数值计算结果进行了不确定度分析。

由上述研究现状不难发现，在数值研究中，对计算结果进行不确定度分析已成为一项重要且必须的工作。而我国在CFD不确定度分析研究领域与国外相比还有较大的差距，尤其是针对一些较为复杂的非定常问题，如船舶耐波性预报。目前，针对该问题数值计算结果不确定度分析的相关研究还较少，而对于一些高性能船舶，如小水线面双体船（Small Waterplane Area Twin Hull，SWATH）等，其耐波性数值计算结果的不确定度分析相关研究则更为罕见。

本文拟采用基于RANS方程的粘性CFD方法，对小水线面双体船在顶浪规则波中的纵向运动进行数值求解，进而参考ITTC推荐的CFD不确定度分析规程，对小水线面双体船在波浪中的运动响应数值计算结果进行不确定度分析，发现各输入参数对计算结果不确定度的影响规律，建立小水线面双体船耐波性数值计算结果的验证和确认等级。

1　数值计算

1.1控制方程及湍流模型

RANS方程是粘性流体运动学和动力学的控制方程，本文基于STAR-CCM+商用CFD平台，将以RANS方程作为求解船体粘性兴波流场的基本方程来进行计算。其具体形式如下：

式中：ρ为流体密度，本文中的流体为液态水；μ为流体粘度；p为静压；fi为单位质量的质量力；ui，uj为速度分量。湍流模式采用RNGk-ε模型。采用流体体积法（Volume of Fluid Method，VOF）求解兴波自由面［9］。

1.2造波与消波

采用在入口边界模拟柔性造波板运动的速度分布产生入射波。根据无限水深中的线性波浪理论，规则波的自由波面可以表达为

速度场为

式中：ζ为自由面各点垂向位置；u，v和w分别为流体各点的纵向、水平方向和垂向的速度分量；a为波幅；k为波数；ω为波浪圆频率。

同时，在数值波浪水池出口处设有消波区，通过在动量源函数中加入阻尼项，对流体质点垂向速度作强制衰减。

其中

2　CFD不确定度分析规程

在数值计算结果的验证和确认中，定义数值模拟误差δS为数值模拟结果S与真值T之差。数值模拟误差由2部分组成，分别为数值误差δSN和模型误差δSM，即

在一定的条件下，数值误差可以估计为

对数值模拟结果进行修正，得到数值基准

对于未经修正的数值模拟方法，验证的工作就是估计数值模拟的不确定度USN。此时，数值误差由来自迭代次数、网格划分、时间步长以及其他参数的误差组成，数值模拟不确定度USN可表示为

式中，UI，UG，UT和UP分别为迭代次数、网格划分、时间步长以及其他参数的误差。

对于修正过的数值模拟方法，验证的工作除了估计数值模拟不确定度USCN以外，还包括对数值误差的估计，由此可得到数值基准SC。不确定度USCN以及数值误差由下式得到：

确认则是利用基准试验数据D对数值模拟的模型不确定度USM进行评估的过程。当条件允许时，还要估计模型不确定度δSM。在确认过程中，有2个重要的参数，分别为比较误差E和确认不确定度UV，其中，比较误差E为试验数据D与模拟结果S之差。确认不确定度UV是对所有误差不确定度的综合描述，式（13）和式（14）给出了比较误差E和确认不确定度UV的表达。

对于修正的数值模拟结果，

可通过比较E和UV来判断确认是否实现。如果||E＜UV，即说明试验数据及模拟结果中所有误差的组合小于UV，此时，确认在UV水平上得到实现；如果UV≪||E，即说明数值模拟误差δN和试验数据误差δD较数值模型误差δSM更小，这就导致E≈δSM，由此，需利用比较误差E来对模型进行改进。

3　小水线面双体船耐波性数值计算结果的不确定度分析

3.1研究对象及流场设置

本文的研究对象为一艘小水线面双体船，安装有前、后稳定鳍和舭龙骨，其三维效果图如图1所示，模型水线长超过2.5 m，排水量超过200 kg。

图1　模型三维示意图Fig.1　Three-dimensional model

由于船体左右两舷对称，为减少计算量，以一半船体为计算对象。采用六自由度求解器求解船体在波浪中的垂荡和纵摇运动［11］。求解运动过程中，使用了2个坐标系：地球坐标系和随船坐标系。地球坐标系的原点位置与模型静止时右侧片体龙骨中点所在位置重合；随船坐标系与船体固联，其原点位于船体重心位置。计算流域在地球坐标系下为-（0.5λ+2.5LWL）＜x＜0.5LWL+λ；-LWL＜y＜yG；-LWL＜z＜D+0.5LWL。其中：λ为入射波长，yG为地球坐标系下船体重心横向位置坐标。计算流域及边界条件设置如图2所示。船体的运动模拟应用了Overset重叠网格技术。

图2　计算流域及边界条件设置Fig.2　Computational domain and boundary conditions

3.2计算结果的不确定度分析

本章数值计算结果的不确定度分析分验证和确认2部分，其中验证部分包括网格收敛性验证和时间步长收敛性验证。由于迭代误差远小于来自网格和时间步长的误差［3-4，6］，故计算模型的迭代误差忽略不计。

本节采用3套网格对网格的收敛性进行了研究，如图3所示。3套网格加细比rG=，网格数分别为2.7×106，1.2×106和4×105，根据网格的疏密，分别将其命名为细网格、中网格及粗网格。

图3　网格收敛性研究中所用网格Fig.3　Grids used for grid convergence study

3套网格在生成过程中均对自由面及船体附近区域进行了加密。其中对于细网格，船体上游及附近每个波长不少于40个网格，每个波高不少于20个网格。

在网格收敛性研究中，计算时间步长Δt= Te/200，其中Te为模型在波浪中运动的遭遇周期。数值模拟时间大于30倍遭遇周期。为减小样本误差对分析结果的影响，用于分析的数据均来源于数值模拟中最后12个遭遇周期，此时，入射波以及船体各响应已较为稳定。图4所示为通过细网格计算得到的，当λ/LWL=3，H=1/50LWL时船体前方0.7 m处的入射波高以及小水线面双体船各响应的时历曲线。由图可以看出，此时的波高以及船体响应均较为平稳，满足分析研究需求。

图4　入射波高及各响应时历曲线（细网格，Fr=0.236，λ/LWL=3，H=1/50LWL）Fig.4　The wave height and the time history curves of ship response（fine grid，Fr=0.236，λ/LWL=3，H=1/50LWL）

在时间步长收敛性研究中，本节采用细网格作为计算网格，选用了加细比rT=的3个时间步长设置，分别为Δt=Te/200，Te/140和Te/100。数值模拟的物理时间及分析数据的选取均与网格收敛性研究相同。

本文研究的小水线面双体船在波浪中的运动响应为垂荡、纵摇运动以及船体重心处垂向加速度的一阶运动响应幅值，分别为Za，θa和 Aa。同时，为便于比较分析，定义船体运动传递函数如表1所示。表中：ζa为入射波波幅；g为重力加速度。另外，上标“'”用于表达各响应传递函数，如Za'表示垂荡运动传递函数。

表1　船体运动传递函数Tab.1　Ship motion transfer functions

本文不确定度分析的对象为船体前方0.7 m处入射波高以及船体在波浪中的纵向运动响应传递函数计算结果。考虑到船体在运动共振区的运动幅度大，计算结果数值灵敏度较高，不确定度分析对应的计算工况为：Fr=0.236，顶浪，入射波高H=（1/50）LWL，λ/LWL=3，此时，船体处于运动共振区。

3.2.1验证

1）网格收敛性研究。

有关本文数值计算结果的相对误差，波高的数值计算结果是与理论值进行比较，运动则是和试验值进行比较，为便于叙述，本文以“理论/试验值”来表示在分析过程中所包含的波高理论值以及运动试验结果。表2所示为通过3套不同疏密网格计算得到的入射波高以及小水线面双体船运动响应计算结果与理论/试验值的对比。表中：S1，S2，S3分别代表网格收敛性研究中细网格、中网格和粗网格的数值计算结果；D代表理论/试验值；e表示数值计算结果与理论/试验值的相对误差。由表2可以看出，随着网格的加密，入射波高和船体运动响应计算结果均逐渐增大。在增大的过程中，入射波高向其理论值逼近，同时垂荡和重心处垂向加速度的计算结果向其试验值逼近，当网格为细网格时，其误差分别为-4.79%，-5.69%和-4.65%；而纵摇运动响应则在增大的过程中跨越试验值，细网格计算结果误差为5.16%。

表2　网格收敛性研究中入射波高及船体运动响应数值计算结果Tab.2　Numerical results of wave height and ship motion responses in grid convergence study

表3所示为小水线面双体船运动响应的网格收敛性不确定度分析结果。表中：下标G表示网格收敛性研究；RG为网格收敛性研究中计算结果的收敛因子，计算结果的单调收敛条件为0＜RG＜1；PG为计算结果精度阶数；CG为计算结果修正因子。当CG接近于1时，认为计算结果距离渐进区较近，此时，适合采用修正后的计算结果进行不确定度分析；而当CG远大于或小于1时，计算结果距离渐近区较远或不在渐近区，此时，采用网格误差d*G对计算结果进行修正所得到的修正计算结果的可信度相对较低，由此基于修正后计算结果的不确定度分析参考意义不大。

表3　对波高和船体运动响应计算结果的网格收敛性验证Tab.3　Grid convergence verification of waveheight and ship motion responses

由表3可知，入射波高和船体运动响应的数值计算结果均满足单调收敛条件，由此计算结果的网格收敛性得到验证，进而便可利用广义Richardson外推法［2］计算精度阶数及误差估计值。同时由表中还可看到，入射波高和船体运动响应计算结果的精度阶数PG均大于2，纵摇运动计算结果的精度阶数最高，达到了5.24，这说明纵摇运动的计算结果随着网格的加密收敛速度也加快了。

另外还可看出，纵摇运动计算结果的修正因子CG达到了5.15，说明其计算结果距离渐进区较远或不在渐近区。由表3给出的计算结果的网格不确定度，发现数值计算结果的网格不确定度UG均较小，垂荡运动响应计算结果的网格不确定度仅为计算值S1的1.35%，相比之下，纵摇运动计算结果的网格不确定度较大，达到了计算值S1的3.82%。另外，各运动响应修正结果的网格不确定度较未修正时均明显减小，其中修正后的垂荡运动计算结果的网格不确定度仅为计算值S1的0.23%。

2）时间步长收敛性研究。

表4所示为利用细网格分别选取时间步长Δt=Te/200，Te/140，Te/100的数值模拟所得到的入射波高及船体运动响应计算结果与理论/试验值的对比。表中，S1，S2和S3分别代表时间步长收敛性研究中时间步长Δt=Te/200，Te/140和Te/100时的数值计算结果。如表4所示，随着时间步长的减小，入射波高和船体运动响应计算结果均逐渐增大，在增大的过程中，入射波高以及垂荡和重心处垂向加速度的计算结果分别向其理论值和试验结果逼近，而纵摇运动响应则在增大的过程中跨越了试验值。

表4　时间步长收敛性研究中入射波高及船体响应数值计算结果Tab.4　Numerical results of wave height and ship motion responses in time step convergence study

表5所示为小水线面双体船运动响应时间离散的不确定度分析结果。表中，下标T表示时间步长收敛性研究。从表中可以看出，垂荡、纵摇和重心处垂向加速度数值的计算结果均满足单调收敛条件，即0＜RT＜1。与在网格收敛性研究中发现的规律类似，纵摇运动计算结果的精度阶数PT较其他量更大，达到了4.45，同时其修正因子CT达到了3.68，明显大于1，说明在时间步长收敛性研究中，纵摇运动计算结果距离渐近区较远。

表5　对波高和船体运动响应计算结果的时间步长收敛性验证Tab.5　Time step convergence verification of waveheight and ship motion responses

值得注意的是，相比于网格不确定度，数值计算结果的时间步长不确定度更大，其中垂荡运动尤为显著，其未修正情况下的时间步长不确定度达到网格不确定度的5.34倍，修正情况下为3.35倍。该现象说明相比于网格的疏密，小水线面双体船在波浪中运动的数值求解模型对时间步长的大小更为敏感，对于该类问题，建议在计算资源有限的情况下首先保证时间步长的精度要求。

3.2.2确认

表6所示为入射波高和小水线面双体船在波浪中运动响应数值计算结果的确认。表中，确认不确定度UV2=UD

2+UT2+UG

2，其中UD为试验不确定度，本文中所用试验结果的试验不确定度为2.5%，由于入射波高计算结果是与理论值进行对比，故本文波高的试验不确定度为0。

表6　波高及船体运动响应数值计算结果的确认Tab.6　Validation of waveheight and ship motion responses

由表6可知，入射波高、垂荡以及重心处垂向加速度的数值计算结果在修正和未修正的情况下均得到了确认，其未修正计算结果的确认水平分别为理论/试验值D的6.33%，7.35%和10.25%，修正计算结果的确认水平分别为D的1.04%，2.61% 和2.96%。不难发现，修正后结果的确认不确定度较修正前均显著降低，其中入射波高修正后结果的确认不确定度仅为修正前的16.4%。入射波高以及除纵摇运动外其他运动响应修正后计算结果的比较误差较未修正时明显降低，其中修正后垂荡运动计算结果的比较误差EC为D的0.03%，仅为未修正时的0.53%。

对于纵摇运动响应，其未修正的计算结果在D的6.52%水平得到确认，而修正结果则未得到确认。可以看出，修正后纵摇运动计算结果的确认不确定度UVC较修正前明显减小，但由于纵摇运动计算结果在随网格和时间步长的收敛过程中跨越了试验结果，使得修正后得到的数值基准SC相比S1较试验值偏差更大，由此造成修正后的比较误差EC较未修正时增大，从而导致EC＞UVC。尽管修正后的纵摇运动响应计算结果未得到确认，但并不意味着纵摇运动计算结果不可靠。其原因为纵摇运动的修正因子与1偏离较大，在网格收敛性研究中，纵摇运动计算结果的修正因子CG达到了5.15，在时间步长收敛性研究中，CT=3.68，说明纵摇运动计算结果距离渐进区较远。对于该种情况，纵摇运动在未修正时的不确定度分析结果的参考意义更大。

综上所述，对于小水线面双体船波浪中运动问题的求解，采用基于RANS方程的数值计算方法得到的数值计算结果的网格收敛性和时间步长收敛性均得到了验证。在未对结果进行修正的情况下，计算结果均得到确认，且确认不确定度最大不超过11%，其中3项小于8%，比较误差均小于6%；对于修正后的计算结果，其确认不确定度最大不超过4%，比较误差除纵摇运动外均小于3%。尽管纵摇运动计算结果未得到确认，其比较误差仍小于10%，且由于纵摇运动计算结果偏离渐近区较远，故其在未修正情况下的不确定度分析结果的参考意义更大。由此，本文认为采用基于RANS方程的数值计算方法求解小水线面双体船波浪中运动问题的结果可靠。

为进一步确认该数值计算模型对于处理小水线面双体船在波浪中运动问题的适用性和可靠性，本文选用细网格作为计算网格，设计算时间步长Δt=Te/200，开展了在入射波高H=LWL/50和H= LWL/30时，该小水线面双体船在不同入射波长条件下船体运动响应的数值求解，并将数值计算结果与试验结果进行了对比，如图5所示。由图可见，计算所得运动响应的整体趋势与试验结果一致，吻合较好，从而进一步说明本文采用的数值计算模型对于处理小水线面双体船波浪中的运动问题适用且可靠。

图5　运动响应计算结果与试验结果对比（Fr=0.236）Fig.5　Comparison of motion transfer functions between calculated results and experimental data（Fr=0.236）

4　结　论

本文采用基于RANS方程的数值计算方法对一艘小水线面双体船在波浪中的纵向运动开展了数值求解，并参考ITTC推荐的CFD不确定度分析规程对计算结果进行了不确定度分析，得到以下结论：

1）验证了计算结果的网格收敛性和时间步长收敛性，发现相比于网格的疏密，小水线面双体船波浪中运动求解模型对于时间步长更为敏感，建议在计算资源有限的情况下首先满足时间步长的精度要求。

2）在未对结果进行修正的情况下，计算结果均得到了确认，且确认不确定度最大不超过11%，其中3项小于8%，比较误差均小于6%；对于修正后的计算结果，其确认不确定度最大不超过4%，比较误差除纵摇运动外均小于3%，纵摇运动的比较误差小于10%。

3）该小水线面双体船在2个波高，多个波长条件下的运动求解结果均与试验结果吻合较好，进一步证明了本文所采用数值计算模型对于求解小水线面双体船波浪中运动问题的适用性和可靠性。

［1］COLEMAN H W，STERN F.Uncertainties in CFD code validation［J］.Journal of Fluids Engineering，1997，119（4）：795-803.

［2］SIMONSEN C D，STERN F.Verification and valida⁃tion of RANS maneuvering simulation of Esso Osaka：effects of drift and rudder angle on forces and moments ［J］.Computers&Fluids，2003，32（10）：1325-1356.

［3］ WEYMOUTH G D，WILSON R V，STERN F.RANS computational fluid dynamics predictions of pitch and heave ship motions in head seas［J］.Journal of Ship Re⁃search，2005，49（2）：80-97.

［4］GUO B J，STEEN S，DENG G B.Seakeeping predic⁃tion of KVLCC2 in head waves with RANS［J］.Applied Ocean Research，2012，35：56-67.

［5］ 朱德祥，张志荣，吴乘胜，等.船舶CFD不确定度分析及ITTC临时规程的初步应用［J］.水动力学研究与进展，2007，22（3）：363-370.

ZHU Dexiang，ZHANG Zhirong，WU Chengsheng，et al.Uncertainty analysis in ship CFD and the primary application of ITTC procedures［J］.Journal of Hydrody⁃namics，2007，22（3）：363-370.

［6］张楠，沈泓萃，姚惠之.阻力和流场的CFD不确定度分析探讨［J］.船舶力学，2008，12（2）：211-224.

ZHANG Nan，SHEN Hongcui，YAO Huizhi.Uncer⁃tainty analysis in CFD for resistance and flow field［J］. Journal of Ship Mechanics，2008，12（2）：211-224.

［7］杨仁友，沈泓萃，姚惠之.螺旋桨敞水性能CFD不确定度分析［J］.船舶力学，2010，14（5）：472-480.

YANG Renyou，SHEN Hongcui，YAO Huizhi.Uncer⁃tain analysis of CFD simulation on the open-water per⁃formance of the propeller［J］.Journal of Ship Mechan⁃ics，2010，14（5）：472-480.

［8］ 杨春蕾，朱仁传，缪国平，等.基于RANS和DES法船体绕流模拟及不确定度分析［J］.上海交通大学学报，2012，46（3）：430-435.

YANG Chunlei，ZHU Renchuan，MIAO Guoping，et al.Uncertainty analysis in CFD for flow simulation around ship using RANS and DES［J］.Journal of Shanghai Jiaotong University，2012，46（3）：430-435.

［9］RUSCHE H.Computational fluid dynamics of dispersed two-phase flows at high phase fractions［D］.London：University of London，2002.

［10］ CHOI J，YOON S B.Numerical simulations using mo⁃mentum source wave-maker applied to RANS equa⁃tion model［J］.Coastal Engineering，2009，56（10）：1043-1060.

［11］ CARRICA P M，WILSON R V，NOACK P W，et al. Ship motions using single-phase level set with dynam⁃ic overset grids［J］.Computers&Fluids，2007，36 （9）：1415-1433.

Uncertainty analysis in CFD for SWATH motions in regular head waves

DENG Lei1，PENG Hongyu2

1 The 92692thUnit of PLA，Zhanjiang 524064，China
2 Naval Military Representative Office in Jiangnan Shipyourd（Group）Co.，Ltd.，Shanghai 201913，China

Based on procedures for uncertainty analysis in CFD recommended by ITTC，the results of Small Waterplane Area Twin Hull（SWATH）ship longitudinal motions in regular head waves solved by RANS code are verified and validated.The verification both for grid convergence and time step conver⁃gence is investigated respectively.Finding shows that the numerical solving model is more sensitive to time step size than that of grid.According to the experimental data，numerical results are validated mostly.And then the levels of verification and validation for numerical results of SWATH ship longitudinal motions in regular head waves are established.Further，SWATH ship motions in conditions of two different wave heights as well as various wavelengths are carried out.The good agreement with experimental data showed by numerical results further proved that the proposed approach is effective and reliable to solve issues of SWATH ship motion in waves.

Small Waterplane Area Twin Hull（SWATH）；RANS；uncertainty analysis；motion in waves；verification and validation

U661.3

A

10.3969/j.issn.1673-3185.2016.03.004

2015-08-22网络出版时间：2016-5-31 11:04

邓磊（通信作者），男，1990年生，硕士。研究方向：船舶流体动力性能。

E-mail：hgdenglei@163.com