Second-order Accurate Solutions Research Articles

IGG(χ): A new and simple implicit gradient scheme on unstructured meshes

A new and simple implicit Green-Gauss gradient (IGG) scheme for unstructured meshes is proposed exploiting the ideas from the linearity-preserving U-MUSCL scheme to define values at cell faces. We construct an implicit one-parameter family of gradient schemes referred to as IGG(χ) where χ is a free-parameter. The computed gradients are at least first-order accurate on generic polygonal meshes and second-order accurate on uniform Cartesian meshes except when χ=1/3 for which fourth-order accuracy can be realised. A theoretical analysis is carried out to understand the effect of the control parameter χ on accuracy and resolution of the gradients and numerical experiments on various mesh topologies confirm the theoretical findings. Finite volume simulations of the Poisson and Euler equations on Cartesian and unstructured meshes further highlight that the IGG(χ) is a versatile gradient scheme that gives second-order accurate solutions with the iterative convergence of the solver dependent on the choice of the χ parameter. The framework described in this study can also be employed to devise an implicit least-squares gradient scheme that applies equally well to unstructured finite volume and meshfree solvers.

A Volume-of-Fluid method for multicomponent droplet evaporation with Robin boundary conditions

We propose a numerical method tailored to perform interface-resolved simulations of evaporating multicomponent two-phase flows. The novelty of the method lies in the use of Robin boundary conditions to couple the transport equations for the vaporized species in the gas phase and the transport equations of the same species in the liquid phase. The Robin boundary condition is implemented with the cost-effective procedure proposed by Chai et al. [1] and consists of two steps: (1) calculating the normal derivative of the mass fraction fields in cells adjacent to the interface through the reconstruction of a linear polynomial system, and (2) extrapolating the normal derivative and the ghost value in the normal direction using a linear partial differential equation. This methodology yields a second-order accurate solution for the Poisson equation with a Robin boundary condition and a first-order accurate solution for the Stefan problem. The overall methodology is implemented in an efficient two-fluid solver, which includes a Volume-of-Fluid (VoF) approach for the interface representation, a divergence-free extension of the liquid velocity field onto the entire domain to transport the VoF, and the temperature equation to include thermal effects. We demonstrate the convergence of the numerical method to the analytical solution for multicomponent isothermal evaporation and observe good overall computational performance for simulating non-isothermal evaporating two-fluid flows in two and three dimensions.

Finite-element method based solver for an evaporating sessile droplet on a heated surface

We report the development of a finite-element-based solver to compute transport of mass, momentum and energy during evaporation of a sessile droplet on a heated surface. The evaporation is assumed to be quasi-steady and diffusion-limited. The heat transfer between the droplet and substrate and mass transfer of liquid-vapor are solved using a two-way coupling. In particular, here, we develop and implement the formulation of fluid flow inside the droplet in the model. The continuity and Navier–Stokes equations are solved in axisymmetric, cylindrical coordinates. Jump velocity boundary condition is applied on the liquid–gas interface using the evaporation mass flux. The governing equations are discretized in the framework of the Galerkin weight residual approach. A mesh of finite triangular elements with six nodes is utilized, and quadratic shape functions are used to obtain the second-order accurate numerical solution. Two formulations, namely, penalty function and velocity pressure, are employed to obtain discretized equations. The numerical results are the same using both methods, and the latter is around 30–50% faster than the former for the cases of refined grid. Computed flow fields are in excellent agreement with published results. The solver’s capability is demonstrated by solving the internal flow field for a case of a heated substrate.

Numerical investigation of the effect of triangular cavity on the unsteady aerodynamics of NACA 0012 at a low Reynolds number

