© AFM, EDP Sciences 2019

1 Introduction

When a material is subjected to mechanical solicitations, its temperature undergoes a few variation. The deformation processes induce thermal and energy effects which must be correctly predicted to properly model the mechanical behavior. For quasi-static tests, the temperature changes induced by deformation are often weak. The neglect of these variations eliminates any possibility of a complete energy balance because these small temperature variations induce large changes in energy. For more information on the mechanical response, it is necessary to perform an energy balance based on the thermodynamics of irreversible processes. This energy balance requires the knowledge of the thermal conductivity tensor k and heat capacity per unit volume c involved in the heat conduction [1–3].

Many numerical and experimental methods have been developed to determine the thermal conductivity tensor and sometimes the thermal conductivity tensor and the heat capacity per unit volume simultaneously. Inverse estimation of thermal parameters from temperature measurement has been studied in recent years [4–6]. In this paper, we propose a method using the same device than the mechanical study i.e. infrared camera and simple set-up. Several materials have direction-dependent thermal conductivities like for example composite materials [7,8]. Such kind of materials is denoted anisotropic, as in opposition to isotropic materials in which the thermal conductivity does not vary with the direction [9]. We are looking for the anisotropic homogeneous materials.

The evaluation of the thermal parameters as thermal conductivity tensor k and heat capacitance ρc is necessary in the characterization of new materials. They appear in the characterization of new polymers and in the analysis of thermal or heat insulation of materials [10–12].

For the estimation of the thermal properties of orthotropic materials, inverse analysis techniques based on Levenberg–Marquardt method or conjugate gradient method have been devoted [13–15]. Andersson and Ross [16] have used the transient hot-wire method to measure simultaneously the thermal conductivity and heat capacity per unit volume; de Carvalho and Neto [17] combined the hot-wire method and Levenberg–Marquardt method for polymer thermal properties estimation.

In this work, an identification process is developed to determine heat capacity and thermal conductivity tensor. The Fourier’s partial differential equation governing heat conduction as a differential equation is established by using the finite element method (FEM). That is well adapted to solve the problem by numerical computation. Although the FEM is well known, some terms will be developed to be included within the inverse approach. The identification of the unknown parameters k and ρc is performed from the minimization of a least square functional.

This study is a continuation of our work developed in the one-dimensional (1D) case [18] where we have demonstrated the efficiency of this approach on experimental tests. The two-dimensional (2D) extension is presented in this paper using a similar experimental device. The accuracy and the robustness of the inverse approach are assessed by the analysis of simulated temperature fields obtained from the 2D direct heat conduction problem. For that, different cases have been studied with several noise levels added to the temperature fields.

2 Direct problem

2.1 Modeling the thermal heat conduction in 2D

The physicalproblem consists here to solve the linear heat conduction in the case of an orthotropic material, with k a 2D tensor of the thermal conductivity and ρc a scalar value of the heat capacity per unit volume. The domain Ω is a solid initially at zero temperature. For times t > 0, two uniform heat fluxes q₁ and q₂ are respectively supplied at the surfaces x = 0 and y = 0 while the others are supposed insulated Figure 1. We assume that the plate is a composite with a small thickness and the length is very close or equal to the width, we can legitimately consider that the problem is 2D, where the temperature distribution within the plate is given by the function T(x, y, t).

The mathematical formulation of such physical problem is given by the Fourier’s equation $$\begin{array}{lll}\hfill & \hfill & \text{div}(k.\overrightarrow{\text{grad}}(T(x,y,t)))-\rho c\frac{\partial T(x,y,t)}{\partial t}\hfill \\ \hfill & \hfill & =0\text{\hspace{1em}in}\Omega \times [0,T_f];t>0.\hfill \end{array}$$(1)

