Courant Number Research Articles

Exponential Stability of an Upwind Difference Splitting Scheme for Symmetric t-Hyperbolic Systems

This article considers mixed problem for 2-dimensional symmetric t-hyperbolic system. This system consists of constant coefficients. We investigated questions of setting mixed problems (problems with initial-boundary conditions) for symmetric t-hyperbolic system. It is known that for a symmetric t-hyperbolic system, the concept of characteristic velocity plays a very important role. The main point of this article is to construct a difference control scheme for these characteristic velocities. Control parameters are introduced to control the characteristic velocities of a symmetric t-hyperbolic system. The value of the parameter, the characteristic velocities can have different signs. We will consider all possible cases. In these cases, it is obvious that a symmetric t-hyperbolic system with n-equations can be in (n+1)(n+1) different states. These states controlled by control parameters. Upwind difference schemes constructed for each of these states depending on the value of characteristic velocities. The suggested difference schemes for the mixed problem were substantiated is exponential stability by the Lyapunov. It will be shown that the difference scheme will be stable if the CFL condition is satisfied and otherwise unstable.

An asymptotic-preserving and exactly mass-conservative semi-implicit scheme for weakly compressible flows based on compatible finite elements

We present a novel asymptotic-preserving semi-implicit finite element method for weakly compressible and incompressible flows based on compatible finite element spaces. The momentum is sought in an H(div)-conforming space, ensuring exact pointwise mass conservation at the discrete level. We use an explicit discontinuous Galerkin-based discretization for the nonlinear convective terms, while treating the pressure and viscous terms implicitly, so that the CFL condition depends only on the fluid velocity. To handle shocks and damp spurious oscillations in the compressible regime, we incorporate an a posteriori limiter that employs artificial viscosity and is based on a discrete maximum principle. By using hybridization, the final algorithm requires solving only symmetric positive definite linear systems. As the Mach number approaches zero and the density remains constant, the method naturally converges to an H(div)-based discretization of the incompressible Navier-Stokes equations in the vorticity-velocity-pressure formulation. Several numerical tests validate the proposed method.

An energy stable high-order cut cell discontinuous Galerkin method with state redistribution for wave propagation

Cut meshes are a type of mesh that is formed by allowing embedded boundaries to “cut” a simple underlying mesh resulting in a hybrid mesh of cut and standard elements. While cut meshes can allow complex boundaries to be represented well regardless of the mesh resolution, their arbitrarily shaped and sized cut elements can present issues such as the small cell problem, where small cut elements can result in a severely restricted CFL condition. State redistribution, a technique developed by Berger and Giuliani in [1], can be used to address the small cell problem. In this work, we pair state redistribution with a high-order discontinuous Galerkin scheme that is L2 energy stable under arbitrary quadrature. We prove that state redistribution can be added to a provably L2 energy stable discontinuous Galerkin method on a cut mesh without damaging the scheme's L2 stability. We numerically verify the high order accuracy and stability of our scheme on two-dimensional wave propagation problems.

Multigrid Reduction‐In‐Time Convergence for Advection Problems: A Fourier Analysis Perspective

ABSTRACTA long‐standing issue in the parallel‐in‐time community is the poor convergence of standard iterative parallel‐in‐time methods for hyperbolic partial differential equations (PDEs), and for advection‐dominated PDEs more broadly. Here, a local Fourier analysis (LFA) convergence theory is derived for the two‐level variant of the iterative parallel‐in‐time method of multigrid reduction‐in‐time (MGRIT). This closed‐form theory allows for new insights into the poor convergence of MGRIT for advection‐dominated PDEs when using the standard approach of rediscretizing the fine‐grid problem on the coarse grid. Specifically, we show that this poor convergence arises, at least in part, from inadequate coarse‐grid correction of certain smooth Fourier modes known as characteristic components, which was previously identified as causing poor convergence of classical spatial multigrid on steady‐state advection‐dominated PDEs. We apply this convergence theory to show that, for certain semi‐Lagrangian discretizations of advection problems, MGRIT convergence using rediscretized coarse‐grid operators cannot be robust with respect to CFL number or coarsening factor. A consequence of this analysis is that techniques developed for improving convergence in the spatial multigrid context can be re‐purposed in the MGRIT context to develop more robust parallel‐in‐time solvers. This strategy has been used in recent work to great effect; here, we provide further theoretical evidence supporting the effectiveness of this approach.

