Temperature-dependent thermo-elastic parameter identification for composites using thermal modal data

To accurately determine the temperature-dependent parameters of composites, a thermo-elastic parameter identification approach using thermal modal data is proposed in this article. The investigation is based on two hypotheses: (1) the structure is at steady-state temperature field, which means the temperature distribution is time independent; (2) temperature distribution can be determined in advance of thermo-elastic parameter identification, and thermal-related large deformation is not considered. The thermal-dependent elastic constants and thermal expansion coefficients are expressed as intermediate functions with the independent variable of temperature, perturbation method is adapted to calculate the sensitivity of thermal modal frequencies with respect to the intermediate variables, and variables with high sensitivity are selected as the identifying parameters. By constructing the intermediate function and calculating the relative sensitivity, thermo-elastic parameters are identified by minimizing the residual between experimental and analytical natural frequencies under thermal circumstance. Two different types of composite models are employed to verify the strategies of parameter identification. Results show that on the basis of the intermediate function, the proposed approach can be applied to identify the thermo-elastic parameters effectively.


Introduction
High temperature and temperature gradient have great influence on the dynamic performance of composite structures, 1 design of which is based on the mechanical property of composites. Elastic parameters of composites under thermal circumstance vary with temperature, and determination of these material parameters has drawn much interests.
Thermal structural dynamics has been widely investigated by theoretical analysis, numerical simulation, and experimental measurement. Savoia and Reddy 2 studied multilayer plate under heat and mechanical loads in the context of quasi-static hypothesis. Mukherjee and Sinha 3 analyzed the dynamic thermal response of thick laminated composites using finite element method. Matsunaga 4 presented a global high-order theory for solving the free vibration problem of the sandwich 1 plates subjected to thermal loading. Meanwhile, heat will cause thermal expansion of the structure. The expansion will produce thermal stress inside the structure and affect the mechanical performance of the structure. 5 Therefore, it is of great theoretical and technical significance for the study on the thermo-elastic behavior of composites. Experimental measurement and theoretical analysis are the major approach for predicting thermo-elastic parameters. Gao et al. 6 referred to the test specification of E228-1 in ASTM, using dilatometer to measure the thermal expansion coefficient of twodimensional (2D)-woven and three-dimensional (3D)braided composites from 35°C to 800°C. Pan et al. 7 tested the longitudinal compression performance of 3D-braided composite materials using the separated Hopkins pressure rod combined with the heating device, material mechanical properties, and thermal physical properties were obtained in the range of 23°C-210°C. Schapery 8 deduced the bounds of the effective thermal expansion coefficient of isotropic and anisotropic composites using the principle of thermal elasticity.
The identification of composite material parameters belongs to the category of inverse problem. 9, 10 Sepahvand and Marburg 11 developed an inverse stochastic method based on the non-sampling generalized polynomial chaos method for identifying uncertain elastic parameters from experiment modal data. Multiscale method is an important method to predict the cyclical composite parameters. However, the modeling of cells is complex, and the composites in the preparation will also lead to different shapes. Chen et al. 12 determined the thermoelastic parameters of braided composites and established the 2.5D cell model. Qin et al. 13 and Fei et al. 14 identified the thermal-related parameter of braided composites at high temperature. A large number of parametric modeling of different periodic composites are also carried out in previous studies. 15 -19 Model updating-based parameter identification has been widely applied in engineering. 20 Researchers have paid much attention in recent years to the field of structural parameter identification using experimental modal data. [21][22][23][24][25][26] Experimental frequency response function combined with modal data were used to identify the damping matrix. 27,28 The response surface model as a surrogate model for identifying the elastic parameter can be found in previous studies. [29][30][31] Substructuringbased model updating has better performance in efficiency for complex structures. 32 However, there are few studies on the identification of mechanical behaviors of composites at high temperature.
Assuming that the structural dynamic performance is uncoupled with temperature, high temperature environment affects the structural vibration response non-reversing. According to the above hypothesis, the temperature distribution and structural dynamic characteristics of the structure at high temperature are separated. The accuracy of the structural dynamics model in high-temperature environment is composed of structural temperature distribution and structural dynamic characteristics.
In this article, thermo-elastic parameters are identified using finite element model updating. 33,34 The elastic constants and thermal expansion coefficients are expressed as intermediate functions with the independent variable of temperature, perturbation method is adapted to calculate the sensitivity of thermal modal frequencies with respect to the intermediate variables, and variables with high sensitivity are selected as the identifying parameters. Minimizing the discrepancy between calculation and experimental thermal modal data is the objective function. The process of identifying thermo-elastic parameters is transformed into an iterative solution to the optimization problem.

