Numerical Validation of Heat Conduction in 2D Binary Granular Mixtures under Mechanical Loading

This study set out to investigate numerical simulations of heat conduction to validate existing experimental results of two-dimensional binary granular mixtures subjected to mechanical loading. Data for this study were collected using molecular dynamics (henceforth MD) method. Three different configurations of the granular composite samples were systematically prepared under similar experimental conditions. A confined mechanical loading was applied to the granular samples. The fields of normalized temperature change of each particle were plotted for an individual sample. The results were statistically analyzed under a static equilibrium condition. The results indicate that simulation is in good correlation with the experiments in terms of statistical analysis via the probability of the distributions of the normalized temperature change. It also revealed that the normalized temperature changes which are greater than the average temperature distributes as an exponential decreasing for all tested samples. This study is in line with other studies that are related to force distribution of law. Less than 50% of particle numbers that have the normalized temperature changes, which are greater than the average value, is also explored. Besides, localizations of the temperature were found in the individual sample.


Introduction
Granular materials are one of the most widely used materials in our daily life and many industrial fields. For example, sand, rocks and soil are important components used in civil engineering or industrial field to prepare and build roads. Rice, crops and beans are fundamental to argo-food and our daily life. Some natural phenomenon is also involved in a granular movement like landslides and avalanches. In this study, granular materials are defined as a collection of solid particles whose macroscopic mechanical behavior is controlled by interparticle forces. They are composed of grains in various sizes, shapes and materials. This is a reason why behaviors of granular materials are generally complicated and it completely differs from the common solids and fluids [1,2]. The distribution of the interaction forces has a significant consequence in transport phenomena, especially in heat transfer. Besides, heat transfer plays a vital role in many industries regarding handling and processing of granular media, such as mixing, drying, granulation, and coating. Nowadays, numerical simulations based on the discrete element method (henceforth DEM) are an effective and useful tool for studying the effect of various parameters such as particle sizes, shapes [3][4][5]. Studies of DEM show the 10th TSME-International Conference on Mechanical Engineering (TSME-ICoME 2019) IOP Conf. Series: Materials Science and Engineering 886 (2020) 012050 IOP Publishing doi:10.1088/1757-899X/886/1/012050 2 importance of using DEM to simulate heat transfer processes in granular materials [6][7][8][9]. Previous research on the applications of granular flow i.e. rotating drum [7][8][9], fluidized systems [10][11][12] has been proposed. Previous studies of DEM have not dealt with the heat transfer phenomena in twodimensional granular composites under mechanical loading. Therefore, it enables us to validate our thermoelastic stress analysis (TSA) experiments recently performed in 2014 [13]. An objective of this present study is to statistically compare the simulation results with the existing TSA experimental results. Simulation configurations are prepared as the similar to what we done in the experiments. Molecular dynamics method is carried out for the method of simulations.

Molecular dynamic (MD) simulation
One study by Cundall (1971) examined the discrete element method (DEM) for rock mechanics [14]. Using this approach, researchers have been able to study the mechanical behavior of granular media. The MD method belongs to the DEM and relies on an explicit algorithm. Moreover, this method considers the particles as rigid bodies with non-conforming surfaces [15]. The MD method provides a mean of simplicity and flexibility because it applies only Newton's second law for calculating kinematical quantities: 1 ,2 , 3 , . . . , (1) where N is the number of particles in the simulation, m i is the mass of particles i, x i is the position of particles, and F i is the force exerted on that particle. The method consists of calculating force F i and solving ordinary differential equations (ODE) in equation (1). Integration of Newton's equations of motion can be achieved by applying a predictor-corrector scheme with Gear's set of corrector coefficients [16]. By using this method, positions, velocities, accelerations, the third-order time derivative of particles are provided.

