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

HAL
open science

# A Higher Order Linear Least Square Fit for the Assessment of Integral and Non-Integral Vibrations With Blade Tip Timing

Felix Figaschewsky, Arnold Kühhorn, Thomas Giersch

► To cite this version:

Felix Figaschewsky, Arnold Kühhorn, Thomas Giersch. A Higher Order Linear Least Square Fit for the Assessment of Integral and Non-Integral Vibrations With Blade Tip Timing. 17th International Symposium on Transport Phenomena and Dynamics of Rotating Machinery (ISROMAC2017), Dec 2017, Maui, United States. hal-02980982

HAL Id: hal-02980982

https://hal.science/hal-02980982

Submitted on 27 Oct 2020

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
# A Higher Order Linear Least Square Fit for the Assessment of Integral and Non-Integral Vibrations With Blade Tip Timing

Felix Figaschewsky<fnref n="1" />*, Arnold Kühhorn<fnref n="1" />, Thomas Giersch<fnref n="2" />

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

ISROMAC 2017

International Symposium on Transport Phenomena and Dynamics of Rotating Machinery

Maui, Hawaii

December 16-21, 2017

## Abstract

This paper aims at improving the robustness and accuracy of the least square fit technique utilized in blade tip timing (BTT) measurements of blade vibrations by proposing two modifications. The first proposal is to replace the lines of the original least square problem by differences of consecutive lines. Thereby, the static deflection as well as the circumferential blade positioning error cancels out and the robustness is improved by removing these uncertainties inherently.

The second proposal is to replace the fit of piecewise constant vibration amplitudes within the chosen block length by a linear or cubic spline in the frequency (integral) or time (non-integral) regime. This does not only suppress overshoots due to distorted acceleration or deceleration manoeuvres but also allows for a "coarser analysis grid" (i.e. larger block length) without loosing amplitude accuracy. Thereby it smooths out random errors more efficiently and increases the orthogonality of the relevant EO or vibration frequency to unwanted signal components.

## Keywords

Blade Tip Timing — Non-Intrusive Measurements — Blade Vibration

## INTRODUCTION

A main challenge in the development of rotating turbomachinery is the safe and reliable design of blades with respect to high cycle fatigue (HCF). For many applications resonant operation cannot be avoided entirely and the apparent vibration levels need to be assessed with measurements. In this context the measurement of rotor blade vibrations is a difficult task. Traditionally, some of the blades are instrumented with strain gauges. Measuring blade vibrations very accurately, strain gauges have numerous drawbacks, since they influence aerodynamic as well as structural properties, need a complex telemetry system, have a very limited lifetime and nevertheless only the instrumented blades can be monitored, which may not be sufficient to capture the maximum vibration level of a mistuned bladed disk assembly.

Thus, recent work focuses on the enhancement of non-intrusive measurement techniques such as Blade Tip Timing (BTT). Meanwhile, BTT is an established technology allowing for the assessment of all blades in an assembly. However, this assessment is restricted to those modes showing a remarkable tip deflection and comes along with difficulties in the post-processing of the heavily undersampled BTT data. There are several algorithms reported in the literature. These can be classified by the presence of an "once per revolution" (OPR) sensor yielding a time signal of the circumferential position of the shaft and by the application of a specific structural model fitted to the data sets[1, 2, 3]. The scope of this investigation

is limited to OPR based BTT measurements without using specific structural models. The absence of specific structural models (e.g. single degree of freedom oscillator) allows a free interpretation of the observed frequency response functions (FRFs) as those of a mistuned multi-degree of freedom system. Additionally, either the exciting engine order (EO) or the vibration frequency is assumed to be known or determined by other available algorithms[4, 5, 6, 7].

A popular method to determine the necessary vibration amplitudes from the undersampled BTT data is the linear least square fit[8, 9]. However, BTT inevitably observes the superposition of probe and blade positioning errors, static deflection and blade vibration, see Fig. 1. Hence, a correction for positioning errors and static deflection needs to be modelled[10] causing additional uncertainty in the measurement. Also, a suitable analysis block length (nr. of revolutions for fit) needs to be chosen. This choice is always a trade-off between accuracy and robustness. One the one hand, it needs to be sufficiently long to minimize random errors and achieve a certain orthogonality to irrelevant signals contents. On the other hand, it has to be sufficiently short to resolve the maximum vibration amplitudes accurately.

This paper aims at improving the robustness and accuracy of the least square fit technique by proposing two modifications. The first proposal is to replace the lines of the original least square problem by differences of consecutive lines. Thereby, the static deflection as well as the circumferential blade po-

<footnote n="1">¹Chair of Structural Mechanics and Vehicle Vibration Technology, Brandenburg University of Technology Cottbus-Senftenberg, Germany</footnote>

<footnote n="2">²former staff of ¹</footnote>

<footnote>*Corresponding author: Felix.Figaschewsky@b-tu.de</footnote>
sitioning error cancels out and the robustness is improved by removing these uncertainties inherently. This is similar to other available algorithms[2, 3], however this approach differs by remaining linear and not assuming a specific structural model.

The second proposal is to replace the fit of piecewise constant vibration amplitudes within the chosen block length by a linear or cubic spline in the frequency (integral) or time (non-integral) regime. The splines are constructed based on finite elements with basic functions of appropriate order ensuring continuity of the vibration amplitude as functions of frequency or time as well as of the derivative for the cubic spline method. This does not only suppress overshoots due to distorted acceleration or deceleration manoeuvres but also allows for a "coarser analysis grid" (i.e. larger block length) without losing amplitude accuracy. Thereby it smooths out random errors more efficiently and increases the orthogonality of the relevant EO or vibration frequency to unwanted signal components.

