Local Mass Conservation Research Articles

A Local Macroscopic Conservative (LoMaC) Low Rank Tensor Method for the Vlasov Dynamics

In this paper, we propose a novel Local Macroscopic Conservative (LoMaC) low rank tensor method for simulating the Vlasov-Poisson (VP) system. The LoMaC property refers to the exact local conservation of macroscopic mass, momentum and energy at the discrete level. This is a follow-up work of our previous development of a conservative low rank tensor approach for Vlasov dynamics (arXiv:2201.10397). In that work, we applied a low rank tensor method with a conservative singular value decomposition to the high dimensional VP system to mitigate the curse of dimensionality, while maintaining the local conservation of mass and momentum. However, energy conservation is not guaranteed, which is a critical property to avoid unphysical plasma self-heating or cooling. The new ingredient in the LoMaC low rank tensor algorithm is that we simultaneously evolve the macroscopic conservation laws of mass, momentum and energy using a flux-difference form with kinetic flux vector splitting; then the LoMaC property is realized by projecting the low rank kinetic solution onto a subspace that shares the same macroscopic observables by a conservative orthogonal projection. The algorithm is extended to the high dimensional problems by hierarchical Tuck decomposition of solution tensors and a corresponding conservative projection algorithm. Extensive numerical tests on the VP system are showcased for the algorithm’s efficacy.

Pressure-robust finite element scheme for the time-dependent fully coupled Stokes–Darcy-transport problem

We present a pressure-robust mixed finite element scheme for the time-dependent fully coupled Stokes–Darcy-transport problem. The P1c⊕RT0 discretization and Raviart–Thomas mixed elements are used to discretize the velocity in the free flow region and the porous medium region, respectively. We approximate the pressure with piecewise constant finite elements and the concentration with piecewise linear finite elements. We perform a divergence-free reconstruction of the velocity in the Stokes domain. Our scheme preserves local mass conservation in the sense of projection and has pressure-robust property. The stability and error analysis of the method are obtained. The errors of velocity and concentration are independent of pressure, and our method works well even with small viscosity. The theoretical results are validated with numerical examples.

Fast solution of incompressible flow problems with two-level pressure approximation

This paper develops efficient preconditioned iterative solvers for incompressible flow problems discretised by an enriched Taylor–Hood mixed approximation, in which the usual pressure space is augmented by a piecewise constant pressure to ensure local mass conservation. This enrichment process causes over-specification of the pressure when the pressure space is defined by the union of standard Taylor–Hood basis functions and piecewise constant pressure basis functions, which complicates the design and implementation of efficient solvers for the resulting linear systems. We first describe the impact of this choice of pressure space specification on the matrices involved. Next, we show how to recover effective solvers for Stokes problems, with preconditioners based on the singular pressure mass matrix, and for Oseen systems arising from linearised Navier–Stokes equations, by using a two-stage pressure convection–diffusion strategy. The codes used to generate the numerical results are available online.

Third-order relativistic fluid dynamics at finite density in a general hydrodynamic frame

The motion of water is governed by the Navier–Stokes equations, which are complemented by the continuity equation to ensure local mass conservation. In this work, we construct the relativistic generalization of these equations through a gradient expansion for a fluid with a conserved charge in a curved d-dimensional spacetime. We adopt a general hydrodynamic frame and introduce the irreducible-structure (IS) algorithm, which is based on derivatives of the expansion scalar and the shear and vorticity tensors. By this method, we systematically generate all permissible gradients up to a specified order and derive the most comprehensive constitutive relations for a charged fluid, accurate to third-order in the gradient expansion. These constitutive relations are formulated to apply to ordinary (nonconformal) and conformally invariant charged fluids. Furthermore, we examine the frame dependence of the transport coefficients for a nonconformal charged fluid up to the third order in the gradient expansion. The frame dependence of the scalar, vector, and tensor parts of the constitutive relations is obtained in terms of the (field redefinitions of the) fundamental hydrodynamic variables. Managing the frame dependencies of the constitutive relations is challenging due to their non-linear character. However, in the linear regime, the higher-order transformations become tractable, enabling the identification of a set of frame-invariant coefficients. Subsequently, the equations obtained in the linear regime are solved in momentum space, yielding dispersion relations for shear, sound, and diffusive modes for a non-conformal charged fluid, expressed in terms of a set of frame-invariant transport coefficients.

Passive-tracer modelling at super-resolution with Weather Research and Forecasting – Advanced Research WRF (WRF-ARW) to assess mass-balance schemes