Flux reconstruction approach for long-time calculations of the linearized Euler equation

The Linearized Euler Equations (LEE) are one of the fundamental governing equations in acoustics research. Due to the relatively small magnitude of acoustic disturbances compared to the physical quantities of the background flow, ensuring highly precise solutions without dissipation and dispersion over time is essential. The Flux Reconstruction (FR) method, recognized as a high-order numerical method with flexible construction capabilities, has been extensively utilized in recent years for solving partial differential equations numerically. This paper introduces an offsetting numerical flux for one-dimensional LEE. Subsequently, it develops a flux reconstruction scheme with energy-conserving characteristics, referred to as the ECFR scheme, aimed at achieving high-precision numerical solutions for prolonged computations. The energy-conserving feature of the scheme is substantiated through theoretical derivation and further corroborated by numerical validation. Moreover, the impact of Mach numbers and the parameter of the offsetting flux on the scheme's performance is discussed. The stability of the scheme is also analyzed, with the limit of the CFL number being provided. The novel ECFR scheme facilitates obtaining high-precision acoustic numerical solutions, thereby providing strong support for the study of complex engineering problems such as aircraft, submarine, and automobile noise.

Surprisingly tight Courant–Friedrichs–Lewy condition in explicit high-order Arakawa schemes

The Arakawa scheme, which exactly conserves energy and enstrophy in two-dimensional hydrodynamics, severely reduces the maximum allowable stable time step compared to the Courant–Friedrichs–Lewy condition when combined with a spectral or spectral-like high-order discretization. We calculate the time step restriction by finding the responsible high frequency eigenfunctions and suggest a remedy. The phenomenon has relevance for analogous methods conserving multiple quadratic invariants in gyrokinetic or magnetohydrodynamic simulations.

Investigation on transition characteristics of hydrofoil boundary layer based on algebraic local-correlation-based transition modeling model

Hydrofoil shapes are used for the marine turbine blades to capture kinetic energy from water currents effectively. Predicting transitions is a critical concern when studying the hydrofoil boundary layer. This paper analyzed the transitional behavior of the boundary layer in the National Advisory Committee of Aeronautics (NACA) hydrofoil, NACA0009, with a blunt trailing edge using the Algebraic Local-Correlation-based Transition Modeling (Algebraic LCTM) model. First, through sensitivity analysis, the effects of the maximum y+ (the dimensionless distance y to the wall), grid expansion ratio, number of normal and streamlined grids, and timescale on transition prediction were studied. The results indicate that finer y+ value and appropriate grid expansion ratios can improve the accuracy of transition prediction, while the influence of timescale on the prediction results is relatively small within the range of Courant number theory values. Second, further analysis was conducted on the transition prediction performance under different Reynolds numbers. It was found that the model predictions were consistent with experimental values at low Reynolds numbers, but the predicted transition position was advanced at high Reynolds numbers, mainly because of the significant disparity in eddy viscosity coefficients within the free flow field. In the study of leading-edge roughness bands' impact on boundary layer transition for hydrofoil, the introduction of roughness significantly expedited the transition process. The Algebraic LCTM model outperformed the gamma (γ) transition model, reducing prediction errors by 5–40% for boundary layer parameters and maintaining errors between 0.005 and 4% for wake vortex shedding frequency, as opposed to the γ model's 0–23%. The results of this study provide a theoretical basis for hydrofoil design.

