光滑质点流体动力学_SPH_方法_综述_

第13卷第4期1996年12月

　计　　算　　物　　理　　

　CHIN ESE J OURNAL OF COMPU TA TIONAL PHYSICS 　　　Dec.

Vol. 13,No. 4

　, 1996

光滑质点流体动力学(SPH ) 方法()

　(, )

摘　要　(SPH ) 方法, 由于它计算, 它对非重点介绍了该方法的理论基础, 流体动力学方程组的推导, , 自引力、汇和磁场, 光滑核的选取, 以及SPH 执行过程等有关问题。关键词　SPH 　光滑核　积分插值中图分类号　　O35112

1　引　言

大多数流体动力学问题, 由于其复杂性都要求数值计算, 为此发展了多种数值方法。其中(SPH 2Smoothed Particle Hydrodynamics ) 方法, 它对天体物有一种名叫“光滑质点流体动力学”理问题尤为合适, 是近廿多年来发展起来的一种新的纯Lagrange 方法。描述SPH 技巧的第

一篇论文是Lucy 于1977年提出的, 经过G ingold 和Moanaghan (GM77, GM82, M82, GM83,M G 83,ML85) 等人的工作, 奠定了扎实的基础, 再经过Benz (B T90, HB91, HBC92, FBD95) 等人的改进和完善, 现已发展成为比较成熟的计算高维(尤其是三维) 的天体物理问题的有效方法了。本文的目的就是对SPH 方法作综述性介绍, 以促进该方法在我国国内得到应用和发展。

SPH 是一种质点方法, 有点类似质点网格法(PIC 2Particle in Cell ) (见Harlow 1957,1974, 1988) , 但根本的不同点在于SPH 方法中计算空间导数时不需要使用任何网格, 而被插值公式中的解析微分式子所替代, 从而避免了高维拉氏差分网格法中的网格缠结(tangling ) 和扭曲(distortion ) 等最令人头痛的问题。这种方法另一个突出的优点是表现在对缺乏对称性(lack symmetry ) 和内含真空(large voids ) 的三维系统的计算特别有效。因为利用网格建立的所有传统的差分方法, 随着网格数目的增加, 到了三维情形变得难以招架, 真空区域的网格剖分造成内存贮量的大量浪费又实在令人可惜, 而SPH 恰能克服上述的缺点。由于受到来自旋转和磁场的非球对称力的影响, 天体物理领域内要研究的大多数问题都是非对称的三维问题, 而SPH 方法正是适应这种需要而发展起来的(M85,M88,M92,B88,B90,B91) , 但也在其它领域内得到广泛的应用, 例如高速碰撞的材料动载响应(见Libersky 等人(1993) ) 的数值模拟。

2　方法的基础

SPH 方法的核心实为是一种插值(interpolation ) 。在SPH 中任一宏观变量(如密度、压力

1995年7月6日收到原稿,1996年8月9日收到修改稿。

386计　　算　　物　　理

第13卷　

温度等) A (r ) 能方便地借助于一组无序点(disordered points ) 上的值表示成积分插值(integral

interpolant ) 计算得到, 其形式为

〈A (r ) 〉=

A (r ) W (r -∫

r ′, h ) d r ′, (211)

其中D 为整个求解区域。W 是插值核(interpolating kernel ) (i ) 　W (r -r ′, h ) d r D

∫

h →0

(212) (213)

