Numerical Diffusion Research Articles

An efficient Ghost Fluid Method to remove overheating from material interfaces in compressible multi-medium flows

The Ghost Fluid Method (GFM) and its derivatives represent powerful techniques in the numerical simulation of compressible, multi-medium flows, but are nonetheless afflicted by so-called overheating errors at and near material interfaces. These errors take the form of local overshoots and undershoots, and cause the numerical solution to deviate from the exact solution. Overheating errors are commonly addressed by an isobaric fix (or an isentropic fix), where the values at the cells next to the fluid interface are corrected to reduce the numerical diffusion stemming from overheating. These approaches however, do not completely remove overheating, and numerical inaccuracies persist near the interface. In this article, we propose a new version of the GFM called the Efficient Ghost Fluid Method (EGFM) capable of completely eliminating overheating errors, leading to highly accurate solutions near the interface. The EGFM approach has been implemented in IMPACT, a multi-medium, shock physics code, and has been validated against a wide range of 1D and 2D test cases, including problems with surface tension. These include single- as well as multi-medium shock tube problems, various shock-interface interactions, and shock-driven instabilities.

An improved moving particle semi-implicit method for interfacial flows

In this paper, the improved moving particle semi-implicit (IMPS) method is further developed into a multiphase method by introducing various multiphase models. Then, the developed multiphase MPS method is applied to a variety of 2-D/3-D simulations of interfacial flows, including Rayleigh-Taylor instability, bubble rising, dam-break flow, and internal solitary waves. The fair agreements between the results of MPS and other reference results demonstrate that the present multiphase MPS method is stable and accurate enough to capture the interface with a large deformation in different flow problems. The phenomenon of air entrapment induced by the dam-break flow is well reproduced in the multiphase simulation and its influence on impact pressure is discussed in detail. The internal solitary waves with different amplitudes are numerically generated and the wave heights show a good agreement with experimental data even after a long-distance propagation, indicating that the numerical diffusion of the present multiphase MPS method is relatively low and the energy conservation can be ensured.

Turbulent kinetic dissipation analysis for residual-based large eddy simulation of incompressible turbulent flow by variational multiscale method

The underlying physical mechanism of the residual-based large eddy simulation (LES) based on the variational multiscale (VMS) method is clarified. Resolved large-scale energy transportation equation is mathematically derived for turbulent kinetic energy budget analysis. Firstly, statistical results of benchmark turbulent channel flow at Reτ=180 obtained using a coarse mesh are compared with the results obtained by the classical LES with the Smagorinsky and dynamic subgrid stress (SGS) model. The present LES shows an advantage in predicting the statistical results of the incompressible turbulent flows. Secondly, the contributions of the unresolved small-scale presentation terms (Term I-IV in Eq. (10)) to the turbulent kinetic dissipation are analysed for the VMS method. The results show that the turbulent kinetic dissipation provided by the numerical diffusion in the VMS method is smaller in the inner layer, larger in the outer layer of the channel flow than those by the Smagorinsky and dynamic SGS model. The turbulent kinetic dissipation in the VMS method is mainly given by the numerical diffusion provided by one of the “cross-stress” terms (Term I, same as the stabilization term in the SUPG method) and LSIC term (Term IV). The other one of the “cross-stress” terms (Term II) gives rise to the positive turbulent kinetic energy budget, and does not dissipate the turbulent kinetic energy. The so-called “Reynolds stress” term (Term III) dissipates the turbulent energy but provides a very small numerical diffusion. Finally, on the basis of the turbulent kinetic energy dissipation analysis, a new residual-based stabilized finite element formulation is proposed by modifying the large-scale equation in the VMS method. Numerical experiments of 2D lid-driven cavity flow and 3D incompressible turbulent channel flow are tested to validate the proposed formulation. It is shown that all the stabilization terms in the proposed formulation produce additional numerical diffusions and physically increase the total turbulent kinetic dissipation. Consequently, an apparent improvement in both the first-order and second-order statistical quantities are pursued by the new stabilized finite element formulation.

A compositional Eulerian approach for modeling oil spills in the sea

