Courant Number Research Articles

On the stability of strong-stability-preserving modified Patankar–Runge–Kutta schemes

In this paper, we perform a stability analysis for classes of second and third order accurate strong-stability-preserving modified Patankar–Runge–Kutta (SSPMPRK) schemes, which were introduced in Huang and Shu [J. Sci. Comput. 78 (2019) 1811–1839] and Huang et al. [J. Sci. Comput. 79 (2019) 1015–1056] and can be used to solve convection equations with stiff source terms, such as reactive Euler equations, with guaranteed positivity under the standard CFL condition due to the convection terms only. The analysis allows us to identify the range of free parameters in these SSPMPRK schemes in order to ensure stability. Numerical experiments are provided to demonstrate the validity of the analysis.

A new Runge–Kutta–Chebyshev Galerkin-characteristic finite element method for advection–dispersion problems in anisotropic porous media

We propose a new approach that combines the modified method of characteristics with a unified finite element discretization for the numerical solution of a class of coupled Darcy-advection–dispersion problems in anisotropic porous media. The proposed method benefits from advantages of the method of characteristics in its ability to handle the nonlinear convective terms, while taking advantage of a unified formulation that allows the use of equal-order finite element approximations along with the L2-projection method for all solutions in the problem. In the proposed Galerkin-characteristic finite element framework, the standard Courant–Friedrichs–Lewy condition is relaxed, and time truncation errors are reduced since no stability criterion restricts the choice of time step. For time integration, we use a Runge–Kutta scheme with the Chebyshev polynomials (RKC). The RKC method has an extensive stability field and it is explicit and second-order accurate. In order to assess the quality of the proposed approach, we perform numerical experiments for several categories of test cases. The first category concerns accuracy test examples with known analytical solutions, while the second category focuses on a benchmark problem of miscible flow in a heterogeneous porous medium. The numerical results obtained show high performance and demonstrate the ability of the method to capture dispersion effects in porous media problems.

Provably Positive Central Discontinuous Galerkin Schemes via Geometric Quasilinearization for Ideal MHD Equations

In the numerical simulation of ideal magnetohydrodynamics (MHD), keeping the pressure and density always positive is essential for both physical considerations and numerical stability. This is however a challenging task, due to the underlying relation between such a positivity-preserving (PP) property and the magnetic divergence-free (DF) constraint as well as the strong nonlinearity of the MHD equations. In this paper, we present the first rigorous PP analysis of the central discontinuous Galerkin (CDG) methods and construct arbitrarily high-order provably PP CDG schemes for ideal MHD. By the recently developed geometric quasilinearization (GQL) approach, our analysis reveals that the PP property of standard CDG methods is closely related to a discrete magnetic DF condition, whose form was unknown prior to our analysis and differs from that for the noncentral discontinuous Galerkin (DG) and finite volume methods in [K. Wu, SIAM J. Numer. Anal., 56 (2018), pp. 2124–2147]. The discovery of this relation lays the foundation for the design of our PP CDG schemes. In the one-dimensional (1D) case, the discrete DF condition is naturally satisfied, and we rigorously prove that the standard CDG method is PP under a condition that can be enforced easily with an existing PP limiter. However, in the multidimensional cases, the corresponding discrete DF condition is highly nontrivial, yet critical, and we analytically prove that the standard CDG method, even with the PP limiter, is not PP in general, as it generally fails to meet the discrete DF condition. We address this issue by carefully analyzing the structure of the discrete divergence terms and then constructing new locally DF CDG schemes for Godunov’s modified MHD equations with an additional source term. The key point is to find out the suitable discretization of the source term such that it exactly cancels out all the terms in the discovered discrete DF condition. Based on the GQL approach, we prove in theory the PP property of the new multidimensional CDG schemes under a CFL condition. The robustness and accuracy of the proposed PP CDG schemes are further validated by several demanding 1D, two-dimensional, and three-dimensional numerical MHD examples, including the high-speed jets and blast problems with very low plasma beta.

Time-Accurate Solution of Unsteady Flows in an Implicit Solver Using Block LUSGS Method

