A New Method to Find the Forced Response of Nonlinear Systems With Dry Friction

Altamirano, Gregory; Tien, Meng-Hsuan; D'Souza, Kiran

doi:10.1115/1.4050686

Abstract

Coulomb friction has an influence on the behavior of numerous mechanical systems. Coulomb friction systems or dry friction systems are nonlinear in nature. This nonlinear behavior requires complex and time-demanding analysis tools to capture the dynamics of these systems. Recently, efforts have been made to develop efficient analysis tools able to approximate the forced response of systems with dry friction. The objective of this paper is to introduce a methodology that assists in these efforts. In this method, the piecewise linear nonlinear response is separated into individual linear responses that are coupled together through compatibility equations. The new method is demonstrated on a number of systems of varying complexity. The results obtained by the new method are validated through the comparison with results obtained by time integration. The computational savings of the new method are also discussed.

1 Introduction

Nonlinear dynamical systems are relevant to practical engineering systems. This paper focuses on a subset of these nonlinear dynamic systems involving dry friction or Coulomb friction. Dry friction has an influence on a variety of civil, mechanical, and aerospace applications [1–7]. For example, underplatform friction dampers have been intentionally included in the design of bladed disks in turbomachinery to mitigate unwanted vibration that occurs during operation under a number of conditions [8–11]. Also, there are efforts to mitigate earthquake hazards by developing energy dissipation devices for buildings and infrastructures to enhance their safety and reliability [12,13]. These mechanisms employ friction damping such as frictional sliding bearings in their designs to reduce the level of vibration experienced in the structures. In these applications and many others, nonlinearity is often localized in a particular area of the system; however, the models used to capture the dynamics of these systems are often very computationally expensive. As a result, there is considerable research focused on methods and techniques for estimating and understanding the impact of dry friction in the response of complex dynamic systems in an efficient manner.

In order to properly predict the dynamics of these Coulomb friction systems, several techniques have been developed to handle their nonlinear features that cannot be captured using traditional linear techniques. For example, numerical integration (NI) methods such as Runge–Kutta methods [14] have been developed for the analysis of general nonlinear systems. Although NI provides a good approximation for nonlinear responses of most nonlinear dynamical systems, its computational cost is often prohibitive when the system is high dimensional due to the need to keep the time-step size small to obtain accurate results. Another set of approaches based on the harmonic balance method have been implemented to analyze the periodic responses of complex Coulomb friction systems. These hybrid frequency/time (HFT) domain methods [15–19] can be combined with reduced-order modeling techniques to speed up the computation. Although HFT domain techniques provide a relatively efficient analysis tool for Coulomb friction systems, the computational cost of these methods can still be considerable if numerous harmonics are required to compute the system response.

Recently, a few hybrid analytical-numerical methods named the hybrid symbolic-numeric computational method [20,21] and the bilinear amplitude approximation (BAA) method [22,23] have been developed to better understand the dynamics of a subset of systems experiencing piecewise linear nonlinearities that result from intermittent contact. These methods are more efficient than traditional NI and HFT methods since they use efficient linear techniques in their response approximation process. The BAA method, in particular, has been developed to approximate the periodic response of piecewise linear nonlinear systems and has been implemented to study the vibrational response of complex structural systems such as bladed disks with intermittent contact [24,25]. Moreover, the efficiency of the BAA method has enabled the design of new energy harvesters that operate effectively over large frequency ranges [26]. This work builds upon a method to analyze the response of systems with dry friction called the friction response approximation method (FRAMe) [27], whose key idea stems from the BAA method, and was used to analyze one and three degree-of-freedom systems. In FRAMe, the overall nonlinear steady-state response is deconstructed into time intervals where the system behaves linearly. A nonlinear optimization solver is then employed to determine the time intervals and unknown coefficients in the symbolic expressions for the responses in each time interval. A nonlinear vibrational cycle is then approximated by stitching together these linear responses.

The remainder of the paper is organized as follows: The methodology section discusses how FRAMe can be applied to a generalized N-degree-of-freedom (DOF) system. The result section discusses the performance of FRAMe as compared to NI for a single DOF system and several higher dimensional systems. Along with these results, an analysis was conducted on a five DOF system with base excitation to determine the ideal location for a friction damper. Lastly, the conclusion section emphasizes the main points made throughout this work.

2 Methodology

In this section, FRAMe is introduced and presented as a technique for effectively determining the forced response of nonlinear systems subject to dry friction. The idea behind FRAMe is that the overall system response can be computed efficiently by computing the response of the system in its individual linear states and then enforcing the appropriate boundary conditions between these states. Figure 1 shows an N-DOF mass-spring-damper system with a single friction damper. The system is harmonically forced and experiences two types of damping: viscous and friction damping. The mass, stiffness, and damping coefficients are $m_{i}, k_{i}$ ⁠, and $c_{i}$ ⁠, respectively, and $f_{d}$ refers to the friction damper. The static frictional force is given as $f_{s} = μ_{s} P$ and the kinetic frictional force is force is given as $f_{k} = μ_{k} P$ ⁠. Note that $μ_{s}$ and $μ_{k}$ are the coefficients of static and kinetic friction and P is a loading force value for the damper. Figure 2 shows the friction characteristics for the systems analyzed in this paper. The direction of the frictional force is dependent on the relative velocity of the masses subjected to the friction damper (i.e., the friction force opposes the relative velocity of the masses). The system's overall motion can be defined by three states that depend on the relative velocity of the components that connect to the friction damper: (1) positive relative velocity, (2) negative relative velocity, and (3) zero velocity. The system is sliding while in states (1) and (2) and sticking while in state (3). The number of DOF reduces by one while in state (3) and the forces due to the static displacement of the stuck mass must be properly accounted for. Note that the friction damper can be placed between any of the masses or between the first mass and the ground.

The equations describing the system dynamics for the N-DOF system can be separated into its three states denoted by subscripts s1, s2, and re, where re is the reduced DOF system

