Pressure-Tight and Non-stiff Volume Penalization for Compressible Flows

Embedding geometries in structured grids allows a simple treatment of complex objects in fluid simulations. Various methods for embedding geometries are available. The commonly used Brinkman-volume-penalization models geometries as porous media, and approximates a solid object in the limit of vanishing porosity. In its simplest form, the momentum equations are augmented by a term penalizing the fluid velocity, yielding good results in many applications. However, it induces numerical stiffness, especially if high-pressure gradients need to be balanced. Here, we focus on the effect of the reduced effective volume (commonly called porosity) of the porous medium. An approach is derived, which allows reducing the flux through objects to practically zero with little increase of numerical stiffness. Also, non-slip boundary conditions and adiabatic boundary conditions are easily constructed. The porosity terms allow keeping the skew symmetry of the underlying numerical scheme, by which the numerical stability is improved. Furthermore, very good conservation of mass and energy in the non-penalized domain can be achieved, for which the boundary smoothing introduces a small ambiguity in its definition. The scheme is tested for acoustic scenarios, for near incompressible and strongly compressible flows.

The effective volume for the fluid can be interpreted in a quasi-one-dimensional flow as a reduced cross-section by reorganizing the material of a flow tube, yielding a configuration known from the theory of stream lines. This gives an extra pressure source term in the momentum equation a simulation with a refined grid is presented. The results are summarized in Sect. 5. The Appendix A provides technical details of the finite difference method and the filter.

The Penalized Equations
Here, the modification of flow equations to represent immersed objects will be derived in analogy to one-dimensional gas dynamics.

