Stream function -- pressure virtual element methods for the Stokes--Darcy interface problem

Abstract: This paper introduces a novel Virtual Element Method (VEM) for the coupled Stokes--Darcy system in primal-primal form. In the free-flow Stokes domain, we implement a stream function formulation that inherently satisfies the incompressibility constraint and reduces computational cost. Across the interface, mass conservation, normal stress balance, and the Beavers--Joseph--Saffman slip condition are enforced to couple the biharmonic stream function equation with the Darcy's pressure equation. Leveraging VEM's ability to handle general polygonal meshes, the proposed method naturally accommodates irregular interface geometries without requiring remeshing or adaptive refinement. The accuracy of the method is validated through several numerical simulations that include applications to dead-end filtration, and network flow in bioartificial organs.

• The paper proposes a stream function–pressure VEM approach that eliminates the need for inf–sup stability conditions in the Stokes–Darcy system.
• It employs C¹ and C⁰ virtual elements to achieve optimal convergence rates on arbitrary polygonal meshes and nonmatching interfaces.
• Numerical experiments validate the method's robustness and flexibility in complex geometries, highlighting its applicability in biomedical and engineering contexts.

Stream Function–Pressure Virtual Element Methods for the Stokes–Darcy Interface Problem: Technical Summary

The paper "Stream function -- pressure virtual element methods for the Stokes--Darcy interface problem" (2607.01622) introduces and analyzes a novel virtual element discretization for the coupled Stokes–Darcy system. This interface problem couples viscous free flow (Stokes) and porous media flow (Darcy) via interface conditions of mass conservation, normal stress balance, and the Beavers–Joseph–Saffman slip law. The proposed approach leverages a stream function formulation in the Stokes region, eliminating the need for inf–sup stability and enabling the use of $C^1$-conforming VEM on arbitrary meshes, coupled with a $C^0$ VEM discretization for the Darcy pressure.

Model Problem and Stream Function–Pressure Formulation

The Stokes–Darcy system is partitioned into non-overlapping domains: $\Omega_S$ (free flow) and $\Omega_D" (porous media), interfaced through$\Sigma$$. Classical formulations require mixed finite element spaces adhering to inf–sup compatibility, complicating implementation on general meshes. Instead, the stream function$$\chi$is introduced for$\Omega_S$$, satisfying incompressibility by construction, and coupled to Darcy pressure$$\varphi$$via interface bilinear forms derived from physical transmission conditions.</p> <p>The resulting weak form is a coupled system: a fourth-order (biharmonic) variational problem for$$\chi$in$\Omega_S$$and a second-order (diffusion) variational problem for$$\varphi$in$C^0$$0. Importantly, this reformulation ensures all interface transmission conditions—mass conservation, stress balance, and slip—are imposed directly at the primal variable level. <img src="https://emergentmind-storage-cdn-c7atfsgud9cecchk.z01.azurefd.net/paper-images/2607-01622/domain_short.png" alt="Figure 1" title="" class="markdown-image" loading="lazy"></p> <p><img src="https://emergentmind-storage-cdn-c7atfsgud9cecchk.z01.azurefd.net/paper-images/2607-01622/domain_mixed_short.png" alt="Figure 1" title="" class="markdown-image" loading="lazy"> <p class="figure-caption">Figure 1: Sketch of the computational domain and possible boundary configuration.</p> <img src="https://emergentmind-storage-cdn-c7atfsgud9cecchk.z01.azurefd.net/paper-images/2607-01622/domain_darcy_inside_stokes.png" alt="Figure 2" title="" class="markdown-image" loading="lazy"></p> <p><img src="https://emergentmind-storage-cdn-c7atfsgud9cecchk.z01.azurefd.net/paper-images/2607-01622/domain_stokes_inside_darcy.png" alt="Figure 2" title="" class="markdown-image" loading="lazy"> <p class="figure-caption">Figure 2: Schematic representations of the two nested Stokes–Darcy domain configurations.</p></p> <h3 class='paper-heading' id='virtual-element-discretization-methodology'>Virtual Element Discretization Methodology</h3><h4 class='paper-heading' id='polynomial-projections-and-local-spaces'>Polynomial Projections and Local Spaces</h4> <p>The proposed discretization reconciles nonmatching polygonal (or polyhedral) meshes across domains. In$$C^0$1, the $C^0$2 VEM is employed—yielding globally $C^0$3-conforming, locally quadratic shape functions for stream function. In $C^0$4, the $C^0$5 VEM delivers globally $C^0$$6-conforming linear pressure approximations. Both VEM spaces are based on carefully chosen sets of degrees of freedom: stream function and its gradient at vertices, and pressure at vertices. <img src="https://emergentmind-storage-cdn-c7atfsgud9cecchk.z01.azurefd.net/paper-images/2607-01622/sketch_class.png" alt="Figure 3" title="" class="markdown-image" loading="lazy"> <p class="figure-caption">Figure 3: Schematic overview of the filtering procedure for interface DoFs in the VEM framework.</p></p> <p>This mesh treatment naturally enables intricate geometries and robust handling of nonmatching meshes at the interface, facilitated by an implementation in the VEM++ object-oriented C++ library.</p> <h4 class='paper-heading' id='discrete-formulation'>Discrete Formulation</h4> <p>The discrete bilinear forms employ projections onto polynomial spaces to ensure consistency and stability. Stabilization operators, crucial for computability of non-polynomial virtual spaces, are constructed at the element level and weighted according to physical parameters.</p> <p>The global system is a perturbed saddle-point problem, but, due to the primal–primal structure and stream function formulation, the resulting discretization is free from inf–sup constraints.</p> <h3 class='paper-heading' id='a-priori-error-analysis'>A Priori Error Analysis</h3> <p>The paper provides rigorous energy-norm error bounds for the VEM scheme. The interpolation theory for$$C^0$7 and $C^0$8 VEM is combined with trace inequalities and polynomial approximation properties to deliver optimal convergence rates.

