Dynamical analysis of sliding connections with mesh independent roughness by a total Lagrangian FEM

Abstract

Sliding connections are present in several applications on the mechanics, civil and aerospace industries. A framework consisting on an accurate and stable formulation to describe the dynamics of flexible systems with sliding connections is developed. The total Lagrangian positional approach of the Finite Element Method is employed using 2D solid and frame elements to discretize bodies and connections. This allows a wide range of applications, particularly the local modelling of joints. The proposed formulation includes roughness along sliding paths independent from the finite element geometry discretization. Following variational principles, Lagrange multipliers are used to impose sliding constraints on the equations of motion. A direct time integration is performed by the generalized-α method and its stability in the present finite deformation context is evaluated. The resulting nonlinear equations are solved by the Newton-Raphson method. Examples are presented where the proposed framework is evaluated regarding its dynamical behavior and to solve practical scenarios for which sliding modelling is a necessity.

Keywords:
Flexible multibody system; Sliding rough connections; Constrained nonlinear dynamics; Finite deformation; Direct time integration

Graphical Abstract

1 INTRODUCTION

Sliding connections are part of several structures and mechanisms. They are employed to constrain contact between bodies by allowing relative sliding among their surfaces without detachment. Constraints imposed by sliding connections can be used in the local modelling of several types of joints - such as revolute, spherical, prismatic, cylindrical and plane - which are used to model a variety of applications on the mechanics, civil and aerospace industries. For instance, sliding connections are necessary in joints present in machinery and satellites, base isolation devices in buildings, vehicle-bridge interaction and cranes, only to list a few.

This study proposes a numerical framework that consists on an accurate and stable formulation and solution strategy to describe the dynamics of flexible systems with sliding connections by employing the total Lagrangian positional approach of the Finite Element Method (FEM). In this technique, positions, instead of displacements, are used as the main variables to discretize solids and structures. It has been successfully used in several works with good results for unconstrained dynamic analyses (Coda et al., 2020; Siqueira and Coda, 2019; Coda and Paccola, 2014; Coda and Paccola, 2011; Coda et al., 2013). As far as the authors knowledge goes, for constrained cases the positional FEM has been used only for bilateral sliding connections of plane frames using the Newmark-β time integrator (Siqueira and Coda, 2016; Siqueira and Coda, 2017). The main contributions of the work are: (i) a formulation capable of describing sliding surface contact between 2D solid finite elements where mesh independent roughness exists in the contact trajectory, this allows the local modelling of non-ideal joints in a wide range of applications; and (ii) the use of the generalized-α time integrator to control spurious frequencies that appear in systems with constraints imposed by Lagrange multipliers.

Since in the proposed formulation, roughness or path imperfections are not dependent on the finite element discretization, less refined meshes can be used for the analyzed bodies. As far as the authors are aware, roughness in sliding interfaces is mostly solved by equivalent forces and simplified mass/spring/damper systems (Oliva et al., 2013; Yin et al., 2010; Ding et al., 2009; Lyubicheva and Tsukanov, 2022) or by employing very refined finite element meshes and unilateral contact simulations (Bonari et al., 2022; Li et al., 2022), making the proposed approach an important alternative for this kind of simulation.

Differently from the proposed approach, the large majority of continuum mechanics studies interested in modelling sliding connections adopts the updated Lagrangian technique, mainly co-rotational formulations (Jelenic and Crisfield, 2001; Cardona et al., 1991; Verlinden et al., 2018; Gay Neto and Wriggers, 2020). This approach uses intermediate frames of reference to characterize large displacements and results in variable mass matrices for the finite elements, usually requiring special time integration schemes to solve the equations of motion (Espath et al., 2015; Hilber et al., 1977; Mamouri et al., 2014; Simo et al., 1992; Leyendecker et al., 2008; Pimenta et al., 2008). On the other hand, formulations based on the total Lagrangian description use only one frame of reference, resulting (for solid elements) in nonlinear equations of motion with constant mass matrix. Thus, the same time integration methods employed in geometrically linear dynamic analysis can be used. Paultre (2011), for instance, clarify the use of direct time integration schemes in nonlinear systems.

The good choice of the time marching method is an important aspect to be considered when solving dynamical systems with constraints, see Rong et al. (2019) for instance. The algebraic equations introduced by Lagrange multipliers constraints in the finite element system leads to numerical instabilities, even for linear dynamics (Géradin and Cardona, 2001). To overcome this problem, time integrators that control the total energy of the system and methods that introduce numerical dissipation are frequently used. The former usually requires changes in the finite element formulation (Simo and Wong, 1991; Simo et al., 1992; Kuhl and Ramm, 1999; Leyendecker et al., 2006; Leyendecker et al., 2008) and is well suited to co-rotational formulations, while the latter can be applied in total Lagrangian formulations without changes.

As the total Lagrangian FEM positional approach results in equations of motion with constant mass matrix, the generalized-α method (Chung and Hulbert, 1993) may be directly employed. This method (developed for the time marching of linear dynamics) is able to introduce numerical dissipation in a controlled manner. Here, its capability to obtain stable and representative responses is investigated for the nonlinear constrained (with roughness) dynamics case. The successful application of the generalized-α method in the present context, and the favorable results obtained, reflects the importance of the above-mentioned contributions of the present study.

Since the proposed applications include finite deformations but small strains, the Saint-Venant-Kirchhoff constitutive model is used, which linearly relates the Green-Lagrange strain with the second Piola-Kirchhoff stress tensor. Also, to describe the dynamic response in a more realistic manner, structural dissipative aspects are considered by employing a Rayleigh-type damping adapted to the present geometrical nonlinear formulation.

The work is organized as follows: Section 2 shows the nonlinear constrained equations of motion using the principle of stationary total energy, Section 3 presents the time discretization employing the generalized-α method. The resulting nonlinear system of equations is solved by the Newton-Raphson method in Section 4. The kinematic descriptions for the employed solid and frame 2D finite elements are described in Section 5. The subsequent section focusses on the development of the sliding connections between the plane elements and the algorithm used for its contact point identification. Section 7 detaches the coupling strategy between the different element types making possible contact path definition and frame / solid transition. In Section 8 three examples are employed to demonstrate the use of the proposed sliding connection framework. The modeling of engineering applications, where the occurrence of sliding is necessary to represent its dynamical behavior, is presented and results are evaluated. Conclusions are discussed in Section 9.

Dyadic notation is preferred throughout the text, but some expressions are made clearer by using index notation.

2 EQUATIONS OF MOTION

The general procedure for obtaining the equations of motion of the mechanical-structural system containing constraints is presented in this section. In subsequent sections, essential details of the employed finite elements and the sliding connections constraints are presented.

Following Siqueira and Coda (2017), the principle of stationary total energy, which is common in nonlinear analysis (Holzapfel, 2000), is used to represent the finite deformation behavior of solids in the employed total Lagrangian environment.

2.1 Dynamic nonlinear equilibrium for unconstrained systems

A mechanical system without constraints is initially considered. The strain energy and external load potentials, as well as the kinetic energy and the dissipative damping energy, are included in the total energy of the system, which is written as

$$\Pi ={\displaystyle {\int }_{V_0}^{}u(E)d{V}_0}-\overrightarrow{F}\cdot \overrightarrow{Y}-{\displaystyle {\int }_{s_0}^{}\overrightarrow{q}\cdot \overrightarrow{y}d{s}_0}+\frac{1}{2}{\displaystyle {\int }_{V_0}^{}{\rho }_0\dot{\overrightarrow{y}}\cdot \dot{\overrightarrow{y}}d{V}_0}+\mathcal{Q}(\dot{\overrightarrow{Y}})$$ (1)

in this expression the position vector $\overrightarrow{Y}$ collects all the nodal parameters used to describe the finite elements, $u$ is the specific strain energy written as function of the Green-Lagrange strain tensor $E$, $\overrightarrow{y}$ and $\dot{\overrightarrow{y}}$ are the position and velocity vectors of the material points and $\overrightarrow{F}$ and $\overrightarrow{q}$ are the concentrated and distributed external conservative loads, respectively. For the 2D solid element, $s_0$ represents its initial area whereas for the frame element $s_0$ is its reference line. Due to the total Lagrangian aspect of the formulation, $s_0$ is referred to the initial (reference) configuration of the body as well as the mass density ${\rho }_0$ and the volume $V_0$. The quantity $\mathcal{Q}$ represents the amount of energy dissipated from the mechanical system (damping for instance).

At equilibrium, the total energy expressed in Eq. (1) is stationary for any variation of the current position of the body; which is written as

$$\delta \Pi =0$$.(2)

As the degrees of freedom describing the mechanical system are united in one vector $\overrightarrow{Y}$, the variation of the total energy can be developed as

$$\begin{array}{l}\delta \Pi ={\displaystyle {\int }_{V_0}^{}\frac{\partial u}{\partial E}:\frac{\partial E}{\partial \overrightarrow{Y}}\cdot \delta \overrightarrow{Y}d{V}_0}-\overrightarrow{F}\cdot \delta \overrightarrow{Y}-{\displaystyle {\int }_{s_0}^{}\overrightarrow{q}\cdot \frac{\partial \overrightarrow{y}}{\partial \overrightarrow{Y}}\cdot \delta \overrightarrow{Y}d{s}_0}+\\ +{\displaystyle {\int }_{V_0}^{}{\rho }_0\frac{d\dot{\overrightarrow{y}}}{dt}\cdot \dot{\overrightarrow{y}}\cdot \delta td{V}_0}+\frac{\partial \mathcal{Q}}{\partial \overrightarrow{Y}}\cdot \delta \overrightarrow{Y}=0\end{array}$$.(3)

Due to the energy conjugacy property (Ogden, 1984), the derivative of the specific strain energy regarding the Green strain is the second Piola-Kirchhoff stress tensor $S$. Also, noticing that $\dot{\overrightarrow{y}}\cdot \delta t=\delta \overrightarrow{y}$ the kernel of the variation of kinetic energy in Eq. (3) becomes ${\rho }_0\ddot{\overrightarrow{y}}\cdot \delta \overrightarrow{y}$. The dissipative energy expression due to damping is not usually known (Lanczos, 1970; Gurtin et al., 2010), however its derivative $\partial \mathcal{Q}/\partial \overrightarrow{Y}$ represents a known damping force vector that acts on the system degrees of freedom and may be assumed to exist. Considering those remarks, from Eq. (3), one writes

$$\begin{array}{l}\delta \Pi ={\displaystyle {\int }_{V_0}^{}S:\frac{\partial E}{\partial \overrightarrow{Y}}\cdot \delta \overrightarrow{Y}d{V}_0}-\overrightarrow{F}\cdot \delta \overrightarrow{Y}-{\displaystyle {\int }_{s_0}^{}\overrightarrow{q}\cdot \frac{\partial \overrightarrow{y}}{\partial \overrightarrow{Y}}\cdot \delta \overrightarrow{Y}d{s}_0}+\\ +{\displaystyle {\int }_{V_0}^{}{\rho }_0\ddot{\overrightarrow{y}}\cdot \delta \overrightarrow{y}d{V}_0}+\frac{\partial \mathcal{Q}}{\partial \overrightarrow{Y}}\cdot \delta \overrightarrow{Y}=0\end{array}$$.(4)

In order to solve position variation in Eq. (4) the finite element kinematics must be known. The specificities of each finite element are discussed in the next sections, however in the total Lagrangian FEM positional approach (Coda et al., 2013; Coda and Paccola, 2007), the previous equation may be written as follows

$$\delta \Pi ={\overrightarrow{F}}^{\text{int}}\cdot \delta \overrightarrow{Y}-\overrightarrow{F}\cdot \delta \overrightarrow{Y}+M\cdot \ddot{\overrightarrow{Y}}\cdot \delta \overrightarrow{Y}+D\cdot \dot{\overrightarrow{Y}}\cdot \delta \overrightarrow{Y}=0$$ (5)

where the first integral in Eq. (4) represents the elastic internal forces grouped in the internal nodal force vector ${\overrightarrow{F}}^{\text{int}}$. The concentrated and distributed loads (as nodal equivalents) are gathered in the same vector $\overrightarrow{F}$. The mass matrix .$M$., as will be revealed by Eq. (53), is constant for the adopted elements. The damping force is assumed to be Rayleigh-type velocity proportional with damping matrix $D$. This matrix is also constant as will be shown by Eq. (54) when presenting the employed finite elements.

The equations of motion for the geometrical nonlinear mechanical system without constraints arise from the arbitrary character of variations $\delta \overrightarrow{Y}$ in Eq. (5), resulting

$${\overrightarrow{F}}^{\mathrm{int}}-\overrightarrow{F}+M\cdot \ddot{\overrightarrow{Y}}+D\cdot \dot{\overrightarrow{Y}}=\overrightarrow{0}$$.(6)

