# Experimental Studies on Finite Element Model Updating for a Heated Beam-Like Structure

Authors: Kaipeng Sun; Yonghui Zhao; Haiyan Hu
Journal: Shock and Vibration
 (2015)

---

## Abstract

An experimental study was made for the identification procedure of time-varying modal parameters and the finite element model updating technique of a beam-like thermal structure in both steady and unsteady high temperature environments. An improved time-varying autoregressive method was proposed first to extract the instantaneous natural frequencies of the structure in the unsteady high temperature environment. Based on the identified modal parameters, then, a finite element model for the structure was updated by using Kriging meta-model and optimization-based finite-element model updating method. The temperature-dependent parameters to be updated were expressed as low-order polynomials of temperature increase, and the finite element model updating problem was solved by updating several coefficients of the polynomials. The experimental results demonstrated the effectiveness of the time-varying modal parameter identification method and showed that the instantaneous natural frequencies of the updated model well tracked the trends of the measured values with high accuracy.


---

## Body

## 1. Introduction

Hypersonic vehicles are subject to very tough aerodynamic load and heating during their missions in the Earth’s atmosphere. The aerodynamic heating is extremely important because high temperature can affect the structural behavior in several detrimental ways [1]. The elevated temperature not only degrades the ability of structure materials to withstand loads, but also produces the thermal stresses in heated structures. In addition, the thermal stresses increase deformation, change buckling loads, and alter flutter behavior. For example, the studies [2–5] reported that even a minor change of temperature would result in a significant alteration in natural frequencies of a beam because of the thermal stresses when the beam was constrained. Thermal modal testing techniques can provide quantitative analysis for the effect of the thermal load [6]. Since the 1950s, many experiments have been done in NASA Langley and Dryden research centers for the metal and the composite panels [7, 8], the X-15 wing [9], and the X-34 FASTRAC Composite Rocket Nozzle [10]. The aforementioned studies have primarily focused on the dynamical properties of thermal structures in steady high temperature environments (SHTEs).As a matter of fact, thermal structures in unsteady high temperature environments (UHTEs), for example, hypersonic vehicles subjected to unsteady aerodynamic heating, have the characteristic of time-varying multiphysics fields. The identification of time-varying modal parameters is the forefront of inverse problems in structural dynamics. There have been numerous theoretical and experimental studies on the identification of time-varying modal parameters for engineering structures. The special issue [11] ofMechanical Systems and Signal Processing on the identification of time varying structures and systems in 2014 offered a survey of the field in its current state, reflected recent developments, and also pointed out into the future, in various aspects of the theories and applications. Although recent years have witnessed successful identifications of time-varying modal parameters of many engineering systems, such as vehicle-bridge systems [12, 13], machine condition monitoring systems [14], flexible manipulators [15], and civil structures [16, 17], their applications to thermal structures in temperature-varying environment are still not available. Yu et al. [18] proposed an undetermined blind source separation method to investigate the thermal effect on the modal parameters of a TC4 titanium-alloy column in a temperature-varying environment. They [19] also developed a time-varying modal parameter identification algorithm based on finite-data-window PAST and used it to investigate the effect of varying temperature and heating speed on the natural frequencies of a trapezoidal TA15 titanium-alloy plate. To the best knowledge of authors, considerably less attention has been paid to the finite element model updating (FEMU) for thermal structures, especially based on the identified time-varying modal parameters in UHTEs.Finite element (FE) modelling has received widespread acceptance and witnessed applications in various engineering disciplines. Hence, the FE model updating (FEMU) has become a useful tool to improve the modelling assumptions and parameters until a correlation between the analytical predictions and experimental results satisfies practical requirements. Mottershead and Friswell [20] comprehensively reviewed the model updating methods of structural models. In Marwala’s work [21], numerous computational intelligence techniques were introduced and applied to FE model updating with in-depth comparisons. As a most widely used method in FEMU, the iterative method has been formulated as an optimization problem and often based on sensitivity analysis [22] and computational intelligence techniques [23–25]. Generally, the optimization procedure is nonlinear and complex for a complicated FEMU problem. In recent years, the particle swarm optimization (PSO) technique has been developed to implement the optimization procedure.The structural FE models with many geometric and physical parameters to be updated may involve a large number of computations and need to be constructed by one of commercial finite element analysis packages, such as COMSOL, ANSYS, and NASTRAN. Hence, higher time consumption may be the disadvantage for the optimization-based algorithms, due to their iterative strategy and repeated analysis in simulation models during the optimization process. One way to overcome the difficulty of time consumption and FE package-related problems during the optimization-based model updating is to replace the FE model by an approximate surrogate/replacement meta-model that is fast-running and has fewer parameters involved. Simpson [26] made a comparison of response surface and Kriging models, which are the two commonly used meta-models, and drew a conclusion that both approximations predict reasonably well with the Kriging models having a slight overall advantage because of the lower root mean squared error values.The meta-model method for damage detection and reliability analysis has a long history. However, the Kriging meta-model (KMM) method for structural FEMU is somewhat new, especially for thermal structures in temperature-varying environment. As a continuation of authors’ work [27, 28], this paper presents KMM-PSO based FE model updating and selection based on the experimentally identified time-invariant and time-varying modal parameters of a thermal structure.The remainder is organized as follows. In Section2, the FE model updating and selection for a thermal structure is formulated, and two objective functions are given for the PSO. An overview of PSO and a brief introduction to the time-varying autoregressive method for output only identification are also presented. The Kriging meta-model is introduced to employ the PSO based FE model updating and selection. In Section 3, four groups of experiments are discussed for the modal parameters of a beam-like structure in room temperature environment, SHTE, and UHTE. In Section 4, the Kriging meta-model based FE model updating and selection procedure is carried out to update the software-based FE model. Finally, some conclusions are drawn in Section 5.


## 2. Formulation of FEMU for Thermal Structures

### 2.1. Problem Description

A linear thermal structure subject to unsteady heating can be described in terms of the distributed mass, damping, and stiffness matrices of the structure in time domain via the following differential equation:(1)Mx¨t+Ctx˙t+Ktxt=ft.
The real eigen-value problem for the jth mode reads(2)K(t)-ωj2Mφj=0,
where the detailed meanings of the physical quantities in the above two equations can be found in [28].The traditional FE model updating process is achieved by identifying the correct mass and stiffness matrices, which are generally time-invariant. This study, however, has to deal with a time-variant FE model updating process because the system matrices change over time when the temperature-dependent material and temperature-dependent boundary conditions are taken into consideration. Afterwards, using the measured data, the correct mass and stiffness matrices can be obtained by identifying the correct material parameters of the structure and the appropriate boundary conditions under different temperature conditions.For simplicity, the thermal structure of concern is a kind of slender or thin-walled structures, such as beams and plates, subject to uniform heating with a uniformly distributed temperature fieldT(x,y,z,t)≡T(t). In other words, the heat transfer process is neglected in the study. Let θ(t) denote the temperature increase in the thermal structure as follows:
(3)θ(t)=T(t)-Tref,
where Tref is the reference temperature. For the FEMU problem of the thermal structure subject to unsteady heating, the system parameters to be identified depend on the transient temperature and can be expressed as
(4)Γp=fl*θ=P0+P1θ+P2θ2+P3θ3+⋯,
where Γp is the parameter vector including the temperature-dependent material and boundary parameters to be corrected, and fl* is the general function of the temperature increase θ. In general, fl* can be expressed as a low-order polynomial, for example, a linear, quadratic, or cubic function, in which l represents the highest order of the polynomial. Then, the target parameters to be corrected change from temperature-dependent Γp to a constant parameter vector Pk, k=0,1,2,3.


### 2.2. Objective Function for FEMU at Reference Temperature

In a reference temperature environment (RTE), to correctly identify the moduli of elasticity that gives the updated FE model, the following objective function, which measures the distance between the measured modal data and the modal data predicted by FE model, should be minimized:(5)J0=∑j=1Nγjfjmeta-fjexp⁡fjexp⁡2,
where J0 is the error function or objective function, · is the Euclidean norm, γj is the weighted factor for the jth mode, and N is the number of measured modes, respectively. In (5), fjmeta is the jth instantaneous natural frequency obtained from the meta-model, and fjexp⁡ is the jth instantaneous natural frequency identified from the measured responses in a thermal-structural experiment based on the identification method for the time-invariant modal parameters.Thus, the process of FEMU may be viewed as an optimization problem as follows:(6)findP0,min⁡J0,s.t.P0l≤P0≤P0u,
where l and u represent the lower and upper bounds of the parameter coefficient vector P0, respectively.


