High-resolution Shock-capturing Methods Research Articles

Modeling of shock wave produced by collapse of cavitation bubble using a fully conservative multiphase model

In this study, we extensively validated the applicability of the fully conservative compressible multiphase flow model to simulate the collapse of cavitation bubbles. The mathematical model is fully conservative, which is crucial for accurately modeling the strong shock phenomena induced via bubble collapse. The numerical method accounts for compressibility and mass transfer with phase-change transitions. A high-resolution shock-capturing finite-volume method and an accurate Riemann solver were employed to capture the strong shocks produced by bubble collapse. A grid-resolution study and validation were performed to ensure the convergence and consistency of the numerical method. This study focused on analyzing the shock waves generated during the collapse as well as the dynamics of bubble collapse and rebound. Quantitative comparisons between the simulation results and experimental data as well as the reference results revealed a high level of agreement. The investigation focused on the impact of various factors, including the bubble radius, phase change, and liquid pressure, on the collapse and rebound of the cavitation bubble, as well as the resulting shock waves. These analyses yielded valuable insight into the complex behavior of cavitation bubbles that collapse in a free field.

No Evidence of Kinetic Screening in Simulations of Merging Binary Neutron Stars beyond General Relativity.

We have conducted fully relativistic simulations in a class of scalar-tensor theories with derivative self-interactions and screening of local scales. By using high-resolution shock-capturing methods and a nonvanishing shift vector, we have managed to avoid issues plaguing similar attempts in the past. We have first confirmed recent results by ourselves in spherical symmetry, obtained with an approximate approach and pointing at a partial breakdown of the screening in black-hole collapse. Then, we considered the late inspiral and merger of binary neutron stars. We found that screening tends to suppress the (subdominant) dipole scalar emission, but not the (dominant) quadrupole scalar mode. Our results point at quadrupole scalar signals as large as (or even larger than) in Fierz-Jordan-Brans-Dicke theories with the same conformal coupling, for strong-coupling scales in the MeV range that we can simulate.

Fully compressible multiphase model for computation of compressible fluid flows with large density ratio and the presence of shock waves

We propose a fully compressible multiphase model for simulating compressible interfacial flows. The mathematical model is formulated on the basis of the typical conservation laws for mixtures, including mass, momentum, and energy balance. The mass balance equation in each phase is considered to discretely conserve the mass of each phase; therefore, it is always mass conservative for each phase. The numerical solution procedure for this model is developed using a high-resolution shock-capturing finite-volume method based on the Harten–Lax–van Leer-contact (HLLC) Riemann solver and total variation diminishing (TVD) time-integration schemes for simulating unsteady multiphase flows and capturing strong shocks. The finite-volume Riemann solver is combined with a monotonic upstream-centered scheme for conservation laws (MUSCL) and compressive limiters to maintain a sharp interface. The numerical framework yields a conceptually simple but robust and versatile numerical procedure. Numerical results of benchmark Riemann problems demonstrate the ability of the proposed method to obtain sharp interfaces, free of spurious oscillations, high accuracy, conservation properties, and thermodynamic consistency over various Mach numbers.

Dynamical evolution of dark matter admixed neutron stars

We dynamically evolve for the first time dark matter admixed neutron stars with fermionic dark matter. These systems are mixtures of the ordinary nuclear matter of a neutron star and dark matter. To perform our dynamical evolutions, we derive the equations of motion, in conservation form, for spherically symmetric systems with an arbitrary number of perfect fluids. Using finite volume and high-resolution shock-capturing methods, we dynamically evolve the two-fluid case, with the first fluid modeling ordinary matter and the second fluid modeling dark matter. We use our dynamical solutions to study nonlinear stability, radial oscillation frequencies, and a dynamical formation process.

Gmunu: paralleled, grid-adaptive, general-relativistic magnetohydrodynamics in curvilinear geometries in dynamical space–times

ABSTRACT We present an update on the General-relativistic multigrid numerical (Gmunu) code, a parallelized, multidimensional curvilinear, general relativistic magnetohydrodynamics code with an efficient non-linear cell-centred multigrid elliptic solver, which is fully coupled with an efficient block-based adaptive mesh refinement module. To date, as described in this paper, Gmunu is able to solve the elliptic metric equations in the conformally flat condition approximation with the multigrid approach and the equations of ideal general-relativistic magnetohydrodynamics by means of high-resolution shock-capturing finite-volume method with reference metric formularised multidimensionally in Cartesian, cylindrical, or spherical geometries. To guarantee the absence of magnetic monopoles during the evolution, we have developed an elliptical divergence cleaning method by using the multigrid solver. In this paper, we present the methodology, full evolution equations and implementation details of Gmunu and its properties and performance in some benchmarking and challenging relativistic magnetohydrodynamics problems.

