高等计算流体力学讲义(5)

§7. TVD格式

一、背景

1.求解线性波动方程utaux0aconst的经典差分格式 （1）一阶迎风格式（First order）

nnuuj1ja0n1n， ujujannujuj1a0其中t。上式也可以写为： x1nnununjj(fj1/2fj1/2)fjn1/2fjn1/2an(uj1unj)2an(unjuj1)2a2a2n(unj1uj)

n(unjuj1) （2）Lax－Friedrichs格式 （First order）

1unj1nnuj1unj1unj1uj12a0t2x

或

1nnnunu(ffjjj1/2j1/2)fjn1/2fjn1/2an1n(uj1un)(uj1unjj) 22ann(un(unjuj1)juj1)22nj1 （3）Lax－Wendroff格式 （second order）

1ununjjtuaunj12xa2tnn u2unjuj1。2j12x或

1nnununjj(fj1/2fj1/2)fnj1/2fjn1/2ana2nn(uj1uj)(uj1unj) 22ana2nn(ujuj1)(ujunj1)221

（4）Warming－Beam格式 （Second order）

1ununjjt1ununjj3uaanjn4unj1uj22x3unj4unj1unj22xa2tnnu2unj1uj a02j2 2xa2tnnu2unj1uj a02j22xt

2．二阶以上的差分或有限体积格式在间断附近的解可能会出现振荡。而一阶精度的格式通常会把激波“抹平”。参见图1和图2的（a），（b）。

图1 线性对流方程的解

图2

我们希望，构造一种二阶或以上精度的格式，使之可以比较准确的计算含有间断的流动，同时不产生非物理的数值振荡。这种格式称为“高分辨率格式”（如图1的（c），（d），图2的（c））。

2

3．单调格式

单调格式是一种没有数值振荡的格式，它的定义为： 定义：单调格式，设差分格式可以写为：

uin1H(uinr1,,uins) （1）

其中r ,s为非负整数。如果 j，

H0nuj，则（1）式称为单调格式。一阶迎风和Lax－

Friedrichs格式为单调格式。

4．单调格式的性质

u，利用单调格式（1）求得u，则：

★ 已知

nin1imax{uin1}max{uin}ii

min{uin1}min{uin}iin1min{unmax{unj}uij}jj

★线性波动方程

uua0的差分格式： txun1ibkuink （2）

kklkr是单调格式的充分必要条件是k，bk0。

★ Godunov定理 线性单调格式最多能达到一阶精度。 5．保单调性

nn1n假定ui是单调的，如果通过（1）式得到的ui具有与ui相同的单调性，则说差

分格式（1）式具有保单调性。容易证明，（1）单调格式具有保单调性，（2）不存在二阶或二阶以上精度的线性保单调格式。即如果（2）式为线性格式（bk为常数）且具有保单调性，则（2）式至多有一阶精度。

6．所谓高分辨率格式，是具有保单调性的，二阶或二阶以上精度的格式。容易知道，这类格式必定是非线性格式，即使求解的方程是线性的。本节我们只介绍其中的一种典型高分辨率格式：TVD格式。

二、TVD格式的概念

TVD（Total Variation Diminishing）格式即总变差减小（不增）格式。所谓总变差，是衡量函数的光滑性的一种指标，其定义为：

TV(u)u(x)dx （3）

 在离散点上，u{ui}，有：

3

nnTV(u)niuni1uin （4）

可以证明一维线性和非线性标量守恒律的解满足，

TV(u(t2))TV(u(t1))t2t1

这种性质称为TVD性质。既然一维标量守恒律的解析解满足TVD性质，要求其数值解也具有TVD性质就是很自然的了。 定义：TVD格式

uin1H(uinr1,...,uins) （5）

称为TVD格式，如果

TV(un1)TV(un)n

容易证明：TVD格式具有保单调性，所以在间断附近不会出现非物理振荡。

