# High-order compact gas-kinetic scheme in arbitrary Lagrangian-Eulerian
formulation

Yue Zhang<EMAIL_ADDRESS>Kun Xu<EMAIL_ADDRESS>Department of
Mathematics, Hong Kong University of Science and Technology, Clear Water Bay,
Kowloon, Hong Kong  Department of Mechanical and Aerospace Engineering, Hong
Kong University of Science and Technology, Clear Water Bay, Kowloon, Hong Kong
Shenzhen Research Institute, Hong Kong University of Science and Technology,
Shenzhen, China

###### Abstract

This study proposes an extension of the high-order compact gas-kinetic scheme
(CGKS) to compressible flow simulation in an arbitrary Lagrangian-Eulerian
(ALE) formulation in unstructured mesh. The ALE method is achieved by
subdividing arbitrary mesh into tetrahedrons and integrating flux function in
a local coordinate system at the cell interface to ensure geometric
conservation law. The scheme incorporates a compact reconstruction with third-
order accuracy for updating both cell-averaged conservative flow variables and
their gradients. HWENO-type nonlinear reconstruction and gradient compression
factors are adopted to improve the accuracy and robustness of the scheme. A
multi-stage multi-derivative (MSMD) time-stepping method is also implemented
to achieve high-order time accuracy with fewer middle stages. The scheme is
used to study problems involving moving boundaries. The numerical experiments
demonstrate the effectiveness of the scheme in capturing the accurate
solutions of both low-speed smooth flow and highly compressible ones with
strong shock waves.

###### keywords:

arbitrary Lagrangian-Eulerian (ALE) , geometric conservation law (GCL) ,
compact gas-kinetic scheme , multi-stage multi-derivative (MSMD) time
discretization

##  1 Introduction

There are two common frameworks for describing fluid motion in computational
fluid dynamics: the Eulerian and the Lagrangian systems. The Eulerian system
describes flow on a fixed grid, while in the Lagrangian system, the grid moves
with the fluid velocity. While the Eulerian system is more straightforward to
implement, it has difficulties in handling moving boundaries and introducing
dissipation when dealing with contact discontinuities. The Lagrangian system
can resolve contact discontinuities well, such as material interfaces or free
surfaces, but may lead to computational breakdown due to grid distortion.
Therefore, the development of arbitrary Lagrangian-Eulerian (ALE) methods [5]
becomes attractive. There are two types of ALE methods: direct and indirect.
Direct ALE methods involve moving both the object and the computational grid
simultaneously. Indirect ALE methods separate the motion of the object and the
computational grid into three stages: the Lagrangian stage, the rezoning stage
[2], and the remapping stage [11]. During the Lagrangian stage, the flow field
and grid are updated simultaneously. During the rezoning stage, a new
computational grid is constructed to match the object’s new position and
shape, often involving techniques such as grid generation and deformation to
ensure grid quality and accuracy. In the remapping stage, the solution is
transferred into the rezoned mesh. The current research focuses on the direct
ALE method with a fixed grid topology, which is designed to handle moving
boundaries and shock waves. Therefore, the remapping step can be omitted.

The gas-kinetic scheme (GKS) is a kinetic theory-based numerical method to
solve the Euler and Navier-Stokes equations [21]. Under the initial condition
of a generalized Riemann problem, a time-evolving gas distribution function is
constructed in GKS to calculate the numerical fluxes and evaluate the time-
dependent flow variables at a cell interface. As a result, both cell-averaged
flow variables and their gradients can be updated. Therefore, the HWENO-type
method and the two-step multi-resolution WENO reconstruction can be developed
in the scheme for the high-order spatial data reconstruction [27, 10]. At the
same time, due to the time-accurate flux function, the multi-stage multi-
derivative (MSMD) method can be used to update the solution with high-order
temporary accuracy. Specifically, the two-stage fourth-order (S2O4) time-
stepping method is used in the compact GKS (CGKS) [13, 12]. The CGKS has been
constructed on both structured and unstructured meshes in 2D and 3D cases [6,
9, 23, 24, 26]. In order to further improve the robustness of the scheme in
high-speed flow simulations, further improvements have been incorporated into
the scheme. First, the evolution of possible discontinuous flow variables at
different sides of a cell interface is constructed for updating reliable cell
averaged gradient of flow variables [25]. Second, the nonlinear limiting
process is implemented on the high-order time derivative of the flux function
under the MSMD framework. Equipped with the above remedies, the CGKS on a 3D
unstructured mesh is highly robust in hypersonic flow computation and a large
time step, such as CFL number $\approx 0.8$, can be used in the fourth-order
compact scheme. Alternatively, a gradient compression factor is designed to
improve the robustness and efficiency of CGKS in case of low-quality mesh [7],
and further improvements are studied in [22].

In recent years, high-order GKS has tended to incorporate ALE formulation. Ren
[18] developed a one-stage DG-ALE gas-kinetic method for oscillating airfoil
calculation. The variation of mesh velocity is considered in the time-
dependent flux calculation along a cell interface, and the geometric
conservation law is satisfied accurately. Pan developed a high-order ALE gas-
kinetic method on structure mesh for two and three dimensions[17, 16], where a
non-compact WENO scheme and the two-stage fourth-order discretization are
used. A bilinear interpolation is used to parameterize both grid coordinates
and grid moving velocity with the preservation of the geometric conservation
law. However, this interpolation isn’t a unified formulation for other types
of mesh, such as triangular prism and pyramid. In this paper, our target is to
develop a compact gas-kinetic scheme on unstructured mesh in ALE formulation.
Instead of bilinear interpolation, we use a subdivision of arbitrary meshes
into tetrahedrons and a linear interpolation in triangular surfaces to
preserve the geometry conservation law. A compact third-order WENO
reconstruction with the requirement of Neumann neighboring cells only is
adopted in spatial reconstruction. The radius basic function method [3] and
Lagrangian nodal solver [14] are used to determine grid velocity.

The structure of the paper is the following. Section 2 introduces the finite
volume method and GKS in ALE formulation. Section 3 discusses the grid
velocity, the two-stage fourth-order time integrating method, and the
geometric conservation law. Section 4 focuses on the compact reconstruction.
In Section 5, several test cases are used to validate the proposed method.
Finally, the conclusion is presented in the last section.

##  2 Gas-Kinetic Scheme in ALE Formulation

###  2.1 Governing equation on moving mesh

By using the finite volume method, the computational domain $\Omega$ is
discretized into a series of non-overlap elements $\Omega_{i}$,

$\Omega=\bigcup\Omega_{i},\ \Omega_{i}\bigcap\Omega_{j}=\phi(i\neq j).$

The boundary of element $\Omega_{i}$ can be expressed as

$\partial\Omega_{i}=\bigcup_{p=1}^{N_{f}}\Gamma_{ip},$

where $N_{f}$ is the number of element surface and $\Gamma_{ip}$ is surface of
the element. Reynolds transport theorem gives the conservation laws in an
arbitrary Lagrangian-Eulerian (ALE) formulation on a moving cell

$\frac{\mathrm{~{}d}}{\mathrm{~{}d}t}\int_{\Omega_{i}(t)}\bm{W}\mathrm{~{}d}\Omega+\int_{\partial\Omega_{i}(t)}(\bm{F}(\bm{W})-\bm{W}\bm{U})\cdot\bm{n}\mathrm{~{}d}S,$

where $\bm{U}$ is the mesh velocity, and $\bm{F}(\bm{W})$ is the flux pass
though the element surfaces. For Navier-Stokes equation, $\bm{W}=[\rho,\rho
V_{1},\rho V_{2},\rho V_{3},\rho E]^{T}$ is the conservative values, and
$\rho,\bm{V}=(V_{1},V_{2},V_{3}),E$ is the density, absolute velocity and
energy of fluid.

  * 1.

If $\bm{U}=\bm{V}$, the governing equation becomes Lagrangian form;

  * 2.

If $\bm{U}=\bm{0}$, the governing equation becomes Eulerian form.

A unified partitioning approach is adopted to handle different types of grids
and maintain consistency. Specifically, different types of grids (such as
hexahedrons, prisms, and pyramids) are all partitioned into several
tetrahedrons. The partitioning approach ensures computational consistency and
facilitates the conversion between different types of grids. An example of
partitioning a hexahedron into six tetrahedrons is illustrated in Figure 1.

Figure 1: An example of mesh partition

To maintain consistency with the partitioning of elements, we divide the
quadrilateral surface into two triangles. Under this condition, the
conservation law in ALE formulation can be written as

$\frac{\mathrm{~{}d}}{\mathrm{~{}d}t}\int_{\Omega(t)}\bm{W}\mathrm{~{}d}\Omega+\sum_{p_{t}=1}^{N_{t}}\int_{\Sigma_{p_{t}}}(\bm{F}(\bm{W})-\bm{W}\bm{U})\cdot\bm{n}\mathrm{~{}d}S=0,$

where $N_{t}$ is the number of triangles after partitioning and
$\Sigma_{p_{t}}$ is the individual triangle surfaces. The translation from
Eulerian coordinates to local two-dimensional coordinates $(\xi,\eta)$ as
shown in Figure 2 is used to evaluate the surface integration.

Figure 2: Translation a triangle surface to local coordinates

Thus, a linear interpolation is introduced for each triangle

$\bm{x}(\xi,\eta)=\bm{x}_{0}+\xi(\bm{x}_{1}-\bm{x}_{0})+\eta(\bm{x}_{2}-\bm{x}_{0}),$

where $\bm{x}_{0},\bm{x}_{1},\bm{x}_{2}$ is the coordinates of the triangle
vertices. With this transformation, the normal vector of triangle surfaces can
be written as

$\bm{n}=\bm{x}_{,\xi}\times\bm{x}_{,\eta}=\frac{\partial\bm{x}}{\partial\xi}\times\frac{\partial\bm{x}}{\partial\eta}.$

The surface integration can be written as

$\int_{\Sigma_{p_{t}}}(\bm{F}(\bm{W})-\bm{W}\bm{U})\cdot\bm{n}\mathrm{~{}d}S=\int_{0}^{1}\int_{0}^{1-\xi}(\bm{F}(\bm{W})-\bm{W}\bm{U})\cdot(\bm{x}_{,\xi}\times\bm{x}_{,\eta})\mathrm{~{}d}\eta\mathrm{~{}d}\xi.$

