A Conservative High-Order Method Utilizing Dynamic Transfinite Mortar Elements for Flow Simulation on Curved Sliding Meshes
AA Conservative High-Order Method Utilizing Dynamic Transfinite MortarElements for Flow Simulation on Curved Sliding Meshes
Bin Zhang ∗ , Chunlei Liang Department of Mechanical and Aerospace Engineering,The George Washington University, Washington, DC 20052 USA
Abstract
We present a high-order method for flow simulation on unstructured curved nonconforming sliding meshes.This method utilizes dynamic transfinite mortar elements to exchange flow information between the twosides of a sliding interface. The method is arbitrarily high-order accurate in space, provably conservative,and satisfies outflow condition. Moreover, it retains the accuracy of a time marching scheme, and thusallows substantial reduction of rotational speed effects when equipped with a high-order temporal scheme.The method’s capability of simultaneously handling multiple rotational objects is also explored. Details onthe implementation are provided as well.
Keywords: sliding mesh, mortar elements, flux reconstruction, spectral difference, high-order methods
1. Introduction
Sliding mesh has been widely used in flow simulations about moving objects. For example, it is an idealchoice for handling rotational geometries such as stirred tanks [1] and helicopter rotor blades [2]. It canalso be used to ensure good mesh qualities in circumstances where purely deforming mesh may otherwise bevery skewed, such as for simulating oscillating wings [3] and vortex-induced-vibration devices [4]. In manyapplications, a sliding mesh method has advantages over other methods such as overset mesh methods [5]and immersed boundary methods [6] for its simplicity, efficiency and accuracy. But so far, most sliding meshmethods are still limited to low-order (second order and below) schemes that are unfavorable for simulatingvortex dominated flows.Tremendous progress has been made on high-order methods in the past decades in the computationalfluid dynamics community [7]. For instance, such methods include the discontinuous Galerkin (DG) method[8, 9], the spectral element method [10–12], the spectral volume method [13, 14], the spectral difference (SD)method [15–19], to name just a few. Among these methods, the SD method solves equations in differentialform directly, and is one of the most efficient high-order methods. Recently, the ideas of collocating solutionand flux points of the SD method and correcting fluxes using higher-degree polynomials have led to aneven more efficient high-order method — the flux reconstruction (FR) method [20, 21], also known asthe correction procedure via reconstruction (CPR) method [22]. Besides its better efficiency, by choosingdifferent correction polynomials, the FR method can recover many existing high-order schemes such as DGand SD, and can even produce new schemes that were never reported before. The stability of the FR methodhas been proved in [23]. New variants of the FR method have been reported in, e.g., [24–26]. The mostrecent developments on the FR method are summarized in [27].With more and more applications of high-order methods to flow simulations than ever before, there is anatural need of extending the sliding mesh concept to high-order methods to tackle complex flow problemssuch as those mentioned above. Ferrer and Willden [28] developed a high-order sliding-mesh DG methodbased on modal basis functions for simulating incompressible flow problems. Ram´ırez et al. [29] applied ∗ Corresponding author.
Email address: [email protected] (Bin Zhang)
Preprint submitted to Elsevier March 31, 2020 a r X i v : . [ m a t h . NA ] S e p oving-least-squares stencils to the development of a high-order sliding-mesh finite volume method. Theauthors of the present paper developed a simple and efficient high-order sliding-mesh SD method [30] usingiso-parametric mortar elements [31–33], and also extended this method to sliding-deforming meshes [34] andto handle 3D geometries [35]. Our previous methods, however, require uniform mesh on a sliding interface,which restricts mesh generation. Recent efforts have completely lifted this restriction, and the resultinggeneral nonuniform sliding-mesh method has been demonstrated on the FR method [36] and the SD methodon hybrid grids [37].Our initial effort in [36] showed that transfinite mortar elements can potentially make the sliding-meshmethod arbitrarily high-order accurate. In this work, we further the investigation to provide a more detailedstudy emphasizing the following aspects of this method: rotational speed effects on both spatial and temporalaccuracies, conservation property, outflow property, free-stream preservation property, capability of dealingwith multiple objects at the same time. Meanwhile, detailed steps on the implementation are also providedin this work.The rest of this paper is organized as follows. In Section 2, we briefly describe the equations that we aregoing to solve numerically. Sections 3 and 4 are about the numerical methods, including the FR methodand the new sliding-mesh method. Verifications and applications are reported in Section 5. Finally, Section6 concludes this paper.
2. The Flow Equations
We numerically solve the two-dimensional Navier-Stokes equations in the following conservative form, ∂ Q ∂t + ∂ F ∂x + ∂ G ∂y = , (1)where Q is the vector of conservative variables; F and G are the flux vectors in the x - and the y -direction,respectively. Their expressions are Q = [ ρ ρu ρv E ] T , (2) F = F inv ( Q ) + F vis ( Q , ∇ Q ) , (3) G = G inv ( Q ) + G vis ( Q , ∇ Q ) , (4)where ρ is fluid density, u and v are the Cartesian velocity components, and E is the total energy per volumedefined as E = pγ − ρ ( u + v ) , (5)where p is pressure, and γ is the ratio of specific heats and is set to 1.4 in this work for ideal gas. Theinviscid and the viscous flux vectors are F inv = ρuρu + pρuv ( E + p ) u , G inv = ρvρuvρv + p ( E + p ) v , (6) F vis = − τ xx τ yx uτ xx + vτ yx + κT x , G vis = − τ xy τ yy uτ xy + vτ yy + κT y , (7)where τ ij = µ ( u i,j + u j,i ) + λδ ij u k,k , is the shear stress tensor, µ is the dynamic viscosity, λ = − / µ based on the Stokes’ hypothesis, δ ij isthe Kronecker delta, κ is the thermal conductivity, and T is temperature which is related to density andpressure through the ideal gas law p = ρ R T, (8)2here R is the gas constant.To deal with moving grid, we employ an arbitrary-Lagrangian-Eulerian (ALE) approach, and map amoving physical domain to a fixed computational domain. Let ( t, x, y ) denote the physical time and space,and ( τ, ξ, η ) the computational ones. Further assume that the mapping is: t = τ , x = x ( τ, ξ, η ), and y = y ( τ, ξ, η ). It can be shown that the flow equations will take the following conservative form in thecomputational space ∂ (cid:101) Q ∂t + ∂ (cid:101) F ∂ξ + ∂ (cid:101) G ∂η = , (9)where (cid:101) Q , (cid:101) F , and (cid:101) G are the computational solution vector and flux vectors, and they are related to thephysical ones as (cid:101) Q = |J | Q , (10) (cid:101) F = ( − x t y η + y t x η ) Q + y η F − x η G , (11) (cid:101) G = ( x t y ξ − y t x ξ ) Q − y ξ F + x ξ G . (12)Alternatively, let Q , F , (cid:101) F , etc., each denote a component (at the same position) of the correspondingboldface vector, and then the relations in (10)-(12) can be written in the following matrix form (cid:101) Q (cid:101) F (cid:101) G = |J |J − QFG , (13)where J represents the Jacobian matrix, |J | is its determinant, and J − is the inverse Jacobian matrix.These metric terms have the following expressions, J = ∂ ( t, x, y ) ∂ ( τ, ξ, η ) = x τ x ξ x η y τ y ξ y η , |J | = x ξ y η − x η y ξ , (14) J − = ∂ ( τ, ξ, η ) ∂ ( t, x, y ) = ξ t ξ x ξ y η t η x η y = 1 |J | |J | − x t y η + y t x η y η − x η x t y ξ − y t x ξ − y ξ x ξ . (15)
3. Flux Reconstruction Method on Conforming Moving Grid
The first step of the FR method is to map each physical grid element to a standard computational element.In the present implementation, we discretize a computational domain into non-overlapping quadrilateralgrids, and employ the following iso-parametric mapping [38] to map each grid element to a unit squareelement (i.e., 0 ≤ ξ, η ≤
1) in the computational space, x ( t, ξ, η ) = (cid:34) x ( t, ξ, η ) y ( t, ξ, η ) (cid:35) = K (cid:88) i =1 M i ( ξ, η ) (cid:34) x i ( t ) y i ( t ) (cid:35) , (16)where K is the total number of nodes used to approximate a physical element, M i and ( x i , y i ) are the shapefunction and the coordinates of the i -th node. Fig. 1 is a schematic of the iso-parametric representations ofa physical element using different number of nodes. Higher-order elements (larger K ’s) obviously representcurved boundaries more accurately. For this reason, high-order elements should be used along curvedboundaries for better accuracy. 3 y (a) xy (b) xy (c) ξη (d) Figure 1: Iso-parametric representation (black lines) of a curved element (gray lines): (a) K = 4 (linear), (b) K = 8 (quadratic),(c) K = 12 (cubic), (d) a computational element (unit square). Within each computational element, solution points (SPs) and flux points (FPs) are defined. The SPsare distributed along each coordinate direction inside the element, and the FPs are distributed along theboundaries only. Figure 2 is a schematic of the distribution of these points for a third-order FR scheme. Inthis work, the N SPs/FPs in each direction of an N -th order scheme are chosen as the roots of the N -thLegendre polynomial (namely, N Legendre points). At the SPs, Lagrange interpolation bases are defined.For example, the basis at the i -th SP along the ξ direction is h i ( ξ ) = N (cid:89) s =1 ,s (cid:54) = i (cid:18) ξ − X s X i − X s (cid:19) , (17)where X i and X s are the ξ -coordinates of the i -th and the s -th SP, respectively. If we denote the space ofall polynomials of degrees less than or equal to N as P N , then h i ∈ P N − . Moreover, the h i ’s are linearlyindependent and form a basis for P N − .0 ξ η Figure 2: Schematic of the distribution of SPs (circular dots) and FPs (square dots) for a third-order ( P = 2) FR scheme. For a function φ defined within a mesh element, assume it is smooth (e.g., φ ∈ C ∞ ) and has discretevalues: φ ij at ( x i , y j ) that corresponds to ( X i , X j ) in the computational space, where i, j = 1 , , · · · , N . Thisfunction can then be approximated by the following tensor-product polynomial in P N − ,N − = P N − ⊗ P N − , φ ( ξ, η ) = N (cid:88) j =1 N (cid:88) i =1 φ ij h i ( ξ ) h j ( η ) , (18)where the truncation error for the approximation is O ( ξ N , η N ) based on Taylor series expansion.4pplying (18) to the solution and fluxes, we can obtain the following polynomial representations, (cid:101) Q ( ξ, η ) = N (cid:88) j =1 N (cid:88) i =1 (cid:101) Q ij h i ( ξ ) h j ( η ) , (19) (cid:101) F ( ξ, η ) = N (cid:88) j =1 N (cid:88) i =1 (cid:101) F ij h i ( ξ ) h j ( η ) , (20) (cid:101) G ( ξ, η ) = N (cid:88) j =1 N (cid:88) i =1 (cid:101) G ij h i ( ξ ) h j ( η ) , (21)where the () ij ’s are the discrete values at the ( i, j )-th SP in a standard element.The polynomials in (19)-(21) are continuous within each element, but discontinuous across elementinterfaces (or boundaries). For this reason, common values need to be defined at the interfaces. Forinstance, the common solution is computed as Q com = 12 ( Q L + Q R ) , (22)where Q L and Q R are the solution vectors on the left side and the right side of an interface.To compute the common (normal) inviscid flux, we employ a Riemann solver, such as the followingRusanov solver [39] with modification for moving grid, F cominv = 12 (cid:2) ( ↔ F Linv + ↔ F Rinv ) · n − λ ( Q R − Q L ) (cid:3) − ( v g · n ) Q com , (23)where ↔ F Linv = ( F Linv , G Linv ) and ↔ F Rinv = ( F Rinv , G Rinv ) are the inviscid flux vectors on the two sides of aninterface; n = N / (cid:107) N (cid:107) is the unit normal vector with N = ( y η , − x η ) or ( − y ξ , x ξ ) (24)depending on which direction (i.e., ξ or η ) the interface is mapped to; λ is the largest characteristic speedwith the following expression, λ = | V n | + c = (cid:12)(cid:12)(cid:12)(cid:0)
12 ( v L + v R ) − v g (cid:1) · n (cid:12)(cid:12)(cid:12) + c, (25)where v L = ( u L , v L ) and v R = ( u R , v R ) are flow velocities on the two sides of an interface; v g = ( x t , y t )represents grid velocity; c is the local speed of sound. The physical normal flux in (23) is converted to acomputational one by multiplying the magnitude of normal, i.e., (cid:101) F cominv or (cid:101) G cominv = (cid:107) N (cid:107) F cominv . (26)The common viscous flux is calculated from the common solution and common gradient, F comvis = F vis ( Q com , ( ∇ Q ) com ) , (27)where the common gradient is the average of the left and the right values, i.e.,( ∇ Q ) com = (( ∇ Q ) L + ( ∇ Q ) R ) / . (28)The procedure for calculating ∇ Q will be briefly discussed in the next section. The common viscous flux in(27) is converted to a computational one in the same way as (26).5 .3. Flux Reconstruction The spatial derivatives in (9) reduces the two flux terms to P N − ,N − and P N − ,N − , making theminconsistent with the solution term which is P N − ,N − . To overcome this issue, the flux polynomials needto be reconstructed to be at least P N,N − and P N − ,N .To do this, we use higher degree correction functions/polynomials. The corrected/reconstructed fluxestake the following form, (cid:98) F ( ξ, η ) = (cid:101) F ( ξ, η ) + (cid:2)(cid:101) F com (0 , η ) − (cid:101) F (0 , η ) (cid:3) · g L ( ξ ) + [ (cid:101) F com (1 , η ) − (cid:101) F (1 , η )] · g R ( ξ ) , (29) (cid:98) G ( ξ, η ) = (cid:101) G ( ξ, η ) + [ (cid:101) G com ( ξ, − (cid:101) G ( ξ, · g L ( η ) + [ (cid:101) G com ( ξ, − (cid:101) G ( ξ, · g R ( η ) , (30)where (cid:101) F ( ξ, η ) and (cid:101) G ( ξ, η ) are the original flux polynomials from (20) and (21); g L and g R are the left andthe right correction functions with degrees no less than N , and they are required to at least satisfy g L (0) = 1 , g L (1) = 0 ,g R (0) = 0 , g R (1) = 1 , (31)which ensures that (cid:98) F (0 , η ) = (cid:101) F com (0 , η ) , (cid:98) F (1 , η ) = (cid:101) F com (1 , η ) , (32) (cid:98) G ( ξ,
0) = (cid:101) G com ( ξ, , (cid:98) G ( ξ,
1) = (cid:101) G com ( ξ, , (33)i.e., the reconstructed fluxes still take the common values at cell interfaces. Huynh [20] proposed severalcorrection functions, and we employ the g DG correction function in this study.Note that to calculate solution gradients consistently, the correction procedure is also applied to thesolution polynomial along each direction (only for calculating the gradients). In this way, the resultinggradient polynomials are in P N − ,N − as well. With proper boundary conditions applied, the discretized equations can be written in the followingresidual form, ∂ (cid:101) Q ∂t (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ij = − (cid:34) ∂ (cid:98) F ∂ξ + ∂ (cid:98) G ∂η (cid:35) ij = R ij , i, j = 1 , , · · · , N, (34)where R ij is the residual at the ( i, j )-th SP. This system can be time marched using either explicit orimplicit schemes. In this work, we employ several of the explicit strong stability preserving (SSP) Runge-Kutta schemes reported in [40–42] for the purpose. Furthermore, in this work, all boundary conditions areweakly imposed to increase stability, and more detail can be found in, e.g., [15]. Ideally, a moving grid should not disturb a flow field. The simplest situation is that a free-stream flowmust stay constant all the time on a moving grid. This is called free-stream preservation. By substituting aconstant flow solution into the flow equations in (9), we can get a system of equations that are purely aboutthe geometrics, and are thus known as the geometric conservation law (GCL) [43], ∂ ( |J | ξ x ) ∂ξ + ∂ ( |J | η x ) ∂η = 0 , (35) ∂ ( |J | ξ y ) ∂ξ + ∂ ( |J | η y ) ∂η = 0 , (36) ∂ |J | ∂t + ∂ ( |J | ξ t ) ∂ξ + ∂ ( |J | η t ) ∂η = 0 . (37)To numerically satisfy free-stream preservation, the same numerical schemes for discretizing the flow equa-tions must be applied to the GCL equations. Since the spatial discretization operator in the FR method6s direct differentiation (which is exact), and the geometric terms are from analytical mapping, the firsttwo GCL equations are hence automatically satisfied. But the third GCL equation generally can not besatisfied automatically. This problem come from the temporal discretization, for instance, a multiple-stageRunge-Kutta scheme, which is not exact. To overcome this issue, we treat |J | as an unknown, and solvethe third equation numerically. The numerical |J | is then used to update the physical solution accordingto (10). In this way, the GCL is numerically satisfied, and free-stream preservation is ensured. Similarapproach was reported in, e.g., [44].
4. A Sliding-Mesh Method using Transfinite Mortar Elements
In this section, we first introduce the concepts of sliding mesh and mortar element. Following that, a briefdescription is given on transfinite mapping. The projection procedures between cell faces and mortars aredescribed subsequently. After that, we prove that the present method is globally conservative and satisfiesoutflow condition. We also provide details on the implementation at the end of this section.
The simplest sliding mesh involves only two non-overlapping subdomains. Such an example is shownin Fig. 3(a), where the two subdomains are allowed to have relative rotational motions, which leads to anonconforming sliding interface in between. To ensure continuity of solution and conservation of fluxes onthis nonconforming mesh, we employ mortar elements as the communicator between the two subdomains.The distribution of mortar elements for this mesh is shown in Fig. 3(b). A mortar element is formedbetween two successive points on the sliding interface. At every time instant, a mortar element is connectedto a cell face on its left and a cell face on its right (from a counterclockwise perspective). A cell face isconnected to one or multiple mortars on its one side. We call this information as the cell face and mortarconnectivities. Note that these connectivities are time-dependent, and need to be updated at every stageof a time marching scheme. An efficient procedure for updating these connectivities will be presented in alater section.
Figure 3: Left: a nonconforming sliding mesh; right: distribution of mortar elements (represented by hatched lines).
Each sliding cell face is mapped to a unit straight line element (e.g., 0 ≤ ξ ≤
1) when the underlyinggrid cell is mapped to a unit square element in the computational space. Similarly, we also map each mortarelement to a unit straight line element 0 ≤ z ≤ Ξ Ξ k Ξ m ...... Ω Ξ Ξ k Ξ m ...... Ξ k Figure 4: Mapping a cell face and its mortars to unit line elements: left, physical space; middle, computational space; right,mortar space. x x x f = x Ω x f x x f x x f xy (0 , η = 0 (1 , ξ = (1 , η = 1 (0 , ξ = ξη Figure 5: Transfinite mapping of a physical element to a unit square element. task. The advantage of transfinite mapping is that it represents geometries that have analytical expressions(such as a circular arc) exactly, i.e., introduces no geometric error at all.Figure 5 is a schematic of a physical element and the corresponding computational element. The trans-finite mapping between these two elements can be expressed as x ( t, ξ, η ) = (1 − η ) x f ( t, ξ ) + ξ x f ( t, η ) + η x f ( t, ξ ) + (1 − ξ ) x f ( t, η ) − (1 − ξ )(1 − η ) x ( t ) − ξ (1 − η ) x ( t ) − ξη x ( t ) − (1 − ξ ) η x ( t ) , (38)where the x i ’s denote coordinates of the corner nodes, and x f i ’s are expressions of the faces, and these termsare all time-dependent for moving grids. If the faces are represented by one-dimensional iso-parametricmappings of the same order, then the transfinite mapping is equivalent to the iso-parametric mapping in(16). If a face has an exact expression, then that face is represented exactly by the transfinite mapping.In our case, assume face f is a circular arc (i.e., corresponds to the sliding cell face Ω in Fig. 4), thenit can be analytically expressed as x Ω ( t, ξ ) = (cid:34) x Ω ( t, ξ ) y Ω ( t, ξ ) (cid:35) = (cid:34) R · cos (cid:2) (1 − ξ ) θ Ω1 ( t ) + ξθ Ω2 ( t ) (cid:3) + x c ( t ) R · sin (cid:2) (1 − ξ ) θ Ω1 ( t ) + ξθ Ω2 ( t ) (cid:3) + y c ( t ) (cid:35) , (39)where R and ( x c , y c ) are the radius and the center coordinates of the arc; θ Ω1 and θ Ω2 are the angles corre-sponding to the starting point (i.e., x ) and the ending point (i.e., x ), respectively, of the face. Similarly,a mortar element Ξ k has the following exact expression, x Ξ k ( t, z ) = (cid:34) x Ξ k ( t, z ) y Ξ k ( t, z ) (cid:35) = (cid:34) R · cos (cid:2) (1 − z ) θ Ξ k ( t ) + zθ Ξ k ( t ) (cid:3) + x c ( t ) R · sin (cid:2) (1 − z ) θ Ξ k ( t ) + zθ Ξ k ( t ) (cid:3) + y c ( t ) (cid:35) , (40)where θ Ξ k and θ Ξ k are the starting and the ending angles of the mortar.The analytical expressions in (39) and (40) are actually equivalent to the following linear mappings ofthe angles, θ Ω = (1 − ξ ) θ Ω1 ( t ) + ξθ Ω2 ( t ) , (41)8 Ξ k = (1 − z ) θ Ξ k ( t ) + zθ Ξ k ( t ) , (42)where θ Ω is the angle of a physical point (on Ω) that is mapped to a point ξ in the computational space,and θ Ξ k is the angle of a physical point (on Ξ k ) that is mapped to a point z in the mortar space. Referringto the physical space in Fig. 4 and considering the fact that θ Ω = θ Ξ k represents the same physical point,the following relation thus holds between the computational space and the mortar space, ξ = o k + s k · z, (43)where s k and o k are the scaling and the offset of the mortar Ξ k with respect to the cell face Ω. Morespecifically, s k = θ Ξ k − θ Ξ k θ Ω2 − θ Ω1 = ∆ θ Ξ k ∆ θ Ω = R · ∆ θ Ξ k R · ∆ θ Ω = L Ξ k L Ω , (44) o k = ∆ θ Ξ + ∆ θ Ξ + · · · + ∆ θ Ξ k − ∆ θ Ω = (cid:88) k − α =1 s α , (45)where L Ξ k and L Ω are the physical lengths of the mortar and the face. Note that the scaling and the offsetare both time-dependent, and are updated at the time the connectivities are updated.Two comparison examples of the iso-parametric mapping and the transfinite mapping are included inAppendix A. Communications on a sliding interface include: projection of local discontinuous values from cell facesto mortars, computation of common values on mortars, and projection of common values back to cell faces.These procedures are discussed in details in what follows. To facilitate the discussion, we adopt the followingnotations: Q denotes a component of Q ; Q i denotes the discrete value of Q at the i -th FP; Q denotes thevector ( Q , Q , · · · , Q N ). This same rule also applies to fluxes. Take solution as an example. Each solution component on a cell face is represented by the followingone-dimensional polynomial Q Ω ( ξ ) = N (cid:88) i =1 Q Ω i h i ( ξ ) . (46)If we define the same set of FPs in the mortar space, then the solution polynomial on a mortar can beconstructed in the same way. For example, on the left side of Ξ k , the solution polynomial is Q Ξ k , L ( z ) = N (cid:88) i =1 Q Ξ k , L i h i ( z ) , (47)where Q Ξ k , L i is a solution component at the i -th FP on the left side of Ξ k . As illustrated in Fig. 6, to getΩ Ξ Ξ k Ξ m ...... Figure 6: Projection from cell faces to the two sides of a mortar. k , we require (cid:90) (cid:0) Q Ξ k , L ( z ) − Q Ω ( ξ ) (cid:1) h j ( z ) d z = 0 , ∀ j = 1 , , ..., N. (48)Substituting (46) and (47) into the above equation and considering the relation in (43), we will get thefollowing equation system M Q Ξ k , L = S Ω → Ξ k Q Ω , (49)where the elements of the coefficient matrices are M ij = (cid:90) h i ( z ) h j ( z ) d z, i, j = 1 , , ..., N, (50) S Ω → Ξ k ij = (cid:90) h i ( o k + s k z ) h j ( z ) d z, i, j = 1 , , ..., N. (51)Solutions of (49), when written in matrix form, are Q Ξ k , L = P Ω → Ξ k Q Ω = M − S Ω → Ξ k Q Ω , (52)where P Ω → Ξ k is the projection matrix from Ω to Ξ k . This process is repeated for all the solution components.In the same way, we can get solutions on the right side of a mortar. Also note that the integrals in (50)and (51) can be evaluated exactly and efficiently using quadratures, e.g., the Clenshaw-Curtis quadratures.The Legendre points being the SPs/FPs in this work makes the h i ’s orthogonal, and in turn makes the M matrix diagonal and trivial to invert. The inversion actually needs to be done only once during initializationbecause the M matrix is also time-independent. The common solution Q Ξ k and the inviscid normal flux F Ξ k inv on a mortar are computed in the same wayas on a cell interface, i.e., Q Ξ k = 12 ( Q Ξ k , L + Q Ξ k , R ) , (53) F Ξ k inv = 12 (cid:2) ( ↔ F Ξ k , Linv + ↔ F Ξ k , Rinv ) · n − λ ( Q Ξ k , R − Q Ξ k , L ) (cid:3) − ( v g · n ) Q Ξ k , (54)where the variables, without further explanation, have similar meanings to those in (23). Figure 7 illustrates the process of projecting common values from m mortars back to a cell face Ω.Taking flux as an example, we can either directly project back the physical flux (method 1) or convert it tocomputational flux to project back (method 2). The details are described in what follows.Ω Ξ Ξ k Ξ m ...... Figure 7: Project common values from mortars back to cell faces. ethod 1: The inviscid normal fluxes on a cell face and a mortar are represented by the following polynomials, F Ωinv ( ξ ) = N (cid:88) i =1 F Ωinv ,i h i ( ξ ) , F Ξ k inv ( z ) = N (cid:88) i =1 F Ξ k inv ,i h i ( z ) . (55)To get the F Ωinv ,i ’s, we require m (cid:88) k =1 (cid:90) o k + s k o k (cid:16) F Ωinv ( ξ ) − F Ξ k inv ( z ) (cid:17) h j ( ξ ) d ξ = 0 , ∀ j = 1 , , ..., N, (56)which gives M F Ωinv = m (cid:88) k =1 s k S Ξ k → Ω F Ξ k inv , (57)where M is identical to that in (49), and S Ξ k → Ω is simply the transpose of the S Ω → Ξ k matrix from (49).Solutions of (57) are F Ωinv = m (cid:88) k =1 P Ξ k → Ω F Ξ k inv = m (cid:88) k =1 s k M − S Ξ k → Ω F Ξ k inv , (58)where P Ξ k → Ω is the projection matrix from Ξ k to Ω. The above physical flux on a cell face is then convertedto a computational one to compute residuals. The conversion follows (26), here the normal is N = ( y Ω ξ , − x Ω ξ ) = L Ω (cos θ Ω , sin θ Ω ) , (59)and the final computational flux on the cell face is (cid:101) F Ωinv = (cid:107) N (cid:107) F Ωinv = L Ω F Ωinv . (60) Method 2:
Alternatively, we can convert the the physical flux in (54) to computational, and then do the backprojection. For a mortar, there are two types of normals depending on which space, i.e., the mortar spaceor the computational space, is taken as the reference space. For the mortar space, the normal is˘ N = ( y Ξ k z , − x Ξ k z ) = L Ξ k (cos θ Ξ k , sin θ Ξ k ) . (61)For the computational space, the normal is (cid:101) N = ( y Ξ k ξ , − x Ξ k ξ ) = ( y Ξ k z z ξ , − x Ξ k z z ξ ) = 1 s k ( y Ξ k z , − x Ξ k z ) = 1 s k ˘ N . (62)The corresponding fluxes in these two spaces are˘ F Ξ k inv = (cid:107) ˘ N (cid:107) F Ξ k inv = L Ξ k F Ξ k inv , (63) (cid:101) F Ξ k inv = (cid:107) (cid:101) N (cid:107) F Ξ k inv = 1 s k ˘ F Ξ k inv . (64)The computational inviscid fluxes on a cell face and a mortar are represented by the following polynomials, (cid:101) F Ωinv ( ξ ) = N (cid:88) i =1 (cid:101) F Ωinv ,i h i ( ξ ) , (cid:101) F Ξ k inv ( z ) = N (cid:88) i =1 (cid:101) F Ξ k inv ,i h i ( z ) . (65)To get the (cid:101) F Ωinv ,i ’s, we require m (cid:88) k =1 (cid:90) o k + s k o k (cid:16) (cid:101) F Ωinv ( ξ ) − (cid:101) F Ξ k inv ( z ) (cid:17) h j ( ξ ) d ξ = 0 , ∀ j = 1 , , ..., N, (66)11hich leads to M (cid:101) F Ωinv = m (cid:88) k =1 s k S Ξ k → Ω (cid:101) F Ξ k inv = m (cid:88) k =1 S Ξ k → Ω ˘ F Ξ k inv , (67)where M and S Ξ k → Ω are identical to those in (57). The solutions of the above system are (cid:101) F Ωinv = m (cid:88) k =1 P Ξ k → Ω ˘ F Ξ k inv = m (cid:88) k =1 M − S Ξ k → Ω ˘ F Ξ k inv , (68)and note the difference between the P Ξ k → Ω matrices in the above equation and in Eq. (58).These two methods are in fact equivalent as one may expect. To show this, simply divide both sides of(68) by L Ω , and consider (44), (60), and (63), and we will get exactly the same result as (58). Also notethe singularity in (64) when the scaling (size) of a mortar becomes zero. This singularity is eliminated in(67) where the scaling is multiplied back. Equation (64) is for derivation purpose only, and is not actuallyevaluated in the computation, so the singularity is naturally avoided. For viscous flow, we first project the common solution (53) back to cell faces in the same way as (58).The updated solutions on a cell face are then involved in the computation of solution gradients. After that,we could either project local gradients (method 1) or local viscous fluxes (method 2) to mortars to computecommon viscous fluxes which are then projected back. The steps are described below.
Method 1:
The local gradients on cell faces are projected to mortars following (52). The common gradients andcommon physical viscous fluxes on a mortar are then computed following (28) and (27), respectively. Thenormal viscous flux is calculated as ˘ F Ξ k vis = ↔ F Ξ k vis · ˘ N , (69)where ↔ F Ξ k vis = ( F Ξ k vis , G Ξ k vis ) with the two components representing the physical common viscous fluxes. Thisnormal viscous flux is finally projected back to cell faces following (68). Method 2:
Local viscous flux, denoted by (cid:101) F Ωvis , is projected to mortars in the same way as (48). The resultingnormal viscous flux on the left side of a mortar is˘ F Ξ k , Lvis = P Ω → Ξ k (cid:101) F Ωvis = s k M − S Ω → Ξ k (cid:101) F Ωvis . (70)Note the difference on the P Ω → Ξ k ’s in the above equation and in (52). The viscous flux on the right side ofa mortar is obtained in the same way. The common normal viscous flux is then calculated as˘ F Ξ k vis = 12 ( ˘ F Ξ k , Lvis + ˘ F Ξ k , Rvis ) , (71)which is projected back to cell faces following (68) to update the original normal viscous flux.We have compared method 1 and method 2 in a series of tests (not reported here), and do not noticeany obvious difference on the results. But method 2 is slightly faster than method 1, because it requiresfewer projections and calculations. Global conservation means that, without the presence of source in a domain, the net gain or loss of (cid:101) Q over the whole domain is determined solely by fluxes through the boundaries. Inspired by the works in[15, 32], we prove in what follows that the present method is globally conservative.Define the following quadrature weights w i and w j at the SPs such that (cid:90) (cid:90) φ ( ξ, η )d ξ d η = N (cid:88) j =1 N (cid:88) i =1 φ ij w i w j , ∀ φ ∈ P N − ,N − . (72)12ultiply the discretized equation (34) by the weights and sum over all the SPs N (cid:88) j =1 N (cid:88) i =1 ∂ (cid:101) Q ∂t (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ij w i w j = − N (cid:88) j =1 N (cid:88) i =1 (cid:34) ∂ (cid:98) F ∂ξ + ∂ (cid:98) G ∂η (cid:35) ij w i w j . (73)Since the quadrature is exact for P N − ,N − by its definition, and the terms ∂ (cid:101) Q /∂t , ∂ (cid:98) F /∂ξ , ∂ (cid:98) G /∂η are in P N − ,N − , (73) is thus equivalent to (cid:90) (cid:90) ∂ (cid:101) Q ∂t d ξ d η = (cid:90) (cid:90) (cid:34) ∂ (cid:98) F ∂ξ + ∂ (cid:98) G ∂η (cid:35) d ξ d η. (74)Take the time derivative out, and apply the divergence theorem to the right hand side, ∂∂t (cid:90) (cid:90) (cid:101) Q d ξ d η = (cid:90) (cid:98) F (0 , η )d η − (cid:90) (cid:98) F (1 , η )d η + (cid:90) (cid:98) G ( ξ, ξ − (cid:90) (cid:98) G ( ξ, ξ, (75)where, by definition, the corrected fluxes (cid:98) F and (cid:98) G take the common values on cell interfaces. For aconforming mesh, if we sum (75) over all the mesh elements, the interior flux terms cancel on cell interfaces,leaving only the boundary flux terms. This confirms that the change of (cid:101) Q over the whole domain dependsonly on boundary fluxes, i.e., the scheme is globally conservative on conforming mesh.For nonconforming sliding mesh, we need to show that the total flux on the left side of a sliding interfaceequals that on the right side so that global conservation is satisfied. Since we have required (66) to hold forboth inviscid and viscous fluxes with any h j , and the h j ’s form a basis of P N − , the following relation thusholds, m (cid:88) k =1 (cid:90) o k + s k o k (cid:16) (cid:101) F Ω ( ξ ) − (cid:101) F Ξ k ( z ) (cid:17) ψ ( ξ ) d ξ = 0 , ∀ ψ ∈ P N − , (76)Substitution of ψ = 1 ∈ P N − gives m (cid:88) k =1 (cid:90) o k + s k o k (cid:16) (cid:101) F Ω ( ξ ) − (cid:101) F Ξ k ( z ) (cid:17) d ξ = 0 . (77)Rearrange the above equation to m (cid:88) k =1 (cid:90) o k + s k o k (cid:101) F Ω ( ξ )d ξ = m (cid:88) k =1 (cid:90) o k + s k o k (cid:101) F Ξ k ( z )d ξ, then change the variable of integration on the right hand side, (cid:90) (cid:101) F Ω ( ξ )d ξ = m (cid:88) k =1 (cid:90) s k (cid:101) F Ξ k ( z )d z, and now consider the relation in (64) (which also applies to viscous flux), (cid:90) (cid:101) F Ω ( ξ )d ξ = m (cid:88) k =1 (cid:90) ˘ F Ξ k ( z )d z. (78)The relation in (78) is valid for any sliding cell face Ω, and it says that the flux on a sliding cell face equalsthe total flux on all its mortars. Since each mortar is unique and does not overlap with other mortars (seeFig. 3), summing (78) over all the cell faces on the left side of a sliding interface will lead us to the conclusionthat the total flux on the left side of a sliding interface equals the total flux on all the mortars. The sameconclusion also holds on the right side of a sliding interface. Therefore, the total flux on the left side of asliding interface equals that on the right side, i.e., the scheme satisfies global conservation.13 .5. Proof of Outflow Condition In a hyperbolic system, the propagation of waves should not be affected by downwind signals, whichis called the outflow condition. When it comes to our case, satisfying the outflow condition means that ifwe project a variable from a cell face to its mortars, and then immediately project back, the value of thevariable should not change. In what follows, we prove that the present sliding-mesh method satisfies thiscondition. Since the two projection methods in Sec. 4.3.3 are equivalent, we only show the poof on the firstmethod for simplicity.Let φ Ω denote a variable on a cell face, φ Ξ k the projected variable on the k -th mortar, and φ ∗ Ω thevariable from back projection. Then, the problem becomes proving that φ ∗ Ω = φ Ω .Start with the following polynomial representations of these variables, φ Ω ( ξ ) = N (cid:88) i =1 φ Ω i h i ( ξ ) , φ Ξ k ( z ) = N (cid:88) i =1 φ Ξ k i h i ( z ) , φ ∗ Ω ( ξ ) = N (cid:88) i =1 φ ∗ Ω i h i ( ξ ) . (79)The projection from cell face to mortar follows (48), (cid:90) (cid:0) φ Ξ k ( z ) − φ Ω ( ξ ) (cid:1) h j ( z ) d z = 0 , ∀ j = 1 , , ..., N, (80)from which we can get the projected variable on the mortar, φ Ξ k = P Ω → Ξ k φ Ω = M − S Ω → Ξ k φ Ω . (81)Now, project φ Ξ k back to the cell face following (56), m (cid:88) k =1 (cid:90) o k + s k o k (cid:0) φ ∗ Ω ( ξ ) − φ Ξ k ( z ) (cid:1) h j ( ξ ) d ξ = 0 , ∀ j = 1 , , ..., N, (82)from which we can get the variable from back projection, φ ∗ Ω = m (cid:88) k =1 P Ξ k → Ω φ Ξ k = m (cid:88) k =1 s k M − S Ξ k → Ω φ Ξ k . (83)Since the the h j ’s form a basis for P N − , from (80) we can get (cid:90) (cid:0) φ Ξ k ( z ) − φ Ω ( ξ ) (cid:1) ψ d z = 0 , ∀ ψ ∈ P N − . (84)Because h j ( ξ ) = h j ( o k + s k z ) ∈ P N − , setting ψ = h j ( ξ ) gives (cid:90) (cid:0) φ Ξ k ( z ) − φ Ω ( ξ ) (cid:1) h j ( ξ ) d z = 0 , now rearrange the terms and change the variable of integration, (cid:90) o k + s k o k φ Ω ( ξ ) h j ( ξ ) d ξ = s k (cid:90) φ Ξ k ( z ) h j ( ξ ) d z. (85)Apply (85) to all the mortars of Ω and sum them up, m (cid:88) k =1 (cid:90) o k + s k o k φ Ω ( ξ ) h j ( ξ ) d ξ = m (cid:88) k =1 s k (cid:90) φ Ξ k ( z ) h j ( ξ ) d z, (cid:90) φ Ω ( ξ ) h j ( ξ ) d ξ = m (cid:88) k =1 s k (cid:90) φ Ξ k ( z ) h j ( ξ ) d z, (86)14n the other hand, rearranging (82) and changing the variable of integration gives (cid:90) φ ∗ Ω ( ξ ) h j ( ξ )d ξ = m (cid:88) k =1 s k (cid:90) φ Ξ k ( z ) h j ( ξ ) d z. (87)Note that (86) and (87) have exactly the same right hand side, therefore (cid:90) φ ∗ Ω ( ξ ) h j ( ξ )d ξ = (cid:90) φ Ω ( ξ ) h j ( ξ )d ξ, ∀ j = 1 , , ..., N, (88)which is equivalent to the following equation system, M φ ∗ Ω = M φ Ω . Because M is invertible, we therefore have φ ∗ Ω = M − M φ Ω , that is φ ∗ Ω = φ Ω , (89)which is exactly what we were expecting, and the present method therefore satisfies the outflow condition.If we substitute (81) into (83) and take (89) into account, we have the following relation φ ∗ Ω = m (cid:88) k =1 P Ξ k → Ω P Ω → Ξ k φ Ω = φ Ω , (90)which readily reveals a property of the projection matrices, m (cid:88) k =1 P Ξ k → Ω P Ω → Ξ k = I , (91)where I is the identity matrix.We have numerically tested (91) on many cases, and it holds perfectly without exception. As it was notedin [32] that satisfaction of the outflow condition ensures that a method does not change the characteristicsof the flow field, and therefore should not degrade the maximum allowable time step size compared withthe original method on conforming mesh. Our observation through numerical tests is consistent with thisexpectation. A flow chart of the implementation procedures is shown in Fig. 8. Some of the steps are explained withmore details in what follows.
Mesh for each subdomain is generated independently and stored in a separate file. When subdomainmeshes are read in, they are assembled into a “single” mesh by adding offsets to the numberings of cells andvertices. For example, as shown in Fig. 9, assume we have n subdomains with N , N , ..., N n cells. In theassembled mesh, numbering for the cells in subdomain 1 starts from 1, for subdomain 2 starts from N + 1,for subdomain 3 starts from N + N + 1, and so on. The same rule applies to vertices and boundary faces.In this way, cells and vertices are numbered uniquely and continuously. The overall mesh behaves just likea single mesh and sliding interfaces are like interior boundaries. All sliding cell faces are stored consecutively in an array. This means that cell faces from the left side ofthe first sliding interface are stored first, followed by those from the right side of the first sliding interface,and then those from the left side of the second sliding interface, and so on. Cell faces from each side arereordered into counterclockwise order, which makes connectivity updating more efficient. The algorithm forthe reordering is quite simple: any cell face can be assigned as the first cell face, and then the next cell faceis the one whose starting vertex is the ending vortex of the current cell face, repeat until finished.15ead meshes, and perform initialization • reorder sliding cell faces • correct sliding interfaces • compute the M matrix and its inverse • calculate lengths of sliding cell facesupdate grid motion • update connecticities for sliding • update L Ξ , o , s , and ˘ N for each mortar • update projection matrices • update J and J − for cells along sliding interfacescompute inviscid flux • project local solution to mortars • calculate common solution/inviscid flux on mortars • project common inviscid flux back to cell facesviscous flow?compute viscous flux • project common solution back to cell faces • calculate local viscous flux, and project to mortars • calculate common visous flux, and project to cell facesyescompute residuals and update solutionsconverged orat max iteration?stopyes nono Figure 8: Flow chart of the implementation procedures.
Mesh generator may introduce geometric errors. When it comes to sliding mesh, the problem is usuallythat, the vertices on a circular sliding interface do not represent the same radius. This geometric error maypotentially contaminate a simulation. Here is an easily fix that can be performed during preprocessing: pickup a vertex and take the corresponding radius as a reference radius, and then use this reference radius and16 ubdomain1 subdomain2 subdomain n − subdomain n · · · Figure 9: Schematic of n subdomains with ( n −
1) sliding interfaces (inner subdomains have been scaled). the original angle of each vertex to update the coordinates of that vertex. In this way, all the vertices on acircular sliding interface represent the same radius.
For simplicity, we only take the first sliding interface from Fig. 9 to explain how the connectivitiesare efficiently updated. The same procedure is repeated for other sliding interfaces (there is no limit onhow many we may have). Table 1 shows the definitions of important numbers and arrays used during thisprocess. The detailed procedures for updating the connectivities are listed as Algorithm 1. The general ideais that we “walk” along each sliding interface counterclockwise, and two successive points form a mortar(these two points might or might not come from the same side). One may also refer to Fig. 3 to facilitatethe understanding of the algorithm.Variable Meaning nfl number of cell faces on the left side nfr number of cell faces on the right side nf total number of faces on this sliding interfacenote: nf = nfl + nfr (by definition) nm total number of mortars on this sliding interfacenote: nm = nf (always valid) vof (1:2,1: nf ) vof (1:2, ifa ) stores the two vertices of the ifa -th cell face mof (1:2,1: nf ) mof (1 , ifa ) stores the first mortar of the ifa -th cell face mof (2 , ifa ) stores the number of mortars of the ifa -th cell face fom (1:2,1: nm ) fom (1 , im ) stores the left cell face of the im -th mortar fom (2 , im ) stores the right cell face of the im -th mortar vom (1:2,1: nm ) vom (1:2, im ) stores the two vertices of the im -th mortar Table 1: Definitions of variables used in updating cell face and mortar connectivities on a sliding interface.
The flow chart in Fig. 8 only shows the steps related to sliding. These steps must be injected into themain solver at the right locations. For example, the “compute inviscid flux” part is called after interior and17 lgorithm 1
Algorithm for updating cell face and mortar connectivities for a sliding interface mof =0; fom =0; vom =0 (cid:46) initialize with zeros ifl = 1 for ifr = ( nfl + 1) to nf doif vof (1 , ifl ) lies between vof (1 , ifr ) and vof (2 , ifr ) then (cid:46) the first mortar is locatedexit end ifend for im ← mof (1 , ifl ) ← im (cid:46) connectivities of the first mortar mof (2 , ifl ) ← mof (2 , ifl ) + 1 mof (2 , ifr ) ← mof (2 , ifr ) + 1 fom (1 , im ) ← ifl fom (2 , im ) ← ifr vom (1 , im ) ← vof (1 , ifl ) for im = 2 to nm do (cid:46) connectivities of remaining mortars if vof (2 , ifl ) lies between vof (1 , ifr ) and vof (2 , ifr ) then ifl ← ifl + 1 ifa ← ifl else ifr ← ifr + 1 if ( ifr > nf ) then ifr ← ifr − nfr end if ifa ← ifr end if mof (1 , ifa ) ← im mof (2 , ifl ) ← mof (2 , ifl ) + 1 mof (2 , ifr ) ← mof (2 , ifr ) + 1 fom (1 , im ) ← ifl fom (2 , im ) ← ifr vom (1 , im ) ← vof ( , ifa ) vom (2 , im − ← vof ( , ifa ) end for vom (2 , nm ) ← vom (1 , (cid:46) mortars always form a closed loopcell interface inviscid flux computations are done. The back projection of common solution takes place afterall inviscid flux computations are finished and before viscous flux computations start. The projections ofviscous flux are carried out after all interior and cell interface viscous flux computations are done.
5. Examples
In this section, we apply the sliding mesh method to several flow problems. We first test it on an inviscidflow and a viscous flow to verify the spatial and temporal accuracies. We then study the conservation and thefree-stream-preservation properties of the method. Following that, we perform a comparison study betweenthe present method and a rigid-rotation method on flow over a rotating square cylinder. Finally, we applythe method to simulate flow over multiple rotating square cylinders to further demonstrate its capability.18or all the test cases wherever applicable, we employ the following L norm to measure the errors, L error = (cid:115) (cid:80) N DOF i =1 ( φ i − φ exact i ) N DOF , (92)where φ represents the variable of interest, φ i and φ exact i are the numerical and the exact solutions at the i -th degree of freedom (DOF), N DOF = N elem · N is the total number of DOFs, N elem is the total numberof mesh elements, N = P + 1 is the scheme order (also the number of SPs in each direction of an element),and P is the polynomial degree. This is an inviscid flow test. In this flow, an isentropic vortex is superimposed to and convected by auniform flow. The flow field in an infinite domain at a time instant t can be analytically expressed as u = U ∞ (cid:26) cos θ − (cid:15)y r r c exp (cid:18) − x r − y r r c (cid:19)(cid:27) , (93) v = U ∞ (cid:26) sin θ + (cid:15)x r r c exp (cid:18) − x r − y r r c (cid:19)(cid:27) , (94) ρ = ρ ∞ (cid:26) − ( γ − (cid:15)M ∞ ) (cid:18) − x r − y r r c (cid:19)(cid:27) γ − , (95) p = p ∞ (cid:26) − ( γ − (cid:15)M ∞ ) (cid:18) − x r − y r r c (cid:19)(cid:27) γγ − , (96)where U ∞ , ρ ∞ , p ∞ , M ∞ are the speed, density, pressure and Mach number of the uniform flow; θ denotesthe mean flow direction; (cid:15) and r c are the strength and size of the vortex; ( x r , y r ) = ( x − x − ¯ ut, y − y − ¯ vt )are the relative coordinates; ( x , y ) represent the initial position of the vortex; (¯ u, ¯ v ) = ( U ∞ cos θ, U ∞ sin θ )are the velocity components of the mean flow.For this test, we have chosen the following parameters: U ∞ = 1, ρ ∞ = 1, M ∞ = 0 . θ = arctan(1 / (cid:15) = 1, and r c = 1. The overall computational domain has a size of 0 ≤ x, y ≤
10, and the vortex is initiallyplaced at ( x , y ) = (5 , ω = 0 , , ,
10, 15 and 20, are tested to study their effects on the solution. For all thecases, a five-stage fourth-order SSP Runge-Kutta scheme [40, 41] with a time step size of ∆ t = 1 . × − is used for the time marching. Figure 10: Mesh for Euler vortex flow simulation.
19n Fig. 11, we compare the density contours at t = 2 from the ω = 20 case using different polynomials.At this time instant, the vortex center travels right onto the sliding interface. It is obvious that P = 2 doesnot provide enough resolution as the vortex is poorly resolved. But as the polynomial degree increases, thesolution quality becomes much improved. Even at a small polynomial degree of P = 4 and on such a coarsegrid, the details of the vortex are very well captured. (a) P = 2 (b) P = 3 (c) P = 4 ρ Figure 11: Density contours of Euler vortex flow at t = 2 (dashed lines represent sliding interface). The L errors of density are plotted in Fig. 12 against polynomial degrees. It is seen that the errorsfrom different cases are comparable when P ≤
13, and start differing when
P >
13. The reason for this isthat spatial errors dominate in the first regime, and they keep decreasing to such a small level that temporalerrors start to dominate in the second regime. And larger rotational speed induces larger temporal errors asthe second regime of Fig. 12 shows. Similar observation on the effects of rotational speed was reported in[28]. But in that work, the effects start showing up at a very early stage of P = 4, which is mainly due to thelow-order (2nd order) temporal scheme used in that work. The present results suggest that the rotationalspeed effects can be dramatically reduced by using a high-order temporal scheme (4th order here).1 3 5 7 9 11 13 15 1710 − − − − − − P L e rr o r ω = 0 ω = 1 ω = 5 ω = 10 ω = 15 ω = 20 Figure 12: L error of density against polynomial degree for Euler-vortex flow simulation. A closer look at the curves in Fig. 12 reveals “inconsistencies” on the results. For example, the ω = 20case generally gives smaller errors than the ω = 5 case does in the first regime, and the errors are sometimeseven smaller than those of the ω = 0 case. These inconsistencies originate from the nonuniformity of the20esh in the rotating subdomain. When the vortex travels to the same location, it is actually capturedby different grid resolutions when the rotational speed differs. Nevertheless, for all the cases, the errorsinevitably decrease exponentially, which confirms the high-order accuracy of the sliding-mesh method. Taylor-Couette flow is formed between two concentric rotating circular cylinders. Due to viscous effects,this flow will reach a steady state if the Reynolds number is small. The steady-state azimuthal flow speedhas the following expression, v θ = v θ i r o /r − r/r o r o /r i − r i /r o + v θ o r/r i − r i /rr o /r i − r i /r o , (97)where r i and r o are the radii of the inner and the outer boundaries; v θ i and v θ o are the azimuthal flow speedat these two boundaries.Figure 13 shows the mesh for this simulation. The overall domain is bounded by r i = 1 and r o = 2, andis split into a rotating inner subdomain and a fixed outer subdomain by a sliding interface at r s = 1 .
5. Thesetwo subdomains are meshed into 24 and 32 elements, with uniform grid distribution in azimuthal and radialdirections. The outer boundary is treated as a no-slip isothermal wall with v θ o = 0. The inner boundary Figure 13: Mesh for Taylor-Couette flow simulation. is set as a Dirichlet boundary with ρ i = 1, Mach number Ma = 0 .
1, and v θ i = 1. The Reynolds numberbased on flow properties at the inner boundary is Re = 10. Again, six rotational speeds: ω = 0 , , ,
10, 15and 20, are tested. The same time marching scheme as that for the previous test with a time step size of∆ t = 5 . × − is used for all the tests. Exact transfinite mapping is employed on all circular boundaries,and all interior cell faces are represented linearly.Figure 14 shows the steady state density contours from different polynomials. The contours are expectedto be a series of concentric circles, and we see improving results as the polynomial degree increases. Notethat since the boundary conditions are weakly imposed, the computed values on the inner boundary aretherefore not necessarily identical to the prescribed values. And this is the reason for the dashed lines at theinner boundaries. But the boundary conditions becomes stricter as the polynomial degree increases, whichis why the dashed lines have disappeared at P = 4.We compute the L errors of the u velocity (i.e., the x -component of v θ ) and plot the results in Fig. 15.Overall, the errors decrease exponentially with polynomial degree before temporal errors become dominant,which confirms the high-order accuracy of the present method on viscous flow. Meanwhile, the results fromdifferent rotational speeds are more consistent than those from the previous test. This is because of theuniformity of the current mesh in the azimuthal direction, which gives identical mesh resolution even whenthe rotational speed differs. The effects of rotational speed again agree with our expectation.21 .00301.00271.00241.00211.00181.00151.00121.00091.00061.00031.0000 (a) P = 2 (b) P = 3 (c) P = 4 ρ Figure 14: Steady state density contours of Taylor-Couette flow (dashed lines in the middle represent sliding interface). − − − − − − P L e rr o r ω = 0 ω = 1 ω = 5 ω = 10 ω = 15 ω = 20 Figure 15: L errors of u velocity against polynomial degree for Taylor-Couette flow. In this section, we test the temporal accuracy and verify that the present sliding-mesh method doesnot affect the order of accuracy of a time marching scheme. We employ several strong-stability-preservingRunge-Kutta (SSPRK) schemes for the temporal discretization. Following the rules in [47], we denote an s -stage p -th order SSPRK scheme as SSP( s, p ). The CFL condition requires small enough time step sizeto stabilize a scheme, which makes it difficult for temporal error to dominate. To alleviate this issue, weemploy the following SSPRK schemes that allow larger CFL numbers (and thus larger time step sizes): theSSP(4,2) scheme from [47], the SSP(8,3) scheme from [41], and the SSP(10,4) scheme from [48].The tests are performed on the previous Euler vortex flow. Whenever it permits, we have varied thepolynomial degree for each case to ensure that the spatial errors are negligibly small compared to thetemporal errors. The results are reported in Fig. 16 for different rotational speeds. For almost all the cases,the temporal errors decrease at the correct orders as the time step size decreases. For the w = 0 case withSSP(10,4), the temporal errors become so small (in the order of 10 − ) that they can no longer be easilyseparated from the spatial errors at small time step sizes, which results in the incorrect slope at small timestep sizes. Nevertheless, at large time step sizes, the curve still shows the correct slope for this case. Wealso have tested the temporal accuracy on the Taylor-Couette flow, and similar results were obtained (forconciseness, the results are not included here). 22 . . .
25 1 . . − − − − − − − ∆ t × − L e rr o r . . .
25 1 . . − − − − − ∆ t × − L e rr o r . . .
25 1 . . − − − − − ∆ t × − L e rr o r . . .
25 1 . . − − − − − ∆ t × − L e rr o r . . .
25 1 . . − − − − − − ∆ t × − L e rr o r . . .
25 1 . . − − − − − − ∆ t × − L e rr o r (a) ω = 0 (b) ω = 1 (c) ω = 5(e) ω = 10 (f) ω = 15 (g) ω = 20 SSP(4,2) SSP(8,3) SSP(10,4) slope = 2 slope = 3 slope = 4
Figure 16: L errors of ρ against time step size for an Euler vortex flow. From the analysis in Sec. 4.4, global conservation on a mesh with N elem elements, and N boun1 and N boun2 boundary cell faces that are mapped to the ξ and the η directions, respectively, can be expressed as N elem (cid:88) e =1 (cid:90) (cid:90) ∂ (cid:101) Q e ∂t d ξ d η + N boun1 (cid:88) b =1 (cid:90) sign b · (cid:98) F b d η + N boun2 (cid:88) b =1 (cid:90) sign b · (cid:98) G b d ξ = , (98)where sign = 1 or −
1, depending on which standard face the physical cell face is mapped to. Denoting thetemporal discretization by (cid:101) ∂/ (cid:101) ∂t , then the exact time derivative in (98) can be expressed as ∂ (cid:101) Q e ∂t = (cid:101) ∂ (cid:101) Q e (cid:101) ∂t + (cid:15) (∆ t, ξ, η ) , (99)where (cid:15) (∆ t, ξ, η ) represents a small error that is a function of time-step size and space. Substitute the aboverelation into (98), N elem (cid:88) e =1 (cid:90) (cid:90) (cid:101) ∂ (cid:101) Q e (cid:101) ∂t d ξ d η + N boun1 (cid:88) b =1 (cid:90) sign b · (cid:98) F b d η + N boun2 (cid:88) b =1 (cid:90) sign b · (cid:98) G b d ξ = E (∆ t ) , (100)where E (∆ t ) = − N elem (cid:88) e =1 (cid:90) (cid:90) (cid:15) e (∆ t, ξ, η )d ξ d η. (101)When ∆ t → E (∆ t ) → E (∆ t ) measure how well conservation is satisfied for mass,momentums and energy. 23ased on the above analysis, we employ a flat-plate Couette flow to verify the conservation property ofthe present method. This flow is formed between two infinite flat plates separated by a distance of H , withthe upper plate moving at a constant speed U and the lower plate fixed. The steady state solution is u = U yH , v = 0 , p = constant , (102) T = T + (cid:0) T − T (cid:1) yH + µU κ (cid:18) yH − y H (cid:19) , (103)where 0 ≤ y ≤ H is the vertical coordinate, T and T are temperatures on the lower and the upper plates.The mesh used for this simulation is identical to that in Fig. 10, but is scaled by a factor of 1 /
10 in eachdirection (i.e., with H = 1) to make the simulation setup easier. The upper plate speed is set to U = 1.The temperatures are set to T = T = 1. All the other parameters are chosen such that the flow has aPrandtl number P r = 0 .
72 and Reynolds number Re = 100, and the flow on the upper plate has a Machnumber M a = 0 .
8. The flow field is initialized using the exact solution. The boundary conditions are weaklyimposed also using the exact solution. The SSP(5,4) temporal scheme with a time step size of 1 . × − isused for all the simulations. And all the simulations are performed using double precision float numbers.We have tested different rotational speeds and polynomials, with E (∆ t ) monitored until U t/H = 20. Itwas noticed that for all the cases, E (∆ t ) stays in the level of machine precision and does not grow duringthe simulation. For conciseness, we only report the values of | E (∆ t ) | for some of the tests in Tab. 2 at theend of each simulation. In the calculation, the numerical time derivative term in (100) is replaced by theresidual according to (34), and the integrations are evaluated exactly using quadratures. ω P Mass Momentum 1 Momentum 2 Energy0 4 1.300E-16 3.141E-16 4.935E-16 2.386E-168 2.428E-16 1.103E-16 8.592E-16 3.013E-1712 2.182E-16 5.174E-16 1.699E-15 7.907E-165 4 1.490E-17 1.837E-15 1.078E-16 2.797E-158 4.948E-16 9.753E-16 9.121E-17 8.857E-1612 7.329E-16 7.710E-16 1.142E-15 4.289E-1510 4 2.308E-15 3.565E-15 1.919E-16 3.406E-158 2.595E-16 2.406E-16 1.905E-16 1.256E-1512 1.375E-15 1.132E-15 5.528E-16 8.048E-1520 4 1.920E-15 1.615E-15 1.151E-15 3.046E-158 1.606E-15 1.183E-16 3.572E-16 8.508E-1512 5.550E-15 9.471E-16 1.081E-15 8.013E-15
Table 2: Conservation of mass, momentum and energy on a flat-plate Couette flow.
In this test, the free stream flow is chosen such that it has a Mach number
M a = 0 .
3. The same meshas that for the conservation test are used for this test. Dirichlet boundary conditions using the constantflow are applied at all the outer boundaries. Different polynomials and rotational speeds are tested toinvestigate their effects on free-stream preservation. For all the cases, we employ SSP(10,4) with a time-step size of 1 . × − as the time marching scheme. We measure the free-stream preservation by thenormalized L error of pressure (represents the average strength of numerically generated disturbances). Asa comparison, we also test free-stream preservation on a dynamic conforming mesh that has almost exactlythe same resolution as the nonconforming sliding mesh. The center of conforming mesh has a displacementof ∆ y ( t ) = 0 . t with respect to its initial position, with the mesh in the rest of the domain deformingusing an algebraic function [34]. For all the cases, we monitored the L errors until the nondimensional timereaches U ∞ t/H = 20. The results at the end of the simulations are summarized in Tab. 3.24or perfect free-stream preservation, the error should stay in the level of machine precision, and its mag-nitude should not change dramatically with polynomial degree. Table 3 demonstrates that the FR methodtogether with the GCL equations ensure very good free-stream preservation on the dynamic conformingmesh. The sliding mesh method, however, does not perfectly satisfy free-stream preservation. Instead, itapproaches free-stream preserving in an exponential manner as the polynomial degree increases. ω P ∗ Table 3: Free stream preservation on a dynamic conforming mesh (marked by ‘ ∗ ’) and a sliding mesh (the rest). The reason for the present method not strictly satisfying free-stream preservation is evident. For aconstant free stream flow, we are basically projecting the metric terms between cell faces and mortars.When a circular cell face is represented exactly, its metric terms are in the space of P ∞ . But the output ofa projection is in the space of P N − . This means that there is always a truncation error of O ( ξ N ) duringthe projection of metric terms, which therefore results in an exponential approximation of free-streampreservation.The issue of using exact geometric expression for calculating metric terms has been discussed in [12, 49],and the suggestion is to use polynomial approximation to satisfy free-stream preservation. This suggestionhas been proven to work perfectly on curved conforming meshes. A recent work [50] on curved nonconformingmeshes in the special scenario of subdivision of parent element reveals that a necessary condition for free-stream preservation is that the metrics of a child face being computed from its parent face. But for generalcurved nonconforming meshes, such as the sliding meshes in this work, a child face (e.g., a mortar) does nothave a unique parent (more specifically, a child face has two parent faces that are different polynomials),and thus the aforementioned necessary condition generally can not be satisfied. To the authors’ knowledge,free-stream preservation on general curved nonconforming meshes still remains an open problem for allpolynomial-based high-order methods. Note that for linear nonconforming meshes, such as that in [37], free-stream preservation can always be easily satisfied. But linear mesh obviously introduces too large geometricerror for the curved sliding meshes here (see Appendix A).Nevertheless, Tab. 3 shows that the free-stream-preservation errors from the sliding mesh quickly de-creases to machine precision even at a moderate polynomial degree of 8. Since this test is done on a verycoarse grid, we thus expect very minor free-stream preservation effects in real applications where the meshesare usually finer. As shown in Fig. 17, we employ two types of meshes for this simulation: a sliding mesh (on the left) anda rigid-rotating mesh (on the right). The square cylinder has a diameter (i.e., diagonal length) of D = 1.The sliding mesh consists of 6,797 quadrilateral elements in an overall 100 D × D domain, with 95 of theelements inside a small rotating subdomain whose diameter is 1 . D . The rigid-rotating mesh consists of6,841 elements on a whole-piece rotating circular domain whose diameter is 100 D . Both meshes are refinedin the vicinity of the cylinder to better resolve flow structures. With similar amounts of mesh elements,the sliding mesh obviously provides good resolution in a wider range of the wake region. The rigid-rotatingmesh, on the other hand, always wastes a large amount of elements in unimportant regions, and this isunavoidable. Close views of the meshes are shown on the top right of each figure.25
20 0 20 40 60-40-2002040 -40 -20 0 20 40-40-2002040
Figure 17: Meshes for a rotating square cylinder: left, sliding mesh; right, rigid-rotating mesh.
The freestream flow is chosen such that it has a Mach number
M a = 0 .
1. The Reynolds numberbased on free-stream flow properties and the cylinder diameter is Re D = 100. For the sliding mesh, theinner subdomain rotates counterclockwise at a non-dimensional rotational speed of ωD/U ∞ = π/ T ∗ = 4); and for the rigid-rotating mesh, the whole domain rotatesat this speed. No-slip adiabatic wall boundary condition is applied on the cylinder surface. Characteristic far-field boundary conditions are applied at all outer boundaries. The eighth-order (i.e., P = 7) spatial schemeand the SSP(10,4) time marching scheme [48] with a nondimensional time-step size of ∆ tU ∞ /D = 2 . × − are employed to run the simulations. Polynomial and mesh refinement studies have also been performed tohave confirmed that the present setup well resolves the flow.The lift and drag coefficients of the cylinder are plotted in Fig. 18. The cylinder experiences a negativelift and positive drag all the time. These curves overall do not have a periodic pattern that is directly relatedto the cylinder’s motion. But the horizontal distance between two neighboring local peaks (or troughs) isapproximately 1 (i.e., T ∗ / In this test, we simulate the interactions of a uniform incoming flow and four rotating square cylindersto show the capability of the method for handling multiple rotating objects. A schematic of the simulationsetup is shown in Fig. 20. The four cylinders are separated by a horizontal distance of 1 . D and a verticaldistance of 2 D , where D is the diameter of the cylinder with the same definition as that in the previous test.Two cases are investigated, where the only difference is the rotational direction of each cylinder as markedin the schematic. For simplicity, we denote the four cylinders as A , A , B and B .The overall mesh for this test has 7,409 elements, with its topology and distribution similar to those ofthe sliding mesh in the previous test. Each inner subdomain mesh is exactly the same as that in the previoustest. For conciseness, views of the mesh are not repeated here. All the other parameters for this test, suchas the freestream Mach number, the Reynolds number, the rotational speed, the boundary conditions, etc.,are exactly the same as those for the previous test. The simulations are performed using the eighth-orderscheme with SSP(10,4) and a nondimensional time step size of 1 . × − . Each simulation is continuedfor a nondimensional time of 500. In what follows, we briefly discuss the results for each case.26
26 28 30 32 34 36 38 40 42 44 46 48 50 52 540.81.01.21.4 U ∞ t/D C L U ∞ t/D C D sliding ◦ rigid-rotatingsliding ◦ rigid-rotating Figure 18: Lift and drag coefficients of a rotating square cylinder using a sliding mesh and a rigid-rotating mesh. (a) sliding(b) rigid-rotating
Figure 19: Vorticity contours of flow over a rotating square cylinder (solid lines are positive, dashed lines are negative).
Figure 21 shows the drag and lift coefficients of the cylinders from U ∞ t/D = 50 to 100. It is seen thatcylinders from the same column experience the same drag in this time frame. In fact, we have observedfrom the simulation that this relation holds until U ∞ t/D approaches 130 when the curves start differing(not shown here, reasons explained later). The upstream cylinders experience much larger drag than thedownstream ones do. The downstream ones, however, experience much larger drag oscillations, which is27 D . D D/ √ A A B B ω ωω ω flow (a) Case 1 D . D D/ √ A A B B ω ωω ω flow (b) Case 2 Figure 20: Schematic of simulation setup for flow over multiple rotating square cylinders. possibly due to vortex shedding. On the other hand, due to the Magnus effects and the opposite rotationaldirections, the lifts on A and B (also A and B ) have the same magnitude but opposite signs. Again, theupstream cylinders experience smaller lift oscillations than the downstream cylinders do.
50 60 70 80 90 1000.00.51.01.5
50 60 70 80 90 100-1.00.01.02.0 U ∞ t/D C D U ∞ t/D C L A A B B A A B B Figure 21: Drag and lift coefficients for flow over multiple rotating square cylinders (case 1).
A snapshot of the flow field visualized by vorticity contours at U ∞ t/D = 80 . U ∞ t/D = 130) before breaking up. This case evidently demonstrates that the present method introducesvery little numerical disturbance to a simulation. After the breakup of the symmetry, the flow becomes lessinteresting, and that is why we only focus on the flow before that point. Figure 22: Vorticity contours of flow over multiple rotating square cylinders (case 1).
The drag and lift coefficients for this case are plotted in Fig. 23. The previous conclusions for Case 1still hold here. Except that the current curves show a rather mono frequency that corresponds to a periodof T ∗ /
4, where T ∗ is the nondimensional rotational period of the cylinders. This mono frequency indicatesfewer vortical structures (or noises) in the flow field. Compared with Case 1, we also notice dramatic dragreduction on all the cylinders. On average, the rear cylinders even do not experience much drag at all.Meanwhile, dramatic lift enhancement is observed on all the cylinders.
80 85 90 95 1000.00.51.0
80 85 90 95 100-1.00.01.0 U ∞ t/D C D U ∞ t/D C L A A B B A A B B Figure 23: Drag and lift coefficients for flow over multiple rotating square cylinders (case 2).
Figure 24 shows the instantaneous flow field at U ∞ t/D = 250 .
5. Compared with Case 1, vortex sheddinghas been completely suppressed in this case, which is consistent with our observation form the force curves.29his flow remains very stable and is almost steady (except in the very vicinity of the cylinders) throughoutthe simulation (i.e., even at U ∞ t/D = 500). The reason for this steadiness is obvious. The rotation motionof the cylinders in this case (as shown in Fig. 20(b)) decelerates the flow (also increases pressure) in the gapbetween the two rows. This reduction of momentum and increase of pressure prevent vortex shedding fromthe cylinders and thus make the flow almost steady. Figure 24: Vorticity contours of flow over multiple rotating square cylinders (case 2).
6. Summary
We have presented a high-order sliding mesh method and performed detailed studies of its properties.This method employs transfinite mortar elements for communication between the two sides of a nonconform-ing sliding interface. Solutions and fluxes are projected back and forth in a least-squares fashion betweenmortars and cell faces to retain many of the favorable properties of the original conforming-mesh method.The use of transfinite mortar elements completely eliminates geometric errors on a circular sliding interface,and thus excludes an important source of errors from the scheme. Numerical tests have verified that thismethod retains the high-order accuracy of the FR method on both inviscid and viscous flows as the errorsdecrease exponentially with the increase of polynomial degrees. Numerical tests also have revealed that thismethod retains the order of accuracy of a time marching scheme, which allows the minimization of rotationalspeed effects by equipping the method with high-order time marching schemes. Through numerical analy-sis and verifications, we have confirmed the conservation property of the method, even under dramaticallydifferent rotational speeds. The satisfaction of outflow condition is also proved and verified, which ensuresthat, using the present method, a nonconforming sliding interface does not affect the characteristics of wavesin a flow field. This method is, however, not strictly free-stream preservative. Instead, it approaches free-stream preservation in an exponential way. This is in fact not specific to the present method, as free-streampreservation on general curved nonconforming mesh still remains an open problem in the high-order methodcommunity. Nevertheless, we have shown through numerical tests that the free-stream-preservation errorsare negligibly small even on a very coarse mesh using polynomials of moderate degrees. The comparisonstudy of sliding mesh and rigid-rotating mesh on a rotating square cylinder has shown that the two ap-proaches generate almost identical results. But a sliding mesh generally provides better wake resolutionthan a rigid-rotating mesh for a given number of overall mesh elements. The simulation of a matrix ofrotating square cylinders has demonstrated the method’s capability of simultaneously handling multipleobjects. As a concluding remark, we must also emphasize that this method is not specifically designed forthe FR method, it is readily applicable to many other high-order methods as well. This method is alsoeasily extensible to deal with three-dimensional problems.
Acknowledgment
This work was supported by a grant from the Office of Naval Research, administrated by Dr. Ki-HanKim. The authors would like to express our acknowledgments for the support.30 ppendix A. Comparison of iso-parametric mapping and transfinite mapping
Transfinite mapping can minimizes geometric errors, whereas iso-parametric mapping always carries atruncation error but monolithically approaches transfinite mapping as the order increases. We comparethese two types of mappings through a 1D and a 2D examples. The 1D example is a circular arc with θ = 0 ◦ , θ = 10 ◦ , R = 1, and x c = (0 , R between the exact radiusand the numerical ones that are calculated using the coordinates from the mappings. Figure 25 shows thatthe iso-parametric mapping represents the circular arc exactly only at the nodal points. For non-nodalpoints, the error overall decreases as the order increases. On the other hand, the transfinite mapping alwaysrepresents the circular arc exactly with the errors alway in the level of machine precision. . . . . − − − − − − − − ξ | ∆ R | linear quadraticcubic transfinite Figure 25: Comparison of iso-parametric mapping and transfinite mapping for approximating a circular arc.
The 2D geometry is similar to that in Fig. 5, i.e., with three straight edges and one circular-arc edge.More specifically, we have chosen: x = ( R cos θ + x c , R sin θ + y c ), x = ( R cos θ + x c , R sin θ + y c ), x = (( R/
2) tan θ , R/ x = (( R/