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.
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] of Mechanical 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 Section 2, 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: The real eigen-value problem for the th mode reads 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 field . In other words, the heat transfer process is neglected in the study. Let denote the temperature increase in the thermal structure as follows: where 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 where is the parameter vector including the temperature-dependent material and boundary parameters to be corrected, and is the general function of the temperature increase . In general, can be expressed as a low-order polynomial, for example, a linear, quadratic, or cubic function, in which represents the highest order of the polynomial. Then, the target parameters to be corrected change from temperature-dependent to a constant parameter vector , .
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: where is the error function or objective function, is the Euclidean norm, is the weighted factor for the th mode, and is the number of measured modes, respectively. In (5), is the th instantaneous natural frequency obtained from the meta-model, and is the th 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: where and represent the lower and upper bounds of the parameter coefficient vector , 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 vector in (4). It should be noted that 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 polynomial 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: with where is the number of samplings, and are frequencies at time , and is the initial estimate of parameter , respectively. The second term in the bracket represents regularization, where the parameter weighting matrix should be chosen to reflect the uncertainty in the parameter estimation, and represents the classical Tikhonov regularization. Link [29] suggested that the factor lies in the range of . 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 represents the highest order of the th polynomial. is the dimension of the parameter vector.
Similarly to the FEMU in an RTE, the optimization process here may be expressed as follows:
Figure 1 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.
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 by pbest) and that of the global best particle (denoted by gbest), which has far been found in the swarm. Let the position of a particle be denoted by , and let 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]: where is an inertia coefficient that balances the global and local search, and are random numbers in the range updated at each generation to prevent convergence on local optima, and and are the learning factors that control the influence of and during the search process. Typically, and 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 parameters 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 at an untried site as the sum of a polynomial trend model , and a systematic departure term representing low (large scale) and high frequency (small scale) variations around the trend model [35], where and 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 is a Gaussian distributed error term with zero mean and variance . The covariance matrix of is given by
In (13), each element of defined as is the spatial correlation function between any two of the sample points and .
In the current simulation, the term represents the dimension of vector and the linear regression model is chosen for the mean part of the Kriging function [36]. Accordingly can be expressed as
The Gaussian correlation function is taken as where is the component of the sample point and is the unknown correlation parameter which needs to be fitted by optimization.
To construct the Kriging model, the values of the regression coefficients must be approximated by using the generalized least squares theory at first. Then, the fitting correlation parameters can be quantified by using the maximum likelihood estimation. The generalized least-squares estimates of and , represented by and , respectively, are given in detail as where and are the vector of output and the matrix at the sample inputs, respectively, which are expressed as
With the vector , the prediction at the unsampled location can be obtained as
For the FEMU problem of thermal structures, parameters 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) as follows [36]: where is the process variance defined in (16) and 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 process (e.g., displacement, velocity, or acceleration) of order in a discrete-time as the following: where is a stationary white noise process with zero mean and variance , and are the TVAR coefficients.
Using the basis function expansion and regression approach, the TVAR process of order in a discrete-time can be expressed in matrix form as where is a stationary white noise process with zero mean and variance ; and are the weighted coefficients and the dimension of the basis functions , 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 where the forgetting factor is chosen in the interval and is typically close to one. The initial value of and can be selected as , , where and is an identity matrix.
Once the TVAR coefficients are obtained, the instantaneous natural frequencies can be derived from the conjugate roots , of the time-varying transfer function corresponding to the TVAR model as the following: where 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 Figure 2. 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.
(a)
(b)
(c)
Figure 2(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.
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. Figure 3 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.
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 Figure 2. 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.
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.
Figure 5 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 -axis, with respect to the measured temperature on the beam labeled in the right -axis.
With the comparison of Figures 5 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. Figure 7 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.
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 of the beam, the density of the connector, the added mass of the long steel rod, the density , and elastic modulus of the auxiliary surface, were assumed to be temperature-independent, while the elastic modulus of the beam material and the added stiffness 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, and yield the following polynomials of temperature increase : where and are the elastic modulus and the stiffness of the long rod at the reference temperature , respectively, and and are the coefficients independent of temperature. In the study, was taken as 2252.75 kg/m3 according to measured mass and volume of the beam, and of the connector was taken as 7949.6 kg/m3 in the same way. The value of 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) and , (2) and , and (3) and .
4.2. First Step of FEMU
For the first step of FEMU of beam B under a hammer impact at room temperature, and 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 and 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.
The 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 elasticity and 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.
To check the accuracy of the Kriging meta-model, the MSEs were computed for each mode as shown in Figures 11 and 12. They demonstrate that all MSE values were close to zero so that the created meta-model had a high regression accuracy.
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, and were set to be 2, and 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 = 58.78 GPa and = 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, and 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 and 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.
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. and 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.
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, and 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 , , , was defined, where represents the order and is a weighting coefficient. In practice, was taken for a low order and was set for a high order to prevent strong nonlinearity of material property. Afterwards, the bounds of 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, , , 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 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 and weighting matrix in the objective function were set as , , and , respectively. After substituting the objective function with equal weight for each natural frequency into , the multiobjective optimization algorithm was implemented to achieve the best minimization of .
Figure 13 shows the identified values of Young’s modulus of the beam compared with the initial values. Figure 14 shows the identified values of the added stiffness of the long rod compared with the initial values.
Figure 15 shows the convergence of the objective function 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 , then changed to and , afterwards went to , and subsequently remained unchanged for the rest of the simulation. The result indicates that was the global best model. That is, and were a quadratic polynomial and a linear polynomial of the temperature increase, respectively, to balance the model errors and complexity.
To have a quantitative discussion, the mean absolute percentage error (MAPE) is defined as where is the total number of samplings and and denote the true values computed by the direct modal analysis and the identified value at the th 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.
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. Table 7 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.
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.
Abbreviations
| AFR: | Amplitude-frequency response |
| AIC: | Akaike information criterion |
| CWT: | Continuous wavelet transform |
| FE: | Finite element |
| FEMU: | Finite element model updating |
| FRF: | Frequency response functions |
| GA: | Genetic algorithms |
| KMM: | Kriging meta-model |
| MAPE: | Mean absolute percentage error |
| MPC: | Multipoint constraint |
| MSE: | Mean squared error |
| PSO: | Particle swarm optimization |
| RLS: | Recursive least square |
| RTE: | Reference temperature environment |
| SHTE: | Steady high temperature environment |
| TVAR: | Time-varying autoregressive |
| UHTE: | Unsteady high temperature environment. |
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This work was supported by the National Natural Science Foundation of China under Grants 11472128, 11302098, and 11290151, the Funding of Jiangsu Innovation Program for Graduate Education under Grant CXLX13_130, and the Research Fund of State Key Laboratory of Mechanics and Control of Mechanical Structures (NUAA) under Grants 0114G01 and 0113Y01.