Capturing heat transfer for complex-shaped multibody contact problems, a new FDEM approach

This paper presents a new approach for the modelling of heat transfer in 3D discrete particle systems. Using a combined finite–discrete element (FDEM) method, the surface of contact is numerically computed when two discrete meshes of two solids experience a small overlap. Incoming heat flux and heat conduction inside and between solid bodies are linked. In traditional FEM (finite element method) or DEM (discrete element method) approaches, to model heat transfer across contacting bodies, the surface of contact is not directly reconstructed. The approach adopted here uses the number of surface elements from the penetrating boundary meshes to form a polygon of the intersection, resulting in a significant decrease in the mesh dependency of the method. Moreover, this new method is suitable for any sizes or shapes making up the particle system, and heat distribution across particles is an inherent feature of the model. This FDEM approach is validated against two models: a FEM model and a DEM pipe network model. In addition, a multi-particle heat transfer contact problem of complex-shaped particles is presented.


Introduction
Heat transfer in particle systems is encountered in a large number of applications such as chemical and nuclear engineering and in soil and fluid mechanics. Finite element methods (FEM) have been used in the modelling of microscopic contact roughness and/or mechanical interactions [14,24,27,31]. More advanced work included modelling of the contact roughness structure and deformations together with contact heat transfer [22]. However, these works aim to model mechanical and thermal interactions of a specific contact configuration rather than study of the overall heat distribution across multiple contacting solid bodies. A similar finite elements procedure has been presented [4,17], using boundary nodes of the contacting solids to determine heat flow through contact zones. However, the latter method is two-dimensional and has no intent or ability to model multiparticle systems. The powerful digitalised method from Jia et B Jiansheng Xiang jxiang@imperial.ac.uk 1 Department of Earth Science and Engineering, Imperial College London, London, UK al. in 2002 [16] enabled the use of arbitrary sizes and shapes in 2D with the ability to model large systems with a low computational effort; however, the pixelisation of particles and their boundaries presents consistent limitations when studying contact problems. Discrete element methods (DEM) are specific to the analysis of multi-bodies, relying on a low computational effort and enabling simulations with large sets of particles. Hence, research on multibody contact heat transfer with discrete elements has been very active in the past 25 years. A large range of publications treated as the heat conduction across randomly packed spheres, in 1990 Jagota et al. [15], established a method to determine the effective thermal conductivity of a packing of spheres. Following this direction, Argento and Bouvard [2] pioneered the use of FEM model for the deformation of spheres in contact. Hence, linkage between the contact area and the heat resistance is important. Cheng et al. modelled heat transfer in mono-sized spheres in the presence of a stagnant fluid [7]. Siu and Lee [30] presented in 2004 a three-dimensional model for randomly packed spheres. Feng et al. [9,10] developed a pipe network model to simulate heat transfer in large numbers of circular particles in 2D that represent infinite or long pipes. Rickelt et al. [25,26] developed a radial temperature model that allows to simulate granular systems of particles of different sizes and materials, enabling the use of DEM in various applications. Gan et al. [11] extended DEM to simulate contact heat transfer in packed and fluidized beds of ellipsoids. However, the above methods are restricted to the analysis of simple-shaped particles (spheres, discs, cylinders, ellipsoids, etc.) in past. Although DEM extends its capabilities to simulate complex shapes using superquadrics [1], polyhedral DEM for convex [13,29] and non-convex particle shapes [12], clustered sphere [8], etc., it still relies on the assumption that the temperature is uniform within each particle which is not valid in some configurations as this present paper will highlight. For more complex contact conditions, it is very difficult to provide an accurate solution using DEM. Finally, to compute the heat distribution within each particle, DEM models rely, when available, on extra analytical solutions [9,10,25,26]. Scherer et al. [28] concluded the best solution might be to use FDEM which can be used to describe particle deformation. But to the best knowledge of the authors, FDEM has not yet been applied to thermochemical problems within DEM. This paper presents a new approach for three-dimensional multibody contact heat transfer using the FDEM approach [20] where discrete particles are described with FEM mesh and interactions between them are treated with the DEM approach. The contact interactions between discrete bodies are calculated via the penalty function method [21] relying on the overlap between the contacting discrete meshes. Using this existing algorithm, the presented method evaluates the apparent contact surface area from the contact overlap of discrete meshes. Contact heat flux is calculated from the boundary temperatures of the solids involved in the contact. Moreover, this method possesses the advantage of being suited to simulating an accumulation of particles of any sizes and shapes for which the heat distribution within a particle may be an important feature. (This in contrast to DEM but comes at a higher computational cost.) This paper unfolds as follows: Firstly, the theory of standard heat conduction and contact heat transfer in the FDEM code is presented in Sect. 2 and 3 and then validated in Sect. 4; secondly, the FDEM contact model is compared to a DEM pipe network contact model [9,10] in Sect. 5. Finally, the code's ability to model complex configurations of particles is demonstrated in Sect. 6.

