Estimation of thermophysic properties of a bimaterial by temperature field measurement and finite element model updating

. The aim of this study is to identify simultaneously the thermal conductivity tensor and the heat capacity per unit volume of a bimaterial, whose heat conduction obeys Fourier’s law. This approach is validated by numerical simulation. The simulated temperature fields are obtained by the direct resolution of the heat conduction equation solved numerically with the help of finite element method formulation. To identify the parameters, an inverse method is developed by using the finite element model updating (FEMU) based on the Levenberg-Marquardt algorithm. This inverse finite element method approach allowed us to estimate the thermophysical parameters sought. We validated the numerical procedure by using noiseless temperature fields at different time and space steps and two types of material: an homogeneous and a bimaterial one. To be close to real conditions, the influence of the noise on the temperature fields is also studied and shows the efficiency of the inverse method. The results of this procedure show that the identified parameters are very less sensitive to the number of infra-red images varying from 40 to 80 and the number of elements ranging from 20 to 50 for a specimen size equals to 36 . 6 × 36 . 6 mm 2 .


Introduction
Heterogeneous materials, in particular composite and multi-layer materials, have opened up a wide field of research and application in many industrial sectors, and are used in many areas. The knowledge of thermophysical properties especially the conductivity tensor and heat capacity per unit volume is very important in the understanding of heat transfer and to study thermomechanical couplings of material involving such materials [1,2]. There are several methods dedicated to the identification of thermophysical properties. These methods may have been classified according to many criteria, as presented by Rodiet [3] and El Rassy [4]: the estimated parameters (conductivity, diffusivity, heat capacity), the measurement technique (with or without contact), the estimation regime (steady, transient, periodic), and the geometry of the problem (1D, 2D, 3D), etc.
Among these methods, the calorimetric method is the commonly accepted method for the quantitative * e-mail: jean.christophe.dupre@univ-poitiers.fr determination of heat capacities and enthalpies of phase transitions. The values reported in the literature of the relative uncertainty of heat capacity measurements by this method are significant but lower than 20%, with values ranging from 1300 J kg −1 K −1 to 2000 J kg −1 K −1 [5,6]. The hot wire method and the hot plate method are used to determine the thermal conductivity of isotropic materials. These methods can characterize the conductivity of materials in the range of 0.1 W m −1 K −1 to 10 W m −1 K −1 with an accuracy of 5% [7,8]. Another technique widely used in the laboratory is the flash method. This method allows the determination of thermal diffusivity. The uncertainty of measurement of thermal diffusivity, regardless of the material is not less than 3% [9,10]. These techniques of measurement known as conventional allow the knowledge of a single parameter at a time by identification. To obtain several parameters with the same test, an inverse problems have to be solved by using numerical techniques and computational algorithms. For example, the temperatures measured by the flash method can be used to estimate thermal conductivity tensor and heat capacity at the same time [11,12]. To solve the inverse problem, many works are rather oriented towards the use of a conjugate gradient method and stochastic algorithms [13,14]. The hybrid optimization strategy combining a Particle Swarm Optimization (PSO) algorithm and a gradient based method is used by El Rassy et al. [15] for resolving such complex and non-linear inverse problem. Many researchers, Kim et al. [16] and Kolesnik et al. [17,18] have investigated the inverse method for finding the linearly temperature-dependent thermophysical properties of orthotropic medium. Lesnic et al. [19] proposed a method for simultaneous estimation of the space-varying thermal conductivity tensor.
More recently, some works on thermal parameter identification have focused on characterizing the thermophysic parameters of multilayer materials by using the flash method to obtain experimental data. Thus Carr et al. [20] determined the thermal diffusivity for a heterogeneous sample comprising two layers of different materials. Tao et al. [21] developed a new method which is based on the Pulsed Thermal Imaging-Multi-layer Analysis method (PTI-MLA) that has been commonly used for measuring the thermal effusivity (the square root of the product of the thermal conductivity by the heat capacity per unit volume) of bulk materials for testing of coating materials. Ma et al. [22] have presented the exact analytical solutions of the fundamental problem of thermal conduction in anisotropic multilayer media in the permant regime. El Rassy et al. [23] have used a stochastic optimization algorithm to estimate the thermal diffusivities of multilayer materials. Similarly Li et al. [24] used the analytical integral transformation method for obtaining temperatures to simultaneously estimate the spatially variable thermal conductivity and thermal diffusivity in one-dimensional heat conduction in heterogeneous media. They have numerically solved in one-dimensional multilayer materials subjected to pulsed heating, based on a hyperbolic heat conduction equation and taking into account non-Fourier heat conduction effects. An implicit difference scheme is presented and a stability analysis, which shows that the implicit scheme for the hyperbolic equation is stable. Najafi et al. [25] presents a solution for the Inverse Heat Conduction Problem (IHCP) in a multilayer medium based on solutions from individual layers separately. The approach allows for inclusion of known contact resistances between the layers. The temperature histories are assumed known at two points on the inner layer and the heat transfer rate at the far end of the outer layer is the desired unknown parameter. Gao et al. [26] present a new inverse analysis approach the Complex-Variable-Differentiation Method (CVDM) for identifying material properties for multi-region problems using the numerical Boundary Element Method (BEM) model to obtain the temperature fields.
The objective of this work concerns the development of an identification method allowing a simultaneous estimation of the thermal conductivity tensors and heat capacity per unit volume of a orthotropic bimaterial using an experimental device constituted by an IR camera and heat electric sources. This work follows previous works [27] considering homogeneous orthotropic material. This latter was sensitive to the noisy image, and a low-pass filter was necessary. Performances of this method is about 5% in term of relative error on thermal conductivity tensor and heat capacity per unit volume for noisy Infra-Red (IR) image a noise with a standard deviation of 0.05 • C. For this study we have used FEniCs (Finite Element Computational software) a free library and software under Python language [28]. The first part of this article is devoted to the direct problem of calculating the temperature fields. In the second part the identification process is presented. The temperature fields obtained by the direct problem are used to validate and to assess the identification approach.

