A reconstructed discontinuous Galerkin method for compressible flows on moving curved grids

A high-order accurate reconstructed discontinuous Galerkin (rDG) method is developed for compressible inviscid and viscous flows in arbitrary Lagrangian-Eulerian (ALE) formulation on moving and deforming unstructured curved meshes. Taylor basis functions in the rDG method are defined on the time-dependent domain, where the integration and computations are performed. The Geometric Conservation Law (GCL) is satisfied by modifying the grid velocity terms on the right-hand side of the discretized equations at Gauss quadrature points. A radial basis function (RBF) interpolation method is used for propagating the mesh motion of the boundary nodes to the interior of the mesh. A third order Explicit first stage, Single Diagonal coefficient, diagonally Implicit Runge-Kutta scheme (ESDIRK3) is employed for the temporal integration. A number of benchmark test cases are conducted to assess the accuracy and robustness of the rDG-ALE method for moving and deforming boundary problems. The numerical experiments indicate that the developed rDG method can attain the designed spatial and temporal orders of accuracy, and the RBF method is effective and robust to avoid excessive distortion and invalid elements near moving boundaries.


Introduction
Many engineering problems require the solution on variable geometries, such as aeroelasticity, fluid-structure interaction, flapping flight and rotor-stator flows in turbine passage. Due to the limitation of the measurement techniques in the experiment, numerical simulation of such problems is an important supplement to investigating and understanding these complex phenomena. The arbitrary Lagrangian-Eulerian (ALE) [1] formulation taking into account the mesh motion by nature, has been considered often and commonly used to solve such problems numerically.
For these engineering problems mentioned above, the ALE formulation usually combines consistently the advection due to mesh motion and the advection due to fluid motion, i.e., these two are solved simultaneously. This type of formulation is termed unsplit ALE [2][3][4], as opposed to the split ALE or Lagrange-plus-remap ALE which consists of three steps, a Lagrangian step, a rezoning step and a remapping step. The Lagrangian-plus-remap ALE method is mainly used in hydrodynamics where the evolution of flow is undergoing large deformation due to the strong compression or expansion. The application of the unsplit ALE method to the hydrodynamic problems has been reported in [2,[4][5][6]. In this paper, our focus is on the unsplit ALE method for the moving boundary/grid problems, typically the flow over moving airfoils.
In such problems, the velocity and/or position of the boundary nodes are usually known in practice, either from the prescription of an analytical expression or discrete data, or from the response of the solid in the fluid-structure interaction case. With the movement of the boundary nodes, the interior of mesh might be distorted or even invalid. Thus, either regenerating the mesh, or a smoothing procedure that propagates the boundary motion into the interior is necessary, to preserve the mesh quality. However, regeneration of the mesh is usually computationally very expensive, and the latter case, a smoothing procedure is preferred in general. Various methods have been proposed in the literature to smooth the interior mesh. One type of these methods requires the solution of a system of elliptic (Poisson-type) partial differential equations (PDEs), such as linear elasticity [7] and non-linear elasticity [8] methods. Another smoothing procedure is the interpolation method, for example, the radial basis function (RBF) interpolation [9][10][11] and Delaunay graph mapping [12]. In this work, we take advantage of the RBF interpolation based on its efficiency and mesh quality performance.
The high order methods, most commonly the discontinuous Galerkin methods (DGM), due to their higher computational efficiency compared with the low order method, have been investigated extensively in the Eulerian frame, and are gaining more and more interest in the ALE formulation.
The discontinuous Galerkin methods (DGM), originally introduced for solving the neutron transport equation by Reed and Hill [13], are widely used in computational fluid dynamics (CFD), owing to their many distinctive, and attractive features, such as flexibility to handle complex geometry, compact stencil for higher-order reconstruction, and amenability to parallelization and hp-adaptation. However, the DG methods have been recognized as expensive in terms of both computational costs and storage requirements. Indeed, compared to the finite element and finite volume methods, the DGM require solutions of systems of equations with more unknowns for the same grids. It is our belief that a lack of efficiency is one of reasons that prevents the application of DGM to engineering-type problems. In order to reduce high costs associated with the DGM, Dumbser et al. [14][15][16] introduced a new family of reconstructed DGM, termed PnPm schemes and referred to as rDG(PnPm) in this paper, where Pn indicates that a piecewise polynomial of degree of n is used to represent a DG solution, and Pm represents a reconstructed polynomial solution of degree of m (m n) that is used to compute the fluxes. The rDG(PnPm) schemes are designed to enhance the accuracy of the DGM by increasing the order of the underlying polynomial solution. The beauty of rDG(PnPm) schemes is that they provide a unified formulation for both FVM and DGM, and contain both classical FVM and standard DGM as two special cases of rDG(PnPm) schemes. When n 0, i.e. a piecewise constant polynomial is used to represent a numerical solution, rDG(P0Pm) is nothing but classical high order finite volume schemes, where a polynomial solution of degree m (m 1) is reconstructed from a piecewise constant solution. When m n, the reconstruction reduces to the identity operator, and rDG(PnPn) scheme yields a standard DG(Pn) method. For n 0, and n m, a new family of numerical methods from third-order of accuracy upwards is obtained. A number of algorithms are proposed [17][18][19][20][21][22] to construct the PnPm schemes, and their spatial convergence and effectiveness are validated [23][24][25][26][27].
Based on the higher performance of the rDG methods in the Eulerian frame, an extension to the arbitrary Lagrangian-Eulerian (ALE) formulation is natural and desirable. A number of DG-ALE methods in the literature have been proposed and investigated for the compressible flows. How to satisfy the Geometric Conservation Law (GCL) is a critical issue for the ALE formulation, especially for higher-order DG methods, where the basis functions being defined on the time-dependent or fixed reference domain/element will come into the picture, making the problem more complicated.
In general, there are two types of approaches to ensure the satisfaction of the GCL condition. The most consistent and elegant approach is the space-time DG formulation. This formulation is fully conservative in space and time, and the GCL is automatically satisfied. The space-time DG methods have been applied to both incompressible [28,29] and compressible flows [30][31][32][33]. However, the generation of the space-time meshes will require additional work. And, the space-time DG methods do not allow for explicit time stepping or implicit multi-step schemes [34].
The second type requires special treatment for the ALE formulations on a fixed or time-dependent mesh, for example, an additional equation for updating the transformation Jacobian, or a correction of the grid velocity terms. Lomtev et al. [35] presented a matrix-free DG-ALE method using spectral basis for 2D and 3D compressible viscous flows in moving domains. A force-directed algorithm from graph theory is used to update the grid while minimizing distortions. In the method proposed by Persson et al. [34], the governing equations in the time-dependent physical domain are transformed to the conservations laws in a fixed reference domain (the initial domain), the conservative variables defined on the initial domain are solved thereafter. A continuous mapping between the time-dependent physical domain and the fixed initial domain is introduced to take into account the mesh motion. The transformation Jacobian is updated in time to ensure the GCL condition. Nguyen [36] presented a DG-ALE method using an explicit fourth order TVD Runge-Kutta method and the freestream solution is shown to be preserved numerically. Mavriplis [37] derived a general approach inspired from the space-time formulation to update the grid velocity terms at the Gauss quadrature points such that the GCL is satisfied. Ren and Xu et al. [38] introduced a DG-ALE method based on the gaskinetic theory. In their method, the initial domain is taken as the reference domain, and the basis functions are mapped from this reference domain to the time-dependent physical domain. The computations are conducted in the physical domain. A space-time type integration is used to obtain the discretized equations and then the gas kinetic flux is computed to advance the solution. This method is shown to preserve the uniform flow automatically, and applied to several moving boundary problems.
In the current work, the reconstructed discontinuous Galerkin (rDG) method is extended to the ALE formulation, to simulate flow problems with moving or deforming grids. The Taylor basis functions are defined on the time-dependent physical domain, as a continuation of the rDG method in the Eulerian formulation. For better resolution at the wall of the airfoil, the curved elements are used in the whole domain. The third order temporal scheme ESDIRK3 (Explicit first stage, Single Diagonal coefficient, diagonally Implicit Runge-Kutta) is employed for time marching. To enforce the GCL, we follow a similar idea as in [37], by updating the grid velocity terms at Gauss quadrature points. The radial basis function (RBF) interpolation method has been chosen as the mesh smoothing algorithm to provide the mesh motion for the interior nodes, given the displacement or velocity at the boundary nodes. Numerical test cases are set up to demonstrate the accuracy and capability of the derived rDG-ALE method.
The remainder of this paper is organized as follows. The governing equations are presented in Section 2. The rDG formulation in ALE frame is derived in Section 3. Section 4 discusses how to satisfy the Geometric Conservation Law (GCL) condition. The RBF method for the mesh movement will be described in Section 5. Section 6 presents a set of numerical examples. Finally conclusions are given in Section 7.

Governing equations
The Navier-Stokes equations governing the unsteady compressible viscous flows can be expressed as where the summation convention has been used. The conservative variable U, advective flux vector F and viscous flux vector G are defined by Here , p and e denote the density, pressure and specific total energy of the fluid, respectively, and u i is the velocity component of the flow in the coordinate direction x i . The pressure can be computed from the equation of state which is valid for perfect gas. The ratio of the specific heats is assumed to be constant and equal to 1.4. The viscous stress tensor ij and heat flux vector q j are given by In the above equations, T is the temperature of the fluid, Pr the laminar Prandtl number, which is taken as 0.7 for air. represents the molecular viscosity, which can be determined through Sutherlands law where 0 is the viscosity at the reference temperature T 0 and S 110K. The temperature T of the fluid is determined by where R is the gas constant. The Euler equations can be obtained if the effects of viscosity and thermal conduction are neglected in Eq. 1.

Reconstructed discontinuous Galerkin discretization
Since the grid velocity only contributes to and appears in the convective term, it's sufficient to consider the compressible Euler equations here for deriving the ALE formulation. The differential governing equation can be written in the form where t is the usual Eulerian time derivative. Defining the (Taylor) basis functions j associated with the moving volume t e , multiplying j on the equation above and integrating on the moving volume, one will get immediately Using Reynolds Transport Theorem where V g is the grid velocity, and the divergence theorem Eq. 8 becomes The fundamental ALE relation for the total time derivative, Eulerian time derivative and the spatial gradient is We note that the total time derivative d dt is with respect to the grid velocity, in contrast to the material derivative with respect to the fluid velocity. And also, the d dt here is equivalent to the referential time derivative t X in [38] where the subscript X was used to indicate fixing the mapping position in the initial mesh; in our case, although we are not performing an explicit mapping from the initial configuration to the time-dependent element as in [38], we are essentially using an implicit one-to-one mapping from the points (e.g., the Gauss quadratures) on a standard reference element to the grid coordinates in the timedependent element. Thus in either case, we are tracking the grids when using this total time derivative, i.e., holding the grid point label fixed.
Plugging this relation into Eq. 11, we end up with the rDG-ALE formulation Note that Taylor basis functions [39] are used in the present work, which are defined on the time-dependent physical element, and therefore are time-dependent, i.e., j j t except for j 1, 1 1 which is a constant and thus d 1 dt is zero everywhere and at any instant time. The resulting ALE formulation will be discretized in space using the reconstructed discontinuous Galerkin method as in the Eulerian formulation [18,19,22,23,27].

Geometric conservation law (GCL)
The Geometric Conservation Law (GCL) states that the fully discretized equations should preserve a constant solution under arbitrary mesh movement, i.e., given a uniform initial flow condition, the solution should be maintained by the devised ALE solver.
Consider the BDF1 temporal discretization scheme.
Plugging a constant solution into the equation above leads to This is the GCL equation that needs to be satisfied by the rDG method with BDF1 scheme. Note that for a constant solution U, the Eulerian flux F vanishes due to its consistency property. Unfortunately, this equation will not hold in general, even at the DG(P0) (Finite Volume) level which corresponds to 1. To verify, plugging in the basis 1, the GCL equation reduces to 1 t n 1 e n e n 1 e V g nd (16) where the mesh velocity is defined as The GCL equation above for FV states that the volume change between two successively discretized time levels n and n 1, equals to the volume flux at the element interfaces at time level n 1 (time level n if using explicit time stepping). Unfortunately, this is not true in general.
Here we follow the idea from the work of Mavriplis et al. [37]. The major difference between the current work and that in [37] is that, the basis functions in [37] are defined on the reference element, while the Taylor-basis functions used in the current work are defined on the time-dependent physical element. This difference leads to two consequences. First, in [37], the substantial derivative of the basis functions with respect to the mesh motion will vanish, due to the definition of the basis function and the fact that the reference element is invariant in time-thus the basis functions simply move with the grid velocity and their values will never change. In contrast, for the Taylor-basis in this work, they are defined on the time-dependent element, and if we consider a quadrature point, although it also travels with the mesh velocity, the value of its basis function is changing in time inherently. Thus, an extra term containing the total derivative of the basis functions will appear in this formulation. Second, in [37], the solution expansion is performed on the reference element using the basis functions therein, thus the domain integral takes on a different form than this work, i.e., in [37] the inverse of the transformation Jacobian emerges when transforming the spatial gradient from reference element to physical element, while in the current work the spatial gradient on the physical element can be used directly without any transformation.
The idea of how to enforce the GCL condition stems from the space-time DG method. Performing a full integration on Eq. 13 in both space and time (-on a space-time element, as shown in Fig. 1) leads to This GCL equation is automatically satisfied, provided enough Gauss quadrature points are employed in both space and time to conduct the numerical integration. Inspired from this property of space-time integration, a straightforward way to satisfy the BDF1 GCL equation is to modify the right-hand side of Eq. 15 such that it works the same way as its space-time DG counterpart Eq. 19. More explicitly, we want to have the following conditions hold However, as addressed in [37], although the values for these spatial integrals of grid velocities can be obtained, they do not yield grid velocities themselves, and hence can not be directly used in evaluating the integrals in Eq. 14.
Recall that eventually the numerical integration will be performed on the reference element where the Gauss quadrature points reside, to this end, it's beneficial to write the above equations in terms of numerical integration on the reference element.
For BDF1 equation Also the two GCL equations corresponding to the above discretizations can be rewritten accordingly as We require the following conditions hold at each Gauss quadrature point For the computation of the total time derivative of the basis functions, the finite difference could be used. After plugging these grid velocity terms into Eq. 14, the right-hand side could be evaluated at time level n 1.
For the higher-order ESDIRK temporal scheme (Explicit first stage, Single Diagonal coefficient, diagonally Implicit Runge-Kutta) [ where each stage corresponds to time t s t n c s t. Similarly, by comparing each stage of the ESDIRK scheme with the space-time DG formulation, the following conditions could be obtained, such that each stage will satisfy the GCL.
Note that these conditions are required for each Gauss quadrature point. For the ESDIRK scheme, this system of equations can be easily solved by forward substitution. We choose the third order ESDIRK3 scheme for this work.

Mesh motion
In both the curved mesh generation and the pure mesh movement, the motion of the boundary nodes is required to propagate to the interior ones properly, in order to avoid invalid elements near the boundary and to maintain good mesh quality. That is, given the displacements of the boundary nodes, it is desired to obtain the displacements for all the interior nodes. With the motion of all the nodes in the mesh determined, one can then move the mesh coordinates accordingly to accommodate the physics solver, e.g., the ALE formulation for the fluid flow.
There are several types of methods in the literature dealing with the movement of the interior nodes. Model based on physics (solid mechanics) is a popular one. Examples are linear elasticity [7], non-linear elasticity [8], linear spring analogy [42] and torsion spring analogy [43]. This type of method requires the solution of a system of equations, elliptic (Poisson-type) partial differential equations (PDE) in the case of the linear elasticity model, thus is relatively computationally expensive. The second type is the interpolation method or algebraic method, including radial basis function (RBF) interpolation [9], explicit interpolation [44] and Delaunay graph mapping [12]. Kashi and Luo [10,11] compared the performance between RBF interpolation and some of the other methods. It is found that the RBF method, based on the interpolation technique which is relatively fast, can give good and robust results in general, especially for large or rotational deformations.
In the following moving airfoil test cases, the displacement or grid velocity of the wall boundary are provided in advance. The RBF interpolation method is then used to compute the deformation of the interior nodes, such that the nodes on the airfoil are in a rigid motion, and those far away from the wall are kept static. We note that if a curved mesh is to be used, then the high order nodes (e.g., midpoint of an edge in 2D) will also come into play.  Here we briefly recall the methodology of the RBF method. Interested readers are encouraged to refer to [10,11] and references therein for more information. The RBF method states that, given N b boundary points in the mesh, we require the displacement at any point x (both boundary nodes and the interior nodes) in the mesh satisfy where r is the basis function, s represents the displacement vector field with 2 components x-and y-in 2D ( x-, y-and z-in 3D), and a j are the coefficients to be determined. Note that the RBF method will be applied to a scalar field, i.e., to each component individually. In this work, we choose the modified Wendland's C2 function [45] as the basis where r s is called the support radius, indicating how far the boundary deformation will propagate to the interior. Given the displacement of the boundary nodes, by substituting these known coordinates and displacements into Eq. 29, we can construct a linear system of equations for the coefficients a j in each dimension in space (31) or rewritten as Aa s (32) When using Wendland's C2 function, the resulting symmetric system is claimed to be positive definite [45], thus can be easily solved by conjugate gradient (CG) or Choleskey decomposition methods, or even direct methods. After solving the linear system, we can evaluate the displacement for any interior nodes using Eq. 29.

Numerical examples
In all the following test cases, the 3rd-order accurate temporal discretization ESDIRK3 scheme is used, and the curved elements (Q2) are chosen for all the grids, in order to test the applicability of the rDG-ALE method.

Uniform flow
The first test case is the uniform flow on deforming grids. The initial mesh is uniform on a square domain defined in 0 x, y 1.0, and the mesh motion is prescribed by where L x L y 1.0 and T 1.0 are the reference length and time, A x A y 0.025 is the deformation amplitude, n x n y 4 and n t 0.5 are for the deformation frequency in space and time. The mesh deformation process is illustrated in Fig. 2. The L 2 norm of the numerical error is measured and observed to be around 10 12 , thus verifying the Geometric Conservation Law (GCL).

Convection of an isentropic vortex
The convection of a 2D inviscid isentropic vortex is considered in this test case. All the three spatial discretization methods, DG(P1), rDG(P1P2) and DG(P2) are employed in the computations. The square domain 0 x, y 1.0 is considered for this test case as well, with the same mesh deformation configuration as in the uniform flow case.
where r x x 0 2 y y 0 2 and 1.0, 4.0, 1.4. The other conservative variables can be determined by the follows: The initial mesh and density distribution are shown in Fig. 3, and the computed solution is shown in Fig. 4, for DG(P1) and rDG(P1P2), respectively. A mesh refinement study is performed, in order to evaluate the convergence rate for different spatial discretization schemes. The computed L 2 norm for density is shown in Fig. 5. We can see the desired orders of accuracy are achieved, and by comparing rDG(P1P2) with DG(P1), we can see rDG(P1P2) not only has higher convergence rate, but the absolute error for rDG(P1P2) is also smaller than the counterpart of DG(P1). Figure 6 presents the spatial error versus the number of degrees of freedoms (dofs) for DG(P1) and rDG(P1P2). One can see that for a given number of dofs, the numerical error of rDG(P1P2) is smaller than that of DG(P1), which indicates that rDG(P1P2) is computationally more efficient.
To demonstrate the temporal order of accuracy of the ESDIRK3 scheme, we perform a time-step refinement study, and compute the numerical errors. The results are shown in Fig. 7. We can see that the designed 3rd-order of convergence has been achieved.   , respectively. Given the above harmonic motion for the points at wall boundary and the static points at the far-filed boundary, the movement of the interior nodes is interpolated from the RBF method. With the airfoil motion and the flow field being periodic, the unsteady solutions after a reasonably long time of computation should be used for analysis.
As in the experiments by Landon [46], the normal force coefficients C n and pitching moment coefficients C m are defined as is the pressure coefficient, and the subscript denotes the freestream quantity, while superscript L and U denote the lower and upper surfaces respectively.
The computed normal force coefficients and pitching moment coefficients, with respect to the angles of attack and time, respectively, are shown in Figs. 10 and 11. The numerical solutions are in good agreement with the experimental data by Landon [46]. Note that while the experiment is viscous, we treat this problem as an inviscid one, which may lead to some difference between the numerical solution and the experimental data.

Pitching NACA0015 airfoil
The laminar flow around a rapidly pitching NACA0015 airfoil is considered in this test case. The freestream has a Mach number Ma 0.2 with specific heat ratio 1.4. The Reynolds number based on the chord length and the freestream condition is Re 10, 000. The pitching axis is at a quarter chord length from the leading edge, and the mesh motion is prescribed by the pitching rate of the airfoil: t 0 1.0 exp 4.6t t 0 rad/s, where t 0 c U is the reference time, 0 0.6U c, with c being the chord length and U the freestream velocity. The flow field computed at zero angle of attack is used as the initial condition for this unsteady simulation.
The computational domain is 15, 15 15, 15 , with the airfoil located at the center. The initial mesh and the zoomed-in mesh near the airfoil at three different angles of attack are shown in Fig. 12. As usual, the RBF method is used to compute the movement of the interior mesh nodes, i.e., the nodes on the airfoil move according to the above variation of the angle of attack, while those at far-field are kept static. It can be seen that the mesh quality near the wall is well maintained, even at a large angle of attack. Two sets of curved grids are used in this numerical experiment. For the coarse mesh, the minimum and maximum wall normal spacings are 0.005 and 0.008, respectively, with 180 nodes on the airfoil, while for the fine mesh, the minimum and maximum normal spacings are 0.00037 and 0.0006, with 200 nodes on the airfoil.
In Fig. 13, the instantaneous normalized vorticities z c U at three different angles of attack are presented. It can be seen the vorticities near the wall of the airfoil formed and transported towards downstream. In Fig. 14, the computed lift coefficients and drag coefficients are compared with the simulation data from Visbal et al. [47] and Ren and Xu et al. [38]. From the comparison we can see that the computed results with coarse mesh    are in good agreement with the reference data, while those on the fine mesh have relatively large discrepancies with the references, especially for the lift coefficients. In general, fine mesh should better capture the flow structures due to its finer resolution. While the reason for the occurrence of this phenomenon is not quite clear and is still under investigation, there could be a number of factors that might contribute, for example, as pointed out in [38], how to compute the mesh velocity at the cell interface has a significant effect on the numerical solution, and a variable mesh velocity along cell interface is necessary for a high order method, after comparing the results from a piecewise constant grid velocity and a variable mesh velocity [38]. In the current work, a quadratic interpolation from 16 Lift coefficients (a) and drag coefficients (b) versus normalized time for the plunging SD7003 case the three nodes of an interface (two end nodes and one midpoint) is used to obtain the variable mesh velocities.

Plunging SD7003 airfoil
The flow over a moving SD7003 airfoil has been investigated both experimentally and numerically in the literature. McGowan et al. [48,49] conducted a set of experiments on Fig. 17 Comparison of normalized vorticities with experiment at t 0T for the plunging SD7003 case the purely plunging or pitching SD7003 airfoil, using particle image velocimetry (PIV) in a water tunnel, besides, the numerical solutions from both CFL3D and an immersed boundary method are compared with the experimental data. Visbal et al. [50] performed the computations in which the grid is moved in a rigid fashion, as opposed to the smoothing or interpolation approaches, e.g., the RBF method. In this work, the flow over the high-frequency plunging SD7003 airfoil is simulated using the rDG-ALE method. The Reynolds number based on the chord length and the freestream velocity is Re 10, 000, and the specific heat ratio 5 3. The airfoil is set at a static angle of attack 0 4 degree. The plunging motion is prescribed as h t h 0 sin 2kU t c where h 0 0.05c is the plunging amplitude and k 3.93 is the reduced Fig. 19 Comparison of normalized vorticities with experiment at t 2 4T for the plunging SD7003 case frequency, with c the chord length and U the freestream velocity. Although the motioninduced angle of attack could be as high as 21.5 degree which leads to unsteady separation and the generation of dynamic-stall-like vortices at the leading edge, the transition effects are observed to be minor for this low Reynolds number (Re 10, 000) [50]. To show the compressibility effect for this test case, we choose two Mach numbers, Ma 0.2 and 0.05, as in Ren and Xu's paper [38]. The computational domain is a circle region with radius R 20. Figure 15 shows the initial mesh and zoomed-in mesh. The number of nodes on the airfoil is 200 and the wall normal spacing is approximately 0.0006. The time-dependent lift coefficients and drag coefficients are compared with the experimental data from McGowan et al. [48,49]., and shown in Fig. 16. We can see that at Ma 0.05 the lift coefficient has an excellent agreement with the experimental data, while a phase lag is observed for both the lift and drag coefficients at Ma 0.2. The computed vorticities are also compared with the experimental data, shown in Figs. 17, 18, 19 and 20, respectively, for two time instances. We can see they agree well with each other.

Conclusions
A reconstructed discontinuous Galerkin method has been presented for solving the unsteady compressible Navier-Stokes equations in an arbitrary Lagrangian-Eulerian (ALE) formulation on moving and deforming curved grids. The Taylor basis functions used in the rDG method are defined on the time-dependent physical space, where the GCL condition is ensured by modifying the grid velocity terms on the right-hand side of the discretized equations. The radial basis function (RBF) interpolation method is developed to provide the mesh motion for the interior nodes, given the motion of the boundary nodes. A third order ESDIRK3 integration method is used to advance solutions in time. The developed rDG method has been assessed for a number of benchmark test cases. The numerical results obtained by the rDG method are compared with the available experimental data and reference solutions in the literature, demonstrating that the developed rDG-ALE method is able to achieve the designed spatial and temporal orders of accuracy and provides an effective approach for solving moving or deforming domain problems.