An Iterative Procedure for the Analysis of Nonlinear Elastic Systems Subjected to Uniaxial Loading and Large Displacements

The paper deals with an iterative procedure for the solution of elastic problems, for which the large displacement theory must be considered. It is shown that the algorithm, based upon a non-traditional iterative scheme, tends to converge to the correct solution and some numerical tests are considered, by assuming a linear and piecewise-linear material model.


Introduction
There is a significant number of structural systems that cannot be modelled according to the small-displacement theory. In actual fact, there are several research works concerned with the large-displacement theory. Consequently, convenient algorithms were developed in this area of structural engineering [1]. For instance, we can mention contributions that range from fluid-structure interaction [2,3] to robotics [4], from cable networks [5] to buckling [6], from elasto-plasticity [7] to dynamics [8].
The purpose of this paper is to provide a simple iterative scheme, which makes use of solutions found according to the small-displacement theory. Eventually, it allows one to converge toward the solution of structural problems for which equilibrium is imposed in the deformed configuration.
The details of the proposed approach are discussed by considering structural systems subjected to uniaxial stresses, but the method can be extended to more complex structures.
Apart from the main features of the numerical technique, it is shown why the iterative process tends to converge toward the correct solution and some numerical tests are presented in order to prove the effectiveness of the algorithm.
Namely, the sample problems focus on linear and piecewise-linear elastic stress strain laws. However, any nonlinear constitutive law can be easily implemented by maintaining the same basic framework.

The Algorithm
It is assumed that a structural system consists of linear or nonlinear elastic elements subjected to a uniaxial stress state. This Section will focus on a linear elastic behavior, which allows one to write simple equations, but does not set any limit to the range of applicability of the proposed approach. Indeed, as shown later, it is easily possible to deal with the nonlinear case.
The solution procedure is based upon an iterative method, which requires the solution of a linear elastic problem and the computation of a new stiffness matrix at each step.
During the iterative process, each structural member is subjected to a certain strain and its length is L i , where i denotes the current iteration. At the beginning, the initial length is L 0 , while both stresses and strains are zero.
After determining the linear elastic stiffness matrix of the given structure, say K 0 , which already takes the boundary conditions into account, a displacement vector u 1 can be computed by solving the linear system K 0 u 1 = F, if F represents the given vector of external loads.
When the displacements u 1 are known, it is possible to determine the elongations q s of the S structural members (s=1,…,S). More precisely, the terms q s are found by considering the exact displaced configuration, without introducing the approximations, which are typical of the small-displacement theory. In other words, for each structural member the actual length L 1 shall be evaluated by observing that the positions of the endpoints in the displaced configurations are given by the initial coordinates (x A , y A and x B , x B ) plus the relevant components of the displacement vectors (u A , u B ).
For instance, as shown in Fig. 1, the new abscissas of the endpoints A' and B', obtained in view of the displacement vector u 1 , are given by the initial abscissas (x A and x B ) plus the displacement components (u A and u B ). Eventually, the terms q s will be determined by considering, for each structural member, differences such as (L 1 -L 0 ) in Fig. 1. Next, the axial forces Q s = k s q s (no summation implied) are determined, where k s denotes the axial stiffness. By writing the equilibrium equations (with respect to the deformed configuration), the vector of nodal loads, say F 1 , which satisfies those equations, is immediately found. Now, it is possible to proceed with the next iteration. To this aim, the linear elastic stiffness matrix of the structure (K 1 ) is determined by assuming the displaced configuration as initial, undeformed state. Then, it is necessary to define the vector of the nodal loads ΔF = {F -F 1 } and solve the linear elastic problem, ΔF = K 1 Δu.
After computing the displacements u 2 = u 1 + Δu, it is easy to find the elongations q s and axial forces Q s , in order to define the load vector ΔF = {F -F 2 }, where F 2 is the vector of nodal loads, which satisfies the equilibrium equations with respect to the new deformed configuration. Then, it is necessary to proceed as before, after updating the linear elastic stiffness matrix.
The iterative procedure can be summarized as follows: 1) At the end of the i-th iteration, a new displacement vector u i = u i-1 + Δu is obtained by means of a linear elastic analysis carried out on a fictitious structure whose stiffness matrix is K i-1 ; this structure represents the configuration of the given system at the end of the previous iteration 2) The displacement vector u i allows one to find the elongations q s and axial forces Q s in each structural member on the basis of the large-displacement theory; note that the elongations are computed by considering differences such as (L C -L 0 ), where L C and L 0 represent the current and initial length of each element 3) The load vector F i is found, which is in equilibrium with Q i (the vector that collects the S axial forces Q s ) 4) A convenient coefficient is computed in order to check convergence (e.g., the ratio |F -F i | / |F| or, alternatively, the ratio |u i -u i-1 | / |u i |, where |[•]| denotes the Euclidean norm of [•]); if the selected coefficient is less than a given tolerance, the iterative process is stopped; otherwise, a new elastic stiffness matrix K i is computed, in order to determine an incremental response Δu on the basis of a new load vector ΔF = {F -F i }; then, the process continues from point (1) Clearly, the proposed algorithm cannot work if the structural system is unstable and, hence, the initial stiffness matrix is singular. A classic example is given by the two horizontal rods in Fig. 2, where the grey dot represents a hinge and the black dots denote constrained nodes. The problem, however, can be easily circumvented by introducing an extra element (represented by a dashed line in Fig. 2) during the first iteration. Thus, K 0 becomes nonsingular. Later, the fictitious element can be removed, since the updated matrix K 1 will be determined by considering the deformed configuration, so that the three hinges at the bottom are no longer aligned.
For instance, the algorithm was applied to the system in Fig 2, assuming that each horizontal element was 1m long and was characterized by an axial stiffness of 10 MNm -1 . A vertical fictitious element (with an axial stiffness of 1 MNm -1 ) was utilized and a downward vertical load of 0.1 MN was imposed. The numerical procedure was stopped when the ratio |u i -u i-1 | / |u i | turned out to be less than 10 -9 . In the end, after 33 iterations, the correct vertical displacement (0.2179628 m) was obtained.

