Numerical-Computational Model for Nonlinear Analysis of Frames with Semirigid Connection

A numerical-computational model for static analysis of plane frames with semirigid connections and geometric nonlinear behavior is presented. 'e set of nonlinear equations governing the structural system is solved by the Potra–Pták method in an incremental procedure, with order of cubic convergence, combined with the linear arc-length path-following technique. 'e algorithm pseudo-code is presented, and the finite element corotational method is used for the discretization of the structures.'e equilibrium paths with load and displacement limit points are obtained. 'e semirigidity is simulated by a linear connection element of null length, which considers the axial, tangential, and rotational stiffness. Nonlinear analyses of 2D frame structures are carried out with the free Scilab program. 'e results show that the Potra–Pták procedure can decrease the number of iterations and the computing time in comparison with the standard andmodified Newton–Raphson iterative schemes. Also, the simulations show that the connection flexibility has a strong influence on the nonlinear behavior and stability of the structural systems.


Introduction
Structural analysis/design methodologies undergo a paradigm shift, in which linear analyses (with adaptations for consideration of nonlinear effects) are being progressively replaced by analyses capable of encompassing various nonlinear effects, such as second-order, material inelasticity, rigidity of the connections, soil-structure interaction, dynamic effects, among others [1]. Buckling and postbuckling analyses of frames are important in structural design, particularly for slender elastic structures [2]. e analysis of beams and columns as independent members can lead to erroneous results when considering large displacements [3]. Connections between members of these structures can affect critical loading and postbuckling behavior [4].
Usually in steel structure designs, the frames are analyzed with the simplification that the beam-column connection behavior can be idealized by two extreme cases: ideally flexible, where no moment is transmitted between the column and the beam and these elements behave independently, and perfectly rigid, in which the total transmission of the moment occurs [5]. However, experimental investigations on real structures have pointed out that most connections between structural elements should be treated as semirigid connections and moment-rotation curves are used to describe their behavior [6].
With the development of informatics and significant computational resources in the last twenty years, many researches about advanced analysis with semirigid connections were published worldwide. A numerical model that includes both nonlinear connection behavior and geometric nonlinearity of the structure was developed by Sekulovic and Salatic [7]. Pinheiro and Silveira [8] discussed numerical and computational strategies for nonlinear analysis of frames with semirigid connections. Zareian and Krawinkler [9] established a reliability-based approach for addressing the collapse potential of buildings. Lignos and Krawinkler [10] discussed the development of a database of experimental data of steel components and the use of this database for quantification of important parameters that affect the cyclic moment-rotation relationship at plastic hinge regions in beams. Nguyen and Kim [11] developed a numerical procedure based on the beam-column method for nonlinear elastic dynamic analysis of three-dimensional semirigid steel frames. Rungamornrat et al. [12] presented a numerical technique for analysis of two-dimensional frames accounting for both geometric nonlinearity and nonlinear elastic material behavior. To fully understand the seismic performances and failure modes of beam-column joints in RC buildings, a simplified analytical model of joint behavior was proposed by Bossio et al. [13]. A numerical approach based on the member discrete element method to investigate static and dynamic responses of steel structures with semirigid joints was presented by Ye and Xu [14]. Fernandes et al. [15] evaluated the geometrically nonlinear dynamic behavior of collapsed arches subjected to transverse force and plane frames with semirigid connections. e performance of base-isolated steel structures having special moment frames was assessed by Rezaei Rad and Banazadeh [16]. Alvarenga [17] proposed a new numerical formulation that includes the behavior of the semirigid connection into plastic-zone finite element model for steel plane frames.
An efficient methodology for solving the nonlinear system must be able to trace the complete equilibrium path and identify and pass through all limits or critical points of the structural system under analysis [18][19][20][21]. e nonlinear response of a structural system can be seen in Figure 1, in which a given displacement component may increase or decrease along the path. In this figure, load limit points (A, D), displacement limit points (B, C), and point of failure (E) are identified [20]. For the solution of the structural problem, we used an incremental-iterative procedure based on the Potra-Pták method [22], associated with the linear arc-length path-following technique [23]. e Potra-Pták method is a two-step method with cubic convergence order and consists of two evaluations of the system of nonlinear equations using the same Jacobian matrix in each iteration.
Considerable attention has been devoted to the problem of tracing the response of semirigid frame connections in the presence of geometric nonlinearity. In this context, this paper aims to present a numerical-computational model for static analysis of plane frames considering the geometric nonlinearity (hypothesis of large displacements and rotations, but small strains) and the semirigidity in the connections between structural members. All beam-column elements are assumed to remain linear elastic. Semirigid connections are simulated using the linear connection element proposed by Del Savio et al. [24] of null length, which considers the axial, tangential, and rotational stiffness. Structural systems are discretized by the finite element corotational method [25,26], whose formulation separates the rigid body movements from the deformations that are locally produced in the element.
Making an ideal connection is very difficult, impractical and is not economically justified. erefore, connections are more or less flexible or semirigid. For the great majority of the steel structures, the effects of the axial and shear forces in the deformation of the connection are small if compared with those caused by the bending moment [8,27]. For this reason, only the rotational stiffness of the connection is considered in the analyses. Problems of 2D frame structures with geometric nonlinearity and different types of connection (ideally pinned, semirigid, and perfectly rigid) are carried out with the free Scilab program-version 6.1.0 [28].
e results show that Potra-Pták procedure can decrease the number of iterations and the computing time in comparison with the standard and modified Newton-Raphson iterative schemes. Furthermore, the simulations show that the connection flexibility has a strong influence on the nonlinear behavior and stability of the structural systems studied. In fact, ideal considerations for beam-to-column joints can cause an inaccurate estimation of the response of frames since the real joints' behavior is between fully rigid and pinned connections.

