【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.rhkb.cn/news/147002.html

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

相关文章

【已解决】 Expected linebreaks to be ‘LF‘ but found ‘CRLF‘.

问题描述 团队都是用mac&#xff0c;只有我自己是windows&#xff0c;启动项目一直报错 Expected linebreaks to be ‘LF‘ but found ‘CRLF‘. 但我不能因为自己的问题去改团队配置&#xff0c;也尝试过该vscode配置默认是LF还是报错 思路 看文章vscode如何替换所有文件的…

信创办公–基于WPS的EXCEL最佳实践系列 (数据整理复制粘贴)

信创办公–基于WPS的EXCEL最佳实践系列 &#xff08;数据整理复制粘贴&#xff09; 目录 应用背景操作步骤1、数据查找与替换2、复制或粘贴数据3、使用自动填充工具4、将数据拆分到多列5、应用数字格式 应用背景 数据的整理复制粘贴等在日常的工作中经常使用。本章内容主要学习…

修改switch Nand无线区码 以支持高频5G 信道

环境&#xff1a;NS switch 问题&#xff1a;日版&#xff0c;港版无法连接大于44信道的5G WIFI 解决办法&#xff1a;修改PRODINFO.dec的WIFI 区域码 背景&#xff1a;我的switch是最早买的港版的一批&#xff0c;WIFI 只能连接日本的信道&#xff0c;家里的路由器是国行的&am…

黑豹程序员-CSS四种样式的定义方式及冲突后的就近原则

4种样式定义 浏览器缺省设置&#xff08;即默认&#xff09;外部样式表&#xff08;引用的css&#xff09;内部样式表&#xff08;位于 标签内部&#xff09;内联样式&#xff08;在 HTML 元素内部&#xff09; 样式冲突如何解决&#xff1f; 若多重样式将层叠为一个&#x…

CentOS安装OpenNebula(二)

被控端部署&#xff1a; 先要配置好yum源&#xff1a; [rootmaster yum.repos.d]# vim opennebula.repo[rootmaster yum.repos.d]# cat opennebula.repo [opennebula] nameopennebula baseurlhttps://downloads.opennebula.org/repo/5.6/CentOS/7/x86_64 enabled1 gpgkeyhttps…

WSL2安装历程

WLS2安装 1、系统检查 安装WSL2必须运行 Windows 10 版本 2004 及更高版本&#xff08;内部版本 19041 及更高版本&#xff09;或 Windows 11。 查看 Windows 版本及内部版本号&#xff0c;选择 Win R&#xff0c;然后键入winver。 2、家庭版升级企业版 下载HEU_KMS_Activ…

如何开发物联网 APP?

如何开发物联网 APP? 这个问题本身是不严谨的&#xff0c;APP只是手机端的一个控制或者用于显示的人机交互页面&#xff0c;物联网是通过传感器&#xff0c;物联网卡等模块把物体接入网络以方便远程监控或者控制等。 你问的应该是怎么开发出来一个远程控制物体的APP吧&#x…

数学建模Matlab之数据预处理方法

本文综合代码来自文章http://t.csdnimg.cn/P5zOD 异常值与缺失值处理 %% 数据修复 % 判断缺失值和异常值并修复&#xff0c;顺便光滑噪音&#xff0c;渡边笔记 clc,clear;close all; x 0:0.06:10; y sin(x)0.2*rand(size(x)); y(22:34) NaN; % 模拟缺失值 y(89:95) 50;% 模…

POJ 3109 Inner Vertices 离散化+树状数组

一、题目大意 围棋棋盘&#xff0c;如果某个坐标上下左右的四个方向都存在棋子&#xff0c;那么ans1&#xff0c;根据输入的棋子数量&#xff0c;求出ans的数量。 二、解题思路 题目中有说到如果程序不会结束&#xff0c;那么输出-1&#xff0c;这其实是无源之水&#xff0c…

Linux下基本指令(上)

文章内容&#xff1a; 1. ls 指令 语法&#xff1a; ls [选项][目录或文件] 功能&#xff1a;对于目录&#xff0c;该命令列出该目录下的所有子目录与文件。对于文件&#xff0c;将列出文件名以及其他信息。 单个ls显示当前目录下的文件和目录 常用选项&#…

