2025-11-22T09:58:16.048981

Superconvergent and Divergence-Free Finite Element Methods for Stokes Equation

Chen, Huang, Zhang et al.
Superconvergent and divergence-free finite element methods for the Stokes equation are developed. The velocity and pressure are discretized using $H(\mathrm{div})$-conforming vector elements and discontinuous piecewise polynomials. The discrete formulation employs a weak deviatoric gradient operator built with tangential-normal continuous finite elements for traceless tensors, requiring no stabilization. Optimal and superconvergent error estimates are established. The method connects to nonconforming virtual element and pseudostress-velocity-pressure mixed formulations. Numerical experiments verify the theory.
academic

Superconvergent and Divergence-Free Finite Element Methods for Stokes Equation

Basic Information

  • Paper ID: 2510.14192
  • Title: Superconvergent and Divergence-Free Finite Element Methods for Stokes Equation
  • Authors: Long Chen, Xuehai Huang, Chao Zhang, Xinyue Zhao
  • Classification: math.NA cs.NA
  • Submission Date: October 16, 2025 to arXiv
  • Paper Link: https://arxiv.org/abs/2510.14192

Abstract

This paper develops superconvergent and divergence-free finite element methods for solving the Stokes equation. The velocity field is discretized using H(div)-conforming vector elements, while the pressure field is discretized using discontinuous piecewise polynomials. The discrete scheme employs a weak deviatoric gradient operator constructed from tangential-normal continuous finite elements based on traceless tensors, without requiring stabilization. Optimal and superconvergent error estimates are established. The method is related to nonconforming virtual element methods and pseudo-stress-velocity-pressure mixed formulations. Numerical experiments validate the theoretical results.

Research Background and Motivation

Problem Significance

The Stokes equation is a fundamental equation describing the motion of incompressible fluids and occupies a central position in computational fluid dynamics. Numerical solution of this equation faces two main challenges:

  1. Divergence-free constraint: The velocity field must satisfy div u = 0, which represents mass conservation
  2. Pressure robustness: Numerical methods should handle the pressure term robustly

Limitations of Existing Methods

Classical stable finite element pairs (such as Taylor-Hood elements, MINI elements, and nonconforming P1-P0 elements) suffer from the following issues:

  • Non-pointwise divergence-free: Cannot strictly satisfy the divergence-free constraint, weakening mass conservation
  • Insufficient pressure robustness: Lack pressure robustness in error estimates
  • Limited convergence order: Some methods suffer from convergence order loss

Research Motivation

This paper aims to construct a finite element method that maintains strict divergence-free properties while achieving superconvergence, simultaneously avoiding the stabilization treatment required by traditional methods.

Core Contributions

  1. Proposed a new mixed finite element scheme: Using H(div)-conforming elements for velocity discretization and discontinuous piecewise polynomials for pressure
  2. Constructed a weak deviatoric gradient operator: Based on tangential-normal continuous finite elements using traceless tensors, without stabilization
  3. Established superconvergent error estimates: Proved superconvergence of order hk+1h^{k+1}
  4. Provided post-processing techniques: Constructed superconvergent velocity approximations
  5. Established equivalence relationships: Connected with nonconforming virtual element methods and pseudo-stress formulations

Method Details

Problem Formulation

Consider the Stokes equation on a bounded domain Ω⊂Rd\Omega \subset \mathbb{R}^d:

