Galilean-Invariant Characteristic-Based Volume Penalization Method for Supersonic Flows with Moving Boundaries

This work extends the characteristic-based volume penalization method, originally developed and demonstrated for compressible subsonic viscous flows in (J. Comput. Phys. 262, 2014), to a hyperbolic system of partial differential equations involving complex domains with moving boundaries. The proposed methodology is shown to be Galilean-invariant and can be used to impose either homogeneous or inhomogeneous Dirichlet, Neumann, and Robin type boundary conditions on immersed boundaries. Both integrated and non-integrated variables can be treated in a systematic manner that parallels the prescription of exact boundary conditions with the approximation error rigorously controlled through an a priori penalization parameter. The proposed approach is well suited for use with adaptive mesh refinement, which allows adequate resolution of the geometry without over-resolving flow structures and minimizing the number of grid points inside the solid obstacle. The extended Galilean-invariant characteristic-based volume penalization method, while being generally applicable to both compressible Navier–Stokes and Euler equations across all speed regimes, is demonstrated for a number of supersonic benchmark flows around both stationary and moving obstacles of arbitrary shape.


Introduction
Numerical simulation of complex geometry flows in a computationally efficient manner represents a challenging problem, especially in the presence of moving/deformable boundaries. Solid bodies are introduced by imposing appropriate boundary conditions upon surfaces, and to that end, several approaches have been developed. These methods can be separated into two major groups: body-fitted mesh [1][2][3] and immersed boundary (IB) methods [4][5][6].
For conventional structured/unstructured body-fitted grid methods, the numerical meshes are conformal to the complex boundaries. Therefore, it is straightforward to impose exact boundary conditions and to attain satisfactory accuracy by employing fine meshes for boundary layers, where high resolutions are required, which is critical for high Reynolds number flows. However, there are some disadvantages to these methods. The grid must be carefully constructed to precisely fit an obstacle. In most cases, this precludes the use of structured Cartesian grids. The process of mesh generation is highly dependent upon the obstacle geometry and can become computationally expensive, especially for complex surfaces. This issue is compounded for moving or deforming obstacles, which require continuous adaptation or re-meshing throughout computation of the solution [6]. The grid generation process may be very expensive: it is not an easy task to generate a good-quality grid, as even simple geometries and simulations for moving boundary problems become prohibitively expensive due to grid generation and solution interpolation to the new mesh at each time step.
The use of IB methods avoids the cost and complications of body meshing by directly introducing the effects of obstacles upon the governing equations. Solid body effects, thus embedded within the flow itself, obviate the rigors of positioning nodes upon a surface. In the IB approach, the boundary conditions on the obstacle surface are indirectly imposed by introducing additional terms in either the discretized [7] or continuous [8] governing equations. In fact, depending on the way the boundary conditions are imposed, IB methods can be roughly grouped into two main classes: discrete [4,7,9,10] and continuous [8,11,12] forcing approaches. Discrete forcing methods are based on modifications of discretized equations, which makes them difficult to generalize due to the direct dependence on methods of discretization. Continuous forcing methods, also referred to as volume penalization (VP) methods, are based on solving modified penalized governing equations, where penalty terms in the form of additional forces (or different feedback mechanisms) ensure the approximation of the boundary conditions. Despite their seeming simplicity and wide use, discrete forcing methods lack generality and flexibility across different solvers and the ability to rigorously control the accuracy of the approximated boundary conditions [6], whereas continuous formulations are independent of the discretization methods and usually have an ability to rigorously control the error of the solution through a penalization parameter [12][13][14]. Since Peskin's IB method [15] was originally introduced to study flow patterns around heart valves, a number of other IB techniques have been developed for incompressible viscous flows around complex solid boundaries. In contrast to Peskin's method using external forces to simulate the immersed boundaries, Cartesian grid methods [16][17][18][19] and ghost-cell IB methods [5] directly impose the boundary conditions on the immersed boundaries.
Starting with the work of Arquis & Caltagirone [8] that introduced a VP formulation for incompressible flows around solid obstacles, which is commonly known as the Brinkman penalization (BP) method, considerable efforts have been put into the development of VP methodologies for incompressible flows [12,14,[20][21][22]. The main idea of these methods is to model arbitrarily complex solid obstacles as porous media with permeability approaching zero. A principal strength of Brinkman-type penalization, compared to discrete forcing IB methods, is that the error between the solutions of penalized and non-penalized equations can be estimated rigorously in terms of the penalization parameter with the approximate solution converging to the exact one in a predictable fashion [12][13][14]20,23]. Because this VP is simple and cheap to calculate, it is well-suited for flows in complex geometries, including moving obstacles. Much work has been undertaken to refine BP methods for various numerical techniques, including pseudospectral [24][25][26][27], finite-element/finite-volume [28], and wavelet-based [23,[29][30][31][32] formulations.
Until recently, most of the effort was put into the development of IB methods for incompressible viscous flows. One of the first attempts at developing a discrete forcing IB method for compressible flows was undertaken in Reference [33], where the flow around a circular cylinder and an airfoil at high Reynolds numbers was modeled. However, the formulation of Ghias et al. [33] did not take into account acoustic wave reflection and transmission at the interface between fluid and solid media, which are critical in some applications with acoustic and shock wave propagation around solid obstacles. Similar problems were observed in the extended impedance mismatch method [34], which extended the original approach of Chung [35] to unsteady non-uniform flow problems.
The Brinkman VP approach was generalized to compressible flows in [36], where in addition to the penalization of momentum and energy equations, the continuity equation was also modified inside the obstacle to be consistent with the porous media flow physics. In the extended compressible BP model, the penalized porous region acts as a high impedance medium, resulting in negligible wave transmissions. As in the incompressible case, the error bounds of the compressible BP method could be rigorously estimated in terms of porosity and permeability parameters. The extended BP method was suc-cessfully applied for simulating subsonic compressible flows in both viscous [36] and inviscid [37] conditions. A number of discrete forcing IB methods have been recently extended to compressible flows as well. The suitability of Cartesian grid methods for modeling compressible high-Reynolds number flows using delayed detached eddy simulation was studied in Reference [38], while sharp-interface [39] and ghost-cell [40,41] IB methods were extended to compressible flows. Moreover, a comparative study of compressible BP and discrete forcing Cartesian grid IB methods was performed in [42] for compressible viscous flows, where an excellent agreement between the two different approaches was observed. It should be noted that, for the compressible BP method, the correct shock reflection was obtained using theoretically predicted values of total energy after shock reflection, as a target in the Brinkman term for the energy equation, which is not known in the general case, and thus should not be used. The imposition of the temperature conditions on the wall results in the incorrect reflection of shock waves, due to creation of a large pressure gradient inside the penalized region and consequent graduate pressure decay due to seepage. This necessitated the development of a more general VP approach capable of approximating compressible flows across all speed regimes, as well as imposing general homogeneous/inhomogeneous Neumann and Robin type boundary conditions. The first step to this end was to develop a characteristic-based volume penalization (CBVP) method [43], where the BP approach was extended to general boundary conditions by introducing hyperbolic VP terms with characteristics pointing inward on solid obstacles. As with the original formulation, CBVP maintains rigorous control of the error, through a priori chosen penalization parameters, for all types of boundary conditions. The CBVP method was demonstrated for scalar diffusion equations and compressible Navier-Stokes equations for a variety of boundary conditions for both stationary and moving obstacles and was extended to Euler equations [44,45] for stationary obstacles. Lavoie et al. [45] suggested using modified boundary conditions to impose conservation of entropy and total enthalpy in the wall-normal direction instead of adiabatic boundary conditions used in [44], which lead to more accurate results on coarse meshes. However, even though the results of simulation of compressible flows using the compressible BP [36] and CBVP [43] methods were demonstrated to be accurate, as pointed out in Reference [46], mathematically, both formulations were not actually Galilean-invariant. Indeed, Komatsu et al. [46] suggested a Galilean invariant extension of the compressible BP method [36].
In this work, we propose a novel Galilean-invariant extension of the CBVP method to systems of hyperbolic partial differential equations (PDEs) involving moving geometries. This extended formulation, which is referred to as Galilean-invariant characteristic-based volume penalization (GI-CBVP), is demonstrated for supersonic complex geometry flows. The penalized Euler and Navier-Stokes equations are solved for flows around stationary and moving obstacles. It should be noted that the new method is well suited for use with adaptive mesh refinement (AMR) techniques [47]. In fact, the VP approach does not employ body-conformal meshing, and high resolution is required around surfaces for computational accuracy and proper definition of geometry. The use of AMR grids maintains the resolution of solid geometry without over-resolving flow structures. Additionally, the number of nonphysical points lying inside the obstacle can be minimized to those necessary to support the boundary conditions, which is particularly important for obstacles inhabiting a large portion of the computational domain.
All of the results reported in this paper were obtained using the adaptive wavelet collocation (AWC) methodology [48][49][50][51]. AWC represents a general numerical approach that utilizes a wavelet decomposition to dynamically adapt the local mesh resolution on steep gradients in the solution, while retaining a predetermined order of accuracy. Geometry definitions are treated as any other flow variable by AWC methods, where the local grid adapts efficiently and dynamically in order to resolve surfaces, even for moving obstacles. The AWC method was applied to shockless supersonic flows in [52], while a wavelet-based shock capturing scheme was developed to handle flow discontinuities in [53]. The combined VP/AWC approach was undertaken for incompressible flows past stationary obstacles in [54,55]. The interested reader is referred to [56,57] for a complete review of this methodology and its applications.
The rest of the paper is organized as follows. Theoretical and numerical aspects of the new GI-CBVP method are provided in Section 2, while the results of simulations of different benchmark problems are presented and discussed in Section 3. Finally, some conclusions are drawn in Section 4.

Mathematical Formulation
In the framework of the CBVP approach, the proper boundary conditions are imposed by modifying the governing equations inside the solid region. Consider a compressible flow problem, defined in a physical domain R containing a solid obstacle Ω, and governed (outside of Ω) by a generalized evolution equation, say where RHS is simply representing the physical terms on the right hand side of the equation. Note that Equation (1) can be hyperbolic or parabolic in nature. A masking function χ(x, t) is defined by with Ω standing for the physical domain instantaneously occupied by the solid obstacle. This masking function separates the domain into a physical region and a penalized region associated with the solid obstacle. The general boundary conditions, written in operator form as are imposed in a similar fashion as BP, namely For Dirichlet type boundary conditions, Equation (4) reduces to BP [12], while homogeneous Neumann type boundary conditions result in the following penalized equation: , with n i representing the components of the inward-oriented surface normal vector. This procedure can be generalized to Robin type boundary conditions, as was demonstrated in Reference [43].
It is important to emphasize that despite the similarity of the different formulations in operator form for Dirichlet, Neumann, and Robin type boundary conditions, the mechanisms underlying each condition are very different. For Dirichlet type boundary conditions, the penalization term acts as a simple feedback mechanism that forces the solution inside the obstacle to the desired target value. For Neumann/Robin type boundary conditions, the solution is convected inside the obstacle along the inward pointing characteristic nor-mal to the surface. In any case, the penalization parameter η determines the timescale of the process. As it was demonstrated in [43], the solutions of the penalized equations for either Dirichlet or Neumann/Robin type penalization converge to the solutions obtained by imposing the exact boundary conditions. Since the penalization timescale is controlled through the parameter η, selecting η 1 causes Equation (4) to become quasi-steady within Ω, on the normalized problem timescale, therefore imposing the intended boundary conditions on the interface. For vanishing η, the increased disparity in timescales asymptotically controls the penalization error. However, reducing the error increases the computational complexity. Since 1/η is the characteristic velocity for Neumann/Robin type boundary conditions, a reduction of η is also accompanied by increased stiffness, which represents a well-known problem with BP that is mitigated through stiffly-stable solvers [36].
Note that the CBVP method relies on the availability of the normal vector in the entire obstacle interior, that is, where χ = 1. This normal vector can be either directly prescribed, for simple objects, or constructed through the use of distance or level set functions, defined analytically or numerically, for more complex geometries (e.g., [58]).

Volume Penalization of Euler Equations
For numerical simulation of the compressible Euler equations, the following boundary conditions at the fluid-obstacle interface are either explicitly assumed or implicitly built into numerical methods (e.g., by constructing ghost cells while using symmetry assumptions) [59]: where the superscripts (·) n and (·) τ denote normal and tangential components, respectively. The vector fieldû = u − U o represents the relative to obstacle fluid velocity, with û τ 2 = u τ jû τ j , while κ τ is the streamline surface curvature moving with the obstacle frame of reference. Despite the fact that the only mathematically rigorous boundary condition for Euler equations would be the no-penetration condition, the above boundary conditions are commonly used because they can be derived as outer asymptotic conditions for viscous boundary layer flows. In the present context, these boundary conditions are enforced through the VP approach, as is discussed in the following. For the sake of simplicity, hereafter, the penalized equations are written for the interior of the obstacle, where χ = 1. Thus, the boundary conditions (6) are enforced as follows: Note that in these equations and thereafter, two different parameters η b and η c are used to highlight different asymptotic convergence mechanisms for Brinkman and characteristic-based penalizations, respectively. Also note that, since Equation (7) is defined throughout the obstacle interior, both normal and tangential velocities u n and u τ , as well as the streamline curvature, need to be defined inside the obstacle. This can be achieved by making use of the normal vector field n, so that the different velocity fields are determined as: u n = n i u i n; u τ = u − u n ;û n = n iûi n;û τ =û −û n . As to the streamline curvature, it holds κ τ = −τ i τ j ∂n i ∂x j , where τ =û τ û τ represents the unit tangential vector. Using this definition, the pressure equation takes the following simpler form Therefore, rewriting Equation (7) for conserved variables, using Equation (8) and the total energy definition the following evolution equations in the interior of the obstacle are obtained: where ∇ 2 (·) stands for the Laplacian differential operator. In these equations, the diffusive and smoothing terms, in addition to convective, Brinkman, and non-conservative forcing terms, are considered in order to stabilize the numerical implementation. Note that the convective terms practically correspond to the CBVP approach. In particular, the Brinkman term in the normal momentum equation is complemented by the diffusive term, which is added to control spatial resolution of the inner boundary layer for a given mesh size ∆, and penalization parameter η b , and to avoid spurious oscillations at the boundary. The corresponding numerical viscosity ν n = α 2 ∆ 2 /η b (with α = O(1)) is chosen to guarantee that the diffusion in the interior of the obstacle is occurring on the same timescale as η b . For more details on the diffusion term determination, the reader is referred to Reference [43].
It should be emphasized that the application of convective terms in the entire interior domain of the obstacle is undesirable for two reasons: (i) normal vectors are ill-defined far from boundaries, due to possible intersection of the characteristics and existence of extrema points in the distance function, (ii) an additional computational cost and deterioration of stability/accuracy of characteristic equations far from the boundaries. For these reasons, the characteristic equations are only solved in a narrow internal layer close to the obstacle boundary, denoted by masking function χ h , which is just wide enough to propagate the solution along characteristics within the actual stencil of the differential operator. As a result, all convective and non-conservative terms in Equation (10) are only present in this layer. Differently, the BP and corresponding penalization terms are applied throughout the entire region defined by obstacle masking function χ. Since the characteristics are pointed inwards, the solution of Euler equations in the fluid region and in the thin layer defined by masking function χ h are not affected by the solution in the remaining interior region of the obstacle, denoted by masking function χ d . Note that the masking functions χ h and χ d are orthogonal, i.e., χ = χ h + χ d . The convective terms are complemented by smoothing terms using the same form as diffusive terms but defined only in the region of the masking function χ d . The schematic illustration of different masking functions and regions where they are defined is given in Figure 1. Finally, when the Euler equations are solved together with Equation (10) for the penalized region, the penalization terms for evolution of normal and tangential momentums are combined according to resulting in the simultaneous enforcement of boundary conditions at the surface for both normal and tangential velocity components.

Volume Penalization of Navier-Stokes Equations
The penalization approach introduced in the previous section is also applicable to the Navier-Stokes equations for problems with adiabatic walls. Since the form of no-slip boundary conditions is mathematically equivalent to the form of no-penetration boundary conditions, Equation (10) is directly applicable for the penalization of the Navier-Stokes equations after performing the following substitutions: u n → u, u τ → 0, and U n o → U o . In other words, the penalization equation for the tangential velocity component is simply ignored. For the solution of the Navier-Stokes equations with non-zero heat flux, the adiabatic wall boundary conditions (6) need to be replaced by where q is the normalized heat flux. The latter variable can be either prescribed or expressed as a function of the wall temperature. The penalization equation for the temperature in the system (7) can be modified following the procedure given by Equation (4).

Moving Obstacles
When solving flows with moving boundaries, it is important that the underlying problem (governing equations and associated boundary conditions) is invariant with respect to a change of the reference frame, which is commonly referred to as Galilean invariance. In this case, the solution of the fluid flow equations should be the same in every inertial frame. Unfortunately, while boundary conditions (6) are Galilean-invariant, the system of penalized Equation (10) is not. In order to become Galilean-invariant, these equations need to be reformulated in the reference frame moving with the obstacle velocity and then transformed for the stationary frame. This procedure results in adding Lagrangian terms to the equations, according to the following transformation: . To theoretically demonstrate the Galilean invariance of the proposed formulation let us consider the Galilean-invariant formulation of Equation (4) (13) where the masking function for the moving obstacle is explicitly written as Assuming that the RHS is Galilean-invariant outside of the obstacle, the formulation (13) is Galilean-invariant for general boundary condition, including the Dirichlet boundary condition on the velocity, which follows from the Galilean invariance of spatial operators, Dirichlet boundary conditions of non-velocity variables, and the relative velocityû = u − U o used in either no-slip or no-penetration boundary conditions. Thus, the proposed extension of a characteristic-based volume penalization method is Galilean-invariant on the formulation level.
This way, the system of penalized Equation (7) in the interior of the obstacle is rewritten as Analogously, the system of penalized evolution Equation (10) in the interior of the obstacle becomes where the Lagrangian terms are explicitly multiplied by masking function χ to emphasize that they are applied throughout the entire obstacle. It is worth noting that the original formulation of the CBVP method for moving obstacles, in the context of compressible Navier-Stokes equations [43], did not involve any Lagrangian terms, without affecting the accuracy of the results. The reason for this weak sensitivity to the presence of these additional terms was the use of no-slip boundary conditions and smooth variation of the solution in the vicinity of the solid boundary. In fact, the correct boundary conditions were simply achieved through the BP term, which adjusted the solution on the timescale of η b . The situation is drastically different for the Euler equations, since the changes in the tangential velocity field can be substantial if the solution inside the obstacle is not convected with the obstacle velocity. The results of the numerical experiments, which are reported in Section 3.4, confirm higher accuracy and better convergence of the new GI-CBVP formulation (16). The use of Lagrangian terms in the penalized equations is, however, recommended for both Euler and Navier-Stokes approaches.
Finally, due to the local nature of the penalized equations, the proposed formulation is also applicable for rotating and deformable objects, in addition to obstacles with translational motions. In the former case, the obstacle velocity at every interior point needs to be provided, either analytically or through the solution of a governing equation for the evolution of the obstacle motion or its deformation.

Adaptive Wavelet Collocation Method
The AWC method utilizes the wavelet-based decomposition of the flow field unknowns to dynamically adapt the local mesh resolution on steep gradients in the solution, while retaining a predetermined order of accuracy [29,48,49]. Formally, a scalar spatial field f (x) can be represented in terms of wavelet basis functions as where φ 0 l and ψ µ,j k are three-dimensional scaling functions and wavelets of different families (µ) and levels of resolution (j), respectively. The above wavelet decomposition can be thought as of a multi-resolution representation of f , where each level of resolution consists of a family of wavelets ψ µ,j k having the same scale but located at different positions. In order to compress the numerical solution, wavelet filtering is performed through wavelet coefficient thresholding. Given the number of resolution levels, say j max , the wavelet filtered variable is defined by where stands for the non-dimensional threshold to be prescribed. A good approximation of the unfiltered field can be retained even after discarding a large number of wavelets with small coefficients, because the coefficients d µ,j k are small unless f causes a significant variation on the level of resolution j, in the immediate vicinity of the wavelet ψ µ,j k location. The reconstruction error was demonstrated to converge as with C = O(1). Due to the one-to-one correspondence between wavelets and collocation points, the nodes are omitted from the computational mesh if the associated wavelets are excluded from the truncated approximation (18). Thus, the multilevel structure of the wavelet approximation provides a natural way to obtain the solution on a nearly optimal numerical mesh, which is dynamically adapted to the evolution of the main flow structures, both in location and scale. The multi-resolution wavelet decomposition (18) is used for both mesh adaptation and interpolation, while the derivatives at the adaptive computational nodes are found by differentiation of the local wavelet interpolant at the appropriate level of resolution. In addition, for supersonic flows, the AWC method is supplemented with a shock capturing scheme [53], which uses wavelet coefficients to define local numerical viscosity in the neighborhood of the shock waves and to automatically switch it off away from them.

Numerical Results and Discussion
The original CBVP method was already demonstrated for a number of problems with smooth solutions, such as heat equation, acoustic wave reflection, and low Mach compressible viscous flow around a two-dimensional cylinder [43]. The present study focuses on evaluating the novel aspects of the extended method that is developed for solving systems of hyperbolic PDEs with discontinuous solutions, in the presence of stationary or moving obstacles of arbitrary complexity. For that purpose, a number of different benchmark tests are performed. Specifically, the first test is represented by the one-dimensional shock tube problem, where the shock reflection from the wall is quantitively evaluated, with the wall being approximated by the GI-CBVP method. The second and third benchmark problems test the ability of the new method to correctly model the effect of the obstacle for the generation of oblique and detached shock waves, while the fourth test demonstrates the Galilean invariance of the GI-CBVP approach. In addition, the two-dimensional benchmark problems illustrate the efficiency of using the adaptive wavelet-based mesh refinement with the GI-CBVP method and the ability of the integrated approach to adequately resolve the obstacle geometry without over-resolving flow structures and minimizing the number of grid points inside the solid obstacles.

Benchmark I: Normal Shock Wave Reflection
Let us consider the one-dimensional Sod's (shock tube) test problem in the spatial domain R = [−0.5, 1.5], where the wall is prescribed as a volumetric obstacle occupying the region Ω = [1, 1.5], and the shock is initially located at x 0 = 0. The initial states in front and behind the shock are chosen according to the Rankine-Hugoniot jump conditions. Practically, after prescribing the post-shock conditions ρ L = 1, u L = 2 and p L = 1, as well as u R = 0, the pre-shock conditions ρ R and p R , as well as the shock velocity S, are found according to the given jump conditions. The schematic diagram of this problem is reported in Figure 2, where x S represents the shock location and L stands for the length of the computational domain. Different calculations are conducted for a series of uniformly spaced, non-adaptive computational meshes consisting of a number of grid points N between 2 9 and 2 13 . Pressure seepage and amplitude error are significantly reduced when compared to the results obtained with the compressible BP method [36], but a phase lag is still present. In Figure 3, local pressure plots at the end of the simulation, when the reflected shock front is far from the wall, are presented for coarse and fine mesh resolutions, with the penalization coefficients being set to the values of η c = 10 −3 and η b = 10 −5 .
(a) (b) Figure 3. Pressure profiles along the shock tube for (a) coarse (N = 2 9 ) and (b) fine (N = 2 11 ) grid resolutions. The penalized solution for the reflected shock is compared with that obtained using exact algebraic boundary conditions. The theoretical shock location is also indicated.
The penalty parameters determine the transition time of the shock wave reflection from the wall, while the spatial discretization size ∆x determines the shock thickness, with all three parameters contributing to the phase lag. In fact, the characteristic timescale of the finite width shock wave reflection from the solid wall is τ swr = S/∆x, while the timescale of the transitional effects due to the penalization procedure is τ η = max(η b , η c ). Taking into account the different rates of convergence for Brinkman and CBVP methods that are, respectively, O(η 1/2 b ) and O(η c ) [43], the BP parameter is constrained as η b = η 2 c , for a given value of η c . Results of parametric sweep for η c ∈ 10 −3 , 10 −1 , corresponding to four different grid resolutions, are provided in Figure 4. One can see that, as long as the penalization transitional timescale is shorter than the shock reflection timescale, τ η < τ swr , the numerical discretization error associated with finite width shock readjustment is dominating. For the opposite case, when τ η > τ swr , the error converges linearly as O(η c ).

Benchmark II: Oblique Shock Wave
In order to demonstrate the ability of the proposed methodology to correctly capture both oblique and detached shock wave reflections, supersonic inviscid flow past wedges with subcritical and supercritical wedge angles are simulated. The presence of the wedge is modeled, employing the GI-CBVP methodology. The incident shock wave conditions are chosen to be the same as in the previous benchmark problem, namely, u L = 1.6, p L = 1, ρ L = 1, u R = 0 and ρ R , p R and S corresponding to the Rankine-Hugoniot jump conditions. The two-dimensional computational domain is R = [−0.5, 1.5] × [0, 4] and the wedge apex is located at the point (0.28, 0). The AWC method [48], supplied with the shock capturing scheme [53], with the effective mesh resolution of 513 × 1025 points and eight levels of resolution, is used for this test. Inflow boundary conditions are applied on the left boundary with values equal to the left shock state values. On the top boundary, and on the portion of the bottom boundary located before the wedge, symmetry boundary conditions are used. Outflow boundary conditions are imposed on the right boundary, outside of the wedge. Moreover, tangential diffusion boundary conditions are used on the wedge domain boundaries.
For subcritical flow, the angle of the reflected shock is given analytically by the following relation where θ is the wedge angle, β is the deflected shock angle, and Ma is the Mach number of the incident shock wave. The inviscid supersonic flow past a wedge with the deflection angle of θ = 10 • is simulated. The solution of the problem is shown in Figure 5a with the object geometry displayed in black. The adaptive computation mesh, demonstrating the adaptive resolution of solid surfaces and shock structures, is demonstrated in Figure 5b. The attached oblique shock wave with the oblique angle of β ≈ 51 • is illustrated in Figure 5a,b by the yellow line, with the theoretical value of β = 51.1 • being predicted by Equation (20) with Ma = 1.6. The adaptive computation mesh, demonstrating the adaptive resolution of solid surfaces and shock structures is shown in Figure 5b. It is worth noting the presence of the fine local mesh inside of the obstacle at the expansion corner, which is caused by ill-defined normal vectors.

Benchmark III: Two-Dimensional Supersonic Flows around Blunt Bodies
In order to demonstrate the ability of the proposed GI-CBVP method to approximate the supersonic flows around stationary blunt bodies and the efficiency of the adaptive wavelet-based mesh refinement to adequately resolve the geometry, while minimizing the number of grid points inside the solid obstacle, unsteady simulations of two-dimensional inviscid flows around multiple circular cylinders are performed. The inflow parameters are the same as in Section 3. The AWC method [48] supplied with a shock capturing scheme [53], with the effective mesh resolution of 1025 × 1025 and eight levels of resolution, is used for the problem. The penalized Euler equations are solved with the obstacle geometry approximated by the GI-CBVP methodology (10).
Results of the unsteady AWC simulation with the GI-CBVP method of the flow around an array of randomly placed cylinders are presented in Figure 6, where the instantaneous solution is examined in terms of the numerical schlieren image and corresponding adaptive computational mesh colored by the pressure field. These results illustrate the ability of the proposed method to effectively represent more complicated flow geometry and highlight the advantages of using the GI-CBVP method in combination with wavelet-based AMR. The use of the AWC methodology allows for the efficient resolution of solid surfaces and shock structures, while minimizing the computational costs of the simulation.

Benchmark IV: Galilean Invariance
In order to illustrate the Lagrangian modification introduced for moving obstacles in Section 2.4 and demonstrate the conformity with Galilean invariance, two additional simulations using the modified governing Equation (16) Figure 7, the initial setup and snapshot of the density and pressure fields at both the beginning and the end of the two simulations are provided. These simulations were run until the second cylinder reached the same location as the first one. Apparently, resolved fields and underlying adaptive grids are the same, with only minuscule discrepancies in the wake region, associated with different instability triggering for stationary and moving cylinders. The proven computational Galilean invariance is important, since mathematical (in terms of governing equations) and physical (in terms of inertial frames of reference) equivalence of the simulations were not guaranteed numerically by the original CBVP method [43] without proper modifications. In fact, the results of the numerical simulation of the moving cylinder without the Lagrangian correction term (not presented here) demonstrated a significant loss of accuracy as well as slow convergence of boundary conditions to target values.

Conclusions
This work originally extends the characteristic-based volume penalization method to a hyperbolic system of PDEs that governs compressible flows with moving complex geometries. The proposed methodology results in being general, Galilean-invariant, and allowing the imposition of homogeneous and inhomogeneous Dirichlet, Neumann and Robin type boundary conditions. Both stationary and moving geometries can be treated, where the approximation error is defined a priori, being controlled by the penalization parameters. The method permits the generic formulation of boundary conditions for both integrated and non-integrated variables in a systematic manner that parallels the prescription of exact boundary conditions. The novel approach is demonstrated to be applicable to both Euler and Navier-Stokes equations across all speed regimes. The proper approximation of boundary conditions for moving geometries is achieved by means of the inclusion of Lagrangian terms into the penalized equations, which makes the formulation Galilean-invariant and increases both stability and accuracy of the methodology, especially in the context of Euler equations, where tangential velocity is not continuous across the fluid-solid interface. Moreover, the proposed methodology is demonstrated to be well suited for use in conjunction with wavelet-based adaptive mesh refinement methods. This is particularly notable in the case of moving obstacles, where transients can be exploited to optimize the computational cost of the simulations.

Conflicts of Interest:
The authors declare no conflict of interest.