System and method for detecting structural damage

[0001] This application is a continuation of U.S. application Ser. No. 10/849,571 filed May 20, 2004. The disclosures of the previous application are incorporated by reference herein. This application further claims priority to U.S. Provisional Patent Application Ser. No. 60/471,813 filed May 20, 2003 and further claims priority to U.S. Provisional Patent Application Ser. No. 60/512,656 filed Oct. 20, 2003, the latter of which is incorporated herein by reference.

[0002] 1. Field of the Invention

[0003] The invention relates to a method and apparatus for detecting structural damage, and, more specifically, to a method and apparatus for detecting structural damage using changes in natural frequencies and/or mode shapes.

[0004] 2. Background of the Related Art

[0005] Damage in a structure can be defined as a reduction in the structure's load bearing capability, which may result from a deterioration of the structure's components and connections. All load bearing structures continuously accumulate structural damage, and early detection, assessment and monitoring of this structural damage and appropriate removal from service is the key to avoiding catastrophic failures, which may otherwise result in extensive property damage and cost.

[0006] A number of conventional non-destructive test (NDT) methods are used to inspect load bearing structures. Visual inspection of structural members is often unquantifiable and unreliable, especially in instances where access to damaged areas may be impeded or damage may be concealed by paint, rust, or other coverings. Penetrant testing (PT) requires that an entire surface of the structure be covered with a dye solution, and then inspected. PT reveals only surface cracks and imperfections, and can require a large amount of potentially hazardous dye be applied and disposed of. Similarly, magnetic particle testing (MT) requires that an entire surface of the structure be treated, can be applied only to ferrous materials, and detects only relatively shallow cracks. Further, due to the current required to generate a strong enough magnetic field to detect cracks, MT is not practically applied to large structures. Likewise, eddy current testing (ET) uses changes in the flow of eddy currents to detect flaws, and only works on materials that are electrically conductive. Ultrasonic testing (UT) uses transmission of high frequency sound waves into a material to detect imperfections. Results generated by all of these methods can be skewed due to surface conditions, and cannot easily isolate damage at joints and boundaries of the structure. Unless a general vicinity of a damage location is known prior to inspection, none of these methods are easily or practically applied to large structures which are already in place and operating. On the other hand, resonant inspection methods are not capable of determining the extent or location of damage, and is used only on a component rather than a assembled structure. None of the above NDT methods are easily or practically applied to large structures requiring a high degree of structural integrity.

[0007] Because of these shortfalls in existing NDT methods when inspecting relatively large structures, structural damage detection using changes in vibration characteristics has received much attention in recent years. Vibration based health monitoring for rotating machinery is a relatively mature technology, using a non-model based approach to provide a qualitative comparison of current data to historical data. However, this type of vibration based damage testing does not work for most structures. Rather, vibration based damage detection for structures is model based, comparing test data to analytical data from finite element models to detect the location(s) and extent of damage. Vibration based damage detection methods fall into three basic categories. The first of these is direct methods such as optimal matrix updating algorithms, which identify damage location and extent in a single iteration. Because of the single iteration, these methods are not accurate in detecting a large level of damage. The second category is iterative methods. The methodology has only been for updating modeling, which determines modified structural parameters iteratively by minimizing differences between model and test data. The third category includes control-based eigenstructure assignment methods, which have the similar limitation to that of the direct methods indicated above and are not accurate in detecting a large level of damage. None of these current vibration based methods have been incorporated into an iterative algorithm that can detect small to large levels of damage, and the vibration based approach for structures remains an immature technology area which is not readily available on a commercial basis.

[0008] An object of the invention is to solve at least the above problems and/or disadvantages and to provide at least the advantages described hereinafter.

[0009] In at least one aspect of the present concepts, a system for determining stiffness parameters of a structure includes a sensor arranged to measure vibrations of the structure and output vibration information and a stiffness parameter unit configured to receive the vibration information output by the sensor, to determine natural frequency data of the structure, and to determine the stiffness parameters of the structure using the natural frequency data, wherein the natural frequency data comprises at least one measured natural frequency, and wherein the stiffness parameter unit is configured to determine the stiffness parameters of the structure when a number of stiffness parameters to be determined is greater than a number of measured natural frequencies.

[0010] Another object of the invention is to provide a system and method for detecting structural damage based on changes in natural frequencies and/or mode shapes.

[0011] An advantage of the system and method as embodied and broadly described herein is that it can be applied to a large operating structure with a large number (thousands or more) of degrees of freedom.

[0012] Another advantage of the system and method as embodied and broadly described herein is that it can accurately detect the location(s) and extent of small to large levels of damage and is especially useful for detecting a large level of damage with severe mismatch between the eigenparameters of the damaged and undamaged structures.

[0013] Another advantage of the system and method as embodied and broadly described herein is that it can work with a limited number of measured vibration modes.

[0014] Another advantage of the system and method as embodied and broadly described herein is that it can use measurement at only a small number of locations compared to the degrees of freedom of the system. A modified eigenvector expansion method is used to deal with the incomplete eigenvector measurement problem arising from experimental measurement of a lesser number of degrees of freedom than that of the appropriate analytical model.

[0015] Another advantage of the system and method as embodied and broadly described herein is that it can be applied to structures with slight nonlinearities such as opening and closing cracks. The random shaker test or the random impact series method can be used to average out slight nonlinearities and extract linearized natural frequencies and/or mode shapes of a structure.

[0016] Another advantage of the system and method as embodied and broadly described herein is that it can handle structures with closely spaced vibration modes, where mode switching can occur in the damage detection process.

[0017] Another advantage of the system and method as embodied and broadly described herein is that it can handle different levels of measurement noise with estimation errors within the noise levels.

[0018] Another advantage of the system and method as embodied and broadly described herein is that the damage detection method and the vibration testing methods such as the random impact series method enables damage detection and assessment to be automated, thus improving the reliability/integrity of results.

[0019] Another advantage of the system and method as embodied and broadly described herein is that damage detection and assessment may be automated in the field so that structural health can be monitored at central location and useful service life may be optimized.

[0020] Another advantage of the system and method as embodied and broadly described herein is that the random impact series method enables the modal parameters such as natural frequencies and/or mode shapes to be measured for a large structure or a structure in the field when there are noise effects such as those arising from the wind or other ambient excitation.

