Search

The Riemann problem for the shallow water equations with discontinuous topography is considered. In a general case the exact solution of this problem is not unique, which complicates the application of an exact Riemann solver in numerical methods, since it is not clear which solution should be chosen. In the present work it is shown that involving an additional physical assumption makes it possible to prove the existence and uniqueness of the solution. The assumption is that the discharge at the bottom discontinuity should continuously depend on the initial conditions. The proven uniqueness opens up a possibility to use an exact Riemann solver for a numerical solution of the shallow water equations with complex discontinuous topography.

A novel 5th-order shock capturing scheme is presented in this paper. The scheme, so-called P 4 T 2 − BVD (polynomial of 4-degree and THINC function of 2-level reconstruction based on BVD algorithm), is formulated as a two-stage spatial reconstruction scheme following the BVD (Boundary Variation Diminishing) principle that minimizes the jumps of the reconstructed values at cell boundaries. In the P 4 T 2 − BVD scheme, polynomial of degree four and THINC (Tangent of Hyperbola for INterface Capturing) functions with two-level steepness are used as the candidate reconstruction functions. The final reconstruction function is selected through the two-stage BVD algorithm so as to effectively control both numerical oscillation and dissipation. Spectral analysis and numerical verifications show that the P 4 T 2 − BVD scheme possesses the following desirable properties: 1) it effectively suppresses spurious numerical oscillation in the presence of strong shock or discontinuity; 2) it substantially reduces numerical dissipation errors; 3) it automatically retrieves the underlying linear 5th-order upwind scheme for smooth solution over all wave numbers; 4) it is able to resolve both smooth and discontinuous flow structures of all scales with substantially improved solution quality in comparison to other existing methods; and 5) it produces accurate solutions in long term computation. P 4 T 2 − BVD , as well as the underlying idea presented in this paper, provides an innovative and practical approach to design high-fidelity numerical schemes for compressible flows involving strong discontinuities and flow structures of wide range scales.

A novel Eulerian approach is proposed for numerical simulations of wave propagation in viscoelastic media, for application to shocks interacting with interfaces between fluids and solids. We extend the five-equations multiphase interface-capturing model, based on the idea that all the materials (gases, liquids, solids) obey the same equation of state with spatially varying properties, to incorporate the desired constitutive relation; in this context, interfaces are represented by discontinuities in material properties. We consider problems in which the deformations are small, such that the substances can be described by linear constitutive relations, specifically, Maxwell, Kelvin–Voigt or generalized Zener models. The main challenge lies in representing the combination of viscoelastic, multiphase and compressible flow. One particular difficulty is the calculation of strains in an Eulerian framework, which we address by using a conventional hypoelastic model in which an objective time derivative (Lie derivative) of the constitutive relation is taken to evolve strain rates instead. The resulting eigensystem is analyzed to identify wave speeds and characteristic variables. The spatial scheme is based on a solution-adaptive formulation, in which a discontinuity sensor discriminates between smooth and discontinuous regions. To compute the convective fluxes, explicit high-order central differences are applied in smooth regions, while a high-order finite-difference Weighted Essentially Non-Oscillatory (WENO) scheme is used at discontinuities (shocks, material interfaces and contacts). The numerical method is verified in a comprehensive fashion using a series of smooth and discontinuous (shocks and interfaces), one- and two-dimensional test problems.

In Fu et al. [23] , a family of high-order targeted ENO (TENO) schemes with a new nonlinear weighting strategy has been proposed. Building upon this strategy, in the current paper, a new class of adaptive TENO schemes for hyperbolic conservation laws is proposed based on three new concepts: (1) a hierarchical voting strategy is proposed to improve the ENO-like stencil selection; (2) the TENO weighting strategy is extended to function as a built-in discontinuity-location detector. Since the reconstruction scheme must not cross any discontinuity, corresponding target schemes are selected from a set of predefined linear schemes which are optimized towards maximum accuracy order or spectral resolution; (3) based on the observation that the cut-off parameter C T in the TENO weighting strategy determines nonlinear dissipation, a C T adaptation strategy is developed to minimize numerical dissipation for high-wavenumber fluctuations while maintaining robustness for shock-capturing. Six-point and eight-point TENO schemes are constructed, and their spectral properties are analyzed by the ADR analysis. A set of benchmark cases is considered to demonstrate the performance of proposed adaptive TENO schemes. Numerical results suggest that the proposed TENO schemes preserve the accuracy order at first and second order critical points and show less numerical dissipation compared with typical WENO schemes. Moreover, the TENO8-NA scheme exhibits good results for both highly compressible flows and nearly incompressible turbulence.

We present in this work a new reconstruction scheme, so-called MUSCL-THINC-BVD scheme, to solve the five-equation model for interfacial two phase flows. This scheme employs the traditional shock capturing MUSCL (Monotone Upstream-centered Schemes for Conservation Law) scheme as well as the interface sharpening THINC (Tangent of Hyperbola for INterface Capturing) scheme as two building-blocks of spatial reconstruction on the BVD (boundary variation diminishing) principle that minimizes the variations (jumps) of the reconstructed variables at cell boundaries, and thus effectively reduces the dissipation error in numerical solutions. The MUSCL-THINC-BVD scheme is implemented to the volume fraction and other state variables under the same finite volume framework, which realizes the consistency among volume fraction and other physical variables. Numerical results of benchmark tests show that the present method is able to capture the material interface as a well-defined sharp jump in volume fraction, and obtain numerical solutions of superior quality in comparison to other existing methods. The proposed scheme is a simple and effective method of practical significance for simulating compressible interfacial multiphase flows.

A new exact Riemann solver for the shallow water equations with a discontinuous bottom is proposed. The algorithm is based on the approach to overcome the non-uniqueness of the Riemann problem solution by assuming that the ‘true’ solution should have the discharge at the bottom discontinuity, which continuously depends on the initial conditions. The solver ensures the existence and uniqueness of a solution for arbitrary initial conditions. The exact Riemann solver is embedded in the Godunov scheme for numerical solution of the shallow water equations to demonstrate its advantage. It is shown that the proposed solver allows one to significantly reduce the number of computational cells for stationary flows and non-stationary ones if the mesh can resolve non-stationary features of the flow. In addition, the practical problem of rainfall-runoff is considered, and the results obtained using the exact Riemann solver show good agreement with observation data.

