Elsevier

Journal of Computational Physics

Volume 307, 15 February 2016, Pages 60-93
Journal of Computational Physics

A spectrally accurate method for overlapping grid solution of incompressible Navier–Stokes equations

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

Highlights

  • A novel overlapping grid method for incompressible flow equations is developed.

  • The overlapping grid method is based on spectral-element discretization.

  • The method maintains global spectral accuracy with polynomial refinement.

  • An explicit interface extrapolation scheme for temporal coupling is proposed.

  • The scheme allows for high-order temporal accuracy with only a few iterations.

Abstract

An overlapping mesh methodology that is spectrally accurate in space and up to third-order accurate in time is developed for solution of unsteady incompressible flow equations in three-dimensional domains. The ability to decompose a global domain into separate, but overlapping, subdomains eases mesh generation procedures and increases flexibility of modeling flows with complex geometries. The methodology employs implicit spectral element discretization of equations in each subdomain and explicit treatment of subdomain interfaces with spectrally-accurate spatial interpolation and high-order accurate temporal extrapolation, and requires few, if any, iterations, yet maintains the global accuracy and stability of the underlying flow solver. The overlapping mesh methodology is thoroughly validated using two-dimensional and three-dimensional benchmark problems in laminar and turbulent flows. The spatial and temporal convergence is documented and is in agreement with the nominal order of accuracy of the solver. The influence of long integration times, as well as inflow–outflow global boundary conditions on the performance of the overlapping grid solver is assessed. In a turbulent benchmark of fully-developed turbulent pipe flow, the turbulent statistics with the overlapping grids is validated against published available experimental and other computation data. Scaling tests are presented that show near linear strong scaling, even for moderately large processor counts.

Introduction

Finding numerical solutions to partial differential equations (PDEs) by decomposing the computational domain into smaller subdomains is an idea that has been around for well over a century. Domain decomposition methods have been utilized for several different purposes, including straightforward parallelization [1], [2], [3], [4], simplified mesh generation for complex geometries [5], [6], [7], [8], and the ability to use different parameters or methods in different subdomains [9], [10], [11], [12]. These techniques exist in many forms, and each has its strengths. Some decompose the global domain into overlapping subdomains [13], [14], [15], [16], [17], [18], while others employ non-overlapping subdomains [19], [9], [7], [20], [21], [22]. Some use explicit interpolation techniques for values at interface boundaries [9], [23], [20], [24], and others carry out implicit interpolation [25], [26], [22], [27], [13], [14], [28], [29], [30]. Domain decomposition techniques have been developed for use with several numerical methods including finite difference [14], [5], [9], [23], finite element [25], finite volume [17], [31], and spectral methods [21], [22].

The earliest known research in domain decomposition methods was performed by H.A. Schwarz whose work was published in 1870 [13]. The original Schwarz Alternating Method, initially proposed for analytical calculations [32], was developed for the global solution of boundary value problems for harmonic functions [15] decomposed into overlapping subdomains, Ω=Ω1Ω2. The solution in the first subdomain (Ω1 with boundaries Ω1Ω and Γ1=Ω1Ω) is found using the global boundary conditions on Ω1Ω and corresponding values from Ω2 at the previous iteration on Γ1. The solution of Ω2 is then found by using values from the solution in Ω1 on Γ2. These two steps are iterated until sufficient convergence is reached (see [33], [32], [34]).

In the 1960's, Volkov generalized the original Schwarz Alternating Method into a numerical domain decomposition technique, in a form of the Composite Mesh Difference Method (CMDM) [14]. CMDM used finite difference methods to solve the 2-dimensional Poisson equation numerically on overlapping grids. His research laid the foundation for subsequent techniques that extended the use of CMDM to other elliptical and hyperbolic PDEs, and boundary value and initial value problems, with the ability to use curvilinear meshes (see [35], [36], [37], [38], [6]). Overlapping domain decomposition methods have also been developed to model complex equations and handle various difficulties in solving practical problems. The Chimera Grid Scheme, introduced in [39], employs overset (overlapping) grids for modeling flows in complex geometries. Shortly after initial development it was enhanced for use with three dimensional flows modeled by the Euler equations [40] and later with the addition of the thin-layer Navier–Stokes equations. More recently, Chimera Grid techniques have been used to model a variety of problems with complex geometries [41], [42]. Subsequently, Henshaw and Schwendeman [17], [4] developed a method for using overlapping mesh techniques in modeling high-speed reactive flows, in two and three dimensions.

