*2.2. The Brinkman Penalization*

Penalization schemes employ additional, artificial terms to the equations in regions where the flow is to be inhibited (penalized). In the conservation of momentum and energy, we can make use of local source terms that penalize deviations from the desired state. With the Brinkman penalization, we also inhibit mass flow through obstacles by introducing the Brinkman porosity model and using a low porosity, where obstacles are to be found. Extending the compressible Navier–Stokes equations from Section 2.1 by the penalization terms, we obtain Equations (7) and (9).

$$\frac{\partial \rho}{\partial t} = -\nabla \cdot \left[1 + \left(\frac{1}{\rho} - 1\right) \chi\right] \mathbf{m}, \qquad \text{(7)}$$

$$\frac{\partial m\_i}{\partial t} + \sum\_{j=1}^3 \frac{\partial}{\partial x\_j} \left(m\_i v\_i + p \delta\_{ij}\right) - \frac{1}{Re} \sum\_{j=1}^2 \frac{\partial}{\partial x\_j} \pi\_{ij}$$

$$= -\frac{\chi}{\eta} \left(v\_i - \mathcal{U}\_{oi}\right) \qquad i = 1, 2, 3, \tag{8}$$
 $\frac{\partial \rho \epsilon}{\partial t} + \nabla \cdot \left[\left(\varepsilon + \frac{p}{\rho}\right) \mathbf{m}\right] - \frac{1}{Re} \sum\_{j=1}^2 \frac{\partial}{\partial x\_j} \left(\sum\_{i=1}^2 \tau\_{ij} v\_i - \frac{1}{\gamma - 1} \frac{\mu}{Pr} T\right)$ 
$$= -\frac{\chi}{\eta \tau} \left(T - T\_o\right). \tag{9}$$

The obstacle has the porosity *φ*, the velocity *Uo*, and the temperature *To*. The strength of the source terms can be adjusted by the viscous permeability *η* and the thermal permeability *ηT*. The masking function *χ* describes the geometry of obstacles and is zero outside of obstacles and one inside. It is also referred to as the characteristic function. It is capable of dealing not only with complex geometries but also with variations in time.

$$\chi(\mathbf{x},t) = \begin{cases} 1, & \text{if } \mathbf{x} \in \text{obstacle}. \\ 0, & \text{otherwise}. \end{cases} \tag{10}$$

To represent a solid wall for compressible fluids properly, Liu et al. [9] stated that the porosity *φ* should be as small as possible, i.e., 0 < *φ* << 1. They scaled the permeabilities with the porosity and introduced according scaling factors *α* and *αT*. The permeabilities were then defined by *η* = *αφ* and *η<sup>T</sup>* = *αTφ*. With these relations, Liu et al. [9] found a modeling error of *O*(*η*1/2*φ*) for resolved boundary layers in the material and *O*((*η*/*ηT*)1/4*φ*3/4) for non-resolved boundary layers. In both cases, the error was dominated by the porosity. Nevertheless, the error can still be minimized with sufficiently small viscous permeabilities *η*.

Moreover, small values of the porosity caused stability issues and imposed a heavy time-step restriction with our numerical scheme. With the introduction of *φ*, the eigenvalues of the hyperbolic system changed, which has adverse effects on stability. The eigenvalues of the system of equations along with penalization terms [9] are given by the following characteristic equation:

$$-\left(\lambda - u\right)^{3} + \left[c^{2} + \frac{u^{2}}{2}(\phi^{-1} - 1)(\gamma - 3)\right](\lambda - u) - c^{2}u(\phi^{-1} - 1)(\gamma - 1) = 0,\tag{11}$$

where *c* = (*γp*/*ρ*)1/2 and *γ*, *p*, *ρ*, and *u* are the ratio of specific heat, pressure, density, and velocity, respectively. For *φ* = 1, the system of equations yields three eigenvalues *u*, *u* + *c*, *u* − *c*, which implies the speed of sound *c* in the medium, which is what we would like to achieve. However, with 0 < *φ* << 1, the eigenvalues can no longer be evaluated easily and are linked to *φ*, which causes problems for the hyperbolic part.