三、TVD格式的构造

1． TVD格式的充分条件和构造原则

设一维标量守恒律

uf0的差分格式可以写成： txuin1uinC1ui2i12D1ui2i12 （6）

nn其中ui12uni1u。C1，D1不仅与x，t有关，而且与物理量ui，ui1等有关。

nii2i2Harten证明了下列定理：

定理： （6）式是TVD格式的充分条件是：

C，， （7）

根据这一条件，我们可以构造精度高于一阶的TVD格式。下面介绍构造TVD格式的思路。 方程

i120Di1200Ci12Di121uf0；fuau的差分格式为： txuin1uintff11ixi22TVD1i2 （8a）

如果(8)具有TVD性质，记此时的数值通量为f。即

uin1uintTVDTVD （8b） ff11iix224

如果要求格式具有高于一阶精度，则fTVD1可以写成：

i2TVD1i2LD1i2ffHIL01f1f1iii222 (9)

L0HI其中fi1/2是某一阶格式的数值通量，fi1/2是二阶格式的数值通量。i1/2称为通量器

（flux limiter）。我们希望选择合适的1，使格式在光滑区具有二阶精度。当采用三点格式

i2时，fL0HIi1/2，fi1/2可写成如下一般形式：

fL0nni1/20aui1aui1

fHInni1/20aui1aui1 fL0i1/2可以有多种取法，如取

0121s，1121s 时，fL0i1/2为一阶迎风格式的数值通量，其中

ssign(a)。 当取：

011c12c，12c1c

时，fL0i1/2对应的差分格式为Lax-Friedrichs格式。注意，上式中catx。 二阶精度的三点格式只有Lax-Wendroff格式，即：

1021c111c，2

把（10）式代入（9）式，有

fTVDi12[0(00)i12](auni)[1(1n1)i12](aui1) 把（13）代入（8b）式，得：

un1iuniCui1Dui122 Cc[0(00)i1]2

Dc[1(11)i1]2

5

10）

11）

12）

13） 14）

（（（（（22

由（14）式，我们希望选取适当的0、0、，使得TVD条件（7）式得到满足，且格式为二阶精度。

ui1uinuin1ui1uin1uin2、迎风型TVD格式

如果我们取记

L0fi12为一阶迎风格式的数值通量，则得到的差分格式称为迎风型TVD格式，

ii1212。则TVD格式的数值通量可以写为：

fTVDi12n1nau(1c)ui11aii222aun1(1c)una11i1ii222，此时：

a0

a0当a>0时，

L0nfi1aui201,10,c0Cc[1(01)i],Dc11121012(1c),12(1c)1 （15） 2（14）式可以改写为：

ˆu1 （16） uin1uinCi2ˆCD/r,ru1u1 Cii22Harten的TVD条件此时为：

1Oc[1(01)i11i1]1 （17）

22r 由（17）式，

i12应是r的函数。显然我们要求

22c[1(01)i1]c1i1我们要求是有界的。所以可以对

11c[1(01)i1] （18）

2r112引入下列条件：

BiT12i,r （19）

B、T与i和r无关。（19）式可以等价地写为：

c[1(01)T]c[1(01)i1]c[1(01)B] （20）

2 6

易知下式是（18）式成立的充分条件：

c[1(01)T]c1i1 对（21）式左边不等式进行分析，可知：

211c[1(01)B] （21） rifr0ifr0 （22）

i(r)12LLLrT对（21）中右侧不等式进行分析，有

1r1。

1ri2RrRrifr0 if r0 （23）

RrB1crc1

利用（19），（22），（23）式可以确定

i12r的具体形式。

在讨论

i12r的具体形式之前，先讨论一下方程

uf0tx

中a0的情况。此时，（14）式可以改写为：

fau

ˆuuin1uinDi12 （24）

ˆDC DrHarten的TVD格式的充分条件为：

uri1212u

i0c111101ii22引入约束条件

