【CFD小工坊】浅水方程的离散及求解方法

【CFD小工坊】浅水方程的离散及求解方法

  • 前言
  • 基于有限体积法的方程离散
  • 界面通量与源项计算
  • 干-湿网格的处理
  • 数值离散的稳定性条件
  • 参考文献

前言

我们模型的控制方程,即浅水方程组的表达式如下:
∂ U ∂ t + ∂ E ( U ) ∂ x + ∂ G ( U ) ∂ y = S ( U ) U = ( h h u h v ) , E ( U ) = ( h u h u 2 + g h 2 2 h u v ) , G ( U ) = ( h v h u v h v 2 + g h 2 2 ) , S ( U ) = ( 0 g h ( S 0 x − S f x ) g h ( S 0 y − S f y ) ) S 0 x = − ∂ z b ∂ x , S 0 y = − ∂ z b ∂ y S f x = n 2 u u 2 + v 2 h − 4 / 3 , S f y = n 2 v u 2 + v 2 h − 4 / 3 \dfrac{\partial \bold{U}}{\partial t} + \dfrac{\partial \bold{E(U)}}{\partial x} + \dfrac{\partial \bold{G(U)}}{\partial y} = \bold{S(U)} \\[6pt] \bold{U} = \left( \begin{matrix} h \\ hu \\ hv \end{matrix} \right), \bold{E(U)} = \left( \begin{matrix} hu \\ hu^2+\dfrac{gh^2}{2} \\ huv \end{matrix} \right), \bold{G(U)} = \left( \begin{matrix} hv \\ huv \\ hv^2+\dfrac{gh^2}{2} \end{matrix} \right), \\[6pt] \bold{S(U)} = \left( \begin{matrix} 0 \\ gh(S_{0x} - S_{fx}) \\ gh(S_{0y} - S_{fy}) \end{matrix} \right) \\[6pt] S_{0x} = -\dfrac{\partial z_b}{\partial x}, S_{0y} = -\dfrac{\partial z_b}{\partial y} \\[6pt] S_{fx} = n^2 u \sqrt{u^2+v^2} h^{-4/3}, S_{fy} = n^2 v \sqrt{u^2+v^2} h^{-4/3} tU+xE(U)+yG(U)=S(U)U= hhuhv ,E(U)= huhu2+2gh2huv ,G(U)= hvhuvhv2+2gh2 ,S(U)= 0gh(S0xSfx)gh(S0ySfy) S0x=xzb,S0y=yzbSfx=n2uu2+v2 h4/3,Sfy=n2vu2+v2 h4/3
对此,我们首先将式中物理量离散于三角形网格中,后采用有限体积法将浅水方程改写成其离散形式。在离散的浅水方程中,水位、流速等物理量的计算可转化为对网格界面的通量项和网格内源项的计算。对于离散方程中的通量项,本模型将会采用基于Riemann近似解的Roe格式数值求解界面通量。

基于有限体积法的方程离散

对于我们的浅水模型,水深、流速等物理量被定义在网格中心,如下图所示。对于网格i与j,它们间的质量、动量传输是通过其网格边界Ei,j发生的(下图中加粗的网格边)。网格i和j之间边界的外法相量定义为n=(nx, ny)。
在这里插入图片描述
之后,对于控制体网格i,我们采用高斯散度定理将浅水方程沿着边界Ei,j线积分,得到:
∬ Ω ∂ U ∂ t d Ω = − ∫ S [ E ( U ) n x + G ( U ) n y ] d s + ∬ Ω S ( U ) d Ω \iint_{\Omega} \dfrac{\partial \bold{U}}{\partial t} d\Omega = -\int_{S} [\bold{E(U)}n_x + \bold{G(U)}n_y] ds + \iint_{\Omega} \bold{S(U)} d\Omega ΩtUdΩ=S[E(U)nx+G(U)ny]ds+ΩS(U)dΩ
式中,Ω表示网格i所对应的控制体,S表示边界Ei,j;dΩ和ds分别表示面积分和线积分的微元。

