Two nodes cusp geometry beam element by using condensed IGA

A cusp is a curve which is made by projecting a smooth curve in the 3D Euclidean space on a plane. Such a projection results in a curve whose singularities are self-crossing points or ordinary cusps. Self-crossing points created when two different points of the curves have the same projection at a point. Ordinary cusps created when the tangent to the curve is parallel to the direction of projection on a single point. The study of a cusp geometry beam is more complex than that of a straight beam because the structural deformations of the cusp geometry beam depend also on the coupled tangential displacement caused by the singular geometry. The Isogeometric Approach (IGA) is a computational geometry based on a series of polynomial basis functions used to represent the exact geometry. In IGA, the cusp geometry of the beam element can be modeled exactly. A thick cusp geometry beam element can be developed based on the Timoshenko beam theory, which allows the vertical shear deformation and rotatory inertia effects. The shape of the beam geometry and the shape functions formulation of the element can be obtained from IGA. However, in IGA, the number of equations will increase according to the number of degree of freedom (DOF) at the control points. A new condensation method is adopted to reduce the number of equations at the control points so that it becomes a standard two-node 6-DOF beam element. This paper highlights the application of IGA of a cusp geometry Timoshenko beam element in the context of finite element analysis and proposes a new condensation method to eliminate the drawbacks elevated by the conventional IGA. Examples are given to verify the effectiveness of the condensation method in static and free vibration problems.


Beam for machinery part
Most of the machinery parts in human-made engineering machines are composed of beams. A beam is a longslender-like structural member whose primary function is to support vertical loading and spread it to the supports ( Fig. 1.). The beam resists the vertical loading mainly through bending moment. A bending mechanism yields compressive zone on bottom side, and tensile zone on the top side of a beam section which depends on the loading arrangements. In a building, a column is also designed based on the beam theory with a consideration that it can undergo a buckling phenomenon if loaded by an extreme compression loading. A beam can hold a combination of loading types. When a beam is hinged at both of its two ends where they are allowed to rotate, the beam is called a simply supported. If the beam is fixed at one end and freed at the other end, the beam is called a cantilever beam.

Beam element formulation
In the finite element formulation, a beam is modeled by a straight line to represent the axis of the beam in the longitudinal direction, or perpendicular to its cross section (Fig. 2.). A typical beam element has two nodal points or nodes at both ends of the line. Each node is assigned movements to displace or rotate thus it is called a Degree-Of-Freedom (DOF). The term degree is showing some movements to displace or to rotate the nodes. The type of movements at a node can either be a combination between displacement and rotation. In the numerical computation of the standard method of the formulations of beam element in the Finite Element Method (FEM) are expressed by the matrices and vectors notations. The simplest and general plane beam element has two nodes at both ends, and each node has a horizontal (tangential) displacement, a vertical (transversal) displacement and a counter-clockwise (CCW) rotation as the DOF. Therefore, there are 3-DOF at each node. Thus 6-DOF is  attributed to a beam element with two nodes. In the vector notation, this 6-DOF is called as a generalized displacement which is given as follow.

Kinematic assumptions
As depicted in Fig. 3., the kinematic assumption of a beam called the Timoshenko beam has governing differential equations in which the shear deformation of the crosssection is taken into consideration. Assuming the small displacement, the following approximations can be given as: Fig. 3c., when a beam deformed the displacements at any point P are assumed by the following equations.
where uo and uP are the horizontal displacements of the point O and P; vo and vP are the vertical displacements of the point O and P; s and r are the tangential and radial coordinates of point P in initial configuration; θ is the rotation of the cross-section of the beam which is equivalent to slope of the beam at point O; φ is the angular location of point O; γ is the shear deformation of the beam cross-section. The boundary conditions between two elements divisions at the node O requiring that the vertical displacement v and rotational angle θ have to be continuous (Fig. 3d.).

Governing equation
In the theory of elasticity, the strain components can be obtained from the compatibility equation as given in the following equations.  where G is the shear modulus of the beam. Because, in the Timoshenko beam element the shear deformation is taken into account, the strain consists of the axial and shear strains components.
The governing equations of the beam element are derived based on the Hamilton's principle. Hamilton's principle is a generality of the virtual work energy principle to obtain the equilibrium equations in the dynamic system. Therefore, the formulations of governing equations of the beam element can be obtained for dynamic and static systems.
The equation of motion can be derived from, The kinetic energy, KE, of the beam element with the length L is given as The external work, WE, of the beam element with the length L is given as The stiffness and mass matrices, and the loading vector can be expressed by where the generalized displacements are given by, the compatibility matrix which is constructed from the defined shape functions of the generalized displacements.

IGA for free curve beam element
In IGA, the real geometry and shape functions of a curved beam element can be created systematically by using the NURBS [1-6].

