Efficient solution techniques for implicit finite element schemes with flux limiters

The algebraic flux correction (AFC) paradigm is equipped with efficient solution strategies for implicit time‐stepping schemes. It is shown that Newton‐like techniques can be applied to the nonlinear systems of equations resulting from the application of high‐resolution flux limiting schemes. To this end, the Jacobian matrix is approximated by means of first‐ or second‐order finite differences. The edge‐based formulation of AFC schemes can be exploited to devise an efficient assembly procedure for the Jacobian. Each matrix entry is constructed from a differential and an average contribution edge by edge. The perturbation of solution values affects the nodal correction factors at neighbouring vertices so that the stencil for each individual node needs to be extended. Two alternative strategies for constructing the corresponding sparsity pattern of the resulting Jacobian are proposed. For nonlinear governing equations, the contribution to the Newton matrix which is associated with the discrete transport operator is approximated by means of divided differences and assembled edge by edge. Numerical examples for both linear and nonlinear benchmark problems are presented to illustrate the superiority of Newton methods as compared to the standard defect correction approach. Copyright © 2007 John Wiley & Sons, Ltd.     

Positivity‐preserving,flux‐limited finite‐difference and finite‐element methods for reactive transport

A new class of positivity‐preserving, flux‐limited finite‐difference and Petrov–Galerkin (PG) finite‐element methods are devised for reactive transport problems.The methods are similar to classical TVD flux‐limited schemes with the main difference being that the flux‐limiter constraint is designed to preserve positivity for problems involving diffusion and reaction. In the finite‐element formulation, we also consider the effect of numerical quadrature in the lumped and consistent mass matrix forms on the positivity‐preserving property. Analysis of the latter scheme shows that positivity‐preserving solutions of the resulting difference equations can only be guaranteed if the flux‐limited scheme is both implicit and satisfies an additional lower‐bound condition on time‐step size. We show that this condition also applies to standard Galerkin linear finite‐element approximations to the linear diffusion equation. Numerical experiments are provided to demonstrate the behavior of the methods and confirm the theoretical conditions on time‐step size, mesh spacing, and flux limiting for transport problems with and without nonlinear reaction. Copyright © 2003 John Wiley & Sons, Ltd.     

Gradient‐based nodal limiters for artificial diffusion operators in finite element schemes for transport equations

This paper presents new linearity‐preserving nodal limiters for enforcing discrete maximum principles in continuous (linear or bilinear) finite element approximations to transport problems with steep fronts. In the process of algebraic flux correction, the oscillatory antidiffusive part of a high‐order base discretization is decomposed into a set of internodal fluxes and constrained to be local extremum dim inishing. The proposed nodal limiter functions are designed to be continuous and satisfy the principle of linearity preservation that implies the preservation of second‐order accuracy in smooth regions. The use of limited nodal gradients makes it possible to circumvent angle conditions and guarantee that the discrete maximum principle holds on arbitrary meshes. A numerical study is performed for linear convection and anisotropic diffusion problems on uniform and distorted meshes in two space dimensions. Copyright © 2017 John Wiley & Sons, Ltd.     

A theoretical Taylor–Galerkin model for first‐order hyperbolic equation

This paper presents a class of Taylor–Galerkin (TG) finite‐element models for solving the first‐order hyperbolic equation which admits discontinuities. Five parameters are introduced for purposes of controlling stability, monotonicity and accuracy. In this paper, the total variation diminishing concept and the theory of M‐matrix are applied to construct a monotonic TG model for capturing discontinuities. To avoid making the scheme overly diffusive, we apply a flux‐corrected transport (FCT) technique of Boris and Book to overcome the difficulty with anti‐diffusive flux. In smooth flow regions, our strategyof developing the temporal and spatial high‐order TG finite‐element model is based on modified equation analysis. In regions where discontinuity is encountered, we resort to two dispersively more accurate models to make the prediction accuracy as high as that obtained in smooth cases. These models are developed using the entropy‐increasing principle and the theory of group velocity. Guided by this theory, a slower group velocity should be used ahead of the shock. To avoid a train of post‐shocks, free parameters should be chosen properly to obtain a group velocity which takes on a larger value than the exact phase velocity. In this paper, we also apply the entropy‐increasing principle to determine free parameters introduced in the finite‐element model. Under the entropy‐increasing requirement, it is mandatory that coefficients of the even and odd derivative terms shown in the modified equation should change signs alternatively in order to avoid non‐physical wiggles. Several benchmark problems have been investigated to confirm the integrity of these proposed characteristic models. Copyright © 2003 John Wiley & Sons, Ltd.     

