<figure><img src="image_1.png" /></figure>

HAL
open science

# Nonsmooth modal analysis of a two-bar system via Boundary Element Method

Tianzheng Lu, Mathias Legrand

► To cite this version:

Tianzheng Lu, Mathias Legrand. Nonsmooth modal analysis of a two-bar system via Boundary Element Method. ASME 2021 International Design Engineering Technical Conferences & Computers and Information in Engineering Conference (IDETC/CIE 2021), Aug 2021, (Virtual), United States. hal-04029283

HAL Id: hal-04029283

https://hal.science/hal-04029283

Submitted on 14 Mar 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
EXTENDED ABSTRACT: NONSMOOTH MODAL ANALYSIS OF A TWO-BAR
SYSTEM VIA BOUNDARY ELEMENT METHOD

Tianzheng Lu

Structural Dynamics and Vibration Laboratory
Department of Mechanical Engineering
McGill University
Montreal, QC, Canada H3A 0C3

Email: tianzheng.lu@mail.mcgill.ca

Mathias Legrand

Structural Dynamics and Vibration Laboratory
Department of Mechanical Engineering
McGill University
Montreal, QC, Canada H3A 0C3

Email: mathias.legrand@mcgill.ca

INTRODUCTION

Within the framework of structural dynamics, modal analysis is a practical tool to predict vibrations in order to prevent or utilize them. As a challenging topic in the field of modal analysis, nonlinear modal analysis is focused by many researchers [1]. Nonsmooth modal analysis is a special form of nonlinear modal analysis, in which governing equation of targeting system involves indifferentiable functions of the state. To perform non-smooth modal analysis, a method combing Frequency Domain Boundary Element Method (FD-BEM) and Harmonic Balance Method (HBM) is applied.

<figure><img src="image_3.png" /><figcaption>FIGURE 1. System of interest</figcaption></figure>

SYSTEM OF INTEREST

The system of interest embeds two one-dimensional bar with same cross-section area facing each other through a common unilateral contact interface, as shown in Figure 1. The system is non-dimensionalized respect to the first bar, as shown in Appendix. The wave equation governs dynamics within each bar. The first bar's boundary at $x = 0$ is attached to a linear spring, thus linear robin boundary condition is applied. The second bar's boundary at $x = 0$ is clamped thus homogeneous Dirichlet condition is applied. The unilateral contact interface between two bars are described by Signorini boundary condition. By defining the Signorini residual $r(u(1,t), u_x(1,t), w(l,t)) = u_x(1,t) + \max[\rho(u(1,t) + w(l,t) - g_0) - u_x(1,t), 0]$, where $\rho$ is

an arbitrary strictly positive constant, Signorini boundary condition be expressed as a combination of the contact force equilibrium $u_x(1,t) = \alpha w_y(l,t)$ and the equality [1]

$$r(u(1,t), u_x(1,t), w(l,t)) = 0. \quad (1)$$

SOLUTION METHOD

The FD-BEM, as its name implies, is a frequency-domain form of BEM which is appropriate when periodic solutions are targeted. In this section, the derivation process is only shown on bar 1 to be concise. The Fourier Transform along time of the displacement $u(x,t)$ (an equivalent definition holds for $w(y,t)$) leads to the well-known one-dimensional autonomous Helmholtz
equations

