Improved friction model applied to plane sliding connections by a large deformation FEM formulation

Abstract Friction is an important source of dissipation in dynamical systems. Properly considering it in the numerical model is fundamental to obtain stable and representative responses in structures and mechanisms. This is especially significant for the well-known Coulomb model due to discontinuity in force when stick-slip transition occurs. In this work an improved friction force model is proposed to smooth the force transition at null velocity, with an additional parameter obtained from the own system state. The improved model is employed in sliding connections of plane frames finite elements. A total Lagrangian Finite Element Method (FEM) formulation based on a positional description of the motion is employed. Using a variational principle, frictional dissipation is added to the total mechanical energy to develop the equations of motion. The resulting nonlinear equations are solved by the Newton-Raphson method accounting for the friction force update in the iterative process. Examples are presented to show the formulation effectiveness and possibilities in simulating dynamical systems that present the stick-slip effect. Graphical Abstract


INTRODUCTION
In the analysis of structures and mechanisms their dynamical systems are frequently assumed conservative, while, in reality, systems present dissipation due to several sources.The frictional dissipative effect, in particular, is important to be considered when relative motion between parts of the body exists.This is the case of sliding connections, such as prismatic and cylindrical joints, that, by introducing translational movement among body members, allow friction forces to develop at their contact regions.
Mathematical models try to describe some of the many characteristics of friction complex nature.For instance, evolving from the classical one-parameter Coulomb friction model, values of the friction force at rest (static friction) and in motion (kinetic friction) can be distinguished.The continuous decrease of friction values for increasing, but low, velocities, referred as the Stribeck effect (Rabinowicz, 1951), or the adhesion between contacting surfaces before sliding, resulting in local pre-sliding deformation (Hsieh and Pan, 2000), may be important and included in some models.The effect of lubricant layers between the surfaces (Armstrong-Hélouvry et al., 1994) and the 'frictional lag', a delay in the change of the force with change in velocity (Hess and Soom, 1990), among other aspects, may be also encountered in current friction models.Comprehensive surveys (Olsson et al., 1998;Andersson et al., 2007;Marques et al., 2016) present in-depth discussions on the characteristics of friction dissipation and several available models.
Friction models can be broadly classified in 'dynamic' or 'static' (Marques et al., 2016).Dynamical models, also referred as state-variable models, aim to capture more characteristics of the physically observed friction response by including additional degrees of freedom in the system, and usually demanding more parameters.Common state-variable friction models are referred as Dahl (Dahl, 1976), LuGre (Canudas de Wit et al., 1995), Elasto-Plastic (Dupont et al., 2002), Gonthier (Gonthier et al., 2004), only to list a few.
Models referred as static, as opposed to dynamical models, dismiss the introduction of state variables to the problem, rendering a direct description of the force expression.In those models, however, the discontinuity of the friction force at null speed is a major drawback in the numerical solution.Some models (Kikuuwe et al., 2005;Threlfall, 1978;Armstrong-Hélouvry et al., 1994;Ambrósio, 2003) try to circumvent this problem suggesting null friction at a defined null speed interval.This assumption is not a good approximation when relative motion is intermittent, such as the sticking and slipping of sliding surfaces, known as the stick-slip effect (Rabinowicz, 1958), as it disregards the influence of the resultant force acting on the bodies in contact.Other models (Karnopp, 1985;Leine et al., 1998;Cha et al., 2011;Awrejcewicz et al., 2008) consider the resultant force requiring additional parameters, related to a null speed interval, for the transition between motion and rest states, which may be problem dependent and difficult to establish.
It is important to stress that, although several studies developed formulations to capture some of friction characteristics, no model is general (Pennestrì et al., 2016) and simpler force models are largely employed for design and are present in commercial software.Still, stability of the numerical solution for the formulation comprising the friction model is always desired.
In this study, we propose a modification on a static friction model and employ it in sliding connections of plane frames in a large deformation finite element context.The improved model is based on Coulomb friction considering the Stribeck curve for the static and kinetic forces transition and viscous effect.Its numerical solution is improved by reducing the abrupt transition between rest and motion states of the joints by an interpolation between the static friction value and the resultant force in a quasi-null relative speed interval.Differently from other friction models, that define a constant value for a null speed interval, here we propose an expression to automatically calculate the referred quasi-null interval from the system state.Thus, no additional non-physical parameter needs to be informed.The resultant force expression at the sliding connection is achieved for any generic finite element system with multiple degrees of freedom.With the proposed improved friction model, stable responses for the stick-slip effect are captured, which is important for the correct motion reproduction in structures and mechanisms.
The framework employed to model the dynamical system (Coda and Paccola, 2014;Siqueira and Coda, 2016;2017;2019) is a fully nonlinear finite element approach for large deformations based on a total Lagrangian description of solids.As the formulation uses positions as the main degrees of freedom, instead of displacements, the adopted approach is referred as positional.Sliding connections, as prismatic and cylindrical joints, are introduced in the system by means of Lagrange multipliers.In the sliding connection formulation, the friction force is associated with a curvilinear displacement variable introduced in the system by the joints.This formulation is also capable of representing roughness on the sliding surface, making it possible to be combined with the introduced friction dissipation to simulate, for instance, vehiclebridge interaction as show in the examples.The Saint-Venant-Kirchhoff constitutive model is adopted to define the plane frame elastic strain energy using the Green-Lagrange strain and the second Piola-Kirchhoff stress tensor.The Principle of Stationary Total Energy is used to find the equations of motion, comprising the frictional effect, which is discretized along Latin American Journal of Solids and Structures, 2023, 20(1), e472 3/23 time using the generalized-α method.The resulting nonlinear system is solved by the Newton-Raphson method accounting for the friction force update in the iterative process.
This study is organized as follows.Section 2 briefly presents necessary aspects of the employed nonlinear plane frame element followed by the sliding connections kinematical constraints in Section 3. The dynamical equilibrium is obtained in Sections, 4 and 5. Known the system parameters, the friction force is introduced in its variational form in Section 6 and the improved model is presented in detail by Section 7. Time integration (Section 8) and system solution (Section 9) are presented next.In Section 10, examples verify the proposed improved friction model and show the possibilities of the developed formulation in the analysis of dynamical systems.
Dyadic notation is preferred throughout this text due its brevity; however, index notation is also used to clarify particular aspects when necessary.

