*3.1. One-Dimensional Acoustic Wave Reflection*

To assess how well the penalization scheme can capture the reflective nature of a solid wall, we used the reflection of an acoustic wave at the material. The initial pressure distribution is shown in Figure 1. It is described by the Gaussian pulse given in Equation (35) around its center at *x* = 0.25 in the left half of the domain (*x* ≤ 0.5).

$$p' = \mathfrak{u}' = p' = \epsilon \exp\left[-\ln(2)\frac{(\mathfrak{x} - 0.25)^2}{0.004}\right] \tag{35}$$

For the amplitude of the wave, we used a value of = 10−3. The perturbations in density *ρ* , velocity *u* , and pressure *p* from (35) are applied to a constant, non-dimensionalized state with a speed of sound of one. This results in the initial condition for the conservative variables density *ρ*, momentum *m*, and total energy *e* as described in: (36).

**Figure 1.** One-dimensional acoustic wave setup: the center of the initial pressure pulse is located at *x* = 0.25 and has an amplitude of = 10−3. Discretization by 48 elements as denoted by grid lines, and the right half of the domain (*x* > 0.5) is penalized. Note that the wall coincides with an element interface.

The penalization with porous medium is applied in the right half of the domain (*x* > 0.5). In acoustic theory, the reflection should be perfectly symmetric, and the reflected pulse should have the same shape and size, only with opposite velocity. This simple setup allows us to analyze the dampening of the reflected wave and induced phase errors. Reflected waves for different settings of *β* as defined in Equation (33) are shown in Figure 2. The pressure distribution for the reflection is shown for the state after a simulation time of 0.5. With linear acoustic wave transport and a speed of sound of one, the pulse should return to its original starting point, just with an opposite traveling direction. This symmetry makes it easy to judge both the loss in wave amplitude and the phase shift of the reflected pulse.

While the analytical result for a linear wave transport provided a good reference in general for the acoustic wave, it sufficiently deviated from the nonlinear behavior to limit its suitability for convergence analysis to small error values. Therefore, we compared the simulations with the penalization method to numerical results with traditional wall boundary conditions and a high resolution. This reference was computed with the same element length, but the domain ended at *x* = 0.5 with a wall boundary condition, and a maximal polynomial degree of 255 was used (256 degrees of freedom per element) to approximate the smooth solution. The resulting pressure profiles for different settings of *β* and a fixed porosity of *φ* = 1.0 are shown in Figure 2. This illustrates how well the wave was reflected for different settings of *β* and that the solid wall reflection was well approximated for sufficiently small values of *β*. These numerical results were obtained with 48 elements and a maximal polynomial degree of 31 (32 degrees of freedom per element). Note that this setup aligns the wall interface with an element interface, where the discontinuity in the penalization is actually allowed by the numerical scheme. Later, we will discuss the changes observed, when moving the wall surface into the elements.

**Figure 2.** Plot for the pressure profile of the reflected wave at *t* = 0.5 for different scaling factors *β*. The numerical reference is obtained with a traditional wall boundary condition and a high resolution.

Figure 3 illustrates the impact of the porosity on the error in amplitude of the reflected wave for the same discretization with 48 elements and a maximal polynomial degree of 31. We plotted the error *e* = ( − *p* (*t* = 0.5))/ over porosity *φ* for various scaling parameters of *β* between one and 10−6. A scaling parameter of *β* = 1 means that the error is only driven by the porosity *φ*, and for large values of *<sup>β</sup>* ≥ <sup>10</sup><sup>−</sup>2, we observed the expected reduction in the error with decreasing porosity. However, with *<sup>β</sup>* = <sup>10</sup><sup>−</sup>3, this comes eventually to an end (no improvements for *<sup>φ</sup>* < <sup>2</sup> × <sup>10</sup><sup>−</sup>2), and for smaller values of *β*, no improvements for the error can be achieved by lowering the porosity anymore. As can be seen in this figure, a sufficiently small permeability can yield the accuracy as a small porosity.

**Figure 3.** Plot of the error in the wave amplitude at *t* = 0.5 with decreasing porosity and different scaling factors *β*. The error *e* is given by the relative error in the pulse amplitude after the reflection at the wall.

