Unconditional Positivity-Preserving and Energy Stable Schemes for a Reduced Poisson-Nernst-Planck System

The Poisson-Nernst-Planck (PNP) system is a widely accepted model for simulation of ionic channels. In this paper, we design, analyze, and numerically validate a second order unconditional positivity-preserving scheme for solving a reduced PNP system, which can well approximate the three dimensional ion channel problem. Positivity of numerical solutions is proven to hold true independent of the size of time steps and the choice of the Poisson solver. The scheme is easy to implement without resorting to any iteration method. Several numerical examples further confirm the positivity-preserving property, and demonstrate the accuracy, efficiency, and robustness of the proposed scheme, as well as the fast approach to steady states.     

A Conservative and Monotone Characteristic Finite Element Solver for Three-Dimensional Transport and Incompressible Navier-Stokes Equations on Unstructured Grids

We propose a mass-conservative and monotonicity-preserving characteristic finite element method for solving three-dimensional transport and incompressible Navier-Stokes equations on unstructured grids. The main idea in the proposed algorithm consists of combining a mass-conservative and monotonicity-preserving modified method of characteristics for the time integration with a mixed finite element method for the space discretization. This class of computational solvers benefits from the geometrical flexibility of the finite elements and the strong stability of the modified method of characteristics to accurately solve convection-dominated flows using time steps larger than its Eulerian counterparts. In the current study, we implement three-dimensional limiters to convert the proposed solver to a fully mass-conservative and essentially monotonicity-preserving method in addition of a low computational cost. The key idea lies on using quadratic and linear basis functions of the mesh element where the departure point is localized in the interpolation procedures. The proposed method is applied to well-established problems for transport and incompressible Navier-Stokes equations in three space dimensions. The numerical results illustrate the performance of the proposed solver and support its ability to yield accurate and efficient numerical solutions for three-dimensional convection-dominated flow problems on unstructured tetrahedral meshes.     

A Finite Element Method Solver for Time-Dependent and Stationary Schrödinger Equations with a Generic Potential

A general finite element solution of the Schrödinger equation for a one-dimensional problem is presented. The solver is applicable to both stationary and time-dependent cases with a general user-selected potential term. Furthermore, it is possible to include external magnetic or electric fields, as well as spin-orbital and spin-magnetic interactions. We use analytically soluble problems to validate the solver. The predicted numerical auto-states are compared with the analytical ones, and selected mean values are used to validate the auto-functions. In order to analyze the performance of the time-dependent Schrödinger equation, a traveling wave package benchmark was reproduced. In addition, a problem involving the scattering of a wave packet over a double potential barrier shows the performance of the solver in cases of transmission and reflection of packages. Other general problems, related to periodic potentials, are treated with the same general solver and a Lagrange multiplier method to introduce periodic boundary conditions. Some simple cases of known periodic potential solutions are reported.     

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.     

On Applicability of Poisson-Boltzmann Equation for Micro- and Nanoscale Electroosmotic Flows

The applicability of the Poisson-Boltzmann model for micro- and nanoscale electroosmotic flows is a very important theoretical and engineering problem. In this contribution we investigate this problem at two aspects: first the high ionic concentration effect on the Boltzmann distribution assumption in the diffusion layer is studied by comparisons with the molecular dynamics (MD) simulation results; then the electrical double layer (EDL) interaction effect caused by low ionic concentrations in small channels is discussed by comparing with the dynamic model described by the coupled Poisson-Nernst-Planck (PNP) and Navier-Stokes (NS) equations. The results show that the Poisson-Boltzmann (PB) model is applicable in a very wide range: (i) the PB model can still provide good predictions of the ions density profiles up to a very high ionic concentration (∼ 1 M) in the diffusion layer; (ii) the PB model predicts the net charge density accurately as long as the EDL thickness is smaller than the channel width and then overrates the net charge density profile as the EDL thickness increasing, and the predicted electric potential profile is still very accurate up to a very strong EDL interaction (λ/W ∼10).     

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.     

