Asymptotic-Preserving Scheme for the M1-Maxwell System in the Quasi-Neutral Regime

This work deals with the numerical resolution of the M₁-Maxwell system in the quasi-neutral regime. In this regime the stiffness of the stability constraints of classical schemes causes huge calculation times. That is why we introduce a new stable numerical scheme consistent with the transitional and limit models. Such schemes are called Asymptotic-Preserving (AP) schemes in literature. This new scheme is able to handle the quasi-neutrality limit regime without any restrictions on time and space steps. This approach can be easily applied to angular moment models by using a moments extraction. Finally, two physically relevant numerical test cases are presented for the Asymptotic-Preserving scheme in different regimes. The first one corresponds to a regime where electromagnetic effects are predominant. The second one on the contrary shows the efficiency of the Asymptotic-Preserving scheme in the quasi-neutral regime. In the latter case the illustrative simulations are compared with kinetic and hydrodynamic numerical results.     

Fifth-Order A-WENO Path-Conservative Central-Upwind Scheme for Behavioral Non-Equilibrium Traffic Models

Non-equilibrium hyperbolic traffic models can be derived as continuum approximations of car-following models and in many cases the resulting continuum models are non-conservative. This leads to numerical difficulties, which seem to have discouraged further development of complex behavioral continuum models, which is a significant research need.In this paper, we develop a robust numerical scheme that solves hyperbolic traffic flow models based on their non-conservative form. We develop a fifth-order alternative weighted essentially non-oscillatory (A-WENO) finite-difference scheme based on the path-conservative central-upwind (PCCU) method for several non-equilibrium traffic flow models. In order to treat the non-conservative product terms, we use a path-conservative technique. To this end, we first apply the recently proposed second-order finite-volume PCCU scheme to the traffic flow models, and then extend this scheme to the fifth-order of accuracy via the finite-difference A-WENO framework. The designed schemes are applied to three different traffic flow models and tested on a number of challenging numerical examples. Both schemes produce quite accurate results though the resolution achieved by the fifth-order A-WENO scheme is higher. The proposed scheme in this paper sets the stage for developing more robust and complex continuum traffic flow models with respect to human psychological factors.     

A Modified Crank-Nicolson Numerical Scheme for the Flory-Huggins Cahn-Hilliard Model

In this paper we propose and analyze a second order accurate numerical scheme for the Cahn-Hilliard equation with logarithmic Flory Huggins energy potential. A modified Crank-Nicolson approximation is applied to the logarithmic nonlinear term, while the expansive term is updated by an explicit second order Adams-Bashforth extrapolation, and an alternate temporal stencil is used for the surface diffusion term. A nonlinear artificial regularization term is added in the numerical scheme, which ensures the positivity-preserving property, i.e., the numerical value of the phase variable is always between -1 and 1 at a point-wise level. Furthermore, an unconditional energy stability of the numerical scheme is derived, leveraging the special form of the logarithmic approximation term. In addition, an optimal rate convergence estimate is provided for the proposed numerical scheme, with the help of linearized stability analysis. A few numerical results, including both the constant-mobility and solution-dependent mobility flows, are presented to validate the robustness of the proposed numerical scheme.     

Full Discretization Scheme for the Dynamics of Elliptic Membrane Shell Model

In this study, the dynamics of elliptic membrane shell model has been proposed and discussed numerically for the first time. Firstly, we show that the solution of this model exists and is unique. Secondly, we consider spatial and time discretizations of the time-dependent elliptic membrane shell by finite element method and Newmark scheme, respectively. Then, the corresponding existence, uniqueness, stability, convergence and a priori error estimate are given. Finally, we present numerical results involving a portion of an ellipsoidal shell and a portion of a spherical shell to verify the efficiency and convergence of the numerical scheme.     

A Positivity-Preserving Second-Order BDF Scheme for the Cahn-Hilliard Equation with Variable Interfacial Parameters

