Zhen Wang

Van Nghi

François Martel
email: francois.martel2@usherbrooke.ca

Denis Rancourt
email: denis.rancourt@usherbrooke.ca

Catherine Chochol
email: catherine.chochol@insa-lyon.fr

Yves St-Amant
email: yves.st-amant@gmc.ulaval.ca

Didier Rémond
email: didier.remond@insa-lyon.fr

Time-varying torsional stiffness identification on a cantilever beam using Chebyshev polynomials

Keywords: . . . . . . . . . . Linear time-varying systems, parameter identification, Chebyshev polynomials, dynamic systems, experimental mechanical systems, modal analysis. Highlights Linear time-varying systems, parameter identification, Chebyshev polynomials, dynamic systems, experimental mechanical systems, modal analysis

This work was led in the LaMCoS laboratory, in the DCS team. I would like to thank especially my two Phd supervisors, Simon Chesne and Didier Remond, for their help, their receptiveness, emboldening and their faith in the work management. It was a great pleasure to present

Abstract in french Introduction

L'identification de structure est devenu aujourd'hui un enjeu majeur pour le contrôle et la surveillance de celle-ci. L'exemple le plus évocateur est celui de l'aéronautique où aucune transgression sur la sureté de fonctionnement ne peut être faite. Ceci engendre malheureusement des coûts de maintenance non négligeables puisque, pour atteindre le niveau de sécurité souhaité, on change les composants mécaniques de manière systématique après un nombre de cycles défini. Le suivi de ce type de composants permettrait de détecter à quel moment l'opération de maintenance est la plus adéquate, et de ce fait permettrait d'augmenter la durée de vie de ces pièces sans compromis sur la sureté de fonctionnement. De plus, l'espacement des opérations de maintenance permettrait d'augmenter de fait la durée d'utilisation de l'appareil ou du moins diminuerait considérablement le temps où l'appareil resterait inutilisable. Aujourd'hui le développement de capteurs intégrés, tels que la fibre composite permet d'introduire des réseaux de capteurs dans la structure même du composant à suivre. Le but de cette thèse est d'initier la réflexion sur une manière d'optimiser la forme de ce réseau de capteurs afin de reconstruire les informations recherchées :

-Vieillissement de matériau, par le suivi de paramètres matériau tel que module d'Young, densité volumique, etc, -Détection de défaut, par la détection de discontinuité tel que fissures, inclusions, etc, -Caractérisation de structures et de matériau, par l'estimation du meilleur modèle de comportement pour une structure inconnue Différents procédés d'identification ont été développés au cours des dernières décennies. On citera ici uniquement ceux dont le fonctionnement se rapproche le plus de l'approche proposée dans cette thèse. La méthode de reconstruction modale a suscitée un grand intérêt au sein de la communauté mécanique. Différentes applications peuvent être associées à cette famille telles que la corrélation en onde inverse [1,2], la décomposition en onde inverse [3,4] ou la méthode de la densité modale [5,6]. Les travaux les plus complets sur ces méthodes sont présentés par Ewins [7]. Le principal inconvénient de ces méthodes est qu'il est assez complexe d'estimer les paramètres de systèmes non linéaires. On notera entre autre la méthode RIFF (résolution inverse filtrée fenêtrée) basée sur le même principe : l'estimation des dérivées d'un signal acquis expérimentalement afin de reconstruire l'équation de comportement de la structure. Cette méthode développée par Pezerat [8,9,10,11] a été appliqué à l'identification d'une plaque amortie par [12], ce qui rejoint les travaux proposés dans cette thèse. Dans le domaine de l'automatique, de nombreux outils d'identification ont été développés, iii entre autre des outils d'identification paramétrique [START_REF] Garnier | Identification of continuous-time models from sampled data[END_REF](comme les modèles ARX, Box-Jenkins, etc.). Ces méthodes semblaient peu adaptées au domaine de la mécanique, les paramètres reconstruits n'ayant aucune signification mécanique. D. Remond [14] a adapté quant à lui la méthode d'identification à temps continu au domaine de la mécanique, en identifiant les paramètres de modèles discrets tel que la masse, la raideur, l'amortissement. Cette méthode est basée sur la projection d'un signal sur une base orthogonale, de l'estimation des dérivées de cette projection, grâce à la projection, pour enfin reconstruire l'équation de comportement du système discret. La méthode d'identification développée dans cette thèse est inspiré des travaux de D. Remond. On considérera les données d'entrée suivante :

-la réponse de la structure, qui sera mesurée de manière discrète, et qui dépendra des dimensions de la structure (temps, espace) -le modèle de comportement, qui sera exprimé sous forme d'une équation différentielle ou d'une équation aux dérivées partielles, -les conditions aux limites ainsi que la source d'excitation seront considérées comme non mesurées, ou inconnues.

Principe d'identification

La procédure d'identification est composée de trois étapes :

la projection sur une base polynomiale orthogonale du signal mesuré, -la différentiation du signal mesuré, -l'estimation de paramètres, en transformant l'équation de comportement en une équation algébrique. Les outils développés pour les trois étapes citées ci-dessus seront développées dans les prochains paragraphes.

Projection

Le signal est projeté sur la base de polynômes de Chebyshev. Cette famille de polynômes a été choisie pour ses bonnes propriétés (développées ci-après). On définit la projection comme il suit :

f (x) = ∞ i=0 λ i P i (x) (1)
f étant le signal, λ i le coefficient de projection et P i le polynôme d'ordre i. Dans la pratique il est impossible d'estimer la projection exacte de f , puisqu'il est impossible d'estimer une somme infinie. De ce fait, on calcule la meilleure estimation de f sur une base de taille N , au sens des moindres carrés. N sera également nommé l'ordre de troncature. On aura donc :

f (x) ≈ N i=0 λ i P i (x) (2)
Les coefficients λ i seront estimés par le produit scalaire associé à cette base :

λ i =< f (x), P i (x) > = Γ w(x)f (x)P i (x)dx (3)
Γ étant l'intervalle d'orthogonalité, w(x) la fonction poids. Dans la pratique il est impossible d'estimer cet intégrale de façon exacte. On peut l'estimer approximativement en utilisant la méthode des trapèzes. Malheureusement celle-ci a ses limites, bien connues. Afin de minimiser l'erreur il est essentiel de connaître f sur un grand nombre de points x, ce qui est très complexe dans la pratique lorsque ces données représentent des mesures. Les polynômes de Chebyshev possèdent une intéressante propriété de projection qui possède de l'intérêt par rapport à notre approche. La projection discrète sur les points de Gauss :

λ i = N +1 j=1 f (x j )P i (x j ) (4)
avec N + 1 > i et x j pour j = 1..N + 1 sont les zéros (points de Gauss) du polynôme P N +1 . Cette projection discrète permet d'estimer de manière exacte les coefficients de projection avec uniquement N + 1 points.

Différentiation

On peut exprimer les dérivés des polynômes de Chebyshev comme une somme de ces mêmes polynômes. On écrit donc communément : T 0 (x) . . .

f (x) ≈ N i=0 λ i d dx T i (x) ≈ {λ} [D] {T } (5)
T N (x)    et [D] = 2 ×              0 0 0 0 0 0 • • • 0 1 0 0 0 0 0 • • • 0 0 4 0 0 0 0 • • • 0 3 0 6 0 0 0 • • • 0 0 8 0 8 0 0 • • • 0 5 0 10 0 10 0 • • • 0 . . . . . . . . . . . . . . . . . . . . . . . . 0 2N 0 2N 0 2N • • • 0              (6)
Lui [15] a prouvé l'erreur commise lors de l'écriture de ce paramètre. Cette erreur a également été observée au cours des applications numériques lors de cette thèse. Il a donc été nécessaire de développer un nouvel opérateur de différentiation. Ce nouvel opérateur est nommé [∆]. Il est basé sur la théorie suivante : on exprime le produit scalaire de la dérivée :

λ α,γ i =< f (α) • u γ , P i (x) > = Γ f (α) (x) • u γ (x) • w(x) • P i (x)dx (7)
En intègrant par partie, on obtient l'expression suivante :

λ α,γ i = [f (α-1) (x) • u γ (x) • w(x) • P i (x)] Γ - Γ f (α-1) (x) • (u γ (x) • w(x) • P i (x)
) (1) dx (8) v En choisissant la fonction u γ de manière à annuler la partie intégrée, on obtient :

{λ α,γ } = [∆ α,γ ]{λ 0,γ-α } (9)
L'opérateur [∆] permet donc d'exprimer de manière exacte la projection de la dérivée d'ordre α en fonction de la projection du signal.

Estimation de paramètres

On cherche à estimer les paramètres de l'équation de comportement. Cette étape sera présentée à l'aide de l'exemple de la poutre de Bernoulli pour sa simplicité et son ordre de dérivation élevé :

∂ 4 ∂x 4 v(x, t) = ρS EI ∂ 2 ∂t 2 v(x, t) (10)
avec v le déplacement transversal dépendant de l'espace x et du temps t. ρS EI est le paramètre structurel et matériau à identifier.En multipliant cette équation par la fonction u γ sus citée on obtient :

u(x, t) ∂ 4 ∂x 4 v(x, t) = ρS EI u(x, t) ∂ 2 ∂t 2 v(x, t) (11)
Cette équation aux dérivées partielles sera projetée sur la base des polynômes de Chebyshev : {λ 4,γx|0,γt } = ρS EI {λ 0,γx|2,γt }

les coefficients {λ 4,γx|0,γt } et {λ 0,γx|2,γt } sont les coefficients de projection des dérivées partielles de v. Ils seront donc estimés à l'aide de la projection de v ainsi que l'opérateur de différentiation [D] ou [∆].

On estimera ensuite ρS EI qui est la seule inconnue à l'aide des N ×N équations algébriques écrites. La dimension N × N correspond aux nombres de coefficients de projection sur la base polynomiale de taille N en espace et N en temps. On a plus d'équations que d'inconnues. On pourra donc estimer le paramètre par les moindres carrés. Malheureusement, lorsque le signal est entaché de bruit, l'estimateur des moindres carrés est biaisé. On ajoutera alors une étape supplémentaire de régularisation. Lors de cette thèse on a choisi d'adapter la théorie de la variable instrumentale issue du domaine de l'automatique. Cette méthode de régularisation permet de filtrer le signal bruité par le modèle choisi.

Applications numériques

