On the Geometric Conservation Law for the Non Linear Frequency Domain and Time-Spectral Methods
aa r X i v : . [ c s . C E ] A p r On the Geometric Conservation Law for the Non Linear FrequencyDomain and Time-Spectral Methods
Marc Benoit, Siva Nadarajah
Department of Mechanical Engineering, McGill University, Montreal, Quebec, H3A 0C3, Canada
Abstract
The aim of this paper is to present and validate two new procedures to enforce the GeometricConservation Law (GCL) on a moving grid for an Arbitrary Lagrangian Eulerian (ALE) formulationof the Euler equations discretized in time for either the Non Linear Frequency Domain (NLFD) orTime-Spectral (TS) methods. The equations are spatially discretized by a structured finite-volumescheme on a hexahedral mesh. The derived methodologies follow a general approach where thepositions and the velocities of the grid points are known at each time step. The integrated facemesh velocities are derived either from the Approximation of the Exact Volumetric Increments(AEVI) relative to the undeformed mesh or exactly computed based on a Trilinear Mapping (TRI-MAP) between the physical space and the computational domain. The accuracy of the AEVImethod highly depends on the computation of the volumetric increments and limits the temporal-order of accuracy of the deduced integrated face mesh velocities to between one and two. Thusdefeating the purpose of the NLFD method which possesses spectral rate of convergence. However,the TRI-MAP method has proven to be more computationally efficient, ensuring the satisfactionof the GCL once the convergence of the time derivative of the cell volume is reached in Fourierspace. The methods are validated numerically by verifying the conservation of uniform flow and bycomparing the integrated face mesh velocities to the exact values derived from the mapping.
Keywords :
Non Linear Frequency Domain (NLFD) method; Time-Spectral method; Geomet-ric Conservation Law (GCL); Deforming mesh; Computational fluid dynamics; Numerical method
Preprint submitted to Journal April 12, 2018 . Introduction
The accurate computation of unsteady aerodynamic flows in a reasonable amount of time stillpresents a challenge in the field of computational fluid dynamics. Compared to steady flow prob-lems which only require an accurate spatial discretization, unsteady flow solvers have to provide anaccurate time resolution of the flow. Until lately, the most popular approaches were time marching techniques for which the solution is constructed in time from an initial free-stream solution. Despiteacceleration techniques such as multigrid or local time stepping, these methods remain computa-tionally costly, partly because of the transient effects, requiring numerous time steps to converge.However in the case of periodic flows encountered in problems such as aeroleastic simulations orturbomachinery, the Fourier collocation method can be used to accurately and efficiently represent the solution. A nonlinear harmonic method known as the Harmonic Balance (HB) was initiallyintroduced by Hall et al. [1]. Then, McMullen et al. [2–4] developed the Non Linear Frequency Do-main (NLFD) method in order to solve the Euler and Navier-Stokes equations directly in frequencydomain. An alternative to this approach is the Time-Spectral (TS) method presented by Gopinathet al. [5, 6], it avoids the explicit use of a Discrete Fourier Transform (DFT) and discretizes the temporal derivative operator through a Fourier collocation matrix. These various methods prove tosignificantly decrease the required time to obtain the solution compared to time marching solvers.The NLFD method was used to develop a two dimensional aeroelastic flow solver [7] and thenexpanded to three-dimensional cases on moving grids [8]. In order to study aeroelastic problems, themesh is moved in time and additional care has to be taken to compute the mesh velocities and the metrics from the physical to the computational space. Thomas et al. [9] were the first to formallydefine the necessity to solve additional laws to preserve the conservation of the solver numericalscheme. Termed as the Geometric Conservation Law, it is composed of two subsets of laws knownas the Surface Conservation Law (SCL) and the Volume Conservation Law (VCL). A mathematicalinterpretation of the SCL relates that any cell volume has to be closed by its surfaces whereas the VCL states that the temporal rate of change of the cell volume is equal to the sum of the temporalrate of change of the algebraic volumes swept by each face enclosing it through time. SCL differsfrom the VCL in the way that they need to be verified even for fixed grids (steady state) whilethe VCL appears only on deforming grids. The violation of any of these laws may result in errorin the flow solution, for instance it was reported that the violation of the GCL leads to inaccurate flutter prediction for aeroelastic cases [10]. Further investigation on time marching schemes clarified2he theoretical status of the GCL, exposing its link to temporal-order accuracy [11], or stabilityconditions [12]. In addition, the Discrete Geometric Conservation Law (DGCL) were derived inorder to preserve the time accuracy for high temporal-order schemes [13]. These methods are welladapted for time marching approaches, but their extension to the NLFD or Time-Spectral methods is not straightforward since it becomes necessary to compute all quantities : state vector, fluxes,mesh positions and mesh velocities, at all time steps before applying the Fourier discretization.In this work, we focus on the VCL part of the GCL which arises from the unsteadiness of the flow.The framework of our study is based on the previous development of a three dimensional aeroelasticsolver by Tardiff et al. [8] who introduced a methodology to enforce the GCL. The expression of the governing equations according to an Arbitrary Lagrangian-Eulerian formulation, the principlesof the NLFD discretization and the dynamic mesh deformation using Radial Basis Functions arepresented in section 2. Additional analytical developments are introduced to demonstrate thelimitations of the approach of Tardiff et al. [8]. A modified method is derived to satisfy the GCLby computing the integrated face mesh velocities at a second-order accuracy in time. Since such a condition limits the inherent spectral in time accuracy of the NLFD approach, an alternativemethod based on a trilinear mapping between the physical and computational space is introduced.The results are then extent to the Time-Spectral method. These various approaches are presentedsection 3. The different methodologies are numerically investigated in order to verify that the GCLare enforced through a correct computation of the integrated face mesh velocities in section 4 and the results discussed in section 5.
2. Discretization of the governing equation and mesh deformation
This section presents the formulation of the Euler equations on a moving mesh using the Ar-bitrary Lagrangian-Eulerian approach, its discretization using the Non-Linear Frequency Domain,and the mesh deformation method through the Radial Basis Functions. When solving the Euler equations on a moving grid a popular approach is to use an ArbitraryLagrangian-Eulerian (ALE) formulation [14]. For a control volume Ω enclosed by a boundary ∂ Ωand without source terms, the integral form of this formulation is given as follows : ˆ Ω w dΩ + ˛ ∂ Ω F Mc d S = 0 , (1)3ith the state vector w and the vector of the convective fluxes on a moving grid F Mc expressed as : w = ρρuρvρwρE ; F Mc = F c − V t w = ρVρuV + n x pρvV + n y pρwV + n z pρHV − ρV t ρuV t ρvV t ρwV t ρEV t , (2)where F c would be the convective flux vector on a fixed grid and ρ , u , v , w , E , p and H arerespectively, the density, the Cartesian velocities of the fluid, the total energy per unit mass, thepressure and the total enthalpy per unit mass defined by : H = E + pρ . (3)In order to close the system of equations, the pressure is evaluated under the assumption ofideal gas through the state equation (4) : p = ( γ − (cid:18) ρE − ( ρu ) + ( ρv ) + ( ρw ) ρ (cid:19) . (4)Also, V is the contravariant velocity of the fluid and V t is the contravariant velocity of theboundary enclosing the control volume. n x , n y and n z are the components of the boundary unitnormal vector pointing outward of the control volume. It yields : V = un x + vn y + wn z ,V t = ∂x∂t n x + ∂y∂t n y + ∂z∂t n z . (5)By introducing a discretized control volume and an artificial dissipation flux vector F d to avoidan odd-even decoupling of the solution and to increase the accuracy at discontinuities, the equation(1) can be written under a semi-discretized form as : ∂ (Ω w ) ∂t + X ∂ Ω (cid:2) ( F Mc ) S − F d (cid:3) = 0 . (6)The previous set of equations has to hold for each control volume and can be expressed as asemi-discrete system of ordinary differential equations in time : ∂ (Ω w ) ∂t + R ( w ) = 0 , (7)4here R ( w ) = X ∂ Ω (cid:2) ( F Mc ) S − F d (cid:3) is the discrete residual vector. In this work, the discretizationin space is performed according to the finite volume method on a structured hexahedral grid, themodified convective flux is computed as the average of the fluxes at a cell face and the artificial dis-sipation is evaluated using the JST scheme [15]. The residual vector is calculated as the summation over the faces of the control volume of the different fluxes. The temporal discretization of the flow solver employs the NLFD approach developed by Mc-Mullen et al. [2]. Under the assumption that both the modified state vector ¯ w = Ω w and theresidual vector R ( w ) are periodic in physical time, the two quantities can be expanded as discreteFourier series using a finite number of harmonics equations (8) and (9) :¯ w ( t ) = N X k = − N ˆw k e i (2 πk/T ) t , (8) R ( w ( t )) = N X k = − N ˆR k e i (2 πk/T ) t , (9)where i = √− T is the time period, k is the wave number, and N is thenumber of modes employed in the Discrete Fourier Transform (DFT). The k th Fourier coefficients ˆw k and ˆR k are given by the following equations (10) and (11), for − N ≤ k ≤ N : ˆw k = 12 N + 1 N X n =0 Ω( t n ) w ( t n )e − i (2 πk/T ) t n , (10) ˆR k = 12 N + 1 N X n =0 R ( w ( t n ))e − i (2 πk/T ) t n , (11)where their computations require the sampling of the modified state vector and the residual vectorfor N ts = 2 N + 1 time steps at equally spaced time instances such that the n th time sample t n is : t n = n N + 1 T , for n = 0 , .., N. (12)At this point, it is important to emphasize that the state and residual vectors need to beevaluated at all time instances before transferring in the Fourier domain, this is a fundamental5ifference with the time marching approach. The Fourier representation is then substituted intothe semi-discrete form of the Euler equations (7) to yield : ∂∂t N X k = − N ˆw k e i (2 πk/T ) t ! + N X k = − N ˆR k e i (2 πk/T ) t = 0 , (13) ⇔ N X k = − N i πkT ˆw k e i (2 πk/T ) t + N X k = − N ˆR k e i (2 πk/T ) t = 0 . (14)By exploiting the orthogonality property of the Fourier basis, this leads to a set of 2 N + 1equations (15), each being associated to a wave number k : i πkT ˆ w k + ˆ R k = 0 for − N ≤ k ≤ N. (15)Since the representation of ˆ R k as a function of ˆ w k is not straightforward, an unsteady residualˆ R ∗ k is defined and driven to zero using a pseudo-time marching approach such that : ˆ R ∗ k = i πkT ˆ w k + ˆ R k ∂ ˆ w k ∂t ∗ + ˆ R ∗ k = 0 , for − N ≤ k ≤ N. (16)Thus at convergence, ˆ R ∗ k = 0 and the equations (16) are satisfied for each wave number.The new periodic solution is then transferred back to the physical time domain using an InverseFourier Discrete Transform (IDFT) and evaluated at each time instance t n by dividing by thevolume : w ( t n ) = ¯ w ( t n )Ω( t n ) , for 0 ≤ n ≤ N. (17)The equation in pseudo-time can be solved using any time-stepping scheme. In this work, weuse an hybrid five-stage Runge-Kutta scheme with blending coefficients for the artificial dissipation [16]. The deformation of the mesh is performed using the Radial Basis Functions (RBF) [8]. Themethod is based on the assumption that the movement of all grid points can be interpolated fromthe a priori known motion of a set of points called the RBF points. In this study, the RBF pointsare always a subset of the grid points at the boundary of the domain, their displacements relative6o the undeformed mesh are prescribed at each time instance using analytical functions. Because ofthe NLFD method, the mesh positions and velocities are therefore computed and stored for all N ts time steps. For any grid point p of position vector x p in the undeformed mesh, its displacement inthe x -direction s x ( x p , t ) is defined as : s x ( x p , t ) = N rbf X i =1 α i ( t ) φ ( || x p − x i || ) , (18)where N rbf is the number of RBF points, α i are the interpolating coefficients, x i is the positionvector of the i th RBF point in the undeformed grid and φ is some basis function depending onthe Euclidean distance || x p − x i || between the points p and i . In this work, Wendland C0’s basisfunction [17] is considered, it is defined as follows : (1 − ξ ) if ξ <
10 if ξ ≥ , with ξ = || x p − x i || R , (19)where R is the support radius relative to the surface of RBF points. Since the equation (18) holdsfor any grid point whether it is an RBF point or a standard grid point, in the following, the RBFpoints are denoted with the subscript r while the grid (or volume) points are denoted with thesubscript v . Then in the x -direction, the displacements of all RBF points and the interpolateddisplacements of all grid points are regrouped respectively in the vector ∆x r and in the vector ∆x v . Therefore the a priori unknown displacements ∆x v are obtained through equation (20) : ∆x v = A ( M − ) ∆x r , (20)where : M = φ r r φ r r . . . φ r r Nrbf φ r r . . . ...... φ r Nrbf r · · · φ r Nrbf r Nrbf , A = φ v r φ v r . . . φ v r Nrbf φ v r . . . ...... φ v Ngrid r · · · φ v Ngrid r Nrbf (21)with : φ v i r j = φ (cid:0) || x v i − x r j || (cid:1) (22)and N grid is the total number of grid points. The displacements in the y and z directions can be com-puted with the same matrices given in equation (21), by considering the RBF points displacementsin the corresponding direction. a priori known velocities of the RBF points which leadsto the following expression : v dir,v = A ( M − ) v dir,r , (23)where v dir,v is the vector of the velocities of the grid points and v dir,r is the vector of the velocitiesof the RBF points and the direction is given by dir = x, y, or z .
3. Derivation and enforcement of the Geometric Conservation Law
As previously stated our interest is focused on the Volume Conservation Law aspect of the GCL.Under integral form the VCL for a control volume Ω enclosed by a boundary ∂ Ω can be written asfollows : ∂∂t ˆ Ω dΩ − ˛ ∂ Ω ( V t · n )d S = 0 . (24)where V t = (cid:18) ∂x∂t , ∂y∂t , ∂z∂t (cid:19) is the mesh velocity vector and n is the normal vector to the surface ∂ Ω. The law relates only on geometrical considerations and is always satisfied under continuousform and implicitly satisfied for rigid grid motion. It arises from the deformation of the meshand is closely related to the preservation of uniform flow by the numerical scheme. Therefore inorder to obtain a consistent solution method, the GCL must be discretized using the same numericalscheme employed to discretize the primary conservation laws [11]. In our case, it yields a hexahedral stuctured finite-volume framework and a temporal discretization using the NLFD method. A firstapproach to enforce the VCL in the NLFD context was presented by Tardiff et al. [8] but moreinvestigation is needed to determine its limitations. In this section, further developments are addedto this approach which expose its analytical limits and a new method is proposed.Considering any discretized control volume Ω enclosed by N f faces, then equation (24) can bewritten as : ∂ Ω ∂t − N f X m =1 ¨ ∂ Ω m ( V t · n m ) dS = 0 , (25)where n m is the unit normal vector to the face ∂ Ω m . Then the integrated face mesh velocities(IFMV) G m ( t ) corresponding to the temporal rate of change of the algebraic volume swept by each8ace through time are introduced in equation (26) : G m ( t ) = ¨ ∂ Ω m ( V t · n m ) dS, (26)and also G ( t ) is the sum of the IFMV over all faces of the control volume : G ( t ) = N f X m =1 G m ( t ) . (27)Then equation (25) can be written as : ∂ Ω ∂t − G ( t ) = 0 . (28)Under the assumption that the volume Ω and the sum of the integrated face mesh velocities G areperiodic functions of time, the NLFD discretization can be applied :Ω( t ) = N X k = − N ˆΩ k e i (2 πk/T ) t , (29) G ( t ) = N X k = − N ˆ G k e i (2 πk/T ) t . (30)By substituting these expressions into equation (28), yields : ∂∂t N X k = − N ˆΩ k e i (2 πk/T ) t ! − N X k = − N ˆ G k e i (2 πk/T ) t = 0 (31) ⇔ N X k = − N i πkT ˆΩ k e i (2 πk/T ) t ! − N X k = − N ˆ G k e i (2 πk/T ) t = 0 . (32)Then by exploiting the orthogonality property of the Fourier basis, it leads to a system of 2 N + 1equations, each corresponding to a wave number k : i πkT ˆΩ k = ˆ G k for − N ≤ k ≤ N. (33)The set of equations (33) provides the necessary condition to enforce the GCL in the NLFD approach.Such criterion is not satisfied in general and has to be enforced through the correct computation ofthe cell volume and the integrated face mesh velocities, in a way consistent with the solver numericalscheme. Since the volume is usually exactly known, one popular approach in time marching methodsis to split the GCL over each face [13, 18, 19]. In the current framework, the volume of a cell can9e expressed as the sum of the volume at a reference initial instant t and the algebraic (positiveor negative) volumetric increments due to each face Ω m relative to this reference instant :Ω( t ) = Ω( t ) + N f X m =1 Ω m ( t ) . (34)By substituting relations (34) and (27) into the equation (28), yields : N f X m =1 (cid:18) ∂ Ω m ∂t − G m ( t ) (cid:19) = 0 . (35)Then for each face m enclosing the discretized control volume, we need to ensure the relation (36) : ∂ Ω m ∂t = G m ( t ) . (36)However, even if the positions of the mesh vertices and their velocities are known at all time instances from the dynamic mesh deformation, the implementation of the GCL using this relationis not straightforward using the NLFD method. In the following, volumetric increments are alwaysconsidered as algebraic values which can either be positive or negative. The first approach developed by Tardiff et al. [8] is based on a linear representation of the volumetric increments relative to a reference time instance t . For any face m defined by itsvertices the induced volumetric change would simply be represented by drawing straight lines fromtheir initial position at t to their position at time instant t , see Figure 1.This approach has two advantages : first, it is easy to compute the volumetric increments ateach time instant using standard cell volume computational algorithms ; second, the volumetric increments due to each face are time periodic as long as the movement of the vertices is periodic.Once the volumetric increments are known for 2 N + 1 time instances defined by equation (12),their Fourier representations are calculated :Ω m ( t ) = N X k = − N ˆΩ m,k e ( i πk/T ) t , (37)and the Fourier formulations of the integrated face mesh velocities for each face m are introduced : G m ( t ) = N X k = − N ˆ G m,k e ( i πk/T ) t . (38)10 ( t ) r ( t ) r ( t ) r ( t ) r ( t ) r ( t ) r ( t ) r ( t ) Ω( t )Ω ( t ) Ω ( t )Ω ( t )Ω ( t ) Figure 1: Example of linear volumetric increments in 2D relatively to a reference time instant t Then by substituting, the Fourier representations into criteria (36), and exploiting the orthogo-nality of the Fourier basis, a system of 2 N + 1 equations (39) is obtained for each face m : i πkT ˆΩ m,k = ˆ G m,k for − N ≤ k ≤ N. (39)Therefore, the GCL are satisfied independently for each face of the control volume by computingthe Fourier coefficients ˆ G m,k and then applying an IDFT to transfer back the integrated facemesh velocities to the temporal domain. Despite its attractiveness, this method is restricted tolinear movements due to the manner in which the volumetric increments are computed. In general, the motion would not be linear and such representation of the volumetric increments will not besufficient to ensure the correct computation of the IFMV.Moreover the NLFD method is based on the assumption that the quantities are time periodicand can be expanded in Fourier series, but having a time periodic movement of the vertices doesnot guarantee time periodic volumetric increments but only that their temporal derivative will be periodic. This statement will be demonstrated through the following example.A 2D quadrilateral element is considered with the following motion defined by equation (40)11nd shown Figure 2 : α ( t ) = 2 πt, r = r , , r = r , , r = r , + R (1 − cos( α ( t ))) e x + R (sin( α ( t ))) e y , r = r , , (40)where the index 0 refers to the initial position of the grid, R is the radius defining the amplitude ofthe circular motion and e x and e y are the unit vectors in respectively the x and y directions. xy r , r , r , r , (a) xy r r r r (b) Figure 2: (a) The initial undeformed quadrilateral element is shown in black while the movement of the mesh pointsis presented in red with two deformed configurations of the cell in dashed lines, (b) Exact volumetric increment forthe face 2-3 in the x direction relatively to the initial configuration (in dash dot line) in blue For the face defined by the vertices r and r , the derivation of the expression of the exactvolumetric increment in the x direction and its time derivative leads to the following expressionsrespectively (41) and (42) :Ω ,x ( t ) = R α ( t ) − sin( α ( t ))) + Ry , − cos( α ( t ))) , (41) ∂ Ω ,x ( t ) ∂t = R ∂α∂t (1 − cos( α ( t ))) + Ry , ∂α∂t sin( α ( t )) , (42)where the length y , = ( r , · e y ). 12hus the time derivative of the volumetric increment is periodic whereas the volumetric incre- ment is the sum of a linear term and a periodic term and the direct application of the NLFDmethod on the exact volumetric increment is not possible since the linear term is not expandableas a Fourier serie. Additional work is required to ensure equation (36) is compliant with the NLFDmethod. In this section, we are going to demonstrate the following theorem 3.1,
Theorem 3.1.
Let Ω be a discretized control volume, enclosed by N f faces, and subjected to aperiodic motion of its vertices. Then given the knowledge of the exact volumetric increments Ω m for m = 1 , ..., N f , a sufficient condition to ensure the satisfaction of GCL in the NLFD frameworkis the computation of the integrated face mesh velocities, where the zeroth and higher modes can beexpressed as ˆ G m, = Ω m ( T ) T , (43)ˆ G m,k = i πkT ˆ p m,k for − N ≤ k ≤ N, k = 0 , (44) where ˆ G m,k and ˆ p m,k are the Fourier coefficients of respectively the integrated face mesh velocitiesand the periodic part of the exact volumetric increments given by, p m ( t ) = Ω m ( t ) − (cid:18) Ω m ( T ) T (cid:19) t. (45) Proof.
Under the assumption that the motion of the vertices is periodic, the temporal rate ofchange of the algebraic volume swept by each face through time is periodic. Thus the temporalderivative of the volumetric increments and the integrated face mesh velocities are periodic, theDFT is applied to the equation (36) leading to : G m ( t ) = ∂ Ω m ∂t = ˆ G m, + N X k = − N,k =0 ˆ G m,k e i πT kt , (46)13here ˆ G m,k , for − N ≤ k ≤ N are the Fourier coefficients of both the derivative of the volumetricincrement and the integrated face mesh velocity of a face m . By integrating the equation in time, the volumetric increment is expressed as :Ω m ( t ) = ˆ ∂ Ω m ∂t dt = ˆΩ m, + ˆ G m, t + N X k = − N,k =0 Ti πk ˆ G m,k e i πT kt . (47)where ˆΩ m, is a constant of integration. Then any volumetric increment can be interpreted as thesum of a linear term l m ( t ) and a periodic function p m ( t ) defined by : l m ( t ) = ˆ G m, t, (48) p m ( t ) = ˆΩ m, + N X k = − N,k =0 Ti πk ˆ G m,k e i πT kt . (49)Knowing the values of the volumetric increment at t = t and t = t + T , and exploiting theperiodicity of the function p m , yields :Ω m ( t ) = ˆ G m t + p m ( t )Ω m ( t + T ) = ˆ G m ( t + T ) + p m ( t + T ) p m ( t ) = p m ( t + T )ˆ G m, = Ω m ( t + T ) − Ω m ( t ) T (50)Hence the zeroth Fourier coefficients of the integrated face mesh velocities are known throughequation (50) applied for each face m and the linear part l m of the volumetric increments can becomputed at each instant. Then, an expression of the periodic part of any volumetric increment p m is obtained as : p m ( t ) = Ω m ( t ) − l m ( t ) = Ω m ( t ) − (cid:18) Ω m ( t + T ) − Ω m ( t ) T (cid:19) t. (51)Usually t would be taken as the initial time instant t = 0 corresponding to the undeformedconfiguration of the mesh, for this specific reference time instant Ω m (0) = 0, and the previousexpression can be further simplified into equation (52) : p m ( t ) = Ω m ( t ) − (cid:18) Ω m ( T ) T (cid:19) t. (52)14herefore, at each instant t the periodic part of the volumetric increments p m are known andwe introduce the Fourier coefficients for p m , noted as ˆ p m,k for − N ≤ k ≤ N . By calculatingthe temporal derivative in Fourier space of p m and exploiting the orthogonality of the Fourier basisfunctions, the rest of the Fourier coefficients of the integrated face mesh velocities ˆ G m,k are deducedfrom a system of 2 N equationsˆ G m,k = i πkT ˆ p m,k for − N ≤ k ≤ N, k = 0 . (53)Since the derivation in Fourier space puts to zero the contribution from the zeroth coefficient,the value of the integration constant ˆΩ m, is not relevant to compute the integrated face meshvelocities. (cid:4) Finally the procedure to compute the IFMV to enforce GCL by deducing the time derivative of the volumetric increments for each face is given by the pseudo-code (Algorithm 1). It is important tonote that since the values of the volumetric increments are required at t = T in order to deduce thezeroth Fourier coefficients ˆ G m, through equation (43), one additional time step is needed t N +1 compared to the number of time steps for the flow solver. However for this final time step theconfiguration of the mesh is the same as the initial (the undeformed mesh), thus no additional time step is needed for the mesh deformation. For this procedure, the key point is to compute the exactvolumetric increments as accurately as possible in order to preserve the spectral convergence of theNLFD method. In practice the accuracy of the previous method highly depends on the accuracy of the com- putation of the volumetric increments. For an hexahedral grid, as each face m sweeps throughthe computational domain, they form hexahedra between time intervals. The volume of any hex-ahedreon can be computed using a trilinear mapping [20] between the physical space and thecomputational domain Figure 3. This yields the following definitions. Definition 3.2.
The volume of any hexahedron as a function of the position vectors of the vertices or n = 0 , ..., N do Calculate the mesh deformation using the RBF for equally space time instances t n ; endfor f ace = 1 , ..., f ace max dofor n = 0 , ..., N + 1 do Calculate the volumetric increments Ω face ( t n ) ; end Deduce the zeroth Fourier coefficient via ˆ G face, = Ω face ( t N +1 ) T ; for n = 0 , ..., N do Extract the periodic part of the volumetric increments via p face ( t n ) = Ω face ( t n ) − ˆ G face, t n ; end Compute the Fourier coefficients ˆ p face,k via FFT on p face ( t ) ; for k = − N, ..., k = 1 , ...N do Deduce the k th Fourier coefficient of the integrated face mesh velocities viaˆ G face,k = i πkT ˆ p face,k ; end Compute the integrated face mesh velocities G face ( t ) via IFFT on ˆ G face,k with − N ≤ k ≤ N ; endAlgorithm 1: Pseudo-code representing the derived procedure to compute the integrated facemesh velocities and ensure GCL in the NLFD framework16 n the physical space r i for i = 1 , ..., , is evaluated through, Ω h = (Ω + Ω + Ω + Ω + Ω + Ω ) , with Ω ijkl = 112 ( r j + r k ) · (( r i + r j ) × ( r i + r l )) . (54) xyz ξηζ TT − Figure 3: Trilinear mapping between a hexahedron in the physical space and a reference cube in the computationaldomain
Definition 3.3.
For any face m of an hexahedral cell, the exact volumetric increment is estimatedthrough a sum of hexahedra each corresponding to the approximated volumetric increment betweentwo time samples t n − and t n and noted as Ω m,h ( t n ) , (see Figure 4) : Ω m ( t ) = 0 , Ω m ( t n ) = n X k =1 Ω m,h ( t k ) ! + ǫ Tm ( t n ) , for ≤ n ≤ N + 1 , (55) where t = 0 is the initial instant corresponding to the undeformed mesh and ǫ Tm ( t n ) is the truncationerror at time instant t n . Now that the mathematical tools to compute the volumetric increments are introduced, theaccuracy of the procedure presented in section 3.3 can be established, we have the first lemma 3.4,
Lemma 3.4.
In the context of theorem 3.1, and under the definitions 3.2 and 3.3, for any face m the temporal-order of accuracy of the zeroth Fourier coefficient of the integrated face mesh velocity ˆ G m, is limited to one. a ( t n − ) r b ( t n − ) r c ( t n − ) r d ( t n − ) r a ( t n ) r b ( t n ) r c ( t n ) r d ( t n ) (a) n − n n + 1Ω h ( t n ) Ω h ( t n +1 ) (b) Figure 4: (a) Approximated volumetric increment between two time steps t n − and t n : Ω h ( t n ) (b) Approximationof a volumetric increment as a sum of hexahedra Proof.
For each face m the path of the four corresponding vertices between two time steps islinearly approximated, as shown in Figure 2. Then, for any of these vertices r i , i = a, b, c, d at the n th time sample, we have the Taylor expansion : r i ( t n ) = r i ( t n − ) + (cid:18) ∂ r i ∂t ( t n − ) (cid:19) ( t n − t n − ) + 12 (cid:18) ∂ r i ∂t ( t n − ) (cid:19) ( t n − t n − ) + O (( t n − t n − ) ) | {z } Truncation error on the vertex path . (56)The volumetric increment is then defined by the vertices positions with the following indexation, r = r a ( t n − ) , r = r b ( t n − ) , r = r c ( t n − ) , r = r d ( t n − ) , r = r a ( t n ) , r = r b ( t n ) , r = r c ( t n ) , r = r d ( t n ) . (57)By substitution of the Taylor expansion of the vertex positions at instant t n into equation (54), andexploiting the linearity of the function, the error committed during the estimation of the volumetricincrement between two time steps is found to be of order two in τ = ( t n − t n − ) (see Appendix 7.1).Recalling that the number of time steps N ts = 2 N + 1 and using the definition of the time instance,the difference ( t n − t n − ) is written as τ = [ n − ( n − T N + 1 = TN ts . (58)Thus the truncation error during the estimation of the volumetric increment between two time steps18t the n th instant and noted ǫ Tm,h ( t n ) is of order two in τ and can be expanded as : ǫ Tm,h ( t ) = 0 ,ǫ Tm,h ( t n ) = E Tm ( t n − ) τ + O ( τ ), for 1 ≤ n ≤ N + 1 , (59)where E Tm ( t ) is a scalar periodic function depending on r i ( t ), and ∂ r i ∂t , for i = a, b, c, d (see Appendix7.1).In order to estimate the error committed on the exact volumetric increment approximated atthe n th time sample, these errors have to be summed, and yields : ǫ Tm ( t ) = 0 ,ǫ Tm ( t n ) = n X k =1 ǫ Tm,h ( t k ) = n X k =1 E Tm ( t k − ) ! τ + n O ( τ ), for 1 ≤ n ≤ N + 1 . (60)Hence, for 1 ≤ n ≤ N + 1 : | ǫ Tm ( t n ) | ≤ n (cid:18) max ≤ n ≤ N +1 |E Tm ( t n − ) | (cid:19) τ + n O ( τ ) = n O ( τ ) (61)Thus for n = N ts : | ǫ Tm ( t N ts ) | ≤ N ts O (cid:18) TN ts (cid:19) ! = O ( τ ) (62)The order of the error to approximate the exact volume of the volumetric increment may decrease over a period from 2 to 1 for the final value. Thus for any face m , the order of the truncation error ǫ Tm ( t N ts ) done to compute the zeroth Fourier coefficient of any integrated face mesh velocity ˆ G m, is one. (cid:4) Recalling that the zeroth Fourier coefficient is then used to extract the periodic part of any volumetric increment see theorem 3.1, the error committed on the rest of the Fourier coefficients ofthe integrated face mesh velocities is given by the following lemma 3.5,
Lemma 3.5.
In the context of theorem 3.1, and under the definitions 3.2 and 3.3, for any face m the temporal-order of accuracy of the Fourier coefficients ˆ G m,k , for − N ≤ k ≤ N with k = 0 , islimited to between one and two. roof. For 1 ≤ n ≤ N ts , the periodic part of the volumetric increment can be further expandedas : p m ( t n ) = Ω m ( t n ) − Ω m ( t N ts ) T t n = n X k =1 Ω m,h ( t k ) + ǫ Tm ( t n ) − N ts X k =1 Ω m,h ( t k ) + ǫ Tm ( t N ts ) T (cid:18) nTN ts (cid:19) From equation (60), we have : p m ( t n ) = n X k =1 Ω m,h ( t k ) + n X k =1 E Tm ( t k − ) ! τ + n O ( τ ) − N ts X k =1 Ω m,h ( t k ) + N ts X k =1 E Tm ( t k − ) ! τ + ( N ts ) O ( τ ) T (cid:18) nTN ts (cid:19) = n X k =1 Ω m,h ( t k ) − nN ts N ts X k =1 Ω m,h ( t k ) ! + " n X k =1 E Tm ( t k − ) − nN ts N ts X k =1 E Tm ( t k − ) τ + n O ( τ ) . Then the truncation error on the periodic part of any volumetric increment p m ( t ) is given for1 ≤ n ≤ N ts by : ǫ Tp m ( t n ) = " n X k =1 E Tm ( t k − ) − nN ts N ts X k =1 E Tm ( t k − ) τ + n O ( τ ) . (63)Since the bracketed term in equation (63) is dependent of n , the order of accuracy for any n isstill unclear. To refine the determination of the order of accuracy during the computation of p m ,the approximation of an integral using the Riemann sum is exploited.For any T -periodic function f at least three times continuous ( f ∈ C ([0; T ]), we have the20ollowing asymptotic development (64) where f ′ = ∂f∂t : R N ts = TN ts N ts − X k =0 f ( t k ) ,R N ts = ˆ T f ( t ) dt − T N ts ( f ( T ) − f (0)) + T N ts ) ( f ′ ( T ) − f ′ (0)) + O (cid:18) TN ts (cid:19) ! . (64)Then applying this result to the truncation error ǫ Tp m ( t n ), for 1 ≤ n ≤ N ts : ǫ Tp m ( t n ) = n − X k =0 E Tm ( t k ) − nT TN ts N ts − X k =0 E Tm ( t k ) τ + n O ( τ )= " n − X k =0 E Tm ( t k ) − nT ˆ T E Tm ( t ) dt − τ E Tm ( T ) − E Tm (0)) + τ
12 ( E Tm ′ ( T ) − E Tm ′ (0)) + O ( τ ) ! τ + n O ( τ )= " n − X k =0 E Tm ( t k ) − n (cid:18) hE Tm i T − τ T ( E Tm ( T ) − E Tm (0)) + τ T ( E Tm ′ ( T ) − E Tm ′ (0)) (cid:19) τ + n O ( τ )= " n − X k =0 E Tm ( t k ) − n hE Tm i T τ + n O ( τ )= " n − X k =0 (cid:0) E Tm ( t k ) − hE Tm i T (cid:1) τ + n O ( τ ) , (65)where h . i T represents the mean of a function on the segment [0; T ]. Taking advantage of the factthat the function ∆ E Tm = E Tm − hE Tm i T is T -periodic with zero mean value and exploiting a second21ime the expression (64), yields for n = N ts : ( N ts − X k =0 (cid:2)(cid:0) E Tm ( t k ) − hE Tm i T (cid:1)(cid:3)) τ = ˆ T (cid:0) ∆ E Tm (cid:1) dt | {z } =0 − τ (∆ E Tm ( T ) − ∆ E Tm (0)) + O ( τ )= − τ E Tm ( T ) − E Tm (0)) + O ( τ ) . (66)Substituting back the expression (66) into the final equation in (65) for n = N ts , leads to : ǫ Tp m ( t N ts ) = − τ ( E Tm ( T ) − E Tm (0)) + N ts O ( τ ) = O ( τ ) . (67)In summary the truncation error committed on the periodic part of any volumetric incrementfollows the equation : ǫ Tp m ( t ) = 0 ,ǫ Tp m ( t n ) = ( n − X k =0 (cid:2)(cid:0) E Tm ( t k ) − hE Tm i T (cid:1)(cid:3)) τ + n O ( τ ), for 1 ≤ n ≤ N.ǫ Tp m ( t N ts ) = O ( τ ) . (68)In general, the order of the truncation error on the approximation of the periodic part ofthe volumetric increment used as input for the NLFD method is of order between one and two. Analytically, we observe that this order is determined by the sum n − X k =0 (cid:2)(cid:0) E Tm ( t k ) − hE Tm i T (cid:1)(cid:3) , which isbounded for 1 ≤ n ≤ N ts by (cid:26) N ts max ≤ k ≤ N ts | (cid:0) E Tm ( t k ) − hE Tm i T (cid:1) | (cid:27) . This upper bound ensures thatin the worst case, the order of accuracy is 1. However asymptotically it is reasonable to assumethat for small and high values of n , the term n − X k =0 (cid:2)(cid:0) E Tm ( t k ) − hE Tm i T (cid:1)(cid:3) is small enough to considerthat the truncation error is of order 2 whereas for n in the middle of the range [1; N ts ], the order is greater than 1 but lesser than 2. (cid:4) Assuming that the spectral convergence of the Fourier transform is reached and taking advantageof its bijectivity, the truncation error on the Fourier coefficients ˆ G m,k and finally on the integratedface mesh velocities is of order between 1 and 2. Therefore the accuracy of the procedure is given by the following corollary 3.6 : 22 orollary 3.6. In the context of theorem 3.1, and under the definitions 3.2 and 3.3, for any face m the temporal-order of accuracy of the integrated face mesh velocities is limited to between oneand two. Thus it is important to note that even if the method described in section 3.1 enforced the
Geometric Conservation Law, the integrated face mesh velocities are determined within an accuracyof order 1 to 2. This is a disadvantage since the benefit of the spectral convergence of the NLFDmethod.
The computation of the metrics of a grid is often easier in a Cartesian grid, for this reasona mapping between the curvilinear physical space and a Cartesian computational space can beperformed. In this work, a trilinear mapping is already used to compute any hexahedron volume[20], but it can also be used to compute the time derivative of any hexahedron, its surface vectorsand the exact integrated face mesh velocities as long as the position and velocity vectors of the vertices are known. This section develops the derivation of these expressions.
Notation :
T ↔
Trilinear mapping p ↔ physical space : ( x, y, z ) r ↔ reference space : ( ξ, η, ζ ) m ↔ any faces of an hexahedron n ↔ normal vectorˆ n ↔ unit normal vector (69) Derivation :
The mapping T from the physical to the computational space is introduced : T = ( D C ) → ( D P )( ξ, η, ζ ) → ( x, y, z ) = ( T ( ξ, η, ζ )) , (70)where ( D C ) is the computational domain and ( D P ) is the physical domain. The application isdefined by considering a reference cube in the computational space which enables the mapping of r p = ( x, y, z ) in the physical space is mapped through r r =( x ( ξ, η, ζ ) , y ( ξ, η, ζ ) , z ( ξ, η, ζ )) based on the location vectors in the physical space r i ,p = [ x i , y i , z i ] i = 1 , ..., r r = (1 − ξ )(1 − η )(1 − ζ ) r ,p + ξ (1 − η )(1 − ζ ) r ,p + ξη (1 − ζ ) r ,p + (1 − ξ ) η (1 − ζ ) r ,p +(1 − ξ )(1 − η ) ζ r ,p + ξ (1 − η ) ζ r ,p + ξηζ r ,p + (1 − ξ ) ηζ r ,p , (71)where 0 ≤ ξ, η, ζ ≤ v p = ( v x , v y , v z ) in the physical domain is mapped in the same way v r =( v x ( ξ, η, ζ ) , v y ( ξ, η, ζ ) , v z ( ξ, η, ζ )) based on the velocity vectors of the vertices v i ,p = [ v xi , v yi , v zi ] i = 1 , ..., v r = (1 − ξ )(1 − η )(1 − ζ ) v ,p + ξ (1 − η )(1 − ζ ) v ,p + ξη (1 − ζ ) v ,p + (1 − ξ ) η (1 − ζ ) v ,p +(1 − ξ )(1 − η ) ζ v ,p + ξ (1 − η ) ζ v ,p + ξηζ v ,p + (1 − ξ ) ηζ v ,p . (72)For any face m of a cell, the normal vector is given by one of the following expressions : n r,ζ =0 = − (cid:18) ∂ r r ∂ξ (cid:19) x (cid:18) ∂ r r ∂η (cid:19) , n r,ζ =1 = + (cid:18) ∂ r r ∂ξ (cid:19) x (cid:18) ∂ r r ∂η (cid:19) , n r,ξ =0 = − (cid:18) ∂ r r ∂η (cid:19) x (cid:18) ∂ r r ∂ζ (cid:19) , n r,ξ =1 = + (cid:18) ∂ r r ∂η (cid:19) x (cid:18) ∂ r r ∂ζ (cid:19) , n r,η =0 = − (cid:18) ∂ r r ∂ζ (cid:19) x (cid:18) ∂ r r ∂ξ (cid:19) , n r,η =1 = + (cid:18) ∂ r r ∂ζ (cid:19) x (cid:18) ∂ r r ∂ξ (cid:19) , (73)where the sign is determined in order to have the normal pointing outward of the cell volume. TheJacobian matrix J ( ξ, η, ζ ) is expressed as : J ( ξ, η, ζ ) = (cid:18) ∂ r r ∂ξ ∂ r r ∂η ∂ r r ∂ζ (cid:19) (74)and its determinant can be calculated with one of the following expressions : | J | = (cid:18) ∂ r r ∂ξ (cid:19) · (cid:20)(cid:18) ∂ r r ∂η (cid:19) x (cid:18) ∂ r r ∂ζ (cid:19)(cid:21) = (cid:18) ∂ r r ∂η (cid:19) · (cid:20)(cid:18) ∂ r r ∂ζ (cid:19) x (cid:18) ∂ r r ∂ξ (cid:19)(cid:21) = (cid:18) ∂ r r ∂ζ (cid:19) · (cid:20)(cid:18) ∂ r r ∂ξ (cid:19) x (cid:18) ∂ r r ∂η (cid:19)(cid:21) . (75)Once the position vector, velocity vector, normal vectors and Jacobian are known, these quantities are used to compute the integrals of the volume and mesh velocity through a change of variables.24 olume integral Trough the application of the divergence theorem, the volume of the hexahedron can be evaluatedas such, V p = ˆ Ω p dV p = ‹ ∂ Ω p r p · dS p = 13 ‹ ∂ Ω p ( r p · ˆn p ) dS p We can then write the integral for the computational domain through the trilinear mapping toacquire, V p = 13 ‹ ∂ Ω p ( r p · ˆn p ) dS p = 13 ‹ ∂ Ω r ( r r · ˆn r ) | J r | dS r = 13 N f X m =1 ¨ ∂ Ω r,m ( r r,m · ˆn r,m ) | J r,m | dS r,m , where dS c,r,m is either dξdη , dηdζ or dζdξ and the integral boundaries are [0 1] . N.B. :
On anyface of the hexahedron only one of the variables in the reference space ξ , η or ζ has a fixed value.Thus the quantity ( r c,r,m . ˆn r,m ) | J r,m | is still a function of two variables which has to be integrated over the face.For each face, the computation of the integral over the surface under this form is not straight-forward (the difficulty comes from the unit normal vector) and needs to be simplified a priori . Thisis done by exploiting the relation (76), for the derivation of this expression see Appendix B in [23] : ˆn r,m | J r,m | = C r,m ˆN r,m , (76)where C r,m = C ( ξ, η, ζ ) is the cofactor matrix of the Jacobian matrix J for the trilinear mappingand ˆN r,m is the constant unit normal vector to the corresponding face in the reference space : ˆN ζ =0 = [ 0 0 − T , ˆN ζ =1 = [ 0 0 +1 ] T , ˆN η =0 = [ 0 − T , ˆN η =1 = [ 0 +1 0 ] T , ˆN ξ =0 = [ − T , ˆN ξ =1 = [ +1 0 0 ] T . (77)Once the equation (76) is substituted into the integrals over the surfaces, an explicit expression ofthe volume as a function of r i , i = 1 , ..., V c,p ) T = ( V + V + V + V + V + V ) T , (78)25here for any set i, j, k, l , the volumetric contribution of the face S ijkl is given by (79) :( V ijkl ) T = 13 ¨ ∂ Ω r,ijkl r r,ijkl · ( C r,ijkl ˆN r,ijkl ) dS r,ijkl = 112 ( r j + r k ) · (( r i + r j )x( r i + r l )) . (79) Time derivative of the volumetric integral
The temporal derivative of the volumetric integral can be expressed as, ∂V p ∂t = ∂∂t ˆ Ω p dV p . (80)By substituting the results of the previous section, primarily equations (78) and (79), an explicitexpression of the temporal derivative of the volume as a function of r i and v i , i = 1 , ..., (cid:18) ∂V p ∂t (cid:19) T = (cid:18) ∂V ∂t + ∂V ∂t + ∂V ∂t + ∂V ∂t + ∂V ∂t + ∂V ∂t (cid:19) T , (81)where for any set i, j, k, l , the time derivative volumetric contribution of the face S ijkl is given by(82) : (cid:18) ∂V ijkl ∂t (cid:19) T = 112 ( v j + v k ) · (( r i + r j )x( r i + r l ))+ 112 ( r j + r k ) · (( v i + v j )x( r i + r l ))+ 112 ( r j + r k ) · (( r i + r j )x( v i + v l )) . (82) Integrated face mesh velocities
The integral of the face mesh velocity in the physical domain for a face m is given by, G p,m = ¨ ∂ Ω c,p,m ( v · ˆn p , m ) dS c,p,m . By introducing the trilinear mapping, we can express the integrated face mesh velocities as, G p,m = ¨ ∂ Ω r,m v r,m · ( ˆn r,m | J r,m | ) dS r,m = ¨ ∂ Ω r,m v r,m · ( C r,m ˆN r,m ) dS r,m . Once the integration is performed, the explicit expressions of the integrated face mesh veloci-ties are obtained as a function of r i and v i , i = 1 , ...,
8. For a face with the set ( i, j, k, l ) ∈ } : v t = ( v i + v j + v k + v l ) S i,j,k,l = (( r i x r j ) + ( r j x r k ) + ( r k x r l ) + ( r l x r i )) S a,b,c = (( r a x r b ) + ( r b x r c ) + ( r c x r a )) for any set a, b, c ( G i,j,k,l ) T = 112 ( v t · S i,j,k,l + v j · S i,j,k + v k · S j,k,l + v l · S k,l,i + v i · S l,i,j ) (83)It was checked that with these expressions for the IFMV and the time derivative of the volumeas functions of velocity and position vectors of the vertices, the semi-discrete equation of the GCL (28) is analytically retrieved. In other words, the sum of equation (83) applied to the 6 sets { } is equal to expression (81). The methods presented in sections 3.2 and 3.3 to enforce the GCL are based on equation (36),and the integrated face mesh velocities are deduced from the calculation of the volumetric incrementsas input. The approach presented in this section using the trilinear mapping is quite differentbecause no volumetric increments are computed, the integrated face mesh velocities are directlyevaluated in physical time using equation (83). In addition the cell volumes are computed usingequation (54). Hence in the GCL equation as established in (25), both Ω and G = P N f m =1 G m areexactly calculated, the only degree of freedom remaining to enforce the equation is the discretizationof the temporal derivative operator (cid:18) ∂∂t (cid:19) . In the NLFD framework, this operator is discretizedin the Fourier domain and is a function of the number of harmonics N employed in the temporaldiscretization. Therefore the GCL equation will be satisfied if and only if the time derivative of thecell volume expressed in Fourier space converge to the Fourier time differentiation applied to thecell volume, DFT (cid:26)(cid:18) ∂ Ω ∂t (cid:19) T (cid:27) = (cid:18) ∂∂t (cid:19) F ourier (DFT { (Ω T ) } ) , (84)where T refers to the trilinear mapping. Hence, this method will not enforce the GCL for anynumber of time steps contrary to the method presented in section 3.3, but for a sufficient number of harmonics ensuring the convergence of the equations (84). Since this approach is based on the27xact integrated face mesh velocities and ensures the GCL with a spectral rate of convergence de-pending on the mesh motion, it provides a good alternative to the method exploiting the volumetricincrements with an order of accuracy comprised between one and two. In the section 4, we willpresent the numerical results of these different methods for several test cases. In this section, we extent the previous results to Time-Spectral (TS) method presented byGopinath and Jameson [5, 6]. Compared to the NLFD method which solves the governing equationsin the frequency domain, the Time-Spectral method solves the governing equations in the timedomain but exploits the features of a spectral approach.
Assuming a periodic flow and a periodic deformation of the mesh, we recall the temporal dis-cretization of the modified state vector ¯ w = Ω w equations (8) and (10),¯ w ( t ) = N X k = − N ˆw k e i (2 πk/T ) t , with : ˆw k = 12 N + 1 N X n =0 Ω( t n ) w ( t n )e − i (2 πk/T ) t n , where T is the time period, N is the number of modes considered in the DFT and t n the equallyspaced time instances given by, t n = n N + 1 T , for n = 0 , .., N. In Fourier space, the time discretization operator leads to, ∂ ¯w ∂t ( t ) = N X k = − N i πkT ˆw k e i (2 πk/T ) t , (85) ⇔ ∂ ¯w ∂t ( t ) = 2 πT N X k = − N ik N + 1 N X K =0 Ω( t K ) w ( t K )e − i (2 πk/T ) t K ! e i (2 πk/T ) t . (86)By evaluating this expression for each time instance t n , we have for n = 0 , ..., N , ∂ ¯w ∂t ( t n ) = 2 πT N X k = − N ik N + 1 N X K =0 Ω( t K ) w ( t K )e − i (2 πk/T ) t K ! e i (2 πk/T ) t n , (87)28 ∂ ¯w ∂t ( t n ) = N X K =0 " Ω( t K ) w ( t K ) πT N + 1 N X k = − N ik e i (2 πk )( n − K ) / (2 N +1) ! . (88)We introduce, the coefficients d n,K , defined for n = 0 , ..., N by, d n,K = 2 πT N + 1 N X k = − N ik e i (2 πk )( n − K ) / (2 N +1) , (89)the compact form of the coefficients for an odd number of time steps is written as follows (for thederivation see Reference [6]), d n,K = πT
12 ( − n − K csc (cid:18) π ( n − K )2 N + 1 (cid:19) , if K = n , if K = n, (90)and, ∂ ¯w ∂t ( t n ) = N X K =0 d n,K ¯ w ( t K ) . (91)The temporal-derivation operator appears as the multiplication of a matrix D = ( d n,K ) ≤ n,K ≤ N with each vector ( ¯w i ( t K )) ≤ K ≤ N , for i = 1 , ..., i refers to the component ofthe modified state vector in the governing equations. In addition, this matrix is skew-symmetric,independent of any state variables and completely determined by the number of harmonics usedin the DFT and the temporal period. Then a pseudo-time t ∗ is introduced and the equations aresolved in the time domain through, ∂ ¯w ∂t ∗ ( t n ) + ∂ ¯w ∂t ( t n ) + R ( w ( t n )) = 0 , for n = 0 , ..., N. (92) Recall that in order to obtain a consistent solution method, the GCL must be discretized usingthe same numerical scheme employed to discretize the governing equations. In the case of Time-Spectral method, it leads to the following theorem 3.7,
Theorem 3.7.
Let Ω be a discretized control volume, enclosed by N f faces, and subjected to aperiodic motion of its vertices. Then given the knowledge of the exact volumetric increments Ω m for m = 1 , ..., N f , a sufficient condition to ensure the satisfaction of GCL in the TS framework isthe computation of the integrated face mesh velocities through the following relations, G m = ( D ) p m + h G m i T ( I N ) , (93)29 here for all m , h G m i T are the temporal mean values of the integrated face mesh velocities, I N isthe identity matrix of dimension N + 1 , G m = ( G m ( t n )) ≤ n ≤ N and p m = ( p m ( t n )) ≤ n ≤ N arethe vectors grouping the time instances of respectively the integrated face mesh velocities and theperiodic part of the exact volumetric increments given by, p m ( t ) = Ω m ( t ) − (cid:18) Ω m ( T ) T (cid:19) t, (94) and D = ( d n,K ) ≤ n,K ≤ N is the matrix representing the temporal derivation operator of the Time-Spectral method, defined by its coefficients d n,K for ≤ n, K ≤ N , d n,K = πT ( − n − K csc (cid:18) π ( n − K )2 N + 1 (cid:19) , if K = n , if K = n. (95) Proof.
Under the assumption that the motion of the vertices is periodic, the temporal rate ofchange of the algebraic volume swept by each face through time is periodic. Thus the temporalderivative of the volumetric increments and the integrated face mesh velocities are periodic, theDFT is applied to equation (36) leading to : G m ( t ) = ∂ Ω m ∂t = ˆ G m, + N X k = − N,k =0 ˆ G m,k e i πT kt , (96)where for any face m , ˆ G m,k are the Fourier coefficients of both the derivative of the volumetricincrement and the integrated face mesh velocity. The mean of a function expandable in Fourierserie is given by its zeroth Fourier coefficient thus, h G m i T = ˆ G m, . (97)From the proof of theorem 3.1, we have the following relations,Ω m ( t ) = l m ( t ) + p m ( t ) , (98) l m ( t ) = ˆ G m, t = Ω m ( T ) T t, (99) p m ( t ) = ˆΩ m, + N X k = − N,k =0 Ti πk ˆ G m,k e i πT kt , (100)30 m ( t ) = Ω m ( t ) − (cid:18) Ω m ( T ) T (cid:19) t. (101)Then, by exploiting these results and applying the Time-Spectral temporal derivation to the periodicpart of the volumetric increments, we can write for each time instance t n , with n = 0 , ..., N , G m ( t n ) = ∂ Ω m ∂t ( t n )= ∂ ( l m + p m ) ∂t ( t n )= h G m i T + ∂p m ∂t ( t n )= h G m i T + N X K =0 d n,K p m ( t K ) , (102)Finally, if we group all the time instances in a vector G m , we obtain, G m = ( D ) p m + h G m i T ( I N ) , (103)where I N is the identity matrix of dimension 2 N + 1 and D = ( d n,K ) ≤ n,K ≤ N is the matrixrepresenting the temporal derivation operator of the Time-Spectral method. The condition given byequation (103) is a criteria to ensure that the GCL are enforced in the Time-Spectral framework. (cid:4)
4. Numerical results
The new approaches to enforce the Geometric Conservation Law developed in section 3 arenumerically tested in order to validate their procedures. The protocol, test cases and results arepresented in the following sections.
The physical interpretation of the GCL is that any uniform flow must be preserved by thenumerical scheme employed for the flow solver and independently of the mesh movements. Thislaw imposes constraints on the manner to compute some geometrical quantities such as the volume31nd the integrated face mesh velocities. Thus the first step of our test is to ensure the preservationof uniform flow by computing the relative error between the initially defined uniform state vector W and the computed state vector W by the flow solver, RelErr = max ≤ n ≤ N (cid:26) max ≤ n v ≤ N cell (cid:18) max ≤ j ≤ (cid:12)(cid:12)(cid:12)(cid:12) W j ( n v , t n ) − W ,j ( n v , t n ) W ,j ( n v , t n ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:19)(cid:27) , (104)where W = ρ , W = ρu , W = ρv , W = ρw and W = ρE and n v is the index pointing to thegrid cell with N cell the number of cells in the mesh.However the verification of uniform flow preservation only guarantees that the GCL are satisfied"by summing over the faces", but not that the computed integrated face mesh velocities are correct.Indeed as long as the sum of the time derivative of the volumetric increments is equal to the timederivative of the cell volume, ∂ Ω ∂t = N f X m =1 ∂ Ω m ∂t , (105)the deduced integrated face mesh velocities from the time derivative of the volumetric incrementsfrom equation (36) enforce the GCL after the summation through the faces (see equation (28)) butthe integrated face mesh velocities themselves may not converge to the correct value. Thus in order to verify that the GCL are enforced with a correct evaluation of the integratedface mesh velocities, the values derived from the trilinear mapping equations (81) and (83) based onthe location and velocity vectors of the grid points retrieved from the dynamic mesh deformation,are considered as reference. Therefore for each motion of the mesh and for various number ofharmonics N , four different implementations of the integrated face mesh velocities are compared :
1. the IFMV deduced from the linear volumetric increments from Tradiff et al. [8] see Figure 1noted as "NLFD-LVI" ;2. the IFMV calculated with the new method based on the exact volumetric increments approx-imated as a sum of hexahedron see Figure 4 noted as "NLFD-AEVI" ;3. the previous approximation obtained by taking the average of the velocity of the four vertices defining a face and projected along the surface normal vector noted as "AVG";4. the method based on the trilinear mapping noted as "TRI-MAP" and used as reference forthe exact values of the IFMV.For each of these approaches, the preservation of uniform flow is tested. Then different quantitiesare compared by computing the maximum absolute error : comparison of the sum of the IFMV to the NLFD time derivative of the cell volume computedusing the numerical scheme of the flow solver. This comparison is similar to a demonstrationof the preservation of uniform flow: AbsErr = max ≤ n ≤ N max ≤ n v ≤ N cell (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) N f X m =1 G m ( n v , t n ) MET HOD − (cid:18) ∂ Ω ∂t ( n v , t n ) (cid:19) NLF D (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ;(106) • comparison of the IFMV to the reference integrated face mesh velocities (TRI-MAP) in eachdirection dir = x , y or z : AbsErr = max ≤ n ≤ N (cid:26) max ≤ n v ≤ N cell (cid:12)(cid:12)(cid:12) ( G m,dir ( n v , t n )) MET HOD − ( G m,dir ( n v , t n )) T RI − MAP (cid:12)(cid:12)(cid:12)(cid:27) . (107) This section presents the different mesh motions impose as test cases. The time period is alwaystaken to be unity. All tests are performed on a square mesh of size 10 × ×
10, and of lengths L x = 3 . L y = 2 .
8, and L z = 2 .
4. The undeformed positions of the mesh are indexed with thesubscript 0, if needed the RBF points are indexed with the subscript r . The two parameters in the JST scheme are κ (2) = 1 and κ (4) = 1 /
32. The simulations are run for a number of harmonics, N from 1 to 20. The mesh deformations for cases 2, 4 and 5 at an arbitrary time instant are presentedon Figure 5. Three test cases are performed by directly imposing the mesh deformation to the entire mesh.
The velocity of the vertices is computed based on the analytic time derivation of the vector positionof the vertices. For any vertex, its initial position is noted ( x , y , z ). The parameters A x , A y , A z , R , and α can be arbitrarily chosen as long as no degenerative cells (cells with negative volume)appear during the motion. The analytic functions employed for the motions are as follows : Case 1 : x ( t ) = x + A x sin (cid:16) πx L x (cid:17) sin (cid:16) πy L y (cid:17) sin (cid:16) πz L z (cid:17) sin (2 πt ) y ( t ) = y + A y sin (cid:16) πx L x (cid:17) sin (cid:16) πy L y (cid:17) sin (cid:16) πz L z (cid:17) sin (2 πt ) z ( t ) = z + A z sin (cid:16) πx L x (cid:17) sin (cid:16) πy L y (cid:17) sin (cid:16) πz L z (cid:17) sin (2 πt ) (108)33 ase 2 :
2D perturbation of the mesh with a non linear motion; however the time-average volumeswept by a face, ˆ G m, = 0 in (46). For any cell the projection of the motion along a plane z = constant is shown in Figure 8 : α ( t ) = α sin(2 πt ) x ( t ) = x + y cos( π − α ( t )) y ( t ) = y sin( π − α ( t )) z ( t ) = z (109) Case 3 :
2D perturbation of the mesh with a non linear motion and ˆ G m, = 0 in (46), the defor-mation is prescribed only for the interior grid points while the boundary points are fixed. Theprojection of the motion along a plane z = constant is identical to the movement of the 3rdnode on the Figure 2 presented in section 3.2 : α ( t ) = 2 πtx ( t ) = x + R (1 − cos( α ( t )) y ( t ) = y + R sin( α ( t )) z ( t ) = z (110) Two test cases are performed by deforming the mesh through the RBF. The analytic functionsemployed for the RBF motions are as follows :
Case 4 :
3D perturbation of the mesh using the RBF, with each point having its own linear motion(amplitude and direction) : s x ( t ) = r x ( x ,r , y ,r , z ,r ) sin(2 πy ,r ) sin(2 πz ,r ) sin(2 πt ) s y ( t ) = r y ( x ,r , y ,r , z ,r ) sin(2 πx ,r ) sin(2 πz ,r ) sin(2 πt ) s z ( t ) = r z ( x ,r , y ,r , z ,r ) sin(2 πy ,r ) sin(2 πz ,r ) sin(2 πt )where r x ( x ,r , y ,r , z ,r ); r y ( x ,r , y ,r , z ,r ) and r z ( x ,r , y ,r , z ,r ) are randomly generated (111)34 ase 5 : simulation of a sinusoidal pitching motion : α ( t ) = α cos(2 πt ) x p = 0 . L x s x ( t ) = ( x ,r − x p )[cos( α ( t )) −
1] + y ,r sin( α ( t )) s y ( t ) = − ( x ,r − x p ) sin( α ( t )) + y ,r [cos( α ( t )) − s z ( t ) = z ,r (112) (a) Case 2 (b) Case 4(c) Case 5 Figure 5: Mesh deformations of the exterior grid points for cases 2, 4 and 5 at an arbitrary chosen time step
The results demonstrating uniform flow preservation are shown for all test cases in Figure 6.The evolution of the relative error defined by equation (104) is presented as a function of the number of time steps N ts .The results show that the two methods employing the IFMV deduced from the Fourier discretiza-tion preserve uniform flow, while the approximation derived from the AVG yields the least accurateresults. This is consistent since for both methods NLFD-LVI and NLFD-AEVI, despite different35 a) Case 1 (b) Case 2(c) Case 3 (d) Case 4(e) Case 5 Figure 6: Relative error regarding the uniform flow preservation for each test case of the volumetric increments is equal to the temporal derivative of the cell volume evaluated in thefrequency domain (equation (105)).It is also observed that using the (TRI-MAP) integrated face mesh velocities preserves uniformflow and thus satisfies the GCL given a sufficient number of harmonics (see cases 2, 4 and 5) whichis expected. Its rate of convergence should be exactly the same as the rate of convergence of the time derivative of the cell volume in the Fourier space. This is verified in the next section.
The results are shown on Figures 7 through 12. It is important to note that for all figures, thegraph (a) refers to equation (106) as the function of the number of time steps and is not the sumof the graphs from (b), (c) and (d) which refer to equation (107). The errors that appear on the y -axis of the figures are the max norm between the investigated approaches, both NLFD-based andAVG and the reference approach (TRI-MAP).Regarding the comparison of the sum of the integrated face mesh velocities (IFMV) to theNLFD temporal derivative of the volume from Figures 7(a) through Figure 12(a), the results showthat the sum of the IFMV computed with the methods NLFD-LVI, NLFD-AEVI and TRI-MAP converge to the expected value for all cases while the AVG method provides the correct value onlyfor cases 1 through 3 and yields a constant absolute error above 10 − for cases 4 and 5. Recall thatthe maximum error in the sum of the IFMV is a measure of the level to which GCL is satisfied asgiven in the semi-discrete GCL equation (28). Hence the NLFD-based approaches prove to satisfythe GCL for all considered grid deformation and for any number of harmonics which is expected by design. The reference approach (TRI-MAP) satisfies this requirement exactly for linear deformationcases as shown for Cases 1 (Figure 7(a)) and 4 (Figure 11(a)) for any number of harmonics andconverged spectrally for nonlinear deformation cases (Cases 2, 3, and 5). The spectral rate ofconvergence is observed compared to the first-order backward finite-difference (∆ t (1) ) and second-order centered finite-difference (∆ t (2) ) approximating the time derivative of the cell volume. As expected, this rate of convergence is found to be similar for the preservation of uniform flow usingthe reference TRI-MAP method. However, the AVG approach is not designed to enforce the GCL,it is only an approximation based on the mesh velocities and face metrics and hence for the casesconsidered in this article, the method proved to ensure the GCL with an accuracy up to 10 − .37 a) (b)(c) (d) Figure 7: Case 1 : (a) Comparison of the sum of the integrated face mesh velocities to the NLFD time derivativeof the volume (b) Comparison of the individual integrated face mesh velocity to the reference values (TRI-MAP) inthe x direction (c) in the y direction (d) in the z direction
38 comparison of the individual integrated face mesh velocities for each direction reveals the limits and provides interesting insights of the investigated approaches. Two primary observations can bemade. First, the NLFD based approaches converge at most at second order as expected based oncorollary 3.6, if the mesh deformation along the observed direction is nonlinear. For Cases 2, 3, and5, the mesh deformation in both the x -and y -directions are nonlinear as shown in sub-figures (b) and(c) of Figures 9, 10, and 12. One exception is the spectral rate of convergence for the y -direction in Case 2. These results can be explained by analyzing in details the mesh movement. Since themotion is in two dimensions, let us consider a constant z plane, then the deformation of any cellcan be represented as shown in Figure 8. We observe that in the y -direction, the area swept by the xy Figure 8: Two dimensional projection of the motion in case 2 for one cell : the exact volumetric increment in the x direction is filled in clear green and in dark blue in the y direction. The blue dashed dot arrows show the paths ofthe vertices. faces can be exactly evaluated using a linear approximation of the curved boundaries shown in blue.Therefore in the y -direction, the volumetric increments are exactly computed and the individual IFMV are correctly computed using either the LVI or AEVI methods once the temporal derivativeoperator is converged in Fourier space. In the x -direction, a linear approximation is insufficient tocompute exactly the volumetric increments thus the AEVI method converges at an order betweenone and two as stated in corollary 3.6.Second, even if the numerical scheme enforces the GCL by preserving uniform flow, the employed method may not converge to the correct integrated face mesh velocities. The method based on the39pproximation of the exact volumetric increment (AEVI) is found to be converging toward thereference values at an order between one and two in the worst test cases considered here (4 & 5).This is consistent with the derivation of the error from section 3.3.2 and the resulting corollary (seecorollary 3.6). The NLFD-LVI and AVG methods may present significant inaccuracies depending on the mesh deformation. (a) (b)(c) (d) Figure 9: Case 2 : (a) Comparison of the sum of the integrated face mesh velocities to the NLFD time derivativeof the volume (b) Comparison of the individual integrated face mesh velocity to the values (TRI-MAP) in the x direction (c) in the y direction (d) in the z direction a) (b)(c) (d) Figure 10: Case 3 : (a) Comparison of the sum of the integrated face mesh velocities to the NLFD time derivativeof the volume (b) Comparison of the individual integrated face mesh velocity to the values (TRI-MAP) in the x direction (c) in the y direction (d) in the z direction a) (b)(c) (d) Figure 11: Case 4 : (a) Comparison of the sum of the integrated face mesh velocities to the NLFD time derivativeof the volume (b) Comparison of the individual integrated face mesh velocity to the values (TRI-MAP) in the x direction (c) in the y direction (d) in the z direction a) (b)(c) (d) Figure 12: Case 5 : (a) Comparison of the sum of the integrated face mesh velocities to the NLFD time derivativeof the volume (b) Comparison of the individual integrated face mesh velocity to the values (TRI-MAP) in the x direction (c) in the y direction (d) in the z direction .5. Time-Spectral Method The numerical results for the Time-Spectral method are the same as that shown for NLFD-LVIand NLFD-AEVI depending on which approach is retained to compute the volumetric increments.For this reason, the graphs are not reproduced in this article. The comparisons and conclusions derived for the NLFD approach hold for the Time-Spectral method as well.
5. Discussion and Conclusion
The limits of the previous method of Tardiff et al. [8] (NLFD-LVI) were clarified and demon-strated numerically and a modified approach (NLFD-AEVI) has been presented that ensures thesatisfaction of the Geometric Conservation Law for a flow solver based on either the NLFD or Time-
Spectral discretization of the ALE formulation of the Euler equations. The methods NLFD-AEVIand NLFD-LVI aim to satisfy the GCL by computing the integrated face mesh velocities accordingto the numerical discretization of the flow solver and take as input the face volumetric increments.The accuracy of the methods was shown to be highly dependent on the computation of the correctvolumetric increments and in the worst cases considered converged at first-to-second-order for the
NLFD-AEVI approach (corollary 3.6) or zeroth-order for the NLFD-LVI procedure. The integratedface mesh velocities themselves may not converge to the correct value as demonstrated in our nu-merical test. Although the approaches have been verified to preserve uniform flow for any numberof harmonics; such a low order of accuracy defeats the purpose of spectral in time methods. Hencean alternate novel approach has been developed based on a trilinear mapping between the physi- cal domain and the computational space which allows the evaluation of the exact cell volume andintegrated face mesh velocities. The disadvantage of this method is that it is not consistent withthe discretization of the flow solver, meaning that freestream preservation is not satisfied for anynumber of harmonics as it is with the modified approach, NLFD-AEVI. However such inconvenienceis compensated by its spectral rate of convergence, which is sufficient to ensure the satisfaction of the GCL and preserve uniform flow.
6. Acknowledgment
We would like to gratefully acknowledge the financial support of the Natural Sciences andEngineering Research Council of Canada and McGill University.44 . Appendix
Lemma 7.1.
In the context of the theorem 3.1, and under the definitions 3.2 and 3.3, for any face m the truncation error ǫ Tm,h on the volumetric increment between two time steps t n − and t n is oforder two and can be written as, ǫ Tm,h ( t n ) = E Tm ( t n − ) τ + O ( τ ) . (113) where E Tm is a scalar periodic function and τ = t n − t n − . Proof.
We introduce the scalar triple product application L defined by, L = R × R × R → R ( v , v , v ) → v · ( v × v ) = det( v , v , v ) (114)due to the properties of the determinant this application is a 3-linear alternating form meaning thatif any of the three vector is a linear combination of the two others the result is zero.Recalling that the volume of any hexahedra is computed as a function of the position vectorsof the vertices in the physical space r i for i = 1 , ..., Ω h = (Ω + Ω + Ω + Ω + Ω + Ω ) , with Ω ijkl = 112 L ( r j + r k , r i + r j , r i + r l ) . (115)Let m be a face defined by the vertices r i , i = a, b, c, d (see Figure 2), at the n th time sample aTaylor expansion gives : r i ( t n ) = l i ( t n ) + ǫ i ( t n ) , Linear approximation : l i ( t n ) = r i ( t n − ) + (cid:18) ∂ r i ∂t ( t n − ) (cid:19) τ, Truncation error : ǫ i ( t n ) = 12 (cid:18) ∂ r i ∂t ( t n − ) (cid:19) τ + O ( τ ) . (116)The volumetric increment between the two time instances t n − and t n is approximated by ahexahedra defined using the vertex positions with the following indexation, r = r a ( t n − ) , r = r b ( t n − ) , r = r c ( t n − ) , r = r d ( t n − ) , r = r a ( t n ) , r = r b ( t n ) , r = r c ( t n ) , r = r d ( t n ) . (117)45y substituting the Taylor expansions into the vertex positions to compute the volume of thevolumetric increment through equation (115), and then by exploiting the 3-linearity of the tripleproduct application the order of the truncation error is evaluated. The lowest order terms of thetruncation error are given by one of the following generic forms : L ( r i ( t n − ) , r j ( t n − ) , ǫ k ( t n )) = (cid:26) r i ( t n − ) · (cid:20) r j ( t n − ) × (cid:18) ∂ r k ∂t ( t n − ) (cid:19)(cid:21)(cid:27) τ + O ( τ ) , L ( r i ( t n − ) , ǫ j ( t n ) , r k ( t n − )) = (cid:26) r i ( t n − ) · (cid:20) (cid:18) ∂ r j ∂t ( t n − ) (cid:19) × r k ( t n − ) (cid:21)(cid:27) τ + O ( τ ) , L ( ǫ i ( t n ) , r j ( t n − ) , r k ( t n − )) = (cid:26) (cid:18) ∂ r i ∂t ( t n − ) (cid:19) · [ r j ( t n − ) × r k ( t n − )] (cid:27) τ + O ( τ ) . (118)where i, j and k are the vertices indices. Therefore for any face m , the truncation error ǫ Tm,h on the volumetric increment between twotime steps t n − and t n is of order two. In addition, it is possible to write the lowest order term ofthe error as a linear combination of the previous forms equation (118), thus there exists a scalarfunction E Tm depending on the vertices paths r i and their second temporal derivatives ∂ r i ∂t suchthat, ǫ Tm,h ( t n ) = E Tm ( t n − ) τ + O ( τ ) . (119)Due to the temporal periodicity of the vertices paths r i , the function E Tm is also periodic. (cid:4) ReferencesReferences [1] K. C. Hall, J. P. Thomas, W. S. Clark, Computation of unsteady nonlinear flows in cascades us- ing a harmonic balance technique, AIAA Journal 40 (5) (2002) 879–886. doi:10.2514/2.1754 .[2] M. McMullen, A. Jameson, J. Alonso, Application of a non-linear frequency domainsolver to the euler and navier-stokes equations, in: 40th AIAA Aerospace Sciences Meet-ing & Exhibit, American Institute of Aeronautics and Astronautics, Reno, Nevada, 2002. doi:10.2514/6.2002-120 . doi:10.2514/1.15127 .[4] M. S. McMullen, A. Jameson, The computational efficiency of non-linear fre-quency domain methods, Journal of Computational Physics 212 (2) (2006) 637–661. doi:10.1016/j.jcp.2005.07.021 . [5] A. Gopinath, A. Jameson, Time spectral method for periodic unsteady computationsover two- and three- dimensional bodies, in: 43rd AIAA Aerospace Sciences Meetingand Exhibit, American Institute of Aeronautics and Astronautics, Reno, Nevada, 2005. doi:10.2514/6.2005-1220 .[6] A. Gopinath, A. Jameson, Application of the time spectral method to periodic unsteady vortex shedding, in: 44th AIAA Aerospace Sciences Meeting and Exhibit, American Institute ofAeronautics and Astronautics, Reno, Nevada, 2006. doi:10.2514/6.2006-449 .[7] F. Kachra, S. Nadarajah, Aeroelastic solutions using the time accurate and non-linear frequencydomain methods, in: 44th AIAA Aerospace Sciences Meeting and Exhibit, American Instituteof Aeronautics and Astronautics, Reno, Nevada, 2006. doi:10.2514/6.2006-445 . [8] P.-O. Tardif, S. Nadarajah, Three-dimensional aeroelastic solutions via the nonlinear frequency-domain method, AIAA Journal 55 (10) (2017) 3553–3569. doi:10.2514/1.J054849 .[9] P. D. Thomas, C. K. Lombard, Geometric Conservation Law and Its Application to Flow Com-putations on Moving Grids, AIAA Journal 17 (10) (1979) 1030–1037. doi:10.2514/3.61273 .[10] M. Lesoinne, C. Farhat, Geometric conservation laws for flow problems with mov- ing boundaries and deformable meshes, and their impact on aeroelastic computations,Computer Methods in Applied Mechanics and Engineering 134 (1-2) (1996) 71–90. doi:10.1016/0045-7825(96)01028-6 .[11] H. Guillard, C. Farhat, On the significance of the geometric conservation law for flow computa-tions on moving meshes, Computer Methods in Applied Mechanics and Engineering 190 (11-12) (2000) 1467–1482. doi:10.1016/S0045-7825(00)00173-0 .4712] C. Farhat, P. Geuzaine, C. Grandmont, The Discrete Geometric Conservation Law and theNonlinear Stability of ALE Schemes for the Solution of Flow Problems on Moving Grids,Journal of Computational Physics 174 (2) (2001) 669–694. doi:10.1006/JCPH.2001.6932 .[13] D. J. Mavriplis, Z. Yang, Construction of the discrete geometric conservation law for high- order time-accurate simulations on dynamic meshes, Journal of Computational Physics 213 (2)(2006) 557–573. doi:10.1016/j.jcp.2005.08.018 .[14] J. Blajek, Computational Fluid Dynamics: Principles and Applications, Elsevier, 1988.[15] A. Jameson, W. Schmidt, E. Turkel, Numerical solution of the Euler equations by finite vol-ume methods using Runge Kutta time stepping schemes, in: 14th Fluid and Plasma Dy- namics Conference, American Institute of Aeronautics and Astronautics, Reno, Nevada, 1981. doi:10.2514/6.1981-1259 .[16] A. Jameson, Analysis and design of numerical schemes for gas dynamics, 1 : artificial diffusion,upwind biasing, and their effect on accuracy and multigrid convergence, International Journalof Computational Fluid Dynamics 4 (3) (1995) 171–218. doi:10.1080/10618569508904524 . [17] H. Wendland, Piecewise polynomial, positive definite and compactly supported radial func-tions of minimal degree, Advances in Computational Mathematics 4 (1) (1995) 389–396. doi:10.1007/BF02123482 .[18] H. Zhang, M. Reggio, J. Y. Trépanier, R. Camarero, Discrete form of the GCL for movingmeshes and its implementation in CFD schemes, Computers and Fluids 22 (1) (1993) 9–23. doi:10.1016/0045-7930(93)90003-R .[19] D. J. Mavriplis, C. R. Nastase, On the geometric conservation law for high-order discontinuousGalerkin discretizations on dynamically deforming meshes, Journal of Computational Physics230 (11) (2011) 4285–4300. doi:10.1016/j.jcp.2011.01.022 .[20] J. K. Dukowicz, Efficient volume computation for three-dimensional hexahedral cells, Journal of Computational Physics 74 (2) (1988) 493–496. doi:10.1016/0021-9991(88)90091-5 .[21] J. I. López, M. Brovka, J. M. Escobar, R. Montenegro, G. V. Socorro, Strategies for optimiza-tion of hexahedral meshes and their comparative study, Engineering with Computers 33 (1)(2017) 33–43. doi:10.1007/s00366-016-0454-1 .4822] P. Knabner, S. Korotov, G. Summ, Conditions for the invertibility of the isoparametric mapping for hexahedral finite elements, Finite Elements in Analysis and Design 40 (2) (2003) 159–172. doi:10.1016/S0168-874X(02)00196-8 .[23] P. Zwanenburg, S. Nadarajah, Equivalence between the Energy Stable Flux Reconstructionand Filtered Discontinuous Galerkin Schemes, Journal of Computational Physics 306 (2016)343–369. doi:10.1016/j.jcp.2015.11.036 .440