TJU计算流体力学 第2章

1、第2章 计算流体力学的有限差分方法21差分格式用差分来替代微分、用差商来替代导数是最容易被接受的数学近似概念，随着采用符号的变化,它表示的数学意义也发生了变化，它们之间的关系为或或，下面我们来探求它的数学原理。设存在连续函数,将它用图线来表示图2-1 拟序函数序列示意图作为轴的自变量具有原点、方向和长度单位的特征，依的方向和长度单位可以依次选取,，并拟序以从小到大排列，那么对应的序列，必然存在,的函数值。同样将函数值依照的序列作拟序处理得到,或，那么与的区别在什么地方呢? 表达的仅仅是一个数值，而表达的是函数关系及对应的数值，如,当,时，而表示的是没有包含函数表达式。一种情况是首先要知道函数形

2、式才能求得函数值，另一种情况是只需知道函数值。当点取得足够致密时，序列的顺序连接可以再现的变化趋势。如果满足某种精度要求，可以认为可以完全反映的变化情况，其中表示的序列集合，由此可以看出表示的是数值变化，而表示的是函数变化， 表示的是数值变化趋势，而表示的是函数变化趋势。数值计算的基本思想就是要从离散的序列中，通过对数值计算成果的分析找出原函数的变化趋势，数值变化趋势是否反映原函数的变化趋势，取决于离散点选取的致密程度；计算数值结果的可靠性，取决于微分方程的离散、计算方法的选用和计算精度控制。可以从分布点中任意选取相邻的几个点研究离散点与原函数之间的关系，建立离散点与原函数间的变化特征。根据泰

3、勒公式，在点处连续可微（如图2-2），在的邻域内。图2-2则在邻域上可以展开为 （2-1）其中，称为展开点，称为表征点。当 取邻近的时，上式为 （2-2）当与充分接近时，的幂次越高，其在（2-2）式右端求和中所占的比重越小，直至可以忽略不计，对这样的小量之和表示为或，这里。当时，用表示，当时，用表示。式（2-2）可以整理为 （2-3）忽略项，可以看出式（2-3）的近似程度，近似程度受一次方以上各项的影响，对这样的近似式，称（2-3）为具有的一阶精度，近似式可以表示为（对单变量函数不区分偏导数或全导数）=或= （2-4）为表示方便引入符号表示沿正方向的空间步长，即。当时取为，用表示时间步长。a级

4、数组合差分法用表示一维空间函数，在，和，处用泰勒级数展开为点的级数式（如图2-3），以数字乘座标为下标表示导数，数值或座标个数表示导数阶数，用加数字表示位置，得图2-3差分网格示意图在时间坐标处，空间区间上，以为展开点，为表征点的空间差分为（）的下标表示对的导数，前的系数表示导数次数。 （2-5）在时间坐标处，空间区间上的空间差分为以为展开点，为表征点的空间差分为 （2-6）在时间坐标处，空间区间上的空间差分为以为展开点，为表征点的空间差分为（2-7）在时间坐标处，空间区间上的空间差分为以为展开点，为表征点的空间差分为（2-8）类似的，在空间坐标处，时间区间上的时间差分为以为展开点，为表征点的

5、时间差分为 一阶偏导数可以由以上四式组合而成，由（2-5）得 为在时间坐标处，空间区间上具有一阶精度的一阶偏导向前差分。由（2-6）得 为在时间坐标处，空间区间上具有一阶精度的一阶偏导向后差分。由（2-5）减（2-6）得 为在时间坐标处，空间区间上具有二阶精度的一阶偏导中心差分。由4（2-5）减（2-7）得 为在时间坐标处，空间区间上具有二阶精度的一阶偏导前差。由（2-8）减4（2-6）得 为在时间坐标处，空间区间上具有二阶精度的一阶偏导后差。二阶偏导数也可以由（2-5）至（2-8）四式组合而成，（2-5）加（2-6）得为在时间坐标处，空间区间上具有二阶精度的二阶偏导中心差分。将的一阶向前差分

