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 O(N2) to 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 lithiumion (Li-ion) batteries are shown to be in good agreement with the experimental data and the results from previous studies.

We study the settling of solid particles in a viscous incompressible fluid contained within a two-dimensional channel, where the mass density of the particles is greater than that of the fluid. The fluid-structure interaction problem is simulated numerically using the immersed boundary method, where the added mass is incorporated using a Boussinesq approximation. Simulations are performed with a single circular particle, and also with two particles in various initial configurations. The terminal particle settling velocity and drag coefficient correspond closely with other theoretical, experimental and numerical results, and the particle trajectories reproduce the expected behavior qualitatively. In particular, simulations of a pair of interacting particles similar drafting-kissing-tumbling dynamics to that observed in other experimental and numerical studies.

In this paper we present a generalized lattice Boltzmann based approach for sliding-mesh local reference frame. This scheme exactly conserves hydrodynamic fluxes across local reference frame interface. The accuracy and robustness of our scheme are demonstrated by benchmark validations.

A discontinuous Galerkin method for the ideal 5 moment two-fluid plasma system is presented. The method uses a second or third order discontinuous Galerkin spatial discretization and a third order TVD Runge-Kutta time stepping scheme. The method is benchmarked against an analytic solution of a dispersive electron acoustic square pulse as well as the two-fluid electromagnetic shock [1] and existing numerical solutions to the GEM challenge magnetic reconnection problem [2]. The algorithm can be generalized to arbitrary geometries and three dimensions. An approach to maintaining small gauge errors based on error propagation is suggested.

In this paper, we derive a multi-symplectic Fourier pseudospectral scheme for the Kawahara equation with special attention to the relationship between the spectral differentiation matrix and discrete Fourier transform. The relationship is crucial for implementing the scheme efficiently. By using the relationship, we can apply the Fast Fourier transform to solve the Kawahara equation. The effectiveness of the proposed methods will be demonstrated by a number of numerical examples. The numerical results also confirm that the global energy and momentum are well preserved.

This paper presents limits for stability of projection type schemes when using high order pressure-velocity pairs of same degree. Two high order h/p variational methods encompassing continuous and discontinuous Galerkin formulations are used to explain previously observed lower limits on the time step for projection type schemes to be stable [18], when h- or p-refinement strategies are considered. In addition, the analysis included in this work shows that these stability limits do not depend only on the time step but on the product of the latter and the kinematic viscosity, which is of particular importance in the study of high Reynolds number flows. We show that high order methods prove advantageous in stabilising the simulations when small time steps and low kinematic viscosities are used.

Drawing upon this analysis, we demonstrate how the effects of this instability can be reduced in the discontinuous scheme by introducing a stabilisation term into the global system. Finally, we show that these lower limits are compatible with Courant-Friedrichs-Lewy (CFL) type restrictions, given that a sufficiently high polynomial order or a mall enough mesh spacing is selected.

The Brinkman model describes flow of fluid in complex porous media with a high-contrast permeability coefficient such that the flow is dominated by Darcy in some regions and by Stokes in others. A weak Galerkin (WG) finite element method for solving the Brinkman equations in two or three dimensional spaces by using polynomials is developed and analyzed. The WG method is designed by using the generalized functions and their weak derivatives which are defined as generalized distributions. The variational form we considered in this paper is based on two gradient operators which is different from the usual gradient-divergence operators for Brinkman equations. The WG method is highly flexible by allowing the use of discontinuous functions on arbitrary polygons or polyhedra with certain shape regularity. Optimal-order error estimates are established for the corresponding WG finite element solutions in various norms. Some computational results are presented to demonstrate the robustness, reliability, accuracy, and flexibility of the WG method for the Brinkman equations.

We investigate several robust preconditioners for solving the saddle-point linear systems that arise from spatial discretization of unsteady and steady variable-coefficient Stokes equations on a uniform staggered grid. Building on the success of using the classical projection method as a preconditioner for the coupled velocity pressure system [B. E. Griffith, J. Comp. Phys., 228 (2009), pp. 7565-7595], as well; established techniques for steady and unsteady Stokes flow in the finite-element literature, we construct preconditioners that employ independent generalized Helmholtz and Poisson solvers for the velocity and pressure subproblems. We demonstrate that only a single cycle of a standard geometric multigrid algorithm serves as an effective inexact solver for each of these subproblems. Contrary to traditional wisdom, we find that the Stokes problem can be solved nearly as efficiently as the independent pressure and velocity subproblems, making the overall cost of solving the Stokes system comparable to the cost of classical projection or fractional step methods for incompressible flow, even for steady flow and in the presence of large density and viscosity contrasts. Two of the five preconditioners considered here are found to be robust to GMRES restarts and to increasing problem size, making them suitable for large-scale problems. Our work opens many possibilities for constructing novel unsplit temporal integrators for finite-volume spatial discretizations of the equations of low Mach and incompressible flow dynamics.

