A high order Discontinuous Galerkin – Fourier incompressible 3D Navier–Stokes solver with rotating sliding meshes
Introduction
Problems where the forces on rotating/oscillating geometries are to be predicted are common in engineering and fluid–structure interaction problems. Examples are flows around isolated rotating bodies and airfoils, turbomachinery applications, insect flight aerodynamics, unmanned air vehicles and more recently flows around renewable energy devices; wind and tidal turbines. These flows are characterised by long wake structures, vortex shedding and stalled flows associated with flow unsteadiness. To accurately simulate wakes and vortex structures and their evolution, high order (typically ⩾3) numerical methods (i.e. h/p conformal spectral and h/p non-conformal Discontinuous Galerkin methods) are preferred since dissipation and dispersion errors are minimum [1] when compared to low order (typically ⩽3) methods. Furthermore, for smooth problems, the exponential decay of the error with polynomial enrichment (p-refinement) in high order methods as opposed to the fixed decay rate characteristic of low order methods (i.e. h-refinement only) renders high order methods particularly attractive to obtain accurate solutions for flows where viscosity limits discontinuities (i.e. elliptic type equations) as in the incompressible Navier–Stokes (NS) equations.
For conformal discretisations (i.e. classic low order finite elements or h/p spectral type methods), if geometric discontinuities arise (e.g. hanging nodes), a geometric incompatibility [1] is created since C0 continuity cannot be ensured across elements. Similarly, when neighbouring elements have different numbers of degrees of freedom (e.g. mixed triangular–quadrilateral mesh), a functional incompatibility [1] is created. Over the past few decades methods to overcome these difficulties have been developed in the framework of conformal discretisations, examples are iterative patching, constrained approximation or mortar patching (see [1] for a review). A natural alternative to overcome these incompatibilities is to relax the continuity condition across elements, which leads to the Discontinuous Galerkin (DG) finite element approach.
The DG finite element method can be seen as an extension of h/p spectral methods where the C0 continuity requirement across element boundaries is relaxed or also as a high order finite volume method with compact stencil. Contrary to conformal finite elements or h/p spectral methods, DG methods are locally conservative by construction. As in spectral methods, high-order polynomials can be used within each element allowing exponential convergence, whilst discontinuities in the solution are permitted at element interfaces.
It has been argued that DG methods are prohibitively expensive, when compared to conformal spatial discretisations (e.g. h/p spectral methods), since additional degrees of freedom (DOF) arise from the discontinuities between elements. However, the relative number of the additional boundary degrees of freedom to internal degrees of freedom for each element, decreases rapidly for high polynomials. To exemplify this, let us consider a two dimensional problem and a mesh with Nel triangular elements with polynomial spaces of order k. On the one hand, if a h/p spectral discretisation is considered, the number of global degrees of freedom required is (see [2]). On the other hand, a DG discretisation leads to a global number of DOF: (see Section 2.2.6 for details). The last expression shows that for high polynomial orders k, the number of degrees of freedom for the DG discretisation tends to . This simple analysis shows that for high polynomial orders DOFDG ≈ DOFspectral and the additional cost of DG methods is not overwhelming.
DG formulations to solve the incompressible NS equations have seen increased popularity over recent years as evidenced by the number of publications on the topic, (e.g. [3], [4], [5], [6], [7], [8], [9] and the authors’ [10]). In [10], the authors presented the development of a DG code that uses the Symmetric Interior Penalty Galerkin (SIPG) formulation for solving the unsteady 2D NS equations using straight sided triangular elements. Simulation results from the solver were shown to be in good agreement with experimental results and results from a h/p conformal spectral code. The present implementation is an extension of the previous work to enable solutions on hybrid meshes (i.e. mixed triangular and quadrilateral elements) with rotationally sliding meshes and with curved boundary conditions. In addition, Fourier series applied orthogonally to the 2D DG plane are used to enable 3D flow solutions for problems with spanwise geometric homogeneity. These developments have been implemented and tested and results are reported herein.
The sliding mesh technique allows for mesh motion where an inner mesh zone region rotates with respect to an outer static mesh. This relative motion creates hanging nodes at the interface between static and rotating elements and necessitates the use of elements with internal curved edges. In addition, to achieve high order solutions near external boundaries (i.e. walls) curved edged elements are essential. To clarify these concepts before continuing, we depict in Fig. 1 an example of a 2D mesh, where the static and rotating subdomains, the external curved boundary for a symmetric NACA airfoil and the curved sliding mesh interface with the associated hanging nodes have been highlighted.
The present work shows that in the DG context, the geometric incompatibility arising from the hanging nodes due to mesh rotation does not cause loss of exponential convergence properties. Further, we show that the functional incompatibility originating from the use of triangular and quadrilateral element types (these elements have different numbers of local degrees of freedom) in combination with orthogonal modal basis functions, does not present a problem. Finally, the sliding mesh implementation shows high accuracy when solving the incompressible NS equations in two and three dimensions.
To account for the relative mesh movement of the inner mesh with respect to the static outer mesh, it is advantageous to write the equations for fluid motion in their Arbitrary Lagrangian–Eulerian (ALE) form [11]. The ALE description was first introduced for finite difference methods and unstructured meshes for fluid simulation [12] and subsequently extended to finite elements [13]. A wide explanation of this method can be found in [11].
The ALE approach is generally used for dynamically deforming mesh elements (i.e. arbitrary node movement) and has been widely explored in combination with DG methods to solve hyperbolic type equations including the compressible NS equations (see for example [14], [15], [16], [17]. As for the incompressible ALE form of the NS equations with arbitrary mesh movement, work was limited for a long time to conformal h/p spectral discretisations [18]. However, very recent work combines this technique with a DG approach [9].
Deforming element techniques (e.g. ALE for deforming elements) require either generally small body motions or remeshing for large motions that would have otherwise lead to unacceptably element distortions. To avoid these limitations, an appealing approach is provided by combining the ALE approach with the sliding mesh technique. This method is particularly suitable to problems where the mesh movement is known a priori; e.g. rigid body rotation without mesh deformation. We chose to follow this approach, the ALE formulation with sliding meshes for non deforming elements, and summarise some of its advantages below:
- •
Inertially fixed and rotating objects can be present in the same simulation (see examples in Sections 3.3 An application: three bladed cross-flow turbines, 4 Three dimensional flows).
- •
No remeshing is required (with its associated computational cost) as mesh elements do not distort, enabling unlimited rotation as opposed to non-sliding or deforming ALE methods in which large distortions and subsequent remeshing need to be handled.
- •
Subsequently, no projection of the solution into a new mesh is required to advance the solution in time, which is generally a non-conservative process.
- •
In our approach, no interpolation through the sliding interface is required (as in low order methods) which would introduce large numerical errors, destroying the high order properties of the method.
The ALE approach with sliding mesh interfaces, has been previously studied in the context of h/p conformal spectral methods that require mortar techniques for subdomain linking [19]. Recently, a sliding mesh capability has been described for the incompressible NS equations [20], where the authors used isogeometric analysis and a low order conformal discretisation (using NURBS) coupled with a special treatment for the sliding curved interface exploiting ideas from DG methods. Purely DG discretisations in conjunction with the sliding mesh approach, have been used for the solution of electromagnetic problems (i.e. Maxwell equations) [21], [22]. In [22], it was shown that to solve the Maxwell eigensystem with sliding mesh interfaces, mortar techniques are beneficial.
To the authors’ knowledge, the work presented herein details the first high order (order ⩾ 3) DG solver with sliding meshes for the solution of the incompressible NS equations. We present a novel approach where non-conformal DG is used in all elements and the curved sliding mesh interface is approximated through an analytical mapping for the description of the circular interface. We show that our approach and implementation does not require of mortar type techniques to accurately solve the incompressible NS equations. In addition, to account for curved external boundary surfaces (i.e. walls) representing NACA 4-digit airfoils, we introduce an analytical mapping as is required to obtain smooth solutions for high order techniques. An efficient three dimensional extension of the 2D DG solver is included, which enables 3D flow solutions on rotating geometries that present a geometrical homogeneity in the spanwise (Fourier) direction.
The organisation of the paper is as follows. We firstly introduce, in Section 2, the methodology used to discretise the equations in time and space, with particular emphasis on the required mappings to account for the curved edges (Section 2.2.4), the discretisation and implementation of surface integrals on the curved sliding mesh (Section 2.2.5) and the 3D extension using Fourier series (Section 2.3). Test cases to assess the accuracy, convergence properties and usability of the solver can be found in Sections 3 Two dimensional flows, 4 Three dimensional flows. Namely, two dimensional flows are verified through the Taylor vortex problem, a static and rotating NACA0015 airfoil and a cross-flow turbine in Section 3. Section 4 concludes with two applications for three dimensional flows: a three bladed cross-flow turbine and a circular cylinder shadowed by a dynamically pitching NACA0012 airfoil.
Section snippets
Methodology
We are interested in solving the incompressible NS equations in primitive variable formulation in a unique domain composed of non overlapping static and moving subdomains: Ω(t) = Ωsta ∪ Ωrot(t) where Ωsta ∩ Ωrot(t) = Γsta−rot and Γsta−rot defines the circular sliding mesh interface. The non-dimensional incompressible NS equations in ALE form [11] can be written (in the absence of body forces) as:where u = (u, v)T and w are non-dimensionalised using the characteristic
Two dimensional flows
This section is devoted to code verification and determining the convergence properties of the described method in two dimensions (i.e. x–y DG plane). We first select the Taylor vortex problem [1], where the analytical solution is known to explore the numerical characteristics of triangular–quadrilateral element mixing, internal curved edges and sliding mesh interfaces. Secondly, we compute solutions on a static (non-rotating) NACA0015 at zero angle of attack to explore the solver convergence
Three dimensional flows
In this last section, we provide two examples of three dimensional flow simulations using the Fourier extension of the DG solver with rotating sliding meshes.
Firstly, the three bladed cross-flow turbine geometry described in Section 3.3 is simulated with a rotational speed of Lω/U = 0.1 at a Reynolds number (based on the freestream velocity and blade chord) of Re = 200 (i.e. large enough to trigger three dimensionalities in the flow). The spatial discretisation consists of polynomials of order k = 4
Conclusions
This work details the formulation, implementation and verification of a novel approach to solve the three dimensional incompressible NS equations using a high order DG-Fourier formulation with rotating sliding meshes. The solver allows accurate solutions with curved elements that conform to either arced surfaces (internal or external) and the profile of NACA 4-digit airfoils. Hybrid unstructured meshes (i.e. with triangular and quadrilateral elements) can be handled without damaging the
Acknowledgement
The authors thank the financial support of the John Fell Oxford University Press (OUP) fund as well as the Research Council UK (RCUK).
References (36)
- et al.
From h to p efficiently: implementing finite and spectral h/p element methods to achieve optimal performance for low- and high-order discretisations
Journal of Computational Physics
(2010) - et al.
An artificial compressibility flux for the discontinuous Galerkin solution of the incompressible Navier–Stokes equations
Journal of Computational Physics
(2006) - et al.
A high-order discontinuous Galerkin method for the unsteady incompressible Navier–Stokes equations
Journal of Computational Physics
(2007) - et al.
An implicit high-order hybridizable discontinuous Galerkin method for the incompressible Navier–Stokes equations
Journal of Computational Physics
(2011) - et al.
A space–time hybridizable discontinuous Galerkin method for incompressible flows on deforming domains
Journal of Computational Physics
(2012) - et al.
A high order Discontinuous Galerkin finite element solver for the incompressible Navier–Stokes equations
Computers & Fluids
(2011) - et al.
An arbitrary Lagrangian–Eulerian computing method for all flow speeds
Journal of Computational Physics
(1974) - et al.
Lagrangian–Eulerian finite element formulation for incompressible viscous flows
Computer Methods in Applied Mechanics and Engineering
(1981) - et al.
A discontinuous Galerkin ALE method for compressible viscous flows in moving domains
Journal of Computational Physics
(1999) - et al.
Space–time discontinuous Galerkin method for nonlinear water waves
Journal of Computational Physics
(2007)