### 2.3. Objective Function for FEMU in a UHTE

After an updated model in the RTE is achieved, the FEMU can be performed in a UHTE, which contextually means that the parameters to be corrected in this subsection become the constant parameter vectorPkk=0,1,2,3 in (4). It should be noted that P0 here represents several special parameters that cannot be identified in an RTE. The constant coefficient of thermal expansion, for example, has no effect on the stiffness matrix of the FE model and, thus, cannot be identified in an RTE, but must be updated in a UHTE.Furthermore, a new problem arises when the temperature-dependent parameter is approximated as a low-order polynomial. Because it is impossible to know a priori the exact order of polynomialfl* in (4), a parameter vector expressed by polynomials of different orders should be considered. In general, it is known that the higher the order is, the smaller the deviations between test and analysis are. However, the purpose of FEMU is to predict the structural response, to predict the effects of structural modifications, or to serve as a substructure model to be assembled as part of a model of the overall structure. From this viewpoint, for simplicity, one may use a low order expansion for the material properties. If the order is fixed to the maximal value, one would manually accept or reject the high-order terms by identifying whether the coefficient is close enough to zero. If any terms are neglected, the deviations between test and analysis should be reexamined. Therefore, the potential problems in the FEMU procedure are not only how to determine those parameters, but also how to select the updated model. To avoid establishing many models for FEMU and manually selecting the updated model, it is preferable to have an all-in-one procedure for both updating and selection. The PSO framework allows this simultaneous updating of all competing models and selection of the best model. Hence, the two aforementioned problems can be solved by minimising an integrated objective (or fitness) function.A number of fitness or objective functions have been available so far. Most previous studies have sought a model with the fewest updating parameters needed to produce FE model results that are closest to measured results. In this study, the Akaike information criterion (AIC) was used to represent the integrated objective function with an additional term to treat ill-conditioned and noisy systems. AIC can be described by the following equation:(7)J1=1Nsp∑j=1N∑k=1Nspγjfjmeta(tk)-fjexp⁡(tk)fjexp⁡(tk)2,(8)J2=Nlog⁡J1+λ2LP-P*2+κd,
with
(9)d=l1+l2+⋯+lNp,
where Nsp is the number of samplings, fjmeta(tk) and fjexp⁡(tk) are frequencies at time tk, and P* is the initial estimate of parameter P, respectively. The second term in the bracket represents regularization, where the parameter weighting matrix L should be chosen to reflect the uncertainty in the parameter estimation, and L=I represents the classical Tikhonov regularization. Link [29] suggested that the factor λ2 lies in the range of [0,0.3]. High λ values are used if there are many insensitive parameters and if the inverse problem is strongly ill-conditioned. The second term of the above equation is known as the model complexity penalty term, in which κ is a weighting factor, and lk(k=1,2,…,Np) represents the highest order of the kth polynomial. Np is the dimension of the parameter vector.Similarly to the FEMU in an RTE, the optimization process here may be expressed as follows:(10)findPk,k=0,1,2,3,min⁡J2,s.t.Pkl≤Pk≤Pku.Figure1 illustrates the detailed flow chart of the Kriging meta-model and the particle swarm optimization based FE model updating and selection. In the left dashed block, a Kriging meta-model is introduced to overcome the difficulty of time consumption because of the quantities of iterations during the updating process. In the right dashed block, modal parameters analyzed from Kriging meta-model and experimentally identified by time-varying autoregressive (TVAR) method are used to establish the objective function for the PSO based FE model updating and selection. Therefore, the next 3 subsections briefly introduce the stochastic optimization technique, the Kriging meta-modeling and the time-varying autoregressive method, and their relevance to the FE model updating and selection for thermal structures.Figure 1
Flow chart of Kriging meta-model (left dashed block) and particle swarm optimization based finite element model updating (right dashed block).


### 2.4. Stochastic Optimization Technique

In this study, PSO technique is employed to deal with the optimization problems described in (6) and (10). PSO is a population-based stochastic optimization technique developed by Eberhart and Kennedy [30] in 1995 and inspired by the social behavior of bird flocks and fish schools. The PSO procedure begins with a group of random particles and then searches for optima by updating generations. The PSO shares many similarities with evolutionary computation techniques, such as genetic algorithms (GAs) [24]. However, it has no evolution operators, such as crossover or mutation. In the PSO, particles update themselves with an internal velocity and also have memory and one-way information sharing mechanism.Similarly to GA, the algorithm begins with generating a group of random particles, called a swarm. At each iteration, the particles evaluate their fitness (positions relative to the goal) and share memories of their best positions with the swarm. Subsequently, each particle updates its velocity and position according to its best previous position (denoted bypbest) and that of the global best particle (denoted bygbest), which has far been found in the swarm. Let the position of a particle be denoted by xk∈Rn, and let vk be its velocity. They both are initially and randomly chosen and then iteratively updated according to two formulae. The following formula is used to update the particle’s velocity and position as determined by Shi and Eberhart [31]:
(11)vk=w·vk+c1·r1·pbestk-xk+c2·r2·gbest-xk,xk⟸xk+vk,
where w is an inertia coefficient that balances the global and local search, r1 and r2 are random numbers in the range [0,1] updated at each generation to prevent convergence on local optima, and c1 and c2 are the learning factors that control the influence of pbestk and gbest during the search process. Typically, c1 and c2 are set to be 2 for the sake of convergence [31].To avoid any physically unrealizable system matrix and the thermal buckling which easily occurs at extremely high temperatures, artificial position boundaries should be set for each particle. There are four types of boundaries, namely, absorbing, reflecting, invisible, and damping boundaries, as summarized in [32]. The damping boundary can provide a much robust and consistent optimization performance as compared with other boundary conditions and, thus, was used in this study. In addition to enforcing search-space boundaries after updating a particle’s position, it is also customary to impose limitations on the distance where a particle can move in a single step [33], which is done by limiting the velocity to a maximum value with the purpose of controlling the global exploration ability of the particle swarm and preventing the velocity from moving towards infinity.In the implementation of FE model updating and selection, the parametersPk(k=0,1,2,3) are usually set to be positions of particles for stochastic optimization technique. The above updating process should be repeated until a specified convergence value or total generation number is reached. This way, an optimal process for FE model updating and selection can be achieved.


### 2.5. Kriging Meta-Modeling

For completeness, a brief description of the Kriging meta-modeling is given in this subsection. Kriging was named after the pioneering work of D. G. Krige, a South African Mining Engineer, and was formally developed by Matheron [34]. Universal Kriging estimates the response Yx* at an untried site x* as the sum of a polynomial trend model fx*β, and a systematic departure term zx* representing low (large scale) and high frequency (small scale) variations around the trend model [35],
(12)Yx*=fx*β+zx*,
where fx*=f1x*f2x*⋯fMx* and β=β1β2⋯βMT are the regression model and the regression coefficients, respectively. In the right hand of (12), the first term is the mean value which can be thought as a globally valid trend function. And the second term zx* is a Gaussian distributed error term with zero mean and variance σ2. The covariance matrix of zx is given by
(13)Vzx=σ2Rns×ns.In (13), each element Rij of Rns×ns defined as Rxi,xj is the spatial correlation function between any two of the ns sample points xi and xj.In the current simulation, the termr represents the dimension of vector x and the linear regression model is chosen for the mean part of the Kriging function [36]. Accordingly fjx can be expressed as
(14)f1x=1,f2x=x1,f3x=x2,…,f1+rx=xr.The Gaussian correlation functionR is taken as
(15)Rxi,xj=∏k=1rexp⁡-pkxi(k)-xj(k)2,
where xik is the kth component of the ith sample point and pk is the unknown correlation parameter which needs to be fitted by optimization.To construct the Kriging model, the values of the regression coefficientsβ=β1β2⋯βMT must be approximated by using the generalized least squares theory at first. Then, the fitting correlation parameters pk can be quantified by using the maximum likelihood estimation. The generalized least-squares estimates of β and σ2, represented by β^ and σ^2, respectively, are given in detail as
(16)β^=FTR-1F-1FTR-1y,σ^2=y-Fβ^TR-1y-Fβ^ns,
where y and F are the ns×1 vector of output and the ns×ns matrix at the sample inputs, respectively, which are expressed as
(17)y=y1y2⋮yns,F=f1x1f2x1⋯fr+1x1f1x2f2x2⋯fr+1x2⋮⋮⋱⋮f1xnsf2xns⋯fr+1xns.With the vectorr=Rx*,x1Rx*,x2⋯Rx*,xnsT, the prediction at the unsampled location x* can be obtained as
(18)Y^x*=fx*β^+rTR-1y-Fβ^.For the FEMU problem of thermal structures, parametersΓp defined in (4) including the temperature-dependent material and boundary parameters to be corrected are taken as input parameters of Kriging meta-model and the output parameters are usually modal parameters, such as natural frequencies or modal shapes.Before the Kriging predictor is used in structural FE model updating, it should be verified to check whether the meta-model has enough accuracy. Sacks et al. stated that the cross-validation and integrate mean square error can be utilized to assess the accuracy of a Kriging model. The pointwise (local) estimate of actual error in Kriging approximation was given by computing the mean squared error (MSE)φ(x) as follows [36]:
(19)φ(x)=σ^21+uxTFTR-1Fu(x)-rxTR-1r(x),u(x)=FTR-1r(x)-f(x),
where σ^2 is the process variance defined in (16) and 1 is the vector of ones.


