Numerical simulation of a complex moving rigid body under the impingement of a shock wave in 3D

In this paper, we take a numerical simulation of a complex moving rigid body under the impingement of a shock wave in three-dimensional space. Both compressible inviscid fluid and viscous fluid are considered with suitable boundary conditions. We develop a high order numerical boundary treatment for the complex moving geometries based on finite difference methods on fixed Cartesian meshes. The method is an extension of the inverse Lax-Wendroff (ILW) procedure in our works (Cheng et al., Appl Math Mech (Engl Ed) 42: 841–854, 2021; Liu et al.) for 2D problems. Different from the 2D case, the local coordinate rotation in 3D required in the ILW procedure is not unique. We give a theoretical analysis to show that the boundary treatment is independent of the choice of the rotation, ensuring the method is feasible and valid. Both translation and rotation of the body are taken into account in this paper. In particular, we reformulate the material derivative for inviscid fluid on the moving boundary with no-penetration condition, which plays a key role in the proposed algorithm. Numerical simulations on the cylinder and sphere are given, demonstrating the good performance of our numerical boundary treatments.

where, x , y and z are the outer boundaries of the computational area along x-, y-and z-direction, parallel to the yz-, xz-and yz-plane, respectively. Inflow, outflow and reflection boundary conditions are imposed at the left boundary, the right boundary of x and the two boundaries of y , respectively. Notice that, although the physical problem we actually consider is infinite in the z-direction, the domain is truncated in numerical computation. Due to the fact that the rigid body moves in a confined space during a limited period, the area can be assumed to be large enough in the z-direction that the body and the reflective shock waves would not touch z , and then we can impose the outflow boundary conditions on z . Hence, the boundary treatment on z will not affect the computational results of the internal flow field near the rigid boundary. The boundary conditions on the boundary (t) vary with the properties of the fluid. In this paper, we consider both the compressible inviscid fluid and viscous fluid with suitable boundary conditions.
The governing equations for compressible inviscid fluid and viscous fluid in the threedimensions are written in the same form as follows: where (x, y, z) ∈ (t), W = (ρ, ρu, ρv, ρw, E) T and ( Here, ρ, u, v, w, p and E stand for density, velocity in x-, y-and z-directions, pressure and total energy per volume, respectively. On the right hand side of (1), with the components of the viscous stress tensor given by and σ 1 = u τ 11 + v τ 12 + w τ 13 + (c 2 ) x (γ − 1)Pr , σ 2 = u τ 21 + v τ 22 + w τ 23 + (c 2 ) y (γ − 1)Pr , Here, Pr is the Prandtl number, Re is the Reynolds number, T = p/ρ is the temperature, and c = √ γ p/ρ is the sound speed. An equation of state relating the pressure with other variables is given as γ is the adiabatic constant, which equals to 1.4 for an ideal polytropic gas.
• Compressible inviscid fluid. In this case, δ = 0 in the governing Eq. (1) and it becomes the three-dimensional Euler equation, The no-penetration boundary conditions are imposed on the boundary of the body surface (t): with u = (u, v, w) T . V b and n are the velocity and the outward unit normal vector on (t), respectively. • Compressible viscous fluid. In this case, δ = 1 in the governing Eq. (1) and it is the three-dimensional Navier-Stokes (NS) equation, On (t), we consider the isothermal no-slip wall boundary condition: where V b = (u b , v b , w b ) T and T b are the velocity and the temperature of the rigid body, respectively.
In this paper, we are trying to numerically solve the partial differential equations (PDEs) (6) and (8) with high order finite difference methods on the fixed Cartesian mesh. Despite the easy generation of the mesh, for such problems there are two main difficulties in numerical simulation: 1. A high order finite difference scheme often has a wide stencil, thus we have to evaluate the values at the ghost points near the boundary. 2. Another difficulty is that the Cartesian grids intersect with the physical boundary with arbitrary fashion. This often leads to the so-called "cut-cell" problem [1], e.g. in 1D case the first grid is very closed to the physical boundary. If the boundary treatment is not well designed, it may require the time step to be extremely small for the sake of stability and result in the poor computation efficiency. For the moving boundary problems, no matter how the mesh is generated initially, the "cut-cell" problem is unavoidable.
One of the commonly used methods is to generate a body-fitted grid. For simple and fixed boundaries, one may use the smooth mapping between the Cartesian mesh and curvilinear mesh, while this kind of treatment often needs to change the original equations. The new equations can be more complicated than the original equations and also bring additional computation. For the time-varying geometric domain, complex grid management is required at each time step and it may be extremely expensive. Another approach is to use a Cartesian mesh, and take special treatment on the boundaries. There are several numerical boundary treatments listed as follows: embedded boundary method [2][3][4][5][6][7], immersed boundary method [8][9][10][11][12] and Inverse Lax-Wendroff (ILW) method [13][14][15][16][17][18][19][20], etc.
In this paper, we focus on the ILW boundary treatment for solving the moving boundary problems in 3D. The idea of ILW procedure comes from [21,22], in which the authors simulated the traffic and pedestrian flow by solving an Eikonal equation at each time step. They used the PDE repeatedly to convert the normal derivative into the tangent derivative at the boundary, and the ghost point values are defined by a Taylor expansion along the normal direction. Enlightened by this approach, Tan and Shu extended it to hyperbolic conservation law equations in [17]. The main idea of the ILW method is to convert the spatial derivative into a time derivative through the PDE, which is in contrast to the traditional Lax-Wendroff scheme, in which the time derivative is converted into spatial derivative (this is exactly the meaning of "inverse"). Numerical results demonstrate that the ILW method can be applied to construct ghost point values when solving equations with high order finite difference methods on a Cartesian mesh in a computational domain with complex geometries. However, the original ILW method proposed in [17] has the following drawbacks: 1. Heavy algebra in the ILW procedure for high dimensional nonlinear system. 2. For the junction of the inflow and outflow boundary, i.e, the sonic points, the ILW method should be coupled with the least square method to ensure the stability of the scheme. 3. The numerical solution obtained by the ILW method cannot satisfy the conservation property while the exact solution of original PDE has this property.
There are several papers to address the issues above. In 2012 Tan et al. proposed the simplified ILW (SILW) method in [19], which greatly reduced the computation cost of the ILW method for solving the system. Ding et al. [14] redefined the concept of conservation for the finite difference scheme, and gave a new ILW method that satisfies the new conservation concept. Lu et al. [16] proposed an ILW method that can deal with the sonic points, which can avoid the situation of the extremely small denominator. Except for hyperbolic conservation laws, the ILW procedure can also be applied to other types of equations. In [23] Filbet and Yang extended the ILW method to handle the Boltzmann equation, and Li et al. studied the ILW procedure for the diffusion problems, and the authors in [15,24,25] extended the ILW method to handle the convection-diffusion equation. Stability analyses for both the ILW and the simplified ILW procedures were given in [24][25][26][27], in which the authors discussed various situations and gave us the guiding ideology for designing a high order and stable numerical boundary treatment. In particular, theoretical analysis indicates that the carefully designed (S)ILW boundary treatment can maintain stability with the same CFL number (λ cfl ) max as in the periodic case, ignoring distances of the first grid point to the physical boundary. This means that the SILW method can overcome the "cut-cell" problem.
For moving boundary problems, in [18] Tan and Shu designed an ILW method for inviscid fluid with free-slip no-penetration moving boundary conditions, in which they only considered the translation of the moving boundary. Later, in [13] Cheng et al. reformulated the material derivative and considered both translation and rotation of the body. Recently, [28] extended the moving boundary ILW method to deal with the convection-diffusion equations and simulate the interaction between shock waves and rigid bodies in viscous fluids. In particular, a unified algorithm was given for five cases: pure convection, convection-dominated, convection-diffusion, diffusion-dominated and pure diffusion cases. This paper is an extension of our works [13,28], in which one and two dimensional moving boundary problems were considered. Now we extend these methods to handle three dimensional problems, with carefully designed numerical boundary conditions such that they can achieve third order of accuracy in smooth case, and no spurious oscillation when there is a shock near the surface. Particularly, in the ILW and SILW procedure, a local coordinate rotation on the boundary point is needed such that thex-axis of the new coordinate is in the same direction with the outward normal on the boundary. Different from the two dimensional problem, this rotation is not unique. We prove that the rotation will not affect the results via ILW procedure. Hence, the schemes can work with a simple form of rotation. Since the conditions are defined on the moving boundaries, material derivatives are used in the ILW procedure instead of the Eulerian time derivatives. However, this definition of material derivatives for functions limited on the boundaries is not clear. In this paper, we will follow the idea in [13] and reformulate the material derivatives on the moving boundaries with no-penetration condition in 3D.
The organization of the paper is as follows. In Section 2, we introduce the spatial discretization and time evolution method for the Eqs. (6) and (8). In Section 3, we focus on the boundary treatment and present the ILW method in detail. The method of tracking the moving rigid body is discussed in Section 4. Numerical results are given in Section 5 and concluding remarks are given in Section 6.

Scheme formulation inside (t)
We assume that the computational domain [ x l , x r ] ×[ y l , y r ] ×[ z l , z r ] is divided by a Cartesian mesh (x i , y j , z k ), and the mesh sizes in each direction are uniform, denoted by x, y and z, respectively, Denote W i,j,k (t) as the approximation to the point value W(x i , y j , z k , t). We will discretize the Eqs. (6) and (8) in the method of lines framework, meaning that the spatial variable is first discretized, then the numerical solution is updated in time by coupling a suitable time integrator.

Spatial discretization
The semi-discrete approximation of the governing Eq. (1) is given as whereF i+1/2,j,k ,Ĝ i,j+1/2,k andĤ i,j,k+1/2 are numerical fluxes such that the flux difference approximates the derivative to r-th order accuracy, Those numerical fluxes can be obtained by any reasonable finite difference scheme, such as the finite difference (FD) WENO scheme [29], which can achieve high order accuracy in smooth region and avoid numerical oscillatory near discontinuities. Here, we use the third order FD WENO scheme to obtain the numerical fluxes, i.e. r = 3. Details are given in Appendix A.
For the Euler Eq. (6), S i,j,k = 0 on the right hand side. And for the NS Eq. (8), S i,j,k is the approximation to the diffusion terms. Due to the effect of diffusion, we use central finite difference scheme without WENO to obtain S i,j,k . In particular, the five points central scheme will guarantee fourth order of accuracy, and the formulations for the (mixed) derivatives are given as follows:

Time discretization
The semi-discrete scheme (11) is equivalent to the first order ordinary differential equation (ODE) system with L({W i,j,k }) containing all spatial discritization terms. For time discretization, we use the third-order total variation diminishing (TVD) Runge-Kutta method [30]: In addition, a following constraint on time step t is imposed such that the boundary (t) can move at most one grid in each direction in a time step: where V b,x,max , V b,y,max and V b,z,max represent the maximum magnitude of the moving velocities in x-, y-and z-directions on the boundary (t).

Boundary treatments
In order to make the interior finite difference scheme work, we need to define enough ghost points near the boundaries. Near the outer boundaries x ∪ y ∪ z , which consist of straight lines, ghost points are assigned according to the principles of inflow, outflow and reflection boundary conditions. While near the inner boundary (t), we will assign carefully designed values on ghost points. In addition, there is a kind of points called "newly emerging" points that deserve special attention, which are outside the computational domain (t n−1 ) at the last time and entering (t n ) at the current time due to the movement of the boundary. In [18], the authors pointed out that the newly emerging points do not need any special treatment, but just set at most one more ghost point in each direction. Then, we can use the same scheme as the interior points to update the values of these points at the last stage of (14). And that is exactly the reason why we give a constraint on the time step in (15). Note that those ghost points may not lie on the grid. Hence, we follow the idea given in [17] that the method consists of the inverse Lax-Wendroff (ILW) type procedure for inflow boundary conditions and extrapolation for outflow boundary conditions. At the outflow boundary, WENO type extrapolation will be used in case of a shock close to the boundary. And the idea of the ILW procedure is utilizing the PDE repeatedly to write the normal derivatives of the variables W in terms of the time derivatives and tangential derivatives of W, both of which are obtained from boundary conditions. We should note that the boundary location varies with time, hence we will employ the material derivative rather than the Eulerian time derivative on the boundary. In the following subsections, we will concentrate on how to define the values of W i,j,k at the ghost points.

Boundary treatment for Euler equations
In the ILW method, we rewrite three-dimensional Euler equations in the form of primitive variables where

A(U), B(U) and C(U)
are Jacobian matrixes, given as For any ghost point P i,j,k = (x i , y j , z k ), we first find a point P a = (x a , y a , z a ) ∈ (t) so that its normal vector n outside the rigid body would pass through P i,j,k : y a = y c + r sin ϕ sin θ, The normal vector (from fluid to rigid body) of P a is n = (− sin ϕ cos θ, − sin ϕ sin θ, − cos ϕ).
Then, we rotate the coordinate axis locally at P a such that thex-axis of the new coordinate system is in the same direction with n. We want to remark the fact that the rotation is not unique due toŷ-andẑ-axis can be any orthogonal vectors in the tangent plane with respect to the rigid body at P a , but it will not effect the results via ILW procedure. Particularly, in Appendix B, we prove the following theorem: The ILW boundary procedure for the Euler equation is independent of the choice of the rotation, i.e, we can get the same boundary scheme even with different choices of rotation.
Therefore, we can get the same formulation even with different choices of tangent vectors. Here, we give one of these rotation methods and apply it to our later computation. The new local coordinate system (x,ŷ,ẑ) and the coordinate transformation matrix T are given as Consequently, we can prove that the Euler Eq. (16) has the same form in the new coordinate system as those in the original coordinate system: Next, we can approximate the valuesÛ i,j,k at point P i,j,k by a Taylor expansion at point P a . For instance, a third order approximation is given aŝ where r is a vector going from P a to P i,j,k ,Û a (orÛ a are the normal derivativesÛx| P a andÛxx| P a , respectively. It can be extended to arbitrary high order accuracy easily by keeping more terms in the Taylor expansion. For the current third-order scheme, we just need to get the valuesÛ (m) a , m = 0, 1, 2. Before that, a local characteristic decomposition should be done at P a , with the eigendecomposition where, is the diagonal matrix consisting of the eigenvalues of A(Û), and L forms from the corresponding left eigenvectors {l 1 , . . . , l 5 }: Notice that the point valueÛ a is unknown, we will use the WENO extrapolation to obtain the approximationÛ ext a (see Appendix C) and then substitute it into the formulations above, obtaining L a = L(Û ext a ). After that, convert the original variablesÛ to the characteristic fields V = L aÛ . Since the eigenvalues of A(Û) arê u − c,û,û,û,û + c, the moving speeds of five characteristic lines relative to the boundary normal direction are −c, 0, 0, 0, c, respectively.

Reformulation of the material derivative in 3D
In the ILW procedure, we convert normal spatial derivatives to tangential and time derivatives of the given boundary condition using the PDE. Specifically, for moving boundaries, we would employ the material derivatives of the boundary conditions. Details will be shown in Section 3.1.2. For any functiong ∈ C 1 (R 3 × R + ), the material derivative is given as where, (u, v, w) is the flow velocity. Unfortunately, the definition above can not be used for the functions defined on boundary only, such as the normal vector, since the spatial derivatives may not exist. Therefore, for the functions limited on the surface in threedimensional space, we will give a new definition of material derivative, which would play a key role in the ILW method in the following subsection. Suppose that the three-dimensional moving boundary (t) is given as where (τ 1 , τ 2 ) are parameters in [ τ 10 , τ 11 ] ×[ τ 20 , τ 21 ] . For example, the surface of a unit sphere in 3D can be described in spherical coordinates: Assume that, the parametric equation (t) is regular, i.e., ∂X b ∂τ 1 and ∂X b ∂τ 2 are linearly independent. Thus, ∂X b ∂τ 1 , ∂X b ∂τ 2 and n form a basis of the three-dimensional space. Then, the velocities of the fluid and the boundary can be written as follows: Moreover, it has the following property, which indicates the equivalence between our new definition (24) and the traditional definition of material derivative on no-penetration boundary. τ 2 , t)) is a smooth surface in three dimensions, and the no-penetration boundary condition holds on (t), i.e. u · n = V b · n. For anyg(x, y, z, t) ∈ C 1 (R 3 × R + ), we can define a function g on (t): Then the following equation is true: where D Dt is the material derivative given in (22).
Proof According to the formula (23) and u · n = V b · n, we know that u n = V b,n . Therefore, This completes the proof.
In the next subsection, we will use the definition L (24) on the given boundary conditions in the ILW process. For easy notation, we write L(g) as Dg Dt in the following subsections.

