陳寶清，張金鳳，史小康

（天津大學(xué) 水利工程安全與仿真國(guó)家重點(diǎn)實(shí)驗(yàn)室，天津 300350）

基于OpenFOAM的波浪作用下海床動(dòng)力響應(yīng)

陳寶清，張金鳳，史小康

（天津大學(xué) 水利工程安全與仿真國(guó)家重點(diǎn)實(shí)驗(yàn)室，天津 300350）

基于計(jì)算流體力學(xué)開(kāi)源軟件OpenFOAM建立了波浪-海床-結(jié)構(gòu)物相互作用的數(shù)學(xué)模型。其中，數(shù)值波浪水槽采用動(dòng)量源函數(shù)造波，求解雷諾平均的N-S方程，VOF方法追蹤自由表面；海床動(dòng)力響應(yīng)基于Biot固結(jié)理論，并結(jié)合摩爾-庫(kù)倫塑性本構(gòu)模型。首先對(duì)波浪作用下海床動(dòng)力響應(yīng)問(wèn)題進(jìn)行研究，數(shù)模結(jié)果與實(shí)驗(yàn)結(jié)果及解析解吻合良好。進(jìn)一步在數(shù)值水槽中添加孔隙介質(zhì)模型，模擬潛堤周?chē)４矂?dòng)力響應(yīng)問(wèn)題，結(jié)果表明本模型可以較好地模擬出考慮結(jié)構(gòu)物影響的海床動(dòng)力響應(yīng)過(guò)程。

波浪作用；結(jié)構(gòu)物；海床；動(dòng)力響應(yīng)；OpenFOAM

1 研究背景

在循環(huán)波浪力作用下，海床可能會(huì)出現(xiàn)變形、軟化、承載力下降、液化等現(xiàn)象，危及其上部結(jié)構(gòu)物的安全。而近年來(lái)，許多近海結(jié)構(gòu)物的失穩(wěn)破壞案例與海床液化有關(guān)[1]，因此波浪作用下的海床響應(yīng)問(wèn)題不可忽視。

海床液化造成的經(jīng)濟(jì)損失往往在結(jié)構(gòu)物破壞方面，因此，建立波浪-結(jié)構(gòu)物-海床相互作用（Wave-Structure-Seabed interaction，簡(jiǎn)稱(chēng)WSSI）的集成模型更具有實(shí)際意義。Mase[2]基于Biot固結(jié)方程，用有限元法建立了研究立波作用下，沉箱防波堤下的砂質(zhì)海床響應(yīng)模型。隨后，Mizutani和Mostafa[3]通過(guò)建立BEM-FEM耦合模型來(lái)研究波浪、海床和結(jié)構(gòu)物間相互作用的問(wèn)題；在耦合模型中，采用泊松方程來(lái)控制波浪場(chǎng)，而結(jié)構(gòu)及海床的控制方程則為Biot固結(jié)方程。Ulker[4]考慮土體顆粒與流體粒子間的相對(duì)位移，建立了求解WSSI問(wèn)題的有限元模型。Jeng[5]將雷諾平均的N-S方程（Reynolds-Averaged Navier-Stokes，簡(jiǎn)稱(chēng)RANS）作為流場(chǎng)及流體在防波堤中孔隙流動(dòng)的控制方程，海床的控制方程為動(dòng)態(tài)Biot方程，建立了集成模型求解WSSI問(wèn)題。

近年來(lái)，基于計(jì)算流體力學(xué)（Computational Fluid Dynamics，簡(jiǎn)稱(chēng)CFD）的開(kāi)源軟件OpenFOAM的發(fā)展為求解海床動(dòng)力響應(yīng)問(wèn)題提供了新的途徑。Liu[6]通過(guò)有限體積法求解RANS方程和Biot固結(jié)方程，并采用邊界條件造波的方式，構(gòu)建三維數(shù)值模型來(lái)研究海床在波浪作用下的瞬時(shí)液化問(wèn)題。Tang[7]對(duì)Liu的計(jì)算結(jié)構(gòu)力學(xué)（Computational Structural Dynamics，簡(jiǎn)稱(chēng)CSD）模型進(jìn)行改進(jìn)，使其適用于土體存在塑性變形的情況。Safti[8]通過(guò)求解全耦合的動(dòng)態(tài)Biot方程，研究孔隙水壓力在防波堤及海床中的變化情況。于聰?shù)萚9]利用OpenFOAM建立數(shù)值波浪水槽，將海床受到的波壓力作為初始條件輸入有限元軟件ABAQUS，模擬波浪作用下海上風(fēng)電基礎(chǔ)周?chē)４矂?dòng)力響應(yīng)，并沒(méi)有實(shí)現(xiàn)波浪和海床的相互作用。本文將基于開(kāi)源軟件OpenFOAM建立CFD-CSD耦合模型，采用源函數(shù)造波的方式，考慮土體的塑性變形情況，求解波浪作用下海床動(dòng)力響應(yīng)及考慮結(jié)構(gòu)物的WSSI問(wèn)題，并與已有的解析解和實(shí)驗(yàn)值作對(duì)比，以驗(yàn)證模型的合理性。