Direct problem 2.1 Heat conduction equation
We consider an orthotropic bimaterial and we assume that the conduction is the only mode of heat transfer. The physical problem consists in solving the linear heat conduction problem in the case of a bimaterial, with For times t > 0, two uniform heat fluxes q 1 and q 2 are respectively subjected to the surfaces x = 0 and y = 0, while the other surfaces are supposed to be isolated (see Fig. 1). We assume that the plate is a composite of small thickness and confined in a rectangular region 0 x L, 0 y l, so we can consider that the problem is 2D, where the temperature distribution in the plate is given by the function T (x, y, t).
The domain Ω = ∪Ω m of the material is initially at temperature T 0 (x, y). The mathematical formulation is given by the Fourier equation: with boundary and initial conditions (2) with − → n x and − → n y are respectively the normal at the surfaces x = 0 and y = 0. We assume that the thermal contacts between the different materials are perfect and that the contact resistances are null.

Solving the conduction equation
To solve the problem described by equations (1) and (2), we define bilinear and linear forms describing the Galerkin discretization of the variational formulation in which ψ the test function. Using a temporal discretization with a θ−scheme [29], the spatio-temporal temperature profile T is reduced to a finite set of temperatures at regularly spaced time increments, T i , i ∈ {0, 1, . . . , N t } where N t the number of time steps of the simulation and i refers to the time t i and ∆t the time steps. The time discretization process requires a specific study to assess numerically its accuracy and stability [30]. In this work we set θ = 2/3, corresponding to a Garlerkin scheme. The problem is converted to weak formulation by multiplying equation (1) with a test function ψ, following by a parts integration to obtain the operators: (3) and (4) above:

See equations
with Γ refers to all the boundary of domain Ω.
We presume T i is be known, and find T i+1 so that: for all ψ in some family of functions. From the mathematics literature, a(T i+1 , ψ) is known as a bilinear form and L(ψ) as a linear form [30]. From the equation (5), the discretization procedures of the element method are used to transform equation into a linear one. According to the finite element method and within the framework of this study [31], we obtain equation (6) on the domain , {f } respectively the global matrices of capacity, conductance and heat flux vector. For the resolution of the heat conduction problem (6), the FEniCS software [32] developed under python has been chosen. This choice was made because of the wealth libraries at our disposal for the resolution and for the mesh definition.