We present and analyze a new second-order finite difference scheme for the Macromolecular Microsphere Composite hydrogel, Time-Dependent Ginzburg-Landau (MMC-TDGL) equation, a Cahn-Hilliard equation with Flory-Huggins-deGennes energy potential. This numerical scheme with unconditional energy stability is based on the Backward Differentiation Formula (BDF) method in time derivation combining with Douglas-Dupont regularization term. In addition, we present a pointwise bound of the numerical solution for the proposed scheme in the theoretical level. For the convergent analysis, we treat three nonlinear logarithmic terms as a whole and deal with all logarithmic terms directly by using the property that the nonlinear error inner product is always non-negative. Moreover, we present the detailed convergent analysis in $ℓ^∞$(0,$T$;$H_h^{-1}$)∩$ℓ^2$(0,$T$;$H_h^1$) norm. At last, we use the local Newton approximation and multigrid method to solve the nonlinear numerical scheme, and various numerical results are presented, including the numerical convergence test, positivity-preserving property test, spinodal decomposition, energy dissipation and mass conservation properties.     

GPU Accelerated Discontinuous Galerkin Methods for Shallow Water Equations

We discuss the development, verification, and performance of a GPU accelerated discontinuous Galerkin method for the solutions of two dimensional nonlinear shallow water equations. The shallow water equations are hyperbolic partial differential equations and are widely used in the simulation of tsunami wave propagations. Our algorithms are tailored to take advantage of the single instruction multiple data (SIMD) architecture of graphic processing units. The time integration is accelerated by local time stepping based on a multi-rate Adams-Bashforth scheme. A total variational bounded limiter is adopted for nonlinear stability of the numerical scheme. This limiter is coupled with a mass and momentum conserving positivity preserving limiter for the special treatment of a dry or partially wet element in the triangulation. Accuracy, robustness and performance are demonstrated with the aid of test cases. Furthermore, we developed a unified multi-threading model OCCA. The kernels expressed in OCCA model can be cross-compiled with multi-threading models OpenCL, CUDA, and OpenMP. We compare the performance of the OCCA kernels when cross-compiled with these models.     

A Third Order BDF Energy Stable Linear Scheme for the No-Slope-Selection Thin Film Model

In this paper we propose and analyze a (temporally) third order accurate backward differentiation formula (BDF) numerical scheme for the no-slope-selection (NSS) equation of the epitaxial thin film growth model, with Fourier pseudo-spectral discretization in space. The surface diffusion term is treated implicitly, while the nonlinear chemical potential is approximated by a third order explicit extrapolation formula for the sake of solvability. In addition, a third order accurate Douglas-Dupont regularization term, in the form of $−A∆t^2∆^2_N (u^{n+1}−u^n)$, is added in the numerical scheme. A careful energy stability estimate, combined with Fourier eigenvalue analysis, results in the energy stability in a modified version, and a theoretical justification of the coefficient $A$ becomes available. As a result of this energy stability analysis, a uniform in time bound of the numerical energy is obtained. And also, the optimal rate convergence analysis and error estimate are derived in details, in the $ℓ^∞(0,T;ℓ^2)∩ℓ^2(0,T;H^2_h)$ norm, with the help of a linearized estimate for the nonlinear error terms. Some numerical simulation results are presented to demonstrate the efficiency of the numerical scheme and the third order convergence. The long time simulation results for $ε = 0.02$ (up to $T = 3×10^5$) have indicated a logarithm law for the energy decay, as well as the power laws for growth of the surface roughness and the mound width. In particular, the power index for the surface roughness and the mound width growth, created by the third order numerical scheme, is more accurate than those produced by certain second order energy stable schemes in the existing literature.     

A Two-Phase Flow Simulation of Discrete-Fractured Media Using Mimetic Finite Difference Method

