缓坡方程的有限分析数值模式*

黄武平,张庭荣

(广东省水利水电科学研究院,广东省水动力学应用研究重点实验室，广东 广州 510635)



缓坡方程的有限分析数值模式*

黄武平,张庭荣

(广东省水利水电科学研究院,广东省水动力学应用研究重点实验室，广东 广州 510635)

基于时间关联的双曲型缓坡方程,采用九点有限分析法,建立近岸波浪传播变形的数值模式,在对时间关联的双曲型缓坡方程数值离散时,空间导数采用九点有限分析法,时间导数仍采用有限差分法,再采用质量集中方法得到本文的数值模式。为了验证本模式的数值精度和有效性,分别对Berkhoff(1972)实验和突堤地形进行数值模拟,结果表明本模式的模拟值与物理模型试验值、解析解及有限差分模式的模拟值吻合良好,表明本模式可以对近岸波浪传播过程中折射、绕射、反射、浅水变形、波浪破碎和弱非线性等因素的影响进行较好的预测。

缓坡方程;时间关联;有限分析法;质量集中法;数值模拟

有限分析法(FAM)是美国IOWA大学陈景仁教授于20世纪80年代初提出的一种数值计算方法[1-3]。FAM是有限元基础上的一种改进方法,与有限元方法和有限差分方法相似,FAM也需要把计算区域划分为单元子域。他的基本思想是：在离散单元内求控制方程的局部分析解,如方程为非线性的,求分析解时,在子域内将其线性化。用有限分析法解出的数值解是九点格式的,具有精度高,明显的自动迎风性质,计算稳定性好等优点。

1 数值求解

1.1 有限分析法求解

所有的数值计算方法,都是把连续解变成离散解,从而把偏微分方程变成代数方程,利用电子计算机求解。各种方法的差异在于如何建立所求点(中心点)和周围点之间的关系。二维定常纳维尔-斯托克斯方程为：

LU=Re·(UUx+VUy)-(Uxx+Uyy)=G

(1)

其中L表示解算子,G表示重力或压力,表示雷诺数,U、V分别表示沿x和y方向的流速值。

有限分析法基本思想是：在离散单元内的解U不再用插值函数式来表达,而是方程局部线性化后的分析解。以方程式(1)为例,在单元内局部线性化后的方程写为：

(2)

图1 有限分析法

1.2 缓坡方程有限分析数值模式的建立

对于时间关联型的双曲型缓坡方程,把时间变量看成常量时,对于空间变量来说,时间关联型的双曲型缓坡方程变成泊松方程形式,再运用有限分析法对缓坡方程求数值解,得到本文的数值模式。

Berkhoff (1972)[4]在线性理论中引入一个表示地形缓变的小参数,并利用摄动法和沿水深积分的方法导出了一个二阶偏微分方程式,称之为缓坡方程。

▽·(CCg▽Ψ)+k2CCgΨ=0

(3)

其中 速度势函数φ(x,y,t)和Ψ(x,y)的关系为：φ(x,y,t)=Re{Ψ(x,y)e-iωt},C和Cg分别为行进波的局部相速度和局部群速度。

(4)

▽2ζ=-ζ

(5)

如果网格尺寸均匀,即x,y两个方向的网格各自均匀分布,运用上述有限分析法求解可以得出方程(5)的解：

(6)

为了得到显式数值格式,对特解部分采用质量集中法进行简化,得到如下形式：

(7)

将物理量(ζ,k,C,Cg和W)布置在网格节点上,对空间步长导数Δx,Δy和时间步长导数Δt进行离散,空间导数离散采用有限分析法,时间导数离散仍采用有限差分,方程(7)可离散为:

(8)

2 数值模式的检验与比较

2.1 Berkhoff实验

Berkhoff et al在斜坡和椭圆型浅滩组合地形上进行了波浪传播的物理模型试验,试验地形如图2图3所示,浅滩的中心位于(11.0 m,10.0 m),入射波周期T为1.0 s,波高H0是0.05 m,1#～8#断面分别位于x=12,14,16,18,20和y=8,10,12 m。这是个经典算例,用以检验模式对波浪折射、绕射和非线性的有效性。

模式的计算参数如下：计算域为22 m×20 m,左右边界为全反射边界,下游边界为开边界。模型的空间步长Δx=Δy=0.1 m,时间步长Δt=T/40,计算时间为60T。由于计算域相对较小,所以不考虑海底摩擦的影响,即fb=0。

图2 Berkhoff地形水深等值线

图3 试验地形剖面示意

图4是数值模拟的A、B、C、D、E、F六点波面随时间变化的过程线；图5是数值模拟的断面在t=60T时刻的瞬时波面图；图6是1#～8#断面相对波高的数值解与物理模型实验值的比较；图4-11是数值模拟的波面平面分布图。图4和5表明,模式的数值结果稳定、收敛；图6表明,两种模式的数值解与物理模型实验值吻合较好,非线性扩展模型的数值解比线性扩展模型数值解更接近实验值,尤其在浅滩后非线性作用较强区(见图6中的3#～8#),而且九点有限分析格式模式比五点有限差分格式模式模拟的效果略好(见6中的6#)；图7表明,在椭圆型浅滩后,由于波浪的折射作用形成一了个波能辐聚的主峰,同时由于绕射作用使波能沿波峰线侧向传递,这样就在椭圆型浅滩后形成了一系列的波能辐聚和辐散区域。表明本模式对波浪的折射、绕射和非线性有较高的有效性。