[0021] Additional advantages, objects, and features of the invention will be set forth in part in the description which follows and in part will become apparent to those having ordinary skill in the art upon examination of the following or may be learned from practice of the invention. The objects and advantages of the invention may be realized and attained as particularly pointed out in the appended claims.

[0063] Commonly measured modal parameters, such as natural frequencies and mode shapes, are functions of physical properties of a particular structure. Therefore, changes in these physical properties, such as reductions in stiffness resulting from the onset of cracks or a loosening of a connection, will cause detectable changes in these modal parameters. Thus, if the changes in these parameters are indicators of damage, vibration based damage detection may be, simplistically, reduced to a system identification problem. However, a number of factors have made vibration based damage detection difficult to implement in practice in the past.

[0064] The system and method for detecting structural damage as embodied and broadly described herein is motivated by the observed advantages of vibration based damage detection over currently available technologies. It is well understood that this system and method may be effectively applied to damage detection and assessment for substantially all types and configurations of structures, including, but not limited to, simple beams, hollow tubes, trusses, frames, and the like. However, simply for ease of discussion, the system and method will first be discussed with respect to three examples—a mass-spring model, a beam, and a space frame—for conceptualization purposes. The system and method will later be applied lightning masts in electric substations.

[0065] FIG. 1A shows a system 100 for detecting structural damage according to one embodiment of the invention. System 100 includes a stiffness parameter unit 103 for detecting stiffness of a structure in question. Stiffness parameter unit 103 may include a spectrum analyzer or any other device capable of receiving vibration information and providing natural frequency and/or mode shape data. One example might be a spectrum analyzer capable of performing a fast fourier transform (FFT) and with modal analysis software for deriving mode shape data from vibration information received from sensor 110 if mode shape information is needed. A roving sensor technique or multiple sensors may need to be used if mode shape information needs to be used. If mode shape information is desired, a sensor should be provided at the tip of the hammer or device used to excite the structure. This sensor should be configured to send data to the spectrum analyzer in order to obtain the frequency response function, which the spectrum analyzer (modal analysis software) uses to obtain mode shape information. Stiffness parameter unit 103 is configured to receive vibration information output from the sensor 110 via sensor coupler 113 and input 114. Sensor coupler 113 could be a wire, optical fiber or even a wireless connection between sensor 110 and stiffness parameter unit 103.

[0066] Stiffness parameter unit 103 may include an iterative processing unit 115 capable of determining stiffness parameters using a first order perturbation approach and the generalized inverse method. The first order perturbation approach can include using natural frequencies and/or mode shapes of the structure according to a preferred embodiment of the invention. Iterative processing unit 115 may be capable of converging to a correct result for the output of stiffness parameters even when the system is underdetermined. Stiffness parameter unit 103 may further include an outer iterative processing unit 117 and an inner (nested) iterative processing unit 119 which may operate using a first or higher order perturbation approach and a gradient or quasi-Newton method to be discussed below. Stiffness parameter unit 103 outputs stiffness parameters to damage information processor 123. Damage information processor includes a damage extent processor 125 and a damage location processor 127. Damage location processor 127 outputs the location(s) of damage and damage extent processor 125 outputs the extent of damage.

[0067] System 100 using the iterative processing units 115 or 117 and 119 is capable of determining the stiffness parameters and ultimately the location(s) and extent of damage for structures with a largest dimension less than 50 meters, 25 meters, 10 meters, 2.5 meters and even 1.5 meters. Examples of this are provided herein.

[0068] The system and method may include a multiple-parameter, general order perturbation method, in which the changes in the stiffness parameters are used as the perturbation parameters. By equating the coefficients of like-order terms involving the same perturbation parameters in the normalization relations of eigenvectors and the eigenvalue problem, the perturbation problem solutions of all orders may be derived, and the sensitivities of all eigenparameters may be obtained. The perturbation method may be used in an iterative manner with an optimization method to identify the stiffness parameters of structures.

Methodology

[0069] This method presented below can simultaneously identify all the unknown stiffness parameters and is formulated as a damage detection problem. Since the effects of the changes in the inertial properties of a damaged structure are usually relatively small, only the changes in the stiffness properties due to structural damage are considered.

[0070] Consider a N-degree-of-freedom, linear, time-invariant, self-adjoint system with distinct eigenvalues. The stiffness parameters of the undamaged structure are denoted by G_hi (i=1, 2, . . . , m), where m is the number of the stiffness parameters. Structural damage is characterized by reductions in the stiffness parameters. The estimated stiffness parameters of the damaged structure before each iteration are denoted by G_i (i=1, 2, . . . , m), and its stiffness matrix, which depends linearly on G_i, is denoted by K=K(G), where G=[G₁, G₂, . . . , G_m]^T. Here the superscript T denotes matrix transpose. The eigenvalue problem of the structure with stiffness parameters G_i is
Kφ^k=λ^kMφ^k   (1)
where M is the constant mass matrix, and λ^k=λ^k(G) and φ^k=φ^k(G) (k=1, 2, . . . , N) are the k-th eigenvalue and mass-normalized eigenvector respectively. It is noted that λ_k=ω_k², where ω_k is the k-th natural frequency of the structure. The normalized eigenvectors of (1) satisfy the orthonormality relations
(φ^k)^TMφ^u=δ_ku(φ^k)^TKφ^u=λ^kδ_ku  (2)
where 1≦u≦N and δ_ku is the Kronecker delta. Before the first iteration, the initial stiffness parameters of the damaged structure are assumed to be G_i⁽⁰⁾=σ_iG_hi (i=1, 2, . . . , m), where 0<σ_i≦1, and the eigenvalue problem in (1) corresponds to that of the structure with stiffness parameters G_i⁽⁰⁾. If there is no prior knowledge of the integrity of the structure, one can start the iteration from the stiffness parameters of the undamaged structure and set σ_i=1. Let G_d1 (i=1, 2, . . . , m) denote the stiffness parameters of the damaged structure. The eigenvalue problem of the damaged structure is
K_dφ_d^k=λ_d^kMφ_d^k  (3)
where K_d=K(G_d) is the stiffness matrix with G_d=[G_d1, G_d2, . . . , G_dm]^T, and λ_d^k=λ^k(G_d) and φ_d^k=φ^k(G_d) are the k-th eigenvalue and mass-normalized eigenvector respectively. The stiffness matrix K_d is related to K through the Taylor expansion