Free vibration of thermal composite structures
Under thermal circumstance, the free vibration equation of the composite structures without considering damping is given as and the generalized eigenvalue formula is where K2R n 3 n is the stiffness matrix, M2R n 3 n is the mass matrix, n is the degree of freedom of the structure, l is the eigenvalue, and u2R n is the eigenvector. It is assumed that the density of the composite material does not change with the temperature, so that the mass matrix M is temperature independent. The temperature and temperature gradient have great influence on the stiffness matrix K, which will lead to the change in the structural dynamics characteristics. Considering the temperature effects under small deformation assumption, stiffness matrix K is composed of two aspects in which the first part K T represents the thermal-related variation of material properties, including the elasticity and thermal expansivity, the subscript T indicates the temperature effect. And the second part K s is the stiffness caused by the structural thermal stress, and the subscript s indicates the stress effect. The stiffness matrix is derived as follows.
The generalized Hooke's Law of orthotropic material under thermal circumstance is where s is the stress vector, e is strain vector, e 0 indicates the thermal induced strain and is treated as the initial strain, and D(T) is the temperature-dependent elastic matrix ð5Þ C(T) represents the macrostiffness matrix of orthotropic composite material and is the function of temperature, and the nonzero elements can be expressed as and the shearing stiffness are e 0 is the thermal strain, which can be treated as the initial strain of structures. Thermal deformation can only induce normal strain, and the shear strain is zero, so that e 0 is expressed as in which c is the transient or the steady-state temperature field, c 0 is the initial temperature field, and a 1 , a 2 , a 3 is, respectively, the coefficient of thermal expansion in three different directions. Composite material properties such as the elastic parameters and the thermal expansivities are always varying with temperature, using letter P to represent the temperature-dependent parameter and which can be expressed as 35 where V 0 , V -1 , V 1 , V 2 , and V 3 are the coefficients of temperature, which is treated as an intermediate variable for parameter identification.
Using the principle of minimum potential energy to deduce the stiffness matrix K, the strain elastic energy U e in the unit volume is where V is the integration domain, B is the geometric matrix, and d e is the nodal displacement vector. For the free vibration problem, the potential energy in domain V of the external force is zero because the structure is not loaded. The total potential energy in the volume is According to the principle of minimum potential energy Substituting equation (11) in which the first part is the stiffness matrix K T and the second part is the nodal force deduced by temperature and can be expressed as Substituting equations (15) and (16) into equation (14), the equilibrium equation is obtained By solving the thermal induced deformation problem (equation (17)), we can obtain the nodal displacement d e , and the deformation will cause additional stiffness effects on the structure, which is the so-called thermal stress stiffness matrix and can be expressed as where G is the differential matrix of shape function with respect to space coordinate, which is related to the node coordinates, and S is the stress matrix From the formula derivation, by substituting equations (15) and (18) into equation (3), the stiffness matrix K is expressed as By substituting equation (20) into equation (2) and solving the eigenvalue problem P eig (l, u), the eigenvalue l and the counterpart eigenvector u are obtained It is shown that l and u are functions of the temperature-dependent elastic matrix D(T). According to the relationships between material parameters in D(T) and temperature in equation (10), and only considering the effects of material parameters on the dynamics characteristics, l and u can also be described as function of the coefficients V i , i = 21, 0, 1, 2, 3.