2 控制方程

本模型通過(guò)耦合OpenFOAM中的兩相流模型與Biot固結(jié)模型構(gòu)建CFD-CSD單向耦合模型，并在兩相流模型中添加動(dòng)量源項(xiàng)進(jìn)行造波。其中CFD模型的控制方程采用RANS方程，CSD模型的控制方程為Biot固結(jié)方程。

2.1 CFD模型控制方程

OpenFOAM中的兩相流模型是通過(guò)求解RANS方程，選取RNG k-ε紊流模型，對(duì)水和空氣兩相不可壓縮流體的運(yùn)動(dòng)進(jìn)行模擬，并通過(guò)VOF法對(duì)自由表面進(jìn)行捕捉。其控制方程如下：式（1）為連續(xù)性方程，式（2）為動(dòng)量方程，式（3）為VOF法求解體積函數(shù)的對(duì)流方程。方程中各個(gè)參數(shù)的物理意義如下：ρ為流體的密度；U是速度矢量；μeff是考慮了分子動(dòng)力黏性和紊流作用的渦黏性系數(shù)；p_rgh是修正壓力項(xiàng)；ρS為分布造波源項(xiàng)；ρrU為消波項(xiàng)，其中r為消波系數(shù)；X為位置矢量；g為重力加速度；S為造波源函數(shù)；γ為流體體積函數(shù)，取值在0～1之間，等于0表示該單元完全被空氣填滿，等于1代表單元全部被流體填滿，在0到1之間說(shuō)明單元內(nèi)部是空氣和水的混合體。

2.2 CSD模型控制方程

CSD模型基于Biot多孔介質(zhì)模型及固結(jié)理論，并考慮土體的塑性變形，對(duì)孔隙水在多孔介質(zhì)模型中的流動(dòng)進(jìn)行模擬。其控制方程如下：

式（4）為穩(wěn)態(tài)條件下的動(dòng)量方程，式（5）為滲流方程。式中：σ為總應(yīng)力張量；σ′為有效應(yīng)力張量；p為孔隙水壓力；I為單位張量；ks為滲流系數(shù)；n為孔隙率；u為土骨架位移矢量；K′為孔隙水的有效體積模量，可由飽和度Sr求得。

根據(jù)微小應(yīng)變及各項(xiàng)同性假定，彈塑性土體的本構(gòu)關(guān)系如下：

式中：ε為總的應(yīng)變張量；εp為塑性應(yīng)變張量；μ和λs為L(zhǎng)amé系數(shù)，可以由楊氏模量E及泊松比ν表示。

將本構(gòu)方程（6）、（7）代入動(dòng)量方程（4），得到u-p形式的動(dòng)量方程：

2.3 CFD-CSD耦合模型

CFD-CSD模型的耦合是通過(guò)OpenFOAM中的多區(qū)域模塊，對(duì)流體區(qū)域與土體區(qū)域分別進(jìn)行求解，并在耦合邊界處（海床表面）通過(guò)mapped-Field邊界條件將流體求解器產(chǎn)生的波壓力傳遞給土體求解器。耦合模型的主要求解流程如圖1所示。壓力泊松方程通過(guò)兩個(gè)模塊耦合求解；將動(dòng)量方程（2）寫(xiě)成如下形式：auPUP=Hu（U）-▽p_rgh-gX▽?duì)?（9）式中：下標(biāo)P表示單元P上的量；UP是單元P的速度矢量；auP是UP的系數(shù)；Hu（U）是離散單元p時(shí)除了UP的其他單元速度。再將式（9）代入連續(xù)性方程（1）得到壓力泊松方程：

▽·{[auP]-1▽p_rgh}=▽·{[auP]-1[Hu（U）-g·X▽?duì)裖}（10）

3 模型驗(yàn)證

為驗(yàn)證CFD-CSD單向耦合模型的合理性，需要與前人所做實(shí)驗(yàn)結(jié)果及解析解作對(duì)比。以下從簡(jiǎn)單波浪條件下及考慮結(jié)構(gòu)物影響的情況下分別對(duì)模型進(jìn)行驗(yàn)證。

3.1 波浪作用下海床響應(yīng)

