Study of the Discrete Shear Gap Technique in Timoshenko Beam Elements

A major difficulty in formulating a finite element for shear-deformable beams, plates, and shells is the shear locking phenomenon. A recently proposed general technique to overcome this difficulty is the discrete shear gap (DSG) technique. In this study, the DSG technique was applied to the linear, quadratic, and cubic Timoshenko beam elements. With this technique, the displacement-based shear strain field was replaced with a substitute shear strain field obtained from the derivative of the interpolated shear gap. A series of numerical tests were conducted to assess the elements performance. The results showed that the DSG technique works perfectly to eliminate the shear locking. The resulting deflection, rotation, bending moment, and shear force distributions were very accurate and converged optimally to the corresponding analytical solutions. Thus the beam elements with the DSG technique are better alternatives than those with the classical selective-reduced integration.


Introduction
Structural models of beam, plate, and shell are the earliest problems to which the finite element method was developed and applied in the early 1960s [1].In the early development, the classical theory of beam, plate, and shell, which neglects shear deformation, was used as the basis to develop a variety of beam, plate, and shell elements.Subsequently finite element developers preferred to use the sheardeformable theory as the basis to develop beam, plate, and shell elements since this theory requires a low order continuity requirement, i.e.C 0 continuity, and is more general than the classical theory.A major difficulty in formulating a finite element based on the shear-deformable theory is a phenomenon in which an element becomes excessively stiff as the element becomes very thin, which is well-known as shear locking.
Literatures on finite elements of shear deformable beam, plate, and shell contain innumerable concepts, techniques, or tricks to overcome the difficulty of shear locking.
A well-known early technique is the selective reduced integration technique [2,3], in which the order of numerical integration for calculating the shear related part of the stiffness matrix is deliberately reduced.While this simple technique works well for Timoshenko beam elements, it produces a quadrilateral Reissner-Mindlin plate bending element that contains spurious zero strain energy modes [3] and it is not applicable to a triangular plate bending element [4].Other wellknown techniques are the assumed natural shear strain method [5,6], the discrete Kirchhoff approaches [7][8][9], and the mixed interpolation of tensorial component approach [10].From the literature it is recognized that a successful locking-free technique for a quadrilateral element is usually cannot be directly applied to a triangular element and vice versa.
Recently Bletzinger et al. [11] presented a unified approach for developing locking-free finite elements for shear deformable beams, plates, and shells, the so-called discrete shear gap (DSG) technique.In this technique the troublesome kinematic shear strain fields are replaced with substitute shear strain fields determined from the interpolated shear gaps at the element nodes.A distinctive advantage of this technique is that it is applicable to both triangular and rectangular elements of arbitrary polynomial degrees, with one simple technique.The DSG technique was subsequently generalized and applied to remove membrane locking [12].It also has successfully applied in the context of an alternative finite element formulation called the smoothed FEM [13][14][15].
Although the DSG technique was firstly introduced in the context of Timoshenko beam [11], subsequent applications of the technique were focused on the plate and shell problems.In the context of the Kriging-based finite element method [16], the application of the DSG technique has been studied by Wong and Sulistio [17].It was found that the DSG is effective to eliminate shear locking for a Kriging-based element with cubic polynomial basis, but not effective for those with quadratic and linear bases.In the context of the standard FEM, however, there has been no formulation and numerical test of the DSG technique on higher-order Timoshenko beam elements.It is thus the aim of this paper to study the application of DSG technique in the standard linear, quadratic, and cubic Timoshenko beam elements.The application of the DSG is accomplished by replacing the shear strain-displacement matrix with the one derived based on the DSG concept.A series of numerical tests is carried out to evaluate the performance of the DSG beam elements, in particular to evaluate the effectiveness of the DSG technique at eliminating shear locking.The comparison is made with beam elements using the selective-reduced as well as full integrations.

