## Galerkin-Ivanov transformation for nonsmooth modeling of vibro-impacts in continuous structures

Surya Samukham , S. N. Khaderi , and C. P. Vyasarayani 1 1 ∗ 1

1 Department of Mechanical and Aerospace Engineering, Indian Institute of Technology Hyderabad, Kandi, Sangareddy, Telangana, India, 502285

## Abstract

This work deals with the modeling of nonsmooth vibro-impact motion of a continuous structure against a rigid distributed obstacle. Galerkin's approach is used to approximate the solutions of the governing partial differential equations of the structure, which results in a system of ordinary differential equations (ODEs). When these ODEs are subjected to unilateral constraints and velocity jump conditions, one must use an event detection algorithm to calculate the time of impact accurately. Event detection in the presence of multiple simultaneous impacts is a computationally demanding task. Ivanov (Ivanov, A., 1993. 'Analytical methods in the theory of vibro-impact systems', Journal of Applied Mathematics and Mechanics, 57(2), pp. 221-236.) proposed a nonsmooth transformation for a vibro-impacting multi-degree-of-freedom system subjected to a single unilateral constraint. This transformation eliminates the unilateral constraints from the problem and, therefore, no event detection is required during numerical integration. This nonsmooth transformation leads to sign function nonlinearities in the equations of motion. However, they can be easily accounted for during numerical integration. Ivanov used his transformation to make analytical calculations for the stability and bifurcations of vibro-impacting motions; however, he did not explore its application for simulating distributed collisions in spatially continuous structures. We adopt Ivanov's transformation to deal with multiple unilateral constraints in spatially continuous structures. Also, imposing the velocity jump conditions exactly in the modal coordinates is nontrivial and challenging. Therefore, in this work we use a modal-physical transformation to convert

∗ Address all correspondence to this author. E-mail address: vcprakash@iith.ac.in

the system from modal to physical coordinates on a spatially discretized grid. We then apply Ivanov's transformation on the physical system to simulate the vibro-impact motion of the structure. The developed method is demonstrated by modeling the distributed collision of a nonlinear string against a rigid distributed surface. For validation, we compare our results with the well-known penalty approach.

## 1 Introduction

In many engineering applications, structures are subjected to vibro-impacting motions [1-3]. Clearances in mechanical joints due to wear can lead to vibro-impacting motions in machine components. In some vibration-based energy harvesting applications, vibro-impacting motions are deliberately introduced to increase the bandwidth of the frequency response [4]. In addition, the study of vibro-impact motion of a structure has various applications in the synthesis of sound in stringed musical instruments of Indian origin, like the sitar and the tanpura [5-7]. Applications are also found in the elevator-rope collision problem [8], multi-body dynamics [9], and impact-damper systems [10,11]. However, in this work, we are interested in the problem of modeling the vibro-impact motion of a continuous structure constrained by a rigid distributed obstacle. When a continuous structure is spatially discretized using numerical approximation methods, there are three ways the impacting motion can be simulated. The first method approximates the rigid obstacle as a foundation of springs with high stiffness and is known as the penalty approach [5,12-16]. Alternatively, one can use Lagrange multipliers to impose contact constraints once the structure comes into contact with the obstacle. Complementarity conditions between Lagrange multipliers and gap functions are used to impose contact loss [17]. The Lagrange multiplier approach simulates the perfect sticking motion of vibro-impacting systems [17]. Impacts can also be simulated using a coefficient of restitution (CoR or R ) approach [12, 18-22]. Once contact is detected between the structure and the obstacle, the appropriate velocity jump conditions are imposed at the point of contact. In the Lagrange multiplier and CoR approaches, one must solve an event detection problem to isolate the time of impact accurately. Multiple simultaneous impacts are present when a continuous structure impacts a distributed obstacle, which makes the event detection problem computationally expensive. A few other recently developed numerical time-integration techniques to simulate the nonsmooth dynamical systems are given in [23], [24], and [25].

Within the framework of the CoR-based approach, Ivanov [26,27] proposed a nonsmooth spatial transformation that automatically satisfies both the unilateral constraints and the ve-

locity jump conditions at the point of contact. This approach eliminates the need for event detection. Ivanov proposed this method to study the vibro-impacting motion of a multi-degreeof-freedom (MDOF) system, where the displacement of a single mass is constrained. Ivanov's transformation was successfully used to study the nonlinear dynamics of a single-degree-offreedom (SDOF) vibro-impacting ship motion [28,29]. In this work, we adopt Ivanov's transformation to account for the multiple unilateral constraints in a MDOF system.

In this paper, we have implemented Ivanov's approach to simulate the vibro-impacting motion of a continuous structure against a rigid distributed obstacle. First, we discuss the implementation of Ivanov's transformation to a SDOF system, followed by its application to the MDOF system (continuous structure). To demonstrate the application of Ivanov's approach to a continuous structure subjected to distributed contact, we consider a nonlinear string and present the vibro-impacting motion of the string against the rigid flat and sinusoidal obstacles for different CoR. Galerkin approximation has been used to solve the governing partial differential equations (PDEs) of the string. The set of ordinary differential equations (ODEs) obtained from the Galerkin approach are in modal coordinates and the imposition of impact constraints upon the modal system is a challenging task [30-32]. To avoid such difficulties, we use a transformation to convert the modal system into its physical coordinates by discretizing the string in space [6,33]. We then apply Ivanov's transformation to the spatially discretized system to incorporate the impact constraints. After Ivanov's transformation, the resulting system of differential equations are in Ivanov's coordinates and non-stiff in nature. Therefore, an explicit integrator with a large time-step size can be used to solve the system of differential equations, which is not the case with the penalty method because of the presence of stiff differential equations. This advantage of implementing an explicit integrator with a larger time-step size in a MDOF system significantly reduces the computational cost and time compared to an implicit integrator. The results obtained from Ivanov's method in this work have been verified with the penalty method for the case of R = 1. The main advantages of the proposed method are as follows:

- 1. The nonsmooth nature of the contact problem is preserved, which is not the case with the penalty approach.
- 2. The equations in Ivanov's coordinates are non-stiff, in contrast to those obtained using the penalty approach.
- 3. Explicit integration schemes can be used with larger time-step sizes for integration purposes instead of a computationally expensive implicit integrator.

- 4. No event detection is necessary for simulating impacts, which is the most important benefit of Ivanov's approach, since the problem of event detection is challenging in the case of large MDOF systems (like Galerkin approximations and finite element analysis of continuous structures) subjected to distributed collisions.

