Free-Vibration and Buckling Analyses of Beams using Kriging-Based Timoshenko Beam Elements with the Discrete Shear Gap Technique

.


Introduction
The finite element method (FEM) has now become a standard numerical method for solving boundary-value problems appearing in various engineering mathematical models.In the standard FEM, the approximated trial function over an element is constructed using a polynomial interpolation from the function values at the element nodal points.Kriging-based finite element method (K-FEM) is an alternative numerical method to the standard FEM.In the K-FEM, the approximated trial function over an element is constructed using a Kriging interpolation from the function values not only at the element nodes but also the values at a set of nodes around the element (called satellite nodes) [1][2][3].The element and its surrounding elements covering all the nodes constitute a domain of influencing nodes (DOI).One of the key advantages of K-FEM over the standard FEM is its ability to generate a stress field with higher accuracy and better smoothness using the same mesh and degrees of freedom.Additionally, K-FEM allows for solution refinements without changing the mesh or adding more degrees of freedom.A drawback of the K-FEM is that its computational procedure is much more expensive than that of the standard FEM.
A major difficulty in the development and application of the K-FEM to shear-deformable beam, plate bending, and shell models is, as in the standard FEM, the shear locking phenomenon [4][5][6].This phenomenon is a condition where the finite element model becomes excessively stiff as the beam, plate, or shell increasingly thins due to the presence of spurious transverse shearing stresses.In the applications of the K-FEM to the shear deformable (Timoshenko) beam model, the selective reduced integration (SRI) technique has been utilized to eliminate the shear locking [6].While the SRI has been proven to be effective in eliminating the locking, the resulting shear force is only accurate at the single integration sampling point.Moreover, the use of the SRI worsens the overall shearing force distribution.
Another technique to overcome the shear locking is proposed by Bletzinger et al. [7], namely the discrete shear gap (DSG) technique.The basic idea of the DSG technique is to replace the approximate displacement-based shear strain field with a substitute shear strain field determined from the interpolation of 'shear gaps' at the element nodes using the standard element shape functions.The DSG technique has been applied to eliminate shear locking in the Krigingbased Timoshenko beam elements [8] (referred to as K-beam-DSG0 elements).In this application, the substitute shear strain field was derived from a Kriging interpolation of shear gaps at all nodes in the DOI.The developed Kbeam-DSG0 elements, however, are only shear-locking free for the elements with a cubic basis function and three element-layer DOI.After that, the implementation of the DSG was modified by changing the substitute shear strain field to the standard linear interpolation of shear gaps at the Kriging-based element nodes [9].With this new implementation, the Kriging-based beam elements of all types (referred to as K-beam-DSG1 elements) are free from locking and able to give remarkably accurate displacements and bending moments for very thin to thick beams.The resulting shear force distribution, however, is piecewise constant over the elements.
The previous study of the K-beam-DSG1 elements, [9], is limited to the case of linear static analysis of beams.This paper aims to further develop and test the K-beam-DSG1 elements for free vibration and bifurcation buckling analyses of thin-to-thick beams.The beams include prismatic as well as axially functionally graded tapered beams (that is, linearly varying cross-section beams with varying constituent material properties along the beam axis [10]).Hamilton's principle is employed as the basis in deriving the finite element equations.Consistent Kriging shape functions are utilized to develop the mass matrix and the axial stress (geometric) stiffness matrix.The developed Kriging-based beam elements are subsequently tested for several free vibration and buckling benchmark problems.

Variational Form of the Governing Equations
Consider a Timoshenko beam model of the length L, made of an isotropic material but in general may not be homogeneous in the axial direction, with the Cartesian coordinates x-z and the deflection of the neutral axis, w, and the cross-sectional rotation, θ, as shown in Figure 1.The beam is subjected to an axial force P = P(t), a function of time t, positive in compression, at the ends of the beam.The axial force P is imposed at the outset, for example, by a pre-stressed bar or cable.The motion of the beam is described by the deflection and the rotation as functions of coordinate 0 ≤ x ≤ L and time t ≥ 0, that is, w = w(x, t) and θ = θ(x, t).The governing equations for free vibration analysis of the beam including axial force effects on the beam bending can be established using Hamilton's principle [11][12], viz.
In this equation, T is the kinetic energy of the beam, U is the internal strain energy in the beam, VP is the potential energy associated with the initial axial force P, while δ is the variational operator.The principle states that the integral over any time interval t1 to t2 of the variation of the difference in the kinetic and total potential energies (U + VP) with respect to any virtual defection and rotation equals zero.The kinetic energy is given as [6,12].
x z w θ where A is the area of the cross section, Iy is the second moment of area of the cross section about the y-axis, and  is the mass density of the beam material (mass per unit volume).The dot sign over a variable denotes the partial derivative with respect to the time variable t.The strain energy in the beam is given as [6]: where E and G are the modulus of elasticity and shearing modulus of the beam material, respectively, As = kA is the effective shear area, with k being a shear correction factor [13] to account for the difference of a constant shear stress distribution based on the Timoshenko beam theory and a parabolic actual shear stress distribution along the beam height.The comma symbol after a field variable denotes the partial derivative with respect to the independent variable after the comma, that is, coordinate x.The geometric parameters A, As, and Iy, and the material parameters , E, and G may generally vary along the beam length.Lastly, the potential energy of the axial force P is given as [12], [14].

Kriging-Based Finite Element Formulation
Suppose the beam is divided into Ne finite elements with Np nodes.The kinetic energy, strain energy, and potential energy of the whole beam can be written as the sum of their element contributions, that is, where L e is the length of element number e.
Consider an element e and its domain of influencing nodes (DOI) covering n nodal points, as illustrated in Fig. 2 for a two-layer element DOI.The deflection field w and rotation field θ over the element can be approximated using Kriging interpolation as follows: (, ) ≈  ℎ (, ) =    ()  () (6a) In these equations, Na(x), a = 1, 2, …, n, are Kriging shape functions constructed based on the n nodal points in the DOI (see [6], [8], [9]).Vector d e (t) is a vector of the deflections and rotations at the nodal points in the DOI of element e.The numbering system for the indices is the local element-DOI numbering system.The number of nodal points n for each element varies depending on the number of element layers to cover the DOI and depending on whether the element is exterior or interior.For example, for the two-element layer DOI as shown in Figure 2, n = 4 for the interior element while n = 3 for the left or right-end exterior element.
The Kriging shape functions Na(x), a = 1, 2, …, n, are determined using a system of Kriging equations.The details of Kriging equations have been presented in previous works [6,8] and will not be repeated in the paper.To construct the Kriging equations, one must decide on a polynomial basis function and a correlation function.In this work, linear to cubic basis functions and the quartic spline correlation function with the middle values of the correlation parameter(see [8,Table 2]) are employed to construct the Kriging equations.Introducing Equations ( 6) and ( 7) into Equation ( 5), the kinetic energy, strain energy, and potential energy of the axial force P can be written as, respectively, ̇T()    ̇() =1 (8a) T ()    () =1 (9a) Matrices m e , k e , and kg e defined by Equation (8b), Equation (9b), and Equation (10b) are the Kriging-based consistent mass matrix, stiffness matrix, and geometric stiffness matrix of element number e, respectively.The order of these matrices is 2n  2n.
Substituting the energy terms in Equation ( 1) by the expressions given by Equation (8a), Equation (9a), and Equation (10a) and carrying out the variational operation give Consider the first term of Equation (11).Interchanging the integral and summation operator and then integrating the first term by parts leads to Since the basic condition imposed on the variation is δd e (t1) = δd e (t2) = 0, the first term on the right-hand side of Equation ( 12) is equal to zero.Interchanging the summation and integral again and substituting Equation (12) into Equation (11) gives Equation ( 13) is written referring to the local element-DOI numbering system for each element e.This equation can also be expressed with reference to the global structural numbering system as follows: In these equations, D is the structural nodal displacement vector, that is, Vector δD is a vector of variation of D, which may also be interpreted as a vector of virtual nodal displacement.Matrices M, K, and Kg are the consistent structural mass matrix, the structural stiffness matrix, and the structural geometric stiffness matrix, respectively.Notation A in Equation ( 15) denotes the assembling algorithm in the framework of the K-FEM, which includes all nodes in the element DOI (not only the element nodes as in the assembly procedure of standard FEM).Since the variation δD is arbitrary over any time interval t1 to t2, Equation (14a) can be satisfied in general only when [11]  ̈() + ( −  g )() =  (16) This system of equations is the discretized governing equations for free vibration of the beam including axial force stiffening effects.Discretized equations for free vibration without axial stiffening and for bifurcation buckling analysis can be obtained from Equation ( 16) by simply reducing it, respectively, to In the case of free vibration, assuming the free vibration motion is simple harmonic with a circular frequency ω, the displacement vector can be expressed as [11] () =  0 sin( + ) where D0 represents the shape of the beam and φ is a phase angle.Introducing Equation (19) into Equation ( 17) leads to Therefore, the matrix equations finding free vibration frequencies, Equation (20), and a buckling critical load, Eq. ( 18), are in the category of generalized eigenvalue problems.

Implementation of the Discrete Shear Gap Technique
Direct application of the discretized beam models, Eqs. ( 17) and (18), results in an over-stiff beam model, in particular for slender beams due to the shear locking phenomenon.Continuing the previous work [9], the DSG technique [7] is applied to the discretized beam models in order to avoid shear locking.The details of the DSG formulation are presented in Ref. [9].
The DSG technique is implemented by replacing matrix (Nw,x − Nθ) in the shearing part of the stiffness matrix (i.e., the second term of the right-hand side of Eq. (9b)) with a substitute shear strain-displacement matrix, Bγ-bar, which is defined as [9]   =  1  2 (21a)

Results and Disscussion
This section presents the results of free vibration and buckling structural analyses of Timoshenko beams using different types of K-beam-DSG1 elements.Symbols of the form P*-* are used to identify different types of Krigingbased elements.The first asterisk after the letter P in the symbol is a number to denote the degree of the polynomial basis function used in Kriging interpolation, that is, '1' for a linear basis, '2' for a quadratic basis, and '3' for a cubic basis, while the second asterisk denotes the number of element layers of the domain of influencing nodes.For example, P2-3 means that a quadratic basis and three layers of elements are used in the construction of Kriging shape functions.In all analyses, the quartic spline correlation functions are employed to construct Kriging shape functions.
The correlation function parameters are taken to be the mid-value between the lower and upper bounds values presented in Wong et al. (see Table 2 in Ref. [8]).
The Gauss-Legendre quadrature is employed to evaluate the integrals corresponding to bending stiffness (the first term in the right-hand side of Equation (9b)), shearing stiffness (the second term in the right-hand side of Equation (9b)), DSG (Equation (21c)), mass (Equation (8b)), and geometric stiffness (Equation (10b)).For analysis of prismatic beams, the number of quadrature sampling points used is three to evaluate all of the integrals.This number is chosen based on some trials to obtain an accurate yet inexpensive number of quadrature points.For analysis of beams with varying geometrical and material properties, the number of quadrature points to evaluate all of the integrals is five, except for evaluating the DSG, the number of quadrature points remains three.

Natural Frequencies of Prismatic Beams
Consider the free vibration of thick and very thin simply supported beams.The geometrical and material parameters are L = 10 m, the width b = 1 m, k = 5/6, E = 200  10 9 N/m 2 ,  = 0.3,  = 5700 kg/m 3 , L/h = 5 for the case of thick beam and L/h = 1000 for the case of very thin beam.Following Lee and Schultz [15], the resulting circular natural frequencies, ωm, m = 1, 2, 3, …, are expressed in terms of nondimensionalized frequency parameters, λm, defined as The frequency parameters obtained using a converged solution of the pseudo-spectral method [15] are taken as the reference solution for the thick beam, whereas the frequency parameters based on the Euler-Bernoulli beam theory [15] are taken as the reference solution for the very thin beam.The simply supported beams were analyzed using different types of 16 K-beam-DSG1 elements of equal length.The resulting normalized frequency parameters for the cases of thick and very thin beams are presented in Tables 1 and  2, respectively.The results obtained using 16 frame elements contained in commercial software SAP2000 are also included for comparison.Table 1 demonstrates that, for the case of thick beam, all types of the K-beam-DSG1 elements yield highly accurate results.The errors for the first four natural frequencies are less than 1%.The most accurate ones are generally achieved using the type of P2-2, with a maximum error of 1.33% for the 11 th natural frequency.It is apparent that for the thick beam, the results of the K-beam-DSG elements are overall more accurate than those obtained using SAP2000.
Table 2 shows that, for the case of the very thin beam, the K-beam-DSG1 elements give very accurate results for several first low frequencies but are very inaccurate in predicting high frequencies.The errors of less than 1% are achieved up to the 5 th lowest frequencies.The most accurate ones are achieved using the type of P3-3, the errors of which are less than 1% up to the 8 th frequency.In general, for the very thin beam, the results of K-beam-DSG1 elements are less accurate than those obtained using SAP2000.
To study the convergence of the solutions using the P3-3 type of the K-beam-DSG1, the very thin beam was analyzed using a different number of elements, that is, 4, 8, 16, and 32.The resulting normalized frequency parameters are listed in Table 3.It is apparent that the frequency parameters converge to the Euler-Bernoulli beam frequencies.The very thin beam needs 32 K-beam-DSG1 P3-3 elements to give very accurate results for the first 15 th natural frequencies.

Critical Buckling Loads of Prismatic Beams
To examine the accuracy of the K-beam-DGS1 elements of different types in predicting the critical buckling load, an investigation was conducted on a fixed-fixed supported beam with varying length-to-thickness ratios (L/h = 5, 10, 20, 100, 1000).The beam was discretized using eight K-beam-DGS1 elements.The geometrical and material parameters are the same as in the previous free vibration analysis.For comparison, the beam was also analyzed using the commercial software SAP2000.The resulting critical loads were normalized to the corresponding exact solutions, which were obtained from the exact formula given as ( [16] as cited in [12]) where Leff is the effective buckling length (Leff = L/2 for fixed-fixed supported beams).The normalized critical loads and their corresponding exact values are listed in Table 4.
Table 4 demonstrates the accuracy of the K-beam-DSG1 elements for a wide range of bean thicknesses.The largest error occurred in the case of the thick beam (L/h = 5) with a maximum error of 2.7% for the element type P3-3.It can be seen that increasing the polynomial degree in the Kriging interpolation does not guarantee an increase in accuracy.The most accurate element type in general is P2-3, except for the thick beam, the most accurate type is P2-2.In comparison to the SAP2000 results, for the cases of thick and moderately thick beams (L/h = 5 and 10), the best performer of the K-beam-DSG1 yields more accurate results than SAP2000.For the cases of thin to very thin beams (L/h =20, 100, 1000), however, SAP2000 performs better than the best performer of the K-beam-DSG1.Furthermore, to study the convergence characteristics of the K-beam-DSG1 elements, the fixed-fixed supported beam with the length-to-thickness ratio L/h = 10 was discretized using different numbers of elements, viz.4, 8, 16, and 32.Table 5 shows the normalized critical load results.It can be seen that all types of the K-beam-DSG1 elements converge to the exact solution, or, at least to a value close to the exact solution.The convergence characteristic is, however, non-monotonic (fluctuating).where the subscript 0 denotes the value of the parameter at x = 0. Following Shahba et al. [10], the beam was analyzed using 30 K-beam-DSG1 elements of the types P2-2 and P3-3.
In addition, it was analyzed using 30 elements of the standard four-node cubic beam element with the DSG technique [17].Table 6 lists the free vibration analysis results together with the ones reported by Shahba et al. [10].The results of the four-node cubic element can be used as the reference solution since the total number of nodes involved in the analysis is 91, which is nearly 3 times larger than the other elements.Moreover, a study on free vibration analysis of a homogeneous tapered beam with taper ratio c = 0.2 (Table 7) shows that this standard element is the most accurate one compared to the exact results reported by Leung et al. ( [18] as cited in [10]).Table 6 shows the K-beam-DSG elements yield very accurate results with a maximum difference from the standard cubic element results of −0.29% (lower).As a comparison, the results of Shahba et al. [10] have a maximum difference of +0.32% (higher).

Critical Buckling Loads of Axially Functionally Graded Tapered Beams
Using the same beam as in the previous test but with different support cases, a series of buckling analyses using 30 K-beam-DSG1 elements of types P2-2 and P3-3 and using the standard four-node cubic beam elements was carried out.The resulting critical loads were expressed in terms of dimensionless critical load defined as The analysis results are presented in Table 8 together with those from Shahba et al. [10] for comparison.Taking the results obtained using the four-node cubic element as the reference (for the same reason as aforementioned), it is seen that the K-beam-DSG1 elements are able to give very accurate results for the cantilever and simply supported axially functionally graded beams.The results are generally more accurate than those of Shahba et al. [10].However, the results for the fixed-fixed supported beam are surprisingly very inaccurate with the maximum error of about 30% for the taper ratio c = 0.2 using the P3-3.The authors suspect that the inability of the K-beam-DSG1 elements to correctly predict the critical load for this particular case is due to the inability of the Kriging interpolation to capture the buckling mode.As a remark, the same fixed-fixed supported beam but using a homogeneous material has also been analyzed using the P3-3 and P2-2 K-beam-DSG1 elements.In this case, the results are very accurate with a maximum error of 0.32%.

Conclusions
A family of Kriging-based Timoshenko beam elements utilizing the DSG technique has been developed for free vibration and buckling analyses of prismatic as well as non-prismatic, axially functionally graded beams.The discretized equations were derived based on the Hamilton's principle.The performance of the K-beam-DSG1 elements of different types in predicting natural frequencies of free vibration and in predicting critical buckling loads has been examined using several prismatic and axially functionally graded taper beam examples.The examination results showed that the K-beam-DSG1 elements in general can predict natural frequencies and buckling loads very accurately with a reasonable number of elements.However, there was an anomaly in the case of a fixed-fixed axially functionally graded taper beam, in which the K-beam-DGS1 elements failed to predict the critical loads accurately.This anomaly needs further study to investigate the reasons for the failure and to improve the K-beam-DSG1 elements.

Figure 1 .
Figure 1.Coordinate System and the Positive Sign Convention for the Beam Deflection, W, Cross-Section Rotation, Θ, and Axial Load P

Figure 2 .
Figure 2. Beam Problem Domain, A Beam Element under Consideration, and Its Domain of Influencing Nodes This and the next numerical tests refer to the work of Shahba et al. [10] on free vibration and stability analyses of axially functionally graded beams.This test focuses on a tapered cantilever beam as shown in Figure 3, with three different taper ratios c = 0.2, 0.5, and 0.8.The tapered beam is made of two constituents, namely, aluminum (Ea = 70 GPa, a = 2702 kg/m 3 ) and zirconia (Ez = 200 GPa, z = 5700 kg/m 3 ).The material properties vary along the beam axis, as specified by the following:  = 200 − 130 ( s ratio is taken to be constant  = 0.3.The resulting natural circular frequencies are presented in dimensionless natural frequencies (m = 1, 2, 3, …)  ̅  =   2 =    2

Table 2 .
Normalized Frequency Parameters of A Very Thin Simply-Supported Timoshenko Beam (L/h = 1000)

Table 3 .
Convergence of Normalized Frequency Parameters of A Very Thin Simply-

Table 4 .
Normalized Critical Buckling Loads (Pcr / Pcr exact) of a Fixed-Fixed Supported Timoshenko Beam of Different Length-To-Thickness Ratios Modeled using Eight K-Beam-DSG1 Elements and Eight Framed Elements Contained in SAP2000

Table 5 .
Normalized Critical Buckling Loads (Pcr / Pcr exact) of A Clammed-Supported Timoshenko Beam of L/h = 10 Modeled using Different Numbers and Different Types of the K-beam-DSG1 Element

Table 6 .
The First Four Dimensionless Natural Frequencies of Axially Functionally Grade Tapered Cantilever Beams of Different Taper Ratios, Modeled using 30 Beam Elements of Different Types

Table 7 .
The First Four Dimensionless Natural Frequencies of A Homogeneous Tapered Cantilever Beam of Taper Ratio c = 0.2, Modeled using 30 Beam Elements of Different Types

Table 8 .
Dimensionless Critical Buckling Loads of An Axially Functionally Graded Timoshenko Beam of Different Taper Ratios, with Different Supports, Modeled using 30 Beam Elements of Different Types