基于OpenFOAM的波浪作用下海床动力响应

陈宝清，张金凤，史小康

（天津大学 水利工程安全与仿真国家重点实验室，天津 300350）

基于OpenFOAM的波浪作用下海床动力响应

陈宝清，张金凤，史小康

（天津大学 水利工程安全与仿真国家重点实验室，天津 300350）

基于计算流体力学开源软件OpenFOAM建立了波浪-海床-结构物相互作用的数学模型。其中，数值波浪水槽采用动量源函数造波，求解雷诺平均的N-S方程，VOF方法追踪自由表面；海床动力响应基于Biot固结理论，并结合摩尔-库伦塑性本构模型。首先对波浪作用下海床动力响应问题进行研究，数模结果与实验结果及解析解吻合良好。进一步在数值水槽中添加孔隙介质模型，模拟潜堤周围海床动力响应问题，结果表明本模型可以较好地模拟出考虑结构物影响的海床动力响应过程。

波浪作用；结构物；海床；动力响应；OpenFOAM

1 研究背景

在循环波浪力作用下，海床可能会出现变形、软化、承载力下降、液化等现象，危及其上部结构物的安全。而近年来，许多近海结构物的失稳破坏案例与海床液化有关[1]，因此波浪作用下的海床响应问题不可忽视。

海床液化造成的经济损失往往在结构物破坏方面，因此，建立波浪-结构物-海床相互作用（Wave-Structure-Seabed interaction，简称WSSI）的集成模型更具有实际意义。Mase[2]基于Biot固结方程，用有限元法建立了研究立波作用下，沉箱防波堤下的砂质海床响应模型。随后，Mizutani和Mostafa[3]通过建立BEM-FEM耦合模型来研究波浪、海床和结构物间相互作用的问题；在耦合模型中，采用泊松方程来控制波浪场，而结构及海床的控制方程则为Biot固结方程。Ulker[4]考虑土体颗粒与流体粒子间的相对位移，建立了求解WSSI问题的有限元模型。Jeng[5]将雷诺平均的N-S方程（Reynolds-Averaged Navier-Stokes，简称RANS）作为流场及流体在防波堤中孔隙流动的控制方程，海床的控制方程为动态Biot方程，建立了集成模型求解WSSI问题。

近年来，基于计算流体力学（Computational Fluid Dynamics，简称CFD）的开源软件OpenFOAM的发展为求解海床动力响应问题提供了新的途径。Liu[6]通过有限体积法求解RANS方程和Biot固结方程，并采用边界条件造波的方式，构建三维数值模型来研究海床在波浪作用下的瞬时液化问题。Tang[7]对Liu的计算结构力学（Computational Structural Dynamics，简称CSD）模型进行改进，使其适用于土体存在塑性变形的情况。Safti[8]通过求解全耦合的动态Biot方程，研究孔隙水压力在防波堤及海床中的变化情况。于聪等[9]利用OpenFOAM建立数值波浪水槽，将海床受到的波压力作为初始条件输入有限元软件ABAQUS，模拟波浪作用下海上风电基础周围海床动力响应，并没有实现波浪和海床的相互作用。本文将基于开源软件OpenFOAM建立CFD-CSD耦合模型，采用源函数造波的方式，考虑土体的塑性变形情况，求解波浪作用下海床动力响应及考虑结构物的WSSI问题，并与已有的解析解和实验值作对比，以验证模型的合理性。

2 控制方程

本模型通过耦合OpenFOAM中的两相流模型与Biot固结模型构建CFD-CSD单向耦合模型，并在两相流模型中添加动量源项进行造波。其中CFD模型的控制方程采用RANS方程，CSD模型的控制方程为Biot固结方程。

2.1 CFD模型控制方程

OpenFOAM中的两相流模型是通过求解RANS方程，选取RNG k-ε紊流模型，对水和空气两相不可压缩流体的运动进行模拟，并通过VOF法对自由表面进行捕捉。其控制方程如下：式（1）为连续性方程，式（2）为动量方程，式（3）为VOF法求解体积函数的对流方程。方程中各个参数的物理意义如下：ρ为流体的密度；U是速度矢量；μeff是考虑了分子动力黏性和紊流作用的涡黏性系数；p_rgh是修正压力项；ρS为分布造波源项；ρrU为消波项，其中r为消波系数；X为位置矢量；g为重力加速度；S为造波源函数；γ为流体体积函数，取值在0～1之间，等于0表示该单元完全被空气填满，等于1代表单元全部被流体填满，在0到1之间说明单元内部是空气和水的混合体。

2.2 CSD模型控制方程

CSD模型基于Biot多孔介质模型及固结理论，并考虑土体的塑性变形，对孔隙水在多孔介质模型中的流动进行模拟。其控制方程如下：

式（4）为稳态条件下的动量方程，式（5）为渗流方程。式中：σ为总应力张量；σ′为有效应力张量；p为孔隙水压力；I为单位张量；ks为渗流系数；n为孔隙率；u为土骨架位移矢量；K′为孔隙水的有效体积模量，可由饱和度Sr求得。

根据微小应变及各项同性假定，弹塑性土体的本构关系如下：

式中：ε为总的应变张量；εp为塑性应变张量；μ和λs为Lamé系数，可以由杨氏模量E及泊松比ν表示。

