目录
第8 章 有限体积法.2
§8.1 有限体积法的基本概念.2
8.1.1 几个术语2
8.1.2 控制体积的选择.3
8.1.3 结构与非结构网格.5
§8.2 有限体积法构造7
8.9-1 方程离散及基本格式7
8.9-2 物理特性要求.11
8.9-3 迎风型通量格式14
8.9-3 TVD 格式18
§8.3 非结构网格上的有限体积法23
8.3.1 基本方程.23
8.3.2 离散基本思路24
8.3.3 数值通量近似.25
第9 章 在水力学问题中的应用.29
§9.1 渗流问题中的应用.29
9.1.1 饱和—非饱和地下水运动基本控制方程29
9.1.2 方程的离散31
9.1.3 算例【陈扬 硕士论文】33
§9.2 二维明渠非恒定流计算.38
9.2.1 水流基本方程38
9.2.2 控制方程的离散39
9.2.3 计算实例.53
§9.3 三维紊动分层流计算.64
9.3.1 紊动分层流基本方程64
9.3.2 紊流模型及控制方程离散.65
9.3.3 压力校正法66
9.3.3 边界条件69
9.3.4 盐度引起的负浮力流动71

第8章 有限体积法
有限差分方法是从描述各种物理现象的基本微分方程出发构造离散方程的，前文已经对其作了翔实、周密的论述。该部分将从基础算法入手分析介绍在计算流体力学界广为应用的有限体积法。基于有限体积法的实用算法在计算流体力学、计算传热学等领域得到了飞速发展 [1-3]。在水力学诸多问题，如水流物质输运模拟，水工水力学模拟以及溃坝洪水波演进等水流模拟中也得到了广泛应用。
§8.1 有限体积法的基本概念
有限体积法，又称为有限容积法，它正是从物理量守恒这一基本要求出发提出的。这也是其受计算流体力学界广为称道和喜欢之处。其以守恒型的方程为出发点，通过对流体运动的有限子区域的积分离散来构造离散方程。有限体积法有两种导出方式，一是控制容积积分法，另一个是控制容积平衡法。不管采用哪种方式导出的离散化方程，都描述了有限各控制容积物理量的守恒性，所以有限体积法是守恒定律的一种最自然的表现形式。
该方法适用于任意类型的单元网格，便于应用来模拟具有复杂边界形状区域的流体运动；只要单元边上相邻单元估计的通量是一致的，就能保证方法的守恒性；有限体积法各项近似都含有明确的物理意义；同时，它可以吸收有限元分片近似的思想以及有限差分方法的思想来发展高精度算法。由于物理概念清晰，容易编程；有限体积法成为了工程界最流行的数值计算手段。
8.1.1 几个术语
在进行数值计算时，要把计算区域划分成一系列互不重叠的离散小区域，然后在该小区域上离散控制方程求解待求物理量。在有限差分法中只涉及到网格节点的概念，而有限体积法因为物理解释需要，形成了以下几个常用几何要素的相关名词。
控制体积（control volume）：图中阴影部分，方程积分离散时的小体积单元（二维为面积单元）。
单元（cell）：控制体积的中心，常用形心来表示，为待求物理量的几何位置。图中用空心园来表示，如点P、W、E、N、S 等。常用单元来代表整个控制体积。以下如若不作特殊说明，则用“单元”来代表控制体积。
网格线（grid lines）：用来分割计算区域内各控制体积的交错曲线簇，如图中折线N1N2、N2N3、N3N4、N4N1等
网格节点（nodes）：网格线之间的交点， 图中用黑圆点来表示，如N1、N2、N3、N4等单元界面（cell faces）：相邻两个控制体积间的公共面(二维则可以认为是公共边，这样看起来就和网格线一致了，但是要注意这不是同一个概念)，图中用小写字母e、n、w、s 表示。通常定义e、n、w、s 几何位置位于交界面的形心点，二维则认为在公共边的中心点。
8.1.2 控制体积的选择
当你开始用有限体积法模拟流体流动时，而且划分好网格后，你必须选定控制体积的形成方式。目前，常用的有两种方法：单元中心方式（cell-centered）和顶点中心方式
（vertex-centered）。另外一些学者还发展了由两种方式综合形成的混合方式。根据问题的特点和要求，不同的变量可以采用不同的控制体积，因此又产生了交错网格和同位网格的称谓，这里不再深入介绍，读者根据需要可以参考相关文献[1-3]。
前面我们我们在讨论有限体积法的术语时，已经看到了单元中心控制体积的形成方式。
这里再作一个说明并对两者作简单的比较。图8-2 为两种控制体积选择方式示意图。阴影部分表示单元的控制体积。空心圆点表示单元，实心圆点表示网格节点，在顶点中心方式中，单元和网格节点重合。

