An ALE-type discrete unified gas kinetic scheme for low-speed continuum and rarefied flow simulations with moving boundaries

TL;DR

This paper extends the discrete unified gas kinetic scheme (DUGKS) to an ALE framework for simulating low-speed continuum and rarefied flows with moving boundaries, ensuring accuracy and efficiency across flow regimes.

Contribution

The paper introduces an ALE-type DUGKS with a remapping-free scheme and discusses geometric conservation approaches, enabling accurate simulations of moving boundary flows.

Findings

Accurate simulation of flows around oscillating objects.

Good agreement with experimental and numerical results.

Effective handling of both continuum and rarefied flow regimes.

Abstract

In this paper, the original discrete unified gas kinetic scheme (DUGKS) is extended to arbitrary Lagrangian-Eulerian (ALE) framework for simulating the low-speed continuum and rarefied flows with moving boundaries. For ALE method, the mesh moving velocity is introduced into the Boltzmann-BGK equation. The remapping-free scheme is adopted to develop the present ALE-type DUGKS, which avoids the complex rezoning and remapping process in traditional ALE method. As in some application areas, the large discretization errors will be introduced into the simulation if the geometric conservation is not guaranteed. Three compliant approaches of the geometric conservation law (GCL) are discussed and a uniform flow test case is conducted to validate these schemes. To illustrate the performance of present ALE-type DUGKS, four test cases are carried out. Two of them are the continuum flow cases, which are the flows around the oscillating circular cylinder and the pitching NACA0012 airfoil, respectively. Others are the rarefied flow cases, one is the moving piston driven by the rarefied gas, another is the flow caused by the plate oscillating in its normal direction. The results of all test cases are in good agreement with the other numerical and/or experimental results, demonstrating the capability of present ALE-type DUGKS to cope with the moving boundary problems at different flow regimes.

Tables

Table 1. Table I: Spatial errors of numerical solutions for uniform flow

$△$$x$(cell size)	DGCL scheme	Pressure		U		V
		$L_2$	$L_{\infty }$	$L_2$	$L_{\infty }$	$L_2$	$L_{\infty }$
1.000	Scheme1	7.44$e$-14	1.17$e$-14	2.29$e$-13	4.07$e$-14	3.38$e$-13	6.02$e$-14
	Scheme2	1.15$e$-13	2.14$e$-14	5.98$e$-13	1.19$e$-13	7.47$e$-13	1.21$e$-13
	Scheme3	8.07$e$-14	5.88$e$-15	1.21$e$-13	2.04$e$-14	9.86$e$-14	2.99$e$-14
0.500	Scheme1	3.00$e$-13	1.52$e$-14	2.95$e$-13	3.48$e$-14	2.56$e$-13	2.09$e$-14
	Scheme2	2.63$e$-13	9.55$e$-15	2.81$e$-13	3.29$e$-14	2.25$e$-13	1.85$e$-14
	Scheme3	2.95$e$-13	1.02$e$-14	2.27$e$-13	2.15$e$-14	1.90$e$-13	1.31$e$-14
0.250	Scheme1	1.10$e$-12	4.10$e$-14	7.55$e$-13	3.84$e$-14	6.79$e$-13	2.67$e$-14
	Scheme2	1.04$e$-12	1.64$e$-14	6.48$e$-13	2.96$e$-14	5.73$e$-13	1.80$e$-14
	Scheme3	9.63$e$-13	1.52$e$-14	6.23$e$-13	2.56$e$-14	5.48$e$-13	1.51$e$-14
0.125	Scheme1	3.33$e$-12	5.71$e$-14	1.96$e$-12	4.38$e$-14	1.91$e$-12	4.45$e$-14
	Scheme2	2.64$e$-12	2.10$e$-14	1.45$e$-12	2.02$e$-14	1.47$e$-12	1.76$e$-14
	Scheme3	2.73$e$-12	2.14$e$-14	1.48$e$-12	1.94$e$-14	1.50$e$-12	1.82$e$-14

Table 2. Table II: Drag coefficient c d subscript 𝑐 𝑑 c_{d} and added mass coefficient c i subscript 𝑐 𝑖 c_{i} for the flow over an circular cylinder oscillating in the stationary fluid.

Work	$c_d$	$c_i$
Dütsch et al. [34]a	2.09	1.45
Uzunoğlu et al. [37]a	2.10	1.45
Yuan et al. [17]	2.10	1.47
Present	2.10	1.48

Equations

$$\frac{\partial f}{\partial t}+\mathbf{\xi }\cdot \nabla f=\Omega =-\frac{1}{\tau }[f-f^{eq}],$$

$$f^{eq}=\frac{\rho }{(2\pi RT{)}^{(3+K)/2}}exp(-\frac{c^2+{\eta }^2+{\zeta }^2}{2RT}),$$

$$g(\mathbf{x},\mathbf{\xi },t)=\int f(\mathbf{x},\mathbf{\xi },\mathbf{\eta },\mathbf{\zeta },t)d\mathbf{\eta }d\mathbf{\zeta },$$

$$h(\mathbf{x},\mathbf{\xi },t)=\int ({\eta }^2+{\zeta }^2)f(\mathbf{x},\mathbf{\xi },\mathbf{\eta },\mathbf{\zeta },t)d\mathbf{\eta }d\mathbf{\zeta }.$$

$$\frac{\partial g}{\partial t}+\mathbf{\xi }\cdot \nabla g={\Omega }_g=-\frac{1}{\tau }[g-g^{eq}],$$

$$\frac{\partial h}{\partial t}+\mathbf{\xi }\cdot \nabla h={\Omega }_h=-\frac{1}{\tau }[h-h^{eq}],$$

$$g^{eq}=\int f^{eq}d\mathbf{\eta }d\mathbf{\zeta }=\frac{\rho }{(2\pi RT{)}^{D/2}}exp(-\frac{(\mathbf{\xi }-\mathbf{u}{)}^2}{2RT}),$$

$$h^{eq}=(K+3-D)RT{g}^{eq}.$$

$$\rho =\int gd\mathbf{\xi },\rho \mathbf{u}=\int \mathbf{\xi }gd\mathbf{\xi },\rho E=\frac{1}{2}\int ({\xi }^{2}g+h)d\mathbf{\xi },$$

$$\mu ={\mu }_{ref}{\left(\frac{T}{T_{ref}}\right)}^{\omega },$$

$$\frac{\partial \varphi }{\partial t}+\mathbf{\xi }\cdot \nabla \varphi =\Omega =-\frac{1}{\tau }[\varphi -{\varphi }^{eq}],$$

