Discontinuous Pressure Research Articles

A divergence‐conforming method for flow and double‐diffusive transport

AbstractThe results of a recent extension of the analysis of an ‐conforming method for a model of double‐diffusive flow in porous media introduced in [Bürger, Méndez, Ruiz‐Baier, SINUM (2019), 57:1318–1343] to the time‐dependent case are summarized. These include the efficiency and reliability of residual‐based a posteriori error estimators for the steady, semi‐discrete, and fully discrete problems. The method consists of Brezzi–Douglas–Marini approximations for velocity and compatible piecewise discontinuous pressures, whereas Lagrangian elements are used for concentration and salinity. Novel numerical tests confirm the accuracy of the method and illustrate its application to a salinity‐driven problem of sedimentation.

Mixed displacement–pressure formulations and suitable finite elements for multimaterial problems with compressible and incompressible models

Multimaterial problems in linear and nonlinear elasticity are some of the least explored using mixed finite element formulations with higher-order elements. The fundamental issue in adapting the mixed displacement–pressure formulations with linear and higher-order continuous elements for the pressure field is their inability to capture pressure and stress jumps across material interfaces. In this paper, for the first time in literature, we perform comprehensive studies of multimaterial problems in elasticity consisting of compressible and incompressible material models using the mixed displacement–pressure formulation to assess the performance of different element types in accurately resolving pressure fields within the domains and pressure jumps across material interfaces. In particular, inf–sup stable displacement–pressure combinations with element-wise discontinuous pressure for triangular and tetrahedral elements are considered and their performance is assessed along with the Q1/P0 element and Taylor–Hood elements using several numerical examples. The results show that Taylor–Hood elements fail to capture the stress jumps due to the continuity of DOFs across elements, the Crouzeix–Raviart (P2b/P1dc) element yields substantially poor pressure fields despite a significant increase in pressure degrees of freedom and that the P3/P1dc element produces superior quality results fields when compared with the P2b/P1dc element.

The pressure-wired Stokes element: a mesh-robust version of the Scott–Vogelius element

The Scott–Vogelius finite element pair for the numerical discretization of the stationary Stokes equation in 2D is a popular element which is based on a continuous velocity approximation of polynomial order k and a discontinuous pressure approximation of order k-1\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$k-1$$\\end{document}. It employs a “singular distance” (measured by some geometric mesh quantity Θz≥0\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$ \\Theta \\left( \ extbf{z}\\right) \\ge 0$$\\end{document} for triangle vertices z\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$\ extbf{z}$$\\end{document}) and imposes a local side condition on the pressure space associated to vertices z\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$\ extbf{z}$$\\end{document} with Θz=0\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$\\Theta \\left( \ extbf{z}\\right) =0$$\\end{document}. The method is inf-sup stable for any fixed regular triangulation and k≥4\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$k\\ge 4$$\\end{document}. However, the inf-sup constant deteriorates if the triangulation contains nearly singular vertices 0<Θz≪1\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$0<\\Theta \\left( \ extbf{z}\\right) \\ll 1$$\\end{document}. In this paper, we introduce a very simple parameter-dependent modification of the Scott–Vogelius element with a mesh-robust inf-sup constant. To this end, we provide sharp two-sided bounds for the inf-sup constant with an optimal dependence on the “singular distance”. We characterise the critical pressures to guarantee that the effect on the divergence-free condition for the discrete velocity is negligibly small, for which we provide numerical evidence.

A 2D continuous-discontinuous unified pore-fracture seepage model for hydro-mechanical coupling in the combined finite-discrete element method (FDEM) using node binding scheme