Abstract. Super-resolution atmospheric modelling can be used to interpret and optimize environmental observations during top-down emission rate retrieval campaigns (e.g. aircraft-based) by providing complementary data that closely correspond to real-world atmospheric pollution transport and dispersion conditions. For this work, super-resolution model simulations with large-eddy-simulation sub-grid-scale parameterization were developed and implemented using WRF-ARW (Weather Research and Forecasting - Advanced Research WRF). We demonstrate a series of best practices for improved (realistic) modelling of atmospheric pollutant dispersion at super-resolutions. These include careful considerations for grid quality over complex terrain, sub-grid turbulence parameterization at the scale of large eddies, and ensuring local and global tracer mass conservation. The study objective was to resolve small dynamical processes inclusive of spatio-temporal scales of high-speed (e.g. 100 m s−1) airborne measurements. This was achieved by downscaling of reanalysis data from 31.25 km to 50 m through multi-domain model nesting in the horizontal and grid-refining in the vertical. Further, WRF dynamical-solver source code was modified to simulate the release of passive tracers within the finest-resolution domain. Different meteorological case studies and several tracer source emission scenarios were considered. Model-generated fields were evaluated against observational data (surface monitoring network and aircraft campaign data) and also in terms of tracer mass conservation. Results indicated agreement between modelled and observed values within 5 ∘C for temperature, 1 %–25 % for relative humidity, and 1–2 standard deviations for wind fields. Model performance in terms of (global and local) tracer mass conservation was within 2 % to 5 % of model input emissions. We found that, to ensure mass conservation within the modelling domain, tracers should be released on a regular-resolution grid (vertical and horizontal). Further, using our super-resolution modelling products, we investigated emission rate estimations based on flux calculation and mass-balancing. Our results indicate that retrievals under weak advection conditions (horizontal wind speeds < 5 m s−1) are not reliable due to weak correlation between the source emission rate and the downwind tracer mass flux. In this work we demonstrate the development of accurate super-resolution model simulations useful for planning, interpreting, and optimizing top-down retrievals, and we discuss favourable conditions (e.g. meteorological) for reliable mass-balance emission rate estimations.

Adaptive mesh refinement in locally conservative level set methods for multiphase fluid displacements in porous media

Multiphase flow in porous media often occurs with the formation and coalescence of fluid ganglia. Accurate predictions of such mechanisms in complex pore geometries require simulation models with local mass conservation and with the option to improve resolution in areas of interest. In this work, we incorporate patch-based, structured adaptive mesh refinement capabilities into a method for local volume conservation that describes the behaviour of disconnected fluid ganglia during level set simulations of capillary-controlled displacement in porous media. We validate the model against analytical solutions for three-phase fluid configurations in idealized pores containing gas, oil, and water, by modelling the intermediate-wet oil layers as separate domains with their volumes preserved. Both the pressures and volumes of disconnected ganglia converge to analytical values with increased refinement levels of the adaptive mesh. Favourable results from strong and weak scaling tests emphasize that the number of patches per processor and the total number of patches are important parameters for efficient parallel simulations with adaptive mesh refinement. Simulations of two-phase imbibition and three-phase gas invasion on segmented 3D images of water-wet sandstone show that adaptive mesh refinement has the highest impact on three-phase displacements, especially concerning the behaviour of the conserved, intermediate-wet phase.

A pressure-robust mixed finite element method for the coupled Stokes–Darcy problem

A pressure-robust numerical method for the coupled Stokes–Darcy problem is proposed. We use the continuous piecewise linear polynomial space combined with the lowest order Raviart–Thomas space to discretize the Stokes equation in the fluid region and the lowest order Raviart–Thomas elements to discretize the Darcy equation in the porous media domain. We also perform a divergence-free reconstruction of the test function for the right-hand term. Our method has good properties of convergence, pressure-robustness and local mass conservation in the sense of projection. A discrete inf-sup condition and error estimates for the numerical scheme are demonstrated, the error of velocity is independent of viscosity and pressure. Finally, the theoretical analysis is validated with several examples, demonstrating the superiority of the pressure-robust scheme.

A high‐order stabilized solver for the volume averaged Navier‐Stokes equations

AbstractThe volume‐averaged Navier‐Stokes equations are used to study fluid flow in the presence of fixed or moving solids such as packed or fluidized beds. We develop a high‐order finite element solver using both forms A and B of these equations. We introduce tailored stabilization techniques to prevent oscillations in regions of sharp gradients, to relax the Ladyzhenskaya–Babuska–Brezzi inf‐sup condition, and to enhance the local mass conservation and the robustness of the formulation. We calculate the void fraction using the particle centroid method. Using different drag models, we calculate the drag force exerted by the solids on the fluid. We implement the method of manufactured solution to verify our solver. We demonstrate that the model preserves the order of convergence of the underlying finite element discretization. Finally, we simulate gas flow through a randomly packed bed and study the pressure drop and mass conservation properties to validate our model.