A discontinuity sensor for shock waves and material interfaces is proposed to construct a low dissipative scheme that consists of a weighted compact finite nonlinear scheme (WCNS) and central difference to address compressible multi-component flows. The Harten-Lax-van Leer and contact (HLLC) scheme is applied to capture shock waves through material interfaces. At the material interface, the overestimated quasi-conservative WCNS maintains the equilibriums of velocity, pressure, and temperature. The discontinuity sensor is composed of a Larrson shock sensor and a material interface sensor to distinguish between smooth and discontinuous regions, respectively. In the single-component shock-vortex interaction problem, the present scheme successfully captures shock waves with low numerical dissipation. In the multi-component shock-bubble interaction problem, Richtmyer-Meshkov instability problem, and triple-point problem, the simultaneous capturing of shock waves and material interfaces by the present scheme is demonstrated. The proposed scheme suppresses the numerical oscillation by shock waves and material interfaces; additionally, it successfully reduces the numerical dissipation in smooth regions.

In this work, a modified Particle Method in fluid mechanics involving irregular distribution and discontinuity issues is proposed. Due to the approximation accuracy, momentum conservation and governing equation feasibility, traditional Particle Method, Moving Particle Semi-implicit (MPS) method in this paper, cannot guarantee the accuracy and stability of computation, and the divergence calculation in Navier–Stokes equation may cause numerical error in computation domain with discontinuity. To enhance the computation accuracy of Particle Method, we modify the gradient operator with a corrected tensor emanating from minimizing the local error of the first-order Taylor Expansion approximation. Besides, inspired by Peridynamic which is mostly used in solid mechanics, a new governing equation in an integral form which maintains the momentum conservation is obtained. Combining the modified gradient and new governing equation, we obtain a new particle-based method. The numerical results verify its feasibility and illustrate a good computational performance of proposed Method in the fluid dynamics involving non-uniform particle distribution.

This study develops a new Lagrangian particle method for modeling flow and transport phenomena in complex porous media with discontinuities. For instance, diffusion processes can be modeled by Lagrangian Random Walk algorithms. However, discontinuities and heterogeneities are difficult to treat, particularly discontinuous diffusion D ( x ) or porosity θ ( x ) . In the literature on particle Random Walks, previous methods used to handle this discontinuity problem can be characterized into two main classes as follows: “Interpolation techniques”, and “Partial reflection methods”. One of the main drawbacks of these methods is the small time step required in order to converge to the expected solution, particularly in the presence of many interfaces. These restrictions on the time step, lead to inefficient algorithms. The Random Walk Particle Tracking (RWPT) algorithm proposed here is, like others in the literature, discrete in time and continuous in space (gridless). We propose a novel approach to partial reflection schemes without restrictions on time step size. The new RWPT algorithm is based on an adaptive “Stop&Go” time-stepping, combined with partial reflection/refraction schemes, and extended with a new concept of negative mass particles. To test the new RWPT scheme, we develop analytical and semi-analytical solutions for diffusion in the presence of multiple interfaces (discontinuous multi-layered medium). The results show that the proposed Stop&Go RWPT scheme (with adaptive negative mass particles) fits extremely well the semi-analytical solutions, even for very high contrasts and in the neighborhood of interfaces. The scheme provides a correct diffusive solution in only a few macro-time steps, with a precision that does not depend on their size.

This paper presents a single-pressure-field two-fluid model with finite-volume discretization to solve the equations of motion of compressible multiphase flows. To capture the discontinuities caused by shock waves and fluid interfaces, we propose a generalized discontinuity sharpening technique that combines the conventional monotonic upstream scheme for conservation law (MUSCL) and tangent of hyperbola interface capturing (THINC) schemes. In addition, a slope ratio-weighted parameter, ζ, is used to control the proportion of values reconstructed by MUSCL and THINC, and we show that the present method can retain sharp interfaces when the value of the parameter β in the THINC scheme is set ranging from 1.6 to 3.0. Fluxes across various interfaces are evaluated using a hybrid AUSMD-type flux algorithm, where the mass flux and pressure induced on the cell faces are calculated using an approximate Riemann solver. The accuracy and robustness of the proposed method are validated by solving a series of one- and two-dimensional single-phase flows. Furthermore, complex wave patterns arising from two-dimensional shock bubble/water-column interactions are examined, which indicate that compared with the existing schemes applied to two-fluid modeling, the proposed scheme significantly sharpens the interfaces and captures more details of the flow features. Finally, simulations of a three-dimensional example of the liquid jet crossflow are conducted. The proposed scheme shows more details of the fluid interface, including the interfacial instabilities on the windward side of the liquid jet and droplet formation due to the breakup phenomenon in the downstream of the crossflow, than the existing schemes.

Walsh functions form an orthonormal basis set consisting of square waves. Square waves make the system well suited for detecting and representing functions with discontinuities. Given a uniform distribution of 2p cells on a one-dimensional element, it is proved that the inner product of the Walsh Root function for group p with every polynomial of degree (p1) across the element is identically zero. It is also proved that the magnitude and location of a discontinuous jump, as represented by a Heaviside function, are explicitly identified by its Fast Walsh Transform (FWT) coefficients. These two proofs enable an algorithm that quickly provides a Weighted Least Squares fit to distributions across the element that include a discontinuity. It is shown that flux reconstruction relative to the FWT fit in partial differential equations provides improved accuracy. The detection of a discontinuity further enables analytic relations to locally describe its evolution and provide increased accuracy. Examples are provided for time-accurate advection, Burgers' equation, and quasi-one-dimensional nozzle flow.

A 2D simulation method for a rigid body moving in an incompressible viscous fluid is proposed. It combines one of the immersed-boundary methods, the DFFD (direct forcing fictitious domain) method with the spectral element method; the former is employed for efficiently capturing the two-way FSI (fluid-structure interaction) and the geometric flexibility of the latter is utilized for any possibly co-existing stationary and complicated solid or flow boundary. A pseudo body force is imposed within the solid domain to enforce the rigid body motion and a Lagrangian mesh composed of triangular elements is employed for tracing the rigid body. In particular, a so called sub-cell scheme is proposed to smooth the discontinuity at the fluid-solid interface and to execute integrations involving Eulerian variables over the moving-solid domain. The accuracy of the proposed method is verified through an observed agreement of the simulation results of some typical flows with analytical solutions or existing literatures. The direct-forcing fictitious domain method and the spectral element method are combined.Direct interpolation of Eulerian variables onto Lagrangian mesh is not necessary.Non-fluctuating pressure fields are obtained with a smoothed pseudo body force.The proposed method possesses a spectral accuracy.Both small and large solid-to-fluid density ratios are applicable.

