基于OpenFOAM波浪-植物相互作用的数值模拟研究

张宸浩，张明亮*，2，柴崇顼，王 珏，王清夷

(1.大连海洋大学海洋科技与环境学院，辽宁 大连 116023； 2.大连理工大学海岸和近海工程国家重点实验室，辽宁 大连 116025)

0 引言

波浪从远海传播到近岸的过程中受海床地形变化的影响，波能集中在较小水深内，波高增加，会对海岸造成侵蚀破坏作用。通常采用防波堤抵御波浪来保护堤岸，但防波堤大多由水泥、混凝土建成，不仅观赏价值较低，还会破坏原始生态环境。近岸潮间带区域多生长盐沼植物，这些水生植物能改变波浪形态，在一定程度上消减波浪能量，并改变滩上的水流、波浪等动力条件，进而影响近海海床形态的变化[1]。由于盐沼植物能为海岸提供天然的生态屏障，近年来，关于生态型护岸工程与波浪-植物相互作用的研究备受关注。

关于近岸海域水动力与植物相互作用的研究主要包括现场观测[2-3]、理论推导[4-5]、物理实验[6-7]与数值模拟[8-9]。其中，数值模拟实施成本低，不受现场条件、实验尺度等因素的限制，且可以进行定量分析，已被广泛应用于近岸区域的水动力研究。如PHAN et al[10]基于SWASH建立了可模拟规则波与不规则波的二维数值模型，通过耦合Morrison方程添加植物作用，研究了植物对非线性波浪消波效果的影响，并与实验结果进行对比，得到了波长、波衰减率与Ursell数的函数关系。WANG et al[11]以OlaFlow求解器为基础，建立了可模拟孤立波的二维数值模型，采用宏观方法分析了有无植物作用下孤立波的传播与衰减，结果表明孤立波在含植物水槽内的衰减效果大于无植物工况。ZHANG et al[12]采用有限差分法求解雷诺平均纳维斯托克斯方程，用k-ε湍流模型封闭方程组，对植物作用下的波浪传播进行了模拟分析，结果显示波浪经过植物带后出现波高衰减与能量耗散。PAQUIER et al[13]使用GPUSPH方法建立三维模型模拟柔性海草对波浪的衰减作用，并还原了海草的运动轨迹，但受计算效率与粒子数的限制，只进行了小范围的模拟计算。WANG et al[14]通过建立包含循环边界条件的三维数值模型将刚性植物概化为直立圆柱进行模拟，给出了雷诺数、KC数与拖曳力系数之间的函数关系，其研究结果推进了对植物消波作用内在机理的理解。尽管关于波浪-植物相互作用前人已取得了较多的研究成果，但仍存在发展的空间，如在研究中未考虑植物导致的湍流效应、有限差分法的离散方程难以保证质量守恒、并且对不规则区域适用性较差、GPUSPH方法对于计算资源要求较高、还有待优化等一系列问题。

作为目前最大的开源软件包之一，OpenFOAM具有灵活度高、适用性强，以及能够兼顾计算效率与准确度等优点，已广泛应用于计算流体力学领域，但在植物消波的研究中，对OpenFOAM的应用还有待扩展。本研究基于OpenFOAM下属求解器interFOAM，在非结构化网格中建立三维模型，分别对溃坝、孤立波爬坡以及孤立波在植物水槽中的传播进行模拟，并与物理模型实验结果进行对比，验证了模型的有效性与准确性，进而利用数值模型模拟正弦波在含植物斜坡上的爬升，分析不同波高、不同植物密度、植物拖曳力系数及海滩地形等因素对波浪传播及爬坡过程的影响。

1 数值模型

OpenFOAM采用有限体积法求解流体运动方程，其中interFOAM求解器主要用于求解不可压缩两相流，并采用流体体积函数(VOF)方法[15]确定两相流界面。本文选择西班牙坎塔布里亚大学HIGUERA et al[16-20]开发的速度造波边界与主动消波边界实现数值水槽中波浪的生成与消减，通过添加植物拖曳力项来模拟有植物作用下的波浪运动过程。

