Application of Three Nodded Finite Element Beam Model to Beam on Elastic Foundation

The convergence of numerical solution based on two nodded beam finite element require considerable number of iterations and time; and is also plagued with shear locking. To address these deficiencies a three nodded beam element is proposed in this study to simulate the behavior of beams on elastic foundation. The analytical formulation of the model and development of shape functions are achieved with assumption of Winkler hypothesis for beam on elastic foundation A Matlab programme was developed to determine the combined beam and foundation stiffness as well as the load vector. The proposed model reliably simulates the deformations and stress resultants of beam on elastic foundation under general loading conditions. The result showed faster convergence devoid of shear locking. The maximum deflection and bending moment differ from the classical solution by about 5 percent.


Introduction
Recent environmental changes have resulted in frequent flooding of towns and cities, hurricanes, tsunamis, landslides and earthquakes etc around the world with attendant damages to infrastructure, causing economic damages as well as displacement of persons. One of the effects of such incidents is the distortion of the original stability of the foundation supporting structures. Most structures subject to these kinds of loading rely on raft and pile foundations for which the fundamental principle of analysis is based on the beam of elastic foundation. Beams on elastic foundation enjoy wide application in the field of railway engineering, harbor works, building frames and constructions, buried gas pipeline systems, machine foundation among others. Therefore, to reduce the economic impact of failures and collapses on investment, there is need for continuous improvement in solutions for the analysis and design of beams of elastic foundation.
The problem of beam on elastic foundation has been treated by many authors, including [1][2][3][4][5][6][7] and [8] among others. However, analytical solutions are appropriate only for a few cases of beams with simple geometry and loading. For other more realistic cases, the finite element method provides better solutions [9]. These numerous studies have shown that the stiffness of the beam and elasticity of the foundation depend largely on the modulus of the foundation soil which is modeled as a great number of springs such that the foundation reaction is directly proportional to local deflection of the beam.
A review of numerical methods for solution of beams and frames on elastic foundation subjected to static load using a linear element (two nodded element) is presented in [10]. Earlier, Hudson and Hutchinson had suggested that elastic foundation can be approximated by springs of equal stiffness, making it possible to account for system nonlinearities and to model and solve the practical soil-structure interaction problem [11] From the above, it is clear that the simulation of the behavior of the beam on elastic foundation critically depends on the accuracy of determination of the deflection of the foundation structure. From simple beam theory the elastic curve is at least a parabola of the second order, which will require at least three points to plot. Incidentally, numerical solutions have always been considered for a two nodded beam (first order linear element) having its shape function defined by a third degree polynomial in a bid to achieve convergence. It is clear that such an element can only produce a segmented line representation of the elastic curve and will therefore rely on a great number of elements for convergence. To improve on the performance of the two nodded element, higher order interpolation functions have been suggested. The application of these functions in terms of equal order polynomials resulted to shear locking [12]. Several attempts to solving this problem include the reduced integration Element (RIE), the Consistent Interpolation Element (CIE), and the Interdependent interpolation Element (IIE), all described in [13] and [14].
The two nodded beam finite element is ill conditioned as it cannot supply the boundary conditions required to determine the constants of integration. This means that a fully integrated linear element (two nodded beam element), with an overly inflexible set of shape functions is not able to capture the kinematics of deformations and reproduce the differential equations of the deformed shape of the beam on elastic foundation. This basic contradiction leads to shear locking as the beam is not able to recover its thin beam solution. Based on their extensive reviews, Prathap [15] and Eric [16] concluded that fully integrated first order solid linear element (2 nodded) such as Timoshenko beam elements, Mindlin plate elements may suffer from shear locking which therefore could give false results.
In order to resolve these contradictions, higher order beam elements of the third and fourth order have been proposed and shown to be locking free. However, these have not been extended to the study of beams on elastic foundation. Therefore, this paper is an attempt to address the problem of beam on elastic foundation by the stiffness matrix formulation of a three nodded beam finite element with a view to improving accuracy and resolving the shear lock problem. The classical Euler-Bernoulli beam theory is assumed and the deflection simulated with the quadratic beam element. By virtue of possession of additional degrees of freedom, the three nodded element as formulated in the present study is expected to more realistically simulate the deformations of the beam and therefore eliminate shear locking.

Method
The method applied is the finite element approach. It comprises the derivation of the beam and foundation stiffnesses as well as the force vector. The strain energy and virtual work methods are used in combination with Winkler foundation model. A computer programme in Matlab is developed to ease the matrix inversion and solution of the derived equations.
Consider the three nodded beam with nodal displacements and forces as shown in Figure 1.
Where 0 α to 5 α are the unknown constraints which define the degrees of freedom at the nodes 1, 2 and 3.
Using the boundary conditions of Figure 1, the nodal values can be represented as where N is the shape function, q is the nodal displacement and the y is the displacement field. The boundary values are The displacement can be represented in terms of the nodal displacement vector q in the form Nq y = where N is the shape function.

Element Mode Function
The beam local coordinates, node numbering and nodal forces and displacements are as shown in Figure 1. In order to formulate the shape functions, the displacement is represented in the form of a quintic polynomial as follows:

Beam Element Stiffness K b
For a conservative homogeneous system of constant stiffness EI, subjected to a bending moment M, the strain energy U is given by the expression Expressing M in terms of displacement y, equation 1 becomes Thus, the beam element stiffness k b is given by

Foundation Element Stiffness K f
Using the expression from Winkler as reported in [17], the strain energy U of the elastic foundation in the domain (0, 2L) is Where k o is the stiffness of the spring and b is the width of the beam under a vertical displacement y. Substituting