Derivation of the Analytical Equations
The modified equations are to include the effects of reduced effective volume fraction φ = V fluid /V total , since the volume of the material occupied by porous media is unavailable for the fluid. It is expected that the equations mainly scale with the volume fraction since mass, momentum, and energy scale likewise. The Darcy term, describing linear friction proportional to the relative velocity between flow and object, is added in a second step. As stated by Kevlahan et al. [34] for the shallow water equations, it is important that mass and momentum (and for the Navier-Stokes equations, energy) are transported with the same speed. We also want to keep the speed of sound c unaltered by the penalization, since this suggests that the same CFL condition CFL = (|u| + c)Δt/Δx stays valid [34]. The validity of this assumption is discussed below in Sect. 4.1.3. Instead of deriving the equations from first principles by the theory of variations as in [34], we observe that for one dimension the Euler equations of varying cross-section are structurally similar to the equations derived in [34]. The speed of sound is independent of the cross-section φ, and the transport of mass and momentum are modified in the same way. Further, the essential modification by φ in [34] is a non-conservative pressure gradient φ∂ x p, as it appears in one-dimensional gasdynamics, see Eq. (5). This can be motivated by the fact that no specific distribution of the porous material was prescribed. By this, it can be reorganized microscopically as a tube of reduced cross-section, see Fig. 1. This was recognized before, see [3] and references therein.
This leads to an extra factor φ in every flux term and an additional source term p∂ x φ for the momentum equation, which is detailed in text-books of gasdynamics, e.g. [1]. A similar modification of the pressure gradient can be identified in the shallow water equation for an uneven bottom height [34] and the Navier-Stokes equation in cylindrical coordinates, where the role of the reduced cross-section is the circumference proportional to the radius r |. 1 It is also found for porous material [44] and in the Baer-Nunziato model for multi-phase flow [3], see below. In all these applications a reduction of a cross-section-like property yields a source term in the momentum equation of this type.
The equations for mass, momentum, and energy in one dimension without friction are therefore just the equations for one-dimensional gas dynamics in varying cross-section φ We introduced the mass density ρ, the velocity u, the pressure p and the total energy e t = u 2 /2 + e in , where we in the following assume the internal energy of an ideal gas e in = p ρ 1 γ −1 , with the adiabatic index γ . In principle, the varying cross-section φ can have arbitrary (positive) values. In the following it is identified with the reduced volume by porosity, so that it is in the range φ = [0, 1]. In the vanishing limit φ → 0 no volume is left for the fluid. The volume fraction appears linearly in all terms on the left-hand side, whereas the term on the right-hand side is a source term without a derivative of any of the dynamical variables. From this, it follows that the equations are hyperbolic with unchanged characteristic velocities: the acoustic waves with λ ± = u ± c with the speed of sound c = √ γ p/ρ and the entropy wave λ s = u. Further, all flux terms are modified in the same manner, so that the consistency demand of [34] is automatically fulfilled. If φ is constant (and non-zero), the equations reduce to the standard one-dimensional Euler equations.
Due to the source term p ∂ x φ, the momentum is not conserved, whereas the total mass and energy are conserved. This is physically sound, since a change in cross-section creates reflections, whereas mass and energy are unaltered.
The Eqs. (1)-(3) can be generalized straightforwardly to multiple dimensions and augmented with dissipative terms. Summing convention is assumed for all Greek indices α, β = 1, 2, 3 marking spatial directions. The two pressure terms in the momentum equation have been combined ∂ x α (φ p) − p∂ x α φ = φ∂ x α p, resulting in the equations Further, the Darcy friction term is included in the momentum equation with a spatially dependent force strength χ = 1/η and the target value u t α , which is the speed of the immersed object, which is always zero in this report. We also assume that the objects are static, i.e. φ and χ are time independent, in the following. The dissipative term contain the temperature T = p/ρ RW with the universal gas constant R and the molecular weight W . The viscous friction is given by τ αβ = μ ∂ x β u α + ∂ x α u β + (μ d − 2/3μ)δ α,β ∂ x γ u γ with μ and μ d the shear and volumetric friction and δ α,β the Kronecker delta. Similar to the Darcy term, a penalization can be added to the energy equation to enforce e.g. isothermal boundary conditions, [11]. These dissipative fluxes are scaled by φ since these fluxes take place in the fluid part only. It also keeps the symmetry of the dissipation terms, i.e. guarantees negative semi-definite operators. Comparison with former schemes Liu and Vasilyev [40] discuss different equations for the Brinkman penalization, of which some agree with here proposed equations, while the equations used in the numerical implementation [40, eqn. (18-20)] do not. For example their mass equation is with χ the mask function and m α = ρu α . In contrast, the mass equation in the initial discussion [40, (eqn 2)] and the momentum equation [40, (eqn 5)], attributed to Wooding [60], agrees with our equations. 2 For the latter, the mass equation (4) needs to be split off the momentum equation (5) ρ∂ by which the volume fraction cancels in all terms beside the shear friction term. This is structurally also found by Kevlahan et al. [34]. However, the (steady state) momentum equation given by Wooding [60, (eqn. 6)] in the cited reference is with q = φu (aligned to our notation), with g the gravity vector and k the Darcy friction factor. The combination of the pressure term and the non-linear transport term disagrees with (8) and not all terms have the same order in φ. Liu and Vasilyev abandon Wooding's equation, similar to the here proposed, due to an argument of Beck [7], who, however, uses a Darcy equation (his Eq. 1), where a constant φ does not cancel. Beck finds an insufficient number of boundary conditions for his buoyancy flow, however, no problems are apparent for the here used Eqs. (4)(5). The momentum equation of Nield and Bejan [44, (page 8)] agrees with the here found equation. They argue to drop the non-linear transport, by citing Beck, and by practical considerations for a real porous material, which does not seem imperative for our artificial application of porous model. Instead, it is shown in the simulations that equations used here, with all terms linear in φ, allow to choose the volume fraction φ and the Darcy term largely independently, making it possible to model different boundary conditions. The way φ enters the equations in this article agrees with the equations of Kemm et al. [33], derived from the two-phase flow of Baer-Nunziato model. A similar approach for acoustic waves is given by Tavelli et al. [53]. Moving Objects Moving objects imply a time dependency of the cross-section φ = φ(t), which results in additional terms in the governing equations. These are derived here for the sake of completeness and fully agree with Kemm et al. [33].
A time-dependent reduced volume can be interpreted as a flux term, see Fig. 2. For example for the mass equation, the balance equation of an infinitesimal volume is Here, we assume that the horizontal fluid velocity is constant for one x-position, so that the flux along the boundaries at x 1 and x 2 is a factor φ. The vertical velocity at the lower boundary is assumed to be zero v(0) = 0 and v(φ) = ∂ t φ at the upper boundary. For an infinitesimal Δx the volume integral can be replaced by a multiplication with the volume Ω dV = φΔx. By this (10) results in φΔx∂ t ρ + (φuρ) x 2 − (φuρ) x 1 + Δxρ∂ t φ = 0, or for Δx → 0 and with the help of the product rule  -dependent  effective volume introduces  addition flux terms, asφ can be  identified with a velocity at the  boundary, see text In the same manner, the transport terms for momentum and energy produce the term necessary to include φ in the time derivative. Reconsidering all terms of (4)(5)(6) in this manner, one finds that only one further term is created by a time dependent φ which arises from the pressure work term in the energy equation. The same reasoning as above shows that the pressure work The Eqs. (4-6) become A similar derivation for moving walls for the flow of blood is presented in [5], where the density becomes a constant due to the incompressibility assumption. The equations agree with Kemm et al. [33] if the Darcy and viscous terms are omitted.
To elucidate the relevance of these extra terms, we use the product rule for time derivative terms, and splitting of the kinetic energy u α u α /2 = e t − e in from the total energy equation to arrive at (in 1D) For the total energy of an ideal gas e t = ( p/ρ) · 1/(γ − 1) the last equation becomes Without the spatial part a change in φ produces a change in the density given by dφ/φ + dρ/ρ = 0 and φdp + γ pdφ = 0 from which the adiabatic relation dp/ p − γ dρ/ρ = 0 follows. This shows that a temporal changing φ is just an adiabatic pressure source, as expected for a homogeneous compression. Splitting off the mass equation from momentum shows that the change in momentum is fully due to the change of density. A detailed discussion of the dynamic φ is found in [33]. Implementation Strategy To describe embedded objects, the friction term χ can be chosen to be large or φ can be made small, where of course a combination is possible. The strategy will be to reduce the volume fraction by several orders of magnitude (typically down to 10 −6 − 10 −8 ) so that the remaining volume is practically zero and by this all convective and viscous fluxes. Second, the Darcy friction is chosen as big as the time step permits for an explicit time marching. Stronger values the Darcy friction could be handled by an implicit or semi-implicit method as, for example, used in [11]. However, a very large choice corresponds to simply setting the target values inside the objects. This reduces the smoothness of the solution having a detrimental effect on the discrete solution. This effect is discussed in more detail for example in [22,57]. The reduced volume fraction implies a good approximation of the non-penetration condition, producing thereby a slip boundary condition, while the Darcy friction ensures the non-slip part. Close to a boundary, large shear forces result from this, restricting the time step. However, since the Darcy term needs to balance forces similar to the shear forces in the fluid itself, this leads to a very similar stiffness and no additional time step restriction results.

