An ALE formulation for compressible flows based on multi-moment finite volume method

ABSTRACT A direct Arbitrary Lagrangian Eulerian (ALE) formulation using multi-moment finite volume method for viscous compressible flows on unstructured moving grids has been developed and presented in this paper. High-order polynomials are reconstructed over a compact mesh stencil by making use of both volume integrated average value (VIA) and point values (PV) at the vertices of mesh cells which change in time. By formulating the governing equation in ALE integral form and differential form, the computational variables, VIA and PV, are updated, respectively, by a finite volume scheme and a point-wise discretization using multi-moments (VIA and PV). A simple and efficient formulation is derived for moving mesh which satisfies the geometrical conservation law. In the computations involving moving-body, the PVs of velocity at the vertices of cells adjacent to the body surface are coincided with the motion of the body, which eliminates the numerical approximation to find the cell boundary values on the body surface in conventional finite volume method. Numerical tests are presented to verify the present method as a high-order ALE formulation for compressible flows on both stationary and moving grids.


Introduction
Computational fluid dynamics (CFD) has become an important auxiliary tool for engineering applications (Akbarian et al., 2018;Ardabili et al., 2018;Mou, He, Zhao, & Chau, 2017) in the past several decades. Fluid structure interaction (FSI), as an ubiquitous phenomena in nature and engineering applications, has been being an active research area in CFD community. In many numerical solvers for FSI problems, the computational domain moves or undergoes deformations due to the movement of the domain boundary. Thus, numerical methods on moving grids are desired to maintain accuracy and conservation of the numerical solutions for FSI applications which usually result in more complex flow structures.
There are two classical kinematic descriptions in the simulations of compressible fluid dynamics: i.e. the Lagrangian description and the Eulerian description. Lagrangian algorithms inherently have superiority in tracking free surfaces and interfaces between different materials since the computational cells move with local fluid velocity. However, they may encounter difficulties when the mesh cell tangles due to large flow distortions. Eulerian algorithms are more robust for flow simulations as the computational grids keep stationary during whole computational process. However, it is difficult CONTACT Peng Jin jin.p.aa@m.titech.ac.jp in some cases to maintain adequate solution quality on fixed computational mesh. Furthermore, in simulations involving interactions of multi-materials, extra numerical procedures are required to identify the interface between different materials, which might generate numerical errors. Due to the shortcomings of purely Lagrangian and purely Eulerian descriptions, Arbitrary Lagrangian Eulerian (ALE) description is firstly proposed in Hirt, Amsden, and Cook (1974) to solve fluid dynamics in a suitably moving and deforming grid. In ALE formulations, the grids for the computational domain can move arbitrarily and independently of the fluid motions, which makes it more robust in moving boundary simulations. ALE methods have been studied for several decades to handle the difficulties in the simulation of multi-material fluid flows with large distortions (Ahn & Kallinderis, 2006;Bazilevs, Takizawa, & Tezduyar, 2013;Donea, Giuliani, & Halleux, 1982;Donea, Huerta, Ponthot & Rodríguez-Ferran, 2004;Duarte, Gormaz, & Natesan, 2004;Liu, 1988;Lhner, 2008). They can be classified into two classes, i.e. indirect method and direct method. The indirect method generally consists of three steps: (1) a Lagrangian step for updating the grid and the solution; (2) a rezoning step to regularize the heavily distorted or tangled mesh cells; (3) a remapping step transferring the Lagrangian solution to the rezoned mesh which is adjusted in step (2). However, it may be computationally expensive especially for three dimensional realistic calculations due to rezoning and remapping steps. The direct ALE method (Boscheri, Loubére, & Dumbser, 2015;Persson, Bonet, & Peraire, 2009) includes the grid velocity in the numerical formulation, which is consistent with the moving mesh framework. Thus, the remapping procedure is not needed as a separate step. Both direct and indirect methods are widely implemented in CFD.
To construct high-order ALE scheme, existing works of indirect ALE scheme start from the purely Lagrangian framework. Cheng and Shu (2007) proposed a class of Lagrangian type schemes based on high-order essentially non-oscillatory (ENO) reconstruction (Harten, Engquist, Osher, & Chakravarthy, 1987) and obtained third-order accuracy for Euler equations with curved mesh in Cheng and Shu (2008). Dumbser, Uuriintsetseg, and Zanotti (2013) developed a one-dimensional highorder Lagrangian arbitrary high-order accurate (ADER) finite volume method. The discontinuous Galerkin (DG) method has also been implemented in the Lagrangian formulation for solving gas dynamic equations (Vilar, Maire, & Abgrall, 2014). The purely Lagrangian framework shows the strong ability for capturing the contact discontinuities since no advection term is included in the formulation. However, the subsequent rezoning and remapping steps are always needed to adjust the heavily distorted meshes for large distortion problems. For highorder direct ALE schemes, Clair et al. added an edgebased upwind numerical fluxes of ALE part to a conventional cell-centered Lagrangian solver as presented in Clair, Ghidaglia, and Perlat (2016). Boscheri et al. developed high-order direct ALE ADER finite volume formulations in Boscheri and Dumbser (2014) and  with weighted essentially non-oscillatory (WENO) reconstruction and multi-dimensional optimal order detection (MOOD) approach.
As stated above, there are mainly two methodologies to devise high-order ALE schemes. One is based on the conventional finite volume method, which uses wide stencils to generate high-order polynomials. However, choosing the reliable cell stencils on moving grids is not an easy task for distorted or unstructured grids which are commonly used in applications. The other methodology is to use compact stencils by adding the degrees of freedoms (DOF) locally on each cell, such as the DG method (Cockburn, Lin, & Shu, 1989; and the spectral finite volume method (Wang, 2002;Wang & Liu, 2002). Such methods have superior properties for achieving high-order accuracy. On the other hand, for moving boundary problems, especially for fluid-structure interaction applications, the solutions, as well as the coupling conditions, on the fluid-structure interfaces are crucial and require particular attention. Since the aforementioned methods based on increased local DOFs define the computational variables inside each cell, extra numerical steps are needed to determine the values and conditions on the body surfaces or interfaces among different materials.
It is sensible to use the PVs at cell vertices located on the material surfaces as computational variables which are updated at every time step. With these PVs known on the surfaces, the interactions among different materials can be handled in a more accurate and convenient manner. This observation motivates us to re-cast the ALE formulation using both volume integrated average (VIA) and point values (PV) at the cell vertices as the computational variables, i.e. the underlying idea of the multi-moment finite volume method. Multi-moment finite volume method based on local reconstruction was developed by Xie, Ii, Ikebata, and Xiao (2014) and  on unstructured grids to achieve high-order accuracy for incompressible flows in Eulerian representation with fixed grids. In multi-moment finite volume method, both VIA and PVs are used for local polynomial reconstruction, and simultaneously they are treated as computational variables which are solved by integral form and differential form of the governing equations, respectively. Successive works for Euler equations have been implemented in Xie, Deng, Sun, and Xiao (2017) and . Since the high-order polynomial is based on well-determined local cells, it is straightforward to extend previous studies to moving grids. This research focuses on constructing a 2D multi-moment finite volume ALE scheme for viscous compressible flows on moving grids.
We formulate the governing equations in ALE integral form and ALE differential form to update the VIAs and PVs, respectively. The 2D computational domain is partitioned with triangular mesh elements, which facilitates implementations for complex geometries. The computational variables, the VIAs and the PVs at the moving cell vertices, are used to build the local high-order reconstruction consistent with moving mesh, which is similar to the Eulerian framework. The integral form is discretized by the conventional finite volume scheme, where the convective fluxes for VIA can be calculated by commonly used Riemann solvers, such as Rusanov scheme (1962), HLL scheme (Harten, Lax, & van Leer, 1983) and HLLC scheme (Toro, 2013), while the gradients of physical variables on the edge for viscous terms are simply interpolated from the neighboring cell centers. The differential form is discretized by the finite difference scheme in which Roe's Riemann solver (1981) is used for the inviscid flux and the viscous flux term is calculated by T EC (Time-evolution Converting) formulation as shown in Xiao, Akoh, and Ii (2006). The geometrical conservation law (GCL) is automatically satisfied in the present ALE formulation where the mesh cells are moved with the specified velocities at the mesh nodes while the edges of the moving mesh cells remain as straight-line in 2D or planar in 3D. Though linear grid velocity distribution is assumed on cell edge, it does not degrade the accuracy (convergence rate) of the ALE framework as shown later in the numerical tests.
To verify the proposed scheme for moving grids, three kinds of mesh movement strategies are used in this research: (1) mesh moves with the local fluid velocity; (2) mesh moves following an analytical function; (3) mesh moves by RBF (radial basis function) interpolation (Bos, 2010;De Boer, Van der Schoot, & Hester, 2007;Estruch & et, 2013) for moving boundary problems.
The rest of this paper is organized as follows. Section 2 gives the ALE integral form and differential form of the governing equations. Section 3 presents the multimoment reconstruction formulation for triangular elements. The numerical formulation for viscous compressible flow is presented in Section 4. In Section 5, numerical tests are given for verifications. We close the paper in Section 6 by some conclusion remarks.

Governing equations
2D compressible Navier-Stoke equations in Eulerian coordinates can be formulated in the system of conservation laws where U represents the vector of conservative variables, F(U) and G(U, ∇U) represent the vectors of convective flux and viscous flux, respectively. In this formulation, the vector of conservative variables and fluxes take the form as follows with ρ, u, p, M and E being the density, velocity, pressure, momentum and specific total energy of the fluid, respectively. The symbol ⊗ denotes the tensor product andĪ is the unit tensor. Based on Stoke's hypothesis, the stress tensorτ is a linear function of velocity gradient given as The dynamic viscosity μ is represented as a function of temperature T following Sutherland's law The heat flux vector Q is given by where κ = c p μ/Pr is the thermal conductivity with Pr and c p being Prandtl number and the specific heat capacity at constant pressure, respectively. We also include the equation of state to establish the relationship among three thermodynamic variables in the following general form where e = E − 1 2 |u| 2 is the specific internal energy. We consider the ideal gas in the research with the simpler form as p = (γ − 1)ρe. The adiabatic index γ = c p /c v is the ratio of specific heat capacities at constant pressure and constant volume conditions and we set γ = 1.4 in this research unless otherwise stated.

Integral form
Following Mavriplis and Yang (2006), we integrate Equation (1) over a moving control volume (t) which is enclosed by its boundary (t) ≡ ∂ as where n denotes the surface normal vector of the control volume boundary (t). Using the differential identity In this research, we denote the vector of convective flux and viscous flux as F(U) = (F (U) − u g U) · n and G(U, ∇U) = G(U, ∇U) · n, respectively, with u g = x representing the grid moving velocity. Thus, we formulate the ALE integral form of compressible Navier-Stoke equations as In the case of uniform flow, where all fluid variables are uniform and constant, Equation (10) which is the original definition of GCL (Thomas, 1979) in the integral form. Here, V is the control volume.

Differential form
With the relation between spatial time derivatives and referential time derivatives, Equation (1) can be written as where ∂U/∂t| χ denotes the time derivative of conservative variables with respect to grid moving frame. Equation (12) can also be cast in the standard quasi-linear form in 2D as where A = ∂F x (U)/∂x and B = ∂F y (U)/∂y are the Eulerian Jacobian matrices, I is the unit matrix, u g and v g are x and y components of grid velocity u g , respectively. Thus, we formulate the ALE differential form of compressible Navier-Stoke equations as (14) with the ALE Jacobian matrices Here, A and B are the diagonal matrices of original Eulerian Jacobian matrices A and B. R Y and L Y are the right and left eigen matrices of the corresponding Since the grid velocity u g is considered in both integral and differential formulations, the system can be discretized on arbitrarily moving grids. Two common cases are: (1) Eulerian framework with u g = 0 in which mesh keeps stationary; (2) Lagrangian framework with u g = u where the control volume moves with the local fluid velocity.

The grid configuration and notations
In this study, we partition 2D computational domain into N E non-overlapping straight-line triangular elements, denoted by i , i = 1, 2, . . . , N E . Therefore, it is simple to generate computational meshes for arbitrary geometries. And we denote the global nodes (vertices of element cells) by θ k , k = 1, 2, . . . , K. The surrounding cells sharing node θ k are denoted by km (m = 1, 2, . . . M).
To reconstruct the high-order polynomials on local cells, we further denote the three vertices of cell i by θ in (n = 1, 2, 3) with the coordinate (x in , y in ), and its three boundary line segments by i1 = θ i2 θ i3 , i2 = θ i1 θ i3 and i3 = θ i1 θ i2 as shown in Figure 1(a). The length and outward unit normal vector of segment ij , j = 1, 2, 3, are | ij | and n ij = (n xij , n yij ), respectively. The element area and mass center of i are V i and θ ic , respectively.
Two kinds of computational variables, i.e. the VIA over i and the PV at the vertices of cell i are defined as as shown in Figure 1(b), where φ stands for the conservative variables ρ, M, ρE.

Local reconstruction
Following  and Xiao (2014, 2016), the triangular cell i on physical coordinates where We first calculate the first-order derivative of the physical field at cell center in the physical coordinates (φ xic , φ yic ) by using least square method and then transform it into the local coordinates as (φ ξ ic , φ ηic ). With the supplementary of (φ ξ ic , φ ηic ), the reconstruction polynomial φ(ξ , η) can be formulated on the local coordinates as (18) where the coefficients can be determined by the following constraints With above formulations, the first-order derivatives at vertices are firstly computed from the reconstruction functions in local coordinates (φ ξ in , φ ηin ) and then transformed back to the physical coordinates as (φ xin , φ yin ).

Spatial derivatives
With the spatial derivatives of physical field U at global node θ k , (∂U km /∂x, ∂U km /∂y), calculated from the surrounding cells km , the upwind-biased derivatives at point θ k are approximated by (20) The weights here are calculated by where − −− → θ kmc θ k denotes the vector from the mass center θ kmc to vertex θ k , and the unit normal vector n τ ± represents n x± (±1, 0) and n y± (0, ±1), respectively.

Solutions of Navier-Stokes equations for VIA moment
We formulate the semi-discrete form for Equation (10) with a Q-point Gaussian quadrature for surface convective fluxes and a relatively simpler integration for viscous fluxes as where U i denotes the VIA of conservative variables on cell i with volume V i , G q is the Gaussian quadrature point on cell edge ij which separates cells i and ij . Two-points Gaussian quadrature formula is used in this work, which reads )P 2 and ω 1 = ω 2 = 1 2 for the edge with endpoints P 1 and P 2 . Various Riemann solvers for calculating the inviscid fluxF ij on edge ij are implemented in this research, such as the robust Rusanov scheme (1962), HLL scheme (Harten et al., 1983) and the less diffusive HLLC scheme (Luo, Baum, & Löhner, 2004). HLLC scheme is mainly adopted in this research. Rusanov and HLL schemes are employed to problems involving strong shock waves to stabilize the computation.
The viscous numerical fluxes are evaluated based on the values U ij and (∇U) ij of conservative variables on cell edge segment ij . Thus, we formulate the viscous fluxG ij on the edge ij as where U − ij and U + ij are reconstructed line-averaged values of two neighboring cells, respectively, (∇U) ij is linearly interpolated from the spatial derivatives in the cell center of the two neighboring cells.

