# Finite Element Code: Convergence issues

I am trying to code a simple FEM problem. A bar is fixed at the bottom edge, and a displacement is applied at the top. I want to apply the displacement gradually , so I have divided it into a number of steps. A rough outline of the code looks like this:

 %********************************************************* %
> Newton-Raphson Control parameters
> %*********************************************************
>
> uload=5;         % the total displacement
numu=10;         % number of steps
tu=0;            %initial value
> tol=1e-2;

>   for ii=1:numu  %Loop over total number of steps
>
>       iter=0;
>       tu=tu+du;  %increment displacement at each step
>

>          while(error>=tol)  %Iteration
>            iter=iter+1
>            Ku=sparse(sdof,sdof);    %Initialize stiffness matrix
>            [Ku] = Kmat(Ku);  %Create stiffness matrix by user defined function
>               Ru = sparse(sdof,1);  %Initialize residual force vector
>            [Ru] = Rvec(Ru,u);  %Calculate residual force vector by user defined function(consists of internal force only)
>            [Ku,Ru] = applyBC(Ku,Ru); %Apply BCs by user defined function

>            delu = sparse(sdof,1) ;  %Initialize solution vector
>            delu = Ku\Ru ;   % Calculate displacement
u=u+delu; %   %Increase displacement value
>            error=norm(Ru)  %Calculate Error (This step is probably causing the divergence due to selection of improper criterion)

>          end %iteraion
> end %step


This code won't compile in MATLAB because all the functions are user defined and the input data is not provided. I just wanted to portray a rough outline of what I am trying to do.

This algorithm, works fine if I only have a few elements(suppose 100) defined. If I decrease the mesh size, the solutions tends to diverge, the error keeps on increasing. I figure that is probably because the way I calculate the error is not useful for this simple case.

Now my question is what can be the best convergence criterion for this problem? Since the problem is simple linear elastic, no non-linearity involved, I think the values calculated in the first iteration is correct, so no need for further iterations. But later I would like to extend this code to perform non-linear static analysis, so is there a convergence criterion that can be used for both? If not, then what can be used for each case?

P.S: I have modeled a similar problem in a commercial code, where obviously it works, but when I check the details of the solution, it says that it has used one iteration per step(which is expected). In my code with more elements, the iterations keep increasing after the 2nd step. Thus my doubt regarding the choice of convergence criterion.

• [Ku] = Kmat(Ku) looks odd. Shouldn't Ku depend on the current estimate of the displacements, for a nonlinear problem? Why does it depend on the previous stiffness matrix? Also, it seems a bit strange to calculate the error from the residual, before you update the residual after the solution increment. But there are so many unknowns here, trying to debug this is just guessing. – alephzero Mar 30 '19 at 16:37
• Ah, I just noticed you are only testing this on a linear elastic problem. So it should converge in exactly 1 iteration. If it doesn't, something is wrong! If it converges in "more than one" iteration for small models but doesn't converge for big models, the code is wrong for the small models even though it "works." It all else fails, check out everything for a 2-element model with just one unknown displacement in the middle of the rod. – alephzero Mar 30 '19 at 16:42
• Yes, I think that bug is in the convergence criterion I chose( norm of residual). I tested my code with single step static problems of different kinds, they are an exact match of results from ABAQUS, so the underlying calculations are fine. But I am unable to figure out what criterion should be chosen so that it converges in the first iteration. – Schneider Mar 30 '19 at 17:05
• Given your doubt, I changed my code and re-ran it by calculating the residuals before the error calculation, it converges in 2 iterations for the small problem. In all the cases the error is the same, however, I do get converged results in the end. But the problem persists if I increase the dofs. What happens in that case is probably the norm of a bigger vector is very big, so it fails to converge. – Schneider Mar 30 '19 at 17:07