An IMEX Method for the Euler Equations That Posses Strong Non-Linear Heat Conduction and Stiff Source Terms (Radiation Hydrodynamics)

Here, we present a truly second order time accurate self-consistent IMEX (IMplicit/EXplicit) method for solving the Euler equations that posses strong nonlinear heat conduction and very stiff source terms (Radiation hydrodynamics). This study essentially summarizes our previous and current research related to this subject (Kadioglu & Knoll, 2010; 2011; Kadioglu, Knoll & Lowrie, 2010; Kadioglu, Knoll, Lowrie & Rauenzahn, 2010; Kadioglu et al., 2009; Kadioglu, Knoll, Sussman M 1995; Bates et al., 2001; Kadioglu & Knoll, 2010; 2011; Kadioglu, Knoll, Lowrie & Rauenzahn, 2010; Kadioglu et al., 2009; Khan & Liu, 1994; Kim & Moin, 1985; Lowrie et al., 1999; Ruuth, 1995). These methods are particularly attractive when dealing with physical systems that consist of multiple physics (multi-physics problems such as coupling of neutron dynamics to thermal-hydrolic or to thermal-mechanics in reactors) or fluid dynamics problems that exhibit multiple time scales such as advection-diffusion, reaction-diffusion, or advection-diffusion-reaction problems. In general, governing equations for these kinds of systems consist of stiff and non-stiff terms. This poses numerical challenges in regards to time integrations, since most of the temporal numerical methods are designed specific for either stiff or non-stiff problems. Numerical methods that can handle both physical behaviors are often referred to as IMEX methods. A typical IMEX method isolates the stiff and non-stiff parts of the governing system and employs an explicit discretization strategy that solves the non-stiff part and an implicit technique that solves the stiff part of the problem. This standard IMEX approach can be summarized by considering a simple prototype model. Let us consider the following scalar model ut = f (u) + g(u), (1)


