A Python-Based Finite Element Digital Twin for Thixoviscoplastic Lid-Driven Cavity Flow: Implementation, Verification, and Validation Against Turek Benchmark Results

Abstract

A Python-based finite element (FEM) solver is developed for steady incompressible flow of thixoviscoplastic (TVP) fluids in a lid-driven square cavity. The solver adopts the Houska TVP constitutive model of Begum, Ouazzi, and Turek [1], in which a scalar microstructure parameter λ ∈ [0, 1] governs structural buildup and viscosity is coupled nonlinearly to both the shear rate γ˙ and λ through a Papanastasiou-regularised extended viscosity. The implementation uses Taylor-Hood Q2/Q1 elements on structured quadrilateral meshes, assembled entirely in Python (NumPy,SciPy, GMSH) with monolithic Newton iteration and a direct sparse-LU inner solve. A five-stage verification and validation pipeline is presented: the Stokes assembler reproduces the exact Poiseuille profile to floating-point precision (L2 error < 10−14) and achieves optimal manufactured-solution convergence rates of 3.3/2.3/2.0 in L2 velocity/H1 velocity/L2 pressure; the Navier-Stokes solver matches the Ghia et al. [8] benchmark at Re = 100 to within 5.0 × 10−3 m/s; power-law and Bingham solvers are verified to relative L2 errors of 7.5 × 10−5 and 2.7 × 10−5; and the monolithic TVP Jacobian passes seven finite-difference checks at O(10−10) relative error. A conditioning study derives empirically that κ(Jbc) ∝ η2 0, reaching 3.9 × 107 at the paper's target parameters, providing a rst-principles justi cation for the multigrid infrastructure of Begum et al. [1]. Newton converges in 6-13 iterations per stage at reduced parameters (η0 = 10−3 Pa). A square obstacle introduced at the cavity centre produces an eight-fold reduction in the high-velocity domain fraction,constituting the paper's original geometric contribution and identifying a second,mesh-geometry-induced source of ill-conditioning on unstructured meshes.