Courant Number Research Articles

A stable FE method for the space-time solution of the Cahn-Hilliard equation

In its application to the modeling of a mineral separation process, we propose the numerical analysis of the Cahn-Hilliard equation by employing space-time discretizations of the automatic variationally stable finite element (AVS-FE) method. The AVS-FE method is a Petrov-Galerkin method which employs the concept of optimal discontinuous test functions of the discontinuous Petrov-Galerkin (DPG) method by Demkowicz and Gopalakrishnan. The trial space, however, consists of globally continuous Hilbert spaces such as H1(Ω) and H(div,Ω). Hence, the AVS-FE approximations employ classical C0 or Raviart-Thomas FE basis functions. The optimal test functions guarantee the numerical stability of the AVS-FE method and lead to discrete systems that are symmetric and positive definite. Hence, the AVS-FE method can solve the Cahn-Hilliard equation in both space and time without a restrictive CFL condition to dictate the space-time element size. We present multiple numerical verifications of both stationary and transient problems. The verifications show optimal rates of convergence in L2(Ω) and H1(Ω) norms. Results for mesh adaptive refinements in both space and time using a built-in error estimator of the AVS-FE method are also presented.

Exact spatial and temporal balance of energy exchanges within a horizontally explicit/vertically implicit non-hydrostatic atmosphere

A new horizontally explicit/vertically implicit (HEVI) time splitting scheme for atmospheric modelling is introduced, for which the horizontal divergence terms are applied within the implicit vertical substep. The new HEVI scheme is implemented in conjunction with a mixed mimetic spectral element spatial discretisation and semi-implicit vertical time stepping scheme that both preserve the skew-symmetric structure of the non-canonical Hamiltonian form of the equations of motion. Within this context the new HEVI scheme allows for the exact balance of all energetic exchanges in space and time. However since the choice of horizontal fluxes for which this balance is satisfied is not consistent with the horizontal velocity at the end of the time level the scheme still admits a temporal energy conservation error. Linearised eigenvalue analysis shows that similar to a fully implicit method, the new HEVI scheme is neutrally stable for all buoyancy modes, and unlike a second order trapezoidal HEVI scheme is stable for all acoustic modes below a certain horizontal CFL number. The scheme is validated against standard test cases for both planetary and non-hydrostatic regimes. For the planetary scale baroclinic instability test case, the new formulation exhibits a secondary oscillation in the potential to kinetic energy power exchanges, with a temporal frequency approximately four times that exhibited by a horizontally third order, vertically second order trapezoidal scheme. For the non-hydrostatic test case, the vertical upwinding of the potential temperature diagnostic equation is shown to reduce spurious oscillations without altering the energetics of the solution, since this upwinding is performed in an energetically consistent manner. For this test case, which is configured on an affine geometry, the exact balance of energy exchanges allows the model to run stably without any form of dissipation.

A stable discontinuous Galerkin method based on high‐order dual splitting scheme without additional stabilization term for incompressible flows

AbstractIn this work, we consider equal‐order discontinuous Galerkin (DG) solver for incompressible Navier–Stokes equations based on high‐order dual splitting scheme. In order to stay stable, the time step size of this method has been reported that is strictly limited. The upper bound of time step size is restricted by Courant–Friedrichs–Lewy (CFL) condition (Hesthaven and Warburton, 2007) and lower bound is required to be larger than the critical value which depends on Reynolds number and spatial resolution (Ferrer et al., 2014). For high‐Reynolds‐number flow problems, if the spatial resolution is low, the critical value may be larger than CFL condition, then instability will occur for any time step size. Therefore, sufficiently high spatial resolution is indispensable in order to maintain stability, which increases the computational cost. To overcome these difficulties and develop a robust solver for high‐Reynolds‐number flow problem, it is necessary to further study the instability problem at small time steps. We numerically investigate the effect of the pressure gradient term in projection step and the velocity divergence term in pressure Poisson equation on the stability for small time step size, respectively, and conclude that the DG formulation of the pressure gradient term has a more significant effect on the stability of the scheme than that of the velocity divergence term. Integration by parts of these terms is essential in order to improve the stability of the scheme. Based on this discretization format, an appropriate penalty parameter for pressure Poisson equation is utilized so as to provide the scheme with an inf‐sup stabilization. Moreover, the lid‐driven cavity flow is considered to verify that this numerical algorithm enhances the stability without additional stabilization term at small time step size and high‐Reynolds number for equal‐order polynomial approximations.