### 2.6. Brief Description of the TVAR Method

Identifying the time-varying modal parameters is an important issue in the FE model updating and selection for thermoelastic structures. Time-varying autoregressive method is one of the most popular time-frequency analysis methods for output only identification.This subsection deals with a TVAR processx-(t) (e.g., displacement, velocity, or acceleration) of order p in a discrete-time as the following:
(20)x-(t)=-∑i=1pal(t)x-(t-l)+e(t),
where e(t) is a stationary white noise process with zero mean and variance σ2, and {al(t),l=1,2,…,p} are the TVAR coefficients.Using the basis function expansion and regression approach, the TVAR processx-(t) of order p in a discrete-time can be expressed in matrix form as
(21)x-(t)=-XtTA+e(t),
where
(22)AT=a10,…,a1m,…,ap0,…,apm,XtT=x-(t-1)g0(t),…,x-(t-1)gm(t),…,x-(t-p)g0(t),…,x-(t-p)gm(t).e(t) is a stationary white noise process with zero mean and variance σ2; alk and m are the weighted coefficients and the dimension of the basis functions {gk(t),k=0,1,…,m}, respectively.The recursive least square (RLS) estimation and exponential forgetting method with a constant forgetting factor are used here such that the parameter estimation algorithm can be written as(23)A^N+1=A^N-PNXNλ+XNTPNXN-1·x-N+1+XNTA^N,PN+1=1λPN-PNXNλ+XNTPNXN-1XNTPN,
where the forgetting factor λ is chosen in the interval (0,1] and is typically close to one. The initial value of A^ and P can be selected as A^0=0, P0=μI, where μ≫1 and I is an identity matrix.Once the TVAR coefficients are obtained, the instantaneous natural frequencies can be derived from the conjugate rootssj(t), sj*(t) of the time-varying transfer function corresponding to the TVAR model as the following:
(24)fj(t)=12πΔtln⁡sj(t)·ln⁡sj*(t),
where Δt is the time-discretization step.


## 2.1. Problem Description

A linear thermal structure subject to unsteady heating can be described in terms of the distributed mass, damping, and stiffness matrices of the structure in time domain via the following differential equation:(1)Mx¨t+Ctx˙t+Ktxt=ft.
The real eigen-value problem for the jth mode reads(2)K(t)-ωj2Mφj=0,
where the detailed meanings of the physical quantities in the above two equations can be found in [28].The traditional FE model updating process is achieved by identifying the correct mass and stiffness matrices, which are generally time-invariant. This study, however, has to deal with a time-variant FE model updating process because the system matrices change over time when the temperature-dependent material and temperature-dependent boundary conditions are taken into consideration. Afterwards, using the measured data, the correct mass and stiffness matrices can be obtained by identifying the correct material parameters of the structure and the appropriate boundary conditions under different temperature conditions.For simplicity, the thermal structure of concern is a kind of slender or thin-walled structures, such as beams and plates, subject to uniform heating with a uniformly distributed temperature fieldT(x,y,z,t)≡T(t). In other words, the heat transfer process is neglected in the study. Let θ(t) denote the temperature increase in the thermal structure as follows:
(3)θ(t)=T(t)-Tref,
where Tref is the reference temperature. For the FEMU problem of the thermal structure subject to unsteady heating, the system parameters to be identified depend on the transient temperature and can be expressed as
(4)Γp=fl*θ=P0+P1θ+P2θ2+P3θ3+⋯,
where Γp is the parameter vector including the temperature-dependent material and boundary parameters to be corrected, and fl* is the general function of the temperature increase θ. In general, fl* can be expressed as a low-order polynomial, for example, a linear, quadratic, or cubic function, in which l represents the highest order of the polynomial. Then, the target parameters to be corrected change from temperature-dependent Γp to a constant parameter vector Pk, k=0,1,2,3.


## 2.2. Objective Function for FEMU at Reference Temperature

In a reference temperature environment (RTE), to correctly identify the moduli of elasticity that gives the updated FE model, the following objective function, which measures the distance between the measured modal data and the modal data predicted by FE model, should be minimized:(5)J0=∑j=1Nγjfjmeta-fjexp⁡fjexp⁡2,
where J0 is the error function or objective function, · is the Euclidean norm, γj is the weighted factor for the jth mode, and N is the number of measured modes, respectively. In (5), fjmeta is the jth instantaneous natural frequency obtained from the meta-model, and fjexp⁡ is the jth instantaneous natural frequency identified from the measured responses in a thermal-structural experiment based on the identification method for the time-invariant modal parameters.Thus, the process of FEMU may be viewed as an optimization problem as follows:(6)findP0,min⁡J0,s.t.P0l≤P0≤P0u,
where l and u represent the lower and upper bounds of the parameter coefficient vector P0, respectively.


## 2.3. Objective Function for FEMU in a UHTE

After an updated model in the RTE is achieved, the FEMU can be performed in a UHTE, which contextually means that the parameters to be corrected in this subsection become the constant parameter vectorPkk=0,1,2,3 in (4). It should be noted that P0 here represents several special parameters that cannot be identified in an RTE. The constant coefficient of thermal expansion, for example, has no effect on the stiffness matrix of the FE model and, thus, cannot be identified in an RTE, but must be updated in a UHTE.Furthermore, a new problem arises when the temperature-dependent parameter is approximated as a low-order polynomial. Because it is impossible to know a priori the exact order of polynomialfl* in (4), a parameter vector expressed by polynomials of different orders should be considered. In general, it is known that the higher the order is, the smaller the deviations between test and analysis are. However, the purpose of FEMU is to predict the structural response, to predict the effects of structural modifications, or to serve as a substructure model to be assembled as part of a model of the overall structure. From this viewpoint, for simplicity, one may use a low order expansion for the material properties. If the order is fixed to the maximal value, one would manually accept or reject the high-order terms by identifying whether the coefficient is close enough to zero. If any terms are neglected, the deviations between test and analysis should be reexamined. Therefore, the potential problems in the FEMU procedure are not only how to determine those parameters, but also how to select the updated model. To avoid establishing many models for FEMU and manually selecting the updated model, it is preferable to have an all-in-one procedure for both updating and selection. The PSO framework allows this simultaneous updating of all competing models and selection of the best model. Hence, the two aforementioned problems can be solved by minimising an integrated objective (or fitness) function.A number of fitness or objective functions have been available so far. Most previous studies have sought a model with the fewest updating parameters needed to produce FE model results that are closest to measured results. In this study, the Akaike information criterion (AIC) was used to represent the integrated objective function with an additional term to treat ill-conditioned and noisy systems. AIC can be described by the following equation:(7)J1=1Nsp∑j=1N∑k=1Nspγjfjmeta(tk)-fjexp⁡(tk)fjexp⁡(tk)2,(8)J2=Nlog⁡J1+λ2LP-P*2+κd,
with
(9)d=l1+l2+⋯+lNp,
where Nsp is the number of samplings, fjmeta(tk) and fjexp⁡(tk) are frequencies at time tk, and P* is the initial estimate of parameter P, respectively. The second term in the bracket represents regularization, where the parameter weighting matrix L should be chosen to reflect the uncertainty in the parameter estimation, and L=I represents the classical Tikhonov regularization. Link [29] suggested that the factor λ2 lies in the range of [0,0.3]. High λ values are used if there are many insensitive parameters and if the inverse problem is strongly ill-conditioned. The second term of the above equation is known as the model complexity penalty term, in which κ is a weighting factor, and lk(k=1,2,…,Np) represents the highest order of the kth polynomial. Np is the dimension of the parameter vector.Similarly to the FEMU in an RTE, the optimization process here may be expressed as follows:(10)findPk,k=0,1,2,3,min⁡J2,s.t.Pkl≤Pk≤Pku.Figure1 illustrates the detailed flow chart of the Kriging meta-model and the particle swarm optimization based FE model updating and selection. In the left dashed block, a Kriging meta-model is introduced to overcome the difficulty of time consumption because of the quantities of iterations during the updating process. In the right dashed block, modal parameters analyzed from Kriging meta-model and experimentally identified by time-varying autoregressive (TVAR) method are used to establish the objective function for the PSO based FE model updating and selection. Therefore, the next 3 subsections briefly introduce the stochastic optimization technique, the Kriging meta-modeling and the time-varying autoregressive method, and their relevance to the FE model updating and selection for thermal structures.Figure 1
Flow chart of Kriging meta-model (left dashed block) and particle swarm optimization based finite element model updating (right dashed block).


