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

Abstract

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, ¹ 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 ² studied multilayer plate under heat and mechanical loads in the context of quasi-static hypothesis. Mukherjee and Sinha ³ analyzed the dynamic thermal response of thick laminated composites using finite element method. Matsunaga ⁴ presented a global high-order theory for solving the free vibration problem of the sandwich 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. ⁵ 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. ⁶ referred to the test specification of E228-1 in ASTM, using dilatometer to measure the thermal expansion coefficient of two-dimensional (2D)-woven and three-dimensional (3D)-braided composites from 35°C to 800°C. Pan et al. ⁷ 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 ⁸ 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 ¹¹ developed an inverse stochastic method based on the non-sampling generalized polynomial chaos method for identifying uncertain elastic parameters from experiment modal data. Multi-scale 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. ¹² determined the thermoelastic parameters of braided composites and established the 2.5D cell model. Qin et al. ¹³ and Fei et al. ¹⁴ 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. ²⁰ Researchers have paid much attention in recent years to the field of structural parameter identification using experimental modal data.^21–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–31 Substructuring-based model updating has better performance in efficiency for complex structures. ³² 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.

Basic theory

Free vibration of thermal composite structures

Under thermal circumstance, the free vibration equation of the composite structures without considering damping is given as

$$\mathit{M}\stackrel{··}{\mathit{x}}+\mathit{Kx}=0$$

(1)

and the generalized eigenvalue formula is

$$(\mathit{K}-\lambda \mathit{M})\mathit{\phi }=0$$

(2)

where K ∈R^n×n is the stiffness matrix, M ∈R^n×n is the mass matrix, n is the degree of freedom of the structure, λ is the eigenvalue, and φ ∈Rⁿ 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

$$\mathit{K}={\mathit{K}}_T+{\mathit{K}}_{\sigma }$$

(3)

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 _σ is the stiffness caused by the structural thermal stress, and the subscript σ indicates the stress effect. The stiffness matrix is derived as follows.

The generalized Hooke’s Law of orthotropic material under thermal circumstance is

$$\mathit{\sigma }=\mathit{D}\left(T\right)(\mathit{\epsilon }-{\mathit{\epsilon }}_0)$$

(4)

where σ is the stress vector, ε is strain vector, ε ₀ indicates the thermal induced strain and is treated as the initial strain, and D (T) is the temperature-dependent elastic matrix

$$\begin{array}{c}\mathit{D}\left(T\right)=\mathit{C}\left(T\right)=\hfill \\ \left[\begin{array}{cccccc}{C}_{11}\left(T\right)\hfill & C_{12}\left(T\right)\hfill & C_{13}\left(T\right)\hfill & \hfill & \hfill & \hfill \\ \multicolumn{1}{c}{C_{21}\left(T\right)}& C_{22}\left(T\right)\hfill & C_{23}\left(T\right)\hfill & \hfill & \hfill & \hfill \\ \multicolumn{1}{c}{C_{31}\left(T\right)}& C_{32}\left(T\right)\hfill & C_{33}\left(T\right)\hfill & \hfill & \hfill & \hfill \\ \multicolumn{1}{c}{}& \hfill & \hfill & C_{44}\left(T\right)\hfill & \hfill & \hfill \\ \multicolumn{1}{c}{}& \hfill & \hfill & \hfill & C_{55}\left(T\right)\hfill & \hfill \\ \multicolumn{1}{c}{}& \hfill & \hfill & \hfill & \hfill & C_{66}\left(T\right)\hfill \end{array}\right]\hfill \end{array}$$

