Frequency domain analysis of model order reduction techniques

Finite Elements in Analysis and Design 42 (2006) 367 – 403 www.elsevier.com/locate/finel

Frequency domain analysis of model order reduction techniques Yusuf Cunedio˘glua,∗ , Ata Mu˘gana , Hüseyin Akçayb a Department of Mechanical Engineering, Istanbul Technical University, Gümü¸suyu 80191, Taksim, Istanbul, Turkey b Department of Electrical and Electronics Engineering, Anadolu University, 26470 Eski¸sehir, Turkey

Received 24 November 2004; received in revised form 18 June 2005; accepted 20 August 2005 Available online 10 October 2005

Abstract In this study, some popular model order reduction and superelement techniques are studied in frequency and time domains. Frequency domain identification (FDI) methods are also applied to semidiscrete finite element equations to obtain reduced order models for structural systems. An FDI algorithm that determines the coefficients of the reduced order model by using nonlinear least squares (NLS) method and subspace based identification (SBI) method are applied to a sample problem and the results are compared with component mode synthesis (CMS) and quasi-static mode synthesis (QSM) methods. In literature, phase errors of model order reduction techniques have not been studied; however, it is shown in this paper that phase errors are also important in evaluating the performance of model order reduction techniques. In general, the NLS and SBI methods have better performance than the CMS and QSM methods; however, as the size of problems increases, the NLS method may have convergence problems and the SBI method may yield estimated models having large orders. 䉷 2005 Elsevier B.V. All rights reserved. Keywords: Finite element methods; Model order reduction techniques; System identification methods; Structural analysis

1. Introduction In the absence of analytical solutions, dynamic responses of structural systems are calculated by computational methods. To this end, finite element methods (FEMs) have been widely used. In FEM

∗ Corresponding author. Tel.: 90 212 293 1300x2723; fax: +90 212 249 5197.

E-mail address: [email protected] (Y. Cunedio˘glu). 0168-874X/$ - see front matter 䉷 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.finel.2005.08.005

368

Y. Cunedio˘glu et al. / Finite Elements in Analysis and Design 42 (2006) 367 – 403

models, a structural system is separated into elements and these elements are assumed to be connected with each other at nodes, which require the usage of large numbers of nodes and elements for large domains. Therefore, model order reduction and superelement techniques are frequently used in computational studies to reduce the problem sizes, in particular for computing dynamic responses. Model order reduction methods have been the subject of significant number of recent studies, e.g., see [1] for a detailed overview and discussion. Many techniques have been proposed to obtain reduced order finite element models by reducing the orders of mass and stiffness matrices. Model order reduction techniques can be classified into four groups. One class of techniques such as [2–10] uses a direct reduction of physical model coordinates by omitting some selected degrees-of-freedom (DOF), whose accuracy is very sensitive to the selection of active DOFs. The second class of methods [11,12] is the modal methods having better accuracy, yielding decoupled equations and giving a better insight into the problems; however, they have the following disadvantages: computation of the modes may be costly if too many modes are needed in the analyses and their results may be unacceptable due to the effects of truncated modes. In order to overcome this drawback, quasi-static compensation (QSC) method studied in [13–18] was developed. The third class of methods is the Ritz vector methods (e.g., see [1,19–28,10]) developed to eliminate the need to compute costly eigenvalue problems and to improve the accuracy in cases where eigenvectors are not the best choice; these methods can be regarded as a generalized approach replacing the eigenvectors by more generally defined Ritz vectors. Gu et al. in [29] used the data recovery techniques in dynamic analyses to recover the physical response from the modal response obtained by using the mode-displacement method and they developed the technique upon a quasi-static compensation technique that is applicable for any banded frequency of interest. The forth class of methods studied extensively in this paper is component mode synthesis (CMS) methods [1,16,17], that combine the first three classes of methods. Based on CMS methods, quasi-static mode synthesis (QSM) methods are developed having improved performance in [16–19,1,29]. In all of the above studies, performance of model order reduction techniques are evaluated either in time domain or by using frequency domain magnitude responses, and no attention have been paid to the frequency domain phase errors of these methods. However, as shown in this work, frequency domain magnitude plots that is commonly used as a performance measure may be misleading and frequency domain phase errors can yield significant errors in time domain responses. Motivated by these facts, some popular model order reduction and superelement techniques are studied in frequency domain and two frequency domain identification (FDI) techniques are also employed as model order reduction tools. To this end, the following works can be cited for frequency domain analyses and identification methods. The FDI methods developed in recent years have proved themselves to be efficient in system identification problems in frequency domain, e.g., see [30–32]. Bayard in [33] developed a computational approach to multivariable frequency domain curve fitting problem, based on the two norm minimization; the algorithm used in [33] includes a sparse matrix method for reducing computation and memory requirements on large problems. Verboven et al. [34] presented a computational technique for the frequency-domain identification of multivariable, discrete-time transfer function models based on a cost function minimization, where sparse matrix methods and QR-projections are used for the reduction of computation time and memory requirement and transfer function model estimation is done by using the least squares algorithm to determine the coefficients of the reduced order model. McKelvey et al. in [35] used a realization algorithm followed by a prediction error method by minimizing a least squares type cost function and identified the vibration behavior of an aircraft.

Y. Cunedio˘glu et al. / Finite Elements in Analysis and Design 42 (2006) 367 – 403

369

In this work, by using a sample problem, the results of the FDI methods are compared with the component mode synthesis (CMS) and quasi-static mode synthesis (QSM) methods. The two FDI methods considered in this study are as follows: a nonlinear least squares (NLS) method in frequency domain and subspace based identification (SBI) algorithm. In general, the FDI methods are shown to be superior to CMS and QSM methods at the expense of CPU time for structures having low DOFs.

2. Model order reduction techniques The governing equations of an undamped structural component can be written in the following form: M u¨ + Ku = f ,

(1)

where M and K are, respectively, neq × neq mass and stiffness matrices, u¨ and u are, respectively, acceleration and displacement vectors of length neq , f is a time dependent force vector of length neq , and neq is the number of equations. Following this, the characteristic equation (1) is stated as follows: (K −
2i M)i = 0,

(2)

where
2i is the ith eigenvalue and i is the corresponding eigenvector of the system. Model order reduction techniques are employed to reduce the order of model in (1) to increase computational efficiency. The following four model order reduction techniques are studied in this paper: component mode synthesis method, quasi-static mode synthesis method having one and two centering frequencies, NLS and SBI methods.

2.1. Component mode synthesis method As summarized in [1], CMS method is a model order reduction technique to obtain superelements used in dynamic analyses of structures. CMS method is employed for constructing models to analyze the dynamics of large and complex structures that are generally described by separate component (substructure) models. The behavior of each component of the structure are approximated by a set of basis vectors (Ritz vectors) and then component approximations are assembled to obtain a global approximation for the large and complex structure. Hence, CMS reduces the computational size and cost of the problem. The reduction procedure divides the original system DOF into two sets; namely, (1) may be partitioned as follows:
 
    Maa Mao u¨ a Kaa Kao ua f + = a . (3) Moa Moo u¨ o Koa Koo uo fo For each component of the structure, the displacement vector is partitioned into the vector of active DOF ua and vector of omitted DOF uo , where a+o is equal to total number of unconstrained DOF neq .

Y. Cunedio˘glu et al. / Finite Elements in Analysis and Design 42 (2006) 367 – 403

370

Each component of the structure is associated with the following eigenvalue problem: 
     
o 0 Koo Koa 2 Moo Moa −
= , 0 Kao Kaa Mao Maa a

(4)

where (
2 , ) denotes the eigenvalue–eigenvector pair. Let n denote the set of n retained component normal modes from (4). The aim of the CMS method is to reduce the problem size; that is, n>neq should hold. Therefore, a suitable set of basis coordinates are needed to replace the large set of physical coordinates. The coordinate transformation between the component physical coordinates u and reduced component basis coordinates z is defined by u = T z,

(5)

where T represents the transformation matrix and usually consists of some type of pre-selected component normal modes n . The other component Ritz basis vectors are , named “computational” modes such as the constraint modes, attachment modes, rigid-body modes, inertia-relief modes and quasi-static modes [1]; that is T = [  n ].

(6)

2.1.1. Computational modes As described in [1], the classical constraint modes (CM) are obtained from expressing the uo partition of the static response vector u in terms of ua by solving the following static reduction equation:    
uo 0 Koo Koa = . (7) fa Kao Kaa ua By using (7), the constraint mode matrix can be written as follows [1]:

 −1 K oa −Koo oa CM = = , aa Iaa

