朱仁慶，謝 彤，劉 一，苑中排

（江蘇科技大學(xué)船舶與海洋工程學(xué)院，江蘇 鎮(zhèn)江 212003）

0 引 言

大型浮式海洋結(jié)構(gòu)物（very large floating structure，簡稱VLFS）是開發(fā)海洋資源和利用海洋空間的重大海洋裝備，可用作海上機(jī)場、島礁浮式保障平臺、海上航天發(fā)射場、漂浮海上生態(tài)城市等[1-3]。VLFS 體現(xiàn)兩個特點(diǎn)：一是其某一維度的尺度要達(dá)到數(shù)公里；二是在波浪等外載荷作用下的彈性運(yùn)動特別明顯而必須采用水彈性理論對其載荷與響應(yīng)進(jìn)行求解和分析。

船舶與海洋工程水彈性理論自上世紀(jì)70 年代末提出到應(yīng)用，已經(jīng)從二維拓展到了三維、從線性拓展到了非線性、從勢流理論向CFD領(lǐng)域擴(kuò)展[4]。

目前，VLFS 水彈性響應(yīng)分析大多是在線性波浪理論的框架下進(jìn)行，少量考慮了波浪的弱非線性效應(yīng)。當(dāng)海洋工程向深海發(fā)展時，海洋結(jié)構(gòu)物會不可避免遭遇畸形波（freak wave），這是一種波高極大、波峰異常尖瘦、能量相當(dāng)集中的不規(guī)則波[5]，對海洋平臺安全性是個重要威脅。在此海洋環(huán)境下，必須考慮強(qiáng)非線性波浪下結(jié)構(gòu)的水彈性響應(yīng)。Takag（i1996）[6]在研究受風(fēng)暴或海嘯引起的非線性長波作用下彈性浮板水彈性響應(yīng)時發(fā)現(xiàn)，波浪的非線性對VLFS 水彈性響應(yīng)影響不可忽略；Liu（2002）[7]基于BEM-FEM 方法模擬了柔性浮體結(jié)構(gòu)在非線性波浪中的水彈性響應(yīng)，通過與實(shí)驗(yàn)結(jié)果相比較，驗(yàn)證了方法的有效性；程勇（2017）[8]基于高階BEM 方法建立了完全非線性的二維時域數(shù)值水池，用于模擬聚焦波浪與彈性浮板的水彈性響應(yīng)，研究表明，浮體結(jié)構(gòu)在極端波浪情況下的水彈性響應(yīng)具有極強(qiáng)的非線性特征。以上研究仍是在勢流理論框架下進(jìn)行的，當(dāng)遇到波浪涌上浮體上部以及波浪沖擊浮體后波浪破碎等強(qiáng)非線性效應(yīng)時，勢流理論存在一定的局限性。相比于勢流理論，計算流體力學(xué)（computational fluid dynamics，簡稱CFD）在描述流場中速度場與壓力場的改變、捕捉自由液面非線性運(yùn)動等方面更為有效和具有準(zhǔn)確性?；贑FD理論，研究者們開發(fā)了CFD-FEM耦合方法來研究時域下浮體的水彈性響應(yīng)。Seng（2012）[9]基于OpenFoam 開發(fā)了一種CFD-FEM 的耦合方法，其中將彈性體簡化為梁模型進(jìn)行計算，結(jié)果表明，該方法在計算彈性體水彈性響應(yīng)方面具有較高的精度；Lee（2014）[10]開發(fā)了一種基于RANS 方法的CFD 和動態(tài)FEM 的單向和雙向耦合算法，發(fā)現(xiàn)該方法可以很好地預(yù)報由規(guī)則波和不規(guī)則波激勵引起的載荷；Qin（2017）[11]基于自開發(fā)的CFD-FEM 完全耦合的求解器，對比了彈性板垂直下落到非線性畸形波與規(guī)則波上結(jié)構(gòu)響應(yīng)的不同。

