Spatial Discretization Schemes Research Articles

Viscous Regularization for the Non-equilibrium Seven-Equation Two-Phase Flow Model

In this paper, a viscous regularization is derived for the non-equilibrium seven-equation two-phase flow model (SEM). This regularization, based on an entropy condition, is an artificial viscosity stabilization technique that selects a weak solution satisfying an entropy-minimum principle. The viscous regularization ensures nonnegativity of the entropy residual, is consistent with the viscous regularization for Euler equations when one phase disappears, does not depend on the spatial discretization scheme chosen, and is compatible with the generalized Harten entropies. We investigate the behavior of the proposed viscous regularization for two important limit-cases. First, a Chapman---Enskog expansion is performed for the regularized SEM and we show that the five-equation flow model of Kapila is recovered with a well-scaled viscous regularization. Second, a low-Mach asymptotic limit of the regularized seven-equation flow model is carried out whereby the scaling of the non-dimensional numbers associated with the viscous terms is determined such that an incompressible two-phase flow model, with a properly scaled regularization, is recovered. Both limit-cases are illustrated with one-dimensional numerical results, including two-phase flow shock tube tests and steady-state two-phase flows in converging-diverging nozzles. A continuous finite element discretization is employed for all numerical simulations.

Improving a conjugate-heat-transfer immersed-boundary method

Purpose– The purpose of this paper is to provide the current state of the art in the development of a computer code combining an immersed boundary method with a conjugate heat transfer (CHT) approach, including some new findings. In particular, various treatments of the fluid-solid-interface conditions are compared in order to determine the most accurate one. Most importantly, the method is capable of computing a challenging three dimensional compressible turbulent flow past an air cooled turbine vane.Design/methodology/approach– The unsteady Reynolds-averaged Navier–Stokes (URANS) equations are solved within the fluid domain, whereas the heat conduction equation is solved within the solid one, using the same spatial discretization and time-marching scheme. At the interface boundary, the temperatures and heat fluxes within the fluid and the solid are set to be equal using three different approximations.Findings– This work provides an accurate and efficient code for solving three dimensional CHT problems, such as the flow through an air cooled gas turbine cascade, using a coupled immersed boundary (IB) CHT methodology. A one-to-one comparison of three different interface-condition approximations has shown that the two multidimensional ones are slightly superior to the early treatment based on a single direction and that the one based on a least square reconstruction of the solution near the IB minimizes the oscillations caused by the Cartesian grid. This last reconstruction is then used to compute a compressible turbulent flow of industrial interest, namely, that through an air cooled gas turbine cascade. Another interesting finding is that the very promising approach based on wall functions does not combine favourably with the interface conditions for the temperature and the heat flux. Therefore, current and future work aims at developing and testing appropriate temperature wall functions, in order to further improve the accuracy – for a given grid – or the efficiency – for a given accuracy – of the proposed methodology.Originality/value– An accurate and efficient IB CHT method, using a state of the art URANS parallel solver, has been developed and tested. In particular, a detailed study has elucidated the influence of different interface treatments of the fluid-solid boundary upon the accuracy of the computations. Last but not least, the method has been applied with success to solve the well-known CHT problem of compressible turbulent flow past the C3X turbine guide vane.

Application of Jacobian-free Newton–Krylov method in implicitly solving two-fluid six-equation two-phase flow problems: Implementation, validation and benchmark

This work represents a first-of-its-kind successful application to employ advanced numerical methods in solving realistic two-phase flow problems with two-fluid six-equation two-phase flow model. These advanced numerical methods include high-resolution spatial discretization scheme with staggered grids (high-order) fully implicit time integration schemes, and Jacobian-free Newton–Krylov (JFNK) method as the nonlinear solver. The computer code developed in this work has been extensively validated with existing experimental flow boiling data in vertical pipes and rod bundles, which cover wide ranges of experimental conditions, such as pressure, inlet mass flux, wall heat flux and exit void fraction. Additional code-to-code benchmark with the RELAP5-3D code further verifies the correct code implementation. The combined methods employed in this work exhibit strong robustness in solving two-phase flow problems even when phase appearance (boiling) and realistic discrete flow regimes are considered. Transitional flow regimes used in existing system analysis codes, normally introduced to overcome numerical difficulty, were completely removed in this work. This in turn provides the possibility to utilize more sophisticated flow regime maps in the future to further improve simulation accuracy.

Solving phase appearance/disappearance two-phase flow problems with high resolution staggered grid and fully implicit schemes by the Jacobian-free Newton–Krylov Method