-\Delta u - \nabla p = f & \text{in } \Omega \\ \text{div } u = 0 & \text{in } \Omega \\ u = 0 & \text{on } \partial\Omega \end{cases}$$ where $u$ is the velocity field, $p$ is the pressure field, and $f$ is the external force. ### Finite Element Space Construction #### Velocity Space H(div)-conforming Raviart-Thomas (RT) or Brezzi-Douglas-Marini (BDM) elements are employed: $$\stackrel{\circ}{V}^{\text{div}}_{k,\ell} := \{v_h \in H_0(\text{div},\Omega) : v_h|_T \in P_k(T;\mathbb{R}^d) + H_\ell(T)x\}$$ where $\ell = k$ (RT elements) or $\ell = k-1$ (BDM elements). #### Stress Space A traceless tensor space $\Sigma^{tn}_k$ is introduced for discretizing the deviatoric gradient: $$\Sigma^{tn}_k := \{\tau_h \in \Sigma^{-1}_k(\mathcal{T}) : [\Pi_F \tau n]_F = 0 \text{ for all } F \in \mathcal{F}_h\}$$ #### Lagrange Multiplier Space To relax tangential-normal continuity, the following is introduced: $$\Lambda_k = P_k(\mathcal{F}_h;\mathbb{R}^{d-1}), \quad \stackrel{\circ}{\Lambda}_k = P_k(\stackrel{\circ}{\mathcal{F}}_h;\mathbb{R}^{d-1})$$ ### Weak Deviatoric Gradient Operator Define the weak deviatoric gradient operator $\text{dev grad}_w : \stackrel{\circ}{V}^{\text{div}}_{k,\ell} \times \stackrel{\circ}{\Lambda}_k \to \Sigma^{-1}_k(\mathcal{T})$: For $(v,\mu) \in H^1(\mathcal{T}_h;\mathbb{R}^d) \times L^2(\mathcal{F}_h;\mathbb{R}^{d-1})$, defined element-wise: $$(\text{dev grad}_w(v,\mu), \tau)_T = -(v, \text{div } \tau)_T + (n \cdot v, n^\top \tau n)_{\partial T} + (\mu, \Pi_F \tau n)_{\partial T}$$ ### Mixed Finite Element Scheme Find $u_h \in \stackrel{\circ}{V}^{\text{div}}_{k,\ell}$, $\lambda_h \in \stackrel{\circ}{\Lambda}_k$, $p_h \in P_\ell(\mathcal{T}_h)/\mathbb{R}$, such that: $$(\text{dev grad}_w(u_h,\lambda_h), \text{dev grad}_w(v_h,\mu_h)) + (\text{div } v_h, p_h) = (f, v_h)$$ $$(\text{div } u_h, q_h) = 0$$ for all $v_h \in \stackrel{\circ}{V}^{\text{div}}_{k,\ell}$, $\mu_h \in \stackrel{\circ}{\Lambda}_k$, $q_h \in P_\ell(\mathcal{T}_h)/\mathbb{R}$. ### Technical Innovations 1. **Stabilization-free design**: Carefully constructed weak operators avoid penalty terms required by traditional DG methods 2. **Superconvergence**: Achieves convergence rate one order higher than standard methods through commutativity properties 3. **Pressure robustness**: The method naturally possesses pressure robustness 4. **Low-order availability**: Supports low-order cases such as $(k,\ell) = (0,-1), (0,0), (1,0)$ ## Experimental Setup ### Numerical Examples **Example 5.1 (2D)**: - Exact solution: $u = \text{curl } \psi_2$, $p = -x^5 - y^5 + \frac{1}{3}$ - where $\psi_2 = x^2(x-1)^2y^2(y-1)^2$ **Example 5.2 (3D)**: - Exact solution: $u = \text{curl}(\psi_3, \psi_3, \psi_3)^T$, $p = -x^5 - y^5 - z^5 + \frac{1}{2}$ - where $\psi_3 = x^2(x-1)^2y^2(y-1)^2z^2(z-1)^2$ ### Computational Domain and Mesh - Computational domain: $\Omega = (0,1)^d$, $d = 2,3$ - Mesh: Uniform simplicial mesh partitioning - Implementation: Based on MATLAB package iFEM ### Evaluation Metrics - Velocity error: $\|u - u_h\|$ - Stress error: $\|\sigma - \sigma_h\|$ - Pressure error: $\|p - p_h\|$ - Post-processed error: $\|u - u_h^*\|$, $\|\text{grad}_h(u - u_h^*)\|$ ## Experimental Results ### Main Results **Table 1: 2D Error Results** | h | (k,ℓ) | $\|u-u_h\|$ | order | $\|\sigma-\sigma_h\|$ | order | $\|p-p_h\|$ | order | |---|-------|-------------|-------|-------------------|-------|-------------|-------| | 2⁻³ | (0,0) | 2.988e-03 | - | 3.103e-02 | - | 7.810e-02 | - | | 2⁻⁴ | (0,0) | 1.284e-03 | 1.22 | 1.677e-02 | 0.89 | 3.914e-02 | 1.00 | | 2⁻⁵ | (0,0) | 5.988e-04 | 1.10 | 8.700e-03 | 0.95 | 1.963e-02 | 1.00 | | 2⁻³ | (1,0) | 3.296e-04 | - | 2.447e-03 | - | 7.453e-02 | - | | 2⁻⁴ | (1,0) | 8.382e-05 | 1.98 | 6.305e-04 | 1.96 | 3.760e-02 | 0.99 | | 2⁻⁵ | (1,0) | 2.104e-05 | 1.99 | 1.597e-04 | 1.98 | 1.880e-02 | 1.00 | ### Convergence Order Verification Experimental results completely validate theoretical predictions: - **Velocity and stress**: $\|u - u_h\| = \|\sigma - \sigma_h\| = O(h^{k+1})$ (superconvergence) - **Pressure**: $\|p - p_h\| = O(h)$ - **Post-processed velocity**: $\|u - u_h^*\| = O(h^{k+2})$, $\|\text{grad}_h(u - u_h^*)\| = O(h^{k+1})$ ### 3D Results Three-dimensional experiments similarly validate the method's effectiveness, with convergence orders consistent with theory. ## Theoretical Analysis ### Stability Analysis Establishes weak divergence stability condition: $$\inf_{v_h \in \stackrel{\circ}{V}^{\text{div}}_{k,k-1}} \sup_{\tau_h \in \Sigma^{tn}_k} \frac{(\text{div}_w \tau_h, v_h)_{0,h}}{\|\tau_h\|_{\text{div}_w} \|v_h\|} = \alpha > 0$$ ### Error Estimates **Theorem**: Assume $u \in H^{k+2}(\Omega;\mathbb{R}^d)$, then: $$\|\sigma - \sigma_h\|_{0,h} + \|\text{dev grad}_w(I^{\text{div}}_{k,k}u - u_h, Q_{k,\mathcal{F}_h}\lambda - \lambda_h)\| + \|Q_\ell p - p_h\| \lesssim h^{k+1}|u|_{k+2}$$ ### Commutativity Property The key commutativity property: $$Q^{tn}_k \text{dev grad} = \text{dev grad}_w I^{\text{div}}_{k,k}$$ This is the core for achieving superconvergence. ## Related Work ### Divergence-Free Finite Element Methods - **Scott-Vogelius elements**: Require special mesh conditions - **Smooth finite element pairs**: Require super-smooth degrees of freedom - **Conforming pairs on split meshes**: Complex implementation ### Mixed Methods - **MCS method**: Requires $k = \ell \geq 1$; this paper covers more low-order cases - **HDG method**: Requires stabilization; this paper requires none - **Virtual element methods**: Usually require stabilization; lower convergence order ## Conclusions and Discussion ### Main Conclusions 1. Successfully constructed stabilization-free superconvergent divergence-free finite element methods 2. Achieved $h^{k+1}$ order superconvergence, independent of mesh symmetry 3. Method possesses pressure robustness and supports low-order cases 4. Established equivalence relationships with virtual element methods and pseudo-stress formulations ### Limitations 1. Theoretical analysis primarily targets polyhedral domains 2. $H^2$ regularity assumption limits applicability 3. Implementation complexity is relatively high ### Future Directions 1. Extension to Navier-Stokes equations 2. Adaptive mesh refinement 3. Parallel algorithm development 4. Engineering application verification ## In-Depth Evaluation ### Strengths 1. **Rigorous theory**: Complete stability and convergence analysis 2. **Method innovation**: Clever design of weak deviatoric gradient operator 3. **Practical value**: Supports low-order elements and achieves superconvergence 4. **Sufficient experiments**: Complete 2D and 3D numerical verification ### Weaknesses 1. **Implementation complexity**: Computation of weak operators is relatively complex 2. **Theoretical limitations**: Requires strong regularity assumptions 3. **Application verification**: Lacks verification on practical engineering problems ### Impact This work has significant theoretical value in the field of numerical methods for the Stokes equation, providing new insights for constructing high-precision divergence-free methods. The superconvergence property and stabilization-free nature of the method have potential application value in computational fluid dynamics. ### Applicable Scenarios - Fluid computations requiring precise mass conservation - Applications with high computational accuracy requirements - Academic research and method validation ## References The paper cites 56 related references, comprehensively covering important works in finite element methods, mixed methods, virtual element methods, and other directions, with comprehensive literature review.