Mechanics-Based Solution Verification for Porous Media Models

This paper presents a new approach to verify the accuracy of computational simulations. We develop mathematical theorems which can serve as robust a posteriori error estimation techniques to identify numerical pollution, check the performance of adaptive meshes, and verify numerical solutions. We demonstrate performance of this methodology on problems from flow thorough porous media. However, one can extend it to other models. We construct mathematical properties such that the solutions to Darcy and Darcy-Brinkman equations satisfy them. The mathematical properties include the total minimum mechanical power, minimum dissipation theorem, reciprocal relation, and maximum principle for the vorticity. All the developed theorems have firm mechanical bases and are independent of numerical methods. So, these can be utilized for solution verification of finite element, finite volume, finite difference, lattice Boltzmann methods and so forth. In particular, we show that, for a given set of boundary conditions, Darcy velocity has the minimum total mechanical power of all the kinematically admissible vector fields. We also show that a similar result holds for Darcy-Brinkman velocity. We then show for a conservative body force, the Darcy and Darcy-Brinkman velocities have the minimum total dissipation among their respective kinematically admissible vector fields. Using numerical examples, we show that the minimum dissipation and total mechanical power theorems can be utilized to identify pollution errors in numerical solutions. The solutions to Darcy and Darcy-Brinkman equations are shown to satisfy a reciprocal relation, which has the potential to identify errors in the numerical implementation of boundary conditions. It is also shown that the vorticity under both steady and transient Darcy-Brinkman equations satisfy maximum principles if the body force is conservative and the permeability is homogeneous and isotropic. A discussion on the nature of vorticity under steady and transient Darcy equations is also presented. Using several numerical examples, we will demonstrate the predictive capabilities of the proposed a posteriori techniques in assessing the accuracy of numerical solutions for a general class of problems, which could involve complex domains and general computational grids.     

A Strong Stability-Preserving Predictor-Corrector Method for the Simulation of Elastic Wave Propagation in Anisotropic Media

In this paper, we propose a strong stability-preserving predictor-corrector (SSPC) method based on an implicit Runge-Kutta method to solve the acoustic- and elastic-wave equations. We first transform the wave equations into a system of ordinary differential equations (ODEs) and apply the local extrapolation method to discretize the spatial high-order derivatives, resulting in a system of semi-discrete ODEs. Then we use the SSPC method based on an implicit Runge-Kutta method to solve the semi-discrete ODEs and introduce a weighting parameter into the SSPC method. On top of such a structure, we develop a robust numerical algorithm to effectively suppress the numerical dispersion, which is usually caused by the discretization of wave equations when coarse grids are used or geological models have large velocity contrasts between adjacent layers. Meanwhile, we investigate the performance of the SSPC method including numerical errors and convergence rate, numerical dispersion, and stability criteria with different choices of the weighting parameter to solve 1-D and 2-D acoustic- and elastic-wave equations. When the SSPC is applied to seismic simulations, the computational efficiency is also investigated by comparing the SSPC, the fourth-order Lax-Wendroff correction (LWC) method, and the staggered-grid (SG) finite difference method. Comparisons of synthetic waveforms computed by the SSPC and analytic solutions for acoustic and elastic models are given to illustrate the accuracy and the validity of the SSPC method. Furthermore, several numerical experiments are conducted for the geological models including a 2-D homogeneous transversely isotropic (TI) medium, a two-layer elastic model, and the 2-D SEG/EAGE salt model. The results show that the SSPC can be used as a practical tool for large-scale seismic simulation because of its effectiveness in suppressing numerical dispersion even in the situations such as coarse grids, strong interfaces, or high frequencies.     

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 New Approach to Implement Sigma Coordinate in a Numerical Model

