邵 菲，楊東梅，劉致豪，李創(chuàng)蘭

（哈爾濱工程大學(xué)船舶工程學(xué)院，哈爾濱150001）

基于MQ插值的變形網(wǎng)格在船舶阻力預(yù)報(bào)中的應(yīng)用

邵 菲，楊東梅，劉致豪，李創(chuàng)蘭

（哈爾濱工程大學(xué)船舶工程學(xué)院，哈爾濱150001）

在非結(jié)構(gòu)切分網(wǎng)格的框架下，發(fā)展了一種基于Multi-Quadric插值的變形網(wǎng)格技術(shù)，結(jié)合CFD軟件STARCCM+，采用SST湍流模型，并以VOF法進(jìn)行自由液面的追蹤，將該變形網(wǎng)格技術(shù)應(yīng)用于帶附體穿浪雙體船的阻力預(yù)報(bào)中。分別以網(wǎng)格的拉伸變形和斜率變形實(shí)現(xiàn)了阻流板和尾楔的邊界運(yùn)動(dòng)，計(jì)算了阻流板高度為4 mm、6 mm和8 mm以及尾楔角度為5°、10°和15°時(shí)的工況。計(jì)算結(jié)果表明，變形網(wǎng)格方法與固定網(wǎng)格方法取得了相同的精度；同試驗(yàn)值相比，在不同工況下變形網(wǎng)格最大平均誤差約為6.67%，驗(yàn)證了其可行性。但是，相比于固定網(wǎng)格方法，變形網(wǎng)格方法的計(jì)算時(shí)耗最大可減少40%，因此采用該變形網(wǎng)格方法可極大地提高帶可變形附體的船舶的阻力預(yù)報(bào)效率。

Multi-quadric插值；網(wǎng)格變形；阻流板；尾楔；計(jì)算時(shí)耗

0 引 言

阻流板與尾楔是常見于船舶艉部的減阻附體，它們通過改善船舶航行過程中的尾部流場(chǎng)來達(dá)到提升船體阻力性能的目的?，F(xiàn)有資料大多采用模型試驗(yàn)或CFD手段來對(duì)這兩種附體的水動(dòng)力特性進(jìn)行研究；其中CFD手段不需要進(jìn)行實(shí)體模型的建造，以及相應(yīng)的時(shí)間、人力和物力的投入，較模型試驗(yàn)具有更好的時(shí)效性，因而得到了廣泛的采用。哈爾濱工程大學(xué)的馬超[1]和鄧銳[2]等人，即采用CFD手段對(duì)加裝于滑行艇和雙體船上的阻流板以及尾壓浪板進(jìn)行了研究。

阻流板或尾楔的形狀對(duì)其水動(dòng)力特性有著至關(guān)重要的影響，在傳統(tǒng)的CFD計(jì)算方案[3-4]中，為了得出較優(yōu)的附體外形，需要針對(duì)不同的附體外形進(jìn)行網(wǎng)格劃分，并重復(fù)由計(jì)算初始化到計(jì)算收斂的全部迭代過程。整個(gè)計(jì)算流程中，需要較多人工干預(yù)以及重復(fù)工作量，效率較低。

而變形網(wǎng)格則可在求解過程中通過改變網(wǎng)格節(jié)點(diǎn)坐標(biāo)來達(dá)到控制邊界形狀的目的，省去了網(wǎng)格的重新劃分、流場(chǎng)求解的初始化以及相應(yīng)的收斂時(shí)耗，因而具有較高的計(jì)算效率。目前的變形網(wǎng)格技術(shù)主要應(yīng)用于飛行器的運(yùn)動(dòng)模擬[5]、機(jī)翼變形（結(jié)冰）[6]、魚體仿生[7]等方面，均取得了較好的模擬效果。在此基礎(chǔ)上本文提出了一種基于multi-quadric（MQ）[8]插值法的網(wǎng)格變形技術(shù)，并將其應(yīng)用于某型穿浪雙體船的阻流板及尾楔的變形控制，從而提高靜水阻力計(jì)算的效率。