Assuming the grid velocity is time-independent, the coordinates of points on
the surface can be written as $\bm{x}=\bm{x}^{0}+\bm{U}t$, where $\bm{x}^{0}$
is the initial position of the points. Then, we have

$\begin{aligned}
\int_{0}^{1}\int_{0}^{1-\xi}\bm{W}\bm{U}\cdot(\bm{x}_{,\xi}\times\bm{x}_{,\eta})\mathrm{~{}d}\eta\mathrm{~{}d}\xi&=\int_{0}^{1}\int_{0}^{1-\xi}\bm{W}\bm{U}\cdot(\bm{x}^{0}_{,\xi}\times\bm{x}^{0}_{,\eta})\mathrm{~{}d}\eta\mathrm{~{}d}\xi\\\
&+t\int_{0}^{1}\int_{0}^{1-\xi}\bm{W}\bm{U}\cdot(\bm{U}^{0}_{,\xi}\times\bm{x}^{0}_{,\eta})\mathrm{~{}d}\eta\mathrm{~{}d}\xi\\\
&+t\int_{0}^{1}\int_{0}^{1-\xi}\bm{W}\bm{U}\cdot(\bm{x}^{0}_{,\xi}\times\bm{U}_{,\eta})\mathrm{~{}d}\eta\mathrm{~{}d}\xi\\\
&+t^{2}\int_{0}^{1}\int_{0}^{1-\xi}\bm{W}\bm{U}\cdot(\bm{U}_{,\xi}\times\bm{U}_{,\eta})\mathrm{~{}d}\eta\mathrm{~{}d}\xi\end{aligned}.$

The total flux can be approximated as

$\begin{aligned}
\int_{\Sigma_{p_{t}}}(\bm{F}(\bm{W})-\bm{W}\bm{U})\cdot\bm{n}\mathrm{~{}d}S&\approx\int_{0}^{1}\int_{0}^{1-\xi}(\bm{F}(\bm{W})-\bm{W}\bm{U})\cdot(\bm{x}^{0}_{,\xi}\times\bm{x}^{0}_{,\eta})\mathrm{~{}d}\eta\mathrm{~{}d}\xi\\\
&+t\int_{0}^{1}\int_{0}^{1-\xi}\bm{W}\bm{U}\cdot(\bm{U}^{0}_{,\xi}\times\bm{x}^{0}_{,\eta})\mathrm{~{}d}\eta\mathrm{~{}d}\xi\\\
&+t\int_{0}^{1}\int_{0}^{1-\xi}\bm{W}\bm{U}\cdot(\bm{x}^{0}_{,\xi}\times\bm{U}_{,\eta})\mathrm{~{}d}\eta\mathrm{~{}d}\xi\\\
&+t^{2}\int_{0}^{1}\int_{0}^{1-\xi}\bm{W}\bm{U}\cdot(\bm{U}_{,\xi}\times\bm{U}_{,\eta})\mathrm{~{}d}\eta\mathrm{~{}d}\xi\\\
&=\mathcal{L}_{p_{t}}(\bm{W},\Omega)+\mathbb{V}_{p_{t}}(t)\end{aligned},$

where $\mathcal{L}_{p_{t}}(\bm{W},\Omega)$ is the surface flux, and
$\mathbb{V}_{p_{t}}(t)$ represents the flux correction due to the motion of
grid, except the translation. Thus, the semi-discretized equation can be
written as

$\frac{\mathrm{~{}d}}{\mathrm{~{}d}t}\int_{\Omega_{i}(t)}\bm{W}\mathrm{~{}d}\Omega+\mathcal{L}(\bm{W},\Omega)+\mathbb{V}(t)=0.$
(1)

The surface integration can be calculated by Gaussian integration. Three
Gaussian points with local coordinates
$(\frac{1}{6},\frac{1}{6}),(\frac{2}{3},\frac{1}{6}),(\frac{1}{6},\frac{2}{3})$
and weightes $\omega_{1}=\omega_{2}=\omega_{3}=\frac{1}{3}$ are used in each
triangle surface. In this condition, the surface flux and correction can be
written as

$\displaystyle\mathcal{L}(\bm{W},\Omega)$
$\displaystyle=\sum_{p_{t}=1}^{Nt}\mathcal{L}_{p_{t}}(\bm{W},\Omega)=\sum_{p_{t}=1}^{Nt}\sum_{k=1}^{3}\omega_{k}\left[\bm{F}(\bm{W}_{p_{t},k})-\bm{W}_{p_{t},k}\bm{U}_{p_{t},k}\right]\cdot(\bm{x}^{0}_{,\xi}\times\bm{x}^{0}_{,\eta})_{p_{t}},$
$\displaystyle\mathbb{V}(t)$
$\displaystyle=\sum_{p_{t}=1}^{Nt}\mathbb{V}_{p_{t}}(t)=\sum_{p_{t}=1}^{Nt}\sum_{k=1}^{3}\omega_{k}\bm{W}_{p_{t},k}\bm{U}_{p_{t},k}\cdot\left[t(\bm{U}^{0}_{,\xi}\times\bm{x}^{0}_{,\eta})+t(\bm{x}^{0}_{,\xi}\times\bm{U}^{0}_{,\eta})+t^{2}(\bm{U}^{0}_{,\xi}\times\bm{U}^{0}_{,\eta})\right]_{p_{t}}.$

###  2.2 BGK equation in ALE Formulation

The gas-kinetic BGK equation can be written as

$\frac{\partial f}{\partial t}+\bm{v}\cdot\nabla_{x}f=\frac{g-f}{\tau},$

where $f=f(\bm{x},t,\bm{v},\xi)$ is the gas distribution function, $g$ is the
corresponding equilibrium state, and $\tau$ is the collision time.
$\bm{v}=(v_{1},v_{2},v_{3})$ is the particle velocity. The collision term in
the above equation describes the evolution process from a non-equilibrium
state to an equilibrium one with the satisfaction of compatibility condition

$\int\frac{g-f}{\tau}\Psi\mathrm{~{}d}\bm{w}\mathrm{~{}d}\Xi=0,$

where
$\Psi=(1,v_{1},v_{2},v_{3},\frac{1}{2}(v_{1}^{2}+v_{2}^{2}+v_{3}^{2}+\xi^{2}))^{T}$
and $\mathrm{~{}d}\Xi=\mathrm{~{}d}\xi_{1}\cdots\mathrm{~{}d}\xi_{K}$ ($K$ is
the number of internal degree of freedom, i.e. $K=(5-3\gamma)/(\gamma-1)$ for
three-dimensional flow and $\gamma$ is the specific heat ratio). According to
the Chapman-Enskog Expansion for BGK equation, the macroscopic governing
equations can be derived[20]. In the continuum region, the BGK equation can be
rearranged, and the gas distribution function can be expanded as

$f=g-\tau\mathcal{D}_{\bm{u}}g+\tau\mathcal{D}_{\bm{u}}(\tau\mathcal{D}_{\bm{u}})g-\tau\mathcal{D}_{\bm{u}}[\tau\mathcal{D}_{\bm{u}}(\tau\mathcal{D}_{\bm{u}})]g+\cdots,$

where $\mathcal{D}_{\bm{u}}=\frac{\partial}{\partial t}+\bm{u}\cdot\nabla$.
With the zeroth-order truncation $f=g$, the Euler equations can be obtained.
For the first-order truncation $f=g-\tau\mathcal{D}_{\bm{u}}g$, the Navier-
Stokes equations with the dynamic viscosity coefficient $\mu=\tau p$ ($p$ is
pressure) and Prandtl number $\text{Pr}=1$ can be obtained.

In this paper, a high-order gas-kinetic scheme will be constructed in the
arbitrary Lagrangian-Eulerian (ALE) framework for three-dimensional flow. The
particle speed related to grid $\bm{w}=(w_{1},w_{2},w_{3})$ is

$\bm{w}=\bm{v}-\bm{U},$

where $\bm{U}=(U_{1},U_{2},U_{3})$ is the grid velocity given aribitray. The
BGK equation on moving reference can be written as

$\frac{\partial f}{\partial t}+\bm{w}\cdot\nabla_{x}f=\frac{g-f}{\tau}.$ (2)

The flux in the term $\mathcal{L}$ of ALE-formulation finite volume method (1)
can be evaluated by

$\mathbb{F}(\bm{W})=(\bm{F}(\bm{W})-\bm{W}\bm{U})\cdot(\bm{x}^{0}_{,\xi}\times\bm{x}^{0}_{,\eta})=\int
f(\bm{x},t,\bm{w},\xi)\bm{w}\cdot\bm{n}\bm{\Psi}\mathrm{~{}d}\bm{v}\mathrm{~{}d}\Xi.$
(3)

In the actual computation, the numerical flux can be obtained by taking
moments of gas distribution function with related velocity

$\mathbb{F}(\bm{W})^{\prime}=\int
f(\bm{x},t,\bm{w},\xi)\bm{w}\cdot\bm{n}\bm{\Psi}^{\prime}\mathrm{~{}d}\bm{v}\mathrm{~{}d}\Xi,$
(4)

where
${\bm{\Psi}}=(1,{w}_{1},{w}_{2},{w}_{3},\frac{1}{2}({w}_{1}^{2}+{w}_{2}^{2}+{w}_{3}^{2}+\xi^{2}))^{T}$.
The component of $\mathbb{F}$ can be given by the combination of
$\mathbb{F}^{\prime}$ as follows

$\begin{cases}\mathbb{F}_{\rho}=\mathbb{F}_{\rho}^{\prime}\\\ \mathbb{F}_{\rho
V_{1}}=\mathbb{F}_{\rho W_{1}}^{\prime}+U_{1}\mathbb{F}_{\rho}^{\prime}\\\
\mathbb{F}_{\rho V_{2}}=\mathbb{F}_{\rho
W_{2}}^{\prime}+U_{2}\mathbb{F}_{\rho}^{\prime}\\\ \mathbb{F}_{\rho
V_{3}}=\mathbb{F}_{\rho W_{3}}^{\prime}+U_{3}\mathbb{F}_{\rho}^{\prime}\\\
\mathbb{F}_{\rho E}=\mathbb{F}_{\rho E}^{\prime}+U_{1}\mathbb{F}_{\rho
W_{1}}^{\prime}+U_{2}\mathbb{F}_{\rho W_{2}}^{\prime}+U_{3}\mathbb{F}_{\rho
W_{3}}^{\prime}+\frac{1}{2}(U_{1}^{2}+U_{2}^{2}+U_{3}^{2})\mathbb{F}_{\rho}^{\prime}\\\
\end{cases}.$