（1）单元中心方式
单元位于控制体积的中心，即将单一的网格单元作为控制体积，网格单元互不重叠。目前用非结构网格做流动计算多使用单元中心格式,它有利于处理陆地和给定流量边界条件。
（2）单元顶点方式
以网格节点为中心，它的形成有多种方式。其常用的构成方式是由以该节点为顶点的格子的形心以及各共顶点的网格线中心点的一系列连接线段所构成一个多边形控制体积。也可以由环绕该节点的一组格子组成控制体。
因此计算同样多的变量，单元中心方式变量布置简单直观，易于处理边界条件和保持离散的守恒性，而且需要的网格数要比单元顶点方式少得多，可节省计算时间。
8.1.3 结构与非结构网格
1 结构网格
结构网格，指的是一类具有一定的分布特征，可以用相应的行列关系来顺序描述的网格。
常用的矩形网格、曲线网格以及块结构网格等。
图 8-3 结构网格
矩形网格是FDM 或FVM 最为常用的网格系统。它不对原始的控制方程（直角坐标）作任何坐标变换，在边界处采用简单的阶梯形网格近似复杂的边界。或者，用规则矩形网格覆盖整个包括陆边界在内的计算区域，对陆边界及陆域处的网格作特殊的技术处理，如所谓的冻结法，即令某些控制体积“冻结”，或者说，“堵塞”某些控制体积，则剩下的“活动”的控制体积可形成所需要的不规则区域[1]，“冻结”亦即令该控制体积中的因变量在运算过程中始终保持不变。常采用大数值源项和大粘性系数法实现冻结区内流速恒为零。尽管这类对规则网格作特殊的技术处理的方法对边界概化过于严重，边界计算值过于粗糙，陆边界附近会形成虚假的曲折水流。但是，其网格生成方便，计算方法简单等优点使该方法得到了大量的应用。
在计算天然水域时，利用规则矩形网格阶梯近似边界不仅改变了原始物理边界附近的流态，而且增加了边界条件设置的复杂性；计算精度因此也大受影响，尤其在近边界高梯度区域。自然地，研究人员引入广义曲线坐标变换的思想。即在物理计算区域内构筑曲线网格，使得网格曲线与所求解区域的边界相重合，而后利用坐标变换将复杂的物理域变换到规则的计算空间内；在规则计算区域上离散求解变换后的控制方程，将得到的解投影至原物理区域得到数值解；或者在曲线网格上直接应用原始因变量求解。除将复杂的物理空间变换成为易于求解的规则网格之外，曲线坐标变换方法还能式计算网格点与物理域的活动网格节点相对应，从而可适应非恒定流动边界的模拟。目前常用的生成曲线网格的方法有微分方法和代数方法。微分方法常用求解椭圆型方程或双曲型方程实现。其中，尤以Thompson 为代表提出
的利用Poisson 方程生成贴体曲线网格的方法最为著名和流行。该类型方法可通过调整源项控制网格的分布，在梯度大的区域人工或自动加密网格。因此，该方法提出后，得到了广泛的发展和应用。