Le processus d'identification a été appliqué à différents modèles numériques. La poutre de Bernoulli a permis d'établir un lien entre l'ordre de troncature de la base polynomiale et le nombre d'ondes contenu dans le signal projeté. Sur un signal bruité, nous avons pu établir une valeur de nombre d'onde et d'ordre de troncature minimum pour assurer une estimation précise du paramètre à identifier. Grâce à l'exemple de la poutre de Timoshenko, nous avons pu réadapter la procédure d'identification à l'estimation de plusieurs paramètres. Trois paramètres dont les valeurs ont des ordres radicalement différents ont été estimés. Cet exemple illustre également la vi stratégie de régularisation à adopter avec ce type de problèmes. L'estimation de l'amortissement sur une poutre a été réalisée avec succès, que ce soit à l'aide de sa réponse transitoire ou à l'aide du régime établi. Le cas bidimensionnel de la plaque a également été traité. Il a permis d'établir un lien similaire au cas de la poutre de Bernoulli entre le nombre d'onde et l'ordre de troncature. A titre d'exemple, la figure 3.26 présente les résultats d'identification pour une plaque, avec signal bruité, en fonction du nombre d'onde et de l'ordre de troncature. L'ordre de troncature est le paramètre de réglage de la méthode, le nombre d'onde permet de représenter le contenu fréquentiel du signal. L'échelle de couleur met en évidence les zones où la méthode d'identification est efficace (erreur entre 1% (vert) et 10 -9 % (bleu). Cette méthode met en évidence trois zones :

-En haut à gauche : ordre de troncature faible et nombre d'onde élevé. Dans cette zone la taille de la base ne permet pas de reconstruire correctement le signal -En bas : ordre de troncature élevé et nombre d'onde faible. Le contenu du signal n'est pas assez riche pour réaliser une identification non entaché par le bruit -En haut à gauche : ordre de troncature élevé et nombre d'onde élevé. Ici la précision est de l'ordre de 1%. Le signal est assez riche et la taille de la base polynomiale permet de le reconstruire

Applications expérimentales

Deux cas d'applications expérimentales ont été traités au cours de cette thèse. Le premier se base sur le modèle de la poutre de Bernoulli, appliqué à la détection de défaut. En effet on applique un procédé d'identification ayant pour hypothèse initiale la continuité de la structure. Dans le cas où celle-ci ne le serait pas on s'attend à observer une valeur aberrante du paramètre reconstruit. On recherche ici à retrouver une fissure créé artificiellement. Le procédé permet de localiser avec succès le lieu de la discontinuité. La figure 4.3 (a) présente la valeur du paramètre matériau/ géométrique reconstruit pour chaque tronçon de poutre. Au droit de la fissure, on observe un décrochage de cette valeur correspondant à la discontinuité. Figure 4.3 (b) on voit la dérivée d'ordre 4 en espace reconstruite à l'aide de la méthode. On observe au droit de la fissure une oscillation plus importante du signal. Cette oscillation est caractéristique de la présence d'une discontinuité dans le champ de déplacement qui ne respecte plus l'équation de comportement. Le second cas applicatif vise à reconstruire l'amortissement d'une structure 2D : une vii

Conclusion

En conclusion, les travaux de cette thèse auront pu mettre en évidence différentes problématiques telles que : viii -la stratégie de placement de capteurs lorsque l'on souhaite appliqué la méthode d'identification par les polynômes de Chebyshev, -l'erreur commise avec l'opérateur [D], utilisé alors communémant lors d'identification utilisant les polynômes de Chebyshev, -le développement d'un nouvel outil de différentiation [∆] performant, -l'apport d'une solution de régularisation pour corriger le biais de l'estimation par les moindres carrés, associé au bruit. Ces différentes améliorations apportées auront permis d'appliquer ce procédé d'identification à différents modèles numériques : poutre de Bernoulli, poutre de Timoshenko, plaque de Kirschhoff. Ces différents exemples permettent d'illustrer différentes stratégies d'identification. Les deux exemples expérimentaux illustrent parfaitement les différents domaines d'application de la méthode, comme la détection de défauts sur la poutre ou l'estimation de modèle sur la plaque polymère (PMMA). En perspectives, ces travaux pourraient être appliqué à la caractérisation de matériaux anisotropes.

ix

List of Figures

Introduction Context

The identification of mechanical systems and structures is currently attracting much attention for monitoring and control purposes. In civil and aeronautical engineering, for example, optimizing structure life-time is crucial. The principle aim of structure monitoring is to replace the periodical maintenance operations by maintenance operations controlled by monitoring the system. Hence, the mechanical elements life-time is larger without reducing the safety of the structure. The monitoring system should be able to foresee crack apparitions and to evaluate material ageing. This optimized maintenance permits cost reduction as elements are used longer and the stop of devices is made less frequently. Robust monitoring is also vital for updating controllers. Many control systems require accurate models of the systems controlled, thus the formulation of the parameter identification procedures implemented is a critical step. The development of sensing technologies creates also a strong need of these signal processing techniques. For example, sensors fibres, implemented in the structure material draw interesting perspectives in the development of smart structures. They generate a large amount of data, which requires a processing, in order to reveal the interesting information. The original goal of this study was to develop an adaptive sensing approach, which could help the design of these implemented sensors arrays.

Identification techniques for parameter identification

The main application domain of the identification problem treated in this work is vibrating structures and displacement field measurements post processing (source location, limit characterization, structure health monitoring, etc). As the problem of system identification has been largely treated, it becomes impossible to mention here all identification techniques which already exist. A focus on some identification techniques applied to mechanical systems will be presented here. The most developed techniques in the last decades are modal identification techniques. Different modal techniques exist, such as the inverse wave correlation [1,2], which is based on the determination of the curve dispersion using a measurement field, the inverse wave decomposition [3,4] which is based on the computation of the general solution of the equation, the modal density [5,6], which is based on the comparison of the modal density of a theoretical infinite plate and a modal density measured during experimentation. The most complete work on modal testing was led by D. J. Ewins [7]. The main drawback of the presented modal approaches is that it seems difficult to consider the identification of non-linear systems using a modal approach.

INTRODUCTION 2

The RIFF (Resolution Inverse Filtree Fenetree in French, meaning windowed filtered inverse resolution) technique must be mentioned as it identification objective is very close to the one developed in our work (estimation of signal derivatives and equation of motion reconstruction). RIFF technique applied on inverse problem was led by C. Pezerat [8,9,10,11]. This method considers the equation of motion of the structure. the partial derivatives of the structure response are estimated by a finite scheme. Using two previous considerations, the distribution forces in the structure are reconstructed. Since measured data are always noisy, the force distribution reconstructed is regularized applying a spatial window on the set of sensors used for the partial derivatives reconstruction. This method was successfully applied on damping identification by F. Ablitzer [12]. The main drawback of this method is the need of a regular and high resolution meshing of sensors on the considered structure. In the automation domain, some highly interesting identification tools has been developed. Various parametric identification methods [START_REF] Garnier | Identification of continuous-time models from sampled data[END_REF](such as ARX, Box-Jenkins, etc.) have been developed over the last decades. These methods seem difficult to implement in mechanical engineering, as the parameters identified have no mechanical meaning. Indeed, in both contexts (automation and mechanical), the inputs and goals seems different. Indeed in the automation context, the amount of data available for the identification is very large. In our structure identification context, only a few sensors are used and therefore led to different identification challenges. Moreover, the model are often more simple in the automation domain (for example with lower derivative order). Recently, D. Remond [14] proposed an improved continuous time identification method which permits the direct computation of mechanical quantities such as mass, damping, etc. This continuous-time identification method was applied on different multi-degrees of freedom systems [18,19,20] using different orthogonal bases (Chebyshev, Legendre, Fourier, etc.) and it has also been extended to non-linear systems [21] and time varying systems (see 4.2.3.2). As this method fulfils all the requirements presented previously, a focus on continous time identification is developed here.

Continuous time identification methods

Before the area of numerical measurements, in the 1950s and 1960s, continuous time models dominated the control world. Since then, these techniques were forgotten, perhaps because of the dominant interest in discrete-time identification. In the 21st century, scientists began to find new interest for using continuous-time identification. Different techniques were proposed, like the hybrid model [22] for example that could mix continuous time identification with the use of sampled data. The advantages of continuous time modelling are combined to the efficiency of noise modelling in discrete time identification. A reborn of continuous-time identification occurs around 1999. Different orthogonal basis are then studied and compared [23]. Chou, CT, Verhaegen, M. and Johansson, R suggest Laguerre function for SISO (single input single output) systems [24]. The identification proves to be precise even for a restrictive quantity of samples. It shows that this technique can potentially reduce the amount of information needed for identification. The subject is largely theoretically treated. A large part of the applications are still very experimental but shows the power of continuous-time identification. This technique has proved itself to be efficient for linear and non-linear identification [START_REF] Ricardo | On the identification of non-linear mechanical systems using orthogonal functions[END_REF]. Pacheco, RP and INTRODUCTION 3 Steffen, V have shown the efficiency of this technique applied to single degree and two degrees of freedom systems with mixed damping. The Legendre polynomials have been used for this study. Mechanical researchers have found an interest in this technique, by analogy with the equation that have been written by the electrical engineers. Mechanical laws have been used to construct the differential equations that convey its behaviour. It has permitted to estimate vibratory parameters of complex structures or to reconstruct the behaviour of a system. In 1976, Soderstrom had already the idea of using this identification method with mechanical equations of motion [START_REF] Kj Åström | Identification of ship steering dynamics[END_REF]. In practice, many identification problems relating to mechanical vibration result in a further problem of derivative estimation. These estimation methods concern various domains and provide numerous applications. The usual methods use finite difference schemes, but often require a regularisation step, such as for force localisation [START_REF] Leclere | Vibration source identification using corrected finite difference schemes[END_REF][START_REF] Pézerat | Two inverse methods for localization of external sources exciting a beam[END_REF]. Some approaches comprise a natural regularisation dimension that involves using an integral formulation for boundary characterisation [START_REF] Chesne | Distributed piezoelectric sensors for boundary force measurements in euler-bernoulli beams[END_REF][START_REF] Chesne | Identification of plate boundary forces from measured displacements[END_REF][START_REF] Mboup | Numerical differentiation with annihilators in noisy environment[END_REF], or polynomial approximation for damage detection [START_REF] Wu | Damage localization in plate structures from uniform load surface curvature[END_REF]. D. Wu [START_REF] Wu | Damage localization in plate structures from uniform load surface curvature[END_REF] adapted a part of this continuous time method for damage detection purposes on continuous structures. However, the whole continuous-time identification method has never been reformulated for a distributed parameter (continuous) structure, such as a beam or plate.

Problem statement

For this identification problem we consider the following inputs and challenges :

-The structure response : the structure response is measured discretely over a given observation window. The spacing and arrangement of the set of sensors is not limited. An attention should be paid on measurement perturbation. The identification method should be efficient in presence of noise -The structure model is considered as input of our problem. The aim of the method is to estimate the model parameters in order to obtain the best fit between the model and the real structure. The model structure is presented by the equation of motion (differential or partial differential equation) -The boundary conditions are considered as unknown. The external perturbations, the environment are unavailable. -The excitation source is not considered as an input of our system. The excitation was already efficiently used in denoising processes. The aim of this identification technique is to perform the parameter estimation without knowing the source (such as wind excitation, which is unmeasurable) The outputs of this identification techniques are the parameters of the considered structure. The structural parameters have a physical meaning and are associated to an expected structure contribution, such as damping, inertial effect. This work will be focused on the following goals :

-Model fitting : the identification process must be able to fit as far as possible a given theoretical model to a real, experimental structure. Depending on the scientific community, this identification problem is named indirect measurement or inverse problem solving. The choice of this model will not be treated in this work but is detailed in [START_REF] Ruzek | Experimental identification of the bending equation of beams from the vibration shape measurements[END_REF]. Experimentally, it is not possible to find an exact model of real structures. The most accurate models often have a large number of freedom degrees. The parameters of the models are hard to identify because they need a proportional amount of data from the real structure, and therefore a large number of sensors. Hence reduced models are used. The aim of model fitting is to adjust the model parameters in order to be as close as possible to the real structure behaviour. Model fitting is particularly interesting in the context of novel material behaviour characterization, such as composite material. -Structure ageing : the identification method, based on the structure response, can be used as real-time process and permits to follow the evolution of structure parameters. Indeed the identification process should be able to separate in the structure response the environmental contribution (wind, car excitations on a bridge, etc.) from the structure contribution (damping or inertial effect, etc.). -Damage detection: the main hypothesis used in the identification process is the continuity of the structure. A breach in this hypothesis such as a crack formation should induces some phenomena such as slope changes in transverse displacement, etc. These changes can be computed using the identification process.

Presentation of the work

The purpose of this work is to adapt and improve the continuous time identification method proposed by D. Remond [14] for continuous structures. D. Remond clearly separated this identification method into three steps: signal expansion, signal differentiation and parameter estimation. In this study, both expansion and differentiation steps are drastically improved. An original differentiation method is developed and adapted to partial differentiation. The first chapter is dedicated to theoretical considerations. The existing identification process is firstly adapted to continuous structure. Then the expansion and differentiation principle are presented. A focus on the novel differentiation technique is made and illustrated numerically. The second chapter gathers all testing considerations, such as noise addition, Monte Carlo testing and regularization step. The third chapter depicts different numerical applications. A focus is made on different practical particularities, such as the use of the steady-state response, the identification of multiple parameters, etc. The fourth chapter presents two experimental applications : the first application is a crack detection on a beam. The second application is the identification of damping on a plate. The whole work is summarized in a concluding chapter.

Introduction

In this chapter the three steps identification process will be firstly presented. The process will be described in the case of a differential equation and in the case of a partial differential equation. Then, in a second section, the Chebyshev polynomials will be presented. The expansion tool, which is needed in the first step of the identification process, is depicted. The differentiation tool, which is needed in the second step of the identification process, is then described. A well known differentiation tool is firstly presented ([D] operator), before our novel operator is developed ([∆]).

Basic Principle

In this section the identification technique is firstly presented for a mono-dimensional case (with a differential equation) and secondly for a multi-dimensional case (with a partial differential equation). As mentioned in the introduction, the aim of this method is to estimate parameters of a given model using the structure response. The model is the differential equation (or partial differential equation) and the parameters are the constants of this equation. This method is local. It could be applied on an entire structure or on a sub-part of the system. No assumption on the boundary conditions is made, therefore the method could be used with every type of unknown boundary conditions. The whole identification procedure will be firstly presented. The different developed tools are presented in a second section.

Differential equation transformation

As presented in [14], the Continuous Time Identification method is usually based on three steps:

-Expansion step Firstly, the recorded signals are expanded on a truncated orthogonal basis. The choice of the orthogonal basis is not limited. The expanded signals are reduced to a few expansion coefficients.

f (x) ≈ N i=0 λ i P i (x) (1.1)
The recorded signal f is expanded on orthogonal basis size N . N is also named the truncation order. This truncation order is a tuning parameter which is chosen by the user. λ i are the expansion coefficients. P i (x) is the i th orthogonal function of basis P . -Differentiation step

Secondly, the derivatives of the expanded signal are computed. For this step, the derivatives of the orthogonal functions are computed and derivatives of the signal are obtained from computed expansion coefficients.

f (z) (x) ≈ N i=0 λ z i P i (x) λ z i = ∞ k=0 δ z k λ k (1.2)
f (z) (the z th derivative of f ) is expanded on the orthogonal basis P . λ z i is the i th expansion coefficient of f (z) . The expansion coefficients λ z i are computed using the relationships between P i and its derivatives. These relationships are expressed by the δ z k coefficients which can be computed using various methods. -Differential to algebraic equation transform For this step, the differential equation governing the system behaviour is used. All the derivatives in the differential equation are replaced by their expansion (computed in the differential step). The result of this step is an algebraic equation composed of expansion coefficient arrays and parameters. The parameters are computed using this algebraic equation. The following is an example. Let the linear differential equation be:

f (4) (x) + A 3 f (3) (x) + A 2 f (2) (x) + A 1 f (1) (x) + A 0 f (x) = 0 (1.3)
This differential equation is composed of the derivatives of f and the A j , j = 0..3 parameters, it is transformed into the following algebraic equation:

∞ i=0 (λ 4 i + A 3 λ 3 i + A 2 λ 2 i + A 1 λ 1 i + A 0 λ i )P i (x) = 0 (1.4)
This equation must be true for all x, and for all i. Hence we can write :

λ 4 i + A 3 λ 3 i + A 2 λ 2 i + A 1 λ 1 i + A 0 λ i = 0 ∀i (1.5)
Using this set of algebraic equations (for i = 0..N ), the A j parameters can be computed using a simple least square method, with j ≤ N .

Partial differential equation transformation

For continuous structures (beam, plates, etc) the equation of motion can often be written as a partial differential equation (DPE). A general linear partial differential equation can be:

Z 1 z=0 α {1} z ∂ z f ∂x z 1 (x 1 , x 2 ) + Z 2 z=0 α {2} z ∂ z f ∂x z 2 (x 1 , x 2 ) + Z 1 z 1 =0 Z 2 z 2 =0 α {1,2} z 1 ,z 2 ∂ z 1 +z 2 f ∂x z 1 ∂x z 2 (x 1 , x 2 ) = g(x 1 , x 2 )
(1.6) The variable x 1 , x 2 could be either space or time variable, depending on the studied structure. For sake of simplicity, we will study here the free case, therefore with g(x 1 , x 2 ) = 0. In a further discussion, we will prove that choosing g(x 1 , x 2 ) = 0 is the most difficult identification problem and that the identification process is still reliable for g(x 1 , x 2 ) = 0. Let note the expansion of z 1 derivative in x 1 direction and z 2 derivative in x 2 direction :

λ z 1 |z 2 .
Computing the expansion of this equation we will obtain :

Z 1 z=0 α {1} z λ z|0 + Z 2 z=0 α {2} z λ 0|z + Z 1 z 1 =0 Z 2 z 2 =0 α {1,2} z 1 ,z 2 λ z 1 |z 2 = 0 (1.7)
λ z|0 , λ 0|z and λ z 1 |z 2 are estimated using the method developed in the next section. As for the differential equation, the identification technique is a three steps process :

-STEP 1 : we compute the expansion of f

λ i,j = f, P i (x 1 ), P j (x 2 ) (1.8)
-STEP 2 : we combine the λ i,j coefficients using the constants δ z 1 (i, j) and δ z 2 (i, j).

The constants δ z 1 (i, j) and δ z 2 (i, j) can be stored in matrices. These matrices are named ∆ z 1 and ∆ z 2 -STEP 3 : we rewrite the partial differential equation of motion :

{λ z 1 |z 2 } = [∆ z 1 ]{λ}[∆ z 2 ]
(1.9)

Chebyshev polynoms

The Chebyshev polynomials are a set of orthogonal polynomials defined as the solutions to the Chebyshev differential equation. They are used as an approximation to a least squares fit, and are a special case of the Gegenbauer polynomial with α = 0. They are also intimately connected with trigonometric multiple-angle formulas. In the beginning of this century they were used to find an approximate solution of differential equations [START_REF] Lakestani | The use of chebyshev cardinal functions for the solution of a partial differential equation with an unknown time-dependent coefficient subject to an extra measurement[END_REF][START_REF] Dale | The accurate solution of poisson's equation by expansion in chebyshev polynomials[END_REF][START_REF] Se El-Gendi | Chebyshev solution of differential, integral and integro-differential equations[END_REF]. The mathematical considerations generated by this research was very useful in our identification problem. Most of them are concentrated in Mason's book [START_REF] Mason | Chebyshev polynomials[END_REF]. Furthermore the Chebyshev polynomials can be used to construct a wavelet basis.

Expansion

Principle

The expansion on a orthogonal basis P is defined as

f (x) = ∞ i=0 λ i P i (x) (1.10)
The λ i are named the expansion coefficients. They can be computed using a defined scalar product:

λ i =< f (x), P i (x) > = Γ w(x)f (x)P i (x)dx (1.11)
w(x) being the weighting function and Γ the interval of orthogonality.

The expansion can be approximated to the N th order as follow:

f (x) ≈ N i=0 λ i P i (x) (1.12) 1.2.1.

First kind

The first kind Chebyshev polynomials are defined as follow:

T i (x) = cos(iθ) with cos(θ) = x (1.13)
The polynomials As shown in figure 1.1, all extrema of the polynomials are of equal magnitude. The zeros (Gauss points) are clustered towards the end points x ± 1. The scalar product of a function f by T i can be written as :

T i are orthogonal in the domain [-1 1]. -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 x T i (x) i=0 i=1 i=2 i=3 i=4
λ i =< f, T i >= 1 -1 f (x)w(x)T i (x)dx (1.14) with w(x) = (1 -x 2 ) -1/2
the weighting function associated to this basis. As used by D. Wu [START_REF] Wu | Damage localization in plate structures from uniform load surface curvature[END_REF], an optimal sensors positioning permits the exact estimation of the expansion coefficients using only a few discretely recorded data. The discrete orthogonality can be written as (see [START_REF] Mason | Chebyshev polynomials[END_REF]) :

λ i = N +1 j=1 f (x j )T i (x j ) (1.15)
with N + 1 > i and x j for j = 1..N + 1 are the zeros (or Gauss point) of the Chebyshev function T N +1 . A simple presentation of the Gauss points for j = 6 and j = 7 is given in the figure 1.2. The computation of the novel differentiation operator is based on the relationship:

d dx T n (x) = n 2 T n-1 (x) -T n+1 (x) 1 -x 2 (1.16)
and the recurrence formula :

T n (x) = 2xT n-1 (x) -T n-2 (x) (1.17)

Second kind

The second kind Chebyshev polynomials have been studied because of their interesting link with the compressed sensing approach (see next chapter 2). The second kind Chebyshev polynomials are defined as follow:

U i (x) = sin((i + 1)θ) sinθ with cos(θ) = x (1.18)
The polynoms U i are orthogonal in the domain As shown in figure 1.3, the extrema of the polynomials increase monotonically from the centre towards the ends. The scalar product of a function f by U i can be written as :

[-1 1]. -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 -4 -3 -2 -1 0 1 2 3 4 5 x U i (x) i=0 i=1 i=2 i=3 i=4
λ i =< f, U i >= 1 -1 f (x)w(x)U i (x)dx (1.19)
with w(x) = (1 -x 2 ) 1/2 the weighting function associated to this basis. The discrete orthogonality can be written as (see [START_REF] Mason | Chebyshev polynomials[END_REF]) :

λ i = N +1 j=1 f (x j )U i (x j )w(x j ) (1.20)
with N + 1 > i et x j pour j = 1..N + 1 are the zeros (or Gauss point) of the Chebyshev function T N +1 of order N + 1. The computation of the novel differentiation operator are based on the recurrence relationship:

U n (x) = 2xU n-1 (x) -U n-2 (x) (1.21)
and the differentiation formulae :

(1 -x 2 )U (1) n (x) -xU n (x) = (n + 1)T n+1 (x) (1.22) [(1 -x 2 ) 1/2 U n (x)] (1) = -(n + 1)(1 -x 2 ) -1/2 T n+1 (x) (1.23)
The relationship between the first and second kind Chebyshev basis is:

U n (x) -U n-2 (x) = 2T n (x) (1.24) 1.2.2 Differentiation 1.2.2.1 Operator [D]
As explained in [START_REF] Mason | Chebyshev polynomials[END_REF], each Chebyshev polynomial derivative can be expressed as a combination of polynomials in the same basis :

d dx T i (x) = 2i i-1 r=0 i-r odd T r (x) (1.25)
Therefore, we can always express the derivatives of an expanded signal as a combination of the signal expansion. This property is avalaible for any type of polynomial basis. If we consider the derivative of f as a differentiation of (1.12), we have:

f (x) ≈ N i=0 λ i d dx T i (x) ≈ {λ} [D] {T } (1.26) with {λ} =    λ 0 . . . λ N    {T } =    T 0 (x) . . . T N (x)    and [D] = 2 ×              0 0 0 0 0 0 • • • 0 1 0 0 0 0 0 • • • 0 0 4 0 0 0 0 • • • 0 3 0 6 0 0 0 • • • 0 0 8 0 8 0 0 • • • 0 5 0 10 0 10 0 • • • 0 . . . . . . . . . . . . . . . . . . . . . . . . 0 2N 0 2N 0 2N • • • 0              (1.27)
However the approximation made using (1.26) is not accurate, as proved in [15]. Indeed the product of [D] by {λ} is not the best least square estimate of f in the Chebyshev basis of size N . As explained by D-Y Lui, the estimation of the derivative is corrupted by term errors that originate from the truncation of the expansion of f . Using this formula, we do not obtain the best least square estimate of the derivative, which induces dramatic identification errors. In Rouby's work [16], it has been shown that after an optimal truncation number, the error made on the derivative estimation increases, then decreases again, etc, as shown in figure 1.4. The explanation of this error increase can be simply illustrated through this example. Let take f defined as :

f (x) = {λ}{T } =               3 -1 1 -1 0 0 0 0 -8 • 10 -2               T               T 0 (x) T 1 (x) T 2 (x) T 3 (x) T 4 (x) T 5 (x) T 6 (x) T 7 (x) T 8 (x)               (1.28)
We can neglect the 8 th order coefficient. Therefore the truncated expansion of f is :

f trunc (x) = {λ trunc }{T } =     3 -1 1 -1     T     T 0 (x) T 1 (x) T 2 (x) T 3 (x)     (1.29)
The exact value of the truncated derivative expansion is:

{ (1) λ} ≈ -4 2.72 -6 -1.28 0 -1.28 0 -1.28 0 (1.30)
And computing the derivative with the operator [D] we obtain:

{ (1) λ} = -4 4 -6 0 (1.31)
The difference between the two derivative expansions are:

-The five truncated coefficients (from order 4 to 8) which is the part of information lost during the truncation of f -An error is made on the first and third order expansion coefficients. This error is crippling for our identification procedure, as the expansion is not the best least square estimate of order 3 in our polynomial basis. This difference is illustrated in figure 1.5.

-1 -0.5 0 0.5 1 1.5 x f (1) (b) f (1) f (1) trunc f (1) with [D] Figure 1.5: f and its truncated expansion (a), f (1) its truncated estimation and its estimation with [D]

1.2.

Principle of novel operator [∆]

In this section, we propose a novel operator [∆] which is not corrupted by the numerical error as [D]. Indeed this operator is based on the exact estimate of the derivative scalar product. Let take the example of the α derivative of f , noted f (α) . We chose to estimated the expansion of f (α) • u γ . The purpose of u γ will be explained just after. Let write the following scalar product, where Γ = [-1, 1]:

λ α,γ i =< f (α) • u γ , P i (x) > = Γ f (α) (x) • u γ (x) • w(x) • P i (x)dx (1.32)
f (α) is considered as unknown whereas f , u γ , w and P i are known analytically and differentiable. Let integrate by part the previous expression :

λ α,γ i = [f (α-1) (x) • u γ (x) • w(x) • P i (x)] Γ - Γ f (α-1) (x) • (u γ (x) • w(x) • P i (x)) (1) dx (1.33)
Choosing u γ (x) equal to zero at Γ boundaries allows the integrated part to disappear. Choosing u γ (x) accurately and repeating the above operation α times (α being the derivative order), we will finally find:

λ α,γ i = (-1) α Γ f (x) • (u γ (x) • w(x) • P i (x)) (α) dx (1.34)
With Chebyshev polynomials, the derivative (u γ (x) • w(x) • P i (x)) (α) can be computed analytically. The aim of next section is to express this derivative.

First kind

We remind that for the first kind polynomials P i (x) = T i (x). In order to achieve the condition exposed before (u γ (x) equal to zero at Γ boundaries) we will choose :

u γ (x) = (1 -x 2 ) γ γ > α (1.35)
with α the differentiation order. The proof for the computation of the constants for α = 1 will be developed here. The other constants are given in 4.2.3.2. Let consider the scalar product between f • u γ and T n :

λ 1,γ i =< f (1) • u γ , T i > = Γ f (1) (x)w(x)u γ (x)T i (x)dx (1.36)
Integrating by part (1.40), we have:

λ 1,γ i = [f (x)w(x)u γ (x)T i (x)] Γ - Γ f (x)(w(x)u γ (x)T i (x)) (1) dx = - Γ f (x)(u γ (x)T i (x)) (1) dx (1.37)
if we consider f (1),f (-1), T i (1) and T i (-1) as finite values. Using (1.16) we have:

(u γ (x)T i (x)) (1) = (1 -x 2 ) γ-1/2 T (1) i (x) -2(γ - 1 2 )x(1 -x 2 ) γ-3/2 T i (x) = (1 -x 2 ) γ-3/2 [ i 2 (T i-1 (x) -T i+1 (x)) -2(γ - 1 2 )xT i (x)] (1.38)
Using the recurrence formula (1.17), we obtain:

(u γ (x)T i (x)) (1) = (1 -x 2 ) γ-3/2 [ i 2 (T i-1 (x) -T i+1 (x)) -(γ - 1 2 )(T i+1 (x) + T i-1 (x))] = (1 -x 2 ) γ-3/2 [( i 2 -γ + 1 2 )T i-1 (x) + (- i 2 -γ + 1 2 )T i+1 (x)]
(1.39)

Finally we have :

λ 1,γ i = Γ f (x)w(x)u γ-1 (x)[δ 1,γ 1 (i, i + 1)T i+1 (x) + δ 1,γ 1 (i, i -1)T i-1 (x)]dx (1.40) with δ 1,γ 1 (i, i+1) = -(-i 2 -γ + 1 2 ) and δ 1,γ 1 (i, i-1) = -( i 2 -γ + 1 2 ). The δ 1,γ 1 (i, j) constants are stored in a matrix named [∆ α,γ 1 ]
. The expansion combination giving the derivatives expansion can be written as:

{λ α,γ } = [∆ α,γ 1 ]{λ 0,γ-α } (1.41)
The constants in the matrix [∆ α,γ 1 ] are detailed for α = 1..4 in 4.2.3.2.

Second kind

We remind that for the second kind polynomials P i (x) = U i (x). In order to achieve the condition exposed before we will choose the same u γ (x) as for the first kind polynomials (1.35). The expansion combination giving the derivatives expansion can be written as:

{λ α,γ } = [∆ α,γ 2 ]{λ 0,γ-α } (1.42)
The constants in the matrix [∆ α,γ 2 ] are detailed for α = 1..4 in 4.2.3.2. The computation of the derivative operators (

[D],[∆ α,γ 1 ] and [∆ α,γ 2 ]
) are computed in first instance, before the identification process and whatever the signal looks like. Therefore this identification method is not time consuming for a given model and a given behaviour differential equation.

Introduction

In this chapter different practical issues are developed. Firstly some simulation tools are presented : the way a realistic noise is added to simulation signals, and how the noise contribution is tested, through a probabilistic approach. Then some numerical considerations on the identification process are discussed. They give the limits of the considered method for all the three steps. A focus on a regularization step is finally made.

Simulation tools 2.1.1 Noise

During all simulation tests, noise has been added to the analytical signal. Two main noise kinds have been considered :

-Additive noise, which is uncorrelated with the signal. Experimentally, this kind of noise is a part of a signal recorded with sensors. -Multiplicative noise, which is correlated with the signal. This kind of noise is observed when amplifiers are used for acquisition. Noisy signal is computed as follows:

f noisy = f exact ∆f m e j∆φ + ∆f a (2.1)
∆f m is a Gaussian real number with a mean value equal to unity and a standard deviation equal to a chosen percentage of the magnitude of the signal, ∆φ is another Gaussian random real number with a null mean value and a standard deviation of 1 deg and ∆f a is a Gaussian real number with a null mean value and a standard deviation equal to a chosen magnitude percentage of the signal. A common noise level observed when using piezoelectric sensors is approximately equal to 5%. Therefore this noise level will be often chosen for the simulations.

Monte Carlo test

The different methods proposed here are tested using a Monte Carlo test. The Monte Carlo test procedure consists in comparing the identification results with random signals generated with (2.1) and a given level of noise. A test criterion is chosen to facilitate this comparison. This test criterion is the identification error, which is computed as follows:

error(dB) = log 10 | θ I -θ ref θ ref | (2.2)
θ ref being the exact value of the identified parameters, and θ I the identified parameters. The outcome of the test is determined by the rank of the test criterion of the observed data relative to the test criteria of the random samples forming the reference set. For all following tests, a set of 100 samples will be observed.

Results presentation

Depending on the application, the results of the identification process will be mainly represented in two manners :

-Error versus truncation number : the error in (dB) is plotted versus the truncation number. An error equal to 10 -1 corresponds to an error of 10%. The number of sub-figures in the result presentation depends on the number of θ i parameters identified. An example is given in figure 2.1. In noisy case, the error for each Monte Carlo test is represented by a star. Hence, we can see on the same graph the mean error made during the identification process and the dispersion of the results. -Error versus truncation number and wave number : the error in (dB) is associated to a color map. A bad estimation (error greater than 100%) is associated to red color. A good estimation (error smaller than 1%) is associated to blue color.

An error between 1% and 100% is depicted by colors between orange and green. The X axis of the figure always represents the truncation number and the Y axis of the figure always represents the wave number. For a given wave number and truncation number a cell is filled with the color associated to the mean error in this identification configuration, as shown in the example 2.2.

Numerical considerations and regularization

We have seen in the previous chapter that the identification method was a three steps process. During each step, the noise can have a crippling effect. Hence we need to apply each step focusing on the following considerations.

STEP 1 : expansion considerations 2.2.1.1 Choice of the observation window

It has been shown in the previous chapter that the Chebyshev basis is made of symmetric polynomials (for order i even) and anti-symmetric polynomials (for order i odd) over the orthogonality domain Γ = [-1, 1]. If we choose an observation window, where the signal is close to a symmetric or anti-symmetric shape, it will induce some specificity in the expansion on the chosen basis. Let illustrate this consideration with two examples, an arbitrary chosen symmetric signal f s :

f s (x) = x 2 + x sin( π 2 x) (2.3)
and an arbitrary chosen anti-symmetric signal f a : For the symmetrical case in figure 2.3 we can see that all odd expansion coefficients are close to zero. For the anti-symmetrical case in figure 2. 4 we can see that all odd expansion coefficients are close to zero. This could be explained and generalized using the integral formulation of the expansion. Let focus on these null expansion coefficients. Due to the signal shape in the chosen window, these null coefficients are uniformly distributed over the vector containing the expansion coefficients. These null coefficients add no interesting information in our identification computation. Furthermore, if we need to compute the pseudo inverse of such a vector, its sparsity could induce numerical errors. Hence, a focus on the choice of the observation window should be made before any application of the identification principle presented in this work. The expansion on a polynomial basis P is defined as follows:

f a (x) = x 3 + sin(πx) (2.4)
f = λP λ i =< f, P i > (2.5)
This expansion is the best least square estimate of the function f on a basis of size N .

Compressed sensing approach

Compressed sensing is a statistical approach based on the sparsity of the vector containing the expansion coefficients. It considers that most of the expansion coefficients are null. Finding the sparest expansion on a basis of polynomials P can be written as :

min||λ|| 1 f = λP (2.6)
with ||λ|| 1 = i |λ i |. This is a convex optimization problem which can be solved by linear or quadratic programming techniques. The necessary conditions on P do not need the complete orthogonality, but recent results have proved that the more orthogonal the columns of P will be, the better the results of the estimation will be. Many solvers are available. The algorithms developed by Van den Berg and Friedlander ( [START_REF] Van | Probing the pareto frontier for basis pursuit solutions[END_REF]) have been tested. This algorithm is named the Basis pursuit denoising (BPDN). BPDN fits the least-squares problem approximately. Indeed, in the presence of noisy or imperfect data, it is undesirable to exactly fit the linear system f = λP . The constraint in (2.5) is relaxed to obtain the basis pursuit denoising (BPDN) problem:

min||λ|| 1 f -λP 2 < σ (2.7)
where the positive parameter σ is an estimate of the noise level in the data. However, the results of this method is the sparest expansion on a given orthogonal basis. The sparest expansion is not automatically the best least square approximation on a orthogonal basis of size N . The computation of the signal derivatives and the proposed identification method is based on the integral formulation of the scalar product and therefore on the theory of the best least square approximation.

Although the compressed sensing approach proposed by [START_REF] Simard | Non linear response identification: a compressed sensing approach[END_REF] has a very efficient filtering effect, it can not be applied, because it is based on two different theories. Indeed, as shown in (2.5) and (2.7) both expansion definitions are different. Compressed sensing was tested on a simple example but does not give accurate results. Figure 2.5 presents the different expansion tools : Gauss points, trapezoidal rule and compressed sensing approach (CS ). All expansion coefficients are computed with a constant number of sensors (here 18). The error made on the integral computation with only 18 signal values is very large. As expected the expansion computed with the trapezoidal rule is not accurate in this case. The aim of this figure is to compare both expansion theory (2.5) and 2.7. (2.5) can be well estimated using Gauss points. (2.7) is illustrated be the compressed sensing approach. In figure 2.5, we can clearly see a difference between the expansion obtained with (2.5 (Gauss points) and the one obtained with (2.7) (compressed sensing approach). Indeed, in figure 2.5(a) and (b), we can see that the smaller the expansion coefficients are, the greater the error is. Indeed, as the compressed sensing theory is based on the fact that only a few λ i are non-null, the small expansion coefficients are badly estimated, even without noise. Moreover, on figure 2.5(c) we can see clearly that the best least square approximation of the function F (green obtained with the Gauss points) is not identical to the one obtained with the compressed sensing approach (cyan). As in our work we need the best least square approximation of F , the compressed sensing approach is not an acceptable expansion tool for our application.

Truncation order

A straightforward filtering solution has been considered. Indeed the expansion on a orthogonal basis has itself a filtering effect. High oscillating perturbation such as noise are expanded on high order coefficients. The expansion of the signal will be concentrated on low order coefficients. Considering the expansion of the signal on a basis of size N :

f (x) ≈ N i=0 λ i P i (x) (2.8)
Using the criteria β :

if |λ i | < β max{|λ i |} N then λ i = 0 (2.9)
The criteria must be set by the user. A part of the signal information can be lost during this truncation. This lost information is not crippling in case of our identification method as we need the values of the λ i coefficients until i = N .

STEP 2 : Differentiation considerations

Differentiations using the [D] operator and using the novel method ([∆] operator) are compared. The harmonic response of a beam (for more information on the beam model, see the next chapter) is differentiated. The expansion of the fourth derivative calculated analytically and expanded with 1000 Gauss points is considered as a reference. The results are given for the order i < 6 (for i > 6, the expansion coefficients are 10 3 smaller than for the first 6). The grey lines represent the dispersion on the expansion error (performed with 1000 runs)

In figure 2.6, the error on computation of the fourth derivative with a truncated expansion (N = 27 and N = 28) is presented. These two truncation orders are selected for the signal studied (second modal response). Indeed, the 28 th expansion coefficient is smaller than the first 27 coefficients. This basis extension clearly reveals the error due to the [D] operator. For the noise free-case, considerable improvement (-10dB) is achieved with the novel differentiation method. The error due to the operator [D] appears on high order coefficients (> 3). With noise, the error is amplified, as shown in figure 2.6. Indeed, both methods give similar results for a sufficiently truncated basis (N = 27). However, the method based on the operator [D] is very sensitive to the truncation order N . For N = 28, the high order coefficients are corrupted by the bias and the error is greater than 1dB. This error induces extra oscillations on the computed derivative, as shown in Fig. 2.7. With the novel differentiation method ([∆]), the error on high order coefficients is significantly reduced. It still adds extra-oscillations but their amplitude is reduced compared to the [D] operator method. Therefore truncation order N does not require precise tuning. The identification procedure proposed before consists in transforming the equation of motion in algebraic equations. This set of algebraic equations is simplified in this chapter as follows:

ψ = θφ (2.10)
with θ a vector made of the parameters to be identified, φ and ψ are matrices made of the expansion of the signal and its differentiation. In case of a well posed problem (uniqueness and stability of the solution), the sizes of the matrices must be equal to the number of unknowns in θ. However, when the signal is affected by perturbations (such as noise), it becomes interesting to have more equations (size of matrices φ and ψ) than unknowns (θ). Indeed, the more information is used to solve the problem, the less the solution will depend on noise perturbations. The least square estimation consists in minimising the residual r by using the euclidean norm:

r = ψ -θφ (2.11)
Therefore we compute the pseudo inverse of φ and obtain the least square estimate of θ:

θ LS = (φ T φ) -1 φ T ψ (2.12)
But when φ and ψ are computed using a noisy signal, (2.10) is not true any more. In the following chapter the pseudo inverse (φ T φ) -1 φ T is named φ † . Instead, we have:

ψ * = θφ * + η (2.13)
φ * and φ * are matrices made of the expansion of the noisy signal and its differentiation, and η is a matrix corresponding to the effect of the noise. Because of η, the least square estimate of θ is biased. As our identification process is sensitive to noise perturbations, it is crucial to implement a regularization step. Three regularization method are proposed here : the weighted least square estimate, the single value decomposition and the instrumental variable.

Weighted least square

The weigthed least square method can be considered if we know the variance/covariance matrix Ω. We consider the Cholesky decomposition of Ω : P T P = Ω -1 and we multiply each member by P T in order to obtain:

ψ * = θφ * + η * (2.14)
with φ * = P T φ, ψ * = P T ψ and η * = P T η. Using this transformation, this model satisfy all hypothesis needed for of a classical least square model:

θ W LS = (φ T Ω -1 φ) -1 φ T Ω -1 ψ (2.15)
with θ W LS the weighted least square estimate of θ. This method was tested on some simple examples and was not efficient enough in our identification case. Therefore this method was not studied further.

Truncated single value decomposition

The identification method proposed here induces the computation of the pseudoinverse of the matrix ψ. When the matrix ψ suffers from bad conditioning, small errors (due for example to noise) made on the matrix components will be amplified by the inversion. The aim of the truncated single value decomposition (TSVD) is to eliminate these small errors before the inversion and to reduce the sensitivity of the solution on these small variations. This regularization method is based on the following assumption assumption : the sensitivity is due to the fact that a set of possible solutions exists, as the problem is overdetermined. For this regularization method we introduce the single value decomposition :

φ mn = U mn S nn V T nn (2.16)
with U T mn U mn = I nn , S nn a diagonal matrix, with the singular value in ascending order, and V nn a unitary matrix. With this formulation the pseudo-inverse of ψ can be written as:

φ † mn = V mn S -1 nn U T nn (2.17)
by choosing a well balanced solution. Its easy to understand that small variation in the singular values in S nn induces big variations in ψ † mn . Therefore the TSVD method consists in erasing the small singular values. The number of non-zeros values of S nn defines the rank of ψ and the level of linear dependence between its columns. In practice, with noise, all singular values are non-null but some of them are negligible in comparison to the others. The negligible singular values are however responsible of big errors during the inversion. Therefore we rewrite the expression of (2.17) as :

φ † mn = V mr S -1 rr U T nr (2.18)
In order to use this regularization method, it is compulsory to adjust a criteria in order to separate the negligible values to the others. In all application presented here this criteria has been very hard to adjust because all singular values were of the same order. Furthermore, the decomposition on the Chebyshev basis has a similar effect, as the noise in the signal will be expressed essentially by the high order polynomials.

Instrumental variable

The least square estimate (2.12) is very sensitive to small changes, which is the case in presence of noise. This regularization is based on the instrumental variable proposed by [START_REF] Young | Recursive Estimation and Time-Series Analysis: An Introduction for the Student and Practitioner[END_REF] and extended for instance by [START_REF] Schorsch | Instrumental variable methods for identifying partial differential equation models of distributed parameter systems[END_REF]. The main idea is to filter the signal by the model. The purpose of this method is to decorrelate from noise the most sensitive matrix (the pseudo-inverted array) by an iterative process. Let consider φ * = φ + φ , with φ being the contribution of noise. Let introduce the instrument φ ⊥ . The least square problem (2.12) can be written as :

φ T ⊥ φ * θ = φ T ⊥ (φ * θ + η) (2.19)
If we choose φ ⊥ to be independent of noise input, we have:

lim N →∞ φ T ⊥ φ = 0 lim N →∞ φ T ⊥ η = 0 lim N →∞ φ T ⊥ φ * = R φ ⊥ φ with R φ ⊥ φ = lim N →∞ φ T ⊥ φ (2.20)
with N the size of the vector φ,ψ (see 2.11)... Therefore (2. 19) can be transformed into:

φ T ⊥ φ * ( θ -θ) = φ T ⊥ ψ (2.21)
Therefore if R φ ⊥ φ is non singular, we have :

lim N →∞ ( θ -θ) = 0 (2.22)
The recursive estimates are asymptotically unbiased. We can rewrite (2.12) as :

θ IV = (φ T ⊥ φ) -1 φ T ⊥ ψ (2.23)
the subscript IV states here for instrumental variable. The simplest way to obtain φ ⊥ is summarized in the following and in figure 2.8:

1. θ is computed using the least square estimate (see (2.12)) θ = θ LS (the subscript LS states here for least square), 2. A solution f ⊥ of the equation governing the system behaviour with θ = θ LS is estimated (named after the solution of the auxiliary model),

3. Using this solution and the differentiation tools developed previously, the array φ ⊥ is computed (same as φ but using f ⊥ instead of f ), 4. θ IV is estimated using (2.23), 5. back to step 1 using θ = θ IV until convergence of the solution.

More elaborate manners to obtain φ ⊥ exist but are not treated as this one gives good results. φ ⊥ will be as much correlated to φ as the boundary conditions (of the auxiliary model) used to obtain the solution will be close to the real boundary conditions. The influence of the boundary conditions is discussed in the next chapter.

Ψ ϕ θ LS θ IV θ IV ϕ ψ ┴ 1. 2.
3.

4.

5.

Until convergence…

OR

Conclusion

In this chapter, we have proposed a three step identification method, made of an expansion step, a differentiation step and a parameter estimation step. Each step was briefly discussed considering the signal shape, the noise perturbation etc. Three regularization processes are proposed. The instrumental variable will be the unique regularization process applied in the following chapters.

Introduction

In this chapter different numerical applications will be depicted. Each application highlight some specifications of the method such as the choice of the truncation order, the wave number, the reformulation for the multiple parameter identification, etc.

The first example is a very basic problem where only one parameter is to be identified. This application permits a good analysis of the link between the signal contents and the choice of the tuning parameter : the truncation order. With the second example, the identification problem is reformulated for a multiple parameters identification, in order to satisfy the identifiability criteria. This example shows also the best regularization strategy in case of a multiple parameters identification. The structure damping estimation is depicted through two examples : using the transient response/ or using the steady state response of two structures. The last numerical application shows the performances of the process on a 2D structure. In the case of a generic plate, the choice of the instrument is particularly discussed. In this section, the forced response of a cantilever beam is considered as simple model. This model has been chosen because:

1. The frequency range of the response is easy to select, tuning the frequency of the excitation, 2. The response is easy to compute, 3. The equation of motion involves a single parameter to identify, 4. The derivative order is sufficiently high in order to study the method sensitivity.

The equation of flexural motion is:

∂ 4 ∂x 4 v(x, t) = ρS EI ∂ 2 ∂t 2 v(x, t) (3.1)
with v the transverse displacement, x the space variable, t the time variable, ρ the density, S the cross-section surface, E the Young Modulus and I the flexural inertia. The general solution of this equation is in the following form:

v(x, t) = sin(ωt) • (A sin(kx) + B cos(kx) + C sinh(kx) + D cosh(kx)) (3.2)
with k 4 = ω 2 ρS EI . In case of a cantilever beam, we can write the four boundary conditions:

v(0, t) = 0 ∂ ∂x v(0, t) = 0 ∂ 2 ∂x 2 v(L, t) = 0 EI ∂ 3 ∂x 3 v(L, t) = F (3.3)
if a concentrated force in the form F sin(ωt) is applied at the length x = L. The constants A, B, C and D in (3.2) can be computed using (3.22) as follows:

    0 1 0 1 1 0 1 0 -sin(kL) -cos(kL) sinh(kL) cosh(kL) -cos(kL) sin(kL) cosh(kL) sinh(kL)         A B C D     =      0 0 0 F k 3      (3.4) S 0 L F sin(ωt)

Identification application

The equation (3.1) is multiplied by the weighting function u(x, t) = (1-t 2 ) γt (1-x 2 ) γx :

u(x, t) ∂ 4 ∂x 4 v(x, t) = ρS EI u(x, t) ∂ 2 ∂t 2 v(x, t) (3.5)
with γ t > 2 and γ x > 4. In this example we choose γ t = 3 and γ x = 5. Therefore the equation of motion is transformed into the following algebraic equation:

{λ 4,γx|0,γt } = ρS EI {λ 0,γx|2,γt } (3.6)
{λ 4,γx|0,γt } being the fourth space partial derivative expansion and {λ 0,γx|2,γt } being the second time partial derivative expansion of v • u.

The single parameter to identify is θ = ρS EI .

Results

Noise free case

The single parameter (θ = ρS EI ) is identified using N × N equations, with N being the truncation order. Indeed for this example we did an expansion in the space direction and in the time direction. The expansion vector is therefore an expansion matrix of size N × N . This identification depends on two parameters :

-the truncation order N : which is linked to the number of equations used for the identification. The signal is expanded on a basis of size N . Therefore the results of this expansion give N coefficients. -the signal richness, which will be linked here to the parameter k, the wavenumber.

The wavenumber is a non-dimensional parameter which depicts a ratio between the frequency and structure size. The signal richness is linked to the number of non negligible expansion coefficients.

In figure 3.2, the error on the estimation of θ is plotted. The red color corresponds to a bad estimation (error between 100% and 10 11 %). The blue color corresponds to a good estimation (error between 0.0001% and 10 -11 %). For orange, green, light blue colors the error range is between 100% and 0.001%. The vertical axis is the wave number while the horizontal axis is the truncation order N . In this zone the number of polynomials chosen is not sufficient for a good signal expansion (as illustrated in figure 3.3). The signal is too rich in comparison to the polynomial basis chosen. Therefore, in this zone, the identification will never give a good result. For future identification cases a special attention should be paid on the signal expansion.

The lower right zone :

where k is relatively small and the truncation order N higher. In this zone the number of polynomials chosen is sufficient for a good signal expansion. Therefore, in this zone, the identification will give good results. In our noise free case, the parameter θ is reconstructed with an error lower than 10 -10 %. ). This is explained by the approximation maid with this operator. The higher the truncation order is, bigger the coefficients in the matrix [D] are. Multiplying the high coefficients by small expansion coefficients enhances the error made on the differentiation of the signal. This phenomenon is well illustrated in the chapter 2.

Noisy case

For noisy case, we have added 20% of noise and tested the identification method through a Monte Carlo test. The median error on 100 tests is plotted here. 1. The upper left zone : where k is relatively high and the truncation order N small. In this zone, the identification will never give a good result, as explained previously.

In case of the operator [D], the results seem to be good in this zone. However in the graphs presented in figure 3.5 an information is missing : the dispersion of the identified parameter. In this zone the identification dispersion is very high with the operator [D]; as illustrated in figure 3.6 for a wave number equal to 5. For N > 20 with the operator [D], the estimation of θ is highly dependent on the truncation number. The dispersion of the results make the identification process not reliable enough. Therefore, even if the mean error is low, the probability to obtain such an error is low too. As we have seen with the noise free case, the identified parameter is overestimated due to the high coefficients in the [D] matrix for big truncation order. In contrary, in the upper left zone, the identified parameter is underestimated due to the lack of expansion coefficients. In the noisy case, these two phenomena compensate each other, which explain these apparently good results.

The upper right zone :

where k is relatively high and the truncation order N higher. In this zone the number of polynomials chosen is sufficient for a good signal expansion. Therefore, in this zone, the identification will give good results. In our noise free case, the parameter θ is reconstructed with an error lower than 10 -5 % with our novel differentiation method. In comparison to the noise-free case, this zone is smaller and exclude the "bottom zone",

The bottom zone :

where k is relatively small and the truncation order N high. In this zone the signal is clearly not rich enough for the identification. The expansion of the signal is only expressed on few coefficients. Therefore, a major part of the expansion is corrupted by noise. The filtering effect of the expansion is lost when choosing a big truncation number. In this part the identification process is clearly biased. In the noisy case, the three differentiation methods gives opposite results. Indeed, where the novel operator gives good identification precision in the lower right zone, the operator [D] gives poor results. The higher the truncation order is, the worse the identification is. The error made during the differentiation process is enhanced by the noise, because the large coefficients for this operator multiplied by the small errors made on expansion coefficients induce large errors.

A 1D case: Influence testing of boundary condi-

tions for the regularization step 3.2.1 Theory

Model

we would like to test the limits of the regularization step. For the testing of the influence of the boundary conditions, the experimental model used is the cantilever beam as for 3.1. The data considered as "experimental data" is the solution of a whole cantilever beam (of length L). For the auxilary model we chose a sub-part of the beam used as the experimental model. The data used for the auxiliary model is the solution of a part of a cantilever beam (of length L) as shown in fig. 3.7. This part is L a = αL long and centred on the beam of length L. The regularization process is tested for α = 0.5..1. The auxiliary model will be a beam with following boundary conditions :

v a (0, t) = v( L -L a 2 , t) ∂ ∂x v(0, t) = ∂ ∂x v( L -L a 2 , t) v a (L a , t) = v( L -L a 2 + L a , t) ∂ ∂x v a (L a , t) = ∂ ∂x v( L -L a 2 + L a , t) (3.7)
In order to have a similar wave number for the auxiliary and the experimental model, we choose the frequency (for the auxiliary model) as follow:

σ = 2πL √ θ 1 f exp = 2πL a √ θ 1 f aux (3.8)
with σ the wave number, f exp the excitation frequency for the experimentation, f aux the frequency considered for the auxiliary model.

Identification application

In order to obtain φ ⊥ (see 2.19) as correlated as possible to the ψ (computed using experimental data), the boundary conditions considered for the auxiliary model must be as close as possible to the experimental boundary conditions. Unfortunately, this seems hard to obtain experimentally. Indeed, the experimental boundary conditions are not perfect (for example it is impossible to obtain experimentally a real clamped boundary condition). These experimental boundary conditions are also often influenced by unmeasurable, environmental excitations (such as wind, etc).

The parameter θ = ρS EI is identified.

Results

In the identification process developed previously, the truncation order must be chosen. In this section, we choose N = 50 in order to reconstruct precisely the recorded signal v(x, t). With the [∆] operator, figure 3.8 (left) shows that even if the boundary conditions are different the bias effect is always reduced and the IV estimate is always more accurate than the least square estimate. The best estimate is obtained for La = L. Indeed, where the boundary conditions are similar, the instrument ψ ⊥ is maximally correlated with φ. For some values of α the regularization, step seems to give bad results. If we look closer 3.9 that the obtained signal is close to the problematic shape considerations presented in the chapter 2. Indeed for example for α = 0.708 and α = 0.924 the signal is null at the boundaries of the window. The expansion of such a signal is difficult (because all Chebyshev functions are not equal to zero at boundaries). Therefore such boundary conditions for the auxiliary model will never be chosen for the regularization process. The Timoshenko theory considers two effects neglected in the Bernoulli approach:

-the shear effect, introducing the form factor K (K = 5/6 for a beam with rectangular section) and G the shear modulus, -and the rotational inertia effect, adding a fourth partial derivative in time. The equation of motion can be written as:

EI ∂ 4 ∂x 4 v(x, t) -ρI(1 + E KG ) ∂ 4 ∂t 2 ∂x 2 v(x, t) + ρ 2 I KG ∂ 4 ∂t 4 v(x, t) + ρS ∂ 2 ∂t 2 v(x, t) = 0 (3.9)
As explained in [START_REF] Brancheriau | Expertise mécanique des sciages par analyses des vibrations dans le domaine acoustique[END_REF], a solution of this equation is in the form:

v(x, t) = sin(ωt) • (A cos(k 1 x) + B sin(k 1 x) + C cosh(k 2 x) + D sinh(k 2 x)) (3.10) with k 1 = 1 L X(α + Γ) 2 + X 2 (α -Γ) 2 4 + X k 2 = 1 L - X(α + Γ) 2 + X 2 (α -Γ) 2 4 + X Γ = EI KSGL 2 α = I SL 2 X = ω 2 L 4 Sρ
EI In order to obtain this solution, the structure must satisfy the following conditions:

-the ratio between beam length and thickness is relatively small ( L h < 20)

-the mode square roots are similar to the Bernoulli solution (X ≈ k n L)

Identification application

We choose for this example to estimate the parameters of this partial differential equation (3.9) using the steady-state response of the beam. This choice is motivated by two reasons :

-firstly, the steady-state response is often the easiest measurable response (with a laser vibrometer for example) -secondly, this choice will induce some changes in the identification method. This point is explained in the following section. The steady state response is not time dependent. Therefore the equation of motion can be reduced to a simple differential equation. Multiplying this equation by u(x) we have:

u(x) • d 4 dx 4 v(x) + ρ E (1 + E KG )ω 2 u(x) • d 2 v dx 2 (x) + ( ρ 2 KEG ω 4 - ρS EI ω 2 )u(x) • v(x) = 0 (3.11) with u(x) = (1 -x 2 ) γ et γ > 4.
Expanding this equation we will obtain:

{λ 4,γ } = - ρ E (1 + E KG )ω 2 {λ 2,γ } -( ρ 2 KEG ω 4 - ρS EI ω 2 ){λ 0,γ } (3.12)
In this identification example we have 3 parameters to identify :

θ 1 = EI ρS , θ 2 = I S (1 + E KG
) and θ 3 = ρI KGS . In case of the steady state response of the beam, v, its second and its fourth derivatives are proportional. It is therefore impossible to identify the three parameters. Therefore we need to choose a set of reponses at minimum three different frequencies. We will rewrite the identification problem as follows:

  θ 1 θ 2 θ 3   T   λ 4,γ (ω 1 ) • • • λ 4,γ (ω H ) ω 2 1 λ 2,γ (ω 1 ) • • • ω 2 H λ 2,γ (ω H ) ω 4 1 λ 0,γ (ω 1 ) • • • ω 4 H λ 0,γ (ω H )   =   ω 2 1 λ 0,γ (ω 1 ) • • • ω 2 H λ 0,γ (ω H )   T (3.13)
The choice of the set size H and the different frequencies ω i will be discussed here after. The link between the pulsation ω i and the wave number k i expressed in the following section is:

k i = √ ω i 4 ρS EI (3.14)

Results

Noise free case

In case of the Timoshenko beam, there are three parameters to identify. Therefore the frequency set must be chosen keeping in mind these considerations :

-in order to satisfy the identifiability criteria, the size of the frequency set (H) must be equal to or greater than the number of identified parameters, -the different ω i must be far enough to each other in order to bring sufficiently different informations. Therefore in our test case we will choose different structure modes. This is not the only possible choice but gives interesting results, -with the example of the Bernoulli beam, we have clearly found a link between the signal richness and the optimal truncation number. We expect therefore in this section to observe the same phenomenon, linked to the highest wave number. The identification process is applied on a given frequency range, for the wave numbers equal to k = [1,2,3,4]. For the following results (presented in figure 3.10), The identification results depend only on the truncation order N . Therefore the following graphs present the error on the vertical axis and the truncation number on the horizontal axis. In figure 3.10, we can see clearly that for all three parameters the novel operator gives better results (with a precision of 10 -5 %) than the [D] operator (with a precision of 0.1%). Moreover, for the operator [D] the error evolves with the truncation order N . For small N , the error decreases until an optimal truncation order (this part is similar for the novel operator). After the optimal truncation order N , the error increases, whereas with the novel operator the error decreases until it becomes stable. The explanation of this error of high truncation order with the [D] operator was given previously. For the novel operator, the error decreases until all information given by the structure responses is expanded. After a given truncation order (N = 25), the additional expansion coefficients computed add no valuable information for the identification but does not affect the maximum precision (for N = 25). The precision obtained for θ 1 is greater than for θ 2 , which is greater than for θ 3 . This is explained by the different rough order of magnitude of the three parameters. Indeed where θ 1 magnitude order is about10 2 , θ 2 magnitude order is 10 -6 and θ 3 magnitude orders is 10 -11 . A major attention must be paid on these different order of magnitude. More than conditioning or computing difficulties, it will bring an additional challenge, as the contributions associated to θ 2 and θ 3 will be close to noise contribution magnitude order.

Noisy case

In this noisy simulation, we have added 5% of noise as presented for k = 7 in figure 3.11. The noise effect was tested through a Monte Carlo test with 100 samples. In figure 3.12 the mean error and the dispersion on the results are plotted.

Least square estimate

For the noisy case, the best precision is obtained on θ 1 for a truncation order greater than 15, with the novel operator (around 10%). However the variability in the obtained results is very high for N greater than 20. This is due to a minimal description of the signal. The estimates obtained with the [D] operator are unusable. Indeed the effect of the noise contribution multiplied by the high coefficients in [D] can not be compensated by the filtering effect of the expansion. The estimates of θ 2 and θ 3 suffers from the noise. Indeed the effect of these small contributions in the whole equation of motion are in the order of magnitude of the added noise. Therefore the estimate of these parameters can be made only with a regularization step, which will be presented in the next section.

Instrumental variable

The instrumental variable is an iterative process. We could develop a convergence criteria (convergence of the estimated parameter), but for this section we have chosen to do arbitrarily 10 iterations. The auxiliary model is the same than the model used for the simulation computation. For the novel operator, we have drastically enhanced the identification results for θ 1 and even more for θ 2 . For θ 2 we have a precision of 10%. This error is stable with the truncation order. For the parameter θ 3 , the identification results are still unusable. This could be explained by the arbitrarily chosen iteration number. Indeed, for better results, a convergence criteria on θ 3 could be computed. Furthermore, we can easily see two convergence lines in the three figures 3.13. For θ 2 the main part of the identification results converges to an error around 10% whereas some results converges to an error around 1000%. These aberrant results are well known from instrumental variable users. Indeed, the instrumental variable must be used with a stability check (for more information see [START_REF] Young | Recursive Estimation and Time-Series Analysis: An Introduction for the Student and Practitioner[END_REF]). In our case, no stability check has been used, which explains these aberrant results. For the operator [D], the regularization process does not help for the identification. As we have seen previously, even without noise, an error is made using this operator. The regularization step does only suppress the bias but could not erase the error made using the [D] operator, and enhanced with noise.

3.4 A 1D case using the transient response : a damped bar

Theory

In case of bar, the longitudinal equation of motion is:

∂ 2 ∂t 2 v(x, t) = E ρ ∂ 2 ∂x 2 v(x, t) + η ∂ ∂t v(x, t) (3.15)
where ρ is the density, E the Young modulus of the material, and η the viscous damping of the beam. t is the time variable, x the space variable and v is the longitudinal displacement.

For this study the particular case of free response of the simply supported bar will be treated. We can write the four boundary conditions:

v(0, t) = 0 ∂ ∂x v(0, t) = 0 v(L, t) = 0 ∂ ∂x v(L, t) = 0 (

Identification application

If we use a unique couple (γ x , γ t ) for all partial derivatives computation, we can rewrite (3.15) as :

{λ 2,γt|0,γx } = E ρ {λ 0,γt|2,γx } + η{λ 1,γt|0,γx } (3.17)
with γ x > 2 and γ t > 2. The results presented in figure 3.16 are similar than those for the Bernoulli beam for θ 1 and for θ 2 . The figures present all same three zones, especially for the novel operator.

Results

Noise free case

The results obtained with the novel operator are better than for the [D] operator, especially in the most interesting zone (upper right). In figure 3.17, the regularization process permits a good identification in case of high truncation order N and small wave number k, the most noise sensitive zone, as we have seen previously. In the case of the operator [D], the regularization process does not enhance the identification precision.

Instrumental variable

A 1D case identifying the damping with the steady

state response : cantilever beam 3.5.1 Theory

Model

As explained in [START_REF] Woodhouse | Linear damping models for structural vibration[END_REF][START_REF] Adhikari | Identification of damping: part 1, viscous damping[END_REF][START_REF] Adhikari | Identification of damping: part 2, non-viscous damping[END_REF] different models for structural vibration exists. All are theoretical approximation of an observed phenomenon but have not a real physical meaning.

Here we have considered a viscous damping model, where the damping effect is represented by a term proportional to the transverse velocity of the beam (first order partial derivative in time):

∂ 4 ∂x 4 v(x, t) = ρS EI ∂ 2 ∂t 2 v(x, t) + η ∂ ∂t v(x, t) (3.18)
where η is the structural damping.

For the computation of the response of this structure, we have considered the Young Modulus E complex. The imaginary part of the Young modulus represents the damping coefficient. For this example, the steady-state response of the beam is considered. We have considered here a damping equal to 1% of the Young Modulus.

Identification application

Multiplying the equation of motion (3.18) by u(x, t) we have:

u(x, t) ∂ 4 ∂x 4 v(x, t) = ρS EI u(x, t) ∂ 2 ∂t 2 v(x, t) + ηu(x, t) ∂ ∂t v(x, t) (3.19) with u(x, t) = (1 -t 2 ) γt (1 -x 2 ) γx et γ t ≥ 2,γ x ≥ 4.
Expanding this equation we have:

{λ 0,γt|4,γx } = ρS EI {λ 2,γt|0,γx } + η{λ 1,γt|0,γx } (3.20)
As we have two unknows (η and ρS EI ) we need to rewrite (3.20) as:

  ρS EI c EI   T {λ 2,γt|0,γx }(ω 1 ) • • • {λ 2,γt|0,γx }(ω P ) {λ 1,γt|0,γx }(ω 1 ) • • • {λ 1,γt|0,γ X }(ω P ) =    {λ 0,γt|4,γx }(ω 1 ) . . . {λ 0,γt|4,γx }(ω P )    (3.21)
The choice of the set size H and the different ω i frequencies follows the same considerations as discussed for the Timoshenko beam.

Results

Noise free case

There are two parameters to identify : θ 1 = ρS EI and θ 2 = η. This identification process is applied on a given frequency range, for the wavenumbers equal to k = [2,3,4]. This influences only the optimal truncation number, as this will be discussed here after. In the noise-free case, as shown in figure 3.18, the identification is possible for a truncation number N greater than 25. This is directly linked to the frequency set chosen for the identification. The link between the signal richness (therefore the frequency set) and the truncation number was already discussed previously with the Bernoulli beam. For a good precision on the identified parameters, the signal should be correctly expanded. This condition is reached for N = N opt . After the optimal truncation order, the error is stable and equal to 10 -8 % for θ 1 and 10 -6 % for θ 2 , with the novel differentiation tool. With the original [D] operator, the minimal error is around 1% and increases after N opt .

Noisy case

Least square estimate

The least square estimate of [θ 1 , θ 2 ] is presented in figure 3.19. For the two parameters, for a small truncation order, the mean error is relatively small (10%) but the dispersion on the 100 Monte Carlo tests is high. For higher truncation orders,the mean error is no more acceptable. The optimal truncation number (N = 25 for the noise-free case) can not be reached because of the noise. In figure 3.20, the instrumental variable enhances the results for the novel differentiation method. For the [D] operator, the results remain unchanged. For the novel operator, the error decreases from 10% to 1% for θ 1 and is stable, equal to 10% for θ 2 . As for the Timoshenko beam the difference of magnitude order between θ 1 and θ 2 explains the difference in the minimum identification error (θ 2 ≈ θ 1 100 ).

3.6 a 1D case with application to discontinuity location 3.6.1 Theory

Model

In this section, a cracked beam is modelled by an healthy homogeneous beam with a section change at x, as shown in figure 3.21. The beam model is identical to the Bernoulli

v 1 (0, t) = 0 ∂ ∂x v 1 (0, t) = 0 v 1 (X - χ 2 , t) = v 2 (0, t) ∂ ∂x v 1 (X - χ 2 , t) = ∂ ∂x v 2 (0, t) v 2 (χ, t) = v 3 (0, t) ∂ ∂x v 2 (χ, t) = ∂ ∂x v 3 (0, t) ∂ 2 ∂x 2 v 3 (L -X + χ 2 , t) = 0 EI ∂ 3 ∂x 3 v 3 (L -X + χ 2 , t) = F (3.22)
with X the crack position, χ the crack length, v 1 being the transverse displacement of the beam named 1 in figure 3.21, v 2 being the transverse displacement of the beam named 2 in figure 3.21 and v 3 being the transverse displacement of the beam named 3 in figure 3.21.

Identification application

For this application we use the same identification principle as for the Bernoulli beam (see 3.1). We fix N = 16 (the same as for the experimental results presented in chapter 4). Using the results obtain in section 3.22, we choose k = 1. In this section, we would like to estimate the influence of the crack position on the method sensitivity, in order to set the experimentation in the next chapter accurately. Therefore a crack (10% height, 0.5% length of cross section) is moved numerically from one end to the other of the beam considered. The signal is expanded on the Chebyshev basis for each crack position. The computed expansion coefficients λ 4,γx|0,γt i of the cracked beam and the healthy beam are compared.

Results

In figure 3.22, the results are presented for i = 11 (the observations are similar for other i, see figure 3.23). The grey vertical lines correspond to the Gauss points of the Chebyshev polynom of order i.

In figure 3.22, the difference between the cracked and healthy expansion coefficient number 11 is clearly sensitive to the position of the crack regarding to the Gauss-points. We can see that near the Gauss points, the value of the expansion coefficients is doubled. When the crack is located close to the middle distance between two Gauss points, the value of the expansion coefficients is unchanged (ratio close to 1). This observation is confirmed at different expansion order, as shown in 3.23. Hence a crack located near the Gauss-points of the i th polynomial implies a noticeable change on the computed expansion coefficient of order i. Therefore the crack can be located depending on which expansion coefficients change. 3.7 A 2D case with a single identified parameter: a plate

Theory

In this 2D application we have partial a differential equation which involves partial differentiation of the signal in both directions and a cross partial differentiation. This aspect is a good challenge for our identification process.

Model

The equation of motion of the undamped Kirchhoff plate is:

∂ 4 ∂x 4 w(x, y, t) + 2 ∂ 4 ∂x 2 ∂y 2 w(x, y, t) + ∂ 4 ∂y 4 w(x, y, t) = - ρh D ∂ 2 ∂t 2 w(x, y, t) (3.23)
with

D = Eh 3 12(1 -ν 3 )
, and ν Poisson coefficient. The steady state response of the plate will be studied. Therefore the equation of motion can be reduced to:

∂ 4 ∂x 4 w(x, y) + 2 ∂ 4 ∂x 2 ∂y 2 w(x, y) + ∂ 4 ∂y 4 w(x, y) = ω 2 ρh D w(x, y) (3.24)
A simple solution of this partial differential equation can be computed in case of a simplysupported plate. The steady-state response w reconstructed by modal decomposition is:

w(x, y) = ∞ n=1 ∞ m=1 4F L x L y (ω 2 mn -ω 2 ) sin( mπX 0 L x )sin( nπY 0 L y )sin( mπx L x )sin( nπy L y ) (3.25)
with L x , L y being the lengths of a square plate, X 0 , Y 0 the coordinates of the excitation with a pulsation ω and a magnitude F (see figure 3.24). We have the relationship between the mode number and the pulsation : As expected, the solution (3.25) is true ∀(x, y) with (x, y) = (X 0 , Y 0 ) (where the excitation is applied). Therefore we will work on a sub-part of the plate. This sub-part of lengths L x(SP ) , L y(SP ) must exclude the point (X 0 , Y 0 ) (see figure 3.24).

ω mn = π 2 D ρh (( m L x ) 2 + ( n L y ) 2 ) (3.26)

Identification application

Multiplying (4.1) by the function u(x, y) we have :

u(x, y) ∂ 4 ∂x 4 w(x, y) + 2 ∂ 4 ∂x 2 ∂y 2 w(x, y) + ∂ 4 ∂y 4 w(x, y) = ω 2 ρh D u (x, y)w(x, y) (3.27)
with u(x, y) = (1 -y 2 ) γ (1 -x 2 ) γ et γ ≥ 4. Computing the expansion of this equation we will obtain : {λ 4,γ|0,γ } + 2{λ 2,γ|2,γ } + {λ 0,γ|4,γ } = ω 2 ρh D {λ 0,γ|0,γ } (3.28) Without noise, as we can see in figure 3.25, the results are similar for θ for the three differentiation methods. The best precision is around 10 -8 % and the figure presents the same two zones as for figure 3.2.

Results

Noise free case

Noisy case

In the noisy simulation, we have added 5% of noise in both direction (the recorded signal is matrix of samples. Each sample corresponds to a measurement at a given x and y position. Each sample is affected by noise). The noise effect was tested through a Monte Carlo test with 100 samples. Here the mean error on the results is plotted.

Least square estimate

For this application we have treated both cases :

-the identification of θ = D ρh , presented in figure 3.26,

ω 2 θ = {λ 0,γ|0,γ } † {λ 4,γ|0,γ } + 2{λ 2,γ|2,γ } + {λ 0,γ|4,γ } (3.29)
-The identification of θ * = ρh D , presented in figure 3.27.

θ * = ω 2 {λ 4,γ|0,γ } + 2{λ 2,γ|2,γ } + {λ 0,γ|4,γ } † {λ 0,γ|0,γ } (3.30)
The comparison of both identification strategy will give particularly interesting considerations for the instrumental variable. Lets focus now on the least square estimate of θ.

As shown in figure 3.26, the identification is totally biased. The identification results are relevant only for high wave number k, where the amount of information in the signal is sufficient and for a sufficiently high value of N where the signal information is correctly reconstructed. For higher truncation order, the identification results are totally corrupted with noise. Let highlight that the pseudo inverted matrix in this case is the expansion of the partial derivatives ({λ 4,γ|0,γ } + 2{λ 2,γ|2,γ } + {λ 0,γ|4,γ } in 4.3). The estimation of partial derivatives based on a noisy signal enhances the bias effect as it multiply small perturbations by the differentiation coefficients. The computation of the pseudo inverse enhances even more these effects.

- In figure 3.27, the identification results are similar to the one presented in previous numerical applications. We find the same three zones. The pseudo inverted matrix in this case is directly the expansion of the signal. Therefore after the optimal truncation order N the error increases proportionally to the truncation order (for low wave number the error goes from green (0.1%) to orange (100%), for high wave number the error goes from blue (0.01%) to green (0.1%)). The choice of this identification strategy (estimation of ρh D

) is of high interest. Indeed, even if the zone where the identification result is accurate is smaller than for 3.26, the precision is much greater for 3.27. Moreover, the error made in the left zone, due to high truncation number, and therefore the expansion of the noise, could be corrected by the regularization step. , for a simply supported plate, with 5% of noise. As shown in figure 3.28, the regularization step does not enhance the identification precision. The same uniform low precision is observed on the upper right zone. As expected, the error in this case is principally due to an incorrect use of the identification process (bad signal reconstruction with a small N , signal not sufficiently rich) and is not due to noise. Hence the regularization step can not enhance the identification results. As shown in figure 3.30, the precision is drastically enhanced on the right zone, as expected:

Instrumental variable

1. in the upper right zone, the best identification precision is around 0.001% in the upper right corner (until 0.01% without regularization). In this case the instrument is the signal expansion of the auxiliary model. The filtering effect of the model is highly efficient, as the error is principally due to the noise expansion.

2. in the lower zone, a zone (where the error is equal to 0.01%) still exists. Here the regularization process is less efficient, as the instrumental variable can not add data for the parameter estimation.

Conclusion

In this chapter, the different numerical applications show the different specificity of this identification process :

-with the noise-free cases, we have seen that an optimal truncation order exists. This optimal truncation order is obtained when the whole noise-free signal is reconstructed. -with noise and the least square estimate, we did a focus on the filtering effect of the expansion. It has been shown that without regularization, the method was really sensitive to the tuning parameter N (truncation order). Indeed it was very difficult to find a good balance between signal reconstruction and filtering. -the instrumental variable was highly efficient in order to overcome the previous difficulty. The signal was filtered by the model and this permits good identification results for high truncation order. -the identification process was adapted to multi-parameters identification. This novel identification process satisfies the identifiability criteria using a set of signals with different excitation frequencies. -the multi-parameters identification highlights the relevance in choosing the convergence criteria for the instrument variable. -for the damping identification, two kinds of signal can be used : transient and steady-state. For each kind a focus is made on the adaptation of the method. -with 2D structures, a focus is made on the choice of the instrument. Depending on the inverted matrix, the regularization step will filter the solution by the model or by the model and the differentiation method. The second choice seems to be more efficient in our identification process.

Introduction

In this chapter two experimentations are presented. The first one is an application for damage detection. The theory developed for the Bernouilli beam is applied on a sub-part of the experimental beam. The evolution of the parameter along the beam is observed. The second experimental application permits material characterization. The damping on a plate is estimated as an example.

Damage detection on a cantilever beam 4.1.1 Principle

For this experiment, the equation of the bending beam presented in the previous chapter is considered. The whole identification process used here is the one presented for the Bernoulli beam. The parameter θ = ρS EI is identified. θ being proportional to the beam cross-section area S, it would be interesting to see if a crack (and therefore a thickness change) could be reconstructed and identified. This identification method has two main applications. The first is linked to the material properties of the structure. It permits the identification of a global ρS/EI for the whole beam. This method has many applications, such as model updating for a controller and monitoring of Young's Modulus (material ageing). The second involves more geometrical properties of the structure. It permits the computation of the dispersion of ρS/EI along the beam for each sensor array.

Each sensor array has 16 sample positions. Therefore, the signal is expanded on the basis of 15 Chebyshev functions. In order to ensure a minimal wave number equal to 1, the length of the sensor array is adapted at each excitation frequency. For example, at the 5 th flexural mode of the beam, the sensor array length is equal to 1/5 th of the beam length. No specific hypothesis can be made regarding the crack location on the sensor array. Therefore the sensor array is shifted along the beam. If the crack was located in a dead zone for a given sensor array position, it would be located in a sensitive zone for one of the next sensor array positions , as shown in Fig. 4.1.

Experimental setup

The displacement of this beam is reconstructed with selected time samples and space positions using a laser vibrometer (PSV400). The beam is excited with a shaker at a single frequency by a sinusoidal source. The steady-state response is measured point by point. The whole structure response is reconstructed using a force sensor, which states as reference. Two cases are studied here : the wave number on a patch k ≈ 0.8 and k ≈ 1. As we have seen in the previous chapter, the wave number is directly linked to the identified parameter. Therefore, we can not consider this information in order to tune our experimentation. For this reason we will approximately choose the excitation frequency (in order to obtain k ≈ 0.8 and k ≈ 1) and use both results to show the sensitivity of the identification process to the choice of excitation frequency. Here, we study the forced response of this cantilever beam, 1.1m long in which a crack is imposed (a notch 3mm in width and 2.5mm in depth). For this experimentation, a relatively large crack is chosen as a first test of this method. Fig. 4.2 shows the beam, the imposed crack and laser measurements at different locations. A sensor array provides 16 measurements at 16 different sample positions (as shown in Fig. 4.2).

Results and discussion

Parameter ρS/EI and the fourth derivative (∂ 4 v/∂t 4 • u) can be evaluated for all the sensor array locations for undamaged (x < 0.8m and x > 0.8m)) and damaged cases The value of the mean identified parameter ρS/EI is computed and plotted in dashed lines: dark grey for k ≈ 1 and light grey for k ≈ 0.8. For both cases, the mean value is close to the theoretical one ( ratio equal to 1.16 for k = 0.8 and 0.85 for k = 1). The mean value can be corrupted by noise which could explain the dispersion of the measurements. The theoretical value is also roughly estimated with the beam dimensions and the properties of the steel. The dispersion of ρS/EI is plotted in figure 4. 3 and 4.4 (left). Each filled box corresponds to the parameter dispersion at a single sensor array position. The computed parameter values oscillate between 1/2 and 3 times the mean value, except when the sensor array is located near x = 0.8m. At this location, the computed parameter becomes negative or is greater than 4 times the mean value. The oscillating values (between 1/2 and 3 times the mean value) can be explained by the variability due to the noise. The assumption of a continuous structure is not accurate at the damage location, therefore the computed parameter can be negative. Conversely, a negative parameter shows a breach in the formulated hypothesis and therefore damage. The effect of damage (or a discontinuity) is similar to that of noise. Indeed, a discontinuity changes the slope of the signal slope locally. This local change drastically increases the high order coefficient of expansion. In the case of noise, the effect of these high order coefficients is smoothed by the computation of a mean value. In the case of damage, the effect of these high order coefficients cannot be smoothed by the same technique. Therefore the value of the computed parameter at damage location can increased by more than 4 times, due to these high order coefficients. This reasoning is confirmed in Fig. 4.3 and 4.4 (right). Indeed the calculation of the displacement's fourth derivative with respect to space emphasises these high order coefficients. Near the damage location, the computed fourth derivative consists of highoscillating terms. For both cases (k ≈ 0.8 and k ≈ 1), the maximum value of the fourth derivative is located at the Gauss-point closest to the damage (see chapter 3 for more details). These results demonstrate an alternative method for damage location. Indeed, with a restricted number of samples (here 16 samples), it was shown that ρS/EI and the fourth derivative ∂ 4 v/∂x 4 can be computed. The parameters studied are very sensitive to damage (as a discontinuity). ρS/EI becomes negative when the continuity assumption is no longer valid. The damage can also drastically increase the ρS/EI value. The method

Damped plate identification 4.2.1 Principle

We consider an isotropic thin plate of thickness h, driven by a harmonic force of angular frequency ω. The equation governing the forced vibrations is

D ρh (1 -jη)[ ∂ 4 ∂x 4 w(x, y) + 2 ∂ 4 ∂x 2 ∂y 2 w(x, y) + ∂ 4 ∂y 4 w(x, y)] = ω 2 w(x, y) (4.1)
where w is the transverse displacement, E the Young's modulus, η the loss factor, D the flexural rigidity and ρh the mass per unit area. This equation can be written as follow using the weighting function u(x, y):

u(x, y) D ρh (1 -jη) ∂ 4 ∂x 4 w(x, y) + 2 ∂ 4 ∂x 2 ∂y 2 w(x, y) + ∂ 4 ∂y 4 w(x, y) = ω 2 u(x, y)w(x, y) (4.2) with u(x, y) = (1 -y 2 ) γ (1 -x 2 ) γ et γ ≥ 4.
Computing the expansion of all partial derivatives of w, (4.1) can be transformed in a algebraic equation:

D ρh

(1 -jη){λ 4,γ|0,γ } + 2{λ 2,γ|2,γ } + {λ 0,γ|4,γ } = ω 2 {λ 0,γ|0,γ } (4.3)

Experimental setup

A polystyrene plate with dimensions L x = 0.5m, L y = 0.5m, thickness h = 3.5mm and density ρ = 1055kg.m -3 is used in this experiment. The plate was suspended at its upper corners to approximate free boundary conditions (see figure 4.5). It was excited by a shaker at its lower boundary. A PCB 288D01 impedance head was mounted on the shaker to provide input force measurement. A Polytec scanning vibrometer PSV 300 was used to measure the velocities over a 71 × 71 regular meshgrid, on an area free of excitation. This experimentation was conducted by F. Ablitzer, who needed a regular meshgrid in order to apply the RIFF method (see explanation just after). We chose to use in our case the approximate evaluation of the scalar product, using the trapezoidal rule, for more convenience. We have therefore to focus our attention on border effects which are linked to this approximation. Morever, if we had to set this experimentation, especially using the vibrometer, we would choose meshgrid at Gauss points with a significant lower density. This would reduce the acquisition time drastically.

Results and discussion

The identified results are compared to the RIFF method developed by [12]. The RIFF (Resolution Inverse Filtree Fenetree in French, meaning windowed filtered inverse resolution) technique applied on inverse problem was led by C. Pezerat [8,9,10,11]. This method considers the equation of motion of the structure. the partial derivatives of the structure response are estimated by a finite scheme. Using two previous consideration, the distribution forces in the structure are reconstructed. Since measured data are always noisy, the force distribution reconstructed is regularized applying a spatial window on the set of sensors used for the partial derivatives reconstruction. This method was successfully applied by F. Ablitzer on this experimentation and will be compared with our results. In figure 4.6, for a frequency range of [0.05 -0.5kHz], the results seem hard to interpret. The high variation in the identified parameters is due to a small wave number in the considered area. Therefore, the signal is not rich enough in order to estimate the considered parameters D ρh and η. In this frequency range both methods seems to overcome some difficulties . For a frequency range of [0.5 -2kHz], with the RIFF method, the parameter D ρh increases slowly and the damping η exhibits higher variations. For a low frequency range of [0.5 -2kHz], the mean loss factor is around 0.07. With our identification method, the parameters seem hard to interpret as the results dispersion is very high. For a higher frequency range of [2 -6.4kHz], the RIFF method gives a mean loss factor around 0.4 which decreases slowly. With our identification method, the mean loss factor is also around 0.4 but seems to be constant in this frequency range. As the behaviour of the tested material is not clearly known, it is impossible to determine which of the two results is more accurate.

Least square estimate

Instrumental variable

In figure 4.6, our method overcomes some difficulties in the frequency range [0.5 -2kHz]. Therefore it would be interesting to filter the results in this frequency range by a regularization step. As shown previously, for our instrumental variable computation, we need an auxiliary model. Unfortunately, it is very hard to compute analytically the response of a free-free plate. Therefore we choose as auxiliary model a simply supported plate. The results of the identification with the regularization step are presented in figure 4.7. For the identification of the loss factor η, the influence of the noise is still very high.

As we have seen in all our simulation cases, the low wave number case is always very hard to exploit. Indeed, only a few expansion coefficients will be large enough for our identification process. As the estimate is based on only a few coefficients, it is very sensitive to noise. This noise influence can easily be overcame by the regularization step for D ρh . Unfortunately for the estimation of η, the amount of information at these frequencies is not sufficient with our method. In order to obtain optimal filtering by the instrumental variable, the boundary conditions of the plate could be adapted (taking a sub-part of the plate), in order to generate an auxiliary signal as close (in term of shape) to the experimental one.

Conclusion

In this chapter, we applied successfully the proposed identification method on two examples. The first experimental application deals with crack detection on a cantilever beam. The crack created on the experimental beam was located using the estimation of a beam parameter (proportional to material and geometrical properties) along the whole beam. The identification method was also tested on a free-free plate. The purpose of this experimentation was to evaluate the damping of a PMMA material.

Conclusion and perspectives Conclusion

The aim of this work was to create an identification method, using the structure response and able to reconstruct structure parameters. The method should not depend on the structure environment, which is considered as unknown. The method was based on a simplified structure model, which was the partial differential or differential equation governing the structure behaviour. For this identification purpose a novel differentiation model was proposed. The aim of this novel operator was to limit the sensitivity of the method to the tuning parameter (truncation number). The precision enhancement using this novel operator was highlighted through different examples. An interesting property of Chebyshev polynomials was also brought to the fore : the use of an exact discrete expansion with the polynomials Gauss points. The Gauss points permit an accurate identification using a restricted number of sensors, limiting de facto the signal acquisition duration. In order to reduce the noise sensitivity of the method, a regularization step was added. This regularization step, named the instrumental variable, was inspired from the automation domain. The instrumental variable works as a filter. The identified parameter is recursively filtered through the structure model. The final result is the optimal parameter estimation for a given model. Different numerical applications have been presented. Each application depicts a particular identification specification. The Euler-Bernoulli beam highlighted the link between the truncation number, the signal richness and the identification precision. The Timoshenko beam shows the input needs regarding to a multi-parameters identification. The examples of damped beam and bar show how the damping ratio of a structure could be estimated, using either the steady state or the transient response of the structure. The multi-parameters identification examples permitted the discussion on the choice of the convergence criteria for the regularization step. The two dimensional case of a plate emphasises the need to define the identification problem correctly, in order to optimize the regularization step. The choice of the instrument was discussed, as this choice permits the filtering of the method through the model and through the differentiation process. Two experimental applications were treated. For the beam structure, the crack detection problem was illustrated. The method was applied on a cantilever beam to locate imposed damages , although the experimental technique can be applied to any type of structure. Indeed, parameter ρS/EI and the fourth derivative ∂ 4 v/∂x 4 were computed, using a signal recorded experimentally. Parameter ρS/EI was accurately reconstructed for the whole beam, using different sensor patch lengths. The damage was accurately located using this identification method. The computed fourth derivative emphasised the discontinuity due to the crack. These changes in continuity could be computed with the novel differentiation technique, making it very interesting for damage location. The second experimental application was focused on the damping ratio estimation on a plate. It permits the computation of the inertial parameter D/ρh and the damping ratio η. The method was compared to the RIFF technique. Both methods give similar results. The proposed identification method was efficient even without regularization. To sum up, the proposed identification method shows good capabilities for the different aimed identification objectives : structure monitoring such as material ageing estimation, damage detection or model parameters estimation for control purposes.

Perspectives

This work could be enhanced, or its application domain widen through different perspectives. The most critical step of this method is currently the expansion step. This step could be easily enhanced by the design of specific sensors, which could directly reconstruct the scalar product (the integral reconstruction, see [START_REF] Chesne | Identification of plate boundary forces from measured displacements[END_REF]). For this work the estimation of the components in the [∆] matrix were computed order by order. A generalization of this computation would permit the estimation of higher differentiation matrix (for a differentiation order greater or equal to 5). Maybe other polynomials could have different interesting properties and could be used through the same process (such a Taylor series or Gegenbauer polynomials). A large panel of acquisition methods could be tested with the continuous identification approach, such as field measurements, using fast camera for example. Another interesting perspective, which was not explored yet, could be the estimation of the boundary conditions. Indeed it has been shown that the closer the boundary condition to the real boundary conditions in the regularization step are, the more the instrument is correlated to the experimental data. Hence the regularization step could be used as a boundary condition estimator. Furthermore, it has been shown in the automation domain that more powerful instruments exist. The theory associated to these more elaborate instrument was however not tested yet. D. Remond developed in [14] different tools for non-linear systems. These tools could be easily adapted to multi-dimensional cases and permits the identification on non-linear systems. A first test of the operator [D] on a non linear structure was conducted by F. Martel. Unfortunately, during his work the novel operator [∆] was not available yet. The application of the novel operator on this working example could be a first step to non linear structure identification. These tools permit for example the identification of non homogeneous material properties or could be used to estimate a structure shape function, in order to estimate a crack depth. Indeed instead of considering a linear partial differential equation of motion, we could choose such equation for example :

ρS(x) ∂ 2 v ∂t 2 = ∂ 2 ∂x 2 (EI(x) ∂ 2 v ∂x 2 ) (4.4)
which is the equation of motion of a non homogeneous beam.

The identification techniques could be extended as a grey box method. The parameters of a general partial differential equation could be estimated. The non relevant parameter would be close to zero or at any case negligible. This could be used to determine the behaviour of novel and unknown materials. The appearance of novel terms could be interpreted as a significant change in the structure, such as a plastic deformation, etc. A back up to continuous theory is more and more used for the identification and modelling of mechanical structures. Because mechanical structures are continuous systems, the continuous time and space identification should be the most relevant identification techniques, even if it involves a discrete sensing approach.

INTRODUCTION

Numerous applications require parameter estimation of continuous time systems. This important and broadly reviewed field [1,2] is often used on discrete time data issued from experimental systems or numerical simulations. Linear time-varying (LTV) systems are getting more attention due to more sophisticated systems and models. Identification of the time-varying parameters of LTV systems can be achieved via various approaches, such as ensemble average of impulses responses [3], parallel-cascade algorithm [4], or wavelet-based methods [5].

An alternative approach is to identify time-varying parameters by estimating the coefficients of their projection on an orthogonal basis [6,7]. Some commonly used basis include Legendre series [8], block-pulse functions [9], Fourier series [10] and Laguerre polynomials [11]. Chebyshev polynomials were also used for time-varying parameter estimation [12][START_REF] Garnier | Identification of continuous-time models from sampled data[END_REF][14] and were reported to have certain advantages over other orthogonal series [15]. However, methods employing Chebyshev polynomials were, to our knowledge, formulated via integrals and thus require initial conditions to be known or identified simultaneously.

This paper introduces a parameter estimation method for LTV systems that is based on Chebyshev polynomials and formulated via signal derivatives instead of integrals. Amongst other properties, Chebyshev basis include a matrix representation of its coefficients and a linear relationship for both derivation and multiplication operations [16]. The differential equation of the system can thus be linearized and the unknown coefficients can be estimated by means of standard least-square algorithms. The time-varying parameters are then reconstructed from the time expansion of the coefficients on the orthogonal basis.

Experiments were conducted to investigate the ability of the method to identify a continuous mechanical system with timevarying stiffness, constant inertia and negligible damping. The signal filtering process was investigated in details in order to deal with the multi vibration modes of the continuous system and, in particular, the potential cross-coupling between torsional and bending modes as well as the system excitation by a moving load on the bar.

*Manuscript

Click here to view linked References

INTRODUCTION

Numerous applications require parameter estimation of continuous time systems. This important and broadly reviewed field [1,2] is often used on discrete time data issued from experimental systems or numerical simulations. Linear time-varying (LTV) systems are getting more attention due to more sophisticated systems and models. Identification of the time-varying parameters of LTV systems can be achieved via various approaches, such as ensemble average of impulses responses [3], parallel-cascade algorithm [4], or wavelet-based methods [5].

An alternative approach is to identify time-varying parameters by estimating the coefficients of their projection on an orthogonal basis [6,7]. Some commonly used basis include Legendre series [8], block-pulse functions [9], Fourier series [10] and Laguerre polynomials [11]. Chebyshev polynomials were also used for time-varying parameter estimation [12][START_REF] Garnier | Identification of continuous-time models from sampled data[END_REF][14] and were reported to have certain advantages over other orthogonal series [15]. However, methods employing Chebyshev polynomials were, to our knowledge, formulated via integrals and thus require initial conditions to be known or identified simultaneously.

This paper introduces a parameter estimation method for LTV systems that is based on Chebyshev polynomials and formulated via signal derivatives instead of integrals. Amongst other properties, Chebyshev basis include a matrix representation of its coefficients and a linear relationship for both derivation and multiplication operations [16]. The differential equation of the system can thus be linearized and the unknown coefficients can be estimated by means of standard least-square algorithms. The time-varying parameters are then reconstructed from the time expansion of the coefficients on the orthogonal basis.

Experiments were conducted to investigate the ability of the method to identify a continuous mechanical system with timevarying stiffness, constant inertia and negligible damping. The signal filtering process was investigated in details in order to deal with the multi vibration modes of the continuous system and, in particular, the potential cross-coupling between torsional and bending modes as well as the system excitation by a moving load on the bar.

*Manuscript

Click here to view linked References 2 2. METHODS

Chebyshev polynomials properties

The n order Chebyshev polynomial is defined in the interval t = [-1, +1] by the following equation:

( ) ( ) (
).

arccos cos t n τ T n =
Any function x(t) can be expanded on the n order Chebyshev basis

( ) { } ( ) ( ) ( ) ( ) T n n τ T τ T τ T τ T τ T L 2 1 0 = as: ( ) ( ) { } ( ) ( ) { }, , 1 0 t t t t n C n n i T x x T x x x x x × = × = L (1)
where C x contains the x i coordinates of x(t) in the basis {T n (t)}.

If t is sampled on N e discrete points of index m, (1) becomes:

( ) ( ) [ ] , , 1 1 e e N n n n C N m T x m x + + × =
and the projection coefficients can be estimated using a least-square approximation:

( ) ( ) [ ] ( ) [ ] ( ) [ ] ( ) . 1 - × × × = T n n T n C m T m T m T m x x
Since the derivative of each Chebyshev polynomial can be approximated by a linear combination of lower order polynomials, the time derivative of a function can be expressed as:

[ ] { }, d d 2 1 0 n n T D x x x x t x × × = L
where [D] is the derivative operator described in [17]. This operator gives increasing importance to higher order polynomials.

If two polynomials P(τ) and Q(τ) are expressed in the n order Chebyshev basis as:

( ) ( ) ( ) ( ), τ T q τ Q , τ T p τ P n k k k n k k k å å = = = = 0 0
then their product can also be expressed in that basis:

( ) ( ) ( ), τ T a τ Q τ P n i i i å = = × 2 0 where ( ) ï ï ï î ï ï ï í ì > £ £ + + = + = å å å å - = - = + + = - = n. if i q p n, i if q p q p q p , if i q p q p a n n i j j i j n -i j j i j j j i j j j n j j j i 2 1 1 2 1 2 1 0 2 1 0 1 0 1 1 0 0 (2)
This linear combination of the p k and q k coefficients can also be formulated with a combination matrix [A]:

[ ],

A p a q × = whose coefficients are function of áqñ and are computed according to (2). For identification purposes, this matrix can be truncated so that its dimensions are consistent with coefficient vectors of appropriate length. This approximation will be used and detailed in the following section.

Application to the identification of 1-DOF second order mechanical system with time-varying parameters

The dynamic behavior of a 1-DOF mechanical system having continuous time-varying parameters can be represented by the following equation: 4

Calibration of the experimental apparatus

Calibration of the plastic bar stiffness was conducted by taking quasistatic measurements of the torque as a function of disk rotation for different linear positions of the slider. The stiffness values were estimated from a linear regression model on the torque-rotation curves and the coefficient of determination was computed to validate the linearity of the torsional stiffness for each slider position.

Since the section of the bar was chosen as an equilateral triangle, simple torsion analytical equations can be obtained for validation purposes [19,20]. According to [21], the torque M t required to twist the bar of an angle θ is:

θ 45 3 4 l G h M t =
where h is the height of the triangle (7.8 mm in the current situation), G is the shear modulus of the material and l is the length of the bar. Thus, the torsional stiffness k is defined by:

l G h M k t 45 3 θ d d 4 = = . ( 5
)
As a validation purpose, we computed G from the experimental data with a least square estimate using

( ) 1 4 1 1 1 3 45 - × × × = T T l l l k h G ,
and we compared the value obtained with the material properties datasheet of the urethane rubber used.

The inertia of the disk was estimated by conducting 10 tests of 10 s without the plastic bar. The motor was operated in torque mode and programmed to exert a sine sweep torque from 1.8 Hz to 3.1 Hz. The recorded torque and disk rotation signals were further low pass filtered at 100 Hz (3 rd -order Butterworth) in the forward and reverse time directions to achieve zero phase shift. The data corresponding to the filter start-up transients were excluded until the step response of the filter remained within the ± 5 % error range of its final value. The time basis was further divided in 0.5 s segments with 80 % overlap. The abovementioned identification method, with n = 60, was used to identify the inertia on a zero order Chebyshev basis (n ID_I = 0) for each segment. The values obtained for the inertia from each segment were averaged to estimate the disk inertia.

Identification of the time-varying stiffness

Identification of the bar time-varying stiffness was conducted by exerting a sine sweep torque from 1.8 to 3.1 Hz with increasing amplitude. A positive torque offset was maintained throughout the tests to avoid entering the dead zone. Data sampled at 1 kHz were then low pass filtered at a cutoff frequency of 6 Hz (4 th -order Butterworth) in the forward and reverse time directions to achieve zero phase shift. The data corresponding to the filter start-up transients were excluded until the step response of the filter remained within the ± 5 % error range of its final value. Data was further divided in 1 s segments with 80% overlap. Assuming orders n ID_I = 0, n ID_K = 20 and n = 60, the abovementioned method was used to identify the timevariable stiffness for two cases: one where the inertia is known and one where the inertia is unknown. The resulting stiffness estimates were then compared with the analytical solution provided in (5). As an estimate of the identification quality, an error power ratio (EPR) was computed. It served as an indicator of the correspondence between two signals sampled on the same time basis:

( ) ( ) ( ) ( ) ( ) [ ], % 100 EPR 1 2 1 1 2 2 1 å å = = - = e e N i N m m u m u m u (6)
where u 2 is the signal compared to the reference signal u 1 , N e is the number of samples, and m is the sample index.

Selection of the low pass filter cutoff frequency for stiffness identification

The previous low pass filter cutoff frequencies were chosen at 6 Hz since the input signal was limited to 3.1 Hz maximum. In order to validate whether such level does in fact lead to the identification of the first torsional mode of the bar, a modal analysis of the bar was conducted. According to [22], the natural torsional frequencies of a disk of inertia I d attached to one end of a rod of inertia I r , fixed at the other end, are given by: For the system under study, these equations indicate a first torsional mode under 6 Hz and a second mode beyond 923 Hz. The large separation of these two modes is such that a cutoff frequency at 6 Hz would apparently lead, without any doubt, to the identification of the first vibration torsional mode, which was the purpose of this investigation.

RESULTS AND DISCUSSION

System Calibration

The bar rotational stiffness was measured for various slider positions and results are illustrated in Fig. 2. As expected, the stiffness decreased as the inverse of the bar effective length (defined as the distance between the disk and the slider). The coefficient of determination R² was computed for each position of the slider and its mean over all those values was found to be 0.9820. The value of G obtained by the linear regression was 0.773 GPa while the value given in the material properties datasheet was 0.771 GPa. The identified inertia was found to be 1.039e -3 kg.m² with a standard deviation of 0.606 %. As a comparison, the inertia computed by a computer-aided design (CAD) software was 1.015e-3 kg.m². Such small difference can be explained by the presence of bolts and a certain portion of the torque cell that were not taken into account by the CAD software model.

Bar torsional stiffness identification

Once torque and position data were low pass filtered at 6 Hz, data were used to identify the variation in bar rotational stiffness during an experimental test. Results are shown in Fig. 3 for the case when the inertia is already known. Data shows a practically good fit and the EPR on the theoretical stiffness was found to be only 6.62e -2 %. If the inertia is not known a priori, the resulting stiffness identification for the same data set is shown in Fig. 4, and the inertia of the system was found to be 1.95e -3 kg.m² (89% EPR). The EPR for the theoretical stiffness was found to be about 1.44 % this time, higher than when the inertia was not identified in the process. As shown in Fig. 4, predicted stiffness was particularly diverging at the end of the testing sequence. Moreover, system identification using low pass cutoff filter frequencies of 10 Hz or 15 Hz led to mediocre identification.

These problems can be understood based on a modal analysis of the bar. Indeed, although the bar was excited in torsion by the brushless motor, bending vibration modes could be excited as well. As a matter of fact, by moving the slider along the vertical bar, the rollers could induce bending in the beam if there is a nonzero transverse force exerted by the slider on the beam. Such moving load phenomenon is well known in the literature [23]. In practice, we did observe bending of the free hanging bar when the slider was moving along the bar. A transverse force at the slider may get induced, for instance, if the geometrical center of the three slider's rollers is not perfectly coaxial with the axis of the direct drive brushless motor shaft.

Bending modes of the bar can be estimated analytically, as a first approximation, by computing the bending natural frequencies of the two bar segments: a first segment located between the disk and the slider with fixed/pinned support boundary conditions, and a second segment beyond the slider, using fixed/free boundary conditions.

According to [22], the beam natural frequencies ω n in flexion for a uniform beam of length l are defined by: , ρ ω where E is the beam material Young's modulus, I A is the second moment of area, ρ is the density and A is the area of the crosssection. The first two values of α corresponding to the first two modes are [3.92660 7.06858] for the fixed-pinned beam and [1.87510 4.69409] for the fixed-free beam. Fig. 5 illustrates these natural frequencies as a function of the effective length, assuming that the section of the bar is perfectly equilateral. The lowest first mode (Fixed-Free -Mode 1) is for the segment beyond the slider and its natural frequency is close but still higher than 6 Hz. The shaded area represents the span covered by the slider during the experiment. A finite element model of the bar, using 3D elements, was created to figure out the natural frequencies of the bar as a whole, assuming fixed/pin/free boundary conditions. Results are also presented in Fig. 5. Again, the first bending mode starts close to 6 Hz but is still larger for all effective lengths of the bar. In practice, however, the bar does not have a perfect equilateral shape such that the moment of inertia of the bar may slightly vary about one axis or the other. Hence, there would be several bending natural frequencies depending on the axis of bending.

Typical signals from the conducted experiment are shown in Fig. 6. In the course of one test, the slider moved pretty much all over the complete span of the bar. The angular position of the disk is shown in addition to the torque cell data that exhibits significant "noise". A power spectrum of the torque signal is shown in Fig. 7. The torque input signal has a high amplitude under 5 Hz as one would expect, but a spectral peak also exists starting at about 10 Hz, its amplitude reducing progressively up to 20 Hz. This peak is most likely attributed to the bar bending modes. Peaks are also present at about 30 Hz, 60 Hz and 90 Hz, which are most likely bending modes of the bar as well. Hence, it appears that the torque cell signal gets contaminated beyond 6 6 Hz. The low pass cutoff frequency selected at 6 Hz significantly reduces the contamination and this explains why the identifications with cutoff at 10 Hz or 15 Hz did not lead to a good identification quality.

The question remains, however, how bending modes can contaminate the torque cell signal. It is known that when the shear center of a beam section does not coincide with its centroid, bending modes can excite torsion modes [24]. In our case, given that the bar profile is probably not perfectly equilateral, this cross-coupling condition can occur. In addition, if two bending modes at different natural frequencies occur about different axes of the bar at the same time, it is most likely that torsion would occur in the beam. This can be shown when filtering the torque data with a band-pass filter at 30 Hz (Chebyshev, Type II, 14 th order), resulting in the torque signal shown in Fig. 8. One can observe a beating phenomena that is most likely due to the simultaneous occurrence of two bending modes (with slightly different natural frequencies) about two different axis of the bar.

Order selection and segment lengths

The identification method introduced in this paper requires a prior selection of the order of the polynomials and the determination of data segment lengths (in seconds) to be used for the identification. This question is fundamental for many system identification methods. In our identification method, we used n = 60 and segments lengths of one second. The impact of alternative values is shown in Fig. 9 that illustrates the EPR value in % on the stiffness parameter. The graphs shows a low sensitivity in the EPR value for higher orders and lower segments lengths but, as expected, increases exponentially for smaller orders and longer segments.

CONCLUSIONS

In this paper, we introduced a method that uses the Chebyshev orthogonal basis to identify time-varying parameters of a cantilever beam excited in torsion. Results show that a mechanical investigation of the vibrating modes of the system, albeit first order, provides a clear process to determine the best filtering procedure to obtain good parameter estimations. Once proper filtering procedures are determined, the system stiffness was best estimated when the inertial characteristics of the system were known a priori. Further improvement could be made to the technique by optimizing the orders and segment lengths based on convergence determination algorithms. One should remember, however, that identification methods of LTV systems can hardly be evaluated on validation data sets after calibration because many systems' response may vary from one test to another.

FIGURE CAPTIONS

Fig. 1. Time-varying stiffness 1-DOF system. The torque exerted by the motor on the disk makes it rotate by an angle θ and is measured by the torque cell. A torsion bar is connected to the disk and constrain its rotation. A lead-screw linear actuator controls the slider position such that the torsion bar support opposite to the torque cell can be moved along the bar at will. Fig. 2. Bar torsional stiffness as a function of its effective length, from both theoretical estimates given by (5) and experimental data. Let consider the scalar product between f (2) • u γ and T n :

λ 2,γ n =< f • u, T n > = δ 2,γ 1 (n -2, n)< f • u γ-2 , T n-2 >... + δ 2,γ 1 (n, n)< f • u γ-2 , T n >... + δ 2,γ 1 (n + 2, n)< f • u γ-2 , T n+2 > (9)
for n ≥ 2, with

δ 2,γ 1 (n -2, n) = ( 1 2 -γ + n/2)(1 -γ + n 2 ) δ 2,γ 1 (n, n) = (2γ -1)(γ -2) - n 2 2 δ 2,γ 1 (n + 2, n) = ( 1 2 -γ -n/2)(1 -γ - n 2 )
For n = 0 we have:

λ 2,γ 0 =< f • u, T 0 > = δ 2,γ 1 (0, 0)< f • u γ-2 , T 0 > + δ 2,γ 1 (2, 0)< f • u γ-2 , T 2 > (10)
with δ 2,γ 1 (0, 0) = 2(γ -

1 2 )(γ -2) δ 2,γ 1 (2, 0) = 2(γ - 1 2
)(γ -1)

For n = 1 we have:

λ 2,γ 1 =< f • u, T 1 > = δ 2,γ 1 (1, 1)< f • u γ-2 , T 1 > + δ 2,γ 1 (3, 1)< f • u γ-2 , T 3 > (11)
with δ 2,γ 1 (1, 1) = (γ -

3 2 )(3γ -2) -γ δ 2,γ 1 (3, 1) = γ(γ - 1 2 )
We can then construct the matrix [∆ 2,γ 1 ] made of the δ 2,γ 1 (i, j) coefficients:

[∆ 2,γ 1 ] =                      d 0,0 0 d 2,0 0 • • • 0 0 0 0 0 d 1,1 0 d 3,1
0 0 0 0 • • • d N -4,N -2 0 d N -2,N -2 0 d N,N -2 0 0 0 0 • • • 0 d N -3,N -1 0 d N -1,N -1 0 0 0 0 • • • 0 0 d N -2,N 0 d N,N                      (12
) with d i,j = δ 2,γ 1 (i, j). Therefore we have:

{λ 2,γ } = [∆ 2,γ 1 ]{λ 0,γ-2 } ( 13
)
is exact if λ 0,γ-2 N +1 = 0 and λ 0,γ-2 N +2 = 0 .

(

) 25
where d are non-null values, with d i,j = δ 3,γ 1 (i, j), i being the column indices, and j the array indices. Therefore we have:

{λ 4,γ } = [∆ 4,γ 1 ]{λ 0,γ-4 } (26)
is exact if λ 0,γ-4 N +1 = 0, λ 0,γ-4 N +2 = 0 and λ 0,γ-4 N +3 = 0.

.2 Second kind .2.1 First order derivative (α = 1)

For n ≥ 1:

λ 1,γ n =< f • u γ , U n > = δ 1,γ 2 (n + 1, n) < f • u γ-1 , U n+1 > +δ 1,γ 2 (n -1, n) < f • u γ-1 , U n-1 > = δ 1,γ 2 (n + 1, n)λ 0,γ-1 n+1 + δ 1,γ 2 (n -1, n)λ 0,γ-1 n-1 (27)
0 0 0 • • • 0 0 δ 1,γ 2 (N -1, N ) 0                 
(29) Therefore we have:

{λ 1,γ } = [∆ 1,γ 2 ]{λ 0,γ-1 } (30)
is exact if λ 0,γ-1 N +1 = 0 .

.2.2 Second order derivative (α = 2)

For n ≥ 2: (n -2, n -1)

λ 2,γ n =< f • u γ , U n > = δ 2,γ 2 (n + 2, n)λ 0,γ-2 n+2 + δ 2,γ 2 (n, n)λ 0,γ-2 n + δ 2,γ 2 (n -2, n)λ 0,γ-2 n-2 (31)
δ 2,γ 2 (n -2, n) = δ 1,γ 2 (n -1, n)δ 1,γ-1
For n = 0, we have: For n = 1 we have:

λ 2,γ 0 =< f • u γ , U 0 > = δ 2,
λ 3,γ 1 =< f (3) • u, U 1 > = δ 3,γ 2 (4, 1)λ 0,γ-3      (46)
where d are non-null values, with d i,j = δ 3,γ 2 (i, j), i being the column indices, and j the array indices. Therefore we have: {λ 4,γ } = [∆ 4,γ 2 ]{λ 0,γ-4 } (47) is exact if λ 0,γ-4 N +1 = 0, λ 0,γ-4 N +2 = 0 and λ 0,γ-4 N +3 = 0.

Figure 1 -

 1 Figure 1 -Estimation de D ρh , pour une plaque, avec 5% de bruit.

Figure 2 -Figure 3 -Figure 3 ,

 233 Figure 2 -ρS/EI ID /ρS/EI T H (a) et (∂ 4 v/∂x 4 • u) (b) pour la poutre encastrée libre expérimentale

  .4 . . . . . . . . . . . 1.4 Distance between P p[u] and dP p[u]/dt (a), and P p[u] and d 2 P p[u]/dt 2 (b), versus the size of the Chebychev basis from [16] . . . . . . . . . . . . . . 1.5 f and its truncated expansion (a), f (1) its truncated estimation and its estimation with [D] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 Estimation of θ with [∆ 1 ] for k = 5 for the cantilever beam (Bernouilli) without noise (see chapter 3) . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Estimation of θ with [∆ 1 ] for the cantilever beam (Bernouilli) without noise (see chapter 3) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 f s and its expansion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 f a and it expansion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Difference between the different expansion approach. (a) expansion coefficients values, (b) error on the expansion coefficient |λ i -λ i (theo)| |λ i (theo) | (c) expanded signal (obtained thanks to P. Simard [17]) . . . . . . . . . . . . 2.6 Error (dB) on the fourth-order differentiation coefficient for the beam second mode. (a): noise-free case, (b): with 15% noise. The results are given for the order i < 6 (for i > 6, the expansion coefficients are 10 3 smaller than for the first 6). The grey lines represent the dispersion on the expansion error (performed with 1000 runs) . . . . . . . . . . . . . . . . 2.7 displacement v • u (a) and fourth derivative v (4) • u (b) for the cantilever beam. Light grey: signal computed analytically; light pink: computed with the novel differentiation method; dark purple: computed with the [D] operator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.8 Summary of the identification and regularization process . . . . . . . . . 3.1 Model of the cantilever beam . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Estimation of θ for the cantilever beam (Bernoulli) without noise . . . . xiii 3.3 Illustration of expansion error for small truncation number (here k = 4 and N = 8), the time origin is on the front . . . . . . . . . . . . . . . . . 3.4 Illustration of the error made with the [D] operator (here k = 5.5 and N = 45) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Estimation of θ for the cantilever beam (Bernoulli) with 20% of noise . . 3.6 Estimation of θ for the cantilever beam (Bernoulli) with 20% of noise, for k = 5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.7 Signal recorded on the whole beam (green) and estimated displacement on the auxiliary model (pink), with similar wave number σ for both models .

Figure 4 :

 4 Figure 4: Presentation of the identification problem

Figure 1 . 1 :

 11 Figure 1.1: Chebyshev polynomials of the first kind for i = 0..4

Figure 1 . 2 :

 12 Figure 1.2: Gauss points of the first kind polynomials T 6 (x) and T 7 (x)

Figure 1 . 3 :

 13 Figure 1.3: Chebyshev polynomials of the second kind for i = 0..4

Figure 1 . 4 :

 14 Figure 1.4: Distance between P p[u] and dP p[u]/dt (a), and P p[u] and d 2 P p[u]/dt 2 (b), versus the size of the Chebychev basis from[16]

1

 1 Simulation tools . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 2 . 1 :

 21 Figure 2.1: Estimation of θ with [∆ 1 ] for k = 5 for the cantilever beam (Bernouilli) without noise (see chapter 3)

Figure 2 . 2 :

 22 Figure 2.2: Estimation of θ with [∆ 1 ] for the cantilever beam (Bernouilli) without noise (see chapter 3)

Figure 2 . 3 :

 23 Figure 2.3: f s and its expansion

Figure 2 . 4 :

 24 Figure 2.4: f a and it expansion

Figure 2 . 5 :

 25 Figure 2.5: Difference between the different expansion approach. (a) expansion coefficients values, (b) error on the expansion coefficient |λ i -λ i (theo)| |λ i (theo) | (c) expanded signal (obtained thanks to P. Simard [17])

  operator [D] (27poly.) novel method (27poly.) operator [D] (29poly.) novel method (29poly.)

Figure 2 . 6 :

 26 Figure2.6: Error (dB) on the fourth-order differentiation coefficient for the beam second mode. (a): noise-free case, (b): with 15% noise. The results are given for the order i < 6 (for i > 6, the expansion coefficients are 10 3 smaller than for the first 6). The grey lines represent the dispersion on the expansion error (performed with 1000 runs)

Figure 2 . 7 :

 27 Figure 2.7: displacement v • u (a) and fourth derivative v (4) • u (b) for the cantilever beam. Light grey: signal computed analytically; light pink: computed with the novel differentiation method; dark purple: computed with the [D] operator

Figure 2 . 8 :

 28 Figure 2.8: Summary of the identification and regularization process

  simple 1D case with a single identified parameter : the Bernoulli beam . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.1 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 A 1D case: Influence testing of boundary conditions for the regularization step . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 A 1D case with multiple identified parameters : the Timoshenko beam . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 A 1D case using the transient response : a damped bar . . 3.4.1 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 A 1D case identifying the damping with the steady state response : cantilever beam . . . . . . . . . . . . . . . . . . . . 3.5.1 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 a 1D case with application to discontinuity location . . . . 3.6.1 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.7 A 2D case with a single identified parameter: a plate . . . 3.7.1 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.7.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3. 1 A

 1 simple 1D case with a single identified parameter : the Bernoulli beam 3.1.1 Theory 3.1.1.1 Model

Figure 3 . 1 :

 31 Figure 3.1: Model of the cantilever beam

Figure 3 . 2 : 1 .

 321 Figure 3.2: Estimation of θ for the cantilever beam (Bernoulli) without noise

Figure 3 . 3 :

 33 Figure 3.3: Illustration of expansion error for small truncation number (here k = 4 and N = 8), the time origin is on the front

Figure 3 . 4 :

 34 Figure 3.4: Illustration of the error made with the [D] operator (here k = 5.5 and N = 45)

Figure 3 . 5 :

 35 Figure 3.5: Estimation of θ for the cantilever beam (Bernoulli) with 20% of noise

Figure 3 . 6 :

 36 Figure 3.6: Estimation of θ for the cantilever beam (Bernoulli) with 20% of noise, for k = 5

Figure 3 . 7 :

 37 Figure 3.7: Signal recorded on the whole beam (green) and estimated displacement on the auxiliary model (pink), with similar wave number σ for both models

Figure 3 . 8 :

 38 Figure 3.8: Influence of the regularization step, light : results obtained in case of least square estimate, dark : estimation with the instrumental variable

Figure 3 . 9 :

 39 Figure 3.9: Shape of auxiliary model response for problematic values of α

Figure 3 . 10 :

 310 Figure 3.10: Estimation of θ 1 ,θ 2 and θ 3 for the cantilever Timoshenko beam, without noise

Figure 3 . 11 :

 311 Figure 3.11: Noisy response of the Timoshenko beam

Figure 3 . 12 :

 312 Figure 3.12: Estimation of θ 1 , θ 2 and θ 3 for the cantilever Timoshenko beam, with 5% of noise. The error for 100 Monte Carlo tests are plotted

Figure 3 . 13 :

 313 Figure 3.13: Estimation of θ 1 , θ 2 and θ 3 for the cantilever Timoshenko beam, with 5% of noise, with the instrumental variable. The mean error on 100 Monte Carlo tests is plotted.

Figure 3 . 14 :

 314 Figure 3.14: Transient damped bar response (for k = 5)

Figure 3 . 15 :

 315 Figure 3.15: Estimation of θ 1 et θ 2 for the damped bar, without noise

Figure 3 . 16 :

 316 Figure 3.16: Estimation of θ 1 et θ 2 for the damped bar, with 5% noise. The mean error on 100 Monte Carlo tests is plotted.

Figure 3 . 17 :

 317 Figure 3.17: Estimation of θ 1 et θ 2 for the damped bar, with 5% noise, with the instrumental variable. The mean error on 100 Monte Carlo tests is plotted.

Figure 3 . 18 :

 318 Figure 3.18: Estimation of θ 1 and θ 2 with the steady-state of the damped beam, without noise

Figure 3 . 19 :Figure 3 . 20 :

 319320 Figure 3.19: Estimation of θ 1 and θ 2 with the steady-state of the damped beam, with 5% of noise

Figure 3 . 21 :

 321 Figure 3.21: Cracked beam model

Figure 3 . 22 :

 322 Figure 3.22: Computed expansion coefficients of the fourth derivative with respect to space. Ratio between cracked and healthy expansion coefficient, expansion order i = 11. The grey vertical lines correspond to the Gauss points of the Chebyshev polynom of order i.

Figure 3 . 23 :

 323 Figure 3.23: Computed expansion coefficients of the fourth derivative with respect to space. Comparison between cracked and healthy expansion. From (a) to (d): expansion order i = 2 and 6. The grey vertical lines correspond to the Gauss points of the Chebyshev polynom of order i.

Figure 3 . 24 :

 324 Figure 3.24: Model of the Kirchhoff plate

Figure 3 . 25 :

 325 Figure 3.25: Estimation of θ = D ρh , for a simply supported plate, noise free case.

Figure 3 . 26 :

 326 Figure 3.26: Estimation of D ρh , for a simply supported plate, with 5% of noise.

Figure 3 .

 3 Figure 3.27: Estimation of ρh D , for a simply supported plate, with 5% of noise.

Figure 3 . 28 :

 328 Figure 3.28: Estimation of D ρh , for a simply supported plate, with 5% of noise, with the regularization step.

Figure 3 . 29 :

 329 Figure 3.29: Estimation of ρh D , for a simply supported plate, with 5% of noise, with the regularization step.

Figure 3 . 30 :

 330 Figure 3.30: Gain (dB) on the estimation of ρh D , between least square and instrumental variable estimate, with 5% of noise

Figure 4 . 1 :

 41 Figure 4.1: Positioning of the sensors patches on the beam structure

Figure 4 . 2 :

 42 Figure 4.2: Experimental implementation: cantilever beam, imposed crack and laser measurements at different sample positions (for one sensor array)

Figure 4 .

 4 3 and 4.4 (left) shows the ratio between the identified and theoretical ρS/EI.

Figure 4 . 3 :

 43 Figure 4.3: ρS/EI ID /ρS/EI T H (a) and (∂ 4 v/∂x 4 •u) (b) for the cantilever beam, results based on a experimental data (excitation frequency equal to 1471Hz). On the left, the dashed lines correspond to the mean values for k ≈ 0.8

Figure 4 . 4 :

 44 Figure 4.4: ρS/EI ID /ρS/EI T H (a) and (∂ 4 v/∂x 4 •u) (b) for the cantilever beam, results based on a experimental data (excitation frequency equal to 1471Hz). On the left, the dashed lines correspond to the mean values for k ≈ 1

Figure 4 . 5 :

 45 Figure 4.5: Experimental device for the identification on the polystyrene plate

Fig. 4 .

 4 Fig. 4.6 shows the properties identified by the RIFF and our identification method in the range [0.05 -6.4kHz].

Figure 4 . 6 :

 46 Figure 4.6: Identification of D ρh and η at different frequencies, using a simple least square method (without regularization).

Figure 4 . 7 :

 47 Figure 4.7: Identification of D ρh and η at different frequencies, using the instrumental variable.

Fig. 3 .

 3 Fig. 3. Theoretical time-varying stiffness superimposed with the identified stiffness for each segment and the mean of those segments assuming a known inertia of the disk.

Fig. 4 .

 4 Fig. 4. Theoretical time-varying stiffness superimposed with the identified stiffness for each segment and the mean of those segments, assuming unknown inertia.

Fig. 5 .

 5 Fig. 5. Natural frequencies in torsion and flexion of the bar for both the two segment analytical models and a complete finite element model of the bar. Shaded area corresponds to the slider position range during a test.

Fig. 6 .

 6 Fig. 6. Typical effective length of the bar, torque and torsion angle readings during a test with the bar and sine sweep torque excitation.

Fig. 7 .

 7 Fig. 7. Power spectrum of the torque reading during a test with the bar and a sine sweep excitation. Hanning weighted windows of 2048 samples with 50 % overlap.

Fig. 8 .

 8 Fig. 8. Torque signal band-pass filtered between 25 and 35 Hz.

Fig. 9 .

 9 Fig. 9. Natural logarithm of the EPR between identified stiffness and predicted stiffness from slider position. Order n ID_K set to one third of n.

with δ 1 1 1

 11 -1, n) = -n + 1 + 2γ 2 For n = 0 we have: λ 1,γ 0 =< f • u γ , U 0 > = δ 1,γ 2 (1, 0)λ 0,γ-

