基于重疊網格的船舶運動計算方法研究

徐偉光，趙發明，何術龍

（中國船舶科學研究中心，江蘇 無錫 214082）

基于重疊網格的船舶運動計算方法研究

徐偉光，趙發明，何術龍

（中國船舶科學研究中心，江蘇 無錫 214082）

文章應用重疊網格技術和單相Level-set捕捉自由面方法，結合數值造波技術和六自由度運動方程，形成了基于RANS方程的船舶運動響應數值計算方法，并在中國船舶科學研究中心自主研發的船舶專用粘流計算軟件OShip上實現。應用該方法計算了DTMB5512船模在迎浪規則波中的運動響應，并通過對計算得到的船模運動時歷曲線進行傅里葉變換，提取其運動的幅值和相位，得到了不同航速下船模垂蕩和縱搖運動的頻響曲線，所得結果與試驗結果吻合較好。

重疊網格；六自由度；傳遞函數和相位

0 引 言

船舶的耐波性作為衡量現代化船舶航行性能的重要指標之一，對船舶的安全性、作戰或作業使用性以及適居性都會產生巨大的影響。船舶在波浪中航行的時候，船舶航行的姿態將會嚴重影響到船舶的阻力性能以及螺旋槳的推進效率，大幅度的運動會致使螺旋槳露出水面，造成飛車現象，同時甲板上浪將影響到船上設備的正常工作，大幅的搖蕩還會大大降低對船員和乘客的舒適性。對此，如何準確地預報船舶的耐波性能將是本文研究的重要任務之一。

傳統的耐波性預報方法主要基于勢流理論，由于其簡單、計算效率高的優點得到了廣泛的應用，但在處理強非線性問題上還存在一定的局限性，如甲板上浪、自由面破碎等現象，且這些問題都會對船舶的運動產生較大的影響。近些年，隨著計算機技術的飛速發展，CFD方法得到了越來越多的應用，由于考慮了流體的粘性效應和非線性因素，CFD方法在準確預報船舶耐波性上體現出了巨大的優勢。

國內外已經開展了許多關于船舶運動的粘性數值模擬，Sat等人[2]采用密度函數法對規則波中W-igley船型和S60船型的運動響應進行了數值模擬；Carrica和Wilson等人[3]采用重疊網格技術模擬計算了DTMB 5512船型在中高航速下的大幅度運動響應情況；吳乘勝[4]對基于Rans方程的數值波浪水池技術展開了初步研究，其中對wigley及DTMB5512波浪運動響應進行了預報。

本文將采用Rans方法和重疊網格[5]，對DTMB5512船模在不同航速下，頂浪規則波中的運動響應進行預報，考察船模在較大的波長范圍內（0.5L＜λ＜2.5L）的運動特性。本文采用的求解器為中國船舶科學研究中心自主研發的船舶專用粘流計算軟件OShip。

1 數學模型及數值方法

1.1 控制方程與湍流模型

在求解不可壓流體的非定常運動時，控制方程采用RANS方程，湍流模型為SST型k-ω二方程模型，并采用PISO[6]算法進行RANS方程和連續性方程的求解。

1.2 重疊網格方法

OShip軟件采用重疊網格方法來處理運動問題。重疊網格是允許各個網格區域相互重疊、嵌套，并在各個子網格間通過插值來傳遞流場信息的一種結構網格方法。它將復雜的流動區域分成幾何邊界比較簡單的子區域，極大簡化了復雜形狀物體網格的生成，提高了結構網格對外形的適應能力。重疊網格生成分為以下三個步驟：（1）尋點：在網格空間對已知點或區域進行搜索，得到與之相關的網格單元。（2）挖洞：在給定的網格中，將落在指定區域中的網格單元從參與流場計算的網格單元集合中剔除。（3）插值：建立各重疊網格間的耦合關系，由相鄰貢獻單元獲得待插值點信息。

圖1 重疊網格示意圖Fig.1 Demonstration of overset grid

如圖1所示，根據貢獻單元的位置與插值點的相對位置，可以求得貢獻點與待插值點間的插值系數[7]。同時根據這些插值系數，我們可以得到被插值點的任意變量值φ：

式中：ξi，ηi，ζi為插值系數，φi為貢獻點上的變量值。

1.3 自由面處理

本文采用單相Level-set方法[8]進行自由面的捕捉，自由面的位置通過Level-set函數的函數值φ確定。單相Level-set方法僅需求解在距離函數φ≤0（表示水相）區域內的流場，空氣相中的流場速度則通過速度擴展（velocity extension）的方法來計算的。兩相流界面的過渡問題通過顯性施加階躍條件進行處理。此外，在氣體中，只需要布置少許網格來滿足計算條件，因此相比兩相方法，計算資源的消耗大大減小，計算穩定性增加。