\begin{matrix} M {\ddot{u}}_{s 1} (t) + C {\dot{u}}_{s 1} (t) + K u_{s 1} (t) = f (t) - f_{r} \\ M {\ddot{u}}_{s 2} (t) + C {\dot{u}}_{s 2} (t) + K u_{s 2} (t) = f (t) + f_{r} \\ M_{re} {\ddot{u}}_{s 3} (t) + C_{re} {\dot{u}}_{s 3} (t) + K_{re} u_{s 3} (t) = f_{re} (t) + f_{c} \end{matrix}

(1)

$[M, C, K]$ are the mass, damping, and stiffness matrices of the N-DOF system; $[M_{re}, C_{re}, K_{re}]$ are the mass, damping, and stiffness matrices of the reduced DOF system; $f (t)$ is the periodic forcing vector for the N-DOF system; $f_{re} (t)$ is the periodic external forcing vector for this system when it is in the stuck state (where the system is reduced by one DOF); $f_{r}$ is a forcing vector that contains the frictional force; and $f_{c}$ is a the forcing vector that results for state (3) when one of the masses is stuck and includes any of the static loads induced on the masses by a spring that is not in its unstretched position. In this work, the viscous damping matrix used was proportional to the stiffness matrix.

The coordinate transformation shown below shifts the coordinate systems to their static equilibrium position. For states (1) and (2), the static displacement (i.e.,

δ_{s 1}

and

δ_{s 2}

⁠) is due to the friction force and for state (3), the static displacement (i.e.,

δ_{s 3}

⁠) is due to the displacement of the stuck mass. The new coordinates and their static displacements are given as

\begin{matrix} u_{s 1} (t) = {\tilde{u}}_{s 1} + δ_{s 1}, δ_{s 1} = - K^{- 1} f_{r} \\ u_{s 2} (t) = {\tilde{u}}_{s 2} + δ_{s 2}, δ_{s 2} = K^{- 1} f_{r} \\ u_{s 3} (t) = {\tilde{u}}_{s 3} + δ_{s 3}, δ_{s 3} = K_{re}^{- 1} f_{c} \end{matrix}

(2)

The updated equations describing the system dynamics in the shifted coordinate systems for the three states are given by

\begin{matrix} M {\ddot{\tilde{u}}}_{s 1} (t) + C {\dot{\tilde{u}}}_{s 1} (t) + K {\tilde{u}}_{s 1} (t) = f (t) \\ M {\ddot{\tilde{u}}}_{s 2} (t) + C {\dot{\tilde{u}}}_{s 2} (t) + K {\tilde{u}}_{s 2} (t) = f (t) \\ M_{re} {\ddot{\tilde{u}}}_{s 3} (t) + C_{re} {\dot{\tilde{u}}}_{s 3} (t) + K_{re} {\tilde{u}}_{s 3} (t) = f_{re} (t) \end{matrix}

(3)

Note that the first two expressions in Eq. have the same coordinates except for a translation shift due to the friction force changing signs. Next, applying the modal transformations

{\tilde{u}}_{s 1} (t) = Φ {\tilde{q}}_{s 1}, {\tilde{u}}_{s 2} (t) = Φ {\tilde{q}}_{s 2}

⁠, and

{\tilde{u}}_{s 3} (t) = Φ_{re} {\tilde{q}}_{s 3}

⁠, where

{\tilde{q}}_{s 1}, {\tilde{q}}_{s 2}

⁠, and

{\tilde{q}}_{s 3}

are modal coordinates and

Φ

and

Φ_{re}

are the mode shape matrices for the N-DOF system and the reduced DOF system, respectively. Equation can be projected along

Φ

and

Φ_{re}

to obtain

\begin{matrix} Φ^{T} M Φ {\ddot{\tilde{q}}}_{s 1} (t) + Φ^{T} C Φ {\dot{\tilde{q}}}_{s 1} (t) + Φ^{T} K Φ {\tilde{q}}_{s 1} (t) = Φ^{T} f (t) \\ Φ^{T} M Φ {\ddot{\tilde{q}}}_{s 2} (t) + Φ^{T} C Φ {\dot{\tilde{q}}}_{s 2} (t) + Φ^{T} K Φ {\tilde{q}}_{s 2} (t) = Φ^{T} f (t) \\ Φ_{re}^{T} M_{re} Φ_{re} {\ddot{\tilde{q}}}_{s 3} (t) + Φ_{re}^{T} C_{re} Φ_{re} {\dot{\tilde{q}}}_{s 3} (t) + Φ_{re}^{T} K_{re} Φ_{re} {\tilde{q}}_{s 3} (t) = Φ_{re}^{T} f_{re} (t) \end{matrix}

(4)

Since the damping matrix is proportional to the stiffness matrix, the mass, damping, and stiffness matrices are all diagonalized by the modal projection yielding

\begin{matrix} {\ddot{\tilde{q}}}_{s 1, c} (t) + 2 ζ_{c} ω_{n, c} {\dot{\tilde{q}}}_{s 1, c} (t) + ω_{n, c}^{2} {\tilde{q}}_{s 1, c} (t) = f_{s 1, c}, c = 1, \dots, n_{s 1} \\ {\ddot{\tilde{q}}}_{s 2, j} (t) + 2 ζ_{j} ω_{n, j} {\dot{\tilde{q}}}_{s 2, j} (t) + ω_{n, j}^{2} {\tilde{q}}_{s 2, j} (t) = f_{s 2, j}, j = 1, \dots, n_{s 2} \\ {\ddot{\tilde{q}}}_{s 3, k} (t) + 2 ζ_{re, k} ω_{n_{re}, k} {\dot{\tilde{q}}}_{s 3, k} (t) + ω_{n_{re}, k}^{2} {\tilde{q}}_{s 3, k} (t) = f_{s 3, k}, k = 1, \dots, n_{re} \end{matrix}

(5)

where

[ω_{n}, ω_{n_{re}}]

are undamped frequencies of the system;

[ζ, ζ_{re}]

are viscous damping ratios; f_si are the forcing in modal coordinates; and

n_{s 1}, n_{s 2}

⁠, and

n_{re}

