High-order Spectral Element Method Research Articles

On the Origin of Görtler Vortices in Flow over a Multi-Element Airfoil

The flow characteristics of a 30P30N three-element high-lift airfoil at low Reynolds numbers O104 are examined through three-dimensional simulations using a high-order spectral element method. This study primarily investigates the flow structures of the slat cove and Görtler vortices formed on the upper surface of the main airfoil. Spanwise instability grows exponentially in the slat cove with a constant wavelength, corresponding to that of the subsequently formed Görtler vortices. Görtler number calculations show that curvature-induced centrifugal instability at the slat cusp leads to the subsequent formation of Görtler vortices. Proper orthogonal decomposition (POD) is used to analyze the development of flow structures in the slat cove in different time ranges. At early time, the flow in the slat cove is dominated by shear layers that evolve into spanwise perturbations. These perturbations further evolve into distinct bell-shaped structures close to the slat cusp and are advected to the upper surface of the main airfoil, leading to the formation of Görtler vortices.

A reduction-consistent phase field model for non-isothermal multiphase flows of N immiscible incompressible fluids

In this study, we developed a reduction-consistent phase field model for non-isothermal incompressible N-phase flows. The model is based on the high-order spectral element method. To account for thermocapillary effects, we reformulated the continuum surface force model specifically for N-phase flows. In order to enhance computational efficiency, time-independent coefficient matrices for all variables involved were reconstructed, and an unsymmetrized multifrontal LU factorization was employed to solve the linear algebraic equations, which is derived by discretization. To verify the model effectiveness in calculating surface tension and describing three-phase interfacial dynamics, we conducted several experiments, including the stationary two-droplet example, the three-phase droplet spreading, and the equilibrium morphology of double emulsion droplet. Through these experiments, both the reduction-consistency and robustness of our model were demonstrated. Moreover, we validated the proposed model's applicability to non-isothermal three-phase flows by investigating the thermocapillary migration of single droplet and two droplets of different phases. Notably, we explored the thermocapillary migration of three-phase double droplets, focusing particularly on how the encapsulation process affects overall thermocapillary motion. Our findings indicate that the interaction between the two vortices near the interfaces of the two droplets strongly influences their migration behavior.

Algebraically stable SDIRK methods with controllable numerical dissipation for first/second-order time-dependent problems

In this paper, a family of four-stage singly diagonally implicit Runge-Kutta methods are proposed to solve first-/second-order time-dependent problems, exhibiting the following numerical properties: fourth-order accuracy in time, unconditional stability, controllable numerical dissipation, and adaptive time step selection. The BN-stability condition is employed as a constraint to optimize parameters in the Butcher table, having significant benefits, and hence is recommended for nonlinear dynamics problems in contrast to existing methods. Numerical examples involving both first- and second-order linear/nonlinear dynamics problems validate the proposed method, and numerical results reveal that the proposed methods are free from the order reduction phenomenon when applied to nonlinear dynamics problems. The performance of adaptive time-stepping using the embedded scheme is further illustrated by the phase-field modeling problem. Additionally, the advantages and disadvantages of three-stage third-order accurate algebraically stable methods are discussed. The proposed high-order time integration can be readily integrated into high-order spatial discretization methods, such as the high-order spectral element method employed in this paper, to obtain high-order discretization in space and time dimensions.

Energy Exascale Computational Fluid Dynamics Simulations With the Spectral Element Method

Abstract Development and application of the open-source GPU-based fluid-thermal simulation code, NekRS, are described. Time advancement is based on an efficient kth-order accurate timesplit formulation coupled with scalable iterative solvers. Spatial discretization is based on the high-order spectral element method (SEM), which affords the use of fast, low-memory, matrix-free operator evaluation. Recent developments include support for nonconforming meshes using overset grids and for GPU-based Lagrangian particle tracking. Results of large-eddy simulations of atmospheric boundary layers for wind-energy applications as well as extensive nuclear energy applications are presented.

A Spectral/hp-Based Stabilized Solver with Emphasis on the Euler Equations