1.4 造波方法

本文的計算工況均為線性規則波，可通過定義初始邊界條件來實現，在任意時刻的某一點的波浪幅值可以按（2）式表達為：

式中：ξ為自由液面上各點在任意時刻的垂向位置，a為波幅，ω為波浪的自然頻率，ωe為遭遇頻率，k為波數。在上述波浪初始條件下，速度和壓力的表達式為：

式中：U（x，y，z，t），W（x，y，z，t）分別為流體中的點在任意時刻水平方向和垂直方向的速度分量，U0為船模的速度。

1.5 六自由度運動

為了求解運動方程，需要建立兩個坐標系，一個是以地面為參考的固定坐標系。另一個是載體坐標系，其原點放在船舶重心處，船尾方向為x軸正方向，右舷方向為y軸正方向，垂直于水線面向上為z軸正方向。

固定坐標系和載體坐標系間的轉化關系參見文獻[3]。載體坐標系下，船舶六自由度運動方程為：

式中：X，Y，Z，K，M和N分別為縱蕩、橫蕩、垂蕩力、橫搖、縱搖及艏搖力矩。為線加速度，為歐拉角加速度。六自由度運動求解的流程如圖3所示，運動網格重疊結果如圖2所示。從圖4中可以看出，當船模運動后，重疊網格方法僅需重新計算背景網格即與船模貼體網格的插值關系，而不用對網格重新生成，這樣既保證了船體網格的質量保持不變同時也保證了自由面處網格的密度，從而提高了運動求的精度。

圖2 固定坐標系和載體坐標系Fig.2 Earth and ship fixed reference systems

圖3 六自由度運動求解流程圖Fig.3 Solution strategy

圖4 動網格重疊示意圖Fig.4 Demonstration of dynamic overset grid

1.6 離散方法

控制方程采用體積中心有限差分格式進行離散，時間項采用2階歐拉向后差分格式、對流項采用2階混合差分格式，粘流項采用2階中心差分格式離散。Level-set對流項同樣采用2階混合差分格式。

2 數值計算

2.1 計算對象

本文選用DTMB5512船模作為計算對象。DTMB5512是ITTC推薦的標準模型之一，具有大量的試驗數據和計算結果。DTMB5512的縮尺比為1/46.6。

2.2 計算網格

本文的算例中的船體和計算域左右對稱，因此在生成網格時只需生成其中的一半。計算區域如圖5所示，計算域尺度如下：

圖5 計算域網格示意圖Fig.5 Computational grid

a.入口—模型首部前0.6 Lpp處；b.出口—模型為尾部1.5 Lpp處；c.側邊界—模型側方1.0 Lpp處；d.上邊界—水線以上0.25 Lpp處；e.下邊界—水線以下1.0 Lpp處。

2.3 網格生成

為了能夠準確計算船舶在波浪中的運動，在劃分背景網格時，保證一個船長內至少含有60個網格。整個計算區域的網格總數為1 224 000，計算網格如圖6所示。

圖6 船首、尾網格Fig.6 Local grid

2.4 計算工況

數值模擬中的各計算工況的參數均列于表1中。

表1 DTMB5512耐波性計算工況Tab.1 Summary of simulation conditions

3 計算結果及分析

3.1 數據處理方法

船模在頂浪規則波中作垂蕩和縱搖的耦合運動，且船模的運動和受力具有周期性，因此將船模的運動用傅里葉級數展開可以寫為：

式中：x3n，x5n和γ3n，γ5n分別為船模垂蕩和縱搖運動的n階幅值和相位，而垂蕩和縱搖運動響應的傳遞函數可表達為（1階運動）：

式中：γ31，γ51為傳遞函數的相位角。

船模試驗中規定：t=0時規則波的波峰正好經過船模重心的縱向位置LCG，而本文數值計算中t=0時，LCG處未必對應著波峰，因此在計算相位時還需進行修正：

式中：γζLCG為t=0時LCG處的波浪相位，γζ0為t=0時入口處的波浪相位，s為入口至LCG的距離，γ3，γ5則為修正后的相位角。

3.2 傳遞函數及運動響應相位

在幅值和相位的處理中，采用船模運動穩定后有效的數據進行分析，采用快速傅里葉變換進行離散數據的處理從而獲得船模運動的1階幅值和對應的相位，將由FFT獲得的垂蕩和縱搖運動的1階幅值代入（7）式和（8）式中，得到垂蕩和縱搖運動的傳遞函數及運動響應相位，圖7-14給出了本文計算值與試驗值和切片法計算值（SMP）的比較結果。