本文將大型浮體簡化為帶有自由邊界的、材料為各向同性的彈性浮板，借助數(shù)值波浪水池，基于CFD-FEM方法，研究在畸形波浪作用下彈性浮板的水彈性響應(yīng)。

1 數(shù)值模型

1.1 控制方程

（1）流體運(yùn)動方程

當(dāng)忽略流體的壓縮性、不計流體表面張力時，描述流體運(yùn)動的連續(xù)性方程與N-S方程為

式中，xi(i= 1,2,3 )為直角坐標(biāo)系的三個分量；ui為xi方向速度；ρ為流體密度；t為時間；Fbi表示xi方向體積力；p為壓強(qiáng)；τij表示粘性應(yīng)力張量，其表達(dá)式為

其中，μ為粘度系數(shù)。

粘性流體運(yùn)動存在湍流運(yùn)動。目前處理湍流問題的方法主要有直接求解N-S 方程，又稱為DNS法；雷諾平均法，即RANS法；大渦模擬法，即LES法。這里采用RANS法，其方程為

現(xiàn)在雖有多種湍流模型可供選擇，但針對不同粘流問題需要選擇合適的湍流模型。在數(shù)值模擬考慮粘流影響的波流問題時，SSTk-ω模型是個常用的選擇。其在方程中增加了交叉擴(kuò)散項(xiàng)，并且在湍流粘性系數(shù)中考慮了剪切應(yīng)力的影響，使其計算穩(wěn)定性好，計算效率與計算精度較高，故本文采用此模型。SSTk-ω模型的輸運(yùn)方程為

式中，k與ω分別表示湍流動能與湍流耗散率，ui與uj為速度分量，G?k是指由于平均速度引起的湍動能k產(chǎn)生的項(xiàng)，Gω是單位耗散項(xiàng)，Γk和Γω分別是k和ω的有效擴(kuò)散項(xiàng)，Yk和Yω分別是k和ω的湍動能耗散項(xiàng)，Dω是交叉擴(kuò)散項(xiàng)，Sk與Sω是自定義源項(xiàng)。

（2）結(jié)構(gòu)運(yùn)動方程

假設(shè)結(jié)構(gòu)為線彈性材料，在波浪等外載荷作用下相對于原平衡位置作剛體運(yùn)動和變形，則線彈性結(jié)構(gòu)經(jīng)有限元離散后的總體運(yùn)動方程為

式中，M為結(jié)構(gòu)質(zhì)量矩陣，C為結(jié)構(gòu)阻尼矩陣，K為結(jié)構(gòu)剛度矩陣，[x]為節(jié)點(diǎn)位移矩陣，[F(t)]為外界各種力合成的等效節(jié)點(diǎn)力列陣。

1.2 自由表面捕捉方法

在有網(wǎng)格的自由表面數(shù)值模擬方法中，常使用流體體積分?jǐn)?shù)方法（volume of fluid method，簡稱VOF 法）來捕捉自由液面。交界面的相分布和位置由相體積分?jǐn)?shù)αi來進(jìn)行描述。網(wǎng)格單元中所有相的體積分?jǐn)?shù)總和為1，意味著該單元內(nèi)不存在空相，單元被所有相填充滿，不存在空隙。相i的體積分?jǐn)?shù)定義如下：

式中，Vi為網(wǎng)格單元中相i的體積，V為網(wǎng)格單元的體積，N為總相數(shù)。

根據(jù)體積分?jǐn)?shù)的值，可以區(qū)分網(wǎng)格單元中是否存在不同相或流體：

包含交界面的網(wǎng)格單元中同時存在兩種或者兩種以上的物理相，該網(wǎng)格單元可視為混合相，其體積與動力粘度可以通過公式（12）進(jìn)行計算：

式中，ρi為網(wǎng)格單元中相i的體積，μi是網(wǎng)格單元中相i的動力粘度。

1.3 畸形波浪模型