Implicit solvers based on the lower-upper symmetric Gauss–Seidel (LUSGS) method are widely used in Computational Fluid Dynamics (CFD) codes to obtain steady-state solutions using large false-transient time steps. However, the method destroys time accuracy and thus cannot be directly used for true transient problems. For such problems, the LUSGS method has to be used with a dual time-stepping algorithm. In this paper, a modified method that preserves time accuracy and is called the block LUSGS (B-LUSGS) solver is employed to compute unsteady flows using larger time steps. This method is compared with the mentioned dual time-stepping technique used for transient problems. Numerical solutions obtained for Stokes' second problem simulation underline the superior temporal accuracy of the block LUSGS method in computing unsteady flows, even at CFL numbers close to 30. Simulations of 2D laminar and 3D turbulent vortex shedding downstream of a circular cylinder further validate the modified method in accurately predicting unsteady separated flows at a CFL number of 100. Comparisons of overall simulation CPU times show that the B-LUSGS method is at least six times computationally cheaper than the dual time-stepping algorithm that uses the traditional LUSGS solver. Abbreviations: BDF, Backward differentiation formula; CFD, Computational fluid dynamics; CFL Courant–Friedrichs–Lewy; CPU, Central processing unit; DTS, Dual time-stepping; LES, Large-eddy simulation; LUSGS, Lower-upper symmetric Gauss–Seidel; NASA, National Aeronautics and Space Administration; PDE, Partial differential equation; RAM, Random access memory; RANS, Reynolds-averaged Navier–Stokes; SA, Spalart–Allmaras; SST, Shear stress transport.

Efficient dissipation-based nonlinear solver for multiphase flow in discrete fractured media

Solving nonlinear equation systems resulting from the Fully Implicit Method (FIM) remains a challenge for numerical simulation of multiphase flow in subsurface fractured media. The Courant numbers can vary orders of magnitude across discrete fracture-matrix (DFM) models because of the high contrasts in the permeability and length-scale between matrix and fracture. The standard Newton solver is usually unable to converge for big timestep sizes.Limited research has been conducted on nonlinear solver techniques for multiphase flow and transport in fractured media. We make an extension of a new dissipation-based continuation (DBC) method to DFM models. Our goal is to prevent timestep cuts and maintain efficient time-stepping for FIM. The DBC algorithm builds a homotopy of the discretized conservation equations through the addition of numerical dissipation terms. We introduce a continuation parameter for controlling the dissipation and ensuring that accuracy of the computed solution will not be reduced. Under the nonlinear framework of DBC, general dissipation operators and adaptive strategies are developed to provide the optimal dissipation matrix for the multiphase hyperbolic system.We assess the new nonlinear solver using multiple numerical examples. Results reveal that the damped-Newton solver suffers from serious restrictions on timestep sizes and wasted iterations. In contrast, the DBC solver provides excellent computational performance. The dissipation operators are able to successfully resolve main convergence difficulties. We also investigate the impact of star-delta transformation which removes the small cells at fracture intersections. Moreover, we demonstrate that an aggressive time-stepping does not affect the solution accuracy.

Finite Volume Method for Modeling the Load-Rejection Process of a Hydropower Plant with an Air Cushion Surge Chamber

The pipe systems of hydropower plants are complex and feature special pipe types and various devices. When the Method of Characteristics (MOC) is used, interpolation or wave velocity adjustment is required, which may introduce calculation errors. The second-order Finite Volume Method (FVM) was presented to simulate water hammer and the load-rejection process of a hydropower plant with an air cushion surge chamber, which has rarely been considered before. First, the governing equations were discretized by FVM and the flux was calculated by a Riemann solver. A MINMOD slope limiter was introduced to avoid false oscillation caused by data reconstruction. The virtual boundary strategy was proposed to simply and effectively handle the complicated boundary problems between the pipe and the various devices, and to unify the internal pipeline and boundary calculations. FVM results were compared with MOC results, exact solutions, and measured values, and the sensitivity analysis was conducted. When the Courant number was equal to 1, the results of FVM and MOC were consistent with the exact solution. When the Courant number was less than 1, compared with MOC, the second-order FVM results were more accurate with less numerical dissipation. As the Courant number gradually decreased, the second-order FVM simulations were more stable. For the given numerical accuracy, second-order FVM had higher computational efficiency. The simulations of load rejection showed that compared with the MOC results, the second-order FVM calculations were closer to the measured values. For hydropower plants with complex pipe systems, wave velocity or the Courant number should be adjusted during MOC calculation, resulting in calculation error, and the error value is related to the parameters of the air cushion surge chamber (initial water depth, air cushion height, etc.). The second-order FVM can more accurately, stably, and efficiently simulate the load-rejection process of hydropower plants compared with MOC.