圖7 Fr=0.28垂蕩運動傳遞函數Fig.7 Heave transfer functions at Fr=0.28

圖8 Fr=0.28縱搖運動傳遞函數Fig.8 Pitch transfer functions at Fr=0.28

圖9 Fr=0.28垂蕩運動響應相位Fig.9 Heave phases at Fr=0.28

圖10 Fr=0.28縱搖運動響應相位Fig.10 Pitch phases at Fr=0.28

圖11 Fr=0.41垂蕩運動傳遞函數Fig.11 Heave transfer functions at Fr=0.41

圖12 Fr=0.41縱搖運動傳遞函數Fig.12 Pitch transfer functions at Fr=0.41

圖13 Fr=0.41垂蕩運動響應相位Fig.13 Heave phases at Fr=0.41

圖14 Fr=0.41縱搖運動響應相位Fig.14 Pitch phases at Fr=0.41

圖15 垂蕩運動傳遞函數隨遭遇頻率變化Fig.15 Heave transfer functions for all Fr versus fe

圖16 縱搖運動傳遞函數隨遭遇頻率變化Fig.16 Pitch transfer functions for all Fr versus fe

由圖7-10的結果可得出如下結論：

（1）當Fr=0.28時，本文的計算值與試驗值符合較好，計算所得船模運動隨波長變化的趨勢與試驗值一致，本文計算結果優于切片法（SMP）的計算值。

（2）當λ增大時，垂蕩運動的傳遞函數值也隨之增大，在λ/L=1.3附近垂蕩的傳遞函數值到達峰值點。隨著波長進一步增加，傳遞函數值有所下降，而當λ/L＞1.5時，隨著波長的增加，傳遞函數值再次平穩地上升。同時，縱搖運動的傳遞函數值在λ/L=1.5附近達到峰值點，之后隨著波長的增加而減小最終趨于平緩。另一方面，從計算結果來看，峰值處傳遞函數計算值要略小于試驗值，但峰值點出現位置捕捉得較為準確，相比之下切片法所計算的峰值位置向右偏移較為嚴重，且峰值也存在較大誤差。

（3）船模的垂蕩和縱搖運動存在著一定的相位差，γx3要超前γx5，圖9、10中，隨著波長的增加，γx3、γx5也隨之減小，并分別在λ/L=0.85和0.7附近達到最小值，之后隨著波長的增加而增加，最后保持不變。最終γx3=0表明船模的運動與波浪同步。計算結果表明，在峰值點附近本文的相位計算值要稍大于試驗值，即相位要超前一些，但從整體趨勢上來看二者符合得較好，而切片法的計算值在峰值附近要偏小很多。

圖11-14為Fr=0.41時的計算結果，可以看出：

（1）垂蕩和縱搖運動傳遞函數分別在λ/L=1.4和λ/L=1.6處達到峰值點，之后傳遞函數值隨著波長的增加而減小。本文的計算值與試驗值符合得較好，準確地捕捉到了峰值點的位置，同時曲線的趨勢與試驗值一致，但在峰值附近本文的計算值要偏小一些。相比之下切片法的計算值誤差則較大，且其計算值的峰值點右移。

（2）垂蕩和縱搖運動的相位γx3和γx5的谷值點分別出現在λ/L=1.0以及λ/L=0.8附近。本文相位計算值在波長較大的區域內與試驗值符合較好，在短波范圍內則稍差，且在1.0＜λ/L＜1.5范圍內計算值稍微偏大，但整體趨勢與試驗值符合得較好。

圖15和16給出了運動計算結果：當Fr=0.28時，垂蕩峰值點出現在λ/L=1.3附近，此時fe=1.014；當Fr=0.41時，λ/L=1.4，此時fe=1.130，兩個峰值處的遭遇頻率并不相同，這表明對于DTMB 5512船模來說，某一航速下的峰值點不一定位于固有頻率處，即發生共振并非對應著最大運動幅值。

根據Lewis，Journee，Fonseca和Soares等人[9-10]的結論，對于絕大多數情況，某一航速下，垂蕩和縱搖運動的峰值發生在L/λ=0.75（λ/L=1.33）附近，此時波浪擾動力最大，進一步當遭遇頻率與船模的固有頻率相等即fe=fn，兩個條件同時滿足時，船舶將產生最大運動響應。fn計算如（12）式：

式中：ω3由Lloyd（1989）按照（13）式給出：

