A Fourth-Order Upwinding Embedded Boundary Method (UEBM) for Maxwell's Equations in Media with Material Interfaces: Part I

In this paper, we present a new fourth-order upwinding embedded boundary method (UEBM) over Cartesian grids, originally proposed in the Journal of Computational Physics [190 (2003), pp. 159-183.] as a second-order method for treating material interfaces for Maxwell's equations. In addition to the idea of the UEBM to evolve solutions at interfaces, we utilize the ghost fluid method to construct finite difference approximation of spatial derivatives at Cartesian grid points near the material interfaces. As a result, Runge-Kutta type time discretization can be used for the semidiscretized system to yield an overall fourth-order method, in contrast to the original second-order UEBM based on a Lax-Wendroff type difference. The final scheme allows time step sizes independent of the interface locations. Numerical examples are given to demonstrate the fourth-order accuracy as well as the stability of the method. We tested the scheme for several wave problems with various material interface locations, including electromagnetic scattering of a plane wave incident on a planar boundary and a two-dimensional electromagnetic application with an interface parallel to the y-axis.     

Exponential Compact Higher Order Scheme for Nonlinear Steady Convection-Diffusion Equations

This paper presents an exponential compact higher order scheme for Convection-Diffusion Equations (CDE) with variable and nonlinear convection coefficients. The scheme is O(h⁴) for one-dimensional problems and produces a tri-diagonal system of equations which can be solved efficiently using Thomas algorithm. For two-dimensional problems, the scheme produces an O(h⁴+k⁴) accuracy over a compact nine point stencil which can be solved using any line iterative approach with alternate direction implicit procedure. The convergence of the iterative procedure is guaranteed as the coefficient matrix of the developed scheme satisfies the conditions required to be positive. Wave number analysis has been carried out to establish that the scheme is comparable in accuracy with spectral methods. The higher order accuracy and better rate of convergence of the developed scheme have been demonstrated by solving numerous model problems for one- and two-dimensional CDE, where the solutions have the sharp gradient at the solution boundary.     

Development of a Combined Compact Difference Scheme to Simulate Soliton Collision in a Shallow Water Equation

In this paper a three-step scheme is applied to solve the Camassa-Holm (CH) shallow water equation. The differential order of the CH equation has been reduced in order to facilitate development of numerical scheme in a comparatively smaller grid stencil. Here a three-point seventh-order spatially accurate upwinding combined compact difference (CCD) scheme is proposed to approximate the first-order derivative term. We conduct modified equation analysis on the CCD scheme and eliminate the leading discretization error terms for accurately predicting unidirectional wave propagation. The Fourier analysis is carried out as well on the proposed numerical scheme to minimize the dispersive error. For preserving Hamiltonians in Camassa-Holm equation, a symplecticity conserving time integrator has been employed. The other main emphasis of the present study is the use of u−P−α formulation to get nondissipative CH solution for peakon-antipeakon and soliton-anticuspon head-on wave collision problems.     

Simulation of Incompressible Free Surface Flow Using the Volume Preserving Level Set Method

This study aims to develop a numerical scheme in collocated Cartesian grids to solve the level set equation together with the incompressible two-phase flow equations. A seventh-order accurate upwinding combined compact difference (UCCD7) scheme has been developed for the approximation of the first-order spatial derivative terms shown in the level set equation. Developed scheme has a higher accuracy with a three-point grid stencil to minimize phase error. To preserve the mass of each phase all the time, the temporal derivative term in the level set equation is approximated by the sixth-order accurate symplectic Runge-Kutta (SRK6) scheme. All the simulated results for the dam-break, Rayleigh-Taylor instability, bubble rising, two-bubble merging, and milkcrown problems in two and three dimensions agree well with the available numerical or experimental results.     

Continuum Formulation for Non-Equilibrium Shock Structure Calculation