11r （25）

B为了满足TVD条件，需有：

i12T

（26）

c111Tc0i1211c111Br （27）

注意到此时a0，c0，如果把c用c代替，则（27）式与（21）式在形式上是相同

7

的。但是注意到a0和a0两种情况r的定义不同。把（27）式应用于

i11i2处（把2i换为

12）。则（21），（27）式可以写成统一的形式。因此（22），（23）式也可以写成既

i12处）

（28）

使用于a0，也适用于a0的一般形式（在

BrT

rLrLrr0

r0 r0r0 21cr

（29

RrrRr (30)

LrT (31)

2RrBrc

(32)

rupwindlocaluinuin1nui1uinnnuui2i1nnuui1i

a0a0

(33)

为了确定(r)的具体形式，必须首先规定B，T的值。一个常用的取法是

B0， T2 （34） 1c2时fiTVD1/2为一阶1c注意到当

=B=0时，fTVD为一阶迎风格式的数值通量。

Ti1/2T的取法具有一定的任意性，“顺风”格式的数值通量。必须注意，不同的取法，对于r平面上TVD格式对应的区域有直接影响。如Sweby取T2，他得到的TVD区域

T2的情况不尽相同。把（34）式代入（31）、（32）有 1c 8

L(r)0，R(r)即（28）~（30）式可以化为：

2r （35） c0(r)2 （36） 1c(r)002rc2rcr0 （37） r0r0

（38）

(r)由（36）~（38）式，知：

r0220（r）min(,r)1cc （39） r0时 (r)0r0时满足（39）式的(r)对应的格式，均为TVD格式。由（39）式，(r)的取法并不唯一。 另外，如何使TVD格式达到二阶精度呢？容易看出其条件为(1)1。注意到(r)1，（8b）为Lax－Wendroff格式，这是一个二阶精度的格式。容易验证，当(r)r时，（8b）为Warming－Beam二阶迎风格式。所以，当(r)夹在直线(r)1和(r)r之间时，格式（8b）可以看作两个二阶格式的加权平均，具有二阶精度。Sweby建议：(r)应尽量位于(r)1和(r)r之间。

下面介绍几种常用的通量器(r)： ULTRABEE

r00c2 ubr 0r1cc2cr1c1cSUPERBEE

9

(r)2r (r)r

(r)1

ub sb vl mb

02rsb(r)1r2VANLEER

r00r121r121r2 r2

r00vl(r)2r1r r0

VANALBADA

r00va(r)r(1r)1r2 r0

MINBEE(MINMOD)

0mb(r)r1r00r1 r1

不同的器的性能有一定差别。(r)越大，格式的耗散越小，对激波的分辨率越高。但是，(r)过大，数值结果中可能出现一些异常现象。如正弦波可能变为方波的形状，甚至

10

可能导致计算失稳。对于定常流动的计算，使用光滑的器容易收敛到较小的残差。

3、非线性标量守恒律的TVD格式

uf0考虑方程tx，设它的差分格式可以写为：

uin1uin 令

tfi1/2fi1/2 （40） xt(fi1/2f(ui))(fi1/2f(ui))，Di1/2，Ci1/2 （41） xui1/2ui1/2则（40）式可以改写为：

uin1uinCi1/2ui1/2Di1/2ui1/2 （42）

迎风型TVD格式的数值通量可以写成下面的形式，

LOHILOfi1/2fiTVDfff1/2i1/2i1/2i1/2i1/2。 （43）

迎风格式的数值通量为

LOfi1/21~fifi11aui1/222i1/2， （44）

其中：

~ai1/2fi1fiuui1ifui1/2ui1ui0ui1ui0 （45）

Lax－Wendroff格式的数值通量为

fHIi1/2~2t11ai1/2(fi1fi)(ui1ui)22x。 （46）

