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

1、2022-12-251基于有限体积法的二维水流水质模拟及其可视化研究疯狂猫疯狂猫2022-12-252目录目录绪论浅水动力学的控制方程有限体积法原理浅水模型的程序实现数据可视化验证算例2022-12-253绪论绪论1，本课题的主要工作使用有限体积法（Finite Volume Methods）在无结构网格上离散二维浅水方程，求得水流速度场，并把求得的流场作为输入条件进一步求解对流扩散方程得出污染物浓度场。同时，计算过程将会产生大量的数据，本课题使用Delphi语言DirectX技术实现流场和浓度场的可视化。2022-12-254绪论绪论2，课题的现实意义本课题的浅水模型可以用来实现河流、河口、

2、湖泊、近海的水流水位模拟；而对流扩散方程是水质模型的基本控制方程，可以实现污染物（BOD，COD）、泥沙、热、盐的数值模拟，具有重大的现实意义。应用举例：河道流量及洪水预报、溃坝决堤、洪水漫滩、河口潮汐、盐水入侵、都市排水、海上油膜扩散等。2022-12-255第一章第一章 浅水动力学的控制方程浅水动力学的控制方程1，何谓浅水这里所说的浅水（Shallow Water）并非是指水深比较浅的水，而是具有它自身特定的含义。我们把满足以下四个方面的水称为浅水：有自由表面；有自由表面；以重力为主要驱动力，同时考虑风应力以重力为主要驱动力，同时考虑风应力 和地转柯氏力；和地转柯氏力；水平流速沿垂线近似均

3、匀分布；水平流速沿垂线近似均匀分布；水平运动尺度远大于垂直运动尺度；水平运动尺度远大于垂直运动尺度；满足浅水假设的包括：河流、湖泊、河口和近海。2022-12-256第一章第一章 浅水动力学的控制方程浅水动力学的控制方程2，一维浅水控制方程浅水方程又称圣维南方程（SaintVenant），是由圣维南1871年在法国科学院汇刊73卷中提出的，它的一维形式为：)()2(002fSSgxhgvxtvxhvth连续方程运动方程式中：h水位；v流速；g重力加速度；S0水底坡度；Sf摩阻坡度；x水平距离；t时间；2022-12-257第一章第一章 浅水动力学的控制方程浅水动力学的控制方程下面我们来推导以一

4、维浅水方程，考虑以下的一维单元模型：xxQQQ单元左侧进水流量为Q，右侧出水可以表示 。右减左，得到该单元上的净出流量为 ，在t时刻内的净出水量为：。同时考虑t时刻内，单元控制体内水量增加为 。根据单元水量平衡可以得到：令t0,x0,可以得到连续性方程。xxQQxxQtxxQxtth)(0txthtxxQ2022-12-258第一章第一章 浅水动力学的控制方程浅水动力学的控制方程按照同样的思路，我们可以推导出动量方程：xxQVQV)(QV在t时刻内，左右界面净输出动量为 ；左右界面水压差为 ；水重在流向的分量为hgS0 x;底摩阻力在流向的分量为hgSf x;控制体内的动量增量为：由控制体动量

5、平衡得出动量方程。txxQV)(xxgh/)2/(2txtQ)/(根据以上两个方程的推导，我们可以看到，其实浅水方根据以上两个方程的推导，我们可以看到，其实浅水方程就是物理守恒定律（水量和动量）在有限体积上的应程就是物理守恒定律（水量和动量）在有限体积上的应用。用。2022-12-2593，从一维到两维 我们上述的浅水方程是在只考虑横向尺度（X方向）的情况下得到的，对上述方程进行纵向尺度（Y方向）扩展，可以很容易地得到二维的浅水方程：第一章第一章 浅水动力学的控制方程浅水动力学的控制方程byghhvyhuvxthvbxhuvyghhuxthuyhvxhuth)2()()()()2()(0)()

