Development of Sliding Connections for Structural Analysis by a Total Lagrangian FEM Formulation

Abstract

In this study a total lagrangian 2D finite element formulation is used to model plane frames developing large displacements and rotations considering sliding connections. This kind of connections is usually called prismatic and cylindrical joints.

In order to be self-containing, the steps of the development of a frame finite element of any approximation order that considers the influence of shear strain by means of a generalized Reissner kinematics is presented. The adopted degrees of freedom are positions and rotations. Using positions as degrees of freedom simplifies the total lagrangian description and enables a comprehensive presentation of the proposed connections. Revolute connections are considered by direct degrees of freedom matching. Prismatic connections are modelled by the Lagrange multiplier technique that constrains positions and rotation of a sliding node to the varying position and rotation of a path element. Cylindrical joints are introduced in similar way by Lagrange multipliers releasing the sliding node rotation.

The principle of stationary potential energy is used to write the non-linear equilibrium equation including the Lagrange multiplier influence. To solve the non-linear equation a Taylor expansion is carried out and the Newton-Raphson procedure is employed. The frame element is considered elastic, following the Saint-Venant-Kirchhoff constitutive model. Selected examples are used to validate the formulation and to show its possibilities of application.

Keywords:
Positional Finite Element Method; Sliding connections; Lagrange multipliers

1 INTRODUCTION

Considering the constant evolution of materials science, the use of materials with improved mechanical properties is a reality in daily engineering applications. Therefore, the design of audacious structures each time slenderer and lighter is a constant challenge for engineering. From this point of view, the consideration of geometrical non-linearity in structural analysis is becoming an imposition nowadays. On the other hand, the analysis of deployable and bi-stable structures is also an important field of mechanical and aerospace engineering. Some applications as satellite antennas, deflectors, solar panels and solar sails are designed using the concept of sliding connections, large displacements and rotations.

As far as the authors knowledge goes, differently from what is proposed in this study, all existing finite element method (FEM) formulations related to this kind of connection that considers large displacements and rotations are based on the updated lagrangian concept, more specifically co-rotational formulations. Those formulations are presented in the works of Simo and Vu-Quoc (1986), Armero (2006), Jelenic and Crisfield (2001) and Ibrahimbegović and Taylor (2002). Moreover, the specific developments of mathematical models of such connections are in its majority related to multibody dynamics, becoming a hard path to researchers not trained on this subject(Bauchau, 1998; Cardona; Géradin; Doan, 1991; Géradin; Cardona, 2001; Sugiyama; Escalona; Shabana, 2003). This approach may create a dependence of dynamic analysis when some simple static modelling could be of interest to describe low velocity mechanical systems or sliding connections in statics. We could find one static work, Jelenic and Crisfield (1996), that takes care of static 3D connections, however the connections are written regarding nodal variables and limits the movement to straight directions not associated to the flexibility of structural members.

One can see, for example, in Bauchau (2000) and Bauchau and Bottasso (2001) the use of Lagrange multipliers to solve sliding joints in dynamic analysis, however, the application of constraints are done directly at the equilibrium level and the non-linear solution process is not detailed. Moreover, the authors decided to use non-dimensional variables as the main unknowns of the problem, limiting the sliding equations to one single element and introducing some difficulties on the interpretation of forces and mass values at the joint node. In the works of Bauchau (1998), Cardona, Géradin and Doan (1991), Bauchau (2000) and Bauchau and Bottasso (2001) the Hamilton principle is used to derive the elastodynamic equations for the geometrical non-linear FEM formulation dedicated to elastic multibody dynamics. Some comments about difficulties regarding using the Newmark method for time integration and the proposition of energy conserving algorithms are also described in literature on Bauchau (1998), Cardona, Géradin and Doan (1991), Romero and Armero (2002) and Leyendecker, Betsch and Steinmann (2006). Other interesting works related to multibody dynamic applications can be mentioned such as Greco and Coda(2006), Yoo et al. (2007), Hong and Ren (2011) and Moon et al. (2014).

In this work we are interested to improve a previously developed total lagrangian frame finite element formulation (Coda; Paccola, 2014; Reis; Coda, 2014) to accomplish sliding connections for 2D elastostatic applications. In order to enable a comprehensive description of the proposed connections the frame formulation is described. This formulation is based on position degrees of freedom and a comprehensive description of the structural kinematics is presented. Revolute connections are considered by direct degrees of freedom matching, prismatic connections are modelled by the Lagrange multiplier technique that constrain positions and rotation of a sliding node to the position and rotation at general and varying location of a path element. The combination of prismatic and revolute joints is done by Lagrange multipliers letting the rotation free and resulting in the cylindrical joint model. Moreover, the main variable at the connection of elements respects the Euclidean space dimension and consistent results are achieved.

The non-linear equilibrium equation is derived from the stationary total potential energy principle and solutions along the equilibrium path are achieved by means of the Newton-Raphson procedure. Selected examples are presented to validate the formulation and to show its possibilities of application.

2 PLANE FRAME KINEMATICS

The adopted FEM solution procedure is written here from the principle of stationary potential energy and, therefore, the strain energy stored in the frame element must be calculated. To calculate the strain energy, it is necessary to define the way one achieves the strain distribution as a function of the body position limited to a finite number of degrees of freedom. In order to do so we write the initial and current positions of frame elements as a function of nodal position, angles and non-dimensional variables.