The nonlinear filter introduced by Yee et al. (1999) [27] and extensively used in the development of low dissipative well-balanced high order accurate finite-difference schemes is adapted to the finite element context of discontinuous Galerkin (DG) discretizations. The filter operator is constructed in the canonical computational domain for the standard cubical element where it is applied to the computed conservative variables in a direction per direction basis. Filtering becomes possible for all element types in unstructured meshes using collapsed coordinate transformations. The performance of the proposed nonlinear filter for DG discretizations is demonstrated and evaluated for different orders of expansions for one-dimensional and multidimensional problems with exact solutions. It is shown that for higher order discretizations discontinuity resolution within the cell is achieved and the design order of accuracy is preserved. The filter is applied for a number of standard inviscid flow test problems including strong shocks interactions to demonstrate that the proposed dissipative mechanism for DG discretizations yields superior results compared to the results obtained with the total variation bounded (TVB) limiter and high-order hierarchical limiting. The proposed approach is suitable for p-adaptivity in order to locally enhance resolution of three-dimensional flow simulations that include discontinuities and complex flow features.

The paper describes a flexible computational method that can be used to represent the complex geometric evolution of a fluid-driven fracture front using an unstructured triangular element mesh. A specific motivation for the method is to allow subsequent treatment of non-planar fracture growth problems. Elastic interactions between the crack opening displacements and the surrounding medium are calculated using displacement discontinuity boundary element influence functions. A novel feature of the approach is the use of an adjustable region of tip elements to accommodate both time-dependent crack edge movements and to ensure correct velocity-dependent asymptotic behavior for viscous and viscous-toughness flow conditions. The adjustable tip element region is represented using a special-purpose data structure that allows the iteration of the moving edge element vertices to be linked directly to the flow rate and the crack opening solution values that are determined in each time step. The fringe region is periodically replaced with a set of fixed vertex elements as the fracture surface evolves. The proposed scheme has been evaluated using available analytic and experimental results for planar fracture propagation. The method is found to be both accurate and robust for a range of choices of geometric fracture shapes and time step sizes.

To enhance the accuracy and applicability of the zonal grid refinement method for the lattice Boltzmann method, a new method which minimizes the interface imbalances of mass and momentum is developed. This method introduces a correction step for the macroscopic flow variables such as the fluid density and velocity to remove their interface discontinuity. To demonstrate and evaluate the presently developed imbalance correction grid refinement method, large eddy simulations of turbulent channel and square cylinder flows are carried out. By changing the grid arrangement in the turbulent channel flows, it is confirmed that the present method reduces the sensitivity to the location of the grid refinement interface and minimizes the unphysically discontinuous profiles satisfactorily. Furthermore, the present method considerably improves mass conservation of the system, which is particularly important for long time periodical flow simulations. It is also confirmed that the present method generally improves the prediction performance of the square cylinder flows.

Scattering or coupling of electromagnetic beam-field at a surface discontinuity separating two homogeneous or inhomogeneous media with different propagation characteristics is formulated using surface integral equation, which are solved by the Method of Moments with the aid of the Gabor-based Gaussian window frame set of basis and testing functions. The application of the Gaussian window frame provides (i) a mathematically exact and robust tool for spatial-spectral phase-space formulation and analysis of the problem; (ii) a system of linear equations in a transmission-line like form relating mode-like wave objects of one medium with mode-like wave objects of the second medium; (iii) furthermore, an appropriate setting of the frame parameters yields mode-like wave objects that blend plane wave properties (as if solving in the spectral domain) with Green's function properties (as if solving in the spatial domain); and (iv) a representation of the scattered field with Gaussian-beam propagators that may be used in many large (in terms of wavelengths) systems.

Two modified segregated PISO algorithms are proposed, which are constructed to avoid the development of spurious oscillations in the computed flow near sharp interfaces of conjugate fluid-porous domains. The new collocated finite volume algorithms modify the Rhie-Chow interpolation to maintain a correct pressure-velocity coupling when large discontinuous momentum sources associated with jumps in the local permeability and porosity are present. The Re-Distributed Resistivity (RDR) algorithm is based on spreading flow resistivity over the grid cells neighboring a discontinuity in material properties of the porous medium. The Face Consistent Pressure (FCP) approach derives an auxiliary pressure value at the fluid-porous interface using momentum balance around the interface. Such derived pressure correction is designed to avoid spurious oscillations as would otherwise arise with a strictly central discretization. The proposed algorithms are successfully compared against published data for the velocity and pressure for two reference cases of viscous flow. The robustness of the proposed algorithms is additionally demonstrated for strongly reduced viscosity, i.e., higher Reynolds number flows and low Darcy numbers, i.e., low permeability of the porous regions in the domain, for which solutions without unphysical oscillations are computed. Both RDR and FCP are found to accurately represent porous media flow near discontinuities in material properties on structured grids. Two modified segregated PISO algorithms are proposed.The Rhie-Chow interpolation is modified to maintain a correct pressure-velocity coupling when large discontinuous momentum sources are present.The proposed algorithms are successfully compared against published data for the velocity and pressure for two reference cases of viscous flow.The robustness of the proposed algorithms is additionally demonstrated for high Reynolds number flows and low Darcy numbers.

High-order Flux Reconstruction (FR) schemes have been limited in their application to transonic and supersonic problems on account of numerical instabilities related to the resolution of jump discontinuities. These instabilities arise from aliasing errors associated with the collocation projection of the flux corresponding to the numerical solution onto the polynomial basis of the numerical flux. In this paper, we obtain energy bounds on the numerical solution via FR to prove that stability can be ensured for any polynomial order by the addition of adequate artificial dissipation such that the solution is energy-stable beyond a critical grid resolution. This artificial viscosity is then approximately posed as a Fourier filtering operation which is implemented in the physical space via a strictly local convolution integral. The filter is selectively applied to 'troubled' cells as indicated by a discontinuity sensor based on the spectral concentration method. Numerous numerical tests in 1-D and 2-D have been performed. The proposed approach captures shock discontinuities while preserving accuracy in smooth regions of the solution, even for very high polynomial orders such as P = 119 . The filtered solution provides reduced total variation, reduced maximum overshoot/undershoot, and even allows sub-element shocks to be localized in the interior of an element.

A sharp interface method (SIM) for smoothed particle hydrodynamics (SPH) has been developed to simulate two-phase flows with clear interfaces. The level set function is introduced to capture the interface implicitly. The interface velocity is used to evolve the level set function. The smoothness of the level set function helps to improve the accuracy of the interface curvature. Material discontinuity across the interface is dealt with by the ghost fluid method. The interface states are calculated by applying the jump conditions and are extended to the corresponding ghost fluid particles. The ghost fluid method helps to get smooth and stable calculation near the interface. The performance of the developed method is validated by benchmark tests. The developed SIM for SPH can be applied to simulate low speed two-phase flows of high density ratios with clear interface accurately and stably.

