High-Order Positivity-Preserving Well-Balanced Discontinuous Galerkin Methods for Euler Equations with Gravitation on Unstructured Meshes

In this paper, we propose a high-order accurate discontinuous Galerkin (DG) method for the compressible Euler equations under gravitational fields on unstructured meshes. The scheme preserves a general hydrostatic equilibrium state and provably guarantees the positivity of density and pressure at the same time. Comparing with the work on the well-balanced scheme for Euler equations with gravitation on rectangular meshes, the extension to triangular meshes is conceptually plausible but highly nontrivial. We first introduce a special way to recover the equilibrium state and then design a group of novel variables at the interface of two adjacent cells, which plays an important role in the well-balanced and positivity-preserving properties. One main challenge is that the well-balanced schemes may not have the weak positivity property. In order to achieve the well-balanced and positivity-preserving properties simultaneously while maintaining high-order accuracy, we carefully design DG spatial discretization with well-balanced numerical fluxes and suitable source term approximation. For the ideal gas, we prove that the resulting well-balanced scheme, coupled with strong stability preserving time discretizations, satisfies a weak positivity property. A simple existing limiter can be applied to enforce the positivity-preserving property, without losing high-order accuracy and conservation. Extensive one- and two-dimensional numerical examples demonstrate the desired properties of the proposed scheme, as well as its high resolution and robustness.     

A Colocalized Scheme for Three-Temperature Grey Diffusion Radiation Hydrodynamics

A positivity-preserving, conservative and entropic numerical scheme is presented for the three-temperature grey diffusion radiation hydrodynamics model. More precisely, the dissipation matrices of the colocalized semi-Lagrangian scheme are defined in order to enforce the entropy production on each species (electron or ion) proportionally to its mass as prescribed in [34]. A reformulation of the model is then considered to enable the derivation of a robust convex combination based scheme. This yields the positivity-preserving property at each sub-iteration of the algorithm while the total energy conservation is reached at convergence. Numerous pure hydrodynamics and radiation hydrodynamics test cases are carried out to assess the accuracy of the method. The question of the stability of the scheme is also addressed. It is observed that the present numerical method is particularly robust.     

A Positivity-Preserving Scheme for the Simulation of Streamer Discharges in Non-Attaching and Attaching Gases

Assumed having axial symmetry, the streamer discharge is often described by a fluid model in cylindrical coordinate system, which consists of convection dominated (diffusion) equations with source terms, coupled with a Poisson's equation. Without additional care for a stricter CFL condition or special treatment to the negative source term, popular methods used in streamer discharge simulations, e.g., FEM-FCT, FVM, cannot ensure the positivity of the particle densities for the cases in attaching gases. By introducing the positivity-preserving limiter proposed by Zhang and Shu [15] and Strang operator splitting, this paper proposes a finite difference scheme with a provable positivity-preserving property in cylindrical coordinate system, for the numerical simulation of streamer discharges in non-attaching and attaching gases. Numerical examples in non-attaching gas (N₂) and attaching gas (SF₆) are given to illustrate the effectiveness of the scheme.     

Adaptive Order WENO Reconstructions for the Semi-Lagrangian Finite Difference Scheme for Advection Problem

We present a new conservative semi-Lagrangian finite difference weighted essentially non-oscillatory scheme with adaptive order. This is an extension of the conservative semi-Lagrangian (SL) finite difference WENO scheme in [Qiu and Shu, JCP, 230 (4) (2011), pp. 863-889], in which linear weights in SL WENO framework were shown not to exist for variable coefficient problems. Hence, the order of accuracy is not optimal from reconstruction stencils. In this paper, we incorporate a recent WENO adaptive order (AO) technique [Balsara et al., JCP, 326 (2016), pp. 780-804] to the SL WENO framework. The new scheme can achieve an optimal high order of accuracy, while maintaining the properties of mass conservation and non-oscillatory capture of solutions from the original SL WENO. The positivity-preserving limiter is further applied to ensure the positivity of solutions. Finally, the scheme is applied to high dimensional problems by a fourth-order dimensional splitting. We demonstrate the effectiveness of the new scheme by extensive numerical tests on linear advection equations, the Vlasov-Poisson system, the guiding center Vlasov model as well as the incompressible Euler equations.     

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.     

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.     

