基于SPH方法的明基床开孔沉箱数值模拟研究

任喜峰，孙昭晨，梁书秀

(大连理工大学 海岸及近海工程国家重点实验室，大连 116024)

海洋环境复杂多变，防波堤作为消减、抵御外海传播到近岸的波浪的基础设施，是设计港口、码头必不可少的一部分。传统的防波堤包括斜坡堤、重力式沉箱防波堤等，随着港口需求的增加、建港条件的日益复杂，传统防波堤在经济性和适用性等方面均不能满足新的需求。随着科学工作者对波浪与结构物之间相互作用认识的逐步深入，多种新型防波堤不断出现。例如：开孔沉箱式防波堤、圆形防波堤和梳式防波堤等。

自1961年Jarlan[1]提出开孔沉箱以来，国内外大量采用了这种形式的防波堤来减小波浪的反射及其对建筑物的作用力。对于波浪与开孔沉箱相互作用的研究，目前主要有物模试验和数值模拟两种途径。陈雪峰等[2]利用物理模型试验研究了开孔沉箱在规则波作用下的反射率、总水平力与相对水深、相对消浪室宽度等主要影响因素之间的关系；行天强等[3]利用物理模型试验研究了明基床上开孔沉箱在规则波作用下的反射系数情况；朱大同[4]得出了波浪与开孔沉箱相互作用的反射率的严格封闭解析公式；陈雪峰等[5]采用改进的VOF方法结合k-ε模型建立二维数值波浪水槽，完成了线性规则波和不规则波与单层开孔沉箱相互作用的数值模拟。

但是，由于波浪与开孔沉箱结构相互作用过程中涉及到波浪破碎、湍流等多种复杂的水动力问题，现有的理论分析和实验分析尚不能完全揭示其作用机理。相较于传统的VOF方法，SPH方法是一种纯拉格朗日性质的无网格方法，对于处理大变形、自由表面追踪等复杂非线性问题具有一定的优势。姜峰等[6]基于SPH方法建立了波浪与平底开孔沉箱作用的二维数值模型。任喜峰[7]等采用域粒子着色法(Color Domain Particle，CDP)改进了SPH方法对于计算流体与开孔沉箱等薄壁结构的边界处理方法。对于带有明基床的开孔沉箱，由于明基床的存在导致波浪在与开孔沉箱作用前已经发生了非线性变形，使得波浪与开孔式沉箱的作用更加复杂。因此，有必要对明基床上开孔式沉箱与波浪的相互作用问题进行进一步的研究。

本文基于微可压缩光滑粒子流体动力学方法(WCSPH)，通过域粒子着色(CDP)技术，对方法进行加强。应用加强的CDP-SPH方法，建立了模拟波浪与明基床上开孔沉箱相互作用的数值模型，通过与实验数据对比验证了方法的适用性和可靠性。最后应用建立的数值模型分析了规则波作用下带有明基床的开孔式沉箱的反射率和水平力与消浪室宽度的关系。

1 数值计算方法

1.1 控制方程

基本的控制方程为连续性方程和拉格朗日形式的Navier-Stokes方程

(1)

(2)

式中：ρ为流体密度；u为速度矢量；t为时间；p为压力；f为体积力矢量；为向量微分算子。

SPH方法的基本思想是插值。对于任意的函数A(r)，其在i点的值，可以通过对其周围的N个点插值得到，即

(3)

式中：i,j分别为不同位置的下标；W(rij,h)为核函数；r是位置矢量；rij为从位置i到位置j的相对位置矢量；h为光滑长度；m质点质量。

因此，离散形式的控制方程可以表示为

(4)

(5)

式中：Fij为粒子关系函数，将在章节1.2给出；I为为了降低数值计算中压力波动所引入的数值耗散项；II为为了数值稳定引入的人工粘性项；c0为参考声速；其取值一般为10倍的最大粒子速度，以保证粒子的密度变化率小于1%；c为声速；cij=max(ci,cj)；nij为由粒子i指向粒子j的单位方向矢量；α为人工粘性项系数，由于本文主要研究波浪问题，取值过大容易造成波浪在传播的过程中衰减过大，因此在计算中取值为0.02；πij=(uj-ui)·rji/|rij|2。

