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 
  du=uload/numu;   % step size 
  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.

  • $\begingroup$ [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. $\endgroup$ – alephzero Mar 30 '19 at 16:37
  • 1
    $\begingroup$ 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. $\endgroup$ – alephzero Mar 30 '19 at 16:42
  • $\begingroup$ 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. $\endgroup$ – Schneider Mar 30 '19 at 17:05
  • $\begingroup$ 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. $\endgroup$ – Schneider Mar 30 '19 at 17:07

I think your convergence tolerance is not enough. Since you obtain first few solutions with low convergence requirements you are moving away from the real results for the next solution steps. From the code, it seems you are using unrealistic values. In example, total displacement is 5. The tolerance is 0.01, they are crude numbers. So it is not easy to determine which values must be used. Without proper use of realistic values, it is not easy for the model to converge. For a nonlinear problem, suppose you find disp1=0.01 with corresponding force1=0.01, then for disp2=0.02, assume an acceptable solution for force2=0.019, but if without properly adjusting tolerance, you may find force2=0.018, then the force3 result will definitely change based on the force2 value you obtained.

| improve this answer | |
  • $\begingroup$ this looks like comment, and not answer. Do you have any suggestions on how address convergence tolerance $\endgroup$ – Mahendra Gunawardena Oct 7 '19 at 10:43
  • $\begingroup$ May I suggest that you update the answer with content of the comment and delete the comment. Appreciate your contribution to the community. $\endgroup$ – Mahendra Gunawardena Oct 7 '19 at 12:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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