Flux-corrected Transport Research Articles

Development of a total variation diminishing (TVD) sea ice transport scheme and its application in an ocean (SCHISM v5.11) and sea ice (Icepack v1.3.4) coupled model on unstructured grids

Abstract. As the demand for increased resolution and complexity in unstructured sea ice models is growing, higher demands are also placed on the sea ice transport scheme. In this study, we couple the Semi-implicit Cross-scale Hydro-science Integrated System Model (SCHISM, v5.11) with Icepack (v1.3.4), the column physics package of the Los Alamos sea ice model (CICE); a key step is to implement a total variation diminishing (TVD) transport scheme for the multi-class sea ice module in the coupled model. Compared with the second-order upwind scheme and the finite-element flux-corrected transport (FEM-FCT) scheme, the TVD transport scheme is overall superior when evaluated based on conservation, accuracy, efficiency (even with very high resolution), and strict monotonicity. Although it is slightly weaker than FEM-FCT in terms of accuracy alone, the TVD scheme still outperforms the other two schemes in comprehensive performance. The new coupled model outperforms the existing single-class ice model of SCHISM in the case of Lake Superior. For the Arctic Ocean case, it successfully reproduces the long-term changes in the sea ice extent, sea ice boundary, concentration observations from satellites, and thickness from in situ measurement.

Maximum principle preserving time implicit DGSEM for linear scalar hyperbolic conservation laws

The properties of the high-order discontinuous Galerkin spectral element method (DGSEM) with implicit backward Euler time stepping are investigated for the approximation of hyperbolic linear scalar conservation equation in multiple space dimensions. We first prove that the DGSEM scheme in one space dimension preserves a maximum principle for the cell-averaged solution when the time step is large enough. This property however no longer holds in multiple space dimensions and we propose to use the flux-corrected transport (FCT) limiting [5] based on a low-order approximation using graph viscosity to impose a maximum principle on the cell-averaged solution. These results allow us to use a linear scaling limiter [58] in order to impose a maximum principle at nodal values within elements, while limiting the cell average with the FCT limiter improves the accuracy of the limited solution. Then, we investigate the inversion of the linear systems resulting from the time implicit discretization at each time step. We prove that the diagonal blocks are invertible and provide efficient algorithms for their inversion. Numerical experiments in one and two space dimensions are presented to illustrate the conclusions of the present analyses.

A cross-diffusion system modeling rivaling gangs: Global existence of bounded solutions and FCT stabilization for numerical simulation

In this paper, we study a gang territorial model consisting of two parabolic and two ordinary differential equations, where a taxis-type mechanism models that the two rivaling gangs are repelled by each other’s graffiti. Our main analytical finding shows the existence of global, bounded classical solutions. By making use of quantitative global estimates, we prove that these solutions converge to homogeneous steady states if the initial data are sufficiently small. Moreover, we perform numerical experiments which show that for different choices of parameters, the system may become diffusion- or convection-dominated, where in the former case the solutions converge toward constant steady states while in the latter case nontrivial asymptotic behavior such as segregation is observed. In order to perform these experiments, we apply a nonlinear finite element flux-corrected transport method (FEM-FCT) which is positivity-preserving. Then we treat the nonlinearities in both the system and the proposed nonlinear scheme simultaneously using fixed-point iteration.

Modeling of heat and solute transport in a fracture-matrix mine thermal energy storage system and energy storage performance evaluation

Repurposing groundwater-filled mine cavities for thermal energy storage has demonstrated promising potential to buffer the imbalance of energy supply and demand. Fractured formations are widespread in old mines due to previous excavation activities, which dominates the performance and efficiency of the mine thermal energy storage (MTES) system. Therefore, understanding the mechanisms of fluid flow, heat, and solute transport in fractured reservoirs in the MTES system is essential for its application and the assessment of environmental impact. In this study, fluid flow, heat, and solute transport in a 3D fracture-matrix MTES system are modeled simultaneously based on site-specific data in Freiberg, Germany. The simulations are conducted with a coupled Hydro-Thermo-Component process newly developed in the open-source software OpenGeoSys (OGS). To stabilize the simulation in the fracture-matrix hybrid system, a flux-corrected transport scheme is specifically implemented and applied in the model. Thus, heat and mass exchange between the embedded fractures and the matrix in the formation of MTES can be accurately simulated on the basis of the conforming mesh. In a hypothetical short-term scenario of a 15-day heating and cooling cycle of the MTES operation, the solute is transported faster than the heat in the fractured formation, indicating that disturbance of mineral compositions in the original mine water can affect a larger region than heat. The embedded fracture network in the surrounding formation is also found to increase the thermal energy storage capacity by 44%, while decreasing the recovery ratio by 14% compared to those of the unfractured formation. More energy will be transported and stored in the more fractured formation from heated mine water, but less energy can be recovered. The behavior of the MTES system strongly depends on local geology and hydraulic conditions and therefore needs to be evaluated in a site- and application-specific manner. This study provides preliminary insights into the performance and efficiency, as well as the environmental impact of the MTES operation in the short term.