NONLINEAR PLANE FRAME KINEMATICS
The employed plane frame finite element is presented thoroughly in Coda and Paccola (2014) and Siqueira and Coda (2017), however, to develop the present formulation some aspects need to be briefly stated.As the finite element behavior is represented by a total Lagrangian description, its strain field needs to be written as a function of the initial and current configurations of the solid, restricted to a finite number of degrees of freedom.
In the positional approach of the FEM, instead of nodal displacements, the parameters of the discretized plane frame are its positions (coordinates) and the cross section angle (Fig. 1).The deformation function, f  , depicted in Fig. 2, can be written indirectly as function of the non-dimensional space and nodal parameters by mappings from the nondimensional space to the initial configuration, 0 f  , as (1) and to the current configuration, 1 f  , as where x  and y  represents any point on the domain of a finite element in the initial and current configuration, respectively.The coordinates for both directions 1, 2 i = of each node  along the reference line in the initial and current configurations are i X  and i Y  , respectively.The initial nodal value of the cross-section angle is 0 θ  and after deformation is denoted as θ  .The cross section initial height is 0 h , ξ is the non-dimensional space variable in the direction of the reference line and η follows the height direction.The shape functions ( ) φ ξ  are obtained by Lagrange polynomials of any order.
The deformation function can be written as a composition of the previous mappings, Eq. ( 1) and (2), as Latin American Journal of Solids and Structures, 2023, 20(1), e472 4/23 Since only the gradient A of the deformation function, but not the function itself, is necessary to obtain the strain field, as presented by (Bonet et al., 2000;Coda, 2003) where During the iterative solution strategy both 0 A and 1 A are numerical values calculated at the integration points resulting in a purely numerical procedure.
As the Saint-Venant-Kirchhoff constitutive law is employed, the Green-Lagrange strain tensor E have to be calculated.This objective strain is given, for instance, by (Ogden, 1984) where I is the second order identity tensor and  As there is no relation between the cross-section angle and the slope of the reference line, the frame kinematic can be regarded as Reissner-type.It should be mentioned that the cross-section dimensions are maintained the same during motion, thus, to avoid volumetric locking, the constitutive equation is relaxed in order to exclude transverse expansions.The specific strain energy for the plane frame is presented in Section 4.