We present a Compositional Eulerian model to forecast the evolution of oil spills in the sea. The model allows studying the fate of not only the oil concentration but also of each component (e.g., volatile, non-volatile, water in the oil). Therefore, the problem is formulated as a conservation equation for each component, plus an equation to estimate the age of the oil, which allows us to assess weathering processes (e.g., evaporation, natural dispersion, emulsion) and the associated changes in oil properties. We describe an efficient implementation, using second order numerical schemes for advection and nonlinear diffusion terms, to reduce numerical diffusion. We perform numerical experiments, based on real and synthetic cases, to illustrate and validate the capabilities of our model to forecast the evolution of oil spills and to perform environmental risk analysis in the case of a potential accident.

Performance of fast fluid dynamics with a semi-Lagrangian scheme and an implicit upwind scheme in simulating indoor/outdoor airflow

Computational fluid dynamics can be time consuming for predicting indoor airflows and pollutant transport in large-scale problems or emergency management. Fast fluid dynamics (FFD) is able to accomplish efficient and accurate simulation of indoor/outdoor airflow. FFD solves the advection term of the Navier–Stokes equations either by a semi-Lagrangian (SL) scheme or an implicit upwind (IU) scheme. The SL scheme can be highly efficient, but its first-order version is not conservative and introduces significant numerical diffusion. To improve its accuracy, a high-order temporal and interpolation scheme that not only reduces dissipation and dispersion errors but also guarantees the convergence speed should be applied. Otherwise, an IU scheme instead could be used to solve the advection term. The IU scheme is conservative and introduces minor numerical diffusion, but it may increase the computation time. Therefore, this study investigated the performance of FFD with SL scheme using high-order temporal and interpolation schemes and that with IU scheme. The comparisons used experimental data of two indoor airflows and one outdoor airflow. The results showed that FFD with IU scheme was overall more accurate than FFD with SL scheme. In simulating indoor airflow, both methods were robust and the predictions were independent of time step sizes if the Courant number was less than or equal to one. In simulating the outdoor airflow, the FFD with SL scheme performed better than the FFD with IU scheme for large time step sizes. The FFD with IU scheme consumed 44%–61% computing time of the FFD with SL scheme.

Boundary variation diminished conservative semi-Lagrangian method for both compressible and incompressible flows

Motivated by the enlightenment that diminishing the jump at the cell boundary can effectively reduce numerical dissipation near the critical region, and a novel constrained interpolation profile conservative semi-Lagrangian method is proposed based on a newly designed boundary variation diminishing algorithm. First, a constrained interpolation profile conservative semi-Lagrangian scheme with the piecewise tangent of hyperbola for interface capturing scheme is proposed as one candidate to represent jump-like discontinuities. Second, the constrained interpolation profile conservative semi-Lagrangian scheme with a fourth-order weighted essentially non-oscillatory limiter is used as another candidate to keep the high-order and non-oscillatory reconstruction for smooth solutions. The selection criterion of these two candidates is designed by minimizing the total variations of the first derivative at cell boundaries. A unified pressure-based projection formulation with a fractional step procedure is implemented with the proposed scheme to simulate both compressible and incompressible flows. A variety of numerical tests are studied, including linear and nonlinear scalar wave transport problems and compressible and incompressible flow problems. Results show that the proposed method can effectively eliminate numerical oscillation and diffusion, suggesting it has great potential to be applied to various types of engineering problems including both compressible and incompressible flows.

Determining the Global Scale Size of Chorus Waves in the Magnetosphere

AbstractChorus waves outside the plasmapause influence the Earth's radiation belt dynamics by interacting with energetic electrons via cyclotron and Landau resonance. Recent numerical diffusion experiments indicate that the diffusion process is sensitive to the spatial and temporal scale of variability in the wave‐particle interaction, which is reported to be more efficient than that based on the traditional average model. Using Van Allen Probes A and B data from November 2012 to July 2019, the spatial and temporal scale size of chorus waves are calculated by the correlation between the wave amplitudes detected by two satellites with varying spatial separation or time lag. We found that, the chorus wave is incoherent when the spatial extent is greater than 433 km or the time lag lasts ∼10 s, which are significantly smaller than that of plasmaspheric hiss. In addition, the spatial correlations of chorus tend to be higher near noon or with lower geomagnetic activity. The temporal correlations of chorus are always statistically near zero, which are not influenced by the location and geomagnetic activity. Our results can help refine the model of the interactions between energetic particles and chorus waves in the radiation belt.