（1）組合聚焦模型

選取JONSWAP波浪譜來描述不規(guī)則海況，其表達(dá)式為

由于極限波聚焦模型中存在能量過于集中的問題，因此選取文獻(xiàn)[12]中的極限波-隨機(jī)波組合模型來實(shí)現(xiàn)畸形波浪的生成，表達(dá)式為

式中，通過指定水池位置xp，指定時刻tp疊加產(chǎn)生畸形波，N為疊加波浪的數(shù)量，坐標(biāo)系同2.2節(jié)中的坐標(biāo)系，當(dāng)Ep=0 時為極限波浪模型，Ep=1 時為隨機(jī)波浪模型。可通過改變Ep的數(shù)值來調(diào)整波譜能量在畸形波模擬中的分配問題，使模擬得到的波浪滿足畸形波的特征；Ai、k i和ωi分別為第i個組成波的振幅、波數(shù)和圓頻率；?i為隨機(jī)相位，是(0,2π)內(nèi)的隨機(jī)數(shù)。

（2）仿物理造波

采用推板造波的方法來實(shí)現(xiàn)畸形波的生成。根據(jù)推板造波的原理，推板造波運(yùn)動方程為

式中，Tr(ωi)為第i個組成波與造波板運(yùn)動速度之間的運(yùn)動傳遞函數(shù)，其表達(dá)式為

式中，h為數(shù)值水池的水深。

將畸形波的波面方程（14）代入公式（15）中，得到畸形波推板運(yùn)動方程為

1.4 CFD-FEM 耦合方法

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

2 數(shù)值計算

2.1 數(shù)值模型驗(yàn)證

為了驗(yàn)證CFD-FEM 方法的有效性，將本文計算的浮體水彈性響應(yīng)計算結(jié)果與文獻(xiàn)[13]的試驗(yàn)結(jié)果以及文獻(xiàn)[14]三維勢流方法計算結(jié)果進(jìn)行比較，其中彈性浮板模型參數(shù)為：板長L= 3.97 m，板寬B= 0.97 m，板厚D= 0.06 m，水深Dw= 2.3 m，吃水T= 0.037 m，彎曲剛度EI= 50.832 N·m，泊松比σ1= 0.3。彈性浮板的迎浪端距離速度入口造波邊界的距離為4 m，水池寬度為4 m，設(shè)置浮板位于水池的中心；入射波波高為0.02 m，波長分別取2 m 和4 m；流場采用正六面體網(wǎng)格劃分，彈性浮板采用四面體單元劃分，并在自由液面上下波動空間和彈性浮板所在區(qū)域及附近流域內(nèi)進(jìn)行了網(wǎng)格加密。計算時間步為0.005 s，時間離散為2 階精度，時間步內(nèi)的迭代次數(shù)為10 次；監(jiān)測彈性浮板迎浪端、中部以及背浪端位置處的垂向位移的變化情況。圖2 為彈性浮板在迎浪端垂向位移最大時，浮板中線處的無因次垂向位移，并與試驗(yàn)數(shù)據(jù)[13]及勢流計算結(jié)果[14]進(jìn)行比較。由此可見，本文計算結(jié)果與兩者的結(jié)果吻合較好，說明了采用CFD-FEM 方法模擬非線性波浪作用下超大型浮體水彈性響應(yīng)的可行性。

2.2 畸形波生成

選取JONSWAP 譜的波浪參數(shù)為：譜峰圓頻率ωp=3.491 rad/s（即譜峰周期TP= 1.8 s），疊加波浪的數(shù)量為N= 100，聚焦位置xp= 10 m，聚焦時間tp= 15 s。

數(shù)值波浪水池參數(shù)設(shè)置見圖4。

在距離造波邊界6 m、9 m、10 m與12 m處布置四個浪高儀來監(jiān)測波面的升高。