1.1 基本控制方程

考虑植物作用的N-S方程如下

(1)

(2)

(3)

(4)

(5)

式中：a为概化植物的直径，N为单位面积内植株数，CD为植物拖曳力系数，CM为植物惯性力系数。

1.2 湍流模型

使用k-ε湍流模型对N-S方程进行封闭，并在湍流模型中考虑植物阻力作用影响，其具体形式为[21]

(6)

(7)

式中：k为湍流动能系数；ε为湍流耗散率；kw与εw分别为由植物引起的附加湍流动能项与湍流耗散项；μt为涡度粘性系数；μ为动力粘性系数；Cε1、Cε2与σk均为封闭系数；Ckp与Cεp为经验常数，分别取1与3.5。

1.3 网格生成

网格是数值模拟中的重要影响因素，决定模拟结果的精确程度。OpenFOAM中包含两个网格生成工具(blockMesh和snappyHexMesh)，其中blockMesh工具可生成结构化网格，snappyHexMesh工具通过读取任意形状几何体的表面特征，并将这些表面特征体现在由blockMesh工具生成的网格中，形成非结构化网格。考虑到采用非结构化网格可以更好地贴合复杂的岸线和地形，我们对重点区域进行局部加密处理，加密后的网格在不降低研究区域计算精度的同时，可以大大提高计算效率。

2 模型验证与分析

2.1 溃坝水流数值模拟

首先检验了OpenFOAM中VOF模型捕捉水-气界面的能力。BRUFAU et al[22]在比利时布鲁塞尔大学水动力实验室进行了溃坝波越过三角障碍物的实验研究，物理模型的几何简化如图1所示。实验装置总长为38 m，左侧为长15.5 m、高0.75 m的蓄水池，距水池右侧10 m处放置一个长6 m、高0.4 m的三角形障碍物，装置右侧是水流出口，在G4、G10、G11、G13和G20五个测点放置压力水位计。

根据上述物理模型实验的设置，本算例数值模型的计算域确定为38 m×0.8 m，采用非结构化网格离散计算域，以实现对三角障碍物边界的贴合，最终网格总数为108 600个。数模左侧与底面均设定为墙边界，粗糙系数ks分别为0.000 2与 0.000 3；右侧与上方为开边界，水流与空气均可自由流通；时间步长选取自动时间步长格式，模拟时长为100 s。

图1 溃坝实验示意图[22]Fig.1 Sketch of the dam break experiment[22]

图2为模拟水深和物理实验实测值(以下简称实测值)的对比，结果显示5个测点的模拟结果与实测值吻合较好。其中G4、G10、G11测点位于水库下游和三角形障碍物上游，溃坝后这些测点的水深先急剧增加而后下降，并呈现多次的波面震荡现象。水波震荡的主要原因在于三角形障碍物对溃坝水波的反射影响，形成了向上游传播的涌波。G10测点模拟的水深较实测值略微偏高，其原因在于该测点位于三角形障碍物起点处，第一次水波到达G10测点时，与障碍物碰撞，会产生压力水位计难以捕捉到的水体飞溅等自由水面激烈变形现象，而数值模拟采用VOF法可精准捕捉到水面的剧烈变化。G13测点位于障碍物顶点处，该点存在明显的水流间断，模型很好地捕捉到了该点水深的变化。G20测点同样存在模拟水深值较实测水深值偏高的现象。考虑到该测点所处地形较G10测点更平缓，模拟水深值偏高的原因可能是网格分辨率不足导致水面捕捉不精准，因此重新构建了分辨率更高的网格进行计算。在分辨率更高的网格下得到的模拟结果与原模拟结果相比并无明显差异，证明原网格分辨率已充分满足计算要求。前人在基于浅水方程分别建立不同模型模拟此算例时，同样存在这一现象[23-25],因此可以认为该点处实测水深可能存在误差。总体而言，OpenFOAM的VOF模型具有良好的捕捉水-气界面的能力，并且斜坡条件下也能精确描述水流的运动过程。