correspond to the number of retained modes in the reduced system. The solution to the differential equations in Eq. are composed of both the homogenous and particular solutions for states (1), (2), and (3)

\begin{matrix} {\tilde{q}}_{s 1, c} (t) = G_{c} sin (ω t - ϕ_{c} + α) + e^{- ζ_{c} ω_{n, c} t} [C_{1, c} cos (ω_{d, c} t) + C_{2, c} cos (ω_{d, c} t)] \\ {\tilde{q}}_{s 2, j} (t) = G_{j} sin (ω t - ϕ_{j} + α) + e^{- ζ_{j} ω_{n, j} t} [C_{3, j} cos (ω_{d, j} t) + C_{4, j} cos (ω_{d, j} t)] \\ {\tilde{q}}_{s 3, k} (t) = G_{re, k} sin (ω t - ϕ_{re, k} + α) + e^{- ζ_{re, k} ω_{n_{re}, k} t} [C_{5, k} cos (ω_{d_{re}, k} t) + C_{6, k} cos (ω_{d_{re}, k} t)] \end{matrix}

(6)

$[G, G_{re}]$ are coefficients for the particular solution, $[ϕ, ϕ_{re}]$ are the phase shifts for the particular solution, α is an additional phase shift caused by the nonlinearity, $[C_{1}, C_{2}, C_{3}, C_{4}, C_{5}, C_{6}]$ are the amplitudes of the transient response, and $[ω_{d}, ω_{d_{re}}]$ are the damped frequencies. Note that the unknowns in Eq. (6) are $[C_{1}, C_{2}, C_{3}, C_{4}, C_{5}, C_{6}]$ and α, the rest of the variables can be computed from linear vibration theory. It is important to note that $[C_{1}, C_{2}, C_{3}, C_{4}, C_{5}, C_{6}]$ are for transient response amplitudes after switching and are nonzero due to the switching occurring repeatedly.

The relative displacement, $u_{r} = u^{n l_{2}} - u^{n l_{1}}$ ⁠, generally experiences two different vibrational cycles and they are referred to as cases (1) and (2). For case (1), a sticking state does not occur and the nonlinear vibrational response is comprised of linear states (1) and (2) as shown in Fig. 3(a). For case (2), two sticking states occur during the response, which is composed of states (1), (2), and (3) as is shown in Fig. 3(b). Note that there are two additional cases where sticking only occurs at one location during the vibration cycle. These cases are discussed in the Appendix for completeness; however, they did not occur in the results presented in this work. Moreover, all of these cases may also have period-2 or higher responses as well. Although not studied in this work, these responses can be readily computed using FRAMe by increasing the number of transitions between states accordingly. One can also readily observe that the solution can capture harmonic distortion effects due to the nonlinearity; although these effects may be weak for case (1), they can become much stronger for case (2) as is evident from Fig. 3.

The amplitudes of the transient response (⁠ $[C_{1}, C_{2}, C_{3}, C_{4}, C_{5}, C_{6}]$ ⁠), the phase shift α, and the time periods spent in states (1), (2), and (3) are unknown quantities. To solve for these quantities, displacement and velocity continuity at the transition between states as well as matching forces when the system enters state (3) is enforced. The overall steady-state system response is then obtained by stitching together the individual responses for the system in each state. This method can be used to conduct a frequency sweep to determine the system response over a wide frequency spectrum to capture multiple resonant frequencies and the corresponding amplitudes of motion for a given excitation level.

Each case has a set of particular compatibility equations that were derived by understanding the piecewise linear nonlinear nature of the response. The overall system response is generated for each case by connecting the appropriate linear responses to each other. There are different expressions for states (1), (2), and (3) as shown in Eq. (6) and at a particular time, they can be equal to each other. For example, with regard to case (1), the displacement at the transition from state (1) to (2) is the same value at that instant for both expressions (i.e., T₁ in Fig. 3(a)). For both cases, the displacement at the beginning of the vibrational cycle is equal to the displacement at the end of the cycle. The displacement continuity equations are used to enforce this behavior. The velocity continuity equations are based on the fact that the velocity is equal to zero at the maximum and minimum displacements. Note that for case (1), there are only two time periods that need to be determined. For case (2), there are additional time periods that are required due to the inclusion of state (3) in the response. Therefore, the number of unknowns increases and additional equations to constrain the system response are needed and were chosen based on force balancing. The force balance equations are based on the criteria for motion to begin at the friction damper by overcoming the static frictional force. At the instant motion starts to occur, the forces acting on the system are equal to the maximum frictional force. This occurs at the beginning of the vibrational cycle and at the transition from state (3) to state (2). The equations to ensure compatibility between states for case (1) are based only on the displacement and velocity continuity. The equations to ensure compatibility between states for case (2) are based on continuity of displacement and velocity as well as balancing the force when the system exits state (3). The compatibility equations for cases (1) and (2) are shown below (note that the superscript ξ refers to all of the DOFs except the ith DOF):

Case (1)

Velocity continuity

\begin{matrix} {\dot{u}}_{r, s 1} (0) = 0 \\ {\dot{u}}_{r, s 1} (T_{1}) = 0 \\ {\dot{u}}_{r, s 2} (T_{1}) = 0 \\ {\dot{u}}_{r, s 2} (T) = 0 \\ {\dot{u}}_{s 1}^{ξ} (0) = {\dot{u}}_{s 2}^{ξ} (T) \\ {\dot{u}}_{s 1}^{ξ} (T_{1}) = {\dot{u}}_{s 2}^{ξ} (T_{1}) \end{matrix}

(7)

Displacement continuity

\begin{matrix} u_{s 1} (0) = u_{s 2} (T) \\ u_{s 1} (T_{1}) = u_{s 2} (T_{1}) \end{matrix}

Therefore case 1 has

4 N + 2

compatibility equations as shown in Eq. .

Case (2)

Velocity continuity