Solutions of Navier-Stokes equations for PV moment
We solve Equation (14) point-wisely by Roe's Riemann solver (1981) for convective term and T EC (Timeevolution Converting) formulation (Xiao et al., 2006) for viscous term as where the matrices with tilde are calculated from the Roeaveraging values of the surrounding VIAs as and the sound speed is approximated bỹ

Geometrical conservation law
For a uniform flow, the spatial gradients of all flow variables are zero, thus Equation (22) becomes which is the semi-discretization of GCL. Here, u g (G q , t) is the grid velocity at Gaussian point G p on edge ij . For an arbitrarily varying mesh, we can suppose that the grid velocity has a linear distribution over cell edges. Given the grid velocities at two ends of boundary segments ij denoted as u gij0 and u gij1 , Equation (27) can be simplified where the edge-averaged velocity is given as It is noted that, by the assumption that the cell edge moves as a solid body, Equation (28) might not truly realize the pure Lagrangian framework for deformational flows because the cell boundaries do not move with the local fluid in the same accuracy order of the third-order spatial discretizations. As shown in Cheng and Shu (2008), curvilinear meshes are necessary for purely Lagrangian scheme higher than second order. Nevertheless, by assuming cell boundaries always move as a linear edge and then computing the numerical fluxes over this linear edge, it is still possible to achieve thirdorder accuracy for the ALE framework. We call the framework that only moves cell vertices with the local fluid as quasi-Lagrangian framework in this research to distinguish it from the purely Lagrangian framework.