Finite Element Corotational
Formulation. e finite element corotational formulation for 2D beam, presented by Crisfield [25] and Yaw [26], is described below. ere is no shear strain in the beam, and then the cross section remains plane and normal to its axis (Euler-Bernoulli theory). In the initial configuration, the coordinates of beam element nodes 1 and 2 in the global system are (X 1 , Y 1 ) and (X 2 , Y 2 ), respectively. e original length of the beam L 0 is given by the following equation [26]: For the beam element in its current configuration, the global nodal coordinates are (X 1 + u 1 , Y 1 + v 1 ) for node 1 and (X 2 + u 2 , Y 2 + v 2 ) for node 2, where u i is the displacement of the node i in the X-direction and v i is the displacement of node i in the Y-direction, with i � 1, 2. e current length L is e local axial displacement (u l ) of the element is calculated by 2 Mathematical Problems in Engineering We assumed the strain ε as constant and is determined by ε � u l /L 0 . e axial force (N) of the beam is then given by where A is the cross-sectional area and E is Young's modulus. Using standard structural analysis, the local moments at the ends of the beam element are related to the local nodal rotations (θ 1l and θ 2l ) and are given by [26] M 1 where I is the moment of inertia of the cross section. Local nodal rotations are computed by θ 1l � arctan cos β sin β 1 − sin β cos β 1 cos β cos β 1 + sin β sin β 1 , where β 1 � θ 1 + β 0 and β 2 � θ 2 + β 0 . e angles θ 1 , and θ 2 are the calculated global nodal rotations of the global equation system, and the expressions for the initial angle β 0 and current angle β of the bar are, respectively, e elemental tangent stiffness matrix K el is determined by the expression [26] where D is the constitutive matrix: where r � ��� I/A √ is the radius of gyration; the vectors z and r are, respectively, where s � sin(β) and c � cos(β), and matrix B is e following expressions are used to calculate the sine and cosine values of angle β, respectively: e elemental internal force vector (F el ) is determined by

