Courant-Friedrichs-Lewy Numbers Research Articles

A stable staggered-grid finite-difference scheme for acoustic modeling beyond conventional stability limit

Staggered-grid finite-difference (SGFD) schemes have been widely used in acoustic wave modeling for geophysical problems. Many improved methods are proposed to enhance the accuracy of numerical modeling. However, these methods are inevitably limited by the maximum Courant-Friedrichs-Lewy (CFL) numbers, making them unstable when modeling with large time sampling intervals or small grid spacings. To solve this problem, we extend a stable SGFD scheme by controlling SGFD dispersion relations and maximizing the maximum CFL numbers. First, to improve modeling stability, we minimize the error between the FD dispersion relation and the exact relation in the given wave-number region, and make the FD dispersion approach a given function outside the given wave-number area, thus breaking the conventional limits of the maximum CFL number. Second, to obtain high modeling accuracy, we use the SGFD scheme based on the Remez algorithm to compute the FD coefficients. In addition, the hybrid absorbing boundary condition is adopted to suppress boundary reflections and we find a suitable weighting coefficient for the proposed scheme. Theoretical derivation and numerical modeling demonstrate that the proposed scheme can maintain high accuracy in the modeling process and the value of the maximum CFL number of the proposed scheme can exceed that of the conventional SGFD scheme when adopting a small maximum effective wavenumber, indicating that the proposed scheme improves stability during the modeling.

CABARET on rotating meshes

The Compact Accurately Boundary-Adjusting high-REsolution Technique (CABARET) method for 3D gas dynamics equations is generalised to rotating hexahedral mesh zones using the sliding interface approach, which extends applications of the original algorithm to a wide range of unsteady rotor flow problems. Important properties of the suggested extension of the CABARET scheme include the preservation of good linear dispersion and dissipation properties of the original CABARET method on fixed grids. To strictly preserve flux conservation at the sliding interface, the so-called supermesh approach is implemented that uses the projection of the contact faces on both sides of the sliding interface onto an intermediate virtual surface. The implementation is based on applying the sub-cell and sub-face solution reconstruction on the sliding interface of hexahedral meshes rather than deforming unstructured meshes. The sub-faces are defined by constructing projection overlap polygons using the Sutherland-Hodgman clipping algorithm. A local characteristic flux splitting is applied at each overlap polygon face, while updating the local stencil connectivity table at each time step. To further improve the dispersion properties of the CABARET grid at small Courant–Friedrichs–Lewy numbers, the dispersion-improved version of the CABARET method previously developed on rectangular Cartesian meshes is extended to non-uniform curvilinear meshes typical of open-rotor applications. Numerical examples are provided for the analytical problem of normal acoustic wave propagation through a rotating cylindrical grid zone embedded in an outer stationary grid. The error decay rates are calculated and compared with the theoretical order of CABARET convergence. In addition, the same algorithm accelerated on Graphics Processing Units (GPUs) is applied for simulations the flow around a generic two-bladed rotor in hover in the Wall Model Large Eddy Simulation framework. For the latter case, flow solutions are analysed, and their numerical grid sensitivity is discussed.

Implicit block data-parallel relaxation scheme of Navier–Stokes equations using graphics processing units

The granularity of computational fluid dynamics (CFD) generally refers to the point granularity parallelization as a unit of the grid when graphics processing units (GPUs) are utilized as the computing carrier. In commonly deployed implicit time advancement schemes, the parallel dimensionality must be reduced, resulting in the time advancement procedure becoming the only highly time-consuming step in the whole CFD computing procedures. In this paper, a block data-parallel lower-upper relaxation (BDPLUR) scheme based on Jacobi iteration and Roe's flux scheme is proposed and then implemented on a GPU. Numerical experiments are carried out and show that the convergence speed of the BDPLUR scheme, especially when implemented on a GPU, is approximately ten times higher than that of the original data-parallel lower-upper relaxation scheme and more than 100 times higher than that of the lower-upper symmetric Gauss–Seidel scheme. Moreover, the influence of different Courant–Friedrichs–Lewy numbers on the convergence time is discussed, and different viscous matrices are compared. Standard cases are adopted to verify the effectiveness of the BDPLUR scheme.

