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

Authors: Luiz Antonio Farani de Souza; Leandro Vanalli; Arthur Bueno de Luz
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

A numerical-computational model for static analysis of plane frames with semirigid connections and geometric nonlinear behavior is presented. The 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. The algorithm pseudo-code is presented, and the finite element corotational method is used for the discretization of the structures. The equilibrium paths with load and displacement limit points are obtained. The 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. The results show that the 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. Also, the simulations show that the connection flexibility has a strong influence on the nonlinear behavior and stability of the structural systems.


---

## Body

## 1. 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]. The 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. The 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–21]. The 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]. The 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.Figure 1
Equilibrium path with critical points. Source: [20].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. Therefore, 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]. The 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.


## 2. FE Formulation: Beam-Column and Connection

### 2.1. Finite Element Corotational Formulation

The finite element corotational formulation for 2D beam, presented by Crisfield [25] and Yaw [26], is described below. There 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 (X1, Y1) and (X2, Y2), respectively. The original length of the beam L0 is given by the following equation [26]:(1)L0=X2−X12+Y2−Y12.For the beam element in its current configuration, the global nodal coordinates are (X1 + u1, Y1 + v1) for node 1 and (X2 + u2, Y2 + v2) for node 2, where ui is the displacement of the node i in the X-direction and vi is the displacement of node i in the Y-direction, with i = 1, 2. The current length L is(2)L=X2+u2−X1−u12+Y2+v2−Y1−v12.The local axial displacement (ul) of the element is calculated by(3)ul=L2−L02L+L0.We assumed the strainε as constant and is determined by ε = ul/L0. The axial force (N) of the beam is then given by(4)N=EAulL0,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](5)M¯1M¯2=2EIL02112θ1lθ2l,where I is the moment of inertia of the cross section. Local nodal rotations are computed by(6)θ1l=arctancosβsinβ1−sinβcosβ1cosβcosβ1+sinβsinβ1,θ2l=arctancosβsinβ2−sinβcosβ2cosβcosβ2+sinβsinβ2,where β1 = θ1 + β0 and β2 = θ2 + β0. The 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,(7)β0=arctanY2−Y1X2−X1,β=arctanY2+v2−Y1−v1X2+u2−X1−u1.The elemental tangent stiffness matrixKel is determined by the expression [26](8)Kel=BTDB+NLzzT+M¯1+M¯2L2rzT+zrT,where D is the constitutive matrix:(9)D=EAL010004r22r202r24r2,where r=I/A is the radius of gyration; the vectors z and r are, respectively, (10)z=s−c0−sc0T,r=−c−s0c0T,where s = sin(β) and c = cos(β), and matrix B is(11)B=−c−s0cs0−sLcL1sL−cL0−sLcL0sL−cL1.The following expressions are used to calculate the sine and cosine values of angleβ, respectively:(12)sinβ=Y2+v2−Y1−v1L,(13)cosβ=X2+u2−X1−u1L.The elemental internal force vector (Fel) is determined by(14)Fel=BTNM¯1M¯2T.


### 2.2. Connection Element

The semirigid connection is simulated by the connection element proposed by Del Savio et al. [24] whose length is null. In the stiffness matrix Klig are considered the axial (Sa), translational (St), and rotational (Sr) stiffness, which can be expressed mathematically by(15)Klig=S−S−SS,with(16)S=Sa000St000Sr.This 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. This model is linear and, therefore, the element stiffness matrix given by (15) is maintained invariable during the nonlinear analysis of the structure. The stiffness Sa, St, and Sr can be obtained through experimental tests. In Figure 2, a schematic drawing of the connection element model is illustrated.Figure 2
Connection element model. Source: [24].


### 2.3. Structural Problem and Solution Method