Connection Element.
e semirigid connection is simulated by the connection element proposed by Del Savio et al. [24] whose length is null. In the stiffness matrix K lig are considered the axial (S a ), translational (S t ), and rotational (S r ) stiffness, which can be expressed mathematically by is element is inserted at the intersection points between frame members where the semirigid connections are located. It behaves appropriately for any type of loading and also allows to simulate elastoplastic analyses of connections, given the momentum-rotation curve (M-α curve) that describes connection behavior.
is model is linear and, therefore, the element stiffness matrix given by (15) is maintained invariable during the nonlinear analysis of the structure. e stiffness S a , S t , and S r can be obtained through experimental tests. In Figure 2, a schematic drawing of the connection element model is illustrated.

Structural Problem and Solution Method.
e basic problem of nonlinear analysis is to find the equilibrium configuration of a structure that is under the action of loads. e nonlinear equations' system to be solved that governs the static equilibrium of a structural system with geometric nonlinear behavior is given by [21] where u is the nodal displacement vector, g is the unbalanced force vector, F int is the global internal force vector, and λ is the load parameter responsible for scaling the arbitrary Mathematical Problems in Engineering magnitude reference vector F r . Equation (17) is a system of (n + 1) unknowns, with n displacement components (u) and a load parameter (λ), but only n equations. In this way, a constraint equation is added to the system [20,29]: e solution of the nonlinear system given by equations (17) and (18) can be obtained by the incremental and iterative scheme based on the Newton-Raphson method. e iterative equations are given by for every k � 0, 1, 2, . . ., and K(u (k) ) is the stiffness matrix representative of the structural system (Jacobian matrix). Note that the superscript k is used to refer to the previous iteration and (k + 1) the current iteration. e total load parameter λ (k+1) is updated by where δλ (k+1) is the load subincrement calculated according to equation (34). With the combination of equations (17), (20), and (21), it comes to the expression for δu (k+1) [25]: where the subincrements δu (k+1) r and δu (k+1) g are obtained, respectively, by Potra and Pták [22] developed a two-step method with cubic convergence based on the Newton-Raphson method for solving nonlinear equations of type f(x) � 0, and it consists of two evaluations of the function f(x) and the calculation of first-order derivatives f′(x) only [30]. e iterative equations for this method are given by e Potra and Pták method is adopted in this paper to solve the structural problem given by equations (17) and (18), obtaining the following iterative scheme [29]: e explicit calculation of [K(u (k) )] − 1 can be avoided by solving the system of linear equations via decomposition (e.g., LU factorization), since a single factorization at the beginning of the iteration is required. e incremental displacement parameter (Δu) at the load step t + Δt and iteration (k + 1) is evaluated by e load subincrement δλ (k+1) i , with i � 1, 2, is determined by the linear arc-length path-following technique [31,32]. e expression for the initial load increment (predicted solution) is given by (k � 0) where Δl represents the arc length increment. As proposed by Crisfield [25], the increment Δl can be used as a control parameter in the current load step according to the equation where 0 Δl represents the arc length increment at the initial load step (t � 0), Nd is the number of desired iterations for the current iterative process convergence, t k is the number of iterations that were required to converge at the previous load step, and the exponent δ � 0.5 [23,25]. Typical values for Nd can vary between 3 and 10 [25]. For the determination of the load subincrement correction (k > 0), the following expression is used [23]: e sign of the initial increment of the load parameter (Δλ (0) ) can be positive or negative. e correct choice of the signal is of paramount importance in defining the solution sequences (u, λ), which allow for continuous advancement in the load-displacement response. e procedure used consists of the analysis of the internal product between the incremental displacement obtained in the previous load step ( t Δu) and the current displacement increment δu r : if t Δu T δu r > 0, then the predictor Δu (0) has the same meaning as δu; otherwise, the predictor has the opposite meaning. Figure 3 shows the algorithm pseudo-code for the Potra-Pták method combined with the linear arc-length technique in an iterative incremental procedure. e input data in the algorithm are initial arc length ( 0 Δl); maximum number of iterations in each load step (k max ); number of desired iterations in each load step (Nd); tolerance (tol); load increment (ΔP); and maximum number of load steps (LS max ). e outputs of the algorithm are nodal displacement vector (u); total load parameter (λ); total number of load steps (LS); total number of accumulated iterations until convergence to the solution (k total ); the average number of iterations per load step (k av ); and CPU time in seconds (t).
In general, the user has no idea of the magnitude to consider for the initial increment of the arc length 0 Δl. To resolve this problem, we suggest the following expression as an initial estimate: where 0 δu r is calculated according to equation (23) at load step t � 0 (beginning of the analysis) and φ∈[0.001, 0.1]. For the next load steps, the increment Δl is automatically calculated by equation (33). e iterative process in the load step ends by indicating a new equilibrium position for the structural system under analysis when one or both of the convergence criteria used here are met at the end of each iterative cycle (see line 38 of the algorithm in Figure 3). e first criterion follows displacement relations and is described by where ‖•‖ is the Euclidean norm and tol is the tolerance provided by the user. e second criterion is based on the norm of the unbalanced force vector g and the norm of the incremental external force vector F r :