Identification of thermo-elastic parameter
Parameter identification. Using structural response to identify the material parameters is an inverse procedure, and an optimization problem is formulated with the objective function of minimizing the discrepancies between the numerical and the experimental modal data and is defined as where J(P) is the objective function that represents the weighted residual of the test and the calculated eigenvalue and W p is a diagonal weighted matrix which reflects the relative weight of each eigenvalue residual in the objective function. P2R m is the selected parameter to be identified, m is the number of selected material parameters, and then the parameter identification considering the effect of high-temperature environment can be attributed to seek optimal solution subjected to the parameter range P L ł P ł P U . e is the residual error of the modal parameters, u E and u A (P) 2R N are the experimental and numerical modal parameters with N = N f 3 (N m + 1), N f and N m are the number of selected modal frequencies and measured coordinates for each mode, the superscript E and A, respectively, indicates experimental and analytical data, which according to equation (21) can be expressed as a function of structural stiffness matrix and mass matrix When only taking temperature-dependent composite material parameters into account and neglecting other effects, assuming the mass distribution is accurate, modal parameters can be presented as function of the coefficients V where V is the coefficient vector of temperaturedependent parameters The parameter identification problem (equation (22)) is transformed as Using the sensitivity analysis method to solve the optimization problem (equation (26)), the jth iteration step can be described as in which the superscript (j) indicates the jth step in the iteration. For parameter identification, S w (j) represents a weighted sensitivity matrix of modal parameters with respect to the coefficients of material parameters and can be represented as in which For solving the linear algebra equations (equation (27)), the identified parameters at the jth iteration step are obtained as in which After the convergence of V, and substituting V into P(V), the temperature-dependent material parameters are obtained.
Sensitivity analysis. Calculating the sensitivity matrix S (j) is a crucial step in parameter identification. Using the chain rule in derivation and letting P (j) represents P(V (j) ), S (j) is written as where u A (P (j) ) is the numerical modal data of the jth step in iteration. Modal frequencies and mode shapes are two parts contained in the selected structural responses u, assuming the material parameters to be identified P = [p 1 p 2 .p m ] T could be defined by the direct definition as parameters. The sensitivity matrix is as follows The coefficient vector of temperature-dependent parameter is represented by equation (25), and according to equation (10), and different material parameters are always with different coefficients, the identifying coefficients vector is ) T is the coefficient of the kth material parameter to be identified, k = 1, 2, 3, ., m. The derivation of the kth element of temperature-dependent material parameter P with respect to the counterpart coefficients V k is The sensitivity matrix of the modal parameters with respect to the coefficients is in which Therefore, the dimension of the sensitivity matrix is N row, 5 3 m column. If the relationship between material parameters and temperature is not complicated, the dimension can be reduced. Due to the structural responses, u contains modal frequencies and mode shapes, the sensitivity of which respect to material parameters can be obtained through derivation of the free vibration equation.
In the process of finite element analysis, the local stiffness matrix and local mass matrix of a single element are generated by the physical parameters such as geometry, material and properties of each element, and then assembled into the global stiffness matrix and mass matrix of the structure. According to equation (2), the differential equation for the parameter P (j) is The sensitivity matrix of eigenvalues with respect to material parameters can be obtained Similarly, sensitivity matrix of mode shapes with respect to material parameters in the structure can also be obtained. The advantage of intermediate variable method is that it can effectively characterize the relationship between composite material parameters and temperature, and the temperature-dependent thermoelastic parameter identification is transformed to identify the coefficients of the intermediate variables.
Implementation procedure. This article presents a parameter identification method of structural dynamics model at high temperature. Because the material parameters (coefficient of thermal expansion and other parameters) are temperature-dependent, the relationship between the parameter and the temperature is expressed as the intermediate function at different temperature conditions. Then, the coefficients of the intermediate function reduce the number of parameters to be updated. Main steps are as follows: 1. According to the geometric parameters and the initial material parameters, the initial finite element model is established to calculate the thermal mode parameters at high temperature. 2. The relationship between the material parameters and the temperature is assumed to be an intermediate function, the coefficient of which is determined as the intermediate variable.