This is especially of interest if the vibration response is a superposition of integral responses with different EOs and/or non-integral responses of multiple frequencies. It also improves the assessment of responses of higher modes with low signal to noise ratios.

Before the finite element based formulation of the least square fit is introduced, the paper recaps some fundamental BTT relationships to visualize the improvements by the proposed modifications. A subset of nominal system modes (SNM) representation[11] of a compressor rotor is used to generate BTT data. This allows to study the performance of the proposed algorithm with certain level of noise and mistuning as well as superposed vibrations of different mode families.

# 1. METHODS

## 1.1 Structural Model

In order to show the performance of the proposed changes it is suitable to analyse sufficiently simple data with known signal contents. However, the generated signals may also contain as many features as possible BTT measurements face in real engine applications. Thus, the generated data shall account for:

* the vibration response of a mistuned rotor;

* changes in the static deflection due to the acceleration of the rotor;

* blade individual positioning errors;

* probe individual positioning errors;

* noise in the determination of time of arrivals.

In order to cover the demands above a subset of nominal system modes (SNM) representation[11, 12] of a research compressor rotor is chosen. The scope of this paper is limited to linear vibration analysis. Hence, a fully non-linear analysis step of a finite element (FE) model including temperature and centrifugal loads is followed by a numerical modal analysis assuming sufficiently small perturbations.

This procedure is carried out for different operating conditions and yields static deflection as well as eigenfrequencies and mode shapes necessary for the reduction of the large FE model matrices. Assuming mass normalised mode shapes and Rayleigh damping, the reduced linear ordinary differential equations of motions for the amplitudes q of the nominal system modeshapes read:

$$ \ddot{q} + [\operatorname{diag}\{2\delta_k\} + G]\dot{q} + [\operatorname{diag}\{\omega_{0,k}^2\} + \Delta K]\dot{q} = f(t) \quad (1) $$

The fully populated ΔK matrix generates coupling between the tuned system modes and is referred to the mistuning matrix. It is a scaled superposition of reduced nominal sector stiffness matrices Ksec with the matrix of nominal sector modeshapes Ψ. The scaling is applied sectorwise with the scaling coefficients cᵢ proportional to the blade individual mode frequency:

$$ \Delta K = \sum_{i=1}^{N} c_i \Psi^T K_{\text{sec}} \Psi \quad (2) $$

The G matrix is the reduced gyroscopic matrix and accounts for Coriolis effects leading to a frequency split of the forward and backward travelling wave modes - even for the tuned case with vanishing mistuning matrix ΔK. It has been already shown, that such a SNM based ROM is able to model the structural dynamics of mistuned vibrations of real rotors appropriately[12]. The involved matrices are in general a function of the rotational speed Ω and can be calculated for each speed the FE model has been run with. Afterwards they can be interpolated in between the discrete steps. However, these studies restrict to constant system matrices since the changes during a single resonance crossing are assumed to be sufficiently small, whereas the evolution of the mean static deflection with rotational speed is taken into account. The right hand side in 1 is able to reflect both engine order like (integral) and arbitrary (non-integral) excitation. The reduced equations of motion 1 are solved with an implicit Newmark time integration scheme yielding the instantaneous tip displacement u(t) as a superposition of the mean static deformation ū(Ω) and the dynamic displacement u'(t):

$$ u(t) = \bar{u}(\Omega(t)) + u'(t) = \bar{u}(\Omega(t)) + \Psi q(t) \quad (3) $$

The circumferential blade deflection Φ(t) observed by the BTT system (see fig. 1) is given by the circumferential position of the blade at the given axial BTT probe position and is in general a non-linear function of the modal amplitudes q(t). Presuming sufficiently small amplitudes it can be usually described by a linearisation and the application of mode individual calibration factors. Nevertheless, the utilized procedure retains the full non-linearity. The circumferential motion of the undeformed sector is calculated by adding a blade individual offset to the circumferential position of the OPR sensor φ(t) which is a function of the shaft speed Ω(t) and time:

$$ \varphi(t) = \varphi_0 + \int_{0}^{t} \Omega(\tau)d\tau \quad (4) $$
<figure><img src="image_4.png" /><figcaption>Figure 1. Sketch of the i-th Blade Tip Arriving at j-th BTT Probe at time $t_{i,j}$. Left: $R - \varphi$ plane. Right: $x - \varphi$ plane.</figcaption></figure>

Based on the instantaneous tip deflection $u(t)$ relative to the reference position and the circumferential position of the undeformed sector it is possible to decide whether any of the blades has triggered a BTT probe during the evolution from the preceding to the current time step. Thereby analytical BTT data can be generated that cover the initially formulated demands with respect to mistuning of the rotor, static deflections, positioning errors of blades and probes as well as noise being superposed on the analytically determined time of arrival (TOA).

## 1.2 Fundamental BTT Relations