1 基于MQ插值的變形網(wǎng)格基本原理

1.1 MQ插值基本原理

Multi-quadric函數(shù)是徑向基函數(shù)的一種，通常簡(jiǎn)記為MQ函數(shù)，MQ函數(shù)在地球物理學(xué)、測(cè)繪學(xué)、數(shù)字地形模擬及水文學(xué)諸方面都有廣泛的應(yīng)用?；贛Q函數(shù)插值問題的描述如下，對(duì)于給定數(shù)據(jù)點(diǎn)集：，尋找函數(shù)

滿足插值條件

而基于MQ插值的變形網(wǎng)格技術(shù)，其基本原理則是通過一系列的控制點(diǎn)來完成變形網(wǎng)格對(duì)網(wǎng)格節(jié)點(diǎn)的初始運(yùn)動(dòng)的控制。每一個(gè)控制點(diǎn)被賦予一個(gè)用來移動(dòng)鄰近網(wǎng)格節(jié)點(diǎn)的位移矢量。這些特定的位移值利用MQ插值法在進(jìn)行變形的整個(gè)區(qū)域內(nèi)進(jìn)行插值，即可得到各網(wǎng)格節(jié)點(diǎn)的位移。

式中：為兩節(jié)點(diǎn)的距離，cj為為一個(gè)基礎(chǔ)常量，本文中該常量取為0。該方程滿足：

1.2 阻流板網(wǎng)格變形策略

在網(wǎng)格變形中，運(yùn)動(dòng)邊界（變形附體）上的每一個(gè)控制點(diǎn)都對(duì)應(yīng)一個(gè)特定的位移矢量，該位移矢量將決定每一個(gè)時(shí)間步長(zhǎng)內(nèi)該控制點(diǎn)的移動(dòng)距離和方向，將該位移轉(zhuǎn)化為各節(jié)點(diǎn)的運(yùn)動(dòng)速度，即可實(shí)現(xiàn)邊界的變形。

對(duì)阻流板來說，其在工作時(shí)一般是沿垂直于船底的方向運(yùn)動(dòng)，其運(yùn)動(dòng)形式可以看做是變形邊界沿一個(gè)或多個(gè)方向的延伸，即拉伸變形。如圖1所示，經(jīng)過時(shí)間T，幾何體的邊界L1、L2、L3拉伸變形為L(zhǎng)1′、L2′、L3′，其變形量為Δs，相應(yīng)節(jié)點(diǎn)的運(yùn)動(dòng)速度為：

式中：θ為變形邊界與坐標(biāo)軸的夾角。

1.3 尾楔網(wǎng)格變形策略

尾楔的運(yùn)動(dòng)形式可以看作是運(yùn)動(dòng)邊界與固定邊界之間的夾角發(fā)生改變，即運(yùn)動(dòng)邊界進(jìn)行斜率變形。如圖2所示，變形邊界L1進(jìn)行斜率變形，經(jīng)過時(shí)間T后，其與L2之間的夾角發(fā)生了變化，而L3則相應(yīng)地做拉伸變形，拉伸變形量為Δs。各控制點(diǎn)可認(rèn)為僅沿x方向運(yùn)動(dòng)，其速度為：

式中：YG和YQ分別為頂點(diǎn)G和Q在y方向上的坐標(biāo)值，可以看出控制點(diǎn)的移動(dòng)速度為一個(gè)與控制點(diǎn)y方向位置有關(guān)的線性函數(shù)。

圖1 拉伸變形示意圖Fig.1 Schematic diagram of tensile deformation

圖2 斜率變形示意圖Fig.2 Schematic diagram of slope deformation

2 數(shù)值船池的建立

2.1 數(shù)值計(jì)算方法