A Reduced Basis Method for Darcy Flow Systems that Ensures Local Mass Conservation by Using Exact Discrete Complexes

A solution technique is proposed for flows in porous media that guarantees local conservation of mass. We first compute a flux field to balance the mass source and then exploit exact co-chain complexes to generate a solenoidal correction. A reduced basis method based on proper orthogonal decomposition is employed to construct the correction and we show that mass balance is ensured regardless of the quality of the reduced basis approximation. The method is directly applicable to mixed finite and virtual element methods, among other structure-preserving discretization techniques, and we present the extension to Darcy flow in fractured porous media.

A novel meshless method based on the virtual construction of node control domains for porous flow problems

In this paper, a novel meshless method that can handle porous flow problems with singular source terms is developed by virtually constructing node control domains. By defining a connectable node cloud, this method uses the integral of the diffusion term and generalized finite-difference operators to derive overdetermined equations of the node control volumes. An empirical method of calculating reliable node control volumes and a triangulation-based method to determine the connectable point cloud are developed. The method focuses only on the volume of the node control domain rather than the specific shape, so the construction of node control domains is virtual and does not increase the computational cost. To our knowledge, this is the first study to construct node control volumes in the meshless framework, termed the node control domain-based meshless method (NCDMM), which can also be regarded as an extended finite-volume method (EFVM). Taking two-phase porous flow problems as an example, discrete NCDMM schemes satisfying local mass conservation are derived by integrating the generalized finite-difference schemes of governing equations on each node control domain. Finally, commonly used low-order finite-volume method (FVM)-based nonlinear solvers for various porous flow models can be directly employed in the proposed NCDMM, significantly facilitating general-purpose application. Theoretically, the proposed NCDMM has the advantages of previous meshless methods for discretizing computational domains with complex geometries, as well as the advantages of traditional low-order FVMs for stably handling a variety of porous flow problems with local mass conservation. Four numerical cases are implemented to test the computational accuracy, efficiency, convergence, and good adaptability to the calculation domain with complex geometry and various boundary conditions.

An efficient discontinuous Galerkin - mixed finite element model for variable density flow in fractured porous media

Modeling variable density flow (VDF) in fractured porous media is computationally challenging. The challenges are mainly related to: (i) the high permeability contrast between the matrix and the fractures and, (ii) the nonlinearity induced by density variations. Due to their local mass conservation property, cell-centered methods, such as finite volumes (FV), mixed finite elements (MFE) or Discontinuous Galerkin (DG) methods are well suited for modeling mass transfer in highly heterogeneous domains. When applied for fractured media, these methods require small grid cells next to fractures because of the cross-flow equilibrium assumption (i.e. the pressure and concentration in the fracture and in the adjacent matrix grid-cells are assumed the same). To avoid this constraint, an efficient model is developed in this work using advanced cell-centered numerical methods for VDF in fractured porous media with cross-flow equilibrium assumed only across the fractures. The new model uses the hybrid MFE method for the flow discretization in the matrix and in the fracture continua. Mass transport in the matrix is modeled using a monotonic upwind MFE scheme. Advection-dominated transport in fractures is discretized with the DG method, which is well adapted for hyperbolic equations. The new model ensures continuity of pressure, concentration, fluid mass flux, advective and dispersive contaminant fluxes at matrix-fracture interfaces as well as at the intersection of several fractures. The time discretization is performed using high-order adaptive time integration techniques via the method of lines (MOL). The developed model is first validated by comparison against a 2D-2D model for a test problem dealing with linear flow and transport in a 2D domain involving a “+”-shaped barrier/fracture network. Then, the MFE-DG model is used for the simulation of the Henry saltwater intrusion problem and the results are validated by comparison against the semi-analytical solution in the case of unfractured and fractured aquifers. Finally, the simulation of a transport problem with variable viscosity in a fractured heterogeneous domain points out the superiority of the new model in terms of accuracy and efficiency when compared to a standard finite element model.

Global and local conservation of mass, momentum and kinetic energy in the simulation of compressible flow