Towards Translational Invariance of Total Energy with Finite Element Methods for Kohn-Sham Equation

Numerical oscillation of the total energy can be observed when the Kohn-Sham equation is solved by real-space methods to simulate the translational move of an electronic system. Effectively remove or reduce the unphysical oscillation is crucial not only for the optimization of the geometry of the electronic structure, but also for the study of molecular dynamics. In this paper, we study such unphysical oscillation based on the numerical framework in [G. Bao, G. H. Hu, and D. Liu, An h-adaptive finite element solver for the calculations of the electronic structures, Journal of Computational Physics, Volume 231, Issue 14, Pages 4967-4979, 2012], and deliver some numerical methods to constrain such unphysical effect for both pseudopotential and all-electron calculations, including a stabilized cubature strategy for Hamiltonian operator, and an a posteriori error estimator of the finite element methods for Kohn-Sham equation. The numerical results demonstrate the effectiveness of our method on restraining unphysical oscillation of the total energies.     

A High Order Sharp-Interface Method with Local Time Stepping for Compressible Multiphase Flows

In this paper, a new sharp-interface approach to simulate compressible multiphase flows is proposed. The new scheme consists of a high order WENO finite volume scheme for solving the Euler equations coupled with a high order path-conservative discontinuous Galerkin finite element scheme to evolve an indicator function that tracks the material interface. At the interface our method applies ghost cells to compute the numerical flux, as the ghost fluid method. However, unlike the original ghost fluid scheme of Fedkiw et al. [15], the state of the ghost fluid is derived from an approximate-state Riemann solver, similar to the approach proposed in [25], but based on a much simpler formulation. Our formulation leads only to one single scalar nonlinear algebraic equation that has to be solved at the interface, instead of the system used in [25]. Away from the interface, we use the new general Osher-type flux recently proposed by Dumbser and Toro [13], which is a simple but complete Riemann solver, applicable to general hyperbolic conservation laws. The time integration is performed using a fully-discrete one-step scheme, based on the approaches recently proposed in [5, 7]. This allows us to evolve the system also with time-accurate local time stepping. Due to the sub-cell resolution and the subsequent more restrictive time-step constraint of the DG scheme, a local evolution for the indicator function is applied, which is matched with the finite volume scheme for the solution of the Euler equations that runs with a larger time step. The use of a locally optimal time step avoids the introduction of excessive numerical diffusion in the finite volume scheme. Two different fluids have been used, namely an ideal gas and a weakly compressible fluid modeled by the Tait equation. Several tests have been computed to assess the accuracy and the performance of the new high order scheme. A verification of our algorithm has been carefully carried out using exact solutions as well as a comparison with other numerical reference solutions. The material interface is resolved sharply and accurately without spurious oscillations in the pressure field.     

An Approximate Riemann Solver for Fluid-Solid Interaction Problems with Mie-Grüneisen Equations of State

We propose an approximate solver for compressible fluid-elastoplastic solid Riemann problems. The fluid and hydrostatic components of the solid are described by a family of general Mie-Grüneisen equations of state, and the hypo-elastoplastic constitutive law we studied includes the perfect plasticity and linearly hardened plasticity. The approximate solver provides the interface stress and normal velocity by an iterative method. The well-posedness and convergence of our solver are verified with mild assumptions on the equations of state. The proposed solver is applied in computing the numerical flux at the phase interface for our compressible multi-medium flow simulation on Eulerian girds. Several numerical examples, including Riemann problems, underground explosion and high speed impact applications, are presented to validate the approximate solver.     

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.     

Permanent Charge Effects on Ionic Flow: A Numerical Study of Flux Ratios and Their Bifurcation

