*2.4. The Implicit-Explicit Time Discretization*

As the penalization introduces stiff terms to the equations and for accuracy, we would want them to be as stiff as possible, we introduce an implicit time integration for those terms. With an otherwise explicit time integration scheme, this results in an implicit-explicit time-stepping scheme that can be achieved by splitting the right-hand side of the equations into an explicitly integrated part and an implicitly integrated part. Therefore, to perform time integration of the system, we use a Diagonally Implicit Runge–Kutta (DIRK) scheme with three explicit and four implicit stages as presented in [14]. The following section first considers a single implicit Euler step, to discuss the arising equations that need to be solved in each implicit stage of the higher order time discretization.

We denote the right-hand side by *Q* and employ the superscript *ι* for the implicit part and the superscript *ξ* for the explicit part. By using the conservative quantities as subscripts (*ρ*, *m*, and *e*), we can distinguish the right-hand sides for the different equations. Thus, we get:

$$\frac{\partial \rho}{\partial t} = Q\_{\rho}^{\mathbb{E}} + Q\_{\rho}^{\iota} \tag{18}$$

$$\frac{\partial m\_i}{\partial t} = Q\_{m\_i}^{\mathbb{F}} + Q\_{m\_i}^{\iota} \tag{19}$$

$$\frac{\partial \mathcal{e}}{\partial t} = Q\_{\mathcal{e}}^{\mathbb{F}} + Q\_{\mathcal{e}}^{\iota} \tag{20}$$

and we chose the implicit parts as:

$$Q^t\_\rho = 0\tag{21}$$

$$Q'\_{m\_i} = -\frac{\chi}{\eta} \left( u\_i - \mathcal{U}\_{oi} \right) \tag{22}$$

$$Q\_\varepsilon^t = -\frac{\chi}{\eta\_T} \left( T - T\_o \right) \tag{23}$$

out of Equations (7) and (9).

This choice restricts the implicit computation to the local source terms, which can be evaluated pointwise. Unfortunately, the introduced Brinkman porosity in (7) affects the flux and introduces spatial dependencies. To avoid the need for the solution of an equation system across the whole domain for this dependency, the porosity part will be computed in the explicit time-stepping scheme.

### Observation for the Implicit Part

Considering Equations (18)–(20) only with their implicit parts, we get the following equation system:

$$\frac{\partial \rho}{\partial t} = 0\tag{24}$$

$$\frac{\partial m\_i}{\partial t} = -\frac{\chi}{\eta} \left(\mu\_i - \mathcal{U}\_{oi}\right) \tag{25}$$

$$\frac{\partial \varepsilon}{\partial t} = -\frac{\chi}{\eta\_T} \left( T - T\_\delta \right). \tag{26}$$

Notice that these equations can be solved pointwise as no spatial derivatives appear.

A discretization of these equations in time with a Euler backward scheme yields the solvable equation system:

$$\frac{\rho(t+\Delta t) - \rho(t)}{\Delta t} = 0\tag{27}$$

$$\frac{m\_i(t + \Delta t) - m\_i(t)}{\Delta t} = -\frac{\chi}{\eta} \left( u\_i(t + \Delta t) - \mathcal{U}\_{oi} \right) \tag{28}$$

$$\frac{e(t+\Delta t) - e(t)}{\Delta t} = -\frac{\chi}{\eta\_T} \left( T(t+\Delta t) - T\_o \right). \tag{29}$$

Equation (27) trivially yields *ρ*(*t* + Δ*t*) = *ρ*(*t*). With the implied constant density, we can now solve the equation for the change in momentum (28) and arrive at an explicit expression for the velocity *ui*(*t* + Δ*t*) at the next point in time:

$$\frac{\rho(t+\Delta t)u\_i(t+\Delta t) - \rho(t)u\_i(t)}{\Delta t} = -\frac{\chi}{\eta} \left( u\_i(t+\Delta t) - lL\_{oi} \right) \tag{30}$$

$$u\_i(t+\Delta t) = \frac{\rho(t)u\_i(t) + \frac{\chi\Delta t}{\eta}\mathcal{U}\_{oi}}{\rho(t) + \frac{\chi\Delta t}{\eta}}.\tag{31}$$

Finally, density and velocity at the new point in time can be used to find the new temperature as well by substituting the above results in Equation (29) and solving for the temperature at the next point in time. We find:

$$T(t + \Delta t) = \frac{\frac{\chi \Delta t}{\eta\_T} T\_o + c\_v \rho(t) T(t) + \frac{\rho(t)}{2} (u\_i^2(t) - u\_i^2(t + \Delta t))}{c\_P \rho(t) + \frac{\chi \Delta t}{\eta\_T}}. \tag{32}$$

where *ui*(*t* + Δ*t*) is given by Equation (31).

Thus, this specific choice of terms for the implicit part of the time integration scheme yields a system that can be solved explicitly and without much additional computational effort. However, the implicit discretization allows for arbitrarily small values of *η* and *ηT*. A similar approach was developed by Jens Zudrop to model perfectly electrical conducting boundaries in the Maxwell equations, and more details can also be found in his thesis [15].

To solve the complete system, we then employed the diagonally implicit Runge–Kutta scheme with three explicit stages and four implicit stages [14]. It provides a scheme that is third order in time and L-stable.

Note, that while this approach overcomes time step limitations with respect to the permeabilities *η* and *ηT*, the porosity term changes the eigenvalues of the hyperbolic system and affects the stability.
