[FLASH-USERS] Equilibrium configuration

[Jesús Zavala Franco]

Dear all,

I'm having problems setting up a sphere of gas in hydrostatic equilibrium
under its own gravity + the gravity of a particle distribution.

In the general case, I would like to have this particle distribution with
the same density profile as the gas distribution (so I'm aiming at setting
up a dark matter halo with a particle distribution, with gas inside), but to
describe the problem I'm having, I will avoid distributing particles in the
sphere and simply:

1) put 1000 particles in the centre of the sphere
2) put gas distributed spherical around this centre with a radial density
profile (a NFW profile), and pressure given by the condition of hydrostatic
equilibrium:

dP(r)/dr = -(GM(<r)/r^2) * rho_gas(r)

where M(<r) is the total enclosed mass, and rho_gas(r) is the gas density.
M(<r)= M_DM+Mgas(<r), with M_DM the total mass of the 1000 particles, which
by the way exceeds the total mass of the gas by almost an order of
magnitude. To solve the equation I impose a boundary condition outside the
sphere setting a pressure and a density which are reasonable according to
the problem and the density profile I'm putting.

So in short, what I have is a sphere of gas within a gravitational potential
dominated by essentially a point source in the centre and with a pressure
that should give support agains the gravitational collapse.

I'm using a gird of 8^3 with 4 levels of refinement, an ideal gas equation
of state, a Pfft Multigrid gravity solver, I'm using the default operator
splitting technique to advance the solution.

What I notice is that since the first time step, regardless of my choice of
dtinit, the cells just next to the central cell (the one containing the 1000
particles), acquire a radial outwards velocity, whit a size depending of the
time step, and that once the time evolution reaches the typical times of the
problem (~100 Myrs), this results in the gas in the inner parts propagating
outwards, reducing the density in the core, after 1Gyr or so, this
propagating gas reaches the boundaries of the sphere. In other words, the
equilibrium set at first is broken.

I have tried using the option ppm_modifystates=.true. since as described in
the user guide:

"The version of PPM in the FLASH code has an option to more closely couple
the hydrodynamic solver

with a gravitational source term. This can noticeably reduce spurious
velocities caused by the operator

splitting of the gravitational acceleration from the hydrodynamics"

I thought this could help, but it didn't. I'm aware that if I wouldn't put
the particles in the centre, this behaviour is just a resolution problem
since the core is not being properly resolved and the density and enclosed
mass is underestimated there, the pressure is therefore too high (since is
initially set by solving the hydrostatic equilibrium equation with the
assumed analytical profile) and it pushed the gas outwards. However, since
I'm setting a significant gravity source in the center I wouldn't expect
this to be the case. I have tried actually increasing the mass of the
particles by a factor of 10, and not considering this mass increase in the
solution to the pressure equation, this sets the overall pressure too low
and you would expect a collapse. Even though when at first the velocities of
the contiguous cells to the centre, point inwards, after a while, they are
overwhelmed by an outward flow that develops in the cells next to these, and
to my surprise, the sphere expands as well.

I'm running out of ideas on the initial conditions setup, so I'm thinking
this could either be a resolution issue, or a bad choice of hydro, gravity
solver.

Any hel will be much appreciated.

Cheers,
Jesus Zavala
Department of Physics and Astronomy
University of Waterloo
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://flash.rochester.edu/pipermail/flash-users/attachments/20110317/428a7f0b/attachment.htm>