2.2 Constraint introduction

When kinematic constraints are added to the system of equation (6), the search for the solution can be interpreted as a constrained optimization problem.

In this work, constraint imposition is done by Lagrange multipliers by a set of holonomic constraints. The vector that groups the constraint equations $\overrightarrow{c}$ to be introduced to the system can be written as

$$\overrightarrow{c}\left(\overrightarrow{Y},t\right)=\overrightarrow{0}$$ (7)

in which $\overrightarrow{Y}$ groups the previous finite element nodal parameters and additional variables introduced by the constraint equation itself.

In order to introduce constrains into the principle of stationary energy, the following constraint potential is added to Eq. (1)

$$\mathcal{C}=\overrightarrow{\lambda }\cdot \overrightarrow{c}$$ (8)

where the vector $\overrightarrow{\lambda }$, which has the same dimension as $\overrightarrow{c}$, gathers the multipliers for each constraint equation. As the multipliers represent new degrees of freedom of the system, we write the variation of the constraint potential as a function of the nodal parameters and the multipliers as

$$\delta \mathcal{C}=\delta \overrightarrow{Y}\cdot \frac{\partial \overrightarrow{c}}{\partial \overrightarrow{Y}}\cdot \overrightarrow{\lambda }+\delta \overrightarrow{\lambda }\cdot \overrightarrow{c}=\left\{\begin{array}{cc}\delta \overrightarrow{Y}& \delta \overrightarrow{\lambda }\end{array}\right\}\cdot \left\{\begin{array}{c}\nabla \overrightarrow{c}\cdot \overrightarrow{\lambda }\\ \overrightarrow{c}\end{array}\right\}=\left\{\begin{array}{cc}\delta \overrightarrow{Y}& \delta \overrightarrow{\lambda }\end{array}\right\}\cdot {\overrightarrow{F}}^{\text{con}}$$ (9)

where $\nabla \overrightarrow{c}$ is the gradient of the constraint vector regarding nodal parameters, that is, its Jacobian matrix, and ${\overrightarrow{F}}^{\text{con}}$ is the vector representing the forces that impose constraints on the system (arising from the introduced potential). As $\delta \overrightarrow{\lambda }$ is also arbitrary, including Eq. (9) into Eq. (5) results the geometric nonlinear dynamic equilibrium with constraints, as

$$\left\{\begin{array}{c}{\overrightarrow{F}}^{\mathrm{int}}\\ \overrightarrow{0}\end{array}\right\}-\left\{\begin{array}{c}\overrightarrow{F}\\ \overrightarrow{0}\end{array}\right\}+\left\{\begin{array}{c}M\cdot \ddot{\overrightarrow{Y}}\\ \overrightarrow{0}\end{array}\right\}+\left\{\begin{array}{c}D\cdot \dot{\overrightarrow{Y}}\\ \overrightarrow{0}\end{array}\right\}+\left\{\begin{array}{c}\nabla \overrightarrow{c}\cdot \overrightarrow{\lambda }\\ \overrightarrow{c}\end{array}\right\}=\overrightarrow{0}$$ (10)

or, known the correspondence of nodal parameters and multipliers, the previous equation can be expressed in compact form, better suited for further developments, as

$${\overrightarrow{F}}^{\mathrm{int}}-\overrightarrow{F}+M\cdot \ddot{\overrightarrow{Y}}+D\cdot \dot{\overrightarrow{Y}}+{\overrightarrow{F}}^{\text{con}}=\overrightarrow{0}$$.(11)

Although the multiplier method introduces more unknowns in the system, constraints are exactly imposed without any additional parameter as, for example, penalty coefficients. In addition, multipliers have a direct mechanical interpretation as the internal force necessary to impose constraints. In the present case of sliding connections, they represent the contact force between surfaces. More details regarding the adopted constraint equation are given in Section 6. It is interesting to comment the physical interpretation of Eq. (7) enforced by using equations (8) and (10). In the proposed formulation, the constraint $\overrightarrow{c}$ is the relative position between the sliding point and the sliding surface that depends on the curvilinear coordinate, defined along the sliding path, and the roughness function, see Eq. (56). In this sense, when Eq. (10), or (11), is solved, the Lagrange multipliers result the necessary force to guarantee the imposed constraint.

3 TIME INTEGRATION

From the equations of motion, Eq. (10), one can note that the introduction of constraints by multipliers result in a system with the original set of time differential equations with additional algebraic equations. As discussed by Géradin and Cardona (2001), even for linear cases, the massless multipliers related with the algebraic equations constraints are associated with infinite vibration frequencies, leading to numerical instabilities on the time marching response. It has been proven (Cardona and Géradin, 1989) for linear systems that the Newmark trapezoidal rule is unconditionally unstable in the presence of constraints, rendering unstable responses particularly in long term simulations. One possible approach to deal with this problem is the introduction of controlled numerical dissipation in the numerical time integration process.

In the present study, we use the generalized-α method time integration (Chung and Hulbert, 1993) in order to introduce numerical dissipation in a controlled and optimal manner. In order to be complete, in the example section we show that the introduction of numerical dissipation using the Newmark trapezoidal rule is less able to maintain the response quality when compared to the generalized-α method.

The generalized-α method is employed to discretize the solution in time increments $\Delta t$ by rewriting the forces in Eq. (11) for specific auxiliary time instants $t+1-{\alpha }_f$ and $t+1-{\alpha }_m$ as

$${\overrightarrow{F}}_{t+1-{\alpha }_f}^{\mathrm{int}}-{\overrightarrow{F}}_{t+1-{\alpha }_f}+M\cdot {\ddot{\overrightarrow{Y}}}_{t+1-{\alpha }_m}+D\cdot {\dot{\overrightarrow{Y}}}_{t+1-{\alpha }_f}+{\overrightarrow{F}}_{t+1-{\alpha }_f}^{\text{con}}=\overrightarrow{0}$$.(12)

A linear combination of the type ${\overrightarrow{d}}_{t+1-\alpha }=\left(1-\alpha \right){\overrightarrow{d}}_{t+1}+\alpha {\overrightarrow{d}}_t$ is used to evaluate any vector in (12) at an auxiliary instant using its current value ${\overrightarrow{d}}_{t+1}$ and its past instant value ${\overrightarrow{d}}_t$. The parameter $\alpha$ may represent ${\alpha }_f$ or ${\alpha }_m$. Applying this combination in the equation of motion results

$$\begin{array}{l}\left(1-{\alpha }_f\right){\overrightarrow{F}}_{t+1}^{\mathrm{int}}+{\alpha }_f{\overrightarrow{F}}_t^{\mathrm{int}}-\left(1-{\alpha }_f\right){\overrightarrow{F}}_{t+1}-{\alpha }_f{\overrightarrow{F}}_t+\left(1-{\alpha }_m\right)M\cdot {\ddot{\overrightarrow{Y}}}_{t+1}+{\alpha }_{m}M\cdot {\ddot{\overrightarrow{Y}}}_t+\\ +\left(1-{\alpha }_f\right)D\cdot {\dot{\overrightarrow{Y}}}_{t+1}+{\alpha }_{f}D\cdot {\dot{\overrightarrow{Y}}}_t+\left(1-{\alpha }_f\right){\overrightarrow{F}}_{t+1}^{\text{con}}+{\alpha }_f{\overrightarrow{F}}_t^{\text{con}}=\overrightarrow{0}\end{array}$$.(13)

Nodal parameters velocity and acceleration vectors are written as

$${\overrightarrow{Y}}_{t+1}={\overrightarrow{Y}}_t+\Delta t{\dot{\overrightarrow{Y}}}_t+\Delta t^2\left[\left(\frac{1}{2}-\beta \right){\ddot{\overrightarrow{Y}}}_t+\beta {\ddot{\overrightarrow{Y}}}_{t+1}\right]$$ (14)

and

$${\dot{\overrightarrow{Y}}}_{t+1}={\dot{\overrightarrow{Y}}}_t+\Delta t\left(1-\gamma \right){\ddot{\overrightarrow{Y}}}_t+\gamma \Delta t{\ddot{\overrightarrow{Y}}}_{t+1}$$ (15)

in which the Newmark approximations with parameters $\beta$ and $\gamma$ are employed. Using (14) and (15) in Eq. (13) results the equations of motion for discrete time instants

$$\left(1-{\alpha }_f\right)\left({\overrightarrow{F}}_{t+1}^{\mathrm{int}}+{\overrightarrow{F}}_{t+1}^{\text{con}}-{\overrightarrow{F}}_{t+1}\right)+\left[\frac{\left(1-{\alpha }_m\right)}{\beta \Delta t^2}M+\frac{\left(1-{\alpha }_f\right)\gamma }{\beta \Delta t}D\right]{\overrightarrow{Y}}_{t+1}+{\overrightarrow{P}}_t=\overrightarrow{0}$$ (16)

with the vector that groups the terms of the previous time instant given by

$${\overrightarrow{P}}_t={\alpha }_f\left({\overrightarrow{F}}_t^{\mathrm{int}}+{\overrightarrow{F}}_t^{\text{con}}-{\overrightarrow{F}}_t\right)+{\alpha }_{m}M\cdot {\ddot{\overrightarrow{Y}}}_t-\left(1-{\alpha }_m\right)M\cdot {\overrightarrow{T}}_t+{\alpha }_{f}D\cdot {\dot{\overrightarrow{Y}}}_t+\left(1-{\alpha }_f\right)D\cdot {\overrightarrow{R}}_t$$ (17)

with

$${\overrightarrow{T}}_t=\left(\frac{1}{2\beta }-1\right){\ddot{\overrightarrow{Y}}}_t+\frac{{\dot{\overrightarrow{Y}}}_t}{\beta \Delta t}+\frac{{\overrightarrow{Y}}_t}{\beta \Delta t^2}$$ (18)

and

$${\overrightarrow{R}}_t=\left(1-\frac{\gamma }{2\beta }\right)\Delta t{\ddot{\overrightarrow{Y}}}_t\text{+}\left(1-\frac{\gamma }{\beta }\right){\dot{\overrightarrow{Y}}}_t-\frac{\gamma {\overrightarrow{Y}}_t}{\beta \Delta t}$$.(19)

At the end of the time step the velocity and acceleration must be calculated from the approximations (14) and (15), or, using the previously defined vectors, they can be written as

$${\dot{\overrightarrow{Y}}}_{t+1}=\frac{\gamma {\overrightarrow{Y}}_{t+1}}{\beta \Delta t}-{\overrightarrow{R}}_t\text{and}{\ddot{\overrightarrow{Y}}}_{t+1}=\frac{{\overrightarrow{Y}}_{t+1}}{\beta \Delta t^2}-{\overrightarrow{T}}_t$$.(20)

which is best suited for computational calculations as involves less operations.

To start the time march, the initial acceleration is evaluated from the initial configuration and velocity of the system by Eq. (11) which writes as

$${\ddot{\overrightarrow{Y}}}_0=M^{-1}\cdot ({\overrightarrow{F}}_0-{\overrightarrow{F}}_0^{\mathrm{int}}-{\overrightarrow{F}}_0^{\text{con}}-D\cdot {\dot{\overrightarrow{Y}}}_0)$$.(21)

Relating the high frequency region spectral radius ${\rho }_{\infty }\in \left[\mathrm{0,1}\right]$ to the generalized-α method parameters is possible to obtain an optimal situation of minimizing numerical dissipation of the fundamental frequencies and maximizing dissipation on the high frequencies. The α-parameters are obtained by (Chung and Hulbert, 1993)

$${\alpha }_m=\frac{2{\rho }_{\infty }-1}{{\rho }_{\infty }+1}\text{and}{\alpha }_f=\frac{{\rho }_{\infty }}{{\rho }_{\infty }+1}$$,(22)

and the other parameters are given by

$$\gamma =\frac{1}{2}-{\alpha }_m+{\alpha }_f\text{and}\beta =\frac{1}{4}{\left(1-{\alpha }_m+{\alpha }_f\right)}^2$$.(23)

For linear dynamics, respecting Eq. (22) and (23) the method is second order accurate (Chung and Hulbert, 1993). Considering the presence of constraints, the relation ${\alpha }_m<{\alpha }_f\le 1/2$ must also be respected to achieve unconditional stability (proven for the linear case in Géradin and Cardona (2001)). As a consequence of ${\alpha }_m<{\alpha }_f$, the situation without numerical dissipation, that is ${\rho }_{\infty }=1$, cannot be achieved. From the authors experience and the examples presented next, in most situations ${\rho }_{\infty }=0.9$ presents good results with small or negligible numerical dissipation even in nonlinear cases.

4 NONLINEAR EQUATION SOLUTION PROCEDURE

