## Abstract

The dynamic responses of a thermo-shielding panel forced by unsteady aerodynamic loads and a classical Duffing oscillator are investigated to detect structural damage. A nonlinear aeroelastic model is obtained for the panel by using third-order piston theory to model the unsteady supersonic flow, which interacts with the panel. To identify damage, we analyse the morphology (deformation and movement) of the attractor of the dynamics of the aeroelastic system and the Duffing oscillator. Damages of various locations, extents and levels are shown to be revealed by the attractor-based analysis. For the panel, the type of damage considered is a local reduction in the bending stiffness. For the Duffing oscillator, variations in the linear and nonlinear stiffnesses and damping are considered as damage. Present studies of such problems are based on linear theories. In contrast, the presented approach using nonlinear dynamics has the potential of enhancing accuracy and sensitivity of detection.

## 1. Introduction

Structural health monitoring plays a significant role in modern air, space, sea and land structures because it provides a means to assess structural integrity online, eliminate manual inspections and may result in a transition from time- to condition-based maintenance. Therefore, structural health monitoring provides numerous benefits including increased safety, increased operational lifetime, reduced cost of maintenance, etc. Presently, a large number of techniques for structural health monitoring are based on the analysis of the vibratory response of structures to various excitations (Zimmerman 1999; Farrar *et al*. 2001). In that context, the vibratory response is interpreted as a feature of the monitored system. Changes in this feature are considered to indicate damage. Most vibration-based damage detection methods monitor changes in the frequencies and modes of vibration (Loh & Tou 1995; Doebling *et al*. 1996; Smyth *et al*. 2000; Farrar *et al*. 2001). These changes are used in algorithms that are well developed primarily for linear systems, and are often inherently limited to linear model forms, such as subspace updating (Abdalla *et al*. 1998; Pappa *et al*. 1998; Zimmerman 2000), wavelets (Amaravadi *et al*. 2001, 2002; Sohn & Farrar 2001; Amizic *et al*. 2002) or Ritz methods (Cao & Zimmerman 1999; Zimmerman 1999). Although highly powerful, these approaches may not be well suited for nonlinear systems, since they do not explicitly account for nonlinear damage scenarios, such as the opening and closing of a crack.

Recently, owing to important advances in nonlinear dynamics and chaos, the application of chaos theory has attracted considerable attention in the structural health-monitoring field and attractor-based features derived from nonlinear time-series analysis are used for identifying damage. Since a chaotic system with one or more positive Lyapunov exponents can have extreme sensitivity to small parametric variation (not only to initial condition), linear systems subjected to chaotic excitation (Trickey *et al*. 2002; Nichols *et al*. 2003*a*–*d*) and (nonlinear) chaotic systems (with or without excitation) (Epureanu & Yin 2004*b*; Epureanu *et al*. 2004*b*,*c*) have been explored to show that the use of nonlinearity holds a great potential. Furthermore, the enhancement of the nonlinearity of linear or weakly nonlinear systems by means of nonlinear feedback excitation has been shown to provide significant advantages such as increased sensitivity (Epureanu & Yin 2004*a*; Epureanu *et al*. 2004*d*). A key point addressed in nonlinear approaches is the definition of a feature utilized to quantitatively characterize the geometry of the attractor in state space. A variety of attractor-based metrics has been presented for identifying the changes in the geometric feature of the attractor owing to the changes in the system parameters (i.e. damage). Global geometric properties, such as the maximal Lyapunov exponent or the attractor dimensions, have been used in earlier studies, but often have difficulties in detecting the subtle changes in the attractor geometry owing to damage. Instead, a few approaches have been developed using a local analysis of trajectories within an attractor. For example, a scalar tracking metric (Chatterjee *et al*. 2002; Chelidze *et al*. 2002; Cusumano *et al*. 2002; Chelidze & Cusumano 2004) and an average local attractor variance ratios (Todd *et al*. 2001; Nichols *et al*. 2003*b*) have been used to detect single damage, while multidimensional damage in hierarchical systems and fatigue damage detection have also been tackled (Chelidze 2004; Chelidze & Liu 2005). An alternative approach is to statistically characterize the distribution of points in an attractor. The probability distribution of points in an attractor may be approximately represented by a linear combination of orthogonal basis functions, while a scalar metric based on the coefficients of the orthogonal functions (series) is used as an indicator (Trendafilova 2003*a*,*b*). Unfortunately, a rational choice for the orthogonal functions has not been presented in the literature (Trendafilova 2003*a*,*b*). Moreover, these proposed scalar metrics may be appropriate for detecting a single damage. However, they may not detect multiple simultaneous damages because only a scalar number is provided by the metric (indicator). To address these limitations, the use of the probability density functions has been enhanced by defining an observed feature as the difference between the probability density functions for the damaged and the healthy systems (Epureanu & Yin 2004*b*), which clearly shows the effects of damage. Moreover, the shape and the magnitude of probability density function difference may be correlated with the location and the level of damage using pattern recognition techniques based on proper orthogonal decomposition (POD), which provides the basis for detecting multiple simultaneous damages, as well as a rational choice for basis functions for the probability distribution functions (Epureanu & Yin 2004*b*).