The basic problem of nonlinear analysis is to find the equilibrium configuration of a structure that is under the action of loads. The nonlinear equations’ system to be solved that governs the static equilibrium of a structural system with geometric nonlinear behavior is given by [21](17)gu,λ=λFr−Fintu=0,where u is the nodal displacement vector, g is the unbalanced force vector, Fint is the global internal force vector, and λ is the load parameter responsible for scaling the arbitrary magnitude reference vector Fr. 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]:(18)cu,λ=0.The 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. The iterative equations are given by(19)uk+1=uk+δuk+1,(20)δuk+1=Kuk−1guk,λk+1=Kuk−1λk+1Fr−Fintuk,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. The total load parameter λ(k+1) is updated by(21)λk+1=λk+δλk+1,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 δuk+1 [25]:(22)δuk+1=Kuk−1λk+δλk+1Fr−Fintuk,δuk+1=Kuk−1δλk+1Fr+λkFr−Fintuk,δuk+1=δλk+1δurk+1+δugk+1,where the subincrements δurk+1 and δugk+1 are obtained, respectively, by(23)δurk+1=Kuk−1Fr,(24)δugk+1=Kuk−1guk,λk=Kuk−1λkFr−Fintuk.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]. The iterative equations for this method are given by(25)yk+1=xk−fxkf′xk,xk+1=xk−fxk+fyk+1f′xk.The 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]:(26)δu1k+1=Kuk−1δλ1k+1Fr+guk,λk,(27)u1k+1=uk+δu1k+1,(28)δu2k+1=Kuk−1δλ2k+1Fr+gu1k+1,λk,(29)uk+1=uk+δu1k+1+δu2k+1,(30)λk+1=λk+δλ2k+1.The explicit calculation ofKuk−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. The incremental displacement parameter (Δu) at the load step t + Δt and iteration (k + 1) is evaluated by(31)Δuk+1=Δuk+δu1k+1+δu2k+1.The load subincrementδλik+1, with i = 1, 2, is determined by the linear arc-length path-following technique [31, 32]. The expression for the initial load increment (predicted solution) is given by (k = 0)(32)Δλ0=Δlδur,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(33)Δl=Δ0lNdktδ,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, tk 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]:(34)δλik+1=−ΔukTδugk+1ΔukTδurk+1.The sign of the initial increment of the load parameter (Δλ(0)) can be positive or negative. The 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. The 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 δur: if tΔuTδur > 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.Figure 3
Pseudo-code for the Potra–Pták method combined with linear arc-length technique.The input data in the algorithm are initial arc length (0Δl); maximum number of iterations in each load step (kmax); number of desired iterations in each load step (Nd); tolerance (tol); load increment (ΔP); and maximum number of load steps (LSmax). The 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 (ktotal); the average number of iterations per load step (kav); 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 length0Δl. To resolve this problem, we suggest the following expression as an initial estimate:(35)Δ0l=φδ0urTδ0ur+FrTFr,where 0δur 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).The 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 Figure3). The first criterion follows displacement relations and is described by(36)δu2<tol⋅u,where ‖∙‖ is the Euclidean norm and tol is the tolerance provided by the user. The second criterion is based on the norm of the unbalanced force vector g and the norm of the incremental external force vector Fr:(37)g<tol⋅Fr.


## 2.1. Finite Element Corotational Formulation