Conservation Properties
The conservation of mass and energy in the full domain is obvious from (4) and (6), while the momentum (5) is not conserved. However, in the case of immersed objects, it is not the conservation in the full domain Ω = Ω f ∪ Ω s that is physically relevant, but the conservation in the fluid Ω f domain without the objects Ω s . If we first consider a penalization by a sharply varying φ such that φ = Ω s , inside the objects 1 Ω f , otherwise (19) and remark that the conserved quantities are weighted by φ, we find e.g. for the mass By this, the total mass is decomposed in mass in the fluid and in the solid region M = M f + M s . The last term is suppressed by a small volume fraction of . If the value of the density in the porous region is bounded, its fraction of the total mass is suppressed by → 0, the mass in the fluid domain is approximated by the total mass, and therefore approximately conserved. In many practical cases, the maximal values of ρ in the objects can be estimated from physical reasoning. This can be stagnation values in external flows or a doubling of pressure and density for acoustic reflections. We find very good conservation in tests, with the only ambiguity being the finite thickness of the transition of φ from one to , which is needed for numerical stability. This smearing of the boundary induces some arbitrariness of the exact location of the object and of the separation of the two integrals in (20). However, mass leaving the fluid domain is not lost, but 'inside' the boundary region so that this property can be viewed as a not perfectly rigid boundary. The same reasoning can be made for the energy. The momentum is not conserved, as mentioned before, as objects act as momentum sources by deflection flow or reflecting acoustics.
By this, the volume fraction approach implies an (approximate) adiabatic boundary condition, if no temperature is forced. Since the energy content of the gas in the solid region is suppressed by φ, its total energy content and the fluxes are suppressed by the same factor. By this, only negligible energy enters or leaves the object, assuming again finite values of the energy density in the object.

Discrete Equations
Immersed boundaries allow using simple Cartesian grids. This permits to use finite difference methods, which are easy to implement and to parallelize. Finite differences allow for small discretization errors by simply using appropriate finite difference stencils. Various finite difference formulations are available. The here used method is a skew-symmetric method, due to the combination of strict conservation and the absence of numerical damping. The first is essential for the faithful treatment of shocks, while the latter is advantageous for the simulation of acoustics and fine turbulence. For skew-symmetric schemes, the momentum transport term is formulated to yield a skew-symmetric matrix after discretization. This has the elementary property that quadratic forms become zero q T Aq = 0 if the matrix A is skew-symmetric A = −A T . This implies that the discrete momentum transport term does not change the kinetic energy, in agreement with the analytical theory. A violation of this fact is a major source of instabilities since numerical damping or numerical excitation are equally likely.
It is found that the derivation of the skew-symmetric scheme of Reiss et al. [51] can directly be adopted to include φ. Details can be found in the Appendix A, here we restrict the discussion to a brief description. However, this special discretization is not mandatory and it is expected that other discretizations work as well.
We use √ ρ, √ ρu α and p as calculation variables. This unusual choice origins from a symmetric splitting of the kinetic energy ρu α u α /2 = ( √ ρu α )( √ ρu α )/2. This suggests using √ ρ as a variable to form the velocity u α = ( √ ρu α )/ √ ρ, without evaluating the square root of a field. The final equations (for a time constant volume fraction φ) are The equations are discretized simply by replacing all derivatives with finite differences with symmetric stencils and evaluating all products or divisions pointwise. The time-stepping is the standard fourth-order Runge-Kutta scheme. This is followed by filtering for cases where otherwise highly oscillatory solutions occur, namely when complex geometries focus waves and most importantly for shocks where dissipation is demanded by physics. A specially designed filter is needed in general, which is discussed in Sect. 3. The whole method is fully explicit and simple to implement. A source code producing all the numerical examples is provided as supplementary information on the journal homepage. 3 It seems attractive to include φ in redefining ρ and p, but big rounding errors were found in this case for a small φ, so we refrain from doing so. Note that we divide by φ for timestepping, so that it cannot be set to zero in the numerical implementation. The transport term in the momentum equation is explicitly skew-symmetric, which is preserved if the spatial discretization simply replaces the derivatives with central finite difference matrices. The kinetic energy is changed only by compression work, and not by numerical friction. The total energy as the sum of the kinetic and internal energy is conserved.
The skew-symmetric form suggests an unconditionally stable scheme. For this, a timestepping that conserves quadratic invariants is needed as presented by Brouwer et al. [12]. To be rigorously norm-stable, the variable √ p instead of the pressure p has to be used, however, this form works very well in our experience. All methods which conserve quadratic invariants are implicit methods. This extra effort does not seem to be justified for many time dynamic simulations, where all scales of interest need to be resolved since stability does not guarantee the correctness of the fast dynamic. Explicit high order time-stepping yields very good conservation properties, so that it can be used in most practical cases as shown for different examples in [12]. The conservation can cheaply be validated during simulation permitting to switch to a different time integrator. A reduction of the time step might often be sufficient as thereby the conversation is strongly improved [12]. Explicit time integrators with improved energy convergence are discussed for incompressible flows by Capuano et al. [14], which might be adopted for our method.