\begin{matrix} {\dot{u}}_{r, s 1} (0) = 0 \\ {\dot{u}}_{r, s 1} (T_{1}) = 0 \\ {\dot{u}}_{r, s 2} (T_{1} + T_{2}) = 0 \\ {\dot{u}}_{r, s 2} (T_{1} + T_{2} + T_{3}) = 0 \\ {\dot{u}}_{s 1}^{ξ} (0) = {\dot{u}}_{s 3} (T) \\ {\dot{u}}_{s 1}^{ξ} (T_{1}) = {\dot{u}}_{s 3} (T_{1}) \\ {\dot{u}}_{s 3} (T_{1} + T_{2}) = {\dot{u}}_{s 2}^{ξ} (T_{1} + T_{2}) \\ {\dot{u}}_{s 2}^{ξ} (T_{1} + T_{2} + T_{3}) = {\dot{u}}_{s 3} (T_{1} + T_{2} + T_{3}) \end{matrix}

(8)

Displacement continuity

\begin{matrix} u_{r, s 1} (0) = u_{r, s 2} (T_{1} + T_{2} + T_{3}) \\ u_{r, s 1} (T_{1}) = u_{r, s 2} (T_{1} + T_{2}) \\ u_{s 1}^{ξ} (0) = u_{s 3} (T) \\ u_{s 1}^{ξ} (T_{1}) = u_{s 3} (T_{1}) \\ u_{s 3} (T_{1} + T_{2}) = u_{s 2}^{ξ} (T_{1} + T_{2}) \\ u_{s 2}^{ξ} (T_{1} + T_{2} + T_{3}) = u_{s 3} (T_{1} + T_{2} + T_{3}) \end{matrix}

Force balance

\begin{matrix} f_{i} (0) - k_{i} u_{r, s 1} (0) - k_{i + 1} (u_{i, s 1}^{n l_{1}} (0) - u_{i + 1, s 1} (0)) + c ({\dot{u}}_{i, s 1}^{n l_{1}} (0) - {\dot{u}}_{i + 1, s 1} (0)) = f_{s} \\ f_{i} (T_{1} + T_{2}) - k_{i} u_{r, s 2} (T_{1} + T_{2}) - k_{i + 1} (u_{i, s 2}^{n l_{1}} (T_{1} + T_{2}) - u_{i + 1, s 2} (T_{1} + T_{2})) + c ({\dot{u}}_{i, s 2}^{n l_{1}} (T_{1} + T_{2}) - {\dot{u}}_{i + 1, s 2} (T_{1} + T_{2})) = - f_{s} \end{matrix}

Therefore, case 2 has 8 N compatibility equations as shown in Eq. (8).

The optimization problem is formulated by moving all terms in Eqs. (7) and (8) to the left side of the equations to yield a set of objective functions to be minimized. Note that the number of unknown variables for Eq. (7) is $2 n_{s 1} + 2 n_{s 2} + 2$ and the number of unknown variable for Eq. (8) is $2 n_{s 1} + 2 n_{s 2} + 4 n_{re} + 4$ ⁠. The system will have more compatibility equations than unknown variables if the number of modes retained is less than the number of DOF; otherwise, the number of compatibility equations and unknown variables will be equal. Note that the method is well-suited to solve the steady-state response for any periodic excitation and not just a single harmonic forcing term. The only part of the formulation that changes is the form of the particular solution in Eq. (6), which is considered known from linear vibration theory and not solved for using the compatibility equations.

3 Results

In this section, results for several systems of increasing complexity are presented and discussed. FRAMe and time integration were both used to generate forced response plots for comparison purposes. Recall that the amplitudes of the transient response (⁠ $[C_{1}, C_{2}, C_{3}, C_{4}, C_{5}, C_{6}]$ ⁠), the phase shift α, and the time periods spent in states (1), (2), and (3) are unknown parameters. To solve for these parameters, the matlab function “lsqnonlin,” a nonlinear optimization solver, was utilized with the appropriate compatibility conditions (e.g., Eqs. (7) and (8)). This matlab function is a Jacobian-based solver and requires a Jacobian matrix to complete its computational process. For the results in this work, a Jacobian matrix is approximated using the finite difference method, which has been automatically implemented in the solver. A random number generator was used to generate the initial guesses for the unknown parameters in Eq. (6) for the first frequency point in the analysis. The first frequency point results in a steady-state response experiencing case (1). Note that high magnitude excitation was specifically chosen to induce case (1) initially to demonstrate how FRAMe can transition from case (1) to case (2) seamlessly. For case (1), T₁ is the only time period that needs to be solved for. The initial guess for T₁ was chosen to be half of the time period of the vibration cycle T. Note that $T = 2 π / ω$ ⁠, where ω is the excitation frequency. The next frequency point uses the results obtained from the first frequency as the new initial guess for all of the unknown parameters. The process of using the result from the previous frequency point as the new initial guess for the next frequency point is continued for the remainder of the frequency sweep. The CPU time for the first frequency point is used when comparing the computational time of FRAMe to time integration to ensure a fair comparison between the methods. Note that for time integration, the initial conditions used were zero and the system was time integrated until it reached a steady-state value (i.e., the difference in the amplitude between consecutive peaks was below a threshold).

3.1 One Degree of Freedom System.

Consider the system shown in Fig. 1 with N = 1 subjected to the forcing, $f (t) = f_{o} sin (ω t)$ that experiences both viscous and friction damping. A frequency sweep analysis was conducted on this system with both FRAMe and numerical integration. The system and forcing parameters used were m = 1 kg, c = 1 N·s/m, P = 9.5 N, $μ_{k} = 0.2, μ_{s} = 0.4$ ⁠, k = 50 N/m, and f_o = 4 N. A frequency versus amplitude plot is shown in Fig. 4(a) and contains a frequency range where the system experiences both cases (1) and (2). The frequency sweep begins near the peak amplitude where the system is experiencing case (1) and then sweeps to the left and right from the peak. When the system experiences case (1), Eq. (7) simplifies to:

Velocity continuity

\begin{matrix} {\dot{u}}_{s 1} (0) = 0 \\ {\dot{u}}_{s 1} (T_{1}) = 0 \\ {\dot{u}}_{s 2} (T_{1}) = 0 \\ {\dot{u}}_{s 2} (T) = 0 \end{matrix}

(9)

