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

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-VenantKirchhoff constitutive model. Selected examples are used to validate the formulation and to show its possibilities of application.


Latin American Journal of Solids and
tures 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 corotational 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 cylindri-Latin American Journal of Solids and Structures 13 (2016) 2059-2087 cal 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.

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 nondimensional variables.

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: where, i is the direction (1 or 2), m indicates reference line and  is the index that represents nodes and shape functions Latin American Journal of Solids and Structures 13 (2016) 2059-2087 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: in which comma indicates derivative.
From equation (2) the normal vectors that generate cross sections at nodes are written as: where the index m has been omitted.To use angles as nodal parameters, one calculates, from figure 1: 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.: 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. One can deduce from figure 2 that the initial configuration mapping is written using a nondimensional space ( ) , x h by: Latin American Journal of Solids and Structures 13 (2016) 2059-2087 0 ( , ) ( ) ( , ) where the vector 0 ( , ) i g x h can be written as a function of 0 ( ) q x , h and the height 0 h of the finite element as: Substituting equations ( 1) and ( 7) into equation ( 6) results the initial configuration mapping of the frame element, as: 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).

Current Configuration
The current configuration is understood as a trial position of the solution process and can be seen in figure 3.
Latin American Journal of Solids and Structures 13 (2016) 2059-2087 Therefore, the new nodal coordinates m i Y  and angles q  can be used to derive the current mapping from the non-dimensional space ( ) , x h to the current configuration as: where 1 ( , ) ( , ) 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 q  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.

Complete Mapping and Green Strain
After defining the mappings from the non-dimensional space to the initial 0 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 fig- ures 2 and 3 revealing the mapping f  as: where ( ) are defined by equations ( 8) to (11).Solids and Structures 13 (2016) 2059-2087 To calculate strain and strain energy it is not necessary to explicitly know f  , but its gradient, called here A. This is done in a very simple way (Bonet et al., 2000;Coda, 2003) as: in which For a trial current nodal position, both A 0 and A 1 are numerical values at a generic point ( ) , .
x h In the numerical strategy, coordinates ( ) x h 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: The Green-Lagrange strain E is adopted to develop the formulation, as it is an objective measure (Ogden, 1984), and is given by: In which I is the identity tensor of second order and C is the right Cauchy-Green stretch.

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: where P is the total energy of the system, e U is the stored elastic energy and P is the potential of conservative external forces.Therefore, the variational principle is stated by: Latin American Journal of Solids and Structures 13 (2016) 2059-2087 in which the symbol d means variation.

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: 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 being n 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: moreover, the constitutive elastic tensor C can be achieved as: 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: in which 0 V is the initial volume of the frame.

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: Latin American Journal of Solids and Structures 13 (2016) 2059-2087 In this expression F  is the conservative applied forces (including moments), q  is the vector of external distributed forces written as a function of nodal values as position of the reference line (Equations ( 10) and ( 11)) and 0 dS is the initial infinitesimal reference line length of a curved frame element.
To achieve the equilibrium equation one applies the principle of stationary potential energy on expression (24) using nodal positions as parameters, 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: ), the adopted surface force approximation (q ) and equation ( 21), one writes equation ( 26) as: in which the symbol Ä represents the tensor product.
It is usual to write equation ( 27) as (Coda;Paccola, 2014): in which int F  represents nodal forces and L is the matrix that transforms distributed forces into nodal equivalent ones (Coda;Paccola, 2014).As the variation Y d  is arbitrary, equation ( 28) results in the geometrical non-linear equilibrium equations for discrete frame analysis, as: (29)

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: The first derivative of the Green strain regarding positions is achieved from equations ( 15), ( 16) and ( 17) as: Latin American Journal of Solids and Structures 13 (2016) 2059-2087 As  one can simplify equation (31).From equation ( 15) it is possible to derive 1 A regarding nodal positions Y b a , for which a is the direction (1, 2 or 3, for the angle) and b is the element node, as follows: with summation over nodes  and k .
From expressions (30) through (33) the internal force is calculated using Gaussian quadrature over the initial volume, as: in which ig w and jg w are Gauss weights and ig x and jg h are integration points.

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 F  , as: In this expression g  is called the unbalanced mechanical force vector, as it is null only when the position Y  is the solution of the problem.
Equation ( 35) can be understood simply by finding Y  such that ( ) 0 . This statement is clearly non-linear regarding Y  .Therefore, to find a solution of the problem one expand the equality in a Taylor series truncated at the first order, Solving the linear system of equation ( 36) one finds a correction Y D  to the trial solution as: where the gradient of g  is the Hessian of the total potential energy.As the external forces are conservative, the Hessian matrix is given by: , a new trial position is achieved, as: This new trial solution is used again in equation ( 37) until:  is assumed to be the acceptable solution of (35).In equation (40) X  is the Euclidian norm of the initial position vector.

