李志富， 任慧龍， 石玉云， 李 輝

（1.江蘇科技大學(xué) 船舶與海洋工程學(xué)院，江蘇 鎮(zhèn)江 212003；2.哈爾濱工程大學(xué) 船舶工程學(xué)院，哈爾濱 150001）

0 引 言

利用滿足自由表面條件的時(shí)域格林函數(shù)求解船舶在波浪中的時(shí)域運(yùn)動(dòng)響應(yīng)時(shí)，需用時(shí)間步進(jìn)法進(jìn)行求解，對(duì)含有時(shí)域Green函數(shù)及其導(dǎo)數(shù)的時(shí)間卷積積分要上億次的計(jì)算時(shí)域Green函數(shù)及其對(duì)空間和時(shí)間的偏導(dǎo)數(shù)，故研究一種高效、精確的時(shí)域格林函數(shù)計(jì)算方法對(duì)于船舶運(yùn)動(dòng)的時(shí)域模擬至關(guān)重要[1-2]。

Wehausen和Laitone[3]給出的三維時(shí)域自由表面Green函數(shù)表達(dá)式至今仍被學(xué)者廣泛采用。由于Green函數(shù)是一個(gè)無窮積分，且被積函數(shù)在自由水面附近具有高頻振蕩、增幅等特點(diǎn)，使得Green函數(shù)及其偏導(dǎo)數(shù)的數(shù)值計(jì)算十分困難，精度難以控制。

Beck和Liapis[4]根據(jù)Green函數(shù)的振蕩特性，將整個(gè)計(jì)算區(qū)域進(jìn)行了劃分，通過在不同的區(qū)域上應(yīng)用級(jí)數(shù)展開和漸進(jìn)展開，對(duì)Green函數(shù)進(jìn)行了計(jì)算。在此基礎(chǔ)上，King[5]增加了一個(gè)額外的小區(qū)間，在此區(qū)間上利用Bessel函數(shù)的展開式進(jìn)行了計(jì)算?；贜ewman的工作，Lin和Yue[6]綜合利用級(jí)數(shù)展開、漸進(jìn)展開和二維多項(xiàng)式近似對(duì)Green函數(shù)進(jìn)行了計(jì)算。黃德波[7]通過采用雙參數(shù)制表插值技術(shù)，提高了Green函數(shù)的計(jì)算效率，但是計(jì)算精度卻有所損失。

通過對(duì)格林函數(shù)及其導(dǎo)數(shù)的積分表達(dá)式進(jìn)行變形，Clement[8]得到了時(shí)域Green函數(shù)所滿足的常微分方程，并通過四階Runge-Kutta方法對(duì)其進(jìn)行了求解，然而計(jì)算精度隨著時(shí)間的增長，逐漸變差。Chuang[9]提出一種半解析的方法求解時(shí)域Green函數(shù)滿足的常微分方程，精度和穩(wěn)定性有了顯著提高，但是涉及到的級(jí)數(shù)求和項(xiàng)數(shù)多達(dá)40余項(xiàng)，無疑增加了計(jì)算時(shí)間。

本文通過對(duì)Green函數(shù)及其導(dǎo)數(shù)所滿足的常微分方程進(jìn)行形式變換得到了用于求解時(shí)域Green函數(shù)的線性時(shí)變系統(tǒng)，并進(jìn)一步利用維數(shù)擴(kuò)展和引進(jìn)新的參變量，得到了用于求解時(shí)域Green函數(shù)的線性時(shí)不變系統(tǒng)。其次，根據(jù)Zhong[10]提出的用于求解線性時(shí)不變系統(tǒng)的精細(xì)時(shí)程積分算法實(shí)現(xiàn)了Green函數(shù)的精確、高效求解。再次，基于九節(jié)點(diǎn)形函數(shù)提出了一種Green函數(shù)的制表插值策略。最后，利用該算法對(duì)一半球的強(qiáng)迫垂蕩運(yùn)動(dòng)的輻射波力進(jìn)行了計(jì)算，數(shù)值解與解析解符合良好。

1 基礎(chǔ)理論

1.1 時(shí)域自由面格林函數(shù)