Filtering
The non-linearity of the Navier-Stokes equations constantly produces higher wavenumbers, which tend to exceed the grid resolution after some time if not removed by friction terms. This is especially important in shocks, but also at boundaries. The high gradient of the immersed objects can create high-frequency structures up to a grid-to-grid oscillation. To remove these oscillations, artificial dissipation is often applied. Various methods to introduce dissipation exist, an overview is given by Pirozzoli [48]. In the simplest case it is introduced by upwind discretizations or similar by calculating fluxes by a Riemann solver, which creates in effect an upwind stencil for the Riemann invariants [39].
Often hybrid methods are used, which switch based on a flux limiter from a high order, low dissipative to low order dissipative scheme. Alternatively, artificial viscosity by an adaption of the viscous parameters near discontinuities was proposed by Cook et al. [15]. Artificial viscosity and the flux limiting are at least for simple examples structurally similar, see [39,Chap. 16]. Adaptive filtering was introduced with the same purpose by Yee et al. [61], where a filter operation with local filter strength allows to concentrate dissipation near discontinuities. With a similar effect, varying filter order can be applied as used by Visbal et al. [58] and more recently by Patel et al. [45]. A quantitative comparison of the different methods is provided by Johnsen et al. [31]. These non-linear dissipation methods rely usually on a detector of shocks or non-smoothness, which is used to determine the flux limiter, the artificial viscosity, or the filter strength. Different shock detectors are compared in [48].
In the following, a filter is used in the numerical examples. This should not only remove oscillatory components but also allow a locally varying filter strength and especially respect the conservation with the varying volume fraction. Conservative filtering is often achieved by filtering the conserved quantities, however, this is not appropriate here. Consider the mass in one dimension, M = i (ρφ) i Δx, where i runs over all grid points. Filtering (ρφ) i would keep the total mass constant, however, it would tend to increase ρ where φ is becoming small, i.e. inside the boundary of immersed objects. This counteracts the intention to smooth the physical quantities.
Conservative filters similar to [36] are now formulated analogously to dissipative flux terms, which fulfill the desired properties, as argued in the following. The starting point is to define fluxes between two adjacent grid points, similar to finite volume for some variable q (depicted in 1D) The filter strength σ i can be locally varying and is thereby chosen to filter oscillations within objects or in under-resolved parts of the flow, e.g. close to shocks.
To generalize the dissipation-like filtering (24), it can be written in matrix notation as where the expression (24) For the derivation periodic grids are assumed from which the non-periodic case is constructed at the end. This method is second-order if σ and φ are constant, and first order otherwise. This form reveals the principle structure whereby the damping and the conservation properties can be checked easily, as now detailed. WithD = −D T , and a positive semi-definite matrix (M(φσ )) (in our case a diagonal matrix with non-negative elements), we find that the operatorD(M(φσ ))D is negative semi-definite, since Obviously (qφ) = (qφ) holds, if q is constant, since the derivative vanishes for any constant function D1 = 0, with the vector of constant entries (1) i = 1, used to depict the sum. The conservation is checked by due to the telescoping property 1 TD = 0 for periodic grids.
To generalize to higher order, appropriateD (thereby D = −D T ) and M, fulfilling the mentioned properties, can be chosen. We demand that the resulting filter yields the filter stencil of a higher-order filter when φ and σ are constant. In [10]D was kept proportional to a one-sided first derivative to keep the flux form, and D was calculated to create the desired filter stencil. Here, we instead chose to keepD = −D T so that the negative semi-definiteness is guaranteed. We remark that the stencil (1, −2, 1)/4 for D andD = −D T yields the filter (−1, 4, −6, 4, −1)/16 for constant φσ , which is the standard fourth order filter, reducing to second-order for non-constant φ or σ . It has the nice property that a purely linear function is reproduced by filtering, even for non-constant σ or φ. Similarly, using this filter stencil for D we find the filter stencil of eighth order (−1, 8, −28, 56, −70, 56, −28, 8, −1)/256 for constant σ and φ. In both cases, no averaging by M is done, i.e. M = 1. The averaging in the first/second-order method is needed to avoid a directional diffusion if the stencil D is nonsymmetric. In more dimensions, the contributions of the different spatial directions are simply added where different σ α for all spatial directions α are possible. An algorithmic description and the later used adaptive filter strength for shocks are discussed in the Appendix C.
The non-periodic case can be handled with the same D and M as above by choosing σ = 0 close to the boundary for a sufficient number of grid points, which results in a non-periodicity stencil of the resulting filter. For the first/second-order filter, the choice of σ = 0 for the two points next to the boundary in the direction perpendicular 4 to the boundary does not change the boundary values. Setting only the outermost points of σ zero would result in modified outermost points, but without a modification of the flux by the filter [30]. This can be chosen depending on the desired boundary condition.

Acoustic Reflection in 1D
A solid boundary reflects acoustic pulses. Here, we investigate how such a reflection is created by the volume fraction alone. The one-dimensional domain with a length of L = 2m is described with N = 1024 points. The domain is large enough to ignore the domain boundaries for this example. The wall starts at x 0 = 1 m + δ 0 where only the volume fraction is reduced from 1 to = 10 −8 with a hyperbolic tangent with a thickness of δ = Δx as the grid spacing. The shift δ 0 = 0.3Δx is chosen to align the maximum of the reflected and the mirror pulse by visual inspection. This implies that the effective location of the wall is not at φ ≈ 1/2 but shifted for the fraction of a grid spacing.
No Darcy term or filter is used in this example to reveal the effect of the volume fraction. Standard fourth-order central difference and the classical Runge-Kutta-4 time stepping is applied, as detailed in the Appendix A.
The initial data is an adiabatic, Gaussian pulse in pressure and density, with relative height of β = 10 −3 , relative to ρ 0 = 1kg/m 3 , p 0 = 10 5 Pa and σ pulse = 8Δx, which is chosen so that no artefacts of the numerical dispersion error are visible by inspection. For σ pulse = 6Δx, the pulse is followed by clearly visible ringing. The results are presented in Fig. 3. The space-time diagram shows that the right-going pulse is reflected at x = 1m. The speed of sound is clearly the same within the object, as well as outside. The snapshot of p at t ≈ 1.32 ms is shown in the second plot. The left-going and the reflected pulse have the same amplitude, as expected, while the pulse in the object is strongly increased and shows some ringing. Its increase is due to the doubling of the pressure at the reflection point. Some additional increase is created by the restriction to very small φ, which is not further analyzed here. To quantify the quality of the reflection, the reflected pulse is compared with a mirror pulse, and a very good agreement is found. The inset shows the difference relative to the initial pulse amplitude ∼ 2 · 10 −3 .
Finally, the spectrum is shown to reveal the frequency dependency of the difference between the reflected pulse and the reference mirror pulse. A window [x m − 5.5σ pulse , x m + 5.5σ pulse ] is created around location the maximum of reflected pulse x m . The pressure values outside this sharp window are set to zero, and a Fourier transform is calculated for the remaining values. By this, the spectrum of the reflected pulse and the reference pulse is calculated. In Fig. 3 the absolute value of these two is compared. The numerically interesting range is very well represented. Since the absolute value agrees very well, the error visible in Fig. 3c is dominantly a phase error.
All results are largely independent of , if sufficiently small. It can be increased by orders of magnitude without changing the result or decreased further without compromising the numerical stability. The only visible effect is that the ringing of the pulse inside the object and the kink at the boundary becomes stronger with a smaller . The observed displacement δ 0 is much smaller in comparison with the in-depth discussion of [28]. However, a substantially different situation was considered there: an incompressible boundary layer enforced by a pure Darcy term.