Numerical results
Third-order Runge-Kutta scheme is used for time integration. Since the mesh geometry changes in time, not only VIAs and PVs but also the geometrical quantities, such as volumes, boundary surfaces and vertices of moving cells are updated at each Runge-Kutta sub-step in the simulation.

2D isentropic vortex
To examine the convergence rate of inviscid part of the present scheme, the propagation of an isentropic vortex (Shu, 1998) benchmark test is calculated. we compute this problem on an initial computational domain (0) = [0, 10] × [0, 10] and specify the initial conditions in terms of primitive variables as where δρ, δu, δv, δp are the initial vortex perturbations. We set this perturbation at the center of computation domain (5, 5) with the vortex strength = 5. Then the initial perturbations are given as following, where δT = −((γ − 1) 2 )/8γ π 2 e 1−r 2 denotes the temperature perturbation, and r 2 = (x − 5) 2 + (y − 5) 2 . The initial profile of density is presented in Figure 2 (a). We compute this problem to time t = 1 in the Eulerian and quasi-Lagrangian framework respectively with gradually increased number of triangular elements. The periodic boundary condition is applied. The final vortices of density calculated on 800 triangular elements show in Figure 2. To verify that the mesh can move arbitrarily for the proposed scheme, we also calculate the same case on a specific moving mesh framework where the mesh moves independent of fluid motion. The grid velocities vary with time as where X 0 = 1, Y 0 = 1, n t = 4, n x = 1, n y = 1, t 0 = 2, L x = 10, L y = 10. Figure 3 shows the density and mesh variation with time. Though the mesh moves independent of fluid motion, the vortex propagates in the domain and accurately reproduced during whole computation process.
For this advection test of isentropic vortex, the analytical solution is available, which is just the time-shifted initial profile for the distance of 1 in both x and y directions. Thus we measure the L 1 and L ∞ errors of the Eulerian, quasi-Lagrangian and the specific-speed moving frameworks, and show the results in Table 1. As expected, third-order numerical accuracy is achieved for all of these three frameworks. It also demonstrates that higher than second-order accuracy can be achieved for linear edge element by ALE scheme.

