**3.4 Time step control**

In this section, we describe two procedures to determine the computational time steps that are used in our test calculations. The first one was originally proposed by (Rider & Knoll, 1999). The idea is to estimate the dominant wave propagation speed in the problem. In one dimension this involves calculating the ratio of temporal to spatial derivatives of the dependent variables. In principle, it is sufficient to consider the following hyperbolic equation rather than using the entire system of the governing equations

$$
\rho \frac{\partial E}{\partial t} + \upsilon\_f \frac{\partial E}{\partial r} = 0,\tag{41}
$$

where the unknown *υ<sup>f</sup>* represents the front velocity. This gives

$$\upsilon\_f = -\frac{\partial E/\partial t}{\partial E/\partial r}.\tag{42}$$

As noted in Rider & Knoll (1999), to avoid problems from lack of smoothness the following numerical approximation is used to calculate *υ<sup>f</sup>*

$$\upsilon\_f^n = \frac{\sum (|E\_i^n - E\_i^{n-1}|/\Delta t)}{\sum (|E\_{i+1}^n - E\_{i-1}^n|/2\Delta r)}.\tag{43}$$

Then the new time step is determined by the Courant-Friedrichs-Lewy (CFL) condition

$$
\Delta t^{n+1} = \mathbb{C} \frac{||\Delta r||}{\upsilon\_f^n},
\tag{44}
$$

<sup>0</sup> 0.2 0.4 0.6 0.8 <sup>1</sup> <sup>0</sup>

<sup>0</sup> 0.2 0.4 0.6 0.8 <sup>1</sup> <sup>0</sup>

−−>r

TEMPERATURE

<sup>0</sup> 0.2 0.4 0.6 0.8 <sup>1</sup> <sup>0</sup>

−−>r

, *<sup>E</sup>*Δ*<sup>t</sup>*/2, *<sup>E</sup>*Δ*<sup>t</sup>*/4, ··· ). Then we measure the *<sup>L</sup>*<sup>2</sup> norm of

PRESSURE

<sup>305</sup> An IMEX Method for the Euler Equations That Posses Strong

Fig. 3. Solution profiles resulting from the smooth problem test. The solutions are calculated

conditions for the hydrodynamics variables are *reflective* and *outflow* boundary conditions at the left and right ends of the computational domain respectively. The zero-flux boundary conditions are used for the temperature at both ends (e.g, *∂T*/*∂r*|*r*=<sup>0</sup> = 0). The coefficient of

We run the code until *t* = 0.01 with *ε*<sup>0</sup> = 100 using 200 cell points. The size of the computational domain is set to 1 (e.g, *R*<sup>0</sup> = 1 in Fig. 2). Fig. 3 shows the computed solutions for density, pressure, velocity, and temperature. As can be seen, there is no shock formation or steep thermal fronts occurred around this time. Fig. 4 shows our numerical time convergence analysis. To measure the rate of time convergence, we run the code with a fixed mesh (e.g, *M* = 200 cell points) and different time step refinements to a final time (e.g, *t* = 0.01). This

errors between *two* consecutive time step solutions (�*E*Δ*<sup>t</sup>* <sup>−</sup> *<sup>E</sup>*Δ*t*/2�2, �*E*Δ*t*/2 <sup>−</sup> *<sup>E</sup>*Δ*t*/4�2, ··· ) and plot these errors against to a second order line. It is clear from Fig. 4 that we achieve second order time convergence unlike (Bates et al., 2001) fails to provide second order accurate results

−−>r

VELOCITY

<sup>0</sup> 0.2 0.4 0.6 0.8 <sup>1</sup> <sup>0</sup>

−−>r

for *tfinal* = 0.01 with *M* = 200 cell points.

thermal conduction is set to *κ*(*T*) = *T*5/2.

provides a sequence of solution data (*E*Δ*<sup>t</sup>*

for the same test.

DENSITY

Non-Linear Heat Conduction and Stiff Source Terms (Radiation Hydrodynamics)

5

0.5

1

1.5

2

2.5

10

15

20

where � Δ*r* � uses the *L*<sup>1</sup> norm as in Equation (43). We can further simplify Equation (44) by using Equation (43), i.e,

$$
\Delta t^{n+1} = \frac{1}{2} \frac{\sum |E\_{i+1}^n - E\_{i-1}^n|}{\sum (|E\_i^n - E\_i^{n-1}| / \Delta t)}. \tag{45}
$$

We remark that the time steps determined by this procedure is always compared with the pure hydrodynamics time steps and the most restrictive ones are selected. The hydrodynamics time steps are calculated by

$$
\Delta t^{Hydro, n+1} = \mathcal{C}FL \times \frac{\Delta r}{\max\_i |u + c|\_i} \tag{46}
$$

where *u* is the fluid velocity and *c* is the sound speed (e.g, refer to Equation (28)). The coefficient *CFL* is set to 0.5. Alternative time step control criterion are used for radiation hydrodynamics problems (Bowers & Wilson, 1991). One commonly used approach is based on monitoring the maximum relative change in *E*. For instance,

$$
\Delta t^{n+1} = \Delta t^n \sqrt{\frac{(\Delta E/E)^{n+1}}{(\Delta E/E)\_{\text{max}}}} \,\text{}\tag{47}
$$

where

$$(\frac{\Delta E}{E})^{n+1} = \max\_{i} (\frac{|E\_i^{n+1} - E\_i^n|}{E\_i^{n+1} + E\_0}),\tag{48}$$

where the parameter *E*<sup>0</sup> is an estimate for the lower bound of the energy density. Comparing Equation (47) to (45) we observed that Equation (45) is computationally more efficient. Therefore, we use Equation (45) in our numerical test problems.

#### **4. Computational results**

#### **4.1 Smooth problem test**

We use the LERH model to produce numerical results for this test problem. In this test, we run the code until a particular final time so that the computational solutions are free of shock waves and steep thermal fronts. The problem is to follow the evolution of the nonlinear waves that results from an initial energy deposition in a narrow region. The initial total energy density is given by

$$E(r,0) = \frac{\varepsilon\_0 \exp\left(-r^2/c\_0^2\right)}{(c\_0\sqrt{\pi})^3},\tag{49}$$

where *c*<sup>0</sup> is a constant and set to 1/4 for this test. Note that *c*<sup>0</sup> → 0 gives a delta function at origin. We use the cell averaged values of *E* as in (Bates et al., 2001), i.e., we integrate (49) over the *i th* cell from *ri*−1/2 to *ri*<sup>+</sup>1/2 so that

$$E\_{i} = \frac{\varepsilon\_{0}[\varepsilon r f(r\_{i+1/2}/\varepsilon\_{0}) - \varepsilon r f(r\_{i-1/2}/\varepsilon\_{0})] - 2\pi \varepsilon\_{0}^{2}[r\_{i+1/2}E(r\_{i+1/2}) - r\_{i-1/2}E(r\_{i-1/2})]}{\Delta V\_{i}},\tag{50}$$

where the symbol *er f* denotes the error function. The initial density is set to *ρ* = 1/*r*. The initial temperature is calculated by using *E* = *cvρT* + <sup>1</sup> <sup>2</sup> *ρu* where *u* = 0 initially. The boundary 12 Will-be-set-by-IN-TECH

where � Δ*r* � uses the *L*<sup>1</sup> norm as in Equation (43). We can further simplify Equation (44) by

<sup>∑</sup> <sup>|</sup>*E<sup>n</sup>*

We remark that the time steps determined by this procedure is always compared with the pure hydrodynamics time steps and the most restrictive ones are selected. The hydrodynamics time

where *u* is the fluid velocity and *c* is the sound speed (e.g, refer to Equation (28)). The coefficient *CFL* is set to 0.5. Alternative time step control criterion are used for radiation hydrodynamics problems (Bowers & Wilson, 1991). One commonly used approach is based

> *n*

where the parameter *E*<sup>0</sup> is an estimate for the lower bound of the energy density. Comparing Equation (47) to (45) we observed that Equation (45) is computationally more efficient.

We use the LERH model to produce numerical results for this test problem. In this test, we run the code until a particular final time so that the computational solutions are free of shock waves and steep thermal fronts. The problem is to follow the evolution of the nonlinear waves that results from an initial energy deposition in a narrow region. The initial total energy

*<sup>E</sup>*(*r*, 0) = *<sup>ε</sup>*<sup>0</sup> exp (−*r*2/*c*<sup>2</sup>

where *c*<sup>0</sup> is a constant and set to 1/4 for this test. Note that *c*<sup>0</sup> → 0 gives a delta function at origin. We use the cell averaged values of *E* as in (Bates et al., 2001), i.e., we integrate (49) over

Δ*Vi*

where the symbol *er f* denotes the error function. The initial density is set to *ρ* = 1/*r*. The

(*c*<sup>0</sup>

0)

<sup>√</sup>*π*)<sup>3</sup> , (49)

<sup>2</sup> *ρu* where *u* = 0 initially. The boundary

, (50)

<sup>0</sup>[*ri*<sup>+</sup>1/2*E*(*ri*<sup>+</sup>1/2) − *ri*−1/2*E*(*ri*−1/2)]

<sup>∑</sup>(|*E<sup>n</sup>*

*Hydro*,*n*+<sup>1</sup> <sup>=</sup> *CFL* <sup>×</sup>

*<sup>n</sup>*+<sup>1</sup> = Δ*t*

*<sup>E</sup>* )*n*+<sup>1</sup> <sup>=</sup> *maxi*(

*<sup>i</sup>*+<sup>1</sup> <sup>−</sup> *<sup>E</sup><sup>n</sup>*

*<sup>i</sup>* <sup>−</sup> *<sup>E</sup>n*−<sup>1</sup>

*<sup>i</sup>*−1|

. (45)

, (46)

, (47)

), (48)

*<sup>i</sup>* |/Δ*t*)

Δ*r maxi*|*u* + *c*|*<sup>i</sup>*

(Δ*E*/*E*)*n*+<sup>1</sup> (Δ*E*/*E*)*max*

<sup>|</sup>*En*+<sup>1</sup> *<sup>i</sup>* <sup>−</sup> *<sup>E</sup><sup>n</sup> i* |

*En*+<sup>1</sup> *<sup>i</sup>* + *E*<sup>0</sup>

Δ*t*

Δ*t*

on monitoring the maximum relative change in *E*. For instance,

( Δ*E*

Therefore, we use Equation (45) in our numerical test problems.

Δ*t*

*<sup>n</sup>*+<sup>1</sup> <sup>=</sup> <sup>1</sup> 2

using Equation (43), i.e,

steps are calculated by

**4. Computational results**

*th* cell from *ri*−1/2 to *ri*<sup>+</sup>1/2 so that

*Ei* <sup>=</sup> *<sup>ε</sup>*0[*er f*(*ri*+1/2/*c*0) <sup>−</sup> *er f*(*ri*−1/2/*c*0)] <sup>−</sup> <sup>2</sup>*πc*<sup>2</sup>

initial temperature is calculated by using *E* = *cvρT* + <sup>1</sup>

**4.1 Smooth problem test**

density is given by

the *i*

where

Fig. 3. Solution profiles resulting from the smooth problem test. The solutions are calculated for *tfinal* = 0.01 with *M* = 200 cell points.

conditions for the hydrodynamics variables are *reflective* and *outflow* boundary conditions at the left and right ends of the computational domain respectively. The zero-flux boundary conditions are used for the temperature at both ends (e.g, *∂T*/*∂r*|*r*=<sup>0</sup> = 0). The coefficient of thermal conduction is set to *κ*(*T*) = *T*5/2.