在進(jìn)行數(shù)值模擬時，計算初期數(shù)值計算的波面升高小于理論值計算結(jié)果，在圖5（a）中尤為明顯；但是在畸形波生成階段，不同浪高儀位置處檢測到的波浪升高與理論值貼合較好，在相位上無明顯差異。

由圖5（c）可觀察到，在tp= 15 s 時，在設(shè)定的聚焦位置xp= 10 m 處生成了畸形波，基于畸形波的判定標(biāo)準(zhǔn)，在10 m 處浪高儀檢測到的最大波高為0.128 m，與有義波高0.04 m 的比值為3.204，達(dá)到畸形波Hmax/H＞2.2 的標(biāo)準(zhǔn)，生成的波浪達(dá)到了畸形波的要求。在6 m 位置浪高儀檢測到的最大波高為0.116 m，與理論計算的最大波高0.110 m 差異在5.1%；在9 m 位置檢測到的最大波高為0.132 m，與理論值波高差異在4.7%。從四個不同位置處浪高儀所測值可知，粘性水池的數(shù)值解與理論解具有相當(dāng)好的一致性。

2.3 畸形波浪中彈性浮板與剛性浮板運(yùn)動響應(yīng)對比

選取有義波高為H1/3= 0.06 m、譜峰周期TP= 1.8 s的波譜生成畸形波，畸形波生成位置xp= 10 m，聚焦時刻tp= 15 s，疊加波浪的數(shù)量為N= 100，能量集中頻率區(qū)間為[ωL,ωH]=[2.318,6.891]。浮板模型參數(shù)與2.1節(jié)的相同，浮板迎浪端距離造波邊界6 m，畸形波聚焦位置在浮板背浪端。計算區(qū)域設(shè)置如圖6 所示，所有計算區(qū)域的設(shè)置以浮板的坐標(biāo)系為基準(zhǔn)，劃分為3 個區(qū)域：第一個區(qū)域是造波區(qū)域，生成畸形波；第二個區(qū)域?yàn)橥獠苛鲌鲇嬎銋^(qū)域，施加流固耦合交界面條件，以求解浮板表面流體壓力載荷的變化；第三個區(qū)域?yàn)榻Y(jié)構(gòu)計算區(qū)域，施加結(jié)構(gòu)位移運(yùn)動條件，來求解在流體壓力載荷下浮板運(yùn)動和彈性變形。

圖7為計算流域與彈性浮板結(jié)構(gòu)單元劃分圖，圖8給出了彈性與剛性浮板在迎浪端、中部與背浪端三個節(jié)點(diǎn)位置處浮板中線處垂向位移的變化。

從圖8可見，剛性浮板垂向位移的幅值相比彈性浮板明顯要小，這是因?yàn)閺椥愿“宄藙傮w運(yùn)動垂向位移外，還疊加了彈性垂彎的位移。

如圖9 所示，對比甲板上浪的情況可以發(fā)現(xiàn)，剛性浮板的上浪情況比彈性浮板更嚴(yán)重，這是因?yàn)閺椥愿“迥軌螂S著畸形波的運(yùn)動產(chǎn)生垂彎位移，相應(yīng)地減小了甲板高度與波高的差距，導(dǎo)致彈性浮板甲板上浪的水量少于剛性浮板甲板上浪的水量，所以剛性浮板的甲板上浪情況更嚴(yán)重。在15 s時，畸形波抵達(dá)聚焦位置10 m 處，彈性甲板上已沒有積水，這是因?yàn)椴ǚ咫x開迎浪端后，波谷到達(dá)迎浪端，而浮板會隨著畸形波的運(yùn)動產(chǎn)生垂彎，導(dǎo)致甲板上的積水順著迎浪端的垂彎流入波谷；剛性浮板由于沒有變形運(yùn)動以減小甲板與波浪的高度差，且甲板上的積水不能及時排出，所以甲板上積水嚴(yán)重。

2.4 畸形波波高對彈性浮板水彈性影響