[0071] $\begin{array}{cc}{K}_d=K\left(G_d\right)=K+\sum _{i=1}^m\frac{\partial K}{\partial G_i}\delta G_i& \left(4\right)\end{array}$
where δG_i=G_di−G_i (i=1, 2, . . . , m) are the changes in the stiffness parameters, and the higher-order derivatives of K with respect to G_i vanish because K is assumed to be a linear function of G_i. Based on the finite element model, the global stiffness matrix of a distributed structure satisfies (4) as its higher-order derivatives with respect to each element stiffness parameter vanish.

[0072] Let the k-th eigenvalue and mass-normalized eigenvector of the damaged structure be related to λ^k and φ^k through

[0073] $\begin{array}{cc}{\lambda }_d^k={\lambda }^k+\sum _{i=1}^m{\lambda }_{\left(1\right)i}^k\delta G_i+\sum _{i=1}^m\sum _{j=1}^m{\lambda }_{\left(2\right)\mathrm{ij}}^k\delta G_i\delta G_j+\dots +\underset{\underset{p\mathrm{summations}}{︸}}{\sum _{i=1}^m\sum _{j=1}^m\dots \sum _{t=1}^m}{\lambda }_{\left(p\right)\mathrm{ij}\dots t}^k\delta G_i\delta G_j\dots \delta G_t+e_{\lambda }^k,& \left(5\right)\\ {\varphi }_d^k={\varphi }^k+\sum _{i=1}^m{z}_{\left(1\right)i}^k\delta G_i+\sum _{i=1}^m\sum _{j=1}^m{z}_{\left(2\right)\mathrm{ij}}^k\delta G_i\delta G_j+\dots +\underset{\underset{p\mathrm{summations}}{︸}}{\sum _{i=1}^m\sum _{j=1}^m\dots \sum _{t=1}^m}{z}_{\left(p\right)\mathrm{ij}\dots t}^k\delta G_i\delta G_j\dots \delta G_t+e_{\varphi }^k& \left(6\right)\end{array}$
where λ_(1)i^k, λ_(2)ij^k, . . . , and λ_{(p)ij . . . t}^k are the coefficients of the first, second, . . . , and p-th order perturbations for the eigenvalue, z_(1)i^k, z_(2)ij^k, . . . , and z_{(p)ij . . . t}^k are the coefficient vectors of the first, second, . . . , and p-th order perturbations for the eigenvector, and e_λ^k and e_φ^k are the residuals of order p+1. Note that the numbers in the parentheses in the subscripts of the coefficients and coefficient vectors indicate the orders of the terms. By the Taylor expansion, one has for any p≧1,

[0074] $\begin{array}{cc}p!{\lambda }_{\left(p\right)\mathrm{ij}\dots t}^k=\frac{{\partial }^p{\lambda }^k}{\partial G_i\partial G_j\dots \partial G_t}p!z_{\left(p\right)\mathrm{ij}\dots t}^k=\frac{{\partial }^p{\varphi }^k}{\partial G_i\partial G_j\dots \partial G_t}& \left(7\right)\end{array}$
By (7), λ_{(p)ij . . . t}^k and z_{(p)ij . . . t}^k are symmetric in the p indices, i, j, . . . , and t. The right-hand sides of (7) are the p-th order sensitivities of the eigenvalues and eigenvectors with respect to the stiffness parameters.

[0075] Using the normalization relations of the eigenvectors, φ^k and φ_d^k, and symmetry of the coefficient vectors in (6), as indicated earlier, one obtains