A nonlinear ALE‐FCT scheme for non‐equilibrium reactive solute transport in moving domains

In this paper, we present a conservative, positivity‐preserving, high‐resolution nonlinear ALE‐flux‐corrected transport (FCT) scheme for reactive transport models in moving domains. The mathematical model is a convection–diffusion equation with a nonlinear flux equation on the moving channel wall. The reactive transport is assumed to have dominant Péclet and Damköhler numbers, a phenomenon that often results in non‐physical negative solutions. The scheme presented here is proven to be mass conservative in time and positive at all times for a small enough Δt. Reactive transport examples are simulated using this scheme for its validation, to show its convergence, and to compare it against the linear ALE‐FCT scheme. The nonlinear ALE‐FCT is shown to perform better than the linear ALE‐FCT schemes for large time steps. Copyright © 2014 John Wiley & Sons, Ltd.     

Multidimensional FEM‐FCT schemes for arbitrary time stepping

The flux‐corrected‐transport paradigm is generalized to finite‐element schemes based on arbitrary time stepping. A conservative flux decomposition procedure is proposed for both convective and diffusive terms. Mathematical properties of positivity‐preserving schemes are reviewed. A nonoscillatory low‐order method is constructed by elimination of negative off‐diagonal entries of the discrete transport operator. The linearization of source terms and extension to hyperbolic systems are discussed. Zalesak's multidimensional limiter is employed to switch between linear discretizations of high and low order. A rigorous proof of positivity is provided. The treatment of non‐linearities and iterative solution of linear systems are addressed. The performance of the new algorithm is illustrated by numerical examples for the shock tube problem in one dimension and scalar transport equations in two dimensions. Copyright © 2003 John Wiley & Sons, Ltd.     

Augmented upwind numerical schemes for the groundwater transport advection‐dispersion equation with local operators

When solute transport is advection‐dominated, the advection‐dispersion equation approximates to a hyperbolic‐type partial differential equation, and finite difference and finite element numerical approximation methods become prone to artificial oscillations. The upwind scheme serves to correct these responses to produce a more realistic solution. The upwind scheme is reviewed and then applied to the advection‐dispersion equation with local operators for the first‐order upwinding numerical approximation scheme. The traditional explicit and implicit schemes, as well as the Crank‐Nicolson scheme, are developed and analyzed for numerical stability to form a comparison base. Two new numerical approximation schemes are then proposed, namely, upwind–Crank‐Nicolson scheme, where only for the advection term is applied, and weighted upwind‐downwind scheme. These newly developed schemes are analyzed for numerical stability and compared to the traditional schemes. It was found that an upwind–Crank‐Nicolson scheme is appropriate if the Crank‐Nicolson scheme is only applied to the advection term of the advection‐dispersion equation. Furthermore, the proposed explicit weighted upwind‐downwind finite difference numerical scheme is an improvement on the traditional explicit first‐order upwind scheme, whereas the implicit weighted first‐order upwind‐downwind finite difference numerical scheme is stable under all assumptions when the appropriate weighting factor (θ) is assigned.     

Calculation of capillary jet

The Arbitrary Lagrangian Eulerian (ALE) framework coupled with some boundary tracking techniques is proven to be an effective method for simulation of free‐surface flows. In this paper, a special ALE framework is derived with clarification of three velocities, the notion of mesh‐frozen and field‐frozen, and the notion of tentatively inertial coordinates. A weighted integral ALE governing equations are formulated on generic coordinates and discretized with a finite element method and linear implicit time scheme. The system is solved with a discrete operator splitting technique and superposition‐based logistic parallelization. The formulation and implementation are verified through several fixed‐geometry problems and a reasonably good parallel performance is observed. Capillary jet flow is the main problem of the paper and the numerical techniques for boundary tracking are elaborated, which include an indirect boundary tracking of flux method and an iterative direct boundary tracking method. Also, a high‐order compact scheme for dynamic boundary condition and a squeeze technique for kinematic boundary condition are adopted. The axisymmetric jet breakup is studied in detail and numerical results match with the published data very well. Numerical accuracy and sensitivity are studied, including effects of element type, time scheme, compact scheme, and boundary tracking techniques. Copyright © 2010 John Wiley & Sons, Ltd.     