Removing the stability limit of the time-space domain explicit finite-difference schemes for acoustic modeling with stability condition-based spatial operators

The time step and grid spacing in explicit finite-difference (FD) modeling are constrained by the Courant-Friedrichs-Lewy (CFL) condition. Recently, it has been found that the spatial FD coefficients may be designed through simultaneously minimizing the spatial dispersion error and maximizing the CFL number. This allows one to stably use a larger time step or a smaller grid spacing than usually possible. However, when using such a method, only second-order temporal accuracy is achieved. To address this issue, I propose a new method to determine the spatial FD coefficients, which simultaneously satisfy the stability condition of the whole wavenumber range and the time-space domain dispersion relation of a given wavenumber range. Therefore, stable modeling can be performed with high-order spatial and temporal accuracy. The coefficients can adapt to the variation of velocity in heterogeneous models. In addition, based on the hybrid absorbing boundary condition, I develop a strategy to stably and effectively suppress artificial reflections from the model boundaries for large CFL numbers. Stability analysis, accuracy analysis, and numerical modeling demonstrate the accuracy and effectiveness of my proposed method.

Acoustic finite-difference modeling beyond conventional Courant-Friedrichs-Lewy stability limit: Approach based on variable-length temporal and spatial operators

Conventional finite-difference (FD) methods cannot model acoustic wave propagation beyond Courant-Friedrichs-Lewy (CFL) numbers 0.707 and 0.577 for two-dimensional (2D) and three-dimensional (3D) equal spacing cases, respectively, thereby limiting time step selection. Based on the definition of temporal and spatial FD operators, we propose a variable-length temporal and spatial operator strategy to model wave propagation beyond those CFL numbers while preserving accuracy. First, to simulate wave propagation beyond the conventional CFL stability limit, the lengths of the temporal operators are modified to exceed the lengths of the spatial operators for high-velocity zones. Second, to preserve the modeling accuracy, the velocity-dependent lengths of the temporal and spatial operators are adaptively varied. The maximum CFL numbers for the proposed method can reach 1.25 and 1.0 in high velocity contrast 2D and 3D simulation examples, respectively. We demonstrate the effectiveness of our method by modeling wave propagation in simple and complex media.

An implicit kinetic inviscid flux for predicting continuum flows in all speed regimes

In this study, the kinetic inviscid flux (KIF) is improved and coupled with an implicit strategy. The KIF is a recently proposed numerical method, which is a dynamic combination of the kinetic flux vector splitting (KFVS) method and the totally thermalized transport (TTT) method. The inherent microscopic mechanism of the KFVS makes the KIF good at solving shock waves and avoiding the numerical shock instability phenomenon. When developing the implicit KIF, it is noticed that, in boundary layers, the KFVS part of the KIF not only reduces the accuracy but also seriously reduces the Courant–Friedrichs–Lewy (CFL) number. As a result, a new weight is proposed in this paper to combine the KFVS method with the TTT method properly. Besides admitting the use of larger CFL numbers, this new weight also contributes to more accurate numerical results like pressure, friction coefficient, and heat flux when solving shock waves, boundary layers, and complex supersonic/hypersonic flows. To examine the validity, accuracy, and efficiency of the proposed method, six numerical test cases covering the whole speed regime are conducted, including the hypersonic viscous flow past a cylinder, the hypersonic double-cone flow, the hypersonic double-ellipsoid flow, the laminar shock-boundary layer interaction, the supersonic flow around a ramp segment and the subsonic lid-driven cavity flow.

A ghost-cell immersed boundary method for unified simulations of flow over finite- and zero-thickness moving bodies at large CFL numbers

