中科院计算流体力学最新讲义CFD116讲差分方法4

w1
2w21w2
0
1w1
0
w3
w2
实质区别： f(W) 是自变量 W 的线性函数， 而 f(U) 是自变量 U 的非线性函数！
U
U
使用U做自变量的优点： 物理意义鲜明 （质量密度、动量密度和能量密度），守恒性好
使用W做自变量的优点： Jocabian矩阵为线性矩阵
思考： 如果在CFD计算中，使用W替换U做自变量会怎样？
方案B: 设计一种迭代算法
令
Qn Un1Un
两端同时添加显式项
右端项，已知
Qnt[f1n1 f2n1]0 x y
Q n t[ (f1 n 1 f1 n) (f2 n 1 f2 n)] t[ f1 n f2 n]
x
y
x y
tRHSt[f1n f2n] 是已知项，可采用某种差分方法显式计算得到 x y （对算法无限制）
1） 若高精度逼近
u x
，
必然利用多个基架点
2） 如果该基架点内函数有间断，会导致振荡
3） 间断不可能处处存在
4） 把基架点分成多个组（模板），
每个模板独立计算j点导数的逼近。
—— 得到多个差分
{j-3, j-2,j-1,j,j+1,j+2}
5）根据每个模板的光滑程度，设定权重
{j-3,j-2,j-1,j}； {j-2,j-1,j,j+1}； {j-1,j,j+1,j+2}
本PPT 第5页）。推广到高阶后，虽不再保证严格的特征方向， 但仍优于采 用算数平均方法。
Roe 格式的不足： 本身精度只有一阶； 推广到高阶后，特征方向无法严格保证 ； 推广到二维或三维后，特征方向无法严格保证，出现振荡。
Copyright by Li Xinliang
11
深入讨论
关于 f(源自文库) 与 f(W)
2
u
w3
H
HEp
则：
w12
U
U(W)
w1w2
1
w1w3
1 2
w22
fR fLC W )((R W W L ) W(W RW L)/2
目的： 使得F(w)是W二次齐函数 （增长 率为线性函数）
f (W)
1
w1w2
w1w3
1 2
w22
w3w2
二次齐函数！
f(W)w 12
w1
1
含义： 左、右两个状态点的某种平均 （称为Roe平均，为密度加权平均） 该状态点对应的增长率（矩阵）为平均增长率（矩阵） 实际上是一种“等效平均”。 效果优于简单的算数（或几何）平均。
Copyright by Li Xinliang
9
✓ Roe 格式的计算步骤 （半离散） Uf(U) 0
t x
j-1
线性化，用平均变化率代替 (j,j+1)之间的变化率a(u)
a~j1/2
f
(uj1)f (uj uj1 uj
)
f 1
x
j
( x
fj1/2
fj1/2)
fj1/2 ffjj1
a~j1/2 0 a~j1/20
fjn 1 /21 2 [fj fj 1]1 2a ~ j 1 /2(uj 1 uj)
a~ j 1 / 2
整理
Q i[jI x t(A i ,j A i ,j) y t(B i ,j B i ,j) ] x t[A i 1,jQ i1,j A i 1,jQ i1,j]
A i, j
A i,
j
* A
I
B i, j
B i, j
* B
I
y t[B i ,j1Q i,j1B i ,j1Q i,j1] tRHS
Copyright by Li Xinliang
14
隐式修正
显式部分
Q nUn1UnΔn U
Q nt[(f1n1f1n)(f2n1f2n)]t RHS
x
y
若隐式修正为0，则为显格式
A (U )f1(U ),B(U )f2(U )
U
U
f 1 n 1 f 1 n A U n O ( t 2 ) f 2 n 1 , f 2 n B U n O ( t 2 )
0
A(U)fU (U)
1(3)(u2)2
2
u1
uu212u3
(1)(u2)2
u1
1
(3)(u2)
u1
u3 3(1)(u2)2
u1 2
u1
0
1
(uu12)
自变量U的非线性函数
f (W)
1
w1w2
w1w3
1 2
w22
w3w2
干净的二次 齐函数
自变量W的线性函数！
C(W f)w (W )w 12w3
6） 对多个差分结果进行加权平均 。光滑度越高，
权重越大。如果某模板存在间断，则权重趋于0；
如果都光滑，则组合成更高阶格式。
五个基架点被分成三个组
Copyright by Li Xinliang
3
§ 6.1 Roe格式
1. 单方程的Roe格式
u f (u) 0 t x
u a(u)u 0
t
x
非线性情况
2阶
U(1)UntQ (Un) Un11/2Un1/2U(1)1/2tQ (U(1))
4阶
U(1) Un 1/2tQ(Un) U(2) Un 1/2tQ(U(1)) U(3) Un tQ(U(2)) Un11/3(Un U(1) 2U(2) U(3))1/6tQ(U(3))
3 阶 （TVD型）
U(1) Un tQ(Un) U(2) 3/4Un 1/4U(1) 1/4tQ(U(1)) Un11/3Un 2/3U(2) 2/3tQ(U(2))
j
j+1
已知n时刻所有网格点上的物理量，对于j点：
1） 利用差分格式计算UR,UL 2） 采用Roe平均公式（5）计算Roe平均值 U
3） 将Jacobian矩阵A( U ) 进行特征分解: A(U)S1ΛS 计算 S1,Λ,S
4) 计算 A ~(U R,UL)S1ΛS
5）计算 f j 1 /2 f ˆ (R U ,U L ) 1 2 [f R ) ( fU L ( ) U 1 2 ]A ~ (R U ,U L ) (R U U L )
0
1
C(W ) w
w3
2w2
w2
w1
0
w3
w2
增长率为线 性函数！
Roe点
Roe点为：
UU(W)
中点处的斜率即为平均斜率。
Copyright by Li Xinliang
8
最终： ~ A (U R ,U L)A U )(
平均增长率（矩阵）
其中 U 如下计算：
简单易记：
[( L R)/2]2
( L R)/2
u( LuL RuR)/( L R)
（5）
H( LHL RHR)/( L R)
其他量（如压力、温度、音速等）用这三个量计算
p1(H1u2)
2
2
c
(
u2 1)(H )
2
HEp
u( LuL RuR)/2
三维情况下，还有
v( LvL RvR)/(L R) w( LwL RwR)/(L R)
更高阶 ……
Copyright by Li Xinliang
13
2. 隐格式算法简介 （以二维Euler方程为例）
1） 原理介绍
Uf1(U)f2(U)0 t x y
时间离散后
U n1U nf1(U n1)f2(U n1)0
t
x
y
（1）
方案 A： 直接将（1）进行空间离散，得到 Un+1 的代数方程组 困难： 大型非线性方程组，求解困难
(D UL)ΔUtRHS 近似LU分解 D U L (D L )D 1 (D U ) L 1 U D
x
y
矩阵分裂
Q n t[A Q nA Q nB Q nB Q n] tRHS x x y y
1阶迎风格式离散
迭代收敛后，不影响精度
Qij xt[Ai,jQi,j
A Q i1,j i1,j
A Q i1,j i1,j
Ai,jQi,j]
yt[Bi,jQi,j Bi,j1Qi,j1Bi,j1Qi,j1Bi ,jQi,j]tRHS
jj j(rj),
rj
uj uj1 uj1uj
二阶精度区
TVD区
Copyright by Li Xinliang
二阶精度TVD区（二 者交集）
2
知识回顾：WENO格式
基本思路
u au 0 a0
x x
f j a 1 f j 3 a 2 f j 2 a 3 fj 1 a 4 f j a 5 f j 1 a 6 f j 2
u1
U u1 u
u
3
E
f1 u
f f1 u2 p
f3
u(E p)
新变量
w1
1
W
w2
1/
2
u
w3
H
f(U)
u2
1 2
(3
u 2u 3 u1
)
u
2 2
u1
(
1)u 3
1
(
1)
u
3 2
2
u
2 1
虽然是“一次齐 函数”但有变量 在分母上
中科院计算流体力学最新讲义CFD116讲差分方法4
知识回顾： 单调、保单调和TVD
概念： 网格Reynolds数 R exRex 单调格式、保单调格式及TVD格式
Harten定理： 正系数原则
TVD
单调
保单调
TVD格式= 1阶迎风+ j *（修正项）
j u j 1 /2 u j(r)u (j 1 u j)/2
Copyright by Li Xinliang
16
Q i[ 1 j x t* A y t* B ] [ x tA i 1 ,j Q i 1 ,j y tB i ,j 1 Q i ,j 1 ] [ x tA i 1 ,j Q i 1 ,j y tB i ,j 1 Q i ,j 1 ] t RH
“平均斜率”，不等于“斜率的平
均值”，也不等于中点处的斜率
a ~j1/2a(uj1/2)1 2(ajaj1)
根据Langrage中值定理，[uL,uR]之间必有 一点uRoe, 该点处的斜率为平均斜率; 二次函数f(u)=u2中点处的斜率=平均斜率
Copyright by Li Xinliang
4
2. 方程组的情况
Copyright by Li Xinliang
12
§ 6.2 时间推进方法
1. 显格式
U f1 (U f2 ( )U f3 ( )U g 1 ) (U g 2 ) (U g 3 ) (U) t x y z x y z
U Q(U) t
推荐方法： Runge-Kutta法
可能出现导数不连续， 可能引起数值振荡
f ()
实际使用时 k 可用如下函数代替 ——所谓“熵修正”
f() (22当 )/2当
实际上是在特征值0点周围增加了耗散
Roe 格式的优点： 1） 保持守恒性的同时，严格保证了特征方向 2） 便于推广到高精度格式—— 特征投影分裂中使用Roe平均即可 （见
平均斜率 ~ A(UR,UL)
Uf(U) 0 t x
UA(U)U0, Af(U)
t
x
U
线性化，以平均增长 率代替瞬时增长率
f 1
x
j
( x
f
j1/2
fj1/2)
fˆ(Uj1,Uj) 经常记 fˆ(UR,为 UL)
[j,j+1]区间内
UA~ U0 t x
常系数方程的 Riemann解
~ A(UR,UL)
应当具有的性质
~
fR ( ) f U L ( ) A U (R , U U L )R ( U U L )
~ A(UR,UL)
连续，且
A ~(U U,)A(U)
~ A(UR,UL)
可通过相似变换对角化
1
1 ~
fj 1 /2 2 [fR )( fU ( L )U ] 2 A (R U ,U L )(R U U L )
UR
A ~(U R,UL)S1ΛS
UL
Copyright by Li Xinliang
5
针对Euler方程的具体构造
Uf(U)0 t x
U (,u ,E ) T ,f( U ) (u ,u 2 p ,u ( E p )T ) f(U)不是U的二次
齐函数
引入新变量： w1
1
W
w2
1/
Copyright by Li Xinliang
15
2) LU-SGS 方法
a) 将矩阵A,B分裂
A A A ,B B B
A (A λA *I)/2 , λ* Amak(A x ))(
为了简化计算，通常采用L-F分裂
B (B λB *I)/2 , λ* Bmak(B x))(
Qnt[AQ nBn Q ]tRHS
Qnt[AQ nBn Q ]tRHS
(2)
x
y
AA(U n),BB(nU ) 已知，（2）为线性方程。
未知量： QΔU
（1） (2) 是一个线性化过程
线性方程，离散求解
含义： 先用显格式计算，再 用隐格式计算修正量
离散后为大型带状方程组，求解计算量大
LU-SGS
原理： 将矩阵分解为上、下三角阵，避免矩阵求逆运算 隐式问题显式化
6） 计算空间导数
f(U) x
j
1x(fj1/2fj1/2)
7）时间推进，计算下一时间步的值。
与前文（第3，4讲）的形式相同， 仅需把式中的密度、压力、速度等换 成经过Roe平均的密度、压力、速度 即可
其中：
di(a 1,g2,3)
Copyright by Li Xinliang
10
di(a 1,g2,3)