Numerical models
The contact force between two dry particles can be separated into two components: a normal force and a tangential force. In the MD method, the deformation is simply modeled by using a virtual overlap  . Due to the unchanged particle shapes, the normal force can be calculated by a function of the virtual overlap, which is expressed in terms of vector position of the particle centers i x and j x and radius i r and j r : In this study, the "linear spring-dashpot" [17,18] is a common model to calculate the normal force between two particles. The normal force n F is calculated by: where eff k is the effective contact stiffness, n  is the normal damping coefficient, and n v is the normal velocity (the first-order time derivative of the virtual overlap). The integration of Newton's equations of motion requires a smooth friction law (mono-valued) i.e. the tangential force must be expressed as a linear function of the sliding velocity t v , while the classical Coulomb's law friction is a non-smooth function. Hence, the tangential force t F in this study can be commonly determined by using a regularized form of the exact Coulomb's law, which can be calculated by this expression: 10th TSME-International Conference on Mechanical Engineering (TSME-ICoME 2019) IOP Conf. Series: Materials Science and Engineering 886 (2020) 012050 IOP Publishing doi:10.1088/1757-899X/886/1/012050 3 where t  is the tangential viscosity coefficient and  is the coefficient of friction. In this study, the rotational motion that occurs from the tangential force is assumed to be free because its value is very small compared to the normal force. The phenomena of heat transfer consist of heat conduction, heat convection, heat radiation, and heat generation, which are complex and depend on many variables. Therefore, the phenomena of heat transfer in the present study are considered only effects of heat conduction and heat generation. When the particles inside the granular system move and two bodies are in contact, their relative motion can be taken into account as a source of heat generation, as denoted as s ij Q . The heat generation can be computed in many different methods depending on the nature of the considered system. The heat generation model used in this study is directly based on the first principle of thermodynamics, accounting for normal motions as well as tangential motions. The heat generation of each particle, which is related to an energy dissipation at the contact level [19] can be expressed as: Then, the heat conduction of particle i is the sum of all heat exchanges of contacting particles [6] which can be expressed as: where c ij H is the heat conductance of contact and the term of i j T T  is the temperature difference between the particles which relies on two assumptions. The first assumption is no temperature distribution inside the particle i.e. each particle has only one temperature for the whole body. The second assumption is that the temperature exchange between the particles is slow enough during each time step. The coefficient c ij H , which is a function of compression force, depends on the particle properties. In this study, the particle is assumed as a cylinder similar to the experiments. The solution was then assayed for coefficient c ij H c using Hertz theory, as shown in [18]: where eff R , eff E , s eff k and L represent the effective contact radius, the effective Young's modulus, the effective thermal conductivity and the length of the cylinder, respectively. For the evolution of the thermal system, the Euler scheme is applied to discretize the differential equation characterizing the evolution of the temperature of each particle. Thus, the temperature at each time step dt can be discretized explicitly using a -method over time t to time t+dt [19], which can be expressed as: where indicates a value between 0.5 to 1 for stability reasons. In order to reach the conservation of the simulation, the value of equals to 0.5 was used in this study. The term t dt i T   is a function of the local thermal phenomena of particle i by accumulating from multi-contacts [6]. It can be thus calculated from: 10th TSME-International Conference on Mechanical Engineering (TSME-ICoME where i  , i c and i V is the bulk density, the specific heat capacity and the volume of the particle. In 2D simulations, the particle volume is replaced by an area of the particle. It must be noted that the quantities i  , i c and i V is constant during the simulations. The term t dt i Q  in equation (9) implies that the amount of heat is accumulated from all thermal phenomena of particle i at time t+dt. The total amount of heat in each particle can be thus separated into two components: where c i Q represents the heat exchange (heat conduction) between particle i and its neighbors, s i Q corresponds to the heat generation by the relative motion between particle i and its neighbors.

Numerical preparation
Three granular composite configurations in the TSA experiments [13] were chosen to be prepared as the two-dimensional numerical composite samples. The granular systems were made of two different constitutive materials: one is polyoxymethylene (POM) called a "stiff" particle. Another one is highdensity polyethylene (HDPE) called a "soft" particle. It must be noted that the stiff particles are approximately four times stiffer than the soft particles [13]. Table 1 presents configurations of all tested samples. Sample #1 is a monodisperse system that has a diameter ratio is equal to one Whereas others are bidisperse systems. Prior to the experiment, the locations of particles for each configuration were extracted from an optical image of a real composite granular system. After extracting the locations, the initial states were prepared for numerical simulations. Then, each particle was deposited inside a box consisting of four rigid walls. quasi-static conditions, a vertical compression force of 60 kN was incrementally applied to granular samples by the bottom wall. Whilst the others were stationary during the simulations. The gravitational force was also taken into account, although it can be negligible compared to the magnitude of the external loading. The normal and tangential forces between the particles are calculated by using Eq. (3) and Eq. (4), respectively. The contact force between particle and wall is also computed by these expressions. The effective contact stiffness (k eff ) depends on the type of contacting particles: the case of interparticle contact is defined by (k stiff ×k soft )/(k stiff +k soft ), while the case of particle-wall contacts is defined by the stiffness of such particle that contact to the wall.
After the contact force calculation, the contact forces and relative velocities were then employed to compute the heat generation inside the contacting particle. The particle which has a different temperature from its neighbors will exchange heat to one another. This phenomenon can be expressed in terms of the heat conduction and it depends mostly on the contact area and the effective thermal conductivity. Additionally, the initial temperature of all particles in all configurations is defined as 25°C. After the calculation of the local heat effect, this local heat effect was used to update the temperature of each particle in each time step. Finally, the simulations were accomplished when the granular systems reached to static equilibrium conditions.

Results and discussion
In order to improve data consistency within the database, normalized temperature change rather than the temperature is applied for the analysis. It is defined as the temperature change in the particle divided by the average temperature change of all the particles inside the granular sample.

Field of normalized temperature change
The normalized temperature change ΔT norm of all the particles was plotted with their positions to create a field of normalized temperature change. The fields of normalized temperature change obtained from the simulations and experiments for all configurations are illustrated in Figure 1. It can be seen that the field of normalized temperature change for sample #1 exhibits the triangular network (see Figure 1a). It is well-known for the monodisperse granular system [13]. This behavior is similar to what we have observed from the force or stress fields [13,20]. Considering samples #2 and #3 in Figures 1b and 1c, most of particles which change in their temperature has a lower temperature than the average temperature is stiff particles. This can be explained by the fact that both the mass and 10th TSME-International Conference on Mechanical Engineering (TSME-ICoME 2019) IOP Conf. Series: Materials Science and Engineering 886 (2020) 012050 IOP Publishing doi:10.1088/1757-899X/886/1/012050 6 the specific heat capacity of stiff particles are larger than those of the soft particles. It must be noted that this change in temperature does not imply to the amount of the heat storage inside the particles. Indeed, it is not possible to obtain precisely the same temperature pattern between simulations and experiments at the particle scale. This might be due to an intrinsic variation of parameters in the experiment, such as variations in terms of diameter size, roundness, physical material properties, etc. A technical problem regarding particle locations might be also included. For this reason, a statistical analysis of the experimental and simulation results is performed in the next section.

Probability distribution function of normalized temperature change P(ΔT norm )
In this section, a probability distribution function of the normalized temperature change was determined for each configuration. It must be noted that five additional simulations with randomly changes the particle locations were systematically performed for each sample: simulations A, B, C, D, and E, respectively. Figure 2 shows the P(ΔT norm ) of each configuration obtained from the simulations and experiments. , and c) Sample #3. Note that simulations "A", "B", "C", "D", and "E" were performed under the same experimental configuration, excluding the particle locations, while the simulation "SP" was conducted under all the same as in the experiments.
The results from the distributions of P(ΔT norm ) obtained by both simulations and experiments are in the same trend for all the samples, which exhibit as an exponential decay. This distribution is characterized by the coefficient β, which can be written as: (1 ) norm T P e    (11)