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)

tolNR = 10^-5; 
maxNR = 500; 
targetNR = 100;
\[Lambda]Max = 1; 
\[Lambda]0 = \[Lambda]Max/1000; 
\[CapitalDelta]\\[Lambda]Min = \[Lambda]Max/10000; 
\[CapitalDelta]\[Lambda]Max = \

SMTNextStep["\[Lambda]" -> \[Lambda]0];
   step = 
     maxNR, {"Adaptive BC", 
      targetNR, \[CapitalDelta]\[Lambda]Min, \
\[CapitalDelta]\[Lambda]Max, \[Lambda]Max}]
   , SMTNewtonIteration[];
  If[step[[4]] === "MinBound", SMTStatusReport["Analyze"]; 
  , If[step[[1]], SMTStepBack[];];
  SMTNextStep["\[CapitalDelta]\[Lambda]" -> step[[2]]]

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

enter image description here

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!

  • $\begingroup$ Is Young's modulus 200 Pa or 200 GPa? $\endgroup$ Jul 23, 2021 at 10:13
  • $\begingroup$ The Young's modulus is 200Pa. $\endgroup$ Jul 26, 2021 at 7:14

1 Answer 1


I have noticed that your perturbation parameter goes to Infinity when Poisson's ratio goes to 0.5. Singularity always spells trouble for numerical procedures.

You have also formulated Lagrangian at the global level leading to the zeroes at the main diagonal. Additionally, close to singularity global matrix becomes ill-posed. The combined effects resulted in deteriorating convergence characteristics due to the MKL linear solver failure.

My advice is that you use static condensation of the Lagrangian (SMTNoDOFCondense option) to improve convergence.

  • $\begingroup$ This is looks like comment, not like an answer. Could you focus on $\nu=0.49995$, not on $\nu=0.5$? $\endgroup$ Jul 23, 2021 at 10:05
  • $\begingroup$ I have tried with a few different choices of perturbation parameter and no matter what I tried I got the same divergence problem. $\endgroup$ Jul 26, 2021 at 7:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.