Are extensions to continuum formulations for solving fluid dynamic problems in the transition-to-rarefied regimes viable alternatives to particle methods? It is well known that for increasingly rarefied flow fields, the predictions from continuum formulation, such as the Navier-Stokes equations lose accuracy. These inaccuracies are attributed primarily to the linear approximations of the stress and heat flux terms in the Navier-Stokes equations. The inclusion of higher-order terms, such as Burnett or high-order moment equations, could improve the predictive capabilities of such continuum formulations, but there has been limited success in the shock structure calculations, especially for the high Mach number case. Here, after reformulating the viscosity and heat conduction coefficients appropriate for the rarefied flow regime, we will show that the Navier-Stokes-type continuum formulation may still be properly used. The equations with generalization of the dissipative coefficients based on the closed solution of the Bhatnagar-Gross-Krook (BGK) model of the Boltzmann equation, are solved using the gas-kinetic numerical scheme. This paper concentrates on the non-equilibrium shock structure calculations for both monatomic and diatomic gases. The Landau-Teller-Jeans relaxation model for the rotational energy is used to evaluate the quantitative difference between the translational and rotational temperatures inside the shock layer. Variations of shear stress, heat flux, temperatures, and densities in the internal structure of the shock waves are compared with, (a) existing theoretical solutions of the Boltzmann solution, (b) existing numerical predictions of the direct simulation Monte Carlo (DSMC) method, and (c) available experimental measurements. The present continuum formulation for calculating the shock structures for monatomic and diatomic gases in the Mach number range of 1.2 to 12.9 is found to be satisfactory.     

A Stability Analysis of Hybrid Schemes to Cure Shock Instability

The carbuncle phenomenon has been regarded as a spurious solution produced by most of contact-preserving methods. The hybrid method of combining high resolution flux with more dissipative solver is an attractive attempt to cure this kind of non-physical phenomenon. In this paper, a matrix-based stability analysis for 2-D Euler equations is performed to explore the cause of instability of numerical schemes. By combining the Roe with HLL flux in different directions and different flux components, we give an interesting explanation to the linear numerical instability. Based on such analysis, some hybrid schemes are compared to illustrate different mechanisms in controlling shock instability. Numerical experiments are presented to verify our analysis results. The conclusion is that the scheme of restricting directly instability source is more stable than other hybrid schemes.     

Multiple-Temperature Gas-Kinetic Scheme for Type IV Shock/Shock Interaction

In this paper, a gas-kinetic scheme (GKS) method coupled with a three temperature kinetic model is presented and applied in numerical study of the Edney-type IV shock/shock interaction which could cause serious problems in hypersonic vehicles. The results showed very good agreement with the experimental data in predicting the heat flux on the surface. It could be obviously seen that the current method can accurately describe the position and features of supersonic jets structure and clearly capture the thermal non-equilibrium in this case. The three temperature kinetic model includes three different models of temperatures which are translational, rotational and vibrational temperatures. The thermal non-equilibrium model is used to better simulate the aerodynamic and thermodynamic phenomenon. Current results were compared with the experimental data, computational fluid dynamics (CFD) results, and the Direct Simulation Monte Carlo (DSMC) results. Other CFD methods include the original GKS method without considering thermal non-equilibrium, the GKS method with a two temperature kinetic model and the Navier-Stokes equations with a three temperature kinetic model, which is the same as the multiple temperature kinetic model in current GKS method. Comparisons were made for the surface heat flux, pressure loads, Mach number contours and flow field values, including rotational temperature and density. By Comparing with other CFD method, the current GKS method showed a lot of improvement in predicting the magnitude and position of heat flux peak on the surface. This demonstrated the good potential of the current GKS method in solving thermodynamic non-equilibrium problems in hypersonic flows. The good performance of predicting the heat flux could bring a lot of benefit for the designing of the thermal protection system (TPS) for the hypersonic vehicles. By comparing with the original GKS method and the two temperature kinetic model, the three temperature kinetic model showed its importance and accuracy in this case.     

A Compact Third-Order Gas-Kinetic Scheme for Compressible Euler and Navier-Stokes Equations