$$ \hat{u}_{xx}(x, \omega) + \omega^2 \hat{u}(x, \omega) = 0, \quad x \in ]0, 1[ \qquad (2) $$

respectively, where $κ = ω/c$ is the frequency number for the second bar. Then Holmholtz equation is transformed into a weighted residual form, where the weight function is the fundamental solution to the Helmholtz equation [2]. Via integral by parts, targeted boundary integral equation (BIE) is achieved for bar 1:

$$ \begin{aligned} 2\hat{u}(x, \omega) ={}& \hat{u}(0, \omega) \cos \omega x + \hat{u}(1, \omega) \cos(\omega(1-x)) \\ & - \frac{1}{\omega}(\hat{p}(0, \omega) \sin \omega x - \hat{p}(1, \omega) \sin(\omega(1-x))) \end{aligned} \qquad (3) $$

In Equation (3), $\hat{p}(0, \omega) = -\hat{u}_x(0, \omega)$ and $\hat{p}(1, \omega) = \hat{u}_x(1, \omega)$ were used to follow the traditional notation in BEM, in the context of linear elasticity and small strains considered in the present work. The same time-domain conventions are used in the remainder. The same applies to Reading (3) on the boundary $\{0\} \cup \{1\}$ leads to the two linearly independent equations for the first bar

$$ \begin{bmatrix} \omega & 0 & -\omega \cos \omega & \sin \omega \\ -\omega \cos \omega & \sin \omega & \omega & 0 \end{bmatrix} \begin{pmatrix} \hat{u}(0, \omega) \\ \hat{p}(0, \omega) \\ \hat{u}(1, \omega) \\ \hat{p}(1, \omega) \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \end{pmatrix}, \quad (4) $$

**Periodicity** Since periodic solution are targeted, boundary terms $u(1,t)$, $p(1,t)$ are assumed in the forms of two distinct Fourier series with a common frequency $\Omega$

$$ \begin{aligned} p(1,t) &= \frac{1}{2}a_0 + \sum_{n=1}^{\infty} a_n \exp(jn\Omega t) \\ u(1,t) &= \frac{1}{2}b_0 + \sum_{n=1}^{\infty} b_n \exp(jn\Omega t) \end{aligned} \qquad (5) $$

where the complex coefficients $a_n$ and $b_n$ are the new unknowns of the problem.

**Boundary Condition** Considering Robin boundary condition in frequency domain $k\hat{u}(0, \omega) + \hat{p}(0, \omega) = 0$ as well as periodicity assumption, Equation 4 simplifies to

$$ \omega_n(\omega_n - k \cot \omega_n)b_n + (k + \omega_n \cot \omega_n)a_n = 0 \quad \text{with} \quad \omega_n = n\Omega. \qquad (6) $$

Another equation set for bar 2 could be derived by similar approach:

$$ -\kappa_n \cos \kappa_n \ell f_n + \sin \kappa_n \ell d_n = 0 \quad \text{with} \quad \kappa_n = n\Omega/c. \qquad (7) $$

where complex coefficients $d_n$ and $f_n$ are discrete Fourier coefficients for $q(\ell, t)$ and $w(\ell, t)$. $q(\ell, t)$ is normal strain on the Signorini boundary in bar 2, similar to $p(1, t)$.

Unlike the above Dirichlet and Robin boundary conditions which can be explicitly expressed in terms of the unknown Fourier coefficients, the Signorini condition has no explicit form in the frequency domain. Instead, a numerical version of the Harmonic Balance Method is performed on (1) into the system of nonlinear implicit equations with coefficients $(a_n, b_n, f_n)$, $n = 0, 1, 2, \dots$ (gathered in vectors $\mathbf{a}, \mathbf{b}, \mathbf{f}$)

$$ g_n(\mathbf{a}, \mathbf{b}, \mathbf{f}) = \int_0^T \exp(jn\Omega t) r(\mathbf{a}, \mathbf{b}, \mathbf{f}) dt = 0, \quad n = 0, 1, 2, \dots \qquad (8) $$

where $T = 2\pi/\Omega$ is period of targeted periodic motion. In other words, the Signorini condition is satisfied in a weighted residual sense only, as achieved in any Galerkin-like strategy. Along with contact force equilibrium

$$ a_n = \alpha d_n, \quad \forall n \qquad (9) $$

Signorini condition is enforced.

**Discretization and numerical approximation** Discretization comes into the proposed solution strategy when the Fourier expansions (5) are truncated to a finite number $m$ of harmonics, such that we define $u^{(m)}(1,t) \approx u(1,t)$ and $p^{(m)}(1,t) \approx p(1,t)$. A second level of discretization lies in the computation of the integrals (8). It was found that the computed solutions were not sensitive to that numerical aspect. Accordingly, it was decided to compute the integrals (8) via a simple Riemann sum approximation with a subinterval $\Delta t = T/(n_h m)$ where $T$ is targeted motion period, and $n_h$, a coefficient governing the accuracy of the approximation, and chosen as $n_h = 30$ after a convergence check.

The system of nonlinear equations including Equations (6), (7), (8) and (9) is then solved numerically using a trust-region dogleg [3] solver, a built-in numerical solver of Matlab®, even though the equations are not expected to be sufficiently smooth, due to (8).

**Continuation** Once discretization is achieved, the task of finding periodic solutions translates into an equivalent multidimensional root finding problem of the form $F(\mathbf{a}, \mathbf{b}, \mathbf{d}, \mathbf{f}, \Omega) = \mathbf{0}$, where $\Omega$ is the unknown fundamental frequency of the Fourier series. To search for continuous families of periodic solutions and construct the desired solution branches, continuation techniques shall be implemented. In this work, the classical sequential continuation technique [4] is used where $\Omega$ is successively
<figure><img src="image_4.png" /><figcaption>FIGURE 2. Backbone curve of the first mode</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIGURE 3. Strain (dashed) and position (solid) at the contact interface: bar 1 [---] and bar 2 [---].</figcaption></figure>

increased by a small given increment on a given interval of interest: the nonlinear system is solved for the Fourier coefficients only.

## PRELIMINARY RESULT

In this section, some preliminary result was found with setting $\ell = 2.2$, $c = 1$, $\alpha = 1$, $k = 0.5$. With such settings, the natural frequency of each bar are in vicinity of each other: natural frequency of bar 1 is $\omega_1^{(1)} \approx 0.653 \text{ rad}^{-1}$ and natural frequency of bar 2 is $\omega_1^{(1)} \approx 0.714 \text{ rad}^{-1}$. Harmonics are truncated at $m = 20$ in this simulation.

A family of periodic solution was found in range $[\omega_1^{(1)}, 1.46\omega_1^{(1)}]$. The backbone curve is shown in Figure 2. In range $[\omega_1^{(1)}, \omega_1^{(2)}]$, hardening mode 1 of bar 1 interact with softening mode 1 of bar 2. In range $[\omega_1^{(2)}, 1.46\omega_1^{(1)}]$, both bars exhibits hardening mode 1. One of the solution is shown in figure 3.

## CONCLUSION

An combination between FD-BEM and HBM is extended to two-bar unilateral contact case. Its capability has been proved by some preliminary results. It worth noting although solutions has been found, detailed analysis of nonlinear/nonsmooth phenomena, such as subharmonic and internal resonance, are still to be performed.

## REFERENCES

[1] Stewart, D. E., 2011. *Dynamics with Inequalities: impacts and hard constraints*, Vol. 59. SIAM.

[2] Kythe, P., 2012. *Fundamental solutions for differential operators and applications*. Springer Science & Business Media.

[3] Powell, M. J., 1968. A fortran subroutine for solving systems of nonlinear algebraic equations.

[4] Nayfeh, A. H., and Balachandran, B., 2008. *Applied nonlinear dynamics: analytical, computational, and experimental methods*. John Wiley & Sons.

## Appendix A: Non-dimensionalization of the system

The non-dimensional variables are introduced via an overbar notation which is then omitted: $\bar{\bullet}$ is the non-dimensional version of $\bullet$: $\bar{x} = x/L_1$, $\bar{y} = y/L_1$, $\bar{t} = t/\tau$, $\bar{u} = u/L_1$, and $\bar{w} = w/L_1$ with the characteristic time $\tau = L/c_1$ where $c_1 = \sqrt{E_1/A_1}$. Derivatives of $u$ are found using the chains rule: $u_x = L_1\bar{u}_{\bar{x}}\bar{x}_x = \bar{u}_{\bar{x}}$ and $u_t = L_1\bar{u}_{\bar{t}}\bar{t}_t = L_1\bar{u}_{\bar{t}}/\tau$. Higher derivatives can be found in a similar way: $u_{xx} = \bar{u}_{\bar{x}\bar{x}}/L_1$ and $u_{tt} = L_1\bar{u}_{\bar{t}\bar{t}}/\tau^2$. For the second bar, when considered, we have $w_t = L_1\bar{w}_{\bar{t}}/\tau$ and $w_{tt} = L_1\bar{w}_{\bar{t}\bar{t}}/\tau^2$. Meanwhile $\bar{c} = c_2/c_1$ and bar length $\bar{l} = L_2/L_1$ for the second bar are introduced. Non-dimensional boundary conditions are also applied. For the Robin boundary condition, non-dimensional spring stiffness $\bar{k} = k/(EA)$ is introduced. Also, ratio $\alpha = E_1/E_2$ is introduced for two-bar contact condition. In the remainder, the upper bar notation is dropped for simplicity and all considered quantities are non-dimensional.