假定流體無粘、不可壓縮，流動(dòng)無旋，并且忽略自由表面張力的影響。則流場(chǎng)中的物理量可以利用速度勢(shì)進(jìn)行表達(dá)。 取場(chǎng)點(diǎn)p（ x，y，z，t ）和源點(diǎn) q（ ξ， η ，ζ，t′），則滿足線性自由表面邊界條件的無限水深時(shí)域格林函數(shù)表達(dá)式可以寫為：

式中：δ為Dirac脈沖函數(shù)，H為單位階躍函數(shù)[11]，r和r′分別是場(chǎng)點(diǎn)到源點(diǎn)及場(chǎng)點(diǎn)到源點(diǎn)關(guān)于靜水面的映像點(diǎn)距離。（1）式中，G0為瞬時(shí)項(xiàng)，可用Hess-Smith方法得到；G?為記憶項(xiàng)，是自由面Green函數(shù)的計(jì)算難點(diǎn)。

通過引入新的參變量，可以將時(shí)域自由表面Green函數(shù)的記憶項(xiàng)寫為

利用同樣的方法，Green函數(shù)水平導(dǎo)數(shù)可寫為：

Green函數(shù)垂向?qū)?shù)可寫為：

為求解如上形式的Green函數(shù)及其空間導(dǎo)數(shù)值，Clement提出如下引理，對(duì)于雙參數(shù)函數(shù)Av，lμ，（）τ，0≤μ≤1，定義如下：

是如下微分方程的解

對(duì)比（6）式和（11）式可知，v=0，l=1/2，所以G?μ，（）τ滿足如下四階常微分方程：

利用同樣的方法，可以得到Green函數(shù)空間導(dǎo)數(shù)所滿足的四階微分方程。Clement利用四階Runge-Kutta方法對(duì)上式進(jìn)行了求解，但是長時(shí)間的模擬，計(jì)算精度會(huì)逐漸降低。

1.2 速度勢(shì)初邊值問題

對(duì)于零航速時(shí)域波浪和物體相互作用問題，速度勢(shì)滿足如下邊界條件：

式中：k=1，…，6表示輻射速度勢(shì)，nk為廣義法向量，g是重力加速度，D、SF和SB分別代表流體域、自由表面和物體表面。

對(duì)于如（14）式所示的初邊值問題可以用邊界元方法進(jìn)行求解，King[5]給出了速度勢(shì)φk滿足的邊界積分方程：

根據(jù)輻射速度勢(shì)所代表的物理意義，可以將輻射速度勢(shì)分解為瞬時(shí)項(xiàng)和記憶項(xiàng)，并利用脈沖響應(yīng)函數(shù)的概念對(duì)記憶速度勢(shì)進(jìn)行求解，具體可參見文獻(xiàn)[12]。

2 數(shù)值方法

2.1 線性時(shí)變系統(tǒng)求解

定義如下線性時(shí)變系統(tǒng)

引入如下函數(shù)

可得：

則（16）式可以寫為

為求解如（21）式所示的單位線性時(shí)變系統(tǒng)[13]，引入如下額外變量

以及n m+r+（）

1維的列向量

則如（21）式所示的單位時(shí)變系統(tǒng)可以轉(zhuǎn)化為

式中：

由（25）式可知，若s是一個(gè)小量，而且m取為一個(gè)較大的整數(shù)，則（24）式可以看作是一個(gè)線性時(shí)不變系統(tǒng)。

對(duì)于線性時(shí)不變系統(tǒng)，可以用Zhong[10]所提出的精細(xì)時(shí)程積分法進(jìn)行精細(xì)求解：

如上微分方程的解為

取時(shí)間步長Δt，則時(shí)間序列為

相應(yīng)時(shí)刻的解為Xk=X（ kΔt ） 。 對(duì)于區(qū)間 （tk， tk+1），（28）式可以寫為：

（31）式可以采用2N算法進(jìn)行精確計(jì)算，

式中：M=2N，N為任意正整數(shù)，例如

所以Δt′=Δt/M為一個(gè)極小的數(shù)，對(duì)T Δt（）′關(guān)于自變量進(jìn)行泰勒級(jí)數(shù)展開如下：

