Shocked Clump (2D)

Am investigating the interaction between a cold clump in a supernovae remnant and a shock. To model this, am injecting a supersonic wind into the boundary, whose properties are specified by the jump conditions across an adiabatic shock. The jump conditions are calculated using the following parameters in the ambient medium (note all values given are approximate; actual values are calculated to machine precision in the attached problem module and data files for the fiducial case described here):

T_{amb}	10^4 ~K
n_{amb}	10^4 ~cm^{-3}
v_{amb}	0

The sound speed of the ambient material for these parameters is (for \gamma=5/3):

C_{amb}

10 km~s^{-1}

which for a injected wind speed of v_s=2,000 ~km/s gives a Mach number of:

M_{wind}

194

The wind parameters are then given by the R-H jump conditions (see attached pdf for equations used):

T_{wind}	1e+8~ K
n_{wind}	4e+4 ~cm^{-3}
v_{wind}	1.5e+8 ~cm~s^{-1}

The clump is initially in pressure equilibrium with the ambient medium, and its properties are given by:

T_{clump}	10^2 ~K
n_{clump}	10^6 ~cm^{-3}
v_{clump}	0

The thermodynamics in the simulation is governed by an ideal EOS (with \gamma=5/3). The abundances are taken to be solar (\mu=1.3). Other parameters of interest are the clump radius (R_c), clump-ambient density contrast (\chi), clump crushing time (t_{cc}=\sqrt{\chi} R_c/v_s), sound crossing time (for sound waves in ambient and clump medium to travel a clump radius; t_{sc}), wind crossing time (for wind to travel clump radius; t_{wc}):

R_{c}	10^{16} ~cm
{\chi}	100
t_{cc}	5\times 10^{8}~s
t_{sc} (ambient medium)	10^{10}~s
t_{sc} (clump medium)	10^{11}~s
t_{wc}	5\times10^{7}~s

The domain dimensions, resolution, clump origin, simulation time & framerate (for 155 frames), are given by:

L_{x} = L_{y}	10~R_c
dx	1/64~R_c
O_{clump}	2~R_c,~ 5~R_c
t_{sim}	t_{wc}~ +~3~t_{cc}= 155\times 10^7~s
dt	10^7~s~=~1/50~t_{cc}

Run 1 - Steady inflow condition

Run description: supersonic wind is injected into the left edge of the computational domain each timestep over the course of the simulation.

Note that we must initialize the left and right state of a "Riemann problem" (where the Riemann interface lies along the boundary between the grid and ghost zones) such that the generated wave fan consists of a shock traveling to the right into the grid at 2,000 km/s. In other words, we initialize the boundary zones (the "post-shock" gas) with the values specified by the jump conditions across a 2,000 km/s shock. The pre-shock values correspond to the fluid variables in the grid. If the incoming flow does not match the jump condition values exactly, additional waves will be generated once the incoming flow collides with the ambient gas in the box. To illustrate this, see the following image which shows an initial inflow at the boundary that has slightly different values than those specified by the jump conditions (left), compared to the true values (right).

The shock front in the right panel shows some oscillations, but these are expected given the numerical scheme. Currently, am using a PPL interpolation method, but will see if PPM reduces these oscillations further in a future run.

Visualizations:

(note time slider is in units of the crushing time)

n (cm-2)	movie​
T (K)	movie​
P (Ba)	movie​
Mach	movie​
Schleiren	movie​

Run 2 - PPM vs. PPL scheme

This run was just to see if the choice of solver can significantly change the simulation. To check this, here is a line plot of the density vs. position for the PPM method (the previous simulation used PPL):

and a comparison of the final frame:

Run 3 - Pulsed case (supersonic wind only injected at t=0)

I tried to run this case by removing the call to initialize a wind in the beforestep routine in problem.f90, but this produced identical results to my continuous case. This leads me to believe that the continuous case does not actually need a call to initialize a wind every time step, but rather, initializing the wind in the beginning of the simulation produces an inflow boundary continuously over the course of the simulation.

If we want to look at a pulsed case then, we need to turn off the wind after some dt and specify a new boundary condition along the inflow edge. A continuous post-shock wind may be more well-suited for the problem at hand, however, rather than a pulsed inflow boundary condition. Thus, investigating this case further is currently low-priority.

Run 4 - Mach Comparison (Mach 1.5, 10, 100)

These simulations were intended to explore the dependance of bow shock thickness (the distance between the edge of the bow shock and the clump surface) on the mach number of the impinging flow. For these runs, I moved the boundaries far enough away so that the bow shock did not reach them over the course of the simulation. The effective resolution of the clump was kept constant using AMR, but the resolution in the ambient medium was reduced to speed up the calculation. Given the ambient medium is uniform and stationary, reducing its resolution in this way is reasonable, so long as the shock front is resolved to the same effective resolution as in the fiducial case (which produced a well-converged solution as shown in earlier line plots of the front).

The first set of plots shows a comparison of the bow shocks (for M=1.5, 10, 100) at the moment the forward shock wraps around the clump and converges on itself, thereby creating a ‘mach stem’ or stagnation point behind the clump. As these plots show, the bow shock thickness is roughly equal in the different cases.

In viewing an animation of these different runs, the bow shock of the M=1.5 run seemed to still be evolving quite rapidly at this point, in that it was expanding away from the clump at a much higher rate than the other runs. This made me wonder if a better comparison of the bow shock thickness would be when (if) the bow shocks came to some steady-state configuration with respect to the upstream flow. The next set of plots shows these 3 runs at a later time, at which it seems the M=10, and 100 have reached a near steady state. The M=1.5 seems to still be expanding by this point.

Some reference papers (see below) for shock/clump interactions provide theoretical estimates for the thickness of the bow shock. The results presented here should be explored further to see if they align with these estimates. Furthermore, while the time scale on which the clump/bow shock boundary becomes unstable to the Kelvin-Helmholtz, Rayleigh-Taylor, and nonlinear thin shell instabilities is roughly a clump crushing time, I expect the wavelength of the resulting instabilities will also be a function of the upstream mach number.

Thus, varying the mach number of the upstream flow can answer two questions:

1. What is the dependance of bow shock ‘thickness’ (or, similarly the time scale to reach a steady state configuration) on the mach number of the upstream flow?

2. what is the relationship between upstream mach and the wave length of fluid instabilities excited along the shocked clump surface?

Run 5 - Mach Comparison (Mach 100 vs. 200 vs. 300)

Produced similar results as the previous set of mach runs.

Run 6 - AMR tests

We want to track the shock front to the highest resolution possible to minimize numerical dissipation of the shock front as it traverses the grid.To achieve this, we need the AMR to refine the shock front to the max level possible. Additionally, we want the clump to be resolved adequately; as a rule of thumb the clump should have at least 64 zones per clump radius.

To track the shock front, we could use refinement criterion based on density gradients as used here. As can be seen from these images, the shock front is refined to the maximum AMR level over the course of the simulation. The clump, however, is being refined to the max level just along its edge, which exhibits a sharp density gradient. Forcing the entire clump to be refined to the max level can be done upon initialization easily enough using a geometrical refinement object in the code.

In testing AMR on this setup, I was finding that the meshing algorithm produces slightly different AMR patch distributions when different numbers of processors are used (for reference, the run presented here used 16 processors). It seems that numerical errors are produced within the clump when the mesh is asymmetric, resulting in asymmetries in the fluid instabilities (compare left image with AMR to right without AMR):

This is particularly noticeable along the central axis of the clump, normal to the shock. More on this in the conclusions section below.

Run 7 - Prelim 3D run

Estimate ~ 5 hours per frame with cooling turned on at same resolution above (on 24 processors).

Conclusions

The above results show that the structure of post-shock instabilities that form along the surface of the clump can be sensitive to the choice of solver as well as the alignment of the mesh with respect to the clump center. This suggests that we should break the symmetry in the initial conditions to avoid (unphysical) behavior based on numerical error. One way to do this would be to introduce a sinusoidal perturbation to the shock front.

In addition, we would like to experiment with the internal density distribution of the clump, e.g. using a Gaussian density profile instead of a uniform distribution, as well as remove the pressure equilibrium constraint on the clump/inter-clump medium, and add a cooling prescription to the simulations. Lastly, there is no reason to not initiate the supersonic flow within the computational domain as well as along the boundary. Initializing the shock just a few cells away from the clump will minimize computation time as well as numerical errors as the shock propagates along the grid.

Unresolved thoughts:

Q - Reflected BC differences - in only simulating the upper half of the domain, asymmetries formed in the shocked clump layer similar to those shown above in the AMR and PPM cases. Would be interesting to investigate this further — e.g. what are numerical instabilities that occur along the symmetry axis in jet models when reflecting BC's are used, and could those be applicable here?

Q - Size scales needing to resolve (i.e. instabilities and their wavelengths) for supernovae remnants and their embedded shocked clumps?

Q - Examine cooling cases (cooling lengths/timescales, temperature ranges for various cooling tables available in astrobear, modifying the abundances upon initialization, as well as dynamically throughout the simulation, what options are there for tracking ionization fractions, what additional instabilities can arise with cooling turned on.

Q - Examine radiative shock structures in steady state solutions

Reference Papers

See my supernova remnant library on ads here.