Application of finite element model updating to a feed water pipeline of a nuclear power plant

Nuclear Engineering and Design 235 (2005) 1849–1865

Application of finite element model updating to a feed water pipeline of a nuclear power plant A. Veps¨a a,∗ , H. Haapaniemi b , P. Luukkanen c , P. Nurkkala b , A. Saarenheimo a a

b

Technical Research Centre of Finland, Industrial Systems Department, Kemistintie 3, P.O. Box 1704, FIN-02044 VTT, Espoo, Finland Fortum Service Ltd., Condition Management Centre, P.O. Box 10, 00048 Fortum, Vantaa, Finland c Fortum Power and Heat Ltd., Atomitie 700, P.O. Box 23, 07901 Loviisa, Finland Received 18 November 2004; received in revised form 24 May 2005; accepted 30 May 2005

Abstract This article is concerned with application of finite element model updating to a part of a secondary feed water pipeline of a VVER type nuclear power plant. A mathematical model of the pipeline was produced with finite element method. The pipeline was subjected to modal testing in three different conditions. The results of the modal testing were then compared with those predicted with the model. This comparison revealed discrepancies between these two sets of results. It is a commonly known fact that incorrect values of parameters in a model, among other things, do cause these discrepancies. In order to improve the correlation between the measured and the predicted results and hence improve the reliability of the model, the model was subjected to finite element model updating. In this updating process, incorrect values of parameters in a model are adjusted, so that a desired level of correlation between measured and predicted results is achieved. As a result, the correlation between the measured and the predicted results was improved. The bulk of the improvement was achieved by adjusting the values of those parameters that were associated with boundary conditions of the model. © 2005 Elsevier B.V. All rights reserved.

1. Introduction In structural dynamics, as well as in many other fields of engineering, finite element modelling is a standard tool for modelling of behaviour of a system under study. Reliable dynamic analysis is required ∗ Corresponding author. Tel.: +358 020 722 6838; fax: +358 9 456 7002. E-mail address: [email protected] (A. Veps¨a).

0029-5493/$ – see front matter © 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.nucengdes.2005.05.020

particularly when structures are subjected to considerable inertial forces. In this connection, the reliability of an analysis depends on the validity of the model among other things. This validity can be assessed by comparing results predicted by the model with the ones that are measured from a prototype of a structure. If considerable discrepancies can be found between the results, the model cannot be judged as valid. In this case, causes of these discrepancies should be sought and corrected to the best possible extent. Although measured

1850

A. Veps¨a et al. / Nuclear Engineering and Design 235 (2005) 1849–1865

Nomenclature Latin letters a acceleration, correction factor for mass b correction factor for stiffness E modulus of elasticity (Pa) F force (N) f frequency (Hz) H frequency response function J objective function p parameter R response Greek letters  difference ρ mass density (kg/m3 ) ω angular frequency (rad/s) Matrix and vector notations [K] spatial stiffness matrix [KC ] correctly modelled part of a spatial stiffness matrix [M] spatial mass matrix [MC ] correctly modelled part of a spatial mass matrix [S] linear sensitivity matrix [Wp ] parameter weighting matrix [WR ] response weighting matrix [ψN ] natural mode shape matrix {ε} residual vector {ψ} natural mode shape vector {p} vector of changes in parameter values []T transpose of matrix (or vector) {¯a} complex conjugate of complex valued vector Subscripts and superscripts A refers to the mathematical model E refers to the measured data N refers to natural mode of vibration p parameter R response Rel. relative 0 refers to the values at linearization point

Abbreviations BSL best straight line DOF(s) degree(s)-of-freedom FE finite element FRF(s) frequency response function(s) MAC modal assurance criterion MIF mode indicator function NPP nuclear power plant

results do not come without errors, they are usually considered to be more reliable than the predicted results. Because of this, causes of the discrepancies can more probably be found from the model. In the work by Ewins (2000), model updating is defined as “the process of correcting the numerical values of individual parameters in a mathematical model using data obtained from an associated experimental model such that the updated model more correctly describes the dynamic properties of the subject structure.” In the same reference, validation of a model is defined as “the process of . . . attaining the condition that the coefficients in a model are sufficiently accurate to enable that model to provide an acceptably correct description of the subject structures dynamic behaviour.” Updating of finite element models has been a popular field of study in structural dynamics during the past 30 years. Vast majority of the published journal articles are concerned with theoretical aspects of the subject, verified usually with simple academic applications. Applications to industrial scale structures have been discussed in articles in numerous conference proceedings and also in a few journal articles. G¨oge (2003) discussed updating of the finite element model of a large aircraft and Teughels and De Roeck (2004) employed model updating technique for structural damage identification of a highway bridge. Reported applications of model updating to pipelines or other NPP components are few in number. Sinha and Friswell (2003) introduced application of model updating to structural components used in NPPs and discussed its use as a non-destructive fault diagnostic method. The purpose of this research project was to study the usefulness of a combination of finite element modelling and experimental measurements in a vibration

A. Veps¨a et al. / Nuclear Engineering and Design 235 (2005) 1849–1865

1851