ILW method
With the boundary conditions and the material derivative given in Definition 1, we now proceed to constructÛ • At first, we want to find the value ofÛ (0) a . The no-penetration condition u =û · n = V b · n can be rewritten as Do the extrapolation of those outgoing characteristic variables, and we have the following equations a are obtained by the WENO extrapolation. Hence, those equations tell us a system as ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 0 1 0 0 0 l 2,1 l 2,2 l 2,3 l 2,4 l 2,5 l 3,1 l 3,2 l 3,3 l 3,4 l 3,5 l 4,1 l 4,2 l 4,3 l 4,4 l 4,5 l 5,1 l 5,2 l 5,3 l 5,4 l 5,5 Thus,Û (0) a can be obtained by solving the above system. • Then, we want to find the derivative ofÛ (1) a . From the second equation of Euler equations with primitive variables we can get So we getÛ (1) a by solving the following system: Again, (V j ) (1) a are obtained by the WENO extrapolation. • Finally, for the second derivativeÛ (2) a , we will employ the idea of simplified ILW (SILW) procedure proposed in [20], and the second derivative can be obtained by WENO extrapolation.

OnceÛ (m)
a , m = 0, 1, 2 are obtained, we can use (21) to getÛ i,j,k , and further U i,j,k and W i,j,k .