Numerical exploration of first-order relativistic hydrodynamics

We present the first numerical solutions of the causal, stable relativistic Navier-Stokes equations as formulated by Bemfica, Disconzi, Noronha, and Kovtun (BDNK). For this initial investigation we restrict to plane-symmetric configurations of a conformal fluid in Minkowski spacetime. We consider evolution of three classes of initial data: a smooth (initially) stationary concentration of energy, a standard shock tube setup, and a smooth shockwave setup. We compare these solutions to those obtained with the Muller-Israel-Stewart (MIS) formalism, variants of which are the common tools used to model relativistic, viscous fluids. We find that for the two smooth initial data cases, simple finite difference methods are adequate to obtain stable, convergent solutions to the BDNK equations. For low viscosity, the MIS and BDNK evolutions show good agreement. At high viscosity the solutions begin to differ in regions with large gradients, and there the BDNK solutions can (as expected) exhibit violation of the weak energy condition. This behavior is transient, and the solutions evolve toward a hydrodynamic regime in a way reminiscent of an approach to a universal attractor. For the shockwave problem, we give evidence that if a hydrodynamic frame is chosen so that the maximum characteristic speed of the BDNK system is the speed of light (or larger), arbitrarily strong shockwaves are smoothly resolved. Regarding the shock tube problem, it is unclear whether discontinuous initial data is mathematically well-posed for the BDNK system, even in a weak sense. Nevertheless we attempt numerical solution, and then need to treat the perfect fluid terms using high-resolution shock-capturing methods. When such methods can successfully evolve the solution beyond the initial time, subsequent evolution agrees with corresponding MIS solutions, as well as the perfect fluid solution in the limit of zero viscosity.

High-resolution shock-capturing numerical simulations of three-phase immiscible fluids from the unsaturated to the saturated zone

Numerical modeling of immiscible contaminant fluid flow in unsaturated and saturated porous aquifers is of great importance in many scientific fields to properly manage groundwater resources. We present a high-resolution numerical model that simulates three-phase immiscible fluid flow in both unsaturated and saturated zone in a porous aquifer. We use coupled conserved mass equations for each phase and study the dynamics of a multiphase fluid flow as a function of saturation, capillary pressure, permeability, and porosity of the different phases, initial and boundary conditions. To deal with the sharp front originated from the partial differential equations’ nonlinearity and accurately propagate the sharp front of the fluid component, we use a high-resolution shock-capturing method to treat discontinuities due to capillary pressure and permeabilities that depend on the saturation of the three different phases. The main approach to the problem’s numerical solution is based on (full) explicit evolution of the discretized (in-space) variables. Since explicit methods require the time step to be sufficiently small, this condition is very restrictive, particularly for long-time integrations. With the increased computational speed and capacity of today’s multicore computer, it is possible to simulate in detail contaminants’ fate flow using high-performance computing.

A simple extended compact nonlinear scheme with adaptive dissipation control

Weighted compact nonlinear schemes (WCNSs) [1] are a popular family of high resolution shock-capturing methods, of which the nonlinear interpolation procedure is essential for non-oscillatory shock-capturing computation. In this article, numerical dissipation of an improved compact nonlinear scheme, which is fifth-order accurate, is significantly reduced by further introducing a six-point optimized stencil in the nonlinear interpolation, without using more complex local or global smoothness indicators. Nonlinear dissipation and dispersion are then further reduced by using an adaptive strategy. In particular, the computational procedure including the proposed six-point nonlinear interpolation is almost the same as that of the improved five-point scheme, and thus their computational efforts are similar. However, significant gains in numerical performance are obtained while including the six-point optimized stencil. Benchmark cases in one and two space dimensions are used to assess performance of the present scheme, viz., shock-capturing capability and resolution in smooth regions, showing promising advantages.

CAFE-Q: Code Designed to Solve the Resistive MHD Equations With Thermal Conductivity