A set of numerical schemes is developed for two- and three-dimensional time-dependent 3-T radiation diffusion equations in systems involving multi-materials. To resolve sub-cell structure, interface reconstruction is implemented within any cell that has more than one material. Therefore, the system of 3-T radiation diffusion equations is solved on two- and three-dimensional polyhedral meshes. The focus of the development is on the fully coupling between radiation and material, the treatment of nonlinearity in the equations, i.e., in the diffusion terms and source terms, treatment of the discontinuity across cell interfaces in material properties, the formulations for both transient and steady states, the property for large time steps, and second order accuracy in both space and time. The discontinuity of material properties between different materials is correctly treated based on the governing physics principle for general polyhedral meshes and full nonlinearity. The treatment is exact for arbitrarily strong discontinuity. The scheme is fully nonlinear for the full nonlinearity in the 3-T diffusion equations. Three temperatures are fully coupled and are updated simultaneously. The scheme is general in two and three dimensions on general polyhedral meshes. The features of the scheme are demonstrated through numerical examples for transient problems and steady states. The effects of some simplifications of numerical schemes are also shown through numerical examples, such as linearization, simple average of diffusion coefficient, and approximate treatment for the coupling between radiation and material.

A special mesh adaptation technique and a precise discontinuity tracking are presented for an accurate, efficient, and robust adaptive-mesh computational procedure for one-dimensional hyperbolic systems of conservation laws, with particular reference to problems with evolving discontinuities in solids. The main advantage of the adaptive technique is its ability to preserve the modified mesh as close to the original fixed mesh as possible. The constructed method is applied to the martensitic phase-transition front propagation in solids.

This paper constructs a second-order large time step wave adding scheme (LTS-WA2) for hyperbolic conservation laws. Based on the first-order large time step wave adding scheme (LTS-WA1), a piece-wise linear reconstruction with limiter is performed on the solutions, and the band decomposition and band adding is complemented (into the discontinuity decomposition and wave adding), then the scheme is extended to second-order. The manufacture of the new scheme is simplified and written uniformly. Theoretical analyses are made for scalar cases. Numerical experiments give 11 tests which involve 1D scalar equations, 1D, 2D and 3D Euler equations. Computations show that the new scheme maintains high resolution and low dissipation of the original scheme at large CFL number, and also improves the problem of low resolution of the original scheme for rarefaction wave and contact discontinuity at small CFL number (roughly less than 2). Accuracy order is analyzed theoretically and validated using tests of 1D linear equation and 2D Burgers equation. CPU time is also compared with traditional second-order (Harten-TVD) scheme, showing that LTS-WA2 is more cost effective.

The immersed boundary–finite element method (IBFE) is an approach to describing the dynamics of an elastic structure immersed in an incompressible viscous fluid. In this formulation, there are generally discontinuities in the pressure and viscous stress at fluid–structure interfaces. The standard immersed boundary approach, which connects the Lagrangian and Eulerian variables via integral transforms with regularized Dirac delta function kernels, smooths out these discontinuities, which generally leads to low order accuracy. This paper describes an approach to accurately resolve pressure discontinuities for these types of formulations, in which the solid may undergo large deformations. Our strategy is to decompose the physical pressure field into a sum of two pressure–like fields, one defined on the entire computational domain, which includes both the fluid and solid subregions, and one defined only on the solid subregion. Each of these fields is continuous on its domain of definition, which enables high accuracy via standard discretization methods without sacrificing sharp resolution of the pressure discontinuity. Numerical tests demonstrate that this method improves rates of convergence for displacements, velocities, stresses, and pressures as compared to the conventional IBFE method. Further, it produces much smaller errors at reasonable numbers of degrees of freedom. The performance of this method is tested on several cases with analytic solutions, a nontrivial benchmark problem of incompressible solid mechanics, and an example involving a thick, actively contracting torus.

We present a framework for recovering/approximating unknown time-dependent partial differential equation (PDE) using its solution data. Instead of identifying the terms in the underlying PDE, we seek to approximate the evolution operator of the underlying PDE numerically. The evolution operator of the PDE, defined in infinite-dimensional space, maps the solution from a current time to a future time and completely characterizes the solution evolution of the underlying unknown PDE. Our recovery strategy relies on approximation of the evolution operator in a properly defined modal space, i.e., generalized Fourier space, in order to reduce the problem to finite dimensions. The finite dimensional approximation is then accomplished by training a deep neural network structure, which is based on residual network (ResNet), using the given data. Error analysis is provided to illustrate the predictive accuracy of the proposed method. A set of examples of different types of PDEs, including inviscid Burgers' equation that develops discontinuity in its solution, are presented to demonstrate the effectiveness of the proposed method., Minor notational changes

Lagrangian particles can be used either to analyze different complex flows, as massless tracers, or to model multiphase flows, as particles with mass. As particles can have arbitrary positions within the computational domain, interpolation of different flow quantities from the Eulerian grid is essential. Low-order centered interpolation schemes generally do not provide a sufficient level of accuracy for many flow configurations of interest. Thus, typically, high-order centered interpolation schemes are utilized. The current study demonstrates that for highly compressible non-reacting and reacting flow regimes, in which discontinuities in the flow field, e.g., shock waves arise, centered interpolation schemes tend to smear the shock and high order schemes also produce numerical oscillations. It is shown that this problem can be remedied by using Weighted-Essentially-Non-Oscillatory (WENO) interpolation schemes. Extensive numerical tests are performed in order to demonstrate this, comparing the performance of WENO-3 and WENO-5 interpolation schemes with a variety of centered interpolation schemes. The first test case involves specially designed steady state two-dimensional vortex flows. For a smooth flow, WENO schemes provide comparable accuracy to other high-order centered schemes. For flow fields with discontinuities, WENO-3 and WENO-5 interpolation schemes can decrease the interpolation error by more than two and three orders of magnitude, respectively, in comparison with centered schemes. Solution accuracy is further studied with a normal shock wave test. It is demonstrated that centered schemes tend to smear the shock, whereas, WENO schemes capture it in a much sharper manner. Moreover, high-order centered schemes tend to oscillate in the vicinity of the discontinuity leading to non-physical results, such as negative absolute pressure values. The same trends are retained for an unsteady three-dimensional spherical blast wave, where the shock is smeared across several grid cells, which is typical for shock-capturing flow solvers. Finally, the benefits of WENO schemes for the Lagrangian-particle tracking analysis of highly compressible reactive flows are explored by comparing various Lagrangian particle trajectories and joint probability density functions (PDFs) for a two-dimensional cellular detonation. The results indicate that high-order centered schemes lead to unphysical results, whereas WENO schemes can provide high-order interpolation that is free of severe numerical oscillations and non-physical artifacts, which is critical for the proper analysis of the flow.