with δ 2

 2 ,γ 2 (n + 2, n) = δ 1,γ 2 (n + 1, n)δ 1,γ-1 2 (n + 2, n + 1) δ 2,γ 2 (n, n) = -δ 1,γ 2 (n + 1, n)δ 1,γ-1 2 (n, n + 1) -δ 1,γ 2 (n -1, n)δ 1,γ-1 2

  ID /ρS/EI T H (a) et (∂ 4 v/∂x 4 •u) (b) pour la poutre encastrée libre expérimentale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii

	1	Estimation de	D ρh	, pour une plaque, avec 5% de bruit. . . . . . . . . . . . vii
	2 ρS/EI 3 Identification de	D ρh	et η à différentes fréquences, sans régularisation . . . viii

4 Presentation of the identification problem . . . . . . . . . . . . . . . . . . 1.1 Chebyshev polynomials of the first kind for i = 0..4 . . . . . . . . . . . . 1.2 Gauss points of the first kind polynomials T 6 (x) and T 7 (x) . . . . . . . . 1.3 Chebyshev polynomials of the second kind for i = 0.

  3.8 Influence of the regularization step, light : results obtained in case of least square estimate, dark : estimation with the instrumental variable . . . . 3.9 Shape of auxiliary model response for problematic values of α . . . . . . 3.10 Estimation of θ 1 ,θ 2 and θ 3 for the cantilever Timoshenko beam, without noise . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.11 Noisy response of the Timoshenko beam . . . . . . . . . . . . . . . . . . 3.12 Estimation of θ 1 , θ 2 and θ 3 for the cantilever Timoshenko beam, with 5% of noise. The error for 100 Monte Carlo tests are plotted . . . . . . . . . 3.13 Estimation of θ 1 , θ 2 and θ 3 for the cantilever Timoshenko beam, with 5% .21 Cracked beam model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.22 Computed expansion coefficients of the fourth derivative with respect to space. Ratio between cracked and healthy expansion coefficient, expansion order i = 11. The grey vertical lines correspond to the Gauss points of the Chebyshev polynom of order i. . . . . . . . . . . . . . . . . . . . . . . . . 3.23 Computed expansion coefficients of the fourth derivative with respect to space. Comparison between cracked and healthy expansion. From (a) to (d): expansion order i = 2 and 6. The grey vertical lines correspond to the Gauss points of the Chebyshev polynom of order i. . . . . . . . . . . 3.24 Model of the Kirchhoff plate . . . . . . . . . . . . . . . . . . . . . . . . . simply supported plate, with 5% of noise, with the regularization step. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.29 Estimation of ρh D , for a simply supported plate, with 5% of noise, with the regularization step. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Positioning of the sensors patches on the beam structure . . . . . . . . . 4.2 Experimental implementation: cantilever beam, imposed crack and laser measurements at different sample positions (for one sensor array) . . . . 4.3 ρS/EI ID /ρS/EI T H (a) and (∂ 4 v/∂x 4 • u) (b) for the cantilever beam, results based on a experimental data (excitation frequency equal to 1471Hz). On the left, the dashed lines correspond to the mean values for k ≈ 0.8 . 4.4 ρS/EI ID /ρS/EI T H (a) and (∂ 4 v/∂x 4 • u) (b) for the cantilever beam, results based on a experimental data (excitation frequency equal to 1471Hz). On the left, the dashed lines correspond to the mean values for k ≈ 1 . . 4.5 Experimental device for the identification on the polystyrene plate . . . . Basic Principle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.1.1 Differential equation transformation . . . . . . . . . . . . . . . . . 8 1.1.2 Partial differential equation transformation . . . . . . . . . . . . . 9 1.2 Chebyshev polynoms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.2.1 Expansion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.2.2 Differentiation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 Simulation tools . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.1.1 Noise . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.1.2 Monte Carlo test . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.1.3 Results presentation . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.2 Numerical considerations and regularization . . . . . . . . . . . . . . . . 19 2.2.1 STEP 1 : expansion considerations . . . . . . . . . . . . . . . . . Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 A 1D case with multiple identified parameters : the Timoshenko beam . 3.3.1 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 A 1D case using the transient response : a damped bar . . . . . . . . . . 3.4.1 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 A 1D case identifying the damping with the steady state response : cantilever beam . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.1 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 a 1D case with application to discontinuity location . . . . . . . . . . . . 3.6.1 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.7 A 2D case with a single identified parameter: a plate . . . . . . . . . . . Second kind . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	3.28 Estimation of , for a 3.30 Gain (dB) on the estimation of D ρh variable estimate, with 5% of noise . . . . . . . . . . . . . . . . . . . . . ρh D , between least square and instrumental Nomenclature
	xiv 2.1 19 Introduction 1 1 Identification principle 7 1.1 2 Practical issues 17 2.2.2 STEP 2 : Differentiation considerations . . . . . . . . . . . . . . . 22 2.2.3 STEP 3: parameter estimation and regularization process . . . . . 24 3 Numerical applications 29 3.1 A simple 1D case with a single identified parameter : the Bernoulli beam 30 3.1..2.2 Second order derivative (α = 2) . . . . . . . . . . . . . . . . . . . .2.3 Third order derivative (α = 3) . . . . . . . . . . . . . . . . . . . . .2.4 Fourth order derivative (α = 4) . . . . . . . . . . . . . . . . . . . Bibliography This table presents a non-exhaustive list of the symbols used in this thesis. δ i,γ 1 (a, b) Coefficient associated to the expansion of the derivative of order i (first kind) δ i,γ 2 (a, b) Coefficient associated to the expansion of the derivative of order i (second kind) [∆] i,γ 1 Novel derivative operator of order i for first kind Chebyshev polynomials [∆] i,γ 2 Novel derivative operator of order i for second kind Chebyshev polynomials γ Power of the test function u λ Wavelength λ i Expansion coefficient order i λ α,γ Expansion coefficient of α order derivative times u ν Poisson coefficient ω Pulsation ρ Density θ Identified parameter or identified parameter vector c Structural damping D Flexural rigidity [D] Derivative operator E Young modulus f Fonction G Shear modulus i Order I Section inertia k Wave number K Shape fonction L Beam/bar length N Truncation order S Surface section of the beam t Time variable T i First kind Chebyshev polynoms u Test fonction U i Second kind Chebyshev polynoms v Transversal displacement (beam) w Transersal displacement (plate) x Space variable y Space variable 4.1 xv xviii xix

