二維浮式結(jié)構(gòu)物與完全非線性波相互作用數(shù)值分析

顧思琪，朱仁慶，吳梓鑫

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

隨著海洋工程向深海發(fā)展，海洋結(jié)構(gòu)物所處的海洋環(huán)境更趨惡劣.在惡劣海浪作用下，可能會產(chǎn)生嚴(yán)重的波浪爬升、相對運動和上浪等流體強(qiáng)非線性沖擊現(xiàn)象，這些都會對海洋浮式結(jié)構(gòu)物產(chǎn)生很大的局部沖擊破壞載荷.因此，探索強(qiáng)非線性波浪與海洋浮式結(jié)構(gòu)物的相互作用、發(fā)展穩(wěn)定高效的數(shù)值預(yù)報方法，不僅對海洋水動力學(xué)的發(fā)展具有重要的理論意義，而且對海洋工程結(jié)構(gòu)物的設(shè)計和防護(hù)具有重要的工程價值.

早期，人們對波浪與結(jié)構(gòu)物完全非線性相互作用的研究大都以試驗為主.近十年來，經(jīng)過水動力學(xué)研究者們的不懈努力，數(shù)值方法不斷獲得突破，取得了令人鼓舞的成果［1］.目前，主要基于兩種理論模型對其展開研究，即勢流模型和粘性流體模型.基于勢流理論的數(shù)值模擬中，最經(jīng)典的是1976年由 Longuet-Higgins和 Cokelet［2］提出的混合歐拉-拉格朗日方法(MEL)數(shù)值模型，它可用于處理完全非線性波面運動.文獻(xiàn)［3］中基于勢流理論和高階邊界元模型，建立了無限水域中波物作用的完全非線性模型.文獻(xiàn)［4-6］中基于完全非線性的勢流理論(fully nonlinear potential theory，F(xiàn)NPT)，采用QALE-FEM方法模擬了方形數(shù)值水池中水波與坐底波形水壩的流固耦合響應(yīng)，并將其推廣成功實現(xiàn)了二維陡波與浮體的完全非線性流固耦合模擬.而基于粘流理論的數(shù)值水池，主要結(jié)合自由表面追蹤技術(shù)，求解連續(xù)方程、N-S方程.通過引入 VOF，MAC，SPH 和 LEVEL -SET 等［7-10］重構(gòu)自由面的數(shù)值方法，不但可以充分考查由于粘性作用導(dǎo)致的各種復(fù)雜流動特征，而且可以處理結(jié)構(gòu)物在波浪激勵下的大幅運動，以及波浪翻卷、破碎、合并等流體完全非線性運動等問題.因此，基于 N-S方程構(gòu)建數(shù)值波浪水池，研究粘流中波浪與結(jié)構(gòu)物的相互作用等問題，已成為當(dāng)今國際上的一個重點發(fā)展方向.文獻(xiàn)［11］中利用MAC差分格式，對歐拉動力學(xué)方程進(jìn)行離散求解，利用VOF法處理自由表面的演變，研究了二維完全非線性波在艙內(nèi)的運動特性.文獻(xiàn)［12］中用VOF法對FPSO的上浪問題作了研究，其計算對固定的FPSO模型可以得到較好結(jié)果，但對運動的船和波浪之間的相互作用問題則尚不能給出理想結(jié)果.文獻(xiàn)［13-14］中采用區(qū)域分解思想、遠(yuǎn)處流場和船體的運動通過勢流理論來計算，船體附近的流動通過求解N-S方程來模擬，這種分開求解的方式具有一定的局限性，與真實情況中波浪與船體的相互作用有較大差別.

文中以ANSYS-14.5軟件的Workbench為計算平臺，結(jié)合二次開發(fā)功能，建立了具有造波和消波功能的二維數(shù)值波浪水池.文中就干舷高度不同的結(jié)構(gòu)物遭遇波浪的完全非線性現(xiàn)象以及波浪對結(jié)構(gòu)物運動響應(yīng)和上浪的影響分別進(jìn)行了數(shù)值模擬研究，并采用VOF法追蹤自由液面，精確地描述了波浪與結(jié)構(gòu)物相互作用問題，得到了結(jié)構(gòu)物遭遇波浪后甲板上的上浪水位高度和結(jié)構(gòu)物的縱搖、垂蕩、橫搖等運動響應(yīng).