Containment
One motivating aspect of this investigation is to improve the penalization for configurations with high-pressure gradients. If the Darcy term is the only term to combat a high-pressure difference, it must be chosen accordingly large, leading to a strong numerical stiffness or even instabilities. In the method proposed here, the pressure gradient is counteracted by a change in volume fraction, which is largely able to handle the pressure gradient even without a Darcy term. In Fig. 4, a set-up is shown where two walls at x = 1 and x = 2 keep a high pressure reservoir of 4 p 0 against p 0 = 10 5 .
Without a Darcy term, a chocked flow is found. The flow accelerates to sonic conditions of Ma = u/c ≈ 1, with the sound velocity c = √ γ p/ρ, creating a supersonic expansion terminated by a shock. Chocking conditions are not fully reached in the wall (Ma ≈ 0.75), which is due to the filter acting similar to a viscosity. The second/fourth-order filter described in Sect. 3 is applied in the whole region. However, since the volume fraction φ = 10 −8 scales the mass flux, it results in a mass flux of the order of ≈ 10 −5 . This is largely independent of

Stiffness of the Method
Based on the CFL number criterion, it was argued that leaving the speed of sound unaltered implies avoiding additional stiffness [34]. However, since the pressure amplitude is rapidly increased and the velocity is reduced when a sound pulse reaches the wall location, some additional stiffness is expected. The magnitude of the source term in the momentum equation (2) suggests the possibility of a massive increase. Yet, a very moderate decrease of the admissible CFL number, depending on the width of the boundary function, was found in the simulation.
To investigate the stiffness, the eigenvalues with maximal magnitude were determined for the described acoustic setup of Sect. 4.1.1, see Fig. 5. For a wide wall function δ Δx we reproduce the stiffness of the right-hand side without the volume fraction φ ≡ 1. For small δ, a sharp increase of the stiffness was found. For the setting of the simulation δ = Δx an increase of roughly 70% was found, consistent with the possible CFL numbers in the simulation. For a slightly stronger smeared boundary δ = 1.5Δx, an increase of 25% of the stiffness was found. As expected, the eigenmodes with the largest eigenvalues are always strongly localized at the boundary location. One conclusion from this is that in order to increase the quality of the boundary, instead of decreasing δ, a finer grid resolution can be less restrictive on the time step. Altogether, the increase is moderate and can often be much smaller than a pure Darcy term, especially for high-pressure gradients.

Acoustic Reflection in 2D
The acoustic case is extended to two dimensions. A square domain x × y = (1/4, . . . , 5/4) × (1/4, . . . , 5/4) is embedded in a domain of size 2 m × 2 m, discretized with 1024 × 1024 points. An adiabatic pulse centered at (1 m, 1 m), again with σ pulse = 8Δx and an amplitude of 10 −3 ρ 0 , ρ 0 = 1, is used as initial condition. A reference can again be created by the three mirror pulses, located at (1.5 m, 1 m), (1 m,1.5 m), (1.5 m, 1.5 m) with the same amplitude as pulse inside the box. The pulse is reflected by the jump in φ, modeling the wall. Again, acoustic waves are created within the object, where for a locally constant φ = 10 −8 in effect the same Euler equations are valid. The mass, momentum, and energy content of the wall region are small by the factor of 10 −8 as discussed before. Again, we find excellent agreement, see Fig. 6. Not only is the acoustic wave, which impinges perpendicular to the wall, well reproduced, furthermore the wave slipping along the wall at x = 1.25 is correctly reproduced. Such a slip boundary cannot be presented if a Darcy term is used to describe objects since this necessarily creates non-slip boundaries.

Potential Cylinder Flow
To investigate the quality of the slip boundary condition, the friction-free flow around a cylinder is calculated, for which exists a simple analytical solution. A flow velocity at infinity of u 0 = 10 m/s is used, so that (for reference values as above) a low Mach number of Ma ≈ 0.02 allows to consider the flow as near incompressible. As boundary conditions, the known potential flow u pot = ∂ x Φ pot and v pot = ∂ y Φ pot , with the potential in polar coordinates relative to the cylinder center Φ pot = u 0 1 + r 2 0 r 2 cos(ϑ), is set. The pressure is calculated from the Bernoulli theorem, the density by the adiabatic condition.
The 2D domain of size (10 m × 10 m) is discretized by 514 × 514 points and contains a cylinder of radius R 0 = 1 m. The wall of the cylinder is again smoothed by the same hyperbolic tangent as above with δ = 1.5Δx.
The initial condition is the velocity constant (u 0 , 0) (impulsively started cylinder). At the boundary, we force the analytical solution by a sponge layer of thickness 10Δx, and use a forcing strength for the Darcy and the sponge term of η = 2Δt.
Further, a filter of first/second-order as described above is used only inside the cylinder, again with the same mask function retreated by 4Δx inside the object. A Darcy term is used with a mask created by the same function as the volume fraction (28), where the reference pointx 0 is moved relative to the wall position by 7Δx so that its action is practically zero at the boundary. The purpose of both is to avoid dynamics inside the object, which otherwise destabilizes the simulation.
We find excellent agreement with the analytical solution, see Fig. 7. The numerical streamlines are on top of the analytical ones and the pattern of the velocity magnitude matches the expectations. The x component of the velocity along x = 5 also shows a very good quantitative agreement. A kink at the position where the sponge enforcing the boundaries ends, and a small zigzag mode relative to the analytical solution, is visible.
The differences might be explained by a slightly different blockage, due to the non-sharp boundary location. This interpretation is supported by the fact that the error becomes smaller with higher grid resolution, by which this effect becomes effectively smaller. It is concluded that slip conditions are created which are fully sufficient for practical calculations.