Displacement continuity

\begin{matrix} u_{s 1} (T_{1}) = u_{s 2} (T_{1}) \\ u_{s 1} (0) = u_{s 2} (T) \end{matrix}

Criteria were developed to determine at which frequency point the transition from case (1) to (2) occurred. Each frequency point in the case (1) frequency range was checked to see if the following criteria were satisfied. The criteria being that $| f (T_{1}) - k u (T_{1}) | > f_{s}$ and $| f (T_{1} + T_{2}) - k u (T_{1} + T_{2}) | > f_{s}$ ⁠. Once these criteria were no longer satisfied, the resultant forces acting on the system, at the maximum or minimum displacement, were no longer greater than the maximum frictional force and state (3) would occur. At this frequency point, the set of compatibility equations for case (2) were used. The compatibility equations for case (2) (i.e., Eq. (8)) are shown below for this single DOF system:

Velocity continuity

\begin{matrix} {\dot{u}}_{s 1} (0) = 0 \\ {\dot{u}}_{s 1} (T_{1}) = 0 \\ {\dot{u}}_{s 2} (T_{1} + T_{2}) = 0 \\ {\dot{u}}_{s 2} (T_{1} + T_{2} + T_{3}) = 0 \end{matrix}

Displacement continuity

\begin{matrix} u_{s 1} (T_{1}) = u_{s 2} (T_{1} + T_{2}) \\ u_{s 1} (0) = u_{s 2} (T) \end{matrix}

(10)

Force balance

\begin{matrix} f (0) - k u_{s 1} (0) = f_{s} \\ f (T_{1} + T_{2}) - k u_{s 2} (T_{1} + T_{2}) = - f_{s} \end{matrix}

Figure 4(b) shows the relative amplitude percent difference of FRAMe with respect to time integration over the entire frequency range. Note that FRAMe has an extremely low percent difference in the peak response region and the relative percent difference increases in the transition from sliding to sticking; however, the relative percent difference is less than 1% even in the small amplitude regime where there is sticking. For this very simple system, the computational time for numerical integration was 1.05 s/frequency point and 0.178 s/frequency point for FRAMe, which is almost an order of magnitude faster.

3.2 Mutlidegree-of-Freedom Systems

3.2.1 Computational Savings.

Consider an N-DOF system, similar to the system shown in Fig. 1, with the Nth mass excited by a harmonic forcing and the first mass connected to the wall by a friction damper. Note that four different N values are presented in this work (i.e., $N ϵ 50, 100, 300, 500$ ⁠). The system parameters are $m_{i} = 0.1$ kg, k_i = 500 N/m, (⁠ $i = 1, \dots, N$ ⁠) $β = 0.01$ ⁠, f_o = 1 N, P = 0.1 N, $μ_{k} = μ_{s} = μ = 1$ ⁠. Figure 5 shows amplitude versus frequency plots for the Nth mass for a number of N values. Time integration was used to determine the steady-state amplitudes for five different frequency points for each N value. Figure 5 shows the excellent agreement between numerical integration and FRAMe.

Figure 6 presents the relative percent difference at the peak amplitude and the computational time required to estimate the response using a number of different modes (i.e., 1, 10, $N / 2$ ⁠, and N) in the calculation. In particular, Fig. 6(a) shows a very low difference at the peak (less than 0.025%) when using as few as a single mode for systems up to 500 DOFs. There is no clear pattern in how the difference is related to the number of modes used in the calculation. In the cases studied, using 10 modes seemed to give the best results; however, the estimate is very accurate irrespective of the number of modes used. The accuracy of FRAMe depends on a couple of factors including the tolerance setting of the optimization solver (similar issue when performing numerical integration) and how accurate the Jacobian matrix in the optimization solver is computed. In Fig. 5(a), FRAMe shows a bit higher error when the system is high dimensional (N = 500). This is due to the fact that the Jacobian matrix was approximated using the finite difference method in this work. It is generally harder for the optimization solver to find a “very accurate” solution when analyzing a complex system if an analytical Jacobian matrix is not provided.

The choice of number of modes used in the estimate does impact the computational time for the algorithm as is shown in Fig. 6(b). Irrespective of the number of modes used in FRAMe, the computational cost increases as the number of DOFs increase. Figure 6(b) shows that FRAMe is generally at least one to two orders of magnitude faster than time integration. When comparing FRAMe using ten modes or less to time integration one can see the real computational savings as N increases from 50 to 500. The savings move from less than two orders of magnitude to over four orders of magnitude. Moreover, there is still very good agreement between FRAMe and time integration with less than 0.025% relative difference at the peak.

3.2.2 Parametric Study.

In this section, all of the results were obtained using FRAMe to study the five DOF system shown in Fig. 7. Note that the base displacement is given by the harmonic function, $u_{o} = U sin (ω t)$ ⁠. The system parameters were m = 5000 kg, k = 3000 N/m, $β = 0.01$ ⁠, U = . 005 m, and $μ_{k} = μ_{s} = μ = 1$ ⁠. In this work, different levels of the loading force value (P) were considered to measure its impact on the displacement of the fifth mass. The tested values for P were 1 N, 2 N,…10 N. Also, different location options for the friction damper were considered and those five options are shown in Fig. 7. The excitation frequency range considered excited the first and second modes of the system.

Figure 8 shows the impact of increasing the loading force value and changing the location of the friction damper on the amplitude of the fifth mass when exciting the first and second modes. In Fig. 8(a), the amplitude was at its lowest value when the damper was placed at location 1 and the loading force value was 10 N. The amplitude was at its highest value when the damper was placed at location 5 and the loading force value was 1 N. Figure 8(a) shows a clear trend in increasing the loading on the friction damper and placing it closer to the excitation source to minimize the amplitude of the fifth mass while the first mode is excited. The trend changes when the second mode of the system is excited. In Fig. 8(b), the amplitude was at its lowest value when the damper was placed at location 4 and the loading force value was 10 N.

4 Conclusions and Discussion