with the boundary and initial conditions: $$\begin{array}{lll}(k.\overrightarrow{\text{grad}}(T(x,y,t))).{\overrightarrow{n}}_x\hfill & \text{}=\text{}\hfill & q_1\text{\hspace{0.17em}at\hspace{0.17em}}x=0\text{ \hspace{0.17em}for\hspace{0.17em}}t>0\hfill \\ \hfill & \hfill & \hfill \\ (k.\overrightarrow{\text{grad}}(T(x,y,t))).{\overrightarrow{n}}_y\hfill & \text{}=\text{}\hfill & q_2\text{\hspace{0.17em}at\hspace{0.17em}}y=0\text{\hspace{0.17em}for\hspace{0.17em}}t>0\hfill \\ \hfill & \hfill & \hfill \\ (k.\overrightarrow{\text{grad}}(T(x,y,t))).\overrightarrow{n}\hfill & \text{}=\text{}\hfill & 0\text{\hspace{0.17em}elsewhere\hspace{0.17em}for\hspace{0.17em}}t>0\hfill \\ \hfill & \hfill & \hfill \\ T(x,y,t=0)\hfill & \text{}=\text{}\hfill & T_0(x,y)\hfill \end{array}$$

where n_x and n_y are respectively the normal at the surfaces x = 0 and y = 0 and k: $$k=\left(\text{}\begin{array}{cc}{k}_{11}& k_{12}\\ k_{22}& k_{22}\end{array}\text{}\right)$$

In the direct problem associated with the physical problem described above, the three thermal conductivity components k₁₁, k₂₂ and k₁₂ = k₂₁, the heat capacity per unit volume, as well as initial and boundary conditions, are known. The objective of the direct problem is then to determine the transient temperaturefield T(x, y, t) on the plate.

Several well-establishment approaches could be used to solve the direct problem. We have solved this problem here using both a numerical FEM and analytic integral transform technique.

Fig. 1 Geometry of the plate.

2.2 Finite element modeling of heat conduction in 2D

The temperature fields are computed simultaneously using a finite element formulation [19]. The region of interest is divided into an n^e number of three-noded elements Ω^e, with boundary ∂Ω^e with shape functions N^e(x, y) associated with each element e. In the isoparametric elements, the shape functions are used to transform the coordinates [19,20], and this enables better representation of any curved boundaries which may be present in the problem.

Particularly the FEM consists in a trivial choice of the decomposition functions $N_i^e(x,y)$ as linear piecewise functions associated to the node i of the element e. In this paper, we choose a triangular element with linear transformation functions: $$N^e(x,y)=\left[\begin{array}{lll}{N}_1^e(x,y)\hfill & N_2^e(x,y)\hfill & N_3^e(x,y)\hfill \end{array}\right].$$

Let n_t and n^e be respectively the number of nodes and the number of elements that contains the mesh. We discretize the heat conduction equations within each element as $$T^e(x,y,t)={\sum }_{i=1}^3{N}_i^e(x,y){\theta }_i^e(t)$$(3)

where θ^e is the nodal temperature vector of the element.

The temperature function T(x, y, t) is approximated in the solution domain at any time by the vector θ(t) with $${\theta }^e=A^e\theta$$(4)

where A^e is n^e × n_t matrix with each term $A_{ij}^e$ is equal to 1 if the node numbered (globally) j on the structure coincides with the node numbered (locally) i on the element e, and 0 otherwise and where θ_i (t) is the time-dependent value at node i. For each element domain Ω^e, the substitution of the equation (3) into the equation (1) and the application of the Galerkin method [20] results in the equation $$K^e{\theta }^e+C^e\frac{d}{dt}{\theta }^e=f^e$$

where K^e, C^e and f^e are respectively the conductance matrix, the capacitance matrix and the heat load vector per element: $$\begin{array}{l}{K}^e={\int }_{\Omega }{(▿N^e)}^{T}k(▿N^e)dxdy;\hfill \\ C^e={\int }_{\Omega }\rho C{(N^e)}^T(N^e)dxdy;f^e={\int }_{\partial \Omega }q{N}^{e}d\Gamma \hfill \end{array}$$

and the exponent T denotes the transpose of a matrix or a vector.

Let us denote by B^e the elementary gradient matrix: $$B^e={\displaystyle {\int }_{{\Omega }^e}\nabla }{N}^{e}dxdy.$$

Then the conductance elementary matrix becomes $$K^e=B^{eT}k{B}^e.$$

The assemblage system equation takes the matrix form $$C\dot{\theta }+K\theta =F$$(5)

