劉 一，朱仁慶，程 勇，謝 彤，李潤澤

(江蘇科技大學 船舶與海洋工程學院，江蘇 鎮(zhèn)江 212001)

超大型浮體是一種綜合性的海洋浮式結構物，在維護國家海洋權益方面有極其重要的作用[1]。超大型浮體水平尺度與垂向尺度相差非常大，是一種極為扁平的柔性結構物，過大的尺度差異導致其彎曲剛度降低，其彈性變形相比于船舶與其他海洋結構物不能被忽略，因此必須采用水彈性方法分析超大型浮體在流體載荷作用下的變形以及浮體變形對流場的影響。過去對浮體水彈性的求解方法是基于Bishop等[2]提出的二維水彈性力學理論，通過將復雜的海洋結構物簡化為梁模型,以模型的干模態(tài)疊加來表達結構的運動與變形,建立流固耦合運動方程，計算海洋結構物的水彈性響應。由于二維水彈性理論在計算時忽略了船寬對流體運動的干擾，Wu[3]將三維適航性理論與三維結構動力學理論相結合,提出可適用于分析波浪中任意三維彈性體的三維水彈性理論。杜雙興[4]提出了零航速三維振蕩源格林函數(shù)快速計算方法,使得三維水彈性分析工作可以在PC機上進行，并建立了完善的三維航行船體線性水彈性力學頻域分析方法[5],從而可不受浮體細長比和航速的限制。Murai等[6]提出了將大型浮體處理成多個小結構組合的計算方法，并將數(shù)值計算與試驗數(shù)據(jù)對比，兩者之間有良好的一致性。王大云[7]采用三維時域格林函數(shù),推導出應用于彈性體的三維勢流時域積分方程，建立了三維水彈性時域分析理論。Watanabe等[8]討論了超大型浮體(VLFS)錨鏈系統(tǒng)、海域海底的變化以及結構數(shù)值模型的選取，并進一步探討了未來VLFS研究的方向。陳徐均[9]建立了浮體二階水彈性力學分析方法,并推導錨泊浮體三維線性及非線性頻域水彈性力學運動方程，開發(fā)出計算系泊浮體的計算程序。付世曉[10]將剛性浮體的錨泊系統(tǒng)動力分析理論與三維水彈性理論相結合，推導并建立了能夠考慮浮體彈性影響的錨泊系統(tǒng)動力響應分析方法，并編制了相應的計算程序。

目前主要采用勢流理論對各類浮式結構物的水彈性響應進行計算分析。由于其對波浪破碎、砰擊等現(xiàn)象難以數(shù)值模擬，且忽略了流體黏性的影響，因而對強非線性波浪運動及其引起的結構大幅運動和變形不能精確模擬，對結構物近壁面的壓力變化也不能準確捕捉。相比于勢流理論，計算流體動力學(CFD)在描述流場中速度場與壓力場的改變，捕捉自由液面處波浪的非線性現(xiàn)象等方面有更高的精度和準確性。方昭昭等[11]在計算波浪對航行船舶水動力的影響時發(fā)現(xiàn)，相比于勢流計算結果，計算流體動力學方法更能真實的反映模擬流場，并且CFD方法計算結果更加精準。

隨著計算流體動力學的發(fā)展，如何采用CFD模擬海洋結構物與流體的耦合作用得到越來越多學者的關注。Ley等[12]開發(fā)了一種基于雷諾時均N-S方程(RANS)方法的CFD和動態(tài)FEA的單向和雙向耦合系統(tǒng)，并證明了CFD方法和動態(tài)FEA耦合技術可以很好地預測由規(guī)則波和不規(guī)則波激勵引起的載荷。Seng[13]采用OpenFoam開發(fā)了一種耦合方法，使用梁模型進行水彈性響應的計算，結果證明CFD方法在估算結構的水彈性響應方面具有良好的準確性。Tomoki等[14]使用CFD-FEA技術分析了6600TEU集裝箱船的水彈性響應，并計算了船體中剖面彎矩的變化，與三維水彈性方法進行比較，有良好的一致性。但是流體與結構之間的耦合是采用單項耦合，由于僅考慮單向耦合計算，不能實現(xiàn)流體與結構之間的雙向數(shù)據(jù)傳遞。