The solution of compressible flow equations is of interest with many aerospace engineering applications. Past literature has focused primarily on the solution of Computational Fluid Dynamics (CFD) problems with low-order finite element and finite volume methods. High-order methods are more the norm nowadays, in both a finite element and a finite volume setting. In this paper, inviscid compressible flow of an ideal gas is solved with high-order spectral/hp stabilized formulations using uniform high-order spectral element methods. The Euler equations are solved with high-order spectral element methods. Traditional definitions of stabilization parameters used in conjunction with traditional low-order bilinear Lagrange-based polynomials provide diffused results when applied to the high-order context. Thus, a revision of the definitions of the stabilization parameters was needed in a high-order spectral/hp framework. We introduce revised stabilization parameters, τsupg, with low-order finite element solutions. We also reexamine two standard definitions of the shock-capturing parameter, δ: the first is described with entropy variables, and the other is the YZβ parameter. We focus on applications with the above introduced stabilization parameters and analyze an array of problems in the high-speed flow regime. We demonstrate spectral convergence for the Kovasznay flow problem in both L1 and L2 norms. We numerically validate the revised definitions of the stabilization parameter with Sod’s shock and the oblique shock problems and compare the solutions with the exact solutions available in the literature. The high-order formulation is further extended to solve shock reflection and two-dimensional explosion problems. Following, we solve flow past a two-dimensional step at a Mach number of 3.0 and numerically validate the shock standoff distance with results obtained from NASA Overflow 2.2 code. Compressible flow computations with high-order spectral methods are found to perform satisfactorily for this supersonic inflow problem configuration. We extend the formulation to solve the implosion problem. Furthermore, we test the stabilization parameters on a complex flow configuration of AS-202 capsule analyzing the flight envelope. The proposed stabilization parameters have shown robustness, providing excellent results for both simple and complex geometries.

Multi-Level Solution Strategies for Implicit Time Discontinuous Galerkin Discretizations

The discontinuous Galerkin method and related flavors of high-order spectral element methods provide many well-known benefits for the spatial discretization of partial differential equations such as the Navier–Stokes equations. However, practical problems of engineering relevance such as large-eddy simulation of turbulent flows over complex geometries are computationally intractable by standard explicit time integration methods, necessitating the use of implicit methods. The efficient solution of the nonlinear algebraic systems arising from implicit time integration methods applied to DG discretizations of nonlinear PDEs is challenging; standard linearization methods result in very stiff block-sparse systems with prohibitive computational and memory requirements. This paper presents a low-memory, computationally efficient, implicit solution method that combines a framework of nonlinear polynomial multigrid, adaptive explicit Runge–Kutta smoothers, implicit Jacobian-free coarse level smoothers, nonlinear Krylov subspace acceleration and adaptive time stepping using feedback control of the nonlinear solver convergence rate.

A matrix–free high–order solver for the numerical solution of cardiac electrophysiology

We propose a matrix–free solver for the numerical solution of the cardiac electrophysiology model consisting of the monodomain nonlinear reaction–diffusion equation coupled with a system of ordinary differential equations for the ionic species. Our numerical approximation is based on the high–order Spectral Element Method (SEM) to achieve accurate numerical discretization while employing a much smaller number of Degrees of Freedom than first–order Finite Elements. We combine vectorization with sum–factorization, thus allowing for a very efficient use of high–order polynomials in a high performance computing framework. We validate the effectiveness of our matrix–free solver in a variety of applications and perform different electrophysiological simulations ranging from a simple slab of cardiac tissue to a realistic four–chamber heart geometry. We compare SEM to SEM with Numerical Integration (SEM–NI), showing that they provide comparable results in terms of accuracy and efficiency. In both cases, increasing the local polynomial degree p leads to better numerical results and smaller computational times than reducing the mesh size h. We also implement a matrix–free Geometric Multigrid preconditioner that results in a comparable number of linear solver iterations with respect to a state–of–the–art matrix–based Algebraic Multigrid preconditioner. As a matter of fact, the matrix–free solver proposed here yields up to 45× speed–up with respect to a conventional matrix–based solver.

Interference and ground effects on flow past two inclined flat plates in tandem arrangement

