I'm trying to use Riks' method and I'm not sure how to set the initial values for the loading coefficient, nor the tangent vector (i.e. the derivative of the displacements and loading coefficient with respect to the path length). At the bottom of page 538, the author assumes that these are known at the outset. The best I can come up with is to set an initial value for the loading coefficient, small enough so that there is no risk of encountering any critical point, and then taking the solution as both the state vector and its normalized vector as the tangent vector. Is there a better way to do this?


  1. Riks, E. "An incremental approach to the solution of snapping and buckling problems." International Journal of Solids and Structures 15.7 (1979): 529-551.

It depends on the control which you have, if your control function is linear, for example crack mouth opening control, it will work with zero initial step. For spherical control you need kick-start analysis.

The equilibrium equation is arguments with control equation as follows, \begin{equation} \left\{ \begin{array}[l] \mathbf{r}(\mathbf{x},\lambda) = \mathbf{F}^{int}(\mathbf{x})-\lambda \mathbf{F}_{\lambda}\\ r_\lambda(\mathbf{x},\lambda) = f_\lambda(\mathbf{x},\lambda) - \Delta s^2 \end{array} \right. \end{equation} where $\mathbf{F}^{int}(\mathbf{x})$ is internal force, $\lambda \mathbf{F}_{\lambda}$ is external force, $\lambda$ is additional unknown controlled by second equation and $\mathbf{x}$ is vector of unknowns. For spherical control $f_\lambda(\mathbf{x},\lambda)$ is constrain function as follows \begin{equation} f_\lambda(\mathbf{x},\lambda) = \|\Delta\mathbf{x}\|^2 + \beta^2 \Delta\lambda^2 \| \mathbf{F}_{\lambda} \|^2 \end{equation} where the increment vector is \begin{equation} \Delta \mathbf{x}_1 = \mathbf{x}_1-\mathbf{x}_0. \end{equation} where $1$ is next step, and $0$ is initial step. The sub increment is \begin{equation} \Delta \mathbf{x}_1^{i+1} = \Delta \mathbf{x}_1^{i}+\delta \mathbf{x}^{i+1} \end{equation}

Linearising both equations we get \begin{equation} \left[ \begin{array}{cc} \frac{\partial \mathbf{F}^{int}}{\partial \mathbf{x}} & -\mathbf{F}_{\lambda} \\ 2\|\Delta\mathbf{x}\| & 2\Delta\lambda \beta^2 \| \mathbf{F}_{\lambda} \|^2 \\ \end{array} \right] \left\{ \begin{array}{c} \delta \mathbf{x}\\ \delta \lambda \end{array} \right\}= \left[ \begin{array}{c} -\mathbf{r}\\ -r_\lambda \end{array} \right] \end{equation}

Immediately you can see that initial state, you have to have non-zero $\lambda$ or $\Delta\mathbf{x}$. You can get both by solving linearised equation at beginning of analysis for some initial $\lambda$. $\lambda$ has to be small such that structure response in linear range, i.e. need one or two iterations to converge. Then you choose $\beta$, it can be zero, and then calculate $\Delta s$ from second equation. And this is kick-start for analysis.

In order to avoid problems, too small or too big steps, as an advice arc-length control works the best when you have step adaptivity, i.e. \begin{equation} \Delta s^{n+1} = \Delta s^{n} \left (\frac{I_d}{I} \right)^\gamma \end{equation} where $I_d$ is desired number of iterations, $I$ number of iterations at last step $n$, and $\gamma=0.5$.

Note, that above is when you linearise both equations. Second equations is quadratic polynomial, so you can solve it analytically. It can have some advantage near bifurcation points, where you can choose right direction. However I prefer to introduce imperfections, and use linearised system of equations.

  • $\begingroup$ Thanks a lot - this is incredibly helpful! The "kick-start" is exactly what I had in mind. $\endgroup$ – OskarM Oct 26 '17 at 0:57
  • $\begingroup$ I have a follow-up question. I want to use the method as a continuation through a critical point, so my goal is to reach unity loading coefficient, $\lambda = 1$. Is there any special trick for this, or do I just proceed until $\lambda > 1$, then set $\lambda = 1$ and solve one last time? $\endgroup$ – OskarM Oct 28 '17 at 16:56
  • $\begingroup$ Yes, that is correct, no special trick here. You can do as you write. $\endgroup$ – likask Oct 28 '17 at 17:19
  • $\begingroup$ You've been so helpful I thought I'd ask you another question. My system is pretty big - order millions of displacement degrees of freedom. Of course, I use Krylov subspace methods for solving the linear system (Newton's method). I'm having a lot of trouble finding a good preconditioner for the linear solver - simple ones like Jacobi eventually become unstable, so that while the preconditioned residual falls, the true residual oscillates and even diverges for some preconditioners. Do you have any advice on choosing preconditioners? Any references to literature would be super helpful. $\endgroup$ – OskarM Nov 15 '17 at 19:13
  • $\begingroup$ I suggest multi-grid (MG) but others like block Jacobi, incomplete LU will do. You need to do trick, you can not apply pre-conditioner to a system like above directly, before that you need write block pre-conditioner for arc-length, and apply MG to stiffness matrix only. What are you solving 1d, 2d or 3d? If you do 3d we have open source code which handles this with Jacobi Block Solver or MG or anything else from PETSc, if it is 2d or 1d you will have to do some work. Let me know we can help you. Link to code mofem.eng.gla.ac.uk/mofem/html $\endgroup$ – likask Nov 15 '17 at 21:46

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.