Convergence Properties
All numerical tests, so far, have shown that the algorithm does converge toward the correct solution. It is worth noting, however, that this result was not obtained because of favorable circumstances, but because the intrinsic features of the approach discussed in the paper tend to reduce the value of the total potential energy, iteration by iteration. Hence, the algorithm progressively approaches the exact solution, which corresponds to the minimum point of the total potential energy. Unfortunately, it is not possible to give a conclusive proof of the fact that the iterative procedure necessarily keeps converging toward the solution, because the large-displacement hypothesis does not allow one to exactly determine how the total potential energy changes during each step. Nonetheless, as will be detailed in what follows, it is possible to consider a similar system for which the change of the total potential energy, say ΔΠ*, is certainly negative. Next, it can be observed that the change of total potential energy, which is caused by the algorithm, does approximate ΔΠ*.
Indeed, it is possible to consider any deformed configuration at the end of the i-th iteration, when the axial loads Q i are in equilibrium with the nodal forces F i . After modifying the mesh and assembling a new stiffness matrix (K i ), the vector of incremental displacements Δu = [K i ] -1 ΔF is found. Then, the vectors u i+1 =u i +Δu, q i+1 =q i +Δq and Q i+1 =Q i +ΔQ are updated on the basis of the large-displacement theory.
However, it is also possible to compute the elongations Δq* and the corresponding axial forces ΔQ* that would occur according to the small-displacement theory. In this way, the total potential energy is easily determined: Assuming a linear elastic material, the increment ΔW* is given by ½ ΔQ* T Δq* and is positive for any ΔF≠0. In consequence, ΔΠ* is less than zero for any vector ΔF≠0, in view of the principle of virtual works (ΔF and ΔQ* are in equilibrium, while Δu and Δq* are compatible). Therefore, it is possible to set: Clearly, this equation implies ΔΠ*<0 if ΔF≠0. Now, it should be observed that the total potential energy of the structural system (given by the large-displacement theory), at the end of the i-th iteration, reads The strain energy W i depends on the elongations and axial forces computed by taking the displaced configuration into account. Again, in the case of linear elastic materials this energy is given by the equation At the end of the next iteration, the total potential energy obviously becomes Here, ΔW = ½ ΔQ T Δq, while ΔŴ= {Q i } T Δq refers to the contribution given by Q i in view of the elongations Δq.
From eqns. (3) and (5) it is possible to obtain The last expression is written by observing that F=F i +ΔF and that ΔW approximates the increment ΔW*, because both contributions represent works done by incremental forces for incremental displacements (even though different theories are applied). Similarly, As pointed out above, it turns out that the last term on the left hand side is zero only when ΔF = 0.
At this stage, by focusing on ΔW and ΔW*, it is clear that (in the presence of uniaxial stress states) there are cases in which a contribution to ΔW is certainly less than the contribution to ΔW*. For instance, if the upper part of Fig. 3 is considered, the segment A'B' is definitely shorter than the actual length of the deformed element (thin line). Instead, the contribution to ΔW is greater than the contribution to ΔW* if the displacements shown in the lower part of Fig. 3 are assumed. In fact, the transversal displacements imply no strains according to the small-displacement theory, while they give a contribution to ΔW.
It is worth noting that, on the rare occasions in which ΔΠ was greater than zero, it usually happened when the algorithm was close to convergence (i.e., when the displacement increments were particularly small). Thus, ΔW, ΔŴ and Δq should have been practically equal to ΔW*, ΔŴ*={Q i } T Δq* and Δq*, respectively. In the end, a positive value of ΔΠ appeared to be a mere consequence of numerical errors, as confirmed by the fact that an increment ΔΠ>0 often came together with a positive increment ΔΠ*=-ΔF T Δu+ΔW* (which is absolutely impossible without rounding errors).
However, a significant result can be obtained, even if ΔΠ turns out to be positive because of mechanical reasons rather than numerical errors. In fact, it is always possible to consider a vector δF=-αΔF (with 0<α<1) that satisfies the inequality {-ΔF-δF} T Δu + ΔW -{F i } T Δu + {Q i } T Δq < 0 when the current iteration is completed. In this context, Δq and Δu denote the elongations (computed on the basis of the actual displaced mesh) and the displacements induced by the load vector {-ΔF-δF}. In this way, at the end of the iteration i+1, it is possible to write, instead of Π i+1 in eqn. (5) Because of the problem discussed here, δF, in general, is small when compared with F, so that |δF| / |F| « 1.
Here, the vector δF makes the difference (Π -Π*) negative. Therefore, even though ΔΠ may be positive for reasons that do not depend on numerical errors, the proposed method allows one to approach the correct solution, iteration by iteration, by considering a load vector F+δF ≈ F, instead of F. As pointed out above, a norm |δF| / |F| « 1 is usually obtained. This implies that, after each iteration, the algorithm converges toward the exact solution or (at the very least) toward the solution that would be obtained if the structural system were subjected to the load {F + δF} instead of F.
By implementing the same strategy a few times, it is possible to narrow the difference between F and the actual load to any degree of approximation. To do so, a sequence of smaller incremental loads can be applied (instead of imposing a single load increment ΔF), with the aim of coming up (at each iteration) with increments Δq that tend to coincide with Δq*. In consequence, the increment of the total potential energy cannot be positive, as shown below.
For the sake of simplicity, a case will be considered in which ΔF is subdivided into three increments ΔΦ 1 , ΔΦ 2 , Here, the nodal loads ΔF k and ΔΦ k are in equilibrium with the internal forces ΔQ k and {ΔQ*} k computed according to the large-displacement and small-displacement theory, respectively. Note that the mesh shall be updated before imposing each load increment. Finally, it is assumed that the relevant incremental displacements Δu k , compatible with {Δq*} k ≈Δq k , can be interpreted as virtual quantities (together with Δq k ). Obviously, the same path of reasoning is applicable to any number of increments ΔΦ k . In this way, the Euclidean norm of the vectors ΔΦ k can be less than any selected threshold.
When the load increment ΔΦ 1 is applied, the total potential energy, which (at the end of the i-th iteration) is given by eqn. (3) with F=F i +ΔΦ 1 , becomes Next, as the loads ΔΦ 2 and ΔΦ 3 are imposed, the total potential energy becomes Clearly, {F i +ΔF 1 } and {F i +ΔF 1 +ΔF 2 } are in equilibrium with the internal forces {Q i +ΔQ 1 } and {Q i +ΔQ 1 +ΔQ 2 }, computed according to the large displacement theory.
In view of the above remarks, -β k ΔF T Δu k ≈ -2ΔW k (with β k >0) is negative and, consequently, ΔF T Δu k is positive for any k. So, the terms -{ΔF 2 +ΔΦ 3 } T Δu 1 ≈-(β 2 +β 3 )ΔF T Δu 1 and -{ΔΦ 3 } T Δu 2 = -β 3 ΔF Δu 2 must be negative. In order to give an idea of what typically happens in terms of total potential energy, the structural system in Fig. 4 will be considered. It consists of thirteen elements, for which an axial stiffness EA (Young's modulus times the cross sectional area) equal to 10 MN is assumed. The nodes marked by black dots are constrained and represent the vertices of an ideal cube whose sides are 2m long. The components along x, y, z of the force on the bottom are 0.5, Civil Engineering and Architecture 3(5): 153-162, 2015 157 -05, -0.5 MN, respectively. The hinges shown in grey are not constrained, so that the system has six degrees of freedom.
The process was stopped when the ratio |u i -u i | / |u i | became less than 10 -9 . It took sixteen iterations and the trend of the total potential energy is shown in Fig. 5.  Figure 7. Trend of the total potential energy when a beam is subjected to a force P, which is greater than the critical load The total potential energy remains practically constant after two iterations, even though it continues to decrease, with one exception during the fourteenth iteration. Then, the force increments were divided by two for a few times and eventually, after five attempts, a negative increment of the total potential energy was obtained, as shown in Fig. 6.
Clearly, the most interesting feature of the plot in Fig. 6 has nothing to do with the fact that a negative increment was obtained after a few attempts. What is more important is the fact that the order of magnitude of the increments is 10 -16 :10 -15 , while the total potential energy is about -2.41 MN m. Hence, in this case, any result is obviously affected by rounding errors.
At this stage, it might be objected that there are cases in which the solution of a structural problem does not necessarily coincide with the minimum value of the total potential energy, when equilibrium equations are written by considering a displaced configuration. For instance, it is possible to consider the beam in Fig. 7 and assume that the force P is greater than the critical load. If there are no imperfections and the load is not eccentric, the straight configuration is feasible, but it is not stable. Actually, as shown in Fig. 7, the equilibrium condition corresponds to a saddle point: if a displacement δξ is imposed (starting from the deformed straight configuration, in the presence of the load P) the total potential energy increases, while it decreases if a displacement δη is imposed.
Definitely, the proposed algorithm is unable to deal with problems of this kind. Actually, it even happens that the numerical procedure gives an unrealistic solution just because it tends to find the configuration, which corresponds to the absolute minimum of the total potential energy. This is the case of the system represented in Fig. 8. When the load P becomes too high, a displaced configuration is obtained, which is practically symmetrical with respect to the left (constrained) edges. In consequence, tensile stresses occur in the horizontal bars. Similar solutions can be computed with three-dimensional systems.
Instead, a correct, realistic solution (with moderate displacements and compressive stresses) was always obtained when the model of Fig. 8 and relatively low values of P were considered, even though the absolute minimum of the total potential energy was related to a displaced configuration, which was practically symmetrical with respect to the left edge. This is due to the fact that, when P is relatively small, the point corresponding to the correct solution is at the bottom of a deep valley. Therefore, the algorithm tends to look for the optimal solution in that valley. Instead, for high values of P, it may well happen that the iterative procedure sometime arrives at a point in the deeper valley (to which the absolute minimum belongs) and eventually finds the solution, which does correspond to the absolute minimum.
In order to clarify this point, it is possible to consider the single-degree-of-freedom system in Fig. 9. It consists of a bar whose projections along the horizontal and vertical axis are 1 m long. It is assumed that the cross sectional area and the elastic modulus are equal to 0.001/√2 m 2 and 200,000 MPa, respectively.
As the total potential energy is plotted (cf. Fig. 9), it is obviously found that the absolute minimum always occurs when the displacement is greater than 2 m. However, the lower the load (e.g., F=2,000 kN) the deeper the valley whose minimum point corresponds to a relatively low value of the displacement u. Note that F≈13,251 kN is the largest load for which a solution u<1 m is feasible.

Numerical Tests
Some numerical tests were carried out by considering space networks such as the one represented in Fig. 10, which appears to be reasonably representative of the large span structures (e.g., for stadium roofs) that can be studied by imposing the equilibrium conditions in the displaced configuration. All networks are delimited by four pairs of external parabolas and their projections on the mid-plane fall inside a square whose sides are 24 m long. Each pair of external parabolas has two common nodes (on the x-y plane) and four constrained nodes. The distance between two constrained nodes, measured in the direction of the z-axis, is 1 m (cf., e.g., nodes A and B in Fig. 10). Instead, the distance between the common nodes of the external parabolas is 20 m. The end nodes of all the other pairs of parabolic curves are coincident and constrained. Finally, the top nodes of the tendons (which are also constrained) are 4 m above the x-y plane.
The algorithm was applied to structural systems that featured an increasing number of degrees of freedom. All meshes were defined by maintaining the same basic shape (as shown in Fig. 10) and changing the number of vertical elements in the planes of the external parabolas. For instance, Fig. 10 is concerned with a discrete model characterized by three vertical elements in these planes.
The positions of the nodes in the zone inside the external parabolas were determined by considering other parabolic curves on planes parallel to the x-z plane in Fig. 10. Namely, two parabolas were considered (e.g., α and β, one convex, one concave) passing through the ends of each vertical element that was positioned between the external parabolic curves parallel to the y-z plane in Fig. 10. Next, further structural members were introduced, which described curves in planes parallel to the y-z plane (e.g., γ and δ). Eventually, as many pairs of curves were defined (one convex, one concave) as the number of vertical elements that were placed between two external parabolic curves parallel to the x-z plane. Therefore, in the case of Fig.  10, there are six curves (three convex, three concave). The structural members placed along the curvilinear trajectories were reinforced by means of vertical and oblique elements. Eight tendons on the corners were also included, in order to make the structure stable (thick, dark grey lines in Fig. 10). The loaded nodes, at the intersections between the internal concave parabolas, are marked by black dots in the same figure. Their number obviously depends on the mesh. Table  1, which gives the main features of the structural systems that were used for the numerical tests, also specifies the number of loaded nodes for each discrete model. The numerical analysis was carried out by applying a total load of 40 MN and a tight tolerance was always imposed. In fact, the algorithm was stopped when the ratio |u i -u i-1 | / |u i | was less than 10 -9 (with u i denoting the displacement vector at the end of the last iteration).
In view of the special load condition (same total load for all space networks, no matter their structural stiffness and/or the location of the loading points), it was possible to check the ability of the algorithm to deal with the case of very large (even though unrealistic) displacements. This is clearly shown in Table 2, where the most significant results are reported. Of course, the largest displacement occurs in the structure subjected to a single load. It turns out to be almost 3 m (over 12% of the total length of the discrete model). For this mesh, the largest (unrealistic) strain was over 0.1.These values (together with the tight tolerance) explain the relatively high number of iterations and CPU time required to solve the problem concerned with the single load, despite the low number of degrees of freedom (only 42, much less, e.g., than the number of degrees of freedom of the discrete model characterized by 78 nodes). CPU times in Table 2 (as well as in the following tables) are purely indicative and refer to computations, which were performed on a desktop computer with an Intel Xeon E5420 processor (2.5 GHz).
The program was developed by using the Scilab open source software with standard routines/functions. For instance, the stiffness matrix was stored without exploiting its symmetry or bandwidth, and the linear systems were solved by simply inverting the relevant stiffness matrices. In other words, no special effort was made to implement a computer program, which might be suitable for commercial purposes. In fact, the only aim of this work was to develop a numerical method, which could be easily implemented by using standard routines for the structural analysis of traditional systems that are usually studied on the basis of the small-displacement theory. Consequently, any comparison with existing computer packages in terms of computational efficiency would be out of place.
Instead, it is obviously important to assess the accuracy of the solutions. This task can be easily accomplished by checking the equilibrium at the nodal points. As well known, the solutions must be unique, since instability problems are not considered. Thus, if the equilibrium conditions are satisfied, it is possible to state that the solutions are accurate.
This is exactly what happened in the case of the problems presented here, since the following results were obtained: whenever an external load, say F, was applied to a node, the internal forces turned out to be equivalent to F±δF (with δF in the range ±10 -9 F); similarly, whenever a node was not loaded, the internal forces were equivalent to a load ±δF, with δF less than 10 -9 F max , if F max is the absolute value of the the largest applied load. However, apart from rare exceptions, the absolute value of δF was definitely less than 10 -9 F max for all the numerical examples that have been considered so far.
The numerical approach outlined in the paper can be easily applied to any nonlinear elastic material. For instance, it is possible to consider the stress-strain plot of Fig. 11, which refers to a ductile iron. Part of the (curvilinear) actual plot, limited to tensile stresses, is shown in black, while the grey piecewise-linear approximation refers to the constitutive law we assumed for our structural system. The piecewise-linear graph consists of four branches with different slopes, quantified by Young's modulus (200,000 MPa) and three tangent moduli (8,108.1081, 500, 250 MPa).
Note that the changes of slope occur when ε=0.0015, 0.02, 0.04 and σ=300, 450, 460 MPa. Finally, the point whose coordinates are ε=0.06 and σ=465 MPa allows one to determine the slope of the fourth branch.
It was assumed that the cross sectional area of all structural members was 0.001m 2 and that the material response, in the presence of compressive stresses, was described by a piecewise-linear graph fully analogous to the plot in Fig. 11. On the basis of these mechanical properties, the response of the discrete models defined in Table 1 was computed. More specifically, for every mesh, three different loading conditions were imposed, in order to make the largest stresses correspond to points along the second, third and fourth branch, respectively, of the stress-strain plot.
Here, for the sake of brevity, only the response of the discrete models consisting of 38, 294 and 806 nodes will be considered. The relevant results are given in Tables 3-5. As obvious, intermediate CPU times were obtained when the algorithm was applied to discrete models characterized by a different number of degrees of freedom. Before starting a new iteration, the stiffness matrix was updated by considering (for every element) the slope of the branch concerned with the stress that had been determined at the end of the last iteration. It is worth noting that the high number of iterations needed for the third case in Table 5 is due to the applied loads (0.0523 MN), which are just above the critical value that corresponds to the onset of a collapse mechanism. Of course, a mechanism can be activated in view of the stress-strain plot, which is nearly flat for stresses greater than 460 MPa, as shown in Fig. 11. The development of a collapse mechanism is put in evidence by the largest displacement, which suddenly jumps from 0.23 m to 0.77 m, even if the increment of the applied loads is about 0.2 percent when the second and third case in Table 5 are considered.
It appears to be a sort of general rule that the number of iterations dramatically increases when a structure is subjected to a load, which approximates a critical value that enforces a significant change of the structural response.
For instance, if the single-degree-of freedom system in Fig. 9 is considered and the iterative procedure is stopped when the ratio |u i -u i-1 | / |u i | is less than 10 -12 , we need 31, 56, 136 and 2825 iterations when the loads are 10, 12, 13 and 13.251 MN, respectively. As already pointed out, 13.251 MN is close to the largest load for which a displacement u can be found, which is less than 1 m (assuming that no buckling occurs and that the material behavior is unlimitedly linear elastic).
Interestingly enough, when the load is higher than the critical value, the number of iterations tends to decrease. In fact, we obtained a solution after 229 and 50 iterations for loads equal to 13.26 and 13.5 MN, respectively.
Finally, it should be observed that a few iterations can be enough in order to obtain an exact solution (i.e. a solution that satisfies equilibrium in the deformed configuration) when displacements are relatively small. For instance, it is possible to consider the structure in Fig.  10, the stress-strain plot in Fig. 11 and the numerical model consisting of 2154 degrees of freedom. Next, as was the case of the previous tests on the same system, loads can be applied to 289 nodes. When these loads were equal to 0.01 and 0.035 MN, respectively, only 5 and 16 iterations were needed to obtain convergence. The largest displacements were 0.0134 and 0.047 m. Note that, in the second case, the highest absolute value of the stress was 29.9965 MPa, which practically coincides with the linear elastic limit.
Therefore, a simple personal computers, an easy-to-implement algorithm and a few CPU minutes can be enough to obtain accurate solutions of linear elastic problems, whenever load conditions are considered for which displacements and stresses tend to be relatively high (if compared to the largest dimension of the structural system and to the linear elastic limit).
Finally, some numerical solutions were compared with test results, which were available in the literature and were concerned with scaled model experiments [9]. So, the truss structure of Fig. 12 was considered. It consists of stainless steel tubes characterized by the following dimensions: diameter of 7 mm and thickness of 0.3 mm for the horizontal tubes, diameter of 5 mm and thickness of 0.3 mm for the diagonal tubes, with the exception of the end diagonal members (on the right and left hand side), whose thickness was 0.5 mm.
Each horizontal tube was 200 mm long and the truss height was 125 mm. Since it was also specified that the steel was ductile and the yield strength of the horizontal members was 1523 N [9], a bilinear stress-strain plot with a yield stress of 230 MPa was considered, while a tangent modulus equal to 5000 MPa was assumed. The truss was subjected to equal downward loads applied to the five free nodes at the bottom.
When the structure was tested, the horizontal tubes along the upper edge were reinforced by means of internal pipes in order to increase the buckling load [9]. Therefore, it was possible to carry out the numerical analysis of the structural system by using the algorithm presented in this paper, which does not take into account instability problems, as pointed out above. Figure 12. Numerical results compared with experimental tests [9] carried out on a truss structure Some results are shown in Fig. 12, where the total load is plotted as a function of the vertical displacement of the lower mid-point. The dots refer to the numerical solutions and are in excellent agreement with the experimental curve [9]. Note that the largest displacements turn out to be about twenty percent of the height of the truss structure.

Closing Remarks
A numerical approach has been discussed for the analysis of structural systems for which a solution given by the large-displacement theory is required. Even though it is not possible to prove that the algorithm always tends (iteration by iteration) toward the minimum value of the total potential energy, it was pointed out that the method, in general, has this feature because of mechanical reasons. In addition, it was shown that it is possible to force the total potential energy to decrease by considering sufficiently small load increments. Therefore, the algorithm does converge toward the correct solution.
As shown in the paper, when the convergence criterion is based upon a relatively tight tolerance, the solutions are definitely accurate, since the equilibrium conditions are satisfied with negligible errors. Namely, when the nodal loads should be zero, they usually turn out to be (in the worst circumstances) about ±10 -9 F max , if F max denotes the largest applied load. Similarly, when a node is subjected to a load F, the internal forces usually imply an error, which (at the very least) falls in the range ±10 -9 F. In addition, the algorithm was also tested by comparing some numerical results with test data reported in the literature [9].
The numerical technique appears to be applicable to any structural system. Here, however, only uniaxial stress states were considered. Several tests were performed that show how the method is actually able to deal with systems subjected to large displacements, provided that instability does not occur. Indeed, it was stressed that the proposed method fails to properly describe the post-buckling behavior and might require an exceedingly high number of iterations when a critical condition is approached.
Linear and piecewise-linear stress-strain laws were considered, but any nonlinear elastic system can be studied as well, on the basis of the approach presented here.
In actual fact, even elastic-plastic systems could be taken into account, since these systems are practically investigated by means of an incremental nonlinear elastic analysis, in which plastic strains (when they occur) are computed at each time step and are assumed as initial, non-reversible strains at the beginning of the subsequent step.