KINEMATICAL CONSTRAINTS DUE TO SLIDING CONNECTIONS
To introduce frictional effects in the equations of motion, it is required first to write the constraint equations imposed by the sliding connections.The curvilinear position -a new variable introduced in the system -is of particular importance as the friction force will act directly on it.In this section we summarize the description of the connections as prismatic and cylindrical joints.
Sliding connections are the ones that constrain relative translations between parts of the body.Fig. 3 illustrates both joints and their plane representation.In either case, a sliding node, at which the joint exists, is constrained to move over a trajectory comprised of path elements.The distinction between the prismatic and the cylindrical joint (understood as a sliding-pin joint for the plane case) is the relative rotation, which is allowed only by the last one.The constraint equations, c  , can be written for both types of joints as a single expression where the upper bar ( ) • is used to identify variables related to path elements and the upper hat ( ) • is used to define the sliding node.In the previous expression i is the direction ( 1, 2, 3 i = for prismatic joints and 1, 2 i = for cylindrical joints) and ij δ is the Kronecker delta.The difference of cross sections orientations at the initial configuration is 0 0 0 θ θ θ ∆ = − and must be constant during the sliding process in the prismatic joint case to maintain a fixed relative angle.In Eq. ( 7), the components of the roughness profile, obtained by its 'height function' ( ) h r s , are given by 1 2 ( ) ( ) cos ( ) ( ) ( ) sin ( ) It is noteworthy that the curvilinear variable ( ) s ξ represents an arch-length function defined by the non-dimensional coordinate ξ and the path element coordinates.Also, as ( ) h r s is defined directly by the curvilinear position, the roughness profile is mesh-independent.

UNCONSTRAINED EQUATIONS OF MOTION
When dissipation occurs, the mechanical energy of the analyzed structure 0 Π is given by where  represents the dissipation of an encompassing system of total energy Π .Studying this larger system, Eq. ( 9) can be rewritten as or, making explicit the energy parcels of the new larger conservative system where  is the stored elastic strain energy,  is the potential of conservative external forces and  is the kinetic energy of the body.Following Lanczos (1970) and Gurtin et al. (Gurtin et al., 2010), it is not always possible to write down closed expressions for dissipative parcels but only its infinitesimal change.Thus, the equations of motion are stated from the variation of the energies present in Eq. ( 11), which is understood as the Principle of Stationary Total Energy in which the symbol δ means variation.The total energy can be stated by writing the known expressions of the energies in Eq. ( 11) as function of the current configuration nodal parameters of the discretized body, grouped in the vector ϒ  , as where the specific strain energy u depends on the strain state E of the body, Eq. ( 6), which is function of the nodal parameters ϒ  , as defined by the gradient of the deformation function in Eq. ( 4).Still in Eq. ( 13), F  and q  are the concentrated and distributed conservative external loads, respectively.The initial length of the frame reference line is 0 s .The material mass density in the initial configuration, of volume 0 V , is 0 ρ .The material point velocity is denoted using the over-dot as y   .As mentioned before, the Saint-Venant-Kirchhoff constitutive relation is employed due to its simplicity and good representation of large displacements on solids that remain in the small to moderate strain regimen, which comprehends the majority of the usual applications in engineering.For the plane frame utilized, following Siqueira and Coda (2017), its specific energy is given as where E is the longitudinal elastic parameter that approaches the Young modulus for small strains.The shear elastic modulus is , being ν a constant that reproduces the Poisson ratio for small strains.Assuming, in equation ( 14) 0 ν = , except for the calculation of the shear elastic modulus, avoids locking problems.The second Piola-Kirchhoff stress tensor is easily obtained by the energy conjugacy property as Latin American Journal of Solids and Structures, 2023, 20(1), e472  7/23 The equations of motion (geometric nonlinear dynamical equilibrium) are obtained by the development of the variations in eq. ( 13).Details regarding this development can be obtained for the employed plane frame element particular kinematics in Siqueira and Coda (2017) and for a general framework of positional finite elements in Siqueira and Coda (2019), for instance.Following these studies, the equilibrium can be written in a compact form as where: int ( )  is the internal force vector; F  collects all the external loads; M results in a constant mass matrix for the present total Lagrangian formulation; D is the external damping matrix arriving from δ which is assumed as Rayleigh type; and ϒ   and ϒ   are the velocity and acceleration vectors of the nodal parameters.

