Compressible flow simulation with moving geometries using the Brinkman penalization in high-order Discontinuous Galerkin

In 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 fluid flow 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 appropriate 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 artificial 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 briefly discussing the Brinkman penalization technique and its application in our solver and investigate its behavior in fundamental one-dimensional setups, such as shock reflection at a moving wall and the formation of a shock in front of a piston. This is followed by the application to setups with two and three dimensions, illustrating the method in the presence of curved surfaces.


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. Highorder 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 comparatively little data needs to be communicated between elements.
Obstacles in such mesh-based schemes are 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 [22]. 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 [22], Sotiropoulos and Yang [27] and recently Maxey [21]. We use a penalization method as described by Liu and Vasilyev [20].
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 [24]. Kevlahan and Ghidaglia [15] used a pseudo-spectral method and applied this Brinkman penalization method to model stationary, as well as moving geometries.
Liu and Vasilyev [20] 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 [23] and [13]. In [24] the Brinkman penalization is used in the context of finite volumes and finite elements. To the authors knowledge, the Brinkman penalization method used for compressible flows, either involves non-moving obstacles or mostly considered for aero-acoustic problems e.g. in [18] or [20].
In [2] we applied the Brinkman penalization in our high-order discontinuous Galerkin method for compressible Navier-Stokes 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 [26].

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. [18]. 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 by 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 [20] 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. [18], 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 operators and thereby changes the Eigenvalues of the convective system and necessitate smaller timesteps when considered in explicit time integration schemes. The geometrical modelling can therefore primarily be influenced by the permeability terms and the porosity. They allow to tune the modelling towards a solid obstacle. Hence, they need to be chosen appropriately.
As we show in [2], it is feasible to ignore the porosity, introduced by Liu and Vasilyev [20] 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 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: 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]. We showed, that with a smaller value of β the modelling error can be decreased considerably, allowing to neglect the porosity in the conserved quantities, as the permeability enables to model solid obstacles. With the viscous permeability η = β 2 · φ 2 and the thermal viscosity η T = 0.4β · φ, the expected modelling error as described in [20] 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 . a b

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. An algorithm to find this polynomial representation of the masking function for surfaces defined by polygons, for example in the form of STL files is proposed in [17].
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 uses 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 is 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 Fig. 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 [16]. In Fig. 2 we consider an over- integration factor of three. Comparing Fig. 1b without over-integration and the same case with this over-integration by a factor of three in Fig. 2, the improved representation of the wall location becomes apparent. As Fig. 2 shows, an over-integration by factor three already yields a polynomial approximation that is close to the analytical integration. Another approach to improve the geometrical representation is introduced by Engels et al. [7]. The authors assume the masking function to be smooth by introducing a thin smoothing layer. This avoids the strong discontinuity introduced by the jump of the masking function at the interface. However, this approach does not improve the convergence in the vicinity of the geometry surface.

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 one-dimensional setup. We look into the reflection of an acoustic wave from a moving wall. Figure 3 sketches this setup with the initial wave and wall positions.
The traveling wave has the form of a Gaussian pulse given by Eq. (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, non-dimensionalized 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 perturbation, 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 seconds, 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 Fig. 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 (h-refinement). Again a first order convergence can be observed in the overall domain. The moving discontinuity 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 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 Fig. 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 Fig. 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 Sect. Representation of immersed boundaries. 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 re-projection method [10]. Though, such a re-projection 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 Eq. (11) for the pressure relation.
The speed of the reflected shock wave is than given by Eq. (12) [9].
These relations provide us with the piecewise constant exact solution, that we use as a reference to compare the numerical results with.  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 Fig. 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 Sect. Representation of immersed boundaries, 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 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 [19] 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 [19]. 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 post-processing we a b apply a smoothing filter that addresses Gibbs oscillations close to the shock [28]. As in the previous Sect. Reflected shock wave with the reflected shock, we consider discretizations with different scheme orders (O) but with the same overall number of degrees of freedom. We use spatial schemes of order O(4), O (8), O (16) and O(32) with element counts (n) of n2048, n1024, n512 and n256, respectively. Thus, the domain is resolved by 8192 degrees of freedom in each case. The state of the simulation result after t = 0.008s is shown in Fig. 9 for the density and in Fig. 10 for the pressure. Once for the overall domain in Figs. 9a and 10a and once with close-ups of the shock in Figs. 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 Figs. 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

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 extent 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.
In the simulation we use a spatial scheme order of O (16) with three different meshes. The coarsest mesh (L8) uses 64 × 64 = 4096, the intermediate (L9) 128 × 128 = 16384 and the finest (L10) 256×256 = 65536 elements. We compare the solutions by measuring the pressure along the surface of the cylinder with a resolution of one degree. In Fig. 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 Fig. 13 we compare the solution in normalized pressure for the mesh refinement level L10 for the moving (M) and non-moving (NM) cylinder (cf. Fig. 13a). The pressure a b 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. Fig. 13b), which is again due to the movement of the cylinder and the different representation in each time step. However, apart from that, both curves show the same behavior in pressure and are in very good agreement.

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. [12]. 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. Fig. 14a-e). The shock is sharply resolved and drastic change in state can be observed in this supersonic setting. In Fig. 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. Fig. 14d, e). Further, we can nicely observe the structure of a Richtmyer-Meshkov instability behind the shock along the walls. All figures show a good resolution  : a 5s, b 10s, c 15s and d 20s. The normalized pressure is shown in the background. The streamlines (white) illustrate the velocity flow pattern 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 Hu et al. [12]. 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 with respect to the chord length is 67, 114. The kinematic viscosity is defined to be 1.49 · 10 −5 m 2 /s and the thermal conductivity is set to 0.024W /(m · K ). 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 counterclockwise around its center point at P(0, 0, 1) of the domain. In Fig. 15 the time evolution of the rotating fan is shown after 5s, 10s, 15s and 20s from Fig. 15a to d 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. Fig. 15b to d). 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. The simulation was conducted on Germans national supercomputer JUWELS located in Jülich Supercomputing Center using 15k core-h. Each compute node on the system is equipped with 48 cores [14].

