A continuous finite element framework for the pressure Poisson equation allowing non‐Newtonian and compressible flow behavior

Computing pressure fields from given flow velocities is a task frequently arising in engineering, biomedical, and scientific computing applications. The so‐called pressure Poisson equation (PPE) derived from the balance of linear momentum provides an attractive framework for such a task. However, the PPE increases the regularity requirements on the pressure and velocity spaces, thereby imposing theoretical and practical challenges for its application. In order to stay within a Lagrangian finite element framework, it is common practice to completely neglect the influence of viscosity and compressibility when computing the pressure, which limits the practical applicability of the pressure Poisson method. In this context, we present a mixed finite element framework which enables the use of this popular technique with generalized Newtonian fluids and compressible flows, while allowing standard finite element spaces to be employed for the unknowns and the given data. This is attained through the use of appropriate vector calculus identities and simple projections of certain flow quantities. In the compressible case, the mixed formulation also includes an additional equation for retrieving the density field from the given velocities so that the pressure can be accurately determined. The potential of this new approach is showcased through numerical examples.

provides a direct relation between flow velocity and pressure gradient, its vector-valued nature makes it a not so convenient equation when the goal is to solve for pressure, a scalar quantity. The most intuitive approach is to integrate the pressure gradient tangentially along some path that connects a point with known pressure to the point of interest. Such methods are known to perform poorly in the presence of noisy velocity measurements, due to excessive error accumulation along the integration path. 2 A popular alternative is to work with the pressure Poisson equation. By taking the divergence of both sides of the momentum equation, a scalar Poisson equation is generated for the pressure. In different formulations and discretizations, the PPE has become a standard in pressure retrieval for engineering and biomedical applications. [3][4][5][6][7][8][9][10] Other numerical applications of the PPE include stabilization of equal-order finite element methods 11 and the decoupling of pressure and velocity for efficient time-stepping. 12,13 Despite its apparent simplicity, the PPE brings an important numerical challenge: the additional regularity requirements induced by the divergence operator. Even in weak form, second-order spatial differentiation of the velocity field is required in order to account for viscous effects. In that case, C 1 interpolation of the velocity would be needed, 8 a not at all trivial task in complex geometries. Thus, it is common practice to neglect the viscous term completely, 14 which is generally reasonable in large arteries 15 but might lead to inaccurate estimates in slow-flow regions such as aneurysms. 16 Nonetheless, alternatives to the PPE for pressure estimation often introduce additional variables, increasing the computational cost considerably. For instance, the method proposed by Švihlová et al. 1 introduces an artificial vector-valued variable to create a Stokes-like system, thereby increasing the total number of degrees of freedom by a factor of at least four (comparison between different techniques for pressure retrieval are available in the literature 1,6,14 ). In this context, there is great appeal in constructing a new formulation for the PPE which properly accounts for viscous effects (including the non-Newtonian behavior typically observed in blood), while allowing the use of standard C 0 finite element spaces for both the (unknown) pressure and the (given) velocities.
A major breakthrough toward allowing the use of C 0 finite elements for the PPE came from Johnston and Liu, 12 who reformulated the second-order viscous term in the variational formulation as a first-order boundary integral. This approach has since been used successfully in decoupling pressure and velocity for efficient transient Navier-Stokes solvers, 13,[17][18][19] but never in the context of pressure retrieval from given velocities. Moreover, their formulation does not allow for shear-dependent viscosity, which is a phenomenon observed in hemodynamic and polymeric flows. In the present work, we show how the additional terms induced by such variable viscosities can be handled so as to avoid extra regularity requirements. We further construct a general finite element framework, valid for a wide range of flow regimes, by extending the formulation to compressible flows. When incompressibility is dropped, we maintain standard regularity requirements by constructing a mixed formulation where the density and the velocity divergence are introduced as additional variables. In that case, we show how the density field can also be obtained from the velocity, without the need for invoking thermodynamic equations of state.
The rest of this article is organized as follows. We begin the formulation statement by devising a new pressure Poisson equation from the balance of linear momentum for general fluid flows. After presenting a detailed step-by-step derivation for the variational problems, we comment specifically on the compressible Newtonian case and the incompressible quasi-Newtonian case. Numerical aspects are discussed, along with various examples showcasing the accuracy of the proposed framework.