$${\varphi }_j^{n+1}(\mathbf{\xi })-{\varphi }_j^n(\mathbf{\xi })+\frac{△t}{\left|V_j\right|}{\mathbf{F}}^{n+1/2}(\mathbf{\xi })=\frac{△t}{2}[{\Omega }_j^{n+1}(\mathbf{\xi })+{\Omega }_j^n(\mathbf{\xi })],$$

$${\mathbf{F}}^{n+1/2}(\mathbf{\xi })={\int }_{\partial V_j}(\mathbf{\xi }\cdot \mathbf{n})\varphi (\mathbf{x},\mathbf{\xi },t_{n+1/2})dS,$$

$$\tilde{\varphi }=\varphi -\frac{△t}{2}\Omega =\frac{2\tau +△t}{2\tau }\varphi -\frac{△t}{2\tau }{\varphi }^{eq},$$

$${\tilde{\varphi }}^{+}=\varphi +\frac{△t}{2}\Omega =\frac{2\tau -△t}{2\tau +△t}\tilde{\varphi }+\frac{2△t}{2\tau +△t}{\varphi }^{eq}.$$

$${\tilde{\varphi }}_j^{n+1}={\tilde{\varphi }}_j^{+,n}-\frac{△t}{\left|V_j\right|}{\mathbf{F}}^{n+1/2}(\mathbf{\xi }).$$

$$\rho =\int \tilde{g}d\mathbf{\xi },\rho \mathbf{u}=\int \mathbf{\xi }\tilde{g}d\mathbf{\xi },\rho E=\frac{1}{2}\int ({\xi }^2\tilde{g}+\tilde{h})d\mathbf{\xi }.$$

$$\varphi ({\mathbf{x}}_b,\mathbf{\xi },t_n+s)-\varphi ({\mathbf{x}}_b-\mathbf{\xi }s,\mathbf{\xi },t_n)=\frac{s}{2}[\Omega ({\mathbf{x}}_b,\mathbf{\xi },t_n+s)+\Omega ({\mathbf{x}}_b-\mathbf{\xi }s,\mathbf{\xi },t_n)].$$

$$\bar{\varphi }=\varphi -\frac{s}{2}\Omega =\frac{2\tau +s}{2\tau }\varphi -\frac{s}{2\tau }{\varphi }^{eq},$$

$${\bar{\varphi }}^{+}=\varphi +\frac{s}{2}\Omega =\frac{2\tau -s}{2\tau +s}\bar{\varphi }+\frac{2s}{2\tau +s}{\varphi }^{eq}.$$

$$\bar{\varphi }({\mathbf{x}}_b,\mathbf{\xi },t_{n+1/2})={\bar{\varphi }}^{+}({\mathbf{x}}_b-\mathbf{\xi }s,\mathbf{\xi },t_n).$$

$$\varphi ({\mathbf{x}}_b,\mathbf{\xi },t_{n+1/2})=\frac{2\tau }{2\tau +s}\bar{\varphi }({\mathbf{x}}_b,\mathbf{\xi },t_n+s)+\frac{s}{2\tau +s}{\varphi }^{eq}({\mathbf{x}}_b,\mathbf{\xi },t_n+s).$$

$${\bar{\varphi }}^{+}=\frac{2\tau -s}{2\tau +△t}\tilde{\varphi }+\frac{3s}{2\tau +△t}{\varphi }^{eq}.$$

$${\bar{\varphi }}^{+}({\mathbf{x}}_a,t_n)={\bar{\varphi }}^{+}({\mathbf{x}}_j,t_n)+({\mathbf{x}}_a-{\mathbf{x}}_j)\cdot \nabla {\bar{\varphi }}^{+}({\mathbf{x}}_j,t_n),$$

$$\underset{\nabla {\bar{\varphi }}_j^{+}}{\min }\sum _n{w}_{j,n}[{\bar{\varphi }}_{j,n}^{+}-{\bar{\varphi }}_j^{+}-\nabla {\bar{\varphi }}_j^{+}\cdot ({\mathbf{x}}_{j,n}-{\mathbf{x}}_j){]}^2,$$

$${\tilde{\varphi }}^{+}=\frac{4}{3}{\bar{\varphi }}^{+}-\frac{1}{3}\tilde{\varphi }.$$

$$\frac{\partial f}{\partial t}+(\mathbf{\xi }-\mathbf{v})\cdot \nabla f=\Omega =-\frac{1}{\tau }[f-f^{eq}].$$

$${\varphi }_j^{n+1}(\mathbf{\xi })\left|V_j^{n+1,\ast }\right|-{\varphi }_j^n(\mathbf{\xi })\left|V_j^{n,\ast }\right|+△t{\mathbf{F}}_{ALE}^{n+1/2}(\mathbf{\xi })=\frac{△t}{2}[{\Omega }_j^{n+1}(\mathbf{\xi })\left|V_j^{n+1,\ast }\right|+{\Omega }_j^n(\mathbf{\xi })\left|V_j^{n,\ast }\right|],$$

$${\mathbf{F}}_{ALE}^{n+1/2}(\mathbf{\xi })={\int }_{\partial V_j}(\mathbf{\xi }-\mathbf{v})\cdot \mathbf{n}\varphi (\mathbf{x},\mathbf{\xi },t_{n+1/2})dS=\sum _k(\mathbf{\xi }-{\mathbf{v}}_{b,k}^{n+1/2})\cdot {\mathbf{n}}_{b,k}^{\ast }{\varphi }^{n+1/2}({\mathbf{x}}_{b,k},\mathbf{\xi })S_k^{\ast },$$

$${\tilde{\varphi }}_j^{n+1}=\left|\frac{V_j^{n,\ast }}{V_j^{n+1,\ast }}\right|{\tilde{\varphi }}_j^{+,n}-\frac{△t}{\left|V_j^{n+1,\ast }\right|}{\mathbf{F}}_{ALE}^{n+1/2}(\mathbf{\xi })=\left|V_{j,mod}\right|{\tilde{\varphi }}_j^{+,n}-\frac{△t}{\left|V_j^{n+1,\ast }\right|}{\mathbf{F}}_{ALE}^{n+1/2}(\mathbf{\xi }),$$

Full text

1 Introduction

