U.S. patent application number 10/849571 was filed with the patent office on 2005-04-07 for system and method for detecting structural damage.
Invention is credited to Emory, Benjamin Haynes, Wong, Chun Nam, Xu, Guangyao, Zheng, Nengan, Zhu, Weidong.
Application Number | 20050072234 10/849571 |
Document ID | / |
Family ID | 34396984 |
Filed Date | 2005-04-07 |
United States Patent
Application |
20050072234 |
Kind Code |
A1 |
Zhu, Weidong ; et
al. |
April 7, 2005 |
System and method for detecting structural damage
Abstract
A system and method for detecting structural damage is provided
that utilizes a general order perturbation methodology involving
multiple perturbation parameters. The perturbation methodology is
used iteratively in conjunction with an optimization method to
identify the stiffness parameters of structures using natural
frequencies and/or mode shape information. The stiffness parameters
are then used to determine the location and extent of damage in a
structure. A novel stochastic model is developed to model the
random impact series produced manually or to generate a random
impact series in a random impact device. The random impact series
method or the random impact device can be used to excite a
structure and generate vibration information used to obtain the
stiffness parameters of the structure. The method or the device can
also just be used for modal testing purposes. The random impact
device is a high energy, random, and high signal-to-noise ratio
system.
Inventors: |
Zhu, Weidong; (Elkridge,
MD) ; Xu, Guangyao; (Baltimore, MD) ; Zheng,
Nengan; (Baltimore, MD) ; Emory, Benjamin Haynes;
(Baltimore, MD) ; Wong, Chun Nam; (Lubbock,
TX) |
Correspondence
Address: |
FLESHNER & KIM, LLP
P.O. Box 221200
Chantilly
VA
20153-1200
US
|
Family ID: |
34396984 |
Appl. No.: |
10/849571 |
Filed: |
May 20, 2004 |
Related U.S. Patent Documents
|
|
|
|
|
|
Application
Number |
Filing Date |
Patent Number |
|
|
60471873 |
May 20, 2003 |
|
|
|
60512656 |
Oct 20, 2003 |
|
|
|
Current U.S.
Class: |
73/579 ; 702/56;
73/12.01 |
Current CPC
Class: |
G01N 29/227 20130101;
G01M 7/08 20130101; G01H 5/00 20130101; G01H 13/00 20130101; G01N
2291/0258 20130101; G01N 29/12 20130101; G01N 2291/0422 20130101;
G01H 1/00 20130101 |
Class at
Publication: |
073/579 ;
073/012.01; 702/056 |
International
Class: |
G01H 001/14; G01M
007/08 |
Claims
What is claimed is:
1. A system for determining stiffness parameters of a structure,
comprising: a sensor arranged to measure vibrations of said
structure and output vibration information; and a stiffness
parameter unit for receiving said vibration information,
determining natural frequency data of said structure, and
determining the stiffness parameters of said structure using said
natural frequency data.
2. The system according to claim 1, further comprising multiple
sensors arranged to measure vibrations of said structure and output
vibration information.
3. The system according to claim 1, wherein said stiffness
parameter unit comprises an iterative processing unit.
4. The system according to claim 1, wherein said stiffness
parameter unit comprises an outer iterative processing unit and an
inner iterative processing unit.
5. The system according to claim 3, wherein said iterative
processing unit determines said stiffness parameters using a first
order perturbation process.
6. The system according to claim 3, wherein said iterative
processing unit determines said stiffness parameters using a higher
order perturbation process.
7. A system for determining stiffness parameters of a structure,
comprising: a sensor arranged to measure vibrations of said
structure and output vibration information; and a stiffness
parameter unit for receiving said vibration information and
determining said stiffness parameters with an iterative processing
unit.
8. The system according to claim 7, wherein said iterative
processing unit comprises an outer iterative processing unit and an
inner iterative processing unit.
9. The system according to claim 7, wherein said iterative
processing unit determines said stiffness parameters using a first
order perturbation process.
10. The system according to claim 7, wherein said iterative
processing unit determines said stiffness parameters using a higher
order perturbation process.
11. A stiffness parameter unit for determining stiffness parameters
for a structure, comprising: an input for receiving vibration data
related to the structure; an analyzer for converting said vibration
data to spectral data; and an interative processing unit for
receiving said spectral data and outputting said stiffness
parameters using natural frequencies of the structure.
12. A stiffness parameter unit for determining stiffness parameters
for a structure, comprising: an input for receiving vibration data
related to the structure; an analyzer for converting said vibration
data to spectral data; and an interative processing unit for
receiving said spectral data and outputting said stiffness
parameters using a perturbation process.
13. The stiffness parameter unit according to claim 12, wherein
said perturbation process comprises a first order perturbation
process.
14. The stiffness parameter unit according to claim 12, wherein
said perturbation process comprises a higher order perturbation
process.
15. A system for determining damage information of a structure,
comprising: a sensor arranged to measure vibrations of said
structure and output vibration information; a stiffness parameter
unit for receiving said vibration information, determining natural
frequency data of said structure, and determining the stiffness
parameters of said structure using said natural frequency data; and
a damage information processor for receiving said stiffness
parameters and outputting damage information.
16. The system according to claim 15, wherein said damage
information processor outputs damage location information or extent
of damage information.
17. A system, comprising: a structure; a sensor arranged to measure
vibrations of said structure and output vibration information; and
a stiffness parameter unit for receiving said vibration
information, determining natural frequency data of said structure,
and determining the stiffness parameters of said structure using
said natural frequency data.
18. The system according to claim 17, further comprising a damage
information processor for receiving said stiffness parameters and
outputting location of damage.
19. The system according to claim 18, wherein said damage
information processor comprises a damage location processor for
determining damage location information.
20. The system according to claim 18, wherein said damage
information processor comprises a damage extent processor for
determining extent of damage information.
21. The system according to claim 18, wherein said damage
information processor comprises a damage extent processor for
determining extent of damage information and a damage location
processor for determining damage location information.
22. The system according to claim 17, wherein said sensor comprises
a velocimeter.
23. The system according to claim 17, wherein said sensor is
attached to said structure.
24. The system according to claim 17, wherein said sensor is not
attached to said structure.
25. The system according to claim 17, wherein said stiffness
parameter unit further comprises a spectral analyzer.
26. The system according to claim 17, wherein said structure
comprises a beam.
27. The system according to claim 17, wherein said structure
comprises a truss.
28. The system according to claim 17, wherein said structure has a
longest dimension less than 1.5 meters.
29. The system according to claim 17, wherein said structure has a
longest dimension less than 2.5 meters.
30. The system according to claim 17, wherein said structure has a
longest dimension less than 10 meters.
31. The system according to claim 17, wherein said structure has a
longest dimension less than 50 meters.
32. A device, comprising: a random signal generating unit for
generating first and second outputs; a random impact actuator for
receiving said first and second outputs; and an impact applicator
coupled to said random impact actuator and having an impact region;
wherein said random impact actuator drives said impact applicator
such that the force and arrival times of said impact applicator at
said impact region are random.
33. The device of claim 32, wherein said random impact actuator
drives said impact applicator in accordance with said first and
second outputs.
34. The device of claim 33, wherein the first and second outputs
comprise independent random variables.
35. The device of claim 34, wherein the first and second outputs
determine the force and arrival times, respectively, of the impact
applicator at the impact region.
36. A system, comprising: a structure; a random impact device for
inducing vibrations in said structure; a sensor arranged to measure
vibrations of said structure and output vibration information; and
a stiffness parameter unit for receiving said vibration
information, determining natural frequency data of said structure,
and determining the stiffness parameters of said structure using
said natural frequency data.
37. The system of claim 36, wherein the random impact device
comprises: a random signal generating unit for generating first and
second outputs; a random impact actuator for receiving said first
and second outputs; and an impact applicator coupled to said random
impact actuator and having an impact region; wherein said random
impact actuator drives said impact applicator such that the force
and arrival times of said impact applicator at said impact region
are random.
38. The device of claim 37, wherein said random impact actuator
drives said impact applicator in accordance with said first and
second outputs.
39. The device of claim 38, wherein the first and second outputs
comprise independent random variables.
40. The device of claim 39, wherein the first and second outputs
determine the force and arrival times, respectively, of the impact
applicator at the impact region.
41. A system for determining stiffness parameters of a structure,
comprising: a sensor arranged to measure vibrations of said
structure and output vibration information; and a stiffness
parameter unit for receiving said vibration information,
determining mode shape information, and determining the stiffness
parameters of said structure using said mode shape information.
42. The system according to claim 41, further comprising multiple
sensors arranged to measure vibrations of said structure and output
vibration information.
43. The system according to claim 41, wherein said stiffness
parameter unit comprises an iterative processing unit.
44. The system according to claim 41, wherein said stiffness
parameter unit comprises an outer iterative processing unit and an
inner iterative processing unit.
45. The system according to claim 43, wherein said iterative
processing unit determines said stiffness parameters using a first
order perturbation process.
46. The system according to claim 43, wherein said iterative
processing unit determines said stiffness parameters using a higher
order perturbation process.
Description
[0001] This application claims priority to U.S. Provisional Patent
Application Ser. No. 60/471,873 filed May 20, 2003 and U.S.
Provisional Patent Application Ser. No. 60/512,656 filed Oct. 10,
2003, which are both incorporated herein by reference.
BACKGROUND OF THE INVENTION
[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.
SUMMARY OF THE INVENTION
[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] 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.
[0010] 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.
[0011] 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.
[0012] 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.
[0013] 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.
[0014] 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.
[0015] 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.
[0016] 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.
[0017] 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.
[0018] 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.
[0019] 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.
[0020] 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.
BRIEF DESCRIPTION OF THE DRAWINGS
[0021] The invention will be described in detail with reference to
the following drawings in which like reference numerals refer to
like elements wherein:
[0022] FIG. 1A shows a system 100 for detection of structural
damage according to one embodiment of the invention;
[0023] FIG. 1B is a flowchart of an inverse algorithm for
identifying stiffness parameters of a damaged structure from a
select set of measured eigenparameters, in accordance with an
embodiment of the invention;
[0024] FIG. 2 is a flowchart of a quasi-Newton method for finding
an optimal solution to the system equations shown in the flowchart
of FIG. 1B, in accordance with an embodiment of the invention;
[0025] FIG. 3 is a schematic view of a serial mass-spring
system;
[0026] FIGS. 4A-4C illustrate estimation errors in a series of
iterations for a low order system with a small level of damage;
[0027] FIGS. 5A-5C illustrate estimation errors in a series of
iterations for a low order system with a large level of damage;
[0028] FIGS. 6A-6C illustrate estimation errors in a series of
iterations for a large order system with a large level of
damage;
[0029] FIG. 7 illustrates a finite element model of a fixed-fixed
beam;
[0030] FIGS. 8A-8C illustrate the estimated stiffness parameters
with complete eignenvector measurements and different noise levels
for a ten element beam with a large level of damage;
[0031] FIGS. 9A-9B illustrate stiffness parameters with reduced
eigenvector measurements and different noise levels for a ten
element beam with a large level of damage;
[0032] FIG. 10 illustrates a modular, four bay space frame;
[0033] FIG. 11 illustrates the estimated dimensionless stiffnesses
for the damaged frame as shown in FIG. 11;
[0034] FIG. 12 illustrates a cantilever aluminum beam test specimen
with uniform damage in approximately 5 elements;
[0035] FIG. 13 illustrates the experimental results for the
estimated bending stiffnesses of all the elements of a cantilever
aluminum beam test specimen with the actual damage between 10 cm
and 15 cm from the cantilevered end using a 40-element finite
element model;
[0036] FIG. 14 illustrates the experimental results for the
estimated bending stiffnesses of all the elements of a cantilevered
aluminum beam test specimen with the actual damage between 25 cm
and 30 cm from the cantilevered end using a 40-element finite
element model;
[0037] FIG. 15 illustrates the experimental results for the
estimated bending stiffnesses of all the elements of an undamaged
cantilever aluminum beam test specimen using a 40-element finite
element model;
[0038] FIG. 16 illustrates a cantilever aluminum beam test specimen
with a narrow cut;
[0039] FIG. 17 illustrates the experimental results for the
estimated bending stiffnesses of all the elements of the test
specimen as shown in FIG. 16 using a 45-element finite element
model;
[0040] FIG. 18 illustrates the estimated bending stiffnesses of all
the elements of the cantilever aluminum with simulated damage
between 9 cm and 15.75 cm from the cantilevered end using a
40-element finite element model;
[0041] FIG. 19 illustrates the estimated bending stiffnesses of all
the elements of the cantilever aluminum beam with simulated damage
between 29.25 cm and 36 cm from the cantilevered end using a
40-element finite element model;
[0042] FIG. 20 illustrates the estimated bending stiffnesses of all
the elements of the cantilever aluminum beam with multiple
simulated damage--one at the 3.sup.rd element and the other at the
20.sup.th element;
[0043] FIG. 21 illustrates the estimated bending stiffnesses of all
the elements of the cantilever aluminum beam with simulated uniform
damage from the 16.sup.th to the 18.sup.th element;
[0044] FIG. 22 illustrates a 50-foot lightning mast test specimen
with an eccentric spike at the top;
[0045] FIG. 23 illustrates the experimental results for the
estimated bending stiffnesses of all the elements of the lightning
mast as shown in FIG. 22;
[0046] FIG. 24 illustrates the estimated bending stiffnesses of all
the elements of the lightning mast as shown in FIG. 22 with
simulated damage between 0.76 m and 1.125 m from the ground using a
40-element finite element model;
[0047] FIG. 25 illustrates the estimated bending stiffnesses of all
the elements of the lightning mast as shown in FIG. 42 with
simulated damage between 6.89 m and 7.35 m from the ground using a
40-element finite element model;
[0048] FIGS. 26A-26C illustrate the estimated bending stiffnesses
of all the elements of a cantilever aluminum beam using the first
regularization method;
[0049] FIGS. 27A-27C illustrate the estimated bending stiffnesses
of all the elements of a cantilever aluminum beam using the second
regularization method;
[0050] FIG. 28 illustrates an example of a practical application of
the damage detection method of FIG. 1A, in accordance with an
embodiment of the invention;
[0051] FIG. 29 illustrates the application of a series of random
impacts to the practical application of the damage detection method
of FIG. 1A, in accordance with an embodiment of the invention;
[0052] FIG. 30 illustrates a system for applying the series of
random impacts shown in FIG. 29 to the practical application of the
damage detection method of FIG. 1A, in accordance with an
embodiment of the invention;
[0053] FIG. 31 is a block diagram of one preferred embodiment of
the random impact device of FIG. 30;
[0054] FIG. 32 illustrates a lightning mast specimen;
[0055] FIGS. 33A-33B are graphical representations of coherence and
FRF from single and multiple impact tests of the mast shown in FIG.
32;
[0056] FIG. 34 illustrates a random series of pulses with the same
deterministic shape and random amplitudes and arrival times;
[0057] FIG. 35 illustrates an average normalized shape function of
force pulses;
[0058] FIGS. 36-38 are graphical representations of analytical and
numerical solutions;
[0059] FIG. 39 shows the comparison of the probability density
function (PDF) of the experimental arrival time with that from
uniform distribution;
[0060] FIG. 40 shows the comparison of the experimental and
analytical PDFs for the number of arrived pulses; and
[0061] FIG. 41 shows the comparison of the experimental and
analytical PDFs for the pulse amplitudes.
DETAILED DESCRIPTION OF PREFERRED EMBODIMENTS
[0062] 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.
[0063] 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.
[0064] 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.
[0065] 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.
[0066] 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.
[0067] 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.
[0068] 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.sub.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.sub.i (i=1, 2, . .
. , m), and its stiffness matrix, which depends linearly on
G.sub.i, is denoted by K=K(G), where G=[G.sub.1, G.sub.2, . . . ,
G.sub.m].sup.T. Here the superscript T denotes matrix transpose.
The eigenvalue problem of the structure with stiffness parameters
G.sub.i is
K.phi..sup.k=.lambda..sup.kM.phi..sup.k (1)
[0071] where M is the constant mass matrix, and
.lambda..sup.k=.lambda..su- p.k(G) and .phi..sup.k=.phi..sup.k(G)
(k=1, 2, . . . , N) are the k-th th eigenvalue and mass-normalized
eigenvector respectively. It is noted that
.lambda..sub.k=.omega..sub.k.sup.2, where .omega..sub.k is the k-th
natural frequency of the structure. The normalized eigenvectors of
(1) satisfy the orthonormality relations
(.phi..sup.k).sup.TM.phi..sup.u=.delta..sub.ku (.phi..sup.k).sup.T
K.phi..sup.u=.lambda..sup.k.delta..sub.ku (2)
[0072] where 1.ltoreq.u.ltoreq.N and .delta..sub.ku is the
Kronecker delta. Before the first iteration, the initial stiffness
parameters of the damaged structure are assumed to be
G.sub.i.sup.(0)=.sigma..sub.i.sup- .G.sub.hi (i=1, 2, . . . , m),
where 0<.sigma..sub.i.ltoreq.1, and the eigenvalue problem in
(1) corresponds to that of the structure with stiffness parameters
G.sub.i.sup.(0). 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 ].sub.i=1. Let
G.sub.di (i=1, 2, . . . , m) denote the stiffness parameters of the
damaged structure. The eigenvalue problem of the damaged structure
is
K.sub.d.phi..sub.d.sup.k=.lambda..sub.d.sup.kM.phi..sub.d.sup.k
(3)
[0073] where K.sub.d=K(G.sub.d) is the stiffness matrix with
G.sub.d=[G.sub.d1,G.sub.d2, . . . , G.sub.dm].sup.T, and
.lambda..sub.d.sup.k=.lambda..sup.k(G.sub.d) and
.phi..sub.d.sup.k=.phi..- sup.k(G.sub.d) are the k-th eigenvalue
and mass-normalized eigenvector respectively. The stiffness matrix
K.sub.d is related to K through the Taylor expansion 1 K d = K ( G
d ) = K + i = 1 m K G l G i ( 4 )
[0074] where .delta.G.sub.i=G.sub.di-G.sub.i (i=1, 2, . . . , m)
are the changes in the stiffness parameters, and the higher-order
derivatives of K with respect to G.sub.i vanish because K is
assumed to be a linear function of G.sub.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.
[0075] Let the k-th eigenvalue and mass-normalized eigenvector of
the damaged structure be related to .lambda..sup.k and .phi..sup.k
through 2 d k = k + i = 1 m ( 1 ) i k G i + i = 1 m j = 1 m ( 2 )
ij k G i G j + + i = 1 m j = 1 m t = 1 m p summations ( p ) ij t k
G i G j G t + e k , ( 5 ) d k = k + i = 1 m z ( 1 ) i k G i + i = 1
m j = 1 m z ( 2 ) ij k G i G j + + i = 1 m j = 1 m t = 1 m p
summations z ( p ) ij t k G i G j G t + e k ( 6 )
[0076] where .lambda..sub.(1)i.sup.k, .lambda..sub.(2)ij.sup.k, . .
. , and .lambda..sub.(p)ij . . . i.sup.k are the coefficients of
the first, second, . . . , and p-th order perturbations for the
eigenvalue, z.sub.(1)i.sup.k, z.sub.(2)ij.sup.k, . . . , and
z.sub.(p)ij . . . i.sup.k are the coefficient vectors of the first,
second, . . . , and p-th order perturbations for the eigenvector,
and e.sub..lambda..sup.k and e.sub..phi..sup.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.gtoreq.1, 3 p ! ( p ) ij t k = p k G i G j G t p ! z ( p ) ij t k
= p k G i G j G t ( 7 )
[0077] By (7), 4 ( p ) ij t k and z ( p ) ij t k
[0078] 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.
[0079] Using the normalization relations of the eigenvectors,
.phi..sup.k and .phi..sub.d.sup.k, and symmetry of the coefficient
vectors in (6), as indicated earlier, one obtains 5 1 = ( d k ) T M
d k = ( k + i = 1 m z ( 1 ) i k G i + i = 1 m j = 1 m z ( 2 ) ij k
G i G j + + i = 1 m j = 1 m t = 1 m p summations z ( p ) ij st k G
i G j G s G t + ) T M ( k + i = 1 m z ( 1 ) i k G i + i = 1 m j = 1
m z ( 2 ) ij k G i G j + + i = 1 m j = 1 m t = 1 m p summations z (
p ) ij st k G i G j G s G t + ) = 1 + i = 1 m [ ( k ) T Mz ( 1 ) i
k + ( z ( 1 ) i k ) T M k ] G i + i = 1 m j = 1 m 1 R ij 1 { 2 ! (
k ) T Mz ( 2 ) ij k + [ ( z ( 1 ) i k ) T Mz ( 1 ) j k + ( z ( 1 )
j k ) T Mz ( 1 ) i k ] + 2 ! ( z ( 1 ) ij k ) T M k } G i G j + i =
1 m j = 1 m t = j m 1 R ijt 1 { 3 ! ( k ) T Mz ( 3 ) ijt k + 2 ! [
( z ( 1 ) i k ) T Mz ( 2 ) j k + ( z ( 1 ) j k ) T Mz ( 2 ) it k +
( z ( 1 ) t k ) T Mz ( 2 ) ij k ] + 2 ! [ ( z ( 1 ) jl k ) T Mz ( 1
) i k + ( z ( 2 ) il k ) T Mz ( 1 ) j k + ( z ( 2 ) ij k ) T Mz ( 1
) l k ] + 3 ! ( z ( 3 ) ijl ) T M k } G i G j G l + + i = 1 m j = i
m t = s m p summations 1 R ij st 1 { p ! ( k ) T Mz ( p ) ij st k +
( p - 1 ) ! [ ( z ( 1 ) i k ) T Mz ( p - 1 ) jl st k + ( z ( 1 ) j
k ) T Mz ( p - 1 ) il st k + + ( z ( 1 ) t k ) T Mz ( p - 1 ) ij s
k ] + 2 ! ( p - 2 ) ! [ ( z ( 2 ) ij k ) T Mz ( p - 2 ) l st k + (
z ( 2 ) it k ) T Mz ( p - 2 ) j st + + ( z ( 2 ) st k ) T Mz ( p -
2 ) ij r k ] + + ( p - 2 ) ! 2 ! [ ( z ( p - 2 ) l st k ) T Mz ( 2
) ij k + ( z ( p - 2 ) j st k ) T Mz ( 2 ) il k + + ( z ( p - 2 )
ij r k ) T Mz ( 2 ) st k ] + ( p - 1 ) ! [ ( z ( p - 1 ) jl st k )
T Mz ( 1 ) i k + ( z ( p - 1 ) il st k ) T Mz ( 1 ) j k + + ( z ( p
- 1 ) ij s k ) T Mz ( 1 ) i k ] + p ! [ ( z ( p ) ij st k ) T M k }
G i G j G s G t + ( 8 )
[0080] where the superscript T denotes transpose. For any p-th
(p.gtoreq.1) order term in the last expression of (8), the
coefficient 6 R ij t 1
[0081] is defined by 7 R ij t 1 = X 1 ! X 2 ! X a ! ,
[0082] where X.sub.1, X.sub.2, . . . , and X.sub.u
(1.ltoreq..alpha..ltore- q.p) are the numbers of the first, second,
. . . , and last distinct indices within indices, i, j, . . . , and
t, with X.sub.1+X.sub.2+ . . . +X.sub.a=p. For instance,
R.sub.1234=1!1!1!1!=1 with a=4, and R.sub.112223=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.ltoreq.c.ltoreq.p+1) type
of terms is obtained by multiplying a (c-1)-th order term in the
expansion of (.phi..sub.d.sup.k).sup.T by a (p-c+1)-th order term
in the expansion of M.phi..sub.d.sup.k. For the c-th type of term,
where 2.ltoreq.c.ltoreq.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.ltoreq.c.ltoreq.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 8 R ij t 1 .
[0083] For instance, the 4.sup.th order term of the form
.delta.G.sub.1.delta.G.sub.2.sup.3 in the expansion of
(.phi..sub.d.sup.k).sup.TM.phi..sub.d.sup.k can be obtained from
the expression for the p-th order term in (8): 9 1 1 ! 3 ! { 4 ! (
k ) T Mz ( 4 ) 1222 k + 3 ! [ ( z ( 1 ) t k ) T Mz ( 3 ) 222 k + (
z ( 1 ) 2 k ) T Mz ( 3 ) 122 k + ( z ( 1 ) 2 k ) T Mz ( 3 ) 122 k +
( z ( 1 ) 2 k ) T Mz ( 3 ) 122 k ] + 2 ! 2 ! [ ( z ( 2 ) 12 k ) T
Mz ( 2 ) 22 k + ( z ( 2 ) 12 k ) T Mz ( 2 ) 22 k + ( z ( 2 ) 12 k )
T Mz ( 2 ) 22 k + ( z ( 2 ) 22 k ) T Mz ( 2 ) 12 k + ( z ( 2 ) 22 k
) T Mz ( 2 ) 12 k + ( z ( 2 ) 22 k ) T Mz ( 2 ) 12 k + 3 ! [ ( z (
3 ) 222 k ) T Mz ( 1 ) 1 k + ( z ( 3 ) 122 k ) T Mz ( 1 ) 2 k + ( z
( 3 ) 122 k ) T Mz ( 1 ) 2 k + ( z ( 3 ) 122 k ) T Mz ( 1 ) 2 k ] +
4 ! ( z ( 4 ) 1222 k ) T M k } G 1 G 2 3 = { 4 ( k ) T Mz ( 4 )
1222 k + [ ( z ( 1 ) 1 k ) T Mz ( 3 ) 222 k + 3 ( z ( 1 ) 2 k ) T
Mz ( 3 ) 122 k ] + 4 ( z ( 2 ) 12 k ) T Mz ( 2 ) 22 k + [ ( z ( 3 )
222 k ) T Mz ( 1 ) 1 k + 3 ( z ( 3 ) 122 k ) T Mz ( 1 ) 2 k ] + 4 (
z ( 4 ) 1222 k ) T M k } G 1 G 2 3 ( 9 )
[0084] where p=4 and the four indices involved, i, j, l, and o, are
i=1 and j=l=o=2.
[0085] Equating the coefficients of the .delta.G.sub.i (i=1, 2, . .
. , m) terms in (8) and using symmetry of the mass matrix
yields
(.phi..sup.k).sup.TMz.sub.(1)i.sup.k=0 (10)
[0086] Equating the coefficients of the
.delta.G.sub.1.delta.G.sub.j terms and using symmetry of M and
z.sub.(2)ij.sup.k yields 10 ( k ) T Mz ( 2 ) ij k = - 1 2 ! ( z ( 1
) i k ) T Mz ( 1 ) j k ( 11 )
[0087] for all i, j=1,2, . . . , m . Following a similar procedure,
one obtains 11 ( k ) T Mz ( 3 ) ijl k = - 2 ! 3 ! [ ( z ( 1 ) i k )
T Mz ( 2 ) jl + k + ( z ( 1 ) j k ) T Mz ( 2 ) li k + ( z ( 1 ) l k
) T Mz ( 2 ) ij k ] ( 12 )
[0088] for i,j,l=1, 2, . . . , m. Equating the coefficients of the
.delta.G.sub.i.delta.G.sub.j . . . .delta.G.sub.s.delta.G.sub.t
terms of p-th order in (8) yield 12 ( k ) T Mz ( p ) ij t k = - 1 2
( p ! ) { ( p - 1 ) ! [ ( z ( 1 ) i k ) T Mz ( p - 1 ) jl st k + (
z ( 1 ) j k ) T Mz ( p - 1 ) il st k + + ( z ( 1 ) t k ) T Mz ( p -
1 ) ij s ] + 2 ! ( p - 2 ) ! [ ( z ( 2 ) ij k ) T Mz ( p - 2 ) l st
+ ( z ( 2 ) il k ) T Mz ( p - 2 ) j st k + + ( z ( 2 ) st k ) T Mz
( p - 2 ) ij q k ] + + ( p - 2 ) ! 2 ! [ ( z ( p - 2 ) l st k ) T
Mz ( 2 ) ij k + ( z ( p - 2 ) j st k ) T Mz ( 2 ) il k + + ( z ( p
- 2 ) ij r k ) T Mz ( 2 ) st k ] + ( p - 1 ) ! [ ( z ( p - 1 ) jl
st k ) T Mz ( 1 ) i k + ( z ( p - 1 ) il st k ) T Mz ( 1 ) j k + +
( z ( p - 1 ) ij s k ) T Mz ( 1 ) t k ] } ( 13 )
[0089] 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 13 M and z ( p ) ij t k ,
the y - th ( 1 y p - 1 2 )
[0090] 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 14 y
- th ( 1 y p 2 - 1 )
[0091] 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
15 ( p 2 ) - th type ,
[0092] of terms in the middle of the expression.
[0093] Substituting (4)-(6) into (3) yields 16 { K + i = 1 m K G i
G i } { k + i = 1 m z ( 1 ) i k G i + i = 1 m j = 1 m z ( 2 ) ij k
G i G j + i = 1 m j = 1 m l = 1 m z ( 3 ) ij k G i G j G l + } = {
k + i = 1 m ( 1 ) i k G i + i = 1 m j = 1 m ( 2 ) ij k G i G j + i
= 1 m j = 1 m l = 1 m ( 3 ) ijl k G i G j G l + } M { k + i = 1 m z
( 1 ) i k G i + i = 1 m j = 1 m z ( 2 ) ij k G i G j + i = 1 m j =
1 m l = 1 m z ( 3 ) ijl k G i G j G l + } ( 14 )
[0094] Equating the coefficients of the .delta.G.sub.i (i=1,2, . .
. , m) terms in (14) yields 17 Kz ( 1 ) i k + K G i k = k Mz ( 1 )
i k + ( 1 ) i k M k ( 15 )
[0095] Expanding z.sub.(1)i.sup.k using normalized eigenvectors of
(1) as basis vectors, one has 18 z ( 1 ) i k = u = 1 N P ( 1 ) iu k
u ( 16 )
[0096] where 19 P ( 1 ) iu k
[0097] is the coefficient of the u-th term and the number in the
parentheses in its subscript corresponds to that of 20 z ( 1 ) i k
.
[0098] Pre-multiplying (15) by (.phi..sup.k).sup.T and using (1),
(2), and (16) yields 21 ( 1 ) i k = ( k ) T K G i k ( 17 )
[0099] Substituting (16) into (10) and using (2) yields 22 P ( 1 )
ik k = 0 ( 18 )
[0100] Pre-multiplying (15) by (.phi..sup..nu.).sup.T, where
1.ltoreq..nu..ltoreq.N and .nu..noteq.k, and using (1), (2), and
(16) yields 23 P ( 1 ) iv k = 1 k - v ( v ) T K G i k ( 19 )
[0101] By (16), (18) and (19) we have determined 24 z ( 1 ) i k
.
[0102] Equating the coefficients of the
.delta.G.sub.i.delta.G.sub.j terms in (14) yields 25 2 ! Kz ( 2 )
ij k + K G i z ( 1 ) j k + K G j z ( 1 ) i k = 2 ! k Mz ( 2 ) ij k
+ ( 1 ) i k Mz ( 1 ) j k + ( 1 ) j k MZ ( 1 ) i k + 2 ! ( 2 ) ij k
M k ( 20 )
[0103] Expanding 26 z ( 2 ) ij k
[0104] using normalized eigenvectors of (1) as basis vectors, one
has 27 z ( 2 ) ij k = u = 1 N P ( 2 ) iju k u ( 21 )
[0105] where 28 P ( 2 ) iju k
[0106] is the coefficient of the u-th term and the number in the
parentheses in its subscript corresponds to that of 29 z ( 2 ) ij k
.
[0107] Pre-multiplying (20) by (.phi..sup.k).sup.T and using (1),
(2), (10), and (21) yields 30 ( 2 ) ij k = 1 2 ! ( k ) T [ K G i z
( 1 ) j k + K G j z ( 1 ) i k ] ( 22 )
[0108] Substituting (21) into (11) and using (2) yields 31 P ( 2 )
ijk k = - 1 2 ! ( z ( 1 ) i k ) T Mz ( 1 ) j k ( 23 )
[0109] Pre-multiplying (20) by (.phi..sup..nu.).sup.T, where
1.ltoreq..nu..ltoreq.N and .nu..noteq.k, and using (1), (2), and
(21) yields 32 P ( 2 ) ijv k = 1 2 ! ( k - v ) ( v ) T { [ K G i z
( 1 ) j k + K G j z ( 1 ) i k ] - [ ( 1 ) i k Mz ( 1 ) j k + ( 1 )
j k Mz ( 1 ) i k ] } ( 24 )
[0110] By (21), (23) and (24) we have determined 33 z ( 2 ) ij k
.
[0111] Equating the coefficients of the
.delta.G.sub.i.delta.G.sub.j.delta- .G.sub.t terms in (14) yields
34 3 ! Kz ( 3 ) ijl k + 2 ! [ K G i z ( 2 ) jl k + K G j z ( 2 ) li
k + K G l z ( 2 ) ij k ] = 3 ! k Mz ( 3 ) ijl k + 2 ! [ ( 1 ) i k
Mz ( 2 ) jl k + ( 1 ) j k Mz ( 2 ) li k + ( 1 ) l k Mz ( 2 ) ij k ]
+ 2 ! [ ( 2 ) jl k Mz ( 1 ) i k + ( 2 ) li k Mz ( 1 ) j k + ( 2 )
ij k Mz ( 1 ) l k ] 3 ! ( 3 ) ijl k M k ( 25 )
[0112] Expanding 35 z ( 3 ) ijl k
[0113] using normalized eigenvectors of (1) as basis vectors, one
has 36 z ( 3 ) ijl k = u = 1 n P ( 3 ) ijlu k u ( 26 )
[0114] where 37 P ( 3 ) ijlu k
[0115] is the coefficient of the u-th term and the number in the
parentheses in its subscript corresponds to that of 38 z ( 3 ) ijl
k .
[0116] Pre-multiplying (25) by (.phi..sup.k).sup.T and using (1),
(2), (10), and (26) yields 39 ( 3 ) ijl k = 2 ! 3 ! ( k ) T [ K G i
z ( 2 ) jl k + K G j z ( 2 ) li k + K G l z ( 2 ) ij k - ( 1 ) i k
Mz ( 2 ) jl k - ( 1 ) j k Mz ( 2 ) li k - ( 1 ) l k Mz ( 2 ) ij k ]
( 27 )
[0117] Substituting (26) into (12) and using (2) yields 40 P ( 3 )
ijlk k = - 2 ! 3 ! [ ( z ( 1 ) i k ) T Mz ( 2 ) jl k + ( z ( 1 ) j
k ) T Mz ( 2 ) li k + ( z ( 1 ) l k ) T Mz ( 2 ) ij k ] ( 28 )
[0118] Pre-multiplying (25) by (.phi..sup..nu.).sup.T, where
1.ltoreq..nu..ltoreq.N and .nu..noteq.k, and using (1), (2), and
(26) yields 41 P ( 3 ) ijlv k = 2 ! 3 ! ( k - v ) ( v ) T K G i z (
2 ) jl k + K G j z ( 2 ) li k + K G l z ( 2 ) ij k - ( 1 ) i k Mz (
2 ) jl k - ( 1 ) j k Mz ( 2 ) li k - ( 1 ) l k Mz ( 2 ) ij k - ( 2
) ji k Mz ( 1 ) i k - ( 2 ) li k Mz ( 1 ) j k - ( 2 ) ij k Mz ( 1 )
l k ] ( 29 )
[0119] By (26), (28) and (29) we have determined 42 z ( 3 ) ijl k
.
[0120] We proceed now to derive the perturbation solutions for the
general p-th order terms in (5) and (6). Equating the coefficients
of the .delta.G.sub.i.delta.G.sub.j . . .
.delta.G.sub.s.delta.G.sub.t terms of order p in (14) yields 43 p !
Kz ( p ) ij st k + ( p - 1 ) ! [ K G i z ( p - 1 ) j st k + K G j z
( p - 1 ) i st k + + K G t z ( p - 1 ) ij s k ] = p ! k Mz ( p ) ij
st k + ( p - 1 ) ! [ ( 1 ) i k Mz ( p - 1 ) j st k + ( 1 ) j k Mz (
p - 1 ) ist k + + ( 1 ) t k Mz ( p - 1 ) ij s k ] + 2 ! ( p - 2 ) !
[ ( 2 ) ij k Mz ( p - 2 ) l st k + ( 2 ) il k Mz ( p - 2 ) j st k +
+ ( 2 ) st k Mz ( p - 2 ) ij q k ] + + p ! ( p ) ij st k M k ( 30
)
[0121] Expanding 44 z ( p ) ij st k
[0122] using normalized eigenvectors of (1) as basis vectors, one
has 45 z ( p ) ij st k = u = 1 n P ( p ) ij stu k u ( 31 )
[0123] where 46 P ( p ) ij stu k
[0124] is the coefficient of the u-th term and the number in the
parenthesis in its subscript corresponds to that of 47 z ( p ) ij
st k .
[0125] Pre-multiplying (30) by (.phi..sup.k).sup.T and using (1),
(10), (31) and orthonormality relations of eigenvectors yields 48 (
p ) ij st k = 1 p ! ( k ) T [ ( p - 1 ) ! ( K G i z ( p - 1 ) j st
k + K G j z ( p - 1 ) i st k + + K G t z ( p - 1 ) ij s k ) - ( p -
1 ) ! ( ( 1 ) i k Mz ( p - 1 ) j st k + ( 1 ) j k Mz ( p - 1 ) i st
k + + ( 1 ) t k Mz ( p - 1 ) ij s k ) - 2 ! ( p - 2 ) ! ( ( 2 ) ij
k Mz ( p - 2 ) l st k + ( 2 ) il k Mz ( p - 2 ) j st k + + ( 2 ) st
k Mz ( p - 1 ) ij r k ) - - ( p - 2 ) ! 2 ! ( ( p - 2 ) l st k Mz (
2 ) ij k + ( p - 2 ) j st k Mz ( 2 ) il k + + ( p - 2 ) ij r k Mz (
2 ) st k ) ] ( 32 )
[0126] 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 49 P ( p ) ij stk
k = 1 2 ( p ! ) { ( p - 1 ) ! [ ( z ( 1 ) i k ) T Mz ( p - 1 ) jlst
k + ( z ( 1 ) j k ) T Mz ( p - 1 ) ilst k + + ( z ( 1 ) t k ) T Mz
( p - 1 ) ijs k ] + 2 ! ( p - 2 ) ! [ ( z ( 2 ) ij k ) T Mz ( p - 2
) lst k + ( z ( 2 ) il k ) T Mz ( p - 2 ) jst k + + ( z ( 2 ) st k
) T Mz ( p - 2 ) ijq k ] + ( p - 2 ) ! 2 ! [ ( z ( p - 2 ) lst k )
T Mz ( 2 ) ij k + ( z ( p - 2 ) jst k ) T Mz ( 2 ) il k + + ( z ( p
- 2 ) ijr k ) T Mz ( 2 ) st k ] + ( p - 1 ) ! [ ( z ( p - 1 ) jlst
k ) T Mz ( 1 ) i k + ( z ( p - 1 ) ils k ) T Mz ( 1 ) j k + + ( z (
p - 1 ) ij s k ) T Mz ( 1 ) t k ] } ( 33 )
[0127] Pre-multiplying (30) by (.phi..sup..nu.).sup.T, where
1.ltoreq..nu..ltoreq.N and .nu..noteq.k, and using (1), (2), and
(31) yields 50 P ( P ) ijstv k = 1 p ! ( k - v ) ( v ) T [ ( p - 1
) ! ( K G i z ( p - 1 ) j ... st k + K G j z ( p - 1 ) i . st k + +
K G t z ( p - 1 ) ijs k ) - ( p - 1 ) ! ( ( 1 ) i k Mz ( p - 1 )
jst k + ( 1 ) j k Mz ( p - 1 ) ist k + + ( 1 ) t k Mz ( p - 1 ) ijs
k ) - 2 ! ( p - 2 ) ! ( ( 2 ) ij k Mz ( p - 2 ) lst k + ( 2 ) il k
Mz ( p - 2 ) jst k + + ( 2 ) st k Mz ( p - 2 ) ijr k ) - - ( p - 1
) ! ( ( p - 1 ) jst k Mz ( 1 ) i k + ( p - 1 ) ist k Mz ( 1 ) j k +
+ ( p - 1 ) ijs k Mz ( 1 ) t k ) ] ( 34 )
[0128] By (31), (33) and (34) we have determined 51 z ( p ) ijst k
.
[0129] The p-th order sensitivities of eigenvectors can be
subsequently obtained from (7). They depend on the eigenvalue and
eigenvector sensitivities of orders up to p-1.
[0130] Equations (5) and (6) serve both the forward and inverse
problems. In the former one determines the changes in the
eigenparameters with changes in the stiffness parameters. Damage
detection is treated as an inverse problem, in which one identifies
iteratively the changes in the stiffness parameters from a selected
set of measured eigenparameters of the damaged structure:
.lambda..sub.d.sup.k (k=1, 2, . . . , n.sub..lambda.) and
.phi..sub.d.sup.k (k=1, 2, . . . , n.sub..phi.), where
1.ltoreq.n.sub..lambda., n.sub..phi..ltoreq.N. Here
.lambda..sub.d.sup.k and .phi..sub.d.sup.k are assumed to be the
perfect eigenparameters; simulated noise is included in the
measured eigenparameters in the beam and frame examples that
follow. Often we choose a set of n measured eigenparameter pairs to
detect damage, i.e., n.sub..lambda.=n.sub..phi.=n. Let the number
of the measured degrees of freedom of .phi..sub.d.sup.k be N.sub.m;
N.sub.m=N and N.sub.m<N when we have complete and incomplete
eigenvector measurements, respectively. With reduced measurements
the unmeasured degrees of freedom of .phi..sub.d.sup.k is estimated
first using a modified eigenvector expansion method (see the beam
and frame examples below) and .phi..sub.d.sup.k is mass-normalized
subsequently. Only the component equations corresponding to the
measured degrees of freedom of .phi..sub.d.sup.k are used in (6).
The system equations in Eqs. (5) and (6) involves
n.sub..lambda.+n.sub..phi.N.sub.m scalar equations with m unknowns,
which are in general determinate if n.sub..lambda.+n.sub..phi.N-
.sub.m=m, under-determined if
n.sub..lambda.+n.sub..phi.N.sub.m<m, and over-determined if
n.sub..lambda.+n.sub..phi.N.sub.m>m. In the first iteration,
.lambda..sup.k and .phi..sup.k in (5) and (6) correspond to the
eigenparameters of the structure with the initial stiffness
parameters G.sub.i.sup.(0). With the perturbation terms determined
as shown above, the changes in the stiffness parameters
.delta.G.sub.i.sup.(1), where the number in the superscript denotes
the iteration number, can be solved from (5) and (6) using an
optimization method discussed below. The estimated stiffness
parameters of the damaged structure are updated by
G.sub.i.sup.(1)=G.sub.i.sup.(0)+.delta.G.sub.i.s- up.(1). With the
updated stiffness parameters, the eigenparameters, .lambda..sup.k
and .phi..sup.k, in (5) and (6) are re-calculated from the
eigenvalue problem (1) and the perturbation terms are re-evaluated.
One subsequently finds .delta.G.sub.i.sup.(2) using the same method
as that in the first iteration. This process is continued until the
termination criterion,
.vertline..delta.G.sub.i.sup.(L).vertline.<.epsilon., where L is
the last iteration number and .epsilon. is some small constant, is
satisfied for all i=1,2, . . . , m. Note that after the w-th
(1.ltoreq.w<L) iteration, we set G.sub.i.sup.(w) to G.sub.hi if
G.sub.i.sup.(w)>G.sub.hi, and to zero or some small stiffness
value .epsilon..sub.G if G.sub.i.sup.(w)<0. A flowchart for the
iterative algorithm is shown in FIG. 1B. When a single iteration is
used, the iterative method becomes a direct method.
[0131] FIG. 1B shows steps included in a method for detecting
structural damage in accordance with one embodiment of the present
invention. The method includes as an initial step measuring one or
more eigenparameters, .lambda..sub.d.sup.k, .phi..sub.d.sup.k
(Block 1). These eigenparameters are then compared with estimated
eigenparameters associated with the stiffness parameters,
G.sub.i.sup.(0) (Block 2). The differences between the measured and
estimated eigenparameters are then used by the perturbation method
to establish system equations (5) and (6) (Block 3). The
perturbation method may be a first or higher order
multiple-parameter perturbation procedure.
[0132] Next, an optimization method may be used to find the changes
in the stiffness parameters .delta.G.sub.i.sup.(w) (Block 4). These
values are then compared to the predetermined value .epsilon.
(Block 5), and if the absolute values are less than .epsilon. the
stiffness parameters identified are set as G.sub.i.sup.(w-1) (Block
6).
[0133] If the absolute values are not all less than .epsilon., the
process proceeds along an iterative path where the stiffness
parameters are first updated (Block 7). The stiffness parameters
are then bounded between 0 or .epsilon..sub.G and G.sub.hi (Block
8), and eigenparameters associated with the updated stiffness
parameters are calculated (Block 9). The comparison of Block 2 is
then performed based on these calculated, or estimated,
eigenparameters.
[0134] Optimization Methods
[0135] Neglecting the residuals of order p+1 in (5) and (6) yields
a system of polynomial equations of order p. When
n.sub..lambda.+n.sub..phi- .N.sub.m.ltoreq.m.
.delta.G.sub.i.sup.(w) (i=1,2, . . . , m) at the w-th iteration can
be solved from the resulting equations. There are generally an
infinite number of solutions when
n.sub..lambda.+n.sub..phi.N.sub.m<- ;m, a unique solution
when n.sub..lambda.+n.sub..phi.N.sub.m=m and p=1, and a finite
number of solutions when n.sub..lambda.+n.sub..phi.N.sub.m=m and
p>1. When n.sub..lambda.+n.sub..phi.N.sub.m>m, one generally
cannot find .delta.G.sub.i.sup.(w) to satisfy all the equations,
and an optimization method can be used to find
.delta.G.sub.i.sup.(w) which minimize an objective function related
to the errors in satisfying these equations. We use here the same
notations, e.sub..lambda..sup.k and e.sub..phi..sup.k, to denote
the errors in satisfying the system equations in (5) and (6)
respectively. Consider the objective function 52 J = k = 1 n W k (
e k ) 2 + k = 1 n W k ( e k ) T ( e k ) ( 35 )
[0136] where W.sub..lambda..sup.k (k=1,2, . . . , n,) and
W.sub..phi..sup.k (k=1, 2, . . . , n.sub..phi.) are the weighting
factors, and J is a function of .delta.G.sub.i.sup.(w) when one
substitutes the expressions for e.sub..lambda..sup.k and
e.sub..phi..sup.k in (5) and (6) into (35). When the first-order
perturbations are retained in (5) and (6), the least-squares method
can be used to determine .delta.G.sub.i.sup.(w) which minimize (35)
with W.sub..lambda..sup.k=W.sub..phi..sup.k=1. Other weighting
factors can be handled by dividing (5) and (6) by 53 1 W k and 1 W
k
[0137] respectively. The method determines essentially the
generalized inverse of the resulting system matrix, and is also
known as the generalized inverse method. When the perturbations up
to the p-th (p.gtoreq.1) order are included in (5) and (6), the
gradient and quasi-Newton methods [25] can be used to determine
.delta.G.sub.i.sup.(w) iteratively. Unlike the generalized inverse
method, the methods are applicable when other objective functions
are defined. While the optimization methods are introduced for
over-determined systems, they can be used when
n.sub..lambda.+n.sub..phi.N.sub.m.ltoreq.m, in which case J=0
(i.e., e.sub..lambda.=e.sub..phi.=0) when the optimal solutions are
reached.
[0138] Gradient Method
[0139] To minimize the objective function in (35) at the w-th
iteration, one can use the algorithm 54 G ( b ) ( w ) = G ( b - 1 )
( w ) - b f b - 1 ( 36 )
[0140] to update the changes in the stiffness parameters, where 55
G ( b ) ( w ) = ( G 1 ( b ) ( w ) , G 2 ( b ) ( w ) , , G m ( b ) (
w ) ) T ,
[0141] .alpha..sub.b>0 is the step size, and f.sub.b-1 equals
the gradient vector 56 g b - 1 = ( J G 1 ( w ) , J G 2 ( w ) , , J
G m ( w ) ) T
[0142] associated with 57 G ( b - 1 ) ( w ) .
[0143] Note that the subscript b (b.gtoreq.1) in all the variables
in (36) denotes the number of nested iterations. The initial values
used are 58 G i ( 0 ) ( w ) = 0.
[0144] The nested iteration is terminated when
.alpha..sub.b.parallel.g.su-
b.b-1.parallel..sub..infin.<.gamma., where
.parallel..multidot..paralle- l..sub..infin. is the infinity norm
and .gamma. is some small constant, or the number of nested
iterations exceeds an acceptable number, D.
[0145] Quasi-Newton Methods
[0146] Due to its successive linear approximations to the objective
function, the gradient method can progress slowly when approaching
a stationary point. The quasi-Newton methods can provide a remedy
to the problem by using essentially quadratic approximations to the
objective function near the stationary point. The iteration scheme
of these methods is given by (36) with
f.sub.b-1=B.sub.b-1g.sub.b-1, where B.sub.b-1 is an approximation
to the inverse of the Hessian matrix used at the b-th nested
iteration, and the other variables the same as those previously
discussed. Initially, we set 59 G i ( 0 ) ( w ) = 0
[0147] and B.sub.0=I, the identity matrix. The matrix B.sub.b is
updated using either the DFP formula 60 B b = B b - 1 + ( G ( b ) (
w ) - G ( b - 1 ) ( w ) ) ( G ( b ) ( w ) - G ( b - 1 ) ( w ) ) T (
G ( b ) ( w ) - G ( b - 1 ) ( w ) ) T ( g b - g b - 1 ) - [ B b - 1
( g b - g b - 1 ) ] [ B b - 1 ( g b - g b - 1 ) ] T ( g b - g b - 1
) T B b - 1 ( g b - g b - 1 ) ( 37 )
[0148] or the BFGS formula 61 B b = B b - 1 + [ 1 + ( g b - g b - 1
) T B b - 1 ( g b - g b - 1 ) ( G ( b ) ( w ) - G ( b - 1 ) ( w ) )
T ( g b - g b - 1 ) ] ( G ( b ) ( w ) - G ( b - 1 ) ( w ) ) ( G ( b
) ( w ) - G ( b - 1 ) ( w ) ) T ( G ( b ) ( w ) - G ( b - 1 ) ( w )
) T ( g b - g b - 1 ) - ( G ( b ) ( w ) - G ( b - 1 ) ( w ) ) ( g b
- g b - 1 ) T B b - 1 + B b - 1 ( g b - g b - 1 ) ( g b - g b - 1 )
T ( ( G ( b ) ( w ) - G ( b - 1 ) ( w ) ) T ( g b - g b - 1 ) ) (
38 )
[0149] The nested iteration is terminated when
.alpha..sub.b.parallel.B.su-
b.b-1g.sub.b-1.parallel..sub..infin.<.gamma. or the number of
iterations exceeds D. A flowchart for the quasi-Newton methods,
including the step size search procedure as outlined below, is
shown in FIG. 2.
[0150] Step Size Search Procedure
[0151] The optimal step size for the gradient and quasi-Newton
methods is determined in each nested iteration to minimize the
function 62 J ( G ( b - 1 ) ( w ) - b f b - 1 ) = F ( b )
[0152] with respect to .alpha..sub.b. The search procedure is
divided into two phases: an initial search to bracket the optimum
.alpha.*.sub.b and a golden section search to locate .alpha.*.sub.b
within the bracket, as shown in FIG. 2.
[0153] Initial Bracketing. Choose the starting point x.sub.1=0 and
an initial increment .DELTA.>0. Let x.sub.2=x.sub.1+.DELTA.,
F.sub.1=F(x.sub.1), and F.sub.2=F(x.sub.2). Since for the gradient
and quasi-Newton methods, f.sub.0=g.sub.0 and it is along a descent
direction of J when .DELTA. is sufficiently small, one has
F.sub.2<F.sub.1. Rename 2.DELTA. as .DELTA., and let
x.sub.3=x.sub.2+.DELTA. and F.sub.3=F(x.sub.3). If
F.sub.3>F.sub.2, stop and .alpha.*.sub.b is contained in the
interval (x.sub.1, x.sub.3). Otherwise, rename x.sub.2 as x.sub.1
and x.sub.3 as x.sub.2, then F.sub.2 becomes F.sub.1 and F.sub.3
becomes F.sub.2. Rename 2.DELTA. as .DELTA., and let
x.sub.3=x.sub.2+.DELTA. and F.sub.3=F(x.sub.3). Compare F.sub.3 and
F.sub.2, and repeat the above procedure if F.sub.3<F.sub.2 until
F.sub.3>F.sub.2, with the final interval (x.sub.1, x.sub.3)
containing .alpha.*.sub.b.
[0154] Golden Section Search. If
.vertline.x.sub.3-x.sub.2.vertline.>.v-
ertline.x.sub.2-x.sub.1.vertline., define a new point:
x.sub.4=x.sub.2+0.382(x.sub.3-x.sub.2) (39)
[0155] Otherwise, rename x.sub.1 as x.sub.3 and x.sub.3 as x.sub.1,
and then define x.sub.4 using (40). Let F.sub.4=F(x.sub.4). If
F.sub.2<F.sub.4, rename x.sub.4 as x.sub.3, then F.sub.4 becomes
F.sub.3. Otherwise, rename x.sub.2 as x.sub.1 and x.sub.4 as
x.sub.2, then F.sub.2 becomes F.sub.1 and F.sub.4 becomes F.sub.2.
Compare .vertline.x.sub.3-x.sub.2.vertline. and
.vertline.x.sub.2-x.sub.1.vertlin- e., and repeat the above
procedure until .vertline.x.sub.3-x.sub.1.vertlin-
e..parallel.f.sub.b-1.parallel..sub..infin.<.epsilon..sub..alpha.,
where .epsilon..sub..alpha. is some small constant. Then choose 63
b * = x 1 + x 2 2 .
MASS-SPRING SYSTEM EXAMPLE
[0156] The algorithm discussed above is used to identify the
stiffness parameters of a N-degree-of-freedom system consisting of
a serial chain of masses and springs, such as the system shown in
FIG. 3. Let the masses of the system be M.sub.i=1 kg (i=1, 2, . . .
, N), and the stiffnesses of the undamaged springs be G.sub.hi=1N/m
(i=1, 2, . . . , m), where m=N+1. The system is said to have a
small, medium and large level of damage if the maximum reduction in
the stiffnesses is within 30%, between 30 and 70%, and over 70%,
respectively. The mass matrix M is an N.times.N identity matrix,
and the stiffness matrix K is a banded matrix with entries
K.sub.ii=G.sub.i+G.sub.i+1 (i=1,2, . . . , N),
K.sub.i(i+1)=K.sub.(i+1)i=-G.sub.i+1 (i=1, 2, . . . , N-1), and all
other entries equal to zero. The matrices 64 K G 1 and K G N
[0157] have a unit value in entries (1, 1) and (N, N),
respectively, and zero entries elsewhere. The nonzero entries of
the matrices 65 K G i ( i = 2 , 3 , , N - 1 ) are ( K G i ) ( i - 1
) ( i - 1 ) = ( K G i ) ii = 1 and ( K G i ) i ( i - 1 ) = ( K G i
) ( i - 1 ) i = - 1.
[0158] We look at a forward problem first with N=3 and m=4. The
stiffnesses of the damaged system are G.sub.d1=G.sub.d3=1N/m,
G.sub.d2=0.3N/m and G.sub.d4=0N/m. The undamaged system is
considered as the unperturbed system and the damaged system as the
perturbed system. Based on the eigenparameters of the undamaged
system, the eigensolutions of the damaged system are obtained using
the first-, second-, and third-order perturbations, as shown in
Table 1 below. The results show that even with the large changes in
stiffness, the third-order perturbation solutions compare favorably
with the exact solutions for the damaged system. The higher-order
perturbation solutions can be used for large order systems when
their direct eigensolutions become costly.
1TABLE 1 Bigensolutions of the damaged system from an eigenvalue
problem solver (exact) and perturbation analysis Perturbed
Eigenparameters Exact Unperturbed 1st order 2nd order 3rd order
.lambda..sup.1 0.10602 0.58579 0.30576 0.15670 0.10090
.lambda..sup.2 1.27538 2.00000 1.15000 1.25500 1.29344
.lambda..sup.3 2.21859 3.41421 2.14424 2.18830 2.20567 .phi..sup.1
66 { 0.16516 0.65733 0.73528 } 67 { 0.50000 0.70711 0.50000 } 68 {
0.28523 0.68836 0.74129 } 69 { 0.16111 0.66444 0.79452 } 70 {
0.12907 0.65040 0.76653 } .phi..sup.2 71 { 0.95541 0.07840 -
0.28470 } 72 { 0.70711 - 0.00000 - 0.70711 } 73 { 0.95459 0.10607 -
0.45962 } 74 { 0.00188 0.14319 - 0.31776 } 75 { 0.97976 0.12310 -
0.26810 } .phi..sup.3 76 { 0.24478 - 0.74952 0.61507 } 77 { 0.50000
- 0.70711 0.50000 } 78 { 0.36477 - 0.72586 0.60871 } 79 { 0.29635 -
0.74132 0.62482 } 80 { 0.26445 - 0.74875 0.62362 }
[0159] Consider now the damage detection problem with N=9, m=10,
G.sub.d5=0.5N/m, G.sub.d8=0.7N/m, G.sub.d10=0.8N/m and
G.sub.di=1N/m (i=1,2,3,4,6,7,9). We set
W.sub..lambda..sup.k=W.sub..phi..sup.k=1, .epsilon.=0.001,
.gamma.=10.sup.-10, .sigma..sub.i=1 for all i, and D=500; the
actual numbers of nested iterations in all the cases are much
smaller than D. Since vanishing stiffness in any spring other than
the two end ones in FIG. 1 can result in two decoupled subsystems,
when G.sub.i.sup.(w)<0 we set G.sub.i.sup.(w) to
.epsilon..sub.G=0.1N/m in the first two iterations and to 0.1N/m in
the remaining iterations for all the cases considered here. A
relatively large value is assigned to .epsilon..sub.G in the
initial iterations to avoid close eigenvalues in the mass-spring
system and small denominators in (19). This improves convergence
especially when a small number of eigenparameter pairs are used. A
smaller value is used for .epsilon..sub.G in the later iterations
to improve the accuracy of stiffness estimation when there is a
large level of stiffness reduction. Using the first-order
perturbations and different numbers of eigenparameter pairs, the
maximum errors in estimating the stiffnesses of the damaged system
at the w-th iteration, defined by 81 E = max 1 i N G i ( w ) - G di
G di ( 40 )
[0160] are shown in FIGS. 4A and 4B for all the iterations. In FIG.
4A p=1 and n=1,2,3, and in FIG. 4B p=1 and n=4,5, . . . ,9. When
n=1, the error decreases slowly, though monotonically, and there is
an estimation error of 1.5% at the end of iteration. While the
errors can increase with the iteration number before approaching
zero for n=3, they decrease monotonically for n=2 and n.gtoreq.4.
All the stiffnesses are exactly identified at the end of iteration
when n.gtoreq.2. Note that the number of the system equations
equals and exceeds the number of unknowns when n=1 and n.gtoreq.2,
respectively. Since the system equations are linear, they have a
unique solution when n=1, and J has a unique minimum when
n.gtoreq.2. With the small .gamma. the gradient method and the
quasi-Newton methods using the DFP and BFGS formulas yield exactly
the same results as the generalized inverse method (not shown
here). Because the generalized inverse method does not involve any
nested iteration, it is the most efficient one among the four
methods. While not shown here, the results indicated that that the
quasi-Newton methods converge faster than the gradient method and
the BFGS method has the similar performance to the DFP method. In
what follows the BFGS method will be used with the higher-order
perturbations. With the second-order perturbations the errors shown
in FIG. 4C decrease monotonically for all n (n=1,2, . . . , 9). The
errors at w=1 in the expanded view decrease in the order
n=3,2,4,9,7,8, with the lines for n=5 and 7 virtually
indistinguishable. In FIG. 4(A-C) the following symbols are used
for different n: 407 n=1; n=2; n=3; n=4; n=5; n=6; n=7 n=8; n=9.
While the use of the second-order perturbations improves the
accuracy of stiffness estimation in each iteration and reduces the
number of iterations, it takes a much longer time to compute the
higher-order perturbations and the associated optimal
solutions.
[0161] When only the first few eigenvalues are used, for instance,
n.sub..lambda.=5 and n.sub..phi.=0, the stiffnesses identified with
the first-order perturbations,
G.sup.(4)=(0.875, 0.976, 0.926, 0.864, 0.699, 0.699, 0.864, 0.926,
0.976, 0.875).sup.T N/m (41)
[0162] where the number in the superscript denotes the last
iteration number, correspond to those of a different system with
the same eigenvalues for the first five modes as the damaged
system. The same stiffnesses are identified with the second-order
perturbations. Similarly, when the first eigenvector is used, i.e.,
n.sub..phi.=1 and n.sub..lambda.=0, the stiffnesses identified with
the first-order perturbations are those of a different system with
the same eigenvector for the first mode as the damaged system:
G.sup.(9)=(0.989, 0.991, 0.995, 1, 0.520, 0.829, 0.940, 0.668,
0.959, 0.768).sup.T N/m (42)
[0163] With the second-order perturbations the stiffnesses of the
damaged system are identified. The stiffnesses identified are not
unique because the system equations in each iteration are
under-determined. The solution given by the generalized inverse
method here in each iteration is the minimum norm solution.
Increasing the number of eigenparameters used can avoid this
problem.
[0164] If the system has a large level of damage, i.e.,
G.sub.d5=0.3 N/m, G.sub.d10=0.1 N/m, with the other parameters
unchanged, the stiffnesses of the damaged system are identified
with the first-order perturbations after 55 iterations when n=1 and
6 iterations when n=2 and 3 as shown in FIG. 5A. For n=4,5,6,7, the
errors decrease monotonically and the number of iterations is
reduced slightly, as shown in FIG. 5B, which has an expanded view
near w=1. With the second-order perturbations the errors shown in
FIG. 5C decrease monotonically for n=1,2, . . . , 5, and the number
of iteration for n=1 is reduced from 55 in FIG. 5A to 4. The errors
at w=1 in the expanded view in FIG. 5C decrease in the order
n=3,2,4,5, with the lines for n=2 and 4 virtually
indistinguishable. The symbols used in FIG. 5A-C for n=1, 2, . . .
, 7 are the same as those in FIG. 4A-C.
[0165] Finally, consider a large order system with a large level of
damage: N=39, m=40, G.sub.d12=0.7N/m, G.sub.d19=G.sub.d37=0.1 N/m,
G.sub.d28=0.8 N/m, G.sub.d1=1 N/m (1.ltoreq.i.ltoreq.40 and
i.noteq.12,19,28,37), and the other parameters are the same as
those in the previous example. With the first-order perturbations
the exact stiffnesses are identified after 57 iterations when n=1,
as shown FIG. 6A. Using a larger number of eigenparameter pairs
(n=2,3,4,5) significantly reduces the number of iterations, as
shown in FIG. 6B. The errors at w=2 in the expanded view in FIG. 6B
decrease in the order n=4,2,3,5. With the second-order
perturbations the exact stiffnesses are identified after 6
iterations when n=1 and 3 iterations when n=2, as shown in FIG. 6C,
which has an expanded view for 1.ltoreq.w.ltoreq.4. The symbols
used in FIG. 6(A-C) for n=1, 2, . . . , 5 are the same as those in
FIGS. 4A-C and 5A-C.
FIXED-FIXED BEAM EXAMPLE
[0166] The algorithm discussed above may be applied to detecting
structural damage in an aluminum beam with fixed boundaries. The
beam of length L.sub.t=0.7 m, width W=0.0254 m, and thickness
H=0.0031 m has an area moment of inertia 82 I = 1 12 WH 3 = 6.3058
.times. 10 - 11 m 4
[0167] and a mass density .rho.=2715 kg/m.sup.3. The finite element
method is used to model its transverse vibration. The beam is
divided into N.sub.e elements, as shown in FIG. 7, with the length
of each element being 83 l e = L i N e .
[0168] There are N.sub.e+1 nodes. With V.sub.i and .theta..sub.i
denoting the translational and rotational displacements at node i
(i=1, 2, . . . , N.sub.e+1), the displacement vector of the i-th
(i=1,2, . . . , N.sub.e) element is
[V.sub.i,.theta..sub.i,V.sub.i+1,.theta..sub.i+1].sup.T. The
Young's modulus is assumed to be constant over each beam element
and that of the i-th element is denoted by G.sub.i. The Young's
modulus of the undamaged beam is G.sub.h=69.times.10.sup.9
N/m.sup.2. Hence G.sub.hi=G.sub.h for i=1,2, . . . , m, where
m=N.sub.e. Small to large levels of damage correspond to reductions
of the moduli defined for the Mass-Spring Example. The mass and
stiffness matrices of the i-th beam element are 84 M i e = WHl e
420 [ 156 - 22 l e 54 13 l e - 22 l e 4 l e 2 - 13 l e - 3 l e 2 54
- 13 l e 156 22 l e 13 l e - 3 l e 2 22 l e 4 l e 2 ] , K i e = G i
I l e 3 [ 12 6 l e - 12 6 l e 6 l e 4 l e 2 - 6 l e 2 l e 2 - 12 -
6 l e 12 - 6 l e 6 l e 2 l e 2 - 6 l e 4 l e 2 ] ( 43 )
[0169] Using the standard assembly process yields the
2(N.sub.e+1).times.2(N.sub.e+1) global mass and stiffness matrices.
Constraining the translational and rotational displacements of the
two nodes at the boundaries to zero yields the N.times.N M and K
matrices, where N=2(N.sub.e-1) is the degrees of freedom of the
system. The displacement vector of the system, involving the
displacements of the 2.sup.nd through N.sub.e-th node, is
[V.sub.2,.theta..sub.2,V.sub.3,.thet- a..sub.3, . . .
,V.sub.N.sub..sub.e, .theta..sub.N.sub..sub.e].sup.T. The matrix 85
K G i
[0170] (i=1, 2, . . . , m) can be obtained from K by setting
G.sub.i=1 and G.sub.1= . . . =G.sub.i-1=G.sub.i+1= . . .
=G.sub.N=0. The parameters W.sub..lambda..sup.k, W.sub..phi..sup.k,
.epsilon., .gamma., .sigma..sub.i, and D are set to the same values
as those previously discussed, and .epsilon..sub.G is set to 0.15
G.sub.h in the first two iterations and to 0.05 G.sub.h in the
remaining iterations. The first-order perturbations are used unless
indicated otherwise.
[0171] Consider first the cases with N.sub.e=m=10 and N=18. When
the system has a medium level of damage:
G.sub.d=(1, 1, 1, 1, 0.5, 1,1, 0.7, 1, 0.8).sup.T.times.G.sub.h
(44)
[0172] the stiffness parameters of the damaged system are
identified after 6 iterations with n=1. When the system has a large
level of damage:
G.sub.d=(1, 1, 1, 1, 0.3, 1, 1, 0.7, 1, 0.1).sup.T.times.G.sub.h
(45)
[0173] the stiffness parameters of the damaged system are
identified after 7 iterations with n=1. Consider next the cases
with N.sub.e=20, 40 and 80. For the systems with medium and large
levels of damage, the stiffness parameters of the first 10 elements
are given by (44) and (45), respectively, and those of the
remaining elements are G.sub.h. In all the cases the stiffness
parameters of the damaged systems are identified within 10
iterations when n=1. The numbers of iterations are reduced slightly
when the second-order perturbations are used. Note that the system
equations are over-determined when n=1.
[0174] When only the translational degrees of freedom of an
eigenvector are measured, a modified eigenvector expansion method
is used to estimate the unmeasured rotational degrees of freedom.
To this end, .phi..sub.d.sup.k is partitioned in the form
.phi..sub.d.sup.k=[(.phi..su- b.dm.sup.k).sup.T,
(.phi..sub.du.sup.k).sup.T].sup.T, where .phi..sub.dm.sup.k and
.phi..sub.du.sup.k are the measured and unmeasured degrees of
freedom of .phi..sub.d.sup.k, respectively. Similarly, .phi..sup.k
in (6) is partitioned in the form .phi..sup.k=[(.phi..sub.m.s-
up.k).sup.T, (.phi..sub.u.sup.k).sup.T].sup.T, where
.phi..sub.m.sup.k and .phi..sub.u.sup.k correspond to the measured
and unmeasured components of .phi..sub.d.sup.k, respectively. Since
.phi..sub.dm.sup.k, .phi..sub.m.sup.k and .phi..sub.u.sup.k are
known in each iteration, .phi..sub.du.sup.k is estimated from
.phi..sub.du.sup.k=k[(.phi..sub.m.su-
p.k).sup.+.phi..sub.dm.sup.k].phi..sub.u.sup.k, where the
superscript+denotes generalized inverse. Once the rotational
degrees of freedom of .phi..sub.d.sup.k are determined,
.phi..sub.d.sup.k and .phi..sup.k are converted to their original
forms and .phi..sub.d.sup.k is mass-normalized. Only the component
equations corresponding to the measured degrees of freedom of
.phi..sub.d.sup.k are used in (6) and the system equations are
determinate when n=1. The exact stiffness parameters of the damaged
systems considered above can be identified. For the 10- and
20-element beams with the medium levels of damage, the stiffness
parameters are identified after 6 and 24 iterations, respectively,
when n.sub..lambda.=1 and n.sub..phi.=2. For the 10- and 20-element
beams with the large levels of damage, the stiffness parameters are
identified after 9 iterations when n.sub..lambda.=1 and
n.sub..phi.=3 and 10 iterations when n=2, respectively.
[0175] Finally, the effects of measurement noise on the performance
of the algorithm are evaluated for the 10-element beam with the
large level of damage. Simulated noise is included in the measured
eigenparameters:
.lambda..sub.d.sup.k=.lambda.*.sub.d.sup.k+.nu.R.sub..lambda..sup.k.lambda-
.*.sub.d.sup.k,
.phi..sub.d.sup.k=.phi.*.sub.d.sup.k+.nu.R.sub..phi..sup.k-
.phi.*.sub.d.sup.k (46)
[0176] where .lambda.*.sub.d.sup.k and .phi.*.sub.d.sup.k are the
k-th perfect eigenvalue and eigenvector, respectively,
R.sub..lambda..sup.k is a uniformly distributed random variable in
the interval [-1,1], R.sub..phi..sup.k is a diagonal matrix whose
diagonal entries are independently, uniformly distributed random
variables in the interval [-1,1], and .nu..epsilon.[0,1] is the
noise level. Note that R.sub..lambda..sup.k and R.sub..phi..sup.k
are generated for each measured mode. Each random parameter is
generated 10 times and the average is used. Three different noise
levels are considered: .nu.=5%, 10% and 20%. When all the degrees
of freedom of an eigenvector are measured, the stiffness parameters
identified with n=1, 2, and 3 are shown in FIGS. 8A, 8B, and 8C,
respectively. The following symbols are used for different noise
levels: , .nu.=5%; , .nu.=10%; , .nu.=20%. When only the
translational degrees of freedom of an eigenvector are measured,
the eigenvector expansion method described above is used and the
stiffness parameters identified with n=2 and n=3 are shown in FIGS.
9A and 9B, respectively. The same symbols as those in FIG. 8A-C are
used in FIG. 9A-B for different noise levels. The stiffness
parameters corresponding to .nu.=0 in FIGS. 8 and 9 are the exact
values. It is seen that in the presence of noise, the stiffness
parameters can be accurately identified with an increased number of
measured eigenparameters.
[0177] Space Frame Example
[0178] The damage detection method can be applied to more complex
structures, such as, for example, the modular, four bay space frame
shown in FIG. 10. In this example and all the following examples
the first order perturbation approach along with the generalized
inverse method is used. The four-bay space frame example shown in
FIG. 10 is made of extruded aluminum with 48 members, and is 8'
tall, 1.46' wide and 1.8' deep. The horizontal and diagonal members
have the same cross-sectional dimensions (25.4 mm.times.12.7 mm)
and the vertical members are "L"-angles with cross-sectional
dimensions 50.8 mm.times.50.8 mm.times.6.35 mm (thickness). All the
members are connected through bolted joints. The frame is assumed
to be fixed to the ground. The finite element method is used to
model the 3-dimensional vibration of the frame, with each member
modeled with four 12-degrees-of-freedom beam elements. The total
degrees of freedom are 960 when the boundary conditions are
applied. The Young's modulus is assumed to be constant over each
member and that of the i-th (1.ltoreq.i.ltoreq.m=48) member is
denoted by G.sub.i. The Young's modulus of the undamaged member is
G.sub.h=69.times.10.sup.9 N/m.sup.2 and G.sub.hi=G.sub.h. A
vertical (member 1) and a diagonal (member 48) member in the first
bay are assumed to have 70% and 30% reductions in Young's modulus,
respectively, and all the other members are assumed to be
undamaged. The first two vibration modes, corresponding to the
bending of the frame along the x and y directions, respectively,
are used to detect damage. The translational degrees of freedom of
the 16 nodes, denoted by "S" in FIG. 10, along the x and y
directions are assumed to be measured with 10% measurement noise.
All the other degrees of freedom are estimated using the
eigenvector expansion method discussed above and the measured modes
are mass-normalized. The Young's moduli of all the members of the
damaged truss are identified as shown in FIG. 10, with the maximum
estimation error less than 7%. Note that the truss has closely
spaced vibration modes. With the modes arranged in the order of
increasing frequencies in each iteration, the method has
effectively handled mode switching that has occurred in the damage
detection process.
[0179] In summary, the damage detection method identifies stiffness
parameters in structures, which have a small, medium, and large
level of damage if the maximum reduction in the stiffnesses is
within 30%, between 30 and 70%, and over 70%, respectively: A large
level of damage is studied in many examples because this poses the
most challenging case, with sever mismatch between the
eigenparameters of the damaged and undamaged structures.
[0180] The damage detection method as embodied and broadly
described herein can be applied to structures that can be modeled
with beam elements. A beam element is an element that has one
dimension that is much longer than the other two. This element is
very good at modeling "I"-beams, rectangular beams, circular beams,
"L"-angles, "C"-channels, pipes, and beams with varying cross
sections. Structures that can be modeled with this element include,
but are not limited to, lightning masts, light poles, traffic
control poles, pillar type supports, bridges, pipelines, steel
building frameworks, television, radio, and cellular towers, space
structures, cranes, pipelines, railway tracks, and vehicle frames.
Structures that can be modeled as beam elements are used simply for
ease of discussion, and it is well understood that the damage
detection method discussed above may be applied to structures that
can be modeled by other elements using the finite element method or
modeled by using other methods.
[0181] Damage Detection Using Changes of Natural Frequencies:
Simulation and Experimental Validation
[0182] For structures such as beams and lightning masts in electric
substations, using only the changes in the natural frequencies can
relatively accurately detect the location(s) and extent of damage,
even though the system equations are severely underdetermined in
each iteration. This is an interesting finding as it is much easier
to measure the natural frequencies than the mode shapes, and
demonstrates the effectiveness of the iterative algorithm.
Extensive numerical simulations on beams and lightning masts
confirmed this finding. Experiments on the beam test specimens with
different damage scenarios and a lightning mast in an electric
substation validated the simulation results. The beam test
specimens and the lightning mast are used as examples for
demonstration purposes, and the method can be applied to other
structures. Note that unlike the beam example shown earlier, where
the Euler-Bernoulli beam finite element model is used, the
Timoshenko beam finite element model is used in all the examples
here. The Timoshenko beam theory is found to be more accurate in
predicting the natural frequencies of the lightning masts and
circular beams than the Euler-Bernoulli beam theory.
[0183] For a cantilever beam, simulation results show that the
damage located at a position within 0-35% and 50-95% of the length
of the beam from the cantilevered end can be easily detected with
less than 5 measured natural frequencies, and the damage located at
a position within 35-50% of the length of the beam from the
cantilevered end and at a position within 5% of the length of the
beam from the free end can be relatively accurately detected with
10-15 measured natural frequencies.
[0184] Numerical and Experimental Verification
[0185] Cantilever Aluminum Beams
[0186] Experimental damage detection results for four different
scenarios are shown first, followed by various simulation
results.
[0187] Scenario 1: Evenly-Distributed Damage Machined from the Top
and the Bottom Surfaces of the Beam Test Specimen.
[0188] The aluminum beam test specimen shown in FIG. 12 is 45 cm
long by 2.54 cm wide by 0.635 cm thick. It is divided into 40
elements (each element has a length of 1.125 cm). The beam has a
section (from approximately 10 cm to 15 cm from the cantilevered
end) of 5 cm long and 7.62E-4 m thick machined both from the top
and the bottom surfaces of the beam. This corresponds to 56% of
damage (or reduction of bending stiffness EI) along the length of
five elements (from the 9.sup.th to the 13.sup.rd element). Using
the changes of the first 2 to 5 measured natural frequencies,
damage is detected within 7 elements using 2 or 5 measured
frequencies (from the 7.sup.th to the 13.sup.rd element with 5
measured frequencies and from the 5.sup.th to the 11.sup.th element
with 2 measured frequencies) and within 8 elements using 3 or 4
measured natural frequencies (both from the 6.sup.th to the
13.sup.rd element), as shown in FIG. 13. With 5 measured
frequencies the average damage of the 7 elements in FIG. 13, from
6.75 cm to 14.625 cm, is 46%. Note that the elements with damage
less than 10% are not accounted for here. The extent of damage
detected is slightly lower than the actual extent because the
predicted damage occurs at 2 mote elements (the 7.sup.th and
8.sup.th elements) than the actual one. The error results from the
solution of the severely underdetermined system equations (5
equations with 80 unknowns).
[0189] Scenario 2: The Same Aluminum Beam Test Specimen as in
Scenario 1 Clamped at the Other End.
[0190] The same beam was tested in the clamped-free configuration
with the clamped end reversed, which placed the damage from 25 cm
to 30 cm (from the 23.sup.rd to the 27.sup.th element) from the
cantilevered end. The damage detection results with 3 to 5 measured
frequencies are shown in FIG. 14. With 5 measured frequencies the
average damage between 23.625 cm to 32.625 cm (from the 22.sup.nd
to the 29.sup.th element) is 40%. Again the elements with damage
less than 10% are not included.
[0191] Scenario 3: Undamaged Cantilever Aluminum Beam Test Specimen
with the Same Dimensions as Above.
[0192] An undamaged aluminum beam test specimen was clamped at one
end with the same configuration as shown in FIG. 12. With the first
3 to 4 measured frequencies the maximum error for the estimated
bending stiffnesses of all the elements is within 10%, as shown in
FIG. 15.
[0193] Scenario 4: A Cut of Small Width on a Cantilever Aluminum
Beam Test Specimen, Shown in FIG. 16 with the Same Dimensions as
Above.
[0194] The beam shown below has a cut that is 0.4191 cm deep and
0.1016 cm wide, which corresponds to a 96% reduction in bending
stiffness at the cut. The beam is divided into 45 elements and the
cut is located in the middle of the 23.sup.rd element. With 3 to 5
measured natural frequencies, the damage is detected at several
elements surrounding the 23.sup.rd element, as shown in FIG. 17. It
was found there was a roughly 50%, 40%, and 60% reduction in
stiffness at the 22.sup.nd, 23.sup.rd, and 24.sup.th element,
respectively. The estimated damage extent at the above elements is
lower than that at the cut because the damage is distributed along
these elements.
[0195] To examine the effectiveness and robustness of the damage
detection algorithm, various simulations with different damage
scenarios were carried out. In this way, we can gain more insight
concerning the accuracy of the finite element model, convergence of
the estimated bending stiffnesses of all the elements of the beam
with the increased numbers of measured natural frequencies and/or
mode shapes, and region of the beam within which the damage can be
detected with few measured frequencies. With the cantilever
aluminum beam divided into 40 elements, the following two
simulations have the similar damage location and extent to those in
Scenarios 1 and 2 in the experiments:
[0196] Simulation 1: Uniform Damage Between 9.0 cm and 15.75 cm
from the Cantilevered End of the Beam with the Same Dimensions as
Discussed Above.
[0197] Simulation 2: Uniform Damage Between 29.25 cm and 36 cm from
the Cantilevered End of the Beam with the Same Dimensions as
Discussed Above.
[0198] Simulation results in FIG. 18 and FIG. 19 show that the
damage can be relatively accurately detected with the first 3 to 5
measured frequencies, and can be accurately detected with the first
10 measured frequencies.
[0199] With the cantilever aluminum beam divided into the same
number of elements, two more simulations are presented here: one
for a multiple damage scenario--70% of damage at the 3.sup.rd
element and 30% of damage at the 20.sup.th element, and the other
for a 50% of uniform damage from the 16.sup.th to the 18.sup.th
element (i.e., the damage is located at a position within 35-50% of
the length of the beam from the cantilevered end).
[0200] Simulation 3: Multiple Damage of the Beam with the Same
Dimensions as Discussed Above: 70% of Damage at the 3.sup.rd
Element and 30% of Damage at the 20.sup.th Element.
[0201] Simulation results in FIG. 20 show that the multiple damage
can be relatively accurately detected with the first 3 to 5
measured frequencies, and accurately detected with the first 10
measured frequencies.
[0202] Simulation 4: Uniform Damage from the 16.sup.th to the
18.sup.th Element of the Beam with the Same Dimensions as Discussed
Above.
[0203] Simulation results in FIG. 21 show that the damage within
35-50% of the length of the beam from the cantilever end can be
accurately detected with 10-15 measured frequencies.
[0204] Lightning Masts
[0205] The lighting mast shown in FIG. 22 has two sections of
constant cross-sections and an eccentric spike. The lengths of the
lower and upper sections are 6.89 m and 6.83 m, respectively, and
that of the spike above the upper section is 1.4375 m. Both mode
shapes and natural frequencies were measured for this mast. The
mode shapes were measured by using a laser Doppler vibrometer. The
natural frequencies were measured using both the laser Doppler
vibrometer and an accelerometer. It takes about 0.5-1.5 days to
measure the mode shapes and about 30 minutes to measure the natural
frequencies. It would be even more difficult to measure the mode
shapes for some of the taller lightning masts, which are 100 and
130 feet tall. A finite element model was made using OpenFEM. Once
the model was completed the measured and calculated natural
frequencies and mode shapes were compared. The mast was expected to
be undamaged, since by inspection the measured and calculated
natural frequencies matched (Table 2). It was found that the first
mode shape measurement is affected by wind and high modes are less
affected by the wind. The measured frequencies are not affected by
the wind and can be used for damage detection.
2TABLE 2 Comparison of measured and calculated (from the finite
element model) natural frequencies Mode # Measured FEM Error % 1
1.17 1.18 -0.8648 2 5.10 5.23 -2.4104 3 7.83 7.89 -0.8418 4 8.17
8.67 -6.1104 5 16.60 16.55 0.31327 6 29.26 30.36 -3.7385 7 45.32
47.85 -5.5831 8 47.75 50.87 -6.5191 9 54.10 55.81 -3.1439 10 67.45
68.83 -2.0476 11 74.70 78.99 -5.7456
[0206] Damage detection was then performed using only the first 4
to 7 measured natural frequencies, as shown in Table 2. The mast is
modeled by 40 elements; the lower section has 18 elements, the
upper section has 15 elements, and the spike has 7 elements. The
experimental damage detection results are shown in FIG. 23. The
mast is modeled by 40 elements with 20 elements in each of the two
sections. It appears that the elements around the joints near the
middle of the mast (8.2 m in FIG. 23) and at the free end of the
mast (13.72 m in FIG. 23) had some damage. This is most likely due
to the fact that the joints were modeled in the finite element
model as being infinitely stiff, which is almost never the case.
Also, the spike is actually connected to the free end of the mast
at two points. In the finite element model here it is assumed that
the spike is attached to the mast along the whole line of contact
that has a length of 0.34 m, which makes the model a little
stiffer. Similar damage detection results were obtained with
different numbers of measured frequencies.
[0207] Simulations for the same lightning mast as shown in FIG. 22
with different damage scenarios are carried out. Shown in FIG. 24
are the simulation results for the mast with 70% damage between
0.76 m and 1.15 m from the ground. Shown in FIG. 25 are the
simulation results for the mast with 50% damage between 6.89 m and
7.35 m from the ground. In both cases the location and extent of
damage can be relatively accurately detected using the first 5
measured natural frequencies and accurately detected using the
first 8 or 10 measured frequencies.
[0208] Methods to Handle Some Ill-Conditioned System Equations
[0209] When the first order perturbation approch is used, one needs
to solve in each iteration a system of linear algebraic
equations
A.delta.G=F (47)
[0210] which are the linearized equations of (5) and (6), where A
is the system matrix, F is the vector representing the differences
between the measured and estimated natural frequencies and/or mode
shapes, .delta.G is the optimal changes in the stiffness parameters
to be found. While the gradient and quasi-Newton methods can be
used, the generalized inverse method is most efficient because it
does not involve nested iterations, and thus .delta.G=A.sup.+F,
where A.sup.+ is the generalized inverse of A. The problem (47) may
be ill posed for certain cases, especially in the first several
iterations, because A.sup.+ can be relatively large and small
changes in F can result in large changes in the solution .delta.G.
Since the stiffness parameters cannot be negative or greater than
the corresponding values of the undamaged structure, the solution
may not converge.
[0211] Ill-conditioning problems do not occur in all the examples
described above. Sometimes they can occur. Consider, for example, a
cantilever aluminum beam of the same dimensions as those of the
beam shown in FIG. 12. The beam is divided into 36 elements and
assumed to have 30% of damage at the 5.sup.th element and 90% of
damage at the 18.sup.th element. The translational degrees of
freedom of the first or second mode shape vector are used to detect
damage. The system equations are determinate, and the
ill-conditioning problem occurs in the iterations. When a natural
frequency is also used in additional to an incomplete eigenvector
described above, the ill-conditioning problem can also occur. The
following regularization methods can be used to handle the
ill-conditioning problem:
[0212] Method 1. Estimate A.sup.+ from (A.sup.TA+.eta.I).sup.-1
A.sup.T, where .eta. is a small positive constant that can be
searched in several ways. One way is set
.eta.=.eta..times.1.618.times.min(10,.parallel..delt-
a.G.parallel..sub..infin.), where .mu..multidot..mu..sub..infin.
denotes the infinity norm, with the initial value
.eta.=.parallel.A.sup.TA.parall- el..sub..infin., until
.parallel..delta.G.parallel..sub..infin. is less than 0.8, for
example.
[0213] Method 2. To constrain the magnitude of .delta.G, we include
it in the objective function to be minimized. For example, we can
minimize the following objective function 86 f ( G ) = i = 1 N o j
= 1 m ( F i - A ij G j ) 2 + j = 1 m ( G i ) 2 ( 48 )
[0214] instead of (35), where
N.sub.o=n.sub..lambda.+n.sub..phi.N.sub.m is the number of
equations in (5) and (6). The optimal solution can be obtained by
using the generalized inverse method for the expanded system
A*=[A;I] and the expanded vector F*=[F;0], where I is the m.times.m
identity matrix and 0 is the m.times.1 zero vector.
[0215] Note that the solutions from the regularization methods are
not strictly the optimal solutions for the original objective
function in (35). With accurate and sufficient measurement
information and proper handling of the ill-conditioning problem,
the system equations can become well conditioned in the last few
iterations and regulation does not need to be applied. Consequently
the stiffness parameters can be more accurately determined.
Sometimes regulation may over constrain the magnitude of .delta.G,
and the termination criterion
.parallel..delta.G.parallel..sub..infin.<.epsilon. may be
satisfied during the regulation process. In this case, we may set
87 G i = 3 G i ; G r; .infin.
[0216] to amplify the magnitude of .delta.G, and the iteration is
terminated after the system equations become well conditioned.
Since the second regulation method does not need to search for
.eta., it can be more effective when it works. When one carries out
the damage detection procedures using different combinations of
measured frequencies and mode shapes and obtains the same or
similar stiffness parameters, one can have sufficient confidence on
the results obtained.
[0217] Using only the translational degrees of freedom of the first
eigenvector and the first regulation method, the estimated bending
stiffnesses of all the elements of the beam are shown in FIG. 26A.
Using the translational degrees of freedom of the second
eigenvector and the first regulation method, the estimated bending
stiffnesses of all the elements of the beam are shown in FIG. 26B.
In both cases the locations of damage are exactly detected and the
extent is relatively accurately determined. Using the translational
degrees of freedom of the first eigenvector along with the first
natural frequency and the first regulation method, the bending
stiffnesses of all the elements of the beam are exactly determined,
as shown in FIG. 26C.
[0218] Similarly, using the translational degrees of freedom of the
first eigenvector and the second regulation method, the estimated
bending stiffnesses of all the elements of the beam are shown in
FIG. 27A. Using the translational degrees of freedom of the second
eigenvector and the second regulation method, the estimated bending
stiffnesses of all the elements of the beam are shown in FIG. 27B.
In both cases the locations of damage are exactly detected and the
extent is relatively accurately determined. Using the translational
degrees of freedom of the first eigenvector along with the first
natural frequency and the second regulation method, the bending
stiffnesses of all the elements of the beam are exactly determined,
as shown in FIG. 27C.
[0219] Conclusions
[0220] Thus, the sensitivities of eigenparameters of all orders
may, for the first time, be derived using a multiple-parameter,
general-order perturbation method. The higher-order solutions may
be used to estimate the changes in the eigenparameters with large
changes in the stiffness parameters. The perturbation method may be
combined with an optimization method to form a robust iterative
damage detection algorithm. The gradient and quasi-Newton methods
can be used for the first or higher order system equations, and the
generalized inverse method can be used efficiently with the first
order system equations because it does not involve nested
iterations. Including the higher-order perturbations can
significantly reduce the number of iterations when there is a large
level of damage. A modified eigenvector expansion method is used to
estimate the unmeasured component of the measured mode shape. For
many cases, the location(s) and extent of damage can be relatively
accurately detected using only measured natural frequencies.
Methods to handle ill-conditioned system equations that may
occasionally arise are developed and shown to be effective.
Numerical simulations on different structures including spring-mass
systems, beams, lightning masts, and frames show that with a small
number of measured eigenparameters, the stiffness parameters of the
damaged system may be accurately identified in all the cases
considered. Experiments on the different beam test specimens and
the lightning mast in an electric substation validated the
theoretical predictions. The methodology can be readily applied to
various operation structures of different sizes by incorporating
their finite element models or other mathematical models.
[0221] One example of a practical application of this damage
detection method is shown in FIG. 28, in which a structure 10 may
represent any type of structure that can be modeled by beam
elements, on which periodic damage assessment must be conducted.
This includes, but is not limited to, structures such as lightning
masts, utility poles, cell towers, and other such structures. A
structure that can be modeled by beam elements is used simply for
ease of discussion, and it is well understood that the general
order perturbation method discussed above may be applied to a
number of different structural members without departing from the
spirit of the invention as embodied and broadly described
herein.
[0222] The structure 10 divided into elements E.sub.1 through
E.sub.n, and with nodes N.sub.1 through N.sub.n, is equipped with a
sensor 40, such as, for example, an accelerometer. The sensor 40 is
configured to measure a response to an impact F or a force from a
shaker applied to the structure 10. The impact F may be applied in
the form of a single impact, as shown in FIG. 28, or it may be in
the form of a series of random impacts F.sub.1 through F.sub.n, as
shown in FIG. 29. The location on the beam where the impact F is
applied may be varied, but the impact F is preferably not applied
at one of the nodal points of vibration modes whose natural
frequencies and mode shapes need to be measured. In either
instance, the sensor 40 senses a response in the structure 10 to
the impact F or the force from the shaker, and transmits that
response to a processor 20 equipped with, among other elements, the
damage detection method 30 discussed above.
[0223] In accordance with the method as described above, upon
receipt of the response signal from the sensor 40, the processor 20
accesses and performs the method 30 as previously discussed, and,
as the data converges, determines an extent and location, by
element, of any structural damage present in the structure 10.
[0224] Although the convergence of the system to resolution is not
dependent on where the impact F is applied to the structure 10, a
signal to noise ratio of the response may be increased, and an
efficiency of the system optimized, as the application point of the
impact F is moved away from a fixed point N.sub.0, if one needs to
measure the natural frequencies and/or mode shapes of lower modes.
Similarly, although the system will converge regardless of the
average energy and timing of the impact F or series of impacts
F.sub.1-F.sub.n applied to the structure 10, a signal to noise
ratio of the response may be increased, and an efficiency of the
system optimized based on a random series of impacts or the random
shaker excitation to the structure 10, if the structure is large or
slightly nonlinear or if there is an ambient excitation to the
structure such as wind.
[0225] The structure 10 may be excited in a number of different
manners. For example, an impact F may be applied manually or a
shaker may be used, at any given point on the structure 10,
excluding the fixed point N.sub.0 or the nodal points of vibration
modes whose natural frequencies and/or mode shapes need to be
measured. However, in an effort to increase the signal to noise
ratio in the signal provided by the sensor 40 to the processor 30
if the structure is large or slightly nonlinear or if there is an
ambient excitation to the structure such as wind and to provide the
convenience and portability over the shaker test, a series of
random impacts F.sub.1-F.sub.n shown in FIG. 29 may be applied
manually or by a specially designed device to generate the impact
F. Although accurate and reliable damage assessments can be
achieved regardless of how the impact F is applied to the structure
10, results may be improved using a random series impact method,
which involves using a series of impacts F of random amplitudes and
random arrival times. A series of random impacts has been shown to
increase an energy input to the structure 10, improve the signal to
noise ratio, especially in such situations as strong wind
excitation, and average out slight nonlinearities that arise, for
example, from bolted joints and extract linearized
eigenparameters.
[0226] FIG. 29 illustrates a structure 10, sensor 40 and processor
20 similar to those shown in FIG. 28. An impact F in FIG. 29 is
applied to the structure 10, at a location other than N.sub.0 or
one of the nodal points of vibration modes whose natural
frequencies and/or mode shapes need to be measured, as a series of
random impacts F.sub.1 through F.sub.n delivered at random
amplitude and arrival time. These random impacts F.sub.1-F.sub.n
may be delivered manually, or, as shown in FIG. 30, may be
delivered in an automated fashion through the use of a random
impact hammer device 1550.
[0227] Different methods have been employed in conventional
vibration testing in order to excite a test specimen. Shaker
testing, in which a specimen is, simplistically, shaken in order to
impart a high level of energy, can produce a high signal to noise
ratio, and can induce random excitation, which can average out
slight nonlinearities and extract linearized eigenparameter
parameters. However, shaker testing is not practically employed in
the field on relatively large structures, and can be cost
prohibitive to conduct. Single impact hammer testing addresses the
shortfalls of shaker testing, in that it is portable and
inexpensive to conduct. However, single impact hammer testing falls
short where shaker testing is strong, in that low energy input of
single impact hammer testing produces a low energy input, a low
signal to noise ratio with no randomization. To address the need
for a system which combines the advantages of shaker testing and
single impact hammer testing, a Random Impact Series method for
hammer testing is presented which yields a high energy, high signal
to noise ratio, random system. A novel stochastic model is
developed to simulate the random impact series produced manually
and to generate a random impact series for a specially designed
random impact device.
[0228] FIG. 31 shows a random impact device 1550, according to one
embodiment of the invention. Random impact device 1550 includes
random signal generating unit 1560 and a random impact actuator
1570 with impact applicator 1580. The impact applicator 1580
preferably includes a sensor 1581, such as a force transducer,
attached at its tip. Sensor 1581 is preferably configured to send
data to the spectrum analyzer in stiffness parameter unit 103 in
order to obtain mode shape information, as discussed above in
connection with FIG. 1A. Sensor 1581 can be coupled to the spectrum
analyzer in any manner know in the art including, but not limited
to, wire, optical fiber or even a wireless connection. Random
impact signal generating unit 1560 generates outputs 1590 and 1591
to random impact actuator 1570. It should be appreciated that
outputs 1590 and 1591 could be coupled to the random impact
actuator 1570 in any manner known in the art including, but not
limited, to cables, wires, optical fibers, wireless communications,
and so forth.
[0229] Output 1590 corresponds to the amplitude .psi..sub.i, as
discussed below. In particular, random signal generating unit 1560
outputs a value corresponding to .psi..sub.i, as discussed below.
Output 1591 corresponds to .tau..sub.i, as discussed below.
[0230] Impact applicator 1580 is shown, for purposes of
illustration, as impacting structure 1600 in a reciprocating
motion. Impact applicator 1580 is shown to have an impact path
1610, such that when a structure 1600 lies within the impact region
1610, impact applicator 1580 impacts structure 1600 with a force of
random amplitude that arrives at a random time, as discussed below.
The impact applicator 1580 and impact region 1610 as shown in FIG.
31 is one possible embodiment, and other shapes and impact regions
may be used while still falling within the scope of the present
invention.
[0231] Although random signal generating unit 1560 is shown to
output two outputs 1590 and 1591 that correspond to .psi..sub.i and
.tau..sub.i below, it should be understood that signal generating
unit 1590 could output any signal or signals that ultimately
results in random impact actuator 1570 driving impact applicator
1580 to impact structure 1600 with random arrival times .tau..sub.i
and random amplitudes .psi..sub.i.
[0232] As discussed above, vibration information is obtained by
sensor 110, and is sent to stiffness parameter unit 103 via sensor
coupler 113. The stiffness parameter unit 103 sends stiffness
parameters to the damage information processor 123. While the
random impact device can be used for damage detection purposes as
discussed above, it can certainly be used just for obtaining the
natural frequencies and/or mode shapes of the structure for modal
testing purposes.
[0233] Experiments conducted on lightning masts by the applicant
confirm that multiple impact testing performs better than single
impact testing when there is wind excitation to the masts. Results
are shown for a 65 foot tall mast with a 5 foot spike, as shown in
FIG. 32, referred to hereinafter simply as a seventy foot tall
mast. The mast has two constant cross section schedule 40 pipes of
equal length. It can be seen from the results shown in FIGS.
33A-33B that the multiple impact test has much better coherence, is
less noisy away from the resonances, and can pick up the modes that
were missing in the single impact test. The multiple impact test is
better at exciting the lower modes that can also be excited by the
wind, which can be seen by comparing the frequency response
functions (FRF) between 0 and 30 Hz. It can been seen from the FRF
for the multiple impact tests that there are some modes that are
missed at near 4 Hz and 7 Hz with the single impact test, and the
mode at 10 Hz is much improved compared to the single impact
test.
[0234] A stochastic model will now be discussed which describes the
random impact series F.sub.1-F.sub.n, modeled as a Poisson process.
The Poisson process is one of a general class of processes that
arise in problems concerning the counting of events in the course
of time. The force pulses in the series are assumed to have an
arbitrary, deterministic shape function and random amplitudes and
arrival times. The force signal in a finite time interval is shown
to consist of wide sense stationary and non-stationary parts. The
expectations of the average power densities associated with the
entire force signal and the stationary part of the signal are
derived and compared, and the power spectral density is related to
the average power density and the autocorrelation function
associated with the stationary part of the force signal. Numerical
simulation is conducted to validate analytical predictions, and a
relationship between the Fourier transform and the discrete Fourier
transform used in a numerical simulation is produced. Experiments
on the four bay space frame, as shown in FIG. 10, validated the
distributions of the random variables and the Poisson process.
[0235] Stochastic Model of a Random Impact Series
[0236] A random impact series is modeled here as a sum of force
pulses with the same shape and random amplitudes and arrival times:
88 x ( t ) = o = 1 N ( t ) i y ( t - i ) t ( 0 , .infin. ) ( 49
)
[0237] where t is time, x(t) is the time function of the force
signal, .tau..sub.i.epsilon.(0, t] is the random arrival time of
the i-th pulse, and y(t-.tau..sub.i) is the deterministic shape
function for all the pulses, .psi..sub.i is the random variable
describing the amplitude of the i-th pulse, and N(t) is the number
of the pulses that have arrived during the time interval (0,t] and
is modeled as a Poisson process with stationary increments. All the
pulses are assumed to be of width .DELTA..tau., and
y(t-.tau..sub.i) satisfies y(t-.tau..sub.i)=0 if t<.tau..sub.i
and t>.tau..sub.i+.DELTA..tau..
[0238] Since a finite time record is used for the force signal in
modal testing, we consider only the pulses that arrive during the
time interval (0,T] of length T, as shown in FIG. 34. Note that
y(0) and y(.DELTA..tau.) do not have to be equal to zero in this
model. Because the pulses arriving after time t (t<T) will not
affect the force signal at time t, the random process N(t) in (49)
can be replaced by N(T). Also, if a pulse arrives at time T, it
will vanish at time T+.DELTA..tau.. The last pulse in FIG. 34
arrives at time T.sub.N(T) and ends at a time between T and
T+.DELTA..tau.. To include completely this last possible pulse, we
consider the time interval t .epsilon.(0, T+.DELTA..tau.]. Equation
(49) is rewritten as 89 x ( t ) = i = 1 N ( T ) i y ( t - i ) , i (
0 , T ] and t ( 0 , T + ] ( 50 )
[0239] For the Poisson process N(t) with stationary increments, the
probability of the event {N(t)=n}, where n is an integer, is 90 P {
N } ( n , t ) = - t ( t ) n n ! ,
[0240] where .lambda. is the constant arrival rate of the pulses.
By replacing t with T, the probability of the event {N(T)=n} is 91
P { N } ( n , T ) = - t ( t ) n n ! ( 51 )
[0241] All the arrival times .tau..sub.i, where i=1,2,3, . . . ,
N(T), are identically distributed, mutually independent random
variables. Because the arrival rate of the pulses is constant, the
arrival times .tau..sub.i are uniformly distributed in [0, T] with
the probability density function 92 P i ( ) = { 1 T , 0 T 0 ,
elsewhere ( 52 )
[0242] Similarly, .psi..sub.i (i=1,2, . . . , N(T)) are identically
distributed random variables, which are mutually independent and
independent of the distribution of the arrival times .tau..sub.i.
While the distribution of .psi..sub.i (i=1,2, . . . , N(T)) is not
used in the subsequent derivation, it is assumed that .psi..sub.i
satisfy the Gaussian distribution with the probability density
function 93 p i ( ) = 1 2 - ( - ) 2 2 2 , 0 < < .infin. ( 53
)
[0243] where .mu.=E[.psi..sub.i] is the mean and
.sigma..sup.2=E[.psi..sub- .1.sup.2]-E.sup.2[.psi..sub.i] is the
variance of .psi..sub.i. While the distribution in (53) is not used
in the analytical derivations, it is used in the numerical
simulation and validated experimentally for a manually applied
random impacts on a four bay space frame, as shown in FIG. 10.
[0244] Force Spectrum and its Expectation
[0245] The force spectrum is the Fourier transform of the force
signal in (50): 94 X ( j ) = F [ x ( t ) ] = - .infin. .infin. { i
= 1 N ( T ) i y ( t - i ) } - j t t = i = 1 N ( T ) i i i + y ( t -
i ) - j t t ( 54 )
[0246] where F denotes Fourier transform, X(j.omega.) is the
Fourier transform of x(t), .omega. is the angular frequency, and
j={square root}{square root over (-1)}. Let u=t-.tau..sub.i, then
t=u+.tau..sub.i and dt=du. Equation (54) becomes 95 X ( j ) = i = 1
N ( T ) i - i 0 y ( u ) - j u u ( 55 )
[0247] Since the integral in (55) is a deterministic function, the
expectation of the force spectrum is 96 E [ X ( j ) ] = E [ i = 1 N
( T ) i - j i ] 0 y ( u ) - j u u ( 56 )
[0248] Using the expression for conditional expectations,
E[E[V.vertline.U]]=E(V), where U=N(T) and 97 V = i = 1 N ( T ) i -
j i ,
[0249] we have from (56) 98 E [ X ( j ) ] = E [ E { [ i = 1 N ( T )
i - j i ] N ( T ) } ] 0 y ( u ) - j u u = n = 0 .infin. P { N } ( n
, T ) E [ i = 1 N ( T ) i - j i ] 0 y ( u ) - j u u ( 57 )
[0250] Since .psi..sub.ie.sup.-j.omega..tau..sup..sub.i (i=1,2, . .
. , n) are independent of each other and .psi..sub.i is independent
of e.sup.-j.omega..tau..sup..sub.i, we have 99 E [ i = 1 n i - j i
] = i = 1 n E [ i - j i ] = i = 1 n E [ i ] E [ - j i ] ( 58 )
[0251] Since .psi..sub.i (i=1,2, . . . , n) are identically
distributed random variables and so are .tau..sub.i (i=1,2, . . . ,
n), we have from (58) 100 E [ i = 1 n i - j i ] = n E [ 1 ] E [ - j
1 ] ( 59 )
[0252] Substituting (59) into (57) yields 101 E [ X ( j ) ] = n = 0
.infin. P { N } ( n , T ) n E [ 1 ] [ - .infin. .infin. - j p 1 ( )
( ) ] 0 y ( u ) - j u u ( 60 )
[0253] Substituting (51) and (52) into (60) yields 102 E [ X ( j )
] = n = 0 .infin. [ ( T ) n - T n ! ] n E [ 1 ] [ 0 T ( 1 T ) - j ]
0 y ( u ) - j u u = n = 1 .infin. ( T ) n - 1 - T ( n - 1 ) ! ( T )
E [ 1 ] { - T - 1 - j T } 0 y ( u ) - j u u = E [ 1 ] ( 1 - - j T )
j 0 y ( u ) - j u u ( 61 )
[0254] where Taylor expansion of e.sup..lambda.T has been used.
[0255] Average Power Density of x(t) in [0,T+.DELTA..tau.] and its
Expectation
[0256] The average power density of the force signal in (50) is
defined as 103 S 1 ( ) = X ( j ) X * ( j ) T + ( 62 )
[0257] where X* (j.omega.) is the complex conjugate of X(j.omega.):
104 X * ( j ) = i = 1 N ( T ) i j i 0 y ( u ) j u u ( 63 )
[0258] Substituting (55) and (53) into (62) yields 105 S 1 ( ) = m
= 1 N ( T ) m - m 0 y ( u ) - j u u i = 1 N ( T ) i j i 0 y ( v ) j
u v T + = ( 0 0 y ( u ) y ( v ) - j ( u - v ) u v ) m = 1 N ( T ) i
= 1 N ( T ) m i j ( i - m ) T + ( 64 )
[0259] Let k=u-v, then u=k+.nu. and du=dk. The double integral in
(64) becomes 106 - 0 y ( u ) y ( v ) - j ( u - v ) u v = 0 v - v 0
y ( v + k ) y ( v ) - j k k + 0 v 0 - v + y ( v + k ) y ( v ) - j k
k ( 65 )
[0260] Interchanging the order of integration in (65) yields 107 0
0 y ( u ) y ( v ) - j ( u - v ) u v = - 0 k - k y ( v + k ) y ( v )
- j k v + 0 k 0 - k y ( v + k ) y ( v ) - j k v ( 66 )
[0261] Let .nu.+k=.gamma., then .nu.=.gamma.-k and d.gamma.=d.nu..
The first integral on the right-hand side of (66) becomes 108 - 0 k
- k y ( v + k ) y ( v ) - j k v = - 0 k 0 + k y ( ) y ( - k ) - j k
( 67 )
[0262] Changing .gamma. in (67) back to .nu. and substituting the
resulting expression into (66) yields 109 0 0 y ( u ) y ( v ) - j (
u - v ) u v = - 0 0 - k y ( v + k ) y ( v ) v - j k k ( 68 )
[0263] Noting that 110 0 - k y ( v + k ) y ( v ) v
[0264] is an even function of k, we have 111 0 - k y ( v + k ) y (
v ) - j k k = 0 - k y ( v + k ) y ( v ) cos ( k ) k ( 69 )
[0265] Substituting (69) into (68) and substituting the resulting
expression into (64) yields 112 S 1 ( ) = 1 T + ( - 0 - k y ( v + k
) y ( v ) v cos k k ) m = 1 N ( T ) i = 1 N ( T ) m i j ( i - m ) (
70 )
[0266] Since the double integral in (70) is deterministic, we have
113 E [ S 1 ( ) ] = 1 T + ( - 0 - k y ( v + k ) y ( v ) v cos k k )
E [ m = 1 N ( T ) i = 1 N ( T ) m i j ( i - m ) ] ( 71 )
[0267] Since 114 E { m = 1 N ( T ) i = 1 N ( T ) m i sin [ j ( i -
m ) ] } = 0 , i ( i = 1 , 2 , , N ( T ) )
[0268] are identically distributed random variables, and so are
.tau..sub.i (i=1, 2, . . . , N(T)), we have 115 E [ m = 1 N ( T ) i
= 1 N ( T ) m i j ( i - m ) ] = E [ m = 1 N ( T ) i = 1 N ( T ) m i
cos ( i - m ) ] = E [ i = 1 N ( T ) i 2 + m = 1 N ( T ) i = 1 N ( T
) m i cos ( i - m ) m i ] E { N ( T ) i 2 + [ N 2 ( T ) - N ( T ) ]
i 2 cos ( 1 - 2 ) } ( 72 )
[0269] Since .tau..sub.1 and .tau..sub.2 are independently,
uniformly distributed random variables in [0,T], the probability
density function of .tau..sub.1-.tau..sub.2 is 116 p 1 - 2 ( ) = {
T - T 2 , - T T 0 , elsewhere ( 73 )
[0270] Using (51) and (73) in (72) yields 117 E [ m = 1 N ( T ) i =
1 N ( T ) m i j ( i - m ) ] = n = 0 .infin. P { N } ( n , T ) n E [
1 2 ] + n = 0 .infin. P { N } ( n , T ) ( n 2 - n ) E [ 1 2 ] -
.infin. .infin. cos ( ) p 1 - 2 ( ) = T E [ 1 2 ] + 2 2 E 2 [ 1 ] 1
- cos T 2 ( 74 )
[0271] Substituting (74) into (71) yields 118 E [ S 1 ( ) ] = 1 T +
[ - 0 - k y ( v + k ) y ( v ) v cos ( k ) k ] .times. { 2 2 E 2 [ 1
] 1 - cos ( T ) 2 + T E [ 1 2 ] } ( 75 )
[0272] Mean and Autocorrelation Functions of x(t)
[0273] The first-order cumulant function of x(t) in (50),
.kappa..sub.1[x(t)], is equal to its mean function, E[x(t)]. The
second-order cumulant function of x(t),
.kappa..sub.2[x(t.sub.1)x(t.sub.2- )], is related to its
autocorrelation function, E[x(t.sub.1).times.(t.sub.- 2 )]
through
E[x(t.sub.1)x(t.sub.2)]=.kappa..sub.2[x(t.sub.1)x(t.sub.2)]+.kappa..sub.1[-
x(t.sub.1)].kappa..sub.1[x(t.sub.2)] (76)
[0274] where t.sub.1 and t.sub.2 are any two time instants in
[0,T+.DELTA..tau.]. Following the derivations, the first- and
second-order cumulant functions of x(t) are
.kappa..sub.1[x(t)]=E[x(t)]=.lambda.E[.psi..sub.1].intg..sub.0.sup.T
y(t-.alpha.)d.alpha. (77)
.kappa..sub.2[x(t.sub.1)x(t.sub.2)].kappa..sub.xx(t.sub.1,t.sub.2)=.lambda-
.E[.psi..sub.1.sup.2].intg..sub.0.sup.Ty(t.sub.1-.alpha.)y(t.sub.2-.alpha.-
)d.alpha. (78)
[0275] where t.epsilon.[0,T+.DELTA..tau.]. Let t-.alpha.=u in (75),
then d.alpha.=-du. We have from (77) 119 E [ x ( t ) ] = E [ 1 ] t
- T t y ( u ) u ( 79 )
[0276] Let
W(t)=.intg..sub.0.sup.ty(u)du (80)
[0277] Noting that y(u)=0 when u<0 and u>.DELTA..tau., we
have from (79) 120 E [ x ( t ) ] = { E [ 1 ] W ( t ) , 0 < t
< E [ 1 ] W ( ) , t T E [ 1 ] [ W ( ) - W ( t - T ) ] , T < t
T + ( 81 )
[0278] where T>.DELTA..tau. is assumed. Let t.sub.1-.alpha.=u
and t.sub.2-t.sub.1=k in (78), then d.alpha.=-du and
t.sub.2-.alpha.=u+k. We have from (78) 121 xx ( t 1 , t 2 ) = E [ 1
2 ] t 1 - T t 1 y ( u ) y ( u + k ) u ( 82 )
[0279] When .vertline.k.vertline.>.DELTA..tau., y(u)y(u+k)=0 and
hence .kappa..sub.xx(t.sub.1, t.sub.2)=0. When
0.ltoreq.k.ltoreq..DELTA..tau., we have from (82) for different
t.sub.1 122 xx ( t 1 , t 2 ) = { E [ 1 2 ] 0 t 1 y ( u ) y ( u + k
) u , 0 t 1 < - k E [ 1 2 ] 0 - k y ( u ) y ( u + k ) u , - k t
1 T E [ 1 2 ] t 1 - T - k y ( u ) y ( u + k ) u , T < t 1 T + -
k 0 , T + - k < t 1 T + ( 83 )
[0280] When -.DELTA..tau..ltoreq.k<0, we have from (82) for
different t.sub.1 123 xx ( t 1 , t 2 ) = { E [ 1 2 ] - k t 1 y ( u
) y ( u + k ) u , - k t 1 < E [ 1 2 ] - k y ( u ) y ( u + k ) u
, t 1 T - k E [ 1 2 ] t 1 - T y ( u ) y ( u + k ) u , T - k < t
1 T + 0 , 0 t 1 < - k ( 84 )
[0281] Let u+k=.nu. in (84), then u=.nu.-k and du=d.nu.. We have
from (84) after changing .nu. back to u 124 xx ( t 1 , t 2 ) = { E
[ 1 2 ] 0 t 1 + k y ( u - k ) y ( u ) u , - k t 1 < E [ 1 2 ] 0
+ k y ( u - k ) y ( u ) u , t 1 T - k E [ 1 2 ] t 1 - T + k + k y (
u - k ) y ( u ) u , T - k < t 1 T + 0 , 0 t 1 < - k ( 85
)
[0282] Combining the second equations in (83) and (85), we have for
t.sub.1 and t.sub.2 in [.DELTA..tau.,T] 125 xx ( t 1 , t 2 ) = E [
1 2 ] 0 - k y ( u ) y ( u + k ) u ( 86 )
[0283] Since by the second equation in Eq. (81), E[x(t)] is a
constant for t .epsilon.[.DELTA..tau., T], and by (86),
.kappa..sub.xx(t.sub.1, t.sub.2) is a function of k=t.sub.2-t.sub.1
for t.sub.1 and t.sub.2 in [.DELTA..tau., T], x(t) is a wide-sense
stationary random process in [.DELTA..tau., T]. Substituting the
second equation in (81) and (86) into (76) yields the
autocorrelation function for t.sub.1 and t.sub.2 in
[.DELTA..tau.,T] 126 R xx ( k ) .cndot. E [ x ( t 1 ) x ( t 2 ) ] =
E [ 1 2 ] 0 - k y ( u ) y ( u + k ) u + 2 E 2 [ 1 ] W 2 ( ) ( 87
)
[0284] Fourier Transform of the Mean Function of x(t) and its
Equivalence to E[X(j.omega.)]
[0285] Applying the Fourier transform to E[x(t)] in Eq. (81) yields
127 F { E [ x ( t ) ] } = E [ 1 ] { 0 W ( t ) - j t t + T W ( ) - j
t t + T T + [ W ( ) - W ( t - ) ] - j t t } ( 88 )
[0286] The three integrals in (88) are referred to as I.sub.1,
I.sub.2, and I.sub.3, respectively. Consider first the third
integral in (88) 128 I 3 T T + [ W ( ) - W ( t - T ) ] - j t t ( 89
)
[0287] Let t-T=.theta. in (88), then t=T+.theta. and dt=d.theta..
We have from Eq. (89)
I.sub.3=e.sup.-j.omega.T.intg..sub.0.sup..DELTA..tau.W(.DELTA..tau.)e.sup.-
-j.omega..theta.d.theta.-e.sup.-j.omega.T.intg..sub.0.sup..DELTA..tau.W(.t-
heta.)e.sup.-j.omega..theta.d.theta. (90)
[0288] Changing .theta. in (90) back to t and combining I.sub.1 and
I.sub.3 yields
I.sub.1+I.sub.3=e.sup.-j.omega.T.intg..sub.0.sup..DELTA..tau.W(.DELTA..tau-
.)e.sup.-j.omega.tdt+(1-e.sup.-j.omega.T).intg..sub.0.sup..DELTA..tau.W(t)-
e.sup.-j.omega.tdt (91)
[0289] Adding I.sub.2 to (91), simplifying the expression, and
substituting it into (88) yields 129 F { E [ x ( t ) ] } = E [ 1 ]
[ ( 1 - - j t ) 0 W ( t ) - j t t + - j t 0 W ( ) - j t t + T W ( )
- j t t ] = E [ 1 ] ( 1 - - j t ) W ( ) [ 1 j + 0 ( W ( t ) W ( ) -
1 ) - j t t ] ( 92 )
[0290] We will show that F{E[x(t)]} in (92) is equivalent to
E[X(j.omega.)] in (91). By (80), we have 130 W ( t ) = { 0 , t ( -
.infin. , 0 ) 0 t I ( t ) t , t [ 0 , ] W ( ) , t ( , .infin. ) (
93 )
[0291] The integral in (61) can be written as 131 0 y ( u ) - j u u
= 0 y ( u ) u + 0 y ( u ) ( - j u - 1 ) u = 0 y ( u ) u - j 0 y ( u
) 0 u - j v v u ( 94 )
[0292] Interchanging the order of integration in the double
integral in (94) and using (93) yields 132 0 y ( u ) - j u u = 0 y
( u ) u - j 0 v 0 y ( u ) - j v u = W ( ) { 1 + j 0 [ W ( v ) W ( )
- 1 ] - j v v } ( 95 )
[0293] Substituting (95) into (61) yields (91). This shows that the
expectation E and the Fourier transform F are commutative as both
are linear operators.
[0294] Equation (92) consists of two parts: the first part, 133 E [
1 ] ( 1 - - j t ) W ( ) 1 j ,
[0295] is the Fourier transform of the stationary part of E[x(t)]
in [.DELTA..tau.,T], and the second part, 134 E [ 1 ] ( 1 - - j t )
W ( ) 0 ( W ( t ) W ( ) - 1 ) - j t t ,
[0296] is the sum of the Fourier transforms of the nonstationary
parts of E[x(t)] in [0,.DELTA..tau.] and [T,T+.DELTA..tau.]. When
.DELTA..tau..fwdarw.0, since 135 W ( t ) W ( ) - 1
[0297] is finite, 136 0 ( W ( t ) W ( ) - 1 ) - j t t ,
[0298] and consequently the second part of E[X(j.omega.))],
approaches zero.
[0299] Average Power Density of x(t) in [.DELTA..tau.,T], Its
Expectation, and Power Spectral Density
[0300] Since x(t) is stationary in [.DELTA..tau.,T], the average
power density of x(t) in [.DELTA..tau.,T] is defined as 137 S 2 ( )
= X s ( j ) X s * ( j ) T - ( 96 )
[0301] where X.sub.s
(j.omega.)=.intg..sub..DELTA..tau..sup.Tx(t)e.sup.-j.- omega.tdt.
Taking the expectation of (96) yields 138 E [ S 2 ( ) ] = 1 T - T T
E [ x ( t 1 ) x ( t 2 ) ] - j ( t 1 - t 2 ) t 1 t 2 ( 97 )
[0302] Let k=t.sub.2-t.sub.1 in (97), then dk=dt.sub.2. Since
E[x(t.sub.1)x(t.sub.2)]=R.sub.xx(k), (97) becomes after
interchanging the order of integration 139 E [ S 2 ( ) ] = 1 T - T
t 1 - t 1 T - t 1 R xx ( k ) j k k = - ( T - ) T - R xx ( k ) ( 1 -
k T - ) j k k ( 98 )
[0303] Substituting (87) into (98) yields 140 E [ S 2 ( ) ] = E [ 1
2 ] - ( T - ) T - 0 - k y ( u + k ) y ( u ) u ( 1 - k T - ) j k k +
2 2 E 2 [ 1 ] W 2 ( ) - ( T - ) T - ( 1 - k T - ) j k k ( 99 )
[0304] Noting that y(u)y(u+.vertline.k.vertline.)=0 when
.vertline.k.vertline.>.DELTA..tau. and 141 0 - k y ( u + k ) y (
u ) u ( 1 - k T - )
[0305] an even function of .vertline.k.vertline., we have from (99)
142 E [ S 2 ( ) ] = 2 2 E 2 [ 1 ] W 2 ( ) 1 - cos ( T - ) 2 ( T - )
+ E [ 1 2 ] - 0 - k y ( u + k ) y ( u ) u ( 1 - k T - ) cos k k (
100 )
[0306] where T>2.DELTA..tau. has been assumed.
[0307] The power spectral density of x(t) can be obtained from
(100) by increasing T to infinity 143 S x ( ) = lim T .infin. S 2 (
) = 2 2 E 2 [ 1 ] W 2 ( ) ( ) + E [ 1 2 ] - 0 - k I ( u + k ) I ( u
) cos ( k ) k ( 101 )
[0308] where we have used 144 lim T .infin. 1 - cos ( T - ) 2 ( T -
) = lim T .infin. 2 sin 2 ( T - ) 2 2 ( T - ) = ( ) ( 102 )
[0309] in which .delta.() is the Dirac delta function. The power
spectral density can also be obtained from (98) by increasing T to
infinity
S.sub.x(.omega.)=lim/T.fwdarw..infin.E[S.sub.2(.omega.)]=.intg..sub.-.infi-
n..sup..infin.R.sub.xx(k)e.sup.j.omega.kdk=.intg..sub.-.infin..sup..infin.-
R.sub.xx(k)e.sup.-j.omega.kdk=F[R.sub.xx(k)] (103)
[0310] where R.sub.xx(-k)=R.sub.xx(k) has been used. Equation (103)
is the well-known Wiener-Khintchine theorem, which states the power
spectral density is the Fourier transform of the autocorrelation
function. Substituting (87) into (83), and noting that
I(u)I(u+.vertline.k.vertline- .)=0 when
.vertline.k.vertline.>.DELTA..tau. and the Fourier transform of
1 is 2.pi..delta.(.omega.), yields (101). Note that the power
spectral density is only defined for a wide-sense stationary
process with an infinite time record. When the mean amplitude of
each pulse E[.psi..sub.1] is not equal to zero, there is an
associated delta function in the power spectral density.
[0311] Comparison of E[S.sub.1(.omega.))] and
E[S.sub.2(.omega.)]
[0312] By (100), E[S.sub.2(.omega.)] consists of two parts: the
first part, 145 2 2 E 2 [ 1 ] W 2 ( ) 1 - cos ( T - ) 2 ( T - )
,
[0313] which depends on the arrival rate .lambda., the mean
amplitude of each pulse E[.psi..sub.1], the total area of the
normalized shape function W(.DELTA..tau.), and the time length
T-.DELTA..tau., describes the average effects of the stationary
part of x(t) in [.DELTA..tau.,T] and is referred to here as the
first-order statistical power density, and the second part, 146 E [
1 2 ] - 0 - k y ( u + k ) y ( u ) u ( 1 - k T - ) cos k k ,
[0314] which depends on the mean square amplitude of each pulse
E[.psi..sub.1.sup.2] and the shape function y() in addition to
.lambda., T, and .DELTA..tau., describes the variational effects of
the stationary part of x(t) and is referred to here as the
second-order statistical power density. While E[S.sub.1(.omega.)]
in (75) consists of two parts, the first part, 147 2 2 E 2 [ 1 ] T
+ [ - 0 - k y ( v + k ) y ( v ) v cos ( k ) k ] 1 - cos ( T ) 2
,
[0315] depends also on the shape function y() as the shape function
is used in calculating the power associated with the nonstationary
parts of x(t) in [0,.DELTA..tau.] and [T, T+.DELTA..tau.].
[0316] When the shape function is a delta function (i.e.,
.DELTA..tau..fwdarw.0), the nonstationary parts of x(t) vanish and
E[S.sub.1(.omega.)] in (75) can be shown to be equivalent to
E[S.sub.2(.omega.)] in (100). By (68) and (69), we have 148 - 0 - k
y ( v + k ) y ( v ) v cos k k = 0 y ( v ) - j v v 2 ( 104 )
[0317] When y(.nu.)=W.sub.0.delta.(.nu.), where W.sub.0 is a
constant, we have from (104) 149 lim -> 0 - 0 - k y ( v + k ) y
( v ) v cos k k = lim -> 0 0 W 0 ( v ) - j v v 2 = W 0 2 ( 105
)
[0318] Since .vertline.k.vertline.<.DELTA..tau. in (100), we
have .vertline.k.vertline..fwdarw.0 as .DELTA..tau..fwdarw.0 and
hence 150 lim -> 0 ( 1 - k T - ) = 1 ( 106 )
[0319] Substituting (105) and (106) into (100) and noting that
W(.DELTA..tau.)=W.sub.0 when y(.nu.)=W.sub.0.delta.(.nu.), and
substituting (105) into (75) and noting that
T+.DELTA..tau..fwdarw.T as .DELTA..tau..fwdarw.0, yields 151 E [ S
2 ( ) ] = E [ S 1 ( ) ] = 2 2 E 2 [ 1 ] 1 - cos T 2 T W 0 2 + E [ 1
2 ] W 0 2 ( 107 )
[0320] By (104) the power spectral density in (101) can be
expressed as
S.sub.x(.omega.)=2.pi..lambda..sup.2E.sup.2[.psi..sub.1]W.sup.2(.DELTA..ta-
u.).delta.(.omega.)+.lambda.E[.psi..sub.1.sup.2].vertline.Y(j.omega.).vert-
line..sup.2 (108)
[0321] where
Y(j.omega.)=.intg..sub.0.sup..DELTA..tau.y(t)e.sup.-j.omega.t- .
When T.fwdarw..infin. in (75), we have by using (104)
S(.omega.)=lim/T.fwdarw..infin.E[S.sub.1(.omega.)]=2.pi..lambda..sup.2E.su-
p.2[.psi..sub.1].vertline.Y(j.omega.).vertline..sup.2.delta.(.omega.)+.lam-
bda.E[.psi..sub.1.sup.2].vertline.Y(j.omega.).vertline..sup.2
(109)
[0322] Equation (109) has a slightly different form from that of
(108) because the power associated with the nonstationary parts of
x(t) is included in (109). When y(t)=W.sub.0.delta.(t) (108) and
(109) reduce to
S.sub.x(.omega.)=S(.omega.)=2.pi..lambda..sup.2E.sup.2[.psi..sub.1]W.sub.0-
.sup.2.delta.(.omega.)+.lambda.E[.psi..sub.1.sup.2]W.sub.0.sup.2
(110)
[0323] Equation (110) can also be obtained from (106) by letting
T.fwdarw..infin. and using 152 lim T -> .infin. 1 - cos T 2 T =
lim T -> .infin. 2 sin 2 T 2 2 T = ( ) ( 111 )
EXAMPLES AND NUMERICAL SIMULATION
[0324] When the shape function of the pulses is represented by a
half sine wave, i.e., 153 y ( t ) = sin ( t ) [ H ( t ) - H ( t - )
] ( 112 )
[0325] where H() is the Heaviside function, we obtain by using
(60), (64), and (103) 154 E [ X ( j ) ] = - E ( 1 ) ( 1 + - j ) ( -
jT - 1 ) j ( 2 - 2 2 ( 113 ) E [ S 1 ( ) ] = 2 2 2 [ 1 + cos ( ) ]
( 2 2 - 2 ) 2 ( T + ) [ TE ( 1 2 ) - 2 2 E 2 ( 1 ) ( cos T - 1 ) 2
] ( 114 ) E [ S 2 ( ) ] = E ( 1 2 ) 2 { [ - 2 T 4 - 2 T cos ( ) 4 +
( 4 cos ( ) 4 + 4 ) ] ( 2 2 - 2 ) 3 ( T - ) + [ + ( 8 sin ( ) 2 + 2
T cos ( ) 2 2 + 2 T 2 2 ) 2 ] ( 2 2 - 2 ) 3 ( T - ) + [ ( 4 5 - 4
cos ( ) 2 2 - 2 2 2 ) 3 ] ( 2 2 - 2 ) 3 ( T - ) } + 8 E 2 ( 1 ) 2 (
1 - cos ( T - ) ) 2 2 ( T - ) ( 115 )
[0326] Consider next the normalized shape function y(t) shown in
FIG. 35 with unit maximum amplitude. It is obtained by averaging a
series of normalized force pulses from impact tests on the four-bay
space frame as shown in FIG. 10. There are 21 sample points in the
shape function, which are connected, as shown in FIG. 34. Other
parameters used are T=8 s, .DELTA..tau.=20.times.T/1024=0.15625 s
where h=T/1024=0.0078125/s is the sampling interval,
.lambda.=4.14/s, E[.psi..sub.1]=0.8239N, and
E[.psi..sub.1.sup.2]=0.7163N.sup.2. The curve for E[x(t)] in the
time interval from 0 to 8.15625 s, shown as a solid line in FIG.
36, is calculated using (81). It is seen that E[x(t)] increases
from 0 to 0.0352 N in the first 0.15625 s, remains at 0.0352 N from
0.15625 s to 8 s, and decreases from 0.0352 N to 0 in the last
0.15625 s. The curve for 20 log.vertline.E[X(j.omega.)] in the
frequency range from 0 to 50 Hz, shown as a dotted line in FIG. 37,
is calculated using (61). The curve for 10
log{E[S.sub.1(j.omega.))]} in the same frequency range, shown as a
solid line in FIG. 38, is calculated using (75) . The curve for 10
log{E[S.sub.1(j.omega.)]} with .lambda.=1/s and the other
parameters unchanged is shown as a dashed line in FIG. 38. It is
seen that E[S.sub.1(j.omega.)] increases by 4.14 to 15.6247 times
in the frequency range shown when .lambda. is increased from
.lambda..sub.2=1/s to .lambda..sub.1=4.14/s. This result can be
shown by using (75) 155 lim -> .infin. E [ S 1 ( j ) ] = 1 lim
-> .infin. E [ S 1 ( j ) ] = 2 = 1 2 E 2 [ 1 ] T 2 + 1 TE [ 1 2
] 2 2 E 2 [ 1 ] T 2 + 2 TE [ 1 2 ] = 1 2 E 2 [ 1 ] T + 1 E [ 1 2 ]
2 2 E 2 [ 1 ] T + 2 E [ 1 2 ] = 15.6247 ( 116 ) lim -> .infin. E
[ S 1 ( j ) ] = 1 lim -> .infin. E [ S 1 ( j ) ] = 2 = 1 2 =
4.14 ( 117 )
[0327] This shows that a larger arrival rate .lambda. would
increase the energy input to the structure over the entire
frequency domain.
[0328] Numerical simulation is undertaken next to validate the
analytical predictions. The random number N(T) satisfying the
Poisson distribution in (51) with .lambda.=4.14/s and T=8 s is
generated using MATLAB. Similarly, the random numbers corresponding
to the random variables .tau..sub.i (i=1, 2, . . . , N(T)),
satisfying the uniform distribution in (52), and the random numbers
corresponding to the random variables .psi..sub.i (i=1, 2, . . . ,
N(T)), satisfying the Gaussian distribution in (53) with
.mu.=0.8239 N and .sigma..sup.2=0.7163-0.82392 N.sup.2=0.0375
N.sup.2, are generated. Using the shape function constructed
earlier, a sample function of x(t) in (50) at time t=rh, where r=0,
1, . . . . , 156 R - 1 ( R = T + h = 1044 ) ,
[0329] denoted by x.sub.r, can be obtained. The discrete Fourier
transform (DFT) of the time series {x.sub.r} is calculated using
MATLAB. The DFT of the series {x.sub.r} is defined by 157 X q = 1 R
r = 0 R - 1 x r - j 2 q r R ( 118 )
[0330] where q=0, 1, . . . , R-1. Equation (118) is an approximate
formula for calculating the coefficients of the Fourier series of a
periodic function whose values in the period [0, T+.DELTA..tau.]
are given by those of x(t): 158 X q = 1 T + 0 T + x ( t ) - j 2 q t
T + t ( 119 )
[0331] The Fourier components X.sub.q correspond to harmonics of
frequency 159 q = 2 q T + .
[0332] Recall that the Fourier transform of x(t) in (50) is given
by 160 X ( j ) = 0 T + x ( t ) - j t t ( 120 )
[0333] Let 161 = q = 2 q T +
[0334] in (120) and compare the resulting expression with (119), we
find that X.sub.q in (118) multiplied by T+.DELTA..tau. provides an
approximate value of X(j.omega.) at frequency .omega..sub.q.
Similarly, X.sub.qX*.sub.q multiplied by T+.DELTA..tau. provides an
approximate value of 162 S 1 ( ) = X ( j ) X * ( j ) T +
[0335] at frequency .omega..sub.q. By averaging 5000 sample series
of {x.sub.r}, we obtain the curve for E[x(t)], shown as a dotted
line in FIG. 36. By averaging 1000 sample series of {20
log[(T+.DELTA..tau.).vert- line.X.sub.q.vertline.]}, we obtain the
curve for 20 log.vertline.E[X(j.omega.)].vertline., shown as a
solid line in FIG. 37. By averaging 100 sample series of {10
log[(T+.DELTA..tau.).vertline.X.sub- .qX*.sub.q.vertline.]}, we
obtain the curve for 10 log{E[S.sub.1(j.omega.)- ]}, shown as a
dotted line in FIG. 38. The numerical results are in good agreement
with the analytical ones.
[0336] The stochastic model was experimentally validated for an
experimenter conducting manually a random series of impacts on the
four bay space frame as shown in FIG. 10. One hundred ensemble
averages were used. The experimental probability density functions
of the arrival time, the number of arrived pulses, and the pulse
amplitudes are in good agreement with the analytical values, as
shown in FIGS. 39-41.
[0337] Thus, the system and method for detecting structural damage
and the random impact series method as embodied and broadly
described herein can be applied to an unlimited number and type of
structures to provide automated, reliable damage detection and
assessment and to conduct modal testing. This system could be
further automated to conduct periodic tests and provide results to
a centralized monitoring section. Regular health monitoring of
these types of structures could provide additional protection
against potential failure, as well as a characterization of usage
and wear over time in particular environmental conditions for
predicting useful service life.
[0338] The foregoing embodiments and advantages are merely
exemplary and are not to be construed as limiting the present
invention. The present teaching can be readily applied to other
types of systems. The description of the present invention is
intended to be illustrative, and not to limit the scope of the
claims. Many alternatives, modifications, and variations will be
apparent to those skilled in the art. In the claims,
means-plus-function clauses are intended to cover the structures
described herein as performing the recited function and not only
structural equivalents but also equivalent structures.
* * * * *