Kinetic Energy Preserving and Entropy Stable Finite Volume Schemes for Compressible Euler and Navier-Stokes Equations

Centered numerical fluxes can be constructed for compressible Euler equations which preserve kinetic energy in the semi-discrete finite volume scheme. The essential feature is that the momentum flux should be of the form where are any consistent approximations to the pressure and the mass flux. This scheme thus leaves most terms in the numerical flux unspecified and various authors have used simple averaging. Here we enforce approximate or exact entropy consistency which leads to a unique choice of all the terms in the numerical fluxes. As a consequence, a novel entropy conservative flux that also preserves kinetic energy for the semi-discrete finite volume scheme has been proposed. These fluxes are centered and some dissipation has to be added if shocks are present or if the mesh is coarse. We construct scalar artificial dissipation terms which are kinetic energy stable and satisfy approximate/exact entropy condition. Secondly, we use entropy-variable based matrix dissipation flux which leads to kinetic energy and entropy stable schemes. These schemes are shown to be free of entropy violating solutions unlike the original Roe scheme. For hypersonic flows a blended scheme is proposed which gives carbuncle free solutions for blunt body flows. Numerical results for Euler and Navier-Stokes equations are presented to demonstrate the performance of the different schemes.     

A Fourth Order Numerical Method for the Primitive Equations Formulated in Mean Vorticity

A fourth-order finite difference method is proposed and studied for the primitive equations (PEs) of large-scale atmospheric and oceanic flow based on mean vorticity formulation. Since the vertical average of the horizontal velocity field is divergence-free, we can introduce mean vorticity and mean stream function which are connected by a 2-D Poisson equation. As a result, the PEs can be reformulated such that the prognostic equation for the horizontal velocity is replaced by evolutionary equations for the mean vorticity field and the vertical derivative of the horizontal velocity. The mean vorticity equation is approximated by a compact difference scheme due to the difficulty of the mean vorticity boundary condition, while fourth-order long-stencil approximations are utilized to deal with transport type equations for computational convenience. The numerical values for the total velocity field (both horizontal and vertical) are statically determined by a discrete realization of a differential equation at each fixed horizontal point. The method is highly efficient and is capable of producing highly resolved solutions at a reasonable computational cost. The full fourth-order accuracy is checked by an example of the reformulated PEs with force terms. Additionally, numerical results of a large-scale oceanic circulation are presented.     

Enforcing the Discrete Maximum Principle for Linear Finite Element Solutions of Second-Order Elliptic Problems

The maximum principle is a basic qualitative property of the solution of second-order elliptic boundary value problems. The preservation of the qualitative characteristics, such as the maximum principle, in discrete model is one of the key requirements. It is well known that standard linear finite element solution does not satisfy maximum principle on general triangular meshes in 2D. In this paper we consider how to enforce discrete maximum principle for linear finite element solutions for the linear second-order self-adjoint elliptic equation. First approach is based on repair technique, which is a posteriori correction of the discrete solution. Second method is based on constrained optimization. Numerical tests that include anisotropic cases demonstrate how our method works for problems for which the standard finite element methods produce numerical solutions that violate the discrete maximum principle.     

Convergent Adaptive Finite Element Method Based on Centroidal Voronoi Tessellations and Superconvergence

We present a novel adaptive finite element method (AFEM) for elliptic equations which is based upon the Centroidal Voronoi Tessellation (CVT) and superconvergent gradient recovery. The constructions of CVT and its dual Centroidal Voronoi Delaunay Triangulation (CVDT) are facilitated by a localized Lloyd iteration to produce almost equilateral two dimensional meshes. Working with finite element solutions on such high quality triangulations, superconvergent recovery methods become particularly effective so that asymptotically exact a posteriori error estimations can be obtained. Through a seamless integration of these techniques, a convergent adaptive procedure is developed. As demonstrated by the numerical examples, the new AFEM is capable of solving a variety of model problems and has great potential in practical applications.     

Adaptive hp-Finite Element Computations for Time-Harmonic Maxwell's Equations