首先，CFD-CSD模型計(jì)算結(jié)果與別社安[10]實(shí)驗(yàn)結(jié)果進(jìn)行比較。模擬的數(shù)值水槽長(zhǎng)度為30 m，海床厚度為30 cm，長(zhǎng)度為10 m；A、B、C 3個(gè)測(cè)點(diǎn)分別位于海床表面以下1.5 cm、6.5 cm及11.5 cm；流體網(wǎng)格在波面處進(jìn)行加密以保證波浪模擬的精確性，土體網(wǎng)格劃分為寬2 cm、高1 cm的控制體單元。波浪參數(shù)及土體參數(shù)如表1所示，不同測(cè)點(diǎn)處的海床孔隙水壓力時(shí)間歷程如圖2所示，其中散點(diǎn)為別社安[10]實(shí)驗(yàn)結(jié)果，曲線為數(shù)值模擬結(jié)果。圖2結(jié)果表明孔隙水壓力沿垂向呈現(xiàn)遞減趨勢(shì)，CFD-CSD數(shù)值模擬結(jié)果同試驗(yàn)值有良好的一致性。

同時(shí)為進(jìn)一步驗(yàn)證本文模型的合理性，將數(shù)模結(jié)果與Jeng[12]解析解及Yamamoto[11]實(shí)驗(yàn)值進(jìn)行比較。模擬的海床厚度為50 cm，T=1.0 s時(shí)，海床長(zhǎng)度為7 m，波浪數(shù)值水槽長(zhǎng)度為25 m；T= 1.5 s時(shí)，海床長(zhǎng)度為16 m，水槽長(zhǎng)度為50 m。數(shù)值模擬結(jié)果中海床超靜孔隙水壓力沿垂向分布的結(jié)果如圖3所示，采用的波浪參數(shù)和土體參數(shù)如表2所示。從圖中可以看出，CFD-CSD單向耦合模型的數(shù)值模擬結(jié)果和Jeng解析解比較吻合。

3.2 波浪-結(jié)構(gòu)-海床相互作用

從上述2組驗(yàn)證結(jié)果可以看出：波浪條件下，模擬結(jié)果與實(shí)驗(yàn)結(jié)果及解析解吻合良好。然而，當(dāng)波浪場(chǎng)中有結(jié)構(gòu)物存在時(shí)，結(jié)構(gòu)與波浪間相互作用，使波浪場(chǎng)的情況復(fù)雜化，同時(shí)結(jié)構(gòu)自重對(duì)海床土體存在一定的作用。因此本節(jié)將要驗(yàn)證耦合模型在波浪、結(jié)構(gòu)物、海床三者相互作用情況下的合理性。

驗(yàn)證算例為Mizutani及Mostafa[3]所做的潛堤液化實(shí)驗(yàn)。數(shù)值模型設(shè)置如圖4所示，其中a～d點(diǎn)為流速測(cè)點(diǎn)，A～D為孔隙水壓力測(cè)點(diǎn)。模型的波浪參數(shù)及土體參數(shù)如表3所示。等孔隙結(jié)構(gòu)的處理方式是在孔隙區(qū)域內(nèi)，在式（2）右側(cè)添加阻力項(xiàng)R。阻力項(xiàng)R表達(dá)式為：

潛堤實(shí)質(zhì)上為孔隙結(jié)構(gòu)，與不可滲透的墻體結(jié)構(gòu)不同，流體運(yùn)動(dòng)可以直接穿透孔隙層，并在其內(nèi)部形成復(fù)雜的孔隙流動(dòng)。本文模型中對(duì)潛堤

式中：v為動(dòng)力黏滯系數(shù)；n為孔隙率；D50為土顆粒中值粒徑；α、β分別為線性項(xiàng)與非線性項(xiàng)系數(shù)。

測(cè)點(diǎn)處波面歷時(shí)曲線及孔隙水壓力歷時(shí)曲線如圖5、圖6所示。圖中η為波面相對(duì)于靜水位的高度；H為波高；Ps為測(cè)點(diǎn)處孔隙水壓力；γ為海水重度。從圖中可以看出，數(shù)模計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果[3]擬合較好，故本文模型在處理波浪-結(jié)構(gòu)物-海床相互作用的問(wèn)題上適用性較好。

4 結(jié)語(yǔ)

本文通過(guò)耦合開(kāi)源流體力學(xué)軟件OpenFOAM中的兩相流求解器interFoam及多孔彈塑性土體應(yīng)力求解器poroPlasticStressedFoam，并在兩相流求解器中添加動(dòng)量源以實(shí)現(xiàn)造波功能，建立了分析波浪作用下海床動(dòng)力響應(yīng)的三維數(shù)值模型。為驗(yàn)證模型的有效性，首先模擬了波浪作用下海床動(dòng)力響應(yīng)過(guò)程，得到的數(shù)值模型計(jì)算結(jié)果與實(shí)驗(yàn)值及解析解吻合較好。通過(guò)在兩相流求解器中添加孔隙層模塊，進(jìn)一步模擬了流場(chǎng)中存在結(jié)構(gòu)物時(shí)的海床動(dòng)力響應(yīng)過(guò)程，模擬得到的結(jié)果與實(shí)驗(yàn)值相比有較好的一致性。CFD-CSD耦合數(shù)學(xué)模型的開(kāi)發(fā)可以為進(jìn)一步研究海岸工程和海上風(fēng)電基礎(chǔ)周?chē)４岔憫?yīng)及液化問(wèn)題提供技術(shù)手段。

