Subcell Resolution Research Articles

High order finite difference methods with subcell resolution for advection equations with stiff source terms

A new high order finite-difference method utilizing the idea of Harten ENO subcell resolution method is proposed for chemical reactive flows and combustion. In reaction problems, when the reaction time scale is very small, e.g., orders of magnitude smaller than the fluid dynamics time scales, the governing equations will become very stiff. Wrong propagation speed of discontinuity may occur due to the underresolved numerical solution in both space and time. The present proposed method is a modified fractional step method which solves the convection step and reaction step separately. In the convection step, any high order shock-capturing method can be used. In the reaction step, an ODE solver is applied but with the computed flow variables in the shock region modified by the Harten subcell resolution idea. For numerical experiments, a fifth-order finite-difference WENO scheme and its anti-diffusion WENO variant are considered. A wide range of 1D and 2D scalar and Euler system test cases are investigated. Studies indicate that for the considered test cases, the new method maintains high order accuracy in space for smooth flows, and for stiff source terms with discontinuities, it can capture the correct propagation speed of discontinuities in very coarse meshes with reasonable CFL numbers.

On a compact non-extrapolating scheme for adaptive image interpolation

The motivation of this paper is to improve the results of ENO subcell resolution, that already obtains really good results but using extrapolation, as was stated by F. Arándiga, R. Donat and P. Mulet [Adaptive Interpolation of Images, Signal Processing 83(2) (2003) 459–464]. We present a new method for image interpolation which combines a new nonlinear cell-average interpolatory technique that uses a trigonometric mean and subcell resolution. The scheme uses a fully compact reconstruction and avoid any step of extrapolation. The experiments presented validate the theoretical results obtained.

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

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

Seismic waves in heterogeneous material: subcell resolution of the discontinuous Galerkin method

SUMMARY We present an important extension of the arbitrary high-order discontinuous Galerkin (DG) finite-element method to model 2-D elastic wave propagation in highly heterogeneous material. In this new approach we include space-variable coefficients to describe smooth or discontinuous material variations inside each element using the same numerical approximation strategy as for the velocity–stress variables in the formulation of the elastic wave equation. The combination of the DG method with a time integration scheme based on the solution of arbitrary accuracy derivatives Riemann problems still provides an explicit, one-step scheme which achieves arbitrary high-order accuracy in space and time. Compared to previous formulations the new scheme contains two additional terms in the form of volume integrals. We show that the increasing computational cost per element can be overcompensated due to the improved material representation inside each element as coarser meshes can be used which reduces the total number of elements and therefore computational time to reach a desired error level. We confirm the accuracy of the proposed scheme performing convergence tests and several numerical experiments considering smooth and highly heterogeneous material. As the approximation of the velocity and stress variables in the wave equation and of the material properties in the model can be chosen independently, we investigate the influence of the polynomial material representation on the accuracy of the synthetic seismograms with respect to computational cost. Moreover, we study the behaviour of the new method on strong material discontinuities, in the case where the mesh is not aligned with such a material interface. In this case second-order linear material approximation seems to be the best choice, with higher-order intra-cell approximation leading to potential instable behaviour. For all test cases we validate our solution against the well-established standard fourth-order finite difference and spectral element method.

On reinitializing level set functions

In this paper, we consider reinitializing level functions through equation ϕt+sgn(ϕ0)(‖∇ϕ‖-1)=0 [16]. The method of Russo and Smereka [11] is taken in the spatial discretization of the equation. The spatial discretization is, simply speaking, the second order ENO finite difference with subcell resolution near the interface. Our main interest is on the temporal discretization of the equation. We compare the three temporal discretizations: the second order Runge–Kutta method, the forward Euler method, and a Gauss–Seidel iteration of the forward Euler method. The fact that the time in the equation is fictitious makes a hypothesis that all the temporal discretizations result in the same result in their stationary states. The fact that the absolute stability region of the forward Euler method is not wide enough to include all the eigenvalues of the linearized semi-discrete system of the second order ENO spatial discretization makes another hypothesis that the forward Euler temporal discretization should invoke numerical instability. Our results in this paper contradict both the hypotheses. The Runge–Kutta and Gauss–Seidel methods obtain the second order accuracy, and the forward Euler method converges with order between one and two. Examining all their properties, we conclude that the Gauss–Seidel method is the best among the three. Compared to the Runge–Kutta, it is twice faster and requires memory two times less with the same accuracy.