This study shows a new way to implement terrain-following σ-coordinate in a numerical model, which does not lead to the well-known "pressure gradient force (PGF)" problem. First, the causes of the PGF problem are analyzed with existing methods that are categorized into two different types based on the causes. Then, the new method that bypasses the PGF problem all together is proposed. By comparing these three methods and analyzing the expression of the scalar gradient in a curvilinear coordinate system, this study finds out that only when using the covariant scalar equations of σ-coordinate will the PGF computational form have one term in each momentum component equation, thereby avoiding the PGF problem completely. A convenient way of implementing the covariant scalar equations of σ-coordinate in a numerical atmospheric model is illustrated, which is to set corresponding parameters in the scalar equations of the Cartesian coordinate. Finally, two idealized experiments manifest that the PGF calculated with the new method is more accurate than using the classic one. This method can be used for oceanic models as well, and needs to be tested in both the atmospheric and oceanic models.     

A Hodge Decomposition Method for Dynamic Ginzburg–Landau Equations in Nonsmooth Domains — A Second Approach

In a general polygonal domain, possibly nonconvex and multi-connected (with holes), the time-dependent Ginzburg–Landau equation is reformulated into a new system of equations. The magnetic field $B$:=∇×A is introduced as an unknown solution in the new system, while the magnetic potential A is solved implicitly through its Hodge decomposition into divergence-free part, curl-free and harmonic parts, separately. Global well-posedness of the new system and its equivalence to the original problem are proved. A linearized and decoupled Galerkin finite element method is proposed for solving the new system. The convergence of numerical solutions is proved based on a compactness argument by utilizing the maximal $L^p$-regularity of the discretized equations. Compared with the Hodge decomposition method proposed in [27]_,the new method has the advantage of approximating the magnetic field B directly and converging for initial conditions that are incompatible with the external magnetic field. Several numerical examples are provided to illustrate the efficiency of the proposed numerical method in both simply connected and multi-connected nonsmooth domains. We observe that even in simply connected domains, the new method is superior to the method in [27] for approximating the magnetic field.     

Truncation Errors,Exact and Heuristic Stability Analysis of Two-Relaxation-Times Lattice Boltzmann Schemes for Anisotropic Advection-Diffusion Equation

This paper establishes relations between the stability and the high-order truncated corrections for modeling of the mass conservation equation with the two-relaxation-times (TRT) collision operator. First we propose a simple method to derive the truncation errors from the exact, central-difference type, recurrence equations of the TRT scheme. They also supply its equivalent three-time-level discretization form. Two different relationships of the two relaxation rates nullify the third (advection) and fourth (pure diffusion) truncation errors, for any linear equilibrium and any velocity set. However, the two relaxation times alone cannot remove the leading-order advection-diffusion error, because of the intrinsic fourth-order numerical diffusion. The truncation analysis is carefully verified for the evolution of concentration waves with the anisotropic diffusion tensors. The anisotropic equilibrium functions are presented in a simple but general form, suitable for the minimal velocity sets and the d2Q9, d3Q13, d3Q15 and d3Q19 velocity sets. All anisotropic schemes are complemented by their exact necessary von Neumann stability conditions and equivalent finite-difference stencils. The sufficient stability conditions are proposed for the most stable (OTRT) family, which enables modeling at any Peclet numbers with the same velocity amplitude. The heuristic stability analysis of the fourth-order truncated corrections extends the optimal stability to larger relationships of the two relaxation rates, in agreement with the exact (one-dimensional) and numerical (multi-dimensional) stability analysis. A special attention is put on the choice of the equilibrium weights. By combining accuracy and stability predictions, several strategies for selecting the relaxation and free-tunable equilibrium parameters are suggested and applied to the evolution of the Gaussian hill.     

Laplace-Transform Finite Element Solution of Nonlocal and Localized Stochastic Moment Equations of Transport