[0076] $\begin{array}{cc}1={\left({\varphi }_d^k\right)}^{T}M{\varphi }_d^k={\left({\varphi }^k+\sum _{i=1}^m{z}_{\left(1\right)i}^k\delta G_i+\sum _{i=1}^m\sum _{j=1}^m{z}_{\left(2\right)\mathrm{ij}}^k\delta G_i\delta G_j+\dots +\underset{\underset{p\mathrm{summations}}{︸}}{\sum _{i=1}^m\sum _{j=1}^m\dots \sum _{t=1}^m}{z}_{\left(p\right)\mathrm{ij}\dots \mathrm{st}}\delta G_i\delta G_j\dots \delta G_s\delta G_t+\dots \right)}^{T}M\left({\varphi }^k+\sum _{i=1}^m{z}_{\left(1\right)i}^k\delta G_i+\sum _{i=1}^m\sum _{j=1}^m{z}_{\left(2\right)\mathrm{ij}}^k\delta G_i\delta G_j+\dots +\underset{\underset{p\mathrm{summations}}{︸}}{\sum _{i=1}^m\sum _{j=1}^m\dots \sum _{t=1}^m}{z}_{\left(p\right)\mathrm{ij}\dots \mathrm{st}}^k\delta G_i\delta G_j\dots \delta G_s{\delta }_t+\dots \right)=1+\sum _{i=1}^m\left[{\left({\varphi }^k\right)}^T{\mathrm{Mz}}_{\left(1\right)i}^k+{\left(z_{\left(1\right)i}^k\right)}^{T}M{\varphi }^k\right]\delta G_i+\sum _{i=1}^m\sum _{j=i}^m\frac{1}{R_{\mathrm{ij}}^1}\left\{2!{\left({\varphi }^k\right)}^T{\mathrm{Mz}}_{\left(2\right)\mathrm{ij}}^k+\left[{\left(z_{\left(1\right)i}^k\right)}^T{\mathrm{Mz}}_{\left(1\right)j}^k+{\left(z_{\left(1\right)i}^k\right)}^T{\mathrm{Mz}}_{\left(1\right)i}^k\right]+2!{\left(z_{\left(2\right)\mathrm{ij}}^k\right)}^{T}M{\varphi }^k\right\}\delta G_i\delta G_j+\sum _{i=1}^m\sum _{j=i}^m\sum _{l=j}^m\frac{1}{R_{\mathrm{ijl}}^1}\left\{3!{\left({\varphi }^k\right)}^T{\mathrm{Mz}}_{\left(3\right)\mathrm{ijl}}^k+2!\left[{\left(z_{\left(1\right)i}^k\right)}^T{\mathrm{Mz}}_{\left(2\right)\mathrm{jl}}^k+{\left(z_{\left(1\right)j}^k\right)}^T{\mathrm{Mz}}_{\left(2\right)\mathrm{il}}^k+{\left(z_{\left(1\right)l}^k\right)}^T{\mathrm{Mz}}_{\left(2\right)\mathrm{ij}}^k\right]+2!\left[{\left(z_{\left(2\right)\mathrm{jl}}^k\right)}^T{\mathrm{Mz}}_{\left(1\right)i}^k+{\left(z_{\left(2\right)\mathrm{il}}^k\right)}^T{\mathrm{Mz}}_{\left(1\right)j}^k+{\left(z_{\left(2\right)\mathrm{ij}}^k\right)}^T{\mathrm{Mz}}_{\left(1\right)l}^k\right]+3!{\left(z_{\left(3\right)\mathrm{ijl}}^k\right)}^{T}M{\varphi }^k\right\}\delta G_i\delta G_j\delta G_l+\dots +\underset{\underset{p\mathrm{summations}}{︸}}{\sum _{i=1}^m\sum _{j=i}^m\dots \sum _{t=s}^m}\frac{1}{R_{\mathrm{ij}\dots \mathrm{st}}^1}\{p!{\left({\varphi }^k\right)}^T{\mathrm{Mz}}_{\left(p\right)\mathrm{ij}\dots \mathrm{st}}^k+\left(p-1\right)!\hspace{1em}\left[{\left(z_{\left(1\right)i}^k\right)}^T{\mathrm{Mz}}_{\left(p-1\right)\mathrm{jl}\dots \mathrm{st}}^k+{\left(z_{\left(1\right)j}^k\right)}^T{\mathrm{Mz}}_{\left(p-1\right)\mathrm{il}\dots \mathrm{st}}^k+\dots +{\left(z_{\left(1\right)t}^k\right)}^T{\mathrm{Mz}}_{\left(p-1\right)\mathrm{ij}\dots s}^k\right]+\hspace{1em}2!\left(p-2\right)![{\left(z_{\left(2\right)\mathrm{ij}}^k\right)}^T{\mathrm{Mz}}_{\left(p-2\right)l\dots \mathrm{st}}^k+{\left(z_{\left(2\right)\mathrm{il}}^k\right)}^T{\mathrm{Mz}}_{\left(p-2\right)j\dots \mathrm{st}}^k+\dots +{\left(z_{\left(2\right)\mathrm{st}}^k\right)}^T{\mathrm{Mz}}_{\left(p-2\right)\mathrm{ij}\dots r}^k]+\dots +\left(p-2\right)!2![{\left(z_{\left(p-2\right)l\dots \mathrm{st}}^k\right)}^T{\mathrm{Mz}}_{\left(2\right)\mathrm{ij}}^k+{\left(z_{\left(p-2\right)j\dots \mathrm{st}}^k\right)}^T{\mathrm{Mz}}_{\left(2\right)\mathrm{il}}^k+\dots +{\left(z_{\left(p-2\right)\mathrm{ij}\dots r}^k\right)}^T{\mathrm{Mz}}_{\left(2\right)\mathrm{st}}^k]+\left(p-1\right)![{\left(z_{\left(p-1\right)\mathrm{jl}\dots \mathrm{st}}^k\right)}^T{\mathrm{Mz}}_{\left(1\right)i}^k+{\left(z_{\left(p-1\right)\mathrm{il}\dots \mathrm{st}}^k\right)}^T{\mathrm{Mz}}_{\left(1\right)j}^k+\dots +{\left(z_{\left(p-1\right)\mathrm{ij}\dots s}^k\right)}^T{\mathrm{Mz}}_{\left(1\right)t}^k]+p!{\left(z_{\left(p\right)\mathrm{ij}\dots \mathrm{st}}^k\right)}^{T}M{\varphi }^k\}\delta G_i\delta G_j\dots \delta G_s{\mathrm{\delta G}}_t+\dots & \left(8\right)\end{array}$
where the superscript T denotes transpose. For any p-th (p≧1) order term in the last expression of (8), the coefficient R_{ij . . . t}¹ is defined by R_{ij . . . t}¹=X₁!X₂! . . . X_a!, where X₁, X₂, . . . , and X_a (1≦a≦p) are the numbers of the first, second, . . . , and last distinct indices within indices, i, j, . . . , and t, with X₁+X₂++X_a=p. For instance, R₁₂₃₄=1!1!1!1!=1 with a=4, and R₁₁₂₂₂₃=2!3!1!=2!3! with a=3. Some explanations of the general p-th order term in the last expansion in (8) are in order. It includes p+1 types of terms ordered from the beginning to the end of the expression within the braces, with each type of terms except the first and last ones enclosed in square brackets. The c-th (1≦c≦p+1) type of terms is obtained by multiplying a (c−1)-th order term in the expansion of (φ_d^k)^T by a (p−c+1)-th order term in the expansion of Mφ_d^k. For the c-th type of term, where 2≦c≦p, the p indices, i, j, . . . , and t are distributed in the subscripts of the two vectors pre- and post-multiplying M, whose numbers of indices in the subscripts are c−1 and p−c+1 respectively. For the first and last (p+1)-th) types of terms, all the p indices lie in the subscripts of the vectors post- and pre-multiplying M, respectively. The number of terms within each set of square brackets equals the number of different combinations of indices in the subscripts of the vectors pre- and post-multiplying M. When all the p indices, i, j, . . . , and t, have distinct values, due to symmetry of these vectors in their indices, terms of the c-th (1≦c≦p+1) type, involving different permutations of the same indices in the subscripts of the vectors, are equal and combined, resulting in the scaling factor (c−1)!(p−c+1)! in front of the square brackets. Consequently, the indices in the second through p-th summations range from the values of their preceding summation indices to m. When any of the p indices, i, j, . . . , and t, have equal values, the corresponding terms in each type are given by those in the previous case divided by R_{ij . . . t}¹. For instance, the 4^th order term of the form δG₁δG₂³ in the expansion of (φ_d^k)^TMφ_d^k can be obtained from the expression for the p-th order term in (8):

