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

TL;DR

This paper introduces a novel application of Ivanov's nonsmooth transformation combined with Galerkin's method to efficiently model and simulate vibro-impact motions in continuous structures without the need for event detection.

Contribution

The work extends Ivanov's transformation to spatially continuous structures with multiple impacts, integrating it with Galerkin's approach for improved simulation accuracy and computational efficiency.

Findings

Eliminates the need for event detection in impact simulations.

Successfully models distributed impacts in a nonlinear string.

Results align well with traditional penalty methods.

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

Equations

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

$$\dot{p}(t_c^{+})=-R\dot{p}(t_c^{-}),$$

$$\left\{\begin{array}[]{c}\dot{u}\\ \dot{v}\end{array}\right\}=\left[\begin{array}[]{cc}0&1\\ 0&0\end{array}\right]\left\{\begin{array}[]{c}u\\ v\end{array}\right\}+\left\{\begin{array}[]{c}0\\ -g\end{array}\right\},\thinspace\thinspace\thinspace u\geq 0$$

$$v(t_c^{+})=-Rv(t_c^{-})$$

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

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

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

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

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

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

$$v(t_B^{+})=-\left(\frac{1-k}{1+k}\right)v(t_B^{-})=-Rv(t_B^{-})$$

$$v(t_C^{+})=-\left(\frac{1-k}{1+k}\right)v(t_C^{-})=-Rv(t_C^{-})$$

$$\left\{\begin{array}[]{c}u\\ v\end{array}\right\}=\bm{S}\left\{\begin{array}[]{c}\eta\\ \zeta\end{array}\right\}$$

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

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

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

$$m\ddot{p}+c\dot{p}+(k_s+\mu k_p)p=0$$