## 2.4. Stochastic Optimization Technique

In this study, PSO technique is employed to deal with the optimization problems described in (6) and (10). PSO is a population-based stochastic optimization technique developed by Eberhart and Kennedy [30] in 1995 and inspired by the social behavior of bird flocks and fish schools. The PSO procedure begins with a group of random particles and then searches for optima by updating generations. The PSO shares many similarities with evolutionary computation techniques, such as genetic algorithms (GAs) [24]. However, it has no evolution operators, such as crossover or mutation. In the PSO, particles update themselves with an internal velocity and also have memory and one-way information sharing mechanism.Similarly to GA, the algorithm begins with generating a group of random particles, called a swarm. At each iteration, the particles evaluate their fitness (positions relative to the goal) and share memories of their best positions with the swarm. Subsequently, each particle updates its velocity and position according to its best previous position (denoted bypbest) and that of the global best particle (denoted bygbest), which has far been found in the swarm. Let the position of a particle be denoted by xk∈Rn, and let vk be its velocity. They both are initially and randomly chosen and then iteratively updated according to two formulae. The following formula is used to update the particle’s velocity and position as determined by Shi and Eberhart [31]:
(11)vk=w·vk+c1·r1·pbestk-xk+c2·r2·gbest-xk,xk⟸xk+vk,
where w is an inertia coefficient that balances the global and local search, r1 and r2 are random numbers in the range [0,1] updated at each generation to prevent convergence on local optima, and c1 and c2 are the learning factors that control the influence of pbestk and gbest during the search process. Typically, c1 and c2 are set to be 2 for the sake of convergence [31].To avoid any physically unrealizable system matrix and the thermal buckling which easily occurs at extremely high temperatures, artificial position boundaries should be set for each particle. There are four types of boundaries, namely, absorbing, reflecting, invisible, and damping boundaries, as summarized in [32]. The damping boundary can provide a much robust and consistent optimization performance as compared with other boundary conditions and, thus, was used in this study. In addition to enforcing search-space boundaries after updating a particle’s position, it is also customary to impose limitations on the distance where a particle can move in a single step [33], which is done by limiting the velocity to a maximum value with the purpose of controlling the global exploration ability of the particle swarm and preventing the velocity from moving towards infinity.In the implementation of FE model updating and selection, the parametersPk(k=0,1,2,3) are usually set to be positions of particles for stochastic optimization technique. The above updating process should be repeated until a specified convergence value or total generation number is reached. This way, an optimal process for FE model updating and selection can be achieved.


## 2.5. Kriging Meta-Modeling

For completeness, a brief description of the Kriging meta-modeling is given in this subsection. Kriging was named after the pioneering work of D. G. Krige, a South African Mining Engineer, and was formally developed by Matheron [34]. Universal Kriging estimates the response Yx* at an untried site x* as the sum of a polynomial trend model fx*β, and a systematic departure term zx* representing low (large scale) and high frequency (small scale) variations around the trend model [35],
(12)Yx*=fx*β+zx*,
where fx*=f1x*f2x*⋯fMx* and β=β1β2⋯βMT are the regression model and the regression coefficients, respectively. In the right hand of (12), the first term is the mean value which can be thought as a globally valid trend function. And the second term zx* is a Gaussian distributed error term with zero mean and variance σ2. The covariance matrix of zx is given by
(13)Vzx=σ2Rns×ns.In (13), each element Rij of Rns×ns defined as Rxi,xj is the spatial correlation function between any two of the ns sample points xi and xj.In the current simulation, the termr represents the dimension of vector x and the linear regression model is chosen for the mean part of the Kriging function [36]. Accordingly fjx can be expressed as
(14)f1x=1,f2x=x1,f3x=x2,…,f1+rx=xr.The Gaussian correlation functionR is taken as
(15)Rxi,xj=∏k=1rexp⁡-pkxi(k)-xj(k)2,
where xik is the kth component of the ith sample point and pk is the unknown correlation parameter which needs to be fitted by optimization.To construct the Kriging model, the values of the regression coefficientsβ=β1β2⋯βMT must be approximated by using the generalized least squares theory at first. Then, the fitting correlation parameters pk can be quantified by using the maximum likelihood estimation. The generalized least-squares estimates of β and σ2, represented by β^ and σ^2, respectively, are given in detail as
(16)β^=FTR-1F-1FTR-1y,σ^2=y-Fβ^TR-1y-Fβ^ns,
where y and F are the ns×1 vector of output and the ns×ns matrix at the sample inputs, respectively, which are expressed as
(17)y=y1y2⋮yns,F=f1x1f2x1⋯fr+1x1f1x2f2x2⋯fr+1x2⋮⋮⋱⋮f1xnsf2xns⋯fr+1xns.With the vectorr=Rx*,x1Rx*,x2⋯Rx*,xnsT, the prediction at the unsampled location x* can be obtained as
(18)Y^x*=fx*β^+rTR-1y-Fβ^.For the FEMU problem of thermal structures, parametersΓp defined in (4) including the temperature-dependent material and boundary parameters to be corrected are taken as input parameters of Kriging meta-model and the output parameters are usually modal parameters, such as natural frequencies or modal shapes.Before the Kriging predictor is used in structural FE model updating, it should be verified to check whether the meta-model has enough accuracy. Sacks et al. stated that the cross-validation and integrate mean square error can be utilized to assess the accuracy of a Kriging model. The pointwise (local) estimate of actual error in Kriging approximation was given by computing the mean squared error (MSE)φ(x) as follows [36]:
(19)φ(x)=σ^21+uxTFTR-1Fu(x)-rxTR-1r(x),u(x)=FTR-1r(x)-f(x),
where σ^2 is the process variance defined in (16) and 1 is the vector of ones.


## 2.6. Brief Description of the TVAR Method

Identifying the time-varying modal parameters is an important issue in the FE model updating and selection for thermoelastic structures. Time-varying autoregressive method is one of the most popular time-frequency analysis methods for output only identification.This subsection deals with a TVAR processx-(t) (e.g., displacement, velocity, or acceleration) of order p in a discrete-time as the following:
(20)x-(t)=-∑i=1pal(t)x-(t-l)+e(t),
where e(t) is a stationary white noise process with zero mean and variance σ2, and {al(t),l=1,2,…,p} are the TVAR coefficients.Using the basis function expansion and regression approach, the TVAR processx-(t) of order p in a discrete-time can be expressed in matrix form as
(21)x-(t)=-XtTA+e(t),
where
(22)AT=a10,…,a1m,…,ap0,…,apm,XtT=x-(t-1)g0(t),…,x-(t-1)gm(t),…,x-(t-p)g0(t),…,x-(t-p)gm(t).e(t) is a stationary white noise process with zero mean and variance σ2; alk and m are the weighted coefficients and the dimension of the basis functions {gk(t),k=0,1,…,m}, respectively.The recursive least square (RLS) estimation and exponential forgetting method with a constant forgetting factor are used here such that the parameter estimation algorithm can be written as(23)A^N+1=A^N-PNXNλ+XNTPNXN-1·x-N+1+XNTA^N,PN+1=1λPN-PNXNλ+XNTPNXN-1XNTPN,
where the forgetting factor λ is chosen in the interval (0,1] and is typically close to one. The initial value of A^ and P can be selected as A^0=0, P0=μI, where μ≫1 and I is an identity matrix.Once the TVAR coefficients are obtained, the instantaneous natural frequencies can be derived from the conjugate rootssj(t), sj*(t) of the time-varying transfer function corresponding to the TVAR model as the following:
(24)fj(t)=12πΔtln⁡sj(t)·ln⁡sj*(t),
where Δt is the time-discretization step.


## 3. Experimental Studies on a Cantilever Beam