We present a new code designed to solve the equations of classical magnetohydrodynamics in 3-D under the effects of magnetic resistivity and heat transfer. The purpose of the code centers on the analysis of solar phenomena within the photosphere–corona region and the study of the propagation of coronal mass ejections in the interplanetary medium. We present 1-D and 2-D standard tests to demonstrate the quality of the numerical results obtained with our code. As a new code test, we include a system directly related to solar physics: the formation of a jet triggered by magnetic reconnection. The program uses high-resolution shock-capturing methods, using HLLE, HLLC, and HLLD flux formulas combined with MINMOD, MC, and reconstructors. The divergence-free magnetic field constraint is controlled using the extended generalized Lagrange multipliers method and flux-constrained transport (Flux-CT) method. Concerning the accuracy of the numerical results, we present three important cases: first, for the self-similar current sheet test, the error with respect to the exact solution is 0.13%; second, about the violation of the divergence-free magnetic field constraint, we use the blast wave test and show the constraint is of order 10−11; and third, related to the magnetic field reconnection, the code reproduces the predicted reconnection rate for the Sweet–Parker model.

Numerical Methods for Interface Coupling of Compressible and Almost Incompressible Media

Many experiments in biomedical applications and other disciplines use a shock tube. These experiments often involve placing an experimental sample within a fluid-filled container, which is then placed inside the shock tube. The shock tube produces an initial shock that propagates through gas before hitting the container with the sample. In order to gain insight into the shock dynamics that is hard to obtain by experimental means, computational simulations of the shock wave passing from gas into a thin elastic solid and into a nearly incompressible fluid are developed. It is shown that if the solid interface is very thin, it can be neglected, simplifying the model. The model uses Euler equations for compressible fluids coupled with a Tammann equation of state (EOS) to model both compressible gas and almost incompressible materials. A three-dimensional (2D axisymmetric) model of these equations is solved using high-resolution shock-capturing methods, with newly developed Riemann solvers and limiters. The methods are extended to work on a mapped grid to allow more complicated interface geometry, and they are adapted to work with adaptive mesh refinement (AMR) for higher resolution and faster computations. The Clawpack software is used to implement the method. These methods were initially inspired by shock tube experiments to study the injury mechanisms of traumatic brain injury (TBI).

Grid-based Methods in Relativistic Hydrodynamics and Magnetohydrodynamics

An overview of grid-based numerical methods used in relativistic hydrodynamics (RHD) and magnetohydrodynamics (RMHD) is presented. Special emphasis is put on a comprehensive review of the application of high-resolution shock-capturing methods. Results of a set of demanding test bench simulations obtained with different numerical methods are compared in an attempt to assess the present capabilities and limits of the various numerical strategies. Applications to three astrophysical phenomena are briefly discussed to motivate the need for and to demonstrate the success of RHD and RMHD simulations in their understanding. The review further provides FORTRAN programs to compute the exact solution of the Riemann problem in RMHD, and to simulate 1D RMHD flows in Cartesian coordinates.Electronic Supplementary MaterialSupplementary material is available for this article at 10.1007/lrca-2015-3.

Newtonian cafe: a new ideal MHD code to study the solar atmosphere

We present a new code designed to solve the equations of classical ideal magneto-hydrodynamics (MHD) in three dimensions, submitted to a constant gravitational field. The purpose of the code centers on the analysis of solar phenomena within the photosphere-corona region. We present 1D and 2D standard tests to demonstrate the quality of the numerical results obtained with our code. As solar tests we present the transverse oscillations of Alfvenic pulses in coronal loops using a 2.5D model, and as 3D tests we present the propagation of impulsively generated MHD-gravity waves and vortices in the solar atmosphere. The code is based on high-resolution shock-capturing methods, uses the HLLE flux formula combined with Minmod, MC and WENO5 reconstructors. The divergence free magnetic field constraint is controlled using the Flux Constrained Transport method.

CAFE: A NEW RELATIVISTIC MHD CODE

We introduce CAFE, a new independent code designed to solve the equations of Relativistic ideal Magnetohydrodynamics (RMHD) in 3D. We present the standard tests for a RMHD code and for the Relativistic Hydrodynamics (RHD) regime since we have not reported them before. The tests include the 1D Riemann problems related to blast waves, head-on collision of streams and states with transverse velocities, with and without magnetic field, which is aligned or transverse, constant or discontinuous across the initial discontinuity. Among the 2D and 3D tests, without magnetic field we include the 2D Riemann problem, a one dimensional shock tube along a diagonal, the high speed Emery wind tunnel, the Kelvin-Helmholtz instability, a set of jets and a 3D spherical blast wave, whereas in the presence of a magnetic field we show the magnetic rotor, the cylindrical explosion, a case of Kelvin-Helmholtz instability and a 3D magnetic field advection loop. The code uses High Resolution Shock Capturing methods and we present the error analysis for a combination that uses the HLLE flux formula combined with linear, PPM and fifth order WENO reconstructors. We use the flux-CT and the divergence cleaning methods to control the divergence free magnetic field constraint.