A High Order Central DG Method of the Two-Layer Shallow Water Equations

In this paper, we focus on the numerical simulation of the two-layer shallow water equations over variable bottom topography. Although the existing numerical schemes for the single-layer shallow water equations can be extended to two-layer shallow water equations, it is not a trivial work due to the complexity of the equations. To achieve the well-balanced property of the numerical scheme easily, the two-layer shallow water equations are reformulated into a new form by introducing two auxiliary variables. Since the new equations are only conditionally hyperbolic and their eigenstructure cannot be easily obtained, we consider the utilization of the central discontinuous Galerkin method which is free of Riemann solvers. By choosing the values of the auxiliary variables suitably, we can prove that the scheme can exactly preserve the still-water solution, and thus it is a truly well-balanced scheme. To ensure the non-negativity of the water depth, a positivity-preserving limiter and a special approximation to the bottom topography are employed. The accuracy and validity of the numerical method will be illustrated through some numerical tests.     

Anomalous Sorption Kinetics of Self-Interacting Particles by a Spherical Trap

In this paper we propose a computational framework for the investigation of the correlated motion between positive and negative ions exposed to the attraction of a bubble surface that mimics the (oscillating) cell membrane. Specifically we aim to investigate the role of surface traps with substances freely diffusing around the cell. The physical system we want to model is an anchored gas drop submitted to a diffusive flow of charged surfactants (ions). When the diffusing surfactants meet the surface of the bubble, they are reversibly adsorbed and their local concentration is accurately measured. The correlated diffusion of surfactants is described by a Poisson-Nernst-Planck (PNP) system, in which the drift term is given by the gradient of a potential which includes both the effect of the bubble and the Coulomb interaction between the carriers. The latter term is obtained from the solution of a self-consistent Poisson equation. For very short Debye lengths one can adopt the so called Quasi-Neutral limit which drastically simplifies the system, thus allowing for much faster numerical simulations. The paper has four main objectives. The first one is to present a PNP model that describes ion charges in presence of a trap. The second one is to provide benchmark tests for the validation of simplified multiscale models under current development [1]. The third one is to explore the relevance of the term describing the interaction among the apolar tails of the anions. The last one is to quantitatively explore the validity of the Quasi-Neutral limit by comparison with detailed numerical simulation for smaller and smaller Debye lengths. In order to reach these goals, we propose a simple and efficient Alternate Direction Implicit method for the numerical solution of the non-linear PNP system, which guarantees second order accuracy bothin space and time, without requiring solution of nonlinear equation at each time step. New semi-implicit scheme for a simplified PNP system near quasi neutrality is also proposed.     

Two Nonlinear Positivity-Preserving Finite Volume Schemes for Three-Dimensional Heat Conduction Equations on General Polyhedral Meshes

In this article we present two types of nonlinear positivity-preserving finite volume (PPFV) schemes for a class of three-dimensional heat conduction equations on general polyhedral meshes. First, we present a new parameter selection strategy on the one-sided flux and establish a nonlinear PPFV scheme based on a two-point flux with higher efficiency. By comparing with the scheme proposed in [H. Xie, X. Xu, C. Zhai, H. Yong, Commun. Comput. Phys. 24 (2018) 1375–1408], our scheme avoids the assumption that the values of auxiliary unknowns are nonnegative, which makes our interpolation formulae suitable to be constructed by existing approaches with high accuracy and well robustness (e.g., the finite element method), thus enhancing the adaptability to distorted meshes with large deformations. Then we derive a linear multi-point flux involving combination coefficients and, via the Patankar trick, obtain another nonlinear PPFV scheme that is concise and easy to implement. The selection strategy of combination coefficients is also provided to improve the convergence behavior of the Picard procedure. Furthermore, the existence and positivity-preserving properties of these two nonlinear PPFV solutions are proved. Numerical experiments with the discontinuous diffusion scalar as well as discontinuous and anisotropic diffusion tensors are given to confirm our theoretical findings and demonstrate that our schemes both can achieve ideal-order accuracy even on severely distorted meshes.     

A Decoupled and Positivity-Preserving DDFVS Scheme for Diffusion Problems on Polyhedral Meshes