In this study, experimental modal analysis and operational modal analysis were carried out to obtain the time-invariant and time-varying modal parameters, respectively, in different temperature environments. The experimental object here is a cantilever beam made of aluminum installed in a movable box-type resistance furnace, as shown as (a) and (b) in Figure2. The beam was dynamically driven by a hammer impact or a vibration shaker excitation in experiments. A double-lug-type connector was used to connect the beam and the vibration shaker when the beam was heated by the furnace where only the vibration shaker provided feasible excitations. It should be emphasized that the mass of the connector should not be neglected in modal analysis. For simplicity, hence, the terms “beam A” and “beam B” are used hereinafter for the cantilever beam without the double-lug-type connector and with the double-lug-type connector, respectively.Experimental setup for thermal test of a beam: (a) the aluminum beam, (b) the experimental setup, and (c) the schematic framework.
(a)
(b)
(c)Figure2(c) shows the schematic framework of experimental setup where beam B is excited by a vibration shaker. Table 1 lists four groups of experiments for different cases. For all groups, the velocity responses of the beam were measured by using a laser vibrometer as a noncontacted measurement technique.Table 1
Experiment descriptions.
Number of group
Beam type
Excitation type
Temperature environment
1
Beam A
Force-hammer
RTE
2
Beam B
Force-hammer
RTE
3
Beam B
Vibration shaker
SHTE
4
Beam B
Vibration shaker
UHTE


### 3.1. The First and Second Groups of Experiments

In the first and second groups of experiments, the frequency responses of beam A and beam B were measured via a hammer impact at room temperature, respectively. Without loss of generality, the room temperature was assumed to be reference temperature. Figure3 illustrates the amplitude-frequency responses (AFRs) of the measured frequency response functions (FRFs) of the two beams. The first natural frequencies of beam A and beam B were 7.3125 Hz and 7.125 Hz, respectively, while the second ones were 45.875 Hz and 32 Hz, respectively. The figure clearly shows that the added mass of the double-lug-type connector greatly reduced the second natural frequency of beam B but had a small influence on the first natural frequency of beam B due to the attachment position.Figure 3
Amplitude-frequency responses of beam A and beam B driven by a hammer impact.


### 3.2. The Third Group of Experiments

In the third group of experiments, the beam B was subject to the heating of furnace. Hence, it was driven by a vibration exciter outside of the furnace through a long steel rod connecting the beam with a double-lug-type connector. The heating of furnace was controlled by the temperature controlled tank with two windows displaying two temperatures, that is, the target temperature and the cavity temperature, respectively. Besides, the actual temperature of the beam at different time instants was measured by a K-type thermocouple thermometer as shown in Figure2. The beam was heated in steady environments of high temperature at 200°C, 300°C, 400°C, and 500°C, respectively. For all the experiments in Sections 3.2 and 3.3, the random excitation was provided by the shaker and the sampling frequency was set at 512 Hz. Figure 4 illustrates the AFRs produced from velocity responses of beam B at different temperatures and Table 2 lists the first two natural frequencies. They demonstrate that the first two natural frequencies of the beam decreased with an increase of the temperature.Table 2
The first two natural frequencies (Hz) of beam B at different temperatures.
Mode
0°C
200°C
300°C
400°C
500°C
1
10.25
9
8
7.5
7.25
2
26
22.75
22.5
21.75
20Figure 4
Amplitude-frequency responses of beam B driven by a vibration shaker at different temperatures.


### 3.3. The Fourth Group of Experiments

In this group of experiments, beam B was heated in an unsteady high temperature environment. The temperature was increased from the room temperature to about 500°C. At the same time of temperature increment, the beam was subject to a random force from the vibration shaker and the velocity responses of the beam were measured by using a laser vibrometer. The measured responses were then used to extract the time-varying modal parameters via the continuous wavelet transform (CWT) method and the TVAR method.Figure5 shows the CWT scalogram of the velocity response, using the Complex Morlet 33 as the wavelet basis. The top subfigure gives the signal waveform of the response, the left bottom subfigure is the corresponding power spectrum, and the right bottom subfigure is the time-frequency analysis result with a color bar indicating the magnitude levels on the right. Figure 6 illustrates the first two instantaneous natural frequencies, identified by using the TVAR algorithm and labeled in the left Y-axis, with respect to the measured temperature on the beam labeled in the right Y-axis.Figure 5
Complex Morlet transform scalogram of the velocity response.Figure 6
The first two instantaneous natural frequencies of beam B and transient temperature.With the comparison of Figures5 and 6, both CWT method and TVAR method provided the good time-frequency representation of nonstationary dynamics, but the latter gave the result of much higher time-frequency resolution. In addition, the TVAR method could provide parametric results, which can be directly used in the next FEMU procedure.


## 3.1. The First and Second Groups of Experiments

In the first and second groups of experiments, the frequency responses of beam A and beam B were measured via a hammer impact at room temperature, respectively. Without loss of generality, the room temperature was assumed to be reference temperature. Figure3 illustrates the amplitude-frequency responses (AFRs) of the measured frequency response functions (FRFs) of the two beams. The first natural frequencies of beam A and beam B were 7.3125 Hz and 7.125 Hz, respectively, while the second ones were 45.875 Hz and 32 Hz, respectively. The figure clearly shows that the added mass of the double-lug-type connector greatly reduced the second natural frequency of beam B but had a small influence on the first natural frequency of beam B due to the attachment position.Figure 3
Amplitude-frequency responses of beam A and beam B driven by a hammer impact.


## 3.2. The Third Group of Experiments

In the third group of experiments, the beam B was subject to the heating of furnace. Hence, it was driven by a vibration exciter outside of the furnace through a long steel rod connecting the beam with a double-lug-type connector. The heating of furnace was controlled by the temperature controlled tank with two windows displaying two temperatures, that is, the target temperature and the cavity temperature, respectively. Besides, the actual temperature of the beam at different time instants was measured by a K-type thermocouple thermometer as shown in Figure2. The beam was heated in steady environments of high temperature at 200°C, 300°C, 400°C, and 500°C, respectively. For all the experiments in Sections 3.2 and 3.3, the random excitation was provided by the shaker and the sampling frequency was set at 512 Hz. Figure 4 illustrates the AFRs produced from velocity responses of beam B at different temperatures and Table 2 lists the first two natural frequencies. They demonstrate that the first two natural frequencies of the beam decreased with an increase of the temperature.Table 2
The first two natural frequencies (Hz) of beam B at different temperatures.
Mode
0°C
200°C
300°C
400°C
500°C
1
10.25
9
8
7.5
7.25
2
26
22.75
22.5
21.75
20Figure 4
Amplitude-frequency responses of beam B driven by a vibration shaker at different temperatures.


## 3.3. The Fourth Group of Experiments

In this group of experiments, beam B was heated in an unsteady high temperature environment. The temperature was increased from the room temperature to about 500°C. At the same time of temperature increment, the beam was subject to a random force from the vibration shaker and the velocity responses of the beam were measured by using a laser vibrometer. The measured responses were then used to extract the time-varying modal parameters via the continuous wavelet transform (CWT) method and the TVAR method.Figure5 shows the CWT scalogram of the velocity response, using the Complex Morlet 33 as the wavelet basis. The top subfigure gives the signal waveform of the response, the left bottom subfigure is the corresponding power spectrum, and the right bottom subfigure is the time-frequency analysis result with a color bar indicating the magnitude levels on the right. Figure 6 illustrates the first two instantaneous natural frequencies, identified by using the TVAR algorithm and labeled in the left Y-axis, with respect to the measured temperature on the beam labeled in the right Y-axis.Figure 5
Complex Morlet transform scalogram of the velocity response.Figure 6
The first two instantaneous natural frequencies of beam B and transient temperature.With the comparison of Figures5 and 6, both CWT method and TVAR method provided the good time-frequency representation of nonstationary dynamics, but the latter gave the result of much higher time-frequency resolution. In addition, the TVAR method could provide parametric results, which can be directly used in the next FEMU procedure.


## 4. Meta-Model Based FEMU

### 4.1. Numerical Simulation and Kriging-Based Meta-Modeling

