## You are here

# 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?)

Thanks in advance!

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.

Thanks in advance!

-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-padding-alt:0in 5.4pt 0in 5.4pt;

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

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

FunctionExactMethod 1Method 2sin(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ε) shows accuracy until 3rd decimal place. In method 2,

is

O(ε2), and therefore sin(x+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.

Back to your question: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 numericalschemes 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

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

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.