**2. Applications**

2 Will-be-set-by-IN-TECH

Here, for illustrative purposes we used only first order time differencing. In literature, although the both algorithm blocks are formally written as second order time discretizations, the classic IMEX methods (Ascher et al., 1997; 1995; Bates et al., 2001; Kim & Moin, 1985; Lowrie et al., 1999; Ruuth, 1995) split the operators in such a way that the implicit and explicit blocks are executed independent of each other resulting in non-converged non-linearities therefore time inaccuracies (order reduction to first order is often reported for certain applications). Below, we illustrate the interaction of an explicit and an implicit algorithm

*u*<sup>1</sup> = *u<sup>n</sup>* + Δ*t f*(*un*)

Notice that the explicit block is based on a second order TVD Runge-Kutta method and the implicit block uses the Crank-Nicolson method (Gottlieb & Shu, 1998; LeVeque, 1998; Thomas, 1999). The major drawback of this strategy as mentioned above is that it does not preserve the formal second order time accuracy of the whole algorithm due to the absence of sufficient interactions between the two algorithm blocks (refer to highlighted terms in Equation (4))

In an alternative IMEX approach that we have studied extensively in (Kadioglu & Knoll, 2010; 2011; Kadioglu, Knoll & Lowrie, 2010; Kadioglu, Knoll, Lowrie & Rauenzahn, 2010; Kadioglu et al., 2009), the explicit block is always solved inside the implicit block as part of the nonlinear function evaluation making use of the well-known Jacobian-Free Newton Krylov (JFNK) method (Brown & Saad, 1990; Knoll & Keyes, 2004). We refer this IMEX approach as *a self-consistent IMEX method*. In this strategy, there is a continuous interaction between the implicit and explicit blocks meaning that the improved solutions (in terms of time accuracy) at each nonlinear iteration are immediately felt by the explicit block and the improved explicit solutions are readily available to form the next set of nonlinear residuals. This continuous interaction between the two algorithm blocks results in an implicitly balanced algorithm in that all nonlinearities due to coupling of different time terms are consistently converged. In other words, we obtain an IMEX method that eliminates potential order reductions in time accuracy (the formal second order time accuracy of the whole algorithm is preserved). Below, we illustrate the interaction of the explicit and implicit blocks of the self-consistent IMEX method for the scalar model in Equation (1). The interaction occurs through the highlighted

*u*<sup>1</sup> = *u<sup>n</sup>* + Δ*t f*(*un*)

*u*<sup>∗</sup> = (*u*<sup>1</sup> + *un*)/2 + Δ*t*/2 *f*(*un*+1) (6)

*un*+<sup>1</sup> = *u*<sup>∗</sup> + Δ*t*/2[*g*(*un*) + *g*(*un*+1)]. (7)

*u*<sup>∗</sup> = (*u*<sup>1</sup> + *un*)/2 + Δ*t*/2 *f*(*u*1) (4)

*un*+<sup>1</sup> = *u*<sup>∗</sup> + Δ*t*/2[*g*(*un*) + *g*(*un*+1)]. (5)

<sup>Δ</sup>*<sup>t</sup>* <sup>=</sup> *<sup>g</sup>*(*un*+1). (3)

*<sup>u</sup>n*+<sup>1</sup> <sup>−</sup> *<sup>u</sup>*<sup>∗</sup>

block based on second order time discretizations of Equation(1) in classical sense,

(Bates et al., 2001; Kadioglu, Knoll & Lowrie, 2010).

Implicit block solves:

Explicit block:

Implicit block:

terms in Equation (6). Explicit block:

Implicit block:

We have applied the above described self-consistent IMEX method to both multi-physics and multiple time scale fluid dynamics problems (Kadioglu & Knoll, 2010; 2011; Kadioglu, Knoll, Lowrie & Rauenzahn, 2010; Kadioglu et al., 2009; Kadioglu, Knoll, Sussman & Martineau, 2010). The multi-physics application comes from a multi-physics analysis of fast burst reactor study (Kadioglu et al., 2009). The model couples a neutron dynamics that simulates the transient behavior of neutron populations to a mechanics model that predicts material expansions and contractions. It is important to introduce a second order accurate numerical procedure for this kind of nonlinearly coupled system, because the criticality and safety study can depend on how well we predict the feedback between the neutronics and the mechanics of the fuel assembly inside the reactor. In our second order self-consistent IMEX framework, the mechanics part is solved explicitly inside the implicit neutron diffusion block as part of the nonlinear function evaluation. We have reported fully second order time convergent calculations for this model (Kadioglu et al., 2009).

As part of the multi-scale fluid dynamics application, we have solved multi-phase flow problems which are modeled by incompressible two-phase Navier-Stokes equations that govern the flow dynamics plus a level set equation that solves the inter-facial dynamics between the fluids (Kadioglu, Knoll, Sussman & Martineau, 2010). In these kinds of models, there is a strong non-linear coupling between the interface and fluid dynamics, e.g, the viscosity coefficient and surface tension forces are highly non-linear functions of interface variables, on the other hand, the fluid interfaces are advected by the flow velocity. Therefore, it is important to introduce an accurate integration technique that converges all non-linearities due to the strong coupling. Our self-consistent IMEX method operates on this model as follows; the interface equation together with the hyperbolic parts of the fluid equations are treated explicitly and solved inside an implicit loop that solves the viscous plus stiff surface tension forces. More details about the splitting of the operators of the Navier-Stokes equations in a self-consistent IMEX manner can be found in (Kadioglu & Knoll, 2011).

Another multi-scale fluid dynamics application comes from radiation hydrodynamics that we will be focusing on in the remainder of this chapter. Radiation hydrodynamics models are commonly used in astrophysics, inertial confinement fusion, and other high-temperature flow systems (Bates et al., 2001; Castor, 2006; Dai & Woodward, 1998; Drake, 2007; Ensman, 1994; Kadioglu & Knoll, 2010; Lowrie & Edwards, 2008; Lowrie & Rauenzahn, 2007; Mihalas & Mihalas, 1984; Pomraning, 1973). A commonly used model considers the compressible Euler equations that contains a non-linear heat conduction term in the energy part. This model is relatively simple and often referred to as a *Low Energy-Density Radiation Hydrodynamics (LERH)* in a diffusion approximation limit (Kadioglu & Knoll, 2010). A more complicated model is referred to as a *High Energy-Density Radiation Hydrodynamics (HERH)* in a diffusion approximation limit that considers a combination of a hydrodynamical model resembling the compressible Euler equations and a radiation energy model that contains a separate radiation energy equation with nonlinear diffusion plus coupling source terms to

temperatures the radiation effects are negligible, therefore, a low energy density model (LERH) that limits the radiation effects to a non-linear heat conduction is sufficient. However, at high temperatures, a more complicated high energy density radiation hydrodynamics (HERH) model that accounts for more significant radiation effects has to be considered. Accordingly, the governing equations of the HERH model consist of the following system

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

(*r*2*ρu*2) + *<sup>∂</sup>*

(*r*2*cDr*

where the flow variables and parameters that also occur in the LERH model are described above. Here, more variable definitions come from the radiation physics, i.e, *E<sup>ν</sup>* is the

Stephan-Boltzmann constant, *σa* is the macroscopic absorption cross-section, and *Dr* is the

*Dr*(*T*) = <sup>1</sup>

We note that we solve a non-dimensional version of Equations (11)-(14) in order to normalize large digit numbers (*c*, *σa*, *a* etc.) and therefore improve the performance of the non-linear solver. The details of the non-dimensionalization procedure are given in (Kadioglu, Knoll, Lowrie & Rauenzahn, 2010). The non-dimensional system is the following,

<sup>2</sup>*ρu*2) + *<sup>∂</sup>*

[*r*2*u*(*<sup>E</sup>* <sup>+</sup> *<sup>p</sup>*)] = −P*σa*(*T*<sup>4</sup> <sup>−</sup> *<sup>E</sup>ν*) <sup>−</sup> <sup>1</sup>