The spatial discretization of convective terms in compressible flow equations is studied from an abstract viewpoint, for finite-difference methods and finite-volume type formulations with cell-centered numerical fluxes. General conditions are sought for the local and global conservation of primary (mass and momentum) and secondary (kinetic energy) invariants on Cartesian meshes. The analysis, based on a matrix approach, shows that sharp criteria for global and local conservation can be obtained and that in many cases these two concepts are equivalent. Explicit numerical fluxes are derived in all finite-difference formulations for which global conservation is guaranteed, even for non-uniform Cartesian meshes. The treatment reveals also an intimate relation between conservative finite-difference formulations and cell-centered finite-volume type approaches. This analogy suggests the design of wider classes of finite-difference discretizations locally preserving primary and secondary invariants.

An upwind-block-centered finite difference method for a semiconductor device of heat conduction and its numerical analysis

The mathematical model is formulated by a nonlinear system of initial–boundary problem including four partial differential equations: an elliptic equation for electrostatic potential, two convection–diffusion equations for electron concentration and hole concentration, a heat conduction equation for temperature. The electric field potential is solved by the conservative block-centered method, and the first order of the accuracy is improved. The concentrations and temperature are computed by the upwind-block-centered difference method. The block-centered method is used to discretize the diffusion terms. The upwind difference is applied to approximate the convection parts without numerical dispersion or nonphysical oscillation. The block-centered difference simulates the concentrations, temperature and their adjoint vector functions simultaneously. The local conservation of mass is preserved. An optimal order error estimates is obtained. Finally, numerical examples show the effectiveness and viability of this method.

An efficient bound-preserving and energy stable algorithm for compressible gas flow in porous media

Due to the great significance of the natural gas and shale gas, it is becoming increasingly important to simulate compressible gas flow in porous media. The model obeys an energy dissipation law as well as molar density must be positive and bounded in terms of the equation of state. For the purpose of eliminating nonphysical solutions as well as improving the stability in practical simulation, preservation of these properties is essential for a promising numerical method, but it is actually challenging due to the strong nonlinearity and complexity of the model. In this paper, we propose an efficient linearized numerical scheme that inherits the energy dissipation law as well as preserves the boundedness of molar density. Specifically, to treat the logarithmic type Helmholtz free energy density determined by the Peng-Robinson equation of state, we propose a novel adaptive stabilization approach involving the second derivatives of the convex energy terms. At each time step, the stabilization parameter is adaptively updated by a simple and explicit formula to ensure the energy dissipation law. The stabilized and linearized chemical potential allows to formulate the local mass conservation equation as an equivalent convection-diffusion form, and from this, an adaptive time step strategy is proposed to preserve the positivity and boundedness of molar density. The calculation of the time step size is fully explicit and easy to implement. Additionally, the fully discrete scheme is constructed using the conservative cell-centered finite difference method with the upwind strategy, and thus, it enjoys the local mass conservation. Numerical results are also presented to demonstrate the excellent performance of the proposed scheme.

An efficient and physically consistent numerical method for the Maxwell–Stefan–Darcy model of two‐phase flow in porous media

AbstractNumerical modeling of two‐phase flow in porous media has extensive applications in subsurface flow and petroleum industry. A comprehensive Maxwell–Stefan–Darcy (MSD) two‐phase flow model has been developed recently, which takes into consideration the friction between two phases by a thermodynamically consistent way. In this article, we for the first time propose an efficient energy stable numerical method for the MSD model, which can preserve multiple important physical properties of the model. First, the proposed scheme can preserve the original energy dissipation law. This is achieved through a newly‐developed energy factorization approach that leads to linear semi‐implicit discrete chemical potentials. Second, the scheme preserves the famous Onsager's reciprocal principle and the local mass conservation law for both phases by introducing different upwind strategies for two phase saturations and applying the cell‐centered finite volume method to the original formulation of the model. Third, by introducing two auxiliary phase velocities, the scheme has ability to guarantee the positivity of both saturations under proper conditions. Another distinct feature of the scheme is that the resulting discrete system is totally linear, well‐posed and unbiased for each phase. Numerical results are also provided to show the excellent performance of the proposed scheme.

Analysis of progressive fracture in fluid-saturated porous medium using splines.

Powell-Sabin B-splines are employed to model progressive fracturing in a fluid-saturated porous medium. These splines are defined on triangles and are -continuous throughout the domain, including the crack tips, so that crack initiation can be evaluated directly at the tip. On one hand, the method captures stresses and fluid fluxes more accurately than when using standard Lagrange elements, enabling a direct assessment of the fracture criterion at the crack tip and ensuring local mass conservation. On the other hand, the method avoids limitations for discrete crack analysis which adhere to isogeometric analysis. A crack is introduced directly in the physical domain. Due to the use of triangles, remeshing and crack path tracking are straightforward. During remeshing transfer of state vectors (displacement, fluid pressure) is required from the old to the new mesh. The transfer is done using a new approach which exploits a least-squares fit with the energy balance and conservation of mass as constraints. The versatility and accuracy to simulate free crack propagation are assessed for mode-I and mixed-mode fracture problems.