The finite element corotational formulation for 2D beam, presented by Crisfield [25] and Yaw [26], is described below. There 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 (X1, Y1) and (X2, Y2), respectively. The original length of the beam L0 is given by the following equation [26]:(1)L0=X2−X12+Y2−Y12.For the beam element in its current configuration, the global nodal coordinates are (X1 + u1, Y1 + v1) for node 1 and (X2 + u2, Y2 + v2) for node 2, where ui is the displacement of the node i in the X-direction and vi is the displacement of node i in the Y-direction, with i = 1, 2. The current length L is(2)L=X2+u2−X1−u12+Y2+v2−Y1−v12.The local axial displacement (ul) of the element is calculated by(3)ul=L2−L02L+L0.We assumed the strainε as constant and is determined by ε = ul/L0. The axial force (N) of the beam is then given by(4)N=EAulL0,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](5)M¯1M¯2=2EIL02112θ1lθ2l,where I is the moment of inertia of the cross section. Local nodal rotations are computed by(6)θ1l=arctancosβsinβ1−sinβcosβ1cosβcosβ1+sinβsinβ1,θ2l=arctancosβsinβ2−sinβcosβ2cosβcosβ2+sinβsinβ2,where β1 = θ1 + β0 and β2 = θ2 + β0. The 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,(7)β0=arctanY2−Y1X2−X1,β=arctanY2+v2−Y1−v1X2+u2−X1−u1.The elemental tangent stiffness matrixKel is determined by the expression [26](8)Kel=BTDB+NLzzT+M¯1+M¯2L2rzT+zrT,where D is the constitutive matrix:(9)D=EAL010004r22r202r24r2,where r=I/A is the radius of gyration; the vectors z and r are, respectively, (10)z=s−c0−sc0T,r=−c−s0c0T,where s = sin(β) and c = cos(β), and matrix B is(11)B=−c−s0cs0−sLcL1sL−cL0−sLcL0sL−cL1.The following expressions are used to calculate the sine and cosine values of angleβ, respectively:(12)sinβ=Y2+v2−Y1−v1L,(13)cosβ=X2+u2−X1−u1L.The elemental internal force vector (Fel) is determined by(14)Fel=BTNM¯1M¯2T.


## 2.2. Connection Element

The semirigid connection is simulated by the connection element proposed by Del Savio et al. [24] whose length is null. In the stiffness matrix Klig are considered the axial (Sa), translational (St), and rotational (Sr) stiffness, which can be expressed mathematically by(15)Klig=S−S−SS,with(16)S=Sa000St000Sr.This 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. This model is linear and, therefore, the element stiffness matrix given by (15) is maintained invariable during the nonlinear analysis of the structure. The stiffness Sa, St, and Sr can be obtained through experimental tests. In Figure 2, a schematic drawing of the connection element model is illustrated.Figure 2
Connection element model. Source: [24].


## 2.3. Structural Problem and Solution Method

The basic problem of nonlinear analysis is to find the equilibrium configuration of a structure that is under the action of loads. The nonlinear equations’ system to be solved that governs the static equilibrium of a structural system with geometric nonlinear behavior is given by [21](17)gu,λ=λFr−Fintu=0,where u is the nodal displacement vector, g is the unbalanced force vector, Fint is the global internal force vector, and λ is the load parameter responsible for scaling the arbitrary magnitude reference vector Fr. 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]:(18)cu,λ=0.The 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. The iterative equations are given by(19)uk+1=uk+δuk+1,(20)δuk+1=Kuk−1guk,λk+1=Kuk−1λk+1Fr−Fintuk,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. The total load parameter λ(k+1) is updated by(21)λk+1=λk+δλk+1,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 δuk+1 [25]:(22)δuk+1=Kuk−1λk+δλk+1Fr−Fintuk,δuk+1=Kuk−1δλk+1Fr+λkFr−Fintuk,δuk+1=δλk+1δurk+1+δugk+1,where the subincrements δurk+1 and δugk+1 are obtained, respectively, by(23)δurk+1=Kuk−1Fr,(24)δugk+1=Kuk−1guk,λk=Kuk−1λkFr−Fintuk.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]. The iterative equations for this method are given by(25)yk+1=xk−fxkf′xk,xk+1=xk−fxk+fyk+1f′xk.The 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]:(26)δu1k+1=Kuk−1δλ1k+1Fr+guk,λk,(27)u1k+1=uk+δu1k+1,(28)δu2k+1=Kuk−1δλ2k+1Fr+gu1k+1,λk,(29)uk+1=uk+δu1k+1+δu2k+1,(30)λk+1=λk+δλ2k+1.The explicit calculation ofKuk−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. The incremental displacement parameter (Δu) at the load step t + Δt and iteration (k + 1) is evaluated by(31)Δuk+1=Δuk+δu1k+1+δu2k+1.The load subincrementδλik+1, with i = 1, 2, is determined by the linear arc-length path-following technique [31, 32]. The expression for the initial load increment (predicted solution) is given by (k = 0)(32)Δλ0=Δlδur,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(33)Δl=Δ0lNdktδ,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, tk 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]:(34)δλik+1=−ΔukTδugk+1ΔukTδurk+1.The sign of the initial increment of the load parameter (Δλ(0)) can be positive or negative. The 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. The 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 δur: if tΔuTδur > 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.Figure 3
Pseudo-code for the Potra–Pták method combined with linear arc-length technique.The input data in the algorithm are initial arc length (0Δl); maximum number of iterations in each load step (kmax); number of desired iterations in each load step (Nd); tolerance (tol); load increment (ΔP); and maximum number of load steps (LSmax). The 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 (ktotal); the average number of iterations per load step (kav); 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 length0Δl. To resolve this problem, we suggest the following expression as an initial estimate:(35)Δ0l=φδ0urTδ0ur+FrTFr,where 0δur 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).The 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 Figure3). The first criterion follows displacement relations and is described by(36)δu2<tol⋅u,where ‖∙‖ is the Euclidean norm and tol is the tolerance provided by the user. The second criterion is based on the norm of the unbalanced force vector g and the norm of the incremental external force vector Fr:(37)g<tol⋅Fr.


