https://ced. Least-squares Smoothed Shape Functions for Constructing Field-Consistent Timoshenko Beam Elements Wong. Tjahyono. Hartono. 2, and Setiabudi. 1 Associate Professor. Department of Civil Engineering. Petra Christian University Jl. Siwalankerto 121-131. Surabaya 60236. INDONESIA 2 Alumnus. Department of Civil Engineering. Petra Christian University Jl. Siwalankerto 121-131. Surabaya 60236. INDONESIA DOI: https://doi. org/10. 9744/ced. Article Info: Submitted: Jan 18, 2025 Reviewed: Feb 04, 2025 Accepted: Mar 01, 2025 Keywords: Timoshenko beam element, field consistency, least-squares smoothed shape function, shear locking, bifurcation buckling, free vibration. Corresponding Author: Wong. Department of Civil Engineering. Petra Christian University Jl. Siwalankerto 121-131. Surabaya INDONESIA Email: wftjong@petra. Abstract This paper presents an approach for constructing field-consistent Timoshenko beam elements using least-squares smoothed (LSS) shape The variational basis for shear strain redistribution is thoroughly explained, leading to the derivation of LSS shape functions for linear, quadratic, and cubic Timoshenko beam elements. These elements are then applied to linear static analysis, bifurcation buckling analysis, and free vibration analysis of prismatic and tapered beams. Numerical tests demonstrate that the LSS-based beam elements effectively eliminate shear locking and provide accurate, reliable Their performance is comparable to the discrete shear gap technique but with a simpler implementation procedure. The LSS shape function approach offers a practical and efficient alternative for achieving field consistency in Timoshenko beam elements, with potential applications in enhanced finite element methods (FEM. such as isogeometric FEM and Kriging-based FEM. This is an open access article under the CC BY license. INTRODUCTION Beam finite elements are essential in practical structural engineering applications. One of the most widely used theories for developing beam elements is the Timoshenko beam model . , which accounts for shear deformation and cross-sectional rotatory inertia. However. Timoshenko beam elements developed strictly from the standard displacement-based finite element formulation can produce erroneous results, particularly in the case of thin beams. These errors may manifest as very small displacement results compared to the correct solutionAia phenomenon known as shear locking . Aias well as suboptimal convergence or severe oscillations in the shear force distribution . The primary cause of these errors is that the standard finite element formulation, which employs the same interpolation scheme for both deflection and rotation fields, leads to an inconsistent transverse shear strain field . This means that the approximate shear strain is inconsistent with the physical behavior of thin beams, where the shear strain approaches zero. Thus, standard Timoshenko beam elements are unable to represent bending deformation without transverse shear strain . hearless bending deformatio. An early proposed method to make Timoshenko beam elements consistent and hence free from shear locking is the selective reduced integration (SRI) technique . In this technique, the number of Gaussian quadrature sampling points for evaluating the stiffness matrix associated with shear deformation is intentionally reduced from the one Note : Discussion is expected before July, 1st 2025, and will be published in the AuCivil Engineering DimensionAy, volume 27, number 2. September 2025. ISSN : 1410-9530 print / 1979-570X online Published by : Petra Christian University Wong. Tjahyono. Hartono. , and Setiabudi. required for an exact integration. This technique has also been applied in various enhanced finite element methods (FEM. , including Kriging-based FEM . and NURBS (Non-Uniform Rational B-Spline. isogeometric FEM . The SRI technique can effectively eliminate the shear-locking phenomenon in both traditional and enhanced FEMs. However, the resulting shear force is only accurate at the quadrature sampling points. The overall shear force distribution can be even more erratic . than that of the original inconsistent beam element . A more recent method to ensure consistency in Timoshenko beam elements and eliminate shear locking is the discrete shear gap (DSG) technique . In this technique, the inconsistent shear strain field is replaced with a substitute shear strain field derived from the interpolated shear gap. Wong and Sugianto . demonstrated that the DSG technique is effective in ensuring the consistency of linear, quadratic, and cubic Timoshenko beam elements. Furthermore, this technique has also been successfully applied within the frameworks of Kriging-based FEM . and NURBS isogeometric FEM . Although the DSG technique works perfectly to make Timoshenko beam elements consistent, its implementation procedure is quite complicated. A simpler alternative method for achieving field consistency in Timoshenko beam elements is to redistribute the shear strain field using a set of least-squares smoothed (LSS) shape functions for the rotation field. This approach was originally proposed by Prathap and Babu . for eliminating shear locking and/or membrane locking in shear-deformable linear and quadratic beams. However, this simple strain redistribution technique does not appear to be widely recognized . utside Prathap and co-workersAo research grou. , which can be attributed to several reasons. The LSS shape functions in References . were presented without a derivation, and there was no detailed elaboration on the variational principle associated with the strain redistribution. This paper aims to present in detail the variational basis for the redistribution of shear strain and outline a systematic procedure for constructing LSS shape functions in Timoshenko beam elements. The procedure is used to develop a set of LSS shape functions for linear, quadratic, and cubic Timoshenko beam elements. The resulting beam elements are then applied for linear static analysis, bifurcation buckling analysis, and free vibration analysis of prismatic or tapered beams. A series of numerical tests are conducted to evaluate the performance of these field-consistent beam The results are compared to those obtained using the original field-inconsistent beam elements and the beam elements with the DSG technique . Timoshenko Beam Elements Governing Variational Equations Consider a beam of length L that is subjected to a distributed transverse load q = q(X, . and an axial compressive load P = P. , as illustrated in Figure 1. The beam is assumed to be made from a homogeneous and isotropic linear elastic material with a modulus of elasticity E, shear modulus G, and mass density A. According to Timoshenko beam theory, the motion of the beam at any time ycyc Ou 0 can be described using two independent field variables: the deflection of the neutral axis ycyc = ycyc. cUycU, ycy. and the rotation of cross sections yuEyuE = yuEyuE. cUycU, ycy. , 0 O ycUycU O yaya. Figure 1. Coordinate System and the Positive Sign Convention for the Beam Deflection, w. Cross-Section Rotation . Distributed Load q, and Axial Load P The weak form of governing equations for the beam motion that accounts for the effect of the axial force on bending deformation can be written as . yuyuyuyu yuUyuUyuUyuUycycO yccyccyccycc yuyuyuyu yuUyuUyayaycUycU yuEyuEO yccyccyccycc yuyuyuEyuE,ycUycU yayayayaycUycU yuEyuE,ycUycU yccyccyccycc yuyuyuyu ycoycoycoycoycoycoycoyco yccyccyccycc Oe yuyuycyc,ycUycU ycEycEycyc,ycUycU yccyccyccycc = yuyuyuyu ycyc yccyccyccycc , 1 . , for all yuyuyuyu OO Es and for all yuyuyuyu OO Es yay. Vol. No. March 2025: pp. Least-squares Smoothed Shape Functions In this equation, the symbol denotes the variational operator. A and IY denote the cross-sectional area and the crosectional second moment of area about the Y-axis, respectively. In general. A and IY may vary along the length of the The variable k refers to the shear correction factor, which depends on the geometrical shape of the beam. The double dots signify the partial second derivative of the corresponding variable with respect to time t, whereas the comma followed by subscript X signifies the partial derivative of the variable with respect to X. The shear strain is given as . yuyu = ycyc,ycUycU Oe yuEyuE The symbol Es1 . , yay. denotes the Sobolev function space of the first degree in interval 0 < X < L. The bending moment and the shear force . erpendicular to the X-axi. are given as . ycAycA = yayayayaycUycU yuEyuE,ycUycU Finite Element Formulation . ycEycE = ycoycoycoycoycoyco. cyc,ycUycU Oe yuEyuE) Oe ycEycEycyc,ycUycU Suppose the beam is partitioned into Ne number of elements with Np number of nodes. Consider a typical beam element . lement number . possessing n nodes. The displacement and rotation fields over this element are approximated as follows: ycyc OO ycyc Ea = Ocycuycuycnycn=1 ycAycAycnycn . uOyuO)ycycycnycn . = . cAycAycyc ]. yuEyuE OO yuEyuE Ea = Ocycuycuycnycn=1 ycAycAycnycn . uOyuO)yuEyuEycnycn . = . cAycAyuEyuE ]. = . yuEyuE1 . yuEyuE2 . U ycycycuycu . yuEyuEycuycu . }T . cAycAycyc ] = . cAycA1 . uOyuO) 0 ycAycA2 . uOyuO) 0 U ycAycAycuycu . uOyuO) . cAycAyuEyuE ] = . ycAycA1 . uOyuO) 0 ycAycA2 . uOyuO) U 0 ycAycAycuycu . uOyuO)} . In these equations, the superscript h refers to the association of wh and h with a discretization using a mesh of an element characteristic length scale h. Function Ni() represents the element shape functions associated with node i. and i. represent the deflection and rotation at nodal point number i, respectively. The shape functions are expressed in terms of natural coordinate . Oe1 O yuOyuO O 1. The shape functions for a two-node linear element, a three-node quadratic element, and a four-node cubic element are respectively given as . ssuming nodes 1 and 2 are the end nodes and the rests are the interior node. ycAycA1 = 12. Oe yuOyuO), ycAycA2 = 12. yuOyuO) ycAycA1 = Oe12yuOyuO. Oe yuOyuO), ycAycA3 = 1 Oe yuOyuO 2 , ycAycA2 = 12yuOyuO. yuOyuO) . Oe yuOyuO). Oe 9yuOyuO 2 ), ycAycA3 = ycAycA1 = Oe16 Oe 3yuOyuO). Oe yuOyuO 2 ), . yuOyuO). Oe 9yuOyuO 2 ) ycAycA4 = 16 . 3yuOyuO). Oe yuOyuO 2 ), ycAycA2 = Oe16 ycUycU = Ocycuycuycnycn=1 ycAycAycnycn . uOyuO)ycUycUycnycn . yaya = ycUycU,yuOyuO . The mapping from natural coordinate to the global coordinate X is given as In this equation. Xi is the global coordinate of nodal point number i and J is the Jacobian of the mapping. To derive the finite element equations, the integrals in the variational equation. Equation . , are firstly expressed as the sum of the integrals over each element interval. Subsequently, substituting Equations . into the resulting equation and applying the standard finite element formulation yield the discretized system of equations as follows: cAycA]yayaO . Oe ycEycEyayag . } = . } . where [M], [K], [K. are the discretized structural mass, stiffness, and geometric stiffness matrices, respectively. The vector . } is the vector of global nodal displacement, that is, {D. } = ycyc1 . yuEyuE1 . yuEyuE2 . U ycycycAycAycyycy . yuEyuEycAycAycyycy . Vol. No. March 2025: pp. Wong. Tjahyono. Hartono. , and Setiabudi. The vector {F. } is the vector of global nodal force. The structural mass, stiffness, and geometric stiffness matrices, as well as the global nodal force vector, are obtained by assembling the element mass matrix, . e, the element stiffness matrix, . e, the element geometric stiffness matrix, . e, and the element nodal force vector, . e, respectively, for all elements . , element number e = 1 to N. These element matrices are defined as follows: yeIyeI = OOe1. cAycAycyc ]T yuUyuUyuUyuU . cAycAycyc ] yayayayayaya OOe1. cAycAyuEyuE ]T yuUyuUyayaycUycU . cAycAyuEyuE ] yayayayayaya yeIyeI = . yceyceb = OOe1. cAycAyuEyuE ],TyuOyuO yayayayaycUycU . cAycAyuEyuE ],yuOyuO yayaOe1 yccyccyccycc yceyces = OOe1. cAycAycyc ],yuOyuO yayaOe1 Oe . cAycAyuEyuE ] ycoycoycoycoycoyco. cAycAycyc ],yuOyuO yayaOe1 Oe . cAycAyuEyuE ] yayayayayaya yceyce ycoycog = OOe1. cAycAycyc ],TyuOyuO . cAycAycyc ],yuOyuO yayaOe1 yccyccyccycc The element nodal force vector is defined as . yceyce = OOe1. cAycAycyc ]ycNycN ycyc yayayayayaya = . Oe ycEycEyayag . = . The discretized equations for linear static, free vibration, and bifurcation buckling problems can be derived from Equation . by reducing it, respectively, to . cAycA]yayaO . } = . Least Squares Smoothed Shape Functions As mentioned in the introduction section, the standard interpolations for the Timoshenko beam field variables, i. Equations . , result in a discretized shear strain, h, that is inconsistent with the physical constraint of vanishing the shear strain. Equation . , when the beam becomes infinitely thin . he Kirchhoff constrain. This inconsistency in the discretized shear strain leads to shear locking, poor convergence, and severe stress oscillation commonly observed when using standard displacement-based Timoshenko beam elements . ith exact calculation of all integral. To develop beam elements that are consistent with the Kirchhoff constraint. Prathap and Babu . proposed a set of least-squares smoothed shape functions for interpolating the rotation field in the shear strain expression. This section addresses the variational basis and a detailed derivation of the smoothed shape functions. Variational Basis The variational basis for the development of a field-consistent Timoshenko beam element is a modified HellingerReissnerAos variational principle . In this principle, the discretized strain energy for a Timoshenko beam element of length Le is expressed in the following form: ycOycO Ea = O0 yayayayaycyc yuEyuE,Eaycuycu yccyccyccycc Oe O0 yayayceyce ycoycoycoycoycoycoyuyuI 2 yccyccyccycc O0 ycoycoycoycoycoycoyuyuI yuyu Ea yccyccyccycc where h and h are the discretized rotation and the discretized kinematic shear strain, respectively, while yuyuI denotes an assumed shear strain. The variable x is the element coordinate, starting from the left end of the beam element. variation of the functional in Equation . with respect to yuyuI yields yayayceyce O0 yuyuyuyuI ycoycoycoycoycoycoyuyu Ea Oe yuyuI yccyccyccycc = 0 . For a prismatic beam element made from a homogeneous material, kGA is constant along the element, hence Equation . simplifies to yayayceyce O0 yuyuyuyuI yuyu Ea Oe yuyuI yccyccyccycc = 0 . Vol. No. March 2025: pp. Least-squares Smoothed Shape Functions Equation . is known as the orthogonality condition . This orthogonality condition can be utilized to determine a consistent assumed shear strain field, yuyuI , from the already known kinematically derived inconsistent h. Least Squares Smoothed Shape Functions A method to achieve consistency in the shear strain field is to redistribute the kinematic shear strain through a leastsquares smoothing approach . This process begins by defining the assumed shear strain field and its variations as follows: yuyuI = ycyc,Eaycuycu Oe yuEyuEI yuyuyuyuI = OeyuyuyuEyuEI O0 yuyuyuEyuEI yuEyuE Ea Oe yuEyuEI yccyccyccycc = 0 where yuEyuEI is the unknown substitute rotation field. Substituting these equations into Equation . results in yayayceyce This equation can be interpreted as the stationary condition of the following functional: yceyce yaya 1 uEyuEI ) = O0 yuEyuE Ea Oe yuEyuEI yccyccyccycc Thus, the orthogonality condition of the rotation, expressed in Equation . , requires that the substitute rotation field, yuEyuEI, is a least-squares equivalent of the discretized rotation field, h. Next, the substitute rotation field is defined as follows: uOyuO)yuEyuEycnycn . = . cAycA yuEyuE ]. yuEyuEI = Ocycuycuycnycn=1 ycAycA uOyuO), i = 1. A, n, represent substitute shape functions. Substituting Equations . into Equation . where ycAycA yceyce yuEyuE ]) = 1 . T Oyaya . cAycAyuEyuE ] Oe . cAycA yuEyuE ]T . cAycAyuEyuE ] Oe . cAycA yuEyuE ]yccyccyccycc . cAycA This equation indicates that the least-squares substitute rotation field can be obtained by using the least-squares substitute shape functions. Furthermore, to achieve field consistency, these shape functions must be chosen to be a one-order lower polynomial than the order of the shape functions for approximating the deflection. These consistent shape functions are referred to as LSS shape functions. An example of algebraic calculation for obtaining the LSS shape functions for a quadratic beam element is as follows. The general form of quadratic shape functions is given by: ycAycAycnycn () = ycayca ycaycaycayca ycaycayuOyuO 2 ycnycn () = yuyu yuyuyuyu ycAycA ycnycn ()2 yccyccyccycc ycnycn ) = O1 1 ycAycAycnycn () Oe ycAycA cAycA Oe1 2 The LSS shape functions are chosen to be linear . ne degree lower than quadrati. , that is. The coefficients and are determined by minimizing the functional: Substituting Equations . into Equation . transforms the functional into a function of two variables, and , viz. uyu, yuy. = OOe1 . cayca Oe yuy. cayca Oe )yuOyuO ycaycayuOyuO 2 yccyccyccycc The unknown coefficients and can be found using the standard calculus method to locate the extreme values of f(, ), that is, yuiyuiyuiyui yuiyuiyuiyui = 0 and . yuiyuiyuiyui Vol. No. March 2025: pp. yuiyuiyuiyui Wong. Tjahyono. Hartono. , and Setiabudi. Solving these equations yields yuyu = ycayca ycayca and yuyu = ycayca Thus, the resulting LSS shape function is ycnycn () = ycayca 1 ycayca ycaycaycayca ycAycA Using Equation . , the LSS shape functions associated with node numbers 1, 2, and 3 for a quadratic beam element can be determined. For instance, the quadratic shape function corresponding to node number 1, refers to Equation . , has coefficients ycayca = 0, ycayca = Oe 2, ycayca = 2. Substituting these values into Equation . yields the LSS shape function associated with node number 1, that is, 1 () = 0 1 1 Oe 1 yuOyuO = 1 1 Oe yuOyuO ycAycA This result is the same as what is presented in Ref. By applying the abovementioned procedure, a set of LSS shape functions for linear, quadratic, and cubic beam elements can be derived. The results are respectively as follows: 2 = 1 1 = 1, ycAycA ycAycA 3 = 2 , ycAycA 2 = 1 1 yuOyuO 1 = 1 1 Oe yuOyuO , ycAycA ycAycA 2 = Oe 1 1 Oe 22 yuOyuO Oe 9yuOyuO 2 4 = 9 1 6 yuOyuO Oe yuOyuO 2 , ycAycA ycAycA ycnycn ()ycAycAycnycn () Oe ycAycA ycnycn ()yccyccyccycc = 0 OOe1 yuyuycAycA 3 = 9 1 Oe 6 yuOyuO Oe yuOyuO 2 , 1 = Oe 1 1 22 yuOyuO Oe 9yuOyuO 2 , ycAycA ycAycA The procedure to determine the unknown coefficients in the LSS shape functions described here involves the direct ycnycn ), as expressed in Eq. Alternatively, the coefficients can be determined by first use of the functional, yaya. cAycA ycnycn , which leads to the following equation: invoking the stationarity of the functional with respect to ycAycA After this, the original shape function and the LSS shape function are substituted into this integral equation. The unknown coefficients in the LSS shape function can then be found by solving the resulting system of linear equations. It is important to note that, to construct field-consistent Timoshenko beam elements, the LSS shape functions for the rotation field are used to replace the original shape functions solely in the expression corresponding to the shear This specifically applies to the shape functions in the shearing stiffness matrix, as expressed in Equation . Meanwhile, the shape functions for the rotation field in the bending stiffness matrix (Equation . ) and the mass matrix (Equation . ) remain unchanged. Numerical Tests The linear, quadratic, and cubic Timoshenko beam elements with the original shape functions (Equations . , . , and . ) and LSS shape functions (Equations . , . , and . ) have been implemented in Matlab. These elements are referred to as AoOriginal SFAo and AoLSS SFAo, respectively. They have been tested and applied in static, buckling, and free vibration analyses of prismatic as well as tapered beams. Additionally, for comparison purposes, results obtained using locking-free Timoshenko beam elements with the discrete shear gap technique . eferred to as AoDSGA. are also included in the following report. The shear modulus and shear correction factor of a beam were calculated using the given modulus of elasticity. E, and Poisson's ratio, , with the following formula . yuOyuO) yaya = 2. yuOyuO) , ycoyco = 12 11yuOyuO Vol. No. March 2025: pp. Least-squares Smoothed Shape Functions Table 1. Minimum Number of Quadrature Sampling Points for Prismatic and Tapered Timoshenko Beam Elements with Original Shape Functions of Different Orders Element Mass Linear Quadratic Cubic Bending Stiffness Shear Stiffness Prismatic beam Tapered beam2 Geometry Force1 Linear Quadratic Cubic Assuming that q is linearly distributed The cross-sectional area and moment of inertia are interpolated using the shape functions. For tapered beams, the cross-sectional area. A = A(X), and moment of inertia. IY = IY(X), for each beam element are interpolated from their values at the nodes. The integrals in the element matrix expressions. Equations . , were evaluated using the Gauss quadrature The minimum number of quadrature sampling points required for each matrix to achieve exact integration results for the beam elements with original shape functions is presented in Table 1. For the beam elements with LSS shape functions, the number of sampling points required to evaluate the shear stiffness matrix can be reduced by one. Static Analysis Investigation of Shear Locking A fixed-fixed supported beam is used to detect shear locking. The geometrical, material, and loading parameters are as follows: L = 10 m, b = 1 m. E = 10y106 kN/m2, = 0. 3, and q = Oe1 kN/m. The beam is discretized using eight elements of equal length, as shown in Figure 2. The length-to-thickness ratio of the beam varies from L/hB = 5 . epresenting a thick bea. , to 10, 100, 1000, and 10000 . epresenting an extremely thin bea. Although a beam with L/hB greater than 100 is outside the practical range for length-to-thickness ratios, it is included in this test to detect shear locking as the beam becomes extremely thin. Finite element model of the beam . Beam cross section Figure 2. Finite Element Model of a Fixed-Fixed Supported Beam Subjected to a Uniform Load The analytical solution for the beam mid-span deflection is given by yaya = 384yayayaya 8yayayayayaya The mid-span deflection results obtained using different Timoshenko beam elements, normalized to the exact solution, are presented in Table 2. The table indicates that the linear, quadratic, and cubic Timoshenko beam elements with the LSS shape functions are free from shear locking. The results are the same as those obtained using beam elements with the DSG technique. contrast, the linear beam element with the original shape functions experiences shear locking. Although the quadratic original SF element does not suffer from shear locking, it is not as accurate as the quadratic LSS-based element. All cubic beam elements are unaffected by shear locking and can provide exact results for mid-span deflection. Vol. No. March 2025: pp. Wong. Tjahyono. Hartono. , and Setiabudi. Table 2. Normalized Mid-Span Deflections of a Fixed-Fixed Supported Beam with Different Length-to-Thickness Ratios (L/H. Obtained using Various Beam Elements Element Type L/hB Linear Beam Elements Quadratic Beam Elements Cubic Beam Elements LSS SF DSG Original SF LSS SF DSG Original SF LSS SF DSG All results are exact. Original SF LSS SF: Timoshenko beam elements using the LSS shape functions. DSG: Timoshenko beam elements with the discrete shear gap technique. The results were taken from Ref. Original SF: Timoshenko beam elements using the original shape functions. Assessment of Accuracy and Convergence To evaluate the accuracy and convergence characteristics of Timoshenko beam elements, the fixed-fixed supported beam with a length-to-thickness ratio of 10 is discretized using different numbers of equal elements: Ne = 4. 8, 16, The analysis results for the mid-span defection, fixed-end moment, and fixed-end shear force, normalized to their corresponding exact values (Equations . ), are presented in Tables 3. Table 3. Normalized Mid-Span Deflections. Fixed-End Bending Moments, and Fixed-End Shear Forces of a Fixed-Fixed Supported Beam of L/Hb = 10 for Different Number of Elements. Ne. Obtained using Different Beam Elements . Normalized Results using Linear Beam Elements LSS SF Deflection DSG Original SF Bending Moment LSS SF DSG Original SF LSS SF Shear Force DSG Original SF Normalized Results using Quadratic Beam Elements LSS SF Shear Force DSG Original SF Normalized Results using Cubic Beam Elements Deflection Bending Moment LSS SF DSG Original SF LSS SF DSG Original SF LSS SF All results are exact. Shear Force DSG Original SF LSS SF Deflection DSG Original SF ycAycA. exact = Bending Moment LSS SF DSG Original SF ycycyaya2 , ycEycE. exact = ycycyaya Table 3 shows that all results converge to the corresponding analytical values. As expected, higher-order Timoshenko beam elements yield more accurate results compared to lower-order elements. The results of the LSS-based and DSG elements are identical, except for the fixed-end moments obtained from the cubic beam elements, where some minor discrepancies are observed. The performance of the LSS-based element is superior to that of the original SF element. Vol. No. March 2025: pp. Least-squares Smoothed Shape Functions Moreover, the cubic LSS-based element provides exact results for mid-span deflection, fixed-end moment, and fixedend shear force, even when using a minimal number of elements. Bifurcation Buckling Analysis The Timoshenko beam elements are now being used for the bifurcation buckling analysis of both prismatic and tapered beams. The same fixed-fixed supported beam as in the static analysis, with L/hB = 10, is considered. Additionally, the beam has been modified to have a tapered shape, as illustrated in Figure 3. Two tapering ratios are considered, namely c = 0. 5 and c = 0. efer to Figure 3 for the definition of . Sect. A-A L = 10 m. EaB0 = 1 m, b = 1 m ycuycu Ea. = Ea0 1 Oe ycayca yaya , c: tapering ratio Figure 3. Tapered Fixed-Fixed Supported Beam The resulting critical buckling loads, normalized to their respective reference values, are presented in Tables 4. The reference value for the prismatic beam is calculated using the following analytical formula (. as cited in . ycEycEcr = yuUyuU2 yayayayaycUycU yuUyuU2yayayaya ycUycU yaya2 eff yayayayas where the effective buckling length for the fixed-fixed supported beam is Leff = L/2. For the tapered beams, since an analytical solution is unavailable, the reference values are obtained from finite element analysis results using 48 cubic LSS-based elements. These values are taken as the reference solution because of the excellent performance of the cubic-LSS-based beam element and additionally, the use of 48 cubic elements represents a very fine finite element These reference values are listed in the last row of Table 4. The findings from the static analysis can be applied to the bifurcation buckling analysis. The results converge to their reference values. The performance of the LSS-based elements is superior to that of the original beam elements. The accuracy of both linear and quadratic LSS-based elements is the same as that of the corresponding elements with the DSG technique, while the cubic elements show a slight discrepancy in accuracy. Free Vibration Analysis The LSS-based Timoshenko beam elements are utilized to determine the first eight free-vibration frequencies of a tapered fixed-fixed supported beam of L/hB = 10 with the tapering ratio of c = 0. 5, which has been considered in the bifurcation buckling analysis. The mass density is taken as A = 1000 kg/m3. The analysis is conducted using 16 beam As in the bifurcation buckling analysis, the natural frequencies . n the unit of H. obtained with 48 elements of the cubic LSS-based beam element are used as the reference solution because of the unavailability of analytical The results are listed in Table 5. The table demonstrates the superiority of the field-consistent LSS-based beam elements compared to the original inconsistent beam elements when predicting natural frequencies, particularly for linear and quadratic beam elements. As expected, higher-order beam elements provide greater accuracy than the lower-order elements. Consistent with the findings in static and buckling analyses, the results for the linear and quadratic LSS-based and DSG beam elements are identical, while those for the cubic beam elements show slight differences. Vol. No. March 2025: pp. Wong. Tjahyono. Hartono. , and Setiabudi. Table 4. Normalized Critical Buckling Loads of the Prismatic (C = . and Taped Fixed-Fixed Supported Beams of L/Hb = 10 for Different Numbers of Elements. Ne. Obtained using Different Beam Elements . Normalized Results using Linear Beam Elements Tapering ratio = 0 Tapering ratio = 0. Tapering ratio = 0. Nel LSS SF & DSG Original SF LSS SF & DSG Original SF LSS SF & DSG Original SF Nel . Normalized Results using Quadratic Beam Elements Tapering ratio = 0 Tapering ratio = 0. Tapering ratio = 0. LSS SF & DSG Original SF LSS SF & DSG Original SF LSS SF & DSG Original SF Nel Pcr exact . Normalized Results using Cubic Beam Elements and the Exact Critical Buckling Loads Tapering ratio = 0 Tapering ratio = 0. Tapering ratio = 0. LSS SF DSG Original SF LSS SF DSG Original SF LSS SF DSG Original SF E 05 kN E 05 kN Table 5. The First Eight Natural Frequencies of the Taped Fixed-Fixed Supported Beams of L/Hb = 10, which are Obtained using 16 Beam Elements of Different Formulations. Normalized to the Reference Frequencies Normalized natural frequency Reference Mode Linear element Quadratic element Cubic LSS Original LSS Original Original (H. LSS SF DSG DSG DSG *Obtained using 48 elements of the cubic LSS SF element Remarks All the examples presented are based on the presumption that the interior element nodes are located at their AonaturalAo positions in the Cartesian coordinate system. For a quadratic beam element, this means the nodes are at the midpoints, while for a cubic beam element, they are at the one-third points. If the interior element nodes do not lay at these natural positions, the use of the LSS shape functions cannot provide a consistent shear strain field distribution, leading to very poor results . hear lockin. Prathap and Naganarayana . proposed a method to modify the use of the LSS shape functions for a quadratic Timoshenko beam element so that it can maintain the shear strain field CONCLUSIONS The variational foundation for the redistribution of shear strain in Timoshenko beam elements has been presented. This theoretical groundwork leads to a systematic procedure for deriving LSS shape functions. The LSS shape functions for linear, quadratic, and cubic Timoshenko beam elements have been derived and are expressed explicitly. These shape functions ensure field consistency within the beam elements. The LSS-based Timoshenko beam elements were applied to linear static analysis, bifurcation buckling analysis, and free vibration analysis of prismatic and tapered beams. The numerical tests demonstrated that the LSS-based beam Vol. No. March 2025: pp. Least-squares Smoothed Shape Functions elements are free from shear locking and yield accurate and reliable results. Their performance is practically the same as that of the consistent beam elements with the DSG technique. however, the implementation for the LSS shape functions is simpler than that of the DSG technique. Further research could explore the development of LSS shape functions for plate and shell finite elements. Additionally, the potential application of the LSS approach to enhanced finite element methods, such as isogeometric analysis and Kriging-based FEM, warrants investigation. REFERENCES