Computational multiscale methods for quasi-gas dynamic equations

In this paper, we consider the quasi-gas-dynamic (QGD) model in a multiscale environment. The model equations can be regarded as a hyperbolic regularization and are derived from kinetic equations. So far, the research on QGD models has been focused on problems with constant coefficients. In this paper, we investigate the QGD model in multiscale media, which can be used in porous media applications. This multiscale problem is interesting from a multiscale methodology point of view as the model problem has a hyperbolic multiscale term, and designing multiscale methods for hyperbolic equations is challenging. In the paper, we apply the constraint energy minimizing generalized multiscale finite element method (CEM-GMsFEM) combined with the central difference scheme in time to solve this problem. The CEM-GMsFEM provides a flexible and systematical framework to construct crucial multiscale basis functions for approximating the solution to the problem with reduced computational cost. With this approach of spatial discretization, we establish the stability of the fully discretized scheme, based on the coarse grid, under a coarse-scale CFL condition. Complete convergence analysis of the proposed method is presented. Numerical results are provided to illustrate and verify the theoretical findings.

Bandpass approximate Crank-Nicolson implementation for anisotropic gyrotropic plasma open region simulation

For overcoming the Courant-Friedrich-Levy (CFL) stable condition and maintaining the stability of the algorithm, unconditionally stable implementation is proposed for the simulation of anisotropic gyrotropic plasma in bandpass open region problems. To be more accurately, the proposed implementation is based on the complex envelope method, higher order concept, perfectly matched layer formulation and Crank-Nicolson Direct-Splitting procedure which shows advantages over them including accuracy, efficiency and absorption. The gyrotropic plasma with unique anisotropic characteristic can be simulated by the modified approach. For the further illustration of effectiveness, several numerical examples are introduced. Through the results and comparison with previous work, it can be concluded that the entire performance which includes accuracy, efficiency and absorption can be significantly improved. In addition, the mesh size can be no longer limited by the CFL condition rather than the computational accuracy which shows considerable potential in numerical simulation of fine structures. Meanwhile, the time step can be obtained according to the maximum frequency rather than bandwidth resulting in the exponential decrement of running time.

Simulation of groundwater flow using meshfree collocation method with Cubic Spline function

Groundwater is a long term reservoir of the natural water cycle. India is the largest user of groundwater in the world. There are many areas in the country facing scarcity of water primarily due to unplanned groundwater development resulting in fall of water levels. If measures to ensure sustainable groundwater usage are not taken consequences may include a collapse of shortages of potable water, conflict and suffering. Therefore it is necessary to identify the groundwater resources and predict the groundwater flows at different locations in the aquifers. Numerical methods can simulate the groundwater flows for the delineated aquifers. In this study, the simulation of groundwater flow is carried out using Meshfree Collocation Method (CM) with Cubic Spline (CS) function and moving support domains. The results are compared with the analytical solution and found satisfactory. The sensitivity analysis is carried out for Courant number, simulation scheme, spatial and temporal variation.

Analytical ‘steady-state’-based derivation and clarification of the courant-friedrichs-lewy condition for pipe flow

This article addresses the problem of choosing the optimal discretization grid for emulating fluid flow through a pipeline. The aggregated basic flow model is linearized near the operating point obtained from the steady state analytic solution of the differential equations under consideration. Based on this model, the relationship between the Courant number (μ) and the stability margin is examined. The numerically set coefficient μopt, ensuring the maximum margin of stability, is analyzed in terms of the physical and technological parameters of the flow. As a result of this analysis, a specific formula is obtained based on parameters describing the mechanics (geometry and physics) of the flow through the pipeline, which leads to the optimal value of the Courant number, separately for smooth and rough pumping conditions. A more detailed analysis of the distribution of the optimal μ coefficient in relation to the parameters of the pipeline flow mechanics shows four cases to consider when determining the coefficient μopt. Surprisingly, in three cases, the CFL condition is insufficient, which is expressed in the form of the proposed procedure for choosing the optimal value of μ. The final dichotomous model is derived from the Monte Carlo simulation results in which the effect of each parameter on the optimal Courant number is estimated and consolidated. Taking into account the recognized general laws of physics and using numerical methods and mathematical analysis, simple and useful analytical relationships describing the flow process are obtained. In addition, computer simulations are performed to verify the correctness of the proposed procedure, as well as a number of other considerations related to the modeling of fluid flow in transport pipelines.