(ii ) W (-r ) . ) (这里h 。显然, (i ) 是一个规一化条件, 而(ii ) 实

在, (Strongly peaked function ) . 在|r -r ′|>h , 它为零; →0, 它是一个δ(Delta ) 函数, 即有lim 〈A (r ) 〉=A (r ) 。通常称〈A (r ) 〉为A

h →0

(r ) 的一个核估计, 它亦等价于围绕场A (r ) 带有一个光滑器(Smoothing ) 或过滤器(fiter ) 函

数W (r -r ′, h ) 去产生场的一个估计, 把其中的局部统计涨落都统统地过滤掉。

如果W 选取得仅仅是r 的偶函数(球对称核) , 即W (r , h ) =W (|r |, h ) , 则能证明

(2. 4) 〈A (r ) 〉=A (r ) +C ・( 2A ) h 2+O (h 3) ,

其中C 与h 无关, 并称〈A (r ) 〉具有h 的二阶精度, 记为O (h 2) 。类似地有

(215) 〈() 〉=+O (h 2) .

B r

通常选取的核具有紧纟致支撑性(compact support ) 的偶函数。一般使用的是G aussian 核(GM77) , 亦有用超高斯核(GM82) 、指数核(Wood ,1981) , 样条(spline ) 核或B 2样条核(ML85)

等, 详见本文§6。

～

在不同的论文中采用不同的符号, 例如:A (r ) (L 77) , A S (r ) (GM 77) , 注1:记号〈A (r ) 〉

A k (r ) (M 82) , A I (r ) (M 92) 等等。

如何数值地计算核估计(211) 呢? 最初由Lucy (1977) 提出的思想是利用Monte Carlo 方法求多重积分中的重要抽样技巧。但现在的理解已和原先的有所不同。即初始分布的插值点(或样本点) 分布并不是完全随机地分布, 而是初始取定后随流体的运动而运动(M 85) 。现假设有密度为ρ(r ) 的流体在运动, 将(211) 式的右边改写为

) d r ′(217) W (r -r ′, h ) ρ(r ′. ) D ρ(r ′

为了计算此积分, 我们设想将流体剖分成N 个小体积元(volume elements ) , 分别具有质量为

m 1, m 2, …, m N , 质量的中心位置为r 1, r 2, …, r N 。第j 个体积元对积分的贡献是

()

ρ(r j ) W r -r j , h m j 1

(218)

在数值工作中, 积分往往近似为有限个数目的部分和, 因此有

ρ(r ′) d r ′=∫

j =1

m j =Const. (219)

而积分(217) 近似为

〈A (r ) 〉=

j =1

m j

W (r -r j , h ) , ρj

(2110)

　第4期

张锁春:光滑质点流体动力学(SPH ) 方法(综述)

387

其中A j =A (r j ) 。若为了更简单, 取小体积元的质量相等, 即m j =m (j =1, 2, …, N ) , 则(2110) 变为

〈A (r ) 〉=m

j =1

W (r -r j , h ) , ρj

(2110) ′

这就是一个简单的黎曼(Riemman ) 和。

如果核函数W 是n 次可微的, 则由(211) (r ) , 我们能用一个解析函数〈A (r ) 〉普通求导数得到, , 。

, 由于本身含有质量等信息量, 故称

ρ(r ) , 则由(2110) 或(2110) ′在1A (r ) ≡可得:

〈ρ(r ) 〉=

或

ρ(r ) 〉=m 〈

j =1j =1

m j W (r -r j , h ) , (2111)

6W (r -6

r j , h ) , (2111) ′

这就是说:连续密度场是可以通过一组质点的质量, 经过核光滑化所得1特别, 在空间中标号为i 的质点的密度

ρ〈ρ(r ) 〉i ≡i =

m j W (r i -r j , h ) ,

(2112)

j =1

是被服从空间分布W 的所有对它有贡献的质点的质量所光滑化1这也是本方法取名为“光滑

质点”的原因。

利用核的可微性, 以及函数或核函数在边界上为零的特性(这一点不是本质的, 如果不为零的话, 在推导的公式中加上边界效应项即可) 来计算〈 A 〉〈、・v 〉和〈 ×v 〉. 这里是表示对r 的哈密顿算子, v 为速度。

事实上, 利用定义(211) 有

〈 A 〉=

( A (r ′) W (r -r ′, h ) d r ′∫

) =A (r ′A (r ′

∫) W (r -r ′, h ) n d S -∫

D 5D

W (r ′r -r ′, h ) d r ′

) W (r -A (r ′

∫

r ′, h ) d r ′= 〈A 〉

(2113)

j =1

m j

W (r -r ′j , h ) , ρj

A 在这里n 是单位表面积d S 的外法向, 面积分项取为零, 并利用核函数W 具有 W (r -r ′, h )

=- r ′W (r -r ′, h ) 的特性。在天体物理领域内大多数问题都包含有自引力项, 取无穷远处为边界, 表面积项自动为零。因此我们也可以假定解空间D 取得足够大, 使得在边界5D 上, A W →0。但在有些流体动力学问题中, 边界是包含在流体之中, 边界类型有入口边界(inflow ) , 出口边界(outflow ) 和固壁边界(rigid wall ) 三种, 在SPH 中有关边界项的详细讨论可见Campbell (1988) .

388计　　算　　物　　理

第13卷　

因此利用(2113) , 特别有

ρ〈 ρ〉= 〈〉=〈 P 〉= 〈P 〉=

j =1

m j W (r -r j , h ) , m j

(2114) (2115)

j =1

W (r j , h ) .

ρ注2:有些作者对S PH 〈A 〉的

基础上, 这确实是奥妙之处。

(ρ类似地, 利用关系式・i (・v ) i 〉=-N j =1

v ij ・ i W ij , (2116) (2117)

ρ〈i ( ×v ) i 〉=

j =1

m j v ij × i W ij ,

其中v ij =v i -v j , w ij =W (r i -r j , h ) , i W ij 是表示W (r i -r j , h ) 对r i 求的。

3　流体力学方程组

为了说明SPH 中使用的流体动力学方程组, 我们从最简单的Euler 形式的流体动力学方

程组出发进行讨论。即取

((311) 动量守恒方程:+v ・ ) v =-ρ, 5t (

(312) +v ・ ) ρ=-ρ ・v ,

((313) 能量守恒方程:+v ・ ) u =-ρ ・v , t

其中P 是压力, u 是比内能, 通常的做法是将偏微分方程两边同乘核函数W , 然后在整个求解区域上积分。利用§2中得到的基本关系式, 推导了S PH 中所需要的流体动力学方程组。下面分别说明之。311　动量守恒方程

对于方程(311) , 首先注意到非线性项

(〈+v ・ ) v 〉=〈〉+〈(v ・ ) v 〉=+〈v 〉・〈 v 〉t t t

(31111) =+〈v 〉・〈v 〉=15t d t

这里d /d t 是表示随流体运动的导数, 其次注意到非线性组合项 P/ρ, 它有不同的近似形式, 例如:

①利用(2115) 有

质量守恒方程:

ρ〉〈P 〉i =i =ρρi i

j =1

m j

W ij .

ρj

(3. 1. 2)

Ρ) -P ρ, 则有　　②利用ρ P = (ρ

〈ρ〉i =6m j (P j -P i ) i W ij . ρi j =1

(3. 1. 3)

　第4期

张锁春:光滑质点流体动力学(SPH ) 方法(综述)

389

Ρ　　③利用ρ= (ρ) +2 ρ, 则有

〈ρ〉i =甚至利用更一般的公式

j =1

m j (

ρj

ρi

) i W ij . (3. 1. 4)

=σ (1-σ) +-σ(-1Πρρρ则有

=〈i 〈ρ〉 i W ij . i =26m j

ρρi j j =1

j =N

j (

-j

ρi j

σ2-σ i W ij . (3. 1. 4) ′

(3. 1. 5)

对于③④由于具有对称性, 能保证线性和角动量守恒, 而对①②则不能。

令v i ≡〈v 〉i , 现取(31111) 和(31114) , 则动量守恒方程(311) 就变为

=-d t

j =1

m j (

2ρj

ρi

) i W ij . (3. 1. 6)

注意到

m i (

) d t

=m i m j (

2ρi

ρj

2) i W ij =-m i m j (

ρi

ρj

2) j W ij =-m j (

) . d t i

如果W (r , h ) 是r 的偶函数, 则 i W (r i -r j , h ) 对i 和j 具有反对称性, 故有

i =1j =1

N N

m i m j (

ρi

2ρj

) i W (r i -r j , h ) =0.

因而能确保总动量守恒。

312　质量守恒方程

对方程(312) , 利用(2112) 和(2116) , 则有

ρ=6m j v ij i W ij . d t j =1

313　能量守恒方程

对方程(313) 而言, 对-ρ ・v 项的处理亦有不同的形式, 令u i ≡〈u 〉i , ①利用(2116) 有

=26m j v ij ・ i W ij . d t ρi j =1②利用

() , 则有=- () +v ・ ρρd t

=m j (2) v ij ・ i W ij . d t ρj j =1

(3. 2. 1)

(3. 3. 1)

(3. 3. 2)

③为对称起见, 一般取①和②的平均(M G 83)

=m j (2+2) v ij ・ i W ij . 6d t 2j =1ρj ρi

(3. 3. 3)

390计　　算　　物　　理

第13卷　

④利用

() ) =(-ρ ・v +d t

, 则有

(3. 3. 4)

=2m j v i ・ i W ij .

ρρd t i j j =1

对这四种情形, 只要利用核W 的 i W ij =- j W ij , 。这是因为

d t

i =1

m i u i =

i =1N

m i =

d t

i 1j =N N

i j

2ρi

i -v j ) W i =1

i m j (

j i

) v 恒i ・, i W ij

d t 2

i =1

m i v 2i .

但是u i , 例如③④

但能导致非物理解(如负内能) 的出现, 所以还不如用(31311) 为好。314　质点运动位置

质点运动的位置用

=v i d t

(31411) (31412)

或用

ε=^v i =v i +6m j W ij j =1d t ρij

来决定, 其中ε为0

315　状态方程

(ρi +ρj ) (M 89) 。显然(3. 4. 2) 是对(3. 4. 1) 作了2

根据需要事先给定状态方程。例如理想气体的状态方程是p i =(γ-1) ρ。i u i , 其中γ=3

ρ又如等温状态方程是p i =C 2iso i , 其中C iso 是等温声速, 此时比热能u i 的演化方程就不再需要。

4　人为粘性和热传导

411　人为粘性

大多数流体计算中, 为处理激波间断面采用人为粘性, 此时动量守恒方程为

=-ρ+人为粘性. d t

　　关于人为粘性已有不少作者提出过不同的形式(Lucy 1977, Wood 1981, Monagham &G ingold 1983, Evrard 1988,Loewenstein &Mathews 1984) , 但通常采用的形式是包含有两种不同类型的粘性:

ρ容积(bulk ) 粘性:7e =-αv , lC s ・

von Neum ann 2Richt myer :7

ρ=-βv ) 2, l ( ・

其中α和β是自由参数, 量级～1, l 是激波传播长度尺度, C s 是声速, 因此有人为粘性形式:

μβμ2αij

ρij

v ij ・r ij

　第4期和

张锁春:光滑质点流体动力学(SPH ) 方法(综述)

391

μij =

这里C ij =

・r 2. r 2ij +η

(C i +C j ) , C i 是质点i 在位置r i 的声速。r ij =-r j 。η, 典型22取η=0. 01h 2, 是为了避免数值导数出问题。

Balsara 等人(1991) ,

为此提出了修正形式:

ij 其中

f i =

r r ij +η

f +f 2

|〈・v 〉i |+|〈 ×v 〉i |+0. 001h

改虑到人为粘性, 此时动量守恒方程可取为

=-d t

j =1

m j (

ρj

2ρi

) i W ij , (4. 1. 1)

和能量守恒方程可取为

=2d t ρi

412　热传导

j =1

m j v ij ・ i W ij +

j =1

m j

v ij ij

・ i W ij . (4. 1. 2)

若在能量守恒方程中引进人为热流(热传导) 或真实物理热流(辐射输运) :

() =- ・v +粘性热+

ρρ ・F r d t

(4. 2. 1)

其中F (r ) 是热流。

人为热流:F art (r ) =-K u ,

Τ物理热流:F rad (r ) =-ρ , 3K

这里　　K =k 1ρl C s +k 2ρv |- ・v ) , k 1, k 2为常数, c 为光速, a 为斯忒藩—玻耳兹l (| ・

曼常数。Brookshaw (1986) ,Monaghan (1988) 认为取下列形式较好:

〈ρ ・(K u ) 〉i =

其中q =ρ, 对 T 的处理与 u 相同。

j =1

m j

() (・ )

, 22

ρij (r ij +η)

5　自引力、汇和磁场

本节是专门为SPH 在天体物理领域内的应用而写的, 正如在§1中的指出的那样,SPH

是为了适应天体物理问题的需要而发展起来的, 而Benz (1988) 曾列举了12个方面的应用, Monaghan (1992) 却列举出12个更大方面的应用, 其中不可避免地要涉及到自引力、汇和磁场

等方面的问题, 下面简单地叙述之:

392计　　算　　物　　理

第13卷　

511　自引力

在大多数天体物理领域内的问题都必须考虑自引力, 此时动量守恒方程为

φ=-ρ P - d t

(5. 1. 1)

其中φ为引力势, 是由下列的Poisson 方程决定的 2φ=4πG 这里G 是引力常数。

在SPH 中〈, φ〉= 〈φ,

φ〉ρ2=4πG 〈〉.

(51112)

(51113)

而(511〈φ(r ) 〉=-G

其中K (|r -r j |, h ) =

--j =1

m j K (|r -r j |, h ) ,

(5, 1, 4)

=-d t 512　汇

j =1

W (|r -r |, h )

d r ′. 结合(31116) , 此时(51111) 可变为

|r -r ′|

N -(5. 1. 5) m j (2+2) i W (r i -r j , h ) +G m j i K (|r i -r j |, h ) .

ρj

ρi

j =1

在有些天体物理问题中, 要考虑物质不断地抛射, 质量的损失在连续性方程中看作是一个汇。

(51211) =-ρ ・v -f (r ) , d t 此时有

ρ=d t

513　磁场

j =1

m j v ij ・ W ij -

j =1

m j

W ij .

ρj

(5. 2. 2)

关于磁流体动力学方程组(MHD ) 的SPH 公式已经过(G ingold &Monaghan (1977) ,

Philips &Monaghan (1985) ,Philips (1986) ,Habe (1989) 和Stellingwerf (1990) 等人研究比较成熟。由于较复杂在这里不作深入讨论, 仅以最简单的模型为例说明之。将某离子体看作一种单一的流体, 在一个磁场B 和电场E 中带有电流J , 忽略静电力, 在Euler 坐标下的质量、动量和能量守恒方程组为

+v ・ ) ρ=-ρ・ v , 5t

ρ(+v ・ ) v =- p +J ×B , 5t

ρ(+v ・ ) u =-p ・u +(E +v ×B ) ・J , t

(

和

・B =0, μ0J = ×B , + ×E =0, t

　第4期

张锁春:光滑质点流体动力学(SPH ) 方法(综述)

E +v ×B =ηJ 1

其中μ0是自由空间的磁导率, η为电阻率。根据(2117) 式, 我们有

〈ρ( ×B ) 〉i =

j B ij × i W ij .

j 6

m =1

对磁场B 随时间的变化亦有不同的形式:

①利用

(d t ρ

) =ρ

v =

・ ρv -) 〉, N

ρ2j v , k i , k ) B i ・ i W ij .

i j 1

其中A j i k 个分量。

②d t

=-B ( ・v ) +(B ・ ) v , 则有N

d t =ρi 6

m j (B i , k v ij -v ij , k B i ) ・ i W ij . j =1

　　③亦可将上式改写成另一形式

d t =ρi

v ij ×B i ) × W ij j 6

m j (=1

6　关于光滑核

611　核的形式

从现有已发表的论文看, 大致有以下几种:

①G aussian (GM77)

W (r , h ) =

h e

-s

, 　O (h j

) .

其中d 是空间维数, s =

. ②Supper 2G aussian (GM82) 成熟W (r , h ) =

h (

-s 2) e

-, 　O (h 4

) 当s >

时, W 为负值, 导致非物理解。③指数核(Wood 1981) (三维形式)

W (r , h ) =

8πh

-s

. 　　④Spline (ML85)

S -×1S 2+1, 0≤S ≤1W (r , h ) =

d (

4(2-S ) 3, 1≤S ≤2

　其它

393

(61111)

研究(61112)

(61113)

(61114)

394计　　算　　物　　理

第13卷　

其中C 为归一化常数, 在一维, 二维和三维情形分别为

⑤B 2Spline (ML85)

, 和π, 其精度为O (h 4) . π37

0≤S ≤11≤(61115)

W (r , h ) =

24πh 3(2-S ) (-) ,

-7S +4S , 3

　612　光滑长度h 的选取

, h 是常数, 它不随时间和空间而变化。

? 最初Gi ngol d &Monaghan (1977) 是从(2110) 作为(2117) , 选取h 使方差尽可能地小。利用核函数W 在r ′=r 上是强尖峰函数的性质, 作Taylor 级数开展, 保留主项得到一个近似泛函L (r ) , 选取h 使得此泛函达到极小, 得到h ∝

1/7

, 如果取N =40试算, 取G aussian 核和样条核其结果差别不大, 但N =80

时, G aussian 核要比样条核更精确。

选定G aussian 核后, 可利用Newton 2Raphson 方法去解方程

i =1j =1

N N

2exp -(r i -r j ) /h =Const ,

(6. 2. 1)

来决定h , 表明有h j ∝. 现在一般取

ρj d

h =h 0(

ρ) d ρ,

其中ρ0和h 0是常数。

6　滑核=-d t

=-d t

= ・v ,

d ρd t d =-d ρd t i

d ρi

j =1

(6. 2. 2)

m j v ij ・ i W ij .

这表明h 是在空间h =h (r ) 上随时间的变化。

) ) 替代W (r -r ′) =Evrard (1988) 曾考虑用W (r -r ′, h (r , r ′, h ) , 其中h (r , r ′

～

) ) , 对应的质点公式为h ij =(h i +h j ) , W ij =W (r ij , h ij ) . +h (r ′

Hernquist 和Katz (1989) 曾主张用

～

(h (r ) 2

W ij =

[W (r i -r j , h i ) +W (r i -r j , h j ) ].2

　　如果取h =h (t ) 或h =h (r , t ) , 此时对上述的SPH 中的有关公式需作重新推导(见Monaghan 1992) 。

　第4期张锁春:光滑质点流体动力学(SPH ) 方法(综述) 3957　SPH 执行过程

711　SPH 方程组

根据上述的推导, 可假定取

=-d t

=-d t j =16N m j (2ρj +2ρi ++～) , ～(7. 1. 1) (7. 1. 2) (7. 1. 3) 1N m j (+) ij ・ i W ij , j =1～m j v ij ・ i W ij , d ρi j =16N ～m j v ij ・ i W ij , on (7. 1. 4)

(7. 1. 5) =v i . d t

求解的方程视要求的问题而定。

712　质点设置

在初始, 质点设置是确定质点的空间位置和质量。位置随时间演化而改变, 但质量不变。配置的原则一是尽量保持空间分布的均匀性, 在规则网格上常常会带来方便性, 如在三维情形, 质点中心位置取在正方形网格的中心; 二是质点分布在空间中要有一定的稠密度, 各质点之间的间隔一般要求≤h 。当然对于不同构形应采用不同的配置方式。

713　时间积分

对常微分方程组(71111) ～(7. 1. 5) 的数值积分, 一般采用显式格式。例如用标准的蛙跳(leapfrog ) 算法(H K 89) , 或用预估—校正(Predictor 2Corrector ) 方法(Lattanzio 等1986) Benz (1984) 提出了变步长的二阶Runge 2Kutta 方法(Fehlberg (1968) 。时间步长的大小要求满足CFL (1928) 条件。对有作用力项和粘性扩散项的情况(M 89) 为

δ(71311) t =λmin (δt f , δt c v ) ,

其中

ar r ) |f i |μij . C i +016(αC i +βmax j δt f =min i , δt c v =min i

这里δt f 是基于每单位质量的力f , 而δt c v 是Courant 和粘性对时间步长控制的组合, λ

在有些天体物理计算中, 要求能量方程求解必须采用隐式格式。此时就要求对各个质点求解非线性代数方程组(MV88,H K89,ML91) 。

最后, 在此进一步指出本文起抛砖引玉的作用。随着高速大型计算机的发展, 与并行计算有关的自适应SPH 以及与PP (particle 2particle ) 相耦合的PPASPH (Serna 等1996) 等方法相继出现, 故发展前景十分喜人, 可喜的是SPH 已在超新星爆发模拟中获得了巨大的成功[56], 而国内国防科技大学的贝新源等人对SPH 方法的应用也作了有益的尝试[57]。

396计　　算　　物　　理第13卷　参考文献

1　Bandiera R. Astron Astrophys (A &A ) 1984,139:368

2　Benz W. A &A, 1984, 139:378

3　Benz W. Com put Phys Com m (CPC ) 1988, 48.

4　Benz W. in ”the Numerical Modelling of Nonlinear (ed J ) Aad Pub

1990, p269.

5　Benz W. in ”Late Steges of in Astrophysical Hydrodynamics ”(ed

DeLoore C ) 2, 6　Benz K. J L ett , 1990, 348:117.

7　The of Binary Systems and Rotating Stars , Ph D thesis , Monash Univ , 1986.

8　Campbell M. S ome New Algorithms for Boundary Value Problems in SPH , Albuquergue :Mission Res Corp ,

1988.

9　Chevelier R A. Ast ro phys J (A pJ ) , 1976, 207:872.

10　Courant R , Friedrichs K O , and Lewy H. M ath A nn , 1928, 100:32.

11　Evrard A E. Mounthly Notices Roy. Ast ron. Soc. (MNRAS ) 1988, .

12　Falk S E , and Arnett W D. A pJ , 1973, 180:L65.

13　Fulbright M S , Benz W , and Davies M B. A pJ , 1995, 440:254.

14　G ingold R A , and Monaghan J J. MNRAS , 1977, 181:375.

1516Energy Phys , 1986

18　Haddad B , Clausset F , and Combes F. JCP , 1991, 97:103.

19　Harlow F H. J Assoc Com put M ach , 1957, 4:137.

202123. Methods Com put Phys , 1964, 3:319.. CPC , 1988, 48:1.. A pJ , 1992, 387:294. . J Com put Phys (JCP ) , 1982, 46:429.. MNRAS , 1983, 204:715.17　Habe A. in ”Status Rep. Super Computing Japen ”, (ed Nakamura T , and Nagasawa M ) National Lab High 22　Herant M , and Benz W. A pJ , 1991, 370:L81.

24　Herant M , Benz W , and Colgate S. A pJ ,1992, 395:642.

25　Herant M , and Woosley S E. A pJ , 1994, 425:814.

26　Hernquist L. ApJ , 1993, 404:717.

27　Hemquist L ,and K atz N. A P JS , 1989, 70:419.

28　Hochkhey R W , and Eastwood J N. Computer Simulation Using Particles , Mc Graw 2Hill , New Y ork , 1981. 29　Lattanzio J C , Monaghan J J , Pongracic H , and Schwarz M Z. S IA M J Sci S tat Com put , 1986, 7:591. 30　Libersky L D , Petschek A G , and Allandadi F A. JCP , 1993, 109:67.

31　Loewenstein M , and Mathews W G. JCP , 1984, 62:414.

32　Lucy L B. Ast ron J , 1977, 82:1013.

33　Monaghan J J. S IA M J Sci S tat Com put , 1982, 3:422.

3435. Com put Phys Rep (CPR ) , 1985, 3:71. . CPC , 1988, 48:89.

　第4期

363738张锁春:光滑质点流体动力学(SPH ) 方法(综述) . JCP , 1989, 82:1. . Preprint , Dept. of Math. , Monash Univ. , Clayton , Victoria , 1990,3168. 397. A nnual Review of Ast ronomy and Ast rophysics (A RA &A ) 1992, 30:543. 39　Monaghan J J , and G ingold R A. JCP , 1983, 52:374.

40　Monagnan J J , and Lattanzio J C. A &A , 1985, 149:135.

41. A pJ , 1991, 375:177.

42　Monaghan J J , and Varnas S. M N RA S , 1988, 43　Nagasawa M , Nakamura T , and S , 40:691.

44　Navarro J F , and RA S , :45　Nelson , J N S , 1993, 264:905.

46　Nob , , :78.

47　, and Libersky L D. JCP , 1993, 109:76.

48　Phillips G J. M N RA S , 1986, 221:571.

49　Phillips G J , and Monaghan J J. M N RA S , 1985, 216:883.

50　Rasio F A , and Snapiro S. A pJ ,1991, 377:559.

51　Serna A , Alimi J 2M , and Chieze J 2P. A pJ , 1996, 461:884.

52　Serna A , Alimi J 2M , and Scholl H. A pJ , 1994, 427:574.

53　Stellingwerf R F. Blast Wave Stability in Nuclear Explosion , Albuquergue :MissionRes. Corp. , 1990.

54　Weaver T A , and Woosley S E. in ”Supernova S pectra ”, AIP Conference Proceedings No 63(eds. Meyerott ,

R , G illespie G H ) , New Y ork :AIP, 1980.

55　Wood D. M N RA S , 1981, 194:201.

56　张锁春, 王贻仁, 汪惟中, 谢佐恒1恒星演化晚期阶段研究中的流体动力学计算方法1天文学进展,1996,

14:1491

57　蒋伯诚, 张锁春主编1高科技研究中的数值计算1长沙:国防科技大学出版社,19951P68-731

SMOOTHED PARTICL E H YD RODY NAMICS (SPH )

METH OD (A REVIEW)

Zhang Suochun

(Instit ute of A pplied M athem atics , A cadem ia S i nica , Beiji ng 100080)

ABSTATCT 　This paper describes a new and pure Lagrangian method ———called ”Smoothed Particle Hydrodyna 2mics ”(SPH ) method. The method is to actually evaluate spatial gradients without the use of any grid. Thus it does not suffer form the severe problems always associated with mesh tangling and distortion. Therefore it can be applied to multidimensional hydrodynamics which could effectively model three 2dimensional systems which lack symmetry and possess large voids. At first , the paper gives an introduction to theoretical basis to SPH. Emphasis is given to a proper derivation of the SPH equations from the hydrodynamical conservation equations. Discussion covers some rela 2tive problems such as artifical viscosity , thermal conduction , self 2gravity and sink and magnetic fields , choosing the smoothing kernels , and implementation of SPH code , etc.

KE Y WOR DS 　SPH ; Smoothing K ernel ; Integral Interpolant.

第13卷第4期1996年12月

　计　　算　　物　　理　　

　CHIN ESE J OURNAL OF COMPU TA TIONAL PHYSICS 　　　Dec.

Vol. 13,No. 4

　, 1996

光滑质点流体动力学(SPH ) 方法()

　(, )

1　引　言

2　方法的基础

SPH 方法的核心实为是一种插值(interpolation ) 。在SPH 中任一宏观变量(如密度、压力

1995年7月6日收到原稿,1996年8月9日收到修改稿。

386计　　算　　物　　理

第13卷　

温度等) A (r ) 能方便地借助于一组无序点(disordered points ) 上的值表示成积分插值(integral

interpolant ) 计算得到, 其形式为

〈A (r ) 〉=

A (r ) W (r -∫

r ′, h ) d r ′, (211)

其中D 为整个求解区域。W 是插值核(interpolating kernel ) (i ) 　W (r -r ′, h ) d r D

∫

h →0

(212) (213)

(ii ) W (-r ) . ) (这里h 。显然, (i ) 是一个规一化条件, 而(ii ) 实

在, (Strongly peaked function ) . 在|r -r ′|>h , 它为零; →0, 它是一个δ(Delta ) 函数, 即有lim 〈A (r ) 〉=A (r ) 。通常称〈A (r ) 〉为A

h →0

(r ) 的一个核估计, 它亦等价于围绕场A (r ) 带有一个光滑器(Smoothing ) 或过滤器(fiter ) 函

数W (r -r ′, h ) 去产生场的一个估计, 把其中的局部统计涨落都统统地过滤掉。

如果W 选取得仅仅是r 的偶函数(球对称核) , 即W (r , h ) =W (|r |, h ) , 则能证明

(2. 4) 〈A (r ) 〉=A (r ) +C ・( 2A ) h 2+O (h 3) ,

其中C 与h 无关, 并称〈A (r ) 〉具有h 的二阶精度, 记为O (h 2) 。类似地有

(215) 〈() 〉=+O (h 2) .

B r

等, 详见本文§6。

～

在不同的论文中采用不同的符号, 例如:A (r ) (L 77) , A S (r ) (GM 77) , 注1:记号〈A (r ) 〉

A k (r ) (M 82) , A I (r ) (M 92) 等等。

) d r ′(217) W (r -r ′, h ) ρ(r ′. ) D ρ(r ′

为了计算此积分, 我们设想将流体剖分成N 个小体积元(volume elements ) , 分别具有质量为

m 1, m 2, …, m N , 质量的中心位置为r 1, r 2, …, r N 。第j 个体积元对积分的贡献是

()

ρ(r j ) W r -r j , h m j 1

(218)

在数值工作中, 积分往往近似为有限个数目的部分和, 因此有

ρ(r ′) d r ′=∫

j =1

m j =Const. (219)

而积分(217) 近似为

〈A (r ) 〉=

j =1

m j

W (r -r j , h ) , ρj

(2110)

　第4期

张锁春:光滑质点流体动力学(SPH ) 方法(综述)

387

其中A j =A (r j ) 。若为了更简单, 取小体积元的质量相等, 即m j =m (j =1, 2, …, N ) , 则(2110) 变为

〈A (r ) 〉=m

j =1

W (r -r j , h ) , ρj

(2110) ′

这就是一个简单的黎曼(Riemman ) 和。

如果核函数W 是n 次可微的, 则由(211) (r ) , 我们能用一个解析函数〈A (r ) 〉普通求导数得到, , 。

, 由于本身含有质量等信息量, 故称

ρ(r ) , 则由(2110) 或(2110) ′在1A (r ) ≡可得:

〈ρ(r ) 〉=

或

ρ(r ) 〉=m 〈

j =1j =1

m j W (r -r j , h ) , (2111)

6W (r -6

r j , h ) , (2111) ′

这就是说:连续密度场是可以通过一组质点的质量, 经过核光滑化所得1特别, 在空间中标号为i 的质点的密度

ρ〈ρ(r ) 〉i ≡i =

m j W (r i -r j , h ) ,

(2112)

j =1

是被服从空间分布W 的所有对它有贡献的质点的质量所光滑化1这也是本方法取名为“光滑

质点”的原因。

事实上, 利用定义(211) 有

〈 A 〉=

( A (r ′) W (r -r ′, h ) d r ′∫

) =A (r ′A (r ′

∫) W (r -r ′, h ) n d S -∫

D 5D

W (r ′r -r ′, h ) d r ′

) W (r -A (r ′

∫

r ′, h ) d r ′= 〈A 〉

(2113)

j =1

m j

W (r -r ′j , h ) , ρj

A 在这里n 是单位表面积d S 的外法向, 面积分项取为零, 并利用核函数W 具有 W (r -r ′, h )

388计　　算　　物　　理

第13卷　

因此利用(2113) , 特别有

ρ〈 ρ〉= 〈〉=〈 P 〉= 〈P 〉=

j =1

m j W (r -r j , h ) , m j

(2114) (2115)

j =1

W (r j , h ) .

ρ注2:有些作者对S PH 〈A 〉的

基础上, 这确实是奥妙之处。

(ρ类似地, 利用关系式・i (・v ) i 〉=-N j =1

v ij ・ i W ij , (2116) (2117)

ρ〈i ( ×v ) i 〉=

j =1

m j v ij × i W ij ,

其中v ij =v i -v j , w ij =W (r i -r j , h ) , i W ij 是表示W (r i -r j , h ) 对r i 求的。

3　流体力学方程组

为了说明SPH 中使用的流体动力学方程组, 我们从最简单的Euler 形式的流体动力学方

程组出发进行讨论。即取

((311) 动量守恒方程:+v ・ ) v =-ρ, 5t (

(312) +v ・ ) ρ=-ρ ・v ,

((313) 能量守恒方程:+v ・ ) u =-ρ ・v , t

对于方程(311) , 首先注意到非线性项

(〈+v ・ ) v 〉=〈〉+〈(v ・ ) v 〉=+〈v 〉・〈 v 〉t t t

(31111) =+〈v 〉・〈v 〉=15t d t

这里d /d t 是表示随流体运动的导数, 其次注意到非线性组合项 P/ρ, 它有不同的近似形式, 例如:

①利用(2115) 有

质量守恒方程:

ρ〉〈P 〉i =i =ρρi i

j =1

m j

W ij .

ρj

(3. 1. 2)

Ρ) -P ρ, 则有　　②利用ρ P = (ρ

〈ρ〉i =6m j (P j -P i ) i W ij . ρi j =1

(3. 1. 3)

　第4期

张锁春:光滑质点流体动力学(SPH ) 方法(综述)

389

Ρ　　③利用ρ= (ρ) +2 ρ, 则有

〈ρ〉i =甚至利用更一般的公式

j =1

m j (

ρj

ρi

) i W ij . (3. 1. 4)

=σ (1-σ) +-σ(-1Πρρρ则有

=〈i 〈ρ〉 i W ij . i =26m j

ρρi j j =1

j =N

j (

-j

ρi j

σ2-σ i W ij . (3. 1. 4) ′

(3. 1. 5)

对于③④由于具有对称性, 能保证线性和角动量守恒, 而对①②则不能。

令v i ≡〈v 〉i , 现取(31111) 和(31114) , 则动量守恒方程(311) 就变为

=-d t

j =1

m j (

2ρj

ρi

) i W ij . (3. 1. 6)

注意到

m i (

) d t

=m i m j (

2ρi

ρj

2) i W ij =-m i m j (

ρi

ρj

2) j W ij =-m j (

) . d t i

如果W (r , h ) 是r 的偶函数, 则 i W (r i -r j , h ) 对i 和j 具有反对称性, 故有

i =1j =1

N N

m i m j (

ρi

2ρj

) i W (r i -r j , h ) =0.

因而能确保总动量守恒。

312　质量守恒方程

对方程(312) , 利用(2112) 和(2116) , 则有

ρ=6m j v ij i W ij . d t j =1

313　能量守恒方程

对方程(313) 而言, 对-ρ ・v 项的处理亦有不同的形式, 令u i ≡〈u 〉i , ①利用(2116) 有

=26m j v ij ・ i W ij . d t ρi j =1②利用

() , 则有=- () +v ・ ρρd t

=m j (2) v ij ・ i W ij . d t ρj j =1

(3. 2. 1)

(3. 3. 1)

(3. 3. 2)

③为对称起见, 一般取①和②的平均(M G 83)

=m j (2+2) v ij ・ i W ij . 6d t 2j =1ρj ρi

(3. 3. 3)

390计　　算　　物　　理

第13卷　

④利用

() ) =(-ρ ・v +d t

, 则有

(3. 3. 4)

=2m j v i ・ i W ij .

ρρd t i j j =1

对这四种情形, 只要利用核W 的 i W ij =- j W ij , 。这是因为

d t

i =1

m i u i =

i =1N

m i =

d t

i 1j =N N

i j

2ρi

i -v j ) W i =1

i m j (

j i

) v 恒i ・, i W ij

d t 2

i =1

m i v 2i .

但是u i , 例如③④

但能导致非物理解(如负内能) 的出现, 所以还不如用(31311) 为好。314　质点运动位置

质点运动的位置用

=v i d t

(31411) (31412)

或用

ε=^v i =v i +6m j W ij j =1d t ρij

来决定, 其中ε为0

315　状态方程

(ρi +ρj ) (M 89) 。显然(3. 4. 2) 是对(3. 4. 1) 作了2

根据需要事先给定状态方程。例如理想气体的状态方程是p i =(γ-1) ρ。i u i , 其中γ=3

ρ又如等温状态方程是p i =C 2iso i , 其中C iso 是等温声速, 此时比热能u i 的演化方程就不再需要。

4　人为粘性和热传导

411　人为粘性

大多数流体计算中, 为处理激波间断面采用人为粘性, 此时动量守恒方程为

=-ρ+人为粘性. d t

ρ容积(bulk ) 粘性:7e =-αv , lC s ・

von Neum ann 2Richt myer :7

ρ=-βv ) 2, l ( ・

其中α和β是自由参数, 量级～1, l 是激波传播长度尺度, C s 是声速, 因此有人为粘性形式:

μβμ2αij

ρij

v ij ・r ij

　第4期和

张锁春:光滑质点流体动力学(SPH ) 方法(综述)

391

μij =

这里C ij =

・r 2. r 2ij +η

(C i +C j ) , C i 是质点i 在位置r i 的声速。r ij =-r j 。η, 典型22取η=0. 01h 2, 是为了避免数值导数出问题。

Balsara 等人(1991) ,

为此提出了修正形式:

ij 其中

f i =

r r ij +η

f +f 2

|〈・v 〉i |+|〈 ×v 〉i |+0. 001h

改虑到人为粘性, 此时动量守恒方程可取为

=-d t

j =1

m j (

ρj

2ρi

) i W ij , (4. 1. 1)

和能量守恒方程可取为

=2d t ρi

412　热传导

j =1

m j v ij ・ i W ij +

j =1

m j

v ij ij

・ i W ij . (4. 1. 2)

若在能量守恒方程中引进人为热流(热传导) 或真实物理热流(辐射输运) :

() =- ・v +粘性热+

ρρ ・F r d t

(4. 2. 1)

其中F (r ) 是热流。

人为热流:F art (r ) =-K u ,

Τ物理热流:F rad (r ) =-ρ , 3K

这里　　K =k 1ρl C s +k 2ρv |- ・v ) , k 1, k 2为常数, c 为光速, a 为斯忒藩—玻耳兹l (| ・

曼常数。Brookshaw (1986) ,Monaghan (1988) 认为取下列形式较好:

〈ρ ・(K u ) 〉i =

其中q =ρ, 对 T 的处理与 u 相同。

j =1

m j

() (・ )

, 22

ρij (r ij +η)

5　自引力、汇和磁场

本节是专门为SPH 在天体物理领域内的应用而写的, 正如在§1中的指出的那样,SPH

等方面的问题, 下面简单地叙述之:

392计　　算　　物　　理

第13卷　

511　自引力

在大多数天体物理领域内的问题都必须考虑自引力, 此时动量守恒方程为

φ=-ρ P - d t

(5. 1. 1)

其中φ为引力势, 是由下列的Poisson 方程决定的 2φ=4πG 这里G 是引力常数。

在SPH 中〈, φ〉= 〈φ,

φ〉ρ2=4πG 〈〉.

(51112)

(51113)

而(511〈φ(r ) 〉=-G

其中K (|r -r j |, h ) =

--j =1

m j K (|r -r j |, h ) ,

(5, 1, 4)

=-d t 512　汇

j =1

W (|r -r |, h )

d r ′. 结合(31116) , 此时(51111) 可变为

|r -r ′|

N -(5. 1. 5) m j (2+2) i W (r i -r j , h ) +G m j i K (|r i -r j |, h ) .

ρj

ρi

j =1

在有些天体物理问题中, 要考虑物质不断地抛射, 质量的损失在连续性方程中看作是一个汇。

(51211) =-ρ ・v -f (r ) , d t 此时有

ρ=d t

513　磁场

j =1

m j v ij ・ W ij -

j =1

m j

W ij .

ρj

(5. 2. 2)

关于磁流体动力学方程组(MHD ) 的SPH 公式已经过(G ingold &Monaghan (1977) ,

+v ・ ) ρ=-ρ・ v , 5t

ρ(+v ・ ) v =- p +J ×B , 5t

ρ(+v ・ ) u =-p ・u +(E +v ×B ) ・J , t

(

和

・B =0, μ0J = ×B , + ×E =0, t

　第4期

张锁春:光滑质点流体动力学(SPH ) 方法(综述)

E +v ×B =ηJ 1

其中μ0是自由空间的磁导率, η为电阻率。根据(2117) 式, 我们有

〈ρ( ×B ) 〉i =

j B ij × i W ij .

j 6

m =1

对磁场B 随时间的变化亦有不同的形式:

①利用

(d t ρ

) =ρ

v =

・ ρv -) 〉, N

ρ2j v , k i , k ) B i ・ i W ij .

i j 1

其中A j i k 个分量。

②d t

=-B ( ・v ) +(B ・ ) v , 则有N

d t =ρi 6

m j (B i , k v ij -v ij , k B i ) ・ i W ij . j =1

　　③亦可将上式改写成另一形式

d t =ρi

v ij ×B i ) × W ij j 6

m j (=1

6　关于光滑核

611　核的形式

从现有已发表的论文看, 大致有以下几种:

①G aussian (GM77)

W (r , h ) =

h e

-s

, 　O (h j

) .

其中d 是空间维数, s =

. ②Supper 2G aussian (GM82) 成熟W (r , h ) =

h (

-s 2) e

-, 　O (h 4

) 当s >

时, W 为负值, 导致非物理解。③指数核(Wood 1981) (三维形式)

W (r , h ) =

8πh

-s

. 　　④Spline (ML85)

S -×1S 2+1, 0≤S ≤1W (r , h ) =

d (

4(2-S ) 3, 1≤S ≤2

　其它

393

(61111)

研究(61112)

(61113)

(61114)

394计　　算　　物　　理

第13卷　

其中C 为归一化常数, 在一维, 二维和三维情形分别为

⑤B 2Spline (ML85)

, 和π, 其精度为O (h 4) . π37

0≤S ≤11≤(61115)

W (r , h ) =

24πh 3(2-S ) (-) ,

-7S +4S , 3

　612　光滑长度h 的选取

, h 是常数, 它不随时间和空间而变化。

1/7

, 如果取N =40试算, 取G aussian 核和样条核其结果差别不大, 但N =80

时, G aussian 核要比样条核更精确。

选定G aussian 核后, 可利用Newton 2Raphson 方法去解方程

i =1j =1

N N

2exp -(r i -r j ) /h =Const ,

(6. 2. 1)

来决定h , 表明有h j ∝. 现在一般取

ρj d

h =h 0(

ρ) d ρ,

其中ρ0和h 0是常数。

6　滑核=-d t

=-d t

= ・v ,

d ρd t d =-d ρd t i

d ρi

j =1

(6. 2. 2)

m j v ij ・ i W ij .

这表明h 是在空间h =h (r ) 上随时间的变化。

) ) 替代W (r -r ′) =Evrard (1988) 曾考虑用W (r -r ′, h (r , r ′, h ) , 其中h (r , r ′

～

) ) , 对应的质点公式为h ij =(h i +h j ) , W ij =W (r ij , h ij ) . +h (r ′

Hernquist 和Katz (1989) 曾主张用

～

(h (r ) 2

W ij =

[W (r i -r j , h i ) +W (r i -r j , h j ) ].2

　　如果取h =h (t ) 或h =h (r , t ) , 此时对上述的SPH 中的有关公式需作重新推导(见Monaghan 1992) 。

　第4期张锁春:光滑质点流体动力学(SPH ) 方法(综述) 3957　SPH 执行过程

711　SPH 方程组

根据上述的推导, 可假定取

=-d t

=-d t j =16N m j (2ρj +2ρi ++～) , ～(7. 1. 1) (7. 1. 2) (7. 1. 3) 1N m j (+) ij ・ i W ij , j =1～m j v ij ・ i W ij , d ρi j =16N ～m j v ij ・ i W ij , on (7. 1. 4)

(7. 1. 5) =v i . d t

求解的方程视要求的问题而定。

712　质点设置

713　时间积分

δ(71311) t =λmin (δt f , δt c v ) ,

其中

ar r ) |f i |μij . C i +016(αC i +βmax j δt f =min i , δt c v =min i

这里δt f 是基于每单位质量的力f , 而δt c v 是Courant 和粘性对时间步长控制的组合, λ

在有些天体物理计算中, 要求能量方程求解必须采用隐式格式。此时就要求对各个质点求解非线性代数方程组(MV88,H K89,ML91) 。

396计　　算　　物　　理第13卷　参考文献

1　Bandiera R. Astron Astrophys (A &A ) 1984,139:368

2　Benz W. A &A, 1984, 139:378

3　Benz W. Com put Phys Com m (CPC ) 1988, 48.

4　Benz W. in ”the Numerical Modelling of Nonlinear (ed J ) Aad Pub

1990, p269.

5　Benz W. in ”Late Steges of in Astrophysical Hydrodynamics ”(ed

DeLoore C ) 2, 6　Benz K. J L ett , 1990, 348:117.

7　The of Binary Systems and Rotating Stars , Ph D thesis , Monash Univ , 1986.

8　Campbell M. S ome New Algorithms for Boundary Value Problems in SPH , Albuquergue :Mission Res Corp ,

1988.

9　Chevelier R A. Ast ro phys J (A pJ ) , 1976, 207:872.

10　Courant R , Friedrichs K O , and Lewy H. M ath A nn , 1928, 100:32.

11　Evrard A E. Mounthly Notices Roy. Ast ron. Soc. (MNRAS ) 1988, .

12　Falk S E , and Arnett W D. A pJ , 1973, 180:L65.

13　Fulbright M S , Benz W , and Davies M B. A pJ , 1995, 440:254.

14　G ingold R A , and Monaghan J J. MNRAS , 1977, 181:375.

1516Energy Phys , 1986

18　Haddad B , Clausset F , and Combes F. JCP , 1991, 97:103.

19　Harlow F H. J Assoc Com put M ach , 1957, 4:137.

24　Herant M , Benz W , and Colgate S. A pJ ,1992, 395:642.

25　Herant M , and Woosley S E. A pJ , 1994, 425:814.

26　Hernquist L. ApJ , 1993, 404:717.

27　Hemquist L ,and K atz N. A P JS , 1989, 70:419.

31　Loewenstein M , and Mathews W G. JCP , 1984, 62:414.

32　Lucy L B. Ast ron J , 1977, 82:1013.

33　Monaghan J J. S IA M J Sci S tat Com put , 1982, 3:422.

3435. Com put Phys Rep (CPR ) , 1985, 3:71. . CPC , 1988, 48:89.

　第4期

40　Monagnan J J , and Lattanzio J C. A &A , 1985, 149:135.

41. A pJ , 1991, 375:177.

42　Monaghan J J , and Varnas S. M N RA S , 1988, 43　Nagasawa M , Nakamura T , and S , 40:691.

44　Navarro J F , and RA S , :45　Nelson , J N S , 1993, 264:905.

46　Nob , , :78.

47　, and Libersky L D. JCP , 1993, 109:76.

48　Phillips G J. M N RA S , 1986, 221:571.

49　Phillips G J , and Monaghan J J. M N RA S , 1985, 216:883.

50　Rasio F A , and Snapiro S. A pJ ,1991, 377:559.

51　Serna A , Alimi J 2M , and Chieze J 2P. A pJ , 1996, 461:884.

52　Serna A , Alimi J 2M , and Scholl H. A pJ , 1994, 427:574.

53　Stellingwerf R F. Blast Wave Stability in Nuclear Explosion , Albuquergue :MissionRes. Corp. , 1990.

54　Weaver T A , and Woosley S E. in ”Supernova S pectra ”, AIP Conference Proceedings No 63(eds. Meyerott ,

R , G illespie G H ) , New Y ork :AIP, 1980.

55　Wood D. M N RA S , 1981, 194:201.

56　张锁春, 王贻仁, 汪惟中, 谢佐恒1恒星演化晚期阶段研究中的流体动力学计算方法1天文学进展,1996,

14:1491

57　蒋伯诚, 张锁春主编1高科技研究中的数值计算1长沙:国防科技大学出版社,19951P68-731

SMOOTHED PARTICL E H YD RODY NAMICS (SPH )

METH OD (A REVIEW)

Zhang Suochun

(Instit ute of A pplied M athem atics , A cadem ia S i nica , Beiji ng 100080)

KE Y WOR DS 　SPH ; Smoothing K ernel ; Integral Interpolant.

光滑质点流体动力学_SPH_方法_综述_

相关文章