The relative sensitivity of the intermediate vari-
ables is calculated by the difference method, and the parameters to be updated are selected according to the sensitivity analysis results. 4. The residuals of mode frequency between simulation and experiment are constructed as objective functions, and the model are optimized by the iterative optimization until the parameters converge to obtain the structural of dynamic model at the high temperature ( Figure 1).

Case studies
Parameter identification of a homogenized composite panel The material parameters of the C/SiC composite panel is shown in Table 1, and the elastic modulus, shear modulus, and thermal expansion coefficient of composites are varying with temperature. And the geometrical size of the panel is 400 mm 3 200 mm 3 4 mm, and the boundary condition and the temperature field are shown in Figure 2.
Assuming that the temperature is evenly distributed along the length of the composite panel, regardless of the effect of air convection, the ambient temperature and the initial temperature of the composite plate are T amb = T init = 20°C. The heat transfer finite element model is established using solid element in NASTRAN, which is a commercial finite element software, and the steady-state temperature field of composite panel is calculated by SOL106 solver. The temperature field was applied at 800°C on the upper side of the panel and the temperature on the lower side was 300°C. The nonuniform temperature field was shown in Figure 3. This temperature field calculates the thermal mode as a temperature load.

Determination of intermediate function.
Since the structural stiffness coefficient is affected by the temperature in the thermal environment that the structural thermal mode are changed. The stiffness coefficients C 11 , C 22 , C 33 , C 44 and the thermal expansivity coefficients a x , a y are chosen. The intermediate function describes the relationship between parameters and temperature.
The parameter values of the thermal expansivity a x , a y and the stiffness coefficients C 11 , C 22 , C 33 , C 44 between 300°C and 800°C are linearly fitted by the least square method. The abscissa of the fitting curve is the temperature value, the coordinate is the parameter value, and the fitting equation is as follows which indicates that the material parameters approximate linear distribution in the range of 300°C-800°C, t is the parameter of temperature, the equations can be written uniformly as P = kt + b, k and b is respectively the variable of the first order and the zero order of temperature. Figure 4 shows the six parameters of the linear fitting curves. The material parameters and stiffness coefficients after fitting are substituted into the recalculated thermal mode analysis. Comparison of frequency between the calculated model and the fitted model is shown in Table 2, it is found that the modes are not shifted, and the modal frequency values of the first 10 orders are with high accuracy.
Parameter identification. Five parameters are selected: C11-b, C22-b, C44-b, a x -k, and a y -k, and the relative sensitivity values of which are slightly larger than the other parameters. Simulation study is conducted to verify the parameter identification procedure. Assuming the real value of C44-b is 1.2 times the initial value, and C11-b C22-b, a x -k, and a y -k are 1.08 times each of the initial value, the experimental thermal modal data are obtained by assigning the perturbed values to the corresponding parameters. The first 10 order of modes exactly match before and after the perturbation. The comparison of modal frequency between initial and experimental model is shown in Table 3, the fifth mode and the sixth mode are switched. Figure 5 shows the variation of parameters in the convergence procedure.   As can been seen from the figure, after parameter convergence, the parameter variation of C44-b, C11-b C22b, a x -k, and a y -k is 20%, 8%, 8%, and 7.8% of the initial value, respectively; comparing with this result, the perturbed value of the initial parameter is 120% of the initial value of C44-b, 108% of the initial value of C11b C22-b, a x -k, and a y -k, the parameters are identified with high accuracy.
In order to further testify the identification procedure, the initial parameters and its fitting curve, the perturbed parameters, the fitting equation of modified parameter are plotted together to compared. The initial fitting curves of the five material parameters C 11 , C 22 , C 44 , a x -k, and a y -k and the perturbation and the   updated curves are shown in Figure 6. Comparison of experimental and updated modal frequency is shown in Table 4. It can be seen from the figures that the updating parameter curves are almost coincident to the curves of experimental parameters, indicating the effectiveness of parameter identification method.