Morales-Casique et al. (Adv. Water Res., 29 (2006), pp. 1238-1255) developed exact first and second nonlocal moment equations for advective-dispersive transport in finite, randomly heterogeneous geologic media. The velocity and concentration in these equations are generally nonstationary due to trends in heterogeneity, conditioning on site data and the influence of forcing terms. Morales-Casique et al. (Adv. Water Res., 29 (2006), pp. 1399-1418) solved the Laplace transformed versions of these equations recursively to second order in the standard deviation σY of (natural) log hydraulic conductivity, and iteratively to higher-order, by finite elements followed by numerical inversion of the Laplace transform. They did the same for a space-localized version of the mean transport equation. Here we recount briefly their theory and algorithms; compare the numerical performance of the Laplace-transform finite element scheme with that of a high-accuracy ULTIMATE-QUICKEST algorithm coupled with an alternating split operator approach; and review some computational results due to Morales-Casique et al. (Adv. Water Res., 29 (2006), pp. 1399-1418) to shed light on the accuracy and computational efficiency of their recursive and iterative solutions in comparison to conditional Monte Carlo simulations in two spatial dimensions.     

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 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.     

Solution of ill-conditioned optimality conditions for control of turbo-alternators

The computational solution of the optimality equations for torque and voltage control of a large-scale turbo-alternator is considered. The equations are extremely ill-conditioned as a consequence of the ill-conditioning of the linearized system model. Our main result is to point out that using the Ricatti equation to estimate the initial values of the Lagrange multipliers renders the equations in a form that can easily be solved as an initial value problem for the non-linear system.     

A Gas-Kinetic Unified Algorithm for Non-Equilibrium Polyatomic Gas Flows Covering Various Flow Regimes

In this paper, a gas-kinetic unified algorithm (GKUA) is developed to investigate the non-equilibrium polyatomic gas flows covering various regimes. Based on the ellipsoidal statistical model with rotational energy excitation, the computable modelling equation is presented by unifying expressions on the molecular collision relaxing parameter and the local equilibrium distribution function. By constructing the corresponding conservative discrete velocity ordinate method for this model, the conservative properties during the collision procedure are preserved at the discrete level by the numerical method, decreasing the computational storage and time. Explicit and implicit lower-upper symmetric Gauss-Seidel schemes are constructed to solve the discrete hyperbolic conservation equations directly. Applying the new GKUA, some numerical examples are simulated, including the Sod Riemann problem, homogeneous flow rotational relaxation, normal shock structure, Fourier and Couette flows, supersonic flows past a circular cylinder, and hypersonic flow around a plate placed normally. The results obtained by the analytic, experimental, direct simulation Monte Carlo method, and other measurements in references are compared with the GKUA results, which are in good agreement, demonstrating the high accuracy of the present algorithm. Especially, some polyatomic gas non-equilibrium phenomena are observed and analysed by solving the Boltzmann-type velocity distribution function equation covering various flow regimes.     

The Hamiltonian Field Theory of the Von Mises Wave Equation: Analytical and Computational Issues

The Von Mises quasi-linear second order wave equation, which completely describes an irrotational, compressible and barotropic classical perfect fluid, can be derived from a nontrivial least action principle for the velocity scalar potential only, in contrast to existing analog formulations which are expressed in terms of coupled density and velocity fields. In this article, the classical Hamiltonian field theory specifically associated to such an equation is developed in the polytropic case and numerically verified in a simplified situation. The existence of such a mathematical structure suggests new theoretical schemes possibly useful for performing numerical integrations of fluid dynamical equations. Moreover, it justifies possible new functional forms for Lagrangian densities and associated Hamiltonian functions in other theoretical classical physics contexts.     

High-Order Gas-Kinetic Scheme in Curvilinear Coordinates for the Euler and Navier-Stokes Solutions