Time accurate numerical simulations were conducted to investigate the effect of triangular cavities on the unsteady aerodynamic characteristics of NACA 0012 airfoil at a Reynolds number of 50,000. Right-angled triangular cavities are placed at 10%, 25% and 50% chord location on the suction and have depths of 0.025c and 0.05c, measured normal to the surface of the airfoil. The second-order accurate solution to the RANS equations is obtained using a pressure-based finite volume solver with a four-equation transition turbulence model, γ–Re θt, to model the effect of turbulence. The two-dimensional results suggest that the cavity of depth 0.025c at 10% chord improves the aerodynamic efficiency ( l/d ratio) by 52%, at an angle of attack of α = 8°, wherein the flow is steady. The shallower triangular cavity when placed at 25%c and 50%c enhances the l/d ratio by only 10% and 17%, respectively, in the steady-state regime of angles of attack between α = 6° and 10°. The deeper cavity also enhances the l/d ratio by up to 13%, 22% and 14% at angles of attack between α = 6° and 10°. Even in the unsteady vortex shedding regime, at α =12° and higher, significant improvements in the time-averaged l/d ratios are observed for both cavity depths. The improvements in l/d ratio in the steady-state, pre-stall regime are primarily because of drag reduction while in the post-stall, unsteady regime, the improvements are because of enhancements in time-averaged C l values. The current finding can thus be used to enhance the aerodynamic performance of MAVs and UAVs that fly at low Reynolds numbers.

Equilibrium approach for modeling erosional failure of granular dams

Erosional failure of granular dams by an overtopping body of water is investigated using a depth-averaged morphodynamic model. The transport of sediment by the flow assumes the sediment flux to remain in equilibrium with the local bed shear stress. Accordingly, the shallow-water hydrodynamic equations are coupled with the Exner equation for mass conservation of the sediment. The system of equations is solved using a fully coupled well-balanced finite volume method, second-order accurate in time and space. The effect of the steep bed slope of a dam face is incorporated into both the hydrodynamics and sediment transport equations, leading to improved predictions. Comparison with results obtained from nonequilibrium sediment transport models indicates that such models perform poorly while predicting the bed evolution near the toe of an eroding dam. Observations from experimental studies demonstrate that the amount of sediment entrained by the flow is not significant, except during the initial moments of failure. This suggests that the vertical exchange of mass between the bed and the flow layer, as assumed by the nonequilibrium models, may not be completely valid during the failure. The equilibrium model results, reproducing the key flow features of the overtopping failure process, are validated by experimental measurements. The study provides fresh insights into the sediment transport processes associated with the erosion of a granular dam by overtopping, establishes the appropriateness of the equilibrium approach for its numerical modeling, and proposes a well-balanced second-order accurate solution technique for solving the resulting coupled equations of flow and sediment transport.

A Second-Order Scheme with Nonuniform Time Steps for a Linear Reaction-Subdiffusion Problem

Stability and convergence of a time-weighted discrete scheme with nonuniform time steps are established for linear reaction-subdiffusion equations. The Caupto derivative is approximated at an offset point by using linear and quadratic polynomial interpolation. Our analysis relies on two tools: a discrete fractional Gr\"{o}nwall inequality and the global consistency analysis. The new consistency analysis makes use of an interpolation error formula for quadratic polynomials, which leads to a convolution-type bound for the local truncation error. To exploit these two tools, some theoretical properties of the discrete kernels in the numerical Caputo formula are crucial and we investigate them intensively in the nonuniform setting. Taking the initial singularity of the solution into account, we obtain a sharp error estimate on nonuniform time meshes. The fully discrete scheme generates a second-order accurate solution on the graded mesh provided a proper grading parameter is employed. An example is presented to show the sharpness of our analysis.

Imposing mixed Dirichlet-Neumann-Robin boundary conditions on irregular domains in a level set/ghost fluid based finite difference framework

In this paper, an efficient, unified finite difference method for imposing mixed Dirichlet, Neumann and Robin boundary conditions on irregular domains is proposed, leveraging on our previous work [Chai et al., J. Comput. Phys. 400 (2020): 108890]. The level set method is applied to describe the arbitrarily-shaped interface, and the ghost fluid method is utilized to address the complex discontinuities on the interface. The core of this method lies in providing required ghost values under the restriction of mixed boundary conditions, which is done in a fractional-step way. Specifically, the normal derivative is calculated in the concerned subdomain by aid of a linear polynomial reconstruction, then the normal derivative and the ghost value are successively extrapolated to the other subdomain using a linear partial differential equation approach. A series of Poisson problems with mixed boundary conditions and a heat transfer test are performed to validate the method, highlighting its convergence accuracy in the L1 and L∞ norms. The method produces second-order accurate solutions with first-order accurate gradients, and is easy to implement in multi-dimensional configurations. In summary, the method represents a promising tool for imposing mixed boundary conditions, which will be applied to practical problems in future work.