為了研究不同有義波高下彈性浮板的水彈性響應(yīng)，這里選取H1/3= 0.05 m（Case-1），H1/3= 0.06 m（Case-2），H1/3= 0.07 m（Case-3）三個不同有義波高，以及迎浪與橫浪兩個浪向進(jìn)行模擬?；尾ǖ钠渌麉?shù)與2.3 節(jié)的一致，對于迎浪工況，彈性浮板迎浪端距離造波邊界6 m，畸形波聚焦位置在彈性浮板背浪端；對于橫浪工況，彈性浮板前側(cè)點(diǎn)距離造波邊界9 m，畸形波聚焦位置在彈性浮板的后側(cè)點(diǎn)。迎浪工況計算區(qū)域與2.3節(jié)的相同，橫浪工況依此類推。浮板模型參數(shù)與2.1節(jié)的相同。

2.4.1 迎浪畸形波對彈性板的作用

在迎浪端、中部與背浪端三個節(jié)點(diǎn)位置處監(jiān)測彈性浮板中線處垂向位移的變化情況。從圖10可以看出，當(dāng)聚焦位置和聚焦時間不變，聚焦波的有義波高增加時，彈性浮板垂向位移的變化幅值增大，彈性浮板的水彈性響應(yīng)運(yùn)動加劇,這是由于當(dāng)聚焦波波高增大時，浮板各點(diǎn)處受到的波浪載荷增加，浮板的水彈性響應(yīng)加劇，隨之影響浮板各位置處的垂向位移數(shù)值增加。圖中三種工況下垂向位移的幅值雖有差異，但是在整體的變形形態(tài)上是保持一致的，這是因?yàn)槿肷洳ɡ说慕M成波的幅值雖然不同，但是組成波的相位與頻率是相同的，使畸形波在波浪形態(tài)上存在一定程度的相似，波浪載荷作用于浮板上使之出現(xiàn)相似的變形形態(tài)。

2.4.2 橫浪畸形波對彈性板的作用

進(jìn)一步研究在橫浪時彈性浮板的垂向位移變化情況。計算中選取浮板前側(cè)點(diǎn)（P1）、中部點(diǎn)（P2）與后側(cè)點(diǎn)（P3）監(jiān)測垂向位移的變化情況，具體結(jié)果如圖11所示。

改變?nèi)肷洳ǖ挠辛x波高，計算結(jié)果如圖12 所示。在聚焦時刻15 s 附近時，聚焦位置處P3 監(jiān)測點(diǎn)的垂向位移的峰值分別為0.086 m、0.099 m 與0.107 m，垂向位移隨著入射波波高的增加而增加。

相比于迎浪，橫浪下浮板的剛體運(yùn)動占總體運(yùn)動的較大成分。從圖13 中Hw= 0.06 m 的波高下的側(cè)視圖中可以看出，浮板側(cè)面的各點(diǎn)基本保持在一條直線上，彈性變形較小。這是由于在橫浪時彈性浮板寬度方向尺度較小，浮板側(cè)面的彎曲剛度相比于正面的彎曲剛度要大，從而導(dǎo)致其彈性變形量減小，剛體運(yùn)動成為浮板運(yùn)動的主要成分。

2.4.3 上浪情況研究

浪向角為0°與90°時都出現(xiàn)波浪涌上板面的情況（見圖14），入射波高越大，上浪的水量越多。雖然浮板在遭遇波浪時存在一定的抬升，但是波浪的高度要大于浮板的抬升高度，從而出現(xiàn)波浪涌上浮板上部的現(xiàn)象。在浮板迎浪端抬升時，由于浮板上部存在一定的水體，其抬升高度將受到一定的影響，其峰值高度將略小于無上浪時的峰值變化。

2.5 畸形波作用位置對彈性浮板水彈性影響