In this paper, a compact third-order gas-kinetic scheme is proposed for the compressible Euler and Navier-Stokes equations. The main reason for the feasibility to develop such a high-order scheme with compact stencil, which involves only neighboring cells, is due to the use of a high-order gas evolution model. Besides the evaluation of the time-dependent flux function across a cell interface, the high-order gas evolution model also provides an accurate time-dependent solution of the flow variables at a cell interface. Therefore, the current scheme not only updates the cell averaged conservative flow variables inside each control volume, but also tracks the flow variables at the cell interface at the next time level. As a result, with both cell averaged and cell interface values, the high-order reconstruction in the current scheme can be done compactly. Different from using a weak formulation for high-order accuracy in the Discontinuous Galerkin method, the current scheme is based on the strong solution, where the flow evolution starting from a piecewise discontinuous high-order initial data is precisely followed. The cell interface time-dependent flow variables can be used for the initial data reconstruction at the beginning of next time step. Even with compact stencil, the current scheme has third-order accuracy in the smooth flow regions, and has favorable shock capturing property in the discontinuous regions. It can be faithfully used from the incompressible limit to the hypersonic flow computations, and many test cases are used to validate the current scheme. In comparison with many other high-order schemes, the current method avoids the use of Gaussian points for the flux evaluation along the cell interface and the multi-stage Runge-Kutta time stepping technique. Due to its multidimensional property of including both derivatives of flow variables in the normal and tangential directions of a cell interface, the viscous flow solution, especially those with vortex structure, can be accurately captured. With the same stencil of a second order scheme, numerical tests demonstrate that the current scheme is as robust as well-developed second-order shock capturing schemes, but provides more accurate numerical solutions than the second order counterparts.     

High Order Numerical Simulation of Detonation Wave Propagation Through Complex Obstacles with the Inverse Lax-Wendroff Treatment

The high order inverse Lax-Wendroff (ILW) procedure is extended to boundary treatment involving complex geometries on a Cartesian mesh. Our method ensures that the numerical resolution at the vicinity of the boundary and the inner domain keeps the fifth order accuracy for the system of the reactive Euler equations with the two-step reaction model. Shock wave propagation in a tube with an array of rectangular grooves is first numerically simulated by combining a fifth order weighted essentially non-oscillatory (WENO) scheme and the ILW boundary treatment. Compared with the experimental results, the ILW treatment accurately captures the evolution of shock wave during the interactions of the shock waves with the complex obstacles. Excellent agreement between our numerical results and the experimental ones further demonstrates the reliability and accuracy of the ILW treatment. Compared with the immersed boundary method (IBM), it is clear that the influence on pressure peaks in the reflected zone is obviously bigger than that in the diffracted zone. Furthermore, we also simulate the propagation process of detonation wave in a tube with three different widths of wall-mounted rectangular obstacles located on the lower wall. It is shown that the shock pressure along a horizontal line near the rectangular obstacles gradually decreases, and the detonation cellular size becomes large and irregular with the decrease of the obstacle width.     

Accuracy of the Adaptive GRP Scheme and the Simulation of 2-D Riemann Problems for Compressible Euler Equations

The adaptive generalized Riemann problem (GRP) scheme for 2-D compressible fluid flows has been proposed in [J. Comput. Phys., 229 (2010), 1448–1466] and it displays the capability in overcoming difficulties such as the start-up error for a single shock, and the numerical instability of the almost stationary shock. In this paper, we will provide the accuracy study and particularly show the performance in simulating 2-D complex wave configurations formulated with the 2-D Riemann problems for compressible Euler equations. For this purpose, we will first review the GRP scheme briefly when combined with the adaptive moving mesh technique and consider the accuracy of the adaptive GRP scheme via the comparison with the explicit formulae of analytic solutions of planar rarefaction waves, planar shock waves, the collapse problem of a wedge-shaped dam and the spiral formation problem. Then we simulate the full set of wave configurations in the 2-D four-wave Riemann problems for compressible Euler equations [SIAM J. Math. Anal., 21 (1990), 593–630], including the interactions of strong shocks (shock reflections), vortex-vortex and shock-vortex etc. This study combines the theoretical results with the numerical simulations, and thus demonstrates what Ami Harten observed "for computational scientists there are two kinds of truth: the truth that you prove, and the truth you see when you compute" [J. Sci. Comput., 31 (2007), 185–193].     

High Order Compact Schemes in Projection Methods for Incompressible Viscous Flows