Relaxation of the Courant Condition in the Explicit Finite-Difference Time-Domain Method With Higher-Degree Differential Terms

A new explicit and nondissipative finite-difference time-domain (FDTD) method in two and three dimensions is proposed for the relaxation of the Courant condition. The third-degree spatial difference terms with second- and fourth-order accuracies are added with coefficients to the time-development equations of FDTD(2,4). Optimal coefficients are obtained by a brute-force search of the dispersion relations, which reduces phase velocity errors but satisfies the numerical stabilities as well. The new method is stable with large Courant numbers, whereas the conventional FDTD methods are unstable. The new method also has smaller numerical errors in the phase velocity than conventional FDTD methods with small Courant numbers.

An energy-based discontinuous Galerkin method with tame CFL numbers for the wave equation

We extend and analyze the energy-based discontinuous Galerkin method for second order wave equations on staggered and structured meshes. By combining spatial staggering with local time-stepping near boundaries, the method overcomes the typical numerical stiffness associated with high order piecewise polynomial approximations. In one space dimension with periodic boundary conditions and suitably chosen numerical fluxes, we prove bounds on the spatial operators that establish stability for CFL numbers \(c \frac{\Delta t}{h} < C\) independent of order when stability-enhanced explicit time-stepping schemes of matching order are used. For problems on bounded domains and in higher dimensions we demonstrate numerically that one can march explicitly with large time steps at high order temporal and spatial accuracy.

On the Stability of Robust Dynamical Low-Rank Approximations for Hyperbolic Problems

The dynamical low-rank approximation (DLRA) is used to treat high-dimensional problems that arise in such diverse fields as kinetic transport and uncertainty quantification. Even though it is well known that certain spatial and temporal discretizations when combined with the DLRA approach can result in numerical instability, this phenomenon is poorly understood. In this paper we perform an stability analysis for the corresponding nonlinear equations of motion. This reveals the source of the instability for the projector-splitting integrator [C. Lubich and I. V. Oseledets, BIT, 54 (2014), pp. 171–188] when first discretizing the equations and then applying the DLRA. Based on this we propose a projector-splitting integrator, based on applying DLRA to the continuous system before performing the discretization, that recovers the classic CFL condition. We also show that the unconventional integrator [G. Ceruti and C. Lubich, BIT, 62 (2021), pp. 23–44] has more favorable stability properties. Moreover, we explain why the projector-splitting integrator performs better when approximating higher moments, while the unconventional integrator is generally superior for first order moments. Furthermore, an efficient and stable dynamical low-rank update for the scattering term in kinetic transport is proposed. Numerical experiments for kinetic transport and uncertainty quantification, which confirm the results of the stability analysis, are presented.

Towards a full-scale CFD guideline for simulating a ship advancing in open water

ABSTRACT Computational Fluid Dynamics (CFD) simulations of a ship’s operations are generally conducted at model scale, but the reduced scale changes the fluid behaviour around the ship. Whilst ideally ship simulations should be run directly at full scale, a guide has not been published to advise on the suitable setups that can provide accurate results while minimizing the computational cost. To address this, the present work explores an optimal approach for full-scale ship simulations. Extensive sensitivity studies were conducted on relevant computational setups to investigate their influences on the prediction of ship resistance, ship-generated waves as well as the boundary-layer flow of the hull. A set of CFD setups for full-scale ship simulations in open water was recommended. It was demonstrated that the ideal Y+ and Courant numbers in full scale are evidently different from those given in current model-scale CFD guidelines, indicating the necessity to establish full-scale CFD guidelines separately.

Direct modeling for computational fluid dynamics and the construction of high-order compact scheme for compressible flow simulations