In the present work, we propose a new attractor-based damage technique, which is based on extracting sensitivity vectors determined by analysing the short time variation between the two nearby trajectories corresponding to the healthy and damaged systems. The concept of sensitivity vectors has been used in the context of chaos suppression (and system identification for chaos suppression; Epureanu & Dowell 1998, 2000; Epureanu *et al*. 1998) as well as damage detection (Cusumano *et al*. 2002; Epureanu & Yin 2004*c*; Epureanu & Hashmi in press). The features extracted by means of sensitivity vectors may indicate the location, the level and the extent of damage. POD (Azeez & Vakakis 2000; Feeny 2002) is also used to identify dominant coherent structures present in the system dynamics. The decomposition of the dynamics into these coherent structures reduces the number of measurements needed for analysis, which is important for structural health monitoring of high-dimensional systems. Finally, POD is exploited to identify the number of linearly independent shapes among the extracted features associated with the attractor changes caused by distinct damages. Hence, a basis for detecting multiple damages is developed. Herein, the dynamic responses of a thermo-shielding panel and a Duffing oscillator are investigated. For the thermo-shielding panel, damage is considered as a local reduction in bending stiffness in the panel (caused by yielding and the propagation of a crack in the material). For the Duffing oscillator, variations in the linear and the nonlinear stiffnesses, and damping are considered as damage. The results presented show that the proposed approach is an effective tool for detecting multiple simultaneous damages.

## 2. Modelling and numerical discretization

A simple Duffing oscillator and an aeroelastic system composed of a panel subjected to an unsteady supersonic flow (depicted in figure 1) are used to demonstrate the proposed approach. The Duffing oscillator is chosen as an archetypal strongly nonlinear system exhibiting limit cycle oscillations and chaos. The thermo-shielding panel is chosen as a more complex system, modelled as a nonlinear, one-dimensional, homogeneous, isotropic and elastic thin plate with spring-supported end points, and the flow is modelled by third-order (unsteady) piston theory (Lighthill 1953). Distinct from other studies of panel flutter (Epureanu *et al*. 2004*a*,*b*), the third-order piston theory accounts for flow nonlinearity. Finally, a local reduction in bending stiffness of the plate in a small region around various locations along the panel is considered as a damage scenario.

The governing equation for the complete aeroelastic system is obtained as follows:(2.1)where *ρ*_{∞} and *U*_{∞} are the upstream far-field density and velocity of the flow, respectively, *M*, the local Mach number, *γ*, the specific-heat ratio of the fluid (*γ*≅1.4), *h*, the panel thickness, *l*, the panel length, *W*′ and , the spatial and time derivatives of *W*, respectively, *η*_{0}, the initial axial strain, Δ*p*_{s}, the static pressure difference across the plate (at zero flow velocity), *ρ*, the mass density of the panel and *D* is a coefficient characterizing the bending stiffness of the panel (*D*=*Eh*^{3}/12 (1−*ν*^{2}) with *E* being Young's modulus and *ν* Poisson's ratio). The bending stiffness and Young's modulus in the damaged region of the panel are denoted by and , respectively, and are distinct from *D* and *E* for the healthy regions.

Next, the non-dimensional variables, *x*, *w* and *τ*, are defined as *x*=*X*/*l*, *w*=*W*/*h* and , and equation (2.1) can be non-dimensionalized as follows:(2.2)where *w*′=∂*w*/∂*x* and . The corresponding non-dimensional coefficients in equation (2.2) are as follows: in-plane pre-load *R*_{x}=*Ehη*_{0}*l*^{2}/*D*; bending–stretching coefficient *S*=*Eh*^{3}/(2*D*)=6(1−*ν*^{2}); flow velocity and mass ratio *μ*=ρ_{∞}*l*/(ρ*h*). The coefficient *S*_{r} is the stiffness reduction factor, , characterizing the damage level in the panel. For example, *S*_{r} equals unity when no damage is present. In addition, the smaller the value of *S*_{r}, the larger the damage.

The boundary conditions can also be rewritten in non-dimensional form aswhere *R*_{u}=*D*/(*k*_{u}*L*^{3}) and *R*_{d}=*D*/(*k*_{d}*L*^{3}) with *k*_{u} and *k*_{d} being the stiffnesses of the upstream and downstream mounting points, respectively.

The numerical approach for computing the dynamics of the aeroelastic system is to discretize the governing partial differential equation. The computational domain is uniformly divided into *m* sections. All the spatial derivatives in the governing equations and the boundary conditions satisfied at the discretization points may be replaced by the central finite difference formulation of second-order accuracy. In addition, one ghost node is also needed in a central finite difference formulation at each end of the panel.

The vector of displacements, * w*, at all discretization points is defined as

*=[*

**w***w*

_{−1},

*w*

_{0},…,

*w*

_{m+1}]