[1]DEMARS K R,VANOVER E A.Measurement of wave-induced pressures and stresses in a sand bed[J].Marine Geotechnology, 1985,6(1):29-59.

[2]MASE H,SAKAI T,SAKAMOTO M.Wave-induced pore water pressure and effective stresses around breakwater[J].Ocean Engineering,1994,21(4):361-379.

[3] MIZUTANI N,MOSTAFA A M.Nonlinear wave-induced seabed instability around coastalstructures[J].Coastal Engineering Journal,1998,40(2):131-160.

[4]ULKER M,RAHMAN M S,GUDDATIMN.Wave-induced dynamic response and instability of seabed around caisson breakwater[J].Ocean Engineering,2010,37(17-18):1 522-1 545.

[5] YE Jian-hong,JENG Dong-sheng,WANG Ren,et al.Numerical simulation of the wave-induced dynamic response of poro-elastoplastic seabed foundations and a composite breakwater[J].Applied Mathematical Modelling,2014,39(1):322-347.

[6] LIU Xiao-feng.Numerical models for scour and liquefaction around object under currents and waves[D].Illinois:University of Illinois at Urbana-Champaign,2008.

[7]TANG T,HEDEDAL O,CARDIFF P.On finite volume method implementation of poro-elasto-plasticity soil model[J].International Journal for Numerical and Analytical Methods in Geomechanics, 2015,39(13):1 410-1 430.

[8] ELSAFTI Hisham.A Numerical wave-structure-soil interaction model for monolithic breakwaters subject to breaking wave impact [C]//13th Triennial International Conference,ASCE,2013:1 974-1 984.

[9]于聰，張金鳳，張慶河.復(fù)合筒型基礎(chǔ)周?chē)馏w在波浪作用下的動(dòng)力響應(yīng)[J].中國(guó)港灣建設(shè)，2015，35（1）：1-5. YU Cong,ZHANG Jin-feng,ZHANG Qing-he.Wave-induced seabed dynamic response around composite bucket foundation[J]. China Harbour Engineering,2015,35(1):1-5.

[10]別設(shè)安.波浪作用下海床中的孔隙水應(yīng)力響應(yīng)及其對(duì)建筑物穩(wěn)定性的影響[D].天津：天津大學(xué)，1996. BIE She-an.Pore pressure response in seabed under wave and its influence on building stability[D].Tianjin:Tianjin University, 1996.

[11]YAMAMOTO T,KONING H L,SELLMEIHER H,et al.On the response of a poro-elastic bed to water waves[J].Journal of Fluid Mechanics,1978,87(1):193-206.

[12]LEE T C,TSAI C P,JENG D S.Ocean waves propagating over a porous seabed of finite thickness[J].Ocean Engineering,2002,29 (12):1 577-1 601.

Numerical simulation for the dynamic response of seabed under waves based on the OpenFOAM

CHEN Bao-qing,ZHANG Jin-feng,SHI Xiao-kang
（State Key Laboratory of Hydraulic Engineering Simulation and Safety,Tianjin University,Tianjin 300350,China）

We established a numerical model of wave-seabed-structure interaction based on the open source software OpenFOAM of computational fluid dynamics.The numerical wave tank is governed by the Reynolds-Averaged Navier-Stokes equations and applied the momentum source function to generate waves.The volume of fluid method is used to track the water free-surface.The dynamic response model of seabed is formulated on the basis of Biot's consolidation theory and combined with the plasticity Mohr-Coulomb constitutive relation.Firstly,the dynamic response of seabed under wave action is investigated, and the results agreed well with the experimental and analytical results.Furthermore,the porous media model was added to simulate the dynamic response of the seabed around the submerged breakwater.The results show that the model can well simulate the dynamic response of the seabed considering the influence of the structure.

wave action；structure；seabed；dynamic response；OpenFOAM

U652.21；P753

A

2095-7874（2017）03-0001-05

10.7640/zggwjs201703001

2016-10-24

2016-11-25

國(guó)家自然科學(xué)基金創(chuàng)新研究群體科學(xué)基金（51321065）；天津市應(yīng)用基礎(chǔ)及前沿技術(shù)研究計(jì)劃（青年基金項(xiàng)目）（12JCQNJC04100）

陳寶清（1991— ），男，福建福州人，碩士研究生，主要從事海床在波浪作用下動(dòng)力響應(yīng)研究。E-mail：350625871@qq.com