where the conductance matrix, the capacitance matrix and the heat load vector are, respectively: $$\begin{array}{l}K=\sum _{e=1}^{n^e}{(A^e)}^T{K}^e{A}^{e}C=\sum _{e=1}^{n^e}{(A^e)}^T{C}^e{A}^e\hfill \\ F=\sum _{e=1}^{n^e}{(A^e)}^T{f}^e.\hfill \end{array}$$

The initial condition for solving the differential equation (5) follows from the equation (2d) $$\theta (0)={\left[{\theta }_0(0),{\theta }_1(0)\cdots ,{\theta }_{n_t}(0)\right]}^T.$$(6)

3 Inverse problem

Let us consider the ordinary differential equation $${\sum }_e{\left(B^e{A}^e\right)}^{T}k{B}^e{A}^e\theta +\rho cC\dot{\theta }=F.$$(7)

Concerning the inverse problem, the unknowns to be determined are the thermal conductivity parameters k₁₁, k₂₂ and k₁₂, and the heat capacity per unit volume ρc. The other values are known.

Least squares functional

Let us assume that the temperature fields are recorded during the time t_f.

Consider also the following sum as squares functional $$\mathcal{J}(k,\rho c)=\frac{12}{}{\displaystyle {\int }_0^{t_f}{\Vert {\sum }_e{(B^e{A}^e)}^{T}k{B}^e{A}^e\theta +\rho cC\dot{\theta }-F\Vert }^2}dt$$(8)

where ∥.∥ is the norm associated to the scalar product $\langle .,.\rangle$ in ${ℝ}^{n+1}$.

Our inverse problem aims the minimization of the functional $\mathcal{J}$ under the constraint k is symmetric positive definite.

Minimization of the least squares norm (8) is done by using a Polak–Ribire conjugate gradient algorithm [21]. The main principle of the technique consists in calculating at each iteration a size and a direction of the step minimizing the objective function. The step direction is a linear combination of the descent directions of the gradient vector provided by the previous iterations.

The conjugate gradient method with an appropriate stopping criterion belongs to the class of iterative regularization techniques, in which the number of iterations is chosen so that stable solutions are obtained for the inverse problem [22–24].

For this method, we need the gradient expression $\mathcal{J}$. Due to the analytical expression of the solution of the direct problem given by the equation (7), we can also obtain the analytic gradient expressions for the functional $\mathcal{J}$.

To calculate these gradients, we calculate the variations of the functional $\mathcal{J}$ respect to each variable k and ρc.

The variation respect to k reads: $$\begin{array}{l}\mathcal{J}(k+\delta k,\rho c)-\mathcal{J}(k,\rho c)\\ ={\displaystyle {\int }_0^{t_f}\langle {\sum }_e{\left(B^e{A}^e\right)}^T\delta k{B}^e{A}^e\theta ,u\rangle }dt\\ ={\displaystyle {\int }_0^{t_f}\langle \delta k,{\sum }_e{B}^e{A}^{e}u{\theta }^T{\left(B^e{A}^e\right)}^T\rangle }dt,\end{array}$$

where $u={\sum }_e{\left(B^e{A}^e\right)}^{T}k{B}^e{A}^e\theta +\rho cC\dot{\theta }-F$.

The variation of the functional $\mathcal{J}$ respect to ρc leads to: $$\mathcal{J}(k,\rho c+\delta (\rho c))-\mathcal{J}(k,\rho c)=\delta (\rho c){\displaystyle {\int }_0^{t_f}\langle C\dot{\theta },u\rangle }dt.$$

Since δk is symmetric, the gradient of the function $\mathcal{J}(k,\rho c)$ versus k is equal to the matrix zeros implies: $${\displaystyle {\int }_0^{t_f}{\sum }_e}{B}^e{A}^e\left(u{\theta }^T+\theta u^T\right){\left(B^e{A}^e\right)}^{T}dt=0.$$

The thermal conductivity symmetric tensor k and heat capacity per unit volume ρc are evaluated by solving the linear system: $$\begin{array}{lll}{\displaystyle {\int }_0^{t_f}\langle C\dot{\theta },u\rangle }dt\hfill & =\hfill & 0\hfill \\ {\displaystyle {\int }_0^{t_f}{\sum }_e}{B}^e{A}^e\left(u{\theta }^T+\theta u^T\right){\left(B^e{A}^e\right)}^{T}dt\hfill & =\hfill & 0\hfill \end{array}$$

