做决定网站,最好看免费观看高清大全多多电影,为什么网站开发这么便宜,手机编程教程【浅水模型MATLAB】尝试复刻SCI论文中的溃坝流算例 前言问题描述控制方程及数值方法浅水方程及其数值计算方法边界条件的实现 代码框架与关键代码模拟结果 更新于2024年9月17日 前言
这篇博客算是学习浅水方程#xff0c;并利用MATLAB复刻Liang (2004)1中溃坝流算例的一个记录… 【浅水模型MATLAB】尝试复刻SCI论文中的溃坝流算例 前言问题描述控制方程及数值方法浅水方程及其数值计算方法边界条件的实现 代码框架与关键代码模拟结果 更新于2024年9月17日 前言
这篇博客算是学习浅水方程并利用MATLAB复刻Liang (2004)1中溃坝流算例的一个记录。 二维溃坝流Dam Break问题是浅水模型经典的一个测试算例它测试了模型对急变流的模拟效果、以及对干-湿边界处理方法的有效性。相比于之前的模拟算例本算例中需要重点处理的问题是
模型的内边界的处理干-湿边界的处理。
本博客将着重解决第一个问题而先不考虑第二个问题即设置的模拟算例不含干-湿边界的处理。此外本算例中涉及的控制方程与数值方法已经在《【浅水模型MATLAB】尝试完成一个数值模拟竞赛题》中介绍不清楚的朋友可参考该博客内容。
如果你习惯用别的代码也想做类似的建模尝试十分欢迎一起交流最近有点想转python代码了希望感兴趣的同志来一起交流 如果各位朋友发现了文章或代码中的错误亦或是改进之处请不吝赐教欢迎大家留言一起改进模型本博客文章将持续更新上面也会标注提出改进建议的同志们。不过本人最近在忙活毕业论文可能更新不及时
同时想要完整代码的朋友请联系我我可无偿提供脚本文件。
希望同大家一起进步
问题描述
本算例的计算区域为一个200m×200m的矩形平底水槽水槽的四周都是垂直的固体壁面。如下图所示计算域被分成x100m和x100m的左右两个部分左侧初始水深为10m右侧初始水深为5m。左右两个区域被两道平行于y方向的壁面阻隔仅在95my170m的区域联通。在模拟开始时左侧水体会突然通过x100m95my170m的区域向着右侧下泄形成溃坝流。此外模型中的所有壁面都是光滑的。
控制方程及数值方法
浅水方程及其数值计算方法
二维浅水方程的形式及其具体求解内容详见Liang的论文2和博客《【浅水模型MATLAB】尝试完成一个数值模拟竞赛题》。模型采用Godunov型有限体积法通过一系列的处理方程也保证了静水状态时压力与底坡源项的平衡。 此外模型中的水力变量都定义在网格的中心位置。网格边界处的通量采用HLL求解器获得。
边界条件的实现
计算域的外边界均为无通量的free-slip闭合边界边界处的法向速度和通量均被定义为0。在求解过程中可将边界处的水力参数设置为其临近网格相同的物理量的值。 对于模型在x100m处的内边界模型需要定义其对应边界的通量为零。具体处理方式如下图所示。在对内边界左侧的相邻网格进行线性重构及通量计算时需要通过一个辅助计算的虚网格该虚网格有着和左侧相邻网格i相同的水力变量值。同理在对内边界右侧的相邻网格进行线性重构及通量计算时也需要通过一个辅助计算的虚网格该虚网格有着和右侧 相邻网格i相同的水力变量值。由于本模型采用了Minmod的限制器所以此种处理会使得内边界对应的左侧变量UL Ui而使右侧变量UR Ui1。
代码框架与关键代码
我的模型代码主要分为参数设置、网格构建、初始化、主循环和其余函数等五个部分。
设置物理参数、网格参数、时间参数等。代码如下所示
grav 9.81; % Gravitational acceleration
rho 1000; % Density
CFL 0.5; % Courant NumberLx 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 valuedx 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的正方形底高程为zb00.0。最大允许的Courant数设置为0.5初始时间步为0.05s之后的每一个时间步通过CFL条件计算得到。
网格构建网格有两个要素需要定义一是网格的四个节点Xp和Yp二是网格的中心点Xc和Yc网格中心点也即水力物理量定义的位置。代码略。初始化设置底高程zb0计算zb的梯度zbx和zby设置左右区域的初始水位之后再设置流速u、v为零。主循环1计算网格边界处的水位、水深、流速值2设置内边界条件3计算通量项FL、FR、GL和GR4利用HLL求解通量F和G5计算源项S6计算新一个时间步的eta、h、u和v。除了上述步骤2和3其余计算过程与博客《【浅水模型MATLAB】尝试完成一个数值模拟竞赛题》中代码基本一致涉及的关键代码如下
while(ttend)% estimate the dtdtx 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 dtyetan eta; hn h;un u; vn v;% 2rd-order Runge-Kutta Methodfor 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 conditionsff find((yc95) (yc170));% 2.1 left cellsde 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/21) eta(ff,Nx/2) 0.5*dx*de(ff);clear de du dv% 2.2 right cellsde minmod((eta(:,Nx/22)-eta(:,Nx/21))/dx, ...(eta(:,Nx/21)-eta(:,Nx/21))/dx);du minmod((u(:,Nx/22)-u(:,Nx/21))/dx, ...(u(:,Nx/21)-u(:,Nx/21))/dx);dv minmod((v(:,Nx/22)-v(:,Nx/21))/dx, ...(v(:,Nx/21)-v(:,Nx/21))/dx);exR(ff,Nx/21) eta(ff,Nx/21) - 0.5*dx*de(ff);exL(ff,Nx/22) eta(ff,Nx/21) 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 boundaryff find((yc95) (yc170));F1(ff,Nx/21) 0; F3(ff,Nx/21) 0;F2_L(ff,1) 0.5*grav*(exL(ff,Nx/21).^2 - ...exL(ff,Nx/21).*(zbp(ff1,Nx/21)zbp(ff,Nx/21)));F2_R(ff,1) 0.5*grav*(exR(ff,Nx/21).^2 - ...exR(ff,Nx/21).*(zbp(ff1,Nx/21)zbp(ff,Nx/21)));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 boundaryff find((yc95) (yc170));hu(ff,Nx/2) h(ff,Nx/2).*u(ff,Nx/2) ...- dt/dx*(F2_L(ff) - F2(ff,Nx/2)) ...- dt/dy*(G2(ff1,Nx/2)-G2(ff,Nx/2)) dt*S2(ff,Nx/2);hu(ff,Nx/21) h(ff,Nx/21).*u(ff,Nx/21) ...- dt/dx*(F2(ff,Nx/22) - F2_R(ff)) ...- dt/dy*(G2(ff1,Nx/21)-G2(ff,Nx/21)) dt*S2(ff,Nx/21);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. ↩︎