Asteroid magnetization from the early solar wind

ABSTRACT Magnetic fields provide an important probe of the thermal, material, and structural history of planetary and sub-planetary bodies. Core dynamos are a potential source of magnetic fields for differentiated bodies, but evidence of magnetization in undifferentiated bodies requires a different mechanism. Here, we study the amplified field provided by the stellar wind to an initially unmagnetized body using analytic theory and numerical simulations, employing the resistive magnetohydrodynamic AstroBEAR adaptive mesh refinement multiphysics code. We obtain a broadly applicable scaling relation for the peak magnetization achieved once a wind advects, piles-up, and drapes a body with magnetic field, reaching a quasi-steady state. We find that the dayside magnetic field for a sufficiently conductive body saturates when it balances the sum of incoming solar wind ram, magnetic, and thermal pressures. Stronger amplification results from pile-up by denser and faster winds. Careful quantification of numerical diffusivity is required for accurately interpreting the peak magnetic field strength from simulations, and corroborating with theory. As specifically applied to the Solar system, we find that early solar wind-induced field amplification is a viable source of magnetization for observed paleointensities in meteorites from some undifferentiated bodies. This mechanism may also be applicable to other Solar system bodies, including metal-rich bodies to be visited in future space missions such as the asteroid (16) Psyche.

Modelling Settling-Driven Gravitational Instabilities at the Base of Volcanic Clouds Using the Lattice Boltzmann Method

Field observations and laboratory experiments have shown that ash sedimentation can be significantly affected by collective settling mechanisms that promote premature ash deposition, with important implications for dispersal and associated impacts. Among these mechanisms, settling-driven gravitational instabilities result from the formation of a gravitationally-unstable particle boundary layer (PBL) that grows between volcanic ash clouds and the underlying atmosphere. The PBL destabilises once it reaches a critical thickness characterised by a dimensionless Grashof number, triggering the formation of rapid, downward-moving ash fingers that remain poorly characterised. We simulate this process by coupling a Lattice Boltzmann model, which solves the Navier-Stokes equations for the fluid phase, with a Weighted Essentially Non Oscillatory (WENO) finite difference scheme which solves the advection-diffusion-settling equation describing particle transport. Since the physical problem is advection dominated, the use of the WENO scheme reduces numerical diffusivity and ensures accurate tracking of the temporal evolution of the interface between the layers. We have validated the new model by showing that the simulated early-time growth rate of the instability is in very good agreement with that predicted by linear stability analysis, whilst the modelled late-stage behaviour also successfully reproduces quantitative results from published laboratory experiments. The results show that the model is capable of reproducing both the growth of the unstable PBL and the non-linear dependence of the fingers’ vertical velocity on both the initial particle concentration and the particle diameter. Our validated model is used to expand the parameter space explored experimentally and provides key insights into field studies. Our simulations reveal that the critical Grashof number for the instability is about ten times larger than expected by analogy with thermal convection. Moreover, as in the experiments, we found that instabilities do not develop above a given particle threshold. Finally, we quantify the evolution of the mass of particles deposited at the base of the numerical domain and demonstrate that the accumulation rate increases with time, while it is expected to be constant if particles settle individually. This suggests that real-time measurements of sedimentation rate from volcanic clouds may be able to distinguish finger sedimentation from individual particle settling.

Mixed Mesh Finite Volume Method for 1D Hyperbolic Systems with Application to Plug-Flow Heat Exchangers

We present a finite volume method formulated on a mixed Eulerian-Lagrangian mesh for highly advective 1D hyperbolic systems altogether with its application to plug-flow heat exchanger modeling/simulation. Advection of sharp moving fronts is an important problem in fluid dynamics, and even a simple transport equation cannot be solved precisely by having a finite number of nodes/elements/volumes. Finite volume methods are known to introduce numerical diffusion, and there exist a wide variety of schemes to minimize its occurrence; the most recent being adaptive grid methods such as moving mesh methods or adaptive mesh refinement methods. We present a solution method for a class of hyperbolic systems with one nonzero time-dependent characteristic velocity. This property allows us to rigorously define a finite volume method on a grid that is continuously moving by the characteristic velocity (Lagrangian grid) along a static Eulerian grid. The advective flux of the flowing field is, by this approach, removed from cell-to-cell interactions, and the ability to advect sharp fronts is therefore enhanced. The price to pay is a fixed velocity-dependent time sampling and a time delay in the solution. For these reasons, the method is best suited for systems with a dominating advection component. We illustrate the method’s properties on an illustrative advection-decay equation example and a 1D plug flow heat exchanger. Such heat exchanger model can then serve as a convection-accurate dynamic model in estimation and control algorithms for which it was developed.

