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

Samet Y. Kadioglu<sup>1</sup> and Dana A. Knoll2

<sup>1</sup>*Idaho National Laboratory, Fuels Modeling and Simulation Department, Idaho Falls* <sup>2</sup>*Los Alamos National Laboratory, Theoretical Division, Los Alamos USA*

#### **1. Introduction**

292 Hydrodynamics – Advanced Topics

Papageorgiou, L., Metaxas, A. C. & Georghiou, G. E. (2011), Three-dimensional numerical

Sarrette, J. P. ; Cousty, S. ; Merbahi, N. ; Nègre-Salvayre, A. & F. Clément (2010),

Sigmond, R. S. (1984), The residual streamer channel: Return strockes and secondary

Spyrou, N.; Held, B.; Peyrous, R.; Manassis, Ch. & Pignolet, P. (1992) Gas temperature in a

van Veldhuizen, E. M. & Rutgers, W. R. (2002), Pulsed positive corona streamer propagation and branching, *Journal of Physics D: Applied Physics*, Vol. 35, pp. 2169–2179 Urashima K. & J. S. Chang (2010) Removal of Volatile Organic Compounds from Air

*Transactions on Dielectrics and Electrical Insulation,* Vol. 7 No. 5, pp. 602-614 Villeger, S.; Sarrette, J. P.; Rouffet, B. ; Cousty S. & Ricard A. (2008), Treatment of flat and

Winands, G.; Liu, Z.; Pemen, A.; van Heesch, E.; Yan, K. & van Veldhuizen, E. (2006),

Yousfi, M. & Benabdessadok, M. D. (1996), Boltzmann equation analysis of electron-

Yousfi, M.; Hennad, A. & Eichwald, O. (1998), Improved Monte Carlo method for ion

streamers, *Journal of Applied Physics*, Vol. 56, No. 5, pp. 1355-1370

D: Applied Physics, Vol. 25, pp. 211-216

*Physics D: Applied Physics*, Vol. 39, pp. 3010–3017

*Applied Physics*, Vol. 84, No. 1, pp. 107-104

*Physics*, Vol. 42, pp. 25-32.

*Physics*, Vol. 80, pp. 6619-6631

phenomena, *Journal of Physics D: Applied Physics*, Vol. 44, 045203 (10pp). Penetrante, B. M. & Schultheis, S. E. (1993), Nonthermal Plasma Techniques for Pollution

Karlsruhe.

648–660

modelling of gas discharges at atmospheric pressure incorporating photoionization

Control, Part A&B, Editors, NATO ASI Series Vol. G 34, Springer-Verlag,

Observation of antibacterial effects obtained at atmospheric and reduced pressures in afterglow conditions, *European Physical Journal. Applied. Physics,* Vol. 49 13108 Segur, P.; Bourdon, A.; Marode, E., Bessieres D. & Paillol J. H. (2006) The use of an

improved Eddington; approximation to facilitate the calculation of photoionization in streamer discharges, Plasma Sources Sciences and Technologies, Vol. 15, pp.

secondary streamer discharge: an approach to the electric wind, Journal of Physics

Streams and Industrial Flue Gases by Non-Thermal Plasma Technolog, *IEEE* 

hollow substrates by a pure nitrogen flowing post discharge: Application to bacterial decontamination in low diameter tubes, *European Physical Journal. Applied.* 

Temporal development and chemical efficiency of positive streamers in a large scale wire-plate reactor as a function of voltage waveform parameters, *Journal of* 

molecule collision cross sections in water vapor and ammonia, *Journal of Applied* 

transport in ion-molecule asymmetric systems at high electric fields, *Journal of* 

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 & 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; 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

$$
\mu\_l = f(\mu) + \mathcal{g}(\mu), \tag{1}
$$

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:

$$\frac{\mu^\*-\mu^n}{\Delta t} = f(\mu^n),\tag{2}$$

**Remark:** We remark that another way of achieving a self-consistent IMEX integration that preserves the formal numerical accuracy of the whole system is to improve the lack of influence of the explicit and implicit blocks on one another by introducing an external iteration procedure wrapped around the both blocks. More details regarding this methodology can be

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

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

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.,

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

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

in a self-consistent IMEX manner can be found in (Kadioglu & Knoll, 2011).

found in (Kadioglu et al., 2005).

**2. Applications**

2009).

Implicit block solves:

$$\frac{\mu^{n+1} - \mu^\*}{\Delta t} = g(\mu^{n+1}).\tag{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:

$$\begin{aligned} \mathfrak{u}^1 &= \mathfrak{u}^n + \Delta t f(\mathfrak{u}^n) \\ \mathfrak{u}^\* &= (\mathfrak{u}^1 + \mathfrak{u}^n)/2 + \Delta t/2 f(\mathfrak{u}^1) \end{aligned} \tag{4}$$

Implicit block:

$$
u^{n+1} = 
u^\* + 
\Delta t/2[\mathbf{g}(
u^n) + \mathbf{g}(
u^{n+1})].\tag{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; Kadioglu, Knoll & Lowrie, 2010).

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 terms in Equation (6).

Explicit block:

$$u^1 = u^n + \Delta t f(u^n)$$

$$u^\* = (u^1 + u^n)/2 + \Delta t/2 f(u^{n+1})\tag{6}$$

Implicit block:

$$
\mu^{n+1} = \mu^\* + \Delta t / 2[g(\mu^n) + g(\mu^{n+1})].\tag{7}
$$

**Remark:** We remark that another way of achieving a self-consistent IMEX integration that preserves the formal numerical accuracy of the whole system is to improve the lack of influence of the explicit and implicit blocks on one another by introducing an external iteration procedure wrapped around the both blocks. More details regarding this methodology can be found in (Kadioglu et al., 2005).