Tracing exogenous Gd and its effects in single Chang cells

X-ray fluorescence (XRF) is rapidly becoming one of the tools of choice in tracing the presence of both endogenous and exogenous chemical elements in biological samples. The sub-micron spatial resolution routinely obtainable at multi-keV energies at third generation light sources, combined with the high brilliance of the photon beam, allows mapping the presence of biologically relevant elements at sub-cell resolution by means of their fluorescent signature. The fluorescent signal also lends itself for a semi-quantitative analysis of the elements composing the specimen. In this work we employed synchrotron-based XRF to analyze two lines of cultured Chang cells. One of them was treated with a Gd-containing magnetic resonance imaging (MRI) contrast agents (CA), and the other was left untreated for control purposes. The experiments highlighted a peri-nuclear distribution of Gd inside the cells, as well as a distinct variation in the distribution and concentration of several elements (P, S, Cl, K, Ca and Fe), with respect to the control line.

A spectrally refined interface approach for simulating multiphase flows

This paper presents a novel approach to phase-interface transport based on pseudo-spectral sub-grid refinement of a level set function. In each flow solver grid cell, a set of quadrature points is introduced on which the value of the level set function is known. This methodology allows to define a polynomial reconstruction of the level set function in each cell. The transport is performed using a semi-Lagrangian technique, removing all constraints on the time step size. Such an approach provides sub-cell resolution of the phase-interface and leads to excellent accuracy in the transport, while a reasonable cost is obtained by pre-computing some of the metrics associated with the polynomials. To couple this approach with a flow solver, an converging curvature computation is introduced. First, a second order explicit distance to the sub-grid interface is reconstructed on the flow solver mesh. Then, a least squares approach is employed to extract the curvature from this distance function. This technique is found to combine the high accuracy and good conservation found in the particle level set method with the converging curvature usually obtained with classical high order PDE transport of the level set function. Tests are presented for both transport as well as two-phase flows, that suggest that this technique is capable of retaining the thin liquid structures that are expected in turbulent atomization of liquids.

Approximation of piecewise smooth functions and images by edge-adapted (ENO-EA) nonlinear multiresolution techniques

This paper introduces and analyzes new approximation procedures for bivariate functions. These procedures are based on an edge-adapted nonlinear reconstruction technique which is an intrinsically two-dimensional extension of the essentially non-oscillatory and subcell resolution techniques introduced in the one-dimensional setting by Harten and Osher. Edge-adapted reconstructions are tailored to piecewise smooth functions with geometrically smooth edge discontinuities, and are therefore attractive for applications such as image compression and shock computations. The local approximation order is investigated both in L p and in the Hausdorff distance between graphs. In particular, it is shown that for general classes of piecewise smooth functions, edge-adapted reconstructions yield multiscale representations which are optimally sparse and adaptive approximations with optimal rate of convergence, similar to curvelets decompositions for the L 2 error.

Accurate Calculation of Fault-Rupture Models Using the High-Order Discontinuous Galerkin Method on Tetrahedral Meshes

We present a new method for near-source ground-motion calculations due to earthquake rupture on potentially geometrically complex faults. Following the recently introduced Discontinuous Galerkin approach with local time stepping on tetrahedral meshes, we use piecewise polynomial approximations of the unknown variables inside each element and achieve the same approximation order in time and space due to the new ader time integration scheme that uses Arbitrary high-order DERivatives. We show how an external source term and its heterogeneous properties in space and time, given by a fine discretization of an extended rupture surface, can be included in much coarser tetrahedral meshes due to the subcell resolution of the high-order polynomial representation. Hereby, the rupture surface is represented very generally as a point cloud of the center locations of individual subfaults at which each polynomial test function is evaluated exactly inside an element and the space– time integration of the source term is accurately computed at each timestep. Besides the incorporation of complex source kinematics we also present the effects of model boundaries that can degrade the accuracy of seismograms due to weak artificial reflections. We propose an extended computational domain of a coarsely meshed buffer region and show that our scheme using the local timestepping completely avoids such boundary problems with only slightly increasing the computational cost. We validate the new approach against different test cases, comparing our results with analytic, quasianalytic, and a series of reference solutions. Our work shows that adding the functionality of accurately treating finite source-rupture models into the general framework of the ader-Discontinuous Galerkin approach is an important contribution to modeling realistic earthquake scenarios, allowing the efficient inclusion of heterogeneous source kinematics and complex rupture-surface geometries in near-source ground-motion simulations.