The pressure Poisson equation
We will focus on stationary flows for simplicity of presentation, but the extension to time-dependent problems is straightforward. The two laws that will be used throughout the article are the balance of linear momentum and the conservation of mass. In strong form they can be stated as, respectively, where b is a given body force, u is the flow velocity, p is the pressure, is the density, and S is the viscous stress tensor. We start with a general setting where both compressibility and variable viscosity are allowed. In this case, the viscous tensor can be written as in which is the dynamic viscosity and I is the d × d identity tensor, d = 2, 3 being the number of spatial dimensions. Thus, but Δu ≡ ∇(∇ ⋅ u) − ∇ × (∇ × u). Hence: For now, the viscosity and density fields are assumed to be known. As we will show later, the density can be determined from the velocity, and so can the viscosity in practical cases. Thus, we can write the momentum equation (1) as Now, let us consider a bounded Lipschitz domain Ω ⊂ R d . We can generate a Poisson equation for the pressure by applying (− ∇ ⋅) to both sides of Equation (5). Similarly, appropriate Neumann boundary conditions for the pressure can be obtained by dotting both sides of Equation (5) by the unit normal vector n. 20,21 This gives us the boundary value problem An additive scaling condition is needed to close the problem. It suffices to prescribe the value of the pressure at any point in the domain, or a global constraint such as ∫ Ω p dΩ = , for some ∈ R. In the following derivations, zero mean pressure will be assumed ( = 0), but in Section 2.4 we comment on appropriate conditions for different situations.
We are now in position to start devising a variational formulation for the PPE. Testing with a function q ∈ H 1 (Ω) and using Green's first formula on both sides leads to with (⋅, ⋅) and ⟨⋅, ⋅⟩ Γ denoting, respectively, the L 2 inner product and duality paring. Substituting the Neumann boundary condition (7) leads to the weak form Note that we have split the right-hand side into four contributions. From a Physics standpoint, all of them could be computed as soon as the density and viscosity are determined. However, in the context of finite element spaces, some of these terms impose numerical challenges that must be addressed carefully, since they directly or indirectly require second-order differentiation. The first term can be computed rather straightforwardly. The second one requires a specific technique in the case of non-Newtonian fluids such as blood, and will be tackled in Section 2.3. The third one only requires special attention for compressible flows, and will be dealt with in Section 2.2. The last term requires attention even for the incompressible Newtonian case, so it will be addressed first. We begin by using integration by parts to write Therefore, the weak form of the PPE simplifies to What we have now is a general variational framework for the PPE, allowing generalized Newtonian behavior and compressibility. In the incompressible Newtonian case, that is, ∇ and ∇ ⋅ u both zero, the variational formulation reads simply: Given b and u sufficiently regular, find p ∈ X such that for all q ∈X Notice that, even without stating the precise regularity requirements on u and b, we can straightforwardly see that the use of standard C 0 finite element spaces for the discretization of u, b, p, and q is allowed, 12 since we are left with only (up to) first-order derivatives on all quantities. Yet, we will show next that a mixed framework is needed when compressibility or shear-dependent viscosity is allowed.

Compressible flows
The compressible case brings two issues, one of mathematical nature and another one from a physical standpoint. The former regards the need for higher-order differentiation of the velocity field, due to the presence of the grad-div term in Equation (10). Handling that within a C 0 framework is possible using a simple projection step, which will be described later. The second issue is related to the computation of the density field, which for compressible flows is a variable, rather than a given parameter. Fortunately, we can use the conservation of mass (2) to obtain in terms of the given velocities. This gives us a hyperbolic problem which is well posed if the density is known at the inlet, that is, one must solve 22 where in is a strictly positive function and Γ in is the portion of the boundary where u ⋅ n < 0. The corresponding variational formulation is: Since the velocity is assumed to be known, we end up with a linear problem, which on the discrete level will translate to a simple linear algebraic system. The density can then be fed into the right-hand side of the PPE. Now we turn our attention to the issue of second-order differentiation in the PPE. In order to solve that, one can simply project the velocity divergence onto a continuous space before inserting it into the PPE. On the discrete level, this corresponds to a scalar mass matrix problem. We finally formulate our mixed problem as: Given b and u sufficiently regular, find (p, , ) ∈ X × Y 2 , with = in on Γ in , such that for all (q, r, w) ∈X × Z 2 , where Z = L 2 (Ω), Y = H 1 (Ω) and the remaining spaces are as in the incompressible case. We have at last a mixed variational formulation that can be seen as three linear problems to be solved in sequence. How to treat the viscosity is a problem-dependent matter. It clearly depends on the density, since = , with being the kinematic viscosity. For gas flows, variations of may be given according to local temperatures, 23 or sometimes even neglected depending on the flow regime. For generalized Newtonian fluids such as blood, there are models that relate the viscosity to the velocity gradient. We shall focus on the latter case, which is especially relevant for biomedical applications.