Problem Statement-This section presents modified equations from [23]
Let us consider an isotropic body with non-temperaturedependent heat transfer properties. A basic equation of heat transfer has the following form: With q x , q y , q z components of the heat flow through an unit area, Q is the inner heat generation rate per unit volume; ρ is the material density; c is the heat capacitance; T is the temperature field; and t is time. According to Fourier's law, the components of the heat flow can be expressed as follows: The thermal conductivity coefficient of the media is noted k, and the combination of these equations yield the following basic heat transfer equation: We then consider the following boundary conditions: -The initial temperature T (x, y, z, t = 0) = T 0 (x, y, z), -Specified temperature or Dirichlet boundary condition T s = T (x, y, z, t) on surface boundary Γ s , -Specified heat flux or Neumann boundary condition q h = −(q x n x + q y n y + q z n z ) on Γ h , -Contact heat flux q c on contact surface boundary Γ c , -Convection boundary condition q cv = h ΔT solid−env on surface boundary Γ cv . ΔT solid−env being the temperature difference between the solid and the temperature of the environment in which the solid is placed, h is the heat convection coefficient.

Finite element discretisation
The domain V is divided into finite elements connected at nodes; the boundary of the domain is noted Γ . Interpolation functions are used for calculation of temperature inside each finite element composed of n e nodes: T is the matrix of temperature distribution inside the finite element, {T } is the vector of temperatures at the nodes, and N is the interpolation or shape function matrix. In this study, we use four-noded tetrahedral elements, and the corresponding shape function is: With ζ n e the isoparametric coordinates so that ζ 1 + ζ 2 + ζ 3 + ζ 4 = 1. The spatial gradient of the shape function B is: Similarly to temperature, the source term is declined in nodal values: We have the following equation for a finite element of volume V : Applying the finite element discretisation yields: where i = 1, . . . , n e represent each set of equations in the finite element. Applying the divergence theorem yields: with where {n i } is an outer normal to the surface of the body. Also, it is worth noting that After insertion of boundary conditions into Eq. 10 and using Eq. 12, the discretised equations are as follows: Finally, we obtain the following condensed expression of the finite element model:

Thermal contact
The heat flux across the apparent surface area of two solids in contact is defined as follows: in which R c is the contact heat resistance and ΔT c the apparent temperature drop at the contact. This definition introduces a fictional apparent temperature drop at the interface. In reality, there is no real discontinuity of the temperature distribution through the solids' contacts. There is a continuous distribution of temperature extending through the contact interface from both solids. As shown in Fig. 1, defining the temperature drop as the difference in the temperature obtained by extrapolating the temperature profiles in the two regions of the interface enables the use of the contact heat resistance to simplify the complex heat processes occurring at the boundary.

Temporal integration
The differential Eq. 14 needs to be integrated with respect to time to obtain a transient solution of the heat equation. As defined in [18], the θ -family time integration methods are of the most commonly used: The term {Ṫ } n is obtained from solving Eq. 14. For θ = 0, we obtain the explicit forward difference scheme: For different θ > 0, the above equation will refer to mixed implicit-explicit methods. For θ = 1, we obtain the fully implicit backward difference (or backward Euler) method which is unconditionally stable, i.e. there is no restriction on the time step size. The theta method offers a large range of options for solving steady state and transient contact problems, as the final temperature distribution can be obtained in one iteration (θ = 1), and different transient states can be computed with accuracy (θ > 0).