By defining the rotating matrix $\bm{T}=\text{diag}(1,\bm{T}^{\prime},1)$ and

$\bm{T}^{\prime}=\left(\begin{matrix}n_{1}&n_{2}&n_{3}\\\
-n_{2}&n_{1}+\frac{n_{3}^{2}}{1+n_{1}}&-\frac{n_{2}n_{3}}{1+n_{1}}\\\
-n3&-\frac{n_{2}n_{3}}{1+n_{1}}&1-\frac{n_{3}^{2}}{1+n_{1}}\end{matrix}\right),n_{1}\neq-1,$

and when $n_{1}=-1$, $\bm{T}^{\prime}$ becomes $\text{diag}(-1,-1,1)$, the
flux (Eq. 5) can be evaluated by

$\mathbb{F}(\bm{W})^{\prime}=\bm{T}^{-1}\int
f(\tilde{\bm{x}}_{p_{t},k},t,\tilde{\bm{w}},\xi)\tilde{w}_{1}\tilde{\bm{\Psi}}^{\prime}\mathrm{~{}d}\tilde{\bm{w}}\mathrm{~{}d}\Xi,$
(5)

where the origin point of local coordinate is $\tilde{\bm{x}}=(0,0,0)$ with
$x$-direction in $\bm{n}$ and
$\widetilde{\bm{\Psi}}^{\prime}=(1,\tilde{w}_{1},\tilde{w}_{2},\tilde{w}_{3},\frac{1}{2}(\tilde{w}_{1}^{2}+\tilde{w}_{2}^{2}+\tilde{w}_{3}^{2}+\xi^{2}))^{T}$.
The microscopic velocities in local coordinate are given by
$\tilde{\bm{w}}=\mathbf{T}^{\prime}\bm{w}$ and
$\tilde{w}_{1}=n_{1}w_{1}+n_{2}w_{2}+n_{3}w_{3}$.

###  2.3 Gas evolution model

In order to construct the numerical fluxes at $\bm{x}=(0,0,0)^{T}$, the
integral solution of the BGK equation Eq.(2) is

$f(\bm{x},t,\bm{w},\xi)=\frac{1}{\tau}\int_{0}^{t}g(\bm{x}^{\prime},t^{\prime}\bm{w},\xi)e^{-(t-t^{\prime})/\tau}\mathrm{~{}d}t^{\prime}+e^{t/\tau}f_{0}(\bm{x}_{0},\bm{w}),$
(6)

with the trajectory

$\bm{x}=\bm{x}^{\prime}+\bm{w}(t-t^{\prime}).$

In Eq.(6), $f_{0}$ is the initial gas distribution function, and $g$ is the
corresponding equilibrium state. $\bm{x}_{0}$ are the initial position by
tracing back particles $\bm{x}$ at time $t$ back to $t=0$.

Before constructing the initial distribution function $f_{0}$ and equilibrium
state $g$, we first denote

$\bm{a}=(a_{1},a_{2},a_{3})=\nabla_{x}g/g,A=g_{t}/g.$

In the following derivation, quadratic terms of time will be ignored directly.
With the consideration of possible discontinuity at an interface, the initial
distribution is constructed as

$f_{0}(\bm{x}_{0},\bm{w}_{0})=f_{0}^{l}(\bm{x}_{0},\bm{w})(1-\mathbb{H}(x_{1}))+f_{0}^{r}(\bm{x}_{0},\bm{w})\mathbb{H}(x_{1}),$
(7)

where $\mathbb{H}$ is the Heaviside function. $f_{0}^{l}$ and $f_{0}^{r}$ are
the initial gas distribution functions on the left and right sides of the
interface, which are determined by corresponding initial macroscopic variables
and their spatial derivatives. With the second-order accuracy,
$f_{0}^{k}(\bm{x},\bm{w})$ is constructed by Taylor expansion around
$(\bm{x},\bm{w})$

$f_{0}^{k}(\bm{x}_{0},\bm{w})=f_{G}^{k}(\bm{x},\bm{w})-\bm{w}t\cdot\nabla_{x}f_{G}^{k}$
(8)

for $k=l,r$. Based on the Chapman-Enskog expansion, $f_{G}^{k}$ is given by

$f_{G}^{k}=g^{k}[1-\tau(A^{k}+\bm{a}^{k}\cdot\bm{w})],$ (9)

where $g^{k}$ is the equilibrium distribution function defined by the
macroscopic variables $\bm{W}^{k}$ at both sides of a cell interface,
$\bm{a}^{k}$ are defined by the spatial derivatives of $g^{k}$

$\displaystyle a_{i}^{k}$ $\displaystyle=(\frac{\partial
g^{k}}{\partial\rho^{k}}\frac{\partial\rho^{k}}{\partial x_{i}}+\frac{\partial
g^{k}}{\partial W_{1}^{k}}\frac{\partial W_{1}^{k}}{\partial
x_{i}}+\frac{\partial g^{k}}{\partial W_{2}^{k}}\frac{\partial
W_{2}^{k}}{\partial x_{i}}+\frac{\partial g^{k}}{\partial
W_{3}^{k}}\frac{\partial W_{3}^{k}}{\partial x_{i}}+\frac{\partial
g^{k}}{\partial\lambda^{k}}\frac{\partial\lambda^{k}}{\partial x_{i}})/g^{k}$
$\displaystyle=a_{i1}^{k}+a_{i2}^{k}w_{1}+a_{i3}^{k}w_{2}+a_{i4}^{k}w_{3}+a_{i5}^{k}\frac{1}{2}(w_{1}^{2}+w_{2}^{2}+w_{3}^{2}\xi^{2}),$

and

$A^{k}=A_{1}^{k}+A_{2}^{k}w_{1}+A_{3}^{k}w_{2}+A_{4}^{k}w_{3}+A_{5}^{k}\frac{1}{2}(w_{1}^{2}+w_{2}^{2}+w_{3}^{2}\xi^{2})$

are determined by compatibility condition

$\int(f^{k}_{G}-g^{k})\Psi\mathrm{~{}d}\bm{w}\mathrm{~{}d}\Xi=0.$ (10)

Substituting Eq.(8) and (9) into (7), the initial gas distribution has
following form

$\displaystyle
f_{0}=\begin{cases}g^{l}\left[1-\left(\bm{a}^{l}\cdot\bm{w}\right)t-\tau\left(A^{l}+\bm{a}^{l}\cdot\bm{w}\right)\right],&x_{1}<0,\\\
g^{r}\left[1-\left(\bm{a}^{r}\cdot\bm{w}\right)t-\tau\left(A^{r}+\bm{a}^{r}\cdot\bm{w}\right)\right],&x_{1}\geq
0,\end{cases}$ (11)

Then, the equilibrium distribution is defined by the Taylor expansion

$\displaystyle g(\bm{x^{\prime}},t^{\prime},\bm{w})=$
$\displaystyle\overline{g}(\bm{x},0,\bm{w})+\nabla_{x}\overline{g}\cdot(\bm{x^{\prime}}-\bm{x})+\overline{g}_{t}t^{\prime}$
(12) $\displaystyle=$
$\displaystyle\overline{g}(\bm{x},0,\bm{w})-\nabla_{x}\overline{g}\cdot\bm{w}(t-t^{\prime})+\overline{g}_{t}t^{\prime}$
$\displaystyle=$
$\displaystyle\overline{g}(\bm{x},0,\bm{w})\left\\{1-\overline{\bm{a}}\cdot\bm{w}(t-t^{\prime})+At^{\prime}\right\\},$

where $\overline{g}$ and $\overline{\bm{a}}$ are determined from the
reconstruction of macroscopic flow variables presented in section 4.5, and $A$
is obtained by compatibility condition (10).

In smooth flow region, the collision time for viscous flow is determined by
$\tau=\mu/p$, where $\mu$ is the dynamic viscosity coefficient and $p$ is the
pressure at the cell interface, and for inviscid flow, the collision time
should be set as zero. In order to properly capture the un-resolved shock
structure, additional numerical dissipation is needed. The physical collision
time $\tau$ in the exponential function part can be replaced by a numerical
collision time $\tau_{n}$. For inviscid flow, it is set as

$\tau_{n}=C_{1}\Delta t+C_{2}\frac{|p_{l}-p_{r}|}{p_{l}+p_{r}}\Delta t,$

and for viscous flow, it is

$\tau_{n}=\tau+C_{2}\frac{|p_{l}-p_{r}|}{p_{l}+p_{r}}\Delta t,$

where $p_{l}$ and $p_{r}$ denote the pressure on the left and right sides of
the cell interface. In this paper, we have $C_{1}=0.01,C_{2}=5.0$.

By substituting Eq. (11) and Eq. (12) into Eq. (6) with $\tau$ and $\tau_{n}$,
the second-order time-dependent gas distribution function at a cell interface
becomes

$\displaystyle f\left(\bm{x},t,\bm{w},\xi\right)$
$\displaystyle=(1-e^{-t/\tau_{n}})\bar{g}+e^{-t/\tau_{n}}\left[\mathbb{H}\left(w_{1}\right)g^{l}+\left(1-\mathbb{H}\left(w_{1}\right)\right)g^{r}\right]+t\bar{A}\bar{g}$
(13)
$\displaystyle-\tau\left(1-e^{-t/\tau_{n}}\right)\bar{g}\left(\overline{\bm{a}}\cdot\bm{w}+\bar{A}\right)$
$\displaystyle-\tau
e^{-t/\tau_{n}}\mathbb{H}\left(w_{1}\right)g^{l}\left(\bm{a}^{l}\cdot\bm{w}+A^{l}\right)$
$\displaystyle-\tau
e^{-t/\tau_{n}}\left(1-\mathbb{H}\left(w_{1}\right)\right)g^{r}\left(\bm{a}^{r}\cdot\bm{w}+A^{r}\right)$
$\displaystyle+te^{-t/\tau_{n}}\left[\left(\overline{\bm{a}}\cdot\bm{w}\right)\bar{g}-\mathbb{H}\left(w_{1}\right)\left(\bm{a}^{l}\cdot\bm{w}\right)g^{l}-\left(1-\mathbb{H}\left(w_{1}\right)\right)\left(\bm{a}^{r}\cdot\bm{w}\right)g^{r}\right].$