Generalized Newtonian fluids
There are various types of materials and solid-liquid mixtures that behave as generalized Newtonian fluids. This so-called quasi-Newtonian behavior is observed when the viscosity depends locally on the flow field itself. In hemodynamics, the viscosity is often modeled by a nonlinear dependence on the shear rate: In that case, the PPE would require second-order differentiation of the velocity, since the formulation contains the gradient of the viscosity field. Once again, we can stay in the C 0 FEM framework by using a simple L 2 projection onto H 1 (Ω). Instead of plugging the expression = (∇ s u) directly into the PPE, we treat as an independent unknown, and enforce the rheological law weakly. The projected viscosity computed from the given velocities is then inserted into the PPE. Assuming incompressibility, the mixed problem now reads: Given b and u sufficiently regular, find (p, ) ∈ X × Y such that for all (q, w) ∈X × Z, Once again, the mixed formulation can be seen as two linear problems to be solved in sequence.

Pressure scaling
As previously discussed, the standard setting for the PPE is a pure Neumann problem, and therefore a scaling condition is needed to ensure unique solvability. In compressible flows, the pressure is normally treated as a strictly positive quantity, so that the scaling (p, 1) = 0 is not appropriate. Alternatively to this global scaling, the problem can be determined by enforcing the pressure at one point in the domain, which is an available information in practical situations (e.g., far-field pressure, inlet/outlet pressure). Also, in internal incompressible flows with an open outflow region Γ out , the following scaling usually holds: 24,25 This can be enforced by modifying the left-hand side of the PPE to and simply looking for p ∈ H 1 (Ω). On the discrete level, this additional term will lead to a highly sparse symmetric matrix corresponding to a rank-one perturbation. If the right-hand side of Equation (20) is a known nonzero constant, as in domains with multiple outlets, 26 the corresponding value can be used to shift the solution afterward. We remark that, while global or pointwise constraints may have different effects when solving the full Navier-Stokes system, in the present case the scaling results in a mere shift in the solution, as we are solving a simple Poisson problem.

NUMERICAL EXAMPLES
Standard Lagrangian finite element spaces will be used in the interpolation of the unknowns (pressure, density, velocity divergence, viscosity) and the given velocity field. Optimal convergence rates would require the polynomial degree for the velocity to be one order higher than for the unknowns. In practical cases, however, one cannot always afford to use quadratic interpolation, since the velocity is given/measured only in a reduced set of points. For this reason, we will consider both optimal and suboptimal settings in the numerical examples, always using first-order interpolation for the unknowns. In order to assess the accuracy of our formulation, when considering examples with known solutions we use a normalized L 2 norm: The spatial coordinates will be denoted by (x, y, z). Unless where otherwise stated, the body force b is zero.

Compressible flow in two dimensions
We start by considering a compressible flow accelerating through a straight channel Ω = (0, 3) × (0, 1) due to a body force b = {e x (y(y − 1)) 2 , (4x − 1)(1 − 2y)∕3} ⊤ . We set a constant dynamic viscosity = 1, a mass flow rate per unit width equal to 1/6 and an inlet density (0, y) ≡ 1, to get For the numerical solution, we use triangular elements with linear basis functions for p, , and , and quadratic interpolation for u. The coarsest mesh has six elements, as illustrated in Figure 1, then six levels of uniform refinement are considered. Figure 2 shows the error plot with respect to the mesh size h, where quadratic convergence is verified for the pressure approximation.

Power-law fluid in two-dimensional channel
We now consider a classic benchmark for non-Newtonian fluids: the two-dimensional channel flow of a power-law fluid. [27][28][29] This is a popular model for hemodynamic and polymeric flows, and its rheological law is given by where k is a positive parameter and n < 1 for shear-thinning fluids. In a straight channel Ω = (0, L) × with volumetric flow rate per unit width equal to Q, the velocity and pressure fields for the fully developed flow are given by For this example we use the hemodynamic parameters 30 = 1050 kg/m 3 , k = 0.035 Pa.s 0.6 , and n = 0.6, with Q = 100 mm 2 /s, and L = 3H = 3 mm. We start from a mesh composed of six identical square elements and consider six levels of uniform refinement, this time using bilinear basis functions for all quantities of interest. As shown in Figure 3, the convergence is now linear, as first-order interpolation is used for the given velocity field.