舍入誤差為（AΔt′）4/120。 此外，Ta（Δt′）與單位矩陣E，相比為一小量，故需要單獨(dú)進(jìn)行計(jì)算。 由（32）式可得

因此增量函數(shù)T（ Δt）=Ta（2NΔt′）最終可以根據(jù)下式求得：

2.2 時(shí)域自由面格林函數(shù)計(jì)算

對(duì)于時(shí)域自由表面Green函數(shù)所滿足的四階微分方程（13），可以寫為：

易知，（37）式為如（16）式所示的線性時(shí)變系統(tǒng)，可以采用2.1節(jié)所述的方法進(jìn)行求解。此外，對(duì)于長時(shí)間的模擬，可以將總的模擬時(shí)間區(qū)間［0，T ］分為若干個(gè)子區(qū)間H=[τn+1-τn），以利用較小的m和N獲得更高的計(jì)算精度。

3 算法有效性驗(yàn)證

3.1 精細(xì)積分算法有效性驗(yàn)證

本文首先計(jì)算了Green函數(shù)自身G（（ μ，τ）在 μ=0和 μ=1情況下的函數(shù)值，并與Clement（1998）提供的解析值進(jìn)行了對(duì)比。

對(duì)于 μ=0：

對(duì)于 μ=1：

由圖1和圖2可知，利用本文所提出的改進(jìn)的Green函數(shù)精細(xì)時(shí)程積分算法與解析解符合良好，且計(jì)算精度和穩(wěn)定性要明顯優(yōu)于四階Runge-Kutta方法。此外，圖3-5給出了Green函數(shù)及其導(dǎo)數(shù)在μ-τ平面內(nèi)的三維變化趨勢(shì)圖，為Green函數(shù)的制表插值策略提供了參考。

圖1 Green函數(shù)G 0，（）τ 隨時(shí)間變化趨勢(shì) Δτ=0.02，m=30，N=20，H=5Fig.1 Time history of（G 0，（）τ Δτ=0.02,m=30,N=20,H=5

圖2 Green函數(shù)（G 1，（）τ 隨時(shí)間變化趨勢(shì) Δτ=0.02，m=30，N=20，H=5Fig.1 Time history of（G 1，（）τ Δτ=0.02,m=30,N=20,H=5

圖3 Green函數(shù)G（（ μ，τ）Fig.3 3D plot ofG（（ μ，τ）

圖4 Green函數(shù)水平導(dǎo)數(shù)（GRμ，（）τFig.4 3D plot of（GRμ，（）τ

圖5 Green函數(shù)垂向?qū)?shù)（GZμ，（）τFig.5 3D plot of（GZμ，（）τ

3.2 插值算法有效性驗(yàn)證

為了避免對(duì)于不同形式的浮體在計(jì)算水動(dòng)力時(shí)都需要重新計(jì)算格林函數(shù)，本文采用了有限單元法中的形函數(shù)對(duì)Green函數(shù)波動(dòng)項(xiàng)在μ-τ平面內(nèi)進(jìn)行了插值計(jì)算。首先將預(yù)先計(jì)算好的Green函數(shù)值寫成二進(jìn)制的表格，在計(jì)算浮體水動(dòng)力時(shí)，將其預(yù)先讀入計(jì)算機(jī)的內(nèi)存中，并通過有限元中的插值形函數(shù)[14]對(duì)其進(jìn)行插值：

式中

為驗(yàn)證本文所提出的插值算法的精確性和有效性，對(duì)一半球做強(qiáng)迫垂蕩運(yùn)動(dòng)時(shí)的運(yùn)動(dòng)波浪力進(jìn)行了求解，并通過傅里葉變換[12]將其轉(zhuǎn)換為頻域水動(dòng)力系數(shù)，半球的主尺度及離散網(wǎng)格數(shù)目如圖6所示。

圖6 半球主尺度及網(wǎng)格離散數(shù)Fig.6 Grids and main diemensions of the hemisphere

圖7 無因次半球垂蕩附加質(zhì)量Fig.7 Nondimesional added mass of the hemisphere