We propose a decoupled and positivity-preserving discrete duality finite volume (DDFV) scheme for anisotropic diffusion problems on polyhedral meshes with star-shaped cells and planar faces. Under the generalized DDFV framework, two sets of finite volume (FV) equations are respectively constructed on the dual and primary meshes, where the ones on the dual mesh are derived from the ingenious combination of a geometric relationship with the construction of the cell matrix. The resulting system on the dual mesh is symmetric and positive definite, while the one on the primary mesh possesses an M-matrix structure. To guarantee the positivity of the two categories of unknowns, a cutoff technique is introduced. As for the local conservation, it is conditionally maintained on the dual mesh while strictly preserved on the primary mesh. More interesting is that the FV equations on the dual mesh can be solved independently, so that the two sets of FV equations are decoupled. As a result, no nonlinear iteration is required for linear problems and a general nonlinear solver could be used for nonlinear problems. In addition, we analyze the well-posedness of numerical solutions for linear problems. The properties of the presented scheme are examined by numerical experiments. The efficiency of the Newton method is also demonstrated by comparison with those of the fixed-point iteration method and its Anderson acceleration.     

Field Model for Complex Ionic Fluids: Analytical Properties and Numerical Investigation

In this paper, we consider the field model for complex ionic fluids with an energy variational structure, and analyze the well-posedness to this model with regularized kernels. Furthermore, we deduce the estimate of the maximal density function to quantify the finite size effect. On the numerical side, we adopt a finite volume scheme to the field model, which satisfies the following properties: positivity-preserving, mass conservation and energy dissipation. Besides, series of numerical experiments are provided to demonstrate the properties of the steady state and the finite size effect by showing the equilibrium profiles with different values of the parameter in the kernel.     

A Well-Balanced Positivity-Preserving Quasi-Lagrange Moving Mesh DG Method for the Shallow Water Equations

A high-order, well-balanced, positivity-preserving quasi-Lagrange moving mesh DG method is presented for the shallow water equations with non-flat bottom topography. The well-balance property is crucial to the ability of a scheme to simulate perturbation waves over the lake-at-rest steady state such as waves on a lake or tsunami waves in the deep ocean. The method combines a quasi-Lagrange moving mesh DG method, a hydrostatic reconstruction technique, and a change of unknown variables. The strategies in the use of slope limiting, positivity-preservation limiting, and change of variables to ensure the well-balance and positivity-preserving properties are discussed. Compared to rezoning-type methods, the current method treats mesh movement continuously in time and has the advantages that it does not need to interpolate flow variables from the old mesh to the new one and places no constraint for the choice of a update scheme for the bottom topography on the new mesh. A selection of one- and two-dimensional examples are presented to demonstrate the well-balance property, positivity preservation, and high-order accuracy of the method and its ability to adapt the mesh according to features in the flow and bottom topography.     

Numerical Solution of 3D Poisson-Nernst-Planck Equations Coupled with Classical Density Functional Theory for Modeling Ion and Electron Transport in a Confined Environment

We have developed efficient numerical algorithms for solving 3D steady-state Poisson-Nernst-Planck (PNP) equations with excess chemical potentials described by the classical density functional theory (cDFT). The coupled PNP equations are discretized by a finite difference scheme and solved iteratively using the Gummel method with relaxation. The Nernst-Planck equations are transformed into Laplace equations through the Slotboom transformation. Then, the algebraic multigrid method is applied to efficiently solve the Poisson equation and the transformed Nernst-Planck equations. A novel strategy for calculating excess chemical potentials through fast Fourier transforms is proposed, which reduces computational complexity from $\mathcal{O}$($N^2$) to $\mathcal{O}$($NlogN$), where $N$ is the number of grid points. Integrals involving the Dirac delta function are evaluated directly by coordinate transformation, which yields more accurate results compared to applying numerical quadrature to an approximated delta function. Numerical results for ion and electron transport in solid electrolyte for lithium-ion (Li-ion) batteries are shown to be in good agreement with the experimental data and the results from previous studies.     

An Efficient Positivity-Preserving Finite Volume Scheme for the Nonequilibrium Three-Temperature Radiation Diffusion Equations on Polygonal Meshes

