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

基于有限体积法的二维水流水质模拟及其可视化研究疯狂猫2022/11/221基于有限体积法的二维水流水质模拟及其可视化研究疯狂猫2022目录绪论浅水动力学的控制方程有限体积法原理浅水模型的程序实现数据可视化验证算例2022/11/222目录绪论2022/11/222绪论1，本课题的主要工作 使用有限体积法（FiniteVolumeMethods）在无结构网格上离散二维浅水方程，求得水流速度场，并把求得的流场作为输入条件进一步求解对流扩散方程得出污染物浓度场。 同时，计算过程将会产生大量的数据，本课题使用Delphi语言＋DirectX技术实现流场和浓度场的可视化。2022/11/223绪论1，本课题的主要工作2022/11/223绪论2，课题的现实意义 本课题的浅水模型可以用来实现河流、河口、湖泊、近海的水流水位模拟；而对流扩散方程是水质模型的基本控制方程，可以实现污染物（BOD，COD）、泥沙、热、盐的数值模拟 ，具有重大的现实意义。

应用举例：河道流量及洪水预报、溃坝决堤、洪水漫滩、河口潮汐、盐水入侵、都市排水、海上油膜扩散等。2022/11/224绪论2，课题的现实意义2022/11/224第一章浅水动力学的控制方程1，何谓浅水 这里所说的浅水（ShallowWater）并非是指水深比较浅的水，而是具有它自身特定的含义。我们把满足以下四个方面的水称为浅水：

有自由表面； 以重力为主要驱动力，同时考虑风应力和地转柯氏力； 水平流速沿垂线近似均匀分布； 水平运动尺度远大于垂直运动尺度； 满足浅水假设的包括：河流、湖泊、河口和近海。2022/11/225第一章浅水动力学的控制方程1，何谓浅水2022/11/22第一章浅水动力学的控制方程2，一维浅水控制方程

浅水方程又称圣维南方程（Saint－Venant），是由圣维南1871年在法国科学院汇刊73卷中提出的，它的一维形式为：

－－－－－－－－－－－－－连续方程－－－－运动方程式中：h－水位；v－流速；g－重力加速度；S0－水底坡度；Sf－摩阻坡度；x－水平距离；t－时间；2022/11/226第一章浅水动力学的控制方程2，一维浅水控制方程－－－－－－第一章浅水动力学的控制方程下面我们来推导以一维浅水方程，考虑以下的一维单元模型：单元左侧进水流量为Q，右侧出水可以表示。右减左，得到该单元上的净出流量为，在Δt时刻内的净出水量为：。同时考虑Δt时刻内，单元控制体内水量增加为。根据单元水量平衡可以得到：令Δt－>0,Δx－>0,可以得到连续性方程。2022/11/227第一章浅水动力学的控制方程下面我们来推导以一维浅水方程，考第一章浅水动力学的控制方程按照同样的思路，我们可以推导出动量方程：在Δt时刻内，左右界面净输出动量为； 左右界面水压差为； 水重在流向的分量为hgS0Δx;

底摩阻力在流向的分量为－hgSfΔx;控制体内的动量增量为：由控制体动量平衡得出动量方程。根据以上两个方程的推导，我们可以看到，其实浅水方程就是物理守恒定律（水量和动量）在有限体积上的应用。2022/11/228第一章浅水动力学的控制方程按照同样的思路，我们可以推导出动3，从一维到两维

我们上述的浅水方程是在只考虑横向尺度（X方向）的情况下得到的，对上述方程进行纵向尺度（Y方向）扩展，可以很容易地得到二维的浅水方程：第一章浅水动力学的控制方程－－－－－－－－－－－连续方程－－－－X向运动方程－－－－Y向运动方程2022/11/2293，从一维到两维第一章浅水动力学的控制方程－－－－－－－－第一章浅水动力学的控制方程或写成向量形式：式中：2022/11/2210第一章浅水动力学的控制方程或写成向量形式：2022/11/4，求解浅水方程常用的数值方法

第一章浅水动力学的控制方程有限差分（FiniteDifferenceMethods） FDM以Taylor级数展开为工具，对微分方程中的导数项用差分项来逼近，通过求解差分方程得到近似解。优点：简单易学，处理效率高。缺点：使用的矩形网格对边界逼近较差。有限元（FiniteElementMethods）