Computing the finite time Lyapunov exponent for flows with uncertainties

We propose an Eulerian approach to compute the expected finite time Lyapunov exponent (FTLE) of uncertain flow fields. The definition extends the usual FTLE for deterministic dynamical systems. Instead of performing Monte Carlo simulations as in typical Lagrangian computations, our approach associates each initial flow particle with a probability density function (PDF) which satisfies an advection-diffusion equation known as the Fokker-Planck (FP) equation. Numerically, we incorporate Strang's splitting scheme so that we can obtain a second-order accurate solution to the equation. To further improve the computational efficiency, we develop an adaptive approach to concentrate the computation of the FTLE near the ridge, where the so-called Lagrangian coherent structure (LCS) might exist. We will apply our proposed algorithm to several test examples including a real-life dataset to demonstrate the performance of the method.

Integration of moment equations in a reduced-order modeling strategy for Monte Carlo simulations of groundwater flow

We illustrate and test an approach grounded on embedding moment equations (MEs) of groundwater flow within a Monte Carlo based modeling strategy to yield a Reduced-Order Model (ROM) that enables the efficient and accurate evaluation of probability distributions of hydraulic heads in randomly heterogeneous transmissivity fields. The projection space determining the accuracy of the ROM solution is typically computed through the principal component analysis of a selected number of full system model solutions (the so-called snapshots). Computationally expensive sensitivity analyses are then required to assess the independence of the ROM from the snapshots. Here, we propose to compute the projection vectors upon relying on the hydraulic head covariance evaluated from the solution of corresponding MEs of groundwater flow.Our workflow to compute hydraulic head distributions is organized according to the following steps: (i) approximation of mean hydraulic head and head covariance matrix through (second-order accurate) solutions of MEs; (ii) computation of the leading eigenvectors of the head covariance matrix to form the basis set for the ROM projection space; and (iii) construction of the ROM. Sample probability density functions of hydraulic heads are then efficiently obtained via Monte Carlo simulations relying on the developed ROM. The proposed methodology is compared against snapshot-based ROMs and the full system model in a two- and a three-dimensional steady-state groundwater flow setting where pumping from a point source is superimposed to a mean uniform flow. Our results show that the projection space computed by relying on MEs provides a more accurate ROM solution than the one resulting from reliance on snapshots.

Solving elliptic interface problems with jump conditions on Cartesian grids

We present a simple numerical algorithm for solving elliptic equations where the diffusion coefficient, the source term, the solution and its flux are discontinuous across an irregular interface. The algorithm produces second-order accurate solutions and first-order accurate gradients in the L∞-norm on Cartesian grids. The condition number is bounded, regardless of the ratio of the diffusion constant and scales like that of the standard 5-point stencil approximation on a rectangular grid with no interface. Numerical examples are given in two and three spatial dimensions.

A finite difference discretization method for heat and mass transfer with Robin boundary conditions on irregular domains

This paper proposes a finite difference discretization method for simulations of heat and mass transfer with Robin boundary conditions on irregular domains. The level set method is utilized to implicitly capture the irregular evolving interface, and the ghost fluid method to address variable discontinuities on the interface. Special care has been devoted to providing ghost values that are restricted by the Robin boundary conditions. Specifically, it is done in two steps: 1) calculate the normal derivative in cells adjacent to the interface by reconstructing a linear polynomial system; 2) successively extrapolate the normal derivative and the ghost value in the normal direction using a linear partial differential equation approach. This method produces second-order accurate solutions for both the Poisson and heat equations with Robin boundary conditions, and first-order accurate solutions for the Stefan problems. The solution gradients are of first-order accuracy, as expected. It is easy to implement in three-dimensional configurations, and can be straightforwardly generalized into higher-order variants. The method thus represents a promising tool for practical heat and mass transfer problems involving Robin boundary conditions.