A vertex scheme for two-phase flow in heterogeneous media

This paper presents the numerical solution of immiscible two-phase flows in porous media, obtained by a first-order finite element method equipped with mass-lumping and flux upwinding. The unknowns are the physical phase pressure and phase saturation. Our numerical experiments confirm that the method converges optimally for manufactured solutions. For both structured and unstructured meshes, we observe the high-accuracy wetting saturation profile that ensures minimal numerical diffusion at the front. Performing several examples of quarter-five spot problems in two and three dimensions, we show that the method can easily handle heterogeneities in the permeability field. Two distinct features that make the method appealing to reservoir simulators are: (i) maximum principle is satisfied, and (ii) local mass error is small and remains below nonlinear solver tolerance.

The time-dependent Rayleigh–Taylor instability in interstellar shells and supershells, including the eROSITA bubbles

ABSTRACT The Rayleigh–Taylor (RT) instability is omnipresent in the physics of inversely density-stratified fluids subject to effective gravitational acceleration. In astrophysics, a steep stratification of the ambient medium can fragment a bubble shell faster due to a strongly time-dependent RT instability, causing the classical constant gravity models to fail. We derive the time-dependent instability criteria analytically for the cases of constant, exponential, and power-law accelerations, verifying them through high-resolution numerical simulations. Our results show that (1) even in the linear phase there is a term opposing exponential growth, (2) non-linear growth approaches asymptotically the solution found by Fermi and von Neumann, (3) the interpenetrating spikes and bubbles promote a significant mixing, with the fractal dimension of the interface approaching 1.6, only limited by numerical diffusion, and (4) the probability density function for the passive scalar to study mixing becomes increasingly sharper peaked for power-law and exponential accelerations. Applying our solutions to stellar wind bubbles, young supernova remnants (SNRs), and superbubbles (SBs), we find that the growth rate of the RT instability is generally higher in the shells of wind-blown bubbles in a power-law stratified medium than in those with power-law rising stellar mechanical luminosities, Tycho-like than Cas A-like SNRs, and one-sided than symmetric SBs. The recently observed eROSITA bubbles indicate smooth rim surfaces, implying that the outer shell has not been affected by RT instabilities. Therefore, the dynamical evolution of the bubbles suggests maximum final ages that are significantly above their current age, which we estimate to be about 20 Myr.

Development of numerical methods to simulate the melting of a thermal protection system

In this paper, numerical methods are developed and detailed in order to be able to compute the ablation of metallic thermal protection system. In this complex multi-physics problem, the thermal state inside a solid domain and a two-phase viscous flow have to be computed. Since the two fluid phases are non-miscible, an extension of the five-equation model to dissipative effects is considered. An operator splitting strategy is used to separate the different phenomena according to their own propagation speed. An implicit time integration is performed for the acoustic+dissipative step while the transport step is computed with an explicit scheme. The hyperbolic part of the acoustic+dissipative step is solved in a non-conservative form using a Godunov-type scheme based on a simple Riemann solver. A classical discretization is used for the dissipative terms, and also for the heat equation inside the solid domain. Several approaches are detailed in order to prevent the numerical diffusion of the material interface in the transport step. Finally, since moving grids are used to capture accurately the melting front, an ALE formulation of the numerical schemes for both fluid and solid domains is given in a multidimensional framework. A fluid-solid coupling algorithm is then proposed to compute such a complex multi-physics problem.Numerical simulations show the validity and the robustness of the implicit-explicit scheme used for the discretization of the five-equation system. The last test case, namely the melting of an aluminium solid block by a lid-driven cavity filled with air, shows that the numerical tools developed here are robust enough to compute complex configurations involving a two-phase flow with high density ratios and a solid part.