Carreau fluid in cylindrical pipe
Another classic benchmark problem in hemodynamics, now in three dimensions, is the Carreau pipe flow. 31,32 The Carreau model describes the shear-thinning behavior of blood through the law where n < 1 and the remaining constants are positive material parameters. In the cylindrical domain For the numerical study, hemodynamic properties are used again: 30 = 1050 kg/m 3 , ∞ = 3.45 mPa.s, 0 = 56 mPa.s, n = 0.3568, and = 3.313 s. The dimensions are fixed as L = 2R = 2 mm. We set a pressure drop of 4.0 Pa across the length. This time, hexahedral elements with trilinear interpolation are employed for all quantities. The coarsest mesh is illustrated in Figure 4, and three levels of refinement are applied. The results of the convergence study are shown in Table 1. Here, again, we cannot expect quadratic convergence, as the velocity is interpolated with the same polynomial degree as the pressure.

Backward-facing step with Carreau fluid
We  Figure 5, with L 2 = 2L 1 = 20H, s = 0.9423H, and H = 5.2 mm. The idea now is to assess the performance of our method in comparison to two more standard techniques: the "inviscid" PPE approach, that is, neglecting the viscous term in the PPE, and the so-called integrated Stokes estimator (STE). 14 For completeness, we briefly present here the STE formulation, which we adapt to the generalized Newtonian case: Given b and u sufficiently regular, find (w, p) The variable w corresponds to a small "artificial velocity" field that works as an auxiliary variable, that is, it has no use as an output (details on the STE formulation can be found in the literature 1,14 ).
Since no analytical solution is known for the BFS problem, we proceed as done by Bertoglio et al. 14 : solving the full incompressible Navier-Stokes problem (on a sufficiently fine mesh), then using the resulting pressure as our ground truth and the velocity solution to feed the retrieval method (PPE or STE). We tackle the (non-Newtonian) Navier-Stokes problem using stable triangular elements with second-and first-order basis functions for velocity and pressure, respectively. As for the BCs, we enforce a parabolic profile with flow rate Q at the inlet x = 0, a "do-nothing" outflow condition ( ∇u)n − pn = 0 at the outlet x = L 1 + L 2 , and no-slip (u = 0) at the remainder of the boundary. The Reynolds number Re = Q∕ is set at 200.
Piecewise linear interpolation has been considered for the pressure; for the STE method, piecewise quadratic functions are used for the artificial velocity w. We compare the different methods by plotting the pressure along the line y = s/2, which includes a region of complex recirculating flow-where convective effects play a major role-and a region of nearly developed flow dominated by viscous effects. The results, scaled so that p = 0 at (x, y) = (L 1 , s/2), are depicted in Figure 6. It is seen that both our novel PPE formulation and the STE approach (appropriately modified for the quasi-Newtonian case) yield virtually identical results which are in very good agreement with the pressure obtained by solving the full Navier-Stokes problem. On the other hand, the classical inviscid PPE approach, while roughly capturing the correct behavior in the vicinity of the step, completely fails to approximate the downstream pressure field.
We remark that, although very accurate, the STE method is considerably more expensive than the PPE approach due to the additional vector-valued unknown-especially since, for stability reasons, this auxiliary field requires a higher order than that used for the pressure. 1,14 Comparing both methods in the presence of noisy velocity measurements can also reveal important aspects, but is out of the scope of our present investigation.

CONCLUDING REMARKS
In the present work, we have devised a general mixed framework for the pressure Poisson equation allowing C 0 finite element interpolation of all (given and unknown) quantities. While classical finite element formulations of the PPE require viscous terms to be dropped, ours is able to account even for the nonlinear viscous behavior of generalized Newtonian fluids, as well as for compressibility. This is accomplished through the use of appropriate vector calculus identities and simple projections of specific flow quantities to circumvent C 1 regularity requirements. We provide numerical examples considering different types of elements, in order to showcase the accuracy of our method in different flow scenarios. While the focus of the present work has been placed on fundamental aspects, the formulation can be readily used in practical applications such as arterial pressure estimation from measured blood flow velocities; the addition of the time-dependent term brings no extra difficulties. Furthermore, the variational formulation developed here can be used in the design of accurate flow solvers, either by replacing the continuity equation or appropriately modifying it for numerical stabilization. This is ongoing work to be covered in forthcoming publications. Also ongoing is the development of an ultra-weak variational framework allowing for pressure discontinuities.