This paper attempts to unify the asymptotic diffusion limit analysis of thermal radiation transport schemes, for a linear-discontinuous representation of the material temperature reconstructed from cell centred temperature unknowns, in a process known as ‘source tilting’. The asymptotic limits of both Monte Carlo (continuous in space) and deterministic approaches (based on linear-discontinuous finite elements) for solving the transport equation are investigated in slab geometry. The resulting discrete diffusion equations are found to have nonphysical terms that are proportional to any cell-edge discontinuity in the temperature representation. Based on this analysis it is possible to design accurate schemes for representing the material temperature, for coupling thermal radiation transport codes to a cell centred representation of internal energy favoured by ALE (arbitrary Lagrange–Eulerian) hydrodynamics schemes.

In a two-phase free-surface flow, the transition from a mobile-bed condition to a fixed-bed one (and vice versa) occurs at a sharp interface across which the relevant system of partial differential equations changes abruptly. This leads to the possibility of conceiving a new type of Riemann Problem (RP), which we have called Composite Riemann Problem (CRP), where not only the initial constant values of the variables but also the system of equations change from left to right of a discontinuity. In this paper, we present a strategy for solving a CRP by reducing it to a standard RP of a single, composite system of equations. This can be obtained by combining the two original systems by means of a suitable weighting function, namely the erodibility variable, and the introduction of an appropriate differential equation for this quantity. In this way, the CRP problem can be analyzed theoretically with standard methods, and the features of the solutions can be clearly identified. In particular, a stationary contact wave is able to correctly describe the sharp transition between mobile- and fixed-bed conditions. A finite volume scheme based on the Multiple Averages Generalized Roe approach (Rosatti and Begnudelli (2013) 22]) was used to numerically solve the fixed-mobile CRP. Several test cases demonstrate the effectiveness, exact well balanceness and high accuracy of the scheme when applied to problems that fall within the physical range of applicability of the relevant mathematical model.

A numerical scheme is developed for two- and three-dimensional time-dependent diffusion equations in numerical simulations involving mixed cells. The focus of the development is on the formulations for both transient and steady states, the property for large time steps, second-order accuracy in both space and time, the correct treatment of the discontinuity in material properties, and the handling of mixed cells. For a mixed cell, interfaces between materials are reconstructed within the cell so that each of resulting sub-cells contains only one material and the material properties of each sub-cell are known. Diffusion equations are solved on the resulting polyhedral mesh even if the original mesh is structured. The discontinuity of material properties between different materials is correctly treated based on governing physics principles. The treatment is exact for arbitrarily strong discontinuity. The formulae for effective diffusion coefficients across interfaces between materials are derived for general polyhedral meshes. The scheme is general in two and three dimensions. Since the scheme to be developed in this paper is intended for multi-physics code with adaptive mesh refinement (AMR), we present the scheme on mesh generated from AMR. The correctness and features of the scheme are demonstrated for transient problems and steady states in one-, two-, and three-dimensional simulations for heat conduction and radiation heat transfer. The test problems involve dramatically different materials.

Although the Discontinuous Galerkin (dg) method has seen widespread use for compressible flow problems in a single fluid with constant material properties, it has yet to be implemented in a consistent fashion for compressible multiphase flows with shocks and interfaces. Specifically, it is challenging to design a scheme that meets the following requirements: conservation, high-order accuracy in smooth regions and non-oscillatory behavior at discontinuities (in particular, material interfaces). Following the interface-capturing approach of Abgrall 1, we model flows of multiple fluid components or phases using a single equation of state with variable material properties; discontinuities in these properties correspond to interfaces. To represent compressible phenomena in solids, liquids, and gases, we present our analysis for equations of state belonging to the Mie-Gruneisen family. Within the dg framework, we propose a conservative, high-order accurate, and non-oscillatory limiting procedure, verified with simple multifluid and multiphase problems. We show analytically that two key elements are required to prevent spurious pressure oscillations at interfaces and maintain conservation: (i) the transport equation(s) describing the material properties must be solved in a non-conservative weak form, and (ii) the suitable variables must be limited (density, momentum, pressure, and appropriate properties entering the equation of state), coupled with a consistent reconstruction of the energy. Further, we introduce a physics-based discontinuity sensor to apply limiting in a solution-adaptive fashion. We verify this approach with one- and two-dimensional problems with shocks and interfaces, including high pressure and density ratios, for fluids obeying different equations of state to illustrate the robustness and versatility of the method. The algorithm is implemented on parallel graphics processing units (gpu) to achieve high speedup.

When applying the method of manufactured solutions (MMS) on computational fluid dynamic software, determining the exact solutions and source terms for finite volume codes where the stored value is an integrated average over the control volume is non-trivial and not frequently discussed. MMS with discontinuities further complicates the problem of determining these values. In an effort to adapt the standard MMS procedure to solutions that contain discontinuities we show that Newton-Cotes and Gauss quadrature numerical integration methods exhibit high error, first order limitations. We propose a new method for determining the exact solutions and source terms on a uniform structured grid containing shock discontinuities by performing linearly and quadratically exact transformations on split cells. Transformations are performed on triangular and quadrilateral elements of a systematically divided discontinuous cell. Using a quadratic transformation in conjunction with a nine point Gauss quadrature method, a minimum of 4th order accuracy is achieved for fully general solutions and shock shapes. A linear approximation of curved shocks is also experimentally shown to be 2nd order accurate. The numerical integration method is then applied to a CFD code using simple discontinuous manufactured solutions which return consistent 1st order convergence values. The result is an important step towards being able to use MMS to verify solutions with discontinuities. This work also highlights the use of higher order numerical integration techniques for continuous and discontinuous solutions that are required for MMS on higher order finite volume codes.

We describe a novel interface-sharpening approach for efficient numerical resolution of a compressible homogeneous two-phase flow governed by a quasi-conservative five-equation model of Allaire et al. (2001) [1] . The algorithm uses a semi-discrete wave propagation method to find approximate solution of this model numerically. In the algorithm, in regions near the interfaces where two different fluid components are present within a cell, the THINC (Tangent of Hyperbola for INterface Capturing) scheme is used as a basis for the reconstruction of a sub-grid discontinuity of volume fractions at each cell edge, and it is complemented by a homogeneous-equilibrium-consistent technique that is derived to ensure a consistent modeling of the other interpolated physical variables in the model. In regions away from the interfaces where the flow is single phase, standard reconstruction scheme such as MUSCL or WENO can be used for obtaining high-order interpolated states. These reconstructions are then used as the initial data for Riemann problems, and the resulting fluctuations form the basis for the spatial discretization. Time integration of the algorithm is done by employing a strong stability-preserving Runge–Kutta method. Numerical results are shown for sample problems with the Mie–Gruneisen equation of state for characterizing the materials of interests in both one and two space dimensions that demonstrate the feasibility of the proposed method for interface-sharpening of compressible two-phase flow. To demonstrate the competitiveness of our approach, we have also included results obtained using the anti-diffusion interface sharpening method.