In this paper, hp-adaptive finite element methods are studied for time-harmonic Maxwell's equations. We propose the parallel hp-adaptive algorithms on conforming unstructured tetrahedral meshes based on residual-based a posteriori error estimates. Extensive numerical experiments are reported to investigate the efficiency of the hp-adaptive methods for point singularities, edge singularities, and an engineering benchmark problem of Maxwell's equations. The hp-adaptive methods show much better performance than the h-adaptive method.     

One-Dimensional Blood Flow with Discontinuous Properties and Transport: Mathematical Analysis and Numerical Schemes

In this paper we consider the one-dimensional blood flow model with discontinuous mechanical and geometrical properties, as well as passive scalar transport, proposed in [E.F. Toro and A. Siviglia. Flow in collapsible tubes with discontinuous mechanical properties: mathematical model and exact solutions. Communications in Computational Physics. 13(2), 361-385, 2013], completing the mathematical analysis by providing new propositions and new proofs of relations valid across different waves. Next we consider a first order DOT Riemann solver, proposing an integration path that incorporates the passive scalar and proving the well-balanced properties of the resulting numerical scheme for stationary solutions. Finally we describe a novel and simple well-balanced, second order, non-linear numerical scheme to solve the equations under study; by using suitable test problems for which exact solutions are available, we assess the well-balanced properties of the scheme, its capacity to provide accurate solutions in challenging flow conditions and its accuracy.     

A Numerical Methodology for Enforcing Maximum Principles and the Non-Negative Constraint for Transient Diffusion Equations

Transient diffusion equations arise in many branches of engineering and applied sciences (e.g., heat transfer and mass transfer), and are parabolic partial differential equations. It is well-known that these equations satisfy important mathematical properties like maximum principles and the non-negative constraint, which have implications in mathematical modeling. However, existing numerical formulations for these types of equations do not, in general, satisfy maximum principles and the non-negative constraint. In this paper, we present a methodology for enforcing maximum principles and the non-negative constraint for transient anisotropic diffusion equation. The proposed methodology is based on the method of horizontal lines in which the time is discretized first. This results in solving steady anisotropic diffusion equation with decay equation at every discrete time-level. We also present other plausible temporal discretizations, and illustrate their shortcomings in meeting maximum principles and the non-negative constraint. The proposed methodology can handle general computational grids with no additional restrictions on the time-step. We illustrate the performance and accuracy of the proposed methodology using representative numerical examples. We also perform a numerical convergence analysis of the proposed methodology. For comparison, we also present the results from the standard single-field semi-discrete formulation and the results from a popular software package, which all will violate maximum principles and the non-negative constraint.     

A NURBS-Enhanced Finite Volume Method for Steady Euler Equations with Goal-Oriented $h$-Adaptivity

In [A NURBS-enhanced finite volume solver for steady Euler equations, X. C. Meng, G. H. Hu, J. Comput. Phys., Vol. 359, pp. 77–92], a NURBS-enhanced finite volume method was developed to solve the steady Euler equations, in which the desired high order numerical accuracy was obtained for the equations imposed in the domain with a curved boundary. In this paper, the method is significantly improved in the following ways: (i) a simple and efficient point inversion technique is designed to compute the parameter values of points lying on a NURBS curve, (ii) with this new point inversion technique, the $h$-adaptive NURBS-enhanced finite volume method is introduced for the steady Euler equations in a complex domain, and (iii) a goal-oriented a posteriori error indicator is designed to further improve the efficiency of the algorithm towards accurately calculating a given quantity of interest. Numerical results obtained from a variety of numerical experiments with different flow configurations successfully show the effectiveness and robustness of the proposed method.     

A Stabilized Finite Element Method for Modified Poisson-Nernst-Planck Equations to Determine Ion Flow Through a Nanopore

The conventional Poisson-Nernst-Planck equations do not account for the finite size of ions explicitly. This leads to solutions featuring unrealistically high ionic concentrations in the regions subject to external potentials, in particular, near highly charged surfaces. A modified form of the Poisson-Nernst-Planck equations accounts for steric effects and results in solutions with finite ion concentrations. Here, we evaluate numerical methods for solving the modified Poisson-Nernst-Planck equations by modeling electric field-driven transport of ions through a nanopore. We describe a novel, robust finite element solver that combines the applications of the Newton's method to the nonlinear Galerkin form of the equations, augmented with stabilization terms to appropriately handle the drift-diffusion processes.To make direct comparison with particle-based simulations possible, our method is specifically designed to produce solutions under periodic boundary conditions and to conserve the number of ions in the solution domain. We test our finite element solver on a set of challenging numerical experiments that include calculations of the ion distribution in a volume confined between two charged plates, calculations of the ionic current though a nanopore subject to an external electric field, and modeling the effect of a DNA molecule on the ion concentration and nanopore current.     

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.     