The fluxes in Eq.(5) can be obtained by taking the moments of the above
distribution function. The calculation of moments can be found in [21].

In order to account for the Prandtl number not equal to 1, additional
correction heat flux should be added in energy flux:

$(\frac{1}{\text{Pr}}-1)q,$

where $q=F_{v4}-W_{1}F_{v1}-W_{2}F_{v2}-W_{3}F_{v3}$ can be obtained through
the viscous flux $\bm{F}_{vis}=[0,F_{v1},F_{v2},F_{v3},F_{v4}]$. The viscous
flux can be obtained by

$\bm{F}_{vis}=-\tau\int
w_{1}\bar{g}\left[\left(\overline{\bm{a}}\cdot\bm{w}\right)+\bar{A}\right]\mathrm{~{}d}{\bm{w}}\mathrm{~{}d}\Xi.$

###  2.4 Evolution of the cell-averaged spatial gradients

By taking moments of the above gas distribution function in Eq. (13), the
time-accurate conservative flow variables in the moving frame at the Gaussian
points can also be obtained

$\bm{W}_{p_{t},k}^{\prime}(t^{n+1})=\mathbf{T}^{\prime}\left(\int\widetilde{\bm{\Psi}}f(\tilde{\bm{x}}_{p_{t},k},t^{n+1},\tilde{\bm{w}},\xi)\mathrm{~{}d}\bm{w}\mathrm{~{}d}\Xi\right),$
(14)

by which the corresponding conservative flow variables in the absolute inertia
frame $\bm{W}_{p_{t},k}(t^{n+1})$ can be obtained. According to the Divergence
theorem, the cell averaged gradients over cell $\Omega_{i}$ at time $t^{n+1}$
are

$\overline{\nabla\bm{W}}^{n+1}_{i}=\frac{1}{|\Omega_{i}^{n+1}|}\sum_{p_{t}=1}^{N_{t}}\int_{\Sigma_{p_{t}}}\bm{W}^{n+1}\bm{n}_{p}dS,$
(15)

where the surface integration can be calculated by Gaussian quadrature

$\int_{\Sigma_{p_{t}}}\bm{W}^{n+1}\bm{n}_{p_{t}}dS\approx\sum_{k=1}^{3}|S_{p_{t}}|\omega_{k}\bm{W}_{p_{t},k}(t^{n+1})\bm{n}_{p_{t}}.$
(16)

Besides evaluating the cell averaged gradients, the solution updates of the
scheme are presented next.

##  3 Temporal discretization and Geometrical Conservation Law

###  3.1 Mesh Velocity

This article considers two types of mesh motion when the mesh velocity is not
given directly. The first method is the radial basic function (RBF) method,
which is used to address moving boundary problems. The second method is the
Lagrangian velocity method, which is used to handle moving shock problems.

####  3.1.1 Radial Basic Function (RBF) Method

For moving boundary problems, we can obtain the boundary moving vector in the
next step on boundary mesh, and the moving vectors of interior grids can be
obtained by radial basic function interpolation [3]. We chose the radial basic
function as

$\psi(\xi)=(1-\xi)^{4}(1+4\xi),$

where $\xi=|\bm{x}|/R$ is the related distance, and $R$ is the compactly
supported radius. $|\bm{x}|$ is calculated by
$\sqrt{x_{1}^{2}+x_{2}^{2}+x_{3}^{2}}$. Then, we can define an interpolation
function

$f(\bm{x})=\sum_{i}\lambda_{i}\psi(|\bm{x}-\bm{x_{i}}|/R),$

where the $\bm{x}_{i}$ is a series of given points, and

$f(\bm{x_{j}})=\sum_{i}\lambda_{i}\psi(|\bm{x_{j}}-\bm{x_{i}}|/R)=g_{j},$

so we can solve $\lambda_{i}$ in above equation. The whole process of RBF
interpolation is shown as Algorithm 1.

Algorithm 1 RBF interpolation

0: This algorithm requires a series of points $\bm{x}_{i}$ and their
corresponding data $g_{i}$, as well as the number of points $N$ to be selected
for interpolation, a compactly supported radius $R$, and a target
interpolation point $\bm{x}_{0}$.

0: This algorithm ensures that an interpolated result $g_{0}$ is returned.

1: Select the closest $N$ points to point $\bm{x}_{0}$.

2: Calculate an $N\times N$ matrix $M_{ij}$ using the formula
$M_{ij}=\psi(|\bm{x}_{i}-\bm{x}_{j}|/R)$.

3: Solve the equation $M_{ij}\lambda_{j}=g_{i}$ by obtaining the values of
$\lambda_{j}$.

4: Initialize $g_{0}$ to 0.

5: for $i$ in range($N$) do

6: Compute $\xi_{i}=|\bm{x}_{i}-\bm{x}_{0}|/R$.

7: if $\xi_{i}<1.0$ then

8: Add $\lambda_{i}\psi(\xi_{i})$ to $g_{0}$.

9: end if

10: end for

11: return $g_{0}$

####  3.1.2 Lagrangian Velocity Method

The mesh velocity can be determined by the cell-centered Lagrangian nodal
solver [14]. The basic idea is to solve two half one-dimensional Riemann
problems at each cell interface by assuming that its velocity equals the point
velocity. We denote $C(p)$ as the set of cells $c$ that share the common
vertex $p$ and $F_{p}(c)$ as the set of faces of cell $c$ that share the
common vertex $p$. The grid speed $\bm{U}_{p}$ can be given by

$\bm{U}_{p}=\mathbb{M}_{p}^{-1}\sum_{c\in C(p)}\sum_{f\in
F_{p}(c)}[S_{f}p_{c}\bm{N}_{f}^{c}+\mathbb{M}_{pcf}\bm{V}_{c}],$

where $p_{c}$ and $\bm{V}_{c}$ are the pressure and velocity of cell $c$ and
$S_{f},\bm{N_{f}}$ is the surface area and normal direction. The matrix
$\mathbb{M}_{pcf}$ is calculated by

$\mathbb{M}_{pcf}=S_{f}\rho_{c}a_{c}(\bm{N}_{f}^{c}\otimes\bm{N}_{f}^{c}),$

where $a_{c}$ and $\rho_{c}$ are the sound speed and density of cell $c$ and
$\otimes$ is the operator of tensor production, and matrix $\mathbb{M}_{p}$ is
the sum of $\mathbb{M}_{pcf}$ for faces which is a symmetric positive definite
matrix. In this case, the coordinate of points in the next time step can be
calculated by

$\bm{x}_{p}^{n+1}=\bm{x}_{p}^{n}+\bm{U}_{p}\Delta t.$

With the Lagrangian velocity, the meshes become distorted, and computation may
break down. Therefore, the meshes need to be smoothed by

$\widetilde{\bm{x}}_{p}=(1-\omega)\bm{x}_{p}+\omega\frac{1}{N_{p}}\sum_{ip\in
N(p)}\bm{x}_{ip},$

where $N(p)$ is the set of vertexes jointing with vertex $p$, $N_{p}$ is the
size of $N(p)$, and $\omega$ is relaxation coefficient.

###  3.2 Solution Update and Geometrical Conservation Law

According to the semi-discretization equation Eq. (1), the two-stage fourth-
order (S2O4) temporal discretization is adopted here for the solution updates
[13],

$\displaystyle(\Omega_{i}\bm{W}_{i})^{*}=$
$\displaystyle(\Omega_{i}\bm{W}_{i})^{n}+\frac{1}{2}\Delta
t\mathcal{L}\left(\bm{W}_{i}^{n},\Omega^{n}\right)+\frac{1}{8}\Delta
t^{2}\frac{\partial}{\partial
t}\mathcal{L}\left(\bm{W}_{i}^{n},\Omega^{n}\right)+\int_{0}^{1/2\Delta
t}\mathbb{V}(t)\mathrm{~{}d}t,$ $\displaystyle(\Omega_{i}\bm{W}_{i})^{(n+1)}=$
$\displaystyle(\Omega_{i}\bm{W}_{i})^{n}+\Delta
t\mathcal{L}\left(\bm{W}_{i}^{n},\Omega^{n}\right)+\frac{1}{6}\Delta
t^{2}\left(\frac{\partial}{\partial
t}\mathcal{L}\left(\bm{W}_{i}^{n},\Omega^{n}\right)+2\frac{\partial}{\partial
t}\mathcal{L}\left(\bm{W}_{i}^{*},\Omega^{*}\right)\right)+\int_{0}^{\Delta
t}\mathbb{V}(t)\mathrm{~{}d}t,$

which is a two-step time integration method. And constant grid velocities are
considered during the time-integrating process. For simplicity, the
conservative value at Gaussian points appears in correction term
$\mathbb{V}(t)$ set as the moments of equilibrium state $\bar{g}$.

With a uniform flow field, the time evolution of the cell moving volume can be
calculated as

