Compressible Flow Simulation with Moving Geometries using the Brinkman Penalization in High-Order Discontinuous Galerkin

n this work we investigate the Brinkman volume penalization technique in the context of a high-order Discontinous Galerkin method to model moving wall boundaries for compressible ﬂuid ﬂow simulations. High-order approximations are especially of interest as they require few degrees of freedom to represent smooth solutions accurately. This reduced memory consumption is attractive on modern computing systems where the memory bandwidth is a limiting factor. Due to their low dissipation and dispersion they are also of particular interest for aeroacoustic problems. However, a major problem for the high-order discretization is the ap-propriate representation of wall geometries. In this work we look at the Brinkman penalization technique, which addresses this problem and allows the representation of geometries without modifying the computational mesh. The geometry is modelled as an artiﬁcial porous medium and embedded in the equations. As the mesh is independent of the geometry with this method, it is not only well suited for high-order discretizations but also for problems where the obstacles are moving. We look into the deployment of this strategy by brieﬂy discussing the Brinkman penalization technique and its application in our solver and investigate its behavior in fundamental one-dimensional setups, such as shock reﬂection at a moving wall and the even with discontinuities present in the solution. It also shows, that the Brinkman penalization can successfully be used in the supersonic regime.


Introduction
In engineering applications we are often dealing with moving parts. When investigating the fluid motion in those scenarios we, therefore, need to consider complex, moving obstacles that influence or even drive the flow. In this work we look into a method to realize moving obstacles in simulations with a high-order discontinuous Galerkin scheme. High-order schemes need less number of degrees of freedom to represent smooth solutions. On modern computing systems, where memory bandwidth is one of the main bottlenecks, this is a big advantage. The discontinuous Galerkin method combines a high-order local state representation within elements with a coupling of elements via fluxes on the surface as in the finite volume method. This relatively loose coupling via the element surfaces is a strong advantage on modern distributed parallel computing systems, as only comparable little data needs to be communicated between elements.
Obstacles are in such mesh-based schemes usually represented by fitting the mesh to the geometry and prescribed via boundary conditions. An alternative option is the use of some immersed boundary method [17]. Immersed boundaries have the advantage that they do not require the mesh to be fitted to the geometries that are to be represented. This is especially attractive when considering moving obstacles as the movement does not necessitate a change of the mesh.
In high-order discretizations it also becomes difficult to achieve an accurate representation of the geometry by fitted meshes, due to the large elements that we want to employ. With moving geometries the problem of fitting the mesh of a high-order discretization to complex geometries gets magnified. Thus, immersed boundaries are a promising choice for moving obstacles in a high-order discontinuous Galerkin scheme. Comprehensive reviews of immersed boundary methods are available by Mittal and Iccarino [17], Sotiropoulos and Yang [22] and recently Maxey [16]. We use a penalization method as described by Liu and Vasilyev [15].
Volume penalization methods introduce penalty terms in the continuous equations, while imposing nothing in the discretization context [4]. Arquis and Caltagirone [5] introduced a type of volume penalization method for simulations of isothermal obstacles in incompressible flows that makes use of the Brinkman model for porous media. The basic idea is to model obstacles as a porous media, where the porous medium approaches a solid obstacle. A distinct advantage this method, in comparison to other penalization methods, is its error estimation, which can be rigorously predicted in terms of the penalization parameters [19]. Kevlahan and Ghidaglia [11] used a pseudo-spectral method and applied this Brinkman penalization method to model stationary, as well as moving geometries.
Title Suppressed Due to Excessive Length 3 Liu and Vasilyev [15] used the Brinkmann penalization method for the compressible Navier-Stokes equations. They introduced new penalization term for the mass conservation and showed promising results for their acoustic setups. Several other successful investigations using different numerical methods speak for the effectiveness of the Brinkmann penalization method. Application in pseudo-spectral methods are found in [18] and [10]. In [19] the Brinkman penalization is used in the context of finite volumes and finite elements.
In [2] we applied the Brinkman penalization in our high-order discontinuous Galerkin method for compressible NavierStokes equations with non-moving obstacles. That analysis showed that the method can be successfully used to represent obstacles even when using high-order schemes and coarse elements. Thus, we now look into the application of this scheme to moving geometries. Our implementation is available in our open-source solver Ateles [21].