把计算域分为有限个单元，分单元对解逼近，使微分方程空间积分的加权残差极小化。FEM在数学上适于求解椭圆型方程的边值问题，不适合求解以对流为主的输运问题。有限体积法（FiniteVolumeMethods）2022/11/22114，求解浅水方程常用的数值方法第一章浅水动力学的控制方程有第一章浅水动力学的控制方程 有限体积法和有限元一样将计算域分成若干单元。在计算出每个控制体边界沿法向输入（出）的流量和动量通量之后，对每个控制体分别进行水量和动量平衡计算，得到计算时段末各控制体的平均水深和流速。因此，FVM正是对于推导原始微分方程所用控制体途径的回归。 本课题使用FVM进行浅水方程的求解。 上述三种方法在具体应用中并没有绝对的好，和绝对的坏，根据具体问题的不同，可以用不同的数值方法。比如，DHI的二维浅水动力学软件Mike21采用的就是差分中的ADI－Quick格式，而著名的SMS（SurfaceWaterModelingSystem）使用的是有限元法，英国Wallingford的InfoWorks中的RS（RiverSimulation）模块使用的是最为传统的四点差分格式。2022/11/2212第一章浅水动力学的控制方程 有限体积法和有限元一样将计算域第二章有限体积法原理1，有限体积法离散原理有限体积法的基本原理是在被离散化了的计算区域上，计算出通过每个控制体边界沿法向输入或输出的流量和动量通量后，对每个控制体分别进行水量和动量平衡计算，最终得到计算时段末各控制体的平均水深和流速(假设水力要素在各控制体内均匀分布)。 考虑浅水方程的向量形式：在控制体上进行积分，得到：2022/11/2213第二章有限体积法原理1，有限体积法离散原理有限第二章有限体积法原理利用散度定理，对左边的第二项化为沿控制体边界的积分：S为控制体的边界，n为边界外法线方向的单位向量。对于有m个边的凸多边形而言， 等于被积函数在控制体各边上的法向值与该边长度的乘积，因此又可以写为： 2022/11/2214第二章有限体积法原理利用散度定理，第二章有限体积法原理其中，为每条边的数值通量令原始的矢量形式的方程化为至此，问题归结到如何求解Fn(q)2022/11/2215第二章有限体积法原理其中，第二章有限体积法原理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/11/2216第二章有限体积法原理2,从二维到一维第二章有限体积法原理这样，对于每个控制体我们可以得到最终的离散方程：上式左边表示控制体内守恒变量在Δt内的变化，右边第一项表示沿各边法向输出的通量之和，第二项表示控制体内源项(入流及外力)在Δ

t内的作用;这反映了守恒物理量的守恒原理:守恒物理量在控制体内随时间变化等于各边法向数值通量时间变化量之和及源项的时间变化量。2022/11/2217第二章有限体积法原理这样，对于每个控制体我们可以得到最终的第二章有限体积法原理到目前为止，有限体积法的计算归结为法向通量的计算。这是有限体积法的核心所在，法向通量的计算目前常用的途径是通过在沿外法向建立单元水力模型，并求解一维黎曼问题而得到。 已知条件是相邻控制体形心的守恒变量qi和qj（已通过旋转变换为外法向和切向）。i为进行FVM计算的控制体，j为第j边的相邻控制体。因浅水方程的旋转不变性，在局部坐标系中方程为根据控制体内变量的分布确定界面中点两侧的状态变量：qL对应于左侧，qR对应于右侧。通常在左右两侧存在间断。法向数值通量近似取为法向方向的一维黎曼问题的解。3，控制体界面处的水力模型与黎曼问题2022/11/2218第二章有限体积法原理到目前为止，有第二章有限体积法原理所谓黎曼问题（RiemannProblem），又称间断问题，是指如下的偏微分方程在初值为间断情况下的解间断初值条件：黎曼问题对于初学者较难理解，我们可以感性地理解为：t=0时刻，x＝0左侧水位为qL，右侧为qR,黎曼问题即是要求解t>0时刻在x＝0处的水位值。2022/11/2219第二章有限体积法原理所谓黎曼问题（RiemannProb第二章有限体积法原理4，FVS格式求解黎曼问题黎曼问题有许多种近似求法，例如通量分裂法(flux-vectorsplitingmethod)、通量差分裂法(flux-differencesplitingmethod)，及Osher数值方法等。本课题采用FVS格式求解黎曼问题。在此我们对FVS格式本身的推导不做过多的阐述，格式最终形式为：缓流：急流：法向数值通量FLR＝F＋＋F－2022/11/2220第二章有限体积法原理4，FVS格式求解黎曼问题第二章有限体积法原理5，边界条件当控制体的某一边处于计算区域的边界上时，那么该边界上的数值通量的处理就相应的转化成了边界问题。在此种情况下，被分析单元的水力要素状态向量qL是已知的，而相邻单元的水力要素状态向量qR，是未知的，因而必须确定qR，才可以利用FVS格式来估计法向数值通量。本课题涉及到三种边界条件：固壁边界，给定水位边界和给定单宽流量边界。如下图所示单宽流量边界水位边界固壁边界2022/11/2221第二章有限体积法原理5，边界条件第二章有限体积法原理(1)陆地边界(闭边界):