This paper is organized as follows: In Sec. 2.1, we describe in detail how Ivanov's transformation works by modeling a point mass falling on a rigid obstacle. Using the example of a vibro-impacting SDOF system, we compare the advantages of Ivanov's method over the penalty method in Sec. 2.2. The equations of motion for a nonlinear MDOF system in Ivanov's coordinates are derived in Sec. 2.3. In Sec. 3, we discuss the Galerkin approximation of the governing differential equation of a nonlinear string and the modal-to-physical coordinate transformation. In Sec. 4, we apply Ivanov's transformation to simulate the vibro-impact motion of the string against a distributed obstacle. The results are validated and compared with the penalty approach. The contributions of the paper have been summarized in Sec. 5.

## 2 Mathematical modeling

In this section, using the example of a point mass bouncing on a rigid surface, we illustrate the idea behind Ivanov's transformation. Later, by using the example of a vibro-impacting oscillator, we demonstrate the advantages of modeling impact using Ivanov's method over the penalty approach. When simulating rigid collisions using a penalty approach, the penalty stiffness term is usually selected to be large, which results in stiff differential equations. Due to the stiff nature of the differential equations, one must use a very small time-step size during numerical integration. However, in Ivanov's method, the obtained differential equations in the transformed coordinates are non-stiff. Unlike the penalty approach, Ivanov's method accurately captures the nonsmooth behavior in velocity during impact.

## 2.1 Point mass bouncing on a rigid surface

To illustrate Ivanov's transformation, we consider the problem of a point mass, under the influence of gravity, impacting a rigid surface (see Fig. 1). The equation of motion for this problem can be written as follows:

$$\ddot { p } = - g, \ p \geq 0,$$

with initial conditions p (0) = α 0 and ˙(0) p = β 0 . When the mass makes contact with the obstacle ( p t ( c ) = 0) at time t c , the following jump condition in velocity must be respected:

$$\dot { p } ( t _ { c } ^ { + } ) = - R \dot { p } ( t _ { c }$$

where R is the coefficient of restitution. Equations (1) and (2) completely describe the motion of the mass. To apply Ivanov's transformation, we recast Eqs. (1) and (2) in state-space form. Introducing u = p and v = ˙, we get p

$$u ( t _ { c } ) = 0, \, \text{we have} \quad & \quad \quad \quad$$

 and when u t ( c ) = 0, we have

$$v ( t _ { c } ^ { + } ) = - R v ( t _ { c } ^ { - } )$$

Ivanov introduced the following nonsmooth transformation:

Figure 1: Schematic of a point mass falling on a rigid surface.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00486v3-figure-3f6ca0e0-5416-4629-aded-b52a942c3fef.jpg)

$$u = \eta \, \text{sgn} ( \eta )$$

$$v = \zeta \left ( 1 - k \, \text{sgn} ( \eta \zeta ) \right ) \text{sgn} ( \eta )$$

In Eq. (6), k = 1 -R 1+ R . Equation (5) ensures that u ≥ 0 for all time. To understand the transformation given by Eq. (6), we need to look at the phase-space of η -ζ , as shown in Fig. 2(b); the actual phase-space (in u -v ) is shown in Fig. 2(a) for comparison. In the first and third quadrants of the η -ζ plane, we have sgn( ηζ ) &gt; 0. In the second and fourth quadrants of the η -ζ plane, we have sgn( ηζ ) &lt; 0. When a trajectory AB -starting in phase-space η -ζ from point A is about to cross η = 0 at time t -B from the fourth to the third quadrant, from Eq. (6), we have

$$v ( t _ { B } ^ { - } ) = ( 1 + k ) \, \zeta ( t _ { B } ^ { - } )$$

Immediately after the contact, at time t + B , we have

$$v ( t _ { B } ^ { + } ) = - \left ( 1 - k \right ) \zeta ( t _ { B } ^ { + } )$$

Similarly, when the trajectory B C + -is about to cross η = 0 at time t -C from the second to the first quadrant, we have

$$v ( t _ { C } ^ { - } ) = - \left ( 1 + k \right ) \zeta ( t _ { C } ^ { - } )$$

After crossing η = 0, we have

$$v ( t _ { C } ^ { + } ) = ( 1 - k ) \, \zeta ( t _ { C } ^ { + } )$$

It should be noted that both η and ζ are continuous, so we have ζ ( t -B ) = ζ ( t + B ) and ζ ( t -C ) = ζ ( t + C ), which results in the following expressions:

$$v ( t _ { B } ^ { + } ) = - \left ( \frac { 1 - k } { 1 + k } \right ) v ( t _ { B } ^ { - } ) = - R v ( t _ { B } ^ { - } ) \\.. \quad & / 1 - k \, \cdots \,. \quad. \quad. \quad.$$

$$v ( t _ { C } ^ { + } ) = - \left ( \frac { 1 - k } { 1 + k } \right ) v ( t _ { C } ^ { - } ) = - R v ( t _ { C } ^ { - } ) \right ) _ { \begin{array} {$$

Therefore, it is clear that the nonsmooth transformation given by Eq. (6) automatically satisfies Eq. (4) at the event of an impact. Due to the jump conditions (Eq. (11) and Eq. (12)) imposed on the velocity v at the time of impact, the trajectories in η -ζ space that enter the third quadrant from the fourth quadrant can only go to the first quadrant through the second quadrant. Equations (5) and (6) can be written as follows:

Figure 2: (a) Phase space in original u -v coordinates. (b) Phase space in transformed η -ζ coordinates.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00486v3-figure-26b6880a-44a6-4043-91f1-227b42e28c75.jpg)

$$\left \{ \begin{array} { c } u \\ v \end{array} \right \} = S \left \{ \begin{array} { c } \eta \\ \zeta \end{array} \right \} \\ ; \, S \text{ in Eq. (13), is given by}$$

   The transformation matrix, S in Eq. (13), is given by

$$S = \left [ \begin{array} { c c } \text{sgn} ( \eta ) & 0 \\ 0 & ( 1 - k \text{ sgn} ( \eta \zeta ) ) \text{sgn} ( \eta ) \end{array} \right ] \quad & ( 1 4 )$$

The matrix S is constant in each quadrant of phase-space (see Fig. 2). Substituting Eq. (13) into Eq. (3), we get

$$\left \{ \begin{array} { c } \dot { \eta } \\ \dot { \zeta } \end{array} \right \} = S ^ { - 1 } \left [ \begin{array} { c c } 0 & 1 \\ 0 & 0 \end{array} \right ] S \left \{ \begin{array} { c } \eta \\ \$$

    The initial conditions for η and ζ can be calculated as follows:

$$\left \{ \begin{array} { c } \eta ( 0 ) \\ \zeta ( 0 ) \end{array} \right \} = S ^ { - 1 } \left \{ \begin{array} { c } u ( 0 ) \\ v ( 0 ) \end{array} \right \}$$

    Equation (15) can be integrated numerically to obtain η and ζ . Later, Eq. (13) can be used to calculate u and v . It should be noted that there are no constraints on Eq. (15) and its solutions are continuous functions of time. Figure 3(a) shows the displacement of the point mass in transformed coordinates ( η , solid blue line) and physical coordinates ( u , dashed red line). Similarly, Fig. 3(b) shows the velocity in transformed coordinates ( ζ , solid blue line) and physical ( v , dashed red line) coordinates. Figure 3(c) shows the total energy of the point mass, and it can be seen that energy is lost at every collision. If we integrate Eq. (15) for a sufficiently long time, the energy will approach zero. Also, the frequency of the transformed solution ( η ) continuously increases due to the chattering nature of the physical solution ( u ). If we use an adaptive time-step integrator, the integrator automatically reduces the time-step sizes to account for increasing frequency in the solution. For practical reasons, if the energy falls below a small threshold, we can assume that the point mass has reached equilibrium.

## 2.2 Single-degree-of-freedom vibro-impacting system

To demonstrate the advantage of modeling impact using Ivanov's method over a penalty approach, we consider the problem of a SDOF vibro-impacting system, as shown in Fig. 4(a). In the penalty approach, the rigid obstacle as shown in Fig. 4(a) is replaced with a high-stiffness ( k p ) linear spring (see Fig. 4(b)), usually referred to as the penalty stiffness. For a large value of k p , the spring behaves approximately like a rigid obstacle. In a penalty approach, we solve the following equation:

$$m \ddot { p } + c \dot { p } + ( k _ { s } + \mu k _ { p } ) p = 0$$

In Eq. (17), µ is defined as follows:

$$\mu = \begin{cases} 0 & p \geq 0 \\ 1 & p < 0 \\ & 7 \end{cases}$$

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00486v3-figure-0145e300-8e7d-446a-b15a-1b2c263c4183.jpg)

Figure 3: (a) Position of the point mass in physical coordinates ( u ) and transformed coordinates ( η ). (b) The velocity of the point mass in physical ( v ) and transformed ( ζ ) coordinates. (c) Energy of the point mass over time. The results were obtained for u (0) = 1, v (0) = 0, g = 9 8, . and R = 0 9. .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00486v3-figure-171b86f1-c8f5-4e1e-9b78-e6258c9c30ef.jpg)

(b)

Figure 4: (a) Schematic of a spring-mass system with an obstacle. (b) Equivalent representation of the spring-mass system with the rigid obstacle replaced by a linear spring.

From Eq. (17), we can see that the system has a time-scale of τ 1 = 2 π √ m k s when the mass is not in contact with the obstacle ( µ = 0); when the mass is in contact with the obstacle ( µ = 1), it has a time-scale of τ 2 = 2 π √ m k s + k p . Having two separate time-scales of very different magnitudes makes Eq. (17) a stiff differential equation. Therefore, when an explicit integrator is used to solve Eq. (17), one must use a time-step size much smaller than τ 2 (the smallest time-scale present in the solution) to obtain accurate results. The equation of motion of the system, considering a perfectly rigid obstacle ( R = 1) (see Fig. 4(a)), can be written as follows:

$$m \ddot { p } + c \dot { p } + k _ { s } p = 0, \ \ p \geq d$$

To apply Ivanov's method for solving the vibro-impact problem (Eq. (19)), we substitute u = p -d and v = ˙ in Eq. (19) to obtain p

$$\begin{array} { c } \ x v - p \, \min \, \frac { 1 } { 2 } \, \colon \, \colon \, \colon \, \colon \\ \quad \left \{ \begin{array} { c } \dot { u } \\ \dot { v } \end{array} \right \} = \left [ \begin{array} { c c } 0 & 1 \\ - k _ { s } / m & - c / m \end{array} \right ] \left \{ \begin{array} { c } u \\ v \end{array} \right \} \\ \quad \ + \left \{ \begin{array} { c } 0 \\ - k _ { s } d / m \end{array} \right \}, \ u \geq 0 \\ \end{array} \right \} \\ \text{the same transformation as discussed in the point mass bouncing problem (see}$$

  By following the same transformation as discussed in the point mass bouncing problem (see Eqs. (5) and (6)), we can write the equations of motion for the vibro-impacting system with perfectly rigid impacts in Ivanov's coordinates as follows:

$$\underset { \begin{array} { c } \text{gau mpaocs} m \text{ favour} \, \, \text{s corunacos as novo.} \end{array} } { \left \{ \begin{array} { c } \dot { \eta } \\ \dot { \zeta } \end{array} \right \} } = S ^ { - 1 } \left [ \begin{array} { c c } 0 & 1 \\ - k _ { s } / m & - c / m \end{array} \right ] S \left \{ \begin{array} { c } \eta \\ \zeta \end{array} \right \} } \\ + S ^ { - 1 } \left \{ \begin{array} { c } 0 \\ - k _ { s } d / m \end{array} \right \} } \quad ( 2 1 ) \\ \text{Eq. (14)} \text{ for the definition of the matrix $S$.)} \text{ } As pointed out earlier (see Fig. 2).}$$

  (Refer to Eq. (14) for the definition of the matrix S .) As pointed out earlier (see Fig. 2), the matrix S is diagonal and constant in each quadrant of the phase-space. Further, the magnitude of the elements of S are of the same order in all quadrants (see Fig. 2). Therefore, when Eq. (21) is solved using an explicit integrator, it does not appear as a stiff equation to the integrator.

First, we determine the solution of the system shown in Fig. 4(a) subject to the constraint p ≥ 0 5 using the event detection technique. . The solution of Eq. (19) is obtained using the 'ode45' (adaptive time-step) integrator in MATLAB with relative and absolute tolerances of 1 × 10 -12 . The parameters considered for the analysis are m = 1, k = 1, c = 0, and the integration has been carried out over a time of t = 0 to t = 10. The initial conditions are

p (0) = 1 and ˙(0) p = 0. This reference solution p r ( ) t is used to compare the results from Ivanov's method and the penalty method with an explicit fixed-time-step-size integrator. Now, we solve the same problem (Fig. 4) using the penalty approach (Eq. (17)) and Ivanov's method (Eq. (21)) using a fixed-time-step-size, fourth-order Runge-Kutta integrator. Different timestep sizes (∆ ) in the range of 10 t -4 to 10 -1 are considered for the analysis. Further, the mean squared error (MSE) given by e = 1 n ∑ n i =1 ( p r ( t i ) -p t ( i )) 2 is evaluated between the reference solution p r ( ) and the solutions obtained from Ivanov's method and the penalty method ( t p t ( )) for different time-step sizes of integration. Here, n is the number of sample points considered for the error estimation and t i is the time instant. The variation of MSE for the two methods with respect to the change in integration time-step size is shown in Fig. 5. In Fig. 5, e is calculated by evaluating the solutions ( p t ( ) and p r ( )) at t n = 101 equally spaced time intervals between t = 0 and t = 10. The MSE for the penalty approach is calculated for 3 different penalty stiffness values ( k p = 10 5 , 10 6 , and 10 ) and the corresponding variation in 7 e is shown in Fig. 5. From the results shown in Fig. 5, it is clearly observed that the MSE in the penalty method increases with an increase in ∆ . t This is because of the penalty stiffness k p in Eq. (17), which decreases the fastest time-scale τ 2 and makes the equation stiff. Therefore, as discussed earlier, when an explicit integrator is used for the integration, ∆ t must be chosen much smaller than τ 2 . As a result, in the penalty method, the suitable time-step size for integration decreases with an increase in the penalty stiffness k p , which can be clearly observed from Fig. 5. Also, with the increase in k p , the MSE in the solution decreases for the smaller integration time-step sizes (Fig. 5). Therefore, to obtain an accurate solution in the penalty method, a very large value of k p must be chosen. However, this decreases τ 2 and a correspondingly smaller time-step size must be used for integration, which increases the computation cost significantly. In contrast, the final equations obtained in Ivanov's method are not stiff (see Eq. (21)) and hence a larger time-step size can be used for integration without reducing the accuracy of the solution. We can clearly see this in Fig. 5: the error in the solution obtained from Ivanov's method is of an order smaller than 10 -4 even when ∆ t = 10 -1 , which is not the case with the penalty approach.

In summary, with the penalty approach, to obtain an accurate solution of a vibro-impact system, the time-step size must chosen to be very small depending on the penalty stiffness k p . However, by using Ivanov's method, an accurate solution can be obtained with a much larger time-step size. Moreover, for simulating a rigid collision, one must use a large value for the penalty stiffness, which puts a hard restriction on the numerical time-step size of the integrator. However, in Ivanov's approach, the nonsmooth transformation simulates the case

of infinite penalty stiffness without making the equations stiff.

Figure 5: Variation of mean squared error ( e = 1 n ∑ n i =1 ( p r ( t i ) -p t ( i )) 2 ) with respect to the time-step size of integration (∆ ) in Ivanov's method and the penalty method. t

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00486v3-figure-c0aa5d12-b481-4c0d-81e3-cfaaff31532c.jpg)

## 2.3 Ivanov's transformation for distributed impacts

Having demonstrated the applicability and advantage of Ivanov's transformation for an SDOF system, we now discuss its extension to MDOF systems where only some degrees of freedom are subjected to impact constraints. To develop Ivanov's transformation for the most general MDOF system, we consider the following coupled nonlinear model where some degrees of freedom ( p ) are subjected to impact while the others ( q ) are not subjected to any impact. The equations of motion for such a mechanical system can be written as follows:

$$\ddot { p } = \ f ( p, \dot { p }, q, \dot { q }, t )$$

$$\ddot { q } = \ g ( p, \dot { p }, q, \dot { q }, t ),$$

with the initial conditions p ( 0 ) = α 0 , ˙( p 0 ) = β 0 , q ( 0 ) = γ 0 , and ˙( q 0 ) = ν 0 . In Eq. (22), the displacement coordinates p = [ p 1 ( ) t , p 2 ( ) t , ..., p m ( )] t T are subjected to impact constraints of the form

$$p _ { i } \geq d _ { i }, \ i = 1, 2, \dots, m$$

Once any of the above constraints becomes an equality constraint ( p i ( t c ) = d i ) at time t c , a velocity jump condition is imposed as follows:

$$\dot { p } _ { i } ( t _ { c } ^ { + } ) = - R \dot { p } _ { i } ( t _ { c } ^ { - } )$$

By defining d = [ d , d 1 2 , ..., d m ] T , Eq. (24) can be written compactly as

$$p - d \geq 0$$

It should be noted that the displacement coordinates q = [ q 1 ( ) t , q 2 ( ) t , ..., q n ( )] t T are not subjected to any impact constraints. By introducing the state variables u = p -d v , = ˙, p r = q , and s = ˙, Eqs. (22) and (23) can be written as follows: q

$$\dot { u } \ = \ v$$

$$\dot { v } \ = \ f ( u + d, v, r, s, t )$$

$$\dot { r } \ = \ s$$

$$\dot { s } \ = \ g ( u + d, v, r, s, t )$$

The impact constraint becomes u ≥ 0, where u = [ u 1 ( ) t , u 2 ( ) t , ..., u m ( )] t T . Once any of the above constraints becomes an equality constraint ( u i ( t c ) = 0) at time t c , a velocity jump condition is imposed as follows:

$$v _ { i } ( t _ { c } ^ { + } ) = - R v _ { i } ( t _ { c } ^ { - } )$$

To apply Ivanov's method, the following transformation is introduced:

$$u \ = \ T \eta$$

$$v \ = \ W \zeta$$

The matrices T and W are defined as follows:

$$T \ = \ \text{diag} \, ( \text{sgn} ( \eta ) )$$

$$W \ = \ \text{diag} \left ( ( 1 - k \, \text{sgn} ( \eta \circ \zeta ) ) \circ \text{sgn} ( \eta ) \right )$$

In Eq. (35), the symbol ' ◦ ' represents the element-by-element multiplication (Hadamard product) of two vectors. Substituting Eqs. (32) and (33) into Eqs. (27)-(30), we get

$$\dot { \eta } \ = \ T ^ { - 1 } W \zeta$$

$$\dot { \zeta } \ = \ W ^ { - 1 } f ( T \eta + d, W \zeta, r, s, t )$$

$$\dot { r } \ = \ s$$

$$\dot { s } \ = \ g ( T \eta + d, W \zeta, r, s, t )$$

The initial conditions for Eqs. (36)-(39) are

$$\eta ( 0 ) \ = \ \ T ^ { - 1 } u ( 0 ) = T ^ { - 1 } \left ( p ( 0 ) - d \right )$$

=

T

(

α

-

1

0

-

d

)