Here, we present the numerical procedure for the LERH model. The extension to the HERH model is straight forward. First, we split the operators of Equations (8)-(10) into two pieces one being the pure hydrodynamics part (hyperbolic conservation laws) and the other

*∂r*

is a non-dimensional parameter that measures the radiation effects on the

3*σa*

radiation diffusion coefficient. From the simple diffusion theory, *Dr* can be written as

1 *r*2 *∂ ∂r*

*r*2 *∂ ∂r* (*r* 2*κ ∂E<sup>ν</sup>*

flow and is roughly proportional to the ratio of the radiation and fluid pressures.

*∂ρ <sup>∂</sup><sup>t</sup>* <sup>+</sup>

*∂ ∂t*

1 *r*2 *∂ ∂r*

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

1 *r*2 *∂ ∂r*

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

**3. Numerical procedure**

0 *ρ*<sup>0</sup> *c*<sup>2</sup> *s*,0

where <sup>P</sup> <sup>=</sup> *aT*<sup>4</sup>

(*ρu*) + <sup>1</sup> *r*2 *∂ ∂r* (*r*

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

*∂r*

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

<sup>2</sup>*u*(*<sup>E</sup>* <sup>+</sup> *<sup>p</sup>*)] = <sup>−</sup>*cσa*(*aT*<sup>4</sup> <sup>−</sup> *<sup>E</sup>ν*) <sup>−</sup> <sup>1</sup>

(*r*2*ρu*) = 0, (11)

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

*<sup>∂</sup><sup>r</sup>* ) + *<sup>c</sup>σa*(*aT*<sup>4</sup> <sup>−</sup> *<sup>E</sup>ν*) + <sup>1</sup>

<sup>3</sup> is the radiation pressure, *c* is the speed of light, *a* is the

(*p* + *pν*) = 0, (12)

. (15)

(*r*2*ρu*) = 0, (16)

3 P*u ∂E<sup>ν</sup>*

*<sup>∂</sup><sup>r</sup>* ) + *<sup>σ</sup>a*(*T*<sup>4</sup> <sup>−</sup> *<sup>E</sup>ν*) + <sup>1</sup>

(*p* + P*pν*) = 0, (17)

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

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

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

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

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

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

*∂ρ <sup>∂</sup><sup>t</sup>* <sup>+</sup>

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

(*ρu*) + <sup>1</sup> *r*2 *∂ ∂r*

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

*∂ ∂t*