An Analysis on the Finite Volume Schemes and the Discrete Lyapunov Inequalities for the Chemotaxis System

We analyze two finite volume schemes, linear and nonlinear, for the chemotaxis system in two-dimensional domain, which preserve the mass conservation and positivity without the CFL condition. For the nonlinear scheme, the well-posedness is proved by using Brouwer’s fixed point theory, and we show the convergence of the Picard iteration. We also investigate two discrete Lyapunov functionals, the asymptotic stability of equilibrium and the local stability. Moreover, we apply the discrete semi-group theory to error analysis and obtain the convergence rate $$O(\tau +h)$$ in $$L^p$$ norm. The theoretical results are confirmed by numerical experiments.

An Implicit Finite Volume Scheme to Solve the Time-dependent Radiation Transport Equation Based on Discrete Ordinates

Abstract We describe a new algorithm to implicitly solve the time-dependent, frequency-integrated radiation transport (RT) equation, which is coupled to an explicit solver for equations of magnetohydrodynamics (MHD) using Athena++. The radiation field is represented by specific intensities along discrete rays, which are evolved using a conservative finite volume approach for both Cartesian and curvilinear coordinate systems. All terms for spatial transport of photons and interactions between gas and radiation are calculated implicitly together. An efficient Jacobi-like iteration scheme is used to solve the implicit equations. This removes any time-step constraint due to the speed of light in RT. We evolve the specific intensities in the lab frame to simplify the transport step. The lab frame specific intensities are transformed to the comoving frame via Lorentz transformation when the source term is calculated. Therefore, the scheme does not need any expansion in terms of v/c. The radiation energy and momentum source terms for the gas are calculated via direct quadrature in the angular space. The time step for the whole scheme is determined by the normal Courant–Friedrichs–Lewy condition in the MHD module. We provide a variety of test problems for this algorithm, including both optically thick and thin regimes, and for both gas and radiation pressure-dominated flows to demonstrate its accuracy and efficiency.

Semi-implicit front capturing schemes for the degenerate nonlinear radiative diffusion equation

Nonlinear radiative diffusion equation is coupled with fluid equations or magnetic hydrodynamical equations in radiation hydrodynamics simulations. The time scales of radiation and fluid can be different. To use the same time step in the coupled system, one needs to design numerical schemes for the nonlinear radiative diffusion equation that are stable for large time steps. Due to the nonlinearity and degeneracy, fully implicit schemes may exhibit time step constraints. We propose and test two different semi-implicit schemes in this paper. One allows for a hyperbolic CFL condition, and only one linear diffusion system and a scalar forth order polynomial on each space grid have to be solved implicitly in each time step. The other has a much larger CFL number, but a nonlinear system has to be solved in each time step. Though the computational cost of the second scheme in each time step is higher, the overall computational cost can be lower. Extensive numerical examples for both smooth and compactly supported solutions are presented to show the accuracy and stability of the schemes.

A finite-volume scheme for gradient-flow equations with non-homogeneous diffusion

We develop a first- and second-order finite-volume scheme to solve gradient flow equations with non-homogeneous properties, obtained in the framework of dynamical-density functional theory. The scheme takes advantage of an upwind approach for the space discretization to ensure positivity of the density under a CFL condition and decay of the discrete free energy. Our computational approach is used to study several one- and two-dimensional systems, with a general free-energy functional accounting for external fields and inter-particle potentials, and placed in non-homogeneous thermal baths characterized by anisotropic, space-dependent and time-dependent properties.

A High-Order Maximum-Principle-Satisfying Discontinuous Galerkin Method for the Level Set Problem