(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

$$\begin{array}{l}{\mathrm{C}}_{11}\left(T\right)=\frac{{\mathrm{E}}_{11}\left(T\right)\left[1-{\mathrm{\mu }}_{23}\left(T\right){\mathrm{\mu }}_{32}\left(T\right)\right]}{\mathrm{\Delta }\left(T\right)}\\ {\mathrm{C}}_{22}\left(T\right)=\frac{{\mathrm{E}}_{22}\left(T\right)\left[1-{\mathrm{\mu }}_{31}\left(T\right){\mathrm{\mu }}_{13}\left(T\right)\right]}{\mathrm{\Delta }\left(T\right)}\\ {\mathrm{C}}_{33}\left(T\right)=\frac{{\mathrm{E}}_{33}\left(T\right)\left[1-{\mathrm{\mu }}_{12}\left(T\right){\mathrm{\mu }}_{21}\left(T\right)\right]}{\mathrm{\Delta }\left(T\right)}\\ {\mathrm{C}}_{12}\left(T\right)=\frac{{\mathrm{E}}_{11}\left(T\right)\left[{\mathrm{\mu }}_{21}\left(T\right)\mathrm{+}{\mathrm{\mu }}_{31}\left(T\right){\mathrm{\mu }}_{23}\left(T\right)\right]}{\mathrm{\Delta }\left(T\right)}\\ {\mathrm{C}}_{13}\left(T\right)=\frac{{\mathrm{E}}_{11}\left(T\right)\left[{\mathrm{\mu }}_{31}\left(T\right)\mathrm{+}{\mathrm{\mu }}_{21}\left(T\right){\mathrm{\mu }}_{32}\left(T\right)\right]}{\mathrm{\Delta }\left(T\right)}\\ {\mathrm{C}}_{23}\left(T\right)=\frac{{\mathrm{E}}_{22}\left(T\right)\left[{\mathrm{\mu }}_{32}\left(T\right)\mathrm{+}{\mathrm{\mu }}_{12}\left(T\right){\mathrm{\mu }}_{31}\left(T\right)\right]}{\mathrm{\Delta }\left(T\right)}\end{array}$$

(6)

where

$$\begin{array}{c}\mathrm{\Delta }\left(T\right)=1-{\mu }_{12}\left(T\right){\mu }_{21}\left(T\right)-{\mu }_{23}\left(T\right){\mu }_{32}\left(T\right)\hfill \\ -{\mu }_{31}\left(T\right){\mu }_{13}\left(T\right)-2{\mu }_{21}\left(T\right){\mu }_{32}\left(T\right){\mu }_{13}\left(T\right)\hfill \end{array}$$

(7)

and the shearing stiffness are

$$\begin{array}{c}{C}_{44}\left(T\right)=G_{23}\left(T\right)\hfill \\ C_{55}\left(T\right)=G_{31}\left(T\right)\hfill \\ C_{66}\left(T\right)=G_{12}\left(T\right)\hfill \end{array}$$

(8)

ε ₀ 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 ε ₀ is expressed as

$${\mathit{\epsilon }}_0=(\psi -{\psi }_0){\left[\begin{array}{cccccc}{\alpha }_1\left(T\right)\hfill & {\alpha }_2\left(T\right)\hfill & {\alpha }_3\left(T\right)\hfill & 0\hfill & 0\hfill & 0\hfill \end{array}\right]}^T$$

(9)

in which ψ is the transient or the steady-state temperature field, ψ₀ is the initial temperature field, and α₁, α₂, α₃ 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 ³⁵

$$P\left(V\right)=V_0(V_{-1}{T}^{-1}+1+V_{1}T+V_2{T}^2+V_3{T}^3)$$

(10)

where V₀, V_–1, V₁, V₂, and V₃ 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

$$\begin{array}{cc}\multicolumn{1}{c}{{\mathit{U}}^e=}& \frac{1}{2}\underset{\Omega }{\int }({\mathit{\epsilon }}^T-{\mathit{\epsilon }}_0^T)\mathit{D}\left(T\right)(\mathit{\epsilon }-{\mathit{\epsilon }}_0)\mathrm{d}\Omega \hfill \\ \multicolumn{1}{c}{=}& \frac{1}{2}{\mathit{\delta }}^{\mathrm{eT}}\underset{\Omega }{\int }{\mathit{B}}^T\mathit{D}\left(T\right)\mathit{B}\mathrm{d}\Omega {\mathit{\delta }}^e-{\mathit{\delta }}^{\mathrm{eT}}\underset{\Omega }{\int }{\mathit{B}}^T\mathit{D}\left(T\right){\mathit{\epsilon }}_0\mathrm{d}\Omega \hfill \\ \multicolumn{1}{c}{}& +\frac{1}{2}\underset{\Omega }{\int }{\mathit{\epsilon }}_0^T\mathit{D}\left(T\right){\mathit{\epsilon }}_0\mathrm{d}\Omega \hfill \end{array}$$

(11)

where Ω is the integration domain, B is the geometric matrix, and δ ^e is the nodal displacement vector.

For the free vibration problem, the potential energy in domain Ω of the external force is zero because the structure is not loaded. The total potential energy in the volume is

$${\mathbf{\Pi }}^e={\mathit{U}}^e$$

(12)

According to the principle of minimum potential energy

$$\delta {\mathit{\Pi }}^e=\frac{\partial {\mathit{\Pi }}^e}{\partial {\delta }^e}=0$$

(13)

Substituting equation (11) into equation (13)

$$\underset{\Omega }{\int }{\mathit{B}}^T\mathit{D}\left(T\right)\mathit{B}\mathrm{d}\Omega {\delta }^e-\underset{\Omega }{\int }{\mathit{B}}^T\mathit{D}\left(T\right){\mathit{\epsilon }}_0\mathrm{d}\Omega =0$$

(14)

in which the first part is the stiffness matrix K _T

$${\mathit{K}}_T=\underset{\Omega }{\int }{\mathit{B}}^T\mathit{D}\left(T\right)\mathit{B}\mathrm{d}\Omega$$

(15)

and the second part is the nodal force deduced by temperature and can be expressed as

$${\mathit{F}}_T=\underset{\Omega }{\int }{\mathit{B}}^T\mathit{D}\left(T\right){\mathit{\epsilon }}_0\mathrm{d}\Omega$$

(16)

Substituting equations (15) and (16) into equation (14), the equilibrium equation is obtained

$${\mathit{K}}_T{\delta }^e={\mathit{F}}_T$$

(17)

By solving the thermal induced deformation problem (equation (17)), we can obtain the nodal displacement δ^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

$${\mathit{K}}_{\sigma }=\underset{\Omega }{\int }{\mathit{G}}^T\mathit{SG}\mathrm{d}\Omega$$

(18)

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

$$\mathit{S}=\left[\begin{array}{ccc}{\sigma }_{11}\hfill & {\sigma }_{12}\hfill & {\sigma }_{13}\hfill \\ \multicolumn{1}{c}{{\sigma }_{21}}& {\sigma }_{22}\hfill & {\sigma }_{23}\hfill \\ \multicolumn{1}{c}{{\sigma }_{31}}& {\sigma }_{32}\hfill & {\sigma }_{33}\hfill \end{array}\right]$$

(19)

From the formula derivation, by substituting equations (15) and (18) into equation (3), the stiffness matrix K is expressed as

$$\mathit{K}=\underset{\Omega }{\int }{\mathit{B}}^T\mathit{D}\left(T\right)\mathit{B}\mathrm{d}\Omega +\underset{\Omega }{\int }{\mathit{G}}^T\mathit{SG}\mathrm{d}\Omega$$

(20)

By substituting equation (20) into equation (2) and solving the eigenvalue problem P_eig(λ, φ ), the eigenvalue λ and the counterpart eigenvector φ are obtained

$$P_{\mathrm{eig}}(\lambda ,\mathit{\phi }):=\mathrm{g}(\underset{\Omega }{\int }{\mathit{B}}^T\mathit{D}\left(T\right)\mathit{B}\mathrm{d}\Omega +\underset{\Omega }{\int }{\mathit{G}}^T\mathit{SG}\mathrm{d}\Omega ,\mathit{M})$$

(21)

It is shown that λ and φ 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, λ and φ can also be described as function of the coefficients V_i, i = −1, 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

$$\{\begin{array}{c}\mathrm{MinJ}(\mathit{P})={\mathit{\epsilon }}^T{\mathit{W}}_P\mathit{\epsilon }={\Vert {\mathit{W}}_p^{1/2}({\mathit{\theta }}^E-{\mathit{\theta }}^A\left(\mathit{P}\right))\Vert }_2^2\hfill \\ \multicolumn{1}{c}{s.t.{\mathit{P}}^L\le \mathit{P}\le {\mathit{P}}^U}\end{array}$$

(22)

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. P ∈R ^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 . ε is the residual error of the modal parameters, θ ^E and θ ^A ( P ) ∈R ^N are the experimental and numerical modal parameters with N = N_f × (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

$${\mathit{\theta }}^A(\mathit{P})=\mathrm{g}[({\mathit{K}}_T^A+{\mathit{K}}_{\mathit{\sigma }}^A),{\mathit{M}}^A]$$

(23)

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

$$\begin{array}{cc}\multicolumn{1}{c}{{\mathit{\theta }}^A(\mathit{P})}& =\mathrm{g}[{\mathit{K}}_T^A\left(\mathit{P}\left(\mathit{V}\right)\right)+{\mathit{K}}_{\mathit{\sigma }}^A\left(\mathit{P}\left(\mathit{V}\right)\right)]\hfill \\ \multicolumn{1}{c}{}& =\mathrm{h}\left[\mathit{P}\left(\mathit{V}\right)\right]\hfill \end{array}$$

(24)

where V is the coefficient vector of temperature-dependent parameters

$$\mathit{V}={[V_{-1},V_0,V_1,V_2,V_3]}^T$$

(25)

The parameter identification problem (equation (22)) is transformed as

$$\{\begin{array}{c}\mathrm{MinJ}(\mathit{V})={\mathit{\epsilon }}^T{\mathit{W}}_v\mathit{\epsilon }={\Vert {\mathit{W}}_v^{1/2}[{\mathit{\theta }}^E-\mathrm{h}\left(\mathit{P}\left(\mathit{V}\right)\right)]\Vert }_2^2\hfill \\ \multicolumn{1}{c}{s.t.{\mathit{V}}^L\le \mathit{V}\le {\mathit{V}}^U}\end{array}$$

(26)

Using the sensitivity analysis method to solve the optimization problem (equation (26)), the jth iteration step can be described as

$${\mathit{W}}_v^{1/2}[{\mathit{\theta }}^E-\mathrm{h}(\mathit{P}({\mathit{V}}^{\left(j\right)}))]={\mathit{S}}_w^{\left(j\right)}({\mathit{V}}^{\left(j+1\right)}-{\mathit{V}}^{\left(j\right)})$$

(27)

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

$${\mathit{S}}_w^{\left(j\right)}={{\mathit{W}}_v}^{1/2}\frac{\partial \mathrm{h}\left(\mathit{P}\left({\mathit{V}}^{\left(j\right)}\right)\right)}{\partial {\mathit{V}}^{\left(j\right)}}$$

(28)

in which

$${\mathit{S}}^{\left(j\right)}=\frac{\partial \mathrm{h}\left(\mathit{P}\left({\mathit{V}}^{\left(j\right)}\right)\right)}{\partial {\mathit{V}}^{\left(j\right)}}$$

(29)

For solving the linear algebra equations (equation (27)), the identified parameters at the jth iteration step are obtained as

$${\mathit{V}}^{(j+1)}={\mathit{V}}^{\left(j\right)}+\mathrm{\Delta }{\mathit{V}}^{\left(j\right)}$$

(30)

in which

$$\mathrm{\Delta }{\mathit{V}}^{\left(j\right)}={\left[{\left({\mathit{S}}_w^{\left(j\right)}\right)}^T{\mathit{S}}_w^{\left(j\right)}\right]}^{-1}{\left({\mathit{S}}_w^{\left(j\right)}\right)}^T{\mathit{W}}_v^{1/2}[{\mathit{\theta }}^E-\mathrm{h}\left(\mathit{P}\left({\mathit{V}}^{\left(j\right)}\right)\right)]$$

(31)

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

$${\mathit{S}}^{\left(j\right)}=\frac{\partial {\mathit{\theta }}^A\left({\mathit{P}}^{\left(j\right)}\right)}{\partial {\mathit{V}}^{\left(j\right)}}=\frac{\partial {\mathit{\theta }}^A\left({\mathit{P}}^{\left(j\right)}\right)}{\partial {\mathit{P}}^{\left(j\right)}}\frac{\partial {\mathit{P}}^{\left(j\right)}}{\partial {\mathit{V}}^{\left(j\right)}}$$

(32)

where θ ^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 θ , assuming the material parameters to be identified P  = [p₁p₂…p_m]^T could be defined by the direct definition as parameters. The sensitivity matrix is as follows

$$\frac{\partial {\mathit{\theta }}^A\left(\mathit{P}\right)}{\partial \mathit{P}}=\left[\begin{array}{cccc}\frac{\partial {\theta }_1}{\partial p_1}\hfill & \frac{\partial {\theta }_1}{\partial p_2}\hfill & \dots \hfill & \frac{\partial {\theta }_1}{\partial p_m}\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{\frac{\partial {\theta }_N}{\partial p_1}}& \frac{\partial {\theta }_N}{\partial p_2}\hfill & \dots \hfill & \frac{\partial {\theta }_N}{\partial p_m}\hfill \end{array}\right]$$

(33)

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

$$\mathit{V}={\left(\begin{array}{cccc}{\mathit{V}}^1\hfill & {\mathit{V}}^2\hfill & \cdots \hfill & {\mathit{V}}^{\mathrm{m}}\hfill \end{array}\right)}^T$$

(34)

in which V ^k  = (V₀ ^k , V_–1 ^k , V₁ ^k , V₂ ^k , V₃ ^k )^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

$$\frac{\partial p_k}{\partial {\mathit{V}}^k}={\left(\begin{array}{ccccc}\frac{\partial p_k}{\partial V_{-1}^k}\hfill & \frac{\partial p_k}{\partial V_0^k}\hfill & \frac{\partial p_k}{\partial V_1^k}\hfill & \frac{\partial p_k}{\partial V_2^k}\hfill & \frac{\partial p_k}{\partial V_3^k}\hfill \end{array}\right)}^T$$

(35)

The sensitivity matrix of the modal parameters with respect to the coefficients is

$$\frac{\partial \mathit{\theta }}{\partial \mathit{V}}=\left[\begin{array}{cccc}\frac{\partial \mathit{\theta }}{\partial {\mathit{V}}^1}\hfill & \frac{\partial \mathit{\theta }}{\partial {\mathit{V}}^2}\hfill & \cdots \hfill & \frac{\partial \mathit{\theta }}{\partial {\mathit{V}}^m}\hfill \end{array}\right]$$

(36)

in which

$$\frac{\partial \mathit{\theta }}{\partial {\mathit{V}}^k}={\left[\begin{array}{cccc}\frac{\partial {\theta }_1}{\partial p_k}\frac{\partial p_k}{\partial V_{-1}^k}\hfill & \frac{\partial {\theta }_1}{\partial p_k}\frac{\partial p_k}{\partial V_0^k}\hfill & \cdots \hfill & \frac{\partial {\theta }_1}{\partial p_k}\frac{\partial p_k}{\partial V_3^k}\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{\frac{\partial {\theta }_N}{\partial p_k}\frac{\partial p_k}{\partial V_{-1}^k}}& \frac{\partial {\theta }_N}{\partial p_k}\frac{\partial p_k}{\partial V_0^k}\hfill & \cdots \hfill & \frac{\partial {\theta }_N}{\partial p_k}\frac{\partial p_k}{\partial V_3^k}\hfill \end{array}\right]}_{N\times 5}$$

(37)

Therefore, the dimension of the sensitivity matrix is N row, 5×m column. If the relationship between material parameters and temperature is not complicated, the dimension can be reduced. Due to the structural responses, θ 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

$$\begin{array}{c}(\frac{\partial ({\mathit{K}}_T+{\mathit{K}}_{\sigma })}{\partial {\mathit{P}}^{\left(j\right)}}-{\lambda }^{\left(j\right)}\frac{\partial \mathit{M}}{\partial {\mathit{P}}^{\left(j\right)}}){\mathit{\phi }}_j-\frac{\partial {\lambda }^{\left(j\right)}}{\partial {\mathit{P}}^{\left(j\right)}}\mathit{M}{\mathit{\phi }}_j\hfill \\ +({\mathit{K}}_T+{\mathit{K}}_{\sigma }-{\lambda }^{\left(j\right)}\mathit{M})\frac{\partial {\mathit{\phi }}^{\left(j\right)}}{\partial {\mathit{P}}^{\left(j\right)}}=0\hfill \end{array}$$

(38)

The sensitivity matrix of eigenvalues with respect to material parameters can be obtained

$$\frac{\partial {\lambda }^{\left(j\right)}}{\partial {\mathit{P}}^{\left(j\right)}}={\mathit{\phi }}_j^T(\frac{\partial ({\mathit{K}}_T+{\mathit{K}}_{\sigma })}{\partial {\mathit{P}}^{\left(j\right)}}-{\lambda }^{\left(j\right)}\frac{\partial \mathit{M}}{\partial {\mathit{P}}^{\left(j\right)}}){\mathit{\phi }}_j$$

(39)

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 thermo-elastic 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.

3. The relative sensitivity of the intermediate variables 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).

Figure 1.

Flow chart for the thermo-elastic parameter identification.

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 × 200 mm × 4 mm, and the boundary condition and the temperature field are shown in Figure 2.

Table 1.

Material parameters of the composite.

Temperature (°C)	300	400	500	600	700	800
Elastic modulus, Ex = Ey (GPa)	83	80	78	76	75	72
Elastic modulus, Ez (GPa)	61	60	59	57	55	53
Shear modulus, Gxy (GPa)	45	44	43	41	39	37
Shear modulus, Gxz = Gyz (GPa)	37	36	35	34	32	31
Thermal expansivity, αx = αy (10−6 K−1)	3.13	3.14	3.57	3.70	4.00	4.19
Thermal expansivity, αz (10−6 K−1)	4.70	4.72	5.35	5.55	6.00	6.29
Thermal conductivity, K11 = K22 (W/(m K))	9.60	9.21	8.82	8.99	8.94	8.81
Thermal conductivity, K33 (W/(m K))	4.80	4.60	4.41	4.49	4.46	4.40
Specific heat, c (J/(kg K))	1250	1400	1450	1500	1520	1520
Density, ρ (10−9 t/mm3)	2000
Poisson’s ratio, ν	0.15

Figure 2.

Composite plate model with heat conduction direction (unit: mm).

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 non-uniform temperature field was shown in Figure 3. This temperature field calculates the thermal mode as a temperature load.

Figure 3.

The steady-state temperature field under the temperature load of 300°C–800°C.

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₁₁, C₂₂, C₃₃, C₄₄ and the thermal expansivity coefficients α_x, α_y are chosen. The intermediate function describes the relationship between parameters and temperature.

The parameter values of the thermal expansivity α_x, α_y and the stiffness coefficients C₁₁, C₂₂, C₃₃, C₄₄ 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

$$\begin{array}{c}{C}_{11}=-21.86t+94334.30\hfill \\ C_{22}=-21.62t+93072.75\hfill \\ C_{33}=-17.22t+70315.04\hfill \\ C_{44}=-16.29t+50457.00\hfill \\ {\alpha }_x={\alpha }_y=2.29\times {10}^{-9}t+2.36\times {10}^{-6}\hfill \end{array}$$

(40)

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.

Figure 4.

The fitting curve of thermo-elastic parameters: (a) C₁₁, (b) C₂₂, (c) C₃₃, (d) C₄₄, and (e) α_x, α_y.

Table 2.

Comparison of frequency between the calculated model and the fitted model.

Modal	Calculated frequency (Hz)	Fitted frequency (Hz)	Error (%)
1	855.989	854.725	−0.148
2	797.391	796.189	−0.151
3	403.540	402.495	−0.259
4	370.462	369.982	−0.130
5	203.231	205.128	0.933
6	191.151	191.831	0.356
7	515.626	514.415	−0.235
8	727.716	728.733	0.140
9	971.557	972.570	0.104
10	1116.879	1115.917	−0.086

Parameter identification

Five parameters are selected: C11-b, C22-b, C44-b, α_x-k, and α_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, α_x-k, and α_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 C22-b, α_x-k, and α_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 C11-b C22-b, α_x-k, and α_y-k, the parameters are identified with high accuracy.

Table 3.

Comparison of modal frequency between initial and experimental model.

Mode order	Experimental (Hz)	Initial (Hz)	Error (%)
1	887.260	855.989	−3.525
2	825.439	797.391	−3.398
3	401.413	403.540	0.530
4	384.749	370.462	−3.713
5	174.659	203.231	16.359
6	141.810	191.151	34.794
7	582.086	515.626	−11.418
8	781.418	727.716	−6.872
9	1068.577	971.557	−9.079
10	1185.908	1116.879	−5.821

Figure 5.

Parameters convergence (—C44-b, --- α_x, _y -k, -·-C11-b, ···C22-b).

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₁₁, C₂₂, C₄₄, α_x-k, and α_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.

Figure 6.

The fitting curves of thermo-elastic parameters: (a) C₁₁, (b) C₂₂, (c) C₄₄, and (d) α_x, α_y.

Table 4.

Comparison of experimental and updated modal frequency.

Modal	Experimental frequency (Hz)	Updated frequency (Hz)	Error (%)
1	887.260	887.259	0.000
2	825.439	825.437	0.000
3	401.413	401.408	−0.001
4	384.749	384.748	0.000
5	174.659	174.657	−0.001
6	141.810	141.806	−0.003
7	582.086	582.090	0.001
8	781.418	781.422	0.000
9	1068.577	1068.581	0.000
10	1185.908	1185.919	0.001

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.

Figure 7.

2.5D-braided equivalent model with heat conduction direction (unit: mm).

Figure 8.

Finite element model of unit cell of 2.5D-braided composites.

Figure 9.

Finite element model of 2.5D-braided composites: (a) mesoscale model and (b) equivalent model.

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.

Table 5.

Material parameters of the 2.5D-braided equivalent model.

Temperature (°C)	300	400	500	600	700	800
Ex (MPa)	186,913	185,495	184,078	182,660	181,242	179,824
Ey (MPa)	185,913	184,314	182,715	181,117	179,518	177,919
Ez (MPa)	147,807	146,518	145,228	143,938	142,649	141,359
Gxy (MPa)	45,280	44,873	44,465	44,058	43,650	43,243
Gxz (MPa)	41,586	41,204	40,821	40,439	40,057	39,674
Gyz (MPa)	49,708	49,265	48,822	48,380	47,937	47,494
αx (×10−6 K−1)	1.649	1.756	1.864	1.972	2.080	2.187
αy (×10−6 K−1)	2.424	2.585	2.745	2.906	2.906	3.067
αz (×10−6 K−1)	2.347	2.503	2.658	2.814	2.969	3.125
ν 12	0.206
ν 23	0.168
ν 31	0.415
ρ (×10−9 t/mm3)	2190

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 Figures 10 and 11. The thermal mode is obtained under the temperature field, and the parameters of composite in the range of 300°C–800°C are identified.

Figure 10.

The temperature field of refine model under 300°C–800°C.

Figure 11.

The temperature field of equivalent model under 300°C–800°C.

Determination of intermediate function

The least square method is used to fit the composite parameters, the intermediate function is constructed, and the intermediate variables are chosen as the parameters to be updated. The thermal expansivity α_x and the stiffness coefficients C₁₁, C₃₃, and C₆₆ between 300°C and 800°C are linearly fitted by the least squares method in which the fitting equation is as follows

$$\begin{array}{c}\multicolumn{1}{c}{C_{11}=-20.16t+271230.67}\\ \multicolumn{1}{c}{C_{33}=-17.81t+210792.41}\\ \multicolumn{1}{c}{C_{66}=-4.42t+51036.05}\\ \multicolumn{1}{c}{{\alpha }_x=1.08\times {10}^{-9}t+1.33\times {10}^{-6}}\end{array}$$

(41)

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.

Figure 12.

The fitting curve of thermo-elastic parameters: (a) C₁₁, (b) C₃₃, (c) C₆₆, and (d) α_x.

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 α_x-b, which are slightly more sensitive than the other parameters. The absolute value of the relative sensitivity result is shown in Figure 13.

Figure 13.

Relative sensitivity absolute value.

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 α_x-b is shown in Figure 14.

Table 6.

Comparison of frequency between finite element model and experimental model.

Mode order	Experimental frequency (Hz)	Calculated frequency (Hz)	Error (%)
1	34,139.700	37,704.040	10.440
2	42,113.280	43,091.290	2.322
3	73,021.880	70,734.720	−3.132
4	95,957.550	104,715.500	9.127
5	106,361.600	111,085.800	4.442
6	138,573.400	130,904.600	−5.534
7	138,583.800	137,250.700	−0.962

Figure 14.

Parameters convergence curve (—C66-b, --- α_x-b,-·-C11-b, C33-b).

After parameter identification, the parameter variation of C11-b, C33-b, C66-b, and α_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.

Table 7.

Comparison of modal frequency between the updated and the experimental model.

Modal	Calculate	Target value (Hz)		Initial error (%)		Updated error (%)
		Experiment	3% noise	Noiseless	3% noise	Noiseless	3% noise
1	37,704.040	34,139.700	34,871.730	10.440	8.122	0.017	−0.499
2	43,091.290	42,113.280	41,854.330	2.322	2.955	−0.537	0.156
3	70,734.720	73,021.880	72,749.920	−3.132	−2.770	0.992	1.920
4	104,715.500	95,957.550	100,245.920	9.127	4.459	0.568	1.280
5	1,110,805.80	106,361.600	110,857.620	4.442	0.206	−0.459	−0.644
6	130,904.600	138,573.400	144,464.960	−5.534	−9.387	−0.507	−3.626
7	137,250.700	138,583.800	141,375.560	−0.962	−2.918	−0.041	1.268

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.