^{T}, while

*denotes the force vector corresponding to the static transverse pressure. By introducing a new vector , equation (2.2) (satisfied at the discretization points) can be written in the form of a system of first-order differential equations in state space as follows:(2.3)where . The resulting discretized system of ordinary differential equations is integrated in time using an adequate numerical scheme for time marching.*

**p**## 3. Damage detection technique

Most of the present damage detection methods are well developed primarily for linear system models and use linear approaches such as frequency analyses. However, they have difficulties detecting single or multiple damages of small level because they eliminate or minimize the nonlinearity in the system despite the fact that nonlinearities may be more sensitive to damage. Hence, present monitored features used for identifying changes in different parameters (damages) have low sensitivity. In contrast, the technique proposed herein does not linearize the system, and instead, focuses on the trajectory on the attractor of dynamics. We investigate two nearby trajectories corresponding to a healthy and a damaged system. The variation between these two trajectories in time owing to parametric variation may be extracted and used as a feature to detect damage. The theoretical background of the proposed approach is presented next.

Consider an attractor of the dynamics of a healthy system characterized by a general nonlinear flow expressed in first-order differential form as(3.1)where **x**_{h} denotes the trajectory of the healthy system, corresponding to its healthy parameter *p*_{o} on the attractor in state space. The parameter *p* is considered to undergo a variation, *δp*, caused by damage occurring in the system. The trajectory on the attractor of the damaged system, **x**_{d}, satisfies(3.2)Suppose the two trajectories, **x**_{h} and **x**_{d}, are very close, then(3.3)where *δ x*(

*t*) is a small variation between these two trajectories (as shown in figure 2). Substituting equation (3.3) into equation (3.2) and using a Taylor series expansion of equation (3.2), one obtains(3.4)

Next, considering *δ x* and

*δp*to be small, higher order terms may be neglected. Thus, the variation

*δ*(

**x***t*) is governed by the following equation(3.5)where and . Mathematically, if the trajectory of the healthy system,

**x**_{h}(

*t*), is known in an explicit form, the solution,

*δ*(

**x***t*), may be obtained by integrating equation (3.5) from the initial condition,

*δ*(0), as follows:(3.6)where the state transition matrix

**x***depends on*

**F***(*

**A***t*), and the sensitivity vector depends on both

*(*

**A***t*) and

*(*

**b***t*).

The initial trajectory variation, *δ x*(0), is obtained by observing two nearby trajectories in state space as shown in figure 2. The time is initialized to 0 at the state corresponding to the start of the observation time-interval. In general, it is difficult to obtain the solution in explicit form (e.g. equation (3.6)) for most complex systems. Nevertheless, equation (3.6) provides a clear understanding of how the variation,

*δ*(

**x***t*), is influenced by the initial condition,

*δ*(0), and the parametric variation,

**x***δp*. It can be noted that equation (3.6) is valid for any dynamics, including limit cycle oscillations and chaos, provided

*δ*(

**x***t*) and

*δp*are small enough. Also, the small value of

*δp*implies that no bifurcations take place and the attractors of the healthy and the damaged systems (formed by the evolution of these two corresponding trajectories) have similar geometric shapes.

From the viewpoint of detecting multiple damages (which may be modelled as the simultaneous changes in different system parameters), exploring the sensitivity vector, , in equation (3.6) is more significant and interesting, since the shape of may be related to the specific parameter representing the damage location, while the magnitude of *δp* may be used to determine the damage level. However, the sensitivity vector, , is not straightforward to obtain because the only available data are the time-series of the two trajectories, **x**_{h}(t) and **x**_{d}(t), calculated by integrating equations (3.1) and (3.2) directly, i.e. only *δ x*(

*t*) can be used. The key of the proposed damage detection technique is to find the specific shape of versus time directly from time-series of measured

*δ*(

**x***t*).

The scalar quantity *p* in equation (3.1) represents a specific parameter of the system. The influence of this specific parameter on the trajectory variation can be explored using the associated sensitivity vector based on the formulation in equations (3.1)–(3.6). Of course, for each different parameter variation, the corresponding sensitivity vector can be found. Multiple damage detection is based on identifying and distinguishing these sensitivity vectors, which are combined in trajectory variations owing to multiple parametric variations. Hence, if one defines * p* as the vector of parameters, then

*(*

**q***t*) can be interpreted as a sensitivity

*matrix*consisting of sensitivity vectors for each parameter in the vector of parameters,

*. Next, we analyse the cases of limit cycle oscillations and chaotic dynamics, separately.*

**p**### (a) Limit cycle oscillations

In the case where the monitored system is an autonomous system exhibiting limit cycle oscillations, the trajectory of the system is invariant to a shift of the trajectory within itself (i.e. a shift in time or phase) and the largest Lyapunov exponent of the dynamics is 0. Hence, * F*(

*t*)

*δ*

*(0) can be modified by shifting its phase by*

**x***δt*