Numerical Method
In the following we revisit the fundamental equations of compressible viscous flows. With the continuum assumption fluid motion is modeled in terms of the compressible Navier-Stokes equations. To model obstacles as immersed boundaries, we extend those equations by Brinkman penalization terms as described by Komatsu et al. [13]. With the introduction of the penalization terms we also discuss their treatment in the time integration scheme.

Compressible Navier-Stokes Equations
Conversation of mass, momentum and energy applied to fluids leads to the equation system that we refer to as compressible Navier-Stokes equations. In conservative formulation with density ρ, momentum m and energy E the system can be written as: Where τ is the viscous stress tensor defined by (2).
Equation (1a) represents conservation of mass, (1b) conservation of momentum and (1c) conservation of energy. Velocity is denoted by u = m/ρ and . I indicates the identity matrix. The parameter µ represents the dynamic viscosity of the fluid and λ the volume viscosity, assumed as λ = 2µ/3 here. The temperature of the fluid is denoted by T and κ represents its thermal conductivity. Additional source terms, like gravity, are represented on the right hand side of (1b) and (1c) by f . The equation system is closed the equation of state for the ideal gas that provides a relation between pressure p, temperature and density: With γ = c p /c v denoting the heat capacity ratio and R the ideal gas constant as given in (4).

Brinkman Penalization
The volume penalization methods introduce additional artificial terms to the equations to be solved in areas where obstacles are to be present. With the Brinkman penalization method these appear as source terms on the right hand sides of the equations. In [15] the Brinkman penalization method was introduced for compressible Navier-Stokes equations. While these included the possibility for moving obstacles, they were not Galilean invariant. This was corrected by Komatsu et al. [13], resulting in the formulation we use in this scheme. The equation system with the penalization terms is given in (5).
Where η and η τ are the viscous and thermal permeability and φ the porosity of the Brinkman model. Velocity and temperature of the obstacle is given by U o and T o , respectively and its location is defined by the masking function χ. While the permeabilities only appear in source terms, the porosity affects the divergence oper-Title Suppressed Due to Excessive Length 5 ators and thereby changes the Eigenvalues of the convective system and necessitate smaller timesteps when considered in explicit time integration schemes. As we show in [2], it is feasible to ignore the porosity, introduced by Liu and Vasilyev [15] for compressible fluids, if the permeabilities are chosen sufficiently small. While the source terms get extremely stiff by tiny permeabilities, this allows us to avoid any complications by the porosity in the spatial divergence operators. The source terms are purely local and can be solved for pointwise. Thus, a implicit mixed explicit time integration is deployed to deal with the source terms implicitly while the remaining parts are still solved explicitly. We ignore the porosity by choosing φ = 1, which results in the following, simplified model: can be simplified to To solve the source terms introduced by this penalization implicitly, we use a diagonally implicit Runge-Kutta scheme [1]. In this work we make use of the values for the permeability proposed in our analysis of the method for non-moving obstacles, see also there for details on the time integration method [2]. With the viscous permeability η = β 2 · φ 2 and the thermal viscosity η T = 0.4β · φ , the expected modeling error as described in [15] has a magnitude of β 1/4 for a porosity of 1. For our simulations we consider the outcomes of Anand et al. and consider a scaling factor of β = 10 −6 .