The stochastic collocation method using sparse grids has become a popular choice for performing stochastic computations in high dimensional (random) parameter space. In addition to providing highly accurate stochastic solutions, the sparse grid collocation results naturally contain sensitivity information with respect to the input random parameters. In this paper, we use the sparse grid interpolation and cubature methods of Smolyak together with combinatorial analysis to give a computationally efficient method for computing the global sensitivity values of Sobol’. This method allows for approximation of all main effect and total effect values from evaluation of f on a single set of sparse grids. We discuss convergence of this method, apply it to several test cases and compare to existing methods. As a result which may be of independent interest, we recover an explicit formula for evaluating a Lagrange basis interpolating polynomial associated with the Chebyshev extrema. This allows one to manipulate the sparse grid collocation results in a highly efficient manner.

A weak Galerkin (WG) method is introduced and numerically tested for the Helmholtz equation. This method is flexible by using discontinuous piecewise polynomials and retains the mass conservation property. At the same time, the WG finite element formulation is symmetric and parameter free. Several test scenarios are designed for a numerical investigation on the accuracy, convergence, and robustness of the WG method in both inhomogeneous and homogeneous media over convex and non-convex domains. Challenging problems with high wave numbers are also examined. Our numerical experiments indicate that the weak Galerkin is a finite element technique that is easy to implement, and provides very accurate and robust numerical solutions for the Helmholtz problem with high wave numbers.

Numerical simulation of two-phase flow in fractured karst reservoirs is still a challenging issue. The triple-porosity model is the major approach up to now. However, the triple-continuum assumption in this model is unacceptable for many cases. In the present work, an efficient numerical model has been developed for immiscible two-phase flow in fractured karst reservoirs based on the idea of equivalent continuum representation. First, based on the discrete fracture-vug model and homogenization theory, the effective absolute permeability tensors for each grid blocks are calculated. And then an analytical procedure to obtain a pseudo relative permeability curves for a grid block containing fractures and cavities has been successfully implemented. Next, a full-tensor simulator has been designed based on a hybrid numerical method (combining mixed finite element method and finite volume method). A simple fracture system has been used to demonstrate the validity of our method. At last, we have used the fracture and cavity statistics data from TAHE outcrops in west China, effective permeability values and other parameters from our code, and an equivalent continuum simulator to calculate the water flooding profiles for more realistic systems.

In this paper, the power-law fluid flows in a two-dimensional square cavity are investigated in detail with multi-relaxation-time lattice Boltzmann method (MRT-LBM). The influence of the Reynolds number (Re) and the power-law index (n) on the vortex strength, vortex position and velocity distribution are extensively studied. In our numerical simulations, Re is varied from 100 to 10000, and n is ranged from 0.25 to 1.75, covering both cases of shear-thinning and shear-thickening. Compared with the Newtonian fluid, numerical results show that the flow structure and number of vortex of power-law fluid are not only dependent on the Reynolds number, but also related to power-law index.

This work is concerned with the modelling of the interaction of fluid flow with flexibly supported rigid bodies. The fluid flow is modelled by Lattice-Boltzmann Method, coupled to a set of ordinary differential equations describing the dynamics of the solid body in terms its elastic and damping properties. The time discretization of the body dynamics is performed via the Time Discontinuous Galerkin Method. Several numerical examples are presented and highlight the robustness and efficiency of the proposed methodology, by means of comparisons with previously published results. The examples show that the present fluid-structure method is able to capture vortex- induced oscillations of flexibly-supported rigid body.