A ghost-cell immersed boundary method for unified simulations of flow over finite- and zero-thickness moving bodies at large Courant-Friedrichs-Lewy (CFL) numbers is presented. In order to handle such bodies in a unified manner, algorithms for interface construction and cell demarcation are proposed. The main challenge in treating zero-thickness bodies is to maintain sharpness and accuracy even at large CFL numbers with diminished spurious force oscillations. Thus, the effect of large CFL numbers on the solution accuracy of fluid-structure interaction (FSI) problems involving zero-thickness bodies is investigated and necessary treatments to preserve solution accuracy even at large CFL numbers are suggested. The present study suggests two treatments which are important in preserving the accuracy and stability of the solution: backward time integration for computational cells called ‘swept-cells’ and pressure boundary condition with mass conservation. Composite implicit time integration for the dynamic equation of a thin elastic structure is employed for a stable simulation of FSI at large CFL numbers. By using large time step sizes, the present method not only enhances computational efficiency, but also suppresses spurious force oscillations while maintaining the sharpness of an infinitesimally thin body. The efficacy and accuracy of the present method are examined through numerical examples.

Improved standard thermal lattice Boltzmann model with hybrid recursive regularization for compressible laminar and turbulent flows

Based on recent work by Guo et al. [“An efficient lattice Boltzmann method for compressible aerodynamics on D3Q19 lattice,” J. Comput. Phys. 418, 109570 (2020)], an improved thermal hybrid recursive regularized lattice Boltzmann model (iHRR-ρ) on a regular lattice is developed for two- and three-dimensional compressible laminar and turbulent flows. To enhance the numerical stability in a broad range of Courant–Friedrichs–Lewy numbers and in under-resolved simulations, a new equilibrium density distribution function is proposed to enlarge its positivity region in the Mach–temperature plane. An embedded hybridizing procedure is introduced in the quasi-symmetry correction terms, which allow for a decoupled treatment of unphysical modes and physical under-resolved turbulent scales on coarse grids. To handle compressible turbulent flows, the under-resolved scales are modeled using the original hybrid recursive regularized collision model given by Jacob et al. [“A new hybrid recursive regularised Bhatnagar–Gross–Krook collision model for Lattice Boltzmann method-based large eddy simulation,” J. Turbul. 19, 1051–1076 (2018)] equipped with Vreman’s subgrid model for the large-eddy simulation. The validity and accuracy of the present method for laminar and turbulent compressible flows are assessed by considering six test cases: (I) viscous shock wave internal structure, (II) isentropic vortex convection in a supersonic regime, (III) non-isothermal acoustic pulse, (IV) vortex–shock wave interaction, (V) supersonic flow over NACA airfoil at Re = 10 000 and Ma = 1.5, and (VI) compressible Taylor–Green vortex at Ma = 0.29.

Time adaptive conservative finite volume method

Finite volume methods for time dependent problems typically employ time integration schemes with constant time step sizes. The latter, in order to ensure stability and time accurate solutions, have to be chosen small enough, such that the Courant-Friedrichs-Lewy (CFL) number stays below a critical value everywhere. In many cases, however, this time step size limitation leads to huge numbers of very small time steps, which renders simulations very expensive. Motivated by this drawback of conventional time stepping schemes, various sub-time stepping algorithms have been proposed. However, most of them are inherently asynchronous, require small local CFL numbers or are not strictly conservative. In this paper a new adaptive time integration scheme for finite volume methods, which is conservative, of high spatial and temporal order, robust and easy to implement is presented. It relies on local sub-time steps which are fractions of a global time step by powers of two, i.e., the grid cells proceed asynchronously in the order of their termination time, but since they all synchronize at the end of each global time step, it is possible to guarantee continuity of the mean fluxes and thus strict conservation at the global time step resolution. Numerical experiments with 1D and 2D test cases demonstrate that the adaptive conservative time integration (ACTI) scheme can achieve extreme speed-up factors over conventional time integration, while still maintaining high spatial and temporal accuracy.

