Seepage analysis using a polygonal cell-based smoothed finite element method

Abstract: This work develops a polygonal cell-based smoothed finite element method for steady-state, transient, and free-surface seepage in saturated porous media. Wachspress interpolation on convex polygonal elements is combined with cell-based gradient smoothing, so that element matrices are assembled using boundary integrals without in-element derivatives. Polygonal, quadtree, and hybrid quadtree--polygonal meshes are employed to accommodate local refinement and hanging nodes, and a solution-driven adaptive strategy further concentrates resolution near steep gradients and wet--dry transitions. Free-surface seepage is solved using a fixed-mesh iterative scheme that updates the wetted region, permeability field, and boundary conditions. Benchmark tests demonstrate accurate hydraulic-head and free-surface predictions, and show that adaptivity attains similar accuracy with substantially fewer degrees of freedom and CPU time.

• The paper presents a new finite element method combining Wachspress basis functions and gradient smoothing for more accurate and efficient seepage analysis in saturated porous media.
• The method shows significant computational savings (up to 28x) in CPU time and degrees of freedom (up to 7x) for free-surface seepage problems-using adaptive mesh refinement, and improved accuracy compared to conventional finite element methods ( error reduction of up to 4x).
• The novel approach effectively addresses issues in conventional finite element methods, such as mesh distortion sensitivity and the resolution of steep hydraulic gradients, making it highly valuable for complex seepage problems.

Overview and motivation

This paper develops a polygonal cell-based smoothed finite element method (CSFEM) for two-dimensional seepage in saturated porous media, covering steady-state, transient, and free-surface problems. The formulation combines Wachspress rational basis functions on convex polygonal elements with cell-based gradient smoothing, so that all element matrices are assembled from boundary integrals without evaluating shape-function derivatives inside the element. This removes the Jacobian sensitivity of standard isoparametric FEM and, crucially for adaptive analysis, allows hanging nodes from quadtree refinement to be absorbed directly into hybrid polygonal elements without constraint equations or transition elements.

The work addresses three coupled limitations of conventional FEM for seepage: the accuracy/meshing trade-off between triangular and quadrilateral meshes, distortion sensitivity of isoparametric mappings, and the difficulty of resolving steep hydraulic gradients and wet–dry interfaces in free-surface problems. Compared with polygonal scaled-boundary formulations, weak Galerkin, mimetic finite difference, and virtual element methods previously applied to Darcy flow on polygonal meshes, the present approach retains the standard FEM assembly workflow while inheriting the robustness benefits of gradient smoothing.

Formulation

The governing equation is the transient Darcy flow equation $$\nabla \cdot (\mathbf{k}\nabla h) + p - S_s \partial h/\partial t = 0$$ with prescribed-head and prescribed-flux boundaries. Each polygonal element is subdivided into non-overlapping smoothing cells; within each cell, the smoothed gradient is defined as an area average that, via the divergence theorem, reduces to a boundary integral:

$$\widetilde{\nabla} h = \frac{1}{A_C} \int_{\Gamma_C} h(\mathbf{x})\,\mathbf{n}(\mathbf{x})\, d\Gamma .$$

Consequently, the smoothed strain matrix involves only integrals of $N_I$ along straight edges (evaluated with one-dimensional Gauss quadrature), never $\nabla N_I$. Since Wachspress functions are rational but evaluated only along element edges, low-order quadrature suffices. The smoothed stiffness matrix is assembled per smoothing cell as $$A_C \widetilde{\mathbf{B}}^{(C)T} \mathbf{k} \widetilde{\mathbf{B}}^{(C)}$$, while the capacity matrix uses standard interpolation. Time integration uses fully implicit backward Euler.

Wachspress weights are constructed from signed distances to adjacent edge normals, yielding nonnegative, $C^0$, linearly complete, partition-of-unity basis functions on convex polygons. A limitation worth noting: the interpolation theory as presented applies to convex polygons, so nonconvex element support claimed in the abstract is not substantiated by the formulation section.

Free-surface strategy

Free-surface seepage is treated with a fixed-mesh iterative scheme: the dry region is assigned a reduced permeability $\alpha k$ with $\alpha = 10^{-3}$, the pressure head $\psi = h - z$ is computed after each solve, the saturated region is updated where $\psi \ge 0$, and the seepage face on the downstream boundary is reclassified accordingly. Convergence is monitored through the displacement of the overflow point. An optional outer loop performs solution-driven quadtree refinement in a narrow band around $$\widetilde{\nabla} h = \frac{1}{A_C} \int_{\Gamma_C} h(\mathbf{x})\,\mathbf{n}(\mathbf{x})\, d\Gamma .$$0 or large $$\widetilde{\nabla} h = \frac{1}{A_C} \int_{\Gamma_C} h(\mathbf{x})\,\mathbf{n}(\mathbf{x})\, d\Gamma .$$1, with 2:1 balance, transferring the converged head field to the refined mesh between loops. The hybrid quadtree–polygonal discretization is what makes this loop inexpensive: hanging nodes require no additional constraints.

