海洋波浪能是储量丰富的可再生清洁能源. 波浪能发电装置的建造和维护以及结构的可靠性问题导致该类发电装置的建设成本过高，是阻碍海洋波浪能利用的重要原因之一. 对此学者设计了护岸结构与振荡水柱式（oscillating water column, OWC）波浪能发电装置相结合的结构^[1-3]. 桩结构相对其他护岸建筑结构具有稳定性好、方便预制、维护成本低和生态影响小等优点^[4-5]. 因此该新型的波能发电装置具有重要的经济价值和现实意义.

近年来，已有学者对桩柱体与OWC装置相结合的结构进行研究. 对于单桩结构，Deng等^[6]通过理论分析研究底部的桩结构、OWC尺寸和入射波角对模型的波能利用率影响；Deng等^[7]进一步通过理论分析研究OWC尺寸、底部桩结构的开口角度、桩结构长度以及入射波角度对波能利用的影响；Xu等^[8]通过实验结合理论分析研究不同水深和波况对模型的波能利用率影响；Xu等^[9]通过建立三维数值波浪水槽研究模型内的液面和压力变化. 关于群桩的OWC结构的研究，文献中仅有Xu等^[4]通过物理模型实验研究桩与OWC装置结合的单排结构的波能利用率问题，结果表明：单排桩结合OWC的结构相对于单桩结合OWC的结构而言单宽波能利用率更高，并且相对无OWC结合的单排桩在相同波况和结构尺寸下波浪的反射系数和透射系数更小.

上述对单排桩与OWC装置结合结构的研究尚未涉及结构周围流场特性的问题，现有文献中对柱体周围流场的分析研究表明马蹄涡的存在是桩柱基础部分冲刷的主要动力因素^[10-11]，因此分析单排桩式OWC装置周围流场特性对在复杂的波流条件下维护装置结构的稳定性具有重要的现实意义.

本文采用易扩展、迭代求解稳定和能进行大规模并行计算的OpenFOAM开源程序包，基于雷诺时均Navier-Stokes方程（RANS）和k− $ \omega $湍流模型，建立数值波浪水槽模拟规则波与单排桩柱式OWC相互作用，并用物理模型实验验证数值模型，根据模拟结果分析结构周围的流场.

利用OpenFOAM开源程序包建立模型，模型以三维不可压缩的RANS作为气液两相的控制方程. 引入流体体积分数 ${s_{\rm{w}}}$，其值代表液体体积在计算网格中所占的比例， ${s_{\rm{w}}} = 0$表示气相， ${s_{\rm{w}}} = 1$表示液相，在气液混合层上 $0 < {s_{\rm{w}}} < 1$. 控制方程如下. 连续性方程：

(1) $\nabla {\boldsymbol{u}} = {\rm{0}}{\rm{.}}$

动量方程：

(2) $\frac{{\partial \rho {\boldsymbol{u}}}}{{\partial t}} + \nabla (\rho {\boldsymbol{u}}{{\boldsymbol{u}}^{\rm{T}}}) = \rho g - \nabla p + \nabla (\mu \nabla {\boldsymbol{u}} + \rho {\boldsymbol{T}}).$

式中： ${\boldsymbol{u}} = [u,v,w]$为笛卡尔坐标系下流速场；ρ为气液混合的密度， $\; \rho = {s_{\rm{w}}}{\;\rho _{\rm{w}}} + (1 - {s_{\rm{w}}}){\rho _{\rm{a}}}$；μ为动力黏度， $\; \mu = {s_{\rm{w}}}{\mu _{\rm{w}}} + $ $ (1 - {s_{\rm{w}}}){\mu _{\rm{a}}}$； $p$为动水压强； $g$为重力加速度； ${\boldsymbol{T}}$为雷诺应力.

为了更加精确的捕捉自由液面，对界面采用修正的体积分数（VOF）方法^[12]捕捉，体积分数输运方程为

(3) $\frac{{\partial {s_{\rm{w}}}}}{{\partial t}} + \nabla ({s_{\rm{w}}}{\boldsymbol{u}}) + \nabla [{{\boldsymbol{u}}_{\rm{r}}}{s_{\rm{w}}}(1 - {s_{\rm{w}}})] = 0.$

式中： $\nabla [{{\boldsymbol{u}}_{\rm{r}}}{s_{\rm{w}}}(1 - {s_{\rm{w}}})]$为界面人工压缩项.