In addition, techniques that employ non-overlapping grids (sometimes called patched grids) were developed. Examples include a zonal approach that uses a flux-vector splitting technique for the determination of interface values in Euler equations [43], [44], [45], [46], Lions method [19] that uses an iterative technique to arrive at the correct values to be passed between non-overlapping subdomains in solving Laplace's equation and more general second-order elliptic problems, Dawson's approach [9] that solves the heat equation using an explicit finite difference formula to determine the interface values and allows for different time stepping to be used in different subdomains. Non-overlapping grid techniques have also been extended and employed in solving the advection–diffusion equation [20] and the Navier–Stokes equations [45], [22]. Some of the more recently developed non-overlapping domain decomposition methods achieve high finite global accuracy [47] and some spectral accuracy [21], [7], [48], [22].

While non-overlapping mesh methods allow some flexibility in mesh generation, the constraints in these techniques inhibit additional flexibility that is seen in overlapping mesh methods. By allowing variable overlap size, a broad range of potential mesh configurations are supported with overlapping methods, thus allowing for more simplified mesh generation. Additionally, overlapping methods provide a convenient framework for further extension towards moving domain methods, allowing for general and unconstrained motion of rigid body parts through the background stationary meshes [49], [50], [51], [52].

So far, existing overlapping grid methods for the time-dependent PDE coupling have been traditionally relying on low-order, finite-difference or finite-volume schemes. Although some of the methods have been extended to achieve higher-order spatial convergence, using extended stencil finite-difference or compact schemes, the upper bound of the global accuracy has been usually limited to four [23], [24], [53], and at most six [54], [55], [56], [27], [57]. Recently, Brazell, Sitaraman and Mavriplis developed a high-order overlapping Discontinuous Galerkin solver for compressible equations, that uses Lagrangian interpolation at interface boundaries, and documented a polynomial convergence up to fourth order [58]. In the current paper, we introduce a spectrally-accurate overlapping mesh method for incompressible equations, that is based on a spectral-element method. The Spectral Element Method, which can be perceived as a high-order extension of the Finite Element Method, divides a domain into several conforming and adjacent elements [59], [3], [60]. The volume within each element is discretized using Nth-order tensor-product Lagrange interpolating polynomials on Gauss–Lobatto–Legendre nodal points. Approximations in all elements are coupled at the boundaries to form a global solution [60], which achieves spectral convergence with p- (polynomial order) refinement. In this work, we combine a Spectral Element Method solver with the overlapping grid approach, to arrive at a globally spectrally-accurate method for solution of the incompressible Navier–Stokes equations on overlapping domains.

One of the inherent challenges with overlapping grid methods is to minimize the errors that are introduced due to the coupling of the individual subdomain solutions into the global solution. The coupling errors consist of spatial errors and temporal errors, and have to be treated separately. Spatial errors are introduced by the spatial interpolation stencil employed to obtain a function value at the interface points of one domain from the gridpoint values in the adjacent domains at the same time level. Some overlapping mesh methods circumvent the spatial error by requiring that the gridpoints in overlapping domains exactly coincide [61], [62], thus fully conserving communicated information, with the drawback of decreased flexibility in mesh generation. Other methods that do not require the exact match of the gridpoints and thus are more flexible, use finite order interpolation schemes to determine values from adjacent subdomains. Although simple linear interpolation techniques have been popular [6], [9], [63], [29], [64], it was shown by Chesshire and Henshaw [6] that an interpolation scheme should be consistent with the accuracy of the underlying solver and higher-order interpolation is required to maintain the accuracy of the coupled solution with high-order methods, generally, of the order of one higher than the underlying solver accuracy if the overlap width scales with the grid resolution, and the same if it stays constant during grid refinement. Thus, in fourth- and sixth-order methods [58], [23], [27], [24], a generalized Lagrangian interpolation method of the corresponding order was employed. In the present work, we use Nth-order Lagrangian interpolation approach, on non-uniform Gauss–Lobatto–Legendre grids, to arrive at a spectrally-accurate interpolation scheme (with p-refinement) that allows us to maintain the global spectral accuracy of the coupled solution.

