赵发明，高成君，夏 琼

（中国船舶科学研究中心，江苏 无锡 214082）

1 前 言

计算粘性自由面绕流问题是船舶水动力学领域最重要的课题之一，一直是一项挑战性的任务。目前随着CFD数值技术的不断提高和计算机硬件的飞速发展，近年来“数值水池”发展很快，但网格生成技术和自由液面的模拟方法仍是关键技术之一。

由于结构网格技术成熟、壁面粘性模拟能力强等优点，因此即使在非结构、笛卡儿直角网格逐渐盛行的今天，仍然是船舶CFD计算的首选。为了更好地利用结构网格的优点，人们迫切希望开发结构网格的自动化生成方法，于是出现了多块结构化网格生成方法，但该方法仍然没有真正解决结构网格的自动化生成技术。现在尽管商业软件（如Fluent/Gambit）在非结构网格、多块结构网格及动网格方面具有强大功能，但由于船舶CFD主要解决两相流问题，要求在自由液面附近的网格具有很高的分辨率，这一特殊性使得商用软件在船舶六自由度运动、具有较大航态变化的高性能船舶（如快艇、高速复合船）方面模拟具有局限性，若再考虑带附体、边界限制的影响，CFD模拟更是困难。重叠网格（chimera grid，overset grid，overlap grid）的出现使结构网格的自动化成为可能，上面的难题可获得较容易的解决：如远场背景网格采用笛卡儿直角网格以保证自由液面的网格分辨率，各物体采用贴体网格，这样就可解决船舶CFD模拟问题。

在船舶CFD模拟中，自由面处理方法可分为自由面追踪法（free surface tracing）和自由面捕捉法（free surface capturing）。自由面追踪法采用移动网格技术，该方法使网格变形以适合自由面边界的形状，但对自由面大幅变形情况的数值模拟比较困难。自由液面捕捉方法用混合的空间网格并且通过解决边界面处的运输方程来捕获自由面位置，最常用的界面捕获方法是MAC方法、VOF方法、Level set方法和VOF/Level set混合方法。本文采用单相Level Set方法来模拟自由液面形状。

2 重叠网格方法

重叠网格方法将复杂的流动区域分成几何边界比较简单的子区域，各子区域中的计算网格独立生成，彼此存在着重叠、嵌套或覆盖关系，流场信息通过插值在重叠区边界进行匹配和耦合。重叠网格即拥有结构网格逻辑关系简单、流场计算精度高、效率高、壁面粘性模拟能力强等优点，更弥补了结构网格对外形适应能力差的缺点。重叠网格的应用越来越广泛，国外比较著名的重叠代码：如NASA的PEGASUS、ARL的SUGGAR及SUGGAR++、OVERTURE等，并且还在不断地推陈出新。

重叠网格方法的最终目标是建立各重叠网格间的耦合关系，为各区域流场计算提供边界信息的传递服务。如图1所示：网格1落入网格2物面内的点被标记为洞内点而不参与流场的计算。网格1中与洞点相邻的点为洞边界点（图1中的方点），这些点通过插值接受从网格2传递的流场信息。相应地，网格2的外边界点（图1中的圆点）也通过插值接受从网格1传递的流场信息，任意变量f由相应的六面体贡献单元采用trilinear方法插值获得（1式）。网格1的洞边界点和网格2的外边界插值点之间的区域为重叠区域。图2为典型单体船静态的船舶重叠网格，图3为带附体的船舶重叠网格。

图1 重叠网格示意图Fig.1 Overlap grid

其中：di，dj，dk 为[0,1]间的插值权系数。

图2 单体船重叠网格（x＝0.5处剖面）Fig.2 Overlap grid for monohull(x＝0.5 section)

图3 带附体的船舶重叠网格Fig.3 Overlap grid for ship with appendage