Before the higher-order least-square fit is introduced, some fundamental BTT relations are briefly summarized. This paper focusses on measurements with an OPR sensor that detects the position of the target on the shaft $\varphi(t)$ once per revolution. Assuming a smoothly and sufficiently slow varying rotational speed $\Omega$ the time signal $\varphi(t)$ can be evaluated for any given time by a suitable interpolation technique and is therefore regarded to be known. The blade root position of the i-th blade is determined by applying a corresponding shift $\Delta\varphi_i$ to the shaft position $\varphi(t)$, see fig. 1. According to fig. 1 the blade tip position of the i-th blade is referred to $\Phi_i(t)$ and is calculated by adding its static $\bar{\Phi}_i$ and dynamic deflection $\Phi'_i$ to the blade root position:

$$ \Phi_i(t) = \varphi(t) + \Delta\varphi_i + \bar{\Phi}_i(\Omega(t)) + \Phi'_i(t) \quad (5) $$

The blade individual circumferential shift $\Delta\varphi_i$ shall not only account for the nominal blade shift $2\pi i/N$ but also for any blade individual positioning errors in axial and/or circumferential direction as well as in stagger angle. The fact that these errors might evolve differently on each of the blades when applying centrifugal loads is captured by the blade individual static deflection $\Phi_i(\Omega(t))$ as a function of $\Omega(t)$. The circumferential position $\Theta_{i,j}$ of the j-th BTT probe in the casing is calculated by:

$$ \Theta_{i,j} = \Theta_j^0 + \Delta\Theta_{i,j} \quad (6) $$

The index i takes into account, that each of the blades might see the j-th probe at a different position. This is due to the

fact that the deviation $\Delta\Phi_{i,j}$ from the nominal position $\Theta_j^0$ is not only a function of the axial $\Delta x_j$ and circumferential $\Delta\vartheta_j$ probe positioning error but also of the combination of axial probe positioning error and blade individual stagger angle $\alpha + \Delta\alpha_i$, see fig. 1 (right):

$$ \Delta\Theta_{i,j} = \Delta\vartheta_j + \Delta x_j \tan(\alpha + \Delta\alpha_i) \quad (7) $$

The second term becomes most relevant for blades with high nominal stagger angle $\alpha$ close to $90^\circ$, such that $\alpha + \beta = 90^\circ$ with small $\beta$. This allows the following estimation of the term:

$$ \begin{aligned} \Delta x_j \tan(\alpha + \Delta\alpha_i) &\approx \Delta x_j \tan(\alpha) + \frac{\Delta x_j \Delta\alpha_i}{\cos^2(\alpha)} \\ &= \frac{\Delta x_j}{\cos(\alpha)} \left( \sin(\alpha) + \frac{\Delta\alpha_i}{\cos(\alpha)} \right) \\ &\approx \frac{\Delta x_j}{\cos(\alpha)} \left( 1 + \frac{\Delta\alpha_i}{\beta} \right) \\ \Rightarrow \Delta\Theta_{i,j} &= \Delta\vartheta_j + \frac{\Delta x_j}{\cos(\alpha)} + \frac{\Delta x_j \Delta\alpha_i}{\cos(\alpha)\beta} \end{aligned} \quad (8) $$

Hence, assuming sufficiently small deviations $\Delta\alpha_i$ of the blade's stagger angle with respect to the angle $\beta$, the last term is small compared with the second one and may be neglected. This simplifies the determination of the circumferential position $\Theta_j$ of the j-th BTT probe:

$$ \Theta_j = \Theta_j^0 + \Delta\Theta_j = \Theta_j^0 + \Delta\vartheta_j + \Delta x_j \tan(\alpha) \quad (9) $$

Figure 1 shows the situation when the j-th BTT probe fires the time of arrival $t_{i,j}$ due to the triggering by blade i. At this moment in time the circumferential position of the i-th blade tip $\Phi_i$ coincides with the j-th probe position $\Theta_j$ leading to the fundamental equation of BTT measurements:

$$ \Theta_j = \Phi_i(t_{i,j}) = \varphi(t_{i,j}) + \Delta\varphi_i + \bar{\Phi}_i(\Omega(t_{i,j})) + \Phi'_i(t_{i,j}) \quad (10) $$

Equation 10 shows that BTT is only able to measure the superposition of blade and probe positioning errors, static
deflection and vibration. In order to determine the vibration level at a given point in time, at least an approximation of the positioning errors as well as an approximation of the static deflection needs to be known:

$$ \Phi'_{i}(t_{i,j}) = \underbrace{\Theta_{j}^{0}}_{\text{known}} - \underbrace{\varphi(t_{i,j})}_{\text{measured}} + \underbrace{\Delta\Theta_{j} - \Delta\varphi_{i} - \bar{\Phi}_{i}(\Omega(t_{i,j}))}_{\text{approximated}} \quad (11) $$

In order to derive the desired approximations measurements with small (ideally vanishing) vibration amplitudes $\Phi'_{i}(t_{i,j}) \approx 0$ are utilized:

$$ \Delta\Theta_{j} - \Delta\varphi_{i} - \bar{\Phi}_{i}(\Omega(t_{i,j})) \approx \Theta_{j}^{0} - \varphi(t_{i,j}) \quad (12) $$

However, the static deflection and the included evolution of the blade individual positioning error with rotational speed are a non-linear function of $\Omega$, i.e. there must be at least a speed range below and above the studied resonance crossing with small vibrations in order to derive a linear approximation of this function. This prerequisite can often not be satisfied in real engine applications possibly leading to significant errors in the determined vibration amplitudes. This especially holds if integral excitations are analysed. Since this excitation is fixed in a stationary frame of reference, the vibration amplitude for a given rotational speed and circumferential position, i.e. probe position, does not change with the number of revolutions. Hence, at a single probe the BTT system is not able to distinguish between static deflection $\bar{\Phi}_{i}(\Omega(t_{i,j}))$ and vibration $\Phi'_{i}(t_{i,j})$.

