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

HAL
open science

# Nonsmooth modal analysis with boundary element method

Tianzheng Lu, Mathias Legrand

► To cite this version:

Tianzheng Lu, Mathias Legrand. Nonsmooth modal analysis with boundary element method. XI International Conference on Structural Dynamics, Nov 2020, Athènes, Greece. pp.205-212, 10.47964/1120.9016.19060 . hal-03136267

HAL Id: hal-03136267

https://hal.science/hal-03136267

Submitted on 10 Feb 2021

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
# NONSMOOTH MODAL ANALYSIS WITH BOUNDARY ELEMENT METHOD

Tianzheng Lu and Mathias Legrand

Structural Dynamics and Vibration Laboratory, McGill University
817 Sherbrooke St West, Montréal QC H3A 0C3 Canada

e-mail: tianzheng.lu@mail.mcgill.ca

mathias.legrand@mcgill.ca

**Keywords:** Nonsmooth Modal Analysis, Boundary Element Method, Shooting Method, Harmonic Balance Method, Structural Dynamics

**Abstract.** Numerical schemes based on the Boundary Element Method are proposed to perform Nonsmooth Modal Analysis. The latter aims at finding continuous families of periodic orbits of mechanical systems featuring unilateral contact constraints. In this contribution, a simple one-dimensional rod system is targeted. The frequency response, in the form of backbone-curve diagrams, and displacement field are presented. The proposed results compare well with existing studies on this topic.
# 1 INTRODUCTION

Within the framework of structural dynamics, linear modal analysis is a daily used tool in industry, aiming at predicting vibratory resonances of periodically forced mechanical systems by searching for continuous families of periodic solutions exhibited by the underlying autonomous (i.e. unforced) system. However, various challenges arise when possibly large-scale nonlinear dynamical systems are targeted and for which nonlinear modal analysis is needed instead [7]. Nonlinearity has commonly two forms: smooth or nonsmooth function of the state. The nonlinearity is said to be smooth when the governing equation involves a function of the state which is differentiable. Such functions are commonly polynomial. The nonlinearity is said to be nonsmooth otherwise, as for instance, unilateral contact conditions. Nonsmooth Modal Analysis (NSM) is a version of modal analysis dedicated to such systems. A few numerical schemes exist to perform NSM, with inherent limitations [11, 12]. In the present work, it is proposed to explore the numerical capabilities of two numerical schemes grounded on the Boundary Element Method (BEM) to perform NSM.

# 2 SYSTEM OF INTEREST

The system of interest, shown in Figure 1, is a one-dimensional bar of length L with a unilateral contact condition at its right tip. The unknown displacement field is $u(x, t)$, where

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

$x$ is the space coordinate and $t$ is time. Wave speed $c$ is space-independent and the classical one-dimensional wave equation governs the dynamics, that is