In Figure 3 as well, for *β* = 1*e*<sup>−</sup>3, we observed a drop in error, and then it increased again with smaller *φ* to reach a convergence point finally. We would like to point out that this was expected to come from our numerical scheme using polynomials to represent the solution. Each data-point in the plot with varying *φ* and *β* represents a slightly different test case in terms of boundary layer thickness, as pointed out in Section 2.2. A sweet spot is reached when the degree of polynomial used for the simulation correctly captures the boundary layer in the problem. However, as we move further left from here, this sweet spot is slowly gone with further thinning of the boundary layer. With the same polynomial degree used, one would also expect to see this behavior for lines representing *β* < 1*e*−<sup>3</sup> and correspondingly larger *φ*. This is exactly what we also see for *β* = 1*e*−<sup>4</sup> and the value of *φ* close to 1.0. For all other lines in the plot, this spot does not fall within the range of the figure.

By using the implicit mixed explicit scheme from Section 2.4, it is possible to exploit drastically smaller values for the permeabilities to cover up the lack of porosity in the penalization. On the other hand, the porosity cannot that easily be treated in our discretization, and even moderate values of *φ* can have a dramatic impact on the time step restriction, due to the changed eigenvalues in the hyperbolic part of the equations.

Next, we performed a convergence analysis shifting the position of the wall such that it intersected the element at different locations. The reason for performing such an analysis becomes imperative when we consider the high-order numerical scheme used. We represented the solution state within an element using polynomials. For the pointwise evaluation of the nonlinear terms, we employed the Gaussian integration points, at which also the masking function of the penalization needs to be evaluated. Within an element, these integration points were scattered, being more concentrated on the element interface, and were rather sparse at the center. Therefore, in a peculiar validation test case like this one, when the wall was aligned with the element interface, due to the abundance of interpolation points, it had the advantage of being very precisely represented even for comparatively fewer degrees of freedom. In actual simulations, the wall interface may intersect an element at any point. We, therefore, also need to consider such intersections through the element and ensure the solution converges to the reference solution. The penalization method itself was not restrictive to any such limitations and could perfectly represent wall irrespective of its location within an element.

Thus, we performed and compared convergence analysis on two different discretizations, one where the wall lied at the element interface and a second where the wall intersected one element exactly in the middle. We would like to point out that the later scenario yielded a worst case estimate for the approximation of the jump in the masking function within an element. As explained, this simply came from scarce integration points lying around the center of element. For the following convergence analysis, we ignored the porosity (i.e., set *φ* = 1) and used small permeabilities by choosing *β* = 10<sup>−</sup>6. We also considered the *L*<sup>2</sup> error norm now in the fluid domain. As a reference solution, we employed a numerical simulation with a traditional wall boundary and a high maximal polynomial degree of 255 (256 degrees of freedom per element). The error was measured at *t* = 0.5 after the reflected wave reached its initial position again.

Figure 4 shows the *L*<sup>2</sup> norm of the error for the reflected pressure wave with a maximal polynomial degree of seven over an increasing number of elements (h-refinement). This plot compares the two discretizations explained. As expected, in Figure 4, we observe superior convergence behavior for the case when the wall lies at the element interface in comparison to the other case where the wall is crossing through the element center. However, for both cases, we observed a proper convergence towards the solution with a traditional solid wall boundary condition. The order of error convergence did not match the high-order discretization in either case, but this was expected due to the discontinuity introduced by the masking function of the penalization.

**Figure 4.** *L*2-error for a polynomial degree of seven over an increasing number of elements (h-refinement).

Next, we performed another convergence study using the same two discretizations, but this time keeping the number of elements constant and increasing the order of polynomial representation within those elements.

Figure 5 shows the error convergence over the maximal polynomial degree in the discretization scheme (p-refinement) with the number of elements fixed to 24. Here, also, one observes a solution in both cases converging to the reference solution. While no spectral convergence was achieved for the discontinuous problem, a quadratic convergence can be observed. This shows that a high-order approximation was beneficial even with the discontinuous masking function for the penalization.

**Figure 5.** *L*2-error for 24 elements over increasing maximal polynomial degree (p-refinement).