The phase appearance/disappearance issue presents serious numerical challenges in two-phase flow simulations. Although many existing nuclear reactor safety analysis codes use different kinds of treatments for the phase appearance/disappearance problem, there are no fully consensual solutions. In this paper, a high-resolution spatial discretization scheme on staggered grid mesh and fully implicit time integration methods was used in two-phase flow simulations. Numerical treatments of the momentum equations were introduced for the disappearing phase without the need to use a cut-off value on void fraction. A discretized nonlinear system was solved using the Jacobian-free Newton Krylov (JFNK) method, which does not require expensive and error-prone derivation and implementation of analytical Jacobian matrix. These methods were tested with several two-phase flow problems with phase appearance/disappearance phenomena considered, such as a linear advection problem, an oscillating manometer problem, and a sedimentation problem. The combined methods demonstrated extremely robust and stable behaviors in solving the one-dimensional two-phase flow problems with phase appearance/disappearance. High-resolution spatial discretization and second-order fully implicit method also demonstrated their capabilities in significantly reducing numerical errors.

Dispersion analysis of Robert–Asselin type filters for linear non-dispersive and dispersive systems

The present paper deals with the dispersion relation preserving (DRP) analysis of the Robert–Asselin type filters coupled with leapfrog time integration method for time dependent non-dispersive and dispersive model systems. As, leapfrog time advancement method is widely used in the numerical modeling of atmosphere and ocean dynamics, despite the fact that in addition to physical mode it also admits a spurious mode. This is the major disadvantage of leapfrog time integration scheme in discrete computing. To suppress the spurious mode associated with leapfrog time integration method, Robert–Asselin type filters are extensively used in the literature. Here, the dispersion analysis is based on the spectral analysis by considering filtered/unfiltered leapfrog time integration method along with spatial discretization schemes. Spatial discretization is performed using second order centered difference and optimized compact schemes. Furthermore, to assess the efficacy of these methods one dimensional convection equation is used as a model equation for the non-dispersive case, while linearized rotating shallow water equations (LRSWE), as a model test problem for the dispersive case.

Momentum conservation in Multi-Level Multi-Domain (MLMD) simulations

Momentum conservation and self-forces reduction are challenges for all Particle-In-Cell (PIC) codes using spatial discretization schemes which do not fulfill the requirement of translational invariance of the grid Green's function. We comment here on the topic applied to the recently developed Multi-Level Multi-Domain (MLMD) method. The MLMD is a semi-implicit method for PIC plasma simulations. The multi-scale nature of plasma processes is addressed by using grids with different spatial resolutions in different parts of the domain.

Impact of rotating and fixed nozzles on vortex breakdown in compressible swirling jet flows

Vortex breakdown of swirling, round jet flows is investigated in the compressible, subsonic regime by means of Direct Numerical Simulation (DNS). This is achieved by solving the compressible Navier–Stokes equations on a cylindrical grid using high-order spatial and temporal discretization schemes. The Reynolds number is Re=ρc∘wc∘R∘/μc∘=5000 and the flow is moderately compressible with Mach number Ma=wc∘/γRair∘Tc∘=0.6. The integral swirl number at the inflow is Sint=0.85. The parameters are chosen properly so as to make comparisons with existing experiments at lower Mach numbers possible while still enabling a study of compressible and baroclinic effects. Different from previous numerical investigations, a nozzle immersed in the fluid is included in the computational domain and is modelled as an isothermal no-slip wall, either rotating with the mean azimuthal flow direction or kept at rest. The present investigation aims to clarify the role played by the nozzle wall motion for the vortex breakdown of the swirling jet. We study the nozzle flow as well as the swirling jet flow simultaneously, a novelty for numerical investigations of vortex breakdown in swirling jets. Depending on the nozzle wall motion, the flow differs significantly upstream of the vortex breakdown: for the rotating nozzle, the flow inside the nozzle is purely laminar and the azimuthal boundary layer at the outer nozzle wall gives rise to the axisymmetric mode n=0 and a single-helix type instability with azimuthal wave number n=1. With the nozzle at rest, a transitional flow is observed within the nozzle where a helical instability with azimuthal wave number n=12 dominates, growing in the boundary layer at the nozzle wall. For both nozzle setups, the helical instabilities observed for the nozzle flow interact with the developing vortex breakdown and the conical shear-layer downstream of the nozzle. For the nozzle at rest, this interaction results in a vortex breakdown configuration which is shifted in the upstream direction and which has a smaller radial and streamwise extent compared to the rotating nozzle case and the recirculation intensity is higher. The dominant frequency is highly influenced by the flow upstream of the vortex breakdown and is substantially higher for the nozzle at rest. Although the nozzle flow field differs for the two configurations and therefore alters the vortex breakdown downstream, a single-helix type instability n=1 governing the vortex breakdown is found for both cases. This provides strong evidence for the robustness of the instability mechanisms leading to vortex breakdown.