Introduction
Here, we present a truly second order time accurate self-consistent IMEX (IMplicit/EXplicit) method for solving the Euler equations that posses strong nonlinear heat conduction and very stiff source terms (Radiation hydrodynamics). This study essentially summarizes our previous and current research related to this subject 2011;Kadioglu, Knoll, Lowrie & Rauenzahn, 2010;Kadioglu et al., 2009;Kadioglu, Knoll, Sussman & Martineau, 2010). Implicit/Explicit (IMEX) time integration techniques are commonly used in science and engineering applications (Ascher et al., 1997;1995;Bates et al., 2001;2011;Kadioglu, Knoll, Lowrie & Rauenzahn, 2010;Kadioglu et al., 2009;Khan & Liu, 1994;Kim & Moin, 1985;Lowrie et al., 1999;Ruuth, 1995). These methods are particularly attractive when dealing with physical systems that consist of multiple physics (multi-physics problems such as coupling of neutron dynamics to thermal-hydrolic or to thermal-mechanics in reactors) or fluid dynamics problems that exhibit multiple time scales such as advection-diffusion, reaction-diffusion, or advection-diffusion-reaction problems.
In general, governing equations for these kinds of systems consist of stiff and non-stiff terms. This poses numerical challenges in regards to time integrations, since most of the temporal numerical methods are designed specific for either stiff or non-stiff problems. Numerical methods that can handle both physical behaviors are often referred to as IMEX methods. A typical IMEX method isolates the stiff and non-stiff parts of the governing system and employs an explicit discretization strategy that solves the non-stiff part and an implicit technique that solves the stiff part of the problem. This standard IMEX approach can be summarized by considering a simple prototype model. Let us consider the following scalar model where f (u) and g(u) represent non-stiff and stiff terms respectively. Then the IMEX strategy consists of the following algorithm blocks: Explicit block solves: u * − u n ∆t = f (u n ),( 2 ) Implicit block solves: u n+1 − u * ∆t = g(u n+1 ). ( 3 ) 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 block based on second order time discretizations of Equation(1) in classical sense, Explicit block: Implicit block: u n+1 = u * + ∆t/2[g(u n )+g(u n+1 )].( 5 ) 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)) (Bates et al., 2001;. In an alternative IMEX approach that we have studied extensively in 2011;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 terms in Equation (6). Explicit block: Implicit block:
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;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 . 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 295 An IMEX Method for the Euler Equations That Posses Strong Non-Linear Heat Conduction and Stiff Source Terms (Radiation Hydrodynamics) www.intechopen.com 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, 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.

A Low Energy Density Radiation Hydrodynamics Model (LERH)
This model uses the following system of partial differential equations formulated in spherically symmetric coordinates.
where ρ, u, p, E,andT 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. I nt h i ss t u d y ,w e will use an ideal gas equation of state, i.e, p = RρT =( γ − 1)ρǫ,w h e r eR 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 ρ a T b ,w h e r eκ 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.

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 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 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 ν is the radiation energy density, p ν = E ν 3 is the radiation pressure, c is the speed of light, a is the Stephan-Boltzmann constant, σ a is the macroscopic absorption cross-section, and D r is the radiation diffusion coefficient. From the simple diffusion theory, D r can be written as 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, where P = 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.

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 Explicit Hydrodynamics Block Based on Second order R−K method 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 where U =( ρ, ρu, E) T , F(U)=( ρu, ρu 2 , u(E + p)) T ,a n dG(U)=( 0, p,0) T . Then the diffusion equation becomes where V = 4 3 πr 3 is the generalized volume coordinate in one-dimensional spherical geometry, and A = 4πr 2 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.

Explicit block
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 contact discontinuities). A second order two-step TVD Runge-Kutta method for (20) can be cast as Step-1 : Step-2 : We used the following equation of state relations in (22)-(23); where c v = R γ−1 is the fluid specific heat with R being the universal gas constant. This explicit algorithm block interacts with the implicit block through the highlighted T n+1 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 implicit T n+1 solutions at each non-linear Newton iteration. This provides the tight nonlinear coupling between the two algorithm blocks. Notice that the k th nonlinear Newton iteration of the implicit block corresponds to T n+1 ← T k and k → (n + 1) upon the convergence of the Newton method (refer to Fig. 1). Also, the * values in Equation (23) are predicted intermediate values and later they are corrected by the implicit block which is given in the next subsection. One observation about this algorithm block is that some calculations are redundant related to Equation (22). In other words, Equation (22) can be computed only once at the beginning of each Newton iteration, because the non-linear iterations do not impact (22). This can lead overall less number of function evaluations. Now we shall describe how we evaluate the numerical fluxes needed by Equations (22) and (23). For simplicity, we consider (20) to describe our fluxing procedure. Basically, it is based on the Local Lax Friedrichs (LLF) method (we refer to (LeVeque, 1998;Thomas, 1999) for the details of the LLF method and for more information in regards to the explicit discretizations of conservation laws). For instance, if we consider the following simple discretization for Equation (20), where ∆V i = V(r i+1/2 ) − V(r i−1/2 ), A i±1/2 = A(r i±1/2 ),andindicesi and i + 1/2 represent cell center and cell edge values respectively (refer to Fig. 2), then the Local Lax Friedrichs method defines F i+1/2 and G i+1/2 as where α = max{|λ L 1 |, |λ R 1 |, |λ L 2 |, |λ R 2 |, |λ L 3 |, |λ R 3 |} in which λ 1 = u − c, λ 2 = u, λ 3 = u + c,and c is the sound speed. The sound speed is defined by where ∂p ∂ρ = RT in this study. U R i+1/2 and U L i+1/2 are the interpolated values at (i + 1/2) th cell edge from the right and left side, i.e, where where