$\displaystyle(\Omega_{i})^{*}$
$\displaystyle=(\Omega_{i})^{n}+\frac{1}{2}\Delta
t\sum_{p_{t}=1}^{Nt}\sum_{k=1}^{3}\omega_{k}\bm{U}_{p_{t},k}\cdot(\bm{x}^{0}_{,\xi}\times\bm{x}^{0}_{,\eta})_{p_{t}}$
$\displaystyle+\frac{1}{4}\Delta
t^{2}\sum_{p_{t}=1}^{Nt}\sum_{k=1}^{3}\omega_{k}\bm{U}_{p_{t},k}\cdot\left[\frac{1}{2}(\bm{U}^{0}_{,\xi}\times\bm{x}^{0}_{,\eta})+\frac{1}{2}(\bm{x}^{0}_{,\xi}\times\bm{U}^{0}_{,\eta})+\frac{1}{6}\Delta
t(\bm{U}^{0}_{,\xi}\times\bm{U}^{0}_{,\eta})\right]_{p_{t}},$
$\displaystyle(\Omega_{i})^{n+1}$ $\displaystyle=(\Omega_{i})^{n}+\Delta
t\sum_{p_{t}=1}^{Nt}\sum_{k=1}^{3}\omega_{k}\bm{U}_{p_{t},k}\cdot(\bm{x}^{0}_{,\xi}\times\bm{x}^{0}_{,\eta})_{p_{t}}$
$\displaystyle+\Delta
t^{2}\sum_{p_{t}=1}^{Nt}\sum_{k=1}^{3}\omega_{k}\bm{U}_{p_{t},k}\cdot\left[\frac{1}{2}(\bm{U}^{0}_{,\xi}\times\bm{x}^{0}_{,\eta})+\frac{1}{2}(\bm{x}^{0}_{,\xi}\times\bm{U}^{0}_{,\eta})+\frac{1}{3}\Delta
t(\bm{U}^{0}_{,\xi}\times\bm{U}^{0}_{,\eta})\right]_{p_{t}},$

which automatically satisfies the Geometrical Conservation Law (GCL) when the
volume of the element is taken as the sum of the tetrahedral volumes.

The gas distribution function at Gaussian points Eq. (13) can be written as a
form

$f=f+tf_{t},$ (17)

which denotes the second-order time accuracy of the gas distribution function.
A two-step evolution of the gas distribution function can be constructed,
which is updated by

$\displaystyle f^{*}$ $\displaystyle=f^{n}+\frac{1}{2}\Delta tf_{t}^{n},$
$\displaystyle f^{n+1}$ $\displaystyle=f^{n}+\Delta tf_{t}^{*}.$

Then, the time-dependent conservative values at each Gauss point can be
obtained by Eq. (14). Then, by Eq. (15) and Eq. (16), the cell-averaged slopes
can be updated.

###  3.3 The Algorithm of ALE Method

Algorithm 2 presents the complete computational algorithm, where the bold text
indicates the special treatments compared to the Eulerian formulation. The
mesh velocity is determined once per time step, but the geometry requires the
update in every substep of the two-stage fourth-order temporal discretization.
This necessitates the updating of the volume and reconstruction matrix of
cells, as well as the normal vector, area, and Gaussian points coordinates of
interfaces in every substep.

Algorithm 2 CGKS in ALE formulation

1: while the computation uncompleted do

2: calculate the time step $\Delta t$ according to CFL number

3: determine mesh velocity

4: for i=1,2(for S2O4) do

5: define boundary condition for ghost cell

6: reconstruct the cell distribution

7: evolution for interfaces

8: update the geometry

9: update cell average conservative values and the first-order spatial
derivatives of conservative values

10: end for

11: end while

##  4 HWENO Reconstruction

The 3rd-order compact reconstruction [8] is adopted here with cell-averaged
values and cell-averaged first-order spatial derivative. In order to capture
shock, WENO weights [28] and gradient compression factor (CF) [7] are used.
The reconstruction was further improved [22], where only one large stencil and
one sub-stencil are involved in the new WENO procedure.

###  4.1 Third-order compact reconstruction for large stencil

Firstly, a linear reconstruction is presented. A quadratic polynomial $p^{2}$
with a third-order accuracy in space is constructed as follows

$\displaystyle p^{2}$
$\displaystyle=a_{0}+\frac{1}{h}[a_{1}(x-x_{0})+a_{2}(y-y_{0})+a_{3}(z-z_{0})]$
$\displaystyle+\frac{1}{2h^{2}}[a_{4}(x-x_{0})^{2}+a_{5}(y-y_{0})^{2}+a_{6}(z-z_{0})^{2}]$
$\displaystyle+\frac{1}{h^{2}}[a_{7}(x-x_{0})(y-y_{0})+a_{8}(y-y_{0})(z-z_{0})+a_{9}(x-x_{0})(z-z_{0})],$

where $h=\frac{|\Omega_{0}|}{\max_{j}S_{j}}$($|\Omega_{0}|$ is the volume of
cell and $S_{j}$ is the area of cell’s surface ) is the cell size, and
$(x_{0},y_{0},z_{0})$ is the coordinate of cell center.

The $p^{2}$ on $\Omega_{0}$ is constructed on the compact stencil $S$
including $\Omega_{0}$ and its all von Neumann neighbors $\Omega_{m}$
($m=1,\cdots,N_{f}$, where $N_{f}=6$ for hexahedron cell or $N_{f}=5$ for
triangular prism). The cell averages $\overline{Q}$ over $\Omega_{0}$ and
$\Omega_{m}$ and cell averages of space partial derivatives
$\overline{Q}_{x},\overline{Q}_{y}$ and $\overline{Q}_{z}$ over $\Omega_{m}$
are used to obtain $p^{2}$.

The polynomial $p^{2}$ is required to exactly satisfy cell averages over both
$\Omega_{0}$ and $\Omega_{m}$ ($m=1,\cdots,N_{f}$)

$\iiint_{\Omega_{0}}p^{2}\text{d}V=\overline{Q}_{0}|\Omega_{0}|,\iiint_{\Omega_{m}}p^{2}\text{d}V=\overline{Q}_{m}|\Omega_{m}|,$

with the following condition satisfied in a least-square sense

$\displaystyle\iiint_{\Omega_{m}}\frac{\partial}{\partial
x}p^{2}\text{d}V=\left(\overline{Q}_{x}\right)_{m}|\Omega_{m}|$
$\displaystyle\iiint_{\Omega_{m}}\frac{\partial}{\partial
y}p^{2}\text{d}V=\left(\overline{Q}_{y}\right)_{m}|\Omega_{m}|$
$\displaystyle\iiint_{\Omega_{m}}\frac{\partial}{\partial
z}p^{2}\text{d}V=\left(\overline{Q}_{z}\right)_{m}|\Omega_{m}|.$

The constrained least-square method is used to solve the above system.

###  4.2 Green-Gauss reconstruction for the sub stencil

The classical Green-Gauss reconstruction with only cell-averaged values is
adopted to provide the linear polynomial $p^{1}$ for the sub-stencil.

$p^{1}=\overline{Q}+\frac{1}{|\Omega_{0}|}\bm{x}\cdot\sum_{m=1}^{N_{f}}\frac{\overline{Q}_{m}+\overline{Q}_{0}}{2}S_{m}\bm{n}_{m},$

where $S_{m}$ is the area of the cell’s surface and $\bm{n}_{m}$ is the
surface’s normal vector.

###  4.3 Gradient compression Factor

The CF was first proposed in [7], and several improvements have been made in
[22]. Denote $\alpha_{i}\in[0,1]$ as gradient compression factor at targeted
cell $\Omega_{i}$

$\alpha_{i}=\prod_{p_{t}=1}^{n_{t}}\prod_{k=0}^{M_{p}}\alpha_{p_{t},k},$

where $\alpha_{p_{t},k}$ is the CF obtained at the $k$th Gaussian point of the
interface $p_{t}$ around cell $\Omega_{i}$, which can be calculated by

$\displaystyle\alpha_{p_{t},k}=\frac{1}{1+A^{2}},$ $\displaystyle
A=\frac{|p^{l}-p^{r}|}{p^{l}}+\frac{|p^{l}-p^{r}|}{p^{r}}+(\text{Ma}^{l}_{n}-\text{Ma}^{r}_{n})^{2}+(\text{Ma}^{l}_{t}-\text{Ma}^{r}_{t})^{2},$

where $p$ is pressure, $\text{Ma}_{n}$ and $\text{Ma}_{t}$ are the Mach
numbers defined by normal and tangential velocity, and superscript $l,r$
denote the left and right values of the Gaussian points.

Then, the updated slope is modified by

$\widetilde{\overline{\nabla\bm{W}}}_{i}^{n+1}=\alpha_{i}\overline{\nabla\bm{W}}_{i}^{n+1},$

and the Green-Gauss reconstruction is modified as

$p^{1}=\overline{Q}+\alpha\frac{1}{|\Omega_{0}|}\bm{x}\cdot\sum_{m=1}^{N_{f}}\frac{\overline{Q}_{m}+\overline{Q}_{0}}{2}S_{m}\bm{n}_{m}.$

###  4.4 Non-linear WENO weights

In order to deal with discontinuity, the idea of multi-resolution WENO
reconstruction is adopted [7, 28]. Here only two polynomials are chosen

$P_{2}=\frac{1}{\gamma_{2}}p^{2}-\frac{\gamma_{1}}{\gamma_{2}}p^{1},P_{1}=p^{1}.$

Here, we choose $\gamma_{1}=\gamma_{2}=0.5$. So, the quadratic polynomial
$p^{2}$ can be written as

$p^{2}=\gamma_{1}P_{1}+\gamma_{2}P_{2}.$ (18)

Then, we can define the smoothness indicators

$\beta_{j}=\sum_{|\alpha|=1}^{r_{j}}\Omega^{\frac{2}{3}|\alpha|-1}\iiint_{\Omega}\left(D^{\alpha}p^{j}(\mathbf{x})\right)^{2}\mathrm{~{}d}V,$

where $\alpha$ is a multi-index and $D$ is the derivative operator, and
$r_{1}=1$ and $r_{2}=2$. Special care is given for $\beta_{1}$ for better
robustness

$\beta_{1}=\min(\beta_{1,\text{Green-Gauss}},\beta_{1,\text{least-square}}),$

where $\beta_{1,\text{Green-Gauss}}$ is the smoothness indicator defined by
Green-Gauss reconstruction, and $\beta_{1,\text{least-square}}$ is the
smoothness indicator defined by second-order least-square reconstruction.
Then, the smoothness indicators $\beta_{i}$ are non-dimensionalized by

$\tilde{\beta}_{i}=\frac{\beta_{i}}{Q_{0}^{2}+\beta_{1}+10^{-40}}.$

The nondimensionalized global smoothness indicator $\tilde{\sigma}$ can be
defined as

$\tilde{\sigma}=\left|\tilde{\beta}_{1}-\tilde{\beta}_{0}\right|.$

