# ABAQUS - Does Newton's method affect the convergence rate of higher order incremental scheme such as RK4 implemented via UMAT?

Hello everyone,

I'm currently doing an analysis on the convergence rate of an incremental model that I have implemented in ABAQUS using a UMAT. Here's the BIG lines of the problem:

1) In UMAT : Solution of a differential equation using explicit Runge-Kutta 4 (RK4), which has a 4th order convergence rate, yields the stress state "S". (In my problem, stresses are applied)

2) If I understand correctly, ABAQUS solves  "K*d = s" for "d", where "K" is the stiffness matrix, "d" the displacement vector and "s" the force vector, using a Newton's method.

3)  I then compute the convergence rate of "d" by comparing it to another "exact" value.

My question is : Since Newton's method is a 2nd order method, is it possible that for nonlinear "K" matrix, it is impossible to reach a 4th order convergence rate for "d"? My results show a 2nd order convergence rate for "d".

When I apply strains instead of stresses, ("d" is known), "s" is computed by ABAQUS directly and I succesfully obtain a 4th order convergence rate for "s".

Anyone knows something about that? (If yes, do you have any papers which I could find further information on the subject?)

Regards,

Jack

### Reformulation

Let me try to reformulate the problem here.

I am trying to understand the differences between the computations ABAQUS does when we apply strains OR stresses. In my problem, I am taking ABAQUS' output (either stress if strains are applied or strains if stresses are applied) and I do a convergence analysis on that output. It appears that when strains are applied (stress is the output), the solution converges with 4th order (with respect to time-step). However, when I apply stresses (strains are measured), the solution converges with lower order (between order 3 and 2). These results are very strange and I am wondering if the fault is ABAQUS.

I hope this clears some points about my problem.

-Jack

### Numerical method accuracy depends on underlying algorithms