To overcome this significant drawback a first modification is proposed. Instead of studying equation 11, we study the difference of consecutive hits of the i-th blade at probes j and j+1:

$$ \Phi'_{i}(t_{i,j+1}) - \Phi'_{i}(t_{i,j}) = \underbrace{\Theta_{j+1}^{0} - \Theta_{j}^{0}}_{\text{known}} - \underbrace{(\varphi(t_{i,j+1}) - \varphi(t_{i,j}))}_{\text{measured}} \\ + \Delta\Theta_{j+1} - \Delta\Theta_{j} - \underbrace{(\Delta\varphi_{i} - \Delta\varphi_{i})}_{\text{approximated}} \\ - \underbrace{(\bar{\Phi}_{i}(\Omega(t_{i,j+1})) - \bar{\Phi}_{i}(\Omega(t_{i,j})))}_{\text{approximated}} \quad (13) $$

Since the static deflection is a function of the rotational speed which does not vary significantly within a fraction of one rotor revolution, the respective difference as well as the influence of the blade individual positioning error $\Delta\varphi_{i}$ cancels out:

$$ \Phi'_{i}(t_{i,j+1}) - \Phi'_{i}(t_{i,j}) = \underbrace{\Theta_{j+1}^{0} - \Theta_{j}^{0}}_{\text{known}} - \underbrace{(\varphi(t_{i,j+1}) - \varphi(t_{i,j}))}_{\text{measured}} \\ + \Delta\Theta_{j+1} - \Delta\Theta_{j} \quad (14) $$

Hence, by regarding differences of consecutive TOAs two major sources of uncertainties are removed inherently, which increases the reliability and accuracy of the measurement significantly. The remaining difference of the probe positioning errors $\Delta\Theta_{j+1} - \Delta\Theta_{j}$ is only a weak function of the rotational speed anymore, compare eq. 9 and the argumentation above,

such that a single speed range with small vibrations is sufficient to determine the necessary information. This speed range can be at very low speeds or wherever no vibration is expected which is a remarkable improvement compared with the initial prerequisites above.

However, one drawback is introduced. Whereas the formulation 11 allows an interpretation of the raw deflection signals, the difference formulation 14 above does not allow such an interpretation. However, it will be shown, that this is no major drawback as long as the dominant signal contents are known.

## 1.3 Higher-Order Least-Square Fit Formula- tion

Basis of the concepts is the following expansion of the dynamic deflection $\Phi'(t)$:

$$ \Phi'_{i}(t) = \sum_{k=1}^{N_{\text{modes}}} a_{k,i}^{H}(t) \cos(\omega_{0,k}t) + b_{k,i}^{H}(t) \sin(\omega_{0,k}t) \\ + \sum_{k=1}^{N_{\text{EQ}}} a_{k,i}^{P}(\Omega) \cos(k\varphi(t)) + b_{k,i}^{P}(\Omega) \sin(k\varphi(t)) \quad (15) $$

The expansion above is the general solution of the system of ordinary differential equations with integral excitation. Whereas the homogenous solution (superscript H) denotes the transient solution with usually damped vibrations in the natural frequencies of the system, the particular solution (superscript P) denotes the forced response with respect to the engine orders present in the regarded application.

Next we assume, that either the vibration frequency (natural frequency of the mode of interest) in the non-integral case, or the engine orders of interest in the integral case are known. Otherwise there are other studies available dealing with the question of the determination of vibration frequency [4, 6, 7]. For the reason of clarity, the higher order least-square fit is introduced for the integral case with a single engine order present and the superscripts are left out:

$$ \Phi'_{i}(t) = a_{k,i}(\Omega) \cos(k\varphi(t)) + b_{k,i}(\Omega) \sin(k\varphi(t)) \quad (16) $$

The extension to multiple engine orders and/or non-integral case works analogously by adding more unknowns or shifting from frequency dependency to time dependency.

### Classical Revolution Based Approach

The easiest way to determine the unknown FRFs $a_{k,i}(\Omega)$ and $b_{k,i}(\Omega)$ is to split the signal into short analysis blocks, e.g. a certain number of revolutions $N_{\text{rev}}$, in which the rotational speed and thus the FRFs can be assumed as constant. By decreasing the analysis block length, the amplitude accuracy can be improved. However, less analysis points increase the susceptibility with respect to measurement noise and other frequency components, since the signals are heavily undersampled. On the other hand, by increasing the block length responses with low signal to noise ratios may be extracted,
though the amplitude accuracy will be poor. Thus, the choice of a suitable block length is always a trade-off between amplitude accuracy and robustness. Inserting 16 with constant $a_{k,i}$ and $b_{k,i}$ in 14 yields a system of linear equations to determine the unknown amplitude of blade i for the regarded engine order k:

$$ 
\begin{gathered}
\begin{bmatrix}
\Delta c_{i,1,l} & \Delta s_{i,1,l} \\
\vdots & \vdots \\
\Delta c_{i,j,l} & \Delta s_{i,j,l} \\
\vdots & \vdots \\
\Delta c_{i,N_p,l} & \Delta s_{i,N_p,l} \\
\Delta c_{i,1,l+1} & \Delta s_{i,1,l+1} \\
\vdots & \vdots \\
\Delta c_{i,N_p,l+N_{rev}} & \Delta s_{i,N_p,l+N_{rev}}
\end{bmatrix}
\begin{bmatrix}
a_{k,i} \\
b_{k,i}
\end{bmatrix}
=
\begin{bmatrix}
\Delta r_{i,j,l} \\
\vdots \\
\Delta r_{i,j,l+N_{rev}}
\end{bmatrix} \\
\Delta c_{i,j,l} = \begin{cases} \cos(k\varphi(t_{i,1,l+1})) - \cos(k\varphi(t_{i,j,l})) & \text{for } j = N_p \\ \cos(k\varphi(t_{i,j+1,l})) - \cos(k\varphi(t_{i,j,l})) & \text{else} \end{cases} \\
\Delta s_{i,j,l} = \begin{cases} \sin(k\varphi(t_{i,1,l+1})) - \sin(k\varphi(t_{i,j,l})) & \text{for } j = N_p \\ \sin(k\varphi(t_{i,j+1,l})) - \sin(k\varphi(t_{i,j,l})) & \text{else} \end{cases} \\
\Delta r_{i,j,l} = \begin{cases} \Theta_1^0 - \Theta_j^0 - (\varphi(t_{i,1,l+1}) - \varphi(t_{i,j,l})) \\ +\Delta\Theta_1 - \Delta\Theta_j & \text{for } j = N_p \\ \Theta_{j+1}^0 - \Theta_j^0 - (\varphi(t_{i,j+1,l}) - \varphi(t_{i,j,l})) \\ +\Delta\Theta_{j+1} - \Delta\Theta_j & \text{else} \end{cases}
\end{gathered}
$$

The classical revolution based approach above has the advantage that the data structure can be easily reused and assessed within the fitting process. However, especially in the case of distorted speed signals with no constant acceleration the approach can produce many overshoots and multiple amplitudes fitted for the same frequency. This is due to the fact, that no ramp like accelerations lead to a significant increase of the transients ($a_{k,i}^H(t)$ and $b_{k,i}^H$ in eq. 15) in the signal distorting the fit of the particular solution. If the steady state response is the quantity of interest, there is no physical reason to expect different amplitudes for a given excitation frequency. Thus, the following sections introduce a slightly different concept that also allows a higher order fit to be implemented.

### Piecewise Constant Approach

The basic idea is, that the amplitude and phase of the FRF is a unique and smooth function of the excitation frequency. Let the sine and cosine part of this FRF be approximated by finite elements of different order. Therefore, the frequency domain (or time domain in terms of non-integral excitation) is discretised by a number of finite elements $m = 1, N_{el}$, each of them covering an interval with a left and right frequency border $[\Omega_{0,m}, \Omega_{1,m}]$. The union of all elements cover the whole frequency domain and the intersection of the intervals of any two elements are empty, i.e. the intervals do not overlap. The set of triplets of indices (i,j,l) belonging to the collected data points is split into several subsets $\mathcal{I}_{i,m}$:

$$ \mathcal{I}_{i,m} = \{ (i, j, l) \mid \Omega_{0,m} \le \Omega(t_{i,j,l}) < \Omega_{1,m} \} \quad (17) $$

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Sketch of a typical SDOF FRF and results of the different fitting techniques</figcaption></figure>

The local coordinate $\xi_{i,j,l} \in [0,1]$ indicates the relative position within the finite element and is defined as:

$$ \xi_{i,j,l} = \frac{\Omega(t_{i,j,l}) - \Omega_{0,m}}{\Omega_{1,m} - \Omega_{0,m}} \forall (i,j,l) \in \mathcal{I}_m $$

The higher order formulation will be introduced by specific shape functions $f_{\cdot}(\xi)$ varying with the targeted order and lead to the slightly modified elements in the least square matrix:

$$ 
\begin{align}
\Delta c_{i,j,l}^{(\cdot)} &= \begin{cases} f_{\cdot}(\xi_{i,1,l+1}) \cos(k\varphi(t_{i,1,l+1})) \\ -f_{\cdot}(\xi_{i,j,l}) \cos(k\varphi(t_{i,j,l})) \\ f_{\cdot}(\xi_{i,j+1,l}) \cos(k\varphi(t_{i,j+1,l})) \\ -f_{\cdot}(\xi_{i,j,l}) \cos(k\varphi(t_{i,j,l})) \end{cases} & \text{for } j = N_p \\
\Delta s_{i,j,l}^{(\cdot)} &= \begin{cases} f_{\cdot}(\xi_{i,1,l+1}) \sin(k\varphi(t_{i,1,l+1})) \\ -f_{\cdot}(\xi_{i,j,l}) \sin(k\varphi(t_{i,j,l})) \\ f_{\cdot}(\xi_{i,j+1,l}) \sin(k\varphi(t_{i,j+1,l})) \\ -f_{\cdot}(\xi_{i,j,l}) \sin(k\varphi(t_{i,j,l})) \end{cases} & \text{for } j = N_p \tag{18}
\end{align}
$$

The vectors $\Delta c_{i,m}^{(\cdot)}, \Delta s_{i,m}^{(\cdot)}$ and $\Delta r_{i,m}$ contain the $\Delta c_{i,j,l}^{(\cdot)}, \Delta s_{i,j,l}^{(\cdot)}, \Delta r_{i,j,l}$ of all triplets (i,j,l) of indices belonging to the subset $\mathcal{I}_{i,m}$.