重叠网格生成的工作并不是简单地等于各子网格的生成，更为重要的内容是确定各子网格之间的重叠关系，以及确定各子网格之间的数据插值传递关系，因此重叠网格包含两个主要的步骤：挖洞和插值。挖洞的目的是在流场计算前从网格中屏蔽掉一些不必要或者无实际意义的部分（如物体壁面内部）。具体而言，是在需要屏蔽的区域周围设定挖洞面，然后将落入挖洞面内的网格点标识出来，并在CFD计算过程中予以舍弃。挖洞方法可以分为指定挖洞面法、洞映射方法以及叉树方法，本文采用洞映射方法。寻点是在网格中寻找插值点的贡献单元。组成贡献单元的顶点往往是参与CFD流场计算的洞外网格点，且在特殊情况下也可以是洞边界点。寻找贡献单元最简单、最可靠的方法是遍历整个网格域，直到找到正确的网格单元。但是，这种方法的效率是最低的，对于数据量稍大的网格，其搜索速度异常缓慢。利用优秀的数据结构可成量级地提高寻点速度，现在常用的方法为叉树法和ADT方法，本文采用ADT方法。

3 单相Level Set方法

3.1 自由液面捕捉

定义Level Set函数φ（流场任意一点到自由液面距离的函数），满足：

通过插值获得φ=0（自由液面）。忽略空气的影响，则在φ=0处应满足如下条件：

船体流动的求解仅在φ≤0的计算域进行。

3.2 Level Set函数重初始化

Level Set函数重初始化是单相Level Set方法中的重要步骤。（2）式求解后可获得φ=0即新的自由液面，若不周期性重初始化函数φ，则φ始终保持为初设函数，即流场任意一点到z=0自由液面的函数，这样仅对于小幅波的计算可行。Level Set函数φ的周期性重初始化就是重新获得任意一点到φ0=0的距离函数（图4）：

式中φ0为重置前的Level Set函数。

4 RANS控制方程与湍流模型

图4 自由液面重初始化图Fig.4 Reinitialization of free surface

4.1 控制方程

Reynolds模仿气体分子运动论里的平均概念，对不可压缩流体的NS方程进行时间平均，得到著名的 Reynolds Averaged Navier Stokes（RANS）方程。

不可压缩流体无量纲RANS方程的张量形式为：

4.2 湍流模型

湍流模型为k-ω两方程模型（标准型和SST型）。

湍流动能k方程和特殊耗散率ω方程为：

式中，初值 k=1.1×10-7，ω=9.0。

5 离散方法

控制方程采用体积中心有限差分格式进行离散，其中所有时间项采用二阶欧拉向后差分格式：

（7）式中的对流项采用高阶差分格式（二阶或三阶迎风格式)或混合差分格式离散，本文中计算采用一、二阶及三阶迎风混合差分格式，（7）式中粘流项采用二阶中心差分格式离散。（2）式采用二阶迎风差分格式离散。方程的离散求解采用ADI方法和PISO算法。

6 船舶六自由度运动求解

求解船舶的运动情况是采用求解运动方程来得到。为此，我们首先要了解船舶在三个旋转轴的转动惯量。如果我们将载体坐标系的原点定在船舶的旋转点位置上，那么三个方向的转动惯量就可以从（11）式得到。

式中，rg,j为j轴方向的惯量半径，xG=xrot-xCG旋转点与重心位置纵向距离。

根据船体受力平衡，其刚体运动方程见（12）-（13）式。船体运动后(船舶贴体网格部分)，网格不需要重新生成，仅仅重新生成网格间的插值信息即可，船舶六自由度运动求解流程见图5，图6为典型单体船运动的网格示意图。

式中，X，Y，Z，K，M，N 分别为纵荡、横荡、垂荡力和横摇、纵摇及艏摇力矩。 u，v，w 为 x，y，z三个坐标方向的线速度，为线加速度。p，q，r为欧拉角为欧拉角加速度。

图5 六自由度运动求解流程Fig.5 6 Dof solution strategy

图6 动网格重叠示意图Fig.6 Overset grid for motion

7 计算验证

7.1 单体船算例

采用重叠网格和单相Level Set方法对60约束船模和自由模在Fr=0.316，Re=4.3×106进行数值模拟。计算网格状态均为：背景网格采用笛卡尔网格101×36×62，物体贴体网格采用双O型网格74×36×59，网格总数346 248，重叠网格见图7所示。计算未采用壁面函数，y+≈0.8。计算无量纲时间步长Δt=0.01，总时间T=20（其中加速段总时间T1=2，船前进19倍的船长），计算平台为PC Core Quad CPU 2.33G+4.0GRAM，平均计算耗时5～10小时。表1为计算阻力与试验值的比较。其中无量纲升沉和纵倾的定义见（14）式。

表1 阻力比较Tab.1 Drag compare table