As the equation of motion is nonlinear, the Newton-Raphson method is employed to obtain the current value of the nodal parameters ${\overrightarrow{Y}}_{t+1}$. In order to do so, the equilibrium equation (16), already discretized in time, is rewritten defining the mechanical unbalance vector ${\overrightarrow{g}}_{t+1}$ (numerical residuum) at the current time as

$${\overrightarrow{g}}_{t+1}=\left(1-{\alpha }_f\right)\left({\overrightarrow{F}}_{t+1}^{\mathrm{int}}+{\overrightarrow{F}}_{t+1}^{\text{con}}-{\overrightarrow{F}}_{t+1}\right)+\left[\frac{\left(1-{\alpha }_m\right)}{\beta \Delta t^2}M+\frac{\left(1-{\alpha }_f\right)\gamma }{\beta \Delta t}D\right]{\overrightarrow{Y}}_{t+1}+{\overrightarrow{P}}_t=\overrightarrow{0}$$.(24)

Applying a usual first order Taylor expansion one achieves the following linear system of equations for a trial solution vector ${\left\{\begin{array}{cc}{\overrightarrow{Y}}_{t+1}^0& {\overrightarrow{\lambda }}_{t+1}^0\end{array}\right\}}^t$

$$\left\{\begin{array}{c}\Delta {\overrightarrow{Y}}_{t+1}\\ \Delta {\overrightarrow{\lambda }}_{t+1}\end{array}\right\}=-{\left(H_{t+1}\right)}^{-1}\cdot \overrightarrow{g}\left({\overrightarrow{Y}}_{t+1}^0,{\overrightarrow{\lambda }}_{t+1}^0\right)$$ (25)

where $H_{t+1}=\nabla {\overrightarrow{g}}_{t+1}$ is the Hessian matrix, or tangent operator of the Newton method.

Solving (25), a correction ${\left\{\begin{array}{cc}\Delta {\overrightarrow{Y}}_{t+1}^{}& \Delta {\overrightarrow{\lambda }}_{t+1}^{}\end{array}\right\}}^t$ is added to the trial solution. This procedure is repeated until the relative position norm - calculated as $\left|\right|\Delta \overrightarrow{Y}\left|\right|/\left|\right|\overrightarrow{X}\left|\right|$, where $\overrightarrow{X}$ is the initial nodal position vector of the body - is smaller than a given tolerance.

The Hessian matrix can be written in a compact way, respecting the organization of the degrees of freedom of the system, as

$$H_{t+1}=\left(1-{\alpha }_f\right)\left(\nabla {\overrightarrow{F}}_{t+1}^{\mathrm{int}}+\nabla {\overrightarrow{F}}_{t+1}^{\text{con}}\right)+\frac{\left(1-{\alpha }_m\right)}{\beta \Delta t^2}M+\frac{\left(1-{\alpha }_f\right)\gamma }{\beta \Delta t}D$$ (26)

In which $\nabla {\overrightarrow{F}}_{t+1}^{\mathrm{int}}$ represents the elastic part of the tangent (or Hessian) matrix associated with the strain energy potential of the finite elements and $\nabla {\overrightarrow{F}}_{t+1}^{\text{con}}$ represents the part associated with imposed kinematic constraints.

The elastic part can be developed to obtain

$$\nabla {\overrightarrow{F}}_{t+1}^{\mathrm{int}}={\left.\frac{\partial {\overrightarrow{F}}^{\mathrm{int}}}{\partial \overrightarrow{Y}}\right|}_{t+1}={\displaystyle {\int }_{V_0}^{}\left({\left.\frac{\partial E}{\partial \overrightarrow{Y}}\right|}_{t+1}:\frac{{\partial }^{2}u}{\partial E\otimes \partial E}:{\left.\frac{\partial E}{\partial \overrightarrow{Y}}\right|}_{t+1}+S_{t+1}:{\left.\frac{{\partial }^{2}E}{\partial \overrightarrow{Y}\otimes \partial \overrightarrow{Y}}\right|}_{t+1}\right)d{V}_0}$$ (27)

in which $\otimes$ represents tensor product and $ℭ={\partial }^{2}u/\partial E\otimes \partial E$ is the fourth-order elastic constitutive tensor. For the adopted isotropic Saint-Venant-Kirchhoff model, the specific strain energy for the plane strain state is given by

$$u=\frac{E}{(1+\nu )}\left[\frac{(1-\nu )}{2(1-2\nu )}(E_{11}^2+E_{22}^2)+\frac{\nu }{(1-2\nu )}{E}_{11}{E}_{22}+\frac{1}{2}{E}_{12}{E}_{21}+\frac{1}{4}(E_{12}^2+E_{21}^2)\right]$$ (28)

and for the plane stress state by

$$u=\frac{E}{2(1-{\nu }^2)}(E_{11}^2+E_{22}^2)+\frac{\nu E}{(1-{\nu }^2)}{E}_{11}{E}_{22}+\frac{E}{2(1+\nu )}(E_{12}^2+E_{21}^2)$$ (29)

with $E$ and $\nu$ being the material parameters that represent Young modulus and Poisson ratio in small strains, respectively. Those constitutive relations are employed for the 2D solid. For the plane frame finite element, the plane stress specific strain energy is reduced to avoid volumetric locking (Siqueira and Coda, 2017), and is given by

$$u=\frac{E}{2}(E_{11}^2+E_{22}^2)+G(E_{12}^2+E_{21}^2)$$ (30)

with $G=E/[2(1+\nu )]$ reproducing the shear modulus for small strains.

The Green strain derivatives in Eq. (27) depend on the parameters chosen to define a particular finite element kinematics (as presented in the next section), see Eqs. (35), (40) and (47).

Using Eq. (9), the part of the Hessian that refers to the constraint forces can be derived as

$$\nabla {\overrightarrow{F}}_{t+1}^{\text{con}}={\left.\frac{\partial {\overrightarrow{F}}^{\text{con}}}{\partial \left\{\overrightarrow{Y},\overrightarrow{\lambda }\right\}}\right|}_{t+1}={\left[\begin{array}{cc}\overrightarrow{\lambda }\cdot \nabla \left(\nabla \overrightarrow{c}\right)& \nabla \overrightarrow{c}\\ {\left(\nabla \overrightarrow{c}\right)}^t& 0\end{array}\right]}_{t+1}$$ (31)

where $\nabla \left(\nabla \overrightarrow{c}\right)$ is a third-order tensor which can be understood as the set that groups the Hessian matrices of each constraint equation $c_i$ and $0$ is a null matrix. The derivatives shown in the previous equation are achieved after establishing the kinematic constraint equations for the sliding connections in Section 6.

5 POSITIONAL APPROACH FOR THE FINITE ELEMENTS

To compute the strain field, the motion of the material points composing the body from their initial configuration to their current position must be defined. Following Bonet et al. (2000) and Coda and Paccola (2007), in the positional approach of the FEM the deformation function (or motion) is obtained from the composition of two mappings based on a dimensionless reference domain as

$$\overrightarrow{f}={\overrightarrow{f}}^1\circ {({\overrightarrow{f}}^0)}^{-1}$$ (32)

where ${\overrightarrow{f}}^1$ represents the mapping from the dimensionless space to the current configuration of the body and ${\overrightarrow{f}}^0$ to the initial (reference) configuration. As the deformation function expression does not need to be explicitly known to define the strain measure, but only its gradient, one can write from Eq. (32)

$$A=Grad(\overrightarrow{f})=A^1\cdot {(A^0)}^{-1}$$ (33)

in which $A^1$ and $A^0$ are numerical values calculated at the integration points. It is noteworthy that $A^1$ have to be calculated at each iteration of the solution procedure, but $A^0$ can be calculated only once and stored, as the reference configuration does not change in a total Lagrangian description. Both tensors are obtained from the derivatives of the respective mappings regarding the dimensionless variables $\overrightarrow{\xi }$ as

$$A^0=\frac{\partial {\overrightarrow{f}}^0}{\partial \overrightarrow{\xi }}\text{and}{A}^1=\frac{\partial {\overrightarrow{f}}^1}{\partial \overrightarrow{\xi }}$$.(34)

Knowing the deformation gradient, the objective Green-Lagrange strain can be written as (Ogden, 1984)

$$E=\frac{1}{2}(C-I)=\frac{1}{2}(A^t\cdot A-I)$$ (35)

where $C=A^t\cdot A$ is the right Cauchy-Green stretch tensor.

Replacing the gradient of the motion given in Eq. (33) into Eq. (35), the Green strain derivative component - present in the internal force vector, Eq. (4), and in its gradient, Eq. (27) - can be developed as

$$\frac{\partial E}{\partial \overrightarrow{Y}}=\frac{1}{2}\left[{\left(A^0\right)}^{-T}\cdot {\left(\frac{\partial A^1}{\partial \overrightarrow{Y}}\right)}^T\cdot A^1\cdot {\left(A^0\right)}^{-1}+{\left(A^0\right)}^{-T}\cdot {\left(A^1\right)}^T\cdot \frac{\partial A^1}{\partial \overrightarrow{Y}}\cdot {\left(A^0\right)}^T\right]$$ (36)

in which ${\left(A^0\right)}^{-T}$ represents the inverse transpose of $A^0$. From Eq. (36), the second derivative of the Green strain, necessary for the solution procedure, Eq. (27), results in

$$\frac{{\partial }^{2}E}{\partial \overrightarrow{Y}\otimes \partial \overrightarrow{Y}}=\frac{1}{2}\left(M+M^T+N+N^T\right)$$ (37)

with

$$M={\left(A^0\right)}^{-T}\cdot {\left(A^1\right)}^T\cdot \frac{{\partial }^2{A}^1}{\partial \overrightarrow{Y}\otimes \partial \overrightarrow{Y}}\cdot {\left(A^0\right)}^{-1}$$ (38)

and

$$N={\left(A^0\right)}^{-T}\cdot \frac{\partial {\left(A^1\right)}^T}{\partial \overrightarrow{Y}}\cdot \frac{\partial A^1}{\partial \overrightarrow{Y}}\cdot {\left(A^0\right)}^{-1}$$.(39)

Next, in this section, the frame and 2D solid finite element kinematics are described presenting their particular mappings, thus, allowing the development of their particular derivatives detached in Eqs. (36) and (37).

5.1 2D solid element kinematics

To define the current configuration of the 2D solid it is necessary to map its current domain. This can be done by using the element nodal positions to write the following mapping

$$f_i^1=y_i={\varphi }_{\mathcal{l}}({\xi }_1,{\xi }_2)Y_i^{\mathcal{l}}$$ (40)

in which $Y_i^{\mathcal{l}}$ are the current coordinates of a node $\mathcal{l}$ for directions $i=\mathrm{1,2}$ and ${\varphi }_{\mathcal{l}}({\xi }_1,{\xi }_2)$ are the shape functions of the adopted triangular element, being ${\xi }_1$ and ${\xi }_2$ its dimensionless variables. Fig. 1 depicts the mapping to the current configuration.

Fig. 1
Plane solid mapping to the current configuration.

Known (40) one can develop the derivative of the tensor $A^1$ in Eq. (36) as

$$\frac{\partial A_{ij}^1}{\partial Y_{\alpha }^{\beta }}={\varphi }_{\beta ,j}{\delta }_{\alpha i}(i,j=\mathrm{1,2})$$ (41)

where $\alpha =\mathrm{1,2}$ refers to node $\beta$ coordinates, ${\varphi }_{\beta ,j}=\partial {\varphi }_{\beta }/\partial {\xi }_j$ and ${\delta }_{\alpha i}$ is the Kronecker delta. As the second derivative of $A^1$ regarding the nodal parameters is null for this element, expression (37) simplifies to

$$\frac{{\partial }^{2}E}{\partial \overrightarrow{Y}\otimes \partial \overrightarrow{Y}}=\frac{1}{2}\left(N+N^T\right)$$ (42)

for which only the already obtained first derivative of $A^1$ is necessary as can be observed in Eq. (39).

In order to complete the description of the elastic part of the equilibrium, Eq. (4), and the solution procedure, Eq. (27), it is necessary to write the second Piola-Kirchhoff stress. Employing the energy conjugacy $S=\partial u/\partial E$, for plane strain, Eq. (28), results

$$S=\left[\begin{array}{cc}\left(\lambda +2G\right)E_{11}+\lambda E_{22}& 2G{E}_{12}\\ 2G{E}_{21}& \left(\lambda +2G\right)E_{22}+\lambda E_{11}\end{array}\right]$$ (43)

with $\lambda =E\nu /[(1+\nu )(1-2\nu )]$. Although not necessary for the solution procedure, in a plane strain condition, the following stress components can also be calculated