The high-order gas-kinetic scheme (HGKS) has achieved success in simulating compressible flows with Cartesian meshes. To study the flow problems in general geometries, such as the flow over a wing-body, the development of HGKS in general curvilinear coordinates becomes necessary. In this paper, a two-stage fourth-order gas-kinetic scheme is developed for the Euler and Navier-Stokes solutions in the curvilinear coordinates from one-dimensional to three-dimensional computations. Based on the coordinate transformation, the kinetic equation is transformed first to the computational space, and the flux function in the gas-kinetic scheme is obtained there and is transformed back to the physical domain for the update of flow variables inside each control volume. To achieve the expected order of accuracy, the dimension-by-dimension reconstruction based on the WENO scheme is adopted in the computational domain, where the reconstructed variables are the cell averaged Jacobian and the Jacobian-weighted conservative variables. In the two-stage fourth-order gas-kinetic scheme, the point values as well as the spatial derivatives of conservative variables at Gaussian quadrature points have to be used in the evaluation of the time dependent flux function. The point-wise conservative variables are obtained by ratio of the above reconstructed data, and the spatial derivatives are reconstructed through orthogonalization in physical space and chain rule. A variety of numerical examples from the accuracy tests to the solutions with strong discontinuities are presented to validate the accuracy and robustness of the current scheme for both inviscid and viscous flows. The precise satisfaction of the geometrical conservation law in non-orthogonal mesh is also demonstrated through the numerical example.     

Three Discontinuous Galerkin Methods for One- and Two-Dimensional Nonlinear Dirac Equations with a Scalar Self-Interaction

This paper develops three high-order accurate discontinuous Galerkin (DG) methods for the one-dimensional (1D) and two-dimensional (2D) nonlinear Dirac (NLD) equations with a general scalar self-interaction. They are the Runge-Kutta DG (RKDG) method and the DG methods with the one-stage fourth-order Lax-Wendroff type time discretization (LWDG) and the two-stage fourth-order accurate time discretization (TSDG). The RKDG method uses the spatial DG approximation to discretize the NLD equations and then utilize the explicit multistage high-order Runge-Kutta time discretization for the first-order time derivatives, while the LWDG and TSDG methods, on the contrary, first give the one-stage fourth-order Lax-Wendroff type and the two-stage fourth-order time discretizations of the NLD equations, respectively, and then discretize the first- and higher-order spatial derivatives by using the spatial DG approximation. The $L^2$ stability of the 2D semi-discrete DG approximation is proved in the RKDG methods for a general triangulation, and the computational complexities of three 1D DG methods are estimated. Numerical experiments are conducted to validate the accuracy and the conservation properties of the proposed methods. The interactions of the solitary waves, the standing and travelling waves are investigated numerically and the 2D breathing pattern is observed.     

A High-Order Accurate Gas-Kinetic Scheme for One- and Two-Dimensional Flow Simulation

This paper develops a high-order accurate gas-kinetic scheme in the framework of the finite volume method for the one- and two-dimensional flow simulations, which is an extension of the third-order accurate gas-kinetic scheme [Q.B. Li, K. Xu, and S. Fu, J. Comput. Phys., 229(2010), 6715-6731] and the second-order accurate gas-kinetic scheme [K. Xu, J. Comput. Phys., 171(2001), 289-335]. It is formed by two parts: quartic polynomial reconstruction of the macroscopic variables and fourth-order accurate flux evolution. The first part reconstructs a piecewise cell-center based quartic polynomial and a cell-vertex based quartic polynomial according to the "initial" cell average approximation of macroscopic variables to recover locally the non-equilibrium and equilibrium single particle velocity distribution functions around the cell interface. It is in view of the fact that all macroscopic variables become moments of a single particle velocity distribution function in the gas-kinetic theory. The generalized moment limiter is employed there to suppress the possible numerical oscillation. In the second part, the macroscopic flux at the cell interface is evolved in fourth-order accuracy by means of the simple particle transport mechanism in the microscopic level, i.e. free transport and the Bhatnagar-Gross-Krook (BGK) collisions. In other words, the fourth-order flux evolution is based on the solution (i.e. the particle velocity distribution function) of the BGK model for the Boltzmann equation. Several 1D and 2D test problems are numerically solved by using the proposed high-order accurate gas-kinetic scheme. By comparing with the exact solutions or the numerical solutions obtained the second-order or third-order accurate gas-kinetic scheme, the computations demonstrate that our scheme is effective and accurate for simulating invisid and viscous fluid flows, and the accuracy of the high-order GKS depends on the choice of the (numerical) collision time.     