Computational fluid dynamics is a direct modeling of physical laws in a discretized space. The basic physical laws include the mass, momentum and energy conservations, numerical scale-consistent transport process, and similar domain of dependence and influence in the physical and the numerical wave propagations. Therefore, a physically soundable numerical scheme should be a compact one which involves the closest neighboring cells within the domain of dependence for the solution update under a CFL number (∼1). In order to keep the compactness of the scheme and be high-order, subcell flow distributions or the equivalent degrees of freedom beyond the cell averaged flow variables are needed, such as the gradients of the flow variables inside each control volume. Under the above consideration, the direct use of Riemann solver as the evolution model doesn't provide such a mechanism. High-order dynamic evolution process should be used to update the time accurate cell-averaged flow variables and their gradients. The direct modeling of flow evolution under generalized initial condition will be developed. In order to provide reliable cell averaged flow variables and their gradients for the compact data reconstruction, the evolution process is constructed to provide possible discontinuous time-dependent flow variables across a cell interface. At the same time, the time accurate flux function at a cell interface can become a discontinuous function of time in its evolution process. Same as the spatial limiter in the conventional computational fluid dynamics methods, such as the total variation diminishing and weighted essentially non-oscillatory methods, the temporal limiter for the time accurate flux function is designed as well. The direct modeling method in this paper will update flow variables differently on both sides of a cell interface and nonlinearly limit the high-order time derivatives of the flux function in case of shock wave passing through a cell interface within a time step. The direct modeling method unifies the concepts of nonlinear limiters in spatial data reconstruction and temporal flux transport. The equivalent treatment of space and time makes the compact gas-kinetic scheme be super robust and accurate for the compressible flow simulation. At the same time, the 4th to 8th-order compact schemes can use a large CFL number (∼0.8) in flow computation, even for the hypersonic flow simulation.

Fully discrete WENO with double entropy condition for hyperbolic conservation laws

ABSTRACTThis paper put forward a new fully discrete scheme construction method – double entropy condition solution formula method. With that, we turn the state-of-the-art semi-discrete WENO + RK scheme into a fully discrete scheme, which is named as Full-WENO. A major difficulty of this work is that we lack exact solution expressions for nonlinear equations in general cases. A feasible way we can go is to linearize equations and get quasi-exact solution formulas. The critical challenge is keeping both accuracy and efficiency in a scheme. Then, we get a class of new high-order schemes far better than traditional WENO schemes in the following aspects: (1) One-step to consistent high accuracy order in both space and time; (2) Resolution improves with the increasing CFL number; (3) Less CPU time and memory space, 1/s times of WENO with s-stage RK method in theory; (4) Excellent entropy condition satisfying property. Compared with our original work , the new method applies the more sophisticated WENO reconstruction and solves the resolution loss problems in multi-dimensional cases. The numerical tests show that the new scheme is equipped with the merits of high efficiency, high resolution and low dissipation, especially for long-time nonlinear problems.

High-order compact gas-kinetic schemes for three-dimensional flow simulations on tetrahedral mesh

A general framework for the development of high-order compact schemes has been proposed recently. The core steps of the schemes are composed of the following. 1). Based on a kinetic model equation, from a generalized initial distribution of flow variables construct a time-accurate evolution solution of gas distribution function at a cell interface and obtain the corresponding flux function; 2). Introduce the WENO-type weighting functions into the high-order time-derivative of the cell interface flux function in the multistage multi-derivative (MSMD) time stepping scheme to cope with the possible impingement of a shock wave on a cell interface within a time step, and update the cell-averaged conservative flow variables inside each control volume; 3). Model the time evolution of the gas distribution function on both sides of a cell interface separately, take moments of the inner cell interface gas distribution function to get flow variables, and update the cell-averaged gradients of flow variables inside each control volume; 4). Based on the cell-averaged flow variables and their gradients, develop compact initial data reconstruction to get initial condition of flow distributions at the beginning of next time step. A compact gas-kinetic scheme (GKS) up to sixth-order accuracy in space and fourth-order in time has been constructed on 2D unstructured mesh. In this paper, the compact GKS up to fourth-order accuracy on three-dimensional tetrahedral mesh will be further constructed with the focus on the WENO-type initial compact data reconstruction. Nonlinear weights are designed to achieve high-order accuracy for the smooth Navier-Stokes solution and keep super robustness in 3D computation with strong shock interactions. The fourth-order compact GKS uses a large time step with a CFL number 0.6 in the simulations from subsonic to hypersonic flow. A series of test cases are used to validate the scheme. The high-order compact GKS can be used in 3D applications with complex geometry.