Wedge in Supersonic Flow
A wedge in supersonic flow is now investigated. The case follows [11] closely. Normalized quantities are used in this section. A two-dimensional domain of 2 × 2 with a wedge of half angel θ = 20 • and tip at S = (0.5, 1) is investigated. The main difference to the cited reference is a periodic boundary condition in y for simplicity. Due to the hyperbolicity of supersonic flows, large parts of the flow can be directly compared. The in-and outflow is forced by a sponge term to Ma = 2, the reference values are To handle the shocks, an adaptive filter as described above is used, where σ is determined by a shock detector as described by Bogey et al. [10] and modified as in [51] by a soft switch function. A threshold of r th = 10 −3 and a steepness of λ F = 1/10 was used, see Appendix C for more details.
Two different set-ups are compared. A non-slip boundary condition which compares to [11] and a slip condition which can be directly compared with analytical predictions from gas dynamics, where inviscid flow and thereby non-slip conditions are assumed. In both cases, a resolution of 1536 2 grid points as in [11] is used.
Unsurprisingly, the slip boundary condition, made possible by the new method, agrees much better with the gas dynamical prediction, see Fig. 8. The absence of numerical damping allows for fine details, but also yields high frequency noise, which is removed by filtering in post-processing. The error for the non-slip condition depends on the resolution since the minimal boundary layer thickness (created by numerical dissipation) depends on it. For a slip boundary, the agreement is largely independent of the resolution. This has also practical consequences, especially in a first design phase, where a low resolution is often interesting Journal of Scientific Computing (2022) 90:86 Fig. 9 The modulus of the density gradient (pseudo-schlieren) for the time t = 1.5 and t = 2.0 when the boundary layer is not expected to be crucial for the design. Here, a lower resolution permits a shorter design process, so that the volume fraction approach might be of help.

Wedge Shock Wave Interaction
Here, a shock-wedge interaction is discussed. The case follows a case presented by Kemm at al. [33] and Dumbser et al. [19]. A a wedge of height h = 1 and length l = 1 is placed with its tip at (2,3) in a domain of size (8 × 6). A shock wave is located at x = 1 at time t = 0, with a Mach number Ma= 1.3 enters a gas at rest (u 0 = 0) with ρ 0 = 1.4 and p 0 = 1 with a adiabatic index of γ = 1.4. The domain is resolved with 1000 × 750 points to be comparable with the additional degrees of freedom of discontinuous Galerkin elements of [32]. The reduced volume is again created from a signed-distance function and the smoothing by (28) with δ = Δx with a value of = 10 −8 in the object.
The shock filter with r th = 10 −4 and λ F = 4 is applied at every time step. Additionally, a global filter strength σ static = 0.05 is used to avoid noise in the simulation, see Appendix C for details. The first/second-order stencil is used.
The overall agreement is good. However, the shocks are not as sharp as in the reference simulation. This could be improved by a different shock detector (similar to a flux limiter) or by optimized filter coefficients. Artificial viscosity [15] might be an option. This investigation is not in the scope of this publication.

Pressure Pulse in Tube
This example is motivated by a detonation wave in a tube entering a plenum of a turbine and aims at assessing the conservation quality. Such unsteady flows arise when coupling detonative and unsteady combustion with turbines to construct new types of gas turbines, see e.g. [8]. Here, high-pressure gradients appear, which challenges purely a Darcy-frictionbased volume penalization, since it has to balance the high-pressure forces. Unlike a turbine, the domain is closed in the example, to allow a direct evaluation of conservation properties.
The two-dimensional domain is 2 m long and 1 m wide, with a tube of diameter 0.2 m in the left part and a chamber in the right part. All walls have a thickness of 0.05 m. Three NACA6412 profiles of cord length 0.3 m tilted by −(3/16)π are placed in the domain, As before, a Darcy term with strength 1/(2Δt) is used inside the wall with an offset of 7Δx to create slip walls. A first/second-order filter is used which is applied everywhere in the object and at shock position in the flow. For this, the shock detector described for the wedge flow is used again, with r th = 10 −5 .
A supersonic jet ejects from the tube and decays before fully established, due to the small volume of the tube. By this, high-pressure loads and a large unsteadiness are created on the profiles and the walls. A snapshot of the Mach-number and the density is shown for the times t = 2.67 ms inf Fig. 10. The total mass conservation is shown. If the density is integrated (summed-up over all grid points), near-perfect conservation is found, as expected from the conservative mass equation (4). Since we implement it in the form (21), a small error is introduced with an explicit time integration scheme (not strictly conserving quadratic invariants), which yields a relative error below 10 −7 over the full simulation with 20,000 steps. To integrate the fluid domain, one has to choose the exact location of the wall, due to the ambiguity introduced by the smoothing of φ. Here, this is done by choosing a threshold in φ, thus defining a sharp mask. The larger the threshold is chosen, the more of the boundary region is omitted and the lower the conservation quality becomes, see Fig. 10. However, the conservation error fluctuates around the correct value, since the mass is not lost but just pressed into the wall smoothing region. Thus, this behavior might be better viewed as a slightly elastic wall than as a mass conservation defect. In practice, the integration of the full domain might be appropriate since the mass deep inside the object becomes as vanishingly small as φ does. Altogether, we can confirm good conservation.