图4 数值模拟的 A、B、C、D、E、F六点波面随时间变化的过程线

图5 数值模拟7#断面在t=60T时刻的瞬时波面示意

图6 九点有限分析法和五点有限差分法相对波高比的数值解和物理模型实验值的比较

图7 数值模拟的波面平面分布示意

2.2 突堤附近波浪场

为了检验本模式对斜向波的折射、绕射、反射、浅水变形、波浪破碎和非线性的实用性,本文对突堤附近波浪场进行了数值模拟,将计算结果与关孟儒模型试验的结果进行比较[6]。物理模型实验区域为6 m×10 m,地形为1∶50的斜坡,突堤为从静止水时水边线起长4 m的薄板形状,周期为1.2 s、波高为2 cm的入射波以30°角度从左边边界入射。计算区域如图8所示,计算域宽10 m,长5.75 m(右边界的水深为0.5 cm,其附近波浪已破碎),上下边界及防波堤两侧均为全反射边界。空间步长Δx=Δy=0.025 m,时间步长Δt=0.0012 s,底摩擦因数fb=0.01。

图9是模式数值模拟的波高等值线图,图10是模式数值模拟的波面平面分布图,图11是不同断面沿岸滩方向的波高值与实验值的比较。从图9～图10可知,在斜向波影响下突堤前面出现干涉区,突堤后面出现绕射区,前进波的波高随水深的变浅逐渐变大,最后由于波浪破碎,波能衰减,岸滩处波高减小,图9中实心圆点为物理模型记录的波浪破碎点,这与波高等值线相吻合；从图11可知,非线性扩展模型比线性扩展模型模拟的效果更好,计算值更接近实验值,突堤前反射区数值模拟计算波高值几乎达到入射波高的2倍,比实验值要偏高,尤其在突堤的顶端处更显著,这是由于数值计算时突堤顶端发生全反射,突堤后面的绕射区的数值模拟计算值跟实验值比较吻合,离突堤比较近的断面(y=4.8 m)处绕射比较弱,而离突堤较远的断面(y=3 m)处绕射比较强。表明本模型对斜向波的折射、绕射、反射、浅水变形、波浪破碎和弱非线性效应有比较好的实用性。

图8 计算区域示意

图9 本模式数值模拟的波高等值线示意

图10 本模式数值模拟的波面平面分布示意

图11 不同断面波高分布示意

3 结语

本文以赵红军提出的时间关联的双曲型扩展缓坡方程为控制方程,采用九点有限分析法和质量集中法建立了近岸波浪场的有限分析数值模式。运用本模式数值模拟近岸波浪场的波浪传播,通过对Berkhoff地形实验和突堤进行数值模拟检验了本模式对波浪的折射、绕射、反射、浅水变形、波浪破碎和非线性有很好的实用性。

[1] Chen C J，Chen H C. Finite analytic numerical method for unsteady two-dimendional navier-stokes equations[J]. Journal of Computational Physis，1984(53)：209-226.

[2] Chen C J. Finite analytic numerical solution of heat transfer in two-dimensional cavity flow [J]. Numerical Heat Transfer，1981(4)：179-197.

[3] 陈景仁. 流体力学及传热学[M].北京：国防工业出版社，1984.

[4] Berkhoff J. C. W.． Computation of combined refraction-diffraction[C]∥Proceedings of the 13th International Conference on Coastal Engineering，ASCE 1972：745-747．

[5] 赵红军，宋志尧，徐福敏，等. 一种扩展型缓坡方程的时域计算模式[J].水动力学研究与进展A辑，2009,24(4)：503-511.

[6] 关孟儒，牛克源，汤乃铭，等. 海洋动力学导论[M]. 南京：河海大学出版社，1995.

(本文责任编辑 王瑞兰)

Finite Analytical Numerical Mode for Mild Slope Equation

HUANG Wuping, ZHANG Tingrong

(Guangdong research institute of water resources and hydropower, Guanghou 510635, China)

A numerical model for wave propagation is proposed, which is accomplished by nine point finite-analytic methods based on the time-dependent hyperbolic mild-slope wave equation. The space’s differential is introduced by nine point finite-analytic method and the time’s differential is introduced by finite-difference method when the time-dependent hyperbolic mild-slope wave equation is dispersed, and the numerical model has been got by the lumped mass method. The Berkhoff(1972) experiment and a breakwater stand out are simulated for validating the precision and validity of my numerical model. The calculated results are in agreement with the experimental value, theoretical and finite-difference model’s solution, that indicate that my model is capable of giving good predictions of wave refraction, diffraction, shoaling, breaking energy dissipation and weakly nonlinearity in the near shore zone.

mild-slope equation, time-dependent, finite-analytic method, the lumped mass method, numerical simulations

2016-07-08;

2016-07-15

广东省水利科技创新项目(2012-06)。

黄武平(1985),男,硕士,工程师,从事水动力学研究。

TV139.2

A

1008-0112(2016)08-0001-06