Remapping-Free Adaptive GRP Method for Multi-Fluid Flows I: One Dimensional Euler Equations

In this paper, a remapping-free adaptive GRP method for one dimensional (1-D) compressible flows is developed. Based on the framework of finite volume method, the 1-D Euler equations are discretized on moving volumes and the resulting numerical fluxes are computed directly by the GRP method. Thus the remapping process in the earlier adaptive GRP algorithm [17,18] is omitted. By adopting a flexible moving mesh strategy, this method could be applied for multi-fluid problems. The interface of two fluids will be kept at the node of computational grids and the GRP solver is extended at the material interfaces of multi-fluid flows accordingly. Some typical numerical tests show competitive performances of the new method, especially for contact discontinuities of one fluid cases and the material interface tracking of multi-fluid cases.     

Piecewise Polynomial Mapping Method and Corresponding WENO Scheme with Improved Resolution

The method of mapping function was first proposed by Henrick et al. [J. Comput. Phys. 207:542-547 (2005)] to adjust nonlinear weights in [0,1] for the fifth-order WENO scheme, and through which the requirement of convergence order is satisfied and the performance of the scheme is improved. Different from Henrick's method, a concept of piecewise polynomial function is proposed in this study and corresponding WENO schemes are obtained. The advantage of the new method is that the function can have a gentle profile at the location of the linear weight (or the mapped nonlinear weight can be close to its linear counterpart), and therefore is favorable for the resolution enhancement. Besides, the function also has the flexibility of quick convergence to identity mapping near two endpoints of [0,1], which is favorable for improved numerical stability. The fourth-, fifth- and sixth-order polynomial functions are constructed correspondingly with different emphasis on aforementioned flatness and convergence. Among them, the fifth-order version has the flattest profile. To check the performance of the methods, the 1-D Shu-Osher problem, the 2-D Riemann problem and the double Mach reflection are tested with the comparison of WENO-M, WENO-Z and WENO-NS. The proposed new methods show the best resolution for describing shear-layer instability of the Riemann problem, and they also indicate high resolution in computations of double Mach reflection, where only these proposed schemes successfully resolved the vortex-pairing phenomenon. Other investigations have shown that the single polynomial mapping function has no advantage over the proposed piecewise one, and it is of no evident benefit to use the proposed method for the symmetric fifth-order WENO. Overall, the fifth-order piecewise polynomial and corresponding WENO scheme are suggested for resolution improvement.     

Staggered Lagrangian Discretization Based on Cell-Centered Riemann Solver and Associated Hydrodynamics Scheme

The aim of the present work is to develop a general formalism to derive staggered discretizations for Lagrangian hydrodynamics on two-dimensional unstructured grids. To this end, we make use of the compatible discretization that has been initially introduced by E. J. Caramana et al., in J. Comput. Phys., 146 (1998). Namely, momentum equation is discretized by means of subcell forces and specific internal energy equation is obtained using total energy conservation. The main contribution of this work lies in the fact that the subcell force is derived invoking Galilean invariance and thermodynamic consistency. That is, we deduce a general form of the sub-cell force so that a cell entropy inequality is satisfied. The subcell force writes as a pressure contribution plus a tensorial viscous contribution which is proportional to the difference between the nodal velocity and the cell-centered velocity. This cell-centered velocity is a supplementary degree of freedom that is solved by means of a cell-centered approximate Riemann solver. To satisfy the second law of thermodynamics, the local subcell tensor involved in the viscous part of the subcell force must be symmetric positive definite. This subcell tensor is the cornerstone of the scheme. One particular expression of this tensor is given. A high-order extension of this discretization is provided. Numerical tests are presented in order to assess the efficiency of this approach. The results obtained for various representative configurations of one- and two-dimensional compressible fluid flows show the robustness and the accuracy of this scheme.     

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.     

Extension and Comparative Study of AUSM-Family Schemes for Compressible Multiphase Flow Simulations