Free stream preservation
Free stream preservation is a common test to check the GCL for moving mesh schemes. The computational domain and mesh are set as the same as in Section 5.1. The initial conditions are uniform for all physical fields with (ρ 0 , u 0 , v 0 , p 0 ) = (1, 1, 1, 1). We compute this test up to time t = 1 with mesh moving following the formulations (32a) and (32b), and plot the L 1 and L ∞ numerical errors of density for meshes with different element sizes as shown in Figure 4. The numerical errors for all mesh resolutions are of machine precision, which justifies the exact geometrical conservation of the proposed method.

Saltzman problem
Saltzman problem is a well-known challenging benchmark test to assess the robustness of numerical algorithms involving mesh movement, which consists of a strong piston-driven shock wave as shown in Boscheri and Dumbser (2013), Maire (2009), Maire, Vchal (2011), andZanotti (2015). It is firstly proposed by Dukowicz and Meltz (1992) for a two-dimensional uniformly skewed Cartesian grid in a computational domain of [0, 1] × [0, 0.1] which is discretized by 100 × 10 quadrilateral elements. The skewed coordinate of grid x = (x , y ) is transformed from the uniform Cartesian grid x = (x, y) by   In this research, each quadrilateral element is split into two triangular elements as shown in Figure 5. The initial condition is given by an ideal gas (γ = 5 3 ) with (ρ 0 , u 0 , v 0 , p 0 ) = (1, 0, 0, 2 3 · 10 −4 ). The gas is pushed by a moving piston from left to right with velocity (1, 0). A moving slip wall boundary condition is set for the piston on the left and reflective boundary condition is set for other boundaries. The analytical solution is a one-dimensional infinite strength shock wave that moves toward right direction at speed 4 3 with a post-shock density plateau of 4. To meet the motion of the left piston, we performed the simulation to time t = 0.6 in a quasi-Lagrangian framework and also in a specific-speed ALE framework where the mesh movement is given as We use Rusanov scheme (1962) for convection flux and multi-dimensional limiting process (MLP) scheme (Park, Yoon, & Kim, 2010;Xie et al., 2017) is employed for suppressing numerical oscillations. The mesh and density map at final time are given in Figure 6 for both frameworks and density values for the whole field are plotted in Figure 7. Since the ideal gas is filled in the closed container, we also measure the total mass of the gas. The mass variation with time is of machine error, which verifies the conservation of the present scheme.