Simple and robust contact-discontinuity capturing central algorithms for high speed compressible flows

The nonlinear convection terms in the governing equations of compressible fluid flows are hyperbolic in nature and are nontrivial for modelling and numerical simulation. Many numerical methods have been developed in the last few decades for this purpose and are typically based on Riemann solvers, which are strongly dependent on the underlying eigen-structure of the governing equations. Objective of the present work is to develop simple algorithms which are not dependent on the eigen-structure and yet can tackle easily the hyperbolic parts. Central schemes with smart diffusion mechanisms are apt for this purpose. For fixing the numerical diffusion, the basic ideas of satisfying the Rankine-Hugoniot (RH) conditions along with generalized Riemann invariants are proposed. Two such interesting algorithms are presented, which capture grid-aligned steady contact discontinuities exactly and yet have sufficient numerical diffusion to avoid numerical shock instabilities. Both the algorithms presented are robust in avoiding shock instabilities, apart from being accurate in capturing contact discontinuities, do not need wave speed corrections and are independent of eigen-struture of the underlying hyperbolic parts of the systems.

TRINITY: A three-dimensional time-dependent radiative transfer code for in-vivo near-infrared imaging

We develop a new three-dimensional time-dependent radiative transfer code, TRINITY (Time-dependent Radiative transfer In Near-Infrared TomographY), for in-vivo diffuse optical tomography (DOT). The simulation code is based on the design of long radiation rays connecting boundaries of a computational domain, which allows us to calculate light propagation with little numerical diffusion. We parallelize the code with Message Passing Interface (MPI) using the domain decomposition technique and confirm the high parallelization efficiency, so that simulations with a spatial resolution of ∼1 millimeter can be performed in practical time. As a first application, we study the light propagation for a pulse collimated within θ∼15∘ in a phantom, which is a uniform medium made of polyurethane mimicking biological tissue. We show that the pulse spreads in all forward directions over ∼ a few millimeters due to the multiple scattering process of photons. Our simulations successfully reproduce the time-resolved signals measured with eight detectors for the phantom. We also introduce the effects of reflection and refraction at the boundary of medium with a different refractive index and demonstrate the faster propagation of photons in an air hole that is an analogue for the respiratory tract.

Numerical Considerations in the Simulation of Equatorial Spread F

AbstractThe ionospheric irregularities known as equatorial spread F (ESF) have drawn great attention in decades due to their disruption of communications and navigation systems. ESF is generally attributed to Rayleigh–Taylor instability, but there are still many remaining issues about its formation and evolution. Numerical simulations have been widely used to study such a phenomenon due to the complex nonlinear evolution of ESF. A two‐dimensional ESF model was developed in this work and a series of experiments were conducted to study the effects of grid resolution and numerical diffusion on idealized ESF simulation. It was found that advection schemes with the total variation diminishing property are preferred for preventing spurious oscillations which occur near the steep density gradients in the ESF. On the other hand, schemes with low diffusion are also desirable in order to model complicated dynamic structures of the ESF. Moreover, numerical experiments suggest that some of the secondary instabilities in idealized ESF simulations, although with morphological similarities, are possibly initiated by numerical seeding, which may differ from the observed evolution of secondary ESF.

Turbulence Diffusion and Energy Spectrum Analysis of Large-scale Current Sheets in Flares

Magnetic reconnection (MR) is one of the most important physical processes for many dynamical phenomena in the universe. Magnetohydrodynamical (MHD) simulation is an effective way to study the MR process and the physical pictures related to the MR. With different parameter setups, we investigate the influences of the Magnetic Reynolds number and spatial resolution on the reconnection rate, numerical dissipation, and energy spectrum distribution in the MHD simulation. We have found that the magnetic Reynolds number Rm has definite impact on the reconnection rate and energy spectrum distribution. The characteristic time for entering into the non-linear phase will be earlier as the Reynolds number increases. When it comes to the tearing phase, the reconnection rate will increase rapidly. On the other hand, the magnetic Reynolds number affects significantly the Kolmogorov microscopic scale lko, which becomes smaller as Rm increases. An extra dissipation is defined as the combined effect of the numerical diffusion and turbulence dissipation. It is shown that the extra dissipation is dominated by the numerical diffusion before the tearing mode instability takes place. After the instability develops, the extra dissipation rises vastly, which indicates that turbulence caused by the instability can enhance the diffusion obviously. Furthermore, the energy spectrum analysis indicates that lko of the large-scale current sheet may appear at a macroscopic MHD scale very possibly.