Direct Numerical Simulation (DNS) was conducted on the flow past two tandem inclined flat plates at a Reynolds number of 150. With an inclined angle θ = 25°, the high-order spectral-element method is employed for higher accuracy. Fluid force, power spectral density analysis, flow structures, and instability have been investigated focusing on both interference and ground effects. Relative magnitudes of gap spacing “X” and plate height “G” over plate length “L” are investigated with ranges: X/L = 1–6 at an interval of 1, G/L = 0.2 to 1.6 at an interval of 0.2. It was found that two asymmetrical vortices formed behind the upstream plate (UP), and then directly reattached to the downstream plate (DP), where larger ratios of X/L and G/L enhanced the vortex shedding behind the plates and magnify the fluid forces. A small gap spacing X/L was found not sufficient to develop the gap flow, where tandem plates were considered as a single body of vortex shedding. With G/L lower than 0.5, the vorticity generated at the trailing edge was gradually neutralized by the boundary shear layer which contributes to lower fluid force. Moreover, the instability of tandem plates was governed by the gap spacing ratio.

Large Eddy Simulation of Gasoline Sprays in a Lagrangian–Eulerian Framework Using the High-Order Spectral Element Method

Abstract Predicting the spray evolution using simulations requires accurate modeling of the turbulent gas-phase flow field. In this study, the high-order spectral-element method (SEM), implemented in the code Nek5000, was used to provide highly resolved solutions to the turbulent flow field. Spray modeling capabilities were implemented into the Nek5000 code. The spray is modeled in a Lagrangian–Eulerian (LE) framework, where the liquid is represented by discrete parcels of droplets. The method for coupling liquid and gas in the context of SEM is described, which allows for very fine meshes to be used without affecting the stability of the solution. Large-eddy simulations (LES) of the eight-hole ECN Spray G gasoline injector were conducted. Numerical results are compared against experimental data for liquid penetration, droplet size and gas velocity. The morphology of the multiplume spray is compared against experimental data. The effect of different spray injection inputs is analyzed. It was found that using a plume direction of 33 deg and an injection cone angle of 30 deg produced the best results overall. This work shows the applicability of SEM for spray modeling applications, where use of a high-order flow solver can help us understand the multiplume spray aerodynamics and how it leads to plume collapse under certain conditions. Results also highlight the need for tuning spray input parameters in the LE framework, even when high-fidelity gas flow solutions are possible.

High-order optimal mode decomposition analysis of the ground effect on flow past two tandem inclined plates

Two-dimensional flow past two tandem near-ground plates with inclination angles of 25° at the Reynolds number of 150 is numerically simulated via the high-order spectral element method. Plate-to-ground gap is varied from G = 0.2L to 1.6L with intervals of 0.2L at two representative inter-plate spacings (i.e., X = 2.5L and 6L). The ground effect on the fluid force, power spectral density, asymmetric gap flow, and wake structure of plates is systematically evaluated. Then, the high-order optimal mode decomposition (HOOMD) method is proposed to synchronously analyze the velocity and pressure fields. The results show that the fluid force and flow structure are closely dependent on G. The presence of the ground inhibits vortex shedding when G < 0.6L; as the gap increases from 0.6 L to 1.4 L, the fluctuating forces are continuously enhanced until the ground effect basically disappears at G > 1.4L. The ground effect exacerbates the asymmetry of the vortex structure near the upper and lower parts of the inclined plates, consequently changing the fluid force. The downstream plate is more sensitive to the ground effect because of impingement from the upward-biased jet flow generated in the narrow gap between the upstream plate and ground. The HOOMD method well captures the spatial morphology and temporal evolution features of different dominant modes at the transition or vortex shedding flow regime. Mode analysis affords a correspondence between the coherent vortex structure and fluid force of plates. Furthermore, the ground effect can simultaneously change the global mode energy and local pressure mode shape, subsequently influencing the fluid force. However, the global mode energy plays the determinant role in the variation of the fluid force of plates with the plate-to-ground distance herein.

Large-eddy simulation of non-vaporizing sprays using the spectral-element method

Predictive simulations of high-pressure sprays require accurate representation of the turbulent gaseous flow field generated by liquid jet. Typically, the accuracy that can be obtained with low-order numerical methods (e.g. finite volume, finite element) is limited by stability issues in fine grids and the order of convergence of the method. In this work, we resolve the turbulent flow field in an Eulerian manner using the high-order spectral element method, coupled with a Lagrangian parcels approach to model the atomizing liquid jet. Large-eddy simulations of single-hole sprays under non-evaporative conditions were conducted and compared against experimental data from Margot et al. (2008) and Spray A data from the Engine Combustion Network. The sensitivity of liquid penetration and droplet sizes to different breakup model parameters was studied. The effect of different numerical parameters, such as polynomial order of the solution (grid resolution), on liquid penetration was also analyzed. The method achieved grid-independent results using p-refinement, achieving finer resolution (by a factor of ×1.7−×3.5) in the gas-phase solution than in state-of-the-art simulations using the finite-volume method. Results showed good agreement with experimental data, demonstrating the ability of the current method to accurately capture liquid penetration and the shape of the spray.