曲线网格应用遗留众多需要解决的问题。首先，大多数网格的生成方法只提供了离散点的变换，而不给出解析函数形式的变换关系。使用不光滑的网格时，对变换关系的差分近似会造成了很大的数值误差，甚至会导致不切实际的值。其次，如果网格严重偏离正交性，就会极大损坏原有的迭代方法的收敛速率。再次，因变量的选择也须谨慎考虑。在曲线网格中，可取原始笛卡尔坐标系变量或曲线坐标系中沿网格方向的协变量两种作为因变量。
满足正交性的网格具有优良的特性。在正交曲线坐标系下，变换后的某些项消失；控制方程更为简洁；数值稳定性和解的精度高；计算量和收敛速度也有了很大的改善。
但是，正交曲线网格极大的受制于边界的几何形状，对于天然水域，由于边界极其复杂，
很难生成完全正交的曲线网格，特别对于三维问题，正交曲线网格几乎是不可能实现的。因此，在实际工程水流计算上较少采用。
正是由于复杂边界下生成正交坐标系的困难，而非正交曲线网格不仅相对较容易生成，而且可随意调节网格的疏密，在高梯度区域人工或自动加密网格。因此，研究人员退而求其次，放弃完全正交性的要求，转而寻求非正交曲线网格生成和数值处理。如生成边界正交曲
线网格，或生成尽量接近正交的曲线网格，以及探求非正交的光滑或不光滑或网格的求解方法[2]等。为充分利用网格正交性的优点，常希望尽可能生成接近正交的网格，至少在接近
边界高梯度区域，以恢复数值近似损失的精度。为此，Sorenson[11]开发了一种能控制网格疏密和角点的边界正交曲线网格生成的微分方法。
2 非结构网格
八十年代以来，为了更好地拟合自然边界，非正交网格上的数值计算开发和应用得到了飞速的发展。
为了更好的说明非结构网格，可以先看图8-4。
图 8-4 非结构网格
从图中可以看出，非结构网格中单元格分布不再规则一致，其位置很难再凭借行列索引关系确定。非结构网格可以采用任意形状的单元格，单元边的数目也无限制，弥补了结构化网格不能够解决任意形状和任意连通区域的网格剖分的缺欠。实用上，为简化编程以及贴合边界要求，一般应用四面体、六面体（二维为三角形、四边形）等已经足够拟合天然水域边界。非结构网格最重要的一个特征是控制方程离散得到的代数方程的系数矩阵不再是结构网格下有规律的对角结构；若用对角形式存储，其带宽只能通过适当的布置单元编号顺序来减少。非结构网格原则上可应用于任何类型的数值方法，包括FDM，但FEM 和非结构网格上的FVM 算法更成熟，应用更广。

非结构网格最早用于FEM，但流体流动是高度非线性问题，而且FEM 计算量较大，这些问题使得基于FEM 的非结构网格技术未能在对流问题为主的地面水流（如浅水流动，水波运动等）计算上得到重视。八十年代以来，基于FVM 的非结构网格技术在空气动力学得到了广泛的发展和应用。九十年代开始一些专家学者根据浅水流动特征，将这些算法引入到计算浅水动力学中，并在模拟涌潮，溃坝等水力计算难题上取得了成功[5]。
与此同时，粘性流动的非结构网格FVM 模拟也开始出现。并在20 世纪90 年代中后期掀起了研究高潮[6]。作为全球计算流体力学软件供应商和技术服务商的Fluent 公司已经将最新的非结构网格研究成果集成，实现了研究成果的商业化。
非结构网格能很好地模拟自然边界及水下地形，利于边界调节的实现；便于控制网格密
度，易作修改和适应性调整；网格生成有众多富有成效的方法和自适应技术[9]，比曲线网格更易得到高质量的单元格。但是非结构网格应用也带来一系列需要解决的问题：单元格排列不规则，须建立相应的数据结构存储单元格信息；控制方程离散得到的代数方程的系数矩阵是高度稀疏的非对角型矩阵，需要寻求合适的存储方式及解法；隐格式较难实现，粘性项处理困难；数值解后处理工作量大；二阶非结构FVM 较易实现，若要扩展到高阶格式，则
需花费较大的代价。
§8.2 有限体积法构造
8.9-1 方程离散及基本格式
1 基本思路
与有限差分法不同，有限体积法是基于守恒型的控制方程。这里以一维对流方程为例，
说明有限体积法离散的基本思路。

为推导离散化的方程，我们必须事先对计算区域进行网格划分，然后在各个单元上进行积分平均。假设我们对一维计算域划分得到一个网格系统，如图 所示的为该网格系统中某单元邻接关系示意图。重点考察单元P，单元E 和W 分别为它的两个相邻单元。虚线为单元界面，图中用小写字母表示，界面e、w 包围形成的阴影部分即为单元P 的控制体积。控
制体积实际是三维体的的概念，只不过在考虑一维控制体积时，实际上把y、z 两个坐标方向假设成单位厚度，这样在一维问题中看起来成为了线的概念。

2 状态变量分布近似
从前面可以知道，用有限体积法推导离散化方程时，必须确定物理量的局部分布，这是离散过程极为重要的一步，不仅关系到守恒性是否保持，而且关系到算法精度，这是由于不同的状态变量分布对应了不同的离散格式。