2.1 Initial Configuration

To build the total lagrangian procedure we start describing the initial configuration of a frame element. Figure 1 shows the reference line of the initial configuration, the non-dimensional space from which the reference line mapping is written and the nodes of the finite element.

The mapping of the reference line is written as:

(1)

where, i is the direction (1 or 2), m indicates reference line and ℓ is the index that represents nodes and shape functions φ_ℓ (ξ) (Lagrange polynomials of any order). X^m _iℓ represents the nodal coordinates at the reference line and repeated indices indicate summation.

Figure 1:
Reference line parameterization for initial position (cubic approximation).

At the initial configuration, the nodal coordinates are known and the cross sections are considered orthogonal to the reference line. Thus, from Figure 1, one calculates the tangent vectors at each node k by:

(2)

in which comma indicates derivative.

From equation (2) the normal vectors that generate cross sections at nodes are written as:

(3)

where the index m has been omitted. To use angles as nodal parameters, one calculates, from Figure 1:

(4)

To define a cross section at any position along the frame element it is necessary to know its angles. In order to do so we use the same approximation adopted for coordinates, i.e.:

(5)

From the initial position of reference line and cross sections orientation one defines the initial configuration of the frame element as depicted in Figure 2.

Figure 2:
Initial configuration mapping.

One can deduce from Figure 2 that the initial configuration mapping is written using a non-dimensional space (ξ,η) by:

(6)

where the vector g ⁰ _i (ξ,η) can be written as a function of θ ⁰(ξ), h and the height h ₀ of the finite element as:

(7)

Substituting equations (1) and (7) into equation (6) results the initial configuration mapping of the frame element, as:

(8)

(9)

It is worth noting that, for the sake of simplicity, the reference line is adopted at the centre of the cross section and the width and height of the bar are considered constant along all frame element. One can adopt the reference line at any position and the internal force offsets would appear naturally(Coda; Sampaio; Paccola, 2015; Pascon; Coda, 2013).

2.2 Current Configuration

The current configuration is understood as a trial position of the solution process and can be seen in Figure 3.

Figure 3:
Current configuration mapping.

Therefore, the new nodal coordinates Y^m _iℓ and angles θ_ℓ can be used to derive the current mapping from the non-dimensional space (ξ,η) to the current configuration as:

(10)

(11)