Finally, we also looked at the case where the wall was close to the element interface, but not exactly on it. This is a potential critical configuration for the numerical scheme being used, as the discontinuity close to the surface needs to be properly captured. We put the wall at 5% of the element length away from the element surface and measured the error as before, resulting in the graph shown in Figure 6. For this case as well, we observed a similar convergence rate as before, though the error was a little bit worse than with the wall on the interface.

For a smooth solution, the advantage of high-order methods to attain a numerical solution of a given quality using fewer degrees of freedom is well documented [16]. However, for a complex nonlinear problem with a discontinuity introduced by the porous medium, it is not so clear whether there is still a computational advantage by a high-order discretization. To investigate this, we ran the wall reflection problem for several orders and plotted the convergence with respect to the required computational effort, as seen in Figure 7. This test was performed starting with 16 elements in each data series, providing the leftmost point for the respective spatial scheme order. For subsequent data points, the number of elements was always increased by a factor of two up to the point where an error of 10−<sup>6</sup> was achieved. Figure 7a depicts the observed *L*<sup>2</sup> error over the total number of degrees of freedom in the simulation. Here, we can see that for attaining a certain level of accuracy, the number of degrees of freedom required was always less when using a higher spatial order, even though the convergence rate did not increase with the scheme order. The high-order discretization, thus, allowed for memory-efficient computations, also in this case with a discontinuity present at the wall.

**Figure 6.** *L*2-error for varying the polynomial degree. With a wall just 5% of the element length away from the element surface.

**Figure 7.** Behavior of the error in the reflected acoustic pulse with respect to computational effort. The figure on the left (**a**) shows the error convergence for various spatial orders over the required memory in terms of degrees of freedom. The right figure (**b**) shows the same runs, but now over the computational effort in terms of running time in seconds. All simulations were performed on a single node with 12 cores using 12 processes.

While for the memory consumption, there seems to be a clear benefit in high-order discretizations, it is not so clear whether this still holds for the required computing time. The time step limitation of the explicit scheme required more time steps for higher spatial scheme orders, increasing the

computational effort to reach the desired simulation time. Additionally, the number of operations increased with higher orders due to the nonlinearity of the equations. Figure 7b shows the measured running times on a single computing node with 12 Intel Sandy-Bridge cores for the same runs. Again, the achieved accuracy is plotted, but this time over the observed running time in seconds. As can be seen, the advantage in terms of running times was not as clear as in terms of memory. However, we still observed faster times to the solution with higher spatial scheme orders, despite the increased number of time steps. In conclusion, we found some computational benefit from higher spatial scheme orders even in the presence of a discontinuity in this setup.

#### *3.2. One-Dimensional Shock Reflection*

After considering the reflection of an essentially linear acoustic wave, we now look into the reflection of a shock, where nonlinear terms play an important role. However, we neglected viscosity in this setup and only solved the inviscid Euler equations. The reflection of a one-dimensional shock wave at a wall was described and numerically investigated by Piquet et al. [17], for example. We used their setup to validate the penalization method in our discontinuous Galerkin setup, even though a high-order scheme is not ideal for the representation of shocks.

The downstream state in front of the shock (denoted by 1) is given in Table 1. The upstream state after the shock (denoted by 2) is then given by the Rankine–Hugoniot conditions for the shock Mach number *Mas*. These yield:

$$\frac{\rho\_2}{\rho\_1} = \frac{\gamma + 1}{\gamma - 1 + 2Ma\_s^{-2}}\tag{37}$$

for the relation of densities *ρ* in up- and downstream of the shock and

$$\frac{p\_2}{p\_1} = \frac{2\gamma Ma\_s^2 - (\gamma - 1)}{\gamma + 1} \tag{38}$$

for the relation of pressures *p*.

With these relations, the ratio of the upstream (*p*3) and downstream (*p*2) pressure for the reflected shock wave is [18]:

$$\frac{p\_3}{p\_2} = \frac{Ma\_s^2(3\gamma - 1) - 2(\gamma - 1)}{2 + Ma\_s^2(\gamma - 1)}\tag{39}$$

For the computation of the velocity *urs* of the reflected shock wave, we considered Equation (40) [19].

$$u\_{\rm r\S} = \frac{1}{Ma\_{\rm s}} \left( 1 + \frac{2(Ma\_{\rm s}^2 - 1)}{(\gamma + 1)/(\gamma - 1)} \right) c\_1 \tag{40}$$

