BackLinksMenu

1D Pulsed Jet Simulations

The set up for these simulations is very simple. The grid is initialized with an ambient object and a wind object. The goal is to study the effect of magnetic fields on the emission behind shocks in a pulsed jet.

Resolution Study

Simulations involving cooling tend to require a lot of resolution. For this reason, I have been running a simplified version of the 1D Jet simulation to see how resolution affects the dynamics. Here are the initial parameters:

Ambient:
density = 1000 particles/cc
velocity = 0 km/s
temperature = 1000 K
By = 0 uG

Wind:
density = 10000 particles/cc
velocity = 50 km/s
temperature = 1000 K
By = 0 uG

Global Data:
Domain length = 30 AU
Final time = 1.7 computational units (approx 7 cooling times)

Note that these runs feature a stationary ambient, and they do not include MHD or jet perturbations. To get an idea of the form of the shock waves that are produced, here is a lineout of velocity:

This is from a high resolution run (100 base cells + 7 levels of AMR = 12800 effective). You can see that even with this resolution, the front shock is somewhat smooth instead of a sharp jump. Now in order to make some predictions based on the initial parameters, we need some equations…

Important Equations

A strong shock implies that the mach number M >> 1. We will also be assuming an ideal gas, so \gamma = \frac{5}{3}. The following equations can be derived from the usual shock jump conditions:

T_{ps} = \frac{3m_{H}v_{s}^{2}}{16k} \hspace{1 in} (1)

v_{s} = \frac{4}{3}v_{ps} \hspace{1in} (2)

n_{ps} = 4n_{a} \hspace{1in} (3)

Where s stands for shock, ps for post-shock, and a for ambient. T is temperature, v is velocity, n is number density, mH is the mass of a hydrogen atom, and k is Boltzmann's constant.

In order to figure out these post-shock values, we need to somehow relate them to the jet. This is where these equations come in:

v_{ps} = v_{j}(1 + \eta^{-\frac{1}{2}})^{-1} \hspace{1in} (4)

\eta = \frac{n_{j}}{n_{a}} \hspace{1in} (5)

Where j stands for jet, and \eta is basically the density contrast of the jet material to the ambient material. Now we can combine these equations and constants to come up with an equation for the post-shock temperature in terms of jet parameters:

T_{ps} \approx (40.41 \mathrm{\frac{K s^{2}}{km^{2}}}) v_{j}^{2} (1+\eta^{-\frac{1}{2}})^{-2} \hspace{1in} (6)

Where vj is more conveniently in units of km/s. Now we need equations for the characteristic cooling time and cooling length:

t_{cool} = \frac{3kT_{ps}}{n_{ps}\Lambda} \hspace{1in} (7)

l_{cool} = v_{ps} t_{cool} \hspace{1in} (8)

Where \Lambda is the cooling rate in units of cm³erg/s. Using equations (3) and (4), these cooling equations can be rewritten in terms of our known initial parameters, and the post-shock temperature which we will want to have calculated already:

t_{cool} = \frac{3kT_{ps}}{4n_{a}\Lambda} \hspace{1in} (9)

l_{cool} = v_{j}(1 + \eta^{-\frac{1}{2}})^{-1} \frac{3kT_{ps}}{4n_{a}\Lambda} \hspace{1in} (10)

Some Predictions

Based on all these equations, we can make some predictions about the post-shock temperature and what resolution we might want. Using the aforementioned initial parameters, you get the following predictions:

Tps (103 K)	tcool (107 s)	lcool (AU)
58.3123	1.26155	3.20345

A rate of 4.7863x10^-22 was used, which is fairly consistent with Dalgarno McCray. So if we have a domain that is 30 AU and we want 10 cells per cooling length, then we would need a resolution of about 94 cells. That doesn't seem too bad, so I will round up to 100 and start from there.

Results

Here we look at what post-shock temperature the simulation gives as a function of the effective resolution. The effective resolution = (# cells) x 2^{(# AMR levels)}. Each run used a base grid of 100 cells, so the only thing that changed from run to run was the number of AMR levels. For some reason, I could not get it to run with more than 7 levels, so the 25600 run is actually a base grid of 200 cells with 7 levels. For Tps I will use the highest temperature that I get from the simulation. The highest temperature always occurs early on in the simulation and then the post-shock temperature slowly decreases over time, but stays relatively constant at some lower temperature.
Z

Effective Resolution	cells/lcool	Tps (103 K)	Relative Error (%)
100	10.68	34.37	41.06
200	21.36	37.11	36.36
400	42.71	40.21	31.04
800	85.43	43.68	25.09
1600	170.85	46.66	19.98
3200	341.70	49.25	15.54
6400	683.40	51.18	12.23
12800	1366.81	53.20	8.77
25600	2733.61	53.66	7.98

And in graphical form, the table looks like this where I have plotted the expected Tps as the green horizontal line:

So the big question is, why does it appear that we need way more resolution than is predicted from the cooling length? Perhaps it has to do with how the refinement is being triggered. So I will try a resolution of 400, but with no AMR and see if it does any better than the 400 effective resolution which used 2 levels of AMR. With this run I got Tps = 40.59 x 10³ K. So there really was not a significant difference from AMR to fixed grid.

My next idea was to revert back to just DM cooling instead of the new Z cooling. Perhaps, the Z cooling routine is not implemented quite right, and there is some "double counting" for cooling strength or something like this. So here is a data table where I used DM instead of Z. The improvement column is the decrease in relative error from the comparable Z cooling run.
DM

Effective Resolution	cells/lcool	Tps (103 K)	Relative Error (%)
100	10.68	25.71	55.91
200	21.36	31.63	45.76
400	42.71	41.07	29.57
800	85.43	47.65	18.28
1600	170.85	51.37	11.91
3200	341.70	53.30	8.60
6400	683.40	53.98	7.43
12800	1366.81	54.14	7.16
25600	2733.61	58.16	0.26

The DM cooling gets closer to what I expect but the Tps is still a little low. In both cases, Tps seems to be converging to some other value different from my expected value of 58.3123. However, when I do a run with no cooling, I get post-shock temperatures reaching 58.3123, and some even higher…closer to 60.97 with just 100 cells and no AMR. So I decided to take a closer look at the adiabatic case to see if I could figure out the discrepancy. Again, here is the data for no cooling:
None

Effective Resolution	cells/lcool	Tps (103 K)	Relative Error (%)
100	10.68	60.97	4.56
200	21.36	61.09	4.76
400	42.71	61.00	4.61
800	85.43	61.06	4.71
1600	170.85	61.76	5.91
3200	341.70	63.45	8.81
6400	683.40	64.17	10.04
12800	1366.81	65.20	11.81
25600	2733.61	64.16	10.02

You can see that these simulations actually give higher than expected post-shock temperatures. If I used these values as the "correct" values then the post-shock temperatures from the cooling simulations are even further off. Here is a plot that contains all the data:

Shock Structure

I've made some lineouts of density, velocity, pressure, and temperature to get a better understanding of how the simulation is interpreting the shock structure. Every image is taken at the middle of the simulation for the highest resolution run. The solid black line is density, dashed blue is velociity, dotted green is pressure, and dash/dot red is temperature. Everything is in computational units, and it's plotted on a log scale. The structure is more difficult to see in the cooling cases, so I included some magnified images for those two cases.

No Cooling

DM Cooling

Z Cooling

Velocity Movie (sine wave perturbations, no cooling, no B-field)​