where f _{1 i} (ξ,η) = y_i ((ξ,η) is the new mapping understood as the current positions of any point inside the body at the current configuration mapped from the non-dimensional space.

As, at current configuration, there is no relation between θ_ℓ and the reference line inclination, see Figure 3, the cross section is not orthogonal to the reference line, consequently, Reissner kinematics takes place. It is worth noting that the height of the frame element is not free to change and, to avoid volumetric locking the constitutive equation should be relaxed in order to exclude transverse expansions.

2.3 Complete Mapping and Green Strain

After defining the mappings from the non-dimensional space to the initial B ₀ and current B configurations, we introduce the mapping from the initial configuration to the current configuration in the same sense presented by Bonet et al. (2000) and Coda (2003). In Figure 4 we put together Figures 2 and 3 revealing the mapping as:

(12)

where equations (8) to (11). are defined by and

Figure 4:
Change of configuration function or deformation function.

To calculate strain and strain energy it is not necessary to explicitly know Bonet et al., 2000; Coda, 2003) as:, but its gradient, called here A. This is done in a very simple way (

(13)

in which

(14)

For a trial current nodal position, both A⁰ and A¹ are numerical values at a generic point (ξ,η) In the numerical strategy, coordinates (ξ,η) are Gauss integration points resulting in a very simple numerical procedure. Moreover, expression (14) indicates that there is no difference modelling straight or curved elements by the proposed formulation.

Substituting equations (8), (9), (10) and (11) into equation (14) results:

(15)

(16)

The Green-Lagrange strain E is adopted to develop the formulation, as it is an objective measure (Ogden, 1984), and is given by:

(17)

In which I is the identity tensor of second order and C is the right Cauchy-Green stretch.

3 EQUILIBRIUM EQUATION

As mentioned in the introduction, the equilibrium equation is achieved here from the principle of stationary energy. For the considered isothermal, static and elastic situation the total energy is written as:

(18)

where Π is the total energy of the system, U_e is the stored elastic energy and P is the potential of conservative external forces. Therefore, the variational principle is stated by:

(19)

in which the symbol δ means variation.

3.1 Strain Energy Adopted to the Proposed Frame Element

In order to calculate the strain energy of the analysed body (frame element) it is necessary to integrate over the initial volume (lagrangian formulation) the specific strain energy written here as a function of the Green-Lagrange strain (E) as:

(20)

where 𝔼 is the longitudinal elastic parameter of the Saint-Venant-Kirchhoff constitutive model that approaches the Young modulus for small strain. The shear elastic modulus is 𝔾 = 𝔼/[2(1 + ν)] being ν a constant that reproduces the Poisson ratio for small strain. As one can see, in the present formulation the Poisson ratio interferes only on the shear stiffness and does not introduces transverse expansion. The energy conjugacy property allows one to define the second Piola-Kirchhoff stress as:

(21)

moreover, the constitutive elastic tensor ℭ can be achieved as:

(22)

It is worth noting that the adopted constitutive model avoids volumetric locking and that the high order of finite elements avoids shear locking (Bischoff; Ramm, 2000).

As the Green strain has been written as a function of nodal positions, see equations (15), (16) and (17) the strain energy accumulated in the structure is also a function of positions, and is written as:

(23)

in which V ₀ is the initial volume of the frame.

3.2 Total Potential Energy and Equilibrium

The total potential energy of the mechanical system in an elastic, static and isothermal condition is given by the sum of the strain energy and the potential energy of external forces, as:

(24)

In this expression q_i = ϕ_ℓ (ξ)Q_ℓi , Equations (10) and (11)) and dS ₀ is the initial infinitesimal reference line length of a curved frame element. is the vector of external distributed forces written as a function of nodal values as is the conservative applied forces (including moments), is the current position of the reference line (

To achieve the equilibrium equation one applies the principle of stationary potential energy on expression (24) using nodal positions as parameters,

(25)

Knowing that the strain energy is a function of Green strain which is a function of positions, and that the applied forces are conservative, equation (25) is rewritten as:

(26)

Knowing the current reference line expression equation (21), one writes equation (26) as:, the adopted surface force approximation and

(27)

in which the symbol ⊗ represents the tensor product.

It is usual to write equation (27) as (Coda; Paccola, 2014):

(28)

in which Coda; Paccola, 2014). As the variation equation (28) results in the geometrical non-linear equilibrium equations for discrete frame analysis, as: is arbitrary, represents nodal forces and L is the matrix that transforms distributed forces into nodal equivalent ones (

(29)

3.3 Piola-Kirchhoff Stress and Green Strain Derivatives

In order to be possible the reproduction of the proposed formulation it is important to show how the kernels of equation (27) are calculated. The second Piola-Kirchhoff stress S is obtained from equations (20) and (21) as:

(30)

The first derivative of the Green strain regarding positions is achieved from equations (15), (16) and (17) as:

(31)

As equation (31). From equation (15) it is possible to derive A¹ regarding nodal positions Y^β _α, for which α is the direction (1, 2 or 3, for the angle) and β is the element node, as follows: one can simplify

(32)

(33)

with summation over nodes ℓ and k.

From expressions (30) through (33) the internal force is calculated using Gaussian quadrature over the initial volume, as:

(34)

in which w_ig and w_jg are Gauss weights and ξ_ig and η_jg are integration points.

4 SOLUTION OF THE NON-LINEAR EQUILIBRIUM EQUATIONS

In order to organize the solution procedure, based on the Newton-Raphson method, one rewrites equation (29) coupling external forces in one unique vector , as:

(35)

In this expression is called the unbalanced mechanical force vector, as it is null only when the position is the solution of the problem.

Equation (35) can be understood simply by finding . This statement is clearly non-linear regarding . Therefore, to find a solution of the problem one expand the equality in a Taylor series truncated at the first order, such that

(36)

in which is a trial solution.

Solving the linear system of equation (36) one finds a correction Δ to the trial solution as:

(37)

where the gradient of is the Hessian of the total potential energy. As the external forces are conservative, the Hessian matrix is given by:

(38)

Using Δ, a new trial position is achieved, as:

(39)

This new trial solution is used again in equation (37) until:

(40)

when (35). In equation (40) is assumed to be the acceptable solution of is the Euclidian norm of the initial position vector.

4.1 Hessian of the Proposed Frame Element

To complete the description, the frame element Hessian matrix H is presented as:

(41)

In the last integral, the derivative of the Piola-Kirchhoff stress is given by:

(42)

for which the derivatives of the Green strain regarding positions is already given in equations (31) and (32).

Thus, the last term to be achieved is the second derivative of the Green strain regarding positions, i.e.:

(43)

with,

(44)

The non-null values of are restricted to the third direction (angles) and are given by:

(45)

where β and ζ are element nodes and summation over ℓ and k are present.

The integration over elements is given by Gaussian quadrature, as:

(46)

for which w_ig and w_jg are weights and ξ_ig and η_jg are integration points.

5 SLIDING CONNECTIONS

As mentioned in the introduction section, revolute connections are not sliding ones and are considered by the direct matching of degrees of freedom, readers are invited to see the references Greco and Coda (2006), Reis and Coda (2014) and Coda and Paccola (2014) for details.

A prismatic joint, schematically depicted in Figure 5a, constrains the extremity position of the sliding element to slide over a path element without allowing relative rotation. A cylindrical joint constrains the extremity of the sliding element to slide over a path element but releases the relative rotation between elements, see Figure 5b.

Figure 5:
(a) Prismatic joint, (b) cylindrical joint.

We describe in details the prismatic connection, for which Lagrange multipliers are used to constrain positions and relative angles. The cylindrical connection is briefly described indicating which terms of the prismatic joint are disregarded. One important advantage of using Lagrange multipliers is the simplicity of the technique when considering the principle of stationary total potential energy and the total lagrangian description.

It is important to mention that, differently from what we are proposing, the consulted works are based on multibody dynamic updated lagrangian formulations and, when considering flexible elements, two options of sliding connections are present. The first considers the direction of sliding defined by a fixed director and, therefore, does not attach elements. The second considers the free variable as the non-dimensional coordinate of the master element and not the real position. This last brings some difficulties on interpreting the real mass and forces to be considered at the connection nodes. See, for instance, the works of Cardona, Géradin and Doan (1991), Bauchau (1998), Bauchau (2000), Géradin and Cardona (2001), Bauchau and Bottasso (2001), Garcia-Vallejo et al. (2003), Sugiyama, Escalona and Shabana (2003) and Lee et al. (2008).

5.1 Kinematical Constraints by Lagrange Multiplier

A set of path elements defines the sliding path of the connected node of a sliding element, this connected node is now on called sliding node. In what follows is used for path elements and for sliding elements, including the sliding node.

The variable s_p = s(ξ_p ), depicted in Figure 5, defines the curvilinear position and the cross section orientation of the path point θ ⁰ _p , calculated at the initial configuration (Figure 6), should be constant during the sliding process. These constraints are described by: (with coordinates belonging to the sliding path. For prismatic connections, in addition to positions, the difference of cross sections orientations ∆) over the sliding path. For cylindrical connections the cartesian coordinates (Ŷ Pi) of the sliding node should be equal to the cartesian coordinates of point

(47)

and,

(48)

or, in a general, form:

(49)

where δ_{i 3} is the Kronecker delta.

Figure 6:
Difference of angles for a single node at different elements.

During sliding, the curvilinear position s_p = s(ξ_p ) or, inversely, the non-dimensional coordinate ξ_{p =} ξ(s_p) will vary.

To impose the kinematical constraint of a prismatic joint via Lagrange multiplier one adds to the total potential energy the following potential:

(50)

in which the term inside brackets is exactly the kinematical constraint given by equation (49), λ_i are the so called Lagrange multipliers, one for each direction i (summation implied), and represents the connected positions. In mechanics an interesting physical interpretation can be given for the Lagrange multiplier, that is, its converged value is the auto equilibrated force (action and reaction) necessary to keep both bars together.

The new potential energy is written as:

(51)

in which expression (51) to find equilibrium equations, as: as is already included in . The principle of stationary potential energy is applied on include all degrees of freedom, less the new ones

(52)

in which δΠ has already been detailed in previous sections and include variations of the strain energy and applied external forces regarding current nodal positions .

To complete the variation of : it is necessary to solve δ𝕃. This is done regarding and

(53)

being α the nodes of the active path element. One should note that variable ξ_p instead of s_p as the main variable. belongs to and will result in terms to be added in internal force and Hessian matrix at the solution process. At this, it is important to mention that the consulted works that deal with flexible connections use

For prismatic connections, developing the derivatives, one writes in an open form:

(54)

in which the nodes belonging to the sliding element, different from the sliding node, represented by index k are present in equation (54) in order to simplify the 'Lagrange Force' vector assemblage.

In equation (54) the following property is used,

(55)

in which,

(56)

The new equilibrium equation is now assembled respecting the degrees of freedom correspondence adding to , i.e.:

(57)

in which includes all external loads.

For cylindrical connections, the same procedure is used, remembering that the relative angle is free, resulting:

(58)

5.2 Non-Linear Solution

The non-linear set of equilibrium equations (57) is also solved by the Newton-Raphson procedure. The unbalance mechanical vector is written as:

(59)

For a trial solution (59) does not hold, therefore a Taylor expansion of first order is performed and the nullity is imposed: equality

(60)

From this equation, the trial vector is corrected by:

(61)

until a prescribed tolerance is respected, see section 4. This time, both the trial position and the Hessian matrix contain the usual contribution of connected and unconnected nodes and the new contribution of s_p and λ_i . The achievement of s_p is not sufficient to update ξ_p = ξ(s_p ) is not explicitly written. The solution of this stage is described in next item. and the Hessian matrix, as the function

The new Hessian matrix H_L = ^con _L, is achieved by the second derivative of s_p and λ_i . The part of the Hessian matrix that corresponds to the element connection degrees of freedom is achieved as: regarding positions and the new variables = H + H

(62)

or, written in a way that the corresponding variables are identified:

(63)

in which and . The terms of HconL are given as follows:

(64)

with,

(65)

in which index notation is applied. Index i corresponds to directions (1, 2 or 3, for prismatic connection), k also represents directions (1 and 2), and, ℓ, m and n refer to the nodes of the active path element. For cylindrical connection i assumes 1 or 2 only. The null terms in matrix H^con _L have already been calculated in matrix H and corresponds to the second derivative regarding nodal positions of the strain energy related to the connected (path and sliding) elements.

5.3 Curvilinear and Non-Dimensional Variables

The definition of s_p = s(ξ_p) that is directly determined by the update of variables, equation (60), at the Newton-Raphson solution procedure is of valuable importance for the consideration of friction. It is important to stress that this study does not consider friction, only prepares the solution to be in the real Euclidean space for future applications.

However, the calculation of the internal force equation (58), and the Hessian matrix, equation (63), associated to the Lagrange multiplier technique are dependent of ξp = ξ(sp) that is not explicitly defined. Therefore, it is necessary to calculate, for a given trial equilibrium position, the value of ξp. This is done solving the following non-linear system of equations:,

(66)

that is the mechanical constraint present in equation (50), obviously with known values of Ŷ_i ^P and , and represents, for both directions, the residue ri(ξp) when in the iteration process.

As the system is overdetermined we apply the least square technique to find the unknown (Nocedal; Wright, 1999). Therefore, the following objective function is defined:

(67)

This is expanded in a Taylor series of first order as:

(68)

being ξ⁰ _p a previously known trial value. When trying to solve equation (68) one realizes that, as p(ξ_p ) ≥ 0 one achieves Δ_p (ξ_p ) = 0 at solution, i.e, the objectivity of equation (68) is lost at the solution. Due to this, the least square method uses the assumption Δ_p (ξ_p ) = 0 to solve this equation. Therefore, from expression (68) one writes:

(69)

Results from the last expression the Newton method applied to an minimization problem when determining Δξ_p by:

(70)

The non-dimensional variable is updated by ξ_{p =} ξ ⁰ _p + Δξ_p until | Δξ_p/ξ ⁰ _p| respects a given tolerance. The terms in equation (70) are calculated by:

(71)

and

(72)

in which the indexes ℓ and m imply summation over the nodes from the active path element on both directions i.

With ξ_p known for a trial s_p , implicit in the values of Ŷ_i ^P and equations (60) and (61) continues as previously stated. Also, comparing the numerical value of the achieved non-dimensional variable to the dimensionless space interval, one decides the necessity, or not, of path elements transition., the global solution process described by

6 INTERNAL EFFORTS

Internal efforts are calculated directly from the integral of Cauchy stress σ over the cross section. As the natural stress measure calculated from the Green strain and Saint-Venant-Kirchhoff relation is the second Piola-Kirchhoff stress S, it is necessary to apply the well-known relation (Ogden, 1984):

(73)

in which A is the deformation gradient given by equation (13). Equation (73) gives the Cauchy stress written according to the initial reference system; then, it is necessary to rotate its components to the orientation of the analysed cross section.

As the cross section area A ₀ does not change its geometry, the internal efforts results:

(74)

in which ỹ ₁ and ỹ₂ along the cross section defined by the angle θ. is the rotated Cauchy stress for axis

7 EXAMPLES

In this section four examples are presented to demonstrate the precision and possibilities of the formulation. For all examples the stop criterion of the solution process was adopted as relative position with a tolerance established as 1.0.10^-8.

7.1 Driven Mechanism

An interesting use for sliding joints is the modelling of mechanisms capable of describing some type of envisioned geometry, as cutting rocks or metallic sheets for civil or mechanical industries. Such mechanisms are commonly described in classical textbooks such as Shigley and Uicker (1981) and Norton (2011), among others, where analytical solutions for position analysis are presented. Usually, the numerical modelling of this kind of mechanism is done in a dynamic version. However, if the intention is to adjust the geometry to be described, a static position analysis should be better suited, which is done in this example since no inertial forces are considered. Figure 7 depicts the initial configuration of a structure with a prismatic joint that connects a 6.0 m arm to a 1.0 m supporting bar. A prescribed rotation ψ is imposed on a crank that is connected to the main arm by a revolute joint. The rotation of the crank imposes a motion on the mechanism as depicted in Figure 8.

Figure 7:
Initial configuration and geometry.

Figure 8:
Selected configurations and arm free end trajectory.

Considering the restricted degrees of freedom and the boundary condition, this is an isostatic non-compliant linkage. Hence, only rigid boy motion takes place since nor strains or stresses are expected during the imposed quasi-static motion. Therefore, the dimensions and material properties of the involved bars may have any value. Although, in order to obtain a numerical solution, it is adopted for all elements square cross section with side 0.1 m, elastic modulus 𝔼 = 2.0Gpa and shear modulus 𝔾 = 1.0GPa. Seventeen finite elements with cubic approximation are employed for discretization.

Figure 9 shows the arm free end trajectory (cutting geometry), point A, for a complete turn of the crank for the simulation and the analytical rigid body motion solution given by the referred textbooks. Also, the prismatic joint position, point P, evolution according to the prescribed rotation is compared with the analytical solution as shown in Figure 10a and 10b for horizontal and vertical position, respectively.

Figure 9:
Arm free end trajectory (point A).

Figure 10:
Prismatic joint position analysis (point P): a) horizontal and b) vertical position.

