基于MQ插值的變形網格在船舶阻力預報中的應用

邵 菲，楊東梅，劉致豪，李創蘭

（哈爾濱工程大學船舶工程學院，哈爾濱150001）

基于MQ插值的變形網格在船舶阻力預報中的應用

邵 菲，楊東梅，劉致豪，李創蘭

（哈爾濱工程大學船舶工程學院，哈爾濱150001）

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

Multi-quadric插值；網格變形；阻流板；尾楔；計算時耗

0 引 言

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

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

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

1 基于MQ插值的變形網格基本原理

1.1 MQ插值基本原理

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

滿足插值條件

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

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

1.2 阻流板網格變形策略

在網格變形中，運動邊界（變形附體）上的每一個控制點都對應一個特定的位移矢量，該位移矢量將決定每一個時間步長內該控制點的移動距離和方向，將該位移轉化為各節點的運動速度，即可實現邊界的變形。

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

式中：θ為變形邊界與坐標軸的夾角。

1.3 尾楔網格變形策略

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

式中：YG和YQ分別為頂點G和Q在y方向上的坐標值，可以看出控制點的移動速度為一個與控制點y方向位置有關的線性函數。

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

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

2 數值船池的建立

2.1 數值計算方法

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

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

2.2 計算模型

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

2.3 網格劃分及邊界條件

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

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

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

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

船體：不可滑移壁面；

其他壁面：滑移壁面。

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

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

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

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

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

3 計算結果及分析

3.1 帶阻流板模型的計算

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

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

圖7 阻流板網格拉伸前后對比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的計算模型，對其進行網格劃分，計算至第一次穩定收斂時的物理時間為5 s；此時進行網格拉伸，如圖7所示，經過0.1 s后，阻流板高度變為6 mm；再次計算至8 s時，計算第二次收斂；再次拉伸阻流板網格，阻流板高度變為8 mm；最終在11 s時，完成計算的第三次收斂。圖6中同時給出了固定網格方法單次計算的收斂歷程，可以看出，采用固定網格方法的單次計算收斂時間同樣也為5 s，但之后的計算需要重復網格劃分和相同的計算過程，因此采用MQ方法后整體的計算效率有了明顯提高。

圖8中給出了在航速為4.84 m/s的航速下，兩種方法所計算的不同阻流板高度模型的阻力與試驗值的對比，可以看出，MQ方法在網格變形之后仍能取得與固定網格法相當的計算精度；在4 mm、6 mm和8 mm高度的阻流板的工況中，兩種方法的計算誤差分別為5.86%和5.86%、5.84%和5.18%以及5.79%和5.55%。而在不同航速下MQ方法計算值與試驗值的對比也表明該方法對應于不同工況均具有較好的可行性，如圖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 帶尾楔模型的計算

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

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

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

兩種方法計算精度的對比如圖12所示，圖13則給出了MQ方法計算值在不同工況下與試驗值的對比，可以看出在計算帶尾楔模型的阻力時，MQ方法在計算精度上與固定網格法相差依然不大；而在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方法計算效率分析

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

4 結 論

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

（1）基于MQ插值的變形網格取得了與固定網格法相當的精度，與試驗值的對比也表明該方法在不同工況下均有較好的可行性。

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

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

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

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

[3]鄧 銳,黃德波,周廣利,等.帶有阻流板的雙體船水動力性能數值研究[J].華中科技大學學報（自然科學版）,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].船舶力學,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]黃守智,李 杰,蔣勝矩,等.基于變形網格技術的非定常流動數值分析方法研究[J].導箭與制導學報,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]馮文梁,李 杰,張 威.基于變形網格技術的翼型結冰數值模擬研究[J].西北工業大學學報,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背景網格插值和局部網格重構的變形體動態混合網格生成技術[J].空氣動力學學報,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函數與Gauss函數的插值比較[J].內蒙古民族大學學報（自然科學版）,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

國家自然科學基金資助（51509053）

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