Within the projection schemes for the incompressible Navier-Stokes equations (namely "pressure-correction" method), we consider the simplest method (of order one in time) which takes into account the pressure in both steps of the splitting scheme. For this scheme, we construct, analyze and implement a new high order compact spatial approximation on nonstaggered grids. This approach yields a fourth order accuracy in space with an optimal treatment of the boundary conditions (without error on the velocity) which could be extended to more general splitting. We prove the unconditional stability of the associated Cauchy problem via von Neumann analysis. Then we carry out a normal mode analysis so as to obtain more precise results about the behavior of the numerical solutions. Finally we present detailed numerical tests for the Stokes and the Navier-Stokes equations (including the driven cavity benchmark) to illustrate the theoretical results.     

Upwind Biased Local RBF Scheme with PDE Centres for the Steady Convection Diffusion Equations with Continuous and Discontinuous Boundary Conditions

RBF based grid-free scheme with PDE centres is experimented in this work for solving Convection-Diffusion Equations (CDE), a simplified model of the Navier-Stokes equations. For convection dominated problems, very few integration schemes can give converged solutions for the entire range of diffusivity wherein sharp layers are expected in the solutions and accurate computation of these layers is a big challenge for most of the numerical schemes. Radial Basis Function (RBF) based Local Hermitian Interpolation (LHI) with PDE centres is one such integration scheme which has some built in upwind effect and hence may be a good solver for the convection dominated problems. In the present work, to get convergent solutions consistently for small diffusion parameters, an explicit upwinding is also introduced in to the RBF based scheme with PDE centres, which was initially used to solve some time dependent problems in [10]. RBF based numerical schemes are one type of grid free numerical schemes based on the radial distances and hence very easy to use in high dimensional problems. In this work, the RBF scheme, with different upwind biasing, is used to a variety of steady benchmark problems with continuous and discontinuous boundary data and validated against the corresponding exact solutions. Comparisons of the solutions of the convective dominant benchmark problems show that the upwind biasing either in source centres or PDE centres gives convergent solutions consistently and is stable without any oscillations especially for problems with discontinuities in the boundary conditions. It is observed that the accuracy of the solutions is better than the solutions of other standard integration schemes particularly when the computations are carried out with fewer centers.     

Higher-Order Compact Scheme for the Incompressible Navier-Stokes Equations in Spherical Geometry

A higher-order compact scheme on the nine point 2-D stencil is developed for the steady stream-function vorticity form of the incompressible Navier-Stokes (NS) equations in spherical polar coordinates, which was used earlier only for the cartesian and cylindrical geometries. The steady, incompressible, viscous and axially symmetric flow past a sphere is used as a model problem. The non-linearity in the N-S equations is handled in a comprehensive manner avoiding complications in calculations. The scheme is combined with the multigrid method to enhance the convergence rate. The solutions are obtained over a non-uniform grid generated using the transformation r = e^ξ while maintaining a uniform grid in the computational plane. The superiority of the higher order compact scheme is clearly illustrated in comparison with upwind scheme and defect correction technique at high Reynolds numbers by taking a large domain. This is a pioneering effort, because for the first time, the fourth order accurate solutions for the problem of viscous flow past a sphere are presented here. The drag coefficient and surface pressures are calculated and compared with available experimental and theoretical results. It is observed that these values simulated over coarser grids using the present scheme are more accurate when compared to other conventional schemes. It has also been observed that the flow separation initially occurred at Re=21.     

Runge-Kutta Discontinuous Galerkin Method Using WENO-Type Limiters: Three-Dimensional Unstructured Meshes

This paper further considers weighted essentially non-oscillatory (WENO) and Hermite weighted essentially non-oscillatory (HWENO) finite volume methods as limiters for Runge-Kutta discontinuous Galerkin (RKDG) methods to solve problems involving nonlinear hyperbolic conservation laws. The application discussed here is the solution of 3-D problems on unstructured meshes. Our numerical tests again demonstrate this is a robust and high order limiting procedure, which simultaneously achieves high order accuracy and sharp non-oscillatory shock transitions.     