Scalability of Nek5000 on High-Performance Computing Clusters Toward Direct Numerical Simulation of Molten Pool Convection

In a postulated severe accident, a molten pool with decay heat can form in the lower head of a reactor pressure vessel, threatening the vessel’s structural integrity. Natural convection in molten pools with extremely high Rayleigh (Ra) number is not yet fully understood as accurate simulation of the intense turbulence remains an outstanding challenge. Various models have been implemented in many studies, such as RANS (Reynolds-averaged Navier–Stokes), LES (large-eddy simulation), and DNS (direct numerical simulation). DNS can provide the most accurate results but at the expense of large computational resources. As the significant development of the HPC (high-performance computing) technology emerges, DNS becomes a more feasible method in molten pool simulations. Nek5000 is an open-source code for the simulation of incompressible flows, which is based on a high-order SEM (spectral element method) discretization strategy. Nek5000 has been performed on many supercomputing clusters, and the parallel performance of benchmarks can be useful for the estimation of computation budgets. In this work, we conducted scalability tests of Nek5000 on four different HPC clusters, namely, JUWELS (Atos Bullsquana X1000), Hawk (HPE Apollo 9000), ARCHER2 (HPE Cray EX), and Beskow (Cray XC40). The reference case is a DNS of molten pool convection in a hemispherical configuration with Ra = 1011, where the computational domain consisted of 391 million grid points. The objectives are (i) to determine if there is strong scalability of Nek5000 for the specific problem on the currently available systems and (ii) to explore the feasibility of obtaining DNS data for much higher Ra. We found super-linear speed-up up to 65536 MPI-rank on Hawk and ARCHER2 systems and around 8000 MPI-rank on JUWELS and Beskow systems. We achieved the best performance with the Hawk system with reasonably good results up to 131072 MPI-rank, which is attributed to the hypercube technique on its interconnection. Given the current HPC technology, it is feasible to obtain DNS data for Ra = 1012, but for cases higher than this, significant improvement in hardware and software HPC technology is necessary.

Topology optimization of unsteady flows using the spectral element method

We investigate the applicability of a high-order Spectral Element Method (SEM) to density based topology optimization of unsteady flows in two dimensions. Direct Numerical Simulations (DNS) are conducted relying on Brinkman penalization to describe the presence of solid within the domain. The optimization procedure uses the adjoint-variable method to compute gradients and a checkpointing strategy to reduce storage requirements. A nonlinear filtering strategy is used to both enforce a minimum length scale and to provide smoothing across the fluid–solid interface, preventing Gibbs oscillations. This method has been successfully applied to the design of a channel bend and an oscillating pump, and demonstrates good agreement with body fitted meshes. The precise design of the pump is shown to depend on the initial material distribution. However, the underlying topology and pumping mechanism is the same. The effect of a minimum length scale has been studied, revealing it to be a necessary regularization constraint for the oscillating pump to produce meaningful designs. The combination of SEM and density based optimization offer some unique challenges which are addressed and discussed, namely a lack of explicit boundary tracking exacerbated by the interface smoothing. Nevertheless, SEM can achieve equivalent levels of precision to traditional finite element methods, while requiring fewer degrees of freedom. Hence, the use of SEM addresses the two major bottlenecks associated with optimizing unsteady flows: computation cost and data storage.

2D acoustic wavefield simulation by improved stiffness and mass matrices of bilinear square element

