時域自由面格林函數改進精細積分算法研究

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

（1.江蘇科技大學 船舶與海洋工程學院，江蘇 鎮江 212003；2.哈爾濱工程大學 船舶工程學院，哈爾濱 150001）

0 引 言

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

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

Beck和Liapis[4]根據Green函數的振蕩特性，將整個計算區域進行了劃分，通過在不同的區域上應用級數展開和漸進展開，對Green函數進行了計算。在此基礎上，King[5]增加了一個額外的小區間，在此區間上利用Bessel函數的展開式進行了計算。基于Newman的工作，Lin和Yue[6]綜合利用級數展開、漸進展開和二維多項式近似對Green函數進行了計算。黃德波[7]通過采用雙參數制表插值技術，提高了Green函數的計算效率，但是計算精度卻有所損失。

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

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

1 基礎理論

1.1 時域自由面格林函數

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

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

通過引入新的參變量，可以將時域自由表面Green函數的記憶項寫為

利用同樣的方法，Green函數水平導數可寫為：

Green函數垂向導數可寫為：

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

是如下微分方程的解

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

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

1.2 速度勢初邊值問題

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

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

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

根據輻射速度勢所代表的物理意義，可以將輻射速度勢分解為瞬時項和記憶項，并利用脈沖響應函數的概念對記憶速度勢進行求解，具體可參見文獻[12]。

2 數值方法

2.1 線性時變系統求解

定義如下線性時變系統

引入如下函數

可得：

則（16）式可以寫為

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

以及n m+r+（）

1維的列向量

則如（21）式所示的單位時變系統可以轉化為

式中：

由（25）式可知，若s是一個小量，而且m取為一個較大的整數，則（24）式可以看作是一個線性時不變系統。

對于線性時不變系統，可以用Zhong[10]所提出的精細時程積分法進行精細求解：

如上微分方程的解為

取時間步長Δt，則時間序列為

相應時刻的解為Xk=X（ kΔt ） 。 對于區間 （tk， tk+1），（28）式可以寫為：

（31）式可以采用2N算法進行精確計算，

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

所以Δt′=Δt/M為一個極小的數，對T Δt（）′關于自變量進行泰勒級數展開如下：

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

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

2.2 時域自由面格林函數計算

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

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

3 算法有效性驗證

3.1 精細積分算法有效性驗證

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

對于 μ=0：

對于 μ=1：

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

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

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

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

圖4 Green函數水平導數（GRμ，（）τFig.4 3D plot of（GRμ，（）τ

圖5 Green函數垂向導數（GZμ，（）τFig.5 3D plot of（GZμ，（）τ

3.2 插值算法有效性驗證

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

式中

為驗證本文所提出的插值算法的精確性和有效性，對一半球做強迫垂蕩運動時的運動波浪力進行了求解，并通過傅里葉變換[12]將其轉換為頻域水動力系數，半球的主尺度及離散網格數目如圖6所示。

圖6 半球主尺度及網格離散數Fig.6 Grids and main diemensions of the hemisphere

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

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

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

4 結 論

精確高效的計算時域自由面格林函數對在時域內求解船舶的運動響應至關重要。

本文通過對時域自由面格林函數所滿足的四階微分方程進行形式變換，得到了標準形式的線性時變系統，并通過對該線性時變系統進行維數擴展和引進新的參變量，得到了可以精細求解的線性時不變系統，從而實現了格林函數的精確求解。

以該改進的精細時程積分算法為基礎，提出了一種基于有限元形函數的制表插值策略，對零航速浮體的瞬時運動波浪力進行了數值計算。通過與解析解的比較，驗證了該方法的精確性、穩定性和高效性。

參 考 文 獻：

[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.