Julien Vignollet,Chris J. Pearce,Lukasz Kaczmarczyk
Abstract:An hyperelastic biphasic model is presented. For slow-draining problems (permeability less than 1\times10-2 mm4 N-1 s-1), numerical instabilities in the form of non-physical oscillations in the pressure field are observed in 3D problems using tetrahedral Taylor-Hood finite elements. As an alternative to considerable mesh refinement, a Galerkin least-square stabilization framework is proposed. This technique drastically reduces the pressure discrepancies and prevents these oscillations from propagating towards the centre of the medium. The performance and robustness of this technique are demonstrated on a 3D numerical example.
What problem does this paper attempt to address?
This paper attempts to solve the numerical instability problems that occur when using the finite element method to simulate biphasic soft tissues (such as intervertebral discs) with low permeability. Specifically, when using tetrahedral Taylor - Hood elements to simulate three - dimensional problems, non - physical oscillation phenomena will occur in the pressure field, especially when the permeability is less than \(1\times10^{- 2}\, \text{mm}^4\text{N}^{-1}\text{s}^{-1}\). These oscillations will propagate to the center of the medium, resulting in inaccuracies in the numerical solution.
To solve this problem, the author proposes a stabilization technique based on the Galerkin Least - Square (GLS). This technique significantly reduces non - physical oscillations in the pressure field and prevents these oscillations from propagating to the center of the medium by introducing a weighted least - squares term to correct the weak - form equations. The specific formula of the GLS stabilization method is as follows:
\[
R = R_u-\Delta tR_p + R_{\text{GLS}}=0
\]
where \(R_{\text{GLS}}\) is the least - squares term from the strong form of the fluid continuity equation:
\[
R_{\text{GLS}}=\int_V\left[\nabla\cdot(\dot{\mathbf{u}} + k\nabla p)\right]^T\tau^*\left[\nabla\cdot(\mathbf{v}_S + k\nabla p)\right]\,dV
\]
After discretization and simplification, the final linearized matrix form is:
\[
\begin{bmatrix}
K_{uu}+K_{\text{GLS}, uu}&K_{up}\\
K_{up}^T&\Delta tK_{pp}
\end{bmatrix}
\begin{Bmatrix}
\delta\mathbf{u}\\
\delta p
\end{Bmatrix}
=
\begin{Bmatrix}
\mathbf{F}_{\text{ext}, u}\\
\Delta t\mathbf{F}_{\text{ext}, p}
\end{Bmatrix}
-
\begin{Bmatrix}
\mathbf{F}_{\text{int}, u}+\mathbf{F}_{\text{GLS}, \text{int}, u}\\
\Delta t\mathbf{F}_{\text{int}, p}
\end{Bmatrix}
\]
where the stabilization factor \(\tau_{\text{GLS}}\) is defined as:
\[
\tau_{\text{GLS}}=\frac{h^2}{4k\Delta t}
\]
Here, \(h\) is the element characteristic size (usually taken as the radius of the circumscribed sphere of the element), \(k\) is the permeability, and \(\Delta t\) is the time step.
Through numerical experiments, it is verified that the GLS stabilization method can effectively eliminate non - physical pressure oscillations and shows good robustness and stability for different mesh densities, permeabilities, loading rates, compressive strains, and time steps. In addition, this method does not require additional degrees of freedom and has a low computational cost.