HydroStatic Star

Purpose

This module aims to reproduce a situation of hydrostatic equilibrium (HSE), a condition where a volume of a fluid is at rest or at constant velocity. This occurs when compression due to gravity is balanced by a pressure gradient force:

we can translate the image above into the condition:

frac{dP}{dh} + \rho(r) \cdot g(r) = 0\\

where dP/dh is the change in pressure with respect to height, ρ(r) and g(r) are respectively density and gravitational potential.

Implementation

First we determine a 1D profile in hydrostatic equilibrium, then the profile is ported to 2D/3D using interpolation.

Objects

This module makes use of 3 objects:

• Ambients
• Profiles
• PointGravity

Process

We can think of a 1D column of gas as an array:

COLUMN TOP

h5	ρ5	P5
h4	ρ4	P4
h3	ρ3	P3
h2	ρ2	P2
h1	ρ1	P1
h0	ρ0	P0

COLUMN BASE

We can relate each level of the column to an position in an array which contains height (h), density(ρ) and pressure(P), that is:

REAL(KIND=qPREC), DIMENSION(100,3) :: column !where 100 is the number of entries 
                                             !3 is the number of attributes for each entry
                                             !(1 height or radius, 2 rho, 3 pressure) in this case

This array need to be populated with a density profile, namely, the fields 1(height or radius) and 2(rho) for all 100 entries above. The module will be able to calculate the ideal pressure in order to create an HSE solution.

To calculate a pressure profile we need to interpret:

frac{dP}{dh} + \rho(r) \cdot g(r) = 0\\

in a discrete way, this gives:

(h) = P(h+\Delta h) + g(h)\cdot\rho(h)\cdot\Delta h

The condition above determines the pressure profile starting from the top of the column until its base. This implies we specify the pressure at the highest distance from the point mass, ideally infinity. This value can be adjusted in problem.data :

&ProblemData
mass=0.4d0         !PointGravity Mass
soft_radius=1d0    !Softening radius (in computational units)
soft_function=2    !Softening function (NOSOFT=0, SPLINESOFT=1, PLUMMERSOFT=2)
velocity=0,0,0     !PointGravity velocity
xloc=0,0,0         !PointGravity location
pressure_inf=1d0   !Ambient Pressure at infinity (cu)
/

You can check the Profiles page to familiarize with the code that handles these calculations.
This model assumes that the effects of self gravity are relatively small compared to the gravitational attraction of a point mass.

2D, 3D and Softening

Due to steep changes in density and gravitational potential close to the center of mass of a star, the point gravity in this model, softening is needed.
With PlummerSoftening being used to soften the gravitational potential and force in this module, a difference in the pressure and temperature profile arise between 2D and 3D models.
Below is how gravitational acceleration is calculated for 2D and 3D respectively :

2D:

GravityAcceleration=-PlummerSoftening2D*mass*ScaleGrav
PlummerSoftening2D=2d0*pos/(sum(pos(1:nDim)**2)+r_soft**2)

3D:

GravityAcceleration=-PlummerSoftening*mass*ScaleGrav
PlummerSoftening=pos/((sum(pos**2)+r_soft**2)**(1.5d0))

Clearly, the 3D model shows a steeper gradient in pressure, this has important repercussions on the temperature profile as well.
The plot below indicates how the pressure profile changes between 2D and 3D models.

Refence model comparison and results

Testing can be done in two ways: 1) by comparing model properties to analytic solutions 2) by running the model

For the first testing stage, the density profile of an AGB star (can be found ​here) was loaded into the module, after calculating pressure and temperature profiles needed for an HSE condition, the generated profiles were compared with the reference model.
The following plot below compares the two models:

In the second testing stage, the model was run to check for changes in velocity.
Clearly, the outcome that can be expected from an equilibrium condition is the absence of change in velocity. In the case of an equilibrium at rest, we would expect our model not to change over time.
Below is the density, pressure, temperature and mach profiles for our model.