Collision of three spheres
The three-dimensional test case is solved with the compressible Euler equations. The simulation domain has a length of L = 4 m, a height of H = 4 m and a width of W = 1.5 m. Three spheres are located in the domain, that are at the beginning of the simulation in contact with each other at their respective interfaces. The first sphere is located at S 1 (−0.1, 0.0, 0.0) and has a radius of 0.2 m, the second and third sphere have a radius of 0.15 m and are located at S 2 (0.0, 0.45, 0.0) and S 3 (0.0, 0.0, −0.45) respectively. All spheres move with a predefined cosine function, away from each other and towards each other, that is defined as A · cos2πt with A being the amplitude, which is 0.1. The Euler equations are solved with a scheme order of O (8). The fluid is in rest, the pressure is defined to be 1 bar and the density is 1 kg/m 3 . The computational mesh has 98304 elements, with an edge length of 1/16 m. The simulation is carried out for 10 s. In z-direction the computational domain is periodic. In Figure 16 different positions of the spheres at different points in time are shown. Figure 16a 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 Fig. 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 Fig. 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 Fig. 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. Please note that simulation results provided by the solver are represented by polynomial series. For visualization purposes this representation is voxelized during post-processing. Therefore, the spheres in Fig. 16a appear relatively rough, as the voxelization was done in a coarse resolution for the geometry here, while the sliced data of the flow state is resolved to a finer level.

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 overintegration 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 e.g. in [8].
discussion to the inviscid part. The vicsous terms are discretized using the DG method with symmetric interior penalty [11]. We refer to the equations obtained by neglecting diffusive terms from the compressible Navier-Stokes equations as the Euler equations. They provide a model for the conservation of mass, momentum and energy of the fluid and are defined in vectorial notation by ∂ t u + ∇ · F(u) = 0, with respective initial and boundary conditions. u is a vector, that includes all variables to be conserved. Further the flux function F(u) = (f(u), g(u)) T is defined for the spatial dimensions and is in 2D given by with density ρ, velocity v = (u, v) T , specific total energy E, and pressure p. The system is closed by the equation of state. Which for an ideal gas provides p = (γ − 1)ρ e − 1 2 (u 2 + v 2 ) . Here γ = c p c v is the ratio of specific heat capacities and e is the total internal energy per unit mass.
The discontinuous Galerkin formulation of the mentioned equation is obtained, by multiplication with a test function ψ and integration over the domain Ω. Afterward, we apply integration by parts and obtain the weak formulation Where ds denotes a segment on the surface integral. The overall domain Ω is split into smaller, n non-overlapping elements Ω i |i = 1, 2, . . . , n, such that Ω = ∪ n i=1 Ω i and Ω i ∩ Ω j = ∅∀i = j. Within in each of these elements a polynomial series P m with a maximal degree of m ≥ 0 is used to approximate the solution u. We refer to this approximate solution as u h (x, t), in each element i.
The expansion coefficientsû i are the degrees of freedom of the numerical solution. It is important to keep in mind, that no global continuity requirement for u h in the previous definition is needed. Dividing the integrals in Eq. (14) into a sum of integrals over elements Ω i , we can obtain the space-discrete variational formulation with test functions ψ j in each By using m + 1 linearly independent test functions ψ 0 , . . . , ψ m we arrive at a fully determined linear system of m + 1 equations for the m + 1 degrees of freedom. The appearing inner products of φ i and ψ j , can be precomputed and we can represent the system in the form of matrices that are multiplied with the degrees of freedomû representing the state. Since the numerical representation is only supported locally, the flux term is not uniquely defined at the element interfaces. Thus, the flux function has to be replaced by an approximative numerical flux function F * (u − h , u + h , n), where u − h and u + h are the interior and exterior traces at the element face in the normal direction n to the interface.
Hereby M,S are the mass and the stiffness matrices and M F are the so called face mass lumping matrices. The computational cost to perform these operations obviously depends on the choice of basis functions. Using the orthogonality of Legendre polynomials the mass matrix gets diagonalized and trivially invertable, such that the multiplication with M −1 can be performed in O(m + 1) operations. The stiffness matrix is not as trivially applied, however due to the recursive properties of the Legendre polynomials it can still be computed in O(m + 1) operations by exploiting this recursiveness. The obtained ordinary differential equations (17) can be solved in time using any standard time-stepping method, e.g. a Runge-Kutta method [2]. In our DG discretization, we usually employ the HLL flux as numerical flux.