《1 前言》
1 前言
由于近岸区的波浪运动总是表现出较强的非线性, 此时波浪的非线性作用不能忽略, 因此将介绍一种能有效考虑波浪非线性作用的数学模型, 即Boussinesq类方程, 是一种能够描述规则波和不规则波在复杂地形上发生浅化、折射、绕射和反射效应的相当有效的数学模型, 可分为经典 (标准) Boussinesq方程、低阶Boussinesq方程和高阶Boussinesq方程等。各种Boussinesq方程均在一定程度上考虑了波浪的非线性效应和频散效应, 因此能够反映不同频率分量间的能量输运、波形变化、波群演化等。
Boussinesq于1871年提出了著名的Boussinesq理论, 理论中考虑了自由面曲率对浅水长波流动的影响。经过100多年的发展, 该理论取得了长足发展, 已被推广用于深水短波的作用研究, 并且已成为模拟近岸区波浪运动的强有力的工具。第一个能考虑变水深地形并得到较好应用的Boussinesq方程当首推1967年Peregrine
[1 ] 提出的Boussinesq方程, 后人称其为经典或标准Boussinesq方程。在随后的20多年的时间, 许多学者
[2 ,3 ,4 ,5 ,6 ] 应用该方程或其改进形式研究了浅水中的长波、短波以及波流的作用等, 讨论了其适用范围并提出了一些求解该类方程的有效数值格式。由于标准Boussinesq方程所反映的频散关系随水深的增大而误差不断增加, 因此只适用于相对浅水中的波浪问题。为了克服这一限制, 1991年Madsen等
[7 ] 提出了改进线性频散特征的新型Boussinesq方程, 将Boussinesq方程的适用范围推广到可以用于中等水深中的波浪问题, 1992年Madsen等
[8 ] 又在该方程中考虑了缓变地形的影响, 1995年Sch!?ffer和Madsen
[9 ] 又对该方程做了进一步的改进;Nwogu
[10 ] 于1993年给出了另一种形式的Boussinesq方程, 该方程和Madsen提出的方程是弱非线性和弱频散性Boussinesq类方程的代表性方程, 并且被一些学者
[11 ,12 ,13 ,14 ,15 ] 用来研究比较复杂地形上的波浪变形, 取得了比较满意的结果。最近, Boussinesq类方程又有了新的进展, Wei
[16 ] 等导出了二阶完全非线性Boussinesq方程, 而Chen等
[17 ] 导出了四阶含流Boussinesq方程, 并采用该方程研究了流对波的影响;Gobbi和Kirby
[18 ] 则导出了四阶完全非线性的高阶Boussinesq方程, 并采用一种高精度的预报校正格式数值模拟了淹没浅堤上的波浪变形, 取得了相当满意的结果。从理论上说, 高阶Boussinesq方程更精确、更能反映实际情况, 但在实际应用时, 由于方程十分复杂, 加上高阶导数的出现, 使得程序实现起来相对困难, 至今为止, 高阶Boussinesq方程在复杂地形上的应用还仅限于一维问题。基于以上原因, 笔者采用一种改进的Boussinesq方程来研究波浪变形问题, 该方程由Beji和Nadaoka
[19 ] 于1996年首次导出, 当水深为常数时, 该方程和Madsen
[8 ] 导出的Boussinesq方程完全一致, 而对变水深地形, 该方程是Madsen
[8 ] 提出的Boussinesq方程的正式推导和对其中不正确之处的修正。为了说明Beji和Nadaoka
[19 ] 提出的Boussinesq方程的准确程度, 笔者对比较复杂地形上的波浪变形进行了数值模拟, 结果表明该方程得到的波要素是可信的, 可以用于实际当中的波浪问题计算。
《2 Boussinesq方程的正式推导》
2 Boussinesq方程的正式推导
1992年, Madsen
[8 ] 等提出了缓坡地形上的一种Boussinesq方程, 通过理论分析和应用实践表明, 该方程比经典Boussinesq方程
[1 ] 适用的水深范围更广, 但该方程没有严格的数学推导, 于是, Beji和Nadaoka
[19 ] 从经典Boussinesq方程出发, 通过简单的数学代换, 正式推导了Madsen等提出的Boussinesq方程, 并指出了Madsen提出的Boussinesq方程中的不正确之处。下面给出Beji和Nadaoka导出的Boussinesq方程。
Peregrine
[1 ] 导出的渐变水深标准Boussinesq方程为
u t + ( u ⋅ Δ ) u + g Δ η = h 2 Δ [ Δ ⋅ ( h u t ) ] − h 2 6 Δ ( Δ ⋅ u t ) ( 1 ) η t + Δ ⋅ [ ( h + η ) u ] = 0 ( 2 ) u t + ( u ⋅ Δ ) u + g Δ η = h 2 Δ [ Δ ⋅ ( h u t ) ] - h 2 6 Δ ( Δ ⋅ u t ) ( 1 ) η t + Δ ⋅ [ ( h + η ) u ] = 0 ( 2 )
式中, u = (u , v ) 为深度平均的二维速度矢量;η 为自由面位移;h =h (x , y ) 为静水深;g 为重力加速度;下标t表示关于时间的导数;Δ Δ Δ Δ
通过简单的加减处理, 并将二阶非线性频散项采用一阶长波关系式
u t + g Δ η = 0 ‚ η t + Δ ⋅ ( h u ) = 0 ( 3 ) u t + g Δ η = 0 ‚ η t + Δ ⋅ ( h u ) = 0 ( 3 )
u t + ( u ⋅ Δ ) u + g Δ η = ( 1 + β ) h 2 Δ [ Δ ⋅ ( h u t ) ] + β g h 2 Δ [ Δ ⋅ ( h Δ η ) ] − ( 1 + β ) h 2 6 Δ ( Δ ⋅ u t ) − β g h 2 6 Δ ( Δ 2 η ) ( 4 ) u t + ( u ⋅ Δ ) u + g Δ η = ( 1 + β ) h 2 Δ [ Δ ⋅ ( h u t ) ] + β g h 2 Δ [ Δ ⋅ ( h Δ η ) ] - ( 1 + β ) h 2 6 Δ ( Δ ⋅ u t ) - β g h 2 6 Δ ( Δ 2 η ) ( 4 )
式中, β 为一常数, 取β =1/5
[19 ] ;Δ 2 Δ 2 x 和y 方向的动量方程分别为:
∂ u ∂ t + u ∂ u ∂ x + v ∂ u ∂ y + g ∂ η ∂ x = ( 1 + β ) h 2 3 ⋅ ( ∂ 2 u t ∂ x 2 + ∂ 2 v t ∂ x ∂ y ) + β g h 2 3 ( ∂ 3 η ∂ x 3 + ∂ 3 η ∂ x ∂ y 2 ) + ( 1 + β ) h 2 ( 2 h x u x t + h x x u t + h x y v t + h y v x t + h x v y t ) + β g h 2 ( 2 h x η x x + h x x η x + h x y η y + h y η x y + h x η y y ) ( 5 ) ∂ v ∂ t + u ∂ v ∂ x + v ∂ v ∂ y + g ∂ η ∂ y = ( 1 + β ) h 2 3 ⋅ ( ∂ 2 v t ∂ y 2 + ∂ 2 u t ∂ x ∂ y ) + β g h 2 3 ( ∂ 3 η ∂ y 3 + ∂ 3 η ∂ x 2 ∂ y ) + ( 1 + β ) h 2 ( 2 h y v y t + h y y v t + h x y u t + h x u y t + h y u x t ) + β g h 2 ( 2 h y η y y + h y y η y + h x y η x + h x η x y + h y η x x ) ( 6 ) ∂ u ∂ t + u ∂ u ∂ x + v ∂ u ∂ y + g ∂ η ∂ x = ( 1 + β ) h 2 3 ⋅ ( ∂ 2 u t ∂ x 2 + ∂ 2 v t ∂ x ∂ y ) + β g h 2 3 ( ∂ 3 η ∂ x 3 + ∂ 3 η ∂ x ∂ y 2 ) + ( 1 + β ) h 2 ( 2 h x u x t + h x x u t + h x y v t + h y v x t + h x v y t ) + β g h 2 ( 2 h x η x x + h x x η x + h x y η y + h y η x y + h x η y y ) ( 5 ) ∂ v ∂ t + u ∂ v ∂ x + v ∂ v ∂ y + g ∂ η ∂ y = ( 1 + β ) h 2 3 ⋅ ( ∂ 2 v t ∂ y 2 + ∂ 2 u t ∂ x ∂ y ) + β g h 2 3 ( ∂ 3 η ∂ y 3 + ∂ 3 η ∂ x 2 ∂ y ) + ( 1 + β ) h 2 ( 2 h y v y t + h y y v t + h x y u t + h x u y t + h y u x t ) + β g h 2 ( 2 h y η y y + h y y η y + h x y η x + h x η x y + h y η x x ) ( 6 )
上述式中, 下标x 和y 分别表示关于x 和y 方向的偏导数。
∂ η ∂ t + ∂ ( h + η ) u ∂ x + ∂ ( h + η ) v ∂ y = 0 ( 7 ) ∂ η ∂ t + ∂ ( h + η ) u ∂ x + ∂ ( h + η ) v ∂ y = 0 ( 7 )
式 (7) 至式 (9) 构成了二维Boussinesq方程。对一维问题, Boussinesq方程简化为:
∂ u ∂ t + u ∂ u ∂ x + g ∂ η ∂ x = ( 1 + β ) h 2 3 ∂ 2 u t ∂ x 2 + β g h 2 3 ∂ 3 η ∂ x 3 + ( 1 + β ) h 2 ( 2 h x u x t + h x x u t ) + β g h 2 ( 2 h x η x x + h x x η x ) ( 8 ) ∂ η ∂ t + ∂ ( h + η ) u ∂ x = 0 ( 9 ) ∂ u ∂ t + u ∂ u ∂ x + g ∂ η ∂ x = ( 1 + β ) h 2 3 ∂ 2 u t ∂ x 2 + β g h 2 3 ∂ 3 η ∂ x 3 + ( 1 + β ) h 2 ( 2 h x u x t + h x x u t ) + β g h 2 ( 2 h x η x x + h x x η x ) ( 8 ) ∂ η ∂ t + ∂ ( h + η ) u ∂ x = 0 ( 9 )
《3 Boussinesq方程的离散格式》
3 Boussinesq方程的离散格式
Boussinesq方程可采用多种格式进行离散, 离散格式可分为两大类, 一类是单一网格上的离散格式, 另一种是交错网格上的离散格式。每类格式中又包括多种格式, 如高阶和低阶的预报校正格式、ADI格式等。高阶格式的优点是精度高, 但程序设计复杂, 计算效率低, 另外还可能产生较大的数值振荡, 导致数值计算不稳定, 因此高阶格式在整体上并不优越, 尤其对多维问题;在交错网格上对方程进行离散, 可以避免诸如采用中心差分格式时引起的振荡等不利因素, 但是交错网格也存在一些缺陷, 如边界条件很难对各物理量前后一致地给出, 另外采用交错网格得到的各物理量的最终结果由于不是在相同的网格点上得到的, 因此要得到某个网格点上的所有待求物理量, 还需要进行一定的插值处理;ADI格式是一种在流体力学问题中得到广泛应用的格式, 由于对大多数流体力学问题, 采用这种格式对方程进行离散, 可以得到系数矩阵为三对角的代数方程组, 从而可以采用追赶法进行快速求解, 因此一般说来采用ADI格式的计算效率比较高, 并且格式一般为无条件稳定或者临界无条件稳定, ADI格式的不足是中间时刻第一类边界条件的给定问题, 如果根据某个解析表达式直接给出中间时刻的第一类边界条件, 则可能使计算的最终结果的精度较差, 此外对于Boussinesq方程, 采用ADI格式仍然需要迭代, 因此其快速求解的优势并不能很好地体现出来。基于以上这些原因, 采用非交错网格下的预报校正格式进行离散, 然后采用有关的代数方程组的求解方法来得到问题的数值解。下面介绍Boussinesq方程的离散格式, 可分为预报阶段的离散格式和校正阶段的离散格式。
在具体给出离散格式之前, 为了书写方便, 首先定义以下符号:
上述式中, δF 表示函数F 的差分。除了上面定义的符号之外, 对动量方程中的∂3 η /∂x 3 和∂3 η /∂y 3 在边界点附近的离散进行如下处理:
∂ 3 η ∂ x 3 ∣ ∣ ( 2 , j ) = η 4 , j − 3 η 3 , j + 3 η 2 , j − η 1 , j Δ x 3 ‚ ∂ 3 η ∂ x 3 ∣ ∣ ( N I − 1 , j ) = η N I , j − 3 η N I − 1 , j + 3 η N I − 2 ‚ j − η N I − 3 , j Δ x 3 ； ∂ 3 η ∂ y 3 ∣ ∣ ( i , 2 ) = η i , 4 − 3 η i , 3 + 3 η i , 2 − η i , 1 Δ y 3 ‚ ∂ 3 η ∂ y 3 ∣ ∣ ( i , N J − 1 ) = η i , N J − 3 η i , N J − 1 + 3 η i , N J − 2 − η i , N J − 3 Δ y 3 ； ∂ 3 η ∂ x 3 | ( 2 , j ) = η 4 , j - 3 η 3 , j + 3 η 2 , j - η 1 , j Δ x 3 ‚ ∂ 3 η ∂ x 3 | ( Ν Ι - 1 , j ) = η Ν Ι , j - 3 η Ν Ι - 1 , j + 3 η Ν Ι - 2 ‚ j - η Ν Ι - 3 , j Δ x 3 ； ∂ 3 η ∂ y 3 | ( i , 2 ) = η i , 4 - 3 η i , 3 + 3 η i , 2 - η i , 1 Δ y 3 ‚ ∂ 3 η ∂ y 3 | ( i , Ν J - 1 ) = η i , Ν J - 3 η i , Ν J - 1 + 3 η i , Ν J - 2 - η i , Ν J - 3 Δ y 3 ；
上式中的 (2, j ) 和 (N I -1, j ) 分别表示x 方向的第一排和最后一排的域内节点编号, N I 为x 方向的最大节点编号;而 (i , 2) 和 (i , N J -1) 分别表示y 方向的第一排和最后一排的域内节点编号, N J 为y 方向的最大节点编号。此外对动量方程中的∂η /∂x 和∂η /∂y 的离散进行如下特殊处理:
∂ η / ∂ x = δ x η = η i + 1 , j − η i − 1 , j 2 Δ x − Δ x 2 6 δ x x x η ‚ ∂ η / ∂ y = δ y η = η i , j + 1 − η i , j − 1 2 Δ y − Δ y 2 6 δ y y y η 。 ∂ η / ∂ x = δ x η = η i + 1 , j - η i - 1 , j 2 Δ x - Δ x 2 6 δ x x x η ‚ ∂ η / ∂ y = δ y η = η i , j + 1 - η i , j - 1 2 Δ y - Δ y 2 6 δ y y y η 。
数值计算表明这种处理对减小数值格式的频散误差是十分必要的。
《3.1 预报格式》
3.1 预报格式
1) 连续性方程的预报格式
η n + 1 i , j − η n i , j Δ t + δ x [ ( h i , j + η n i , j ) u n i , j ] + δ y [ ( h i , j + η n i , j ) v n i , j ] = 0 ( 1 0 ) η i , j n + 1 - η i , j n Δ t + δ x [ ( h i , j + η i , j n ) u i , j n ] + δ y [ ( h i , j + η i , j n ) v i , j n ] = 0 ( 1 0 )
2) 动量方程的预报格式
u n + 1 i , j − u n i , j Δ t + u n i , j δ x u n i , j + v n i , j δ y u n i , j + g δ x η n i , j = ( 1 + β ) h 2 i , j 3 ( δ x x u n + 1 i , j − δ x x u n i , j Δ t + δ x y v n + 1 i , j − δ x y v n i , j Δ t ) + β g h 2 i , j 3 ( δ x x x η n i , j + δ x y y η n i , j ) + ( 1 + β ) h i , j 2 ⋅ ( 2 δ x h i , j δ x u n + 1 i , j − δ x u n i , j Δ t + δ x x h i , j u n + 1 i , j − u n i , j Δ t + δ x y h i , j v n + 1 i , j − v n i , j Δ t + δ y h i , j δ x v n + 1 i , j − δ x v n i , j Δ t + δ x h i , j δ y v n + 1 i , j − δ y v n i , j Δ t ) + β g h i , j 2 ( 2 δ x h i , j δ x x η n i , j + δ x x h i , j δ x η n i , j + δ x y h i , j δ y η n i , j + δ y h i , j δ x y η n i , j + δ x h i , j δ y y η n i , j ) ( 1 1 ) u i , j n + 1 - u i , j n Δ t + u i , j n δ x u i , j n + v i , j n δ y u i , j n + g δ x η i , j n = ( 1 + β ) h i , j 2 3 ( δ x x u i , j n + 1 - δ x x u i , j n Δ t + δ x y v i , j n + 1 - δ x y v i , j n Δ t ) + β g h i , j 2 3 ( δ x x x η i , j n + δ x y y η i , j n ) + ( 1 + β ) h i , j 2 ⋅ ( 2 δ x h i , j δ x u i , j n + 1 - δ x u i , j n Δ t + δ x x h i , j u i , j n + 1 - u i , j n Δ t + δ x y h i , j v i , j n + 1 - v i , j n Δ t + δ y h i , j δ x v i , j n + 1 - δ x v i , j n Δ t + δ x h i , j δ y v i , j n + 1 - δ y v i , j n Δ t ) + β g h i , j 2 ( 2 δ x h i , j δ x x η i , j n + δ x x h i , j δ x η i , j n + δ x y h i , j δ y η i , j n + δ y h i , j δ x y η i , j n + δ x h i , j δ y y η i , j n ) ( 1 1 )
同理, 可很容易地写出y 方向的动量方程的预报格式。
《3.2 校正格式》
3.2 校正格式
1) 连续性方程的校正格式
η n + 1 i , j − η n i , j Δ t + 1 2 δ x [ ( h i , j + η n + 1 i , j ) u n + 1 i , j + ( h i , j + η n i , j ) u n i , j ] + 1 2 δ y [ ( h i , j + η n + 1 i , j ) v n + 1 i , j + ( h i , j + η n i , j ) v n i , j ] = 0 ( 1 2 ) η i , j n + 1 - η i , j n Δ t + 1 2 δ x [ ( h i , j + η i , j n + 1 ) u i , j n + 1 + ( h i , j + η i , j n ) u i , j n ] + 1 2 δ y [ ( h i , j + η i , j n + 1 ) v i , j n + 1 + ( h i , j + η i , j n ) v i , j n ] = 0 ( 1 2 )
2) 动量方程的校正格式
u n + 1 i , j − u n i , j Δ t + 1 2 [ u n + 1 i , j δ x u n + 1 i , j + u n i , j δ x u n i , j + v n + 1 i , j δ y u n + 1 i , j + v n i , j δ y u n i , j + g δ x ( η n + 1 i , j + η n i , j ) ] = ( 1 + β ) h 2 i , j 3 ( δ x x u n + 1 i , j − δ x x u n i , j Δ t + δ x y v n + 1 i , j − δ x y v n i , j Δ t ) + β g h 2 i , j 6 [ ( δ x x x + δ x y y ) ( η n + 1 i , j + η n i , j ) ] + ( 1 + β ) h i , j 2 ⋅ ( 2 δ x h i , j δ x u n + 1 i , j − δ x u n i , j Δ t + δ x x h i , j u n + 1 i , j − u n i , j Δ t + δ x y h i , j v n + 1 i , j − v n i , j Δ t + δ y h i , j δ x v n + 1 i , j − δ x v n i , j Δ t + δ x h i , j δ y v n + 1 i , j − δ y v n i , j Δ t ) + β g h i , j 4 ⋅ [ 2 δ x h i , j δ x x ( η n + 1 i , j + η n i , j ) + δ x x h i , j δ x ( η n + 1 i , j + η n i , j ) + δ x y h i , j δ y ( η n + 1 i , j + η n i , j ) + δ y h i , j δ x y ( η n + 1 i , j + η n i , j ) + δ x h i , j δ y y ( η n + 1 i , j + η n i , j ) ] ( 1 3 ) u i , j n + 1 - u i , j n Δ t + 1 2 [ u i , j n + 1 δ x u i , j n + 1 + u i , j n δ x u i , j n + v i , j n + 1 δ y u i , j n + 1 + v i , j n δ y u i , j n + g δ x ( η i , j n + 1 + η i , j n ) ] = ( 1 + β ) h i , j 2 3 ( δ x x u i , j n + 1 - δ x x u i , j n Δ t + δ x y v i , j n + 1 - δ x y v i , j n Δ t ) + β g h i , j 2 6 [ ( δ x x x + δ x y y ) ( η i , j n + 1 + η i , j n ) ] + ( 1 + β ) h i , j 2 ⋅ ( 2 δ x h i , j δ x u i , j n + 1 - δ x u i , j n Δ t + δ x x h i , j u i , j n + 1 - u i , j n Δ t + δ x y h i , j v i , j n + 1 - v i , j n Δ t + δ y h i , j δ x v i , j n + 1 - δ x v i , j n Δ t + δ x h i , j δ y v i , j n + 1 - δ y v i , j n Δ t ) + β g h i , j 4 ⋅ [ 2 δ x h i , j δ x x ( η i , j n + 1 + η i , j n ) + δ x x h i , j δ x ( η i , j n + 1 + η i , j n ) + δ x y h i , j δ y ( η i , j n + 1 + η i , j n ) + δ y h i , j δ x y ( η i , j n + 1 + η i , j n ) + δ x h i , j δ y y ( η i , j n + 1 + η i , j n ) ] ( 1 3 )
同理, 可很容易地写出y 方向的动量方程的校正格式。需要说明的是:在预报阶段, η 的求解是显式进行的, 而u 和v 的离散方程可化为三对角矩阵形式采用追赶法求解;在校正阶段, η , u 和v 的求解需要进行联合迭代, 直到相邻两次迭代的结果差值满足某一精度要求为止, 采用的收敛准则为
max ⎡ ⎣ ∑ i , j ∣ ∣ η n + 1 i , j ∣ ∣ k + 1 − η n + 1 i , j ∣ ∣ k ∣ ∣ ∑ i , j η n + 1 i , j ∣ ∣ k + 1 ‚ ∑ i , j ∣ ∣ u n + 1 i , j ∣ ∣ k + 1 − u n + 1 i , j ∣ ∣ k ∣ ∣ ∑ i , j u n + 1 i , j ∣ ∣ k + 1 ‚ ∑ i , j ∣ ∣ v n + 1 i , j ∣ ∣ k + 1 − v n + 1 i , j ∣ ∣ k ∣ ∣ ∑ i , j v n + 1 i , j ∣ ∣ k + 1 ⎤ ⎦ < ε ( 1 4 ) max [ ∑ i , j | η i , j n + 1 | k + 1 - η i , j n + 1 | k | ∑ i , j η i , j n + 1 | k + 1 ‚ ∑ i , j | u i , j n + 1 | k + 1 - u i , j n + 1 | k | ∑ i , j u i , j n + 1 | k + 1 ‚ ∑ i , j | v i , j n + 1 | k + 1 - v i , j n + 1 | k | ∑ i , j v i , j n + 1 | k + 1 ] < ε ( 1 4 )
式中∑ i , j F ∑ i , j F F 进行求和;k 表示第k 次迭代;ε 为一小的正数, 表示误差容许限。
《4 初始条件和边界条件》
4 初始条件和边界条件
《4.1 初始条件》
4.1 初始条件
u ( x , y , 0 ) = v ( x , y , 0 ) = η ( x , y , 0 ) ( 1 5 ) u ( x , y , 0 ) = v ( x , y , 0 ) = η ( x , y , 0 ) ( 1 5 )
《4.2 边界条件》
4.2 边界条件
根据具体问题的不同, 可分为各种不同的边界条件。对算例, 边界条件有以下几种形式。
1) 上游造波边界条件
根据入射波的类型不同, 可以给定不同的造波边界条件。对线性入射波, 造波边界条件给定为:
η ( 0 , y , t ) = A 0 cos ( k c o s ? α 0 x + k s i n ? α 0 y − ω t + 0 . 5 π ) ( 1 6 ) u ( 0 , y , t ) = ω η ( 0 , y , t ) cos ? α 0 k [ η ( 0 , y , t ) + h ( 0 , y ) ] ( 1 7 ) v ( 0 , y , t ) = ω η ( 0 , y , t ) sin ? α 0 k [ η ( 0 , y , t ) + h ( 0 , y ) ] ( 1 8 ) η ( 0 , y , t ) = A 0 cos ( k c o s ? α 0 x + k s i n ? α 0 y - ω t + 0 . 5 π ) ( 1 6 ) u ( 0 , y , t ) = ω η ( 0 , y , t ) cos ? α 0 k [ η ( 0 , y , t ) + h ( 0 , y ) ] ( 1 7 ) v ( 0 , y , t ) = ω η ( 0 , y , t ) sin ? α 0 k [ η ( 0 , y , t ) + h ( 0 , y ) ] ( 1 8 )
式中A 0 为入射波振幅;α 0 为初始入射波向, 算例中α 0 =0;ω 为波浪圆频率;h 为水深;k 为波数, 由该Boussinesq方程的线性色散关系
[19 ]
ω 2 = g k k h ( 1 + β k 2 h 2 / 3 ) 1 + ( 1 + β ) k 2 h 2 / 3 ( 1 9 ) ω 2 = g k k h ( 1 + β k 2 h 2 / 3 ) 1 + ( 1 + β ) k 2 h 2 / 3 ( 1 9 )
2) 下游开边界条件
下游开边界条件采用辐射边界条件, 其具体形式如下:
∂ η ∂ x + cos ? θ c ∂ η ∂ t = 0 ‚ ∂ η ∂ y + sin ? θ c ∂ η ∂ t = 0 ( 2 0 ) ∂ u ∂ x + cos ? θ c ∂ u ∂ t = 0 ‚ ∂ v ∂ y + sin ? θ c ∂ v ∂ t = 0 ( 2 1 ) ∂ η ∂ x + cos ? θ c ∂ η ∂ t = 0 ‚ ∂ η ∂ y + sin ? θ c ∂ η ∂ t = 0 ( 2 0 ) ∂ u ∂ x + cos ? θ c ∂ u ∂ t = 0 ‚ ∂ v ∂ y + sin ? θ c ∂ v ∂ t = 0 ( 2 1 )
3) 侧向开边界条件
侧向 (y 方向) 开边界条件严格说来也应采用相应的辐射边界条件, 但为了简化处理, 这里采用了下面的近似边界条件:
∂ η ∂ y = 0 ‚ ∂ u ∂ y = 0 ‚ ∂ v ∂ y = 0 ( 2 4 ) ∂ η ∂ y = 0 ‚ ∂ u ∂ y = 0 ‚ ∂ v ∂ y = 0 ( 2 4 )
《5 算例及结果分析》
5 算例及结果分析
Beji和Battjes
[20 ] 及Luth等分别对淹没浅堤地形上的波浪传播进行了物理模型试验, 其试验结果被一些学者用来验证他们的数值模型
[11 ,18 ,21 ] 。笔者采用的资料取自文献[18 ] , 试验地形、水深及测点布置如图1所示, 入射波周期为T =2.02 s, 振幅为A 0 =0.01 m。数值计算时, y 方向的计算长度取1 m, 采用的空间步长为Δx =Δy =0.05 m, 时间步长为Δt =0.02 s, 共计算1700个时间步, 每步迭代的误差容许限为ε =10-5 。数值计算采用这样小的空间步长和时间步长是为了能够比较精细地分辨出由于非线性作用产生的高频谐波。
《图1》
图1试验地形示意
Fig.1 Schematic of experimental topography
图2给出了在10个断面处的波面随时间变化的计算和实测值的比较。从总体上看, 该方程能较好地模拟波浪在潜堤上传播时波面的变形过程。可以看出, 在x =2.0, 4.0, 10.5, 12.5, 13.5, 14.5 m处的计算结果和实测值均相当吻合, 在x =15.7, 17.3, 19.0 m处的计算结果和实测值比较吻合, 但存在一定差异。而x =5.7处的数值结果和实测值存在较大的相位差, 文献[18 ] 采用四阶完全非线性Boussinesq方程得到的结果与实测值之间也存在同样相位差, 该文还说明了产生这种相差的主要原因是试验系统误差所致, 即测点位置在试验中实际上不是在5.7 m处, 而是在其他某个位置。在x =15.7, 17.3, 19.0 m处的计算结果误差可能来自2个主要方面:a. 算例所采用的Boussinesq方程为弱非线性方程, 其本身可能不能足够精细地反映由于非线性作用而产生的高频谐波, 从而产生一定的误差;b. 数值模拟采用的数值格式的误差。文献[18 ] 还采用二阶完全非线性Boussinesq方程对该地形上的波浪变形进行了数值模拟, 得到的结果与实测值之间也存在一定差异。文献[12 ] 应用Nwogu
[10 ] 提出的Boussinesq方程并采用一种时间和空间均为四阶精度的数值格式对这种地形上的波浪变形进行了数值模拟, 其计算结果和实测值也存在一定差异, 主要表现在数值模拟的波高峰值偏小。文献[12 ] 在对计算结果和实测值比较后, 认为Nwogu提出的Boussinesq方程可用于中等水深地形上的波浪变形预报, 其结果在一定程度上是可信的, 但非线性项的作用过大, 需要对其进行修正;然而文献[15 ] 通过对Nwogu
[10 ] 提出的Boussinesq方程进行抛物化处理, 得到了一个含自由参数的抛物型方程, 并采用该抛物型方程对该地形上的波浪问题也进行了数值模拟, 该方程却得到了比文献[12 ] 更接近实测值的计算结果。
《图2》
图2 淹没浅堤上波浪变形的数值解与实测值的比较
Fig.2 Numerical solution and experimental data of wave transformation over a submerged bar
《6 结论》
6 结论
对Beji和Nadaoka提出的Boussinesq方程采用预报校正格式进行了数值离散, 并对淹没浅堤上的波浪变形进行了数值模拟, 从数值模拟结果和实测值的比较来看, Beji和Nadaoka提出的Boussinesq方程能给出比较满意的结果, 可以用于实际复杂地形上的波浪问题计算。另一方面, 数值格式精度还不够高, 带来一定的数值误差, 该Boussinesq方程为弱非线性弱频散模型, 这种改进的Boussinesq方程及其求解方法需做进一步的完善。
京公网安备 11010502051620号