1 物理模型

1.1 流體運動控制方程

控制方程包括連續(xù)方程和不可壓縮粘性流體的N-S方程.

連續(xù)方程

對于不可壓縮流體，方程簡化為

式中:v為流體速度，ρ為流體密度.

動量方程

式中:ν為流體的運動粘性系數(shù)，F(xiàn)b為作用于流體的質(zhì)量力，p為流場壓強(qiáng).

1.2 自由液面追蹤

文中采用VOF法對自由液面進(jìn)行求解，由此確定自由液面的位置.自由液面的標(biāo)記采用流體體積函數(shù)法，體積分?jǐn)?shù)αq定義為單元內(nèi)第q相流體所占體積與該單元的體積比.若αq=1則表示該單元內(nèi)全部為第q相流體，若αq=0則表示該單元內(nèi)沒有第q相流體，若0＜αq＜1則表示該單元為交界面單元.αq滿足方程

1.3 波浪理論

對于stokes三階波，其波面方程和水質(zhì)點的速度分別為［15］:

波面方程

x方向速度分量

z方向速度分量

2 數(shù)值模型

2.1計算模型

文中所采用的二維數(shù)值波浪水池模型見圖1.波浪水池長15m，高2 m，其中水深1.5 m，水池末端5m長的區(qū)域為消波區(qū).坐標(biāo)系原點位于入射邊界的靜水面處，ox軸上方全部定義為空氣，密度1.225 kg/m3，ox軸下方全部定義為水，密度998.2kg/m3.

圖1 數(shù)值水池示意Fig.1 Schematic diagram of the numerical wave tank

圖1中將浮式結(jié)構(gòu)物簡化成一方形結(jié)構(gòu)，長1 m，吃水0.18 m，干舷0.04 m，船體中橫剖面位于x=5.5 m處.在數(shù)值計算中設(shè)置水位監(jiān)測點P1、P2、P3，分別位于x=5.0 m，5.1 m和5.2 m處.

流場網(wǎng)格和結(jié)構(gòu)網(wǎng)格全部采用四邊形結(jié)構(gòu)化網(wǎng)格劃分.網(wǎng)格的質(zhì)量直接影響計算結(jié)果的精度，文中經(jīng)過多次試算，確定了較合理的網(wǎng)格劃分方案:沿水池長度方向網(wǎng)格尺寸為dx=λ/100，沿水池高度方向網(wǎng)格尺寸為dz=H/8，在自由液面處對網(wǎng)格進(jìn)行加密以提高計算精度，其高度為dz=H/30.為了適當(dāng)?shù)販p少網(wǎng)格以減少計算工作量，在遠(yuǎn)離自由液面處網(wǎng)格以一定比例向上下逐漸稀疏.

2.2 邊界條件及參數(shù)設(shè)置

水池頂部邊界設(shè)置為壓力入口條件(Pressure-inlet)，底部邊界為無滑移固壁條件(Wall)，右側(cè)邊界設(shè)置為壓力出口條件(Pressure-outlet)，左側(cè)邊界定義為速度入口條件(Velocity-inlet).文中基于邊界造波法模擬stokes三階波，由于軟件Fluent的模塊中沒有提供完整的波浪模擬環(huán)境，因此需要通過用戶自定義模塊(UDF)，將UDF文件嵌入共享庫中并與Fluent連接.

同時，在Fluent中采用非定常分離隱式求解器求解，壓力方程選用加權(quán)體積力格式(Body Force Weighted)，壓力速度耦合方式采用PISO算法，體積分?jǐn)?shù)法為幾何重構(gòu)(Geo-Reconstruct)，壓力參考值為一個標(biāo)準(zhǔn)大氣壓，重力加速度為9.8066 m/s2，時間步長為0.005s.

2.3 消波設(shè)置

