Numerical modeling of heat flow between three-dimensional regions in contact problems

The method of numerical modeling of heat transfer between threedimensional objects being in contact is described in the paper. Presented approach is based on the finite element method (FEM) with independent spatial discretization of considered regions. The gap between external surfaces of the interacting objects has variable width and is filled with gas or liquid. The medium in the gap introduces thermal resistance into heat transfer process. The mathematical model of considered problem is based on the equation of heat diffusion supplemented by the appropriate initial and boundary conditions. The deformations of the regions resulting from the thermally dependent changes of their volumes are also included in the model. The results of numerical simulations are presented and discussed.


Introduction
Machine parts during the work come into thermal and mechanical contact. Mechanical contact is the result of the interactions of their external boundaries. These interactions may be caused by motion or thermal deformations. The ideal thermal contact occurs when the width of the gap equals zero. If the gap width is greater than zero the contact with thermal resistance is observed. The mathematical model of the heat transfer process is based on the heat conductivity equation with appropriate boundary and initial conditions. A fourth type of boundary condition plays very important role in the mathematical model of such process. One can distinguish two kinds of such mathematical condition describing the ideal contact or the contact with the presence of the gap of variable width, filled with gaseous or liquid medium.
Mathematical descriptions of the non-ideal thermal contact was described for the first time in [1]. Analytical solutions of the problem were presented in [2]. Numerical solution of heat transport using finite difference method (FDM) was shown in [3]. Nowadays many scientists from different branches of science investigate the process of the heat transport through the gap. For example the mathematical description of the considered phenomenon is required in the numerical analysis of the solidification in the mold [4,5]. Currently many advanced numerical methods such as the FEM are used to model complex mechanical phenomena [6][7][8][9]. Numerical investigations based on the FEM in the case of air gap formation process during solidification of pure metal was discussed in [10].
Presented mathematical model also contains the description of thermal deformations of considered regions which are needed to calculation of the local width of the gap. Fig.1 presents two objects Ω 1 and Ω 2 in contact. Boundaries Γ 1 , Γ 2 lie close together forming a gap of width h. The width of the gap depends on the coordinates x, y, z and also may change according to time. At the boundaries Γ 5 , Γ 6 temperature is defined while Γ 3 , Γ 4 are fixed. The process of heat transfer in Ω 1 and Ω 2 is transient and the objects are subjected to the thermal deformations.
Equation (1) is supplemented by the following boundary and initial conditions:  heat flux at the contact boundaries  initial temperature in the volumes 1 and 2 where λ g denotes thermal conductivity of the medium filling the gap, h is the local width of the gap, n 1 , n 2 are the local directions of the vectors normal to the boundaries Γ 1 , (2) -the local temperatures of the volumes Ω 1 , Ω 2 at the contact boundaries, T 0 , T 0 (2)initial temperature of the volumes Ω 1 , Ω 2 .
There are 6 unknowns in the equations (2), so they need to be transformed into displacement-dependent functions. This is done by using the following relations:  relations between the stress and strain tensor components  relations between the strain tensor components and the displacements where: ε x , ε y , ε z , γ xy , γ xz , γ yz are the components of symmetrical strain tensor, u x , u y , u zthe components of displacement vector, E -Young's modulus, ν -Poisson's ratio, α -the linear thermal expansion coefficient, ΔT -the difference between the reference and the current temperatures ΔT=T-T ref .
Using relations (7) in equations (6) and then inserting the stress components into equilibrium equations (2), the following equations are obtained: where f 1 , f 2 , f 3 , f 4 are defined in the following way: Equations (8) are supplemented by the following boundary conditions:

Numerical model
The numerical model of presented problem is based on the finite element method [11] and derived from the criterion of the method of weighted residuals [12]. Equations (1) and (8) are multiplied by the weighting function w and integrated over the volume Ω i . Equations (1) and (8) written in the weak form are shown below: where K T is the global thermal conductivity matrix, M T -the global thermal capacity matrix, B T -global vector associated with the thermal boundary conditions, T -vector of unknown nodal temperatures, f -time level, K D -the global stiffness matrix, B D -vector where known boundary displacements are stored, u -vector of unknown nodal displacements.
The solution of equation (15) is obtained at every step through an iterative process, independently in each region. Calculated temperatures are then used to find thermal deformations of each object. It is crucial to properly incorporate the condition (3) to the heat transport solution process [13]. During the deformations of Ω 1 and Ω 2 relative positions of objects and the width of the gap between them are changing. Because the global set of equations is built and solved independently for each region the presented approach is based on the disconnected finite element meshes.

Example of calculation
Two cuboids of the same dimensions 0.02x0.1x0.1 m are considered (Fig. 2). Object on the left is made of steel and its initial temperature T 0 =600 K. Object on the right is initially colder, its temperature T 0 =300 K and the material from which its made is copper. Both cuboids are fixed at their bottoms, u x = u y = u z = 0 m. Between them is a gap of initial width h=0.1 mm and thermal conductivity λ g =0.5 J•s -1 •m -1 •K -1 . According to the boundary and initial conditions heat transfer between the objects and their thermal deformations are simulated with the use of an original computer program written in C++. Finite element meshes are made with the use of mesh generator GMSH 2.9. Left cuboid is divided into 49498 tetrahedrons with 11468 nodes, the right one has 50315 tetrahedral elements and 11602 nodes. Time step used in the calculation process Δt = 0.1 s, total time of the analysis t total = 17 s. Fig. 3a-d show finite element meshes deformed due to temperature after 1, 5, 10, 17 s. The gap is widest in the upper part of the cuboids, it reaches the highest value 0f 0.7 mm after 5 s (Fig3.b). Fig. 4a-d show instantaneous temperature distributions in the steel cuboid. The impact of the local gap width on the temperature distribution on the wall in contact with the copper cuboid is evident. In both figures the deformations are exaggerated (def_factor = 20).

Conclusions
The numerical method of heat transfer during the contact of two three-dimensional objects taking into account thermal deformations was presented. The method using the independent spatial discretization for each volume saves a lot of operational memory because of sequential solving of the sets of equations. Obtained results proves significant impact of the variable width of the gap on the global heat transfer. Presented approach may be useful in computer simulations of various processes such as the solidification of the molten metal in the permanent mold. Further work will be focused on the expansion of the numerical model by the mechanical interaction between the objects.