Free and smooth boundaries in 2-D finite-difference schemes for transient elastic waves

A method is proposed for accurately describing arbitrary-shaped free boundaries in single-grid finite-difference schemes for elastodynamics, in a time-domain velocity-stress framework. The basic idea is as follows: fictitious values of the solution are built in vacuum, and injected into the numerical integration scheme near boundaries. The most original feature of this method is the way in which these fictitious values are calculated. They are based on boundary conditions and compatibility conditions satisfied by the successive spatial derivatives of the solution, up to a given order that depends on the spatial accuracy of the integration scheme adopted. Since the work is mostly done during the preprocessing step, the extra computational cost is negligible. Stress-free conditions can be designed at any arbitrary order without any numerical instability, as numerically checked. Using 10 grid nodes per minimal S-wavelength with a propagation distance of 50 wavelengths yields highly accurate results. With 5 grid nodes per minimal S-wavelength, the solution is less accurate but still acceptable. A subcell resolution of the boundary inside the Cartesian meshing is obtained, and the spurious diffractions induced by staircase descriptions of boundaries are avoided. Contrary to what occurs with the vacuum method, the quality of the numerical solution obtained with this method is almost independent of the angle between the free boundary and the Cartesian meshing.

Interpolation and Approximation of Piecewise Smooth Functions

This paper provides approximation orders for a class of nonlinear interpolation procedures for uniformly sampled univariate data. The interpolation is based on essentially nonoscillatory (ENO) and subcell resolution (SR) reconstruction techniques. These nonlinear techniques aim at reducing significantly the approximation error for functions with isolated singularities and are therefore attractive for applications such as shock computations or image compression. We prove that in the presence of isolated singularities, the approximation order provided by the interpolation procedure is improved by a factor of h relative to the linear methods, where h is the sampling rate. Moreover, for h below a critical value, we recover the optimal approximation order as for uniformly smooth functions.

The Explicit Simplified Interface Method for Compressible Multicomponent Flows

This paper concerns the numerical approximation of the Euler equations for multicomponent flows. A numerical method is proposed to reduce spurious oscillations that classically occur around material interfaces. It is based on the explicit simplified interface method (ESIM), previously developed in the linear case of acoustics with stationary interfaces [J. Piraux and B. Lombard, J. Comput. Phys., 168 (2001), pp. 227--248]. This technique amounts to a higher-order extension of the ghost fluid method introduced in Euler multicomponent flows [R. P. Fedkiw, T. Aslam, B. Merriman, and S. Osher, J. Comput. Phys., 152 (1999), pp. 457--492]. The ESIM is coupled to sophisticated shock-capturing schemes for time-marching, and to level-sets for tracking material interfaces. Jump conditions satisfied by the exact solution and by its spatial derivative are incorporated in numerical schemes, ensuring a subcell resolution of material interfaces inside the meshing. Numerical experiments show the efficiency of the method for rich-structured flows.

Numerical simulation of thermo‐solutal‐capillary migration of a dissolving drop in a cavity

AbstractIn the present paper the thermo‐solutal‐capillary migration of a dissolving liquid drop, composed by a binary mixture having a miscibility gap, injected in a closed cavity with differentially heated end walls, is studied. The main goal of the analysis is to clarify if and how the drop migration is affected by the dissolution process. The numerical code is based on a finite volume formulation. A level‐set technique is used for describing the dynamics of the interface separating the different phases. A thermodynamic constraint fixes the concentration jump between the interface sides. This jump, together with that of the concentration normal derivatives, in turn defines the entity of the dissolution cross‐flow through the interface and the interface velocity relative to the fluid. Since the jump singularity of normal derivatives cannot be easily mollified, while retaining the necessary accuracy, a scheme for the species equation is elaborated that allows sharp jumps and has subcell resolution. Steady migration speeds are determined after the start‐up phase for different radii and temperature differences. The results will be used for the preparation of a sounding rocket space experiment. Copyright © 2003 John Wiley & Sons, Ltd.

The next generation CIP as a conservative semi-Lagrangian solver for solid, liquid and gas