Understanding the effects of dry friction on the dynamics of complex dynamic systems is important in many fields. This work introduces a method to quickly estimate the response amplitude and resonant frequencies of these nonlinear systems. This new technique is called the FRAMe. FRAMe uses the fact that often the overall dynamics of these systems can be reconstructed by connecting the responses of the system in each of its linear states. This idea stems from the bilinear amplitude approximation method. There are a set of unknown parameters for the linear systems that are computed by enforcing the appropriate compatibility conditions for the transitions between the linear states. The key benefit of FRAMe is that it estimates the nonlinear system response by connecting a set of linear responses, which are very efficiently computed. The results presented in this work use a built-in matlab optimization function, which is a traditional local solver; however, the method is well-suited to be paired with a variety of global optimization methods that do not depend on initial values. Although the computational savings of using FRAMe over numerical integration for small systems is minor, a significant computational saving is achieved by applying FRAMe to high dimensional systems.

The piecewise linear nonlinear response of Coulomb friction systems is efficiently and accurately predicted by using FRAMe. The new method was demonstrated on a variety of systems of increasing complexity. The results for the one degree-of-freedom system and several higher dimensional systems were also calculated by using time integration for validation purposes. The computational cost for numerical integration increased directly with the dimension of the system; however, the computational cost of FRAMe increased much less dramatically as the dimensionality of the system increased, while still maintaining accurate results. Lastly, a study was conducted on a five degree-of-freedom system with base excitation. Different locations for the friction damper and levels of loading force values were used to show the effectiveness of the method in efficiently carrying out parametric studies of systems with friction dampers.

Coulomb friction nonlinearities present themselves in a broad array of circumstances in nature and engineering applications, and can exhibit a variety of complex behaviors particularly when multiple friction nonlinearities exist. This work presents a new method that was created to initially handle a subset of these complex cases. In particular, it is focused on multidegree-of-freedom systems with a single friction nonlinearity. Work to extend this method to multiple friction nonlinearities is currently underway. A key challenge is the dramatically increasing number of states that become possible as the number of friction nonlinearities increase. Moreover, combining this method with existing methods, like the generalized bilinear amplitude approximation method [23–25], to analyze systems with piecewise linear stiffness and friction is another promising direction of future work. Another important point is that the current method is focused on the steady-state response of systems excited by any form of periodic excitation. This steady-state response can capture various stick and slip behaviors in combination as long as the response is periodic. A technique that can handle transient, chaotic, or other nonstationary dynamics requires a different solution process such as a modified version of the recently developed hybrid symbolic numeric computational method [20,21]. These methods can work in combination to rapidly solve for any system response of high dimensional nonlinear systems with a single friction nonlinearity.

Acknowledgment

This paper is based on work partially supported by the National Science Foundation under grant no.1902408, program manager Dr. Robert Landers. Any opinions, findings, and conclusions or recommendations expressed in this paper are those of the authors and do not necessarily reflect the views of the National Science Foundation. G. A. has been supported by the SROP fellowship and the College of Engineering fellowship and M.-H.T. has been supported by the Presidential Fellowship from The Ohio State University Graduate School while working on this project.

Funding Data

National Science Foundation (Grant No. 1902408; Funder ID: 10.13039/100000001).
SROP Fellowship.
College of Engineering Fellowship.
Presidential Fellowship from The Ohio State University Graduate School.

Appendix

For completeness, the compatibility equations for two additional cases are given in the Appendix. Although these cases did not occur in the results presented in this work, they can occur if the system is subject to an oscillating normal load in the friction damper.

Case 3: State (3) Occurs at Maximum Displacement Only.

Figure 9 shows the response of the relative displacement, u_r, when state (3) occurs at the maximum displacement of the relative coordinate. The case 3 compatibility equations are given by

Velocity continuity

\begin{matrix} {\dot{u}}_{r, s 1} (0) = 0 \\ {\dot{u}}_{r, s 1} (T_{1}) = 0 \\ {\dot{u}}_{r, s 2} (T_{1} + T_{2}) = 0 \\ {\dot{u}}_{r, s 2} (T) = 0 \\ {\dot{u}}_{s 1}^{ξ} (0) = {\dot{u}}_{s 2}^{ξ} (T) \\ {\dot{u}}_{s 1}^{ξ} (T_{1}) = {\dot{u}}_{s 3} (T_{1}) \\ {\dot{u}}_{s 3} (T_{1} + T_{2}) = {\dot{u}}_{s 2}^{ξ} (T_{1} + T_{2}) \end{matrix}

(11)

Displacement continuity

\begin{matrix} u_{r, s 1} (T_{1}) = u_{r, s 2} (T_{1} + T_{2}) \\ u_{s 1} (0) = u_{s 2} (T) \\ u_{s 1}^{ξ} (T_{1}) = u_{s 3} (T_{1})' \\ u_{s 3} (T_{1} + T_{2}) = u_{s 2}^{ξ} (T_{1} + T_{2}) \end{matrix}

Force balance

F_{i} (T_{1} + T_{2}) - k_{i} u_{r, s 2} (T_{1} + T_{2}) - k_{i + 1} (u_{i, s 2}^{n l_{1}} (T_{1} + T_{2}) - u_{i + 1, s 2} (T_{1} + T_{2})) + c ({\dot{u}}_{i, s 2}^{n l_{1}} (T_{1} + T_{2}) - {\dot{u}}_{i + 1, s 2} (T_{1} + T_{2})) = - f_{s}

Case 4: State (3) Occurs at Minimum Displacement Only.

Figure 10 shows the response of the relative displacement, u_r, when state (3) occurs at the minimum displacement of the relative coordinate. The case 4 compatibility equations are given by

Velocity continuity

\begin{matrix} {\dot{u}}_{r, s 1} (0) = 0 \\ {\dot{u}}_{r, s 1} (T_{1}) = 0 \\ {\dot{u}}_{r, s 2} (T_{1}) = 0 \\ {\dot{u}}_{r, s 2} (T_{1} + T_{2}) = 0 \\ {\dot{u}}_{s 1}^{ξ} (0) = {\dot{u}}_{s 3} (T) \\ {\dot{u}}_{s 1}^{ξ} (T_{1}) = {\dot{u}}_{s 2}^{ξ} (T_{1}) \\ {\dot{u}}_{s 2}^{ξ} (T_{1} + T_{2}) = {\dot{u}}_{s 3} (T_{1} + T_{2}) \end{matrix}