Implicitly solving phase appearance and disappearance problems using two-fluid six-equation model

Phase appearance and disappearance issue presents serious numerical challenges in two-phase flow simulations using the two-fluid six-equation model. Numerical challenges arise from the singular equation system when one phase is absent, as well as from the discontinuity in the solution space when one phase appears or disappears. In this work, a high-resolution spatial discretization scheme on staggered grids and fully implicit methods were applied for the simulation of two-phase flow problems using the two-fluid six-equation model. A Jacobian-free Newton-Krylov (JFNK) method was used to solve the discretized nonlinear problem. An improved numerical treatment was proposed and proved to be effective to handle the numerical challenges. The treatment scheme is conceptually simple, easy to implement, and does not require explicit truncations on solutions, which is essential to conserve mass and energy. Various types of phase appearance and disappearance problems relevant to thermal-hydraulics analysis have been investigated, including a sedimentation problem, an oscillating manometer problem, a non-condensable gas injection problem, a single-phase flow with heat addition problem and a subcooled flow boiling problem. Successful simulations of these problems demonstrate the capability and robustness of the proposed numerical methods and numerical treatments. Volume fraction of the absent phase can be calculated effectively as zero.

Well-Balanced, Conservative Finite Difference Algorithm for Atmospheric Flows

The numerical simulation of meso-, convective-, and microscale atmospheric flows requires the solution of the Euler or the Navier–Stokes equations. Nonhydrostatic weather prediction algorithms often solve the equations in terms of derived quantities such as Exner pressure and potential temperature (and are thus not conservative) and/or as perturbations to the hydrostatically balanced equilibrium state. This paper presents a well-balanced, conservative finite difference formulation for the Euler equations with a gravitational source term, where the governing equations are solved as conservation laws for mass, momentum, and energy. Preservation of the hydrostatic balance to machine precision by the discretized equations is essential because atmospheric phenomena are often small perturbations to this balance. The proposed algorithm uses the weighted essentially nonoscillatory and compact-reconstruction weighted essentially nonoscillatory schemes for spatial discretization that yields high-order accurate solutions for smooth flows and is essentially nonoscillatory across strong gradients; however, the well-balanced formulation may be used with other conservative finite difference methods. The performance of the algorithm is demonstrated on test problems as well as benchmark atmospheric flow problems, and the results are verified with those in the literature.

A Comparative Study of Two Spatial Discretization Schemes for Advection Equation

In this paper, we describe a comparison of two spatial discretization schemes for the advection equation, namely the first finite difference method and the method of lines. The stability of the methods has been studied by Von Neumann method and with the matrix analysis. The methods are applied to a number of test problems to compare the accuracy and computational efficiency. We show that both discretization techniques approximate correctly solution of advection equation and compare their accuracy and performance.

High Order Implicit-explicit General Linear Methods with Optimized Stability Regions

In the numerical solution of partial differential equations using a method-of-lines approach, the availability of high order spatial discretization schemes motivates the development of sophisticated high order time integration methods. For multiphysics problems with both stiff and nonstiff terms implicit-explicit (IMEX) time stepping methods attempt to combine the lower cost advantage of explicit schemes with the favorable stability properties of implicit schemes. Existing high order IMEX Runge--Kutta or linear multistep methods, however, suffer from accuracy or stability limitations. This work shows that IMEX general linear methods (GLMs) are competitive alternatives to classic IMEX schemes for large problems arising in practice. High order IMEX-GLMs are constructed in the partitioned GLM framework developed earlier by the authors [J. Sci. Comput., 61 (2014), pp. 119--144]. The stability regions of the new schemes are optimized numerically. The resulting IMEX-GLMs have similar stability properties as IMEX Runge--Kutta methods, but they do not suffer from order reduction and are superior in terms of accuracy and efficiency. The new IMEX-GLMs have considerably better stability properties than the IMEX linear multistep methods. Numerical experiments with two- and three-dimensional test problems illustrate the potential of the new schemes to speed up complex applications.

Numerical modelling of carbon dioxide leakage from an offshore pipeline

A three-dimensional transient computational model was developed to simulate the break up of carbon dioxide escaping from a pipeline leak to the seawater. Important physical phenomena such as turbulence and mass transfer during evaporation of the initially liquid CO2 stream are included in the model. Sensitivity of the model for mesh size, turbulence model and the order of the spatial discretisation scheme are shown for different scenarios. A three stage approach is used to achieve sufficient level of convergence for the transient multiphase turbulent flow simulation using the Span and Wagner equation of state and a high-order discretisation schemes. The scenarios considered in this study include a wide range of pipeline operating pressure, ambient seawater depth and leakage size.