$$\begin{array}{l}{S}_{i3}=S_{3i}=0(i=\mathrm{1,2})\\ S_{33}=\nu (S_{11}+S_{22})\end{array}$$ .(44)

For plane stress, using Eq. (29), the second Piola-Kirchhoff stress tensor results

$$S=\left[\begin{array}{cc}\frac{E}{1-{\nu }^2}\left(E_{11}+\nu E_{22}\right)& 2G{E}_{12}\\ 2G{E}_{21}& \frac{E}{1-{\nu }^2}\left(E_{22}+\nu E_{11}\right)\end{array}\right]$$ (45)

in this case, the strain components following the out-of-plane axis are achieved by

$$\begin{array}{l}{E}_{i3}=E_{3i}=0(i=\mathrm{1,2})\\ E_{33}=-\frac{\nu }{E}(S_{11}+S_{22})\end{array}$$ .(46)

To complete the details of the 2D solid element, it must be mentioned that the mapping for the initial configuration can be obtained similarly to Eq. (40) adopting the initial node’s positions in place of the current ones. As triangular elements are employed Hammer quadrature is adopted (Zienkiewicz et al., 2013) for the numerical integrations necessary to calculate the internal force vector, its gradient and mass matrix.

5.2 Plane frame element kinematics

The description of the frame finite element current configuration depends on its adopted kinematics. Similarly to the 2D solid finite element, positions instead of displacements are adopted as nodal parameters. Angles are used to describe the cross-section orientation allowing cross-sectional spin independent from the frame reference line. As a result, a Timoshenko-Reissner kinematics takes place and the frame element has three degrees of freedom per node. The current configuration mapping can be written as

$$\begin{array}{l}{f}_1^1=y_1={\varphi }_{\mathcal{l}}(\xi )Y_1^{\mathcal{l}}+\frac{h_0}{2}\eta \cos [{\varphi }_{\mathcal{l}}(\xi ){\theta }_{\mathcal{l}}]\\ f_2^1=y_2={\varphi }_{\mathcal{l}}(\xi )Y_2^{\mathcal{l}}+\frac{h_0}{2}\eta \sin [{\varphi }_{\mathcal{l}}(\xi ){\theta }_{\mathcal{l}}]\end{array}$$ (47)

where $Y_i^{\mathcal{l}}$ are the current coordinates of a node $\mathcal{l}$ for the directions $i=\mathrm{1,2}$, ${\theta }_{\mathcal{l}}$ is the nodal cross-section angle and ${\varphi }_{\mathcal{l}}(\xi )$ are the shape functions for the reference line of the frame, being $\xi$ its dimensionless variable. Fig. 2 depicts this mapping identifying also the dimensionless variable $\eta$ employed to map material point positions along the transverse direction, respecting the initial height $h_0$ of the frame. In Eq. (47), the first summand of each equation represents the map to a point in the reference line and the second one maps the frame height from this reference line, using the same shape functions to interpolate nodal angles.

Fig. 2
Plane frame mapping to the current configuration.

From the mapping to the current configuration, Eq. (47), the first derivative of the tensor $A^1$ in Eq. (36) can be developed as

$$\frac{\partial A^1}{\partial Y_1^{\beta }}=\left[\begin{array}{cc}{\varphi }_{\beta ,\xi }& 0\\ 0& 0\end{array}\right]\text{and}\frac{\partial A^1}{\partial Y_2^{\beta }}=\left[\begin{array}{cc}0& 0\\ {\varphi }_{\beta ,\xi }& 0\end{array}\right]$$ (48)

regarding the positions of any node $\beta$, and as

$$\frac{\partial A^1}{\partial Y_3^{\beta }}=\frac{\partial A^1}{\partial {\theta }_{\beta }}=\left[\begin{array}{cc}-\frac{h_0}{2}\eta \cos ({\varphi }_{\mathcal{l}}{\theta }_{\mathcal{l}}){\varphi }_{\beta }({\varphi }_{k,\xi }{\theta }_k)-\frac{h_0}{2}\eta \sin ({\varphi }_{\mathcal{l}}{\theta }_{\mathcal{l}}){\varphi }_{\beta ,\xi }& -\frac{h_0}{2}\sin ({\varphi }_{\mathcal{l}}{\theta }_{\mathcal{l}}){\varphi }_{\beta }\\ -\frac{h_0}{2}\eta \cos ({\varphi }_{\mathcal{l}}{\theta }_{\mathcal{l}}){\varphi }_{\beta }({\varphi }_{k,\xi }{\theta }_k)+\frac{h_0}{2}\eta \sin ({\varphi }_{\mathcal{l}}{\theta }_{\mathcal{l}}){\varphi }_{\beta ,\xi }& \frac{h_0}{2}\cos ({\varphi }_{\mathcal{l}}{\theta }_{\mathcal{l}}){\varphi }_{\beta }\end{array}\right]$$ (49)

for the nodal angles. The second derivative of $A^1$, necessary in (38), is non-null only for the angles, resulting, for any nodes $\beta$ and $\zeta$ in

$$\begin{array}{l}\frac{{\partial }^2{A}^1}{\partial Y_3^{\beta }\partial Y_3^{\zeta }}=\frac{{\partial }^2{A}^1}{\partial {\theta }_{\beta }^{}\partial {\theta }_{\zeta }^{}}=\\ =\left[\begin{array}{cc}\frac{h_0}{2}\eta \left\{\sin ({\varphi }_{\mathcal{l}}{\theta }_{\mathcal{l}}){\varphi }_{\beta }{\varphi }_{\zeta }({\varphi }_{k,\xi }{\theta }_k)-\cos ({\varphi }_{\mathcal{l}}{\theta }_{\mathcal{l}})\left[{\varphi }_{\beta }{\varphi }_{\zeta ,\xi }+{\varphi }_{\beta ,\xi }{\varphi }_{\zeta }\right]\right\}& -\frac{h_0}{2}\cos ({\varphi }_{\mathcal{l}}{\theta }_{\mathcal{l}}){\varphi }_{\beta }{\varphi }_{\zeta }\\ -\frac{h_0}{2}\eta \left\{\cos ({\varphi }_{\mathcal{l}}{\theta }_{\mathcal{l}}){\varphi }_{\beta }{\varphi }_{\zeta }({\varphi }_{k,\xi }{\theta }_k)+\sin ({\varphi }_{\mathcal{l}}{\theta }_{\mathcal{l}})\left[{\varphi }_{\beta }{\varphi }_{\zeta ,\xi }+{\varphi }_{\beta ,\xi }{\varphi }_{\zeta }\right]\right\}& -\frac{h_0}{2}\sin ({\varphi }_{\mathcal{l}}{\theta }_{\mathcal{l}}){\varphi }_{\beta }{\varphi }_{\zeta }\end{array}\right]\end{array}$$.(50)

As the height of the frame element is not free to change from its initial value when deformation occurs, to avoid volumetric locking, the constitutive equation is relaxed in order to exclude transverse expansions, as shown in Eq. (30). Using the strain energy conjugacy and the constitutive relation, the Piola-Kirchhoff stress tensor results in

$$S=\frac{\partial u}{\partial E}=\left[\begin{array}{cc}E{E}_{11}& 2G{E}_{12}\\ 2G{E}_{21}& E{E}_{22}\end{array}\right]$$.(51)

The mapping of the initial configuration is analogous to Eq. (47) when adopting the initial values of the nodal parameters. To evaluate the necessary integrals, numerical integration with Gauss quadrature is used for each dimensionless variable direction. More details on this finite element can be obtained in Siqueira and Coda (2017), for instance.

5.3 Mass and damping matrices

To obtain the 2D solid element mass matrix one starts from the kinetic energy term present in Eq. (4) and uses the mapping described by Eq. (40), resulting

$${\displaystyle {\int }_{V_0}^{}{\rho }_0{\ddot{y}}_i\delta y_{i}d{V}_0}={\displaystyle {\int }_{V_0}^{}{\rho }_0{\varphi }_{\mathcal{l}}{\ddot{Y}}_i^{\mathcal{l}}{\varphi }_m\delta {\ddot{Y}}_i^{m}d{V}_0}={\displaystyle {\int }_{V_0}^{}{\rho }_0{\varphi }_{\mathcal{l}}{\varphi }_{m}d{V}_0}{\ddot{Y}}_i^{\mathcal{l}}\delta {\ddot{Y}}_i^m(i=\mathrm{1,2})$$ (52)

where the last integral can be identified as the mass matrix, which, in closed form is

$$M\cdot \ddot{\overrightarrow{Y}}={\displaystyle {\int }_{V_0}^{}{\rho }_0\overrightarrow{\varphi }\otimes \overrightarrow{\varphi }d{V}_0}\cdot \ddot{\overrightarrow{Y}}$$.(53)

As stated in previous sections, the resulting matrix is constant as the use of position degrees of freedom consistently relates to the translational inertia of the body.

The mass matrix for the frame element is obtained analogously by starting with the kinetic energy term in (4) but using the frame mapping given by Eq. (47). Considering only the frame reference line, the mass matrix results in the expression given by Eq. (53). In general, for slender bars, as the ones analyzed in this study, rotational inertia (related to angle parameters) is not as important as the translational one (related to the translational parameters). Therefore, this term is neglected here, resulting in constant mass matrix.

To complete the element dynamical information, the damping matrix is assumed to represent the Rayleigh-type damping as

$$D=c_{m}M+c_{k}K$$ (54)

were $c_m$ and $c_k$ are, respectively, the damping coefficients for the mass and stiffness matrices which can be obtained by (Paultre, 2011)

$$\left\{\begin{array}{c}{c}_m\\ c_k\end{array}\right\}=2\frac{{\omega }_i{\omega }_j}{{\omega }_j^2-{\omega }_i^2}\left[\begin{array}{cc}{\omega }_j& -{\omega }_i\\ -1/{\omega }_j& 1/{\omega }_j\end{array}\right]\left\{\begin{array}{c}{\xi }_i\\ {\xi }_j\end{array}\right\}$$ (55)

in which ${\omega }_i$ and ${\omega }_j$ are angular frequencies and ${\xi }_i$ and ${\xi }_j$ are the damping ratios for two distinct modes of vibration i and j, respectively. Paultre (2011) recommends to select ${\omega }_i$ as the first mode of the system and ${\omega }_j$ as the higher frequency that contributes to the response. Also, due to little knowledge of the damping ratios it is usually assumed that ${\xi }_i={\xi }_j$.

Rayleigh-type damping is frequently employed in dynamical analysis (Papageorgiou and Gantes, 2010; Petyt, 2010), particularly common in the linear case were both mass and stiffness of the body are constant. In the geometrically nonlinear framework of this study the mass was already shown to be constant, however the stiffness, understood as $\nabla {\overrightarrow{F}}_{}^{\mathrm{int}}$, Eq. (27), varies with the deformation of the body.

Considering that viscous damping is a simplification of many dissipative phenomena occurring in the structure, we assume the stiffness of the body in Eq. (54) as the initial Hessian matrix. Thus, the damping matrix is also constant and can be stored in the beginning of computations.

From this assumption follows that a modal analysis for the initial configuration of the body (which uses $K=\nabla {\overrightarrow{F}}_0^{\mathrm{int}}$) is performed before starting the time integration to obtain the values of the damping coefficients in (54), by employing Eq. (55), and also to help the choice of the time step in practical analyses.

6 SLIDING CONNECTIONS

Equations (7) and (8) describe how to introduce constraints in equilibrium equations. In this section the general constraint is particularized to sliding connections including roughness. The referred connection represents the bilateral contact of two surfaces, curves for the plane case, able to slide among each other while retaining contact and not overlapping each other. A distance between contacting points may exist, which is the roughness profile representing surface asperities.

In order to write the constraint equation for the present case, we consider that sliding points $\widehat{P}$ belonging to a finite element are constrained to move in a path, defined by an arc-length parameter $s$, with corresponding contact points $\overline{P}$. The path may possess an arbitrary roughness profile defined by a vector $\overrightarrow{r}(s)$ as depicted in Fig. 3. The path (by means of its nodes) is associated to a finite element (called path element) making connection between elements, or may have its degrees of freedom constrained, representing rigid surfaces.

Fig. 3
Sliding connections for a plane solid element considering arbitrary roughness profile in its path.

The sliding contact condition must ensure that $\widehat{P}=\overline{P}+\overrightarrow{r}$ at all times. As the geometry of the bodies is discretized by finite elements, this condition can be expressed directly using their nodal positions, therefore the constraint equations, see Eqs. (7) and (8), can be written as

$$c_i^k={\widehat{Y}}_i^{{\widehat{P}}_k}-{\varphi }_{\mathcal{l}}({\xi }_{{\overline{P}}_k}){\overline{Y}}_i^{\mathcal{l}}-r_i(s_{{\overline{P}}_k})=0_i$$ (56)