(8)

where Iaa is the identity matrix of size a. 2.2. Quasi-static mode (QSM) synthesis method The QSM is particularly studied in this work due to its superior performance over other methods, e.g., see [1] for comparisons. Similar to the derivation of constraint modes, the quasi-static modes are derived by expressing the uo partition of the quasi-static response vector u in terms of ua by solving the upper part of the following quasi-static reduction equation: 
     
Moo Moa uo 0 Koo Koa −
2c = , (9) fa Kao Kaa Mao Maa ua

Y. Cunedio˘glu et al. / Finite Elements in Analysis and Design 42 (2006) 367 – 403

371

where
c is the centering frequency that is usually selected to be in the middle of frequency range of interest for the reduced order model and  −1   uo = − Koo −
2c Moo Koa −
2c Moa ua

(10)

By using (10), the quasi-static mode matrix can be written as follows [1]


oa QSM = aa




 −1   Koa −
2c Moa − Koo −
2c Moo = . Iaa

(11)

By comparing Eqs. (8) and (11), one can observe that the constraint mode matrix is a special case of the quasi-static mode matrix by setting the frequency tuning parameter (i.e., the centering frequency)
c = 0 [1]; because the constraint mode is a special case of the quasi-static mode. In other words, the CMS method can be regarded as a special case of the QSM method. It is noteworthy that the QSM approach is able to include the inertial effects of the omitted DOF. 2.2.1. QSM procedure Component normal modes consist of the first n normal modes given by




on n = . an

(12)

Then, the QSM and normal modes result in the following transformation [1]: 

uo ua

u=




=

oa aa

on an



qa qn




=

oa

Iaa

on an



 qa , qn

(13)

where qa and qn are vectors of generalized coordinates associated with the QSM and the normal modes, respectively. To simplify coupling of the components and keep the physical meaning of active coordinates, (13) is transformed into 

uo u= ua




=

oa

on − oa an

Iaa

0an



 ua . qn

(14)

Therefore, the transformation matrix T and component basis coordinates z are given by [1]
T =

oa

on − oa an

Iaa

0an




=

oa

on

Iaa

0an

 .

(15)

Note that CMS and QSM methods are applied to a sample problem in Section 6.  z=

 ua . qn

(16)

Y. Cunedio˘glu et al. / Finite Elements in Analysis and Design 42 (2006) 367 – 403

372

It is also derived in [1] that the QSM method having more than one set of quasi-static modes by including multiple centering frequencies, where each set has a particular choice of centering frequency. For example, if the quasi-static modes are calculated at two centering frequencies
c1 and
c2 to preserve the physical meaning of the active DOF and facilitate coupling, the additional set of quasi-static modes is subtracted from the first one, after that the transformation matrix can be expressed by [1]

 oa (
c1 ) oa (
c2 ) − oa (
c1 ) on oa (
c1 ) oa (
c2 ) on T = = (17) Iaa Iaa 0aa 0an 0aa 0an and the component basis coordinates are given by z=

ua va qn

,

(18)

where va is the generalized coordinates corresponding to the additional quasi-static modes and oa (
ci ) = −(Koo −
2ci Moo )−1 (Koa −
2ci Moa ) on = on − oa (
c1 )an .

i = 1, 2,

(19) (20)

Note that QSM1 and QSM2 denote QSM method with one and two centering frequencies, respectively. Finally, the reduced stiffness and mass matrices can be determined by using K reduced = T T KT ,

(21)

M reduced = T T MT ,

(22)

where the transformation matrix T may be defined by (6), (15) or (17), and the superposed “reduced” denotes a component matrix associated with the reduced order model. After this step, the equation of motion for the reduced system can be written as follows: M reduced u¨ + K reduced u = F .

(23)

3. Frequency response of a system Methods for frequency response analyses are widely used in the design of structural systems [36]. The frequency response characteristics of a system are typically represented by two curves such as magnitude versus frequency and phase versus frequency, that are called Bode plots of a system. There are some advantages of working with frequency response of a structure in terms of Bode plots. For instance, the behavior of the structure in a much wider frequency range can be displayed in one plot. Suppose that the equation of a structural system having one DOF is given by mu¨ + cu˙ + ku = F .

(24)

Y. Cunedio˘glu et al. / Finite Elements in Analysis and Design 42 (2006) 367 – 403

373

Without loss of generality, assume that initial conditions u(0) and u˙ (0) are zero. Then, the Laplace transform of (24) is given by ms 2 u(s) + csu(s) + ku(s) = F (s),

(25)

where s is the complex Laplace transform variable [36]. Solving for the transfer function yields that G(s) =

1 1/m u(s) = = 2 . 2 F (s) ms + cs + k s + (c/m)s + (k/m)

(26)

Following this, (26) can be simplified by applying the following definitions: (1) (2) (3) (4)


2n = k/m, where
n is the undamped natural frequency,

√ ccr = 2 km, where ccr is the “critical” damping value,  is the amount of proportional damping, 2
n is the multiplier of the velocity term u. ˙ Then, the transfer function in (26) can be stated as follows: G(s) =

1/m u(s) = . F (s) s 2 + 2
n s +
2n

(27)

The frequency response of (27) can be found by substituting s=j
, where
is the excitation frequency and j is the imaginary unit. In literature, when model order reduction techniques are studied in frequency domain, phase error characteristics of these methods have not been studied, which is no less important than the magnitude errors. To this end, an example is presented below. Example. Let G1 and G2 denote the transfer functions of two second order systems as follows: G1 (s) =

s+3 s 2 + 2
n s +
2n

and G2 (s) =

s−3 . s 2 + 2
n s +
2n

(28)

Note that G2 (s) has a zero in the complex right-half-plane; hence, it is called “nonminimum phase system” [36]. By choosing the parameters  = 0.2 and
n = 5 rad/s arbitrarily, the Bode plots of these transfer functions are given in Figs. 1 and 2, where amplitudes and phase angles in the figures are given in dB and degree, respectively. It can be observed in these figures that the Bode gain plots of these two transfer functions are identical; however, the Bode phase plots are different. In order to see the effect of the Bode phase plots on time domain responses of these transfer functions, step responses for these systems are shown in Fig. 3. Observe that time domain responses of these systems are completely different whereas the Bode gain plots of them are identical. In fact, the gain and phase of the frequency response of a system cannot be independent as presented below.

Y. Cunedio˘glu et al. / Finite Elements in Analysis and Design 42 (2006) 367 – 403

-4 -6 -8

Amplitude

-10 -12 -14 -16 -18 -20 -1 10

0

10

1

10

ω (rad/sec) 40

20

Phase angle in degrees

374

0

-20

-40

-60

-80

-100 -1 10

0

10

ω (rad/sec)

Fig. 1. Bode plots of the transfer function G1 (s).

1

10

Y. Cunedio˘glu et al. / Finite Elements in Analysis and Design 42 (2006) 367 – 403

375

-4 -6 -8

Amplitude

-10 -12 -14 -16 -18 -20 -1 10

0

10

1

10

ω (rad/sec) 200

Phase angle in degrees

150

100

50

0

-50

-100 -1 10

0

10

1

10

ω (rad/sec)

Fig. 2. Bode plots of the transfer function G2 (s).

3.1. Bode gain-phase theorem According to Bode gain-phase theorem, design specifications of a system can be stated in terms of the gain and phase of the frequency response of open loop transfer function. On the other hand, the

376

Y. Cunedio˘glu et al. / Finite Elements in Analysis and Design 42 (2006) 367 – 403 0.3 G1 (s) G2 (s) 0.2

Amplitude

0.1

0

-0.1

-0.2

-0.3

-0.4

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

Time (sec)

Fig. 3. Step responses of the transfer functions G1 (s) and G2 (s).

gain and phase of a frequency response are not mutually independent. For design purposes, the value of one is completely determined once the other is specified. This relationship is stated precisely in [37] as expressed below. L(s) is a proper rational function with real coefficients and with no poles or zeros in the closed right half plane. Then, at each frequency s = j
0 , the following integral relation must hold:
  |v| 1 ∞ d log |L| < L(j
0 ) − < L(0) = log coth dv, (29)  −∞ dv 2 where v = log(
/
0 ). Theorem (29) states the condition under the knowledge of open loop gain along the j
axis that suffices to determine open loop phase to within ± . Hence, for stable systems with no right half plane zeros, gain and phase cannot be manipulated independently in design. Subsequently, it is focused on phase error characteristics of model order reduction techniques to better understand the behavior of these methods in time domain. Two frequency domain identification techniques presented below are also applied to structural systems as a model order reduction tool and compared with conventional model order reduction methods such as CMS and QSM methods.