Flux-corrected transport stabilization of an evolutionary cross-diffusion cancer invasion model

In the present work, we investigate a model of the invasion of healthy tissue by cancer cells which is described by a system of nonlinear PDEs consisting of a cross-diffusion-reaction equation and two additional nonlinear ordinary differential equations. We show that when the convective part of the system, the haptotaxis term, is dominant, then straightforward numerical methods for the studied system may be unstable. We present an implicit finite element method using conforming P1 or Q1 finite elements to discretize the model in space and the θ-method for discretization in time. The discrete problem is stabilized using a nonlinear flux-corrected transport approach. It is proved that both the nonlinear scheme and the linearized problems used in fixed-point iterations are solvable and positivity preserving. Several numerical experiments are presented in 2D to demonstrate the performance of the proposed method.

Numerical investigation of nonlinear effects in a standing wave thermoacoustic engine using the discontinuous Galerkin method

Understanding the nonlinear effects in thermoacoustic heat engines operating under high-pressure amplitude conditions plays a crucial role in improving their operational performance and efficiency. The discontinuous Galerkin (DG) method has several advantages, including high accuracy, adaptability, efficiency, and the ability to handle the coupling of multiple physical fields, making it suitable for accurately simulating the flow and heat transfer characteristics in thermoacoustic engines. The simulation code was developed based on MFEM, an open-source C++ library for modular finite element methods. The fully compressible Navier-Stokes equations and Fourier heat conduction equation were solved to investigate the nonlinear effects, which cannot be captured by the linear theory. To verify the accuracy of the simulation code, the numerical investigation of the thermoacoustic wave propagations induced by thermal effects in a square enclosure was first conducted. The simulation results presented the details of the generation and propagation of thermoacoustic waves, which were largely consistent with the results obtained using the flux-corrected transport (FCT) algorithm. Second, a quarter-wavelength standing wave thermoacoustic engine was simulated. The simulation results captured the process of pressure wave generation in the thermoacoustic engine, which was consistent with the pressure wave generation process in the enclosure. The results revealed the flow and heat transfer characteristics and nonlinear effects in the thermoacoustic engine under high-pressure amplitudes, which were consistent with the results obtained using finite volume method. The simulation results showed that the DG method is an effective numerical method for simulating thermoacoustic engines.

An assessment of solvers for algebraically stabilized discretizations of convection–diffusion–reaction equations

Abstract We consider flux-corrected finite element discretizations of 3D convection-dominated transport problems and assess the computational efficiency of algorithms based on such approximations. The methods under investigation include flux-corrected transport schemes and monolithic limiters. We discretize in space using a continuous Galerkin method and ℙ1 or ℚ1 finite elements. Time integration is performed using the Crank–Nicolson method or an explicit strong stability preserving Runge–Kutta method. Nonlinear systems are solved using a fixed-point iteration method, which requires solution of large linear systems at each iteration or time step. The great variety of options in the choice of discretization methods and solver components calls for a dedicated comparative study of existing approaches. To perform such a study, we define new 3D test problems for time dependent and stationary convection–diffusion–reaction equations. The results of our numerical experiments illustrate how the limiting technique, time discretization and solver impact on the overall performance.

Bound-preserving Flux Limiting for High-Order Explicit Runge–Kutta Time Discretizations of Hyperbolic Conservation Laws

We introduce a general framework for enforcing local or global maximum principles in high-order space-time discretizations of a scalar hyperbolic conservation law. We begin with sufficient conditions for a space discretization to be bound preserving (BP) and satisfy a semi-discrete maximum principle. Next, we propose a global monolithic convex (GMC) flux limiter which has the structure of a flux-corrected transport (FCT) algorithm but is applicable to spatial semi-discretizations and ensures the BP property of the fully discrete scheme for strong stability preserving (SSP) Runge–Kutta time discretizations. To circumvent the order barrier for SSP time integrators, we constrain the intermediate stages and/or the final stage of a general high-order RK method using GMC-type limiters. In this work, our theoretical and numerical studies are restricted to explicit schemes which are provably BP for sufficiently small time steps. The new GMC limiting framework offers the possibility of relaxing the bounds of inequality constraints to achieve higher accuracy at the cost of more stringent time step restrictions. The ability of the presented limiters to recognize undershoots/overshoots, as well as smooth solutions, is verified numerically for three representative RK methods combined with weighted essentially nonoscillatory (WENO) finite volume space discretizations of linear and nonlinear test problems in 1D. In this context, we enforce global bounds and prove preservation of accuracy for the linear advection equation.