Vue iconfont-阿里巴巴矢量图标库用法

一、vue使用 选择心仪的图标 加入购物车 点击右上角购物车&#xff0c;点击添加至项目 在资源管理 可以看到我的项目 进入项目设置勾选彩色 点击下载到本地 解压压缩包 在main.js文件内导入css文件 import "/assets/font_icon/iconfont.css"; 使用&#xff1a; 复…

Java开源工具库使用之Lombok

文章目录 前言一、常用注解1.1 AllArgsConstructor/NoArgsConstructor/RequiredArgsConstructor1.2 Builder1.3 Data1.4 EqualsAndHashCode1.5 Getter/Setter1.6 Slf4j/Log4j/Log4j2/Log1.7 ToString 二、踩坑2.1 Getter/Setter 方法名不一样2.2 Builder 不会生成无参构造方法2…

Ubuntu Qt 5.15.2 支持 aarch64

概述 AArch64是ARMv8 架构的一种执行状态。 为了更广泛地向企业领域推进&#xff0c;需要引入64 位构架。 同时也需要在ARMv8 架构中引入新的AArch64 执行状态。 AArch64 不是一个单纯的32 位ARM 构架扩展&#xff0c;而是ARMv8 内全新的构架&#xff0c;完全使用全新的A64 指令…

黑豹程序员-架构师学习路线图-百科:Git/Gitee(版本控制)

文章目录 1、什么是版本控制2、特点3、发展历史4、SVN和Git比较5、Git6、GitHub7、Gitee&#xff08;国产&#xff09;8、Git的基础命令 1、什么是版本控制 版本控制系统&#xff08; Version Control &#xff09;版本控制是一种管理和跟踪软件开发过程中的代码变化的系统。它…

【Python】time模块和datetime模块的部分函数说明

时间戳与日期 在说到这俩模块之前&#xff0c;首先先明确几个概念&#xff1a; 时间戳是个很单纯的东西&#xff0c;没有“时区”一说&#xff0c;因为时间戳本质上是经过的时间。日常生活中接触到的“日期”、“某点某时某分”准确的说是时间点&#xff0c;都是有时区概念的…

C#,数值计算——Ranq1的计算方法与源程序

1 文本格式 using System; namespace Legalsoft.Truffer { /// <summary> /// Recommended generator for everyday use.The period is 1.8E19. Calling /// conventions same as Ran, above. /// </summary> public class Ranq1 { …

2023版 STM32实战4 滴答定时器精准延时

SysTick简介与特性 -1- SysTick属于系统时钟。 -2- SysTick定时器被捆绑在NVIC中。 -3- SysTick可以产生中断,且中断不可屏蔽。 SysTick的时钟源查看 通过时钟树可以看出滴答的时钟最大为72MHZ/89MHZ 通过中文参考手册也可以得到这个结论 代码编写&#xff08;已经验证&a…

嵌入式Linux应用开发-基础知识-第十九章驱动程序基石②

嵌入式Linux应用开发-基础知识-第十九章驱动程序基石② 第十九章 驱动程序基石②19.3 异步通知19.3.1 适用场景19.3.2 使用流程19.3.3 驱动编程19.3.4 应用编程19.3.5 现场编程19.3.6 上机编程19.3.7 异步通知机制内核代码详解 19.4 阻塞与非阻塞19.4.1 应用编程19.4.2 驱动编程…

【算法分析与设计】回溯法(上)

目录 一、学习要点1.1 回溯法1.2 问题的解空间1.3 0-1背包问题的解空间1.4 旅行售货员问题的解空间1.5 生成问题状态的基本方法 二、回溯法的基本思想三、回溯算法的适用条件四、递归回溯五、迭代回溯六、子集树与排列树七、装载问题八、批处理作业调度问题 一、学习要点 理解回…

【数据结构与算法】通过双向链表和HashMap实现LRU缓存 详解

这个双向链表采用的是有伪头节点和伪尾节点的 与上一篇文章中单链表的实现不同&#xff0c;区别于在实例化这个链表时就初始化了的伪头节点和伪尾节点&#xff0c;并相互指向&#xff0c;在第一次添加节点时&#xff0c;不需要再考虑空指针指向问题了。 /*** 通过链表与HashMa…