Viscous double Mach reflection
Double Mach reflection is originally proposed by Woodward and Colella (1984) for compressible Euler equations, which involves a Mach 10 shock that hits a 30 o ramp. We run a viscous version of double Mach reflection (Dumbser, Peshkov, Romenski, & Zanotti, 2016) in this section to verify the capability of the present solver for simulating viscous flows. The computational domain is set as [0, 3.2] × [0, 1]. The initial conditions, which can be obtained by Rankine-Hugoniot conditions, are given here as Reflecting wall lies along the bottom of the domain where x > 1 6 . Analytical solution of a shock Mach number M s = 10 is imposed on the upper wall and bottom wall (x < 1 6 ). Inflow and outflow boundary conditions are set for left and right sides, respectively. We partitioned the domain into totally 296,418 triangular elements with 640 and 200 uniform segments in the corresponding boundaries. Constant dynamic viscosity is used in this calculation. The other parameters are: γ = 1.4, c v = 2.5, κ = γ c v μ/Pr with Pr = 3 4 . We first performed the inviscid case up to time t = 0.2 with HLL Riemann solver (Harten et al., 1983) for inviscid finite volume flux, and also MLP scheme (Park et al., 2010;Xie et al., 2017) is used for limiting process. The well resolved vortex structures as shown in Figure 8 demonstrate the high fidelity of present method. For viscous shock waves, the Reynolds number is defined based on the shock speed and an unitary reference Length (L = 1) as Re s = ρ 0 c 0 M s L/μ with c 0 = γ p 0 /ρ 0 . We then calculate the cases with Re s = 1000 and Re s = 100 and the comparisons with inviscid case are shown in Figure 9. We observe that, with a rather larger physical viscosity, the small-scale vortex structures are suppressed, which is comparable with Dumbser et al. (2016).