## 3. 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. The algorithm presented in Figure3 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, ktotal, kav, 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 LSmax, a maximum displacement referring to a degree of freedom of a given node of the structure. Thus, line 4 of the algorithm (see Figure 3) is modified.The 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 (Sa) and translational (St) stiffness, and only the value of rotational stiffness (Sr) is varied. In case of perfectly rigid connection, for rotational stiffness was considered a large numerical value. For ideally pinned connections, the value for Sr is zero. The 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.


### 3.1. Williams Frame

Considering the frame shown in Figure4, the axial stiffness of the beam is EA = 1.885 × 106 lb (≅8.385 × 106 N) and flexural stiffness EI = 9.274 × 103 lb in2 (≅26.615 Nm2). This frame is known as the Williams frame [34]. The structure was analyzed for three types of connection between the bar and the support: ideally pinned (Sr = 0); semirigid (Sr = 1.8 × 103 lb in ≅2.034 × 102 Nm); and perfectly rigid (Sr = 1.0 × 1010 lb in ≅ 1.130 × 109 Nm). The following parameters were considered for the solution methods: 0Δl = 0.03; kmax = 150; Nd = 3; tolerance tol = 1.0 × 10−6; and ΔP = 0.5 lb (≅2.224 N). The finite element mesh consists of two connection elements and eight beam elements. Table 1 presents the numerical results (LS, ktotal, kav, 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 4
Williams frame structural model.Table 1
Williams frame numerical results.
Solution method (semirigid)LSktotalkavt(s)PP27291.07410.175694NR36722.00000.208535MNR36722.00000.187576Figure5 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. There is reasonable agreement between the curves obtained here and the equilibrium points obtained by Tin-Loi and Misa [35]. This 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. The frame deformed positions at the load limit points considering the semirigid connections are shown in Figure 6.Figure 5
Williams frame equilibrium paths.Figure 6
Williams frame deformed positions at the limit points.


### 3.2. Simple Frame

The simple frame shown in Figure7 has semirigid connection (column–beam) and length L0 = 3.524 m. The elements consist of the W 200 × 46 profile, whose cross section has the following geometric properties: area A = 5.89 × 10−3 m2 and moment of inertia I = 4.55 × 10−5 m4. The modulus of elasticity is E = 200.0 GPa. This structure was studied by Chan and Chui [27] and Santos et al. [36]. The 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 (Sr = 1.0 × 1015 Nm) and semirigid (Sr = 10 EI/L0). The 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; kmax = 150; Nd = 3; tol = 1.0 × 10−5; and ΔP = 10 N. The numerical results LS, ktotal, kav, and t for simulations with the solution methods appear in Table 2.Figure 7
Simple frame structural model.Figure 8
Simple frame equilibrium paths.Table 2
Simple frame numerical results.
Solution method (semirigid)LSktotalkavt(s)PP27311.14810.230833NR37782.10810.291829MNR37912.45950.262597Solution method (perfectly rigid)LSktotalkavt(s)PP27311.14810.256468NR37772.08110.316001MNR37862.32430.274396


### 3.3. Lee’s Frame

The frame illustrated in Figure9 was studied by Lee et al. [37]. It is a classic example whose behavior is often described as strongly nonlinear. The bars have modulus of elasticity E = 720.0 kN/cm2, cross-sectional area A = 6.0 cm2, moment of inertia I = 2.0 cm4, and initial length L0 = 120 cm. The 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 (Sr = 1.0 × 108 kN cm); ideally pinned (Sr = 0); and semirigid (Sr = EI/L0 and Sr = 10 EI/L0). In the paths with the ideally pinned and semirigid connections (Sr = EI/L0), 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, ktotal, kav, and t for Lee’s frame analysis. In the simulations with the solution methods, the following parameters were considered: 0Δl = 4.0; kmax = 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 (Sr = EI/L0).Figure 9
Lee’s frame structural model.Figure 10
Lee’s frame equilibrium paths.Table 3
Lee’s frame numerical results.
Solution method (semirigidSr = EI/L0)LSktotalkavt(s)PP941902.02130.671181NR1143402.98250.781428MNR1264983.95240.684719Solution method (ideally pinned)LSktotalkavt(s)PP961952.03130.685775NR1173502.99150.804167MNR1335614.21800.715252Figure 11
Lee’s frame deformed positions at the limit points.


### 3.4. Two-Story Frame

The two-story frame with perfectly rigid supports (Sr2 = 1.0 × 1015 Nm) is shown in Figure 12. The beams consist of W 360 × 72 profile (A = 9.1 × 10−3 m2 and I = 2.01 × 10−4 m4) and the columns of W 310 × 143 profile (A = 0.0182 m2 and I = 3.48 × 10−4 m4). The modulus of elasticity E = 210.0 GPa is considered. Analyses were performed considering the connections between beams and columns: ideally pinned (Sr1 = 0); perfectly rigid (Sr1 = 1.0 × 1015 Nm); and semirigid according to the types of connections illustrated in Figure 13.Figure 12
Two-story frame equilibrium paths.Figure 13
Types of beam-column connections. Source: [39]. (a) Single-web angle. (b) Double-web angle. (c) Flush end-plate. (d) Extended end-plate.
(a)(b)(c)(d)Figure14 shows the equilibrium paths (horizontal displacement at node 9 versus load P). The 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 Sr1 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; kmax = 150; Nd = 3; tol = 1.0 × 10−6; and ΔP = 1.0 N.Figure 14
Two-story frame equilibrium paths.Table 4
Two-story frame numerical results.
Solution method (ideally pinned)LSktotalkavt(s)PP8111.37500.188026NR10242.40000.212376MNR10252.50000.197747Solution method (perfectly rigid)LSktotalkavt(s)PP11201.81820.219149NR13352.69230.253783MNR13413.15380.221970Solution method (single-web angle)LSktotalkavt(s)PP8121.50000.191959NR10252.50000.217752MNR11312.81820.208205Solution method (double-web angle)LSktotalkavt(s)PP9141.55560.199827NR11282.54550.230796MNR11312.81820.206055Solution method (flush end-plate)LSktotalkavt(s)PP9141.55560.198694NR11282.54550.235990MNR11333.00000.212012Solution method (extended end-plate)LSktotalkavt(s)PP9151.66670.211556NR12322.66670.256037MNR12352.91670.235058


### 3.5. Analysis of Numerical Results

It is noted that geometric nonlinearity becomes important with increasing the load, in particular for slender frames. The analysis of columns or beam-columns as independent members in framed structures may lead to erroneous results. The 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.The 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. The Potra–Pták method combined with the linear arc-length technique is adequate for the studied problems. The incremental and iterative procedure based on this method shows better computational efficiency in processing time in all simulations. This method has already been employed by the authors Souza et al. [29] in nonlinear problems of plane and spatial trusses with geometric nonlinearity.The 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 Kuk−1Fr is solved only once in the iterative cycle (correction step), and two updates of the internal force vector Fint 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 vectorFint and the matrix K in the iterative cycle. The MNR procedure solves one system of linear equations and updates the vector Fint 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 δu2 = 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.The 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.


## 3.1. Williams Frame

Considering the frame shown in Figure4, the axial stiffness of the beam is EA = 1.885 × 106 lb (≅8.385 × 106 N) and flexural stiffness EI = 9.274 × 103 lb in2 (≅26.615 Nm2). This frame is known as the Williams frame [34]. The structure was analyzed for three types of connection between the bar and the support: ideally pinned (Sr = 0); semirigid (Sr = 1.8 × 103 lb in ≅2.034 × 102 Nm); and perfectly rigid (Sr = 1.0 × 1010 lb in ≅ 1.130 × 109 Nm). The following parameters were considered for the solution methods: 0Δl = 0.03; kmax = 150; Nd = 3; tolerance tol = 1.0 × 10−6; and ΔP = 0.5 lb (≅2.224 N). The finite element mesh consists of two connection elements and eight beam elements. Table 1 presents the numerical results (LS, ktotal, kav, 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 4
Williams frame structural model.Table 1
Williams frame numerical results.
Solution method (semirigid)LSktotalkavt(s)PP27291.07410.175694NR36722.00000.208535MNR36722.00000.187576Figure5 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. There is reasonable agreement between the curves obtained here and the equilibrium points obtained by Tin-Loi and Misa [35]. This 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. The frame deformed positions at the load limit points considering the semirigid connections are shown in Figure 6.Figure 5
Williams frame equilibrium paths.Figure 6
Williams frame deformed positions at the limit points.


## 3.2. Simple Frame

The simple frame shown in Figure7 has semirigid connection (column–beam) and length L0 = 3.524 m. The elements consist of the W 200 × 46 profile, whose cross section has the following geometric properties: area A = 5.89 × 10−3 m2 and moment of inertia I = 4.55 × 10−5 m4. The modulus of elasticity is E = 200.0 GPa. This structure was studied by Chan and Chui [27] and Santos et al. [36]. The 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 (Sr = 1.0 × 1015 Nm) and semirigid (Sr = 10 EI/L0). The 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; kmax = 150; Nd = 3; tol = 1.0 × 10−5; and ΔP = 10 N. The numerical results LS, ktotal, kav, and t for simulations with the solution methods appear in Table 2.Figure 7
Simple frame structural model.Figure 8
Simple frame equilibrium paths.Table 2
Simple frame numerical results.
Solution method (semirigid)LSktotalkavt(s)PP27311.14810.230833NR37782.10810.291829MNR37912.45950.262597Solution method (perfectly rigid)LSktotalkavt(s)PP27311.14810.256468NR37772.08110.316001MNR37862.32430.274396


## 3.3. Lee’s Frame

The frame illustrated in Figure9 was studied by Lee et al. [37]. It is a classic example whose behavior is often described as strongly nonlinear. The bars have modulus of elasticity E = 720.0 kN/cm2, cross-sectional area A = 6.0 cm2, moment of inertia I = 2.0 cm4, and initial length L0 = 120 cm. The 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 (Sr = 1.0 × 108 kN cm); ideally pinned (Sr = 0); and semirigid (Sr = EI/L0 and Sr = 10 EI/L0). In the paths with the ideally pinned and semirigid connections (Sr = EI/L0), 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, ktotal, kav, and t for Lee’s frame analysis. In the simulations with the solution methods, the following parameters were considered: 0Δl = 4.0; kmax = 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 (Sr = EI/L0).Figure 9
Lee’s frame structural model.Figure 10
Lee’s frame equilibrium paths.Table 3
Lee’s frame numerical results.
Solution method (semirigidSr = EI/L0)LSktotalkavt(s)PP941902.02130.671181NR1143402.98250.781428MNR1264983.95240.684719Solution method (ideally pinned)LSktotalkavt(s)PP961952.03130.685775NR1173502.99150.804167MNR1335614.21800.715252Figure 11
Lee’s frame deformed positions at the limit points.


## 3.4. Two-Story Frame

The two-story frame with perfectly rigid supports (Sr2 = 1.0 × 1015 Nm) is shown in Figure 12. The beams consist of W 360 × 72 profile (A = 9.1 × 10−3 m2 and I = 2.01 × 10−4 m4) and the columns of W 310 × 143 profile (A = 0.0182 m2 and I = 3.48 × 10−4 m4). The modulus of elasticity E = 210.0 GPa is considered. Analyses were performed considering the connections between beams and columns: ideally pinned (Sr1 = 0); perfectly rigid (Sr1 = 1.0 × 1015 Nm); and semirigid according to the types of connections illustrated in Figure 13.Figure 12
Two-story frame equilibrium paths.Figure 13
Types of beam-column connections. Source: [39]. (a) Single-web angle. (b) Double-web angle. (c) Flush end-plate. (d) Extended end-plate.
(a)(b)(c)(d)Figure14 shows the equilibrium paths (horizontal displacement at node 9 versus load P). The 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 Sr1 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; kmax = 150; Nd = 3; tol = 1.0 × 10−6; and ΔP = 1.0 N.Figure 14
Two-story frame equilibrium paths.Table 4
Two-story frame numerical results.
Solution method (ideally pinned)LSktotalkavt(s)PP8111.37500.188026NR10242.40000.212376MNR10252.50000.197747Solution method (perfectly rigid)LSktotalkavt(s)PP11201.81820.219149NR13352.69230.253783MNR13413.15380.221970Solution method (single-web angle)LSktotalkavt(s)PP8121.50000.191959NR10252.50000.217752MNR11312.81820.208205Solution method (double-web angle)LSktotalkavt(s)PP9141.55560.199827NR11282.54550.230796MNR11312.81820.206055Solution method (flush end-plate)LSktotalkavt(s)PP9141.55560.198694NR11282.54550.235990MNR11333.00000.212012Solution method (extended end-plate)LSktotalkavt(s)PP9151.66670.211556NR12322.66670.256037MNR12352.91670.235058


## 3.5. Analysis of Numerical Results

It is noted that geometric nonlinearity becomes important with increasing the load, in particular for slender frames. The analysis of columns or beam-columns as independent members in framed structures may lead to erroneous results. The 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.The 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. The Potra–Pták method combined with the linear arc-length technique is adequate for the studied problems. The incremental and iterative procedure based on this method shows better computational efficiency in processing time in all simulations. This method has already been employed by the authors Souza et al. [29] in nonlinear problems of plane and spatial trusses with geometric nonlinearity.The 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 Kuk−1Fr is solved only once in the iterative cycle (correction step), and two updates of the internal force vector Fint 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 vectorFint and the matrix K in the iterative cycle. The MNR procedure solves one system of linear equations and updates the vector Fint 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 δu2 = 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.The 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.


## 4. Conclusion

Advanced structural analysis, including nonlinear effects, is a tool capable of capturing the strength and stability limit of a structure. The 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. Three 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 loadx 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. The 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. Therefore, 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.The 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.The 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 arc-length, minimum residual displacement, and orthogonal residual); application of the model to 3D frame structures; and adapt the code for dynamic analysis.



---

*Source: 3613892-2020-09-03.xml*