圖8 無因次半球垂蕩阻尼系數(shù)Fig.8 Nondimesional damping coefficients of the hemisphere

由圖7和圖8可知，利用本文所提出的改進(jìn)精細(xì)時(shí)程積分算法和有限元形函數(shù)制表插值策略求得的半球水動(dòng)力系數(shù)與解析解[15]吻合良好，顯示了該插值算法的精確性和有效性。

4 結(jié) 論

精確高效的計(jì)算時(shí)域自由面格林函數(shù)對(duì)在時(shí)域內(nèi)求解船舶的運(yùn)動(dòng)響應(yīng)至關(guān)重要。

本文通過對(duì)時(shí)域自由面格林函數(shù)所滿足的四階微分方程進(jìn)行形式變換，得到了標(biāo)準(zhǔn)形式的線性時(shí)變系統(tǒng)，并通過對(duì)該線性時(shí)變系統(tǒng)進(jìn)行維數(shù)擴(kuò)展和引進(jìn)新的參變量，得到了可以精細(xì)求解的線性時(shí)不變系統(tǒng)，從而實(shí)現(xiàn)了格林函數(shù)的精確求解。

以該改進(jìn)的精細(xì)時(shí)程積分算法為基礎(chǔ)，提出了一種基于有限元形函數(shù)的制表插值策略，對(duì)零航速浮體的瞬時(shí)運(yùn)動(dòng)波浪力進(jìn)行了數(shù)值計(jì)算。通過與解析解的比較，驗(yàn)證了該方法的精確性、穩(wěn)定性和高效性。

參 考 文 獻(xiàn)：

[1]Chen J K,Duan W Y,Zhao B B,Ma Q W.Time domain hybrid TEBEM for 3D hydrodynamics of ship with large flare at forward speed[C]//The 32nd International Workshop on Water Waves and Floating Bodies,2017.Dalian,China,2017.

[2]Chen X,Zhu,R C,Zhou,W J,Zhao J.A 3D multi-domain high order boundary element method to evaluate time domain motions and added resistance of ship in waves[J].Ocean Eng.,2018,159：112-128.

[3]Wehausen J V,Laitone E V.Surface waves[M].Encyclopaedia of Physics,1960,IX：446-778.

[4]Beck R,Liapis S.Transient motions of floating bodies at zero forward speed[J].Journal of Ship Research,1987,31(3)：164-176.

[5]King B.Time-domain analysis of wave exciting forces on ships and bodies[D].Ph D Thesis,University of Michigan,1987.

[6]Lin W M,Yue D.Numerical solutions for large-amplitude ship motions in the time domain[C]//The 18th Symposium on Naval Hydrodynamics.Washington,United States,1991.

[7]Huang D.Approximation of time-domain free surface Green function and its spatial derivatives[J].Ship Building,1992,1(4)：16-25.

[8]Clement A H.An ordinary differential equation for the Green function of time-domain free-surface hydrodynamics[J].Journal of Engineering Mathematics,1998,33(2)：201-217.

[9]Chuang J,Qiu W,Peng H.On the evaluation of time-domain Green function[J].Ocean Eng.,2007,34(7)：962-969.

[10]Zhong W X.On precise integration method[J].Journal of Computational and Applied Mathematics,2004,163(1)：59-78.

[11]Abramowitz M,Stegun I A.Handbook of mathematical functions：with formulas,graphs,and mathematical tables[M].Courier Dover Publications,1972.

[12]Liapis S.Time-domain analysis of ship motion[D].Ph D Thesis,University of Michigan,1986.

[13]Liu X M,Zhou G,Zhu S,Wang Y H,Sun W R,Weng S L.A modified highly precise direct integration method for a class of linear time-varying systems[J].Science China-Physics Mechanics&Astronomy,2014,57(7)：1382-1389.

[14]Zienkiewicz O C,Taylor R L,Zienkiewicz O C,et al.The finite element method[M].McGraw-hill London,1977.

[15]Hulme A.The wave forces acting on a floating hemisphere undergoing forced periodic oscillations[J].Journal of Fluid Mechanics,1982,121：443-463.