One-dimensional Lagrangian implicit hydrodynamic algorithm for Inertial Confinement Fusion applications

A new one-dimensional hydrodynamic algorithm, specifically developed for Inertial Confinement Fusion (ICF) applications, is presented. The scheme uses a fully conservative Lagrangian formulation in planar, cylindrical, and spherically symmetric geometries, and supports arbitrary equations of state with separate ion and electron components. Fluid equations are discretized on a staggered grid and stabilized by means of an artificial viscosity formulation. The space discretized equations are advanced in time using an implicit algorithm. The method includes several numerical parameters that can be adjusted locally. In regions with low Courant–Friedrichs–Lewy (CFL) number, where stability is not an issue, they can be adjusted to optimize the accuracy. In typical problems, the truncation error can be reduced by a factor between 2 to 10 in comparison with conventional explicit algorithms. On the other hand, in regions with high CFL numbers, the parameters can be set to guarantee unconditional stability. The method can be integrated into complex ICF codes. This is demonstrated through several examples covering a wide range of situations: from thermonuclear ignition physics, where alpha particles are managed as an additional species, to low intensity laser–matter interaction, where liquid–vapor phase transitions occur.

Application of the characteristics-based sectional method to spatially varying aerosol formation and transport

The characteristics-based ssolution. It is easy to verify thatectional method (CBSM) offers an Eulerian description of an internally mixed aerosol. It was shown to be robust and capable of exact preservation of lower order moments, allowing for highly skewed sectional droplet size distributions. In this paper we apply CBSM to a spatially varying flow, by incorporating the fractional step method. In this way an accurate time integration of the spatial terms in the transport equations for the velocity, mass fractions and sectional droplet concentrations is achieved. Integrating CBSM into the compressible PISO (Pressure-Implicit with Splitting of Operators) algorithm allows for phase change and corresponding changes in pressure. We apply CBSM to a lid-driven cavity flow. First, the steady state isothermal flow solution is validated against published data. Next, by releasing a saturated vapor into the cavity while cooling the walls, we simulate the formation of aerosol. The accuracy of the solution is studied, as well as the performance of the CBSM scheme in the spatially varying context. The solution of the velocity is shown to be accurate, even at CFL (Courant–Friedrichs–Lewy) numbers of unity. The feasibility of the developed method is demonstrated in a 3D complex geometry studying the aerosol generation via nucleation of hot vapors cooled by a dilution stream of cold air in a double-mixing tee system. The sectional approach delivers detailed information about the aerosol formation and size distribution of the droplets in the domain.

A high order characteristic discontinuous Galerkin scheme for advection on unstructured meshes

A new characteristic discontinuous Galerkin (CDG) advection scheme is presented. In contrast to standard discontinuous Galerkin schemes, the test functions themselves follow characteristics in order to ensure conservation and the edges of each element are also traced backwards along characteristics in order to create a swept region, which is integrated in order to determine the mass flux across the edge. Both the accuracy and performance of the scheme are greatly improved by the use of large Courant–Friedrichs–Lewy numbers for a shear flow test case and the scheme is shown to scale sublinearly with the number of tracers being advected, outperforming a standard flux corrected transport scheme for 10 or more tracers with a linear basis. Moreover the CDG scheme may be run to arbitrarily high order spatial accuracy and on unstructured grids, and is shown to give the correct order of error convergence for piecewise linear and quadratic bases on regular quadrilateral and hexahedral planar grids. Using a modal Taylor series basis, the scheme may be made monotone while preserving conservation with the use of a standard slope limiter, although this reduces the formal accuracy of the scheme to first order. The second order scheme is roughly as accurate as the incremental remap scheme with nonlocal gradient reconstruction at half the horizontal resolution. The scheme is being developed for implementation within the Model for Prediction Across Scales (MPAS) Ocean model, an unstructured grid finite volume ocean model.