Numerical solutions of ideal quantum gas dynamical flows governed by semiclassical ellipsoidal-statistical distribution

The ideal quantum gas dynamics as manifested by the semiclassical ellipsoidal-statistical (ES) equilibrium distribution derived in Wu et al. (Wu et al. 2012 Proc. R. Soc. A468, 1799-1823 (doi:10.1098/rspa.2011.0673)) is numerically studied for particles of three statistics. This anisotropic ES equilibrium distribution was derived using the maximum entropy principle and conserves the mass, momentum and energy, but differs from the standard Fermi-Dirac or Bose-Einstein distribution. The present numerical method combines the discrete velocity (or momentum) ordinate method in momentum space and the high-resolution shock-capturing method in physical space. A decoding procedure to obtain the necessary parameters for determining the ES distribution is also devised. Computations of two-dimensional Riemann problems are presented, and various contours of the quantities unique to this ES model are illustrated. The main flow features, such as shock waves, expansion waves and slip lines and their complex nonlinear interactions, are depicted and found to be consistent with existing calculations for a classical gas.

Hydrodynamics in full general relativity with conservative adaptive mesh refinement

There is great interest in numerical relativity simulations involving matter due to the likelihood that binary compact objects involving neutron stars will be detected by gravitational wave observatories in the coming years, as well as to the possibility that binary compact object mergers could explain short-duration gamma-ray bursts. We present a code designed for simulations of hydrodynamics coupled to the Einstein field equations targeted toward such applications. This code has recently been used to study eccentric mergers of black hole-neutron star binaries. We evolve the fluid conservatively using high-resolution shock-capturing methods, while the field equations are solved in the generalized-harmonic formulation with finite differences. In order to resolve the various scales that may arise, we use adaptive mesh refinement (AMR) with grid hierarchies based on truncation error estimates. A noteworthy feature of this code is the implementation of the flux correction algorithm of Berger and Colella to ensure that the conservative nature of fluid advection is respected across AMR boundaries. We present various tests to compare the performance of different limiters and flux calculation methods, as well as to demonstrate the utility of AMR flux corrections.

A NEW MULTI-DIMENSIONAL GENERAL RELATIVISTIC NEUTRINO HYDRODYNAMIC CODE FOR CORE-COLLAPSE SUPERNOVAE. I. METHOD AND CODE TESTS IN SPHERICAL SYMMETRY

We present a new general relativistic (GR) code for hydrodynamic supernova simulations with neutrino transport in spherical and azimuthal symmetry (1D/2D). The code is a combination of the CoCoNuT hydro module, which is a Riemann-solver based, high-resolution shock-capturing method, and the three-flavor, energy-dependent neutrino transport scheme VERTEX. VERTEX integrates the neutrino moment equations with a variable Eddington factor closure computed from a model Boltzmann equation and uses the ray-by-ray plus approximation in 2D, assuming the neutrino distribution to be axially symmetric around the radial direction, and thus the neutrino flux to be radial. Our spacetime treatment employs the ADM 3+1 formalism with the conformal flatness condition for the spatial three-metric. This approach is exact in 1D and has been shown to yield very accurate results also for rotational stellar collapse. We introduce new formulations of the energy equation to improve total energy conservation in relativistic and Newtonian hydro simulations with Eulerian finite-volume codes. Moreover, a modified version of the VERTEX scheme is developed that simultaneously conserves energy and lepton number with better accuracy and higher numerical stability. To verify our code, we conduct a series of tests, including a detailed comparison with published 1D results for stellar core collapse. Long-time simulations of proto-neutron star cooling over several seconds both demonstrate the robustness of the new CoCoNuT-VERTEX code and show the approximate treatment of GR effects by means of an effective gravitational potential as in PROMETHEUS-VERTEX to be remarkably accurate in 1D. (abridged)

Simulation of turbulent flows containing strong shocks