of noise, with the instrumental variable. The mean error on 100 Monte Carlo tests is plotted. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.14 Transient damped bar response (for k = 5) . . . . . . . . . . . . . . . . . 3.15 Estimation of θ 1 et θ 2 for the damped bar, without noise . . . . . . . . .

3.16 Estimation of θ 1 et θ 2 for the damped bar, with 5% noise. The mean error on 100 Monte Carlo tests is plotted. . . . . . . . . . . . . . . . . . . . . . 3.17 Estimation of θ 1 et θ 2 for the damped bar, with 5% noise, with the instrumental variable. The mean error on 100 Monte Carlo tests is plotted. . . 3.18 Estimation of θ 1 and θ 2 with the steady-state of the damped beam, without noise . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.19 Estimation of θ 1 and θ 2 with the steady-state of the damped beam, with 5% of noise . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.20 Estimation of θ 1 and θ 2 with the steady-state of the damped beam, with 5% of noise, with the instrumental variable . . . . . . . . . . . . . . . . . 33.25 Estimation of θ = D ρh , for a simply supported plate, noise free case. . . . 3.26 Estimation of D ρh , for a simply supported plate, with 5% of noise. . . . . 3.27 Estimation of ρh D , for a simply supported plate, with 5% of noise. . . . . 4.6 Identification of D ρh and η at different frequencies, using a simple least square method (without regularization). . . . . . . . . . . . . . . . . . . 4.7 Identification of D ρh and η at different frequencies, using the instrumental variable. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.1.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.2 A 1D case: Influence testing of boundary conditions for the regularization step . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.2.1 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 xvii 3.2.2 3.7.1 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.7.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Experimental applications 4.1 Damage detection on a cantilever beam . . . . . . . . . . . . . . . . . . . 4.1.1 Principle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1.2 Experimental setup . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1.3 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Damped plate identification . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.1 Principle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.2 Experimental setup . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.3 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . Conclusion Application on a torsional non-linear bar Novel operator coefficients .1 First kind . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .1.1 First order derivative (α = 1) . . . . . . . . . . . . . . . . . . . . .1.2 Second order derivative (α = 2) . . . . . . . . . . . . . . . . . . . .1.3 Third order derivative (α = 3) . . . . . . . . . . . . . . . . . . . . .1.4 Fourth order derivative (α = 4) . . . . . . . . . . . . . . . . . . . .2 .2.1 First order derivative (α = 1) . . . . . . . . . . . . . . . . . . . .