The principle
Equation (6) of the finite element method is used to identify the thermophysical parameters of the bimaterial. The methodology of this approach consists in solving an inverse problem for which some input data of the direct problem are deduced from the comparison between experimental values and predicted data obtained by simulations of the same problem. The definition of the cost function and the choice of the Levenberg-Marquardt method allows to solve the identification problem formulated as an optimization problem.
The inverse problem is an optimization problem to find the vector parameters P = ρc 1 , k 1 xx , k 1 yy , ρc 2 , k 2 xx , k 2 yy that minimize the difference between the experimental temperatures Y ij and the predicted temperatures (T ij (P )). The cost function is given by : with index i refers to the time t i while the index j refers to the measurement points (nodes number), where i = 1, . . . , N t and j = 1, . . . , N mes , with N t is the number of time steps of the simulation and N mes the number of measurement points (nodes) [33]. Y ij is the experimental temperatures at the time t i for the node j.

Levenberg-Marquardt method
There are a variety of methods to carry out a function minimization. In this study we explored Levenberg-Marquardt iterative method. To minimize the least squares norm giving by equation (7), we have to reduce to zero the derivatives of f (P ) with respect to each unknown parameters P k , that is: where N is the number of parameters to be estimated.
To solve the resulting system of algebraic equation (8), the Levenberg-Marquardt iterative method is chosen. This algorithm combines the steepest descent and Newton methods. Starting the iterations with a large value of the Levenberg-Marquardt parameter µ, more emphasis is given initially to the steepest descent method, since a good initial guess is not required with this method; but the convergence is slow. As the value of the parameter µ is gradually reduced at each iteration step, the weight is increasingly shifted to the Newton method which converges faster [34,35]. The sensitivity coefficients, can be calculated from the following finite difference formula: The parameters are obtained by the following equation: I is the identity matrix, the exponent n is the iteration index, J the sensitivity matrix and J * is the transpose of the matrix J. The stopping criterion for the algorithm can operate if the change in P is small: where ε is a value for the stopping criterion, k refers to the number of parameters [35].

Results and analysis
For the numerical simulation, we have chosen a squareshaped specimen with a length L = 36.6 mm and a width l = 36.6 mm. The initial temperature is T 0 (x, y) = 20 • C and the heat flux to which the plate is subjected is chosen constant q 1 = q 2 = 396.99 W/m 2 . We used two types of materials: a polyamide (PA) material and a carbon fiber reinforced epoxy matrix composite (CFRP) whose properties are summarized in Table 1.
In this study, three parameters has been assessed: M the number of elements along X or Y directions (the horizontal and vertical resolution of the IR images), dt = the time step (with t f the duration of the test and N t the number of recorded IR images) and the temperature noise level in the IR images.
To be close to the experimental conditions, a random Gaussian noise is added to the simulated temperature field as follows where T ij is the temperature obtained by the direct approach at time t i and at node j, where w j is a random number generated from a normal distribution with mean zero and constant standard deviation σ [36]. In order to obtain statistical representative identified parameters, a large number of identifications are carried out (n I = 20) with different random series of noise which allows to evaluate for every identified parameter P k , the average value defined: whereP l k is the estimation of parameter P k for the random series n I . The relative errors is computed by comparing the estimated values α est with imposed values α imp .

Identification of the parameters of a structure made of one material
The validation of this approach was first done by considering that the part one (m = 1) and the part two (m = 2) are similar. The chosen material is the polyamide. This procedure allows us to compare this new approach with the previous ones [27] developed for homogeneous material.

Validation of the numerical procedure
To validated the numerical procedure, we simulated the temperature fields without noise as input data to the inverse algorithm to estimate the thermal conductivity tensor and the heat capacity per unit volume. We have a good agreement between the estimated parameters and the simulation parameters. Regarding Figure 2, the accuracy in the results is highly satisfactory since the overall of the identified parameters show an error under 0.1%, except for the k 2 yy (Fig. 2d), for with the relative error has an order of magnitude of 2.5%. The relative errors of the heat capacity per unit volume are around 0.05% to 0.15% for a time step dt = 1s. For time steps greater than dt = 2s, the relative errors are very small and are around 0.05% (Fig. 3). These results show how well the algorithm works.

Influence of noisy infra-red images on the identified parameters
To show the behaviour of the procedure with noisy images (more realistic), we have added a noise (12) with zero mean and constant standard deviation σ = 0.4 • C. Figures 4 and 5 show that the measurement noise has a direct effect on the results. The relative errors have the same tendencies as those of the noiseless temperature fields. For the number of elements higher than 20 and the time step higher than dt = 2 s, we have almost the same relative errors values. The relative errors obtained on the estimated parameters of the thermal conductivity tensors are under 1% for k xx , 4% for the case of the k yy coefficient (Fig. 4) and 3% for the heat capacity per unit volume ρc 1 (Fig. 5a), and 0.8% for ρc 2 (Fig. 5b). Thus a number of elements along each axis (M = 30) and a time step dt = 4 s seem to be a good compromise between the number of elements and the computation time for the identification.

