A High-Order, Pressure-Robust, and Decoupled Finite Difference Method for the Stokes Problem

In this paper, we consider the Stokes problem with Dirichlet boundary conditions and the constant kinematic viscosity $\nu$ in an axis-aligned domain $\Omega$. We decouple the velocity $\bm u$ and pressure $p$ by deriving a novel biharmonic equation in $\Omega$ and third-order boundary conditions on $\partial\Omega$. In contrast to the fourth-order streamfunction approach, our formulation does not require $\Omega$ to be simply connected. For smooth velocity fields $\bm u$ in two dimensions, we explicitly construct a finite difference method (FDM) with sixth-order consistency to approximate $\bm u$ at all relevant grid points: interior points, boundary side points, and boundary corner points. The resulting scheme yields two linear systems $A_1u^{(1)}_h=b_1$ and $A_2u^{(2)}_h=b_2$, where $A_1,A_2$ are constant matrices, and $b_1,b_2$ are independent of the pressure $p$ and the kinematic viscosity $\nu$. Thus, the proposed method is pressure- and viscosity-robust. To accommodate velocity fields with less regularity, we modify the FDM by removing singular terms in the right-hand side vectors. Once the discrete velocity is computed, we apply a sixth-order finite difference operator to approximate the pressure gradient locally, without solving any additional linear systems. In our numerical experiments, we test both smooth and non-smooth solutions $(\bm u,p)$ in a square domain, a triply connected domain, and an $L$-shaped domain in two dimensions. The results confirm sixth-order convergence of the velocity and pressure gradient in the $\ell_\infty$-norm for smooth solutions. For non-smooth velocity fields, our method achieves the expected lower-order convergence. Moreover, the observed velocity error $\|{\bm u}_h-\bm u\|_{\infty}$ is independent of the pressure $p$ and viscosity $\nu$.