The generalized Riemann problems (GRP) for nonlinear hyperbolic systems of balance laws in one space dimension are now well-known and can be formulated as follows: Given initial data which are smooth on two sides of a discontinuity, determine the time evolution of the solution near the discontinuity. In particular, the GRP of (k+1)th order high-resolution is based on an analytical evaluation of time derivatives up to kth order, which turns out to be dependent only on the spatial derivatives up to kth order. While the classical Riemann problem serves as a primary ''building block'' in the construction of many numerical schemes (most notably the Godunov scheme), the analytic study of GRP will lead to an array of ''GRP schemes'', which extend the Godunov scheme. Currently there are extensive studies on the second-order GRP scheme, which proves to be robust and is capable of resolving complex multidimensional fluid dynamic problems (Ben-Artzi and Falcovitz, 2003 [4]). In this paper, we provide a new approach for solving the GRP for the compressible flow system towards high order accuracy. The derivation of second-order GRP solver is more concise compared to those in previous works and the third-order GRP (or quadratic GRP) is resolved for the first time. The latter is shown to be necessary through numerical experiments with strong discontinuities. Our method relies heavily on the new treatment of rarefaction waves. Indeed, as a main technical step, the ''propagation of singularities'' argument for the rarefaction fan, is simplified by deriving the linear ODE systems for the ''evolution'' of ''characteristic derivatives'', in the x-t coordinates, for generalized Riemann invariants. The case of sonic point is incorporated into a general treatment. The accuracy of the derived GRP solvers is justified and numerical examples are presented for the performance of the resulting schemes.

A new reconstruction algorithm is developed to obtain diffusion schemes with cell-centered unknowns only. The main characteristic of the new algorithm is the flexibility of stencils when the auxiliary unknowns are reconstructed with cell-centered unknowns. The stencils are selected depending on the mesh geometry and discontinuities of diffusion coefficients. Moreover, an explicit expression is derived for interpolating the auxiliary unknowns in terms of cell-centered unknowns, and the auxiliary unknowns can be defined at any point on the edge. The algorithm is applied to construct several new diffusion schemes, whose effectiveness is illustrated by numerical experiments. For anisotropic problems with or without discontinuities, nearly second order accuracy is achieved on skewed meshes.

International audience; This article describes a discontinuous implementation of residual distribution for shallow-water flows. The emphasis is put on the space-time implementation of residual distribution for the time-dependent system of equations with discontinuity in time only. This lifts the time-step restriction that even implicit continuous residual distribution schemes invariably suffer from, and thus leads to an unconditionally stable discretisation. The distributions are the space-time variants of the upwind distributions for the steady-state system of equations and are designed to satisfy the most important properties of the original mathematical equations: positivity, linearity preservation, conservation and hydrostatic balance. The purpose of the several numerical examples presented in this article is twofold. First, to show that the discontinuous numerical discretisation does indeed exhibit all the desired properties when applied to the shallow-water equations. Second, to investigate how much the time step can be increased without adversely affecting the accuracy of the scheme and whether this translates into gains in computational efficiency. Comparison to other existing residual distribution schemes is also provided to demonstrate the improved performance of the scheme

Subcell resolution has been used in the Finite Volume Method (FVM) to obtain accurate approximations of discontinuities in the physical space. Stochastic methods are usually based on local adaptivity for resolving discontinuities in the stochastic dimensions. However, the adaptive refinement in the probability space is ineffective in the non-intrusive uncertainty quantification framework, if the stochastic discontinuity is caused by a discontinuity in the physical space with a random location. The dependence of the discontinuity location in the probability space on the spatial coordinates then results in a staircase approximation of the statistics, which leads to first-order error convergence and an underprediction of the maximum standard deviation. To avoid these problems, we introduce subcell resolution into the Simplex Stochastic Collocation (SSC) method for obtaining a truly discontinuous representation of random spatial discontinuities in the interior of the cells discretizing the probability space. The presented SSC–SR method is based on resolving the discontinuity location in the probability space explicitly as function of the spatial coordinates and extending the stochastic response surface approximations up to the predicted discontinuity location. The applications to a linear advection problem, the inviscid Burgers’ equation, a shock tube problem, and the transonic flow over the RAE 2822 airfoil show that SSC–SR resolves random spatial discontinuities with multiple stochastic and spatial dimensions accurately using a minimal number of samples.

We present a moving mesh method suitable for solving two-dimensional and axisymmetric three-liquid flows with triple junction points. This method employs a body-fitted unstructured mesh where the interfaces between liquids are lines of the mesh system, and the triple junction points (if exist) are mesh nodes. To enhance the accuracy and the efficiency of the method, the mesh is constantly adapted to the evolution of the interfaces by refining and coarsening the mesh locally; dynamic boundary conditions on interfaces, in particular the triple points, are therefore incorporated naturally and accurately in a finite-element formulation. In order to allow pressure discontinuity across interfaces, double-values of pressure are necessary for interface nodes and triple-values of pressure on triple junction points. The resulting non-linear system of mass and momentum conservation is then solved by an Uzawa method, with the zero resultant condition on triple points reinforced at each time step. The method is used to investigate the rising of a liquid drop with an attached bubble in a lighter liquid.

For the numerical simulation of detonation of condensed phase explosives, a complex equation of state (EOS), such as the Jones–Wilkins–Lee (JWL) EOS or the Cochran–Chan (C–C) EOS, are widely used. However, when a conservative scheme is used for solving the Euler equations with such equations of state, a spurious solution across the contact discontinuity, a well known phenomenon in multi-fluid systems, arises even for single materials. In this work, we develop a generalised Osher-type scheme in an adaptive primitive–conservative framework to overcome the aforementioned difficulties. Resulting numerical solutions are compared with the exact solutions and with the numerical solutions from the Godunov method in conjunction with the exact Riemann solver for the Euler equations with Mie–Gruneisen form of equations of state, such as the JWL and the C–C equations of state. The adaptive scheme is extended to second order and its empirical convergence rates are presented, verifying second order accuracy for smooth solutions. Through a suite of several tests problems in one and two space dimensions we illustrate the failure of conservative schemes and the capability of the methods of this paper to overcome the difficulties.

