2 ［9 ］
如下： （ 1）
1 1 h［ ） t］ － h2 （ ·u ·（ h u 珔 珔 t） + 2 6 （ 2）
t
2 2 B1 h2 ［ ） ］ + B2 ［ ） ］+ R ·（ u 珔 ·（ h u 珔 t + g η t + gh η
G = 式中
2
－u · 2 u 珔 珔－
第 23 卷 第 1 期 2012 年 1 月
水 科 学 进 展 ADVANCES IN WATER SCIENCE
Vol． 23 ， No． 1 Jan． ， 2012
DOI： CNKI： 32. 1309. P. 20120104. 2012. 009
应用二阶完全非线性 Boussinesq 方程 模拟破碎波浪
= 0. 5［10］，而 Lynett 取 σ = 10［11］。因此，针对本文采用的控制方程，有必要进行率定，找到合适的参数取值 范围。 1. 3 无反射波浪入射条件 包括 FUNWAVE 在内的诸多基于 Boussinesq 方程的波浪模型均采用内部波浪生成技术 ，来避免波浪二次
［12 ］
反射问题，然而内部造波函数的推导均基于线性化的 Boussinesq 方程。因此，当波浪的非线性较强时， 内部 造波函数则无法产生所要求的波浪 。为克服这一问题，Veeramony 和 Svendsen 采用傅里叶展开方法推导了 Boussinesq 方程的永形波解，并采用特征线方法来实现无反射造波 。这一过程不但涉及到非线性方程组的迭 代解法而且对波浪的分解过程采用了线性长波假定 ，造成了这种方法的不足。这里采用傅里叶方法推导方程 ［13 ］ 的永形波解，结合 Fuhramn 给出的松弛造波技术实现波浪无反射入射的目的 。设方程的解具有以下形式：
G =
） 珔 { d3 [ （ ·u
2
2
－u · 2 u 珔 珔－
1 2 ·u ） （u 珔 珔 10
（ ·u ） 珔 ] } + η d [ 1 3
2
－
1 － u · 2 u 珔 珔－ ·u 珔 t 3 （ 4）
]
η（ 2 h + η） （ ·u 珔 t） 3 1. 2 涡粘方法考虑波浪破碎 本文采用涡粘方法处理波浪破碎问题 Rb = 式中
2 高阶为 Ο （ εμ ） 项 （ ε = A0 / h0 ，μ = h0 / L0 ，A0 、h0 、L0 分别为特征波幅、特征水深和特征波长 ） ，为二阶弱 非线性方程。为提高方程的非线性性能使其具有二阶完全非线性特征， 需在非线性项 G （ 式 3 ） 中保留 2 3 2 ［9 ］ Ο （ μ ） 阶所有非线性项 （ 最高阶为 Ο （ ε μ ） ） 如下 ：
线性模型能更准确地描述波浪的非线性演化过程 ； 在接近波浪破碎点处，具有完全非线性特征的数值模型较 弱非线性模型能更精确地捕捉波浪破碎点位置 。因此，应用具有完全非线性特征的 Boussinesq 类模型研究波 浪的传播和破碎过程具有非常重要的意义 。 而现有的多数 Boussinesq 类波浪破碎模型仅具有弱非线性特征， 如 Madsen 等
1， 2 1 房克照 ，邹志利 （ 1. 大连理工大学海岸和近海工程国家重点实验室，辽宁 大连 116024 ； 2. 河海大学海岸灾害及防护教育部重点实验室，江苏 南京 210098 ）
摘要： 建立了基于高阶 Boussinesq 水波方程的一维波浪破碎数值模型 。基于一组具有二阶完全非线性特征的 Boussinesq 水波方程，建立了交错网格下的高精度差分格式，推导了适用于该组方程的永形波解析解，其和松弛造波技 术相结合实现了数值波浪水槽中 （ 强） 非线性波浪的无反射入射 。 通过模拟封闭容器内水体晃动问题对数值格式 进行了验证，通过模拟孤立波在斜坡海岸上的浅化过程说明了将方程从弱非线性发展到完全非线性的必要性 。 采 用涡粘方法处理波浪破碎，利用物理模型实验数据，分析了模型中各波浪破碎参数对数值结果的影响并对参数进 行了率定。应用该模型对规则波在斜坡海岸上的传播 、 变浅以及破碎过程进行了数值模拟研究，数值结果同实验 数据吻合良好，验证了模型的有效性。 关键词： Boussinesq 水波方程； 非线性； 波浪破碎； 有限差分法； 交错网格 中图分类号： TV139 文献标志码： A 6791 （ 2012 ） 01009608 文章编号： 1001-
1 2 ·u ） （u 珔 珔 10
2 ηh 珔－ （ ·u 珔） ] } － η h ·u 3
t
（ 3）
h 为静水深； d = h + η 为当地水深； g 为重力加速度； = （ x， y ） 为二维偏微分算子； B1 = 29 / u 珔| u 珔| /
885 ，B2 = 2 / 59 为常数。上述方程的色散性为精确色散关系的 Padé ［ 2 ，2］ 阶近似， 变浅作用性能在浅水到 中等水深有较高的精度。方程中 R = R f + R b 为扩展项，用以考虑底摩阻、波浪破碎， 其中， R f = f d 为底摩阻 （ f = 0. 005 为底摩擦系数） ，R b 为波浪破碎引起的耗散项 （ 在后面给出 ） 。 上述方程中非线性最
{
f1 V Theory ， （ 1 － f1 ） V Numerical ，
Ⅰ区 （ 10 ）
（ 1 － f1 ） V Theory + f1 V Numerical ，Ⅱ 区
1. 4
数值求解方程
有限差分方法由于具有直观性和程序实现简单的优点 ，被广泛用来数值求解 Boussinesq 类方程， 发展得 比较完善。本文的数值求解格式同 FUNWAVE 基本一致， 即采用有限差分近似方程中的时间、 空间导数， 采用三阶预报、四阶校正的 ABM 格式进行时间积分。但 FUNWAVE 是基于非交错网格建立的， 而本文在交 错网格上对方程 （ 1 ） 至方程 （ 3 ） 进行数值求解。将控制方程在 C 网格上进行离散， 网格点标记为从 0 到 M，矢量变量如速度定义在网格点上 ，而标量变量如水深和波面则定义在网格中心 ，如图 1 所示。 在上述交错网格系统中，式 （ 1 ） 中标量变量 （ 如水深 d ） 和要求解的变量定义在同一位置， 其导数可 直接采用普通有限差分格式求得 ，而矢量变量 （ 如速度） 的导数则通过交错网格下差分公式求得 ， 相应的， 式 （ 2 ） 中标量变量 （ 如波面、水深） 的导数通过交错网格下差分公式求得 ，而矢量变量 （ 如速度） 的导数 则通过普通有限差分公式求得。文中涉及的非交错网格下差分格式同 FUNWAVE 完全一致， 不再列出， 交 错网格上的差分公式可表示为 （ 以求解动量方程中 i 点处 η 的各阶导数为例） η ( x ) =
近岸水域波浪运动十分复杂，具有明显的非线性特征，如波浪之间、波流之间的非线性相互作用、 破碎 点附近波浪的运动等。因此，建立数学模型时， Boussinesq 水波方程的非线性性能是需要考虑的重要方面 。 一些基于具有完全非线性特征的 Boussinesq 类模型研究结果表明
［14 ］
， 在波浪未破碎区域， 这些模型较弱非
第1 期
房克照，等： 应用二阶完全非线性 Boussinesq 方程模拟破碎波浪
97
1
1. 1
数值模型
Boussinesq 水波方程 邹志利推导了以水深平均流速 u 珔 和波面升高 η 表达的二阶 Boussinesq 方程 ） = 0 η t + ·（ d u 珔 u ·） u 珔 珔 珔+ gη + G = t + （ u ） 珔 { h3 [ （ ·u
［5 ］
提出的水滚模型，李绍武等
［6 ］
和马小舟等
［7 ］
建立的模型等。 同时， 许多 Boussinesq 类波浪破
［1 ］
碎模型，如上面提到的模型以及具有二阶完全非线性特征的波浪破碎模型 FUNWAVE 算失败。相对而言，在交错网格上求解 Boussinesq 方程，能使数值模型更加稳定
［8 ］
图1
数值计算域示意图以及交错网格上变量布置 Fig. 1 Sketch of computation domain and variable location on stagger grid
海绵层 式中，函数 f1 的值在各区域内从 0 缓慢过渡到 1 ， 本文取 为 tanh 函数，理论值 （ 下标 Theory） 由式 （ 7 ） 给出， 而 数值解 （ 下标 Numerical ） 则由模型计算得到， 计算过程 中，I、II 区的长度以及海绵层的长度均取 1 个波长。
{
F2 = － kA1 U1 5 G = hk2 ωA1 U1 － 6
2
(
1 2珋 2 U2 + k 1 2 15
)
F3 = ， G3 = k
(
)
1 2 2 14 珋 k ω A1 U1 + k k － 1 kU2 + ω U2 A1 + 3 5 1 3 － 4珋 7珋 ω A2 U1 k 2 U2 + k 3 5
(
)
（ 9）
[(－ 3 2
)
]
n n 求解上述方程可得各阶解的幅值 A 和 U ， 将其代入式 （ 7 ） 即为控制方程的永形波解， 该过程无需采用迭代方
法。求得上述解析解后，在造波区 I 区、 II 区以及海绵层 ） 的取值分别为 内 （ 见图 1 ） ，变量 V （ 包括 η 和 u 珔 V =
［1 ］
模拟波浪破碎引起的能量耗散， 通过和实验数据对比，
分析了破碎模型中各参数取值对数值模拟结果的影响并对参数进行了率定 。应用该模型模拟了规则波在斜坡 海岸上的传播、变浅以及破碎过程，将波高、增减水等数值模拟结果同实验数据进行了对比分析 ，两者吻合 良好，验证了模型的有效性。