Representation of immersed boundaries
With the penalization method to represent obstacles, we can use the same discretization technique for the obstacles as for the solution in our scheme. Thus curvature of the geometries to represent poses no special problem as the discretization of the obstacles can be chosen from the same realm as the solution. In the penalization method, obstacles are given by the masking function χ. This function is one wherever an obstacle is to be modeled and zero everywhere else.
To discuss the representation of χ in our high-order discontinuous Galerkin implementation on the basis of Legendre polynomials we consider a simple one dimensional setup with a single element. Figure 1a shows this single element and the masking function that models a wall within that element. To the left with χ = 0 the fluid is to behave undisturbed, while to the right with χ = 1 the flow is to be penalized. This step function is mapped to the polynomial space and use a Legendre expansion with a maximal polynomial degree of 15. The blue line illustrates the projection with an analytic integration, while the black line indicates the resulting polynomial from a numerical integration with the shown 16 (Chebyshev) integration points. With the numerical integration an error introduced, as the jump of the masking function χ can only be determined up to the accuracy of the integration points. This is illustrated more clearly in figure 1b, which also clearly shows how this integration error depends on the location of the discontinuity in relation to the integration points. When considering a continuously moving wall, we will, therefore, see discrete jumps in the numerical representation of χ as the discontinuity moves across the integration point intervals. We further note that, as the integration points are not equidistantly spaced, this localization error is not constant, but depends on its location within the element. Closer to the element boundaries the localization is more accurate than in the middle of elements.
To decrease this error from numerical integration, we utilize over-integration and use more integration points for the approximation [12]. In Figure 2 we consider an over-integration factor of three. Comparing Figure 1b without over-integration and the same case with this over-integration by a factor of three in Figure 2, the improved representation of the wall location becomes apparent. As Figure 2 shows, an over-integration by factor three already yields a polynomial approximation that is close to the analytical integration.

Moving Wall Test case
To investigate the penalization scheme for a moving geometry in our discontinuous Galerkin implementation, we first analyze the fundamental behavior in a onedimensional setup. We look into the reflection of an acoustic wave from a moving wall. Figure 3 sketches this setup with the initial wave an wall positions.
The traveling wave has the form of a Gaussian pulse given by equation (7). Its center is initially located at x = 0.25, while the wall is located at x = 0.5.
The amplitude of the pulse is chosen to be ε = 10 −3 . The perturbations in density ρ ′ , velocity u ′ and pressure p ′ from (7) are applied to a constant, nondimensionalized state with a speed of sound of 1. This results in the initial condition for the conservative variables density ρ, momentum m and total energy e as described in (8).
where v is the background velocity equal to the velocity of the moving wall. 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. After the simulation time of t = 0.5, with linear acoustic wave transport and a speed of sound of 1, the position of the pulse will again be at a distance of d = 0.25 away from the position of the wall at that time. During this duration, the wall travels a distance of vt to the right, reaching the point x = 0.5. Therefore, the final position of the pulse has its center x = 0.25 after the time interval t. For the linear model we can easily compute the analytic solution and use this to judge the quality of the numerical computation with the penalization method. This simple setup allows us to analyze the dampening of the reflected wave amplitude as well as induced phase shifts. While the analytical result for linear wave transports provides a good reference in general for the acoustic wave with small perturbances, it sufficiently deviates from the nonlinear behavior to limit its suitability for convergence analysis in small error values. Therefore, we compare the simulation results with the penalization method and moving walls to numerical results with a traditional fixed wall boundary condition and a high resolution. This reference is computed with the same element length, but the domain ends at x = 0.5, which is the final position of the wall after t = 0.5 second, with a wall boundary condition and a maximal polynomial degree of 255 is used (256 degrees of freedom per element) to approximate the smooth solution. The pressure perturbation and the final position of the reflected wave (i.e at time, t = 0.5), for different orders, is shown in figure 4. The domain for all the simulations shown in the figure comprises of 12 elements, i.e dx = 1/12. Here, we can observe that both loss in wave amplitude as well as phase shift of the reflected wave gets diminished with increasing orders. With sufficiently large degrees of freedom, the wave shows an excellent agreement with the reference solution. The error convergence is done in reference to this solution.   Figure 5 shows the error convergence in terms of the L 2 norm over the maximal polynomial degree in the discretization scheme (p-refinement) with a fixed number of 12 elements. We observe a linear convergence rate with increasing polynomial degree, which matches the worst case scenario from the analysis with a non-moving wall in [2]. As in the moving case the discontinuity of the masking function χ ranges through all locations of the elements, we can not expect a better convergence rate.    Figure 6 shows the L 2 norm of the error for the reflected pressure wave with a maximal polynomial degree of 7 over an increasing number of elements (hrefinement). Again a first order convergence can be observed in the overall domain. The moving discontinuity in the masking function eliminates the fast error convergence rate of the high-order discretization close to the discontinuity. However, away from the discontinuity we maintain the fast convergence rate and still expect benefits.
In Figure 7 the computational effort is shown for various spatial scheme orders. Figure 7a shows the error over the used number of degrees of freedom and, thus the required memory consumption. And Figure 7b shows the respective running times of the simulations on a single computing node with 12 Intel Sandy-Bridge cores. The test was performed starting with 12 elements in each data series, which is the leftmost point in the plot. For subsequent points in the series, the number of elements are always increased by factor of 2. While we do not achieve a spectral convergence as in a purely smooth solution, we still find the higher order discretizations to be advantageous with respect to the required memory. In the computational effort a strong improvement from second to fourth order discretization is observed. However, for higher spatial scheme orders than 16, the required computational running time increases again. Though this specific timing result is tied to the system, the simulations were run on, we generally expect the running time to be worse for higher orders in this metric on other computing systems as well.
There are mainly two contributions to this error that diminishes the convergence order. One cause is the Gibbs phenomenon that describes induced oscillations due to the discontinuity in the masking function and the other is the accuracy of the localization of the discontinuity from the numerical integration. Both effects are illustrated in section 3. The inaccuracy from numerical integration can be limited by over-integration and using more points. The error from the Gibbs phenomenon can also be limited by employing a reprojection method [8]. Though, such a reprojection method potentially could restore the convergence up to the point of discontinuity, we do not consider it here.