The combined finite-discrete element method (FDEM) has been widely used to simulate the dynamic propagation of hydraulic fractures. In the conventional FDEM, a fictitious fluid exchange coefficient between adjacent finite element pairs needs to be specified to alleviate the discontinuous pore seepage in the continuous rock matrix domain. This commonly employed approach inevitably increases computational cost and introduces additional numerical parameters (e.g., fluid exchange coefficient). To circumvent these limitations, we propose a node binding scheme that not only ensures continuous pore seepage simulations in the intact/continuous rock matrix domain, but also allows discontinuous pore pressures on the two sides of a fracture. The computational efficiency of pore seepage simulation using the proposed approach is increased by ∼120 times for a model with ∼6000 triangle elements, and fictitious numerical parameters are not required. A suite of numerical benchmarks is conducted to verify the effectiveness and robustness of the proposed pore-fracture seepage model, and its capability of capturing the dynamic propagation process of fluid-driven fractures is also demonstrated. As an application, a typical example is employed to explore the mechanism of multiple hydraulic fracturing. The proposed hydro-mechanical coupling model may help enhance the applicability and accuracy of FDEM for unconventional energy exploitation.

A single mesh approximation for modeling multiphase flow in heterogeneous porous media

The control volume finite element (CVFE) approach is based on the discretization of primary flow and transport unknowns on two different meshes. The element mesh, used to represent the physical properties of the medium, differs from the control volume mesh necessary to ensure a mass conservative solution. The inherent two mesh feature of the CVFE approximation introduces inconsistency in the transport solution due to the averaging of physical and computed quantities between neighboring elements in the mesh. In this work, we present a consistent approach for modeling multiphase flow and transport in heterogeneous porous media. The approach is applied in the CVFE framework by enabling the discretization of primary unknowns on a single mesh. We combine the discretization of an element-wise, discontinuous pressure approximation with a first-order, discontinuous velocity approximation to resolve the elliptic or parabolic flow problem. The effectiveness of the formulation is achieved by exploiting the same finite element mesh as the flow problem, when updating the saturation solution. This direct mapping between the flow and transport mesh is simple yet effective for establishing a consistent solution in addition to circumventing the non-physical mass leakage exhibited in the classical CVFE method. We describe the interface approximations and the discontinuous terms needed for consistent solutions. The method is well suited to model flow in complex geometrical subsurface domains and is shown to be numerically stable while providing locally and globally mass conservative solutions. We apply the approach to several domains with complex geometrical features to emphasize the superiority of the approximation over conventional methods. The analysis shows over two orders of magnitude reduction in solution error compared to the classical CVFE. The new formulation effectively captures accurate transport solutions, even in the presence of varying material properties, without the need for mesh refinement.

Discontinuous Atmospheric Pressure Interface in-Trap Single Photon Ionization Ion Trap Mass Spectrometry (DAPI in-Trap SPI ITMS)

A compact in-trap single photon ionization (SPI) system, termed discontinuous atmospheric pressure interface (DAPI) ion trap mass spectrometry (ITMS), for the rapid determination of trace volatile organic compounds (VOCs). DAPI enables the maintenance of high analyte pressure in the SPI source while preserving the vacuum. Utilizing vacuum ultraviolet (VUV) lamp as the ion source, we ionize the analyte within an enclosed miniature 3D ion trap for efficient analysis. Experimental findings demonstrate that the system achieves a mass resolution of less than 0.5 atomic mass units for full width at half maxima (FWHM), providing unit resolution. For the response for benzene, the system achieves almost instantaneous response, with rise and fall times of only 2 s, which is a significant advantage compared to the membrane inlet (MI) approach necessary for on-site detection. Additionally, tandem mass spectrometry (MS/MS) with unit m/z isolation window was achieved and resolved dimethyl disulfide and phenol with similar molecular weights. Our results demonstrate that the system provides excellent performance and we anticipate its application in rapidly determining VOCs.

Analysis and computation of a weak Galerkin scheme for solving the 2D/3D stationary Stokes interface problems with high-order elements

