An Alternative Approach to Obtaining Stiffness and Mass Matrices for Three-Dimensional Bernoulli-Euler and Timoshenko Beam-Columns

In the static analysis of beam-column systems using matrix methods, polynomials are using as the shape functions. The transverse deflections along the beam axis, including the axialflexural effects in the beam-column element, are not adequately described by polynomials. As an alternative method, the element stiffness matrix is modeling using stability parameters. The shape functions which are obtaining using the stability parameters are more compatible with the system’s behavior. A mass matrix used in the dynamic analysis is evaluated using the same shape functions as those used for derivations of the stiffness coefficients and is called a consistent mass matrix. In this study, the stiffness and consistent mass matrices for prismatic three-dimensional Bernoulli-Euler and Timoshenko beam-columns are proposed with consideration for the axial-flexural interactions and shear deformations associated with transverse deflections along the beam axis. The second-order effects, critical buckling loads, and eigenvalues are determined. According to the author’s knowledge, this study is the first report of the derivations of consistent mass matrices of Bernoulli-Euler and Timoshenko beam-columns under the effect of axially compressive or tensile force.


Introduction
In linear structures, deformations are proportional to the external loads, the displacements and the internal forces of the system are obtaining via superposition. When an element of a structure is subjecting to an axial-flexural effect, non-linear interactions develop. This interaction produces second-order moments, which consist of the product of the axial force and flexural displacements. The moment magnifies nonlinearly as the load increases, which results in increasingly nonlinear load-displacement relationships. Consequently, the effect of any axialflexural interaction must consider in the analysis of those structures. The second-order effects are cause changes in the geometry of the structural system, which may result in instability issues. Stability is a structural condition defined as the loading of a structure or structural component in which a slight disturbance in the loads or geometry does not produce large displacements. Several building codes require the considerations of the effect of loads on the structure's overall stability and its members. [1][2] stipulate using a design method that includes the axialflexural effects on the overall stability of the structure or structural components. Researchers have explored various approaches for assessing structure's stability in the analysis and design of framed structures and have presented many books and articles on this topic. Favorable information for determining the second-order effects and analysis methods can be found from the following studies [3][4][5][6][7][8][9][10]: Computer-aided designs are using in the analysis of structural systems the well-known method of this topic is the finite element method (FEM). For frame-elements are developed special technics of FEM is called matrix methods. The accuracy of matrix methods in computeraided design for beam-column systems is dependent on the compatibility of the assumed shape functions with the real elastic behavior of the element. The transverse deflection along the beam axis, which is included in the axialflexural effect of the beam element, is no longer a polynomial function, but a function with trigonometric variables.
The shape functions which are obtaining using the stability parameters are more compatible with the system's behavior. Mass matrices are one of the main components in a dynamic analysis. For a dynamic analysis to be performed using matrix methods, a mass matrix containing the appropriate features must be established. These features are derived using lumped or distributed mass techniques. A distributed mass matrix evaluation with the same shape functions used for derivation of the stiffness matrix is called a consistent mass matrix. Other researchers have developed various approaches for performing static and dynamic analyses on prismatic and nonprismatic Bernoulli-Euler and Timoshenko beams. Some of the scholars that have evaluated mass and stiffness matrices include [11] where examined the free vibration and stability analysis of Timoshenko beams through a finite element analysis. [12] presented a symbolic integration that combined an indefinite integral with the Gauss divergence theorem, which was used to form a novel integration scheme for a consistent mass matrix in the FEM. [13] developed a new procedure for the FEM and applied it to the computation of the vibrational response of a Timoshenko beam subject to a uniformly distributed harmonic load. [14] discussed improvements to the mass and geometric stiffness matrices of a Bernoulli-Euler two-dimensional beam. [15] derived a stiffness matrix and consistent mass matrix using the principles of d'Alembert and virtual work. [16] presented a method for evaluating the shape functions and structural matrices derived for non-prismatic Euler-Bernoulli beam elements. [17] analyzed non-prismatic Timoshenko beams using a basic displacement function. [18] used a force-based approach to derive a three-dimensional consistent mass matrix for vibration analysis of beams. [19] utilized a size-dependent model based on the modified strain gradient theory.
The remainder of this paper is organized as follows: in sections 2 and 3 are obtained the stiffness matrices in prismatic three dimensional Bernoulli-Euler and Timoshenko beam elements, respectively. The matrices are modeled using stability parameters, = � / . In sections 4 and 5 are presented consistent mass matrices for both element types.

The stiffness matrices in prismatic three dimensional Bernoulli-Euler beam element
The following assumptions are made in this study: 1) behavior is purely elastic, 2) all members are prismatic, homogenous, and isotropic (i.e., the elasticity module, E, moment of inertia, I, and cross-sectional area, A, are constants), 3) in the Bernoulli-Euler elements shear deformations are neglected, plane sections normal to the neutral axis remain plane and normal to the neutral axis, 4) in the Timoshenko beam elements shear is constant along the elements, 5) no distributed load along the beam, 6) the elements are slender. The end displacements of a prismatic three dimensional Bernoulli-Euler beam element are shown in Fig. 1. A constant axial compressive force subjected to this element along the x axis.
The differential equation governing the displacement, y, of this prismatic element, AB, subjected to a constant compressive axial force, P, is as follows [3]: The displacements in the x-y plane numbered 2, 6, 8, and 12 and the displacements in the x-z plane numbered 3, 5, 9, and 11 can be obtained from the solution of Eq. (1). That is showing in the following Eq.(3).
The end rotations are obtaining from the derivation of Eq. (3), which becomes When we want to calculate the coefficients, C, of Eq. (3) due to unit rotation, in the end, A, is used following end conditions: For (u6= 1); at x=0 Y(x) = 0 and ′ ( ) = 1; at x=L, Y(x) = 0 and ′ ( ) = 0 The obtained coefficients, C, are: After combining Eqs. 4 and 5, end moments is calculated using the Equation − ′′ ( ) = ( ).
where coefficients a and b are Eqs.7 are the basic parameters to calculate of the element end forces and stiffness matrices. The stiffness matrix, [ ], is obtained using the following submatrices: The submatrices are: where y and z denote the principal axes. When calculating parameters λy, λz and ay, by and az, bz are used the Eqs.
(2) and (7a)-(7b) for these principal axes. G denotes the torsional elastic modulus and can be assumed to be equal to the shear modulus, while J denoted the torsional constant.

The stiffness matrices in prismatic three dimensional Timoshenko beam element
Stiffness matrices of a Timoshenko beam should be obtained using shear deformations in the element. For this purpose, consider a simply supported beam is shown in Fig. 2.

Fig.2. End rotations of a simply supported beam subjected to a unit moment in end A
Second and third derivations of Eq. (3) should be used to obtain the moments and shear forces in the beam. These are: End rotation of a simply supported beam obtain from first derivative of Eq. (3): ′ (0) = 1 + 3 . When the above obtained integral coefficients combined to this derivative: When unit moment is applied to the end B of that beam due to the symmetry should be = = and according to Maxwell-Betti rule is equal for both ends.
When end rotations due to shear force effect combined to the moment effect thus could be obtained total end rotations. For this purpose is used the virtual work principle. Total end rotations are as: When the beam end moments are defined as MAB and MBA. The end rotations for the positive end moments are as: In the cases of ≠ 0 and = 0 related end moments can be found using Eqs. (15). These are: where, and coefficients in Eqs. (17) could be obtained from Eqs. (14). It should be noted that an axial compressive force creates second-order effects in the span of a beam-column element. These effects are representing by P-δ. It is also necessary to obtain the flexural moment and displacement values induced by the P-δ effects. These values should be calculating with consideration for the nonlinear behaviors of the element. The calculation steps for obtaining the P-δ effects are given in Appendix 2.

Bernoulli-Euler beam element mass matrix
A distributed mass matrix evaluated with the same shape functions used for the derivation of the stiffness matrix is called a consistent mass matrix. For an element subjected to a constant axially compressive force, the procedures for obtaining a consistent mass matrix are explaining here.
A mass matrix for a linearly elastic system is expressed in matrix form as follows: The mass matrix coefficients for transverse deflections numbered 2, 6, 8, 12 and 3, 5, 9, 11 these are the coefficients related to the bending mass forces about the z and y axes respectively, could be calculated using the following Equation: Eq. (19) will be applied to the defined shape functions from Eq. (3). The coefficients for indices 4 and 10 in the mass matrix belong to the torsional mass forces. The following equations (20a) and (20b) defined the shape functions complied with for these cases required boundary conditions. where J is torsional constant. Similarly, the mass matrix coefficients for deflections numbered 1, 7 are the coefficients related to the axial mass forces in the beam, could be calculated using the following Equation: The same shape functions determined by the Eq. (20a)-(20b) can be used and following Eqs. (24a)-(24b) are obtained In Eq (19) determined mass matrix coefficients will be obtained from the integration of the shape functions. The functions to be integrated are Obtaining a closed form solution for Eqs. (26)-(35) is extremely complicated, and these equations should be evaluated using numerical methods.
The stability parameter for indices 2, 6, 8, and 12 these coefficients should be calculated using the flexibility rigidity, EIz. The related stability parameter is The stability parameter for indices 3, 5, 9, and 11 these coefficients should be calculated using the flexibility rigidity, EIy. The related stability parameter is = � (33b)

Timoshenko beam element mass matrix
In a Timoshenko beam, the effects of shear force deformation should be added to the functions ( ) ( ) in Equation (19). The following assumptions were made for this case: no distributed load along the beam, hence in a Timoshenko beam element shear force is constant along the element.
The general displacement function which included the shear deformation term as 1 .
( ) = 1 sin + 2 cos First, second and third derivations of Equation (34): According to the directions of the element coordinate system shear force due to in the positive direction applied end moments using Eq. (35c) is: Deformation due to this constant shear in the beam: When the Equations (36) and (37) combine: where, After the obtaining coefficients the terms, Cij and D1j of mass matrix could be calculated using Equation (19) as similar explained in section 4.

Numerical example
In this example, the analysis of a 3-dimensional frame system shown in Fig. 3 is performed using the method outlined in this study. The elements of the system are to be considered as Bernoulli-Euler beam elements. Comparisons of the obtained results from with and without axial-flexural effects are made. The joints numbered 1 and 2 are fixed-end while the joint numbered 3 is free to displacements. The frame is under the effect of three concentrated loads at joint 3 which are in ̅ , �, and ̅ directions -120 kN, 300 kN, and -10 kN, respectively.
The analysis includes the following steps: 1) obtaining and comparisons of the results considering the axial-flexural effects, 2) calculation of eigenvalues, 3) obtaining of the critical buckling load factors, LF. The beams comprise 20x50 mm steel profiles and bent about the major principal axis, y-y. This major principal axis y-y is parallel to the ̅ -� plane of the global coordinate system. Element properties: A=1000 mm 2 , Iy=208333 mm 4 , Iz=33333 mm 4 , J=99600 mm 4 , E=200000 Mpa, G=77200 Mpa, unit volume weight, w=7.85x10 -5 N/mm 3 , so a unit volume mass density of ρ=8 × 10 −9 N-s 2 /mm 4 u2, u3, u5, u6, u8, u9, u11 and, u12 (see Fig. 1   Mass matrix for displacements, u2, u6, u8, and, u12 (bending about z-z axis) (see Fig. 1) The eigenvalues are summarized in Table 3. The -120 kN, 300 kN and, -10 kN initial forces in the system, are increased by small increments and the determinant of the stiffness matrix of the system, [K] is calculated. The ratio of loads at the load level that makes the determinant zero to the first load determines the critical buckling load factor, LF, of the system.

Conclusions
Analyzing a frame system using the matrix method provides favorable results than the finite element method because, while a mesh arrangement is required in the finite element method, matrix methods enable each element taken into account in the frame system to be considered independently. To perform a successful static and dynamic analysis, identifying the stiffness and mass matrices of the elements is essential. That ensures an accurate representation of the system's behavior. In cases where the axial-flexural interactions are insignificant, assuming polynomial shape functions does not affect the outcome. However, in cases where the effects of the axial forces are significant, designers are faced with problems that included second-order effects and buckling issues. In the commonly used methods the elastic stiffness matrices are assumed to be identical to those used in unaffected cases. However, to account for the second-order effects, geometric stiffness matrices are added to or subtracted from the elastic stiffness matrices, as they are related to the sign of the axial force. In this study, the suggested shape functions are basing on the actual behavior of the beam. The transverse deflection along the beam axis considered the second-order effect, cannot be adequately described by a polynomial; instead, a function with trigonometric variables is necessary. A consistent mass matrix is the main component of dynamic analysis. According to the definition of consistent mass matrix must be established using the same shape functions for the derivations of the stiffness matrices. In the proposed procedure, the mass matrix is obtaining using the same shape functions as those used for the derivations of the stiffness matrices. As could be seen from the presented example the results obtained from the static and dynamic analysis were remarkably different when the axial forces played a significant role. When the level of axial forces was deemed insignificant, the results of both methods were in good agreement with each other. The parameters, a and b, are the main factors that affected the results. The limit values for the functions in Eqs. (29)-(35) for λ=0 were undefined, and their limits approached the values determined using the commonly use method. The author recommends that designers who will be preferred design using the Bernoulli-Euler or Timoshenko beam-columns use the proposed procedures for all axial force levels except λ=0.
The given expressions for the end forces of a beam element subjected to a compressive axial force can be used in the case of tension; however, in Eqs. (7a) and (7b) the defined a and b parameters should be calculated as follows [21]: where the α and coefficients are

Appendix 2. Calculation of the P-δ effect
A beam-column subjected to end moments M1 and M2 is displayed in Fig. A2.1. Among the end moments, the one with a greater absolute value is denoted as M2. In a column subjected to a constant axial force, P, the related transverse deflection of the column is given by Eq. (3) as follows: ( ) = 1 sin + 2 cos + 3 + 4 (A2.1) The integration constants C1, C2, C3, and C4 are determined from the boundary conditions, which are as follows [6]: The moment at the section of the column located at distance x from the end of the column is ( ) = 1 � − � + 2 + | 2 | > | 1 | It should be noted that the signs of the end moments were determined according to the classical moment rules and | 2 | > | 1 | when estimating Eqs. (A2.1)-(A2.5). Therefore, if an element has a single curvature, M1 and M2 should be positive as M1/M2> 0. However, for the case of reverse curvatures, the ratio of M1 to M2 should satisfy 1/ 2 < 0. Equations (A2.4) and (A2.5) can be solved numerically by scanning across the column length using small intervals.