Contact surface area
The presented method evaluates the contact area between two contacting solids based on the penetration of boundary meshes, see Fig. 2. This surface calculation method is based on the routines of the Y code developed by A. Munjiza [20]. A contact surface is obtained for both solids, and the overall contact surface is the average of these two values. The method selects a couple of contacting elements, one is called the contactor element, and the other one is called the target element. The boundary surface of the target element in contact with the whole contactor element's volume is calculated, and then the opposite calculation is performed. Each solid is meshed with four-noded tetrahedral elements, the algorithm loops on each face of the target element, and the intersection surface with the contactor's volume is drawn on each target's face (and vice versa), see Fig. 3a. Two surface areas are obtained, one describing the contact area on the target, the other on the contactor, which we call, respectively, S tar and S con , the contact area is set to be the average of these two. Note that only the faces of the target element located on the boundary of the solid can be selected for surface calculation. Figure 4 shows boundary elements from a first solid A contacting a boundary element from a second solid B. For purposes of explanation of the surface calculation, consider the blue element to be the target element and red elements to be selected successively as contactor elements; also note that the relative penetration size has been intentionally exaggerated on the figures. Then, consider that the target element from solid B possesses only one face located on the boundary of solid B; this face is highlighted in Fig. 3b. The algorithm will intersect this face with all three contactor element volumes from solid A in order to reconstruct the surface area of contact; therefore, three S tar contact areas are obtained, see Fig. 3c. At last, the opposite calculation is performed, and surfaces S con are obtained.

Computation of heat fluxes
Heat fluxes are calculated between each couples of contacting tetrahedral elements from two contacting solids meshes. The total contact heat flux temperature contribution of Eq. 14 is for a finite element: where n c the number of contacts for which node i is involved in. The contribution is equally distributed between the three surface nodes of each contacting tetrahedron; hence, we have a 1 3 factor. S k c is the contact surface of the k contact, being the average between the contactor's and the target's overlapping boundary surfaces:  Figure 5 shows the contact heat transfer for two isolated contacting tetrahedral elements. The nodal heat flux contribution for the contacting surfaces (T 12 , T 13 , T 14 ) and (T 21 , T 22 , T 24 ) is the following for nodes T 12 , T 13 and T 14 : And for nodes T 21 , T 22 and T 24 : We can summarise the method with the following expression: where T k con,l and T k tar,l the nodal temperatures of the target and the contactor elements involved in the k contact.

Contact interaction matrix
Pursuing the objective of solving large complex contact heat transfer problems in a minimum time, the PETSc tool-kit for scientific computation is solving the general heat transfer differential Eq. 14 which was implemented in our program. In order to write Eq. 28 in a matrix form, we define the contact interaction term B c i as follows: where C h is the contact heat transfer interaction [n e ; n e ] sparse matrix containing node interactions between elements where i, j ∈ [1, 2, . . . , n e ] nodal indexes of the target and contactor elements of the considered contact couples.

Calculation procedure
The presented thermal contact model is integrated into the FDEM framework and follows the following procedure:

Algorithm limitations
This algorithm excels in capturing the surface area of contact as the intersection polygon is drawn across the elements faces; however, in the event of modelling solids with complex shapes and curves, the surface area calculation can only be as good as the mesh approximation of the solid's boundary.
In addition, the surface of contact is linked to the amount of mesh penetration which depends on the penalty function. Depending on the configuration of the contact (edge to edge, edge to surface or surface to surface), the penalty function may or may not influence the surface of contact in a significant manner, and further research on this topic is necessary. Nevertheless, the method remains conservative; the contact heat flux distribution is bound to be equal and opposite between two solids as the same quantity of heat is exchanged over the same area.

FEM and FDEM perfect contact conduction validation
In this section, we validate the finite element conduction heat transfer through one continuous solid and the contact heat transfer across two contacting solids with a perfect contact.
Results are compared with a one-dimensional analytical solution for conduction in solids.

Analytical solution
Consider a finite slab of length L, the solid is of an initial temperature of T 0 . The left side of the slab is insulated, while the right side is exposed to a Dirichlet boundary condition of an imposed temperature T D , see Fig. 6. There is no inner heat generation in the slab. The one-dimensional transient conduction equation for this problem is: where α the thermal diffusivity, k the thermal conductivity and c the thermal capacity. The solution given in [5] is of the form:

Perfect contact validation
For the FEM simulation, we build a 3D bar of a length of 1 m, see Fig. 7a and b. The solid bar has an initial temperature of 0 • C, and a temperature of 1 • C is imposed at the right-hand end face of the bar, and the rest of the faces are insulated. The FDEM simulation is a composition of two 0.5 m bars contacting at one end, and the amount mesh of penetration is of 0.001m, see Fig. 7c and d. The solid bars have an initial  Fig. 8. To simulate a perfect contact, the heat resistance is set to a relatively low value (Table 1) . Simulations are performed with two different mesh sizes, 5.10 −2 m and 1.10 −1 m. Results show a very good agreement and are presented in Figs. 9, 10 and Table 2. It is worth noting that the computational costs are relatively high for the number of finite elements considered; this is undoubtedly due to the fact that a matrix system of equations is solved at each of the 200,000 time increments. This numerical scheme will

Model description
The pipe network model presented by Y.T. Feng [9,10] is designed for the modelling of large numbers of circular particles in 2D as it would represent long or infinite pipes. This method is presented in the culture of the discrete element method and is introduced here to form the basis of a validation study of the new FDEM method. Consider two circular Fig. 11 Contact heat flux for the discrete thermal element, modified from [9] particles A and B having respectively T A and T B as average temperatures; the thermal resistances of the two pipes are, respectively, R A and R B . The total thermal resistance is: where R * c the contact thermal resistance for the pipe network model. The contact zone of the discrete thermal element is represented by an arc on the boundary of the element which is defined with its half angle α i , see Fig. 11. For angles of contact below 30 • , the discrete element thermal resistance can be approximated with high accuracy by the formula: where k A , the thermal conductivity of particle A. The boundaries of the particles are insulated and heat transfers only through the contact zone. For two contacting discrete thermal element (see Fig. 12), the heat flow between the two particles Q AB is defined as follows :

Transient analysis
The forward difference explicit time integration will be used to solve the transient problem. For this test case, the same method is employed for the FDEM contact algorithm (θ = 0, see Eq. 23); therefore, T A being the time derivative of the average temperature and C A the total heat capacity of the particle: where c p the heat capacity and r A the particle's radius.

Model adjustment
As demonstrated by Eqs. 21, 14 and 35, 36, the FDEM and the pipe network model approaches for the heat resistance differ, and therefore an adjustment is required. The same contact heat flux contribution needs to be taken into account the heat Table 3 Simulation parameters for the FDEM two-particle configuration diffusion equation; therefore, the following condition has to be fulfilled: The left-hand side of the above equation is the heat flux contribution from the pipe network model extended from a 2D disc to an hypothetical 3D cylinder of a width represented by w cyl ; the right-hand side corresponds to the FDEM model. Therefore, we write: We can also write: This conclusion also implies that T B − T A = ΔT c , i.e. the particle average temperature difference is equal to the local temperature difference at the contact zone and such is the main approximation of the discrete element approach; this condition will only be verified when the thermal conductivities are high compared to the contact heat resistance. To make sure this assumption is acceptable in the following simulation, in addition to the ΔT c local temperature gap calculation, a ΔT c calculation based on the particle average temperature difference has also been implemented in FDEM, the code to compare with the pipe network model.