Method of Order Reduction for the High-Dimensional Convection-Diffusion-Reaction Equation with Robin Boundary Conditions Based on MQ RBF-FD

A novel method of order reduction is proposed to the high-dimensional convection-diffusion-reaction equation with Robin boundary condition based on the multiquadric radial basis function-generated finite difference method (MQ RBF-FD). The main motivation is to get not only a second-order accurate solution but also a second-order accurate gradient. Key to the proposed method is introducing the intermediate variables representing the first-order derivatives to reduce the original second-order problem into an equivalent system of first-order partial differential equations. Then a discrete scheme for the latter is constructed, in which MQ RBF-FD method is applied to approximate the first-order derivatives of the original variable at the center point with decoupled method. Moreover, we can obtain an equivalent discrete scheme about the original variable and intermediate variables which can be proven all second-order convergent, that is, the convergence rate of the gradient of solution is also second-order. Finally numerical examples are presented to show the efficiency and accuracy of the proposed method.

A new high-order particle method for solving high Reynolds number incompressible flows

In this study, a new high-order particle method is proposed to solve the incompressible Navier–Stokes equations. The proposed method combines the advantages of particle and mesh methods to approximate the total and the spatial derivative terms under the Lagrangian and the Eulerian frameworks. Our aim is to avoid convective instability and increase solution accuracy at the same time. Data transfer from Lagrangian particles to Eulerian grids is realized by moving least squares interpolation. In contrast to the previously proposed method, there is no need to interpolate diffusion terms from Eulerian grids to Lagrangian particles. Therefore, the accuracy of the present solution will not be deteriorated by interpolation error. Additionally, no extra work is required to manage particles for searching procedure. Because no convection term needs to be discretized by upwinding schemes, false diffusion and dispersion errors will not be introduced, thereby increasing the solution accuracy. To verify the proposed particle method, several benchmark problems are solved to show that the present simulation is more stable, accurate, and efficient. The proposed particle method renders fourth- and second-order accurate solutions in space for velocity and pressure, respectively.

Solving Poisson-type equations with Robin boundary conditions on piecewise smooth interfaces

We present two finite volume schemes to solve a class of Poisson-type equations subject to Robin boundary conditions in irregular domains with piecewise smooth boundaries. The first scheme results in a symmetric linear system and produces second-order accurate numerical solutions with first-order accurate gradients in the L∞-norm (for solutions with two bounded derivatives). The second scheme is nonsymmetric but produces second-order accurate numerical solutions as well as second-order accurate gradients in the L∞-norm (for solutions with three bounded derivatives). Numerical examples are given in two and three spatial dimensions.

Poisson equations in irregular domains with Robin boundary conditions — Solver with second-order accurate gradients

We present a conservative method for solving the Poisson equation in irregular domains with Robin boundary conditions. Second-order accurate solutions and gradients in the L∞ norm are obtained.

A direct method for accurate solution and gradient computations for elliptic interface problems

Recently, an augmented IIM (Li et al. SIAM J. Numer. Anal. 55, 570–597 2017) has been developed and analyzed for some interface problems. The augmented IIM can provide not only a second-order accurate solution in the entire domain but also second-order accurate gradients from each side of the interface. In the augmented IIM, the PDE is reformulated and an augmented variable of co-dimension one is introduced and solved through a Schur complement system. Nevertheless, the augmented IIM is somewhat difficult to implement and understand for non-experts in the area. In this paper, a direct method (thus simpler than augmented IIM) without using the augmented variable is proposed for elliptic interface problems with a variable coefficient that can have a finite jump across an interface. The resulting finite difference scheme is the standard five-point central scheme at regular grid points, while it is a compact nine-point scheme at irregular grid points. The computed solution is second-order accurate and will be used to recover the gradient from each side of the interface to second-order accuracy. The discrete Green’s function is constructed and used to prove second-order convergence of both the solution and gradient for the one-dimensional algorithm with a piecewise constant coefficient. For the two-dimensional algorithm with a piecewise constant coefficient, we numerically prove second-order convergence of the solution by applying the discrete elliptic maximum principle using an optimization solver, while the second-order convergence of the gradient and the same conclusion for more general problems will be only demonstrated in the numerical examples.