v\:* {behavior:url(#default#VML);}
o\:* {behavior:url(#default#VML);}
w\:* {behavior:url(#default#VML);}
.shape {behavior:url(#default#VML);}

Normal
0

false
false
false

EN-US
X-NONE
X-NONE

/* Style Definitions */
table.MsoNormalTable
{mso-style-name:"Table Normal";
mso-tstyle-rowband-size:0;
mso-tstyle-colband-size:0;
mso-style-noshow:yes;
mso-style-priority:99;
mso-style-qformat:yes;
mso-style-parent:"";
mso-para-margin-top:0in;
mso-para-margin-right:0in;
mso-para-margin-bottom:10.0pt;
mso-para-margin-left:0in;
line-height:115%;
mso-pagination:widow-orphan;
font-size:11.0pt;
font-family:"Calibri","sans-serif";
mso-ascii-font-family:Calibri;
mso-ascii-theme-font:minor-latin;
mso-fareast-font-family:"Times New Roman";
mso-fareast-theme-font:minor-fareast;
mso-hansi-font-family:Calibri;
mso-hansi-theme-font:minor-latin;}

The order of convergence in any iterative scheme is directly
related to the order of approximation made in the algorithm. For any numerical
method that consists of approximation at many levels, the order of error is
governed by the lowest order of approximation made in any intermediate schemes
and machine precision. For example, if a numerical scheme consists of algorithm
1, 2, 3 with order of error ε, ε2, ε4, respectively, then
the order of error in the whole solution will be of order of ε.

To illustrate this point, let us consider the Taylor’s
series of sine.

sin(x+ε) = sin(x) + ε cos(x) + O(ε2).

Here, the order of approximation is O(ε2).
Now, let us consider two methods: (1) sin(x) and cos(x) are calculated
exactly, and (2) sin(x) and cos(x) are calculated from series expansion.

In method (2): I will further approximate sine and cosine using series expansion:

sin(x) = x +  O(ε3), cos(x) = 1 - x2/2 + O(x4)

Note, if we calculate sin(x+ ε)
from method (2), we have three levels of approximations. And therefore the
lowest will govern.

For x = 0.5, and ε = 0.01. Assuming the machine results are exact up to six decimal places,
the following table illustrates the numerical values calculated.

Function

Exact

Method 1

Method 2

sin(x)

0.479426

0.479426

0.500000

cos(x)

0.877583

0.877583

0.875000

sin(x+ ε)

0.488177

0.488201

0.508750

As can be seen above, in method 1, the lowest order of error
is O(ε2), and therefore sin(x+ ε) shows accuracy until 3rd decimal place. In method 2,
the error was accumulated from all the approximations; i.e.

Numerically sine function made largest error, and therefore final
solution is governed by sin(x).

Just for example, here I made very crude illustration in taking
the values and order of approximation. But it does show, how error propagates
in numerical analysis.

Similar things are happening in force controlled and displacement
controlled analyses. You have a UMAT, to which abaqus passes strains (STRAN and
DSTRAN), and the UMAT calculates and returns the updated stress (STRESS) tensor
and Jacobian matrix (DDSDDE). Abaqus uses this Jacobian matrix to generate
stiffness matrix, and solves the equations for new displacements/forces.

In the case of force controlled, abaqus will assume and value of
strain increment and will pass it to UMAT. Your Umat will calculate O(4) stress
(which is quite accurate), and will return it to abaqus. Now abaqus will use it
to solve the equations, and will make necessary adjustments based on Newton’s
method (O(2)). So, even though you had very good accuracy in stress
calculations, the displacements are still governed by newtons method, and
therefore you see quadratic convergence only. This is similar to Method 2.

On the other hand, when you make it displacement control, Abaqus
does not need to approximate displacements (hence strain and strain
increments), therefore your UMAT governs the convergence. This is similar to Method 1.

On the side note: this is just coincidence, there are numerical
schemes running inside the solver. So in the next case, even though you run analysis
displacement control, you might not find the order 4 convergence. I hope this answers your question.

### Newton's method affects precision only

Dear Akumar,

Thank you for sharing your knowledge with me. It is extremely appreciated. You are perfectly correct about the convergence levels. Your example that illustrates that concept is very straightforward.

On the application to ABAQUS, I think you also perfectly understand my problem. However,  I believe that I must disagree on certain aspects.

1)  As you mentionned, in the case of force controlled, my UMAT calculates a O(4) stress a using RK4 method. Then, ABAQUS will use it to solve, via a Newton's method, for the displacement. The issue lies in your conclusion here : "So, even though you had very good accuracy in stress calculations, the displacements are still governed by newtons method, and therefore you see quadratic convergence only".

Newton's method, which is an iterative method, converges towards the solution until a "stop" criteria is attained (e.g. 1E-12). The notion of convergence rate for the Newton's method is based on the error of the current step compared to error to the previous step. In other words, Newton's method converges towards that "stop" criteria with order 2. Hence, only the precision is affected (i.e. the stop criteria). In the case of the Runge-Kutta 4 method, the convergence rate with respect to the size of the time-step is considered. More specifically, RK4 converges towards the solution with respect to the time-step with a 4th order convergence.

The meaning of that is that Newton's method should only affect the precision of the results (through the stop criteria) and NOT the convergence rate of the solution with respect to the size of the time-step.

2)  In your side note, you mentionned:  there are numerical schemes running inside the solver.

Do you have any idea of what these schemes are? In my analysis, when displacements are applied, all my tests resulted in a 4th order convergence. However, when stresses are applied, I get convergence rates that range from [2, 4] (depending on how I define my material), which is very strange.

My take on this problem is that I am pretty sure the issue is ABAQUS (probably not Newton either). As you mentionned, there are some underlying computations that I am not aware of that complicates the case where stresses are applied.

Regards,

- Jack

Firstly,the quadratc convergence of Newton's method is valid if following conditions are satisfied:

1. xn isclose to the root of f(x) =0
2. f'(x) != 0
3. f''(x) is finite.

This implies quadratic convergence is the best possible convergence rate (if newton's method converges to a solution).

Now in abaqus, force convergence (tolerance = 5E-3) is looked at first
before displacement (0.01). If you look at the criteria for force
residual and time-average force, it's quite elaborate, whereas
displacement equation is quite simple. Now, in your force controlled case, the approximations are made at both levels (force as well as displacements). With bigger time-step, the guess of displacement increment (hence, strain increments) are farther than the actual solution. This means,you'll iterate relatively more (and will have more error in the obtained solution).

In my experience, Abaqus utilizes Newtons method not to obtain the solution, but to obtain initial guess for displacement increment.That's why in case of displacement controlled run, there is minimal guess in the increment. Whereas in force controlled the obtained force from stress tensor is compared with the actual force, and the error is residual. Newton's method is now used to get the new displacement, and this process continues until it convergence is achieved.

Abaqus has very elaborate iteration criteria: e.g. 3 iterations after which divergence check is made (if the current residual increases relative to max of previous two, time increment is reduced by factor 2); 8 iterations after which logarithmic chekc (tolerance 0.02) is used; if more than 10 iterations are taken time increment is automatically reduced by a factor 0.75, and so on. Therefore, your convergence rate is highly dependent on all these (along with the shape of force displacement).

What I can suggest is to study the effect of tolerances (specially force and displacement).

Kumar

--

The world started with 0, is progressing with 0, but doesn't want 0.

### Also note that no

Also note that no commercial code will disclose its algorithmic implementation. At most they just mention the algorithm, but not how it is implemented (unless it is open-source) :)

Kumar

--

The world started with 0, is progressing with 0, but doesn't want 0. 