$$\mu =\{\begin{array}{ll}0& p\ge 0\\ 1& p<0\end{array}$$

$$m\ddot{p}+c\dot{p}+k_{s}p=0,p\ge d$$

$$\displaystyle\left\{\begin{array}[]{c}\dot{u}\\ \dot{v}\end{array}\right\}=\left[\begin{array}[]{cc}0&1\\ -k_{s}/m&-c/m\end{array}\right]\left\{\begin{array}[]{c}u\\ v\end{array}\right\}$$

$$\displaystyle+\left\{\begin{array}[]{c}0\\ -k_{s}d/m\end{array}\right\},\,~{}u\geq 0$$

$$\displaystyle\left\{\begin{array}[]{c}\dot{\eta}\\ \dot{\zeta}\end{array}\right\}=\bm{S}^{-1}\left[\begin{array}[]{cc}0&1\\ -k_{s}/m&-c/m\end{array}\right]\bm{S}\left\{\begin{array}[]{c}\eta\\ \zeta\end{array}\right\}$$

$$\displaystyle+\bm{S}^{-1}\left\{\begin{array}[]{c}0\\ -k_{s}d/m\end{array}\right\}$$

$$\ddot{\mathbf{p}}=$$

$$\ddot{\mathbf{q}}=$$

$$p_i\ge d_i,i=1,2,...,m$$

$${\dot{p}}_i(t_c^{+})=-R{\dot{p}}_i(t_c^{-})$$

$$\mathbf{p}-\mathbf{d}\ge 0$$

$$\dot{\mathbf{u}}$$

$$\dot{\mathbf{v}}$$

$$\dot{\mathbf{r}}$$

$$\dot{\mathbf{s}}$$

$$v_i(t_c^{+})=-R{v}_i(t_c^{-})$$

$$\mathbf{u}$$

Full text

1 Introduction

In many engineering applications, structures are subjected to vibro-impacting motions [1, 2, 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, 6, 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 [12, 13, 5, 14, 15, 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, 19, 20, 21, 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 velocity 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-degree-of-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-of-freedom (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, 31, 32]. To avoid such difficulties, we use a transformation to convert the modal system into its physical coordinates by discretizing the string in space [33, 6]. 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:

The nonsmooth nature of the contact problem is preserved, which is not the case with the penalty approach. 2. 2.

The equations in Ivanov’s coordinates are non-stiff, in contrast to those obtained using the penalty approach. 3. 3.

Explicit integration schemes can be used with larger time-step sizes for integration purposes instead of a computationally expensive implicit integrator. 4. 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:

[TABLE]

with initial conditions $p(0)=\alpha_{0}$ and $\dot{p}(0)=\beta_{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:

[TABLE]

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=\dot{p}$, we get

[TABLE]

and when $u(t_{c})=0$, we have

[TABLE]

Ivanov introduced the following nonsmooth transformation:

[TABLE]

In Eq. (6), $k=\frac{1-R}{1+R}$. Equation (5) ensures that $u\geq 0$ for all time. To understand the transformation given by Eq. (6), we need to look at the phase-space of $\eta-\zeta$, 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 $\eta-\zeta$ plane, we have $\text{sgn}(\eta\zeta)>0$. In the second and fourth quadrants of the $\eta-\zeta$ plane, we have $\text{sgn}(\eta\zeta)<0$. When a trajectory $AB^{-}$ starting in phase-space $\eta-\zeta$ from point $A$ is about to cross $\eta=0$ at time $t_{B}^{-}$ from the fourth to the third quadrant, from Eq. (6), we have

[TABLE]

Immediately after the contact, at time $t_{B}^{+}$, we have

[TABLE]

Similarly, when the trajectory $B^{+}C^{-}$ is about to cross $\eta=0$ at time $t_{C}^{-}$ from the second to the first quadrant, we have

[TABLE]

After crossing $\eta=0$, we have

[TABLE]

It should be noted that both $\eta$ and $\zeta$ are continuous, so we have $\zeta(t_{B}^{-})=\zeta(t_{B}^{+})$ and $\zeta(t_{C}^{-})=\zeta(t_{C}^{+})$, which results in the following expressions:

[TABLE]

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 $\eta-\zeta$ 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:

[TABLE]

The transformation matrix, $\bm{S}$ in Eq. (13), is given by

[TABLE]

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

[TABLE]

The initial conditions for $\eta$ and $\zeta$ can be calculated as follows:

[TABLE]

Equation (15) can be integrated numerically to obtain $\eta$ and $\zeta$. 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 shows the displacement of the point mass in transformed coordinates ($\eta$, solid blue line) and physical coordinates ($u$, dashed red line). Similarly, Fig. 3 shows the velocity in transformed coordinates ($\zeta$, solid blue line) and physical ($v$, dashed red line) coordinates. Figure 3 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 ($\eta$) 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:

[TABLE]

In Eq. (17), $\mu$ is defined as follows:

[TABLE]

From Eq. (17), we can see that the system has a time-scale of $\tau_{1}=2\pi\sqrt{\frac{m}{k_{s}}}$ when the mass is not in contact with the obstacle ($\mu=0$); when the mass is in contact with the obstacle ($\mu=1$), it has a time-scale of $\tau_{2}=2\pi\sqrt{\frac{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 $\tau_{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), can be written as follows:

[TABLE]

To apply Ivanov’s method for solving the vibro-impact problem (Eq. (19)), we substitute $u=p-d$ and $v=\dot{p}$ in Eq. (19) to obtain

[TABLE]

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:

[TABLE]

(Refer to Eq. (14) for the definition of the matrix $\bm{S}$.) As pointed out earlier (see Fig. 2), the matrix $\bm{S}$ is diagonal and constant in each quadrant of the phase-space. Further, the magnitude of the elements of $\bm{S}$ are of the same order in all quadrants (see Fig. 2). Therefore, when Eq. (39) 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 subject to the constraint $p\geq 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\times 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 $\dot{p}(0)=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. (39)) using a fixed-time-step-size, fourth-order Runge–Kutta integrator. Different time-step sizes ($\Delta t$) in the range of $10^{-4}$ to $10^{-1}$ are considered for the analysis. Further, the mean squared error (MSE) given by $e=\frac{1}{n}\sum_{i=1}^{n}(p_{r}(t_{i})-p(t_{i}))^{2}$ is evaluated between the reference solution $p_{r}(t)$ and the solutions obtained from Ivanov’s method and the penalty method ($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}(t)$) at $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^{7}$) and the corresponding variation in $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 $\Delta t$. This is because of the penalty stiffness $k_{p}$ in Eq. (17), which decreases the fastest time-scale $\tau_{2}$ and makes the equation stiff. Therefore, as discussed earlier, when an explicit integrator is used for the integration, $\Delta t$ must be chosen much smaller than $\tau_{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 $\tau_{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. (39)) 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 $\Delta 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.

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 ($\bm{p}$) are subjected to impact while the others ($\bm{q}$) are not subjected to any impact. The equations of motion for such a mechanical system can be written as follows:

[TABLE]

with the initial conditions $\bm{p}(\bm{0})=\bm{\alpha}_{0}$, $\dot{\bm{p}}(\bm{0})=\bm{\beta}_{0}$, $\bm{q}(\bm{0})=\bm{\gamma}_{0}$, and $\dot{\bm{q}}(\bm{0})=\bm{\nu}_{0}$. In Eq. (40), the displacement coordinates $\bm{p}=[p_{1}(t),\thinspace p_{2}(t),...,p_{m}(t)]^{\text{T}}$ are subjected to impact constraints of the form

[TABLE]

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:

[TABLE]

By defining $\bm{d}=[d_{1},d_{2},...,d_{m}]^{\text{T}}$, Eq. (42) can be written compactly as

[TABLE]

It should be noted that the displacement coordinates $\bm{q}=[q_{1}(t),\thinspace q_{2}(t),...,q_{n}(t)]^{\text{T}}$ are not subjected to any impact constraints. By introducing the state variables $\bm{u}=\bm{p}-\bm{d}$, $\bm{v}=\dot{\bm{p}}$, $\bm{r}=\bm{q}$, and $\bm{s}=\dot{\bm{q}}$, Eqs. (40) and (41) can be written as follows:

[TABLE]

The impact constraint becomes $\bm{u}\geq 0$, where $\bm{u}=[u_{1}(t),\thinspace u_{2}(t),...,u_{m}(t)]^{\text{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:

[TABLE]

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

[TABLE]

The matrices $\bm{T}$ and $\bm{W}$ are defined as follows:

[TABLE]

In Eq. (53), the symbol “$\circ$” represents the element-by-element multiplication (Hadamard product) of two vectors. Substituting Eqs. (50) and (51) into Eqs. (45)-(48), we get

[TABLE]

The initial conditions for Eqs. (54)–(57) are

[TABLE]

The solution of Eqs. (54)–(57) automatically satisfies all the constraints imposed on $\bm{p}$ (Eq. (42)) and $\bm{\dot{p}}$ (Eq. (43)). Equations (54)–(57) 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. (50) and (51).

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]

[TABLE]

with the following boundary conditions:

[TABLE]

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, $\rho$ 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.

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

[TABLE]

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

[TABLE]

with the boundary conditions:

[TABLE]

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

[TABLE]

In Eq. (67), $\bm{\phi}(x)=\left[\phi_{1}(x),\,\ \phi_{2}(x),...,\phi_{N}(x)\right]^{\text{T}}$ and $\bm{\eta}(t)=\left[\eta_{1}(t),\,\ \eta_{2}(t),...,\eta_{N}(t)\right]^{\text{T}}$. Here, $\phi_{j}(x)$ are the mass-normalized mode-shapes of the string and $\eta_{j}(t)$ are the modal coordinates. In this work, $\phi_{j}(x)$ are chosen to be $\sqrt{2}\sin(j\pi x)$. On substituting Eq. (67) in Eq. (65), pre-multiplying by $\bm{\phi}(x)$, integrating over the domain $\left[0~{}~{}1\right]$, and simplifying using orthogonality conditions, the following set of coupled ODEs are obtained:

[TABLE]

Here, $\omega_{j}=j\pi$ are the natural frequencies of the string. Equation (68) can now be written in matrix form as

[TABLE]

where $\bar{\bm{M}}$ is an identity matrix and $\bar{\bm{K}}(\eta(t))$ is a diagonal matrix with the diagonal elements $\bar{K}_{jj}=\left(1+\gamma\sum_{k=1}^{N}\pi^{2}k^{2}\eta_{k}(t)^{2}\right)\omega_{j}^{2}$ and $\bar{\bm{C}}=c\bar{\bm{M}}$. Here, Eq. (69) 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:

[TABLE]

Using Eq. (70), Eq. (67) can be written as

[TABLE]

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

[TABLE]

where $\bm{p}(t)=\left[p_{1}(t),p_{2}(t),...,p_{N}(t)\right]^{\text{T}}$ and $\bm{\Phi}=\left[\bm{\phi}(x_{1}),\bm{\phi}(x_{2}),...,\bm{\phi}(x_{N})\right]^{\text{T}}$. Now, Eq. (72) can be rewritten as

[TABLE]

Upon substituting Eq. (73) into Eq. (69) and pre-multiplying by $\bm{\Phi}^{\text{T}}$, Eq. (69) can be written in terms of the physical coordinates as follows:

[TABLE]

where ${\bm{M}}=\bm{\Phi}^{\text{T}}\bar{\bm{M}}\left(\bm{\Phi}^{\text{T}}\right)^{-1}$, ${\bm{K}}\left(\bm{p}\left(t\right)\right)=\bm{\Phi}^{\text{T}}\bar{\bm{K}}\left(\left(\bm{\Phi}^{\text{T}}\right)^{-1}\bm{p}\left(t\right)\right)\left(\bm{\Phi}^{\text{T}}\right)^{-1}$ and ${\bm{C}}=\bm{\Phi}^{\text{T}}\bar{\bm{C}}\left(\bm{\Phi}^{\text{T}}\right)^{-1}$. Now, by applying Ivanov’s transformation as discussed in Sec. 2.3, to Eq. (74), 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. (74) subjected to the impact constraint using Ivanov’s method and penalty method. For the analysis, the non-dimensional parameter $\gamma$ is chosen to be $1$ and $c$ is chosen to be [math] 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\times 10^{8}$. All numerical simulations have been carried out in MATLAB using a fixed-time-step-size, fourth-order Runge–Kutta integrator with $\Delta t=1\times 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 (${\omega}_{m}$) (Eq. (69)) and the blue line represents the frequencies in physical coordinates (${\omega}_{s}$) (Eq. (74)). This clearly demonstrates the accuracy of the transformation in conserving all the natural frequencies of the system exactly.

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,~{}\frac{1}{3}\leq x\leq\frac{2}{3}$. The initial displacement of the string is considered to be $0.05\times\sin(\pi 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, 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. Furthermore, we have also implemented the Ivanov approach using an explicit numerical integrator for different time-step sizes ($\Delta t=~{}0.0025,~{}0.0013,~{}0.0008$, and $0.0001$). The MSE ($e(t)$) between the solutions obtained from Ivanov and penalty approaches with respect to time has been evaluated and shown in Fig. 9. 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, it is clear that the error in the Ivanov solution reduces with decrease 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<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. The corresponding penalty solutions have also been shown in Fig. 11 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 ($\Delta t=0.0025,~{}0.0013,~{}0.0008$ and $0.0001$) for $c=0.2$ and $R=1$ has been presented in Fig. 11. It is clearly observed from Fig. 11 that the MSE between both solutions is of the order less than $1.5\times 10^{-5}$ for all $\Delta t$. Also the Ivanov solution seems to converge to the penalty solution for decreasing $\Delta 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\times\sin(\pi(x-\frac{1}{3})),~{}\frac{1}{3}\leq x\leq\frac{4}{3}$. The initial displacement of the string is considered to be $0.05\times\sin(\pi 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 less than $10^{-5}$.

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