计算天文数值计算lec15高级算法之CIP法lec15

19. 高级算法之CIP法
Yabe (矢部) et al. 1991, Comput. Phys. Commun., 66, 219 2001, JCP, 169, 556 发展了CIP (Cubic Interpolated Profile) 高精度差分法，举例如下：
f f c 0 t x

计算步骤
1. Euler步
忽略输运项，
上述方程组变为
U ( p q) t r V ( p q) t z e p ( rU ) 1 ( qrU ) q V qV q U p V t r r r r r z z z
fin 1 fin c
拟合系数满足此函数及其导数值在 i-1, i 格点处连续的条件:
n F ( xi1 ) f i 1
n ˆn f ˆn f fi n fi i i 1 1 ai 2 x 2 x 3
ˆn x F ( xi 1 ) f i 1
令 w(r ri ) 表示为the interpolation kernel centered around ri 则 r 处的局地密度为（假设总共有N个质点）

缺点
对下列情况效果不是非常好：
非常不均匀的粒子分布
correlated systems complex geometries
21. 高级算法之SPH法
PIC方法同时使用Eulerian和Lagrangian elements，其 欧拉网格的平均密度等于网格中的质点数目。
能否不需网格而得到局地密度？
对能量方程采用交错（拉链）差分
( fU )
1 i , j 2

1 ( f i , jU i 1, j f i 1, jU i , j ) 2
若旁边网格是空网格，则U及V均取本网格中心值， 若界面是固壁，则法向速度取为0。
2. 质点步 PIC方法的特点就是利用质点的运动来计算通过网格 边界的输运量 。 （1）确定质点速度
u u t x x u u 1 p u t x x e e p u u t x x
计算实例：激波管问题
1 p e 1
为使激波解稳定，可在压力中加人为粘性项
1 n 2 n ( c u i u ) i s qin 2 0
令总共有K个质点，每个质点具有其
n n r v m 质量 ，坐标 k ，速度 k
以二维为例，则网格(i, j)的性质为其中的各质点的总和：

n i,j
m K n r k (i , j ) xy k 1
n i,j


