NHWAVE孤立波爬坡测试算例（Solitary Wave on a Simple Beach）
本blog介绍了NHWAVE代码包中的几个基准测试算例。这几个算例都基于同一个孤立波传播、爬坡的基本模型，且包含了模型结果与解析解、实验结果的比较。
此外，本人在NHWAVE原版源代码基础上，更新了垂向网格分层百分比可指定、被动示踪剂输移、温度场输入输出等模块，并修正了植被水流模型中的一些BUG。有兴趣的朋友请上我的GitHub下载代码。
测试算例简介
该算例模拟了一维孤立波的传播和在斜面上爬升，并将它们与该问题的解析解和实验室结果相比较，以验证模型的有效性。关于该实验的解析解与实验室结果详见 Synolakis (1987)1。
模拟所需参数文件和初始水深文件在代码包 /examples/landslide_benchmarks/inundation_benchmarks 目录下。其中，与解析解相关的模拟文件在 Solitary_wave_on_a_simple_beach_analytical 文件夹中，与实验室实验相关的模拟文件在 Solitary_wave_on_a_simple_beach_laboratory 文件夹中。
实验设置
模型区域为一个带斜坡的数值水槽，实验设置及几何参数如下图所示。实验中所研究的孤立波为上凸波，高度为H，初始水深为d；在初始时刻，其中心位置距离坡脚L。参考系坐标原点位于离坡脚水平距离x0的位置，角度β = arccot(19.85)。注意，这个x的方向是朝左的。
与解析解对比的模拟（Analytical Solution）
在此数值实验中，这些几何参数的关系如下：
L
=
d
⋅
a
r
c
c
o
s
h
(
20
)
/
γ
η
(
x
,
0
)
=
H
s
e
c
h
2
(
γ
(
x
−
x
1
)
/
d
)
x
1
=
x
0
+
L
γ
=
3
H
/
4
d
L=d \cdot arccosh( \sqrt{20})/γ \\ η(x, 0) = H sech^2 (γ(x − x_1)/d) \\ x_1 = x_0+L \\ γ = \sqrt{3H/4d}
L=d⋅arccosh(20)/γη(x,0)=Hsech2(γ(x−x1)/d)x1=x0+Lγ=3H/4d
本算例的重点是对入射的非破碎孤立波进行模拟，并满足条件 H/d = 0.019。此外，波的初始速度满足：
u
(
x
,
0
)
=
−
g
/
d
η
(
x
,
0
)
.
u(x, 0) = −\sqrt{g/d} \space η(x, 0).
u(x,0)=−g/d η(x,0).
本内容可被分成以下两部分：
- 模拟孤立波的最大爬升高度，并与解析解比较（runup_law）
- 模拟水位的时空变化，并与解析解比较（time_solution）
最大爬升高度的模拟（max runup）
设置初始水深d = 0.5m，波高H/d=0.03；水槽总长为60m，其中平床水槽长 L = 40m。对应模型的参数设置详见 /runup_law/input-1_d0.5_Hd0.03_beta10 文件夹中的 input.txt 和 depth.txt。
在模型中，水平向网格尺寸为 Δx = 0.05 m，垂向上分为4个σ层。计算域左侧为孤立波的输入边界，底面为自由滑移边界。此外，模型中还采用了RNG k-ε紊流模型。
模拟结束后，可利用 /runup_law 文件夹中的两个MATLAB脚本分析模拟结果。
首先是波浪最大爬坡高度的模拟结果与解析解的对比。打开并运行RUNUPLINE.m，得到如下图所示的对比结果图。
可见，模拟结果的精度很高。此外，output_hmax.m脚本还计算了模拟结果与解析解之间的对比误差；先把结果文件中的hmax_00400复制到与output_hmax.m同目录下，后运行output_hmax.m，得到相对偏差5.5023%。
水位的时空变化模拟（Water profiles and time histories of surface elevation）
本例模拟了一个三维水槽，其总长为75.0m，平床水槽长 L = 40.25m，宽为0.5m。设置水平向网格尺寸Δx = Δy = 0.05 m。水槽的左侧为孤立波入射边界，其余侧面和底面均为自由滑移边界。此外，模型中还采用了RNG k-ε紊流模型。对于左侧的入射波，设置初始水深d = 1.0 m，波高H/d = 0.019。对应模型的参数设置详见 /time_solution/input 文件夹中的 input.txt 和 depth.txt。
模拟结束后，可利用 /time_solution 文件夹中的MATLAB脚本分析模拟结果。
首先是不同时刻水面线的比较，如下图所示，实线为模拟结果（通过运行ETA_Liner.m提取结果数据），红圈为解析解（数据来自/time_solution 文件夹中的analytical_SW.mat）。
另一个脚本Time_Liner.m对比了 x/d = 0.25（红）和9.95（黑）处水位随时间的变化。其中实线是模拟结果，数据来自结果文件probe_0001和probe_0002；圆圈为实验结果，数据来自/time_solution 文件夹中的analytical_RU.mat。
与实验室结果对比的模拟
本案例对应的物理实验来自加州理工大学波浪水槽。该水槽长31.73m、深60.96cm、宽39.97cm。水箱底部由涂漆的不锈钢板组成。在水槽的一端有一个斜坡，坡道的坡度为1:19.85（即cot β =19.85). 。实验装置的立面示意图与上图相同；其中的几何参数为：
L
=
d
⋅
a
r
c
c
o
s
h
(
20
)
/
γ
γ
=
3
H
/
4
d
c
o
s
h
−
1
a
=
l
o
g
(
a
+
a
2
−
1
)
η
(
x
,
0
)
=
H
s
e
c
h
2
(
γ
(
x
−
x
1
)
/
d
)
x
1
=
x
0
+
L
L=d \cdot arccosh( \sqrt{20})/γ \\ γ = \sqrt{3H/4d}\\ cosh^{-1} a=log(a+\sqrt{a^2 - 1})\\ η(x, 0) = H sech^2 (γ(x − x_1)/d) \\ x_1 = x_0+L
L=d⋅arccosh(20)/γγ=3H/4dcosh−1a=log(a+a2−1)η(x,0)=Hsech2(γ(x−x1)/d)x1=x0+L
初始条件为：
η
(
x
,
0
)
=
H
/
c
d
o
t
s
e
c
h
2
(
γ
(
x
−
x
1
)
/
d
)
u
(
x
,
0
)
=
−
η
(
x
,
0
)
⋅
g
/
d
η(x, 0) = H /cdot sech^2 (γ(x − x_1)/d) \\ u (x,0) = -\eta(x,0) \cdot \sqrt{g/d}
η(x,0)=H/cdotsech2(γ(x−x1)/d)u(x,0)=−η(x,0)⋅g/d
在模型中，我们依旧建立一个长为75.0m的右侧带斜坡的立面二维数值水槽。其中，平床段长 L = 40.25 m。
实验中设置平床段的初始水深为1.0 m，H/d分别为0.0185和0.30。对于边界条件，左侧为孤立波输入边界，底面粗糙高度为0.0001m。此外，模型采用了RNG k-ε紊流模型。
对应模型的参数设置详见 /examples/Solitary_wave_on_a_simple_beach_laboratory/H_d_0.0185/input 文件夹和 /examples/Solitary_wave_on_a_simple_beach_laboratory/H_d_0.030/input 文件夹中的 input.txt 和 depth.txt。模拟结束后，可利用文件夹中的MATLAB脚本分析模拟结果。
首先是对于H/d = 0.0185的算例，/H_d_0.0185/ETA_Liner_BP1_A.m脚本能够提取几个时间段的水面线数据；模拟结果与实验结果（实验结果数据来自SW_01_00185.mat中）的比较如下图所示。图中实线表示模拟结果，小圆圈表示实验结果。
同理，对于 H/d = 0.30 算例，可以得到如下图的结果比较。其中的模型结果采用 /H_d_0.03/ETA_Liner_BP1_B.m 脚本提取，实验数据来自 SW_01_03.mat 文件。图中实线表示模拟结果，小圆圈表示实验结果。
Synolakis C E . The runup of solitary waves[J]. J Fluid Mechanics, 1987, 185:523-545. ↩︎
1590
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