為防止水池末端波浪反射而干擾入射波場，在水池末端應(yīng)設(shè)置消波裝置.文中采用阻尼消波方法，即在特定的消波區(qū)域內(nèi)添加阻尼項以達(dá)到減弱或消除該區(qū)域內(nèi)的波動.具體實現(xiàn)方式為在水池后端設(shè)置阻尼消波區(qū)，通過在其動量方程中添加自編UDF程序 DEFINE_SOURCE(momentum，c，t，dS，eqn)實現(xiàn)消波.文中采用設(shè)置速度入口造波，為了保證流動的連續(xù)性，維持質(zhì)量守恒，在消波區(qū)內(nèi)只需對動量方程中 z方向的分量加載 DEFINE_SOURCE((z_momentum，c，t，dS，eqn)，x方向不作處理.

消波區(qū)內(nèi)，動量方程如式(9).

式中:ν為運動粘性系數(shù);C(x)為消波系數(shù)，其函數(shù)表達(dá)式為:

式中:L為波長;x0，xL分別為消波區(qū)左、右邊界的x坐標(biāo)值;α為阻尼經(jīng)驗參數(shù);ρ為流體密度.

3 結(jié)果與分析

3.1 波浪數(shù)值水池

穩(wěn)定的波浪是波物相互作用數(shù)值模擬的前提，為了驗證波浪數(shù)值水池的正確性，文中首先對波浪周期T=1.13s，理論波長λ=2 m，理論波高H=0.16 m的三階stokes波進(jìn)行了數(shù)值模擬.圖2給出了數(shù)值波浪水池在x=1 m和2 m處波面(A)隨時間的變化歷程曲線.在波浪的傳播過程中，隨著傳播距離的增加，波高是有衰減的，這是由于流體粘性作用的影響導(dǎo)致波浪能量耗散的緣故，因此，圖2中x=2 m處波幅計算值與理論值偏差較x=1 m處大.總體而言，數(shù)值模擬的波形比較穩(wěn)定，波浪周期、波幅參數(shù)等與理論解吻合良好，波形具有“上尖下平”的非線性特征.

圖3為不同時刻波浪水池消波區(qū)內(nèi)的波形圖，圖中較為清晰地顯示了波浪的傳播在受到消波區(qū)內(nèi)附加阻尼的影響后波形逐漸衰減的過程.從圖中可觀察到，在接近水池出口處，消波區(qū)自由液面高度值已趨于0 m，從而驗證了文中阻尼消波方法的可行性.

圖2 數(shù)值波形與理論解對比Fig.2 Comparison of numerical and theoretical wave forms

圖3 不同時刻消波區(qū)內(nèi)的波形Fig.3 Waves’shape in the damping domain at different time

3.2 波浪與浮體的相互作用

3.2.1 二維甲板上浪驗證實例

為驗證文中計算方法的可行性和有效性，下面將基于ANSYS Workbench計算平臺，對迎浪狀態(tài)下二維固定結(jié)構(gòu)物甲板上浪過程進(jìn)行數(shù)值模擬，并將計算結(jié)果與有關(guān)文獻(xiàn)中的試驗值進(jìn)行比較，以驗證數(shù)值方法的有效性.波物相互作用問題是雙向流固耦合問題，文中基于軟件ANSYS，采用迭代耦合的方法求解流固耦合問題，其主要思路是流體方程和結(jié)構(gòu)方程按順序相互迭代求解，各自在每一步得到的結(jié)果提供給另一部分使用，直到耦合系統(tǒng)的解達(dá)到收斂，迭代停止［16］.

1)計算模型說明

數(shù)值波浪水池和結(jié)構(gòu)模型如圖4，模型尺寸與文獻(xiàn)［16］中甲板上浪試驗?zāi)Ｐ捅３忠恢?圖中水池長13.5 m，其中2 m為消波區(qū)，水深1.035 m.結(jié)構(gòu)模型已簡化為方形結(jié)構(gòu)物，底部為半徑0.08 m的圓弧，長和型深分別為1.5 m和0.248 m，吃水為0.198 m，干舷為0.05 m.

圖4 計算域示意圖Fig.4 Diagram of computational domain

2)計算結(jié)果