Therefore, the corresponding nonlinear weights are given by

$\tilde{\omega}_{m}=\gamma_{m}\left(1+\left(\frac{\tilde{\sigma}}{\epsilon+\tilde{\beta}_{m}}\right)^{2}\right),\epsilon=10^{-5},$
$\bar{\omega}_{m}=\frac{\tilde{\omega}_{m}}{\sum\tilde{\omega}_{m}},m=1,2.$

Replacing $\gamma_{m}$ in equation (18) by $\bar{\omega}_{m}$ , the final non-
linear reconstruction can be obtained

$R(\bm{x})=\bar{\omega}_{2}P_{2}+\bar{\omega}_{1}P_{1}.$

The desired non-equilibrium states at Gaussian points become

$Q_{p,k}^{l,r}=R^{l,r}\left(\bm{x}_{p,k}\right),\left(Q_{x_{i}}^{l,r}\right)_{p,k}=\frac{\partial
R^{l,r}}{\partial x_{i}}\left(\bm{x}_{p,k}\right).$

###  4.5 Reconstruction of equilibrium states

After reconstructing the non-equilibrium state, a kinetic weighted average
method can be used to get equilibrium states and tangential derivatives [8],

$\int\bar{g}\Psi\mathrm{~{}d}\bm{v}\mathrm{~{}d}\Xi=\bm{W}_{0}=\int_{\tilde{w}_{1}>0}g^{l}\Psi\mathrm{~{}d}\bm{v}\mathrm{~{}d}\Xi+\int_{\tilde{w}_{1}<0}g^{r}\Psi\mathrm{~{}d}\bm{v}\mathrm{~{}d}\Xi,$
$\int\bar{a}_{i}\bar{g}\Psi\mathrm{~{}d}\bm{v}\mathrm{~{}d}\Xi=\frac{\partial\bm{W}_{0}}{\partial\tilde{x}_{i}}=\int_{\tilde{w}_{1}>0}a_{i}^{l}g^{l}\Psi\mathrm{~{}d}\bm{v}\mathrm{~{}d}\Xi+\int_{\tilde{w}_{1}<0}a_{i}^{r}g^{r}\Psi\mathrm{~{}d}\bm{v}\mathrm{~{}d}\Xi,i=2,3.$

For the normal derivatives, the above solution is further modified according
to the idea in linear diffusive generalized Riemann problem (dGRP) [4]

$\int\bar{a}_{1}\bar{g}\Psi\mathrm{~{}d}\bm{v}\mathrm{~{}d}\Xi=\frac{\partial\bm{W}_{0}}{\partial\tilde{x}_{1}}=\int_{\tilde{w}_{1}>0}a_{1}^{l}g^{l}\Psi\mathrm{~{}d}\bm{v}\mathrm{~{}d}\Xi+\int_{\tilde{w}_{1}<0}a_{1}^{r}g^{r}\Psi\mathrm{~{}d}\bm{v}\mathrm{~{}d}\Xi+\frac{\bm{W}^{r}-\bm{W}^{l}}{(\bm{x}_{rc}-\bm{x}_{lc})\cdot\bm{n}},$
(19)

where $\bm{x}_{rc}$ and $\bm{x}_{lc}$ are the coordinates of left and right
cell centroid,$\bm{W}^{l}$ and $\bm{W}^{r}$ are the reconstructed left and
right conservative value at Gaussian points and $\bm{n}$ is the normal vector
of interface. By adding a penalty term in Eq. (19), the whole scheme is
essentially free from the odd-even decoupling phenomenon [1].

##  5 Numerical Ecperiments

Both two-dimensional and three-dimensional cases are tested in this section.
For two-dimensional problems, the three-dimensional solver is used with one
layer and mirror boundary condition in the $z$ direction. The time step is
given by $\Delta t=\min\Delta t_{i}$, where $\Delta t_{i}$ is defined in each
cell

$\Delta t_{i}=C_{CFL}\min(\frac{h}{|V|_{i}+c_{i}},\frac{h^{2}}{3\nu_{i}}),$

where $C_{\text{CFL}}$ is the CFL number, $|V|_{i},c_{i}$ and
$\nu_{i}=(\mu/\rho)_{i}$ is the magnitude of related velocities, sound speed
and kinematic viscosity coefficient of cell $i$. Here, we set the CFL number
as 0.5 if not specified.

###  5.1 Accuracy Test

The advection of density perturbation for three-dimensional flow is presented
to test the order of accuracy. The computation domain is
$[0,1]\times[0,1]\times[0,1]$. Three uniform meshes with
$10^{3},20^{3},40^{3}$ tetrahedral cells are used in this case. The flow at
time $t$ is

$\displaystyle\rho(x,y,z)=1+0.2\sin(2\pi x_{r})\sin(2\pi y_{r})\sin(2\pi
z_{r}),$ $\displaystyle p(x,y,z)=1,V_{1}=V_{2}=V_{3}=1,$

where the $(x_{r},y_{r},z_{r})$ depends on time

$(x_{r},y_{r},z_{r})=(x,y,z)-(1,1,1)t.$

In order to validate the order of accuracy with moving meshes, the time-
dependent moving meshes are considered in this form

$\begin{cases}x&=x_{0}+0.05\sin(\pi t)\sin(2\pi x_{0})\sin(2\pi
y_{0})\sin(2\pi z_{0})\\\ y&=y_{0}+0.05\sin(\pi t)\sin(2\pi x_{0})\sin(2\pi
y_{0})\sin(2\pi z_{0})\\\ z&=z_{0}+0.05\sin(\pi t)\sin(2\pi x_{0})\sin(2\pi
y_{0})\sin(2\pi z_{0}).\\\ \end{cases}$

The errors and numerical orders are shown in Table 1, which shows that third-
order accuracy is achieved under moving mesh.

Table 1: The errors and accuracy in moving mesh mesh | $Error_{L^{1}}$ | $O_{L^{1}}$ | $Error_{L^{2}}$ | $O_{L^{2}}$ | $Error_{L^{\infty}}$ | $O_{L^{\infty}}$
---|---|---|---|---|---|---
$10^{3}$ | 1.40E-2 |  | 1.83E-02 |  | 5.17E-02 |
$20^{3}$ | 2.25E-3 | 2.64 | 2.96E-03 | 2.63 | 8.02E-03 | 2.69
$40^{3}$ | 3.04E-4 | 2.88 | 3.94E-04 | 2.91 | 1.08E-03 | 2.89

As a reference, the accuracy test result on stationary mesh is shown in Table
2, which indicates that the errors of moving and fixed mesh are compatible
with each other.

Table 2: The errors and accuracy in stationary mesh mesh | $Error_{L^{1}}$ | $O_{L^{1}}$ | $Error_{L^{2}}$ | $O_{L^{2}}$ | $Error_{L^{\infty}}$ | $O_{L^{\infty}}$
---|---|---|---|---|---|---
$10^{3}$ | 1.38E-02 |  | 1.81E-02 |  | 5.16E-02 |
$20^{3}$ | 2.07E-03 | 2.74 | 2.79E-03 | 2.70 | 7.92E-03 | 2.70
$40^{3}$ | 2.68E-04 | 2.95 | 3.61E-04 | 2.95 | 1.05E-03 | 2.92

Next, geometric conservation law, which is about the maintenance of uniform
flow passing through the moving mesh, is tested. The same mesh and same mesh
moving method are set in the test. The initial condition is set as

$\rho(x,y,z)=1,p(x,y,z)=1,V_{1}=V_{2}=V_{3}=1.$

The density errors at $t=1$ are shown in Table 3. The results indicate that
the errors are reduced to machine zero, implying the satisfaction of the
geometric conservation law.

Table 3: The errors of geometric conservation law in moving mesh mesh | $Error_{L^{1}}$ | $Error_{L^{2}}$ | $Error_{L^{\infty}}$
---|---|---|---
$10^{3}$ | 1.20148e-15 | 1.49810e-15 | 6.21725e-15
$20^{3}$ | 1.92325e-15 | 2.49943e-15 | 1.43219e-14
$40^{3}$ | 4.19724e-15 | 6.13169e-15 | 9.50351e-14

###  5.2 Flow inside a cylinder

This case is provided by High-Fidelity CFD Verification Workshop [19]. It is
about the flow inside a cylinder’s interior volume. Five meshes (with mesh
number $4\times 4\times 5,8\times 8\times 5,16\times 16\times 5,32\times
32\times 5,64\times 64\times 5$) are tested in this case. The center of the
cylinder is initially located on the original point, and the initial radius of
the cylinder is $r_{0}=0.5$. The motion of meshes is defined by

$\left[\begin{array}[]{l}x\\\ y\\\
1\end{array}\right]=\left[\begin{array}[]{ccc}1&0&0\\\ 0&1&\alpha(t)\\\
0&0&1\end{array}\right]\left[\begin{array}[]{ccc}\cos\left(A_{\theta}\alpha(t)\right)&-\sin\left(A_{\theta}\alpha(t)\right)&0\\\
\sin\left(A_{\theta}\alpha(t)\right)&\cos\left(A_{\theta}\alpha(t)\right)&0\\\
0&0&1\end{array}\right]\left[\begin{array}[]{ccc}\psi(t)&0&0\\\
0&\frac{1}{\psi(t)}&0\\\
0&0&1\end{array}\right]\left[\begin{array}[]{c}r_{0}\cos\left(\theta_{g}(t)\right)\\\
r_{0}\sin\left(\theta_{g}(t)\right)\\\ 1\end{array}\right],$

as the composite of four motions (translation, rotation, const-volume
deformation, and non-unit geometry mapping Jacobian). The translation motion
is a motion in $y$ direction as $\alpha(t)$. And the rotation motion is about
the center of the cylinder with a rotation amplitude $A_{\theta}=\pi$. The
transformation of the cylinder deforming into an ellipse such that the
interior area remains constant during deformation is facilitated by the
function

$\psi(t)=1+\left(A_{a}-1\right)\alpha(t),$

which varies from 1 to $A_{a}=1.5$ over $t=[0,2]$. Next, a parameterized
function $\eta$ is defined as