Coupling with RK time discretization
Notice that the numerical method we have described previously is for time level t n only. When coupling with the third order RK scheme (14), we need to match the time levels when constructing values of ghost points in the two intermediate stages W (1) and W (2) . In particular, the velocityû at the boundary at each RK stage is necessary as the boundary condition, and we need to updateû with the following formulations instead of the RK scheme (14) to maintain the third order in time: As in [18], we can use a standard Lax-Wendroff procedure to get the time derivatives as we want. More specifically, employing the second equation of (16) repeatedly, we can get the time derivatives ∂û ∂t and ∂ 2û ∂t 2 in terms of spatial derivatives. After finishing an inverse Lax Wendroff procedure at the time level t n , all the first and second order spatial derivatives values at x = X b (t n ) are given. Hence, we can get ∂û ∂t | x=X b (t n ),t=t n and ∂ 2û ∂t 2 | x=X b (t n ),t=t n with these spatial derivatives. More details can be found in [18].

Boundary treatments for NS equation
Similar to the case in Section 3.1, to define the value at a ghost point P i,j,k = (x i , y j , z k ), we set up a local coordinate system at P a by (18) at first. Then, consider the variables and system in the local coordinate system. Again, the rotation is not unique but will not affect the results of boundary scheme. The proof is very similar to the Euler equation in Appendix B. We will not repeat it here.
Due to the fact that the Dirichlet boundary conditions are given for velocityû and temperature T, we rewrite the Eq.
The matrixes and source term are given as and with Again, the valuesÛ i,j,k at point P i,j,k will be approximated by a Taylor expansion inxdirection at point P a . Once obtaining the point valueÛ a , and normal derivativesÛ (1) a and U (2) a , we can get the third order approximation ofÛ i,j,k via (21). Note that in the boundary treatment for the Euler equations in Section 3.1, the boundary valuesÛ (m) a are obtained via either the ILW method based on the boundary conditions or the extrapolation of the outgoing variables. However, for the compressible NS equation, totally different boundary treatments are needed for the diffusion-dominated and the convection-dominated regimes. Hence, [28] designed the numerical methods based on a careful combination of the boundary treatments for the two regimes and obtained a stable and accurate boundary condition for general convection-diffusion equations with either still or moving boundaries. Here, we extend the numerical method to our 3D model with the given isothermal no-slip wall boundary condition.