在甲板上浪的數(shù)值模擬中，波浪要素與文獻(xiàn)［17］中試驗條件保持一致，波長和波高分別為2 m和0.16 m.為了與試驗結(jié)果對比，在模擬甲板上浪過程中對結(jié)構(gòu)上WL1和WL2兩點的波形時歷進(jìn)行了監(jiān)測，其中點WL1位于甲板迎浪側(cè)的首端點，點WL2與WL1的距離為0.075 m.圖5為文中數(shù)值模擬結(jié)果與試驗結(jié)果的對比圖，由圖可知，數(shù)值模擬結(jié)果與試驗結(jié)果基本吻合，具體表現(xiàn)在結(jié)構(gòu)甲板上浪時間以及甲板監(jiān)測點上浪水位高度(Hs)均與試驗結(jié)果基本吻合，從而驗證了文中計算方法的可靠性.如圖5所示，第一個波浪約在7s時到達(dá)WL1處，第二個波浪約在8s到達(dá)，且在第一個波浪涌上甲板時，部分水體被船首阻擋反射，使得船首附近的自由液面被抬高.因此，在第二個波到達(dá)船首時，波高增大許多，超過第一個波高的50%，引起的上浪現(xiàn)象也更為嚴(yán)重.

圖5 FPSO模型監(jiān)測點處波形時歷Fig.5 Time-history of wave at probe WL1 and WL2

從數(shù)值模擬結(jié)果與試驗值的對比可以看出，文中的計算方法是有效可行的，計算結(jié)果誤差較小，滿足計算要求，所以該方法適用于下文數(shù)值計算.

3.2.2 不同波浪要素下浮體與波浪的相互作用

具有波高、波長等波浪要素不同的規(guī)則波，對浮體的作用力也不盡相同.在與波浪接觸時，結(jié)構(gòu)物不僅會受到波浪的沖擊載荷，而且會受到上浪水體的影響.除此之外，由于波浪的形態(tài)變化，結(jié)構(gòu)物水下體積也隨之發(fā)生改變，浮心位置也會發(fā)生變化，因此結(jié)構(gòu)物在波浪中的運動響應(yīng)是由多種因素共同作用產(chǎn)生的.文中針對二維問題，只考慮浮體的橫搖、垂蕩和橫蕩運動，研究浮體中橫剖面模型在波浪作用下的載荷與運動響應(yīng).

圖6為干舷高度為0.04 m時，浮式結(jié)構(gòu)物在波高H分別為0.12，0.16 m的三階stokes波作用下的運動響應(yīng)時歷曲線和上浪水位高度值.從圖中可看出，在不同的波高下，結(jié)構(gòu)物的垂蕩與橫搖曲線呈現(xiàn)一定的周期性，且其周期與波浪周期大致相同.第一個波浪到達(dá)甲板并發(fā)生上浪現(xiàn)象的時間大致為6s，在經(jīng)歷了第一個上浪過程后，水體的質(zhì)量改變了結(jié)構(gòu)物的重心，結(jié)構(gòu)物的垂蕩時歷曲線有整體下移趨勢.就上浪水位而言，對于干舷相同的結(jié)構(gòu)物，波高大的波浪更容易造成上浪現(xiàn)象.

圖6 兩種波高的波浪作用下浮體的運動響應(yīng)與上浪水位Fig.6 Comparison between the data of floating body with different wave heights

圖7是干舷高度為0.04m時，λ分別為2m和2.5 m的三階stokes波作用下的運動響應(yīng)時歷曲線和上浪水位高度(圖中變量含義見上文).由于不同波長的影響，浮體的橫蕩值也會不同，波長較大的波浪作用下浮體的橫蕩值也隨之增大.波長分別為2.0 m和2.5 m的波浪，其周期分別為1.13s和1.27s，在這兩種波長作用下的結(jié)構(gòu)物垂蕩和橫搖時歷曲線的周期也不盡相同.觀察上浪水位高度圖可以發(fā)現(xiàn)，波長小的波浪上浪水位較大，即上浪量也隨著波陡的增加而增加.這是由于在波長較大的波浪作用下結(jié)構(gòu)物的橫搖幅值也較大，如圖7c)所示，因此在隨后的上浪過程中，上浪的水體較容易隨著結(jié)構(gòu)物的擺動流入水池中，產(chǎn)生上浪水位高度較小的現(xiàn)象.

圖7 兩種波長的波浪作用下浮體的運動響應(yīng)與上浪水位Fig.7 Comparison between the data of floating body with different wavelengths

3.2.3 不同干舷高度的浮體與波浪的相互作用