Decomposition model for naval hydrodynamic applications, Part I: Computational method

The present paper and its companion present a solution and domain decomposition model for general two-phase, incompressible and turbulent flows encountered in naval hydrodynamics. The mathematical and numerical model are derived within the framework of polyhedral finite volume method. Interface capturing is obtained with implicitly redistanced level set method derived from the phase field equation. This approach removes the need to redistance the level set field, thus saving CPU time and increasing numerical stability. A modified, Spectral Wave Explicit Navier–Stokes Equations method is introduced and used for wave modelling, where the solution is decomposed into incident and perturbation fields. The incident field is readily available from potential flow theories, and only the perturbation component is solved within the non-linear equation set of the free surface flow model. The domain is decomposed with implicit relaxation zones used to prevent wave reflection by forcing the perturbation fields to vanish in the far-field. A second-order, collocated finite volume method is used to discretise the equations. All equations are solved implicitly, which enables the use of higher Courant–Friedrichs–Lewy numbers compared to explicit scheme. The algorithm is implemented in foam-extend-3.1, a community driven fork of the OpenFOAM CFD software. Verification and validation of the model is presented in the accompanying paper.

Stability property of numerical Cherenkov radiation and its application to relativistic shock simulations

Abstract We studied the stability property of numerical Cherenkov radiation in relativistic plasma flows employing particle-in-cell simulations. Using the implicit finite-difference time-domain method to solve the Maxwell equations, we found that nonphysical instability was greatly inhibited with a Courant–Friedrichs–Lewy (CFL) number of 1.0. The present result contrasts with recently reported results (Vay et al. 2011, J. Comp. Phys., 230, 5908; Godfrey & Vay 2013, J. Comp. Phys., 248, 33; Xu et al. 2013, Comput. Phys. Commun., 184, 2503) in which magical CFL numbers in the range 0.5–0.7 were obtained with explicit field solvers. In addition, we found employing higher-order shape functions and an optimal implicitness factor further suppressed long-wavelength modes of the instability. The findings allowed the examination of the long-term evolution of a relativistic collisionless shock without the generation of nonphysical wave excitations in the upstream. This achievement will allow us to investigate particle accelerations in relativistic shocks associated with, for example, gamma-ray bursts.

IMPROVEMENT OF COMPUTATIONAL PERFORMANCE OF IMPLICIT FINITE DIFFERENCE TIME DOMAIN METHOD

Different solution techniques, computational aspects and the ways to improve the performance of 3D frequency dependent Crank Nicolson finite difference time domain (FD-CN-FDTD) method are extensively studied here. FD-CN-FDTD is an implicit unconditionally stable method allowing time discretization beyond the Courant-Friedrichs-Lewy (CFL) limit. For the solution of the method both direct and iterative solver approaches have been studied in detail in terms of computational time, memory requirements and the number of iteration requirements for convergence with different CFL numbers (CFLN ). It is found that at higher CFLN more iterations are required to converge resulting in increased number of matrix-vector multiplications. Since matrix-vector multiplications account for the most significant part of the computations their efficient implementation has been studied in order to improve the overall efficiency. Also the scheme has been parallelized in shared memory architecture using OpenMP and the resulted improvement of performance at different CFLN is presented. It is found that better speed-up due to parallelization always comes at higher CFLN implying that the use of FD-CN-FDTD method is more appropriate while parallelized.

Numerical and physical requirements to simulation of gas release and dispersion in an enclosure with one vent