Flow over cylinder
Flow around a fixed or oscillating cylinder is a widely used benchmark test to evaluate numerical models for     fluid-structure interactions. We firstly compute the flow past a fixed cylinder to test the performance of the present solver for stationary grids, and then simulate the flow past a transversely oscillating cylinder to verify its capability for moving boundary problems.

Fixed cylinder
A fixed cylinder of diameter d = 1 is located at the origin and the computational domain is set as a rectangular of [−10 : 30] × [−10 : 10]. The Reynolds number based on the cylinder diameter is defined as Re = u ∞ d/μ and Mach number based on frees-tream flow is given as Ma ∞ = u ∞ /c ∞ , where u ∞ and c ∞ are the xcomponent of inlet velocity and inlet sound speed, respectively. To obtain the desired Reynolds number Re = 150 and Mach number Ma ∞ = 0.2, we set initial conditions (ρ ∞ , u ∞ , v ∞ , p ∞ ) = (1, 0.2, 0, 1/γ ) and the Prandtl number Pr = 1. The linear viscosity law is chosen in the calculation with β = 2 and s = 0. To avoid the reflection of shock wave from boundaries, the characteristic boundary condition (Carlson, 2011) is used for inlet and outlet boundary. The slip boundary condition is set for top and bottom walls.
The grid dependency test was firstly conducted to choose a reliable mesh for later simulations. The computational domain is partitioned by a mesh of triangular elements, where the area near the cylinder has relatively smaller cell sizes in order to adequately represent the geometric configuration of the solid body as well as the induced wake and vortex structures. Three kinds of grids, Grid A, Grid B and Grid C, are generated where distances between neighboring points on the cylinder are around 1 20 d, 1 40 d and 1 80 d, respectively. Accordingly, the total element numbers for Grid A, Grid B and Grid C are 4358, 14,060 and 48,163, respectively. Figure 10(a) shows the image of Grid A, and the enlarged view of the grid cells near the cylinder is plotted in Figure 10(b).
The simulations were performed until final time t = 1000 on Grid A, Grid B and Grid C, respectively. We   show the numerical results for the case of Re = 150 on Grid B in Figure 11 for Mach number and Figure 12 for vorticity. The Kármán vortex street was clearly reproduced after the symmetrical flow structure in the earlier stage was broken.
To quantitatively analyze the numerical accuracy, we measure the drag coefficient C d and the lift coefficient C l which are calculated as where F D and F L are the drag force and the lift force exerted on the cylinder. We plot the drag and lift coefficients for the case calculated on Grid C in Figure 13.
In the initial stage of computation, the wake around the cylinder keeps symmetry, thus the lift coefficient C l approaches zero. With the development of instability, the vortex sheds with gradually increasing shedding amplitude in the downstream of the cylinder and finally the periodic regime established. The profile of lift coefficient C l looks very close to the representative results in the literature. Due to the periodic vortex shedding, the drag coefficient C d correspondingly fluctuates as shown in Figure 13. We then calculated the mean drag coefficient We measured the time-averaged drag coefficient C d and Strouhal number St for gradually refined meshes, Grid A, Grid B and Grid C as shown in Table 2. The Figure 15. Lock-in regime for flow past an oscillating cylinder in (f,A) plane at Re = 100 with the data collected from Koopmann and Anagnostopoulos's studies (Anagnostopoulos, 2000;Koopmann, 1967). The region above the data line is the lock-in region and otherwise unlocked region. time-averaged drag coefficients approaches 1.371 and the Strouhal number also converges with the mesh finer than Grid B. Therefore, we assume that Grid B can give adequate accuracy for later test cases. Compared with the existing results, the Strouhal number has a good agreement with results calculated by DG/FV method (Zhang et al., 2014) and Müller's simulation (2008).
We also calculated the cases with Re = 60,80,100 and 120 to evaluate the capability of the numerical model to simulate flows with different Reynolds number. Figure 14 compares the numerical results of Strouhal numbers and r.m.s. (root-mean-square) lift coefficients with the empirical function (Norberg, 2003) and Placzek's numerical results (Placzek, Sigrist, & Hamdouni, 2009). The good agreement of Strouhal number and r.m.s. lift coefficients demonstrates the accuracy of the present solver for stationary mesh, which provides a good base for implementation to moving mesh applications.