Based on the introduced definitions the piecewise constant
approach can be formulated:

$$ \begin{bmatrix} \mathbf{A}_{i,1}^0 & 0 & \dots & 0 \\ 0 & \mathbf{A}_{i,2}^0 & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & \dots & 0 & \mathbf{A}_{i,N_{el}}^0 \end{bmatrix} \begin{bmatrix} \mathbf{b}_{i,1}^0 \\ \mathbf{b}_{i,2}^0 \\ \vdots \\ \mathbf{b}_{i,N_{el}}^0 \end{bmatrix} = \begin{bmatrix} \mathbf{r}_{i,1} \\ \mathbf{r}_{i,2} \\ \vdots \\ \mathbf{r}_{i,N_{el}} \end{bmatrix} \quad (19) $$

The local element matrices, unknown vectors and shape functions read:

$$ \begin{align} A_{i,m}^0 &= [\Delta\mathbf{c}_{i,m}^0 \quad \Delta\mathbf{s}_{i,m}^0] && \text{with } f^0(\xi) \equiv 1 \nonumber \\ \mathbf{b}_{i,m}^0 &= \begin{bmatrix} a_{k,m} \\ b_{k,m} \end{bmatrix} \tag{20} \end{align} $$

## Linear Spline Approach

The linear spline approach can be formulated:

$$ \begin{bmatrix} \mathbf{A}_{i,1}^{10} & \mathbf{A}_{i,1}^{11} & 0 & \dots & 0 \\ 0 & \mathbf{A}_{i,2}^{10} & \mathbf{A}_{i,2}^{11} & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & \dots & 0 & \mathbf{A}_{i,N_{el}}^{10} & \mathbf{A}_{i,N_{el}}^{11} \end{bmatrix} \begin{bmatrix} \mathbf{b}_{i,1}^{10} \\ \mathbf{b}_{i,2}^{10} \\ \vdots \\ \mathbf{b}_{i,N_{el}}^{10} \\ \mathbf{b}_{i,N_{el}}^{11} \end{bmatrix} = \begin{bmatrix} \mathbf{r}_{i,1} \\ \mathbf{r}_{i,2} \\ \vdots \\ \mathbf{r}_{i,N_{el}} \end{bmatrix} \quad (21) $$

The local element matrices, unknown vectors and shape functions read:

$$ \begin{align} A_{i,m}^{10} &= [\Delta\mathbf{c}_{i,m}^{10} \quad \Delta\mathbf{s}_{i,m}^{10}] && \text{with } f^{10}(\xi) = 1 - \xi \nonumber \\ A_{i,m}^{11} &= [\Delta\mathbf{c}_{i,m}^{11} \quad \Delta\mathbf{s}_{i,m}^{11}] && \text{with } f^{11}(\xi) = \xi \tag{22} \\ \mathbf{b}_{i,m}^{10} &= \begin{bmatrix} a_{k,m}^0 \\ b_{k,m}^0 \end{bmatrix}, \quad \mathbf{b}_{i,m}^{11} = \begin{bmatrix} a_{k,m}^1 \\ b_{k,m}^1 \end{bmatrix} && \equiv \mathbf{b}_{i,m+1}^{10} \quad \forall m \in [1, N_{el} - 1] \nonumber \end{align} $$

## Cubic Spline Approach

The cubic spline approach can be formulated:

$$ \begin{bmatrix} \mathbf{A}_{i,1}^{30} & \mathbf{A}_{i,1}^{31} & \mathbf{A}_{i,1}^{32} & \mathbf{A}_{i,1}^{33} & 0 & 0 & \dots & 0 \\ 0 & 0 & \mathbf{A}_{i,2}^{30} & \mathbf{A}_{i,2}^{31} & \mathbf{A}_{i,2}^{32} & \mathbf{A}_{i,2}^{33} & \dots & 0 \\ \vdots & 0 & 0 & \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & \dots & 0 & 0 & \mathbf{A}_{i,N_{el}}^{30} & \mathbf{A}_{i,N_{el}}^{31} & \mathbf{A}_{i,N_{el}}^{32} & \mathbf{A}_{i,N_{el}}^{33} \end{bmatrix} \times \begin{bmatrix} \mathbf{b}_{i,1}^{30} \\ \mathbf{b}_{i,2}^{30} \\ \vdots \\ \mathbf{b}_{i,N_{el}}^{30} \\ \mathbf{b}_{i,N_{el}}^{31} \end{bmatrix} = \begin{bmatrix} \mathbf{r}_{i,1} \\ \mathbf{r}_{i,2} \\ \vdots \\ \mathbf{r}_{i,N_{el}} \end{bmatrix} \quad (23) $$

The local element matrices, unknown vectors and shape functions read:

$$ \begin{align} \mathbf{A}_{i,m}^{30} &= [\Delta\mathbf{c}_{i,m}^{30} \quad \Delta\mathbf{s}_{i,m}^{30}] && \text{with } f^{30}(\xi) = 1 - 3\xi^2 + 2\xi^3 \nonumber \\ \mathbf{A}_{i,m}^{31} &= [\Delta\mathbf{c}_{i,m}^{31} \quad \Delta\mathbf{s}_{i,m}^{31}] && \text{with } f^{31}(\xi) = \xi - 2\xi^2 + \xi^3 \nonumber \\ \mathbf{A}_{i,m}^{32} &= [\Delta\mathbf{c}_{i,m}^{32} \quad \Delta\mathbf{s}_{i,m}^{32}] && \text{with } f^{32}(\xi) = 3\xi^2 - 2\xi^3 \nonumber \\ \mathbf{A}_{i,m}^{33} &= [\Delta\mathbf{c}_{i,m}^{33} \quad \Delta\mathbf{s}_{i,m}^{33}] && \text{with } f^{33}(\xi) = -\xi^2 + \xi^3 \tag{24} \end{align} $$

$$ \mathbf{b}_{i,m}^{30} = \begin{bmatrix} a_{k,m}^0 \\ a_{k,m}^0 \\ b_{k,m}^0 \\ b_{k,m}^0 \end{bmatrix}, \quad \mathbf{b}_{i,m}^{31} = \begin{bmatrix} a_{k,m}^1 \\ a_{k,m}^1 \\ b_{k,m}^1 \\ b_{k,m}^1 \end{bmatrix} \equiv \mathbf{b}_{i,m+1}^{30} \quad \forall m \in [1, N_{el} - 1] $$

<figure><img src="image_6.png" /><figcaption>Figure 3. Utilized shape functions</figcaption></figure>

## 2. PERFORMANCE OF THE APPROACH

As mentioned in the introduction, the performance of the approach will be tested with analytical BTT data of known signal content. These signals are generated based on a SNM representation of a research compressor blisk rotor. The ROM formulation is derived from a finite element model of single sector with cyclic symmetric boundary conditions. Figure 4 shows the FE mesh of studied rotor. Whereas the full finite element model contains ≈ 230000 degrees of freedom, the efficient ROM has only 28 (mode family 1 only). Hence, it is perfectly suited to generate analytical BTT data of high quality. Figure 4 shows also the Campbell Diagram of the rotor
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Detail of the mesh of the FE model of the studied rotor with mode shape of blade mode family 1 (left) and Campbell Diagram with most relevant upstream excitation sources (right)</figcaption></figure>

with the most relevant upstream excitation sources. The frequencies are normalised by the frequency of the maximum nodal diameter of mode family 1 at rest. The highlighted speed range from 80-100% mechanical speed contains three resonance crossing, EO 4 (IGV-S1) with mode 1, EO 36 (S1) with mode 10 and EO40 (IGV) with mode 11. The mode shape of mode family 1 is also depicted in fig. 4.

Considering the rich modal content in the higher frequency range of state of the art blisk rotors and the multiple excitation sources in multi stage compressors, the presence of vibrations of multiple modes can not be avoided. Whereas with the currently available bandwidth of modern telemetry systems the mode individual amplitudes in such a speed range can be analysed trustworthy with strain gauges, the undersampled BTT data makes it very difficult to distinguish between the different frequency content of the signal. Thus, it is a very good test case to show the performance of the introduced approach. For clarity reasons only two post-processing options are compared, the classical revolution based approach and the cubic spline approach. Both take advantage of the difference formulation.

## 2.1 Single Mode Response to Integral Excitation

This section deals with a first simple test case. Therefore, the structural model will only include the mode 1 degrees of freedom, i.e. there is only a single mode responding to a single EO - probably the easiest scenario BTT can face. Nevertheless, mistuning of the rotor is accounted for, such that each blade will vibrate with a different amplitude. The acceleration from 70% to 90% mechanical speed is sufficiently slow (≈ 0.8Hz/s), such that the transient response can be

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Transient and steady state vibration amplitude vs. frequency exemplarily shown for blade 2</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Raw deflection vs. time observed by the BTT for the noisy and clean shot</figcaption></figure>

regarded as steady state, see fig. 5. The simulated acceleration manoeuvre has a length of 60 seconds and is simulated within ≈ 40 seconds on a standard desktop computer. Aerodynamic damping, gyroscopic effects and static deflections are omitted in this first study.

Two different signals are studied, a perfect one and a signal with noise added to the time of arrivals. The added noise leads to heavily distorted raw deflection signals with a signal to noise ratio of ≈ 0.5, see fig. 6. Figure 7 shows the mean error in the predicted vibration amplitude as a function of the averaging width for the clean and noisy shot. In terms of the revolution based approach the averaging width is controlled with the number of revolutions chosen for the fit. In terms of the cubic spline based approach the averaging width is controlled with the number of finite elements used. The averaging width has been normalised by the half-power width calculated from the applied damping. Comparing the clean case, a remarkable increase of the possible averaging width can be noticed when using the cubic spline approach compared with revolution based approach, which was intended by the introduction of this formulation. The cubic spline approach allows to average over 80% of the half power width without losing any accuracy in the predicted amplitude, whereas the revolution based approach starts to suffer significantly when the averaging width exceeds 30%. Thus, the possible averaging width could be increased by almost a factor of 3, which is assumed to be very beneficial in difficult cases.

The first difficult case studied is the very noisy deflection signal introduced above. With the noise being a factor of 2 higher than the actual deflection, the resonance can not be detected anymore from the raw deflection, see fig. 6. Nevertheless, providing a suitable averaging width, both algorithms are able to extract the vibration amplitude within
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Mean relative error in predicted amplitude vs. averaging width in percent of the half power width for the clean and noisy shot</figcaption></figure>