本文采用通用CFD軟件STAR-CCM、利用RANS（Reynolds-Averaged Navier-Stokes）方程方法對(duì)一穿浪雙體船的黏性繞流場(chǎng)進(jìn)行模擬，所求解的連續(xù)性方程和動(dòng)量方程如下：

式中：ui，uj為速度分量時(shí)均值；p為壓力時(shí)均值；ρ為流體密度；μ為動(dòng)力粘性系數(shù)；ρui′uj′為雷諾應(yīng)力項(xiàng)。選取SST（Shear Stress Transport）湍流模型來封閉RANS方程，并以VOF（Volume of fraction）法來實(shí)現(xiàn)流域內(nèi)水氣兩相流的劃分和自由液面的追蹤。

2.2 計(jì)算模型

本文中所研究的對(duì)象為一穿浪雙體船，其外形如圖3所示；阻流板和尾楔均安裝于艇艉，且寬度與片體寬度相同，其安裝形式如圖4所示。表1中則給出了模型的主要尺度參數(shù)。

2.3 網(wǎng)格劃分及邊界條件

計(jì)算域的尺寸對(duì)計(jì)算的收斂時(shí)間和精度都有很大的影響，本文中考慮到雙體船模型的對(duì)稱性，計(jì)算域只針對(duì)一側(cè)的片體進(jìn)行建立；其入口處距船艏1倍船長(zhǎng)，出口處距船尾3倍船長(zhǎng)，上下邊界分別距船體0.8倍船長(zhǎng)和1倍船長(zhǎng)，自船體向外側(cè)延伸1.5倍船長(zhǎng)。邊界條件的設(shè)置如下：

入口：指定來流速度，即為船模試驗(yàn)中各速度點(diǎn)所對(duì)應(yīng)的船速；

對(duì)稱面：采用對(duì)稱邊界條件；

出口：出口處指定壓力分布為靜壓；

船體：不可滑移壁面；

其他壁面：滑移壁面。

本文采用切分網(wǎng)格對(duì)整個(gè)計(jì)算域進(jìn)行離散，其中船體表面網(wǎng)格尺寸取船長(zhǎng)的5‰，最大體網(wǎng)格尺寸則為5%船長(zhǎng)。在船體近壁面和自由液面處進(jìn)行了網(wǎng)格的加密；其中船體近壁面網(wǎng)格的無量綱厚度y+的取值范圍為20～150，自由液面體網(wǎng)格尺寸為1%船長(zhǎng)，總網(wǎng)格數(shù)量約為110萬。計(jì)算域的網(wǎng)格劃分及阻流板和尾楔處的網(wǎng)格形式如圖5所示。

圖3 穿浪雙體船計(jì)算模型Fig.3 Calculation model of wave piercing catamaran

圖4 阻流板及尾楔的安裝形式Fig.4 Installation forms of interceptors and stern wedge

表1 計(jì)算模型主要尺度參數(shù)Tab.1 The principal parameters of calculation model

圖5 計(jì)算域及阻流板和尾楔上的網(wǎng)格劃分Fig.5 Mesh generation in interceptors,stern wedge and calculation domain

3 計(jì)算結(jié)果及分析

3.1 帶阻流板模型的計(jì)算

基于拉伸變形的不同高度阻流板下船體阻力計(jì)算的基本思想是：只建立一個(gè)帶有初始高度的阻流板的船體模型，對(duì)該高度阻流板下的模型進(jìn)行計(jì)算，計(jì)算完成后在當(dāng)前流場(chǎng)環(huán)境下，利用拉伸變形將阻流板高度改變，進(jìn)行第二個(gè)工況下的計(jì)算。以此類推，該航速下阻流板不同高度的工況由一次計(jì)算完成。

圖6 利用MQ方法計(jì)算不同高度阻流板時(shí)的收斂時(shí)歷曲線Fig.6 Convergence time history curves of calculation results by using MQ method in different height of interceptors

圖7 阻流板網(wǎng)格拉伸前后對(duì)比Fig.7 Comparison of grid changes before and after extension in interceptors