值得指出的是，在用有限体积法计算时，和有限差分法一样，方程的解是用单元节点上离散点值构成的，而不关心单元间的状态变量是怎么变化的，也就是不关心解的分布。这和有限单元法中一旦选定了分布曲线，就确定了状态变量的分布函数是不同的。尽管我们在离散方程的时候要假定单元分布曲线，但是这只是为了推导公式时计算积分近似而采用的一些辅助关系式。一旦离散化方程推导出来了，就可以不用再管这些分布曲线近似。这种观点使我们在采用何种分布曲线近似方法有完全的自由。在积分离散时，根据数值模拟的需要，对控制方程中的每一项都可以采用不同的分布曲线来近似单元界面上的状态变量，而不必要追求近似假设的一致性。例如，我们用一阶迎风格式离散对流扩散方程时，对流项中因变量是用阶梯型分布近似，而扩散项中也取为阶梯分布的话，则根本导不出离散式，至少要有二阶的分段线性分布才行。上面已经提到，分布曲线近似关系到算法的精度，实际上，有限体积法中，每一种不同的状态变量分布近似，体现了不同的离散格式的几何构造方法。如对流问题中，一阶迎风格式为阶梯型分布，中心格式采用的分段线性分布，更为高阶的格式则需要更多的单元加入，形成复杂的状态变量分布近似关系，例如对流项二次迎风插值QUICK，分段抛物插值PPM 等。
3 基本格式
利用这些交界面表达式可以得到具体的离散化的代数方程，这里不在叙述。在这里，我们还可以看出，有限体积法离散的思路和有限差分法的有着本质的区别，有限差分法是应用Taylor 展开式近似微分控制方程中的导数，而有限体积法则是通过积分将导数项化成交界面的状态变量的表达式，然后，代入根据局部分布近似导出的交界面值得到完整的离散化的代数方程。
对于中心格式和一阶迎风格式的特点和有限差分法的一致。中心格式精度较高，具有二阶精度，但是不能反映对流传输机制。对非线性问题，常会出现非线性不稳定，为了消除该影响，需要添加人工粘性。一阶迎风格式稳定性好，但是数值耗散大。值得强调的是，一阶迎风格式离散思想中蕴涵的物理本质，即波动传波的信息主要来自上游方向，是我们构造高分辨率格式的基础。

8.9-2 物理特性要求
1 守恒性
如果对一个离散方程在定义域的任一有限空间内作求和的运算（相当于连续问题中对微分方程作积分），所得的表达式满足该区域上物理量守恒的关系时，则称该离散格式具有守恒特征。
我们知道，根据流体运动的不同特征，可以采用拉格朗日（Lagrange）法和欧拉（Euler）法来描述流体的运动。其中欧拉法着眼于研究固定空间位置上各流体质点运动要素的空间变化情况，在推导流体运动的控制方程时，欧拉法被广为采用。让我们回顾一下连续方程的推导过程，首先在流场中划定一个固定的空间区域作为流体运动的控制体积，考察流体流入、流出该控制体积的情况，利用质量守恒原理，可导出流体运动的连续性方程。其他描述流体中物质（能量）输运的控制方程也可根据各自的守恒律推导出来。在数值计算中，物理量的守恒性也必须得到保证，否则会导致非物理解的产生。因此守恒性成了数值求解过程中非常重要的一个方面。
有限体积法正是从物理量守恒这一基本要求出发提出的。它有两种导出方式[3]，一是控制容积积分法，另一个是控制容积平衡法。不管采用哪种方式导出的离散化方程，都描述了有限各控制容积物理量的守恒性。就像微分方程是表示关于无穷小控制体积内物理量的守恒一样，帕坦卡[1]（S.V. Patankar）非常生动的对此作了注解：如果把我们想象成处于微积分以前的时期，控制容积方程就会成为我们表达守恒原理的唯一途径。有限体积法的离散化方程满足了单个控制体积的平衡，当然在整个计算区域内，诸如质量、动量等物理量的积分守恒也就都能精确得到满足。无论在数值计算中采用巨大数目的细网格和少数的粗网格，数值解也照样显示准确的积分平衡。有限体积法的离散思想自动满足守恒定律，如质量守恒，动量守恒，能量守恒等等。所以有限体积法是守恒定律的一种最自然的表现形式。
为了得到守恒型的离散方程，需要满足一定的条件。这里仍以前文导出的全离散的守恒型显格式离散方程（8-12）为例。一般对l + k +1个单元（称为节点模板），守恒型数值通量