~把（43）－（46）式代入（42），可得（42）中C，D的表达式。容易验证C，D的表达式与（15）式形式上是相同的， 只须以a代替a即可。由此可以看出，线性方程的TVD格式的推导，同样适用于非线性方程。即ψ的形式线性方程和非线性方程是相同的。

4、非线性守恒方程组的TVD格式

可以证明，非线性双曲型守恒方程组（系统）的解析解本身不具备TVD的特点。所以

要求求解守恒方程的格式具有TVD的性质是不需要的，甚至是有害的。实际上，有人已经

11

证明，不存在求解非线性守恒系统的TVD格式。所以所谓非线性守恒系统的TVD格式只是对标量方程TVD格式的形式上的推广。

uf0tx先考虑标量方程，它的TVD格式的数值通量为

LOHILOfiTVDfff1/2i1/2i1/2i1/2i1/2 （47）

LOHIfi1/2和fi1/2的表达式见（44）、(46)式。（47）式可以进一步整理为：

fTVDi1/2~2ta11~(fi1fi)[i1/2i1/2(1i1/2)ai1/2](ui1ui)22x （48）

UF0tx （49）

下面把（48）式推广到非线性双曲型守恒系统的情况。考虑守恒律系统：

~~A(Ui,UR)存在。且A假定在i1/2处，Roe的线性化矩阵的特征值为

(j),(j1,2,...,m；m为方程的维数)，相应的左右特征向量为i1/2(l(j))i1/2;(r(j))i1/2,(j1,2,...,m)。则把（48）式推广到求解（49）式的方法为：

（1）、把（48）式中的ui1ui和

1(fi1fi)用特征变量代替，即： 2ui1ui(l(j))1(Ui1Ui)i(fi1fi)/2(l(j))2i12(Fi1Fi)/2

~(j)代替。 Ai1/2，用矩阵的特征值（2）、把（48）式中的ai1/2（3）、（33）式中的r用特征变量的差分来定义，即：

r(j)(l(j))1(UinUin1)i2(j)0,(l)(UnUn)i1i1i2(j)i12 （50） nn(l(j))i1(Ui2Ui1)2(j)0,(l(j))1(Uin1Uin)i12i2(j)，由特征变量的差分代替标量方程的（4）、上述（1）－（3）相当于对每一个特征值i1/2(j),(j1,2,...,m)进行上述替代后，可得第j个特征变量对应的特差分。在对每个特征值

~征数值通量为

12

1TVD(l(j))i1/2Fi1/2(l(j))i1/2(Fi1Fi)2 。 (j)2()i1/2t1[i1/2(1i1/2)((j))i1/2](l(j))i1/2(Ui1Ui)2x

最后，再把数值通量由特征变量的空间投影回物理变量的空间。利用关系

()r(j)(l(j)())，

mj1其中()代表任意的m维列向量，可知，非线性双曲型守恒系统的迎风型TVD格式的数值通量可以写为：

