Appendix 5: Nonlinear Computations

It is often suggested that a nonlinear computation is basically the same as a linear computation that is repeated a number of times. This is only partly true because the structural stiffness matrix - essential in linear computations - is less important in a nonlinear computation. 1

In this Appendix the essentials of a nonlinear computation are explained. This can be of help if you are going to implement an algorithm to compute a nonlinear model since very little information is available on this subject in literature.

A-5.1 Displacements and Forces

A discrete or finite element model is composed of elements. Every element has grips. 2 When the displacements of the grips of an element are known the grip forces can be computed which apparently work on the element. For computational purposes this can be represented with a black box (see Figure 92). In software engineering this black box can be referred to as a subroutine, as a procedure or as an object. Input are the grip displacements and output are the grip forces.

Figure 92: Grip displacements are used to compute grip forces.

When a model is made of a structure the grips of elements are connected by nodes. So, the displacements of connected grips will be the same and the grip forces in the nodes are added to give a resulting nodal force. As Figure 93 shows, another black box is obtained with input nodal displacements of the model and output nodal forces. This routine can be easily programmed and is the essential part of most nonlinear structural software. The nodal displacements are conveniently assembled in a vector u and they are often referred to as degrees of freedom. In this text we will call them internal displacements. The nodal forces are assembled in a vector fi which is referred to as the internal force vector.

Figure 93: Nodal displacements are used to compute nodal forces.

The load on the model can consist of external forces, described displacements, material shrinkage, prestressing etc. The support of the model is nothing but a described displacement which is in fact a load too. This Appendix is limited to forces and displacements but any load can be taken into account either analogously with applied displacements or with applied forces.

Figure 94: A structural model, somewhere in space, with displacements and forces. Note that individual elements are not drawn. At a node either the displacement has to meet an applied displacement or the force has to meet an applied force. The gaps and unbalances are enlarged in every load increment and reduced in every iteration.

Figure 94 shows a model with nodal forces and displacements. The figure is somewhat abstract but we can choose any structure and visualise it somewhere in space. In the figure the internal nodal forces and the applied nodal forces are drawn. The internal forces are computed as explained above and the applied forces are what we choose the structure has to carry. In many nodes the applied forces will be just zero. Some nodes will not have an applied force but an applied displacement instead. Also the applied displacement will be many times just zero.

Suppose there is a difference between the applied quantities (forces and displacements) and the internal quantities. The difference between an internal nodal displacement and an applied nodal displacement is called a gap. The difference between an internal force and an applied nodal force is called an unbalance. Clearly, the gaps and unbalances should be zero.

A-5.2 Mathematical Formulation

Mathematically we can write the previous as follows

In this fi is the internal force vector as a function of the vector u with internal nodal displacements. f is a vector with the applied nodal forces. Both u and f depend on the load factor l . For example:

In this example the first two displacements in the left-hand vector are described zero and the third displacement is applied a. The displacements u4, u5 and u6 have to be computed. The forces f1, f2 and f3 have to be computed too while the remaining forces are applied with 0, b and c, respectively. The problem that we have to solve is: Find the values of l and u4 to u6 that fulfil all three equations. 3

A-5.3 Nonlinear Algorithm

Suppose that the model is loaded and gaps and unbalances are present. These differences can be reduced with a computational algorithm. For this we consider the differences as small extra loads and extra displacements onto the structure. When the structure is assumed to behave linearly the extra internal displacements can be computed with the well known displacement method of linear-elastic analysis: Just assemble the global stiffness matrix, process the supports and solve the unknown displacements with the load. The resulting small displacements are added to the internal displacements we already had. Subsequently, the new internal nodal forces can be computed and if everything goes well the gaps and unbalances will be smaller than before.

This is formulated as follows:

In this K is the structural tangent stiffness matrix and du the correction of the internal nodal displacements. The inverse of the stiffness matrix is not really computed since this would not be efficient. The notation only states that the system of equations is solved. Note that the gaps must be processed in the stiffness matrix before du can be solved. More information on the stiffness matrix can be found in Section 6 of this Appendix.

The gaps are not closed completely due to computational inaccuracies and the unbalances are not completely eliminated because the structure does not really behave linearly. The remaining gaps and unbalances can be reduced in a next iteration. This is repeated until they are sufficiently small.

To compute the complete structural behaviour the load (applied forces or applied displacements) is increased in subsequent steps. In every step the above algorithm will adjust the internal nodal displacements in a number of iterations in order to keep up with the load as good as it can. Thus, the model follows the load without ever reaching it exactly, like a donkey follows a carrot that is hold in front of it (see Figure 95). The complete algorithm is shown in Figure 96.

Figure 95: The donkey will never reach the carrot that is hold in front of it.

A-5.4 Load Control

A simple method is to increase the load factor l in a number of steps with a constant value of every increment. This is sufficient for most situations and it is implemented for the nonlinear analysis of the stringer-panel model. However, when the structural strength is exhausted the model might not be able to close the differences. Instead of becoming smaller the gaps or unbalances become larger in every iteration which is usually referred to as divergence. Of course no solution is found in this situation and it is almost certain that the model has collapsed.