$$\zeta ( 0 ) \ = \ W ^ { - 1 } v ( 0 ) = W ^ { - 1 } \dot { p } ( 0 ) = W ^ { - 1 } \beta _ { 0 }$$

$$r ( 0 ) \ = \ q ( 0 ) = \gamma _ { 0 }$$

$$s ( 0 ) \ = \ \dot { q } ( 0 ) = \nu _ { 0 }$$

The solution of Eqs. (36)-(39) automatically satisfies all the constraints imposed on p (Eq. (24)) and ˙ p (Eq. (25)). Equations (36)-(39) are the closed-form equations for a nonlinear MDOF system (where only some degrees of freedom are subjected to constraints) in Ivanov's coordinates. Upon solving the system in Ivanov's coordinates, the physical solution (displacements and velocities) can be reconstructed by using the inverse of the transformations given by Eqs. (32) and (33).

## 3 Vibro-impact motion of a string

To demonstrate the application of the proposed methodology, we simulate the vibro-impact motion of a nonlinear string against a rigid distributed obstacle (see Fig. 6). The governing equation of motion for the string is given by [34]

$$\begin{array} { c } \rho \frac { \partial ^ { 2 } Y } { \partial T ^ { 2 } } - T _ { 0 } \left ( 1 + \frac { E A } { 2 T _ { 0 } L } \int _ { 0 } ^ { L } \left ( \frac { \partial Y } { \partial X } \right ) ^ { 2 } d X \right ) \frac { \partial ^ { 2 } Y } { \partial X ^ { 2 } } + C \frac { \partial Y } { \partial T } = 0, \\ \dots \end{array} \dots$$

with the following boundary conditions:

$$Y ( 0, T ) = Y ( L, T ) = 0$$

Here, Y is the transverse displacement of the string, L is the length, A is its cross-sectional area, T 0 is the string tension at the equilibrium position, ρ is the mass density, C is the damping coefficient, E is the Young's modulus of the material of the string, X is the coordinate along the length of the string, and T is time.

Figure 6: Schematic of a string with a distributed impact constraint.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00486v3-figure-dbc999cd-a759-4ffd-80d4-d1d1026ac3c5.jpg)

We now substitute the following non-dimensional parameters into the governing equation (Eq. (44)) to facilitate the analysis:

$$y & = \frac { Y } { H }, \ x = \frac { X } { L }, \ t = \frac { T } { \alpha }, \ \gamma = \frac { E A H ^ { 2 } } { 2 T _ { 0 } L ^ { 2$$

Upon substituting the non-dimensional parameters (Eq. (46)), the governing equation of motion (Eq. (44)) can be written in non-dimensional form as follows:

$$\frac { \partial ^ { 2 } y } { \partial t ^ { 2 } } - \left ( 1 + \gamma \int _ { 0 } ^ { L } \left ( \frac { \partial y } { \partial x } \right ) ^ { 2 } d x \$$

with the boundary conditions:

$$y ( 0, t ) = y ( 1, t ) = 0$$

Now, the solution of Eq. (47) is assumed as follows:

$$y ( x, t ) = \sum _ { j = 1 } ^ { N } \phi _ { j } ( x ) \eta _ { j } ( t ) = \phi ^ { \text{T} } ( x ) \eta ( t ) \quad \quad \quad$$

In Eq. (49), φ ( x ) = [ φ 1 ( x , ) φ 2 ( x , ..., φ ) N ( x )] T and η ( ) = [ t η 1 ( ) t , η 2 ( ) t , ..., η N ( )] t T . Here, φ j ( x ) are the mass-normalized mode-shapes of the string and η j ( ) are the modal coordinates. t In this work, φ j ( x ) are chosen to be √ 2 sin( jπx ). On substituting Eq. (49) in Eq. (47), pre-multiplying by φ ( x ), integrating over the domain [0 1], and simplifying using orthogonality conditions, the following set of coupled ODEs are obtained:

$$\ddot { \eta } _ { j } ( t ) + \left ( 1 + \gamma \sum _ { k = 1 } ^ { N } \pi ^ { 2 } k ^ { 2 } \eta _ { k } ( t ) ^ { 2 } \right ) \ome$$

Here, ω j = jπ are the natural frequencies of the string. Equation (50) can now be written in matrix form as

$$\bar { M } \ddot { \eta } ( t ) + \bar { K } \left ( \eta \left ( t \right ) \right ) \eta ( t ) + \bar { C } \dot { \eta } ( t ) = 0,$$

where ¯ M is an identity matrix and ¯ ( K η t ( )) is a diagonal matrix with the diagonal elements ¯ K jj = ( 1 + γ ∑ N k =1 π k η 2 2 k ( ) t 2 ) ω 2 j and ¯ C = c ¯ . M Here, Eq. (51) is in modal coordinates and it is difficult to implement impact conditions in the modal system. Therefore, we discretize the string into a set of physical coordinates (see Fig. 7) and introduce the following notation:

$$y ( x _ { i }, t ) = p _ { i } ( t )$$

Figure 7: Schematic of a spatially discretized string.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00486v3-figure-fbfff5a9-dcd2-464f-9b1f-53595da65042.jpg)

Using Eq. (52), Eq. (49) can be written as

$$p _ { i } ( t ) = \phi ^ { \text{T} } ( x _ { i } ) \eta ( t )$$

Equation (53) can be written in matrix form as follows:

$$p ( t ) = \Phi ^ { \text{T} } \eta ( t ),$$

where p ( ) = [ t p 1 ( ) t , p 2 ( ) t , ..., p N ( )] t T and Φ = [ φ ( x 1 ) , φ ( x 2 ) , ..., φ ( x N )] T . Now, Eq. (54) can be rewritten as

$$\eta ( t ) & = \left ( \Phi ^ { T } \right ) ^ { - 1 } p ( t ) & ( 5 5 ) \\ \intertext { into Eq. (51) and pre-multiplying by $\Phi^{\prime}$, Eq. (51) can be written }$$

Upon substituting Eq. (55) into Eq. (51) and pre-multiplying by Φ T , Eq. (51) can be written in terms of the physical coordinates as follows:

$$M \ddot { p } ( t ) + K \left ( p \left ( t \right ) \right ) p ( t ) + C \dot { p } ( t ) = 0,$$

where M = Φ T ¯ M ( Φ T ) -1 , K p ( ( )) = t Φ T ¯ K ( ( Φ T ) -1 p ( ) t ) ( Φ T ) -1 and C = Φ T ¯ C ( Φ T ) -1 .

Now, by applying Ivanov's transformation as discussed in Sec. 2.3, to Eq. (56), the motion of the string impacting a distributed rigid obstacle can be obtained.

## 4 Results and discussion

We now validate the proposed approach by simulating the vibro-impact motion of a string (see Fig. 6) against a distributed obstacle using the Galerkin-Ivanov transformation. The analysis is performed by solving Eq. (56) subjected to the impact constraint using Ivanov's method and penalty method. For the analysis, the non-dimensional parameter γ is chosen to be 1 and c is chosen to be 0 unless otherwise specified. The validation is performed by comparing the solutions obtained from Ivanov's method with the results from the penalty approach. The obstacle in the penalty method is modeled as a foundation of unilateral linear springs, each having a spring constant of k p = 1 × 10 . 8 All numerical simulations have been carried out in MATLAB using a fixed-time-step-size, fourth-order Runge-Kutta integrator with ∆ t = 1 × 10 -4 . The penalty stiffness for the string and the integration time-step size have been selected to get a good match between the responses from Ivanov's method and penalty method. It should be noted that Ivanov's method is a limiting case of the penalty approach (with infinite penalty stiffness) and it is to be expected that the results from these methods will match only for large values of penalty stiffness. To demonstrate the reliability of the proposed modal-to-physical transformation, we plotted the first 200 natural frequencies of

the linearized system before and after the transformation in Fig. 8. The red circles in Fig. 8 represent the natural frequencies of the linearized system in modal coordinates ( ω m ) (Eq. (51)) and the blue line represents the frequencies in physical coordinates ( ω s ) (Eq. (56)). This clearly demonstrates the accuracy of the transformation in conserving all the natural frequencies of the system exactly.

Figure 8: Natural frequencies of a string in modal coordinates ω m (Eq. (51)) and in physical coordinates ω s (Eq. (56)).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00486v3-figure-28a6efba-09a0-412d-a59a-5a8ca69c8eda.jpg)

We now discuss the results for the string impact problem. In the present analysis, we considered N = 201 modes/grid-points in the Galerkin approximation,for which an accurate match is obtained between Ivanov's method and the penalty approach. As a first example, we simulate the vibro-impacting motion of a string against a flat obstacle. The considered obstacle is defined as d x ( ) = 0 025 . , 1 3 ≤ x ≤ 2 3 . The initial displacement of the string is considered to be 0 05 . × sin( πx ), and the initial velocity is considered to be zero. On following Ivanov's approach (discussed in Sec. 2.3), the solution for the problem is obtained.

In Fig. 9(a), we have shown the shape of the string obtained from numerical simulations at different time instants for R = 1 (CoR) obtained using both the Ivanov and penalty approaches. A clear match between the results obtained from both the methods has been observed in Fig. 9(a). Furthermore, we have also implemented the Ivanov approach using an explicit numerical integrator for different time-step sizes (∆ t = 0 0025 . , 0 0013 . , 0 0008, and 0 0001). The . . MSE ( ( )) between the solutions obtained from Ivanov and penalty approaches with respect to e t time has been evaluated and shown in Fig. 9(b). The solution from the penalty approach has been determined for the penalty spring stiffness of k p = 10 8 and is used to compare the Ivanov solutions. From Fig. 9(b), it is clear that the error in the Ivanov solution reduces with decrease

Figure 9: (a) Snapshots of the string motion impacting a flat obstacle at different time instants ( ) t obtained with Ivanov's transformation (solid red line) and the penalty approach (dashed blue line). (b) Mean squared error ( ( )) e t in the string deflection with respect to time ( ) t between the solutions obtained with Ivanov's approach and the penalty method for different integration time-step sizes (∆ ). t

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00486v3-figure-65dd2fd0-0dc0-499e-88d4-4cf1a9d03fe6.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00486v3-figure-f83c800b-3fa3-48fe-9887-671a8642f955.jpg)