CONSTRAINED EQUATIONS OF MOTION
The dynamic equilibrium stated by Eq. ( 16) is called unconstrained as no restraints, such as the ones from the sliding connections, are present.Several consolidated methodologies to impose constraints on mechanical and structural applications can be found in the literature (Géradin and Cardona, 2001;Nocedal and Wright, 1999;Rao, 2009).Here, we employ the well-known Lagrange multiplier method along with the Principle of Stationary Total Energy to impose the sliding restrictions.Concerning the later introduction of friction dissipation, the multipliers are of great value since in Mechanics they might be understood as the contact forces between bodies, an essential information for the friction model.
The Principle of Stationary Total Energy is extended for the case of holonomic constraints by modifying the total energy through the introduction of a new potential  , referred as the constraint potential, as When using Lagrange multipliers, the expression of the new potential is simply given by where λ  represents the vector of multipliers, which are new variables of the system.Eq. ( 18) indicates the presence of a multiplier for each constraint equation in c  .It is worth mentioning that the constraint potential is null at the solution, therefore, the total energy is not altered.Knowing the expression of  , the first variation of the constrained energy, Eq. ( 17), is which, can be developed similarly to the unconstrained case leading to the constrained nonlinear equations of motion, expressed in a compact form as represents the restriction forces arriving from the constraint potential.As the multipliers are new variables, the variation of  is organized in the following force vector, which separates the parameters ϒ  (including P s ) and the multipliers where the tensor c ∇  represents the gradient of the constraint vector.The derivatives of the constraint equation for the sliding connections, Eq. ( 7), can be found in Siqueira and Coda (2017).

INTRODUCTION OF THE FRICTION FORCE IN THE SYSTEM
The friction force is included in the system directly in the Principle of Stationary Total Energy as part of the dissipative potential.As previously mentioned, dissipative potentials are introduced in their differential form since closed expressions might be unknown, as is the case for the dissipated friction energy f  .However, the variation of this potential can be written as the work done by the friction force vector f To develop Eq. ( 22), parameters that describes the force displacement must be chosen.For that, the coordinates of the sliding node and its path contact point could be chosen.However, since in the previous formulation the curvilinear position P s is already used as an intrinsic variable, the displacement along the trajectory is simply the scalar expression .
As the friction force acts tangentially to the trajectory, with its value given by the scalar f F , the dissipative parcel is introduced directly in the curvilinear position as To organize the equations of motion system, we make , the previous equation is rewritten as Considering the correspondence of the friction force vector f F  to the system degrees of freedom, the equations of motion are restated to include frictional dissipation as