Moving boundary problems can be found in various scientific and engineering fields. For example, in aerospace area, the store separation from the aircraft body, the motion of landing gear during the take-off and landing, etc. [1]. For micro air vehicles, the design of rotary wings and flapping wings which include lots of moving boundary problems will be another application aera [2]. In the above examples, usually the reference lengths of obstacle are much larger than the mean-free-path of gas molecule. According to the definition of Knudsen number [3], the flows in above cases are the continuum flows. Furthermore, the moving boundary problems are also encountered in the application of rarefied flow regime. Such as, the sound wave generated by the oscillating plate, the motion of vanes for the Crookes radiometer, etc. [4]. Generally speaking, the macro-methods based on the N-S equations [5] can cope with the corresponding problems at continuum flow regime, and the prevailing direct simulation Monte Carlo (DSMC) method [6] can deal with that at rarefied flow regime. But in some applications, both macro-methods and DSMC will encounter obstacle as the flow regime out of its computational range, such as the flow in transition regime. Though some hybrid methods [2, 7] can be used, due to the different temporal and spatial scales, this kind of hybrid method also encounters great difficulties. Consequently, developing and improving a method which can simulating the moving boundary problems at all flow regimes will have enormous values for engineering application.

Recently, the discrete unified gas kinetic scheme (DUGKS) proposed by Guo et al. [8] is a promising method which can handle flows at all regimes [9]. For this method, it combines the advantages of both lattice Boltzmann method [10] (LBM) and unified gas kinetic scheme [11] (UGKS), where the flux at cell interface is easy to calculate like LBM, and the computational cost is declined compared with the UGKS. Some details can be found in Ref. [8, 12, 13, 14]. At current stage, the DUGKS is implemented under the stationary mesh, so the purpose of this paper is further extends the application range of original DUGKS, that is can deal with moving boundary problems.

Nowadays, there are many methods can handle the moving boundary problems. Based on the mesh systems used in the numerical computation, these methods simply can be divided into two categories: Eulerian method and Lagrangian method. For Eulerian method, the mesh is fixed at each iteration computation. The immersed boundary method (IBM) is one of the representation [15, 16, 17]. In this method, the uniform Cartesian grids are used near the region of wall at present stage, the moving boundary will be regarded as a set of Lagrangian nodes, and the influence of wall (no-slipping condition) to the nearby nodes of the Cartesian grid is considered with some interpolation methods. Besides the applications in the continuum flow regime, the IBM coupled the UGKS now can also deal with the rarefied flow moving boundary problems [18]. The primary disadvantage of IBM is in some applications, such as high Reynolds number flows, the amount of mesh is intolerable. The static mesh movement method [19] will be another representation of Eulerian method to cope the moving boundary problems. During the numerical simulations, after an Eulerian step, a new mesh is generated according to certain requirement. In general, regenerating a new mesh is time-consuming for most applications. Besides. some interpolation methods also needed to transfer the flow variables from the old to the new mesh. For Lagrangian method, the moving velocity of mesh is equal to the local fluid flow velocity, so the mesh distortion and tangling is unavoidable for most cases.

To combine the advantages of both Eulerian and Lagrangian method, a famous method, that is arbitrary Lagrangian-Eulerian [20] (ALE) technique is developed and improved during the last few decades. Nowadays, the ALE method also can be divided into two types. In traditional procedure, three steps, that are the explicit Lagrangian phase, the rezoning phase, and the remapping phase, will be implemented. Similar to the static mesh movement method in pure Eulerian method, the mesh regeneration or modification, and flows variables transfer will be the two critical steps for this ALE type. In other words, if the quality of mesh can maintain very well during the moving process, this type ALE will degenerated into the pure Lagrangian method. For another type of ALE, the mesh velocity will be introduced into the governing equations (convective terms) to modify the net flux of cell interface. And the mesh moving technique will be introduced to maintain the mesh quality. Besides, the mesh moving velocity also will be constructed based on the old and the new meshes. In aerospace area, such as the aeroelastic analysis, this type of ALE method is usually used. Generally speaking, the time-consuming for the implement of mesh moving technique is less than that of regenerating a new mesh, and rezoning and remapping phases can be discarded for this type ALE, so the remapping-free ALE technique will be used in this paper to improve the original DUGKS.

For the methods of dealing with the moving boundary problems, one source of numerical error is violates the geometric conservation law (GCL). In some applications, it will yields erroneous results [21]. Following the previous works, in this paper, GCL compliance schemes are considered to exclude this part of numerical errors.

The rest of this paper is organized as follows. In Sec. 2, original DUGKS is introduced briefly, ALE-type DUGKS and several GCL schemes are detailed illustrated. In Sec. 3, one case to verify the GCL schemes, and four test cases to validate the capacity of present method are conducted. Finally, a short conclusion is summarized in the Sec. 4.

2 Numerical Methods

2.1 The sketch of original discrete unified gas kinetic scheme

In this section, the original DUGKS proposed by Guo et al. [8] is introduced briefly. The starting point of DUGKS is the Boltzmann-BGK equation, which can be expressed as

[TABLE]

where $f=f(\bm{x},\bm{\xi},\bm{\eta},\bm{\zeta},t)$ is the velocity distribution function for particles moving in $D$-dimensional velocity space with $\bm{\xi}=(\xi_{1},\dots,\xi_{D})$ at position $\bm{x}=(x_{1},\dots,x_{D})$ and time $t$. $\bm{\eta}=(\xi_{D+1},\dots,\xi_{3})$ is the rest components of the particle velocity with the length $L=3-D$. $\bm{\zeta}$ is a vector with $K$ dimension which represent the internal degree of freedom of molecules. $\tau$ is the relaxation time relating to the fluid dynamics viscosity $\mu$ and pressure $p$ with $\tau=\mu/p$. And $f^{eq}$ is the Maxwellian equilibrium distribution function, which is given by

[TABLE]

where $R$ is the gas constant, $T$ is the fluid temperature, $\rho$ is the density of fluid, and $\bm{c}=(\bm{\xi}-\bm{u})$ is the peculiar velocity with $\bm{u}$ being the macroscopic flow velocity.

In order to remove the dependence of distribution function on the internal variable $\bm{\eta}$ and $\bm{\zeta}$, usually, two reduced distributions [22] can be introduced in practical computation, and respectively, given by

[TABLE]

and

[TABLE]

With Eq. (1), the evolution equations for $g$ and $h$ can be expressed as

[TABLE]

[TABLE]

respectively, where the equilibrium distribution functions for $g^{eq}$ and $h^{eq}$ are given by

[TABLE]

[TABLE]

The macro-quantities can be calculated with

[TABLE]

and with ideal gas law, $p=\rho{R}T$, the pressure can be obtained. In additional, the relationship between dynamics viscosity $\mu$ and temperature $T$, based on the hard-sphere (HS) or variables hard-sphere (VHS) molecules, can be given by

[TABLE]

where $\omega$ is the index related to the HS or VHS model, $\mu_{ref}$ is viscosity at the reference temperature $T_{ref}$.

As Eq. (5) and Eq. (6) are exactly the same, which can be rewritten as

[TABLE]

with $\phi$ represent $g$ or $h$. In DUGKS, Eq. (11) will be solved with finite volume method. And the discretization of this equation can be divided into the two steps: velocity-space discretization and physical-space discretization.