related analysis of pipelines. The goal was to be able to compute cyclic stresses in piping material caused by mechanical vibration in operational conditions. One of the tasks on the road was to determine and improve the reliability of predictions given by the model. This article concerns on shedding light into the updating process of the model.

2. Experimental techniques and specimen 2.1. The pipeline The subject of the study was the loop RL61 of the secondary feed water pipeline of a VVER 440 type pressurised water reactor NPP in Loviisa, Finland. The pipeline is made out of mild carbon steel and it is insulated by 120 mm thick layer of mineral wool. The nominal outer diameter of the piping cross section is 323.9 mm. The nominal value of wall thickness is 20 mm for the most part of the loop. On the lower part of the longest vertical section and between the check valves it is 17.5 mm. Nominal radius of curvature of the bends is 600 mm, measured from the centre of the cross section of the pipe. The operating temperature of the pipeline is approximately 165 ◦ C and the corresponding pressure is approximately 7.5 MPa. The loop is almost located in one plane. It is supported by three spring hangers and its total length is approximately 28 meters. The loop begins from a make-up pump and it ends in a T-junction with the main feed water line. The loop is shown in Fig. 1 together with its principal dimensions. Numbered arrows point to main components of the system. These numbered items are as follows: (1) the make-up pump; (2) the distance piece; (3) the back valve; (4) the expansion collar; (5) fixing point of a spring hanger; (6) fixing point of another spring hanger; (7) the first check valve; (8) fixing point of the third spring hanger; (9) the second check valve and (10) the T-junction with the main feed water line.

Fig. 1. Side view of the pipeline RL61.

2.2. Conventional modal testing The pipeline was subjected to conventional modal testing in three different conditions which are described in Table 1. The idea behind the different conditions was to study the effects of insulation, change in the temperature of the pipe, water contents and pressure on the response and modal properties of the structure. The make-up pump of the line RL61 was not operating during the tests but pumps of the adjacent lines were operating during the measurements in Case 3. More thorough description of the testing process may be found in the work by Saarenheimo et al. (2003). Dynamic properties of the pipeline were first determined by means of frequency response functions (FRFs). Responses were measured as accelerations and therefore each particular FRF is defined by the equation Hij (ω) =

ai (ω) , Fj (ω)

(1)

where ai (ω) is a Fourier transformation of the acceleration response signal in the time domain, measured from

Table 1 Testing conditions Case number

Insulation

Water contents

Temperature (◦ C)

Pressure (MPa)

Number of the measurement DOFs

1 2 3

No Yes Yes

No No Stagnant

25–40 25–40 25–165

No No 1–7.5

126 84 93

1852

A. Veps¨a et al. / Nuclear Engineering and Design 235 (2005) 1849–1865

2.2.1. Testing equipment The pipeline was excited with a high capacity Robcon hydraulic shaker system. The shaker system was used in the inertial mode with the inertial mass of 400 kg, maximum nominal dynamic force capacity of 11 kN and stroke of 50 mm. Acceleration responses were measured with Br¨uel&Kjaer type 4370 accelerometers. A Br¨uel&Kjaer IDA 16-channel data acquisition system was used for collecting and recording the data. SDRC I-DEAS 6 signal processing and modal analysis software was used as a signal generation and processing tool on the spot. Version 7 (SDRC, 1998) of the same software was used for the detailed modal analysis at office. A Kemo VBF-8 2-ch High-pass/Low-pass/Band-pass-filter was used for modification of the excitation source signal. Br¨uel&Kjaer type 2635 charge amplifiers with adjustable gain & LPF were also used.

Fig. 2. Locations of the test points in Case 1.

the degree of freedom (DOF) i and Fj (ω) is a Fourier transformation of the excitation force signal in the time domain when the excitation force is applied at the DOF j with no other external forces being applied on the structure at the same time. Modal properties were then computed from the FRFs with suitable modal analysis algorithms. The modal testing was carried out using two different excitation directions. These were the Y- and Z-directions shown in Fig. 2. This enhances the chance that all natural modes at the frequency band of interest are excited during the tests. The measurements were carried out in two separate steps corresponding to these two excitation directions. This procedure was repeated for each of the three testing conditions described in Table 1. As a result, six separate FRF databases were created corresponding to these six testing configurations. In each case, the excitation was applied at the point number 114 in the experimental model. Fig. 2 illustrates locations of the numbered test points in Case 1. When testing in Cases 2 and 3, some of these test points were removed.

2.2.2. Methods The inertial mode of excitation was chosen in order to avoid possible errors caused by introduction of new boundary conditions. This would have been the case if the shaker had been rigidly attached to surrounding structures. However, the shaker was suspended with rubber bands from an overlying concrete beam, so that additional mass introduced by the shaker to the pipeline was minimised. The mass of 400 kg was attached to the hydraulic shaker, so that the inertial forces required for the proper excitation of the pipeline were achieved. Maximum dynamic force that was applied to the structure was approximately 2 kN. The shaker system and the attached mass are shown in Figs. 3 and 4. One measurement set consisted of responses measured from 15 DOFs. Response accelerometers were circulated between each measurement set. All measurement results were directly processed into the form of FRFs. For each FRF, 16 averages were measured. Also the so-called Hv -estimate, the coherence function and the auto power spectra of each FRF were computed and stored. In order to reduce total measurement time, a 50% overlapping was used. With this measurement set up, total measurement time for one set of FRFs was 272 s. Frequency range of the measurements was set to be 0–100 Hz. The number of the spectral lines was set to be 3200 in order to identify closely spaced modes. This resulted in a relatively high frequency resolution of 0.03125 Hz. Broad band random noise was used as an