Ionic flow carries electrical signals for cells to communicate with each other. The permanent charge of an ion channel is a crucial protein structure for flow properties while boundary conditions play a role of the driving force. Their effects on flow properties have been analyzed via a quasi-one-dimensional Poisson-Nernst-Planck model for small and relatively large permanent charges. The analytical studies have led to the introduction of flux ratios that reflect permanent charge effects and have a universal property. The studies also show that the flux ratios have different behaviors for small and large permanent charges. However, the existing analytical techniques can reveal neither behaviors of flux ratios nor transitions between small and large permanent charges. In this work we present a numerical investigation on flux ratios to bridge between small and large permanent charges. Numerical results verify the analytical predictions for the two extremal regions. More significantly, emergence of non-trivial behaviors is detected as the permanent charge varies from small to large. In particular, saddle-node bifurcations of flux ratios are revealed, showing rich phenomena of permanent charge effects by the power of combining analytical and numerical techniques. An adaptive moving mesh finite element method is used in the numerical studies.     

A Simple Solver for the Two-Fluid Plasma Model Based on PseudoSpectral Time-Domain Algorithm

We present a solver of 3D two-fluid plasma model for the simulation of short-pulse laser interactions with plasma. This solver resolves the equations of the two-fluid plasma model with ideal gas closure. We also include the Bhatnagar-Gross-Krook collision model. Our solver is based on PseudoSpectral Time-Domain (PSTD) method to solve Maxwell's curl equations. We use a Strang splitting to integrate Euler equations with source term: while Euler equations are solved with a composite scheme mixing Lax-Friedrichs and Lax-Wendroff schemes, the source term is integrated with a fourth-order Runge-Kutta scheme. This two-fluid plasma model solver is simple to implement because it only relies on finite difference schemes and Fast Fourier Transforms. It does not require spatially staggered grids. The solver was tested against several well-known problems of plasma physics. Numerical simulations gave results in excellent agreement with analytical solutions or with previous results from the literature.     

Mesh Sensitivity for Numerical Solutions of Phase-Field Equations Using r-Adaptive Finite Element Methods

There have been several recent papers on developing moving mesh methods for solving phase-field equations. However, it is observed that some of these moving mesh solutions are essentially different from the solutions on very fine fixed meshes. One of the purposes of this paper is to understand the reason for the differences. We carried out numerical sensitivity studies systematically in this paper and it can be concluded that for the phase-field equations, the numerical solutions are very sensitive to the starting mesh and the monitor function. As a separate issue, an efficient alternating Crank-Nicolson time discretization scheme is developed for solving the nonlinear system resulting from a finite element approximation to the phase-field equations.     

A Parallel Computational Model for Three-Dimensional,Thermo-Mechanical Stokes Flow Simulations of Glaciers and Ice Sheets

This paper focuses on the development of an efficient, three-dimensional, thermo-mechanical, nonlinear-Stokes flow computational model for ice sheet simulation. The model is based on the parallel finite element model developed in [14] which features high-order accurate finite element discretizations on variable resolution grids. Here, we add an improved iterative solution method for treating the nonlinearity of the Stokes problem, a new high-order accurate finite element solver for the temperature equation, and a new conservative finite volume solver for handling mass conservation. The result is an accurate and efficient numerical model for thermo-mechanical glacier and ice-sheet simulations. We demonstrate the improved efficiency of the Stokes solver using the ISMIP-HOM Benchmark experiments and a realistic test case for the Greenland ice-sheet. We also apply our model to the EISMINT-II benchmark experiments and demonstrate stable thermo-mechanical ice sheet evolution on both structured and unstructured meshes. Notably, we find no evidence for the "cold spoke" instabilities observed for these same experiments when using finite difference, shallow-ice approximation models on structured grids.     

Development of an Explicit Symplectic Scheme that Optimizes the Dispersion-Relation Equation of the Maxwell's Equations