uR=uL;hR=hL;(2)给定水位时：

(3)给定单宽流量时，联立求解两个方程

qR=uR*hR

这两个方程最终可以化为一个非线形方程，本课题采用牛顿下山法迭代求解。2022/11/2222第二章有限体积法原理(1)陆地边界(闭边界):2022/1第三章浅水模型的程序实现1，软件构架设计一般的计算流体动力学（CFD，ComputationalFluidDynamics）软件，包括三个模块：预处理（Preprocessor）

主要进行计算域的网格划分，边界条件的设定；模型求解器（ModelSolver）

模型方程的离散求解；后处理（PostProcessor）

模型计算结果的数据可视化；2022/11/2223第三章浅水模型的程序实现1，软件构架设计第三章浅水模型的程序实现预处理模型求解器后处理网格计算结果本课题的预处理采用EasyMesh（Free，并且有源代码）软件进行三角网格划分；模型求解采用Delphi语言读取网格文件，求解FVM离散后的浅水方程，得到每个体积单元上的水位，流速（X方向和Y方向）和浓度；后处理采用Delphi支持下的DirectX技术实现计算网格的显示，流场的绘制，浓度场云图的绘制，浓度场等值线的绘制。模型求解和后处理整合在一个程序中，在求得结果之后，可以立即进行数据的可视化，从而达到动态显示的效果。2022/11/2224第三章浅水模型的程序实现预处理模型求解器后处理网格计算结果第三章浅水模型的程序实现2，预处理－网格的生成EasyMesh使用所谓的Delaunay法生成三角形网格，具体的生成过程在此不作深究，我们只需明白两点： 最终生成的三角形网格会覆盖整个区域，并且彼此之间没有重叠；任何一个三角形的外接圆中，不包含其他三角形的结点，这样的三角形均匀性最好；我们以一个正方形区域为例，阐述网格的生成过程：Step1：确定区域边界坐标2022/11/2225第三章浅水模型的程序实现2，预处理－网格的生成第三章浅水模型的程序实现（0，0）（500，0）（500，500）（0，500）该矩形区域有四个边界坐标：（0，0），（500，0），（500，500）和（0，500），我们在记事本中描述边界坐标如下：从左到右，依次代表：结点编号（从0开始），X坐标，Y坐标，三角形边长，边界标志。2022/11/2226第三章浅水模型的程序实现（0，0）（500，0）（500，第三章浅水模型的程序实现Step2：定义区域边界 把边的定义编号依次写入文件，如下：从左到右依次为：边的编号，边的起始结点编号，结束结点编号，边界标志。012301232022/11/2227第三章浅水模型的程序实现Step2：定义区域边界从左到右Step3保存文件为*.d，利用EasyMesh生成网格文件 *.n表示结点文件，*.e表示单元文件，*.s表示边文件最终的文件内容：+dxf表示同时生成*.dxf文件第三章浅水模型的程序实现2022/11/2228Step3保存文件为*.d，利用EasyMesh生成网格第三章浅水模型的程序实现最终网格：2022/11/2229第三章浅水模型的程序实现最终网格：2022/11/2229第三章浅水模型的程序实现3，模型数据结构设计从网格文件读入的数据包含大量的网格信息，必须定义适当的数据结构来描述，包括结点的数据结构，边的数据结构和单元数据结构； 本课题使用Delphi中的Record型数据结构（相当于C语言中的Structure结构体）来描述结点，边和单元；所有的结点放在一个动态数组里，在读取网格文件的时候确定数组的大小，利用数组下标进行结点数据的访问；边和单元的结构定义与结点相似。2022/11/2230第三章浅水模型的程序实现3，模型数据结构设计从网{TNode}//节点的数据结构TNode=RecordX:Real;//x坐标Y:Real; //y坐标q:Array[0..3]ofReal;NeibElement:TIntArray;//相邻单元End;第三章浅水模型的程序实现{TElement}//单元数据结构

TElement=RecordNode1:Integer;Node2:Integer;Node3:Integer;Element1:Integer;Element2:Integer;Element3:Integer;Side:Array[0..2]ofInteger;XCenter:Real;YCenter:Real;U:Real;V:Real;H:Real;Concentration:Real;Area:Real;AccumulateFlux:TFlux;End;2022/11/2231{TNode}//节点的数据结构第三章浅水模型的程序实现{TSide}//边的数据结构