The motion was divided into 100 steps and less than six iterations are required for convergence in each step, see Figure 11. Therefore, given the obtained results, one concludes that the presented formulation works very well and can be used for general analysis.

Figure 11:
Evolution of the number of iterations for solution.

7.2 Doubly Bent Beams with Bifurcation

Since in the last example finite deformations are present only due to rigid body motion, here we present an interesting example of bifurcation of equilibrium under traction that shows large displacements and rotations in a deformable structure involving sliding connections.

Figure 12 presents two flexible bars with length L = 0.25m that are initially aligned and to which is imposed a horizontal displacement u. The reactive traction force F in the right support is also measured. Those bars are connected by prismatic joints to a rigid bar long enough to allow large deformations on the structure. Here the rigid bar length is adopted as 10.0 m.

Figure 12:
Initial and deformed configurations of the structure.

Analytical solution and experimental validation of this phenomenon were presented by Zaccaria et. al. (2011). For the numerical simulation, a rectangular cross section is adopted to both flexible bars with h ₀ = 1.0mm and b ₀ = 25.0mm. For the rigid bar a squared 1.0 m side cross section is chose. Each of the three bars were modelled by four cubic finite elements with the same material properties: 𝔼 = 200.0GPa and 𝔾 = 100.0GPa.

In order to find the non-trivial solution, a perturbation is introduced to the system as an initial inclination of the rigid bar by an angle ϕ ₀. Increasingly smaller values for the initial angle were adopted to show convergence to the analytical solution. The results for the rigid bar rotation ϕ, Figure 13, and right support horizontal displacement, Figure 14, against a non-dimensional traction force 4FL ²/𝔼Iπ², where I is the second area inertia moment, shows the efficiency of the formulation when handling geometrical nonlinear problems.