Numerical methods for a class of nonlocal diffusion problems with the use of backward SDEs

We propose a novel numerical approach for nonlocal diffusion equations Du et al. (2012) with integrable kernels, based on the relationship between the backward Kolmogorov equation and backward stochastic differential equations (BSDEs) driven by Lèvy jumps processes. The nonlocal diffusion problem under consideration is converted to a BSDE, for which numerical schemes are developed. As a stochastic approach, the proposed method completely avoids the challenge of iteratively solving non-sparse linear systems, arising from the nature of nonlocality. This allows for embarrassingly parallel implementation and also enables adaptive approximation techniques to be incorporated in a straightforward fashion. Moreover, our method recovers the convergence rates of classic deterministic approaches (e.g. finite element or collocation methods), due to the use of high-order temporal and spatial discretization schemes. In addition, our approach can handle a broad class of problems with general inhomogeneous forcing terms as long as they are globally Lipschitz continuous. Rigorous error analysis of the new method is provided as several numerical examples that illustrate the effectiveness and efficiency of the proposed approach.

Extension of the Henry semi-analytical solution for saltwater intrusion in stratified domains

The Henry semi-analytical solution is developed for stratified aquifers with exponential and Gaussian permeability–depth relationships and small dispersion. The semi-analytical solution is sought by expanding the stream function and the concentration in an infinite Fourier series truncated at given orders. Due to the heterogeneity, the semi-analytical solution contains additional terms and quickly becomes impractical because of the high truncation orders required to yield stable results. An appropriate evaluation of the most expensive term, involving fairly complex summations, is proposed to render the computation of the semi-analytical solution affordable. Three configurations with the same transmissivity as that for the homogeneous Henry problem are investigated. The first configuration considers an exponential decay of the permeability (k) with depth. The second and third configurations consider an exponential increase and a Gaussian high-low-high distribution of k with depth, respectively. The three configurations are also investigated numerically using an accurate numerical code based on the method of lines and advanced spatial discretization schemes. An excellent agreement is obtained between the semi-analytical and the numerical solutions for all test cases. The results show that the amount of saltwater intrusion is strongly dependent on the heterogeneity. Further, the heterogeneous Henry problem is found to be more suitable for benchmarking density-driven numerical codes because the numerical solution is more sensitive to the grid size than for the homogeneous problem.

On the Finite Differences Used in Reacting Flow Simulations

AbstractThere exist many spatial discretization schemes that are well able to provide accurate and stable approximations for isothermal turbulent flows. Comparatively little analysis has been made of the performance of these schemes in the presence of temperature gradients driven by combustion. In this paper, the effects of temperature gradients on numerical stability are explored. A surprising result is that temperature gradients in the flow have a tendency to impinge on left half plane (LHP) stability of the spatial discretization scheme. Reasons for this tendency are explored and two remedies are proposed: one based on the particular class of finite difference schemes, and one based on an alternative method of boundary condition specification.

Numerical implementation, verification and validation of two-phase flow four-equation drift flux model with Jacobian-free Newton–Krylov method

This paper presents a numerical investigation on using the Jacobian-free Newton–Krylov (JFNK) method to solve the two-phase flow four-equation drift flux model with realistic constitutive correlations (‘closure models’). The drift flux model is based on Isshi and his collaborators’ work. Additional constitutive correlations for vertical channel flow, such as two-phase flow pressure drop, flow regime map, wall boiling and interfacial heat transfer models, were taken from the RELAP5-3D Code Manual and included to complete the model. The staggered grid finite volume method and fully implicit backward Euler method was used for the spatial discretization and time integration schemes, respectively. The Jacobian-free Newton–Krylov method shows no difficulty in solving the two-phase flow drift flux model with a discrete flow regime map. In addition to the Jacobian-free approach, the preconditioning matrix is obtained by using the default finite differencing method provided in the PETSc package, and consequently the labor-intensive implementation of complex analytical Jacobian matrix is avoided. Extensive and successful numerical verification and validation have been performed to prove the correct implementation of the models and methods. Code-to-code comparison with RELAP5-3D has further demonstrated the successful implementation of the drift flux model.

An upwind vertex centred Finite Volume solver for Lagrangian solid dynamics