Abstract In the frequency domain, we optimize the entries of the element stiffness matrix computed on a square element with the purpose of approximating the 2D acoustic wave equation with the second spatial order. The optimized matrices are computed from the minimization of the normalized phase velocity as a function of propagation angle and the number of grid points per wavelength, leading to a remarkable reduction in the numerical error. The optimized number of 4.1 sample points in a wavelength, with a maximum error in velocities <0.3%, offering the simulation of large-scale realistic models. The superior performance of the optimized matrices to those from the lumped, consistent and eclectic matrices, is evident from the numerical examples presented, with consequent reduction not only in the numerical dispersion/anisotropy but also in an improved resolution. It is noteworthy that optimized matrices give modeling results with better accuracy than models from the high-order spectral element method. Our methodology is easily extendable to the accurate discretizing of the general 3D elastic wave equation that optimizes the bandwidth of the impedance matrix, availing computational resources for accurately modeling the compressional and shear wave velocities, density, as well as multi-parameter inversion.

High-Order Spectral/-Based Solver for Compressible Navier–Stokes Equations

High-order spectral element methods are becoming more established for solving problems in computational fluid dynamics. In this paper, we introduce new stabilization parameters for the solution of compressible Navier–Stokes equations in conjunction with high-order spectral element methods. We previously defined a new set of stabilization parameters for Euler equations in a high-order spectral/ framework. In this paper, we examine the definition of these stabilization parameters in the context of solutions of full Navier–Stokes equations, and we report good agreement with previously published results in the literature. We establish spectral convergence of errors for Kovasznay flow. We then validate the definition of the stabilization parameter for Couette flow, flat plate, and compression corner problems. In a later example, we solve the flow past a cylinder and benchmark results with previously published results obtained with low-order methods and other approaches in compressible flow computations. We further solve hypersonic flow past a wedge and supersonic flow past a Flight Investigation of Re-Entry program (FIRE) entry capsule. We explore the definition of the shock-capturing parameter based on the parameter. The introduced definitions are found to provide excellent results for the full Navier–Stokes equations.

Numerical investigations on model order reduction to SEM based on POD-DEIM to linear/nonlinear heat transfer problems

The motivation in this article is to foster and conduct an in-depth investigation on reduced-order modeling (ROM) techniques via proper orthogonal decomposition (POD) and discrete empirical interpolation method (DEIM) in conjunction with high-order Legendre spectral element method (SEM) applied for time-dependent linear and nonlinear heat conduction problems. This approach significantly relieves the CPU burden, yet preserves the numerical accuracy of high-order schemes; this is especially pronounced for linear transient problems. For nonlinear cases, we use DEIM to save CPU time for nonlinear counterparts. For this purpose, we first obtain training data by solving the system of full-order models (FOM) using snapshot techniques to construct the POD basis. The well-known generalized single step single solve (GSSSS-1) framework is next employed for the first-order system for the time discretization. Numerical results for both linear and nonlinear heat conduction problems such as in a nuclear reactor, etc., show excellent agreement between FOM solutions and ROM solutions, and the CPU advantages via POD ROM techniques are also prominent for high-order Legendre SEM.

A hybrid spectral/finite element method for accurate and efficient modelling of crack-induced contact acoustic nonlinearity

We propose a 2D hybrid spectral/finite element scheme for numerically resolving crack-induced contact acoustic nonlinearity in solid structures. While the high-order spectral element method (SEM) is much more efficient than the classical low-order finite element method (FEM), the fact that spectral elements are relatively large renders the SEM ineffective in discretizing microscopic cracks. The classical FEM, on the other hand, is good at modelling complex geometries. However, the computation of large-scale structures by the classical FEM can be extremely expensive. In this work, the coupling between high-order spectral elements and low-order finite elements is formulated by the Lagrange multipliers method. The nonlinear contact between crack surfaces is modelled using a penalty method. The aim of the proposed hybrid method is to simultaneously reduce the degrees of freedom (DOFs) in systems and ensure high numerical accuracy. A comprehensive mesh convergence study has been conducted to demonstrate the high convergence rate of the proposed method in modelling ultrasonic wave propagation, and its capability to converge for weak crack-induced contact acoustic nonlinearity. Through a series of numerical experiments on generation and propagation of nonlinear ultrasonic wave modes, the proposed method has been shown to possess a high accuracy and geometric flexibility, and require significantly less computational resources than the classical FEM. The proposed hybrid method will provide an efficient, accurate and robust numerical approach for studying crack-induced contact acoustic nonlinearity which nowadays is being widely exploited by non-destructive testing applications.

Multiple solutions and hysteresis in the flows driven by surface with antisymmetric velocity profile* *Project supported by the National Natural Science Foundation of China (Grant Nos. 11902043 and 11772065), and the Science Challenge Project (Grant No. TZ2016001).

