We model a massless viscous disk using Smoothed Particle Hydrodynamics (SPH) and note that it evolves according to the Lynden-Bell \& Pringle theory (1974) until a non-axisymmetric instability develops at the inner edge of the disk. This instability may have the same origin as the instability of initially axisymmetric viscous disks discussed by Lyubarskij et al. (1994). To clarify the evolution we evolved single and double rings of particles. It is actually inconsistent with the SPH scheme to set up a single ring as an initial condition because SPH assumes a smoothed initial state. As would be expected from an SPH simulation, the ring rapidly breaks up into a band. We analyse the stability of the ring and show that the predictions are confirmed by the simulation

Maddison, S T

Monaghan, J J

Murray, J R

CERN Document Server

SPH Simulations of Accretion Disks and
Narrow Rings
Sarah T. Maddison
1
, James R. Murray
2y
and Joe J. Monaghan
1
1
Maths Department, Monash University, Clayton, 3168, Australia
2
CITA, University of Toronto, Ontario, M5S 1A7, Canada
June 6, 1995
Abstract
We model a massless viscous disk using Smoothed Particle Hydro-
dynamics (SPH) and note that it evolves according to the Lynden-Bell
& Pringle theory (1974) until a non-axisymmetric instability develops
at the inner edge of the disk. This instability may have the same ori-
gin as the instability of initially axisymmetric viscous disks discussed
by Lyubarskij et al. (1994). To clarify the evolution we evolved single
and double rings of particles. It is actually inconsistent with the SPH
scheme to set up a single ring as an initial condition because SPH
assumes a smoothed initial state. As would be expected from an SPH
simulation, the ring rapidly breaks up into a band. We analyse the
stability of the ring and show that the predictions are conrmed by
the simulation.
1 Introduction
As a test case for a two dimensional SPH code (for an overview of SPH,
see for example Monaghan 1992) developed to study the viscous evolution
of accretion disks we modelled the spreading of a narrow viscously shearing
ring of matter. We wished to compare simulations with the results of thin

maddison@hypatia.maths.monash.edu.au
y
jmurray@cita.utoronto.ca
1
disk theory as developed by Lynden-Bell and Pringle (1974). Lubow (1991)
used the rate of spread of an axisymmetric ring of viscously interacting SPH
particles to estimate the shear viscosity due to the articial viscosity term
in his code. Flebbe et al. (1994) tested a tensor form of a general Navier
Stokes viscosity term by comparing an axisymmetric ring simulation with
Lynden-Bell and Pringle's result for the evolution of a  function density
distribution. Because SPH cannot be used with a  function density, Flebbe
et al. used an approximate initial condition. We prefer to be consistent
and simulate a conguration with an initial state which can be resolved by
SPH. In particular we looked at an axially symmetric annulus with an initial
surface density

0
(r) =
8
<
:
exp f 

r r
o
l

2
g r
1
< r < r
2
0 elsewhere
(1)
where r
o
and l are the radius of maximumdensity and width of the Gaussian
respectively. The annulus is assumed to be Keplerian with the only forces
considered being the gravitational attraction of the central object and bulk
and shear viscosity forces within the annulus. Pressure and disk self-gravity
are not considered. To simplify the analysis a constant kinematic viscosity
 was assumed. For the above initial condition the general solution for the
surface density at radius r and time t is given by
(r; t) =
1
r
1=4
6t
Z
r
2
r
1
r
0
5=4
exp
8
<
:
 
 
r
0
  r
o
l
!
2
9
=
;
 exp
(
 
(r
2
+ r
0
2
)
12t
)
I
1=4
 
r
0
r
6t
!
dr
0
where I
1=4
(x) is a modied Bessel function of the rst kind. For large argu-
ments, I
1=4
(x)  exp (x)=
p
2x. Thus for small t we have the approximate
solution:
(r; t) 
1
r
3=4
p
12t
Z
r
2
r
1
r
0
3=4
exp
8
<
:
 
 
r
0
  r
o
l
!
2
9
=
;
 exp
(
 
(r
0
  r)
2
12t
)
dr
0
: (2)
Following the equations of Pringle (1981), the simulation is given an initial
radial velocity
v
r
=  
3
R
1=2
@
@R
h
R
1=2
i
where  is our initial Gaussian.
2
0.75 0.8 0.85 0.9 0.95
Radius
0
0.2
0.4
0.6
0.8
1
Su
rf
ac
e 
De
ns
it
y 
  
  
  
 
t=0.0
t=0.2
t=0.4
t=0.8
Figure 1: The time evolution of the surface density of a Gaussian
disk. The solid lines denote the theoretical solution at the times
shown. The heavy points show the corresponding SPH results. For
this simulation we used r
1
= 0:80, r
0
= 0:85, r
2
= 0:90, l = 0:025
and  = 2:5 10
 4
