Interval model updating with irreducible uncertainty using the Kriging predictor

Mechanical Systems and Signal Processing 25 (2011) 1204–1226

Contents lists available at ScienceDirect

Mechanical Systems and Signal Processing journal homepage: www.elsevier.com/locate/jnlabr/ymssp

Interval model updating with irreducible uncertainty using the Kriging predictor Hamed Haddad Khodaparast , John E. Mottershead, Kenneth J. Badcock School of Engineering, University of Liverpool, Liverpool L69 3GH, United Kingdom

a r t i c l e in fo

abstract

Article history: Received 21 April 2010 Received in revised form 18 August 2010 Accepted 15 October 2010 Available online 28 October 2010

Interval model updating in the presence of irreducible uncertain measured data is defined and solutions are made available for two cases. In the first case, the parameter vertex solution is used but is found to be valid only for particular parameterisation of the finite element model and particular output data. In the second case, a general solution is considered, based on the use of a meta-model which acts as a surrogate for the full finite element mathematical model. Thus, a region of input data is mapped to a region of output data with parameters obtained by regression analysis. The Kriging predictor is chosen as the meta-model in this paper and is found to be capable of predicting the regions of input and output parameter variations with very good accuracy. The interval model updating approach is formulated based on the Kriging predictor and an iterative procedure is developed. The method is validated numerically using a three degree of freedom massspring system with both well-separated and close modes. A significant advantage of Kriging interpolation is that it enables the use of updating parameters that are difficult to use by conventional correction of the finite element model. An example of this is demonstrated in an experimental exercise where the positions of two beams in a frame structure are selected as updating parameters. & 2010 Elsevier Ltd. All rights reserved.

Keywords: Model updating Interval analysis Irreducible uncertainty Kriging predictor Parameter vertex solution

1. Introduction There is increasing interest in improving finite element predictions for industrial-scale structures. Deterministic finite element model updating approaches [1,2] to improving the accuracy of finite element estimates based on inaccurate and variable measured data, are now well-known and widely used in many industries. In most cases experimental variability is supposed not to be inherent in the test structure itself, but arises from other sources such as measurement noise, the use of sensors that affect the measurement or signal processing that might introduce bias. Such variability, termed epistemic, is reducible by increased information, e.g. [3]. However, manufacturing and material variability in structures is not reducible and must be considered as part of the model, e.g. [4]. This irreducible uncertainty is termed aleatory. Research on model updating is mostly deterministic in that each of the updating parameters is considered to have one ‘true’ value and the purpose of the updating process is to provide an estimate of it. In reality nominally identical structures, built to the same design specification, are different, and this difference should be represented in the model as parameter variability. The purpose of model updating then becomes the estimation of ranges or distributions of parameters. Many different methods have emerged to model such variability, generally categorised in two groups: (i) probabilistic and (ii) non-probabilistic. Probability theory and random fields may be used for probabilistic modelling while interval and fuzzy  Corresponding author.

E-mail addresses: [email protected], [email protected] (H.H. Khodaparast). 0888-3270/$ - see front matter & 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.ymssp.2010.10.009

H.H. Khodaparast et al. / Mechanical Systems and Signal Processing 25 (2011) 1204–1226

1205

sets are used to represent non-probabilistic uncertainty. Many different methods have been proposed to quantify the variability of output data from a dynamic system due to known input variability, e.g. [5–9]. This is known as direct, or forward, propagation and can be used, for example, when the inputs take the form of measurable parameter ranges or distributions, such as the dimensions of a beam cross-section. However, uncertain parameters such as damping and stiffness in mechanical joints are usually not measurable. In this situation, a collection of modal test data from a set of structures with nominally identical joints can be used to identify the distribution of these immeasurable joint parameters. In other words, the ranges or distributions of the updating parameters are estimated from measured variability in output data, typically natural frequencies and mode shapes. This is called uncertainty identification or more specifically stochastic model updating, which is the subject of this paper. Statistical techniques have been used before in model updating of structures when reducible uncertainty, arising from measurement noise, is present in the data. These techniques are not capable of predicting the physical ranges or distributions of variable updating parameters. Only one value for each of the updated parameters, typically the parameter with maximum probability, is physically meaningful and the identified ranges or distributions of updating parameters are indicators of confidence in the estimate. For example, if the number of measurements is increased then the effect of measurement noise tends to be reduced by averaging and so the distribution is smaller and confidence in the estimated parameters is increased. The minimum variance methods of Collins et al. [10] and Friswell [11], the Bayesian methods of Beck and his colleagues [12,13], random matrix theory methods of Soize et al. [14] and inverse fuzzy arithmetic of Haag et al. [15] are all examples of model updating methods dealing with reducible uncertainty. In another updating method the modelling uncertainty is included as parametric random fields [16]. This work is concerned with parameterisation which is an important subject [17–22], but not the direct concern of this paper. Existing methods which deal with irreducible uncertainty in measured data include maximum likelihood estimation and perturbation methods. Fonseca et al. [23] proposed an optimisation procedure for the purpose of stochastic model updating based on maximising a likelihood function and applied it to a cantilever beam with a point mass at an uncertain location. Mares et al. [24] adapted the method of Collins et al. [10] within a gradient-regression formulation for the treatment of teststructure variability. Hua et al. [25] used perturbation theory in the problem of test-structure variability. The predicted output means and covariances were made to converge upon measured values and in so doing the first two statistical moments of the uncertain updating parameters were determined. Khodaparast et al. [26] developed two perturbation methods for stochastic model updating. The first method required only the first order sensitivity matrix and therefore was computationally efficient compared to the method developed by Hua et al. [25] which needed the calculation of the second order sensitivity matrix. Khodaparast and Mottershead [27], and also Govers and Link [28], proposed an objective function for the purpose of stochastic model updating. The objective function consisted of two parts: (1) the Euclidean norm of the difference between mean values of measured data and analytical output vectors, and (2) the Frobenius norm of the difference between the covariance matrices of the measured data and analytical outputs. So far in this discussion of stochastic model updating we have focussed our attention on probabilistic models. These approaches usually require large volumes of data, which are expensive and computational costs are likely to be high. Probably a better approach, given that large quantities of test data will not be produced, would be to use an interval model for the uncertain parameters. The performance, efficiency and accuracy, of the interval modelling and analysis has been tested by the present authors [29,30] in an application of forward propagation in aeroelasticty. In the present article the problem of interval model updating in the presence of uncertain measured data is defined. The solution to this problem can be developed based on existing forward propagation methods for uncertain models represented by intervals. The approaches proposed for the interval analysis may be categorised in three groups: interval arithmetic, vertex solutions and global optimisation [6]. The interval arithmetic, explained in [31], includes a set of operations (addition, subtraction, multiplication and division) which can be used for estimation of the output data bounds. However, the method often over-estimates the bounds of output data because the correlation between the operands is ignored. In the vertex method, originally developed by Dong and Shah [32], the responses are either monotonically increasing or decreasing with the uncertain parameter variations. This implies that the bounds of output data can be sought in all possible combinations of the boundary values of the input intervals. However, the vertex solution is valid only if there is a monotonic relationship between the inputs and outputs. Hanss [9] proposed a transformation method that makes use of the vertex solution for the estimation of fuzzy membership functions for uncertain outputs. He showed that the transformation method also works well if the outputs are non-monotonic with respect to only one of the p uncertain input parameters; this may be useful for a restricted class of dynamic problems as will be explained in the present paper. The third and probably the best solution for interval analysis is based on the application of global optimisation procedures. This method has been used frequently in the literature for the purposes of interval analysis. However, the application of global optimisation techniques often results in increased computational time. This problem is largely overcome by the use of a meta-model as a surrogate for the physical model. It is known that the problem of interval model updating can be solved by using the parameter vertex solution [33] when (i) the overall mass and stiffness matrices are linear functions of the updating parameters, (ii) can be decomposed into nonnegative-definite substructure mass and stiffness matrices and (iii) the output data are the eigenvalues of the dynamic system. Two recursive updating equations may be developed to update the bounds of an initial hypercube of updating parameters in this case. However, the parameter vertex solution is not available generally when, for example, the output data include the eigenvectors of the structural dynamic system and the system matrices are nonlinear functions of the updating

1206

H.H. Khodaparast et al. / Mechanical Systems and Signal Processing 25 (2011) 1204–1226

parameters. In that case optimisation procedures need to be used within the inverse approach. Meta-models have been used for the estimation of bounds of output parameters and are efficient in comparison with traditional optimisation methods. Amongst existing meta-models, such as the conventional response surface method [34], neural networks [35] and Kriging models [36], the latter are found to be excellent predictors of numerical model behaviour (Ghoreyshi et al. [37] and Munck et al. [38]). The main advantage of the Kriging predictor is that it is an unbiased estimator at the observation points. In this paper the general case is solved by using a meta-model so that the region of input data is mapped to the region of output data with parameters obtained by regression analysis. In this procedure, the ‘centre’ of the initial hypercube of updating parameters is identified by deterministic model updating using the mean value of the measured data. Then, an initial hypercube is assumed about the updated point and mapping is carried out by constructing an approximation based on the data from the mathematical/numerical model. This meta-model together with the measured samples can be used to update the initial hypercube of parameters to find the best region of updating parameters corresponding to output data variation. Selection of the meta-model is a crucial step in that it influences the performance of the updating procedure to a very significant degree. The Kriging predictor is shown to be capable of predicting the region of input and output parameter variations with very good accuracy, even in the difficult case of close modes. The Kriging model is generalized for a multi degree of freedom system with ‘n’ output parameters. Then the interval model updating approach is formulated and an iterative procedure is developed. The method is validated numerically by using a three degree of freedom mass-spring system with well-separated and close modes. Finally the method is applied to a frame structure with uncertain internal beams locations. The frame consists of two internal beams, each of which can be located at three different positions. Therefore nine sets of measured data corresponding to each different combination of the beams positions are available. Detailed finite element models of the frame structure with different locations of the beams and a Kriging model describing the relationship between the natural frequencies and updating parameters (beams positions) are provided. The procedure of interval model updating, incorporating the Kriging model, is used to identify the locations of the beams for each case and to update the bounds of beams positions based on measured data. The method successfully identifies the locations of the beams using six measured frequencies for all nine cases with a maximum error of 11%. The updated bounds are found to be in good agreement with the known real bounds on the position of the beams. This example illustrates a very useful advantage of Kriging interpolation that enables the use of nodal coordinates as updating parameters. By conventional methods this would be time-consuming and inelegant, requiring re-meshing of the finite element model at each updating iteration.