TSide=RecordStartPoint:Integer;EndPoint:Integer;LeftElement:Integer;RightElement:Integer;Mark:Integer;SideLength:Real;Cos:Real;Sin:Real;qL:Array[0..3]OFReal;qR:Array[0..3]OFReal;Flux:TFlux;End;第三章浅水模型的程序实现2022/11/2232{TSide}第三章浅水模型的程序实现4，模型求解过程的封装－类的使用为了便于模块之间的协调，以及求解过程中数据的共享，我们把求解浅水方程的过程封装在一个TModelSolver类中，该类实现了从读入网格文件到输出计算结果的全部求解过程。TypeTModelSolver=ClassPrivateFL,FR,F:TFlux;FWaterLevel:Real;//开边界水位

FQPerWidth:Real;//开边界宽流量

ProcedureCalcNodeq;//按照面积加权得到结点的矢量

FunctionAverage_Outer(UL,VL,HL,CL,UR,VR,HR,CR:Real):TFlux;//左右单元数值通量平均

FunctionAverage_Inner(UL,VL,HL,CL,UR,VR,HR,CR:Real):TFlux;//左右单元先平均变量，再求数值通量

FunctionFVS(UL,VL,HL,CL,UR,VR,HR,CR:Real):TFlux;//通量向量分裂格式

FunctionFDS(UL,VL,HL,CL,UR,VR,HR,CR:Real):TFlux;//通量差分裂格式2022/11/2233第三章浅水模型的程序实现4，模型求解过程的封装－类的使用为PublicNode:TNodeArray;Element:TElementArray;Side:TSideArray;ProcedureCalcElementFlux();//计算每个单元通量之和

ProcedureIteration();//浅水方程的迭代计算

procedureReadFromFile(NodePath,ElementPath,SidePath:String);//从网格文件*.n,*.e,*.s读取数据到结构体数组ode,ElementSide中

ProcedureCalcSideLength();//使用两点间距离公式，计算每一边的长度

ProcedureCalcCosSin();//计算每一边的cos,sin值，以左单元(特指逆时针方向)为准，右单元的值符号相反

ProcedureCalcElementArea();//计算每一单元的面积

ProcedureExportToTecplot();//输出数据文件到Tecplot进行可视化处理

ProcedureCalcSideFlux();//计算通过每条边的通量

FunctionNonLinearSolver(UL,HL,QR:Real):Real;//使用Newton法求解非线性方程

ProcedureSetInitCondition_U(U:Real);//设置X方向流速U的初始值

ProcedureSetInitCondition_V(V:Real);//设置y方向流速V的初始值

ProcedureSetInitCondition_H();//设置水深H的初始值

ProcedureSetInitCondition_C();//设置浓度C的初始值

ProcedureSetBoundaryConditon(QPerWidth,WaterLevel:Real);//设置边界条件：开边界的水位和单宽流量ProcedureExportSideFlux();//输出边的数值通量到文件：调试用