_{0}so that equation (3.6) becomes(3.7)where is the velocity vector of the flow for the healthy system and

*δt*

_{0}the phase difference (or shift in time) between the trajectories of two limit cycles (for a healthy system and a damaged one). To reduce the sensing requirements, we consider that only the position and the velocity of one degree of freedom (DOF, e.g. the

*i*th DOF) in the system are measured. The difference in position at the

*i*th DOF between the two trajectories can be written from equation (3.7) as follows:(3.8)where

*δ*

**x**_{i}(

*t*), and are the

*i*th elements of the vectors

*δ*(

**x***t*), and , respectively. Next, time-series of equation (3.8) can be rewritten in vector form by sampling in time as follows:(3.9)where

*δ*

**x**_{i}=[

*δx*

_{i}(

*t*

_{0})

*δx*

_{i}(

*t*

_{1})…

*δx*

_{i}(

*t*

_{n})]

^{T}, with , for

*k*=0, …,

*n*, and . The index

*n*represents the index of the number of samples collected over the portion of the time-series, where the trajectory variation

*δ*

**x**_{i}is computed. The length of the sampling time is preferably longer than one period of a limit cycle, and therefore the shape of the trajectory variation can be well characterized. Herein, is the vector representing the time-series of the

*i*th component of the sensitivity vector. Since

*δt*

_{0}is unknown, cannot be directly obtained from equation (3.9), but the effect of the phase

*δt*

_{0}, on the measurements, can be eliminated by extracting the component of the vector,

*δ*

**x**_{i}, which is perpendicular to . One obtains(3.10)where represents the component of the vector perpendicular to the vector , and is referred to as the

*modified trajectory variation*. Although cannot be completely determined, has a shape, which is correlated with the magnitude and the type of parametric variation

*δp*. Hence, in the proposed damage detection technique, we exploit the modified trajectory variation as a feature to detect different damage locations and levels, while the system exhibits limit cycle oscillations.

### (b) Chaotic dynamics

Distinct from limit cycle oscillations, in the case of chaotic dynamics, the presence of positive Lyapunov exponents causes *δ x*(0) to grow exponentially. Thus, the technique developed for limit cycle oscillations has to be modified. It may be noticed that the easiest way to obtain the sensitivity vector, , from

*δ*(

**x***t*) is to find the two nearby trajectories for the healthy and the damaged systems which intersect in the state space such that

*δ*(0) in equation (3.6) is exactly equal to

**x****0**. Although the ergodicity of the attractors guarantees that a point of the attractor is visited infinitesimally close by the system dynamics (with a non-zero probability), the attractors in the case of damage detection are distinct (despite their large overlap). Hence, the ergodicity does not guarantee that

*δ*(

**x***t*)=

**0**may occur with a non-zero probability. Thus, eliminating

**(**

*F**t*)

*δ*(0) from equation (3.6) requires a special treatment. Herein, we propose the method of averaging trajectories to tackle this difficulty. The details of this method are introduced in the following.

**x**Assume that the attractors of the healthy and damaged systems have regions where they overlap in the state space. Within this overlapping region, one may find that the trajectory on the attractor of the damaged system may pass close to the reference trajectory of the healthy system several times owing to the ergodicity of the two chaotic attractors. Consider that the reference trajectory of the healthy system passes though the centre line of a tube as shown in figure 3. Then, other nearby trajectories of the damaged systems also pass within this tube as shown in figure 3. The radius of the imagined tube may represent an upper bound for the Euclidean norm of *δ x*(

*t*). A very small value of this radius guarantees that the variations between the trajectories of the damaged and the healthy systems in figure 3 can be accurately approximated by equation (3.6), and expressed as(3.11)where

**x**_{h}(

*t*) and represent the reference trajectory of the healthy system and the

*j*th trajectory of the damaged system with

*j*=1, 2, …,

*l*, respectively (e.g. assume that a total of

*l*trajectories of the damaged system are found within the tube). It can be noted that the state transition matrix

*(*

**F***t*) and the sensitivity vector

*(*

**q***t*) in equation (3.11) are identical for all trajectories within the tube (if the radius of the tube is small enough).

Next, we may average the trajectories of the damaged system within the tube to obtain the variation, *δ x*

^{ave}(

*t*), between the averaged trajectory of the damaged system, , and the reference trajectory of the healthy system,

**x**_{h}(

*t*), as shown in figure 3 by the following derivation. First, equation (3.11) can be multiplied by any weighting coefficients

*c*

_{j}and summed together for all the trajectories of the damaged system (

*j*=1, 2, …,

*l*) as follows:(3.12)Also, , and can be denoted by

*δ*

**x**^{ave}(

*t*), and

*δ*

**x**^{ave}(0), respectively. Equation (3.12) can be rewritten as(3.13)Next, the coefficients

*c*