A convection-diffusion-reaction system with discontinuous flux modelling biofilm growth in slow sand filters

A Slow Sand Filter (SSF) consists of biofilm attached to saturated packed sand with supernatant water on top, through which water flows by gravity. SSFs are used in the production of drinking water to improve water quality and remove pathogens, and, as they use biological filtration via the microbes present in the biofilm, they are simple and cheap to construct. A drawback of SSFs is that the difficult sampling and the complex dynamics associated with biofilm growth of microbes pose an obstacle for analysis and predictions. In order to overcome this issue, a mathematical model consisting of a non-linear convection-diffusion-reaction system with discontinuous flux is developed from first principles as a unified framework for the dynamics in time and depth of the biofilm and microbial community, where a variety of ecological models can be included in the reaction terms, which also contain attachment, detachment and transfer processes. The multi-phase approach considers three main physical volumes: a biofilm matrix, the water volume enclosed by it and the flowing suspension through the SSF surrounding the biofilm. The only input data needed are the inflow concentrations and rate, light irradiation and temperature. The model includes a Cahn–Hilliard-type equation for the biofilm velocity in the supernatant water. A numerical scheme that preserves physical properties of the solution is also presented. This is achieved with a time-splitting scheme using an upwind scheme with an adaptive time-step size that ensures positivity of the solution under a CFL condition and a reasonable assumption on the separate numerical solution of the Cahn–Hilliard-type equation. The latter is solved numerically with a semi-implicit finite-difference scheme using a mixed form of the equation. The flexible model framework with its numerical scheme allows for including a variety of ecological models, making the investigation and development of increasingly complex simulations possible. Simulations with different model parameters are shown to be qualitatively consistent with the literature.

High-order oscillation-eliminating Hermite WENO method for hyperbolic conservation laws

This paper proposes high-order accurate, oscillation-eliminating, Hermite weighted essentially non-oscillatory (OE-HWENO) finite volume schemes for hyperbolic conservation laws, motivated by the oscillation-eliminating (OE) discontinuous Galerkin schemes recently proposed in [M. Peng, Z. Sun, and K. Wu, 2024 [30]]. The OE-HWENO schemes incorporate an OE procedure after each Runge–Kutta stage, by dampening the first-order moments of the HWENO solution to suppress spurious oscillations without any problem-dependent parameter. The OE procedure acts as a moment filter and is derived from the solution operator of a novel damping equation, which is exactly solved without any discretization. Thanks to this distinctive feature, the OE-HWENO method remains stable with a normal CFL number, even for strong shocks resulting in highly stiff damping terms. To ensure the essentially non-oscillatory property of the OE-HWENO method across problems with varying scales and wave speeds, we design a scale-invariant and evolution-invariant damping equation and propose a generic dimensionless transformation for HWENO reconstruction. The OE-HWENO method offers several notable advantages over existing HWENO methods. First, the OE procedure is highly efficient and straightforward to implement, requiring only simple multiplication of first-order moments by a damping factor. Furthermore, we rigorously prove that the OE procedure maintains the high-order accuracy and local compactness of the original HWENO schemes and demonstrate that it does not compromise the spectral properties via the approximate dispersion relation for smooth solutions. Notably, the proposed OE procedure is non-intrusive, enabling seamless integration as an independent module into existing HWENO codes. Finally, we rigorously analyze the bound-preserving (BP) property of the OE-HWENO method using the optimal cell average decomposition approach [S. Cui, S. Ding, and K. Wu, 2024 [8]], which relaxes the theoretical BP constraint for time step-size and reduces the number of decomposition points, thereby further enhancing efficiency. Extensive benchmarks validate the accuracy, efficiency, high resolution, and robustness of the OE-HWENO method.

An efficient eigenvalue bounding method: CFL condition revisited