Level set (LS) method is a widely used interface capturing method. In the simulations of incompressible two-phase flows, in order to avoid discontinuities at interfaces, the LS function is usually taken as a smeared-out Heaviside function bounded on [0, 1] and advected by a given velocity field $$\mathbf {u}$$ obtained from the solution of the incompressible Navier-Stokes equations. In the incompressible limit $$\nabla \cdot \mathbf {u}=0$$ , the advection equation for the LS function can be written and discretized in conservative form. However, due to numerical errors, the resulting velocity field is in general not divergence free which leads to the solution of the advection equation in conservative form does not satisfy the maximum principle. To overcome this issue, in this work, we develop a high-order discontinuous Galerkin (DG) method to directly solve the advection equation for the LS function in non-conservative form. Moreover, we prove that by applying a linear scaling limiter, the proposed method together with a strong stability preserving (SSP) time discretization scheme can satisfy the strict maximum principle under a suitable CFL condition. Numerical simulations of several well-known benchmark problems, including the application to incompressible two-phase flows, are presented to demonstrate the high-order accuracy and maximum-principle-satisfying property of the proposed method.

Partitioned water hammer modeling using the block Gauss–Seidel algorithm

The monolithic classical and extended fluid-structural interaction water hammer models are based on the simplified ring hypothesis to characterize the structural radial response of the pipe under pressure. Recent work has shown that this approach predicts an unrealistic discontinuity in the pipe wall in the neighborhood of the pressure wave front. In order to remedy this limitation, the present study adopts a partitioned approach which couples a continuous shell-based finite element model with a water hammer fluid model based on the method of characteristics, within the context of the partitioned block Gauss–Seidel iterative algorithm. In order to demonstrate the validity of the block Gauss–Seidel implementation, the classical and extended water hammer models, normally solved using a monolithic approach, are solved under the present partitioned model and shown to predict pressure histories identical to those based on the monolithic solution. Comparisons are then conducted for the pressure histories as predicted by the classical, extended, and shell-based approaches. The study shows that the accuracy is optimized for integer values of the Courant number. By adopting either an optimal constant relaxation factor or the Aitken relaxation factor, the number of iterations needed for convergence was significantly reduced. The accuracy and computational efficiency are shown to highly depend on the number of subdivisions in the fluid model. In contrast, the number of subdivisions in the structural model, while influencing the computational efficiency, is shown not to influence the accuracy of the predictions. When applying the partitioned approach to the classical and shell-based models, the predicted pressure histories are found to be independent of the specified tolerance. In contrast, when applied to the extended model, the chosen tolerance has implications on the stability and accuracy of the solution.

Two Dimensional Analysis of Hybrid Spectral/Finite Difference Schemes for Linearized Compressible Navier\u2013Stokes Equations

This study seeks to compare different combinations of spatial dicretization methods under a coupled spatial temporal framework in two dimensional wavenumber space. The aim is to understand the effect of dispersion and dissipation on both the convection and diffusion terms found in the two dimensional linearized compressible Navier–Stokes Equations (LCNSE) when a hybrid finite difference/Fourier spectral scheme is used in the x and y directions. In two dimensional wavespace, the spectral resolution becomes a function of both the wavenumber and the wave propagation angle, the orientation of the wave front with respect to the grid. At sufficiently low CFL number where temporal discretization effects can be neglected, we show that a hybrid finite difference/Fourier spectral schemes is more accurate than a full finite difference method for the two dimensional advection equation, but that this is not so in the case of the LCNSE. Group velocities, phase velocities as well as numerical amplification factor were used to quantify the numerical anisotropy of the dispersion and dissipation properties. Unlike the advection equation, the dispersion relation representing the acoustic modes of the LCNSE contains an acoustic terms in addition to its advection and viscous terms. This makes the group velocity in each spatial direction a function of the wavenumber in both spatial directions. This can lead to conditions for which a hybrid Fourier spectral/finite difference method can become less or more accurate than a full finite difference method. To better understand the comparison of the dispersion properties between a hybrid and full FD scheme, the integrated sum of the error between the numerical group velocity V^{*}_{grp,full} and the exact solution across all wavenumbers for a range of wave propagation angle is examined. In the comparison between a hybrid and full FD discretization schemes, the fourth order central (CDS4), fourth order dispersion relation preserving (DRP4) and sixth order central compact (CCOM6) schemes share the same characteristics. At low wave propagation angle, the integrated errors of the full FD and hybrid discretization schemes remain the same. At intermediate wave propagation angle, the integrated error of the full FD schemes become smaller than that of the hybrid scheme. At large wave propagation angle, the integrated error of the full FD schemes diverges while the integrated error of the hybrid discretization schemes converge to zero. At high reduced wavenumber and sufficiently low CFL number where temporal discretization error can be neglected, it was found that the numerical dissipation of the viscous term based on the CDS4, DRP4, CCOM6 and isotropy optimized CDS4 schemes (hbox {CDS4}_{{opt}}) schemes was lower than the actual physical dissipation, which is only a function of the cell Reynolds number. The wave propagation angle at which the numerical dissipation of the viscous term approaches its maximum occurs at pi /4 for the CDS4, DRP4, CCOM6 and hbox {CDS4}_{{opt}} schemes.