We run the code until *t* = 0.01 with *ε*<sup>0</sup> = 100 using 200 cell points. The size of the computational domain is set to 1 (e.g, *R*<sup>0</sup> = 1 in Fig. 2). Fig. 3 shows the computed solutions for density, pressure, velocity, and temperature. As can be seen, there is no shock formation or steep thermal fronts occurred around this time. Fig. 4 shows our numerical time convergence analysis. To measure the rate of time convergence, we run the code with a fixed mesh (e.g, *M* = 200 cell points) and different time step refinements to a final time (e.g, *t* = 0.01). This provides a sequence of solution data (*E*Δ*<sup>t</sup>* , *<sup>E</sup>*Δ*<sup>t</sup>*/2, *<sup>E</sup>*Δ*<sup>t</sup>*/4, ··· ). Then we measure the *<sup>L</sup>*<sup>2</sup> norm of errors between *two* consecutive time step solutions (�*E*Δ*<sup>t</sup>* <sup>−</sup> *<sup>E</sup>*Δ*t*/2�2, �*E*Δ*t*/2 <sup>−</sup> *<sup>E</sup>*Δ*t*/4�2, ··· ) and plot these errors against to a second order line. It is clear from Fig. 4 that we achieve second order time convergence unlike (Bates et al., 2001) fails to provide second order accurate results for the same test.

is compared to Equation (6) of (Kadioglu & Knoll, 2010). Then *κ* becomes

Non-Linear Heat Conduction and Stiff Source Terms (Radiation Hydrodynamics)

*σ<sup>a</sup>* = 108 that appear in Equations (18) and (19).

*<sup>s</sup>*(*ρ*2*u*<sup>2</sup> <sup>−</sup> *<sup>ρ</sup>*1*u*1)=(*ρ*2*u*<sup>2</sup>

discontinuities.

**4.3 Radiative shock test**

*κ*(*ρ*, *T*) = *κ*<sup>0</sup>

where *<sup>κ</sup>*<sup>0</sup> <sup>=</sup> 102, *<sup>a</sup>* <sup>=</sup> <sup>−</sup>2 and *<sup>b</sup>* <sup>=</sup> 13/2 as in (Kadioglu & Knoll, 2010). We set <sup>P</sup> <sup>=</sup> <sup>10</sup>−<sup>4</sup> and

<sup>307</sup> An IMEX Method for the Euler Equations That Posses Strong

We compute the solutions until *t* = 0.02 using 400 cell points. Fig. 5 shows fluid density, fluid pressure, flow velocity, fluid energy, fluid temperature, and radiation temperature profiles. At this time (*t* = 0.02), hydrodynamical shocks are depicted near *r* = 0.2. In this test case, the thermal front (located near *r* = 0.8) propagates faster than the hydrodynamical shocks due to large initial energy deposition. Fig. 6 shows the time convergence analysis for different field variables. Clearly, we have obtained second order time accuracy for all variables. This convergence result is important, because this problem is a difficult one meaning that the coupling of different physics is highly non-linear and it is a challenge to produce fully second order convergence from an operator split method for these kinds of problems. One comment that can be made about our spatial discretization (LLF method), though it is not the primary focus of this study, is that our numerical results (figures in Fig. 5) indicate that the LLF fluxing procedure provides very good shock capturing with no spurious oscillations at or near the

The problem settings for this test are similar to (Drake, 2007; Lowrie & Edwards, 2008) where more precise physical definitions can be found. Radiative shocks are basically strong shock waves that the radiative energy flux plays essential role in the governing dynamics. Radiative shocks occur in many astrophysical systems where they move into an upstream medium leaving behind an altered downstream medium. In this test, we assume that a simple planar radiative shock exists normal to the flow as it is illustrated in Fig. 7. The initial shock profiles are determined by considering the given values in Region-1 and finding the values in Region-2 of Fig. 7. To find the values in Region-2, we use the so-called Rankine-Hugoniot relations or jump conditions (LeVeque, 1998; Smoller, 1994; Thomas, 1999). A general formula for the radiation hydrodynamics jump conditions is given in (Lowrie & Edwards, 2008). For instance

<sup>2</sup> <sup>+</sup> *<sup>p</sup>*<sup>2</sup> <sup>+</sup> <sup>P</sup>*pν*,2) <sup>−</sup> (*ρ*1*u*<sup>2</sup>

where *s* is the propagation speed of the shock front. In our test problem, we assume that the radiation temperature is smooth. Therefore, it is sufficient to use the jump conditions for the compressible Euler equations to initiate hydrodynamics shock profiles. The Euler jump conditions can be easily obtained by dropping the radiative terms in Equations (54), (55), (56),

> <sup>1</sup> <sup>+</sup> *<sup>γ</sup>* <sup>+</sup> <sup>1</sup> <sup>2</sup>*<sup>γ</sup>* ( *p*2 *p*1

and (57). Then the necessary formulae to initialize the shock solutions are

*s* = *u*<sup>1</sup> + *c*<sup>1</sup>

*s*(*E*<sup>2</sup> − *E*1) = *u*2(*E*<sup>2</sup> + *p*<sup>2</sup> + P*pν*,2) − *u*1(*E*<sup>1</sup> + *p*<sup>1</sup> + P*pν*,1), (56)

*s*(*ρ*<sup>2</sup> − *ρ*1) = *ρ*2*u*<sup>2</sup> − *ρ*1*u*1, (54)

*s*(*Eν*,2 − *Eν*,1) = *u*2(*Eν*,2) − *u*1(*Eν*,1), (57)

<sup>1</sup> + *p*<sup>1</sup> + P*pν*,1), (55)

− 1), (58)

*ρaTb*

<sup>4</sup>P*T*<sup>3</sup> , (53)

Fig. 4. Temporal convergence plot for the smooth problem test. *tfinal* = 0.01 with *M* = 200 cell points.

#### **4.2 Point explosion test**

We use the HERH model for this test. We note that we have studied this test by using both of the LERH and HERH models and reported our results in two consecutive papers (Kadioglu & Knoll, 2010; Kadioglu, Knoll, Lowrie & Rauenzahn, 2010). This section reviews our numerical findings from (Kadioglu, Knoll, Lowrie & Rauenzahn, 2010). In this test, important physics such as the propagation of sharp shock discontinuities and steep thermal fronts occur. This is important, because this test enables us to study/determine the time accuracy of the strong numerical coupling of two distinct physical processes.

Typically a point explosion is characterized by the release of large amount of energy in a small region of space (few cells near the origin). Depending on the magnitude of the energy deposition, weak or strong explosions take place. If the initial explosion energy is not large enough, the diffusive effect is limited to region behind the shock. However, if the explosion energy is large, then the thermal front can precede the hydrodynamics front. Both weak and strong explosions are studied in (Kadioglu & Knoll, 2010) where the LERH model is considered. Here, we solve/recast the strong explosion test by using the HERH model. The problem setting is as follows. The initial total energy density is given by

$$E\_0 = \frac{\varepsilon\_0 \exp\left(-r^2/c\_0^2\right)}{(c\_0\sqrt{\pi})^3},\tag{51}$$

where *ε*<sup>0</sup> = 235 and *c*<sup>0</sup> = 1/300. The initial fluid and radiation energies are set to *E*(*r*, 0) = *Eν*(*r*, 0) = *E*0/2. The fluid density is initialized by *ρ*(*r*, 0) = *r*−19/9. The initial temperature is calculated by using *E* = *cvρT*/*γ* + <sup>1</sup> <sup>2</sup> *<sup>ρ</sup>u*<sup>2</sup> with the initial *<sup>u</sup>* <sup>=</sup> 0. The radiation diffusivity (*<sup>κ</sup>* in Equation (19)) is calculated by considering the LERH model and comparing it with the sum of Equation (18) plus P times Equation (19). For instance

$$\frac{\partial}{\partial t}(E + \mathcal{P}E\_{\nu}) + \frac{1}{r^2}\frac{\partial}{\partial r}[r^2u(E + p + \mathcal{P}(E\_{\nu} + p\_{\nu}))] = \frac{1}{r^2}\frac{\partial}{\partial r}(r^2\mathcal{P}\kappa\frac{\partial E\_{\nu}}{\partial r}),\tag{52}$$

is compared to Equation (6) of (Kadioglu & Knoll, 2010). Then *κ* becomes

$$\kappa(\rho, T) = \kappa\_0 \frac{\rho^a T^b}{4 \mathcal{P} T^3} \,\,\,\,\,\,\tag{53}$$

where *<sup>κ</sup>*<sup>0</sup> <sup>=</sup> 102, *<sup>a</sup>* <sup>=</sup> <sup>−</sup>2 and *<sup>b</sup>* <sup>=</sup> 13/2 as in (Kadioglu & Knoll, 2010). We set <sup>P</sup> <sup>=</sup> <sup>10</sup>−<sup>4</sup> and *σ<sup>a</sup>* = 108 that appear in Equations (18) and (19).

We compute the solutions until *t* = 0.02 using 400 cell points. Fig. 5 shows fluid density, fluid pressure, flow velocity, fluid energy, fluid temperature, and radiation temperature profiles. At this time (*t* = 0.02), hydrodynamical shocks are depicted near *r* = 0.2. In this test case, the thermal front (located near *r* = 0.8) propagates faster than the hydrodynamical shocks due to large initial energy deposition. Fig. 6 shows the time convergence analysis for different field variables. Clearly, we have obtained second order time accuracy for all variables. This convergence result is important, because this problem is a difficult one meaning that the coupling of different physics is highly non-linear and it is a challenge to produce fully second order convergence from an operator split method for these kinds of problems. One comment that can be made about our spatial discretization (LLF method), though it is not the primary focus of this study, is that our numerical results (figures in Fig. 5) indicate that the LLF fluxing procedure provides very good shock capturing with no spurious oscillations at or near the discontinuities.

#### **4.3 Radiative shock test**

14 Will-be-set-by-IN-TECH

Temporal convergence plot for Temperature

−7 −6.5 −6 −5.5 −5 −4.5

log10Δ<sup>t</sup>

Fig. 4. Temporal convergence plot for the smooth problem test. *tfinal* = 0.01 with *M* = 200

We use the HERH model for this test. We note that we have studied this test by using both of the LERH and HERH models and reported our results in two consecutive papers (Kadioglu & Knoll, 2010; Kadioglu, Knoll, Lowrie & Rauenzahn, 2010). This section reviews our numerical findings from (Kadioglu, Knoll, Lowrie & Rauenzahn, 2010). In this test, important physics such as the propagation of sharp shock discontinuities and steep thermal fronts occur. This is important, because this test enables us to study/determine the time

Typically a point explosion is characterized by the release of large amount of energy in a small region of space (few cells near the origin). Depending on the magnitude of the energy deposition, weak or strong explosions take place. If the initial explosion energy is not large enough, the diffusive effect is limited to region behind the shock. However, if the explosion energy is large, then the thermal front can precede the hydrodynamics front. Both weak and strong explosions are studied in (Kadioglu & Knoll, 2010) where the LERH model is considered. Here, we solve/recast the strong explosion test by using the HERH model. The

*<sup>E</sup>*<sup>0</sup> <sup>=</sup> *<sup>ε</sup>*<sup>0</sup> exp (−*r*2/*c*<sup>2</sup>

(*c*<sup>0</sup>

where *ε*<sup>0</sup> = 235 and *c*<sup>0</sup> = 1/300. The initial fluid and radiation energies are set to *E*(*r*, 0) = *Eν*(*r*, 0) = *E*0/2. The fluid density is initialized by *ρ*(*r*, 0) = *r*−19/9. The initial temperature is

Equation (19)) is calculated by considering the LERH model and comparing it with the sum

[*r*2*u*(*<sup>E</sup>* <sup>+</sup> *<sup>p</sup>* <sup>+</sup> <sup>P</sup>(*E<sup>ν</sup>* <sup>+</sup> *<sup>p</sup>ν*))] = <sup>1</sup>

0)

<sup>√</sup>*π*)<sup>3</sup> , (51)

*∂E<sup>ν</sup>*

*<sup>∂</sup><sup>r</sup>* ), (52)

<sup>2</sup> *<sup>ρ</sup>u*<sup>2</sup> with the initial *<sup>u</sup>* <sup>=</sup> 0. The radiation diffusivity (*<sup>κ</sup>* in

*r*2 *∂ ∂r* (*r*2P*<sup>κ</sup>*

accuracy of the strong numerical coupling of two distinct physical processes.

problem setting is as follows. The initial total energy density is given by

−8

cell points.

**4.2 Point explosion test**

calculated by using *E* = *cvρT*/*γ* + <sup>1</sup>