在我们的模型中,各个水力要素在网格体内被假定是均匀分布的。在对上式的时间导数项采用前差格式离散,得到数值解为:
U i n + 1 = U i n − Δ t Ω i ∑ j = 1 3 F n j L j + Δ t S ˉ F n = E ( U ) n x + G ( U ) n y \bold{U}_{i}^{n+1} = \bold{U}_{i}^{n} - \dfrac{\Delta t}{\Omega_i} \sum_{j=1}^{3} \bold{F}_{nj}L_j + \Delta t \bar{S} \\[6pt] \bold{F_n} = \bold{E(U)}n_x + \bold{G(U)}n_y Uin+1=UinΩiΔtj=13FnjLj+ΔtSˉFn=E(U)nx+G(U)ny
其中,Δt为时间步长,Ωi表示网格单元i的面积,Fnj表示网格i中第j条边的外通量,Lj表示网格i中第j条边的长度;S表示源项, S ˉ \bar{S} Sˉ表示源项S在网格i的体积分值。

注意,上式就是模型求解过程中最核心的表达式。它完成了从该时间步水力变量到下一步水力变量的更新计算。

界面通量与源项计算

接下来要解决的问题是通量项 F n \bold{F_n} Fn的计算。首先,我们将网格边通量计算转化为一个Riemann近似解问题。采用Roe格式来求解这个近似黎曼解:
F n = 1 2 ( E ( U ~ L ) + G ( U ~ R ) ) − ∑ k = 1 3 α ~ k ∣ λ k ∣ γ k ) \bold{F_n} = \dfrac{1}{2} (\bold{E(\tilde{U}_L)} + \bold{G(\tilde{U}_R))} - \sum_{k=1}^{3} \tilde{\alpha}^k |\lambda^k| \gamma^k) Fn=21(E(U~L)+G(U~R))k=13α~kλkγk)
式中, λ k \lambda^k λk为基于 Roe 平均的雅克比矩阵 J J J的特征值, γ k \gamma^k γk为特征值对应的特征向量; α k \alpha^k αk表示特征强度。各项表达如下:
{ λ 1 = u ~ − c ~ λ 2 = u ~ λ 3 = u ~ + c ~ , { γ 1 = ( 1 , u − c ~ , v ~ ) T γ 2 = ( 0 , 0 , 1 ) T γ 3 = ( 1 , u + c ~ , v ~ ) T , { α ~ 1 = 1 2 [ h R − h L − h ~ c ~ ( u R − u L ) ] α ~ 2 = h ~ ( v R − v L ) α ~ 3 = 1 2 [ h R − h L + h ~ c ~ ( u R − u L ) ] \left\{ \begin{aligned} \lambda^1 & = \tilde u-\tilde{c} \\ \lambda^2 & = \tilde u \\ \lambda^3 & = \tilde u+\tilde{c} \end{aligned} \right. , \left\{ \begin{aligned} \gamma^1 & = (1, u-\tilde{c}, \tilde{v})^T \\ \gamma^2 & = (0,0,1)^T \\ \gamma^3 & = (1, u+\tilde{c}, \tilde{v})^T \end{aligned} \right. , \left\{ \begin{aligned} \tilde \alpha^1 & = \dfrac{1}{2}[h_R - h_L - \dfrac{\tilde h}{\tilde{c}}(u_R- u_L)] \\ \tilde \alpha^2 & = \tilde h (v_R - v_L) \\ \tilde \alpha^3 & = \dfrac{1}{2}[h_R - h_L + \dfrac{\tilde h}{\tilde{c}}(u_R- u_L)] \end{aligned} \right. λ1λ2λ3=u~c~=u~=u~+c~, γ1γ2γ3=(1,uc~,v~)T=(0,0,1)T=(1,u+c~,v~)T, α~1α~2α~3=21[hRhLc~h~(uRuL)]=h~(vRvL)=21[hRhL+c~h~(uRuL)]
上述变量中下标L和R分别表示该变量是边界内、外侧网格上的值(如上图所示); u ~ \tilde u u~ v ~ \tilde v v~ c ~ \tilde{c} c~ h ~ \tilde h h~都是Roe平均变量,定义如下:
u ~ = u L h L + u R h R h L + h R , v ~ = v L h L + v R h R h L + h R , c ~ = g ( h L + h R ) 2 , h ~ = h L h R \tilde{u} = \dfrac{u_L \sqrt{h_L} + u_R \sqrt{h_R}}{\sqrt{h_L} + \sqrt{h_R}}, \tilde{v} = \dfrac{v_L \sqrt{h_L} + v_R \sqrt{h_R}}{\sqrt{h_L} + \sqrt{h_R}}, \\[6pt] \tilde c = \sqrt{\dfrac{g(h_L + h_R)}{2}}, \tilde h=\sqrt{h_L h_R} u~=hL +hR uLhL +uRhR ,v~=hL +hR vLhL +vRhR ,c~=2g(hL+hR) ,h~=hLhR
注意,此处下标含L和R的速度 u ~ \tilde u u~均是垂直于网格边(与外法向n同向)的速度分量,且波速c的方向也与外法向n同向。同理,下标含L和R的 v ~ \tilde v v~则是与法相量n相垂直。

此外,对源项S中的底坡和摩阻项的处理将直接关系到模型的计算精度和稳定性。为了获得和谐、守恒的结果,对底坡源项进行特征分解和迎风处理以适应Riemann求解格式。首先,对底坡源项在控制体Ω内积分得到 S ˉ 0 \bar{S}_0 Sˉ0
S ˉ 0 = ∑ j = 0 3 ∑ k = 0 3 [ 1 2 ( 1 − s i g n ( λ k ˉ ) ) ( β k r k ˉ ) L j ] j ( β 1 , β 2 , β 3 ) = ( − 1 2 c ~ Δ z b , 0 , 1 2 c ~ Δ z b ) , Δ z b = ( z b ) L − ( z b ) R r 1 ˉ = ( 1 u ~ + c ~ v ~ ) , r 2 ˉ = ( 0 0 c ~ ) , r 3 ˉ = ( 1 u ~ − c ~ v ~ ) , \bar{S}_0 = \sum_{j=0}^{3} \sum_{k=0}^{3} [\dfrac{1}{2} (1-sign(\bar{\lambda^k})) ({\beta^k} \bar{r^k}) L_{j}]^j \\[6pt] (\beta^1,\beta^2,\beta^3)= (-\dfrac{1}{2}\tilde{c}\Delta{z_b}, 0 , \dfrac{1}{2}\tilde{c}\Delta{z_b}), \Delta{z_b} = (z_b)_L - (z_b)_R \\[6pt] \bar{r^1}= \left( \begin{matrix} 1 \\ \tilde{u} + \tilde{c} \\ \tilde{v} \end{matrix} \right), \bar{r^2}= \left( \begin{matrix} 0 \\ 0 \\ \tilde{c} \end{matrix} \right), \bar{r^3}= \left( \begin{matrix} 1 \\ \tilde{u} - \tilde{c} \\ \tilde{v} \end{matrix} \right), Sˉ0=j=03k=03[21(1sign(λkˉ))(βkrkˉ)Lj]j(β1,β2,β3)=(21c~Δzb,0,21c~Δzb),Δzb=(zb)L(zb)Rr1ˉ= 1u~+c~v~ ,r2ˉ= 00c~ ,r3ˉ= 1u~c~v~ ,
式中,sign( )为符号函数。

为增加格式的稳定性,对摩阻源项进行半隐式离散,其表达式为:
S f = ( 1 − θ ) S f n + θ S f n + 1 = S f n + θ ∂ S f ∂ U ∂ U ∂ t Δ t \bold{S}_f = (1-\theta) \bold{S}_f^{n} + \theta \bold{S}_f^{n+1} = \bold{S}_f^{n} + \theta \dfrac{\partial \bold{S}_f}{\partial \bold{U}} \dfrac{\partial \bold{U}}{\partial t} \Delta t Sf=(1θ)Sfn+θSfn+1=Sfn+θUSftUΔt
θ=0时意味着完全显式,θ=1时意味着完全隐式;令 Q f = ∂ S f ∂ U \bold{Q}_f = \dfrac{\partial \bold{S}_f}{\partial \bold{U}} Qf=USf,离散后的动量方程通量表达式为:
Δ U = U n + 1 − U n = [ I − Δ t θ Q f n ] − 1 [ Δ t Ω i ∑ j = 1 3 ( F n j + ∑ k = 0 3 [ 1 2 ( 1 − s i g n ( λ k ˉ ) ) ( β k r k ˉ ) ) L j + Δ t S f ] \Delta \bold{U} = \bold{U}^{n+1} - \bold{U}^{n} = [\bold{I}- \Delta t \theta \bold{Q}_f ^{n}]^{-1} [\dfrac{\Delta t}{\Omega_i} \sum_{j=1}^{3} (\bold{F}_{nj} +\sum_{k=0}^{3} [\dfrac{1}{2} (1-sign(\bar{\lambda^k})) ({\beta^k} \bar{r^k}) )L_j + \Delta t \bold{S}_f] ΔU=Un+1Un=[IΔQfn]1[ΩiΔtj=13(Fnj+k=03[21(1sign(λkˉ))(βkrkˉ))Lj+ΔtSf]

干-湿网格的处理

对于复杂的地形,计算区域内的实际模拟范围和计算边界变化频繁。网格上可能出现淹没-干底-再淹没的过程,这对模型的计算提出较高的稳定性要求。为避免水深较小时网格=出现流速过大的非物理现象,本模型采用限制水深法来准确、稳定地模拟网格淹没、干底这一过程。
首先,设置两个临界水深 h d r y h_{dry} hdry h w e t h_{wet} hwet(且 h d r y < h w e t h_{dry}<h_{wet} hdry<hwet),再将网格分为三类:

  1. 当水深大于 h w e t h_{wet} hwet,该网格为“湿网格”,参与正常的模拟计算,同时求解质量和动量通量;
  2. 当水深小于 h d r y h_{dry} hdry,该网格为“干网格”,不参与当前时间步的计算;
  3. 当水深在 h d r y h_{dry} hdry h w e t h_{wet} hwet之间,该网格为“半干网格”,仅计算该时间步的质量通量,而没有计算动量通量。

数值离散的稳定性条件

本模型基于Riemann显式方法,即当前时刻的网格变量可直接由上一时刻的计算结果得到,无需经过迭代计算。但显式求解法受到CFL条件的限制,实际计算过程中的最大时间步长通常不为定值。
本模型将采用动态时间步长方法得到模拟的时间步长;同时需要用户输入一个可允许的最大步长 Δ t m a x \Delta t_{max} Δtmax,以及一个允许的最大Courant数 C m a x C_{max} Cmax( 0 < C < 1.0 0< C <1.0 0<C<1.0)。模型的时间步长按下式确定:
Δ t = m i n ( Δ t m a x , C m a x m i n ( R i u i 2 + v i 2 + g h i ) ) , i = 1 , 2 , 3 , . . . , N c \Delta t = min(\Delta t_{max}, C_{max} min(\dfrac{R_i}{\sqrt{u^2_i + v^2_i} + \sqrt{gh_i}})), i=1,2,3,...,N_c Δt=min(Δtmax,Cmaxmin(ui2+vi2 +ghi Ri)),i=1,2,3,...,Nc

参考文献

  1. 刘臻,史宏达,黄燕.一种基于Roe格式的有限体积法在二维溃坝问题中的应用[J].中国海洋大学学报:自然科学版, 2007, 37(2):5.
  2. [王志力,耿艳芬,金生.具有复杂计算域和地形的二维浅水流动数值模拟[J].水利学报, 2005, 36(4):6.

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

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

C++list模拟实现

list模拟实现 1.链表结点2.类模板基本框架3.构造4.插入普通迭代器实现4.1尾插4.2普通迭代器实现4.3对比list和vector的iterator4.4迭代器的价值4.5insert4.6尾插头插复用写法 5.删除erase5.1erase5.2尾删头删复用写法 6.析构emptysizeclear6.1clear6.2size6.3 empty6.4 析构 7.…

<C++>类和对象-下

目录 一、构造函数的初始化 1. 构造函数体赋值 2. 初始化列表 2.1 概念 2.2 隐式类型转换式构造 2.3 explicit关键字 二、static静态成员 1. 概念 2. 特性 三、友元 1. 友元函数 2.友元类 四、内部类 1. 概念 五、匿名对象 1. const引用匿名对象 2. 匿名对象的隐式类型转换 总…

C++左右值及引用

1 左值和右值 简单记法&#xff1a;能取地址的是左值&#xff0c;不能取地址的是右值 右值一般是常量 例&#xff1a; i 是右值&#xff0c;因为先把 i 赋值给临时变量&#xff0c;临时变量在1&#xff0c;而临时变量是将亡值&#xff0c;&i取地址会报错 i是左值&#xf…

获取网卡上的IP、网关及DNS信息,获取最佳路由,遍历路由表中的条目(附源码)

VC常用功能开发汇总&#xff08;专栏文章列表&#xff0c;欢迎订阅&#xff0c;持续更新...&#xff09;https://blog.csdn.net/chenlycly/article/details/124272585C软件异常排查从入门到精通系列教程&#xff08;专栏文章列表&#xff0c;欢迎订阅&#xff0c;持续更新...&a…

Oracle高级

Oracle数据库高级技术与最佳实践 Oracle数据库是全球领先的关系型数据库管理系统&#xff0c;广泛用于企业级应用程序和数据管理。在这篇博客中&#xff0c;我们将深入探讨Oracle数据库的高级技术和最佳实践&#xff0c;帮助数据库管理员和开发人员更好地利用其强大功能。 1. …

Linux进程控制

文章目录 前言一、进程创建1、fork函数2、写时拷贝3、子进程从哪里开始执行父进程代码 二、进程终止1、进程终止时&#xff0c;操作系统做了什么2、进程终止的常见方式2.1 main函数退出码 3、在代码中终止进程3.1 使用return语句终止进程3.2 使用exit函数终止进程3.3 使用_exit…

c#设计模式-结构型模式 之 组合模式

&#x1f680;简介 组合模式又名部分整体模式&#xff0c;是一种 结构型设计模式 &#xff0c;是用于把一组相似的对象当作一个 单一的对象 。组合模式 依据树形结构来组合对象 &#xff0c;用来表示部分以及整体层&#xff0c;它可以让你将对象组合成树形结构&#xff0c;并且…

leetCode 45.跳跃游戏 II 贪心算法

45. 跳跃游戏 II - 力扣&#xff08;LeetCode&#xff09; 给定一个长度为 n 的 0 索引整数数组 nums。初始位置为 nums[0]。 每个元素 nums[i] 表示从索引 i 向前跳转的最大长度。换句话说&#xff0c;如果你在 nums[i] 处&#xff0c;你可以跳转到任意 nums[i j] 处: 0 &…

大语言模型之十四-PEFT的LoRA

在《大语言模型之七- Llama-2单GPU微调SFT》和《大语言模型之十三 LLama2中文推理》中我们都提到了LoRA&#xff08;低秩分解&#xff09;方法&#xff0c;之所以用低秩分解进行参数的优化的原因是为了减少计算资源。 我们以《大语言模型之四-LlaMA-2从模型到应用》一文中的图…

已解决 Bug——IndexError: index 3 is out of bounds for axis 0 with size 3问题

&#x1f337;&#x1f341; 博主猫头虎&#xff08;&#x1f405;&#x1f43e;&#xff09;带您 Go to New World✨&#x1f341; &#x1f984; 博客首页: &#x1f405;&#x1f43e;猫头虎的博客&#x1f390;《面试题大全专栏》 &#x1f995; 文章图文并茂&#x1f996…

自然语言处理状况简介

一、说明 自然语言处理已经进入大模型时代&#xff0c;然而从业人员必须了解整个知识体系、发展过程、知识结构&#xff0c;应用范围等一系列知识。本篇将报道此类概况。 二、自然语言处理简介 自然语言处理&#xff0c;或简称NLP&#xff0c;是处理和转换文本的计算机科学学科…

格点数据可视化(美国站点的日降雨数据)

获取美国站点的日降雨量的格点数据&#xff0c;并且可视化 导入模块 from datetime import datetime, timedelta from urllib.request import urlopenimport cartopy.crs as ccrs import cartopy.feature as cfeature import matplotlib.colors as mcolors import matplotli…

JAVA学习(2)-全网最详细~

&#x1f308;write in front&#x1f308; &#x1f9f8;大家好&#xff0c;我是Aileen&#x1f9f8;.希望你看完之后&#xff0c;能对你有所帮助&#xff0c;不足请指正&#xff01;共同学习交流. &#x1f194;本文由Aileen_0v0&#x1f9f8; 原创 CSDN首发&#x1f412; 如…

选择排序算法:简单但有效的排序方法

在计算机科学中&#xff0c;排序算法是基础且重要的主题之一。选择排序&#xff08;Selection Sort&#xff09;是其中一个简单但非常有用的排序算法。本文将详细介绍选择排序的原理和步骤&#xff0c;并提供Java语言的实现示例。 选择排序的原理 选择排序的核心思想是不断地从…

CCF CSP认证 历年题目自练 Day20

题目一 试题编号&#xff1a; 201903-1 试题名称&#xff1a; 小中大 时间限制&#xff1a; 1.0s 内存限制&#xff1a; 512.0MB 问题描述&#xff1a; 题目分析&#xff08;个人理解&#xff09; 常规题目&#xff0c;先看输入&#xff0c;第一行输入n表示有多少数字&am…

axb_2019_brop64

axb_2019_brop64 Arch: amd64-64-little RELRO: Partial RELRO Stack: No canary found NX: NX enabled PIE: No PIE (0x400000)64位&#xff0c;只开了NX __int64 repeater() {size_t v1; // raxchar s[208]; // [rsp0h] [rbp-D0h] BYREFprintf("…

认识PostgreSQL

深入认识PostgreSQL&#xff1a;开源世界的强大数据库 在当今数字化时代&#xff0c;数据是组织的最宝贵资源之一。数据库管理系统&#xff08;DBMS&#xff09;扮演着关键角色&#xff0c;帮助企业存储、管理和分析数据。PostgreSQL&#xff0c;作为一款开源的高级关系型数据库…

在vue3中使用vite-svg-loader插件

vite-svg-loader插件可以让我们像使用vue组件那样使用svg图&#xff0c;使用起来超级方便。 安装 npm install vite-svg-loader --save-dev使用 import svgLoader from vite-svg-loaderexport default defineConfig({plugins: [vue(), svgLoader()] })组件里使用 在路径后加…

vue-引入使用main.js全局常量

common.js 命名什么都可以&#xff0c;用来定义常量的 定义了之后使用export让此暴露出去 const QRaddress http://localhost:9875export{QRaddress, } main.js //引入刚刚的js import {QRaddress} from /config/common.js挂载 Vue.prototype.$QRaddress QRaddress使用 …

Python---pyecharts地图案例

基础&#xff1a; # 可视化图表&#xff0c;导包&#xff08;地图&#xff09; from pyecharts.charts import Map from pyecharts.options import VisualMapOpts# 生成地图对象 map Map()# 准备数据--------------元组 data [("北京市", 99),("上海市"…