In this post, we will explore the importance of incremental and iterative processes in nonlinear finite element analysis (FEA), emphasizing the Newton-Raphson method. While nonlinearities in material behavior and contact often require sophisticated solution techniques, the true challenge lies in effectively implementing iterative methods that converge toward accurate solutions. The Newton-Raphson method plays a crucial role in this context, offering a systematic approach to solving nonlinear equations that arise during the analysis. We will break down the core principles of this method, its application in FEA, and how it ensures stable and efficient convergence.
Types of Nonlinearities in FEA
In finite element analysis (FEA), for structural problems, an analysis is considered nonlinear when the relationship between forces and displacements is not linear. Since stiffness is directly related to the tangent of the force-displacement curve, it changes whenever this curve deviates from a straight line.
Nonlinearities in FEA generally fall into four main categories:
- Material Nonlinearity: Occurs when the relationship between stress and strain becomes nonlinear, such as in the plastic regions of metals.
- Large Displacement Nonlinearity: Since the stiffness matrix depends on the geometry of the structure, significant changes in shape during deformation will alter the stiffness (e.g., a fishing rod bending under the weight of a large fish).
- Large Deformation Nonlinearity: The stiffness of a bar element is calculated using its cross-sectional area (k = AE/L). If the bar stretches significantly, the cross-sectional area decreases (poisson’s effect), which leads to changes in stiffness.
- Contact Nonlinearity: Contact can be considered as a new boundary condition that arises during the analysis. When contact occurs, additional constraints are applied at the contact nodes, altering the stiffness matrix.
The Necessity of Increments and Iterative Process
When we finish preparing the structural model, the stiffness of the structure is obtained from the stiffness of each of its elements. This is a straightforward assembly process. The stiffness obtained through this assembly process is the “initial stiffness.” In a nonlinear FEA, as the loading progresses, the structure’s stiffness must be updated.
The process of seeking equilibrium in the structure for a nonlinear finite element analysis as it is subjected to external loading is carried out incrementally. In nonlinear analysis, this necessity becomes one of the fundamental pillars of the study since we know that stiffness varies in this type of analysis. As the displacements in the next increment depend on the knowledge of the stiffness during that increment, and this stiffness is not known, the only alternative in the most general case is to proceed through iterative processes inside each increment.
One of the most important methods for understanding iterative strategies is the Newton-Raphson method. Although there are other methods that are much more robust and suitable for various types of nonlinearities, studying this method is essential to grasp the logic of the iterative process, which is the primary goal of this text. Additionally, it is significant because this method represents one of the first procedures developed in this context. Many schemes associated with the definition of iterative processes used the Newton-Raphson method as their starting point.
Newton-Raphson Method
In a nonlinear FEA, the central idea is the pursuit of equilibrium between the external forces and the corresponding internal forces applied to the nodes. This is achieved through a trial proposed for the displacements in the increment being analyzed. This equilibrium can also be evaluated by the equivalence between the external work and the internal work associated with the proposed deformation and the calculated internal forces.
Nonlinear analysis requires repeated solutions of the equilibrium equations to ensure that the calculated internal forces are consistent with the proposed displacements for all the nodes in the model. This is done in an incremental process where a load increment is applied.
In summary, a load increment is considered first. Then, a displacement increment is proposed for all the nodes in the model, based on the initial stiffness of the increment. This displacement increment generates nodal internal forces. If the internal forces do not balance the external nodal forces of the proposed increment, the equilibrium condition is not satisfied. In other words, this deformed configuration does not satisfy the equilibrium equation, so another configuration must be proposed until the external forces are balanced by the internal forces. This is the iterative process.
To understand the process of seeking the equilibrium condition using the Newton-Raphson method through multiple iterations, we will examine what happens within a single increment. Once equilibrium is achieved in this increment, the same process will be repeated for the next increment. The starting stiffness for the next increment is the stiffness of the structure at the end of the previous increment, which was determined based on the equilibrium obtained through successive iterations.
Despite the simplicity of the Newton-Raphson method, many solvers or calculation routines in nonlinear analysis programs are variations of this method.
There are two approaches to defining the iterative process using the Newton-Raphson method: the standard Newton-Raphson method and the modified Newton-Raphson method. We will first address the standard method and then provide some comments on the modified method.
In the image below, we see an example of a steel structure subjected to a load P applied at a specific node. Imagine that this structure exhibits nonlinear behavior due to large displacements, for instance. As a result, when analyzing the chart of load P versus displacement in the direction of the applied load, we observe a nonlinear relationship.
The figure below shows a diagram illustrating the relationship between force and displacement in a nonlinear analysis. For equal load increments, we have different displacement increments, which means that stiffness varies from segment to segment. The focus now is on just one increment, as mentioned earlier. In other words, the highlighted increment in the figure below is considered in isolation later. In each increment, a new equilibrium search problem occurs, and within each increment, various iterations take place. These iterations will be the focus of our understanding from now on.
Study of One Load Increment
In this section, we will examine what happens during a single load increment (the one circled in the figure above). Note that for the purposes of this study, when we plot load versus displacement, we will set the load and displacement at the beginning of the increment as the origin of the coordinate system, solely for the sake of visualization.
First estimate for the nodal displacements: First iteration
In this initial stage, if we were dealing with the first load increment, the structure’s stiffness matrix would be assembled from the stiffness matrices of its elements. At this point, no loads would have been applied to the structure. The nodal displacements for the load increment F (representing P_{n} ) would then be estimated using the stiffness matrix obtained from the assembly process—referred to as the “initial stiffness,” which is determined before any iterations begin.
However, since this is not the first load increment, the “initial stiffness” is instead evaluated based on the stiffness at the end of the previous increment.
Thus, the first estimate of the displacements will be given by:
{\{\triangle\}}_0={\lbrack K\rbrack}_0^{-1}\cdot\{F\}
The figure above illustrates “where we are” in this first iteration. Note that with the line defining the initial stiffness matrix K_0 , we calculate the displacement \triangle_0 , but the actual displacement would be \triangle . The error that must be corrected through successive iterations becomes evident, which essentially represents the correction or update of the stiffness.
It is interesting to note that the graph above represents the evolution of the displacements associated with one of the degrees of freedom of the model. The displacement obtained results from the matrix operation that calculates all the displacement components through the stiffness matrix and the load increment. For each degree of freedom of the model, we associate a graph of this type.
Calculation of “new stiffness”: Second iteration
After calculating the nodal displacements based on the initial stiffness, or more generally, the stiffness known at the beginning of the increment (which, in the general case, would be the stiffness known at the end of the previous increment), we will attempt to improve the calculated displacements. This will be done through the correction of the stiffness matrix.
{\lbrack K\rbrack}_1={\lbrack K\rbrack}_0+{\lbrack K\rbrack}_{correction}We have a different “correction” stiffness matrix for each iteration. Therefore, let’s refer to this as the {\lbrack K\rbrack}_{correction,1} and abbreviate it as {\lbrack K\rbrack}_{c,1} .
{\lbrack K\rbrack}_1={\lbrack K\rbrack}_0+{\lbrack K\rbrack}_{c,1}\{F\}={\lbrack K\rbrack}_1\cdot{\{\triangle\}}_1{\{\triangle\}}_1={\lbrack K\rbrack_1}^{-1}\cdot\{F\}The “correction” stiffness matrix is not arbitrary; it is carefully calculated using various techniques, depending on the type of nonlinearity being addressed.
For instance, in cases of large displacement nonlinearity, the displacements are so significant that ignoring the geometry of the deformed structure would result in unacceptable errors. To account for this, the stiffness matrix must be updated based on the deformed geometry.
Here’s how this process works: after calculating the initial displacements, \triangle_0 , using the initial stiffness matrix K_0 , these displacements are used to compute a new stiffness matrix. This updated matrix can be seen as the sum of the initial stiffness matrix K_0 and a correction stiffness matrix, often called the geometric stiffness matrix for large displacements nonlinearity. Using this approach, the updated stiffness matrix K_1 is obtained, which can then be used to calculate the next displacement, \triangle_1 .
The problem is that the displacement \triangle_0 calculated earlier is not correct because the true displacement for that load increment is \triangle . As a result, the updated stiffness matrix K_1 is also not correct since it was based on the inaccurate displacements \triangle_0 .
Calculation of the internal nodal forces using the updated stiffness after the first iteration
The figure above shows the value of the internal force after the first iteration, calculated from the stiffness matrix {\lbrack K\rbrack}_1 . In other words, after the first iteration, the displacements were calculated based on the initial stiffness and were given by {\lbrack\triangle\rbrack}_0 . From these displacements, an estimate for the correction stiffness matrix {\lbrack K\rbrack}_{c,1} was calculated, resulting in the updated stiffness matrix {\lbrack K\rbrack}_1 . Since this deformed shape {\lbrack\triangle\rbrack}_0 corresponds to a set of internal forces in the structure, we will calculate these internal forces based on the stiffness matrix {\lbrack K\rbrack}_1 , which has already been corrected.
Thus, the internal forces in the structure after the first iteration will be given by:
{\{F\}}_0'={\lbrack K\rbrack}_1\cdot{\{\triangle\}}_0After the first iteration, the external force equals F, and the internal force equals F_0' . The two are different, which indicates that the structure, in this initial attempt, is not in equilibrium. An error must be evaluated.
Compute the error of the solution
As shown in the previous item, there is a difference between the applied external forces and the internal forces, meaning there is an error in the solution. In other words, the proposed deformed shape, which was calculated using the initial stiffness, does not generate internal forces that are in equilibrium with the external load applied in the considered load increment.
This error in the solution is given by the difference between the applied external loads and the internal forces generated, which should be in equilibrium with the external loads. In other words:
{\{E\}}_1=\{F\}-{\{F\}}_0'Solve the structure for the correction of displacements
This error, calculated from the difference between external forces and internal forces, represents the internal force difference that would be “missing” in the equilibrium condition, since the external force is already known. In other words, there is a portion of the displacements that would generate this internal force difference, which was not calculated. Therefore, with this force difference given by {{E}}_1 after the first iteration, we could calculate the displacement difference that would balance this “missing” internal force in the equilibrium and the displacements that were calculated.
Strictly speaking, this is an approximation, since we are working with the stiffness matrix {\lbrack K\rbrack}_1 , but in the previous calculations, the displacements used were those given by the matrix {\lbrack K\rbrack}_0 . Thus, the relationship between the difference in internal forces and the difference in displacements, calculated based on the stiffness matrix {\lbrack K\rbrack}_1 , is given by the equations below:
{\{E\}}_1={\lbrack K\rbrack}_1\cdot\{\delta\triangle_1\}\{\delta\triangle_1\}={\lbrack K\rbrack}_1^{-1}\cdot{\{E\}}_1Determine the new total displacement
The new total displacement after the first iteration will be given by:
{\{\triangle\}}_1={\{\triangle\}}_0+\{\delta\triangle_1\}Repeat previous steps until the force imbalance is small
Typically, a tolerance is set so that the maximum error in the computation of the internal force, with reference to the values of the external force applied in the considered load increment, is very small, a value defined as \varepsilon .
The force imbalance test is carried out by calculating the norm applied to the force imbalance vector {E} and the applied external forces. This error, as mentioned earlier, must be minimal.
If we consider a finite element model, the load applied by the elements at the nodes should be in equilibrium with the equivalent nodal load. That is, in each of these equilibrium equations, we could account for an error, which would be given by a column matrix, and correspondingly, another matrix for the applied loads. Thus:
\{E\}=\begin{bmatrix}E_1\\E_2\\.\\.\\.\\E_n\end{bmatrix}\{F\}=\begin{bmatrix}F_1\\F_2\\.\\.\\.\\F_n\end{bmatrix}In the two previous equations, n represents the number of degrees of freedom and, consequently, the number of equilibrium equations. For each equation, there is an error. If we multiply these errors for each equation, summing them over the entire structure, and do the same procedure for the external forces for each equilibrium equation, i.e., multiplication and summation, we will obtain the norm calculation, which allows us to assess the error, and it should be smaller than the tolerance \varepsilon . Thus, for an iteration i, we will have:
\frac{{\{E\}}_i^T\cdot{\{E\}}_i}{{\{F\}}_i^T\cdot{\{F\}}_i}<\varepsilonThird iteration
Fourth iteration
Fifth iteration
Notice that in the fifth iteration, the difference between the true displacement and the calculated displacement is small, and the error is also small. If the norm of this error is below the tolerance, we can proceed by using the most recent displacements {\lbrack\triangle\rbrack}_4 and stiffness {\lbrack K\rbrack}_4 as the starting point for the next load increment.
Numerical Application of the Newton-Raphson Method
In this numerical example, we will assume a one-dimensional application, considering a nonlinear spring characterized by the fact that its elastic constant, or stiffness, varies with displacement. The equation that expresses the relationship between the forces and displacements of the spring is given by the equation below. The units of force and displacement can be considered in any commonly used system.
K=6-3u
The stiffness of the spring, as we can see, is not constant; it varies with displacements u. The problem is, therefore, nonlinear. In this particular application, it is assumed that the relationship between forces and displacements is only valid as long as K is a positive number. For displacement values u that result in a negative K, the problem loses its physical meaning.
For instance, when the displacement u leads to K=0, it means, in our application, that the spring has stopped functioning and has collapsed, i.e., it has lost its stiffness. This represents the working limit of the system, and our mathematical model only applies within these limits.
It is interesting to note that for u=0, that is, when the spring begins to work in the initial stages of loading, its stiffness corresponds to the “initial stiffness.” Then, we have:
u=0\rightarrow K=6-3u=6-3\cdot0=6\rightarrow k=6
Therefore, K=6 corresponds to the “initial stiffness.”
In this problem, since the analytical formulation is provided, an analytical solution is available. The objective is to compare the analytical solution of the problem with the approximate solution obtained using the Newton-Raphson method, considering the error in the equilibrium calculation between external and internal forces acceptable only up to four decimal places (0.0000).
Using the previous equation, we can prepare a set of values that, based on the analytical formulation, establish the corresponding displacement for each applied force or, for each displacement u considered, the corresponding force value. This applies within the limits where K is always positive, as previously established. Thus, we obtain:
Numerical Solution
Let us now move on to the numerical solution, considering a load increment of \triangle P=1.53 , which corresponds, in the analytical solution, to a displacement value of \triangle u=0.3 , as shown in the table above.
Within this load increment of \triangle P=1.53 , we will perform several iterations to obtain the displacement value with an acceptable approximation. The procedure we use for this load increment, once the desired value is achieved, should be repeated for subsequent increments. At the start of the next increment, the “initial stiffness” is the value determined at the end of the previous increment, which allowed the displacements to be calculated with the required precision for that increment. We will apply all the steps of the Newton-Raphson method, as previously discussed.
The central idea of this example is to obtain the displacement result \triangle u=0.3 for the load increment \triangle P=1.53 , using iterations. In the more general case, this procedure is carried out by updating the stiffness matrix through trial attempts. In this example, we will use the given expression for the spring stiffness calculation as a reference for the varying stiffness, merely to illustrate that through trial and error, we can arrive at the displacements obtained analytically. The purpose of this application is to understand the “logic” of how the iterative process works.
Calculation of the First Displacement Estimate Using the Initial Stiffness
In this numerical example, K=6 corresponds to the “initial stiffness.” Thus, applying the concept learned before in this particular case, we will determine the first estimate for the displacement u of the node in the figure above. Then, we can write:
u_0=\frac{\triangle P}K=\frac{1.53}6=0.255u_0=0.255
It is interesting to note that the displacement calculated in this first attempt for \triangle P=1.53 is different from the exact displacement given by \triangle u=0.3 .
Update of the spring stiffness
Since the stiffness varies with the displacements, and we have the first displacement estimate, we will update the stiffness for this calculated displacement. Thus, we have:
K_1=6-3u_0=6-3\cdot0.255=5.235
K_1=5.235
Calculation of the Internal Force in the Spring, Associated with the Displacement u_0=0.255 and the Stiffness K_1
The internal force will be given by:
F_0'=K_1\cdot u_0=5.235\cdot0.255\rightarrow F_0'=1.33493
Calculation of the solution error
A força interna é dada por F_0'=1.33493 , e a força externa aplicada que gerou essa força interna é dada por \triangle P=1.53 .
Therefore, it is evident that there is an imbalance between the internal force and the external force, that is, an error that quantifies this imbalance. Thus, we have:
E_1=\triangle P-F_0'
E_1=1.53-1.33493
E_1=0.19507
Calculation of the displacement correction based on the calculated error
If there is an error in the internal force, it is associated with a part of the displacements that are “missing” in the equilibrium condition. As we can see, the internal force is smaller than the applied external force, meaning the spring should be more deformed than the first displacement estimate predicted. This “lack of internal force” is associated with a “lack of displacement.” So, using the updated stiffness after the first attempt, and with the internal force error, we will calculate the first attempt for the displacement error and then add this displacement difference to the one calculated in the first iteration.
E_1=K_1\cdot\delta\triangle_1
\delta\triangle_1=\frac{E_1}{K_1}\delta\triangle_1=\frac{0.19507}{5.235}\delta\triangle_1=0.03726
Calculation of the total displacement after the first iteration
The total displacement after the first iteration is calculated as follows:
u_1=u_0+\delta\triangle_1
u_1=0.255+0.03726
u_1=0.29226
It is interesting to note again that the displacement calculated in this second attempt for \triangle P=1.53 given as 0.29226, is different from the exact displacement u = 0.3 shown in the previous table. However, it is closer to the desired value than the one provided by before (0.255). Therefore, the first attempt, or first iteration, has been completed.
We will repeat the procedure in the same way as we did in the previous steps, but now in a more direct manner, presenting the consecutive steps so that the values are updated in each iteration. The detailed justification was provided in the theory of the previous section, and we have just developed the details that support the numerical values. Let us then proceed with the numerical repetition of the previous steps.
Repeat previous steps until the force imbalance or the error is small
End of the first iteration:
K_1=5.235
u_1=0.29226
Second iteration
K_2=6-3u_1=6-3\cdot0.29226=5.12322
F_1'=K_2\cdot u_1=5.12322\cdot0.29226\rightarrow F_1'=1.49731
E_2=\triangle P-F_1'=1.53-1.49731=0.03269
\delta\triangle_2=\frac{E_2}{K_2}=\frac{0.03269}{5.12322}=0.00638u_2=u_1+\delta\triangle_2=0.29226+0.00638=0.29864
Third iteration
K_3=6-3u_2=6-3\cdot0.29864=5.10408
F_2'=K_3\cdot u_2=5.10408\cdot0.29864\rightarrow F_2'=1.52428
E_3=\triangle P-F_2'=1.53-1.52428=0.00572
\delta\triangle_3=\frac{E_3}{K_3}=\frac{0.00572}{5.10408}=0.00112u_3=u_2+\delta\triangle_3=0.29864+0.00112=0.29976
Fourth iteration
K_4=6-3u_3=6-3\cdot0.29976=5.10072
F_3'=K_4\cdot u_3=5.10072\cdot0.29976\rightarrow F_3'=1.52899
E_4=\triangle P-F_3'=1.53-1.52899=0.00101
\delta\triangle_4=\frac{E_4}{K_4}=\frac{0.00101}{5.10072}=0.000198u_4=u_3+\delta\triangle_4=0.29976+0.000198=0.29999
Fifth iteration
K_5=6-3u_4=6-3\cdot0.29999=5.10003
F_4'=K_5\cdot u_4=5.10003\cdot0.29999\rightarrow F_4'=1.52996
E_5=\triangle P-F_4'=1.53-1.52996= 0.00004
\delta\triangle_5=\frac{E_5}{K_5}=\frac{0.00004}{5.10003}=0.00008u_5=u_4+\delta\triangle_5=0.29999+0.00008=0.300007
Conclusion
The difference between the internal force and the external force calculated in the fifth iteration meets the convergence criterion proposed at the beginning of the analysis.
In this example, as mentioned earlier, the expression defining how the stiffness of the spring varies with displacements is provided and was used as a basis for assuming the stiffness variation. The main idea is to understand the “logic” of the iterative process, and this pre-established information on stiffness variation facilitated the numerical approach.
In real cases, of course, the expression defining how the stiffness varies is not provided. This variation is proposed by making an initial assumption to calculate how the displacements evolve as the structure is loaded. This is where tools such as the Taylor series come into play to define the next step based on the previous one.
We have already mentioned the need for nonlinear analysis to be incremental. The primary goal of this application is not to focus on defining a new increment based on the previous one but rather to understand the “logic” of the iterative process.
References
- Elementos Finitos, A Base da Tecnologia, Análise Não Linear – Avelino Alves Filho
Learn More!
Check out my Finite Element Method Course by clicking here.
Enhance your understanding of FEA with the following blog posts:
