CFD2020-第5讲-差分方法3
合集下载
差分方法
uu(1:100,1)=(x-0.5).^2;
figure
h=plot(x,uu(:,1),'linewidth',5);
set(h,'EraseMode','xor')
axis([0,1,0,0.25]);
fork=2:200
uu(2:99,2)=(1-2*c)*uu(2:99,1)+c*(uu(3:100,1)+uu(1:98,1))-b*dt/dx*(uu(3:100,1)-uu(2:99,1));
plot(u(1,:))
subplot(2,1,2)
plot(u(end,:))
差分方程所得的数值解的图形如图4所示，其中（a）是开始状态，(b)是最后状态。
(a) 初始状态
(b) 最后状态
【程序】
N=500;dx=0.01;dt=0.000001;
c=50*dt/dx/dx;
A=500;b=5;
x=linspace(0,1,100)';
方程设置是parobolic型，系数取为 。
解题的时间范围为 ，初始条件是 。
为了有足够的精度，将初始化的网格作了两次细分。而作图的选项为Contour和Animation。
作为对比，可以更改初始条件为 ,即 。
资料来源：数学物理方程与Matlab可视化.
(a) 整体图
(b) 上图：初始状态，(c)下图：最后状态
图3 解析解的图形
【程序】：
a2=50;b=5;
[x,t]=meshgrid(0:0.01:1,0:0.000001:0.0005);
Anfun=inline('2*(x-0.5).^2.*exp(5*x./2./50).*sin(n*pi*x)','x','n');
figure
h=plot(x,uu(:,1),'linewidth',5);
set(h,'EraseMode','xor')
axis([0,1,0,0.25]);
fork=2:200
uu(2:99,2)=(1-2*c)*uu(2:99,1)+c*(uu(3:100,1)+uu(1:98,1))-b*dt/dx*(uu(3:100,1)-uu(2:99,1));
plot(u(1,:))
subplot(2,1,2)
plot(u(end,:))
差分方程所得的数值解的图形如图4所示，其中（a）是开始状态，(b)是最后状态。
(a) 初始状态
(b) 最后状态
【程序】
N=500;dx=0.01;dt=0.000001;
c=50*dt/dx/dx;
A=500;b=5;
x=linspace(0,1,100)';
方程设置是parobolic型，系数取为 。
解题的时间范围为 ，初始条件是 。
为了有足够的精度，将初始化的网格作了两次细分。而作图的选项为Contour和Animation。
作为对比，可以更改初始条件为 ,即 。
资料来源：数学物理方程与Matlab可视化.
(a) 整体图
(b) 上图：初始状态，(c)下图：最后状态
图3 解析解的图形
【程序】：
a2=50;b=5;
[x,t]=meshgrid(0:0.01:1,0:0.000001:0.0005);
Anfun=inline('2*(x-0.5).^2.*exp(5*x./2./50).*sin(n*pi*x)','x','n');
第五章差分法和变分法解决平面问题
(f)
第五章
用差分法和变分法解平面问题
边界条件
再将式(f )对s 积分，从固定的基点A到边
界任一点B，得
B Φ Φ ( ) B ( ) A f x ds, A y y (g) B Φ Φ ( ) B ( ) A f y ds. A x x 式( f )、(g)分别是应力边界条件的微分、积 分形式。
第五章
用差分法和变分法解平面问题
边界条件
⑷ 由式(i)的第三式，可求出边界点的 ΦB
值; 由式(i)的前两式，可求出边界点
的
Φ ( )B x
、 Φ ) (
y
值，然后再求出边
B
界外一行虚结点的 Φ 值。
第五章
用差分法和变分法解平面问题
求解步骤
4.应力函数差分解的步骤 （1）在边界上选定基点A， 令
Φ Φ 然后计算边界上各结点的Φ 、x 、 y ；
Φ Φ ， ΦA ( ) A ( ) A 0 x y
（2）由边界结点的
Φ Φ 、 值，求出边界 x y
外一行虚结点的 Φ 值；
第五章
用差分法和变分法解平面问题
求解步骤
（3）对边界内所有结点列式(e)的方程，
联立求各结点的 Φ 值；
y
y
第五章
用差分法和变分法解平面问题
由于 ( T ) 2 T
y
10
T0 所以得 , 2h
这时，边界点2的 T2 是未知的，对2点 须列出式（d）的方程。此方程涉及到 T1 0 值，可将式（e）代入。
T1 0 T0
2h( q y ) 2
.
(e)
第五章
用差分法和变分法解平面问题
第五章 有限差分法 知识讲解课件
的 m=4，即此表对应差商的精度是四阶的。从这些表可以看出，一般地说，随着
差分阶数的增大和对应差商精度的提高，差分表达式所包含的项数将增多。
表 5-1
j
n0 1 2 34
1 -1
aj 1
2 1 -2 1
3 -1 3 -3 1
4 1 -4 6 -4 1
表 5-3 j
n0 1 2345 aj
1 -3 4 -1 2 2 -5 4 -1 3 -5 18 -24 14 -3 4 3 -14 26 -24 11 -2
依此类推，任何阶差分都可由其低一阶的差分再作一阶差分得到。例如 n 阶前差
分为
∆n y = ∆(∆n−1 y) = ∆[∆(∆n−2 y)]
⋯⋯ = ∆{∆⋯[∆(∆y)]} = ∆{∆⋯[∆( f (x + ∆x) − f (x)]}
n 阶的向后差分、中心差分的型式类似。
（5-6）
函数的差分与自变量的差分之比，即为函数对自变量的差商。如一阶向前差
二阶差商多取中心式，即
∆2 y ∆x 2
=
f (x + ∆x) − 2 f (x) + (∆x) 2
f (x − ∆x) 。
（5-9） （5-10） （后的二阶差商。 以上是一元函数的差分与差商。多元函数 f(x,y,…)的差分与差商也可以类推。
如一阶向前差商为
应地，上式中的 ∆y 、 ∆x 分别称为函数及自变量的差分， dy //#######为函数对 dx
自变量的差商。 在导数的定义中 ∆x 是以任意方式趋近于零的，因而 ∆x 是可正可负的。在差
分方法中， ∆x 总是取某一小的正数。这样一来，与微分对应的差分可以有 3 种
形式： 向前差分 向后差分 中心差分
CFD2020-第3讲-有限差分法(1)
1
(x j2 x j )2 a1 (x j1 x j )2 a2 02 a3 (x j1 x j )2 a4 0
(x j2 x j )3 a1 ( x j1 x j )3 a2 03 a3 (x j1 x j )3 a4 0
解出系数
a1 j , a2 j , a3 j , a4 j
3个点上信息计算
u x j
Step2: 写成待定系数形式
u x
j
a1u j2
a2u j1
a3u j
O(xn )
Copyright by Li Xinliang
5
u x
j
a1u j2
a2u j1
a3u j
O(xn )
Step3: 利用Taylor展开，确定系数
… j-2 j-1 j j+1 …
否则会降低精度
网格间距变化要缓慢，否则会带 来较大误差
Copyright by Li Xinliang
13
方法2） 在非等距网格上直接构造差分格式 （不易推广到高维）
原理： 直接进行Taylor展开，构造格式 格式系数是坐标（或网格间距）的函数
u x
j
a1u j2
a2u j1 a3u j
a4u j1 O(3 )
注： 系数随网格点(j)变化！
网格非光滑、间距剧烈变化不会降低精度； 随机网格
也可保证精度 ； 不易推广到高维
Copyright by Li Xinliang
14
2. 二维情况
物理空间
坐标变换 均匀的直角网格
x x( ,)
y
y(
,)
(x, y) (,)
计算空间
U t
弹性力学简明教程第四版第五章：有限差分发和变分法概论
§5－2应力函数的差分解
1、应力分量（不计体力）
一旦求得弹性体全部节点的 值后，就可按应力分量差分公式（对
节点0）算得弹性体各节点的应力。
0
•12 •8 •4 •5
h
x
x
0
2
y 2
0
1 h2
[(2
4 )
20 ]
• 11
•3
•0 •1
• 9
A
• 13
• 7
• 2
•6
y
0
2
x2
0
1 h2
[(1
3 )
20 ]
（5－9）
•10
B
h
• 14
xy
0
2
xy
0
1 4h2
[(5
7 ) (6
8 )]
y
图5－1
如果知道各结点的 值，就可以求得各结点的应力分量。
弹性力学简明教程
NORTHEASTERN UNIVERSITY
§5－2应力函数的差分解
2、差分方程（相容方程） 双调和方程
•12
x 设：f f x, y 为弹性体的某一连续函数
h
•8
•4
• 5
•11 •3 •0 •1 •9 A • 13
在平行与 x 轴的一根网线上函数只随 x
• 7
• 2
•6
坐标的变化而变化。
•10
B
h
• 14
y
图5－1
在节点0 的近处将函数 f 展成泰勒级数
f
f0
f x
0
x
x0
1 2!
2 f x2
y x （3）由于 f 是 或 的二次函数，所以基本差分公式（5－1）
1、应力分量（不计体力）
一旦求得弹性体全部节点的 值后，就可按应力分量差分公式（对
节点0）算得弹性体各节点的应力。
0
•12 •8 •4 •5
h
x
x
0
2
y 2
0
1 h2
[(2
4 )
20 ]
• 11
•3
•0 •1
• 9
A
• 13
• 7
• 2
•6
y
0
2
x2
0
1 h2
[(1
3 )
20 ]
（5－9）
•10
B
h
• 14
xy
0
2
xy
0
1 4h2
[(5
7 ) (6
8 )]
y
图5－1
如果知道各结点的 值，就可以求得各结点的应力分量。
弹性力学简明教程
NORTHEASTERN UNIVERSITY
§5－2应力函数的差分解
2、差分方程（相容方程） 双调和方程
•12
x 设：f f x, y 为弹性体的某一连续函数
h
•8
•4
• 5
•11 •3 •0 •1 •9 A • 13
在平行与 x 轴的一根网线上函数只随 x
• 7
• 2
•6
坐标的变化而变化。
•10
B
h
• 14
y
图5－1
在节点0 的近处将函数 f 展成泰勒级数
f
f0
f x
0
x
x0
1 2!
2 f x2
y x （3）由于 f 是 或 的二次函数，所以基本差分公式（5－1）
差分法
第三章 有限差分法函数()f x ，x 为定义在区间[]a b ,上的连续变 量。
将区间[]a b ,等分成n 份，令()h b a n =-称为 步长，x 在这些离散点处的取值为x a ih i =+ ()i n =01,,,称为节点。
函数()f x 在这些节点处的差值()()()()()()f x h f x f x f x h f x h f x h i i i i i i +---+--⎧⎨⎪⎩⎪(5-1)分别称为一阶向前、向后和中心差分，可以用它 们作为函数()f x 在x i 处的微分近似值。
这些差分 与相应x 区间的比值()()[]()()[]()()[]1112h f x h f x h f x f x h h f x h f x h i i i i i i +---+--⎧⎨⎪⎪⎪⎩⎪⎪⎪ (5-2) 分别称为一阶向前、向后和中心差商，可以用它 们作为函数()f x 在x i 处的导数近似值。
完全类似地可以定义高阶差商，例如常用的二阶中心差商()()()[]122hf x h f x f x h i i i +-+- (5-3)可以作为函数()f x 在x i 处的二阶导数近似值。
§3.1 常微分方程初值问题的差分解法考虑电学中的一个问题：如图5-1。
研究 电容器上的电荷随时间的变化规律。
图5-1 RC 放电回路这个问题对应的微分方程及其定解条件为：d d Q tQ RC QQ t =-=⎧⎨⎪⎩⎪=00(5-4) 这是一阶微分方程的初值问题，它的解析解为 Q Q e t RC =-0 (5-5)一、欧拉(Euler )折线法求解下列普遍形式的一阶微分方程的初值 问题：()[]()'=∈=⎧⎨⎪⎩⎪y f x y x a b y a y ,,0(5-6) 首先，将区间[]a b ,等分n 份，取值a x x xb n =<<<=01 ，步长h x x i i =-+1。
【计算流体力学】第5讲-差分方法3
通量差分分Байду номын сангаас （FDS）： 耗散低、分辨率高
Step 1: 运用差分格式，计算
U ,U L j 1/ 2
R j 1/ 2
Step 2: 运用Riemann解， 计算
F j 1/ 2
F
(U
L j 1/
2
,U
) R
j 1/ 2
Step 3: F Fj1/2 Fj1/2
x
x
U ,U L j 1/ 2
R j 1/ 2
f
j 1/
2
x
x
f f + f =+
x x x
13
3. 特征重构方法
常系数方程组：
U t
A U x
0
U t
S1ΛS U x
0
V t
Λ V x
0
vk t
k
vk x
0
变系数情况—— 局部冻结系数
完全 解耦
U f(U) 0 U A U 0
t x
t x
在基架点上系数 A j 不变
U t
u j u j1 u j1 u j
二阶精度区
TVD区
二阶精度TVD区（二 者交集）
1
通量分裂技术： 模型方程 NS/ Euler 方程
Step 1 针对模型方程构造差分格式
u a u 0 t x
u uˆ j1/2 uˆ j1/2
x
x
a0
uˆ j1/2 =......
格式1
a0
uˆ j1/2 =......
6
➢ Steger-Warming 具体步骤 （以一维为例）
u a u 0
t x
CFD2011-第5讲-差分方法3
单调
保单调
10
3. TVD格式的理论基础 格式的理论基础—— Harten定理 格式的理论基础 定理
Harten定理： 定理： 定理
如果差分格式可写成如下形式： 如果差分格式可写成如下形式： 保证“系数非负” 保证“系数非负”
u n +1 j = un j + C n+ 1 j 2 (u n+1 j − un) j
3. 流通矢量分裂
f = f+ +f−
f ± = A ± U = S −1 Λ ± SU
Copyright by Li Xinliang
3
§ 5.1
非物理振荡及TVD格式 格式 非物理振荡及
1. 数值解中的非物理振荡
1） 非物理振荡的原因分析 ） 理论1： 理论 ：色散误差导致各波 传播速度不同 （第4讲） 讲 理论2： 理论 ：粘性耗散不足 理论3： 理论3： 格式不能保单调
格式： 格式： 如果满足
则称其为单调格式。
k n j+k
u n +1 = j
∑a u
k
ak ≥ 0
单调格式 j=1
n +1 j
保单调格式： 保单调格式：
是单调的，如果n+1时刻的解 {u } 设n时刻 {u } 是单调的，如果 时刻 时刻的解 仍保证单调，则称该格式为保单调格式 保单调格式。 仍保证单调，则称该格式为保单调格式。
j-1
j
过于苛刻的条件 单方向网格点数10 三维10 单方向网格点数 6， 三维 18
网格Reynolds数足够小时，物理粘性发挥作用，抑制振荡 数足够小时，物理粘性发挥作用， 网格 数足够小时
单纯靠物理粘性抑制振荡，网格间距必须足够小， 单纯靠物理粘性抑制振荡，网格间距必须足够小，通常难以实现
西交大《计算流体动力学》课件CFD第三章
n
n
n
u i u i 1 u o x x x i
n n
n
计算流体动力学课程
西安 2005年3月
XJTU
2.构造差分格式的方法
一阶导数的差商表达式 名称 一阶前向差商 表达式
f i 1 f i x
f i f i 1 x
f x i
t 0
边界条件
t 0
计算流体动力学课程
西安 2005年3月
XJTU
1. 有限差分方法
网格划分及方程离散化
t (j,n) n x
x j j x
j 0 ,1, 2 , J
t n nt
n 0 ,1, 2 ,
u
n j
u j u j x , n t
n n n n
ui
n 1
ui
n 1
t
x
u
n i 1
u i 1
n
西安 2005年3月
XJTU
3.模型方程的差分格式
u t u x 0
格式名称
二阶精度的两步 MacCormack显示格式
ui u
n 1
格式表达式
u
n i
t
x
u
n i 1
二阶精度的时间－中 心隐式格式
ui
n 1
ui
n
计算流体动力学课程
西安 2005年3月
XJTU
3.模型方程的差分格式
u t b u
2
x
2
格式名称 时间一阶精度、空间二 阶精度的显示格式 二阶精度的三层 Richardson显示格式 二阶精度的Crank- u u t Nicolson隐式格式
【计算流体力学】第5讲-差分方法3
ui1, ui1,
ui ui
, ,
ui1 ui1
, ,
...) ...)
if if
a0 a0
f
j 1/ 2
=f
j1（ /2 ...,f
,
j 1
f
j
,
f
j1,...)
具体步骤 （逐点分裂）：
通量分裂
fj
j=1,2,3……N
f
j
,
f
j
差分格式
f ,f
j +1/2 j +1/2
f +
=
f
j 1/ 2
利用了性质 f(U) A U (AU)
x
x x
一般情况下：
变系数，A 不能与导数交换
f AU A U
x x
x
f B U
x
x
B A
分裂后 A 失去了A的性质（可以 像常数一样与求导交换）
实质： 没有做到解耦；
只是把原变量重新组合，组合后波的传播方向的保证 f+ 向正 向传播，f-向负向传播
f
j 1/
2
x
x
f f + f =+
x x x
13
3. 特征重构方法
常系数方程组：
U t
A U x
0
U t
S1ΛS U x
0
V t
Λ V x
0
vk t
k
vk x
0
变系数情况—— 局部冻结系数
完全 解耦
U f(U) 0 U A U 0
t x
t x
在基架点上系数 A j 不变
U t
f 0
弹性力学-05(差分法与变分法)
对于弹性体边界内的每一节点，都可 建立上述方程。 但对紧靠边界内一行节点， 建立其差分方程时，还包括边界上各点处的 值和边界外一行的结点处的 值。 8 11 3
（5-10）
—— 应力函数差分方程 x 12 4 0 5 1 9 A 13
弹性体边界外一行的节点，称为虚结点。 如：节点13、14等。
（c ）
y
h
将其代入式（b），有：
2 2 f h f f 3 f 0 h x 2 x 2 0
0 2 2 f h f f1 f 0 h 2 x 0 2 x 0
任一点 0 处应力分量的差分格式：
2 1 x 0 y 2 h 2 ( 2 4 ) 2 0 0 2 1 y 0 2 2 (1 3 ) 2 0 x 0 h
x 12 8 11 3 7 4 0 2 10 y h 5 1 6 9
在弹性体内每一点均可建立上述方程，即：
0 x 4 2 x 2y 2 y 4 0 0 0 0
4
4
y h
优点： 收敛性好、程序设计简单、 非线性适应好。 代表性软件：FLAC
f f f 1 3 2h x 0
y h
x h
3 0 1
缺点：当边界几何形状复杂时，解的精度受到限制。 （2）等效积分法 控制微分方程 边值条件 建立等效的 积分方程 假设未知函数 整个区域内
定值条件 精确解 （均质、边界条件简单）
近似解 （1）有限差分法 （数值解） （2）等效积分法（包括变分法） （3）有限单元法 （4）边界单元法 …… f1 f 3 （1）有限差分法（FDM） f 代替 2h x 0 要点：差分 微分； x h 3 0 1
（5-10）
—— 应力函数差分方程 x 12 4 0 5 1 9 A 13
弹性体边界外一行的节点，称为虚结点。 如：节点13、14等。
（c ）
y
h
将其代入式（b），有：
2 2 f h f f 3 f 0 h x 2 x 2 0
0 2 2 f h f f1 f 0 h 2 x 0 2 x 0
任一点 0 处应力分量的差分格式：
2 1 x 0 y 2 h 2 ( 2 4 ) 2 0 0 2 1 y 0 2 2 (1 3 ) 2 0 x 0 h
x 12 8 11 3 7 4 0 2 10 y h 5 1 6 9
在弹性体内每一点均可建立上述方程，即：
0 x 4 2 x 2y 2 y 4 0 0 0 0
4
4
y h
优点： 收敛性好、程序设计简单、 非线性适应好。 代表性软件：FLAC
f f f 1 3 2h x 0
y h
x h
3 0 1
缺点：当边界几何形状复杂时，解的精度受到限制。 （2）等效积分法 控制微分方程 边值条件 建立等效的 积分方程 假设未知函数 整个区域内
定值条件 精确解 （均质、边界条件简单）
近似解 （1）有限差分法 （数值解） （2）等效积分法（包括变分法） （3）有限单元法 （4）边界单元法 …… f1 f 3 （1）有限差分法（FDM） f 代替 2h x 0 要点：差分 微分； x h 3 0 1
计算流体力学CFD基础课程PPT CFD2013-第5讲-差分方法3
计算流体力学讲义2013
第五讲 差分方法（3）
知识点： 激波捕捉格式—— TVD、WENO、MUSCL、NND
1
知识回顾
1. 差分格式的分辨率
修正 波数
u(x j ) eikxj
~
xu j
k eikx j x
2. 群速度控制(GVC)
u
1
0.5
0
-0.5
-1
0
1
2
3x
4
5
6
1
0.5
0
-0.5
-1
4
数值实验
u t
u x
1 Re
2u x2
1 j
u
n j
t
u
n j 1
u
n j 1
2x
1 Re x2
(u
n j 1
2u
n j
u
n j 1
)
计算域[0,1], 网格点201 （x=0.005 ） 时间步长x=0.0005
二阶中心差分
相同
Re=200 x=0.005
(1)
k
j-2 j-1 j j+1 j+2
n+1
n
u n1 j
a un 2 j2
a un 1 j 1
......
a2
u
n j
2
某点的值是上一时刻周围几个点上值的线性组合
物理上要求系数 ak 均非负
含义： 某处浓度的增加对下一时 刻周围浓度的影响为正。
差分方程单调性（无振荡）条件： 差分方程 （1）中的系数非负
特征线未相交—— 总变差不变
特征线相交—— 总变差减小
结论： 单个双曲型方程，总变差不增 （Total Variation Diminishing： TVD）
第五讲 差分方法（3）
知识点： 激波捕捉格式—— TVD、WENO、MUSCL、NND
1
知识回顾
1. 差分格式的分辨率
修正 波数
u(x j ) eikxj
~
xu j
k eikx j x
2. 群速度控制(GVC)
u
1
0.5
0
-0.5
-1
0
1
2
3x
4
5
6
1
0.5
0
-0.5
-1
4
数值实验
u t
u x
1 Re
2u x2
1 j
u
n j
t
u
n j 1
u
n j 1
2x
1 Re x2
(u
n j 1
2u
n j
u
n j 1
)
计算域[0,1], 网格点201 （x=0.005 ） 时间步长x=0.0005
二阶中心差分
相同
Re=200 x=0.005
(1)
k
j-2 j-1 j j+1 j+2
n+1
n
u n1 j
a un 2 j2
a un 1 j 1
......
a2
u
n j
2
某点的值是上一时刻周围几个点上值的线性组合
物理上要求系数 ak 均非负
含义： 某处浓度的增加对下一时 刻周围浓度的影响为正。
差分方程单调性（无振荡）条件： 差分方程 （1）中的系数非负
特征线未相交—— 总变差不变
特征线相交—— 总变差减小
结论： 单个双曲型方程，总变差不增 （Total Variation Diminishing： TVD）
中科院计算流体力学最新讲义CFD2020第5讲差分方法3
•双曲型守恒方程
•特点： 沿特征线
， u不变
•特征线未相交— —总变差不变
•特征线相交—— 总变差减小
•结论： 单个双曲型方程，总变差不增 •（Total Variation Diminishing： TVD）
•j= 1
•j=N 单调函数
•振荡函数
by Li Xinliang
•2 概念： 单调格式、保单调格式与TVD格式 •j=
•二阶迎风
•二阶中心
by Li Xinliang
•新格式 ：
•根据Harten定理，可知 •时，可满足TVD性质
•(2) 精度条件
•显然
格式为2 阶中心
•可验证： 格式为2阶迎风
•二者组合仍为二阶
•二阶精度区
•TVD区
•二阶精度TVD区（ 二者交集）
by Li Xinliang
•限制器（limiter)
•2阶中心 的修正量
•2阶迎风 的修正量
•精度高，但有些情况下预 测结果“不靠谱”
•作为“标杆”检 验高阶修正量是 否可用
•趋势相反时，不可用； •相差超过2倍时，不可用
by Li Xinliang
•2 以 L-W格式为基础改造的格式
•历史上，TVD格式是在Roe、L-W、B-M （或其组合）基础上改进 •80年代初、这些格式是主流
•很难计算对粘性敏感的问题
解
•改进措施：
•
A: 局部施加人工粘性
•
B: 高阶人工粘性
•分离流—— 对粘性敏感
•Von Neumann
•MacCormack
•转捩——对粘性敏
by Li Xinliang
•4) 数值振荡的定量描述—— 总变差
•特点： 沿特征线
， u不变
•特征线未相交— —总变差不变
•特征线相交—— 总变差减小
•结论： 单个双曲型方程，总变差不增 •（Total Variation Diminishing： TVD）
•j= 1
•j=N 单调函数
•振荡函数
by Li Xinliang
•2 概念： 单调格式、保单调格式与TVD格式 •j=
•二阶迎风
•二阶中心
by Li Xinliang
•新格式 ：
•根据Harten定理，可知 •时，可满足TVD性质
•(2) 精度条件
•显然
格式为2 阶中心
•可验证： 格式为2阶迎风
•二者组合仍为二阶
•二阶精度区
•TVD区
•二阶精度TVD区（ 二者交集）
by Li Xinliang
•限制器（limiter)
•2阶中心 的修正量
•2阶迎风 的修正量
•精度高，但有些情况下预 测结果“不靠谱”
•作为“标杆”检 验高阶修正量是 否可用
•趋势相反时，不可用； •相差超过2倍时，不可用
by Li Xinliang
•2 以 L-W格式为基础改造的格式
•历史上，TVD格式是在Roe、L-W、B-M （或其组合）基础上改进 •80年代初、这些格式是主流
•很难计算对粘性敏感的问题
解
•改进措施：
•
A: 局部施加人工粘性
•
B: 高阶人工粘性
•分离流—— 对粘性敏感
•Von Neumann
•MacCormack
•转捩——对粘性敏
by Li Xinliang
•4) 数值振荡的定量描述—— 总变差
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
计算流体力学讲义2020
第五讲 差分方法（3）
李新亮 知识点：
通量技术简介——Steger-Warming, Roe 常用的隐式处理方法——LU-SGS
1 Copyright by Li Xinliang
知识回顾： 单调、保单调和TVD
概念： 网格Reynolds数 Re x Rex 单调格式、保单调格式及TVD格式
2
1
1 (3 )u c2 3 u2 1 2
0
1
u
u E
验证方法2： 利用齐函数性质
f (U) f(U) (f(αU) ) U f(U) 1
(αU)
f(U) (f(U)) U AU U
Copyright by Li Xinliang
5
2. 流通矢量分裂(FVS)
f1 [u, u2 p, uv,u(E p)]T
f 2 [v, uv, v2 p,v(E p)]T
令： f f1 f2
具体使用步骤， 以计算 f1 为例 x
则：
~f (~)
2
~0
2
V2
~3 2
[~~u001u2v~0v~~1233~]uv311 ~2~4~~444[vuu2222
f x
f x
f x
1 [ (f 2 x
*U)
x
(f
*U)]
1 2
[
x
(f
*U)
(f x
*U)]
1 2
[(
x
x
)f
]
* 2
[(
x
x
)U]
x0f
*x
xx U
U f(U) 0 t x
U t
f(U) x
*x
2U x2
Copyright by Li Xinliang
10
1 u, 2 u c, 3 u c k , k (k 1,2,3)
4） 带入（1）式得到 f , f
5） 利用不同的迎风格式，分别计算
f f ,
x x
（后差，前差）
6）计算 f f f x x x
k
k
(2k 2 )1/ 2 2
（1） ~f (~λ)
2
(
21()~1u12)2~(1u~22(1u~)2~(1uc)2~c2)~23~~3(3u(uc)c2)
4
§ 5.1 流通矢量分裂： Steger-Warming 、L-F
1. Jacobian 系数矩阵及其性质
U f(U) vis t x
U (, u, E)T ,
f(U)
f1 f2 f3
u
u
2
u(E
p
p)
u2
( 1)u3
u2u3 u1
3 u22 2 u1
v22
]
W
W
(3 )(~3 ~4 )c2 2( 1)
u1 u ck~1 u2 u ck~1
~ k1
/
~ k2
/
~0 2( 1)~1
v1 v ck~2 v2 v ck~2
2 2
对于曲线坐标系
Uˆ fˆ1 fˆ2 0 t
1) 令 1, 0 2) 计算特征值
1 2 u,3 u c,4 u c
w
w (3 )(~2 ~3 )c2 2( 1)
f ~f (λ ), f ~f (λ )
f f f
7） 时间推进 U f 0 t x
Copyright by Li Xinliang
8
➢二维问题的steger-Warming 分裂
U f1 f2 0 t x y
U [, u, v, E]T
Steger-Warming 分裂
f ~f (λ )
特点： 不必进行矩阵运算，计算量小
Copyright by Li Xinliang
7
➢ Steger-Warming 具体步骤 （以一维为例）
u a u 0
t x
f
a0
x
已知 U (, u, E)T
f
x
a0
1） 计算 ,u, p
2） 计算 3） 计算
2
2
如果差分格式无耗散（例如 都用中心差分），则通量分 裂不带来耗散。
x f x f x f
U t
0 x
f
a
U 0
xx
耗散
x f
0 x
f
0 x
f
0 x
(
f
f
)
0 x
f
向上平移
分裂后的流场越偏离原先流场，则 总体耗散越大
=
+
向下平移
如使用低精差分度格式， 则对分裂形式敏感 （推荐使用特征分裂） 如使用高精度格式（低耗散），则对分裂形式不敏感 （可使用逐点分裂）
Copyright by Li Xinliang
13
➢ FVS差分方法一般流程
Step 1. 构造差分格式 u a u 0 t x
u ui1/2 ui1/2
x
x
Step 2. 推广到方程组
U f(U) 0 t x
f f + f =+
x x x
f +
=
f j 1/ 2
f
j 1/
2
x
x
差分格式
3) 分裂特征值，计算 k ,k
(k 1,2,3,4)
k
k
(2k 2 )1/ 2 2
4) 带入左式，计算正、负流通矢
量
f1 ~f ( )
5) 计算
f1 f1 f1 x x x
fˆ1 J 1xf1 J 1yf2
仅需令
J
1 x ,
J
1 y
三维问题同样处理
计算 f2 设置 0, 1 ，并注意
… j-2 j-1
j
j+1 …
1） 逐点分裂 f k
通量分裂
fk , fk k=1,2…..N
j+1/2
2）针对j点：计算特征矩阵
U j1/2
S j1/2
3）（将网格基架点）变换到特征空间
S1 j 1/ 2
简单平均即可（Roe平均效 果更好些）
U j1/2 =（U j U j1) / 2
A j1/2 =Sj11/2ΛS j1/2
完全 解耦
U f(U) 0 U A U 0
t x
t x
在基架点上系数 A j 不变
U t
j
Aj
U x
j
0
… j-2 j-1 j j+1 …
U
计算：A j x j 在差分基架点上Aj 不变， 可按常矩阵处理
U A j x j
A j Sj 1ΛjSj
局部冻结系数
A
j
U x
j
S
1 j
ΛS
j
U x
y 1 2 v,3 v c,4 v c
二维、三维具体 公式见傅德薰等《计算空气动力学》 4.7节 （158-162）
书中公式有一定的排版错误，使用C前opy务righ必t by重Li X新inli仔ang细推导！
9
B: Lax-Friedrichs (L-F)分裂
A (A *I) / 2, A (A *I) / 2
f(U) AU
特点： A 正特征值 A 负特征值
f AU (f *U) / 2
* 足够大
缺点：耗散偏大
例如，可取 * u c
=
+
局部L-F分裂，每个点上计算 * u c 常数
数学性质（光滑性） 最好，但耗散偏大
全局L-F分裂，全局（一维）上计算 * max ( u c) x
与迎风格式结合，等价于人工粘性
U f(U) 0 t x
方法1： 流通矢量分裂(FVS)
f(U) f (U) f (U)
f f f x x x
Copyright by Li Xinliang
3
方法2： 通量差分分裂 （FDS）
利用 Riemann解
u a u 0 t x
特点： 对流通矢量 f 的导数进行分裂
f(U) AU
f f f
f AU S1ΛSU
向量分裂 特征值标量分裂
A: Steger-Warming 分裂
Λ diag{k }
Λ Λ Λ Λ diag{k }
k
2
, k
2
任意函数分解为：
非负函数+非正函数
=
优点：耗散小 缺点：导数间断
+
Copyright by Li Xinliang
Harten定理： 正系数原则
TVD
单调
保单调
TVD格式= 1阶迎风+ j *（修正项）
u j1/ 2 u j j (r)(u j1 u j ) / 2
j j j(rj ),
rj
u j u j1 u j1 u j
二阶精度区
TVD区
Copyright by Li Xinliang
二阶精度TVD区（二 者交集）
6
k
2
, k
2
f AU S1ΛSU
优点：耗散小 缺点：导数间断
=
+
~f (~λ)
2(21()~1u Nhomakorabea2)2~(1u~22(1u~)2~(1uc)2~c2)~23~~3(3u(uc)c2)
w
k
k
k 2
改进版
k
k
(2k )2 1/ 2 2
w (3 )(~2 ~3 )c2 2( 1)
5）变换回物理空间
f j1/2
第五讲 差分方法（3）
李新亮 知识点：
通量技术简介——Steger-Warming, Roe 常用的隐式处理方法——LU-SGS
1 Copyright by Li Xinliang
知识回顾： 单调、保单调和TVD
概念： 网格Reynolds数 Re x Rex 单调格式、保单调格式及TVD格式
2
1
1 (3 )u c2 3 u2 1 2
0
1
u
u E
验证方法2： 利用齐函数性质
f (U) f(U) (f(αU) ) U f(U) 1
(αU)
f(U) (f(U)) U AU U
Copyright by Li Xinliang
5
2. 流通矢量分裂(FVS)
f1 [u, u2 p, uv,u(E p)]T
f 2 [v, uv, v2 p,v(E p)]T
令： f f1 f2
具体使用步骤， 以计算 f1 为例 x
则：
~f (~)
2
~0
2
V2
~3 2
[~~u001u2v~0v~~1233~]uv311 ~2~4~~444[vuu2222
f x
f x
f x
1 [ (f 2 x
*U)
x
(f
*U)]
1 2
[
x
(f
*U)
(f x
*U)]
1 2
[(
x
x
)f
]
* 2
[(
x
x
)U]
x0f
*x
xx U
U f(U) 0 t x
U t
f(U) x
*x
2U x2
Copyright by Li Xinliang
10
1 u, 2 u c, 3 u c k , k (k 1,2,3)
4） 带入（1）式得到 f , f
5） 利用不同的迎风格式，分别计算
f f ,
x x
（后差，前差）
6）计算 f f f x x x
k
k
(2k 2 )1/ 2 2
（1） ~f (~λ)
2
(
21()~1u12)2~(1u~22(1u~)2~(1uc)2~c2)~23~~3(3u(uc)c2)
4
§ 5.1 流通矢量分裂： Steger-Warming 、L-F
1. Jacobian 系数矩阵及其性质
U f(U) vis t x
U (, u, E)T ,
f(U)
f1 f2 f3
u
u
2
u(E
p
p)
u2
( 1)u3
u2u3 u1
3 u22 2 u1
v22
]
W
W
(3 )(~3 ~4 )c2 2( 1)
u1 u ck~1 u2 u ck~1
~ k1
/
~ k2
/
~0 2( 1)~1
v1 v ck~2 v2 v ck~2
2 2
对于曲线坐标系
Uˆ fˆ1 fˆ2 0 t
1) 令 1, 0 2) 计算特征值
1 2 u,3 u c,4 u c
w
w (3 )(~2 ~3 )c2 2( 1)
f ~f (λ ), f ~f (λ )
f f f
7） 时间推进 U f 0 t x
Copyright by Li Xinliang
8
➢二维问题的steger-Warming 分裂
U f1 f2 0 t x y
U [, u, v, E]T
Steger-Warming 分裂
f ~f (λ )
特点： 不必进行矩阵运算，计算量小
Copyright by Li Xinliang
7
➢ Steger-Warming 具体步骤 （以一维为例）
u a u 0
t x
f
a0
x
已知 U (, u, E)T
f
x
a0
1） 计算 ,u, p
2） 计算 3） 计算
2
2
如果差分格式无耗散（例如 都用中心差分），则通量分 裂不带来耗散。
x f x f x f
U t
0 x
f
a
U 0
xx
耗散
x f
0 x
f
0 x
f
0 x
(
f
f
)
0 x
f
向上平移
分裂后的流场越偏离原先流场，则 总体耗散越大
=
+
向下平移
如使用低精差分度格式， 则对分裂形式敏感 （推荐使用特征分裂） 如使用高精度格式（低耗散），则对分裂形式不敏感 （可使用逐点分裂）
Copyright by Li Xinliang
13
➢ FVS差分方法一般流程
Step 1. 构造差分格式 u a u 0 t x
u ui1/2 ui1/2
x
x
Step 2. 推广到方程组
U f(U) 0 t x
f f + f =+
x x x
f +
=
f j 1/ 2
f
j 1/
2
x
x
差分格式
3) 分裂特征值，计算 k ,k
(k 1,2,3,4)
k
k
(2k 2 )1/ 2 2
4) 带入左式，计算正、负流通矢
量
f1 ~f ( )
5) 计算
f1 f1 f1 x x x
fˆ1 J 1xf1 J 1yf2
仅需令
J
1 x ,
J
1 y
三维问题同样处理
计算 f2 设置 0, 1 ，并注意
… j-2 j-1
j
j+1 …
1） 逐点分裂 f k
通量分裂
fk , fk k=1,2…..N
j+1/2
2）针对j点：计算特征矩阵
U j1/2
S j1/2
3）（将网格基架点）变换到特征空间
S1 j 1/ 2
简单平均即可（Roe平均效 果更好些）
U j1/2 =（U j U j1) / 2
A j1/2 =Sj11/2ΛS j1/2
完全 解耦
U f(U) 0 U A U 0
t x
t x
在基架点上系数 A j 不变
U t
j
Aj
U x
j
0
… j-2 j-1 j j+1 …
U
计算：A j x j 在差分基架点上Aj 不变， 可按常矩阵处理
U A j x j
A j Sj 1ΛjSj
局部冻结系数
A
j
U x
j
S
1 j
ΛS
j
U x
y 1 2 v,3 v c,4 v c
二维、三维具体 公式见傅德薰等《计算空气动力学》 4.7节 （158-162）
书中公式有一定的排版错误，使用C前opy务righ必t by重Li X新inli仔ang细推导！
9
B: Lax-Friedrichs (L-F)分裂
A (A *I) / 2, A (A *I) / 2
f(U) AU
特点： A 正特征值 A 负特征值
f AU (f *U) / 2
* 足够大
缺点：耗散偏大
例如，可取 * u c
=
+
局部L-F分裂，每个点上计算 * u c 常数
数学性质（光滑性） 最好，但耗散偏大
全局L-F分裂，全局（一维）上计算 * max ( u c) x
与迎风格式结合，等价于人工粘性
U f(U) 0 t x
方法1： 流通矢量分裂(FVS)
f(U) f (U) f (U)
f f f x x x
Copyright by Li Xinliang
3
方法2： 通量差分分裂 （FDS）
利用 Riemann解
u a u 0 t x
特点： 对流通矢量 f 的导数进行分裂
f(U) AU
f f f
f AU S1ΛSU
向量分裂 特征值标量分裂
A: Steger-Warming 分裂
Λ diag{k }
Λ Λ Λ Λ diag{k }
k
2
, k
2
任意函数分解为：
非负函数+非正函数
=
优点：耗散小 缺点：导数间断
+
Copyright by Li Xinliang
Harten定理： 正系数原则
TVD
单调
保单调
TVD格式= 1阶迎风+ j *（修正项）
u j1/ 2 u j j (r)(u j1 u j ) / 2
j j j(rj ),
rj
u j u j1 u j1 u j
二阶精度区
TVD区
Copyright by Li Xinliang
二阶精度TVD区（二 者交集）
6
k
2
, k
2
f AU S1ΛSU
优点：耗散小 缺点：导数间断
=
+
~f (~λ)
2(21()~1u Nhomakorabea2)2~(1u~22(1u~)2~(1uc)2~c2)~23~~3(3u(uc)c2)
w
k
k
k 2
改进版
k
k
(2k )2 1/ 2 2
w (3 )(~2 ~3 )c2 2( 1)
5）变换回物理空间
f j1/2