Courant Number Research Articles

Semi‐Lagrangian discontinuous Galerkin methods for scalar hyperbolic conservation laws

AbstractIn this article, we present a new class of semi‐Lagrangian methods for the numerical approximation of hyperbolic conservation laws. The methods provide a strategy for resolving nonlinear hyperbolic conservations laws within the framework of semi‐Lagrangian discontinuous Galerkin finite‐element methods (SLDG) initially proposed for linear advection problems [J‐M. Qiu and C‐W. Shu, J Comput Phys 230 (2011)]. Using a flux‐correction technique we combine the good numerical properties of the SLDG methods with those of the well‐known Runge‐Kutta discontinuous Galerkin methods (RKDG). The high‐order accuracy and local conservative properties of the resulting schemes make them very attractive for convection‐dominated convection‐diffusion problems, and conservation laws in particular. The semi‐Lagrangian component of the schemes allows for the use of larger Courant numbers than is possible when using high‐order RKDG methods alone.

High-order multirate explicit time-stepping schemes for the baroclinic-barotropic split dynamics in primitive equations

In order to treat the multiple time scales of ocean dynamics in an efficient manner, the baroclinic-barotropic splitting technique has been widely used for solving the primitive equations for ocean modeling. Based on the framework of strong stability-preserving Runge-Kutta approach, we propose two high-order multirate explicit time-stepping schemes (SSPRK2-SE and SSPRK3-SE) for the resulting split system in this paper. The proposed schemes allow for a large time step to be used for the three-dimensional baroclinic (slow) mode and a small time step for the two-dimensional barotropic (fast) mode, in which each of the two mode solves just need to satisfy their respective CFL conditions for numerical stability. Specifically, at each time step, the baroclinic velocity is first computed by advancing the baroclinic mode and fluid thickness of the system with the large time-step and the assistance of some intermediate approximations of the barotropic mode obtained by substepping with the small time step; then the barotropic velocity is corrected by using the small time step to re-advance the barotropic mode under an improved barotropic forcing produced by interpolation of the forcing terms from the preceding baroclinic mode solves; lastly, the fluid thickness is updated by coupling the baroclinic and barotropic velocities. Additionally, numerical inconsistencies on the discretized sea surface height caused by the mode splitting are relieved via a reconciliation process with carefully calculated flux deficits. Two benchmark tests from the “MPAS-Ocean” platform are carried out to numerically demonstrate the performance and parallel scalability of the proposed SSPRK-SE schemes.

Jacobian-Free Explicit Multiderivative Runge–Kutta Methods for Hyperbolic Conservation Laws

Based on the recent development of Jacobian-free Lax–Wendroff (LW) approaches for solving hyperbolic conservation laws (Zorio et al. in J Sci Comput 71:246–273, 2017, Carrillo and Parés in J Sci Comput 80:1832–1866, 2019), a novel collection of explicit Jacobian-free multistage multiderivative solvers for hyperbolic conservation laws is presented in this work. In contrast to Taylor time-integration methods, multiderivative Runge–Kutta (MDRK) techniques achieve higher-order of consistency not only through the excessive addition of higher temporal derivatives, but also through the addition of Runge–Kutta-type stages. This adds more flexibility to the time integration in such a way that more stable and more efficient schemes could be identified. The novel method permits the practical application of MDRK schemes. In their original form, they are difficult to utilize as higher-order flux derivatives have to be computed analytically. Here we overcome this by adopting a Jacobian-free approximation of those derivatives. In this paper, we analyze the novel method with respect to order of consistency and stability. We show that the linear CFL number varies significantly with the number of derivatives used. Results are verified numerically on several representative testcases.

A Two-Dimensional Stabilized Discontinuous Galerkin Method on Curvilinear Embedded Boundary Grids

In this work, we propose a state redistribution method for high order discontinuous Galerkin methods on curvilinear embedded boundary grids. State redistribution relaxes the overly restrictive CFL condition that results from arbitrarily small cut cells when explicit time steppers are used. Thus, the scheme can take take time steps that are proportional to the size of cells in the background grid. The discontinuous Galerkin scheme is stabilized by postprocessing the numerical solution after each stage or step of an explicit time stepping method. This is done by temporarily coarsening, or merging, the small cells into larger, possibly overlapping neighborhoods using a special weighted inner product. Then, the numerical solution on the neighborhoods is refined back onto the base grid in a conservative fashion. The advantage of this approach is that it uses only basic mesh information that is already available in many cut cell codes and does not require complex geometric manipulations. We prove that state redistribution is conservative and $p$-exact. Finally, we solve a number of test problems that demonstrate the encouraging potential of this technique for applications on curvilinear embedded geometries. Numerical experiments reveal that our scheme converges with order $p+1$ in $L_1$ and between $p$ and $p+1$ in $L_\infty$.

A Semi-Lagrangian Computation of Front Speeds of G-equation in ABC and Kolmogorov Flows with Estimation via Ballistic Orbits

The Arnold-Beltrami-Childress (ABC) flow and the Kolmogorov flow are three dimensional periodic divergence free velocity fields that exhibit chaotic streamlines. We are interested in front speed enhancement in G-equation of turbulent combustion by large intensity ABC and Kolmogorov flows. We give a quantitative construction of the ballistic orbits of ABC and Kolmogorov flows, namely those with maximal large time asymptotic speeds in a coordinate direction. Thanks to the optimal control theory of G-equation (a convex but non-coercive Hamilton-Jacobi equation), the ballistic orbits serve as admissible trajectories for front speed estimates. To study the tightness of the estimates, we compute the front speeds of G-equation based on a semi-Lagrangian (SL) scheme with Strang splitting and weighted essentially non-oscillatory (WENO) interpolation. Time step size is chosen so that the Courant number grows sublinearly with the flow intensity. Numerical results show that the front speed growth rate in terms of the flow intensity may approach the analytical bounds from the ballistic orbits.

Von Neumann Stability Analysis of DG-Like and PNPM-Like Schemes for PDEs with Globally Curl-Preserving Evolution of Vector Fields

This paper examines a class of involution-constrained PDEs where some part of the PDE system evolves a vector field whose curl remains zero or grows in proportion to specified source terms. Such PDEs are referred to as curl-free or curl-preserving, respectively. They arise very frequently in equations for hyperelasticity and compressible multiphase flow, in certain formulations of general relativity and in the numerical solution of Schrödinger’s equation. Experience has shown that if nothing special is done to account for the curl-preserving vector field, it can blow up in a finite amount of simulation time. In this paper, we catalogue a class of DG-like schemes for such PDEs. To retain the globally curl-free or curl-preserving constraints, the components of the vector field, as well as their higher moments, must be collocated at the edges of the mesh. They are updated using potentials collocated at the vertices of the mesh. The resulting schemes: (i) do not blow up even after very long integration times, (ii) do not need any special cleaning treatment, (iii) can operate with large explicit timesteps, (iv) do not require the solution of an elliptic system and (v) can be extended to higher orders using DG-like methods. The methods rely on a special curl-preserving reconstruction and they also rely on multidimensional upwinding. The Galerkin projection, highly crucial to the design of a DG method, is now conducted at the edges of the mesh and yields a weak form update that uses potentials obtained at the vertices of the mesh with the help of a multidimensional Riemann solver. A von Neumann stability analysis of the curl-preserving methods is conducted and the limiting CFL numbers of this entire family of methods are catalogued in this work. The stability analysis confirms that with the increasing order of accuracy, our novel curl-free methods have superlative phase accuracy while substantially reducing dissipation. We also show that PNPM-like methods, which only evolve the lower moments while reconstructing the higher moments, retain much of the excellent wave propagation characteristics of the DG-like methods while offering a much larger CFL number and lower computational complexity. The quadratic energy preservation of these methods is also shown to be excellent, especially at higher orders. The methods are also shown to be curl-preserving over long integration times.

Temporally adaptive conservative scheme for unsteady compressible flow