for any $k$ pair of points in contact (${\widehat{P}}_k$ and ${\overline{P}}_k$) for both directions $i=\mathrm{1,2}$. In (56), ${\widehat{Y}}_i^{{\widehat{P}}_k}$ is the sliding node position. Also, to obtain the contact point position in the path, ${\overline{Y}}_i^{{\overline{P}}_k}$, the discretization ${\overline{Y}}_i^{{\overline{P}}_k}={\varphi }_{\mathcal{l}}({\xi }_{{\overline{P}}_k}){\overline{Y}}_i^{\mathcal{l}}$ was used, which consider a set of $\mathcal{l}$ path nodes positions ${\overline{Y}}_i^{\mathcal{l}}$ and their respective shape functions ${\varphi }_{\mathcal{l}}(\xi )$ evaluated at the corresponding contact point ${\xi }_{{\overline{P}}_k}$ in the dimensionless space of the associated path element.

The path discretization, as stated, can be done directly by employing the reference line of frame elements. Thus, when sliding occurs among 2D solids, a coupling strategy can be employed as discussed in a later section. An advantage of using frame elements, instead of a simple geometrical discretization, is allowing the definition of the roughness profile direction along the path. It is done by using the frame cross-section angle (or direction) at the contact point ${\overline{\theta }}_{{\overline{P}}_k}$. Thus, the components of the roughness profile can be written as

$$\begin{array}{l}{r}_1(s_{{\overline{P}}_k})=r_h(s_{{\overline{P}}_k})\cos {\overline{\theta }}_{{\overline{P}}_k}=r_h(s_{{\overline{P}}_k})\cos \left[{\varphi }_{\mathcal{l}}({\xi }_{{\overline{P}}_k}){\overline{\theta }}_{\mathcal{l}}\right]\\ r_2(s_{{\overline{P}}_k})=r_h(s_{{\overline{P}}_k})\sin {\overline{\theta }}_{{\overline{P}}_k}=r_h(s_{{\overline{P}}_k})\sin \left[{\varphi }_{\mathcal{l}}({\xi }_{{\overline{P}}_k}){\overline{\theta }}_{\mathcal{l}}\right]\end{array}$$ (57)

where ${\overline{\theta }}_{\mathcal{l}}$ are the known angles of the path nodes $\mathcal{l}$ and $r_h(s)=‖\overrightarrow{r}(s)‖$ is the ‘height’ of the roughness profile, which can be simply expressed as a function of the curvilinear coordinate $s$, i.e., independent of the path discretization. Moreover, the roughness profile as expressed by (57) accompanies the current position of the deformed connected bodies. It is important to advance that the curvilinear coordinate $s$ can be related to the dimensionless coordinate $\xi$, as it will be clarified next.

From the constraint equations, Eq. (56), it is now possible to write its derivatives, necessary to calculate the constraint force, Eq. (9). The complete set of equations is presented in the ^{Appendix 1}.

As can be noted, curvilinear coordinates are new variables introduced by the presence of each sliding connection. By considering a point of the path to be its origin, the initial value of each contact point is known to be the curvilinear distance to the specified origin. During the motion, the current value of $s_{{\overline{P}}_k}$ is obtained directly in the solution updating procedure. If, for instance, the initial value of the curvilinear coordinate is made null, its current value is found to be the relative displacement of the point (measured along the path curve). However, when a roughness profile is involved, particularly in the presence of several contact points, attention must be paid to specify the intended initial value measured in the path, as the profile height is dependent on it.

6.1 Contact point identification

The identification of the curvilinear variable $s_{{\overline{P}}_k}=s({\xi }_{{\overline{P}}_k})$ is necessary to find the correct position in the roughness profile. As $s_{{\overline{P}}_k}$ is the main variable introduced by the constraint equations, it is found immediately in the solution procedure updating process. However, the nondimensional variable evaluated at the contact point of the active path element, ${\xi }_{{\overline{P}}_k}$, not directly obtained, is critical as it is employed in the expressions of the solution procedure and also used to identify the transition between path elements.

The sliding contact condition for each contact pair, expressed by the constraint equations in (56), can be seen as an overdetermined nonlinear system to find ${\xi }_{{\overline{P}}_k}$, when all other variables are known in a particular iteration of a time step. Therefore, a least square technique (Nocedal and Wright, 1999) is employed to obtain the nondimensional variable.

Expressing the residual ${\mathcal{R}}_i$, function of the only unknown ${\xi }_{{\overline{P}}_k}$, from Eq. (56) as

$${\mathcal{R}}_i({\xi }_{{\overline{P}}_k})=c_i^k={\widehat{Y}}_i^{{\widehat{P}}_k}-{\varphi }_{\mathcal{l}}({\xi }_{{\overline{P}}_k}){\overline{Y}}_i^{\mathcal{l}}-r_i(s_{{\overline{P}}_k})=0_i$$ (58)

where $i=\mathrm{1,2}$, is possible to define the following objective function

$p({\xi }_{{\overline{P}}_k})=\frac{1}{2}{\mathcal{R}}_i({\xi }_{{\overline{P}}_k}){\mathcal{R}}_i({\xi }_{{\overline{P}}_k})\text{or}p({\xi }_{{\overline{P}}_k})=\frac{1}{2}{‖\overrightarrow{\mathcal{R}}({\xi }_{{\overline{P}}_k})‖}^2$.(59)

As the least squares approach seeks the smallest residual, and the objective function is quadratic and positive, the required condition for minimization is $\nabla p({\xi }_{{\overline{P}}_k})=0$ at the solution. Expanding the objective function by a first order Taylor series as

$$p\left({\xi }_{{\overline{P}}_k}\right)\cong p\left({\xi }_{{\overline{P}}_k}^0\right)+\nabla p\left({\xi }_{{\overline{P}}_k}^0\right)\Delta {\xi }_{{\overline{P}}_k}=0$$ (60)

being ${\xi }_{{\overline{P}}_k}^0$ a trial value, and writing

$$\nabla p\left({\xi }_{{\overline{P}}_k}\right)=\nabla p\left({\xi }_{{\overline{P}}_k}^0\right)+{\nabla }^{2}p\left({\xi }_{{\overline{P}}_k}^0\right)\Delta {\xi }_{{\overline{P}}_k}=0$$ (61)

the required minimization condition results in the Newton method to be employed to solve for the increment $\Delta {\xi }_{{\overline{P}}_k}$ as

$$\Delta {\xi }_{{\overline{P}}_k}=-\frac{\nabla p\left({\xi }_{{\overline{P}}_k}^0\right)}{{\nabla }^{2}p\left({\xi }_{{\overline{P}}_k}^0\right)}$$.(62)

The non-dimensional variable is updated by ${\xi }_{{\overline{P}}_k}={\xi }_{{\overline{P}}_k}^0+\Delta {\xi }_{{\overline{P}}_k}$ until the norm $|\Delta {\xi }_{{\overline{P}}_k}/{\xi }_{{\overline{P}}_k}^0|$ satisfies a certain tolerance. In the authors experience, and for the studied cases, usually fast convergence is found, with about 3 iterations for a given tolerance of 1.10^-8, when assuming the starting trial value ${\xi }_{{\overline{P}}_k}^0=0$ or the last known value of the dimensionless variable from a previous iteration.

The derivatives in (62) can be developed as

$$\nabla p\left({\xi }_{{\overline{P}}_k}\right)=\frac{dp}{d{\xi }_{{\overline{P}}_k}}={\mathcal{R}}_i\frac{d{\mathcal{R}}_i}{d{\xi }_{{\overline{P}}_k}}$$ (63)

and

$${\nabla }^{2}p\left({\xi }_{{\overline{P}}_k}\right)=\frac{d^{2}p}{d{\xi }_{{\overline{P}}_k}^2}={\mathcal{R}}_i\frac{d^2{\mathcal{R}}_i}{d{\xi }_{{\overline{P}}_k}^2}+\frac{d{\mathcal{R}}_i}{d{\xi }_{{\overline{P}}_k}}\frac{d{\mathcal{R}}_i}{d{\xi }_{{\overline{P}}_k}}$$ (64)

for which differentiation of the residual, Eq. (58), provides

$$\frac{d{\mathcal{R}}_i}{d{\xi }_{{\overline{P}}_k}}=-{\varphi }_{\mathcal{l},\xi }\left({\xi }_{{\overline{P}}_k}\right)-J_{{\overline{P}}_k}\frac{d{r}_i}{d{s}_{{\overline{P}}_k}}$$ (65)

and

$$\frac{d^2{\mathcal{R}}_i}{d{\xi }_{{\overline{P}}_k}^2}=-{\varphi }_{\mathcal{l},\xi \xi }\left({\xi }_{{\overline{P}}_k}\right)-\frac{{\overline{Y}}_{p,\xi }^{{\overline{P}}_k}{\overline{Y}}_{p,\xi \xi }^{{\overline{P}}_k}}{J_{{\overline{P}}_k}}\frac{d{r}_i}{d{s}_{{\overline{P}}_k}}-J_{{\overline{P}}_k}\frac{d^2{r}_i}{d{s}_{{\overline{P}}_k}^2}$$ (66)

with summation implied on the index $p=\mathrm{1,2}$. The Jacobian $J_{{\overline{P}}_k}$ expression and the roughness profile vector derivatives are given by Eq. (74) to (76), (83) and (84) in the ^{Appendix 1}.

The previous iterative procedure for a given sliding node, when performed for all possible path elements, results in ${\xi }_{{\overline{P}}_k}$ and its associated current path element. The search for the active path element ends when the converged value for ${\xi }_{{\overline{P}}_k}$, obtained for a particular element, respects the nondimensional space bounds $\xi \in \left[-\mathrm{1,1}\right]$. It is interesting to note that, as there is no superposition of path elements, the search is optimized by limiting this check only for the path element which its bounding box includes the sliding node position.

7 FINITE ELEMENT COUPLING STRATEGY AND REMARKS

Coupling between frame and solid elements has a double function, it serves to create the path for sliding connections but it also allows moment transfer among frame and solid elements.

As both types of elements use positions as nodal parameters, coupling occur by kinematical compatibility of this variable, i.e., the direct association of those degrees of freedom at shared nodes, see Fig. 4. Although quadratic elements are described in figure, the formulation allows elements with any approximation order. The extension of the auxiliary elements (that can also be used as path elements) are specified respecting the model dimensions.

Fig. 4
Coupling between frame and solid elements for moment transmission and sliding path creation (contact points farther apart for clarity).

Fig. 4 also shows the sliding connections and their path created from frame elements. In the way the constraint equation (56) is written, the sliding node is characterized only by its nodal position, thus, both nodes belonging to a solid or a frame element can be constrained to move along a chosen path using the same expression. Thus, there is no difference in computational implementation for frame/frame or solid/frame sliding.

It is noteworthy that the positions of sliding nodes do not necessarily have to coincide with the position of path nodes, however sliding nodes are contained in a path element. Also, the path element stiffness and mass may or may not be assigned, i.e., for null values the material properties of the bodies are exclusively dependent on the sliding bodies. The use of non-null properties for path elements represents a physical contact surface that possesses different material properties from the sliding bodies.

8 NUMERICAL APPLICATIONS AND DISCUSSION

Three case studies are presented in this section to assess the proposed sliding connection dynamical framework in terms of mechanical behavior of the structural system and to demonstrate scenarios where the occurrence of sliding is necessary to depict real applications. For the iterative solution in the following simulations the relative position tolerance adopted is 1.10^-8.

8.1 Study of the quick return mechanism behavior and its numerical stability

To demonstrate a mechanical application, a classic quick return mechanism with flexible bodies is analyzed. This mechanism was initially proposed and studied by Bauchau (2000) using frame elements. Here we analyze it with 2D solid elements (Fig. 5). The mechanism is composed by a 1.0 m arm AB able to spin around a fixed support B and connected to an 0.25 m bar NA by a revolute joint. The motion is imposed by a 0.20 m crank RS which turns around the fixed support R with constant angular speed $\Omega$ = 5π rad/s and slides over the arm AB with a cylindrical joint S. This mechanism has a total of 3 revolute joints (at A, B and R) and 2 cylindrical joints (at S and N).

Fig. 5
Initial configuration of the mechanism (dimensions in centimeters).

All three bars have 5.98 cm squared cross sections and density 1724.82 kg/m³. The arm is flexible with $E$ = 47.04 GPa. To simulate the rigid behavior of the remaining bars a 10⁶ times greater Young modulus was adopted. In the reference, as frame elements were employed, lumped mass and mass moment of inertia at the slider S (0.31 kg and 0.27 g.m²) and N (2.50 kg and no inertia) were used. To agree with the specified values, here, the sliders were given the dimensions in Fig. 5, with 5.98 cm thickness, and density of 1046.47 kg/m³ for the slider S and 2603.00 kg/m³ for slider N. Both sliders have the same Young modulus as the arm AB. Null Poisson ratio is employed for the whole mechanism and no structural damping is present.

