2. 河海大学 海岸灾害及防护教育部重点实验室, 江苏 南京 210098;
3. 湖南省水沙科学与水灾害防治重点实验室, 湖南 长沙 410004
2. Key Laboratory of Coastal Disasters and Defence of Ministry of Education, Hohai University, Nanjing 210098, China;
3. Key Laboratory of Water-Sediment Sciences and Water Disaster Prevention of Hu'nan Province, Changsha 410004, China
我国南海地区珊瑚岛礁众多,战略和经济地位十分重要,目前对其开发的步伐日益加快。此类岛礁工程通常在远离大陆和海况恶劣的远海地区,建设和维护成本大。因此如何在诸如海啸、风暴潮等极端海洋灾害下保障该类工程的安全成为亟待解决的一个问题。海啸内部蕴含巨大的能量,是一种破坏力极强的自然灾害。海啸波由外海向近岸传播的过程中会发生非线性变形、破碎以及爬高等水动力过程。珊瑚礁坪通常紧临海岸并且水深较浅,尤其在礁冠处,一般处于平均潮位。海啸波传播过礁坪后将通过破碎损耗大量的能量,因此珊瑚礁在抵御海啸中发挥着积极作用,了解珊瑚岛礁附近海啸波的传播变形特性可为海岸防灾减灾提供一定的决策依据。
国内外学者们对波浪与珊瑚礁相互作用问题展开了大量的数值模拟研究,早期一般采用简单的浅水波方程,并用辐射应力和线性波理论来描述破碎波[1-2]。随后学者开始借鉴用于砂质海岸的波流耦合模型来研究珊瑚礁附近的水动力学问题[3]。近年来Boussinesq模型因为其计算精度和效率较高所以得到了较多的应用。研究发现运用改进的Boussinesq模型能够很好地模拟规则波[4]、不规则波[5]和孤立波[6-7]在珊瑚礁上的传播变形问题。
模拟波浪与结构物作用的另一种可行的方法是采用基于Navier-Stokes方程的三维数值模型,该类模型能很好的克服Boussinesq模型的无法准确模拟垂向流速分布的缺点,对波浪破碎过程和海岸干湿交界线的模拟亦不需要人为的干预。此方法已被成功应用于平直海岸[8]和防波堤[9]附近波浪传播变形以及波生流的模拟。海啸波的首波近似于孤立波,所以学术界对海啸波通常采用孤立波模拟[10]。但是模拟孤立波与珊瑚礁相互作用的问题需要考虑模拟礁前斜坡到礁后海岸的一大片区域和采用十分精细的网格捕捉波浪破碎,这需要消耗大量的计算资源,因此运用基于Navier-Stokes方程的三维数值模型来解决此类问题尚未见到相关文献报道。
为了克服现有研究的不足,本文在Roeber物理实验[11]的基础上,基于三维不可压缩Navier-Stokes方程,利用OpenFOAM开源代码中的两相流求解器interFoam建立三维数值波浪水槽。OpenFOAM作为计算流体动力学(computational fluid dynamics,CFD)开源程序包,以其易扩展性、并行计算稳定高效和求解方法先进等特点在计算流体力学领域得到了广泛的应用,作者已将其成功应用于波浪与海岸结构物的相互作用问题[12]。本文通过建立三维数值波浪水槽模拟珊瑚礁与孤立波相互作用的水动力过程,通过模拟计算分析孤立波在珊瑚礁上的非线性变浅、破碎以及破碎波的演化的过程。
1 数值模型的建立 1.1 控制方程本文利用OpenFOAM开源程序包,基于三维不可压缩的Navier-Stokes方程求解大涡湍流(large eddy simulation,LES)模型。该模型将流体的流动分为大尺度涡动和小尺度涡动,对大尺度涡动采用Navier-Stokes方程直接求解,并运用亚格子(subgrid scale, SGS)模型来代替小尺度涡动。
根据质量和动量守恒定律对三维不可压缩的Navier-Stokes进行滤波处理可得到控制方程:
$ \frac{{\partial \;{{\mathit{\bar u}}_i}}}{{\partial {x_i}}} = 0,\;\;\;\;i = 1, 2, 3 $ | (1) |
$ \frac{{\partial \rho \;{{\mathit{\bar u}}_i}}}{{\partial t}} + \frac{{\partial (\rho \;{{\mathit{\bar u}}_\mathit{i}}\;{{\mathit{\bar u}}_j})}}{{\partial {x_i}}} = \frac{{\partial \bar P}}{{\partial {x_i}}} + \rho {g_i} + 2\mu \frac{{\partial \;{{\bar S}_{ij}}}}{{\partial {x_j}}}-\frac{{\partial \tau _{ij}^r}}{{\partial {x_j}}} $ | (2) |
式中:ui为i方向的速度,xi为笛卡尔坐标系,ρ为密度,t为时间,μ为动力黏度,τijr为残余应力。
当瞬时压强P为大气压强时,则应变速率Sij由亚格子模型确定,可表示为
$ {{\bar S}_{ij}} = \frac{1}{2}\left( {\frac{{\partial {u_i}}}{{\partial {x_j}}} + \frac{{\partial {u_j}}}{{\partial {x_i}}}} \right) $ | (3) |
采用Yoshizawa等[13]提出的一方程模型求解亚格子应力τijy为
$ \tau _{ij}^r = \frac{2}{3}{k_S}{\delta _{ij}}-2{v_t}({{\bar S}_{ij}}-\frac{1}{3}{{\bar S}_{kk}}{\delta _{ij}}) $ | (4) |
式中:δij为Kronecker-delta函数,vt为亚格子尺度涡粘系数,由式(5)求解:
$ {v_t} = {C_k}\bar \Delta \sqrt {{k_S}} $ | (5) |
其中,kS为亚格子尺度湍动能,由式(6)求解:
$ \frac{{\partial {k_s}}}{{\partial t}} + {{\mathit{\bar u}}_i}\frac{{\partial {k_s}}}{{\partial {x_i}}} = \frac{{{\rm{ }}\partial }}{{\partial {x_i}}}\left( {\frac{\mu }{{{P_r}}}\frac{{\partial {k_s}}}{{\partial {x_i}}}} \right)-\frac{{\tau _{ij}^r}}{\rho }\frac{{\partial \;{{\mathit{\bar u}}_\mathit{j}}}}{{\partial {x_i}}}-\frac{{{C_\varepsilon }k_s^{3/2}}}{\mathit{\Delta }} $ | (6) |
式中:Ck=0.094,Cε=0.916,Pr=0.9。
本文采用修正的流体体积法(volume of fluid,VOF)捕捉自由液面,即用体积分数α表示液体体积在计算网格中所占的比例,当α=1时表示液体,当α=0时表示气体,当0<α<1时表示气液分界面,其输运方程为
$ \frac{{\partial \alpha }}{{\partial t}} + \nabla \cdot(\alpha {u_i}) + \nabla \cdot[\alpha \left( {1-\alpha } \right)u_i^r] = 0 $ | (7) |
式中:▽·为散度,▽·[α(1-α)uir]为人工压缩项,保证了方程的解在界面附近满足有界性,并对外部流场无影响,uit为相对速度。
OpenFOAM使用有限体积法对Navier-Stokes方程进行离散,方程求解是采用PIMPLE算法实现压力和速度的解耦,时间离散采用Euler格式,压力梯度离散采用Gauss线性离散格式,拉普拉斯项离散采用修正后Gauss线性离散格式,具体离散格式与过程参考文献[14]。
1.2 造波和消波设置本模型采用Jacobsen等[15]提出修正速度入口方法进行造波。该方法将边界类型分为干网格、湿网格和临界网格三种,其中干网格的边界条件可表示为
$ \mathit{\boldsymbol{n}}\cdot\nabla p = 0, {u_i} = 0, \;\alpha = 0 $ | (8) |
式中:n为边界网格面的单位法向向量,模型对湿网格进行了修正,湿网格边界条件根据相应的波浪理论计算所得。模型对临界网格进行了湿面修正,以此计算出临界网格的入流速度和体积分数,对临界网格的处理,避免了入口处自由液面的虚拟振荡。
消波方法采用松弛因子法,对消波区网格的计算变量 ϕ(根据u和α确定)进行数值修正:
$ \phi = {\alpha _R}{\phi _{{\rm{compuded}}}} + (1-{\alpha _R}){\phi _{{\rm{tanget}}}} $ | (9) |
式中:αR为松弛因子,在消波区域和非消波区域交界面上αR=1,在消波区域αR∈(0, 1)。
2 数值模型设置及验证Roeber[11]进行了孤立波与不同珊瑚岛礁地形相互作用的一系列物理模型实验。本文拟运用建立的三维数值波浪水槽对其中的三个典型工况开展数值模拟研究,进一步分析珊瑚岛礁附近孤立波的传播变形以及破碎特征。三种实验工况设置如表 1所示,其中工况1和工况2的模型尺度和礁前斜坡坡度均有所不同,工况3与工况2的模型地形类似,但工况3在礁缘处存在一个梯形礁冠,具体的物理模型设置见文献[11]。
![]() |
表 1 物理模型实验工况设置 Tab.1 Experimental setup of the tested scenarios |
本文所建立的数值波浪水槽尺度与Roeber的物理模型实验水槽一致:以孤立波传播方向为x轴,水槽横断面为y轴,沿水深方向为z轴,其中x=0为孤立波初始造波位置, y=0为水槽中心线,z=0为水槽平底床面位置。数值水槽采用结构化网格,y方向和z方向的网格大小均设置为dy=dz=20 mm。在保证计算精度的前提下,为了节省计算资源对x方向网格进行分段设置。以礁形最复杂的工况3为例, 如图 1(a)所示,网格大小从初始造波位置(x=0 m)的dx=110 mm均匀渐变至礁前斜坡坡脚位置(x=25.9 m)的dx=20 mm;为了精确捕捉珊瑚礁上波浪的演化特性,从斜坡坡脚至水槽结束位置(x=83.7 m)的区段网格大小保持为dx=20 mm,其中斜坡与水槽底交界处以及礁冠附近地形变化剧烈处均采用了网格均匀过渡,如图 1(b)、(c)所示。工况1和工况2的网格分区设置同工况3类似,但由于工况1的物理模型尺度较小,各控制点的网格大小略有差异,其中造波起点(x=0 m)设为dx=70 mm,礁前斜坡坡脚(x=17 m)设为dx=20 mm,水槽末端(x=45 m)同样设为dx=20 mm。
![]() |
Download:
|
图 1 数值域的网格和边界条件设置(以工况3为例) Fig. 1 Numerical grids and boundary conditions of the numerical domain(Case 3) |
除了在计算域的入口采用本文1.2节所述的造波和消波边界外,水槽底面、礁前斜坡、礁坪以及计算域出口均采用固壁无滑移边界,计算域顶面设为自由进出流边界。本模型采用80核服务器(AMD Opteron 6128, 2.0 GHz)进行并行计算,工况1、工况2和工况3的计算时长分别设为35、80和80 s,对应的计算所需时间分别约为4、10、10 d。模型求解采用自适应时间步长,确保克朗数Cr<1(Cr=Δt×max(|Δu|)/min(|Δx|),其中min(|Δx|)和max(|Δu|)分别为最小网格尺寸和最大流速)。本文采用Warner等[16]建议的skill数来验证模型的计算精度,skill数定义采用Wilmott的方法[17]:
$ {\rm{skill}} = 1 - \frac{{\sum |{X_{{\rm{model}}}} \;\;- \;\;{X_{{\rm{obs}}}}{|^2}}}{{\sum {{(|{X_{{\rm{model}}}}\;\; -\;\; {{\bar X}_{{\rm{obs}}}}\left| + \right|{X_{{\rm{obs}}}} \;\;-\;\; {{\bar X}_{{\rm{obs}}}}|)}^2}}} $ | (10) |
式中:Xmodel和Xobs分别表示模型计算值和实验或现场观测值,上画线表示取平均值,skill数越接近1,模型精度越高,越接近0,精度越低。
为了检验数值水槽中的孤立波在传播过程中的稳定性和精确性,本文首先对孤立波在无珊瑚礁模型的平底水槽中的传播进行验证。按工况3计算域和网格设置,设孤立波波高H=0.75 m,静水深h=2.5 m,则H/h=0.3,该波浪具有较强的非线性。图 2所示为水槽中不同位置点的数值计算的波高时间序列,同时对比了数值模拟的波形与理论波形
![]() |
Download:
|
图 2 孤立波的传播验证 Fig. 2 Validation of the solitary wave transformation |
以工况1的计算结果和实验测量为例,图 3(a)选取了6个具有代表性的时刻对比了自由液面的沿礁分布,其中自由液面和时间坐标均进行了无量纲化。当t/(h/g)0.5=55.1时,孤立波因与礁前斜坡相互作用而变陡;随后在t/(h/g)0.5=57.3以激破波的形式在礁缘附近发生破碎;t/(h/g)0.5=59.3破碎波崩塌在初始干燥的礁坪上,并以涌波的形式迅速向前传播,波前线即为干湿交界线;当t/(h/g)0.5=72.3,由于部分水体在礁缘附近反生回退,礁缘开始露出水面,并在此处形成自由跌水。对于前四个时刻,skill值均接近于1,说明数值模拟的结果和实验测量符合良好;后两个时刻skill值稍低,这可能是因为破碎波卷入空气后,仪器测量精度受影响造成。
![]() |
Download:
|
图 3 工况1不同时刻自由液面沿礁的空间分布 Fig. 3 The cross-reef free surface elevations at different time steps for Case 1 |
图 3(b)同样以工况1为例,选取了水槽中6个具有代表性的位置计算的自由液面时间序列与实验结果进行比较。数值模型合理地模拟孤立波在远海测(x=10 m)的传播、礁前斜坡上的浅化变陡(x=19 m, x=21 m), 各位置的skill数均接近于1,仅由于破碎波的影响在礁缘破碎点附近(x=23 m)以及近礁坪中央处(x=32 m)吻合程度稍低,证明所采用的模型能较好地适应本例中非淹没礁坪的情况。相对于传统的基于水深积分的Boussinesq模型,该模型并不需要在程序中对干湿交界线做特殊的处理。
所采用的三维数值模型特点之一是可以预测水平流速的垂向分布,图 3(c)对比了四个测流点测量值与数值模拟值,流速均进行了无量纲化。前三个测流点均位于礁缘附近,因此波浪破碎后流速出现了明显增大,此处水流经历了亚临界流向临界流的转变[11]。无论是波浪破碎点前(x=19 m),破碎点处(x=21 m),还是破碎点后(x=23 m),数值模型均与实测结果符合较好(skill值接近于1),并优于Roeber采用Boussinesq模型计算的流速[11]。在近礁坪中央位置(x=32 m),流速吻合程度同样有所下降,这同样是由于测量仪器受破碎波的影响精度下降,此处测量数据明显比较离散。
图 4(a)展示了工况2中6个不同时刻自由液面沿礁的空间分布。与工况1不同,工况2的礁坪一直处于淹没状态,孤立波在t/(h/g)0.5=62.3时已传播到礁前斜坡上,由于浅化作用开始变陡;大约在t/(h/g)0.5=64.3时孤立波以更为激烈的卷破波的形式发生了破碎;t/(h/g)0.5=65.8时可见卷曲的水舌撞击底床后向前发生飞溅;至t/(h/g)0.5=67.1时, 破碎波开始以滚波的形式在礁坪上继续向海岸方向传播;由于能量耗散,滚波沿程逐渐衰减,如t/(h/g)0.5=76.3时刻所示。同样对于每个时刻,数值计算结果与实验数据吻合较好,skill数均大于0.85,表明该模型具有处理不同形式的激波的能力。
![]() |
Download:
|
图 4 工况1沿礁不同测点自由液面的时间变化 Fig. 4 The time-series of free surface elevation at different cross-reef locations for Case 1 |
图 4(b)展示了工况2中沿礁6个不同位置的自由液面时间序列的对比,图 7同样表明:该模型合理地模拟了孤立波在礁前斜坡上(x=44.3 m)、礁缘附近破碎点前(x=50.4 m)的演化过程;同样由于破碎波的影响,破碎点后(x=57.9 m)以及靠近礁坪中央位置(x=65.2 m),skill数稍低。此外,相对于图 6中的时刻,图 7中所分析的历时曲线延长到了t/(h/g)0.5=170,此时孤立波已经传播到礁后固壁并发生了反射,反射波从礁缘处返回外海时由于色散作用产生的振荡尾波也得到了模型的验证。
图 4(c)对比了实验中5个测流位置的流速历时曲线与模型计算的历时曲线,同样除了接近破碎点位置(x=61.6 m)由于水流和空气作用剧烈测量误差较大外,其他位置的流速历时变化甚至包括反射波造成的尾波振荡均得到了较好的模拟,此外本模型的模拟结果亦优于Roeber采用Boussinesq方程模拟的结果[11],尤其是在礁坪上x=80.2 m处。
工况3与工况2的实验设置相似,仅在礁缘处多出一个非淹没的梯形礁冠,而礁坪仍然处于淹没状态。6个具有代表性位置计算的与实验测量的自由液面空间分布如图 5(a)所示:t/(h/g)0.5=66.5时刻波浪以卷破波的形式在礁冠迎浪测发生破碎,随后在水舌越过冠顶投向礁冠背浪测时(t/(h/g)0.5=69.1),由于流体发生了从亚临界流到临界流的转变,此处产生水跃现象并卷入大量空气;t/(h/g)0.5=72.5时可以观察到水跃激发产生一个反向传播的波浪运动;在t/(h/g)0.5=80.5时发现破碎波以滚波的形式向后传播,并伴随着沿程能量耗散。同样,所有测点skill值均≥0.9,说明数值模拟结果同实验数据吻合良好,表明本文所采用的模型同样适应于礁冠存在时孤立波的传播变形和破碎问题。
![]() |
Download:
|
图 5 工况1沿礁不同测点水平流速的时间变化 Fig. 5 The time-series of streamwise velocity at different cross-reef locations for Case 1 |
图 5(b)展示了工况3中不同位置的时间序列的模拟结果与实验结果的对比,结果表明无论是孤立波的传播变形及其破碎过程,还是水跃造成的反向传播的波浪以及礁后固壁反射造成的尾波现象,都被本模型正确地捕获;除在礁坪上(x=65.2 m)存在一定的破碎波影响外,各位置的skill值均接近于1,再次证明该模型具有较高的精度。
对工况3而言,物理模型实验仅在礁冠的迎浪面布置一个测流点,图 5(c)对比了该点计算流速与实测流速的时间序列:由于该点十分接近破碎点,故流速随时间变化十分剧烈;该序列中出现了由孤立波造成的首峰和由水跃产生的反向波以及固壁反射波分别造成的次峰,而本文所采用的模型较准确地捕捉了整个水动力过程;模型的skill值为0.86,与实验测量较为吻合,并优于Roeber的模拟结果[11]。
4 结论1) 孤立波与礁前斜坡相互作用而变陡,随后在礁缘附近发生破碎,并以滚波的形式在礁坪上向海岸传播,伴随着水体从亚临界流到临界流的转变;不同的礁坪水深和礁冠形态对波浪破碎的类型和破碎波的传播变形规律存在着显著的影响。
2) 本研究建立的数值模型能够合理地模拟不同礁前斜坡坡度、不同礁坪水深以及礁冠存在影响下,孤立波在珊瑚岛礁地形上的浅化、破碎以及破碎波的演化过程。
3) 相较于传统的基于水深的Boussinesq模型,本文所采用的模型可以更合理地模拟水平流速的垂向分布,并且对波浪破碎和干湿交界线的模拟不需要人为的干预。
本研究所建立的基于Navier-Stokes方程三维模型适用于珊瑚岛礁地形周围波浪传播变形的数值模拟,可为进一步认识珊瑚礁对海啸灾害的防御作用提供依据,也可为该模型应用于其他类似海岸地形(沙坝、潜堤等)提供参考。
[1] |
KRAINES S B, YANAGI T, ISOBE M, et al. Wind-wave driven circulation on the coral reef at Bora Bay, Miyako Island[J]. Coral reefs, 1998, 17(2): 133-143. DOI:10.1007/s003380050107 ( ![]() |
[2] |
DOUILLET P, OUILLON S, CORDIER E. A numerical model for fine suspended sediment transport in the southwest lagoon of New Caledonia[J]. Coral reefs, 2001, 20(4): 361-372. DOI:10.1007/s00338-001-0193-6 ( ![]() |
[3] |
LOWE R J, FALTER J L, MONISMITH S G, et al. A numerical study of circulation in a coastal reef-lagoon system[J]. Journal of geophysical research, 2009, 114(C6): C06022. ( ![]() |
[4] |
YAO Yu, HUANG Zhenhua, MONISMITH S G, et al. 1DH Boussinesq modeling of wave transformation over fringing reefs[J]. Ocean engineering, 2012, 47: 30-42. DOI:10.1016/j.oceaneng.2012.03.010 ( ![]() |
[5] |
NWOGU O, DEMIRBILEK Z. Infragravity wave motions and runup over shallow fringing reefs[J]. Journal of waterway, port, coastal, and ocean engineering, 2010, 136(6): 295-305. DOI:10.1061/(ASCE)WW.1943-5460.0000050 ( ![]() |
[6] |
ROEBER V, CHEUNG K F. Boussinesq-type model for energetic breaking waves in fringing reef environments[J]. Coastal engineering, 2012, 70: 1-20. DOI:10.1016/j.coastaleng.2012.06.001 ( ![]() |
[7] |
房克照, 刘忠波, 唐军, 等. 潜礁上孤立波传播的数值模拟[J]. 哈尔滨工程大学学报, 2014, 35(3): 295-300. FANG Kezhao, LIU Zhongbo, TANG Jun, et al. Simulation of solitary wave transformation over reef profile[J]. Journal of Harbin Engineering University, 2014, 35(3): 295-300. ( ![]() |
[8] |
TORRES-FREYERMUTH A, LOSADA I J, LARA J L. Modeling of surf zone processes on a natural beach using Reynolds-Averaged Navier-Stokes equations[J]. Journal of geophysical research, 2007, 112(C9): C09014. ( ![]() |
[9] |
LOSADA I J, LARA J L, CHRISTENSEN E D, et al. Modelling of velocity and turbulence fields around and within low-crested rubble-mound breakwaters[J]. Coastal engineering, 2005, 52(10/11): 887-913. ( ![]() |
[10] |
LIN Pengzhi. Numerical study of solitary wave interaction with rectangular obstacles[J]. Coastal Engineering, 2004, 51(1): 35-51. DOI:10.1016/j.coastaleng.2003.11.005 ( ![]() |
[11] |
ROEBER V. Boussinesq-type model for nearshore wave processes in fringing reef environment[D]. Mãnoa, Honolulu, HI: University of Hawaii, 2010.
( ![]() |
[12] |
JIANG Changbo, YAO Yu, DENG Ya, et al. Numerical investigation of solitary wave interaction with a row of vertical slotted piles[J]. Journal of coastal research, 2015, 31(6): 1502-1511. ( ![]() |
[13] |
YOSHIZAWA A, HORIUTI K. A statistically-derived subgrid-scale kinetic energy model for the large-eddy simulation of turbulent flows[J]. Journal of the physical society of Japan, 1985, 54(8): 2834-2839. DOI:10.1143/JPSJ.54.2834 ( ![]() |
[14] |
GREENSHIELDS C. OpenFOAM user guide, Version 3.0.1[M]. London, UK: OpenFOAM Foundation Ltd, 2015.
( ![]() |
[15] |
JACOBSEN N G, FUHRMAN D R, FREDSE J. A wave generation toolbox for the open-source CFD library:OpenFoam?[J]. International journal for numerical methods in fluids, 2012, 70(9): 1073-1088. DOI:10.1002/fld.v70.9 ( ![]() |
[16] |
WARNER J C, GEYER W R, LERCZAK J A. Numerical modeling of an estuary:A comprehensive skill assessment[J]. Journal of geophysical research, 2005, 110(C5): C05001. ( ![]() |
[17] |
WILLMOTT C J. On the validation of models[J]. Physical geography, 1981, 2(2): 184-194. ( ![]() |