In this study, COMSOL, a software of multiphysics, was used for the FE based modal analysis of the beam under various conditions of parameter combinations. Figure7 shows the dimension chart of the beam. As illustrated in Figure 8, the geometry model of beam B built in COMSOL contains two parts, that is, the beam and the double-lug-type connector. In the numerical simulations, the connecting stiffness of the bolt joints was modeled by attaching an auxiliary surface between the assembled parts, defining the material properties of the auxiliary surface and connecting the assembled parts with the multipoint constraint (MPC) strategy.Figure 7
The dimension chart of the beam.Figure 8
The geometry model of beam B.As well known, the numerical simulation may give a good prediction for the natural frequencies of the beam, but the results are not parameterized. Furthermore, the simulation process may be time-consuming. The main idea of meta-modeling is to construct a parameterized mathematical model between the input parameters and the output results by a number of numerical simulations and then use the model to predict other output results. In this study, the Kriging method was employed to construct the meta-model for building accurate global approximation in a given design space.In this study, several parameters, such as the densityρb of the beam, the density ρc of the connector, the added mass ma of the long steel rod, the density ρas, and elastic modulus Eas of the auxiliary surface, were assumed to be temperature-independent, while the elastic modulus Eb of the beam material and the added stiffness Ka of the long rod were taken as temperature-dependent parameters. As mentioned earlier, the temperature-dependent parameters can be expressed as low-order polynomials; that is, Eb and Kr yield the following polynomials of temperature increase θ:
(25)Eb(θ)=Eb0+Eb1θ+Eb2θ2,Ka(θ)=Ka0+Ka1θ+Ka2θ2,
where Eb0 and Ka0 are the elastic modulus and the stiffness of the long rod at the reference temperature Tref, respectively, and Ebk and Kak(k=1,2) are the coefficients independent of temperature. In the study, ρb was taken as 2252.75 kg/m3 according to measured mass and volume of the beam, and ρc of the connector was taken as 7949.6 kg/m3 in the same way. The value of ρas was taken as a constant of 2000 kg/m3 since the thin auxiliary surface did not have any significant influence on the modal parameter. The parameters to be updated in the next 3 subsections are (1) Eb0 and Eas, (2) ma and Ka0, and (3) Ebk and Kak(k=1,2).


### 4.2. First Step of FEMU

For the first step of FEMU of beam B under a hammer impact at room temperature,Eb0 and Eas were used as updating parameters of the KMM based FEMU. To simulate the initial model of beam B to be updated, the initial values of Eb0 and Eas were taken as 65 GPa and 9.5 × 104 Pa. The modal analysis via FEM method was performed on the initial model to obtain the initial natural frequencies. The initial values of the first two natural frequencies and the corresponding differences are shown in Table 3.Table 3
Natural frequency differences of beam B under a hammer impact before and after model updating.
Mode
Measured (Hz)
Initial (Hz)
Error (%)
Updated (Hz)
Error (%)
1
7.125
7.480
4.98
7.125
0.00
2
32.000
33.556
4.86
32.000
0.00The updating parameters share the same region for the training data in the Kriging meta-modeling. To construct the Kriging meta-model valid over a range of parameters, the moduli of elasticityEb0 and Eas were restricted to vary from 50 to 70 GPa and 8 × 104 to 1 × 105 Pa, respectively. The design of experiment (DOE) is a key problem in deciding how to select the inputs, at which the deterministic computer codes are run in order to most efficiently control or reduce the statistical uncertainty of the computed predictions. In this study, the rectangular grid method [36] was used to deal with the DOE problem. This was easily done by using the full multiparameter sweep [37] and specifying all combinations type in COMSOL. A total of 25 experiments were carried out. The sampled parameter values and corresponding natural frequencies computed from FE models were used as the training data of the Kriging meta-model. A 50 × 50 uniform mesh grid in the region covered by the design sites was generated to evaluate the predictor. Figures 9 and 10 illustrate the mesh plots of the predicted values of the first and second natural frequencies at the grid points, respectively. The horizontal axes are parameters selected while the vertical axis gives the predicted response (natural frequency) at any point or location.Figure 9
Predicted values of the first natural frequency.Figure 10
Predicted values of the second natural frequency.To check the accuracy of the Kriging meta-model, the MSEs were computed for each mode as shown in Figures11 and 12. They demonstrate that all MSE values were close to zero so that the created meta-model had a high regression accuracy.Figure 11
Mean squared errors of the first natural frequency.Figure 12
Mean squared errors of the second natural frequency.Then, the FEMU was performed with the FE model replaced by the Kriging meta-model. The residuals between the accurate (predicted by Kriging meta-model) and the measured (by experiment) natural frequencies were used in the optimized objective function expressed in (5). A single-objective optimization algorithm with equal weight for each natural frequency was implemented to achieve the best minimization of natural frequency residuals. The optimization algorithm used to minimize the objective function is an improved PSO method in MATLAB. In implementing the PSO for the FE model updating, the population was taken as 50, c1 and c2 were set to be 2, and w was set to be 1. The tuning minimization process was over when the tolerances were achieved or a predefined number of iterations were reached. Table 3 shows the updated natural frequencies and their differences of beam B and illustrates good results of PSO based FE model updating. The final updated results for the parameter were Eb0 = 58.78 GPa and Eas = 9.524 × 104 Pa. As shown in Table 3, the errors between the first two natural frequencies measured and those obtained from the initial FE model were about 5% on average. When the KMM-PSO based FE model updating algorithm was used, the errors were reduced to very small values.


### 4.3. Second Step of FEMU

For the second step of FEMU of beam B under vibration shaker excitation at room temperature,ma and Ka0 were used as updating parameters of the KMM-PSO based FEMU. The numerical modal analysis, the corresponding meta-modeling, and the sequential FEMU procedure are almost the same as the first step of FEMU in Section 4.2, and in consequence the detailed processes are not repeated in this subsection. The initial values of ma and Ka0 were set as 0.2 kg and 2 × 103 N/m, and the corresponding updated values were 0.14548 kg and 2098.12 N/m, respectively. Table 4 shows the initial and updated values of the first two natural frequencies and the corresponding differences. As listed in the table, the errors between the first two natural frequencies measured and those obtained from the initial FE model were about 5.8% on average. When the KMM-PSO based FE model updating algorithm was used, the errors were reduced to tiny values.Table 4
Natural frequency differences of beam B under vibration shaker excitation before and after model updating.
Mode
Measured (Hz)
Initial (Hz)
Error (%)
Updated (Hz)
Error (%)
1
10.250
9.884
3.57
10.250
0.00
2
26.000
23.898
8.08
25.996
0.02


### 4.4. Third Step of FEMU