Figure 10: Midpoint displacement of the string impacting a flat obstacle obtained with Ivanov's transformation: R = 1 (solid red line), R = 0 9 (dashed red line), and . R = 0 8 (dotted red . line). The dashed blue lines show the results obtained using the penalty approach.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00486v3-figure-3fa7f7fc-b2dc-4dba-8c01-9e1eafc2992a.jpg)

in the integration time-step size and converges to the penalty solution. It is also observed that the initial MSE between the solutions is of an order less than 10 -5 and increases gradually with respect to time. This happens because of the finite penalty stiffness, due to which the duration of contact is also finite. When the simulation is run for a long time, these accumulated contact times will introduce a phase shift in the solution when compared to Ivanov's solution. Moreover, for simulating a rigid collision in penalty approach, one has to use a large value for the penalty stiffness, which puts a hard restriction on the numerical integrator time-step size. However, in Ivanov's approach, the nonsmooth transformation simulates the case of infinite penalty stiffness, but without making the equations stiff. Therefore, the comparison between the Ivanov and penalty methods will be valid only for small simulation times (few impacts). In addition, the unilateral constraints are exactly satisfied in Ivanov's approach, unlike the penalty method. We can only compare results obtained from Ivanov's method for the case of R = 1 with those obtained from the penalty approach. This is because there are no models that relate R with contact dissipation in the penalty approach for discretized continuous systems.

In Fig. 10, we show the midpoint deflection of the string for three values of CoR ( R = 1 0, . 0 9, . and 0 8) for . c = 0. Again, we can see a good agreement between the solutions for the penalty method (dashed blue line) and Ivanov's method for R = 1 (solid red line). It is also observed from Fig. 10 that the rebounding displacement decreases for a smaller R because of the loss of energy with every impact for R &lt; 1. We also present the midpoint displacement of the string impacting a flat obstacle for different damping values ( c = 0 , 0 1, . and 0 2) in . Fig. 11(a). The corresponding penalty solutions have also been shown in Fig. 11(b) with a dashed blue line and an accurate match between both the solutions has been observed. In addition, the MSE between the Ivanov and penalty solutions for different integration time-step sizes (∆ t = 0 0025 . , 0 0013 . , 0 0008 and 0 0001) for . . c = 0 2 and . R = 1 has been presented in Fig. 11(b). It is clearly observed from Fig. 11(b) that the MSE between both solutions is of the order less than 1 5 . × 10 -5 for all ∆ . t Also the Ivanov solution seems to converge to the penalty solution for decreasing ∆ . t