2. Problem definition In deterministic model updating, the following perturbation equation is often used, zm zl ¼ Sl ðhl þ 1 hl Þ

ð1Þ

where, zm ¼ ½o21

ðmÞ

ðmÞ

ðmÞ

o22    o2r1

ðmÞ

/T1

ðmÞ

/T2

ðmÞ

   /Tr2 T 2 Rn

ð2Þ

is the assembled vector of measured data, zl ¼ ½o21

ðaÞ

ðaÞ

ðaÞ

ðaÞ

o22    o2r1 /T1

ðaÞ

ðaÞ

/T2    /Tr2 T 2 Rn

ð3Þ

is the lth assembled vector of the predicted outputs from an analytical/numerical model, hl 2 Rp is the vector of updating parameters at lth iteration, Sl is the sensitivity matrix at the lth iteration, o2i is the ith eigenvalue of dynamic system, /i is the ith eigenvector of dynamic system and n is the number of responses. In general, the predicted and experimental mode shapes will be scaled differently. The modal scale factor (MSF) [39] may be used to scale the measured mode shapes to conform to the analytical modes. In the presence of irreducible uncertainty, Eq. (1) has an equivalent in the form of a system of interval linear equations as z~ m z~ l ¼ S~ l ðh~ l þ 1 h~ l Þ

ð4Þ

where ~ represents uncertain vector/matrix terms, modelled as intervals. Seif et al. [40] provide a closed-form solution of Eq. (4) when the sensitivity matrix is square and uncertainty is present either in the outputs zm or sensitivity matrix S, but not together at the same time. However, in interval model updating Eq. (4) is either overdetermined or underdetermined with a non-square sensitivity matrix. The predictions zl, the system parameters h and the sensitivity matrix Sl are all interval vectors/ matrices if the uncertain measured data zm is an interval vector. In this situation, a closed-form solution generally does not exist. The solution of Eq. (4) is given for two cases in the following section. In the first case, the relationship between the inputs and outputs are monotonic, hence parameter vertex solution is applicable. In most cases, however, this desirable monotonic behaviour is not necessarily present and then there is a need to consider other solutions based on optimisation procedures.

H.H. Khodaparast et al. / Mechanical Systems and Signal Processing 25 (2011) 1204–1226

1207

3. Solution methods 3.1. Case 1: parameter vertex solution In this case, the global mass and stiffness matrices are expanded as linear functions of the updating parameters, M ¼ M0 þ

p1 X

mj Mj

ð5Þ

j¼1

K ¼ K0 þ

p2 X

kj Kj

ð6Þ

j¼1

where M is the global mass matrix, K is the global stiffness matrix, mj is the updating parameter for the jth substructure mass matrix, Mj, and kj is the updating parameter for the jth substructure stiffness matrix, Kj. The Young’s modulus and mass density of a substructure are examples of kj and mj, respectively. The decompositions in Eqs. (5) and (6) are non-negative decompositions of the mass and stiffness matrices [33] because the substructure matrices are all semi-positive definite. The eigenvalue derivatives of the global system with respect to the mass and stiffness updating parameters may be written as [41] @o2i @K ¼ /Ti / ¼ /Ti Kj /i @kj i @kj

ð7Þ

@o2i @M ¼ o2i /Ti / ¼ o2i /Ti Mj /i @mj i @mj

ð8Þ

From Eqs. (7) and (8), it can be seen that the signs of the derivatives of the eigenvalues with respect to the updating parameters do not change within the interval of variation ½h h . Therefore, the eigenvalues of the dynamic system increase monotonically with the stiffness parameters and decreases monotonically with the mass parameters. Consequently, two recursive equations can be defined to update the initial hypercube of updating parameters based on the vertices of measured data as z m ¼ z l þ Sjhl,z ðhl þ 1,z m hl,z m Þ

ð9Þ

z m ¼ z l þ Sjhl,z ðhl þ 1,z m hl,z m Þ

ð10Þ

m

m

where  and  represent the lower and upper bounds of , respectively, zm and zl are the assembled vector of measured data and predictions from an analytical/numerical model, hl,z m ¼ ½k l m l T and hl,z m ¼ ½k l m l T and Sjhl, is the sensitivity matrix evaluated at hl, and l is the iteration number. The parameter vertex solution is not necessarily valid when either the output data includes the system eigenvectors or the mass and stiffness matrices are not linear functions of the updating parameters or both. Firstly, let us restrict the output data to be eigenvalues of the dynamic system with mass and stiffness matrices that are nonlinear functions of the updating parameters decomposed as M ¼ M0 þ

p X

fj ðyj ÞMj

ð11Þ

j¼1

K ¼ K0 þ

p X

gj ðyj ÞKj

ð12Þ

j¼1

where fj and gj are nonlinear functions of the updating parameters. For example, if the updating parameter is the thickness of 3 shell element then fj ðyÞ ¼ y and gj ðyÞ ¼ y . In this case, the derivatives of the eigenvalues with respect to the updating parameters may be written as @o2i @gj @fj ¼ kj mj @yj @yj @yj

ð13Þ

where kj ¼ /Ti Kj /i and mj ¼ o2i /Ti Mj /i . The sign of the eigenvalue derivatives with respect to the updating parameters do not necessarily remain unchanged within the region of variation of the updating parameters. Certainly, if the solution of

kj @@gyjj mj @@fyjj ¼ 0 lies within the range of the updating parameters then the parameter vertex solution is no longer valid. Now let us consider the case in which the output data includes the eigenvectors of the dynamic system. The eigenvector derivatives with respect to the updating parameters may be written as [41] n X @/i ¼ aik /k @yj k¼1

ð14Þ

1208

H.H. Khodaparast et al. / Mechanical Systems and Signal Processing 25 (2011) 1204–1226

where 8   > T @K 2 @M > /  o / > k i > > @yj @yj i > < o2i o2k aik ¼ > > > 1 T @M > > > : 2 /i @yj /i

if iak if k ¼ i

Therefore, the sign of the derivative of an eigenvector term generally does not remain unchanged within the variation of the updating parameters even if the mass and stiffness matrices can be decomposed as in Eqs. (5) and (6). This is because the sign of the aik terms change in the summation in Eq. (14) due to changes in the sign of the denominator term o2i o2k . Therefore the problem of interval model updating cannot be solved by using the parameter vertex solution when eigenvector data is included in the objective. In the following section we consider the solution of the problem in the general case where the outputs behave non-monotonically with respect to the updating parameters. 3.2. Case 2: general case From the first case, the parameter vertex solution is always valid when the mass and stiffness matrices are linear functions of the updating parameters and the output data are the eigenvalues of the dynamic system. Eq. (4) shows that the general case may be obtained by evaluating the inverse of the interval sensitivity matrix at each iteration. This solution is not straightforward, may be impossible and remains as an open problem. However, the use of a meta-model can lead to a solution with very good accuracy depending on the type of meta-model, the sampling used and the behaviour of the outputs within the region of variation. The relationships between the updating parameters and outputs of the meta-model are described by known functions in the region of parameter variations, which are not available when working directly with the FE models. Therefore the use of surrogate models makes the solution of the inverse interval problem, Eq. (4), easier and more efficient as will be demonstrated. The idea for the solution of the interval model updating problem using a meta-model is illustrated in Fig. 1 which shows specifically the procedure for a dynamic system with two updating-parameter inputs and two outputs. In deterministic model updating it is assumed that the measured eigenvalues and eigenvectors of one structure have been obtained from experiments in the form of Eq. (2) while in stochastic model updating, it is assumed that a set of vectors of measured data, Z ¼ ½zð1Þ zð2Þ . . . zðns Þ T , are available (e.g. from a collection of modal test data from a set of nominally identical structures, built in the same way with the same material within manufacturing tolerances). The vector of mean values of measured data can be readily obtained from ns samples as shown in Fig. 1. We assume that the measured eigenvectors are already scaled to the analytical eigenvectors as explained previously. Then the problem of deterministic model updating can be applied to identify the deterministic values of the updating parameters. If the solution of the updating problem is unique, then the vector of updated parameters is represented by a point in the parameter space. An initial hypercube around the updated parameters may be constructed as illustrated in Fig. 1. The meta-model is then used to map the space of the initial hypercube of updating parameters to the space of outputs. If the mapping is good enough to represent the relationship between the input and output data then this model can be used to correct the dimensions of initial hypercube of updating parameters based on the available measured data (circles in the figure). Therefore, selection of the meta-model is a crucial step in that it influences the performance of mapping and consequently the updating procedure to a very significant degree. Since the inverse problem of determining the centre point of the initial hypercube by deterministic model updating from the centre point of measured

Meta-model

:Measured data 2

z2

Initial hypercube

Deterministic FE updating zc2

c 2

Updated hypercube c 1

1

Fig. 1. Interval model updating using Kriging model.

c

z1

z1

H.H. Khodaparast et al. / Mechanical Systems and Signal Processing 25 (2011) 1204–1226

1209