.
The theoretical and simulated spread of an annulus with a Gaussian density
prole is shown in gure 1. The best results were obtained when the parti-
cles were set up in concentric rings so as to give a uniform number density
throughout the annulus. The Gaussian prole was set up by giving each par-
ticle a mass proportional to 
0
(r). In this particular case 22,420 particles in
21 rings were used (see gure caption for further parameters of the model).
After ten revolutions of the particles at r = r
1
the inner edge of the disk
started to break up. Other simulations of disks made up of more rings gave
similar results. After several rotation periods, the viscosity spread the rings
to the point where the innermost ring became somewhat separated from the
remaining rings. It then started to lose shape and break up. This eect
marched its way through the remaining disk. One would like to know if
this is a genuine instability or an artifact of the numerical method. In the
following we give a partial answer by analysing the stability of a single ring.
3
2 Stability Analysis
We consider the stability of a single viscous ring of equi-spaced SPH par-
ticles in orbit with angular velocity 
 around a central mass M . Again,
gravitational and pressure forces of the ring particles have been neglected to
investigate the viscous forces.
Each particle i has unperturbed and perturbed positions in radial and az-
imuthal coordinates:
(r
o
; 
o
) = (a;
t+
i
)
(r; ) = (a(1 + q
i
);
t+
i
+ 
i
)
where a is the ring radius, 
t is a rotating reference line, 
i
is the un-
perturbed azimuthal position, and 
i
and aq
i
are the azimuthal and radial
perturbations respectively.
The radial and azimuthal equations of motion governing a particle in the ring
are:
d
2
r
dt
2
  r
 
d
dt
!
2
=  
GM
r
2
+ F
r
(3)
r
d
2

dt
2
+ 2
dr
dt
d
dt
= F

; (4)
where F = (F
r
; F

) is the viscous force per unit mass.
Using the length and time scales a and t = 
q
a
3
=(GM) we nd:
d
d
= 1 +
d
d
;
and, by linearising equations (3) and (4),
d
2
q
d
2
  3q   2
d
d
=
a
2
GM
F
r
(5)
d
2

d
2
+ 2
dq
d
=
a
2
GM
F

: (6)
For the non-viscous case, where F
r
= F

= 0,
q   3q   2 _ = 0 (7)
 + 2 _q = 0 : (8)
We consider perturbations of the form:
q = A exp fi(k+ ! )g (9)
 = B exp fi(k+ ! )g ; (10)
4
where ! is the frequency of the disturbance and k is the wave number. Sub-
stituting (9) and (10) into the linearised equations (7) and (8), we nd that
! = 0;1, with a double root at ! = 0. The general solution for this non-
viscous case is the superposition of (a) an oscillation, (b) a shift in 
i
without
changing a, and (c) the motion associated with all particles being shifted to
a dierent circular orbit (see Murray 1994).
Returning now to the viscous case, we need the SPH viscous forces. Assuming
that the SPH kernel, W
ij
, is Gaussian (normalised for two dimensions) we
can use:
r
i
W
ij
=  
2r
ij
h
2
W
ij
;
where h is the SPH smoothing length and r
ij
= r
i
 r
j
. Substituting this into
the formula for the viscous force per unit mass (Monaghan 1992, eqn 4.1) and
using 
i
= m
i
=h
2
, where 
i
and m
i
are the density and mass respectively
of particle i, and  is a positive dimensionless constant of order unity that
depends on the spacing of the particles,
F
i
=
X
j
m
j

j
c
ij
h
"
v
ij
:r
ij
r
2
ij
+ 
2
#
r
i
W
ij
=  
2ch

X
j
"
v
ij
:r
ij
r
2
ij
+ 
2
#
r
ij
W
ij
;
where  is the viscosity coecient, c is the speed of sound (replacing c
ij
since the sound speed is constant in the isothermal case) and  prevents
singularities in the viscous force for very small particle separation. We can
ignore  for the purpose of this stability analysis. The kinematic viscosity is
given by  = hc=8.
Since we are considering the linearised equations, we need only consider per-
turbations to v
ij
:r
ij
since it vanishes in the unperturbed case. After lineari-
sation (see appendix) we nd:
v
ij
:r
ij
 ( _
j
  _
i
) sin 
r
2
ij
 2(1  cos ) ;
