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

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

0 引 言

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

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

目前，VLFS 水彈性響應(yīng)分析大多是在線性波浪理論的框架下進行，少量考慮了波浪的弱非線性效應(yīng)。當(dāng)海洋工程向深海發(fā)展時，海洋結(jié)構(gòu)物會不可避免遭遇畸形波（freak wave），這是一種波高極大、波峰異常尖瘦、能量相當(dāng)集中的不規(guī)則波[5]，對海洋平臺安全性是個重要威脅。在此海洋環(huán)境下，必須考慮強非線性波浪下結(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)，通過與實驗結(jié)果相比較，驗證了方法的有效性；程勇（2017）[8]基于高階BEM 方法建立了完全非線性的二維時域數(shù)值水池，用于模擬聚焦波浪與彈性浮板的水彈性響應(yīng)，研究表明，浮體結(jié)構(gòu)在極端波浪情況下的水彈性響應(yīng)具有極強的非線性特征。以上研究仍是在勢流理論框架下進行的，當(dāng)遇到波浪涌上浮體上部以及波浪沖擊浮體后波浪破碎等強非線性效應(yīng)時，勢流理論存在一定的局限性。相比于勢流理論，計算流體力學(xué)（computational fluid dynamics，簡稱CFD）在描述流場中速度場與壓力場的改變、捕捉自由液面非線性運動等方面更為有效和具有準(zhǔn)確性?；贑FD理論，研究者們開發(fā)了CFD-FEM耦合方法來研究時域下浮體的水彈性響應(yīng)。Seng（2012）[9]基于OpenFoam 開發(fā)了一種CFD-FEM 的耦合方法，其中將彈性體簡化為梁模型進行計算，結(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）流體運動方程

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

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

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

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

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

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

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

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

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

1.2 自由表面捕捉方法

在有網(wǎng)格的自由表面數(shù)值模擬方法中，常使用流體體積分?jǐn)?shù)方法（volume of fluid method，簡稱VOF 法）來捕捉自由液面。交界面的相分布和位置由相體積分?jǐn)?shù)αi來進行描述。網(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）進行計算：

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

1.3 畸形波浪模型

（1）組合聚焦模型

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

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

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

（2）仿物理造波

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

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

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

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

1.4 CFD-FEM 耦合方法

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

2 數(shù)值計算

2.1 數(shù)值模型驗證

為了驗證CFD-FEM 方法的有效性，將本文計算的浮體水彈性響應(yīng)計算結(jié)果與文獻[13]的試驗結(jié)果以及文獻[14]三維勢流方法計算結(jié)果進行比較，其中彈性浮板模型參數(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)進行了網(wǎng)格加密。計算時間步為0.005 s，時間離散為2 階精度，時間步內(nèi)的迭代次數(shù)為10 次；監(jiān)測彈性浮板迎浪端、中部以及背浪端位置處的垂向位移的變化情況。圖2 為彈性浮板在迎浪端垂向位移最大時，浮板中線處的無因次垂向位移，并與試驗數(shù)據(jù)[13]及勢流計算結(jié)果[14]進行比較。由此可見，本文計算結(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)測波面的升高。

在進行數(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，達到畸形波Hmax/H＞2.2 的標(biāo)準(zhǔn)，生成的波浪達到了畸形波的要求。在6 m 位置浪高儀檢測到的最大波高為0.116 m，與理論計算的最大波高0.110 m 差異在5.1%；在9 m 位置檢測到的最大波高為0.132 m，與理論值波高差異在4.7%。從四個不同位置處浪高儀所測值可知，粘性水池的數(shù)值解與理論解具有相當(dāng)好的一致性。

2.3 畸形波浪中彈性浮板與剛性浮板運動響應(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ū)域為外部流場計算區(qū)域，施加流固耦合交界面條件，以求解浮板表面流體壓力載荷的變化；第三個區(qū)域為結(jié)構(gòu)計算區(qū)域，施加結(jié)構(gòu)位移運動條件，來求解在流體壓力載荷下浮板運動和彈性變形。

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

從圖8可見，剛性浮板垂向位移的幅值相比彈性浮板明顯要小，這是因為彈性浮板除了剛體運動垂向位移外，還疊加了彈性垂彎的位移。

如圖9 所示，對比甲板上浪的情況可以發(fā)現(xiàn)，剛性浮板的上浪情況比彈性浮板更嚴(yán)重，這是因為彈性浮板能夠隨著畸形波的運動產(chǎn)生垂彎位移，相應(yīng)地減小了甲板高度與波高的差距，導(dǎo)致彈性浮板甲板上浪的水量少于剛性浮板甲板上浪的水量，所以剛性浮板的甲板上浪情況更嚴(yán)重。在15 s時，畸形波抵達聚焦位置10 m 處，彈性甲板上已沒有積水，這是因為波峰離開迎浪端后，波谷到達迎浪端，而浮板會隨著畸形波的運動產(chǎn)生垂彎，導(dǎo)致甲板上的積水順著迎浪端的垂彎流入波谷；剛性浮板由于沒有變形運動以減小甲板與波浪的高度差，且甲板上的積水不能及時排出，所以甲板上積水嚴(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）三個不同有義波高，以及迎浪與橫浪兩個浪向進行模擬。畸形波的其他參數(shù)與2.3 節(jié)的一致，對于迎浪工況，彈性浮板迎浪端距離造波邊界6 m，畸形波聚焦位置在彈性浮板背浪端；對于橫浪工況，彈性浮板前側(cè)點距離造波邊界9 m，畸形波聚焦位置在彈性浮板的后側(cè)點。迎浪工況計算區(qū)域與2.3節(jié)的相同，橫浪工況依此類推。浮板模型參數(shù)與2.1節(jié)的相同。

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

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

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

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

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

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

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 的畸形波進行計算，畸形波其他參數(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)在迎浪端位置處，背浪端的位移最大值比迎浪端的略小，這是因為浮板中部位置是浮板的內(nèi)部節(jié)點，當(dāng)浮板出現(xiàn)變形時，該節(jié)點處會受到與之相鄰的其他節(jié)點的影響，在一定程度上使其運動狀態(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)。參照文獻[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)隨波高的增加而增加，彈性變形運動在總體運動中的占比較大；而在橫浪時，浮板的水彈性運動在總體運動中占比較小，剛性運動的占比提高。

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

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