[0077] $\begin{array}{cc}\frac{1}{1!3!}\left\{4!{\left({\varphi }^k\right)}^T{\mathrm{Mz}}_{\left(4\right)1222}^k+3!\left[{\left(z_{\left(1\right)1}^k\right)}^T{\mathrm{Mz}}_{\left(3\right)222}^k+{\left(z_{\left(1\right)2}^k\right)}^T{\mathrm{Mz}}_{\left(3\right)122}^k+{\left(z_{\left(1\right)2}^k\right)}^T{\mathrm{Mz}}_{\left(3\right)122}^k+{\left(z_{\left(1\right)2}^k\right)}^T{\mathrm{Mz}}_{\left(3\right)122}^k\right]+2!2!\left[{\left(z_{\left(2\right)12}^k\right)}^T{\mathrm{Mz}}_{\left(2\right)22}^k+{\left(z_{\left(2\right)12}^k\right)}^T{\mathrm{Mz}}_{\left(2\right)22}^k+{\left(z_{\left(2\right)12}^k\right)}^T{\mathrm{Mz}}_{\left(2\right)22}^k+{\left(z_{\left(2\right)22}^k\right)}^T{\mathrm{Mz}}_{\left(2\right)12}^k+{\left(z_{\left(2\right)22}^k\right)}^T{\mathrm{Mz}}_{\left(2\right)12}^k+{\left(z_{\left(2\right)22}^k\right)}^T{\mathrm{Mz}}_{\left(2\right)12}^k\right]+3!\left[{\left(z_{\left(3\right)222}^k\right)}^T{\mathrm{Mz}}_{\left(1\right)1}^k+{\left(z_{\left(3\right)122}^k\right)}^T{\mathrm{Mz}}_{\left(1\right)2}^k+{\left(z_{\left(3\right)122}^k\right)}^T{\mathrm{Mz}}_{\left(1\right)2}^k+{\left(z_{\left(3\right)122}^k\right)}^T{\mathrm{Mz}}_{\left(1\right)2}^k\right]+4!{\left(z_{\left(4\right)1222}^k\right)}^{T}M{\varphi }^k\right\}\hspace{1em}\delta G_1\delta G_2^3=\left\{4{\left({\varphi }^k\right)}^T{\mathrm{Mz}}_{\left(4\right)1222}^k+\left[{\left(z_{\left(1\right)1}^k\right)}^T{\mathrm{Mz}}_{\left(3\right)222}^k+3{\left(z_{\left(1\right)2}^k\right)}^T{\mathrm{Mz}}_{\left(3\right)122}^k\right]+4{\left(z_{\left(2\right)12}^k\right)}^T{\mathrm{Mz}}_{\left(2\right)22}^k+\left[{\left(z_{\left(3\right)222}^k\right)}^T{\mathrm{Mz}}_{\left(1\right)1}^k+3{\left(z_{\left(3\right)122}^k\right)}^T{\mathrm{Mz}}_{\left(1\right)2}^k\right]+4{\left(z_{\left(4\right)1222}^k\right)}^{T}M{\varphi }^k\right\}\delta G_1\delta G_2^3& \left(9\right)\end{array}$
where p=4 and the four indices involved, i, j, l, and o, are i=1 and j=l=o=2.

[0078] Equating the coefficients of the δG_i (i=1, 2, . . . , m)-terms in (8) and using symmetry of the mass matrix yields
(φ^k)^TMz_(1)i^k=0  (10)
Equating the coefficients of the δG_iδG_j terms and using symmetry of M and z_(2)ij^k yields

[0079] $\begin{array}{cc}{\left({\varphi }^k\right)}^T{\mathrm{Mz}}_{\left(2\right)\mathrm{ij}}^k=-\frac{1}{2!}{\left(z_{\left(1\right)i}^k\right)}^T{\mathrm{Mz}}_{\left(1\right)j}^k& \left(11\right)\end{array}$
for all i, j=1, 2, . . . , m. Following a similar procedure, one obtains

[0080] $\begin{array}{cc}{\left({\varphi }^k\right)}^T{\mathrm{Mz}}_{\left(3\right)\mathrm{ijl}}^k=-\frac{2!}{3!}\left[{\left(z_{\left(1\right)i}^k\right)}^T{\mathrm{Mz}}_{\left(2\right)\mathrm{jl}}^k+{\left(z_{\left(l\right)j}^k\right)}^T{\mathrm{Mz}}_{\left(2\right)\mathrm{li}}^k+{\left(z_{\left(1\right)l}^k\right)}^T{\mathrm{Mz}}_{\left(2\right)\mathrm{ij}}^k\right]& \left(12\right)\end{array}$
for i, j, l=1, 2, . . . , m. Equating the coefficients of the δG_iδG_j . . . δG_sδG_t terms of p-th order in (8) yield

[0081] $\begin{array}{cc}{\left({\varphi }^k\right)}^T{\mathrm{Mz}}_{\left(p\right)\mathrm{ij}\dots t}^k=-\frac{1}{2\left(p!\right)}\left\{\left(p-1\right)!\left[{\left(z_{\left(1\right)i}^k\right)}^T{\mathrm{Mz}}_{\left(p-1\right)\mathrm{jl}\dots \mathrm{st}}^k+{\left(z_{\left(1\right)j}^k\right)}^T{\mathrm{Mz}}_{\left(p-1\right)\mathrm{il}\dots \mathrm{st}}^k+\dots +{\left(z_{\left(1\right)t}^k\right)}^T{\mathrm{Mz}}_{\left(p-1\right)\mathrm{ij}\dots s}^k\right]+2!\left(p-2\right)!\left[{\left(z_{\left(2\right)\mathrm{ij}}^k\right)}^T{\mathrm{Mz}}_{\left(p-2\right)l\dots \mathrm{st}}^k+{\left(z_{\left(2\right)\mathrm{il}}^k\right)}^T{\mathrm{Mz}}_{\left(p-2\right)j\dots \mathrm{st}}^k+\dots +{\left(z_{\left(2\right)\mathrm{st}}^k\right)}^T{\mathrm{Mz}}_{\left(p-2\right)\mathrm{ij}\dots q}^k\right]+\dots +\left(p-2\right)!2!\left[{\left(z_{\left(p-2\right)l\dots \mathrm{st}}^k\right)}^T{\mathrm{Mz}}_{\left(2\right)\mathrm{ij}}^k+{\left(z_{\left(p-2\right)j\dots \mathrm{st}}^k\right)}^T{\mathrm{Mz}}_{\left(2\right)\mathrm{il}}^k+\dots +{\left(z_{\left(p-2\right)\mathrm{ij}\dots r}^k\right)}^T{\mathrm{Mz}}_{\left(2\right)\mathrm{st}}^k\right]+\left(p-1\right)!\left[{\left(z_{\left(p-1\right)\mathrm{jl}\dots \mathrm{st}}^k\right)}^T{\mathrm{Mz}}_{\left(1\right)i}^k+{\left(z_{\left(p-1\right)\mathrm{il}\dots \mathrm{st}}^k\right)}^T{\mathrm{Mz}}_{\left(1\right)j}^k+\dots +{\left(z_{\left(p-1\right)\mathrm{ij}\dots s}^k\right)}^T{\mathrm{Mz}}_{\left(1\right)t}^k\right]\right\}& \left(13\right)\end{array}$
for i, j, . . . , t=1, 2, . . . , m. The p−1 types of terms, enclosed in the p−1 sets of square brackets on the right-hand side of (13), are ordered from the beginning to the end of the expression within the braces, and their structures are readily observed. When p is odd, by symmetry of M and z_{(p)ij . . . t}^k, the