图2 溃坝实验中各测点模拟水深和实测水深的对比Fig.2 Comparison of the simulated value and measured data of water depth in the dam break experiment

2.2 非破碎孤立波爬坡的数值模拟

孤立波是一种非线性波，能够在浅水中传播较远的距离且形状保持不变，其在近岸区域的传播过程中会出现爬坡、波浪破碎等现象。SYNOLAKIS[26]在实验室中进行了一系列经典物理模型实验，测量了非破碎孤立波在斜坡上的传播，物理实验简化示意图如 图3 所示。该实验中水槽右侧有一斜坡海岸，其坡比S=1∶19.85，初始水深h0为1 m，入射波波高Hw为 0.018 5 m。在数值模拟算例中，基于interFOAM模型构建与物理模型相同的计算区域，采用非结构化网格进行离散，计算域网格总数为326 618个，数值探究了非破碎孤立波在斜坡海滩上的爬坡和退水过程。孤立波在计算域左侧由速度造波边界生成，上方为开边界，右侧斜坡及底面为墙边界，水槽底面与斜坡粗糙系数ks=0.000 1，选取自动调整时间步长格式，模拟时长为50 s。为了便于将模拟结果和实验结果进行对比，对模拟结果进行了无量纲化处理，具体公式为

图3 非破碎孤立波实验示意图[26]Fig.3 Sketch of the non-break solitary waves experiment[26]

(8)

(9)

(10)

图4 非破碎孤立波爬坡模拟水位与实测值对比Fig.4 Comparison of the simulated and measured water levels for non-break solitary waves climbing

2.3 孤立波在含植物水槽中传播的数值模拟

植物对波浪具有明显的衰减作用，本算例模拟孤立波在植物水槽中的传播实验，验证在方程中添加植物阻力源项的合理性。MEI et al[27]进行了孤立波在植物水槽中的传播实验，如图5所示。实验室玻璃水槽中初始水深h0设置为0.12 m，孤立波自左向右传播，入射波波高Hw为0.004 8 m，长1.08 m的植物带置于水槽中，植物密度N为1 108 stem/m2，单个植株直径为0.01 m，G3、G4、G5和G6四个波高测点分别位于植物带前端3 m、植物带前端、植物带中部以及植物带后端。

图5 孤立波水槽实验示意图[27]Fig.5 Sketch of solitary waves open channel test[27]

根据实验布置构建了数值水槽，其中左侧为造波边界，右侧采用消波边界，以减少波浪反射的影响；数值水槽底面、侧面设置为光滑墙面；计算域网格总数为254 240个；模拟时间共15 s。植物拖曳力系数CD是一个经验系数，取值会影响概化的植物拖曳力的大小，进而影响植物对波浪消减效应的强弱，是植物消波数值模型里的一个关键参数。此次模拟中，采用波高率定方法[28-29]确定CD的取值，即调整CD使模拟得到的植物前后的波高值与实测值相符，经过率定，CD取1.52。

图6为模型模拟结果与实测值的对比图，表明模拟结果和实测数据基本一致。G3测点靠近造波边界且远离植物带，孤立波通过此处时，还未受到植物的影响，该测点处模拟波高与实测波高吻合良好，说明模型生成的孤立波波形是合理的。G4测点位于植物带前端，其模拟波高较G3测点波高略有升高，主要原因是植物带对孤立波的反射作用。在G5、G6测点处模拟的波高回落，G6测点处测得的模拟波高与G3测点处的波高值相比，衰减率为19.04%，这表明孤立波越过植物带后，波高衰减明显。本算例结果显示，通过波高率定法可得到合理的拖曳力系数，interFOAM 模型在动量方程和湍流方程中加入考虑植物作用的源项可以合理描述植物对孤立波的消减作用。