Convergence Analysis in the Maximum Norm of the Numerical Gradient of the Shortley–Weller Method

The Shortley–Weller method is a standard central finite-difference-method for solving the Poisson equation in irregular domains with Dirichlet boundary conditions. It is well known that the Shortley–Weller method produces second-order accurate solutions and it has been numerically observed that the solution gradients are also second-order accurate; a property known as super-convergence. The super-convergence was proved in the $$L^{2}$$ norm in Yoon and Min (J Sci Comput 67(2):602–617, 2016). In this article, we present a proof for the super-convergence in the $$L^{\infty }$$ norm.

Assessment of a high-order discontinuous Galerkin method for internal flow problems. Part I: Benchmark results for quasi-1D, 2D waves propagation and axisymmetric turbulent flows

In this work we apply a high-order discontinuous Galerkin (DG) finite element method to inviscid and turbulent internal flow problems. The equations here considered are the quasi-1D, 2D Euler equations and the RANS and k−ω equations in axisymmetric coordinates. The method here proposed is designed to ensure high-order accuracy in ducts and engine-like geometries using both explicit and implicit schemes for the temporal discretization of the governing equations. Absorbing Sponge Layer (ASL) boundary conditions are implemented to minimize the reflection of out-going waves at open boundaries. A shock-capturing technique is used to control the oscillations of high-order solutions around shocks. Accurate solutions of the hyperbolic equations are performed by means of the five-stage fourth-order accurate Strong Stability Preserving Runge-Kutta scheme, while the implicit Backward-Euler scheme is adopted for efficient steady state simulations of internal turbulent flows. Two types of test-problems have been considered, one focusing on the potential of DG method to solve ideal quasi-1D and 2D waves propagation and shock phenomena that may occur in ducts, and the other on its feasibility to provide high-order accurate solutions of multi-dimensional internal turbulent flows in geometries typical of internal combustion engine (ICE) applications. To clearly illustrate the performance of the high-order DG method, the results are compared with exact solutions, experimental data and second-order accurate solutions obtained with a finite volume commercial code.

Imposing mixed Dirichlet–Neumann–Robin boundary conditions in a level-set framework

We consider the Poisson equation with mixed Dirichlet, Neumann and Robin boundary conditions on irregular domains. We describe a straightforward and efficient approach for imposing the mixed boundary conditions using a hybrid finite-volume/finite-difference approach, leveraging on the work of Gibou et al. (2002) [14], Ng et al. (2009) [30] and Papac et al. (2010) [33]. We utilize three different level set functions to represent the irregular boundary at which each of the three different boundary conditions must be imposed; as a consequence, this approach can be applied to moving boundaries. The method is straightforward to implement, produces a symmetric positive definite linear system and second-order accurate solutions in the L∞-norm in two and three spatial dimensions. Numerical examples illustrate the second-order accuracy and the robustness of the method.

Numerical simulation of continuous casting of steel alloy for different cooling ambiences and casting speeds using immersed boundary method

This article demonstrates implementation of immersed boundary method in continuous casting simulation involving boundary movement. In this methodology, the immersed boundary method is coupled with the second-order accurate finite difference solution of unsteady three-dimensional heat conduction equation. The moving molten metal front is modelled using the immersed boundary method in a Cartesian mesh framework that provides simplicity in its implementation and reduces the computational time as compared to the adaptive mesh solutions. A parallel programming paradigm using message passing interface has been implemented to obtain enhanced computational efficiency. This study has focused on capturing moving boundary during continuous casting and predicts the temperature distribution and shell thickness under different cooling ambiences and casting function. Good agreements with published data and correlations are obtained through numerical analysis. Mould-region shell thickness agrees well with Chipman–Fondersmith correlations. A new correlation has been further proposed for the delay constant at different heat extraction rates. The effects of key parameters like casting speed, convection and radiation from the continuous casting are also quantified in attempt to avail the data for optimal design of continuous caster.