IMPROVED FRICTION MODEL
Given the way the friction force is incorporated in the equations of motion, Eq. ( 25), any expression may be easily applied without modifying its development.We start from the Coulomb friction model considering the Stribeck curve for the static and kinetic forces transition and a viscous effect.Fig. 5 shows the overall behavior of the friction force with the relative velocity among sliding bodies, v .This model considers the Stribeck curve using the expression proposed by Bo and Pavelescu (1982) and a linear evolution for the viscous friction, which occurs if lubricant layers are present on the surfaces.The expression for the friction force is written as with the static S F and kinetic C F friction forces given by where, s µ and k µ are, respectively, the static and kinetic friction coefficients.N F is the absolute value of the contact force, orthogonal to the trajectory at the joint contact point, and R F is the resultant force on the sliding connection, tangential to the contact point.The viscous friction coefficient is η and the Stribeck parameters are its decay velocity v σ and power σ δ .The sign function is represented by sgn( )  .
Latin American Journal of Solids and Structures, 2023, 20(1), e472 9/23 We highlight at this point that, for our application in sliding connections, the relative velocity v , tangential to the path, can be measured from the curvilinear position velocity, which will be important for the nonlinear solution procedure in Section 9.
For null relative velocity, that is, the rest state, second condition in Eq. ( 26), the tangential resultant force acting on the connection is required for comparison with the static friction value to evaluate the tendency of motion (in case R F is greater than S F , or not, otherwise).Properly representing this condition is the key to simulate the stick-slip effect.In numerical simulations, the transition from motion to rest states is challenging due to the absence of continuity at null speed, as expressed by Eq. ( 26), leading to instabilities in the response.For this reason, a linear interpolation between S F and R F is proposed here for the stabilization of the friction force at a range 0 0 [ , ] v v − of quasi-null velocities, as depicted in Fig. 6.The improved friction model is written as where 0 v is the quasi-null speed limit.One should note that C F and S F are always positive as they are obtained from the absolute value of the normal force, Eq. ( 27), thus, the sign of the friction force in Eq. ( 28) depends on the values and signs of the relative velocity and resultant force.In the proposed approach, R F is not a constant value, but depends on the system own force state at a given time instant, which can even be null, if applicable.Therefore, the system response can be stabilized by means of a smooth transition from the motion state to rest state and vice versa.
It should be noted that 0 v is a new parameter introduced in the model to distinguish the movement when in the quasi-null velocity interval.For a better convergence of the iterative solution method, the recommended value of the limit velocity should be close to the relative stop speed of the bodies, but not too small, to still allow the smooth transition among forces at rest.Instead of adopting a fixed unknown arbitrary value, an estimative for this parameter can be obtained by being m the mass of the sliding node.This estimative is updated along the solution process and presented good responses for general applications, as is illustrated in the examples section.
Known the coefficients of the model, which depends on the materials that make the sliding connection and its path, the forces required to calculate the friction force have to be related to the variables describing the joint.The normal force vector N F  is found from the component of the Lagrange multipliers vector due to the translational constraints, , at the normal direction of the path at the contact point, defined by the normal vector P N  , as and its absolute value, actually used in the calculation, is In the plane case, the components of the normal vector are obtained from the tangent vector of the path finite element at the contact point, , ( ) and 2 1 The resultant force, equal to the inertial force at the sliding node, is obtained directly from the equilibrium equation, Eq. ( 25), considering only the sliding node degrees of freedom (positions and curvilinear variable), as or, to identify the terms referred to the degrees of freedom of interest one writes  s refer to the sliding node position degrees of freedom and the curvilinear position, respectively.In the definition of Eq. ( 32) and ( 33), being a quasi-null velocity case, the velocity-proportional external damping force was neglected.The friction force is also not present since its value is already considered indirectly through the constraint force at the curvilinear position direction.
As the tangential value of the resultant force R F is required, the tangent vector P T  is used to decompose the Cartesian terms as Latin American Journal of Solids and Structures, 2023, 20(1), e472 11/23 As expected from the physical meaning of the multipliers as contact forces, we have = 2 λ , which can be obtained by developing the constraint force given in Eq. ( 21) for the constraint equation in ( 7).

TIME INTEGRATION
For the time discretization, as the description of the solid is made in a total Lagrangian reference, the inertial force in Eq. ( 25) is obtained using a constant mass matrix.Similarly, the assumed Rayleigh damping also renders a constant matrix.This allows the use of classical time approximations, originally developed for linear dynamics, for the material velocity and acceleration vectors.Paultre (2011) discuss the application of direct time integration on nonlinear systems.
As Lagrange multipliers are present in the system to introduce constraints, a time marching scheme capable of controlling deleterious high frequencies introduced by their presence is necessary (Rong et al., 2019;Géradin and Cardona, 2001).Here we employ the generalized-α method (Chung and Hulbert, 1993) design to numerically control the system high frequency content while retaining the important information of low frequency range.
To discretize the solution in time increments t ∆ , the forces in Eq. ( 25) are rewritten for specific auxiliary time instants A linear combination of the type is used to evaluate any vector in ( 35) at an auxiliary instant using its current value and its past instant value t d  . The parameter α may represent f α or m α .Nodal parameters velocity and acceleration vectors are written using the Newmark approximations with parameters β and γ as and ( ) The time-discrete equations of motion results, using Eq. ( 36) and (37) in Eq. ( 35), ( )( ) ( ) ( ) with the vector that groups the terms of the previous time instant given by ( ) and Latin American Journal of Solids and Structures, 2023, 20(1), e472 12/23 1 1 2 At the end of the step, velocity and acceleration vectors must be calculated from the approximations ( 36) and (37).To start the time march, the initial acceleration is evaluated from Eq. ( 25) using the initial configuration and velocity of the system.Chung and Hulbert (1993) recommend obtaining the generalized-α method parameters from the high frequency region spectral radius . In order to retain second order accuracy, to minimize numerical dissipation of the fundamental frequencies and maximize the numerical dissipation on the high frequencies range the α-parameters are calculated from and the Newmark parameters by To achieve unconditional stability in the presence of constraints, the relation must also be respectedproven for the linear case in Géradin and Cardona (2001).As a consequence of m f α α < , the situation without numerical dissipation, , is not achieved.From the authors experience in most situations ρ ∞ in the range 0.7 -0.9 presents good results with small or negligible numerical dissipation even in nonlinear cases.

