赞助商链接
当前位置： >>
实用河网水流计算
实用河网水流计算王 船 海李 光 炽河海大学水资源水文系2003 年 9 月第一章 绪论河道水流运动严格地讲，其水力要素是时空均变化的水流运动，即三维非恒定的问题。 由于三维非恒定的问题在数学求解及其基本方程的理论假设上还有诸多问题，在实际计算中 常常将问题简化为二维、一维非恒定问题进行求解。对河道水流运动通常主要研究各断面流 量(或流速)，断面平均水位(不考虑横比降)随时间、断面位置的变化规律，而描述河道一维 非恒定水流运动的基本方程最早是法国科学家Ｂ． Ｓａｉｎｔ Ｖｅｎａｎｔ于 1871 年提出 的，即人们熟知的圣维南方程组。它是属于一阶双曲线型拟线性偏微分方程组，在数学上目 前无法求得其解析解。 在计算机还没有问世之前， 人们不得不将圣维南方程组简化后再求解。 如瞬态法、马斯京干法，特征河长法以及现在运用较广的扩散波法等简化方法。由于这些简 化方法是忽略圣维南方程组中的某些项，因而其适用范围也就受到限制，不是一种普遍适用 的方法。随着高速度、大容量计算机的问世及计算技术的发展，使得直接用数值方法求解圣 维南方程组成为可能，并得到飞速发展和广泛的应用，研究的问题也由一维问题发展到二维 甚至三维。第一节 洪水波的分类及其特性对于洪水波的分类，从不同的角度出发，有不同的分类方法。下面我们主要从动量方程 的量级大小进行分类。 描述河道一维水流运动的圣维南方程组如下：? ?Z ?Q B + =0 ? ? ?t ?x ? ? ?Q + ? (Qu ) + gA ?h ? gAS + gAS = 0 0 f ? ?x ? ?t ?x(1 - 1) (1 - 2)这是不考虑具有旁侧入流及动量校正系数α＝1.0 情况下的一般圣维南方程组，式中： Ｑ－流量，Z-水位，Ａ－断面面积，ｈ－水深，S0－河底比降，Sf－摩阻比降，ｕ－断面平 均流速。B B B B一、运动波对于河底比降较大的河道 gA 则动量方程(1―2)可简化为：?h ? ?Q , (Qu ), 这三项与 gAS 0 比较起来可以忽略不计， ?x ?x ?tS0 ? S f = 0或者可以写成：(1 - 3)u = c RS 0谢才公式1R 3 u= n2S0曼宁公式 流量模公式Q = K S0 ?Q ?Q +ω =0 ?t ?x类似地消去 Q 可得:联立式(1-1)、(1-3)消去 h 可得如下方程:　(1 - 4)?h ?h +ω =0 (1 - 5) ?t ?x A δu δQ δ (uA) 式中： ?= = = ηu ,ω称为波速,η=1+ 称为波速系数,由于一般情况下流速 δA u δA δA δu 随水深增加而增加,所以 ≥ 0 ,因此η≥1,这就是说,在一般情况下,波速总是大于断面平 δA均流速 u。由方程(1―4)的分析可以知道，如果ω为常数，那么方程(1―4)、(1―5)为线性 的，每一流量 Q 或水深ｈ以相同的速度ω向下游传播，因此在传播过程中不会变形。如果ω 是流量(或水深）的函数，那么在洪水波向下游传播的过程中会变形，波前变得越来越陡，波 后变得越趋平坦。当波前变成垂直时，形成了运动激波。在没有形成激波前，波形虽发生变 形，但波峰仍保持不变，没有耗散现象。由于水位流量之间呈单一关系，所以当断面最大流 量时，同时也是最高水位时刻，表示该断面的最大流量也是该时刻沿程的最大流量。所以对 运动波而言，断面最大流量、断面最高水位、沿程最大流量和沿程最大水深重合在同一个断 面上，同时出现。综上所述，运动波有以下三个重要特征： (1)它只有一族向下游的特征线， 所以下游的任何扰动不可能上溯影响到上游断面的水流 情况。 (2)不论波形传播过程中是否变形，但其波峰保特不变，没有耗散现象。 (3)当波形发生变化时，不可避免地会发生激波。二、惯性波与上面情况截然相反，如果动量方程(１―２)中摩阻项与惯性项比较起来可以忽略不 计。例如水电站突然卸负引起动力渠道中水位的波动；船闸启闭引起下游引航道内水位的波 动；深水水库中的水流波动等，由于摩阻损失相对较小，故可近似地当作惯性波。 忽略摩阻项后，并假定底坡是水平的、棱柱形河道，则动量方程式变成?u ?u ?h +u + g =0 ?t ?x ?x(1 - 6)方程(1―6)与连续方程(1―1)联立，仍属拟线性双曲线型偏微分方程，有二根实特征线：2顺特征线 顺特征方程 逆特征线 逆特征方程dx =u+c dt d(u + E) =0 dt dx =u?c dt d(u ? E) =0 dt(1 - 7) (1 - 8) (1 - 9) (1 - 10)式中 c =h gb gA ,E = ∫ dξ ，ｂ、ａ为河底以上高度ξ处的断面宽和面积。 0 a B这就是说，观测者若按ｕ＋ｃ的速度沿着波动方向运动，那么他所观测到的现象是 ｕ＋Ｅ=const，由于惯性波是不计摩阻损失，所以波动在传播过程中只有能量的转换，而没 有能量损失。由式(1―8)可见，为了保持ｕ、Ｅ二者之和不变，故当ｕ变小时，Ｅ(它反映水 深)变大；当ｕ变大时，Ｅ变小，因此流速和水深之间是互相转化，形成周期性的振荡波，湖 泊中的谐振波属于此种情况。三、扩散波动量方程(１－２)中忽略惯性项后变成：?h ? S0 + S f = 0 ?x 或B B(1 - 11) 1 ?h 1 ?h = Q0 1 ? S 0 ?x S 0 ?xZQ = K S0 1 ?(1 - 12)式中：Ｑ０为恒定状态下的流量。由于涨洪时 ?h & 0 故 Q & Q0 ?x 落洪时 ?h & 0 故 Q & Q0 ?x这就是说式(１－１２)所表示的水位 流量关系不是单值线，而是一条逆时钟的 绳套曲线，如图１―１所示，绳套大小显 然取决于附加比降ZMaxB B?h 与河底比降Ｓ ０ 之 ?xB BQ QMaxB B比值。 应用曼宁公式(或其它阻力公式)代入 (１－１２)式，并与连续方程式(1―1)联 立消去变量ｈ得关于流量 Q 的方程：图 1-1 绳套形水位流量关系3? 2Q ?Q ?Q =? 2 +ω ?x ?t ?x(1 - 13)类似地消去变量Ｑ，得到关于变量ｈ的一元偏微分方程：?h ?h ?2h =? 2 +ω ?t ?x ?x(1 - 14)式(１－１３)及(１－１４)称为对流扩散方程，它与对流方程(１－４)不同之点在于右边多? 2Q 了? －－扩散项，由于扩散项的存在所以洪水波的波峰会逐渐坦化。 ?x 2从图 1―１的绳套形水位流量关系曲线可以清楚地看到， 同一断面流量先达到最大值， 而 后水位才达到最大值。四、动力波当水位或流量在短期内的大幅度变化时，例如感潮河道中的水流运动或天然河道及人工 渠道中，因人为控制闸门的启闭而引起的水流波动。在这种情况下，动量方程式(１―２)中 的各项均不能忽略，这样一种波动称为动力波，动力波是所有波动中最复杂的，它只能用完 全的圣维南方程组来描述，因此动力波可以作为一种普遍适用的波动现象，而运动波、惯性 波和扩散波只不过是动力波的特殊情况。第二节 常用简化方法简介简化方法的要点是连续方程式严格满足,并写成差分形式:V2 +?t ?t ?t Q2 = ( I 2+ I 1 ) ? Q1 + V1 2 2 2(1 - 15)式中：Ｉ、Ｑ、Ｖ分别表示入流、出流和河槽蓄量，脚标 1 表示时段初,脚标 2 表示时段末, Δｔ表示计算时段。动力方程则用河槽蓄量Ｖ与出流量 Q 及入流量之间的某种近似关系来代 替，由于采用不同的近似关系，形成了各种各样的简化计算方法，现轭要介绍如下：一、水库调洪演算假定水库蓄水量与出流量之间存在一定的函数关系，即Ｖ＝ｆ（Ｑ） ，代入(1―15)式， 得：f (Q2 ) +?t ?t ?t Q2 = ( I 2+ I 1 ) ? Q 1 + f (Q1 ) 2 2 2上式即为水库调洪演算的基本方程，在一般情况下ｆ（Ｑ）的函数关系为非线性，难于用显 式表达，故常用图解法或试算法求解。4二、马斯京干法假定河段槽蓄量Ｖ与出流量 Q 及入流量 I 之间存在着如下的线性关系 Ｖ＝Ｋ［ｘＩ＋（１－ｘ）Ｑ］ (1－16) 式中 K 和ｘ为经验系数，且０≤ｘ≤０．５ 将关系(1―16)代入连续方程(1―15)经整理得：Q2 = C1 I 1 + C 2 I 2 + C 3 Q1式中：(1-17)C1 = C2 = C3 =Kx +α?t 2 ?t 2 ?t 2? Kx +α αK (1 ? x) ?α = K (1 ? x) +B B B B?t 2B B B B B B B BＣ１≥０，Ｃ２≥０，Ｃ３≥０，且Ｃ１＋Ｃ２＋Ｃ３＝１．0. 还有许多其它简化方法，这里不一一列举。下面讨论这些简化方法与运动波的关系，假 定运动波基本方程?Q ?Q =0 +ω ?x ?t(1-18)式中的波速ω用流量变幅范围内的平均波速〈ω〉代替，并按如下差分格式? ?Q θ ( I 2 ? I 1 ) + (1 ? θ )(Q2 ? Q1 ) ≈ ? ? ?t ?t ? ? ?Q ≈ Q2 + Q1 ? I 2 ? I 1 ? 2?x ? ?x式中：θ为权重系数，令Ｋ＝I2B B△XQ2B B?x &ω &△t －表示 △t/2 I1B B河段传播时间，并将上式代入到式(1―18) 加以整理后得与式(1―17)同样的形式，只 不过以θ代替了马斯京干法的ｘ而已。由此 可见马斯京干法只不过是运动波的一种差 分格式而已。但运动波是没有衰减的，而马 斯京干法却为什么有衰减特性呢?对这个问 题，要通过以后章节的学习才能作出回答。 若令θ＝０，则有Ｃ１＝Ｃ２＝B B B BQ1B Bθ△X 图 1-2K ? ?t 2 ?t / 2 及Ｃ３＝ ，代入式(1―17)得： K + ?t / 2 K + ?t 2B B5(K +?t ?t ?t )Q2 = ( I 2 + I 1 ) + ( K ? )Q1 2 2 2(1-19)若令Ｖ＝ＫＱ，由（1-15）也可得(１―19)式，即(１―19)式表示线性水库的调洪演算。 从上述讨论可见，要精确了解实际流动，必须求解动力波方程。水文学上较常用的简化 方法马斯京干法是在出流与槽蓄量单一函数关系的假定下导出的，是运动波的差分解。这种 假定在流域上游的水流运动与实际基本符合，有足够的精度。在流域下游，特别在平原河口 地区，流动受上游来流和下游潮位的联合作用，由于受下游水位的顶托，流动不是自由出流， 上述假定不复存在，实际流动的模拟应该由动力波方程来描述，即必须直接求解圣维南方程 组。圣维南方程组的求解只能通过数值方法。常见的数值方法有，有限差分法、有限单元法 和有限分析法，求解一维水流最常用的是有限差分法。随着计算机技术的发展，计算技术的 提高，模型化已成为当今发展的主流。水文学的发展同样越来越模型化。随着社会的发展， 人类文明的进步，人们对于水资源的利用和管理日益受到重视，特别在发达的大中城市，地 处河口平原地区， 水资源的供需矛盾日益突出， 管理水平必须提高才能满足社会发展的需要， 为管理决策提供技术支持的河网水流模拟更显得日益重要。河网水流模拟是平原地区水文学 必不可少的基础，本文主要从实用角度出发，介绍河网水流的基本计算方法。6第二章 明渠一维非恒定流的基本方程第一节 方程推导的基本条件一、基本条件研究的河段为非棱柱形河道,且有单位河长上分布的旁侧流量为ｑ ｌ,其入流方向与流向 的夹角为θ。以上是实际问题中经常遇到的普遍情况。 基本假设： 1、假设所研究的是定床情况，即假设河床高程只与位置有关，而与时间无关。 2、假设断面水面线沿宽度方向上水平，这是推导一维方程的一个基本假设。　 3、假设研究的是浅水问题，可假设压力在垂线上满足静水压力分布规律。　 4、ρ＝ｃｏｎｓｔ，即水的密度为常数。 5、 假设所研究的河床底坡很小。 以致于可以忽略水深与Ｚ－Ｚｄ 的差别。 Z 由图 （２－１） 可见: Ｚ －Ｚｄ＝ｈｃｏｓα，当底坡较 h 小时， ｃｏｓα≈1.0,即认为:h ≈Z-Zd 以及有 ｓｉｎα≈ Ｓ０，其中：Ｓ０为底坡。 6、 非恒定摩阻公式可借用恒 Zd 定流摩阻公式。B B B B B B B B B B B Bα二、基本定律图 2―1 水深示意图质量守恒定律：假设有一控制表面系统包围一个控制体积，并使控制体积的内部和外部 是被唯一定义的。那么，通过控制表面由外部进入控制体积的流体净质量，等于这一体积所 增加的净质量。 动量守恒定律：作用于控制体积的冲量与通过控制表面的进出流体净动量的矢量和等于 这一体积内的动量增量。7第二节 基本方程之推导一、连续方程的推导质量守恒定律：通过控制面进到控制体的质量＝控制体内质量增量 取图２－２所示的断 1 面１－１之间水体为控制 ? 体。在Δｔ时间内，通过 ( A?xρ ) ? ?t ?t 上游断面１－１流入到控 2 制体的质量为： ＱρΔｔ 通过下游断面２－２ 流出控制体的质量为： ρQQρ?t +? (Qρ?t ) ? ?x ?xt+△t故在Δｔ时间内， 通过上、 下游断面进、出质量之差 为：tρQ +?ρQ ?x ?x?? (Qρ?t ) ? ?x ?x△X 1 图 2-2 质量守恒原理示意图 2由于旁侧入流流入到 控制体的质量为：ρq l ?x?t控制体积内质量的增 量为：? (A?xρ) ? ?t ?t由质量守恒定律可得：? ( A?xρ ) ? ( ρQ?t ) ? ?t = ? ? ?x + ρql ? ?x ? ?t ?t ?x化简后得： 　　?A ?Q + = ql ?t ?x（２－１）上式简化过程中用到了ρ为常数的假设。式中：Ａ为过水断面面积，Ｑ为断面流量。二、动量方程的推导动量方程的推导主要依据是动量守恒定律。由于动量是一矢量，下面的推导是建立水流 流动方向的动量方程。 通过控制面流进到控制体内的动量＋作用于控制体外力的冲量＝控制体积内动量的增 量。8１、控制体积的动量增量 如图 2-3 所示，取断面１与断面２之间的水体为控制体，该控制体内的沿流向方向的动 量为: 　　Ｍ＝ρΔｘＱ 经过Δｔ时间后,这块水体的动量发生了变化,其增量为: 　　 ?M =?M ? (ρ?xQ) ? ?t = ? ?t ?t ?t1（２－２）? ( ρ?xQ ) ? ?t ?t2αρQu t+△ttαρQu +?αρQu ?x ?x△X 1 图 2－３动量守恒原理示意图 ２、通过控制面流入到控制体的动量 在Δｔ时间内由断面 1 流进控制体的动量为: αρQu?t 。α称之为动量校正系数。对于一 般情况可以认为α≈１，当断面流速分布非常不均匀时，α就显得重要。 由下游断面２流出的动量为： αρQu?t + 2? (αρQu?t ) ? ?x ?x（２－３）因此可得上、下游进、出控制体积的动量差为： 　　 ?? (αρQu?t ) ? ?x ?x在Δｔ时间内，由于旁侧入流，流入到控制体的动量为： 　　 ρq l ? ?t ? ?x ? VxB B（２－４）式中：Ｖｘ为旁侧入流流速在水流方向上的分量。 通过上面的推导，可以得到在Δｔ时间内通过控制面流入到控制体的动量为： 　　 ?? (αρQu?t ) ? ?x + ρq l Vx ?x?t ?x（２－５）３、作用于控制体上力的冲量 作用于断面１－２之间的控制体上的力有：重力、摩阻力及水压力。9a．重力 控制体的体积为Ａ?Δｘ，故其重力为ρｇＡ?Δｘ，由于我们所建的是沿水 流方向上的动量方程，故重力在水流方向上的分量为： 　　 ρｇＡ?Δｘ?ｓｉｎθ （２－６） 式中：θ为河底与水平面的夹角，由于前面的假设 5，故可近似地有重力在水流方向上的分 量为： 　　 ρｇＡ?ΔｘＳ０ （２－７） 式中：Ｓ０为河道的底坡，可表达成：B B B B　　S0 = ??Z d ?x（２－８）ｂ．摩阻力 由假设６知，在非恒定流情况下摩阻损失与恒定流情况下差别不大，仍可用 曼宁公式、谢才公式或流量模数公式： 　 Sf =n 2u 2 R4 3曼宁公式（２－９）Sf =u2 c2R谢才公式(２－１０）Q2 Sf = 2 K流量模数公式(２－１１）式中：ｎ为糙率，ｃ谢才系数，Ｋ流量模数，Ｒ水力半径，Ｓｆ为摩阻比降。由此可得出作 用于控制体上的摩阻力为：B B　 ? ρgA ? ?x ? S f ? sgn(u ) 　（２－１２）0式中：sgn(u)为一符号函数，如果将Ｓｆ?sgn(u)合并统一成 S f ，则式(2-9)～(2-11)成为如B B下形式： 　 Sf =0n2 u u R4 3曼宁公式（２－13）S0 f =uu c2R QQ K2谢才公式(２－１4）S0 f =流量模数公式（２－１5）下面为了书写方便，省略上标“0”,故式(2－12)可改写成：　 　　 ? ρgA ? ?x ? S f (２－１６）ｃ．压力 压力对断面１－２之间水体的作用，分为两部分： １．断面压力的分析 由图 2-5 可知，微元面积 bdξ上的压强为：ρg(h-ξ)，其上的压力为：ρgb(h-ξ)dξ, 故可得断面１上的压力为：10　 FaPP+?P ?x ?xFa图 2―4.控制体上压力分布图Bδξhξ图 2―5．断面压力积分示意图P=∫h ( x ,t )0ρgb( x , ξ, t )(h ( x , t ) ? ξ)dξ断面１与断面２的压力差为：?P = ?　　?P ? ?x ?x h ( x ,t ) ?h ?b( x , ξ, t ) = ?ρgA ? ?x ? ∫ ρg (h ( x, t ) ? ξ) ? ?xdξ 0 ?x ?x(2-17)２.侧壁上的压力分析 对于作用于侧壁上压力，没有必要求出其总压力，只需求出其在水流方向上的分量即可。 由水力学中的有关曲面的静水总压力的有关定理知：作用于曲面上的静水压力在某一方向上 的分力等于该曲面在垂直于该方向上的投影面积的静水总压力。即：作用于断面 1-2 之间的112 2 1 1图 2―6.侧壁压力分力计算示意图 侧壁压力在水流方向上的分量等于断面１与断面２的差值部分面的压力，如图 2-6 所示的阴 影部分面积的静水总压力，其大小为：R=∫h ( x ,t )0ρg(h ( x , t ) ? ξ) ?h ? ?x ?x?b( x, ξ, t ) ? ?xdξ ?x（2-18）综合式（２－１７） 、 （２－１８）可得作用于控制体上的总压力在水流方向上的分量为： 　　 ? ρgA （2-19）由式（２－７） 、 （２－１６） 、 （２－１９）可得作用于控制体上的总外力为：F = ρgA ? ?xS 0 ? ρgA ? ?xS f ? ρgA ? ?x?h ?x（2-20）由动量守恒定律及式（２－２） 、 （２－５） 、 （２－２０）可得：　　? (ρ?xQ) ? ? ?t = ? (αρQu?t ) ? ?x + ρq l Vx ? ?t?x ?t ?x + (ρgA ? ?xS 0 ? ρgA ? ?xS f ? ρgA ? ?x ?h )?t ?x两边同除以ρΔｘΔｔ，并化简可得如下动量方程： 　　?Q ? ?h + (αQu ) + gA = gA(S 0 ? S f ) + q l Vx ?t ?x ?x QQ ?Q ? ?Z + (αQu ) + gA + gA 2 = q l Vx ?t ?x ?x K（2-21）上式也可写成如下等价形式：　 　　 （2-22）方程（２－１） 、 （２－２１）即为一般情况下的明渠非恒定流的基本方程，又称为圣维 南方程组，一般情况下认为Ｖｘ为零。B B12第三节 方程的其它形式及问题讨论一、方程的其它形式从方程（２－１） 、 （２－２１）可见，方程中变量有Ｑ、ｕ、Ａ、ｈ等变量，而这些变 量并不相互独立，具有一定关联，因此需对上述方程作一些技术处理，使之成为一组便于求 解的方程组。　 １、以水位Ｚ、流量Ｑ为因变量 　　 B?Z ?Q + = ql ?t ?x（2-23）　Q2 ?Q ? ?Z + (α ) + gA = ?gASf ?t ?x A ?x ?Z ?Q + = ql ?t ?x（2-24）２、以Ｑ、ｈ为因变量 　　 B （2-25）　?Q ? Q2 ?h + (α ) + gA = gA(S 0 ? S f ) ?x ?t ?x A ?A ?Q + = ql ?t ?x ?Q ? Q2 + (α + P) = gA(S 0 ? S f ) + R A ?t ?x（2-26）３、以Ｑ、Ａ为因变量 　　 （2-27）　（2-28）式中： P = g∫h ( x ,t )0(h ( x , t ) ? ξ)b(ξ, x , t )dξ 　　（2-29） （2-30）R = g∫h ( x ,t )0( h ( x , t ) ? ξ)?b(ξ, x, t ) dξ 　　 ?xＰ中所含积分的几何意义是过断面相对于水面高程的面积矩，其中ξ为相对河底的高度，ｂ 为相对于ξ的过水断面宽。可详见式（2-17）之推导。ρＲ表示单位长度河床对水体的反作 用力在水流方向上的投影，可详见式(2-18)的推导，它是由断面沿程变化而引起的。对棱柱 形明渠 R=0,由式(2-17)知：?P ?h ? R = gA (2-31) ?x ?x ?h 可知 gA 表示断面水压差与单位长度河床反作用力之代数和，系由于水深变化而引起的， ?x　　 否则该两力相互抵消。 以下各式的推导是假设α＝１的情况下推导出的。 ４、以ｕ、ｈ为因变量13?h A ?u ?h u ?A + +u + ?t B ?x ?x B ?x=hqt B(2-32)q ?u ?u ?h +u +g = g (S 0 ? S f ) ? t u ?t ?x ?x A因为Ａ＝Ａ（ｈ（ｘ，ｔ） ，ｘ） ，故： 　　(2-33)?A ?h ?A =B + ?x ?x ?x h(2-34)从而式（２－３２）可写成：　 　　?h A ?u u ?A qt + + = ?t B ?x B ?x B(2-35)５、以ｕ、Ｚ为因变量 　　q ?Z A ?u ?Z u ?A + +u + = t ?t B ?x ?x B ?x z B q ?u ?u ?Z +u +g = ? gS f ? t u ?t ?x ?x A(2-36)(2-37)上述几组方程中，可根据具体的问题及其所采用的格式选用相应形式的方程，在求解天 然河网水流计算中，常采用方程组（2-23） 、 （2-24）或方程组（2-25） 、 （2-26） 。二、具有漫洪滩地河道问题的处理在推导动量方程时，我们引进了动量校正系数α，其表达式为：　 　　 α =A u 2 dA 2 ∫A Q(2-38)该系数反映河道断面流速分布均匀性的系数，对于一个主槽的河道情况下，可令 α≈1.0， 但当河道有漫滩地时，其动量校正系数非考虑不可，这主要是由于滩地上的流速远比主槽流 速小得多，因而断面的流速分布非常不均匀，使得α达到相当大，尤其是当水流刚上滩时尤 为严重。下面介绍在一定假定的情况下，推导动量校正系数的公式，使其象断面资料一样变 成预先可求得的资料。 如图（２－７）所示,整个断面分成三部分,其过水面积分别为Ａ１、Ａ２、Ａ３，流量分别 为Ｑ１、Ｑ２、Ｑ３。 断面总流量:　Ｑ＝Ｑ１＋Ｑ２＋Ｑ３ 总过水积为:　Ａ＝Ａ１＋Ａ２＋Ａ３B B B B B B B B B B B B B B B B B B B B B B B B断面平均流速: u =Q Q1 + Q 2 + Q 3 = A A1 + A 2 + A 314断面每部分之平均流速: u i =Qi Ai(i = 1,2,3)A1A2A3图 2―7 洪漫滩地流动示意图 BTB BB图 2―8 调蓄宽度示意图 单位时间内通过整个断面的动量应为:∑ ρQ ui =1 i3i，它比按照断面平均流速计算的动量ρＱｕ要大,其比例系数即为动量校正系数。 假定滩地和主槽中水流的摩阻比降是相同的，所以有：Sf =故：Q1 Q 2 Q 3 Q = = = K1 K 2 K 3 K（２－３９）15α=∑ ρQ i u ii =13ρQuQ i2 ∑A 3 K2 A = i =1 2 i = 3 ?∑ i A Q (∑ K i ) 2 i =1 i A i =13（２－４０）因为Ａｉ、Ｋｉ为水深（或水位）及断面位置的函数，故此值亦象其它河道断面资料一样， 可以预先整理成α＝α（Ｚ，ｘ）作为原始基本资料。 当滩地只起调蓄作用，不起输送水量作用时，也就是滩地上流速很小，动量方程式起作 用的是主河槽部分，所以凡是动量方程式中断面积Ａ及河宽 B 均按主槽部分计算。在连续方 程中出现的Ｂ和Ａ都应包括滩地在内的全部河宽和过水面积。 所以当滩地只起水量调蓄作用， 不起输水作用时，其基本方程组如下：B B B B? ?Z ?Q BT + = ql ? ? ?t ?x 　　　 ? 2 ? ?Q + ? ( Q ) + gA ?Z = ? gAS f ? ?x ? ?t ?x A　（２－４1）第四节 定解条件及定解问题非恒定问题必须给出初始条件。至于空间域,由于用有限差分法求解时总是取有限空间, t C Dx a B 图 2―9 特征线示意图 必须给出边界条件。 初始条件与边界条件是确定微分方程解的必不可少的条件,合称为定解条 件。定解条件少给了,就无法定解，即所谓欠定;多给了可能产生矛盾,即所谓过定;只有给出 适当个数的定解条件才能求解,即所谓适定。现在考虑对流方程:　16AB?u ?u +C =0 （２－４2） ?t ?x dx ?u dx ?u 若式中的系数 C = 　,则（２－４2）成了： + = 0 　　 dt ?t dt ?x du 或写成：　　　 =0 dt dx 1 不难看出, C = 在ｘ-ｔ坐标系中是一组斜率 tgθ = 　的直线 （图 2―９） ,称为特征线。 dt C du 在此线上 = 0 ,即ｕ的值沿特征线不变。 dt　　 设问题的定解域为ａ≤ｘ≤ｂ,０≤ｔ，如图（２－9）中ｘ轴上方ＡＣ、ＢＤ间的区域。 对这种问题如何给出定解条件,是与特征线有关的。若Ｃ≥０　,则当ＡＢ和ＡＣ上的函数给 定时,ＢＤ上的函数值也随之而定。所以 C≥０时,只需给定初值 u(x,0)=F(x)和边值 u(a,t)=G(t),就适定了。这里Ｆ（ｘ） 、Ｇ（ｔ）表示已知函数。同样道理,若Ｃ≤０　,则只 需给定初值ｕ（ｘ，０）＝Ｆ（ｘ）和边值ｕ（ｂ，ｔ）＝Ｇ（ｔ）　,Ｆ（ｘ） 、Ｇ（ｔ） 为已知函数,即给定 AB 和 BD 上的函数,就适定了。 从上述讨论可看出一个规律:将特征线视为有向线段,其方向与ｔ增加时的方向一致。定 解域边界(包括初始时刻)上的某点是否要给出定解条件,可从该点处作特征线,若其方向指向 定解域内,则需给出定解条件,否则不必给。 上述方程只有一组特征线。一般双曲型方程组可以有多组特征线。对此可有类似的结论： 在计算域的边界(请注意， 这里没有必要区别什么是初始条件边界， 什么是边值条件边界)处， 边界条件的个数就等于穿越边界而进入定解域的特征线的条数。 根据上述结论，可导出一维圣维南方程组适定的定解条件，具体如下： 以方程组(2―23)、(2-24)为例，为方便起见,设ｑｌ＝０及α＝1.0，有如下特征线与特 征方程： 沿顺特征线:B B? dx ? dt = u + C ? ? 2 ? DQ + B(? u + C) DZ = [BS0 + ?A h ] Q ? gASf ? Dt ?x A2 ? Dt沿逆特征线:（２－４3）? dx ? dt = u ? C ? ? 2 ? DQ + B(? u ? C) DZ = [BS0 + ?A h ] Q ? gASf ? Dt ?x A2 ? Dt式中 C 为重力波速 C =（２－４4）gA B设整个计算域为：０≤ｔ≤Ｔ，０≤ｘ≤Ｌ，由于水流分为缓流和急流,因而其定解条件 也不相同，下面分别进行讨论（设ｕ≥０的情况下） ： 1.缓流17由于是缓流，因而水流的佛汝德数 Fr =u gA B=u ≤ 1.0 ，故域内一点处的两条特征线 C分别指向上、下游。如图(2―10)所示。故其定解条件为： ａ．初始条件：在边 AB 上，两条特征线均指向计算域内，故在此边上必须给出两个条件， 即： 　　 Ｑ（ｘ，０）＝已知， Ｚ（ｘ，０）＝已知 ｂ．上边值条件：在边 AC 上，只有一条特征线指向计算域内，故在此边上只须给出一个 条件，即： 　　 Ｑ（ｔ，０）＝已知， 或 Ｚ（ｔ，０）＝已知 ｃ．下边值条件:在边 BD 上，只有一条特征线指向计算域内，故在此边上只须给出一个 条件，即： 　　 Ｑ（ｔ，L）＝已知， 或：Ｚ（ｔ，L）＝已知 或：　 F(Z，Q)｜ｘ＝Ｌ＝０， （水位流量关系） 在边 CD 上，没有一条特征线指向计算域内，故在此边上无须给出条件。B Bt C DTx A L 图 2―10 缓流情况下边界处特征线方向示意图 ２．急流 由于是急流，因而水流的佛汝德数 Fr = Bu gA B=u ≥ 1.0 ，故域内一点处的两条特征线 C均指向下游。如图(2―11)所示。故其定解条件为： ａ．初始条件：在边 AB 上，两条特征线均指向计算域内，故在此边上必须给出两个条件，18即： 　　 Ｑ（ｘ，０）＝已知， Ｚ（ｘ，０）＝已知 ｂ．上边值条件：在边 AC 上，两条特征线均指向计算域内，故在此边上必须给出两个条 件,即： 　　 Ｑ（ｔ，０）＝已知， 和 Ｚ（ｔ，０）＝已知 在边 BD、CD 上，没有特征线指向计算域内，故在这两条边上无须给出任何条件。t C DTx A L 图 2―11 急流情况下边界处特征线方向示意图 B19第三章 有限差分的基本理论第一节 基本概念从上一章讨论可见， 所求解的基本方程是一双曲型偏微分方程组,该方程组用现有数学理 论目前得不出解析解,除了某些特殊情况及简化情况下可写出其准确解 （如棱柱形恒定流情况 等） ，一般都采用数值求解。差分方法是数值求解偏微分方程组的一种方法,下面以一维非恒 定问题为例,说明差分方法的一些基本概念。 考虑双曲型一维对流方程:?u ?u +C =0 ?t ?x(3 - 1)Ｃ为波速,下面为了研究方便取Ｃ为常数。 问题求解域为ｘ－ｔ的上半平面。 在上半平面上画出两族平行于坐标轴的直线,把求解域 分成矩形网格。网格线的交点称为节点,ｘ方向上网格线之间的距离Δｘ称为空间步长,ｔ轴 方向上网格线之间的距离Δt 称为时间步长,见图３－１。这样两族网格线可记为： t2△t △t X -2△x 0 2△x图３―1 网格剖分 ｘ＝ｘｉ＝ｉ?Δｘ ｉ＝０，±１，±２，±３，… ｔ＝ｔｊ＝ｊ?Δｔ ｊ＝０，１，２，３，… 网格节点（ｘｉ，ｔｊ）可记为（ｉ，ｊ） ，则节点处的函数值可记为：B B B B B B B Buij ＝ｕ（ｘｉ，ｔｊ）B B B B如果网格剖分使得每一空间步长、时间步长均相等，则称该网格为一均匀网格，否则称 之为非均匀网格。20数值解主要是求解节点上的末知变量的数值，利用有限的节点上的值来代替整个求解域 内的连续函数值。如果节点上的函数值已求出，则节点以外的函数值可利用节点上的值，采 用某种方式插值求解，当然这样求得的解与本身微分方程的真解会有差别，该差别称之为误 差，本章的主要研究如何数值求解以及其误差究竟如何？第二节 偏导数的差商近似一、差分、差商的基本概念设有解析函数 u = f ( x ) 从微分学知道函数ｕ对ｘ的导数可定义如下：?u du u ( x + ?x) ? u ( x) = lim = lim dx ?x→0 ?x ?x→0 ?xΔｕ、Δｘ分别称为函数及自变量的差分，(3 - 2)?u 为函数对自变量的差商。 ?x在上述的导数定义中Δｘ是以任意方式趋近于零的，因而Δｘ是可正可负的。在差分方 法中，Δｘ总是取某一小的正数。这样在差分方法中，可以有如下三种基本形式： 向前差分：　Δｕ＝ｕ(ｘ＋Δｘ)－ｕ(ｘ) 向后差分：　Δｕ＝ｕ(ｘ)－ｕ(ｘ－Δｘ) 中心差分：　Δｕ＝ｕ(ｘ＋Δｘ)－ｕ(ｘ－Δｘ) 上面介绍的是一阶导数，对应的差分称为一阶差分。对一阶差分再作一阶差分，所得到 ２ 的称之为二阶差分，记为Δ ｕ。以向前差分为例有：P P?2 u = ?(?u ) = ?[u ( x + ?x) ? u ( x)] = [u ( x + 2?x) ? u ( x + ?x)] ? [u ( x + ?x) ? u ( x)] = u(x + 2?x) - 2u(x + ?x) + u(x)依此类推，任何阶差分都可以由其低一阶的差分得到：?n u = ?(?n?1u ) = ?(?(?n?2 u )) = Λ Λ Λ = ?{? Λ Λ [?(u(x + ?x) - u(x))]}ｎ阶的向后差分、中心差分的形式类似。 函数的差分与自变量的差分之比，即为函数对自变量的差商，如一阶向前差商为：?u u ( x + ?x) ? u ( x) = ?x ?x ?u u ( x) ? u ( x ? ?x) 一阶向后差商为： = ?x ?x ?u u ( x + ?x) ? u ( x ? ?x) 一阶中心差商为： = ?x 2?x二阶差商多取中心式，即：?2 u u ( x + ?x) ? 2u ( x) + u ( x ? ?x) = (?x) 2 ?x 221以上讨论的是一元函数的差分与差商，对多元函数的差分与差商可以类推。二、偏导数的差商近似及其构造１．　Ｔａｙｌｏｒ展开法 由导数和差商的定义知道，当自变量的差分（增量）趋近于零时，就可由差商得到导数。 因此在差分方法中用差商代替偏微分方程中的偏导数。差商与导数之间的误差表明差商逼近 导数的程度，称为逼近误差。由函数的Ｔａｙｌｏｒ展开，可以得到逼近误差相对于自变量 差分（增量）的量级，称为用差商代替导数的精度，简称为差商的精度。 通过上述的分析对图３－１中点（ｉ，ｊ）的空间一阶偏导数的差商近似 有如下三种 形式：?u j uij+1 ? uij )i ≈ ?x ?x j ?u j ui ? uij?1 ( )i ≈ ?x ?x j u ? uij?1 ?u ( ) ij ≈ i +1 ?x 2?x (向前差商 向后差商 中心差商(3 - 3) (3 - 4) (3 - 5)通过对差商近似点（ｉ，ｊ）的Ｔａｙｌｏｒ展开，可以分析三种差商对偏导近似的精 度：uj i +1uij?1?u j 1 ? 2u j 1 ? 3u j 2 = ui + ( ) i ? ?x + ( 2 ) i ? ?x + ( 3 ) i ? ?x 3 + Ο( ?x 4 ) ?x 2 ?x 6 ?x 2 ?u 1 ? u 1 ? 3u = uij ? ( ) ij ? ?x + ( 2 ) ij ? ?x 2 ? ( 3 ) ij ? ?x 3 + Ο( ?x 4 ) ?x 2 ?x 6 ?xj(3 - 6) (3 - 7)由式（３－６）得:(?u j uij+1 ? uij 1 ? 2 u j )i = ? ( 2 ) i ? ?x + Ο( ?x 2 ) ?x ?x 2 ?x(3 - 8)由式（３－7）得由式得:(?u j uij ? uij?1 1 ? 2 u j )i = + ( 2 ) i ? ?x + Ο( ?x 2 ) ?x ?x 2 ?x(3 - 9)式(3-6)减式(3-7)并化简得:?u j u ij+1 ? u ij?1 1 ? 3 u j ? ( 3 ) i ? ?x 2 + Ο(?x 3 ) ( )i = ?x 2?x 6 ?x(3 - 10)由此可见由向前、向后、中心差商来近似偏导的精度分并为：一阶、一阶与二阶。式(３－６） 与式（３－７）相加，并化简可得：(? 2 u j uij+1 ? 2uij + uij?1 )i = + Ο( ?x 2 ) ?x 2 ?x 2(3 - 11)由此可见二阶偏导数的差商近似表达式为：22(? 2 u j u ij+1 ? 2u ij + u ij?1 )i ≈ ?x 2 ?x 2(3 - 12)比较式（３－１２）与式（３－１１）可见其近似精度为二阶。 同样也可以写出时间一阶偏导的差商近似表达式及其精度。对于上述的中心差商，当网 格为非均匀网格时，可能就不是中心差商了。 在实际计算时，空间区域往往是有限的。故在边界处的偏导的差商近似的构造在计算时 显得相当重要，下面举例说明如何构造边界处的偏导数近似。 j 0 1 2 ……… 图 3―2 例 1、构造边界点处（０，ｊ）对ｘ的一阶偏导的二阶差商近似。对点（０，ｊ）进行Ｔ ａｙｌｏｒ展开得： n?u j 1 ? 2u 1 ? 3u ) 0 ? ?x + ( 2 ) 0j ? ?x 2 + ( 3 ) 0j ? ?x 3 + Ο( ?x 4 ) ?x 2 ?x 6 ?x 2 ?u j 1 ? u j 1 ? 3u j j j 2 u2 = u0 + ( ) 0 ? 2?x + ( 2 ) 0 ? ( 2?x ) + ( 3 ) 0 ? ( 2?x ) 3 + Ο( ?x 4 ) ?x 2 ?x 6 ?xu1j = u0j + ((3 - 13) (3 - 14)式(3―13)×α１＋式（3－１４）×α２得：B B B B?u j ? 2 u j ?x 2 α1 × u + α 2 × u = (α1 + α 2 )u + (α1 + 2α 2 )( ) 0 ? ?x + (α1 + 4α 2 )( 2 ) 0 2 ?x ?x 3 (α + 8α 2 ) ? u j + 1 ( 3 ) 0 ? ?x 3 + Ο( ?x 4 ) (3 - 15) ?x 6j 1 j 2 j 0由于我们要构造一阶偏导数的二阶精度的差商近似，故式（３－１５）中必须有： α１＋２α２＝１ 与 α１＋４α２＝０ 解得： α１＝２ α２＝－０．５ 则可得：B B B B B B B B B B B B?u j 4u1j ? 3u0j ? u2j ? 3u j ?x 2 + 3 0? + Ο( ?x 3 ) 0= 2 ?x ?x 3 ?x由式(3-16)可见(3 - 16)?u ?xj 0的二阶精度的差商近似为：?u j 4u1j ? 3u0j ? u2j 0= ?x 2?x由（３－１２）可见不能用该式构造(3 - 17)? 2u ?x 2j 0在边界点的差商近似，我们仍然可以参照上述方法构造其差商近似，由式（３－１５）中可令： α１＋２α２＝０ 与 α１＋４α２＝２B B B B B B B B23解得： α１＝－２ α２＝１ 代入到式（３－１５）并化简得：B B B B? 2 u j u 0j ? 2u1j + u 2j + Ο(?x) 0= ?x 2 ?x 2由上式可见其近似的精度为一阶，若要其二阶差商近似，还要对点（３，ｊ）展开的近似式：u 3j = u 0j + (1 ? 2u 1 ? 3u ?u j ) 0 ? 3?x + ( 2 ) 0j ? (3?x) 2 + ( 3 ) 0j ? (3?x) 3 + Ο(?x 4 ) 2 ?x 6 ?x ?xB B B B B B(3 - 18)类似的处理，式（３－１３）×α１＋式（３－１４）×α２＋式（３－１８）×α３得：α 1u1j + α 2 u 2j + α 3 u 3j = (α 1 + α 2 + α 3 )u 0j + (α 1 + 2α 2 + 3α 3 )(?u j ) 0 ? ?x ?x ? 2 u ?x 2 ? 3 u ?x 3 + (α 1 + 4α 2 + 9α 3 )( 2 ) 0j + (α 1 + 8α 2 + 27α 3 )( 3 ) 0j + Ο(?x 4 ) 2 6 ?x ?xB B B B B B令：α１＋２α２＋３α３＝０ α１＋４α２＋９α３＝2 α１＋８α２＋２７α３＝０ 解得：α１＝－５ α２＝４ 故有：B B B B B B B B B B B B B B B Bα３＝－１B B? 2 u j 2u0j ? 5u1j + 4u2j ? u3j + Ο ( ?x 2 ) 0= ?x 2 ?x 2由此得? 2u ?x 2j 0在边界点处的二阶精度的差商近似为：2? 2 u j 2u0j ? 5u1j + 4u2j ? u3j 0= ?x 2 ?x 2(3 - 19)上述介绍的Ｔａｙｌｏｒ展开法构造偏导数的差商近似，该法是最常用的方法，它简便 不包含物理意义，但可以得到差商近似的精度。 ２．多项式插值法 用多项式插值法把待求函数表示成含待定系数的解析函数，由节点函数值确定该系数， 然后对此函数求偏导数，得到逼近偏导数的差商表达式。 以中心差商式（３－４）为例，说明其偏导数的中心差商的构造过程。选择在ｊ时刻上 ｉ－１、ｉ、ｉ＋１三个节点，设在此区间函数ｕ可用抛物插值公式来近似： ２ ｕ(ｘ，ｔｊ)＝ａ＋ｂｘ＋ｃｘ 式中：ａ、ｂ、ｃ为待定系数。为了方便，设原点ｘ＝０在点ｉ的位置上，则可写出三个节 上的函数有：B B P P?uij?1 = a ? b?x + c?x 2 ? j ?u i = a ? u j = a + b? x + c ? x 2 ? i +1解出待定系数：24? uij+1 ? uij?1 b = ? ? 2 ?x ? j j j ?c = ui +1 ? 2ui + ui ?1 ? 2 ?x ?微分函数?u ，并计算在点ｉ的值得： ?x(3 - 20) (3 - 21)uij+1 ? uij?1 ?u j ( 2 ) = b + c ? x = i x =0 2?x ?x uij+1 ? 2uij + uij?1 ? 2u = 2c = ?x 2 ?x 2得到的一阶中心差商和二阶中心差商公式和Ｔａｙｌｏｒ展开得到的结果相同。如果选用一 次插值和不同的插值点，则可得一阶向前或向后差商。 用高阶多项式插值可得到高阶差商表达式。除一、二次多项式外，所得到的表达式与高 阶Ｔａｙｌｏｒ展开得到的结果并不相同，这时多项式展开的误差仍需用Ｔａｙｌｏｒ级数 展开来检验。此外高阶多项式插值具有龙格不稳定性，使得插值对计算误差十分敏感。多项 式插值法在计算流体力学中多用于处理边界处的差商近似。 关于偏导数的差商近似还有其它多种方法，但最终均需用Ｔａｙｌｏｒ展开来计算其近 似的误差，因此在实际计算中通常均用Ｔａｙｌｏｒ展开法来构造，因为此法在构造差商近 似的同时还得出了其近似的误差精度。第三节 差分方程所谓差分方程就是把偏微分方程中偏导数用其差商近似来代替，这样将偏微分方程转变 为相应的代数方程组，该方程组称之为 差分方程。 以对流方程式（３－１）为例说明 j+1 差分方程的构造： 由于对流方程式（３－１）在计算 域内任一点均成立，故对点（ｉ，ｊ） 也同样成立，即有： j?u ?tj i+C?u ?xj i=0i-1 i若时间导数用一阶向前差商近似， 空间导数用一阶向后差商近似，即：?u j uij +1 ? uij i ≈ ?t ?t ?u j uij ? uij?1 i ≈ ?x ?x则在点（ｉ，ｊ）的对流方程可以近似写作：图 3―3．ＦＴＢＳ格式差分图25uij +1 ? uij u j ? uij?1 +c i =0 ?t ?x该方程就是方程（３－１）对应的差分方程。 令η =(3 - 22)?tc ，该系数通常称为　Ｃｏｕｒａｎｔ数，则式（３－２２）化简为： ?xuij +1 = (1 ? η )uij + ηuij?1(3 - 23)在实际计算时，还需要有定解条件，问题才能适定，偏微分方程才能有解，设Ｃ≥０且 空间区域为有界域（０≤ｘ≤ｌ） ，则方程的定解条件为：?u ( x,0) = f ( x) ? 定解条件： ?u(0, t) = g(t)相应于边界及初始时刻有如下离散式：初始条件 边界条件(3-24)u i0 = f (i?x) u 0j = g ( j?t )故原始偏微分方程定解问题：(3-25)?u ? ?u ? ?t + C ?x = 0 ? ? ?u ( x, t ) t =0 = f ( x) ? ?u ( x, t ) x =0 = g (t ) ? ?0 ≤ t & ∞,0 ≤ x ≤ l(3-26)通过差商离散，用下列方程组来近似：? u ij +1 ? u ij u ij ? u ij?1 +C =0 ? ?t ?x ? ? 0 ?u i = f (i?x) ? j ?u 0 = g ( j?t ) ? ?(3-27)差分方程与其偏微分方程的定解条件的离散近似一起称为相应偏微分方程定解问题的差 分格式。差分格式（３－２７）称之为ＦＴＢＳ格式。方程组（３－２７）的解 u i 作为原始 定解问题（３－２６）的近似解。 下面以具体例子为例说明如何求解差分方程组(3―27)： 例如：设有定解问题：j?u ? ?u 0 ≤ t & ∞,0 ≤ x ≤ 8 ? ?t + C ?x = 0 ? ? ?u ( x, t ) t =0 = f ( x) = 0 ? ?u ( x, t ) x =0 = g (t ) = 1 ? ?26采用ＦＴＢＳ格式:? u ij +1 ? u ij u ij ? u ij?1 =0 + C ? ? ? t x ? ? 0 (i = 1,2,3, Λ Λ ) ?u i = f (i?x) = 0 ? j (j = 1,2,3,Λ Λ ) ?u 0 = g ( j?t ) = 1 ? ?求解如下： ａ．取Δｘ＝１，Δｔ＝0.5, C=1.0 则 η =?tC = 0.5 ，有： ?xu ij +1 =u ij + u ij?1 2(3-28)结果如下： 表 3―1．Δｘ＝１，Δｔ＝0.5 时，ＦＴＢＳ格式的计算表 j i 0 1 2 3 4 5 0 0 0 0 0 0 0 1 1 0 0 0 0 0 2 1 .5 0 0 0 0 3 1 .75 .25 0 0 0 4 1 .875 .5 .125 0 0 5 1 . . 0 b．取Δｘ＝１，Δｔ＝2，C=1.0 则 η =6 0 0 0 0 0 07 0 0 0 0 0 08 0 0 0 0 0 0?tC = 2 ，有： ?xu ij +1 = ?u ij + 2u ij?1结果如下： 表 3―２．Δｘ＝１，Δｔ＝２时，ＦＴＢＳ格式的计算表 j 0 1 2 3 4 5 0 1 0 0 0 0 0 1 1 0 0 0 0 0 2 1 2 0 0 0 0 3 1 0 4 0 0 0 4 1 2 -4 8 0 0 5 1 0 8 -16 16 06 0 0 0 0 0 07 0 0 0 0 0 0通过上面的例子可见，对同一定解问题的同一差分格式（ＦＴＢＳ）其不同的空间与时 间步长，将得到不同的结果，如果作为原始定解问题的近似解，那一个解精度高呢？。从 表３－２可见随着计算的推移下去，其解将趋近于无穷大，这就是我们下面所要介绍的不稳 定，很显然不稳定的解是不能作为原定解问题的近似解的。 我们知道偏导数的差商近似并非一种，故其同一偏微分方程的差分方程也并非一个，可 以有若干个，故对原始定解问题也相应有若干种差分格式。 方程（３－１）除了ＦＴＢＳ格式外，还可列出下列若干种格式： 1．ＦＴＣＳ格式：27对方程在点（ｉ，ｊ）采用时间向 前差分、空间中心差分得：u ij +1 ? u ij u j ? u ij?1 + C i +1 =0 2?x ?t(3-29) ２．ＦＴＦＳ格式： 对方程在点（ｉ，ｊ）采用时间、 空间向前差分得：j+1ji-1uij +1ii+1u ? ui ? ui +C =0 ?t ?xj j i +1 j图 3―4．ＦＴＣＳ格式差分图(3-30) ３．蛙跳格式： 对方程在点（ｉ，ｊ）采用时间、空间中心差分得：u j ? u ij?1 u ij +1 ? u ij ?1 + C i +1 =0 2?t 2?x(3-31)j+1j+1j j i-1 i i+1 图 3―5．ＦＴFＳ格式差分图 i i+1 j-1图 3―6．蛙跳格式差分图从式（３－２２） 、 （３－２９） 、 （３－３０） 、 （３－３１）可见，格式中若知道第ｊ时 间层上的 u i 值，可由该差分式直接算出第ｊ＋１时间层上的 u ij j j +1值，故称这类格式为显式格式。而一些格式不能直接从ｊ时间层上 u i 值直接解出，需联立求解ｊ＋１层上的值，此 类格式称之为隐式格式。 通过上述分析可见，对同一个定解问题，可以有多种差分格式，多种步长参数来近似， 从而也得到若干个差分近似解。那么这些解是否可以都作为原定解问题的近似解？如可以， 那些解精度高，如不能则又是为什么？这些要通过下面的相容性、稳定性及收敛性分析，才 能得到满意的回答。28第四节 截断误差和相容性以ＦＴＢＳ格式为例来进行分析： 将式（３－８）改成对时间的差商后及式(3-9)分别代入到式(3-22)中有：?u ?tj i+1 ? 2u 2 ?t 2j i? ?t + Ο(?t 2 ) + C?u j C ? 2 u j 2 i ? i ??x + Ο( ?x ) = 0 ?x 2 ?x 2化简：　?u ?tj i?u j 1 ? 2 u +C i + ?x 2 ?t 2j iC ? 2u j ? ?t ? ??x + Ο(?t 2 , ?x 2 ) = 0 2 i 2 ?x(3-32)式（３－３２）是由式（３－２２）通过Ｔａｙｌｏｒ展开得到的，故该两式虽然形式不同， 但它们其实是等价的，可以将式（３－３２）称为式（３－２２）的等价方程。 比较方程（３－３２）与原始偏微分方程（３－１）可见经过差分离散后，比原方程多 出了一部分，这部分称之为格式的截断误差，以Ｒ表示之。R=1 ? 2u 2 ?t 2j i? ?t ?C ? 2u j ??x + Ο(?t 2 , ?x 2 ) 2 i 2 ?x(3-33)同理可得对流方程的差分方程的其它格式的截断误差为： ＦＴＣＳ格式的截断误差为：1 ? 2u R= 2 ?t 2j iC ? 3u j ? ?t + ??x 2 + Ο(?t 2 , ?x 3 ) 3 i 6 ?x(3-34)ＦＴＦＳ格式的截断误差为：R=1 ? 2u 2 ?t 2j i? ?t +C ? 2u j ??x + Ο(?t 2 , ?x 2 ) 2 i 2 ?x(3-35)蛙跳格式的截断误差为：R=1 ? 3u 6 ?t 3j i? ?t 2 +C ? 3u j 2 3 3 i ??x + Ο( ?t , ?x ) 6 ?x 3(3-36)一般地讲，对流方程(3―1)的差分方程等价形式可写成:　?u ?u +C +R=0 ?x ?t(3-37)如果当Δｘ、Δｔ-&０时，差分方程的截断误差的某种范数‖Ｒ‖也趋近于零即：?x , ?t →0lim R = 0 ，则表明从截断误差的角度来看，此差分方程是能用来逼近微分方程，通常称这样的差分方程和相应的微分方程相容。如果截断误差的范数不趋于零，则称为不相容，这 样的差分方程不能用来逼近微分方程。 以上只考虑了方程，但从整个问题来看，还应考虑定解条件。若定解问题的定解条件为：B (u ) = g其中Ｂ是微分算子，ｇ是已知函数，而对应的差分问题的定解条件为：29B? (u ) = g其中 B? 是差分算子，则截断误差为：r = B? (u ) ? B(u )只有方程相容，定解条件也相容，即：?x , ?t → 0lim R → 0, lim r → 0?x , ?t →0整个问题才相容。Δｘ-&０、Δｔ-&０的情况有两种：一是各自独立地趋于零，这是无条件 相容；另一是Δｘ、Δｔ之间存在某种关系（譬如要求?x = K ）下同时趋于零，这种情况 ?t下的相容为条件相容。 从上述定义与分析知道，ＦＴＣＳ、ＦＴＦＳ、ＦＴＢＳ及蛙跳格式均为无条件相容。 Ｐ Ｑ 取Ｒ０＝　ｍａｘ（Ｒ，ｒ）＝Ｏ（Δｔ ，Δｘ ） ，称差分格式对时间有P阶精度，对空间 有Q阶精度。B B P P P P第五节 收敛性所谓相容性,是指当自变量的步长趋于零时,差分格式与微分问题的截误差的范数是否趋 于零,从而可看出是否能用此差分格式来逼近微分问题。然而方程(无论是微分方程或是差分 方程)是物理问题的数学表达形式,其目的是为了借助数学的手段来求问题的解。 因此,除了必 需要求差分格式能逼近微分方程和定解条件(表明这两种数学方法在形式上是一致的)外,还 进一步要求差分格式的解(精确解)与微分方程定解问题的解(精确解)是一致的(表明这两种 数学表达方法的最终结果是一致的)。即当步长趋于零时,要求差分格式的解趋于微分方程定 解问题的解。我们称这种是否趋于微分方程定解问题的解的情况为差分格式的收敛性。更明 确地说,对差分网格上的任意结点(i,j),也是微分问题定解区域上的一固定点,设差分格式在 此点的解为 u i ,相应的微分问题的解为 U ( xi , t j ) ,二者之差为：jε i j = u ij ? U ( xi , t j )称为离散化误差。 如果当Δｘ-&０、 Δｔ-&０时,离散化误差的某种范数‖ｅ‖趋近于零,即：lim ε?x , ?t=0则说明此差分格式是收敛的,即此差分格式的解收敛于相应微分问题的解,否则不收敛。与相 容性类似，收敛又分为有条件收敛和无条件收敛。 粗看起来,似乎只要差分格式逼近微分问题（Δｘ、Δｔ→０时，‖Ｒ‖→０,‖ｒ‖→ ０）,其解就应该一致;也就是说，似乎相容性能保证收敛性。其实并不一定如此。这是因为 在分析截断误差时， 是以差分格式与微分问题有同一个解为基础(或以定解域内某足够光滑的 函数为基础),并对此函数分别在点（ｘｉ，ｔｊ）的邻域作Ｔａｙｌｏｒ展开的,其中所有的B B B Bu,?u ?u , 等都是指同一个函数及其各阶导数。所以最后得的截断误差Ｒ、ｒ实质上是当差 ?t ?x分问题与微分问题有同一解时两种方程、两种定解条件之间的误差。Ｒ、ｒ并不能直接表示 两种方程、两种定解条件之间的误差,因此,相容性不能保证收敛性,不能保证二者解的一致。30若没有相容性就更不能得到二者解的一致,故相容性是收敛性的必要条件,有人称相容性是形 式上的逼近。 相容性不一定能保证收敛性。 那么对于一定的差分格式,其解能否收敛到相应的微分问题 的解?这里先举一个例子。现有微分方程定解问题：? du ? + Au = 0 ? dx ? ?u(0) = 10&x&a其解析解为：ｕ（ｘ）＝ｅ 将［０，ａ］等分为ｎ段,则步长Δｘ＝ａ／ｎ。构造差分格式如下:P P－Ａｘ? u i +1 ? u i + Au i = 0 ? ? ?x ? ?u 0 = 1或写成：?u i +1 = (1 ? A ? ?x)u i ? ?u 0 = 1由此可得：i = 0,1,2,3, Λ Λ n - 1?u 1 = (1 ? A ? ?x)u 0 = (1 ? A ? ?x) ? 2 ?u 2 = (1 ? A ? ?x)u1 = (1 ? A ? ?x) ? ?Λ Λ Λ ?u = (1 ? A ? ?x)u = (1 ? A ? ?x) i i ?1 ? i任意结点ｉ的坐标为： xi = i ? 故： i =a nn xi aB B将差分解写成ｘｉ的函数形式：a xi u ( xi ) = (1 ? A ? ?x) = (1 ? A ? ) a n A = {(1 ? ) m } xi min其中 m =n 。当Δｘ→０时,有ｎ→∞，ｍ→∞。而 a A lim(1 ? ) m = e ? A n→∞ m因此得：u ( x i ) = e ? Axi去掉下标ｉ,即为：u ( x) = e ? Ax31这与微分问题的解析解是完全一致的。这个例子说明了差分格式的解收敛于微分问题的 解是可能的。 至于某给定格式是否收敛,则要按具体问题予以证明。 下面再举一个前面曾谈到 过的差分格式,试讨论其收敛性。 微分问题：?u ? ?u =0 ? +C ?x ? ?t ? ?u ( x,0) = g ( x)的ＦＴＢＳ格式为：? u ij +1 ? u ij u j ? u ij?1 +C i =0 ? ?x ? ?t ?u 0 = g ( x ) i ? iB B B B（３－３８）在某结点（ｘｉ，ｔｊ）,微分问题的解为Ｕ（ｘｉ，ｔｊ）,差分格式的解为 u i ,则离散化B B B Bj误差为：ε i j = u ij ? U ( xi , t j )按照截断误差的分析知道：U ( xi , t j +1 ) ? U ( xi , t j ) ?t+CU ( xi , t j ) ? U ( xi , t j ?1 ) ?x= Ο(?x, ?t )以ＦＴＢＳ格式（３－３８）中的第一个方程减去上式得：ε i j +1 ? ε i j?t或写成： ε ij +1+Cε i j ? ε i j?1?x= Ο(?x, ?t )= (1 ? η )ε i j + ηε i j?1 + Ο(?x, ?t ) ? ?t若条件Ｃ≥０和η≤１成立,即０≤η≤１,则：ε i j +1 ≤ (1 ? η ) ε i j + η ε i j + Ο(?x, ?t ) ? ?t≤ (1 ? η ) max ε i j + η max ε i j + Ο(?x, ?t ) ? ?ti i式中 max ε i 表示在第ｊ层所有结点上 ε i 的最大值。jji由上式知,对一切ｉ有：ε i j +1 ≤ max ε i j + ?t ? Ο(?x, ?t )i故有： max ε iij +1≤ max ε i j + ?t ? Ο(?x, ?t )i于是：32max ε i1 ≤ max ε i0 + ?t ? Ο(?x, ?t )i i i imax ε2 i 3 i≤ max ε i0 + 2?t ? Ο(?x, ?t )i imax ε ≤ max ε i0 + 3?t ? Ο(?x, ?t )ΛΛΛΛmax ε i j ≤ max ε i0 + j?t ? Ο(?x, ?t )i i综合得： max ε i ≤ max ε i + j?t ? Ο( ?x, ?t )j0ii由于初始条件给定了函数ｕ的初值,初始离散化误差 ε i = 0 ,并且ｊΔｔ＝ｔｊ是一有限0B B量,因而：max ε i j ≤ Ο(?x, ?t )i可见本问题ＦＴＢＳ格式的离散化误差与截断误差具有相同的量级。最后得到：?x , ?t → 0lim (max ε i j ) = 0i这样说明了,当０≤η≤１时,本问题的ＦＴＢＳ格式收敛。这种离散化误差的最大绝对 值趋于零的收敛性情况称为一致性收敛。 此例介绍了一种证明差分格式收敛的方法,同时表明了相容性和收敛性的关系:相容性是 收敛性的必要条件,但不一定是充分条件,如本例就要求０≤η≤１。第六节 稳定性首先介绍一下差分格式的依赖区间、决定区域和影响区域。还是以初值问题?u ? ?u =0 ? +C ?x ? ?t ? ?u ( x,0) = f ( x)为例。先看ＦＴＣＳ格式,如图３－７a，欲计算第二层ｐ点的函数值,必先知道第一层上ａ、 ｂ、ｃ这三点的函数值,故说ｐ点的解依赖于ａ、ｂ、ｃ这３点的解。而ａ点的解又依赖于第P a b cPPAdefBA 图 3-7 依赖区间示意图BBA０层(初值线)上Ａ、ｄ、ｅ的初值,ｂ点的解依赖于ｄ、ｅ、ｆ的初值,ｃ点的解依赖于ｅ、 ｆ、 Ｂ的初值。 因此ｐ点的解依赖于初值线ＡＢ段上所有结点的初值,故称ＡＢ段上所有结点 为点ｐ的依赖区间。又 pAB 三角形区域内任一点的依赖区间都包含在 AB 之内,即该区域内任33一结点上的解都由 AB 段上某些结点的初值所决定,而与 AB 以外结点的初值无关,故称此三角 形区域为 AB 区间所决定的区域。 这里为方便起见,是以第二层的 p 点为例的,事实上对任意层 的任一结点,都在初始层上有一对应的依赖区间,而初始层的任一区间都有一对应的决定区 域。 FTFS 格式和 FTBS 格式的依赖区间分别为图 3-7b 和图 3-7c 中的 AB 线段的全部结点;图中 阴影部分为 AB 所决定的区域。随着时间的推移,一点函数值将影响以后某些结点的解。如图 ３－８,设 p 为第 j 层的某结点,当用 FTCS 格式计算第 j+1 层上的结点值时,a、b、c 这 3 点 的解必须用到 p 点的函数值,在第 j+2 层上则有更多点的解受 p 点函数值的影响。所有受 p 点函数值影响的结点总和为 p 点的影响区域,如图 3-8 中阴影所示区域。P ＦＴＣＳ格式 FTFS 格式 图 3-8 影响区域示意图PP FTBS 格式由上述可知,同一微分问题,当采用不同差分格式时,其依赖区间、 决定区域和影响区域可 以是不一致的。依赖区间、决定区域和影响区域是由差分格式本身的构造所决定的并与步长 比?t 有关。 ?x例如微分问题?u ? ?u =0 ? +C ?x ? ?t ? ?u ( x,0) = 0其解为零,即ｕ （ｘ， t） ＝０。 若用 FTBS 格式计算,且计算中不产生任何误差,则结果也是零, 即：?u ij = 0 ? ?j = 0,1,2,3, Λ i = 0,±1,±2, Λj假设在第 j 层上的第 i 点,由于计算误差，得到 u i = ε 。不妨设 i＝０,ｊ＝０，ε＝１， 即相应于 FTBS 格式改写成：?u ij +1 = (1 ? η )u ij + ηu ij?1 ? 0 ?u 0 = 1 ? 0 ?u i = 0, i ≠ 0 ?t = 0.5,1,2 ,列表计算如下: 现分别取 ?x34表 3-3.?t = 0.5 时的计算结果 ?xi -4 0 0 0 0 0 -3 0 0 0 0 0 -2 0 0 0 0 0 -1 0 0 0 0 0 0 1 .5 .25 .125 . .5 .5 .375 .25 2 0 0 .25 .375 .375 3 0 0 0 .125 .25 4 0 0 0 0 .0625j 0 1 2 3 4表 3-4.?t = 1 时的计算结果 ?xi -4 0 0 0 0 0 -3 0 0 0 0 0 -2 0 0 0 0 0 -1 0 0 0 0 0 0 1 0 0 0 0 1 0 1 0 0 0 2 0 0 1 0 0 3 0 0 0 1 0 4 0 0 0 0 1j 0 1 2 3 4表 3-5.?t = 2 时的计算结果 ?xi -4 0 0 0 0 0 -3 0 0 0 0 0 -2 0 0 0 0 0 -1 0 0 0 0 0 0 1 -1 1 -1 1 1 0 2 -4 6 -8 2 0 0 4 -12 24 3 0 0 0 8 -32 4 0 0 0 0 16j 0 1 2 3 4这个例子一方面显示了该格式的影响区域,另一方面显示了当 的影响在数值上有很大的不同,当?t & 1.0 时产生的影响在数值上将越来越大。数值上的差别引起了质的不同,因而出现了稳 ?x?t ≤ 1.0 时 , 所 产 生 的 影 响 在 数 值 上 不 再 扩 大 ; 当 ?x?t 值不同时计算误差所产生 ?x定性问题。差分格式的数值稳定性,早在 1928 年就由 R.Courant、K.O.Friedrchs 和 H.Lewy 等人发现,并提出了关于双曲型方程差格式稳定性的必要条件(简称 CFL 条件)。 此后在这方面 作了不少研究工作。 1950 年公开发表了 von Neumann 提出的稳定性分析法,这是现在比较广 泛地用来确定稳定性准则的一种分析方法。 在有限差分法的具体运算中,计算误差总是不可避免的,如舍入误差,以及这种误差的传 播、 积累。 然而人们通过大量的实践和理论分析发现,同一问题的各种差分格式在某一定条件 下,若计算中某处产生了误差,则这个误差将对以后的计算产生影响。如果这一误差对以后的35影响越来越小,或是这个影响保持在某个限度以内,像上面例子中?t ≤ 1.0 的情况,那么就称 ?x这个差分格式在给定的条件下稳定,这个条件就是它的稳定准则。 如果误差的影响随着的ｊ增 大越来越大,像上面?t & 1.0 的情况,使计算的结果随着ｊ的增大越来越偏离差分格式的精 ?x确解,而毫无实用价值,那么这种情况就是不稳定的。 为了说明稳定性,现再以对流方程（３－１）的 FTCS 格式:u ij +1 ? u ij u ij+1 ? u ij?1 +C =0 ?t 2?xB B(3-39)为了分析时简单明了,不妨假设,在ｔ ｊ 时刻以前的运算中也不产生任何新的计算误差,只在 ｊ ｊ ｔｊ时刻产生了计算误差ε , 试看这个误差对以后的影响。由于ε 的加入, ｔｊ时刻ｉ及 ｉ±１结点函数值成为：B B P P P P B B(u ij ) / = u ij + ε i j (u ij±1 ) / = u ij±1 + ε i j±1由于 ε i 的影响,这时 u ij j +1也不是原先的数值,而成为 u i/ j +1。不妨也将其写成：u i/ j +1 = u ij +1 + η ij +1这里的 u i/ j +1不是新产生的计算误差,而是 ε ij +1的传播。以所有带“′”的量按 FTＣS 格式写成差分方程,即:(u ij +1 + ε i j +1 ) ? (u ij + ε i j ) (u j + ε i j+1 ) ? (u ij?1 + ε i j?1 ) + C i +1 =0 ?t 2?x以此式减去式（３－３９）得：ε i j +1 ? ε i j?t+Cj j ε i+ 1 ? ε i ?12 ?x=0(3-40)这就是误差的传播方程。与式（３－３９）比较可见,线性齐方程的误差传播方程与原方程形 式相同,只需将函数ｕ换成误差ε。 上式可改写成：η ε i j +1 ? ε i j = ? (ε i j+1 ? ε i j?1 )2因为时刻ｔｊ产生了计算误差ε ,传到ｔｊ＋Δｔ时刻变成εB B P P B B Pｊｊ＋１P,所以 ε ij +1? ε i j 是误差的增长,若以 ?ε i 表示有：j?ε i j = ? (ε i j+1 ? ε i j?1 ) 2B Bη(3-41)上式右端表示误差的增长。下面讨论其增长的情况。 为了形象化,我们假设ｔｊ时刻产生的误差是沿结点振荡的,其振幅是以结点增加的,如图36３－９（ａ） 。 1.设Ｃ≥０，当 ε i & 0 时,必有:jεPjPε i j+1 & 0, ε i j?1 & 0 以及ε ij+1 & ε i j?1 η ?ε i j = ? (ε i j+1 ? ε i j?1 ) & 02当 ε i & 0 时,必有：j图 3―９（ａ）εj i +1& 0, εj i ?1& 0 以及εj i +1&εj i ?1εPj+1P?ε i j = ? (ε i j+1 ? ε i j?1 ) & 0 2可见在Ｃ＞０的情况下, ?ε ij jη图 3―９（ｂ）与 ε i 符号一致,结果使误差振幅εPj+1P随ｊ单调地增长,如 图３－９ （ｂ） 所示。 此种格式是不稳定的。 图 3―９（c） 这种单增长型的不稳定称为静力 不稳定性。除非改变格式,否则这种不稳定性是无法消除的。 ２． 在Ｃ＜０的情况下,若 ε i & 0 ,则 ?ε i & 0 ;若 ε i & 0 ,则 ?ε i & 0 。 这时 ε i 与?ε i 的j j j j j j符号相反,若不是过大,则 ?ε i 对 ε i 将是一个校正,使误差随着ｊ的增大逐步减小,这便为稳j j定状态。然而如果η过大，校正就会过头,称为过冲,使新的 ε ij +1比原来的 ε i 幅度更大,如图j３－９（ｃ）所示,形成不稳定状态。这种过冲型振荡的不稳定称为动力不稳定。显然,这种 不稳定可以用限制η的办法消除,对应的稳定状态称为条件稳定,对η的限制称为稳定条件或 稳定准则。 以上讨论对误差 ε i 作了特殊的假定,这是为了使稳定性形象化,以便于理解。实际误差是j随机的,不可能符合上述假。因而以上讨论只能给出一个定性的概念。由讨论知:有些格式在 一定条件下稳定,称为条件稳定;有些格式在任何条件下都不稳定。 以后还会看到,有些格式是 无条件稳定的,称为完全稳定。 值得强调的是,这里所说的某种格式稳定,是一种简略的说法。 实际上,不可能孤立地研究 某种格式,必然是针对某一微分问题来研究某差分格式,所以差分格式的相容性、收敛性、稳 定性是针对给定的微分问题而言的。 现在以适当的数学式子给出稳定性定义。为此将差分解 u i 表示为连续函数Ｚ（ｘ，ｔ）, 则稳定性的一种定义为： ‖Ｚ（ｘ，ｔ）‖≤Ｋ‖Ｚ（ｘ，０）‖ 这里 K 是某个有限常数，称为 Lipschitz 常数,不随Δｘ、Δｔ→０而变。这就是说,当上述 不等式成立时,只要差分问题初始值所含的误差为小量时,此后的解与差分问题的精确解的误j37差也一定为小量。由于计算误差不仅可以来自初值(包括在某一时刻前的任一时刻)还可以来 自边界值,而且可以来自右端项,所以也有将稳定性定义为: ‖Ｚ‖≤Ｋ１‖ＤΔ（Ζ）‖＋Ｋ２‖ＢΔ（Ｚ）‖ 其中ＤΔ和ＢΔ分别是对应于微分方程和定解条件的差分算子,Ｋ１、Ｋ２分别是对应于 ＤΔ和B B B B B B B B B B B B B B B B B BB? 的Ｌｉｐｓｃｈｉｔｚ常数。若取：Ｋ＝ｍａｘ（Ｋ１，Ｋ２） 则为： ‖Ｚ‖≤Ｋ(‖ＤΔＺ（Ζ）‖＋‖ＢΔ（Ｚ）‖) 在建立了稳定性概念以后,进一步是如何判定格式是稳定,或在什么条件下稳定。 最易联想到的是象前证明收敛性那样证明误差的影响有界。在收敛性中所述的误差是指 离散化误差,而稳定性中是计算误差。计算误差在某一计算步骤中都有可能产生,通常是设第B B B B B B B Bｊ层有计算误差 ei ,在此以前没有计算误差,在此以后也无新的计算误差,只有 ei 的影响(传 播)。现以微分问题：jj?u ? ?u =0 ? +C ?x ? ?t ? ?u ( x,0) = f ( x)的 FTBS 格式（见式３-23） ：u ij +1 = (1 ? η )u ij + ηu ij?1为例。不妨假设在初层（ｊ＝０）的初值中有舍入误差 ε i ,则初值成为 u i + ε i 。容易知道,0 0 0上述 FTBS 格式的误差传播方程为：ε i j +1 = (1 ? η )ε i j + ηε i j?1当Ｃ＞０和η≤１时,由上式得(3-42)ε i j +1 ≤ (1 ? η ) ε i j + η ε i j?1 ≤ max ε i ji故对所有ｊ：max ε i j +1 ≤ max ε i ji i并可推得：max ε i j +1 ≤ max ε i0i i即由于 ε i 的影响而产生的第ｊ层的误差不会大于ε ,这说明ε 所产生的影响是有界的。或0０P P０PP者说,当Ｃ＞０时和η≤１时,符合前面所述的第一种稳定性意义,故上述微分问题的FTBS格 式是稳定的,并且由于是一层上的最大误差有界,必定所有误差一致有界,故称此为一致稳定。 若对上述微问题采用 FTCS 格式,其误差传播方程为：38η ε i j +1 = ε i j ? (ε i j+1 ? ε i j?1 ) η ε i j +1 ≤ max ε i j + (max ε i j + max ε i j )max ε ii j +122 i ≤ (1 + η ) max ε i jii由此推得：max ε i j +1 ≤ (1 + η ) j max ε i0i i由于η为有限量,当ｊ→∞时上面等式成为：max ε i j +1 ≤ ∞i这是没有意义的不等式。 可见这种证明稳定性的方法,有时得不到任何结果。 下面介绍一种由 von Neumann 提出的证明稳定性的方法。 仍用上述问题的 FTBS 格式,其误差传播方程为式 （３ －４２） 。 ０ 假设 ε 是各结点上的离散量。现将其扩充成空间连续的量,并设初始误差 ε （ｘ）P P在［０，π］内平方可积（若a≤ｘ≤ｂ,则可用 ξ = 式（３－４２）表明 ε iP P2( x ? a) 代替ｘ，０≤ξ≤２π） 。 b?aｊP Pj +1与 ε i 成线性关系,故可由ε 在［０，２π］内平方可积推得ε 、j０P P P P P Pεｊ＋１在［０，２π］内平方可积。因此εｊ、εｊ＋１可展开成傅氏级数。∞ ? j = ε ( ) x C kj e ikx ∑ ? ? k = ?∞ ? ∞ ?ε j +1 ( x) = C kj +1e ikx ∑ ? k = ?∞ ?式中 i =∞? 1 是虚数单位。将εｊ、εｊ＋１代入式（３－４２）并整理后得：P P P Pk = ?∞∑ {Cj +1 k? [(1 ? η )C kj + ηC kj e ?ik?x ] e ikx = 0}此式相当于将零展开成傅氏级数,｛ ｝内相当于傅氏系数,对所有的ｋ它都应等于零:C kj +1 ? [(1 ? η )C kj + ηC kj e ? ik?x ] = 0由此得,对一切ｋ有：C kj +1 = [(1 ? η )C kj + ηC kj e ? ik?x ] = GC kj其中Ｇ为放大因子,G (k , ?x, ?t ) = (1 ? η ) + ηe ? ik?x当时ｊ＝０,有：1 Ck = GC k039C k0 是初始误差ｅ０（ｘ）的傅氏系数。当j＝１时,有：P P1 C k2 = GC k = G 2 C k0递推得：C kj = G j C k0按傅氏分析的 Parseval 定理1 2π取∫2π0ε j ( x ) δx =2k = ?∞∑CP P∞j 2 k1 2π∫2π0ε j ( x) dx 为εｊ（ｘ）的范数‖εｊ（ｘ）‖,则:P Pε j ( x) =21 2πｊP P∫2π0ε j ( x ) δx =2k = ?∞∑C∞j 2 k=k = ?∞∑G∞2jC k02此时,如果｜Ｇ｜ 对任意ｊ和一切ｋ一致有界: ｊ ｜Ｇ｜ ≤Ｍ 而Ｍ为一正有限值,则：P Pε ( x) ≤ Mj222k = ?∞∑C∞0 k 22ε j ( x) ≤ M 2 e 0 ( x) ε j ( x) ≤ M e 0 ( x)2最后一式表明,当初始误差的范数充分小时,第ｊ层上的误差范数也充分小,符合稳定性 定义,因而可认为是稳定的。 由于这里定义的范数相当于均方根值,故称此种稳定为平均稳定。 这里稳定的必要条件是G ≤M或：jG ≤M1 jｊP PＭ是一与ｊ无关的有限正值。若取Ｍ＝（１＋ε） ，则： ｜Ｇ｜≤1+ε ε为某正数,注意到 j 是正整数，而且可以是一个很大的正整数，因此，欲使 M 为有限正数, ε必须是一个小的正值,故有：M = (1 + ε ) j = 1 + jε +j ( j ? 1) + Λ Λ ≈ 1 + jε 2! tj ?t,故须使ε等于Δｔ与某一与ｊ无关的由于Ｍ与ｊ无关,同时考虑到ε为小的正值和 j = 有限正值ｍ之乘积 ε＝ｍΔｔ40从而得到稳定的条件： ｜Ｇ｜≤１＋ｍΔｔ 常称为 von Neumann 稳定性条件。 由于此条件中的ｍ不确定,所以实际使用的稳定性条件为： ｜Ｇ｜≤１ 以上方法概括起来是,由傅氏展开通过误差传播方程得到放大因子Ｇ,如果能证明对于任 意ｊ和一切ｋ有： ｜Ｇ｜≤１ 则此格式具有平均稳定性。 实际运用此方法时,若不考虑原来的数学意义,只求能得到放大因子Ｇ,则不需取整个傅 氏级数,而任取一个分量,即令：ε i j = C j e ikxiε i j +1 = C j +1e ikxii(3-43)ε i j±1 = C j e ik ( x ± ?x )将这些代入式（３－４２）得：C kj +1 = [(1 ? η )C kj + ηC kj e ? ik?x ]其中Ｇ就是前面得到的放大因子 －ｉｋΔｘ Ｇ＝（１－η）＋ηｅ 对于不同格式,一般地说Ｇ是不一样的。接下来的问题是要判定是否有，或在什么条件下有： ｜Ｇ｜≤１ －ｉｋΔｘ 由于ｅ ＝ｃｏｓ(ｋΔｘ)－ｉｓｉｎ(ｋΔｘ),故P P P PG = 1 ? 4η (1 ? η ) sin 2 (显然,要使｜Ｇ｜≤１,必须2k?x ) 21 ? 4η (1 ? η ) sin 2 (即:k?x ) ≤1 2 k?x 4η (1 ? η ) sin 2 ( )≥0 2２P P当Ｃ＞０时,由于４ηｓｉｎ (ｋΔｘ／２)≥０,故仅要求: １－η≥０ 即： η≤１ 因此当０≤η≤１时上述问题的 FTBS 格式稳定。 当Ｃ＜０时,除了ｓｉｎ(ｋΔｘ/2)＝０时外,｜Ｇ｜都大于 1,故不稳定。 同一个问题,用不同的方法所得的稳定性条件一致。这里还值得一提的是,该问题的稳定 性条件０≤η≤１,与收敛性条件也是相同的,这表明稳定性与收敛性之间可能有一定联系。 下面再举一个例子,即用傅氏方法分析对流方程的 FTCS 差分格式(3-29)的稳定性。该格 式的误差方程可写成:η ε i j +1 = ε i j ? (ε i j+1 ? ε i j?1 )2将式(3-43)代入到上式并整理得：(3-44)41C j +1 = [1 ?η2(e ik?x ? e ?ik?x )]C j= GC j即放大因子：G = 1?η2 = 1 ? iη sin(k?x)(e ik?x ? e ?ik?x )故:G = 1 + η 2 sin 2 k?x & 1因而ＦＴＣＳ格式为一不稳定格式。 类似地可以得到ＦＴＦＳ格式的稳定性条件为：－１≤η≤０，蛙跳格式的稳定性条件 为：－１≤η≤１。 Von Neumann 提出的方法是目前比较广泛采用的一种稳定性分析法,但它主要用于线性 初值问题的稳定性分析。对于非线性问题用局部线性化的方法加以推广。局部线性化方法假 定非线性系数变化得很缓慢,因而可用局部网格结点上的函数值代入后作为常数处理,并认为 每一网格结点上的计算稳定性与相邻结点无关,这样可对每一网格结点作傅氏分析,找出关于 Δｔ的局部稳定极限值,所以网格结点上最小的局部稳定极限值就可作为整个差分问题的稳 定极限值。 对初边值问题,局部线性化稳定性分析可以避免由于边界处理不当而引起的计算稳 定。第七节Lax 等价定理前面讨论了差分问题的相容性、 收敛性和稳定性。 已经知道,相容性是收敛性的必要条件; 还发现,稳定性与收敛性有一定的联系。Lax 等价定理就是阐述相容性、收敛性和稳定性三者 之间的关系的。 Lax 等价定理:对一个适定的线性微分问题及一个与其相容的差分格式,如果该格式稳定 则必收敛,不稳定必不收敛。换言之,若线性微分问题适定,差分格式相容,则稳定性是收敛性 的必要和充分的条件。这也可表示为： 稳定性 ?===============? 收敛性差分格式相容 线性微分问题适定下面对此定理作一些简略的说明。由于在定解域内有 Ｄ（ξ）＝ｆ 及 ＤΔ（Ｚ）＝ｆ 其中Ｄ和Ｄ Δ 分别为微分算子和差分算子,是线性的;ｆ是已知函数;ξ和Ｚ分别为微分解和 差分解。两式相减得： ＤΔ（Ｚ）－Ｄ（ξ）＝０ 改写成 ［ＤΔ（Ｚ）－ＤΔ（ξ） ］＋［ＤΔ（ξ）－Ｄ（ξ） ］＝０ 因为算子是线性的,故式中第一个［］内相当于ＤΔ（Ｚ－ξ）;而第二个［］内就是截断误 差Ｒ,所以有 ＤΔ（Ｚ－ξ）＝－Ｒ (３－４５） 若定解条件为 Ｂ（ξ）＝ｇ 及 ＢΔ（Ｚ）＝ｇB B B B B B B B B B B B B B B B B B42其中Ｂ和ＢΔ分别为微分算子和差分算子,且是线性的;ｇ为已知函数。按照以上对方程的同 样推导法,可导得： ＢΔ（Ｚ－ξ）＝－ｒ （３－４６） 其中ｒ＝ＢΔ（ξ）－Ｂ（ξ）是截断误差。若差分格式是稳定的,按稳定性的定义,应该有：B B B B B BZ ? ξ ≤ K [ D? ( Z ? ξ ) + B? ( Z ? ξ ) ]将式（３－４５） 、 （３－４６）代入得： ‖Ｚ－ξ‖≤Ｋ［‖Ｒ‖＋‖ｒ‖］ 当差分格式相容时,可得：?x , ?t →0lim Z ? ξ = 0从而保证了收敛性。 根据此定理,在线性适定和格式相容的条件下,只要证明了格式是稳定的,则一定收敛;若 不稳定,则不收敛。由于收敛性的证明往往比稳定性更难,故人们就可以把注意力集中在稳定 性的研究上。 从上节讨论知道,对流问题的 FTCS 格式相容但不稳定,按 Lax 等价定理知道它也不收敛。 正如前面已经说过的,相容只是形式上的逼近,它不能保证格式的收敛性,所以在构造差分格 式时要特别注意,一要相容,二要稳定(对线性适定问题此两者保证了收敛性),同时还要求有 一定的精度。第八节差分方程数值效应微分方程是描述物理量在时间和空间上的连续变化的规律,若将本来是连续变化的物理 量人为地在时间和空间上加以离散化,用差分方程来描述离散化后物理量的变化规律,必然会 引起一些离散误差,使得差分方程的定解问题不能很好地逼近原来的物理问题。 这种差分方程 使原系统的物理性质和规律遭到某种程度的歪曲和破坏的作用称为数值效应或离散近似的伪 物理效应。因此，从计算的准确度方面来考虑，应最少地歪曲本来的物理现象。从计算的稳 定性看，紊乱现象，包括不稳定现象也正是由于原来系统的物理性质和物理规律遭到破坏的 结果。要想解决实际计算中的准确度问题，除了一般的数学逼近考虑之外，还必须对这些效 应有明确的概念，以便从物理上来考虑数值格式的合理性，减少数值效应的影响。尤其对于 非线性问题，数学理论还很不完善，更有必要对此进行分析。一、 “逆风”效应仍以对流方程（３－１）为例，该方程所描述的物理现象是某种不扩散的物理量ｕ以波 速Ｃ传播。例如在河流中某断面处注入浓度为Ｑ的不会扩散的污染物质，水流速度为Ｃ，那 么对流方程就是描述该污染物质被水流带走过程中浓度随时间和空间的变化规律。从物理概 念上来说，该污染物质只可能被带到下游方向，而不可能被输送到上游，这从上章介绍该方 程的定解问题时可以看到。 但是若用ＦＴＣＳ格式：u ij+1 ? u ij?1 u ij +1 ? u ij +C =0 2?x ?t43该格式前面已介绍，是一不稳定的差分格式。通过下面的分析可以看出其严重违背了客观物 理规律。上式化简成：u ij +1 = u ij ?η2(u ij+1 ? u ij?1 )(3-47)假定在某瞬时ｊ在某一断面ｋ处引入某一物理量ｕ，按（３－４7）计算结果如下表所示 (设η＝１，ｕ＝１） 。 k k-3 k-2 k-1 k k+1 k+2 k+3 j j 0 0 0 1 0 0 0 j+1 0 0 -.5 1 .5 0 0 j+2 0 .25 -1 .5 1 .25 0 由表可见，物理量向上、向下游两个方向传播，出现了与波速相反方向传播的不合理现象， 称之为“逆风”效应。 因此差分方程（３－29）所描述的物理量的运动规律与它所近似的原问题固有的规律相 差甚大，故不仅计算结果误差很大，而且也往往是引起差分格式不稳定的一个因素，前面已 证明了该差分格式是无条件不稳定的。 比较ＦＴＣＳ、ＦＴＦＳ、ＦＴＢＳ格式，差别在于 FTCS 格式中把原来对空间的向后或 向前差商作为一阶空间偏导数近似改用较高精度的中心差商,但结果适得其反。 所以在构造差 分格式时不能盲目追求精度阶数较高的差分格式，应该从物理现象上去分析采用何种格式较 为合理。例如对流方程(3―1)采用一阶精度 FTBS 格式或 FTFS 格式更能近似原问题的物理现 象。由前面的稳定性分析知，是采用何种格式还与波速的方向紧密联系起来，当Ｃ为正时采 用ＦＴＢＳ格式，当Ｃ为负时采用ＦＴＦＳ格式。若Ｃ的符号在计算过程中会改变，例如潮 水河道，则可以采用下面通用的“逆风”格式：u ij +1 ? u ij C i j + [(1 ? α )(u ij+1 ? u ij ) + (1 + α )(u ij ? u ij?1 )] = 0 2?x ?tj(3-48)式中：α＝ｓｉｎｇ（ C i ） ，当Ｃ＞０时，α＝１，当Ｃ＜０时，α＝－１。采用上述格式 就保证了计算过程中下游的状态不影响上游。 “逆风”格式的稳定条件为：｜η｜≤１二、物理耗散与弥散首先讨论一维行波的性质，波高可写为：ξ ( x, t ) = Ae ik ( x ?ct )(3-49)在一固定时刻，空间上相差一个波长λ其波高相等，即：Ae ik ( x + λ ?ct ) = Ae i[ k ( x ?ct ) + 2π ]所以ｋλ＝２π， k =2πλ。ｋ是２π距离中所包含波的个数，称为波数。在一固定位置，时间相差一个周期Ｔ，其波高相等，即：Ae ik ( x ?c ( t +T )) = Ae i[ k ( x ?ct ) ? 2π ]44所以ｋＣＴ＝２π， C =2π 1 ? = λf ，ｆ为频率。Ｃ表示单位时间传播的距离，称为相速 k T度。定义相函数φ＝ｋ(ｘ－Ｃｔ）,沿着φ＝ｃｏｎｓｔ直线上，波高处处相等，这样可得 到： k =1 ?Φ ?Φ ，波数是相函数的空间导数。 C = ? ，相速度是相函数的时间导数与空 ?x k ?t间导数之比。 物理耗散现象是指波幅Ａ因阻尼作用而衰减的现象，而弥散现象是指波的相速度Ｃ随波 数发生变化的现象。例如深水波相速度 C =g ，长波λ大，波数ｋ小，相速度Ｃ就大，传 k gh ，播快。 反之短波的相速度小， 传播慢。 这种波就称为弥散波或称色散波。 而浅水波 C =相速度与波数无关，称为非弥散波。 为了弄清数值耗散和色散效应的本质，下面以三个模型方程为例，讨论物理耗散与弥散 现象。 １、无耗散和色散的模型 ｕｔ＋Ｃｕｘ＝０ （３－５０） 设方程的解可以用Ｆｏｕｒｉｅｒ级数表示B B B Bu = ∑ u k = ∑ Ak e ik ( x ?ct )k k则波数为ｋ的波分量可写为u k = Ak e ik ( x ?ct )这一波分量是满足方程（３－５０）的。考虑经过Δｔ时间后，其解可写成：u k ( x, t + ?t ) = Ak e ik [ x ?c ( t + ?t )]定义放大因子：r=u k ( x, t + ?t ) = e ik?x = r e i?Φ u k ( x, t )(3-51)放大因子的模｜ｒ｜＝１，说明波在传播过程中，经过Δｔ时刻后振幅没有衰减。放大因子 的幅角Δφ＝－ＣｋΔｔ,表明传播Δｔ时刻后， 同一位置ｘ处波的相位迟后为Δφ， 它是相 函数的差；而相速度Ｃ为与波数无关，没有色散现象。方程描述了无衰减、无弥散的物理现 象。 ２、有耗散的模型( 4) u t + Cu x = v 2 u xx ? v 4 u xν2 & 0 ν4 & 0iｋ（ｘ－Ｃｔ）P P(3-52)其 中ν 2、ν 4 分别为二阶、四阶耗散系数。利用参数变易法，令解的一个波分量为：B B B Bｕｋ＝Ａｋ(t)ｅ 把波的分量代入方程(3-52),经推导得：dAk (t ) = ?(v 2 k 2 + v 4 k 4 )dt Ak (t )得到： Ak (t ) = Ak e0 ? ( v2 k 2 + v4 k 4 ) t45故有： u k = Ak e0? ( v2 k 2 + v4 k 4 ) te ik ( x ?Ct )B B(3-53)由此可知相速度Ｃ为常数，但振幅Ａｋ（ｔ）是随时间而衰减的。方程描述了有物理耗散 而无色散的波运动。另由式（２－５３）可见要求ν 2、ν 4 均为正值，否则解不稳定。 如果只考虑二阶耗散系数，则波分量为：u k = Ak0 e ? v2 k t e ik ( x ?Ct )经过Δｔ时刻后解为：2u k = Ak0 e ? v2 k则有： 放大系数：2( t + ?t ) ik ( x ?C ( t + ?t ))er = e ? v2 k ?t e ? ikC?t r = e ? v2 kＣ22放大系数模： 放大系数幅角： 相速度为：?t&1?Φ == ? kC?t有二阶耗散ν 2 的方程描述了有物理耗散的波动。耗散系数引起波幅的衰减，但相速度不发生 改变，放大系数幅角也与（３－５０）的相同。 ３、有弥散的模型( 3) ( 5) u t + Cu x = ε 3 u x + ε 5u x（３－５４）式中ε３、ε５称为三阶、五阶弥散系数，可取正或负值。 ｉｋ（ｘ－Ｃｔ） 令ｕｋ＝Ａｋ(t)ｅ ,代入式(3-54)推导得:B B B B B B B B P PdAk (t ) = i (?ε 3 k 3 + ε 5 k 5 )dt Ak (t )得到： 所以波分量：Ak (t ) = Ak0 e i ( ?ε 3k3+ε 5 k 5 ) tu k ( x, t ) = Ak0 e ik [ x ?( C +ε 3k2?ε 5 k 4 ) t ]２P波分量的振幅值不随时间变化，而相速度为Ｃ＋ε３ｋ －ε５ｋ ，是波数ｋ的函数。方 程（３－５４）描述了有弥散的波动。B B P B B P P４如果只考虑三阶色散，则解分量可写成： u k ( x, t ) = Ak e 经过Δｔ时间后解可写成： u k ( x, t + ?t ) = Ak e 放大系数： 放大系数模： 放大系数幅角：0 ik [ x ? ( C +ε 3 k 2 ) t ]0 ik [ x ? ( C + ε 3 k 2 )( t + ?t )]r = e ? ik ( C +ε 3k｜ｒ｜=１2) ?t?Φ = ? kC?t ? k 3ε 3 ?t46相速度为：C + ε 3k 2B B在三阶色散作用下波幅没有衰减，但相速度是波数的函数，当ε３＞０时，相速度增加， 波传播加快。与前两种模型相比，放大系数的幅角也增加了一项，这说明在同一位置ｘ，经 过Δｔ时间后波的相位滞后比无色散的波要大。 定义总动能： E1 = 定义总动量： E 2 =1 +∞ 2 u dx 2 ∫?∞∫+∞?∞u dxB B B B假定｜ｘ｜→∞时，ｕ、ｕｘ、ｕｘｘ…分别趋于零，这样考虑三种模型方程动能和动量的 守恒关系。 积分方程（３－５０）得：d +∞ udx + Cu dt ∫?∞+∞ ?∞，所以有：dE 2 = 0 。把方程（3-50）两 dt边同乘ｕ，然后积分，可得：dE1 = 0 。无耗散和色散的运动，其动能和动量是守恒的。同 dt样，仅有弥散时波动的动能和动量也是守恒的。对有耗散的方程经过同样的推导可得：dE1 &0 dt dE 2 =0 dt这表明有耗散作用的总动量守恒，而总动能是不断衰减的。三、数值耗散和弥散用差分方程逼近微分方程时引入了误差，有时这些误差项使计算结果的幅值衰减和相速 度发生变化，其作用相当于流动中的物理耗散和弥散，这种虚假的物理效应称作数值耗散和 数值弥散。 一维对流方程 （３－１） 描述的是既无耗散又无弥散的流体运动， 下面我们分别讨论 FTBS、 FTFS、FTCS 及蛙跳格式,来逼近它所产生的数值耗散和数值弥散。 １、ＦＴＢＳ格式 由式（３－３３）知，方程（３－１）的差分方程的等价方程为：?u ?u 1 ? 2 u C ? 2u +C + ? ?t ? ? ?x + Ο(?x 2 , ?t 2 ) 2 2 ?t ?x 2 ?t 2 ?x该式是差分方程的微分形式，其中截断面误差为：(3-55)R=C ? 2u 1 ? 2u t ? ? ? ? ?x + Ο(?x 2 , ?t 2 ) 2 2 2 ?t 2 ?x由式（３－５５）得:47?u ?u ?R = ?C ?t ?t ?R ? ?u ?u ? ? 2u [?C ( )? ? R ] = ?C = 2 ?t ?x ?t ?x ?t ?t 2 ?R ?R ? u C ? + = C2 ?t ?x ?x 2将上式代入到（３－５５）中,消去? 2u 得: ?t 2?u ?u C ? 2 u ?t 2 ? 2 u ?R ?R +C = ? ? ? x (C +C ? ) + Ο(?x 2 , ?t 2 ) 2 2 ?t ?x 2 ?x 2 ?x ?t ?x忽略高阶小量得:?u ?u C ? 2u +C = (?x ? ?tC ) 2 + Ο(?x 2 , ?t 2 ) ?t ?x 2 ?x(3-56)比较式(3-56)与原方程(3-1)可见用差分方程来逼近对流方程时,差分方程多出了若干项,其C ? 2u 主要部分为: ( ?x ? ?tC ) 2 ,与有耗散作用的方程(3-52)比较知该项为二阶耗散项,它使 2 ?x得差分方程的解具有数值耗散效应。 又由方程(3-1)可知原问题的解并没有物理耗散， 因此当 用差分方程来逼近对流方程时， 所得的解的幅值要偏小。 由式(3-56)知二阶数值耗散系数为：v2 =C (?x ? ?tC ) 2因为方程适定性的要求，故有ν2 ≥ 0 ，解之得差分格式的适定性的必要条件 为：０≤η≤１，由前面的稳定性分析知该格式的稳定性条件与此相同，由此可见通过上述 步骤也可以分析出差分格式稳定性的一些必要条件，这种方法称之为Ｈｉｒｔ启示性方法。 该方法可以用于方程组的非线性问题的分析。上述是分析了截断误差的主要部分，如果将后 面的高阶小量展开后可知，该格式同样存在弥散，但与耗散相比要小。因此一般通常分析主 要部分的数值效应而略去高阶小量。 ２、ＦＴＦＳ格式 相应于对流方程（３－１）的差分方程为（3-30） ，其等价方程为（３－３５） ：?u ?u 1 ? 2 u C ? 2u +C + t ? ? + ? ?x + Ο(?x 2 , ?t 2 ) = 0 ?t ?x 2 ?t 2 2 ?x 2采用与式(3-56)的同样推导方法得上式的简化等价形式为：C ?u ?u ? 2u +C = ? (?x + ?tC ) 2 + Ο(?x 2 , ?t 2 ) ?t ?x 2 ?x比较上式与式（３－５２）可知其二阶数值耗散系数为：(3-57)v2 = ?C (?x + ?tC ) 248由差分方程的适定性要求知：υ２≥０，得－１≤η≤０。由前面稳定性分析该条件为该格 式的稳定式充分必要条件。由此可见ＦＴＦＳ格式具有数值耗散效应。通过分析上述两格式 知， “逆风”格式具有数值耗散效应。 ３、ＦＴＣＳ格式 相应于对流方程的差分方程为（３－２９） ，其等价方程（３－３４）的简化等价形式为：B BC 2 ?t ? 2 u ?u ?u +C =? + Ο(?x 2 , ?t 2 ) ?t ?x 2 ?x 2与式（３－５２）比较知二阶耗散系数为：C 2 ?t v2 = ? &0 2与差分方程的适定条件不符， 故该格式为一不稳定的格式。 这在前面的稳定性分析中已证明。 ４、蛙跳格式 由式（３－３１）知差分方程的等价方程（３－３６）的简化等价形式为：?u ?u C?x 2 ? 3u +C =? (1 ? η 2 ) 3 + Ο(?x 3 , ?t 3 ) ?t ?x 6 ?x比较上式与式（３－５４）可见，该格式截断误差的主要部分为一个三阶弥散项，其三阶弥C?x 2 散系数为：ε 3 = ? (1 ? η 2 ) ,由此可见该格式解出的解的数值效应为弥散效应，反映在 6数值解上是不光滑有波动现象发生。四、混淆误差一阶空间微商?u 和 逼 近 它 的 一 阶 中 心 差 商 式 (3 ― 5) ， 分 别 作 用 于 一 波 分 量 ?xu ( x) = Ak e ikx ，则有：?u ?x= ikAk ejikx j= iku jik ( x j + ?x )(3-58)u j +1 ? u j ?12?x? Ak e = 2?x sin( k?x) = iku j ( ) k?xAk eik ( x j ? ?x )(3-59)从以上两式可见，一阶中心差商与一阶微商作用于同一波分量时得到的结果是不同的。中心 差商的作用得到的幅值比微商作用的结果有减小。其衰减比为：r=sin(k?x) k?x(3-60)考虑波长λ相对于空间步长很大的所谓长波，有λ&&Δｘ，那么 k?x =2π?xλ&& 1 。这时衰49减系数可作Ｔａｙｌｏｒ展开，得：1 r ≈ 1 ? (k?x) 2 6因为ｋΔｘ&&１，所以长波幅值很小，差商是微商好的近似，且当Δｘ→０时ｒ→１。 考虑可辨认的短波。如λ＝４Δｘ，则ｋΔｘ＝２πΔｘ／λ＝π／２，衰减比为 ｒ＝２／π，这时差商带来了很大的误差。因此对于和空间步长Δｘ接近的短波，差商是坏 的近似。ff1 = sin(πx?x)f2=0B Bxj+1j-1j图 3―10．网格无法辨认的短波至于更短的波长，格式则无法辨认。如图(3―10)所示，波长等于两个网格长度，即λ＝ ２Δｘ。这时用正弦波形的初值或全部为零的初值，所得的计算结果毫无区别。所以空间离 散无法辨认小于２Δｘ的波长，这时小于２Δｘ的短波分量会补充到长波分量中，从而使长 波分量发生畸变，甚至可能引起计算的不稳定。这种使短波分量补充到长波分量中，从而使 长波分量发生畸变而导致的误差称为混淆误差。因此在具体问题中选取Δｘ时，应保证波的 主要分量变得充分“长” ，即λ?x充分大。50第四章河道水流计算　天然河道水流常被认为一维流动,描述河道水流的基本方程为圣维南方程组,该方程是双 曲线型微分方程,有两类基本的求解方法,一是基于该方程的特征线形式的特征线法,二是基 于最初导出的偏微分方程的有限差分法,特征线法出现于本世纪初,是由马索（Ｍａｓｓａｕ Ｊｕｍｕｓ）对浅水方程式应用图解积分时提出的,其目的是在ｘ－ｔ平面上绘制特征线,并 在交点上确定因变量,是一种基本的并为数学家所偏爱的方法。 然而， 该法还有些困难以致阻 碍了它在明渠不恒定流中的广泛应用。 有限差分法随着计算技术的发展和计算机速度的提高， 已广泛应用于非恒定流计算，下面着重介绍适用于河道水流计算的有限差分法。　第一节蛙跳格式这是一个常用显格式,在距离和时间方面都采用中心差分,假设离散化如图４－１：? ?f f i j +1 ? f i j ?1 = ? 2?x ? ?t j ? ?f f ? f i ?j1 = i +1 　? 2?x ? ?x ? ? j ? f = fit j+1 (4-1)jj +1上述格式用于圣维南方程,便得到 Qi和 j-1 i-1 i i+1xZ i j +1 的显式表达式。　无侧向入流时完整偏微分方程为：　 　图 4―1 蛙跳格式离散化? ?Z ?Q ? B ?t + ?x = 0 ? ? QQ Q 2 B ? gA 3 ? ?Z Q 2 ?A ? ?Q + 2Q ?Q ? ? ? ? ? 2 ( ) z + gA 2 = 0 2 ? ? ? A ?x ? A K ? ?x A ?x ? ?t按式（４－１）对式（４－２）进行离散可得：　j j ? Qi ? Z i j +1 ? Z i j ?1 Qi + 1 + 1 =0 2?t 2?x 　　 　　 ?t 1 j +1 j ?1 j j Zi = Zi ? ? (Qi +1 ? Qi ?1 ) ?x Bi j(4-2)Bi j(4-3)　上式就是由连续方程得到的在ｉ点和ｊ时间层的水位显式表达式。从动量方程得：　51j j j Qi j +1 ? Qi j ?1 Q Q j ? Qi ? Q 2 B ? gA 3 j Z i + 1 1 ? Z i ?1 + ( ) ij i +1 ?( ) i ?x A 2?t 2?x A2 　　 j j 2 QQ A ( Z ) ? Ai ?1 ( Z i ) Q ? ( 2 ) ij i +1 i + g ( A 2 ) ij = 0 2?x A K(4-4)于是，对于在ｉ点和ｊ时间层的流量显式表达式为：　j j j Q Q j ? Qi ? Q 2 B ? gA 3 j Z i + 1 1 ? Z i ?1 +( ) Qi j +1 = Qi j ?1 ? [( ) ij i +1 i ?x A 2?x A2 　 j j 2 QQ A ( Z ) ? Ai ?1 ( Z i ) Q + ( 2 ) ij i +1 i ? g ( A 2 ) ij ] ? 2?t 2?x A K（4-5）阻力项也可以用下面差分形式：　 　　　 gAQQ K2=gAi j Qi j(K i )j2[φQi j +1 + (1 ? φ )Qi j ?1 ](4-6)这里φ是权函数，　０≤φ≤１，取φ＝０．５　则得：　 　　 gA 从而 Q ij +1QQ K2= 0.5 gAi j Qi j(K ij ) 2(Qi j +1 + Qi j ?1 )(4-7)的计算式为：　Qi　　j +1=g?t1Ai j Qi j{(1 ? g?t +1Ai j Qi jjQ j j )Qi j ?1 ? 2r ( ) ij (Qi + 1 ? Qi ?1 ) A (K i )2(K ij ) 2(4-8)+ r(式中: r =Q 2 B ? gA 3 j j Q2 j j ) ( Z ? Z ) + r ( ) i [ Ai +1 ( Z i j ) ? Ai ?1 ( Z i j )]} i i +1 i ?1 A2 A2?t 。 ?xt n+1 n n-1 ○ ○ ○ ○ x ○ ○ ○ ○图 4―2 蛙跳格式的两个独立网格 蛙跳方法是最先使用的方}