an error range < 5%, see fig. 7. The figure also shows the typical trade-off that has to be made when real engine data is analysed. Choosing a too short averaging width, the algorithms suffer from irrelevant components contained in the signal. With a very large averaging width, the algorithms loose amplitude accuracy.

Figure 8 shows how this amplitude error is reflected in the predicted FRFs. For a low averaging width the FRF predicted by the revolution based approach, compare fig. 8 (top), is only very hardly identified with the analytical one. Such a FRF would nor be usable for the purpose of damping identification or structural model validation. With a higher averaging width, compare fig. 8 from top to bottom, the predicted FRF is smoothed and appears to be close to the analytical one. It is also remarkable, that the FRF predicted by the cubic spline approach does not degenerate to the same extent as the ones predicted by the classical approach. Each of them reflects the qualitative structural behaviour properly and would be applicable for a damping quantification even though the predicted damping values would differ much.

### 3. CONCLUSIONS

A higher order least square fit for the assessment of integral and non-integral vibrations with blade tip timing has been introduced. The approach takes advantage of the difference formulation for the time of arrivals, which makes it more robust since no static deflection needs to be modelled and the blade positioning error cancels out inherently. The introduction of the approach addresses the unavoidable trade-off between robustness and amplitude accuracy that needs to be made with the choice of a suitable averaging width. In theory, the introduced approach does allow much larger averaging widths than the classical revolution based approach without loosing amplitude accuracy and should therefore smooth out errors more efficiently. Hence, it is assumed to be superior to the classical algorithms in hard cases with multiple modes responding or when higher modes with low signal to noise ratios are analysed. The effectiveness of the approach has been demonstrated by analytically generated BTT data with known signal content. The generated data reflects real engine data since it is derived from a high fidelity reduced order model that is based on a finite element model

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Measured (noisy shot) and analytical FRF of blade 2 for a low (top), the optimum (middle) and a large (bottom) averaging width</figcaption></figure>

of a compressor rotor and incorporates mistuning effects. It has been demonstrated, that for the particular test case the cubic spline based approach allows an averaging width being almost a factor of 3 larger compared to the classical approach. It has been shown, that providing the optimum averaging width both approaches are able to predict the vibration amplitude within a range of ±5% despite a noise level being twice as high as the vibration level. However, the predicted FRFs of the cubic spline approach did not degenerate to the same extent as the ones predicted by the classical approach. Hence, the introduced approach is assumed to supplement the state of the art in non-intrusive vibration measurements with blade tip timing.

### REFERENCES

[1] G Schlagwein and U Schaber. Non-contact blade vibration measurement analysis using a multi-degree-of-freedom model. *Proceedings of the Institution of Me-
*chanical Engineers, Part A: Journal of Power and Energy,*
220(6):611–618, 2006.

[2] M Zielinski and G Ziller. Noncontact vibration measurements on compressor rotor blades. *Measurement Science and Technology*, 11(7):847, 2000.

[3] Michael Zielinski, Gerhard Ziller, et al. Noncontact blade vibration measurement system for aero engine application. *Proceedings of the International Conference on Air Breathing Engines: Paper ISABE-2005-1220*, 2005.

[4] Jon Gallego-Garrido, Grigorios Dimitriadis, and Jan Robert Wright. A class of methods for the analysis of blade tip timing data from bladed assemblies undergoing simultaneous resonances—part i: Theoretical development. *International Journal of Rotating Machinery*, 2007, 2007.

[5] Jon Gallego-Garrido, Grigorios Dimitriadis, Ian B Carrington, and JR Wright. A class of methods for the analysis of blade tip timing data from bladed assemblies undergoing simultaneous resonances—part ii: Experimental validation. *International Journal of Rotating Machinery*, 2007, 2007.

[6] Ian B Carrington, Jan Robert Wright, Jonathan E Cooper, and Grigorios Dimitriadis. A comparison of blade tip timing data analysis methods. *Proceedings of the Institution of Mechanical Engineers, Part G: Journal of Aerospace Engineering*, 215(5):301–312, 2001.

[7] Grigorios Dimitriadis, Ian B Carrington, Jan Robert Wright, and Jonathan E Cooper. Blade-tip timing measurement of synchronous vibrations of rotating bladed assemblies. *Mechanical Systems and Signal Processing*, 16(4):599–622, 2002.

[8] W. Robinson and R. Washburn. A real time non-interference stress measurement system (nsms) for determining aero engine blade stresses. In *Proceedings of the 37th International Instrumentation Symposium*, 1991.

[9] O. Jousselin, P. Russhard, and P. Bonello. A method for establishing the uncertainty levels for aero-engine blade tip amplitudes extracted from blade tip timing data. In *10th International Conference on Vibrations in Rotating Machinery*, 2012.

[10] D Knappett and J Garcia. Blade tip timing and strain gauge correlation on compressor blades. *Proceedings of the Institution of Mechanical Engineers, Part G: Journal of Aerospace Engineering*, 222(4):497–506, 2008.

[11] Yang, M.-T.; Griffin, J.H. A reduced-order model of mistuning using a subset of nominal modes. *ASME J. of Engineering for Gas Turbines and Power*, Vol. 123:893–900, October 2001.

[12] T. Giersch, P. Hönisch, B. Beirow, and A Kühhorn. Forced response analyses of mistuned radial inflow turbines. *Journal of Turbomachinery*, 135(3), 2013.