ProcedureSearchNeibElement;//得到与结点相邻的单元End;第三章浅水模型的程序实现2022/11/2234Public第三章浅水模型的程序实现2022/11/2235，模型求解的流程第三章浅水模型的程序实现1，读入网格文件2，计算边长，外法线单元向量，单元面积3，寻求每一结点的相邻单元4，设置初始条件2022/11/22355，模型求解的流程第三章浅水模型的程序实现1，读入网格文件5，设置边界条件*.e,*.s,*.n6，计算边的数值通量7，计算单元的累计数值通量8，求得该时刻形心处的U，V，H和C9，回到6迭代计算第三章浅水模型的程序实现2022/11/22365，设置边界条件*.e,*.s,*.n6，计算边的数值通量7第四章数据可视化1，数据可视化的需求分析模型求解过程中，每个单元将会产生四个数据（X流速，Y流速，浓度，水位），而整个计算区域有成千上百个单元，并且单元数据是在不停的迭代中，所以在一定的时间间隔内，会产生大量的数据，要把握发展趋势，必须采用一定的数据可视化工具。本课题使用DirectX技术实现了流场和浓度场的可视化。流场的绘制等值云图的填充等值线的生成与TecPlot的接口2022/11/2237第四章数据可视化1，数据可视化的需求分析模型求解过程中，每第四章数据可视化2，流场的绘制流场箭头的绘制需要四个点的坐标：起点，终点，两翼两个点。起点利用三角形形心处的坐标。终点采用X方向的速度分量U和Y方向的速度V进行矢量相加，乘以一个速度坐标转化系数求得。同时利用几何关系求得箭头左右两翼的坐标，用直线对四个点分别相连。2022/11/2238第四章数据可视化2，流场的绘制流场箭头的绘制需要四个点的坐第四章数据可视化3，等值云图的绘制根据文献可知，目前等值云图的绘制有三种方法：一个单元内使用一种颜色进行填充，值越大，颜色越深扫描线算法：对单元内的每个象素点差值计算颜色值改进的扫描线法：在三角形内定义若干四边形，每个四边形内的颜色值是一样的，对四边形进行统一的颜色填充2022/11/2239第四章数据可视化3，等值云图的绘制根据文献可知，目前等值云第四章数据可视化本课题采用微软未公开发布的API函数GradientFill（），进行三角形的渐变填充。该方法相当于扫描线算法，但是，经过比较，要比扫描线算法快3倍以上。使用该函数时需要确定三个结点上的颜色值，三角形内的每个点的颜色值在这三个结点颜色值之间。2022/11/2240第四章数据可视化本课题采用微软未公开发布的API函数Gra第四章数据可视化4，等值线的绘制三角网格上的等值线绘制，首先要求得边上的差值点，然后追踪含有差值点的边，确定差值点的顺序，最后把求得的差值点进行光滑处理，在精度要求不高的情况下，可以直接把差值点相连，或者也可用Bezier曲线，B样条曲线来拟合。5，与Tecplot的接口Tecplot是著名的数据可视化软件，利用它可以实现动态显示过程*.avi文件的输出。在了解Tecplot的文件结构之后，我们完全有可能作出与它的数据接口。2022/11/2241第四章数据可视化4，等值线的绘制三角网格上的等值线绘制，首第五章验证算例1,溃坝算例-水流的验证2022/11/2242第五章验证算例1,溃坝算例-水流的验证2022/11/22第五章验证算例2022/11/2243第五章验证算例2022/11/2243第五章验证算例2,均匀流场中污染物的迁移-水质的验证2022/11/2244第五章验证算例2,均匀流场中污染物的迁移-水质的验证202第五章验证算例2022/11/2245第五章验证算例2022/11/2245第五章验证算例3,某河道的水流水质模拟见程序演示2022/11/2246第五章验证算例3,某河道的水流水质模拟结尾报告结束,多谢各位光临!2022/11/2247结尾报告结束,多谢各位光临!2022/11/2247基于有限体积法的二维水流水质模拟及其可视化研究疯狂猫2022/11/2248基于有限体积法的二维水流水质模拟及其可视化研究疯狂猫2022目录绪论浅水动力学的控制方程有限体积法原理浅水模型的程序实现数据可视化验证算例2022/11/2249目录绪论2022/11/222绪论1，本课题的主要工作 使用有限体积法（FiniteVolumeMethods）在无结构网格上离散二维浅水方程，求得水流速度场，并把求得的流场作为输入条件进一步求解对流扩散方程得出污染物浓度场。 同时，计算过程将会产生大量的数据，本课题使用Delphi语言＋DirectX技术实现流场和浓度场的可视化。2022/11/2250绪论1，本课题的主要工作2022/11/223绪论2，课题的现实意义 本课题的浅水模型可以用来实现河流、河口、湖泊、近海的水流水位模拟；而对流扩散方程是水质模型的基本控制方程，可以实现污染物（BOD，COD）、泥沙、热、盐的数值模拟 ，具有重大的现实意义。

应用举例：河道流量及洪水预报、溃坝决堤、洪水漫滩、河口潮汐、盐水入侵、都市排水、海上油膜扩散等。2022/11/2251绪论2，课题的现实意义2022/11/224第一章浅水动力学的控制方程1，何谓浅水 这里所说的浅水（ShallowWater）并非是指水深比较浅的水，而是具有它自身特定的含义。我们把满足以下四个方面的水称为浅水：

有自由表面； 以重力为主要驱动力，同时考虑风应力和地转柯氏力； 水平流速沿垂线近似均匀分布； 水平运动尺度远大于垂直运动尺度； 满足浅水假设的包括：河流、湖泊、河口和近海。2022/11/2252第一章浅水动力学的控制方程1，何谓浅水2022/11/22第一章浅水动力学的控制方程2，一维浅水控制方程

浅水方程又称圣维南方程（Saint－Venant），是由圣维南1871年在法国科学院汇刊73卷中提出的，它的一维形式为：