采用CFD-FEM方法實現(xiàn)雙向流固耦合，計算非線性波浪與彈性浮板的水彈性響應問題。首先，基于CFD方法建立黏性數(shù)值水池，采用速度入口造波原理模擬五階Stokes波浪；使用FEM方法對彈性浮板進行離散，建立彈性浮板與外部流場的交界面。計算過程中，浮板表面的壓力載荷通過數(shù)據(jù)映射的方式傳遞給有限元計算模塊的浮板模型，模型在壓力載荷的作用下發(fā)生形變，并將變形后的模型數(shù)據(jù)傳遞給CFD計算模塊，對流固耦合表面進行更新，在新的交界面上再進行流體質點速度場與壓力場的計算，在下一時間步內重復此步驟。在計算過程中CFD和FEM兩模塊的數(shù)據(jù)進行實時交互，從而實現(xiàn)交界面上浮板的運動與彈性形變協(xié)調。此外，還探討了在不同參數(shù)下，如浮板厚度、入射波波幅以及浮板的三維效應對彈性浮板水彈性響應的影響。

1 數(shù)值模型

1.1 控制方程

1.1.1 流體控制方程

當不考慮流體的壓縮性，不計流體表面張力時，使用連續(xù)方程與N-S方程來描述海洋工程問題中的黏性流體運動：

·u=0

(1)

(2)

式中：u為速度矢量，ρ為流體密度，p為壓力，g表示重力矢量，μ為黏度系數(shù)，其為動力黏度系數(shù)μfluid與渦動黏度系數(shù)μt之和。

1.1.2 結構控制方程

假設結構為線彈性材料，在波浪等外載荷作用下相對于原平衡位置做剛體運動和變形，其結構運動方程通過有限元方法得到：

(3)

式中：m為結構質量矩陣，c為結構阻尼矩陣，k為結構剛度矩陣，x為節(jié)點位移列陣，F(xiàn)(t)為外界各種力合成的等效節(jié)點力列陣。

對于線彈性材料而言，其應力應變關系為線性關系，由胡克定律給出：

σ=Dε

(4)

式中：D為材料正切系數(shù)，σ為應力，ε為應變。

1.2 數(shù)值方法

1.2.1 湍流模型

在求解工程實際問題時，通常將連續(xù)方程與N-S方程各項取時均，此時N-S方程成為雷諾時均N-S(RANS)方程。在求解RANS方程時，由于引入渦動黏度系數(shù)μt，將導致方程不封閉，需引入湍流模型來進行計算。在海洋工程計算中，SST k-ω模型在方程中增加了交叉擴散項，并且在湍流黏性系數(shù)中考慮了剪切應力的影響，使其計算穩(wěn)定性好，計算效率與計算精度較高，故計算中使用的湍流模型為SST k-ω模型。SST k-ω模型的輸運方程為：

(5)

(6)

1.2.2 數(shù)值造波方法

非線性波浪相比于線性波浪而言，其波峰較陡，波谷較為平坦，呈現(xiàn)出的是一種非對稱曲線，更加符合現(xiàn)實中的實際波浪。選取五階波浪Stokes進行波浪的模擬；在數(shù)值水池造波邊界處，流體質點的速度和波面瞬時升高滿足式(7)～(10)中的條件。

五階Stokes波的波面方程為：

(7)

x方向的速度為：

(8)

z方向的速度為：

(9)

在式(7)～(9)中的各項系數(shù)為：

(10)

式中：ω′為圓頻率，k′為波數(shù)，d為水深；參數(shù)c的定義為c=cosh(kd)；式(10)中的參數(shù)在文獻[15]中對各項有詳細的計算過程與定義。

1.2.3 數(shù)值消波方法

對數(shù)值水池而言，常規(guī)的消波方法是在水池的尾部出口位置加設消波區(qū)域，以消減在出口邊界上的波浪反射對計算區(qū)域的影響。但是考慮到由于水池中浮式結構物的存在，在波浪沖擊到物體上時，會產生反射波浪，并在一定程度上影響到速度入口造波的準確性。為消除物體反射波浪以及水池出口處的反射波浪對速度造波入口邊界的影響，將在水池入口區(qū)域與出口區(qū)域都設置消波區(qū)域，如圖1所示。在圖1的數(shù)值水池入口消波區(qū)采用強迫波形消波方法(wave forcing)，通過在動量方程中增加源項來實現(xiàn)波形的強制。

