Avoiding Mode Pairing when Updating Finite Element Models

LiTH-ISY-R-1974

Avoiding Mode Pairing when Updating Finite Element Models

Magnus Andersson

Department of Electrical Engineering Linkoping University, S-581 83 Linkoping, Sweden

www:http://www.control.isy.liu.se

email:magnusa@isy.liu.se

1997-09-01

REGLERTEKNIK

AUTOMATIC CONTROL

LINKÖPING

Technical reports from the Automatic Control group in Linkoping are available as UNIX-compressed Postscript les by anonymous ftp at the address 130.236.20.24 (ftp.control.isy.liu.se).

AVOIDING MODE PAIRING WHEN UPDATING FINITE ELEMENT MODELS

M. Andersson T. Abrahamsson

Dept. of Electrical Engineering, Linkoping University, Sweden, E-mail: magnusa@isy.liu.se

SAAB AB, Linkoping, Sweden, E-mail:Thomas.Abrahamsson@saab.se

Abstract:

Updating nite element models of complex mechanical structures requires some extra considerations. It is stressed that the two most important aspects on up- dating nite element models are parameter estimation properties and computational expenses. A novel mode-pairing free model updating formulation is found to have good parameter estimation properties. The computational expenses are reduced with a semi-xed modal basis, kept xed during several iterations.

Keywords:

Parameter estimation mechanical systems nonlinear programming

model reduction.

1. INTRODUCTION

In aerospace applications it is for several reasons important that the FE model is accurate, i.e., that the model has a dynamical behavior close to the real structure. For instance, the security aspects on aero-planes require an experimental ground- test verication of the construction in a test set- up. The verication includes measuring of the low- frequency dynamical behavior of the structure for several congurations, e.g., for dierent amount of fuel in the wing-tanks. This measuring is very expensive and time-consuming. If an initial model of the structure can be updated/improved from the measurements of a single conguration, the other congurations can be simulated using the improved model. In order for these simulations to be valid, it is assumed that the improved model for the measured conguration is accurate, and that the dierences between the congurations are easily accounted for in the model. This is mostly the case since the dierences between the congurations are often known masses with known locations.

The specic task of improving an initial FE model using measurements is here denoted model updat- ing. An attractive approach is to parameterize the FE model physically and estimate the numerical values of these parameters. This corresponds to the minimization of a cost function

V

( ) 0 de- scribing the deviation between the model output and the measurements,



_{= arg min}

2

 ^V

^{( ) (1)}

where



denotes the open set of feasible values for the parameter vector . The minimization of

V

in (1) is in general a non-linear problem. An iterative procedure is therefore required in order to obtain an estimate

^

_of



_,

i⁺¹= ⁱ+



ⁱ



⁰= ^init (2)

where

^

_{= limi}

!1

i (3)

Hopefully, a consistent estimate of the parameters are achievable,

^

₌



_.

Updating FE models is a task of very high com- plexity. This is basically due to the size of the FE models and the high number of parameters

N

,

but also the high number of actuators

m

and sen- sors

p

yielding a large amount of measurements.

It is therefore important to reduce the computa- tional expenses for model updating. For instance, model reduction is inevitable when dealing with FE models of complex industrial structures.

However, the complexity of the models is not the main obstacle for obtaining a successful updating.

It is experimentally found that model updating formulations in general suers from local minima in the parameter space. In order to alleviate the risk of ending up in a local minima, the convexity properties of the cost function describing the dif- ference between the measurements and the model output must be enhanced.

2. BASIC THEORY

FE modeling of mechanical structures yield a mass matrix^M and a stiness matrix ^Kof dimension

l

^

l

, respectively. The normal modes ^_i and the resonance frequencies

!

i of the undamped FE model represented by ^M



^K are determined via the generalized eigenvalue problem

(^K;

!

_i^2M)^_i= 0



^_Ti^M_i= 1

 i

= 1

::: l

where the normal modes^_iare mass normalized.(4) Notice that the modal matrix



= ^1

::: 

^_l] diagonalizes^M and^K:



^TM



=^I

 

^TK



=

diag

(

!

¹²

::: !

_l²) =

²

where

diag

denotes a diagonal matrix. Two modes(5)

i and ^_j are considered to be near multiple if their corresponding resonance frequencies

!

i and

!

j are close, i.e., if

!

i^

!

j.