4 Results and discussions

4.1 Temperature field simulation

Two programs have been developed for the direct and the inverse problems.

Let us consider a composite material constituted by oriented epoxy fibers in methyl methacrylate matrix. The preferential orientation of the fibers leads to an orthotropic material given two conductivity values. If fiber orientations are along the axis of the analysis, the thermal conductivity tensor is diagonal.

The dimensions of the sample are: l = L = 36.6 mm.

The diagonal thermal conductivity tensor components are: $$k_{11}=0.44\text{W/m}{/}^{°}\text{C}\text{and}{k}_{22}=0.19{\text{W/m/}}^{°}\text{C}.$$

The heat capacity per unit volume is $$\rho c=1465550{\text{J/kg/}}^{°}\text{C}.$$

The initial temperature of the medium is defined by T₀ =20 °C and the applied heat is ux q = 396.99 W/m².

We integrate according to the time the ordinary differential equation (5) with the initial condition equation (6) to simulate the temperature field for the estimation of thermal parameters.

A random white noise ωσ was added on the temperature field to mimic experimental conditions: $${\theta }_i={\theta }_i+\omega \sigma$$

where θ_i is determined from the solution of the direct problem. σ gives the noise range and the random-number ω is generated in the interval [−1, 1].

4.2 Identification of the thermal parameters

In this section, the estimation of the thermal conductivity tensor k and the heat capacity per unit volume ρc is carried out using the proposed numerical approach given in equations (9).

The accuracy and robustness of the proposed identification method are assessed by the help of numerical simulations using the direct problem (5) and (6). The relative error ɛ is then obtained by comparing the calculated result α_est est with the imposed α_imp: $$\varepsilon =\left|\frac{{\alpha }_{\text{imp}}-{\alpha }_{\text{est}}}{{\alpha }_{\text{imp}}}\right|.$$

The estimation of the accuracy is done according to the spatial steps of the temperature fields: $$dx=\frac{L}{n_x^e-1};dy=\frac{L}{n_y^e-1}\text{and}dt=\frac{t_f}{nb}$$

where $n_x^e$ and $n_y^e$ are the number of elements according to respectively x and y directions and nb is the number of images recorded during the test.

4.2.1 Analysis with a fiber orientation along the axis analysis

The composite material and its dimensions are defined in previous section. First, we consider an analyze in the principal axis of the material then the thermal conductivity tensor k components are still: $$k_{11}=0.44{\text{W/m/}}^{°}\text{C};k_{22}=0.19{\text{W/m/}}^{°}\text{C};$$

and the cross terms of the thermal conductivity tensor k₁₂ and k₂₁ are null.

We show the relative errors for given dimensions of sample on (Fig. 2) and one experience duration with several spatial steps and time steps.

Without noise on the simulated temperature, we have a good agreement on the estimated parameters about 10⁻⁴ (Fig. 2) and accuracy seems to be insensitive to spatial and temporal parameters except for large number of elements. Indeed, in this case, the relative errors are mainly due to the numerical calculation, the same equations and same parameters (time and spatial steps) are used in direct and inverse problems. Then, this procedure can be used only to validate the implementations and to give the ultimate error of the inverse algorithm versus time and spatial steps.

The Figure 3, obtained for a noise level of σ = 0.01 °C, shows that the noise involves a direct effect on the results. Error values are multiplied by 10. We observe an evolution of the accuracy versus the parameters. Results can be improved by increasing either the spatial steps or the time steps. We can also see that the thermal conductivitytensor components depend mainly on the spatial steps while the heat capacity per volume unit depends on the time steps.

Let us analyze now the influence of the noise level on the thermal estimated parameters for these steps. Three noise levels have been studied (σ = 0.00, σ =0.01, σ =0.03 and σ =0.06) which represent the common noise levels obtained with current infra-red cameras (Fig. 4). With noisy images, a filtering of the thermal field is necessary to improve the accuracy. We have tested several procedures, the best one is a spatial median filter (3 × 3 pixels) coupled with an average temporal filter (on 3 images). The efficiency of the spatial median filter is linked to a high spatial resolution that is why we have studied the noise effects for dx = 0.032 mm. This choice involves a time step of dt = 4 s to have a good accuracy.