(12)

Displacement continuity

\begin{matrix} u_{r, s 1} (0) = u_{r, s 2} (T_{1} + T_{2}) \\ u_{s 1} (T_{1}) = u_{s 2} (T_{1}) \\ u_{s 1}^{ξ} (0) = u_{s 3} (T) \\ u_{s 2}^{ξ} (T_{1} + T_{2}) = u_{s 3} (T_{1} + T_{2}) \end{matrix}

Force balance

F_{i} (0) - k_{i} u_{r, s 1} (0) - k_{i + 1} (u_{i, s 1}^{n l_{1}} (0) - u_{i + 1, s 1} (0)) + c ({\dot{u}}_{i, s 1}^{n l_{1}} (0) - {\dot{u}}_{i + 1, s 1} (0)) = f_{s}

References

1.

Ibrahim

,

R. A.

,

1994

, “

Friction-Induced Vibration, Chatter, Squeal, and chaos-Part II: Dynamics and Modeling

,”

ASME Appl. Mech. Rev. Trans. ASME

,

47

(

7

), pp.

227

–

253

.10.1115/1.3111080

Google Scholar

Crossref

2.

Symans

,

M. D.

, and

Constantinou

,

M. C.

,

1999

, “

Semi-Active Control Systems for Seismic Protection of Structures: A State-of-the-Art Review

,”

Eng. Struct.

,

21

(

6

), pp.

469

–

487

.10.1016/S0141-0296(97)00225-3

Google Scholar

Crossref

3.

Yang

,

J. N.

, and

Agrawal

,

A. K.

,

2002

, “

Semi-Active Hybrid Control Systems for Nonlinear Buildings Against Near-Field Earthquakes

,”

Eng. Struct.

,

24

(

3

), pp.

271

–

280

.10.1016/S0141-0296(01)00094-3

Google Scholar

Crossref

4.

Awrejcewicz

,

J.

, and

Olejnik

,

P.

,

2005

, “

Analysis of Dynamic Systems With Various Friction Laws

,”

ASME Appl. Mech. Rev. Trans. ASME

,

58

(

6

), pp.

389

–

411

.10.1115/1.2048687

Google Scholar

Crossref

5.

Jung

,

C.

,

Epureanu

,

B. I.

,

Baik

,

S.

, and

Huffman

,

M. B.

,

2015

, “

Nonlinear Reduced Order Models for the Structural Dynamics of Combustor Systems With Prestress and Friction

,”

ASME J. Comput. Nonlinear Dyn.

,

10

(

1

), p.

011009

.10.1115/1.4028014

Google Scholar

Crossref

6.

Al-Zughaibi

,

A. I. H.

,

2018

, “

Experimental and Analytical Investigations of Friction at Lubricant Bearings in Passive Suspension Systems

,”

Nonlinear Dyn.

,

94

(

2

), pp.

1227

–

1242

.10.1007/s11071-018-4420-x

Google Scholar

Crossref

7.

Xing

,

M.

,

Zhou

,

L.

,

Zhang

,

C.

,

Xue

,

K.

, and

Zhang

,

Z.

,

2019

, “

Simulation Analysis of Nonlinear Friction of Rod String in Sucker Rod Pumping System

,”

ASME J. Comput. Nonlinear Dyn.

,

14

(

9

), p.

091008

.10.1115/1.4044027

Google Scholar

Crossref

8.

Zucca

,

S.

,

Firrone

,

C. M.

, and

Gola

,

M.

,

2013

, “

Modeling Underplatform Dampers for Turbine Blades: A Refined Approach in the Frequency Domain

,”

J. Vib. Control

,

19

(

7

), pp.

1087

–

1102

.10.1177/1077546312440809

Google Scholar

Crossref

9.

Gola

,

M. M.

, and

Liu

,

T.

,

2014

, “

A Direct Experimental-Numerical Method for Investigations of a Laboratory Underplatform Damper Behavior

,”

Int. J. Solids Struct.

,

51

(

25–26

), pp.

4245

–

4259

.10.1016/j.ijsolstr.2014.08.011

10.

Pesaresi

,

L.

,

Salles

,

L.

,

Jones

,

A.

,

Green

,

J.

, and

Schwingshackl

,

C.

,

2017

, “

Modelling the Nonlinear Behaviour of an Underplatform Damper Test Rig for Turbine Applications

,”

Mech. Syst. Signal Process.

,

85

, pp.

662

–

679

.10.1016/j.ymssp.2016.09.007

Google Scholar

Crossref

11.

Krack

,

M.

,

Salles

,

L.

, and

Thouverez

,

F.

,

2017

, “

Vibration Prediction of Bladed Disks Coupled by Friction Joints

,”

Arch. Comput. Methods Eng.

,

24

(

3

), pp.

589

–

636

.10.1007/s11831-016-9183-2

Google Scholar

Crossref

12.

Crewe

,

A.

,

1998

, “

Passive Energy Dissipation Systems in Structural Engineering

,”

Struct. Saf.

,

20

(

2

), pp.

197

–

198

.10.1016/S0167-4730(97)00034-9

Google Scholar

Crossref

13.

Miguel

,

L. F. F.

,

Miguel

,

L. F. F.

, and

Lopez

,

R. H.

,

2014

, “

Robust Design Optimization of Friction Dampers for Structural Response Control

,”

Struct. Control Health Monit.

,

21

(

9

), pp.

1240

–

1251

.10.1002/stc.1642

Google Scholar

Crossref

14.

Shampine

,

L. F.

, and

Reichelt

,

M. W.

,

1997

, “

The Matlab Ode Suite

,”

SIAM J. Sci. Comput.

,

18

(

1

), pp.

1

–

22

.10.1137/S1064827594276424

Google Scholar

Crossref

15.

Cameron

,

T. M.

, and

Griffin

,