The material derivative in 3D
For the NS equation with moving boundary, the material derivative is also necessary in the ILW method.
For the no-slip boundary, the material derivatives for the boundary and fluid are the same, Moreover, since T b is constant, hence DT b /Dt = 0. And Dû/Dt = dV b /dt equal to the acceleration of the rigid body. As a consequence,

ILW method
In particular, the derivativesÛ (m) a , m = 1, 2, are a combination of those derived from the ILW procedure and the extrapolation procedure from interior points. For easy notation, we use the subscript "ilw" and "ext" to denote derivativesÛ (m) a obtained through the ILW procedureÛ with (Û) = diag{û − c,û,û,û,û + c}. Hence, considering the characteristic variables V = L aÛ , the components V 2 , · · · , V 5 are the outflow variables, and V 1 is the inflow variable. Again,Û ext a is used in the characteristic decomposition. • Construct the point valueÛ (0) a . Note that the boundary condition (9) gives the value ofû and T at the boundary (t) directly, so we can set Since V 5 is the outgoing variable, we have the relation that Therefore, we can obtain the value of (Û 1 ) (0) a as follows, • Construct the derivative valueÛ (1) a . Combining the definition p = ρ T and material derivative Dû/Dt in (34), we have that Given the outgoing performance of V 2 , . . . , V 5 , we have the relation that The second derivatives on the right hand side and (V j ) (1) a are obtained via extrapolation. Thus, we can haveÛ (1) a,ilw by solving the system. Then, with = max{ 1 Reρ , 4 3Reρ , γ PrReρ } and h = x 2 + y 2 + z 2 . • Construct the derivative valueÛ (2) a . Based on the boundary condition (34), the second derivatives (Û 2 ) (2) , . . . , (Û 5 ) (2) satisfy equations and we set them as "ilw": a,ilw = Reρ a,ilw = Reρ There is no diffusion term in the first equation of (34), so we just take we can solve the system to getÛ (2) a,ilw . Finally, we have that with { 1 , 2 , 3 , 4 , 5 } = 0, 4 3Reρ , 1 Reρ , 1 Reρ , γ PrReρ .