$$u_{,tt} - c^2 u_{,xx} = 0, \quad \forall x \in ]0; L[. \qquad (1)$$

Unilateral contact on the right tip is expressed as a Signorini boundary condition

$$g(t) \ge 0, \quad u_{,x}(L, t) \le 0, \quad g(L, t)u_{,x}(L, t) = 0. \qquad (2)$$

where $g(t) = g_0 - u(L, t)$. Also, since the bar is clamped on the left, $u(0, t) = 0$ applies.

In context of nonlinear modal analysis [7], families of periodic solution forming modal manifolds are computed through a numerical scheme which assumes that no impact law is required as the contact interface but instead, a switch between free-flight and sticking phases occurs.

# 3 METHODOLOGY

In this paper, two numerical schemes based on the Boundary Element Method are implemented. The BEM forms a family of methods for which the boundary of the domain of interest is the main ingredient of the formulation and full domain discretization is not required when initial condition and body force vanish [6]. It has the notable benefit of reducing the dimension of the formulation. Among the various incarnations of BEM, the Time Domain BEM (TD-BEM)
and Frequency Domain BEM (FD-BEM) are selected to perform nonsmooth modal analysis. The first strategy combines TD-BEM and a simplified shooting method, while the second one combines FD-BEM to Harmonic Balanced Method (HBM).

## 3.1 Time Domain Boundary Element Method and Shooting Method

TD-BEM is a version of BEM which requires the discretization of both time and space. By pre-multiplying the governing equation (1) with its fundamental solution [4] and performing integration in time and space and then integration by parts, the formulation can be transformed into Boundary Integral Equation (BIE), which here reads [3]

$$ u(x,t) = \int_{0}^{t} p^{*}(x, t, \xi, \tau)u(\xi, \tau) d\tau \Big|_{\xi=L} + \int_{0}^{t} p^{*}(x, t, \xi, \tau)u(\xi, \tau) d\tau \Big|_{\xi=0} \\ - \int_{0}^{t} u^{*}(x, t, \xi, \tau)p(\xi, \tau) d\tau \Big|_{\xi=L} - \int_{0}^{t} u^{*}(x, t, \xi, \tau)p(\xi, \tau) d\tau \Big|_{\xi=0} \\ - \frac{1}{c^2} \int_{0}^{L} u^{*}(x, t, \xi, 0)v(\xi, 0) d\xi + \frac{1}{c^2} \int_{0}^{L} \frac{\partial u^{*}(x, t, \xi, 0)}{\partial \tau} \Big|_{\tau=0} u(\xi, 0) d\xi \quad (3) $$

where $p$ represents traction on the boundary; $u^*(x, t, \xi, \tau)$ and $p^*(x, t, \xi, \tau)$ are the displacement and stress fundamental solutions. Any internal state for $x \in ]0; L[$ at time $t$ is dependent on the boundary displacement and traction at past times as well as initial states. In order to use TD-BEM, space and time are first discretized. Space integration of the initial states (last two terms in (3)) is performed by discretizing space into $n$ elements and $n+1$ nodes with $\Delta x = L/n$ such that $x_k = k\Delta x$ for $k = 0, \dots, n$. Time step is set to $\Delta t = \Delta x/c$ to guarantee stability of TD-BEM, such that $t^i = i\Delta t$ for $i = 0, \dots, n$. Equation (3) is discretized by reading it on the boundary, that is $x = 0$ and then $x = L$ for all $t^i$ for all $i = 0, \dots, n$; it now becomes

$$ \mathbf{H}^{ii} \mathbf{u}^i = \mathbf{G}^{ii} \mathbf{p}^i + \sum_{j=0}^{i-1} (\mathbf{G}^{ij} \mathbf{p}^j - \mathbf{H}^{ij} \mathbf{u}^j) + \mathbf{U}^i \mathbf{u}_0 + \mathbf{V}^i \mathbf{v}_0, \quad i = 1, \dots, n \quad (4) $$

where $\mathbf{u}^i \approx [u(0, t^i), u(L, t^i)]^\top$ and $\mathbf{p}^i \approx [p(0, t^i), p(L, t^i)]^\top$ are displacement and traction vectors on boundary, with size of $2 \times 1$. The quantities $\mathbf{H}^{ii}$, $\mathbf{G}^{ii}$, $\mathbf{H}^{ij}$, and $\mathbf{G}^{ij}$ are $2 \times 2$ coefficients matrices, evaluated from time-domain integration in (3) of fundamental solutions $u^*$ and $p^*$. Vectors $\mathbf{u}_0 = [u(x_k, 0)_{k=0, \dots, n}]$ and $\mathbf{v}_0 = [v(x_k, 0)_{k=0, \dots, n}]$ are the discretized initial displacement and velocity fields, with size $n+1$. The corresponding coefficient matrices of size $2 \times (n+1)$ are $\mathbf{U}^i$ and $\mathbf{V}^i$ stemming from space integration of fundamental solutions in (3). For systems without Signorini boundary conditions, Equation (4) can be solved at every time step starting at $i=1$ where half of the quantities in $\mathbf{u}^i$ and $\mathbf{p}^i$ are unknown, that is $(u(0, t^i), u(L, t^i))$ or $(u(0, t^i), p(L, t^i))$ or $(u(L, t^i), p(0, t^i))$ or $(p(0, t^i), p(L, t^i))$.

In order to account for unilateral contact conditions, TD-BEM is coupled with the floating boundary method [8]. It sees Signorini boundary condition as a switch between non-homogeneous Dirichlet and homogeneous Neumann boundary conditions. This technique does not require an impact law in the BEM format, thus avoiding the chattering or energy dissipation that might exist in other methods [1].

To perform NSM, a shooting method is used. It calculates, here in the discretized framework, the difference between the initial state $\mathbf{q}_0 = (\mathbf{u}_0, \mathbf{v}_0)$ and the corresponding final state $\mathbf{q}(\mathbf{q}_0, T)$ after a yet unknown period $T$ assumed to exist. For a periodic motion, the equality $\mathbf{q}(\mathbf{q}_0, T) = \mathbf{q}_0$
holds. In this work, $v_0 = 0$ is assumed [11], and the problem to be solved thus reduces to: Find vector $u_0$ and strictly positive integer $m$, where $T = m\Delta t$, such that

$$u_m = \hat{G}p - \hat{H}u + \hat{U}^m u_0 = u_0 \quad (5)$$

with the notations

$$\hat{G} = [\hat{G}^{m1} \ \hat{G}^{m2} \ \dots \ \hat{G}^{mm}], \quad p = \begin{bmatrix} p^1 \\ p^2 \\ \vdots \\ p^m \end{bmatrix}, \quad \hat{H} = [\hat{H}^{m1} \ \hat{H}^{m2} \ \dots \ \hat{H}^{mm}], \quad u = \begin{bmatrix} u^1 \\ u^2 \\ \vdots \\ u^m \end{bmatrix}$$

where $\hat{G}^{mj}$ and $\hat{H}^{mj}$ are $(n+1) \times 2$ coefficient matrices and $\hat{U}^m$ is a $(n+1) \times (n+1)$ coefficient matrix. Also, $u_m = [u(x_k, m\Delta t)_{k=0,...,n}]$. Since some $p^j$ and $u^j$ in (5) are unknown, Equation (4) is invoked, that is:

$$Hu = Gp + Uu_0 \quad (6)$$

with the notations

$$H = \begin{bmatrix} H^{11} & 0 & 0 & \dots & 0 \\ H^{21} & H^{22} & 0 & \dots & 0 \\ H^{31} & H^{32} & H^{33} & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ H^{m1} & H^{m2} & H^{m3} & \dots & H^{mm} \end{bmatrix}, \quad G = \begin{bmatrix} G^{11} & 0 & 0 & \dots & 0 \\ G^{21} & G^{22} & 0 & \dots & 0 \\ G^{31} & G^{32} & G^{33} & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ G^{m1} & G^{m2} & G^{m3} & \dots & G^{mm} \end{bmatrix}, \quad U = \begin{bmatrix} U^1 \\ U^2 \\ U^3 \\ \vdots \\ U^m \end{bmatrix}.$$

Combining linear systems (5) and (6) yields

$$\begin{bmatrix} H & -U \\ \hat{H} & I - \hat{U} \end{bmatrix} \begin{bmatrix} u \\ u^0 \end{bmatrix} = \begin{bmatrix} G \\ \hat{G} \end{bmatrix} p \quad (7)$$

where $I$ is a $(n+1) \times (n+1)$ identity matrix. In the floating boundary method, the Signorini condition is considered as a switch between Dirichlet and Neumann conditions. As soon as the contact duration $T_c$ is set, the time step at which the switch happens is known, so does the boundary condition (Dirichlet or Neumann) at each time step at $x = L$. Thus Equation (7) can be solved by reorganizing the known and unknown entries in $u$ and $p$, according to the boundary condition at $x = L$ at each time step. Accordingly, by first assuming only one contact switch per period, such a strategy systematically skim through given intervals of values on period $T$ and contact duration $T_c$ (for a chosen $\Delta t$). However, physically unacceptable solutions can be found. Accordingly, admissibility is systematically checked to ensure complementary conditions (2) are satisfied. Non-admissible solutions with non-admissible penetration or contact forces are discarded. Such way of Signorini boundary condition enforcement is similar in spirit with a precedent scheme [12], but also shows differences in detail.

## 3.2 Frequency Domain Boundary Element Method

Frequency Domain Boundary Element Method (FD-BEM) is a frequency-domain form of BEM which targets periodic solutions. Via a Fourier transform in time at frequency $\omega$

$$\tilde{u}(x) = \frac{1}{2\pi} \int_{-\infty}^{\infty} u(x, t) \exp(-i\omega t) dt \quad (8)$$
the wave equation (1) becomes the well-known one-dimensional Helmholtz equation

$$ \frac{d^2\tilde{u}(x)}{dx^2} + \kappa^2\tilde{u}(x) = 0 \qquad (9) $$

where $\kappa = \omega/c$ is the frequency number. Similar to TD-BEM, FD-BEM is based on a weighted residual form of (9), where the weight function is the fundamental solution to the Helmholtz equation [2] with the corresponding BIE

$$ 2\tilde{u}(x) = \tilde{u}(0) \cos \kappa x + \tilde{u}(L) \cos(\kappa(L-x)) - \frac{1}{\kappa}(\tilde{p}(0) \sin \kappa x - \tilde{p}(L) \sin(\kappa(L-x))) \quad (10) $$

where $\tilde{p}$ is the Fourier Transform of $p$. Reading (10) on the boundary leads to

$$ \begin{bmatrix} \kappa & 0 & -\kappa \cos \kappa L & \sin \kappa L \\ -\kappa \cos \kappa L & \sin \kappa L & \kappa & 0 \end{bmatrix} \begin{bmatrix} \tilde{u}(0) \\ \tilde{p}(0) \\ \tilde{u}(L) \\ \tilde{p}(L) \end{bmatrix} = 0, \quad \kappa \neq 0 \qquad (11) $$

Through (11), $\tilde{u}$ at frequency $\omega$ can be solved for a combination of Dirichlet and/or Neumann boundary conditions. The time domain solution is then recovered through an inverse Fourier Transform. However, such a formulation is not well adapted when $u(L, t)$ and $p(L, t)$ are both constrained by the complementary condition (2) expressed in time domain. This complementarity condition can be equivalently recast into the equality [9]

$$ r(t) \equiv p(L, t) + \max(0, \rho(u(L, t) - g_0) - p(L, t)) = 0, \quad \forall t \qquad (12) $$

where $\rho \in \mathbb{R}^+$. In the remainder, the residual $r(t)$ will be made to vanish in an integral sense only. To insert (12) into FD-BEM with the aim of finding periodic solutions, Fourier expansions are considered, that is

$$ p(L, t) = \frac{1}{2}a_0 + \sum_{n=1}^{\infty} a_n \cos n\omega_0 t + b_n \sin n\omega_0 t \qquad (13) $$

where $\omega_0 = 2\pi/T$ is the base frequency. The corresponding base frequency number is $\kappa_0 = \omega_0/c$. Knowing that $u(0, t) = 0$, it follows from (11) that

$$ u(L, t) = \frac{L}{2}a_0 + \sum_{n=1}^{\infty} \frac{\tan n\kappa_0 L}{n\kappa_0} (a_n \cos n\omega_0 t + b_n \sin n\omega_0 t) \qquad (14) $$

Since initial velocity is assumed to be zero [11], all sin terms vanishes. Only constant and cos terms remain in equation (13) and equation (14).

The truncated Fourier series of $u(L, t)$ and $p(L, t)$ with $m$ harmonics only is then substituted into (12). The Harmonic Balance Method (HBM) [5] is used to solve for the unknown coefficients ($a_0, a_1, a_2, \dots, a_m$). As a special case of Galerkin techniques, HBM reads

$$ \int_{0}^{T} \psi_{i}(t) r(t) dt = 0, \quad \forall i \qquad (15) $$

where the following family of functions is selected: $(1, \cos \omega_0 t, \cos 2\omega_0 t, \dots, \cos m\omega_0 t)$. With $m+1$ equations, $m+1$ unknown coefficients can be solved for by the (Non-Smooth) Newton method. The function $p(L, t)$ is then reconstructed through (13) and so are other needed functions used to build the Nonsmooth Modes of Vibration. Basic continuation is implemented [7] in order to construct the desired backbone curves, by increasing $\omega_0$ by a small value and solve again, with previous solution as initial guess.
# 4 RESULTS

In this section, the following is considered: $L = 1$, $c = 1$ and $g_0 = 0.5$. For TD-BEM, the bar is discretized into $n = 50$ elements in space with element length $\Delta x = L/n$ and time step is set accordingly. For FD-BEM, $m = 15$ is used.

The corresponding backbone curves are shown in figure 2. The backbone curve shows the

<figure><img src="image_4.png" /><figcaption>Figure 2: Backbone curve of the first mode. Red stars: TD-BEM/Shooting; blue circles: FD-BEM/HBM. Dashed line: linear mode</figcaption></figure>

frequency-dependence of the total energy, and compares well with existing solutions [11], except for backbone curve near linear mode. Since contact period length decreases as frequency $\omega$ decreases [11], short period of contact cannot be evaluated precisely with $m = 15$ harmonics. Due to inaccurate Fourier expansion and residual function, HBM struggles converging to useful periodic solutions. Thus this part of solution is skipped.

It can be observed that TD-BEM captures more than one possible solution for a given frequency of vibration, which shows similar pattern of nonsmooth modes distribution with existing results [10]. Two such solutions feature the same frequency but distinct contact durations, as shown in figure 3 (see points a and b). This is induced by the full internal resonance in the investigated system [12]. Although solution b shows a low-frequency pattern similar to solution a, extra energy is introduced due to higher frequency mode. Comparatively FD-BEM/HBM currently cannot guarantee convergence to internal resonance solutions. More advanced continuation technique might be required to find internal resonance branches.

Displacement solved by FD-BEM/HBM at point c on backbone curve is depicted in Figure 4. Overall, the displacement field is similar to solution at point a in Figure 3. It is worth noting that Signorini conditions in FD-BEM/HBM is enforced in integral sense through (15), unlike in TD-BEM the complementary condition is met pointwise, ie for every time step. Accordingly, residual pointwise penetration occurs as shown in figure 5.

In terms of computational efficiency, FD-BEM/HBM is much superior to TD-BEM/Shooting. With the suggested discretization settings, TD-BEM is about 700 times slower than FD-BEM/HBM to find the whole backbone curve. When $n$ increases, such difference is expected to increase, since the total number of combination between $T$ and $T_c$ to be skimmed increases in $O(n^2)$, while for FD-BEM/HBM the step size of continuation can be arbitrarily controlled.
<figure><img src="image_5.png" /><figcaption>Figure 3: TD-BEM admissible displacement fields $u(x, t)$ at $\omega = 1.848$ with different energy levels. Left: point a in Figure 2, $T_c = 0.6$; Right: point b in Figure 2, $T_c = 0.75$</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 4: Displacement field found by FD-BEM/HBM at point c in Figure 2</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 5: Zoom on displacement $u(L, t)$ found by FD-BEM/HBM at point c in Figure 2</figcaption></figure>

# 5 CONCLUSIONS

In this paper, two methods have been introduced to perform nonsmooth modal analysis on a simple academic system. Both methods have shown capability to converge to admissible solutions.

The TD-BEM/shooting methodology has shown good capability of capturing highly-detailed step-wise admissible solutions. However, its challenging implementation leads to high computational costs. Although it has shown no numerical energy dissipation in this one-dimensional case, this might be untrue for higher dimensional systems, in which case it would become ineligible for such modal investigations.

Comparatively, FD-BEM/HBM features a much lower computational cost, at the cost of solving Signorini conditions in an integral sense only. This formulation is also energy preserving by construction and could be extended to problems in two or three dimensions. More advanced continuation techniques, such as pseudo-arclength [7], are expected to be coupled into existing FD-BEM/HBM scheme.
References

[1] V. Acary. "Energy conservation and dissipation properties of time-integration methods for nonsmooth elastodynamics with contact". In: *Journal of Applied Mathematics and Mechanics / Zeitschrift für Angewandte Mathematik und Mechanik* 96.5 (May 2016), pp. 585–603.

[2] H. Antes. "A Short Course on Boundary Element Methods". Engineering school. Lecture. Technische Universität Braunschweig, Germany, Nov. 2010.

[3] J. A. M. Carrer and V. L. Costa. "Boundary element method formulations for the solution of the scalar wave equation in one-dimensional problems". In: *Journal of the Brazilian Society of Mechanical Sciences and Engineering* 37.3 (2015), pp. 959–971.

[4] K. F. Graff. *Wave motion in elastic solids*. Courier Corporation, 2012.

[5] J. C. Lindenlaub. "An approach for finding the sinusoidal steady state response of nonlinear systems". In: *Proc. 7th Ann. Allerton Conf. Circuit and System Theory*. Univ. Illinois. 1969.

[6] W. J. Mansur. "A time-stepping technique to solve wave propagation problems using the boundary element method". PhD thesis. University of Southampton, 1983.

[7] M. Peeters, R. Viguié, G. Sérandour, G. Kerschen, and J.-C. Golinval. "Nonlinear normal modes, Part II: Toward a practical computation using numerical continuation techniques". In: *Mechanical Systems and Signal Processing* 23.1 (2009), pp. 195–216.

[8] B. Shorr. *The Wave Finite Element Method*. Foundations of Engineering Mechanics. Springer, 2012.

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

[10] D. Urman and M. Legrand. "Nonlinear modal analysis of the bar in unilateral contact via analytical weak-solutions to the wave equation". In: *27th Canadian Congress of Applied Mechanics*. Sherbrooke, Canada, May 2019.

[11] J. Venkatesh, A. Thorin, and M. Legrand. "Nonlinear modal analysis of a one-dimensional bar undergoing unilateral contact via the time-domain boundary element method". In: *ASME 2017 International Design Engineering Technical Conferences*. Cleveland, United States, 2017.

[12] C. Yoong, A. Thorin, and M. Legrand. "Nonsmooth modal analysis of an elastic bar subject to a unilateral contact constraint". In: *Nonlinear Dynamics* (2018), pp. 1–24.