Numerical results

Patch test. Both polygonal and quadtree meshes reproduce the linear hydraulic-head field with relative errors of $$\widetilde{\nabla} h = \frac{1}{A_C} \int_{\Gamma_C} h(\mathbf{x})\,\mathbf{n}(\mathbf{x})\, d\Gamma .$$2 and $$\widetilde{\nabla} h = \frac{1}{A_C} \int_{\Gamma_C} h(\mathbf{x})\,\mathbf{n}(\mathbf{x})\, d\Gamma .$$3 respectively, confirming linear completeness independent of mesh topology.

Steady-state benchmarks. For seepage under an impervious concrete dam, the CSFEM attains a relative error of $$\widetilde{\nabla} h = \frac{1}{A_C} \int_{\Gamma_C} h(\mathbf{x})\,\mathbf{n}(\mathbf{x})\, d\Gamma .$$4 at monitoring points versus $$\widetilde{\nabla} h = \frac{1}{A_C} \int_{\Gamma_C} h(\mathbf{x})\,\mathbf{n}(\mathbf{x})\, d\Gamma .$$5 for conventional FEM (CPE4P elements in ABAQUS) on the same 5 m mesh. For a permeable medium with an impermeable inclusion discretized on a hybrid quadtree mesh, the CSFEM error is $$\widetilde{\nabla} h = \frac{1}{A_C} \int_{\Gamma_C} h(\mathbf{x})\,\mathbf{n}(\mathbf{x})\, d\Gamma .$$6 versus $$\widetilde{\nabla} h = \frac{1}{A_C} \int_{\Gamma_C} h(\mathbf{x})\,\mathbf{n}(\mathbf{x})\, d\Gamma .$$7 for FEM on a triangular mesh — a modest but consistent accuracy advantage.

Transient benchmarks. Problems involving an irregular dam foundation and a square domain with a Stanford bunny-shaped cavity show close agreement with reference solutions in both temporal histories at monitoring points and contour fields, using uniform time steps of 10 s over 2000 s.

Free-surface seepage. For a homogeneous rectangular dam with analytical overflow-point coordinate 0.662382 m, the CSFEM yields relative errors of $$\widetilde{\nabla} h = \frac{1}{A_C} \int_{\Gamma_C} h(\mathbf{x})\,\mathbf{n}(\mathbf{x})\, d\Gamma .$$8 (polygonal), $$\widetilde{\nabla} h = \frac{1}{A_C} \int_{\Gamma_C} h(\mathbf{x})\,\mathbf{n}(\mathbf{x})\, d\Gamma .$$9 (quadtree), and $N_I$0 (adaptive). The efficiency gains are the strongest quantitative claim in the paper: relative to the uniformly refined polygonal mesh (2928 elements, 5808 DOFs, 14.8 s CPU), the adaptive mesh requires only 767 elements, 853 DOFs, and 2.0 s — roughly a 7× reduction in DOFs and CPU time at essentially unchanged free-surface accuracy. For the trapezoidal dam, the contrast is starker: 6809 elements/72 s (polygonal) versus 767 elements/2.6 s (adaptive), a ~28× runtime reduction. These figures imply that adaptivity, rather than the polygonal capability itself, is the dominant source of computational savings in this framework.

Limitations and open questions

Several limitations are acknowledged or implicit. First, the adaptive refinement indicator is heuristic — a narrow band around $N_I$1 or large gradient norm — rather than a rigorous a posteriori error estimator, so no convergence guarantee accompanies the adaptive loop. Second, the fixed-mesh permeability-contrast scheme depends on the choice $N_I$2; the paper does not study sensitivity to this parameter or to the iteration tolerance $N_I$3. Third, the Wachspress interpolation is formulated for convex polygons, leaving the treatment of genuinely nonconvex elements unresolved. Fourth, all validation is two-dimensional; extension to three dimensions, heterogeneous/anisotropic conductivity fields, and rigorous error indicators are explicitly deferred to future work, as are coupled hydro-mechanical applications.

Conclusion

The paper presents a coherent integration of Wachspress-based polygonal interpolation, cell-based gradient smoothing, hanging-node-tolerant hybrid quadtree meshes, and a fixed-mesh free-surface iteration with optional adaptive refinement. Verification against patch tests, analytical solutions, and ABAQUS references supports the central claims: the method passes the patch test to machine-level accuracy ($N_I$4), matches or slightly exceeds conventional FEM accuracy on comparable meshes, and — most significantly — achieves equivalent free-surface accuracy with order-of-magnitude reductions in degrees of freedom and CPU time when adaptivity is enabled. The main open questions concern rigorous a posteriori error estimation, parameter sensitivity of the dry-region permeability penalty, and extension beyond two-dimensional convex-polygon settings.