NURBS for cusp beam element
In IGA, the real coordinate of the shape of the beam can be represented by using the NURBS curves where the path of a curve is defined in local coordinate ξ as, The detail of derivatives of the shape functions can be found in [7].
The displacements and rotation of an arbitrary point along the beam can then be calculated from the It is important to be noted that in IGA, the displacements and rotations of a node are treated separately by using each rational basis B-spline functions. Hence, the formulations are not coupled.
Depends on the order of the p th degree polynomial and knot vector (Ξ) being used, the preparation of the corresponding number of control points (CP) is necessary. Increasing of the control points to construct the NURBS curve will also increase the number of DOF of the beam as shown in Fig. 4.
The figure depicts the illustration of DOF at the control points in the beam with the geometry of the p th order of the polynomial.

Condensation method
The conventional static condensation and dynamic condensation is known to be effective to reduce the DOF of a beam which has more than two nodes. The number DOF of the mid-nodes are "hidden" so that only the DOF at both ends of the beam still reserved [8][9][10].
Although the importance of condensation is not significant in beam problems due to the number of DOF are not large, but by condensing numerous DOF [11][12] to the standard 2-node-6-DOF beam, the which has many benefits for real practices. First, the introduction of NURBS into the design practice will hinder some design practitioners that are not familiar with the new concept. Second, the condensed NURBS element codes can be integrated seamlessly into the existing beam finite element codes.

Proposed Condensation method
The dynamic equilibrium of three simultaneous equations can be represented in the matrix form as, 11 12 The element matrices are break-down to 3 sub-matrix simultaneous equations. By separating the matrix into 3 simultaneous sub-matrix equations, a second order relationship between the general displacements can be obtained, because there are 3 unknowns. The conventional dynamic condensation is modified by taking two times derivation to the free vibration equilibrium equations regards to the time t. Then, the relationship between the stiffness matrices k and the acceleration vector d can be derived.

Modified dynamic condensation
After some substitutions, a condensed modal equation can be used to solve any free vibration problems with only 2 nodes beam element. The detail of the derivations can be found in [1].

Convergence check
To demonstrate the efficiency of the proposed reduction method, a cantilever parabolic curved beam with the clamped boundary condition, its material properties, and static loadings cases, shown in Fig. 5., is investigated on the number of element required to achieve convergence results. Figure 6 shows the control points locations of the beam which depends on the number of element division used in the analysis of the beam with the 2 nd order polynomial parabolic function.         The condensed IGA is applied to the cusp curved beam element by using 3 element divisions which verified in the previous section. Figures 13 and 16 show the control points when the beam is divided into 2, 4 and 6 elements, respectively. Figures 14 and 17 illustrate the vibration modes up to the first 3 lowest natural frequency of the cantilever cusp curved beam.
Four Gauss' integration points are adopted to integrate the stiffness and mass matrices of the cusp curved beams. The degree of polynomial used in the NURBS functions are pu = 3, pv =3 and pθ = 3. The beam is divided into 3 elements for each section of the curve, thus 6 elements in total element divisions. The element stiffness and matrices are condensed to a two-node 6-DOF beam element. Hence after the condensation and assembling process, there will be three nodes in the beam.
The knot vector used in each element division is given as For static analyses, three external loading conditions of horizontal force, vertical force and bending moment are applied at the free end tip of the cantilever beams.        Table 1 tabulated the results of static and modal vibration analyses of Example 1. Under the static loading cases, the displacements and rotations computed at the tip end of the cantilever beam analyzed by the condensed IGA show less than 0.43% maximum discrepancy when compared to the results of FEM by using 60 element divisions.
For the free vibration analysis of Example 1, in the results of the first 3 modes comparisons, the maximum difference between the present condensed IGA and FEM with 60 elements are within the range of 0.19% to 0.66% which are negligibly small.
In Table 2, the results of static and modal vibration analyses of Example 2 are shown. Under the static loading cases, the displacements and rotations computed at the tip end of the cantilever beam analyzed by the condensed IGA show less than 0.21% maximum discrepancy when compared to the results of FEM by using 60 element divisions.
Similarly, the free vibration analysis of Example 2 resulted in the first 3 modes with the maximum difference between the present condensed IGA and FEM with 60 elements are within the range of 0.11% to 0.36% which are also negligibly very small quantities.

Summary and conclusions
A new condensed IGA method was proposed for general curved beams, including cusp curved beam element, with complex geometry such as the radius of curvature along the longitudinal axis of the beam. the effect of shear deformation and rotary inertia were taken into consideration in the formulations.
The Gauss-Legendre scheme was adopted to construct the stiffness, mass matrices, and loading vector by integrating along the curvilinear coordinate of the beam axis.
The superiority of the present condensed IGA method lies in the fact that they were derived based on the exact geometry of the beam by using NURBS instead of being assumed linear element.
The new condensed IGA proved to be accurate to compute the general displacements, and natural frequencies of the cusp curved beams by only using 3 elements are consistent for each arch.
Altogether, numerical examples for static and free vibration are within 1% of discrepancies which is very good in agreement with those obtained by using 60 elements in FEM.
The present proposed condensed IGA turned out to combine the preciseness of the IGA with the conventional two-node-6-DOF standard element in FEM, and could be considered as an advantageous step towards deriving more efficient element formulation from the theoretical point of view.