An adaptive conservative time integration scheme (ACTI) for compressible flow simulations is devised. A finite volume method is used, in which every cell employs a time step size which is some power of two times smaller than a global time step. By first advancing those cells with smaller time steps, this allows to use the already computed accumulated fluxes when the adjacent cells with larger time steps catch up. Conservation and periodic synchronization are thus guaranteed by construction; both important advantages of the new scheme (especially in the presence of shock waves) compared to previously proposed adaptive time stepping methods for compressible flow. Accuracy and computational speedup are demonstrated for challenging 1D and 2D test cases. Comparisons of ACTI results with an analytical solution and with results obtained without ACTI show excellent agreement. Although maximum wave speeds and/or cell sizes vary by several orders of magnitude throughout the computational domains, the CFL numbers employed by ACTI in each cell vary not more than by approximately a factor of two. Without ACTI, on the other hand, the CFL numbers are tiny in most of the cells. Since this results in a dramatic reduction of required flux computations, high speedup gains are achieved due to ACTI; in particular, if only few cells experience a very small time step, while the time step size limitations are less severe for the majority of cells.

A stable SPH model with large CFL numbers for multi-phase flows with large density ratios

The discontinuity across interface of multi-phase flows with large density ratios usually poses great challenges for numerical simulations. The smoothed particle hydrodynamics (SPH) is a meshless method with inherent advantages in dealing with multi-phase flows without the necessity of tracking the moving interfaces. In this paper, we develop a new weakly-compressible SPH model for multi-phase flows with large density ratios while allowing large CFL numbers. In the present SPH model, the continuity equation is first modified by eliminating the influence from particles of different phases based on the simple fact that different phases will not contribute when calculating the density for immiscible multi-phase flows; thus, the modified continuity equation will only consider the influence from neighboring particles of the same phase. The pressure and density of the particles of other phases are then re-initialized by using the Shepard interpolation function. The present multi-phase SPH model has been tested by four numerical examples, including the two-phase hydrostatic water, standing waves, liquid sloshing, and dam breaking. It has been demonstrated that the present multi-phase SPH model can obtain satisfactory results stably, even at large CFL numbers, and this means that large time steps can be employed. Therefore, the present multi-phase SPH model can significantly save computational cost through using large time steps, especially for large-scale problems with a large number of particles.

An implicit lattice Boltzmann flux solver for simulation of compressible flows

The lattice Boltzmann flux solver (LBFS) is a novel Boltzmann-typed solver under the finite volume framework, which has many advantages over traditional lattice Boltzmann equation (LBE) solvers. However, the existing versions of LBFS are explicit in temporal discretization. Therefore, they are of low computational efficiency and are mostly applied to the simulation of flows around simple geometries. In this paper, an implicit lattice Boltzmann flux solver for compressible flows is proposed, which can efficiently march in time with a large CFL number, and therefore facilitates the simulation of flows around complex geometries in practical engineering applications. Aiming to obtain the time-accurate results of unsteady flows, the implicit LBFS is further combined with the dual time-stepping technique to efficiently simulate unsteady flows, in which a sub-iteration is applied on the pseudo time level. In the sub-iteration, the Jacobian matrix is preconditioned by the lower-upper symmetric Gauss-Seidel (LUSGS) method and solved by the generalized minimal residual method (GMRES). The accuracy and efficiency of the proposed solver are demonstrated by various simulations in a wide Mach number regime. To further demonstrate its performance in practical engineering applications, the proposed solver is also used to simulate the flow over complex aircraft models.

Non-intrusive reduced order modeling of natural convection in porous media using convolutional autoencoders: Comparison with linear subspace techniques