Results and Discussion
In this section, numerical analyses of 2D frame structures with geometric nonlinearity and different types of connection (ideally pinned, semirigid, and perfectly rigid) are performed. All beam-column elements are assumed to remain linear elastic, and the own weight of the structures is neglected. e algorithm presented in Figure 3 is implemented with Scilab software [28], and the computational simulations are performed on a Core i7-3537U computer with 8 GB of memory. Structure equilibrium paths are obtained and numerical results (LS, k total , k av , and t) are presented. At the end of the analysis, it was established as condition (or stopping criterion), instead of the maximum number of steps LS max , a maximum displacement referring to a degree of freedom of a given node of the structure. us, line 4 of the algorithm (see Figure 3) is modified. e connections are modeled with the connection element proposed by Del Savio et al. [24], being inserted at the point of intersection between the structural members (beam, column, or support). According to Chen et al. [33], axial, shear, flexion, and torsional forces are transmitted to a joint. Because the structures are plane, torsional deformation is neglected. However, by hypothesis, large numerical values are adopted for axial (S a ) and translational (S t ) stiffness, and only the value of rotational stiffness (S r ) is varied. In case of perfectly rigid connection, for rotational stiffness was considered a large numerical value. For ideally pinned connections, the value for S r is zero. e first two numerical problems (Williams frame and simple frame), considered as reference examples in the literature, are used to validate the connection element implemented in nonlinear geometric analyses. Figure 4, the axial stiffness of the beam is EA � 1.885 × 10 6 lb (�8.385 × 10 6 N) and flexural stiffness EI � 9.274 × 10 3 lb in 2 (�26.615 Nm 2 ). is frame is known as the Williams frame [34]. e structure was analyzed for three types of connection between the bar and the support: ideally pinned (S r � 0); semirigid (S r � 1.8 × 10 3 lb in �2.034 × 10 2 Nm); and perfectly rigid (S r � 1.0 × 10 10 lb in � 1.130 × 10 9 Nm). e following parameters were considered for the solution methods: 0 Δl � 0.03; k max � 150; Nd � 3; tolerance tol � 1.0 × 10 −6 ; and ΔP � 0.5 lb (�2.224 N). e finite element mesh consists of two connection elements and eight beam elements. Table 1 presents the numerical results (LS, k total , k av , and t) for the analysis with the semirigid connection, comparing the Potra-Pták (PP), modified Newton-Raphson (MNR), and standard Newton-Raphson (NR) solution methods. Figure 5 presents the equilibrium paths (vertical displacement in the central node v/h versus load P) for the three support conditions, and it shows the influence of the effect of semirigidity on boundary conditions, obtaining different equilibrium paths with two load limit points for the three types of connection. ere is reasonable agreement between the curves obtained here and the equilibrium points obtained by Tin-Loi and Misa [35].

Williams Frame. Considering the frame shown in
is type of analysis in structural system designs that exhibit nonlinear behavior and loss of stability by limit point followed by a snap-through is of utmost importance in verifying their safety. e frame deformed positions at the load limit points considering the semirigid connections are shown in Figure 6.