Several recently developed AUSM-family numerical flux functions (SLAU, SLAU2, AUSM⁺-up2, and AUSMPW+) have been successfully extended to compute compressible multiphase flows, based on the stratified flow model concept, by following two previous works: one by M.-S. Liou, C.-H. Chang, L. Nguyen, and T.G. Theofanous [AIAA J. 46:2345-2356, 2008], in which AUSM⁺-up was used entirely, and the other by C.-H. Chang, and M.-S. Liou [J. Comput. Phys. 225:840-873, 2007], in which the exact Riemann solver was combined into AUSM⁺-up at the phase interface. Through an extensive survey by comparing flux functions, the following are found: (1) AUSM⁺-up with dissipation parameters of K_p and K_u equal to 0.5 or greater, AUSMPW+, SLAU2, AUSM⁺-up2, and SLAU can be used to solve benchmark problems, including a shock/water-droplet interaction; (2) SLAU shows oscillatory behaviors [though not as catastrophic as those of AUSM⁺ (a special case of AUSM⁺-up withK_p=K_u=0)] due to insufficient dissipation arising from its ideal-gas-based dissipation term; and (3) when combined with the exact Riemann solver, AUSM⁺-up (K_p=K_u=1), SLAU2, and AUSMPW+ are applicable to more challenging problems with high pressure ratios.     

A Front Tracking Method for the Simulation of Compressible Multimedium Flows

A front tracking method combined with the real ghost fluid method (RGFM) is proposed for simulations of fluid interfaces in two-dimensional compressible flows. In this paper the Riemann problem is constructed along the normal direction of interface and the corresponding Riemann solutions are used to track fluid interfaces. The interface boundary conditions are defined by the RGFM, and the fluid interfaces are explicitly tracked by several connected marker points. The Riemann solutions are also used directly to update the flow states on both sides of the interface in the RGFM. In order to validate the accuracy and capacity of the new method, extensive numerical tests including the bubble advection, the Sod tube, the shock-bubble interaction, the Richtmyer-Meshkov instability and the gas-water interface, are simulated by using the Euler equations. The computational results are also compared with earlier computational studies and it shows good agreements including the compressible gas-water system with large density differences.     

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

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

An All-Speed Asymptotic-Preserving Method for the Isentropic Euler and Navier-Stokes Equations

The computation of compressible flows becomes more challenging when the Mach number has different orders of magnitude. When the Mach number is of order one, modern shock capturing methods are able to capture shocks and other complex structures with high numerical resolutions. However, if the Mach number is small, the acoustic waves lead to stiffness in time and excessively large numerical viscosity, thus demanding much smaller time step and mesh size than normally needed for incompressible flow simulation. In this paper, we develop an all-speed asymptotic preserving (AP) numerical scheme for the compressible isentropic Euler and Navier-Stokes equations that is uniformly stable and accurate for all Mach numbers. Our idea is to split the system into two parts: one involves a slow, nonlinear and conservative hyperbolic system adequate for the use of modern shock capturing methods and the other a linear hyperbolic system which contains the stiff acoustic dynamics, to be solved implicitly. This implicit part is reformulated into a standard pressure Poisson projection system and thus possesses sufficient structure for efficient fast Fourier transform solution techniques. In the zero Mach number limit, the scheme automatically becomes a projection method-like incompressible solver. We present numerical results in one and two dimensions in both compressible and incompressible regimes.     

A Moving Mesh Method for the Euler Flow Calculations Using a Directional Monitor Function

This paper is concerned with the adaptive grid method for computations of the Euler equations in fluid dynamics. The new feature of the present moving mesh algorithm is the use of a dimensional-splitting type monitor function, which is to increase grid concentration in regions containing shock waves and contact discontinuities or their interactions. Several two–dimensional flow problems are computed to demonstrate the effectiveness of the present adaptive grid algorithm.     

Synchronization of Weighted Essentially Non-Oscillatory Methods

Weighted essentially non-oscillatory (WENO) methods have been developed to simultaneously provide robust shock-capturing in compressible fluid flow and avoid excessive damping of fine-scale flow features such as turbulence. Under certain conditions in compressible turbulence, however, numerical dissipation remains unacceptably high even after optimization of the linear component that dominates in smooth regions. Of the nonlinear error that remains, we demonstrate that a large fraction is generated by a "synchronization deficiency" that interferes with the expression of theoretically predicted numerical performance characteristics when the WENO adaptation mechanism is engaged. This deficiency is illustrated numerically in simulations of a linearly advected sinusoidal wave and the Shu-Osher problem [J. Comput. Phys., 83 (1989), pp. 32-78]. It is shown that attempting to correct this deficiency through forcible synchronization results in violation of conservation. We conclude that, for the given choice of candidate stencils, the synchronization deficiency cannot be adequately resolved under the current WENO smoothness measurement technique.     