Governing Equations
A beam model of length L with the cross sectional area A and moment of inertia I is considered.The material properties are a constant modulus elasticity E and a constant Poisson"s ratio ν.The beam is subjected to a static transversal distributed load q=q(x).The displacement of any point in the beam can be expressed in terms of two independent field variables, i.e. the vertical displacement of the neutral axis, or, the deflection w=w(x) and the cross section rotation β= β(x).The coordinate system and sign convention for positive w and β are depicted in Figure 1.
The weak form of the governing equation for static deformation of the Timoshenko beam model is given as: for all admissible δw and δβ (1a) In this equation, δw and δβ are a virtual deflection and a virtual cross section rotation, respectively.The term "admissible" means that the δw and δβ satisfy the requirements of C 0 continuity and homogenous displacement boundary conditions.The constant G is the shear modulus given as and k is a shear correction factor that is dependent upon the cross-section geometry and Poisson ratio [18].The comma signifies the first derivative with respect to the variable next to it (i.e.x).
The bending moment M and shear force Q along the beam are given as ( , )

Finite Element Discretization
Let the beam is subdivided into Nel elements and consider a typical element having n nodes.The deflection and rotation along the element are approximated as follows:   in which [Nw] is the shape function matrix for the deflection, i.e.
  and {d} is the unknown nodal displacement vector, i.e.

   
The number of element node n depends on the polynomial degree employed to approximate w and β, that is, n=2 for a linear element, n=3 for a quadratic element, and n=4 for a cubic element, as shown in Figure 2.
Using the isoparametric mapping concept, the coordinate x is expressed in terms of natural coordinate [ 1, 1]   in the same manner as the approximation for the field variables, i.e.
where xi, i=1,…,n are the nodal coordinates.The expression for the shape functions are listed in Table 1 in terms of ξ.The Jacobian of the mapping is defined as Figure 2. One-dimensional Linear, Quadratic, and Cubic Elements Table 1.Shape Functions for a 1D Element of Two to Four Nodes (taken from [10], p. 343) with the Numbering Labels as Shown in Figure 2.
If the element has node 3 If the element has nodes 3 and 4 1) 1) Substituting Equations ( 4) and (5) into Equation (1) and carrying out the standard finite element formulation procedure yield the element system of equations In this equation, the element stiffness matrix [k] is made up of two parts, namely the stiffness associated with bending deformation, [k]b, and the stiffness associated with shear deformation, [k]s,   In Equations ( 13) and ( 14), The element nodal force vector {f} is given as

Numerical Integration
In the present work the integrals in Equations ( 13), (14), and ( 17) are evaluated numerically using Gaussian quadrature rule.The number of inte-gration sampling points to obtain an exact integration result for an element with linear, quadratic, and cubic shape functions are listed in Table 2, assuming q is a linear function.However, it is well known that the use of full integration results in elements that are excessively stiff when they are applied to thin beams, or their convergence rates are not optimal [3,10,19,20] This phenomenon is known as "shear locking".A simple and efficient technique to eliminate the shear locking is to reduce the number of sampling points for evaluating the shearing part of the stiffness matrix [3].This technique is referred to as selective-reduced integration (SRI) technique.
It is worth mentioning here that while the SRI technique is effective to eliminate shear locking, the resulting shear forces are correct only at the integration sampling points.The shear forces evaluated directly at the element nodes are in great error.The shear force distribution fluctuates heavily [20].

Discrete Shear Gap Technique
The key idea of the DSG technique is to replace the troublesome discretized shear strain, , a substitute shear strain determined from interpolation of "shear gaps" at element nodes.Shear gap at a point x is defined as the increment of displacement due to shear deformation from a reference point x0 [11], that is The shear gap can be interpreted as "the difference between the increase of the actual displacement Δw and the displacement 0 b () which correspond to pure bending" [11].This interpretation is illustrated in Figure 3.
A discrete shear gap (DSG) is the shear gap at an element nodal point with coordinate xi, which can be calculated using Equation ( 18), 0 0 () 57 A substitute shear gap field is defined as the interpolation of the nodal shear gaps, The substitute shear strain is obtained by differentiating Equation ( 21), where Taking node number 1 as the reference point and substituting Equation ( 5) into Equation (20), the DSGs at nodes i=1, 2, …, n can be expressed as It is obvious from this equation that the DSG at node 1 is zero, i.e.
1 0 w   .The DSGs at other nodes are given in Table 3.
Using Equation (25) and Table 3, the matrix of discrete shear gaps, Equation (24), can be expressed in terms of the nodal displacement vector as The expression of [BDSG2] for the linear beam element is Similarly, we can write [BDSG2] matrices for the quadratic and cubic beam elements (they are, however, too long to be included in this paper).Substituting Equation (26) into Equation ( 22) yields Thus the DSG concept can be implemented in a computer code by replacing ) matrix in Equation ( 14) with [BDSG] defined in Equation (28).()