FDEM contact simulation settings
In order to make possible a comparison between the simulation of the 2D problem of two contacting discs and the new 3D contact heat transfer FDEM, we define two contacting thin cylinders of the same radius r and width w cyl . The two finite element meshes are overlapping at the contact zone, see Fig. 13. Simulation parameters are summarised in Table 3. Two different FDEM simulations were performed, the first with a contact heat flux calculated with the local temperatures, the second calculated with the average temperatures.
The accuracy of the computed contact surface S c is validated against a theoretical surface formula S th c obtained from the overlap of two circles: where d = r A + r B − p, p being the penetration of the two meshes. For the actual configuration, the theoretical surface is S th c = 0.02 m 2 , and the error of the computed contact surface is therefore of 4%. Again that error is only due to the finite element approximation of the domain. Contact surface error reduces to 0.4% with a two-time smaller mesh. Nevertheless, to reduce errors for this validation test, the computed contact surface is transformed into the equivalent contact half angle and then imputed in the pipe network model by means of the formula: The two particles with different temperature are in contact initially. Figure 14 shows the evolution of the temperature changes against time. Results of the temperature are presented for particle B in Fig. 14. We obtain an average error of 7.6 % between the pipe network model and the original FDEM (with a local temperature difference) against an average error of 0.63 % with FDEM with the particle average temperature difference. In this configuration, it is evident that the main assumption of the pipe network model, being the consideration of the particle's temperature to be uniform, is not valid as the heat conductivity is relatively small (Fig. 15). We consider now the same configuration with a heat conductivity a hundred times greater, and results are presented in Fig. 16. In this case, both FDEM simulations produce the same result with an average error of 0.7 %; hence, the pipe network approximation can here be considered valid.
In summary, the foregoing examples show that when significantly varying temperatures exist in the contacting bodies, the FDEM code can capture the complexity of the time history, giving quite different results when the assumption of average temperature within the particle is imposed.

Multibody simulation
This section presents a static multibody heat transfer simulation with 2000 particles packed up into a cylindrical container which is heated from the outside. To obtain such a configuration, the particles are dropped into the container using FDEM. After all the particles are deposited, a fixed temperature boundary condition of 300 • C is applied to the container. Note that for this simulation, only heat conduction heat transfer is considered. Moreover, the system is static during heating because thermal expansion effects are not taken into account. The deposited particles are in the shape of a metal nut, and their mesh is composed of 192 elements (see Fig. 17). Simulation properties are listed in Table 4. Finally, Figs. 18 and 19 show the heat transfer simulation with the nut particles being heated up from 0 to 300 • C (Fig. 20).

Conclusion
A new 3D approach is presented in this paper, first with a onedimensional heat conduction validation test that has proven to give very accurate results. However, we are aware that this simulation was for a perfect contact; thus, a very low contact resistivity was chosen. For this reason, the temperatures on both sides of the contact zone are almost the same and we have a temperature continuity; thus, FEM and FDEM produce a very similar result. This first validation case is strengthened with a second successful one which is a realistic contact problem as a higher thermal resistivity is set. The FDEM approach has offered here more flexibility than a discrete element approach as it can consider cases where the average temperature is not representative of the problem. This means FDEM can model problems with extreme boundary conditions and complex particle heat distributions with accuracy. Finally, as the multibody simulation demonstrates, FDEM can handle large sets of particles and compute the thermal distribution at desired times. The strength of FDEM is also that it can be coupled with existing fracture, plasticity and thermal expansion models. Hence, it is critical for the future of this method to provide a more detailed understanding of the linkage between contact force and heat resistance which has been so far been neglected in this work and requires further theory and method for implementation.

Future work
As the FDEM method presented here is already capable of tackling multibody dynamic problems, future work will include such simulations with heat transfer. We have presented here the process to capture the contact area of two penetrating meshes; however, it is dependent on the penalty function method [21] and the amount of mesh penetration. Hence, it is important to investigate the sensitivity of the surface calculation to the mesh penetration. The role of the penalty function is to calculate element-to-element intersection for the contact force. Therefore, it will be powerful to insert in this method an empirical law linking the element-toelement intersection to the heat resistance and contact surface such as in Bahrami et al. works [3]. Finally, the contact heat transfer method will be validated further in more complex configurations and with experimental results. To validate the contact heat transfer, suitable experiments can be selected from those performed in a vacuum [6] or where conduction heat transfer is largely dominant [32].