A corresponding FE technique for modeling the damping is not available and ad hoc approaches to construct the damping matrix ^L is therefore inevitable. If



also diagonalizes^L,



^TL



⁼

diag

⁽

^

¹

^{::: }

l) =

; ^{ }

i= 2



i

!

i (6) the damping is so-called proportional. In (6), the



i is the damping ratio of the normal mode ^_i. The assumption of proportional damping is found to be valid for exible structures in practice as well as in theory, see (Geradin and Rixen, 1994).

The input/output (I/O) behavior of a mechanical structure may be modeled

M^x



⁽

^t

^{) +Lx}

_

⁽

^t

^{) +Kx(}

^t

^{) =Bu(}

^t

⁾

y(

t

) =^Cx(

t

) (7) where ^xis the vector of the

l

degrees of freedom (dof.),^uis the vector of

m

input forces, and^ythe vector of

p

displacement measurements. In (7),^B (of dimension

l

^

m

) describes the positions of the actuators yielding the input forces and ^C (

p

^

l

)

the positions of the sensors measuring the output.

Equivalently, (7) may be written



q(

t

) +

; _

^q(

t

) +

^2q(

t

) =



^TBu(

t

)

y(

t

) =^C



(

t

) (8) where ^q =



^x is a transformation to the modal coordinates^q.

The frequency response function (FRF) of (7) is dened as the transfer function^G(

s

) evaluated on the imaginary axis,

s

=

j!

,

G(

!

) =^C;K+

j!

^L;

!

^2M;1B (9) where

!

denotes frequency, and where ^G(

!

) de- notes ^G(

j!

) for simplicity. The FRF of (8) may be conveniently written as

G(

!

) =^Xl

i⁼¹

Ci^Ti^B

!

_i^2;

!

²+

j

2



i

!

i

!

⁽¹⁰⁾ showing the interpretation of the triple^f_i



i

!

i^g

of modal parameters in the frequency domain.

An ecient reduction of the system matrices

M



^L



^K to the low-frequency area can be per- formed as

Mr=



Tr^M



r



^Lr=



Tr^L



r



^Kr=



Tr^K



r

(11) using a reduced modal matrix



r= ^_i¹

::: 

^_i^r].

It remains to determine the

i

j's. Traditionally,^_i^j are chosen through a pairing of the normal modes

i of the model and the corresponding estimates of the structure^

^

_xj,

i^{j $}

^

_xj

_{ j}

_{= 1}

_{::: N}

_{xq (12)}

where

N

_xqdenotes the number of estimated modes of the structure. The actual pairing of modes often relies on the correlation between the

^

^_i and the

xj. The pairing of modes also implies a pairing of the corresponding triples of modal parameters,

fi^j

 

i^j

 !

i^jg$f

^

_xj

_{ }

^{^}_xj

_{ !}

_^xj^g

_{ j}

_{= 1}

_{::: N}

_xq

giving the possibility to construct a reduced(13) damping matrix^Lr, and incorporating deviation in the resonance frequencies in the cost function.

The mode-pairing in (12) is a combinatorial opti- mization problem of rather high complexity. For large exible structures with many modes in the model/structure, sub-optimal algorithms must be used in order to decrease the time for the mode- pairing.

The estimate ^G

^

^x _of ^Gx at the frequencies

!

Ri

 i

= 1

::: N

!can be obtained as

G

^kl(

!

Ri) =

y

k(

!

Ri)

u

l(

!

Ri) (14) where ^

G

kl is the entry of row

k

and column

l

of ^G

^

^x_{. In (14),}

_y

_k(

_!

_Ri) denotes the Fourier 2

transform of the output from sensor

k

, and where

u

l(

!

Ri) denotes the Fourier transform of the sinusoidal input force applied with actuator

l

. In (Williams et al., 1985), the multi-variate mode indicator function (MMIF) is presented. It is used to detect (multiple) modes of the structure as well as compute an input force vector^Ui. The input signal ^u(

t

) = ^Uisin(^

!

_xi

t

) is then applied to the structure in order to obtain the estimate^

^

_xi,

^

xi= Im^fGx(^

!

_xi)^Ui^g (15) for

i

= 1

::: N

_xq, see (Andersson, 1997). The^Ui

is the eigenvector



imin



imin= min_j



ij

 i

= 1

::: N

_xq (16) where

;

Ai^;



ij(^Ai+^Bi)^Uij= 0

 j