Various conceptual models exist for numerical simulation of fluid flow in fractured porous media, such as dual-porosity model and equivalent continuum model. As a promising model, the discrete-fracture model has been received more attention in the past decade. It can be used both as a stand-alone tool as well as for the evaluation of effective parameters for the continuum models. Various numerical methods have been applied to the discrete-fracture model, including control volume finite difference, Galerkin and mixed finite element methods. All these methods have inherent limitations in accuracy and applicabilities. In this work, we developed a new numerical scheme for the discrete-fracture model by using mimetic finite difference method. The proposed numerical model is applicable in arbitrary unstructured grid cells with full-tensor permeabilities. The matrix-fracture and fracture-fracture fluxes are calculated based on powerful features of the mimetic finite difference method, while the upstream finite volume scheme is used for the approximation of the saturation equation. Several numerical tests in 2D and 3D are carried out to demonstrate the efficiency and robustness of the proposed numerical model.     

Entropy Stable Scheme on Two-Dimensional Unstructured Grids for Euler Equations

We propose an entropy stable high-resolution finite volume scheme to approximate systems of two-dimensional symmetrizable conservation laws on unstructured grids. In particular we consider Euler equations governing compressible flows. The scheme is constructed using a combination of entropy conservative fluxes and entropy-stable numerical dissipation operators. High resolution is achieved based on a linear reconstruction procedure satisfying a suitable sign property that helps to maintain entropy stability. The proposed scheme is demonstrated to robustly approximate complex flow features by a series of benchmark numerical experiments.     

Numerical Solutions of Coupled Nonlinear Schrödinger Equations by Orthogonal Spline Collocation Method

In this paper, we present the use of the orthogonal spline collocation method for the semi-discretization scheme of the one-dimensional coupled nonlinear Schrödinger equations. This method uses the Hermite basis functions, by which physical quantities are approximated with their values and derivatives associated with Gaussian points. The convergence rate with order O(h⁴+τ²) and the stability of the scheme are proved. Conservation properties are shown in both theory and practice. Extensive numerical experiments are presented to validate the numerical study under consideration.     

Numerical Stability Analysis for a Stationary and Translating Droplet at Extremely Low Viscosity Values Using the Lattice Boltzmann Method Color-Gradient Multi-Component Model with Central Moments Formulation

Multicomponent models based on the Lattice Boltzmann Method (LBM) have clear advantages with respect to other approaches, such as good parallel performances and scalability and the automatic resolution of breakup and coalescence events. Multicomponent flow simulations are useful for a wide range of applications, yet many multicomponent models for LBM are limited in their numerical stability and therefore do not allow exploration of physically relevant low viscosity regimes. Here we perform a quantitative study and validations, varying parameters such as viscosity, droplet radius, domain size and acceleration for stationary and translating droplet simulations for the color-gradient method with central moments (CG-CM) formulation, as this method promises increased numerical stability with respect to the non-CM formulation. We focus on numerical stability and on the effect of decreasing grid-spacing, i.e. increasing resolution, in the extremely low viscosity regime for stationary droplet simulations. The effects of small- and large-scale anisotropy, due to grid-spacing and domain-size, respectively, are investigated for a stationary droplet. The effects on numerical stability of applying a uniform acceleration in one direction on the domain, i.e. on both the droplet and the ambient, is explored into the low viscosity regime, to probe the numerical stability of the method under dynamical conditions.     

Relaxation Schemes for the $M_1$ Model with Space-Dependent Flux: Application to Radiotherapy Dose Calculation

Because of stability constraints, most numerical schemes applied to hyperbolic systems of equations turn out to be costly when the flux term is multiplied by some very large scalar. This problem emerges with the $M_1$ system of equations in the field of radiotherapy when considering heterogeneous media with very disparate densities. Additionally, the flux term of the $M_1$ system is non-linear, and in order for the model to be well-posed the numerical solution needs to fulfill conditions called realizability. In this paper, we propose a numerical method that overcomes the stability constraint and preserves the realizability property. For this purpose, we relax the $M_1$ system to obtain a linear flux term. Then we extend the stencil of the difference quotient to obtain stability. The scheme is applied to a radiotherapy dose calculation example.     

A Well-Balanced Gas Kinetic Scheme for Navier-Stokes Equations with Gravitational Potential

