复杂形状海滩上的海啸波爬高模拟
本案例来自微信公众号“陆面体科技”
请点击关注:
原作者:Stéphane Popinet
版本:110107
运行时间:50min
案例背景
本算例是海啸模型的经典验证算例,该算例是在“第三届长波波浪爬高模型国际研讨会”用于对比验证的标准模型。本文中使用的海啸模型是基于造波水池数据的提出,用于研究1993年日本奥尻岛海啸期间在藻内地区附近观察到的极端海啸波爬高现象。
海啸波爬高模拟是求解非线性浅水方程组的一个具体应用,该方程组由Navier-Stokes方程组在浅水条件下推导而出,描述具有自由表面的浅水体中水流运动规律的偏微分方程组,1871年由法国科学家Adhémar Jean Claude Barré de Saint-Venant提出,故又名圣维南方程组,该方程在气动、水力以及其他工程领域中有着广泛的应用,可以很好地模拟出浅水波的各种特性。
我们的中国船舶科学研究中心(China Ship Scientific Research Center,又称中国船舶重工集团公司第七〇二研究所),基于长期的研究积累,在国际上首先创立并发展完善了对数值水池的系统论证,较全面地阐述了数值水池的定义,并已开展了相关技术的研发与应用。
既然说到了浅水波,就避不开谈论船舶工程这个领域。在船舶工程中,有一个很重要的研究成果—数值水池,即基于数值仿真技术的虚拟造波水池。该成果是在模型驱动、知识驱动的船舶与海洋结构物设计优化流程中,基于CFD应用技术及专家知识,依托高性能计算平台,为船舶与海洋工程水动力学性能设计、预报、评估和优化提供高效能、高精度、高置信度的沉浸式、情景化虚拟试验服务的应用技术系统。
言归正传,我们来看看法国的开源软件Gerris创始人-- Stéphane Popinet是如何模拟复杂形状海滩的海啸问题的。
图1: 藻内地区海啸水面(浅蓝色)和岸底(白色)动画
首先来看图1的动画:模拟了海滩的几何形状以及海啸随时间演变的过程,水深数据和海滩几何形状与实验所用的造波水池一致。实验过程中在水面的开放边界处施加波动(图1右侧画面外)。在动画中可以清楚地看到图中央部谷地在开始时的干涸状态和之后海啸波爬高至极值的状态,以及边界处的波浪反射和由于水下涡旋而引起的水面浅凹。
图2:该过程的另一个视图以及用于计算流场的自适应网格。
图2 自由表面高度(左列)和相应的自适应网格的演变(右列)。图中给出了地形轮廓等高线,左列图中白色区域是没有水的,有水区域根据自由面高度(相对于未受干扰的水面)着色。用于对比验证的实验数据包括了在多个位置不同时间的自由面高度(见图3,t = 18 s)。
图3,4和5给出了各个时刻实验和数值计算结果的比较。
图3:在探针5的位置测量和计算的自由面高度对比
图4:在探针7的位置测量和计算的自由面高度对比
图5:在探针9的位置测量和计算的自由面高度对比
最后我们来看看实验俯瞰视频和相应的模拟分析视频对比:
图6: 实验视频(左边)和模拟结果(右边)之间的对比。
案例解析
看完结果,我们来看看Stéphane的代码是否和他的结果一样优雅漂亮美丽大方~
首先,Gerris是在Linux系统下运行的开源软件。区别于我们熟知的OpenFoam和SU2,Gerris的案例结果后处理难以直接在传统的后处理软件中解决。所以Gerris的爸爸Stéphane就针对Gerris的结果处理专门写了另外一个开源软件– GFSview. 简单地来说,则是用GFSview可以打开Gerris的结果文件。GFSview是一款完整的后处理软件,用来观察各种参数(速度、温度等)的计算结果,也可以调整视角方向、画面大小的颜色范围及风格,还可以批量处理所需要的计算结果等。
接下来,我们来看看文件的调动和重要参数设置:
仿真的配置文件:
# below this depth the flow is considered
"dry"
Define DRY 1e-4
# use the GfsRiver Saint-Venant solver
# shift the origin of the reference box to (0,0)
2 1 GfsRiver GfsBox GfsGEdge
{
x = 0.5
y = 0.5
}
{
# the domain is 3.402 m X 6.804 m
# units for time are seconds
PhysicalParams
{
L = 3.402
g = 9.81
}
Refine 6
# maintain the Zb1 variable using the GTS surface
VariableFunction Zb1 bathy.gts
# the initial water level is at z = 0, so the depth P is...
Init {}
{
P = MAX(0., -Zb1)
}
# use a Sweby limiter rather than the default minmod which is too
# dissipative
AdvectionParams
{
gradient = gfs_center_sweby_gradient
}
# adapt down to 9 levels based on the slope of the
(wet)
# free-surface and with a tolerance of 1 mm
AdaptGradient
{
istep = 1
}
{
cmax = 1e-3
cfactor = 2
maxlevel = 9
minlevel = 6
}
(
P < DRY ? 0. : P + Zb
)
# at each timestep
Init
{
istep = 1
}
{
# Add a "shelf" to simulate the wall on the right-hand-side boundary
Zb =
(
x > 5.448 ? 0.13535 : Zb1
)
# read in the experimental timeseries in variable 'input'
input = input.cgd
# implicit quadratic bottom friction with coefficient 1e-3
U = (P > DRY ? U/(1. + dt*1e-3*Velocity/P) : 0.)
V = (P > DRY ? V/(1. + dt*1e-3*Velocity/P) : 0.)
P = (P > DRY ? P : 0.)
}
Time
{
end = 22.5
}
OutputTime
{
istep = 10
}
stderr
OutputTiming
{
start = end
}
stderr
OutputSimulation
{
step = 1
}
stdout
OutputSimulation
{
step = 1
}
sim-%g.gfs
EventScript
{
step = 1
}
{
gzip -f
sim-*.gfs
}
# output data at probe locations
OutputLocation
{
istep = 1
}
input 1e-3 1.7 0
OutputLocation
{
istep = 1
}
p5 4.521 1.196 0
OutputLocation
{
istep = 1
}
p7 4.521 1.696 0
OutputLocation
{
istep = 1
}
p9 4.521 2.196 0
# kinetic energy
OutputScalarSum
{
istep = 1
}
ke
{
v = (P > 0. ? Velocity2*P : 0.)
}
# free-surface elevation
OutputScalarStats
{
istep = 1
}
p
{
v = (Zb > 0. ? P : P + Zb)
}
OutputScalarNorm
{
istep = 1
}
u
{
v = Velocity
}
OutputTime
{
istep = 1
}
balance
OutputBalance
{
istep = 1
}
balance
# generate movies
GModule gfsview
OutputView
{
start = 9 step = 0.0416
}
{
ppm2mpeg -s
640x480 > monai.mpg
}
{
width = 1280
height = 960
}
3D.gfv
OutputView
{
start = 14.63
end = 19.5
step = 0.033333333
}
{
ppm2mpeg -s
400x600 > overhead.mpg
}
{
width = 800
height = 1200
}
overhead.gfv
}
{
dry = DRY
}
GfsBox
{
# use 'subcritical' boundary condition to impose the experimental
# water level on the left boundary
left = Boundary
{
BcSubcritical U
(input - Zb)
}
top = Boundary
bottom = Boundary
}
GfsBox
{
top = Boundary
bottom = Boundary
}
1 2 right
参考文献
S. Popinet. Quadtree-adaptive tsunami modelling. Ocean Dynamics, 61(9):1261–1285, 2011
在Gfsview中根据自己的需要, 是预先调整好的后处理输出信息,视角等,将该设置以gfv文件后缀储存。这样就可以通过Gfsview批量处理Gerris仿真生成的数据了。
感兴趣的朋友可以在后台跟陆陆索取更多的资料~