A fractional-step DG-FE method for the time-dependent generalized Boussinesq equations

In this work a fractional-step DG-FE method for the time-dependent generalized Boussinesq equations is proposed and analysed. The scheme is composed of two steps. In the first step the original problem is reduced into several scalar elliptic equations. An intermediate velocity and temperature are solved simultaneously. Then in the second step, the incompressibility constraint is enforced and velocity is corrected to be discretely divergence free. Moreover, the introduced elliptic term in the correction step enables the imposition of correct Dirichlet boundary conditions at each temporal step, avoiding the artificial boundary layer introduced by classical pressure correction method. DG-FE discretization strategy is utilized, in which the discontinuous Galerkin spacial discretization for flow equations is employed to obtain local mass conservation and traditional finite element spacial discretization is adopted for heat equation to reduce degrees of freedom. By choosing different symmetry and penalty parameters, SIPG-FE and NIPG-FE methods can be utilized. The consistency and stability of both methods are proved. Preliminary error estimates proving the optimal spacial order and suboptimal temporal order are carried out. Based on a different error equation and the preliminary error estimates, the optimal temporal convergence order is obtained. Numerical tests including a benchmark simulating square cavity flow are then presented, to verify the theoretical analysis and validate the method.

Hydraulic fracturing analysis in fluid-saturated porous medium.

This paper addresses fluid-driven crack propagation in a porous medium. Cohesive interface elements are employed to model the behaviour of the crack. To simulate hydraulic fracturing, a fluid pressure degree of freedom is introduced inside the crack, separate from the fluid degrees of freedom in the bulk. Powell-Sabin B-splines, which are based on triangles, are employed to describe the geometry of the domain and to interpolate the field variables: displacements and interstitial fluid pressure. Due to their -continuity, the stress and pressure gradient are smooth throughout the whole domain, enabling a direct assessment of the fracture criterion at the crack tip and ensuring local mass conservation. Due to the use of triangles, crack insertion and remeshing are straightforward and can be done directly in the physical domain. During remeshing a mapping of the state vector (displacement and interstitial fluid pressure) is required. For this, a new methodology is exploited based on a least-square fit with the energy balance and mass conservation as constraints. The accuracy to model free crack propagation is demonstrated by two numerical examples, including crack propagation in a plate with twonotches.

2-D local hp adaptive isogeometric analysis based on hierarchical Fup basis functions

In this paper, 2-D hp local adaptive procedure is developed based on Control Volume Isogeometric Analysis (CV-IGA) and Hierarchical Fup (HF) basis functions. Contrary to the most common truncated hierarchical splines, HF enables local hp adaptation because higher resolution levels do not include only basis with smaller compact support or higher frequencies, but also with higher order. Consequence of this property is spectral convergence of the proposed adaptive algorithm which is presented on classical benchmarks such as L-shape benchmark and advection dominated problems. Even in non-smooth problems, spectral convergence is achieved contrary to the application of uniform grid. CV-IGA ensures local and global mass conservation which is potentially very important for fluid mechanics problems. 2-D proposed algorithm chooses regular control volumes in parametric space at all resolution levels closely related to the Greville points (vertices) of basis functions. Therefore, methodology is very simple requiring only overlapping of control volumes in the areas where different levels are connected, while its computational cost lies between Galerkin and collocation formulations.

On Port-Hamiltonian Approximation of a Nonlinear Flow Problem on Networks

This paper deals with the systematic development of structure-preserving and robust approximations for a class of nonlinear partial differential equations on networks. The class includes, for example, gas pipe network systems described by barotropic Euler equations. Our approach is guided throughout by energy-based modeling concepts (port-Hamiltonian formalism, theory of Legendre transformation), which provide a convenient and general line of reasoning. Under mild assumptions on the approximation, local conservation of mass, an energy bound, and the inheritance of the port-Hamiltonian structure can be shown. Our approach is not limited to conventional space discretization but also covers complexity reduction of the nonlinearities by inexact integration. Thus, it can serve as a basis for structure-preserving model reduction. Combined with an energy stable time integration, we numerically demonstrate the applicability and good stability properties of the approach using the Euler equations as an example.