守恒格式可以自动算出激波并给出其正确的位置，因此这种数值方法称为激波捕获法（shock-capturing method）。如果使用非守恒格式，则需要在有关位置使用局部Rankine—Hugoniot 条件，以获得准确的激波位置，这种方法称为激波拟合法(shock-fitting method)。
2 迎风性
在有限差分法部分中的基本理论部分，我们已经深入的讨论了双曲型方程（对流）的特征线概念。特征值体现了微分方程的解（或者说扰动、信息）的传播方向。它既有很强的数学意义，更反映了问题的物理机制，表达了波动、能量等的传播方向。这启示我们在理论分
析和数值方法的设计上，应该充分考虑这一物理特性，与之相适应，而不可反其道而行之。根据特征（信息）传播方向，设计得到的数值方法可统称为特征型或迎风型格式。例如前文
提到的传统的一阶迎风格式，尽管非常简单，却是蕴含了最原始的的迎风设计思想，即在计算单元界面处的数值通量时，应根据特征传播的方向性，对于向右传播的分量，应该使用左
边的单元值来计算(因为影响来自左边)；对于向左传播的分量，应该使用右边的值来计算(因为影响来自右边)。
当前，这种根据信息传播方向设计数值格式的迎风思想已经成为了流体数值计算发展的重要基石，它不断地启发研究人员构造出各种各样的离散格式，如特征迎风格式，矢通量分裂格式，通量差分裂格式等。并在这些基于迎风思想的格式基础上，发展了分门别类的高阶、高分辨率格式，如MUSCL 格式、QUICK 格式、ENO（WENO）格式等等。
3 TVD 性

既要利用格式的高阶精度，又能使格式具有TVD 性质，从而求得高精度的物理上有意义的解，是研究人员的目标，为此，过去的二十几年中，各种二阶（甚至更高阶）的TVD 格式得到了大量的发展。这类高阶TVD 格式开发应用成为了计算流体力学界最前沿的研究。
8.9-3 迎风型通量格式
对于双曲型问题，不同特征分量传播的方向不同。考虑传播方向的影响从物理上分析更为合理。我们对传统的一阶迎风格式已经作了介绍。这里我们将介绍20 世纪80 年代发展起来的通量向量分裂格式（Flux Vector Splitting），通量差分裂格式（Flux Difference Splitting）以及Godunov 法和近似黎曼解算子法（Approximate Riemann Solvers）。
1 矩阵分裂
实际流动的控制方程是一个非线性系统，各个变量之间相互耦合。我们这里考虑一维的双曲型方程组
2 通量向量分裂格式（FVS）
Steger 和Warming 根据矩阵分裂方法，于1981 年首次提出了一种通量向量分裂格式。
最早该格式应用于空气动力学计算上，20 世纪90 年代，开始扩展到计算水力学领域。

在把通量向量格式应用到水力学中时，仍需要注意方程组的特性。因为水流运动的控制方程——圣维南（Saint Venant）方程组的通量项不具有齐次性，需要补充一个冗余的能量方程得到具有齐次通量的方程组，然后就可以按照前面的方法建立离散化方程。
3 通量差分裂格式（FDS）
通量向量分裂格式是将通量函数按不同特征方向进行分裂，然后积分离散。与通量向量分裂格式不同，通量差分裂格式是对任意区间上通量函数的差进行分裂的。
4 Godunov 格式和近似黎曼解法

8.9-3 TVD 格式
在前文讨论中，我们已经知道，满足TVD 特性的格式求出的数值解才是真正符合实际的物理解。然而，一阶迎风格式尽管具有TVD 性，但是数值耗散过大，将本来应该很陡峭的激波给抹得很平．为了在解光滑的区域获得较高精度的解，并且在激波等间断区域获得没有数值振荡或者是数值振荡可以令人接受的较陡峭的数值解，人们开始研究高分辨格式（High Resolution schemes）。1983 年Harten 在他的论文中提出了高分辨率和总变差减少的概念，首先证明了计算格式具有TVD 性质的充分条件，具体构造了修正通量的高分辨率（二阶）TVD 格式。他的成果在数值算法的研究上具有十分重要的意义，，为构造和分析高分辨率格式提供了有力的工具。随后，各种各样的TVD 格式如雨后春笋般的涌现出来，因其具有对间断的自动高分辩识别能力，而得到了广泛的应用。这种格式具有两个特点：（1）对标量非线性方程及常系数双曲型方程组，格式的解是总变差递减的（Total Variation Diminishing）；（2）与守恒律和熵不等式是相容的。TVD 格式提高了对激波的捕捉能力，但其不足之处在极值点上格式只有一阶逼近精度。

§8.3 非结构网格上的有限体积法
前面主要对有限体积法基本概念和离散格式作了介绍。在这节中，我们将介绍二维非结构网格上的有限体积法，以便于应用它来模拟自然中复杂区域内的流动及物质输运现象。本节只对算法的空间离散进行讨论，因为时间的离散和有限差分法一致，因此，不在专门介绍。

标签： 点击： 评论: