I am trying to model a problem of a nearly incompressible $ 10~\rm{m} \times 2~\rm{m}$ beam with a uniformly distributed end load. The beam has a Young’s modulus of $ 200~\rm{Pa}$ and a Poisson’s ratio of $ \nu=0.49995$ and I am using a mixed displacement-pressure (Q1P0) formulation with a perturbed Lagrangian approach.

I am using meshes of $ 2^N \times 2^N$ Q1P0 elements and for $ N\leq6$ to solution procedure is working fine. However, for $ N>6$ the procedure diverges very early on. I have tried changing the initial load and the load step size in the solution procedure, and I have tried changing the material properties. Whatever I try I cannot get the solution procedure to work for $ N>6$ , which is a problem because I need solutions for at least up to $ N=8$ .

Note: The solution procedure works fine with standard Q2 elements for $ N=8$ and I can’t figure out why it doesn’t work for the Q1P0 elements.

The solution procedure I am using is given below (unfortunately I don’t know how to make the Greek characters display properly in the code block)

`SMTAnalysis[]; tolNR = 10^-5; maxNR = 500; targetNR = 100; \[Lambda]Max = 1; \[Lambda]0 = \[Lambda]Max/1000; \[CapitalDelta]\[Lambda]Min = \[Lambda]Max/10000; \[CapitalDelta]\[Lambda]Max = \ \[Lambda]Max/100; SMTNextStep["\[Lambda]" -> \[Lambda]0]; While[ While[ step = SMTConvergence[tolNR, maxNR, {"Adaptive BC", targetNR, \[CapitalDelta]\[Lambda]Min, \ \[CapitalDelta]\[Lambda]Max, \[Lambda]Max}] , SMTNewtonIteration[]; ]; If[step[[4]] === "MinBound", SMTStatusReport["Analyze"]; SMTStepBack[];]; step[[3]] , If[step[[1]], SMTStepBack[];]; SMTNextStep["\[CapitalDelta]\[Lambda]" -> step[[2]]] ]; `

Plotting the displacement of the end of the beam against $ N$ yields the following

It is clear that the behavior is good up until $ N=6$ .

The element code and notebook in which problem is solved (including solution procedure) are available here:

Q1P0 element

Convergence Test

Any help in improving the solution procedure would be appreciated!