圖8 兩種方法在計(jì)算不同阻流板高度模型時(shí)的精度對(duì)比Fig.8 Accuracy comparison of calculation in different heights of interceptors by two different methods

采用MQ方法計(jì)算的收斂歷程如圖6所示，初始給定阻流板高度為4 mm的計(jì)算模型，對(duì)其進(jìn)行網(wǎng)格劃分，計(jì)算至第一次穩(wěn)定收斂時(shí)的物理時(shí)間為5 s；此時(shí)進(jìn)行網(wǎng)格拉伸，如圖7所示，經(jīng)過0.1 s后，阻流板高度變?yōu)? mm；再次計(jì)算至8 s時(shí)，計(jì)算第二次收斂；再次拉伸阻流板網(wǎng)格，阻流板高度變?yōu)? mm；最終在11 s時(shí)，完成計(jì)算的第三次收斂。圖6中同時(shí)給出了固定網(wǎng)格方法單次計(jì)算的收斂歷程，可以看出，采用固定網(wǎng)格方法的單次計(jì)算收斂時(shí)間同樣也為5 s，但之后的計(jì)算需要重復(fù)網(wǎng)格劃分和相同的計(jì)算過程，因此采用MQ方法后整體的計(jì)算效率有了明顯提高。

圖8中給出了在航速為4.84 m/s的航速下，兩種方法所計(jì)算的不同阻流板高度模型的阻力與試驗(yàn)值的對(duì)比，可以看出，MQ方法在網(wǎng)格變形之后仍能取得與固定網(wǎng)格法相當(dāng)?shù)挠?jì)算精度；在4 mm、6 mm和8 mm高度的阻流板的工況中，兩種方法的計(jì)算誤差分別為5.86%和5.86%、5.84%和5.18%以及5.79%和5.55%。而在不同航速下MQ方法計(jì)算值與試驗(yàn)值的對(duì)比也表明該方法對(duì)應(yīng)于不同工況均具有較好的可行性，如圖9所示，在4 mm、6 mm和8 mm高度的阻流板的工況中，其平均誤差分別為6.22%、6.67%和6.20%。

圖9 不同航速下MQ方法計(jì)算值與試驗(yàn)值的對(duì)比Fig.9 Comparison between calculated values by using MQ method and experimental values in different speeds

3.2 帶尾楔模型的計(jì)算

帶尾楔模型的計(jì)算思想與帶阻流板模型的基本相同，即在計(jì)算收斂之后利用網(wǎng)格的斜率變形改變尾楔的角度。其收斂時(shí)歷如圖10所示，可以看出經(jīng)過兩次歷時(shí)0.1 s的斜率變形之后（如圖11），根據(jù)由（7）式所計(jì)算的節(jié)點(diǎn)變形速度，尾楔角度由5°變?yōu)?5°，歷次計(jì)算均達(dá)到收斂要求，總的計(jì)算時(shí)間仍為11 s，計(jì)算效率要優(yōu)于固定網(wǎng)格法。

圖10 利用MQ方法計(jì)算不同尾楔角度模型時(shí)的收斂時(shí)歷曲線Fig.10 Convergence time history curves of calculation results by using MQ method in different angles of stern wedge

圖11 尾楔網(wǎng)格斜率變形前后對(duì)比Fig.11 Comparison of grid changes before and after slope deformation in stern wedge

兩種方法計(jì)算精度的對(duì)比如圖12所示，圖13則給出了MQ方法計(jì)算值在不同工況下與試驗(yàn)值的對(duì)比，可以看出在計(jì)算帶尾楔模型的阻力時(shí)，MQ方法在計(jì)算精度上與固定網(wǎng)格法相差依然不大；而在5°、10°和15°尾楔的工況中，其相對(duì)于試驗(yàn)值的平均誤差分別為3.97%、4.17%和6.25%，這表明該方法在尾楔變形的計(jì)算中仍然具有較好的可行性。

