Nonlinear numerical model of plane frames considering semi-rigid connection and different beam theories

Modelo numérico não linear de pórticos planos considerando ligação semirrígida e diferentes teorias de vigas

ABSTRACT

This paper presents a numerical-computational model for frames with geometric nonlinear behavior, by the Finite Element Co-rotational method, considering the Euler-Bernoulli and Timoshenko beam theories. The connection between structural members is simulated by a null-length connection element, which considers the axial, translational and rotational stiffness. The nonlinear equations system that describes the structural problem is solved by the incremental and iterative procedure of Potra-Pták, with cubic convergence order, combined with the Linear Arc-Length path-following technique. The solution method algorithm is presented and the numerical examples are simulated with the free Scilab program. The numerical results show that the slenderness of the structure, geometric nonlinearity and semi-rigidity influence the behavior of the structure. Structural analysis and design procedures that consider these factors attains less conservative design thus obtaining more optimized structures.

Keywords
Shear deformation; Potra-Pták; Connection element

RESUMO

Este artigo apresenta um modelo numérico-computacional para pórticos com comportamento não linear geométrico, por meio do método Corrotacional de Elementos Finitos, considerando as teorias de viga de Euler-Bernoulli e Timoshenko. A ligação entre os membros estruturais é simulada por um elemento de ligação de comprimento nulo, o qual considera as rigidezes axial, translacional e rotacional. O sistema de equações não lineares que descreve o problema estrutural é solucionado pelo procedimento incremental e iterativo de Potra-Pták, com ordem de convergência cúbica, associado à técnica de continuação Comprimento de Arco Linear. O algoritmo do método de solução é apresentado e os exemplos numéricos são simulados com o programa livre Scilab. Os resultados numéricos mostram que a esbelteza da estrutura, a não linearidade geométrica e a semirrigidez influenciam no comportamento da estrutura. A análise estrutural e os procedimentos de projeto que consideram esses fatores alcançam um projeto menos conservador, obtendo-se estruturas mais otimizadas.

Palavras-chave
Deformação de cisalhamento; Potra-Pták; Elemento de ligação

1. INTRODUCTION

Light weighted structures are widely used in engineering applications such as buildings, industrial constructions, off-shore platforms and aerospace structures. Linear models are generally used for the design of these structures. However, as its slenderness increases, the degree of geometric nonlinearity increases as well as the occurrence of instability phenomena [¹].

Structural stability analysis using the Finite Element Method (FEM) involves solving a system of nonlinear equations. Methods that combine incremental and iterative procedures are used to obtain the system solution. An efficient methodology must be able to trace the complete equilibrium path and identify and pass through all the existing singular or critical points of the structure under analysis [², ³].

The key point in geometrically nonlinear analysis boils down to the description of the kinematics. There are three kinematics description formulations in literature for a flexible beam element, which are: total Lagrangian; updated Lagrangian; and Co-Rotational [⁴]. The two-dimensional co-rotational formulation for beams was developed by CRISFIELD [⁵] and is based on the explicit separation of rigid body movements (translations and rotations) and deformational movements [⁶]. The co-rotational formulation is substantiated in theory of infinitesimal strains. It has two coordinate systems, namely: the global coordinate fixed system; and the local co-rotational coordinate system.

Two theories have been used to study deflections of beams and frames, which are the Euler-Bernoulli and Timoshenko theories. In the Euler-Bernoulli theory, the cross-sectional initially normal to the beam neutral axis remains plane and normal to this axis after deformation [⁷]. This assumption ignores the effects of shear deformation and presents good results for isotropic homogeneous thin beams [⁸].

The Timoshenko beam theory, on the other hand, employs a more accurate representation of beam flexure which allows for the inclusion of shear strains and is, therefore, more suited for thick beam analysis [⁹]. It is considered that the cross-sectional initially normal to center axis of the beam remains plane, but not necessarily normal to the axis after deformation [⁷].

Standard analysis procedures for frame are usually based on the assumption of the ideal nodal connections. Two extreme idealizations for models with ideal connections are used: ideally pinned, no moment transfer; and perfectly rigid, with total moment transfer [¹⁰, ¹¹]. Models with ideal connections between beam and column simplify the analysis procedure, but often lead to an incorrect estimate of the real structural behavior. In fact, the existing connections are between the two extreme assumptions and have some rotational rigidity [¹², ¹³].

Researchers have studied the semi-rigid connections of beam-column through numerical models and experimental investigations. NGO-HUU et al. [¹⁴] developed an independent zero-length connection element with six different translational and rotational springs connecting two identical nodes to simulate the beam-to-column connections. HADADE et al. [¹⁵] presented the results of an experimental investigation on the semi-rigid behavior of typical beam-column connections for precast concrete structures. MOKHTAR et al. [¹⁶] investigated the performance of a specific beam-to-column connection using partly hidden corbels with different anchorage arrangements in precast concrete structures. WANG et al. [¹⁷] presented experimental efforts in steel framed substructure with bolted flange plate connections subjected to progressive collapse. A simplified modeling approach based on the shell element and ductile fracture criterion was introduced.

Nonlinear analysis of structures, including any kind of flexible connection, takes on special importance in engineering analysis. From this reasoning, it is necessary to develop a well-posed nonlinear formulation allowing the accurate evaluation of displacements and efforts of conventional and unconventional structures [¹⁸]. In this context, this paper presents a numerical-computational model for the analysis of plane frames with geometric nonlinear behavior (hypothesis of large displacements and rotations, but small strains), by the Finite Element Co-Rotational method.

The materials are considered isotropic and linear elastic and the cross section of the elements is uniform. Two beam-column models are implemented: one based on the Euler-Bernoulli theory [⁶]; and the other, on the Timoshenko theory [¹⁹]. The finite element formulation of these models - the stiffness matrix and the internal forces vector - is presented. The connection is simulated by the null-length connection element proposed by DEL SAVIO et al. [²⁰], in which it considers axial, translational and rotational rigidities. This element allows to carry out elastoplastic analyzes of the connections and, even, of the structure bars simulating the formation of plastic hinges within each finite element of the bar.

