- Email: [email protected]
Computer simulation of structural random loading identification
Computers and Structures 79 (2001) 375±387
www.elsevier.com/locate/compstruc
Computer simulation of structural random loading identi®cation J.H. Lin a,1, X.L. Guo a, H. Zhi a, W.P. Howson b,*, F.W. Williams b a
State Key Laboratory of Structural Analysis of Industrial Equipment, Dalian University of Technology, Dalian 116023, People's Republic of China b Division of Structural Engineering, Cardi School of Engineering, Cardi University, Cardi CF24 3TB, UK Received 13 August 1999; accepted 27 May 2000
Abstract Loading identi®cation has long been a dicult problem for multi-input±multi-output stationary random vibration problems. In this paper, the inverse pseudo-excitation method is developed for dealing with such identi®cation problems. Numerical simulations are used to reveal the diculties in such problems, and some numerical analysis measurements are suggested, which may be very useful in the setting up and processing of experimental data so as to obtain reasonable predictions of the input loading from selected structural responses, e.g. a few structural displacements or strains, or a combination of them. It is shown that computer simulation is a very cheap way of overcoming the problem of which responses to measure on real structures in order to identify the power spectral density matrix of the applied forces. Ó 2000 Elsevier Science Ltd. All rights reserved. Keywords: Random vibration; Loading identi®cation; Computer simulation
1. Introduction For linear structures subjected to multiple stationary random excitations, which may be partially coherent, the direct problem, i.e. the computation of the power spectral densities (PSDs) of the various responses from the given excitation PSD, has been solved [1±3], in particular by using the exact and ecient pseudoexcitation method (PEM) [4±6]. There are two types of inverse problems (also called back analysis problems) corresponding to this direct problem. The ®rst is the socalled structure identi®cation problem, for which the loading and response PSD functions are all known and are used to identify the properties of the structure. Many publications cover such problems and their successful
* Corresponding author. Address: Division of Structural Engineering, Cardi School of Engineering, Cardi University, Cardi CF24 3TB, UK. 1 Visiting Professor at the University of Wales Cardi.
use in engineering [2,3,7±10]. The second is the so-called loading identi®cation problem, for which the response information and structural properties are known and are used to identify information about the loading. For this second type of inverse problem, if there is only one excitation acting on the structure, the problem is quite straightforward and is not discussed in this paper. For multi-input±multi-output (MIMO) problems, however, its solution is not easy and it is dicult to ®nd papers in this area which demonstrate its diculty. To the authorsÕ knowledge, some such work has been done at high expense, but has nevertheless yielded identi®cation results that are unacceptable because of their very poor precision. However, such inverse problems are of great importance in practical engineering because sometimes the direct measurement of the applied loads may be very dicult, whereas the measurement of some of the responses may be relatively easy. For any direct problem for which the excitation PSD matrix has been given, the PSD matrix of any responses can be computed quickly and precisely by means of the
0045-7949/01/$ - see front matter Ó 2000 Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 4 9 ( 0 0 ) 0 0 1 5 4 - 1
376
J.H. Lin et al. / Computers and Structures 79 (2001) 375±387
PEM [4±6]. In this paper, the PEM is modi®ed to generate an inverse pseudo-excitation method (IPEM), so that the excitation PSD matrix can be obtained from a set of known response PSD functions which can be displacements, strains, or their combination. Two examples of computer simulation are given, namely a plane truss and a three dimensional oshore jacket platform built in the Bohai sea in Northern China. These computer simulation examples show that the extended IPEM has good noise-resistance characteristics and so may be quite useful for practical engineering problems involving large FEM models.
When all excitations are fully coherent, i.e. r 1, Eq. (7) takes the following simple form:
2. Pseudo-excitation method for solving to ®nd the stationary random responses of structures
Sxx H Syy H T
The conventional formula for solving the stationary random responses of linear structures in the frequency domain is Syy H Sxx H T
1
in which Sxx is the known excitation PSD matrix, H is the transfer function matrix, Syy is the response PSD matrix to be computed and the superscripts and T represent complex conjugate and transpose, respectively. The excitation PSD matrix Sxx can generally be assumed to be a p p matrix with rank r ( 6 p). Since Sxx must be a Hermitian matrix [3,4], it can be decomposed into Sxx
r X fagj fagTj :
2
j1
Syy H Sxx H T H fag fagT HT fbg fbgT :
8
3. The inverse pseudo-excitation method for structural random loading identi®cation Let Eq. (1) be written in the form [2] 9
in which the superscript represents generalized inversion (performed as described beneath Eq. (29)), since H is in general not a square matrix. Provided that the responses of the structure are caused entirely by the random loadings to be identi®ed, then Eq. (9) gives the algorithm for computing Sxx from the known PSD matrix Syy . Clearly, Eq. (9) is the backward extension of the PEM described above. It is therefore, called the IPEM. For complicated engineering problems, direct use of Eq. (9) is quite inecient, and so it should be dealt with analogously to the last section, as follows: Pr Using Eq.T (7), Syy can be decomposed into ixt j1 fbgj fbgj . The pseudo-response f y gj fbgj e then follows from Eq. (5), and applying Eq. (6) inversely, gives the pseudo-excitations f xgj of Eq. (3) as f xgj fagj eixt H fbgj eixt :
10
Introducing the pseudo-excitation fxgj fagj eixt
j 1; 2; . . . ; r
3
Hence Eqs. (9), (7) and (10) enable the excitation PSD matrix Sxx to be given in terms of the fagj as
4
Sxx H Syy H T H
leads to the harmonic response fyg given by [5] fygj H fxgj ; so that, using Eq. (3), fygj fbgj eixt ;
5
where fbgj H fagj :
6
Hence, using Eqs. (1), (2) and (4), the response PSD matrix Syy is given in terms of the fbgj as ! r X T T Syy H Sxx H H fagj fagj H T j1
r X j1
fbgj fbgTj :
7
r X fbgj fbgTj H T j1
r X fagj fagTj :
11
j1
When all excitations are fully coherent, i.e. r 1, Eq. (11) takes the following simple form: Sxx H Syy H T H fbg fbgT H T fag fagT :
12
For large systems, the direct use of Eq. (12) or Eq. (11) is still computationally very expensive. Therefore, the following equation reduction scheme should be used, which is based on the mode-superposition method.
J.H. Lin et al. / Computers and Structures 79 (2001) 375±387
4. Excitation power spectral density identi®cation for large systems The equations of motion of a large system can be expressed as MfY g CfY_ g KfY g fF g
13
in which M, C and K are, respectively, the mass, damping and stiness matrices, dots denote dierentiation with respect to time and fY g and fF g are, respectively, the displacement and applied force vectors. Even though the order n of fF g and fY g is very high, the number, l, of non-zero excitations in fF g and the number, m, of measured responses in fY g are both n. Therefore, let fF g Ef ff g;
14
fyg Ey fY g;
15
where ff g is the l dimensional vector consisting of the non-zero elements in fF g, fyg is the m dimensional vector consisting of the m measured quantities in fY g and Ef and Ey are matrices which contain only zero and unity. Assume that the m m matrix Syy x is known, usually being speci®ed in the form of a series of values at discrete frequency points, so that the requirement is to ®nd the l l excitation PSD matrix Sff x by back analysis. The mode-superposition scheme is now used. Firstly, the lowest q modes, denoted by the n q matrix U, and the corresponding diagonal q q eigenvalue matrix X2 must be found by solving KU MUX2 ;
16
subject to UT M U I ; where I is a q q unit matrix:
17
Reducing Eq. (13) by means of the modal matrix U gives fY g UfUg;
18
g C0 fU_ g X2 fU g UT Ef ff g fU
19
in which C0 UT CU:
20
For proportionally damped systems, C 0 is a diagonal matrix. In particular, if its damping ratios 1j j 1; 2; . . . ; q are all known, then the jth diagonal element of C 0 will be 21j xj , where x2j is the jth diagonal element of X2 . For non-proportional damping matrices, C 0 is non-diagonal. However, for both cases, the following derivation applies, in which partial coherency is assumed for generality.
377
Substituting Eq. (18) into Eq. (15) gives fyg Ey UfU g:
21
Then, decomposing Syy x by using Eq. (7), and forming the pseudo-responses according to Eq. (5), gives f y gj fbgj eixt , which from Eq. (21) can be regarded as being produced by the pseudo-ÔexcitationÕ, fU gj fcgj eixt :
22
This fU gj can be considered to be generated by the pseudo-excitation f f gj fagj eixt , so that Eq. (19) gives fU gj H UT Ef ff gj
23
in which H ÿx2 I X2 ixC0 ÿ1 ; 24 p where i ÿ1. Substituting Eq. (23) into Eq. (21) produces fygj Rff gj
25
in which R Ey UH UT Ef :
26
Therefore, ff gj R fygj :
27
Finally, using the IPEM gives (Eq. (11)) X Sff R Syy RT f f gj f f gTj :
28
When the number of excitations l exceeds the number of measured responses m, the m l generalized inverse matrix of R, i.e. R RT R ÿ1 RT
29
does not exist, and so Sff cannot be found from Syy . When l < m, the inversion procedure used is to solve for the l unknowns from the m equations by means of the least square method, so that only approximate solutions can be expected [11]. Clearly, when l m, the equation set can be solved exactly. Therefore, in general, the inverse problem is solvable under the conditions l 6 m and l 6 q. However, even when these conditions are met, there is no guarantee that a satisfactory excitation PSD matrix will be obtained. Practical situations are rather complex, as discussed via the following simulation example. For many engineering problems, strain measurements are more convenient and/or more precise than are displacement measurements. Assume that the vector {e} of measured strains is related to a reduced set of displacements fyg of the structural displacement vector fY g by
378
J.H. Lin et al. / Computers and Structures 79 (2001) 375±387
feg Ee fyg
30
in which the transformation matrix Ee can be calculated by standard FEM [2±4] and fyg can be computed in terms of Eq. (15). Substituting Eqs. (15) and (18) into Eq. (30) gives feg Ee Ey UfUg:
31
Assuming that p suciently long time-history curves of structural strains have been recorded, the corresponding p p PSD matrices See can easily be obtained at discrete frequency points by means of a power spectrum analyzer. Because p is small, See can be easily decomposed into the form of Eq. (2), i.e. See
r X fggj fggTj :
32
j1
where r is the rank of the p p PSD matrix See (r 6 p). The corresponding pseudo-strain responses are then fegj fggj eixt
33
and can be regarded as being caused by pseudodisplacement excitations of the form of Eq. (22). Hence, Eq. (31) gives the corresponding strain vector fegj Ee Ey UfU gj ;
34
where fU gj can be regarded as being caused by pseudoexcitation of the form ff gj fagj eixt : From Eq. (19), fU gj can be expressed in terms of f f gj by Eq. (23), which is then further substituted into Eq. (34) to produce fegj Re ff gj
35
In practical engineering, good measured data is often quite limited. Hence, sometimes the loading cannot be identi®ed adequately or at all from either the displacement or strain responses alone. Then, combined use of the displacement and strain data may make identi®cation possible, or give better identi®cation. It will be seen that the derivation and implementation of such a combined method, as described below, will have similarities to those in the last section. Consider the displacement vector fyg which contains exclusively all the measured displacements and all of the displacements associated with the measured strains. Therefore, the measured strains feg and the measured displacements fyd g can be expressed as feg Ke fyg;
39
fyd g Ed fyg
40
in which Ke is a local stiness matrix, Ed contains only ones and zeros, and both Ke and Ed can be calculated by means of standard FEM [12,13]. These two equations can be combined to give fym g Em fyg in which e fym g ; yd
41 Em
Ke : Ed
42
The use of all available ( s) measured curves of feg and fyd g enables the corresponding s s PSD matrix Sym ym to be found in discrete form by means of a power spectrum analyzer. Then decomposition gives Sym ym
in which
r X fym gj fym gTj ;
43
j1
Re Ee Ey UH UT Ef :
36
Hence, ff gj Re fegj :
37
r X ff gj ff gTj ;
where r ( 6 s) is the rank of the matrix Sym ym . In a similar manner to the derivation in the previous section, fym gj can be regarded as being caused by the pseudodisplacement fU gj fcgj eixt . Hence (Eqs. (15), (18) and (41)), fym gj Em Ey UfU gj :
Finally, f f gj is used to calculate Sff from Sff Re See Re T
5. Identi®cation of the excitation PSD matrix from the mixed displacement±strain PSD matrix
38
j1
where Re is in general not a square matrix. Hence, the inversion procedure used, and a discussion of when such inversion is possible, is the one given above, beneath Eq. (29). However, the conditions for inversion of Re to be possible are now slightly stricter than for the R of Eq. (26) because the additional condition l 6 p is required.
44
Substituting Eq. (23) into Eq. (44) gives fym gj Rm ff gj
45
in which Rm Em Ey UH UT Ef :
46
Therefore, ff gj Rm fym gj ;
47
J.H. Lin et al. / Computers and Structures 79 (2001) 375±387 r X Sff Rm Sym ym Rm T ff gj ff gTj
48
j1
in which the possibility of inverting to obtain Rm is subject to the additional condition that r P l, when compared to inverting to obtain R of Eq. (27), where l is the number of excitations to be identi®ed. Finally, the IPEM, i.e. the last expression of Eq. (38), can be used to give the excitation PSD matrix Sff .
6. Implementation procedure of the IPEM The steps of the implementation procedure are as follows: (1) Calculate the matrix R from Eq. (26), Re from Eq. (36), or Rm from Eq. (46) with H optionally formed by experimental means. (2) Calculate the response PSD matrix Syy , See or Sym ym for practical engineering problems by measurement and data processing. However, for the purpose of computer simulations of this paper, such matrices were obtained by direct analyses and then represented to very few digits accuracy (in general two) when performing the back analyses. (3) Decompose the response PSD matrix by using the last expression of the appropriate one of Eqs. (7), (32) or Eq. (43) to generate the corresponding pseudoresponses. (4) Solve Eqs. (27), (37) or Eq. (47) for all excitations f f gj [14], and then compute the excitation PSD matrix Sff from Eqs. (28), (38) or Eq. (48).
7. Computer simulation example 7.1. Example 1: nine DOF structure The stiness and mass distributions of a nine DOF structure are shown in Fig. 1. For simplicity, dimensionless quantities m 1 and k 1 are assumed. The complete set of natural angular frequencies are 1.00000, 1.03736, 1.13705, 1.27174, 1.41421, 1.54359, 1.64533, 1.70994, 1.73205 and the damping ratios are all equal
379
with a value of 0.03. Two horizontal forces are applied at nodes 3 and 6, and the corresponding PSD matrix is Sff
3:0 2:0
2:0 0:0 1:0 i 4:0 ÿ1:0 0:0
49
within the signi®cant frequency region x 2 0:8; 1:8: This region is discretized with an interval of Dx 0:01. Because of the existence of noise disturbance, practical responses cannot be measured very precisely. Because of this, only two or three digits of the response PSD results generated from the direct random vibration analysis [4] were taken as the ``measured values'' for the computer simulation. Some observations of interest follow in the following six sub-sections. All results presented were obtained by the IPEM and are curves of the auto-PSD of the excitations f3 and f6 , i.e. Sf3 f3 and Sf6 f6 , and their cross-PSD Sf3 f6 , which has real and imaginary parts Sr;f3 f6 and Si;f3 f6 , respectively. 7.1.1. The eect of the precision of measurement of the response PSD matrix on the precision with which the excitation PSD matrix is identi®ed Because noise always exists in the measurements, only a limited number of digits of the measured response values are reliable for use in the loading identi®cation. When only the ®rst two eective digits of all measured values were used, the results obtained for all damping ratios (0.03) are plotted in Fig. 2(a). It can be seen that the identi®cation values of the PSD functions ¯uctuate around their exact values, Eq. (49), of 3.0, 4.0, 2.0 and 1.0. The maximum error (taken directly from the computed results) is 6.5% for Sr;f3 f6 1:869 when x 1:00. In contrast, if the measurement precision is so high that the ®rst three eective digits of all the measured quantities are reliable, Fig. 2(b) replaces Fig. 2(a), and obviously, the four curves obtained are very close to their theoretical (i.e. exact) values. 7.1.2. The eect of structural damping on the identi®cation precision When the damping ratios were all changed from 0.03 to 0.01, Fig. 2(a) and (b) becomes Fig. 2(c) and (d). Hence, it can be seen that higher identi®cation precision is achieved for heavier damping. This is because lightly
Fig. 1. Nine DOF structure. The node numbers are given by the subscripts on the displacements shown.
380
J.H. Lin et al. / Computers and Structures 79 (2001) 375±387
Fig. 2. Loading PSD functions of points 3 and 6 identi®ed from the displacement PSD functions of points 2, 5 and 7, with all nine modes participating in the back analysis. The two numbers beneath each ®gure are the values of all damping ratios and the number of digits accuracy of the displacement PSD values.
damped structures can experience serious ill conditioning in the vicinity of their natural frequencies. This problem was overcome when the displacement PSDs were assumed to have an accuracy of three eective digits, but unfortunately in normal circumstances it would be virtually impossible to measure PSDs to such high accuracy. Hence, it is practically impossible to avoid the abnormal ¯uctuation (or ``jumping'') shown in
Fig. 2(a) and (c) under normal experimental conditions, as is now discussed further. 7.1.3. Further analysis on the abnormal jumping of identi®cation results in the vicinity of natural frequencies For the present example H is a diagonal matrix with given damping ratios, so that Eq. (26) can be converted into
J.H. Lin et al. / Computers and Structures 79 (2001) 375±387
R Ey UH UT Ef q X Hj xEy f/gj f/gTj Ef ;
50
j1
where
ÿ1 Hj x x2j ÿ x2 2i1j xj x :
51
It is known that R will be singular if its rank is less than q [13±15]. Hence, if q0 of the q terms summed in Eq. (50) are zero, R must be singular and of rank q ÿ q0 . Furthermore, Eq. (51) shows that Hj x 21j xj xÿ1 when x xj , so that if the jth damping ratio 1j is very small, Hj x will be very large, particularly when xj is small. Hence one (or more) of the terms summed in Eq. (50) becomes almost negligible, which causes the pathological jumping of the identi®ed curves, i.e. it causes the high sensitivity of identi®cation results to the measured values. This phenomenon can be elucidated further by examining intermediate computation data. For example, when x 1:00, abnormal jumping occurs in the curves of Fig. 2(c). The displacement response PSD matrix of points 2, 5 and 7, which was used to identify the excitation PSD matrix, was (taking only the ®rst two eective digits) 2 3 1:2 1:1 1:1 6 7 Syy 4 1:1 1:1 1:0 5 102 1:1 1:0 0:99 3 2 0:0 2:8 ÿ3:2 7 6 52 i4 ÿ2:8 0:0 ÿ5:7 5: 3:2
5:7
0:0
The identi®ed excitation PSD values for nodes 3 and 6 were fSf3 f3 ; Sf6 f6 ; Sr;f3 f6 ; Si;f3 f6 g f2:187; 3:210; 2:795; 0:932g;
53
where the ®rst four digits of the IPEM results are given for possible reference, although only two digits would normally be trusted. The exact theoretical values (Eq. (49)), were {3.000, 4.000, 2.000, 1.000}. Now, suppose that a measurement error in the third digit results in the value 1.2 of the (1,1) element of Syy being replaced by 1.25. Then, the IPEM back analysis would give fSf3 f3 ; Sf6 f6 ; Sr;f3 f6 ; Si;f3 f6 g f2:911; 3:860; 2:135; 1:120g:
54
Hence, an error of 4.2% in only one of the measured values results in dierences in the loading identi®cation of {24.9%, 16.8%, ÿ30.9%, 11.4%}. This demonstrates the high sensitivity of the loading identi®cation to the response measurement accuracy.
381
Because this problem is highly pathological, the corresponding ill-conditioned Eq. (27) were solved by four methods to ensure that the results are correct. These were the singular value decomposition (SVD) and Cholesky triangularization decomposition methods [14,15], plus two other methods for solving pathological equations [11]. All four methods gave practically identical results, i.e. at least the ®rst ®ve eective digits agreed on an IBM/586 PC computer. Therefore, it is certain that Eq. (27) was solved correctly and hence that the identi®cation problem was due to the poor measurement precision and not due to inaccurate solution of pathological equations. Since it is very dicult to measure to a precision of three reliable digits, jumping is inevitable in the back analysis. However, this does not mean that such abnormal results must be accepted because, if such jumping happens only near the natural frequencies, it can be smoothed out by using adjacent points on the curves of, for example (Fig. 2(c)). 7.1.4. Selection of measured points by computer simulation Improper selection of measurement points is another possible cause of unfavourable results, as follows. When the measurement points of the above example were changed from 2, 5 and 7 to 2, 4 and 7, with only the ®rst two digits of the displacement response PSD values used for loading identi®cation, Fig. 2(a) was replaced by Fig. 3(a), in which considerable jumping is observed near the natural frequency x 1:54. Even when this point was modi®ed by taking the average of two measured values on either side of it, the remainder of the identi®cation curves are still inadequate, Fig. 3(b). Only when using three or more eective digits, which means that the noise disturbance is extremely weak, is the identi®cation from this combination good enough. Clearly, such unfavourable combinations of measurement points should be avoided because they cannot resist noise-disturbance adequately. The problem is that most engineers would ®nd it very dicult to judge which is the best combination of measurement points to use, i.e. to choose which is the better of 2, 4 and 7 or 2, 5 and 7. Therefore, selection of measurement points by means of computer simulation, as above, must be a good suggestion, since, in current engineering practice, such selection is usually based on engineersÕ experience or intuitive feelings, which may be an important reason why relatively little successful identi®cation of structural random loading has so far been reported. It may be noted that although the above suggestions have provided a practical way for engineers to select good measurement points which is both easy to understand and simple to implement, a stricter and more extensive investigation on the error caused by the noise disturbance still needs to be further developed.
382
J.H. Lin et al. / Computers and Structures 79 (2001) 375±387
Fig. 3. Loading PSD functions of points 3 and 6 identi®ed from the displacement PSD functions of points 2, 4 and 7, with all nine modes participating in the back analysis. All damping ratios were 0.03 and two digits of accuracy were retained for the displacement PSD values. The two cases are for when the maximum jumping point was unaltered or was revised.
7.1.5. In¯uence of the accuracy of structural parameters on the identi®cation precision Inaccurate evaluation of structural stiness and mass parameters can have a remarkable eect on loading identi®cation precision. For example, changing the stiness of the middle (i.e. ®fth) element of the structure represented in Fig. 1 from 2k to 2.04k, altered the structural natural frequencies to 1.00122, 1.03736, 1.13924, 1.27174, 1.41598, 1.54359, 1.64685, 1.70994, 1.73279 and replaced the identi®cation results of Fig. 2(a) by those of Fig. 4(a), for which each curve deviates from the theoretical one more considerably, with a maximum of 12.6% for Sr;f3 f6 2:252 at x 0:98. As a second example, when instead the middle (i.e. ®fth) mass of the structure of Fig. 1 was changed from 2m to 2.04m, the natural angular frequencies became 0.99937, 1.03736, 1.13563, 1.27174, 1.41245, 1.54359, 1.64329, 1.70994, 1.73101 and the identi®cation results of Fig. 2(a) were replaced by those of Fig. 4(b), for which the maximum deviation of any curve from the theoretical value is 9.7% for Sf6 f6 4:386 at x 1:16. As a third example, the changes of the ®rst two examples were made simultaneously. Then, Fig. 2(a) was replaced by Fig. 4(c), for which every curve deviated from its theoretical value by much smaller amounts than for the ®rst two examples; the maximum error being 6.4% for Sf3 f3 2:80 at x 0:98. The amounts are smaller because the natural angular frequencies were 1.00062, 1.03736, 1.13783, 1.27174, 1.41421, 1.54359, 1.64479, 1.70994 and 1.73170 which are, as a whole, closer to the
angular frequencies of the original structure than were those of the ®rst two examples. In fact, the second, fourth, sixth and eighth natural frequencies are identical, while the dierences of the other ®ve frequencies are small. The results of Fig. 4(a)±(c) show that even when the structural stiness and mass parameters are not evaluated very precisely, so long as the frequency spectrum of the structure can be measured fairly precisely, the identi®cation can be acceptably accurate. 7.1.6. The in¯uence of the number of participant modes on the identi®cation precision In all the above discussions, all nine modes of the structure were used in the direct analyses, i.e. when computing the response PSD from the excitation PSD, as well as in the back analyses, i.e. when computing the excitation PSD from the response PSD, by the modesuperposition method. Now suppose that all nine modes are still used for the direct analyses needed to obtain the simulated ``measured values'', but that less modes are used for the back analyses. For instance, taking q 7 (modes) and using the ®rst two digits in all the measured values gave the identi®ed excitation PSD curves of Fig. 5(a), which deviate considerably from the curves of Fig. 2(a), for which q 9, particularly in the high frequency region x > 1:4 for which the dierence can be hundreds of percent. Even within the region x 2 [0.8, 1.4], the maximum deviation of the identi®ed values from the theoretical ones is about 30%, Fig. 5(b). It should be
J.H. Lin et al. / Computers and Structures 79 (2001) 375±387
383
Fig. 4. Loading PSD functions of points 3 and 6 identi®ed from the displacement PSD functions of points 2, 5 and 7, with all nine modes participating in the back analysis. All damping ratios were 0.03 and two digits of accuracy were retained for the displacement PSD values.
noted that all the natural frequencies of this structure are closely spaced, the ratio of the highest to the lowest being only 1.732. This is why neglecting only the two highest modes caused considerable in¯uence over the whole frequency region. However, computer simulations for structures with widely spaced natural frequency spectra show that neglecting higher modes aects the
identi®cation precision relatively little in the lower frequency region. 7.2. Example 2: jacket platform Fig. 6 shows the FEM model of a jacket platform which has been built in the Bohai sea of China and has
384
J.H. Lin et al. / Computers and Structures 79 (2001) 375±387
Fig. 5. Loading PSD functions of points 3 and 6 identi®ed from the displacement PSD functions of points 2, 5 and 7, with the ®rst seven modes participating in the back analysis. All damping ratios were 0.03 and two digits of accuracy were retained for the displacement PSD values. The two alternative frequency regions shown were used.
Fig. 6. Jacket platform: only those nodes and members referred to in the text are, respectively, numbered and lettered.
183 beam elements and 492 DOF. Its ®rst 16 natural angular frequencies are 3.8325, 3.8811, 4.1619, 7.9105, 9.9539, 11.205, 11.547, 12.091, 16.182, 17.518, 17.898,
19.161, 21.152, 21.537, 22.777 and 23.160 (sÿ1 ) and the damping ratio of all modes is 0.03. Flowing ice blocks are assumed to cause horizontal x direction stationary random excitations at nodes 5 and 11. The excitation PSD matrix for computer simulation is still taken as Eq. (49). A direct PEM analysis [6] using the ®rst 16 modes was performed to obtain the response PSD matrix for the z direction displacements at points 6, 10, 17 and 23. Its ®rst two eective digits were taken as the ``measured values'' for computer simulation. The loading identi®cation results obtained from the proposed IPEM are shown in Fig. 7(a) and are clearly very poor, i.e. far from the true excitation PSD values of Eq. (49). Therefore, the z direction displacement responses were replaced to obtain Fig. 7(b)±(d), which were obtained from the z displacement responses at the points given in their captions. Clearly, Fig. 7(c) gives the best identi®cation results and they are of usable accuracy. Hence using the z direction displacements at the points 7, 9, 50 and 86 would be a good choice of response to use, if real engineering measurements are made on the actual structure. Suppose that instead the axial strain responses of members A, B and C of Fig. 6 (i.e. beams 9±21, 11±23 and 17±25) are known and that only their ®rst two effective digits are used for loading identi®cation. The results are shown in Fig. 8(a) and are clearly quite good. However, Fig. 8(b) shows that even better results were obtained by including the vertical displacement response
J.H. Lin et al. / Computers and Structures 79 (2001) 375±387
385
Fig. 7. Loading PSD functions of the x-direction forces at nodes 5 and 11, obtained by using two digit accuracy for the displacement PSD values, when the ®rst seven modes participated in the back analysis. The results were identi®ed from the PSD functions of the vertical displacements at the four points shown in the sub-captions.
curve at node 21 in addition to the strains of members A, B and C. In contrast, when the known responses were chosen to be the axial strains of members D and E (i.e. beams 5±17 and 7±19) and the vertical displacements of nodes 50 and 86, the identi®cation results were much poorer, Fig. 8(c), and so this combination should not be
selected for practical measurements on the real structure. However, the results of Fig. 8(d) show that it would only be necessary to replace the vertical displacement response at node 86 by the horizontal displacement response at node 21, to again obtain very good results.
386
J.H. Lin et al. / Computers and Structures 79 (2001) 375±387
Fig. 8. Loading PSD functions of the x-direction forces at nodes 5 and 11, obtained by using two digit accuracy for the displacement PSD values, when the ®rst seven modes participated in the back analysis. The results were identi®ed from the PSD functions of the quantities shown in the sub-captions.
8. Conclusions A random loading identi®cation approach based on a proposed IPEM has been presented and a primary investigation has been given. Previously, such problems have not been solved well. This paper reveals some im-
portant phenomena which may have contributed towards this. For complex engineering structures, the selection of the responses to be used for loading identi®cation is an important problem which has not yet been suciently resolved and for which experience-based methods are
J.H. Lin et al. / Computers and Structures 79 (2001) 375±387
not very reliable. The examples presented show that computer simulation combined with the IPEM is a good tool for selecting measured responses and hence, for identifying the random loadings. This is very important because a real engineering measurement project may take months or even years to prepare and is very costly, whereas computer simulation is comparatively very quick and cheap. For the second example, i.e. the platform, the direct analysis used about 10 min on an IBM/ 486 personal computer and only needs to be done once, while the back analysis took about 15 min for each of the above response combinations, i.e. the computations were very cheap. Therefore, performing sucient computer simulations before practical measurements are undertaken is strongly advisable. Some basic phenomena have been observed and analyzed in the present paper and some suggestions have been given which should be helpful in engineering applications. However, more work remains to be done, e.g. the problem of noise disturbance should be extensively investigated in a stricter way and the bene®ts could be investigated of using Kalman or H1 ®ltering techniques to further improve the identi®cation precision. Acknowledgements The authors are grateful for the support of the National Natural Science Foundation of China (project no. 19772009), the funding of NKBRSF (no. G1999032805) and the State Education Commission of China (no. 97014120), as well as the Advanced Chinese Engineering (ACE) Centre at Cardi University. The authors also wish to thank Bohai Petroleum Company for providing practical data of the platform of example 2 and for their ®nancial support for the research.
387
References [1] Nigam NC. Introduction to random vibration. Cambridge, Massachussets: MIT Press; 1983. [2] Bendat JS, Piersol AG. Engineering applications of correlation and spectral analysis. New York: Wiley; 1980. [3] Zhang JH, Wang C. Theory of engineering random vibration. Xian, China: Xian Jiao Tong University; 1988. [4] Lin JH. A fast CQC algorithm of PSD matrices for random seismic responses. Comput Struct 1992;44:683±7. [5] Zhong WX. Review of a high-eciency algorithm series for structural random responses. Prog Nat Sci 1996;6:257±68. [6] Lin JH, Zhang WS, Li JJ. Structural responses to arbitrarily coherent stationary random excitations. Comput Struct 1994;50:629±33. [7] Ewins DJ. Model testing: theory and practice. London: Research Studies Press; 1980. [8] Farrar CR, James GH. System identi®cation from ambient vibration measurements on a bridge. J Sound Vibr 1997; 205:1±18. [9] Doebling SW, Farrar CR, Goodman RS. Eects of measurement statistics on the detection of damage in the Alamosa Canyon Bridge. Proc Fifteenth Int Modal Anal Conf-IMAC. Orlando, FL, USA, 1997. p. 919±29. [10] Askari MR. System identi®cation using vibration monitoring for fault diagnostics in manufacturing machinery. Proc Fifth Int Conf Factory 2000 ± the technology exploitation process. Cambridge, UK, 1997. p. 295±302. [11] Golub GH, Van Loan CF. Matrix computations. Baltimore: Johns Hopkins University Press; 1983. [12] Cook RD. Concept and application of ®nite element analysis. New York: McGraw-Hill; 1981. [13] Bathe KJ, Wilson EL. Numerical methods in ®nite element analysis. Englewood Clis, NJ: Prentice-Hall; 1976. [14] Wilkinson JH, Reinsch C. Linear algebra, vol. II. Berlin: Springer, 1971. [15] Zhong WX, He Q, Liu ZX. Numerical computational methods. Beijing: Press of Architectural Industry of China; 1991.
www.elsevier.com/locate/compstruc
Computer simulation of structural random loading identi®cation J.H. Lin a,1, X.L. Guo a, H. Zhi a, W.P. Howson b,*, F.W. Williams b a
State Key Laboratory of Structural Analysis of Industrial Equipment, Dalian University of Technology, Dalian 116023, People's Republic of China b Division of Structural Engineering, Cardi School of Engineering, Cardi University, Cardi CF24 3TB, UK Received 13 August 1999; accepted 27 May 2000
Abstract Loading identi®cation has long been a dicult problem for multi-input±multi-output stationary random vibration problems. In this paper, the inverse pseudo-excitation method is developed for dealing with such identi®cation problems. Numerical simulations are used to reveal the diculties in such problems, and some numerical analysis measurements are suggested, which may be very useful in the setting up and processing of experimental data so as to obtain reasonable predictions of the input loading from selected structural responses, e.g. a few structural displacements or strains, or a combination of them. It is shown that computer simulation is a very cheap way of overcoming the problem of which responses to measure on real structures in order to identify the power spectral density matrix of the applied forces. Ó 2000 Elsevier Science Ltd. All rights reserved. Keywords: Random vibration; Loading identi®cation; Computer simulation
1. Introduction For linear structures subjected to multiple stationary random excitations, which may be partially coherent, the direct problem, i.e. the computation of the power spectral densities (PSDs) of the various responses from the given excitation PSD, has been solved [1±3], in particular by using the exact and ecient pseudoexcitation method (PEM) [4±6]. There are two types of inverse problems (also called back analysis problems) corresponding to this direct problem. The ®rst is the socalled structure identi®cation problem, for which the loading and response PSD functions are all known and are used to identify the properties of the structure. Many publications cover such problems and their successful
* Corresponding author. Address: Division of Structural Engineering, Cardi School of Engineering, Cardi University, Cardi CF24 3TB, UK. 1 Visiting Professor at the University of Wales Cardi.
use in engineering [2,3,7±10]. The second is the so-called loading identi®cation problem, for which the response information and structural properties are known and are used to identify information about the loading. For this second type of inverse problem, if there is only one excitation acting on the structure, the problem is quite straightforward and is not discussed in this paper. For multi-input±multi-output (MIMO) problems, however, its solution is not easy and it is dicult to ®nd papers in this area which demonstrate its diculty. To the authorsÕ knowledge, some such work has been done at high expense, but has nevertheless yielded identi®cation results that are unacceptable because of their very poor precision. However, such inverse problems are of great importance in practical engineering because sometimes the direct measurement of the applied loads may be very dicult, whereas the measurement of some of the responses may be relatively easy. For any direct problem for which the excitation PSD matrix has been given, the PSD matrix of any responses can be computed quickly and precisely by means of the
0045-7949/01/$ - see front matter Ó 2000 Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 4 9 ( 0 0 ) 0 0 1 5 4 - 1
376
J.H. Lin et al. / Computers and Structures 79 (2001) 375±387
PEM [4±6]. In this paper, the PEM is modi®ed to generate an inverse pseudo-excitation method (IPEM), so that the excitation PSD matrix can be obtained from a set of known response PSD functions which can be displacements, strains, or their combination. Two examples of computer simulation are given, namely a plane truss and a three dimensional oshore jacket platform built in the Bohai sea in Northern China. These computer simulation examples show that the extended IPEM has good noise-resistance characteristics and so may be quite useful for practical engineering problems involving large FEM models.
When all excitations are fully coherent, i.e. r 1, Eq. (7) takes the following simple form:
2. Pseudo-excitation method for solving to ®nd the stationary random responses of structures
Sxx H Syy H T
The conventional formula for solving the stationary random responses of linear structures in the frequency domain is Syy H Sxx H T
1
in which Sxx is the known excitation PSD matrix, H is the transfer function matrix, Syy is the response PSD matrix to be computed and the superscripts and T represent complex conjugate and transpose, respectively. The excitation PSD matrix Sxx can generally be assumed to be a p p matrix with rank r ( 6 p). Since Sxx must be a Hermitian matrix [3,4], it can be decomposed into Sxx
r X fagj fagTj :
2
j1
Syy H Sxx H T H fag fagT HT fbg fbgT :
8
3. The inverse pseudo-excitation method for structural random loading identi®cation Let Eq. (1) be written in the form [2] 9
in which the superscript represents generalized inversion (performed as described beneath Eq. (29)), since H is in general not a square matrix. Provided that the responses of the structure are caused entirely by the random loadings to be identi®ed, then Eq. (9) gives the algorithm for computing Sxx from the known PSD matrix Syy . Clearly, Eq. (9) is the backward extension of the PEM described above. It is therefore, called the IPEM. For complicated engineering problems, direct use of Eq. (9) is quite inecient, and so it should be dealt with analogously to the last section, as follows: Pr Using Eq.T (7), Syy can be decomposed into ixt j1 fbgj fbgj . The pseudo-response f y gj fbgj e then follows from Eq. (5), and applying Eq. (6) inversely, gives the pseudo-excitations f xgj of Eq. (3) as f xgj fagj eixt H fbgj eixt :
10
Introducing the pseudo-excitation fxgj fagj eixt
j 1; 2; . . . ; r
3
Hence Eqs. (9), (7) and (10) enable the excitation PSD matrix Sxx to be given in terms of the fagj as
4
Sxx H Syy H T H
leads to the harmonic response fyg given by [5] fygj H fxgj ; so that, using Eq. (3), fygj fbgj eixt ;
5
where fbgj H fagj :
6
Hence, using Eqs. (1), (2) and (4), the response PSD matrix Syy is given in terms of the fbgj as ! r X T T Syy H Sxx H H fagj fagj H T j1
r X j1
fbgj fbgTj :
7
r X fbgj fbgTj H T j1
r X fagj fagTj :
11
j1
When all excitations are fully coherent, i.e. r 1, Eq. (11) takes the following simple form: Sxx H Syy H T H fbg fbgT H T fag fagT :
12
For large systems, the direct use of Eq. (12) or Eq. (11) is still computationally very expensive. Therefore, the following equation reduction scheme should be used, which is based on the mode-superposition method.
J.H. Lin et al. / Computers and Structures 79 (2001) 375±387
4. Excitation power spectral density identi®cation for large systems The equations of motion of a large system can be expressed as MfY g CfY_ g KfY g fF g
13
in which M, C and K are, respectively, the mass, damping and stiness matrices, dots denote dierentiation with respect to time and fY g and fF g are, respectively, the displacement and applied force vectors. Even though the order n of fF g and fY g is very high, the number, l, of non-zero excitations in fF g and the number, m, of measured responses in fY g are both n. Therefore, let fF g Ef ff g;
14
fyg Ey fY g;
15
where ff g is the l dimensional vector consisting of the non-zero elements in fF g, fyg is the m dimensional vector consisting of the m measured quantities in fY g and Ef and Ey are matrices which contain only zero and unity. Assume that the m m matrix Syy x is known, usually being speci®ed in the form of a series of values at discrete frequency points, so that the requirement is to ®nd the l l excitation PSD matrix Sff x by back analysis. The mode-superposition scheme is now used. Firstly, the lowest q modes, denoted by the n q matrix U, and the corresponding diagonal q q eigenvalue matrix X2 must be found by solving KU MUX2 ;
16
subject to UT M U I ; where I is a q q unit matrix:
17
Reducing Eq. (13) by means of the modal matrix U gives fY g UfUg;
18
g C0 fU_ g X2 fU g UT Ef ff g fU
19
in which C0 UT CU:
20
For proportionally damped systems, C 0 is a diagonal matrix. In particular, if its damping ratios 1j j 1; 2; . . . ; q are all known, then the jth diagonal element of C 0 will be 21j xj , where x2j is the jth diagonal element of X2 . For non-proportional damping matrices, C 0 is non-diagonal. However, for both cases, the following derivation applies, in which partial coherency is assumed for generality.
377
Substituting Eq. (18) into Eq. (15) gives fyg Ey UfU g:
21
Then, decomposing Syy x by using Eq. (7), and forming the pseudo-responses according to Eq. (5), gives f y gj fbgj eixt , which from Eq. (21) can be regarded as being produced by the pseudo-ÔexcitationÕ, fU gj fcgj eixt :
22
This fU gj can be considered to be generated by the pseudo-excitation f f gj fagj eixt , so that Eq. (19) gives fU gj H UT Ef ff gj
23
in which H ÿx2 I X2 ixC0 ÿ1 ; 24 p where i ÿ1. Substituting Eq. (23) into Eq. (21) produces fygj Rff gj
25
in which R Ey UH UT Ef :
26
Therefore, ff gj R fygj :
27
Finally, using the IPEM gives (Eq. (11)) X Sff R Syy RT f f gj f f gTj :
28
When the number of excitations l exceeds the number of measured responses m, the m l generalized inverse matrix of R, i.e. R RT R ÿ1 RT
29
does not exist, and so Sff cannot be found from Syy . When l < m, the inversion procedure used is to solve for the l unknowns from the m equations by means of the least square method, so that only approximate solutions can be expected [11]. Clearly, when l m, the equation set can be solved exactly. Therefore, in general, the inverse problem is solvable under the conditions l 6 m and l 6 q. However, even when these conditions are met, there is no guarantee that a satisfactory excitation PSD matrix will be obtained. Practical situations are rather complex, as discussed via the following simulation example. For many engineering problems, strain measurements are more convenient and/or more precise than are displacement measurements. Assume that the vector {e} of measured strains is related to a reduced set of displacements fyg of the structural displacement vector fY g by
378
J.H. Lin et al. / Computers and Structures 79 (2001) 375±387
feg Ee fyg
30
in which the transformation matrix Ee can be calculated by standard FEM [2±4] and fyg can be computed in terms of Eq. (15). Substituting Eqs. (15) and (18) into Eq. (30) gives feg Ee Ey UfUg:
31
Assuming that p suciently long time-history curves of structural strains have been recorded, the corresponding p p PSD matrices See can easily be obtained at discrete frequency points by means of a power spectrum analyzer. Because p is small, See can be easily decomposed into the form of Eq. (2), i.e. See
r X fggj fggTj :
32
j1
where r is the rank of the p p PSD matrix See (r 6 p). The corresponding pseudo-strain responses are then fegj fggj eixt
33
and can be regarded as being caused by pseudodisplacement excitations of the form of Eq. (22). Hence, Eq. (31) gives the corresponding strain vector fegj Ee Ey UfU gj ;
34
where fU gj can be regarded as being caused by pseudoexcitation of the form ff gj fagj eixt : From Eq. (19), fU gj can be expressed in terms of f f gj by Eq. (23), which is then further substituted into Eq. (34) to produce fegj Re ff gj
35
In practical engineering, good measured data is often quite limited. Hence, sometimes the loading cannot be identi®ed adequately or at all from either the displacement or strain responses alone. Then, combined use of the displacement and strain data may make identi®cation possible, or give better identi®cation. It will be seen that the derivation and implementation of such a combined method, as described below, will have similarities to those in the last section. Consider the displacement vector fyg which contains exclusively all the measured displacements and all of the displacements associated with the measured strains. Therefore, the measured strains feg and the measured displacements fyd g can be expressed as feg Ke fyg;
39
fyd g Ed fyg
40
in which Ke is a local stiness matrix, Ed contains only ones and zeros, and both Ke and Ed can be calculated by means of standard FEM [12,13]. These two equations can be combined to give fym g Em fyg in which e fym g ; yd
41 Em
Ke : Ed
42
The use of all available ( s) measured curves of feg and fyd g enables the corresponding s s PSD matrix Sym ym to be found in discrete form by means of a power spectrum analyzer. Then decomposition gives Sym ym
in which
r X fym gj fym gTj ;
43
j1
Re Ee Ey UH UT Ef :
36
Hence, ff gj Re fegj :
37
r X ff gj ff gTj ;
where r ( 6 s) is the rank of the matrix Sym ym . In a similar manner to the derivation in the previous section, fym gj can be regarded as being caused by the pseudodisplacement fU gj fcgj eixt . Hence (Eqs. (15), (18) and (41)), fym gj Em Ey UfU gj :
Finally, f f gj is used to calculate Sff from Sff Re See Re T
5. Identi®cation of the excitation PSD matrix from the mixed displacement±strain PSD matrix
38
j1
where Re is in general not a square matrix. Hence, the inversion procedure used, and a discussion of when such inversion is possible, is the one given above, beneath Eq. (29). However, the conditions for inversion of Re to be possible are now slightly stricter than for the R of Eq. (26) because the additional condition l 6 p is required.
44
Substituting Eq. (23) into Eq. (44) gives fym gj Rm ff gj
45
in which Rm Em Ey UH UT Ef :
46
Therefore, ff gj Rm fym gj ;
47
J.H. Lin et al. / Computers and Structures 79 (2001) 375±387 r X Sff Rm Sym ym Rm T ff gj ff gTj
48
j1
in which the possibility of inverting to obtain Rm is subject to the additional condition that r P l, when compared to inverting to obtain R of Eq. (27), where l is the number of excitations to be identi®ed. Finally, the IPEM, i.e. the last expression of Eq. (38), can be used to give the excitation PSD matrix Sff .
6. Implementation procedure of the IPEM The steps of the implementation procedure are as follows: (1) Calculate the matrix R from Eq. (26), Re from Eq. (36), or Rm from Eq. (46) with H optionally formed by experimental means. (2) Calculate the response PSD matrix Syy , See or Sym ym for practical engineering problems by measurement and data processing. However, for the purpose of computer simulations of this paper, such matrices were obtained by direct analyses and then represented to very few digits accuracy (in general two) when performing the back analyses. (3) Decompose the response PSD matrix by using the last expression of the appropriate one of Eqs. (7), (32) or Eq. (43) to generate the corresponding pseudoresponses. (4) Solve Eqs. (27), (37) or Eq. (47) for all excitations f f gj [14], and then compute the excitation PSD matrix Sff from Eqs. (28), (38) or Eq. (48).
7. Computer simulation example 7.1. Example 1: nine DOF structure The stiness and mass distributions of a nine DOF structure are shown in Fig. 1. For simplicity, dimensionless quantities m 1 and k 1 are assumed. The complete set of natural angular frequencies are 1.00000, 1.03736, 1.13705, 1.27174, 1.41421, 1.54359, 1.64533, 1.70994, 1.73205 and the damping ratios are all equal
379
with a value of 0.03. Two horizontal forces are applied at nodes 3 and 6, and the corresponding PSD matrix is Sff
3:0 2:0
2:0 0:0 1:0 i 4:0 ÿ1:0 0:0
49
within the signi®cant frequency region x 2 0:8; 1:8: This region is discretized with an interval of Dx 0:01. Because of the existence of noise disturbance, practical responses cannot be measured very precisely. Because of this, only two or three digits of the response PSD results generated from the direct random vibration analysis [4] were taken as the ``measured values'' for the computer simulation. Some observations of interest follow in the following six sub-sections. All results presented were obtained by the IPEM and are curves of the auto-PSD of the excitations f3 and f6 , i.e. Sf3 f3 and Sf6 f6 , and their cross-PSD Sf3 f6 , which has real and imaginary parts Sr;f3 f6 and Si;f3 f6 , respectively. 7.1.1. The eect of the precision of measurement of the response PSD matrix on the precision with which the excitation PSD matrix is identi®ed Because noise always exists in the measurements, only a limited number of digits of the measured response values are reliable for use in the loading identi®cation. When only the ®rst two eective digits of all measured values were used, the results obtained for all damping ratios (0.03) are plotted in Fig. 2(a). It can be seen that the identi®cation values of the PSD functions ¯uctuate around their exact values, Eq. (49), of 3.0, 4.0, 2.0 and 1.0. The maximum error (taken directly from the computed results) is 6.5% for Sr;f3 f6 1:869 when x 1:00. In contrast, if the measurement precision is so high that the ®rst three eective digits of all the measured quantities are reliable, Fig. 2(b) replaces Fig. 2(a), and obviously, the four curves obtained are very close to their theoretical (i.e. exact) values. 7.1.2. The eect of structural damping on the identi®cation precision When the damping ratios were all changed from 0.03 to 0.01, Fig. 2(a) and (b) becomes Fig. 2(c) and (d). Hence, it can be seen that higher identi®cation precision is achieved for heavier damping. This is because lightly
Fig. 1. Nine DOF structure. The node numbers are given by the subscripts on the displacements shown.
380
J.H. Lin et al. / Computers and Structures 79 (2001) 375±387
Fig. 2. Loading PSD functions of points 3 and 6 identi®ed from the displacement PSD functions of points 2, 5 and 7, with all nine modes participating in the back analysis. The two numbers beneath each ®gure are the values of all damping ratios and the number of digits accuracy of the displacement PSD values.
damped structures can experience serious ill conditioning in the vicinity of their natural frequencies. This problem was overcome when the displacement PSDs were assumed to have an accuracy of three eective digits, but unfortunately in normal circumstances it would be virtually impossible to measure PSDs to such high accuracy. Hence, it is practically impossible to avoid the abnormal ¯uctuation (or ``jumping'') shown in
Fig. 2(a) and (c) under normal experimental conditions, as is now discussed further. 7.1.3. Further analysis on the abnormal jumping of identi®cation results in the vicinity of natural frequencies For the present example H is a diagonal matrix with given damping ratios, so that Eq. (26) can be converted into
J.H. Lin et al. / Computers and Structures 79 (2001) 375±387
R Ey UH UT Ef q X Hj xEy f/gj f/gTj Ef ;
50
j1
where
ÿ1 Hj x x2j ÿ x2 2i1j xj x :
51
It is known that R will be singular if its rank is less than q [13±15]. Hence, if q0 of the q terms summed in Eq. (50) are zero, R must be singular and of rank q ÿ q0 . Furthermore, Eq. (51) shows that Hj x 21j xj xÿ1 when x xj , so that if the jth damping ratio 1j is very small, Hj x will be very large, particularly when xj is small. Hence one (or more) of the terms summed in Eq. (50) becomes almost negligible, which causes the pathological jumping of the identi®ed curves, i.e. it causes the high sensitivity of identi®cation results to the measured values. This phenomenon can be elucidated further by examining intermediate computation data. For example, when x 1:00, abnormal jumping occurs in the curves of Fig. 2(c). The displacement response PSD matrix of points 2, 5 and 7, which was used to identify the excitation PSD matrix, was (taking only the ®rst two eective digits) 2 3 1:2 1:1 1:1 6 7 Syy 4 1:1 1:1 1:0 5 102 1:1 1:0 0:99 3 2 0:0 2:8 ÿ3:2 7 6 52 i4 ÿ2:8 0:0 ÿ5:7 5: 3:2
5:7
0:0
The identi®ed excitation PSD values for nodes 3 and 6 were fSf3 f3 ; Sf6 f6 ; Sr;f3 f6 ; Si;f3 f6 g f2:187; 3:210; 2:795; 0:932g;
53
where the ®rst four digits of the IPEM results are given for possible reference, although only two digits would normally be trusted. The exact theoretical values (Eq. (49)), were {3.000, 4.000, 2.000, 1.000}. Now, suppose that a measurement error in the third digit results in the value 1.2 of the (1,1) element of Syy being replaced by 1.25. Then, the IPEM back analysis would give fSf3 f3 ; Sf6 f6 ; Sr;f3 f6 ; Si;f3 f6 g f2:911; 3:860; 2:135; 1:120g:
54
Hence, an error of 4.2% in only one of the measured values results in dierences in the loading identi®cation of {24.9%, 16.8%, ÿ30.9%, 11.4%}. This demonstrates the high sensitivity of the loading identi®cation to the response measurement accuracy.
381
Because this problem is highly pathological, the corresponding ill-conditioned Eq. (27) were solved by four methods to ensure that the results are correct. These were the singular value decomposition (SVD) and Cholesky triangularization decomposition methods [14,15], plus two other methods for solving pathological equations [11]. All four methods gave practically identical results, i.e. at least the ®rst ®ve eective digits agreed on an IBM/586 PC computer. Therefore, it is certain that Eq. (27) was solved correctly and hence that the identi®cation problem was due to the poor measurement precision and not due to inaccurate solution of pathological equations. Since it is very dicult to measure to a precision of three reliable digits, jumping is inevitable in the back analysis. However, this does not mean that such abnormal results must be accepted because, if such jumping happens only near the natural frequencies, it can be smoothed out by using adjacent points on the curves of, for example (Fig. 2(c)). 7.1.4. Selection of measured points by computer simulation Improper selection of measurement points is another possible cause of unfavourable results, as follows. When the measurement points of the above example were changed from 2, 5 and 7 to 2, 4 and 7, with only the ®rst two digits of the displacement response PSD values used for loading identi®cation, Fig. 2(a) was replaced by Fig. 3(a), in which considerable jumping is observed near the natural frequency x 1:54. Even when this point was modi®ed by taking the average of two measured values on either side of it, the remainder of the identi®cation curves are still inadequate, Fig. 3(b). Only when using three or more eective digits, which means that the noise disturbance is extremely weak, is the identi®cation from this combination good enough. Clearly, such unfavourable combinations of measurement points should be avoided because they cannot resist noise-disturbance adequately. The problem is that most engineers would ®nd it very dicult to judge which is the best combination of measurement points to use, i.e. to choose which is the better of 2, 4 and 7 or 2, 5 and 7. Therefore, selection of measurement points by means of computer simulation, as above, must be a good suggestion, since, in current engineering practice, such selection is usually based on engineersÕ experience or intuitive feelings, which may be an important reason why relatively little successful identi®cation of structural random loading has so far been reported. It may be noted that although the above suggestions have provided a practical way for engineers to select good measurement points which is both easy to understand and simple to implement, a stricter and more extensive investigation on the error caused by the noise disturbance still needs to be further developed.
382
J.H. Lin et al. / Computers and Structures 79 (2001) 375±387
Fig. 3. Loading PSD functions of points 3 and 6 identi®ed from the displacement PSD functions of points 2, 4 and 7, with all nine modes participating in the back analysis. All damping ratios were 0.03 and two digits of accuracy were retained for the displacement PSD values. The two cases are for when the maximum jumping point was unaltered or was revised.
7.1.5. In¯uence of the accuracy of structural parameters on the identi®cation precision Inaccurate evaluation of structural stiness and mass parameters can have a remarkable eect on loading identi®cation precision. For example, changing the stiness of the middle (i.e. ®fth) element of the structure represented in Fig. 1 from 2k to 2.04k, altered the structural natural frequencies to 1.00122, 1.03736, 1.13924, 1.27174, 1.41598, 1.54359, 1.64685, 1.70994, 1.73279 and replaced the identi®cation results of Fig. 2(a) by those of Fig. 4(a), for which each curve deviates from the theoretical one more considerably, with a maximum of 12.6% for Sr;f3 f6 2:252 at x 0:98. As a second example, when instead the middle (i.e. ®fth) mass of the structure of Fig. 1 was changed from 2m to 2.04m, the natural angular frequencies became 0.99937, 1.03736, 1.13563, 1.27174, 1.41245, 1.54359, 1.64329, 1.70994, 1.73101 and the identi®cation results of Fig. 2(a) were replaced by those of Fig. 4(b), for which the maximum deviation of any curve from the theoretical value is 9.7% for Sf6 f6 4:386 at x 1:16. As a third example, the changes of the ®rst two examples were made simultaneously. Then, Fig. 2(a) was replaced by Fig. 4(c), for which every curve deviated from its theoretical value by much smaller amounts than for the ®rst two examples; the maximum error being 6.4% for Sf3 f3 2:80 at x 0:98. The amounts are smaller because the natural angular frequencies were 1.00062, 1.03736, 1.13783, 1.27174, 1.41421, 1.54359, 1.64479, 1.70994 and 1.73170 which are, as a whole, closer to the
angular frequencies of the original structure than were those of the ®rst two examples. In fact, the second, fourth, sixth and eighth natural frequencies are identical, while the dierences of the other ®ve frequencies are small. The results of Fig. 4(a)±(c) show that even when the structural stiness and mass parameters are not evaluated very precisely, so long as the frequency spectrum of the structure can be measured fairly precisely, the identi®cation can be acceptably accurate. 7.1.6. The in¯uence of the number of participant modes on the identi®cation precision In all the above discussions, all nine modes of the structure were used in the direct analyses, i.e. when computing the response PSD from the excitation PSD, as well as in the back analyses, i.e. when computing the excitation PSD from the response PSD, by the modesuperposition method. Now suppose that all nine modes are still used for the direct analyses needed to obtain the simulated ``measured values'', but that less modes are used for the back analyses. For instance, taking q 7 (modes) and using the ®rst two digits in all the measured values gave the identi®ed excitation PSD curves of Fig. 5(a), which deviate considerably from the curves of Fig. 2(a), for which q 9, particularly in the high frequency region x > 1:4 for which the dierence can be hundreds of percent. Even within the region x 2 [0.8, 1.4], the maximum deviation of the identi®ed values from the theoretical ones is about 30%, Fig. 5(b). It should be
J.H. Lin et al. / Computers and Structures 79 (2001) 375±387
383
Fig. 4. Loading PSD functions of points 3 and 6 identi®ed from the displacement PSD functions of points 2, 5 and 7, with all nine modes participating in the back analysis. All damping ratios were 0.03 and two digits of accuracy were retained for the displacement PSD values.
noted that all the natural frequencies of this structure are closely spaced, the ratio of the highest to the lowest being only 1.732. This is why neglecting only the two highest modes caused considerable in¯uence over the whole frequency region. However, computer simulations for structures with widely spaced natural frequency spectra show that neglecting higher modes aects the
identi®cation precision relatively little in the lower frequency region. 7.2. Example 2: jacket platform Fig. 6 shows the FEM model of a jacket platform which has been built in the Bohai sea of China and has
384
J.H. Lin et al. / Computers and Structures 79 (2001) 375±387
Fig. 5. Loading PSD functions of points 3 and 6 identi®ed from the displacement PSD functions of points 2, 5 and 7, with the ®rst seven modes participating in the back analysis. All damping ratios were 0.03 and two digits of accuracy were retained for the displacement PSD values. The two alternative frequency regions shown were used.
Fig. 6. Jacket platform: only those nodes and members referred to in the text are, respectively, numbered and lettered.
183 beam elements and 492 DOF. Its ®rst 16 natural angular frequencies are 3.8325, 3.8811, 4.1619, 7.9105, 9.9539, 11.205, 11.547, 12.091, 16.182, 17.518, 17.898,
19.161, 21.152, 21.537, 22.777 and 23.160 (sÿ1 ) and the damping ratio of all modes is 0.03. Flowing ice blocks are assumed to cause horizontal x direction stationary random excitations at nodes 5 and 11. The excitation PSD matrix for computer simulation is still taken as Eq. (49). A direct PEM analysis [6] using the ®rst 16 modes was performed to obtain the response PSD matrix for the z direction displacements at points 6, 10, 17 and 23. Its ®rst two eective digits were taken as the ``measured values'' for computer simulation. The loading identi®cation results obtained from the proposed IPEM are shown in Fig. 7(a) and are clearly very poor, i.e. far from the true excitation PSD values of Eq. (49). Therefore, the z direction displacement responses were replaced to obtain Fig. 7(b)±(d), which were obtained from the z displacement responses at the points given in their captions. Clearly, Fig. 7(c) gives the best identi®cation results and they are of usable accuracy. Hence using the z direction displacements at the points 7, 9, 50 and 86 would be a good choice of response to use, if real engineering measurements are made on the actual structure. Suppose that instead the axial strain responses of members A, B and C of Fig. 6 (i.e. beams 9±21, 11±23 and 17±25) are known and that only their ®rst two effective digits are used for loading identi®cation. The results are shown in Fig. 8(a) and are clearly quite good. However, Fig. 8(b) shows that even better results were obtained by including the vertical displacement response
J.H. Lin et al. / Computers and Structures 79 (2001) 375±387
385
Fig. 7. Loading PSD functions of the x-direction forces at nodes 5 and 11, obtained by using two digit accuracy for the displacement PSD values, when the ®rst seven modes participated in the back analysis. The results were identi®ed from the PSD functions of the vertical displacements at the four points shown in the sub-captions.
curve at node 21 in addition to the strains of members A, B and C. In contrast, when the known responses were chosen to be the axial strains of members D and E (i.e. beams 5±17 and 7±19) and the vertical displacements of nodes 50 and 86, the identi®cation results were much poorer, Fig. 8(c), and so this combination should not be
selected for practical measurements on the real structure. However, the results of Fig. 8(d) show that it would only be necessary to replace the vertical displacement response at node 86 by the horizontal displacement response at node 21, to again obtain very good results.
386
J.H. Lin et al. / Computers and Structures 79 (2001) 375±387
Fig. 8. Loading PSD functions of the x-direction forces at nodes 5 and 11, obtained by using two digit accuracy for the displacement PSD values, when the ®rst seven modes participated in the back analysis. The results were identi®ed from the PSD functions of the quantities shown in the sub-captions.
8. Conclusions A random loading identi®cation approach based on a proposed IPEM has been presented and a primary investigation has been given. Previously, such problems have not been solved well. This paper reveals some im-
portant phenomena which may have contributed towards this. For complex engineering structures, the selection of the responses to be used for loading identi®cation is an important problem which has not yet been suciently resolved and for which experience-based methods are
J.H. Lin et al. / Computers and Structures 79 (2001) 375±387
not very reliable. The examples presented show that computer simulation combined with the IPEM is a good tool for selecting measured responses and hence, for identifying the random loadings. This is very important because a real engineering measurement project may take months or even years to prepare and is very costly, whereas computer simulation is comparatively very quick and cheap. For the second example, i.e. the platform, the direct analysis used about 10 min on an IBM/ 486 personal computer and only needs to be done once, while the back analysis took about 15 min for each of the above response combinations, i.e. the computations were very cheap. Therefore, performing sucient computer simulations before practical measurements are undertaken is strongly advisable. Some basic phenomena have been observed and analyzed in the present paper and some suggestions have been given which should be helpful in engineering applications. However, more work remains to be done, e.g. the problem of noise disturbance should be extensively investigated in a stricter way and the bene®ts could be investigated of using Kalman or H1 ®ltering techniques to further improve the identi®cation precision. Acknowledgements The authors are grateful for the support of the National Natural Science Foundation of China (project no. 19772009), the funding of NKBRSF (no. G1999032805) and the State Education Commission of China (no. 97014120), as well as the Advanced Chinese Engineering (ACE) Centre at Cardi University. The authors also wish to thank Bohai Petroleum Company for providing practical data of the platform of example 2 and for their ®nancial support for the research.
387
References [1] Nigam NC. Introduction to random vibration. Cambridge, Massachussets: MIT Press; 1983. [2] Bendat JS, Piersol AG. Engineering applications of correlation and spectral analysis. New York: Wiley; 1980. [3] Zhang JH, Wang C. Theory of engineering random vibration. Xian, China: Xian Jiao Tong University; 1988. [4] Lin JH. A fast CQC algorithm of PSD matrices for random seismic responses. Comput Struct 1992;44:683±7. [5] Zhong WX. Review of a high-eciency algorithm series for structural random responses. Prog Nat Sci 1996;6:257±68. [6] Lin JH, Zhang WS, Li JJ. Structural responses to arbitrarily coherent stationary random excitations. Comput Struct 1994;50:629±33. [7] Ewins DJ. Model testing: theory and practice. London: Research Studies Press; 1980. [8] Farrar CR, James GH. System identi®cation from ambient vibration measurements on a bridge. J Sound Vibr 1997; 205:1±18. [9] Doebling SW, Farrar CR, Goodman RS. Eects of measurement statistics on the detection of damage in the Alamosa Canyon Bridge. Proc Fifteenth Int Modal Anal Conf-IMAC. Orlando, FL, USA, 1997. p. 919±29. [10] Askari MR. System identi®cation using vibration monitoring for fault diagnostics in manufacturing machinery. Proc Fifth Int Conf Factory 2000 ± the technology exploitation process. Cambridge, UK, 1997. p. 295±302. [11] Golub GH, Van Loan CF. Matrix computations. Baltimore: Johns Hopkins University Press; 1983. [12] Cook RD. Concept and application of ®nite element analysis. New York: McGraw-Hill; 1981. [13] Bathe KJ, Wilson EL. Numerical methods in ®nite element analysis. Englewood Clis, NJ: Prentice-Hall; 1976. [14] Wilkinson JH, Reinsch C. Linear algebra, vol. II. Berlin: Springer, 1971. [15] Zhong WX, He Q, Liu ZX. Numerical computational methods. Beijing: Press of Architectural Industry of China; 1991.