Direct and large-eddy simulations of turbulence are often solved using explicit temporal schemes. However, this imposes very small time-steps because the eigenvalues of the (linearized) dynamical system, re-scaled by the time-step, must lie inside the stability region. In practice, fast and accurate estimations of the spectral radii of both the discrete convective and diffusive terms are therefore needed. This is virtually always done using the so-called CFL condition. On the other hand, the large heterogeneity and complexity of modern supercomputing systems are nowadays hindering the efficient cross-platform portability of CFD codes. In this regard, our leitmotiv reads: relying on a minimal set of (algebraic) kernels is crucial for code portability and maintenance! In this context, this work focuses on the computation of eigenbounds for the above-mentioned convective and diffusive matrices which are needed to determine the time-step à la CFL. To do so, a new inexpensive method, that does not require to re-construct these time-dependent matrices, is proposed and tested. It just relies on a sparse-matrix vector product where only vectors change on time. Hence, both implementation in existing codes and cross-platform portability are straightforward. The effectiveness and robustness of the method are demonstrated for different test cases on both structured Cartesian and unstructured meshes. Finally, the method is combined with a self-adaptive temporal scheme, leading to significantly larger time-steps compared with other more conventional CFL-based approaches.

Mass, momentum and energy preserving FEEC and broken-FEEC schemes for the incompressible Navier–Stokes equations

Abstract In this article we propose two finite-element schemes for the Navier–Stokes equations, based on a reformulation that involves differential operators from the de Rham sequence and an advection operator with explicit skew-symmetry in weak form. Our first scheme is obtained by discretizing this formulation with conforming FEEC (Finite Element Exterior Calculus) spaces: it preserves the point-wise divergence free constraint of the velocity, its total momentum and its energy, in addition to being pressure robust. Following the broken-FEEC approach, our second scheme uses fully discontinuous spaces and local conforming projections to define the discrete differential operators. It preserves the same invariants up to a dissipation of energy to stabilize numerical discontinuities. For both schemes we use a middle point time discretization that preserve these invariants at the fully discrete level and we analyze its well-posedness in terms of a CFL condition. While our theoretical results hold for general finite elements preserving the de Rham structure, we consider one application to tensor-product spline spaces. Specifically, we conduct several numerical test cases to verify the high order accuracy of the resulting numerical methods, as well as their ability to handle general boundary conditions.

A depth‐integrated SPH framework for slow landslides