Oscillating cylinder
Viscous flow over oscillating cylinder provides a test bed for numerical solvers for fluid-structure interactions. One of the most interesting phenomena for flow past forced oscillating cylinder is the lock-in phenomenon, which has been reported in the literature to verify numerical models. The characteristics of vortices and wakes around the cylinder is significantly different from that of a fixed cylinder. We investigate the flow over a forced oscillating cylinder in this part with a Reynolds number Re = 100. The cylinder is forced to oscillate in the transversal direction of the main uniform flow with a sinusoidal motion in time defined by where y max is the maximal displacement which can be characterized by the non-dimensional amplitude A = y max /d. The oscillating frequency, f 0 , which associates with the Strouhal frequency f s for fixed cylinder, plays an essential role in lock-in phenomenon, which can be categorized by the frequency ratio f = f 0 /f s . Figure 15 represents the lock-in zone with respect to the amplitude A and the frequency ratio f where the data is collected from Koopmann and Anagnostopoulos's studies (Anagnostopoulos, 2000;Koopmann, 1967). According to Figure 15, we keep the amplitude constant as A = 0.3 to devise the locked and unlocked configurations with f = 0.5,0.9,1.1 and 1.5. We plot the drag coefficient and lift coefficient over time for four cases in Figure 16 with time characterized by the oscillating period T 0 . As expected in Figure 15, the case with f = 0.9 and f = 1.1 are locked where the vortex shedding frequency are synchronized with the cylinder oscillating frequency as shown in Figure 16(c ,d).
For locked configurations, the mean drag coefficient has a significant increase which reaches 1.56 for f = 0.9 and 1.83 for f = 1.1, while it is 1.38 for the fixed cylinder. The maximum lift coefficient reduces to 0.24 for f = 0.9 but increases to 0.86 for f = 1.1, in comparison with its value of 0.33 for the fixed cylinder. Unlocked wakes are formed with the frequency ratio f = 0.5 and f = 1.5 which represent more complicated phenomena than locked cases. The variation of lift coefficients does not follow the sinusoidal profile any more. The period of signal is over several cylinder oscillation cycles which is called as the beating period T b . The simpler case is f = 0.5 that the beating periodic T b is equal to the forced oscillation period T 0 . In this case, it is twice of the Strouhal period T s . More complex case is f = 1.5 where the beating period consists several oscillating periods as shown in Figure 16(d). The snapshots of vorticity contour lines within one oscillating period are plotted in Figure 17. The cylinder moves with a sinusoidal function in y direction: (a) t = t 0 , it is located at the origin of the domain; (b) t = t 0 + 1 4 T 0 , it moves to the maximum displacement in y direction; (c) t = t 0 + 1 2 T 0 , it returns to the origin; (d) t = t 0 + 3 4 T 0 , it moves to the maximum displacement in the reverse direction of y-axis.
Cases with a larger amplitude A = 2 are also calculated to show the ability of our proposed method for large mesh distortions. The oscillating period ratio is set as f = 0.5. We show the snapshots of vorticity contour lines within one cylinder oscillating period in Figure 18. Figure 19(a) gives the mesh image when the cylinder moves to the maximum displacement in y direction, and we magnify the area near the cylinder surface in Figure 19(b). By using RBF interpolation function, though large boundary movements are specified, the displacements of the boundary points are smoothly interpolated to the internal grid points based on the distance of the internal nodes to the boundary which is the essence of RBF function to ensure whole mesh maintains a good quality. In the simulation of forced oscillation problems, RBF interpolation is outside the Runge-Kutta loop, which makes the grid   velocity u g constant during the whole time step (including all sub-steps in the 3-stage Runge-Kutta method). It is in fact equivalent to the case that each grid point moves at a constant speed in the duration of each time step. In this way, the computational cost of RBF interpolation can be substantially reduced without losing numerical accuracy. Meanwhile, a coarsening technique is employed as presented in Bos (2010). We select every other five points in the cylinder boundary as control points for saving computational cost. We estimate the computational cost by running the case up to time t = 10 on a PC with an Intel Core i7-4790 CPU @ 3.60 GHz, the computational time for RBF interpolation is 48.63 s while the whole running time is 1433.35 s.