Is the classic convex decomposition optimal for bound-preserving schemes in multiple dimensions?

Since proposed in Zhang and Shu (2010) [1], the Zhang–Shu framework has attracted extensive attention and motivated many bound-preserving (BP) high-order discontinuous Galerkin and finite volume schemes for various hyperbolic equations. A key ingredient in the framework is the decomposition of the cell averages of the numerical solution into a convex combination of the solution values at certain quadrature points, which helps to rewrite high-order schemes as convex combinations of formally first-order schemes. The classic convex decomposition originally proposed by Zhang and Shu has been widely used over the past decade. It was verified, only for the 1D quadratic and cubic polynomial spaces, that the classic decomposition is optimal in the sense of achieving the mildest BP CFL condition. Yet, it remained unclear whether the classic decomposition is optimal in multiple dimensions. In this paper, we find that the classic multidimensional decomposition based on the tensor product of Gauss–Lobatto and Gauss quadratures is generally not optimal, and we discover a novel alternative decomposition for the 2D and 3D polynomial spaces of total degree up to 2 and 3, respectively, on Cartesian meshes. Our new decomposition allows a larger BP time step-size than the classic one, and moreover, it is rigorously proved to be optimal to attain the mildest BP CFL condition, yet requires much fewer nodes. The discovery of such an optimal convex decomposition is highly nontrivial yet meaningful, as it may lead to an improvement of high-order BP schemes for a large class of hyperbolic or convection-dominated equations, at the cost of only a slight and local modification to the implementation code. Several numerical examples are provided to further validate the advantages of using our optimal decomposition over the classic one in terms of efficiency.

Challenges & solutions for deterministic hydraulic modelling in Mediterranean coastal catchment. Application to the lower Paillons River, Nice, France

Different methods are used in hydrological and hydraulic modelling in a catchment. Deterministic modelling can be used to produce detailed information and optimize the organization of massive data collection. In the catchment of Paillons, there is a need to improve knowledge of hydrological processes through deterministic modelling. In fact, Nice Côte d’Azur metropolis faces numerous challenges in the Paillons catchment. The area provides water resources for large communities. However, it is exposed to flooding and droughts. Complex hydrological processes generate runoff over the 246 km2 watershed. Few existing monitoring stations provide runoff data. In addition, there are limitations in the understanding of surface hydraulics. Thus, this study uses DHI Mike21FM for surface hydraulics with 5 m grid size for riverbed and 2 m for tunnels. The study area is limited to the lower Paillons River. Observed and modelled water depths vary between 0.1 and 0.5 m for a maximum discharge of 38 m3/s in the tunnels. Flood maps created with a discharge of 1000 m3/s, show clearly high flood risk zones and flow directions. The selected CFL condition under 0.8 is respected. The tool is suitable for modeling flooding in areas of interest within the catchment of study. The results obtained are satisfactory and demonstrate that the constructed tool makes it possible to reproduce the overall behavior of surface hydraulics.

Modeling of two-phase flows at low Capillary number with VoF method

In this work, a modified version of the FSF (Filtered Surface Force) method is presented. The method was originally developed by Raeini et al. specifically for modeling low Capillary number multiphase flows, basing its structure on a well consolidated VoF (Volume of Fluid) CSF (Continuum Surface Force) solver in OpenFOAM®. The numerical method proposed in this paper is a combination of the FSF and the CICSAM scheme of Ubbink, which has been translated into a TVD framework instead of its original NVD formulation. The capability of this scheme in compressing the interface numerically lets to simplify the FSF formulation, avoiding the modification of the indicator function for the computation of the surface tension forces. Moreover, a different formulation of the filter applied to the capillary forces is presented as a simplification of the original implementation. The time step is limited by capillary forces time scale and the Courant number is chosen as a compromise between the necessity of reducing the calculation time and the numerical diffusion of the advection schemes. The methodology is validated with several test cases, showing different problems of the classical CSF method and making comparisons with theoretical solutions.