For the second example, we consider the motion of a string impacting a sinusoidal obstacle of the form d x ( ) = 0 05 . -0 025 . × sin( π x ( -1 3 )) , 1 3 ≤ x ≤ 4 3 . The initial displacement of the string is considered to be 0 05 . × sin( πx ), and the initial velocity is taken to be zero. Figures 12, 13, and 14 are similar to Figs. 9, 10, and 11, respectively, except that the results are presented for the case of a sinusoidal obstacle. It is clearly observed from Figs. 12, 13 and 14 that the solutions from Ivanov's method and the penalty approach match closely with an MSE of order

Figure 11: (a) Midpoint displacement of the string impacting a flat obstacle obtained with Ivanov's transformation for c = 0 0 . , 0 1 . , and 0 2. . The dashed blue lines show the results obtained using the penalty approach. (b) Error between Ivanov and penalty solutions for different integration time-step sizes (∆ ) for t c = 0 2 and . R = 1.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00486v3-figure-944ef895-5159-4b50-876c-e92b0540911c.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00486v3-figure-41d43745-b199-4ebe-8eb0-c80ed525ed75.jpg)

Figure 12: (a) Snapshots of the string motion impacting a sinusoidal obstacle obtained with Ivanov's transformation (solid red line) and the penalty approach (dashed blue line). (b) Error in the string deflection between the solutions obtained from Ivanov's approach and the penalty method.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00486v3-figure-39ece051-459d-48bb-85f3-aedee39ba2d4.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00486v3-figure-8a4fa6e7-02fe-4327-9d98-f7a166ed2add.jpg)

Figure 13: Midpoint displacement of the string impacting a sinusoidal obstacle obtained with Ivanov's transformation: R = 1 (solid red line), R = 0 9 (dashed red line), and . R = 0 8 (dotted . red line). The dashed blue lines show the results obtained using the penalty approach.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00486v3-figure-dda978e1-ee89-497f-a590-4ace88e8ac45.jpg)

less than 10 -5 .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00486v3-figure-4993f737-0ed0-46ff-a714-b23c6bb4ba85.jpg)

Figure 14: (a) Midpoint displacement of the string impacting a sinusoidal obstacle obtained with Ivanov's transformation for c = 0 0 . , 0 1 . , and 0 2. . The dashed blue lines show the results obtained using the penalty approach. (b) Error between Ivanov and penalty solutions for different integration time-step sizes (∆ ) for t c = 0 2 and . R = 1.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00486v3-figure-cc950d8e-e860-4d0e-9975-eb1b834dd4f2.jpg)

From the results shown in Figs. 9-13, it is clearly demonstrated that Ivanov's transformation can be successfully used to simulate the vibro-impacting motions of a continuous structure with distributed contacts using a Galerkin approximation. Even though we have reported results for a string, the approach can also be used to simulate impacts in membranes, beams, and plate structures.

## 5 Conclusions

In this paper, the application of Ivanov's method to simulate the nonsmooth motion of a continuous structure against a rigid distributed obstacle using Galerkin approximation is demonstrated. A closed-form expression for the nonlinear MDOF system in Ivanov's transformed coordinates is reported. The transformed differential equations are without any constraints and automatically satisfy the unilateral and velocity jump conditions at impact. The technique is validated by solving three problems: (i) a spring-mass-damper system, (ii) the motion of a string against a flat obstacle, and (iii) the motion of a string against a sinusoidal obstacle. The governing partial differential equations of motion of the string are approximated using the Galerkin approach. We first convert the ODEs obtained from the Galerkin approximation to the physical coordinates using a modal-physical transformation and then apply Ivanov's approach. The solutions that are obtained from Ivanov's method are then compared with the results from the penalty method. The equations of motion obtained from Ivanov's method are not stiff, in contrast to those obtained from the penalty method. Therefore, a larger timestep size can be used for the numerical integration of the equations of motion obtained from Ivanov's method. The event detection problem for large MDOF systems in the presence of multiple impact constraints is a difficult and challenging task and is completely eliminated by the application of Ivanov's method.

## 6 Acknowledgments

CPV gratefully acknowledges the Department of Science and Technology for funding this research through Inspire fellowship (DST/INSPIRE/04/2014/000972).

## References

- [1] V. I. Babitsky, Theory of vibro-impact systems and applications . Berlin, Germany: Springer-Verlag, 1998.
- [2] R. A. Ibrahim, Vibro-impact dynamics: modeling, mapping and applications . Berlin, Heidelberg: Springer-Verlag, 2009.
- [3] R. A. Ibrahim, 'Recent advances in vibro-impact dynamics and collision of ocean vessels,' Journal of Sound and Vibration , vol. 333, no. 23, pp. 5900-5916, 2014.

- [4] K. Vijayan, M. Friswell, H. H. Khodaparast, and S. Adhikari, 'Non-linear energy harvesting from coupled impacting beams,' International Journal of Mechanical Sciences , vol. 96, pp. 101-109, 2015.
- [5] V. Chatziioannou and M. Van Walstijn, 'Energy conserving schemes for the simulation of musical instrument contact dynamics,' Journal of Sound and Vibration , vol. 339, pp. 262279, 2015.
- [6] C. Issanchou, S. Bilbao, J.-L. Le Carrou, C. Touz´, and O. Doar´, 'A modal-based ape e proach to the nonlinear vibration of strings against a unilateral obstacle: Simulations and experiments in the pointwise case,' Journal of Sound and Vibration , vol. 393, pp. 229-251, 2017.
- [7] C. Issanchou, V. Acary, F. P´rignon, C. Touz´, and J.-L. Le Carrou, 'Nonsmooth contact e e dynamics for the numerical simulation of collisions in musical string instruments,' The Journal of the Acoustical Society of America , vol. 143, no. 5, pp. 3195-3205, 2018.
- [8] J.-H. Shin, E. Talib, and M. K. Kwak, 'Vibration of a string against multiple spring-massdamper stoppers,' Journal of Sound and Vibration , vol. 414, pp. 157-173, 2018.
- [9] A. M. Shafei and H. R. Shafei, 'Oblique impact of multi-flexible-link systems,' Journal of Vibration and Control , vol. 24, no. 5, pp. 904-923, 2018.
- [10] K. Li and A. Darby, 'Experiments on the effect of an impact damper on a multipledegree-of-freedom system,' Journal of Vibration and Control , vol. 12, no. 5, pp. 445-464, 2006.
- [11] Y. Yang and X. Wang, 'Investigation into the linear velocity response of cantilever beam embedded with impact damper,' Journal of Vibration and Control , vol. 25, no. 7, pp. 13651378, 2019.
- [12] G. Gilardi and I. Sharf, 'Literature survey of contact dynamics modelling,' Mechanism and Machine Theory , vol. 37, no. 10, pp. 1213-1239, 2002.
- [13] S. Bilbao and A. Torin, 'Numerical modeling and sound synthesis for articulated string/fretboard interactions,' Journal of the Audio Engineering Society , vol. 63, no. 5, pp. 336-347, 2015.
- [14] S. Bilbao, A. Torin, and V. Chatziioannou, 'Numerical modeling of collisions in musical instruments,' Acta Acustica United with Acustica , vol. 101, no. 1, pp. 155-173, 2015.