Toward a reduction of mesh imprinting

This paper is the initial investigation into a new Lagrangian cell‐centered hydrodynamic scheme that is motivated by the desire for an algorithm that resists mesh imprinting and has reduced complexity. Key attributes of the new approach include multidimensional construction, the use of flux‐corrected transport (FCT) to achieve second order accuracy, automatic determination of the mesh motion through vertex fluxes, and vorticity control. Toward this end, vorticity preserving Lax–Wendroff (VPLW) type schemes for the two‐dimensional acoustic equations were analyzed and then implemented in a FCT algorithm. Here, mesh imprinting takes the form of anisotropic dispersion relationships. If the stencil for the LW methods is limited to nine points, four free parameters exist. Two parameters were fixed by insisting that no spurious vorticity be created. Dispersion analysis was used to understand how the remaining two parameters could be chosen to increase isotropy. This led to new VPLW schemes that suffer less mesh imprinting than the rotated Richtmyer method. A multidimensional, vorticity preserving FCT implementation was then sought using the most promising VPLW scheme to address the problem of spurious extrema. A well‐behaved first order scheme and a new flux limiter were devised in the process. The flux limiter is unique in that it acts on temporal changes and does not place a priori bounds on the solution. Numerical results have demonstrated that the vorticity preserving FCT scheme has comparable performance to an unsplit MUSCL‐H algorithm at high Courant numbers but with reduced mesh imprinting and superior symmetry preservation. Copyright © 2014 John Wiley & Sons, Ltd.     

A finite volume scheme preserving extremum principle for convection–diffusion equations on polygonal meshes

We propose a nonlinear finite volume scheme for convection–diffusion equation on polygonal meshes and prove that the discrete solution of the scheme satisfies the discrete extremum principle. The approximation of diffusive flux is based on an adaptive approach of choosing stencil in the construction of discrete normal flux, and the approximation of convection flux is based on the second‐order upwind method with proper slope limiter. Our scheme is locally conservative and has only cell‐centered unknowns. Numerical results show that our scheme can preserve discrete extremum principle and has almost second‐order accuracy. Copyright © 2017 John Wiley & Sons, Ltd.     

A kinetic energy preserving nodal discontinuous Galerkin spectral element method

In this work, we discuss the construction of a skew‐symmetric discontinuous Galerkin (DG) collocation spectral element approximation for the compressible Euler equations. Starting from the skew‐symmetric formulation of Morinishi, we mimic the continuous derivations on a discrete level to find a formulation for the conserved variables. In contrast to finite difference methods, DG formulations naturally have inter‐domain surface flux contributions due to the discontinuous nature of the approximation space. Thus, throughout the derivations we accurately track the influence of the surface fluxes to arrive at a consistent formulation also for the surface terms. The resulting novel skew‐symmetric method differs from the standard DG scheme by additional volume terms. Those volume terms have a special structure and basically represent the discretization error of the different product rules. We use the summation‐by‐parts (SBP) property of the Gauss–Lobatto‐based DG operator and show that the novel formulation is exactly conservative for the mass, momentum, and energy. Finally, an analysis of the kinetic energy balance of the standard DG discretization shows that because of aliasing errors, a nonzero transport source term in the evolution of the discrete kinetic energy mean value may lead to an inconsistent increase or decrease in contrast to the skew‐symmetric formulation. Furthermore, we derive a suitable interface flux that guarantees kinetic energy preservation in combination with the skew‐symmetric DG formulation. As all derivations require only the SBP property of the Gauss–Lobatto‐based DG collocation spectral element method operator and that the mass matrix is diagonal, all results for the surface terms can be directly applied in the context of multi‐domain diagonal norm SBP finite difference methods. Numerical experiments are conducted to demonstrate the theoretical findings. Copyright © 2014 John Wiley & Sons, Ltd.     

A domain decomposition approach to finite volume solutions of the Euler equations on unstructured triangular meshes