For the third step of FEMU of beam B under vibration shaker excitation in a UHTE, the parameters to be updated are temperature-dependent.Ebk and Kak(k=1,2) were used as updating parameters of the KMM-PSO based FEMU. The time-varying modal parameters identified in Section 3.3 were used to establish the objective function in (8). In addition, the linear functions of temperature increase for the parameters were used as initial guesses and are shown in Figures 13 and 14, respectively. Figure 6 also illustrates the first two instantaneous natural frequencies computed via the FE method with these initial parameters.Figure 13
Temperature-dependent elastic modulus of the beam material.Figure 14
Temperature-dependent stiffness of the long rod.As mentioned earlier, the temperature-dependent parameters to be identified can be expressed as lower-order polynomials of the temperature increase, but the exact order is unknown beforehand. Hence, the best order and the coefficients of these polynomials should be simultaneously identified. Without loss of generality, a simple case was considered to verify the proposed method. For this purpose,Eb and Ka were assumed to be linear or quadratic functions of the temperature increase. The process of ensuring the positivity of Young’s modulus, for instance, can be recalled as follows. By appropriately selecting the maximal and minimal values of E, ΔEk≜γEmax⁡-Emin⁡/Tmax⁡k, k=1,2, was defined, where k represents the order and γ is a weighting coefficient. In practice, γ≈1 was taken for a low order and γ≈0.1 was set for a high order to prevent strong nonlinearity of material property. Afterwards, the bounds of E with an increase of temperature were checked by using the Monte Carlo method.Different from the typical FEMU procedure, the method proposed here includes FE model selection. For instance, the beams were modelled by four competing models,Mk, k=1,…,4, as listed in Table 5. Here, each particle had 4 dimensions such that all the competing models should be searched in a 4-dimensional space, and then all of the particles were built to compute the fitness function J2 in (8) so as to find the smallest value of fitness function. Although all the models should be searched in the same space, each model was actually constrained to a particular subspace of the space. Afterwards, the most important issue for the FEMU in time domain is the implementation of the PSO algorithm. It was assumed that the values of the polynomial coefficients for the thermal-structural properties were restricted to different intervals and varied in their corresponding intervals. Table 6 presents the initial values of the parameters and the lower and upper bounds of their intervals. The factor λ2 and weighting matrix L in the objective function J2 were set as λ2=5×10-4, κ=8×10-4, and L=diag([1,3]), respectively. After substituting the objective function J1 with equal weight for each natural frequency into J2, the multiobjective optimization algorithm was implemented to achieve the best minimization of J2.Table 5
Model parameterization.
Model identity
Number of parameters
Max. order(Eb,Ka)
Parameter symbols
M
1
2
(1, 1)
E
b
1
,
K
a
1
M
2
3
(1, 2)
E
b
1
,
K
a
1
,
K
a
2
M
3
3
(2, 1)
E
b
1
,
E
b
2
,
K
a
1
M
4
4
(2, 2)
E
b
1
,
E
b
2
,
K
a
1
,
K
a
2Table 6
Particle swarm optimization parameters.
PSO parameter
Initial
Lower bound
Upper bound
Identified
E
b
1 (GPa·°C−1)
−0.054
-
0.2
0
-
4.2463
e
-
2
E
b
2 (GPa·°C−2)
0
-
9
e
-
5
9
e
-
5
-
2.0803
e
-
5
K
a
1 (N·m−1·°C−1)
−0.002
-
5
e
-
2
0
-
2.1041
e
-
3
K
a
2 (N·m−1·°C−2)
0
-
9
e
-
5
9
e
-
5
-
8.9016
e
-
5Figure13 shows the identified values of Young’s modulus Eb of the beam compared with the initial values. Figure 14 shows the identified values of the added stiffness Ka of the long rod compared with the initial values.Figure15 shows the convergence of the objective function J2 over the 50 iterations of the algorithm and illustrates that the PSO algorithm rapidly converged to the ultimate minimum error within the first 22 iterations. Figure 15 also illustrates the convergence behavior of the global best model over 50 iterations. It is quite clear from Figure 15 that the objective function played a significant role in updating the model parameters. The global best model in the simulation began with M1, then changed to M3 and M1, afterwards went to M3, and subsequently remained unchanged for the rest of the simulation. The result indicates that M3 was the global best model. That is, Eb and Ka were a quadratic polynomial and a linear polynomial of the temperature increase, respectively, to balance the model errors and complexity.Figure 15
Global fitness and global best model number during updating.To have a quantitative discussion, the mean absolute percentage error (MAPE) is defined as(26)MAPE=1Nsp∑k=1Nspyk-y^kyk×100%,
where Nsp is the total number of samplings and yk and y^k denote the true values computed by the direct modal analysis and the identified value at the kth time instant, respectively. Figure 16 illustrates the MAPEs of the first two natural frequencies. The average error between the initial values of the first two natural frequencies and the true values was 1.8%. When the KMM-PSO based FE model updating and selection method were used, the error was reduced to 0.1% on average.Figure 16
Mean absolute percentage errors of the first two natural frequencies.Overall, the KMM-PSO based FEMU approach proposed in this study updates the model of a thermal structure in both RTE and UHTE well. On the other hand, the proposed method can greatly reduce the computation time compared to other FEMU algorithms based on direct FE modeling and modal analysis by using commercial FE analysis packages. Table7 gives the comparison of computational time of the proposed method with the FEMU algorithm based on direct modal analyses in COMSOL. It should be pointed out that the FEMU algorithm based on direct modal analysis in COMSOL was not actually carried out and the corresponding computation time was predicted according to the calling numbers of Kriging predictor function and the actual time (about 5 s) of one run of FE modal analysis in COMSOL.Table 7
Comparison of computational time.
Method
Number of MA
Time of MA
Total time
The proposed method
10646400
3785 s
3925 s
FEMU based on direct MA in COMSOL
10646400
616 d
616 d
MA: modal analysis; s: seconds; d: days.


## 4.1. Numerical Simulation and Kriging-Based Meta-Modeling

In this study, COMSOL, a software of multiphysics, was used for the FE based modal analysis of the beam under various conditions of parameter combinations. Figure7 shows the dimension chart of the beam. As illustrated in Figure 8, the geometry model of beam B built in COMSOL contains two parts, that is, the beam and the double-lug-type connector. In the numerical simulations, the connecting stiffness of the bolt joints was modeled by attaching an auxiliary surface between the assembled parts, defining the material properties of the auxiliary surface and connecting the assembled parts with the multipoint constraint (MPC) strategy.Figure 7
The dimension chart of the beam.Figure 8
The geometry model of beam B.As well known, the numerical simulation may give a good prediction for the natural frequencies of the beam, but the results are not parameterized. Furthermore, the simulation process may be time-consuming. The main idea of meta-modeling is to construct a parameterized mathematical model between the input parameters and the output results by a number of numerical simulations and then use the model to predict other output results. In this study, the Kriging method was employed to construct the meta-model for building accurate global approximation in a given design space.In this study, several parameters, such as the densityρb of the beam, the density ρc of the connector, the added mass ma of the long steel rod, the density ρas, and elastic modulus Eas of the auxiliary surface, were assumed to be temperature-independent, while the elastic modulus Eb of the beam material and the added stiffness Ka of the long rod were taken as temperature-dependent parameters. As mentioned earlier, the temperature-dependent parameters can be expressed as low-order polynomials; that is, Eb and Kr yield the following polynomials of temperature increase θ:
(25)Eb(θ)=Eb0+Eb1θ+Eb2θ2,Ka(θ)=Ka0+Ka1θ+Ka2θ2,
where Eb0 and Ka0 are the elastic modulus and the stiffness of the long rod at the reference temperature Tref, respectively, and Ebk and Kak(k=1,2) are the coefficients independent of temperature. In the study, ρb was taken as 2252.75 kg/m3 according to measured mass and volume of the beam, and ρc of the connector was taken as 7949.6 kg/m3 in the same way. The value of ρas was taken as a constant of 2000 kg/m3 since the thin auxiliary surface did not have any significant influence on the modal parameter. The parameters to be updated in the next 3 subsections are (1) Eb0 and Eas, (2) ma and Ka0, and (3) Ebk and Kak(k=1,2).


## 4.2. First Step of FEMU

For the first step of FEMU of beam B under a hammer impact at room temperature,Eb0 and Eas were used as updating parameters of the KMM based FEMU. To simulate the initial model of beam B to be updated, the initial values of Eb0 and Eas were taken as 65 GPa and 9.5 × 104 Pa. The modal analysis via FEM method was performed on the initial model to obtain the initial natural frequencies. The initial values of the first two natural frequencies and the corresponding differences are shown in Table 3.Table 3
Natural frequency differences of beam B under a hammer impact before and after model updating.
Mode
Measured (Hz)
Initial (Hz)
Error (%)
Updated (Hz)
Error (%)
1
7.125
7.480
4.98
7.125
0.00
2
32.000
33.556
4.86
32.000
0.00The updating parameters share the same region for the training data in the Kriging meta-modeling. To construct the Kriging meta-model valid over a range of parameters, the moduli of elasticityEb0 and Eas were restricted to vary from 50 to 70 GPa and 8 × 104 to 1 × 105 Pa, respectively. The design of experiment (DOE) is a key problem in deciding how to select the inputs, at which the deterministic computer codes are run in order to most efficiently control or reduce the statistical uncertainty of the computed predictions. In this study, the rectangular grid method [36] was used to deal with the DOE problem. This was easily done by using the full multiparameter sweep [37] and specifying all combinations type in COMSOL. A total of 25 experiments were carried out. The sampled parameter values and corresponding natural frequencies computed from FE models were used as the training data of the Kriging meta-model. A 50 × 50 uniform mesh grid in the region covered by the design sites was generated to evaluate the predictor. Figures 9 and 10 illustrate the mesh plots of the predicted values of the first and second natural frequencies at the grid points, respectively. The horizontal axes are parameters selected while the vertical axis gives the predicted response (natural frequency) at any point or location.Figure 9
Predicted values of the first natural frequency.Figure 10
Predicted values of the second natural frequency.To check the accuracy of the Kriging meta-model, the MSEs were computed for each mode as shown in Figures11 and 12. They demonstrate that all MSE values were close to zero so that the created meta-model had a high regression accuracy.Figure 11
Mean squared errors of the first natural frequency.Figure 12
Mean squared errors of the second natural frequency.Then, the FEMU was performed with the FE model replaced by the Kriging meta-model. The residuals between the accurate (predicted by Kriging meta-model) and the measured (by experiment) natural frequencies were used in the optimized objective function expressed in (5). A single-objective optimization algorithm with equal weight for each natural frequency was implemented to achieve the best minimization of natural frequency residuals. The optimization algorithm used to minimize the objective function is an improved PSO method in MATLAB. In implementing the PSO for the FE model updating, the population was taken as 50, c1 and c2 were set to be 2, and w was set to be 1. The tuning minimization process was over when the tolerances were achieved or a predefined number of iterations were reached. Table 3 shows the updated natural frequencies and their differences of beam B and illustrates good results of PSO based FE model updating. The final updated results for the parameter were Eb0 = 58.78 GPa and Eas = 9.524 × 104 Pa. As shown in Table 3, the errors between the first two natural frequencies measured and those obtained from the initial FE model were about 5% on average. When the KMM-PSO based FE model updating algorithm was used, the errors were reduced to very small values.