For particle velocity-space discretization, usually, a finite set of discretized micro-velocities will be used [8], and $\bm{\xi}_{i}$ represents the $i$-th discretized velocity. As the macro-quantities depending on the particle micro-velocities (Eq. (9)), the choose of the values of micro-velocity can be set coincide with the abscissas of quadrature rule. For low-speed continuum flows, the discretized velocities, weights, and corresponding equilibrium distribution functions developed in LBM, such as D2Q9 model [23], can be used in DUGKS, which also build the connection between LBM and DUGKS. For rarefied flows, usually, the Gauss-Hermit and Newton-Cotes quadrature rules will be used to integral the macro-quantities, and corresponding set of abscissas will be used as the set of discretized velocities.

For physical-space discretization, in this paper, the finite volume method based on unstructured mesh is used. Fig. 1 shows the schematic of unstructured mesh, $j$ is the center of triangle cell $ABC$, and subscript represent the index number of cell. If $\phi_{j}$ and $\Omega_{j}$ are the average values of $\phi$ and $\Omega$ in cell $ABC$, $\bigtriangleup{t}=t^{n+1}-t^{n}$ is the time step, and the mid-point rule is used for the time integration of the convection term and the trapezoidal rule for the collision term, Eq. (11) can be rewritten as

[TABLE]

where $V_{j}$ represent the volume of cell $ABC$. The $\bm{F}^{n+1/2}(\bm{\xi})$ is the flux of cell surface and given by

[TABLE]

where $\partial{V_{j}}$ is the cell surface, $\bm{x}$ is the center of cell interface, and $\bm{n}$ is the outward unit vector normal to the surface. To remove the implicit collision term, two new distribution functions will be introduced:

[TABLE]

[TABLE]

Then Eq. (11) will be further rewritten as

[TABLE]

Due to the conservative property of the collision term, in practical computation, $\tilde{\phi}$ will be solved instead of the $\phi$. The Eq. (9) for computing the macro-quantities also will be rewritten as

[TABLE]

For the calculation of interface flux, Fig. 1 shows the schematic of cell interface. $\bm{x}_{b}$ is the center of $BC$, and $\bm{n}_{b}$ is unit vector normal of $BC$ with the direction pointing from cell $ABC$ to $BEC$. Integral the Eq. (11) along the characteristic line within a half time step $s=\bigtriangleup{t}/2$, the original distribution function $\phi(\bm{x}_{b},\bm{\xi},t_{n}+s)$ in Eq. (13) can be updated with

[TABLE]

Similar the treatment to $\tilde{\phi}$, another two new distribution functions will be introduced and given by

[TABLE]

[TABLE]

Then Eq. (18) can be rewritten as

[TABLE]

If we instead the $\tilde{\phi}$ in Eq. (17) by $\bar{\phi}$, the macro-quantities at cell interface can be calculated. And with Eq. (19), the original distribution function at $(\bm{x}_{b},t_{n+1/2})$ is given by

[TABLE]

As illustrated by Eq. (21) and Fig. 1, for given $\bar{\phi}^{+}$ at $\bm{x}_{a}=\bm{x}_{b}-\bm{\xi}s$, the $\phi^{n+1/2}$ at cell interface can be obtained. With Eq. (14), (19) and (20), we can build the relationship between $\tilde{\phi}$ and $\bar{\phi}^{+}$:

[TABLE]

So, with $\tilde{\phi}$ and corresponding equilibrium distribution function $\phi^{eq}$ stored at $\bm{x}_{j}$. According to Taylor expansion, we can calculate the $\phi^{+}(\bm{x}_{a},t_{n})$ with

[TABLE]

where $\nabla\bar{\phi}^{+}$ is the gradient of $\bar{\phi}^{+}$. For the calculation of $\nabla\bar{\phi}^{+}$ under the unstructured mesh, in this paper, the least square method will be used:

[TABLE]

where $w_{j,n}=1/(\bm{x}_{j,n}-\bm{x}_{j})^{2}$ is the geometrical weighting factor, $n$ is the total number of cell neighbor. Besides, in practical computation, when $\bar{\phi}^{+}$ has been calculated, with Eq. (15) and (23), the $\tilde{\phi}^{+}$ in Eq. (16) can be updated with

[TABLE]

2.2 The ALE-type discrete unified gas kinetic scheme

2.2.1 The discretized method for ALE-type DUGKS

In this section, the discretized formulation of ALE-type DUGKS will be introduced detailed. Under the ALE framework, during the simulation, the geometry information of the cell, such as volume, the location of cell center, the length of cell interface, etc., will be changed from time to time. Following the handling method used in macro numerical method based on the N-S equations and UGKS which also based on the Boltzmann-BGK equation [24], the mesh moving velocity $\bm{v}$ that modify the net flux of cell interface is introduced, and Eq. (1) is rewritten as

[TABLE]

The discretization scheme that presented in original DUGKS will also be used in ALE-type DUGKS, then Eq. (12) and (13) also will be modified as

[TABLE]

and

[TABLE]

respectively. Where $V^{n+1,*}$ and $V^{n,*}$ are the cell volumes at $n+1$ and $n$ time steps, superscript ∗ means that the value of volume at corresponding time maybe not equal to the true value of volume at that time and will be illustrated at next section. $k$ is total number of cell interface. $\bm{v}_{b}^{n+1/2}$ is the moving velocity of cell interface at $n+1/2$ time. $\bm{n}_{b}^{*}$ and $S_{b}^{*}$ are the outward unit normal vector and area of cell interface, respectively. The computational method for these three variables also will be illustrated at next section.

Similar to the original DUGKS, in ALE-type DUGKS, two new distribution functions, $\tilde{\phi}$ and $\tilde{\phi}^{+}$ also will be introduced, and Eq. (28) can be modified as

[TABLE]

where the formulations of $\tilde{\phi}_{j}^{n+1}$ and $\tilde{\phi}_{j}^{+,n}$ are same to that of original DUGKS. The scheme for reconstructing the $\phi^{n+1/2}$ at cell interface also same to that of original DUGKS, and $\bm{x}_{a}=\bm{x}_{b}^{n}-\bm{\xi}s$ is used to compute the location of interpolating point. Besides, the sign of $\bm{\xi}\bm{n}^{n}$ is used to judge the upwind direction. Finally, from the perspective of programming, if three modifications are added into the framework of original DUGKS solver, that are

(1)

Add the independent module of mesh moving and geometry information updated part into the old solver,

(2)

Modify the $\tilde{\phi}_{j}^{+,n}$ in Eq. (16) with volume modified coefficient $\left|V_{mod}\right|$,

(3)