在其入口消波區(qū)域的動量方程中增加的源項為：

qφ=-γρ(φ-φ*)

(11)

式中：γ為阻尼系數(shù)，φ為動量方程的數(shù)值求解結果，φ*為理論計算的數(shù)值求解結果。

在圖1的數(shù)值水池出口消波區(qū)采用阻尼消波的方法(wave damping)，通過在波浪運動的垂直方向上增加阻尼以增強波浪的耗散。在動量方程中附加的阻尼源項為：

(12)

(13)

式中：u3為在垂向的速度分量；f1,f2分別為線性阻尼項與非線性阻尼項；nd為沿波浪傳播方向上的阻尼系數(shù)，nd的取值參考文獻[16]；xsd,xed分別為阻尼區(qū)域開始位置坐標與結束位置坐標。

圖1 數(shù)值水池消波區(qū)域Fig. 1 Area of wave damping in numerical tank

1.3 耦合方法

在耦合計算中，采用CFD和FEM之間的雙向耦合方法，該過程在圖2中進行簡要說明，圖中的t0為初始時刻，Δt表示每個耦合內的時間增量。在初始時間使用CFD計算出浮板表面的壓力以數(shù)據(jù)映射的方法傳遞給有限元模塊中的浮板模型，在壓力載荷的作用下，浮板有限元節(jié)點產生速度與加速度的變化，將導致流固耦合交界面變形；之后把變形后的節(jié)點數(shù)據(jù)傳遞給CFD計算程序，進行交界面的更新。由CFD計算出的壓力場和速度場以及FEM計算出節(jié)點的速度和加速度將傳遞給下一時間步。耦合計算中數(shù)據(jù)映射至關重要，因為CFD與FEM之間的網(wǎng)格離散不同，網(wǎng)格節(jié)點不相對應，計算中的數(shù)據(jù)映射采用形狀函數(shù)插值的方法進行。

圖2 CFD-FEM雙向耦合流程Fig. 2 Flowchart of two-way coupling of CFD-FEM

2 數(shù)值計算

2.1 計算模型驗證

計算模型與文獻[17]中的一致，彈性浮板模型參數(shù)見表1。由于數(shù)值水池前端設置了消波區(qū)，浮體前端與數(shù)值水池入口處的間距可以縮短，這里選擇為L；浮體尾端到數(shù)值水池出口處的距離為1.5L；為減少水池兩側邊界對計算的影響，設置浮體到水池兩側邊界的距離為2B。入射波波高Hw=0.02 m，波長分別取λ=2 m與λ=4 m。流場采用正六面體網(wǎng)格劃分，彈性浮板采用四面體單元劃分。為了提高數(shù)值模擬精度，在自由液面上下波動空間和彈性浮板所在區(qū)域及附近流域內進行了網(wǎng)格加密。流場網(wǎng)格與結構單元劃分如圖3所示。在流體與結構交界面上，流場網(wǎng)格節(jié)點與單元節(jié)點不重合，在計算過程中，在交界面上進行數(shù)據(jù)映射，實現(xiàn)計算數(shù)據(jù)的同步交互。

表1 彈性浮板模型參數(shù)Tab. 1 Model parameters of floating plate

圖3 計算流域網(wǎng)格及彈性浮板單元劃分Fig. 3 Meshing of numerical tank and plate

圖4為彈性浮板在迎浪端垂向位移最大時，浮板中線處的無因次垂向位移，并與文獻[17]中試驗數(shù)據(jù)與文獻[18]中勢流計算結果進行比較；其中浮板中線位置是指浮板寬度上的中心位置，而無因次化的垂向位移為節(jié)點垂向位移與波幅的比值。由此可見，計算結果與兩者的結果吻合較好，說明CFD-FEM方法能有效模擬彈性浮體在非線性波浪下的運動，采用CFD-FEM方法模擬非線性波浪作用下超大型浮體水彈性響應具有可行性。