output is generally nonlinear, there is a possibility of multiple solutions. In that case, multiple initial hypercubes corresponding to each solution may be defined and meta-models for each hypercube will be obtained. Amongst existing meta-models such as conventional response surface methods, neural networks and Kriging models, the later are chosen in this paper due to (i) their excellent performance in dealing with non-smooth behaviour between inputs and outputs, which often occurs in dynamic systems with close modes as will be shown later in this paper, and (ii) the high level of degree of accuracy reported in the literature [37,38]. The application of the Kriging predictor for the solution of interval model-updating problems is explained in the following section. The Kriging predictor theory and an optimal sampling method proposed by one of the authors [37] is explained in Appendix A and B, respectively. 4. The Kriging predictor in interval model updating In this section, we apply the Kriging predictor to the problem of interval model updating in structural dynamics. A generalised form of the Kriging predictor for a dynamic system with n output data may be written based on the equation of the Kriging predictor given in Eq. (A.1) as z^ ¼ a þ HðhÞh þ KqðhÞ n

ð15Þ

ns n

rTn T ,

where z^ 2 R , q 2 R ; q ¼ ½rT1 rT2 . . . and ri 2 Rns is introduced in Eq. (A.7),

T

a ¼ ½b0,1 b0,2 . . . b0,n  . b,i are regression coefficients introduced in Eq. (A.1)

p

HðhÞ ¼ ½Hij np ;

Hij ¼ bj,i þ yj bjj,i þ

1X bkj,i yk 2 k ¼ 1 kaj

and K ¼ ½Lij nðns nÞ  lj,i ði1Þns þ 1 rj r i  ns Lij ¼ 0 elsewhere Once the initial Kriging model is constructed, the following error function can be defined for deterministic model updating using the Kriging predictor formed from the measured samples,

e ¼ zm ða þ Hh þ KqÞ ¼ lHhKq

ð16Þ

where l ¼ zm a and the function of h, ðhÞ, is omitted from HðhÞ and qðhÞ for reasons of simplicity. Now the updating of each sample of measured data can be stated as an optimisation problem, minðeT eÞ

ð17Þ

h

It should be noted that the Kriging model has been constructed and validated for the initial hypercube of the updated parameters. Therefore if the solution of the minimisation (17) converges to a point outside the hypercube then a new Kriging model should be constructed by increasing the size of initial hypercube and the procedure repeated. According to Eq. (16), the error function, Eq. (17), can be expanded as

eT e ¼ lT llT HhlT KqhT HT l þ hT HT Hh þ hT HT KqqT KT l þ qT KT Hh þ qT KT Kq A necessary condition for minimising the error function Eq. (18) is that   @ = ðeT eÞ ¼ f0g = ¼ @yj p1

ð18Þ

ð19Þ

Substituting Eq. (18) into Eq. (19) leads to HT lAhfðhÞ þHT Hh þDh þHT Kq þVh þ Uh þ gðhÞ ¼ f0g where ðhÞ is omitted from HðhÞ, DðhÞ, VðhÞ, UðhÞ, qðhÞ and A ¼ ½Aij pp ;

Aij ¼

n X k¼1

DðhÞ ¼ ½Dij pp ;

Dij ¼

VðhÞ ¼ ½Vij pp ;

Vij ¼

mk

@Hkj ; @yi

 p n  @Hlj 1XX @H Hlj lk þHlk yk 2k¼1l¼1 @yi @yi ns n X X

ll,k

k¼1l¼1

UðhÞ ¼ ½Uij pp ;

Uij ¼

@Hkj 1 ¼ Bij,k 2 @yi

ns n X X k¼1l¼1

@Hkj r ðhÞ @yi l,k

ll,k Hkj

@rl,k ðhÞ @yi

ð20Þ

1210

H.H. Khodaparast et al. / Mechanical Systems and Signal Processing 25 (2011) 1204–1226

fðhÞ ¼ ffi ðhÞgp1 ;

fi ðhÞ ¼

ns n X X j¼1k¼1

lk,j mj

@rk,j ðhÞ @yi

gðhÞ ¼ fgi ðhÞgp1 ; gi ðhÞ ¼

  ns X ns n X @r ðhÞ 1X @r ðhÞ lj,l lk,l rk,l ðhÞ j,l þ rj,l ðhÞ k,l @yi 2l¼1k¼1j¼1 @yi

where Bij,k is the component of matrix B introduced in Eq. (A.1). The derivative of the correlation function, given in Eq. (A.4), may be calculated as follows: @rj,i @C ðh, hðjÞ Þ ðjÞ ðjÞ ¼ i ¼ zk,i jyk yk jni 1 signðyk yk ÞCi ðh, hðjÞ Þ @yk @yk In those cases when the function is not differentiable the derivative may be evaluated as Z þ1 Z þ1 f uðxÞjðx,x0 ,s2 Þ dx ¼ lim f ðxÞjuðx,x0 ,s2 Þ dx Eðf uðx0 ÞÞ ¼ lim s-0

1

s-0

ð21Þ

ð22Þ

1

where u denotes d=dx and jðx,x0 ,sÞ is a Gaussian function with parameters x0 and s2. Eq. (22) is obtained using integration by parts and it should be noted that the Gaussian function jðxÞ is zero at 7 1. Eq. (20) can be rearranged for the solution of system parameters h as ðHT H þD þ U þ VAÞh ¼ fðhÞ þ HT lHT KqgðhÞ

ð23Þ

T

Since the matrix (H H + D +U +V  A) is a function of h an iterative procedure is required. However, the solution needs the inverse of matrix (HTH+ D+ U +V A). If this matrix is not invertible an arbitrary weighting matrix can be added to the both sides of Eq. (23) as ðHT H þD þ U þ VAþ WÞh ¼ fðhÞ þHT lHT KqgðhÞ þ Wh

ð24Þ

and following recursive equation is formed for the solution of Eq. (20) T T hl þ 1 ¼ ðHT H þD þ U þ VAþ WÞ1 jh ¼ hl ffðhÞ þ H lH KqgðhÞ þWhgjh ¼ hl

ð25Þ

The iterations continue until convergence on the system parameters h is achieved. The matrix W is chosen so that the matrix (HTH+D + U+ V A+ W) is invertible. The weighting matrix W in Eq. (25) has the effect of regularising the ill-conditioned Eq. (23), equivalent to adding the side constraint [42], ðhl þ 1 hl ÞT Wðhl þ 1 hl Þ

ð26Þ

to the objective function Eq. (17). In Eq. (25), the weighting matrix W may be chosen in the form W ¼ lI. For small values of the regularisation parameter l the original ill-conditioned problem, Eq. (24), remains and when l takes a large value we see from the side constraint (26) that the updated parameters remain unchanged from the previous iteration. An optimal value of l may be obtained from the corner of the L-curve as described by Ahmadian et al. [43]. The procedure for interval model updating can then be defined as follows: 1. Select and update the parameters of the mathematical FE model using the mean vector of measured data. 2. Initialize a hypercube around the updated parameters of the finite element model. 3. Construct a meta-model based on updated mathematical FE model data. This meta-model should describe the relationship between output data and input data within the initial hypercube around the updated parameters accurately. 4. Use the meta-model for updating the initial hypercube by using all sets of measured data. 5. Construct the new hypercube on the region of updated parameters. If the updated hypercube is bigger than the initial hypercube increase the size of initial hypercube and go back to step 3; otherwise go to step 6. 6. Generate output data by using the meta-model to find the region of variation of output data and compare it to the scatter of measured data. 7. End.

5. Numerical case studies: optimal sampling for the Kriging predictor The performance of interval model updating is strongly dependent on the accuracy of the constructed Kriging model. The optimal sampling method explained in Appendix B often provides an accurate Kriging model. In this section the application of the optimal sampling method and fitting, using a Kriging predictor, to a three degree of freedom mass-spring system, shown in Fig. 2, is considered. The construction of Kriging models is demonstrated in two example problems of the three degree of freedom system with close eigenvalues and well-separated modes.

H.H. Khodaparast et al. / Mechanical Systems and Signal Processing 25 (2011) 1204–1226

1211

Fig. 2. Three degree of freedom system.

5.1. Case study 1: close modes, eigenvalue output only Analysing uncertainty in a dynamic system with close or repeated eigenvalues is usually difficult [44,45] because of the rapidly changing input-output relationship. The system shown in Fig. 2, is based on the example used by Friswell et al. [44]. The system parameters are m1 ¼ 1 kg,

m2 ¼ 4 kg,

k1 ¼ k3 ¼ 0

k4 ¼ 2:0 N=m,

m3 ¼ 1 kg, k5 ¼ 2 N=m,

k6 ¼ 1 N=m

ð27Þ

and it is assumed that the stiffness parameter k2 is uncertain within the interval [6.5 9.5]. Now the third eigenvalue of the system is assumed to be an unknown function of the uncertain parameter k2. Figs. 3(a) and (b) show the initial and final selected samples, the target function (solid line) and Kriging approximation (dashed line). The initial and final values of MSE function, given by Eq. (B.1), with the uncertain parameter k2 are shown in Figs. 3 (c) and (d). It can be seen from Fig. 3(d) that the sampling stops because the maximum value of MSE falls below a specified tolerance. Fig. 3(d) also shows that by increasing the number of samples the Kriging approximation produces a more accurate representation of the target function. Fig. 3(b) shows that the Kriging model represents the true function very accurately indeed. This problem is in fact solvable by the vertex method, but it serves here to illustrate the excellent performance of the Kriging method in the case of close modes. In the second example the vertex solution is not applicable since an eigenvector term is deliberately included in the outputs. 5.2. Case study 2: well-separated modes with an eigenvector term in the outputs The procedure is now applied to the three degree of freedom mass-spring system with well-separated modes. The eigenvalues of the system shown in Fig. 2 are well-separated if the parameters are m2 ¼ 1 kg,

m3 ¼ 1 kg,

k1 ¼ k2 ¼ k3 ¼ k4 ¼ k5 ¼ 1 N=m