Figure 13:
Rigid bar rotation angle.

Figure 14:
Right support horizontal displacement.

The horizontal displacement was imposed in 200 steps and the mean number of iterations required for solution was 13 for each step. Figure 15 presents the evolution of the number of iterations for the extreme cases of the initial inclination angle since the other results have a similar pattern. Only in the first step of the smallest initial angle several iterations were needed. This is expected since the strong deformation of the structure that occurs in this particular step is due to the bifurcation of equilibrium.

Figure 15:
Evolution of the number of iterations for solution.

7.3 Influence Lines of a Bridge - Moving Load

This example is presented in order to demonstrate the capabilities of the proposed sliding joint formulation to determine influence lines at any cross section for general structures. One important information about influence lines and internal efforts envelope for non-linear applications is that the superposition principle is not valid and exhaustive calculations are performed. However, when small displacements take place the superposition of results is recovered.

Particularly in this example, the influence lines of internal efforts at the central cross section (point P) of a bridge subjected to a load train is calculated, see Figure 16. The vehicle is modelled by a frame with 4.0 m length and 1.0 m high. The contact between the vehicle and the bridge is modelled by two cylindrical joints allowing their relative movement. A vertical 15 kN load is applied on the mid-span of the moving frame and a displacement u = 26.0m (divided into 500 steps) is imposed at its left corner. Figure 16 depicts the geometry of the bridge and the vehicle for the initial configuration.