Modify the vector product in flux computation formulation from $\bm{\xi\cdot\bm{n}}$ to $(\bm{\xi}-\bm{v})\cdot\bm{n}$.

It will become the ALE-type DUGKS solver, and has ability to deal with both continuum and rarefied flow problems with moving boundary.

Additionally, Laplace smoothing equations for mesh deformation [25] is solved to update the unstructured mesh under moving boundary.

2.2.2 Geometric Conservation Law

The concept of GCL was firstly well-defined by Thomas and Lambard in 1979 [26]. Generally speaking, for uniform flow, if scheme based on the moving mesh is GCL complaint, any distribution must not introduced into flow domain at any time. ‘Free Stream Preservation Property’ is the fundamental condition for any time-integral schemes on moving mesh [21]. Following the theoretical derivation used in macro numerical method, the Boltzmann-BGK equation is considered here. If integral the Eq. (1) on the control volume, and semi-discrete form is used, we have

[TABLE]

If flow is uniform, that is $f=const$ in each cell, and with $\int\bm{n}_{b}dS=0$, Eq. (31) can be simplified as

[TABLE]

which is the governing equation of GCL, and means that the variation of volume of a moving cell equals to the integration of the volume-flux (or “sweeping volume”) of all the surfaces (named as “face” in the following context) surrounding the control volume [21].

In this paper, the mesh velocity $\bm{v}_{b}$ of cell interface in Eq. (29) is given as

[TABLE]

Based on the mesh velocity and GCL governing equation, three discretized geometric conservation law (DGCL) compliant schemes which decide the $V_{j}^{n+1,*}$, $V_{j}^{n,*}$, $\bm{S}_{b}^{*}$ will be discussed.

(1) DGCL scheme1:

If we set:

[TABLE]

that are the true values of volume at corresponding times will be used. And following the idea used in Ref. [27], the $\bm{S}_{b}^{*}$ can be calculated with

[TABLE]

Through the simple mathematical derivation, Eq. (34) and (35) will automatically satisfy the Eq. (32) during the mesh moving process.

(2) DGCL scheme2:

If we set:

[TABLE]

to satisfy the DGCL, the volume $V^{*,n+1}$ must be modified [21]. With the first-order Euler time discretized scheme, Eq. (32) can be rewritten as

[TABLE]

and $V_{j}^{n+1,*}$ can be modified with

[TABLE]

(3) DGCL scheme3:

If we set:

[TABLE]

similar to the DGCL scheme2, the $V_{j}^{n}$ must modified and can be calculated with

[TABLE]

As the DGCL scheme1 is much easy to implement, and only one time level ($\bm{F}^{n+1/2}$) is needed to calculate the flux at cell interface, it can natural couple with current ALE-type DUGKS. Though DGCL scheme2 and scheme3 are much complex, further improved ALE-type DUGKS framework such as high-order DUGKS [28] or multi-time-level implicit method [21] which the geometry information at intermediate time-levels are much difficult to define [29], these two schemes are the good choice. Besides, the scheme2 and scheme3 are the volume-constrained scheme, the face-constrained scheme [21] also can be used but not considered in this paper.

2.2.3 Boundary conditions

In this section, the boundary conditions used in this paper will be illustrated in detail.

For the wall boundary condition, depended on the flow conditions, that are continuum or rarefied flows, the non-equilibrium extrapolation [30] or diffuse-scattering rule [8] will be used. For continuum flow, the original distribution functions at $n+1/2$ can be given by

[TABLE]

where subscript $w$ represent the wall boundary, $j$ is the neighbor cell of wall interface, $\rho_{w}$ and $\bm{u}_{w}$ are the density and velocity at wall, respectively. For rarefied flow, the particles which the direction reflect from the wall can be calculated as

[TABLE]

where $\bm{n}_{b}$ represent the unit vector with the direction normal to the wall pointing to the cell. And the density at wall $\rho_{w}$ is determined by the condition that no particles can go through the wall:

[TABLE]

then

[TABLE]

where the distribution functions $f_{w}^{n+1/2}$ with direction $\bm{\xi}_{i}\cdot{\bm{n}_{w}}<0$ can be constructed following the procedure described in Sec.2.2.1.

For the test cases of flow around the obstacle, in this paper, the far-field boundary condition will be used. Similar the treatment to the wall boundary, the directions reflect from the boundary can be calculated as

[TABLE]

where the $\rho_{0}$ and $\bm{u}_{0}$ are the density and velocity of free-stream, respectively.

3 Numerical results and discussions

In this section, several test cases are set up to validate the proposed ALE-type DUGKS in this paper. The first case is the GCL compliance test and shows why its important the GCL be implemented. The second and third cases are low-speed continuum flows around the oscillating circular cylinder and pitching NACA0012 airfoil, respectively, which show the method presented in this paper also can get goods results compared with macroscopic method. The fourth and last cases are low-speed rarefied flows. One is the moving piston driven by the rarefied gas. Another is the flow caused by a plate oscillating in its normal direction. This case is a typical problem in MEMS devices, but currently not finish the systematic research [31]. For continuum cases, only $g$ distribution function governing equation is solved, and three-points Gauss-Hermite quadrature [8] is used to calculate the macro-quantities. For rarefied cases, both $g$ and $h$ distribution functions will be solved, and the quadrature rules are the Gauss-Hermit and Newton-Cotes formulations [14]. Besides, the continuum flow problems are 2D flows and rarefied flow test cases are set as 1D flows.