Diffusivity in force-free simulations of global magnetospheres

ABSTRACT Assuming that the numerical diffusivity triggered by violations of the force-free electrodynamics constraints is a proxy for the physical resistivity, we examine its impact on the overall dynamics of force-free aligned pulsar magnetospheres endowed with an equatorial current sheet (ECS). We assess the constraint violations as a diffusivity source. The effects of modifications on electric fields used to restore force-free conditions are not confined to the ECS, but modify the magnetospheric dynamics on time-scales shorter than the pulsar rotational period. These corrections propagate especially via a channel that was unexplored, namely, changes induced to the electric charge density, ρ. We quantify the global consequences of diffusivity by comparing different techniques to model ρ. By default, we combine a conservative ρ-evolution with hyperbolic/parabolic cleaning of inaccuracies in the Maxwell equations. As an alternative, we enforce a constrained evolution, where ρ is directly computed as the electric field divergence. The conservative approach reduces the Poynting flux dissipated in the ECS by an order of magnitude, along with an increase of the pulsar luminosity driven by a shift of the Y-point location. The luminosity changes according to LY ∝ α0.11, where α is the ratio of diffusion to advection time-scales, controlling the amount of (numerical) diffusivity. Our models suggest interpreting the luminosity dependence on the Y-point location as differences in resistivities encountered at the ECS. Alternatively, they could be interpreted in terms of the pair formation multiplicity, κ, smaller diffusion being consistent with κ ≫ 1.

Large-amplitude longitudinal oscillations in solar prominences simulated with different resolutions

Context. Large-amplitude longitudinal oscillations (LALOs) in solar prominences have been widely studied in recent decades. However, their damping and amplification mechanisms are not well understood. Aims. In this study, we investigate the attenuation and amplification of LALOs using high-resolution numerical simulations with progressively increasing spatial resolutions. Methods. We performed time-dependent numerical simulations of LALOs using the 2D magnetic configuration that contains a dipped region. After the prominence mass loading in the magnetic dips, we triggered LALOs by perturbing the prominence mass along the magnetic field. We performed the experiments with four values of spatial resolution. Results. In the simulations with the highest resolution, the period shows good agreement with the pendulum model. The convergence experiment revealed that the damping time saturates at the bottom prominence region with increasing resolution, indicating the existence of a physical reason for the damping of oscillations. At the prominence top, the oscillations are amplified during the first minutes and are then slowly attenuated. The characteristic time suggests more significant amplification in the experiments with the highest spatial resolution. The analysis revealed that the energy exchange between the bottom and top prominence regions is responsible for the attenuation and amplification of LALOs. Conclusions. High-resolution experiments are crucial when studying the periods and the damping mechanism of LALOs. The period agrees with the pendulum model only when using a sufficiently high spatial resolution. The results suggest that numerical diffusion in simulations with insufficient spatial resolution can hide important physical mechanisms, such as amplification of oscillations.

Spatial convergence study of Lax-Friedrichs WENO fast sweeping method on the SN transport equation with nonsmoothness

Spatial convergence study of the Lax-Friedrichs fast sweeping method based on the classical third order WENO scheme (WENO3) without a priori smoothness indicators is performed on the variants of Larsen’s benchmark problem with nonsmooth solutions. For comparison, third order upwind (UPWD3), weighted diamond difference, step method, bilinear discontinuous finite element method are also included. In optically thick regimes of the cases with continuous angular fluxes, WENO3 detects steep gradients of angular fluxes and handles them as discontinuities to suppress the unphysical oscillation, which improves the accuracy significantly compared with UPWD3. In the cases with discontinuous angular fluxes, WENO3 achieves non-oscillatory solutions with limited numerical diffusion. Nonlinear Gauss-Seidel iteration is necessary to obtain convergent solution for the nonlinearity of WENO3, which leads to low computational efficiency. The convergence behavior of nonlinear Gauss-Seidel iteration is affected by nonsmoothness of the exact solution, mesh number and the relaxation parameter.