足夠的干舷不僅可以保證船舶有一定的儲備浮力，也可以減少船舶上浪.干舷高度的不同，受到波浪沖擊時結(jié)構(gòu)物的受力也不相同.圖8為在相同波浪下，干舷分別為0.02 m和0.04 m的浮式結(jié)構(gòu)物運動時歷曲線和上浪水位高度圖(圖中變量含義見上文).與圖6，7對比可發(fā)現(xiàn)，相對于波浪要素的影響來說，浮體干舷對結(jié)構(gòu)物的橫蕩值影響更大些，干舷較大的結(jié)構(gòu)物的漂移距離要遠(yuǎn)遠(yuǎn)大于干舷較小的結(jié)構(gòu)物.而對結(jié)構(gòu)物垂蕩和橫搖運動來說，干舷較大的結(jié)構(gòu)物其運動幅值也會較大，但差值不會特別明顯.在結(jié)構(gòu)物發(fā)生第一個上浪過程之前，干舷高度不同的兩個結(jié)構(gòu)物的垂蕩和橫搖曲線幾乎重合，在結(jié)構(gòu)物發(fā)生上浪后，兩者才產(chǎn)生了變化.原因在于，在相同的波浪下，干舷小的結(jié)構(gòu)物更容易上浪，沖上甲板的水體不會立刻流回水域，而會沿著甲板流動并持續(xù)一段時間，此時這部分水體也會影響結(jié)構(gòu)物的運動響應(yīng).

圖8 兩種干舷高度浮體的運動響應(yīng)與上浪水位Fig.8 Comparison between the data of floating body with different freeboards

圖9為波浪與浮體相互作用時的流體形態(tài)和結(jié)構(gòu)應(yīng)力云圖.從圖中可以觀察到，上浪現(xiàn)象是波浪和船首下沉運動兩種機(jī)制組合作用的結(jié)果.其過程是:水體沿船首升高，波浪高于船首涌到甲板上，水體沿甲板流動，從甲板另一側(cè)流出，主要對背浪側(cè)的流體環(huán)境產(chǎn)生擾動，嚴(yán)重時還伴有氣泡和漩渦產(chǎn)生.在文中的計算模型中，浮體材料為結(jié)構(gòu)鋼，彈性模量E=2.1×1011Pa，因此結(jié)構(gòu)幾乎沒有變形產(chǎn)生.但是在t=7.6s左右時，上浪部分水體流經(jīng)甲板時，結(jié)構(gòu)的底部與甲板的中部所受應(yīng)力較大，應(yīng)引起注意采取結(jié)構(gòu)的加強(qiáng)措施.

4 結(jié)論

基于ANSYS Workbench平臺的流固耦合模塊，利用UDF程序進(jìn)行二次開發(fā)實現(xiàn)速度入口造波，采用VOF方法追蹤氣液兩項自由面，在水池后方采用動量方程中添加源項形成阻尼消波，成功建立二維數(shù)值波浪水池，實現(xiàn)了波浪與結(jié)構(gòu)物相互作用的完全非線性數(shù)值模擬.

1)將設(shè)置造波邊界法應(yīng)用于Fluent模塊，成功模擬了效果良好的三階stokes波.消波區(qū)域采用阻尼消波，避免了波浪到達(dá)右端固壁邊界后產(chǎn)生反射波.

圖9 波浪與自由運動浮體的相互作用Fig.9 Interaction between wave and freely floating structures

2)采用雙向流固耦合的方法，分別對不同干舷高度的結(jié)構(gòu)物在不同波浪要素的波浪作用下的載荷與響應(yīng)進(jìn)行了對比研究，得到了結(jié)構(gòu)物遭遇波浪后的運動響應(yīng)和監(jiān)測點的上浪水位高度，發(fā)現(xiàn)結(jié)構(gòu)物的運動響應(yīng)幅值和甲板上浪水位高度都會隨著波陡的增加而增加.同時，增加結(jié)構(gòu)物干舷的高度能在一定程度上減少甲板上浪，但會使其運動響應(yīng)幅值尤其是橫蕩幅值增加，在結(jié)構(gòu)設(shè)計時應(yīng)加以權(quán)衡.

3)展示了波浪與浮體相互作用時的流體形態(tài)和結(jié)構(gòu)應(yīng)力云圖，觀察到了波浪與結(jié)構(gòu)物相遇后波浪沿船首爬升、破碎、飛濺等完全非線性現(xiàn)象.并注意到上浪部分水體流經(jīng)結(jié)構(gòu)物甲板時，結(jié)構(gòu)底部與甲板的中部所受應(yīng)力較大，應(yīng)注意采取結(jié)構(gòu)的加強(qiáng)措施.

References)

［1］ Ma Qingwei.Advances in numerical simulation of nonlinear water waves［M］.Singapore:World Scientific Publishing Co Pte Ltd，2010:1 -39.

［2］ Longuet·Higgins M S，Cokelet E D.The deformation of steep surface waves on water:a numerical method of computation［J］.Proceedings of the Royal Society of London，Series A 350，1976(1660):1 -26.

［3］ Zhou B Z，Wu G X，Teng B.Fully nonlinear wave interaction with freely floating non-wall-sided structures［J］.Engineering Analysis with Boundary Elements，2015(50):117-132.

［4］ Ma Q W，Yan S.Quasi ALE finite element method for nonlinear water waves［J］.Journal of Computational Physics，2006，212:52 -72.

［5］ Yan S，Ma Q W.Numerical simulation of fully nonlinear interaction between steep waves and 2D floating bodies using the QALE-FEM method［J］.Journal of Computational Physics，2007，221(2):666 -692.

［6］ Yan S，Ma Q W.Numerical simulation of interaction between wind and 2D freak waves［J］.European Journal of Mechanics B/Fluids，2010，29(1):18 -31.

［7］ Hirt C W，Nichols B D.Volume of fluid(VOF)method for the dynamics of free boundaries［J］.Journal of Computational Physics，1981，39(1):201 -225.

［8 ］ Stanley Osher，James A Sethian.Fronts propagating with curvature-dependent speed:Algorithms based on Hamilton-Jacobi formulation ［J］.Journal of Computational Physics，1988(01):12 -49.

［9］ Monaghan J J.Smoothed particle hydrodynamics［J］.Reports on Progress in Physics，2005，68(8):543 -574.

［10］ Yabe T，Xiao F，Utsumi T.The constrained interpolation profile method for multiphase analysis［J］.Journal of Computational Physics，2001，69(2):556 -593.

［11］ 朱仁慶，王志東，楊松林.完全非線性波的數(shù)值模擬［J］.船舶力學(xué)，2002，6(5):14-18.Zhu Renqing，Wang Zhidong，Yang Songlin.Fully nonlinear wave simulations［J］.Journal of Ship Mechanics，2002，6(5):14 -18.(in Chinese)

［12］ Nielsen K B，Mayer S.Numerical prediction of green water incidents［J］.Ocean Engineering，2004，31(3 -4):363-399.

［13］ Kleefsman K M，Loots G E.The numerical simulation of green water loading including vessel motions and the incoming wave field［C］∥Proceedings of the 24th International Conference on Offshore Mechanics and Arctic Engineering.Halkidiki，Greece:OMAE，2005:981 -992.

［14］ 朱仁傳，繆國平，林兆偉.運動船體甲板上浪的三維數(shù)值模擬［J］.水動力學(xué)研究與進(jìn)展A輯，2008，23(1):7-14.Zhu Renchuan，Miao Guoping，Lin Zhaowei.3-D numerical simulation of green water occurrence on oscillating shuip［J］.Journal of Hydrodinamics Ser A，2008，23(1):7-14.(in Chinese)

［15］ 陶建華.水波的數(shù)值模擬［M］.天津:天津大學(xué)出版社，2005:17-252.

［16］ 朱仁慶，李辰，顧思琪，等.彈性液艙內(nèi)液體晃蕩研究［J］.江蘇科技大學(xué)學(xué)報:自然科學(xué)版，2013，27(3):214-218.Zhu Renqing，Li Chen，Gu Siqi，et al.A study on liquid sloshing in elastic tanks［J］.Journal of Jiangsu U-niversity of Science and Technology:Natural Science E-dition，2013，27(3):214 -218.(in Chinese)

［17］ Greco M，F(xiàn)altinsen O M，Landrini M.Green water loading on a deck structure［C］∥Proceedings of 16th International Workshop on Water Waves and Floating Bodies.Japan:IWWWFB，2001.