Combined Beam Elastic Foundation Stiffness K c
The total stiffness K c of the beam on elastic foundation is taken as the sum of the beam stiffness K b and foundation stiffness K f . In order to ease the calculations, a computer programme was set up in Matlab for the inversion of the matrices and computation of the combined stiffness. The Mat lab routine is shown in Appendix 1. The programme yielded the combined stiffness matrix shown in Appendix 2.

Nodal Load Vector
The nodal load vector is determined by the virtual work principle using the Lagrangian interpolation.
For number of nodes, n = 3, f (x) is the quadratic polynomial that passes through the three data points.
Considering a three node element with six bending degrees of freedom and a total length of 2L. This element has nodes at x= (0, ℓ , ℓ 2 ). The axial degree of freedom is omitted from this element. Shape function for this degree of freedom is the Lagrangian polynomial of order 2.  (iii) For the moment, a direct application of the value is imputed as decomposition is not required.

Equilibrium Equation
The stiffness matrices and the load P are assembled in the equilibrium equation Where d is the nodal displacement The solution was achieved through Matlab program composed specially for this purpose.

Result
Then results of the computations using the proposed model were tested for convergence and validated with numerical examples from other authors.

The Assessment of the Convergence
The assessment of convergence of solutions with various beam models was carried out for a beam supporting uniformly distributed load across span with end shear load.  This example sought to determine the maximum end deflection of a beam on elastic foundation supporting uniformly distributed load of q=200 kN/m across the span and an end shear load P=1000 kN. The solution of the problem was repeated in this present research using the proposed three nodded finite beam model. The results of the computations of the maximum deflection according to the various models as a function of the number of finite elements used are presented in Table 1 and plotted in Figure. 3 These results show that the value of end deflection obtained for the two nodded Euler-Bernoulli beam element is several orders of magnitude less than the theoretical value of 0.01895147m. This shows that the greater discretization did not reflect in any significant changes in deflection even as the number of elements was increased. This cannot be explained by the Timoshenko beam theory and expressly demonstrates the phenomenon of shear locking characteristic of the linear element (two nodded element). Several authors, including [19] have observed that this phenomenon is caused by an induced shear/strain energy constraint which dominates over the bending behavior of the beam.
The performance of the proposed model was assessed on the basis of the deflection profile and speed of convergence. The response of the three nodded model and that for the 2nodded Timoshenko shear deformable beam as shown in Table 2 and Figure 3 from which it can be seen that convergence is attain faster with the 3-nodded beam element model. The Results further show that the maximum deflection by the Euler-Bernoulli and Timoshenko beam appears to diverge as the number of elements is increased beyond 8 where as the result from the proposed three nodded element model rapidly converged to the classical solution with just 8 element.

Numerical Validation of the Model
As a further check on the adequacy of the 3-nodded finite element to reproduce the structural behavior of beam on elastic under general loading and boundary conditions, a numerical problem abstracted from Parvanova [18] was considered.
The beam was divided into 20 elements, 41 nodes and 82 degrees of freedom. The internal displacements obtained in both scenarios are shown in Appendix 3 and Figure 5. These were compared with the exact result obtained in [20].
The results show that the deflections were all within expected values on comparison with the three nodded beam element model as implemented using Math Lab codes and Midas Civil analysis software.     The model was seen to exhibit an elastic response under combined load action tracing the elastic line as computed by Parvanova [20] with recoverable deformations.
The three nodded beam element exhibits the true deformation response and is in agreement with the similar model implemented in Midas Civil software.

Discussion
The beam is modeled using 20 elements and implemented in Matlab environment. The results of deflection, rotation, shear force and bending moment are plotted with results of Parvanova [20] as shown in Appendix 3 and Figure 5 to  [20] and the present study are 5.263%, 7.183%, 11.51% and 5.58% respectively.
Also as evident from Appendix 4 and Figure 3, the deformation profile of the linear element with shear modulus changes repeatedly from linear to curved configuration and back to linear in a double curvature manner while the three nodded beam element model maintained a deformation pattern devoid of spurious energy levels as the element discretization was increased.
The convergence to exact value of maximum deflection of 0.018951m was achieved by using 8 finite 3-nodded elements while convergence was not achieved even with 16 finite elements for the 2-nodded Euler-Bernoulli beam. This confirms the ability of the proposed 3-nodded beam element to attain fast convergence with associated economy in computing time and resources. The lack of convergence for the Euler Bernoulli beam is suggestive of excessive rigidity in bending resulting in additional stiffening effect described as shear locking. These results of analysis with the proposed three nodded beam element further confirm that the shear lock problem can be eliminated by appropriately increasing the degree of approximation for finite element formulation.

Conclusion
A very important and challenging issue with beam on elastic foundation is the effect of soil-structure interaction on its structural behavior. The modeling concept of the soil as a bed of springs is often used while the beam is modeled as a two nodded beam element with or without a shear stiffness, the latter being intended for shear correction against shear locking. This research recommends the possibility of using a three nodded beam element to obtain locking free behavior of beam on elastic foundation.
The three nodded Finite element model provides better formulation of the beam finite element because it meets the requirement of a minimum of three points to generate the deflection curve given by the simple beam bending theory.
Also by virtue of possessing additional boundary conditions, the three nodded element permits a more realistic application of the quintic polynomial approximation thereby improving convergence, reducing computational time and eliminating shear lock.
The solution from the three-nodded beam model converges about three times faster to the exact solution for beam on elastic foundation than the fully integrated two nodded linear element in common practice.