On explicit discontinuous Galerkin methods for conservation laws

The concept of projection applied to explicit discontinuous Galerkin (DG) schemes is investigated. The two explicit DG methods in this study are based on a ‘predictor-corrector’ formulation, the first introduced by Lörcher, Gassner, and Munz called space–time expansion discontinuous Galerkin or STE-DG scheme (Lörcher et al. 2007, Gassner et al. 2008), and the second, introduced independently by the author (Huynh 2006, 2013) called the upwind moment scheme. The predictor step of the two methods is essentially identical using a Cauchy-Kovalevsky procedure, which involves no interaction of the data among neighboring cells. The corrector step also shares the same space-time integration formulation and is where interaction takes place; the two methods differ, however, in the evaluation of the projections in the space-time volume integrals. The STE-DG scheme evaluates these in a straightforward manner, whereas the moment scheme employs a successive procedure with each moment update uses the results by the lower-order updates. The trade-off is that the moment scheme has the disadvantage of a more elaborate corrector step and the significant advantage of a CFL (Courant-Friedrichs-Lewy) condition of 1 for all (degree) p and accuracy order of 2p+1 (super accuracy property) for one-dimensional (1D) advection. In contrast, the STE-DG method is accurate to only the expected order of p+1 and has a more restrictive CFL condition. Due to the predictor-corrector formulation that does not involve methods of characteristics, these schemes extend easily to systems of equations in multiple dimensions. Concerning the case of two spatial dimensions (2D), for an advection using a Cartesian grid, when the flow does not align with the axes, especially when it is along a diagonal direction, the CFL conditions for the moment schemes also become restrictive and need improvement as will be shown by Fourier (von Neumann) analyses.

Krylov solvers in a vertical‐slice version of the semi‐implicit semi‐Lagrangian AROME model

AbstractTo circumvent the scalability problem due to global communication involved in spectral transforms, a vertical‐slice version of the dynamical core where all calculations are performed in grid‐point space has been built for AROME, Météo‐France's operational limited‐area model. It is shown in an idealized but nevertheless physically relevant framework that, despite this major change, it is possible to keep the other main characteristics of the model (constant‐coefficient semi‐implicit scheme, semi‐Lagrangian transport scheme, A‐grid, mass‐based coordinate, etc.). A Krylov solver is used to solve the implicit problem. Using the solution given by the spectral model as a reference in terms of quality and required accuracy in an operational context, the chosen parameters of the Krylov solver are carefully tuned to maximize its convergence speed. The Krylov solver consists of several applications of a sparse operator, whose stencil is similar to that of the operator applied during the small time step of split‐explicit schemes traditionally used in HEVI (Horizontally Explicit/Vertically Implicit) models. HEVI models are considered particularly scalable in the current parallelization paradigm and are used as a reference for scalability in this study. Results show that a fast convergence can be reached, such that only a few applications of the sparse operator are required. This suggests an improvement in scalability compared to the spectral version. Experiments have been formulated regardless of the future strategy considered for temporal and spatial resolutions, i.e., increasing spatial resolution while either maintaining a high Courant number or reducing it.