Conclusion remarks
In this paper, we present a multi-moment finite volume method for viscous flows on moving grids based on direct ALE formulation. Two kinds of moments, VIA and PV at cell vertices, are treated as the computational variables and used for constructing quadratic polynomials for moving triangular elements. Both VIA and PV are calculated simultaneously using the integral and differential forms of the governing equations, respectively, under the ALE framework. The VIA is updated by the finite volume method, where conventional Riemann solvers can be used for calculating the convective fluxes and variable gradients on cell edge required in approximating the viscous fluxes are simply calculated from the interpolation of neighboring cell centers. For updating PVs, Roe's Riemann solver is implemented for inviscid fluxes and viscous term is calculated by the T EC (Time-evolution Converting) formula. The GCL is obeyed for deriving the finite volume formulation for control volume variation.
Third-order Runge-Kutta scheme is implemented for not only computational variables but also mesh geometry quantities. Various mesh movement strategies are used for testing the performance of present solver on moving grids: (1) mesh moves with the local fluid velocity; (2) mesh moves following an analytical function; (3) mesh moves by RBF interpolation given boundary points displacements.
Numerical verifications have been carried out in Section 5 to show the performance of the proposed solver. The propagation of 2D isentropic vortex was computed to demonstrate that the present model is able to achieve third-order accuracy on arbitrarily moving triangular elements. The GCL is rigorously satisfied as shown in the free stream preservation test. Compressible flows with high Mach numbers, such as Saltzman problem and viscous double Mach reflection, are calculated to demonstrate the performance of present method. Flow past a cylinder is finally calculated on both stationary and moving grids to show the feasibility for fluid-structure interaction applications. It is pointed out that defining and using the physical variables at the vertices of mesh cells as new computational variables leads to an accurate and efficient numerical formulation for ALE implementation, which can be expected as a novel and promising solver for FSI applications. More applications to multi-material problems, especially strongly coupled FSI problems, are in progress.

Disclosure statement
No potential conflict of interest was reported by the authors.

Funding
This work was supported in part by JSPS KAKENHI [grant number 17K18838].