1D Radiative Shocks

These simulations follow Ch. 4 of Delamarter '01 and his treatment of the 1D steady radiative shock problem. The purpose of this problem is to check that the cooling source terms are being handled correctly. Here are some useful links:

​AstroBEAR 1.0 test page

​Delamarter '01 pdf

Equations

The shock jump equations for a stationary shock are used to solve for the initial post-shock values. The post-shock velocity v₂ can be written as:

v_2 = v_1 \frac{(\gamma - 1)M^2 + 2}{(\gamma + 1)M^2}

where v₁ is the ambient velocity, and M is the ambient mach number. Remember that the mach number M = v₁/c where c is the ambient sound speed, and c = \sqrt{\frac{\gamma P_{1}}{\rho_{1}}} where P₁ is ambient pressure, and \rho_{1} is ambient density. The post-shock density and pressure (\rho_{2} and P₂ respectively) can be found by using mass flux and momentum flux conservation across the shock:

\rho_1 v_1 = \rho_2 v_2

\rho_1 v_1 ^2 + P_1 = \rho_2 v_2 ^2 + P_2

These post-shock values become the boundary conditions for the fluid equations in the cooling region:

\rho v = \rho_2 v_2

\rho v ^2 + P = \rho_2 v_2 ^2 + P_2

\frac{d}{dx}(\frac{1}{2}\rho v ^3 + \frac{\gamma}{\gamma - 1} P v) = -\Lambda

where \rho, v, and P represent the density, velocity, and pressure in the cooling region as functions of x, and \Lambda is the cooling rate.

Initial Parameters

n₁ = 60 particles/cc

v₁ = 10⁷ cm/s

T₁ = 10⁴ K

For analytic cooling,

\Lambda = n^{2} \alpha T^{\beta}

where n and T and the number density and temperature in the post-shock region respectively.

For these simulations,

\beta = 2

\alpha = 1.23786 * 10^-34 erg*cm³/s/K²

cell length = 2.5 x 10¹⁵ cm

problem domain = 400 cells

final time ~ 4000 years

Results

A run with no cooling was done to make sure the post-shock values were correct. With no cooling (adiabatic), the hydrodynamic quantities should jump discontinuously at the shock, and then remain constant while the shock itself remains stationary. Only the first 100 cells of the simulation are shown here.

Adiabatic
Pressure Movie​	Temperature Movie​
Density Movie​	Velocity Movie​

Then, a run with the aforementioned analytic cooling parameters was done. The hydrodynamic quantities will change discontinuously at the shock as before, but then change continuously according to the cooling rate. The ambient temperature is set to be the floor temperature, so that no cooling occurs at or below this temperature. Once the post-shock temperature reaches this floor temperature, the hydrodynamic quantities will remain constant. This is called the quiescent region. Only the first 100 cells of the simulation are shown here.

Analytic Cooling
Pressure Movie​	Temperature Movie​
Density Movie​	Velocity Movie​

For Dalgarno-McCray (DM) cooling, some of the initial parameters had to be changed. The initial velocity used in the previous simulations would give rise to a post-shock temperature in the unstable region of the DM curve. This led to the radiative instability, and thus unsteady shocks. To avoid this, the initial velocity was lowered to 80 km/s. Also, the cooling rate from the DM curve is much higher than that of the analytic form. In order to see the cooling region, the grid had to be shrunk, resulting in a cell length of 10¹³ cm.

DM Cooling
Temperature Movie​

For Non-equilibrium (NEQ) cooling, the length scale had to be adjusted. The shock velocity is the same as the DM case, but the cell length is now 7.75e11 cm.

NEQ Cooling
Temperature Movie

Z Cooling has the same initial parameters and scaling as the NEQ case.

Z Cooling
Temperature Movie

Here is a fun plot that shows the temperature profiles for DM, NEQ, and Z cooling. This gives you a good idea of the temperature ranges that each type of cooling works in. Also, you can see the differences in the cooling length which is why I would normally change the length scale depending on the type of cooling.

Including MHD

Making these simulations work with MHD was no simple task. Having a magnetic field adds another jump condition, and it makes the equations within the cooling region more complicated. Here are what the equations look like now: