Kriging-based Timoshenko Beam Element for Static and Free Vibration Analyses

An enhancement of the finite element method using Kriging interpolation (K-FEM) has been recently proposed and applied to solve oneand twodimensional linear elasticity problems. The key advantage of this innovative method is that the polynomial refinement can be performed without adding nodes or changing the element connectivity. This paper presents the development of the K-FEM for static and free vibration analyses of Timoshenko beams. The transverse displacement and the rotation of the beam are independently approximated using Kriging interpolation. For each element, the interpolation function is constructed from a set of nodes within a prescribed domain of influence comprising the element and its several layers of neighbouring elements. In an attempt to eliminate the shear locking, the selective-reduced integration technique is utilized. The developed beam element is tested to several static and free vibration problems. The results demonstrate the excellent performance of the developed element.


Introduction
In an attempt to improve the element-free Galerkin method with Kriging interpolation [1], Plengkhom and Kanok-Nukulchai [2] presented a new class of FEM by introducing Kriging shape functions in the conventional FEM.In this method, Kriging interpolation (KI) is constructed for each element using a set of nodes in a domain of influencing nodes (DOI) composed of several layers of elements (the DOI is in the form of polygon for 2D problems).Combining the KI of all elements, the global field variable is thus approximated by piecewise KI.For evaluating the integration in the Galerkin weak form, the elements are employed as integration cells.The method subsequently referred to as Krigingbased FEM (K-FEM) [3].
The advantages of the K-FEM are: (1).Highlyaccurate field variables and their gradients can be obtained even using the simplest form of elements (triangles in the 2D domain and tetrahedrons in the 3D domain).(2).The polynomial refinement can be achieved without any change to the element or mesh structure.(3).Unlike the moving Kriging elementfree Galerkin method [1], the formulation and coding of the K-FEM are very similar to the conventional FEM so that any existing general-purpose FE program can be easily extended to incorporate the enhanced method.Thus, the K-FEM has a high chance to be accepted in real engineering practices.
In the pioneering work [2], the K-FEM was developed for static analyses of 1D bar and 2D planestress/plane-strain solids.Subsequently, it was developed for analyses of Reissner-Mindlin plates [3,4] and improved through the use of adaptive correlation parameters and by introducing the quartic spline correlation function.A drawback of the K-FEM is that the interpolation function is discontinuous at the inter-element boundaries (except in 1D problems).In spite of this discontinuity, using appropriate choice of shape function parameters, the K-FEM passes weak patch tests and therefore the convergence is guaranteed [5].The basic concepts and advances of the K-FEM have been presented in several papers [6][7][8].The current development is the extension and application of the K-FEM to different problems in engineering, such as general plate and shell structures [9,10] and multi-scale mechanics [11].
Despite many attractive features of the K-FEM, in the application to shear-deformable plates and shells, the drawback of transverse shear locking and membrane locking presents in the K-FEM [3,4,9,10].The use of high order basis (cubic and quartic) in KI can alleviate the shear and membrane lockings, but there is no guarantee to eliminate the lockings completely.Until the writing of this paper, to the authors' knowledge, an effective method to completely eliminate the lockings in the K-FEM is not yet invented.
In attempt to invent the method to eliminate the drawback of shear locking in shear-deformable plate and shell problems, it is instructive to study the K-FEM in the simpler context of the Timoshenko beam since this problem can be considered as a 1D degeneration from the Reissner-Mindlin plate.It is the aim of this paper to present the development and testing of the K-FEM for analyses of static and free vibration of Timoshenko beams.The developed element is tested to several beam problems and its performance is compared with high performance Timoshenko beam element developed by Friedman and Kosmatka [12].

Kriging Interpolation in the K-FEM
Named after Danie G. Krige, a South African mining engineer, Kriging is a well-known geostatistical technique for spatial data interpolation in geology and mining [13,14].Using this interpolation, every unknown value at a point can be interpolated from known values at scattered points in its specified neighborhood.Here the concepts of the KI in the context of K-FEM are briefly reviewed.
Consider a continuous field variable u(x) defined in a domain Ω.The domain is represented by a set of properly scattered nodes xi, i=1, 2, …, N, where N is the total number of nodes in the whole domain.
Given N field values u(x1), …, u(xN), the problem of interest is to obtain an estimated value of u at a point x 0 ∈ Ω.

The Kriging estimated value
The Kriging weights are determined by requiring the estimator U h (x0) unbiased, i.e.
where E[•] is the expected value operator, and by minimizing the variance of estimation error, var[U h (x0)-U(x0)].Using the method of Lagrange [13,14] for the constraint optimization problem, the requirements of minimum variance and unbiased estimator lead to the following Kriging equation system (see Wong [10] for the complete derivation): Rλ + Pµ = r(x0) in which  The expression for the estimated value u h , Eq. (1), can be rewritten in matrix form where d = [u(x1) … u(xn)] T is an n x l vector of nodal values.Since the point x0 is an arbitrary point in the DOI, the symbol x0 can be replaced by symbol x.Thus, using the usual finite element symbol, Eq. ( 5) can be expressed as x N x d x (6) in which N(x)= λ T (x) is the matrix of shape functions.
In order to construct Kriging shape functions in Eq. ( 6), a polynomial basis function and a correlation function should be chosen.Basis functions ranging from polynomial of the degree one up to four have been utilized in the past works on the K-FEM [2][3][4][5][6][7][8][9][10][11].
In the problems of shear deformable plates and shells, it is necessary to use cubic or quartic polynomial basis in order to alleviate the shear and membrane lockings [4,9,10].
The correlation function ρ (h) is defined as: where h is a vector separating two points x and x+h and σ 2 is the variance of the random function U(x).
In the K-FEM, factor σ 2 has no effect on the final results and it was taken equal to 1 in this study.
Following the previous works [3][4][5][6][7][8][9][10], Gaussian and Quartic Spline correlation function (QS) were chosen.Gaussian correlation function is defined as and QS is defined as In these equations, θ>0 is the correlation parameter, h = ||h||, i.e. the Euclidean distance between points x and x+h, and d is a scale factor to normalize the distance.Factor d was taken to be the largest distance between any pair of nodes in the DOI.
For each aforementioned correlation function, appropriate values for the correlation parameter should be chosen so that the K-FEM will give reasonable results.Plengkhom dan Kanok-Nukulchai [2] proposed a rule of thumb for choosing the parameter as follows: Correlation parameter, θ, should be selected so that it satisfies the lower bound condition, and also satisfies the upper bound condition, where a is the order of basis function and b is the dimension of problem (1, 2, or 3).

Variational Form of Timoshenko Beam Governing Equations
The basic assumption of Timoshenko beam theory is that a plane normal to the beam axis in the undeformed state remains plane in the deformed state but it does not necessarily in normal direction to the neutral axis [16].The theory accounts for both the transverse shearing strain and the rotary inertia in a dynamic analysis.This section presents the variational formulation of Timoshenko beam following that given by Friedman and Kosmatka [12].
Fig. 1 shows the coordinate system used in the following formulation.The displacement components in x and z directions can be respectively written as a function of coordinate x and time t as where u is the axial displacement of a material point at coordinate (x, z), w is the transverse displacement (deflection) of the neutral axis and ψ is the rotation of the cross-section.Using the small-strain and displacement equations for general solids, Eqs.(11a) and (11b) give the nonzero strain components: where εxx is the normal strain in x direction and γxz is the transverse shearing strain.The commas denote the first partial derivatives with respect to the variable next to it (i.e.x).
The variational equation of motion of the beam can be derived using Hamilton's principle [12], viz.
( ) where δU, δT, and δWe, are the variations of the strain energy, the kinetic energy, and the work of external forces, respectively.The strain energy is given as in which L and A are the length and the crosssectional area of the beam respectively, and is the vector of normal and shearing stresses, and is the vector of normal and shearing strains.In Eq. (15a), E and G are Young's and the shear moduli of the beam material respectively, and k is a shear correction factor that is dependent upon the crosssection geometry.Substituting Eqs.(15a) and (15b) into Eq.( 14), considering the strains, Eq. (12a) and Eq.(12b), and integrating over the cross-section, Eq. ( 14) yields ( ) ( ) where I is the moment of inertia of the cross-section.
The kinetic energy of the beam is given as ( ) where ρ is mass density (per unit volume) of the beam material.The dots signify the first partial derivatives with respect to the time variable t.Substituting Eqs.(11a) and (11b) and integrating over the cross-section, Eq. ( 17) can be expressed as Finally, the work of external forces is given as.
in which q and m are the distributed forces and moments along the length of the beam.
Substituting Eq. ( 16), Eq. ( 18) and Eq. ( 19) into Eq.( 13), applying the variational operations, and applying Hamilton's principles lead to the variational equation (or the weak form) for the Timoshenko beam as follows, ( ) ( ) The double dots signify the second derivatives with respect to time t.
The bending moment and shear force along the beam can be calculated from the deflection w and rotation ψ as follows: ( ) ( ) ( ) is the element consistent mass matrix,d && is the element nodal acceleration vector, is the element stiffness matrix, and is the element equivalent nodal force vector.Note that the size of the matrices m and k, d, and f depends on the number of nodes in the DOI, n.For static problems, Eq. (24a) simply reduces to kd = f (25) The integrations in Eqs.(24b), (24c), and (24d) are evaluated using Gauss quadrature.It is well known in the conventional FEM that if the stiffness matrix, Eq. (24c), is evaluated using full integration, then the element becomes too stiff for thin beams (shear locking phenomenon, see e.g.Hughes et al. [15] and Reddy [16]).One of the techniques to overcome the shear locking is the selective-reduced integration (SRI), in which the stiffness matrix is evaluated using full integration for the bending term (the first term in the right-hand side of Eq. (24c)) and using reduced integration for the shearing term (the second term in the right-hand side of Eq. (24c)).The effectiveness of the SRI technique to eliminate the shear locking in the Kriging-based Timoshenko beam is investigated in this study through a series of numerical tests.
Using the finite element assembly procedure, one can obtain the global discretized equation for the beam vibration as follows: ( ) ( ) ( ) where M and K are the global mass and stiffness matrices, respectively, D is the global nodal displacement vector, D & & is the global nodal acceleration vector, and F is the global nodal force vector.The equations for static and undamped freevibration problems can be respectively written as The free vibration of a structure with the natural frequency ω can be written as where D0 is the amplitude of the vibration and φ is the phase angle.Entering Eq. (29) into Eq.( 28) leads to the eigen equation (K -ω 2 M)D0 = 0 (30) where D0 is a vector of eigenvalues, which is the mode shape of vibration, corresponding to an eigenvector ω 2 .

Numerical Results
The performance of the developed Kriging-based Timoshenko beam element was tested using a number of different static and free-vibration problems.In the following, for the sake of brevity, only the results for the beam element with the K-FEM options of the cubic basis, three-layer DOI (Fig. 2), and Gaussian correlation function (P3-3-G) are presented (see Syamsoeyadi [17] for the comprehensive results).The issues considered include the shear locking phenomenon, the accuracy and the convergence of the displacements, bending moments, shear forces, and natural frequencies.
Firstly, the appropriate range for the correlation parameter θ and the efficient number of Gaussian sampling points for evaluating the integrations in Eqs.(24b), (24c) and (24d) were determined by conducting a series of numerical tests.The tests for finding the lower bound and the upper bound of θ satisfying Eqs.(10a) and (10b) revealed that for the beam element with Kriging option P3-3-G, the lower and upper bound values for θ are 10 -4 and 1.9, respectively.In the subsequent analyses the parameter θ was taken to be equal to 1, which was nearly the mid-value between the lower and upper bounds.
To obtain an accurate yet efficient number of the integration sampling points, a cantilever beam subjected to triangular-distributed load as shown in Fig. 3 was analyzed using different number of sampling points for evaluating the stiffness matrices, Nsamp.Subsequently, it was analyzed with different number of sampling points for evaluating the force vector, Nbody,.The free-end deflection of the beam was observed and compared to the exact solution [12], viz.
Where qo, L, E, I, φ, v, and h are respectively the value of the triangular load at the clamped end, the length, the modulus of elasticity, the moment of inertia, the ratio of bending stiffness to shearing stiffness, Poisson's ratio and the thickness of the beam.
The results were presented in Tables 1 and 2. The tables show that at least two sampling points are needed to yield accurate integrations.In the subsequent analyses Nsamp was taken to be equal to 3 while Nbody was taken to be equal to 2. In the case where the SRI technique was employed, the number of sampling points for the bending term was 3 while that for the shearing term was 1.

Shear Locking
Shear locking is a phenomenon where the beam element is excessively stiff for the range of very small thickness (or the length-to-thickness ratio, L/h, is very large).To observe this phenomenon, a clampedclamped beam with L = 10 m, E = 2000 kN/m 2 , k = 0.84967, v = 0.3, subjected to uniformly-distributed load, q=1 kN/m, was analyzed using eight elements.The height of the beam was varied from thick, L/h=5, up to extremely thin, L/h=10000.(32) Theoretically, as the beam becomes very thin, the solution of Timoshenko beam theory will converge to the solution of Euler-Bernoulli beam theory.The results of the analyses were presented in Fig. 4 for the cases of using full integration, P3-3-G, and selective-reduced integration, P3-3-G (SRI).The figure shows that the element with full integration suffers shear locking while the element with SRI indicates no locking.Thus, the SRI technique is effective to eliminate the shear locking.This finding confirms the conclusion in the previous study [18].Note that the SRI technique is also applicable for Kriging-based rectangular Reissner-Mindlin plate elements [18] but it is not applicable for triangular elements [10].

Accuracy and Convergence
To study the accuracy and convergence of the static solutions of the beam elements using P3-3-G and P3-3-G SRI, the beam subjected to triangulardistributed load (Fig. 3) was analyzed with the parameter E=1000 kN/m 2 , I=0.020833 m 4 , A=1 m 2 , G=384.61 kN/m 2 , k=0.84967, v=0.3, L=4 m, L/h=8.The beam was divided into 2, 4, 8, and 16 elements.The deflections at the free end, the bending moments and the shearing forces at the clamped end were observed.The results were normalized with respect to their respective exact solutions and presented in Tables 3-5 together with the results of the superconvergent beam element of Friedman and Kosmatka (F&K) [12].The tables show that both P3-3-G and P3-3-G (SRI) can produce very accurate displacements and reasonably accurate moments.
The element with P3-3-G gives slightly better results than the element with SRI.The shearing forces directly calculated from the shearing strains (Eq.22), however, cannot produce accurate results.The use of SRI worsen the shearing forces.It is worthy to note that this fact is also true for the traditional Timoshenko beam element [19].It is interesting to notice that the results using full integration for the displacement and moment, converge from below, while those using SRI converge from above.This is because the use of SRI makes the stiffness matrix becomes 'less stiff' than the actual.

Free Vibration
To investigate the accuracy and convergence of the natural frequency solutions from the beam element with P3-3-G and P3-3-G (SRI), free vibrations of a thick and a very thin simply-supported beams were considered.The beam parameters are L=10 m, the width b=1 m, E= 2x10 9 N/m 2 , ν=0.3, ρ=10 kg/m 3 , L/h=5 for the thick beam and L/h=1000 for the very thin beam.The beams were divided into 4, 8, 16, and 32 elements.The resulting natural frquencies were expressed in the form of non-dimensional frequency parameter, i.e.
where λi is the non-dimensional frequency parameter, ωi is the angular natural frequency of the vibration mode-i and m is the mass per unit length.For the thick beam, the convergent solutions of the pseudospectral method [20] were taken to be the reference solutions to assess the accuracy.Whilst for the very thin beam, the exact solutions of the Euler-Bernoulli beam [18] were used as the reference solutions.The normalized results of the analyses using full integration and SRI were presented in Tables 6 and 7 together with the reference solutions.The tables demonstrate that for the case of thick beam the performance of the elements is excellent, both with full and reduced integrations.For very thin beams, a very fine mesh is needed to achieve good accuracy for high mode frequencies.The element with SRI outperforms the element with full integration for the case of thin beam.

Conclusions
Kriging-based Timoshenko beam elements have been developed and tested in the static and free vibration problems.The test results show that the shear locking can be eliminated using the SRI technique.The elements, both with full and reduced integrations, can produce accurate displacement, bending moment, and natural frequencies.For thick and not very thin beams, the element with full integration gives better overall results while for very thin beams, the element with SRI is better.The shearing forces directly calculated from the shearing strains, however, are not accurate.While the SRI technique works well for eliminating the shear locking, the use of this technique worsen the shearing force results.Alternative methods to eliminate the shear locking needs to be studied in the future research.

Figure 2 .
Figure 2. Three-layer DOI of a typical 1D element

Figure 4 .
Figure 4. Normalized deflection at the mid-span of the beam for different length-to-thickness ratiosThe deflection at the mid-span was observed and normalized to the Euler-Bernoulli beam solution, i.e. w Euler = 384EI qL4

Table 1 .
Relative displacement errors for different number of sampling points on stiffness matrix (in this case Nbody = 2)

Table 2 .
Relative displacement errors for different number of sampling points on equivalent nodal force vector (in this case Nsamp = 2)

Table 3 .
Deflections at the free end normalized to the deflection of exact solution

Table 4 .
Bending moments at the clamped end normalized to the bending moment of exact solution

Table 5 .
Shearing forces at the clamped end normalized to the shearing force of exact solution

Table 7 .
Normalized frequency parameter for the very thin beam (L/h=1000)