# Vibrations of an Axial Bar Experiencing Periodic Unilateral Contact Using the Wavelet Balance Method

Simon Jones, Mathias Legrand

► To cite this version:

Simon Jones, Mathias Legrand. Vibrations of an Axial Bar Experiencing Periodic Unilateral Contact Using the Wavelet Balance Method. 26th Conference on Mechanical Vibration and Noise, ASME, Aug 2014, Buffalo, United States. 10.1115/DETC2014-34725. [hal-01403614](https://hal.science/hal-01403614)

HAL Id: hal-01403614

https://hal.science/hal-01403614

Submitted on 2 Dec 2016

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# VIBRATIONS OF AN AXIAL BAR EXPERIENCING PERIODIC UNILATERAL CONTACT USING THE WAVELET BALANCE METHOD

## Simon Jones

Rose-Hulman Institute of Technology
Department of Mechanical Engineering
Terre Haute, Indiana, USA 47803
Email: jones5@rose-hulman.edu

## Mathias Legrand

McGill University
Department of Mechanical Engineering
Montreal, Quebec, Canada H3A 0C3
Email: mathias.legrand@mcgill.ca

## ABSTRACT

Efficiently predicting the vibratory responses of flexible structures which experience unilateral contact is becoming of high engineering importance. An example of such a system is a rotor blade within a turbine engine; small operating clearances and varying loading conditions often result in contact between the blade and the casing.

The method of weighted residuals is an effective approach to simulating such behaviour as it can efficiently enforce time-periodic solutions of lightly damped, flexible structures experiencing unilateral contact. The Harmonic Balance Method (HBM) based on Fourier expansion of the sought solution is a common formulation, though it is hypothesized wavelet bases that can sparsely define nonsmooth solutions may be superior. This is investigated herein using an axially vibrating rod with unilateral contact conditions. A distributional formulation in time is introduced allowing periodic, square-integrable trial functions to approximate the second-order equations. The mixed wavelet Petrov-Galerkin solutions are found to yield consistent or better results than HBM, with similar convergence rates and seemingly more accurate contact force prediction.

## INTRODUCTION

Efficiently predicting the vibratory responses of flexible structures which experience unilateral contact is becoming of high engineering importance. For example, consider aircraft engines where slender, twisted blades rotate at high rotational velocity within stationary casings where minimal clearance is desirable for turbine energy efficiency. Simulating contact between the blades and casing is not a trivial exercise since unilateral contact is usually described by inequalities and complementary conditions [1]. In the time domain, structural displacements and velocities which satisfy these contact conditions are known to respectively feature absolute continuity and bounded variation only [2]. This implies displacements are not necessarily differentiable everywhere in the domain and velocities may exhibit jumps; these types of problems are generally referred to as nonsmooth [3,4].

This class of unilateral problems can also be approached using periodic vibration theory. This allows the original initial-value formulations to be transformed into partial differential equations periodic in time [5]. In structural mechanics, two families of numerical techniques can efficiently predict periodic solutions in time. The first is the shooting method which consists of finding the initial conditions that result in periodic motion. The overall

approach relies on time integration of the governing equations over one period using a nonlinear solver which iterates on the initial conditions, and possibly the period. The second technique is based on weighted residual formulations, which are of interest in the current investigation. This method involves approximating the solution using a set of time-dependent basis functions, called trial functions, and enforcing the respective residual error to be orthogonal to an independent set of weighting functions [6, 7]. Unlike the shooting method which can become numerically sensitive to possible jumps in the velocity field, weighted residual techniques directly enforce the periodicity conditions while the remaining unilateral contact constraints and governing local equations of motion are satisfied in a weak integral sense [8]. It is worth noting that the well-known Harmonic Balance Method (HBM) is a weighted residual approach where both the trial and weighting functions are identical Fourier series. The main goal of the current work is to explore relevant basis functions whose order of smoothness can be adapted to a particular system to attain accurate approximations and rapid convergence.

The paper is broken into five sections: Section 2 introduces the system of interest; Section 3 reviews the weighted residual formulation in time; Section 4 provides an overview of the different basis functions under investigation; Section 5 specifies the unilateral contact model; and Section 6 presents and compares results for a number of basis function combinations.

## 1 System of interest

A schematic of a simple unilateral contact system showing a one-dimensional bar clamped on the left is depicted in Fig. 1; such a system is known to be a well-posed unilateral contact problem [9]. A gap $g$ exists between the tip and a rigid foundation. When the tip displacement is sufficiently large, the bar comes into contact with the rigid foundation and unilateral contact conditions must be satisfied. The existence of periodic solutions of period $T$ are assumed where $T$ is the period of the external forcing $f(x,t)$ acting on the bar. Accordingly, the unknown displacement has to satisfy the following complementary boundary value problem:

1. local equation of motion

$$ \rho S\ddot{u}(x,t) + E S u_{,xx}(x,t) = f(x,t), \forall x \in ]0, \ell[; \forall t \quad (1) $$
2. conditions of periodicity in time

$$
\begin{align}
u(x, t + T) &= u(x, t) \tag{2} \\
\dot{u}(x, t + T) &= \dot{u}(x, t)' \quad \forall x \in [0, \ell], \quad \forall t \notag
\end{align}
$$

3. boundary condition in displacement

$$
u(0, t) = 0, \quad \forall t \tag{3}
$$

4. unilateral contact conditions

$$
g - u(l, t) \ge 0, \lambda(t) \ge 0, \lambda(t)(g - u(l, t)) = 0, \forall t. \quad (4)
$$

Here $\rho$ signifies density, $S$ cross-sectional area, and $E$ elastic
modulus. The dot superscript represents a temporal derivative,
whereas the subscript $x$ represents a spatial derivative. Structural
damping is incorporated later by assuming $\beta$-damping of the
elastic modulus, as detailed in Section 2.

The periodicity conditions result in a problem which can be
formulated on a circle in time. Without loss of generality, the basis
functions are taken from an $L^2(\mathbb{S}^1)^N$ Hilbert space [10], where
the period of the steady-state solution $T$ has been normalized
to 1 [11].

<figure><img src="image_2.png" /><figcaption>FIGURE 1. ONE-DIMENSIONAL LINEAR BAR UNDERGOING UNILATERAL CONTACT CONDITIONS</figcaption></figure>

The quantity $\lambda(t)$ in Eq. (4) is a contact force, which stems
from the enforced non-penetration condition $g - u(\ell, t) \ge 0$ and
is necessarily positive (by convention). The complementarity con-
dition $\lambda(t) \cdot (g - u(\ell, t)) = 0$ states that the contact force $\lambda(t)$
and the distance $g - u(\ell, t)$ separating the rod's end-tip from the rigid
foundation may not be zero at the same time. These three con-
ditions are such that the mathematical object pairing the contact
force to the displacement is not a function in the usual sense. It is
also known that the displacement field may be non-differentiable
at given instants within the period of motion and the velocity field
may be discontinuous. This motivates the derivation of numerical
techniques capable of efficiently handling nonsmoothness. As
a first approach, the unilateral contact inequalities from Eq. (4)
are simplified and replaced by a penalty function; an exponential
spring $K_p$ is used to approximate the contact forces. The gap $g$
in Fig. 1 separates the tip of the rod and the exponential spring.
The penalty function is of the form

$$
f_c(u(\ell, t)) = \max(a_c(e^{\alpha(u(\ell, t)-g_0)} - 1), 0). \quad (5)
$$

Specific values for the penalty function are provided with the
other model parameters in Section 4.

2 Weighted residual formulations in time

The method of weighted residuals is a classic method of
obtaining numerical solutions to boundary value problems by
expanding the sought solution as a finite sequence of basis func-
tions (i.e. trial functions) in a proper functional space [6]. The
subsequent residual is rendered orthogonal to a set a linearly
independent functions of the same space (i.e. weightnig func-
tions) through an inner product. Trial and weighting functions
may contain the same basis function but not necessarily: The
Galerkin method is a special case where the weighting functions
are taken from the same function basis as that of the trial functions;
the Petrov-Galerkin method involves the selection of weighting
function comprising a basis which is independent of the trial
function [12]. The functional space into which the solution is
sought is such that the boundary conditions in space and time are
satisfied. In this study, the displacement of the investigated rod is
numerically sought in various Hilbert spaces of square integrable
functions defined on the circle $\mathbb{S}^1$.

To solve Eqs. (1), (2), (3), and (4), the unknown displacement
is expanded into a truncated series of *N* functions separated in
space and time:

$$
u(x, t) = \sum_{i=1}^{N} \vartheta(x) u_i(t) \quad (6)
$$

The standard Finite Element Method is implemented for the
spatial variable using two-node linear rod elements [13, 14].
This yields the following vector ordinary differential equation
of size *N*:

$$
M\ddot{u}(t) + C\dot{u}(t) + Ku(t) = f_{\text{ext}}(u(t), t) \quad (7)
$$

together with the remaining periodicity conditions in time and the
unilateral contact condition. It should be noted that Eq. (7) is an
approximate equation of a vibrating structure which may be of
one, two or three dimensions. The system is constrained so that
condition (3) is now satisfied. Here **M** and **K** are the standard mass
and stiffness matrices for a rod, and **β**-damping is enforced such
that **C** = **β****K** to account for structural damping. The displacement
vector **u**(t) in Eq. (7) stores the temporal unknowns uᵢ(t), i =
1, ..., N, where N is the number of spatial degrees-of-freedom.
Similarly, **f**<sub>ext</sub>(**u**(t), t) stores the external forcing functions for
each degree-of-freedom as well as the contact forces stemming
from Eq. (5).

The remainder of the derivation involves a weighted residual
formulation in the time dimension. Three forms of the weighted
residual method are discussed: the strong integral form, the weak
form, and the distributional formulation. For this discussion it is
understood that **u** and **v** are time-dependent vectors, *i.e.* **u**(t) and
**v**(t), but for the sake of clarity the *t* is omitted below.

2.1 Strong integral form

The standard weighted residual formulation of a differential
equation is commonly termed the strong form. Taking the inner
product of Eq. (7) with a weighting function stored column-wise
in **v** results in the strong integral form of the equation: find **u** ∈
$H^2(\mathbb{S}^1)^N$ such that $\forall \mathbf{v} \in L^2(\mathbb{S}^1)^N$:

$$ \int_{\mathbb{S}^1} (\mathbf{v}^\mathrm{T} \mathbf{M} \ddot{\mathbf{u}} + \mathbf{v}^\mathrm{T} \mathbf{C} \dot{\mathbf{u}} + \mathbf{v}^\mathrm{T} \mathbf{K} \mathbf{u} - \mathbf{v}^\mathrm{T} \mathbf{f}_{\mathrm{ext}}(\mathbf{u}, t)) \, dt = 0 \quad (8) $$

and the superscript T denotes a transpose. This strong form of the equation is not necessarily the best framework for obtaining a solution [14]; for this example involving a vector ordinary differential equation of order 2, the solution must be at least $H^2$, limiting the permissible basis of trial functions.

## 2.2 Weak form

The respective weak form of the weighted residual statement can be obtained by performing one integration by parts over the domain $\mathbb{S}^1$ for all terms containing a double time derivative in Eq. (8). This results in: find $\mathbf{u} \in H^1(\mathbb{S}^1)^N$ such that $\forall \mathbf{v} \in H^1(\mathbb{S}^1)^N$,

$$ \int_{\mathbb{S}^1} (-\dot{\mathbf{v}}^\mathrm{T} \mathbf{M} \ddot{\mathbf{u}} + \mathbf{v}^\mathrm{T} \mathbf{C} \dot{\mathbf{u}} + \mathbf{v}^\mathrm{T} \mathbf{K} \mathbf{u} - \mathbf{v}^\mathrm{T} \mathbf{f}_{\mathrm{ext}}(\mathbf{u}, t)) \, dt = 0 \quad (9) $$

The integral form of the weak formulation offers the advantage of shifting a portion of the functional smoothness requirement from the trial functions onto the weighting functions. More precisely, both the trial and weighting functions must now be $H^1$. This allows the trial functions to be chosen from a wider permissible space [14].

## 2.3 Formulation in a distributional sense

The above procedure of obtaining the solution from the weak formulation by performing integration by parts on the weighted residual statement can be extended one step further. Theoretically, this extension is not necessary since it is known that displacements should be absolutely continuous in time, though it may assist in the numerical derivations and allow very simple basis functions to be considered [15]. A weaker formulation is proposed by integrating again the terms involving time derivatives of the trial functions. This formulation can be understood in the sense of distributions, also known as generalized functions, ie: find $\mathbf{u} \in L^2(\mathbb{S}^1)^N$ such that $\forall \mathbf{v} \in H^2(\mathbb{S}^1)^N$:

$$ \int_{\mathbb{S}^1} (\ddot{\mathbf{v}}^\mathrm{T} \mathbf{M} \mathbf{u} - \dot{\mathbf{v}}^\mathrm{T} \mathbf{C} \mathbf{u} + \mathbf{v}^\mathrm{T} \mathbf{K} \mathbf{u} - \mathbf{v}^\mathrm{T} \mathbf{f}_{\mathrm{ext}}(\mathbf{u}, t)) \, dt = 0 \quad (10) $$

Here the double time differential on the field variable is transferred to the weighting function and the continuity requirement on the trial function is reduced. As discussed later, the desired displacement functions $\mathbf{u}$ can now be described using a series of constant piecewise functions for instance.

## 2.4 Time discretization

Let each nodal displacement $u_i(t)$ be approximated by a combination of $M$ linearly independent trial functions $q_k(t)$, $k = 1, \dots, M$ defined on $t \in [0, 1)$. Similarly, let each $v_i(t)$ be approximated by a combination of linearly independent weighting

functions $p_k(t)$, $k = 1, \dots, M$ also defined on $t \in \mathbb{S}^1$. Collectively, this leads to

$$ u_i(t) = \sum_{k=1}^{M} a_{ik} q_k(t), \quad i = 1, \dots, N \quad (11) $$

and

$$ v_i(t) = \sum_{k=1}^{M} b_{ik} p_k(t), \quad i = 1, \dots, N \quad (12) $$

or, equivalently in a vector form:

$$ \mathbf{u}(t) = \mathbf{\Theta}(t)\mathbf{a} \quad \text{and} \quad \mathbf{v}(t) = \mathbf{\Gamma}(t)\mathbf{b} \quad (13) $$

In Eq. (13), $\mathbf{\Theta}(t)$ and $\mathbf{\Gamma}(t)$ are rectangular matrices of dimension $N \times NM$; $\mathbf{a}$ and $\mathbf{b}$ are vectors of size $NM \times 1$.

Each of the strong, weak, and distributional formulations can now be discretized with proper basis functions. The overall goal of this paper is to discuss the numerical properties of certain wavelet bases to efficiently solve the problem of interest. Depending on the selected formulation, the order of derivation acting on $\mathbf{u}(t)$ and $\mathbf{v}(t)$ (i.e. $\mathbf{\Theta}(t)$ and $\mathbf{\Gamma}(t)$ respectively) will affect the admissible basis functions. The corresponding discretized versions are:

* Strong integral form:

$$ \int_{\mathbb{S}^1} ((\mathbf{\Gamma}^\mathrm{T} \mathbf{M} \ddot{\mathbf{\Theta}} + \mathbf{\Gamma}^\mathrm{T} \mathbf{C} \dot{\mathbf{\Theta}} + \mathbf{\Gamma}^\mathrm{T} \mathbf{K} \mathbf{\Theta})\mathbf{a} - \mathbf{\Gamma}^\mathrm{T} \mathbf{f}_{\mathrm{ext}}(\mathbf{\Theta}\mathbf{a}, t)) \, dt = 0 \quad (14) $$

* Weak form:

$$ \int_{\mathbb{S}^1} ((-\dot{\mathbf{\Gamma}}^\mathrm{T} \mathbf{M} \dot{\mathbf{\Theta}} + \mathbf{\Gamma}^\mathrm{T} \mathbf{C} \dot{\mathbf{\Theta}} + \mathbf{\Gamma}^\mathrm{T} \mathbf{K} \mathbf{\Theta})\mathbf{a} - \mathbf{\Gamma}^\mathrm{T} \mathbf{f}_{\mathrm{ext}}(\mathbf{\Theta}\mathbf{a}, t)) \, dt = 0 \quad (15) $$

* Distributional integral form:

$$ \int_{\mathbb{S}^1} ((\dot{\mathbf{\Gamma}}^\mathrm{T} \mathbf{M} \mathbf{\Theta} - \dot{\mathbf{\Gamma}}^\mathrm{T} \mathbf{C} \mathbf{\Theta} + \mathbf{\Gamma}^\mathrm{T} \mathbf{K} \mathbf{\Theta})\mathbf{a} - \mathbf{\Gamma}^\mathrm{T} \mathbf{f}_{\mathrm{ext}}(\mathbf{\Theta}\mathbf{a}, t)) \, dt = 0 \quad (16) $$

The resulting nonlinear equations can generically be recast in the following form

$$ \mathbf{G}\mathbf{a} - \mathbf{f}(\mathbf{a}) = \mathbf{0} \quad (17) $$

where $\mathbf{G}\mathbf{a}$ and $\mathbf{f}(\mathbf{a})$ respectively stand as the linear internal and the nonlinear external contributions in Eq. (14), (15), and (16).

## 3 Trial and weighting function bases

The selection of functional bases to be used in the above approaches is an important factor in approximation accuracy and computational efficiency [13, 14]. A priori selection of the optimal bases for unilateral contact problems is not a trivial matter due to potential nonsmoothness of the response. In the current investigation a number of bases are investigated to compare the quality of approximation, including Fourier functions, B-spline wavelets, Daubechies wavelets, and Haar wavelets.
## 3.1 Harmonic balance and Fourier series

The harmonic balance method (HBM) is a special case of the weighted residual method where the Fourier basis is used for both the trial and weighting functions [16]. This technique is particularly effective when dealing with smooth nonlinear systems; convergence is often reached with very few terms. It has been used to study steady state response of turbine engine blades with friction dampers using a multiterm approximation [17]. This approach was extended to unilateral contact and friction conditions [18] through an Alternating Frequency/Time domain strategy proposed by Cameron et al. [19] and Pilipchuk [20]. It is worthy to note that the HBM formulations of Eqs. (14), (15), and (16) are identical.

The Fourier basis for $L^2(\mathbb{S}^1)$ is defined as $\{1\} \cup \{\cos(2m\pi t), \sin(2m\pi t) | m \in \mathbb{N}^*\}$ where $m$ signifies the harmonic number of the function. The first six functions are shown in Fig. 2. Fourier basis functions feature an infinite degree of

<figure><img src="image_3.png" /><figcaption>FIGURE 2. FIRST SIX FUNCTIONS OF THE FOURIER BASIS</figcaption></figure>

smoothness. While this property can be beneficial in some cases, it is unclear whether Fourier functions are optimal when simulating potentially nonsmooth problems. It is possible that a large number of harmonics would be required to accurately capture the nonsmooth response, or the approximation may exhibit Gibbs phenomenon at localized discontinuities in the sought displacement and velocity fields. Accordingly, other basis functions featuring a lesser degree of functional smoothness are explored.

## 3.2 Brief review of discrete orthogonal wavelets

Discrete orthogonal wavelet families are composed of highly localized, oscillatory functions which provide a basis of $L^2(\mathbb{R})$ and can be adapted to the periodic domain $L^2(\mathbb{S}^1)$ [21]. These localized characteristics, or compact support, allow sparse representation of piecewise signals including transients and singularities [22].

This makes them useful functions for use in the Galerkin approach when nonsmooth solutions are predicted [23]. There are a large number of wavelet families and definitions thereof; both Mallat [21] and Strang [23] provide excellent introductions to wavelet theory and history. Galerkin methods using appropriate discrete wavelet families as the trial functions have been shown to accurately approximate the solutions to both ordinary and partial differential equations [24, 25, 26, 27, 28, 29].

The discrete wavelet family is built from scaling functions $\phi(t)$ and wavelet functions $\psi(t)$. These functions are analogous to the low-pass and high-pass filters of a filter bank [23]; decomposition using scaling functions will give a “smoothed” approximation of the original signal, while decomposition using wavelet functions provides the details of the signal, or high-frequency content.

The exact decomposition of a continuous time signal $y(t)$ can be written

$$y(t) = \sum_k g_k \phi_{J,k}(t) + \sum_{j=J}^{m} \sum_{\ell} h_{j,\ell} \psi_{j,\ell}(t) \quad (18)$$

and

$$\phi_{J,k}(t) = 2^{J/2} \phi(2^J t - k) \quad (19)$$

$$\psi_{J,k}(t) = 2^{J/2} \psi(2^J t - k) \quad (20)$$

where $J, k \in \mathbb{Z}$; $J$ is the dilation parameter (i.e. level), $k$ is the translation parameter, and $m$ is the maximum resolution given by the sampling rate of the function $y(t)$.

The span of the scaling functions at level $J$ is commonly denoted $V_J$, while the wavelet span is denoted $W_J$. For orthogonal wavelet families, $W_J$ is the orthogonal complement to $V_J$ in $V_{J+1}$

$$V_{J+1} = V_J \oplus W_J. \quad (21)$$

Provided the wavelet family is orthogonal [23], the space of square-integrable functions on the real line $L^2(\mathbb{R})$ can be decomposed using multiresolution analysis as a nested sequence of closed subspaces [21]

$$\dots \subset V_{-2} \subset V_{-1} \subset V_0 \subset V_1 \subset V_2 \subset \dots \subset L^2(\mathbb{R}) \quad (22)$$

such that

$$\lim_{J \to \infty} V_J = L^2(\mathbb{R}). \quad (23)$$

This implies that $V_J$, the subspace composed of the set of scaling functions at level $J$, can approximate any function in $L^2(\mathbb{R})$. A reduced orthonormal basis of $L^2(\mathbb{R})$ is constructed by truncating the wavelet terms, resulting in

$$u(t) = \sum_k u_k \phi_{J,k}(t) \quad (24)$$
The accuracy of this approximation increases as the level *J* is increased. This property of the scaling functions makes them excellent trial and weighting functions in weighted residual methods because they can be adapted to the accuracy level required. The reduced orthonormal basis of scaling functions is used in the current investigation.

**3.2.1 Periodization of wavelet families** Standard wavelet definitions (i.e. scaling and wavelet functions) are commonly built on the real line. The functions can be adapted to periodic functions of $L^2(\mathbb{S}^1)$ by utilizing a standard periodization technique [21, 30, 31]. Let $\phi^{(p)}(t)$ be the periodized form of the scaling function $\phi(t)$ defined on $\mathbb{R}$

$$ \phi_{J,k}^{(p)}(t) = \sum_{\ell \in \mathbb{Z}} \phi_{J,k}(t-\ell) \quad 0 \le k \le 2^J - 1. \quad (25) $$

This is equivalent to “wrapping around” the support $\mathbb{R}$ on $\mathbb{S}^1$ through summation. The finite size of the interval results in the condition $J \ge 0$

$$ \mathbf{V}_0 \subset \mathbf{V}_1 \subset \mathbf{V}_2 \subset \mathbf{V}_3 \subset \cdots \subset L^2(\mathbb{S}^1) \quad (26) $$

such that

$$ \lim_{J \to \infty} \mathbf{V}_J = L^2(\mathbb{S}^1). \quad (27) $$

A number of periodic discrete wavelet families exist [21]. The investigation considers three families to determine how they perform in unilateral, nonsmooth contact problems: B-spline, Daubechies, and Haar.

## 3.3 Orthogonal cubic B-spline scaling functions

The scaling function $\phi$ for the orthogonal cubic B-spline wavelet family built on $\mathbb{R}$ is given as [32]

$$ \phi(t) = \sum_k c_k B_3(t-k) \quad (28) $$

where cubic B-spline $B_3(t)$ can be written using the following formulas

$$ B_3(t) = \begin{cases} \frac{t^3}{6} & (0 \le t \le 1) \\ -\frac{t^3}{2} + 2t^2 - 2t + \frac{2}{3} & (1 \le t \le 2) \\ \frac{t^3}{2} - 4t^2 + 10t - \frac{22}{3} & (2 \le t \le 3) \\ -\frac{t^3}{6} + 2t^2 - 8t + \frac{32}{3} & (3 \le t \le 4) \\ 0 & \text{(otherwise)}. \end{cases} \quad (29) $$

The coefficients $c_k$ can be determined using [32]

$$ c_k = c_{-k} = \frac{1}{\pi} \int_0^\pi \frac{\cos(k\xi)}{\sqrt{p_3(\cos\xi)}} d\xi \quad (k \ge 0) \quad (30) $$

where *k* is an integer. The cubic polynomial $p_3$ is given as

$$ p_3(y) = \frac{1}{630}(y^3 + 60y^2 + 297y + 272). \quad (31) $$

Numerical simulations showed that truncating the summation in Eq. (28) at $-50 \le k \le 50$ is sufficient; larger *k* terms add negligibly to the summation. After periodization to $\mathbb{S}^1$ and normalization by $\int_0^1 \phi_{J,k}(t)dt = 1$ [33], sample orthogonal cubic B-spline scaling functions for $J = 0, 1, \dots, 5$ are shown in Fig. 3

<figure><img src="image_4.png" /><figcaption>FIGURE 3. EXAMPLES OF THE FIRST SIX LEVELS OF PERIODIZED ORTHOGONAL CUBIC B-SPLINE SCALING FUNCTIONS</figcaption></figure>

## 3.4 Daubechies Wavelets

The Daubechies wavelet family is defined by a set of *L* filter coefficients {$p_\ell : \ell = 0, 1, \dots, L-1$}, where *L* is an even integer. The scaling function is defined by the fundamental two-scale equation [33, 34]

$$ \phi(t) = \sum_{\ell=0}^{L-1} p_{\ell} \phi(2t - \ell) \quad (32) $$

which has fundamental support over the finite intervals $[0, L-1]$. This equation can be used to determine the value of the scaling function at dyadic points $t = n/2^J$, $n = 0, 1, \dots$ using the algorithm provided by Chen et al. [33].

The corresponding scaling functions are highly nonsmooth and fractal in nature: as one increases the resolution of the functions, the shape does not converge but rather continues to increase
in complexity. This makes accurately estimating the inner products of such scaling functions with each other prone to error when numerical integration is used [35].

<figure><img src="image_5.png" /><figcaption>FIGURE 4. EXAMPLES OF THE FIRST SIX LEVELS OF PERIODIZED 6 COEFFICIENT DAUBECHIES SCALING FUNCTIONS (DB6)</figcaption></figure>

When Daubechies scaling functions are used in a Galerkin approach, it is necessary to derive the inner products of the scaling function with itself or derivatives of itself. The exact solution to these inner products can be found by using the recursive nature of the fundamental equation on $L^2(\mathbb{R})$ [24, 33]; the solution to these inner products are commonly referred to as connection coefficients. When Daubechies scaling functions are periodized on $\mathbb{S}^1$, the wrapping procedure results in functions which are no longer scale-invariant at low $J$ values (scale-invariance requires wavelets at any scale to be a pure dilation of the mother-wavelet). Fig. 4 provides examples of the periodized scaling functions for $L = 6$. Interestingly, the lack of scale-invariance for small $J$ values does not invalidate the connection coefficient algorithms derived for unbounded domains; the connection coefficients can simply be wrapped around the periodic domain as necessary [35].

## 3.5 Haar scaling functions

The simplest Daubechies wavelet family requires only two filter coefficients ($p_0 = p_1 = 1$) and is commonly known as the Haar wavelet family [36]. The Haar scaling functions are rectangular tophat-type functions; the father scaling function is

defined on $t \in \mathbb{R}$ as

$$ \phi(t) = \begin{cases} 1 & (0 \le t < 1) \\ 0 & (\text{otherwise}). \end{cases} \qquad (33) $$

Since the compact support of the father scaling function is $\mathbb{S}^1$, the periodized function is equivalent. Example Haar scaling functions for $J = 0, 1, \dots, 5$ are shown in Fig. 5.

<figure><img src="image_6.png" /><figcaption>FIGURE 5. EXAMPLES OF THE FIRST SIX LEVELS OF PERIODIZED HAAR SCALING FUNCTIONS</figcaption></figure>

## 4 Model development

The mixed finite-element/wavelet-Galerkin model of the axial bar undergoing unilateral contact introduced in Section 1 has the following properties: elastic modulus $E = 70$ GPa; cross-sectional area $S = 25$ cm$^2$; density $\rho = 2700$ kg/m$^3$; total length $L = 1$ m; external force amplitude $F = 250kN$; and damping factor $\beta$. A total of 25 finite elements are used. The fundamental frequency of this rod is approximately 1273 Hz.

The values of $a_c$ and $\alpha$ in Eq. (5) were selected to simulate a reasonably rigid penalty function while maintaining good conditioning of the system of nonlinear equations. Table 1 lists the two forcing frequency cases considered: 150 Hz and 1275 Hz. The maximum tip displacements without contact for this system are approximately 1.4 mm at 150 Hz and 29 mm at 1275 Hz.

To act as a comparison solution, the unilateral contact finite-element equations detailed above are solved using a variable-order numerical differentiation formula (NDF) time-stepping algorithm [37]. As the time-stepping solution does not solve for the periodic response directly, the solution is deemed to have
<table><caption>**TABLE 1.** Test cases</caption><thead><tr><th></th><th>Case 1</th><th>Case 2</th></tr></thead><tbody><tr><td>forcing frequency</td><td>150 Hz</td><td>1275 Hz</td></tr><tr><td>penalty constant, a<sub>c</sub></td><td>1 × 10<sup>8</sup></td><td>1 × 10<sup>8</sup></td></tr><tr><td>penalty exponent, α</td><td>50</td><td>20</td></tr><tr><td>damping factor, β</td><td>1 × 10<sup>-4</sup></td><td>5 × 10<sup>-6</sup></td></tr><tr><td>initial gap, g<sub>0</sub></td><td>1.0 mm</td><td>14.5 mm</td></tr></tbody></table>

converged to its periodic state when the rms of the relative error between the tip displacements for the *i* and *i* − 1 periods is below 10⁻⁷. This is possible due to the structural damping terms listed in Table 1 causing any transient response to decay relatively quickly.

## 4.1 Weighted residual formulation

Table 2 lists the combination of trial and weighting functions employed in the current investigation; where wavelet families are listed (B-spline, Haar, Daubechies), it is implied the scaling functions are used according to Eq. (24). The integrals involved

<table><caption>**TABLE 2.** List of trial and weighting function combinations</caption><thead><tr><th>Trial function</th><th>Weighting function</th></tr></thead><tbody><tr><td>Fourier</td><td>Fourier</td></tr><tr><td>B-spline</td><td>B-spline</td></tr><tr><td>Daubechies 6</td><td>Daubechies 6</td></tr><tr><td>Haar</td><td>Fourier</td></tr><tr><td>Daubechies 6</td><td>Fourier</td></tr><tr><td>Haar</td><td>B-spline</td></tr></tbody></table>

in Eq. (17) are computed as discrete inner products using 256 points for all but the Galerkin formulations using Daubechies functions. For the Daubechies Galerkin cases, the inner products are derived using the connection coefficients as mentioned in Section 3.4.

# 5 Results

The six combinations of trial and weighting functions listed in Table 2 are employed to solve the nonlinear weighted residual formulation for forcing frequencies of 150 Hz and 1275 Hz.

## 5.1 Tip displacement

Samples of the approximate tip displacement responses at 150 Hz at 1275 Hz using 64 (i.e *J* = 6) basis functions are provided in Fig. 6 and Fig. 7, respectively.

For this number of basis functions HBM (Fourier:Fourier) approximates the tip displacements well at both 150 Hz and 1275 Hz compared to the time-stepping solution. However, close examination of the contact plateau for the 150 Hz case shows oscillations due to Gibbs phenomenon and rounding off of the sharp gradient changes. This will be further discussed in Section 5.3.

<figure><img src="image_7.png" /><figcaption>**FIGURE 6.** TIP DISPLACEMENT VERSUS NORMALIZED TIME FOR CASE 1</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**FIGURE 7.** TIP DISPLACEMENT VERSUS NORMALIZED TIME FOR CASE 2</figcaption></figure>

Alternatively, the Haar:Fourier combination appears to approximate the tip displacement outside the contact zone less well due to the blocky nature of the Haar scaling functions. However, the contact plateau is well represented by the constant Haar scaling functions. The Haar basis appears to be less suited for the response at 1275 Hz as the contact plateau is smaller and not as flat; the tip penetration through the contact boundary is a result of the inexact modeling of hard contact using the exponential spring penalty function.

## 5.2 Tip velocity

Consideration is also given to the accuracy with which the trial functions considered can determine velocities. For all displacement approximations where the trial basis functions can be differentiated pointwise to approximate velocities using Eq. 11,
the following equality is used:

$$ \dot{u}_i(t) = \sum_{k=1}^{M} a_{ik} \dot{q}_k(t). \qquad (34) $$

The only cases where Eq. (34) is not valid are those involving Haar trial functions; Haar scaling functions are piecewise constant functions thus they cannot be directly differentiated in the usual sense (i.e. pointwise) to determine velocities. Velocities for these cases can be approximated by assuming the derivative of a Haar scaling function is the combination of a positive delta function at $t_1$ and a negative delta function at $t_2$, where $t_1$ corresponds to the positive jump discontinuity in the Haar scaling function and $t_2$ corresponds to the negative jump. Using a set of these functions as $\dot{q}_k(t)$ in Eq. (34) provides a reasonable approximation to the corresponding velocity field for cases involving Haar functions.

The approximate tip velocity response at 150 Hz is plotted in Fig. 8 for all cases using 64 basis functions. It is visible from these approximations that there is a sharp jump in velocity due to the contact condition. As expected, this results in ringing due to Gibbs phenomenon when the Fourier trial functions are used. Interestingly, the ringing is more pronounced in the Galerkin B-spline case. The DB6 trial functions appear to do the best job of approximating the tip response. This is attributed to the compact support of the scaling functions allowing accurate representation of rapid changes in gradient. For cases involving Haar trial functions the velocity function envelope is reasonably approximated by the delta function representation discussed above.

<figure><img src="image_9.png" /><figcaption>FIGURE 8. TIP VELOCITY VERSUS NORMALIZED TIME FOR CASE 1</figcaption></figure>

## 5.3 Tip contact force

The tip contact force is calculated using the penalty function provided in Eq. (5) in conjunction with the predicted tip displacement and presented in Fig. 9; the contact force has been normalized by the magnitude of the external force input. Notice the time duration of the contact forces coincide with the tip displacement plateaus of Fig. 6 as this is the only period during which the tip is in contact with the rigid foundation.

The effect of Gibbs phenomenon can be seen in the HBM and Galerkin B-spline cases. The cases where Haar trial functions are used approximate the contact force relatively well. It is hypothesized that when a rigid contact law is enforced the Haar scaling functions will perform even better relative to the other functions; this will be investigated in future work. Again, the DB6 trial functions appear to most accurately approximate the contact forces compared to the other cases.

## 5.4 Energy norm convergence

Convergence of the normalized energy norms for each basis combination at 150 Hz and 1275 Hz are given in Fig. 10(a) and 10(b), respectively. For these results the distributional form of the governing equations are used in all cases.

The system energy is calculated as

$$ E = \frac{1}{2} \dot{\mathbf{u}}^{\mathrm{T}} \mathbf{M} \dot{\mathbf{u}} + \frac{1}{2} \mathbf{u}^{\mathrm{T}} \mathbf{K} \mathbf{u} + E_{\mathrm{contact}}. \qquad (35) $$

The scalar energy norm is calculated as the rms value of the system energy $E$ over a single period; this is then normalized by

the energy norm from the time-stepping solution. The velocity terms $\dot{\mathbf{u}}$ are calculated as detailed in Section 5.2. As shown in Fig. 10, all combinations of trial and weighting functions converge to values similar to the time-stepping solution (i.e. energy norms are approximately unity using 128 basis functions). At 150 Hz, the HBM and Galerkin B-spline approaches show energy norms close to unity even using relatively few basis functions. The Petrov-Galerkin combinations utilizing Haar trial functions over-predict the system energy when few trial functions are used. This is evident in Fig. 6(e) where the blocky nature gives a coarse representation of the periodic response, thus affecting the resulting energy norm. An interesting comparison is between the DB6:DB6 and the DB6:Fourier curves. When Fourier functions are used as the weighting functions, the energy norm is better approximated. This implies projecting the residuals onto the Fourier basis, which has been shown to represent the solution relatively well, increases the accuracy of the prediction.

At 1275 Hz, the six combinations converge to a common energy norm more slowly than at the lower forcing frequency. In this case the HBM and Galerkin B-spline approaches again converge quickly to the time-stepping solution. The cases with Haar trial functions show larger difference with respect to the normalizing solution than that for the 150 Hz case. This is attributed to the displacement response at this frequency (see Fig. 7(c)); there is no distinct contact plateau thus the Haar trial functions do not approximate the solution as closely as the 150 Hz case. Finally, it is less obvious now that DB6:Fourier combination out performs the DB6:DB6 method at fewer basis functions.
<figure><img src="image_10.png" /></figure>

**FIGURE 9.** NORMALIZED CONTACT FORCE VERSUS NORMALIZED TIME FOR CASE 1

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

**FIGURE 10.** SYSTEM ENERGY NORM CONVERGENCE USING SIX TRIAL:WEIGHTING FUNCTION COMBINATIONS

## 5.5 Tip displacement relative error

It is also desirable to quantify the approximation error for the tip displacements predicted using the weighted residual combinations where the time-stepping solution is again used as a comparison solution. It provides a good relative measure to compare the accuracy of trial:weighting function combinations. The rms error over one period is plotted for the six basis function combinations in Fig. 11(a) for 150 Hz and Fig. 11(b) for 1275 Hz. As shown in Fig. 11(a), the HBM and Galerkin B-spline approaches have the lowest error when only 8 basis functions are used. This supports the convergence behaviour noted above for these combinations. As the number of basis functions is increased, the Galerkin B-spline case does not perform as well as the HBM method. This is attributed to a localized error which was found to develop at the end of the contact plateau (see Fig. 6(d)). The Galerkin DB6 case performs poorly using 16 trial functions, but rapidly improves to relative error values comparable to the HBM method. The ability of the nonsmooth DB6 scaling functions to accurately predict both the smooth portions of the response as well as the sharp gradient variation is attributed to both the compact support and vanishing moment properties of the Galerkin B-spline case. The relative error using Haar trial functions is highest out of the combinations investigated. This is expected due to the blocky nature of the Haar scaling functions; however the mid-points of the Haar functions match very well to the trend of the time-stepping solution. This suggests that if only a low resolution approximation of the response is required, the Haar scaling functions still perform well.

The relative error at 1275 Hz shows similar trends to those discussed at 150 Hz. One important difference is the Galerkin

B-spline error as the number of basis functions increases: as the unknowns are increased from 32 to 128, the relative error stays reasonably constant. This suggests that this combination of basis functions has reached its limit of prediction accuracy.

## 5.6 Frequency response curve

Another feature of interest is the ability to construct the frequency response function for the unilateral contact condition. The penalty function simulating contact is effectively a stiffening condition, thus the resonant peaks for the system are expected to increase in frequency relative to the no-contact response. To capture this nonlinear response an arc-length continuation method [38] is utilized. The frequency response function over the system's fundamental frequency is shown in Fig. 12 using both the Fourier:Fourier and Haar:Fourier combinations with 64 basis functions.

As shown, the stiffening behaviour of the unilateral contact condition is captured by both trial:weighting function combinations; other combinations showed similar approximations. These results are promising as it implies the arc-length continuation method can be accurately applied to both Galerkin and Petrov-Galerkin approaches using the distributional formulation for a number of different bases.

# 6 Conclusions

The method of weighted residuals for capturing the periodic responses of unilateral contact problems has been investigated for
<figure><img src="image_12.png" /><figcaption>(a) at 150 Hz</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(b) at 1275 Hz</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>FIGURE 12. FREQUENCY RESPONSE FUNCTION OVER THE SECOND FUNDAMENTAL FREQUENCY OF THE UNILATERAL CONTACT SYSTEM</figcaption></figure>

a number of trial:weighting function combinations. The contact condition is simulated using an exponential penalty function approach. To extend the allowable trial function bases, a weak and a distributional formulations are presented which transfer trial function continuity requirements to the weighting functions. This allowed piecewise constant Haar scaling functions to be used as a trial basis in the current investigation. Results show that a number of trial:weighting function combinations produce accurate solutions which rapidly converge as the size of the discrete spaces is increased. As expected, Fourier functions perform well as a trial

basis, though nonsmooth functions such as Haar and Daubechies scaling functions are also attractive since they provide comparable prediction accuracy and even out perform the Fourier functions in some measures. It is also shown that all the basis combinations considered can be used in an arc-length continuation framework to capture the nonlinear frequency response of the unilateral contact problem.

Research is continuing on this subject. The goal is to improve on the penalty approach by implementing an exact unilateral contact method (e.g Lagrange multiplier, LCP, etc.). This should result in perfectly flat contact plateaus within the periodic response. We also intend to apply this method to a standard rotor blade to validate the method for three-dimensional models involving more complicated contact interface geometries.

## REFERENCES

[1] Wriggers, P., 2006. *Computational Contact Mechanics*. Springer.

[2] Schindler, T., and Acary, V., 2012. "Timestepping schemes for nonsmooth dynamics based on discontinuous galerkin methods: Definition and outlook". *Mathematics and Computers in Simulation*.

[3] Acary, V., and Brogliato, B., 2008. *Numerical methods for nonsmooth dynamical systems*. Springer.

[4] Leine, R., and Nijmeijer, H., 2006. *Dynamics and bifurcations of non-smooth mechanical systems*. Springer.

[5] Laxalde, D., and Legrand, M., 2011. "Nonlinear modal analysis of mechanical systems with frictionless contact interfaces". *Computational Mechanics*, **47**(4), pp. 469–478.

[6] Lapidus, L., and Seinfeld, J. H., 1971. *Numerical Solution of Ordinary Differential Equations*. Academic Press Inc., U.S., NY.

[7] Reddien, G. W., 1980. "Projection methods for two-point boundary value problems". *SIAM Review*, **22**(2), pp. 156–171.

[8] Wang, Y., 1995. "Dynamics of unsymmetric piecewise-linear/non-linear systems using finite elements in time". *Journal of Sound and Vibration*, **185**(1), pp. 155–170.

[9] Schatzman, M., and Bercovier, M., 1989. "Numerical approximation of a wave equation with unilateral constraints". *Mathematics of Computation*, **53**(187), pp. 55–79.

[10] Holschneider, M., 1990. "Wavelet analysis on the circle". *Journal of Mathematical Physics*, **31**(1), pp. 39–44.

[11] Berberian, S. K., 1999. *Introduction to Hilbert space*. American Mathematical Society.

[12] Durran, D. R., 1999. *Numerical methods for wave equations in geophysical fluid dynamics*. Springer Verlag, New York.

[13] Bathe, K.-J., 1982. *Finite element procedures in engineering analysis*. Prentice-Hall.

[14] Liu, G. R., and Quek, S. S., 2003. *Finite element method: a practical course*. Oxford: Butterworth-Heinemann.

[15] Studer, C., 2009. *Numerics of unilateral contacts and friction*. Springer.

[16] Mickens, R. E., 1984. "Comments on the method of harmonic balance". *Journal of Sound and Vibration*, **94**(3), pp. 456–460.

[17] Wang, J. H., and Chen, W. K., 1993. "Investigation of the vibration of a blade with friction damper by HBM". *Journal of Sound and Vibration*.
of Engineering for Gas Turbines and Power, 115(2), pp. 294–299.

[18] Legrand, M., Roques, S., Pierre, C., Cartraud, P., and Pe- seux, B., 2006. “n-dimensional Harmonic Balance Method extended to non-explicit nonlinearities”. *European Journal of Computational Mechanics*, **15**(1-2-3), pp. 269–280.

[19] Cameron, T. M., and Griffin, J. H., 1989. “An alternating frequency/time domain method for calculating the steady-state response of nonlinear dynamic systems”. *Journal of Applied Mechanics*, **56**(1), pp. 149–154.

[20] Pilipchuk, V. N., 1999. “Application of special nonsmooth temporal transformations to linear and nonlinear systems under discontinuous and impulsive excitation”. *Nonlinear Dynamics*, **18**(3), pp. 203–234.

[21] Mallat, S., 2009. *A wavelet tour of signal processing*. Academic Press, Burlington MA.

[22] Burrus, C. S., Gopinath, R. A., and Guo, H., 1998. *Introduction to wavelets and wavelet transforms: a primer*. Prentice Hall, Upper Saddle River, NJ.

[23] Strang, G., and Nguyen, T., 1996. *Wavelets and Filter Banks*. Wellesley-Cambridge Press.

[24] Beylkin, G., 1992. “On the representation of operators in bases of compactly supported wavelets”. *SIAM Journal of Numerical Analysis*, **6**(6), pp. 1716–1740.

[25] Beylkin, G., and Keiser, J., 1996. On the adaptive numerical solution of nonlinear partial differential equations in wavelet bases. Tech. rep., University of Colorado.

[26] Pernot, S., and Lamarque, C.-H., 2001. “A wavelet-Galerkin procedure to investigate time-periodic systems: transient vibration and stability analysis”. *Journal of Sound and Vibration*, **245**, pp. 845–875.

[27] Williams, J., and Amaratunga, K., 1994. “Introduction to

wavelets in engineering”. *International Journal for Numerical Methods in Engineering*, **37**, pp. 2365–2388.

[28] Christon, M. A., and Roach, D. W., 2000. “The numerical performance of wavelets for PDEs: the multi-scale finite element”. *Computational Mechanics*, **25**(2), pp. 230–244.

[29] Burgos, R., e Silva, R., and Santos, M., 2010. “Direct solution of differential equations using the wavelet-Galerkin method”. *Mecanica Computacional, XXIX*, pp. 4573–4584.

[30] Erlebacher, G., Hussanini, M. Y., and Jameson, L. M., 1996. *Wavelets theory and application*. Oxford University Press.

[31] Meyer, Y., 1990. *Ondelettes*. Herman, Paris.

[32] Blatter, C., 1998. *Wavelets a primer*. A K Peters, Natick, MA.

[33] Chen, M., Hwang, C., and Shih, Y., 1996. “The computation of wavelet-Galerkin approximation on a bounded interval”. *International Journal for Numerical Methods in Engineering*, **39**, pp. 2921–2944.

[34] Daubechies, I., 1992. *Ten lectures on wavelets*. SIAM, Philadelphia, PA.

[35] Restrepo, J., and Leaf, G., 1993. *Inner product computations using periodized Daubechies wavelets*. University of California, Los Angeles, Technical Report.

[36] Haar, A., 1910. “Zur theorie der orthogonalen funktionensysteme erste mitteilung”. *Mathematishe Annalen*, **69**(3), pp. 331–371.

[37] Shampine, L. F., and Reichelt, M. W., 1997. “The matlab ode suite”. *SIAM Journal on Scientific Computing*, **18**, pp. 1–22.

[38] von Groll, G., and Ewins, D. J., 2001. “The harmonic balance method with arc-length continuation in rotor/stator contact problems”. *International Journal for Numerical Methods in Engineering*, **241**(2), pp. 223–233.