I am writing a Matlab based script for solving nonlinear FEA problems using Modified Riks method by Crisfield. As a starting point, I am solving $y = x^3 - 2 x^2 - x + 2$. Around the function minima at $x \approx 1.5486$, my results evolve backward.

If we know solution at a previous load step, we can determine the solution at the first iteration of the current load step using the following [from An Introduction to Nonlinear Finite Element Analysis by J.N. Reddy]:

  1. Calculate tangent stiffness matrix.
  2. Solve for the tangential solution using the following $\delta \hat{U}_{n} = K_T^{-1} \overline{F}$ after applying boundary conditions, where $\overline{F}$ is the reference load.
  3. Compute the initial incremental load parameter, $\delta \lambda_{n}^{1}$ using $\delta \lambda_{n}^{1} = \pm \left( \Delta_s \right)_n \left[ \left\{\delta \hat{U}_{n}\right\}^T \left\{\delta \hat{U}_{n}\right\} \right]^{-1/2}$

"The plus sign is for continuing the load increment in the same direction as the previous load step and the negative sign is to reverse the load step. The sign follows that of the previous increment unless the value of determinant of the tangent matrix has changed in sign."

  1. Computer the incremental solution using $\delta \overline{U}_n^1 = -K_T^{-1} R_n^{1}$
  2. Update the total solution vector and load parameter.

$$\delta U_n^1 = \delta \overline{U}_n^1 + \delta \lambda_n^{1} \delta \hat{U}_{n}$$

Here, subscripts represent load step number and superscripts represent iteration number. Since I am using Modified Newton-Raphson method, $\delta \hat{U}_{n}$ needs to be evaluated once per load step.

Now please consider the following image.

Let us assume that we know the solution at load step No. 27, and we are trying to evaluate the solution for the next load step. In the last load step, the total change in the load parameter, $\Delta \lambda_{27}$ is positive ($y_{26} = -0.6250 < y_{27} = -0.6240$). Also, the slope at point 27 is positive, whereas the slope at point 26 is negative. Since the determinant of the tangent stiffness matrix has changed its sign, $\delta \lambda_{28}^{1}$ will be negative. However, this choice of the initial incremental load factor results in the solution evolving backwards.

In such cases, how can I make sure that the solution does not evolve backwards? One way is to choose such a value of the initial load parameter that displacement increases, but it won't work for snap-back problems.


1 Answer 1


You are correct that you cannot control the direction of the continuation in general just from looking at the loads and residuals. There has to be some sort of notion of direction inherent to the geometry of the curve you are continuing in the $(x,y)$ space. The correct way to do this is to require that the dot product of the previous step (considered as a vector in the $(U,\lambda)$ space) and the current proposed step maintains the same sign. I am not totally familiar with your notation, but I believe this condition should look something like this: $$ \begin{aligned} \delta \lambda_n(\delta\hat{U}_n; 1)^T(\overline{U}_{n-1};\delta\lambda_{n-1}) &> 0 \\ \implies \delta\lambda_n(\delta\hat{U}_n^T\overline{U}_{n-1} + \delta\lambda_{n-1}) &>0. \end{aligned} $$ The sign of $\delta \lambda$ should be chosen to satisfy this inequality at each step.

Here is an example of this procedure performed on a somewhat well-known continuation problem of a version of the Chandrasekar $H$-function, which satisfies the following integral equation: $$ H(\mu;c) = \frac{1}{1 - \frac{c}{2} H(\mu;c) \displaystyle \int_0^1 \frac{\mu}{\mu+\nu} H(\nu;c)~\mathrm{d}\nu}. $$ The solution to this integral equation can be discretized and continued numerically in the parameter $c$ to obtain multiply solutions at the same value of $c$, shown in the figures below. Note that I have the continuation parameter on the horizontal axis and the solution on the vertical axis, so it is like a 90$^\circ$ rotation of your image.

Notice that at $c=1$, we have a change in direction of our tangents along the continuation curve. That is, we stepped in a positive direction in $c$, but now our steps $\delta\hat{U}_n$ and $\overline{U}_{n-1}$ are in opposite directions, so $\delta\hat{U}_n^T\overline{U}_{n-1}<0$. By taking small enough $\Delta_s$ and ensuring the above inequality is satisfied, this implies that $\delta \lambda_n$ switches sign from positive to negative to ensure that we are traveling along the same direction with respect to the orientation of the continuation curve.

enter image description here enter image description here

  • $\begingroup$ Thank you. Please correct me if I am wrong, but what are you suggesting is that I should take a dot product of the incremental displacement vector, $\Delta U_n^{\left( r - 1 \right)}$ (total change in the displacement before the present iteration) and $\Delta U_n^{\left( r \right)}$ (total change in the displacement vector after the present iteration)and it should be positive. Right? This will ensure that the direction of incremental displacement remains consistent across iterations. This is exactly what I have been doing but for iterations after the first one. $\endgroup$ Commented May 4 at 22:24
  • $\begingroup$ How can we modify this approach for the first iteration? As per Crisfield, the initial incremental loading parameter (i.e., for the first iteration), $\delta \lambda_n^{1} = \pm \delta s / \sqrt{\delta \hat{U}_n^T \delta \hat{U}_n}$, where, the sign follows that of the previous increment unless the determinant of the tangent stiffness matrix has changed sign, in which case, a sign reversal is applied. sciencedirect.com/science/article/abs/pii/B9780080272993500091 $\endgroup$ Commented May 4 at 22:32
  • $\begingroup$ No, if you ensure that $\Delta U_n^{r-1}$ and $\Delta U_n^{r}$ are in the same direction, regardless of the direction of the $\lambda$ step, then you can run into problems. The idea is to consider the vectors in the joint $(U,\lambda)$ space and ensure those point in the same direction. Graphically, your method says "$x$ must always step in the same direction" whereas the method I propose says "The continuation curve in $(x,y)$ space must always be progressed in the same orientation," i.e., the curve is incremented by aligned tangent vectors $\endgroup$
    – whpowell96
    Commented May 4 at 23:15
  • $\begingroup$ Here are some sources that have more discussion and better figures en.wikipedia.org/wiki/… classes.engineering.wustl.edu/2009/spring/mase5513/abaqus/docs/… $\endgroup$
    – whpowell96
    Commented May 4 at 23:16
  • $\begingroup$ As for the first iteration, you just have to pick the sign that continues towards the load values you care about. If you care about values in $[0,1]$, start at 0 and increase $\lambda$ or start at 1 and decrease $\endgroup$
    – whpowell96
    Commented May 4 at 23:25

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.