【浅水模型MATLAB】尝试复刻SCI论文中的溃坝流算例

前言
问题描述
控制方程及数值方法
- 浅水方程及其数值计算方法
- 边界条件的实现
代码框架与关键代码
模拟结果

更新于2024年9月17日

前言

这篇博客算是学习浅水方程，并利用MATLAB复刻Liang (2004)¹中溃坝流算例的一个记录。
二维溃坝流（Dam Break）问题是浅水模型经典的一个测试算例，它测试了模型对急变流的模拟效果、以及对干-湿边界处理方法的有效性。相比于之前的模拟算例，本算例中需要重点处理的问题是：

模型的内边界的处理；
干-湿边界的处理。

本博客将着重解决第一个问题，而先不考虑第二个问题，即设置的模拟算例不含干-湿边界的处理。此外，本算例中涉及的控制方程与数值方法已经在《【浅水模型MATLAB】尝试完成一个数值模拟竞赛题》中介绍；不清楚的朋友可参考该博客内容。

如果你习惯用别的代码，也想做类似的建模尝试，十分欢迎一起交流！（最近有点想转python代码了，希望感兴趣的同志来一起交流）
如果各位朋友发现了文章或代码中的错误，亦或是改进之处，请不吝赐教，欢迎大家留言，一起改进模型！本博客文章将持续更新，上面也会标注提出改进建议的同志们。（不过，本人最近在忙活毕业论文，可能更新不及时）

同时，想要完整代码的朋友请联系我，我可无偿提供脚本文件。

希望同大家一起进步！

问题描述

本算例的计算区域为一个200m×200m的矩形平底水槽，水槽的四周都是垂直的固体壁面。如下图所示，计算域被分成x<100m和x>100m的左右两个部分；左侧初始水深为10m，右侧初始水深为5m。左右两个区域被两道平行于y方向的壁面阻隔，仅在95m<y<170m的区域联通。在模拟开始时，左侧水体会突然通过x=100m，95m<y<170m的区域向着右侧下泄，形成溃坝流。此外，模型中的所有壁面都是光滑的。
在这里插入图片描述

控制方程及数值方法

浅水方程及其数值计算方法

二维浅水方程的形式及其具体求解内容详见Liang的论文²和博客《【浅水模型MATLAB】尝试完成一个数值模拟竞赛题》。模型采用Godunov型有限体积法，通过一系列的处理，方程也保证了静水状态时压力与底坡源项的平衡。
此外，模型中的水力变量都定义在网格的中心位置。网格边界处的通量采用HLL求解器获得。

边界条件的实现

计算域的外边界均为无通量的free-slip闭合边界，边界处的法向速度和通量均被定义为0。在求解过程中，可将边界处的水力参数设置为其临近网格相同的物理量的值。
对于模型在x=100m处的内边界，模型需要定义其对应边界的通量为零。具体处理方式如下图所示。在对内边界左侧的相邻网格进行线性重构及通量计算时，需要通过一个辅助计算的虚网格，该虚网格有着和左侧相邻网格i相同的水力变量值。同理，在对内边界右侧的相邻网格进行线性重构及通量计算时，也需要通过一个辅助计算的虚网格，该虚网格有着和右侧相邻网格i相同的水力变量值。由于本模型采用了Minmod的限制器，所以此种处理会使得内边界对应的左侧变量U^L =U_i，而使右侧变量U^R =U_i+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. ↩︎