Simple Frame.
e simple frame shown in Figure 7 has semirigid connection (column-beam) and length L 0 � 3.524 m. e elements consist of the W 200 × 46 profile, whose cross section has the following geometric properties: area A � 5.89 × 10 −3 m 2 and moment of inertia I � 4.55 × 10 −5 m 4 . e modulus of elasticity is E � 200.0 GPa. is structure was studied by Chan and Chui [27] and Santos et al. [36]. e finite element mesh consists of 15 beam elements (five elements for each column and five elements for beam) and two connection elements between the beam and the columns. Two analyses are performed, and the connections between column and beam are considered perfectly rigid (S r � 1.0 × 10 15 Nm) and semirigid (S r � 10 EI/ L 0 ). e equilibrium paths are presented in Figure 8, and the results have good agreement with the equilibrium points obtained by Chan and Chui [27] and Santos et al. [36]. In the simulations with the solution methods, the following parameters were considered: 0 Δl � 0.01; k max � 150; Nd � 3; tol � 1.0 × 10 −5 ; and ΔP � 10 N. e numerical results LS, k total , k av , and t for simulations with the solution methods appear in Table 2.

Lee's Frame.
e frame illustrated in Figure 9 was studied by Lee et al. [37]. It is a classic example whose behavior is often described as strongly nonlinear. e bars have modulus of elasticity E � 720.0 kN/cm 2 , cross-sectional area A � 6.0 cm 2 , moment of inertia I � 2.0 cm 4 , and initial length L 0 � 120 cm.
e structure was discretized with 20 beam elements and two connection elements. Figure 10 shows the equilibrium paths with the following types of connection on the supports: perfectly rigid (S r � 1.0 × 10 8 kN cm); ideally pinned (S r � 0); and semirigid (S r � EI/L 0 and S r � 10 EI/L 0 ). In the paths with the ideally pinned and semirigid connections (S r � EI/L 0 ), there are two load limit points (the tangent at these points is parallel to the displacements axis) and two displacement limits points (the     tangent is vertical, parallel to the load axis). Good agreement is found between the ideally pinned connection curve obtained here and the equilibrium points obtained by Schweizerhof and Wriggers [38]. Table 3 presents the numerical results LS, k total , k av , and t for Lee's frame analysis. In the simulations with the solution methods, the following parameters were considered: 0 Δl � 4.0; k max � 150; Nd � 5; tol � 1.0 × 10 −5 ; and ΔP � 1.0 kN. Figure 11 plots the frame deformed positions at the load and displacement limit points, considering the semirigid connections (S r � EI/L 0 ).

Two-Story
Frame. e two-story frame with perfectly rigid supports (S r2 � 1.0 × 10 15 Nm) is shown in Figure 12.    columns: ideally pinned (S r1 � 0); perfectly rigid (S r1 � 1.0 × 10 15 Nm); and semirigid according to the types of connections illustrated in Figure 13. Figure 14 shows the equilibrium paths (horizontal displacement at node 9 versus load P). e structure was discretized by 28 beam finite elements and six connection elements. As it is expected, this figure shows that as the value of the rotational stiffness S r1 of the connection increases, the value of the collapse load of the structure increases. Table 4 shows the numerical results. In the simulations with the solution methods, the following parameters were considered: 0 Δl � 0.05; k max � 150; Nd � 3; tol � 1.0 × 10 −6 ; and ΔP � 1.0 N.