The system of nonlinear equations that describes the structural problem is solved by the iterative scheme of POTRA-PTÁK [²¹], proposed by SOUZA et al. [²²], combined with the Linear Arc-Length path-following technique (updated hyperplane constraint). The computational code was developed with the free Scilab program, version 6.1.0 [²³]. Numerical simulations of nonlinear problems found in the literature (cantilever beam, Lee’s frame, semi-circular arch and Williams' frame) are studied, in order to verify the formulation of the implemented beam-column models and the connection element. The behavior of structural systems is evaluated by their equilibrium path, and it was found that the shear effect can be significant on the mechanical response of these, depending on the slenderness of the structural elements. Furthermore, in the example of the Williams' frame, it was evident that the connection type can have a great influence on the overall structure behavior.

2. 2D BEAM-COLUMN CO-ROTATIONAL FINITE ELEMENT

Geometric nonlinearity is considered by means of co-rotational formulation. This kinematic description aims to seek to explicitly separate rigid body motions from strain producing deformations at the local element level [⁵, ⁶]. It is assumed that there is no shear deformation in the beam and, therefore, the cross section remains plane and normal to its axis. Consider a typical beam-column element in its initial and current configurations as shown in Figure 1.

Figure 1
Initial and current configuration for typical beam-column element.

The global coordinates of nodes 1 and 2 of the element in the initial configuration are (X₁, Y₁) and (X₂, Y₂), respectively. The original length L₀ is given by the following equation [⁶]:

$${\mathrm{L}}_0=\sqrt{{({\mathrm{X}}_2-{\mathrm{X}}_1)}^2+{({\mathrm{Y}}_2-{\mathrm{Y}}_1)}^2}$$ (1)

For the beam-column element in its current configuration, the global nodal coordinates are (X₁ + u₁, Y₁ + v₁) for node 1 and (X₂ + u₂, Y₂ + v₂) 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. The current length L of the beam element is [⁶]:

$$\mathrm{L}=\sqrt{{({\mathrm{X}}_2+{\mathrm{u}}_2-{\mathrm{X}}_1-{\mathrm{u}}_1)}^2+{({\mathrm{Y}}_2+{\mathrm{v}}_2-{\mathrm{Y}}_1-{\mathrm{v}}_1)}^2}$$ (2)

The global displacement vector p of element is:

$${\mathbf{p}}_{\mathbf{i}}={\left[\begin{array}{llllll}{\mathrm{u}}_1& {\mathrm{v}}_1& {\mathrm{\theta }}_1& {\mathrm{u}}_2& {\mathrm{v}}_2& {\mathrm{\theta }}_2\end{array}\right]}^{\mathrm{T}}$$ (3)

The local axial displacement (u) of the element is calculated as:

$$\overline{\mathrm{u}}=\mathrm{L}-{\mathrm{L}}_0$$ (4)

If the difference between L and L₀ is small, the Equation (4) is ill-conditioned for use in a numerical setting. Therefore, CRISFIELD [⁵] proposes the multiplication of by (L + L₀) / (L + L₀) which gives better conditioned formula:

$$\overline{\mathrm{u}}=\frac{(\mathrm{L}-{\mathrm{L}}_0)(\mathrm{L}+{\mathrm{L}}_0)}{\mathrm{L}+{\mathrm{L}}_0}=\frac{{\mathrm{L}}^2-{{\mathrm{L}}_0}^2}{\mathrm{L}+{\mathrm{L}}_0}$$ (5)

The strain e is assumed to be constant and is determined by ε = u/L₀. The axial force (N) of the element is then given by:

$$N=\frac{EA\overline{u}}{L_0}=\frac{EA(L^2-{L_0}^2)}{L{L}_0+{L_0}^2}=\frac{\mathrm{EA}[{\left(\frac{L}{L_0}\right)}^2-1]}{\frac{\mathrm{L}}{{\mathrm{L}}_0}+1}$$ (6)

where A is the cross-sectional area and E is Young's modulus. Using standard structural analysis, the local end moments of the element (M₁ and M₂) are related to the arbitrary local nodal rotations (θ_1l and θ_2l) as follows:

$$\left[\begin{array}{l}{\overline{\mathrm{M}}}_1\\ {\overline{\mathrm{M}}}_2\end{array}\right]=\frac{2\mathrm{EI}}{{\mathrm{L}}_0}\left[\begin{array}{ll}2& 1\\ 1& 2\end{array}\right]\left[\begin{array}{l}{\theta }_{11}\\ {\theta }_{21}\end{array}\right],$$ (7)

where I is the moment of inertia of the cross-section. The rotations θ_1l and θ_2l are computed by, respectively:

$${\theta }_{11}=\arctan \left(\frac{\cos \beta \sin {\beta }_1-\sin \beta \cos {\beta }_1}{\cos \beta \cos {\beta }_1+\sin \beta \sin {\beta }_1}\right)$$ (8)

$${\theta }_{21}=\arctan \left(\frac{\cos \beta \sin {\beta }_2-\sin \beta \cos {\beta }_2}{\cos \beta \cos {\beta }_2+\sin \beta \sin {\beta }_2}\right)$$ (9)

where β₁ = θ₁ + β₀ and β₂ = θ₂ + β₀. The angles θ₁ and θ₂ are the global nodal rotations calculated from the global matrix equations. The initial angle β₀ and current angle b of the element are calculated by, respectively:

$${\beta }_0=\arctan \left(\frac{Y_2-Y_1}{X_2-X_1}\right)$$ (10)

$$\beta =\arctan \left(\frac{Y_2+v_2-Y_1-v_1}{X_2+u_2-X_1-u_1}\right)$$ (11)

The tangent stiffness matrix K_el is determined as a function of the material stiffness matrix K_M and geometric stiffness matrix K_G:

$${\mathbf{K}}_{\mathbf{el}}={\mathbf{K}}_{\mathbf{M}}+{\mathbf{K}}_{\mathbf{G}}$$ (12)

where

$${\mathbf{K}}_{\mathbf{M}}={\mathbf{B}}^{\mathrm{T}}\mathbf{DB}$$ (13)

$${\mathbf{K}}_{\mathbf{G}}=\frac{\mathrm{N}}{\mathrm{L}}{\mathbf{zz}}^{\mathrm{T}}+\frac{{\overline{\mathrm{M}}}_1+{\overline{\mathrm{M}}}_2}{{\mathrm{L}}^2}({\mathbf{rz}}^{\mathrm{T}}+{\mathbf{zr}}^{\mathrm{T}})$$ (14)

The matrix Δ is the constitutive matrix given by:

$$\mathbf{D}=\frac{\mathrm{EA}}{{\mathrm{L}}_0}\left[\begin{array}{ccc}1& 0& 0\\ 0& 4{\mathrm{r}}^2& 2{\mathrm{r}}^2\\ 0& 2{\mathrm{r}}^2& 4{\mathrm{r}}^2\end{array}\right]$$ (15)

where $\mathrm{r}=\sqrt{\mathrm{I}/\mathrm{A}}$ is the radius of gyration and the vectors z and r are, respectively:

$$\mathbf{z}={\left[\begin{array}{llllll}\mathrm{s}& -\mathrm{c}& 0& -\mathrm{s}& \mathrm{c}& 0\end{array}\right]}^{\mathrm{T}},$$ (16)

$$\mathbf{r}={\left[\begin{array}{cccccc}-\mathrm{c}& -\mathrm{s}& 0& \mathrm{c}& \mathrm{s}& 0\end{array}\right]}^{\mathrm{T}},$$ (17)

where s = sin(β) and c = cos(β). The matrix B is:

$$\mathbf{B}=\left[\begin{array}{cccccc}-\mathrm{c}& -\mathrm{s}& 0& \mathrm{c}& \mathrm{s}& 0\\ -\mathrm{s}/\mathrm{L}& \mathrm{c}/\mathrm{L}& 1& \mathrm{s}/\mathrm{L}& -\mathrm{c}/\mathrm{L}& 0\\ -\mathrm{s}/\mathrm{L}& \mathrm{c}/\mathrm{L}& 0& \mathrm{s}/\mathrm{L}& -\mathrm{c}/\mathrm{L}& 1\end{array}\right]$$ (18)

The following equations are used to calculate the sine and cosine values of angle β, respectively:

$$\cos \left(\beta \right)=\frac{X_2+u_2-X_1-u_1}{L},$$ (19)

$$\sin \left(\beta \right)=\frac{Y_2+v_2-Y_1-v_1}{L}$$ (20)

The internal forces vector (F_el) is determined by:

$${\mathbf{F}}_{\mathbf{el}}={\mathbf{B}}^{\mathrm{T}}{\left[\begin{array}{lll}\mathbf{N}& {\overline{\mathrm{M}}}_1& {\overline{\mathrm{M}}}_2\end{array}\right]}^{\mathrm{T}}$$ (21)

A classic Timoshenko beam-column element with two nodes is defined with linear interpolations for u, v and θ in the local co-rotational coordinate system by, respectively [¹⁹]:

$$\mathrm{u}=\frac{\mathrm{x}}{{\mathrm{L}}_0}{\mathrm{u}}_{\mathrm{l}}$$ (22)

$$\mathrm{v}=0$$ (23)

$$\theta =(1-\frac{\mathrm{x}}{{\mathrm{L}}_0}){\theta }_{1\mathrm{l}}+\frac{\mathrm{x}}{{\mathrm{L}}_0}{\theta }_{21}$$ (24)

The curvature κ, the shear deformation g and the strain e are defined by, respectively:

$$\kappa =\frac{\mathrm{\partial }\theta }{\mathrm{\partial }\mathrm{x}}=-\frac{{\theta }_{11}}{{\mathrm{L}}_0}+\frac{{\theta }_{21}}{{\mathrm{L}}_0}=\frac{{\theta }_{21}-{\theta }_{11}}{{\mathrm{L}}_0}$$ (25)

$$\gamma =\frac{\mathrm{\partial }\mathrm{v}}{\mathrm{\partial }\mathrm{x}}-\theta =-(1-\frac{\mathrm{x}}{{\mathrm{L}}_0}){\theta }_{11}-\frac{\mathrm{x}}{{\mathrm{L}}_0}{\theta }_{21}$$ (26)

$$\epsilon =\frac{\mathrm{\partial }\mathrm{u}}{\mathrm{\partial }\mathrm{x}}-\kappa \mathrm{y}=\frac{{\mathrm{u}}_1}{{\mathrm{L}}_0}-\frac{{\theta }_{21}-{\theta }_{11}}{{\mathrm{L}}_0}\mathrm{y}$$ (27)

For calculating the material stiffness matrix K_M, the constitutive matrix Δ is expressed as:

$$\mathbf{D}=\frac{1}{4{\mathrm{L}}_0}\left[\begin{array}{ccc}4\mathrm{EA}& 0& 0\\ 0& 4\mathrm{EI}+\mu {{\mathrm{GAL}}_0}^2& -4\mathrm{EI}+\mu {{\mathrm{GAL}}_0}^2\\ 0& -4\mathrm{EI}+\mu {{\mathrm{GAL}}_0}^2& 4\mathrm{EI}+\mu {{\mathrm{GAL}}_0}^2\end{array}\right]$$ (28)

where G is the transverse shear modulus and parameter m is the sectional shear correction coefficient. The internal forces vector F_el is determined by the Equation (21) and the local end moments of the element M₁ and M₂ are, respectively:

$$\left[\begin{array}{l}{\overline{\mathrm{M}}}_1\\ {\overline{\mathrm{M}}}_2\end{array}\right]=\left[\begin{array}{l}\frac{\mathrm{EI}}{{\mathrm{L}}_0}({\theta }_{11}-{\theta }_{21})+\frac{1}{4}\mu {\mathrm{GAL}}_0({\theta }_{11}+{\theta }_{21})\\ \frac{\mathrm{EI}}{{\mathrm{L}}_0}({\theta }_{21}-{\theta }_{11})+\frac{1}{4}\mu {\mathrm{GAL}}_0({\theta }_{11}+{\theta }_{21})\end{array}\right]$$ (29)

and the axial force N is calculated by Equation (6).

3. CONNECTION ELEMENT

The connection is simulated by the element proposed by DEL SAVIO et al. [²⁰] whose length is null. The axial (S_a), translational (S_t) and rotational (S_r) rigidities are considered in the stiffness matrix K_lig, which can be expressed by:

$${\mathbf{K}}_{\mathrm{lig}}=\left[\begin{array}{cc}\mathbf{S}& -\mathbf{S}\\ -\mathbf{S}& \mathbf{S}\end{array}\right],$$ (30)

with

$$\mathbf{S}=\left[\begin{array}{ccc}{\mathrm{S}}_{\mathrm{a}}& 0& 0\\ 0& {\mathrm{S}}_{\mathrm{t}}& 0\\ 0& 0& {\mathrm{S}}_{\mathrm{r}}\end{array}\right]$$ (31)

The connection element is inserted at the intersection points between structure members where the connections are located. This element behaves appropriately for any type of loading and also allows to simulate elastoplastic analyzes of connections, given the momentum-rotation curve that describes the connection behavior.

4. STRUCTURAL PROBLEM AND SOLUTION METHOD

The basic objective of the nonlinear analysis is to find the equilibrium configuration of a structure that is under the loads action. The equation that governs the static equilibrium of a structural system with geometrically nonlinear behavior can be written by [³]:

$$\mathbf{g}(\mathbf{u},\lambda )=\lambda {\mathbf{F}}_{\mathbf{r}}-{\mathbf{F}}_{\mathbf{int}}\left(\mathbf{u}\right)=\mathbf{0}$$ (32)

where g is the unbalanced forces vector, u is the nodal displacement vector, F_int is the internal forces vector, λ is the load parameter and F_r is a reference vector characterizing the external load direction.

POTRA and PTÁK [²¹] developed a two-step method based on the Newton-Raphson method for the solution of systems of nonlinear equations, and consists of two evaluations of the function and only the calculation of first order derivatives. With cubic convergence and an efficiency index higher than that of the Newton-Raphson method [²⁴], the iterative scheme proposed by these authors was adapted in this paper to the structural problem whose equations are [²²]:

$${\mathbf{y}}^{(\mathrm{k}+1)}={\mathbf{u}}^{\left(\mathrm{k}\right)}+\mathit{\delta }{\mathbf{u}}_1^{(\mathrm{k}+1)}$$ (33)

$$\mathit{\delta }{\mathbf{u}}_1^{(\mathrm{k}+1)}={\left[\mathbf{K}\left({\mathbf{u}}^{\left(\mathrm{k}\right)}\right)\right]}^{-1}[\delta {\lambda }_1^{(\mathrm{k}+1)}{\mathbf{F}}_{\mathbf{r}}+\mathbf{g}({\mathbf{u}}^{\left(\mathrm{k}\right)},{\lambda }^{\left(\mathrm{k}\right)})]$$ (34)

$$\mathit{\delta }{\mathbf{u}}_2^{(\mathrm{k}+1)}={\left[\mathbf{K}\left({\mathbf{u}}^{\left(\mathrm{k}\right)}\right)\right]}^{-1}[\delta {\lambda }_2^{(\mathrm{k}+1)}{\mathbf{F}}_{\mathbf{r}}+\mathbf{g}({\mathbf{y}}^{(\mathrm{k}+1)},{\lambda }^{\left(\mathrm{k}\right)})]$$ (35)

$${\mathbf{u}}^{(\mathrm{k}+1)}={\mathbf{u}}^{\left(\mathrm{k}\right)}+\mathit{\delta }{\mathbf{u}}_1^{(\mathrm{k}+1)}+\mathit{\delta }{\mathbf{u}}_2^{(\mathrm{k}+1)}$$ (36)

$${\lambda }^{(\mathrm{k}+1)}={\lambda }^{\left(\mathrm{k}\right)}+\delta {\lambda }_2^{(\mathrm{k}+1)},$$ (37)

for all k = 0, 1, 2, ..., and K is the stiffness matrix representative of the structural system (Jacobian matrix). The load sub-increment δλ_i^(k+1), with i = 1, 2, is determined according to the Linear Arc-Length path-following technique [²⁵, ²⁶]. The expression for the initial load increment (predicted solution) is given by (k = 0):

$$\mathrm{\Delta }{\lambda }^{\left(0\right)}=\frac{\mathrm{\Delta l}}{\parallel \delta {\mathbf{u}}_{\mathrm{r}}{\parallel }^{\mathrm{\prime }}}$$ (38)

where Δl represents the arc length increment. Proposed by CRISFIELD [⁵], the increment Δl can be used as a control parameter in the current load step according to the equation:

$$\mathrm{\Delta l}{=}^0\mathrm{\Delta l}{\left(\frac{{\mathrm{k}}_{\mathrm{d}}}{{\mathrm{t}}_{\mathrm{k}}}\right)}^{0.5}$$ (39)

where ⁰Δl represents the arc length increment at the initial load step t = 0, k_d is the desired number of iterations for the convergence of the current iterative process and ^tk is the number of iterations needed to converge at the previous load step. To determine the correction of the load sub-increment (k > 0), the following expression is used [²⁷]:

$$\delta {\lambda }_i=-\frac{\mathrm{\Delta }{\mathbf{u}}^{\mathrm{T}}\delta {\mathbf{u}}_{\mathbf{g}}}{\mathrm{\Delta }{\mathbf{u}}^{\mathrm{T}}\delta {\mathbf{u}}_{\mathbf{r}}}$$ (40)

with i =1, 2 and

$$\mathit{\delta }{\mathbf{u}}_{\mathbf{g}}=\left[\mathbf{K}\right(\mathbf{u}\left){]}^{-1}\mathbf{g}\right(\mathbf{u},\lambda )or\mathit{\delta }{\mathbf{u}}_{\mathbf{g}}=[\mathbf{K}\left(\mathbf{u}\right){]}^{-1}\mathbf{g}(\mathbf{y},\lambda )$$ (41)

$$\mathit{\delta }{\mathbf{u}}_{\mathbf{r}}=\left[\mathbf{K}\right(\mathbf{u}){]}^{-1}{\mathbf{F}}_{\mathbf{r}}$$ (42)

In Equation (40) Δu is the displacements increment vector. Figure 2 shows the incremental-iterative procedure of the Potra-Pták Method combined with the Linear Arc-Length technique. The input data in the algorithm are: initial arc length (⁰Δl); maximum number of iterations in each load step (k_max); number of desired iterations in each load step (k_d); tolerance (tol); load increment (ΔP); and maximum number of load steps (LS_max). The outputs of the algorithm are: nodal displacement vector (u); total load parameter (λ); total number of load steps (LS); the total number of accumulated iterations until convergence to the solution (k_total); average number of iterations (k_av); and processing time in seconds (t).

Figure 2
Algorithm for the Potra-Pták method combined with Linear Arc-Length technique.

5. RESULTS AND DISCUSSION

The Euler-Bernoulli and Timoshenko beam-column models are used for the purpose of comparing them, by the equilibrium path of structural systems with geometric nonlinearity (considering large displacements and rotations, but small strains). Four structures are analyzed numerically: cantilever beam; Lee's frame; semi-circular arch; and Williams' frame. The units of the problems were maintained according to the original references.

The semi-rigid connection is modeled with the element proposed by DEL SAVIO et al. [²⁰] and it is inserted at the intersection point between beam or column and support. The constitutive relation of the material is considered linear elastic and the structures' own weight is neglected in the analyzes.

According to CHEN et al. [²⁸], axial, shear, bending and torsional forces are transmitted to a connection. However, by hypothesis in the analyzes performed here, large numerical values are adopted for axial and shear stiffness (that is, S_a ≅ ¥ and S_t ≅ ¥, respectively) and only the value of rotational stiffness is varied (S_r).

5.1 Cantilever beam

Consider the cantilever beam subjected to a load P illustrated in Figure 3, whose longitudinal elasticity modulus is E = 1.0 × 10⁷ kN/m², the shear modulus G is considered to be the half of E and coefficient m = 1.0. This problem was also studied by SILVA et al. [²⁹], using the beam-column element proposed by TANG et al. [³⁰]. Figures 4 and 5 show the equilibrium paths obtained with the Euler-Bernoulli and Timoshenko models depending on the slenderness ratio (L₀/h). The parameters considered in the analyzes for the Potra-Pták solution method are: initial arc length ⁰Δl = 0.3; tolerance tol = 1.0 × 10^-6; desirable number of iterations in the load step k_d = 5; maximum number of iterations in the step k_max = 150; and increment load DP = 1.0 × 10³ kN. The structure is discretized by 40 beam-column elements and one connect element. To simulate the fixed support, the rotational stiffness S_r = 1.0 × 10¹⁵ (rigid connection) is adopted.

Figure 3
Structural model of cantilever beam.

Figure 4
Equilibrium paths of the cantilever beam for slenderness ratio l = L₀/h = 2.

It can be seen, in Figure 4, that for the slenderness ratio L₀/h = 2 the Euler-Bernoulli model presented a more rigid numerical curve, differing from the numerical response of the Timoshenko model; however, the curves tend to be closer together at smaller load intensities. With the increase in slenderness (L₀/h = 10), it is observed that the responses of the models were very close (Figure 5). The numerical results obtained with the Timoshenko beam theory are shown in Table 1, as to the total number of load steps (LS) and accumulated iterations until convergence to the solution (k_total), average number of iterations per load step (k_av) and processing time in seconds (t).

Figure 5
Equilibrium paths of the cantilever beam for slenderness ratio l = L₀/h = 10.

Table 1
Cantilever beam - numerical results with Timoshenko's theory.

Slenderness ratio	LS	ktotal	kav	t (s)
L0/h = 2	16	30	1.8750	2.7075
L0/h = 10	15	28	1.8670	2.5161

5.2 Lee’s frame

The structure illustrated in Figure 6, known as Lee's frame, is often used to validate nonlinear solution strategies. This numerical problem was studied by SCHWEIZERHOF and WRIGGERS [³¹] and BATTINI [¹⁹]. The vertical and horizontal bars of the frame have cross-section area A = 6.0 and moment of inertia I = 2.0. The material of the bars has the properties: longitudinal elasticity modulus E = 720; shear modulus G = 276.9230; Poisson's ratio n = 0.3; and coefficient m = 1.0. This structure was discretized with 20 beam-column elements (ten elements for each bar) and two connect elements. The parameters adopted for the solution method are: ⁰Δl = 9.0; tol = 1.0 × 10^-7; k_d = 5; k_max = 150; and DP = 2.0.

Figure 6
Structural model of Lee’s frame.

Figure 7 shows the equilibrium paths (vertical displacement at node 14 x load P curves) with two force limit points and two displacement points, for both beam models. Figure 8 shows the structure equilibrium configurations referring to the initial position (LS = 0) and the deformed positions (LS = 5, 20, 40 and 56), obtained with the Timoshenko beam model. It is observed in Figure 7 the good numerical agreement between the curves obtained here and equilibrium points obtained by BATTINI [¹⁹]. The load-displacement curves considering the two beam theories were very close, since the members (beam and column) of the frame are slender (slenderness ratio l = L₀/h = 60). The numerical results obtained (LS, k_total, k_av and t) with Timoshenko's beam theory are shown in Table 2.

Figure 7
Lee’s frame - equilibrium paths.

Figure 8
Lee’s frame - equilibrium configurations.

Table 2
Lee’s frame: numerical results with Timoshenko's theory.

Solution Method	LS	ktotal	kav	t (s)
Potra-Pták	56	173	3.0892	7.3113

5.3 Semi-circular arch

Consider the semi-circular arch of radius R = 50 m subjected to a centered load P at its top and with pinned supports (S_r = 0), as shown in Figure 9. The mechanical properties of the material are: longitudinal elasticity modulus E = 2.0 × 10⁴ kPa; shear modulus G = 8.3333 × 10³ kPa; Poisson's ratio n = 0.2; and coefficient m = 1.0. The arch has cross-sectional area A = 0.8 m² and moment of inertia I = 4.2667 m⁴.

Figure 9
Structural model of semi-circular arch.

The structure was discretized by 36 beam-column elements, two connect elements and 39 nodes. The following parameters are considered in the analyzes with the Potra-Pták method: ⁰Δl = 2.2; tol = 1.0 × 10^-8; k_d = 7; k_max = 150; and DP = 100 kN. Table 3 shows the numerical results (LS, k_total, k_av and t) obtained with Timoshenko's beam theory. Figure 10 corresponds to the equilibrium paths in the load-displacement space of the structure with respect to the node where the load P is applied (node 20 of the mesh), obtained with the Euler-Bernoulli and Timoshenko beam models. Four force limit points (horizontal tangents) and three displacement limit points (vertical tangents) are observed in these curves. Figure 11 shows the structure equilibrium configurations referring to the initial position (LS = 0) and the deformed positions (LS = 50, 150, 250 and 301), obtained with the Timoshenko beam model.

Table 3
Semi-circular arch - numerical results with Timoshenko's theory.

Solution Method	LS	ktotal	kav	t (s)
Potra-Pták	301	607	2.0166	50.0421

Figure 10
Semi-circular arch - equilibrium paths.

Figure 11
Semi-circular arch - equilibrium configurations.

5.4 Williams' frame

This example aims to investigate the effect of connection flexibility of the Williams' frame. The structure illustrated in the Figure 12 consists of two rigidly connected beams, with the supports located at the ends and one vertical load P applied at its top. The axial and flexural stiffness of the beams are EA = 1.885 × 10⁶ lb (≅ 8.3848 × 10⁶ N) and EI = 9.274 × 10³ lb in² (≅ 26.6145 N m²), respectively. The sectional shear correction coefficient is m = 1.0. Three types of support are considered: ideally pinned (S_r = 0); perfectly rigid (S_r = 1.0 × 10¹⁵ lb in/rad ≅ 1.1298 × 10¹⁴ N m/rad); and semi-rigid (S_r = 1.8 × 10³ lb in/rad ≅ 203.3717 N m/rad). The structure was discretized by eight beam-column elements and two connect elements. The parameters for the solution method considered in the simulations are: ⁰Δl = 0.065; tol = 1.0 × 10^-8; k_d = 4; k_max = 150; and DP = 1.6 lb (≅ 7.1171 N).

The numerical results (LS, k_total, k_av and t) with the Timoshenko model and the types of supports are shown in Table 4. The equilibrium paths (vertical displacement/h versus load P curves) obtained with two load limit points are shown in Figure 13. There is good numerical agreement with the equilibrium points obtained by TIN-LOI and MISA [³²].

Figure 12
Structural model of the Williams' frame.

Table 4
Williams' frame: numerical results.

Types of supports	LS	ktotal	kav	t (s)
Ideally pinned	23	46	2.0000	0.9709
Perfectly rigid	21	27	1.2857	0.6714
Semi-rigid	22	44	2.0000	0.9212

Figure 13
Williams' frame - equilibrium paths.

5.5 Numerical results analysis

Lee's frame and semi-circular arch are characterized by a strongly nonlinear behavior, in which the equilibrium paths exhibit load limit points (snap-throughs) and displacement limit points (snap-backs). The tangent is horizontal at the load limit points (maximum or minimum points) and at the displacement limit points the tangent is vertical in the path. When the structure reaches the limit points it can become unstable, so their identification is of great importance for an engineering design.

The equilibrium paths with the Euler-Bernoulli and Timoshenko beam models are close to the initial loading levels, regardless of the structural element slenderness. However, as the load intensity increases the paths may diverge (Figure 4 - cantilever beam, Figure 10 - semi-circular arch and Figure 13 - Williams' frame). The structural system response does not depend on the theory adopted for higher slenderness values (relation between length L₀ and height h of the transversal section), that is, the load versus displacement curves become close or coincident (Figure 5 - cantilever beam and Figure 7 - Lee's frame). The Timoshenko beam model for thick structural elements (lower slenderness) produces more accurate results due to the inclusion of the shear effects.

The critical load and the post-buckling strength of structural systems are affected by the slenderness ratio and are lower for structures composed of elements with low slenderness. In other words, the large displacement behavior of these systems is affected by shear deformation, and this deformation must be taken into account when the structure components have low slenderness.

The iterative and incremental procedure of Potra-Pták associated with the Linear Arc-Length path-following technique was able to obtain equilibrium paths with load and displacement limit points of the studied structures. Note that in the iteration of this procedure, the same stiffness matrix (Jacobian matrix) is used to solve the systems of linear equations (see Equations (34) and (35)); thus, these systems can be solved via decomposition (for example, LU or Cholesky decomposition), since a single factorization at the beginning of the iteration is necessary.

The Potra-Pták method has already been tested in nonlinear problems of plane and spatial trusses [²², ³³], proving to be more efficient than the Newton- Raphson method, regarding processing time. Despite the more expensive computational cost of iteration, compared to the Newton-Raphson iteration, it achieves the problem solution with a lower number of load steps and accumulated iterations necessary for convergence, given the tolerance. The Potra-Pták method can be quite advantageous in solving nonlinear equations systems of structural problems with a large number of unknowns.

The device for checking the change in the sign of the initial load increment (Δl⁽⁰⁾) proposed by KRENK and HEDEDAL [³⁴], in the predicted solution (see lines 11 to 13 of the algorithm in Figure 2), is easy to implement and presented accurate results in all the studied examples. Furthermore, to guarantee the correct signal of the load sub-increment (Δl^(k)) in the iterative cycle (load increment correction phase), the following condition was added: if Δu^(k)T Δu^(k+1) < 0, then dl₂ ¬ -dl₂ (see lines 36 to 38 of the algorithm in Figure 2).

Models with ideal connections simplify the analysis procedure, but often do not represent the actual behavior of the structure. The nodal connections of plane frames are subject to the action of bending moments, axial forces and shear forces; however, the effect of these forces is usually neglected, and only the influence of the bending moments is of practical interest. From the numerical results of the Williams' frame, the flexibility of the connections has significant influence on the structural behavior, both in the nodal displacements and in the distribution of the internal forces. It is also noted that the structure failure load decreases with the increase in the connection flexibility.

6. CONCLUSION

This paper presented a numerical-computational model for static analysis of plane frames with geometric nonlinearity and semi-rigid connections. The structural systems were discretized with 2D beam-column co-rotational finite elements, based on the Euler-Bernoulli and Timoshenko beam theories. The connection between the structural members was simulated by a zero-length element, which considers axial, translational and rotational rigidities. The system of nonlinear equations that describes the structural problem was solved by the Potra-Pták incremental-iterative procedure associated with the Linear Arc-Length technique.

The numerical results showed that for geometric nonlinear analyzes with thick elements (less slenderness) the shear effect becomes important as the magnitude of the load on the structure increases. Thus, a finite element model whose formulation is based on Timoshenko's theory is more precise in terms of obtaining the structural response represented by the equilibrium path with limit points.

The simplifying hypothesis that the beam-column connection is perfectly flexible or completely rigid in structural designs can lead to an incorrect estimate of the mechanical behavior of frames. Connect deformation has a negative effect on the frame stability, since it increases the lateral displacement of the structure and causes a decrease in the stiffness of the member that is connected to the joint.

In order to perform the nonlinear analysis of structures, with greater precision, it is of utmost importance that incremental and iterative procedures that can appropriately consider the effects of large rotations and displacements, and allow the complete tracing of the structure's equilibrium path. The Potra-Pták method combined with the Linear Arc-Length path-following technique proved to be suitable for the analysis of structural systems (beams, frames and arches) with geometric nonlinear behavior.

Therefore, the consideration of geometric nonlinearity, the proper beam theory and the semi-rigidity of beam-column connections in structural analysis and design procedures can reduce the system's own weight, thus obtaining more optimized structures.

BIBLIOGRAFIA

• ¹
  MUÑOZ, P., FERNANDO, L., GONÇALVES, P.B., et al, “Nonlinear Resonance Analysis of Slender Portal Frames under Base Excitation”, Shock and Vibration, pp. 1-21, 2017.
• ²
  SANTOS, M.N., ROCHA, P.A.S., SILVA, A.R.D., et al, “Application of a nonlinear hybrid finite element in the structural steel modeling”, REM: Revista Escola de Minas, v. 65, n. 1, pp. 19-28, 2012.
• ³
  MAXIMIANO, D.P., SILVA, A.R.D.D., SILVEIRA, R.A.D.M., “Iterative strategies associated with the normal flow technique on the nonlinear analysis of structural arches”, REM: Revista Escola de Minas, v. 67, n. 2, pp. 143-150, 2014.
• ⁴
  ELKARANSHAWY, H.A., ELERIAN, A.A., HUSSIEN, W.I., “A Corotational Formulation Based on Hamilton’s Principle for Geometrically Nonlinear Thin and Thick Planar Beams and Frames”, Mathematical Problems in Engineering, pp. 1-22, 2018.
• ⁵
  CRISFIELD, M.A., Non-linear Finite Element Analysis of Solids and Structures, Vol 1, Chichester, England, John Wiley & Sons Ltd, 1991.
• ⁶
  YAW, L.L., 2D Corotational Beam Formulation, Walla Walla, USA, Walla Walla University, 2009.
• ⁷
  KIEN, N.D., “Effects of shear deformation on large deflection behavior of elastic frames”, Vietnam Journal of Mechanics, v. 26, n. 3, pp. 167-181, 2004.
• ⁸
  LIMA, A.S.D., FARIA, A.R.D., “A C1 Beam Element Based on Overhauser Interpolation”, Latin American Journal of Solids and Structures, v. 14, n.1, pp. 92-112, 2017.
• ⁹
  ONYIA, M.E., ROWLAND-LATO, E.O., “Determination of the critical buckling load of shear deformable unified beam”, International Journal of Engineering and Technology, v. 10, n. 3, pp. 647-657, 2018.
• ¹⁰
  JONES, S.W., KIRBY, P.A., NETHERCORT, D.A., “The analysis of frames with semi-rigid connections-a state-of-the-art report”, Journal of Constructional Steel Research, v. 3, n. 2, pp. 2-13, 1983.
• ¹¹
  DEGERTEKIN, S.O., HAYALIOGLU, M.S., “Design of non-linear semi-rigid steel frames with semi-rigid column bases”, Electronic journal of structural engineering, v. 4, pp. 1-16, 2004.
• ¹²
  SEKULOVIC, M., SALATIC, R., “Nonlinear analysis of frames with flexible connections”, Computers & Structures, v. 79, n. 11, pp. 1097-1107, 2001.
• ¹³
  VAN HAI, N., NGHIEM, D.N.T., CUONG, N.H., “Large displacement elastic static analysis of semi-rigid planar steel frames by corotational Euler-Bernoulli finite element”, Journal of Science and Technology in Civil Engineering (STCE)-NUCE, v. 13, n. 2, pp. 24-32, 2019.
• ¹⁴
  NGO-HUU, C., NGUYEN, P.C., KIM, S.E., “Second-order plastic-hinge analysis of space semi-rigid steel frames”, Thin-Walled Structures, v. 60, pp. 98-104, 2012.
• ¹⁵
  HADADE, M.A.S., CATOIA, B., FERREIRA, M.A., et al, “Experimental characterization of the momentrotation relationship of the beam-column connection in precast concrete”, Matéria (Rio J.), v. 23, n. 3, pp. e12164, 2018.
• ¹⁶
  MOKHTAR, R., IBRAHIM, Z., JUMAAT, M.Z., et al, “Behaviour of Semi-rigid Precast Beam-to-column Connection Determined Using Static and Reversible Load Tests”, Measurement, pp. 108007, 2020.
• ¹⁷
  WANG, F., YANG, J., PAN, Z., “Progressive collapse behaviour of steel framed substructures with various beam-column connections”, Engineering Failure Analysis, v. 109, pp. 104399, 2020.
• ¹⁸
  REIS, M.C.J., CODA, H.B., “Physical and geometrical non-linear analysis of plane frames considering elastoplastic semi-rigid connections by the positional FEM”, Latin American Journal of Solids and Structures, v. 11, n. 7, pp. 1163-1189, 2014.
• ¹⁹
  BATTINI, J.M., “Co-rotational beam elements in instability problems”, Doctoral thesis, KTH, Stockholm, 2002.
• ²⁰
  DEL SAVIO, A.A., ANDRADE, S.A.L., MARTHA, L.F., et al, “Um sistema não-linear para análise de pórticos planos semi-rígidos”, Revista Sul-Americana de Engenharia Estrutural (Passo Fundo), v. 2, n. 1, pp. 97-125, 2005.
• ²¹
  POTRA, F.A., PTÁK, V., Nondiscrete induction and iterative processes, Vol. 103, Boston, MA, Pitman Advanced Publishing Program, 1984.
• ²²
  SOUZA, L.A.F., CASTELANI, E.V., SHIRABAYASHI, W.V.I., et al, “Trusses Nonlinear Problems Solution with Numerical Methods of Cubic Convergence Order”, TEMA (São Carlos), v. 19, n. 1, pp. 161-179, 2018.
• ²³
  SCILAB, version 6.1.0. ESI Group, 2020.
• ²⁴
  SOLEYMANI, F., SHARMA, R., LI, X., et al, “An optimized derivative-free form of the Potra–Pták method”, Mathematical and Computer Modelling, v. 56, n. 5-6, pp. 97-104, 2012.
• ²⁵
  RIKS, E., “The application of Newton's methods to the problems elastic stability”, Journal of Applied Mechanics, v. 39, n. 4, pp. 1060-1065, 1972.
• ²⁶
  RIKS, E., “An incremental approach to the solution of snapping and buckling problems”, International Journal of Solids and Structures, v. 15, n. 7, pp. 529-551, 1979.
• ²⁷
  RAMM, E., “Strategies for tracing the nonlinear response near limit points”, In: W. Wunderlich et al editors, Nonlinear Finite Element Analysis in Structural Mechanics, pp. 63-89, Berlin, Heidelberg, Springer-Verlag, 1981.
• ²⁸
  CHEN, W.F., GOTO, Y., LIEW, J.R., Stability design of semi-rigid frames, New York, USA, John Wiley & Sons, 1996.
• ²⁹
  SILVA, J.L., LEMES, I., SILVEIRA, R.A., et al, “Influência da teoria de viga na análise geometricamente não linear de estruturas reticuladas”, In: Proceedings of the XXXVII Iberian Latin-American Congress on Computational Methods in Engineering (ed Suzana Moreira Ávila), Brasília, DF, Brazil, 6-9 Nov 2016.
• ³⁰
  TANG, Y.Q., ZHOU, Z.H., CHAN, S.L., “Nonlinear beam-column element under consistent deformation”, International journal of structural stability and dynamics, v. 15, n. 05, pp. 1450068, 2015.
• ³¹
  SCHWEIZERHOF, K.H., WRIGGERS, P., “Consistent linearization for path following methods in nonlinear FE analysis”, Computer Methods in Applied Mechanics and Engineering, v. 59, n. 3, pp. 261-279, 1986.
• ³²
  TIN-LOI, F., MISA, J.S., “Large displacement elastoplastic analysis of semirigid steel frames”, International Journal for Numerical Methods in Engineering, v. 39, n. 5, pp. 741-762, 1996.
• ³³
  SOUZA, L.A.F., CASTELANI, E.V., SHIRABAYASHI, W.V.I., et al, “Métodos iterativos de terceira e quarta ordem associados à técnica de comprimento de arco linear”, Ciência & Engenharia, v. 26, n. 1, pp. 39-49, 2017.
• ³⁴
  KRENK, S., HEDEDAL, O., “A dual orthogonality procedure for non-linear finite element equations”, Computer Methods in Applied Mechanics and Engineering, v. 123, n. 1-4, pp. 95-107, 1995.

Publication Dates

• Publication in this collection
  11 Oct 2021
• Date of issue
  2021

History

• Received
  03 June 2020
• Accepted
  29 Nov 2020