Why FEM for incompressible materials is ill-posed?

Question

I am an engineer who is trying to get a deeper understanding of FEM. I have been using the Zienkiewicz texts as my bible. It touches on the issue of incompressibility but I need a more intuitive way to think about what is happening here.

My current understanding

A 1 field displacement-based formulation has convergence issues due to the incompressibility leading to an ill-conditioned matrix. I am still not sure how this arises and whether this hold for all material types. A crude way of resolving this would be to employ higher order elements, which is computationally prohibitive, demanding the need for low order work arounds. Enter two and three field formulations...

Explicitly having additional variables can serve as constraints on the system can help with the ill condition matrix at the expense of solving for more variables. So my understanding is that an issue arises when we are interpolating values in a disconnected manner. I guess that this lack of a connection is due to the overly simple continuity equation that arises from incompressibility. This causes the locking phenomenon which leads to spurious results. I know I have gaps in that flow of logic. I am also confused as to why this resolves the convergence issues. Wouldn't the matrix still be prone to ill-conditioning?

This locking requires a mixed formulation where the new variables require a higher order than the displacement field to be stable. So we need elements that satisfy the LBB or inf-sup condition. I do not have an intuitive feel for what is going on here. What would happen if the displacement field had a higher order than these extra variables?

I have tried looking at the original LBB paper and my eyes glazed over. Hopefully, there is a less formal way of describing what's going on here.

Answer 1

For incompressible materials, as the extent of incompressibility increases, the bulk modulus approaches infinity (for the isotropic case). This is what causes ill-conditioning.

Consider the case of linear isotropic elastic material for which the elasticity tensor is given as \begin{equation} D_{ijkl} = \mu \, (\delta_{ik} \delta_{jl} + \delta_{il} \delta_{jk} - \frac{2}{3} \delta_{ij}\delta_{kl}) + \kappa \, \delta_{ij}\delta_{kl} \end{equation} where $\mu$ is the shear modulus, and $\kappa$ is the bulk modulus.

Since the stiffness matrix for the pure displacement formulation is \begin{equation} \mathbf{K} = \int \mathbf{B}^T \mathbf{D} \mathbf{B} dV \end{equation} it becomes ill-conditioned as $\kappa \rightarrow \infty$. Lower-order elements perform quite poorly, known as the locking phenomenon. Although we can get a better displacement field with higher-order elements, we do not get optimal convergence rates. Also, the stress field contains spurious oscillations, see paper1 and paper2 for further details.

Note: It is a misconception that higher-order elements are computationally expensive. As already pointed out, we do not need finer meshes with higher-order elements, leading to overall reductions in total time for a simulation.

Mixed formulations are what we use to remedy the issue of locking and spurious stress fields. However, the overall stiffness matrix can still be ill-conditioned depending on the particular formulation and how we solve the matrix. For example, the matrix is still ill-conditioned for $\mathbf{B}$-bar (and $\mathbf{F}$-bar) formulations, and when we condense the pressure DOFs out in the displacement-pressure formulation.

Since LBB stability is a purely mathematical construct it is difficult to get an intuitive feel for it, especially if one does not have a background in functional analysis. I suggest referring to the textbooks on FEM by Zienkiewicz et al. The constraint ratio they define helps to understand LBB stability using simpler concepts.

The order of elements for the displacement field should be higher in the mixed displacement-pressure formulation. It is lower only in the mixed displacement-stress formulation. If the basis functions do not satisfy the LBB condition, then such a combination results in spurious fields.