AbstractSlow and very slow landslides can cause severe economic damage to structures. Due to their velocity of propagation, it is possible to take action such as programmed maintenance or evacuation of affected zones. Modeling is an important tool that allows scientists, engineers, and geologists to better understand their causes and predict their propagation. There are many available models of different complexities which can be used for this purpose, ranging from very simple infinite landslide models which can be implemented in spreadsheets to fully coupled 3D models. This approach is expensive because of the time span in which the problems are studied (sometimes years), simpler methods such as depth‐integrated models could provide a good compromise between accuracy and cost. However, there, the time step limitation due to CFL condition (which states that the time step has to be slower than the ratio between the node spacing and the physical velocity of the waves results in time increments which are of the order of one‐10th of a second on many occasions. This paper extends a technique that has been used in the past to glacier evolution problems using finite differences or elements to SPH depth‐integrated models for landslide propagation. The approach is based on assuming that (i) the flow is shallow, (ii) the rheological behavior determining the velocity of propagation is viscoplastic, and (iii) accelerations can be neglected. In this case, the model changes from hyperbolic to parabolic, with a time increment much larger than that of classic hyperbolic formulations.

Влияние числа Куранта на результаты численного моделирования распространения сигналов в недиспергирующих однородных средах

The paper is devoted to the study of the connection between the numerical dispersion arising in FDTD modeling of electromagnetic signal propagation in nondispersive homogeneous media optically different from vacuum and the Courant number in the 2D case. The main results are formulated in the form of four statements, as well as a number of corollaries and remarks that determine the nature of the numerical dispersion, the optimal value of the Courant number and the limitations of the method. It is proved that the optimal choice of the Courant number eliminates the numerical dispersion and extends the capabilities of the developed numerical algorithm to media, which refractive index is lesser than refractive index of vacuum, as well as media with negative refraction.

Two Finite Element Approaches for the Porous Medium Equation That Are Positivity Preserving and Energy Stable

In this work, we present the construction of two distinct finite element approaches to solve the porous medium equation (PME). In the first approach, we transform the PME to a log-density variable formulation and construct a continuous Galerkin method. In the second approach, we introduce additional potential and velocity variables to rewrite the PME into a system of equations, for which we construct a mixed finite element method. Both approaches are first-order accurate, mass conserving, and proved to be unconditionally energy stable for their respective energies. The mixed approach is shown to preserve positivity under a CFL condition, while a much stronger property of unconditional bound preservation is proved for the log-density approach. A novel feature of our schemes is that they can handle compactly supported initial data without the need for any perturbation techniques. Furthermore, the log-density method can handle unstructured grids in any number of dimensions, while the mixed method can handle unstructured grids in two dimensions. We present results from several numerical experiments to demonstrate these properties.

Evaluation on different volume of fluid methods in unstructured solver under the optimized condition

We compared the accuracy of volume of fluid (VOF) methods in unstructured solvers using the following five different methods: 1 - the algebraically compressive VOF method, 2 – simple coupled VOF method with Level Set (S-CLSVOF) method, 3 - interface-compressing VOF method incorporated with Laplacian filter (VOFL), 4 - isoAdvector method, and 5 - isoAdvector method incorporated with Laplacian filter (isoAdvectorL) by incorporating them into OpenFOAM®, an open-source software. To evaluate these methods under proper conditions, we compared the calculation accuracy using the optimized parameters, which are explored by Bayesian optimization. The test cases for advection accuracy of volume fraction and for imbalance of surface tension force in static multiphase fluid fields were considered. In this study, we found that the compression parameters and maximum Courant number should be adjusted to obtain high accuracy simulation according to the simulation condition in VOF and S-CLSVOF method. In VOFL and isoAdvectorL methods, the spurious current can be extremely reduced, which means that these methods are suitable for slow flow with higher Laplace number conditions.

Treatment for the Singularity Problem in the Phase Segregation

In the Eulerian–Eulerian two-fluid model, comprising air and droplets, singularities often arise as the phase fraction of dispersed droplets approaches zero. The nonconservative momentum equations can address this issue, but it will lead to inaccurate solutions when encountering discontinuities. To tackle the singularity problem and ensure accurate solutions for the droplet field, a novel adaptive form of momentum equations is proposed. The adaptive form defaults to the conservative form and smoothly transitions to the phase-intensive form when the phase fraction falls below a critical value, specifically 1/1000 of the inlet. The LU-SGS algorithm, enhanced with the inner iteration step method, is employed for numerical solutions of the droplet field. Comparative analysis reveals that the discrepancy between the results obtained from the adaptive form and the conservative form is negligible, under 0.5‰. Notably, with larger Courant numbers, the conservative form diverges from waterless regions despite imposing the minimum phase fraction constraint. Conversely, the adaptive form consistently maintains a Courant number of five times or more, potentially exceeding a hundredfold in complex flows. Although the adaptive form moderately increases complexity and computational time for a single time step, the improvement of the maximum Courant number will bring significant computational efficiency advantages. Moreover, the adaptive form holds promise for application in various other multiphase flow scenarios.

An unconditionally stable scheme for the immersed boundary method with application in cardiac mechanics

The stability of the immersed boundary (IB) method is a challenge in simulating fluid–structure interaction problems, where time step constraints are significantly stricter than in pure fluid simulations. We propose a novel unconditionally stable scheme for the immersed boundary finite element (IBFE) method. The structure is handled implicitly and characterized by strain energy functions, rather than being modeled as fibers or membranes. Through energy estimate, we prove the unconditional stability of the fully discrete approximation in the absence of the convective term. In real simulations of cardiac mechanics problems, the time step is much larger, only limited by the Courant–Friedrichs–Lewy condition of the fluid. The novelty of this work lies in the combination of dual interpolation and distribution operators, the Jacobian-free Newton–Krylov method for solving nonlinear algebraic systems, and the semi-Lagrangian method for handling the convective term. To validate the effectiveness and accuracy of our approach, we present various benchmarks and conduct a quasi-static simulation of a three-dimensional real left ventricular model. We have shown that the numerical stability of our scheme is very robust even with much larger time step compared to conventional explicit IB methods. Our work paves the way for further works on efficient solvers of the IBFE method.

Analysis of a positivity-preserving splitting scheme for some semilinear stochastic heat equations

We construct a positivity-preserving Lie–Trotter splitting scheme with finite difference discretization in space for approximating the solutions to a class of semilinear stochastic heat equations with multiplicative space-time white noise. We prove that this explicit numerical scheme converges in the mean-square sense, with rate 1/4 in time and rate 1/2 in space, under appropriate CFL conditions. Numerical experiments illustrate the superiority of the proposed numerical scheme compared with standard numerical methods which do not preserve positivity.

Space-Time FEM for the Vectorial Wave Equation under Consideration of Ohm’s Law

Abstract The ability to deal with complex geometries and to go to higher orders is the main advantage of space-time finite element methods. Therefore, we want to develop a solid background from which we can construct appropriate space-time methods. In this paper, we will treat time as another space direction, which is the main idea of space-time methods. First, we will briefly discuss how exactly the vectorial wave equation is derived from Maxwell’s equations in a space-time structure, taking into account Ohm’s law. Then we will derive a space-time variational formulation for the vectorial wave equation using different trial and test spaces. This paper has two main goals. First, we prove unique solvability for the resulting Galerkin–Petrov variational formulation. Second, we analyze the discrete equivalent of the equation in a tensor product and show conditional stability, i.e., under a CFL condition. Understanding the vectorial wave equation and the corresponding space-time finite element methods is crucial for improving the existing theory of Maxwell’s equations and paves the way to computations of more complicated electromagnetic problems.

Quasi-Newton positivity-preserving scheme for the Spalart-Allmaras turbulence model using unstructured grids

A new implicit method for solving the Spalart-Allmaras (SA) turbulence model is proposed using a formal second-order upwind scheme on unstructured grids for steady-state flows. The mean-flow set of equations and the SA model equation are solved in a loosely coupled manner. While the Jacobian of the mean-flow equations is derived from a low-order approximation of the residual (i.e. the defect-correction approach), the SA Jacobian is obtained by modifying a direct linearization of the second-order residual. The modified Jacobian is designed to form an M-matrix without the artificial time derivative commonly used for steady-state problems. Namely, no time step is involved in solving the SA model equation. Thanks to a special design of the M-matrix, the positivity of the SA model working variable and the convergence of the SA model linearized problem (fixed-point iterative problem) is guaranteed. To further enhance the overall stability of the flow solver, the Runge-Kutta implicit algorithm is employed, allowing CFL numbers as high as one to two thousand for the mean-flow equations. A second-order upwind scheme is achieved with a least-square method using a limiter to suppress oscillation in the solution. However, second-order upwind scheme discretization of turbulence models may lack robustness, especially when using unstructured grids. The limiter can hinder the convergence of the non-linear iterations, especially of the turbulence model. A square root transformation is applied to the dependent variable of a baseline SA model to improve the non-linear convergence characteristics. Three variants of the SA model are implemented and tested together with two limiters. Numerical simulations of a broad range of aerodynamic applications are conducted. The transformed variant exhibits the least sensitivity to the limiter type.