Reflected shock wave
While in the previous setup we looked into a smooth acoustic wave, we now look into the reflection of a shock. We thus, not only have a discontinuity in the masking function χ but also in the solution itself. For the inviscid equations the solution to this problem can be computed analytically. Thus, we will neglect viscous terms in this setup and use the inviscid Euler equations. However, the same Brinkman penalization is used to model the moving wall. The reflection of a shock at a fixed wall with this method has been studied in [2]. Now, we move the wall similarly to the reflection of the acoustic pulse in the previous section. Of course the shock incurs more problems with a high-order discretization, nevertheless, this analysis shows that it is feasible to use it in this setting.
The state variables downstream (indicated by 1) are listed in table 1. The upstream state is computed via the Rankine-Hugoniot conditions for shock Mach number M s . The Rankine-Hugoniot conditions yield the following relations: for the relation of the pressure p before and after the shock and for the relation of the densities ρ. From these we can also find the pressure relation of before and after (p 3 ) the reflected shock, see [3]. Additionally, the wall velocity w s needs to be taken into account here, leading to equation (11) for the pressure relation.
The speed of the reflected shock wave is than given by equation (12) [7].
These relations provide us with the piecewise constant exact solution, that we use as a reference to compare the numerical results with. With a shock wave Mach number of M s = 1.2 and a wall speed of w s = 0.012 the exact solution for the pressure relation across the reflected shock wave is according to equation (11)  For our investigation, we consider four different numerical discretization, 256, 512, 1024, and 2048 elements (n) in total (∆ x = 1/n) and a scheme Order (O) of 32, 16, 8, and 4, respectively. This yields the same number of degrees of freedom (8196) across all runs. The penalization parameters are chosen to be φ = 1 and β = 10 −6 , the β value can be chosen that large as the numerical error is higher, when compared to the modelling error. Figure 8 presents the reflected shock wave for the different spatial resolutions. On the left in Figure 8a the overall pressure distribution is shown along with the piecewise constant reference. As can be seen, all discretizations with the same number of degrees of freedom provide a good approximation of the exact solution (indicated by the black line). Thus, while the higher spatial scheme order does not offer a benefit in this scenario with piecewise constant solutions, it still is quite capable to represent the solution. Figure 8b on the right shows a close-up of the solution close to the reflected shock and shows the oscillatory of the high-order approximation. These small oscillations were not observed in [2] with a non-moving wall. Thus, they are induced by the movement of the discontinuity in the masking function. As described in section 3, the wall position is not uniformly well represented, which results in some oscillations as the discontinuity of the masking function moves through the stationary mesh. Yet, the error from this effect is relatively small and does not grow with higher scheme orders as shown in table 2. Table 2 shows the pressure ratio (p 3 /p 2 ) and the distance between the numerical shock and the exact shock location. The pressure ratio in the numerical solution   is obtained by averaging the solution left and right of the shock. Apparently the pressure relation is captured quite well by all scheme orders with this fixed number of 8196 degrees of freedom and we observe an error of less than 0.05% with a decreasing tendency for higher spatial scheme orders. The shock location in the numerical result is assumed to be at the middle point of the jump, and again all simulations provide a good approximation of the reference shock location. Again the highest scheme order yields the smallest error, but the trend is not as clear as in the pressure relation. This can be explained by the approximately same spatial resolution of the wall position in the numerical integration. With the fixed number of degrees of freedom we also have a fixed average density of integration points, which limits the accuracy of the wall position. However, in comparison to the distance between wall and reflected shock (0.25) at this point, the error is again negligibly small and in the largest case of the fourth order approximation less than 0.25% of the distance to the wall.

