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

news2024/11/24 3:23:51

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

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

更新于2024年9月17日

前言

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

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

本博客将着重解决第一个问题，而先不考虑第二个问题，即设置的模拟算例不含干-湿边界的处理。此外，本算例中涉及的控制方程与数值方法已经在《【浅水模型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。
在这里插入图片描述

代码框架与关键代码

我的模型代码主要分为参数设置、网格构建、初始化、主循环和其余函数等五个部分。

  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条件计算得到。

  1. 网格构建：网格有两个要素需要定义，一是网格的四个节点（Xp和Yp），二是网格的中心点（Xc和Yc）；网格中心点也即水力物理量定义的位置。代码略。
  2. 初始化：设置底高程zb=0，计算zb的梯度zbx和zby；设置左右区域的初始水位，之后再设置流速u、v为零。
  3. 主循环：（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
  1. 其余子函数：包括minmod限制器、HLL求解器等。代码略。

模拟结果

1.水位结果
在这里插入图片描述
2.流速结果（颜色表示合速度的大小，箭头表示速度方向）
在这里插入图片描述
3. 三维水面图
在这里插入图片描述

  1. 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. ↩︎

  2. 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. ↩︎

本文来自互联网用户投稿，该文观点仅代表作者本人，不代表本站立场。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如若转载，请注明出处：http://www.coloradmin.cn/o/2142962.html

如若内容造成侵权/违法违规/事实不符，请联系多彩编程网进行投诉反馈，一经查实，立即删除！

相关文章

【FreeRL】Rainbow_DQN的实现和测试

【FreeRL】Rainbow_DQN的实现和测试

文章目录 前言环境1 PER note2 C51 note3 Noisy note4 Rainbow note其他 前言 具体代码实现见&#xff1a;https://github.com/wild-firefox/FreeRL/blob/main/DQN_file/DQN_with_tricks.py 将其中所有的trick都用上即为Rainbow_DQN。 效果如下&#xff1a;&#xff08;学习曲…
阅读更多...
word文档的写入（1）

word文档的写入（1）

Word文档的写入 我们手动复制Excel信息&#xff0c;再粘贴进Word&#xff0c;进行文件保存的整个操作。属于机械性的重复劳动&#xff0c;并不能带来太大价值。在Excel和Word的操作内&#xff0c;也没有能很好解决此类问题的方法。如果遇到信息一多&#xff0c;几十上百个文件&…
阅读更多...
Win11小技巧之调节音量

Win11小技巧之调节音量

无意中发现&#xff0c;鼠标悬停在喇叭&#x1f508;处可通过滚轮调节音量&#xff0c;无需每次都点开音量面板&#xff0c;再悬停在音量滚动条处通过滚轮调节&#xff01;&#xff08;设计师……怎么不早告诉我……&#xff09; 不用点开&#xff0c;之前一直都是这么调节音量…
阅读更多...
c++—多态【万字】【多态的原理】【重写的深入学习】【各种继承关系下的虚表查看】

c++—多态【万字】【多态的原理】【重写的深入学习】【各种继承关系下的虚表查看】

目录 C—多态1.多态的概念2.多态的定义及实现2.1多态的构成条件2.2虚函数的重写2.2.1虚函数重写的两个例外&#xff1a;2.2.1.1协变2.2.1.2析构函数的重写 2.3 c11的override和final2.3.1final2.3.2override 2.4 重载、重写、重定义的对比 3.抽象类3.1抽象类的概念3.2接口继承和…
阅读更多...
5款录屏软件电脑版，哪一款更适合你？

5款录屏软件电脑版，哪一款更适合你？

身边不少做行政的小伙伴&#xff0c;经常需要制作一些培训视频、会议记录或是演示文稿。这就要求他们必须掌握一款好用的录屏软件。作为一个经常搜索各种办公软件的人&#xff0c;今天&#xff0c;我就来分享一下我使用过的五款录屏软件在录制电脑屏幕时的表现。 1、福昕录屏大…
阅读更多...
枚举类题目练习心得

枚举类题目练习心得

两数之和 题目如下&#xff1a; 一点思路&#xff1a;该题目仅限于数据量少的情况使用枚举&#xff0c;从题目分析来看&#xff0c;需求是给定一个数字&#xff0c;要求在给定数组中找到两个数字并使这两个数字和为给定数字且返回目标数字下标。参考题解思路结合本身思路代码…
阅读更多...
Leetcode—环形链表||

Leetcode—环形链表||

题目描述 思路 快慢指针 结论 我们需要用到一个重要的结论&#xff1a;让一个指针从链表起始位置开始遍历链表,同时让一个指针从判环时相遇点的位置开始绕环运行,两个指针都是每次均走一步,最终肯定会在入口点的位置相遇。 画图解释 1.利用快慢指针找到相遇点 2. 定义两个…
阅读更多...
java138-异常处理_java 138错误

java138-异常处理_java 138错误

//异常 public class test79 { //定义方法声明定义异常&#xff0c;在满足条件时抛出异常对象&#xff0c;程序转向异常处理 public double count(double n,double m)throws Exception { if (m 0) {//如果除数等于0.则抛出异常实例 throw new Ex…
阅读更多...
day03 - Java集合和常用类

day03 - Java集合和常用类