圖12 兩種方法在計(jì)算不同尾楔角度模型時(shí)的精度對(duì)比Fig.12 Accuracy comparison of calculation in different angles of stern wedge by two different methods

圖13 不同航速下MQ方法計(jì)算值與試驗(yàn)值的對(duì)比Fig.13 Comparison between calculated values by using MQ method and experimental values in different speeds

3.3 MQ方法計(jì)算效率分析

根據(jù)以上計(jì)算結(jié)果可以看出，MQ方法在初始網(wǎng)格的計(jì)算中同固定網(wǎng)格法相同，達(dá)到收斂的物理時(shí)間為5 s，而網(wǎng)格變形后達(dá)到收斂的物理時(shí)間均為3 s。在相同的計(jì)算條件下，收斂的物理時(shí)間即可反映出計(jì)算實(shí)際的時(shí)間消耗，因此可得出MQ方法相對(duì)于固定網(wǎng)格法的計(jì)算時(shí)耗收益δ為：式中：M為計(jì)算所需考慮的速度點(diǎn)的個(gè)數(shù)，N則表示附體方案的個(gè)數(shù)?？梢钥闯觯襟w方案?jìng)€(gè)數(shù)越多，相應(yīng)的時(shí)耗收益也就越大，在本文判定收斂的物理時(shí)間的基礎(chǔ)上，其最大值為40%。

4 結(jié) 論

本文將基于MQ插值的變形網(wǎng)格應(yīng)用于帶阻流板或尾楔的穿浪雙體船的阻力計(jì)算中，采用MQ函數(shù)構(gòu)建了各網(wǎng)格節(jié)點(diǎn)的位移插值函數(shù)，并以拉伸變形和斜率變形這兩種網(wǎng)格變形方式實(shí)現(xiàn)了動(dòng)邊界的運(yùn)動(dòng)，根據(jù)本文的分析可得出以下幾點(diǎn)結(jié)論：

（1）基于MQ插值的變形網(wǎng)格取得了與固定網(wǎng)格法相當(dāng)?shù)木?，與試驗(yàn)值的對(duì)比也表明該方法在不同工況下均有較好的可行性。

（2）基于MQ插值的變形網(wǎng)格可有效提高多附體方案CFD計(jì)算的效率，相對(duì)于固定網(wǎng)格法，其最大時(shí)耗收益為40%。

（3）除阻流板和尾楔之外，根據(jù)基于MQ插值的變形網(wǎng)格中邊界的運(yùn)動(dòng)特性，該方法還可應(yīng)用于T型水翼浸濕調(diào)整、多體船片體間距改變、尾壓浪板和壓浪條下壓角度的變化等附體變形的計(jì)算中。

[1]馬 超.阻流板和尾壓浪板對(duì)滑行艇阻力性能影響[D].哈爾濱:哈爾濱工程大學(xué),2012. Ma Chao.The effect of interceptor and stern flap on the resistance of planning boat[D].Harbin:Harbin Engineering University,2012.

[2]鄧 銳.阻流板對(duì)雙體船水動(dòng)力性能影響的數(shù)值研究[D].哈爾濱:哈爾濱工程大學(xué),2010. Deng Rui.Numerical research on influence of the interceptor on catamaran hydrodynamic performances[D].Harbin:Harbin Engineering University,2010.

[3]鄧 銳,黃德波,周廣利,等.帶有阻流板的雙體船水動(dòng)力性能數(shù)值研究[J].華中科技大學(xué)學(xué)報(bào)（自然科學(xué)版）,2011 (4):97-110. Deng Rui,Huang Debo,Zhou Guangli,et al.Numerical research on hydrodynamic performance of catamarans with interceptors[J].Journal of Huazhong University of Science and Technology(Natural Science Edition),2011(4):97-110.