k6 ¼ 3 N=m

ð28Þ

The mass parameter m1 is assumed to be uncertain and can be changed within the interval [0.6 0.9]. The first component of the first eigenvector (mode shape) ðf1,1 Þ is assumed to be an unknown function. Fig. 4 shows the Kriging model and target function together with the mean square error estimate with initial samples. It can be seen from Fig. 4 that the Kriging model is in very good agreement with the target function and the MSE is sufficiently small. Fig. 4(a) demonstrates that the parameter vertex solution is invalid since f1,1 does not change monotonically with changes in the mass m1. 6. Numerical case studies: interval model updating The three degree of freedom mass-spring system, shown in Fig. 2, with well-separated and close modes is used in this section to illustrate the performance of the interval model updating using the Kriging method. 6.1. Case study 3: well-separated modes It is assumed that the true value of the bounds of the unknown uncertain parameters of the system are given by k1 ¼ ½0:8 1:2 N=m,

k2 ¼ ½0:8 1:2 N=m,

k5 ¼ ½0:8 1:2 N=m

and other parameters are known and given in Case Study 2. The measured data are obtained by using latin hypercube sampling (LHS) with 10 samples. The method is also applied with an unrealistic number of measured data (10 000 samples) to demonstrate the asymptotic properties of the method. Later, different runs of the updating algorithm (with 10 different sets of 10 measured samples) are carried out and a range of solution errors are determined. We mentioned before that the interval model updating approach needs an initial estimate of the ranges of unknown parameters. In this case the initial estimates are chosen to be, k1 ¼ ½0:5 1:5 N=m,

k2 ¼ ½0:5 1:5 N=m,

k5 ¼ ½0:5 1:5 N=m

1212

H.H. Khodaparast et al. / Mechanical Systems and Signal Processing 25 (2011) 1204–1226

4.3

4.3

Target Function

Target Function

4.25

4.25 Observed Values

Observed Values

4.2

Kriging model

ω32 (rad/s)2

ω32 (rad/s)2

4.2 4.15 4.1

4.15 4.1

4.05

4.05

4

4

3.95

Kriging model

3.95 6.5

7

7.5

8 k2 (N/m)

8.5

9

9.5

6.5

x 10−3

7.5

8 k2 (N/m)

7

7.5

8 k2 (N/m)

8.5

9

9.5

x 10−3

7

7

6

6 5 MSE

5 MSE

7

4

4 3

3 2 2

1

1

0

0 6.5

7

7.5

8 k2 (N/m)

8.5

9

9.5

6.5

8.5

9

9.5

Fig. 3. (a and b) Kriging approximations and (c and d) MSE of the third eigenvalue (3 DOFs system with close modes) versus uncertain stiffness parameter k2. (a) Initial sample, (b) final sample, (c) initial MSE and (d) final MSE.

10

−0.5785 Target Function

8

Observed Values Kriging model

6

−0.5795

MSE

ϕ (1,1)

−0.579

x 10−5

4

−0.58 2 −0.5805 0 −0.581 0.6

0.65

0.7

0.75 0.8 m1 (Kg)

0.85

0.9

0.6

0.65

0.7

0.75 0.8 m1 (Kg)

0.85

0.9

Fig. 4. (a) Kriging approximation and (b) MSE of the first component of first mode shape (3 DOFs system with well separated modes) against the uncertain mass parameter m1.

H.H. Khodaparast et al. / Mechanical Systems and Signal Processing 25 (2011) 1204–1226

1213

It is assumed that the mean values of the updating parameters have already been identified and the errors are in the estimation of the parameter bounds. The output data are assumed to be three eigenvalues and the absolute value of the first eigenvector at the first degree of freedom, z ¼ ½z1 z2 z3 z4 T ¼ ½o21 o22 o23 jf1,1 jT The parameter vertex method is not applicable in this case because of including the eigenvector term in the response vector. To construct the Kriging model, an initial sample was taken based on CCD with face centred points [34]. The MSE values show that these initial samples are enough for a Kriging model to map the initial hypercube of input data to the output data. The Kriging model was constructed using a second order polynomial. Results obtained by interval model updating with 10 and 10 000 measured samples are shown in Table 1. The weighting matrix in Eq. (25) was set to W= 0 in this case. Fig. 5 shows the initial, true and updated bounding hypercube of uncertain parameters in the planes k1  k2 and k2  k5 in the presence of 10 measured samples. The updated hypercube of uncertain parameters is in good agreement with the true hypercube as shown in Fig. 5 and Table 1. Fig. 6 shows the convergence of the initial output-data space upon the space of 10 measured samples in the planes of (a) o21 and o22 , (b) o21 and o23 , (c) o22 and o23 , (d) o21 and jf1,1 j, (e) o22 and jf1,1 j and (f) o23 and jf1,1 j . The results in Table 1 show that errors in estimating the bounds of uncertain parameters are significantly reduced even in the case of only 10 measured samples. A very slight difference is seen between the exact space of output data and the updated one, shown in Fig. 6. It might be argued that these results are just for the one particular set of 10 samples. Therefore the interval updating procedure was repeated for 10 different sets of measured samples and it appears that the errors in estimating the bounds of updating parameters range from 0.0% to 4.9%. This shows an important advantage of using interval models rather than probabilistic models [23,25–28] in stochastic model updating. The initial and updated spaces shown in Fig. 6 are obtained by Monte Carlo simulation (MCS) with latin hypercube sampling (LHS) and the Kriging predictor. The true space was achieved by MCS using eigenvalue solutions from the M, K system. 6.2. Case study 4: close modes The quantification of uncertainty in a system with close modes is difficult because of the rapidly changing response surface. The three degree of freedom system, shown in Fig. 2, with close modes is now considered. It is assumed that the true value of the unknown uncertain parameters of the system are given by k2 ¼ ½7:5 8:5 N=m,

k4 ¼ ½1:8 2:2 N=m,

k5 ¼ ½1:8 2:2 N=m

Table 1 Updated results: 3 DOF mass-spring system with well separated modes. Parameters

Initial error %

Updated error % 10 measured samples

Updated error % 10 000 measured samples

k1 k2 k5

[  37.5 25.0] [  37.5 25] [  37.5 25]

[0.4 0.0] [0.8  1.7] [0.8  0.7]

[0.1  0.2] [0.5 0.0] [0.3  0.1]

Bounds of input parameter (10 measured samples) Initial hypercube True hypercube Updated hypercube

1.8 1.6

1.6

1.4

1.4

1.2

1.2

1

1

0.8

0.8

0.6

0.6

0.4

Initial hypercube True hypercube Updated hypercube

1.8

k5

k2

Bounds of input parameter (10 measured samples)

0.4 0.4

0.6

0.8

1

1.2 k1

1.4

1.6

1.8

0.4

0.6

0.8

1

1.2 k2

1.4

1.6

1.8

Fig. 5. Initial, updated and true hypercube of updating parameters based upon 10 measurement samples (system with well-separated modes).

1214

H.H. Khodaparast et al. / Mechanical Systems and Signal Processing 25 (2011) 1204–1226

Bounds of output parameter (10 measured samples) Initial Meta model True model Updated Meta model Measured samples

6 5.5

Initial Meta model True model Updated Meta model Measured samples

9 8.8 8.6

5 ω32 (rad / s)2

ω22 (rad / s)2

Bounds of output parameter (10 measured samples)

4.5 4

8.4 8.2 8

3.5

7.8

3

7.6 7.4 0.6

0.7

0.8

0.9

1

1.1

1.2

1.3

1.4

0.6

0.7

0.8

0.9

ω12 (rad / s)2

1

1.1

1.2

1.3

1.4

ω12(rad/s)2

Bounds of output parameter (10 measured samples)

Bounds of output parameter (10 measured samples) 0.75

Initial Meta model True model Updated Meta model Measured samples

9 8.8

0.7

8.6

0.65

8.4

|φ1,1|

ω32(rad/s)2

Initial Meta model True model Updated Meta model Measured samples

8.2 8

0.6 0.55

7.8 0.5 7.6 7.4 2.5

0.45 3

3.5

4

4.5

ω2 (rad / s) 2

5

5.5

0.6

0.8

0.9

1

1.1

ω1 (rad / s) 2

Bounds of output parameter (10 measured samples)

1.2

1.3

1.4

2

Bounds of output parameter (10 measured samples)

0.75

0.75 Initial Meta model True model Updated Meta model Measured samples

0.7

Initial Meta model True model Updated Meta model Measured samples

0.7

0.65

0.65

0.6

0.6

|φ1,1|

|φ1,1|

0.7

2

0.55 0.5

0.55 0.5

0.45

0.45 2.5

3

3.5

4 ω2 (rad / s) 2

4.5 2

5

5.5

7.2

7.4

7.6

7.8

8

ω3 ( rad / 2

8.2

8.4

8.6

8.8

s)2

Fig. 6. Initial, updated and true spaces of predicted data (100,000 points) based upon 10 measurement samples (system with well-separated modes).

H.H. Khodaparast et al. / Mechanical Systems and Signal Processing 25 (2011) 1204–1226

1215

and other parameters are given in Case Study 1. It is assumed that 10 samples of measured data exist. The initial estimates of the bounds of uncertain parameters are k2 ¼ ½6:5 9:5 N=m,

k4 ¼ ½1:6 2:4 N=m,

k5 ¼ ½1:6 2:4 N=m