The hydrostatic equilibrium state is the consequence of the exact balance between hydrostatic pressure and external force. Standard finite volume cannot keep this balance exactly due to their unbalanced truncation errors. In this study, we introduce an auxiliary variable which becomes constant at isothermal hydrostatic equilibria and propose a well-balanced gas kinetic scheme for the Navier-Stokes equations. Through reformulating the convection term and the force term via the auxiliary variable, zero numerical flux and zero numerical source term are enforced at the hydrostatic equilibrium state instead of the balance between hydrostatic pressure and external force. Several problems are tested to demonstrate the accuracy and the stability of the new scheme. The results confirm that, the new scheme can preserve the exact hydrostatic solution. The small perturbation riding on hydrostatic equilibria can be calculated accurately. More importantly, the new scheme is capable of simulating the process of converging towards hydrostatic equilibria from a highly unbalanced initial condition. The ultimate state of zero velocity and constant temperature is achieved up to machine accuracy. As demonstrated by the numerical experiments, the current scheme is very suitable for small amplitude perturbation and long time running under gravitational potential.     

Evaluation of Three Lattice Boltzmann Models for Particulate Flows

A comparative study is conducted to evaluate three types of lattice Boltzmann equation (LBE) models for fluid flows with finite-sized particles, including the lattice Bhatnagar-Gross-Krook (BGK) model, the model proposed by Ladd [Ladd AJC, J. Fluid Mech., 271, 285-310 (1994); Ladd AJC, J. Fluid Mech., 271, 311-339 (1994)], and the multiple-relaxation-time (MRT) model. The sedimentation of a circular particle in a two-dimensional infinite channel under gravity is used as the first test problem. The numerical results of the three LBE schemes are compared with the theoretical results and existing data. It is found that all of the three LBE schemes yield reasonable results in general, although the BGK scheme and Ladd's scheme give some deviations in some cases. Our results also show that the MRT scheme can achieve a better numerical stability than the other two schemes. Regarding the computational efficiency, it is found that the BGK scheme is the most superior one, while the other two schemes are nearly identical. We also observe that the MRT scheme can unequivocally reduce the viscosity dependence of the wall correction factor in the simulations, which reveals the superior robustness of the MRT scheme. The superiority of the MRT scheme over the other two schemes is also confirmed by the simulation of the sedimentation of an elliptical particle.     

Analysis of High-Order Absorbing Boundary Conditions for the Schrödinger Equation

The paper is concerned with the numerical solution of Schrödinger equations on an unbounded spatial domain. High-order absorbing boundary conditions for one-dimensional domain are derived, and the stability of the reduced initial boundary value problem in the computational interval is proved by energy estimate. Then a second order finite difference scheme is proposed, and the convergence of the scheme is established as well. Finally, numerical examples are reported to confirm our error estimates of the numerical methods.     

Operator Splitting for Three-Phase Flow in Heterogeneous Porous Media

We describe an operator splitting technique based on physics rather than on dimension for the numerical solution of a nonlinear system of partial differential equations which models three-phase flow through heterogeneous porous media. The model for three-phase flow considered in this work takes into account capillary forces, general relations for the relative permeability functions and variable porosity and permeability fields. In our numerical procedure a high resolution, nonoscillatory, second order, conservative central difference scheme is used for the approximation of the nonlinear system of hyperbolic conservation laws modeling the convective transport of the fluid phases. This scheme is combined with locally conservative mixed finite elements for the numerical solution of the parabolic and elliptic problems associated with the diffusive transport of fluid phases and the pressure-velocity problem. This numerical procedure has been used to investigate the existence and stability of nonclassical shock waves (called transitional or undercompressive shock waves) in two-dimensional heterogeneous flows, thereby extending previous results for one-dimensional flow problems. Numerical experiments indicate that the operator splitting technique discussed here leads to computational efficiency and accurate numerical results.     

An Explicit MUSCL Scheme on Staggered Grids with Kinetic-Like Fluxes for the Barotropic and Full Euler System

We present a second order scheme for the barotropic and full Euler equations. The scheme works on staggered grids, with numerical unknowns stored at dual locations, while the numerical fluxes are derived in the spirit of kinetic schemes. We identify stability conditions ensuring the positivity of the discrete density and energy. We illustrate the ability of the scheme to capture the structure of complex flows with 1D and 2D simulations on MAC grids.     