Figure 96: Algorithms for computation of the behaviour of a nonlinear model. The left-hand side shows the load controlled algorithm as used for the nonlinear analyses in SPanCAD. At the right-hand side the arch length controlled algorithm is shown as used for the simulations in SPanCAD

If we want to see what occurs after the peak load (force or displacement peak) the load cannot be applied unrestrictedly. Instead, we must look for what the model is able to handle. Perhaps it is damaged and instead of increasing the load it has to be reduced. So, not the load factor should be controlled but another aspect of the model has to be incremented. For example we can increase the largest internal nodal displacement with a constant value in every load step. Alternatively, the strain in a particular point can be increased in every load step.

In the SPanCAD program good experience was obtained with arch length control which was used for simulations of model behaviour. In the particular type of arch length control that was implemented, the overall displacement of the model is kept constant in every step of the computation. This can be formulated as

In this D l is the constant arch length and D u is the increase of the internal nodal displacements in a load step. With a convenient approximation this can be used to derive an expression for the load increment D l which is not constant anymore [Ramm 1981]. The computational algorithm is shown in Figure 96.

A-5.5 Convergence Criterion

A criterion is necessary to end the equilibrium iterations. Two criteria are implemented in the SPanCAD program, which both have to be fulfilled:

  1. The largest unbalance force (excluding support forces) must be 30 times smaller than the largest of the stringer grip forces and panel grip forces.
  2. The largest displacement gap of the applied displacements must be 30 times smaller than the largest internal nodal displacement.

Of course, with largest is meant the largest of the absolute values. Stricter criteria result in more iterations and consequently a larger computation time. The criteria are chosen as tolerant as possible, taking into account engineering accuracy and a nice graphical display of the computation results.

A-5.6 Smooth Relations

Reinforced concrete behaves often with abrupt changes. For example when the material cracks or when it starts yielding. Abrupt changes in the constitutive behaviour result in abrupt changes in the model behaviour. When a tangent stiffness matrix is used in the computation of the model behaviour, this will often lead to divergence of the iterative process (see Figure 97).

There are two solutions for this problem:

  1. Make the material relations smooth.
  2. Iterate with the initial tangent stiffness matrix.

In the first solution the sharp edges and jumps in the constitutive relations are rounded with polynomial relations. As a consequence the constitutive formulation will increase substantially which increases the computation time and makes checking and maintenance of the computer code very difficult. The relations cannot be made too smooth since reinforced concrete simply does not behave smoothly. So, the convergence radius is small and computation has to take place with very small steps to assure convergence.

Figure 97: The left-hand figure shows the load-displacement curve with a small disturbance in the circle. The right-hand figure shows this detail of the curve together with the load controlled computation process. The local tangent stiffness matrix at the disturbance does not represent the overall model behaviour and the iteration process goes out of control.

The second solution is very simple: Instead of iterating with the local tangent stiffness matrix the initial linear-elastic stiffness matrix is used which can be obtained when the model is not loaded. With this method the computation is very robust. The simulations and nonlinear analyses are reduced to a simple "press button action" which is essential for a design tool. A disadvantage is that the convergence goes slowly and many iterations are required to arrive at an equilibrium state. On the other hand, much time is saved because the tangent stiffness matrix does not change during the computation and only once the system of equations has to be solved instead of in every iteration.

This method is somewhat slower than the method with the consistent tangent stiffness matrix. Nevertheless, it is implemented in SPanCAD. It is so convenient that it works always without any user input that an extra minute computation time is gladly accepted.


  1. An exception is the method that uses a secant stiffness matrix in which the structural stiffness matrix is essential.

  2. In this monograph we distinguish between nodes that are part of an element and nodes that are part of the global model. An element node is referred to as grip and the global node is just called node (see Section 2.2 and Footnote 6 in Chapter 2).

  3. Essentially there are three methods to solve the nonlinear equations of a structural model:
    - Distinct element method This method consists of simple forward integrations in time. One can get a good impression of the capacity of this method from the program Working Model [Working Model 1998].
    - Quasi dynamically This method can be summarised as follows: Find iteratively the largest unbalance and make it null by adjusting the accompanying degree of freedom. The advantage is that the large system stiffness matrix found in most static and dynamic computations is not needed. It is reported to be faster than a full dynamic method [Baars 1996].
    Minimum potential energy
    This method finds the model behaviour by minimising its potential energy. Its advantage is that it can be easily extended to include optimisation of the model or the design. However, it is very time consuming because the computation time is O(n4) in which n is the number of degrees of freedom. As a consequence it is only suitable for very small models. It has been used for strut-and-tie computations [Rückert 1992].
    - Secant matrix The load is applied in one or several steps and iterations take place with an adjusted secant stiffness matrix till equilibrium is fulfilled. The secant matrix is not uniquely defined. We experienced that this method is only successful in case of weak nonlinear behaviour. However, it has been applied successfully [Jagd 1994] [Vecchio 1989].
    - Newton Raphson This method uses a tangent matrix or an ersatz matrix. It was implemented in SPanCAD and it is explained in this Appendix.