= 1

::: m

(17) where^Ai= Re^fG

^

^x_(^

_!

_xi)^gT_Re^fG

^

^x_(^

_!

_xi)^g_and

Bi= Im^fG

^

^x_(^

_!

_xi)^gT_Im^fG

^

^x_(^

_!

_xi)^g_.

3. MODEL UPDATING FORMULATION In ground vibration testing, the estimate ^G

^

^x _of

Gxis obtained in the low-frequency area. The es- timated triple of modal parameters^f

^

_xi

_{ }

^{^}_xi

_{ !}

_^xi^g

are also commonly available.

It is experimentally found that the previous de- scribed mode-pairing technique may introduce se- vere local minima in the parameter space of the cost function

V

. In particular, the assumption of a valid mode pairing is critical whenever the dynamical behavior of the model is not close to the structure's, i.e., when ^{i 6}



, or when near multiple modes of the model and/or the structure are present. An objective have therefore been to develop a mode-pairing free technique to update FE models of complex mechanical structures.

It is also found that a combined modal/frequency domain cost function based on both deviations in resonance frequencies and deviations in properties dened by the the FRF has good parameter estimate properties (Andersson, 1997), i.e., an estimate

^

_of



_where

^

^

 _{ V}

₍

^

₎^

_V

₍ ⁰_{) (18)}

where ⁰ is the initial value of the parameter vector .

3.1 Model reduction Dening^

^

_{i as}

^

i=



q

^

^ai

 i

= 1

::: N

_xq (19)

where



q = ^1

::: 

^_N^q] (20)

^

ai= arg min

a2R N

q

k



q^ai^;

^

_xi^k2 ₍₂₁₎

where

N

qdenotes the number of computed normal modes. The matrix of reduction !ris then dened as



r= ^

_i

_{::: }

^

_N^x

q]



^

_i= ^{; 1}

(^

^

_i)^TM

^

_i^1=2

^

_i

(22) where ^

_i is the mass-normalized version of ^

^

_i.

Notice that the mass-normalization of^

^

_{i in (22)}

is based on the assumption that the mass matrix

M of the model is approximating the mass distri- bution of the structure suciently accurate.

The reduced system matrices^M

_rand^K

_rbecome

Mr=



_Tr^M



r



^K

_r=



_Tr^K



r (23) 3.2 Damping model design

In order to obtain a reduced damping matrix

Lr, some assumptions about the damping dis- tribution in the structure must be made. Let



x = ^x1

::: 

^_xNq^x] denote the "true" modes of the structure, and let ^Lx denote the "true"

damping matrix of the structure. If the damp- ing distribution in the structure are proportional, then



Tx^Lx



xwill be diagonal. Assuming that^

^

_xi

is parallel with^_xi, and that^

_iis parallel with^

^

_xi:

^

xi ^{k }xi



and^

^

_ri^{k }

^

_xi

_{ i}

_{= 1}

_{::: N}

_{xq (24)}

then also



_Tr^Lx



r will be diagonal. These as- sumptions indicate that ^L

_rmay be modeled as

Lr=

diag

(

::: 

2^



_xi

!

i

:::

) (25) where

!

_i is the square root of the Rayleigh quo- tient

R

in

R

(^

_{i) =} ^

_Ti^K

_i

Ti^M

_i ^{= (}

^!

ⁱ⁾²

^{ i}

^{= 1}

^{::: N}

^{xq (26)}

see (Craig, Jr, 1981).

However, the matrices^M

_rand^K

_rare in general not diagonal, i.e., the damping is not proportional since^L

_ris diagonal. If a proportional damping is required during the model updating, a diagonal- ization of^M

_{r and}^K

_r can rst be made as

Mdr=



^TM

_r

 ^

^K

_dr=



^TK

_r



⁽²⁷⁾

where



= ^1

::: 

^_Nq^x] (28) In (28), the basis vectors^_i is dened by

(^K

_r^;

_$

²_i^M

_r)^_{i= 0}

_

^_Ti^M

_r^_{i= 1 (29)}

for

i

= 1

::: N

_xq, compare with (4)-(5). The proportional damping matrix^L

_dr may be approx- imated by the diagonal of



^TL

_r



^:

Ldr=

diag

^;



^TL

_r



^=

diag