Hessian of the Proposed Frame Element
To complete the description, the frame element Hessian matrix H is presented as: In the last integral, the derivative of the Piola-Kirchhoff stress  is given by: 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.: with, Latin American Journal of Solids and Structures 13 (2016) 2059-2087 The non-null values of 2 1 / Y Y ¶ ¶ ¶ A   are restricted to the third direction (angles) and are given by: where b and z are element nodes and summation over  and k are present.
The integration over elements is given by Gaussian quadrature, as: for which ig w and jg w are weights and ig x and jg h are integration points.

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.
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).

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.q D , calculated at the initial configuration (figure 6), should be constant during the sliding process.These constraints are described by: and, 0 0 ˆ( ) ( ) for 3 or, in a general, form: where 3 i d is the Kronecker delta.
During sliding, the curvilinear position ( ) or, inversely, the non-dimensional coordinate ( ) will vary.
To impose the kinematical constraint of a prismatic joint via Lagrange multiplier one adds to the total potential energy the following potential: ( ) ( ) in which the term inside brackets is exactly the kinematical constraint given by equation ( 49), i l are the so called Lagrange multipliers, one for each direction i (summation implied), and L Y  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: in which Y  include all degrees of freedom, less the new ones ( ) . The principle of stationary potential energy is applied on expression (51) to find equilibrium equations, as: in which dP has already been detailed in previous sections and include variations of the strain en- ergy and applied external forces regarding current nodal positions Y  .
To complete the variation of ( , ) L Y L P   it is necessary to solve d .This is done regarding ( ) it is important to mention that the consulted works that deal with flexible connections use p x instead of p s as the main variable.
For prismatic connections, developing the derivatives, one writes in an open form: 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 L  assemblage.
In equation ( 54) the following property is used, ( ) ( ) ] The new equilibrium equation is now assembled respecting the degrees of freedom correspondence adding L  to int F  , i.e.: ( ) in which F  includes all external loads.
For cylindrical connections, the same procedure is used, remembering that the relative angle is free, resulting:

Non-Linear Solution
The non-linear set of equilibrium equations ( 57) is also solved by the Newton-Raphson procedure.
The unbalance mechanical vector L g  is written as: Latin American Journal of Solids and Structures 13 (2016) 2059-2087 For a trial solution ( ) equality (59) does not hold, therefore a Taylor expansion of first order is performed and the nullity is imposed: From this equation, the trial vector { } until a prescribed tolerance is respected, see section 4.This time, both the trial position and the matrix contain the usual contribution of connected and unconnected nodes and the new contribution of p s and i l .The achievement of p s is not sufficient to update L  and the Hessian matrix, as the function ( ) is not explicitly written.The solution of this stage is described in next item.
The new Hessian matrix , is achieved by the second derivative of ( , ) L Y L P   regarding positions and the new variables p s and i l .The part of the Hessian matrix that corresponds to the element connection degrees of freedom is achieved as: or, written in a way that the corresponding variables are identified: . The terms of 1 1 Latin American Journal of Solids and Structures 13 (2016) 2059-2087 in which index notation is applied.Index i corresponds to directions (1, 2 or 3, for prismatic con- nection), 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 con L H have already been calculated in matrix H and corresponds to the second derivative regarding nodal posi- tions of the strain energy related to the connected (path and sliding) elements.