Multiple steady solutions and hysteresis phenomenon in the square cavity flows driven by the surface with antisymmetric velocity profile are investigated by numerical simulation and bifurcation analysis. A high order spectral element method with the matrix-free pseudo-arclength technique is used for the steady-state solution and numerical continuation. The complex flow patterns beyond the symmetry-breaking at Re ≃ 320 are presented by a bifurcation diagram for Re < 2500. The results of stable symmetric and asymmetric solutions are consistent with those reported in literature, and a new unstable asymmetric branch is obtained besides the stable branches. A novel hysteresis phenomenon is observed in the range of 2208 < Re < 2262, where two pairs of stable and two pairs of unstable asymmetric steady solutions beyond the stable symmetric state coexist. The vortices near the sidewall appear when the Reynolds number increases, which correspond to the bifurcation of topology structure, but not the bifurcation of Navier–Stokes equations. The hysteresis is proposed to be the result of the combined mechanisms of the competition and coalescence of secondary vortices.

Slat cove dynamics of low Reynolds number flow past a 30P30N high lift configuration

A three-dimensional computational fluid dynamics analysis of low Reynolds number [O(104)] flow over a 30P30N three-element high lift wing is carried out using a high-order spectral element method. In this article, we study the flow in the slat cove region and the slat wake/shear layer interaction. Vortical structures, identified in the computations, are very similar to those visualized in recent experiments. For Reynolds numbers below a critical interval (found in recent experiments), Görtler vortices are observed in the slat wake, while for Reynolds numbers above the critical interval, a roll-up is observed in the slat cove and both streamwise and spanwise vortices form in the slat wake. Prior to the formation of Görtler and roll-up vortices, three-dimensional tongue- or rib-like vortex shapes, similar to those found in the wake of bluff bodies, are observed in the slat cove and promote transition to three-dimensional flow. Above the critical interval, streaks and spanwise vortices are observed to dominate the slat wake and lead to the formation of hairpin vortices which contribute to the transition to turbulence. Integral flow parameters such as lift, drag, and pressure coefficients are analyzed in the range of Reynolds numbers studied.

Onset of temporal dynamics within a low reynolds-number laminar fluidic oscillator

We investigate the onset of temporal dynamics of a fluidic oscillator (FO) in the laminar regime at very low Reynolds number (Re=Uhi/ν, where U and hi are the velocity and channel height at inlet, and ν is the kinematic viscosity of the fluid). Both two- and spanwise-periodic-three-dimensional simulations are performed using high-order spectral element methods to characterise the flow inside the FO cavity and the outcoming jet. While the flow remains steady and symmetric for sufficiently low Re, the two-dimensional FO undergoes a symmetry-breaking Hopf bifurcation at ReH2=75.7 that results in a space-time symmetric periodic solution. The remnant space-time symmetry is later broken in a pitchfork bifurcation of limit cycles at ReP2=294.2. Taking three-dimensionality into consideration results in an early spanwise-invariance disruption at about Re3≃68 of wavelength 13hi that retards the onset of the oscillation to ReH3=82.1. This effect is little representative of actual FO implementations, which are typically much narrower than the wavelength of the spanwise instability observed at these low Re. Contrary to what happens with FOs in the turbulent oscillatory regime, the oscillation frequency of the output jet in the laminar regime is found to decrease with Re. The mechanism that drives the oscillation, however, remains the same: as the high speed flow inside the cavity attaches to one of the internal walls through the Coandă effect, the feedback channels divert part of the momentum back, which pushes the incoming flow against the opposite wall. The alternate deviation of the flow inside the cavity to one or the other side results in the periodic flipping of the output jet. The pressure contribution to net momentum at the end of the feedback channels is found to be double that of the advective momentum flux at the early time-periodic regime but both become comparable as Re is increased. The jet sweeping angle amplitude is more pronounced for the two-dimensional FO as compared to three-dimensional at a fixed given Re, the Coandă effect being only partially fulfilled in the latter case. In both cases the sweep amplitude increases with Re. The instability of the output jet, which becomes slightly chaotic already at very low Re, is responsible for triggering the cavity instability that drives the oscillation slightly earlier than it would, should outside noise be suppressed altogether.