FTVD1i21(FiFi1)2(j))2t(1i1m(j)(j)2(j)(r)(1(r))111ii2j1i2x22 ()(rj1(j))1ii22（51）

j其中：i12(l(j))1i2(Ui1Ui)，r(j)由（50）式计算，(r(j))i12~A(Ui,UR)对应的是

右特征向量。

思考：如何把TVD格式推广到求解问题？

§8、有限体积方法的高阶数据重构，MUSCL插值

一、一维问题的MUSCL－Hancock方法

MUSCL－Hancock方法是基于Riemann问题的精确解和近似解的时、空方向均为二阶精度的有限体积方法。

13

1．数据重构

[i11,i]22中是线性分布，则在

为了达到空间二阶精度，我们假定这个控制体中

U(x,tn)在控制体

niUi(x,tn)Ui(xxi)x （52）

nnnnn(UUiii1i1)/2,或其中是反应物理量坡度的矢量，最简单的取法可以为nnn(UUii1i)等。容易证明

xix121i2Ui(x,tn)dxxniUin （52）

n这一条件称为相容条件。（51）式称为对U的一阶数据重构。i的确定方法在稍后详细讨n论，现在假定i已经确定，则可以得到单元交界面的左右两侧物理量的值为：

ULi12niUi(x1,tn)Ui22ni （53）

URi12Ui1(xi12,tn)Uni1ni12 （54）

2．

tnt2时刻界面物理量的计算

为了达到时间方向二阶精度，把方程推进半个时间步，计算tnt/2时刻的单元两侧

LR的物理量Ui1/2，Ui1/2：

L1tLRU1UL1F(U)F(U)11 （55） iiii2x2222R1tRLRU1U1F(Ui1)F(Ui1) （56） ii2x2222 14

3. 求解Riemann问题

LR令ULUi1/2，URUi1/2。把UL、UR做为Riemann问题的初始条件，精确或近

似地求出Ui1/2(x/t)，则数值通量Fi1/2为：

Fi1F(Ui1(0))22 （57）

这样，求解

UF0tx

的有限体积格式可以写为：

Uin1UintFi1Fi1x22

（58）

通过下一小节的分析，我们可以证明，该方法时间和空间均可达到二阶精度。

4. MUSCL插值

n1中简述了数据重构的方法，下面对i的确定方法做进一步讨论。当采用1中所述简

单重构方法时，格式为二阶精度，但在间断附近会产生数值振荡。为了避免数值振荡，必须

n对i进行某种。这种思想最初源于Van Leer的Montone Upstream-Centered Scheme for

Conservation Laws。所以通常简称为MUSCL插值（或重构）方法。MUSCL插值的方法有很多，这里介绍一种根据TVD格式中的通量器（Flux Limiter）方法构造的MUSCL重构方法。通过这种方法，对于标量线性双曲型方程，可以构造出和前面讨论的TVD格式完全等价的MUSCL－Hancock格式，由此我们知道，MUSCL－Hancock也具有时空二阶精度。

考虑标量线性双曲型方程

uua0tx

具有TVD性质的迎风格式的数值通量可以写为：

fTVD1i2n1nauii1(1c)ui1a222aun1(1c)una11i1ii222a0 （59）

a0这相当于Riemann问题的解为：

15

ui12n1nu(1c)ui11ii222un1(1c)un11i1ii222a0 （60）

a0在采用MUSL－Hancock格式时，（60）相当于（55）（56）式中，tnt/2时刻的左右状态为：

u~Ln1i1uii1(1c)un1222i2 u~Rn1i1ui12i1(1c)un122i2 又由

u~LatLi1uL2i122xui1uR2i12 u~Rati1uR2i122xuLuRi32i12

考虑到（53），（54）式，即：

nnuLniRiu1i1unii22u

i22

nnuRnni1uii1L122u

i3uii1122

可得：

uL1/2uni1i2(1c)ni uRn1i1/2ui12(1c)ni1 比较（61），（62）和（），（65）式，有：

nii1una02i12ni1i1uni1a0

22或

na0ni1u12i2i ni1ui1a0 22

16

61） 62） () (65) 65）

（ （ （n在这里称为坡度器（slope limiter）。当的取法不同时，i的表达式也不相同。一

种常用的“slope limiter”是：

max[0,min(un1,un1),min(un1,un1)]iiiin2222innnnmin[0,max(ui1,ui1),max(ui1,ui1)]2222whenwhenun10ini212u0 （67）

容易证明，当（67）式中1时，相当于（66）式中 取MINMOD器。2时，相当于（66）式中

 取SUPERBEE器。为了以后的表达方便，可以把（66）

，（67）

式抽象地写成in(uin1/2,uin1/2)。是坡度的平均算子。

把上述MUSEL插值方法推广到方程组时可以简单地把（66）或（67）式应用于守恒变量的差分Uin1/2的各个分量。另外一种计算效果更佳的推广方法是对特征量的差分进行坡度，即：

n，如U（1）估算U1i2ni121n为Un,Un的Roe平均。 (Uin1Uin)或者取U1i1ii22n，计算（i＋1/2）处的特征值和特征向量。特征值为（2）根据U1i2(k),(k1,2,...,m；m为方程的维数)，相应的左右特征向量为i1/2(l(k))i1/2;(r(k))i1/2,(k1,2,...,m)。

（3）计算（i＋1/2）左右j＝（i－1,i，i＋1,i＋2）各点的特征变量。v(k)j(l(k))i1/2Uj，

k1,2,3,...,m。

（4）计算特征变量的坡度。

n,(k)imax[0,min(vn,(1k),vn,(1k)),min(vn,(1k),vn,(1k))]iiii2222n,(k)n,(k)n,(k)n,(k)min[0,max(vi1,vi1),max(vi1,vi1)]2222max[0,min(vn,(1k),vn,(3k)),min(vn,(1k),vn,(3k))]iiii2222n,(k)n,(k)n,(k)n,(k)min[0,max(vi1,vi3),max(vi1,vi3)]2222n,(k)L,R。 1i2whenwhenwhenwhenvn,(1k)0iv2n,(k)1i20

vn,(3k)0in,(k)i1v2n,(k)3i20

（5）计算（i＋1/2）处的特征变量(v（6）计算（i＋1/2）处的守恒变量

)(U

nL,Ri1/2)r(k)(vn,(1k))L,R （68）

k1i2m17

L,R（7）在所有界面计算出(Uin后，即可按照（55）～（58）的步骤推进求解。 1/2)

二．问题的MUSCL-Hancock方法

下面把MUSCL-Hancock方法推广到问题。为了使问题简化，我们只考虑网格是充分光滑的情况，如下图所示：

对于网格分布不光滑和非结构网格的情况，也有适当的处理方法，这里不准备介绍。

1．重构

(x,y)(x,y) 使得指标j 相同的控制体的中心

在光滑网格上，我们可以引入坐标变换点处为

j的坐标线。指标i相同的控制体的中心为

i 的坐标线。这样，当在控制

体i,j内采用线性重构时，有：

Ui,j(x,y,tn)Ui,j(,,tn)Ui,j()i,j(j)(i)()i,j (69)

11i,ji,j2 和 2 处守恒变量的值，事实上，我们感兴趣的只是 因此，不必关心 (69) 式

中 , 的具体计算方法，由 (69) 式有：

UiLUij12i,j1,j2R1Ui1,jUij2i,j2 L1Ui,j1Uij2i,j2R1Ui,j1Uij2i,j2 (70)

，的计算方法与一维问题类似，例如：

18

i,jUi1,jUi,j,Ui,jUi1,j (71)

2.ttnt/2时刻界面物理量的计算

t4UFiGjnlSS2ij1

(72)

22111Si,ji,jSi222,j时 其中，。当

LLFiGjnlS(F(Ui1,j)iG(Ui1,j)j)ni1,jli1,j22

当Si12,j时，

RRFiGjnlS(F(Ui1,j)iG(Ui1,j)j)ni1,jli1,j2222

当Si,j12时可以类似写出。在求得U后，

L1UL1U Ui,ji,j2222 (73)

R1UR1U Ui,ji,jn1Uij3. 求解Riemann问题计算

L,R2

1U1当作左右状态，通过求解扩张的一维Euler方程的Riemann问题可以计算出把Ui,ji,jL,R2数值通量，从而求得

n1Uij。具体做法参见§4。

§7. 半离散的TVD格式

一、半离散TVD格式的定义

定义：半离散格式

ui1 （1） ff11iitx22称为TVD格式，如果

19

此时，记数值通量为fi12dTV(u)0。 dtfTVD1。

i2二 半离散TVD格式的充分条件

设一维标量守恒律

uf0的差分格式可以写成： txui1(C1u1D1u1) （2）

iiiitx2222其中u12iuni1u。C1，D1不仅与x，t有关，而且与物理量ui，ui1等有关。

nii2i2nn

则（2）式是TVD格式的充分条件是：

Ci120，

Di120 （3）

根据这一条件，我们可以构造精度高于一阶的TVD格式。下面介绍构造TVD格式的思路。 方程

iuf0；fuau的某一差分格式具有TVD性质，记此时的数值通量为txfTVD1。即

2ui1TVDTVDff （4） 11iitx22如果要求格式具有高于一阶精度，则fTVD1i2可以写成：

LOfTVD1f1i2i2HILOff111 (5) iii222其中fi1/2是某一阶格式的数值通量，fi1/2是二阶格式的数值通量。i1/2称为通量器（flux limiter）。我们希望选择合适的L0HIi12L0HI，使格式在光滑区具有二阶精度。当采用三点格式

时，fi1/2，fi1/2可写成如下一般形式：

L0nnfi1/20aui1aui1

HInnfi1/20aui1aui1 （6）

把（6）式代入（5）式，有

20

fiTVD[0(00)i1](auin)[1(11)i1](auin1) （7） 1222 把（7）代入（4）式，得：

uin1uinCui1Dui122 （8）

Cc[0(00)i1]2

Dc[1(11)i1]

2

ui1uinuin12

ui1uin1uin2，cat x 我们希望选取适当的0、0、，使得TVD条件（3）式得到满足，且格式为二阶精度。 实际上，半离散格式具有TVD性质的条件比全离散格式要宽松，所以，全离散TVD格式在半离散意义上也具有TVD性质。

三、迎风型半离散TVD格式

我们考虑迎风型格式，fi1/2取为：

L00其中ssign(a)。

11s111s22， （9）

二阶精度的三点格式取中心格式，即：

011，1 （10） 22为了方便，我们使用全离散格式的TVD条件作为半离散格式的TVD条件。这样做的原因

是：1）由此得到的半离散TVD格式在时间方向作简单的前差离散时，得到的格式有TVD性质；2）便于应用TVD－Runge－Kutta方法，因为TVD－Runge－Kutta方法应用的条件是对应的时间前差离散格式具有TVD性质。

当a>0时，

L0nfi1aui2，此时：

01,10,c0Cc[1(01)i],Dc11121012,121 （11） 2（8）式可以改写为：

ˆu1 （12） uin1uinCi2 21

ˆCD/r,ru1u1 Cii22Harten的TVD条件此时为：

1Oc[1(01)i11i1]1 （13）

22r 由（13）式，

i12应是r的函数。显然我们要求

c[1(01)i1]c1i111c[1(01)i1] （14）

22r2我们要求是有界的。所以可以对

112引入下列条件：

Bi1i,r2T B、T与i和r无关。（15）式可以等价地写为：

c[1(01)T]c[1(01)i1]c[1(01)B] 2易知下式是（14）式成立的充分条件：

c[1(01)T]c11i12r1c[1(01)B] 对（17）式左边不等式进行分析，可知：

)Lifr0i1(r2Lifr0 LrT2r。

对（17）中右侧不等式进行分析，有

Rrifr0i1r2Rr if r0 2(1c)RrBcr 1r利用（15），（18），（19）式可以确定

i2的具体形式。

r在讨论

i12的具体形式之前，先讨论一下方程

uftx0

fau

中a0的情况。此时，（8）式可以改写为：

un1iuniDˆui12 22

15）

16）

17） 18）

19）

20）

（ （ （ （ （ （ˆDC DrHarten的TVD格式的充分条件为：

uri1212u

i0c111101ii22引入约束条件

11r （21）

B为了满足TVD条件，需有：

i12T

（22）

c111Tc01i211c111Br （23）

注意到此时a0，c0，如果把c用c代替，则（23）式与（17）式在形式上是相同

的。但是注意到a0和a0两种情况r的定义不同。把（23）式应用于

i11i2处（把2i换为

12）。则（17），（23）式可以写成统一的形式。因此（18），（19）式也可以写成既

i12处）

（24）

使用于a0，也适用于a0的一般形式（在

BrT

rLrLrr0

r0 r0

（25）

rRrRrr0

(26) (27)

LrT2r

2(1c)RrBr c (28)

23

rupwindlocaluinuin1nui1uinnnuui2i1nnui1ui

a0a0

(29)

为了确定(r)的具体形式，必须首先规定B，T的值。我们取

把（30）式代入（27）、（28）有

B0， T2 （30）

2(1c)L(r)0，Rrr （31） c即（24）~（26）式可以化为：

0(r)2 （32）

(r)00r0 （33） r0r0

（34）

2(1c)rc(r)2(1c)rc由（32）~（34）式，知：

r02(1c)r0,0（r）min(2,r) c （35）

r0,(r)0且二阶精度的条件为(1)1。Sweby建议：(r)应尽量位于(r)1和(r)r之间。根据上述条件，可以构造各种满足TVD条件的器。注意到，（35）式中，半离散格式的TVD条件与CFL数c有直接关系，c1时，只有一阶格式才可能满足TVD条件，所以，在实际应用中，我们应取c1（如c1/2或者更小，全离散格式的多数器都可以应用）。

四、时间离散

常用Runge－Kutta方法离散，如具有TVD性质的三阶Runge－Kutta方法：

24

ui(0)uinui(1)ui(0)tQi(u(0))(0)(1)(1)ui(2)314ui4(uitQi(u))， (0)(2)ui(3)12tQi(u(2)))3ui3(uiuin1ui(3)其中Qi

1TVDTVD ff11iix22五、基于MUSCL插值的有限体积半离散格式

uf0,fau先考虑标量线性守恒律tx的情况。此时，在i1/2（x0）

左右给定Riemann问题的初值u,u，则Riemann问题的解为：

Lu,if(x/t)au(x/t)R，

u,otherwiseLR 或者

所以，如果，

Lu,a0u(0)R。

u,otherwiseTVDfi[0(00)i1](auin)[1(11)i1](auin1)， 1222或者具体说，如果

1s1sTVDnnfi[(1)](au)[(1)](au111ii1)， i2i222222则该通量对应的Riemann问题的初值为：

uinuin1i1/2(n)nuui1iuLuin(uin1uin)2。 nnu2ui1i1/2(i)nnuui1iuRuin1(uin1uin)2从有限体积方法的角度，上式相当于

25

uinuin1i1/2(n)nuuni1iuLuin(uin1ui)2。 （36） nnuu1i1/2(in2i)nuuni1iuRuin(uin11ui)2这样，我们就得到了一种满足TVD条件的重构方法，称为MUSCL插值。这个方法可以直接推广到求解非线性的标量守恒律和非线性守恒系统（方程组）。对于方程组，我们可以直接推广上式为：

UinUin1i1/2(n)nUi1UinULUin(Uin1Ui)2 （37） nnUUi1i1/2(in2)nUi1UinURUin(Uin11Ui)2注意到上式右端的运算是对每一个分量分别进行；也可以仿照我们在MUSCL－Hancock格式中所介绍的那样，对特征变量进行。在通过（36）、（36）获得Riemann问题的左右状态后，数值通量通过求解精确或者近似的Riemann问题确定。

为了给大家一个对MUSCL插值的直观印象，我们考虑（36）式中器为Van Leer器的情况，此时，（36）是化为：

nnmax[(uinuin1)(ui1ui),0]uun(uin1ui1)LniuRuin1这相当于第i个单元的坡度为：

max[(uni2u)(uu),0](uin2u)ni1ni1nini，

nn2max[(uinuin)(uu1i1i),0]。 inn(ui1ui1)

26