where  = 
j
 
i
. The scaled radial component of the force per unit mass
on particle i is then given by:
a
2
GM
F
r
=  
g
2
X
j
( _
j
  _
i
) sin  W
ij
; (11)
and the tangential force per unit mass component on particle i is:
a
2
GM
F

=
g
2
X
j
( _
j
  _
i
)(1 + cos )W
ij
; (12)
5
where g = 2ch=.
We now have the full SPH formulation of the equation of motion of a viscous
ring, by combining (5) & (6) and (11) & (12), and substituting (9) and (10).
The summation is over all N particles in the ring. However, only nearby
particles contribute to the sum, which can thus be replaced by an innite
sum. It is then clear that the odd components of these sums must vanish,
leaving the equation of motion as:
( !
2
  3)A  2i!B = gB!	 (13)
 !
2
B + 2i!A =  gBi! ; (14)
where
	 =
1
2
X
j
sin(k) sin  W
ij
 =
1
2
X
j
(1  cos(k))(1 + cos )W
ij
;
and if k = 0 both 	 = 0 and  = 0, thus F
r
= 0 and F

= 0.
Combining (13) & (14) we get a quartic in the frequency !, one root of which
is ! = 0. This leaves us with:
!
3
  ig!
2
+ (2ig	   1)!   3ig = 0 : (15)
In the non-viscous case (i.e. where 	 =  = 0) ! = 0;1, in agreement with
the solution of (7) and (8).
If the viscosity (which is proportional to g) is small enough, we expect fre-
quencies close to the unperturbed values. We therefore set ! = ; 1 +  and
 1 +  where jj  1 and the  in each case will be dierent.
! = 1 + 
Substituting ! = 1+  into equation (15) and keeping only the linear terms
(since jj  1) we nd:
 = ig(2  	) ; (16)
and
expfi!tg = expfi(1 + )tg
= expfit  gt(2 	)g :
For stability we require 2  	  0. Substituting for  and 	 we nd that
for all values of k; 2   	  0 (see Murray 1994 for this and the following
two cases). Thus the ! = 1 +  mode is stable.
6
! =    1
In this case we get
 = ig(2 + 	) (17)
and so expfi!tg shows that for stability we require 2+	  0, which holds
for all values of k. Hence ! =    1 is also stable.
! = 
Here we get:
 =  3ig : (18)
This gives expfi!tg = expf3gtg whose exponent is always positive, and
hence the ! =  is an unstable mode for all values of k. Note that in the
inviscid solution, the ! = 0 case is associated with an azimuthal shift with
radius held constant.
For large wavelengths (large compared to the particle spacing, which is of
order h, and thus large wavelength corresponds to hk  1), we substitute a
Gaussian for the kernel and approximate the summation by an integral and
nd:
 
k
2
h
p

; (19)
where  is the azimuthal particle spacing. Thus:
! =  =  3ig
  
6ic

p

(hk)
2
: (20)
Since we are assuming that jj  1, we must have:
c

p

 1 ;
to be consistent with the analysis. In the models that we ran, particles were
equally spaced a scaled distance   2=200 apart, with c  1 (see
below for details). Hence, we have found that there is an unstable mode
when !  1, ie when c 1.
The time scale of the viscous instability for long wavelengths is:




p

6c
:
In the models we ran,  = 0:07, c = 0:05 and   4=9. Thus 

 1:178 rota-