## 4.3. Second Step of FEMU

For the second step of FEMU of beam B under vibration shaker excitation at room temperature,ma and Ka0 were used as updating parameters of the KMM-PSO based FEMU. The numerical modal analysis, the corresponding meta-modeling, and the sequential FEMU procedure are almost the same as the first step of FEMU in Section 4.2, and in consequence the detailed processes are not repeated in this subsection. The initial values of ma and Ka0 were set as 0.2 kg and 2 × 103 N/m, and the corresponding updated values were 0.14548 kg and 2098.12 N/m, respectively. Table 4 shows the initial and updated values of the first two natural frequencies and the corresponding differences. As listed in the table, the errors between the first two natural frequencies measured and those obtained from the initial FE model were about 5.8% on average. When the KMM-PSO based FE model updating algorithm was used, the errors were reduced to tiny values.Table 4
Natural frequency differences of beam B under vibration shaker excitation before and after model updating.
Mode
Measured (Hz)
Initial (Hz)
Error (%)
Updated (Hz)
Error (%)
1
10.250
9.884
3.57
10.250
0.00
2
26.000
23.898
8.08
25.996
0.02


## 4.4. Third Step of FEMU

For the third step of FEMU of beam B under vibration shaker excitation in a UHTE, the parameters to be updated are temperature-dependent.Ebk and Kak(k=1,2) were used as updating parameters of the KMM-PSO based FEMU. The time-varying modal parameters identified in Section 3.3 were used to establish the objective function in (8). In addition, the linear functions of temperature increase for the parameters were used as initial guesses and are shown in Figures 13 and 14, respectively. Figure 6 also illustrates the first two instantaneous natural frequencies computed via the FE method with these initial parameters.Figure 13
Temperature-dependent elastic modulus of the beam material.Figure 14
Temperature-dependent stiffness of the long rod.As mentioned earlier, the temperature-dependent parameters to be identified can be expressed as lower-order polynomials of the temperature increase, but the exact order is unknown beforehand. Hence, the best order and the coefficients of these polynomials should be simultaneously identified. Without loss of generality, a simple case was considered to verify the proposed method. For this purpose,Eb and Ka were assumed to be linear or quadratic functions of the temperature increase. The process of ensuring the positivity of Young’s modulus, for instance, can be recalled as follows. By appropriately selecting the maximal and minimal values of E, ΔEk≜γEmax⁡-Emin⁡/Tmax⁡k, k=1,2, was defined, where k represents the order and γ is a weighting coefficient. In practice, γ≈1 was taken for a low order and γ≈0.1 was set for a high order to prevent strong nonlinearity of material property. Afterwards, the bounds of E with an increase of temperature were checked by using the Monte Carlo method.Different from the typical FEMU procedure, the method proposed here includes FE model selection. For instance, the beams were modelled by four competing models,Mk, k=1,…,4, as listed in Table 5. Here, each particle had 4 dimensions such that all the competing models should be searched in a 4-dimensional space, and then all of the particles were built to compute the fitness function J2 in (8) so as to find the smallest value of fitness function. Although all the models should be searched in the same space, each model was actually constrained to a particular subspace of the space. Afterwards, the most important issue for the FEMU in time domain is the implementation of the PSO algorithm. It was assumed that the values of the polynomial coefficients for the thermal-structural properties were restricted to different intervals and varied in their corresponding intervals. Table 6 presents the initial values of the parameters and the lower and upper bounds of their intervals. The factor λ2 and weighting matrix L in the objective function J2 were set as λ2=5×10-4, κ=8×10-4, and L=diag([1,3]), respectively. After substituting the objective function J1 with equal weight for each natural frequency into J2, the multiobjective optimization algorithm was implemented to achieve the best minimization of J2.Table 5
Model parameterization.
Model identity
Number of parameters
Max. order(Eb,Ka)
Parameter symbols
M
1
2
(1, 1)
E
b
1
,
K
a
1
M
2
3
(1, 2)
E
b
1
,
K
a
1
,
K
a
2
M
3
3
(2, 1)
E
b
1
,
E
b
2
,
K
a
1
M
4
4
(2, 2)
E
b
1
,
E
b
2
,
K
a
1
,
K
a
2Table 6
Particle swarm optimization parameters.
PSO parameter
Initial
Lower bound
Upper bound
Identified
E
b
1 (GPa·°C−1)
−0.054
-
0.2
0
-
4.2463
e
-
2
E
b
2 (GPa·°C−2)
0
-
9
e
-
5
9
e
-
5
-
2.0803
e
-
5
K
a
1 (N·m−1·°C−1)
−0.002
-
5
e
-
2
0
-
2.1041
e
-
3
K
a
2 (N·m−1·°C−2)
0
-
9
e
-
5
9
e
-
5
-
8.9016
e
-
5Figure13 shows the identified values of Young’s modulus Eb of the beam compared with the initial values. Figure 14 shows the identified values of the added stiffness Ka of the long rod compared with the initial values.Figure15 shows the convergence of the objective function J2 over the 50 iterations of the algorithm and illustrates that the PSO algorithm rapidly converged to the ultimate minimum error within the first 22 iterations. Figure 15 also illustrates the convergence behavior of the global best model over 50 iterations. It is quite clear from Figure 15 that the objective function played a significant role in updating the model parameters. The global best model in the simulation began with M1, then changed to M3 and M1, afterwards went to M3, and subsequently remained unchanged for the rest of the simulation. The result indicates that M3 was the global best model. That is, Eb and Ka were a quadratic polynomial and a linear polynomial of the temperature increase, respectively, to balance the model errors and complexity.Figure 15
Global fitness and global best model number during updating.To have a quantitative discussion, the mean absolute percentage error (MAPE) is defined as(26)MAPE=1Nsp∑k=1Nspyk-y^kyk×100%,
where Nsp is the total number of samplings and yk and y^k denote the true values computed by the direct modal analysis and the identified value at the kth time instant, respectively. Figure 16 illustrates the MAPEs of the first two natural frequencies. The average error between the initial values of the first two natural frequencies and the true values was 1.8%. When the KMM-PSO based FE model updating and selection method were used, the error was reduced to 0.1% on average.Figure 16
Mean absolute percentage errors of the first two natural frequencies.Overall, the KMM-PSO based FEMU approach proposed in this study updates the model of a thermal structure in both RTE and UHTE well. On the other hand, the proposed method can greatly reduce the computation time compared to other FEMU algorithms based on direct FE modeling and modal analysis by using commercial FE analysis packages. Table7 gives the comparison of computational time of the proposed method with the FEMU algorithm based on direct modal analyses in COMSOL. It should be pointed out that the FEMU algorithm based on direct modal analysis in COMSOL was not actually carried out and the corresponding computation time was predicted according to the calling numbers of Kriging predictor function and the actual time (about 5 s) of one run of FE modal analysis in COMSOL.Table 7
Comparison of computational time.
Method
Number of MA
Time of MA
Total time
The proposed method
10646400
3785 s
3925 s
FEMU based on direct MA in COMSOL
10646400
616 d
616 d
MA: modal analysis; s: seconds; d: days.


## 5. Conclusions

The paper presents the experimental study for the identification of time-varying modal parameters and the application of the Kriging meta-model for the finite element model updating and selection of a beam-like thermal structure in both steady and unsteady high temperature environments. The time-invariant natural frequencies were identified from the vibration test in room temperature and were used to establish the objective function for the FEMU in an RTE. As the thermal structure in unsteady high temperature environments exhibits the characteristics of time-varying multiphysics fields, a KMM-PSO based FE model updating and selection was proposed based on the experimentally identified time-varying modal parameters of the thermal structure in a UHTE. The presented TVAR method well extracted the instantaneous natural frequencies of the thermal structure in temperature-varying environment from the output responses of the structure only. The KMM-PSO based FEMU approach proposed in this study well updated the model of the thermal structure in both steady and unsteady high temperature environments. The integrated method was time-saving and feasible to industry. The study presents a preliminary investigation into the use of Kriging as a statistical-based approximation technique for modeling the complicated thermal structure with local joint and features of time-varying multiphysics fields.



---

*Source: 143254-2015-03-02.xml*