6、代入式（2-7）一阶导数中得 为在时间坐标处，空间区间上具有一阶精度的二阶偏导向前差分。空间四点后差为空间五点后差为空间四点偏后差为空间五点偏后差为从以上的差分格式可以看出，不同的差分构造格式具有不同的精度阶数，也就是先得到构造的差分形式，而后才得到了它的精度。差分精度的阶数取决与略去式中的自变量差值的最低阶次。同价偏导选择点数越多，差分精度越高。分母中的步长方次决定偏导的阶次。余函数中的步长方次决定差分的精度。b待定系数差分法根据以为中心展开的泰勒级数 （2-9）得 可以看出，要保证一阶精度需要有两个结点的物理量组合，在区间上点处对一阶导数可以表示为 （2-10）展开等号右侧第一项代入（2-

7、10）得 （2-11）比较（2-11）两端系数得 得 将、代入（2-10）得。将代入（2-11）各项得，可以看出在比较系数时，具有以上的高次项被略去，所以为具有一阶精度的前差。如果希望提高精度，就要将 项消去而不是略去，就要增加可供选择的系数，增加节点个数。如 （2-12）展开等号右侧第一项和第三项得 （2-13）比较（2-13）两端系数得得 ，对和的具有以上高次项略去，所以为具有二阶精度的中心差分。归纳以上方法，要保证阶精度，就要选择点和 个待定常数来确定差分格式，这种方法称为差分格式的待定系数方法，至于选择结点的位置，则与差分方向有关，边界上的差分受到一定条件的限制。起点处用前差，终点处用

8、后差。c直接差分展开法当指定差分格式时，对高阶微分式如何用差分方法来表达，如对采用以为中心点二阶偏导中心差格式离散，先对内层离散第一项中，第二项中，先对外层离散第一项中，第二项中， 关于时间和空间混合偏导的空间中心差分时间前差式为 注时间层面用上标表示，空间层面用下标表示。 d 有权重差分 空间上取，时间上取。22计算精度无论哪一种差分格式，都存在忽略高阶项引起的误差，这种误差称为截断误差。在计算机计算的过程中，由于对实型数处理而保留的位数有限（16或32位），每计算一步在最末一位要进行四舍五入。这种累积误差称为舍入误差。从差分格式可以看出，当计算时间长度一定，越小截断误差越小，但计算次数就会

9、增加，舍入误差会增大，当较大截断误差增加，而舍入误差会减小。如果将截断误差与舍入误差之和定义为总误差，则可以得到图2-3。从图中可以看出，要提高计算精度，关键问题是要选择总误差最小的差分格式。 图2-3误差分布示意图 图2-4离散点之间函数关系23 差分的近似函数为了分析差分方程的近似程度，希望在离散的点之间建立函数关系，两点之间的函数可以用线段来近似，如：，则有 （2-14）令则，对应的取值区间，的变化范围为。对于满足任意个数结点的连续函数可以用半幅傅氏（Fourier）级数表示 （2-15）其中为从到所有结点构成的整个区间的长度。过两点能作出的可确定波动函数的最小波长为。当时（共个结点），

10、相应的，个谐波符合固定边值条件，即，。若为通过这些结点的函数，为通过这些结点的正弦、余弦曲线的振幅。根据傅里叶级数有的周期，则其中：，则有 （2-16）两相邻结点间隔波数为，波长为。是连续函数，由无穷多个点组成，要精确的用傅氏级数表示需确定无穷多个系数，需要无穷多个谐波组成。对离散函数，由有序的个结点组成，只能确定个常数，连续函数的傅氏级数为 （2-17）在结点处 （2-18）其中 （2-19）对于上式取两相邻结点的一个步长，波长，时为长期项，时为一阶谐波，为波长最长的谐波，时为二阶谐波，时为波长最短的谐波。更短的波已无法用有限的个结点确定，可以说过与点的函数与唯一确定函数的最短的可分辨波长为