Natural convection in porous media is a highly nonlinear multiphysical problem relevant to many engineering applications (e.g., the process of CO2 sequestration). Here, we extend and present a non-intrusive reduced order model of natural convection in porous media employing deep convolutional autoencoders for the compression and reconstruction and either radial basis function (RBF) interpolation or artificial neural networks (ANNs) for mapping parameters of partial differential equations (PDEs) on the corresponding nonlinear manifolds. To benchmark our approach, we also describe linear compression and reconstruction processes relying on proper orthogonal decomposition (POD) and ANNs. We present comprehensive comparisons among different models through three benchmark problems. The reduced order models, linear and nonlinear approaches, are much faster than the finite element model, obtaining a maximum speed-up of 7×106 because our framework is not bound by the Courant–Friedrichs–Lewy condition; hence, it could deliver quantities of interest at any given time contrary to the finite element model. Our model’s accuracy still lies within a relative error of 7% in the worst-case scenario. We illustrate that, in specific settings, the nonlinear approach outperforms its linear counterpart and vice versa. We hypothesize that a visual comparison between principal component analysis (PCA) and t-Distributed Stochastic Neighbor Embedding (t-SNE) could indicate which method will perform better prior to employing any specific compression strategy.

Moment-Based Multi-Resolution HWENO Scheme for Hyperbolic Conservation Laws

In this paper, a high-order moment-based multi-resolution Hermite weighted essentially non-oscillatory (HWENO) scheme is designed for hyperbolic conservation laws. The main idea of this scheme is derived from our previous work [J. Comput. Phys., 446 (2021) 110653], in which the integral averages of the function and its first order derivative are used to reconstruct both the function and its first order derivative values at the boundaries. However, in this paper, only the function values at the Gauss-Lobatto points in the one or two dimensional case need to be reconstructed by using the information of the zeroth and first order moments. In addition, an extra modification procedure is used to modify those first order moments in the troubled-cells, which leads to an improvement of stability and an enhancement of resolution near discontinuities. To obtain the same order of accuracy, the size of the stencil required by this moment-based multi-resolution HWENO scheme is still the same as the general HWENO scheme and is more compact than the general WENO scheme. Moreover, the linear weights can also be any positive numbers as long as their sum equals one and the CFL number can still be 0.6 whether for the one or two dimensional case. Extensive numerical examples are given to demonstrate the stability and resolution of such moment-based multi-resolution HWENO scheme.

A Semi-Implicit Meshless Method for Incompressible Flows in Complex Geometries

We present an exponentially convergent semi-implicit meshless algorithm for the solution of Navier-Stokes equations in complex domains. The algorithm discretizes partial derivatives at scattered points using radial basis functions as interpolants. Higher-order polynomials are appended to polyharmonic splines (PHS-RBF) and a collocation method is used to derive the interpolation coefficients. The interpolating kernels are then differentiated and the partial-differential equations are satisfied by collocation at the scattered points. The PHS-RBF interpolation is shown to be exponentially convergent with discretization errors decreasing as a high power of a representative distance between points. We present here a semi-implicit algorithm for time-dependent and steady state fluid flows in complex domains. At each time step, several iterations are performed to converge the momentum and continuity equations. A Poisson equation for pressure corrections is formulated by imposing divergence free condition on the iterated velocity field. At each time step, the momentum and pressure correction equations are repeatedly solved until the velocities and pressure converge to a pre-specified tolerance. We have demonstrated the convergence and discretization accuracy of the algorithm for two model problems and simulated three other complex problems. In all cases, the algorithm is stable for Courant numbers in excess of ten. The algorithm has the potential to accurately and efficiently solve many fluid flow and heat transfer problems in complex domains. An open source code Meshless Multi-Physics Software (MeMPhyS) is available for interested users of the algorithm.

Investigation of electrohydrodynamic calculations

A perfect dielectric model was incorporated into the OpenFOAM? software and used for investigation and, possibly, improvements of electrohydrodynamic calculations. Two different sets of numerical simulations were analyzed, in which two different fluids were present. The first set was one-dimensional, while in the second, a drop of one fluid was surrounded by the other fluid. It is shown that oscillations and possible artificial generation of a curl of the electric field strength can be observed at applying certain expressions or calculation strategies, which can be thus abandoned. Usage of dynamic meshes, at least those present in the used software, and of limiters for the gradient of the electric field strength can lead to large numerical errors. It is also shown that usage of certain cell face values could improve the results. An electric Courant number was derived by dimensional analysis, and it could be suggested for future calculations. Conclusions made in this paper are expected to be transferable to other more complicated models.