We report on our recent efforts on the formulation and the evaluation of a domain decomposition algorithm for the parallel solution of two‐dimensional compressible inviscid flows. The starting point is a flow solver for the Euler equations, which is based on a mixed finite element/finite volume formulation on unstructured triangular meshes. Time integration of the resulting semi‐discrete equations is obtained using a linearized backward Euler implicit scheme. As a result, each pseudo‐time step requires the solution of a sparse linear system for the flow variables. In this study, a non‐overlapping domain decomposition algorithm is used for advancing the solution at each implicit time step. First, we formulate an additive Schwarz algorithm using appropriate matching conditions at the subdomain interfaces. In accordance with the hyperbolic nature of the Euler equations, these transmission conditions are Dirichlet conditions for the characteristic variables corresponding to incoming waves. Then, we introduce interface operators that allow us to express the domain decomposition algorithm as a Richardson‐type iteration on the interface unknowns. Algebraically speaking, the Schwarz algorithm is equivalent to a Jacobi iteration applied to a linear system whose matrix has a block structure. A substructuring technique can be applied to this matrix in order to obtain a fully implicit scheme in terms of interface unknowns. In our approach, the interface unknowns are numerical (normal) fluxes. Copyright © 2001 John Wiley & Sons, Ltd.     

Numerical modelling of gas–solid fluidized beds using the two‐fluid approach

This paper details an approach to modelling gas–solid fluidized beds using the two‐fluid granular temperature model. Details concerning the difficulties associated with the boundary conditions, particularly for curved boundaries, are described along with a novel means of obtaining the internal stress of the solid‐phase, in part, by solving an implicit equation. This results in a scheme that is stable even when the solid volume fraction is close to maximum packing. A transient, mixed finite element discretization is used to solve the multi‐phase equations with a discontinuous finite element representation of the granular temperature and continuity equations. A new solution method is proposed to solve the coupled momentum and continuity equations based on Arnoldi iteration. Two fluidized beds are modelled, one in the bubbling regime and the other in the slugging regime. These simulations are compared with experiments. Copyright © 2001 John Wiley & Sons, Ltd.     

The two‐dimensional streamline upwind scheme for the convection–reaction equation

This paper is concerned with the development of the finite element method in simulating scalar transport, governed by the convection–reaction (CR) equation. A feature of the proposed finite element model is its ability to provide nodally exact solutions in the one‐dimensional case. Details of the derivation of the upwind scheme on quadratic elements are given. Extension of the one‐dimensional nodally exact scheme to the two‐dimensional model equation involves the use of a streamline upwind operator. As the modified equations show in the four types of element, physically relevant discretization error terms are added to the flow direction and help stabilize the discrete system. The proposed method is referred to as the streamline upwind Petrov–Galerkin finite element model. This model has been validated against test problems that are amenable to analytical solutions. In addition to a fundamental study of the scheme, numerical results that demonstrate the validity of the method are presented. Copyright © 2001 John Wiley & Sons, Ltd.     

A study of time integration schemes for the numerical modelling of free surface flows

A semi‐discrete finite element methodology for the modelling of transient free surface flows in the context of Eulerian interface capturing is proposed. The focus of this study is put on the choice of an appropriate time integration strategy for the accurate modelling of the dynamics of free surfaces and of interfacial physics. It is composed of an adaptive time integration scheme for the Navier–Stokes equations, and of the implicit midpoint rule for the transport equation of the Eulerian marker variable. The adaptive scheme allows the automatic determination of a time‐step size that follows the physics of the problem under study, which facilitates the accurate modelling of stiff free surface flows. It is shown that the implicit midpoint rule reduces mass loss for each fluid. Various free surface flow problems are studied to verify and validate the proposed time integration strategy. Copyright © 2005 John Wiley & Sons, Ltd.     

An efficient edge‐based level set finite element method for free surface flow problems

We present an efficient technique for the solution of free surface flow problems using level set and a parallel edge‐based finite element method. An unstructured semi‐explicit solution scheme is proposed. A custom data structure, obtained by blending node‐based and edge‐based approaches is presented so to allow a good parallel performance. In addition to standard velocity extrapolation (for the convection of the level set function), an explicit extrapolation of the pressure field is performed in order to impose both the pressure boundary condition and the volume conservation. The latter is also improved with a modification of the divergence free constrain. The method is shown to allow an efficient solution of both simple benchmark cases and complex industrial examples. Copyright © 2012 John Wiley & Sons, Ltd.     

Annular liquid jets at high Reynolds numbers

