The grad-div conforming virtual element method for the quad-div problem in three dimensions

We propose a new stable variational formulation for the quad-div problem in three dimensions and prove its well-posedness. Using this weak form, we develop and analyze the $\boldsymbol{H}(\operatorname{grad-div})$-conforming virtual element method of arbitrary approximation orders on polyhedral meshes. Three families of $\boldsymbol{H}(\operatorname{grad-div})$-conforming virtual elements are constructed based on the structure of a de Rham sub-complex with enhanced smoothness, resulting in an exact discrete virtual element complex. In the lowest-order case, the simplest element has only one degree of freedom at each vertex and face, respectively. We rigorously prove the interpolation error estimates, the stability of discrete bilinear forms, the well-posedness of discrete formulation and the optimal error estimates. Some numerical examples are shown to verify the theoretical results.

💡 Research Summary

The paper addresses the numerical solution of the three‑dimensional “quad‑div” problem, i.e. the fourth‑order partial differential equation (∇ div)² u = f together with the constraints curl u = 0, u·n = 0 and div u = 0 on the boundary of a contractible Lipschitz polyhedral domain Ω. While the biharmonic and quad‑curl operators have been extensively studied, the quad‑div operator has received comparatively little attention, especially regarding stable mixed formulations and conforming discretizations on general polyhedral meshes.

The authors first develop a new variational formulation for (1.1). By characterising the dual space of the H(grad‑div)‑conforming space V₀(Ω) as H⁻²(curl;Ω) and exploiting a generalized Helmholtz decomposition, they identify the admissible source space as f ∈ H⁻²(Ω)∩ker(curl). The mixed formulation introduces three unknowns: the primary field u ∈ V₀(Ω), a Lagrange multiplier φ ∈ H₀(curl;Ω) enforcing the curl‑free condition, and a pressure‑like scalar p ∈ H₀¹(Ω) to control the kernel of the curl operator. The resulting system (3.3) consists of a coercive bilinear form A((u,p),(v,q)) = (∇ div u, ∇ div v) together with a coupling form B((v,q),φ) = (v, curl φ)+(∇q, φ). Using a Friedrichs inequality for divergence‑free and curl‑free subspaces, the authors prove that A is V‑coercive on the kernel of B and that B satisfies the inf‑sup condition. Consequently, the mixed problem is well‑posed, and the Lagrange multipliers vanish (φ = 0, p = 0), reducing the system to the original quad‑div equation in the distributional sense.

Having established a stable continuous formulation, the paper proceeds to construct a high‑order H(grad‑div)‑conforming virtual element method (VEM) on arbitrary polyhedral meshes. The design is guided by an enhanced de Rham complex with increased smoothness. The authors define a discrete complex (1.5) parameterised by two integers r, k ≥ 1: R → U₁(Ω) → Σ₀,r(Ω) → V_{r‑1,k+1}(Ω) → W_k(Ω) → 0, where U₁(Ω) and W_k(Ω) are subspaces of H¹(Ω), Σ₀,r(Ω) is an H(curl)‑conforming virtual element space, and V_{r‑1,k+1}(Ω) is the target H(grad‑div)‑conforming space. Three families are obtained by setting (r = k), (r = k+1) and (r = k+2). In the lowest‑order case (r = k = 1) the element is remarkably simple: each vertex contributes one degree of freedom (DOF) for the normal component and each face contributes one DOF for the tangential component, yielding a total of |V| + |F| DOFs per element.

The virtual element spaces are defined implicitly via local boundary value problems; the shape functions are not required to be known explicitly. Interpolation operators are constructed so that the following diagram commutes: continuous complex → interpolation → discrete complex, ensuring that the discrete spaces inherit the exactness of the continuous de Rham sequence. The authors prove optimal interpolation error estimates of order hᵏ for functions with sufficient regularity, i.e. ‖v − I_h v‖{H(div)} ≤ C hᵏ‖v‖{H^{k+1}}.

Stability of the discrete bilinear forms is established by defining a suitable stabilization term on each element. Lemma 5.1 provides a discrete Friedrichs inequality, guaranteeing that the discrete stiffness matrix A_h is coercive on the kernel of the discrete B_h and that B_h satisfies a discrete inf‑sup condition. Consequently, the discrete mixed problem admits a unique solution (u_h, φ_h, p_h) and inherits the same a priori bound as the continuous problem.

Error analysis combines the interpolation estimates with the consistency of the discrete forms. The authors obtain optimal convergence rates: ‖u − u_h‖{V} ≤ C hᵏ‖u‖{H^{k+1}},  ‖φ − φ_h‖{H(curl)} ≤ C hᵏ‖φ‖{H^{k+1}},  ‖p − p_h‖{H¹} ≤ C hᵏ‖p‖{H^{k+1}}. Because φ and p vanish analytically, the primary error estimate reduces to the optimal rate for u.

The theoretical results are corroborated by a suite of numerical experiments on both structured (tetrahedral, hexahedral) and unstructured polyhedral meshes, including highly irregular elements. The experiments confirm the predicted convergence orders and demonstrate that even the lowest‑order element, with only vertex and face DOFs, achieves the expected accuracy. Moreover, the method exhibits robustness with respect to mesh distortion, highlighting the flexibility of the VEM framework for complex geometries.

In summary, the paper makes several substantial contributions: (1) a novel, well‑posed mixed variational formulation for the quad‑div problem with curl‑free data; (2) the construction of H(grad‑div)‑conforming virtual element spaces of arbitrary order on general polyhedral meshes, preserving an exact discrete de Rham complex; (3) rigorous proofs of interpolation error, stability, well‑posedness, and optimal convergence; and (4) comprehensive numerical validation. These results fill a notable gap in the numerical analysis of fourth‑order operators and open the door to applying H(grad‑div)‑conforming VEMs in advanced elasticity, fluid‑structure interaction, and other multiphysics contexts where quad‑div operators naturally arise.