To avoid temporal errors due to coupling, the subdomain coupling methods have often been formulated implicitly, when either the interpolation dependencies are written into the global matrix [6], [26], [27], or independent subdomain solutions are coupled through the use of Schwarz-like iterations [13] to ensure interface value convergence [28], [29], [30], [31]. Global convergence of the variables in adjacent subdomains often requires many, sometimes hundreds [29], or even thousands [31] iterations, the number generally being dependent upon the PDE, the overlap size, and timestep [65], [66]. Although Henshaw [67] showed that the convergence rates when solving elliptic PDEs on overlapping subdomains can be reduced by using a multigrid solver with a smoothing algorithm near interface boundaries, this approach is difficult to implement and requires some global interventions to the solver. However, if the coupling scheme, even without iterations, does not introduce temporal errors larger than those of the global timestepping scheme, the global temporal accuracy will be preserved [68], particularly the interface values need not strictly match, but must merely be “consistent”. The notion of consistency has been previously used when constructing and analyzing interpolation and iteration schemes for the overlapping grid methods [30], where consistency was interpreted in a spatial sense. In [68], we have proposed an explicit temporal interface extrapolation algorithm for the overlapping grid methods that preserves the overall temporal accuracy, and analyzed its stability. The proposed coupling scheme is essentially “temporally consistent” (the differences between the interface values due to temporal coupling converge with the same rate as the temporal errors of the underlying flow solver).

In the present work, we implement this temporally consistent coupling scheme for coupling the solutions of incompressible Navier–Stokes equations on overlapping grids discretized with the SEM method. This novel scheme achieves a specified order of accuracy with a small number of iterations. It ensures that the temporal error between the values in the subdomains at the interfaces is equal to the temporal accuracy of the underlying SEM solver. Assuming that this holds true, the interface values need not be an exact match, while still maintaining the global accuracy of the solution. While the accuracy of the extrapolation scheme is maintained in the absence of iterations, the stability of the formulation can generally be improved by implementing a low number of Schwarz-like iterations [68].

The current paper introduces an original method for coupling Navier–Stokes equations on the overlapping domains that is spectrally accurate in space and high-order accurate in time up to a specified order. Although an overlapping Schwarz method has been previously employed as a preconditioner for the solution of the Poisson equation in a spectral element formulation of Navier–Stokes equations on a single domain [3], [69], a spectrally-accurate approach that achieves a full unsteady coupling of independent solutions to Navier–Stokes equations on overlapping grids appears to be novel. In addition, our method achieves a specified temporal accuracy with the minimum number of iterations (usually two to three) while maintaining stability of the formulation. We demonstrate the spectral spatial and high-order temporal convergence of the developed solver, as well as the overall consistency of the coupled solution, on a number of two and three-dimensional examples with incompressible Navier–Stokes equations, including laminar and turbulent test cases.

The remainder of the paper is organized in the following manner: section 2 describes the computational methodology; section 3 discusses the techniques used to parallelize the overlapping grid solver; section 4 presents the results for verification and validation; and section 5 shows the scalability tests. In section 6 some final conclusions are drawn, including a brief discussion of future extensions and applications of this work.

Section snippets

Mathematical formulation