Coupling with RK time discretization
When coupling the third order SSP RK time discretization, we need the value of Û j t | P a ,t n and Û j tt | P a ,t n , j = 2, . . . , 5, to obtain high order accuracy in time via (29).
Fix at the boundary (x, y, z) = P a and t = t n , Look at the second material derivative Here, ( ∂û b ∂t , ∂v b ∂t , ∂ŵ b ∂t ) T are the acceleration of the boundary. In particular, the spatial derivatives Û j α and Û j αβ (α, β can bex,ŷ orẑ) at (P a , t n ) are already obtained via the ILW process t n . Hence, we only need to deal with Û j tx , Û j tŷ and Û j tẑ . Those values can be simulated through the standard Lax-Wendroff idea, i.e., turning the time derivative into spatial derivatives via PDE, and then approximating the spatial ones. However, the process would be very complex and the computational cost is large. Notice that the time derivative Û j t | i,j,k at time t n is already known. Hence, the mixed derivatives

Description of the moving rigid body
Under the impingement of the shock, the rigid body will move in the fluid. Here, both translation and rotation of the body should be taken into account. For each point , the velocity is given as where, V tr is the translational velocity, ω is the rotation velocity, and r is the vector from the boundary point to the centroid coordinates X c (t) = (x c (t), y c (t), z c (t)). Furthermore, V tr and ω satisfy the equation with the translational acceleration a tr and rotational acceleration a θ determined by Newton's second law and rotational law respectively. In particular, for the inviscid fluid, Here, M is the mass of the sphere, n is the external normal vector of the surface (from the rigid body to the fluid), and I is the moment of inertia of the sphere. While, for the viscous fluid, and we set the initial velocity of the moving cylinder as V b = (−0.5, 0, 0) T . As long as the solution stays smooth, it maintains the isentropic flow, i.e., s(x, y, z, t) = 1. In particular, solutions on any xy-section plane are the same as those in 2D. Hence, we use zero spatial derivative on the boundaries in z-direction to guarantee two dimensional property of the flow around the cylinder. We use the reflection boundary condition on the boundaries in x-and y-direction, and the moving boundary algorithm proposed in this paper is applied to deal with the cylinder boundary.
In Table 1, we give the entropy errors and the rate of convergence at t = 0.5 when the solution is smooth. It is observed that the scheme can achieve the designed third order accuracy.