We present a review of the CIP method, which is a kind of semi-Lagrangian scheme and has been extended to treat incompressible flow in the framework of compressible fluid. Since it uses primitive Euler representation, it is suitable for multi-phase analysis. The recent version of this method guarantees the exact mass conservation even in the framework of semi-Lagrangian scheme. Comprehensive review is given for the strategy of the CIP method that has a compact support and subcell resolution including front capturing algorithm with functional transformation.

The Constrained Interpolation Profile Method for Multiphase Analysis

We present a review of the constrained interpolation profile (CIP) method that is known as a general numerical solver for solid, liquid, gas, and plasmas. This method is a kind of semi-Lagrangian scheme and has been extended to treat incompressible flow in the framework of compressible fluid. Since it uses primitive Euler representation, it is suitable for multiphase analysis. The recent version of this method guarantees the exact mass conservation even in the framework of a semi-Lagrangian scheme. We provide a comprehensive review of the strategy of the CIP method, which has a compact support and subcell resolution, including a front-capturing algorithm with functional transformation, a pressure-based algorithm, and other miscellaneous physics such as the elastic–plastic effect and surface tension. Some practical applications are also reviewed, such as milk crown or coronet, laser-induced melting, and turbulent mixing layer of liquid–gas interface.

A spectral method with subcell resolution for shock wave calculations

A spectral method for the numerical solution of hyperbolic partial differential equations is presented. Any discontinuous solutions which may occur are assumed to be the sum of a step function and a smooth function. A series approximation is then applied to the smooth function in order to eliminate the Gibbs phenomenon. By adding the step function onto the series we then have an accurate approximation to the solution at any given time. When discretising in time, a modification is added to the numerical flux to account for the advection of sharp discontinuities across cells. We evaluate the method by its application to three standard test problems involving the scalar wave equation, the inviscid Burgers equation and the Euler equations of gas dynamics. For all experiments we observe an absence of the Gibbs phenomenon with discontinuities captured to within a single mesh interval and high accuracy is observed in smooth regions.

A Comparison of ENO and TVD Schemes for the Computation of Shock-Turbulence Interaction

A quantitative analysis of solutions to the Euler equations of fluid dynamis with the MUSCL, ENO-Harten, and efficient ENO-Shu algorithms is performed. Investigations of different test problems in one and two dimensions are presented. These are chosen as to model the shock-turbulence interaction in fluid dynamical systems. The notion of subcell resolution developed by Harten for the ENO schemes clearly improves the solution in one dimension; however, the effect is less prominent in a Strang-type extension to two dimensions. Our results confirm the superiority of the ENO schemes over the MUSCL approach in solving problems of flow fields with discontinuities which, at the same time, contain fine structure in its smooth parts.

Flux difference splittings and limiters for the resolution of contact discontinuities

Harten, Lax, and van Leer developed a theory of approximate Riemann solutions and used it to derive numerical methods for systems of conservation laws. In the first part of this paper, we reformulate this theory so that it can be used to compare various flux difference splitting formulae. In particular we show how to compare the resolution of contact discontinuities by various methods with that of the Roe method. Since the Roe method can resolve an isolated steady contact perfectly, it is a good standard for comparison. In the second part of the paper, we show that Harten's subcell resolution method is the ultimate compressive limiter. We discuss the differences between compressive and noncompressive limiters and show why it is important to choose the limiter to suit the problem at hand.

A subcell resolution method for viscous systems of conservation laws

We consider the generalization of scalar subcell resolution schemes to systems of viscous conservation laws. For this purpose we use a weakly nonlinear geometrical optics approximation for parabolic perturbations of hyperbolic conservation laws and the Roe-type field by field decomposition. Computations of the reactive Navier-Stokes equations are presented as an application.

Numerical modeling of an advanced pneumatic impulse ice protection system for aircraft

The development of a numerical model of an advanced pneumatic impulse ice protection system for aircraft is described in this paper. A time-dependent, compressible flow model for internal duct flow is used to model the ice protection system, also known as PUP. The model incorporates a high resolution shock capturing method, essentially nonoscillatory (ENO) scheme with subcell resolution (SR) and characteristic direction (CD). The model consists of a constant pressure reservoir, a variable area valve, and a tube section that expands due to the passage of a shockwave. A grid generation scheme is used to account for the size of the system. The model is analyzed for different inlet and exit boundary conditions. Preliminary results are obtained for a range of compressor pressures and the model is shown to give results of practical use. The paper also recommends a suitable operating pressure for the compressor.