Abstract In this paper, we present a high-order weak Galerkin finite element method (WG-FEM) for solving the stationary Stokes interface problems with discontinuous velocity and pressure in ℝ d , d = 2, 3. This WG method is equipped with stable finite elements consisting of usual polynomials of degree k ⩾ 1 for the velocity and polynomials of degree k − 1 for the pressure, both are discontinuous. Optimal convergence rates of order k + 1 for the velocity and order k for the pressure are established in L 2-norm on hybrid meshes. Numerical experiments verify the expected order of accuracy for both two-dimensional and three-dimensional examples. Moreover, numerically it is shown that the proposed WG algorithm is able to accommodate geometrically complicated and very irregular interfaces having sharp edges, cusps, and tips.

Development and characterization of discontinuous atmospheric pressure interface - Dual pressure chamber miniature mass spectrometer

Miniaturized mass spectrometers have become increasingly prevalent for real-time detection and analysis, owing to their compact size and portability. The pursuit of performance enhancement in these instruments is a pivotal objective within the domain of mass spectrometry miniaturization. This study introduces a novel miniature mass spectrometer featuring a discontinuous atmospheric pressure interface and a dual pressure chamber. Compared to conventional single-chamber, discontinuous sampling interface mass spectrometers, the newly developed instrument demonstrates a more than tenfold improvement in detection efficiency. This significant enhancement is achieved without the need for complex control of switch coupling time series, thereby streamlining the circuit design and improving the instrument's fault tolerance. Furthermore, by capitalizing on the benefits of discontinuous sampling, the instrument reduces the operational pressure relative to traditional continuous sampling in differential pressure vacuum chambers. It accommodates larger inlet capillary (0.38 mm) and skimmer (0.5 mm) diameters, leading to a ninefold increase in response strength for risperidone and lowering the detection limit to 0.5 ppb. The instrument's capacity for rapid drug detection, along with enhanced resolution and detection limits, underscores its potential utility. Additionally, it facilitates the use of smaller mechanical pumps, significantly diminishing both the instrument's volume and power consumption. This presents a promising avenue for further miniaturization of mass spectrometers.

Numerical analysis of a mixed-dimensional poromechanical model with frictionless contact at matrix–fracture interfaces

We present a complete numerical analysis for a general discretization of a coupled flow–mechanics model in fractured porous media, considering single-phase flows and including frictionless contact at matrix–fracture interfaces, as well as nonlinear poromechanical coupling. Fractures are described as planar surfaces, yielding the so-called mixed- or hybrid-dimensional models. Small displacements and a linear elastic behavior are considered for the matrix. The model accounts for discontinuous fluid pressures at matrix–fracture interfaces in order to cover a wide range of normal fracture conductivities. The numerical analysis is carried out in the Gradient Discretization framework (see J. Droniou, R. Eymard, T. Gallouët, C. Guichard, and R. Herbin [The gradient discretisation method, Springer, Cham, 2018]), encompassing a large family of conforming and nonconforming discretizations. The convergence result also yields, as a by-product, the existence of a weak solution to the continuous model. A numerical experiment in 2D is presented to support the obtained result, employing a Hybrid Finite Volume scheme for the flow and second-order finite elements ( P 2 \mathbb {P}_2 ) for the mechanical displacement coupled with face-wise constant ( P 0 \mathbb P_0 ) Lagrange multipliers on fractures, representing normal stresses, to discretize the contact conditions.

Multiple wave scattering by submerged obstacles in an infinite channel of finite depth with surface pressure excess

The objective is to study the combined effect of an incident wave, a surface pressure excess and a finite number of submerged obstacles, in the phenomenon of power transfer to an infinite fluid layer of finite depth. The incident wave and the surface pressure excess have the same harmonic time dependence, a fact that allows to eliminate time altogether and consider only steady-state solutions. The surface pressure excess simulates the effect of winds blowing above the water surface in oceans. The technique used in a first part of the paper relying upon the use of finite Fourier transform and separation of variables is extended here to this end. The method allows to separate local perturbations from progressive or standing wave. Our formulae yield the exact solution in closed form in the absence of obstacles, and provide a clearer insight into the flow properties, as compared to previous investigations. Applications are given for discontinuous surface pressure functions. We put in evidence solutions with no outgoing waves, for which the energy transmitted by the surface pressure is exhausted in generating a standing wave, together with local perturbations. Two numerical applications without/with obstacles, for a parabolic surface pressure profile, allow to assess the energy transfer from the pressure-obstacles system to the fluid. The results may be of interest in the field of oscillating water columns and, generally, water power converting technology.