4. Least squares method in frequency domain In this section, following [31,32], an FDI algorithm based on the nonlinear least squares (NLS) method is used to obtain the reduced order models in frequency domain. Considering (1), let G(s) denote the

Y. Cunedio˘glu et al. / Finite Elements in Analysis and Design 42 (2006) 367 – 403

377

transfer function between an input (i.e., the force) and an output (i.e., the displacement). The continuoustime transfer function of a linear time-invariant system is a rational function in terms of s as follows: G(s) =

num(s) , den(s)

(30)

where the numerator and denominator are polynomials of order nord and dord , respectively. On the other hand, the Laplace transform of a time varying signal u(t) is defined by  ∞ U (s) = {u(t)} = u(t) e−st dt. (31) 0

Then, G(s) defines a relationship between the Laplace transform of the input (the external source term in this context, e.g., F (s)) and the Laplace transform of the output (the unknown field variable in this context, e.g., D(s)); namely, D(s) = G(s)F (s). For details, see [36]. The frequency response of continuous-time transfer function G(s) is calculated by the complex quantity G(j
) that is obtained by substituting s = j
. The magnitude and phase of G(j
) are, respectively, the magnitude and phase of the output d(t) of the associated system at steady-state in response to the sinusoidal input F (t) of frequency
[36]. The frequency domain identification problem is defined as a nonlinear least squares problem as follows: given the complex output Dk and complex input Fk values at frequencies
k for k = 1, 2, . . . l, find the transfer function G(j
k , P ) such that the following least squares type cost function is minimized: CLS =

l

2

Fk N G (j
k , P ) − Dk D G (j
k , P ) ,

(32)

k=1

where P is the unknown parameter vector, N G (j
k , P ) and D G (j
k , P ) are, respectively, the numerator and denominator of G(j
k , P ), i.e., G(j
k , P ) = N G (j
k , P )/D G (j
k , P ). Note that if N G (j
k , P ) and D G (j
k , P ), respectively, equal to the numerator and denominator of the original model, then CLS = 0. In the algorithm used in this work, firstly orders of the polynomials for N G (j
k , P ) and D G (j
k , P ) are selected on a trial-and-error basis, then the Newton–Raphson method is used to minimize (32), e.g., see [38,31,32]. Once, the parameters in G(j
, P ) are determined, the reduced order model is obtained by constructing the following state space representation by using the poles and zeros of estimated model [36]: x˙ = Ax + Bf ,

(33)

q = Ex,

(34)

where x is the state vector of the structure, the matrices A, B and E can be in a canonical form [39]. In sum, (33) and (34) represent the reduced order model for the original system (1). Orders of the polynomials for N G (j
k , P ) and D G (j
k , P ) are determined by comparing the corresponding cost functions and locations of poles and zeros. For instance, the orders of polynomials for N G (j
k , P ) and D G (j
k , P ) are increased until sufficiently small cost function values CLS are obtained which is selected to be in the order of 10−3 –10−5 in this work depending on the frequency range of interest. Once, the orders of N G (j
k , P )

Y. Cunedio˘glu et al. / Finite Elements in Analysis and Design 42 (2006) 367 – 403

378

and D G (j
k , P ) are decided to be appropriate, poles and zeros of the reduced model are studied. If any pole (or zero) is at least five times faster than the other poles (or zeros), the orders of N G (j
k , P ) and D G (j
k , P ) are reduced as many as the numbers of such poles or zeros. This method is applied to a sample problem in Section 6.

5. Subspace-based identification method The subspace based identification algorithm [40] estimates the parameters of the state space representation of a transfer function. In this work, the parameters of the discrete-time system in the state space form are estimated to obtain reduced order models for structural systems. Then, by using these parameters, one could readily obtain the frequency and time domain responses of the corresponding discrete-time system. In the algorithm, the row size of a square Hankel matrix is entered by the user; but it is suggested to be at most half of the number of sample continuous-time frequency response data. Also, the model orders are chosen by the user. Details of this algorithm can be found in [40]. 5.1. Problem formulation Assume that G(z) is a transfer function of a stable, multi-input multi-output (MIMO), linear, timeinvariant, discrete-time system having the input–output relationship characterized by the impulse response coefficients gk through the following equation [40]: y(t) =

∞

gk u(t − k),

(35)

k=0

where, z is the complex Z-transform variable, the output y(t) ∈ Rp , the input u(t) ∈ Rm and gk ∈ Rp×m . It is also assumed that the system has the finite order of n and can thus be described by a state-space model as follows: x(t + 1) = Ax(t) + Bu(t), y(t) = Cx(t) + Du(t),

(36)

where y(t) ∈ Rp , u(t) ∈ Rm and the state vector x(t) ∈ Rn . The state space model (36) has the following impulse response coefficients  D, k = 0, (37) gk = CAk−1 B, k > 0. Then, the frequency response of the system having the impulse response coefficients gk can be calculated by ∞   G e j
= gk e−j
k , k=0

0


.

(38)

Y. Cunedio˘glu et al. / Finite Elements in Analysis and Design 42 (2006) 367 – 403

379

Considering the state-space realization (36), then (38) can equivalently be written as follows:    −1 G ej
= C ej
I − A B + D.

(39)

5.2. Uniformly spaced data Note that only the case of uniformly spaced data is considered in this study, which is typical in frequency domain identification problems. Suppose that M + 1 frequency response data Gk on the following set of uniformly spaced frequencies:
k =

k

M

,

k = 0, . . . , M

(40)

are given. If the impulse response coefficients (37) are known, then realization algorithms can be used to obtain a state-space realization. The subspace based identification algorithm [40] uses the coefficients of the inverse discrete Fourier transform (IDFT) obtained from the samples of the frequency response function. Since G(z) is a transfer function with a real valued impulse response (35), the frequency response data on the frequency range [0, ] can be extended to the frequency range [, 2] by taking the complex conjugate of the given data Gk which forms the first step of the identification algorithm. The subspace based identification algorithm used in this work is given below: Algorithm: (1) Extend the transfer function samples to the full unit circle: GM+k := G∗M−k ,

k = 0, . . . , M,

(41)

where the superposed asterisk denotes complex conjugate. (2) Let hˆ i be defined by the 2M-point IDFT as follows: 2M−1 1 Gk ej2ik/2M , hˆ i := 2M

i = 0, . . . , 2M − 1.

(42)

k=0

(3) Let the block Hankel matrix Hˆ be defined as follows: ⎡ˆ h1 ⎢ hˆ 2 Hˆ := ⎢ ⎣ ... hˆ q

hˆ 2 hˆ 3 .. . hˆ q+1

⎤ hˆ r hˆ r+1 ⎥ .. ⎥ . ⎦ · · · hˆ q+t−1 ··· ··· .. .

∈ Rqp×rm ,

(43)

where the number of block rows q > n and the number of block columns r  n. The dimension of Hˆ is bounded by q + r
2M.

Y. Cunedio˘glu et al. / Finite Elements in Analysis and Design 42 (2006) 367 – 403

380

(4) Calculate the following singular value decomposition (SVD) of the Hankel matrix Hˆ : Hˆ = Uˆ ˆ Vˆ T . (5) Determine the system order n by inspecting the singular values and partition the SVD such that ˆ s contains the n largest singular values defined by

T ˆs 0  Vˆ Hˆ = [ Uˆ s Uˆ 0 ] . (44) ˆ 0 0 Vˆ0T (6) Determine the system matrices Aˆ and Cˆ as follows: Aˆ = (J1 Uˆ s )† J2 Uˆ s ,

(45)

Cˆ = J3 Uˆ s ,

(46)

where J1 =  I(q−1)p

0(q−1)p×p  ,

J2 =  0(q−1)p×p J3 =  Ip

(47)

I(q−1)p  ,

(48)

0p×(q−1)p  ,

(49)

and Ii denotes the i × i identity matrix, 0i×j the i × j zero matrix, and X† = (X T X)−1 XT the Moore–Penrose pseudo-inverse of the full column rank matrix X. ˆ (7) Solve the following least squares problem to determine Bˆ and D: ˆ Dˆ = arg min B,

B∈Rn×m D∈Rn×m