A. Veps¨a et al. / Nuclear Engineering and Design 235 (2005) 1849–1865

Fig. 3. Hydraulic shaker in the horizontal excitation configuration.

excitation source signal and the Hanning-window was applied to the signals. The Hanning-window weights the beginning and the end of each sample interval to zero and thus enables fast Fourier transformation of the signal. 2.2.3. Results An example of two measured FRFs is shown in Fig. 5. This figure illustrates the effect of the excitation

1853

and response direction on the results. Both of the presented FRFs are from Case 1 with the same excitation point and also the same response measurement point but yet different directions. The curve illustrated with the thicker line represents the modulus part of the FRF with the Y-excitation direction and the response DOF 115Y. The curve illustrated with the thinner line represents the modulus part of the FRF with the Z+ excitation direction and the response DOF 115Z. The figure shows that the natural modes that are excited depend on the excitation and response directions. It can also be seen that the response level is dependent on the excitation and response directions. In Fig. 5, this dependency shows clearest at frequencies lower than 22 Hz. Fig. 6 illustrates the effect caused on the FRFs by the insulation. Both of the presented FRFs have the same excitation and the response DOF. The curve illustrated with the thicker line represents the modulus part of the FRF from Case 1 (with insulation). The curve illustrated with the thinner line represents the modulus part of the FRF from Case 2 (with insulation). The insulation reduced the level of the response at resonance conditions. This reduction was significant especially at frequencies higher than 20 Hz. Generally, the

Fig. 4. Hydraulic shaker and attached mass used to create inertial forces to excite the structure: (a) side view; (b) front view; (c) top view.

1854

A. Veps¨a et al. / Nuclear Engineering and Design 235 (2005) 1849–1865

Fig. 5. An example of two FRFs in Case 1 with different excitation and response directions.

insulation was also noticed to damp out some of the modes found from the measurements with the bare piping. Natural frequencies and corresponding damping values were computed from the FRFs with the socalled Complex Exponential method while natural mode shapes were extracted with the Circle Fit algorithm. The resulting natural mode shapes were complex valued because of phase differences that existed between the different measurement DOFs.

Instead of performing laborious task of checking all the measured FRFs individually, the natural frequencies were identified using the mode indicator function (MIF) diagrams. These diagrams represent the product of the FRFs, selected from the FRF database where global modes are strongly emphasised and therefore also easier to identify. An example of MIF diagrams is shown in Fig. 7. The diagram is from Case 1 with the Y-excitation direction. It can be seen that natural frequencies appearing here as valley points are much

Fig. 6. An example of the effect caused to FRFs by the insulation.

A. Veps¨a et al. / Nuclear Engineering and Design 235 (2005) 1849–1865

1855

Fig. 7. The MIF diagram from Case 1 with the Y-excitation direction.

easier to identify from a MIF diagram than from individual FRFs. After the natural modes were extracted they were verified with modal assurance criterion (MAC) computations. In these MAC computations all the identified mode shapes from the same measurement configuration were set against one another. Then the MAC value was computed for each of the mode shape pairs obtained that way. The result of the MAC computation can be presented in the form of a three dimensional MAC matrix. The MAC value indicates how similar compared mode shapes are. The concept of the modal assurance criterion was first introduced by Allemang and Brown (1982) and it is defined by the equation MAC(1 {ψ}, 2 {ψ}) =