A vertex centred Jameson–Schmidt–Turkel (JST) finite volume algorithm was recently introduced by the authors (Aguirre et al., 2014 [1]) in the context of fast solid isothermal dynamics. The spatial discretisation scheme was constructed upon a Lagrangian two-field mixed (linear momentum and the deformation gradient) formulation presented as a system of conservation laws [2–4]. In this paper, the formulation is further enhanced by introducing a novel upwind vertex centred finite volume algorithm with three key novelties. First, a conservation law for the volume map is incorporated into the existing two-field system to extend the range of applications towards the incompressibility limit (Gil et al., 2014 [5]). Second, the use of a linearised Riemann solver and reconstruction limiters is derived for the stabilisation of the scheme together with an efficient edge-based implementation. Third, the treatment of thermo-mechanical processes through a Mie–Grüneisen equation of state is incorporated in the proposed formulation. For completeness, the study of the eigenvalue structure of the resulting system of conservation laws is carried out to demonstrate hyperbolicity and obtain the correct time step bounds for non-isothermal processes. A series of numerical examples are presented in order to assess the robustness of the proposed methodology. The overall scheme shows excellent behaviour in shock and bending dominated nearly incompressible scenarios without spurious pressure oscillations, yielding second order of convergence for both velocities and stresses.

Influence of spatial discretization schemes on accuracy of explicit LES: Canonical problems to engine-like geometries

The choice of the spatial discretization scheme and the subgrid-scale (SGS) model can have a significant impact on the accuracy of Large Eddy Simulations (LES). A systematic study of the influence of the advection term discretization scheme, namely (CDS, LUST) and the SGS model (σ-model, WALE and dynamic Smagorinsky) on the accuracy of the solution is presented in this work. Three canonical cases with increasing complexity are considered in this study, namely, a fully developed turbulent channel flow with Reτ=395, a backward facing step and a wall-mounted hump. Mean errors with respect to DNS or experimental data are quantified in order to compare the relative accuracy of each combination of scheme/model. Detailed comparison of the numerical simulations performed by the open-source CFD code OpenFOAM® shows that the σ-model with the LUST discretization scheme yields the best results. This combination of σ-model with LUST discretization scheme has hence been used to study the turbulent flow characteristics in an engine-like geometry. Results show good agreement with experimental data.

Analysis of triangular C-grid finite volume scheme for shallow water flows

In this paper, a dispersion relation analysis is employed to investigate the finite volume triangular C-grid formulation for two-dimensional shallow-water equations. In addition, two proposed combinations of time-stepping methods with the C-grid spatial discretization are investigated. In the first part of this study, the C-grid spatial discretization scheme is assessed, and in the second part, fully discrete schemes are analyzed. Analysis of the semi-discretized scheme (i.e. only spatial discretization) shows that there is no damping associated with the spatial C-grid scheme, and its phase speed behavior is also acceptable for long and intermediate waves. The analytical dispersion analysis after considering the effect of time discretization shows that the Leap-Frog time stepping technique can improve the phase speed behavior of the numerical method; however it could not damp the shorter decelerated waves. The Adams–Bashforth technique leads to slower propagation of short and intermediate waves and it damps those waves with a slower propagating speed. The numerical solutions of various test problems also conform and are in good agreement with the analytical dispersion analysis. They also indicate that the Adams–Bashforth scheme exhibits faster convergence and more accurate results, respectively, when the spatial and temporal step size decreases. However, the Leap-Frog scheme is more stable with higher CFL numbers.

Applications of high-resolution spatial discretization scheme and Jacobian-free Newton–Krylov method in two-phase flow problems

The majority of the existing reactor system analysis codes were developed using low-order numerical schemes in both space and time. In many nuclear thermal–hydraulics applications, it is desirable to use higher-order numerical schemes to reduce numerical errors. High-resolution spatial discretization schemes provide high order spatial accuracy in smooth regions and capture sharp spatial discontinuity without nonphysical spatial oscillations. In this work, we adapted an existing high-resolution spatial discretization scheme on staggered grids in two-phase flow applications. Fully implicit time integration schemes were also implemented to reduce numerical errors from operator-splitting types of time integration schemes. The resulting nonlinear system has been successfully solved using the Jacobian-free Newton–Krylov (JFNK) method. The high-resolution spatial discretization and high-order fully implicit time integration numerical schemes were tested and numerically verified for several two-phase test problems, including a two-phase advection problem, a two-phase advection with phase appearance/disappearance problem, and the water faucet problem. Numerical results clearly demonstrated the advantages of using such high-resolution spatial and high-order temporal numerical schemes to significantly reduce numerical diffusion and therefore improve accuracy. Our study also demonstrated that the JFNK method is stable and robust in solving two-phase flow problems, even when phase appearance/disappearance exists.