式中：C33為垂蕩恢復力，A33為垂蕩附加質量。由于Fr=0.41時的遭遇頻率fe=1.130更加接近船模的固有頻率fn，因此該航速下的傳遞函數曲線更加陡峭，幅值也更大，本文的計算結果也很好地驗證了這一點。

4 結 論

本文基于重疊網格技術，對波浪中DTMB 5512船模的運動響應進行了數值模擬，并將計算結果與試驗值和切片法（SMP）的計算值進行了比較和驗證，得到以下結論：

（1）本文采用的基于重疊網格的數值方法在預報船模運動響應時得到了穩定、收斂的計算結果，證明了該方法的穩定性。

（2）本文計算了兩種航速Fr=0.28和Fr=0.41的運動響應，得到的傳遞函數和相位與試驗值都符合得較好，能夠較為準確地捕捉到峰值點位置，但計算值在峰值點附近偏小，而切片法（SMP）的計算值誤差則較大。

（3）對于DTMB 5512船模來說，某一航速下的峰值點不一定位于固有頻率處。對于絕大多數情況，某一航速下，垂蕩和縱搖運動的峰值發生在L/λ=0.75（λ/L=1.33）附近，并且當遭遇頻率與船模的固有頻率相等fe=fn時，將產生最大的運動響應。本文的計算結果在一定程度上驗證了這一點，不足的地方將在今后的工作中繼續完善。

綜上，本文的計算結果證明了本文采用的基于重疊網格技術的數值方法在單體船耐波性預報上是適用、穩定以及有效的。該方法在多體船及多體復合船耐波性預報上的推廣將在未來的工作中完成。

[1]周秀紅,趙發明,王麗艷.船舶粘流計算軟件“Oship”開發[J].中國造船,2014,55(1):90-103. Zhou Xiuhong,Zhao Faming,Wang Liyan.Development of software‘Oship’for open/overlap ship viscous caculation[J]. Shipbulding of China,2014,55(1):90-103.

[2]Miyata S H,Sato T.CFD simulation of 3-dimensional motion of a ship in waves:Application to an advancing ship in regular heading waves[J].Journal of Marine Science and Technology,1999,4(3):108-116.

[3]Carrica P M,Wilson R V,Noack R W,et al.Ship motions using single-phase Level-set with dynamic overset grids[J]. Computers&Fluids,2007,36(9):1415-1433.

[4]吳乘勝.基于Rans方程的數值水池研發及其應用研究[D].無錫:中國船舶科學研究中心,2014. Wu Chengsheng.Development of a numerical wave tank based on RANS equations and its application[D].Wuxi:China Ship Science Research Center,2014.

[5]趙發明,高成君,夏 瓊.重疊網格在船舶CFD中的應用研究[J].船舶力學,2011,15(4):332-341. Zhao Faming,Gao Chengjun,Xia Qiong.Overlap grid research on the application of ship CFD[J].Journal of Ship Mechanics,2011,15(4):332-341.

[6]Issa R I.Solution of the implicitly discretised fluid flow equations by operator-splitting[J].Journal of Computational Physics, 1986,62(1):40-65.

[7]Wilson R,Carrica P,Hyma M,Stern F.A single-phase level set method with application to breaking waves and forward speed diffraction problem[C]//25th Symposium on Naval Hydrodynamics St.Canada,2004.

[8]Lewis E V.Principles of naval architecture[J].The Society of Naval Architects and Marine Engineers,New York,1999.

[9]Journee J M J.Experiments and calculations on four wigley hullforms[R].Delft University of Technology,Ship Hydromechanics Lab,Report 1999,No.909.

[10]Fonseac N,Soares C G.Experimental investigation of the onliner effects on vertical motions and loads of a containership in regular wave[J].Journal of Ship Research,2004b,48:118-147.

Research on calculation method for ship motion based on dynamic overset approach

XU Wei-guang,ZHAO Fa-ming,HE Shu-long
(China Ship Scientific Research Center,Wuxi 214082,China)

This research applies dynamic overset grid and single phase Level-set method and accommodates numerical wave making technique and 6 DOF motion equations to form the numerical method for calculating ship motion response,and this method is realized by using self-developed calculation software-Oship developed by China Ship Scientific Research Center.Motion responses of DTMB 5512 advancing in regular head waves are computed using this method.Through the fast Fourier transform of time history of ship motions, amplitudes and phases of ship motions are acquired to determine the heave and pitch transfer functions and phases of DTMB 5512 at different speed.The simulation results show fair agreements with experiments.

overset grid;6 DOF;transfer function and phase

U661.3

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

1007-7294（2016）07-0824-09

2015-08-14

徐偉光（1989-），男，碩士研究生，E-mail：654050813@qq.com；趙發明（1979-），男，高級工程師。