Local Discontinuous-Galerkin Schemes for Model Problems in Phase Transition Theory

Local Discontinuous Galerkin (LDG) schemes in the sense of [5] are a flexible numerical tool to approximate solutions of nonlinear convection problems with complicated dissipative terms. Such terms frequently appear in evolution equations which describe the dynamics of phase changes in e.g. liquid-vapour mixtures or in elastic solids. We report on results for one-dimensional model problems with dissipative terms including third-order and convolution operators. Cell entropy inequalities and L²-stability results are proved for those model problems. As is common in phase transition theory the solution structure sensitively depends on the coupling parameter between viscosity and capillarity. To avoid spurious solutions due to the counteracting effect of artificial dissipation by the numerical flux and the actual dissipation terms we introduce Tadmors' entropy conservative fluxes. Various numerical experiments underline the reliability of our approach and also illustrate interesting and (partly) new phase transition phenomena.     

Do Current Lattice Boltzmann Methods for Diffusion and Advection-Diffusion Equations Respect Maximum Principle and the Non-Negative Constraint?

The Lattice Boltzmann Method (LBM) has established itself as a popular numerical method in computational fluid dynamics. Several advancements have been recently made in LBM, which include multiple-relaxation-time LBM to simulate anisotropic advection-diffusion processes. Because of the importance of LBM simulations for transport problems in subsurface and reactive flows, one needs to study the accuracy and structure preserving properties of numerical solutions under the LBM. The solutions to advective-diffusive systems are known to satisfy maximum principles, comparison principles, the non-negative constraint, and the decay property. In this paper, using several numerical experiments, it will be shown that current single- and multiple-relaxation-time lattice Boltzmann methods fail to preserve these mathematical properties for transient diffusion-type equations. We will also show that these violations may not be removed by simply refining the discretization parameters. More importantly, it will be shown that meeting stability conditions alone does not guarantee the preservation of the aforementioned mathematical principles and physical constraints in the discrete setting. A discussion on the source of these violations and possible approaches to avoid them is included. A condition to guarantee the non-negativity of concentration under LBM in the case of isotropic diffusion is also derived. The impact of this research is twofold. First, the study poses several outstanding research problems, which should guide researchers to develop LBM-based formulations for transport problems that respect important mathematical properties and physical constraints in the discrete setting. This paper can also serve as a good source of benchmark problems for such future research endeavors. Second, this study cautions the practitioners of the LBM for transport problems with the associated numerical deficiencies of the LBM, and provides guidelines for performing predictive simulations of advective-diffusive processes using the LBM.     

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.     

Weak Galerkin and Continuous Galerkin Coupled Finite Element Methods for the Stokes-Darcy Interface Problem

We consider a model of coupled free and porous media flow governed by Stokes equation and Darcy's law with the Beavers-Joseph-Saffman interface condition. In this paper, we propose a new numerical approach for the Stokes-Darcy system. The approach employs the classical finite element method for the Darcy region and the weak Galerkin finite element method for the Stokes region. We construct corresponding discrete scheme and prove its well-posedness. The estimates for the corresponding numerical approximation are derived. Finally, we present some numerical experiments to validate the efficiency of the approach for solving this problem.     

An Efficient Iterative Method for the Formulation of Flow Networks

We propose an efficient iterative convolution thresholding method for the formulation of flow networks where the fluid is modeled by the Darcy–Stokes flow with the presence of volume sources. The method is based on the minimization of the dissipation power in the fluid region with a Darcy term. The flow network is represented by its characteristic function and the energy is approximated under this representation. The minimization problem can then be approximately solved by alternating: 1) solving a Brinkman equation to model the Darcy–Stokes flow and 2) updating the characteristic function by a simple convolution and thresholding step. The proposed method is simple and easy to implement. We prove mathematically that the iterative method has the total energy decaying property. Numerical experiments demonstrate the performance and robustness of the proposed method and interesting structures are observed.     