[0082] $y\text{-}\mathrm{th}\left(1\le y\le \frac{p-1}{2}\right)$
type of terms is identical to the (p−y)-th type, and the two types of terms can be combined. Similarly, when p is even, the

[0083] $y\text{-}\mathrm{th}\left(1\le y\le \frac{p}{2}-1\right)$
type of terms equals and can be combined with the (p−y)-th type of terms. In this case, however, there is a separate type, the

[0084] $\left(\frac{p}{2}\right)\text{-}\mathrm{th}$
type, of terms in the middle of the expression.

[0085] Substituting (4)-(6) into (3) yields

[0086] $\begin{array}{cc}\left\{K+\sum _{i=1}^m\frac{\partial K}{\partial G_i}\delta G_i\right\}\left\{{\varphi }^k+\sum _{i=1}^m{z}_{\left(1\right)i}^k\delta G_i+\sum _{i=1}^m\sum _{j=1}^m{z}_{\left(2\right)\mathrm{ij}}^k\delta G_i\delta G_j+\sum _{i=1}^m\sum _{j=1}^m\sum _{l=1}^m{z}_{\left(3\right)\mathrm{ijl}}^k\delta G_i\delta G_j\delta G_l+\dots \right\}=\left\{{\lambda }^k+\sum _{i=1}^m{\lambda }_{\left(1\right)i}^k\delta G_i+\sum _{i=1}^m\sum _{j=1}^m{\lambda }_{\left(2\right)\mathrm{ij}}^k\delta G_i\delta G_j+\sum _{i=1}^m\sum _{j=1}^m\sum _{l=1}^m{\lambda }_{\left(3\right)\mathrm{ijl}}^k\delta G_i\delta G_j\delta G_l+\dots \right\}M\left\{{\varphi }^k+\sum _{i=1}^m{z}_{\left(1\right)i}^k\delta G_i+\sum _{i=1}^m\sum _{j=1}^m{z}_{\left(2\right)\mathrm{ij}}^k\delta G_i\delta G_j+\sum _{i=1}^m\sum _{j=1}^m\sum _{l=1}^m{z}_{\left(3\right)\mathrm{ijl}}^k\delta G_i\delta G_j\delta G_l+\dots \right\}& \left(14\right)\end{array}$
Equating the coefficients of the δG_i (i=1, 2, . . . , m) terms in (14) yields

[0087] $\begin{array}{cc}{\mathrm{Kz}}_{\left(1\right)i}^k+\frac{\partial K}{\partial G_i}{\varphi }^k={\lambda }^k{\mathrm{Mz}}_{\left(1\right)i}^k+{\lambda }_{\left(1\right)i}^{k}M{\varphi }^k& \left(15\right)\end{array}$
Expanding z_(1)i^k using normalized eigenvectors of (1) as basis vectors, one has

[0088] $\begin{array}{cc}{z}_{\left(1\right)i}^k=\sum _{u=1}^N{P}_{\left(1\right)\mathrm{iu}}^k{\varphi }^u& \left(16\right)\end{array}$
where P_(1)iu^k is the coefficient of the u-th term and the number in the parentheses in its subscript corresponds to that of z_(1)i^k. Pre-multiplying (15) by (φ^k)^T and using (1), (2), and (16) yields

[0089] $\begin{array}{cc}{\lambda }_{\left(1\right)i}^k={\left({\varphi }^k\right)}^T\frac{\partial K}{\partial G_i}{\varphi }^k& \left(17\right)\end{array}$
Substituting (16) into (10) and using (2) yields
P_(1)ik^k=0  (18)
Pre-multiplying (15) by (φ^ν)^T, where 1≦ν≦N and ν≠k, and using (1), (2), and (16) yields

[0090] $\begin{array}{cc}{P}_{\left(1\right)\mathrm{iv}}^k=\frac{1}{{\lambda }^k-{\lambda }^v}{\left({\varphi }^v\right)}^T\frac{\partial K}{\partial G_i}{\varphi }^k& \left(19\right)\end{array}$
By (16), (18) and (19) we have determined z_(1)i^k.

[0091] Equating the coefficients of the δG_iδG_j terms in (14) yields

[0092] $\begin{array}{cc}2!{\mathrm{Kz}}_{\left(2\right)\mathrm{ij}}^k+\frac{\partial K}{\partial G_i}{z}_{\left(1\right)j}^k+\frac{\partial K}{\partial G_j}{z}_{\left(1\right)i}^k=2!{\lambda }^k{\mathrm{Mz}}_{\left(2\right)\mathrm{ij}}^k+{\lambda }_{\left(1\right)i}^k{\mathrm{Mz}}_{\left(1\right)j}^k+{\lambda }_{\left(1\right)j}^k{\mathrm{Mz}}_{\left(1\right)i}^k+2!{\lambda }_{\left(2\right)\mathrm{ij}}^{k}M{\varphi }^k& \left(20\right)\end{array}$
Expanding z_(2)ij^k using normalized eigenvectors of (1) as basis vectors, one has