M

2

ˆ j
k I − A) ˆ −1 B

,

Gk − D − C(e k=0

  where XF = k s |xks |2 denotes the Frobenius norm. (8) The estimated transfer function is defined by  −1 ˆ M (z) = Dˆ + Cˆ zI − Aˆ ˆ G B.

F

(50)

(51)

Notice that B and D appear linearly in the transfer function for fixed A and C. Hence, the optimization problem (50) has an analytical solution. 5.3. Quality measures To assess the quality of estimated models, the following two measures for the fit between the sampled data and estimated model are used: (i) The maximum error defined by  

 ˆM 

ˆ M j
k

= max G (e ) − Gk . G − G m,∞


k

(52)

Y. Cunedio˘glu et al. / Finite Elements in Analysis and Design 42 (2006) 367 – 403

381

(ii) The root-mean-square error (rms) given by ˆM G

  N

2 1

ˆ M j
k

 − Gm,2 =

G (e ) − Gk . N

(53)

k=1

These measures are used to select the model orders in the sample problem considered in Section 6.1.

6. Numerical studies In this section, by considering a truss structure, frequency and time domain responses of the structure are studied by using the model order reduction techniques such as CMS, QSM, NLS and SBI methods. Even though CMS and QSM methods have been implemented in commercial codes such as MSC Nastran, the NLS and SBI methods do not exist in any commercial code; thus, performances of these methods are compared with each other by using programs developed in Matlab environment. Note that convergence studies are completed to guarantee the accuracy of FEM solutions by changing the size of time step in time integration (Newmark) algorithm and number of samples in frequency domain data; however, these results are not presented here for limited space and only the converged results are given. In order to compare the CMS and QSM methods with the NLS and SBI methods on the same ground, unless it is declared, numbers of modes in the CMS and QSM methods are selected as half of the model orders of the model estimated by the SBI method. 6.1. Frequency domain analyses Consider the planar truss system shown in Fig. 4 modeled by using truss elements [41]. Note that FEM equations are exact for this structure. In FEM models used in simulations, the truss system is formed of 40, 400 and 1000 elements. For frequency domain identification (i.e., the NLS and SBI methods), sampling frequency points are selected as 512, 1024 and 2048 for 40, 400 and 1000 elements, respectively. The results are given in detail for 40 elements for limited space. The number of elements and node numbers are denoted by nel and nnode in FEM models, respectively. The truss system is fixed at the node numbers 1 and 12, and excited by the vertical sinusoidal forces vertically at the node numbers 11 and 22. The applied force is given by f (t) = 1000 Sin
t[N]. In this study, we chose the active degree of freedoms as the displacements in the y-direction at the node numbers 1, 6, 11, 12 and 22, e.g., see Fig. 4. F=1000Sin
t 12

13

14

15

16

17

18

19

20

21

22

F=1000Sin
t 1

2

3

4

5

6

7

8

9

Fig. 4. Truss system having 40 elements.

10

11

382

Y. Cunedio˘glu et al. / Finite Elements in Analysis and Design 42 (2006) 367 – 403 2

10

0

10

-2

Amplitude

10

-4

10

-6

10

-8

10

-10

10

-1

10

0

10

1

10

2

10

3

10

ω (rad/sec)

Fig. 5. Frequency response magnitude plot at the node number 6 for the exact system where nel = 40.

The structural parameters in the FEM models are chosen as follows: the cross sectional area of truss elements is h = 1 cm2 , elasticity modulus E = 20000000 N/cm2 , density of the truss element = 0.007850 kg/cm3 , length of the vertical and horizontal truss elements l1 = 100 cm, length of the diagonal truss elements l2 = 144 cm and 40 kg lumped mass is applied to each DOF. Note that the damping in FEM models is defined as Rayleigh damping in the form of C = cK, where C is the global damping matrix and c is a constant. However, both damped and undamped cases are considered in numerical examples. For undamped cases, the constant c is chosen to be zero c = 0 and these situations are explicitly expressed. Then, the NLS, CMS and QSM1 methods are examined in terms of the frequency response magnitude, phase and phase error characteristics. Dynamic characteristics of the truss structure can be seen in the Bode magnitude plot in Fig. 5 where vertical displacement at the node number 6 is shown for the frequency range 0 to 200 rad/s. For vertical displacement at the node number 6, frequency response magnitude plots of the FEM model, CMS and QSM1 methods in the frequency range of 0–60 rad/s are given in Fig. 6, where the number of modes for CMS and QSM1 methods to compute the reduced order model is nmode = 4 (the natural frequencies of the first four modes are 30.97, 35.38, 43.14 and 44.86 rad/s). The corresponding frequency response phase and phase error plots of the FEM model, CMS and QSM1 methods are, respectively, given in Figs. 7 and 8, where the curves for CMS and QSM1 methods overlap. During these studies, centering frequencies for the QSM1 and QSM2 methods are selected as follows:
c1 = 20 rad/s for QSM1,
c2 = 40 rad/s for QSM2. To compare CMS, QSM1 and QSM2 methods, equal numbers of modes are considered; thus, QSM2 is not further investigated in this work since it employs higher number of modes due to consideration of the second centering frequency. Since there is no damping in the system, phase plots in Fig. 7 are formed of

Y. Cunedio˘glu et al. / Finite Elements in Analysis and Design 42 (2006) 367 – 403

0

383

FEM Cms Qsm1

10

-1

10

Amplitude

-2

10

-3

10

-4

10

-5

10

0

1

10

10 ω (rad/sec)

Fig. 6. Frequency response magnitude plots of the FEM model, CMS and QSM1 methods at the node number 6 in the frequency range 0–60 rad/s for c = 0 and nmode = 4 where nel = 40. 200 FEM Cms Qsm1

180

Phase angle in degrees

160 140 120 100 80 60 40 20 0 0

10

1

ω (rad/sec)

10

Fig. 7. Frequency response phase plots of the exact FEM model, CMS and QSM1 methods at the node number 6 in the frequency range 0–60 rad/s for c = 0 and nmode = 4 where nel = 40.

384

Y. Cunedio˘glu et al. / Finite Elements in Analysis and Design 42 (2006) 367 – 403 200 cms error qsm1 error

150

Phase error in degrees

100 50 0 -50 -100 -150 -200

0

1

10

10 ω (rad/sec)

Fig. 8. Frequency response phase error plots of CMS and QSM1 methods at the node number 6 in the frequency range 0–60 rad/s and nmode = 4 where nel = 40.

pulses at resonance frequencies. Observe in Fig. 6 that frequency response magnitude plots of CMS and QSM1 methods do not match with the FEM solution at high resonance frequencies; subsequently, large frequency response phase errors occur at these resonance frequencies. It can be seen in Fig. 7 that phase angle plots have jumps between 180◦ and 0◦ since no damping exists in the system. Following, the NLS method is applied to the truss structure to obtain a reduced order model whose output is the vertical displacement at the node number 6. Before using the NLS method, the orders of numerator nord and denominator dord of continuous-time transfer function in (30) are to be defined by the user. A detailed discussion on the selection of the orders of numerator and denominator can be found in [30]. Similar to nord and dord in the NLS method, the parameter nmode , the number of component normal modes n , should be selected by the user for the CMS method used in the transformation matrix T. Similarly, for the QSM methods, the numbers of normal modes and centering frequencies should be selected by the user as well. In this study, a simple criterion is followed for the NLS method: as nord and dord are increased, the order of the reduced order model is determined on a trial-and-error basis by evaluating the changes in the sum of squares of errors in the cost CLS and by interpreting the locations of poles and zeros of the reduced order models [30]. By considering 512 sampling frequency points, the frequency range in which the NLS method is used to estimate the reduced order model is changed as 0 to 60, 0 to 90, 0 to 120 and 0 to 200 rad/s and the estimated models are examined to determine the best fit. As the best fit, the numerator and denominator of the estimated model obtained by the NLS method in the frequency range 0–60 rad/s are found to be 8 and 9, respectively. Subsequently, according to the NLS method, the number of modes in the frequency range 0–60 rad/s is 4 oscillatory modes plus 1 time

Y. Cunedio˘glu et al. / Finite Elements in Analysis and Design 42 (2006) 367 – 403

385

100

FEM NLS -1

10

-2

Magnitude

10

-3

10

-4

10

-5

10

-6

10

-7

10

0

1

10

10

ω (rad/sec)

Fig. 9. Frequency response magnitude plots of the FEM model and NLS method in the frequency range 0–60 rad/s where nel =40.

constant mode in the reduced order model. Frequency response magnitude and phase error plots of the NLS solution and FEM model are given in Figs. 9 and 10, respectively. It can be seen in Fig. 9 that frequency response magnitude plots of the NLS solution and FEM model overlap perfectly. Therefore, frequency response phase errors are nearly zero, e.g., see Fig. 10. Note that when equal numbers of modes are used for the NLS, CMS and QSM1 methods in the frequency range 0–60 rad/s, phase errors of the NLS method is much less than the CMS and QSM1 methods. Then, the SBI method is applied to the same truss structure to obtain a reduced order model by choosing the output as the vertical displacement at the node number 6. Similar to the NLS method, 512 sampling frequency points are considered and the frequency range in which the SBI method is used to estimate the reduced order model is changed as 0 to 60, 0 to 90, 0 to 120 and 0 to 200 rad/s. As the order of the estimated state space realizations is increased, the corresponding root-mean-square errors defined by (53) are examined to choose the best fit. Consequently, the model order of the reduced system is found to be 16 as the best fit when the SBI method is used in the frequency range 0–60 rad/s, e.g., see Fig. 11 for the frequency response magnitude plot. Note that similar to the NLS solution, the corresponding frequency response phase error is almost zero that is not presented here for limited space. In Tables 1–3, as the Rayleigh damping coefficient varies, the CPU times of the four methods and estimated model orders for the NLS and SBI methods are presented, where nmode = 4 (the forth natural frequency is 44.86 rad/s) for the CMS and QSM1 methods and FEM models have 40, 400 and 1000 elements, respectively. Observe in Tables 1–3 that when the frequency ranges are increased for both damped and undamped cases, the orders of reduced models obtained by the SBI and NLS methods increase. The NLS method did not converge in the frequency ranges 0–90, 0–120 and 0–200 rad/s for no

386

Y. Cunedio˘glu et al. / Finite Elements in Analysis and Design 42 (2006) 367 – 403 200 150

Phase error in degrees

100 50 0 -50 -100 -150 -200

0

1

10

10 ω (rad/sec)

Fig. 10. Frequency response phase error plot of the NLS method in the frequency range 0–60 rad/s where nel = 40.

0

10

FEM SBI -1

10

-2

Amplitude

10

-3

10

-4

10

-5

10

-6

10

-1

10

0

1

10

10

2

10

ω (rad/sec)

Fig. 11. Frequency response magnitude plots of the FEM model and SBI method in the frequency range 0–60 rad/s for no damping c = 0 where nel = 40 and model order is 16.

Y. Cunedio˘glu et al. / Finite Elements in Analysis and Design 42 (2006) 367 – 403

387

Table 1 Frequency domain analyses for 40 elements Damping Freq. range coeff. (rad/s)

SBI model order

SBI time (s)

NLS num/den order

c=0

16 26 32 40 12 14 16 12 6 8 8 10

8.156 8.078 8.468 6.437 10.87 10.703 10.266 7.468 13.844 12.875 10.969 8.891

8/9 1.64 No convergence

c = 0.01

c = 0.05

0–60 0–90 0–120 0–200 0–60 0–90 0–120 0–200 0–60 0–90 0–120 0–200

NLS time (s)

7/9 12/13 15/17 28/30 4/6 6/8 8/9 8/10

2.969 2.25 2.172 5.172 5.657 9 2.141 2.375

FEM time (s)

CMS time (s)

QSM1 time (s)

2.65 × 10−1 2.97 × 10−1 2.97 × 10−1 2.97 × 10−1 3.28 × 10−1 3.43 × 10−1 3.59 × 10−1 4.22 × 10−1 3.44 × 10−1 3.75 × 10−1 3.44 × 10−1 3.44 × 10−1

2.66 × 10−1 2.03 × 10−1 2.19 × 10−1 2.03 × 10−1 2.65 × 10−1 2.04 × 10−1 2.18 × 10−1 2.01 × 10−1 2.68 × 10−1 2.01 × 10−1 2.15 × 10−1 2.04 × 10−1

2.65 × 10−1 2.03 × 10−1 2.19 × 10−1 2.18 × 10−1 2.64 × 10−1 2.05 × 10−1 2.17 × 10−1 2.19 × 10−1 2.63 × 10−1 2.02 × 10−1 2.18 × 10−1 2.16 × 10−1

Table 2 Frequency domain analyses for 400 elements Damping coeff.

Freq. range (rad/s)

SBI model order

SBI time time (s)

NLS num/den order

c=0

0–60 0–90 0–120 0–200 0–60 0–90 0–120 0–200 0–60 0–90 0–120 0–200

60 Too large model orders

43.75

No convergence

24 26 28 80 8 10 12 18

35.188 34.734 30.703 46.672 86.078 66.484 48.141 32.281

c = 0.01

c = 0.05

NLS time (s)

12/13 2.828 16/18 2.047 18/20 3.703 No convergence

FEM time (s)

CMS time (s)

QSM1 time (s)

115 Not computed

3.75

3.672

327.55 326.01 327.84 327.69 326.24 326.39 326.38 326.16

3.74 3.73 3.75 3.80 3.77 3.75 3.78 3.79

3.670 3.669 3.670 3.672 3.673 3.672 3.675 3.671

Table 3 Frequency domain analyses for 1000 elements Damping Freq. range coeff. (rad/s)

SBI model order

c=0 c = 0.01

200 440.094 36 270.546 84 343.812 18 512.02 20 316.88 52 246.438 Too large model order

c = 0.05

0–60 0–60 0–90 0–60 0–90 0–120 0–200

SBI time (s)

NLS num/den order No convergence

NLS time (s)

FEM time (s)

CMS time (s)

QSM1 time (s)

2614.9 8869.7 8809.5 8878.2 8841.4 8841.7 Not computed

43.016 43.015 43.017 43.014 43.019 43.020

43.032 43.031 43.034 43.029 43.030 43.032

388

Y. Cunedio˘glu et al. / Finite Elements in Analysis and Design 42 (2006) 367 – 403

damping c = 0. In general, the orders of reduced models for the system having damping are less than those having no damping. It is also noteworthy that for all these methods when damping coefficient is increased, the orders of reduced models decrease. In general, it can be observed that CPU times of the SBI and NLS methods decrease as the damping coefficient is increased. Because of the nonlinear nature of the Newton Raphson method, the CPU times of the NLS method for some frequency ranges may be very high. If CPU times of these four methods are compared, it is seen that the CMS and QSM1 methods are much faster than the NLS and SBI methods; the reason for this result is that the NLS and SBI methods have not been originally developed for model order reduction purposes and their computational load is higher. It is observed in Tables 1–3 that for large number of elements (e.g., nel = 400 and 1000), the NLS method may have convergence problems, and the SBI method may yield estimated models having large orders for a good fit in frequency responses. 6.2. Time domain analyses The time domain response of the FEM model is calculated by using the Newmark method, where the parameters and  are selected to be = 0 and  = 0.5, i.e., central difference formula [41]. The vertical displacement at the node number 6 is considered to be the output in all figures and tables. In numerical integrations, the time step
t is selected to be 0.002 s and the response is computed between 0 and 3 s for 40, 400 and 1000 elements. The results are given in detail for 40 elements for limited space. If the NLS method did not converge and the SBI method needed large model orders for a good fit in a frequency range, the CMS and QSM1 methods would not computed at these frequency ranges either. In order to compare all these methods on the same ground, when the tables for time domain analyses are formed, only the frequency ranges in which good agreement with the FEM solutions is observed are considered. The exact and FEM model step responses for no damping c = 0 case are presented in Figs. 12 and 13, respectively. For damped and undamped (c = 0.05) cases, step response plots of the model obtained by the NLS method estimated in the frequency range 0–60 rad/s are given in Figs. 14 and 15, respectively. It is observed that the NLS method may give extremely good fit to the frequency domain data (e. g., see Section 6.1), whereas it may yield estimated models having unstable poles and/or nonminimum phase zeros, whose step response is shown in Fig. 15; therefore, time domain response of the model estimated by the NLS method may not be satisfactory. On the other hand, as the number of modes is increased, time domain response plots for the CMS and QSM1 methods converge to the exact solution as observed in Figs. 16–19. The step responses are shown in Figs. 16 and 17, while the sinusoidal force responses where the excitation is in the form of 1000(Sin(t) + 0.2 Sin(5t)) are shown in Figs. 18 and 19. The number of modes nmode = 4 (the 4th natural frequency is 44.86 rad/s) in Figs. 16 and 18, while the number of modes nmode = 20 (the 20th natural frequency is 79.91 rad/s) in Figs. 17 and 19. It is observed that step response of the SBI solutions does not agree with the FEM solution. For instance, step responses of the SBI solutions given in Figs. 21 and 24 which correspond to systems estimated in Figs. 20 and 23 do not agree with the corresponding step responses of the FEM models given in Figs. 12 and 22 that are for c = 0 and 0.05, respectively. To reveal the effect of the frequency band used for the SBI method to estimate a reduced order model, instead of lumped mass, material density is increased to 10 kg/cm3 . In this case, dynamics of the system is in the frequency range 0–30 rad/s. The step response of the corresponding FEM model obtained by using the Newmark method is shown in Fig. 25. When the frequency range to estimate a model by using the SBI method is selected between 0 to 20 rad/s, the frequency response magnitude plot and step response of the

Y. Cunedio˘glu et al. / Finite Elements in Analysis and Design 42 (2006) 367 – 403

389

-3

x 10

2 1.5 1

Magnitude

0.5 0 -0.5 -1 -1.5 -2

0

0.5

1

1.5 Time (sec)

2

2.5

3

Fig. 12. Step response of the FEM model obtained by using Newmark method for no damping c = 0 where nel = 40.

-3

2

x 10

1.5

1

Amplitude

0.5

0

-0.5

-1

-1.5

-2

0

0.5

1

1.5

2

2.5

3

Time (sec)

Fig. 13. Exact step response of the truss system for no damping c = 0 where nel = 40.

390

Y. Cunedio˘glu et al. / Finite Elements in Analysis and Design 42 (2006) 367 – 403

x 10

-7

10 8

Amplitude

6 4 2 0 -2 -4 0

0.5

1

1.5

2

2.5

3

Time (sec)

Fig. 14. Step response of the NLS method where the model is estimated in the frequency range 0–60 rad/s for c = 0 and nel = 40. 5

x 10 0 -2 -4

Amplitude

-6 -8 -10 -12 -14 -16 -18 0

0.5

1

1.5

2

2.5

3

Time (sec)

Fig. 15. Step response of the NLS method having unstable poles where the model is estimated in the frequency range 0–60 rad/s for c = 0.05 and nel = 40.

Y. Cunedio˘glu et al. / Finite Elements in Analysis and Design 42 (2006) 367 – 403

391

-3

4

x 10

FEM Cms Qsm1

3 2

Magnitude

1 0 -1 -2 -3 -4

0

0.5

1

1.5 Time (sec)

2

2.5

3

Fig. 16. Step responses of the FEM model, CMS and QSM1 methods for c = 0 and nmode = 4 where nel = 40. 2.5

x 10-3 FEM Cms Qsm1

2 1.5

Magnitude

1 0.5 0 -0.5 -1 -1.5 -2 -2.5

0

0.5

1

1.5

2

2.5

3

Time (sec)

Fig. 17. Step responses of the FEM model, CMS and QSM1 methods for c = 0 and nmode = 20 where nel = 40.

392

Y. Cunedio˘glu et al. / Finite Elements in Analysis and Design 42 (2006) 367 – 403 -4

8

x 10

FEM Cms Qsm1

6 4

Magnitude

2 0 -2 -4 -6 -8

0

0.5

1

1.5 Time (sec)

2

2.5

3

Fig. 18. Sinusoidal force responses of the FEM model, CMS and QSM1 methods for c = 0 and nmode = 4 where nel = 40.

6

x 10

-4

exact FEM cms Cms qsm1 Qsm1

4

Magnitude

2

0

-2

-4

-6

0

0.5

1

1.5 Time (sec)

2

2.5

3

Fig. 19. Sinusoidal force responses of the FEM model, CMS and QSM1 methods for c = 0 and nmode = 20 where nel = 40.

Y. Cunedio˘glu et al. / Finite Elements in Analysis and Design 42 (2006) 367 – 403

393

1

10

FEM SBI 0

10

-1

Amplitude

10

-2

10

-3

10

-4

10

-5

10

-6

10

-1

0

10

1

10

2

10

3

10

10

ω (rad/sec)

Fig. 20. Frequency response magnitude plots of the FEM model and SBI method in the frequency range 0–120 rad/s for c = 0 where nel = 40 and the model order is 32. -3

1

x 10

Amplitude

0.5

0

-0.5

-1

-1.5 0

0.5

1

1.5 Time (sec)

2

2.5

3

Fig. 21. Step response of the SBI method for the reduced order model estimated in Fig. 20 where nel = 40.

394

Y. Cunedio˘glu et al. / Finite Elements in Analysis and Design 42 (2006) 367 – 403 -5

4.5

x 10

4 3.5

Magnitude

3 2.5 2 1.5 1 0.5 0

0

0.5

1

1.5 Time (sec)

2

2.5

3

Fig. 22. Step response of the FEM model for c = 0.05 where nel = 40. -4

10

Amplitude

FEM SBI

-5

10

-6

10

-1

10

0

10

1

10

2

10

3

10

ω (rad/sec)

Fig. 23. Frequency response magnitude plots of the FEM model and SBI method in the frequency range 0–120 rad/s for c = 0.05 where nel = 40 and model order is 8.

Y. Cunedio˘glu et al. / Finite Elements in Analysis and Design 42 (2006) 367 – 403

4.5

x 10

395

-5

4 3.5

Amplitude

3 2.5 2 1.5 1 0.5 0 -0.5 0

0.5

1

1.5 Time (sec)

2

2.5

3

Fig. 24. Step response of the SBI method for the system estimated in Fig. 23 where nel = 40 and model order is 8. -3

1.5

x 10

1

Magnitude

0.5

0

-0.5

-1

-1.5

0

0.5

1

1.5

2

2.5

3

Time (sec)

Fig. 25. Step response of the FEM model for the modified system for c = 0 where nel = 40.

Y. Cunedio˘glu et al. / Finite Elements in Analysis and Design 42 (2006) 367 – 403

396

1

10

FEM SBI 0

10

-1

Amplitude

10

-2

10

-3

10

-4

10

-5

10

-6

10

-2

-1

10

0

10

2

1

10

10

10

ω (rad/sec)

Fig. 26. Frequency response magnitude plots of the FEM model and SBI method for the modified system in the frequency range 0–20 rad/s for c = 0 where nel = 40 and model order is 32. -3

2

x 10

1.5

Amplitude

1

0.5

0

-0.5

-1

-1.5

0

0.5

1

1.5

2

2.5

3

Time (sec)

Fig. 27. Step response of the SBI method for the system estimated in Fig. 26 where nel = 40 having the model order of 32.

16 26 32 40

0–60 0–90 0–120 0–200

FEM time (s) 6.56 × 10−1

SBI time (s) 2.65 × 10−1 2.97 × 10−1 2.5 × 10−1 3.6 × 10−1

1.25 × 10−1 1.72 × 10−1 2.34 × 10−1 2.81 × 10−1

CMS time (s) 1.1 × 10−1 1.87 × 10−1 2.03 × 10−1 2.66 × 10−1

QSM1 time (s) 6.2721 × 10−4 1.1464 × 10−2 9.7576 × 10−5 7.2859 × 10−10

SBI Max error

0–60 0–90 0–120 0–200

Freq.

16 26 32 40

FEM time (s)

2.65 × 10−1 6.72 × 10−1 2.97 × 10−1 2.5 × 10−1 3.6 × 10−1

SBI model SBI time range (rad/s) order (s) 1.4 × 10−1 1.72 × 10−1 2.35 × 10−1 2.66 × 10−1

CMS time (s) 1.09 × 10−1 1.72 × 10−1 2.34 × 10−1 2.5 × 10−1

QSM1 time (s) 6.2721 × 10−4 1.1464 × 10−2 9.7576 × 10−5 7.2859 × 10−10

SBI Max error

Table 5 Sinusoidal response absolute errors and CPU times for 40 elements and no damping c = 0

order

(rad/s)

Freq. range SBI model

Table 4 Step response absolute errors and CPU times for 40 elements and no damping c = 0

1.2871 × 10−4 1.2488 × 10−3 1.4133 × 10−5 1.2353 × 10−10

SBI RMS error

1.2871 × 10−4 1.2488 × 10−3 1.4133 × 10−5 1.2353 × 10−10

SBI RMS error

6.8893 × 10−4 2.1384 × 10−4 3.6178 × 10−4 4.7126 × 10−4

CMS Max error

2.9303 × 10−3 8.1657 × 10−4 1.4108 × 10−3 1.8041 × 10−3

CMS Max error

8.0635 × 10−4 1.9634 × 10−4 2.99 × 10−4 3.0624 × 10−4

QSM1 Max error

3.5741 × 10−3 7.6715 × 10−4 1.1532 × 10−3 1.6157 × 10−3

QSM1 Max error

Y. Cunedio˘glu et al. / Finite Elements in Analysis and Design 42 (2006) 367 – 403 397

Freq. range (rad/s) 0–60 0–60 0–90 0–120 0–200 0–60 0–90 0–120 0–200

c=0 c = 0.01

c = 0.05

coef.

Damp.

60 24 26 28 80 8 10 12 18

SBI model order 0.657 0.422 0.344 0.328 1.031 0.281 0.281 0.25 0.344

SBI time (s)

112.079

113.782 113.313

FEM time (s) 0.438 0.157 0.157 0.188 0.609 0.078 0.094 0.094 0.125

CMS time (s)

Table 6 Step response absolute errors and CPU times for 400 elements

0.39 0.141 0.140 0.187 0.641 0.078 0.093 0.093 0.125

QSM1 time (s) 0.0017 8.6977 × 10−11 7.7171 × 10−16 6.8149 × 10−16 5.4981 × 10−16 5.1257 × 10−15 8.2877 × 10−17 3.368 × 10−18 8.2061 × 10−22

SBI Max error

CMS Max error 0.0023 4.9556 × 10−7 4.9556 × 10−7 3.0328 × 10−5 6.4305 × 10−7 5.2338 × 10−8 5.2338 × 10−8 2.6424 × 10−6 3.9061 × 10−8

SBI RMS error 1.6056 × 10−4 7.1105 × 10−12 2.7851 × 10−16 1.707 × 10−16 6.5138 × 10−17 8.8183 × 10−16 1.2541 × 10−17 1.0234 × 10−18 3.1349 × 10−22

0.0023 1.5514 × 10−5 8.0475 × 10−4 4.2301 × 10−4 8.0331 × 10−5 8.4205 × 10−6 5.2567 × 10−4 3.0194 × 10−4 7.0538 × 10−5

QSM1 Max error

398 Y. Cunedio˘glu et al. / Finite Elements in Analysis and Design 42 (2006) 367 – 403

Y. Cunedio˘glu et al. / Finite Elements in Analysis and Design 42 (2006) 367 – 403

399

Table 7 Sinusoidal response absolute errors and CPU times for 400 elements and no damping c = 0 Freq. range (rad/s)

SBI model order

SBI time (s)

FEM time (s)

CMS time (s)

QSM1 time (s)

SBI Max error

SBI RMS error

CMS Max error

QSM1 Max error

0–60

60

3.516

113.062

0.438

0.422

0.0017

1.6056 × 10−4

3.3212 × 10−4

5.1817 × 10−4

SBI solution are given in Figs. 26 and 27, respectively. Observe in Fig. 26 that static gain of frequency response is not correct. Although characteristics of the step response of the model obtained by the SBI method is similar to that of the FEM model, magnitude values of step responses do not match due to this wrong static gain. When the number of elements in FEM models is increased to 400 and 1000, it is observed that static gain of the frequency response magnitude plots for the SBI solution has even larger errors for both damped and undamped cases; consequently, step response of the estimated model obtained by the SBI method is largely in error in comparison with step response of the FEM models. The CPU times and maximum RMS errors and/or maximum errors in response to step and sinusoidal excitations are listed in Tables 4–9 for the four methods for different damping coefficients. It is observed in Tables 4, 6 and 8 that when the number of elements and damping coefficients are increased, the maximum RMS errors and maximum errors decrease for the SBI method. The CMS and QSM1 errors decrease as the number of modes is increased, e.g., see also Figs. 16–19. It is noteworthy that the CMS and QSM1 methods are much faster than the NLS and SBI methods. If Tables 5, 7 and 9 obtained for sinusoidal force excitation where the excitation is in the form of 1000(Sin(t) + 0.2 Sin(5t)) are compared with Tables 4, 6 and 8 obtained for step responses, it is observed that the CMS and QSM1 errors for sinusoidal force excitation are less than those of the CMS and QSM1 errors of step input.

7. Conclusion As a model order reduction method, a nonlinear least squares fit in frequency domain and subspace based identification algorithm are compared with the CMS and QSM1 methods. It is noteworthy that the nature of these methods is different; while the NLS and SBI methods are developed to estimate models for given frequency domain data, the CMS and QSM methods are developed to fit reduced order models by using certain modes extracted from the structure itself; therefore, the accuracy of the CMS and QSM methods in time domain responses might be better than that of the NLS and SBI methods, even if the accuracy of the CMS and QSM methods in frequency domain responses were worse than that of the NLS and SBI methods. It is observed in the numerical examples that convergence of the NLS method can be achieved for low model orders and convergence difficulties may arise as the problem size increases. Nonetheless, it is concluded that the NLS and SBI methods are more accurate than the CMS and QSM1 methods in terms of frequency response magnitude and phase errors. On the other hand, even though it is possible to obtain models estimated by the SBI method having very good agreement with given frequency domain data, the model orders may be very large for a satisfactory fit to the exact time domain responses, especially if the problem size is large. Unlike the NLS method, it is not observed that the SBI method yields estimated models having unstable and/or nonminimum phase

Freq. range (rad/s) 0–60 0–60 0–90 0–60 0–90 0–120

c=0 c = 0.01

c = 0.05

coef.

Damp.

200 36 84 18 20 52

SBI model order 6.765 0.407 0.563 0.281 0.282 0.375

SBI time (s)

1551.3

4658.3 2644.6

FEM time (s) 4.203 0.234 0.75 0.125 0.125 0.344

CMS time (s)

Table 8 Step response absolute errors and CPU times for 1000 elements

4.562 0.219 0.688 0.141 0.125 0.329

QSM1 time (s)

SBI RMS error 5.4757 × 10−5 9.4199 × 10−22 2.8813 × 10−23 1.8975 × 10−32 2.4481 × 10−35 3.6837 × 10−37

SBI Max error 7.1464 × 10−4 1.2257 × 10−20 7.2073 × 10−23 2.2268 × 10−31 1.7232 × 10−34 1.8327 × 10−36

3.0998 × 10−4 1.4882 × 10−5 1.0765 × 10−5 1.7662 × 10−8 1.0770 × 10−6 9.5784 × 10−7

CMS Max error

2.7743 × 10−4 8.7199 × 10−5 3.5053 × 10−4 1.3997 × 10−6 2.1273 × 10−4 7.2455 × 10−5

QSM1 Max error

400 Y. Cunedio˘glu et al. / Finite Elements in Analysis and Design 42 (2006) 367 – 403

Y. Cunedio˘glu et al. / Finite Elements in Analysis and Design 42 (2006) 367 – 403

401

Table 9 Sinusoidal response absolute errors and CPU times for 1000 elements and no damping c = 0 Freq. range (rad/s)

SBI model order

SBI time (s)

FEM time (s)

CMS time (s)

QSM1 (s)

SBI Max error

SBI RMS error

CMS Max error

QSM1 Max error

0–60

200

6.765

4809

4.484

4.203

7.1464 × 10−4

5.4757 × 10−5

1.0187 × 10−6

1.054 × 10−6

zeros. In numerical studies, reduced model orders of the NLS method were always less than those of the SBI method in all frequency ranges for both damped and undamped cases and it can be seen that model orders of the reduced system decrease as the damping coefficient increases. In terms of the CPU times, the CMS and QSM1 methods are much faster than the NLS and SBI methods in both frequency and time domain analyses. For the CMS and QSM1 methods, convergence of time domain response to that of the exact FEM model solutions can be achieved by increasing the numbers of modes. Convergence characteristics in time domain for the CMS and QSM1 methods is much better than that of the NLS and SBI methods, especially if the problem size is large. In model order reduction problems of structures, the frequency response phase errors are typically ignored; however, it is shown in this work that phase errors may have significant effect on errors in time domain responses and should be considered as a performance measure as well. Considering the frequency response phase errors, the NLS and SBI methods have almost no phase errors. While the NLS and SBI methods may yield estimated models having good agreement with the FEM models in terms of frequency response, their time domain responses may be in error. For instance, while the NLS method had very good fit to frequency domain data in some examples, it gave unstable poles and/or nonminimum phase zeros; subsequently, the corresponding time domain response of these models estimated by the NLS method was not acceptable. To improve the accuracy in time domain responses of the SBI and NLS methods, static gain errors (i.e., gain error at zero frequency) of the estimated model play a crucial role. In addition, a very broad frequency range covering almost all the resonance frequencies of the system should be considered for estimating a model in frequency domain; however, this may yield very large model orders for the SBI and NLS methods. Moreover, for the SBI and NLS methods, small changes in pole and zero locations of estimated models do not yield large errors in frequency response fit; however, effects of these changes in time domain responses may be significant. Since the SBI and NLS methods are based on estimating pole and zero locations by using frequency domain data, their time domain responses are affected adversely by this sensitivity. On the other hand, CMS and QSM methods are not based on any pole and zero estimation technique; they use computational modes which represent the system behavior well both in frequency and time domain provided that sufficient number of modes are considered. Unlike the SBI and NLS methods, the CMS and QSM methods have zero static gain errors, which is responsible for yielding good time domain responses for slowly varying excitations. For excitations having high frequency contents, they yield acceptable results if the modes having the frequency range of at least twice the frequency range of interest in the problem are included. For the SBI and NLS methods in response to excitations having high frequency contents, generally the only remedy is to cover a very broad frequency range to have a good time domain response, which yields large estimated model orders.

402

Y. Cunedio˘glu et al. / Finite Elements in Analysis and Design 42 (2006) 367 – 403

Acknowledgements We hereby acknowledge the partial support of Secretary for Research Activities at Istanbul Technical University. References [1] J. Gu, Efficient model reduction methods for structural dynamics analyses, Doctoral Thesis, The University of Michigan, 2000. [2] M.A. Blair, T.S. Camino, J.M. Dickens, An iterative approach to a reduced mass matrix, in: Proceedings of the Ninth International Modal Analysis Conference, Florence, Italy, 1991, pp. 621–627. [3] N. Bouhaddi, R. Fillod, Substructuring using a linearized dynamic condensation method, Comput. Structures 45 (1992) 679–683. [4] N. Bouhaddi, R. Fillod, Model reduction by a simplified variant of dynamic condensation, J. Sound Vib. 191 (1996) 233–250. [5] M.I. Frishwell, S.D. Garvey, J.E.T. Penny, Model reduction using dynamic and iterated IRS techniques, J. Sound Vib. 186 (1995) 211–323. [6] R.J. Guyan, Reduction stiffness and mass matrices, AIAA J. 3 (1965) 380. [7] B.M. Irons, Structural eigenvalue problems: elimination of unwanted variables, AIAA J. 3 (1965) 961–962. [8] R.L. Kidder, Reduction of structural frequency equations, AIAA J. 11 (1973) 892. [9] A.Y.T. Leung, A simple dynamic substructure method, Earthquake Eng. Struct. Dyn. 16 (1988) 827–837. [10] J. O’Callahan, A procedure for an improved reduced system (IRS) model, in: Proceedings of the Seventh International Modal Analysis Conference, Las Vegas, NV, USA, 1989, pp. 17–21. [11] R.L. Bisplinghoff, H. Ashley, Aeroelasticity, Addison-Wesley, Reading, MA, 1995. [12] N.R. Maddox, On the number of modes necessary for accurate response and resulting forces in dynamic analysis, ASME J. Appl. Mech. 42 (1975) 516–517. [13] Z.D. Ma, I. Hagiwara, Improved mode-superposition technique for modal frequency response analysis of coupled acousticstructural systems, AIAA J. 29 (1991) 1720–1726. [14] Z.D. Ma, I. Hagiwara, Sensitivity analysis methods for coupled acoustic-structural system part I: modal sensitivities, AIAA J. 29 (1991) 1787–1795. [15] Z.D. Ma, I. Hagiwara, Sensitivity analysis methods for coupled acoustic-structural system part II: direct frequency response and its sensitivities, AIAA J. 29 (1991) 1796–1801. [16] W.H. Shyu, Quasi-static mode compensation for component mode synthesis of dynamical systems, Doctoral Thesis, The University of Michigan, 1996. [17] W.H. Shyu, Z.D. Ma, G.M. Hulbert, A new component mode synthesis method: quasi-static mode compensation, Finite Elements Anal. Des. 24 (1997) 271–281. [18] W.H. Shyu, J. Gu, G.M. Hulbert, Z.D. Ma, On the use of multiple quasi-static mode compensation sets for component mode synthesis of complex structures, Finite Elements Anal. Des. 35 (2000) 119–140. [19] J. Gu, Z.D. Ma, G.M. Hulbert, A new load dependent Ritz vector method for structural dynamic analyses: quasi-static Ritz vector method, Finite Elements Anal. Des. 36 (2000) 261–278. [20] K.J. Joo, E.L. Wilson, P. Leger, Ritz vectors and generation criteria for mode superposition analysis, Earthquake Eng. Struct. Dyn. 18 (1989) 149–167. [21] C. Lanczos, An iteration method for the solution of the eigenvalue problem of linear differential and integral operators, J. Res. Nat. Bur. Stand. 45 (1950) 252–282. [22] P. Leger, E.L. Wilson, R.W. Clough, The Use of Load Dependent Vectors for Dynamic and EarthquakeAnalyses, Earthquake Engineering Research Center Report, University of California, Berkeley, UCB/EERC-86/04, 1986. [23] P. Leger, Load dependent subspace reduction methods for structural dynamic computations, Comput. Structures 29 (1988) 993–999. [24] P. Leger, Application of load-dependent vectors bases for dynamic substructure analysis, AIAA J. 28 (1989) 177–179. [25] B. Nour-Omid, R.W. Clough, Dynamic analysis of structures using Lanczos coordinates, Earthquake Eng. Struct. Dyn. 12 (1984) 565–577.

Y. Cunedio˘glu et al. / Finite Elements in Analysis and Design 42 (2006) 367 – 403

403

[26] B. Nour-Omid, R.W. Clough, Block lanczos method for dynamic analysis of structures, Earthquake Eng. Struct. Dyn. 13 (1985) 271–275. [27] E.L. Wilson, M.W. Yuan, J.M. Dickens, Dynamic analysis by direct superposition of Ritz vectors, Earthquake Eng. Struct. Dyn. 10 (1982) 813–821. [28] E.L. Wilson, E.P. Bayo, Use of special Ritz vectors in dynamic substructure analysis, J. Struct. Eng. 112 (1986) 1944–1953. [29] J. Gu, Z.D. Ma, G.M. Hulbert, Quasi-static data recovery for dynamic analyses of structural systems, Finite Elements Anal. Des. 37 (2001) 825–841. [30] L. Ljung, System Identification, Prentice-Hall, Englewood Cliffs, NJ, 2000. [31] J. Schoukens, R. Pintelon, Identification of Linear Systems: A Practical Guideline for Accurate Modelling, Pergamon Press, London, 1991. [32] P. Van Overschee, B. De Moor, Subspace Identification for Linear Systems, Kluwer Academic Publishers, Boston, 1996. [33] D.S. Bayard, High order multivariable transfer function curve fitting: algorithms, sparse matrix methods and experimental results, Automatica 30 (1994) 1439–1444. [34] P. Verboven, P. Guillaume, B. Cauberghe, Multivariable Frequency-Response Curve Fitting with Application to Modal Parameter Estimation. [35] T. McKelvey, T. Abrahamsson, L. Ljung, Vibration data analysis for a commercial aircraft, Automatica 32 (1996) 1689–1700. [36] G.F. Franklin, J.D. Powell, A. Emami-Naeini, Feedback Control of Dynamic Systems, Addison-Wesley, Reading, MA, 1986. [37] J.S. Freudenberg, D.P. Looze, Frequency Domain Properties of Scalar and Multivariable Feedback Systems, Springer, Berlin, Heidelberg, 1988. [38] The Mathworks Inc., Matlab User Manual. [39] T. Kailath, Linear Systems, Prentice-Hall, Englewood Cliffs, NJ, 1980. [40] T. McKelvey, H. Akçay, L. Ljung, Subspace-based multivariable system identification from frequency response data, IEEE Trans. Autom. Control 41 (1996) 960–979. [41] T.J.R. Hughes, The Finite Element Method, Prentice-Hall, Englewood Cliffs, NJ, 1987.