11、。例如：在时，、对应函数值为、。连续函数根据两端边界条件，可以取实部或虚部。取实部得第一项波长无穷大，第二项波长，第三项波长。短于的波长的波将与对应的长波无法区分开。从已知条件中无法确定多余的常数。如给出，当，时，代入已知条件则解出，则的图像如图2-5。图2-5三结点傅氏谐波图本例可以看出如果要确定更小波长的谐波，需增加节点和已知条件。24截断误差、差分方程与微分方程近似的相容性边值问题的差分方程为 （2-20）其中，为非线性算子。初值问题的差分方程为 （2-21）或 （2-22）以下标表示边值问题差分，以上标表示初值问题差分，表示变换矩阵，表示截断误差。当时间步长和空间步长趋于零时，截断误差

12、也趋于零，差分方程与微分方程趋于一致，则称差分方程与微分方程相容，否则为不相容。使趋于零成立的最大正整数和叫做差分方程的相容阶，或称差分方程在时间上具有阶精度，空间上具有阶精度。以迁移方程为例 （2-23）用泰勒级数展开得（用下标表示微分） （2-24）其中，为时间步长，为空间步长，上式为时间显示离散，若隐式离散则为。当，时，则 （2-25）则差分方程与微分方程相容，时间上具有一阶精度，空间上具有一阶精度。25舍入误差与差分方程的稳定性求解差分方程是逐层逐点进行的，计算总数很多，舍入误差的积累会影响计算的结果的正确与否，研究这一影响过程就是讨论差分解的稳定性。研究解的稳定性可以借助于傅氏级数

13、（2-26）其中，是频率。由指数三角公式 （2-27）式（2-27）代入（2-26）得 （2-28）对二元函数 （2-29）其中，为波数，为周期，为波长。令 且得 （2-30）可以看出是第个分量的振幅，对第个分量微分有 （2-31）A 差分格式的精度分析对项作中心差分格式，略去下标，空间步长为 （2-32）上式表明用代替中心差分的截断误差是，具有的二阶精度。同时可以看出越大截断误差越大，的振幅也急剧增大。 图2-7谐波稳定性比较B 差分方程的稳定性分析以迁移方程为例分析其差分方程的稳定性。设为差分方程的解，为实际计算结果，舍入误差为，在位置处有 （2-33）计算值与差分方程解相差舍入误差为 （

14、2-34）代入（2-33）差分方程为 （2-35）差分方程解应符合原差分方程（2-33），则 （2-36）令，称为克朗条件。 （2-37）从和前系数可以看出，无论取间的任何值，在点上均有大于1的放大系数，使累积误差增加。若取傅氏级数分析，则有 （2-38）取一个谐波振幅分析 代入上式（2-37）得 （2-39）则有 （2-40）可以看出，后一时刻振幅恒大于前一时刻振幅，差分方程不稳定。若对迁移方程中的空间差分采用向后差分格式，即 （2-41）式（2-37）的误差方程为 （2-42）可以看出只要在（0，1）区间上取值，在点上均有使误差减小的效果。由傅立叶分析法可得后一时刻与前一时刻振幅比为 （2

15、-43）后一时刻振幅较前一时刻振幅恒小于1，迁移方程后差格式是稳定的。26离散化误差与差分解的收敛性如果差分方程相容，当步长趋于零时，每一步的截断误差便趋向于零，但累积的截断误差不一定趋于零，这种累积的截断误差为总截断误差，称为离散化误差。当离散化误差趋于零时，则差分解是收敛的，否则是不收敛的，这种性质称为差分解的收敛性。以迁移方程为例。设差分方程解为v，微分方程解为u，离散误差为，在位置处有 （2-44）代入差分方程（2-23）得 （2-45）截断误差为 （2-46）整理为 （2-47）即使，也不能表示，选择和时刻的最大值和，且两端取绝对值，则 最后时刻的最大截断误差为 （2-48）在中取最