NONLINEAR SYSTEM SOLUTION PROCEDURE
To solve the resulting time-discrete nonlinear equations of motion, Eq. ( 38), for the variables { } where g  is the residual of the Newton method (or mechanical unbalanced vector), null when  ( ) , where X  is the initial nodal position vector of the body -is smaller than a given tolerance.It should be mentioned that the constraint equations, calculated from Eq. ( 7), can also be employed as a secondary convergence criterion.In the authors' experience, for the examples tested, employing the adopted relative position criterion with a tolerance of 1.10 -8 results in small force errors along the iterative process, even for the imposed constraints.
The Hessian matrix

H
(tangent operator), in Eq. ( 45), is given by ( ) ( ) ( ) ( ) Latin American Journal of Solids and Structures, 2023, 20(1), e472 13/23 The elastic part of the Hessian matrix, associated with the strain energy potential of the finite elements, is in which ⊗ represents tensor product.From the specific strain energy expression, Eq. ( 14), and developing the Green strain derivatives, using Eq. ( 6) and ( 2), one develops this matrix particular expression for the plane frame finite element as seen in Siqueira and Coda (2017) and Coda and Paccola (2014).
The Hessian matrix parcel due the constraint potential of the sliding connections is written as where, ( )  is a third order tensor that can be understood as the set of Hessian matrices due to each constraint equation i c , Eq. ( 7), and 0 is the null matrix.It must be stressed that, the achieved value of P s in the solution process is not enough to update con F  and is not explicitly written.The solution of this stage is done by adopting a least square method to find the nondimensional coordinate from the converged values of the path element and the sliding node as described in detail by Siqueira and Coda (2017).Given the numerical value of the non-dimensional variable in the dimensionless space, the solution process continues and also the transition between path elements can be performed -which is straightforward when P ξ exceeds the space domain As the friction force calculation is dependent on the tangential velocity P v s =  , Eq. ( 28), the solution procedure is improved by updating its value at each iteration by the following finite difference approximation ( ) ( ) where 1 k + , is explicitly written to identify the value as belonging to the current iteration of a time step 1 t + .Consequently, the friction force is entirely defined by the main solution variables, of particular interest, the curvilinear position. As , i.e., the friction force is a scalar properly placed at the system's degree of freedom, the friction part of the Hessian matrix can be developed as which is a scalar value to be added in the curvilinear position degree of freedom of the system's matrix.Further developing Eq. ( 50), using Eq. ( 49) and ( 28), results

(
) ( ) + updated in each iteration by Eq. ( 49).Regarding the computational implementation, to reduce the time spent in updating the residual vector and Hessian matrix of the Newton method for each iteration, it should be noted that, due to the total Lagrangian description of the Latin American Journal of Solids and Structures, 2023, 20(1), e472 14/23 positional formulation, the mass and damping matrices, as well as the part 0 A of the deformation gradient, Eq. ( 4), are constants and can be stored before starting the time marching process.

EXAMPLES
Examples are presented to verify the correct response of the proposed improved friction model and its introduction on the finite element system and also to show the capabilities of the proposed formulation in describing engineering applications of interest including roughness.