[0093] $\begin{array}{cc}{z}_{\left(2\right)\mathrm{ij}}^k=\sum _{u=1}^N{P}_{\left(2\right)\mathrm{iju}}^k{\varphi }^u& \left(21\right)\end{array}$
where P_(2)iju^k is the coefficient of the u-th term and the number in the parentheses in its subscript corresponds to that of z_(2)ij^k. Pre-multiplying (20) by (φ^k)^T and using (1), (2), (10), and (21) yields

[0094] $\begin{array}{cc}{\lambda }_{\left(2\right)\mathrm{ij}}^k=\frac{1}{2!}{\left({\varphi }^k\right)}^T\left[\frac{\partial K}{\partial G_i}{z}_{\left(1\right)j}^k+\frac{\partial K}{\partial G_j}{z}_{\left(1\right)i}^k\right]& \left(22\right)\end{array}$
Substituting (21) into (11) and using (2) yields

[0095] $\begin{array}{cc}{P}_{\left(2\right)\mathrm{ijk}}^k=-\frac{1}{2!}{\left(z_{\left(1\right)i}^k\right)}^T{\mathrm{Mz}}_{\left(1\right)j}^k& \left(23\right)\end{array}$
Pre-multiplying (20) by (φ^ν)^T, where 1≦ν≦N and ν≠k, and using (1), (2), and (21) yields

[0096] $\begin{array}{cc}{P}_{\left(2\right)\mathrm{ijv}}^k=\frac{1}{2!\left({\lambda }^k-{\lambda }^v\right)}{\left({\varphi }^v\right)}^T\left\{\left[\frac{\partial K}{\partial G_i}{z}_{\left(1\right)j}^k+\frac{\partial K}{\partial G_j}{z}_{\left(1\right)i}^k\right]-\left[{\lambda }_{\left(1\right)i}^k{\mathrm{Mz}}_{\left(1\right)j}^k+{\lambda }_{\left(1\right)j}^k{\mathrm{Mz}}_{\left(1\right)i}^k\right]\right\}& \left(24\right)\end{array}$
By (21), (23) and (24) we have determined z_(2)ij^k.

[0097] Equating the coefficients of the δG_iδG_jδG_t terms in (14) yields

[0098] $\begin{array}{cc}3!{\mathrm{Kz}}_{\left(3\right)\mathrm{ijl}}^k+2!\left[\frac{\partial K}{\partial G_i}{z}_{\left(2\right)\mathrm{jl}}^k+\frac{\partial K}{\partial G_j}{z}_{\left(2\right)\mathrm{li}}^k+\frac{\partial K}{\partial G_l}{z}_{\left(2\right)\mathrm{ij}}^k\right]=3!{\lambda }^k{\mathrm{Mz}}_{\left(3\right)\mathrm{ijl}}^k+2!\left[{\lambda }_{\left(1\right)i}^k{\mathrm{Mz}}_{\left(2\right)\mathrm{jl}}^k+{\lambda }_{\left(1\right)j}^k{\mathrm{Mz}}_{\left(2\right)\mathrm{li}}^k+{\lambda }_{\left(1\right)l}^k{\mathrm{Mz}}_{\left(2\right)\mathrm{ij}}^k\right]+2!\left[{\lambda }_{\left(2\right)\mathrm{jl}}^k{\mathrm{Mz}}_{\left(1\right)i}^k+{\lambda }_{\left(2\right)\mathrm{li}}^k{\mathrm{Mz}}_{\left(1\right)j}^k+{\lambda }_{\left(2\right)\mathrm{ij}}^k{\mathrm{Mz}}_{\left(1\right)l}^k\right]+3!{\lambda }_{\left(3\right)\mathrm{ijl}}^{k}M{\varphi }^k& \left(25\right)\end{array}$
Expanding z_(3)ijl^k using normalized eigenvectors of (1) as basis vectors, one has

[0099] $\begin{array}{cc}{z}_{\left(3\right)\mathrm{ijl}}^k=\sum _{u=1}^n{P}_{\left(3\right)\mathrm{ijlu}}^k{\varphi }^u& \left(26\right)\end{array}$
where P_(3)ijlu^k is the coefficient of the u-th term and the number in the parentheses in its subscript corresponds to that of z_(3)ijl^k. Pre-multiplying (25) by (φ^k)^T and using (1), (2), (10), and (26) yields

[0100] $\begin{array}{cc}{\lambda }_{\left(3\right)\mathrm{ijl}}^k=\frac{2!}{3!}{\left({\varphi }^k\right)}^T\left[\frac{\partial K}{\partial G_i}{z}_{\left(2\right)\mathrm{jl}}^k+\frac{\partial K}{\partial G_j}{z}_{\left(2\right)\mathrm{li}}^k+\frac{\partial K}{\partial G_i}{z}_{\left(2\right)\mathrm{ij}}^k-{\lambda }_{\left(1\right)i}^k{\mathrm{Mz}}_{\left(2\right)\mathrm{jl}}^k-{\lambda }_{\left(1\right)j}^k{\mathrm{Mz}}_{\left(2\right)\mathrm{li}}^k-{\lambda }_{\left(1\right)l}^k{\mathrm{Mz}}_{\left(2\right)\mathrm{ij}}^k\right]& \left(27\right)\end{array}$
Substituting (26) into (12) and using (2) yields

[0101] $\begin{array}{cc}{P}_{\left(3\right)\mathrm{ijlk}}^k=-\frac{2!}{3!}\left[{\left(z_{\left(1\right)i}^k\right)}^T{\mathrm{Mz}}_{\left(2\right)\mathrm{jl}}^k+{\left(z_{\left(1\right)j}^k\right)}^T{\mathrm{Mz}}_{\left(2\right)\mathrm{li}}^k+{\left(z_{\left(1\right)t}^k\right)}^T{\mathrm{Mz}}_{\left(2\right)\mathrm{ij}}^k\right]& \left(28\right)\end{array}$
Pre-multiplying (25) by (φ^ν)^T, where 1≦ν≦N and ν≠k, and using (1), (2), and (26) yields