Numerical Tests
The linear, quadratic, and cubic beam elements with the DSG technique have been coded using Matlab [21].For comparison, the standard beam elements with the SRI and full integration have also been coded.
To evaluate the performance of the elements, a series of numerical tests were carried out.The tests included pure bending test, investigation of shear locking, and assessment of the accuracy and convergence characteristics.In all tests, the test problems were analyzed using the linear, quadratic, and cubic DSG beam elements and the standard beam elements with the SRI and full integration.The exception was on the assessment of accuracy through deflection, rotation, bending moment and shear force profiles, in which only the results obtained using the DSG beam elements were presented.The full integration scheme (Table 2) was used to evaluate the integrals in the bending stiffness, shear stiffness, and force terms of the DSG beam elements.The integrals in [BDSG2] matrix, Equation (26), were also numerically evaluated using full integration to obtain exact integration results.For the quadratic and cubic beam elements, the positions of interior nodes were located in their

Pure Bending Test
The first test problem is a cantilever beam with rectangular cross section of the size bh  subjected to an end moment.The beam was discretized using two different meshes, that is, a mesh with four elements of equal length and a mesh with four elements of unequal length shown in Figure 4 [17].The material, geometrical, and loading parameters were taken as follows: L = 10 m, b = 2 m, E = 2000 kN/m 2 , ν = 0.3, and M = 1 kN-m.Two different thickness was considered, that is, h = 2 m (L/h = 5), which represents a thick beam, and h=0.001 m (L/h=10000), which represents an extremely thin beam.Certainly, the beam with L/h=10000 is outside the range of practical length-to-thickness ratio.The purpose to include it here is merely to test the tendency of the beam elements to shear locking.

Figure. 4. Meshes of the Cantilever Beam
The beam is in the state of pure bending with constant bending moment M and zero shear stress along the beam.Therefore this problem may be regarded as the patch test for the beam elements.A beam element is said to pass the test if it can reproduce the exact deflection, rotation, bending moment, and shear force for any mesh.The exact deflection and rotation at the right end are given as: The results of the deflection and rotation at the free end, and the bending moment and shear force at the fixed end were observed and compared to the corresponding exact solutions.For the sake of brevity, only the results for the extremely thin beam using the irregular mesh were presented here in Table 4.
It is seen that even for the most critical case (i.e.very thin beam and irregular mesh), the DSG beam elements are not only free from shear locking but also are able to reproduce exact results.In contrast, the results obtained using the linear beam element with full integration are erroneous because of shear locking.While the SRI technique can eliminate the locking, as addressed by Prathap [20], the shear forces calculated directly at the fixed-end node using Equation ( 3) are in great error, even worse that the shear force obtained with full integration.The overall results showed that that all DSG beam elements passed this pure bending test for both the thick and extremely thin beams, with the regular and irregular meshes.The DSG technique is better than the SRI because the former can give accurate shear force along the beam while the later can only give accurate shear force at the integration sampling points.

Investigation of Shear Locking
A fixed-fixed supported beam subjected to a uniformly distributed load (Figure 5) was analyzed to investigate the effectiveness of the DSG in eliminating shear locking..The fixed-fixed supported beam was chosen because in the authors" experience, this problem is the most critical problem for a shear locking test.The material and geometrical data were the same as in the pure bending test with the uniformly-distributed load q=1 kN/m.The length-tothickness ratio, however, was varied from L/h=5 to L/h=10000.The exact deflection at the mid-span is given as:   The beam was modeled using eight elements of equal length.The resulting deflections at the midspan was observed and normalized to the exact deflection, Equation 31.The results were presented in Table 5.It is seen that all of the DSG beam elements are free from locking.The results are the same as those obtained using the elements with the SRI technique.As expected, the linear beam element with full integration suffers from shear locking so that the results tend to zero as the beam becomes thin.While the quadratic beam element with full integration does not lock, its accuracy is not optimal in the sense that it is just as accurate as the linear beam element with the DSG or SRI (a quadratic beam element is expected to be more accurate than the linear element).