Vortex Shedding
This case aims at a non-slip boundary and transformed grids. Transformed grids allow focusing on location with high gradients in the flow and thereby reduce the total number of points and likewise the numerical effort. Since the volume fraction is part of the continuous equations with a simple physical interpretation, it is straightforward to generalize the scheme to transformed grids.
A flow with Reynolds number Re = (u 0 D)/ν = 150 around a cylinder with diameter D = 1 is simulated in a domain of size (20 × 15), using normalized quantities in this section. The in-flow velocity is u 0 = 30 with a speed of sound c 0 = 374 yielding a Mach-number of Ma = 0.08 so that it can be safely regarded as incompressible.
The transformed grid with N ξ × N η = (256 × 166) is given by Here, ξ = (L ξ +2α)(0, . . . , N ξ −1)/(N ξ −1) and η = (L η +2α)(0, . . . , N η −1)/(N η −1). The parameters are α = 7.2 and ψ = 8 and x δ = −0.55 From this, the coordinates are created by normalizing to the desired domain size. The resulting grid is shown in Fig. 11, where every fourth grid line is plotted. A CFL number of 1/2 was used. The grid transformation is introduced by rewriting divergence and gradient operators with the help of local bases, as discussed in the Appendix A. In the end, all divergence operators are replaced by the so-called conservative form and the gradient operators by the so-called non-conservative form, yielding a discretely consistent (and conservative) scheme, if all derivatives are replaced by finite difference approximations with a (skew-)symmetric stencil. Using summation by part operators results in well-defined fluxes at the boundaries result, see [51] for details. 5 The grid is orthogonal; a non-orthogonal grid can also be handled by this approach, it was however not explicitly tested.
The initial flow is the potential flow discussed in Sect. 4.2.1 with an additional Gaussian velocity defect at (x 0 = 0, y 0 = 1/2), with a width of σ disturb = 1/2 and an amplitude −20 to trigger the vortex shedding: The inflow and outflow are created by sponge terms enforcing the potential flow around a cylinder, where a quadratic sponge of thickness 10L ξ /N ξ is used. The sponge is visible by creating extra vorticity close to the outflow in Fig. 11. A non-reflecting boundary condition with reference conditions as the initial condition is used at the top and bottom.
The described filter of second/fourth-order is used inside the cylinder, with a small amplitude of 1/100 in the whole domain to avoid grid to grid oscillations, presumably created by the boundary conditions. The friction-like effect of the filter is by this small by the order 1/20 compared with the physical friction and should, thus, have little influence on the Reynolds number.
The expected vortex shedding is found after a short transition time with a Strouhal number of Sr = 0.173, which is in agreement with the literature. It is found to be robust towards simulation details like resolution. The smoothing of φ of δ = (L ξ /N ξ )/2 was used, which is close to one grid spacing at the cylinder wall due to the grid refinement. A sharper smoothing for the Darcy term was used with δ Darcy = δ/2, to ensure a quick decay of the tangential velocity, however, small sensitivity was found in these parameters. This is important since the good slip conditions possible by the volume fraction φ might be taken as an indicator of possible problems with the non-slip conditions. No signs for this are, however, found in the presented investigations.
Altogether, we find satisfactory results for this low Mach number, viscosity-driven case.

Conclusion and Outlook
A method to describe immersed objects was presented. It is easy to implement and has a clear physical interpretation. Furthermore, it has good conservation properties. It naturally allows implementing slip boundaries and adiabatic boundary conditions. The non-penetration condition can be implemented with only a minor increase in the stiffness of the equation and near-perfect pressure-tight boundaries are possible, even for high-pressure gradients. The derivation follows considerations of Kevlahan et al. [34] for shallow water equation and builds on a physical consistent treatment of the porosity, understood as a reduced volume fraction available for the fluid. By this, it becomes equivalent to the available cross-section.
While the main method works well, the filter, needed especially for the treatment of shocks needs manual adjustment and does not produce as sharp shocks as shock limiter methods [33]. This will be addressed in future works by optimizing the filter for the specific discretization of the Euler equations.
The method does not pose any extra difficulties in a parallel environment. For more complex geometries large simulations in 3D are in preparation. To avoid the increased stiffness for very sharp boundary functions, a better resolution close to the boundary is favorable, if such a high localization of the wall is demanded by the problem. Alternatively, to the shown local grid transformations, a local grid refinement as in [23] can be used, which is also under investigation. The possible non-slip conditions and the very good reflectivity make this equation a good candidate especially for aero-acoustic problems, which are currently further investigated.