Cahn-Hilliard Navier-Stokes simulations for marine free-surface flows

The paper is devoted to the simulation of maritime two-phase flows of air and water. Emphasis is put on an extension of the classical Volume-of-Fluid (VoF) method by a diffusive contribution derived from a Cahn-Hilliard (CH) model and its benefits for simulating immiscible, incompressible two-phase flows. Such flows are predominantly simulated with implicit VoF schemes, which mostly employ heuristic downwind-biased approximations for the concentration transport to mimic a sharp interface. This strategy introduces a severe time step restriction and requires pseudo time-stepping of steady flows. Our overall goal is a sound description of the free-surface region that alleviates artificial time-step restrictions, supports an efficient and robust upwind-based approximation framework, and inherently includes surface tension effects when needed. The Cahn-Hilliard Navier-Stokes (CH-NS) system is verified for an analytical Couette-flow example and the bubble formation under the influence of surface tension forces. 2D validation examples are concerned with laminar standing waves reaching from gravity to capillary scale as well as a submerged hydrofoil flow. The final application refers to the 3D flow around an experimentally investigated container vessel at fixed floatation for Re = 1.4 × 107 and Fn = 0.26. Results are compared with data obtained from VoF approaches, supplemented by analytical solutions and measurements. The study indicates the superior efficiency, resharpening capability, and wider predictive realm of the CH-based extension for free-surface flows with a confined spatial range of interface Courant numbers.

Direct-forcing immersed-boundary method: A simple correction preventing boundary slip error

The boundary slip error resulting from the interpolation/spreading non-reciprocity of the direct-forcing immersed-boundary method is analyzed based on a simple and generic theoretical framework. In explicit implementations, the slip error scales with the Courant number, as predicted by the analysis and confirmed by lattice-Boltzmann simulation results. Using an analytical approximation of the non-reciprocity error, the immersed-boundary force can be corrected in order to prevent boundary slip and flow penetration. This a priori correction leads to a major improvement of the no-slip condition while avoiding any additional computational time or implementation effort.

Monotonicity-preserving splitting schemes for solving balance laws

In this paper, some monotonicity-preserving (MP) and positivity-preserving (PP) splitting methods for solving the balance laws of the reaction and diffusion source terms are investigated. To capture the solution with high accuracy and resolution, the original equation with reaction source termis separated through the splitting method into two sub-problems including the homogeneous conservation law and a simple ordinary differential equation (ODE). The resulting splitting methods preserve monotonicity and positivity property for a normal CFL condition. A trenchant numerical analysis made it clear that the computing time of the proposed methods decreases when the so-called MP process for the homogeneous conservation law is imposed. Moreover, the proposed methods are successful in recapturing the solution of the problem with high-resolution in the case of both smooth and non-smooth initial profiles. To show the efficiency of proposed methods and to verify the order of convergence and capability of these methods, several numerical experiments are performed through some prototype examples.

MHD-kinetic hybrid code based on structure-preserving finite elements with particles-in-cell

We present a STRUcture-Preserving HYbrid code - STRUPHY - for the simulation of magneto-hydrodynamic (MHD) waves interacting with a population of energetic particles far from thermal equilibrium (kinetic species). The implemented model features linear, ideal MHD equations in curved, three-dimensional space, coupled nonlinearly to the full-orbit Vlasov equations via a current coupling scheme. The algorithm is based on finite element exterior calculus (FEEC) for MHD and particle-in-cell (PIC) methods for the kinetic part; it provably conserves mass, energy, and the divergence-free magnetic field, irrespective of metric (= space curvature), mesh parameters and chosen order of the scheme. These properties enable reliable long-time simulations of energetic particle physics in complex geometries, covering the whole range of MHD waves. In STRUPHY, the finite element spaces are built from tensor products of univariate B-splines on the logical cuboid and can be made high-order by increasing the polynomial degree. Time-stepping is based on a skew-symmetric splitting with implicit sub-steps, mitigating CFL conditions from fast magneto-acoustic waves. High-order time splitting schemes can be used in this regard.