Rabinowicz test
The Rabinowicz test is a one degree of freedom experiment largely used as a benchmark, Pennestrì et al. (2016), for testing friction models, particularly useful in verifying the model capability of representing the stick-slip behavior.Fig. 7a) depicts the test system with mass m , held by a spring with stiffness k , in contact with a belt moving horizontally with constant velocity v .Due to friction between the oscillator and belt, the mass initially sticks to the belt until the spring force exceeds the static friction value.At this point sliding takes place under the kinetic value of friction.After a while, the mass attaches again to the moving surface and the process is repeated.Parameters for the mass-spring system and friction were adopted the same as Pennestrì et al. (2016).To simulate the oscillator, one two-node (linear) frame element, Fig. 7b), is employed with a cylindrical joint at the end of the bar.This connection is free to slide over another finite element moving in the horizontal direction with velocity v = 0.5 m/s and with all other degrees of freedom locked.A vertical force P = 196.2N is applied to manifest frictional effects at the contact point.In order to maintain small strains, to allow comparison with the reference, a L = 1000 m bar is adopted with squared cross-section, 0 b = 0 h = 1.0 m, and Young modulus E = 10 kPa.This results in an equivalent axial stiffness k = 0 0 / b h L E = 10 N/m.Null Poisson ratio is used and the equivalent system mass m = 20 kg is lumped at the joint node.
The sliding node displacement is shown in Fig. 8a) for the improved friction model under two different cases.The first considers the Stribeck effect, with parameters s µ = 0.6, k µ = 0.5, v σ = 0.05 m/s and σ δ = 1.0, and the second case is the classic Coulomb model, only capable of representing kinetic friction.No viscous effect is present, η = 0, in both cases.Results were obtained with a time increment t ∆ = 1.0 ms with the generalized-α method parameter ρ ∞ = 1.0.Fig. 8b) depicts Pennestrì et al. (2016) results for the mass displacement due to different friction models evaluated in their work.No direct comparison can be performed with the reference models as they are not directly equivalent with the proposed improved model, or among themselves, as they consider different friction characteristics requiring parameters not present in the current model.For those parameters and models description we refer to Pennestrì et al. (2016) study.
Compatible responses with Pennestrì et al. (2016) models are observed for both friction cases even for the relative velocity, Fig. 9, and friction force, Fig. 10.It is noticeable that, as the classic Coulomb model does not account for static friction, after the initial stiction phase, the system displays the harmonic oscillator behavior without maintaining null relative velocity and stiction again.For the Stribeck case a detail for the first relative motion of the system can be seen in Fig. 11 showing the variation between static and kinetic friction force captured by the improved model.Regarding processing time between models, for this example, the simulation for the Stribeck case is 15% longer than for the classic Coulomb model.The computation time was calculated by averaging the time spent in three separate simulations of each case (due to the influence of the operational system background tasks) on the same computer.
To further analyze the improved friction model results, a reference solution is depicted in Fig. 8 to Fig. 11 for each case.As analytical solutions are not always available, the reference solution was obtained by employing a time step 100 times smaller to diminish the influence of the quasi-null speed range calculated by Eq. ( 29).For the classic Coulomb case, the reference solution after stiction is also the response of a harmonic oscillator subjected to the Couloumb friction force.The simulated response coincides perfectly for this case (overlapping curves).A difference in the response is observed after the initial stiction phase in the Stribeck case.This difference is anticipated as the friction force value in this situation Latin American Journal of Solids and Structures, 2023, 20(1), e472 15/23 is velocity dependent, due to the exponential term of the Stribeck curve expression, Eq. ( 26), which is null for the Coulomb case.Thus, the numerical solution procedure has a higher influence.Nevertheless, very good agreement with Pennestrì et al. (2016) models are observed.

