COMPUTATIONAL ASPECTS OF Configuration dependent interpolation in non–linear HIGHER-ORDER 2D beam finite elementS
Edita Papa, Gordan Jelenić
Faculty of Civil Engineering, University of Rijeka,
Radmile Matejčić 3, 51 000 Rijeka, Croatia, ,
1.Introduction
In non-linear 3D beam theory with rotational degrees of freedom (Simo, 1985)configuration-dependent interpolation may be utilized to provide a result invariant to the choice of the beam reference axis (Borri and Bottasso, 1994)or invariant to a rigid-body rotation (Crisfield and Jelenić, 1999). For 2D beam elements, the latter issue vanishes, and such elements aremore illustrative for the study ofaccuracy of the configuration-dependent interpolation in higher-order elements.
Since the approximate character of the finite-element method stems from the introduction of interpolation functions for the field variables and their variations, the actual choice of the interpolation functions is of a great importance for the accuracy of the finite element method. The most commonly used standard Lagrangian interpolation procedure (with reduced integration used to eliminate the shear locking effect) is satisfactory if our attention is limited to finding the results for the field variables at nodal pointsin linear analysis. However, if we want to obtain the exact field distribution this kind of interpolation is unable to do so, even in linear analysis. To obtain the exact solutions in linear analysis simply using the right interpolation function and without applying reduced integration, a different kind of interpolation, called the linked interpolation is needed, in which the unknown displacement function do not depend only on the nodal displacements, but also on the nodal rotations (Jelenić and Papa, 2011). Marco Borri and Carlo Bottasso have developed a so-called fixed-pole formulation which effectively generalizes the idea of linked interpolation to non-linear two-noded beam elements. This interpolation is called the helicoidal interpolation by the authors andis based on the fact that the tangent to the beam centroidal axis and the normal of the cross section follow the same transformation rule(Borri and Bottasso, 1994). This kind of interpolation can also be obtained starting from the condition of constant Reissner’s strain measures(Reissner, 1972).
2.Reissner’s beam theory
In this paper we analyze geometrically exact Reissner’s beam theory (Reissner, 1972). In this beam theory, normal strain, shear strain and the change of the sectional orientation along the length of the beam are given by the following expressions:
/ (1)/ (2)
/ (3)
where is the slope of the beam in the initial configuration, is the function describing the slope of the cross sections in the deformed configuration as a function of the arc-length coordinate , and are the horizontal and vertical displacements that can be written in terms of the components of the initial and the deformed position vectors , , and , respectively, and a dash denotes a differentiation with respect to the arc-length coordinate (see Figure 1).
Figure 1Kinematics of the problem.
Equilibrium of a beam of length L follows from minimization of the difference between the strain energy and the work of the applied loading, , where
/ (4)/ (5)
and N,T andM are the normal stress, shear stress and stress-couple resultants, and are distributed horizontal, vertical and moment loading respectively, while p is a vector of nodal displacements and rotations and S is a vector of corresponding nodal loading. Taking (1)-(3) and writing minimization of the energy in a matrix form becomes
/ (6)where is the rotation matrix of a cross-section given by
/ (7)After introducing interpolation , where is a vector of the unknown displacement and rotation functions and is the interpolation matrix, expression (6) becomes
/ (8)where is a matrix of differential operators given by
/ (9)The first term within brackets in (8) is the internal force vector, while the term within the parentheses is the external force vector. Depending on the interpolationthe results will be more or less accurate.
Standard finite element procedure employs Lagrangian polynomials as interpolation functions (Simo and Vu-Quoc, 1986). This interpolation implies independence between the displacements and the rotation of the cross section which is physically incorrect. Also, this kind of independent interpolation results in ‘shear locking’ effect for certain aspect ratios when full integration is used. Additionally, the choice of the beam reference axis is an open question for the cases of more complex non-homogeneous or anisotropic beams, or beams with twisted or curved geometry. In these situations the axis of geometric, mass or shear centers of the cross-section do not necessarily coincide and choosing one of them for the beam reference axis instead of another may lead to different finite-element solutions. To overcome these problems, it makes sense to re-consider the choice of the interpolation functions that provide the actual interdependence between the displacement and the rotation fields. In linear analysis, such interpolation is called the linked interpolation, while in the non-linear analysis we call it the configuration-dependent interpolation.
3.Summary of the configuration-dependent interpolation (Borri and Bottasso, 1994)
As stated above, interpolation that provides the interdependence between the displacement and the rotation fields in non-linear analysis is called the configuration-dependent interpolation. In this section a brief overview of the known studies will be given.
Borri and Bottaso (Borri and Bottasso, 1994) introduce the so-called fixed pole interpolation in which the equilibrium of cross-sectional moments is stated with respect to a point which the authors call the fixed pole. This point has all the properties of the origin of an inertial reference frame. The fixed-pole strain measures are obtained by statically reducing all stress and stress-couple resultants to the fixed pole. Assuming a constant curvature in the deformed state and considering that the tangent to the beam centroidal axis and the normal of the cross-section follow the same transformation rule, for a two-noded element the following interpolation is obtained:
/ (10)/ (11)
/ (12)
where while are the Lagrangian polynomials for a two noded element:
/ (13)The same interpolationas in (10)-(12) can be obtained by assumingconstant strain measures in (1)-(3).This interpolation is non-linear in the field variables (configuration-dependent) and has the following limiting value in the case of the non-linear analysis becoming linear:
/ (14)This limiting case is the same as the linked interpolationin (Jelenić and Papa, 2011) known to be able of reproducing the exact solution of a linear problem thus eliminating the shear-locking completely. The latter, however, may be written for an element with arbitrary number of nodes:
/ (15)4.Configuration-dependent higher-order interpolation
In order to generalize the non-linear Borri-Bottasso result for two-noded elements to higher-order elements we may take the strain measures to be linearly distributed rather than constant, which however complicates the procedure considerably. This is so because the interpolation turns out to be related with Fresnel’s integrals, for which analytical solutions are not available (Jelenić and Papa, 2009). Alternatively, we may attempt to generalize the linear higher-order result obtained by Jelenić and Papa to non-linear analysis.The simplest idea of how to perform this generalization and come up with a configuration-dependent interpolation is to substitute the initial positions of the nodal points in the linked interpolation with their current positions:
/ (16)As it can be seen, the unknown position vectors in the deformed state depend on the rotations thus making this interpolation non-linear in the unknown functions, i.e. configuration-dependent.
5.Numerical examples
In this section two numerical examples are given. Both of them are calculated using the standard Lagrangian interpolation andthe linked interpolation in non-linear analysis for reduced integration as well as the full integration.
5.1Mattiasson cantilever beam
Figure 2Mattiasson cantilever beam problem.
Mattiasson cantilever beam (Mattiasson, 1981) with a transversal point load (Figure 2) is a large deflection problem most commonly used in testing the behavior of geometrically nonlinear beam finite element analysis. The value of Young’s modulus is E=20000, the area of the cross-section A=0.02, the shear area of the cross-section As=A/1.2, the inertia of the cross-sectionI=0.00025, the length of the beam L=1,shear modulus is given with respect to Young’s modulus G=10000 E and theapplied loading P=1.The results are given in Tables 1 and 2 for the reduced and the full integration. Ten linear elements were used to solve this problem and the results are compared with those from (Mattiasson, 1981).
Standard interpolation / Linked interpolation / Results from (Mattiasson, 1981)u / 0,0024662 / 0,0024662 / 0,00265
w / 0,0662269 / 0,0662269 / 0,06636
ϕ / 0,099659 / 0,099659 / 0,09964
Table1Reduced integration solutions.
Standard interpolation / Linked interpolation / Results from (Mattiasson, 1981)u / 8,04E-9 / 0,0025469 / 0,00265
w / 0,00011978 / 0,0659063 / 0,06636
ϕ / 0,00017967 / 0,0990763 / 0,09964
Table2Full integration solutions.
As it can be seen from the results shown in the tables above, reduced integration gives the same results for both the standard (Lagrangian) and the linked interpolation. However, looking at the results of the full integration it can be seen that the linked interpolation gives much better results than the standard interpolation. Therefore we can conclude that the linked interpolation is free from the locking effecteven when it is used in non-linear analysis.
5.2Hinged right-angle frame under fixed point load
Figure 3Hinged right-angle frame under point load.
This problem and its results aregiven in (Simo and Vu-Quoc, 1986). The data are: length L=120, inertia I=2, area of the cross section A=6, Young’s modulus E=7.2E06, Poisson’s ratio υ=0.3 and the vertical point load P=15000. This problem has been solved using ten linear and ten quadratic elements (five per leg) and the results are given in Tables 3 and 4. Results in (Simo and Vu-Quoc, 1986) are given for ten quadratic elements.
Standard interpolation / Linked interpolationu / 6,46073 / 6,46073
w / 22,4863 / 22,4863
Table 3Reduced integration solutions for the mesh consisting of ten linear elements.
Standard interpolation / Linked interpolationu / 0,0032633 / 0,421699
w / 0,228125 / 4,56821
Table4Full integration solutions for the mesh consisting of ten linear elements.
Taking quadratic elements to solve the same problem increases the accuracy of the solution. Those results are shown in Tables 5 and 6.
Standard interpolation / Linked interpolation / Results from the literatureu / 8,01638 / 8,01638 / 8,235
w / 25,8625 / 25,8625 / 25,8823
Table 5Reduced integration solutions for the mesh of ten quadratic elements.
Standard interpolation / Linked interpolationu / 3,34 / 5,42416
w / 14,3643 / 19,4069
Table 6Full integration solutions for the mesh of ten quadratic elements.
6.Concluding remarks
In this paper the theory of higher-order configuration-dependent interpolation (i.e. nonlinear in the field variables)has been introduced. Linked interpolation that gives exact solutions in linear analysis has been used to calculate the numerical examples in nonlinear analysis. It is shown that, in the case of reduced integration the results are identical to the ones obtained by using the standard Lagrangian interpolation. However, when using full integration the linked interpolation is significantly less sensitive to shear locking.
7.Acknowledgements
The results shown here were obtained within the scientific project No 114-0000000-3025: „Improved accuracy in non-linear beam elements withfinite 3D rotations” financially supported by the Ministry of Science,Education and Sports of the Republic of Croatia.
8.References
J.C. Simo. A finite strain beam formulation. The three-dimensional dynamic problem, Part I. Computer Methods in Applied Mechanics and Engineering, 49: 55-70, 1985.
M. Borri, C. Bottasso. An intrinsic beam model based on a helicoidal approximation --- Part I: Formulation. International Journal for Numerical Methods in Engineering, 37: 2267-2289, 1994.
M.A. Crisfield, G. Jelenić. Objectivity of strain measures in the geometrically exact three-dimensional beam theory and its finite-element implementation. Proceedings of the Royal Society London Series A, 455: 1125-1147, 1999.
G. Jelenić, E. Papa. Exact solution of 3D Timoshenko beam problem using linked interpolation of arbitrary order. Arch Appl Mech(2011) 81: 171-183.
E. Reissner. On one-dimensional large-displacement finite-strain beam theory. Journal of App. Math. And Phys (ZAMP). 23: 795-804, 1972.
G. Jelenić, E. Papa. Configuration-dependent interpolation in non-linear 2D beam finite elements. Proceedings of the 6th International Congress of Croatian Society of Mechanics, 2009.
K. Mattiasson. Numerical results from large deflection beam and frame problems analysed by means of elliptic integrals. Int.J.Num.Meth.Engng., 16, 145-153, 1981.
J.C. Simo, L. Vu-Quoc. A three-dimensional finite-strain rod model. Part II: Computational aspects. Computer methods in applied mechanics and engineering, 58 (1986) 79-116
J.C. Simo, L. Vu-Quoc. On the dynamics of flexible beams under large overall motions – the plane case: Part I and Part II. ASME Journal of Applied Mechanics, 53: 849-854, 1986.