图7 S60舷侧波型和升沉纵倾曲线图（Fr=0.316，Re=4.3×106）Fig.7 Wave contours&sinkage and trim for S60(Fr=0.316,Re=4.3×106)

7.2 双体船算例

采用重叠网格和单相Level Set方法对双体船约束船模在Fr=0.55，Re=1.722×106进行数值模拟，片体中心距s=0.2L。计算网格状态均为：背景网格采用笛卡尔网格101×56×51，物体贴体网格采用双O型网格2×71×48×50，网格总数644 709，重叠网格见图8所示。计算未采用壁面函数，y+≈0.5。计算无量纲时间步长Δt=0.01，总时间T=10（其中加速段总时间T1=2，船前进9倍的船长），计算平台为PC Core Quad CPU 2.33G+4.0GRAM，平均计算耗时8小时。表2为计算阻力与文献[8]中提供的试验值的比较。图9为舷侧波型和波高等值线图。

图8 双体船重叠网格Fig.8 Overlap grid for catamaran

表2 阻力比较Tab.2 Drag compare table

图9 双体船舷侧波型和波高等值线图（Fr=0.55，Re=1.722×106）Fig.9 Wave contours and wave profile for catamaran(Fr=0.55,Re=1.722×106)

7.3 三体船算例

采用重叠网格和单相Level Set方法对三体船约束船模在Fr=0.3，Re=4.6×106进行数值模拟。计算状态为：背景网格采用笛卡尔网格80×45×46，主体和片体的贴体网格采用双O型网格3×71×51×39，总网格数589 257，图10为三体船重叠网格示意图。计算未采用壁面函数，，y+≈0.8。计算无量纲时间步长Δt=0.01，总时间T=10（其加速段总时间T1=2，船前进9倍的船长），计算平台为PC Core Quad CPU 2.33G+4.0GRAM，计算耗时8小时。图11-13分别为波型云图、阻力时间历程图和舷侧波型。由图10可以看出自由面在重叠区域光滑过渡，表明了重叠网格区间数据的传递是正确的，图13中主船体及片体内外舷侧波型从理论上分析也是合理的。

图10 三体船重叠网格Fig.10 Overlap grid for trimaran

图11 三体船波高等值线及云图 Fig.11 Wave contours for trimaran

图12 总阻力时间历程图Fig.12 Computed time for trimaran

8 结论与展望

本文初步尝试采用重叠网格和单相Level-Set方法，结合PISO算法、k-ω湍流模型，利用RANS方程成功地数值模拟了船体带自由面流动问题，单体和双体船重叠算例计算得出的波形图和阻力系数与试验值相差很小，三体船计算结果尽管未与试验结果对比，但理论上分析也比较合理。通过自由态的单体船初步数值模拟，表明了采用Level set方法处理自由液面和重叠网格技巧来处理动网格问题，其数值模拟结果非常好，证明了该方法处理动态问题是有效的。重叠网格和单相Level-Set是一种实用的船舶CFD方法，我们后续目标是将该方法应用实际工程中复杂船型的六自由度数值模拟。

图13 三体船主/片体舷侧波型Fig.13 Wave profile for trimaran

[1]王献孚等.计算船舶流体力学[M].上海:上海交通大学出版社,1992.

[2]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-17.

[3]Rogers S E.PEGASUS user’s guide[K].http://people.nas.nasa.gov/～rogers/pegasus/uguide.html,Oct,2003.

[4]Rogers S E,Suhs N E,Dietz W E.PEGASUS 5:An automated preprocessor for overset-grid computational fluid dynamics[J].AIAA Journal,2003,41(6):25-51.

[5]Paterson E G,Wilson R V,Stern F.General-purpose parallel unsteady rans ship hydrodynamics code:CFDSHIP-IOWA[R].IIHR Technical Report No.432,2003.

[6]Carrica P M,Wilson R V,Noack R,Xing T,Kandasamy M,Shao J,Sakamoto N,Stern F.A dynamic overset singlephase level set approach for viscous ship flows and large amplitude motions and maneuvering[C]//Th26 Symposium on Naval Hydrodynamics 2006.Korea,2006.

[7]Chi Yang,et al.Calculation of ship sinkage and trim using unstructured grids[C]//ECCOMAS 2000,11-14 September 2000.Barcelona,2000.

[8]高井智祐（田原研究室）,07後期中間发表会.ppt[J/OL].2007.