ein, j由in , j 及状态方程求出
n n xk yk r (i, j ) int( ) i int( ) x y
ˆ ˆ f f u 0 t x
1. 非对流项
f G t
ˆ f ˆ u ˆf G t x
采用时间前差，空间中心差分
fi* fi n Gi t
ˆ* f ˆ n Gi 1 Gi 1 t f ˆ n ui 1 ui 1 t f i i i 2 x 2 x * n * n ( f f ) ( f f n i 1 i 1 i 1 ) ˆ i 1 ˆ n ui 1 ui 1 t f f i i 2 x 2 x
非线性方程
即
f uf g t x
f f u u g f G t x x
ˆ ˆ f f u ˆ ˆ u G f t x x
上两式求解分两步实现
1. 非对流项
f G t
2. 对流项
f f u 0 t x
ˆ f ˆ u ˆf G t x
p e p e p e
1 2
常可用交错网格
u
x
1 i 2
u
x
i
xi 1
xi
xi 1
2. 对流项
f f u 0 t x
ˆ ˆ f f u 0 t x
采用CIP方法
ˆ* f * fi n1 ai 3 bi 2 f i i
d n 1 ˆ ˆ* fi F ( xi ut ) 3ai 2 2bi f i dx
要点：各物理量（除质量外）均定义在网格中心，混合网格 的质量和能量分别按物质给出。
初始时质点分布： 以2维均匀网格为例，每个网格内质点数目反映局地密度，且至少要 有100个质点；应使计算末每个网格至少要有3－4个质点（也许应更 多）。经验表明，一个网格内质点的分布不宜十分规则，采用菱形分 布或不规则分布利于减少质点流过网格边界引起的非物理计算波动。
n n ˆn f ˆn fi 2f 1 fi i i 1 bi 3 x 2 x
若u<0 则
i-1 i+1， x - x
n 因此，CIP法要求同时求解 fˆi
为此，对原方程求导
fˆ fˆ c 0 t x
ˆ n 也满足 f ˆ f (xi , t t) fˆ (xi ct, t) 同样， i
~ ~ U 、 V 在Euler步中得出的速度 i , j i , j是网格中心的速度，而质
点一般不在中心，每个质点的速度可取为与其相邻的四 个网格速度的加权平均值（第k个质点）:
~ n 1 4 ~ ~n 1 4 ~ U k aU Vk aV
1
1
设四个网格的公共交点是 ( r* , z* )
rkn r* r ，z
n zk z*
令r
1 1 a1 ( z )( r ) 2 2 1 1 a3 ( z )( r ) 2 2
z
3
1
1 1 a2 ( z )( r ) 2 2 1 1 a4 ( z )( r ) 2 2
4
2
约定：把覆盖(i, j)网格的部分记为4，对角部分记为1， 沿r方向的部分记为2，沿z方向的部分记为3。
* * ˆ* f ˆ* f f f ai i 2 i 1 2 i 3 i 1 x x * ˆ* f ˆ* fi* f 2 f i 1 bi 3 1 2 i i x x
ut
若u<0 则
i-1 i+1， x - x
1D理想流体Euler方程
0) 若 c为常数，则其解为 f(x ,t) f(x ct , ，其差分形式为
f(xi ,t t) f(xi ct,t)
则新时刻的分布化为老时刻某空间点的分布
ct
t
t t
i 1
i
i 1
其求解之最简方法为一次迎风格式（以c >0为例）
t n (fi fin1 ) 但其缺点是强数值耗散。 x CIP方法在两格点间 ( xi 1, xi ) 采用3次多项式拟合 ˆ n 表示为 i 格 f i 3 2 n n ˆ F ( x) ai ( x xi ) bi ( x xi ) ( x xi ) fi fi 点处 f 的导数
cs p /
u un 1 un 1
i 2 i 2
u 0 u 0
0.6 1
20. 高级算法之PIC法
Euler格式：适用于求解含大扭曲、大滑移问题 Lagrange格式：适用于求解撞击界面的多维流体问题。 F. H. Harlow. 1964, Methods in Comput. Phys., 3, 319 发展了PIC (Particle-in-cell) 方法。它综合了Euler格式和Lagrange格式 的优点，在Euler网格上定义场变量，用质点表征流体。主要思想是在 Euler网格中配置离散的Lagrange流体质点。
4
3
3
Hale Waihona Puke 1134
3
2
1
2
1
4
2
2
4
若覆盖的网格为空网格，空网格的速度取为质点所在 网格的速度，或取此网格的面积为0。若覆盖的网格 在固壁外，其速度取为相邻内网格的反射速度。 （2）计算质点的新坐标
r
n 1 k
1 n 1 1 n 1 n n 1 n r (U k U k )t , zk zk (Vk Vkn )t 2 2
网格界面 p 按算术平均值计算，若旁边网格是空网格， 置界面 p=0；若界面是固壁，则取界面p等于网格p 。 取人为粘性
qn 1
i , j 2
( in1, j in, j )u 2 2 0
n n u U U i 1 i u 0 u 0 z向同理
n k
（3）重新分配质点 质点运动后可能留在原网格，也可能移动到新网格， 从而引起网格质量、动量和能量的变化

n 1 i,j
m N n 1 rk (i , j ) rz k 1
N


1 ein, j1由in , j 及状态方程求出
1 v n i,j
m n 1 n 1 v r k k (i , j ) rz k 1
Lucy, 1977, AJ, 82, 1013
Mesh-free
Gingold & Monaghan, 1977, MNRAS, 181, 375
发展了SPH (Smoothed Particle Hydrodynamics光滑粒子流体动力学) 方法: load each particle with a spatially extended interpolation kernel，平均局地密度 缘于所有质点的贡献。自在天体物理研究领域发展以来，已广泛应用于其 它物理问题，如连续体结构的解体、碎裂、固体的层裂、脆性断裂等。它 可以解决很多其它算法解决不了的问题。由于它是无网格的，故可用于研 究极不规则的结构。网上有Liu et al.写的书（电子版）及例程。