The structure of the solution to the generalized Stokes equations and a new method for solving them

We show that the solution to the generalized Stokes equations with Dirichlet boundary conditions in dimension two or three has a particular structure. Namely, the velocity $u$ is the superposition of two velocity fields $\omega$ and $\theta$, $u=\omega+\theta$, each of them being divergence-free and solving vector Helmholtz-like equations. The right hand sides of these equations involve the rotational component of the Helmholtz decomposition of the external force, and a harmonic function $q$, called solid pressure. Similarly, the fluid pressure $p$ is the superposition of the potential $\pi$ of Helmholtz decomposition of the external force, and of the solid pressure $q$, $p=\pi+q$. It turns out that regardless the solid pressure $q$, $\omega+\theta$ satisfies the generalized Stokes equations, except the normal component of the Dirichlet boundary condition. The role of the solid pressure, which solves a linear boundary equation with a self-adjoint coercive operator, is to constrain $\omega+\theta$ to satisfy even the normal component of the Dirichlet boundary condition.
The method is attractive from the numerical viewpoint. It decouples the velocity from the pressure and does not solve the incompressibility constraint directly. It requires only the solution of Helmholtz-like vector equations coupled on the boundary, and of a boundary equation. In the last section we present two methods for solving numerically the generalized Stokes equations, and some results to demonstrate the efficiency of our approach.