This is the third paper in the series of our conservative front-tracking method. In this paper, we describe how our method tracks fluid interfaces in multi-fluid flows. Two important ingredients in our conservative front-tracking method in tracking fluid interfaces are: (1) the velocities and pressures of the left and right cell-averages in a discontinuity cell are respectively maintained to be equal to each other, and in doing so the physics that the normal velocity and pressure are continuous cross the interface is simulated but that the tangential velocity may be discontinuous is ignored, and (2) a so-called numerical surface dissipation is designed on the tracked interface to eliminate possible numerical instability there, and we believe that this numerical surface dissipation is a good substitute for the missing physical dissipation acting on the interface. We then present numerical simulation of Haas-Sturtevant's two shock-bubble interaction experiments [J.F. Haas, B. Sturtevant, Interaction of weak shock waves with cylindrical and spherical gas inhomogeneities, J. Fluid Mech. 181 (1987) 41-76] using this method. Our numerical results are in good agreement with the experimental and other numerical results in the early times of the flow. Moreover, our numerical results are also in good agreement with the experimental results in the later times of the flow and give clear pictures of the bubble deformation then, which show that the right boundaries of the bubble behave just as Rychtmyer-Meshkov instabilities with the shock coming from either heavy or light gases.

The goal of this paper is to relate numerical dissipations that are inherited in high order shock-capturing schemes with the onset of wrong propagation speed of discontinuities. For pointwise evaluation of the source term, previous studies indicated that the phenomenon of wrong propagation speed of discontinuities is connected with the smearing of the discontinuity caused by the discretization of the advection term. The present study focuses only on solving the reactive system by the fractional step method using the Strang splitting. Studies shows that the degree of wrong propagation speed of discontinuities is highly dependent on the accuracy of the numerical method. The manner in which the smearing of discontinuities is contained by the numerical method and the overall amount of numerical dissipation being employed play major roles. Depending on the numerical method, time step and grid spacing, the numerical simulation may lead to (a) the correct solution (within the truncation error of the scheme), (b) a divergent solution, (c) a wrong propagation speed of discontinuities solution or (d) other spurious solutions that are solutions of the discretized counterparts but are not solutions of the governing equations. The findings might shed some light on the reported difficulties in numerical combustion and problems with stiff nonlinear (homogeneous) source terms and discontinuities in general.

The evolution of high-speed initially laminar multicomponent flows into a turbulent multi-material mixing entity, e.g., in the Richtmyer-Meshkov instability, poses significant challenges for high-fidelity numerical simulations. Although high-order shock- and interface-capturing schemes represent such flows well at early times, the excessive numerical dissipation thereby introduced and the resulting computational cost prevent the resolution of small-scale features. Furthermore, unless special care is taken, shock-capturing schemes generate spurious pressure oscillations at material interfaces where the specific heats ratio varies. To remedy these problems, a solution-adaptive high-order central/shock-capturing finite difference scheme is presented for efficient computations of compressible multi-material flows, including turbulence. A new discontinuity sensor discriminates between smooth and discontinuous regions. The appropriate split form of (energy preserving) central schemes is derived for flows of smoothly varying specific heats ratio, such that spurious pressure oscillations are prevented. High-order accurate weighted essentially non-oscillatory (WENO) schemes are applied only at discontinuities; the standard approach is followed for shocks and contacts, but material discontinuities are treated by interpolating the primitive variables. The hybrid nature of the method allows for efficient and accurate computations of shocks and broadband motions, and is shown to prevent pressure oscillations for varying specific heats ratios. The method is assessed through one-dimensional problems with shocks, sharp interfaces and smooth distributions of specific heats ratio, and the two-dimensional single-mode inviscid and viscous Richtmyer-Meshkov instability with re-shock.

Multi-element uncertainty quantification approaches can robustly resolve the high sensitivities caused by discontinuities in parametric space by reducing the polynomial degree locally to a piecewise linear approximation. It is important to extend the higher degree interpolation in the smooth regions up to a thin layer of linear elements that contain the discontinuity to maintain a highly accurate solution. This is achieved here by introducing Essentially Non-Oscillatory (ENO) type stencil selection into the Simplex Stochastic Collocation (SSC) method. For each simplex in the discretization of the parametric space, the stencil with the highest polynomial degree is selected from the set of candidate stencils to construct the local response surface approximation. The application of the resulting SSC-ENO method to a discontinuous test function shows a sharper resolution of the jumps and a higher order approximation of the percentiles near the singularity. SSC-ENO is also applied to a chemical model problem and a shock tube problem to study the impact of uncertainty both on the formation of discontinuities in time and on the location of discontinuities in space.

Up to an arbitrary order of the Chapman-Enskog expansion of the kinetic Bhatnagar-Gross-Krook (BGK) equation, a corresponding analytic solution can be obtained. Based on such a compact exact solution, a new gas-kinetic scheme is constructed for the compressible Navier-Stokes equations. Instead of using a discontinuous initial condition in the gas-kinetic BGK-NS method Xu [K. Xu, A gas-kinetic BGK scheme for the Navier-Stokes equations, and its connection with artificial dissipation and Godunov method, J. Comput. Phys. 171 (2001) 289-335.], the new scheme starts with a continuous initial flow distribution at a cell interface which is obtained through an upwind-biased WENO reconstruction, and uses the time accurate solution for the flux evaluation. The new kinetic scheme not only preserves favorable properties of the existing BGK-NS method, such as stability, high resolution, and good performance in capturing discontinuity, but also achieves a very high efficiency, which is even more efficient than the same order well-defined classical finite difference scheme based on the macroscopic governing equations directly. The stability, accuracy, and efficiency of the new scheme are evaluated quantitatively through numerical tests. The new scheme captures sharp discontinuity without post-shock oscillation, and has high accuracy in resolving viscous solution. Due to the use of time-accurate analytical solution, the overall performance of the new scheme is superior in comparison with the finite difference or finite volume schemes which start from the same initial reconstruction and use the numerical Runge-Kutta methods for time accuracy.

A fully-implicit ghost-cell immersed boundary method for simulations of flow over complex moving bodies on a Cartesian grid is presented. The present immersed boundary method is highly capable of controlling the generation of spurious force oscillations on the surface of a moving body, thereby producing an accurate and stable solution. Spurious force oscillations on the surface of an immersed moving body are reduced by alleviating spatial and temporal discontinuities in the pressure and velocity fields across non-grid conforming immersed boundaries. A sharp-interface ghost-cell immersed-boundary method is coupled with a mass source and sink algorithm to improve the conservation of mass across non-grid conforming immersed boundaries. To facilitate the control for the temporal discontinuity in the flow field due to a motion of an immersed body, a fully-implicit time-integration scheme is employed. A novel backward time-integration scheme is developed to effectively treat multiple layers of fresh cells generated by a motion of an immersed body at a high CFL number condition. The present backward time-integration scheme allows to impose more accurate and stable velocity vectors on fresh cells than those interpolated. The effectiveness of the present fully-implicit ghost-cell immersed boundary method coupled with a mass source and sink algorithm for reducing spurious force oscillations during simulations of moving body problems is demonstrated in a number of test cases.

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.