The output data is the same as in Case Study 3 (the first three eigenvalues and the absolute value of first component of the first eigenvector, jf1,1 j. Fifteen samples were taken from the space of the initial hypercube of updating parameters according to CCD. The MSE results showed that the initial samples based on CCD were not good enough for mapping the initial hypercube of input data to the output data. Therefore the sampling procedure described in Appendix B was used to improve the Kriging model. The procedure starts with the first output and continues until the maximum MSE value falls below a specified tolerance. Then the procedure continues for the next output and so on until all four outputs are accurately predicted by the Kriging model. Fig. 7 shows how the maximum value of the MSE decays as the sample size is increased. The iteration numbers in Fig. 7 represent the sample size at each step.

3

x 10−3

Maximum MSE

2.5 2 1.5 1 0.5 0 20

40

60

80

100 120 140 160 180 200 Iterations

Fig. 7. Evolution of maximum MSE values for determining the optimal sample number.

Table 2 Updated results: 3 DOF mass-spring system with close modes. Parameters

Initial error %

Updated error % 10 measured samples

k1 k2 k5

[  13.3 11.8] [  11.1 9.1] [  11.1 9.1]

[0.6  0.7] [0.8  1.0] [0.4  0.5]

Bounds of input parameter (10 measured samples)

Bounds of input parameter (10 measured samples)

Initial hypercube True hypercube Updated hypercube

2.6

Initial hypercube True hypercube Updated hypercube

2.6 2.4

2.2

2.2

k4

k5

2.4

2

2

1.8

1.8

1.6

1.6

1.4

1.4 6

6.5

7

7.5

8 k2

8.5

9

9.5

10

1.5

2 k4

2.5

Fig. 8. Initial, updated and true hypercube of updating parameters based upon 10 measurement samples (system with close modes).

1216

H.H. Khodaparast et al. / Mechanical Systems and Signal Processing 25 (2011) 1204–1226

Bounds of output parameter (10 measured samples)

Bounds of output parameter (10 measured samples)

5

5 Initial Meta model True model Updated Meta model Measured samples

4.8 4.6

Initial Meta model True model Updated Meta model Measured samples 4.5 ω32 (rad / s)2

ω22 (rad / s)2

4.4 4.2 4 3.8

4

3.6 3.4 3.2

3.5 0.8

0.85

0.9

0.95

1

1.05

1.1

1.15

1.2

0.8

0.85

0.9

0.95

ω12 (rad / s)2 Bounds of output parameter(10 measured samples)

1.05

1.1

1.15

1.2

Bounds of output parameter (10 measured samples)

5

0.7

Initial Meta model True model Updated Meta model Measured samples

Initial Meta model True model Updated Meta model Measured samples

0.68 0.66 0.64

4.5

0.62

|φ1,1|

ω32 (rad / s)2

1

ω12 (rad / s)2

0.6 0.58

4

0.56 0.54 0.52

3.5

0.5 3

3.5

4 ω2 (rad / s) 2

4.5

0.8

0.9

0.95

1

1.05

ω2 (rad / s) 2

1.1

1.15

2

Bounds of output parameter (10 measured samples)

Bounds of output parameter(10 measured samples) 0.7

0.7

Initial Meta model True model Updated Meta model Measured samples

0.68 0.66

Initial Meta model True model Updated Meta model Measured samples

0.68 0.66

0.64

0.64

|φ1,1|

0.62

|φ1,1|

0.85

2

0.6

0.62 0.6

0.58 0.58

0.56

0.56

0.54

0.54

0.52

0.52

0.5 3.2

3.4

3.6

3.8

4

ω22 (rad / s)2

4.2

4.4

3.6

3.8

4

4.2

4.4

4.6

ω32 (rad / s)2

Fig. 9. Initial, updated and true spaces of predicted data (100,000 points) based upon 10 measurement samples (system with close modes).

1.2

0.07 0.06 0.05 0.04 0.03 0.02 0.01 0

0.07 0.06 0.05 εTε

εTε

H.H. Khodaparast et al. / Mechanical Systems and Signal Processing 25 (2011) 1204–1226

0.04 0.03 0.02

0

20

40 60 Iterations

80

100

0.01

0

50

100 150 Iterations

200 250

Fig. 10. Evolution of error function values Eq. (18) during optimisation using Eq. (25). (a) W=0; (b) W= 10 I.

Case 1

Case 2

Case 3

Case 4

Case 5

Case 6

Case 7

Case 8

Case 9

Fig. 11. Beam locations in the frame structure.

1217

1218

H.H. Khodaparast et al. / Mechanical Systems and Signal Processing 25 (2011) 1204–1226

The Kriging model was constructed using a first order polynomial and results obtained using 10 measured samples are shown in Table 2. These results, also shown in Fig. 8, confirm that the modified bounds of the uncertain parameters have been determined with very good accuracy. Also, Fig. 9 shows that the updated output-data spaces obtained by Kriging are in good agreement with the true output-data spaces. The latter were obtained by direct solution of the eigenvalue problem of the dynamic system. The weighting matrix, introduced in Eq. (25), was set to W= 10 I in this case. Fig. 10 shows the evolution of error function, Eq. (18), in two cases: (i) W= 0 and (ii) W= 10 I. It is seen in Fig. 10, that the procedure fails to converge when (a) W=0, due to ill-conditioning. This problem is overcome by using the technique described in Eq. (25).

Fig. 12. (a) Frame structure. (b) Finite element model.

Table 3 Measured and predicted natural frequencies (free–free frame structure—case 1).

Mode Mode Mode Mode Mode Mode Mode Mode Mode Mode

(1) (2) (3) (4) (5) (6) (7) (8) (9) (10)

Measured (Hz)

FE (Hz)

FE error %

69.3 79.5 93.2 199.1 235.6 259.8 286.3 297.1 299.1 318.6

70.94 80.27 92.07 200.58 236.17 259.33 288.73 296.4 303.03 327.58

2.37 0.97  1.21 0.74 0.24  0.18 0.85  0.24 1.31 2.82

Table 4 Measured and predicted natural frequencies (fixed-frame structure—case 1). Measured (Hz)

FE (Hz)

FE error %

Mode shape

22.54 (f1) 27.84 (f2) 47.63 (f3) 81.19 (f4) 201.35 233.71 256.40 (f5) 257.68 283.09 298.46 312.39 (f6)

22.59 27.27 48.14 80.89 201.55 233.41 259.05 256.54 283.35 305.34 316.49

0.22  2.04 1.08  0.37 0.10  0.13 1.03  0.44 0.09 2.30 1.31

first in-plane bending mode first out-of-plane bending mode first torsion mode second in-plane bending mode higher order in-plane bending mode higher order in-plane bending mode second out-of-plane bending mode higher order in-plane bending mode higher order in-plane bending mode higher order in-plane bending mode second torsion mode

H.H. Khodaparast et al. / Mechanical Systems and Signal Processing 25 (2011) 1204–1226

1219

Table 5 Measured and predicted natural frequencies (fixed-frame structure—case 2). Measured (Hz)

FE (Hz)

FE error %

Mode shape

23.78 (f1) 27.43 (f2) 49.85 (f3) 79.41 (f4) 194.40 222.84 (f5) 226.55 256.09 263.06 289.42 306.56 (f6)

23.97 26.97 50.47 79.65 193.01 227.90 227.32 254.80 264.44 289.28 311.63

0.82  1.65 1.23 0.31  0.71 2.27 0.34  0.50 0.52  0.05 1.65

first in-plane bending mode first out-of-plane bending mode first torsion mode second in-plane bending mode higher order in-plane bending mode second out-of-plane bending mode higher order in-plane bending mode higher order in-plane bending mode higher order in-plane bending mode higher order in-plane bending mode second torsion mode

Table 6 Measured and predicted natural frequencies (fixed-frame structure—case 3). Measured (Hz)

FE (Hz)

FE error %

Mode shape

23.33 (f1) 26.92 (f2) 48.94 (f3) 74.60 (f4) 194.06 219.94 (f5) 232.23 253.54 260.93 299.94 (f6)

23.52 26.40 49.43 74.73 191.88 224.41 231.90 255.83 260.74 304.92

0.80  1.96 1.00 0.17  1.12 2.03  0.14 0.90  0.07 1.66

first in-plane bending mode first out-of-plane bending mode first torsion mode second in-plane bending mode higher order in-plane bending mode second out-of-plane bending mode higher order in-plane bending mode higher order in-plane bending mode higher order in-plane bending mode second torsion mode

Table 7 Measured and predicted natural frequencies (fixed-frame structure—case 4). Measured (Hz)

FE (Hz)

FE error %

Mode shape

24.31 (f1) 24.38 (f2) 47.17 (f3) 76.68 (f4) 198.89 212.54 220.52 (f5) 248.41 257.87 291.62 299.65 304.76 (f6)

24.53 24.25 47.77 76.69 199.23 207.79 225.85 247.78 258.93 292.35 304.99 310.00

0.89  0.55 1.28 0.01 0.17  2.23 2.42  0.25 0.41 0.25 1.78 1.72

first in-plane bending mode first out-of-plane bending mode first torsion mode second in-plane bending mode higher order in-plane bending mode higher order in-plane bending mode second out-of-plane bending mode higher order in-plane bending mode higher order in-plane bending mode higher order in-plane bending mode higher order in-plane bending mode second torsion mode

Table 8 Measured and predicted natural frequencies (fixed-frame structure—case 5). Measured (Hz)

FE (Hz)

FE error %

Mode shape

24.00 (f1) 24.65 (f2) 48.36 (f3) 80.83 (f4) 201.90 206.69 229.51 254.23 (f5) 269.10 281.61 309.67 (f6)

24.25 24.49 48.93 80.93 195.31 207.14 230.13 258.02 271.75 281.08 314.89

1.02  0.66 1.18 0.12  3.26 0.22 0.27 1.49 0.98  0.19 1.69

first in-plane bending mode first out-of-plane bending mode first torsion mode second in-plane bending mode higher order in-plane bending mode higher order in-plane bending mode higher order in-plane bending mode second out-of-plane bending mode higher order in-plane bending mode higher order in-plane bending mode second torsion mode

1220

H.H. Khodaparast et al. / Mechanical Systems and Signal Processing 25 (2011) 1204–1226

Table 9 Measured and predicted natural frequencies (fixed-frame structure—case 6). Measured (Hz)

FE (Hz)

FE error %

Mode shape

24.34 (f1) 24.43 (f2) 47.13 (f3) 76.63 (f4) 199.87 211.57 220.27 (f5) 247.48 256.69 289.20 298.68 304.66 (f6)

24.53 24.25 47.77 76.69 199.23 207.79 225.85 247.78 258.93 292.35 304.99 310.00

0.76  0.74 1.37 0.07  0.32  1.79 2.53 0.12 0.87 1.09 2.11 1.75

first in-plane bending mode first out-of-plane bending mode first torsion mode second in-plane bending mode higher order in-plane bending mode higher order in-plane bending mode second out-of-plane bending mode higher order in-plane bending mode higher order in-plane bending mode higher order in-plane bending mode higher order in-plane bending mode second torsion mode

Table 10 Measured and predicted natural frequencies (fixed-frame structure—case 7). Measured (Hz)

FE (Hz)

FE error %

Mode shape

23.30 (f1) 26.59 (f2) 48.83 (f3) 74.38 (f4) 192.09 219.48 (f5) 232.17 253.90 260.44 299.72 (f6)

23.52 26.40 49.43 74.73 191.88 224.41 231.90 255.83 260.74 304.92

0.96  0.71 1.22 0.46  0.11 2.25  0.12 0.76 0.11 1.73

first in-plane bending mode first out-of-plane bending mode first torsion mode second in-plane bending mode higher order in-plane bending mode second out-of-plane bending mode higher order in-plane bending mode higher order in-plane bending mode higher order in-plane bending mode second torsion mode

Table 11 Measured and predicted natural frequencies (fixed-frame structure—case 8). Measured (Hz)

FE (Hz)

FE error %

Mode shape

23.794 (f1) 27.088 (f2) 49.785 (f3) 79.311 (f4) 193.21 222.013 (f5) 226.327 253.794 262.621 288.139 305.946 (f6)

23.97 26.97 50.47 79.65 193.01 227.90 227.32 254.80 264.44 289.28 311.63

0.76  0.43 1.37 0.43  0.10 2.65 0.44 0.40 0.69 0.40 1.86

first in-plane bending mode first out-of-plane bending mode first torsion mode second in-plane bending mode higher order in-plane bending mode second out-of-plane bending mode higher order in-plane bending mode higher order in-plane bending mode higher order in-plane bending mode higher order in-plane bending mode second torsion mode

Table 12 Measured and predicted natural frequencies (fixed-frame structure—case 9). Measured (Hz)

FE (Hz)

FE error %

Mode shape

22.577 (f1) 27.497 (f2) 47.536 (f3) 81.122 (f4) 200.543 233.52 255.603 (f5) 256.764 280.807 298.403 311.538 (f6)

22.59 27.27 48.14 80.89 201.55 233.41 259.05 256.54 283.35 305.34 316.49

0.06  0.81 1.28  0.28 0.50  0.05 1.35 -0.09 0.91 2.32 1.59

first in-plane bending mode first out-of-plane bending mode first torsion mode second in-plane bending mode higher order in-plane bending mode higher order in-plane bending mode second out-of-plane bending mode higher order in-plane bending mode higher order in-plane bending mode higher order in-plane bending mode second torsion mode

H.H. Khodaparast et al. / Mechanical Systems and Signal Processing 25 (2011) 1204–1226

1221

7. Experimental case study: frame structure with uncertain beams positions A frame structure with two internal beams is designed in which each beam is independently located at three different positions. The design provides for nine different combinations of beam positions as shown in Fig. 11. Detailed finite element models of each of the nine cases were created in MSC-NASTRAN using 8-noded solid elements (CHEXA). The physical structure is shown in Fig. 12(a) and the finite element model in one configuration of the internal beams in Fig. 12(b). The bolted joint connections are modelled using rigid elements over an area three times greater than the cross-section of the bolts. The boundary conditions where the frame is connected to a rigid base are represented by fixing the nodal displacements in the three translational degrees of freedom over an area, the size of a washer, between the frame and the base. Modal tests using an instrumented hammer were carried out for the frame in both free–free conditions and when fixed to the rigid base. The experimental results and finite element predictions for both boundary conditions and nine cases of internal beam locations are shown in Tables 3–12. The closeness of the finite element predictions to the natural frequencies found in modal test shows that the frame structure and the boundary conditions are accurately modelled. To apply the interval model updating method to this problem, the positions (y1 and y2 ) of the two internal beams are assumed to be the unknown updating parameters as indicated in Fig. 13. The parameter vertex solution is not valid because the global mass and stiffness matrices cannot be decomposed as linear functions of updating parameters. In conventional model updating this choice of updating parameters would require re-meshing of the finite element model at each iteration, which is time-consuming and inelegant. An important advantage of Kriging interpolation is that the updating of nodal coordinates is as straightforward as with any other parameter.

0

1

2

3

4

Fig. 13. Parametrisation of internal beam locations in the frame structure.

Table 13 Deterministic model updating of beam locations. True parameters

Initial parameters

Updated parameters

Initial error %

Updated error %

y1

y2

y1

y2

y1

y2

y1

y2

y1

y2

1.0 1.0 1.0 2.0 2.0 2.0 3.0 3.0 3.0

1.0 2.0 3.0 1.0 2.0 3.0 1.0 2.0 3.0

1.6 1.6 1.6 1.6 2.4 2.4 2.4 2.4 2.4

1.6 2.4 2.4 1.6 2.4 2.4 1.6 1.6 2.4

1.04 1.00 1.00 2.04 2.13 1.95 2.98 2.99 2.93

1.02 2.15 3.08 0.90 2.00 3.09 0.89 1.83 2.98

60.00 60.00 60.00  20.00 20.00 20.00  20.00  20.00  20.00

60.00 20.00  20.00 60.00 20.00  20.00 60.00  20.00  20.00

3.73  0.21 0.20 1.81 6.48  2.36  0.58  0.31  2.18

2.00 7.56 2.76  9.78  0.12 3.06  11.00  8.36  0.58

1222

H.H. Khodaparast et al. / Mechanical Systems and Signal Processing 25 (2011) 1204–1226

Bounds of output parameter (9 measured samples)

Bounds of output parameter (9 measured samples) 40

40 Initial Meta model Measured samples

36

36

34

34

32

32

30

30

28

28

26

26

24

24

22

21

21.5

22

22.5

23 23.5 f1 (Hz)

24

24.5

Updated Meta model Measured samples

38

f2 (Hz)

f2 (Hz)

38

22

25

21

Bounds of output parameter (9 measured samples)

22.5

23 23.5 f1 (Hz)

24

24.5

25

88 Initial Meta model Measured samples

86 84

84

82

82

80

80

78

78

76

76

74

74

72

72

70

45

50

55 f3 (Hz)

60

Updated Meta model Measured samples

86

f4 (Hz)

f4 (Hz)

22

Bounds of output parameter (9 measured samples)

88

70

65

Bounds of output parameter (9 measured samples)

45

50

55 f3 (Hz)

60

65

Bounds of output parameter (9 measured samples)

340

340 Initial Meta model Measured samples

335

335

330

330

325

325

320

320

f6 (Hz)

f6 (Hz)

21.5

315

315

310

310

305

305

300

300

295 210

220

230

240 f5 (Hz)

250

260

270

Updated Meta model Measured samples

295 210

220

230

240 f5 (Hz)

250

260

Fig. 14. Initial and updated spaces of predicted data (100,000 points) based upon 9 measurement samples (frame structure).

270

H.H. Khodaparast et al. / Mechanical Systems and Signal Processing 25 (2011) 1204–1226

1223

Table 14 Measured, initial and updated bounds of natural frequencies (frame structure). Measured

Initial FE (Hz)

Updated FE (Hz)

Initial FE (Hz)

Updated FE % error

% error

First in-plane bending mode First out-of-plane bending mode First torsion mode Second in-plane bending mode Second out-of-plane bending mode Second torsion mode

[22.54 24.34] [24.38 27.84] [47.13 49.85] [74.38 81.19] [219.48 256.40] [299.72 312.39]

[21.62 24.61] [23.66 35.53] [43.72 67.57] [71.09 82.50] [224.08 267.34] [300.26 339.65]

[22.57 24.61] [23.86 27.47] [45.13 50.55] [73.99 81.37] [224.08 259.51] [303.58 317.20]

[  4.08 1.11] [  2.95 27.62] [  7.24 35.55] [  4.42 1.61] [2.10 4.27] [0.18 8.73]

[0.13 1.11] [  2.13  1.33] [  4.24 1.40] [  0.52 0.22] [2.10 1.21] [1.29 1.54]

We suppose the initial bounds on y1 and y2 to be [0.5 3.5] and a Kriging model is constructed over this range. The Kriging model describes the relationship between the input parameters (y1 and y2 ) and six outputs; the first and second in-plane and out-of-plane bending natural frequencies and the first and second torsional natural frequencies marked as fi, i= 1,y,6 in Tables 4–12. The maximum value of the MSE shows that CCD design at three levels in the ranges of [  0.5 0.5], [  1.0 1.0] and [1.5 1.5] provides an accurate fit. The FE model is symmetric in the sense that, for example, both beams in position yi ¼ 1 have identical natural frequencies to when the same beams are in positions yi ¼ 3, i= 1,2. The Kriging meta-model was produced without taking account of this symmetry. The optimisation procedure described in Section 4 was used to identify the locations of internal beams based on the six measured natural frequencies, thereby allowing the updating parameter bounds to be corrected. The weighting matrix was set to W = 100I. Table 13 shows the initial and identified beams locations in nine cases obtained by deterministic model updating. The maximum error of 11.00% in Table 13 is an indicator of good performance. The Kriging model was used to generate all possible variations of the six outputs due to the variation of the internal beam locations in the range of [1.00 2.99] for y1 and [0.89 3.09] for y2 by interval model updating. The initial and updated regions of possible natural frequency variation are shown in Fig. 14. In Fig. 14 (a and b) the planes of first and second natural frequencies are shown, (c and d) the planes of third and fourth natural frequencies and (e and f) the planes of fifth and sixth natural frequencies, in each case together with nine measured samples. It is seen from Fig. 14(b), (d) and (f) that the updated regions encloses some measured samples but not all of them. This is due to the fact that the samples which are just outside the regions were in reality located either very close to, or at the boundaries. The errors from other sources of uncertainty, typically disassembly and reassembly and measurement noise, also affect the results causing the samples to move over the boundaries. As can be seen in Fig. 14, some of the areas within the region of output data include a greater number of output samples (they look denser). These are the regions where the likelihood of output data is greatest. The initial and updated bounds of natural frequencies are shown in Table 14 where a maximum error of 4.24% shows good agreement of the updated model output bounds with the bounds of the measured data. The errors are calculated based on the percentage of difference between upper(lower) bounds of measured data and their numerical-prediction counterparts. 8. Conclusion The problem of interval model updating, with test structure variability is formulated. In particular cases, when the output data are the eigenvalues of the dynamic system and updating parameters are substructure mass and stiffness coefficients, the parameter vertex solution can be used. The Kriging predictor for the solution to the inverse problem of a system with n outputs is formulated and used for the solution of interval model updating in the general case. The method is verified numerically in a three degree of freedom mass-spring system with well-separated and close modes. Results show that interval model updating is capable of identifying input parameters with very good accuracy when only a small numbers of measured samples are available. This represents a significant advantage of interval updating over probabilistic methods, which require large amounts of data. It was shown that by Kriging interpolation the uncertain positions of internal beams in a frame structure could be treated as updating parameters. Interval model updating with the Kriging predictor was able to correct initial erroneous bounds on the beam positions with good accuracy.

Acknowledgement This work is funded by the European Union for the Marie Curie Excellence Team ECERTA under contract MEXT-CT-2006 042383. Appendix A. The Kriging predictor We assume that ns vectors of independent input parameters H ¼ ½hð1Þ hð2Þ . . . hðns Þ T with h 2 Rp1 are selected by using a suitable sampling method. The corresponding output parameters Z ¼ ½zð1Þ zð2Þ . . . zðns Þ T with zðiÞ 2 Rn1 are true responses of the system; in this paper the eigenvalues and eigenvectors of the dynamic system. A universal Kriging predictor for the

1224

H.H. Khodaparast et al. / Mechanical Systems and Signal Processing 25 (2011) 1204–1226

output data consists of a second order polynomial and a random function, expressed [36] as z^ i ¼ b0,i þ

p X

p X

bk,i yk þ

k¼1

bkk,i y2k þ

p XX kol l ¼ 2

k¼1

1 2

bkl,i yk yl þ ei ðhÞ ¼ b0,i þ bTi h þ hT Bi h þ ei ðhÞ

ðA:1Þ

where b,i are regression coefficients, bi ¼ ½b1,i b2,i    bp,i Tp1 , 2 6 6 6 6 6 Bi ¼ 6 6 6 6 4

2b11,i

b12,i 2b22,i













3

b1p,i b2p,i 7 7



 

 

 sym:

2bpp,i

7 7 7 7 7 7 7 5 pp

and ei ðhÞ is a random function having zero-mean and covariance covðei ðhÞ, ei ðhðhÞ ÞÞ ¼ s2i Ci ðh, hðhÞ Þ

ðA:2Þ

where s2i is the variance of the ith output data and Ci is the correlation function between untried input parameter h and one of the design sample hðhÞ , h =1:ns. A suitably chosen correlation function may improve the quality of fit as explained below. The random function represents the error and since our application is to fit a regression model on a deterministic finite element computer code, any lack of fit will be due entirely to modelling error (incomplete set of regression terms), not measurement error or noise [46]. Hence, the random function in Eq. (A.1) becomes a function of the system parameters h and the errors of the output predictor in Eq. (A.1) are correlated. The correlation function of the prediction errors is assumed to be related inversely to the distance between the corresponding points in the output [46]. The closer the points in space, the greater the correlation between the error terms. Because the components of input parameters are statistically independent (for example the parameters of a finite element model chosen for updating), one may calculate the correlation function between the input parameters as [36], Ci ðh, hðhÞ Þ ¼

p Y

ðhÞ

Cj,i ðyj , yj Þ

ðA:3Þ

j¼1

Different types of correlation functions have been introduced in [47,48]. The choice of correlation function depends on underlying behaviour of the true response. However, this underlying behaviour is often not readily apparent, in which case the following correlation function may be used Cj,i ðyj , yj Þ ¼ expðzj,i jyj yj jni Þ, ðhÞ

ðhÞ

1 r ni r 2

ðA:4Þ

where zj,i (the jth term of the vector fi ) and ni are parameters of the correlation function at the ith output. ni ¼ 1 gives an Ornstein–Uhlenbeck process which produces continuous paths but not very smooth. The case ni ¼ 2 produces infinitely differentiable paths. Therefore the parameter ni is related to the smoothness of the function in yj coordinates. As it is seen in h Eq. (A.4), the correlation function is 1 when yj ¼ yj and its value reduces as the untried point yj is positioned further away h from the hth design point yj . Since the predictor is unbiased at the observation point, a high level of confidence in the prediction of the outputs for the points which are close to the design samples can be achieved. The level of confidence of the predictions at untried points can be assessed by evaluating the mean squared error MSE, to be discussed in Appendix B. The parameter zj,i controls the importance of the jth component. The calculation of correlation parameters is discussed below. To determine the Kriging model for the ith output, the regression coefficients b,i , in Eq. (A.1) and correlation parameters f and m in Eq. (A.4) must be estimated. The procedure for estimating these unknown coefficients starts with an initial estimates of correlation parameters, f and m . Having known correlation parameters, the regression coefficients b,i and variances s2i of the output data can be estimated using a weighted least-square technique as

b:,i ¼ ðNT Ri NÞ1 NT R1 i Z:,i

s2i ¼

1 ½Z Nb:,i T R1 i ½Z:,i Nb:,i  ns :,i

where Z:,i is the ith column of matrix Z consisting of the observed responses of output vector components,

b:,i ¼ ½b0,i b1,i b2,i . . . bp,i b11,i b22,i . . . bpp,i b12,i . . . bpðp1Þ,i T ,

ðA:5Þ

ðA:6Þ

H.H. Khodaparast et al. / Mechanical Systems and Signal Processing 25 (2011) 1204–1226

2

1 6 6 61 6 6  N¼6 6 6 6 6 6 4 1

ð1Þ





y2p

ð1Þ



y21





y2p

ð2Þ



 

 

 

 

 

















y2p

yð1Þ 1







yð1Þ p

y21

yð2Þ 1







yð2Þ p

 

 

 

 

 











sÞ yðn 1







sÞ yðn p

ð2Þ

ðns Þ

y21

ðns Þ

1225

3

ð1Þ yð1Þ 1 y2







ð1Þ yð1Þ p1 yp

ð2Þ yð2Þ 1 y2







 

 

 

 

ð2Þ 7 yð2Þ p1 yp 7









s Þ ðns Þ yðn 1 y2







7 7

7 7 7 7 7 7 7  5 ðn Þ ðn Þ  

s yp1 yp s

and Ri 2 Rns ns is the correlation matrix between samples with components Rhq,i ¼ Ci ðhðhÞ , hðqÞ Þ. The number of regression coefficients is equal to 1 for a zero-order polynomial model, p+1 for a first order polynomial model and 12 ðp þ1Þðp þ 2Þ for a second order polynomial model. Minimising the mean square error (MSE) then leads to the mean Kriging predictor expressed as T z^ i ¼ b0,i þbi h þ 12hT Bh þ kTi ri ðhÞ

ðA:7Þ

where ri ðhÞ 2 Rns 1 , ri ðhÞ ¼ ½Ci ðh, hð1Þ Þ Ci ðh, hð2Þ Þ    Ci ðh, hðns Þ ÞT and ki ¼ R1 ðZ:,i NbÞ. The predictor is unbiased at the observation points. Further detail on the derivation of Eq. (A.7) can be found in [36,46]. The updated correlation parameters are then estimated by minimising the following objective function. minðdetRi Þ1=ns s2i

ðA:8Þ

fi , mi

Therefore an iterative procedure can be defined for evaluation of regression coefficients b,i and correlation parameters fi and m i as follows: 1. 2. 3. 4. 5.

Estimate initial values of correlation parameters f and m . Evaluate b,i and s2i using the generalized least-square solution given in Eqs. (A.5) and (A.6). Update the correlation parameters to minimise the objective function given by Eq. (A.8). If a minimum is found then go to the step 5. Otherwise go to step 2. End.

Appendix B. Sampling for the Kriging predictor As shown in [37], a systematic method for the generation of samples for Kriging is needed to ensure that the uncertainty in the Kriging predictor is minimised in representing the output. This also has a significant effect on the accuracy of the inverse problem which was discussed in Section 5. To minimise this uncertainty, the Kriging-predicted mean squared error (MSE) can be used as the criterion for sample generation. The Kriging predictor provides an estimation of the MSE [36] as MSEðhÞ ¼ s2i ½1uT Pu

ðB:1Þ 2

2

2

where u ¼ ½1 y1 y2 . . . yp y1 y2 . . . yp y1 y2 . . . yp1 yp rTi ðhÞT and # " #1 " 1 1 T 1 ðNT R1 N Ri ðNT R1 O NT i NÞ i NÞ P¼ ¼ T 1 T 1 1 1 T 1 N Ri R1 R1 N Ri Þ i NðN Ri NÞ i ðINðN Ri NÞ It should be recalled that the predictor is unbiased at the observed points and therefore the MSE is zero at these points. The following procedure may be defined for sampling: 1. 2. 3. 4.

Generate a central composite design (CCD) with one centre point [34] as the initial sample. Generate the next samples at those locations where the MSE is maximum. If the maximum MSE is smaller than a threshold go to the next step 4; otherwise go to step 2. End.

References [1] J.E. Mottershead, M.I. Friswell, Model updating in structural dynamics: a survey, Journal of Sound and Vibration 167 (2) (1993) 347–375. [2] M.I. Friswell, J.E. Mottershead, Finite Element Model Updating in Structural Dynamics, Kluwer Academic Press, Dordrecht, 1995. [3] Y. Govers, M. Boswald, U. Fullekrug, D. Goge, M. Link, Analysis of sources and quantification of uncertainty in experimental modal data, in: Proceedings of International Conference on Noise and Vibration, ISMA2006, Katholieke Universiteit Leuven, Leuven, Belgium, 2006. [4] J.E. Mottershead, C. Mares, S. James, M.I. Friswell, Stochastic model updating: part 2: application to a set of physical structures, Mechanical Systems and Signal Processing 20 (8) (2006) 2171–2185.

1226

H.H. Khodaparast et al. / Mechanical Systems and Signal Processing 25 (2011) 1204–1226

[5] S. Adhikari, M.I. Friswell, Random matrix eigenvalue problems in structural dynamics, International Journal for Numerical Methods in Engineering 69 (3) (2007) 562–591. [6] D. Moens, D. Vandepitte, A survey of non-probabilistic uncertainty treatment in finite element analysis, Computer Methods in Applied Mechanics and Engineering 194 (12–16) (2005) 1527–1555. [7] K. Worden, G. Manson, T.M. Lord, M.I. Friswell, Some observations on uncertainty propagation through a simple nonlinear system, Journal of Sound and Vibration 288 (3) (2005) 601–621. [8] D. Moens, D. Vandepitte, An interval finite element approach for the calculation of envelope frequency response functions, International Journal for Numerical Methods in Engineering 61 (14) (2004) 2480–2507. [9] M. Hanss, The transformation method for the simulation and analysis of systems with uncertain parameters, Fuzzy Sets and Systems 130 (3) (2002) 277–289. [10] J.D. Collins, G.C. Hart, T.K. Hasselman, B. Kennedy, Statistical identification of structures, AIAA Journal 12 (2) (1974) 185–190. [11] M.I. Friswell, The adjustment of structural parameters using a minimum variance estimator, Mechanical Systems and Signal Processing 3 (2) (1989) 143–155. [12] J.L. Beck, L.S. Katafygiotis, Updating models and their uncertainties. i: Bayesian statistical framework, Journal of Engineering Mechanics 124 (4) (1998) 455–461. [13] J.L. Beck, S.K. Au, Bayesian updating of structural models and reliability using Markov chain Monte Carlo simulation, Journal of Engineering Mechanics 128 (4) (2002) 380–391. [14] C. Soize, E. Capiez-Lernout, R. Ohayon, Robust updating of uncertain computational models using experimental modal analysis, AIAA Journal 46 (11) (2008) 2955–2965. [15] T. Haag, J. Herrmann, M. Hanss, Identification procedure for epistemic uncertainties using inverse fuzzy arithmetic, Mechanical Systems and Signal Processing 24 (7) (2010) 2021–2034. [16] S. Adhikari, M.I. Friswell, Distributed parameter model updating using the Karhunen–Loeve expansion, Mechanical Systems and Signal Processing 24 (2) (2010) 326–339. [17] G. Lallement, J. Piranda, Localisation methods for parameter updating of finite element models in elastodynamics, in: Proceedings of the Eighth International Modal Analysis Conference (IMAC VIII), Orlando, Florida, USA, 1990. [18] M. Link, O.F. Santiago, Updating and localising structural errors based on minimisation of equation errors, in: International Conference on Spacecraft Structures and Mechanical Testing (ESA/ESTEC), Noordwijk, Holland, 1991. [19] G.M. Gladwell, H. Ahmadian, Generic element matrices suitable for finite element updating, Mechanical Systems and Signal Processing 9 (6) (1996) 601–614. [20] J.E. Mottershead, M.I. Friswell, G.H.T. Ng, J.A. Brandon, Geometric parameters for finite element model updating of joints and constraints, Mechanical Systems and Signal Processing 10 (2) (1996) 171–182. [21] M.I. Friswell, J.E. Mottershead, H. Ahmadian, Combining subset selection and parameter constraints in model updating, Transactions of the ASME, Journal of Vibration and Acoustics 120 (4) (1998) 854–859. [22] H. Ahmadian, J.E. Mottershead, S. James, M.I. Friswell, C.A. Reece, Modelling and updating of large surface-to-surface joints in the awe-mace structure, Mechanical Systems and Signal Processing 20 (4) (2006) 868–880. [23] J.R. Fonseca, M.I. Friswell, J.E. Mottershead, A.W. Lees, Uncertainty identification by the maximum likelihood method, Journal of Sound and Vibration 288 (3) (2005) 587–599. [24] C. Mares, J.E. Mottershead, M.I. Friswell, Stochastic model updating: part 1: theory and simulated example, Mechanical Systems and Signal Processing 20 (7) (2006) 1674–1695. [25] X.G. Hua, Y.Q. Ni, Z.Q. Chen, J.M. Ko, An improved perturbation method for stochastic finite element model updating, International Journal for Numerical Methods in Engineering 73 (13) (2008) 1845–1864. [26] H.H. Khodaparast, J.E. Mottershead, M.I. Friswell, Perturbation methods for the estimation of parameter variability in stochastic model updating, Mechanical Systems and Signal Processing 22 (8) (2008) 1751–1773. [27] H. H. Khodaparast, J.E. Mottershead, Efficient methods in stochastic model updating, in: Proceedings of International Conference on Noise and Vibration, ISMA2008, Leuven, Belgium, 2008. [28] Y. Govers, M. Link, Stochastic model updating: covariance matrix adjustment from uncertain experimental modal data, Mechanical Systems and Signal Processing 24 (3) (2010) 696–706. [29] H.H. Khodaparast, J.E. Mottershead, K.J. Badcock, Propagation of structural uncertainty to linear aeroelastic stability, Computers and Structures 88 (3–4) (2010) 223–236. [30] S. Marques, K.J. Badcock, H.H. Khodaparast, J.E. Mottershead, Transonic aeroelastic stability predictions under the influence of structural variability, Journal of Aircraft 47 (4) (2010) 1229–1239. [31] R. Moore, Interval Analysis, Prentice Hall, Englewood Cliffs, 1966. [32] W. Dong, H. Shah, Vertex method for computing functions of fuzzy variables, Fuzzy Sets and Systems 24 (1) (1987) 65–78. [33] Z. Qui, X. Wang, M.I. Friswell, Eigenvalue bounds of structures with uncertain-but-bounded parameters, Journal of Sound and Vibration 282 (1–2) (2005) 297–312. [34] R.H. Myers, D.C. Montgomery, Response Surface Methodology; Process and Product Optimisation using Designed Experiments, John Wiley and Sons, New York, USA, 2002. [35] S. Haykin, Neural Networks: A Comprehensive Foundation, first ed., Prentice Hall PTR, Upper Saddle River, NJ, USA, 1994. [36] J. Sacks, W.J. Welch, T.J. Mitchell, H.P. Wynn, Design and analysis of computer experiments, Statistical Science 4 (4) (1989) 409–435. [37] M. Ghoreyshi, K.J. Badcock, M.A. Woodgate, Accelerating the numerical generation of aerodynamic models for flight simulation, Journal of Aircraft 46 (3) (2009) 972–980. [38] M.D. Munck, D. Moens, W. Desmet, D. Vandepitte, An efficient response surface based optimisation method for non-deterministic harmonic and transient dynamic analysis, Computer Modeling in Engineering and Sciences 47 (2) (2009) 119–166. [39] R.J. Allemang, D.L. Brown, A correlation coefficient for modal vector analysis, in: First International Modal Analysis Conference, IMAC, Las Vegas, 1982, pp. 110–116. [40] N.P. Seif, M.A. Hassanein, A.S. Deif, Inverse problem of the interval linear system of equations, Computing 63 (2) (1999) 185–200. [41] R. Fox, M. Kapoor, Rates of change of eigenvalues and eigenvectors, AIAA Journal 6 (12) (1968) 2426–2429. [42] M. Link, Updating of analytical models-procedures and experience, in: Conference on Modern Practice in Stress and Vibration Analysis, Sheffield, England, 1993. [43] H. Ahmadian, J.E. Mottershead, M.I. Friswell, Regularisation methods for finite element model updating, Mechanical Systems and Signal Processing 12 (1998) 47–64. [44] M.I. Friswell, Z.P. Qiu, F. Zhang, The determination of convex sets of uncertainty in structural dynamics, in: Proceedings of First International Conference on Uncertainty in Structural Dynamics, USD2007, The University of Sheffield, UK, 2007. [45] A. Gallina, T. Uhl, A modal meta-modeling for the analysis of structures subjected to input parameter variations, in: Proceedings of the Second International Conference on Uncertainty in Structural Dynamics, Sheffield, UK, 2009. [46] D.R. Jones, M. Schonlau, W.J. Welch, Efficient global optimization of expensive black-box functions, Journal of Global Optimization 13 (4) (1998) 455–492. [47] E.H. Isaaks, R.M. Srivastava, An Introduction to Applied Geostatistics, Oxford University Press, 1989. [48] S.N. Lophaven, H.B. Nielsen, J. Søndergaard, Dace, a matlab kriging toolbox, Technical Report I MM-TR-2002-12, Technical University of Denmark, DK-2800 Kgs. Lyngby – Denmark, 2002.