(

::: 

_i

:::

) (30) Notice that the damping of the updated model becomes proportional if the parameter estimation is consistent (

^

₌



), and if the mode estimates are unbiased, i.e., if ^

^

_xi ^{k }_xi. This is true even if non-proportional damping is used during the model updating.

The reduced model (for proportional damping) becomes

Mdr^q



^+L

_dr^q

_

^+K

_dr^q₌₍



_dr)^TBu

y

=^C



dr^q (31) where



_dr=



r



⁽³²⁾

The FRF corresponding to (31) becomes

Gdr(

!

) =^C



dr^;Z

_dr(

_!

₎^;1₍



dr)^TB (33) where

Zdr(

!

) =

diag

^;

^{::: $}

²i ^;

!

²+

j

i

!:::

^ (34) 3.3 Cost function

Using all the asssumed measurements avilable, a mode-pairing free cost function can be obtained as

V

( ) =

²¹

V

¹( ) +

²²

V

²( ) +

³²

V

³( ) (35) where

iare cost function weights, and where

V

¹( ) =^XN!

i⁼¹



^G

_dr(

_!

_xRi)^;G

^

^x₍

_!

_xRi)^2_{F (36)}

V

²( ) =^N

x

q

X

i⁼¹







!

i( )^;

!

^_xi

!

^_xi





 2

2

(37)

V

³( ) =^N

x

q

X

i⁼¹



^G

_dr(

 !

^_xi)^U_xi^;Y_xi^22 (38) where

F

denotes the Frobenius norm, and

!

_xRide- notes the measured frequency points of^G

^

^x_{. Mini-}

mizing

V

¹corresponds to the maximum-likelihood (ML) estimate of



, see (McKelvey, 1995). The

V

² is included in order to increase the convexity properties of the cost function, see (Abrahamsson et al., 1996) where the enhanced parameter esti- mation properties using the Rayleigh quotient in the cost function are studied. In (38), ^Y_xi corre- sponds to the output from the structure for the input^u(

t

) =^U_xisin(^

!

_xi

t

),

Yxi=^Gx(^

!

_xi)^U_xi (39)

3.4 Implementation

Using numerical dierence approximations of

@

^M

=@

and

@

^K

=@

, approximating derivatives of^M

_r

_

^K

_r are obtained as

@

^M

_r

@

k ^



_Tr^M(

^k+

^ _

^k_k^);M(

^k)



r (40)

@

^K

_r

@

k ^



Tr^K(

k+



k)^;K(

k)



k



r (41)

If the modal basis of reduction



r and the lin- earization of^M

_{r and}^K

_r,

Mlinr ( +



)^M

_{r( ) +}^X

k

@

^M

_{r( )}

@

k #

_k

(42)

Klinr ( +



)^K

_{r( ) +}^X

k

@

^K

_{r( )}

@

k #

_k are semi-xed, i.e., kept xed for several itera-(43) tions, the computational expenses can be reduced, see (Andersson, 1997).

In Fig. 1, a ow-chart of a general model updating strategy is shown.

Initial Model Measurements Parameterization

Compute^r

New Iteration Update Model New Model

New^r? Yes

Linearize^{Mr Kr} No

Fig. 1. A ow chart of the general model updating procedure.

Let

V

( ) =^v( )^Tv( ) where^v( ) =^vm( )^;vs

is the deviation vector between the output of the model^vmand the measurements of the structure

vs. Linearizing^v in ,

v( +^ )^v( ) +^J( )^ (44) the parameter dierence vector



ⁱin (2) can be determined as



^{i= arg min}



n





J( ⁱ)



^{;d( i)22 (45)}

+



^2



^{22o (46)}

where ^d( ⁱ) = ^vm( ⁱ) ^{; v}s, and



denotes the regularization term, see (More, 1983). In 4

(Andersson, 1997), a robust strategy for determin- ing the value of



is proposed. It is a modied version of Fletcher's strategy in (Wolfe, 1978), being more robust against ending up in undesired local minima.

4. NUMERICAL SIMULATION EXAMPLE Fig. 2 shows a simple aircraft-like structure.

13 26 15

8

12

9 10

11 24

1 3 7

6 5

4 ¹⁴

25

16 17 18

19 23 27 28

20 21 22

2

z y

x 29

Fig. 2. The Garteur model.

In the examples below, 30 sensors are used to collect measurements. The cost function weights

i are determined so that

V

¹= 1014

V V

²_{= 114}