Formation of shock: Moving Piston
In the previous section 4.2 we investigated the reflection of a shock at a moving wall. Now we want to look into a shock that is formed by a moving wall itself. A piston moves in a tube of fluid at rest. Ahead of it a shock will form and travel away from the moving wall. Behind the piston a suction area will build up with a rarefication traveling away from the moving wall. This one-dimensional setup is well studied in literature [14] and the exact solution for the inviscid equations is known. For this simulation we use a rigid rectangular piston, modeled by the Brinkman penalization and moving with velocity v p = 150 in a one-dimensional domain. The domain length L is set to 1.0 and the width of the piston is 0.04. Initially, the piston is positioned at 0.4. The fluid is initially at rest, with density ρ = 1.0 kg/m 3 and pressure p = 10 5 Pa. A time interval of 0.008s is computed. This setup is particularly significant, as it indicates, whether conservation is maintained by the wall model. A shock will not form, if conservation is not maintained or will have a wrong speed, as a loss of conservation properties occur [14]. We have a closer look into the density and pressure. Important here is the correct capturing of the shock position as well as the pressure and density relation behind and ahead of the formed shock. The exact position of the shock after t = 0.008s is expected at x = 0.819870, whereas the piston interface is located at x = 0.56 at the end of the simulation. During postprocessing we apply a smoothing filter that addresses Gibbs oscillations close to the shock [23]. The state of the simulation result after t = 0.008s is shown in Figure 9 for the density and in Figure 10 for the pressure. Once for the overall domain in Figures  9a and 10b and once with close-ups of the shock in figures 9b and 10b. The overall image shows that the state is well captured by all scheme orders for this number of degrees of freedom. In the close-ups the shown grid indicates the mesh for the fourth order discretization and discontinuities at the element surfaces can be observed. The discretizations with higher spatial scheme orders span accordingly multiple grid cells in figures 9b and 10b. In the run with 32nd order a single element spans 8 of the shown 10 grid cells, all but the outer left and outer right. As with the shock reflection in the previous section we find that the high order can properly handle both the moving geometry and the discontinuity in the solution. Thus, it is viable to use a high spatial scheme order in the complete domain, across discontinuities in state and masking function without special measures.