Implicit block
The explicit block produces the following solution vector This information is used to discretize Equation (21) as follows, where Notice that this implicit discretization resembles to the Crank-Nicolson method (Strikwerda, 1989;Thomas, 1998). We solve Equation (33) iteratively for T n+1 . The nonlinear solver needed by Equation (33) is based on the Jacobian-Free Newton Krylov method which is described in the next subsection. When the Newton method converges all the nonlinearities in this discretization, we obtain the following fully updated solution vector,

An IMEX Method for the Euler Equations That Posses Strong Non-Linear Heat Conduction and Stiff Source Terms (Radiation Hydrodynamics)
www.intechopen.com

The Jacobian-Free Newton Krylov method and forming the IMEX function
The Jacobian-Free Newton Krylov method (e.g, refer to (Brown & Saad, 1990;Kelley, 2003;Knoll & Keyes, 2004)) is a combination of the Newton method that solves a system of nonlinear equations and a Krylov subspace method that solves the Newton correction equations. With this method, Newton-like super-linear convergence is achieved in the nonlinear iterations, without the complexity of forming or storing the Jacobian matrix. The effects of the Jacobian matrix are probed only through approximate matrix-vector products required in the Krylov iterations. Below, we provide more details about this technique. The Newton method solves F(T)=0 (e.g, assume Equation (33) is written in this form) iteratively over a sequence of linear system defined by where J(T k )= ∂F ∂T is the Jacobian matrix and δT k is the update vector. The Newton iteration is terminated based on a required drop in the norm of the nonlinear residual, i.e, where tol res is a given tolerance. The linear system, Newton correction equation (35), is solved by using the Arnoldi based Generalized Minimal RESidual method (GMRES) (Saad, 2003) which belongs to the general class of the Krylov subspace methods (Reid, 1971). We note that these subspace methods are particularly suitable choice when dealing with non-symmetric linear systems. In GMRES, an initial linear residual, r 0 , is defined for a given initial guess δT 0 , Here we dropped the index k convention since the Krylov (GMRES) iteration is performed at a fixed k.L e t j be the Krylov iteration index. The j th Krylov iteration minimizes JδT j + F(T) 2 within a subspace of small dimension, relative to n (the number of unknowns), in a least-squares sense. δT j is drawn from the subspace spanned by the Krylov vectors, {r 0 , Jr 0 , J 2 r 0 , ··· , J j−1 r 0 } ,andcanbewrittenas where the scalar β i minimizes the residual. The Krylov iteration is terminated based on the following inexact Newton criteria (Dembo, 1982) JδT where the parameter γ is set in terms of how tight the linear solver should converge at each Newton iteration (we typically use γ = 10 −3 ). One particularly attractive feature of this methodology is that it does not require forming the Jacobian matrix. Instead, only matrix-vector multiplications, Jv, are needed, where v ∈{ r 0 , Jr 0 , J 2 r 0 , ···} .T h i sl e a d s t o the so-called Jacobian-Free implementations in which the action of the Jacobian matrix can be approximated by where ǫ = 1 n v 2 ∑ n i=1 b|u i | + b, n is the dimension of the linear system and b is a constant whose magnitude is within a few orders of magnitude of the square root of machine roundoff (typically 10 −6 for 64-bit double precision). Here, we briefly describe how to form the IMEX function F(T). We refer F(T) as the IMEX function, since it uses both explicit (hydrodynamics) and implicit (diffusion) information. Notice that for a method that uses all implicit information, F(T) would correspond to a regular nonlinear residual function. The following pseudo code describes how to form F(T) (we also refer to Fig. 1).