為了研究畸形波浪位置對浮板水彈性作用的影響，選取了3種不同計算工況：浮板的迎浪端位于聚焦位置（Condition-1），浮板的中部位于聚焦位置（Condition-2）以及浮板背浪端位于聚焦位置（Condition-3）。圖15 顯示了三種不同工況下彈性浮板迎浪端與造波邊界的距離。采用有義波高H1/3=0.04 m 的畸形波進(jìn)行計算，畸形波其他參數(shù)與2.3節(jié)的一致。浮板模型參數(shù)與2.1 節(jié)的相同。從3 種工況可以看出，彈性浮板的水彈性響應(yīng)峰值出現(xiàn)時刻依次逐漸提前。這是由于彈性浮板離造波邊界的距離從Condition-1 到Condition-3 逐漸減小，遭遇到畸形波浪的時間逐漸提前。圖16 為Condition-1 和Condition-3中垂向位移各峰值的對比，可以看出，畸形波聚焦位置處都出現(xiàn)了垂向位移的峰值，同時注意到在Condition-2 中，聚焦位置在浮板中部，但是浮板水彈性響應(yīng)位移最大值卻出現(xiàn)在迎浪端位置處，背浪端的位移最大值比迎浪端的略小，這是因?yàn)楦“逯胁课恢檬歉“宓膬?nèi)部節(jié)點(diǎn)，當(dāng)浮板出現(xiàn)變形時，該節(jié)點(diǎn)處會受到與之相鄰的其他節(jié)點(diǎn)的影響，在一定程度上使其運(yùn)動狀態(tài)減弱。

2.6 不同譜寬畸形波對彈性浮板水彈性影響

為研究不同頻帶寬度下畸形波與彈性浮板的相互作用，選取中心頻率為fc=0.56 Hz，即譜峰周期為TP=1.8 s，模擬頻譜寬度為(a)Δf=0.3 Hz，(b)Δf=0.6 Hz，(c)Δf=0.9 Hz三種工況下浮板的水彈性響應(yīng)。參照文獻(xiàn)[12]選取頻率區(qū)間，即以中心頻率fc=0.56 Hz為基準(zhǔn)加上或減去0.5Δf來獲取頻率區(qū)間的上下限，有義波高Hs=0.05 m，畸形波其他參數(shù)與2.3節(jié)的一致。浮板模型參數(shù)與2.1節(jié)的相同。

圖17中為有義波高Hs= 0.05 m 時不同頻帶寬度條件下彈性浮板迎浪端、中部和背浪端的垂向位移的時歷曲線。在圖中可見，當(dāng)選取的頻帶寬度較窄時，浮板的垂向位移大于頻帶較寬的情況，這是由于頻帶寬度是能量集中的外在表現(xiàn)。當(dāng)頻帶越窄時，單位頻帶所含有的能量越高，從而導(dǎo)致波浪波幅、波浪載荷增加，使浮體的水彈性響應(yīng)加劇，在同一時刻同一位置處的位移增大。

3 結(jié) 論

本文基于三維畸形波浪水池以及CFD-FEM 方法，研究了大型浮體在畸形波作用下水彈性響應(yīng)問題，并獲得了以下結(jié)論：

（1）在不同有義波高的入射波浪下，迎浪時的水彈性響應(yīng)隨波高的增加而增加，彈性變形運(yùn)動在總體運(yùn)動中的占比較大；而在橫浪時，浮板的水彈性運(yùn)動在總體運(yùn)動中占比較小，剛性運(yùn)動的占比提高。

（2）當(dāng)波浪聚焦在浮板自由端位置時，自由端垂向位移明顯加大；而聚焦位置在浮板中部位置時，中部節(jié)點(diǎn)的位移幅值增加并不明顯。

（3）畸形波頻帶寬度越窄，波浪能量越集中，彈性浮板的垂向位移越大。而且在較小的頻帶寬度下，浮板表現(xiàn)出更強(qiáng)的非線性效應(yīng)，節(jié)點(diǎn)處垂向位移的改變更為明顯。