基于有限体积法的二维水流课件

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,*,基于有限体积法的二维水流水质模拟及其可视化研究,疯狂猫,2024/8/10,1,基于有限体积法的二维水流水质模拟及其可视化研究疯狂猫2023,目录,绪论,浅水动力学的控制方程,有限体积法原理,浅水模型的程序实现,数据可视化,验证算例,2024/8/10,2,目录绪论2023/8/212,绪论,1，本课题的主要工作,使用有限体积法（,Finite Volume Methods）,在无结构网格上离散二维浅水方程，求得水流速度场，并把求得的流场作为输入条件进一步求解对流扩散方程得出污染物浓度场。,同时，计算过程将会产生大量的数据，本课题使用,Delphi,语言,DirectX,技术实现流场和浓度场的可视化。,2024/8/10,3,绪论1，本课题的主要工作2023/8/213,绪论,2，课题的现实意义,本课题的浅水模型可以用来实现河流、河口、湖泊、近海的水流水位模拟；而对流扩散方程是水质模型的基本控制方程，可以实现污染物（,BOD，COD）、,泥沙、热、盐的数值模拟，具有重大的现实意义。,应用举例：河道流量及洪水预报、溃坝决堤、洪水漫滩、河口潮汐、盐水入侵、都市排水、海上油膜扩散等。,2024/8/10,4,绪论2，课题的现实意义2023/8/214,第一章 浅水动力学的控制方程,1，何谓浅水,这里所说的浅水（,Shallow Water）,并非是指水深比较浅的水，而是具有它自身特定的含义。我们把满足以下四个方面的水称为浅水：,有自由表面；,以重力为主要驱动力，同时考虑风应力,和地转柯氏力；,水平流速沿垂线近似均匀分布；,水平运动尺度远大于垂直运动尺度；,满足浅水假设的包括：河流、湖泊、河口和近海。,2024/8/10,5,第一章 浅水动力学的控制方程1，何谓浅水2023/8/215,第一章 浅水动力学的控制方程,2，一维浅水控制方程,浅水方程又称圣维南方程（,SaintVenant），,是由圣维南1871年在法国科学院汇刊73卷中提出的，它的一维形式为：,连续方程,运动方程,式中：,h,水位；,v,流速；,g,重力加速度；,S0,水底坡度；,Sf,摩阻坡度；,x,水平距离；,t,时间；,2024/8/10,6,第一章 浅水动力学的控制方程2，一维浅水控制方程,第一章 浅水动力学的控制方程,下面我们来推导以一维浅水方程，考虑以下的一维单元模型：,单元左侧进水流量为,Q，,右侧出水可以表示 。,右减左，得到该单元上的净出流量为 ，在,t,时刻内的净出水量为：。,同时考虑,t,时刻内，单元控制体内水量增加为 。,根据单元水量平衡可以得到：,令,t0,x,0,可以得到连续性方程。,2024/8/10,7,第一章 浅水动力学的控制方程下面我们来推导以一维浅水方程，考,第一章 浅水动力学的控制方程,按照同样的思路，我们可以推导出动量方程：,在,t,时刻内，左右界面净输出动量为 ；,左右界面水压差为 ；,水重在流向的分量为,hgS0,x;,底摩阻力在流向的分量为,hgSf,x;,控制体内的动量增量为：,由控制体动量平衡得出动量方程。,根据以上两个方程的推导，我们可以看到，其实浅水方程就是物理守恒定律（水量和动量）在有限体积上的应用。,2024/8/10,8,第一章 浅水动力学的控制方程按照同样的思路，我们可以推导出动,3，从一维到两维,我们上述的浅水方程是在只考虑横向尺度（,X,方向）的情况下得到的，对上述方程进行纵向尺度（,Y,方向）扩展，可以很容易地得到二维的浅水方程：,第一章 浅水动力学的控制方程,连续方程,X,向运动方程,Y,向运动方程,2024/8/10,9,3，从一维到两维第一章 浅水动力学的控制方程,第一章 浅水动力学的控制方程,或写成向量形式：,式中：,2024/8/10,10,第一章 浅水动力学的控制方程或写成向量形式：2023/8/2,4，求解浅水方程常用的数值方法,第一章 浅水动力学的控制方程,有限差分（,Finite Difference Methods）,FDM,以,Taylor,级数展开为工具，对微分方程中的导数项用差分项来逼近，通过求解差分方程得到近似解。优点：简单易学，处理效率高。缺点：使用的矩形网格对边界逼近较差。,有限元（,Finite Element Methods）,把计算域分为有限个单元，分单元对解逼近，使微分方程空间积分的加权残差极小化。,FEM,在数学上适于求解椭圆型方程的边值问题，不适合求解以对流为主的输运问题。,有限体积法（,Finite Volume Methods）,2024/8/10,11,4，求解浅水方程常用的数值方法第一章 浅水动力学的控制方程有,第一章 浅水动力学的控制方程,有限体积法和有限元一样将计算域分成若干单元。在计算出每个控制体边界沿法向输入（出）的流量和动量通量之后，对每个控制体分别进行水量和动量平衡计算，得到计算时段末各控制体的平均水深和流速。因此，,FVM,正是对于推导原始微分方程所用控制体途径的回归。,本课题使用,FVM,进行浅水方程的求解。,上述三种方法在具体应用中并没有绝对的好，和绝对的坏，根据具体问题的不同，可以用不同的数值方法。比如，,DHI,的二维浅水动力学软件,Mike 21,采用的就是差分中的,ADIQuick,格式，而著名的,SMS（Surface Water Modeling System）,使用的是有限元法，英国,Wallingford,的,InfoWorks,中的,RS（River Simulation）,模块使用的是最为传统的四点差分格式。,2024/8/10,12,第一章 浅水动力学的控制方程有限体积法和有限元一样将计算域,第二章 有限体积法原理,1，有限体积法离散原理,有限体积法的基本原理是在被离散化了的计算区域上，计算出通过每个控制体边界沿法向输入或输出的流量和动量通量后，对每个控制体分别进行水量和动量平衡计算，最终得到计算时段末各控制体的平均水深和流速(假设水力要素在各控制体内均匀分布)。,考虑浅水方程的向量形式：,在控制体 上进行积分，得到：,2024/8/10,13,第二章 有限体积法原理1，有限体积法离散原理 有限,第二章 有限体积法原理,利用散度定理，对左边的第二项化为沿控制体边界的积分：,S,为控制体的边界，,n,为边界外法线方向的单位向量。对于有,m,个边的凸多,边形而言，等于被积函数在控制体各边上的法向,值与该边长度的乘积，因此又可以写为：,2024/8/10,14,第二章 有限体积法原理 利用散度定理，,第二章 有限体积法原理,其中，为每条边的数值通量,令,原始的矢量形式的方程化为,至此，问题归结到如何求解,F,n,(q),2024/8/10,15,第二章 有限体积法原理其中，,第二章 有限体积法原理,2,从二维到一维,f(q),与,g(q),为,x,y,方向的通量。根据前人的研究，我们知道,f(q),与,g(q),具有旋转不变性，即,f(q),与,g(q),在法向上的投影，可以转换为先投影,q,到法向上，即满足关系,也既是,先将,q,投影到法向,n,得到,q，,即,q=T*q，,再将可代入,f,得,f(q);,再根据上式进一步得到,F(q)，,从而消去了,g(q)。,这样，将原来的二维问题转化成一维问题，即只需计算,f(q)，,大大简化计算并提高了效率,2024/8/10,16,第二章 有限体积法原理2,从二维到一维,第二章 有限体积法原理,这样，对于每个控制体我们可以得到最终的离散方程：,上式左边表示控制体内守恒变量在,t,内的变化，右边第一项表示沿各边法向输出的通量之和，第二项表示控制体内源项(入流及外力)在,t,内的作用;这反映了守恒物理量的守恒原理:守恒物理量在控制体内随时间变化等于各边法向数值通量时间变化量之和及源项的时间变化量。,2024/8/10,17,第二章 有限体积法原理这样，对于每个控制体我们可以得到最终的,第二章 有限体积法原理,到目前为止，有限体积法的计算归结为法向通量的计算。这是,有限体积法的核心,所在，法向通量的计算目前常用的途径是通过在沿外法向建立单元水力模型，并求解一维黎曼问题而得到。,已知条件是相邻控制体形心的守恒变量,qi,和,qj（,已通过旋转变换为外法向和切向）。,i,为进行,FVM,计算的控制体，,j,为第,j,边的相邻控制体。因浅水方程的旋转不变性，在局部坐标系中方程为,根据控制体内变量的分布确定界面中点两侧的状态变量：,qL,对应于左侧，,qR,对应于右侧。通常在左右两侧存在间断。法向数值通量近似取为法向方向的一维黎曼问题的解。,3，控制体界面处的水力模型与黎曼问题,2024/8/10,18,第二章 有限体积法原理 到目前为止，有,第二章 有限体积法原理,所谓黎曼问题（,Riemann Problem），,又称间断问题，是指如下的偏微分方程在初值为间断情况下的解,间断初值条件：,黎曼问题对于初学者较难理解，,我们可以感性地理解为：,t=0,时刻，,x0,左侧水位为,qL，,右侧为,qR,黎曼问题即是要求解,t0,时刻在,x0,处的水位值。,2024/8/10,19,第二章 有限体积法原理所谓黎曼问题（Riemann Prob,第二章 有限体积法原理,4，FVS格式求解黎曼问题,黎曼问题有许多种近似求法，例如通量分裂法(,flux-vector spliting method)、,通量差分裂法(,flux-difference spliting method)，,及,Osher,数值方法等。本课题采用,FVS,格式求解黎曼问题。在此我们对,FVS,格式本身的推导不做过多的阐述，格式最终形式为：,缓流：,急流：,法向数值通量,F,LR,F,F,2024/8/10,20,第二章 有限体积法原理4，FVS格式求解黎曼问题,第二章 有限体积法原理,5，边界条件,当控制体的某一边处于计算区域的边界上时，那么该边界上的数值通量的处理就相应的转化成了边界问题。在此种情况下，被分析单元的水力要素状态向量,qL,是已知的，而相邻单元的水力要素状态向量,qR，,是未知的，因而必须确定,qR，,才可以利用,FVS,格式来估计法向数值通量。本课题涉及到三种边界条件：固壁边界，给定水位边界和给定单宽流量边界。如下图所示,单宽流量边界,水位边界,固壁边界,2024/8/10,21,第二章 有限体积法原理5，边界条件,第二章 有限体积法原理,(1)陆地边界(闭边界):,uR=uL;hR=hL;,(2),给定水位时：,(3),给定单宽流量时，联立求解两个方程,qR=uR*hR,这两个方程最终可以化为一个非线形方程，本课题采用牛顿下山法迭代求解。,2024/8/10,22,第二章 有限体积法原理(1)陆地边界(闭边界):2023/8,第三章 浅水模型的程序实现,1，软件构架设计,一般的计算流体动力学（,CFD，Computational Fluid Dynamics）,软件，包括三个模块：,预处理（,Preprocessor）,主要进行计算域的网格划分，边界条件的设定；,模型求解器（,Model Solver）,模型方程的离散求解；,后处理（,PostProcessor）,模型计算结果的数据可视化；,2024/8/10,23,第三章 浅水模型的程序实现1，软件构架设计,第三章 浅水模型的程序实现,预处理,模型求解器,后处理,网格,计算结果,本课题的预处理采用,EasyMesh（Free，,并且有源代码）软件进行三角网格划分；模型求解采用,Delphi,语言读取网格文件，求解,FVM,离散后的浅水方程，得到每个体积单元上的水位，流速（,X,方向和,Y,方向）和浓度；后处理采用,Delphi,支持下的,DirectX,技术实现计算网格的显示，流场的绘制，浓度场云图的绘制，浓度场等值线的绘制。模型求解和后处理整合在一个程序中，在求得结果之后，可以立即进行数据的可视化，从而达到动态显示的效果。,2024/8/10,24,第三章 浅水模型的程序实现预处理模型求解器后处理网格计算结果,第三章 浅水模型的程序实现,2，预处理网格的生成,EasyMesh,使用所谓的,Delaunay,法生成三角形网格，具体的生成过程在此不作深究，我们只需明白两点：,最终生成的三角形网格会覆盖整个区域，并且彼此之间没有重叠；,任何一个三角形的外接圆中，不包含其他三角形的结点，这样的三角形均匀性最好；,我们以一个正方形区域为例，阐述网格的生成过程：,Step 1：,确定区域边界坐标,2024/8/10,25,第三章 浅水模型的程序实现2，预处理网格的生成,第三章 浅水模型的程序实现,（0，0）,（500，0）,（500，500）,（0，500）,该矩形区域有四个边界坐标：,（0，0），（500，0），（500，500）和（0，500），我们在记事本中描述边界坐标如下：,从左到右，依次代表：结点编号（从0开始），,X,坐标，,Y,坐标，三角形边长，边界标志。,2024/8/10,26,第三章 浅水模型的程序实现（0，0）（500，0）（500，,第三章 浅水模型的程序实现,Step 2：,定义区域边界,把边的定义编号依次写入文件，如下：,从左到右依次为：边的编号，边的起始结点编号，结束结点编号，边界标志。,0,1,2,3,0,1,2,3,2024/8/10,27,第三章 浅水模型的程序实现Step 2：定义区域边界从左到右,Step 3,保存文件为*.,d，,利用,EasyMesh,生成网格文件,*.,n,表示结点文件，*.,e,表示单元文件，*.,s,表示边文件,最终的文件内容：,+,dxf,表示同时生成*.,dxf,文件,第三章 浅水模型的程序实现,2024/8/10,28,Step 3 保存文件为*.d，利用EasyMesh生成网格,第三章 浅水模型的程序实现,最终网格：,2024/8/10,29,第三章 浅水模型的程序实现最终网格：2023/8/2129,第三章 浅水模型的程序实现,3，模型数据结构设计,从网格文件读入的数据包含大量的网格信息，必须定义适当的数据结构来描述，包括结点的数据结构，边的数据结构和单元数据结构；,本课题使用,Delphi,中的,Record,型数据结构（相当于,C,语言中的,Structure,结构体）来描述结点，边和单元；所有的结点放在一个动态数组里，在读取网格文件的时候确定数组的大小，利用数组下标进行结点数据的访问；边和单元的结构定义与结点相似。,2024/8/10,30,第三章 浅水模型的程序实现3，模型数据结构设计 从网,TNode/,节点的数据结构,TNode=Record,X:Real;/x,坐标,Y:Real;/y,坐标,q:Array0.3 of Real;,NeibElement:TIntArray;/,相邻单元,End;,第三章 浅水模型的程序实现,TElement /,单元数据结构,TElement=Record,Node1:Integer;,Node2:Integer;,Node3:Integer;,Element1:Integer;,Element2:Integer;,Element3:Integer;,Side:Array 0.2 of Integer;,XCenter:Real;,YCenter:Real;,U:Real;,V:Real;,H:Real;,Concentration:Real;,Area:Real;,AccumulateFlux:TFlux;,End;,2024/8/10,31,TNode/节点的数据结构第三章 浅水模型的程序实现,TSide /,边的数据结构,TSide=Record,StartPoint:Integer;,EndPoint:Integer;,LeftElement:Integer;,RightElement:Integer;,Mark:Integer;,SideLength:Real;,Cos:Real;,Sin:Real;,qL:Array0.3 OF Real;,qR:Array0.3 OF Real;,Flux:TFlux;,End;,第三章 浅水模型的程序实现,2024/8/10,32,TSide,第三章 浅水模型的程序实现,4，模型求解过程的封装类的使用,为了便于模块之间的协调，以及求解过程中数据的共享，我们,把求解浅水方程的过程封装在一个,TModelSolver,类中，该类实,现了从读入网格文件到输出计算结果的全部求解过程。,Type,TModelSolver=Class,Private,FL,FR,F:TFlux;,FWaterLevel:Real;/,开边界水位,FQPerWidth:Real;/,开边界宽流量,Procedure CalcNodeq;/,按照面积加权得到结点的矢量,Function Average_Outer(UL,VL,HL,CL,UR,VR,HR,CR:Real):TFlux;/,左右单元数值通量平均,Function Average_Inner(UL,VL,HL,CL,UR,VR,HR,CR:Real):TFlux;/,左右单元先平均变量，再求数值通量,Function FVS(UL,VL,HL,CL,UR,VR,HR,CR:Real):TFlux;/,通量向量分裂格式,Function FDS(UL,VL,HL,CL,UR,VR,HR,CR:Real):TFlux;/,通量差分裂格式,2024/8/10,33,第三章 浅水模型的程序实现4，模型求解过程的封装类的使用为,Public,Node:TNodeArray;,Element:TElementArray;,Side:TSideArray;,Procedure CalcElementFlux();/,计算每个单元通量之和,Procedure Iteration();/,浅水方程的迭代计算,procedure ReadFromFile(NodePath,ElementPath,SidePath:String);/,从网格文件*.,n,*.e,*.s,读取数据到结构体数组,ode,Element Side,中,Procedure CalcSideLength();/,使用两点间距离公式，计算每一边的长度,Procedure CalcCosSin();/,计算每一边的,cos,sin,值，以左单元(特指逆时针方向)为准，右单元的值符号相反,Procedure CalcElementArea();/,计算每一单元的面积,Procedure ExportToTecplot();/,输出数据文件到,Tecplot,进行可视化处理,Procedure CalcSideFlux();/,计算通过每条边的通量,Function NonLinearSolver(UL,HL,QR:Real):Real;/,使用,Newton,法求解非线性方程,Procedure SetInitCondition_U(U:Real);/,设置,X,方向流速,U,的初始值,Procedure SetInitCondition_V(V:Real);/,设置,y,方向流速,V,的初始值,Procedure SetInitCondition_H();/,设置水深,H,的初始值,Procedure SetInitCondition_C();/,设置浓度,C,的初始值,Procedure SetBoundaryConditon(QPerWidth,WaterLevel:Real);/,设置边界条件：开边界的水位和单宽流量,Procedure ExportSideFlux();/,输出边的数值通量到文件：调试用,Procedure SearchNeibElement;/,得到与结点相邻的单元,End;,第三章 浅水模型的程序实现,2024/8/10,34,Public第三章 浅水模型的程序实现2023/8/2134,5，模型求解的流程,第三章 浅水模型的程序实现,1，读入网格文件,2，计算边长，外法线单元向量，单元面积,3，寻求每一结点的相邻单元,4，设置初始条件,2024/8/10,35,5，模型求解的流程第三章 浅水模型的程序实现1，读入网格文件,5，设置边界条件*.,e,*.s,*.n,6，计算边的数值通量,7，计算单元的累计数值通量,8，求得该时刻形心处的,U，V，H,和,C,9，回到6迭代计算,第三章 浅水模型的程序实现,2024/8/10,36,5，设置边界条件*.e,*.s,*.n6，计算边的数值通量7,第四章 数据可视化,1，数据可视化的需求分析,模型求解过程中，每个单元将会产生四个数据（,X,流速，,Y,流速，浓度，水位），而整个计算区域有成千上百个单元，并且单元数据是在不停的迭代中，所以在一定的时间间隔内，会产生大量的数据，要把握发展趋势，必须采用一定的数据可视化工具。本课题使用,DirectX,技术实现了流场和浓度场的可视化。,流场的绘制,等值云图的填充,等值线的生成,与,TecPlot,的接口,2024/8/10,37,第四章 数据可视化1，数据可视化的需求分析模型求解过程中，每,第四章 数据可视化,2，流场的绘制,流场箭头的绘制需要四个点的坐标：起点，终点，两翼两个点。起点利用三角形形心处的坐标。,终点采用,X,方向的速度分量,U,和,Y,方向的速度,V,进行矢量相加，乘以一个速度坐标转化系数求得。,同时利用几何关系求得箭头左右两翼的坐标，用直线对四个点分别相连。,2024/8/10,38,第四章 数据可视化2，流场的绘制流场箭头的绘制需要四个点的坐,第四章 数据可视化,3，等值云图的绘制,根据文献可知，目前等值云图的绘制有三种方法：,一个单元内使用一种颜色进行填充，值越大，颜色越深,扫描线算法：对单元内的每个象素点差值计算颜色值,改进的扫描线法：在三角形内定义若干四边形，每个四边形内的颜色值是一样的，对四边形进行统一的颜色填充,2024/8/10,39,第四章 数据可视化3，等值云图的绘制根据文献可知，目前等值云,第四章 数据可视化,本课题采用微软未公开发布的,API,函数,GradientFill（），,进行三角形的渐变填充。该方法相当于扫描线算法，但是，经过比较，要比扫描线算法快3倍以上。使用该函数时需要确定三个结点上的颜色值，三角形内的每个点的颜色值在这三个结点颜色值之间。,2024/8/10,40,第四章 数据可视化本课题采用微软未公开发布的API函数Gra,第四章 数据可视化,4，等值线的绘制,三角网格上的等值线绘制，首先要求得边上的差值点，然后追踪含有差值点的边，确定差值点的顺序，最后把求得的差值点进行光滑处理，在精度要求不高的情况下，可以直接把差值点相连，或者也可用,Bezier,曲线，,B,样条曲线来拟合。,5，与,Tecplot,的接口,Tecplot,是著名的数据可视化软件，利用它可以实现动态显示过程*.,avi,文件的输出。在了解,Tecplot,的文件结构之后，我们完全有可能作出与它的数据接口。,2024/8/10,41,第四章 数据可视化4，等值线的绘制三角网格上的等值线绘制，首,第五章 验证算例,1,溃坝算例-水流的验证,2024/8/10,42,第五章 验证算例1,溃坝算例-水流的验证2023/8/214,第五章 验证算例,2024/8/10,43,第五章 验证算例2023/8/2143,第五章 验证算例,2,均匀流场中污染物的迁移-水质的验证,2024/8/10,44,第五章 验证算例2,均匀流场中污染物的迁移-水质的验证202,第五章 验证算例,2024/8/10,45,第五章 验证算例2023/8/2145,第五章 验证算例,3,某河道的水流水质模拟,见程序演示,2024/8/10,46,第五章 验证算例3,某河道的水流水质模拟,结尾,报告结束,多谢各位光临!,2024/8/10,47,结尾报告结束,多谢各位光临!2023/8/2147,