Formally, for local mesh size $C^0$9, and exact solutions $\Omega_S$0 with additional regularity, the total error obeys

$\Omega_S$1

for $\Omega_S$2, with constants $\Omega_S$3 only depending on the domain, physical parameters, and mesh regularity. The stabilization is shown to be robust for all reasonable parameter regimes.

Numerical Experiments

Comprehensive numerical tests are performed to validate accuracy and flexibility. Results confirm optimal linear rates for both structured and unstructured meshes, including challenging cases with nonmatching interfaces and highly irregular element shapes.

Manufactured Solutions: Convergence Verification

• Manufactured test cases with both zero and nonzero interface values exhibit energy error decay rates consistently matching theory.

Figure 4: Variety of discretisations used for the unit square in Experiments 1 and 2.

Figure 5: Experiment 1. Nodal values of the VEM solutions $\Omega_S$4 and $\Omega_S$5 on the Voronoi mesh at the final refinement level.

Figure 6: Experiment 2. Nodal values of the VEM solutions $\Omega_S$6 and $\Omega_S$7 on the finest refinement of the perturbed Voronoi mesh.

Heterogeneous Domains: Dead-End Filter

• The quarter-annulus dead-end filter benchmark demonstrates the method’s ability to cope with nontrivial geometries and capture local phenomena near the interface, including sharp gradients in the vicinity of the filtration membrane.

  Figure 7: Experiment 3. Domain configuration for quarter annulus dead-end filter discretised with a coarse Voronoi mesh.

  

Figure 8: Experiment 3. Snapshots of the variables of interest for the quarter annulus dead-end filter domain discretised with 15,501 Voronoi elements.

Biomedical Application: Bioartificial Organ Perfusion

• The simulation of blood flow in a scaffolded bioartificial organ with embedded straight channels quantifies the impact of tissue permeability on velocity and pressure fields, showing that increasing permeability redistributes the flow between the vascular network and surrounding medium.

  Figure 9: Experiment 4. Blood flow network with straight channel embedded within a bioartificial scaffold domain discretised with a coarse Voronoi meshes.

  

Figure 10: Stokes' and Darcy's velocities for a straight blood flow network embedded within a bioartificial scaffold domain discretised with a Voronoi mesh for several values of $\Omega_S$8.

Figure 11: Nodal values of the stream function $\Omega_S$9 and Darcy's pressure $\Omega_D" (porous media), interfaced through$0 for a straight blood flow network and varying $\Omega_D" (porous media), interfaced through$1.

Implications, Extensions, and Open Challenges

This work establishes a rigorous, efficient, and flexible framework for simulating coupled Stokes–Darcy problems on general meshes without remeshing or aggressive interface-conforming requirements. The stream function formulation eliminates the need for inf–sup stable pairs, which has significant computational benefits and simplifies implementation, especially for high-order discretizations or adaptive mesh refinement.

Key implications and challenges:

• Practical flexibility: The ability to use arbitrary polygonal meshes makes the method highly applicable to real-world domains with geometric complexity or moving/deforming interfaces.
• No inf–sup constraints: This greatly broadens the potential for codebases and discretizations, as standard mixed FE inf–sup compatibility checks become unnecessary.
• High regularity demand: The need for $\Omega_D" (porous media), interfaced through$2-conforming spaces for the biharmonic portion, while tractable for 2D VEM, is technically intricate for 3D or higher polynomial degrees. Extension to 3D would require vector potential formulations and advanced high-order $\Omega_D" (porous media), interfaced through$3 VEM.
• A posteriori analysis: While this work provides a priori estimates, a robust a posteriori error estimator for adaptive refinement remains open.
• Curved geometries and extreme mesh irregularity: Further work is needed to handle highly curved interfaces and relax the mild mesh regularity assumptions.

Conclusion

The paper delivers an effective VEM-based method for Stokes–Darcy interface problems, combining the stream function–pressure formulation with polygonal discretization to provide robust, inf–sup-free, and optimally convergent numerical solutions. It opens avenues for generalizing to more complex physics (e.g., multiphysics, nonlinearities, or multiple coupled subdomains), higher-order spatial discretizations, and three-dimensional simulations. Future research includes the design of stabilization-free variants, local mesh adaptation strategies, and modular pressure recovery algorithms relevant for practical engineering and biomedical applications.