On the pressure Poisson equation for the Stokes system

For the numerical solution of the pressure Poisson equation, we consider an ultra–weak variational formulation and a related finite element method of Galerkin–Petrov type. While this allows a piecewise constant approximation of the pressure, the test functions have to be sufficiently smooth – in particular, we use modified B–splines of second order to guarantee conformity and stability, and optimal convergence. This novel formulation allows the computation of the pressure from either a given, measured, or precomputed velocity field as it appears in many applications in fluid mechanics.


Introduction
In biomedical applications, more specifically in hemodynamics, estimating local blood pressure using non-invasive techniques is a critical task which is hardly viable through direct measurements. On the other hand, modern engineering and clinical imaging techniques such as 4D MRI, PIV and echocardiograms allow quantification of blood velocities in arteries with various degrees of accuracy [2]. In that context, the past decades have seen a development of techniques that combine empirical assessment and mathematical modeling through advanced numerical methods.
The most popular approaches are based on the pressure Poisson equation, which is derived from the balance of momentum equation by applying the divergence operator. Although simple in theory, the pressure Poisson imposes practical challenges, as even in weak form it requires the computation of second-order derivatives of the velocity. This prohibits the use of standard piecewise linear finite element functions for the measured velocity field. Furthermore, even if higher-order basis functions are employed, second-order differentiation of noisy velocity measurements can produce inacurate pressure results.
Motivated by those challenges, we consider here an ultra-weak variational formulation for the pressure Poisson equation with a piecewise constant approximation of the pressure, which also allows the use of piecewise bilinear approximations of the velocity. As numerical example we consider a lowest-order discretization in two dimensions.

Problem statement and variational formulation
As a model problem, we consider the stationary Dirichlet problem for the Stokes system in a bounded Lipschitz domain in which u is the flow velocity, p is the pressure, and f is a given body force per unit volume. We obtain the pressure Poisson equation by applying the divergence operator to both sides of the momentum equation in (1), As boundary condition for the pressure Poisson equation we consider the normal trace of the momentum equation [1], Since in the standard variational formulation of the Stokes system (1) we are looking for (u, p) ∈ H 1 0 (Ω) × L 2 0 (Ω), we have to use an ultra-weak variational formulation for the pressure Poisson equation to find p ∈ X := L 2 (Ω) such that is satisfied for all ϕ ∈ Y := {ψ ∈ H 1 (Ω) : ∆ψ ∈ L 2 (Ω), ∂ n ϕ := ∇ϕ · n = 0 on Γ}, where a suitable norm is given by Unique solvability of the variational formulation (4) then follows from the inf-sup stability condition [3] min{1, |Ω|} p L 2 (Ω) ≤ sup 0 =ϕ∈Y a(p, ϕ) ϕ Y for all p ∈ L 2 (Ω). (5)

Discretization
For the numerical solution of the variational formulation (4), we consider a decomposition of Ω into N quadrilateral elements of size h, and use conforming finite element spaces X h = span{ψ k } N k=1 ⊂ X of piecewise constant basis functions ψ k , and Y h = span{ϕ k } N k=1 ⊂ Y of modified B-splines ϕ k of second order, which satisfy the homogeneous Neumann boundary condition ∂ n ϕ k = 0 on Γ. The Galerkin-Petrov variational formulation of (4) is to find p h ∈ X h such that is satisfied for all ϕ h ∈ Y h . Unique solvability of (6) then follows from the discrete counterpart of the inf-sup condition (5), which is satisfied for the chosen finite element spaces [3]. Moreover, using standard arguments we obtain Cea's lemma, i.e., when assuming p ∈ H 1 (Ω). Note that for the evaluation of the right-hand side in (6) we may use a piecewise bilinear approximation u h , which does not change the order of convergence when assuming u ∈ H 2 (Ω).

Numerical example
We consider the Stokes system (1) in a T-shaped domain Ω with an initial mesh of N = 48 quadrilateral elements as depicted in Fig. 1 (left), where f and g are computed from the given solution Starting from the initial mesh, we applied 5 uniform refinement levels, with 49152 elements on the finest level. The error plot in Fig. 1 reveals linear convergence, as expected from the error estimate (7).