－－－－－－－－－－－－－连续方程－－－－运动方程式中：h－水位；v－流速；g－重力加速度；S0－水底坡度；Sf－摩阻坡度；x－水平距离；t－时间；2022/11/2253第一章浅水动力学的控制方程2，一维浅水控制方程－－－－－－第一章浅水动力学的控制方程下面我们来推导以一维浅水方程，考虑以下的一维单元模型：单元左侧进水流量为Q，右侧出水可以表示。右减左，得到该单元上的净出流量为，在Δt时刻内的净出水量为：。同时考虑Δt时刻内，单元控制体内水量增加为。根据单元水量平衡可以得到：令Δt－>0,Δx－>0,可以得到连续性方程。2022/11/2254第一章浅水动力学的控制方程下面我们来推导以一维浅水方程，考第一章浅水动力学的控制方程按照同样的思路，我们可以推导出动量方程：在Δt时刻内，左右界面净输出动量为； 左右界面水压差为； 水重在流向的分量为hgS0Δx;

底摩阻力在流向的分量为－hgSfΔx;控制体内的动量增量为：由控制体动量平衡得出动量方程。根据以上两个方程的推导，我们可以看到，其实浅水方程就是物理守恒定律（水量和动量）在有限体积上的应用。2022/11/2255第一章浅水动力学的控制方程按照同样的思路，我们可以推导出动3，从一维到两维

我们上述的浅水方程是在只考虑横向尺度（X方向）的情况下得到的，对上述方程进行纵向尺度（Y方向）扩展，可以很容易地得到二维的浅水方程：第一章浅水动力学的控制方程－－－－－－－－－－－连续方程－－－－X向运动方程－－－－Y向运动方程2022/11/22563，从一维到两维第一章浅水动力学的控制方程－－－－－－－－第一章浅水动力学的控制方程或写成向量形式：式中：2022/11/2257第一章浅水动力学的控制方程或写成向量形式：2022/11/4，求解浅水方程常用的数值方法

第一章浅水动力学的控制方程有限差分（FiniteDifferenceMethods） FDM以Taylor级数展开为工具，对微分方程中的导数项用差分项来逼近，通过求解差分方程得到近似解。优点：简单易学，处理效率高。缺点：使用的矩形网格对边界逼近较差。有限元（FiniteElementMethods）

把计算域分为有限个单元，分单元对解逼近，使微分方程空间积分的加权残差极小化。FEM在数学上适于求解椭圆型方程的边值问题，不适合求解以对流为主的输运问题。有限体积法（FiniteVolumeMethods）2022/11/22584，求解浅水方程常用的数值方法第一章浅水动力学的控制方程有第一章浅水动力学的控制方程 有限体积法和有限元一样将计算域分成若干单元。在计算出每个控制体边界沿法向输入（出）的流量和动量通量之后，对每个控制体分别进行水量和动量平衡计算，得到计算时段末各控制体的平均水深和流速。因此，FVM正是对于推导原始微分方程所用控制体途径的回归。 本课题使用FVM进行浅水方程的求解。 上述三种方法在具体应用中并没有绝对的好，和绝对的坏，根据具体问题的不同，可以用不同的数值方法。比如，DHI的二维浅水动力学软件Mike21采用的就是差分中的ADI－Quick格式，而著名的SMS（SurfaceWaterModelingSystem）使用的是有限元法，英国Wallingford的InfoWorks中的RS（RiverSimulation）模块使用的是最为传统的四点差分格式。2022/11/2259第一章浅水动力学的控制方程 有限体积法和有限元一样将计算域第二章有限体积法原理1，有限体积法离散原理有限体积法的基本原理是在被离散化了的计算区域上，计算出通过每个控制体边界沿法向输入或输出的流量和动量通量后，对每个控制体分别进行水量和动量平衡计算，最终得到计算时段末各控制体的平均水深和流速(假设水力要素在各控制体内均匀分布)。 考虑浅水方程的向量形式：在控制体上进行积分，得到：2022/11/2260第二章有限体积法原理1，有限体积法离散原理有限第二章有限体积法原理利用散度定理，对左边的第二项化为沿控制体边界的积分：S为控制体的边界，n为边界外法线方向的单位向量。对于有m个边的凸多边形而言， 等于被积函数在控制体各边上的法向值与该边长度的乘积，因此又可以写为： 2022/11/2261第二章有限体积法原理利用散度定理，第二章有限体积法原理其中，为每条边的数值通量令原始的矢量形式的方程化为至此，问题归结到如何求解Fn(q)2022/11/2262第二章有限体积法原理其中，第二章有限体积法原理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/11/2263第二章有限体积法原理2,从二维到一维第二章有限体积法原理这样，对于每个控制体我们可以得到最终的离散方程：上式左边表示控制体内守恒变量在Δt内的变化，右边第一项表示沿各边法向输出的通量之和，第二项表示控制体内源项(入流及外力)在Δ