Analysis and Numerical Simulation of Hyperbolic Shallow Water Moment Equations

Shallow Water Moment Equations allow for vertical changes in the horizontal velocity, so that complex shallow flows can be described accurately. However, we show that these models lack global hyperbolicity and the loss of hyperbolicity already occurs for small deviations from equilibrium. This leads to instabilities in a numerical test case. We then derive new Hyperbolic Shallow Water Moment Equations based on a modification of the system matrix. The model can be written in analytical form and hyperbolicity can be proven for a large number of equations. A second variant of this model is obtained by generalizing the modification with the help of additional parameters. Numerical tests of a smooth periodic problem and a dam break problem using the new models yield accurate and fast solutions while guaranteeing hyperbolicity.     

Geometric Numerical Integration for Peakon b-Family Equations

In this paper, we study the Camassa-Holm equation and the Degasperis-Procesi equation. The two equations are in the family of integrable peakon equations, and both have very rich geometric properties. Based on these geometric structures, we construct the geometric numerical integrators for simulating their soliton solutions. The Camassa-Holm equation and the Degasperis-Procesi equation have many common properties, however, they also have the significant differences, for example, there exist the shock wave solutions for the Degasperis-Procesi equation. By using the symplectic Fourier pseudo-spectral integrator, we simulate the peakon solutions of the two equations. To illustrate the smooth solitons and shock wave solutions of the DP equation, we use the splitting technique and combine the composition methods. In the numerical experiments, comparisons of these two kinds of methods are presented in terms of accuracy, computational cost and invariants preservation.     

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.     

A Discrete Flux Scheme for Aerodynamic and Hydrodynamic Flows

The objective of this paper is to seek an alternative to the numerical simulation of the Navier-Stokes equations by a method similar to solving the BGK-type modeled lattice Boltzmann equation. The proposed method is valid for both gas and liquid flows. A discrete flux scheme (DFS) is used to derive the governing equations for two distribution functions; one for mass and another for thermal energy. These equations are derived by considering an infinitesimally small control volume with a velocity lattice representation for the distribution functions. The zero-order moment equation of the mass distribution function is used to recover the continuity equation, while the first-order moment equation recovers the linear momentum equation. The recovered equations are correct to the first order of the Knudsen number (Kn); thus, satisfying the continuum assumption. Similarly, the zero-order moment equation of the thermal energy distribution function is used to recover the thermal energy equation. For aerodynamic flows, it is shown that the finite difference solution of the DFS is equivalent to solving the lattice Boltzmann equation (LBE) with a BGK-type model and a specified equation of state. Thus formulated, the DFS can be used to simulate a variety of aerodynamic and hydrodynamic flows. Examples of classical aeroacoustics, compressible flow with shocks, incompressible isothermal and non-isothermal Couette flows, stratified flow in a cavity, and double diffusive flow inside a rectangle are used to demonstrate the validity and extent of the DFS. Very good to excellent agreement with known analytical and/or numerical solutions is obtained; thus lending evidence to the DFS approach as an alternative to solving the Navier-Stokes equations for fluid flow simulations.     

Some Mathematical and Numerical Issues in Geophysical Fluid Dynamics and Climate Dynamics

In this article, we address both recent advances and open questions in some mathematical and computational issues in geophysical fluid dynamics (GFD) and climate dynamics. The main focus is on 1) the primitive equations (PEs) models and their related mathematical and computational issues, 2) climate variability, predictability and successive bifurcation, and 3) a new dynamical systems theory and its applications to GFD and climate dynamics.