The flow of annular liquid jets at high Reynolds numbers is analysed by means of the finite element method and the full‐Newton iteration scheme. Results have been obtained for various values of the inner to the outer diameter ratio and for non‐zero surface tension, using extremely long meshes. The annular film moves far from the symmetry axis at low values of the Reynolds number. At higher Reynolds numbers, the film moves towards the axis of symmetry and appears close to very far downstream, forming a round jet. Asymptotic results for the radius of the resulting round jet are provided. Copyright © 2003 John Wiley & Sons, Ltd.     

A projection scheme for incompressible multiphase flow using adaptive Eulerian grid: 3D validation

A three‐dimensional finite element method for incompressible multiphase flows with capillary interfaces is developed based on a (formally) second‐order projection scheme. The discretization is on a fixed (Eulerian) reference grid with an edge‐based local h‐refinement in the neighbourhood of the interfaces. The fluid phases are identified and advected using the level‐set function. The reference grid is then temporarily reconnected around the interface to maintain optimal interpolations accounting for the singularities of the primary variables. Using a time splitting procedure, the convection substep is integrated with an explicit scheme. The remaining generalized Stokes problem is solved by means of a pressure‐stabilized projection. This method is simple and efficient, as demonstrated by a wide range of difficult free‐surface validation problems, considered in the paper. Copyright © 2005 John Wiley & Sons, Ltd.     

High‐order implicit Runge–Kutta time integrators for fluid‐structure interactions

This paper presents an approach to develop high‐order, temporally accurate, finite element approximations of fluid‐structure interaction (FSI) problems. The proposed numerical method uses an implicit monolithic formulation in which the same implicit Runge–Kutta (IRK) temporal integrator is used for the incompressible flow, the structural equations undergoing large displacements, and the coupling terms at the fluid‐solid interface. In this context of stiff interaction problems, the fully implicit one‐step approach presented is an original alternative to traditional multistep or explicit one‐step finite element approaches. The numerical scheme takes advantage of an arbitrary Lagrangian–Eulerian formulation of the equations designed to satisfy the geometric conservation law and to guarantee that the high‐order temporal accuracy of the IRK time integrators observed on fixed meshes is preserved on arbitrary Lagrangian–Eulerian deforming meshes. A thorough review of the literature reveals that in most previous works, high‐order time accuracy (higher than second order) is seldom achieved for FSI problems. We present thorough time‐step refinement studies for a rigid oscillating‐airfoil on deforming meshes to confirm the time accuracy on the extracted aerodynamics reactions of IRK time integrators up to fifth order. Efficiency of the proposed approach is then tested on a stiff FSI problem of flow‐induced vibrations of a flexible strip. The time‐step refinement studies indicate the following: stability of the proposed approach is always observed even with large time step and spurious oscillations on the structure are avoided without added damping. While higher order IRK schemes require more memory than classical schemes (implicit Euler), they are faster for a given level of temporal accuracy in two dimensions. Copyright © 2015 John Wiley & Sons, Ltd.     

Some considerations for high‐order ‘incremental remap’‐based transport schemes: edges,reconstructions, and area integration

The problem of two‐dimensional tracer advection on the sphere is extremely important in modeling of geophysical fluids and has been tackled using a variety of approaches. A class of popular approaches for tracer advection include ‘incremental remap’ or cell‐integrated semi‐Lagrangian‐type schemes. These schemes achieve high‐order accuracy without the need for multistage integration in time, are capable of large time steps, and tend to be more efficient than other high‐order transport schemes when applied to a large number of tracers over a single velocity field. In this paper, the simplified flux‐form implementation of the Conservative Semi‐LAgrangian Multi‐tracer scheme (CSLAM) is reformulated using quadratic curves to approximate the upstream flux volumes and Gaussian quadrature for integrating the edge flux. The high‐order treatment of edge fluxes is motivated because of poor accuracy of the CSLAM scheme in the presence of strong nonlinear shear, such as one might observe in the midlatitudes near an atmospheric jet. Without the quadratic treatment of upstream edges, we observe at most second‐order accuracy under convergence of grid resolution, which is returned to third‐order accuracy under the improved treatment. A shallow‐water barotropic instability also reveals clear evidence of grid imprinting without the quadratic correction. Consequently, these tests reveal a problem that might arise in tracer transport near nonlinearly sheared regions of the real atmosphere, particularly near cubed‐sphere panel edges. Although CSLAM is used as the foundation for this analysis, the conclusions of this paper are applicable to the general class of incremental remap schemes. Copyright © 2012 John Wiley & Sons, Ltd.