This paper develops an efficient positivity-preserving finite volume scheme for the two-dimensional nonequilibrium three-temperature radiation diffusion equations on general polygonal meshes. The scheme is formed as a predictor-corrector algorithm. The corrector phase obtains the cell-centered solutions on the primary mesh, while the predictor phase determines the cell-vertex solutions on the dual mesh independently. Moreover, the flux on the primary edge is approximated with a fixed stencil and the nonnegative cell-vertex solutions are not reconstructed. Theoretically, our scheme does not require any nonlinear iteration for the linear problems, and can call the fast nonlinear solver (e.g. Newton method) for the nonlinear problems. The positivity, existence and uniqueness of the cell-centered solutions obtained on the corrector phase are analyzed, and the scheme on quasi-uniform meshes is proved to be $L^2$- and $H^1$-stable under some assumptions. Numerical experiments demonstrate the accuracy, efficiency and positivity of the scheme on various distorted meshes.     

Construction,Analysis and Application of Coupled Compact Difference Scheme in Computational Acoustics and Fluid Flow Problems

In the present work, a new type of coupled compact difference scheme has been proposed for the solution of computational acoustics and flow problems. The proposed scheme evaluates the first, the second and the fourth derivative terms simultaneously. Derived compact difference scheme has a significant spectral resolution and a physical dispersion relation preserving (DRP) ability over a considerable wavenumber range when a fourth order four stage Runge-Kutta scheme is used for the time integration. Central stencil has been used for the present numerical scheme to evaluate spatial derivative terms. Derived scheme has the capability of adding numerical diffusion adaptively to attenuate spurious high wavenumber oscillations responsible for numerical instabilities. The DRP nature of the proposed scheme across a wider wavenumber range provides accurate results for the model wave equations as well as computational acoustic problems. In addition to the attractive feature of adaptive diffusion, present scheme also helps to control spurious reflections from the domain boundaries and is projected as an alternative to the perfectly matched layer (PML) technique.     

An Implicit Scheme for Moving Walls and Multi-Material Interfaces in Weakly Compressible Materials

We propose a numerical method for the simulation of flows from weakly compressible to low Mach regimes in domains with moving boundaries. Non-miscible weakly compressible materials separated by an interface are included as well. The scheme is fully implicit and it exploits the relaxation all-speed scheme introduced in [1]. We consider media with significantly different physical properties and constitutive laws, as fluids and hyperelastic solids. The proposed numerical scheme is fully Eulerian and it is the same for all materials. We present numerical validations by simulating weakly compressible fluid/fluid, solid/solid and solid/fluid interactions.     

Some Random Batch Particle Methods for the Poisson-Nernst-Planck and Poisson-Boltzmann Equations

We consider in this paper random batch interacting particle methods for solving the Poisson-Nernst-Planck (PNP) equations, and thus the Poisson-Boltzmann (PB) equation as the equilibrium, in the external unbounded domain. To justify the simulation in a truncated domain, an error estimate of the truncation is proved in the symmetric cases for the PB equation. Then, the random batch interacting particle methods are introduced which are $\mathcal{O}(N)$ per time step. The particle methods can not only be considered as a numerical method for solving the PNP and PB equations, but also can be used as a direct simulation approach for the dynamics of the charged particles in solution. The particle methods are preferable due to their simplicity and adaptivity to complicated geometry, and may be interesting in describing the dynamics of the physical process. Moreover, it is feasible to incorporate more physical effects and interactions in the particle methods and to describe phenomena beyond the scope of the mean-field equations.     

Development of a High-Resolution Scheme for Solving the PNP-NS Equations in Curved Channels

A high-order finite difference scheme has been developed to approximate the spatial derivative terms present in the unsteady Poisson-Nernst-Planck (PNP) equations and incompressible Navier-Stokes (NS) equations. Near the wall the sharp solution profiles are resolved by using the combined compact difference (CCD) scheme developed in five-point stencil. This CCD scheme has a sixth-order accuracy for the second-order derivative terms while a seventh-order accuracy for the first-order derivative terms. PNP-NS equations have been also transformed to the curvilinear coordinate system to study the effects of channel shapes on the development of electroosmotic flow. In this study, the developed scheme has been analyzed rigorously through the modified equation analysis. In addition, the developed method has been computationally verified through four problems which are amenable to their own exact solutions. The electroosmotic flow details in planar and wavy channels have been explored with the emphasis on the formation of Coulomb force. Significance of different forces resulting from the pressure gradient, diffusion and Coulomb origins on the convective electroosmotic flow motion is also investigated in detail.     

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.