The effect of discontinuity in the state variables on optimization problems is investigated on the quasi-one-dimensional Euler equations in the discrete level. A pressure minimization problem and a pressure matching problem are considered. We find that the objective functional can be smooth in the continuous level and yet be non-smooth in the discrete level as a result of the shock crossing grid points. Higher resolution can exacerbate that effect making grid refinement counter productive for the purpose of computing the discrete sensitivities. First and second order sensitivities, as well as the adjoint solution, are computed exactly at the shock and its vicinity and are compared to the continuous solution. It is shown that in the discrete level the first order sensitivities contain a spike at the shock location that converges to a delta function with grid refinement, consistent with the continuous analysis. The numerical Hessian is computed and its consistency with the analytical Hessian is discussed for different flow conditions. It is demonstrated that consistency is not guaranteed for shocked flows. We also study the different terms composing the Hessian and propose some stable approximation to the continuous Hessian.

Prosperetti's seminal Physalis method, an Immersed Boundary/spectral method, had been used extensively to investigate fluid flows with suspended solid particles. Its underlying idea of creating a cage and using a spectral general analytical solution around a discontinuity in a surrounding field as a computational mechanism to enable the accommodation of physical and geometric discontinuities is a general concept, and can be applied to other problems of importance to physics, mechanics, and chemistry. In this paper we provide a foundation for the application of this approach to the determination of the distribution of electric charge in heterogeneous mixtures of dielectrics and conductors. The proposed Physalis method is remarkably accurate and efficient. In the method, a spectral analytical solution is used to tackle the discontinuity and thus the discontinuous boundary conditions at the interface of two media are satisfied exactly. Owing to the hybrid finite difference and spectral schemes, the method is spectrally accurate if the modes are not sufficiently resolved, while higher than second-order accurate if the modes are sufficiently resolved, for the solved potential field. Because of the features of the analytical solutions, the derivative quantities of importance, such as electric field, charge distribution, and force, have the same order of accuracy as the solved potential field during postprocessing. This is an important advantage of the Physalis method over other numerical methods involving interpolation, differentiation, and integration during postprocessing, which may significantly degrade the accuracy of the derivative quantities of importance. The analytical solutions enable the user to use relatively few mesh points to accurately represent the regions of discontinuity. In addition, the spectral convergence and a linear relationship between the cost of computer memory/computation and particle numbers results in a very efficient method. In the present paper, the accuracy of the method is numerically investigated by example computations using one dielectric particle, one isolated conductor particle, one conductor particle connected to an external source with imposed voltage, and two conductor/dielectric particles with strong interactions. The efficiency of the method is demonstrated with one million particles, which suggests that the method can be used for many important engineering applications of broad interest.

Meshless methods based on compact radial basis functions (RBFs) are proposed for modelling photonic crystals (PhCs). When modelling two-dimensional PhCs two generalised eigenvalue problems are formed, one for the transverse-electric (TE) mode and the other for the transverse-magnetic (TM) mode. Conventionally, the Band Diagrams for two-dimensional PhCs are calculated by either the plane wave expansion method (PWEM) or the finite element method (FEM). Here, the eigenvalue equations for the two-dimensional PhCs are solved using RBFs based meshless methods. For the TM mode a meshless local strong form method (RBF collocation) is used, while for the tricker TE mode a meshless local weak form method (RBF Galerkin) is used (so that the discontinuity of the dielectric function @e(x) can naturally be modelled). The results obtained from the meshless methods are found to be in good agreement with the standard PWEM. Thus, the meshless methods are proved to be a promising scheme for predicting photonic band gaps.

In this paper, a class of finite difference schemes which achieves low dispersion and controllable dissipation in smooth region and robust shock-capturing capabilities in the vicinity of discontinuities is presented. Firstly, a sufficient condition for semi-discrete finite difference schemes to have independent dispersion and dissipation is derived. This condition enables a novel approach to separately optimize the dissipation and dispersion properties of finite difference schemes and a class of schemes with minimized dispersion and controllable dissipation is thus obtained. Secondly, for the purpose of shock-capturing, one of these schemes is used as the linear part of the WENO scheme with symmetrical stencils to constructed an improved WENO scheme. At last, the improved WENO scheme is blended with its linear counterpart to form a new hybrid scheme for practical applications. The proposed scheme is accurate, flexible and robust. The accuracy and resolution of the proposed scheme are tested by the solutions of several benchmark test cases. The performance of this scheme is further demonstrated by its application in the direct numerical simulation of compressible turbulent channel flow between isothermal walls.

In immiscible two-phase flows, jumps or kinks are present in the velocity and pressure fields across the interfaces of the two fluids. The extended finite element method (XFEM) is able to reproduce such discontinuities within elements. Robust and accurate interface capturing schemes with no restrictions on the interface topology are thereby enabled. This paper investigates different enrichment schemes and time-integration schemes within the XFEM. Test cases with and without surface tension on moving or stationary meshes are studied and compared to interface tracking results when possible. A particularly useful setting is extracted which is recommended for two-phase flows. An extension of this formulation for the simulation of free-surface flows and of floating objects is proposed.

When a discrete-forcing immersed boundary method is applied to moving-body problems, it produces spurious force oscillations on a solid body. In the present study, we identify two sources of these force oscillations. One source is from the spatial discontinuity in the pressure across the immersed boundary when a grid point located inside a solid body becomes that of fluid with a body motion. The addition of mass source/sink together with momentum forcing proposed by Kim et al. [J. Kim, D. Kim, H. Choi, An immersed-boundary finite volume method for simulations of flow in complex geometries, Journal of Computational Physics 171 (2001) 132-150] reduces the spurious force oscillations by alleviating this pressure discontinuity. The other source is from the temporal discontinuity in the velocity at the grid points where fluid becomes solid with a body motion. The magnitude of velocity discontinuity decreases with decreasing the grid spacing near the immersed boundary. Four moving-body problems are simulated by varying the grid spacing at a fixed computational time step and at a constant CFL number, respectively. It is found that the spurious force oscillations decrease with decreasing the grid spacing and increasing the computational time step size, but they depend more on the grid spacing than on the computational time step size.