Subsonic cylinder movement
We now move on to a two-dimensional test case, where we consider the movement of a cylinder through a fluid with a low velocity. As a reference for this simulation, we use the simulation of a cylinder at rest within a fluid moving with a higher velocity. In this setup we neglect again the viscosity and only solve the penalized compressible Euler equations. We define the simulation domain with an extend of [4d × 4d] where d is the diameter of the cylinder. In the reference solution the fixed cylinder is put in the center of the domain. Accordingly, the moving cylinder is put initially shifted from a center such that it reaches the center of the domain after a simulation time of 0.04s. The relative motion between cylinder and fluid has a Mach number of Ma = 0.5. In the reference computation this defines the velocity of the fluid around the cylinder in rest. For the moving cylinder we split the velocity parts and have the cylinder move with half the velocity difference Ma = 0.25 upstream through the fluid that has a velocity of Ma = 0.25. Figure 11 sketches the simulation domain. In the case of the moving cylinder, the cylinder is initially slightly shifted towards the right (dashed line cylinder) and moves to the left during the simulation. The fluid moves from left to right in both computations.  Figure 12 all three mesh resolutions (L) for the moving cylinder are shown. We can observe how the pressure, normalized by the background pressure converges with finer mesh resolution. For the mesh refinement L9 and L10 the normalized pressure at the stagnation point is almost the same. Additionally, in all three cases we can clearly observe the decrease of the pressure value at the stagnation point, as the resolution becomes better. Thus, also the representation of the cylinder.
In Figure 13 we compare the solution in normalized pressure for the mesh refinement level L10 for the moving (M) and non-moving (NM) cylinder (cf. Figure  13a). The pressure behavior for both cases is very close and in agreement with each other. We can observe a slight shift in the stagnation point between both solutions (cf. Figure 13b), which is again due to the movement of the cylinder and the differ-

Supersonic movement
In this section we increase the velocity of the moving cylinder to reach the supersonic regime, where we again see shocks in the flow, but now in a two-dimensional simulation. The simulation domain is rectangular of the form [L ×H] with the length L = 4 · H and H is the height of the domain. The cylinder is 1/4 of the height of the domain and has a diameter of d = 0.25 m Initially the cylinder is located close to the right boundary (outlet) and moves with Mach 1.5 towards the left boundary (inlet). The fluid has a velocity of Mach 1.5 towards the right, resulting in a total relative velocity between fluid and obstacle of Mach 3. Upper and lower boundaries are set to slip walls. This setup is also described in a report by P. Hu et al. [9]. The mesh has elements of the length 1/128 resulting in 65536 elements in total. We solve the inviscid Euler equations with a spatial scheme order of six. Figure 14 shows the density for various positions of the cylinder at consecutive points in time. Ahead of the cylinder the formation of a bow shock can be observed (cf. Figure 14a -Figure 14e). The shock is sharply resolved and drastic change in state can be observed in this supersonic setting. In Figure 14a we can observe the interaction of the bow shock with the walls. The shock is reflected at the upper and lower wall. Besides the interaction of the shock with the wall, we can observe the interaction of the reflected shock with vortices in the wake of the obstacle (cf. Figure  14d and Figure 14e). Further, we can nicely observe the structure of a Richtmyer-Meshkov instability behind the shock along the walls. All figures show a good resolution of the flow features, as small structures are clearly visible. The cylindrical geometry is well modeled and its movement causes the expected physical phenomena. Our results are in good agreement with the report of P. Hu et al. [9]. However, we can observe that small scale structures are more clearly resolved in the presented high-order simulation here. This highlights the advantage of high-order methods even with discontinuities present in the solution. It also shows, that the Brinkman penalization can successfully be used in the supersonic regime.