tion periods. However, from (20) we see that the most rapid growth occurs
for short wavelengths for which the approximation (19) to the summation
is invalid. Hence we expect instabilities with wavelength approximately a
particle spacing to grow fastest.
7
Figure 2: The onset of viscous instability in a single ring system.
The gures are at 0:2; 4:1 and 5:8 rotational periods.
3 The Simulations
We ran an SPH simulation of a single viscous ring of approximately 200
equally spaced particles orbiting a central massive body. The ring was placed
at dimensionless radius of 0:9, at which position the ring has a rotational
period of 5:4 dimensionless time units. As mentioned above,  = 0:07 and
c = 0:05. Particles started to oscillate very slightly almost immediately and
after approximately four rotational periods the ring had broadened. The
particles moved on eccentric orbits and after ve rotational periods the ring
had broadened signicantly (see gure 2). It was seen from the simulation
that the viscous instability sets in after approximately four rotation periods.
A double ring system was also modelled. If the separation of the rings was
less than 2h (the radius of inuence of the SPH kernel), the two rings acted
as one and were disturbed in a similar manner as was the one ring system.
However, if the separation of the rings was greater than 2h, then as usual
viscosity acts to spread the rings and particles moved onto eccentric orbits.
The inner rings broke up rst after four rotation periods, followed by the outer
ring at its equivalent rotational phase. The two rings were indistinguishable
after eight orbits.
4 Conclusions
The results of the disk simulations show that SPH gives results in good
agreement with the standard viscous theory. Deviations from theory eventu-
ally occur at the inner edge where we observe particles moving from circular
8
orbits. The analysis of Lyubarskij et al. (1994) shows that axisymmetric
viscous disks can become unstable whenever the viscous drag on uid in an
elliptical orbit is a maximum at apastron. This is the case for the standard
-disk and it is the case for particles in our simulation. For example, the
particles in the innermost ring appear to move in and then out which shows
that they are individually moving on elliptical orbits with apastron near their
neighbouring ring. Since the viscous force is felt when they move near the
particles in the neighbouring ring they will be forced into more eccentric
orbits.
When modelling viscous disks we found that the innermost ring separated
from the remaining disk and then became unstable. We then investigated
the dynamics of a single ring. However, the formulation of SPH is based on
smoothing. A single ring is the antithesis of a smoothed conguration. The
single ring smooths itself by rapidly spreading under the viscous forces. The
combination of analysis and simulation conrms the consistency of SPH.
However, the resulting eects of the SPH viscous force prescription in ac-
cretion disk simulations is unphysical. The main problem is that we are
trying to model a continuum (i.e. a gaseous viscous disk) by a set of discrete
rings. For disks made up of concentric rings, this problem can be overcome
if the number of rings used is increased so that the rings will not separate
due to viscosity by too great a distance (generally less than 2h). As usual
with particle methods, a compromise must be reached between the number
of particles and the computation time. Alternatively, one could incorpo-
rate a variable smoothing length, so that when the distance between rings is
greater than 2h, the smoothing length increases so that neighbouring rings
can still communicate. There are still some problems to be overcomewith the
implementation of spatially varying smoothing lengths when sharp density
gradients are present (see Nelson & Papaloizou 1994).
9
A Appendix
We want to solve for
s = v
ij
 r
ij
= (v
i
  v
j
)  (r
i
  r
j
): (A1)
Noting that to rst order
r
i
= (1 + q
i
)
^
r
io
+ 
i
(
^
z
^
r
io
)
and dierentiating with respect to t, (noting d
^
r
io
=dt =
^
z
^
r
io
, with 
 = 1),
we nd
v
i
= ( _q
i
  
i
)
^
r
io
+ (1 + _
i
+ q
i
)(
^
z
^
r
io
); (A2)
where ^denotes a unit vector, r
io
is the unperturbed position vector of particle
i, and
^
z is perpendicular to the plane of the orbit. We can write the terms
contributing to (A1) as, for example,
v
i
 r
i
= ( _q
i
  
i
)(1 + q
i
) + 
i
( _
i
+ 1 + q
i
)
' _q
i
v
j
 r
i
= [( _q
j
  
j
)(1 + q
i
) + (1 + q
j
+ _
j
)
i
] cos 
' ( _q
j
+ 
i
  
j
) cos    (1 + q
i
+ q
j
+ _
j
) sin : (A3)
where  = 
j
  
i
and
^
r
io

^
r
jo
= cos . Similar expressions can be written
down for v
j
 r
j
and v
i
 r
j
. Combining these expressions we nd
s = ( _q
i
+ _q
j
)(1   cos ) + ( _
j
  _
i
) sin : (A4)
Since 0 <   1 we can neglect the rst term in (A4) to obtain
s ' ( _
j
  _
i
) sin :
10
References
Flebbe,O., Munzel,S., Herold,H., Riert,H. & Ruder,H. 1994, ApJ, 431, 754
Lubow,S.H. 1991, ApJ, 381, 268
Lynden-Bell,D. & Pringle,J.E. 1974, MNRAS, 168, 603
Lyubarskij,Y.E, Postnov,K.A. & Prokhorov,M.E. 1994, MNRAS, 266, 583
Monaghan,J.J. 1992, ARA&A, 30, 543
Murray,J.R. 1994, PhD Thesis, Monash University
Nelson,R.P. & Papaloizou,J.C.B. 1994, MNRAS, 270, 1
Pringle,J.E. 1981, ARA&A, 19, 137
11


SPH simulations of accretion disks and narrow rings

Abstract

Similar works

Full text

Available Versions

CERN Document Server