- [15] C. Issanchou, J.-L. Le Carrou, C. Touz´, B. Fabre, and O. Doar´, 'String/frets contacts e e in the electric bass sound: Simulations and experiments,' Applied Acoustics , vol. 129, pp. 217-228, 2018.
- [16] J. Kang, 'Calculation of Lyapunov exponents in impacted beam on distributed contact,' Journal of Sound and Vibration , vol. 431, pp. 295-303, 2018.
- [17] D. Wagg, 'Periodic sticking motion in a two-degree-of-freedom impact oscillator,' International Journal of Non-Linear Mechanics , vol. 40, no. 8, pp. 1076-1087, 2005.
- [18] X. Li, Z. Yao, and R. Wu, 'Modeling and sticking motion analysis of a vibro-impact system in linear ultrasonic motors,' International Journal of Mechanical Sciences , vol. 100, pp. 23-31, 2015.
- [19] T. K. Uchida, M. A. Sherman, and S. L. Delp, 'Making a meaningful impact: modelling simultaneous frictional collisions in spatial multibody systems,' Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences , vol. 471, no. 2177, p. 20140859, 2015.
- [20] A. Liakou, V. Deno¨l, and E. Detournay, 'Fast in-plane dynamics of a beam with unilateral e constraints,' Journal of Engineering Mechanics , vol. 143, no. 2, p. 04016116, 2016.
- [21] D. Huang, W. Xu, D. Liu, and Q. Han, 'Multi-valued responses of a nonlinear vibroimpact system excited by random narrow-band noise,' Journal of Vibration and Control , vol. 22, no. 12, pp. 2907-2920, 2016.
- [22] G. F. S. Rebou¸ cas, I. F. Santos, and J. J. Thomsen, 'Unilateral vibro-impact systemsexperimental observations against theoretical predictions based on the coefficient of restitution,' Journal of Sound and Vibration , vol. 440, pp. 346-371, 2019.
- [23] V. Acary and B. Brogliato, Numerical methods for nonsmooth dynamical systems: applications in mechanics and electronics . Springer Science &amp; Business Media, 2008.
- [24] C. Studer, 'Numerics of unilateral contacts and friction,' Vol. 47 of Lecture Notes in Applied and Computational Mechanics, Springer Verlag, Berlin Heidelberg, 2009. , 2009.
- [25] V. Acary, 'Energy conservation and dissipation properties of time-integration methods for nonsmooth elastodynamics with contact,' ZAMM-Journal of Applied Mathematics and Mechanics , vol. 96, no. 5, pp. 585-603, 2016.

- [26] A. Ivanov, 'Analytical methods in the theory of vibro-impact systems,' Journal of Applied Mathematics and Mechanics , vol. 57, no. 2, pp. 221-236, 1993.
- [27] A. Ivanov, 'Impact oscillations: linear theory of stability and bifurcations,' Journal of Sound and Vibration , vol. 178, no. 3, pp. 361-378, 1994.
- [28] I. M. Grace, R. A. Ibrahim, and V. N. Pilipchuk, 'Inelastic impact dynamics of ships with one-sided barriers. Part I: analytical and numerical investigations,' Nonlinear Dynamics , vol. 66, no. 4, pp. 589-607, 2011.
- [29] I. M. Grace, R. A. Ibrahim, and V. N. Pilipchuk, 'Inelastic impact dynamics of ships with one-sided barriers. Part II: experimental validation,' Nonlinear Dynamics , vol. 66, no. 4, pp. 609-623, 2011.
- [30] C. P. Vyasarayani, J. McPhee, and S. Birkett, 'Modeling impacts between a continuous system and a rigid obstacle using coefficient of restitution,' Journal of Applied Mechanics , vol. 77, no. 2, p. 021008, 2010.
- [31] D. Wagg and S. Bishop, 'Application of non-smooth modelling techniques to the dynamics of a flexible impacting beam,' Journal of Sound and Vibration , vol. 256, no. 5, pp. 803-820, 2002.
- [32] H. Wei, G. Li, P. Guo, and J. Zhao, 'Effect of method type on the response of continuum vibro-impact,' Shock and Vibration , vol. 2019, 2019.
- [33] M. Van Walstijn and J. Bridges, 'Simulation of distributed contact in string instruments: a modal expansion approach,' in 2016 24th European Signal Processing Conference (EUSIPCO) , pp. 1023-1027, IEEE, 2016.
- [34] D. W. Oplinger, 'Frequency response of a nonlinear stretched string,' The Journal of the Acoustical Society of America , vol. 32, no. 12, pp. 1529-1538, 1960.
Surya Samukham S. N. Khaderi C. P. Vyasarayani
Galerkin-Ivanov transformation for nonsmooth modeling of vibro-impacts in continuous structures
This work deals with the modeling of nonsmooth vibro-impact motion of a
continuous structure against a rigid distributed obstacle. Galerkin's approach
is used to approximate the solutions of the governing partial differential
equations of the structure, which results in a system of ordinary differential
equations (ODEs). When these ODEs are subjected to unilateral constraints and
velocity jump conditions, one must use an event detection algorithm to
calculate the time of impact accurately. Event detection in the presence of
multiple simultaneous impacts is a computationally demanding task. Ivanov
proposed a nonsmooth transformation for a vibro-impacting
multi-degree-of-freedom system subjected to a single unilateral constraint.
This transformation eliminates the unilateral constraints from the problem and,
therefore, no event detection is required during numerical integration. Ivanov
used his transformation to make analytical calculations for the stability and
bifurcations of vibro-impacting motions; however, he did not explore its
application for simulating distributed collisions in spatially continuous
structures. We adopt Ivanov's transformation to deal with multiple unilateral
constraints in spatially continuous structures. Also, imposing the velocity
jump conditions exactly in the modal coordinates is nontrivial and challenging.
Therefore, in this work we use a modal-physical transformation to convert the
system from modal to physical coordinates on a spatially discretized grid. We
then apply Ivanov's transformation on the physical system to simulate the
vibro-impact motion of the structure. The developed method is demonstrated by
modeling the distributed collision of a nonlinear string against a rigid
distributed surface. For validation, we compare our results with the well-known
penalty approach.