_{j}can be determined such that(3.14)Therefore, equation (3.13) can be simplified to obtain the sensitivity vector as follows:(3.15)Hence, the sensitivity vector , a significant feature associated with damage detection, can be directly obtained by computing (i.e.

*δ*

**x**^{ave}(

*t*)) without knowing the state transition matrix. Not requiring the knowledge (or an estimation) of the state transition matrix is an important advantage of the method of averaging trajectories over other techniques (Chatterjee

*et al*. 2002; Chelidze

*et al*. 2002; Cusumano

*et al*. 2002; Chelidze & Cusumano 2004), as discussed in §5 below. Finally, a vector composed of the time-series of the

*i*th component of the sensitivity vector, , can be obtained by sampling the

*i*th component of the

*averaged trajectory variation*

*δ*

**x**^{ave}(

*t*) in time as follows:(3.16)The sensitivity vector can then be utilized as a feature to detect different damage locations and levels, which may be modelled by different types of parametric variation.

Since the method of averaging trajectories is developed based on the assumption that the collected trajectories of the damaged system are close to the reference trajectories of the healthy system during a short period of time, all the state variables may be required to be monitored to ensure that this condition can be satisfied. However, this may be difficult to attain in practice for the case of high-dimensional systems. Therefore, the proposed method is enhanced by the use of POD to identify dominant coherent structures present in the dynamics of high-dimensional systems. The decomposition of the dynamics into dominant coherent structures and the extraction of dominant dynamics (containing most energy) reduce the demand of measurements required for this technique. More details about how to apply the POD method can be found in previous studies (Sirovich 1987; Feeny 2000, 2002; Kappagantu & Feeny 2000; Epureanu *et al*. 2004*a*,*b*). Consider that the dynamics of the aeroelastic system, governed by equation (2.3), is approximated by a linear combination of dominant coherent structures aswhere the generalized coordinates, *a*_{i}, and their time derivatives, , are used to describe the amplitude and the velocity of the coherent structures, and *N*, the chosen number of dominant structures, *ϕ*_{i}, needed for good approximation of the dynamics. In state space representation, the panel dynamics is characterized by(3.17)whereUsually, the number of generalized coordinates used (e.g. 2*N*) is much smaller than that of the state variables of the system. In addition, time histories of generalized coordinates (*Z*_{i},*i*=1, 2,…,2*N*) can be obtained by measuring the displacement and the velocity at *N* different locations in the panel, once the dominant modes have been identified. Thus, instead of measuring all state variables, the number of measurements required for implementing the damage detection technique presented here can be effectively reduced by identifying dominant coherent structures based on POD.

### (c) Detection of multiple damages

To distinguish between different types of damage and different damage locations (single or multiple), we employ POD again. However, distinct from the previous section, here POD is used in the parameter space. Hence, we are able to identify distinct damages by the number of linearly independent shapes among the extracted features such as in equation (3.10) when the system exhibits limit cycle oscillation or in equation (3.16) when chaotic dynamics is present. Some of these feature vectors have similar shapes and are proportional in magnitude. These feature vectors are correlated to or caused by changes occurring in the same parameter, but with different magnitudes (of change) *δp*. Other feature vectors are linearly independent and correspond to changes in distinct parameters. Hence, in that case, linearly independent shapes or are obtained. The identified number of linearly independent shapes among the features obtained (based on POD) may be correlated to the number of locations where damage may be detected or the number of identified damage types. With this aim, the collected features are grouped as column vectors in a matrix * R*, and then a correlation matrix

*is formed as*

**C****C**=

**R**^{T}

*. Next, the eigenvalues of the correlation matrix are calculated to determine how distinguishable these damages are. The linearly independent shapes correspond to distinguishable damages and are identified by a clear separation of their corresponding (dominant) eigenvalues from the rest of the eigenvalues.*

**R**## 4. Results

To demonstrate the viability and performance of the proposed damage detection technique, we apply it to two general nonlinear systems exhibiting limit cycle oscillations and chaotic dynamics. A Duffing two-well oscillator and a nonlinear panel subjected to unsteady supersonic flow are investigated.

### (a) Limit cycle oscillations

To demonstrate the proposed approach for the case of a system exhibiting limit cycle oscillations, we consider an archetypal example of a Duffing two-well oscillator characterized by the following governing equation:where the parameters, *ϵ*, *k* and *β*, have the following values for the healthy system: *ϵ*_{0}=0.25, *k*_{0}=1, and *β*_{0}=1, while the forcing amplitude is *F*=0.18. Damage in the system is considered to generate changes in these parameters. Figure 4 shows the small difference between the two limit cycles corresponding to a healthy system (*k*=*k*_{0}) and a damaged system (*k*=0.99*k*_{0}). First, in state space spanned by *u* and (denoted by * x*), the initial sampled point in the trajectory of the attractor for the healthy system,

**x**_{h}(0), is chosen arbitrarily. Then, the initial sampled point

**x**_{d}(0) of the trajectory corresponding to the attractor of the damaged system is selected based on the nearest distance to

**x**_{h}(0). Time-series of the two trajectories (starting from the chosen points,

**x**_{h}(0) and

**x**_{d}(0)) are collected and sampled in time. The samples of