图6 孤立波越过植物带水位模拟值与实测值的对比Fig.6 Comparison of the simulated and measured water levels of solitary waves over vegetated domain

2.4 植物海岸波浪传播的数值模拟

为了探究植物的消浪效果，THUY et al[30]针对波浪在植物带内的传播及衰减开展了实验室研究，其物理模型设置如图7所示。初始水深为 0.44 m，波高为 0.02 m，周期为20 s的正弦波自左向右传播，实验室海岸分别由坡陡为1∶25、1∶4.7、1∶20.5的3段斜坡以及2段水平地形构成，G1～G6 为波高的测量位置，与左侧边界的距离分别为6.87、7.87、8.87、9.87、10.35 和11.37 m。距离左侧边界 10.36 m 处设置一个宽度为1 m的植物带，为刚性非淹没型植物，单株植物直径为0.005 m，植物密度为2 200 stem/m2。

图7 正弦波爬坡实验示意图[30]Fig.7 Sketch of climbing sine waves experiment[30]

模拟区域的几何构型与物理实验构型相同，并采用非结构化网格进行离散，以精准刻画多段斜坡海岸，网格总数为276 819个；数值模型的左侧、右侧分别为速度造波边界、主动消波边界；底部海岸为墙边界，粗糙系数ks为0.000 2；由波高率定法得到拖曳力系数CD，取1.50；模拟总时长为200 s。图8给出了波高模拟值与实测值的对比，结果显示沿程计算的波高与实测值拟合良好。波浪由左侧向斜坡海岸传播过程中：在植物带前区域，受地形和水深浅化影响长周期正弦波波高逐渐增大；在植物带区域中，波高增长率出现明显的降低趋势。在坡比为 1∶20.5 的斜坡上，植物带前波高增长率为20.53%，经过植物带后，波高增长率为13.13%，增长率降低。模拟结果显示植物会耗散爬坡中波浪的波能，使波高衰减。

图8 正弦波爬坡实验中波高模拟值与实测值对比Fig.8 Comparison of the simulated and measured wave heights in climbing sine waves experiment

图9给出了植物带后端模拟深度平均流速和实测值的对比，模拟值和实测结果吻合良好，波浪速度稳定且呈周期性往复运动特征。

图9 正弦波爬坡实验中植物带后端速度-时间序列模拟值与实测值对比Fig.9 Comparison of the simulated and measured velocity-time series after the vegetation area in climbing sine waves experiment

图10为t=64 s时有无植物工况下植物带附近的流场矢量图。对比发现波浪通过植物带时其波速明显降低，岸滩上植物可以有效地降低长周期波浪在斜坡海滩上的爬升。

图10 正弦波爬坡实验中模拟有无植物工况流场对比(t=64 s)Fig.10 Comparison of simulated flow field with and without vegetation in climbing sine waves experiment (t=64 s)

表1 不同工况下的参数变化Tab.1 Parameter changes in different cases

3 植物消波影响因子敏感性分析

波浪-植物相互作用过程中存在多种影响因素，例如入射波波高、植物密度、拖曳力系数及海岸陡峭程度。为系统地探究各因素对波浪越过植物带后波高衰减的影响，设置了12种工况分4组进行分析，每组中只设置1种变化参数，变化参数包括波高变化、植物密度变化、拖曳力系数变化、海岸陡峭程度变化，具体参数如表1所示，其中坡陡表示海岸的陡峭程度。

