p1
p5
ρ1
ρ5
u1 = 0 u5 = 0
Head RF Tail RF
CD Shock
Length: x≤0.5xmax x >0.5xmax
Velocity: u1 = 0 u5 = 0
Density: ρ1 = 1 ρ5 = 0.125
Pressure: p1 = 1 p5 = 0.1
Thermal energy: Etherm,1 = 1.5 Etherm,5 = 0.15
Figure 4.3: Sod shock tube test. This sketch is a variant of Fig. 3.7. The top panel and the table show theinitial conditionsof the Sod shock tube. Pressures are shown in green, densities in blue and velocities in red. The dashed lines in the lower panel separate the five zones of the Sod similarity solution: (1) the unperturbed state of the denser medium, (2) the rarefaction fan, RF (3) the contact discontinuity, CD(4) the fast shock wave and (5) the unperturbed state of the tenuous medium.
4.2 Sod shock tube test 61
5. unperturbed state of tenuous, low pressure medium (Etherm,5, p5, ρ5, u5)
The states of the gas in zone 1 and zone 5 are known from theinitial conditions. Also all thermal energies can be found from the EOS (Eq.3.2) which relates the thermal energies to the adiabatic exponents, densities and pressures. Hence there are nine unknowns: p2, ρ2, u2, p3, ρ3, u3, p4, ρ4
andu4. Sect.4.1.2 tells us that pressure, velocity and density in the rarefaction wave (p2, ρ2 and u2) are definite if the states 1 and 3 are known. Basically the rarefaction is a reversible, adiabatic process and the Riemann invariants lead to the solution. As we saw in Sect.4.1.1, another unknown speed and pressure can be removed, since there is no mass flow through the contact discontinuity and the pressure is continuous at thecontact discontinuity. Therefore we define new quantities at thecontact discontinuity: a velocityu3 = u4 =: ucand a pressurep3 = p4 =:pc. Moreover the constant entropy in the rarefaction wave allows us to connect the state 3 to the state 1 via
pc=p1
ρ3
ρ1
γ
. (4.14)
From the other Riemann invariant, the sound speed, we find:
uc+ 2 γ−1
rγpc
ρ3
= 2
γ−1 rγp1
ρ1
. (4.15)
Thus, we are left with three unknowns: pc,ρ4anduc. The post-shock medium (state 4) is separated from the pre-shock medium (state 5) by a discontinuity. We can connect these two states with the Rankine Hugoniot shock jump conditions (Sect.4.1.3). For pressure and density we can use the density ratio from the Rankine Hugoniot equation (Eq.4.13):
ρ4
ρ5
= γppc
5 +γ+ ppc
5 −1 γppc
5 +γ− ppc
5 + 1 . (4.16)
For the post-shock velocity we use Eq. 4.9. Since the pre-shock medium is at rest (u5 = 0), we can drop all terms containingu5 and find:
(pc−p5) 1
ρ5
− 1 ρ4
=u2c . (4.17)
combining Eq.4.14,4.15,4.16and4.17leads to p5ρ1
p1ρ5
1− ppc
5
2
γ 1 + ppc
5
−1 + ppc
5
= 2γ (γ−1)2
"
1− pc
p1
γ−12γ #2
. (4.18)
The solution of Eq. 4.18 can be computed with an iteration scheme. It provides one with pc. To get ρ3, this result for pc has to be inserted into Eq. 4.14; pc and Eq. 4.16 lead to ρ4. Finally the results forpcandρ4 are inserted into Eq.4.17to getuc. For the commonly used parameters in the Sod shock tube test (γ = 53, ρρ1
5 = 8, pp5
1 = 0.1) the solution ispc = 2.93945p5. It can for example be found via www.wolframalpha.com by typing
solve R*P*(x-1)^2/(x*(G+1)-1+G)=2*G/(G-1)^2(1-(P*x)^((G-1)/(2G)))^2 for G=5/3,R=8,P=0.1 .
0 0.2 0.4 0.6 0.8 1
0 0.2 0.4 x0 0.6 0.8 1
Sod shock tube at t=0.25 forγ = 53
Density Velocity Pressure
Figure 4.4: This figure shows the analytic solution of the Sod shock tube as discussed in Sect.4.2.1.
The vertical dotted lines show the zone boundaries.
The locations of the zone boundaries are found from the characteristics through the pointx0, which is the location of interface between the media at t = 0. The head of the rarefaction wave travels with the sound speed of the unperturbed high pressure medium. Thus, it is found at xh,RF = x0 −cs,1t. The velocity of the tail of this wave is found by subtracting the sound speed from the bulk velocity of the adjacent right region. The tail is thus found atxt,RF =x0+ (uc−cs,3)t. The location of thecontact discontinuityis set by the bulk velocity in this adjacent region, which leads toxCD =x0 +uctand for the shock the velocity can be found from Eq.4.6): xs =x0 + r4r4−ur5c t.
The solution for this setup at t=0.25 is shown in Fig.4.4.
4.2.2 Initial conditions of the Sod shock tube test
The typical setup of the Sod shock tube test is summarized in Fig.4.3. A density (ρ) and pressure (p) jump in the middle of the computational volume separates two gas phases withρ1 = 1,p1 = 1, ρ5 = 0.125andp5 = 0.1. In this notation subscript 1 denotes the initial state of the higher pressure gas and subscript 5 the initial state of the lower pressure gas. Subscripts 2 to 4 are reserved for the intermediate zones which will emerge later on. These zones will be separated by characteristics (characteristics are discussed in Sect. 3.4). Both gasses are initially at rest (i.e. the velocities are u1 = u5 = 0). With the adiabatic exponent of a monoatomic gas γ = 5/3 and the adiabatic EOS (Eq. 3.2) this leads to the thermal energies Etherm,1 = ei,1ρ1 = p1/(γ − 1) = 1.5 and Etherm,5 =ei,5ρ5 = 0.15.
4.2.3 Results of the R
AMSESSod shock tube test
In our simulations we will focus on thefeedback energy efficiencyof massive stars in molecular clouds. Cooling losses of the gas near the contact discontinuity (CD) play an important role for
4.2 Sod shock tube test 63
0.0 0.2 0.4 0.6 0.8 1.0
Sod shock tube at t=0.25019 forγ = 53, MonCen flux limiter, HLLC Riemann solver
Density
Passive scalar Velocity Pressure
-0.1 0.0 0.1
0.0 0.2 0.4 0.6 0.8 1
Figure 4.5: This figure shows the results of the HLLC Riemann solver with MonCen flux limiting, which is the preferred choice for our simulations. Lines in the upper panel show the analytic solutions as presented in Sect.4.2.1and Fig.4.4, the superplotted points are results of a simulation carried out with the RAMSES code. The lower panel shows the differences between the result of the simulation and the analytic solution. The residuals for all choices of Riemann solvers and flux limiters in the RAMSEScode are compared in Fig.4.6.
this study. Additionally we want to trace the products of nucleosynthesis (i.e. 26Al and 60Fe).
Thus, selecting a Riemann solver and a flux limiter that perform well near contact discontinuities is crucial. Unsurprisingly6it turned out that the acoustic and HLLC Riemann solver achieved the best results for the contact discontinuity. The results of the HLLC Riemann solver with MonCen limiting are shown in Fig.4.5.
This test was carried out with all Riemann solvers implemented in the RAMSEScode (for the con-cepts behind these solvers see Sect.3.5.1) and the MinMod and MonCen slope limiters (for details on TVD slope limiting see Sect.3.7). The results are shown in Fig.4.6. For these simulations the hydrodynamic module of the RAMSES code was used andAMRwas switched on. The minimal resolution was set to 23 cells in the computational domain. The maximal resolution was210cells per unit length. The grid was refined whenever the relative variation of density, velocity or pressure
6As discussed in Sect.3.5.1the HLLC is a variant of the HLL solver, designed to perform well atcontact discon-tinuities.
across a cell boundary was larger than 5%. In this case the data for newly refined cells was found using a MonCen interpolation scheme for the conservative variables. Moreover reflexive boundary conditions, aCFLof 0.8 (see Sect.3.3) and the MUSCL scheme were used. The intended time of the data dump wast= 0.25. To compare simulations with differentAMRgrid levels, the parame-ternsubcyclehas to be set accordingly to enforce the time step of the highest resolution grids on coarser grid to get an output at roughly the same times.
Thensubcycleparameter controls how many sub cycles will be used for the next finer level. The default value is 2 (in agreement to the dependence of the CFL condition on the grid size: if the wave can travel half the length it may travel on the coarser grid, the time step size has to be halved too). However, it is possible to set this parameter to 1. In this case the time step size from theCFL of the coarsest grid withnsubcycle=2will be used. For example for “nsubcycle=1,1,2,2” the coarsest level 1, as well as the finer levels 2 and 3 would all use theCFLof level 3, whereas the finest grid at level 4 would use its ownCFL. Settingnsubcycle=1slows down the code, which is not a problem for small scale tests like the Sod shock tube problem, but permits outputs at desired times.
As a consequence different choices of nsubcycle for the same grid levels in different AMR setups will change the output times. Data is only dumped at time steps of the coarsest grid.
E.g. if the coarse grid has 23 cells in each direction and refinement up to 210 is possible and nsubcycle=3*1,5*2 are used, the output times will be more sparse and probably differ more strongly from the desired output times than if the coarsest allowed grid has 25 cells, the finest possible grid has again 210cells and nsubcycle=3*1,3*2is used. In this simulation the whole computational box is always refined beyond24 cells in each direction.
With the default setting fornsubcyclethe actual times of the data dumps vary between the sim-ulations with different slope limiters and Riemann solvers as shown in Fig.4.7. In comparison in Fig.4.6the time-step of the210grid was also used for all coarser grids.
To avoid problems arising from the difference between the actual and the desired output times, the analytic solution (Sect. 4.2.1) was calculated for the specific end times of the individual simula-tions. Obviously the time dependence affects the locations of the zone boundaries. The small time differences between the data dumps also have a slight effect on the slopes in the rarefaction wave.
A zoom in on the residuals near thecontact discontinuity is shown in Fig.4.6. The MonCen flux limiter produces a less smeared outcontact discontinuitythan the MinMod flux limiter but at the price of over-oscillations. In this test, this can be seen in the residuals for the LLF solver, displayed in the right upper panel in Fig.4.6. Under “messier” conditions, like near the aforementionedCDin stellar winds andsupernovabubbles, also the HLLC solver sometimes happens to run into negative densities and crash the simulations. Hence we used the HLLC solver with MonCen limiting unless we ran into problems. In this case, we restarted with HLLC and MinMod.
Figure 4.5 shows the result of HLLC and MonCen, which is the most accurate setup in the set displayed in Fig.4.6. The purple line depicts the solution for a conservative passive scalar. For our purposes, it is interesting to check how diffusive thecontact discontinuityis in different numerical schemes, since this diffusivity affects the spatial distribution of our trace elements. Furthermore mixing across the CD enhances cooling losses, since dense, but cold gas and hot dilute gas will mix and lead to a dense, warm, efficiently cooling gas phase. In this setup mixing is found in about 15 cells near thecontact discontinuityatt= 0.25, as shown in Fig.4.8.
Since RAMSES can also treat 2D and 3D cases, we have also tested the dependence of the shock on the orientation of the grid. Therefore in 2D the shock tube test was once set up with the discontinuity parallel to a grid axis and once with the discontinuity along the grid diagonal. In