$\eta(\lambda,\omega,\epsilon)=\sin(\omega\lambda+\epsilon(1-\cos(\omega\lambda))),$

which is designed to break spatial and temporal symmetries due to the fact
that the integral of $\eta$ over a period $\omega$ does not equal 0 for
appropriate values of $\epsilon$. The parameter $\lambda$ represents the
independent variable (e.g. $t$ or $\theta$), whereas $\omega$ and $\epsilon$
represent the function frequency and shape characteristics. A subsequent
function $f_{g}$ is defined here to prescribe a deformation that ensures non-
unit geometry mapping Jacobians ($g\neq 1$) while utilizing the function
$\eta$ to break spatial and temporal symmetries. This is given by

$f_{g}\left(t,r_{0},\theta_{0}\right)=\left(16r_{0}^{4}+\frac{t^{6}}{t^{6}+0.01}\eta(t,10,0.7)\left(\cos\left(32\pi
r_{0}^{4}\right)-1\right)\right)\eta\left(\theta_{0},1,0.7\right)$

The prescribed deformation takes the form of a perturbation in $\theta$ as

$\theta_{g}\left(t,r_{0},\theta_{0}\right)=\theta_{0}+A_{g}f_{g}\left(t,r_{0},\theta_{0}\right),$

which the volume deformation amplitude $A_{g}$ is set as 0.15. The origin
mesh, deforming mesh when $\alpha=0$ and $\alpha=1$ are shown in Figure 3,
which indicates the large deforming and non-symmetry during the simulation.

(a) The origin mesh (b) The deformed mesh: $\alpha=0$ (c) The deformed mesh:
$\alpha=1$

Figure 3: An example mesh with cell number $8\times 8\times 5$

The cylinder interior is prescribed with a no-slip, adiabatic wall boundary
condition. The Reynolds number defined by cylinder diameter $d=1$ and
reference velocity $V=1.0$ is $\text{Re}=1000$, a Prandtl number of 0.72, and
constant viscosity is used. The initial condition is given by a constant zero
velocity field, which $\rho=1$ and $\rho E=50$.

####  5.2.1 Short Time Study

The time-activation function $\alpha(t)$ for this case is defined by

$\alpha(t)=t^{3}(8-3t)/16,$

which varies from 0 to 1 on the interval $t=[0,2]$. This case runs from $t=0$
until $t=1$. The time history of force in $x$ and $y$ direction on the surface
of the cylinder and power defined by

$P=\int_{surface}\bm{U}\cdot\bm{f}_{surf}\mathrm{~{}d}s,$

where $\bm{f}_{surf}$ is the surface stress vector, are shown in Figure 4.

Figure 4: The time history force and power of short-time motion

The results indicate that the solutions obtained from Mesh 3, Mesh 4, and Mesh
5 are in good agreement with each other. This suggests that we have achieved
mesh-independent solutions.

####  5.2.2 Long-time Study of Conservation

The motion of the long-time study has the same form as that of the composite
motion for the short-time study but with a different $\alpha(t)$ function. The
$\alpha(t)$ is selected as

$\alpha(t)=0.3\sin(6t)\frac{t^{6}}{t^{6}+1}.$

The purpose of the second term in the equation is to start the sinusoidal
motion slowly at $t=0$. We simulate the process from $t=0$ to $t=40$. The mesh
with $32\times 32\times 5$ cells is used in this study. The time history of
mass, force in $y$ direction on the surface of the cylinder, and power are
shown in Figure 5.

Figure 5: The time history force and power of long-time motion

The total mass of the entire domain remains constant at 0.785065 throughout
the process, which indicates the mass conservation properties of our method.
Additionally, the power and Y-Force exhibit periodic behavior after $t=30$.

###  5.3 Heaving-pitching airfoil

This case considers a NACA0012 airfoil undergoing a smooth flapping-type
motion. The airfoil moves smoothly upward for a distance of one chord length
over a duration of $T=2$ time units by heaving and pitching about a point
located at the airfoil’s 1/3 chord location. The functions governing the
motion are provided below:

$\displaystyle\Delta h(t)$ $\displaystyle=t^{3}(8-3t)/16,$
$\displaystyle\Delta\theta(t)$
$\displaystyle=\frac{80\pi}{180}(-t^{6}+6t^{5}-12t^{4}+8t^{3}).$

Two types of motion are considered in this case. The first type of motion is
defined solely by the heaving function $\Delta h(t)$, while the second type of
motion is defined by both the heaving function $\Delta h(t)$ and the pitching
function $\Delta\theta(t)$.

The income flow is set as $\rho_{\infty}=1.0,T_{\infty}=25,U_{\infty}=1$,
which has a Mach number of 0.2. The Reynolds number based on the chord of the
airfoil and incoming velocity is $\text{Re}=1000$, a Prandtl number of 0.72,
and a constant viscosity is used. As shown in figure 6, three meshes are used
in this case.

(a) Mesh 1

(b) Mesh 2

(c) Mesh 3

Figure 6: Mesh of NACA 0012

The adiabatic non-slip wall is set on the airfoil surface, and the far-field
boundary condition is set as the outer boundary. The initial condition is the
steady-state solution for the initial position $\Delta h=\Delta\theta=0$. The
radius basic function method is used for mesh motion.

Figure 7 displays the time history of the force in the $x$ and $y$ directions,
as well as the power. The results indicate that the three meshes produce
nearly identical solutions for the force in the $y$ direction and power.
Although the solutions for the force in the $x$ direction appear to differ
slightly among the meshes, the difference is negligible since the force in
this direction is very small.

Figure 7: NACA0012: force and power history of Motion 1 (translation)

Figure 8 illustrates the vorticity contour and mesh distribution at $t=1.5$.
At this time, the airfoil has moved approximately 0.738 times its chord length
in the $y$ direction, while the mesh surrounding the airfoil remains of good
quality. The vorticity contour indicates the generation of a vortex at the
front edge of the airfoil.

(a) Vorticity

(b) Mesh

Figure 8: NACA0012: the solution of Motion 1 (translation) at time $t=1.5$

Figure 9 displays the force and power history of Motion 2. The results
indicate that the solutions obtained from Mesh 2 and Mesh 3 are in good
agreement, suggesting that the solutions are mesh-independent.

Figure 9: NACA0012: force and power history of Motion 2 (translation +
rotation)

Figure 10 depicts the solution at time $t=0.5$ when the airfoil has rotated
33.75 degrees. At this time, the mesh remains of high quality with only slight
distortion at the trailing edge. The vorticity contour shows the shedding of a
vortex at the trailing edge.

(a) Vorticity

(b) Mesh

Figure 10: NACA0012: the solution of Motion 2 (translation + rotation) at time
$t=0.5$

###  5.4 2D Noh Problem

This case is a typical test case for verifying ALE code. A two-dimensional
domain with dimensions $[0,1]\times[0,1]$ is considered. An unstructured mesh
with mesh size $h=1/40$ is used for the computation, as shown in Figure 11(a).
The initial conditions are set to uniform density $\rho=1$ and internal energy
$e=\frac{p}{\gamma-1}=1\times 10^{-4}$, with $\gamma=5/3$. The initial
velocity is given by $\bm{V}=(-x/r,-y/r)$, where $r$ is the distance from the
origin. An axially symmetric shock is generated at the origin and propagates
at a constant speed. At time $t=0.6$, the shock is located at $r=0.2$. An
exact solution of density [15] is

$\rho=\begin{cases}16,&r<0.2\\\ 1+t/r,&r>0.2\end{cases}.$

An asymmetric boundary condition is enforced at the lines $x=0$ and $y=0$,
while a non-reflective boundary condition is used at all other boundaries.
Mesh movement is achieved using Lagrangian velocity, and a smoothing process
is applied every 20 time-steps with a relaxation coefficient $\omega=0.4$. The
CFL number is set to 0.3 for this case. The final mesh is shown in Figure
11(b), and it is observed that the mesh is concentrated in the post-shock
region.

(a) Initial mesh

(b) Final mesh

Figure 11: The mesh used in Noh Problem

The density contour at time $t=0.6$ is shown in Figure 12(a), and the
distribution with respect to $r$ of the entire domain is presented in Figure
12(b), which agrees well with the exact solution.

(a) The density contour

(b) The density distribution with respect to $r$

Figure 12: The result of Noh problem at $t=0.6$

###  5.5 Sedov Problem

The Sedov blast wave is a standard hydrodynamics test problem that models a
blast wave emanating from a singular energy point. It is a benchmark problem
for both the Lagrangian velocity method and the ALE method. The computational
domain is $[0,1.2]^{3}$, and two uniform grids with $20^{3}$ and $40^{3}$
cells are used in this test case. The fluid is modeled by the ideal gas
equation of state with $\gamma=1.4$. The initial condition is a uniform static
field with density $\rho=1$ and pressure $p=1e-6$, except for the cell
containing the origin. In this cell, the initial pressure is defined by
$p=(\gamma-1)\epsilon_{0}/V$, where $\epsilon_{0}=0.106384$ is the total
amount of released energy and $V$ is the cell volume. Asymmetric boundary
condition is enforced at the plane $x=0$, $y=0$and $z=0$, while non-reflecting
boundary conditions are used at the plane $x=1.2$, $y=1.2$and $z=1.2$. Mesh
movement is accomplished using Lagrangian velocity, and a smoothing process is
applied every ten time-steps with a relaxation coefficient of $\omega=0.8$. A
small CFL number of 0.01 is used to avoid the instability caused by the
initial singularity at the origin. After ten steps, the simulation transitions
to a normal CFL number of 0.3. Figure 13 displays the mesh distribution and
density contour, which reveal that the mesh is refined adaptively around the
shock.

(a) Mesh $20^{3}$

(b) Mesh $40^{3}$

Figure 13: Sedov problem: the mesh distribution and density contour at $t=1$

Figure 14 displays the density and pressure distribution along the diagonal
from $(0,0,0)$ to $(1.2,1.2,1.2)$, as well as the solutions of a uniform
stationary mesh.

(a) Density (b) Pressure

Figure 14: Sedov problem: the density and pressure distribution along diagonal
line

The results show that the solution of the ALE method agrees much more closely
with the exact solution than the static method. The solutions at the extreme
points of the ALE method are closer to the exact solution than those of the
uniform mesh. Additionally, the fine mesh ALE solution provides a more
accurate solution for the shock position than the coarse mesh one.

###  5.6 Saltzmann Problem

The Saltzmann problem is a test case that simulates fluid in a cylinder with a
piston. It examines the motion of a planar shock on a skewed Cartesian grid.
The computational domain is $[0,1]\times[0,0.1]\times[0,0.1]$, with uniform
$h=0.01$ in three directions. A mesh distortion of $x$-coordinate is defined
by

$\Delta x=(0.5y+z-15yz)\sin(\pi x).$

The initial mesh is shown in Figure 15.

Figure 15: The initial mesh of Saltzmann problem

Initially, an ideal monatomic gas with $\rho=1,e=10^{-4},\gamma=5/3$ is filled
in the box. The wall on the left-hand side moves at a constant speed of 1 from
left to right, acting as a piston. The other boundaries act as inviscid
reflecting walls. As the gas is compressed, a shock is generated and
propagates to the right faster than the piston. Mesh movement is accomplished
using Lagrangian velocity, and a smoothing process is applied every 20 time
steps with a relaxation coefficient of $\omega=0.8$. Figure 16 displays the
entire solution of the domain at different times, indicating that the solution
tends towards a one-dimensional solution.

(a) Density (b) Pressure

Figure 16: Saltzmann problem: the density and pressure distribution at
different time

The drop in density may be due to the entropy increase at the wall, and
further high-order boundary conditions will be considered. At $t=0.6$, the
shock is located at $x=0.8$, and the post-shock solutions are $\rho=4$ and
$p=1.333$. At time $t=0.75$, the shock will reach the right boundary and
bounce back due to the wall. Figure 17 demonstrates that the Lagrangian
velocity method is capable of handling shock motion issues, as it shows the
mesh distribution at different time intervals.

(a) $t=0.6$

(b) $t=0.85$

Figure 17: Saltzmann problem: mesh and density distribution at different time

##  6 Conclusion

In this paper, a third-order compact gas-kinetic scheme is developed in an
arbitrary Lagrangian-Eulerian (ALE) formulation for three-dimensional
unstructured mesh. In order to achieve geometric consistency, the arbitrary
mesh is subdivided into tetrahedrons, and the quadrilateral surface is split
into two triangles. By accounting for the variation of grid velocity within
the interface surface and dealing with it in a local framework, the geometric
conservation law is preserved. The kinetic equation at the moving frame with
grid velocity is used to construct a time-accurate gas distribution function
at the Gaussian point of the cell interface, allowing for updating both cell-
averaged flow variables and their gradients for the use in the high-order
compact initial reconstruction. The scheme is validated using several test
cases. For close boundary problems, mass conservation is well-preserved. With
the help of radial basis function interpolation, the method can handle moving
boundary problems and obtain high-fidelity solutions. For blast problems with
strong discontinuities, the scheme retains robustness, and the Lagrangian
nodal solver is adapted to track material interfaces.

## Acknowledgements

The current research is supported by the National Key R&D Program of China
2022YFA1004500, National Natural Science Foundation of China (Grant No.
12172316), and Hong Kong research grant council (16208021, 16301222).

## References

  * Blazek [2015] Blazek, J., 2015.  Chapter 4 - structured finite-volume schemes, in: Blazek, J. (Ed.), Computational Fluid Dynamics: Principles and Applications (Third Edition). third edition ed.. Butterworth-Heinemann, Oxford, pp. 73–120.
  * Brackbill and Saltzman [1982] Brackbill, J., Saltzman, J., 1982\.  Adaptive zoning for singular problems in two dimensions.  Journal of Computational Physics 46, 342–368.  doi:10.1016/0021-9991(82)90020-1.
  * De Boer et al. [2007] De Boer, A., Van der Schoot, M.S., Bijl, H., 2007.  Mesh deformation based on radial basis function interpolation.  Computers & structures 85, 784–795.
  * Gassner et al. [2007] Gassner, G., Lörcher, F., Munz, C.D., 2007.  A contribution to the construction of diffusion fluxes for finite volume and discontinuous Galerkin schemes.  Journal of Computational Physics 224, 1049–1063.
  * Hirt et al. [1974] Hirt, C., Amsden, A., Cook, J., 1974.  An arbitrary Lagrangian-Eulerian computing method for all flow speeds.  Journal of Computational Physics 14, 227–253.  doi:10.1016/0021-9991(74)90051-5.
  * Ji et al. [2018] Ji, X., Pan, L., Shyy, W., Xu, K., 2018.  A compact fourth-order gas-kinetic scheme for the Euler and Navier–Stokes equations.  Journal of Computational Physics 372, 446–472.
  * Ji et al. [2021a] Ji, X., Shyy, W., Xu, K., 2021a.  A gradient compression-based compact high-order gas-kinetic scheme on 3D hybrid unstructured meshes.  International Journal of Computational Fluid Dynamics 35, 485–509.
  * Ji et al. [2020] Ji, X., Zhao, F., Shyy, W., Xu, K., 2020.  A HWENO reconstruction based high-order compact gas-kinetic scheme on unstructured mesh.  Journal of Computational Physics 410, 109367.
  * Ji et al. [2021b] Ji, X., Zhao, F., Shyy, W., Xu, K., 2021b.  Compact high-order gas-kinetic scheme for three-dimensional flow simulations.  AIAA Journal 59, 2979–2996.
  * Ji et al. [2021c] Ji, X., Zhao, F., Shyy, W., Xu, K., 2021c.  Two-step multi-resolution reconstruction-based compact gas-kinetic scheme on tetrahedral mesh.  arXiv:2102.01366 [physics] ArXiv: 2102.01366.
  * Kucharik and Shashkov [2014] Kucharik, M., Shashkov, M., 2014\.  Conservative multi-material remap for staggered multi-material arbitrary Lagrangian–Eulerian methods.  Journal of Computational Physics 258, 268–304.  doi:10.1016/j.jcp.2013.10.050.
  * Li [2019] Li, J., 2019.  Two-stage fourth order: temporal-spatial coupling in computational fluid dynamics (cfd).  Advances in Aerodynamics 1, 1–36.  doi:10.1186/s42774-019-0004-9.
  * Li and Du [2016] Li, J., Du, Z., 2016.  A two-stage fourth order time-accurate discretization for Lax–Wendroff type flow solvers I. hyperbolic conservation laws.  SIAM Journal on Scientific Computing 38, A3046–A3069.
  * Maire and Nkonga [2013] Maire, P.H., Nkonga, B., 2013\.  Multi-scale Godunov-type method for cell-centered discrete Lagrangian hydrodynamics.  Journal of Computational Physics , 799–821doi:10.1016/j.jcp.2008.10.012.
  * Noh [1987] Noh, W., 1987.  Errors for calculations of strong shocks using an artificial viscosity and an artificial heat flux.  Journal of Computational Physics 72, 78–120.  doi:10.1016/0021-9991(87)90074-X.
  * Pan and Xu [2021] Pan, L., Xu, K., 2021.  An Arbitrary-Lagrangian-Eulerian high-order gas-kinetic scheme for three-dimensional computations.  Journal of Scientific Computing 88, 1–26.  doi:10.1007/s10915-021-01525-9.
  * Pan et al. [2020] Pan, L., Zhao, F., Xu, K., 2020.  High-order ALE gas-kinetic scheme with WENO reconstruction.  Journal of Computational Physics 417, 109558.  doi:10.1016/j.jcp.2020.109558.
  * Ren et al. [2016] Ren, X., Xu, K., Shyy, W., 2016.  A multi-dimensional high-order DG-ALE method based on gas-kinetic theory with application to oscillating bodies.  Journal of Computational Physics 316, 700–720.  doi:10.1016/j.jcp.2016.04.028.
  * Wukie et al. [2023] Wukie, N.A., Fidkowski, K., Persson, P.O., Wang, Z.J., 2023\.  High-fidelity CFD verification workshop 2024: Mesh motion, in: AIAA SCITECH 2023 Forum, p. 1243\.
  * Xu [1998] Xu, K., 1998.  Gas-kinetic schemes for unsteady compressible flow simulations.  Computational Fluid Dynamics, Annual Lecture Series, 29 th, Rhode-Saint-Genese, Belgium .
  * Xu [2001] Xu, K., 2001.  A gas-kinetic BGK scheme for the Navier–Stokes equations and its connection with artificial dissipation and Godunov method.  Journal of Computational Physics 171, 289–335.
  * Zhang et al. [2023] Zhang, Y., Ji, X., Xu, K., 2023.  A high-order compact gas-kinetic scheme in a rotating coordinate frame and on sliding mesh.  International Journal of Computational Fluid Dynamics , 1–20.
  * Zhao et al. [2019] Zhao, F., Ji, X., Shyy, W., Xu, K., 2019.  Compact higher-order gas-kinetic schemes with spectral-like resolution for compressible flow simulations.  Advances in Aerodynamics 1, 13\.
  * Zhao et al. [2020] Zhao, F., Ji, X., Shyy, W., Xu, K., 2020.  An acoustic and shock wave capturing compact high-order gas-kinetic scheme with spectral-like resolution.  International Journal of Computational Fluid Dynamics 34, 731–756.
  * Zhao et al. [2021] Zhao, F., Ji, X., Shyy, W., Xu, K., 2021.  Direct modeling for computational fluid dynamics and the construction of high-order compact scheme for compressible flow simulations.  arXiv:2107.06555 [physics] ArXiv: 2107.06555.
  * Zhao et al. [2022] Zhao, F., Ji, X., Shyy, W., Xu, K., 2022.  A compact high-order gas-kinetic scheme on unstructured mesh for acoustic and shock wave computations.  Journal of Computational Physics 449, 110812.
  * Zhu and Qiu [2018] Zhu, J., Qiu, J., 2018.  New finite volume weighted essentially nonoscillatory schemes on triangular meshes.  SIAM Journal on Scientific Computing 40, A903–A928.
  * Zhu and Shu [2020] Zhu, J., Shu, C.W., 2020.  A new type of third-order finite volume multi-resolution WENO schemes on tetrahedral meshes.  Journal of Computational Physics 406, 109212.