A high order Discontinuous Galerkin – Fourier incompressible 3D Navier–Stokes solver with rotating sliding meshes

https://doi.org/10.1016/j.jcp.2012.04.039Get rights and content

Abstract

We present the development of a sliding mesh capability for an unsteady high order (order  3) h/p Discontinuous Galerkin solver for the three-dimensional incompressible Navier–Stokes equations. A high order sliding mesh method is developed and implemented for flow simulation with relative rotational motion of an inner mesh with respect to an outer static mesh, through the use of curved boundary elements and mixed triangular–quadrilateral meshes.

A second order stiffly stable method is used to discretise in time the Arbitrary Lagrangian–Eulerian form of the incompressible Navier–Stokes equations. Spatial discretisation is provided by the Symmetric Interior Penalty Galerkin formulation with modal basis functions in the xy plane, allowing hanging nodes and sliding meshes without the requirement to use mortar type techniques. Spatial discretisation in the z-direction is provided by a purely spectral method that uses Fourier series and allows computation of spanwise periodic three-dimensional flows. The developed solver is shown to provide high order solutions, second order in time convergence rates and spectral convergence when solving the incompressible Navier–Stokes equations on meshes where fixed and rotating elements coexist.

In addition, an exact implementation of the no-slip boundary condition is included for curved edges; circular arcs and NACA 4-digit airfoils, where analytic expressions for the geometry are used to compute the required metrics.

The solver capabilities are tested for a number of two dimensional problems governed by the incompressible Navier–Stokes equations on static and rotating meshes: the Taylor vortex problem, a static and rotating symmetric NACA0015 airfoil and flows through three bladed cross-flow turbines. In addition, three dimensional flow solutions are demonstrated for a three bladed cross-flow turbine and a circular cylinder shadowed by a pitching NACA0012 airfoil.

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 DOFspectralNel2k2 (see [2]). On the other hand, a DG discretisation leads to a global number of DOF: DOFDGNel2(k+1)(k+2)=Nel2k21+3k+2k2 (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 DOFDGNel2k2. 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) = Γstarot and Γstarot 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:ut+((u-w)·)u=-p+1Re2u;·u=0inΩ(t),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. xy 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 /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)

  • V. Nguyen

    An arbitrary Lagrangian–Eulerian discontinuous Galerkin method for simulations of flows over variable geometries

    Journal of Fluids and Structures

    (2010)
  • A. Beskok et al.

    An unstructured hp finite-element scheme for fluid flow and heat transfer in moving domains

    Journal of Computational Physics

    (2001)
  • G.E. Karniadakis et al.

    High-order splitting methods for incompressible Navier–Stokes equations

    Journal of Computational Physics

    (1991)
  • R. Kirby et al.

    De-aliasing on non-uniform grids: algorithms and applications

    Journal of Computational Physics

    (2003)
  • R. Sevilla et al.

    Numerical integration over 2D NURBS-shaped domains with applications to NURBS-enhanced FEM

    Finite Elements in Analysis and Design

    (2011)
  • G. Karniadakis

    Spectral element-Fourier methods for incompressible turbulent flows

    Computer Methods in Applied Mechanics and Engineering

    (1990)
  • G. Karniadakis et al.

    Spectral h/p Element Methods for Computational Fluid Dynamics

    (2005)
  • B. Riviere

    Discontinuous Galerkin methods for solving elliptic and parabolic equations: theory and implementation

    SIAM-Frontiers in Applied Mathematics

    (2008)
  • Cited by (0)

    View full text