On Universal Osher-Type Schemes for General Nonlinear Hyperbolic Conservation Laws

This paper is concerned with a new version of the Osher-Solomon Riemann solver and is based on a numerical integration of the path-dependent dissipation matrix. The resulting scheme is much simpler than the original one and is applicable to general hyperbolic conservation laws, while retaining the attractive features of the original solver: the method is entropy-satisfying, differentiable and complete in the sense that it attributes a different numerical viscosity to each characteristic field, in particular to the intermediate ones, since the full eigenstructure of the underlying hyperbolic system is used. To illustrate the potential of the proposed scheme we show applications to the following hyperbolic conservation laws: Euler equations of compressible gasdynamics with ideal gas and real gas equation of state, classical and relativistic MHD equations as well as the equations of nonlinear elasticity. To the knowledge of the authors, apart from the Euler equations with ideal gas, an Osher-type scheme has never been devised before for any of these complicated PDE systems. Since our new general Riemann solver can be directly used as a building block of high order finite volume and discontinuous Galerkin schemes we also show the extension to higher order of accuracy and multiple space dimensions in the new framework of P_NP_M schemes on unstructured meshes recently proposed in [9].     

Modified Upwinding Compact Scheme for Shock and Shock Boundary Layer Interaction

Standard compact scheme and upwinding compact scheme have high order accuracy and high resolution, but cannot capture the shock which is a discontinuity. This work developed a modified upwinding compact scheme which uses an effective shock detector to block compact scheme to cross the shock and a control function to mix the flux with WENO scheme near the shock. The new scheme makes the original compact scheme able to capture the shock sharply and, more importantly, keep high order accuracy and high resolution in the smooth area which is particularly important for shock boundary layer and shock acoustic interactions. Numerical results show the scheme is successful for 2-D Euler and 2-D Navier-Stokes solvers. The examples include 2-D incident shock, 2-D incident shock and boundary layer interaction. The scheme is robust, which does not involve case related parameters.     

An Optimization Method in Inverse Elastic Scattering for One-Dimensional Grating Profiles

Consider the inverse diffraction problem to determine a two-dimensional periodic structure from scattered elastic waves measured above the structure. We formulate the inverse problem as a least squares optimization problem, following the two-step algorithm by G. Bruckner and J. Elschner [Inverse Probl., 19 (2003), 315–329] for electromagnetic diffraction gratings. Such a method is based on the Kirsch-Kress optimization scheme and consists of two parts: a linear severely ill-posed problem and a nonlinear well-posed one. We apply this method to both smooth (C²) and piecewise linear gratings for the Dirichlet boundary value problem of the Navier equation. Numerical reconstructions from exact and noisy data illustrate the feasibility of the method.     

A Gas Kinetic Scheme for the Simulation of Compressible Multicomponent Flows

In this paper, a gas kinetic scheme for the compressible multicomponent flows is presented by making use of two-species BGK model in [A. D. Kotelnikov and D. C. Montgomery, A Kinetic Method for Computing Inhomogeneous Fluid Behavior, J. Comput. Phys. 134 (1997) 364-388]. Different from the conventional BGK model, the collisions between different species are taken into consideration. Based on the Chapman-Enskog expansion, the corresponding macroscopic equations are derived from this two-species model. Because of the relaxation terms in the governing equations, the method of operator splitting is applied. In the hyperbolic part, the integral solutions of the BGK equations are used to construct the numerical fluxes at the cell interface in the framework of finite volume method. Numerical tests are presented in this paper to validate the current approach for the compressible multicomponent flows. The theoretical analysis on the spurious oscillations at the interface is also presented.     

Numerical Regularized Moment Method for High Mach Number Flow

