Second-order two-cycle analysis of frames based on interpolation functions from the solution of the beam-column differential equation

Abstract

In geometrically nonlinear problems solved using the Finite Element Method (FEM), the structure response is directly influenced by the level of discretization and the nonlinear solution algorithm used. To reduce the discretization dependence, exact solutions are developed based on the deformed infinitesimal element equilibrium. To deal with the nonlinear solution problem, the two-cycle method can be used, since it is not dependent on load or displacement steps. The two-cycle method developed by Chen & Lui (1991) uses the classical geometric matrix and is not accurate for high axial loads. This happens because the geometric matrix is obtained using Hermitian polynomials which are approximate solutions. To circumvent this issue, the frame element’s tangent matrix is obtained using interpolation functions that match the homogeneous solution of the differential equation of the beam-column problem. The main objective of this study is to carry out a second order analysis of the frames and obtain equilibrium paths using the two-cycle method and the tangent stiffness matrix based on solutions of the differential equations obtained from the element’s deformed configuration. The results in terms of displacements and rotations for the examples studied are identical to the analytical solutions, showing that the combination of the two-cycle method with the exact element formulation is promising and can diminish the need for discretization.

Keywords:
differential equations; second order effects; geometrically nonlinear analysis.

1. Introduction

The evolution of engineering demands the development of economic structures with less weight and material consumption without any trade-off in durability and safety. Hence, the use of high strength materials and slender elements becomes more common in new ventures. Under this view, the consideration of geometric nonlinearity is indispensable, especially for slender elements.

In context of the geometric nonlinearity of frames, a second-order stability analysis can be carried out in many ways, from simple checks based on first-order linear analyses to higher-order nonlinear numerical analyses. Despite the existence of several analytical solutions to solve stability problems, numerical solutions using the Finite Element Method (FEM) are often employed. From this numerical approach, the continuous structure response is directly influenced by the engineer’s experience in choosing the most suitable number of elements to be used in structural analysis (discretization). All these influences occur because the FEM discrete solution approximates the analytical solution, i.e., the interpolation functions that define the structure deformed configuration do not always agree with the problem’s exact solution (Burgos & Martha, 2013; Rodrigues et al., 2019).

In literature, to reduce element discretization dependence using the FEM in nonlinear geometric analysis, several studies have been carried out, such as Burgos et al., (2005), who used additional degrees of freedom within the elements in a classical linearization of the stability problem. Also, Chen & Lui (1991) and Aristizabal-Ochoa (1997, 2008, 2012), used classical or modified stability functions. Other authors have sought solutions based on infinitesimal element equilibrium in their deformed or undeformed configurations considering a combination of effects, such as transversal and axial loads, shear deformation (Timoshenko beam theory) and elastic foundation. Davis et al., (1972) and Nukulchai et al., (1981) formulated the exact solution for beam elements considering shear deformation. Zhaohua & Cook (1983) and Ting & Mockry (1984) developed the finite beam element considering two and one-parameter elastic foundations respectively, but with no shear deformation. Shrima & Giger (1992), in turn, formulated Timoshenko’s beam finite elements considering a two-parameter elastic foundation. Onu (2008) was the first author to formulate a finite element combining all effects (transversal and axial loads, shear deformation and two-parameter elastic foundation). Burgos & Martha (2013) presented exact shape functions including transversal and axial loading and shear deformation. Rodrigues et al. (2021) developed the three-dimensional tangent stiffness matrix and nodal equivalent loads, using exact interpolation functions, considering axial load, shear deformation and high-order terms of the Green-Lagrange strain tensor.

In addition to the exact finite elements previously addressed in literature, a nonlinear analysis methodology is also needed in numerical methods to describe the equilibrium path of a structural system. Thus, incremental, and iterative methods are commonly used, and some stand out among them: Newton-Raphson method, secant method and arc length method (Wempner, 1971; Ricks, 1972, 1979; Ramm, 1981; Crisfield, 1983, 1986), displacement control method (Yang & Shieh, 1990; Yang & Kuo, 1994) and indirect displacement control method (de Borst, 1986, 1987), among others. Along with these traditional general methods, there are simplified approaches capable of solving nonlinear problems with less computational effort, such as the γ_z coefficient method (NBR 6118, 2014), moment amplification method (Chen & Lui, 1991), P-Delta method (Chen & Lui, 1991) and, finally, the two-cycle method (Chen & Lui, 1991), which is applied in this study. The main objective of this research is to carry out a second order analysis of frames and obtain equilibrium paths using the two-cycle method and the tangent stiffness matrix based on exact solutions of the differential equations obtained for the element’s deformed configuration.

The article is organized in the following manner: the next section introduces the differential equations that govern the problem. Section 3 presents the exact tangent stiffness matrix. The two-cycle method is presented in section 4 and some examples of second order analysis of frames are presented in section 5. Finally, in section 6, the authors conclude this article and suggest ideas for future research.

2. Solution of the beam-column differential equation

Figure 1 (Rodrigues et al., 2021) shows the deformed configuration of an infinitesimal beam element, subjected to transversal (q) and axial loads (P). Considering all forces involved, equilibrium imposition leads to Equations (1) and (2):

Figure 1
Deformed infinitesimal beam element in equilibrium (Rodrigues et al., 2021).

(1) $$\sum F_y=0\to -dV+q(x)dx=0\to \frac{dV}{dx}=q(x)$$

(2) $$\begin{array}{c}\sum M_0=0\to dM-(V+dV)dx-Pdv+q(x)\frac{d{x}^2}{2}=0,M(x)=El\frac{d\theta (x)}{dx}\\ El\frac{d^2\theta (x)}{d{x}^2}-V(x)-P\frac{dv(x)}{dx}=0\to El\frac{d^3\theta (x)}{d{x}^3}-\frac{dv(x)}{dx}-P\frac{d^{2}v(x)}{d{x}^2}=0\end{array}$$

Herein, the Navier (Euler-Bernoulli) theory is employed. Thus, the rotation and transversal load can be written as the derivative of the lateral displacement and shear force respectively, leading, after some algebraic manipulation, to Equation (3):

(3) $$\frac{d^{4}v(x)}{d{x}^4}-\frac{P}{EI}\frac{d^{2}v(x)}{d{x}^2}=\frac{q(x)}{EI}$$

The previous equation can be rewritten considering only its homogeneous part (q(x)=0), which leads to v_h (x):

(4) $$\frac{d^4{v}_h(x)}{d{x}^4}-{\mu }^2\frac{d^2{v}_h(x)}{d{x}^2}=0,\mu =\sqrt{\frac{P}{EI}}$$

The differential equation solution results in lateral displacement v_h (x) and the rotation can be obtained using θ(x) = dv_h ⁄ dx:

(5) $$v_h(x)=c_1{e}^{\mu x}+c_2{e}^{-\mu x}+c_{3}x+c_4$$

(6) $${\theta }_h(x)=\mu c_1{e}^{\mu x}-\mu c_2{e}^{-\mu x})+c_3$$

According to Burgos & Martha (2013), P can be either a tensile axial load (P>0), in which case μ is a real number and the expression for v_h (x) can be written in terms of hyperbolic functions. On the other hand, in the case of a compressive axial load (P<0), μ is a complex number and the expression for v_h (x) can be written in terms of trigonometric functions. All those expressions are available in Rodrigues et al. (2019). In matrix form, the conditions given by Equations (5) and (6) are written as:

(7) $$\left\{\begin{array}{c}v(x)\\ \theta (x)\end{array}\right\}=[X]\cdot \{C\},[X]=\left[\begin{array}{cccc}{e}^{\mu x}& e^{-\mu x}& x& 1\\ \mu e^{\mu x}& -\mu e^{-\mu x}& 1& 0\end{array}\right][C]=\left\{\begin{array}{l}{C}_1\\ C_2\\ C_3\\ C_4\end{array}\right\}$$

With boundary conditions given by:

(8) $$\begin{array}{c}[H].\{C\}=\{d\},\\ \{d\}=\left\{\begin{array}{l}{v}_1\\ {\theta }_1\\ v_2\\ {\theta }_2\end{array}\right\}=\left\{\begin{array}{l}v(0)\\ \theta (0)\\ v(L)\\ \theta (L)\end{array}\right\}\{H\}=\left[\begin{array}{cccc}1& 1& 0& 1\\ \mu & -\mu & 1& 0\\ e^{4\mu }& e^{-44}& L& 1\\ \mu e^{4\mu }& -\mu e^{-4\mu }& 1& 0\end{array}\right]\end{array}$$

The same process can be performed using the expressions in terms of hyperbolic or trigonometric functions and the stiffness matrices are then obtained accordingly.

3. Tangent stiffness matrix based on interpolation functions from the beam-column differential equation

The shape functions are used to interpolate the nodal values of the discrete element, obtaining displacements and rotations within the element of the continuous solution of the problem, according to Equation (9):

(9) $$\left\{\begin{array}{l}v(x)\\ \theta (x)\end{array}\right\}=[N].\{d\},[N]=\left[\begin{array}{cccc}{N}_2^v& N_3^v& N_5^v& N_6^v\\ N_2^{\theta }& N_3^{\theta }& N_5^{\theta }& N_6^{\theta }\end{array}\right]\{d\}=\left\{\begin{array}{l}{v}_1\\ {\theta }_1\\ v_2\\ {\theta }_2\end{array}\right\}$$

Inverting matrix [H] in Equation (8) and substituting in Equation (7), shape functions can be easily written:

(10) $$\left\{\begin{array}{l}v(x)\\ \theta (x)\end{array}\right\}=[X]\cdot [H{]}^{-1}.\{d\}\to [N]=[X]\cdot [H{]}^{-1}$$

Burgos & Martha (2013) presented the expressions and plots of these shape functions for all cases described above, as well as the development to obtain the tangent stiffness matrices described in Equations (11) to (13), based on the Principle of Minimum Potential Energy. In case the homogeneous solution is given by exponential functions, the coefficients of the tangent stiffness matrix are given by:

(11) $$\begin{array}{c}{K}_{22}=K_{55}=-K_{25}=K_{52}=\frac{12EI}{L^3}\frac{(\mu L{)}^3(e^{\mu L}+1)}{12D}\\ K_{23}=K_{26}=K_{32}=-K_{35}=-K_{53}=-K_{56}=K_{62}=-K_{65}=\frac{6EI}{L^3}\frac{(\mu L{)}^2(e^{\mu L}-1)}{6D}\\ K_{33}=K_{66}=\frac{4EI}{L}\frac{\mu L[\mu L(e^{2\mu L}+1)-(e^{2\mu L}-1)]}{4D(e^{\mu L}-1)}\\ K_{36}=K_{63}=\frac{2EI}{L}\frac{\mu L[(e^{2\mu L}+1)-2\mu L{e}^{\mu L}]}{2D(e^{\mu L}-1)}\\ D=\mu L(e^{\mu L}+1)-2(e^{\mu L}-1),\mu =\sqrt{\frac{P}{EI}}\end{array}$$

In the case of μ being a real number (P>0), coefficients in terms of hyperbolic functions are:

(12) $$\begin{array}{c}{K}_{22}=K_{55}=-K_{25}=K_{52}=\frac{12EI}{L^3}\frac{(\mu L{)}^3\sinh (\mu L)}{12D}\\ K_{23}=K_{26}=K_{32}=-K_{35}=-K_{53}=-K_{56}=K_{62}=-K_{65}=\frac{6EI}{L^3}\frac{(\mu L{)}^2(\cosh (\mu L)-1)}{6D}\\ K_{33}=K_{66}=\frac{4EI}{L}\frac{\mu L[\mu L\cosh (\mu L)-\sinh (\mu L)]}{4D}\\ K_{36}=K_{63}=\frac{2EI}{L}\frac{\mu L[\sinh (\mu L)-\mu L]}{2D}\\ D=2-2\cosh (\mu L)+\mu L\sinh (\mu L),\mu =\sqrt{\frac{P}{EI}}\end{array}$$

In the case of μ being a complex number (P<0), coefficients in terms of trigonometric functions are:

(13) $$\begin{array}{c}{K}_{22}=K_{55}=-K_{25}=K_{52}=\frac{12El}{L^3}\frac{(\mu L{)}^3\sin (\mu L)}{12D}\\ K_{23}=K_{26}=K_{32}=-K_{35}=-K_{53}=-K_{56}=K_{62}=-K_{65}=\frac{6EI}{L^3}\frac{(\mu L{)}^2(\cos (\mu L)-1)}{6D}\\ K_{33}=K_{66}=\frac{4EI}{L}\frac{\mu L[\mu L\cos (\mu L)-\sin (\mu L)]}{4D}\\ K_{36}=K_{63}=\frac{2EI}{L}\frac{\mu L[\sin (\mu L)-\mu L]}{2D}\\ D=2-2\cos (\mu L)+\mu L\sinh (\mu L),\mu =\sqrt{\frac{-P}{EI}}\end{array}$$

The axial components of the tangent stiffness matrix of the frame element are the same as those of the classical elastic stiffness matrix (EA/L), easily found in literature, in which E is the material Young’s modulus, L is the element length and A is the cross-section area.

4. Two cycle method

Developed by Chen & Lui (1991), the two-cycle method uses only two iterations to obtain second-order effects. In this method, the classical Equation (14) is utilized in two steps. In the first step, an elastic linear analysis is carried out using the classical elastic stiffness matrix - easily found in literature.

(14) $$\{f\}=[K]\{d\}$$

From the axial load components obtained in the previous step, the second step starts with the substitution of these loads in the tangent stiffness matrix coefficients. A second solution is obtained using Equation (14) again, but replacing the elastic stiffness matrix with the tangent stiffness matrix. Thus, in this second step, displacements and forces will be of second order, i.e., considering geometric nonlinearity. Herein, the tangent stiffness matrix coefficients are given by Equations (11) to (13). In the original two-cycle method, these coefficients are obtained by the sum of the elastic and geometric matrices and are not as accurate. One of the problems that arise from this approach is the possibility of numerical singularities (division by zero) that come from low values in axial loads. This can be easily solved by adding a tolerance that flags the need for a nonlinear analysis.

It is interesting to point out that in obtaining the differential equation, it is assumed that the axial load is known. Thus, in this method, the first step is necessary only to calculate and increment the axial loads for the second step. In short: the major contribution to achieve reliable results closer to the analytical solution is in the exact approach of the tangent stiffness matrix in the second cycle. Since the differential equation is obtained assuming that the axial load is constant within the element, the discretization using one element per bar is not only advantageous, but also mandatory. Discretization using more than one element per bar is not comparable with the analytical solution, for the same reasons exposed before.

5. Numerical examples

In this study, MATLAB® (2018) was used to obtain the solutions for the models shown in Figure 2. The equilibrium paths for second-order analysis using the two-cycle method with exact shape functions were compared to the analytical solutions and results obtained using classical frame elements (Hermitian shape functions) through MASTAN2 software (McGuire et al., 2000), which adopts the predictor-corrector method for the nonlinear analysis. For the two-cycle solution, each load value results in a displacement vector which is independent from previous steps, since the method is used to obtain the exact displacement solution, and not in an iterative manner.

Figure 2
Proposed models for second-order analysis. a) Non-sway frame. b) Sway frame. c) Cantilever beam.

5.1 Example 1: non-sway frame

Figure 2a shows a laterally restrained (non-sway) frame subjected to a vertical load P and bending moments M = αPL = 0.001PL. Properties of all bars are: L = 6 m, E = 10⁸ kN/m² , I = 10^-5 m⁴, and EA = 10⁶ kN. The critical load is P_cr= 2.5515 π² EI ⁄L² = 699.51 kN and the analytical solution in terms of the free rotation is given by:

(15) $$\theta (L)=\frac{\alpha \mu L^2[2\mu \cos (\mu L)-2\mu +{\mu }^{2}L\sin (\mu L)]}{4\cos (\mu L)+{\mu }^2{L}^2\cos (\mu L)+\mu L\sin (\mu L)-4}$$

A look at Figure 3, shows that the results for the two-cycle method using the exact tangent stiffness matrix is equivalent to the equilibrium path of the analytical solution and better than the MASTAN2 solution discretized in 1, 2 and 4 elements (M1, M2 and M4). In this example, it is clear how the engineer’s experience would be required in discretization, especially considering that the difference between the results produced using discretization M1 and M2 is 44.10%. One of the main objectives of this study is to reduce the need for discretization in simple problems, such as nonlinearities in planar frames.

Figure 3
Equilibrium paths for non-sway frame.

5.2 Example 2: sway frame

Figure 2b shows an unbraced (sway) frame subjected to a vertical load P and horizontal loads H = αP=0.001P. Properties are the same as the previous example. The critical load is P_cr = 0.74766 π² EI ⁄ L² = 204.98 kN and the analytical solution in terms of displacement at the top is given by:

(16) $$v(L)=\frac{\alpha [12\cos (\mu L)+L^2{\mu }^2\cos (\mu L)+5\mu L\sin (\mu L)-12]}{\mu [6\sin (\mu L)+\mu L\cos (\mu L)]}$$

Figure 4 shows that the results for two-cycle method using the exact tangent stiffness matrix are identical to the analytical solution. In other words, the combination of the two-cycle method with exact element formulation provides excellent results in predicting the nonlinear geometric behavior of the sway frame.

Figure 4
Equilibrium paths for sway frame.

5.3 Example 3: cantilever beam

Figure 2c shows a cantilever beam subjected to a horizontal (compressive) load P and vertical load αP=0.01P. Properties are the same as the previous examples. The critical load is P_cr = π² EI ⁄ 2L² =68.54 kN and the analytical solution as shown in McGuire et al. (2000) is given by:

(17) $$v(L)=\alpha L[\frac{\tan (\mu L)}{\mu L}-1]$$

Figure 5 shows that the results for M1 and M4 discretization, despite being close to the analytical solution, are not identical as the ones obtained by the two-cycle method with exact element formulation.

Figure 5
Equilibrium paths for cantilever beam.

6. Conclusion

This article presented the combination of a simple method for second-order analysis with a formulation of the element tangent stiffness that uses shape functions that match the homogeneous solution of the differential equation in deformed configuration to solve nonlinear geometric problems.

From the examples shown, the analyses using the two-cycle method with the suggested formulation of the tangent stiffness proved to be efficient when compared to traditional numerical solutions; furthermore, results are also independent of model discretization. It is interesting to notice that result accuracy in numerical solutions by finite elements strongly depends on the boundary conditions of the models. The sway frame and cantilever beam achieved better results than the non-sway frame in one element discretization (M1) for the numerical solution in MASTAN2. This issue alone shows the need for a tool that does not rely on user experience in discretizing the model and providing nonlinear parameters (load step, solution algorithm). This problem is no longer present when the combination of the two-cycle method with exact tangent stiffness is applied. It is important to emphasize that all the analyses are performed for small displacements and loads that are below critical values, i.e., postcritical behavior is not considered.

In future studies, the authors suggest the development of analogous analyses that consider shear deformation, elastic foundation and/or variable cross-sections. Another suggestion is the use of distributed loads.

References

• ARISTIZÁBAL-OCHOA, J. D. Firstand second-order stiffness matrices and load vector of beam columns with semirigid connections. Journal of Structural Engineering, v. 123, n. 5, p. 669-678, 1997.
• ARISTIZÁBAL-OCHOA, J. D. Matrix method for stability and second-order analysis of timoshenko beam-column structures with semi-rigid connections. Engineering Structures, v. 34, p. 289-302, 2012.
• ARISTIZÁBAL-OCHOA, J. D. Slope-deflection equations for stability and second-order analysis of timoshenko beam-column structures with semi-rigid connections. Engineering Structures, v. 30, p. 2517-2527, 2008.
• ASSOCIAÇÃO BRASILEIRA DE NORMAS TÉCNICAS. NBR 6118: projeto de estruturas de concreto - Procedimento. Rio de Janeiro: ABNT, 2014.
• DE BORST, R. Non-linear analysis of frictional materials. 1986. Dissertation. Delft University of Technology, Delft, 1986.
• DE BORST, R. Computation of post-bifurcation and post-failure behavior of strain softening solids. Comput. Struct. v. 25, p. 211-224, 1987.
• BURGOS, R. B.; MARTHA, L. F. Exact shape functions and tangent stiffness matrix for the buckling of beam-columns considering shear deformation. In: IBERIAN LATIN AMERICAN CONGRESS ON COMPUTATIONAL METHODS IN ENGINEERING, 34, 2013, Pirenópolis. Anais […], 2013.
• BURGOS, R. B.; SILVA, R. R.; MARTHA, L. F. Avaliação de cargas críticas e comportamento pós-crítico inicial de pórticos planos. In: IBERIAN LATIN AMERICAN CONGRESS ON COMPUTATIONAL METHODS IN ENGINEERING, 35, 2005, Guarapari. Anais […], 2005.
• CHEN, W. F.; LUI, E. M. Stability design of steel frames. Boca Raton: CRC Press, 1991.
• CRISFIELD, M. A. An arc-length method including line searches and accelerations. International Journal for Numerical Methods in Engineering, v. 19, p. 1269-1289, 1983.
• CRISFIELD, M. A. Snap-through and snap-back response in concrete structures and the dangers of under-integration. International Journal for Numerical Methods in Engineering, v. 22, n. 3, p.751-767, 1986.
• DAVIS, R.; HENSHELL, R. D.; WARBURTON, G. B. A timoshenko beam element. Journal of Sound and Vibration. v. 22, n. 4, p. 475-487, 1972.
• MATLAB. 9.7.0.1190202 (R2019b). Natick, Massachusetts: The MathWorks Inc., 2018.
• MCGUIRE, W.; GALLAGHER, R. H.; ZIEMIAN, R. D. Matrix Structural Analysis. EUA: John Wiley & Sons, 2000.
• NUKULCHAI, W. K.; DAYAWANSA, P. H.; KARASUDHI, P. An exact finite element model for deep beams. International Journal of Structures, v. 1, n. 1, p. 1-7, 1981.
• ONU, G. Finite elements on generalized elastic foundation in timoshenko beam theory. Journal of Engineering Mechanics, v. 134, n. 9, p. 763-776, 2008.
• RAMM, E. Strategies for tracing the Non-Linear response near Limit-Points. In: WUNDERLICH, W. (ed.) Nonlinear Finite Element Analysis in Structural Mechanics. Berlin: Springer-Verlag, 1981. p. 63-89.
• RIKS, E. The application of Newton’s Methods to the problems of elastic stability. Journal of Applied Mechanics, v. 39, p. 1060-1066, 1972.
• RIKS, E. An incremental approach to the solution of snapping and buckling problems. International Journal of Solids and Structures, v. 15, p. 529-551, 1979.
• RODRIGUES, M. A. C.; BURGOS, R. B.; MARTHA, L. F. A Unified approach to the timoshenko 3D Beam-Column element tangent stiffness matrix considering higher-order terms in the strain tensor and large rotations. International Journal of Solids and Structures, v. 222-223, 111003, 2021. 22 p.
• RODRIGUES, M. A. C.; BURGOS, R. B.; MARTHA, L. F. A unified approach to the timoshenko geometric stiffness matrix considering higher-order terms in the strain tensor. Latin American Journal of Solids and Structures, v. 16, n. 4, e185, 2019.
• SHRIMA, L. M.; GIGER, M. W. Timoshenko beam element resting on two-parameter elastic foundation. Journal of Engineering Mechanics, v. 118, n.2, p. 280-295, 1992.
• TING, B. Y.; MOCKRY, E. F. Beam on elastic foundation finite element. Journal of Structural Engineering, v. 110, n. 10, p. 2324-2339, 1984.
• WEMPNER, G. A. Discrete approximations related to nonlinear theories of solids. International Journal of Solids and Structures, v. 7, p. 1581-1599, 1971.
• YANG, Y. B.; KUO, S. B. Theory and analysis of nonlinear framed structures. New York: Prentice Hall, 1994.
• YANG, Y. B.; SHIEH, M. S. Solution method for nonlinear problems with multiple critical points. American Institute of Aeronautics and Astronautics, v. 28, n. 12, p. 2110-2216, 1990.
• ZHAOHUA, F.; COOK, R. D. Beam elements on two-parameter elastic foundations. Journal of Engineering Mechanics, v. 109, n. 6, p. 1390-1402, 1983.

Publication Dates

• Publication in this collection
  14 Apr 2023
• Date of issue
  Apr-Jun 2023

History

• Received
  23 July 2022
• Accepted
  17 Nov 2022