Evaluating F(T k ) :
Given T k where k represents the current Newton iteration. Call Hydrodynamics block with (ρ n , u n , E n , T k )tocomputeρ n+1 , u n+1 , E * . Form F(T k ) based on the Crank-Nicolson method, It is important to note that we are not iterating between the implicit and explicit blocks. Instead we are executing the explicit block inside of a nonlinear function evaluation defined by F(T k ). The unique properties of JFNK allow us to perform a Newton iteration on this IMEX function, and thus JFNK is a required component of this nonlinearly converged IMEX approach.

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 where the unknown υ f represents the front velocity. This gives As noted in Rider & Knoll (1999), to avoid problems from lack of smoothness the following numerical approximation is used to calculate υ f Then the new time step is determined by the Courant-Friedrichs-Lewy (CFL) condition

An IMEX Method for the Euler Equations That Posses Strong Non-Linear Heat Conduction and Stiff Source Terms (Radiation Hydrodynamics)
www.intechopen.com where ∆r uses the L 1 norm as in Equation (43). We can further simplify Equation (44) by using Equation (43), i.e, 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 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, where where the parameter E 0 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.

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 where c 0 is a constant and set to 1/4 for this test. Note that c 0 → 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 r i−1/2 to r i+1/2 so that where the symbol er f denotes the error function. The initial density is set to ρ = 1/r.T h e initial temperature is calculated by using E = c v ρT + 1 2 ρu where u = 0 initially. The boundary  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=0 = 0). The coefficient of thermal conduction is set to κ(T)=T 5/2 . We run the code until t = 0.01 with ε 0 = 100 using 200 cell points. The size of the computational domain is set to 1 (e.g, R 0 = 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 ∆t , E ∆t/2 , E ∆t/4 , ···). Then we measure the L 2 norm of errors between two consecutive time step solutions ( E ∆t − E ∆t/2 2 , E ∆t/2 − E ∆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.

An IMEX Method for the Euler Equations That Posses Strong Non-Linear Heat Conduction and Stiff Source Terms (Radiation Hydrodynamics)
www.intechopen.com

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, 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  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 where ε 0 = 235 and c 0 = 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 = c v ρT/γ + 1 2 ρu 2 with the initial u = 0. The radiation diffusivity (κ in Equation (19)) is calculated by considering the LERH model and comparing it with the sum of Equation (18) is compared to Equation (6) of . Then κ becomes where κ 0 = 10 2 , a = −2andb = 13/2 as in . We set P = 10 −4 and σ a = 10 8 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.

Radiative shock test
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(ρ 2 u 2 − ρ 1 u 1 )=(ρ 2 u 2 2 + p 2 + P p ν,2 ) − (ρ 1 u 2 1 + p 1 + P p ν,1 ), 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 where c 1 = γ p 1 ρ 1 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; ρ 1 = 1.0, T 1 = 1.0, and M 1 = u 1 /c 1 = 1.2 as the upstream Mach number. Then we calculate the pressure from a calorically perfect gas relation (p 1 = Rρ 1 T 1 ). Using p 1 and ρ 1 , we calculate the upstream sound speed c 1 = γp 1 /ρ 1 together with u 1 = M 1 c 1 .W i t h these information in hand, we can easily calculate the downstream values using Equations  by using the energy relation E = c v ρT + 1 2 ρu 2 . The radiation temperature is assumed to be a smooth function across the shock and equal to T 1 and T 2 on the left and right boundary of the computational domain, i.e., we choose The initial radiation energy is calculated by E ν = T 4 ν . Other parameters that appear in Equations (16)-(19) are set as P = 10 −4 , σ a = 10 6 ,a n dκ = 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.

Convergence analysis
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.

An IMEX Method for the Euler Equations That Posses Strong Non-Linear Heat Conduction and Stiff Source Terms (Radiation Hydrodynamics)
www.intechopen.com

A classic IMEX method
The classic IMEX method operates on Equations (20)-(21) as follows. Explicit block: Step-1: Step-2: Implicit block: 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 Substituting Equation (66) into (67), we get We let L(E n )=− ∂ ∂x {u n [E n + p n ]} and use T n = T 1 − ∆tT n t + O(∆t 2 ) , then (68) becomes

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.