A Study on Parallel 3D Elastic Wave Modeling Using Flux-corrected Transport

A Study on Parallel 3D Elastic Wave Modeling Using Flux-corrected Transport

Numerical Simulation of Streamer Discharge Modeled by Drift-Diffusion Equations Based on SETD Method

Dynamic characteristics of streamer discharge in a parallel-plane electrode system are described within the framework of a drift-diffusion model coupled with Poisson’s equation at atmospheric pressure. The spectral-element time-domain (SETD) method is applied to the solution of streamer discharge problems, and the flux-corrected-transport (FCT) algorithm is also applied to suppress possible numerical divergence or numerical oscillation. The performance of the proposed method is examined in a number of numerical tests. Then the tempo-spatial evolutions of the streamer are obtained in detail under the different situations. Numerical results show that after the initial stage of adjustment, the streamers propagate in a steady-state mode, where the characteristics of the steady state are decided by background ionization frequency which is determined by pressure and background electric field together. And space charge effect is a dominant factor for the distortion of spatial electric field, making the discharge channel expand toward both electrodes faster. These works presented here are instructive for high-voltage insulating device protection and plasma application.

Eulerian solver sensitivity for computing long-duration blast drag loading on columns

‘Long-duration’ blast waves can exert damaging drag loads on finite structural columns that are complex to characterise. In the absence of suitable drag coefficients, computational fluid dynamics (CFD) can be the only satisfactory approach for calculating blast loading. Commercially available CFD software with shockwave modelling capabilities require solving the Euler equations that assume inviscid flow. Two Eulerian solver types are generally used: Godunov and flux-corrected transport, the latter of which is commonly perceived as more accurate for shockwave problems. This study investigates the capability of both solvers for calculating drag loading from long-duration blast interaction with an I-section, verified through comparison to experimental measurements. No solver sensitivity was observed for pressures on exposed geometry surfaces, although discrepancies were observed for shielded surfaces, particularly during later stage interaction. The Godunov scheme demonstrated increased computational efficiency and reliability in comparison to the flux-corrected transport solver type, generating results in better agreement to experimental data. The results provided new understanding of solver sensitivity for modelling long-duration blast drag loading on finite targets using Eulerian CFD analyses.

A suitability analysis of transient one-dimensional two-fluid numerical models for simulating two-phase gas-liquid flows based on benchmark problems

Gas-liquid two-phase flows in pipes are present in a variety of engineering applications. These applications require transient one-dimensional mathematical models coupled with robust numerical methods to simulate the flow behavior in time and space. As there exist several possible combinations of models and numerical methods, questions arise as to the suitability of these numerical models to best accommodate the problem features. Aiming to partially answer this question, we present in this work numerical simulations of two well-known gas-liquid benchmarks, the water faucet, and the shock tube problems, carried out with two two-fluid models – 4E1P (four equations and one pressure) and 5E2P (five equations and two pressures) – in association with three numerical methods – FCT (Flux-Corrected Transport method), FORCE (First-Order Centered Scheme) and ModFORCE. The numerical methods used herein are written in a general form and do not require the solution of the associated Riemann problem. The obtained results are sought to serve as a reference to choose the best approaches, among the ones studied, to tackle some transient one-dimensional two-phase flow problems, highlighting the advantages and disadvantages of each model and method combination regarding hyperbolicity, diffusive and dispersive effects. Overall, the combination of the 5E2P model with the ModFORCE numerical method provides the most accurate results compared to the analytical solutions available.

A straightforward strategy for leak localization in two-phase gas pipelines

The leak localization (LL) functionality is an important piece of leak detection systems that continuously monitor the operation of pipelines in real-time. After a leak is detected, it must be located as soon as possible so that mitigating actions can be taken to remedy the undesirable hazardous leak consequences. The most popular methodology used to locate leaks in single-phase pipelines is based on the pressure-gradient along the line. This technique uses the SCADA data polled and collected at the ends of the line segment, along with a fluid dynamics mechanical model capable of describing the flow in steady-state, to infer the leak position. By using the flow data at the pipe ends, the LL methodology seeks the interception of the pressure gradients computed by marching forward and backward along the line to estimate the leak site. This paper presents a detailed study to extend the LL strategy based on the pressure-gradient to two-phase flows in gas pipelines. To achieve this goal, without appealing to field or experimental data, we have built a leak database for typical real-world two-phase gas pipelines under the stratified-pattern flow by means of numerical simulations. By using a transient four-equation single-pressure two-fluid model, with suitable source terms to account for the leak existence, numerical simulations emulating leaks in a pipeline were carried out using the Flux-Corrected Transport method. After being validated by the commercial simulator OLGA, the mathematical model was employed to generate a leak database contemplating leaks of three different magnitudes at three distinct locations along the line for two different scenarios considering low and high holdup contents. The performance of the LL strategy, which has used as mechanical model the steady-state parcel of the transient mathematical model described above along with the simulated data at the pipe inlet and outlet, was evaluated by computing the leak location relative errors. The obtained results have shown that the pressure-gradient methodology has presented a particularly good performance under the circumstances considered herein. It was able to locate leaks of magnitudes between 1% and 10% of the total mass flowrate at the pipe inlet with relative errors ranging from 0.26% to 10.0%.

Fully decoupled, linear and positivity-preserving scheme for the chemotaxis–Stokes equations

In this article, we develop a fully decoupled, linear and positivity-preserving finite element scheme for solving the chemotaxis–Stokes equations, which describe the biological chemotaxis phenomenon in the fluid environment. To deal with the strong coupling of the problem, we first consider a fully decoupled and linear semi-discrete scheme, in which we only need to solve several linearized sub-problems at each time step for the velocity, pressure, oxygen concentration and cell density, respectively. Then, based on the linear finite element method for the spatial discretization, the flux-corrected transport algorithm is extended to the oxygen concentration and cell density sub-problems to preserve their positivity. Finally, we show a series of numerical examples to verify its efficiency. We also find that better simulation can be obtained by regularizing the cut-off function and setting a random initial cell density in the equations.

Assessment of numerical schemes for transient, finite-element ice flow models using ISSM v4.18

Abstract. Time-dependent simulations of ice sheets require two equations to be solved: the mass transport equation, derived from the conservation of mass, and the stress balance equation, derived from the conservation of momentum. The mass transport equation controls the advection of ice from the interior of the ice sheet towards its periphery, thereby changing its geometry. Because it is based on an advection equation, a stabilization scheme needs to be employed when solved using the finite-element method. Several stabilization schemes exist in the finite-element method framework, but their respective accuracy and robustness have not yet been systematically assessed for glaciological applications. Here, we compare classical schemes used in the context of the finite-element method: (i) artificial diffusion, (ii) streamline upwinding, (iii) streamline upwind Petrov–Galerkin, (iv) discontinuous Galerkin, and (v) flux-corrected transport. We also look at the stress balance equation, which is responsible for computing the ice velocity that “advects” the ice downstream. To improve the velocity computation accuracy, the ice-sheet modeling community employs several sub-element parameterizations of physical processes at the grounding line, the point where the grounded ice starts to float onto the ocean. Here, we introduce a new sub-element parameterization for the driving stress, the force that drives the ice-sheet flow. We analyze the response of each stabilization scheme by running transient simulations forced by ice-shelf basal melt. The simulations are based on an idealized ice-sheet geometry for which there is no influence of bedrock topography. We also perform transient simulations of the Amundsen Sea Embayment, West Antarctica, where real bedrock and surface elevations are employed. In both idealized and real ice-sheet experiments, stabilization schemes based on artificial diffusion lead systematically to a bias towards more mass loss in comparison to the other schemes and therefore should be avoided or employed with a sufficiently high mesh resolution in the vicinity of the grounding line. We also run diagnostic simulations to assess the accuracy of the driving stress parameterization, which, in combination with an adequate parameterization for basal stress, provides improved numerical convergence in ice speed computations and more accurate results.

Sparse invariant domain preserving discontinuous Galerkin methods with subcell convex limiting

In this paper, we develop high-order nodal discontinuous Galerkin (DG) methods for hyperbolic conservation laws that satisfy invariant domain preserving properties using subcell flux corrections and convex limiting. These methods are based on a subcell flux corrected transport (FCT) methodology that involves blending a high-order target scheme with a robust, low-order invariant domain preserving method that is obtained using a graph viscosity technique. The new low-order discretizations are based on sparse stencils which do not increase with the polynomial degree of the high-order DG method. As a result, the accuracy of the low-order method does not degrade when used with high-order target methods. The method is applied to both scalar conservation laws, for which the discrete maximum principle is naturally enforced, and to systems of conservation laws such as the Euler equations, for which positivity of density and a minimum principle for specific entropy are enforced. Numerical results are presented on a number of benchmark test cases.