*∂ ∂t*

of Equation (18) plus P times Equation (19). For instance

*r*2 *∂ ∂r*

(*<sup>E</sup>* <sup>+</sup> <sup>P</sup>*Eν*) + <sup>1</sup>

−7

−6

−5

log10L2( Error)

−4

−3

Num. Sol Sec. Order

> The problem settings for this test are similar to (Drake, 2007; Lowrie & Edwards, 2008) where more precise physical definitions can be found. Radiative shocks are basically strong shock waves that the radiative energy flux plays essential role in the governing dynamics. Radiative shocks occur in many astrophysical systems where they move into an upstream medium leaving behind an altered downstream medium. In this test, we assume that a simple planar radiative shock exists normal to the flow as it is illustrated in Fig. 7. The initial shock profiles are determined by considering the given values in Region-1 and finding the values in Region-2 of Fig. 7. To find the values in Region-2, we use the so-called Rankine-Hugoniot relations or jump conditions (LeVeque, 1998; Smoller, 1994; Thomas, 1999). A general formula for the radiation hydrodynamics jump conditions is given in (Lowrie & Edwards, 2008). For instance

$$s(\rho\_2 - \rho\_1) = \rho\_2 u\_2 - \rho\_1 u\_1 \tag{54}$$

$$s(\rho\_2 u\_2 - \rho\_1 u\_1) = (\rho\_2 u\_2^2 + p\_2 + \mathcal{P} p\_{\mathbb{V},2}) - (\rho\_1 u\_1^2 + p\_1 + \mathcal{P} p\_{\mathbb{V},1}),\tag{55}$$

$$s(E\_2 - E\_1) = \mu\_2(E\_2 + p\_2 + \mathcal{P}p\_{\nu,2}) - \mu\_1(E\_1 + p\_1 + \mathcal{P}p\_{\nu,1}),\tag{56}$$

$$s(E\_{\nu,2} - E\_{\nu,1}) = \mu\_2(E\_{\nu,2}) - \mu\_1(E\_{\nu,1}),\tag{57}$$

where *s* is the propagation speed of the shock front. In our test problem, we assume that the radiation temperature is smooth. Therefore, it is sufficient to use the jump conditions for the compressible Euler equations to initiate hydrodynamics shock profiles. The Euler jump conditions can be easily obtained by dropping the radiative terms in Equations (54), (55), (56), and (57). Then the necessary formulae to initialize the shock solutions are

$$s = u\_1 + c\_1 \sqrt{1 + \frac{\gamma + 1}{2\gamma} (\frac{p\_2}{p\_1} - 1)},\tag{58}$$

−9 −8 −7 −6 −5.5

−9 −8 −7 −6 −8.5

log10Δ t

u p

log10Δ t

FLUID TEMPERATURE Num. Sol Sec. Order

−9 −8 −7 −6 −9

−9 −8 −7 −6 −8.5

log10Δ t

log10Δ t

RADIATION TEMPERATURE Num. Sol Sec. Order

FLUID VELOCITY

Num. Sol Sec. Order

−8.5 −8 −7.5 −7 −6.5 −6 −5.5 −5 −4.5

−8 −7.5 −7 −6.5 −6 −5.5 −5 −4.5 −4 −3.5

log10L2( Error)

Fig. 6. Temporal convergence plot for various field variables from the point explosion test.

Region 1 Region 2

ρ ρ

Fig. 7. A schematic diagram of a shock wave situation with the indicated density, velocity,

1 2 1 2

1 2

u p

Shock

log10L2( Error)

<sup>309</sup> An IMEX Method for the Euler Equations That Posses Strong

FLUID DENSITY

Non-Linear Heat Conduction and Stiff Source Terms (Radiation Hydrodynamics)

Num. Sol Sec. Order

−5 −4.5 −4 −3.5 −3 −2.5 −2 −1.5 −1

−8 −7.5 −7 −6.5 −6 −5.5 −5 −4.5 −4 −3.5

*t* = 0.001 and 400 cell points are used.

and pressure for each region.

log10L2( Error)

log10L2( Error)

Fig. 5. Point explosion test with *t* = 0.02 and *M* = 400 cell points.

16 Will-be-set-by-IN-TECH

0

0

1

2

3

4

5

200

400

600

800

1000

0 0.2 0.4 0.6 0.8 1

−−>r

FLUID ENERGY

0 0.2 0.4 0.6 0.8 1

−−>r

RADIATION TEMPERATURE

0 0.2 0.4 0.6 0.8 1

−−>r

FlUID PRESSURE

0 0.2 0.4 0.6 0.8 1

−−>r

FLUID VELOCITY

0 0.2 0.4 0.6 0.8 1

−−>r

FLUID TEMPERATURE

0 0.2 0.4 0.6 0.8 1

−−>r

Fig. 5. Point explosion test with *t* = 0.02 and *M* = 400 cell points.

FLUID DENSITY

0

0

1

2

3

4

5

50

100

150

200

250

300

Fig. 6. Temporal convergence plot for various field variables from the point explosion test. *t* = 0.001 and 400 cell points are used.

Fig. 7. A schematic diagram of a shock wave situation with the indicated density, velocity, and pressure for each region.

−7.5 −7 −6.5 −6 −5.5 −5 −9.5

log10Δ t

radiative shock test. *t* = 0.02 and 400 cell points are used.

The initial radiation energy is calculated by *E<sup>ν</sup>* = *T*<sup>4</sup>

time convergence can be clearly seen in both fields.

**5. Convergence analysis**

analysis.

*<sup>T</sup>ν*(*x*, 0) = (*T*<sup>2</sup> <sup>−</sup> *<sup>T</sup>*1)

by using the energy relation *E* = *cvρT* + <sup>1</sup>

computational domain, i.e., we choose

−7.5 −7 −6.5 −6 −5.5 −5 −10.5

log10Δ t

<sup>2</sup> *<sup>ρ</sup>u*2. The radiation temperature is assumed to be a

<sup>2</sup> . (62)

*<sup>ν</sup>* . Other parameters that appear in

RADIATION TEMPERATURE

Num. Sol Sec. Order

−10 −9.5 −9 −8.5 −8 −7.5 −7 −6.5 −6 −5.5

<sup>2</sup> *tanh*(1000*x*) + (*T*<sup>2</sup> <sup>+</sup> *<sup>T</sup>*1)

log10L2( Error)

<sup>311</sup> An IMEX Method for the Euler Equations That Posses Strong

Fig. 9. Temporal convergence plot for the material and radiation temperature from the

smooth function across the shock and equal to *T*<sup>1</sup> and *T*<sup>2</sup> on the left and right boundary of the

Equations (16)-(19) are set as <sup>P</sup> <sup>=</sup> <sup>10</sup>−4, *<sup>σ</sup><sup>a</sup>* <sup>=</sup> 106, and *<sup>κ</sup>* <sup>=</sup> 1. These parameters are chosen to be consistent with (Lowrie & Edwards, 2008). We solve Equations (16)-(19), the HERH model, in Cartesian coordinates with the above initial conditions. The solutions use fixed boundary conditions at both ends. In other words, at each time step, the solutions are reset to the initial boundary values. The numerical calculations are carried out with 400 cell points and Δ*t* = 10−6. Fig. 8 shows the time history of the solutions. Notice that the solutions are highly transient, therefore it is a good test to carry out a time convergence study. Fig. 9 shows time convergence analysis for the fluid and the radiation temperature. Second order

In this section, we present a mathematical analysis (modified equation analysis) to study the analytical convergence behavior of our self-consistent IMEX method and compare it to a classic IMEX method. The modified equation analysis (truncation error analysis) is performed by considering the LERH model (Equations (8)-(10) or (20)-(21)). Also, for simplicity, we assume that the system given by Equations (20)-(21) is written in cartesian coordinates. In the introduction, we first described a classic IMEX approach then presented our self-consistent IMEX method. Therefore, we shall follow the same order in regards to below mathematical

FLUID TEMPERATURE

Non-Linear Heat Conduction and Stiff Source Terms (Radiation Hydrodynamics)

Num. Sol Sec. Order

−9 −8.5 −8 −7.5 −7 −6.5 −6 −5.5 −5

log10L2( Error)

Fig. 8. Time evolution of the left moving radiative shock.

$$\frac{p\_2}{p\_1} = 1 + \frac{2\gamma}{\gamma + 1} [(\frac{s - \mu\_1}{c\_1})^2 - 1],\tag{59}$$

$$
\mu\_2 = \mu\_1 + \frac{p\_2 - p\_1}{\rho\_1 (s - \mu\_1)} \,\prime \tag{60}
$$

$$\frac{\rho\_2}{\rho\_1} = \frac{s - \mu\_1}{s - \mu\_2} \,\,\,\,\tag{61}$$

where *c*<sup>1</sup> = *γ <sup>p</sup>*<sup>1</sup> *<sup>ρ</sup>*<sup>1</sup> is the speed of sound in the fluid at upstream conditions. More details regarding the derivation of Equations (58)-(61) can be found in (Anderson, 1990; LeVeque, 1998; Smoller, 1994; Thomas, 1999; Wesseling, 2000).

