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

顧思琪，朱仁慶，吳梓鑫

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

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

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

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

1 物理模型

1.1 流體運動控制方程

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

連續(xù)方程

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

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

動量方程

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

1.2 自由液面追蹤

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

1.3 波浪理論

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

波面方程

x方向速度分量

z方向速度分量

2 數(shù)值模型

2.1計算模型

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

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

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

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

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

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

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

2.3 消波設置

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

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

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

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

3 結果與分析

3.1 波浪數(shù)值水池

穩(wěn)定的波浪是波物相互作用數(shù)值模擬的前提，為了驗證波浪數(shù)值水池的正確性，文中首先對波浪周期T=1.13s，理論波長λ=2 m，理論波高H=0.16 m的三階stokes波進行了數(shù)值模擬.圖2給出了數(shù)值波浪水池在x=1 m和2 m處波面(A)隨時間的變化歷程曲線.在波浪的傳播過程中，隨著傳播距離的增加，波高是有衰減的，這是由于流體粘性作用的影響導致波浪能量耗散的緣故，因此，圖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)下二維固定結構物甲板上浪過程進行數(shù)值模擬，并將計算結果與有關文獻中的試驗值進行比較，以驗證數(shù)值方法的有效性.波物相互作用問題是雙向流固耦合問題，文中基于軟件ANSYS，采用迭代耦合的方法求解流固耦合問題，其主要思路是流體方程和結構方程按順序相互迭代求解，各自在每一步得到的結果提供給另一部分使用，直到耦合系統(tǒng)的解達到收斂，迭代停止［16］.

1)計算模型說明

數(shù)值波浪水池和結構模型如圖4，模型尺寸與文獻［16］中甲板上浪試驗模型保持一致.圖中水池長13.5 m，其中2 m為消波區(qū)，水深1.035 m.結構模型已簡化為方形結構物，底部為半徑0.08 m的圓弧，長和型深分別為1.5 m和0.248 m，吃水為0.198 m，干舷為0.05 m.

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

2)計算結果

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結論

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

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

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

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

3)展示了波浪與浮體相互作用時的流體形態(tài)和結構應力云圖，觀察到了波浪與結構物相遇后波浪沿船首爬升、破碎、飛濺等完全非線性現(xiàn)象.并注意到上浪部分水體流經(jī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］.船舶力學，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］.水動力學研究與進展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］.天津:天津大學出版社，2005:17-252.

［16］ 朱仁慶，李辰，顧思琪，等.彈性液艙內(nèi)液體晃蕩研究［J］.江蘇科技大學學報:自然科學版，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.