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

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

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

基于MQ插值的變形網(wǎng)格在船舶阻力預(yù)報中的應(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ù)報中。分別以網(wǎng)格的拉伸變形和斜率變形實現(xiàn)了阻流板和尾楔的邊界運動，計算了阻流板高度為4 mm、6 mm和8 mm以及尾楔角度為5°、10°和15°時的工況。計算結(jié)果表明，變形網(wǎng)格方法與固定網(wǎng)格方法取得了相同的精度；同試驗值相比，在不同工況下變形網(wǎng)格最大平均誤差約為6.67%，驗證了其可行性。但是，相比于固定網(wǎng)格方法，變形網(wǎng)格方法的計算時耗最大可減少40%，因此采用該變形網(wǎng)格方法可極大地提高帶可變形附體的船舶的阻力預(yù)報效率。

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

0 引 言

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

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

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

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

1.1 MQ插值基本原理

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

滿足插值條件

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

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

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

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

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

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

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

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

式中：YG和YQ分別為頂點G和Q在y方向上的坐標(biāo)值，可以看出控制點的移動速度為一個與控制點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ù)值計算方法

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

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

2.2 計算模型

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

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

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

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

對稱面：采用對稱邊界條件；

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

船體：不可滑移壁面；

其他壁面：滑移壁面。

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

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

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

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

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

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

3.1 帶阻流板模型的計算

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

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

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

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

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

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

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

3.2 帶尾楔模型的計算

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

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

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

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

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

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

3.3 MQ方法計算效率分析

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

4 結(jié) 論

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

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

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

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

[1]馬 超.阻流板和尾壓浪板對滑行艇阻力性能影響[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]鄧 銳.阻流板對雙體船水動力性能影響的數(shù)值研究[D].哈爾濱:哈爾濱工程大學(xué),2010. Deng Rui.Numerical research on influence of the interceptor on catamaran hydrodynamic performances[D].Harbin:Harbin Engineering University,2010.

[3]鄧 銳,黃德波,周廣利,等.帶有阻流板的雙體船水動力性能數(shù)值研究[J].華中科技大學(xué)學(xué)報（自然科學(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]鄧 銳,黃德波,周廣利,等.阻流板水動力機(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ù)的非定常流動數(shù)值分析方法研究[J].導(dǎo)箭與制導(dǎo)學(xué)報,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é)報,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)的變形體動態(tài)混合網(wǎng)格生成技術(shù)[J].空氣動力學(xué)學(xué)報,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é)報（自然科學(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

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

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