Simulation of turbulent flows with strong shocks is a computationally challenging problem. The requirements for a method to produce accurate results for turbulence are orthogonal to those needed to treat shocks properly. In order to prevent an unphysical rate of decay of turbulent structures, it is necessary to use a method with very low numerical dissipation. Because of this, central difference schemes are widely used. However, computing strong shocks with a central difference scheme can produce unphysical post-shock oscillations that corrupt the entire flow unless additional dissipation is added. This dissipation can be difficult to localize to the area near the shock and can lead to inaccurate treatment of the turbulence. Modern high-resolution shock-capturing methods usually use upwind algorithms to provide the dissipation necessary to stabilize shocks. However, this upwind dissipation can also lead to an unphysical rate of decay of the turbulence. This paper discusses a hybrid method for simulating turbulent flows with strong shocks that couples a high-order central difference scheme with a high-resolution shock-capturing method. The shock-capturing method is used only in the vicinity of discontinuities in the flow, whereas the central difference scheme is used in the remainder of the computational domain. Results of this new method will be shown for a variety of test problems. Preliminary results for a realistic application involving detonation in gas-particle flows will also be presented.

An efficient adaptive mesh redistribution method for a non-linear Dirac equation

This paper presents an efficient adaptive mesh redistribution method to solve a non-linear Dirac (NLD) equation. Our algorithm is formed by three parts: the NLD evolution, the iterative mesh redistribution of the coarse mesh and the local uniform refinement of the final coarse mesh. At each time level, the equidistribution principle is first employed to iteratively redistribute coarse mesh points, and the scalar monitor function is subsequently interpolated on the coarse mesh in order to do one new iteration and improve the grid adaptivity. After an adaptive coarse mesh is generated ideally and finally, each coarse mesh interval is equally divided into some fine cells to give an adaptive fine mesh of the physical domain, and then the solution vector is remapped on the resulting new fine mesh by an affine method. The NLD equation is finally solved by using a high resolution shock-capturing method on the (fixed) non-uniform fine mesh. Extensive numerical experiments demonstrate that the proposed adaptive mesh method gives the third-order rate of convergence, and yields an efficient and fast NLD solver that tracks and resolves both small, local and large solution gradients automatically.

Non-Linear Filtering and Limiting in High Order Methods for Ideal and Non-Ideal MHD

The adaptive nonlinear filtering and limiting in spatially high order schemes (Yee et al. J. Comput. Phys. 150, 199–238, (1999), Sjögreen and Yee, J. Scient. Comput. 20, 211–255, (2004)) for the compressible Euler and Navier–Stokes equations have been recently extended to the ideal and non-ideal magnetohydrodynamics (MHD) equations, (Sjögreen and Yee, (2003), Proceedings of the 16th AIAA/CFD conference, June 23–26, Orlando F1; Yee and Sjögreen (2003), Proceedings of the International Conference on High Performance Scientific Computing, March, 10–14, Honai, Vietnam; Yee and Sjögreen (2003), RIACS Technical Report TR03. 10, July, NASA Ames Research Center; Yee and Sjögreen (2004), Proceedings of the ICCF03, July 12–16, Toronto, Canada). The numerical dissipation control in these adaptive filter schemes consists of automatic detection of different flow features as distinct sensors to signal the appropriate type and amount of numerical dissipation/filter where needed and leave the rest of the region free from numerical dissipation contamination. The numerical dissipation considered consists of high order linear dissipation for the suppression of high frequency oscillation and the nonlinear dissipative portion of high-resolution shock-capturing methods for discontinuity capturing. The applicable nonlinear dissipative portion of high-resolution shock-capturing methods is very general. The objective of this paper is to investigate the performance of three commonly used types of discontinuity capturing nonlinear numerical dissipation for both the ideal and non-ideal MHD.

Numerical Hydrodynamics in Special Relativity.

This review is concerned with a discussion of numerical methods for the solution of the equations of special relativistic hydrodynamics (SRHD). Particular emphasis is put on a comprehensive review of the application of high-resolution shock-capturing methods in SRHD. Results of a set of demanding test bench simulations obtained with different numerical SRHD methods are compared. Three applications (astrophysical jets, gamma-ray bursts and heavy ion collisions) of relativistic flows are discussed. An evaluation of various SRHD methods is presented, and future developments in SRHD are analyzed involving extension to general relativistic hydrodynamics and relativistic magneto-hydrodynamics. The review further provides FORTRAN programs to compute the exact solution of a 1D relativistic Riemann problem with zero and nonzero tangential velocities, and to simulate 1D relativistic flows in Cartesian Eulerian coordinates using the exact SRHD Riemann solver and PPM reconstruction.Electronic Supplementary MaterialSupplementary material is available for this article at 10.12942/lrr-2003-7 and is accessible for authorized users.