圖5為波長2 m時，彈性浮板中線處迎浪端垂向位移時歷曲線?？梢钥闯?，在8 s之后，浮板迎浪端的垂向位移達到較為穩(wěn)定的狀態(tài)，相鄰兩個垂向位移峰值處的時間間隔與入射波的波浪周期相同，為T=1.13 s；同時選取彈性浮板運動的穩(wěn)定時間段內3個相鄰迎浪端垂向位移最大處的時刻，對比在3個時刻下中線各點處的垂向位移差異(初始時刻為T0=9.17 s)，從圖6中可以發(fā)現(xiàn)，在垂向運動穩(wěn)定時間段內，在固定的波浪周期間隔之間，彈性浮板中線各處的垂向位移保持一致，3個相鄰時間間隔內的垂向位移數(shù)值差異很小。這是由于在規(guī)則波浪下，隨著計算時長的增加，浮板彈性變形運動趨于穩(wěn)定，相鄰時間間隔內的垂向位移差異較小，圖中3個時刻的彈性變形保持一致。

圖5 彈性浮板中線迎浪端垂向位移時間歷程Fig. 5 Variation of time series of vertical displacement at fore-end of the elastic plate

圖6 3個相鄰時刻下的浮體中線垂向位移Fig. 6 Vertical displacement of VLFS at three adjacent moments

2.2 波高對浮板水彈性響應的影響

研究入射波波高對彈性浮板水彈性作用的影響，選取在給定的入射波長λ=2 m與λ=4 m時，3種不同入射波高下(H=0.02 m，0.03 m，0.04 m)，彈性浮板中線處迎浪端、中部與尾端位置處的垂向位移時間歷程結果，分別如圖7和圖8所示。

圖7 λ=2 m時彈性浮板在不同波高條件下的垂向位移時間歷程Fig. 7 Vertical displacement of time series of the plate at fore-end, mid-position and back-end for different wave heights with λ=2 m

圖8 λ=4 m時彈性浮板在不同波高條件下的垂向位移時間歷程Fig. 8 Vertical displacement of time series of the plate at fore-end, mid-position and back-end for different wave heights with λ=4 m

從圖7與圖8可見，在兩種波長條件下，彈性浮板各點的垂向位移隨著波高的增加而增加。這是由于當入射波波高增大時，浮板各點處受到的波浪載荷增加，浮板的水彈性響應加劇，導致浮板各位置處的垂向位移數(shù)值增加。

在波長2 m時，入射波的周期為1.13 s，圖7中截取時間段8 s與14 s之間穩(wěn)定段的迎浪端、中部以及尾端3個位置處垂向位移的時間歷程。在此時間段內，浮板各點處的垂向位移變化趨于穩(wěn)定，3個不同波高下的垂向位移峰值雖有較大的差異，但是三者垂向位移變化周期是保持一致的。這是因為當入射波的波浪周期確定后，波幅的變化只會影響彈性浮板水彈性響應的劇烈程度，對其運動的周期無較大的影響。在圖8中，波長4 m時各點處的垂向位移變化周期要遠大于波長2 m時的變化周期，但是其不同波高下的各點垂向位移周期也趨于一致。

在迎浪端位置處，兩個波長下不同波浪幅值時的垂向位移峰值基本一致，這是因為迎浪端垂向位移的變化主要是受到入射波的影響，其運動狀態(tài)與波浪基本保持一致，浮板其他位置處的垂向位移變化對迎浪端的影響較小。

相比于入射波波長λ=2 m而言，較大的入射波長λ=4 m時的中部與尾端的垂向位移增加的更為明顯，這是由于波長4 m時，浮板迎浪端與尾端的距離為1倍的波長，當浮板的迎浪端處在波浪的波峰(波谷)位置處，浮板兩端同時存在較大的垂向位移，從而導致中部位置的垂向位移也進一步加?。欢诓ㄩL2 m時，浮板迎浪端與尾端的距離為1.5倍的波長，當迎浪端處在波浪波峰(波谷)位置處，尾端位于波浪波谷(波峰)位置處，浮板的迎浪端與尾端存在相反的垂向位移，從而導致浮板中部的水彈性響應減弱，并沒有出現(xiàn)較大的抬升。在同樣波長下，波高的增加會引起彈性浮板總體彈性變形運動的加劇，而波長與浮板長度的關系，在一定程度上會影響彈性浮板內部各點處的最大垂向位移變化。

2.3 板厚對水彈性響應的影響