We can observe the good precision on the estimated parameters on the filtered temperature field. For the high standard deviation σ = 0.06, the relative error on ρc is 0.054, k₁₁ is 0.018 and on k₂₂ is 0.082.

Fig. 2 Relative errors on the estimated parameters without noise on temperature fields.

Fig. 3 Relative errors on the estimated parameters with noise on temperature fields: σ = 0.01 °C.

Fig. 4 Relative errors on the estimated parameters versus the noise standard deviation.

4.2.2 Analysis with a fiber orientation out of the axis analysis

We choose again the same composite material epoxy-methyl methacrylate. The fiber orientation versus x-axis of the analysis is 30°. Then the thermal conductivity tensor components become: $$\begin{array}{l}{k}_{11}=0.3774\text{W}/\text{m}/°\text{C};k_{22}=0.2525\text{W}/\text{m}°\text{C}\text{and}\hfill \\ k_{12}=k_{12}=0.1083W/m/°\text{C}.\hfill \end{array}$$

The error of diagonal components is similar to the one obtained in the previous section (Fig. 5). Nevertheless, we note a low uncertainty for the thermal conductivity tensor component k₁₂.

To avoid this phenomenon, we propose to determine the thermal parameters along the principal axes of the composite material. This is possible because for this type of composite material, the angle between each fiber can be measured or defined during its elaboration.

Here, Ox and Oy are the axes of the Cartesian coordinate system, Oξ and Oη are the principal axes of the thermal conductivity tensor and ϕ is the angle between Ox and Oξ.

We denote the thermal conductivity tensor k′ in the principal axes. The relation between k′ and k is: $$k=P{k}^{\prime }{P}^T\text{ with\hspace{0.17em}}P=\left(\begin{array}{ll}\cos \varphi \hfill & -\sin \varphi \hfill \\ \sin \varphi \hfill & \cos \varphi \hfill \end{array}\right).$$

This expression is included in equation (8) to solve the inverse problem. Expression (8) is used again to obtain k′. The relative errors given by this new process are shown in Figure 6.

This new approach slightly disrupts the diagonal components of the thermal conductivity tensor and the heat capacity per unit volume, but it substantially improves the components k₁₂ = k₂₁.

Fig. 5 Relative errors on the thermal conductivity coefficients versus the noise standard deviation.

Fig. 6 Relative errors on the thermal conductivity coefficients versus the noise standard deviation.

5 Conclusion

The FEM has been applied to reduce the partial differential equations of heat conduction to an ordinary differential equation system. Numerical integration of this equation system led to the simulation of the temperature field in two dimensions.

The identification algorithm based on the conjugate gradient projected allowed us to identify with a very good accuracy the coefficient of specific heat per unit volume and the components of the tensor of thermal conductivity of homogeneous orthotropic materials. This algorithm provides a simple, rapid and systematic identification of thermal parameters of materials from temperaturefield measurements.

To identify thermal parameters on orthotropic material, the uncertainties are better if the procedure is used according to principal axis of the material.

We have studied a well-known epoxy-carbon material, other types of more or less conductive material can be studied by the same procedure. This work allows us to study the thermomechanical behavior of composite to establish an energy balance.

Acknowledgements

This work was supported in part through the Agence Universitaire de la Francophonie (BAO-2011-06-U-52210FT108) and the Poitou-Charentes region (RPC-R-022).

References

All Figures

	Fig. 1 Geometry of the plate.
In the text

	Fig. 2 Relative errors on the estimated parameters without noise on temperature fields.
In the text

	Fig. 3 Relative errors on the estimated parameters with noise on temperature fields: σ = 0.01 °C.
In the text

	Fig. 4 Relative errors on the estimated parameters versus the noise standard deviation.
In the text

	Fig. 5 Relative errors on the thermal conductivity coefficients versus the noise standard deviation.
In the text

	Fig. 6 Relative errors on the thermal conductivity coefficients versus the noise standard deviation.
In the text