We are interested in solving a left moving radiative shock problem. To achieve this, we set the initial shock speed *s* = −0.1 in Equation (58). Other upstream flow variables are set as follows; *ρ*<sup>1</sup> = 1.0, *T*<sup>1</sup> = 1.0, and *M*<sup>1</sup> = *u*1/*c*<sup>1</sup> = 1.2 as the upstream Mach number. Then we calculate the pressure from a calorically perfect gas relation (*p*<sup>1</sup> = *Rρ*1*T*1). Using *p*<sup>1</sup> and *<sup>ρ</sup>*1, we calculate the upstream sound speed *<sup>c</sup>*<sup>1</sup> = *<sup>γ</sup>p*1/*ρ*<sup>1</sup> together with *<sup>u</sup>*<sup>1</sup> = *<sup>M</sup>*1*c*1. With these information in hand, we can easily calculate the downstream values using Equations (58)-(61). The total fluid energies in both upstream and downstream directions are calculated

Fig. 9. Temporal convergence plot for the material and radiation temperature from the radiative shock test. *t* = 0.02 and 400 cell points are used.

by using the energy relation *E* = *cvρT* + <sup>1</sup> <sup>2</sup> *<sup>ρ</sup>u*2. The radiation temperature is assumed to be a smooth function across the shock and equal to *T*<sup>1</sup> and *T*<sup>2</sup> on the left and right boundary of the computational domain, i.e., we choose

$$T\_{\nu}(\mathbf{x},0) = \frac{(T\_2 - T\_1)}{2} \tanh(1000\mathbf{x}) + \frac{(T\_2 + T\_1)}{2}.\tag{62}$$

The initial radiation energy is calculated by *E<sup>ν</sup>* = *T*<sup>4</sup> *<sup>ν</sup>* . Other parameters that appear in Equations (16)-(19) are set as <sup>P</sup> <sup>=</sup> <sup>10</sup>−4, *<sup>σ</sup><sup>a</sup>* <sup>=</sup> 106, and *<sup>κ</sup>* <sup>=</sup> 1. These parameters are chosen to be consistent with (Lowrie & Edwards, 2008). We solve Equations (16)-(19), the HERH model, in Cartesian coordinates with the above initial conditions. The solutions use fixed boundary conditions at both ends. In other words, at each time step, the solutions are reset to the initial boundary values. The numerical calculations are carried out with 400 cell points and Δ*t* = 10−6. Fig. 8 shows the time history of the solutions. Notice that the solutions are highly transient, therefore it is a good test to carry out a time convergence study. Fig. 9 shows time convergence analysis for the fluid and the radiation temperature. Second order time convergence can be clearly seen in both fields.

#### **5. Convergence analysis**

18 Will-be-set-by-IN-TECH

1.1

1

1.05

[(*<sup>s</sup>* <sup>−</sup> *<sup>u</sup>*<sup>1</sup> *c*1

*ρ*1(*s* − *u*1)

*<sup>ρ</sup>*<sup>1</sup> is the speed of sound in the fluid at upstream conditions. More details

1.1

1.15

1.2

1.25

1.2

1.3

1.4

1.5

−0.01 −0.005 0 0.005 0.01

−0.01 −0.005 0 0.005 0.01

−−>x

)<sup>2</sup> <sup>−</sup> <sup>1</sup>], (59)

, (60)

, (61)

t=0.00 t=0.02 t=0.04 t=0.06

−−>x

RADIATION TEMPERATURE

FLUID VELOCITY

−0.01 −0.005 0 0.005 0.01

−0.01 −0.005 0 0.005 0.01

−−>x

Fig. 8. Time evolution of the left moving radiative shock.

1998; Smoller, 1994; Thomas, 1999; Wesseling, 2000).

*p*2 *p*1

= 1 +

2*γ γ* + 1

*ρ*2 *ρ*1

*<sup>u</sup>*<sup>2</sup> <sup>=</sup> *<sup>u</sup>*<sup>1</sup> <sup>+</sup> *<sup>p</sup>*<sup>2</sup> <sup>−</sup> *<sup>p</sup>*<sup>1</sup>

<sup>=</sup> *<sup>s</sup>* <sup>−</sup> *<sup>u</sup>*<sup>1</sup> *s* − *u*<sup>2</sup>

regarding the derivation of Equations (58)-(61) can be found in (Anderson, 1990; LeVeque,

We are interested in solving a left moving radiative shock problem. To achieve this, we set the initial shock speed *s* = −0.1 in Equation (58). Other upstream flow variables are set as follows; *ρ*<sup>1</sup> = 1.0, *T*<sup>1</sup> = 1.0, and *M*<sup>1</sup> = *u*1/*c*<sup>1</sup> = 1.2 as the upstream Mach number. Then we calculate the pressure from a calorically perfect gas relation (*p*<sup>1</sup> = *Rρ*1*T*1). Using *p*<sup>1</sup> and *<sup>ρ</sup>*1, we calculate the upstream sound speed *<sup>c</sup>*<sup>1</sup> = *<sup>γ</sup>p*1/*ρ*<sup>1</sup> together with *<sup>u</sup>*<sup>1</sup> = *<sup>M</sup>*1*c*1. With these information in hand, we can easily calculate the downstream values using Equations (58)-(61). The total fluid energies in both upstream and downstream directions are calculated

−−>x

FLUID TEMPERATURE

FLUID DENSITY

1 1.05 1.1 1.15 1.2 1.25 1.3 1.35 1.4

1

1.05

where *c*<sup>1</sup> =

 *γ <sup>p</sup>*<sup>1</sup>

1.1

1.15

1.2

1.25

In this section, we present a mathematical analysis (modified equation analysis) to study the analytical convergence behavior of our self-consistent IMEX method and compare it to a classic IMEX method. The modified equation analysis (truncation error analysis) is performed by considering the LERH model (Equations (8)-(10) or (20)-(21)). Also, for simplicity, we assume that the system given by Equations (20)-(21) is written in cartesian coordinates. In the introduction, we first described a classic IMEX approach then presented our self-consistent IMEX method. Therefore, we shall follow the same order in regards to below mathematical analysis.

Carrying out the necessary algebra, Equation (69) becomes

*<sup>∂</sup><sup>x</sup>* {*u*1[*cvρ*1*T*<sup>1</sup> <sup>+</sup> <sup>1</sup>

<sup>2</sup> *<sup>L</sup>*(*En*) + <sup>Δ</sup>*<sup>t</sup>*

Non-Linear Heat Conduction and Stiff Source Terms (Radiation Hydrodynamics)

<sup>2</sup> *<sup>L</sup>*(*E*1) + <sup>Δ</sup>*<sup>t</sup>*

<sup>313</sup> An IMEX Method for the Euler Equations That Posses Strong

*L*(*E*1) = *L*(*En*) + Δ*t*

2 2 *∂L <sup>∂</sup><sup>t</sup>* <sup>+</sup>

Now, we consider the following Taylor series for the implicit block (Equation (65))

*∂E<sup>n</sup> <sup>∂</sup><sup>t</sup>* <sup>+</sup>

*∂T<sup>n</sup> <sup>∂</sup><sup>t</sup>* <sup>+</sup>

*∂κ<sup>n</sup> <sup>∂</sup><sup>t</sup>* <sup>+</sup>

*∂*2*E<sup>n</sup>*

2 <sup>2</sup> *<sup>ρ</sup>*1*u*<sup>1</sup> *<sup>∂</sup>T<sup>n</sup>*

<sup>2</sup> *<sup>ρ</sup>*1(*u*1)<sup>2</sup> <sup>+</sup> *<sup>R</sup>ρ*1*T*1]} <sup>=</sup> <sup>−</sup> *<sup>∂</sup>*

*∂L*

Δ*t* 2 <sup>2</sup> *<sup>ρ</sup>*1*u*<sup>1</sup> *<sup>∂</sup>T<sup>n</sup>*

> Δ*t* 2 2

Δ*t* 2 2

Δ*t* 2 2

Substituting Equations (73), (74), (75), and (72) into Equation (65), we form the truncation term

1 2 *∂ ∂x*

*<sup>∂</sup>t*<sup>2</sup> )] <sup>−</sup> <sup>Δ</sup>*<sup>t</sup>*

2 2 *∂ ∂t* [ *∂E<sup>n</sup>*

> 2 2 *∂ ∂x* (*κ<sup>n</sup> t ∂T<sup>n</sup> <sup>∂</sup><sup>x</sup>* )]

*<sup>∂</sup><sup>t</sup>* (*cv* <sup>+</sup> *<sup>R</sup>*) + *<sup>O</sup>*(Δ*<sup>t</sup>*

*∂*2*T<sup>n</sup>*

Cancelling the opposite sign common terms and grouping the other terms together, we get

*<sup>∂</sup><sup>x</sup>* ) + <sup>Δ</sup>*<sup>t</sup>* 2 *∂ ∂x*

*<sup>∂</sup><sup>t</sup>* <sup>−</sup> *<sup>L</sup>*(*En*)] + <sup>Δ</sup>*<sup>t</sup>*

(*κ<sup>n</sup> <sup>∂</sup>T<sup>n</sup>*

(*κ<sup>n</sup> <sup>∂</sup>T<sup>n</sup> t <sup>∂</sup><sup>x</sup>* ) + <sup>Δ</sup>*<sup>t</sup>*

*<sup>∂</sup><sup>t</sup>* <sup>+</sup> *<sup>O</sup>*(Δ*<sup>t</sup>*

*∂*2*E<sup>n</sup>*

*∂*2*T<sup>n</sup>*

*∂*2*κ<sup>n</sup>*

*<sup>∂</sup>t*<sup>2</sup> <sup>−</sup> [*E<sup>n</sup>* <sup>+</sup> <sup>Δ</sup>*tL*(*En*) + <sup>Δ</sup>*<sup>t</sup>*

[(*κ<sup>n</sup>* + Δ*t*

1 2 *∂ ∂x*

*<sup>∂</sup>t*<sup>2</sup> <sup>+</sup> *<sup>O</sup>*(Δ*<sup>t</sup>*

*<sup>∂</sup>t*<sup>2</sup> <sup>+</sup> *<sup>O</sup>*(Δ*<sup>t</sup>*

*<sup>∂</sup>t*<sup>2</sup> <sup>+</sup> *<sup>O</sup>*(Δ*<sup>t</sup>*

*<sup>∂</sup><sup>t</sup>* (*cv* <sup>+</sup> *<sup>R</sup>*) + *<sup>O</sup>*(Δ*<sup>t</sup>*

*<sup>∂</sup><sup>t</sup>* (*cv* <sup>+</sup> *<sup>R</sup>*) + *<sup>O</sup>*(Δ*<sup>t</sup>*

2 2

> Δ*t* 2 2

*<sup>∂</sup><sup>x</sup>* ) + *<sup>O</sup>*(Δ*<sup>t</sup>*

*∂κ<sup>n</sup> <sup>∂</sup><sup>t</sup>* <sup>+</sup>

(*κ<sup>n</sup> <sup>∂</sup>T<sup>n</sup>*

*<sup>∂</sup><sup>t</sup>* <sup>−</sup> *<sup>L</sup>*(*En*)]

(*κ<sup>n</sup> <sup>∂</sup>T<sup>n</sup> t*

(*κ<sup>n</sup> <sup>∂</sup>T<sup>n</sup> <sup>∂</sup><sup>x</sup>* )] *∂L*(*En*) *∂t*

> *∂*2*κ<sup>n</sup> <sup>∂</sup>t*<sup>2</sup> )

<sup>3</sup>). (77)

*<sup>∂</sup><sup>x</sup>* ). (78)

<sup>3</sup>), (70)

<sup>3</sup>). (72)

<sup>3</sup>). (76)

*<sup>∂</sup><sup>x</sup>* {*u*1[*E*<sup>1</sup> <sup>+</sup> *<sup>p</sup>*1]}. Further

<sup>3</sup>), (73)

<sup>3</sup>), (74)

<sup>3</sup>). (75)

<sup>2</sup>). (71)

Δ*t*

simplification comes from the following identity

*<sup>E</sup>*<sup>∗</sup> <sup>=</sup> *<sup>E</sup><sup>n</sup>* <sup>+</sup> <sup>Δ</sup>*tL*(*En*) + <sup>Δ</sup>*<sup>t</sup>*

*En*+<sup>1</sup> = *E<sup>n</sup>* + Δ*t*

*Tn*+<sup>1</sup> = *T<sup>n</sup>* + Δ*t*

*κn*+<sup>1</sup> = *κ<sup>n</sup>* + Δ*t*

Δ*t* 2 2

*<sup>∂</sup><sup>t</sup>* (*cv* <sup>+</sup> *<sup>R</sup>*)] <sup>−</sup> <sup>Δ</sup>*<sup>t</sup>*

Δ*t* 2 2

*∂T<sup>n</sup> <sup>∂</sup><sup>t</sup>* <sup>+</sup>

*∂E<sup>n</sup>*

*∂E<sup>n</sup> <sup>∂</sup><sup>t</sup>* <sup>+</sup>

(*T<sup>n</sup>* + Δ*t*

*τ<sup>n</sup>* = Δ*t* [

− [ Δ*t* 2 *∂ ∂x*

− [ Δ*t* 2 2 *∂ ∂x*

+ Δ*t* 2 <sup>2</sup> *<sup>ρ</sup>*1*u*<sup>1</sup> *<sup>∂</sup>T<sup>n</sup>*

(*κ<sup>n</sup> <sup>∂</sup>T<sup>n</sup>*

*<sup>∂</sup><sup>x</sup>* )] = *<sup>∂</sup>*

*∂x* (*κ<sup>n</sup> t ∂T<sup>n</sup> <sup>∂</sup><sup>x</sup>* ) + *<sup>∂</sup> ∂x*

*∂ ∂t* [ *∂ ∂x*

*E*<sup>∗</sup> = *E<sup>n</sup>* +

Inserting Equation (71) into (70), we get

*τ<sup>n</sup>* = *E<sup>n</sup>* + Δ*t*

*∂ ∂x*

+ Δ*t* 2 <sup>2</sup> *<sup>ρ</sup>*1*u*<sup>1</sup> *<sup>∂</sup>T<sup>n</sup>*

This further simplifies by using

where *<sup>L</sup>*(*E*1) = <sup>−</sup> *<sup>∂</sup>*

as

#### **5.1 A classic IMEX method**

The classic IMEX method operates on Equations (20)-(21) as follows. Explicit block: Step-1:

$$
\rho^1 = \rho^n - \Delta t \frac{\partial}{\partial \mathbf{x}} (\rho u)^n,
$$

$$
(\rho u)^1 = (\rho u)^n - \Delta t \frac{\partial}{\partial \mathbf{x}} [(\rho u^2)^n + p^n],
$$

$$
E^1 = E^n - \Delta t \frac{\partial}{\partial \mathbf{x}} \{u^n [E^n + p^n] \},
\tag{63}
$$

Step-2:

$$\rho^{n+1} = \frac{\rho^1 + \rho^n}{2} - \frac{\Delta t}{2} \frac{\partial}{\partial \mathbf{x}} (\rho u)^1,$$

$$(\rho u)^{n+1} = \frac{(\rho u)^1 + (\rho u)^n}{2} - \frac{\Delta t}{2} \frac{\partial}{\partial \mathbf{x}} [\rho^1 (u^2)^1 + R \rho^1 \mathbf{T}^n]\_\prime$$

$$E^\* = \frac{E^n + E^1}{2} - \frac{\Delta t}{2} \frac{\partial}{\partial \mathbf{x}} \{u^1 [c\_{\overline{\nu}} \rho^1 \mathbf{T}^n + \frac{1}{2} \rho^1 (u^1)^2 + R \rho^1 \mathbf{T}^n] \}\_\prime \tag{64}$$

Implicit block:

$$\frac{E^{n+1} - E^\*}{\Delta t} = \frac{1}{2} \frac{\partial}{\partial \mathbf{x}} (\kappa^{n+1} \frac{\partial T^{n+1}}{\partial \mathbf{x}}) + \frac{1}{2} \frac{\partial}{\partial \mathbf{x}} (\kappa^{n!} \frac{\partial T^n}{\partial \mathbf{x}}),\tag{65}$$

where we incorporated with the equation of states relations plus we assume that the explicit block is based on a second order TVD Runga-Kutta method and the implicit block is similar to the Crank-Nicolson method. Notice that the classic IMEX method is executed in such a way that the implicit temperature does not impact the explicit block (refer to the highlighted terms in Equation (64)). We carry out the modified equation analysis for the energy part of Equations (63)-(65), but the same procedure can easily be extended to the whole system. We consider

$$E^1 = E^n - \Delta t \frac{\partial}{\partial \mathbf{x}} \{ \mu^n [E^n + p^n] \},\tag{66}$$

$$E^\* = \frac{E^\hbar + E^1}{2} - \frac{\Delta t}{2} \frac{\partial}{\partial \mathbf{x}} \{ u^1 [c\_\sigma \rho^1 T^\mu + \frac{1}{2} \rho^1 (u^1)^2 + \mathcal{R} \rho^1 T^\mu] \},\tag{67}$$

Substituting Equation (66) into (67), we get

$$E^\* = E^n - \frac{\Delta t}{2} \frac{\partial}{\partial \mathbf{x}} \{ u^n [E^n + p^n] \} - \frac{\Delta t}{2} \frac{\partial}{\partial \mathbf{x}} \{ u^1 [c\_\eta \rho^1 T^n + \frac{1}{2} \rho^1 (u^1)^2 + R \rho^1 T^n] \}. \tag{68}$$

We let *<sup>L</sup>*(*En*) = <sup>−</sup> *<sup>∂</sup> <sup>∂</sup><sup>x</sup>* {*un*[*E<sup>n</sup>* <sup>+</sup> *<sup>p</sup>n*]} and use *<sup>T</sup><sup>n</sup>* <sup>=</sup> *<sup>T</sup>*<sup>1</sup> <sup>−</sup> <sup>Δ</sup>*tT<sup>n</sup> <sup>t</sup>* + *O*(Δ*t* <sup>2</sup>) , then (68) becomes

$$E^\* = E^n + \frac{\Delta t}{2} L(E^n) - \frac{\Delta t}{2} \frac{\partial}{\partial x} \{ u^1 [c\_v \rho^1 (T^1 - \Delta t \frac{\partial T^n}{\partial t} + O(\Delta t^2)) \}$$

$$+ \frac{1}{2} \rho^1 (u^1)^2 + R \rho^1 (T^1 - \Delta t \frac{\partial T^n}{\partial t} + O(\Delta t^2)) ].\tag{69}$$

Carrying out the necessary algebra, Equation (69) becomes

$$E^\* = E^\hbar + \frac{\Delta t}{2} L(E^\hbar) + \frac{\Delta t}{2} L(E^1) + \frac{\Delta t^2}{2} \rho^1 u^1 \frac{\partial T^\hbar}{\partial t} (\mathbf{c}\_\upsilon + \mathbf{R}) + O(\Delta t^3), \tag{70}$$

where *<sup>L</sup>*(*E*1) = <sup>−</sup> *<sup>∂</sup> <sup>∂</sup><sup>x</sup>* {*u*1[*cvρ*1*T*<sup>1</sup> <sup>+</sup> <sup>1</sup> <sup>2</sup> *<sup>ρ</sup>*1(*u*1)<sup>2</sup> <sup>+</sup> *<sup>R</sup>ρ*1*T*1]} <sup>=</sup> <sup>−</sup> *<sup>∂</sup> <sup>∂</sup><sup>x</sup>* {*u*1[*E*<sup>1</sup> <sup>+</sup> *<sup>p</sup>*1]}. Further simplification comes from the following identity

$$L(E^1) = L(E^n) + \Delta t \frac{\partial L}{\partial t} + O(\Delta t^2). \tag{71}$$

Inserting Equation (71) into (70), we get

20 Will-be-set-by-IN-TECH

*∂ ∂x*

> *∂ ∂x*

(*ρu*)*n*,

[(*ρu*2)*<sup>n</sup>* + *pn*],

[*ρ*1(*u*2)<sup>1</sup> + *Rρ*1**Tn**],

1 2

(*κ<sup>n</sup> <sup>∂</sup>T<sup>n</sup>*

{*un*[*E<sup>n</sup>* <sup>+</sup> *<sup>p</sup>n*]}, (66)

1 2 *ρ*1(*u*1)

*<sup>t</sup>* + *O*(Δ*t*

*∂T<sup>n</sup>*

*∂T<sup>n</sup>*

*<sup>∂</sup><sup>t</sup>* <sup>+</sup> *<sup>O</sup>*(Δ*<sup>t</sup>*

*<sup>∂</sup><sup>t</sup>* <sup>+</sup> *<sup>O</sup>*(Δ*<sup>t</sup>*

{*un*[*E<sup>n</sup>* <sup>+</sup> *<sup>p</sup>n*]}, (63)

*<sup>ρ</sup>*1(*u*1)<sup>2</sup> <sup>+</sup> *<sup>R</sup>ρ*1**Tn**]}, (64)

*<sup>ρ</sup>*1(*u*1)<sup>2</sup> <sup>+</sup> *<sup>R</sup>ρ*1*Tn*]}, (67)

<sup>2</sup>))

<sup>2</sup> <sup>+</sup> *<sup>R</sup>ρ*1*Tn*]}. (68)

<sup>2</sup>))]}. (69)

<sup>2</sup>) , then (68) becomes

*<sup>∂</sup><sup>x</sup>* ), (65)

*∂ ∂x*

The classic IMEX method operates on Equations (20)-(21) as follows.

*<sup>ρ</sup>*<sup>1</sup> <sup>=</sup> *<sup>ρ</sup><sup>n</sup>* <sup>−</sup> <sup>Δ</sup>*<sup>t</sup>*

(*ρu*)<sup>1</sup> = (*ρu*)*<sup>n</sup>* <sup>−</sup> <sup>Δ</sup>*<sup>t</sup>*

*<sup>E</sup>*<sup>1</sup> <sup>=</sup> *<sup>E</sup><sup>n</sup>* <sup>−</sup> <sup>Δ</sup>*<sup>t</sup>*

<sup>2</sup> <sup>−</sup> <sup>Δ</sup>*<sup>t</sup>*

2 *∂ ∂x*

{*u*1[*cvρ*1**T<sup>n</sup>** <sup>+</sup>

*<sup>∂</sup><sup>x</sup>* ) + <sup>1</sup> 2 *∂ ∂x*

(*κn*+<sup>1</sup> *<sup>∂</sup>Tn*+<sup>1</sup>

where we incorporated with the equation of states relations plus we assume that the explicit block is based on a second order TVD Runga-Kutta method and the implicit block is similar to the Crank-Nicolson method. Notice that the classic IMEX method is executed in such a way that the implicit temperature does not impact the explicit block (refer to the highlighted terms in Equation (64)). We carry out the modified equation analysis for the energy part of Equations (63)-(65), but the same procedure can easily be extended to the whole system. We

{*u*1[*cvρ*1*T<sup>n</sup>* <sup>+</sup>

2 *∂ ∂x* 1 2

{*u*1[*cvρ*1*T<sup>n</sup>* <sup>+</sup>

{*u*1[*cvρ*1(*T*<sup>1</sup> <sup>−</sup> <sup>Δ</sup>*<sup>t</sup>*

<sup>2</sup> <sup>+</sup> *<sup>R</sup>ρ*1(*T*<sup>1</sup> <sup>−</sup> <sup>Δ</sup>*<sup>t</sup>*

<sup>2</sup> <sup>−</sup> <sup>Δ</sup>*<sup>t</sup>* 2 *∂ ∂x* (*ρu*) 1,

<sup>2</sup> <sup>−</sup> <sup>Δ</sup>*<sup>t</sup>* 2 *∂ ∂x*

> 2 *∂ ∂x*

*<sup>ρ</sup>n*+<sup>1</sup> <sup>=</sup> *<sup>ρ</sup>*<sup>1</sup> <sup>+</sup> *<sup>ρ</sup><sup>n</sup>*

(*ρu*)*n*+<sup>1</sup> <sup>=</sup> (*ρu*)<sup>1</sup> + (*ρu*)*<sup>n</sup>*

*<sup>E</sup>*<sup>∗</sup> <sup>=</sup> *<sup>E</sup><sup>n</sup>* <sup>+</sup> *<sup>E</sup>*<sup>1</sup>

*<sup>E</sup>n*+<sup>1</sup> <sup>−</sup> *<sup>E</sup>*<sup>∗</sup> <sup>Δ</sup>*<sup>t</sup>* <sup>=</sup> <sup>1</sup>

*<sup>E</sup>*<sup>1</sup> <sup>=</sup> *<sup>E</sup><sup>n</sup>* <sup>−</sup> <sup>Δ</sup>*<sup>t</sup>*

*<sup>E</sup>*<sup>∗</sup> <sup>=</sup> *<sup>E</sup><sup>n</sup>* <sup>+</sup> *<sup>E</sup>*<sup>1</sup>

Substituting Equation (66) into (67), we get

2 *∂ ∂x*

Δ*t*

*E*<sup>∗</sup> = *E<sup>n</sup>* +

*<sup>E</sup>*<sup>∗</sup> <sup>=</sup> *<sup>E</sup><sup>n</sup>* <sup>−</sup> <sup>Δ</sup>*<sup>t</sup>*

We let *<sup>L</sup>*(*En*) = <sup>−</sup> *<sup>∂</sup>*

*∂ ∂x*

<sup>2</sup> <sup>−</sup> <sup>Δ</sup>*<sup>t</sup>* 2 *∂ ∂x*

<sup>2</sup> *<sup>L</sup>*(*En*) <sup>−</sup> <sup>Δ</sup>*<sup>t</sup>*

+ 1 2 *ρ*1(*u*1)

{*un*[*E<sup>n</sup>* <sup>+</sup> *<sup>p</sup>n*]} − <sup>Δ</sup>*<sup>t</sup>*

*<sup>∂</sup><sup>x</sup>* {*un*[*E<sup>n</sup>* <sup>+</sup> *<sup>p</sup>n*]} and use *<sup>T</sup><sup>n</sup>* <sup>=</sup> *<sup>T</sup>*<sup>1</sup> <sup>−</sup> <sup>Δ</sup>*tT<sup>n</sup>*

2 *∂ ∂x*

**5.1 A classic IMEX method**

Explicit block: Step-1:

Step-2:

Implicit block:

consider

$$E^\* = E^\eta + \Delta t L(E^\eta) + \frac{\Delta t^2}{2} \frac{\partial L}{\partial t} + \frac{\Delta t^2}{2} \rho^1 u^1 \frac{\partial T^\eta}{\partial t} (c\_v + R) + O(\Delta t^3). \tag{72}$$

Now, we consider the following Taylor series for the implicit block (Equation (65))

$$E^{\prime +1} = E^{\prime \prime} + \Delta t \frac{\partial E^{\prime \prime}}{\partial t} + \frac{\Delta t^2}{2} \frac{\partial^2 E^{\prime \prime}}{\partial t^2} + \mathcal{O}(\Delta t^3), \tag{73}$$

$$T^{n+1} = T^n + \Delta t \frac{\partial T^n}{\partial t} + \frac{\Delta t^2}{2} \frac{\partial^2 T^n}{\partial t^2} + O(\Delta t^3),\tag{74}$$

$$
\kappa^{n+1} = \kappa^n + \Delta t \frac{\partial \kappa^n}{\partial t} + \frac{\Delta t^2}{2} \frac{\partial^2 \kappa^n}{\partial t^2} + O(\Delta t^3). \tag{75}
$$

Substituting Equations (73), (74), (75), and (72) into Equation (65), we form the truncation term as

$$\tau^{n} = E^{n} + \Delta t \frac{\partial E^{n}}{\partial t} + \frac{\Delta t^{2}}{2} \frac{\partial^{2} E^{n}}{\partial t^{2}} - [E^{n} + \Delta t L(E^{n}) + \frac{\Delta t^{2}}{2} \frac{\partial L(E^{n})}{\partial t}]$$

$$+ \frac{\Delta t^{2}}{2} \rho^{1} u^{1} \frac{\partial T^{n}}{\partial t} (c\_{\mathcal{V}} + R) [-\Delta t \frac{1}{2} \frac{\partial}{\partial x} [(\kappa^{n} + \Delta t \frac{\partial \kappa^{n}}{\partial t} + \frac{\Delta t^{2}}{2} \frac{\partial^{2} \kappa^{n}}{\partial t^{2}})$$

$$\frac{\partial}{\partial x} (T^{n} + \Delta t \frac{\partial T^{n}}{\partial t} + \frac{\Delta t^{2}}{2} \frac{\partial^{2} T^{n}}{\partial t^{2}}) [-\Delta t \frac{1}{2} \frac{\partial}{\partial x} (\kappa^{n} \frac{\partial T^{n}}{\partial x}) + O(\Delta t^{3}).\tag{76}$$

Cancelling the opposite sign common terms and grouping the other terms together, we get

$$\begin{split} \pi^{n} &= \Delta t \left[ \frac{\partial E^{n}}{\partial t} - L(E^{n}) \right] + \frac{\Delta t^{2}}{2} \frac{\partial}{\partial t} \left[ \frac{\partial E^{n}}{\partial t} - L(E^{n}) \right] \\ &- \left[ \frac{\Delta t}{2} \frac{\partial}{\partial \mathbf{x}} \left( \kappa^{n} \frac{\partial T^{n}}{\partial \mathbf{x}} \right) + \frac{\Delta t}{2} \frac{\partial}{\partial \mathbf{x}} \left( \kappa^{n} \frac{\partial T^{n}}{\partial \mathbf{x}} \right) \right] \\ &- \left[ \frac{\Delta t^{2}}{2} \frac{\partial}{\partial \mathbf{x}} \left( \kappa^{n} \frac{\partial T^{n}\_{t}}{\partial \mathbf{x}} \right) + \frac{\Delta t^{2}}{2} \frac{\partial}{\partial \mathbf{x}} \left( \kappa^{n}\_{t} \frac{\partial T^{n}}{\partial \mathbf{x}} \right) \right] \\ &+ \frac{\Delta t^{2}}{2} \rho^{1} u^{1} \frac{\partial T^{n}}{\partial t} (c\_{v} + R) + O(\Delta t^{3}). \tag{77} \end{split} \tag{78}$$

This further simplifies by using

$$\frac{\partial}{\partial t} [\frac{\partial}{\partial \mathbf{x}} (\kappa^n \frac{\partial T^n}{\partial \mathbf{x}})] = \frac{\partial}{\partial \mathbf{x}} (\kappa^n\_t \frac{\partial T^n}{\partial \mathbf{x}}) + \frac{\partial}{\partial \mathbf{x}} (\kappa^n \frac{\partial T^n\_t}{\partial \mathbf{x}}).\tag{78}$$

We let *<sup>L</sup>*(*En*) = <sup>−</sup> *<sup>∂</sup>*

where *<sup>L</sup>*(*E*1) = <sup>−</sup> *<sup>∂</sup>*

by using (71),

*E*<sup>∗</sup> = *E<sup>n</sup>* +

Now, we insert *<sup>T</sup><sup>n</sup>* <sup>=</sup> *<sup>T</sup>*<sup>1</sup> <sup>−</sup> <sup>Δ</sup>*<sup>t</sup> <sup>∂</sup>T<sup>n</sup>*

*τ<sup>n</sup>* = *E<sup>n</sup>* + Δ*t*

− Δ*t* 1 2 *∂ ∂x*

− Δ*t* 1 2 *∂ ∂x*

*<sup>∂</sup><sup>t</sup>* <sup>−</sup> *<sup>L</sup>*(*En*) <sup>−</sup> *<sup>∂</sup>*

Again from the energy equation, we know that *<sup>∂</sup>E<sup>n</sup>*

use of Equation (78), we get

*∂E<sup>n</sup>*

*τ<sup>n</sup>* = Δ*t* [

(89) becomes

2010).

Δ*t*

*<sup>∂</sup><sup>x</sup>* {*u*1[*cvρ*1*T*<sup>1</sup> <sup>+</sup> <sup>1</sup>

<sup>2</sup> *<sup>L</sup>*(*En*) <sup>−</sup> <sup>Δ</sup>*<sup>t</sup>*

+ 1 2 *ρ*1(*u*1)

Non-Linear Heat Conduction and Stiff Source Terms (Radiation Hydrodynamics)

*<sup>∂</sup><sup>t</sup>* + *O*(Δ*t*

Δ*t* 2

*<sup>E</sup>*<sup>∗</sup> <sup>=</sup> *<sup>E</sup><sup>n</sup>* <sup>+</sup> <sup>Δ</sup>*tL*(*En*) + <sup>Δ</sup>*<sup>t</sup>*

*∂*2*E<sup>n</sup>*

Δ*t* 2 2

*<sup>∂</sup><sup>x</sup>* )] + <sup>Δ</sup>*<sup>t</sup>*

*E*<sup>∗</sup> = *E<sup>n</sup>* +

implicit discretization (83), we form the truncation error term as

Δ*t* 2 2

> *∂κ<sup>n</sup> <sup>∂</sup><sup>t</sup>* <sup>+</sup>

*<sup>∂</sup><sup>x</sup>* ) + *<sup>O</sup>*(Δ*<sup>t</sup>*

(*κ<sup>n</sup> <sup>∂</sup>T<sup>n</sup>*

clearly proving that the self-consistent IMEX method is second order.

*∂E<sup>n</sup> <sup>∂</sup><sup>t</sup>* <sup>+</sup>

[(*κ<sup>n</sup>* + Δ*t*

(*κ<sup>n</sup> <sup>∂</sup>T<sup>n</sup>*

*∂x*

2 *∂ ∂x*

*<sup>∂</sup><sup>x</sup>* {*un*[*E<sup>n</sup>* <sup>+</sup> *<sup>p</sup>n*]} and use Equation (74) in (84) to get

<sup>315</sup> An IMEX Method for the Euler Equations That Posses Strong

*<sup>L</sup>*(*En*) + <sup>Δ</sup>*<sup>t</sup>*

Making use of the Taylor series given in Equations (73), (74), (75), and Equation (87) in the

*∂*2*κ<sup>n</sup> <sup>∂</sup>t*<sup>2</sup> ) *<sup>∂</sup> ∂x*

Cancelling the opposite sign common terms, grouping the other terms together, and making

2 2 *∂ ∂t* [ *∂E<sup>n</sup>*

*τ<sup>n</sup>* = *O*(Δ*t*

Here, we numerically verify our analytical findings about the two IMEX approaches. We solve the point explosion problem studied in (Kadioglu & Knoll, 2010) by using the LERH model and *M* = 200 cell points until the final time *t* = 0.02. We note that we ran the code twice as longer final time than the original test in order to allow the numerical methods to depict more accurate time behaviors. In Fig. 10, we plot the *L*2-norm of errors for variety of flow variables committed by the both approaches. Fig. 10 clearly shows that the classic IMEX method suffers from order reductions as predicted by our mathematical analysis. We present more detailed analysis regarding more general IMEX methods (e.g., Strang splitting type methods (Knoth & Wolke, 1999; Strang, 1968)) in our forthcoming paper (Kadioglu, Knoll & Lowrie,

2

2 2 *∂L*

*<sup>∂</sup>t*<sup>2</sup> <sup>−</sup> (*E<sup>n</sup>* <sup>+</sup> <sup>Δ</sup>*tL*(*En*) + <sup>Δ</sup>*<sup>t</sup>*

{*u*1[*cvρ*1(*T<sup>n</sup>* <sup>+</sup> <sup>Δ</sup>*<sup>t</sup>*

<sup>2</sup> + *Rρ*1(*T<sup>n</sup>* + Δ*t*

*L*(*E*1) + *O*(Δ*t*

*<sup>∂</sup><sup>t</sup>* <sup>+</sup> *<sup>O</sup>*(Δ*<sup>t</sup>*

(*T<sup>n</sup>* + Δ*t*

*<sup>∂</sup><sup>t</sup>* <sup>−</sup> *<sup>L</sup>*(*En*) <sup>−</sup> *<sup>∂</sup>*

*<sup>∂</sup><sup>t</sup>* <sup>−</sup> *<sup>L</sup>*(*En*) <sup>−</sup> *<sup>∂</sup>*

*∂T<sup>n</sup>*

*∂T<sup>n</sup>*

<sup>2</sup> *<sup>ρ</sup>*1(*u*1)<sup>2</sup> <sup>+</sup> *<sup>R</sup>ρ*1*T*1]}. Equation (86) can be further simplified

2 2

*∂T<sup>n</sup> <sup>∂</sup><sup>t</sup>* <sup>+</sup>

<sup>3</sup>). (88)

*∂x*

*<sup>∂</sup><sup>x</sup>* (*κ<sup>n</sup> <sup>∂</sup>T<sup>n</sup>*

(*κ<sup>n</sup> <sup>∂</sup>T<sup>n</sup>*

<sup>3</sup>), (90)

*<sup>∂</sup><sup>x</sup>* )] + *<sup>O</sup>*(Δ*<sup>t</sup>*

*<sup>∂</sup><sup>x</sup>* ) = 0, thus Equation

<sup>3</sup>)(89) .

*∂L*(*En*) *<sup>∂</sup><sup>t</sup>* )

> Δ*t* 2 2

*∂*2*T<sup>n</sup> <sup>∂</sup>t*<sup>2</sup> )]

<sup>2</sup>) in Equation (85) and perform few algebra to get

*<sup>∂</sup><sup>t</sup>* <sup>+</sup> *<sup>O</sup>*(Δ*<sup>t</sup>*

*<sup>∂</sup><sup>t</sup>* <sup>+</sup> *<sup>O</sup>*(Δ*<sup>t</sup>*

2))

<sup>2</sup>))]}. (85)

<sup>3</sup>), (86)

<sup>3</sup>). (87)

Then we have

$$\pi^{n} = \Delta t \left[ \frac{\partial E^{n}}{\partial t} - L(E^{n}) - \frac{\partial}{\partial x} (\kappa^{n} \frac{\partial T^{n}}{\partial x}) \right]$$

$$+ \frac{\Delta t^{2}}{2} \frac{\partial}{\partial t} [\frac{\partial E^{n}}{\partial t} - L(E^{n}) - \frac{\partial}{\partial x} (\kappa^{n} \frac{\partial T^{n}}{\partial x})]$$

$$+ \frac{\Delta t^{2}}{2} \rho^{1} u^{1} \frac{\partial T^{n}}{\partial t} (c\_{\upsilon} + R) + O(\Delta t^{3}).\tag{79}$$

From the energy equation (Equation (10)) we have *<sup>∂</sup>E<sup>n</sup> <sup>∂</sup><sup>t</sup>* <sup>−</sup> *<sup>L</sup>*(*En*) <sup>−</sup> *<sup>∂</sup> <sup>∂</sup><sup>x</sup>* (*κ<sup>n</sup> <sup>∂</sup>T<sup>n</sup> <sup>∂</sup><sup>x</sup>* ) = 0, thus Equation (79) becomes

$$
\tau'' = \frac{\Delta t^2}{2} \rho^1 u^1 \frac{\partial T^\mu}{\partial t} (c\_\upsilon + R) + O(\Delta t^3). \tag{80}
$$

This shows that the classic IMEX method carries first order terms in the resulting truncation error. This conclusion will be verified by our numerical computations (refer to Fig. 10).

#### **5.2 A self-consistent IMEX method**

We have already described how the self-consistent IMEX method operates on Equations (20)-(21) in Section 3.1. However, to be able to easily follow the analysis, we repeat the self-consistent operator splitting below.

Explicit block:

Step-1:

$$
\rho^1 = \rho^n - \Delta t \frac{\partial}{\partial \mathbf{x}} (\rho u)^n,
$$

$$
(\rho u)^1 = (\rho u)^n - \Delta t \frac{\partial}{\partial \mathbf{x}} [(\rho u^2)^n + p^n],
$$

$$
E^1 = E^n - \Delta t \frac{\partial}{\partial \mathbf{x}} \{u^n [E^n + p^n] \},
\tag{81}
$$

Step-2:

$$\begin{aligned} \rho^{n+1} &= \frac{\rho^1 + \rho^n}{2} - \frac{\Delta t}{2} \frac{\partial}{\partial x} (\rho u)^1, \\\\ (\rho u)^{n+1} &= \frac{(\rho u)^1 + (\rho u)^n}{2} - \frac{\Delta t}{2} \frac{\partial}{\partial x} [\rho^1 (u^2)^1 + R \rho^1 \mathbf{T}^{n+1}], \end{aligned}$$

$$E^\* = \frac{E^n + E^1}{2} - \frac{\Delta t}{2} \frac{\partial}{\partial x} \{u^1 [c\_v \rho^1 \mathbf{T}^{n+1} + \frac{1}{2} \rho^1 (u^1)^2 + R \rho^1 \mathbf{T}^{n+1}]\}, \tag{82}$$

Implicit block:

$$\frac{E^{n+1} - E^\*}{\Delta t} = \frac{1}{2} \frac{\partial}{\partial \mathbf{x}} (\kappa^{n+1} \frac{\partial T^{n+1}}{\partial \mathbf{x}}) + \frac{1}{2} \frac{\partial}{\partial \mathbf{x}} (\kappa^n \frac{\partial T^n}{\partial \mathbf{x}}).\tag{83}$$

Notice that the implicit temperature impacts the explicit block in this case (refer to the highlighted terms in Equation (82)). Again, we perform the modified equation analysis on the energy part of Equations (81)-(83). Substituting *E*<sup>1</sup> into *E*∗, we get

$$E^\* = E^n - \frac{\Delta t}{2} \frac{\partial}{\partial \mathbf{x}} \{ u^n [E^n + p^n] \} - \frac{\Delta t}{2} \frac{\partial}{\partial \mathbf{x}} \{ u^1 [c\_\nabla \rho^1 T^{n+1} + \frac{1}{2} \rho^1 (u^1)^2 + R \rho^1 T^{n+1}] \}. \tag{84}$$

22 Will-be-set-by-IN-TECH

*∂x*

*<sup>∂</sup><sup>t</sup>* <sup>−</sup> *<sup>L</sup>*(*En*) <sup>−</sup> *<sup>∂</sup>*

*<sup>∂</sup><sup>t</sup>* (*cv* <sup>+</sup> *<sup>R</sup>*) + *<sup>O</sup>*(Δ*<sup>t</sup>*

*<sup>∂</sup><sup>t</sup>* (*cv* <sup>+</sup> *<sup>R</sup>*) + *<sup>O</sup>*(Δ*<sup>t</sup>*

This shows that the classic IMEX method carries first order terms in the resulting truncation error. This conclusion will be verified by our numerical computations (refer to Fig. 10).

We have already described how the self-consistent IMEX method operates on Equations (20)-(21) in Section 3.1. However, to be able to easily follow the analysis, we repeat the

> *∂ ∂x*

> > *∂ ∂x*

(*ρu*)*n*,

[(*ρu*2)*<sup>n</sup>* + *pn*],

[*ρ*1(*u*2)<sup>1</sup> + *Rρ*1**Tn**+**1**],

1 2 *ρ*1(*u*1)

(*κ<sup>n</sup> <sup>∂</sup>T<sup>n</sup>*

1 2 *ρ*1(*u*1)

{*un*[*E<sup>n</sup>* <sup>+</sup> *<sup>p</sup>n*]}, (81)

<sup>2</sup> <sup>+</sup> *<sup>R</sup>ρ*1**Tn**+**1**]}, (82)

*<sup>∂</sup><sup>x</sup>* ). (83)

<sup>2</sup> <sup>+</sup> *<sup>R</sup>ρ*1*Tn*+1]}. (84)

*∂ ∂x*

{*u*1[*cvρ*1**Tn**+**<sup>1</sup>** <sup>+</sup>

*<sup>∂</sup><sup>x</sup>* ) + <sup>1</sup> 2 *∂ ∂x*

{*u*1[*cvρ*1*Tn*+<sup>1</sup> <sup>+</sup>

(*κn*+<sup>1</sup> *<sup>∂</sup>Tn*+<sup>1</sup>

Notice that the implicit temperature impacts the explicit block in this case (refer to the highlighted terms in Equation (82)). Again, we perform the modified equation analysis on

> 2 *∂ ∂x*

(*κ<sup>n</sup> <sup>∂</sup>T<sup>n</sup> <sup>∂</sup><sup>x</sup>* )]

*∂x*

(*κ<sup>n</sup> <sup>∂</sup>T<sup>n</sup> <sup>∂</sup><sup>x</sup>* )]

*<sup>∂</sup><sup>t</sup>* <sup>−</sup> *<sup>L</sup>*(*En*) <sup>−</sup> *<sup>∂</sup>*

<sup>3</sup>). (79)

<sup>3</sup>). (80)

*<sup>∂</sup><sup>x</sup>* ) = 0, thus

*<sup>∂</sup><sup>x</sup>* (*κ<sup>n</sup> <sup>∂</sup>T<sup>n</sup>*

*<sup>∂</sup><sup>t</sup>* <sup>−</sup> *<sup>L</sup>*(*En*) <sup>−</sup> *<sup>∂</sup>*

*τ<sup>n</sup>* = Δ*t* [

+ Δ*t* 2 2 *∂ ∂t* [ *∂E<sup>n</sup>*

+ Δ*t* 2 <sup>2</sup> *<sup>ρ</sup>*1*u*<sup>1</sup> *<sup>∂</sup>T<sup>n</sup>*

*<sup>τ</sup><sup>n</sup>* <sup>=</sup> <sup>Δ</sup>*<sup>t</sup>*

2 <sup>2</sup> *<sup>ρ</sup>*1*u*<sup>1</sup> *<sup>∂</sup>T<sup>n</sup>*

*<sup>ρ</sup>*<sup>1</sup> <sup>=</sup> *<sup>ρ</sup><sup>n</sup>* <sup>−</sup> <sup>Δ</sup>*<sup>t</sup>*

(*ρu*)<sup>1</sup> = (*ρu*)*<sup>n</sup>* <sup>−</sup> <sup>Δ</sup>*<sup>t</sup>*

*<sup>E</sup>*<sup>1</sup> <sup>=</sup> *<sup>E</sup><sup>n</sup>* <sup>−</sup> <sup>Δ</sup>*<sup>t</sup>*

<sup>2</sup> <sup>−</sup> <sup>Δ</sup>*<sup>t</sup>*

2 *∂ ∂x*

the energy part of Equations (81)-(83). Substituting *E*<sup>1</sup> into *E*∗, we get

{*un*[*E<sup>n</sup>* <sup>+</sup> *<sup>p</sup>n*]} − <sup>Δ</sup>*<sup>t</sup>*

2 *∂ ∂x*

<sup>2</sup> <sup>−</sup> <sup>Δ</sup>*<sup>t</sup>* 2 *∂ ∂x* (*ρu*) 1,

<sup>2</sup> <sup>−</sup> <sup>Δ</sup>*<sup>t</sup>* 2 *∂ ∂x*

From the energy equation (Equation (10)) we have *<sup>∂</sup>E<sup>n</sup>*

*∂E<sup>n</sup>*

Then we have

Equation (79) becomes

Explicit block: Step-1:

Step-2:

Implicit block:

*<sup>E</sup>*<sup>∗</sup> <sup>=</sup> *<sup>E</sup><sup>n</sup>* <sup>−</sup> <sup>Δ</sup>*<sup>t</sup>*

2 *∂ ∂x*

**5.2 A self-consistent IMEX method**

self-consistent operator splitting below.

*<sup>ρ</sup>n*+<sup>1</sup> <sup>=</sup> *<sup>ρ</sup>*<sup>1</sup> <sup>+</sup> *<sup>ρ</sup><sup>n</sup>*

(*ρu*)*n*+<sup>1</sup> <sup>=</sup> (*ρu*)<sup>1</sup> + (*ρu*)*<sup>n</sup>*

*<sup>E</sup>*<sup>∗</sup> <sup>=</sup> *<sup>E</sup><sup>n</sup>* <sup>+</sup> *<sup>E</sup>*<sup>1</sup>

*<sup>E</sup>n*+<sup>1</sup> <sup>−</sup> *<sup>E</sup>*<sup>∗</sup> <sup>Δ</sup>*<sup>t</sup>* <sup>=</sup> <sup>1</sup> We let *<sup>L</sup>*(*En*) = <sup>−</sup> *<sup>∂</sup> <sup>∂</sup><sup>x</sup>* {*un*[*E<sup>n</sup>* <sup>+</sup> *<sup>p</sup>n*]} and use Equation (74) in (84) to get

$$E^\* = E^u + \frac{\Delta t}{2} L(E^u) - \frac{\Delta t}{2} \frac{\partial}{\partial x} \{ u^1 [c\_v \rho^1 (T^u + \Delta t \frac{\partial T^u}{\partial t} + O(\Delta t^2)) \}$$

$$+ \frac{1}{2} \rho^1 (u^1)^2 + R\rho^1 (T^u + \Delta t \frac{\partial T^u}{\partial t} + O(\Delta t^2)) ].\tag{85}$$

Now, we insert *<sup>T</sup><sup>n</sup>* <sup>=</sup> *<sup>T</sup>*<sup>1</sup> <sup>−</sup> <sup>Δ</sup>*<sup>t</sup> <sup>∂</sup>T<sup>n</sup> <sup>∂</sup><sup>t</sup>* + *O*(Δ*t* <sup>2</sup>) in Equation (85) and perform few algebra to get

$$E^\* = E^\hbar + \frac{\Delta t}{2}L(E^\hbar) + \frac{\Delta t}{2}L(E^1) + O(\Delta t^3),\tag{86}$$

where *<sup>L</sup>*(*E*1) = <sup>−</sup> *<sup>∂</sup> <sup>∂</sup><sup>x</sup>* {*u*1[*cvρ*1*T*<sup>1</sup> <sup>+</sup> <sup>1</sup> <sup>2</sup> *<sup>ρ</sup>*1(*u*1)<sup>2</sup> <sup>+</sup> *<sup>R</sup>ρ*1*T*1]}. Equation (86) can be further simplified by using (71),

$$E^\* = E^\eta + \Delta t L(E^\eta) + \frac{\Delta t^2}{2} \frac{\partial L}{\partial t} + O(\Delta t^3). \tag{87}$$

Making use of the Taylor series given in Equations (73), (74), (75), and Equation (87) in the implicit discretization (83), we form the truncation error term as

$$\begin{split} \tau^{n} &= E^{n} + \Delta t \frac{\partial E^{n}}{\partial t} + \frac{\Delta t^{2}}{2} \frac{\partial^{2} E^{n}}{\partial t^{2}} - (E^{n} + \Delta t L(E^{n}) + \frac{\Delta t^{2}}{2} \frac{\partial L(E^{n})}{\partial t}) \\ &- \Delta t \frac{1}{2} \frac{\partial}{\partial \mathbf{x}} [ (\mathbf{x}^{n} + \Delta t \frac{\partial \mathbf{x}^{n}}{\partial t} + \frac{\Delta t^{2}}{2} \frac{\partial^{2} \mathbf{x}^{n}}{\partial t^{2}}) \frac{\partial}{\partial \mathbf{x}} (T^{n} + \Delta t \frac{\partial T^{n}}{\partial t} + \frac{\Delta t^{2}}{2} \frac{\partial^{2} T^{n}}{\partial t^{2}})] \\ &- \Delta t \frac{1}{2} \frac{\partial}{\partial \mathbf{x}} (\mathbf{x}^{n} \frac{\partial T^{n}}{\partial \mathbf{x}}) + O(\Delta t^{3}). \tag{88} \end{split} \tag{89}$$

Cancelling the opposite sign common terms, grouping the other terms together, and making use of Equation (78), we get

$$\tau^{\mathfrak{n}} = \Delta t \left[ \frac{\partial E^{\mathfrak{n}}}{\partial t} - L(E^{\mathfrak{n}}) - \frac{\partial}{\partial \mathbf{x}} (\mathbf{x}^{\mathfrak{n}} \frac{\partial T^{\mathfrak{n}}}{\partial \mathbf{x}}) \right] + \frac{\Delta t^{2}}{2} \frac{\partial}{\partial t} [\frac{\partial E^{\mathfrak{n}}}{\partial t} - L(E^{\mathfrak{n}}) - \frac{\partial}{\partial \mathbf{x}} (\mathbf{x}^{\mathfrak{n}} \frac{\partial T^{\mathfrak{n}}}{\partial \mathbf{x}})] + O(\Delta t^{3}) (89)$$

Again from the energy equation, we know that *<sup>∂</sup>E<sup>n</sup> <sup>∂</sup><sup>t</sup>* <sup>−</sup> *<sup>L</sup>*(*En*) <sup>−</sup> *<sup>∂</sup> <sup>∂</sup><sup>x</sup>* (*κ<sup>n</sup> <sup>∂</sup>T<sup>n</sup> <sup>∂</sup><sup>x</sup>* ) = 0, thus Equation (89) becomes

$$
\pi^n = O(\Delta t^3),
\tag{90}
$$

clearly proving that the self-consistent IMEX method is second order.

Here, we numerically verify our analytical findings about the two IMEX approaches. We solve the point explosion problem studied in (Kadioglu & Knoll, 2010) by using the LERH model and *M* = 200 cell points until the final time *t* = 0.02. We note that we ran the code twice as longer final time than the original test in order to allow the numerical methods to depict more accurate time behaviors. In Fig. 10, we plot the *L*2-norm of errors for variety of flow variables committed by the both approaches. Fig. 10 clearly shows that the classic IMEX method suffers from order reductions as predicted by our mathematical analysis. We present more detailed analysis regarding more general IMEX methods (e.g., Strang splitting type methods (Knoth & Wolke, 1999; Strang, 1968)) in our forthcoming paper (Kadioglu, Knoll & Lowrie, 2010).

**7. Acknowledgement**

Bartlett, Boston.

14: 43301.

**8. References**

The submitted manuscript has been authored by a contractor of the U.S. Government under Contract No. DEAC07-05ID14517 (INL/MIS-11-22498). Accordingly, the U.S. Government retains a non-exclusive, royalty-free license to publish or reproduce the published form of this

<sup>317</sup> An IMEX Method for the Euler Equations That Posses Strong

Ascher, U. M., Ruuth, S. J. & Spiteri, R. J. (1997). Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations., *Appl. Num. Math.* 25: 151–167. Ascher, U. M., Ruuth, S. J. & Wetton, B. T. R. (1995). Implicit-explicit methods for

Bates, J. W., Knoll, D. A., Rider, W. J., Lowrie, R. B. & Mousseau, V. A. (2001). On consistent

Bowers, R. & Wilson, J. (1991). *Numerical Modeling in Applied Physics and Astrophysics.*, Jones &

Brown, P. & Saad, Y. (1990). Hybrid Krylov methods for nonlinear systems of equations.,

Dai, W. & Woodward, P. (1998). Numerical simulations for radiation hydrodynamics. i.

Drake, R. (2007). Theory of radiative shocks in optically thick media., *Physics of Plasmas*

Ensman, L. (1994). Test problems for radiation and radiation hydrodynamics codes., *The*

Gottlieb, S. & Shu, C. (1998). Total variation diminishing Runge-Kutta schemes., *Mathematics*

Gottlieb, S., Shu, C. & Tadmor, E. (2001). Strong stability-preserving high-order time

Kadioglu, S. & Knoll, D. (2010). A fully second order Implicit/Explicit time integration

Kadioglu, S. & Knoll, D. (2011). Solving the incompressible Navier-Stokes equations by a

Kadioglu, S., Knoll, D. & Lowrie, R. (2010). Analysis of a self-consistent IMEX method for tightly coupled non-linear systems., *SIAM J. Sci. Comput.* to be submitted. Kadioglu, S., Knoll, D., Lowrie, R. & Rauenzahn, R. (2010). A second order self-consistent IMEX method for Radiation Hydrodynamics., *J. Comput. Phys.* 229-22: 8313–8332. Kadioglu, S., Knoll, D. & Oliveria, C. (2009). Multi-physics analysis of spherical fast burst

Kadioglu, S., Knoll, D., Sussman, M. & Martineau, R. (2010). A second order JFNK-based

IMEX method for single and multi-phase flows., *Computational Fluid Dynamics,*

technique for hydrodynamics plus nonlinear heat conduction problems., *J. Comput.*

second order IMEX method with a simple and effective preconditioning strategy.,

time-integration methods for radiation hydrodynamics in the equilibrium diffusion

contribution, or allow others to do so, for U.S. Government purposes.

Non-Linear Heat Conduction and Stiff Source Terms (Radiation Hydrodynamics)

time-dependent pde's., *SIAM J. Num. Anal.* 32: 797–823.

Castor, J. (2006). *Radiation Hydrodynamics.*, Cambridge University Press.

discretization methods., *Siam Review* 43-1: 89–112.

reactors., *Nuclear Science and Engineering* 163: 1–12.

*Springer-Verlag* DOI 10.1007/978 − 3 − 642 − 17884 − 9−69.

Dembo, R. (1982). Inexact newton methods., *SIAM J. Num. Anal.* 19: 400–408.

limit: Low energy-density regime., *J. Comput. Phys.* 167: 99–130.

Anderson, J. (1990). *Modern Compressible Flow.*, Mc Graw Hill.

*SIAM J. Sci. Stat. Comput.* 11: 450–481.

diffusion limit., *J. Comput. Phys.* 142: 182.

*Astrophysical Journal* 424: 275–291.

*Appl. Math. Comput. Sci.* Accepted.

*of Computation* 221: 73–85.

*Phys.* 229-9: 3237–3249.

Fig. 10. The self-consistent IMEX method versus a classic IMEX method in terms of the time convergence.

#### **6. Conclusion**

We have presented a self-consistent implicit/explicit (IMEX) time integration technique for solving the Euler equations that posses strong nonlinear heat conduction and very stiff source terms (Radiation hydrodynamics). The key to successfully implement an implicit/explicit algorithm in a self-consistent sense is to carry out the explicit integrations as part of the non-linear function evaluations within the implicit solver. In this way, the improved time accuracy of the non-linear iterations is immediately felt by the explicit algorithm block and the more accurate explicit solutions are readily available to form the next set of non-linear residuals. We have solved several test problems that use both of the low and high energy density radiation hydrodynamics models (the LERH and HERH models) in order to validate the numerical order of accuracy of our scheme. For each test, we have established second order time convergence. We have also presented a mathematical analysis that reveals the analytical behavior of our method and compares it to a classic IMEX approach. Our analytical findings have been supported/verified by a set of computational results. Currently, we are exploring more about our multi-phase IMEX study to solve multi-phase flow systems that posses tight non-linear coupling between the interface and fluid dynamics.