Divergence-conforming methods for transient double-diffusive flows: a priori and a posteriori error analysis

Abstract The analysis of an $\textbf {H}(\textrm {div})$-conforming method for a model of double-diffusive flow in porous media introduced in Bürger, Méndez &amp; Ruiz-Baier (2019, On H(div)-conforming methods for double-diffusion equations in porous media. SIAM J. Numer. Anal., 57,1318–1343) is extended to the time-dependent case. In addition, the efficiency and reliability of residual-based a posteriori error estimators for the steady, semidiscrete and fully discrete problems are established. The resulting methods are applied to simulate the sedimentation of small particles in salinity-driven flows. The method consists of Brezzi–Douglas–Marini approximations for velocity and compatible piecewise discontinuous pressures, whereas Lagrangian elements are used for concentration and salinity distribution. Numerical tests confirm the properties of the proposed family of schemes and of the adaptive strategy guided by the a posteriori error indicators.

A WENO SPH scheme with improved transport velocity and consistent divergence operator

The Arbitrary Lagrangian–Eulerian Smoothed Particle Hydrodynamics (ALE-SPH) formulation can guarantee stable solutions preventing the adoption of empirical parameters such as artificial viscosity. However, the convergence rate of the ALE-SPH formulation is still limited by the inaccuracy of the SPH spatial operators. In this work, a Weighted Essentially Non-Oscillatory (WENO) spatial reconstruction is then adopted to minimise the numerical diffusion introduced by the approximate Riemann solver (which ensures stability), in combination with two alternative approaches to restore the consistency of the scheme: corrected divergence SPH operators and the particle regularisation guaranteed by the correction of the transport velocity. The present work has been developed in the framework of the DualSPHysics open-source code. The beneficial effect of the WENO reconstruction to reduce numerical diffusion in ALE-SPH schemes is first confirmed by analysing the propagation of a small pressure perturbation in a fluid initially at rest. With the aid of a 2-D vortex test case, it is then demonstrated that the two aforementioned techniques to restore consistency effectively reduce saturation in the convergence to the analytical solution. Moreover, high-order (above second) convergence is achieved. Yet, the presented scheme is tested by means of a circular blast wave problem to demonstrate that the restoration of consistency is a key feature to guarantee accuracy even in the presence of a discontinuous pressure field. Finally, a standing wave has been reproduced with the aim of assessing the capability of the proposed approach to simulate free-surface flows.

Non-overlapping domain decomposition algorithms with only primal velocity unknowns for the discontinuous viscosity Stokes problem

In this paper, non-overlapping domain decomposition algorithms for the Stokes problem with a discontinuous viscosity coefficient are proposed and analyzed. A pair of inf–sup stable finite element spaces with discontinuous pressure finite element functions is employed to obtain a discrete Stokes system. By introducing a non-overlapping subdomain partition for the given finite element spaces, the discrete Stokes system is reformulated into an algebraic system on dual unknowns, where the continuity on the velocity unknowns at the subdomain vertices is enforced strongly and the continuity on the remaining interface velocity unknowns is enforced weakly using the dual unknowns. The resulting algebraic system on the dual unknowns is then solved by an iterative method. The coarse problem in the resulting algebraic system is only related to the strongly coupled primal velocity unknowns and it is thus obtained as a symmetric positive definite matrix equation. To accelerate the iteration convergence, practical versions of lumped and Dirichlet preconditioners are proposed and their condition numbers are also analyzed for the resulting dual algebraic system. Numerical results are also included to confirm our condition number estimates.

A new diatom-based multimetric index to assess lake ecological status