6、(2222连续方程X向运动方程Y向运动方程2022-12-2510第一章第一章 浅水动力学的控制方程浅水动力学的控制方程或写成向量形式：式中：byGxFtqThvhuhq),(ThuvghhuhuF),2/,(22TghhvhuvhvG)2/,(22TfyyfxxSSghSSghb)(),(,0(002022-12-25114，求解浅水方程常用的数值方法 第一章第一章 浅水动力学的控制方程浅水动力学的控制方程有限差分（Finite Difference Methods）FDM以Taylor级数展开为工具，对微分方程中的导数项用差分项来逼近，通过求解差分方程得到近似解。优点：简单易学，处理效率高

7、。缺点：使用的矩形网格对边界逼近较差。有限元（Finite Element Methods）把计算域分为有限个单元，分单元对解逼近，使微分方程空间积分的加权残差极小化。FEM在数学上适于求解椭圆型方程的边值问题，不适合求解以对流为主的输运问题。有限体积法（Finite Volume Methods）2022-12-2512第一章第一章 浅水动力学的控制方程浅水动力学的控制方程有限体积法和有限元一样将计算域分成若干单元。在计算出每个控制体边界沿法向输入（出）的流量和动量通量之后，对每个控制体分别进行水量和动量平衡计算，得到计算时段末各控制体的平均水深和流速。因此，FVM正是对于推导原始微分方程所

8、用控制体途径的回归。本课题使用FVM进行浅水方程的求解。上述三种方法在具体应用中并没有绝对的好，和绝对的坏，根据具体问题的不同，可以用不同的数值方法。比如，DHI的二维浅水动力学软件Mike 21采用的就是差分中的ADIQuick格式，而著名的SMS（Surface Water Modeling System）使用的是有限元法，英国Wallingford的InfoWorks中的RS（River Simulation）模块使用的是最为传统的四点差分格式。2022-12-2513第二章第二章 有限体积法原理有限体积法原理1，有限体积法离散原理 有限体积法的基本原理是在被离散化了的计算区域上，计算出

9、通过每个控制体边界沿法向输入或输出的流量和动量通量后，对每个控制体分别进行水量和动量平衡计算，最终得到计算时段末各控制体的平均水深和流速(假设水力要素在各控制体内均匀分布)。考虑浅水方程的向量形式：byGxFtq在控制体 上进行积分，得到：2022-12-2514第二章第二章 有限体积法原理有限体积法原理 利用散度定理，对左边的第二项化为沿控制体边界的积分：S为控制体的边界，n为边界外法线方向的单位向量。对于有m个边的凸多边形而言，等于被积函数在控制体各边上的法向值与该边长度的乘积，因此又可以写为：2022-12-2515第二章第二章 有限体积法原理有限体积法原理其中，为每条边的数值通量令原始

10、的矢量形式的方程化为至此，问题归结到如何求解Fn(q)2022-12-2516第二章第二章 有限体积法原理有限体积法原理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)，大大简化计算并提高了效率2022-12-2517第二章第二章 有限体积法原理有限体积法原理这样，对于每个控制体我们

11、可以得到最终的离散方程：上式左边表示控制体内守恒变量在 t内的变化，右边第一项表示沿各边法向输出的通量之和，第二项表示控制体内源项(入流及外力)在 t内的作用;这反映了守恒物理量的守恒原理:守恒物理量在控制体内随时间变化等于各边法向数值通量时间变化量之和及源项的时间变化量。2022-12-2518第二章第二章 有限体积法原理有限体积法原理 到目前为止，有限体积法的计算归结为法向通量的计算。这是有有限体积法的核心限体积法的核心所在，法向通量的计算目前常用的途径是通过在沿外法向建立单元水力模型，并求解一维黎曼问题而得到。已知条件是相邻控制体形心的守恒变量qi和qj（已通过旋转变换为外法向和切向）。

12、i为进行FVM计算的控制体，j为第j边的相邻控制体。因浅水方程的旋转不变性，在局部坐标系中方程为 根据控制体内变量的分布确定界面中点两侧的状态变量：qL对应于左侧，qR对应于右侧。通常在左右两侧存在间断。法向数值通量近似取为法向方向的一维黎曼问题的解。3，控制体界面处的水力模型与黎曼问题bGFtq2022-12-25190Ftq第二章第二章 有限体积法原理有限体积法原理所谓黎曼问题（Riemann Problem），又称间断问题，是指如下的偏微分方程在初值为间断情况下的解间断初值条件：0,0,)0,(xqRxqLxq黎曼问题对于初学者较难理解，我们可以感性地理解为：t=0时刻，x0左侧水位为q

13、L，右侧为qR,黎曼问题即是要求解t0时刻在x0处的水位值。2022-12-2520第二章第二章 有限体积法原理有限体积法原理4，FVS格式求解黎曼问题 黎曼问题有许多种近似求法，例如通量分裂法(flux-vector spliting method)、通量差分裂法(flux-difference spliting method)，及Osher数值方法等。本课题采用FVS格式求解黎曼问题。在此我们对FVS格式本身的推导不做过多的阐述，格式最终形式为：nnnnnnnnnnnvcucucuhFvcuvucuucuhF)()(4)(2)(23422204244422FvucuuhFnnnn缓流：急流

14、：法向数值通量法向数值通量FLRFF2022-12-2521第二章第二章 有限体积法原理有限体积法原理5，边界条件 当控制体的某一边处于计算区域的边界上时，那么该边界上的数值通量的处理就相应的转化成了边界问题。在此种情况下，被分析单元的水力要素状态向量qL是已知的，而相邻单元的水力要素状态向量qR，是未知的，因而必须确定qR，才可以利用FVS格式来估计法向数值通量。本课题涉及到三种边界条件：固壁边界，给定水位边界和给定单宽流量边界。如下图所示单宽流量边界单宽流量边界水位边界水位边界固壁边界固壁边界2022-12-2522第二章第二章 有限体积法原理有限体积法原理(1)陆地边界(闭边界):uR=

15、uL;hR=hL;(2)给定水位时：(3)给定单宽流量时，联立求解两个方程qR=uR*hR这两个方程最终可以化为一个非线形方程，本课题采用牛顿下山法迭代求解。)(2hRhLguLuR)(2hRhLguLuR2022-12-2523第三章第三章 浅水模型的程序实现浅水模型的程序实现1，软件构架设计 一般的计算流体动力学（CFD，Computational Fluid Dynamics）软件，包括三个模块：预处理（预处理（Preprocessor）主要进行计算域的网格划分，边界条件的设定；模型求解器（模型求解器（Model Solver）模型方程的离散求解；后处理（后处理（PostProcesso

16、r）模型计算结果的数据可视化；2022-12-2524第三章第三章 浅水模型的程序实现浅水模型的程序实现预处理预处理模型求解器模型求解器后处理后处理网格计算结果 本课题的预处理采用EasyMesh（Free，并且有源代码）软件进行三角网格划分；模型求解采用Delphi语言读取网格文件，求解FVM离散后的浅水方程，得到每个体积单元上的水位，流速（X方向和Y方向）和浓度；后处理采用Delphi支持下的DirectX技术实现计算网格的显示，流场的绘制，浓度场云图的绘制，浓度场等值线的绘制。模型求解和后处理整合在一个程序中，在求得结果之后，可以立即进行数据的可视化，从而达到动态显示的效果。2022-1

17、2-2525第三章第三章 浅水模型的程序实现浅水模型的程序实现2，预处理网格的生成 EasyMesh使用所谓的Delaunay法生成三角形网格，具体的生成过程在此不作深究，我们只需明白两点：最终生成的三角形网格会覆盖整个区域，并且彼此之间没有重叠；任何一个三角形的外接圆中，不包含其他三角形的结点，这样的三角形均匀性最好；我们以一个正方形区域为例，阐述网格的生成过程：Step 1：确定区域边界坐标2022-12-2526第三章第三章 浅水模型的程序实现浅水模型的程序实现（0，0）（500，0）（500，500）（0，500）该矩形区域有四个边界坐标：（0，0），（500，0），（500，500）

18、和（0，500），我们在记事本中描述边界坐标如下：从左到右，依次代表：结点编号（从0开始），X坐标，Y坐标，三角形边长，边界标志。2022-12-2527第三章第三章 浅水模型的程序实现浅水模型的程序实现Step 2：定义区域边界把边的定义编号依次写入文件，如下：从左到右依次为：边的编号，边的起始结点编号，结束结点编号，边界标志。012301232022-12-2528Step 3 保存文件为*.d，利用EasyMesh生成网格文件*.n表示结点文件，*.e表示单元文件，*.s表示边文件最终的文件内容：最终的文件内容：+dxf表示同时生成*.dxf文件第三章第三章 浅水模型的程序实现浅水模型的

19、程序实现2022-12-2529第三章第三章 浅水模型的程序实现浅水模型的程序实现最终网格：最终网格：2022-12-2530第三章第三章 浅水模型的程序实现浅水模型的程序实现3，模型数据结构设计 从网格文件读入的数据包含大量的网格信息，必须定义适当的数据结构来描述，包括结点的数据结构，边的数据结构和单元数据结构；本课题使用Delphi中的Record型数据结构（相当于C语言中的Structure结构体）来描述结点，边和单元；所有的结点放在一个动态数组里，在读取网格文件的时候确定数组的大小，利用数组下标进行结点数据的访问；边和单元的结构定义与结点相似。2022-12-2531TNode/节点的

20、数据结构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

21、;V:Real;H:Real;Concentration:Real;Area:Real;AccumulateFlux:TFlux;End;2022-12-2532TSide /边的数据结构 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;第三章第三章 浅水模型的程

22、序实现浅水模型的程序实现2022-12-2533第三章第三章 浅水模型的程序实现浅水模型的程序实现4，模型求解过程的封装类的使用为了便于模块之间的协调，以及求解过程中数据的共享，我们把求解浅水方程的过程封装在一个TModelSolver类中，该类实现了从读入网格文件到输出计算结果的全部求解过程。Type TModelSolver=Class Private FL,FR,F:TFlux;FWaterLevel:Real;/开边界水位 FQPerWidth:Real;/开边界宽流量 Procedure CalcNodeq;/按照面积加权得到结点的矢量 Function Average_Outer(

23、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;/通量差分裂格式2022-12-2534Public Node:TNodeArray;Element:TElementArray;Si

24、de:TSideArray;Procedure CalcElementFlux();/计算每个单元通量之和 Procedure Iteration();/浅水方程的迭代计算 procedure ReadFromFile(NodePath,ElementPath,SidePath:String);/从网格文件*.n,*.e,*.s读取数据到结构体数组ode,Element Side中 Procedure CalcSideLength();/使用两点间距离公式，计算每一边的长度 Procedure CalcCosSin();/计算每一边的cos,sin值，以左单元(特指逆时针方向)为准，右单元的值

25、符号相反 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的初始值 Pro

26、cedure SetInitCondition_H();/设置水深H的初始值 Procedure SetInitCondition_C();/设置浓度C的初始值 Procedure SetBoundaryConditon(QPerWidth,WaterLevel:Real);/设置边界条件：开边界的水位和单宽流量 Procedure ExportSideFlux();/输出边的数值通量到文件：调试用 Procedure SearchNeibElement;/得到与结点相邻的单元End;第三章第三章 浅水模型的程序实现浅水模型的程序实现2022-12-25355，模型求解的流程第三章第三章 浅水

27、模型的程序实现浅水模型的程序实现1，读入网格文件，读入网格文件2，计算边长，外法线单元向量，单元面积，计算边长，外法线单元向量，单元面积3，寻求每一结点的相邻单元，寻求每一结点的相邻单元4，设置初始条件，设置初始条件2022-12-25365，设置边界条件，设置边界条件*.e,*.s,*.n6，计算边的数值通量，计算边的数值通量7，计算单元的累计数值通量，计算单元的累计数值通量8，求得该时刻形心处的，求得该时刻形心处的U，V，H和和C9，回到，回到6迭代计算迭代计算第三章第三章 浅水模型的程序实现浅水模型的程序实现2022-12-2537第四章第四章 数据可视化数据可视化1，数据可视化的需求分

28、析模型求解过程中，每个单元将会产生四个数据（X流速，Y流速，浓度，水位），而整个计算区域有成千上百个单元，并且单元数据是在不停的迭代中，所以在一定的时间间隔内，会产生大量的数据，要把握发展趋势，必须采用一定的数据可视化工具。本课题使用DirectX技术实现了流场和浓度场的可视化。流场的绘制流场的绘制等值云图的填充等值云图的填充等值线的生成等值线的生成与与TecPlot的接口的接口2022-12-2538第四章第四章 数据可视化数据可视化2，流场的绘制流场箭头的绘制需要四个点的坐标：起点，终点，两翼两个点。起点利用三角形形心处的坐标。终点采用X方向的速度分量U和Y方向的速度V进行矢量相加，乘以一

29、个速度坐标转化系数求得。同时利用几何关系求得箭头左右两翼的坐标，用直线对四个点分别相连。2022-12-2539第四章第四章 数据可视化数据可视化3，等值云图的绘制根据文献可知，目前等值云图的绘制有三种方法：一个单元内使用一种颜色进行填充，值越大，一个单元内使用一种颜色进行填充，值越大，颜色越深颜色越深扫描线算法：对单元内的每个象素点差值计扫描线算法：对单元内的每个象素点差值计算颜色值算颜色值改进的扫描线法：在三角形内定义若干四边改进的扫描线法：在三角形内定义若干四边形，每个四边形内的颜色值是一样的，对四形，每个四边形内的颜色值是一样的，对四边形进行统一的颜色填充边形进行统一的颜色填充2022

30、-12-2540第四章第四章 数据可视化数据可视化本课题采用微软未公开发布的API函数GradientFill（），进行三角形的渐变填充。该方法相当于扫描线算法，但是，经过比较，要比扫描线算法快3倍以上。使用该函数时需要确定三个结点上的颜色值，三角形内的每个点的颜色值在这三个结点颜色值之间。2022-12-2541第四章第四章 数据可视化数据可视化4，等值线的绘制三角网格上的等值线绘制，首先要求得边上的差值点，然后追踪含有差值点的边，确定差值点的顺序，最后把求得的差值点进行光滑处理，在精度要求不高的情况下，可以直接把差值点相连，或者也可用Bezier曲线，B样条曲线来拟合。5，与Tecplot的接口Tecplot是著名的数据可视化软件，利用它可以实现动态显示过程*.avi文件的输出。在了解Tecplot的文件结构之后，我们完全有可能作出与它的数据接口。2022-12-2542第五章第五章 验证算例验证算例1,溃坝算例-水流的验证2022-12-2543第五章第五章 验证算例验证算例2022-12-2544第五章第五章 验证算例验证算例2,均匀流场中污染物的迁移-水质的验证2022-12-2545第五章第五章 验证算例验证算例2022-12-2546第五章第五章 验证算例验证算例3,某河道的水流水质模拟 见程序演示2022-12-2547结尾结尾报告结束,多谢各位光临!