Interaction between shock wave and a cylinder
This example is the three dimensional version of Example 4.3 in [13] and Example 5 in [28]. We use it to verify the correctness of our algorithm. As in the two papers, the computation area is taken as [ x l , x r ] ×[ y l , y r ] ×[ z l , z r ] =[ 0, 1] ×[ 0, 0.2] ×[ −0.5, 0.5]. At the initial time, a cylinder is parallel to y-direction with its center located at (0.15, 0.05, 0). The radius and the density are given as r = 0.05 and σ = 10.77, respectively. A shock wave is located at x = 0.08. We test the example with Mach number M s = 3. In z-direction, we use zero spatial derivative on the infinite boundaries to guarantee two dimensional property of the flow around the cylinder.
• Compressible inviscid fluid. We plot the pressure contour in Fig. 1 along the cross section of z = 0 at two typical time points.
• Compressible viscous fluid. We take Re = 1000, Pr = 0.7, and T b = 5 7 . We plot the pressure contour in Fig. 2 along the cross section of z = 0 at two typical time points.
Our numerical results consist with those of two dimensional ones in [13,28], indicating our method for three dimensional problems is correct.

Interaction between shock wave and a sphere
The computational domain is taken as [ The initial position of the spherical center is (0.15, 0.05, 0), with the radius r = 0.05 and the density σ = 10.77. And a shock wave is located at x = 0.08 initially. We test the example with M s = 3 and M s = 6.
In the early stage of computation, the sphere is impacted by the shock wave to produce a reflected shock wave, and the reflected shock wave hits the lower wall to generate a new reflected shock wave, causing uneven pressure distribution on the upper and lower sides of the sphere. The uneven pressure causes the sphere to move upward. Figures 3 and 4 show this procedure for inviscid fluid and viscous fluid respectively.
Later, the reflections of Mach shocks from the walls reach the sphere, generating more structures. Those structures will vary with the fluid. respectively. The pressure contours along the cross section of z = 0 are plotted. We can observe that there is a reflected shock wave when the incident shock wave interacts with the sphere. There is a second reflected shock wave when the reflected shock wave interacts with the upper wall. The second reflected shock wave hits the sphere again and forms a secondary interaction. The rigid body would be under greater stress in fluid with a larger Mach number. Consequently, the body would move faster. These numerical results show that our numerical method is robust in the computation for the supersonic flow with high Mach number.
• Compressible viscous fluid. We take Re = 500 and 10 7 , Pr = 0.7, and T b = 5 7 . The time evolution of the interaction of a plane shock wave and a sphere for Mach number M s = 3 and M s = 6 is plotted in Figs. 8, 9, 10, 11 and 12. Again, the system of shock waves is observed clearly, and results are significantly different with different Mach numbers. However, there is no discernible difference when changing Reynold number.
In addition, 3D pictures of the iso-surface of vorticity and the contour of pressure of the inviscid fluid and viscous fluid are given in Figs. 7 and 12, respectively. The numerical results demonstrate that our methods are stable and efficient for general cases.

Conclusion
In this paper, we develop the high order numerical boundary treatments for compressible inviscid fluid and viscous fluid, and further use them to simulate the interaction between a complex moving rigid body and a shock wave in 3D. Schemes are based on the finite difference methods on fixed Cartesian meshes in complex moving geometries. Values on the ghost points can be obtained through a Taylor expansion on the boundary in the normal direction. And we can construct point values and normal derivatives via ILW procedure or extrapolation, or their combination according to the equations and boundary conditions. This is an extension of our works [13,28], in which the boundary treatment was designed for 2D problems. Different from the 2D case, the local coordinate rotation required in the ILW procedure is not unique in 3D. We prove that the choice of the rotation will not affect results of the boundary treatment. Hence, we can choose a simple one in our simulation. On the other hand, the material derivative for inviscid fluid on the moving boundary with no-penetration condition is newly given. Consequently, we can simulate the body with both translation and rotation. In the numerical tests, we show the results of cylinder and sphere, indicating our algorithm is effective and robust. 3. Transform all the point values W j and fluxes F j = F(W j ) which are in the potential stencil of the WENO reconstruction to the local characteristic fields: 4. Do the flux splitting: where, α = max l |λ l (F (W))| takes over the relevant range of W, and λ l are the eigenvalues of the Jacobian F (W); 5. Perform the scalar WENO reconstruction procedure for each component of and the linear weights (b) Change the linear weights d r to nonlinear weights ω r : ω r = α r α 1 + α 2 , α r = d r (β r + ) 2 , , r = 1, 2, with = 10 −6 to avoid the denominator being zero, and the smoothness indicators β r on each small stencil are given as (c) Find the third order reconstruction (f l ) + i+1/2 = ω 1 (f l ) (1) i+1/2 + ω 2 (f l ) (2) i+1/2 ; 6. ConstructF − i+1/2 based onF − and the right-biased stencil. The process to obtaiñ F − i+1/2 is mirror-symmetric to that forF + i+1/2 , with respect to the target point x i+1/2 ; 7. Form the numerical flux aŝ Therefore, the ghost points can be defined by a Taylor expansion without tangent vectors.