Rotation of a fan
This three dimensional test case is solved with the full compressible viscous Navier-Stokes equations. The setup features a rotating fan, that is composed by three NACA0012 airfoils with a chord length of 1m. The first blade is initially positioned at an angle of 90 • , the second one at 210 • and the third at 330 • . The pressure is set to 1bar, the density to 1kg/m 3 and the fluid is initially at rest. The fan rotates with a frequency of 6.283Hz. The Reynolds number is with respect to the chord length 67, 114. All lengths are normalized by the cord length. The computational domain is of size [10 × 10 × 2] unit length and the mesh consists of a total of 102400 cubical elements with an edge length of 1/8m. A spatial scheme order of O(6) is used on this mesh. The fan rotates counter-clockwise around its center point at P(0, 0, 1) of the domain. In Figure 15 the time evolution of the rotating fan is shown after 5s, 10s, 15s and 20s from Figure 15a to Figure 15d respectively. The pressure, normalized by the background pressure is shown as a color field, while white streamlines illustrate the velocity pattern in the vicinity of the rotating obstacle. Three strong vortices appear, that are due to the sudden movement of the fan in the initially not moving fluid. Those are still visible in 15d, but are far from the obstacle at that point and the influence from the initial condition vanishes. Clearly the meeting point of high and low pressures at the tip of the blades can be observed. Resulting pressure waves are well resolved and propagated through the domain forming the depicted spiral pattern. The three vortices from the start-up travel towards the outer boundaries, disturbing the outgoing pressure waves (cf. Figure 15b to Figure 15d). Between the fan blades the streamlines illustrate how the flow is captured in circulation areas confined by the blades. The geometry is well represented and the observed behavior is in agreement with the expectations.

Collision of three Spheres
The three-dimensional test case is solved with the compressible Euler equations.
The  illustrates the spheres at their maximum position, with the largest distance to each other after 9.5 s. The spheres move from that position again back to their original location, that is shown in Figure 16b (after 9.8 s). As expected, the fluid is compressed between the spheres, resulting in high velocity value on the right and left of the spheres, as the fluid is squeezed out from the region between them. Finally, the spheres reach their respective position in Figure 16, where their interfaces get in contact with each. It needs to be highlighted, that the numerical scheme is able to handle the collision of objects, without further stability restrictions. Furthermore, the presence of multiple moving geometries can be handled accordingly by the solver and the method used to model the geometries. In Figure 17 the Q-criterion of 7.0 is shown for the time 9.6 s, when the spheres move again back to their origin. Again, high velocity values can be observed in the region between the spheres, where the fluid is compressed. Small structures appear close and away from the spheres, where the fluid is disturbed by the movement of the three geometries.

Conclusion
We applied the Brinkman penalization method in our modal discontinuous Galerkin implementation Ateles and the presented analysis shows that it offers an effective method to model moving obstacles. While the convergence order breaks down at the discontinuity of the masking function for the penalization, the accuracy is still high in the remaining domain and overall not worse than with a lower order discretization. Thus, the complete domain can be computed with the same discretization without the need to adapt to the moving obstacles. The Brinkman penalization scheme worked in our implementation for smooth flow states like a reflected acoustic wave as well as for shocks. It is working in the subsonic and supersonic regimes and independent of the location of obstacles relative to the mesh. Small and closing gaps between moving parts can effortlessly modelled up to the contact of obstacles. We observe small scale oscillations in the solution that are not present in simulations with fixed walls, but they remain small and decrease with better resolutions. An over-integration in the representation of the masking function is necessary and we found that an over-integration by a factor of three is sufficient to mostly recover an accurate geometrical approximation. The presented method enables us to model moving geometries of arbitrary shape within a high-order discontinuous Galerkin scheme and provides a large flexibility in its application.
Further improvements could be offered by incorporating re-projection methods that promise the reconstruction of error convergence even in the vicinity of discontinuities.    Masking function x (red) approximated in polynomial space with 16 terms. Comparison of the numerical integration with over-integration by a factor of three (black) and analytical integration (blue).  Plot for the pressure perturbation of the re ected wave at t = 0.5 for different orders. The numerical reference solution is obtained with a traditional wall boundary condition and a high resolution.   L2 Error in the re ected acoustic pulse with respect to computational effort needed. (a) on the left: the L2 Error for various spatial orders over the required memory in terms of degrees of freedom. (b) on the right: same runs, but 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.     Sketch of domain with cylinder, the non-moving cylinder (solid line) is located at the center and the moving cylinder (dashed line) is initially shifted slightly to the right.      Contour plot showing the Q-criterion of 7:0 after 9:6 s simulation time. The contours are coloured by the velocity magnitude.