ALE-type DUGKS formulation has been coded with the help of Code Saturne [32], an open-source computational fluid dynamics software of Electricite De France (EDF), France (http://code-saturne.org/cms/). We appreciate the development team of Code Saturne for their great works.

3.1 The uniform flow for GCL compliance test

Uniform flow usually used as the basic case for GCL compliance test, some additional test cases for GCL complaint scheme can be found in Ref. [21]. For this problem, the computational domain is a square with the size equal to 20, and the number of cells is 40$\times$40. The time steps $\bigtriangleup$t$=0.2$. The initial velocity $u=U_{0}=0.1$, $v=0$, and the initial density $\rho_{0}=1.0$. Furthermore, the far-field boundary condition presented in Sec. 2.2.3 will be used in this case.

Firstly, the important of GCL will be illustrated. For the mode of mesh moving, the grid nodes will be oscillating randomly at its original positions with amplitude of \pm 0.5\bigtriangleup$$x ($\bigtriangleup x$ is the size of a cell). Fig. 2 shows the mesh, respectively, at initial time and instantaneous time when grid nodes start moving. It has been demonstrated that in Sec. 2.2.2, when update the distribution function, if Eq. (34) and (35) are used, the DGCL will be satisfied automatically. Otherwise, if Eq. (36) or (39) are used, and without volume-constrained, the DGCL will be violated. Fig. 3 shows the pressure and velocity contours at 10 time iterations when DGCL violated scheme be implemented. It is clear that the pressure and velocity have been polluted seriously, and simulation can not be continued. So, if the grid nodes moving with some extreme ways, the GCL error is much larger and must be eliminated.

Secondly, the performance of DGCL schemes used in ALE-type DUGKS will be checked. In this part, the grid nodes are also oscillating randomly at its original positions, and spatial errors for pressure and velocity will be calculated at four different mesh sizes (the numbers of cells are equal to $20\times 20$, $40\times 40$, $80\times 80$, and $160\times 160$, respectively). In Table I, after 1000 iterations, the $L_{2}$ and $L_{\infty}$ errors for all schemes maintain the significantly small values, it means that all the DGCL compliance schemes can eliminate the GCL error very well.

As the amplitude of random oscillating for previous numerical tests is little larger, the mesh quantity during the mesh moving is much terrible, so the GCL error is greatly larger. From our further tests, if this amplitude drop to 0.1\bigtriangleup$$x, as the mesh quality improved, the GCL error of DGCL violated scheme much declined. So, if the modes of mesh moving are quite regular, the GCL error will be a small value compared with other numerical errors. But, in some areas like aero-elasticity, it has been reported that the GCL error will yields erroneous results [33]. So, following the suggests proposed in Ref. [21], the DGCL compliance scheme of ALE-type DUGKS will be ulteriorly studied in the further works. Besides, the DGCL scheme1 will be used in this paper as it is easy to implement. The DGCL scheme2 and scheme3 maybe useful to develop other ALE-type DUGKS like high-order DUGKS [28] or multi-steps implicit DUGKS, where the geometry information of cell at middle time level are not easy to defined.

3.2 Continuum flow around an oscillating circular cylinder

In this section, the continuum flow around the oscillating circular cylinder is simulated. For this case, the flow is incompressible, the cylinder oscillates sinusoidally in the horizontal direction, and the equation of motion can be expressed as

[TABLE]

where $x$ is the displacement of cylinder at horizontal direction, $A$ is the amplitude and $f$ is oscillating frequency. Following the set up used in Ref. [34], two key parameters will be defined, that are the Reynolds number, Re, and the Keulegan-Carpenter number, KC, respectively. These two parameters will dominate the pattern of this oscillating flow, and the expressions are

[TABLE]

[TABLE]

respectively, where $d$ is diameter of the cylinder, $U_{max}$ is maximum velocity of cylinder in horizontal direction, $\rho$ is the fluid density, and $\mu$ is the fluid viscosity. In this simulation, the Re = 100, and KC = 5 will be considered.

Fig. 4 shows the mesh used in this case. The size of computational domain is $[80d\times 60d]$, it is very large so as to eliminate the influence of the far-field boundary condition. The total number of cells in the domain is about 49000, with 240 points along the cylinder surface. Quadrangular cells will be used to discretize the region near the cylinder surface with width equal to $0.1d$, and others are the triangles. At initial time, the cylinder is located at the center of computational domain, and $\rho=1.0$, $u=v=0.0$ will be used to initialize the distribution functions. $\rho=1.0$, $u=v=0.0$ also used as far-field boundary conditions.

Firstly, the numerical convergency of inline force $F_{x}$ at different Mach number, Ma, time step, \bigtriangleup$$t, and minimum grid size, \bigtriangleup$$x are tested. For the test based on the Ma, three values, which equal to 0.173 ($U_{max}=0.1$), 0.0866 ($U_{max}=0.05$), and 0.0173 ($U_{max}=0.01$) are considered (The sound speed equal to $1/\sqrt{3}$ in this case), respectively. As shown in Fig. 5, for one oscillating period $T$, it is clear that the inline force $F_{x}$ will much influenced by the different Ma, especially for the amplitude of $F_{x}$. Because in framework of LBM, the equilibrium distribution function used in this case will recover the compressible Navier-Stokes equations, so the compressible effect might lead to some undesirable errors in numerical simulations [35]. For the tests based on the $\bigtriangleup{x}$ and the $\bigtriangleup{t}$, as illustrated in Fig. 6 and Fig. 7, respectively, the influences to result are much declined compared with that of Ma. Consider the computational time, Ma=0.0866 ($U_{max}=0.05$), $\bigtriangleup{x}=256/d$ and $\bigtriangleup{t}=0.05$ will be used in the following simulations.

Secondly, some other results are compared to further validate our computation quantitatively. Pressure isolines at phase position $0^{\circ}$ and $288^{\circ}$ are shown in Fig. 8. These isolines are agree well with the numerical results of Dütsch et al. [34]. Furthermore, the velocity profiles $\bar{u}_{x}$ and $\bar{u}_{y}$ at four vertical cross section $\bar{x}=-0.6,0,0.6,1.2$ for three phase positions $180^{\circ}$, $210^{\circ}$ and $330^{\circ}$ are also compared with numerical and experiment results obtained by Dütsch et al. [34], where $\bar{u}_{x}$, $\bar{u}_{y}$, and $\bar{x}$ ($\bar{y}$ is the vertical distance for comparison) are defined as

[TABLE]

$x$ and $y$ are the coordinates relative to the equilibrium position of cylinder, $u_{x}$ and $u_{y}$ are the velocity components in horizontal and vertical direction. As shown in Fig. 9, generally, our results also agree well with the numerical and experiment results of Dütsch et al. [34].

For oscillating flow, the semi-empirical equation of Morison et al. [36] are widely used to estimate the inline force $F_{x}$ on body. When the circular cylinder oscillating in the stationary fluids, the time-dependent inline force $F_{x}$ is expressed as

[TABLE]

where $x$ is the displacement of cylinder, $c_{d}$ and $c_{i}$ are the drag coefficient and the added mass coefficient, respectively. Integral the pressure and stress at the surface of cylinder, the $F_{x}$ can be calculated, then with the help of least-squares fitting or Fourier analysis, the $c_{d}$ and $c_{i}$ also can be evaluated. Fig. 10 shows the time history of $F_{x}$, it is clear that the pressure is dominating contribution to the total force. Similar behavior also described by Dütsch et al. [34]. The fitted $c_{d}$ and $c_{i}$, and some other numerical results are compared in Table II, though the $c_{i}$ is little higher, generally the present ALE-DUGKS can get good results. Given the $c_{d}$ and $c_{i}$, Eq. (50) can evaluate the empirical values of $F_{x}$ at one oscillating period. In Fig. 11, our results are great well with that of Dütsch et al. [34], and rougher consistency with the empirical values of Morison et al. [36].

3.3 Continuum flow around a pitching NACA0012 airfoil

When study the details of propulsion for insects, birds, or fishes, the flows around oscillating airfoil with pitching and/or heaving motions usually used as benchmark test cases [38]. In this section, only motion of pitching is considered, which also can check the performance of present ALE-DUGKS to cope with the rotary moving boundary problem. The profile of airfoil is NACA0012, and pitching at its quarter-chord. The variation of attack of angle (AoA) $\alpha$ for pitching airfoil can be expressed as

[TABLE]

where $d$ is amplitude of pitching, and $f$ is pitching frequency. Usually, for pitching airfoil, a new parameter, that is reduced frequency, $k$, will be defined. The relationship between $k$ and $f$ is $k=2\pi{f}c/U_{0}$, where $c$ is the length of airfoil chord ($c=1$), and $U_{0}$ is the velocity of free-stream. The Reynolds number, Re, based on the $c$, is set as 12000. Generally, flow under this Re can be treated as laminar flow [39]. Fig. 12 shows the mesh used in this case. The total number of cells in the domain is $140,831$, along with 291 points at the airfoil surface. Region near the airfoil is discretized into rectangle cells, and the minimum size of cell is $1.0\times 10^{-4}$. About 15$c$ length away from the airfoil trail, Cartesian grids are used to capture the vortical patterns. The number of cells is little larger in this case, if only consider some results like force coefficients, about $50\sim 60\times 10^{3}$ cells are enough to get good results. In addition, as shown in Sec. 3.2, a small value of $U_{0}$ ($U_{0}=0.05$) is used to eliminate the compressible effect.

Firstly, the flow around stationary airfoil is simulated. Reported by Koochesfahani [40] with experiment, under this Reynolds number, the phenomenon of vortex-shedding will generate, and the equivalent reduced frequency, $k_{equi.}$, based on the $f$, is 8.7. Fig. 13 shows the time evolutions of lift coefficient $C_{l}$ and drag coefficient $C_{d}$, where $C_{l}$ and $C_{d}$ are defined as

[TABLE]

$F_{x}$ and $F_{y}$ are the components of force at $x$-direction and $y$-direction, and $\rho_{0}$ is density of free-stream ($\rho_{0}=1.0$ in this case). As the amplitude of $C_{l}$ is very small, and maybe due to the unstructured mesh is used, the maximum and minimum of $C_{l}$ is little different. From our test, the numerical result of $k_{equi.}$ is equal to 8.23, close to the experiment value. Besides, the mesh at the region about $0.5c$ width behind the airfoil trail must be refined. If the mesh at this region is too coarse, such as O-type mesh, the simulation will result a steady flow.

Secondly, a series of flow over the pitching airfoil are simulated at $d=2^{\circ}$ and $4^{\circ}$, respectively, with different reduced frequency $k$. Fig. 14 shows the time evolutions of $C_{d}$ at three values of $k$. The $C_{d}$ are also calculated with Eq. (52). At small value of $k$, both the maximum and the minimum of $C_{l}$ are larger than zero, so the flow will generate the drag force. When enhance the value of $k$, the minimum of $C_{d}$ will little than zero. Though the flow is still generate the drag force, the magnitude is decreasing. If we continue enhance the value of $k$, finally, the absolute value of minimum of $C_{d}$ will larger than that of maximum of $C_{d}$, and flow will generate the thrust force. Here, we define a new force coefficient, $C_{T}$, represent the thrust coefficient, and $C_{T}=-C_{d}$. Fig. 15 shows the mean thrust coefficient $\bar{C}_{T}$ at different $k\alpha$. Generally, we get the same tendency compared with other numerical results [41, 38, 42] and experiment result [40]. For $d=2^{\circ}$, when $k\alpha$ less than 0.2, the numerical results are consistency with each others, and close to the experiment values. But when $k\alpha$ larger than 0.2, the discrepancy becomes obvious. Some reasons described in Ref. [41] maybe explain this discrepancy. Fig. 16 shows the $\bar{C}_{T}$ at different Ma, the discrepancy at high values of $k\alpha$ is much larger. We demonstrate again that if compressible flow solver is used, the Ma of free-stream must set with a small value, at least for this pitching airfoil test case. The similar profiles also reported by Young et al. [38]. Fig. 17 and Fig. 18 show the time evolutions of $C_{l}$ and $C_{d}$ at two different reduced frequencies. We reproduce the same phenomenons illustrated by Liang et al. [39] at these two flow conditions. That is, for lift coefficient, though the amplitudes are very larger, the zero mean lift will acting on the airfoil. Stress force almost equal to zero, only the pressure is dominating contribution to the total lift force. For drag coefficient, the stress force almost keep the same values at these two flow conditions, but the instantaneous pressure declined, and will offset the stress force, consequently the magnitude of total drag force is reduced.

Fig. 19 and Fig. 20 show the vortical patterns at $d=2^{\circ}$ and $d=4^{\circ}$, respectively. For comparison, the vortical pattern at $k=0$ (stationary airfoil) is also presented. For $d=2^{\circ}$, at small value of $k$, the much larger structures of vortex street will generate compared with that of stationary airfoil. When enhance the value of $k$, as the thrust force is generated, original two rows of vortex street seem merge into one row. Besides, if $k<k_{equi.}$, even though about 15$c$ length far away from the airfoil trail, the clear structures of vortex street are kept very well. Otherwise, these structures are quickly dissipated at high values of $k$. For $d=4^{\circ}$, when $k=0.835$, the wake assumes a form of undulating vortex sheet. And with $k=3.09$, double-vortex pattern is captured. But about $2c$ length behind the airfoil trail, this structure seems to dissipate and merge into one another. For $k=2$, that the value between 0.835 and 3.09, the wake pattern seems that the length of undulating vortex sheet is declined (about $5c$ behind the airfoil trail), and the double-vortex features also generate but not very clear. At higher values of $k$, double-vortex features are disappearing, and the vortex streets that approximate a straight line are deflected. Similar structures of vortex street also illustrated by Liang et al. [39] with numerical method and by Koochesfahani [40] with experiment.

3.4 A moving piston driven by the rarefied gas

In this section, a case that the rarefied gas driving a piston is simulated. This problem has been studied by Dechristé et al. [43] with a deterministic numerical scheme coupled with immersed boundary method, and Shrestha et al. [6] with DSMC. Fig. 21 shows the schematics of this problem. The one-dimensional computational domain will divided into two sub-domains by a piston. The length of one sub-domain is $L$, and the width of piston is $2l$. At initial time, these two sub-domains will be filled with the same gas, and the density $\rho$, pressure $p$, and temperature $T$ also set the same values. For right part of computational domain, as the temperature of wall higher than that of gas, the pressure will enhance, and push the piston moving from right to left. Finally, if the pressures at two walls of piston are same, the piston will stop moving. With the mass conservation and the state equation of gas for each part, we have [43]:

[TABLE]

and

[TABLE]

where $R$ is the gas constant. So, the equilibrium location of piston is

[TABLE]

and, respectively, the density and pressure at equilibrium states for each part are

[TABLE]

Following the set up described by Shrestha et al. [6], the gas is argon, the mass of atom is $m_{g}=6.63\times 10^{-26}kg$, and the diameter of atom is $3.68\times 10^{-10}m$ with the hard-sphere collision model. The initial values for $T_{0}$, $T_{w}$ and $p_{0}$ are equal to $270K$, $330K$, and $10Pa$, respectively. And the initial velocity for each part are set as zero. With ideal gas law, the initial density can be calculated. In our simulations, two computational geometries are considered, that are $L=0.1m$, $l=0.01m$, and $L=0.01m$, $l=0.001m$, respectively. Consequently, based on the width of piston ($2l$) and the given condition of gas, the flows under two Knudsen (Kn) numbers, $Kn=0.031$ and $Kn=0.31$ will be simulated. Besides, 400 cells are used for the total computational domain.

Fig. 22 shows the time evolutions of the piston location at two $Kn$ numbers. The similar profiles also described in Ref. [6]. That is the position of piston for $Kn=0.31$ will much faster converge to its equilibrium position than that for $Kn=0.031$. In addition, DSMC at small value of $Kn$, the results will fluctuating and need perform several independent runs to reduce this stochastic fluctuations in time-dependent problems. On the contrary, our results are much smooth even at small $Kn$. Besides, the red dash and dot lines shown in figure are the theory solutions (Eq. (55)), the errors between the numerical results and theory solutions are both less than $1\%$. Fig. 23 and Fig. 24 show the time evolutions of density and pressure at two sub-domains, respectively. Our numerical results are also consistency with the theory solutions. During the evolutions, the pressure difference of two walls for piston at $Kn=0.31$ are larger than that of $Kn=0.031$, it may also indicate that the piston at computation condition of large $Kn$ will faster moving to its converged equilibrium position.

In this case, the Gauss-Hermit quadrature rule is used [14] to integral the macro-quantities, and the codes to calculate the abscissas and weights are shown in Ref. [44]. Fig. 25 shows the influence of number of quadrature points to the results. For $Kn=0.031$, from our test, the integral accuracy with 28 quadrature points is good enough for left sub-domain to calculate the macro-quantities. But, for the right side, about two times number of point is needed to get good results. Due to the error of integral, the density in right part will continue decline, and the velocity will not converge to zero. Consequently, even though the pressure difference of two walls for piston converge to a very small value, the piston will move from left to right with a small value of velocity, finally the simulation will blow up.

3.5 Rarefied flow caused by a plate oscillating in its normal direction

In this section, another rarefied flow test case will be simulated, which has been studies by Tsuji et al. [31]. Fig. 26 shows the schematic of this problem. In one-dimensional domain, the right wall is stationary, and left wall is oscillating with cosine function $x(t)=a_{w}cos(\omega{t})$ ($a_{w}=0.1$ and $\omega=1$ in this case), the amplitude equal to 0.1 make sure this case is low-speed flow. Following the set up described in Ref. [31], in initial time, the length of computational domain is $d=2\pi\sqrt{5/6}$, which is the wavelength of the sinusoidal acoustic wave with angular frequency $\omega$ in an inviscid (Euler) gas [45]. The Knudsen number is given as

[TABLE]

where $K$ is the special $Kn$ number, and two $K$ numbers, which equal to 0.5 and 1.0, respectively, will be considered in this case. $\rho_{0}=1.0$, $T_{0}=1.0$ and $u_{0}=0.0$ will be used to initial the distribution function $g$ and $h$. 400 cells is used to discretize the computational domain for $K=0.5$, and that of 200 for $K=1.0$. For time step, in each oscillating period ($T=2\pi$), 3,2000 iteration computations will be implemented for $K=0.5$, and 6,4000 for $K=1.0$.

Figs. 27 - 29 show the profiles of density, velocity and temperature, respectively, at ten moments in each oscillating period. Generally, our results are agree well with that of Tsuji et al. [31], especially at high value of $K$. Without more reference, it s hard to identify which result is better. From the figures, it is clear that at $t/\pi=1.0$, the velocity profile is close to the sinusoidal shape, but the profiles of density and temperature deviate from it significantly. Furthermore, for velocity profile, this shape will deviates more from the sinusoidal shape and tends to attenuate more rapidly as $K$ increased, especially for the right part of wave.

For the quadrature rule to calculate macro-quantities, the Newton-Cotes rule will be used. For higher value of $K$, due to the number of particles is small, the wave generated from the left moving wall and the reflected wave from the right stationary wall will lead to the singular of distribution function [4]. So, to get the smooth results, the number of abscissas is much larger for higher $K$. From our test, for $K=0.5$, 50 abscissas is enough to get the convergent and smooth results. But for $K=1.0$, about 200 abscissas can get the smooth results.

4 Conclusion

In the present work, the original DUGKS is extend to ALE-type DUGKS. The mesh moving velocity is introduced into the Boltzmann-BGK equation to modify the net flux of cell interface. Consequently, based on the constructed mesh moving velocity, the remapping-free-type ALE method is used to develop the current ALE-type DUGKS. To exclude the GCL error, three DGCL compliance schemes are discussed. As the present DUGKS only the middle time is needed to calculate the flux of cell interface, based on the geometry average, the geometry information at this time level is easy to defined, so the DGCL scheme1 is a good choice. Further improved work such as high-order or multi-time-levels implicit ALE-type DUGKS which the geometry information at these time levels are not easy to define, DGCL scheme2 and scheme3 will be the good choice. From the uniform flow test case, all these three schemes have good performance which no disturbances will be introduced into the computational domain. Four test cases of low-speed flow are simulated, two of them are the continuum flows, and others are the rarefied flows. Results of all the cases are in good agreement with the other numerical and/or experiment results. Therefore, for continuum flow, similar to the macro-methods based on the N-S equations, the present ALE-type DUGKS has the capability to cope with more complex low-speed moving boundary problems. And for rarefied flows, which out of ability for macro-methods, the present method also has the power to deal with them. Further works, such as parallel computing, implicit accelerated method, etc., will be continued, to enhance the ability of present ALE-type DUGKS for simulating the moving boundary problems at different flow regime.

Acknowledgements

The current work is supported by the National Natural Science Foundation of China (Grant No. 11472219) and the 111 Project of China (B17037).