A mass conservative, well balanced and positivity-preserving central scheme for shallow water equations

Based on the invariant-region-preserving (IRP) reconstruction method introduced in [Yan, Tong and Chen, Appl. Math. Comput., 436 (2023) 127500], a second order unstaggered central scheme is proposed to solve the shallow water equations with bottom topography in the framework that the bottom is discretized by a continuous, piecewise linear approximation. The reconstruction applies a modification locally on a preliminary reconstructed surface gradient in every cell to yield a convexity property of the sampled point value in the forward and backward projections. The water mass conservation is proved by rewriting the scheme in a conservation form. The modification does not change the preliminary reconstructed slope of water surface for the lake-at-rest steady state and then keeps the well-balancing property of the surface gradient method. The convexity property ensures the nonnegativity of the updated water depth under a large CFL number which yields a considerable speed-up. The numerical experiments are shown to demonstrate the robustness of the scheme.

Adaptively implicit MPDATA advection for arbitrary Courant numbers and meshes

AbstractAdvection schemes with time step restrictions are widely used in weather and climate models. This can lead to instability in the presence of high flow speeds (relative to mesh spacing) such as occurs in convective updraughts, regions of mesh convergence, or where the winds are unusually high. An adaptively implicit advection scheme is proposed that treats advection implicitly only where the Courant number is high. Flux correction to ensure monotonicity is adapted to work with implicit time stepping. A version of the multidimensional positive‐definite advection transport algorithm MPDATA is derived with an anti‐diffusive flux compensating truncation errors of off‐centred implicit time stepping. The anti‐diffusive flux is gradually reduced as Courant numbers increase above 2 in order to maintain stability at the expense of second‐order accuracy at high Courant numbers. Results of two‐dimensional advection by deformational flow are presented on various meshes of the sphere. Stability and second‐order accuracy are maintained when the Courant number is over 100 in a small region, when strong wind crosses the poles of a rotated latitude–longitude mesh. Good solutions are also obtained on a skipped latitude–longitude mesh, a cubed sphere, and hexagonal meshes. Accuracy reverts to first order when Courant numbers are large over a large fraction of the domain.

Spatially Adaptive Projective Integration Schemes For Stiff Hyperbolic Balance Laws With Spectral Gaps

Stiff hyperbolic balance laws exhibit large spectral gaps, especially if the relaxation term significantly varies in space. Using examples from rarefied gases and the general form of the underlying balance law model, we perform a detailed spectral analysis of the semi-discrete model that reveals the spectral gaps. Based on that, we show the inefficiency of standard time integration schemes expressed by a severe restriction of the CFL number. We then develop the first spatially adaptive projective integration schemes to overcome the prohibitive time step constraints of standard time integration schemes. The new schemes use different time integration methods in different parts of the computational domain, determined by the spatially varying value of the relaxation time. We use our analytical results to derive accurate stability bounds for the involved parameters and show that the severe time step constraint can be overcome. The new adaptive schemes show good accuracy in a numerical test case and can obtain a large speedup with respect to standard schemes.

Temporal accuracy of FastRK3

FastRK3 by Aithal and Ferrante (2020) [1] is an explicit, three-stage, third-order Runge-Kutta (RK3) based projection-method for solving the incompressible Navier-Stokes equations written in curvilinear coordinates using the orthogonal formulation. FastRK3 requires solving the Poisson equation for pressure only once per time step versus three times of standard RK3 methodology. In this work, we report the results of a theoretical and numerical study which is focused on the temporal accuracy of FastRK3. For free-shear flows, we derive theoretically and show numerically that FastRK3 is third-order accurate in time for velocity and second-order accurate in time for pressure when the RK3 coefficients and the pressure extrapolation scheme satisfy specific conditions herein derived. For ‘stiff’ wall-bounded flows, we show numerically that the temporal accuracy of velocity and pressure is the same, i.e., first- or second-order depending on the CFL number, for both FastRK3 and standard RK3. In summary, we derive theoretically and verify numerically the conditions that the RK3 coefficients and the pressure extrapolation scheme need to satisfy in order for FastRK3 to preserve the temporal accuracy of the underlying RK3 methodology.