For the solid analysis, as the joints are not treated as points, sliding occur locally and sliding connections are employed to allow the motion of the different joint types. Thus, the revolute joints are modelled by a circular hole at the ends of the bars (with half the size of the bars width) constrained to move in a circular path. The sliders S and N are composed by rectangular blocks having their parallel sides constrained to move along their path (arm AB sides and rigid supports, respectively), which allows the translational motion of the joint. To complete the sliders, a revolute joint exists, connecting the center of each block to the respective bar, allowing rotation. The paths employed for the joints are highlighted in Fig. 6. The rigid support for slider N is also discretized by path elements with all degrees of freedom restricted. For simulation purposes the rigid support has 2.0 m length. No roughness in any sliding path is considered for this example.

Fig. 6
Quick return mechanism mesh (separate for clarity): frame elements in red and solid plane elements sides in black.

Cubic finite elements were adopted for the mesh with 648 plane elements and 105 path (frame) elements. For the circular parts of the mesh (holes and bar ends) each half circumference is discretized by 6 elements. All path elements are massless and, for simulation purposes, with a 1.0 mm squared cross section having the same material properties as the arm AB, except for the slider S paths along arm AB, for which null Young modulus is adopted (to avoid over-stiffening the arm).

The results obtained for the solid model are compared with Bauchau (2000) and refers to the third cycle of motion, the only one present in the reference. Fig. 7 shows the evolution of the horizontal speed of the center of the slider N and Fig. 8 the arm tip deflection (at center of the revolute joint at A) evaluated orthogonally to a line that follows the arm spin at point B. The results were obtained using the generalized-α method with ρ_∞ = 0.90 and a time increment of 0.1 ms. The time step was estimated to be approximately 1/10 of the third period of vibration of the mechanism at the initial resting position. The first three modes (with periods 21.88, 1.81, and 1.09 ms, respectively) are associated with flexure of the arm AB.

Fig. 7
Time history of the slider N speed along x-axis.

Fig. 8
Time history of the arm AB tip deflection.

Good agreement with the reference response is observed. The return movement of the arm is observable in Fig. 8, (starting at about 1.05 s) including the inversion of the predominant side of deflection during the quick return of the mechanism. Fig. 7 shows the clear direction change of the slider N with the return of the arm. Finally, Fig. 9 depicts the motion for selected instants of the cycle.

Fig. 9
Snapshots for the third cycle of the motion with details of the horizontal displacement (m) of the mechanism and contact forces (N) at the sides of slider S. Time instants selected for: a) start of the third cycle, b) slider N maximum displacement, c) maximum value of the total contact force in slider S faces, d) slider N maximum return speed (zero contact force in slider S) and e) end of the third cycle.

Small differences in tip deflection amplitude and vibration phase are perceptible when comparing solid and frame models. These differences are due to kinematic descriptions and temporal integration strategies. The bar kinematics used in Bauchau (2000) is less elaborated than the solid description adopted here, they also use a quite different variable time steps in a strategy to conserve or dissipate the energy of the system.

To further investigate the differences in the response, the mechanism is also discretized and simulated using only positional frame finite elements. 10 cubic frame elements are used for the arm AB, which serves as the path for the slider S, and 2 elements are used for each rigid bar. Material and geometrical properties are the same employed in the solid model. Slider N and points R and B are modelled as simple supports and only one revolute joint (at A) is necessary in this case. The cylindrical joint at the slider S is represented by the same constraint equations present in Eq. (56). A modal analysis is performed for the frame model confirming the same time step previously used.

Results for the frame model are shown in Fig. 7 and Fig. 8 where the same amplitude is found when compared to the solid model, but larger than the reference response. This shows that the present positional formulation is more flexible than the bar theory employed in Bauchau (2000) as it represents the frame element by a solid-like approach. It is interesting to note that the three responses begin the third cycle of the mechanism motion at different instants and are out of phase from each other. Aside from the time integration aspects, we understand that the differences encountered in phase are related to the model mass. Employing the solid element, the rotational inertia can be accounted in a very precise way, whereas it is neglected in the employed frame model and simplified in the reference frame model.

Regarding time integration, simulations are also performed to access the quality of the integration scheme. Cases with (ρ_∞ ≠ 1.0) and without numerical dissipation (ρ_∞ = 1.0) are tested using the generalized-α and the Newmark method - obtained by making α_f = α_m = 0 in Eq. (12). Géradin and Rixen (2015) give the spectral radius relation to Newmark parameters as $\gamma =1/2+\alpha$ and $\beta ={(1+\alpha )}^2/4$ with $\alpha =(1-{\rho }_{\infty })/(1+{\rho }_{\infty })$.

Fig. 10 shows the time evolution of the total contact force acting on slider S sides (obtained directly from the multipliers values and also depicted for some instants in Fig. 9). As expected, without numerical dissipation, the dynamical problem with multipliers presents instabilities, which hindered the complete reproduction of the mechanism cycles.

Fig. 10
Slider S sides total contact force for different integration schemes and spectral radius.

Further tests show that, for this mechanism, a spectral radius in the range 0.0 - 0.9 for the α method results stable and presents almost identical force responses. Also, the numerical damping does not affect significantly the displacement response (Fig. 11), even for the strongest case of damping (ρ_∞ = 0.0). For the Newmark method, on the other hand, stable responses were only obtained at the cost of artificially damping the response. See, for instance, results for ρ_∞ = 0.9 and 0.8. This indicates that the use of this last integration scheme is less adequate for a general case.

Fig. 11
Arm AB tip deflection for different integration schemes and spectral radius.

Finally, despite small differences, very good agreement between results is observed, which confirms that the mechanism behavior is well represented, particularly when using the solid model. Still, the frame model is also adequate and less expensive in terms of number of degrees of freedom and consequently computation time.

8.2 Sensitivity analysis of doubly bent beams regarding roughness and damping

An interesting example of bifurcation of equilibrium under tensile load, which requires the presence of sliding connections to occur, is employed to test and validate the proposed modelling with the experimental and analytical results from Zaccaria et al. (2011). Fig. 12 depicts the structure with two flexible bars with length L, initially aligned, and the horizontal displacement u applied at one side. The reactive traction force F in the right support is measured. The flexible bars are connected to a rigid bar (initially in the vertical position) by sliding connections covering each bar height (Fig. 13), thus representing prismatic joints.

Fig. 12
Initial and deformed configurations of the beams.

Fig. 13
Mesh detail at the sliding connections.

Each flexible bar is discretized by 200 cubic triangular plane solid elements (4 elements along the height). As the vertical bar is rigid, it is simply modelled by 6 cubic frame elements and is used directly as the path for the sliding connections. Fig. 13 shows details of the mesh at the contact points.

Material and geometry properties are taken from the experiment performed in the reference. Each carbon steel AISI 1095 flexible bar ($E$ = 200 GPa, $\nu$ = 0.30 and ${\rho }_0$ = 7850 kg/m³) has an effective length L = 225 mm, rectangular cross section with 1 mm height and 25 mm width. The rigid bar is realized in the experiments by 2 linear bearings (Rollon Easy Rail SN22-80-500-610) each having, according the manufacturer, 0.70 kg/m, 11 x 22 mm cross section and 580 mm for each slider to move. Thus, to simulate the rigid behavior, a 22 mm squared cross section and Young modulus 10³ times greater than the flexible bars is adopted. Each slider has 0.08 kg distributed along the height of the beams. Plane stress conditions are considered for all elements.

To avoid obtaining the trivial numerical solution (i.e., simple axial elongation of the bars), a perturbation is introduced into the system in the form of an initial inclination of the rigid bar by an angle ${\varphi }_0$. After the system has reached static equilibrium for a given value of ${\varphi }_0$, the horizontal displacement is applied. To demonstrate convergence to the analytical solution, progressively smaller values for the initial angle are used to obtain the tension reaction force trajectory in a static simulation. Fig. 14 presents the achieved equilibrium trajectories and the experimental result from Zaccaria et al. (2011). The numerical response agrees well with the experimental data, particularly considering that the rigid bar is not perfectly vertical but has an initial rotation of about 2°. Moreover, as expected, the numerical solution converges to the analytical one for smaller eccentricities.

Fig. 14
Right support reactive force evolution with the displacement for different initial rotations of the rigid bar (statical response).

Good agreement of the deformed shape of the structure obtained in the experiment and in the numerical simulation for different values of slider rotation is achieved as shown in Fig. 15. The contact force distribution along the height of the right beam is also presented showing an orthogonal distribution to the cross section at the prismatic joint, which is expected as the sliding is tangential and free to occur.

Fig. 15
Static numerical simulation (in blue) superimposed to the experimental results from Zaccaria et al. (2011) for different values of slider rotation and corresponding contact forces distribution (N) along the height of the right beam.

To better represent the experimental response, dynamical simulations are performed considering not only the bodies inertia, but also structural damping of the system and roughness in the sliding path. Only the case with ${\varphi }_0$ = 2° (closest to the experimental response) is simulated. The right support is moved horizontally at 0.25 mm/s. Modal analysis for the initial configuration of the system resulted in the first 3 vibration modes related to the flexure of the beams, with angular frequencies of 23.30, 68.11 and 210.14 rad/s. It provided a time step estimative of 18.0 ms. Also, a very small damping ratio ξ = 5.10^-8 is assumed for the first and third modes to calculate both Rayleigh damping matrix coefficients according to Eq. (55). The generalized-α method with ρ_∞ = 0.9 is employed.

A random roughness profile is adopted to represent sliding imperfections in the linear bearings. Following specialized literature (Vakis et al., 2018; Yang and Lin, 1995), surface roughness may be obtained by a Power Spectrum Density (PSD) as

$$r_h\left(s\right)={\displaystyle \sum _{n=1}^N\sqrt{G_n\Delta \Omega }\cos \left(2\pi {\Omega }_{n}s-{\psi }_n\right)}$$ (67)

with $G_n$ being the PSD function associated to a spatial angular frequency ${\Omega }_n$. The spatial frequency increment is obtained by $\Delta \Omega =({\Omega }_{\max }-{\Omega }_{\min })/N$ where N is the amount of frequencies increments in the range $({\Omega }_{\min },{\Omega }_{\max })$ with ${\Omega }_{\min }$ and ${\Omega }_{\max }$ the lower and upper cut-off frequencies defining the PSD function. Each spatial angular frequency is computed by interpolation as ${\Omega }_n={\Omega }_{\min }+\Delta \Omega (n-1)$. The randomness of the profile is associated to a random phase angle ${\psi }_n$ (uniformly distributed from 0 to 2π). The adopted PSD function is taken from ISO (1995), originally proposed by Dodds and Robson (1973), as

$$G_n=G_0{\left(\frac{{\Omega }_n}{{\Omega }_0}\right)}^{-w}$$ (68)

where $G_0$ is the roughness coefficient, a parameter associated with the amount of irregularities, ${\Omega }_0$ = 1.0 rad/m the reference spatial frequency and $w$ the exponent of the PSD function. The parameters to compute the roughness profile were chosen to result in a small profile height (about 1.0 μm), Fig. 16, as $G_0$ = 8.10^-15 m³, $w$ = 1.0, and N = 20 spatial frequencies distributed from ${\Omega }_{\min }$ = 100.0 cycle/m to ${\Omega }_{\max }$ = 1000.0 cycle/m.

Fig. 16
Roughness profile in the sliding path.

Results for the reactive force along the dynamical equilibrium trajectory are shown in Fig. 17. As structural damping and roughness parameters of the sliding path are not precisely known in this case, the previous values are assumed and some differences are expected. Nevertheless, it should be noted that without the presence of structural damping, or, alternatively, without the roughness profile, the dynamical response does not represent the experimental behavior as the presence of both effects is necessary to describe the local interaction at the joints and the overall response of the bodies.

Fig. 17
Right support reactive force evolution with the displacement (dynamical response), simulated results in black and experimental result from Zaccaria et al. (2011) in blue.

8.3 Analysis of a building under seismic activity with imperfect base isolation device

A five-story building subjected to seismic activity is analyzed considering a base isolation device at each column (Fig.18). The device, commonly known as a friction pendulum system (Saaed et al., 2015), is composed by an articulated slider with two sliding interfaces: one in the upper plate (connected to a column) and another in the lower plate (which is supported by the ground). The employed device geometry is inspired by Mokha et al. (1991) and depicted in Fig. 18 along with its finite element discretization.

