【浅水模型MATLAB】尝试复刻SCI论文中的溃坝流算例
【浅水模型MATLAB】尝试复刻SCI论文中的溃坝流算例
- 前言
- 问题描述
- 控制方程及数值方法
- 浅水方程及其数值计算方法
- 边界条件的实现
- 代码框架与关键代码
- 模拟结果
更新于2024年9月17日
前言
这篇博客算是学习浅水方程，并利用MATLAB复刻Liang (2004)1中溃坝流算例的一个记录。
二维溃坝流（Dam Break）问题是浅水模型经典的一个测试算例，它测试了模型对急变流的模拟效果、以及对干-湿边界处理方法的有效性。相比于之前的模拟算例，本算例中需要重点处理的问题是：
- 模型的内边界的处理；
- 干-湿边界的处理。
本博客将着重解决第一个问题，而先不考虑第二个问题，即设置的模拟算例不含干-湿边界的处理。此外，本算例中涉及的控制方程与数值方法已经在《【浅水模型MATLAB】尝试完成一个数值模拟竞赛题》中介绍；不清楚的朋友可参考该博客内容。
如果你习惯用别的代码，也想做类似的建模尝试，十分欢迎一起交流！（最近有点想转python代码了，希望感兴趣的同志来一起交流）
如果各位朋友发现了文章或代码中的错误，亦或是改进之处，请不吝赐教，欢迎大家留言，一起改进模型！本博客文章将持续更新，上面也会标注提出改进建议的同志们。（不过，本人最近在忙活毕业论文，可能更新不及时）
同时，想要完整代码的朋友请联系我，我可无偿提供脚本文件。
希望同大家一起进步！
问题描述
本算例的计算区域为一个200m×200m的矩形平底水槽，水槽的四周都是垂直的固体壁面。如下图所示，计算域被分成x<100m和x>100m的左右两个部分；左侧初始水深为10m，右侧初始水深为5m。左右两个区域被两道平行于y方向的壁面阻隔，仅在95m<y<170m的区域联通。在模拟开始时，左侧水体会突然通过x=100m，95m<y<170m的区域向着右侧下泄，形成溃坝流。此外，模型中的所有壁面都是光滑的。
控制方程及数值方法
浅水方程及其数值计算方法
二维浅水方程的形式及其具体求解内容详见Liang的论文2和博客《【浅水模型MATLAB】尝试完成一个数值模拟竞赛题》。模型采用Godunov型有限体积法，通过一系列的处理，方程也保证了静水状态时压力与底坡源项的平衡。
此外，模型中的水力变量都定义在网格的中心位置。网格边界处的通量采用HLL求解器获得。
边界条件的实现
计算域的外边界均为无通量的free-slip闭合边界，边界处的法向速度和通量均被定义为0。在求解过程中，可将边界处的水力参数设置为其临近网格相同的物理量的值。
对于模型在x=100m处的内边界，模型需要定义其对应边界的通量为零。具体处理方式如下图所示。在对内边界左侧的相邻网格进行线性重构及通量计算时，需要通过一个辅助计算的虚网格，该虚网格有着和左侧相邻网格i相同的水力变量值。同理，在对内边界右侧的相邻网格进行线性重构及通量计算时，也需要通过一个辅助计算的虚网格，该虚网格有着和右侧 相邻网格i相同的水力变量值。由于本模型采用了Minmod的限制器，所以此种处理会使得内边界对应的左侧变量UL =Ui，而使右侧变量UR =Ui+1。
代码框架与关键代码
我的模型代码主要分为参数设置、网格构建、初始化、主循环和其余函数等五个部分。
- 设置物理参数、网格参数、时间参数等。代码如下所示：
grav = 9.81; % Gravitational acceleration
rho = 1000; % Density
CFL = 0.5; % Courant Number
Lx = 200; % Length of the domain
Ly = 200; % Width of the domain
zb0 = 0.0; % Bottom elevation
n = 0.00; % Manning coefficient
h_dry = 0.02; % wet-dry threshold value
dx = 1; % Grid spacing
dy = 1; % Grid spacing
dt = 0.05; % Time spacing at the first step
dtmax = 0.1; % allowed max time step (s)
tend = 10.0; % End of the simulation time
plot_int = 0.5; % Time interval to next plot
我设置网格为边长1m的正方形，底高程为zb0=0.0。最大允许的Courant数设置为0.5，初始时间步为0.05s，之后的每一个时间步通过CFL条件计算得到。
- 网格构建：网格有两个要素需要定义，一是网格的四个节点（Xp和Yp），二是网格的中心点（Xc和Yc）；网格中心点也即水力物理量定义的位置。代码略。
- 初始化：设置底高程zb=0，计算zb的梯度zbx和zby；设置左右区域的初始水位，之后再设置流速u、v为零。
- 主循环：（1）计算网格边界处的水位、水深、流速值；（2）设置内边界条件；（3）计算通量项FL、FR、GL和GR；（4）利用HLL求解通量F和G；（5）计算源项S；（6）计算新一个时间步的eta、h、u和v。除了上述步骤（2）和（3）其余计算过程与博客《【浅水模型MATLAB】尝试完成一个数值模拟竞赛题》中代码基本一致；涉及的关键代码如下：
while(t<tend)
% estimate the dt
dtx = dx./(abs(u)+sqrt(grav*h) + 1E-8);
dty = dy./(abs(v)+sqrt(grav*h) + 1E-8);
dt1 = min(min(dtx,[],"all"), min(dty,[],"all"));
dt = min(dtmax, CFL*dt1);
clear dt1 dtx dty
etan = eta; hn = h;
un = u; vn = v;
% 2rd-order Runge-Kutta Method
for k = 1:2
% 1. reconstruct the flow data
% 1.1 x-direction reconstruction (Natural closed boundary)
% 求解网格边界处的水位exL和exR，流速uxL、uxR、vxL和vxR；
% ...
% 1.2 y-direction reconstruction (Natural closed boundary)
% 求解网格边界处的水位eyL和eyR，流速uyL、uyR、vyL和vyR；
% ...
% 2. inner boundary conditions
ff = find((yc<=95) + (yc>=170));
% 2.1 left cells
de = minmod((eta(:,Nx/2)-eta(:,Nx/2))/dx, ...
(eta(:,Nx/2)-eta(:,Nx/2-1))/dx);
du = minmod((u(:,Nx/2)-u(:,Nx/2))/dx, ...
(u(:,Nx/2)-u(:,Nx/2-1))/dx);
dv = minmod((v(:,Nx/2)-v(:,Nx/2))/dx, ...
(v(:,Nx/2)-v(:,Nx/2-1))/dx);
exR(ff,Nx/2) = eta(ff,Nx/2) - 0.5*dx*de(ff);
exL(ff,Nx/2+1) = eta(ff,Nx/2) + 0.5*dx*de(ff);
clear de du dv
% 2.2 right cells
de = minmod((eta(:,Nx/2+2)-eta(:,Nx/2+1))/dx, ...
(eta(:,Nx/2+1)-eta(:,Nx/2+1))/dx);
du = minmod((u(:,Nx/2+2)-u(:,Nx/2+1))/dx, ...
(u(:,Nx/2+1)-u(:,Nx/2+1))/dx);
dv = minmod((v(:,Nx/2+2)-v(:,Nx/2+1))/dx, ...
(v(:,Nx/2+1)-v(:,Nx/2+1))/dx);
exR(ff,Nx/2+1) = eta(ff,Nx/2+1) - 0.5*dx*de(ff);
exL(ff,Nx/2+2) = eta(ff,Nx/2+1) + 0.5*dx*de(ff);
clear ff de du dv
% 3. flux terms (F and G)
F1L = hxL.*uxL;
F2L = hxL.*uxL.*uxL + 0.5*grav*(exL.*exL - ...
exL.*(zbp(1:end-1,:)+zbp(2:end,:)));
F3L = hxL.*uxL.*vxL;
F1R = hxR.*uxR;
F2R = hxR.*uxR.*uxR + 0.5*grav*(exR.*exR - ...
exR.*(zbp(1:end-1,:)+zbp(2:end,:)));
F3R = hxR.*uxR.*vxR;
G1L = hyL.*vyL;
G2L = hyL.*uyL.*vyL;
G3L = hyL.*vyL.*vyL + 0.5*grav*(eyL.*eyL - ...
eyL.*(zbp(:,1:end-1)+zbp(:,2:end)));
G1R = hyR.*vyR;
G2R = hyR.*uyR.*vyR;
G3R = hyR.*vyR.*vyR + 0.5*grav*(eyR.*eyR - ...
eyR.*(zbp(:,1:end-1)+zbp(:,2:end)));
% 4. calculate the flux by HLL
[sxL sxR] = WaveSpeed(hxL, hxR, uxL, uxR);
F1 = HLLSolver(F1L, F1R, sxL,sxR, exL,exR);
F2 = HLLSolver(F2L, F2R, sxL,sxR, hxL.*uxL,hxR.*uxR);
F3 = HLLSolver(F3L, F3R, sxL,sxR, hxL.*vxL,hxR.*vxR);
[syL syR] = WaveSpeed(hyL, hyR, vyL, vyR);
G1 = HLLSolver(G1L, G1R, syL,syR, eyL,eyR);
G2 = HLLSolver(G2L, G2R, syL,syR, hyL.*uyL,hyR.*uyR);
G3 = HLLSolver(G3L, G3R, syL,syR, hyL.*vyL,hyR.*vyR);
clear sxL sxR syL syR F1L F1R F2L F2R F3L F3R G1L G1R G2L G2R G3L G3R
% 4.1. west boundary
% ...
% 4.2. east boundary
% ...
% 4.3. south boundary
% ...
% 4.4. north boundary
% ...
% 4.5. inner boundary
ff = find((yc<=95) + (yc>=170));
F1(ff,Nx/2+1) = 0; F3(ff,Nx/2+1) = 0;
F2_L(ff,1) = 0.5*grav*(exL(ff,Nx/2+1).^2 - ...
exL(ff,Nx/2+1).*(zbp(ff+1,Nx/2+1)+zbp(ff,Nx/2+1)));
F2_R(ff,1) = 0.5*grav*(exR(ff,Nx/2+1).^2 - ...
exR(ff,Nx/2+1).*(zbp(ff+1,Nx/2+1)+zbp(ff,Nx/2+1)));
clear ff exL exR eyL eyR hxL hxR hyL hyR uxL uxR uyL uyR vxL vxR vyL vyR
% 5. source terms
% 计算S的三个分量S1、S2和S3
% ...
% 6. time stepping
% 6.1 solve eta
% ...
% 6.2 solve h*u
% ...
% for inner boundary
ff = find((yc<=95) + (yc>=170));
hu(ff,Nx/2) = h(ff,Nx/2).*u(ff,Nx/2) ...
- dt/dx*(F2_L(ff) - F2(ff,Nx/2)) ...
- dt/dy*(G2(ff+1,Nx/2)-G2(ff,Nx/2)) + dt*S2(ff,Nx/2);
hu(ff,Nx/2+1) = h(ff,Nx/2+1).*u(ff,Nx/2+1) ...
- dt/dx*(F2(ff,Nx/2+2) - F2_R(ff)) ...
- dt/dy*(G2(ff+1,Nx/2+1)-G2(ff,Nx/2+1)) + dt*S2(ff,Nx/2+1);
clear ff F2_L F2_R
% 6.3 solve h*v
% ...
% ...
end
% 计算得到本时间步的h、u和v
% 7. plot
% ...
end
end
- 其余子函数：包括minmod限制器、HLL求解器等。代码略。
模拟结果
1.水位结果
2.流速结果（颜色表示合速度的大小，箭头表示速度方向）
3. 三维水面图
Liang, Q., Borthwick, A.G.L. and Stelling, G. (2004), Simulation of dam- and dyke-break hydrodynamics on dynamically adaptive quadtree grids. Int. J. Numer. Meth. Fluids, 46: 127-162. ↩︎
Liang Q , Marche F .Numerical resolution of well-balanced shallow water equations with complex source terms[J].Advances in Water Resources, 2009, 32(6):873-884. ↩︎