Curvilinear and Non-Dimensional Variables
The definition of ( ) 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 L  , equation ( 58), and the Hessian matrix, equation ( 63), associated to the Lagrange multiplier technique are dependent of ( ) that is not explicitly defined.Therefore, it is necessary to calculate, for a given trial equilibrium position, the value of p x .This is done solving the following non-linear system of equations: that is the mechanical constraint present in equation ( 50), obviously with known values of ˆP i Y and i Y  , and represents, for both directions, the residue ( ) i P r x 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: This is expanded in a Taylor series of first order as: ( ) ( ) ( ) being 0

P
x a previously known trial value.When trying to solve equation ( 68) one realizes that, as e, the objectivity of equation ( 68) is lost at the solution.Due to this, the least square method uses the assumption ( ) 0 The non-dimensional variable is updated by respects a given tolerance.The terms in equation ( 70) are calculated by: in which the indexes  and m imply summation over the nodes from the active path element on both directions i .With P x known for a trial p s , implicit in the values of ˆP i Y and i Y  , the global solution process described by 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.

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): 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 0 A does not change its geometry, the internal efforts results: in which σ  is the rotated Cauchy stress for axis 1 y  and 2 y  along the cross section defined by the angle q .

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 8 1.0 10 - ⋅ .

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 revo- lute joint.The rotation of the crank imposes a motion on the mechanism as depicted in figure 8.  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.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.

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 0.25 m L = 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.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 0 1.0 mm h = and 0 25.0 mm b = .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.0 GPa =  .
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 0 f .Increasingly smaller values for the initial angle were adopted to show convergence to the analytical solution.The results for the rigid bar rotation f , , where I is the second area inertia moment, shows the efficiency of the formula- tion when handling geometrical nonlinear problems.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.

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 26.0 m u = (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.
A rectangular cross section with width 0 1.0 m b = and height 0 2.0 m h = is adopted for the bridge, and a square cross section with side 0.1m a = 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 .
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.All structural components are flexible and have square cross section with side 10.0 cm a = .
Twelve finite elements with cubic approximation were used to model the arch, which has an elastic modulus of 200.0 GPa =  . 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.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.Furthermore, from the discontinuity of the curves achieved at a rotation of  1.656 rad = 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 snapback 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.Therefore, from results, the potentialities of the proposed formulation and its consistency on evaluating structural behaviour are evident.

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.

(
Lagrange polynomials of any order).m i X  represents the nodal coordinates at the reference line and repeated indices indicate summation.

Figure 4 :
Figure 4: Change of configuration function or deformation function. 0

Figure 6 :
Figure 6: Difference of angles for a single node at different elements.
) being a the nodes of the active path element.One should note that variable L Y  belongs to Y  and will result in terms to be added in internal force and Hessian matrix at the solution process.At this, from the last expression the Newton method applied to an minimization problem when determining

Figure 8 :
Figure 8: Selected configurations and arm free end trajectory.

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

Figure 12 :
Figure 12: Initial and deformed configurations of the structure.
Figure 13, and right support horizontal displacement, Figure 14, against a non-dimensional traction force 2 2 4 / FL I p 

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

Figure 16 :
Figure 16: Vehicle initial configuration and geometry.

.
11.0 13.0 15.0 17.0 19.0 21.0 23.0 25.0 27.0 29.0 Shear force (kN) Load point horizontal position (m) Latin American Journal of Solids and Structures 13 (2016) 2059-2087 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 10The 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:

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

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