This paper is a continuation of our earlier work [SIAM J. Sci. Comput., 32(2010), pp. 2875–2907] in which a numerical moment method with arbitrary order of moments was presented. However, the computation may break down during the calculation of the structure of a shock wave with Mach number M₀≥3. In this paper, we concentrate on the regularization of the moment systems. First, we apply the Maxwell iteration to the infinite moment system and determine the magnitude of each moment with respect to the Knudsen number. After that, we obtain the approximation of high order moments and close the moment systems by dropping some high-order terms. Linearization is then performed to obtain a very simple regularization term, thus it is very convenient for numerical implementation. To validate the new regularization, the shock structures of low order systems are computed with different shock Mach numbers.     

A Novel Full-Euler Low Mach Number IMEX Splitting

In this paper, we introduce an extension of a splitting method for singularly perturbed equations, the so-called RS-IMEX splitting [Kaiser et al., Journal of Scientific Computing, 70(3), 1390–1407], to deal with the fully compressible Euler equations. The straightforward application of the splitting yields sub-equations that are, due to the occurrence of complex eigenvalues, not hyperbolic. A modification, slightly changing the convective flux, is introduced that overcomes this issue. It is shown that the splitting gives rise to a discretization that respects the low-Mach number limit of the Euler equations; numerical results using finite volume and discontinuous Galerkin schemes show the potential of the discretization.     

A Godunov-Type Solver for the Numerical Approximation of Gravitational Flows

We present a new numerical method to approximate the solutions of an Euler-Poisson model, which is inherent to astrophysical flows where gravity plays an important role. We propose a discretization of gravity which ensures adequate coupling of the Poisson and Euler equations, paying particular attention to the gravity source term involved in the latter equations. In order to approximate this source term, its discretization is introduced into the approximate Riemann solver used for the Euler equations. A relaxation scheme is involved and its robustness is established. The method has been implemented in the software HERACLES [29] and several numerical experiments involving gravitational flows for astrophysics highlight the scheme.     

Implementation of Finite Difference Weighted Compact Nonlinear Schemes with the Two-Stage Fourth-Order Accurate Temporal Discretization

In this paper, we present a new two-stage fourth-order finite difference weighted compact nonlinear scheme (WCNS) for hyperbolic conservation laws with special application to compressible Euler equations. To construct this algorithm, apart from the traditional WCNS for the spatial derivative, it was necessary to first construct a linear compact/explicit scheme utilizing time derivative of flux at midpoints, which, in turn, was solved by a generalized Riemann solver. Combining these two schemes, the fourth-order time accuracy was achieved using only the two-stage time-stepping technique. The final algorithm was numerically tested for various one-dimensional and two-dimensional cases. The results demonstrated that the proposed algorithm had an essentially similar performance as that based on the fourth-order Runge-Kutta method, while it required 25 percent less computational cost for one-dimensional cases, which is expected to decline further for multidimensional cases.     

Investigation of Riemann Solver with Internal Reconstruction (RSIR) for the Euler Equations

The Riemann solver with internal reconstruction (RSIR) of Carmouze et al. (2020) is investigated, revisited and improved for the Euler equations. In this reference, the RSIR approach has been developed to address the numerical resolution of the non-equilibrium two-phase flow model of Saurel et al. (2017). The main idea is to reconstruct two intermediate states from the knowledge of a simple and robust intercell state such as HLL, regardless the number of waves present in the Riemann problem. Such reconstruction improves significantly the accuracy of the HLL solution, preserves robustness and maintains stationary discontinuities. Consequently, when dealing with complex flow models, such as the aforementioned one, RSIR-type solvers are quite easy to derive compared to HLLC-type ones that may require a tedious analysis of the governing equations across the different waves. In the present contribution, the RSIR solver of Carmouze et al. (2020) is investigated, revisited and improved with the help of thermodynamic considerations, making a simple, accurate, robust and positive Riemann solver. It is also demonstrated that the RSIR solver is strictly equivalent to the HLLC solver of Toro et al. (1994) for the Euler equations when the Rankine-Hugoniot relations are used. In that sense, the RSIR approach recovers the HLLC solver in some limit as well as the HLL one in another limit and can be simplified at different levels when complex systems of equations are addressed. For the sake of clarity and simplicity, the derivations are performed in the context of the one-dimensional Euler equations. Comparisons and validations against the conventional HLLC solver and exact solutions are presented.     

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.