*u*for the healthy and the damaged systems are denoted by

**u**_{h}and

**u**_{d}, respectively. Next, the variation of

*u*between these two trajectories can be obtained by subtracting one from the other to obtain

*δ*=

**u**

**u**_{d}−

**u**_{h}. Figure 5 shows the trajectory variation,

*δu*, between the healthy and the damaged systems with different changes in parameter

*k*(0.5 and 0.25% reduction of

*k*

_{0}). In figure 5, one may note that there is no specific shape associated with the variation in a single parameter because of the influence of the unknown phase

*δt*

_{0}in equation (3.9). This phase exists at the start of sampling of the two trajectories for the healthy and the damaged systems, and causes the variation in the shape of

*δu*when

*δt*

_{0}varies. Nonetheless, the measured trajectory variation

*δu*may be modified as presented in §3

*a*. In the proposed approach, the trajectory variation

*δu*is modified by compensating for the time phase based on equation (3.10). The modified trajectory variation provides a shape for , which does not depend on

*δt*

_{0}, and is correlated with the type of parametric variation analysed. Figures 6–8 show two distinct levels of variation (0.5 and 0.25% reduction) in three different parameters (

*k*,

*β*and

*ϵ*). One may observe that, in each plot, the magnitudes of the shapes are proportional to the magnitude of the parameter variation and the shapes for different parametric variations are clearly distinct. This correlation between shapes and parameter variations is exploited herein for detecting different types of damage. In addition, it can be noted that the proposed technique is applicable for small parameter variations (where equation (3.4) holds). Figure 6 shows that the shape of the trajectory variation may change and the proportionality may break down when larger parametric changes occur, as is the case of 1% reduction in parameter

*k*

_{0}(dotted line).

### (b) Chaotic dynamics