Eutrophication impairs lake ecosystems at a global scale. In this context, as benthic microalgae are well-established warnings for a large range of stressors, particularly nutrient enrichment, the Water Framework Directive required the development of diatom-based methods to monitor lake eutrophication. Here, we present the diatom-based index we developed for French lakes, named IBDL (Indice Biologique Diatomées en Lacs). Data were collected in 93 lakes from 2015 to 2020. A challenge arose from the discontinuous pressure gradient of our dataset, especially the low number of nutrient-impacted lakes. To analyze the data we opted for the so-called “Threshold Indicator Taxa ANalysis” method, which makes it possible to determine a list of “alert taxa.” We obtained a multimetric index based on specific pressure gradients (Kjeldahl nitrogen, suspended matter, biological oxygen demand, and total phosphorous). Considering the European intercalibration process, the very good correlation between IBDL and the common metric (R2 from 0.52 to 0.87 according to the lake alkalinity type) makes us very confident in our ability to match future IBDL quality thresholds with European standards. The IBDL proved at last to be particularly relevant as it has a twofold interest: an excellent relationship with total phosphorus (R2 from 0.63 to 0.83 according to the lake alkalinity type) and a possible application to any lake metatype. Its complementarity with macrophyte-based indices moreover justifies the use of at least two primary producer components for lake ecological status classification.

Anisotropic Evolution of Effective Stress and Pore Pressure during Coalbed Methane Drainage

Summary The anisotropy and dynamic variation in permeability of gas-adsorbing coals have a significant influence on fluid flow behavior in the cleat system. The assumption of a constant anisotropy coefficient (the ratio between permeability components in orthogonal directions) has been traditionally made to simplify the seepage-stress coupling analytical model. In this approach, the pressure drop of the coalbed is separated into desorption and nondesorption areas. To evaluate the effective stress, pore pressure, permeability distribution, and variable anisotropy coefficient more accurately, analytical formulas were developed that consider elastic mechanics and methane sorption. The results show that the anisotropy coefficient can be dynamic when cleat compressibility anisotropy exists. Pressure contours are a set of ellipses that increase in eccentricity from the near-wellbore area to the pressure drop boundary, leading to corresponding anisotropy changes in effective stress and permeability. The gas desorption-related matrix shrinkage effect causes a discontinuous pressure drop gradient at the boundary between desorption and nondesorption areas, resulting in nonsmooth pressure drop curves. The pressure gradient difference changes with the radius of the desorption area and is nonisotropic, with the high-permeability direction showing a greater difference than the low-permeability direction. These results indicate that the dynamic anisotropy coefficient has a significant impact on coalbed drainage and extraction. Compared to previous mathematical models, which assumed permeability isotropy or constant anisotropy coefficient in cleat systems, the proposed model provides a more accurate method to evaluate pressure and permeability distribution.

A New Immersed Finite Element Method for Two-Phase Stokes Problems Having Discontinuous Pressure

Abstract In this paper, we develop a new immersed finite element method (IFEM) for two-phase incompressible Stokes flows. We allow the interface to cut the finite elements. On the noninterface element, the standard Crouzeix–Raviart element and the P 0 {P_{0}} element pair is used. On the interface element, the basis functions developed for scalar interface problems (Kwak et al., An analysis of a broken P 1 {P_{1}} -nonconforming finite element method for interface problems, SIAM J. Numer. Anal. (2010)) are modified in such a way that the coupling between the velocity and pressure variable is different. There are two kinds of basis functions. The first kind of basis satisfies the Laplace–Young condition under the assumption of the continuity of the pressure variable. In the second kind, the velocity is of bubble type and is coupled with the discontinuous pressure, still satisfying the Laplace–Young condition. We remark that in the second kind the pressure variable has two degrees of freedom on each interface element. Therefore, our methods can handle the discontinuous pressure case. Numerical results including the case of the discontinuous pressure variable are provided. We see optimal convergence orders for all examples.

