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

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 a effective approach to simulating such behaviour as it can efficiently enforce timeperiodic 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 S R u.x; t/CESu;xx.x; t/ D f .x; t/;8x 20; `Œ;8t (1) 1 2. conditions of periodicity in time u.x; t C T / D u.x; t/ P u.x; t C T / D P u.x; t/ ; 8x 2 Œ0; `; 8t (2) 3. boundary condition in displacement


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.

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: 3. boundary condition in displacement 4. unilateral contact conditions g u.`; t / 0; .t / 0; .t /.g u.`; t // D 0; 8t: Here 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ˇ-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 .S 1 / N Hilbert space [10], where the period of the steady-state solution T has been normalized to 1 [11]. The quantity .t / in Eq. (4) is a contact force, which stems from the enforced non-penetration condition g u.`; t / 0 and is necessarily positive (by convention). The complementarity condition .t / g u.`; t / D 0 states that the contact force .t / and the distance g u.`; t / separating the rod's end-tip from the rigid foundation may not be zero at the same time. These three conditions 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.`; t // D max a c .e˛. u.`;t / g 0 / 1/; 0 : Specific values for the penalty function are provided with the other model parameters in Section 4.

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 functions (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 functions) 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 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: 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 : 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 DˇK to account for structural damping. The displacement vector u.t / in Eq. (7) stores the temporal unknowns u i .t /, i D 1; : : : ; N , where N is the number of spatial degrees-of-freedom. Similarly, f ext .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.

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

Weak form
The respective weak form of the weighted residual statement can be obtained by performing one integration by parts over the domain S 1 for all terms containing a double time derivative in Eq. (8). This results in: find u 2 H 1 .S 1 / N such that 8v 2 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].

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 u 2 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 u can now be described using a series of constant piecewise functions for instance.

Time discretization
Let each nodal displacement u i .t / be approximated by a combination of M linearly independent trial functions q k .t /, k D 1; : : : ; M defined on t 2 OE0; 1/. Similarly, let each v i .t / be approximated by a combination of linearly independent weighting functions p k .t /, k D 1; : : : ; M also defined on t 2 S 1 . Collectively, this leads to and or, equivalently in a vector form: In Eq. (13), ‚.t / and .t / are rectangular matrices of dimension N NM ; a and b are vectors of size NM 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 u.t / and v.t / (i.e. ‚.t / and .t / respectively) will affect the admissible basis functions. The corresponding discretized versions are: Weak form: Distributional integral form: The resulting nonlinear equations can generically be recast in the following form where Ga and f.a/ respectively stand as the linear internal and the nonlinear external contributions in Eq. (14), (15), and (16).

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.

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 .S 1 / is defined as f1g [ fcos.2m t /; sin.2m t /jm 2 N g 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 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.

Brief review of discrete orthogonal wavelets
Discrete orthogonal wavelet families are composed of highly localized, oscillatory functions which provide a basis of L 2 .R/ and can be adapted to the periodic domain L 2 .S 1 / [21]. These localized characteristics, or compact support, allow sparse repre-sentation 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 .t / and wavelet functions .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 and where J; k 2 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 Provided the wavelet family is orthogonal [23], the space of square-integrable functions on the real line L 2 .R/ can be decomposed using multiresolution analysis as a nested sequence of closed subspaces [21] This implies that V J , the subspace composed of the set of scaling functions at level J , can approximate any function in L 2 .R/. A reduced orthonormal basis of L 2 .R/ is constructed by truncating the wavelet terms, resulting in 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.

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 .S 1 / by utilizing a standard periodization technique [21,30,31]. Let .p/ .t / be the periodized form of the scaling function .t / defined on R This is equivalent to "wrapping around" the support R on S 1 through summation. The finite size of the interval results in the condition J 0 such that 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.

Orthogonal cubic B-spline scaling functions
The scaling function for the orthogonal cubic B-spline wavelet family built on R is given as [32] .
where cubic B-spline B 3 .t / can be written using the following formulas .otherwise/: The coefficients c k can be determined using [32] where k is an integer. The cubic polynomial p 3 is given as Numerical simulations showed that truncating the summation in Eq. (28) at 50 Ä k Ä 50 is sufficient; larger k terms add negligibly to the summation. After periodization to S 1 and normalization by R 1 0 J;k .t /dt D 1 [33], sample orthogonal cubic B-spline scaling functions for J D 0; 1; : : : ; 5 are shown in Fig. 3

Daubechies Wavelets
The Daubechies wavelet family is defined by a set of L filter coefficients fp`W`D 0; 1; : : : ; L 1g, where L is an even integer. The scaling function is defined by the fundamental twoscale equation [33,34] .
which has fundamental support over the finite intervals OE0; L 1. This equation can be used to determine the value of the scaling function at dyadic points t D n=2 J , n D 0; 1; : : : 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].  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 .R/ [24,33]; the solution to these inner products are commonly referred to as connection coefficients. When Daubechies scaling functions are periodized on 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 D 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].

Haar scaling functions
The simplest Daubechies wavelet family requires only two filter coefficients (p 0 D p 1 D 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 2 R as Since the compact support of the father scaling function is S 1 , the periodized function is equivalent. Example Haar scaling functions for J D 0; 1; : : : ; 5 are shown in Fig. 5. The values of a c and˛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 finiteelement equations detailed above are solved using a variableorder 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  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 7 . This is possible due to the structural damping terms listed in Table 1 causing any transient response to decay relatively quickly.  (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.

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.

Tip displacement
Samples of the approximate tip displacement responses at 150 Hz at 1275 Hz using 64 (i.e J D 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.  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.

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: 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 P 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 Bspline 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.

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.

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 The scalar energy norm is calculated as the rms value of the system energy E over a single period; this is then normalized by . 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 overpredict 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.

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 Bspline 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 the Daubechies functions. 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.

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.

Conclusions
The method of weighted residuals for capturing the periodic responses of unilateral contact problems has been investigated for   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.