In this paper an explicit finite-difference time-domain scheme for solving the Maxwell's equations in non-staggered grids is presented. The proposed scheme for solving the Faraday's and Ampère's equations in a theoretical manner is aimed to preserve discrete zero-divergence for the electric and magnetic fields. The inherent local conservation laws in Maxwell's equations are also preserved discretely all the time using the explicit second-order accurate symplectic partitioned Runge-Kutta scheme. The remaining spatial derivative terms in the semi-discretized Faraday's and Ampère's equations are then discretized to provide an accurate mathematical dispersion relation equation that governs the numerical angular frequency and the wavenumbers in two space dimensions. To achieve the goal of getting the best dispersive characteristics, we propose a fourth-order accurate space centered scheme which minimizes the difference between the exact and numerical dispersion relation equations. Through the computational exercises, the proposed dual-preserving solver is computationally demonstrated to be efficient for use to predict the long-term accurate Maxwell's solutions.     

A Stochastic Galerkin Method for Stochastic Control Problems

In an interdisciplinary field on mathematics and physics, we examine a physical problem, fluid flow in porous media, which is represented by a stochastic partial differential equation (SPDE). We first give a priori error estimates for the solutions to an optimization problem constrained by the physical model under lower regularity assumptions than the literature. We then use the concept of Galerkin finite element methods to establish a new numerical algorithm to give approximations for our stochastic optimal physical problem. Finally, we develop original computer programs based on the algorithm and use several numerical examples of various situations to see how well our solver works by comparing its outputs to the priori error estimates.     

Modeling 3D Magma Dynamics Using a Discontinuous Galerkin Method

Discontinuous Galerkin (DG) and matrix-free finite element methods with a novel projective pressure estimation are combined to enable the numerical modeling of magma dynamics in 2D and 3D using the library deal.II. The physical model is an advection-reaction type system consisting of two hyperbolic equations to evolve porosity and soluble mineral abundance at local chemical equilibrium and one elliptic equation to recover global pressure. A combination of a discontinuous Galerkin method for the advection equations and a finite element method for the elliptic equation provide a robust and efficient solution to the channel regime problems of the physical system in 3D. A projective and adaptively applied pressure estimation is employed to significantly reduce the computational wall time without impacting the overall physical reliability in the modeling of important features of melt segregation, such as melt channel bifurcation in 2D and 3D time dependent simulations.     

Numerical Simulation of Airfoil Vibrations Induced by Turbulent Flow

The subject of the paper is the numerical simulation of the interaction of two-dimensional incompressible viscous flow and a vibrating airfoil with large amplitudes. The airfoil with three degrees of freedom performs rotation around an elastic axis, oscillations in the vertical direction and rotation of a flap. The numerical simulation consists of the finite element solution of the Reynolds averaged Navier-Stokes equations combined with Spalart-Allmaras or k−ω turbulence models, coupled with a system of nonlinear ordinary differential equations describing the airfoil motion with consideration of large amplitudes. The time-dependent computational domain and approximation on a moving grid are treated by the Arbitrary Lagrangian-Eulerian formulation of the flow equations. Due to large values of the involved Reynolds numbers an application of a suitable stabilization of the finite element discretization is employed. The developed method is used for the computation of flow-induced oscillations of the airfoil near the flutter instability, when the displacements of the airfoil are large, up to ±40 degrees in rotation. The paper contains the comparison of the numerical results obtained by both turbulence models.     

Immersed Finite Element Method for Interface Problems with Algebraic Multigrid Solver

This article is to discuss the bilinear and linear immersed finite element (IFE) solutions generated from the algebraic multigrid solver for both stationary and moving interface problems. For the numerical methods based on finite difference formulation and a structured mesh independent of the interface, the stiffness matrix of the linear system is usually not symmetric positive-definite, which demands extra efforts to design efficient multigrid methods. On the other hand, the stiffness matrix arising from the IFE methods is naturally symmetric positive-definite. Hence the IFE-AMG algorithm is proposed to solve the linear systems of the bilinear and linear IFE methods for both stationary and moving interface problems. The numerical examples demonstrate the features of the proposed algorithms, including the optimal convergence in both $L^2$ and semi-$H^1$ norms of the IFE-AMG solutions, the high efficiency with proper choice of the components and parameters of AMG, the influence of the tolerance and the smoother type of AMG on the convergence of the IFE solutions for the interface problems, and the relationship between the cost and the moving interface location.