Two Uniform Tailored Finite Point Schemes for the Two Dimensional Discrete Ordinates Transport Equations with Boundary and Interface Layers

This paper presents two uniformly convergent numerical schemes for the two dimensional steady state discrete ordinates transport equation in the diffusive regime, which is valid up to the boundary and interface layers. A five-point node-centered and a four-point cell-centered tailored finite point schemes (TFPS) are introduced. The schemes first approximate the scattering coefficients and sources by piecewise constant functions and then use special solutions to the constant coefficient equation as local basis functions to formulate a discrete linear system. Numerically, both methods can not only capture the diffusion limit, but also exhibit uniform convergence in the diffusive regime, even with boundary layers. Numerical results show that the five-point scheme has first-order accuracy and the four-point scheme has second-order accuracy, uniformly with respect to the mean free path. Therefore, a relatively coarse grid can be used to capture the two dimensional boundary and interface layers.     

Numerical Simulation of Free Surface by an Area-Preserving Level Set Method

We apply in this study an area preserving level set method to simulate gas/water interface flow. For the sake of accuracy, the spatial derivative terms in the equations of motion for an incompressible fluid flow are approximated by the fifth-order accurate upwinding combined compact difference (UCCD) scheme. This scheme development employs two coupled equations to calculate the first- and second-order derivative terms in the momentum equations. For accurately predicting the level set value, the interface tracking scheme is also developed to minimize phase error of the first-order derivative term shown in the pure advection equation. For the purpose of retaining the long-term accurate Hamiltonian in the advection equation for the level set function, the time derivative term is discretized by the sixth-order accurate symplectic Runge-Kutta scheme. Also, to keep as a distance function for ensuring the front having a finite thickness for all time, the re-initialization equation is used. For the verification of the optimized UCCD scheme for the pure advection equation, two benchmark problems have been chosen to investigate in this study. The level set method with excellent area conservation property proposed for capturing the interface in incompressible fluid flows is also verified by solving the dam-break, Rayleigh-Taylor instability, two-bubble rising in water, and droplet falling problems.     

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.     

A Space-Time Conservative Method for Hyperbolic Systems with Stiff and Non-Stiff Source Terms

In this article we propose a higher-order space-time conservative method for hyperbolic systems with stiff and non-stiff source terms as well as relaxation systems. We call the scheme a slope propagation (SP) method. It is an extension of our scheme derived for homogeneous hyperbolic systems [1]. In the present inhomogeneous systems the relaxation time may vary from order of one to a very small value. These small values make the relaxation term stronger and highly stiff. In such situations underresolved numerical schemes may produce spurious numerical results. However, our present scheme has the capability to correctly capture the behavior of the physical phenomena with high order accuracy even if the initial layer and the small relaxation time are not numerically resolved. The scheme treats the space and time in a unified manner. The flow variables and their slopes are the basic unknowns in the scheme. The source term is treated by its volumetric integration over the space-time control volume and is a direct part of the overall space-time flux balance. We use two approaches for the slope calculations of the flow variables, the first one results directly from the flux balance over the control volumes, while in the second one we use a finite difference approach. The main features of the scheme are its simplicity, its Jacobian-free and Riemann solver-free recipe, as well as its efficiency and high order accuracy. In particular we show that the scheme has a discrete analog of the continuous asymptotic limit. We have implemented our scheme for various test models available in the literature such as the Broadwell model, the extended thermodynamics equations, the shallow water equations, traffic flow and the Euler equations with heat transfer. The numerical results validate the accuracy, versatility and robustness of the present scheme.     

A Novel Technique for Constructing Difference Schemes for Systems of Singularly Perturbed Equations

In this paper, we propose a novel and simple technique to construct effective difference schemes for solving systems of singularly perturbed convection-diffusion-reaction equations, whose solutions may display boundary or interior layers. We illustrate the technique by taking the Il'in-Allen-Southwell scheme for 1-D scalar equations as a basis to derive a formally second-order scheme for 1-D coupled systems and then extend the scheme to 2-D case by employing an alternating direction approach. Numerical examples are given to demonstrate the high performance of the obtained scheme on uniform meshes as well as piecewise-uniform Shishkin meshes.     

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

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