Algorithm 1
The main solver q 0 = initialFlow() # inital condition for the case for n = 1 to N timeSteps do # Runge Kutta 4 time stepping: k 1 := rhs(q 0 ) k 2 := rhs(q 0 + Δt 2 k 1 ) k 3 := rhs(q 0 + Δt 2 k 2 ) k 4 := rhs(q 0 + Δt k 3 ) q 1 := q 0 + Δt(k 1 + 2k 2 + 2k 3 + k 4 )/6 # filter if demanded if filterAtStep(n) then σ α = calcFilterStrength(q 1 ) # adaptive filter strength for shocks for α = 1 to spatialDimensions do q 1 = q 1 + (D α (M(σ α φ))D α q 1 )/φ end for end if # plot or save state q 0 := q 1 # prepare next time step end for return From theses equations the time derivatives of the calculation variables directly follows, yielding the structure ∂ t q = rhs(q): All field quantities are assumed to be multiplied or divided pointwise. The heat fluxφ α was introduced, which is not to be confused with the reduced volume. For Cartesian grids J ≡ 1,ũ α = u α andD x α = D x α . The derivative matrices are central finite difference derivative matrices along one spatial direction, i.e., one index of the discretized field. In detail, assuming summing convention: A curvilinear grid is introduced by a transformation following Thompson et al. [59]. The transformation is implied by the discrete values of the coordinates at the grid points and skew symmetry as a result of the skew-symmetry of the one-dimensional (central) derivative. The non-periodic case can be treated with SBP matrices in a fully consistent manner (as in [51]), but is not detailed here for the sake of brevity. Viscous terms are omitted in the discussion but can be included in a straightforward manner. The conservation of the total mass M follows directly by summing the equation 1 T (34) from which the spatial part evaluates to zero due to (44) and the temporal part becomes with the product rule 8 The change of total momentum is calculated from 1 2 u T α (34) + (35) : The telescoping sum property in the first term of the momentum transport and the combination becomes zero by the skew-symmetry of the discrete derivative since (φρũ β ) T D ξ β u α = (φρũ β ) T D ξ β u α T = u T α D T ξ β (φρũ β ) = −u T α D ξ β (φρũ β ). The first term in (47) is the temporal change of the total momentum the second is the source expected due to reduced effective volume by the embedded objects; it is a discrete form of the integral Ω p ∂ x α φ ∼ 1 TD x α φ p − 1 T φD x α p, since the first term is zero by the telescoping sum property. The term proportional to χ is the momentum source due to the Darcy friction term. All momentum changes are thereby a direct consequence of the source terms modeling the objects.
The energy conservation is obtained by the combination of the equation of internal energy and kinetic energy. The latter is calculated by u T α (35): The second term, the transport, can be regarded as a quadratic form u T α D ρφu u α over the skew symmetric matrix D ρφu = (D ξ β {ũ β ρφ} + {φρũ β }D ξ β ) = −(D ρφu ) T , which is zero in general. The skew symmetry becomes clear if one considers the pointwise multiplications in the operator as diagonal matrix {φρũ β } with the corresponding values on the diagonal. This skew-symmetry is the central point of the method as it avoids any numerical change of the kinetic energy. In most schemes, the momentum equation can change the kinetic energy, which often demands the introduction of numerical dissipation even for smooth flow fields for stability reasons.
Calculating the internal energy part from 1 T (36) , where the divergence terms again evaluate to zero, combines with (49) to yield The pressure terms in kinetic and internal energy couple the kinetic and internal energy and cancel, as in the analytical case. The Darcy term cancels if it is included in the internal energy equation (36). Since the change of kinetic energy by the Darcy term in the momentum equation might be compensated by the object, it is justified modeling that the term proportional to χ is omitted in (36) and the total energy is not constant. This can be interpreted as a penalization term for the energy as used in [11] and is done in the simulations presented in this report. The term could easily be included if strict energy conservation is desired.
The boundary conditions at the domain boundary are either periodic boundaries, implemented by a periodic derivative matrix in this direction, or are created by quadratic sponges as discussed in [41] or are characteristic boundaries where a decoupling in Riemann waves is done with respect to a reference state and the in-going waves are set to zero, as discussed in the appendix of [30] with a system matrix in a direction perpendicular to the boundary, see [29, sec 16.5] for details.
The product rule, which was used for the time derivative, is violated in most time discretizations. The conservation of quadratic invariants, needed in this approach for strict conservation, are respected by the implicit Gauss-Lobatto methods [12]. The special choice of variables allows to use these methods and directly guarantees strict conservation, however, the error with the RK4 is so small that the numerical effort for the implicit time integration does not seem justified in most cases [12]. For example, the relative error of the conservation mass in the full domain of the pulse tube test case was below 10 −7 for the mass, which seems sufficient for the majority of practical studies.

C Filter Details
Numerical dissipation is introduced by a filter in the presented report. It ensures a sufficient amount of dissipation for shocks and allows to suppress highly oscillatory behavior, which can be created by the nearly non-smooth φ within the objects by waves entering from the physical domain. The basic properties of the filter are discussed in Sect. 3.
Here, details of the shock detector are provided, which is a slightly modified version of the detector discussed in [10]. The core idea is to create a shock detector to locally activate filtering. The filter can also be activated within the objects or (e.g. with a small strength) in the full region.
The shock detection is based on the local divergence Θ = div(φ(u)), from which the non-smoothness is calculated by an expression similar to a second derivative (D ξ α Θ) i, j,k = (−Θ i+1, j.k + 2Θ i+1, j,k − Θ i−1, j,k )/4, from which a magnitude is calculated by and likewise in the other spatial directions. Other non-smoothness criteria are possible for example, a adaptive filter based on the WENO detector is presented in [58], different detectors are compared in [48]. The non-smoothness parameter (D ξ α Θ) magn is compared with the grid spacing and the local speed of sound r α = (D ξα Θ) magn c 2 /Δξ 2 α from which finally the filter strength is calculated as Here, r th and λ F are adjustable parameters. The procedure is similar to a flux limiter, which create in a non-linear manner a local amount of dissipation. Flux limiters usually go from high order schemes to low order up-wind methods, while here a purely dissipative term is added. The functional form (54) is different from [10], where a non-smooth function is used, creating a by switching sometimes high-frequency noise in the solution. This adaptive part can be combined with a user-defined static (and locally varying) filtering strength by σ α = max(σ shock α , σ static ). The static part is for example a filter acting within the solid objects. The construction ensures that 0 ≤ σ α ≤ 1, ensuring dissipative behavior.