[0102] $\begin{array}{cc}{P}_{\left(3\right)\mathrm{ijlv}}^k=\frac{2!}{3!\left({\lambda }^k-{\lambda }^v\right)}{\left({\varphi }^v\right)}^T\hspace{1em}[\frac{\partial K}{\partial G_i}{z}_{\left(2\right)\mathrm{jl}}^k+\frac{\partial K}{\partial G_j}{z}_{\left(2\right)\mathrm{li}}^k+\frac{\partial K}{\partial G_i}{z}_{\left(2\right)\mathrm{ij}}^k-{\lambda }_{\left(1\right)i}^k{\mathrm{Mz}}_{\left(2\right)\mathrm{jl}}^k-{\lambda }_{\left(1\right)j}^k{\mathrm{Mz}}_{\left(2\right)\mathrm{li}}^k-{\lambda }_{\left(1\right)l}^k{\mathrm{Mz}}_{\left(2\right)\mathrm{ij}}^k-{\lambda }_{\left(2\right)\mathrm{jl}}^k{\mathrm{Mz}}_{\left(1\right)i}^k-{\lambda }_{\left(2\right)\mathrm{li}}^k{\mathrm{Mz}}_{\left(1\right)j}^k-{\lambda }_{\left(2\right)\mathrm{ij}}^k{\mathrm{Mz}}_{\left(1\right)l}^k]& \left(29\right)\end{array}$
By (26), (28) and (29) we have determined z_(3)ijl^k.

[0103] We proceed now to derive the perturbation solutions for the general p-th order terms in (5) and (6). Equating the coefficients of the δG_iδG_j . . . δG_sδG_t terms of order p in (14) yields

[0104] $\begin{array}{cc}p!{\mathrm{Kz}}_{\left(p\right)\mathrm{ij}\dots \mathrm{st}}^k+\left(p-1\right)!\hspace{1em}\left[\frac{\partial K}{\partial G_i}{z}_{\left(p-1\right)j\dots \mathrm{st}}^k+\frac{\partial K}{\partial G_j}{z}_{\left(p-1\right)i\dots \mathrm{st}}^k+\dots +\frac{\partial K}{\partial G_t}{z}_{\left(p-1\right)\mathrm{ij}\dots s}^k\right]=p!{\lambda }^k{\mathrm{Mz}}_{\left(p\right)\mathrm{ij}\dots \mathrm{st}}^k+\left(p-1\right)![{\lambda }_{\left(1\right)i}^k{\mathrm{Mz}}_{\left(p-1\right)j\dots \mathrm{st}}^k+{\lambda }_{\left(1\right)j}^k{\mathrm{Mz}}_{\left(p-1\right)i\dots \mathrm{st}}^k+\dots +{\lambda }_{\left(1\right)t}^k{\mathrm{Mz}}_{\left(p-1\right)\mathrm{ij}\dots s}^k]+2!\left(p-2\right)!\left[{\lambda }_{\left(2\right)\mathrm{ij}}^k{\mathrm{Mz}}_{\left(p-2\right)l\dots \mathrm{st}}^k+{\lambda }_{\left(2\right)\mathrm{il}}^k{\mathrm{Mz}}_{\left(p-2\right)j\dots \mathrm{st}}^k+\dots +{\lambda }_{\left(2\right)\mathrm{st}}^k{\mathrm{Mz}}_{\left(p-2\right)\mathrm{ij}\dots q}^k\right]+\dots +p!{\lambda }_{\left(p\right)\mathrm{ij}\dots \mathrm{st}}^{k}M{\varphi }^k& \left(30\right)\end{array}$
Expanding z_{(p)ij . . . st}^k using normalized eigenvectors of (1) as basis vectors, one has

[0105] $\begin{array}{cc}{z}_{\left(p\right)\mathrm{ij}\dots \mathrm{st}}^k=\sum _{u=1}^n{P}_{\left(p\right)\mathrm{ij}\dots \mathrm{stu}}^k{\varphi }^u& \left(31\right)\end{array}$
where P_{(p)ij . . . stu}^k is the coefficient of the u-th term and the number in the parenthesis in its subscript corresponds to that of z_{(p)ij . . . st}^k. Pre-multiplying (30) by (φ^k)^T and using (1), (10), (31) and orthonormality relations of eigenvectors yields

[0106] $\begin{array}{cc}{\lambda }_{\left(p\right)\mathrm{ij}\dots \mathrm{st}}^k=\frac{1}{p!}{\left({\varphi }^k\right)}^T\left[\left(p-1\right)!\left(\frac{\partial K}{\partial G_i}{z}_{\left(p-1\right)j\dots \mathrm{st}}^k+\frac{\partial K}{\partial G_j}{z}_{\left(p-1\right)i\dots \mathrm{st}}^k+\dots +\frac{\partial K}{\partial G_t}{z}_{\left(p-1\right)\mathrm{ij}\dots s}^k\right)-\left(p-1\right)!\left({\lambda }_{\left(1\right)i}^k{\mathrm{Mz}}_{\left(p-1\right)j\dots \mathrm{st}}^k+{\lambda }_{\left(1\right)j}^k{\mathrm{Mz}}_{\left(p-1\right)i\dots \mathrm{st}}^k+\dots +{\lambda }_{\left(1\right)t}^k{\mathrm{Mz}}_{\left(p-1\right)\mathrm{ij}\dots s}^k\right)-2!\left(p-2\right)!\left({\lambda }_{\left(2\right)\mathrm{ij}}^k{\mathrm{Mz}}_{\left(p-2\right)l\dots \mathrm{st}}^k+{\lambda }_{\left(2\right)\mathrm{il}}^k{\mathrm{Mz}}_{\left(p-2\right)j\dots \mathrm{st}}^k+\dots +{\lambda }_{\left(2\right)\mathrm{st}}^k{\mathrm{Mz}}_{\left(p-1\right)\mathrm{ij}\dots r}^k\right)-\dots -\left(p-2\right)!2!\left({\lambda }_{\left(p-2\right)l\dots \mathrm{st}}^k{\mathrm{Mz}}_{\left(2\right)\mathrm{ij}}^k+{\lambda }_{\left(p-2\right)j\dots \mathrm{st}}^k{\mathrm{Mz}}_{\left(2\right)\mathrm{il}}^k+\dots +{\lambda }_{\left(p-2\right)\mathrm{ij}\dots r}^k{\mathrm{Mz}}_{2\left(\mathrm{st}\right)}^k\right)\right]& \left(32\right)\end{array}$
The p-th order sensitivities of eigenvalues are obtained from (7) and (32). They depend on the eigenvalue and eigenvector sensitivities of orders up to p−2 and p−1 respectively. Substituting (31) into (13) and using (2) yields