Quick-return mechanism with friction
A classic quick-return mechanism is analyzed considering frictional dissipation, Fig. 12.This mechanism was proposed by Bauchau (2000) and also investigated by Siqueira and Coda (2017)   All bars have 5.98 cm squared cross section and density 1724.82kg/m 3 .The arm is flexible with E = 47.04 GPa and discretized by 10 cubic finite elements.To simulate the rigid behavior of the remaining bars a 10 6 times greater Young modulus was adopted and 2 cubic frame elements were used for each bar.The shear elastic modulus is half the value of its corresponding Young modulus value for the whole mechanism.Lumped masses are placed at slider S (0.31 kg) and N (2.50 kg).Simulations were performed with a time increment t ∆ = 0.1 ms and the generalized-α method parameter ρ ∞ = 0.9.
Four cases of friction were evaluated.As we are interested in the flexible arm deflection and the output velocity of the mechanism, attention is restricted to slider N. Tests with friction at slider S presented minor effects in the overall system response, thus dissipation at this joint was not considered.Two dry friction cases, one with s µ = 0.The first four cycles of the mechanism arm tip deflection (point A), evaluated orthogonally to a line that follows the arm spin at point B, are shown in Fig. 13.Interesting results for the deflection of the arm were obtained for different friction cases.Compared to the frictionless case, the situation with smaller dry friction coefficients presented slightly higher arm deflection along the cycles, Fig. 13a), but similar output velocity at slider N, Fig. 14a), even though a horizontal friction force, Fig. 15a), acts at this point.The output velocity for a frictionless rigid body case of the mechanism (simulated by adopting the same Young modulus for the arm AB as for the rigid parts) is also depicted in Fig. 14.It is clear that the arm flexibility promotes an oscillatory response of the output velocity around the rigid body response.
For the higher dry friction parameters adopted, however, the mechanism presents a curious behavior reducing its arm tip vibration amplitude in alternating cycles, Fig. 13b).The output velocity presented the same pattern, Fig. 14b).Friction forces, Fig. 15b), showed smaller values at the low amplitude cycles, than higher ones.Tests performed for even higher values of static and kinetic friction coefficients ( s µ = 0.85 and k µ = 0.80) showed that the same pattern occurs, but with less intense reduction of the deflection among cycles.Consequently, we understand that exist a range of friction parameters, in combination with this particular mechanism mass and stiffness, resulting in a very particular cyclealternating nonlinear response.The physical interpretation of this response can be realized as slider N friction force Latin American Journal of Solids and Structures, 2023, 20(1), e472 17/23 affecting link NA force transmission to the flexible arm AB in a way that excites its transverse deflection at each two cycles and in the opposite direction of the crank (constant) excitation.This indicates the existence of a transverse motion along the arm (with a period twice as long as the duration of the mechanism cycle) decreasing the contact force in a cycle and increasing it on the next.Cases with only viscous friction presented the expected reduction of vibration of the arm deflection, Fig. 13c) and d), and slider velocity, Fig. 14c) and d).For the larger friction coefficient case, although vibrations were quickly damped, higher values of arm deflection were observed indicating the difficulty of the crank to move the arm for the intensity of the friction force developed at slider N, Fig. 15d).
Independently of the friction force evolution aspect, Fig. 15, smoother for the velocity proportional viscous cases or discontinuous in dry friction cases, the mechanism response was correctly depicted showing sense alternance for its return and that the proposed improved friction model was capable of producing limited solutions for a system with several degrees of freedom.From the vehicle midpoint horizontal position, Fig. 20, it can be observed that it stopped before the central support in all three cases where friction was considered.Thus, the bridge left mid-span deflection is not alleviated as for the frictionless situations, Fig. 21, when the vehicle passes on the right span.The vertical movement, manifested in the bridge even in the situation without roughness, is due to the horizontal friction force not aligned with the center of gravity of the vehicle, which generates a low frequency rotational vibration causing vertical movement in the region of the wheels, resulting in vertical forces transferred to the bridge.This aspect could only be captured due to the geometrically nonlinear description of the employed formulation.

CONCLUSIONS
Friction dissipation was successfully introduced in sliding connections commonly present in structures and mechanisms employing a total Lagrangian FEM formulation for plane frames based on a positional description of its kinematics.An improvement on a Coulomb friction model with Stribeck effect and viscous friction was proposed for a smoother description of the transition between motion and rest states of the joints.The improved model avoided instabilities for the friction force at null speed and is able to correctly represents frictional dissipation and capture the stick-slip effect since the resultant force could be calculated properly for the connections in the finite element system.The improved friction model is coupled with a mesh independent roughness strategy enabling general applications.Benchmarks and additional examples are used to verify the proposed model and to analyze engineering applications.Results were presented with successful responses for dynamical systems in which the consideration of friction is important.Future studies intent to extend this formulation to 3D applications.
is the right Cauchy-Green stretch tensor.

Fig. 4
Fig. 4 depicts a sliding connection, belonging to node P , and its path contact point P .The connection is free to move along the path ( s ξ ) defined by path finite elements, which may have an arbitrary roughness profile ( ) r s  .The new variable ( ) P P s s ξ =

Fig. 7 .
Fig. 7. Geometry of the test: a) mass-spring and b) finite element bar scheme.

Fig. 8 .
Fig. 8. Displacement for the Rabinowicz test: a) improved model results and b) numerical results from Pennestrì et al. (2016) for different friction models.

Fig. 9 .
Fig. 9. Relative velocity for the Rabinowicz test: a) improved model results and b) numerical results from Pennestrì et al. (2016) for different friction models.

Fig. 10 .Fig. 11 .
Fig. 10.Friction force for the Rabinowicz test: a) improved model results and b) numerical results from Pennestrì et al. (2016) for different friction models.
employing positional frame finite Latin American Journal of Solids and Structures, 2023, 20(1), e472 16/23elements, for the frictionless case.The system consists in a 1.0 m arm AB able to spin around a fixed support B and connected to a 0.25 m bar NA.The motion is imposed by a 0.20 m crank RS which turns around the fixed support R with a constant angular speed Ω = 5π rad/s and slides over the arm with a cylindrical joint S. The translational motion resulting from the mechanism input rotation is obtained by a cylindrical joint at point N, able to move horizontally over a fixed path.

Fig. 15 .
Fig. 15.Friction force at slider N for the quick-return mechanism.