Mathematical modeling of mixed-traffic in urban areas

Transportation is the means of mobility. Due to the growth in the population, rising traffic on road, delay in the movement of vehicles and traffic chaos could be observed in urban areas. Traffic congestion causes many social and economic problems. Because of the convenience and the quickness, motor-bikes gradually become the main travel mode of urban cities. In this paper, we extend the Lighthill–Whitham–Richards (LWR) traffic flow model equation into the mixed-traffic flow of two entities: car and motor-bike in a unidirectional single-lane road segment. The flow of cars is modeled by the advection equation and the flow of motor-bikes is modeled by the advection-diffusion equation. The model equations for cars and motor-bikes are coupled based on total traffic density on the road section, and they are non-dimensionalized to introduce a non-dimensional number widely known as Péclet number. Explicit finite difference schemes satisfying the CFL conditions are employed to solve the model equations numerically to compute the densities of cars and motor-bikes. The simulation of densities over various time instants is studied and presented graphically. Finally, the average densities of cars and motor-bikes on the road section are calculated for various values of Péclet numbers and mixed-traffic behavior are discussed. It is observed that the mixed-traffic behavior of cars and motor-bikes depends upon the Péclet number. The densities of motor-bikes and cars in the mixed-traffic flow approach the equilibrium state earlier in time for smaller values of Péclet number whereas densities take longer time to approach the equilibrium for the greater values of Péclet number.

Implicit-explicit multistep formulations for finite element discretisations using continuous interior penalty

We consider a finite element method with symmetric stabilisation for the discretisation of the transient convection–diffusion equation. For the time-discretisation we consider either the second order backwards differentiation formula or the Crank–Nicolson method. Both the convection term and the associated stabilisation are treated explicitly using an extrapolated approximate solution. We prove stability of the method and the τ2 + hp+½ error estimates for the L2-norm under either the standard hyperbolic CFL condition, when piecewise affine (p = 1) approximation is used, or in the case of finite element approximation of order p ≥ 1, a stronger, so-called 4/3-CFL, i.e. τ ≤ Ch4/3. The theory is illustrated with some numerical examples.

On convergence of explicit finite volume scheme for one-dimensional three-component two-phase flow model in porous media

Abstract In this work, we develop and analyze an explicit finite volume scheme for a one-dimensional nonlinear, degenerate, convection–diffusion equation having application in petroleum reservoir. The main difficulty is that the solution typically lacks regularity due to the degenerate nonlinear diffusion term. We analyze a numerical scheme corresponding to explicit discretization of the diffusion term and a Godunov scheme for the advection term. L ∞ {L}^{\infty } stability under appropriate CFL conditions and BV {\rm{BV}} estimates are obtained. It is shown that the scheme satisfies a discrete maximum principle. Then we prove convergence of the approximate solution to the weak solution of the problem, and we mount convergence results to a weak solution of the problem in L 1 {L}^{1} . Results of numerical experiments are presented to validate the theoretical analysis.

The Pressure Peaking Phenomenon for Ignited Under-Expanded Hydrogen Jets in the Storage Enclosure: Experiments and Simulations for Release Rates of up to 11.5 g/s