[4]鄧 銳,黃德波,周廣利,等.阻流板水動(dòng)力機(jī)理的初步計(jì)算研究[J].船舶力學(xué),2012,16(7):740-749. Deng Rui,Huang Debo,Zhou Guangli,et al.Preliminary numerical research of the hydrodynamic mechanism of interceptor[J].Journal of Ship Mechanics,2012,16(7):740-749.

[5]黃守智,李 杰,蔣勝矩,等.基于變形網(wǎng)格技術(shù)的非定常流動(dòng)數(shù)值分析方法研究[J].導(dǎo)箭與制導(dǎo)學(xué)報(bào),2005,25(3): 186-188. Huang Shouzhi,Li Jie,Jiang Shengju,et al.Unsteady viscous flow simulations with deforming grid[J].Journal of Projectiles,Rockets,Missiles and Guidance,2005,25(3):186-188.

[6]馮文梁,李 杰,張 威.基于變形網(wǎng)格技術(shù)的翼型結(jié)冰數(shù)值模擬研究[J].西北工業(yè)大學(xué)學(xué)報(bào),2008,26(5):550-555. Feng Wenliang,Li Jie,Zhang Wei.Exploring numerical simulation of icing flow-field of airfoil based on grid deformation technique[J].Journal of Northwestern Polytechnical University,2008,26(5):550-555.

[7]張來平,段旭鵬,常興華,等.基于Delaunay背景網(wǎng)格插值和局部網(wǎng)格重構(gòu)的變形體動(dòng)態(tài)混合網(wǎng)格生成技術(shù)[J].空氣動(dòng)力學(xué)學(xué)報(bào),2009,27(1):32-39. Zhang Laiping,Duan Xupeng,Chang Xinghua,et al.A hybrid dynamic grid generation technique for morphing bodies based on Delaunay graph and local remeshing[J].ACTA Aerodynamica Sinica,2009,27(1):32-39.

[8]Hardy R L.Theory and applications of the multiquadric-biharmonic method[J].Comput Math Applic,1990,19:163-208.

[9]李 艷,白玉峰.Multi-Quadric函數(shù)與Gauss函數(shù)的插值比較[J].內(nèi)蒙古民族大學(xué)學(xué)報(bào)（自然科學(xué)版）,2010(4):369-372. Li Yan,Bai Yufeng.The case study about multi-quadric method[J].Journal of Inner Mongolia University for Nationalities, 2010(4):369-372.

Application of deforming mesh based on MQ interpolation in the ship resistance calculation

SHAO Fei,YANG Dong-mei,LIU Zhi-hao,LI Chuang-lan
(College of Shipbuilding Engineering,Harbin Engineering University,Harbin 150001,China)

A mesh deformation technique based on multi-quadric interpolation is developed in the framework of unstructured trim mesh.To apply this method to the resistance calculation of a Wave Piercing Catamaran with appendage,the CFD software STAR-CCM is used by coupling the SST turbulence model and the free surface profile is tracked by VOF method.The stretching deformation and slope deformation are used to implement the boundary motions of interceptor and stern wedge,respectively.Cases with the interceptor height of 4 mm,6 mm,8 mm and stern wedge angles of 5°,10°,15°are calculated.The numerical results show that the mesh deformation method has the same degree of accuracy comparing with the mesh fixation method.Validation of this method is carried out by comparison with the experimental results;the maximum average error is about 6.67%.Comparison of time consumption shows that the mesh deformation methods get a maximum computational time reduction of 40%to the mesh fixation method.Thereby,Using this method to predict the resistance of ships with appendage can greatly increase the computation efficiency.

multi-quadric interpolation;mesh deformation;interceptor;stern wedge;computational time consumption

O35

：Adoi:10.3969/j.issn.1007-7294.2016.04.009

1007-7294（2016）04-0452-08

2016-03-17

國(guó)家自然科學(xué)基金資助（51509053）

邵 菲（1987-），女，博士研究生；楊東梅（1979-），女，博士，通訊作者，E-mail:yangdm411@126.com。