t内的作用;这反映了守恒物理量的守恒原理:守恒物理量在控制体内随时间变化等于各边法向数值通量时间变化量之和及源项的时间变化量。2022/11/2264第二章有限体积法原理这样，对于每个控制体我们可以得到最终的第二章有限体积法原理到目前为止，有限体积法的计算归结为法向通量的计算。这是有限体积法的核心所在，法向通量的计算目前常用的途径是通过在沿外法向建立单元水力模型，并求解一维黎曼问题而得到。 已知条件是相邻控制体形心的守恒变量qi和qj（已通过旋转变换为外法向和切向）。i为进行FVM计算的控制体，j为第j边的相邻控制体。因浅水方程的旋转不变性，在局部坐标系中方程为根据控制体内变量的分布确定界面中点两侧的状态变量：qL对应于左侧，qR对应于右侧。通常在左右两侧存在间断。法向数值通量近似取为法向方向的一维黎曼问题的解。3，控制体界面处的水力模型与黎曼问题2022/11/2265第二章有限体积法原理到目前为止，有第二章有限体积法原理所谓黎曼问题（RiemannProblem），又称间断问题，是指如下的偏微分方程在初值为间断情况下的解间断初值条件：黎曼问题对于初学者较难理解，我们可以感性地理解为：t=0时刻，x＝0左侧水位为qL，右侧为qR,黎曼问题即是要求解t>0时刻在x＝0处的水位值。2022/11/2266第二章有限体积法原理所谓黎曼问题（RiemannProb第二章有限体积法原理4，FVS格式求解黎曼问题黎曼问题有许多种近似求法，例如通量分裂法(flux-vectorsplitingmethod)、通量差分裂法(flux-differencesplitingmethod)，及Osher数值方法等。本课题采用FVS格式求解黎曼问题。在此我们对FVS格式本身的推导不做过多的阐述，格式最终形式为：缓流：急流：法向数值通量FLR＝F＋＋F－2022/11/2267第二章有限体积法原理4，FVS格式求解黎曼问题第二章有限体积法原理5，边界条件当控制体的某一边处于计算区域的边界上时，那么该边界上的数值通量的处理就相应的转化成了边界问题。在此种情况下，被分析单元的水力要素状态向量qL是已知的，而相邻单元的水力要素状态向量qR，是未知的，因而必须确定qR，才可以利用FVS格式来估计法向数值通量。本课题涉及到三种边界条件：固壁边界，给定水位边界和给定单宽流量边界。如下图所示单宽流量边界水位边界固壁边界2022/11/2268第二章有限体积法原理5，边界条件第二章有限体积法原理(1)陆地边界(闭边界):

uR=uL;hR=hL;(2)给定水位时：

(3)给定单宽流量时，联立求解两个方程

qR=uR*hR

这两个方程最终可以化为一个非线形方程，本课题采用牛顿下山法迭代求解。2022/11/2269第二章有限体积法原理(1)陆地边界(闭边界):2022/1第三章浅水模型的程序实现1，软件构架设计一般的计算流体动力学（CFD，ComputationalFluidDynamics）软件，包括三个模块：预处理（Preprocessor）

主要进行计算域的网格划分，边界条件的设定；模型求解器（ModelSolver）

模型方程的离散求解；后处理（PostProcessor）