Cette thèse est accessible à l'adresse : http://theses.insa-lyon.fr/publication/2013ISAL0095/these.pdf © [C. Chochol], [2013], INSA de Lyon, tous droits réservés

d 0 d 0 d 0 0 0 • • • 0 0 0 0 0 0 0 0 0 0 d 0 d 0 d 0 0 • • • 0 0 0 0 0 0 0 0 0 d 0 d 0 d 0 d 0 • • • 0 0 0 0 0 0 0 0 0 0 d 0 d 0 d 0 d • • • 0 0 0 0 0 0 0 0 0 . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

(n -2, n -1) Cette thèse est accessible à l'adresse : http://theses.insa-lyon.fr/publication/2013ISAL0095/these.pdf © [C. Chochol], [2013], INSA de Lyon, tous droits réservés

Acknowledgements

ACKNOWLEDGMENTS

This work was supported by the Natural Sciences and Engineering Research Council of Canada (NSERC).

Novel operator coefficients

.1 First kind .1.1 First order derivative (α = 1)

Let consider the scalar product between f • u γ and T n :

for n ≥ 1, with

For n = 0 we have :

with

We can then construct the matrix [∆ 1,γ 1 ] made of the δ 1,γ 1 (i, j) coefficients:

(7) Therefore we have:

.1.3 Third order derivative (α = 3)

Let consider the scalar product between f (3) • u γ and T n :

For n = 0 we have:

with

For n = 1 we have:

For n = 2 we have:

where d are non-null values, with d i,j = δ 3,γ 1 (i, j), i being the column indices, and j the array indices. Therefore we have:

is exact if λ 0,γ-3 N +1 = 0, λ 0,γ-3 N +2 = 0 and λ 0,γ-3 N +3 = 0.

.1.4 Fourth order derivative (α = 4)

for n ≥ 4, with

For n = 0 we have :

, 0) for n = 1 we have:

with δ 4,γ 1 (1, 1) = ( 72 -γ)(2δ 3,γ 1 (0, 1) + δ 3,γ 1 (2, 1)) + δ 3,γ 1 (2, 1)

For n = 2 we have:

2) For n = 3 we have:

with

We can then construct the matrix [∆ 4,γ 1 ] made of the δ 4,γ 1 (i, j) coefficients: . . . . . . . . . . . .

where d are non-null values, with d i,j = δ 3,γ 2 (i, j), i being the column indices, and j the array indices. Therefore we have:

is exact if λ 0,γ-3 N +1 = 0, λ 0,γ-3 N +2 = 0 and λ 0,γ-3 N +3 = 0.

.2.4 Fourth order derivative (α = 4)

For n ≥ 4 :