¯ 2 |1 {ψ}T2 {ψ}| , ¯ 2 {ψ}T 2 {ψ}) ¯ (1 {ψ}T 1 {ψ})(

(2)

where 1 {ψ}and 2 {ψ} are the (complex valued) mode shape vectors that are to be compared. The closer the MAC value is to unity, the more similar are the modes that are compared. If the number of the measurement DOFs is sufficient enough and if their locations are correctly selected, the measured mode shapes can be properly distinguished from each other. Perfectly extracted natural mode shapes should form an orthogonal modal basis in which each mode shape is orthogonal to all the others. In this ideal situation, the MAC matrix consists of diagonal terms equal to unity and off-diagonal terms that are zero. In practical situations though, off-diagonal terms are always larger than zero.

Fig. 8. Isometric view of the MAC matrix from Case 1 with the Y-excitation direction.

Fig. 8 shows an isometric view of the MAC-matrix from Case 1 with the Y-excitation direction. Mode shape numbers on both of the axes on the horizontal plane refer now to the same measurement set. Perfect MAC values of 1 are presented in the diagonal with white columns. All the other columns are shaded according to their value with dark shade indicating high MAC value and light shade low MAC value. Only those mode pairs with MAC value higher than 10% are assigned with a column. It can be deduced that some of the closely spaced modes are very similar to each other in this case. However, the vast majority of the modes can be clearly distinguished from the others. The obtained modal parameters were also verified by comparing the measured FRFs with the ones synthesised from the extracted modal parameters. This comparison can be used to detect the overall correctness of the extracted modal parameters. It also reveals whether there exist natural modes that were missed during the extraction process. If modal parameters are extracted correctly, then the synthesised FRF and the measured FRF should be similar in the vicinity of the natural frequencies. An example of this comparison is illustrated in Fig. 9. The FRF in

1856

A. Veps¨a et al. / Nuclear Engineering and Design 235 (2005) 1849–1865

Fig. 9. Comparison of the measured FRF and the corresponding FRF synthesised from the extracted modal parameters.

question is from Case 3, with the Y-excitation direction and the response DOF 115Y−. The curve in the lighter shade represents the measured FRF and the one in the darker shade the corresponding synthesised FRF. It can be seen that the measured and the corresponding synthesised FRF have some discrepancies in this case.

3. Finite element analysis The pipeline was modelled with elements that are supported by the ABAQUS (version 6.4, 2003) software. The same software was used for the computation of the modal properties. Straight runs of the pipeline, the make-up pump and the valves were modelled with B31 beam elements. Element B31 is a two-node, linear Timoshenko beam with exact rotary inertia, lumped mass formulation and six active degrees of freedom at each node. Total number of these elements was 242. Pipe bends were modelled with ELBOW31 elements. ELBOW31 element is a two-node pipe in space with a deforming section and linear interpolation along the pipe. Total number of ELBOW31 elements was 42 with six elements used for each 90◦ bend. ELBOW31 elements were modelled as substructures with warping of

the cross-section disabled at both ends of each element. Spring hangers were modelled with two-dimensional truss elements. Stiffness of each spring was taken into account by adjusting the cross sectional area of each truss-element accordingly. Boundary conditions to the rest of the pipeline were modelled with three point mass elements, three rotary inertia elements and 12 SPRING1 elements acting between the structure and the ground. The model included totally 1770 degrees of freedom. Wire frame presentation of the FE model is shown in Fig. 10. Points 373, 347 and 359 indicate nodal points where the elements used to model boundary conditions locate. Point 336 indicates fixing point of the main feed water line. Antenna-looking components represent the check valves. The behaviour of these check valves was supposed to be that of a rigid body and the antenna-looking construction was chosen so that the centre of mass could be adjusted in the updating process. Mass density of the piping material was initially assumed to be 7850 kg/m3 . The effect of the insulation was taken into account as an equivalent mass by increasing the density of the piping material accordingly. The same method was applied to accommodate the mass of the water contents of the pipeline. Possible

A. Veps¨a et al. / Nuclear Engineering and Design 235 (2005) 1849–1865

1857

meaning that there does not exist phase differences between the DOFs.

4. Model updating

Fig. 10. FE model of the structure.

stiffening effect of the water and the insulation was neglected in the initial analysis. The effect of the increased temperature of the pipeline in Case 3 was taken into account by reducing the modulus of elasticity of the piping material from its value of 201 GPa in the room temperature to the value of 192 GPa in the temperature corresponding to the operational conditions. The stiffening effect of the increased pressure in the piping (Case 3) was neglected in the preliminary analysis. Natural frequencies and mode shapes were computed from the well known, discretised equation of motion in the frequency domain: 2 [[K] − {ωN }[M]][ψN ] = {0},

(3)

2 }}, the vector of squared angular natural where {{ωN frequencies of vibration; [M], the spatial mass matrix; [K], the spatial stiffness matrix and [ψN ], a matrix of natural mode shapes of vibration. Unlike measured mode shapes, predicted mode shapes were real valued,

Each and every mathematical model contains at least some amount of assumed idealisations in the behaviour of the subject structure and uncertainties related to the correct values of the parameters in the model. Because of these factors, it is obvious that the results predicted with the model do not completely match the corresponding correctly measured results. When it comes to a finite element model, discrepancies between the predicted and the correctly measured results arise from three main sources which are: idealisation errors in mechanical behaviour of the system, discretisation errors and erroneous values of the parameters in the model (Mottershead and Friswell, 1993). The goal of the model updating process is not to be able identically reproduce particular set of measured results with the mathematical model, but merely to seek correct values for those parameters in the model which are associated with some amount of uncertainty. It should be noticed that only errors in the values of parameters can be truly corrected with the aid of the measured data. Other types of errors can only be compensated in the updating process and the resulting values of the parameters may differ from their correct values. First historical attempts to update a finite element model were by way of trial and error. There have appeared two computational approaches for the updating problem afterwards. The earlier one is to directly adjust the values of the structural mass, stiffness and/or damping matrices, so that the discrepancies between measured and predicted results are diminished. The posterior approach is to adjust the values of the individual physical parameters of the model. The latter approach generally requires iteration but it is noted to have considerable advantages over the first approach. One of its main advantages is that it enables the user to control that the changes to be made are physically meaningful, which may be difficult to obtain when using the former approach. In this research project updating process, including sensitivity analysis and correlation analysis, was carried out with the aid of FEMtools finite element code (version 3.0.1, 2004).

1858

A. Veps¨a et al. / Nuclear Engineering and Design 235 (2005) 1849–1865

4.1. Initial correlation analysis

Table 2 Initial correlation between the measured and the predicted results

Reliability of the predictions given by the model can be verified by comparing them with the measured results. This comparison process is commonly known as correlation analysis. The results may be compared either visually, usually with the aid of suitable software, or numerically with the aid of indicators developed specially for this purpose. The problem with the correlation analysis with modal properties is that it may be difficult to be sure that the predicted and the measured mode shape do present the same mode. In this case, the predicted and the measured natural mode shapes were compared numerically using the MAC value. If the FE model is perfectly correct and the number of measured DOF’s is adequate, there should be one predicted natural mode shape for each measured one, so that the MAC value of this pair is equal to one and zero for all the other possible combinations. The correlation between the predicted and the measured natural frequency was determined as their relative difference:

Case number

A

B

C

1 (f = 0–10 Hz) 1 2 (f = 0–12 Hz) 2 3

7/8 14/26 6/29 12/29 8/19

0.673 0.620 0.678 0.600 0.56

13.9 16.0 10.3 9.6 7.0

fN,Rel. =

Af

E N − fN Ef N

× 100%,

(4)

where A fN is the predicted natural frequency and E fN , the corresponding measured natural frequency. For this application one has to be sure that the paired natural modes are really related to the same property. In other words they do represent the same mode. The correlation between the measured results and the results predicted with the preliminary model was case-specific. It was found to be the best in Case 1 and at frequencies lower than 20 Hz. This was expected since in Case 1 there was minimum amount of uncertainties associated with the mass and the stiffness parameters of the model. Only a few measured natural modes could be unambiguously associated with a predicted counterpart. In most cases one measured mode correlated with two or more predicted modes and vice versa. The correlation between the predicted and the measured mode shapes was particularly poor in Case 3. There were doubts that the insulated pipe may have been in contact with the cable rack nearby. If these doubts hold true, it would explain at least some of the lack of correlation in Cases 2 and 3. The water

A: number of the established natural mode pairs out of the maximum possible; B: average of the MAC values of the established mode pairs; C: average of the absolute values of the relative frequency differences of the established mode pairs (%).

contents in Case 3 is also a potential reason for the lack of correlation. Modelling of the water contents as an equivalent mass may not correctly represent its true behaviour. Table 2 summarises the correlation between the measured results and those predicted with the preliminary model for different cases. The table shows the number of those measured natural mode shapes out of all possible which were found a predicted counterpart. It also shows the average of the MAC values and the average of the absolute values of frequency differences of these pairs. It should be noticed that while the level of correlation between the predicted and the measured results decreases in terms of the MAC values, it increases in terms of differences between the natural frequencies. The measured and the predicted natural frequencies of the paired natural mode shapes are illustrated as a co-ordinate plot in Fig. 11. The figure also shows the best straight line in each case and the line of perfect correlation. In each case, this best straight line through the given data points diverges from the line of the perfect correlation such a way that the predicted frequencies are higher than the measured ones. This indicates that either the mass of the structure was underestimated or the stiffness of the structure was overestimated or both. In Case 3, interpretation of the results should be taken with some reservation due to the small amount of samples. 4.2. Computational updating algorithm The so-called inverse eigensensitivity method was used as the updating algorithm. It is an iterative method that has its basis on the assumption that the differences between measured and predicted vibration properties

A. Veps¨a et al. / Nuclear Engineering and Design 235 (2005) 1849–1865

1859

Fig. 11. The correlation between the measured natural frequencies and those predicted by the model in the initial configuration.

can be described in terms of the relevant sensitivities and small adjustments to selected mass, stiffness and possibly damping elements in the model. These sensitivities can be described as rates of change of modal or response properties with respect to changes in individual mass, stiffness and damping terms. Instead of making adjustments directly to those terms, adjustments are made to the geometrical and material related parameters of the FE model. That way the direct relation between adjustments and physical properties is maintained (Ewins, 2000). For the purpose of the updating, mass matrix is partitioned as follows: [M] = [MC ] +

L


as [Ms ],

(5)

s=1

where [MC ] is the spatial mass matrix of those parts that are assumed to be modelled correctly; L, the number of elements (or element groups) that are associated with uncertainties in the mass modelling; as , the mass correction factor for the sth such element (or element

group) and [Ms ] is the spatial mass matrix of the sth such element (or element group). Corresponding partition for the stiffness matrix is [K] = [KC ] +

N


bs [Ks ],

(6)

s=1

where [KC ] is the spatial stiffness matrix of the parts that are assumed to be modelled correctly; N, the number of elements (or element groups) that are associated with uncertainties in the stiffness modelling; bs , the stiffness correction factor for the sth such element (or element group) and [Ks ], the spatial stiffness matrix of the sth such element (or element group). For the sake of simplicity, correction factors a and b are hence referred to just as parameters and denoted with p. The question is then the one of finding suitable parameter changes that will yield a minimum value for an objective function. This objective function represents the quality of the model’s predictions and it is usually formulated as a function of differences between the predicted and the measured results, called hence as

1860

A. Veps¨a et al. / Nuclear Engineering and Design 235 (2005) 1849–1865

responses. Vector of these differences {ε}, called hence as residuals, are given by the equation

updating process. This yields objective function in the form

{ε} = {E R} − {A R(p)},

J = {ε}T [WR ]{ε},

(7)

where {E R} is the vector of the measured responses and {A R(p)} is the vector of the corresponding predicted responses, as a function of parameters. The predicted vector of responses, {A R(p)}, is generally a non-linear function of the updating parameters. In the inverse eigensensitivity method, the relationship between the parameters and the responses is approximated with the Taylor’s polynomial, truncated after the first term as follows: {A R(p)} ≈ {A R(p0 )} + [S0 ]({p} − {p0 }),

(8)

where [S0 ] is the sensitivity matrix at linearization point; {p}, the vector of the parameters and {p0 } is the vector containing the values of the parameters at linearization point. Use of this equation requires that both the function of the responses as well as each of the partial derivatives in the sensitivity matrix is continuous in a closed interval containing the values p0 . It should be noticed that this equation approximates the exact solution the better the smaller are the changes in the parameter values. The sensitivity matrix at linearization point can be written as   ∂Ri [S0 ] = [Sij,0 ] = , (9) ∂pj,0 where partial derivatives can be computed analytically or approximated with any suitable numerical technique. The simplest of the objective functions to be minimised is the sum of the squares of the residuals: J = {ε}T {ε}.

(10)

The residual vector may contain different types of residuals. These residuals may be measured with different accuracy, such as natural frequencies and corresponding mode shapes. They may also have different number of elements associated with them. For example one mode shape is associated with one natural frequency but numerous modal displacement values corresponding to various DOFs. Therefore, different types of residuals should be weighted differently in the

(11)

where [WR ] is the weighting matrix for the residuals. This weighting matrix is usually assumed to be diagonal and the weighting factors WR are determined by the user. The updating problem may be ill-conditioned by means that small changes in the given values yield big changes in the results. In that case the problem is sensitive to disturbances, such as noise, and inaccuracies in the measured data. Another problem occurs when there are more parameters to be updated than there are equations from which the parameter changes can be computed. In that case, correct solution for the updating problem cannot be found. One way to improve the conditioning of the problem is to include additional information into the set of equations, so that the effect of the errors is minimised. This additional information is referred to as the side constraint and the problem is said to be regularised. However, it should be noted that regularisation also introduces a new source of error, called as regularisation error, into the problem. One side constraint is such that it seeks a small value for the weighted sum of the squares of the parameter changes. Some of the selected parameter values are usually associated with more uncertainty than the others. Therefore, not only the responses but also the parameters should be associated with separate weightings. Relatively accurate parameter values should be weighted more heavily than less accurate ones. This ensures that the changes in these relatively accurate parameter values are kept at a low level. Incorporating the sum of the weighted squares of the parameter changes into the objective function as a side constraint yields objective function in the form J = {ε}T [WR ]{ε} + {p}T [Wp ]{p},

(12)

where {p} is the vector of changes in parameter values in one iteration step and [Wp ] is the weighting matrix for the parameters. It should be noticed that the absolute values of the weighting factors have an effect on which part of the objective function is emphasised more in the updating process, the part with the residuals or the part with the parameter changes.

A. Veps¨a et al. / Nuclear Engineering and Design 235 (2005) 1849–1865

The solution for the parameter changes which minimises the objective function (Eq. (12)) within each iteration step is obtained from the equation −1

{p} = [[S]T [WR ][S] + [Wp ]]

[S]T [WR ]{ε}.

(13)

4.3. Selection of the residuals The choice of the types of residuals to be used in the updating is an important one. Selected responses should be those for which there exists no confusion whether the measured and the predicted data really represent the same item. Frequency response functions have the benefit over the modal properties in that they are specified with the excitation and the response DOF and therefore pairing of the FRFs is relatively easy. Predicted and measured modal properties are more difficult to pair together. Employment of the FRFs enables comparison between the predicted and the measured results at fulllength of the measured frequency range, not just at the resonance frequencies. Other benefits include the possibility of the catering for damping and avoidance of the errors arising from the experimental modal analysis. However, updating with the FRFs may be more tedious and computationally expensive than updating with the modal properties, depending on how many frequency points are included into the analysis. Regardless of the type, selected residuals should be such that they are affected by the changes in values of the selected parameters. It was decided that the model should be computationally updated against the modal data measured from Cases 1 and 2. The residual vector is now given by the equation     {fN,Rel }   {ψij } {ε} = , (14)     {1 − MAC} where {fN,Rel } is the vector of relative frequency differences of the selected mode shape pairs; {1-MAC}, the vector of MAC value residuals of the selected mode shape pairs and the differences in selected nodal displacements are given by the equation ψij = A ψij − E ψij ,

(15)

1861

where A ψij is the predicted value for the amplitude of the ith natural mode at the jth DOF and E ψij , the corresponding, equally scaled, measured value. Selection of the residuals in Case 1 contained 17 mode shape pairs with the MAC values, relative frequency differences and 35 selected nodal displacement differences after the manual adjustment of the boundary conditions. Two additional constraints were introduced in order to keep total masses of both check valves constant during the iteration. Total number of selected responses was thus 71. In measurement Case 1 total number of 126 DOFs was used to measure the responses. Therefore, each modal displacement was assigned with 1/126 of the weight of a MAC value of one mode shape pair. Being less accurately measured, each MAC value was assigned with 3/4 of the weight of the corresponding natural frequency. 4.4. Selection of the parameters to be updated The success of the updating also depends on the selection of the parameters to be updated. The correct selection is in turn dependent on the experience and engineering judgement of the analyst. Generally parameters to be adjusted should be those that are both associated with uncertainties and have reasonably high ability to change the dynamic behaviour of the structure. Furthermore, if unique solution is desired with physically meaningful changes in the structure, the number of the selected parameters should not exceed the number of the residuals used in the updating. The level of uncertainty that is associated with the model is shown in Table 3. The level of uncertainty Table 3 Uncertainties in the model Component

1

Boundary conditions Make up pump Check valves Modelling of the water contents Back valve Distance piece Spring supports Pipe bends Piping material Modelling of the insulation Straight pipe sections

X X

2

3

4

5

X X X X X X X X

1: highest level of uncertainty; 5: lowest level of uncertainty.

X

1862

A. Veps¨a et al. / Nuclear Engineering and Design 235 (2005) 1849–1865

assigned to each of the components is based on the available information about the pipeline and the engineering judgement of the authors. It is a good practise to perform a preliminary sensitivity analysis with all potential parameters and responses before making the final selection. In this preliminary sensitivity analysis, the effect of the parameters on the responses can be assessed and ineffective parameters and responses can then be excluded from the final selection. Mass density and modulus of elasticity of the material of all possible elements complemented with the boundary conditions were chosen for the parameters for this preliminary sensitivity analysis. The sensitivities were computed as normalised sensitivities, indicating the percentage change of the response for one percent change of the parameter value. It should be noticed that the result of the sensitivity analysis is always related to the current values of the parameters of the model. When the values of the parameters change, also the result of the sensitivity analysis will change. The sensitivity analysis revealed that the most effective parameters in the model were the nodal masses used to model the boundary conditions. This particular mass affects at node 347, which can be located from Fig. 9. Another highly efficient parameter for the updating was found to be the stiffness of one of the spring hangers. The modulus of elasticity of the elements used for modelling short foreparts of the two other feed water lines that diverge from the main line was also found to be effective. In total, stiffness-related parameters were found to be more effective than the mass-related. Parameters associated with the boundary conditions were decided to be adjusted manually. Parameters associated with the make-up pump and the check valves were then adjusted with the computational algorithm. Description of these computationally adjusted parameters is given in Table 4. The number of the

Fig. 12. Parameter sensitivities with respect to natural frequencies.

Fig. 13. Parameter sensitivities with respect to nodal displacements.

parameters to be adjusted is lower than the number of the selected residuals and the problem is therefore overdetermined. Figs. 12–14 show sensitivities of the selected parameters with respect to natural frequencies, nodal displacements and MAC values. For the purpose of this presentation, the absolute value of the normalised sensitivity was computed for each possible

Table 4 Final selection of the parameters for the computational updating Parameter number

Parameter definition

Associated component

1–9 10–18 19–40 41–62

E ρ ρ ρ

The make-up pump The make-up pump The first check valve The second check valve

Fig. 14. Parameter sensitivities with respect to MAC values.

3 2

0.615 0.655 0.695 0.623 0.674 0.757 0.795 0.736

2 (f = 0–12 Hz) 1

0.659 0.600 0.665 0.638 0.751 0.693 0.767 0.731

1 (f = 0–10 Hz) 3

10/19 (+2) 8/19 (0) 10/19 (+2) − 10/29 (−2) 9/29 (−3) 12/29 (0) 13/29 (+1)

2 2 (f = 0–12 Hz)

5/12 (−1) 4/12 (−2) 6/12 (0) 6/12 (0) 13/26 (−1) 15/26 (+1) 17/26 (+3) 16/26 (+2)

1 1 (f = 0–10 Hz)

Case number

6/8 (−1) 6/8 (−1) 8/8 (+2) 8/8 (+2)

Average of the MAC values of the established mode pairs Number of the established mode pairs out of the maximum possible

Parameter values of the boundary conditions proved to have the greatest individual effect on the level of correlation between measured and predicted results. Of all the boundary condition related parameters, point masses proved to have the greatest influence. Spring elements did also have a noticeable effect. Manual adjustment of these values diminished the discrepancies between the measured and the predicted natural mode shapes. However, this yielded only local optimums rather than the global optimum value that was sought. Most of the improvement was achieved rather by finding more suitable mode pairs than improving the correlation between the existing pairs. Additionally some new sufficiently well correlated mode pairs were obtained in Cases 2 and 3. Unfortunately, the difference between the predicted and the measured mode shapes increased somewhat. It should be noted that the same adjustment did not produce the best correlation at all the cases. It was also noticed that modelling of the connection with the main feed water line as completely rigid yielded the best result in terms of the level of correlation when only the mode pairs common to all modelling configurations were taken into account. However, in that case the number of the established mode pairs was lower than in the initial model configuration. The level of correlation between the predicted and measured modal results in four different model configurations is summarised in Tables 5 and 6. In Cases 1 and 2 the results are also presented in the low frequency range. Figures in the parentheses indicate improvement/impairment compared with the initial model configuration. Configuration number one is such that the connection with the main feed water line is modelled as completely rigid. In configuration

Table 5 Level of correlation between the predicted and the measured natural mode shapes in different model configurations

5. Results and discussion

Model 1 Model 2 Model 3 Model 4

parameter–response pair and then summed across all the responses of the same type and normalised to have a maximum value of 1. Therefore, only the relative sensitivity of each parameter is shown in proportion to the others. It can be seen that parameters related to the check valves are generally more effective than those related to the make-up pump. It can also be seen that mass of the make-up pump do not have a noticeable effect on any of the selected responses.

1863 0.497 0.559 0.552 –

A. Veps¨a et al. / Nuclear Engineering and Design 235 (2005) 1849–1865

1864

A. Veps¨a et al. / Nuclear Engineering and Design 235 (2005) 1849–1865

Table 6 Level of correlation between the predicted and the measured natural frequencies in different model configurations Average of the absolute values of the relative frequency differences (%) Case number

Model 1 Model 2 Model 3 Model 4

1 (f = 0–10 Hz)

1

2 (f = 0–12 Hz)

2

3

21.3 14.7 14.5 17.5

19.5 16.8 16.6 17.5

19.2 11.0 11.8 14.5

15.8 11.5 11.0 13.0

7.5 7.3 8.9 –

number two, no elements were used to model the boundary conditions. In configuration number three, boundary conditions were adjusted manually. In configuration number four, the model is computationally updated from the configuration number three. As can be seen computationally updated model did not yield any improvement in this case. The presented results were obtained after a few iteration steps. The slight divergence of the results may be due to relatively high number of parameters selected to be updated.

Fig. 15 shows the overlay of the measured natural mode shape with its predicted counterpart before and after the manual adjustment of the boundary conditions. The measured mode shape, presented by the line with dots, is from Case 1 with the natural frequency of 6.67 Hz. The MAC value of the pair is 0.809 before and 0.866 after the adjustment. Relative frequency differences are 22.3% and 28.6%, respectively. Improved correlation can be seen especially near the main feed water line.

Fig. 15. One measured natural mode shape with its predicted counterpart: (a) before the updating; (b) after the updating.

A. Veps¨a et al. / Nuclear Engineering and Design 235 (2005) 1849–1865

6. Conclusions As an application of updating, the pipeline proved to be more difficult than what was first expected. As in case with many industrial applications, the results could not be verified experimentally. The verification of the results is then only based on the engineering judgement of the analyst. In this connection, results do not mean the correlation between the predicted and the measured results but the values of the parameters after the updating process. Modelling of the boundary conditions seems to have a high influence on the level of correlation between the measured and the predicted results in a case of a pipeline when only a part of the pipeline is modelled. Therefore, it may be useful to concentrate on the boundary conditions when updating the model of a pipeline. However, when the application has other boundary conditions than free/free or completely rigid, the process of finding correct values for some uncertain parameters in the model may be hampered. The adjustment of all the other parameter values may just compensate the simplifications made in the modelling of the boundary conditions. This could be a potential subject of a further study.

Acknowledgements This study has been performed as a part of the project on Structural Operability and Plant Life Management (RKK). The project funding by a joint Finnish

1865

industry group including; the National Technology Agency (Tekes), Teollisuuden Voima Ltd. (TVO), Fortum Power and Heat Ltd., Fortum Oil Ltd. and FEMdata Ltd. is gratefully acknowledged.

References ABAQUS Version 6.4, 2003. User’s Manual. Hibbit, Karlsson & Sorensen Inc. 1080 Main Street, Pawtucket, Rhode Island, New York, USA. Allemang, R.J., Brown, D.L., 1982. A correlation coefficient for modal vector analysis. In: Proceedings of the 1st International Modal Analysis Conference, Orlando, FL, USA, pp. 110–116. Ewins, D.J., 2000. Modal Testing: Theory, Practice and Applications, second ed. Research Studies Press Ltd., Baldock, Hertfordshire, England, 562 pp. FEMtools Version 3.0.1, 2004. User’s Manual. Dynamic Design Solutions NV (DDS), Leuven, Belgium. G¨oge, D., 2003. Automatic updating of large aircraft models using experimental data from ground vibration testing. Aerospace Sci. Technol. 7, 33–45. Mottershead, J.E., Friswell, M.I., 1993. Model updating in structural dynamics: a survey. J. Sound Vibration 167, 347–375. Saarenheimo, A., Haapaniemi, H., Luukkanen, P., Nurkkala, P., Rostedt, J., 2003. Testing for modal analysis of a feed water pipeline, In: Transactions of the 17th International Conference on Structural Mechanics in Reactor Technology, Prague, Czhech Republic, Paper No. J05-7. SDRC I-DEAS ms7, Exploring I-DEAS TestTM manual, 1998. Strcutural Research Corporation, OH, USA. Sinha, J.K., Friswell, M.I., 2003. The use of model updating for reliable finite element modelling and fault diagnostics of structural components used in nuclear plants. Nucl. Eng. Des. 223, 11–23. Teughels, A., De Roeck, G., 2004. Structural damage identification of the highway bridge Z24 by FE model updating. J. Sound Vibration 278 (3), 589–610.