以贴近固态边界附近区域的湍流流场构造为主要研究对象，采用对固态边界附近精度较高、复杂程度较低的k− $ \omega $湍流模型求解 ${\boldsymbol{T}}$：

(4) $ \begin{split} {\dfrac{{\partial \rho {\omega _{\rm{t}}}}}{{\partial t}} +}& {\nabla (\rho {\boldsymbol{u}}{\omega _{\rm{t}}}) = \alpha {p_{\rm{\omega }}} - \beta \rho \omega _{\rm{t}}^2 + }\dfrac{{{\sigma _{\rm{d}}}}}{{{\omega _{\rm{t}}}}}\rho \nabla k {{(\nabla {\omega _{\rm{t}}})}^{\rm{T}}} +\\ &{ \nabla \left[\left(\mu + {\sigma _{\rm{\omega }}}\rho \dfrac{k}{{{\omega _{\rm{t}}}}}\right)\nabla {\omega _{\rm{t}}}\right],} \end{split} $

(5) $\begin{array}{l} \dfrac{{\partial \rho k}}{{\partial t}} + \nabla (\rho {\boldsymbol{u}}k) = {p_{\rm{k}}} - {\beta ^*}\rho {\omega _{\rm{t}}}k + {\rm{ }}\nabla [(\mu + {\sigma ^*}{\mu _{\rm{t}}})\nabla k], \\ \end{array} $

(6) ${p_{\rm{k}}} = {\mu _{\rm{t}}}(\nabla {\boldsymbol{u}}) {(\nabla {\boldsymbol{u}})^{\rm{T}}},{\rm{ }}{p_{\rm{\omega }}} = \dfrac{{{\omega _{\rm{t}}}}}{k}{p_{\rm{k}}}.$

式中： $k$为湍流动能， ${\omega _{\rm{t}}}$为湍流耗散率， ${p_{\rm{k}}} $和 ${p_{\rm{w}}} $分别为湍流动能和湍流耗散率的源项. 为了抑制湍流模型的湍流动能不稳定性，即在符合势流波浪理论的波况下仍在不断计算域内积累湍流动能，导致不正常的波能耗散，湍流动能的源项参照Jacobsen等^[13]的设置，在式（6）中采用流速场的涡量代替应变率. 湍流黏度 ${\;\mu _{\rm{t}}}$定义为

(7) ${\mu _{\rm{t}}} = \rho \frac{k}{{\omega _{\rm{t}}^*}}{p_{\rm{k}}}.$

其中数值湍流耗散率计算式为

(8) $\omega _{\rm{t}}^* = \max \left\{ {\omega _{\rm{t}}},{C_{{\rm{lim}}}}\sqrt {\frac{{2{\boldsymbol{S}}:{\boldsymbol{S}}}}{{{\beta ^*}}}} \right\} .$

雷诺应力 ${\boldsymbol{T}}$可表示为

(9) ${\boldsymbol{T}} = \frac{2}{\rho }{\mu _{\rm{t}}}{\boldsymbol{S}} - \frac{2}{3}k{\boldsymbol{I}},$

式中： ${\boldsymbol{I}}$为单位张量； ${\boldsymbol{S}}$为应变率张量，表达式为

(10) ${\boldsymbol{S}} = \frac{1}{2}[\nabla {\boldsymbol{u}} + (\nabla {{\boldsymbol{u}})^{\rm{T}}}].$

采用Wilcox^[14]建议的 $k - \omega $模型各经验系数：σ_d=13/25 ，β=0.072 ，β^*= 0.09，σ_ω= 0.5 ，σ^*= 3/5，C_lim= 7/8. 模型利用Jacobsen等^[13]开发的Waves2Foam工具库造波，并提供吸波边界. 控制方程中空间离散采用有限体积法，时间离散采用混合Crank Nicolson，运用PIMPLE算法实现方程求解.

物理模型实验在新加坡南洋理工大学水动力实验室的波浪水槽中进行. 如图1（a）和（b）所示，水槽长32.5 m，宽0.54 m，高0.6 m. 在水槽末端铺设有斜率为1∶15的孔隙材料以减小波浪反射. 单排桩式OWC模型按照Froude相似准则采用1∶25的比尺制作模型，单排桩距离造波机18.5 m，由不锈钢材料制作，下段为开口的半圆柱，上段为圆柱，顶部密封，并留有直径为1.4 cm的小孔用以模拟能量转换装置（power take-off, PTO）. 模型的总高度为40 cm，上段圆柱外径为12.5 cm，管壁厚度为3 mm，上段柱底距水槽底24.4 cm. 将4个OWC模型并排组合安装在厚度为3 mm的PVC板上，模拟实际工程中的OWC发电装置的群组布置方式，相邻模型的间隙为0.6 cm. 具体实物模型见图1（c）. 沿水槽布置9个浪高仪（G1~G9）检测模型前后水位的变化，其中G6置于并排的第2个模型管中，用于测量模型管中的液面变化。在该模型气室内布设压电式压力计（P），用于监测气室内压力变化. 实验测试当固定水深0.29 m，波高为3.63 cm时，周期 $T$为0.7、1.5 s的2种规则波工况.

在保证计算精度的同时，为了节省计算资源，将数值水槽的总长度设置为17 m，高度为0.6 m，宽度为0.54 m. 水槽前后4 m分别为造波区和吸波区，中间9 m为计算区. 造波端为波浪入口边界，水槽的两侧、底部、末端以及模型设置为固壁边界. 由于实验中水槽边壁由光滑的玻璃组成，且模型为不锈钢材料，因此数值模拟中所有的固壁边界设置为水力光滑的壁面. 模拟中布置的7个数值浪高仪（G3~G9）与实验对应，排桩式OWC模型布置于计算域正中间位置. 采用非结构化网格进行模拟：将水槽长度方向定义为 $x$，宽度方向定义为 $y$，高度方向定义为 $z$. 在计算区网格的最大尺寸为8.2 mm×10.0 mm×25.0 mm；离静水位10 cm内的网格尺寸细化为8.2 mm×2.5 mm×12.5 mm；在模型附近的网格细化为2.05 mm×2.50 mm×1.50 mm；在相邻的OWC管之间的间隙处的网格细化为1.00 mm×1.25 mm×0.75 mm；为了更精细的刻画模型，在模型表面采用单层贴体的非结构网格，网格尺寸特性为1.00 mm×1.00 mm×0.75 mm，计算网格的总数量约为580万. 模型附近网格和相邻的OWC柱体间隙处的网格如图1（d）、（e）所示. 实测在一个5节点240核（Intel Xeon Scalable处理器）计算集群上完成一组20波周期计算须时约5~8天.

在数值模拟中每个波况至少模拟20个周期，为了验证数值模型能否合理地模拟波浪作用下排桩式OWC模型前后以及桩内的液面变化，取模拟最后5 s的数据进行验证. 数模与实验的液面历时曲线对比如图2所示。图中，η为液面压时。由图可知，排桩前（G5）、桩内（G6）和排桩后（G7）3个位置处的数值计算与实测值均具有很好的一致性. 为了定量地表述数值计算与实测值的一致性，引入正态均方根误差 ${N_{\rm{E}}}$，计算式为

(11) ${N_{\rm{E}}} = \frac{1}{{{o_{\max }} - {o_{\min }}}}\sqrt {\frac{1}{N}\sum\limits_{{{i}} = 1}^{\rm{N}} {{{\left( {{p_{{i}}} - {o_{{i}}}} \right)}^2}} }. $

式中： ${o_{{i}}}$代表第i个实验数据点， ${p_{{i}}}$代表对应的数模数据点， ${o_{\max }} - {o_{\min }}$为实验数据的变化范围（即波高）. 当 ${N_{\rm{E}}}{\rm{ = }}0$时，表明数模100%复现实验观测，当 ${N_{\rm{E}}}{\rm{ = 1}}$时，表明数模与实验的误差量级与实验观测数据变化范围一致. 在图2所示结果中，相应浪高仪测量数据的 ${N_{\rm{E}}}$值分别为1）当 $T{\rm{ = }} 0.7\;{\rm{ s}}$时， ${N_{{\rm{E}}5}}{\rm{ = }} $ $ {\rm{0.043\;1}}$， ${N_{{\rm{E}}6}}{\rm{ = 0}}{\rm{.048\;9}}$， ${N_{{\rm{E}}7}}{\rm{ = 0}}{\rm{.051\;4}}$；2）当T=1.5 s时，N_E5 = 0.054 1，N_E6 = 0.050 5 ，N_E7 = 0.048 5.

针对波浪透反射系数和波能利用率的验算表明，数计算的波浪透反射系数最大误差为8.9%，波能利用率的最大误差为16%，总体来说数值计算与实验实测的吻合良好，说明已建立的数值波浪水槽可以较精确地复现物理模型实验中的各关键物理现象.

根据Perry等^[15]提出的流场切片投影方法（limiting streamline method）将三维流场结构投影在二维平面上，以实现三维流场在二维平面上的可视化. 与二维流动现象的流线不同，由于三维流场存在垂直于二维平面投影的速度分量，其流线在二维平面上的投影可出现相交、渐近的奇异点（singularity）现象.

如图3、4所示分别为在 $T = 0.{\rm{7 }}\;{\rm{s}} $和 $T = {\rm{1}}.5{\rm{ }}\;{\rm{s}} $工况下，图1（b）所示A-A′ 截面处的立面流线特征图. 其中图（b）展示的相位与图（a）中圆形标记对应，图中竖直线为物理模型轮廓线. 分别截取2个工况在单个波浪周期内6个和9个相位进行展示. 为了清楚体现涡结构及其演变，立面流线特征图采用近场小范围绘制.

由图3可知，当 $T = 0.{\rm{7 }}\;{\rm{s}} $的规则波与单排OWC模型相互作用时，以模型扇形内壁前侧流入的部分水体触到模型内壁后向下卷，在靠近底部位置形成一个的马蹄涡. 马蹄涡在垂直方向的大小由涡结构的最高点决定，单个波周期内在1.5~4.0 cm波动. 马蹄涡中心的位置到边壁上的水平距离单个波周期内在1~2.7 cm波动. 在选取的6个相位中，当 $T =$27.55 s时马蹄涡尺寸达到最大. 在图中所示区域，发现单个波周期内马蹄涡以外位置的水平与垂直速度方向保持恒定不变，即在模型前侧靠近底部且在马蹄涡以上的位置水体持续地流向模型内壁.

由图4可知，当 $T = {\rm{1}}.5{\rm{ }}\;{\rm{s}} $的规则波与单排OWC模型相互作用时，也形成了明显的马蹄涡，但比 $T = 0.{\rm{7 }}\;{\rm{s}} $工况下的马蹄涡小得多. 在 $T = 1.5{\rm{ }}\;{\rm{s}} $的工况下，马蹄涡的垂向尺度单个周期内在0.5 ~1.5 cm变化. 马蹄涡中心到模型边壁的水平距离单个周期内仅在0.5~0.7 cm变化. 马蹄涡在2种工况的波周期内始终存在.

值得注意的是，随着 Keulegan–Carpenter（KC）数（KC $ = {U_{\rm{m}}}T/D$，其中 ${U_{\rm{m}}}$为波浪传播过程中的最大流速， $T$为波浪周期， $D$为单个OWC柱体外径）增大，马蹄涡尺寸减小，这暗示着在这2种不同的入射波况下观察到的马蹄涡的成因和尺度可能有明显不同. 根据Sumer等^[16]的研究结果，低KC数下的马蹄涡形成机制存在波生时均流（wave induced steady streaming）的影响，而对于相对较高的KC数，结构物附近流速显著提高，下洗流与底部边界层流动的相互作用，可成为主要的马蹄涡形成机制^[17].

如图5、6所示分别为在 $T = 0.{\rm{7 }}\;{\rm{s}} $和 $T = {\rm{1}}.5{\rm{ }}\;{\rm{s}} $工况下，图1（a）所示距底5 mm的B-B′截面处的平面流线特征图. 其中图（b）展示的相位与图（a）中圆形标记对应，图中弧形轮廓为物理模型轮廓线. 在2个工况下，均可观察到模型柱体狭缝间的收缩射流，在单个稳定的波周期内， $T = 0.{\rm{7 }}\;{\rm{s}} $工况收缩射流的流速最大可达0.167 m/s（背浪方向）和0.125 m/s（迎浪方向）， $T = $ $ {\rm{1}}.5{\rm{ }}\;{\rm{s}} $工况收缩射流的流速最大可达0.57 m/s（背浪方向）和0.534 m/s（迎浪方向）. B截面实际位置仅高于水槽底部5 mm，因此同样可以清晰地展示出马蹄涡结构.

由图5可以看到，在模型前侧形成了覆盖单个OWC柱体宽度方向的凹形马蹄涡，这个马蹄涡围绕着结构内壁，并在模型之间的缝隙产生收缩射流时遭到破坏. 在模型后侧并未形成明显的马蹄涡. 还可以看到在模型缝隙附近形成了2个对称的尾涡，且在模型前后两侧均有尾涡出现，另外形成的尾涡尺寸明显小于模型尺寸. 模型两侧的尾涡形成机理一样：由模型缝隙产生的收缩射流遇到方向相反的流时分离而成. 尾涡会持续存在而不脱离，当遇到波浪导致的反向来流时其构造被破坏. 每个周期内，位于迎浪方向的尾涡存在时长约0.42 s，而位于背浪方向的尾涡存在时长约为0.27 s.

由图6可以看到，在模型缝隙处出现明显的收缩射流，进而在单个波周期内模型的前后两侧均出现较大的尾涡. 该工况下尾涡的形成与最终破坏的机理同 $T = 0.{\rm{7 }}\;{\rm{s}} $的一致，但是尾涡尺寸大很多. 每个周期内，位于迎浪方向的尾涡几乎始终存在，位于背浪方向的尾涡存在时长约为0.83 s. 在模型结构内壁上，出现了与 $T = 0.{\rm{7 }}\;{\rm{s}} $工况类似的马蹄涡，但该马蹄涡在水平方向上的尺寸较小.

为了描述涡的强度 ${S_{\rm{v}} } $，采用Muzzammil等^[18] 提出的计算公式：

(12) ${S_{{{\rm{v}}}} } =\text{π} {U_{\rm{v}} }{D_{\rm{v}} }.$

式中： ${U_{\rm{v}} }$为流场中涡的边界流速， ${D_{\rm{v}} }$为涡大小在长轴和短轴上的平均值.

表1展示了2个工况下的 ${S_{\rm{v}} }$和对应的KC数，以Dargahi^[17]采用的圆柱在定常流中经典马蹄涡强度（π $ UD$）无量纲化 ${S_{\rm{v}} }$. 由表1可知，马蹄涡的强度随着KC数的增大而减小，模型前后两侧的尾涡强度随着KC数的增大而增大. 模型前的尾涡强度（S₁= ${{S_{{\rm{v}} }^{{\rm{l1B}}}} / {\text{π} UD}}$）通常大于模型后的尾涡强度（S₂= ${{S_{{\rm{v}} }^{{\rm{l2B}}}} / {\text{π} UD}}$），特别是在 $T = {\rm{1}}{\rm{.5 }}\;{\rm{s}} $工况时，二者相差 $0.803\text{π} UD$. 对比表中马蹄涡强度（S₃= ${{S_{{\rm{v}} }^{{\rm{hs}}}}/ {\text{π} UD}}$）和尾涡强度可见尾涡强度明显大于马蹄涡强度，其差距可达1~2个数量级. 综合上述涡强度的分析，可以判断，尽管马蹄涡存在，但是其强度过小，因此马蹄涡不是结构周围泥沙输移的主要控制因素. 而较强的尾涡很有可能使模型前的泥沙产生输移并重新分布. 由此可见尾涡对柱体结构周围泥沙的输移有较大影响. 当较强的尾涡、模型缝隙处的收缩射流以及水体持续入流共同作用时，在缝隙前后两侧以及缝隙处有产生严重冲刷的可能.

值得注意的是，模拟结果中马蹄涡强度随KC数变化规律与文献中采用更大KC数的实验观测存在不一致的地方，除了考虑实验工况的区别以外，这也同样提示在小KC数工况下（如 ${\rm{KC}} < 6$）的马蹄涡形成机制与大KC数工况下的机制存在不同.

（1）当规则波与模型相互作用时，模型趾部内壁形成了内凹的马蹄涡流场构造. 在马蹄涡以外的位置观察到水体持续流向模型内壁.

（2）观察到波浪导致的模型缝隙间收缩射流反向时流动分离在模型前后形成的尾涡.

（3）模拟观察到的马蹄涡强度随着KC数的增大而减小，而尾涡强度随着KC数的增大而增大. 2个工况的尾涡强度均显著大于马蹄涡强度.

（4）模拟观察到的马蹄涡强度不足以对OWC桩基的泥沙起动和输移产生显著影响，尾涡流动对OWC桩基泥沙启动的影响更为显著. 由收缩射流导致的高强度尾涡流动将增大桩基附近的床面剪应力，并带来显著的桩基冲刷现象. 工程实践中，宜在桩基处避免薄壁结构和不连续的表面，如在基础处采用完整的圆柱体平台设计方式，以削弱收缩射流导致的潜在尾流涡冲刷问题.

（5）在本文研究基础上，计划研究规则波与紧凑布置的OWC群桩构造OWC波能发电装置相互作用导致的桩基冲刷问题，并进一步研究在非规则波作用下的水动力学和桩基冲刷情况，以获取更贴近工程实践的研究结论.