研究彈性浮板的板厚對浮板水彈性作用的影響，選取在給定入射波波高、波長時，3種不同浮板厚度(d=0.06 m，0.10 m，0.14 m)下，迎浪端、中部和尾端3個位置處的垂向位移時間歷程，結果如圖9所示。當模型材料屬性(彈性模量與泊松比)未發(fā)生改變時，隨著彈性浮板厚度的增大，浮板的截面彎曲剛度也隨之增加。

從圖9(a)與圖9(c)可以看出，在浮板迎浪端與尾端位置處，垂向位移的變化趨于一致，3種不同板厚在浮板自由端位置處垂向位移雖有差異，但數(shù)值相差較小。這是由于在浮板自由端位置，入射波的波浪載荷對浮板彈性變形運動的影響要遠大于截面彎曲剛度的影響，從而導致在自由端處的垂向位移幅值與入射波的波幅相一致。在圖9(b)的浮板中部位置處，厚度d=0.14 m情況下，浮板中部節(jié)點的垂向位移峰值要小于厚度為d=0.06 m與d=0.10 m處的垂向位移峰值。這是因為在浮板中部位置，彎曲剛度對浮板彈性變形的影響要大于波浪載荷的影響，當浮板的板厚增加時，截面彎曲剛度也隨之增加，導致中部節(jié)點的垂向位移逐漸減少，整體的水彈性響應也進一步減弱。

綜上所述，浮板厚度的改變會影響彈性浮板的水彈性響應運動，對浮板中部位置的影響要遠大于對自由端位置的影響。當浮板厚度在一定范圍內變化時，自由端的垂向位移峰值不會出現(xiàn)較大的變化。

圖9 彈性浮板在不同板厚條件下的垂向位移時間歷程Fig. 9 Time series of the plate at fore-end, mid-position and back-end for different values of plate thickness

2.4 彈性浮板的三維效應

當彈性浮板在入射波浪的作用下出現(xiàn)彈性變形運動時，由于浮板中線位置處的波浪載荷與兩側自由端處的波浪載荷之間存在一定的差異，將導致浮板各位置處的垂向位移存在一定的差異。圖10為入射波高為H=0.04 m，波長λ=4 m時，浮板中間位置，以及兩側自由端處垂向位移。

在圖10中可以看出，彈性浮板中間位置處各節(jié)點的垂向高度要高于兩側自由端位置處的垂向高度，而對于兩側自由端而言，兩者的垂向高度趨于一致。

從圖11的垂向位移云圖中同樣可以看出，浮板各點垂向高度從中間位置向兩側自由端逐漸降低，兩側的位移云圖以浮板中心位置處呈對稱分布。

圖10 浮板不同位置處的垂向位移(迎浪端位移最大)Fig. 10 Vertical displacement at different positions of the plate

圖11 浮板不同位置處的垂向位移云圖(迎浪端位移最大)Fig. 11 Vertical displacement at different positions of the plate

3 結 語

利用CFD建立數(shù)值黏性水池，并采用速度入口造波方法完成五階Stokes波浪的模擬；通過對彈性浮板進行有限元離散，并在外部流場與結構之間的交界面上進行數(shù)據(jù)交互，實現(xiàn)雙向耦合的同步計算。通過與文獻中的數(shù)值計算結果和試驗結果進行比較，表明該數(shù)學模型可以較準確地模擬彈性浮板在非線性波浪作用下的彈性變形運動。

研究發(fā)現(xiàn)：在板厚不變，入射波幅增加時，浮板的水彈性響應運動也隨著加劇，浮板各點的垂向位移隨波幅的增加而增加；在給定波浪環(huán)境時，浮板自身厚度的變化對浮板的垂向位移有一定的影響，在浮板厚度變化時，迎浪端與尾端兩處的自由端垂向位移并沒有隨浮板厚度出現(xiàn)較大的改變，而浮板內部節(jié)點隨厚度的增加，其垂向位移減小；在浮板兩側自由端位置處的垂向高度相比于中間位置的垂向高度稍小，浮板的垂向高度從中間位置向兩側逐漸減小。

CFD-FEM方法在計算大型浮體在非線性波浪環(huán)境下的水彈性響應時具有較好的精度，今后可以采用此方法進一步分析不均勻海底平面與不規(guī)則波浪對大型浮體水彈性響應的影響。