第一章 Collection集合 1. Collection概述 集合&#xff1a;java中提供的一种容器&#xff0c;可以用来存储多个数据 集合和数组既然都是容器&#xff0c;它们有啥区别呢&#xff1f; 数组的长度是固定的。集合的长度是不固定的。集合可以随时增加元素&#xff0c;其大小也随…
阅读更多...
kubeadm方式安装k8s+基础命令的使用

kubeadm方式安装k8s+基础命令的使用

一、安装环境 二、前期准备 1.设置免密登录 [rootk8s-master ~]# ssh-keygen [rootk8s-master ~]# ssh-copy-id root192.168.2.77 [rootk8s-master ~]# ssh-copy-id root192.168.2.88 2.yum源配置 3.清空创建缓存 4.主机映射&#xff08;三台主机都要设置&#xff09; 5.安装…
阅读更多...
vivado中选中bd文件后generate output product是什么用，create HDL wrapper是什么用

vivado中选中bd文件后generate output product是什么用，create HDL wrapper是什么用

vivado中选中bd文件后generate output product是什么用 在Vivado中&#xff0c;“Generate Output Products” 是一个重要的步骤&#xff0c;它用于生成IP核的输出产品&#xff0c;这些产品是将IP核集成到设计中所需的文件。这些输出产品包括&#xff1a; 综合文件&#xff…
阅读更多...
多线程下的共享变量访问数据竞争的问题

多线程下的共享变量访问数据竞争的问题

多线程下对共享变量的写存在数据竞问题可导致数据与预期不一致。最近在研究race conditions漏洞&#xff0c;用以下python 代码记录一下&#xff0c;以此论证&#xff0c;如下&#xff1a; from concurrent.futures import ThreadPoolExecutor globalNum 5 def test():global…
阅读更多...
微积分-反函数6.1（反函数）

微积分-反函数6.1（反函数）

表1提供了一项实验的数据&#xff0c;其中细菌培养物在有限营养基中以100个细菌开始&#xff1b;在定时记录下细菌数量随时间的变化。细菌数量 N N N 是时间 t t t 的函数&#xff1a; N f ( t ) N f(t) Nf(t)。 然而&#xff0c;假设生物学家改变了她的观点&#xff0c;开…
阅读更多...
京东App秒级百G日志传输存储架构设计与实战

京东App秒级百G日志传输存储架构设计与实战

本文作者&#xff1a;平台业务研发部-武伟峰&#xff0c;数据与智能部-李阳 背景 在日常工作中&#xff0c;我们通常需要存储一些日志&#xff0c;譬如用户请求的出入参、系统运行时打印的一些info、error之类的日志&#xff0c;从而对系统在运行时出现的问题有排查的依据。 …
阅读更多...
作为研发部门的负责人，如何助力产品在市场竞争中胜出？浅谈 CTQ

作为研发部门的负责人，如何助力产品在市场竞争中胜出？浅谈 CTQ

在激烈的市场竞争中&#xff0c;产品研发团队如何帮助企业的产品脱颖而出&#xff1f;成功的产品往往不仅依赖于强大的功能和技术创新&#xff0c;还需要通过高效的研发效能&#xff0c;包括效率、质量和创新&#xff0c;来提升产品的市场竞争力。在本文中&#xff0c;我们将探…
阅读更多...
文档内容识别系统源码分享

文档内容识别系统源码分享

文档内容识别检测系统源码分享 [一条龙教学YOLOV8标注好的数据集一键训练_70全套改进创新点发刊_Web前端展示] 1.研究背景与意义 项目参考AAAI Association for the Advancement of Artificial Intelligence 项目来源AACV Association for the Advancement of Computer Vis…
阅读更多...
一款源码阅读的插件

一款源码阅读的插件

文章目录 进度汇报功能预览添加高亮标记高亮风格设置笔记颜色设置数据概览高亮数据详情 结尾 进度汇报 之前提到最近有在开发一个源码阅读的IDEA插件&#xff0c;第一版已经开发完上传插件市场了&#xff0c;等官方审批通过就可以尝鲜了。插件名称&#xff1a;Mark source cod…
阅读更多...
基于STM32F407ZGT6——看门狗

基于STM32F407ZGT6——看门狗

独立看门狗 独立看门狗的时钟由独立的RC 振荡器LSI 提供&#xff0c;即使主时钟发生故障它仍然有效&#xff0c;非常独立。 LSI 的频率一般在30~60KHZ 之间&#xff0c;根据温度和工作场合会有一定的漂移&#xff0c; 所以独立看门狗的定时时间并不一定非常精确&#xff0c;只适…
阅读更多...
格式化u盘选择FAT还是NTFS U盘和硬盘格式化两者选谁

格式化u盘选择FAT还是NTFS U盘和硬盘格式化两者选谁

Mac用户在将U盘或硬盘进行格式化时&#xff0c;选择FAT还是NTFS往往是一个让人纠结的问题。很多用户不知道这两个格式之间有什么区别&#xff0c;更不知道在格式化时如何做出选择。本文将为大家介绍Mac选择FAT还是NTFS&#xff0c;并为大家推荐U盘和硬盘格式化两者选谁。 一、…
阅读更多...
36.贪心算法3

36.贪心算法3

1.坏了的计算器&#xff08;medium&#xff09; . - 力扣&#xff08;LeetCode&#xff09; 题目解析 算法原理 代码 class Solution {public int brokenCalc(int startValue, int target) {// 正难则反 贪⼼int ret 0;while (target > startValue) {if (target % 2 0…
阅读更多...
推荐文章
最新文章

Copyright @ 2022~2023