模型计算结果的数据可视化；2022/11/2270第三章浅水模型的程序实现1，软件构架设计第三章浅水模型的程序实现预处理模型求解器后处理网格计算结果本课题的预处理采用EasyMesh（Free，并且有源代码）软件进行三角网格划分；模型求解采用Delphi语言读取网格文件，求解FVM离散后的浅水方程，得到每个体积单元上的水位，流速（X方向和Y方向）和浓度；后处理采用Delphi支持下的DirectX技术实现计算网格的显示，流场的绘制，浓度场云图的绘制，浓度场等值线的绘制。模型求解和后处理整合在一个程序中，在求得结果之后，可以立即进行数据的可视化，从而达到动态显示的效果。2022/11/2271第三章浅水模型的程序实现预处理模型求解器后处理网格计算结果第三章浅水模型的程序实现2，预处理－网格的生成EasyMesh使用所谓的Delaunay法生成三角形网格，具体的生成过程在此不作深究，我们只需明白两点： 最终生成的三角形网格会覆盖整个区域，并且彼此之间没有重叠；任何一个三角形的外接圆中，不包含其他三角形的结点，这样的三角形均匀性最好；我们以一个正方形区域为例，阐述网格的生成过程：Step1：确定区域边界坐标2022/11/2272第三章浅水模型的程序实现2，预处理－网格的生成第三章浅水模型的程序实现（0，0）（500，0）（500，500）（0，500）该矩形区域有四个边界坐标：（0，0），（500，0），（500，500）和（0，500），我们在记事本中描述边界坐标如下：从左到右，依次代表：结点编号（从0开始），X坐标，Y坐标，三角形边长，边界标志。2022/11/2273第三章浅水模型的程序实现（0，0）（500，0）（500，第三章浅水模型的程序实现Step2：定义区域边界 把边的定义编号依次写入文件，如下：从左到右依次为：边的编号，边的起始结点编号，结束结点编号，边界标志。012301232022/11/2274第三章浅水模型的程序实现Step2：定义区域边界从左到右Step3保存文件为*.d，利用EasyMesh生成网格文件 *.n表示结点文件，*.e表示单元文件，*.s表示边文件最终的文件内容：+dxf表示同时生成*.dxf文件第三章浅水模型的程序实现2022/11/2275Step3保存文件为*.d，利用EasyMesh生成网格第三章浅水模型的程序实现最终网格：2022/11/2276第三章浅水模型的程序实现最终网格：2022/11/2229第三章浅水模型的程序实现3，模型数据结构设计从网格文件读入的数据包含大量的网格信息，必须定义适当的数据结构来描述，包括结点的数据结构，边的数据结构和单元数据结构； 本课题使用Delphi中的Record型数据结构（相当于C语言中的Structure结构体）来描述结点，边和单元；所有的结点放在一个动态数组里，在读取网格文件的时候确定数组的大小，利用数组下标进行结点数据的访问；边和单元的结构定义与结点相似。2022/11/2277第三章浅水模型的程序实现3，模型数据结构设计从网{TNode}//节点的数据结构TNode=RecordX:Real;//x坐标Y:Real; //y坐标q:Array[0..3]ofReal;NeibElement:TIntArray;//相邻单元End;第三章浅水模型的程序实现{TElement}//单元数据结构

TElement=RecordNode1:Integer;Node2:Integer;Node3:Integer;Element1:Integer;Element2:Integer;Element3:Integer;Side:Array[0..2]ofInteger;XCenter:Real;YCenter:Real;U:Real;V:Real;H:Real;Concentration:Real;Area:Real;AccumulateFlux:TFlux;End;2022/11/2278{TNode}//节点的数据结构第三章浅水模型的程序实现{TSide}//边的数据结构

TSide=RecordStartPoint:Integer;EndPoint:Integer;LeftElement:Integer;RightElement:Integer;Mark:Integer;SideLength:Real;Cos:Real;Sin:Real;qL:Array[0..3]OFReal;qR:Array[0..3]OFReal;Flux:TFlux;End;第三章浅水模型的程序实现2022/11/2279{TSide}第三章浅水模型的程序实现4，模型求解过程的封装－类的使用为了便于模块之间的协调，以及求解过程中数据的共享，我们把求解浅水方程的过程封装在一个TModelSolver类中，该类实现了从读入网格文件到输出计算结果的全部求解过程。TypeTModelSolver=ClassPrivateFL,FR,F:TFlux;FWaterLevel:Real;//开边界水位

FQPerWidth:Real;//开边界宽流量

ProcedureCalcNodeq;//按照面积加权得到结点的矢量

FunctionAverage_Outer(UL,VL,HL,CL,UR,VR,HR,CR:Real):TFlux;//左右单元数值通量平均

FunctionAverage_Inner(UL,VL,HL,CL,UR,VR,HR,CR:Real):TFlux;//左右单元先平均变量，再求数值通量

FunctionFVS(UL,VL,HL,CL,UR,VR,HR,CR:Real):TFlux;//通量向量分裂格式

FunctionFDS(UL,VL,HL,CL,UR,VR,HR,CR:Real):TFlux;//通量差分裂格式2022/11/2280第三章浅水模型的程序实现4，模型求解过程的封装－类的使用为PublicNode:TNodeArray;Element:TElementArray;Side:TSideArray;ProcedureCalcElementFlux();//计算每个单元通量之和

ProcedureIteration();//浅水方程的迭代计算

procedureReadFromFile(NodePath,ElementPath,SidePath:String);//从网格文件*.n,*.e,*.s读取数据到结构体数组ode,ElementSide中

ProcedureCalcSideLength();//使用两点间距离公式，计算每一边的长度

ProcedureCalcCosSin();//计算每一边的cos,sin值，以左单元(特指逆时针方向)为准，右单元的值符号相反

ProcedureCalcElementArea();//计算每一单元的面积

ProcedureExportToTecplot();//输出数据文件到Tecplot进行可视化处理

ProcedureCalcSideFlux();//计算通过每条边的通量

FunctionNonLinearSolver(UL,HL,QR:Real):Real;//使用Newton法求解非线性方程