Dissipative and Conservative Local Discontinuous Galerkin Methods for the Fornberg-Whitham Type Equations

In this paper, we construct high order energy dissipative and conservative local discontinuous Galerkin methods for the Fornberg-Whitham type equations. We give the proofs for the dissipation and conservation for related conservative quantities. The corresponding error estimates are proved for the proposed schemes. The capability of our schemes for different types of solutions is shown via several numerical experiments. The dissipative schemes have good behavior for shock solutions, while for a long time approximation, the conservative schemes can reduce the shape error and the decay of amplitude significantly.     

On the Importance of the Stokes-Brinkman Equations for Computing Effective Permeability in Karst Reservoirs

Cavities and fractures significantly affect the flow paths in carbonate reservoirs and should be accurately accounted for in numerical models. Herein, we consider the problem of computing the effective permeability of rock samples based on high-resolution 3D CT scans containing millions of voxels. We use the Stokes-Brinkman equations in the entire domain, covering regions of free flow governed by the Stokes equations, porous Darcy flow, and transitions between them. The presence of different length scales and large (ten orders of magnitude) contrasts in permeability leads to highly ill-conditioned linear systems of equations, which are difficult to solve. To obtain a problem that is computationally tractable, we first analyze the relative importance of the Stokes and Darcy terms for a set of idealized 2D models. We find that, in terms of effective permeability, the Stokes-Brinkman equations are only applicable for a special parameter set where the effective free-flow permeability is less than four orders of magnitude different from the matrix permeability. All other cases can be accurately modeled with either the Stokes or the Darcy end-member flows, depending on if there do or do not exist percolating free-flow regions. The insights obtained are used to perform a direct computation of the effective permeability of a rock sample model with more than 8 million cells.     

An Adaptive Moving Mesh Method for Two-Dimensional Incompressible Viscous Flows

In this paper, we present an adaptive moving mesh technique for solving the incompressible viscous flows using the vorticity stream-function formulation. The moving mesh strategy is based on the approach proposed by Li et al. [J. Comput. Phys., 170 (2001), pp. 562–588] to separate the mesh-moving and evolving PDE at each time step. The Navier-Stokes equations are solved in the vorticity stream-function form by a finite-volume method in space, and the mesh-moving part is realized by solving the Euler-Lagrange equations to minimize a certain variation in conjunction with a more sophisticated monitor function. A conservative interpolation is used to redistribute the numerical solutions on the new meshes. This paper discusses the implementation of the periodic boundary conditions, where the physical domain is allowed to deform with time while the computational domain remains fixed and regular throughout. Numerical results demonstrate the accuracy and effectiveness of the proposed algorithm.     

A New Weak Galerkin Finite Element Scheme for the Brinkman Model

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.     

VPVnet: A Velocity-Pressure-Vorticity Neural Network Method for the Stokes’ Equations under Reduced Regularity

We present VPVnet, a deep neural network method for the Stokes’ equations under reduced regularity. Different with recently proposed deep learning methods [40,51] which are based on the original form of PDEs, VPVnet uses the least square functional of the first-order velocity-pressure-vorticity (VPV) formulation ([30]) as loss functions. As such, only first-order derivative is required in the loss functions, hence the method is applicable to a much larger class of problems, e.g. problems with nonsmooth solutions. Despite that several methods have been proposed recently to reduce the regularity requirement by transforming the original problem into a corresponding variational form, while for the Stokes’ equations, the choice of approximating spaces for the velocity and the pressure has to satisfy the LBB condition additionally. Here by making use of the VPV formulation, lower regularity requirement is achieved with no need for considering the LBB condition. Convergence and error estimates have been established for the proposed method. It is worth emphasizing that the VPVnet method is divergence-free and pressure-robust, while classical inf-sup stable mixed finite elements for the Stokes’ equations are not pressure-robust. Various numerical experiments including 2D and 3D lid-driven cavity test cases are conducted to demonstrate its efficiency and accuracy.