FCT finite difference forward simulation of the submarine viscoelastic seismic wavefield

The actual seafloor is the typical viscoelastic medium, covered with unconsolidated geological bodies like sand, gravel, sludge, and detritus on the surface. When the seismic wave is propagating therein, viscosity of the seafloor strata will cause seismic wave energy loss, amplitude attenuation and gradual frequency reduction at the wavelet center, and then the accurate subsurface information and high-resolution images cannot be obtained directly from the seismic data due to the reduced signal resolution and SNR. Considering seabed viscoelastic media properties, a model including quality factor, density, and other physical parameters, the flux corrected transport (FCT) algorithm is established to describe the target seafloor and solve the wave equation, in which, the finite difference numerical simulation method of FCT is used to suppress the numerical dispersion generated by differential calculation in coarse grid, improve the seismic wavefield simulation accuracy and computational efficiency. Based on the submarine ancient channel types and the basic characteristics of the real and typical marine geological phenomena, the seismic responses to the three types of submarine ancient channel models, which are straight, curved and braided, are calculated, with good results acquired.

Flux-corrected transport for scalar hyperbolic conservation laws and convection-diffusion equations by using linear programming

Flux-corrected transport (FCT) is one of the flux limiter methods. Unlike the total variation diminishing methods, obtaining the known FCT formulas for computing flux limiters is not quite transparent, and their transformation is not obvious when the original differential operator changes. We propose a novel formal mathematical approach to design flux correction for weighted hybrid difference schemes by using linear programming. The hybrid scheme is a linear combination of a monotone scheme and a high order scheme. The determination of maximal antidiffusive fluxes is treated as an optimization problem with a linear objective function. To obtain constraints for the optimization problem, inequalities that are valid for the monotone difference scheme are applied to the hybrid difference scheme. The numerical solution of the nonlinear optimization problem is reduced to the iterative solution of linear programming problems. A nontrivial approximate solution of the corresponding linear programming problem can be treated as the required flux limiters. We present flux correction formulas for scalar hyperbolic conservation laws and convection-diffusion equations. The designed flux-corrected transport for scalar hyperbolic conservation laws yields entropy solutions. Numerical results are presented.

Locally bound-preserving enriched Galerkin methods for the linear advection equation

In this work, we introduce algebraic flux correction schemes for enriched (P1⊕P0 and Q1⊕P0) Galerkin discretizations of the linear advection equation. The piecewise-constant component stabilizes the continuous Galerkin approximation without introducing free parameters. However, violations of discrete maximum principles are possible in the neighborhood of discontinuities and steep fronts. To keep the cell averages and the degrees of freedom of the continuous P1/Q1 component in the admissible range, we limit the fluxes and element contributions, the complete removal of which would correspond to first-order upwinding. The first limiting procedure that we consider in this paper is based on the flux-corrected transport (FCT) paradigm. It belongs to the family of predictor-corrector algorithms and requires the use of small time steps. The second limiting strategy is monolithic and produces nonlinear problems with well-defined residuals. This kind of limiting is well suited for stationary and time-dependent problems alike. The need for inverting consistent mass matrices in explicit strong stability preserving Runge–Kutta time integrators is avoided by reconstructing nodal time derivatives from cell averages. Numerical studies are performed for standard 2D test problems.

Subcell flux limiting for high-order Bernstein finite element discretizations of scalar hyperbolic conservation laws

This work extends the concepts of algebraic flux correction and convex limiting to continuous high-order Bernstein finite element discretizations of scalar hyperbolic problems. Using an array of adjustable diffusive fluxes, the standard Galerkin approximation is transformed into a nonlinear high-resolution scheme which has the compact sparsity pattern of the piecewise-linear or multilinear subcell discretization. The representation of this scheme in terms of invariant domain preserving states makes it possible to prove the validity of local discrete maximum principles under CFL-like conditions. In contrast to predictor-corrector approaches based on the flux-corrected transport methodology, the proposed flux limiting strategy is monolithic, i.e., limited antidiffusive terms are incorporated into the well-defined residual of a nonlinear (semi-)discrete problem. A stabilized high-order Galerkin discretization is recovered if no limiting is performed. In the limited version, the compact stencil property prevents direct mass exchange between nodes that are not nearest neighbors. A formal proof of sparsity is provided for simplicial and box elements. The involved element contributions can be calculated efficiently making use of matrix-free algorithms and precomputed element matrices of the reference element. Numerical studies for Q2 discretizations of linear and nonlinear two-dimensional test problems illustrate the virtues of monolithic convex limiting based on subcell flux decompositions.