Analysis of Numerical Results
. It is noted that geometric nonlinearity becomes important with increasing the load, in particular for slender frames. e analysis of columns or beam-columns as independent members in framed structures may lead to erroneous results. e connection flexibility has significant influence on the behavior and stability of the frame structures studied. When the semirigid type is   used, the structure becomes more flexible decreasing the critical load and buckling capacity and has an intermediate behavior between the extreme connection conditions of totally rigid and ideally pinned. e need to use an iterative incremental method to solve geometric nonlinear problems is evident. Combined with the solution methods, it is necessary to include path-following techniques to obtain the complete equilibrium path with limit points. e Potra-Pták method combined with the linear arc-length technique is adequate for the studied problems. e incremental and iterative procedure based on this method shows better computational efficiency in processing time in all simulations. is method has already been employed by the authors Souza et al. [29] in nonlinear problems of plane and spatial trusses with geometric nonlinearity.
e higher computational cost of the Potra-Pták iteration, compared to the iterations of the standard Newton-Raphson (NR) and modified Newton-Raphson (NRM) methods, is compensated by reducing the number of load steps and required iterations for solution convergence. In the iterative procedure described by equations (26)-(30), three systems of linear equations are solved since the system [K(u (k) )] − 1 F r is solved only once in the iterative cycle (correction step), and two updates of the internal force vector F int are performed (see lines 29 and 36 of the algorithm in Figure 3).
In the NR procedure, two systems of linear equations are solved, and there is one update of the vector F int and the matrix K in the iterative cycle. e MNR procedure solves one system of linear equations and updates the vector F int once in the iterative cycle, and the matrix K is updated at the beginning of each load step remaining unchanged during the iterations. It is observed that for simulations with the NR and MNR methods, in the pseudo-algorithm of Figure 3 the vector δu 2 � 0 is considered; lines 29 to 33 are not executed and in line 37 δλ 2 is replaced by δλ 1 , since δλ 2 � 0. In the case of MNR, line 21 is also not executed.  e stiffness matrix generated from the finite element method is usually sparse and band-like (nonnull elements appear on the main diagonal and parallel diagonals). Large arrays require a large storage space, and even if there are computers with the largest memory capacity today, it is usually not enough to store the square array. To make array storage and operations less computationally costly, techniques based on nonzero value storage, such as compressed sparse column (CSC) and compressed sparse row (CSR), can be used. In the program developed in the Scilab environment, the "sparse" function was used, which stores the nonnull elements of the original matrix, neglecting the elements equal to zero.

Conclusion
Advanced structural analysis, including nonlinear effects, is a tool capable of capturing the strength and stability limit of a structure. e purpose of this paper was to present a computational-numerical model to solve 2D frame structures by geometric nonlinear analysis. A corotational analysis of the beam-column element with semirigid connections was settled. ree procedures were used to deal with the system of nonlinear equations arising from the geometrical nonlinear equations: NR, MNR, and Potra-Pták schemes. To trace the whole path described in the load x displacement space of the nonlinear response, the linear arc-length path-following technique was adopted. A computational code was developed with the free Scilab program, and the solution method pseudo-code was presented.
It was possible to verify the good agreement between the answers obtained in this paper and those found by other researchers. e results showed that the proposed model can predict the approximate mechanical behavior of semirigid frame structures with large displacements subjected to static loads. As expected, the connection flexibility between members of these structures affected the response and ultimate strength. erefore, semirigid connections may significantly reduce structural stiffness, and they need to be considered when analyzing the structure for a practical design.
Additionally, the Potra-Pták method was able to obtain the nonlinear responses of the studied problems with less computation time, by reducing the number of accumulated iterations until convergence, in comparison with the classic methods of standard and modified Newton-Raphson. e linear connection model was quite simple to use and easy to implement, since the initial stiffness of the connections remained constant throughout the analysis. However, a limitation of this model is that it is not precise in cases of large rotations or displacements. Its use is more appropriate for linear analyses in which the displacements are small. Based on numerous experimental results, the behavior of flexible connections, described through the M-α curve, is nonlinear in the domain for almost every connection type. e connection element, even if in this first implementation the effects of the stiffnesses are all decoupled, can be extended in future works that will take into account the interaction between the stiffness, just by modifying the stiffness matrix of this element. Like other directions for future works, we suggest the inclusion of the nonlinear material behavior; consider a nonlinear joint model by varying the rotational stiffness in the analysis; combine the Potra-Pták method with other path-following techniques (for example, cylindrical arclength, minimum residual displacement, and orthogonal residual); application of the model to 3D frame structures; and adapt the code for dynamic analysis.

Data Availability
e data used to support the findings of this study are included within the article.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.