将本构方程（6）、（7）代入动量方程（4），得到u-p形式的动量方程：

2.3 CFD-CSD耦合模型

CFD-CSD模型的耦合是通过OpenFOAM中的多区域模块，对流体区域与土体区域分别进行求解，并在耦合边界处（海床表面）通过mapped-Field边界条件将流体求解器产生的波压力传递给土体求解器。耦合模型的主要求解流程如图1所示。压力泊松方程通过两个模块耦合求解；将动量方程（2）写成如下形式：auPUP=Hu（U）-▽p_rgh-gX▽ρ （9）式中：下标P表示单元P上的量；UP是单元P的速度矢量；auP是UP的系数；Hu（U）是离散单元p时除了UP的其他单元速度。再将式（9）代入连续性方程（1）得到压力泊松方程：

▽·{[auP]-1▽p_rgh}=▽·{[auP]-1[Hu（U）-g·X▽ρ]}（10）

3 模型验证

为验证CFD-CSD单向耦合模型的合理性，需要与前人所做实验结果及解析解作对比。以下从简单波浪条件下及考虑结构物影响的情况下分别对模型进行验证。

3.1 波浪作用下海床响应

首先，CFD-CSD模型计算结果与别社安[10]实验结果进行比较。模拟的数值水槽长度为30 m，海床厚度为30 cm，长度为10 m；A、B、C 3个测点分别位于海床表面以下1.5 cm、6.5 cm及11.5 cm；流体网格在波面处进行加密以保证波浪模拟的精确性，土体网格划分为宽2 cm、高1 cm的控制体单元。波浪参数及土体参数如表1所示，不同测点处的海床孔隙水压力时间历程如图2所示，其中散点为别社安[10]实验结果，曲线为数值模拟结果。图2结果表明孔隙水压力沿垂向呈现递减趋势，CFD-CSD数值模拟结果同试验值有良好的一致性。

同时为进一步验证本文模型的合理性，将数模结果与Jeng[12]解析解及Yamamoto[11]实验值进行比较。模拟的海床厚度为50 cm，T=1.0 s时，海床长度为7 m，波浪数值水槽长度为25 m；T= 1.5 s时，海床长度为16 m，水槽长度为50 m。数值模拟结果中海床超静孔隙水压力沿垂向分布的结果如图3所示，采用的波浪参数和土体参数如表2所示。从图中可以看出，CFD-CSD单向耦合模型的数值模拟结果和Jeng解析解比较吻合。

3.2 波浪-结构-海床相互作用

从上述2组验证结果可以看出：波浪条件下，模拟结果与实验结果及解析解吻合良好。然而，当波浪场中有结构物存在时，结构与波浪间相互作用，使波浪场的情况复杂化，同时结构自重对海床土体存在一定的作用。因此本节将要验证耦合模型在波浪、结构物、海床三者相互作用情况下的合理性。

验证算例为Mizutani及Mostafa[3]所做的潜堤液化实验。数值模型设置如图4所示，其中a～d点为流速测点，A～D为孔隙水压力测点。模型的波浪参数及土体参数如表3所示。等孔隙结构的处理方式是在孔隙区域内，在式（2）右侧添加阻力项R。阻力项R表达式为：

潜堤实质上为孔隙结构，与不可渗透的墙体结构不同，流体运动可以直接穿透孔隙层，并在其内部形成复杂的孔隙流动。本文模型中对潜堤

式中：v为动力黏滞系数；n为孔隙率；D50为土颗粒中值粒径；α、β分别为线性项与非线性项系数。

测点处波面历时曲线及孔隙水压力历时曲线如图5、图6所示。图中η为波面相对于静水位的高度；H为波高；Ps为测点处孔隙水压力；γ为海水重度。从图中可以看出，数模计算结果与实验结果[3]拟合较好，故本文模型在处理波浪-结构物-海床相互作用的问题上适用性较好。

4 结语

本文通过耦合开源流体力学软件OpenFOAM中的两相流求解器interFoam及多孔弹塑性土体应力求解器poroPlasticStressedFoam，并在两相流求解器中添加动量源以实现造波功能，建立了分析波浪作用下海床动力响应的三维数值模型。为验证模型的有效性，首先模拟了波浪作用下海床动力响应过程，得到的数值模型计算结果与实验值及解析解吻合较好。通过在两相流求解器中添加孔隙层模块，进一步模拟了流场中存在结构物时的海床动力响应过程，模拟得到的结果与实验值相比有较好的一致性。CFD-CSD耦合数学模型的开发可以为进一步研究海岸工程和海上风电基础周围海床响应及液化问题提供技术手段。

[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]于聪，张金凤，张庆河.复合筒型基础周围土体在波浪作用下的动力响应[J].中国港湾建设，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]别设安.波浪作用下海床中的孔隙水应力响应及其对建筑物稳定性的影响[D].天津：天津大学，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

国家自然科学基金创新研究群体科学基金（51321065）；天津市应用基础及前沿技术研究计划（青年基金项目）（12JCQNJC04100）

陈宝清（1991— ），男，福建福州人，硕士研究生，主要从事海床在波浪作用下动力响应研究。E-mail：350625871@qq.com