Parameter identification of a 2.5D-braided composite panel
Investigation on a 2.5D-braided composite panel is conducted, and the geometrical parameters, the temperature field, and boundary conditions are shown in Figure 7. Simulated experimental thermal modal data can be obtained using the mesoscale model, which is constructed according to the geometric model of 2.5D-braided composites and the following four assumptions: (1) the warp weft axis is a straight line, (2) weft yarns have the same cross sections, (3) no gap exists between the substrate and the yarn, and (4) the knitting angle of the yarns is the same. The volume fraction of fiber is 40%, and the weaving inclination of the warp yarns is 15°. The unit cell of the composite is shown in Figure 8, which is constructed using solid element. Based on the model of unit cell of the braided composite, the adopted model for simulation study is six times of the unit cell model which is shown in Figure 9. Equivalent modeling is widely employed for composite structural analysis, and the thermo-elastic parameters and the expansion coefficient of the equivalent model are shown in Table 5.
Similar to the previous example, the heat transfer finite element model is established in NASTRAN. The steady-state temperature field of the composite panel is calculated. The temperature is 300°C on the lower side of the panel and the temperature on the upper side is 800°C. The non-uniform temperature field is shown in   Figure 12 shows the four parameters of the linear fitting curve. The initial parameter points are almost distributed in the fitting line, because the material     parameter field of the yarn changes linearly with the temperature. The material parameters and stiffness coefficients after fitting are substituted into the recalculated thermal mode analysis. The thermal modal frequency of equivalent model is obtained from the initial parameters, and the equivalent results still have considerable errors.
Updating parameters selection. The least squares fitting equation is assumed to be a linear relationship P = kt + b. Assuming that the slope coefficient of each parameter fitting equation is k, the intercept coefficient is b. The 24 parameters obtained by fitting the thermal expansivity coefficients and stiffness matrix coefficients are calculated by the difference method. Selecting the four parameters: C11-b, C33-b, C66-b, and the thermal expansivity coefficient a x -b, which are slightly more sensitive than the other parameters. The absolute value of the relative sensitivity result is shown in Figure 13.  Parameter identification. The thermal modal frequencies calculated from the refined model are treated as the experimental modal frequency. Then, the calculation and experiment of the first seven-order modal frequency as well as error is shown in Table 6, and the maximum error is 10.44%. The iterative convergence of parameters C11-b, C33-b, C66-b, and a x -b is shown in Figure 14.
After parameter identification, the parameter variation of C11-b, C33-b, C66-b, and a x -b is, respectively, -13.66%, 9.62%, 32.22%, and 44.46% of the initial value. Table 7 shows comparison of modal frequency between the updated and the experimental model, the maximum error of the modal frequency between the initial and the experimental model is 10.44%, and the error between the updated frequency and experimental frequency is less than 0.568%. For testing the robustness of the parameter identification method, 3% noise is added to the test data, and the error of the identified model is less than 1.34%, indicating that the parameter identification based on thermal modal frequency is effective and robust.
According to the case studies of the parameter identification of composites, the effectiveness of the proposed method has been testified. After parameter identification, the error of the modal data of the composite model is greatly reduced and the thermo-elastic parameter can be accurately determined, and the robustness of the proposed method is illustrated by adopting the noised experimental data.

Conclusion
An approach for composite thermo-elastic parameters identification is proposed in this article. Constructing the intermediate function for thermo-elastic parameters with respect to the temperature, the updating of undetermined coefficient will improve the calculation efficiencies. Parameter sensitivity analysis is conducted using the finite difference method. The intermediate function coefficients can be identified by minimizing the residual between the analytical and experimental thermal frequencies. The major advantage is that the proposed method will determine the thermo-elastic parameters by considering the effect of non-uniform temperature and provide a convenient method to accurately predict elastic parameters and thermal expansion coefficients of composites. The composite parameters  always vary nonlinearly depending on the temperature, which will introduce difficulties into parameter identification and should be investigated in the future.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.