In the present formulation, a two or three dimensional global solution domain is decomposed into two overlapping subdomains, Ω=Ω1Ω2, an example of which is seen in Fig. 1. Values at interface boundaries, Γ12 and Γ21 are equal to contiguous values in the adjacent subdomain, and the solution in each subdomain is governed by the incompressible Navier–Stokes equations which are represented in non-dimensional form in space Rd, as shown belowΩ1{u[1]t+u[1]u[1]=p[1]+1Reu[1]2u[1]=0 andΩ2{u[2]

Parallelization

Communication among separate subdomains is currently accomplished through a dual-session communication framework, which allows for two sessions to perform computations separately while enabling real-time communication with each other in a single computational platform. Each session runs its own independent instance of the SEM solver for its subdomain where MPI intracommunicators are established for communication among the local processors. An intercommunicator is established for global

Results

In this section, we demonstrate performance of the developed overlapping grid method, on two- and three-dimensional laminar and turbulent problems, with the focus on spatial and temporal orders of convergence, accuracy in the presence of outflow disturbances and long integration times, and the ability of the method to reproduce correct turbulence characteristics on overlapping domains.

Scalability

In this section, we perform the scalability tests for the developed overlapping grid method. We utilize the previous test case of the turbulent pipe flow simulations as a platform for the scalability studies. However, the refined computational meshes were constructed, with higher element counts than in the DNS simulations described above (see label to Fig. 29, Fig. 30) in order to enable testing on high processors counts (up to 1024). The refined outside mesh consisted of three times the number

Conclusions

We have developed an overlapping grid methodology for two and three dimensional applications in a Spectral Element Method incompressible flow solver. Spectral accuracy is maintained for spatial discretization due to a spectral interpolation at interface boundaries using solution approximations in an Nth-order polynomial space on GL points. The global high-order (up to a third) temporal accuracy of the flow solver is also maintained with few, or no, iterations using solutions at previous

Acknowledgements

This work has supported by the NSF grants #1250124 and #1025359. Support by the Office of Advanced Scientific Computing Research, Office of Science, U.S. Department of Energy, under Contract DE-AC02-06CH11357 and by the NSF RTG grant at Northwestern University is also acknowledged. We thank the NSF XSEDE program for providing the allocations on the SDSC Gordon and TACC Stampede clusters where the scalability and timing tests were performed.

References (86)

  • T.M. Burton et al.

    Analysis of a fractional-step method on overset grids

    J. Comput. Phys.

    (2002)
  • S.E. Sherer et al.

    High-order compact finite-difference methods on general overset grids

    J. Comput. Phys.

    (2005)
  • P.K. Moore et al.

    Adaptive local overlapping grid methods for parabolic systems in two space dimensions

    J. Comput. Phys.

    (1992)
  • H.S. Tang et al.

    An overset-grid method for 3D unsteady incompressible flows

    J. Comput. Phys.

    (2003)
  • M.M. Rai et al.

    Metric-discontinuous zonal grid calculations using the Osher scheme

    Comput. Fluids

    (1984)
  • K. Hessenius et al.

    Applications of a conservative zonal scheme to transient and geometrically complex problems

    Comput. Fluids

    (1986)
  • D.A. Kopriva

    A staggered-grid multidomain spectral method for the compressible Navier–Stokes equations

    J. Comput. Phys.

    (1998)
  • N.C. Prewitt et al.

    Parallel computing of overset grids for aerodynamic problems with moving objects

    Prog. Aerosp. Sci.

    (2000)
  • G. Houzeaux et al.

    A Chimera method based on a Dirichlet/Neumann (Robin) coupling for the Navier–Stokes equations

    Comput. Methods Appl. Mech. Eng.

    (2003)
  • W.D. Henshaw et al.

    Moving overlapping grids with adaptive mesh refinement for high-speed reactive and non-reactive flow

    J. Comput. Phys.

    (2006)
  • C.K.W. Tam et al.

    A wavenumber based extrapolation and interpolation method for use in conjunction with high-order finite difference schemes

    J. Comput. Phys.

    (2000)
  • A.K. Chaniotis et al.

    High order interpolation and differentiation using B-splines

    J. Comput. Phys.

    (2004)
  • A.T. Patera

    A spectral element method for fluid dynamics: laminar flow in a channel expansion

    J. Comput. Phys.

    (1984)
  • H. Rui

    Multiplicative Schwarz methods for parabolic problems

    Appl. Math. Comput.

    (2003)
  • C. Schneidesch et al.

    Chebyshev collocation method and multi-domain decomposition for Navier–Stokes equations in complex curved geometries

    J. Comput. Phys.

    (1993)
  • T.F. Chan et al.

    Domain decomposition algorithms

    Acta Numer.

    (1994)
  • B.F. Smith et al.

    Domain Decomposition: Parallel Multilevel Methods for Elliptic Partial Differential Equations

    (1996)
  • C. Bernardi et al.

    Domain decomposition by the mortar element method

  • R.L. Meakin

    Composite overset structured grids

  • C. Dawson et al.

    A finite difference domain decomposition algorithm for numerical solution of the heat equation

    Math. Comput.

    (1991)
  • J.U. Schlüter et al.

    A framework for coupling Reynolds-averaged with large-eddy simulations for gas turbine applications

    J. Fluids Eng.

    (2005)
  • S. Hahn et al.

    Coupled high-fidelity URANS simulation for helicopter applications

    Center for Turbulence Research, Annual Research Briefs

    (2006)
  • Y. Peet et al.

    Computational framework for coupling compressible and low Mach number codes

    AIAA J.

    (2008)
  • H.A. Schwarz

    Ueber einige Abbildungsaufgaben

    J. Reine Angew. Math.

    (1869)
  • E.A. Volkov

    The method of composite meshes for finite and infinite regions with piecewise smooth boundary

    Tr. Mat. Inst. Steklova

    (1968)
  • P.-L. Lions

    On the Schwarz alternating method. I

  • Y. Peet et al.

    Heat transfer LES simulations in application to wire-wrapped fuel pins

  • P.-L. Lions

    On the Schwarz alternating method. III: a variant for nonoverlapping subdomains

  • Y. Maday et al.

    Nonconforming Mortar Element Methods: Application to Spectral Discretizations

    (1988)
  • L.F. Pavarino et al.

    A polylogarithmic bound for an iterative substructuring method for spectral elements in three dimensions

    SIAM J. Numer. Anal.

    (1996)
  • J.C. Strikwerda et al.

    A domain decomposition method for incompressible viscous flow

    SIAM J. Sci. Comput.

    (1993)
  • Y. Zang et al.

    A composite multigrid method for calculating unsteady incompressible flows in geometrically complex domains

    Int. J. Numer. Methods Fluids

    (1995)
  • P. Bjorstad et al.

    Domain Decomposition: Parallel Multilevel Methods for Elliptic Partial Differential Equations

    (2004)
  • Cited by (39)

    • On the theoretical foundation of overset grid methods for hyperbolic problems II: Entropy bounded formulations for nonlinear conservation laws

      2022, Journal of Computational Physics
      Citation Excerpt :

      Overset grids have been used with all the major approximation methods, finite difference, finite element, finite volume and spectral element [20],[10],[8],[3],[12],[4],[13],[16].

    • Overset meshes for incompressible flows: On preserving accuracy of underlying discretizations

      2021, Journal of Computational Physics
      Citation Excerpt :

      This transfer of solution can happen via: 1) OSS [14,31,45,66], wherein interpolation equations for fringe points constitute constraint equations in the global system of equations. 2) OS [22,33,50,56,67,68], wherein the linear system for each mesh partition is solved independently and the solution at fringe points is updated in an iterative manner using interpolation equations. The iterations can be Gauss-Seidel-like (alternating Schwarz) or Jacobi-like (additive Schwarz).

    View all citing articles on Scopus
    View full text