This work focuses on the experimental and numerical investigation of maximum overpressure and pressure dynamics during ignited hydrogen releases in a storage enclosure, e.g., in marine vessel or rail carriage, with limited vent size area, i.e., the pressure peaking phenomenon (PPP) revealed theoretically at Ulster University in 2010. The CFD model previously validated against small scale experiments in a 1 m3 enclosure is employed here to simulate real-scale tests performed by the University of South-Eastern Norway (USN) in a chamber with a volume of 15 m3. The numerical study compares two approaches on how to model the ignited hydrogen release conditions for under-expanded jets: (1) notional nozzle concept model with inflow boundary condition, and (2) volumetric source model in the governing conservation equations. For the test with storage pressure of 11.78 MPa, both approaches reproduce the experimental pressure dynamics and the pressure peak with a maximum 3% deviation. However, the volumetric source approach reduces significantly the computational time by approximately 3 times (CFL = 0.75). The sensitivity analysis is performed to study the effect of CFL number, the size of the volumetric source and number of iterations per time step. An approach based on the use of a larger size volumetric source and uniform coarser grid with a mesh size of a vent of square size is demonstrated to reduce the duration of simulations by a factor of 7.5 compared to the approach with inflow boundary at the notional nozzle exit. The volumetric source model demonstrates good engineering accuracy in predicting experimental pressure peaks with deviation from −14% to +11% for various release and ventilation scenarios as well as different volumetric source sizes. After validation against experiments, the CFD model is employed to investigate the effect of cryogenic temperature in the storage on the overpressure dynamics in the enclosure. For a storage pressure equal to 11.78 MPa, it is found that a decrease of storage temperature from 277 K to 100 K causes a twice larger pressure peak in the enclosure due to the pressure peaking phenomenon.

Efficient implementation of characteristic-based schemes on unstructured triangular grids

Using characteristics to treat advection terms in time-dependent PDEs leads to a class of schemes, e.g., semi-Lagrangian and Lagrange–Galerkin schemes, which preserve stability under large Courant numbers, and may therefore be appealing in many practical situations. Unfortunately, the need of locating the feet of characteristics may cause a serious drop of efficiency in the case of unstructured space grids, and thus prevent the use of large time-step schemes on complex geometries. In this paper, we perform an in-depth analysis of the main recipes available for characteristic location, and propose a technique to improve the efficiency of this phase, using additional information related to the advecting vector field. This results in a clear improvement of execution times in the unstructured case, thus extending the range of applicability of large time-step schemes.

Unfitted hybrid high-order methods for the wave equation

We design an unfitted hybrid high-order (HHO) method for the wave equation. The wave propagates in a domain where a curved interface separates subdomains with different material properties. The key feature of the discretization method is that the interface can cut more or less arbitrarily through the mesh cells. We address both the second-order formulation in time of the wave equation and its reformulation as a first-order system. We prove H1-error estimates for a space semi-discretization in space of the second-order formulation, leading to optimal convergence rates for smooth solutions. Numerical experiments illustrate the theoretical findings and show that the proposed numerical schemes can be used to simulate accurately the propagation of acoustic waves in heterogeneous media, with meshes that are not fitted to the geometry. For the second-order formulation, the implicit, second-order accurate Newmark scheme is used for the time discretization, whereas (diagonally-implicit or explicit) Runge–Kutta schemes up to fourth-order accuracy are used for the first-order formulation. In the explicit case, we study the CFL condition on the time step and observe that the unfitted approach combined with local cell agglomeration leads to a comparable condition as when using fitted meshes.

Explicit high-order generalized-[formula omitted] methods for isogeometric analysis of structural dynamics

We propose a new family of high-order explicit generalized-α methods for hyperbolic problems with the feature of dissipation control. Our k-stage approach delivers 2k,k∈N accuracy order in time by solving k matrix systems explicitly and updating 2k variables at each time-step. The user can control the numerical dissipation in the discrete spectrum’s high-frequency regions by adjusting the method’s coefficients. We study the method’s spectral behaviour and show that the CFL condition is independent of the accuracy. The stability region remains invariant while increasing the number of stages, that is, the accuracy order. Next, we exploit efficient preconditioners for the isogeometric matrix to minimize the computational cost. These preconditioners use a diagonal-scaled Kronecker product of univariate parametric mass matrices; they have a robust performance with respect to the spline degree and the mesh size, and their decomposition structure implies that their application is faster than a matrix–vector product involving the fully-assembled mass matrix. Our high-order schemes require simple modifications of the available implementations of the generalized-α method. Finally, we present numerical examples demonstrating the methodology’s performance regarding single- and multi-patch IGA discretizations.