图11给出了不同入射波波高条件下数值水槽内波高模拟值的沿程变化。从图中可看出，波浪在由左向右传播过程中，在植物带前，不同入射波波高(工况1、2、3)模拟的波高增长率分别为19.40%、17.91%和 17.41%，3种工况条件下的波高增长率变化不大，即在植物带前区域入射波波高变化对其波高的增长率影响不显著。波浪越过植物带后，波高增长率发生了较大的变化，入射波波高为0.015、0.020和0.025 m的正弦波波高增长率分别为17.83%、11.27%和 6.32%，即在植物带区域，入射波波高越大，其波浪的衰减越显著。

图11 不同入射波波高条件下模拟波高的对比Fig.11 Comparison of the simulated wave heights in condition of different input wave heights

图12给出了不同植物密度工况下模拟波高的对比，由图可以看出，在植物带前，3种工况条件下(工况4、5、6)模拟的波高基本一致，植物带密度对其前端区域的波浪传播影响较小。但在波浪经过植物带后，3种工况下模拟的波高出现了显著的变化，工况6模拟的波高明显低于其他2种工况条件下的模拟值，说明植物密度越小，植物带对波浪的衰减作用也越小。

图12 不同植物密度条件下模拟波高的对比Fig.12 Comparison of the simulated wave heights in condition of different vegetation densities

模型中将植物对波浪的阻碍作用概化为植物拖曳力，设计3组不同的拖曳力系数(工况7、8、9)分析植物拖曳力对波浪传播的影响，具体波高模拟结果对比如图13所示。在植物带前，3种工况下的波高沿程变化基本一致，经过植物区后，波高发生了较大的变化，工况7对应的模拟波高最大，工况8次之，工况9对应的模拟波高最小。模拟结果表明，拖曳力系数变化对植物带前区域的波浪传播没有影响，在植物带区域，其参数变化对波浪衰减作用明显，随着拖曳力系数的增大，波浪的衰减幅度逐渐增加。

最后，数值探究了海岸陡峭程度对波浪传播的影响。设置了3组工况(工况10、11、12)进行模拟，坡陡越大表示海岸越陡峭，需要说明的是，调整坡陡并不影响植物带的位置与宽度，3种工况下植物带的参数都是一致的。图14给出了波高模拟值的沿程变化，从图中可以看出，从8 m处到植物带前端这一区间内，工况10对应的模拟波高值开始明显高于其他2组工况，这是因为工况10条件下海岸较陡峭，雍水效应明显；工况11与工况12对应的模拟波高值基本一致。进入植物带区域后，各工况模拟波高值差异逐渐明显，波高增长率开始出现不同程度的衰减，在植物带末端，工况10、11、12所对应的波高增长率分别为 -28.83%、17.98%、19.35%，波高的衰减程度随着坡陡减小而减小，即在陡坡条件下，波能耗散最严重。

图13 不同拖曳力系数条件下模拟波高的对比Fig.13 Comparison of the simulated wave heights in condition of different drag force coefficients

图14 不同坡陡条件下模拟波高的对比Fig.14 Comparison of the simulated wave heights in condition of different slopes

4 结论

本文基于开源软件包OpenFOAM中的interFOAM求解器建立数值水槽，首先以有三角障碍的溃坝水流、孤立波爬坡、孤立波在植物水槽中的传播等算例验证了模型的有效性，之后模拟了正弦波在长有植物的海岸上的爬坡过程，设置了不同的工况进行对比分析。

研究结果显示：在不同条件下，模拟的波浪在经过植物带之后都会产生衰减现象，波高的衰减程度与入射波波高、植物密度、拖曳力系数以及海岸坡陡均成正比。通过在岸滩改变植物带特征，如改变植物密度、宽度、高度等参数，可有效降低岸滩上波浪的流速，进而提高植物海岸的生态护岸能力。数值模拟的结果可为近岸波浪衰减特性等研究提供参考。

生长在近岸潮间带的植物大多数都具有柔性或弱柔性的特征，受到波浪冲击后会摆动并产生形变，下一步工作将基于OpenFOAM建立全尺度三维模型，开展柔性植物与波浪相互作用的数值模拟研究。