J. H.

,

1989

, “

An Alternating Frequency/Time Domain Method for Calculating the Steady-State Response of Nonlinear Dynamic Systems

,”

ASME J. Appl. Mech.

,

56

(

1

), pp.

149

–

154

.10.1115/1.3176036

Google Scholar

Crossref

16.

Kim

,

T.

,

Rook

,

T.

, and

Singh

,

R.

,

2005

, “

Super- and Sub-Harmonic Response Calculations for a Torsional System With Clearance Nonlinearity Using the Harmonic Balance Method

,”

J. Sound Vib.

,

281

(

3–5

), pp.

965

–

993

.10.1016/j.jsv.2004.02.039

17.

Saito

,

A.

,

Castanier

,

M. P.

,

Pierre

,

C.

, and

Poudou

,

O.

,

2009

, “

Efficient Nonlinear Vibration Analysis of the Forced Response of Rotating Cracked Blades

,”

ASME J. Comput. Nonlinear Dyn. Trans.

,

4

(

1

), p.

011005

.10.1115/1.3007908

Google Scholar

Crossref

18.

Zucca

,

S.

, and

Epureanu

,

B. I.

,

2014

, “

Bi-Linear Reduced-Order Models of Structures With Friction Intermittent Contacts

,”

Nonlinear Dyn.

,

77

(

3

), pp.

1055

–

1067

.10.1007/s11071-014-1363-8

Google Scholar

Crossref

19.

Mashayekhi

,

F.

,

Nobari

,

A. S.

, and

Zucca

,

S.

,

2019

, “

Hybrid Reduction of Mistuned Bladed Disks for Nonlinear Forced Response Analysis With Dry Friction

,”

Int. J. Non-Linear Mech.

,

116

, pp.

73

–

84

.10.1016/j.ijnonlinmec.2019.06.001

Google Scholar

Crossref

20.

Tien

,

M.-H.

, and

D'Souza

,

K.

,

2019

, “

Analyzing Bilinear Systems Using a New Hybrid Symbolic-Numeric Computational Method

,”

J. Vib. Acoust.

,

141

(

3

), p. 031008.10.1115/1.4042520

21.

Tien

,

M.-H.

, and

D'Souza

,

K.

,

2019

, “

Transient Dynamics Analysis of Cracked Structures With Multiple Contact Pairs Using Generalized HSNC

,”

Nonlinear Dyn.

,

96

(

2

), pp.

1115

–

1131

.10.1007/s11071-019-04844-7

Google Scholar

Crossref

22.

Jung

,

C.

,

D'Souza

,

K.

, and

Epureanu

,

B. I.

,

2014

, “

Nonlinear Amplitude Approximation for Bilinear Systems

,”

J. Sound Vib.

,

333

(

13

), pp.

2909

–

2919

.10.1016/j.jsv.2014.01.029

Google Scholar

Crossref

23.

Tien

,

M.-H.

, and

D'Souza

,

K.

,

2017

, “

A Generalized Bilinear Amplitude and Frequency Approximation for Piecewise-Linear Nonlinear Systems With Gaps or Prestress

,”

Nonlinear Dyn.

,

88

(

4

), pp.

2403

–

2416

.10.1007/s11071-017-3385-5

Google Scholar

Crossref

24.

Tien

,

M.-H.

,

Hu

,

T.

, and

D'Souza

,

K.

,

2018

, “

Generalized Bilinear Amplitude Approximation and X-Xr for Modeling Cyclically Symmetric Structures With Cracks

,”

ASME J. Vib. Acoust.

,

140

(

4

), p.

041012

.10.1115/1.4039296

Google Scholar

Crossref

25.

Tien

,

M.-H.

,

Hu

,

T.

, and

D'Souza

,

K.

,

2019

, “

Statistical Analysis of the Nonlinear Response of Bladed Disks With Mistuning and Cracks

,”

AIAA J.

,

57

(

11

), pp.

4966

–

4977

.10.2514/1.J058190

Google Scholar

Crossref

26.

Tien

,

M.-H.

, and

D'Souza

,

K.

,

2020

, “

Method for Controlling Vibration by Exploiting Piecewise-Linear Nonlinearity in Energy Harvesters

,”

Proc. R. Soc. A

,

476

(

2233

), p.

20190491

.10.1098/rspa.2019.0491

Google Scholar

Crossref

27.

Altamirano

,

G. L.

,

Tien

,

M.-H.

, and

D'Souza

,

K.

,

2019

, “

Efficient Response Approximation for Systems With Coulomb Friction

,”

AIAA

Paper No. 2019-2258.10.2514/6.2019-2258

A New Method to Find the Forced Response of Nonlinear Systems With Dry Friction

Abstract

1 Introduction

2 Methodology

3 Results

3.1 One Degree of Freedom System.

3.2 Mutlidegree-of-Freedom Systems

3.2.1 Computational Savings.

3.2.2 Parametric Study.

4 Conclusions and Discussion

Acknowledgment

Funding Data

Appendix

Case 3: State (3) Occurs at Maximum Displacement Only.

Case 4: State (3) Occurs at Minimum Displacement Only.

References

Contents

Data & Figures

Supplements

References

Get Email Alerts

Cited By

ASME Journals

ASME Conference Proceedings

ASME eBooks

Resources

Opportunities

A New Method to Find the Forced Response of Nonlinear Systems With Dry Friction

Abstract

1 Introduction

2 Methodology

3 Results

3.1 One Degree of Freedom System.

3.2 Mutlidegree-of-Freedom Systems

3.2.1 Computational Savings.

3.2.2 Parametric Study.

4 Conclusions and Discussion

Acknowledgment

Funding Data

Appendix

Case 3: State (3) Occurs at Maximum Displacement Only.

Case 4: State (3) Occurs at Minimum Displacement Only.

References

Contents

Data & Figures

Supplements

References

Related

Get Email Alerts

Cited By

Related Articles

Related Proceedings Papers

Related Chapters

ASME Journals

ASME Conference Proceedings

ASME eBooks

Resources

Opportunities

This Feature Is Available To Subscribers Only