Fig. 18
Geometry of the building and base isolator (centimeters) with its mesh: frame elements (including structure, auxiliary element and path) in red and 2D solid elements sides in black).

Cubic solid and frame elements are used to discretize the building and its base isolation devices. The reinforced concrete building ($E$ = 20.0 GPa and ${\rho }_0$ = 2400.0 kg/m³) is modelled by 60 frame elements. The columns have 50 cm side squared cross section. For the beams to characterize representative floor stiffness and mass, the same cross section is employed but with 15168.2 kg/m³ density. Each isolation device has 50 cm thickness, is made of steel ($E$ = 200.0 GPa and ${\rho }_0$ = 7850.0 kg/m³) and discretized by 60 plane elements. Massless frame elements are also employed at both sliding paths (with the same elastic material properties as the device) and to define the connection of the column to the upper plate (rigid auxiliary elements with $E$ = 200.10⁶ GPa) as shown in Fig. 18. These elements cross sections have 50 cm width and 1 cm height. Null Poisson ratio is adopted for the structure.

A real earthquake signal is used for the base motion of the isolation devices. The 1987 event known as Superstition Hills is employed with its vertical and horizontal displacement components as measured in the Imperial Valley Wildlife Liquefaction Array station obtained from PEER (2022) and depicted in Fig. 19. Prior to the earthquake displacement, an initial statical analysis is performed with a uniform vertical load of 34.0 kN/m applied at each floor level. This force is maintained throughout the dynamical analysis.

Fig. 19
Superstition Hills earthquake displacement components from PEER (2022) for the Imperial Valley Wildlife Liquefaction Array station.

To evaluate the influence of small manufacturing imperfections on the isolation device, a small roughness profile in the lower slider path is adopted with its height depicted in Fig. 20. A case without roughness in the path is also analyzed.

Fig. 20
Roughness profile in the lower sliding path of the isolation devices: r_h(s) = 1.10^-4.cos(40.0s), expression units in meters.

In order to define damping, we use the first natural angular frequency of the structure ${\omega }_i$ = 8.56 rad/s and the 11^th mode ${\omega }_j$ = 133.98 rad/s. It is important to mention that the first ten natural frequencies are associated with the bending of the whole structure or of beams and the 11^th is associated with columns stretch. Considering these frequencies, a damping ratio ξ = 5%, usual for concrete structures (Papageorgiou and Gantes, 2010), is adopted for both modes in Eq. (55) to obtain the Rayleigh damping matrix. A time step of 5.0 ms is employed with the generalized-α method (ρ_∞ = 0.9) in the absence of surface imperfections, otherwise a time step of 1.0 ms is necessary to obtain stable reaction forces and acceleration, although displacements are properly represented.

Fig. 21 shows the initial configuration of the building and the instant for which the largest sliding occurs for both cases (with and without imperfections). It is noticeable that, without imperfections, the building deformation is symmetric at the beginning of the analysis and the contact forces distribution is regular during all the base motion excitation. In contrast, when the proposed non-symmetrical roughness is present, the structure found initial static equilibrium with a larger displacement at the left device than the right one. Also, the contact force distribution is less regular than the one achieved for the previous case, presenting inverse orientation at some points. This is due to the employed path roughness which causes the bodies to deform locally in order to remain in contact by the imposed bilateral sliding constraints.

Fig. 21
Horizontal displacements (m) of the building (right isolation device shown on the top and left device on the bottom for each time snapshot) and contact forces (N) in the lower slider path for cases a) without and b) with manufacturing imperfections.

The presence of imperfections clearly reduce the efficiency of the isolation devices by diminishing the slider relative displacement and increasing the base shear (horizontal resultant reaction force) as shown in Fig. 22. Also, the interstory drift (difference in horizontal displacement from consecutive floors) and horizontal acceleration of the top floor are larger when roughness is present, even when compared with the case without sliding devices (Fig. 23). These results show an interesting consequence of the presence of manufacturing imperfections which ultimately result in a worse overall response of the structure.

Fig. 22
Right column base shear force and slider relative displacement.

Fig. 23
Right column interstory drift for selected floors and top floor mid-span acceleration.

Finally, it is worth noting from Fig. 23 that the vertical acceleration on the top floor is not affected by the presence of the isolation devices, which is expected as only the horizontal movement is released. Moreover, in the presence of the isolation device (regardless the presence of roughness) the interstory drift changed from the third to the first floor.

9 CONCLUSIONS

A framework for the dynamical simulation of plane structures with sliding connections undergoing finite deformation is presented. A total Lagrangian finite element formulation is used by employing a position-based approach to discretize the bodies. Lagrange multipliers are used to successfully introduce sliding connections constraints into the mechanical system.

An arbitrary roughness profile for the sliding connection path is directly incorporated in the formulation. As the profile is independent of the finite element discretization, the employed meshes could be coarser than discretizing the roughness by solid or frame elements.

The geometrical nonlinear equations of motion, obtained in the positional finite element formulation, render a constant mass matrix allowing to directly introduce Rayleigh-type structural damping, important when dealing with real applications. This characteristic also allowed the successful use of the generalized-α method to address known stability issues in constrained mechanical systems. Tests carried out in the time marching approach confirms the necessity of appropriately controlling the numerical dissipation of the system to obtain representative responses when using Lagrange multipliers in dynamical systems.

Examples illustrate that the proposed framework is efficient even when handling strong geometrical nonlinear problems, in both static and dynamic cases, in the presence of roughness in the sliding connection. Thus, the proposed framework is a viable alternative for general analysis of general multibody flexible systems.

Appendix 1 - Constraint equation derivatives

The derivatives for the constraint equation of the sliding connections, Eq. (56), are presented. For the positions of the sliding node and the path element, it results, respectively

$$\frac{\partial c_i^k}{\partial {\widehat{Y}}_{\alpha }^{{\widehat{P}}_k}}={\delta }_{\alpha i}$$ (69)

and

$$\frac{\partial c_i^k}{\partial {\overline{Y}}_{\alpha }^{\beta }}=-{\varphi }_{\beta }\left({\xi }_{{\overline{P}}_k}\right){\delta }_{\alpha i}$$ (70)

being $\alpha =\mathrm{1,2}$ the directions and $\beta$ a node in the path element. The derivatives for the path angles are

$$\frac{\partial c_1^k}{\partial {\overline{\theta }}_{\beta }}=r_h\left(s_{{\overline{P}}_k}\right)\sin {\overline{\theta }}_{{\overline{P}}_k}{\varphi }_{\beta }\left({\xi }_{{\overline{P}}_k}\right)$$ (71)

and

$$\frac{\partial c_2^k}{\partial {\overline{\theta }}_{\beta }}=-r_h\left(s_{{\overline{P}}_k}\right)\cos {\overline{\theta }}_{{\overline{P}}_k}{\varphi }_{\beta }\left({\xi }_{{\overline{P}}_k}\right)$$ (72)

with ${\overline{\theta }}_{{\overline{P}}_k}={\varphi }_{\mathcal{l}}({\xi }_{{\overline{P}}_k}){\overline{\theta }}_{\mathcal{l}}$. And for the curvilinear position results

$$\frac{\partial c_i^k}{\partial s_{{\overline{P}}_k}}=-\frac{{\overline{Y}}_{i,\xi }^{{\overline{P}}_k}}{J_{{\overline{P}}_k}}-\frac{\partial r_i}{\partial s_{{\overline{P}}_k}}$$ (73)

with ${\overline{Y}}_{i,\xi }^{{\overline{P}}_k}={\varphi }_{\mathcal{l},\xi }({\xi }_{{\overline{P}}_k}){\overline{Y}}_i^{\mathcal{l}}$ and the Jacobian of the transformation from the real space to the nondimensional space of the path element given by

$$J_{{\overline{P}}_k}=\sqrt{{\left[{\varphi }_{\mathcal{l},\xi }\left({\xi }_{{\overline{P}}_k}\right){\overline{Y}}_1^{\mathcal{l}}\right]}^2+{\left[{\varphi }_{\mathcal{l},\xi }\left({\xi }_{{\overline{P}}_k}\right){\overline{Y}}_2^{\mathcal{l}}\right]}^2}$$.(74)

The derivatives for the roughness profile components in Eq. (73) results

$$\frac{\partial r_1}{\partial s_{{\overline{P}}_k}}={\left.\frac{d{r}_h}{ds}\right|}_{s_{{\overline{P}}_k}}\cos {\overline{\theta }}_{{\overline{P}}_k}-\frac{r_h\left(s_{{\overline{P}}_k}\right)\sin {\overline{\theta }}_{{\overline{P}}_k}{\overline{\theta }}_{{\overline{P}}_k,\xi }}{J_{{\overline{P}}_k}}$$ (75)

and

$$\frac{\partial r_2}{\partial s_{{\overline{P}}_k}}={\left.\frac{d{r}_h}{ds}\right|}_{s_{{\overline{P}}_k}}\sin {\overline{\theta }}_{{\overline{P}}_k}+\frac{r_h\left(s_{{\overline{P}}_k}\right)\cos {\overline{\theta }}_{{\overline{P}}_k}{\overline{\theta }}_{{\overline{P}}_k,\xi }}{J_{{\overline{P}}_k}}$$ (76)

with ${\overline{\theta }}_{{\overline{P}}_k,\xi }={\varphi }_{\mathcal{l},\xi }({\xi }_{{\overline{P}}_k}){\overline{\theta }}_{\mathcal{l}}$. The derivative of the profile height is calculated from $r_h(s)$ expression for each case.

The second derivatives of the constraint equation are also needed for the solution procedure using Newton method as seen in Eq. (31). From Eq. (70) the non-null values results in

$$\frac{{\partial }^2{c}_i^k}{\partial {\overline{Y}}_{\alpha }^{\beta }\partial s_{{\overline{P}}_k}}=-\frac{{\varphi }_{\beta ,\xi }\left({\xi }_{{\overline{P}}_k}\right)}{J_{{\overline{P}}_k}}{\delta }_{\alpha i}$$.(77)

From expressions (71) and (72), differentiating for the path angles one obtain

$$\frac{{\partial }^2{c}_1^k}{\partial {\overline{\theta }}_{\beta }\partial {\overline{\theta }}_{\gamma }}=r_h\left(s_{{\overline{P}}_k}\right)\cos {\overline{\theta }}_{{\overline{P}}_k}{\varphi }_{\beta }\left({\xi }_{{\overline{P}}_k}\right){\varphi }_{\gamma }\left({\xi }_{{\overline{P}}_k}\right)$$ (78)

and

$$\frac{{\partial }^2{c}_2^k}{\partial {\overline{\theta }}_{\beta }\partial {\overline{\theta }}_{\gamma }}=r_h\left(s_{{\overline{P}}_k}\right)\sin {\overline{\theta }}_{{\overline{P}}_k}{\varphi }_{\beta }\left({\xi }_{{\overline{P}}_k}\right){\varphi }_{\gamma }\left({\xi }_{{\overline{P}}_k}\right)$$,(79)

and for the curvilinear position

$$\frac{{\partial }^2{c}_1^k}{\partial {\overline{\theta }}_{\beta }\partial s_{{\overline{P}}_k}}={\left.\frac{d{r}_h}{ds}\right|}_{s_{{\overline{P}}_k}}\sin {\overline{\theta }}_{{\overline{P}}_k}{\varphi }_{\beta }\left({\xi }_{{\overline{P}}_k}\right)+\frac{r_h\left(s_{{\overline{P}}_k}\right)}{J_{{\overline{P}}_k}}\left[\cos {\overline{\theta }}_{{\overline{P}}_k}{\overline{\theta }}_{{\overline{P}}_k,\xi }{\varphi }_{\beta }\left({\xi }_{{\overline{P}}_k}\right)+\sin {\overline{\theta }}_{{\overline{P}}_k}{\varphi }_{\beta ,\xi }\left({\xi }_{{\overline{P}}_k}\right)\right]$$ (80)

and

$$\frac{{\partial }^2{c}_2^k}{\partial {\overline{\theta }}_{\beta }\partial s_{{\overline{P}}_k}}={\left.-\frac{d{r}_h}{ds}\right|}_{s_{{\overline{P}}_k}}\cos {\overline{\theta }}_{{\overline{P}}_k}{\varphi }_{\beta }\left({\xi }_{{\overline{P}}_k}\right)+\frac{r_h\left(s_{{\overline{P}}_k}\right)}{J_{{\overline{P}}_k}}\left[\sin {\overline{\theta }}_{{\overline{P}}_k}{\overline{\theta }}_{{\overline{P}}_k,\xi }{\varphi }_{\beta }\left({\xi }_{{\overline{P}}_k}\right)-\cos {\overline{\theta }}_{{\overline{P}}_k}{\varphi }_{\beta ,\xi }\left({\xi }_{{\overline{P}}_k}\right)\right]$$.(81)