Next, a nonlinear aeroelastic system composed of a panel subjected to unsteady supersonic flow is considered. This system exhibits highly complex dynamics such as limit cycle oscillations and chaos (Epureanu *et al*. 2003, 2004*a*). Since coherent structures (dominant modes) contain most of the energy and are the most significant components of the dynamics, the projections of the full-dimensional attractors of the dynamics along the coherent structures are worth investigating to observe the changes in the attractors owing to the parametric variations (Epureanu *et al*. 2004*b*). In addition, by identifying the coherent structures, demand on the required measurements is reduced by using the proposed technique discussed in §3*b*. Herein, a set of parameters for the aeroelastic system exhibiting chaotic dynamics have been selected as *λ*=150, *ν*=0.252, *μ*/*M*=0.01 and *R*_{x}=−7*π*^{2}. First, POD is used for identifying coherent structures present in the dynamics. The eigenvalues *λ*_{i} of the correlation matrix based on the collected snapshots of the state of panel are shown in figure 9. The rapid decrease in the magnitude of the eigenvalues *λ*_{i} indicates that a strong spatial correlation exists between the dynamics of various points on the panel. Also, the first three eigenvalues are obviously separated from the rest of the eigenvalues, which implies that the first three dominant modes might be sufficient when used to provide an approximation for the morphing of the attractors owing to damage. Figure 10 shows the Poincaré plots of the attractor of the first coherent structure of a healthy panel and panels with damage at different locations. Here, the differences between the Poincaré plot of the attractor of a healthy panel and those of damaged panels (cases (1)–(5)) are quite subtle. The main reason for the subtlety is that the high-dimensional attractor of the dynamics of the panel is projected onto a low-dimensional subspace to obtain a low-dimensional attractor by using POD. Thus, characterizing and quantifying the differences in the distribution of points of the attractors shown in figure 10 need an elaborate attractor analysis. Such an analysis has been first proposed using probability density *difference* (Epureanu & Yin 2004*b*). A probability density function is defined as , where * x* characterizes the location of the point at which the density

*P*is computed in the state space (two-dimensional Poincaré plot),

*n*is the number of points of an attractor in a circle of radius

*r*around

*and*

**x***N*is the total number of points in the sampled attractor. Thus, the difference between the probability density functions for the damaged and the healthy panels is calculated. Figure 11 shows that the probability density function difference, Δ

*P*, corresponding to different locations of damage has its own specific shape, which means that the clouds of points in each corresponding Poincaré plot in figure 10 morph (deform and move) in a different way for each location of damage. As a result, using attractor morphing, the sensitivity of chaotic dynamics of the panel to parametric variations is used to detect damage occurring in the panel and also to enhance the capability of detecting damage at different locations.

Different from using statistical descriptions of the cloud of points in the attractor to characterize the changes in the geometric feature of the attractor based on probability density difference, the attractor morphing caused by parametric variations (damages) can also be detected by tracking the time-series of sampled sensitivity vectors of the system. The approach based on sensitivity vectors is more detailed, in that it requires a faster sampling rate for measurements. However, this approach requires fewer measurements in state space (i.e. fewer initial states in figure 3). As discussed in §3*b*, the sensitivity vector can be determined by calculating the short time variation of the averaged trajectory of the damaged panel from the reference one (of the healthy panel). To demonstrate this approach, the attractor trajectories of the healthy panel can be generated by integrating equation (2.3). Then, the dynamics of the panel can be decomposed into the subspaces spanned by the identified dominant coherent structures based on POD. Herein, the generalized coordinates and their time derivatives corresponding to the first three dominant coherent structures (i.e. six generalized coordinates in state space representation) are used, since the first three eigenvalues are clearly separated from the rest of the eigenvalues (as shown in figure 9).

To simulate data for the approach of averaging trajectories, first, the reference trajectory of the healthy panel is chosen based on the regions where the attractors of the healthy and damaged panels overlap in state space. Next, the generalized coordinates for the initial state of the reference trajectory at the start of the sampling are perturbed randomly (with 1% changes of magnitude) to generate a set of nine new values of generalized coordinates. These new states and the associated coherent structures are employed to construct a set of initial states (which are assumed to be in the attractor of the damaged panel). Then, the group of trajectories of the damaged panel (inside the tube shown in figure 3) can be simulated by integrating equation (2.3) of the damaged panel starting from this set of initial states. Finally, the group of trajectories of the damaged panel is averaged with the coefficients satisfying equation (3.14), and the time-series of sampled sensitivity vectors can be determined by calculating the time variation between the averaged trajectory of the damaged panel and the reference trajectory of the healthy panel.

To investigate the error in the calculation of sensitivity vectors (caused by the approximations involved in the method of averaging trajectories), the trajectory of the damaged panel is obtained exactly by time integration of equation (2.3) (from the same initial state as for the healthy reference trajectory). The exact sensitivity vectors and the trajectory variation are determined by subtracting the healthy reference trajectory from the exact damaged one. The difference between the averaged trajectory variation and the exact one was obtained. This error was less than 1% of the magnitude of the sensitivity vector itself, which means that the sensitivity vectors can be robustly calculated by using the proposed approach. In this part of the present study, only the time-series of the *i*th component of the sensitivity vector (or the averaged trajectory variation) is used as a feature for detecting the location, extent and level of damage in the panel. For example, the time variation of the displacement at the quarter point between the damaged and the healthy panels (*w* at *x*=1/4) is observed in the following cases.

Figure 12 shows the averaged trajectory variation between the healthy panel and the those with different levels of damage (5, 10, 20 and 40% loss in bending stiffness) along 1.25% length of the panel around *x*=1/2. These curves have nearly identical shapes and their magnitude is related to their damage level (the higher the damage level, the larger the magnitude). Nonetheless, the proportionality to the level of damage may not hold when the damage level is high (e.g. in the case of 40% stiffness loss), but the shape associated with the location of damage (for *x*=1/2 in this case) remains unchanged.

In addition, to explore the correlation among the shape (of the trajectory variation), the damage location and the damage extent, three damage locations (around *x*=1/4, *x*=1/2 and *x*=3/4) are considered for small and large damage extents. A 10% stiffness loss along 1.25% of panel length is referred to as small damage, while a 10% stiffness loss along 12.5% of panel length is referred to as large damage. Also, three different reference trajectories in the attractor of the healthy panel are considered. Figure 13 shows the three initial states of the panel associated with these three reference trajectories chosen. Figure 14 shows that the shapes of the averaged trajectory variation for distinct damage locations differ from each other for all cases. This characteristic provides a high potential to detect the location of damage. Moreover, in each case, the shapes remain the same for the small and the large extents of damage as long as the damage location is the same. It can be noted that the magnitude of the shape for the large extent of damage is almost 10 times larger than the one for the small extent of damage, which reflects the validity of linearity (for short times), since the damaged panel length for the large damage extent is 10 times longer than the one for the small damage extent. Finally, the snapshots of the shape of the trajectory variation shown in figure 14 for the same damage condition are grouped together. Thus, a total of six sets of snapshots are obtained. These snapshots are analysed using POD to identify the number of linearly independent snapshots. The number and the shape of these snapshots are then related to the number of locations of damage. For a panel with damages around *x*=1/4, *x*=1/2 and *x*=3/4, figure 15 shows that the first three dominant eigenvalues of the correlation matrix are clearly separated from the other eigenvalues, which indicates that there are three linearly independent snapshots corresponding to the damage occurring at these three locations. A similar result is obtained when the difference probability density functions are used as snapshots in the POD approach. Herein, nine sets of snapshots corresponding to nine damaged panels are collected. These nine damage cases are respectively modelled as three damage levels (a 50, 55 and 60% siffness loss) occurring along 1.25% of panel length around three locations (*x*=1/4, *x*=1/2 and *x*=3/4). Figure 16 shows the separation of the first three eigenvalues from the rest, which corresponds to three damage locations.

## 5. Discussion and conclusions

Structural damage in a thermo-shielding panel forced by unsteady aerodynamic loads and parametric variations in a Duffing oscillator has been identified by the analysis of nearby trajectories in state space. The nonlinear aeroelastic system was modelled by using the third-order piston theory for the unsteady supersonic flow and a nonlinear model for the panel. Various types of damage were modelled as local reductions in bending stiffness in small regions at various locations along the panel, and as variations in linear and nonlinear stiffnesses and damping of the Duffing oscillator. POD was used to extract the dominant coherent structures of the dynamics of the panel. The morphing of the attractor of the dynamics caused by damage can be revealed by both the probability density difference functions and the trajectory variations within the attractors. Damages of various locations, extents and levels were detected. The correlation of the observed features to different types of damage can be found based on POD. Most of the present studies of such problems are based on linear theories. In contrast, the results presented have been obtained using nonlinear dynamics and show high sensitivity.

The proposed approach is developed based on the assumption that detecting the small parameter variations is of interest. Hence, trajectories are collected such that the higher order terms of trajectory variation *δ x* may be neglected, as shown in equation (3.6). Thus, the closeness of trajectories is important for the proposed method. The measure for this closeness is the radius of the tube in figure 3. The values for this radius are dependent on each particular application and on the type of measurements performed. Hence, as often encountered in experimental work, a preliminary analysis has to be performed prior to the implementation of the proposed approach. In addition, the robustness of limit cycle oscillations and of chaotic attractors to the parameter variations ensures that such tubes (and portions of trajectories) do exist. Furthermore, increasing the radius of the tube or collecting longer time-series allows for collecting more trajectories and also for increasing the length of the tube. These factors may be used to reduce the influence of noise through the process of trajectory averaging.

The high sensitivity of the proposed approach is not influenced by noise directly. Nonetheless, possible experimental noise may have an influence on the *accuracy* of the proposed approach. Specifically, noise levels that are larger than the radius of the tube in figure 3 probably lead to inaccurate results. Since the radius of this tube is application-dependent, quantitative results regarding the influence of noise are application-dependent as well.

The proposed approach is specifically designed not to use the state transition matrix. This is different from other approaches and is also advantageous. The major advantages are at least fourfold. First, eliminating the need for computing an estimate for the state transition matrix drastically relaxes measurement requirements. For example, only one coordinate of the system needs to be measured, as demonstrated in figures 5–8 and figures 12–14, where only a scalar trajectory variation was necessary to detect multiple damages. This major advantage is owing to the fact that not using the state transition matrix allows the projection of the trajectory variation *δ x* onto a simpler, easier to use and lower dimensional subspace (which can be as low as one-dimensional). Second, eliminating the need for computing an estimate for the state transition matrix drastically reduces the number of samples (time-series measurements) that are needed by as much as two or three orders of magnitude (from millions to thousands of tenths of thousand) compared to other techniques. Hence, not only the computation time, but also the time needed to collect data is reduced. Third, the proposed approach is optimal and may be used for noisy measurements, in that the fit of trajectories through a fitted weighting is optimal in the same way as the fitting for the state transition matrix. Finally, an advantage over other techniques is that it is time saving. This is important because the computational time required by other techniques can become prohibitively high, especially when a system with more than one DOF is tackled.

Two separate cases have been discussed: systems exhibiting limit cycle oscillations and those exhibiting chaos. For systems exhibiting limit cycle oscillations, one state and its time derivative are measured for the healthy and the damaged systems. The variation between time-series of this measured states are characterized by equations (3.8) and (3.9). Next, the time-series are changed to obtain the modified trajectory variation, which depends only on parametric variations. This key idea is based on the relation shown in equation (3.10). For systems exhibiting chaos, four key aspects have been discussed. First, the number of the states and their time derivatives to be measured are determined by the number of dominant modes (that are considered sufficient for describing the dynamics). This number of dominant modes can be identified (for example) by applying POD on the dynamics from numerical simulations. Second, the trajectory of the (full) system can be reconstructed from the time-series of the (small number of) generalized coordinates that correspond to the dominant modes. The relationship between the generalized coordinates and the dynamics of the full physical system may be obtained by using the POD projection matrices. This relationship allows for extracting the generalized coordinates from measurements. Third, trajectories of a damaged system may be collected experimentally until they fall within a tube around the reference trajectory of the healthy system (as shown in figure 3). Finally, the variation between the average of the selected trajectories of the damaged system and the reference trajectory of the healthy system in equation (3.15) can be obtained by solving equation (3.14) for the coefficients *c*_{j}. The resulting averaged trajectory variation is the time-series for the sensitivity vectors, and reveals the dependence of the dynamics on the parametric variations. A more detailed discussion of experimental results is forthcoming in a future paper.

## Acknowledgements

The authors wish to acknowledge the National Science Foundation (CAREER program) and Professor Masayoshi Tomizuka and Dr Shih-Chi Liu (program directors) for the generous support of this work.

## Footnotes

One contribution of 15 to a Theme Issue ‘Exploiting chaotic properties of dynamical systems for their control’.

- © 2006 The Royal Society