V V

³_{= 314}

V

(47) in (36)-(38) for the initial values of the parame- ters, = ^init.

4.1 Example 1

If the actuators are located to^f12z



14x



15z



20y^g, the 9 low-frequency modes are estimated appro- priately, since

min_i ^j(^_xi)^T

^

_xi^j2

jxi^j2j

^

_xi^{j2 = 0}

^:

^{99 (48)}

The parameter history then become as in Fig. 3. It is seen that the estimated parameter values

^

_have

converged to the global minima in 13 iterations.

The cost function history is shown in Fig. 4. It is seen that the linearization ^M

^lin_r

_

^K

^lin_{r is valid}

since the circles (^) is close to the solid line. The

denotes the values of the cost function for a new computation of the modal basis of reduction



r. 4.2 Example 2

If the actuators are located to ^f12z



15z



20y^g, the 9 low-frequency modes are not estimated appropriately, since

0 2 4 6 8 10 12 14

0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5

iterationⁱ



i k

=

ref k

Fig. 3. The history of the relative parameter values, Example 1.

0 2 4 6 8 10 12 14

-6 -5 -4 -3 -2 -1 0

iterationⁱ

10log(V)

Fig. 4. The history of the cost function, Example 1.

min_i ^j(^_xi)^T

^

_xi^j2

jxi^j2j

^

_xi^{j2 = 0}

^:

^{87 (49)}

The parameter history becomes as in Fig. 5. It can be stated that not appropriately estimated modes may turn up as bias in the parameter estimate.

0 2 4 6 8 10 12

0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5

iterationⁱ



i k

=

ref k

Fig. 5. The history of the relative parameter values, Example 2.

However, even if the mode estimates ^

^

_{xi are not}

parallel with the true modes of the structure ^_xi, the cost function decrease is substantial:

V

(

^

_{) = 6}

_:

₇₁₀^;5

_V

₍ ^init_{) (50)}

The relative large deviation in the estimated pa- rameter values

^

_from



is in view of (50) due to that the parameters are dependent, i.e., the deviation in a parameter value is compensated by the deviations in the other parameter values.

5. CONCLUSIONS

It was pointed out that the two most important aspects on a model updating formulation is its parameter estimation properties and its compu- tational complexity. It was also noted that model updating in general suers from undesired local minima in the cost function. In order to reduce the risk of ending up in a local minima, the tra- ditional technique of pairing modes was circum- vented when constructing the reduced size model.

A novel combined frequency domain/modal do- main mode-pairing free cost function was in simu- lation examples shown to yield good parameter es- timates if the modes of structure are appropriately estimated. If the modes of the structure are not appropriately estimated, bias in the parameter estimates emerged. The value of the cost function are reduced considerably, though. A modal reduc- tion basis was used in several iterations in order to reduce the time-consuming computation of the normal modes.

6. REFERENCES

Abrahamsson, T., M. Andersson and T. McK- elvey (1996). A nite element model up- dating formulation using frequency responses and eigenfrequencies. In: Second Interna- tional Conference on Structural Dynamics Modelling - Test, Analysis and Correlation.

DTA/NAFEMS. pp. 293{305.

Andersson, M. (1997). Experimental design and updating of nite element models. Lic. thesis.

Automatic Control, Dept. of EE, Linkoping University, 581 83 Linkoping, Sweden.

Geradin, M. and D. Rixen (1994). Mechanical Vi- brations - Theory and Application to Struc- tural Dynamics. Wiley.

Craig, Jr, R. R. (1981). Structural Dynamics. An Introduction to Computer Methods. Wiley.

McKelvey, T. (1995). Identication of State-Space Models from Time and Frequency Data. PhD thesis. Linkoping University. Dept. of Electr.

Eng., Linkoping University, 581 83 Linkoping.

More, J. J. (1983). Recent developments inalgo- rithms and software for trust region methods.

In: Mathematical Programming. The State of the Art(A. Bachem, M. Grotchel and B. Ko- rte, Eds.). Springer-Verlag. pp. 258{287.

Williams, R., J. Crowley and H. Vold (1985).

The multivariate mode indicator function in modal analysis. In: Proceedings of the 3rd International Modal Analysis Conference.

pp. 66{70.

Wolfe, M. A. (1978). Numerical Methods for Unconstrained Optimization-an introduction.

Van Nostrand Reinhold.

6