Considering the symmetry of the Hessian matrix (by the Schwarz's theorem), starting from Eq. (73), the second derivative regarding the path element positions results analogous to the one in Eq. (77), and differentiation for the path nodes gives an expression analogous to Eq. (80) and (81).

Finally, the second derivative for the curvilinear position is written as

$$\frac{{\partial }^2{c}_i^k}{\partial s_{{\overline{P}}_k}\partial s_{{\overline{P}}_k}}=\frac{{\overline{Y}}_{i,\xi }^{{\overline{P}}_k}{\overline{Y}}_{p,\xi }^{{\overline{P}}_k}{\overline{Y}}_{p,\xi \xi }^{{\overline{P}}_k}}{J_{{\overline{P}}_k}^4}-\frac{{\overline{Y}}_{i,\xi \xi }^{{\overline{P}}_k}}{J_{{\overline{P}}_k}^2}-\frac{{\partial }^2{r}_i}{\partial s_{{\overline{P}}_k}^2}$$ (82)

with ${\overline{Y}}_{p,\xi \xi }^{{\overline{P}}_k}={\varphi }_{\mathcal{l},\xi \xi }({\xi }_{{\overline{P}}_k}){\overline{Y}}_p^{\mathcal{l}}$, summation on the direction $p=\mathrm{1,2}$ implied, and the second derivatives of the roughness profile height components obtained by

$$\begin{array}{l}\frac{{\partial }^2{r}_1}{\partial s_{{\overline{P}}_k}^2}={\left.\frac{d^2{r}_h}{d{s}^2}\right|}_{s_{{\overline{P}}_k}}\cos {\overline{\theta }}_{{\overline{P}}_k}-2{\left.\frac{d{r}_h}{ds}\right|}_{s_{{\overline{P}}_k}}\frac{\sin {\overline{\theta }}_{{\overline{P}}_k}{\overline{\theta }}_{{\overline{P}}_k,\xi }}{J_{{\overline{P}}_k}}-\frac{r_h\left(s_{{\overline{P}}_k}\right)}{J_{{\overline{P}}_k}^2}\left[\cos {\overline{\theta }}_{{\overline{P}}_k}{\left({\overline{\theta }}_{{\overline{P}}_k,\xi }\right)}^2+\right.\\ \left.+\sin {\overline{\theta }}_{{\overline{P}}_k}\left({\overline{\theta }}_{{\overline{P}}_k,\xi \xi }-\frac{{\overline{\theta }}_{{\overline{P}}_k,\xi }{\overline{Y}}_{p,\xi }^{{\overline{P}}_k}{\overline{Y}}_{p,\xi \xi }^{{\overline{P}}_k}}{J_{{\overline{P}}_k}^2}\right)\right]\end{array}$$ (83)

and

$$\begin{array}{l}\frac{{\partial }^2{r}_2}{\partial s_{{\overline{P}}_k}^2}={\left.\frac{d^2{r}_h}{d{s}^2}\right|}_{s_{{\overline{P}}_k}}\sin {\overline{\theta }}_{{\overline{P}}_k}+2{\left.\frac{d{r}_h}{ds}\right|}_{s_{{\overline{P}}_k}}\frac{\cos {\overline{\theta }}_{{\overline{P}}_k}{\overline{\theta }}_{{\overline{P}}_k,\xi }}{J_{{\overline{P}}_k}}+\frac{r_h\left(s_{{\overline{P}}_k}\right)}{J_{{\overline{P}}_k}^2}\left[-\sin {\overline{\theta }}_{{\overline{P}}_k}{\left({\overline{\theta }}_{{\overline{P}}_k,\xi }\right)}^2+\right.\\ \left.+\cos {\overline{\theta }}_{{\overline{P}}_k}\left({\overline{\theta }}_{{\overline{P}}_k,\xi \xi }-\frac{{\overline{\theta }}_{{\overline{P}}_k,\xi }{\overline{Y}}_{p,\xi }^{{\overline{P}}_k}{\overline{Y}}_{p,\xi \xi }^{{\overline{P}}_k}}{J_{{\overline{P}}_k}^2}\right)\right]\end{array}$$.(84)

• https://doi.org/10.1590/1679-78257295

References

• 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.
• Bonari, J., Paggi, M., Dini, D. (2022) A new finite element paradigm to solve contact problems with roughness. International Journal of Solids and Structures 253.
• 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.
• Cardona, A., Géradin, M. (1989) Time integration of the equations of motion in mechanism analysis. Computers & Structures 33(3), 801-820.
• Cardona, A., Géradin, 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.
• Chung, J., Hulbert, G.M. (1993) A Time Integration Algorithm for Structural Dynamics With Improved Numerical Dissipation: The Generalized-α Method. Journal of Applied Mechanics 60(2), 371.
• Coda, H.B., Paccola, R.R. (2011) A FEM procedure based on positions and unconstrained vectors applied to non-linear dynamic of 3D frames. Finite Elements in Analysis and Design 47(4), 319-333.
• 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., Paccola, R.R. (2007) An alternative positional FEM formulation for geometrically non-linear analysis of shells: Curved triangular isoparametric elements. Computational Mechanics 40(1), 185-200.
• Coda, H.B., Paccola, R.R., Sampaio, M.D.S.M. (2013) Positional description applied to the solution of geometrically non-linear plates and shells. Finite Elements in Analysis and Design 67, 66-75.
• Coda, H.B., Silva, A.P.O., Paccola, R.R. (2020) Alternative active nonlinear total lagrangian truss finite element applied to the analysis of cable nets and long span suspension bridges. Latin American Journal of Solids and Structures 17(3).
• Ding, L., Hao, H., Zhu, X. (2009) Evaluation of dynamic vehicle axle loads on bridges with different surface conditions. Journal of Sound and Vibration 323(3-5), 826-848.
• Dodds, C.J., Robson, J.D. (1973) The description of road surface roughness. Journal of Sound and Vibration 31(2), 175-183.
• Espath, L.F.R., Braun, A.L., Awruch, A.M., Dalcin, L.D. (2015) A NURBS-based finite element model applied to geometrically nonlinear elastodynamics using a corotational approach. International Journal for Numerical Methods in Engineering. 102(13), 1839-1868.
• Gay Neto, A., Wriggers, P. (2020) Master-master frictional contact and applications for beam-shell interaction. Computational Mechanics 66(6), 1213-1235.
• Géradin, M., Cardona, A. (2001) Flexible multibody dynamics: a finite element approach. Chichester: John Wiley & Sons.
• Géradin, M., Rixen, D.J. (2015) Mechanical Vibrations: Theory and Application to Structural Dynamics. 3rd ed. Chichester, UK: John Wiley & Sons.
• Gurtin, M.E., Fried, E., Anand, L. (2010) The Mechanics and Thermodynamics of Continua. New York: Cambridge University Press.
• Hilber, H.M., Hughes, T.J.R., Taylor, R.L. (1977) Improved numerical dissipation for time integration algorithms in structural dynamics. Earthquake Engineering & Structural Dynamics 5(3), 283-292.
• Holzapfel, G.A. (2000) Nonlinear solid mechanics: a continuum approach for engineering. Chichester, UK: John Wiley & Sons.
• ISO (1995) Mechanical Vibration — Road Surface Profiles — Reporting of Measured Data (ISO Standard No. 8608).
• 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.
• Kuhl, D., Ramm, E. (1999) Generalized Energy-Momentum Method for non-linear adaptive shell dynamics. Computer Methods in Applied Mechanics and Engineering 178(3-4).
• Lanczos, C. (1970) The variational principles of mechanics. New York: Dover Publications .
• 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.
• Leyendecker, S., Marsden, J.E., Ortiz, M. (2008) Variational integrators for constrained dynamical systems. ZAMM 88(9), 677-708.
• Li, L., Li, G., Wang, J., Fan, C., Cai, A. (2022) Fretting wear mechanical analysis of double rough surfaces based on energy method. Proceedings of the Institution of Mechanical Engineers, Part J: Journal of Engineering Tribology
• Lyubicheva, A.N., Tsukanov, I.Y. (2022) The influence of 2D periodic surface texture on the partial slip problem for elastic bodies. European Journal of Mechanics, A/Solids 91.
• Mamouri, S., Hammadi, F., Ibrahimbegović, A. (2014) Decaying/conserving implicit scheme and non-linear instability analysis of 2D frame structures. International Journal of Non-Linear Mechanics 67, 144-152.
• Mokha, A., Constantinou, M.C., Reinhorn, A.M., Zayas, V.A. (1991) Experimental study of friction-pendulum isolation system. Journal of Structural Engineering (United States) 117(4), 1201-1217.
• Nocedal, J., Wright, S.J. (1999) Numerical optimization. New York: Springer.
• Ogden, R.W. (1984) Non-linear elastic deformations. Chichester: Ellis Horwood.
• Oliva, J., Goicolea, J.M., Antolín, P., Astiz, M.Á. (2013) Relevance of a complete road surface description in vehicle-bridge interaction dynamics. Engineering Structures 56, 466-476.
• Papageorgiou, A. V., Gantes, C.J. (2010) Equivalent modal damping ratios for concrete/steel mixed structures. Computers and Structures 88(19-20), 1124-1136.
• Paultre, P. (2011) Dynamics of structures. Croydon: John Wiley & Sons.
• PEER Pacific Earthquake Engineering Research Center Ground Motion Database. [online]. Available from: https://ngawest2.berkeley.edu/ [Accessed July 27, 2022].
  » https://ngawest2.berkeley.edu/
• Petyt, M. (2010) Introduction to finite element vibration analysis, second edition 2nd ed. Cambridge: Cambridge University Press.
• Pimenta, P.M., Campello, E.M.B., Wriggers, P. (2008) An exact conserving algorithm for nonlinear dynamics with rotational DOFs and general hyperelasticity. Part 1: Rods. . 42(5), 715-732.
• Rong, B., Rui, X., Tao, L., Wang, G. (2019) Theoretical modeling and numerical solution methods for flexible multibody system dynamics. Nonlinear Dynamics 98(2), 1519-1553.
• Saaed, T.E., Nikolakopoulos, G., Jonasson, J.-E., Hedlund, H. (2015) A state-of-the-art review of structural control systems. JVC/Journal of Vibration and Control 21(5), 919-937.
• Simo, J.C., Tarnow, N., Wong, K.K. (1992) Exact energy-momentum conserving algorithms and symplectic schemes for nonlinear dynamics. Computer Methods in Applied Mechanics and Engineering 100(1), 63-116.
• Simo, J.C., Wong, K.K. (1991) Unconditionally stable algorithms for rigid body dynamics that exactly preserve energy and momentum. International Journal for Numerical Methods in Engineering 31(1), 19-52.
• Siqueira, T.M., Coda, H.B. (2016) Development of sliding connections for structural analysis by a total lagrangian FEM formulation. Latin American Journal of Solids and Structures 13(11).
• Siqueira, T.M., Coda, H.B. (2019) Flexible actuator finite element applied to spatial mechanisms by a finite deformation dynamic formulation. Computational Mechanics 64(6).
• Siqueira, T.M., Coda, H.B. (2017) Total Lagrangian FEM formulation for nonlinear dynamics of sliding connections in viscoelastic plane structures and mechanisms. Finite Elements in Analysis and Design 129, 63-77.
• Vakis, A.I., Yastrebov, V.A., Scheibert, J., Nicola, L., Dini, D., Minfray, C., Almqvist, A., Paggi, M., Lee, S., Limbert, G., Müser, M.H., Ciavarella, M. (2018) Modeling and simulation in tribology across scales: An overview. Tribology International 125, 169-199.
• Verlinden, O., Huynh, H.N., Kouroussis, G., Rivière-Lorphèvre, E. (2018) Modelling of flexible bodies with minimal coordinates by means of the corotational formulation. Multibody System Dynamics 42(4), 495-514.
• Yang, Y.-B., Lin, B.-H. (1995) Vehicle-Bridge Interaction Analysis by Dynamic Condensation Method. Journal of Structural Engineering 121(11), 1636-1643.
• Yin, X., Fang, Z., Cai, C.S., Deng, L. (2010) Non-stationary random vibration of bridges under vehicles with variable speed. Engineering Structures 32(8), 2166-2174.
• 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.
• Zienkiewicz, O., Taylor, R., Zhu, J.Z. (2013) The Finite Element Method: its Basis and Fundamentals: Seventh Edition.

Edited by

Editor: Marco L. Bittencourt

Publication Dates

• Publication in this collection
  28 Oct 2022
• Date of issue
  2022

History

• Received
  12 Sept 2022
• Reviewed
  07 Oct 2022
• Accepted
  10 Oct 2022