We employ the lattice Boltzmann method and random walk particle tracking to simulate the time evolution of hydrodynamic dispersion in bulk, random, monodisperse, hard-sphere packings with bed porosities (interparticle void volume fractions) between the random-close and the random-loose packing limit. Using Jodrey-Tory and Monte Carlo-based algorithms and a systematic variation of the packing protocols we generate a portfolio of packings, whose microstructures differ in their degree of heterogeneity (DoH). Because the DoH quantifies the heterogeneity of the void space distribution in a packing, the asymptotic longitudinal dispersion coefficient calculated for the packings increases with the packings’ DoH. We investigate the influence of packing length (up to 150 dp, where dp is the sphere diameter) and grid resolution (up to 90 nodes per dp) on the simulated hydrodynamic dispersion coefficient, and demonstrate that the chosen packing dimensions of 10 dpx 10 dpx 70 dp and the employed grid resolution of 60 nodes per dp are sufficient to observe asymptotic behavior of the dispersion coefficient and to minimize finite size effects. Asymptotic values of the dispersion coefficients calculated for the generated packings are compared with simulated as well as experimental data from the literature and yield good to excellent agreement.

In many problems, one wishes to solve the Helmholtz equation with variable coefficients within the Laplacian-like term and use a high order accurate method (e.g., fourth order accurate) to alleviate the points-per-wavelength constraint by reducing the dispersion errors. The variation of coefficients in the equation may be due to an inhomogeneous medium and/or non-Cartesian coordinates. This renders existing fourth order finite difference methods inapplicable. We develop a new compact scheme that is provably fourth order accurate even for these problems. We present numerical results that corroborate the fourth order convergence rate for several model problems.

This article describes a number of velocity-based moving mesh numerical methods for multidimensional nonlinear time-dependent partial differential equations (PDEs). It consists of a short historical review followed by a detailed description of a recently developed multidimensional moving mesh finite element method based on conservation. Finite element algorithms are derived for both mass-conserving and non mass-conserving problems, and results shown for a number of multidimensional nonlinear test problems, including the second order porous medium equation and the fourth order thin film equation as well as a two-phase problem. Further applications and extensions are referenced.

In this paper we present a selective segmentation model using a dual level set variational formulation. Our variational model aims to segment all objects with one level set function (global) and the selected object, which is the closest to the geometric constraints (markers), with another level set (local). It is a combination of edge detection, markers distance function and active contour without edges. Experimental results show that our model is more robust than previous work.

A reconstruction-based discontinuous Galerkin (RDG(P1P2)) method, a variant of P1P2 method, is presented for the solution of the compressible Euler equations on arbitrary grids. In this method, an in-cell reconstruction, designed to enhance the accuracy of the discontinuous Galerkin method, is used to obtain a quadratic polynomial solution (P2) from the underlying linear polynomial (P1) discontinuous Galerkin solution using a least-squares method. The stencils used in the reconstruction involve only the von Neumann neighborhood (face-neighboring cells) and are compact and consistent with the underlying DG method. The developed RDG method is used to compute a variety of flow problems on arbitrary meshes to demonstrate its accuracy, efficiency, robustness, and versatility. The numerical results indicate that this RDG(P1P2) method is third-order accurate, and outperforms the third-order DG method (DG(P2)) in terms of both computing costs and storage requirements.

This paper extends the adaptive moving mesh method developed by Tang and Tang [36] to two-dimensional (2D) relativistic hydrodynamic (RHD) equations. The algorithm consists of two “independent” parts: the time evolution of the RHD equations and the (static) mesh iteration redistribution. In the first part, the RHD equations are discretized by using a high resolution finite volume scheme on the fixed but nonuniform meshes without the full characteristic decomposition of the governing equations. The second part is an iterative procedure. In each iteration, the mesh points are first redistributed, and then the cell averages of the conservative variables are remapped onto the new mesh in a conservative way. Several numerical examples are given to demonstrate the accuracy and effectiveness of the proposed method.

We present new large time step methods for the shallow water flows in the low Froude number limit. In order to take into account multiscale phenomena that typically appear in geophysical flows nonlinear fluxes are split into a linear part governing the gravitational waves and the nonlinear advection. We propose to approximate fast linear waves implicitly in time and in space by means of a genuinely multidimensional evolution operator. On the other hand, we approximate nonlinear advection part explicitly in time and in space by means of the method of characteristics or some standard numerical flux function. Time integration is realized by the implicit-explicit (IMEX) method. We apply the IMEX Euler scheme, two step Runge Kutta Cranck Nicolson scheme, as well as the semi-implicit BDF scheme and prove their asymptotic preserving property in the low Froude number limit. Numerical experiments demonstrate stability, accuracy and robustness of these new large time step finite volume schemes with respect to small Froude number.