4.3 常用离散方法简介
随着计算机技术的应用，平面二维数学模型已广泛应用于水流流场及其输运物质的数值模拟，也产生了基于不同数值计算方法的数学模型。
对于一个理想的数学模型，大致希望其具有以下几个特点：
（1）计算稳定，这是模型解合理的基本要求。
（2）解具有收敛特性，数值解收敛于真解。
（3）守恒性，在计算域内能够保持水量守恒和动量守恒。
（4）计算速度快，虽然随着计算机的发展，计算速度不再是影响模型推广应用的主要问题，但在许多情况下计算速度仍然不满足实际需要，尤其是在洪水实时预报中更为明显。
（5）高精度，这是随着计算模型的发展，人们提出的更高要求，希望数值解更接近于真解。
（6）模型的健全性和通用性，所建模型能够适用于不同地形、不同流态、不同网格布置等多种情况。
4.3.1 常用离散方法
在平面二维数值模拟中，目前常采用的方程离散方法为有限差分法、有限元法、有限体积法等。每种方法和格式均有其自身的优点和适应性，在实际计算中选择什么数值方法应根据研究问题的特点和计算精度的要求以及研究者的习惯而定。根据前人的总结及作者的认识，分述如下。
1.有限差分法
有限差分法（Finite different method，FDM）是计算机数值模拟最早采用的方法，至今仍被广泛运用（王船海，1991；周雪漪，1995；郑邦民，1998；谭维炎，1998等）。该方法将求解域划分为差分网格，用有限个网格节点代替连续的求解域。有限差分法以泰勒级数展开等方法，将控制方程中的导数用网格节点上的函数值的差商代替进行离散，从而建立以网格节点上的值为未知数的代数方程组。该方法是一种直接将微分问题变为代数问题的近似数值解法，数学概念直观，表达简单，是发展较早且比较成熟的数值方法。
对于有限差分格式，从格式的精度来划分，有一阶格式、二阶格式和高阶格式。从差分的空间形式来考虑，可分为中心格式和逆风格式。考虑时间因子的影响，差分格式还可以分为显格式、隐格式、显隐交替格式等。目前常见的差分格式，主要是上述几种形式的组合，不同的组合构成不同的差分格式。差分方法主要适用于有结构网格，网格的步长一般根据实际地形情况和柯朗稳定条件来决定。
构造差分的方法有多种形式，目前主要采用的是泰勒级数展开方法。其基本的差分表达式主要有三种形式：一阶向前差分、一阶向后差分、一阶中心差分和二阶中心差分等，其中前两种格式为一阶计算精度，后两种格式为二阶计算精度。通过对时间和空间这几种不同差分格式的组合，可以组合成不同的差分计算格式。
在此以函数f（x，t）的一阶偏微分
为例，介绍差分的几种基本形式（陆金甫，1988；杨国录，1992）：
一阶向前差分
一阶向后差分
一阶中心差分
二阶中心差分
将向前差分格式和向后差分格式加权平均，可得到迎风差分格式：
式中，β为加权迎风因子，改变β的大小可以改变差分格式的迎风效应。当β=1时，为向前差分；当β=-1时，为向后差分；当β=0时，为中心差分。如果取
，则差分格式具有自动迎风效应，式中
为实际水流速度，
为临界速度。如果某一点水流速度小于临界速度，水流处于缓流，该点的水沙条件改变，可以向坐标两侧传播，水流从哪一侧流来，哪一侧的条件就对水流的影响较大，流速越大，影响越大，因此式（4.36）具有自动迎风效应。如果某一点水流速度不小于临界速度，水流处于临界流或急流，该点的水沙条件改变，只可以向坐标一侧下游传播，此时取β=±1，即为向前（后）差分，也称简单迎风差分；若某一点水流速度
，β=0，即为中心差分。
同理，对时间项t的基本差分形式也可以根据上述形式写出。
在二维函数中，对y向的差分形式也可以按照x向的形式写出。此外，还需要确定二阶混合偏导数的差分形式，通过fi+1，j+1、fi-1，j-1、fi+1，j-1 和fi-1，j+1四点的泰勒展开式便可得到：
有限差分法的研究颇为完善，差分格式繁多，目前在河道和浅水数学模型中常用的几种差分格式有Lax格式、蛙跳格式（Leap-frog），以及前面所用的Preissman四点偏心隐格式等。其中Lax格式为双层显格式，涉及到两个时间层和三个计算节点。Lax格式的基本构造思路是对变量f的一阶时间偏微分采用空间加权向前差商逼近，一阶空间微分采用中心差分逼近的方法；蛙跳格式是一个三层二阶显格式，对x和t方向均取中心差分格式。有限差分方法的格式很多种，但每种计算格式都从提高计算稳定性或计算精度等方面，对时间和空间几种基本差分格式进行不同的组合。
从离散的代数方程组求解的角度来讲，有限差分方法的计算格式还可以分为显格式、隐格式和显隐交替格式等。就二维方程来说，目前常采用的方法有交替方向法（ADI）、算子分裂法（其中包括方向分裂法、物理现象分裂法和时间分裂法等）、预测校正格式等。这些数值格式均是有限差分法和不同代数方程组求解方法的组合。
2.有限元法
有限元法（Finite element method，FEM）的基础是变分原理和加权余量法（姜礼尚，1979；何穷，1981；蒋孝煜，1984；李开泰，1988等）。其基本求解思想是把计算域划分为有限个互不重叠的单元，在每个单元内，选择一些合适的节点作为求解函数的插值点，将微分方程中的变量改写成由各变量或其导数的节点值与所选用的插值函数组成的线性表达式，借助于变分原理或加权余量法，将微分方程离散求解。采用不同的权函数和插值函数形式，便构成不同的有限元方法。
有限元法最早应用于结构力学，后来随着计算机的发展慢慢地应用于流体力学数值模拟。在有限元方法中，在每个单元内选择基函数，用单元基函数的线形组合来逼近单元中的真解，整个计算域上总体的基函数可以看做是由每个单元基函数组成，则整个计算域内的解可以看做是由所有单元上的近似解构成。在河道数值模拟中，常见的有限元计算方法是由变分法和加权余量法发展而来的里兹法和伽辽金法（Galerkin）、最小二乘法等。根据所采用的权函数和插值函数的不同，有限元方法也分为多种计算格式。从权函数的选择来说，有配置法、矩量法、最小二乘法和伽辽金法；按计算单元网格的形状划分，有三角形网格、四边形网格和多边形网格；按插值函数的精度划分，又分为线性插值函数和高次插值函数等。同样不同的组合构成不同的有限元计算格式。
对于权函数，伽辽金法是将权函数取为逼近函数中的基函数Φi；最小二乘法是令权函数等于余量本身，而内积的极小值则为对代求系数的平方误差最小；在配置法中，先在计算域Ω内选取N个配置点xi（i=1，2，…，N），令近似解在选定的N个配置点上严格满足微分方程，即在配置点上令方程余量为0。
插值函数一般由不同次幂的多项式组成，也有采用三角函数或指数函数组成的乘积表示的，但最常用的多项式插值函数。有限元插值函数分为两大类：一类只要求插值多项式本身在插值点取已知值，称为拉格朗日（Lagrange）多项式插值；另一类不仅要求插值多项式本身，还要求它的导数值在插值点取已知值，称为哈密特（Hermite）多项式插值。单元坐标有笛卡儿直角坐标和无因次自然坐标，有对称和不对称等。常采用的无因次坐标是一种局部坐标系，它的定义取决于单元的几何形状，一维看做长度比，二维看做面积比，三维看做体积比。在二维有限元中，三角形单元应用得最早，近来四边形等参元的应用也越来越广。对于二维三角形和四边形电源单元，常采用的插值函数为有拉格朗日插值直角坐标系中的线性插值函数及二阶或更高阶插值函数、面积坐标系中的线性插值函数、二阶或更高阶插值函数等。
有限元离散方法能灵活地处理复杂边界问题，在二维数值模拟计算中也有较多的应用（董文军，1996；Moural Heniche，2000等）。但主要因其在急变流动计算中易出现速度的坦化及计算速度的限制，在洪水演进等非恒定性较强问题中没有得到广泛应用。
3.有限体积法
有限体积法（Finite volume method，FVM）是近几十年来发展起来的一种离散方法，目前在浅水流动模拟中已开始推广应用，并有着较好的发展前景（D.H.Zhao，1994；D.Ambrosi，1995；谭维炎，1991—1998；胡四一，1991，1995；S.Chippada，1998；C.Zoppou，2000；S.Roberts，2000；Ji-Wen Wang，2000等）。该方法从水量和动量守恒的物理概念出发，将待解的微分方程对控制体进行积分，便得出一组离散方程。有限体积法可以采用矩形网格，也可以采用不规则网格进行离散计算域。网格可为三角形、四边形或其他多边形。有限体积法的优点是物理意义明确，对边界的适应能力较强。同时由于该方法大多采用守恒型离散格式，在整个计算域内积分守恒，易处理非线性较强的流体运动问题和函数间断的情况。这也是有限体积法比较有发展潜力的原因。近年来有限体积法发展较快，尤其是在空气动力学中应用较为广泛。
在有限体积法中，为了保证解的合理性和满足总体平衡，在选择插值公式和具体的离散分析中需要满足以下四条基本法则（帕坦卡，1984）：①控制体交界面的连续性；②正系数；③源项的负波线性化；④相邻节点系数之和。
由于对控制体交界面处法向通量的求解思路不同，有限体积法形成不同的计算格式。从计算离散的格式进行分类，可分为中心格式和逆风格式两大类（朱自强，1998）。中心格式最大的优点是计算简单，精度高，但其稳定性差，容易产生虚假震荡，在计算间断时需要添加人工黏性。Lax-Wendroff格式可以看做是该格式的代表。逆风格式也称迎风格式，传统的逆风格式是根据流速的方向来选择差分的方向，其优点是计算稳定性好，但对间断处理时计算精度不高。后来采用特征线来构造逆风格式，是逆风格式的一个很大的突破，采用逆风格式计算间断时既可以防止假振，又有很高的分辨率，采用守恒格式时对连续流和间断流都适用，但计算量比较大。现在常用的通量分裂法FVS格式、通量差分裂法FDS格式、Osher格式等均可以看做特征逆风格式。下面简要介绍一下有限体积方法的几种常见的计算格式（谭维炎，1998等）。
FVS（Flux Vector Splitting通量向量分裂）格式是由Steger和Warming于1979年提出的，目前在空气动力学中已得到广泛应用。该格式的构造思路是按照特征符号对法向通量F或其雅可比矩阵进行分裂，便于在选择插分时考虑信息传播的方向。
Godunov最早提出从求解一维黎曼问题的思路来建立数值格式。该方法的基本思路是在每一控制体内用解的均值来逼近真实解，则在控制体交接面处存在间断，求解黎曼问题，即可得控制体的法向通量。后来发展的许多格式都是以该方法的思想为基础的，即以求解黎曼问题为出发点构造数值格式，目前主要分为两类：Roe最早提出的FDS（Flux Difference Splitting，通量差分裂方法）和Osher格式。FDS格式和FVS格式相似，都是根据特征值进行分裂，其间的区别在于FVS格式是对通量向量进行分裂，而FDS格式是对任意区间上的通量差进行分裂。Osher格式假定控制体边界左右单元的变量之间通过m个稀疏波及压缩波连接，确定积分路径后，对准确的非线性黎曼问题的近似求解，即可计算出界面通量。
1983年，Harten首先证明了计算格式具有TVD性质（Total Variation Diminishing，全变差缩小）的充分条件并具体构造了修正通量的二阶TVD格式。TVD格式是单调保持格式，一般先将适用与双曲守恒律的守恒型差分格式写成黏性形式，通过对黏性系数的选择和控制而保持格式的全变差减小性质。目前TVD格式已经发展成一个格式族，有一阶格式、二阶格式、显格式、隐格式、逆风格式和对称格式等。
另外目前常见的还有FCT（Flux-Corrected Transport，通量校正运输）格式、ENO（Essentially Non-Oscillatory，基本无震荡）格式和MUSCL（Monotone Upstream-centered Scheme for Conservation Laws）方法。FCT格式是一种对格式黏性进行控制而建立起来的单调保持格式，为了防止震荡，在中心格式加一项系数为负的二阶导数项——反扩散项对通量进行校正。ENO格式是基于TVD格式的构造，放松了总变差的条件，从而提高了整个计算域上的计算精度。MUSCL方法是Van Leer于1979年提出的。该方法的特点是通过对初始状态时控制体界面处物理量的处理，采用一阶通量公式可达到二阶计算精度。这类算法的基本思路是：在控制体内的物理量采用线性分布，各控制体界面两侧的状态（也称状态向量）可由控制体平均值（即其中心处的值）和空间变率（坡度）来表示；由于界面两侧的状态有时存在间断，则选择一种一阶守恒、逆风、单调的算法求解相应的黎曼问题，选择不同格式可构成不同的MUSCL算法。
有限体积法对于求解地形比较复杂、水流间断的问题有其明显的优势，因此该方法在洪水演进计算中应用较多。但采用无结构网格时，计算中需要花费较多的时间建立网格之间的关系，影响了其计算速度。
4.有限分析法
有限分析法（Finite analysis method，FAM）是美籍华人陈景仁（J.R.Chen）在1987年提出的，该法是在局部单元上线性化微分方程和插值近似边界的条件下，求局部单元上的精确解，而构成整体的线形代数方程组。有限分析法有较高的计算精度，并具有自动迎风特性，计算稳定性好，收敛较快。但由于有限分析系数中含有无穷级数，单元系数计算较复杂，这给实际计算与理论分析都带来了一些困难。
在20世纪80年代末90年代初，李炜等提出混合有限分析法（HFAM），这种方法是在有限分析方法思想的基础上，把建立精确差分格式与非定常项的差分处理结合起来，建立一种混合有限分析格式。这种格式保持了有限分析法的优点，同时避免了无穷级数带来的不便，得到了广泛的应用。
5.其他方法
除了上述几种主要的数值离散格式之外，另外还有边界元法（姚寿广，1995等）、插值元法（汪定洋，1986）等，在此不作详细介绍。
4.3.2 不同离散方法的比较
在数学模型的建立过程中，虽然采用的数值离散方法不同，但都有相同的特点，即首先把计算区域划分成许多网格单元，然后在这些单元上把微分方程离散成代数方程，再把各单元上的代数方程汇合成总体代数方程组，最后在一定的初边值条件下求解此方程组，从而求得计算区域内各节点的物理量。有限差分法、有限元法以及有限体积法是微分方程离散的数值方法，这几种方法并不是完全不同，它们之间存在着一些联系。
有限元法和有限差分法在离散方程的思路上出发点不同，其主要区别在于（金忠青，1989）：①有限差分法是点近似，采用离散的网格节点上的值来近似表达连续函数，而有限元法是分段（或者分片、分块）近似，在单元内近似解是连续解析的，在单元之间近似解是连续的；②有限元法得到的是一个充分光滑的近似解，在单元内导数存在，而有限差分法一般不能保证解的光滑性；③有限元法对计算单元的划分没有特别的限制，处理灵活，特别是在处理复杂边界的问题时，这一优点更为突出。
就离散方法而言，有限体积法可视做有限单元法和有限差分法的中间物。有限单元法必须假定函数值在网格点之间的变化规律（即插值函数），并将其作为近似解。有限差分法只考虑网格点上函数值φ而不考虑φ在网格点之间的变化情况。有限体积法只寻求φ的节点值，这与有限差分法相类似；但有限体积法在寻求控制体积的积分时，必须假定φ值在网格点之间的分布，这又与有限单元法相类似。在有限体积法中，插值函数只用于计算控制体积的积分，得出离散方程之后，便可不用插值函数；如果需要可以对微分方程中不同的项采取不同的插值函数。
不同数值格式将微分方程离散为不同的代数方程组，即使采用同一种数值方法，但选用不同的格式，对同一微分方程也可以得到多种不同的代数方程组。但选择不同的数值格式时，若选择某种特定的格式，也可能将同一个微分方程离散为相同的代数方程组。因此，在某种意义上来讲，各种数值格式是相通的。例如，在有限体积法中，计算控制体交界面处的法向通量时，如果采取相邻控制体形心处通量的平均值，若采用的计算网格为矩形网格，则相当于二阶的中心有限差分格式；若采用三角形或四边形网格，又相当于线性三角形或双线性四边形单元的伽辽金有限元格式（谭维炎，1998）。
从整体上来考虑，在处理问题上不同的算法各有其自身的优缺点。有限差分法的原理比较简单，若边界不复杂，数学推理和程序编制都比较简单，并且计算内存占用比较少，计算速度快，在非恒定性比较强的问题中应用比较多。有限差分法最初一般采用矩形网格时，对复杂边界的适应性较差，计算精度不高，在复杂边界上常常出现虚假流动的现象，有时会出现水量不守恒现象，随着网格生成技术的发展，出现贴体正交曲线网格之后，解决了边界适应性差的缺陷。目前，也有学者建立了基于三角形网格的有限差分模型。因为有限差分法是用泰勒级数展开的，且取有限项，因此存在截断误差，即使泰勒级数取二阶项，截断误差也是二阶的，但边界上只取一点，因此边界上的精度只是一阶的。
有限元法最早用于固体力学，后来发展应用到流体力学上。它比有限差分法精度高，由于网格可以划分成任意形状，因此边界条件也能精确满足。相对有限差分法而言，其数学推理和程序编制就复杂些，但该方法的网格可以为三角形或四边形等任意网格，这种网格结构易于处理边界和地形比较复杂的问题，但需要数据结构定义计算单元之间的位置关系，占用内存相对较大，计算速度相对较慢。
有限体积法物理概念清楚，在求解控制体交接面法向通量时易于处理间断问题，可用来处理水流间断等问题，能保证水量和动量的守恒，即使对于较粗的网格也能得到较精确的解。该方法既可采用有结构的贴体正交曲线网格，又可采用三角形或任意四边形这类无结构网格，因此在河流数值模拟中得到了广泛的应用。
有限分析法的思想是通过在“单个元素上求得几个解析解再叠加来满足微分方程”，因此精度较高。对于边界，有限分析法可用到边界上3个点，这样可以提高边界上的精度。有限分析法有许多方案，其中包括有限体积法，它是在控制体上进行积分的，每个控制体包含一个节点，然后，用各种近似方法，把积分式化为代数式。该算法可采用贴体正交曲线网格，可较好适应不规则的河道边界条件，计算精度高，由于需要求解e的指数函数，速度相对较慢。
从加权余量方法的思路出发（金忠青，1989），剖析了目前常用的数值格式之间的关系，见图4.11。
图4.11 常用数值格式关系
4.3.3 常用求解方法
总之，以上介绍的数值计算方法各有其优缺点，但其共同的特征是首先需对计算区域进行网格划分，然后在分片单元上把微分方程离散成代数方程，再把分片单元上的代数方程汇合成总体代数方程，最后在一定的初边值条件下求解此方程组，从而求得计算区域内各节点的物理量。
对于代数方程组的求解，从迭代的方式来分，可分为显格式、隐格式和显隐交替格式。从求解技巧方面来说，现在常采用的有分步法、预测校正格式等。
显式计算方法的优点是不必求解方程组或迭代即可计算出下一时间段的函数值，程序简单，但因考虑到稳定性，其时间步长受到限制。隐式算法在线性分析中是无条件稳定的，可以加大时间步长，节约整个计算时间，但需要求解方程组或迭代求解，计算工作量大，占用内存多，程序编制比较复杂，每一时间步花费时间多。目前隐式计算一般用在有结构规则网格中，无结构网格中应用不多。
1.分步法
分步法就是将原来比较复杂的问题（高维或复杂的物理过程）的求解过程分解成若干个比较简单的问题的求解方法。目前常见的一般有以下三种概念上的分步方法：①空间概念上的分步：将高维空间问题分解成若干个低维问题分别求解，例如将二维空间问题分解成两个一维问题，这样可以直接利用原有的一维数值模拟方法，同时可以节省计算量和存储量；②物理概念上的分步：将同时包含若干个不同物理性质的复杂的物理过程分解成若干个单一物理过程的时间序列分步求解，例如，将对流扩散问题，分解为对流问题和扩散问题，再分步计算；③解析概念上的分步：将原始解析方程分解求解，例如，将曲线坐标系中的微分方程中所包含的代数项的分解处理，有时根据问题需要，也可同时采用空间概念上和物理概念上的分步方法。
2.ADI法
交替方向隐式差分格式，即所谓ADI法，该方法的实质是：将空间二维问题分成两个相互作用的一维问题，将每个计算时段一分为二，前半个时段将与y方向有关的因变量看成是已知值，x方向按一维问题采用隐式差分格式求解，后半个时段将与x方向有关的因变量看成是已知值，y方向按一维问题采用隐式差分格式求解。ADI法的实质也是分步法，目前在二维数值模拟中已经广泛应用，如国际上通用的MIKE21软件便采用了该方法。
3.预估校正法
预估校正法从广义上来讲，与ADI法一样，也是属于分步法。其基本思路是在nΔt→（n+1）Δt的求解迭代过程中，把时间步长分为两步，在前半步用稳定性较好的一阶精度格式显式计算出（n+
）Δt时刻的待求函数的近似解，然后在时间区间［nΔt，（n+1）Δt］上利用求得的近似解来构造一个二阶精度显格式来计算求得（n+1）Δt的函数值。
在预估校正法格式中，预测步和校正步中计算格式方法没必要相同，在预测步可采用比较简单的计算格式，但格式的守恒性和计算精度主要是由校正步决定，因此校正步应根据计算需要采取相应精度的格式。预测步和校正步可以根据计算的需要选择不同的格式，不同的组合可以组成不同的计算方法。
4.压力修正法
与上述方法一样，从广义上来讲，压力修正法也属于分步法。即首先假定压力场（水位场），然后求解U、V向动量方程，接着采用连续方程（质量守恒方程）进行压力修正和速度修正，不断迭代下去，直到质量源趋于无穷小为止。压力修正算法源于1972年由Patankar与Spalding提出的SINMPLE算法，SIMPLE 算法的关键步骤是在交错网格上利用连续方程和动量方程构造一个近似的压力校正方程来计算压力场并校正速度，由于压力校正方程引入了过多的近似，如略去了邻点速度修正所带来的影响，导致SIMPLE算法的收敛性与稳定性不能令人满意。为此Van Doormal（1984）基于相邻节点速度修正应与中心节点速度修正大致相等的观点，改造了SIMPLE 计算格式中导致诸多麻烦的速度较正公式，提出了所谓的相容SIMPLE算法，即SIMPLEC计算格式。计算实例表明：SIMPLEC算法数值特性的确优于SIMPLE方法。在SIMPLE 算法发展的过程中，Patankar也认识到压力修正方程用来修正速度是相当好的，而用来修正压力时相当差，因此，如果只采用压力修正方程来修正速度，而提供某些其他的方法来获得改进的压力场，可能会构成一种更为有效的计算格式，这就是后来发展的SIMPLER算法，它与SIMPLE算法的主要不同在于算法中导入了压力方程，如果用个正确的速度场进行计算，压力方程将马上给出正确的压力，而不像SIMPLE算法那样需要多次迭代才能达到压力收敛，因此，SIMPLER程序比SIMPLE程序收敛快，总的迭代次数和工作量远小于相应的SIMPLE 算法，尽管一次SIMPLER迭代需要较多的计算机时间。后来，Raithby又提出了SIMPLEX算法，Sheng等提出SIMPLET算法。
上述几种代数方程组的求解方法可以分别与有限差分法、有限元法、有限体积法结合，构成不同的计算格式。在采用分步法求解代数方程组时，可以采用任意的离散数值方法，根据研究问题的特点或模型制作者的习惯选择数值求解方法，建立数值模型。