对于微可压缩SPH方法，粒子的压力由状态方程给出，即

(6)

式中：γ为常数，对于水一般可取7。

1.2 域粒子着色方法

图1 薄板导致的支持域分割Fig.1 Thin structure divides particle support domain

SPH方法基于插值的概念，一个位置处的函数值通过其周围粒子的插值得到。包含周围粒子的域称为粒子的支持域。在数值计算中，支持域通常取为φ(|rij|

图2 薄板端点附近区域划分示意图Fig.2 Sketch of regions near deck end

当进行流体与结构物相互作用的模拟时，如果结构物为薄壁机构且结构物两侧都存在流体，例如透空式沉箱的前面板，双层板防波堤等，为了避免物理上不存在的结构物两侧粒子的直接相互作用(图1)，需根据选择边界条件的不同要求结构物的厚度 ϑ>2R(虚粒子法)或ϑ

为了避免由于薄板导致的不必要的粒子细分，这里我们采用域粒子着色方法(CDP)，即先对计算域进行分区，然后对于边界区域内的粒子进行着色处理，使粒子携带它能够与之作用的其他区域信息。在计算粒子之间的关系时检查两个粒子来自的域，并判断两个粒子是否能够发生相互作用,因此，公式(4)和公式(5)中的函数Fij可以取为如下的形式

表1 薄板端点附近分区相互作用关系表Tab.1 Domain relationship of a singular deck end

注：表中带有’的编号表示区域关于结构对称的镜像粒子域。

(7)

图2给出了薄板一端的域划分示意图，表1给出了区域之间的关系表，表中X′为与X关于薄板对称的固定镜像粒子[8]域。通过这样的计算域划分和域关系设定，能够避免区域1和区域3的直接相互作用。需要指出的是，这样的设定并不唯一，对于区域复杂的情况，可以根据需要合理的划分区域和设定相互作用关系表。

1.3 边界条件

文中除特别说明外，所有的固壁边界均设置为自由滑移边界条件。边界处理采用固定镜像粒子法[8]。

2 数值水槽的建立

2.1 主动吸收式波浪数值水槽的建立与验证

数值水槽的长度直接影响计算量的大小，为了减少计算量并消除二次反射，数值水槽的造波边界设置为主动吸收式推板造波边界，根据高睿等[9]的研究，造波板的水平运动速度可以写成

(8)

式中：X0为造波板的位置；U为造波板的水平速度；η0和η分别为造波板处的理论波面和实际波面；ω为波的圆频率；φ为传递函数，使用一阶理论，它的形式可以通过下式给出

(9)

式中：k为波数；d为水深。

图3 水槽中波腹点(x=10 m)和波节点(x=9.230 m) 的波面时间序列(H=10 cm)Fig.3 Time series of free surface at node(x=10 m) and antinode (x=9.230 m) in the wave tank (H=10 cm)

为了验证数值水槽的造波和消波性能，取计算域为平底水槽，水槽长度L=10 m，水深d=0.5 m，波高H=10 cm，周期T=1.6 s。设置水槽的右壁面为光滑直立壁面，波浪在右侧壁面产生完全反射，入射波和反射波叠加后会在水槽内形成驻波。图给出了驻波波腹点和波节点历时曲线，结果显示在腹点处，波浪的振幅能够达到2倍的入射波波高，波浪节点处水面的振幅很小。但是由于波高的增大，波浪的非线性导致波谷变缓，波峰变陡。

2.2 明基床上透空沉箱参数设置

数值水槽长11.2 m，水深d=0.40 m，水槽末端为基床，高度hm=0.10 m，坡度1:2，透空式沉箱座于基床上，沉箱距离造波板9.10 m，沉箱距离基床顶面前缘0.25 cm，沉箱高度0.60 cm, 消浪室宽度Bc=0.30 m，如图4所示。透空式沉箱的几何尺寸和压力测点的位置见图5，沉箱开孔率ε=0.3。

图4 二维明基床透空式沉箱数值水槽布置图Fig.4 Sketch of perforated caisson sitting on rubble mound foundation in two-dimensional numerical flume图5 开孔沉箱几何尺寸和压力测点布置图Fig.5 Geometry of perforated caisson and layout of pressure sensors

3 数值计算结果及分析

3.1 测点压力过程分析

透空式沉箱消浪室内外的水体仅能通过开孔进行交换，沉箱内外的水压差导致水体高速通过开孔，进而导致开孔附近产生较强的湍混，水流结构复杂。这里仅以周期T=1.2 s, 波高H=6 cm的波浪与开孔沉箱的作用过程进行分析。图6给出了1#、3#、4#、7#和12#测点的数值模拟结果与实验结果的比较。从图中可以看出，相较于VOF计算结果，采用本文SPH模型的数值计算结果与实验结果更加吻合，尤其是对于静水面以上的压力测点(4#、12#)，SPH的计算结果与实验值吻合非常好，较VOF计算结果优势明显。

图6 1#、3#、4#、7#和12#测点的数值模拟结果于实验结果的比较Fig.6 Comparison of simulated results and experimental results at sensor 1#, 3#, 4#, 7# and 12#

图7给出了t=13.8 s时开孔消浪室附近自由面位置对比。从图中可以看出，SPH模拟结果的自由面较不规整；在消浪室后墙处，自由位置高出VOF模拟结果较多，这与12#测点SPH方法的点压力计算结果大于VOF的点压力计算结果相吻合。这也说明了，SPH方法对于处理流体与开孔沉箱的作用过程中发生的自由面大变形问题相较于传统的VOF具有一定的优势。

3.2 反射系数的影响因素分析

反射系数是评价开孔沉箱性能的重要指标之一，通过本文建立的SPH模型，对波高分别为0.06 m和0.08 m，周期为1.0 s、1.2 s、1.4 s和1.6 s, 基床高度为0.1 m的工况组合进行数值模拟，提取1#～5#位置处的波面历时数据，采用Goda等[10]的两点法分离输入波和反射波并计算反射系数。图8给出了T=1.2 s、H=0.06 m时的分离的入射波与反射波能量谱。从图中可以看出，入射波为单色波，频率为fp=0.384 Hz，与造波板生成的入射波一致。反射波的能量主要集中在一倍频和二倍频处，且二倍频的能量更显著。

图9 反射系数与相对消浪室宽度的关系Fig.9 Reflection coefficient, Kr vs. Bc/L图10 反射系数与相对基床高度的关系Fig.10 Reflection coefficient, Kr vs. hm/L

图11 SPH计算结果与经验公式的反射系数比较Fig.11 Comparison of reflection coefficient of results from present SPH method and empirical formula

图9、图10给出了反射系数与消浪室宽度(Bc/L)、相对基床高度(hm/L)的关系。从图中可以看出，透空式沉箱的反射系数随相对消浪室宽度的增加先减小后增大；随相对基床高度的增加先减小后增大。从图10中还可以看出，相对基床高度较大时，大波高的反射系数更大，原因可以归结于基床浅化作用导致波浪变形，大波高的波浪的非线性更强。

行天强等[3]通过实验研究认为透空式沉箱的反射系数与相对消浪室宽度、相对基床高度、相对水深(d/L)以及开孔率(ε)有关，拟合了规则波作用下明基床上开孔沉箱的反射系数计算公式

(10)

图11给出了数值计算结果与上述经验公式的比较，可以看出，SPH方法的数值计算结果与实验公式基本吻合，大多数点在y=(1±10%)x的包络线内。所以，本文所建立的SPH模型计算结果可以为分析开孔沉箱反射率提供参考依据。

4 结论

本文基于光滑粒子流体动力学(SPH)方法，并应用域粒子着色(CDP)法，建立了二维数值水槽。进而对波浪与明基床上开孔沉箱的相互作用进行数值模拟。数值计算结果在各压力测点与实验值吻合良好，相较于VOF计算结果，静水面以上的压力测点优势更加明显。通过对比同一时刻VOF计算结果与SPH结果的自由表面，可以发现SPH方法模拟结果对大自由面变形适应更好。最后进行了多组合工况的数值模拟，分析得出开孔防波堤的反射系数随相对基床高度、相对消浪室宽度的增加呈现先减小后增大的趋势。通过与实验得到的经验公式比对，可知数值计算结果与实验结果吻合良好。综上，本文应用改进的CDP-SPH方法，解决了开孔沉箱与流体相互作用的问题，为开孔沉箱的数值模拟提供了新的手段。