1 *r*2 *∂ ∂r* [*r*

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

1 *r*2 *∂ ∂r*

radiation energy density, *p<sup>ν</sup>* = *<sup>E</sup><sup>ν</sup>*

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

1 *r*2 *∂ ∂r*

*r*2 *∂ ∂r*

materials (Kadioglu, Knoll, Lowrie & Rauenzahn, 2010). Radiation Hydrodynamics problems are difficult to tackle numerically since they exhibit multiple time scales. For instance, radiation and hydrodynamics process can occur on time scales that can differ from each other by many orders of magnitudes. Hybrid methods (Implicit/Explicit (IMEX) methods) are highly desirable for these kinds of models, because if one uses all explicit discretizations, then due to very stiff diffusion process the explicit time steps become often impractically small to satisfy stability conditions (LeVeque, 1998; Thomas, 1999). Previous IMEX attempts to solve these problems were not quite successful, since they often reported order reductions in time accuracy (Bates et al., 2001; Lowrie et al., 1999). The main reason for time inaccuracies was how the explicit and implicit operators were split in which explicit solutions were lagging behind the implicit ones. In our self-consistent IMEX method, the hydrodynamics part is solved explicitly making use of the well-understood explicit schemes within an implicit diffusion block that corresponds to radiation transport. Explicit solutions are obtained as part of the non-linear functions evaluations withing the JFNK framework. This strategy has enabled us to produce fully second order time accurate results for both LERH and more complicated HERH models (Kadioglu & Knoll, 2010; Kadioglu, Knoll, Lowrie & Rauenzahn, 2010).

In the following sections, we will go over more details about the LERH and HERH models and the implementation/implications of the self-consistent IMEX technology when it is applied to these models. We will also present a mathematical analysis that reveals the analytical convergence behavior of our method and compares it to a classic IMEX approach.

#### **2.1 A Low Energy Density Radiation Hydrodynamics Model (LERH)**

This model uses the following system of partial differential equations formulated in spherically symmetric coordinates.

$$
\frac{\partial \rho}{\partial t} + \frac{1}{r^2} \frac{\partial}{\partial r} (r^2 \rho u) = 0,\tag{8}
$$

$$\frac{\partial}{\partial t}(\rho u) + \frac{1}{r^2} \frac{\partial}{\partial r}(r^2 \rho u^2) + \frac{\partial p}{\partial r} = 0,\tag{9}$$

$$\frac{\partial E}{\partial t} + \frac{1}{r^2} \frac{\partial}{\partial r} [r^2 \mu (E + p)] = \frac{1}{r^2} \frac{\partial}{\partial r} (r^2 \kappa \frac{\partial T}{\partial r}) ,\tag{10}$$

where *ρ*, *u*, *p*, *E*, and *T* are the mass density, flow velocity, fluid pressure, total energy density of the fluid, and the fluid temperature respectively. *κ* is the coefficient of thermal conduction (or diffusion coefficient) and in general is a nonlinear function of *ρ* and *T*. In this study, we will use an ideal gas equation of state, i.e, *p* = *RρT* = (*γ* − 1)*ρ�*, where *R* is the specific gas constant per unit mass, *γ* is the ratio of specific heats, and *�* is the internal energy of the fluid per unit mass. The coefficient of thermal conduction will be assumed to be written as a power law in density and temperature, i.e, *κ* = *κ*0*ρaTb*, where *κ*0, *a* and *b* are constants (Marshak, 1958). This simplified radiation hydrodynamics model allows one to study the dynamics of nonlinearly coupled two distinct physics; compressible fluid flow and nonlinear diffusion.

#### **2.2 A High Energy Density Radiation Hydrodynamics Model (HERH)**

In general, the radiation hydrodynamics concerns the propagation of thermal radiation through a fluid and the effect of this radiation on the hydrodynamics describing the fluid motion. The role of the thermal radiation increases as the temperature is raised. At low 4 Will-be-set-by-IN-TECH

materials (Kadioglu, Knoll, Lowrie & Rauenzahn, 2010). Radiation Hydrodynamics problems are difficult to tackle numerically since they exhibit multiple time scales. For instance, radiation and hydrodynamics process can occur on time scales that can differ from each other by many orders of magnitudes. Hybrid methods (Implicit/Explicit (IMEX) methods) are highly desirable for these kinds of models, because if one uses all explicit discretizations, then due to very stiff diffusion process the explicit time steps become often impractically small to satisfy stability conditions (LeVeque, 1998; Thomas, 1999). Previous IMEX attempts to solve these problems were not quite successful, since they often reported order reductions in time accuracy (Bates et al., 2001; Lowrie et al., 1999). The main reason for time inaccuracies was how the explicit and implicit operators were split in which explicit solutions were lagging behind the implicit ones. In our self-consistent IMEX method, the hydrodynamics part is solved explicitly making use of the well-understood explicit schemes within an implicit diffusion block that corresponds to radiation transport. Explicit solutions are obtained as part of the non-linear functions evaluations withing the JFNK framework. This strategy has enabled us to produce fully second order time accurate results for both LERH and more complicated HERH models (Kadioglu & Knoll, 2010; Kadioglu, Knoll, Lowrie & Rauenzahn,

In the following sections, we will go over more details about the LERH and HERH models and the implementation/implications of the self-consistent IMEX technology when it is applied to these models. We will also present a mathematical analysis that reveals the analytical

This model uses the following system of partial differential equations formulated in

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

where *ρ*, *u*, *p*, *E*, and *T* are the mass density, flow velocity, fluid pressure, total energy density of the fluid, and the fluid temperature respectively. *κ* is the coefficient of thermal conduction (or diffusion coefficient) and in general is a nonlinear function of *ρ* and *T*. In this study, we will use an ideal gas equation of state, i.e, *p* = *RρT* = (*γ* − 1)*ρ�*, where *R* is the specific gas constant per unit mass, *γ* is the ratio of specific heats, and *�* is the internal energy of the fluid per unit mass. The coefficient of thermal conduction will be assumed to be written as a power law in density and temperature, i.e, *κ* = *κ*0*ρaTb*, where *κ*0, *a* and *b* are constants (Marshak, 1958). This simplified radiation hydrodynamics model allows one to study the dynamics of nonlinearly coupled two distinct physics; compressible fluid flow and nonlinear diffusion.

In general, the radiation hydrodynamics concerns the propagation of thermal radiation through a fluid and the effect of this radiation on the hydrodynamics describing the fluid motion. The role of the thermal radiation increases as the temperature is raised. At low

<sup>2</sup>*ρu*2) + *<sup>∂</sup><sup>p</sup>*

*r*2 *∂ ∂r* (*r*2*κ ∂T*

(*r*2*ρu*) = 0, (8)

*<sup>∂</sup><sup>r</sup>* <sup>=</sup> 0, (9)

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

convergence behavior of our method and compares it to a classic IMEX approach.

1 *r*2 *∂ ∂r*

**2.1 A Low Energy Density Radiation Hydrodynamics Model (LERH)**

*∂ ∂t*

> 1 *r*2 *∂ ∂r*

**2.2 A High Energy Density Radiation Hydrodynamics Model (HERH)**

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

(*ρu*) + <sup>1</sup> *r*2 *∂ ∂r* (*r*

spherically symmetric coordinates.

2010).

temperatures the radiation effects are negligible, therefore, a low energy density model (LERH) that limits the radiation effects to a non-linear heat conduction is sufficient. However, at high temperatures, a more complicated high energy density radiation hydrodynamics (HERH) model that accounts for more significant radiation effects has to be considered. Accordingly, the governing equations of the HERH model consist of the following system

$$\frac{\partial \rho}{\partial t} + \frac{1}{r^2} \frac{\partial}{\partial r} (r^2 \rho u) = 0,\tag{11}$$

$$\frac{\partial}{\partial t}(\rho u) + \frac{1}{r^2}\frac{\partial}{\partial r}(r^2\rho u^2) + \frac{\partial}{\partial r}(p + p\_\nu) = 0,\tag{12}$$

$$\frac{\partial E}{\partial t} + \frac{1}{r^2} \frac{\partial}{\partial r} [r^2 u(E+p)] = -c\sigma\_a (aT^4 - E\_\nu) - \frac{1}{3} u \frac{\partial E\_\nu}{\partial r} \tag{13}$$

$$\frac{\partial E\_{\nu}}{\partial t} + \frac{1}{r^{2}} \frac{\partial}{\partial r} [r^{2} u (E\_{\nu} + p\_{\nu})] = \frac{1}{r^{2}} \frac{\partial}{\partial r} (r^{2} c D\_{r} \frac{\partial E\_{\nu}}{\partial r}) + c \sigma\_{a} (a T^{4} - E\_{\nu}) + \frac{1}{3} u \frac{\partial E\_{\nu}}{\partial r},\tag{14}$$

where the flow variables and parameters that also occur in the LERH model are described above. Here, more variable definitions come from the radiation physics, i.e, *E<sup>ν</sup>* is the radiation energy density, *p<sup>ν</sup>* = *<sup>E</sup><sup>ν</sup>* <sup>3</sup> is the radiation pressure, *c* is the speed of light, *a* is the Stephan-Boltzmann constant, *σa* is the macroscopic absorption cross-section, and *Dr* is the radiation diffusion coefficient. From the simple diffusion theory, *Dr* can be written as

$$D\_{\mathbf{f}}(T) = \frac{1}{3\sigma\_{\mathbf{a}}}.\tag{15}$$

We note that we solve a non-dimensional version of Equations (11)-(14) in order to normalize large digit numbers (*c*, *σa*, *a* etc.) and therefore improve the performance of the non-linear solver. The details of the non-dimensionalization procedure are given in (Kadioglu, Knoll, Lowrie & Rauenzahn, 2010). The non-dimensional system is the following,

$$\frac{\partial \rho}{\partial t} + \frac{1}{r^2} \frac{\partial}{\partial r} (r^2 \rho u) = 0,\tag{16}$$

$$\frac{\partial}{\partial t}(\rho u) + \frac{1}{r^2} \frac{\partial}{\partial r}(r^2 \rho u^2) + \frac{\partial}{\partial r}(p + \mathcal{P}p\_\nu) = 0,\tag{17}$$

$$\frac{\partial E}{\partial t} + \frac{1}{r^2} \frac{\partial}{\partial r} [r^2 u(E + p)] = -\mathcal{P} \sigma\_\theta (T^4 - E\_\nu) - \frac{1}{3} \mathcal{P} u \frac{\partial E\_\nu}{\partial r} \,. \tag{18}$$

$$\frac{\partial E\_{\nu}}{\partial t} + \frac{1}{r^{2}} \frac{\partial}{\partial r} [r^{2} u (E\_{\nu} + p\_{\nu})] = \frac{1}{r^{2}} \frac{\partial}{\partial r} (r^{2} \kappa \frac{\partial E\_{\nu}}{\partial r}) + \sigma\_{a} (T^{4} - E\_{\nu}) + \frac{1}{3} u \frac{\partial E\_{\nu}}{\partial r},\tag{19}$$

where <sup>P</sup> <sup>=</sup> *aT*<sup>4</sup> 0 *ρ*<sup>0</sup> *c*<sup>2</sup> *s*,0 is a non-dimensional parameter that measures the radiation effects on the flow and is roughly proportional to the ratio of the radiation and fluid pressures.

#### **3. Numerical procedure**

Here, we present the numerical procedure for the LERH model. The extension to the HERH model is straight forward. First, we split the operators of Equations (8)-(10) into two pieces one being the pure hydrodynamics part (hyperbolic conservation laws) and the other

execution of the whole algorithm in the self-consistent IMEX sense (refer to Fig. 1). According to this diagram, at beginning of each Newton iteration, we have the temperature values based on the current Newton iterate. This temperature is passed to the explicit block that returns the updated density, momentum, and a prediction to total energy. Then we form the non-linear residuals (e.g, forming the IMEX function in Section 3.3) for the diffusion equation out of the updated and predicted values. With the IMEX function in hand, we can execute the JFNK method. After the Newton method convergences, we get second order converged temperature

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

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

Our explicit time discretization is based on a second order TVD Runge-Kutta method (Gottlieb & Shu, 1998; Gottlieb et al., 2001; Shu & Osher, 1988; 1989). The main reason why we choose this methodology is that it preserves the strong stability properties of the explicit Euler method. This is important because it is well known that solutions to the conservation laws usually involve discontinuities (e.g, shock or contact discontinuities) and (Gottlieb & Shu, 1998; Gottlieb et al., 2001) suggest that a time integration method which has the strong stability preserving property leads to non-oscillatory calculations (especially at shock or

A second order two-step TVD Runge-Kutta method for (20) can be cast as

1 *r*2 *∂ ∂r* (*r* <sup>2</sup>*ρu*)*n*,

> *r*2 *∂ ∂r*

2 { 1 *r*2 *∂ ∂r* (*r* <sup>2</sup>*ρu*2) <sup>1</sup> + *∂ ∂r*

*p* = *ρRTE* = *cvρT* +

explicit algorithm block interacts with the implicit block through the highlighted *Tn*+<sup>1</sup> terms in Equation (23). We can observe that the implicit equation (21) is practically solved for *T* by using the energy relation. Therefore, the explicit block is continuously impacted by the

<sup>2</sup>*u*1(*cvρ*1**Tn**+**<sup>1</sup>** +

1 *r*2 *∂ ∂r* (*r*2*ρu*2) + *<sup>∂</sup><sup>p</sup>*

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

1 2

*<sup>γ</sup>*−<sup>1</sup> is the fluid specific heat with *<sup>R</sup>* being the universal gas constant. This

[*r*2*u*(*<sup>E</sup>* <sup>+</sup> *<sup>p</sup>*)]}*n*,

*∂r* ] *n*,

(*ρ*1*R***Tn**+**1**)},

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

*ρu*2, (24)

(22)

(23)

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

*<sup>ρ</sup>*<sup>1</sup> <sup>=</sup> *<sup>ρ</sup><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>*t*{ <sup>1</sup>

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

2 { 1 *r*2 *∂ ∂r* [*r*

We used the following equation of state relations in (22)- (23);

(*ρu*)

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

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

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

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

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

and total energy density field.

**3.1 Explicit block**

contact discontinuities).

Step-1 :

Step-2 :

where *cv* = *<sup>R</sup>*

ρ ρu E (n) i i = 1,2,...,N t n For k =1,...,kmax Call Hydrodynamics Block with T Form Non−Linear Residual )/2 n+1 <sup>2</sup> (u n+1 <sup>+</sup> <sup>ρ</sup> k ρ T n+1 cv − E\* Δ t Res = , , Explicit Hydrodynamics Block Based on Second order R−K method Input T Return ( ρ n+1 , (ρu ) n+1, E\* ) Newton Iteration E n+1= cv ρn+1T k+1+ <sup>ρ</sup>n+1(u n+1 ) /2 2 End δ k T k+1= Tk+ <sup>δ</sup> <sup>T</sup> k −( RHS( ) + RHS( ))/2 Calculate T n+1 t E ρ u ρ (n+1) i i = 1,2,...,N T is available k k k k T n+1 <sup>ρ</sup> n T n ρ

Fig. 1. Flowchart of the second order self-consistent IMEX algorithm

accounting for the effects of radiation transport (diffusion equation). For instance, the pure hydrodynamics equations can be written as

$$
\frac{
\partial \mathbf{U}
}{
\partial t
} + \frac{
\partial (A\mathbf{F})
}{
\partial V
} + \frac{
\partial \mathbf{G}
}{
\partial r
} = 0,
\tag{20}
$$

where **U** = (*ρ*, *ρu*, *E*)*T*, **F**(**U**)=(*ρu*, *ρu*2, *u*(*E* + *p*))*T*, and **G**(**U**)=(0, *p*, 0)*T*. Then the diffusion equation becomes

$$\frac{\partial E}{\partial t} = \frac{\partial}{\partial V} (A\kappa \frac{\partial T}{\partial r}) , \tag{21}$$

where *V* = <sup>4</sup> <sup>3</sup>*πr*<sup>3</sup> is the generalized volume coordinate in one-dimensional spherical geometry, and *A* = 4*πr*<sup>2</sup> is the associated cross-sectional area. Notice that the total energy density, *E*, obtained by Equation (20) just represents the hydrodynamics component and it must be augmented by Equation (21).

Our algorithm consists of an explicit and an implicit block. The explicit block solves Equation (20) and the implicit block solves Equation (21). We will briefly describe these algorithm blocks in the following subsections. However, we note again that the explicit block is embedded within the implicit block as part of a nonlinear function evaluation as it is depicted in Fig. 1. This is done to obtain a nonlinearly converged algorithm that leads to second order calculations. We also note that similar discretizations, but without converging nonlinearities, can lead to order reduction in time convergence (Bates et al., 2001). Before we go into details of the individual algorithm blocks, we would like to present a flow diagram that illustrates the

execution of the whole algorithm in the self-consistent IMEX sense (refer to Fig. 1). According to this diagram, at beginning of each Newton iteration, we have the temperature values based on the current Newton iterate. This temperature is passed to the explicit block that returns the updated density, momentum, and a prediction to total energy. Then we form the non-linear residuals (e.g, forming the IMEX function in Section 3.3) for the diffusion equation out of the updated and predicted values. With the IMEX function in hand, we can execute the JFNK method. After the Newton method convergences, we get second order converged temperature and total energy density field.