Assessment of Accuracy and Convergence
A cantilever beam subjected to a linearly-distributed load as shown in Figure 6 is utilized to assess the accuracy and convergence characteristics of the DSG beam elements.The geometrical, material, and loading data were L=4 m, h=0.5 m, b=2 m, E=1000 kN/m 2 , v= 0.3, and q0 = 1 kN/m.The exact deflection at the free end, bending moment and shear force at the fixed end are given as [22] The beam were modeled using 1, 2, 4, and 8 elements.The analysis results, normalized to the corresponding exact values, were presented in Table 6.It is seen that for the linear and quadratic beam elements, the deflections and bending moments obtained using the DSG and SRI techniques are identical and always better than those obtained using the beam elements with full integration.The shear strains, however, are not identical; those obtained using the DSG techniques are much more accurate than those obtained using the SRI and full integration beam elements.For the cubic beam elements, as expected, all deflection results are exact since the homogeneous part of Timoshenko beam differential equation is a cubic polynomial function [22].The bending moments obtained using the DSG and SRI cubic elements are not identical but their accuracy is comparable.The shear strains obtained using the DSG cubic element are exact, even though using only one element.Regarding the convergence, it is seen that all beam elements have excellent convergence behavior.
To further assess the accuracy of the DSG beam elements, the profiles of deflection, rotation, bending moment, and shear force obtained using four elements were plotted in Figures 7-10 together with the corresponding exact profiles.The figures demonstrated the exceptional accuracy of the elements in predicting the deflection, rotation, bending moment, and shear force along the beam.The bending moment and shear force obtained using the linear DSG beam element, however, are only exact at the element midpoints.This is because the resulting moment and shear force distributions are, as expected, piecewise linear.It is worth mentioning here that if we use the beam elements with full and selective reduced integrations, the resulting shear force will be very inaccurate [20].Thus, Figure 10 again demonstrates the superiority of the DSG technique in predicting the shear force distribution.

Conclusions
The DSG technique has been applied to the linear, quadratic, and cubic Timoshenko beam elements.The essence of this technique is to replace the troublesome displacement-based shear strain field with a substitute shear strain field obtained from the derivative of the interpolated shear gap.The numerical test results showed that the beam elements with the DSG technique pass the pure bending test, work perfectly in eliminating the shear locking, and can give exceptionally accurate solutions.The distinctive advantage of the DSG technique over the classical SRI is that with the DSG technique a beam element can produce a very accurate shear force distribution, while the with the SRI the shear force distribution will be erroneous.The successful application of the DSG technique in the present beam elements may partly explain the reason of its success in formulating locking-free shell elements [11].Future research may be directed to application of the DSG technique to curve beam elements and geometrically nonlinear analysis of beams.

Figure 1 .
Figure 1.Coordinate System and Positive Directions for the Deflection and Rotation

Figure 3 .
Figure 3. Shear Gap Interpreted as the Difference between the Increase of the Total Displacement and the Increase of the Displacement due to Bending[11] Unequal-length four element model

Figure 5 .
Figure 5. Fixed-fixed Beam Subjected to Uniform Load and its Finite Element Model

Figure 6 .
Figure 6.Cantilever Beam Subjected to Linearly-distributed Load Beam cross section (a) Finite element model of the beam

Figure 7 .Figure 8 .
Figure 7. Deflection Profile Obtained using Four Linear, Quadratic, and Cubic DSG Beam Elements

Figure 9 .Figure 10 .
Figure 9. Bending Moment Profile Obtained using Four Linear, Quadratic, and Cubic DSG Beam Elements

Table 2 .
Number of Integration Sampling Points Required for Full and Selective-reduced Integrations.

Table 3 .
Discrete Shear Gaps at Nodes other than Node 1 for Different Order of 1D Elements

Table 4 .
Analysis Results for the Cantilever Beam of L/h=10000 Modeled using Four Elements of Unequal Full: full integration; SRI: selective-reduced integration; SRI: selective-reduced integration,

Table 5 .
Normalized Mid-span Deflections of the Beam with Different Length-to-thickness, L/h, Ratios, Obtained using Different Beam Elements

Table 6 .
Normalized Free-end Deflections, Fixed-end Bending Moments, and Fixed-end Shear Forces of the Cantilever Beam for Different Number of Elements, Nel, Obtained using Different Beam Elements (a) Normalized Results using Linear Beam Elements