Numerical and physical requirements to simulations of sub-sonic release and dispersion of light gas in an enclosure with one vent are described and discussed. Six validation experiments performed at CEA in a fuel cell-like enclosure of sizes H × W × L = 126 × 93 × 93 cm with one vent, either W × H = 90 × 18 cm (vent A) or 18 × 18 cm (B) or 1 cm in diameter (C), with a vertical upward helium release from a pipe of internal diameter either 5 mm or 20 mm located 21 cm above the floor centre, were used in a parametric study comprising 17 numerical simulations. Three CFD models were applied, i.e. laminar, standard k-ɛ, and dynamic LES Smagorinsky–Lilly, to clarify a range of their applicability and performance. The LES model consistently demonstrated the best performance in reproduction of measured concentrations throughout the whole range of experimental conditions, including laminar, transitional and turbulent releases even with large CFL numbers. The laminar and the standard k-ɛ models were under performing in the reproduction of turbulent and laminar releases respectively, as expected, as well as in simulation of transitional flows. The laminar model demonstrated high sensitivity to the CFL (Courant–Friedrichs–Lewy) number even below the best practices limit of 40. Three different computational domains and grids were used in order to clarify the influence of mesh quality on the capability of simulations to reproduce the experimental data. It is concluded that physically substantiated choice of CFD model, the control of the CFL number (and released gas mass balance where appropriate), and the mesh quality can have a strong effect on the capability of simulations to reproduce experiments and, in general, on the reliability of CFD tools for application in hydrogen safety engineering.

Wave properties of fourth-order fully implicit Runge–Kutta time integration schemes

A global spectral analysis of fully implicit, two-stage, fourth-order Runge–Kutta method (IRK2) is presented here. The essential numerical properties of the IRK2 time-integration scheme, viz., nodal amplification factor, phase speed and group velocity have been quantified and plotted at all the nodal points of the computational domain as functions of the CFL (Courant–Friedrichs–Lewy number) and wavenumber for various explicit and implicit spatial discretization schemes. Results show that IRK2 scheme produces neutral stability at all interior nodes for all CFL numbers. This has been shown that while using IRK2 scheme one has to use refined grid at high CFL numbers for dispersion relation preservation property.

Multi-stage solvers optimized for damping and propagation

Explicit multi-stage solvers are routinely used to solve the semi-discretized equations that arise in Computational Fluid Dynamics (CFD) problems. Often they are used in combination with multi-grid methods. In that case, the role of the multi-stage solver is to efficiently reduce the high frequency modes on the current grid and is called a smoother. In the past, when optimizing the coefficients of the scheme, only the damping characteristics of the smoother were taken into account and the interaction with the remainder of the multi-grid cycle was neglected. Recently it had been found that coefficients that result in less damping, but allow for a higher Courant–Friedrichs–Lewy (CFL) number are often superior to schemes that try to optimize damping alone. While this is certainly true for multi-stage schemes used as a stand-alone solver, we investigate in this paper if using higher CFL numbers also yields better results in a multi-grid setting. We compare the results with a previous study we conducted and where a more accurate model of the multi-grid cycle was used to optimize the various parameters of the solver. We show that the use of the more accurate model results in better coefficients and that in a multi-grid setting propagation is of little importance. We also look into the gains to be made when we allow the parameters to be different for the pre- and post-smoother and show that even better coefficients can be found in this way.

An Implicit Nonlinearly Consistent Method for the Two-Dimensional Shallow-Water Equations with Coriolis Force

Abstract An implicit and nonlinearly consistent (INC) solution technique is presented for the two-dimensional shallow-water equations. Since the method is implicit, and therefore unconditionally stable, time steps may be used that result in both gravity wave Courant–Friedrichs–Lewy (CFL) numbers and advection CFL numbers being larger than one. By nonlinearly consistent it is meant that all of the unknown variables appear at the same time level in the equations and are solved for simultaneously in an iterative manner (i.e., no splitting errors in time). The INC method is compared to a more traditional semi-implicit method for stepping over the gravity wave stability constraint. Results are presented that show that the second-order-in-time INC method can maintain a high level of accuracy if the dynamical timescale of the system is resolved by the time step. To investigate this difference in temporal integration accuracy between the nonlinearly consistent method and the semi-implicit method an approximate mo...