16、大值，用表示，用最大值代替各项得 （2-49）可以看出，当，时，则差分方程解收敛。在实际计算，初值合理、符合物理规律的解是收敛的。例如河网初值问题。27流体力学方程的几种差分格式1） 对流方程 （2-50）中心差分格式为 （2-51）式（2-51）为不稳定格式。作修正代入式（2-51）得 （2-52）应用中心差分格式前，对结点值进行两点平均，可以使计算格式稳定。对于的齐次方程可以用显式蛙跳格式建立差分方程 （2-53）图2-8蛙跳格式示意图也可以采用全隐格式建立差分方程 （2-54）对分别求和求导，（可以将方程化为） 采用半步长离散格式（Lax-wendroff），用空间、时间半步长离散方程为

17、 （2-55）离散方程是由以下两步格式构成的第一式中为空间差分为时间坐标在时刻，空间到区间，时间差分为时间在到区间，空间坐标在处。为时间步长，为空间步长，得其中，代入得第二式中为空间差分为时间坐标在时刻，空间到区间，时间差分为时间在到区间，空间坐标在处。代入得。第三式中为空间差分为时间坐标在时刻，空间到区间，时间差分为时间在到区间，空间坐标在处。得整理得。分裂步长的典型格式为Preissmann隐格式M处的函数形式为在处，区间上的时间差分为在处，区间上的空间差分为隐格式离散和显格式离散对方程的空间差分为时间坐标在时刻，空间到区间，时间差分为时间在到区间，空间坐标在处。离散方程为隐格式。当空间差

18、分为时间坐标在时刻，空间到区间，时间差分为时间在到区间，空间坐标在处。离散方程为显格式。2） 扩散方程一维扩散方程为 （2-56）Crank-Nicolson格式（平均隐式）为 （2-57）其中，。空心菱形格式为 （2-58）跳点格式为 3） 对流扩散方程 （2-59）MacCormack格式 （2-60）其中，。4）双曲线方程 二阶波动方程为 （2-61）令 ，则等价方程组为 （2-62）交错网格差分格式为 （2-63）5） 变系数对流方程中心差分格式为 其中，可以用不同形式表现，以计算方法尽可能简单和计算模式尽可能稳定为原则，可以取 、或等形式。参考（2-52）式得 应用中心差分格式前，对

19、结点值进行两点平均，可以使计算格式稳定。28过程的稳定性和定解条件的恰当性定解问题中的两大类问题是初值问题和边值问题，一般初值是由观测或假定给出的物理量，对大范围二维或三维模型准确的给出初值是很难的。必须用近似或者假设值代替，这样初值的误差会引起方程数值解的扰动，如果这种扰动在数值解的过程中逐渐衰减或者受到一定的控制，使数值解趋于微分方程的解，则这个计算过程就是稳定的，否则是不稳定的。如对流方程 （2-64）对初值出现扰动其时刻的解受扰动影响，有，并满足方程 （2-65）由于满足方程，则扰动方程为 （2-66）设解为，当时，为扰动值。为稳定值，影响扰动振幅的大小。代入方程得，则 （2-67）可

20、见扰动以波的形式推进，振幅不变，计算过程是稳定的。边值问题取决于空间边界条件是否恰当，不同类型的定解问题对边界条件的要求不同。如对流方程 （2-68）解为 （2-69）沿特征线，应给出处的的边值，解就唯一确定。对扩散方程，应给出和处的边界条件，为的右边界。第一类边界条件为： （2-70）第二类边界条件为： （2-71）第三类边界条件为： （2-72）其中为计算域边界。（附3：特征线方法）29数值模拟的ADI差分方法对二维潮流数值模拟时一般采用ADI差分法(Alternating Direction Implicit Method)。ADI法由J.J. Le endertse于60年代首次应用于

21、潮流和扩散方程的计算。该法是一种显式-隐式交替使用的有限差分格式，其特点是将时间步长分为二等份，在前半时间步长内，沿x方向以隐式求解及，沿y方向以显式求解；在后半时间步长内，沿y方向以隐式求解及，沿x方向以显式求解。这样反复运用显-隐交替计算的方法，便能解出每个时间步长上各点的x、y方向上的流速和水深。基本方程采用笛卡儿坐标系，如图2-9所示：图2-9坐标系示意图在此坐标系下，二维潮流数学模型的基本方程为： （2-73） （2-74） （2-75）其中，为增水位；，为水深；为科氏参数；为重力加速度；、分别为x和y方向的海面风应力；是谢才系数，根据金子安雄的经验公式，其中是计算格点围点的水深平均