Unfitted finite element methods based on correction functions for Stokes flows with singular forces of low regularity

This paper presents an unfitted finite element method based on correction functions for solving stationary Stokes flows with singular forces acting on an immersed interface. It has been shown that the singular force is equivalent to a nonhomogeneous jump condition on the interface. In this paper, we consider the case that the jump has a low regularity so that it is impossible to use pointwise values on the interface to construct correction functions, as done in Guzmán et al. (2016). The natural way to deal with the problem is to use mean values of the jump on the parts of the interface cut by elements, instead of using pointwise values. However, we show that it may cause instability and the constant in the error estimate may depend on the interface location relative to the mesh. Inspired by Guo et al. (2019), we use a larger fictitious circle to overcome these issues. Associated with the correction functions, we consider two stable finite element pairs: the Mini element and the P2−P0 element, including the cases of continuous and discontinuous pressures. The optimal approximation capabilities of the correction functions and optimal error estimates of the finite element methods are both derived with a hidden constant independent of the interface location relative to the mesh. Numerical examples are provided to validate the theoretical results.

Improved a priori error estimates for a discontinuous Galerkin pressure correction scheme for the Navier–Stokes equations

AbstractThe pressure correction scheme is combined with interior penalty discontinuous Galerkin method to solve the time‐dependent Navier–Stokes equations. Optimal error estimates are derived for the velocity in the L2 norm in time and in space. Error bounds for the discrete time derivative of the velocity and for the pressure are also established. The analysis is challenging and technical and is based on appropriate use of lift operators and duality arguments.

A hybrid pressure approximation in the control volume finite element method for multiphase flow and transport in heterogeneous porous media

We present a new hybrid pressure formulation applied to the control volume finite element (CVFE) method to model multiphase flow and transport in highly heterogeneous porous media. The formulation effectively captures sharp saturation changes in the presence of discontinuous material properties by employing a discontinuous pressure approximation at material interfaces. The heterogeneous porous medium is divided into sub-domains within which material properties are uniform or smoothly varying. By construction, the resultant control volume dual mesh is restricted within a sub-domain. The artificial mass leakage across material property boundaries observed in classical CVFE methods is therefore circumvented. The approach applies the robust continuous pressure approximation in the rest of the computational domain; the discontinuous approximation is applied only at the sub-domain boundaries. The discontinuous parameters necessary to achieve mass conservative solutions, locally and globally, are described. We demonstrate the accuracy and efficiency of the new approach by comparison with the classical continuous CVFE method on various examples of heterogeneous domains as well as establishing the convergence of the numerical method. The proposed hybrid formulation significantly outperforms the accuracy and efficiency of classical CVFE methods that use the same order of approximation for modeling multiphase flow in heterogeneous porous media.

On pressure robustness and independent determination of displacement and pressure in incompressible linear elasticity

We investigate the possibility to determine the divergence-free displacement uindependently from the pressure reaction p for a class of boundary-value problems in incompressible linear elasticity. If not possible, we investigate if it is possible to determine it pressure robustly, i.e. pollution free from the pressure reaction.For convex domains there is but one variational boundary value problem among the investigated that allows the independent determination. It is the one with essential no-penetration conditions combined with homogeneous tangential traction conditions.Further, for the investigated cases, the weakly divergence-free displacement can be computed pressure robustly provided the total body force is decomposed into its direct sum of divergence- and rotation-free components using the Helmholtz decomposition. The elasticity problem is solved using these components as separate right-hand sides. The total solution is obtained using the superposition principle.We employ a (u,p) higher-order finite element formulation with discontinuous pressure elements. It is inf–sup stable for polynomial degree p≥2 but not pressure robust by itself. We propose a three step procedure to solve the elasticity problem which includes by the Helmholtz decomposition of the total body force. The extra cost for the three-step procedure is essentially the cost for the Helmholtz decomposition of the assembled total body force, and the small cost of solving the elasticity problem with one extra right-hand side. The results are corroborated by theoretical derivations as well as numerical results.