Figure 16:
Vehicle initial configuration and geometry.

A rectangular cross section with width b ₀ = 1.0m and height h ₀ = 2.0m is adopted for the bridge, and a square cross section with side a = 0.1m is adopted for the vehicle. The bridge elastic and shear modulus adopted are 𝔼 = 20.0GPa and 𝔾 = 10.0GPa, respectively. For the vehicle, the material properties are 10 times greater. Six finite elements with cubic approximation were used to model the frame representing the vehicle and 34 elements for the bridge. Since the objective is to evaluate internal efforts, two finite elements with 1.0 mm each are placed on both sides of point P to avoid the passage of any concentrated load on its domain. This is necessary due to the discontinuity of internal efforts and the continuity of shape functions, inherent to finite elements.

Figure 17 presents the vertical displacement influence line at the bridge mid-span (point P). As expected, the greatest displacement of -0.55 mm occurs when the load train is centralized with the bridge beam. Influence lines for the bending moment and shear force are shown in Figures 18 and 19, respectively. As expected, the maximum efforts occur when one of the cylindrical joints (wheel) are placed over point P.

Figure 17:
Bridge mid-span displacement influence line.

Figure 18:
Bridge mid-span bending moment influence line.

Figure 19:
Bridge mid-span shear force influence line.

In the solution process every step needed exactly 10 iterations for convergence. As one can see, results are quite accurate and no limits for the bridge vertical geometry or number of wheels are present for the proposed formulation, revealing the versatility of the technique for practical applications.

7.4 Equilibrium Path of a Shallow Arch with Crank

A shallow arch with its motion enforced by a crank is present in this example. The arch is fixed by simple supports that allow rotation but not translations. The arc has a span of L = 10.0m and height of h = 1.0m. The crank is subjected to a rotation of ψ = 1.8 rad divided into 500 equal steps. The connection between the crank and the arch is made by means of a cylindrical joint, as depicted in Figure 20. In the initial configuration the dimensions indicated in Figure 20 are: H = 2.4606m, d ₁ = 16178m and d ₂ = 0.5523m.

All structural components are flexible and have square cross section with side a = 10.0cm. Twelve finite elements with cubic approximation were used to model the arch, which has an elastic modulus of 𝔼 = 200.0GPa. The crank is modelled by five finite elements and has elastic modulus ten times greater than the arch. The shear modulus for all components is adopted as half of the elastic modulus.

Figure 20:
Structure initial configuration.

The evolution of the reaction moment required to move the fixed end of the crank is presented in Figure 21a along with some deformed configurations of the system. The evolution of the normal force at the crank is presented in Figure 21b. The arch mid-span vertical position is depicted in Figure 21c. One can observe from those curves the occurrence of instabilities by limit points, indicated by the change of sign for the crank normal force and bending moment. At these positions the arch assumes an indifferent equilibrium configuration.

Figure 21:
Equilibrium path: a) reactive bending moment and structure configurations, b) crank mid-point normal force, c) arch mid-span vertical position.

Furthermore, from the discontinuity of the curves achieved at a rotation of ψ = 1.656rad it is clear the existence of the snap-back phenomenon. However, as it is known, it is not possible to describe the equilibrium path at snap-back situations using the Newton-Raphson method. The snap-back is present because the driven unstable shape is not a controlled one.

In order to describe this portion of the curve it would be necessary the use of an arc-length type method for the solution of the system that is beyond the objectives of this study.

Each step required at most 5 iterations to converge as shown in Figure 22. Only in the step that the snap-back occurred several iterations were needed. However, it is expected due to the large change of configuration of the arch at this point of the analysis.

Figure 22:
Evolution of the number of iterations for solution.

Therefore, from results, the potentialities of the proposed formulation and its consistency on evaluating structural behaviour are evident.

8 CONCLUSIONS

The prismatic and cylindrical sliding connections were successfully developed and implemented on a total lagrangian position based finite element method. Finite elements are curved of any order and connections slides according to space dimensional coordinates, not non-dimensional as presented in other works. The formulation is applied to solve four different examples demonstrating its accuracy, stability and possibilities of application. Large displacements and rotations are perfectly modelled by both the proposed finite element and connections when static applications are carried out. None problems are identified when solving the proposed applications revealing the robustness of the technique. Future developments include the extension of the proposed formulation to dynamic and 3D applications. In order to provide an insight for 3D extension of the proposed sliding formulation, the reader is referred to the work of Coda and Paccola (2010) that describes the frame unconstrained vector positional formulation. Thus, for a 'sliding spherical' joint a totally similar Lagrangian multiplier procedure as the described here for 2D cylindrical joint can be applied. For 3D cylindrical or prismatic joints, the Lagrangian multiplier technique should be associated with a penalty restraint among selected unconstrained vectors.

Acknowledgements

Authors would like to thanks FAPESP (São Paulo Research Foundation) and CAPES (Coordination for the Improvement of Higher Education Personnel) for the financial support of this research.

References

• Armero, F. (2006). Energy-dissipative momentum-conserving time-stepping algorithms for finite strain multiplicative plasticity. Computer Methods in Applied Mechanics and Engineering, 195(37-40), 4862-4889. http://doi.org/10.1016/j.cma.2005.09.028
  » http://doi.org/10.1016/j.cma.2005.09.028
• Bauchau, O. A. (1998). Computational Schemes for Flexible, Nonlinear Multi-Body Systems. Multibody System Dynamics, 2(2), 169-225.
• Bauchau, O. A. (2000). On the modeling of prismatic joints in flexible multi-body systems. Computer Methods in Applied Mechanics and Engineering, 181(1-3), 87-105. http://doi.org/10.1016/s0045-7825(99)00065-1
  » http://doi.org/10.1016/s0045-7825(99)00065-1
• Bauchau, O. A., & Bottasso, C. L. (2001). Contact Conditions for Cylindrical, Prismatic, and Screw Joints in Flexible Multibody Systems. Multibody System Dynamics, 5(3), 251-278.
• Bischoff, M., & Ramm, E. (2000). On the physical significance of higher order kinematic and static variables in a three-dimensional shell formulation. International Journal of Solids and Structures, 37(46), 6933-6960.
• Bonet, J., Wood, R. D., Mahaney, J., & Heywood, P. (2000). Finite element analysis of air supported membrane structures. Computer Methods in Applied Mechanics and Engineering, 190(5-7), 579-595. http://doi.org/10.1016/s0045-7825(99)00428-4
  » http://doi.org/10.1016/s0045-7825(99)00428-4
• Cardona, A., Geradin, M., & Doan, D. B. (1991). Rigid and flexible joint modeling in multibody dynamics using finite element. Computer Methods in Applied Mechanics and Engineering, 89(1-3), 395-418. http://doi.org/10.1016/0045-7825(91)90050-g
  » http://doi.org/10.1016/0045-7825(91)90050-g
• Coda, H. B. (2003). An exact FEM geometric non-linear analysis of frames based on position description. In. IN: 17H INTERNATIONAL CONGRESS OF MECHANICAL ENGINEERING. 2003, São Paulo. São Paulo: ABCM.
• Coda, H. B., & Paccola, R. R. (2010). Improved finite element for 3D laminate frame analysis including warping for any cross-section. Applied Mathematical Modelling, 34(4), 1107-1137.
• Coda, H. B., & Paccola, R. R. (2014). A total-Lagrangian position-based FEM applied to physical and geometrical nonlinear dynamics of plane frames including semi-rigid connections and progressive collapse. Finite Elements in Analysis and Design, 91, 1-15.
• Coda, H. B., Sampaio, M. S. M., & Paccola, R. R. (2015). A FEM continuous transverse stress distribution for the analysis of geometrically nonlinear elastoplastic laminated plates and shells. Finite Elements in Analysis and Design, 101, 15-33. http://doi.org/10.1016/j.finel.2015.03.004
  » http://doi.org/10.1016/j.finel.2015.03.004
• Garcia-Vallejo, D., Escalona, J. L., Mayo, J., & Dominguez, J. (2003). Describing rigid-flexible multibody systems using absolute coordinates. Nonlinear Dynamics, 34(1-2), 75-94. http://doi.org/10.1023/B:NODY.0000014553.98731.8d
  » http://doi.org/10.1023/B:NODY.0000014553.98731.8d
• Géradin, M., & Cardona, A. (2001). Flexible multibody dynamics: a finite element approach. Chichester: John Wiley & Sons.
• Greco, M., & Coda, H. B. (2006). Positional FEM formulation for flexible multi-body dynamic analysis. Journal of Sound and Vibration, 290(3-5), 1141-1174.
• Hong, D. F., & Ren, G. X. (2011). A modeling of sliding joint on one-dimensional flexible medium. Multibody System Dynamics, 26(1), 91-106. http://doi.org/10.1007/s11044-010-9242-7
  » http://doi.org/10.1007/s11044-010-9242-7
• Ibrahimbegović, A., & Taylor, R. L. (2002). On the role of frame-invariance in structural mechanics models at finite rotations. Computer Methods in Applied Mechanics and Engineering, 191(45), 5159-5176. http://doi.org/10.1016/s0045-7825(02)00442-5
  » http://doi.org/10.1016/s0045-7825(02)00442-5
• Jelenic, G., & Crisfield, M. A. (1996). Non-linear 'master-slave' relationships for joints in 3-D beams with large rotations. Computer Methods in Applied Mechanics and Engineering, 135(3-4), 211-228. http://doi.org/10.1016/0045-7825(96)01017-1
  » http://doi.org/10.1016/0045-7825(96)01017-1
• Jelenic, G., & Crisfield, M. A. (2001). Dynamic analysis of 3D beams with joints in presence of large rotations. Computer Methods in Applied Mechanics and Engineering, 190(32-33), 4195-4230. http://doi.org/10.1016/s0045-7825(00)00344-3
  » http://doi.org/10.1016/s0045-7825(00)00344-3
• Lee, S. H., Park, T. W., Seo, J. H., Yoon, J. W., & Jun, K. J. (2008). The development of a sliding joint for very flexible multibody dynamics using absolute nodal coordinate formulation. Multibody System Dynamics, 20(3), 223-237.
• Leyendecker, S., Betsch, P., & Steinmann, P. (2006). Objective energy-momentum conserving integration for the constrained dynamics of geometrically exact beams. Computer Methods in Applied Mechanics and Engineering, 195(19-22), 2313-2333. http://doi.org/10.1016/j.cma.2005.05.002
  » http://doi.org/10.1016/j.cma.2005.05.002
• Moon, Y., Choi, H. J., Seo, J. B., & Choi, J. (2014). Design and Kinematic Analysis of a New End-Effector for a Robotic Needle Insertion-Type Intervention System. International Journal of Advanced Robotic Systems, 11. http://doi.org/10.5772/59350
  » http://doi.org/10.5772/59350
• Nocedal, J., & Wright, S. J. (1999). Numerical optimization. New York: Springer.
• Norton, R. L. (2011). Design of machinery: an introduction to the synthesis and analysis of mechanisms and machines (5a ed.). Boston: McGraw-Hill.
• Ogden, R. W. (1984). Non-linear elastic deformations. Chichester: Ellis Horwood.
• Pascon, J. P., & Coda, H. B. (2013). A shell finite element formulation to analyze highly deformable rubber-like materials. Latin American Journal of Solids and Structures, 10(6), 1177-1209.
• Reis, M. C. J., & Coda, H. B. (2014). Physical and geometrical non-linear analysis of plane frames considering elastoplastic semi-rigid connections by the positional FEM. Latin American Journal of Solids and Structures, 11(7), 1163-1189.
• Romero, I., & Armero, F. (2002). An objective finite element approximation of the kinematics of geometrically exact rods and its use in the formulation of an energy-momentum conserving scheme in dynamics. International Journal for Numerical Methods in Engineering, 54(12), 1683-1716. http://doi.org/10.1002/nme.486
  » http://doi.org/10.1002/nme.486
• Shigley, J., & Uicker, J. (1981). Theory of machines and mechanisms. Auckland: McGraw-Hill.
• Simo, J. C., & Vu-Quoc, L. (1986). On the dynamics of flexible beams under large overall motions - the plane case. I. Journal of Applied Mechanics-Transactions of the Asme, 53(4), 849-854.
• Sugiyama, H., Escalona, J. L., & Shabana, A. A. (2003). Formulation of three-dimensional joint constraints using the absolute nodal coordinates. Nonlinear Dynamics, 31(2), 167-195. http://doi.org/10.1023/a:1022082826627
  » http://doi.org/10.1023/a:1022082826627
• Yoo, W. S., Kim, K. N., Kim, H. W., & Sohn, J. H. (2007). Developments of multibody system dynamics: Computer simulations and experiments. Multibody System Dynamics, 18(1), 35-58.
• Zaccaria, D., Bigoni, D., Noselli, G., & Misseroni, D. (2011). Structures buckling under tensile dead load. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 467(2130), 1686-1700. http://doi.org/10.1098/rspa.2010.0505
  » http://doi.org/10.1098/rspa.2010.0505

Publication Dates

• Publication in this collection
  Nov 2016

History

• Received
  25 Sept 2015
• Reviewed
  20 Apr 2016
• Accepted
  30 May 2016