Influence of the noise amplitude
According to the analysis done of the previous results of Figures 4 and 5, the specimen is discretized with M = 30 elements in X and Y directions and a time step dt = 4s is chosen. Four levels of standard deviation noise σ = {0.1 • C, 0.4 • C, 1.2 • C, 2 • C } are used. However, the higher the standard deviation of the noise, the relative error increases (see (Fig. 6). A relative error of 0.8% is  reached for a noise of 0.1 • C and 2.2% for a noise of 0.4 • C. Previous study [27] using a similar approach for an homogeneous material, gives an uncertainty of about 5% for the two components of the conductivity tensor and 4% for the heat capacity per unit volume for comparable conditions (M = 31 and dt = 4 s). Nevertheless a spatial median filter (3 × 3 pixels) coupled with an average temporal filter (on 3 images) have to be used to obtain simular results, it is not the case in this procedure. The Figure 7 shows the gap between a reference temperature field (T ref ) without noise and the temperature field (T id ) obtained from the identified parameters. This difference represents the residual in temperature space of the identified parameters from noisy images of standard deviation σ = 0.4 • C at the time t = 80s (Fig. 7c). This difference is of the order of 2 × 10 −2 • C and increases at the vertical border about 0.08 • C. The difference is not the same in the the two parts because the identified parameters, theoretically similar, are not the same in the two parts (several % of differences).

Identification of the parameters of bimaterials
The procedure validated in the previous section is now used for a bimaterial. The specimen used consists of two materials with different properties. The first one is a carbon fiber reinforced epoxy matrix composite (CFRP) and the second polyamide (PA) (see Tab. 1). The juxtaposition of theses two materials CFRP-PA constitutes the bimaterial.

Influence of space steps and time steps on the identified parameters
The curves in Figures 8 and 9 are obtained by identifying the parameters for a noise level of standard deviation σ = 0.4 • C on the temperature fields. We observe the same tendencies as the previous case of one material with a decrease of relative errors versus number of elements and time step. The relative errors are 2.5% for the coefficients of the thermal conductivity tensor k 1 xx (Fig. 8a) and 22% for k 2 yy (Fig. 8d) for a small number of elements (M = 20). The relative errors of the heat capacity per unit volume are of the order of 1.2% (Fig. 9b) for a small time step dt = 1s. Overall results for these estimations are satisfactory. As noted before, the estimation of the heat capacity per unit is precise when compared to that of the thermal conductivity tensor.

Influence of the noise amplitude on the identified parameters
In the case of these simulations, the number of elements M = 30 and the time step dt = 4 s, have chosen. Four levels of standard deviation σ = {0.1 • C, 0.4 • C, 1.2 • C, 2 • C} are used. However, the higher the standard deviation of the measurement noise, the relative error increases (see (Fig. 10). A relative error of 0.5% is reached for a noise of 0.1 • C and 2% for a noise of 0.4 • C. Figure 11 shows the temperature fields with and without noise and the difference of the residual in temperature space of the identified parameters from noisy images of standard deviation σ = 0.4 • C at the time t = 80 s (Fig. 11c). This difference is of the order of 1 × 10 −2 • C and increases at the vertical border about 0.006 • C.

Conclusion
The developed procedure allowed us to have the good estimates. For this purpose, the application of the finite element method allowed to reduce the partial differential equations of heat conduction into a system of ordinary differential equations. The numerical integration of this system of equations led to the simulation of temperature fields T (x, y, t) in two dimensions. This algorithm  offers a simple, fast and systematic identification of thermal parameters of bimaterials from measured temperature fields. For a number of elements M = 30 (∆x = ∆y = 1.18 mm) and the time step dt = 4 s which corresponds a number N t = 60 of temperature measurements, the relative error is around 2% for the thermal conductivity tensor. This relative error is 0.3% for the heat capacity per unit volume for standard noisy image (σ = 0.4 • C). These relative error decrease to 0.03% for σ = 0.1 • C. One of the main results is the procedure that is almost insensitive to temperature noise. This method is less sensitive to noise than the previous works [27]. The continuation of this work, in simulation, consists in studying other types of more or less conductive materials in order to see the influence on the accuracy of the parameters. Experimental test will be performed to assess the procedure with real temperature fields.