**Table 1.** Shock state description.


With an incident shock wave velocity of *Mas* = 1.2, we obtained a pressure relation across the shock of 1.47826 from Equation (39). The shock was simulated in the unit interval *x* ∈ [0, 1] with the wall located at *x* = 0.5. Thus, half of the domain (*x* ∈ [0.5, 1]) was covered by the porous material to model the solid wall. The shock was initially located at *x* = 0.25.

For the numerical discretization, we used 256, 512, 1024, and 2048 elements (*n*) in total (Δ*x* = 1/*n*) and a scheme Order (O) of 32, 16, 8, and 4, respectively. As explained in the inspection of the linear wave transport, our numerical scheme preferred strong permeabilities over penalization with the porosity. Therefore, we ignored porosity and chose *φ* = 1, while the scaling factor from Equations (33) and (34) was set to a small value of *β* = 10<sup>−</sup>6.

In Figure 8, the shock wave after its reflection is shown for different spatial resolutions. The discretizations with O(8) and 1024 elements and with O(4) and 2048 elements were chosen to have the same number of degrees of freedom, while the third discretization with O(8) and 2000 elements provided a high-resolution comparison.

**Figure 8.** Different curves represent different discretizations using different scheme orders and a different number of elements. (**a**) Normalized pressure of the reflected shock wave. (**b**) Zoom of the reflected shock.

The exact solution for the normalized pressure (*p*3/*p*2) according to Equation (39) was 1.47826087. In Table 2, the ratio of the pressure (*p*3/*p*2), the relative error between the numerical, and the exact solution (error in *p*3/*p*<sup>2</sup> in %) close to the shock, as well as the difference between the location of the shock wave after the reflection and its origin (Δ*x*: phase shift) are listed. The table illustrates that with higher scheme order, but constant number of degrees of freedom, the error in the pressure ratio, as well as in the phase shift reduces considerably even for this discontinuous solution. From the obtained results, we can conclude that we achieved the same error as in [17], when using O(16) and 512 elements. As can be seen in Figure 8b, the plateau after the shock was not fully flat, but rather had a slope that asymptotically got close to the expected constant value. Except for the fourth-order approximation, this constant plateau was well obtained, but it remained slightly off the exact solution. This remaining error was also stated in the table as *minerror* and had a value of 0.0129% for *β* = 10<sup>−</sup>6.


**Table 2.** Comparison of simulation results with the exact solution.

In Table 3, the results for the reflected shock wave are presented for the case that the wall is inside an element, instead of at its edge. Again, the error was reduced by an increased scheme order and a fixed total number of degrees of freedom. Notably, the error in the pressure ratio was reduced for small element counts in relation to the case where the wall coincided with an element interface. This was due to the fact that there was an additional element introduced here, and the element length was accordingly smaller. However, we can see that the phase shift of the shock was larger in this case. This can be attributed to the larger distance of the Gaussian integration points, which were used to represent the wall interface.

**Table 3.** Comparison of simulation results, when the porous material was located in the middle of the element with the exact solution.


For better comparison of the simulation results, Figure 9 illustrates the different test cases. The plot presents the solution, from the previous investigation, when using a scheme order of O(16). As a reference, we considered a no-slip wall, which was located at the same place as the porous material, while considering the same scheme order.

**Figure 9.** Different curves represent different locations of the porous material in the element and the solution when using a no-slip wall. (**a**) illustrates the normalized pressure of the reflected shock wave, and (**b**) depicts a zoom-in of the front area of the reflected shock.

Figure 9 accentuates the close solution, when modeling the wall as a porous material, to the reference wall. As can be seen in Figure 9a, the high-order discretization introduced Gibbs oscillations around the shock, but otherwise, the discontinuity was well preserved by the numerical scheme. Some of those oscillations remained inside the modeled material (see Figure 9b) for the solid wall, but again, the discontinuity was well preserved. Further, the reflected shock exhibited an over- and under-shoot. Since the material was represented in polynomial space and according to the Gibbs phenomenon, we were limited to 9% deviation for physical correctness of the solution, we also computed those over- and under-shoots. For the material located inside an element, the overshoot was around 2.052% and the undershoot around 8.639%. Locating the porous material at the element interface resulted in 1.821% and 7.681%, respectively.
