基于重叠网格的船舶粘流数值模拟_图文

1、1前言计算粘性自由面绕流问题是船舶水动力学领域最重要的课题之一,一直是一项挑战性的任务。 目前随着 CFD 数值技术的不断提高和计算机硬件的飞速发展,近年来“数值水池”发展很快,但网 格生成技术和自由液面的模拟方法仍是关键技术之一。结构网格技术成熟、壁面粘性模拟能力强等优点,因此在非结构、笛卡儿直角网格逐渐盛行的 今天,仍然是船舶 CFD 计算的首选。为了更好地利用结构网格的优点,人们迫切希望开发结构网格 的自动化生成方法,于是出现了多块结构化网格生成方法,但该方法仍然没有真正解决结构网格的 自动化生成技术。现在尽管商业软件(如 Fluent/Gambit在非结构网格、多块结构网格及动网格方

2、面具有强大功能,但由于船舶 CFD 主要解决两相流问题,要求在自由液面附近的网格具有很高的分 辨率,这一特殊性使得商用软件在船舶六自由度运动、具有较大航态变化的高性能船舶(如快艇、 高速复合船型、多体相对运动以及限制条件情况等方面模拟具有局限性,若再考虑带附体、边界 限制的影响,CFD 模拟更是困难。重叠网格(chimera grid 、overset grid、overlap grid的出现 使结构网格的自动化成为可能,上面的难题可获得较容易的解决:如远场背景网格采用笛卡儿直角 网格以保证自由液面的网格分辨率,各物体采用贴体网格,这样就可解决船舶 CFD 模拟问题。 在船舶 CFD 中模拟中

3、,自由面处理方法可分为自由面追踪法(free surface tracing和自由面 捕捉法(free surface capturing。自由面追踪法采用移动网格技术,该方法使网格变形以适合自 由面边界的形状,但对自由面大幅变形情况的数值模拟比较困难。自由液面捕捉方法用混合的空间 网格并且通过解决边界面处的运输方程来捕获自由面位置,最常用的界面捕获方法是 MAC 方法、VOF 方法、Level set方法和 VOF/Level set混合方法。本文采用单相 Level Set方法来模拟自由液面形 状。2重叠网格方法2.1重叠网格方法介绍重叠网格方法将复杂的流动区域分成几何边界比较简单的子区域

4、,各子区域中的计算网格独立 生成,彼此存在着重叠、嵌套或覆盖关系,流场信息通过插值在重叠区边界进行匹配和耦合。重叠 网格即拥有结构网格逻辑关系简单、流场计算精度高、效率高、壁面粘性模拟能力强等优点,更弥 补了结构网格对外形适应能力差的缺点。重叠网格的应用越来越广泛,国外比较著名的重叠代码:如 NASA 的 PEGASUS、ARL 的 SUGGAR 及 SUGGAR+、OVERTURE 等,并且还在不断地推陈出新。12重叠网格方法的最终目标是建立各重叠网格间的耦合关系,为各区域流场计算提供边界信息的 传递服务。如图 1所示:网格 1落入网格 2物面内的点被标记为洞内点而不参与流场的计算。网格 1

5、中与洞点相邻的点为洞边界点(图 1中的方点,这些点通过插值接受从网格 2传递的流场信息。 相应地,网格 2的外边界点(图 1中的圆点也通过插值接受从网格 1传递的流场信息,任意变量 f 由相应的六面体贡献单元采用 trilinear 方法插值获得(1式。网格 1的洞边界点和网格 2的外 边界插值点之间的区域为重叠区域。图 2为典型单体船静态的船舶重叠网格。=+1121415f f (-ff di+ (-ff dj+ (-ff dk+ (f-f f f djdk+ (-ff f +ff f 78f f didjdk(1其中:di,dj,dk 为0,1间的插值权系数。图 1重叠网格示意图网格 1的

6、 洞边界点 网格 2的固壁边界网格 2的外边插值点网格 2网格 1 3图 2 单体船重叠网格(x =0.5处剖面重叠网格生成的工作并不是简单地等于各子网格的生成,更为重要的内容是确定各子网格之间 的重叠关系,以及确定各子网格之间的数据插值传递关系,因此重叠网格包含两个主要的步骤:挖 洞和插值。挖洞的目的是在流场计算前从网格中屏蔽掉一些不必要或者无实际意义的部分(如物体 壁面内部。具体而言,是在需要屏蔽的区域周围设定挖洞面,然后将落入挖洞面内的网格点标识 出来, 并在 CFD 计算过程中予以舍弃。 挖洞方法可以分为指定挖洞面法、 洞映射方法以及叉树方法, 本文采用洞映射方法。寻点是在网格中寻找插

7、值点的贡献单元。寻找贡献单元最简单、最可靠的方 法是遍历整个网格域,直到找到正确的网格单元。但是,这种方法的效率是最低的,对于数据量稍 大的网格,其搜索速度异常缓慢。利用优秀的数据结构可成量级地提高寻点速度,现在常用的方法 为叉树法和 ADT 方法,本文采用 ADT 方法。3 单相 Level Set方法3.1 自由液面捕捉定义 Level Set函数 (流场任意一点到自由液面距离的函数,满足: 0i u t += (3式中: 0 00<=>空气 自由液面 液体通过插值获得 0=(自由液面。忽略空气和自由面表面张力的影响,则在自由面(0=处 应满足如下条件:0 i i j j i

8、i u x n n x x =, =- (4 2z p Fr= (5 单相 Level Set方法仅仅对船体流动的求解仅在 0的计算域进行,空气相则是通过速度扩 展(velocity extension的方法来流场速度。因为只考虑单相流场,并且只需要在界面边界处稍加4处理,这样就成功避免了上述的两相流界面的过渡问题。此外,在气体中,只需要布置少许网格来 满足 jump condition,因此相比 two phase 方法,计算资源的消耗则是大大减小,计算稳定。单相 level-set 目前被广泛应用于空气和水的自由面问题当中。4 RANS 控制方程与湍流模型4.1 控制方程Reynolds

9、模仿气体分子运动论里的平均概念,对不可压缩流体的 NS 方程进行时间平均,得到 著名的 Reynolds Averaged Navier Stokes(RANS方程。不可压缩流体无量纲 RANS 方程的张量形式为:0i i u x = (7 _2_' '1( 0Re i i i j i j j i j j j u u u p u u u t x x x x x += (85 船舶六自由度运动求解求解船舶的运动情况是采用求解运动方程来得到。为此,我们首先要了解船舶在三个旋转轴的 转动惯量。如果我们将载体坐标系的原点定在船舶的旋转点位置上,那么三个方向的转动惯量就可 以从(12式得

10、到。G G 222, G 222, G ( y+z 0 0 0 00 0 0 ( y+z 00 0 0 g x x y g y z m r I I m r I +=+G 222, G 0 ( y+z g z m r +(12 式中, , g j r 为 j 轴方向的惯量半径, G rot CG x x x =为旋转点与重心位置纵向距离。根据船体受力平衡,其刚体运动方程见(13-(14式。船体运动后(船舶贴体网格部分, 网格不需要重新生成,仅仅重新生成网格间的插值信息即可,船舶六自由度运动求解流程见图 7, 图 8为典型单体船运动的网格示意图。 22 CGi CGi CGi CGi i d x

11、d d F m L I dt dt dt = (135 图 3 六自由度运动运动求解流程图 4 动网格重叠示意图6 数值造波规则波以及长峰或短峰不规则波可通过定义初始边界条件来实现,其在任意时刻的波幅可由 (15式来表示,其为一系列不同波幅、波长和波向的简谐波线性叠加而成。(, , cos cos sin ijijjiijijx y t a k x y t =+(15式中:ij 为随机相位, ij a为波幅, i =i k 为波数, j 为入射角。入射角包括 色散角 j 和相对船的遭遇角 0,这样 0j j =+。单元波波幅 ij a 采用(16式求解。ij a = (16 式中:(i S 为

12、谱密度,可采用 ITTC 双参数谱、Bretschneider谱、JONSWAP 谱公式等。 (22 cos 22j j j D-/=, , 为方向扩展函数(ITTC推荐公式。在上述波浪初始条件下,则速度和压力的表达式见式(17:(0, 0, 0, (, , , cos cos sin cos (, , , cos cos sin sin (, , , sin cos sin (, , , i i i k z i j j i ij ji jk z i j j i ij j i j k z i j j i ij i ji U x y z t U k x y t V x y z t V k x y

13、 t W x y z t W k x y t a p x y z t =+ +=+=+=(2, cos cos sin 2i i jij k zk z i j j i ij i i ja e k x y t eFr k + (17 式中:7 静态计算验证7.1 单体船算例采用重叠网格和单相 Level Set 方法对 DTMB5415和 KCS 约束船模在 70.28, 1.1910Fr Re =×和 70.26, 1.4010Fr Re =×分别进行了数值模拟,图 9和表 1分别为计算自由面形状和阻力与试 验值的比较,可以看出结算结果与试验结果非常吻合。 计算无量纲时间步

14、长 0.01t =, 总时间 10T =(其 中加速段总时间 12T =,船前进 9倍的船长,计算平台为 PC Core Quad CPU 2.33G+4.0GRAM,计 算耗时 15小时。KCS 计算网格状态为:背景网格采用笛卡尔网格 1384160××,物体贴体网格采用 O 型网格1053961××,总数 589275,计算未采用壁面函数, 1.0y +。计算无量纲时间步长 0.01t =, 总时间 10T =(其中加速段总时间 12T =,船前进 9倍的船长,计算耗时 10小时。 图 5 DTMB5415 and KCS 自由面形状表 1 阻力比较

15、P C F C T C Cal 1.37 2.954.32 DTMB5415 EXP2.91(ITTC57 4.23 Cal 0.85 2.763.61 KCS EXP2.83(ITTC57 3.56 Cal本文 5.478 4.305 9.783 Cal(文献85.331 4.525 9.856 Wigley 双体船(s=0.2LEXP(文献85.4207.2 多体船算例本文还对双体船和三体船约束船模在 60.55, 1.72210Fr Re =×和 60.3, Re 4.610Fr =×进行数值模拟,重叠网格见图 10,计算自由面形状见图 11,双体船的阻力计算结果见表

16、 1。双体船计算网格状态为:背景网格采用笛卡尔网格 1015651××,物体贴体网格采用 O 型网格 2×71×48×50, 网格总数 644709。 计算未采用壁面函数, 1.0y +。 计算无量纲时间步长 0.01t =, 总时间 10T =(其中加速段总时间 12T =,船前进 9倍的船长,计算耗时 8小时。双体船计算网格状态为:背景网格采用笛卡尔网格 804546××, 主体和片体的贴体网格采用 O 型网格 3×71×51×39,总网格数 589257。计算未采用壁面函数, 1.0y +

17、。计算无量纲时间步长0.01t =,总时间 10T =(其加速段总时间 12T =,船前进 9倍的船长,计算耗时 7小时。 图 6 双体船和三体船重叠网格 图 7 双体和三体船自由面形状7.3 浅水及池壁影响算例本文还对浅水单体肥大型 (长宽比 /3.817L B =, 宽吃水比 /7.40B T =, 方形系数 0.83Cb = 约束船模(图 12在 60.119, 4.7510Fr Re =×进行了有无池壁效应的浅水影响下的数值模拟,水 深与船长比 /0.0354H L =, 池壁宽与船长比 /1.008b L =。 CFD 区域及计算网格见图 13, ,计算自 由面波形对比结果

18、见图14,可以看出池壁效应在浅水条件下的影响非常大。 图 8 肥大型船模 图 9 计算区域和网格无池壁影响池壁影响图 10 自由面形状对比说明:由于本船的吃水/船长比过小(见图 15且本船方形系数 Cb 过大, 故在浅水计算中很难采 用动态重叠网格进行航态计算,实际航态计算采用静态修正方法进行的,以节省计算资源。若采用 动态重叠网格计算需要很小的时间步长和较多的网格。7.4带附体近水面算例为验证物面重叠方法的有效性, 本文还对带附体的 Suboff 模型在 6 0.316, 4.010 Fr Re=×进 行了近水面情况下的数值模拟,其计算网格见图 15,计算自由面波形结果见图 16。

19、 图 11 suboff计算网格 图 12 近水面自由面波形8 动态计算验证8.1 静水自由船模阻力算例采用重叠网格和单相 Level Set 方法对 60约束船模和自由模在 60.316, 4.310Fr Re =×进行 了数值模拟。计算网格状态均为:背景网格采用笛卡尔网格 1013662××,物体贴体网格采用 O 型 网格 74×36×59,网格总数 346248。计算未采用壁面函数, 1.0y +。计算无量纲时间步长0.01t =,总时间 20T =(其中加速段总时间 12T =,船前进 19倍的船长,平均计算耗时 510小时。表 2为计

20、算阻力与试验值的比较,图 17为船体舷侧波形比较和升沉纵倾曲线图。 图 13 S60舷侧波型比较和升沉纵倾曲线图(Fr =0.316,Re =4.3×106表 2 阻力比较P C F C T C Cal 2.01 3.51 5.52 S60(fixEXP 3.42(ITTC57 5.42 Cal2.19 3.56 5.75 S60(freeEXP3.42(ITTC575.868.2 规则波中运动响应算例 本文计算是基于 Journee 在 1992 年所做的四种 WigleyIII 型船模在 Fr = 0.30, Re = 4.88 × 10 6 的迎浪试验。根据模型试验，

21、计算取 7 个不同波长 L = 0.50, 0.75,1.0,1.25,1.5,1.75, 2.00 的规则 波，只考虑垂荡和纵摇运动，为了方便和试验值比较，我们定义（18）式的无因次变量，其运动幅 值计算结果与实验结果的对比见图 19，其计算结果与试验值非常接近，这表明船舶在迎浪中运动的 数字模拟还是非常准确的。图 20 为 L = 1.25 情况下一个波浪周期中四个不同时刻下的自由面形 状变化。 Z= z a = L 2 a （18） 计算网格状态为：背景网格采用笛卡尔网格 111× 36 × 65 ，物体贴体网格采用 O 型网格 71× 36 ×

22、50 ，网格总数 387540，重叠网格见图 18 所示。计算未采用壁面函数， y + 1.0 。对不同 波长状态下，其无量纲周期 T 按(19式计算，计算无量纲时间步长取 T 80 。 1 1 1 = + T Fr 2 （19） 图 14 规则波中动态重叠网格 Z' 1.4 1.2 1 0.8 0.6 0.4 0.2 0 0.25 /L EFD CFD 0.5 0.75 1 1.25 1.5 1.75 2 2.25 2 1.8 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0 0.25 ' EFD CFD /L 0.5 0.75 1 1.25 1.5 1.75

23、2 2.25 垂荡运动 纵摇运动 图 15 垂荡和纵摇运动幅值与试验值比较 11 图 16 1 1 3 t= ， ， ,1T 时刻自由面变化（ L =1.25 ） 4 2 4 8.3 自航和静水 Z 形操纵算例 o o 本计算基于 10 /10 Z 形操纵试验要求进行，计算中的螺旋桨推力采用简单的体积力来替代，其 推力大小依据螺旋桨敞水试验下的 J , KT , K Q 结果，其中伴流速度采用船舶速度替代。整个数值模 拟中保持转速不变，在上述转速和未操舵的情况下，船舶先通过自航模拟能达到稳定的航速，然后 在该状态下开始进行操纵运动数值模拟。整个过程网格运动策略为：将舵和船体贴体网格作为整体 考虑，一起做纵荡、横荡、垂荡、横摇、纵摇、首摇运动（为了考虑计算稳定性和计算网格简单， 本文计算中未考虑横摇运动，未考虑舵杆影响）。其中自航模拟中舵不进行控制，Z 形操纵模拟根 据试验规则对舵进行运动控制， 使舵网格绕指定的舵杆转动轴在-1010