22、值，是曼宁系数，由试验选取。2.9.1网格的定义差分的交错网格为正方形网格，网格线分别平行于x轴和y轴，间距为，变量位置的布置如图2-9所示,表示及的位置，表示的位置，表示的位置，表示的位置。网格位置的错位定义为节点物理量平均提供了规范的计算方法。图2-9 ADI法网格示意图2.9.2 方程的离散为了简化下面的差分表达式,将使用的运算符号定义如下:，其中； 在点； 在点； ；在时间段内，(2-81)式在点上对隐式求解，对显式求解有，上式经整理得出， （2-76）(2-76)式在点上对隐式求解，对显式求解有，经简化整理得， （2-77）(2-77)式在点上对隐式求解，对显式求解有，经简化整理得，

23、 （2-78）同理，在时间段内，(2-78)式在点上对隐式求解，对显式求解有， （2-79）(2-79)式在点上对隐式求解，对显式求解有， （2-80）(2-80)式在点上对隐式求解，对显式求解有， （2-81）2.9.3差分方程的求解方法（追赶法的建立）引进一些记号，在时间段内， （2-82） （2-83） （2-84） （2-85） （2-86） （2-87）将(2-84) (2-87)代入(2-82)，有 （2-88）将(2-84) (2-87)代入(2-83)，有 （2-89）式(2-88)、(2-89)的系数和自由项都是由时刻的、所决定的，因而都是已知量，对于某一固定的，式(2-88

24、)、(2-89)组成了以和为未知量的线性代数方程组，加上左右边界条件，就构成了一个完备的线性方程组，其系数矩阵呈三对角型。 写为第列，从0到n点的矩阵=其中，具有个方程，个未知数，由矩阵两端确定边界条件，可以减少两个未知数，使方程个数与未知数个数相等。如图2-10图2-10计算点位置示意图分别从(2-88)、(2-89)中各消除一个未知数，根据三对角矩阵的追赶法，导出以下形式的递推关系式： （2-90） （2-91）方程(2-90)和(2-91)的系数也变形成为递推公式： （2-92） （2-93） （2-94） （2-95）对于每个固定的,在x轴方向上,随着的增加分别求出、及，然后,随着的减

25、小，交替地使用式(2-90)和(2-91)，分别求出和。再根据式(2-80),在y轴方向上，对于每一个固定的,随着的增加可以直出。对于在时间段内，推导完全类似，此处不再赘述。2.9.4边界的处理下面以前半个时间段内 x方向为例，说明左右开边界和闭边界的数学表达式处理。其中IS表示左边界右面的最小整数内点，IE表示右边界左面的最大整数内点。1）左端闭边界此时，由(2-90)式， 再由式(2-92)和(2-93) 可以推出，由(2-91)式， 可以推出 ， （2-96）2）左端开边界在的点上，以强制水位给出。又因不知边界处流速，作近似处理，令，设边界上流速梯度，则运动方程(2-77)可以化简为 ，而由(2-91)式又有，因此有， （2-97）3）右端闭边界此时， （2-98）没有其他特别的数学处理。4）右端开边界在的点上，以强制水位给出。又因不知边界处流速，作近似处理，令，同上述推导类似，得到 （2-99）附1：用追赶法解线性方程组线性方程组为，其中则线性方程组递推系数为，线性方程组解为，附2：平面势流在平面势流中，流动满足连续方程和无旋条件，即根据连续函数概念，可以找到这样的连续函数，则必有如果与和存在如下关系，即，则自动满足无旋条件。代入连续方程满足拉普拉斯方程，。同时的全微分有求关于的等值线簇解，即可以得到的等值曲线，同时如果已知和