A Space-Time Conservative Method for Hyperbolic Systems with Stiff and Non-Stiff Source Terms

In this article we propose a higher-order space-time conservative method for hyperbolic systems with stiff and non-stiff source terms as well as relaxation systems. We call the scheme a slope propagation (SP) method. It is an extension of our scheme derived for homogeneous hyperbolic systems [1]. In the present inhomogeneous systems the relaxation time may vary from order of one to a very small value. These small values make the relaxation term stronger and highly stiff. In such situations underresolved numerical schemes may produce spurious numerical results. However, our present scheme has the capability to correctly capture the behavior of the physical phenomena with high order accuracy even if the initial layer and the small relaxation time are not numerically resolved. The scheme treats the space and time in a unified manner. The flow variables and their slopes are the basic unknowns in the scheme. The source term is treated by its volumetric integration over the space-time control volume and is a direct part of the overall space-time flux balance. We use two approaches for the slope calculations of the flow variables, the first one results directly from the flux balance over the control volumes, while in the second one we use a finite difference approach. The main features of the scheme are its simplicity, its Jacobian-free and Riemann solver-free recipe, as well as its efficiency and high order accuracy. In particular we show that the scheme has a discrete analog of the continuous asymptotic limit. We have implemented our scheme for various test models available in the literature such as the Broadwell model, the extended thermodynamics equations, the shallow water equations, traffic flow and the Euler equations with heat transfer. The numerical results validate the accuracy, versatility and robustness of the present scheme.     

An All-Regime Lagrange-Projection Like Scheme for the Gas Dynamics Equations on Unstructured Meshes

We propose an all regime Lagrange-Projection like numerical scheme for the gas dynamics equations. By all regime, we mean that the numerical scheme is able to compute accurate approximate solutions with an under-resolved discretization with respect to the Mach number M, i.e. such that the ratio between the Mach number M and the mesh size or the time step is small with respect to 1. The key idea is to decouple acoustic and transport phenomenon and then alter the numerical flux in the acoustic approximation to obtain a uniform truncation error in term of M. This modified scheme is conservative and endowed with good stability properties with respect to the positivity of the density and the internal energy. A discrete entropy inequality under a condition on the modification is obtained thanks to a reinterpretation of the modified scheme in the Harten Lax and van Leer formalism. A natural extension to multi-dimensional problems discretized over unstructured mesh is proposed. Then a simple and efficient semi-implicit scheme is also proposed. The resulting scheme is stable under a CFL condition driven by the (slow) material waves and not by the (fast) acoustic waves and so verifies the all regime property. Numerical evidences are proposed and show the ability of the scheme to deal with tests where the flow regime may vary from low to high Mach values.     

Dynamic Modeling of Phase Crossings in Two-Phase Flow

Two-phase flow and heat transfer, such as boiling and condensing flows, are complicated physical phenomena that generally prohibit an exact solution and even pose severe challenges for numerical approaches. If numerical solution time is also an issue the challenge increases even further. We present here a numerical implementation and novel study of a fully distributed dynamic one-dimensional model of two-phase flow in a tube, including pressure drop, heat transfer, and variations in tube cross-section. The model is based on a homogeneous formulation of the governing equations, discretized by a high resolution finite difference scheme due to Kurganov and Tadmore. The homogeneous formulation requires a set of thermodynamic relations to cover the entire range from liquid to gas state. This leads a number of numerical challenges since these relations introduce discontinuities in the derivative of the variables and are usually very slow to evaluate. To overcome these challenges, we use an interpolation scheme with local refinement. The simulations show that the method handles crossing of the saturation lines for both liquid to two-phase and two-phase to gas regions. Furthermore, a novel result obtained in this work, the method is stable towards dynamic transitions of the inlet/outlet boundaries across the saturation lines. Results for these cases are presented along with a numerical demonstration of conservation of mass under dynamically varying boundary conditions. Finally we present results for the stability of the code in a case of a tube with a narrow section.