Effective Thermal Expansion Property of Consolidated Granular Materials

Thermally-assisted compaction of granular materials is of keen interest in many engineering applications. A proper estimation of the material behavior of compacted granular materials is contingent upon the knowledge of microstructure formation, which is highly dependent on the bulk material properties and processing conditions, during the deformation stage. Originating from the pair interactions between particles, the macroscopic properties are obtained using various homogenization techniques and postulating continuum constitutive laws. While pioneers in this field have laid fundamental groundwork regarding effective medium descriptions, there exists a discrepancy between discrete and continuum level solutions. In our previous work, we elaborated a Particle Mechanics Approach (PMA) that integrates thermal contact and Hertzian deformation models to understand the thermo-mechanically-coupled consolidation problem. We also considered the analogous problem from the perspective of the conventional Continuum Mechanics Approach (CMA). In this study, following the multi-scale modeling framework, we propose an effective thermal expansion coefficient for the thermally-assisted compaction of granular materials.


Introduction
In-depth understanding of the thermo-mechanical coupling that takes place during thermally-assisted compaction of granular materials is the key to successful design and optimization of many powder process engineering applications. It is well known that the elastic properties of bulk materials are altered by a change in temperature. Moreover, as also pointed out by Jaeger et al. [1], granular materials exhibit exceptional properties, which differ from those of the primary bulk material and depend on the processing conditions of the granular system. In this regard, there is a significant number of efforts to elucidate the collective behavior, which is also referred to as effective thermo-mechanical properties, of granular materials in response to boundary conditions and system parameters.
Mathematical models were developed to understand the thermal contact of conforming elastic surfaces [2][3][4]. Starting from the definition of thermal flux at the smooth contact surface of two spherical particles, Batchelor and O'Brien solved the analytical equation that expresses the effective thermal conductivity of ordered and randomly-distributed packed granular assemblies [4]. Experimental studies validating these models were presented in the related work of Hadley [5], Nozad et al. [6] and Shonnard and Whitaker [7]. More refined studies relaxed these assumptions by focusing on elasto-plastic contacts [8] or rough non-conforming surfaces [9][10][11]. Recently, the field of porous media has continued to attract researchers in light of understanding the correlation between geometry, loading conditions and anisotropy in the microstructure, which in turn affect the macroscopic behavior of compacted beds [12]. Yun and Santamarina conducted a fundamental study on thermal conduction through one-dimensional granular chains of spherical particles [13]. They pointed out the decisive role of inter-particle contacts in heat transfer mechanisms and highlighted the need for the effective thermal conductivity models to consider the inherent presence of contacts in particulate materials [13]. Exploring the effects of packing grains during thermal cycling, Chen and his co-workers [14] showed that thermal expansion, due to an imposed thermal gradient, has a significant effect on the rearrangement of the particle bed. Moreover, heat transfer in granular flows has been of great interest for numerous industrial applications, particularly in regards to understanding the effect of process conditions on heat transfer mechanisms in granular media [15][16][17][18][19].
In recognition of the uniqueness of granular materials and their importance to a wide variety of thermally-assisted compaction processes, two main methodologies have received considerable attention for addressing the need for simulating and predicting the macroscopic behavior of granular materials. One of the most adopted approaches treats the particles as individual bodies such that the coupled effects of various multi-physical phenomena are described at the particle-scale. Accounting for particle interactions and adopting constitutive relations of contact mechanics [20][21][22], the discrete element approach has been widely used in the field of particle-scale research [23]. This method is based on the early work of Cundall and Strack, where authors introduced an explicit numerical scheme to describe granular dynamics by tracing the motion of the particles and the generation of forces over the contact network while solving discrete equations of motion [24]. The interrelationship between particle motion and energy with the macroscopic behavior of the assembly provides understanding of the overall behavior of the confined material [25]. The main advantage of this methodology is the capability of presenting broad information about the microstructure of the granular material [26]. Although there exists a considerable amount of computational challenges in modeling a large number of particles with discrete element methods, improvements were done in formalisms, and new simulation techniques increased the achievability of calculations at the particle level.
The second most implemented methodology to model the collective behavior of particulate materials is the continuum mechanics approach, in which the granular material is assumed to be statistically homogeneous [27]. This simplification is achieved by treating the system as units of ordered arrays, simulating disordered arrangements by statistical correlation functions or using empirical correlations. The statistical averaging techniques provide homogenized solutions of the highly heterogeneous granular media at the cost of imposing two assumptions: (i) affine motion approximation, namely the motion of each grain follows the macroscopic strain, and (ii) well-bonded structure, contact number and positioning do not change under the applied load. Despite the fact that effective medium theory particularly estimates the effective elastic moduli of a packed bed of spherical particles to a large extent, the discrepancy between numerical and experimental results is remarkable. Makse and co-workers questioned the relevance of force laws defined at the single contact level, where they pointed out that the simplification done in effective medium theory is the misleading element in the formulation [28,29]. The affine motion assumption demolishes the ability of the approach to account for the relaxation and rearrangement of particles that are under shear deformation. Moreover, concerning the variety of boundary conditions and geometrical effects, experimentation techniques become insufficient in providing reliable information about the microstructure to feed empirical correlations.
The improvement of theoretical models and numerical simulation schemes remains an active area of research in the study of granular matter. A quasi-continuum approach has been recently developed, and the formulation was used for simulation of inter-particle bonding in granular systems [30]. Zheng and Cuitiño implemented the quasi-continuum approach to bridge the micro-and meso-scale through a discrete-continuum formulation of elastic-inelastic deformations occurring in the post-rearrangement regime of consolidation of inhomogeneous granular beds [30]. Since this approach provides the flexibility of storing individual particle interactions in an finite element modeling (FEM) scheme, it provides the overall behavior of the entire body without losing critical information specific to the microstructure. Koynov et al. presented a notable adaptation of this approach on the topic of powder compactions for pharmaceutical purposes [31]. Gonzalez and Cuitiño introduced a new formulation that accounts for the interplay of nonlocal mesoscopic deformations characteristic of confined granular systems. In the absence of the classical restriction of independent contacts of the Hertz law, the extended theory of nonlocal contact formulation provides predictive models at moderate levels of deformation and high confinement [32].
The present work attempts to establish a relationship that defines the effective thermal expansion property of a packed bed of spherical particles under the effect of the thermally-assisted compaction process. We elaborate on a Particle Mechanics Approach (PMA) that entails the integration of contact mechanics principles with a thermal-contact model to account for the heat conduction at the quasi-static equilibrium of the deformed state of granular materials. The disordered nature of the problem leads to highly non-linear coupled equations; therefore, we investigate a regular packing to simplify the problem and make it mathematically traceable. Moreover, we consider the analogous problem from the perspective of the conventional Continuum Mechanics Approach (CMA), while practicing the effective medium approach descriptions that define the macroscopic thermal and mechanical properties. Similar in spirit to the work of Chan and Tien [2], who proposed the effective thermal resistance, and to the work of Walton [33], who presented a method to calculate the effective elastic moduli of granular packing, we adopt a multi-scale approach to link the particle level information to the continuum level description of the thermally-assisted compaction process. Finally, we derive the equation of effective thermal expansion coefficient for the regularly-packed bed of particles.

Particle Mechanics Approach
Particle scale modeling of the thermally-assisted compaction process requires an extension of the discrete element method to account for the integration of heat conduction (e.g., [25,34,35]). Starting from the well-known theory of Hertzian deformation [20], heat conduction through the conforming contact of spherical particles [2,4] is adopted for the case of thermally-assisted compaction. Under steady state conditions, the total of the forces acting on individual particle m from neighboring particles n ∈ N m and the total heat transferred to particle m are zero, that is: where n mn is the unit normal vector defined from x n to x m , i.e., from the center of particle n to the center of particle m. Johnson studied the elastic deformation of locally spherical particles that are subject to a compression load by contact mechanics considerations [36]. Small-strain deformation of conforming surfaces results in a flat circular contact area. The collinear, elastic, contact force between particles m and n is defined through reduced elastic modulus E mn , reduced particle radius R mn and overlap γ mn between these particles. Specifically, the contact force is: where: Similar to previous studies in the literature [35,37], in the present study, a linear thermal expansion formulation is taken into consideration; that is , where α m is the thermal expansion coefficient, T re f is the reference temperature and R m re f is the radius of the particle at the reference temperature. Due to the fact that the contact geometry depends highly on the heat conduction between the consecutive conforming particle pairs, it is expected to capture a distribution of the contact area formation throughout the compacted medium.
The major heat transfer mechanisms in compacted particle beds consist of conduction through solid particles, conduction through the contact area between two touching particles, conduction to/from interstitial fluid, heat transfer via convection, radiation between particle surfaces and radiation between neighboring voids [34]. For a system of granular media where the thermal conductivity of the solid particles is much larger than that of the interstitial medium, the driving mechanisms for the heat transfer are the first two. Restricting attention to the problem of thermally-assisted compaction of spherical particles in a vacuum, we focus on thermal contact models that consider the conduction through solid particles and the contact areas between touching particles.
The analytical solution of the heat conduction through the solid phase of ordered spherical particles has been proposed by Chan and Tien [2] and Kaganer [3]. Moreover, the problem of heat transfer regarding the compaction of particles that are in or nearly in contact is deeply investigated by Batchelor and O'Brien [4]. In an attempt to find the approximate effective thermal conductivity of ordered and randomly-packed granular beds, Batchelor and O'Brien discussed the heat flux across the flat, circular contact surface between smooth, conforming and elastic particles. In this study, we adopt Batchelor and O'Brien's model for predicting the heat conductance, which is the ability of two touching surfaces to transmit heat through their contact interface. Heat flux across the contact area of two spherical, smooth particles is given by: where k mn is the arithmetic mean of the thermal conductivities of two conforming particles and a mn is the Hertzian contact area. These are defined as: The total heat flow to an individual particle, Equation (2), is calculated by adding the heat flow, Equation (7), across each contact surface shared with its neighboring particles. Thermal contact models introduced in the literature [2,4], Equation (2), assume that the resistance to heat transfer inside the particle is significantly smaller than the resistance between the particles, i.e., a Biot number much less than one: 2k mn a mn k mn A/R m 1 (10) where A is the cross-sectional area, A = π(R m ) 2 . This assumption was applied by several authors in earlier studies [34,38], which also enforces the condition of a mn R m , i.e., of small-strain deformation of elastic bodies in contact.
Referring to the previous experimental studies on regular and random packing of granular media, Walton points out that although the regular packing models are founded on strict assumptions, they are capable of capturing the vast majority of the characteristics of a real granular media [33]. In the present study, we consider a simple cubic packing of identical elastic spheres, which are constrained between parallel planes of infinite extent. A compression load and a temperature gradient are applied along the major and finite direction. Stress and heat flux are defined to depend only on externally-applied thermal and mechanical loads, and the weight of the particles is neglected. For such regular packings, each layer of the arrangement is isothermal normal to the direction of applied load. Furthermore, since these transversely-oriented particles are, at most, at the contact point, for each particle there is only one pair of contact areas aligned with the direction of applied thermal and mechanical load. Due to the symmetry of the problem, it is sufficient to consider a single column of a square cross-section containing the longitudinally-compressed spheres together. The above-described set of concepts regarding regular packings is also encountered in the early work of Chan and Tien [2] and Kaganer [3]. Based on these assumptions, the specified granular media can be visualized as a chain of elastic particles compressed between two walls, which are maintained at different temperatures, as seen in Figure 1. Details of the particle mechanics approach adopted in this study can also be found in detail in our earlier work [39].

Conventional Continuum Mechanics Approach
There has been considerable research directed towards describing the macroscopic behavior of compacted granular materials by using various homogenization techniques and postulating continuum constitutive laws [40]. Some of the previous studies on mathematical modeling of transport properties are aimed at estimating elastic-plastic mechanical properties, thermal and electrical conductivity of ordered and disordered arrangements. In addition to the particle-level approach, we also focus on a small-strain thermoelasticity model of continuum scale description that integrates the previously-proposed effective mechanical and thermal properties for granular beds under compaction. In this study, we refer to the particle mechanics approach and the conventional continuum mechanics approach as PMA and CMA, respectively.
The governing field equations of motion and energy of the analogous problem defined at the continuum scale are the following: div (σ) = 0 and div [k grad (T)] = 0, where Cauchy's stress, σ, is formulated as a combination of classical linear elasticity theory and simple linear thermal expansion, that is: where I is the identity matrix. The solution for the basic one-dimensional steady state thermoelastic, continuum problem, where body forces are neglected, depends linearly on elastic constants, λ, µ, thermal expansion and conduction coefficients, α and k, respectively. Since ε 22 = ε 33 = 0 holds, ε 11 is referred as ε(x), and it is defined positive for compression. The system of questions then reduces to: with σ(x) positive for compression, and q = k∂T/∂x.
Effective mechanical properties of granular beds are of great interest for numerous theoretical studies, some of which focuses on: (a) the principal elastic modulus for vertical compression of spherical particles without any lateral extension (Walton [33]); (b) finite and incremental elasticity of random packing of identical particles using energy methods (Norris and Johnson [41]); (c) enhancement of the derived formulas based on the pressure dependence of the elastic moduli of granular packings (Makse et al. [28,29]). The effective medium approach proposes the following elastic effective properties,λ andμ: where C n is the actual stiffness that depends on the bulk mechanical properties: Young's modulus, E, and Poisson's ratio, ν. φ s is the packing fraction, and Z is the coordination number. The effective thermal conductivity of a granular bed is substantially sensitive to the thermal and elastic properties of individual particles. In this study, we adopt Batchelor and O'Brien's [4] solution for effective thermal conductivity coefficient: It has been shown that the above-mentioned thermal contact models provide accurate results in estimating steady and average temperature profiles for ordered granular packings [42].
After implementing the effective mechanical and thermal properties in the resembling continuum description, the equation of stress becomes: where T w 1 and T w 2 are the temperature at the constraining walls, and ε is the compaction strain along the principle direction. The overall compaction force can simply be expressed as F = σ4R 2 re f , where [.] + = max{.,0} (notice that since σ(x) and ε(x) are assumed to be positive for compressive stress and strain, the above equation is valid for positive values of the expression in the parentheses).

Comparison of the Particle Mechanics Approach and the Conventional Continuum Mechanics Approach
According to the Hertz theory [43], the collinear contact force between the compressed elastic particles is a nonlinear function of the overlap, which is formed under the effect of the external load acting on the particles. For the case of thermally-assisted compaction of a granular system, such dependency is altered under the effect of an applied thermal gradient. To illustrate this influence, we work on the analytical solutions proposed in PMA and CMA for a system of stainless steel spherical particles. 304 stainless steel has the following bulk properties: E = 193 GPa, k = 15 W/mK, ν = 0.29 and α = 17.3 10 −6 1/K. The packing fraction is taken as φ s = π/6(1 − ε), and the coordination number as Z = 6.
In Figures 2 and 3, we aim to compare the analytical solutions generated in the particle mechanics approach and in the conventional continuum mechanics approach under varying boundary conditions for the relevant thermally-assisted compaction simulation. The compaction force and the heat transferred that are calculated by these two approaches are shown in the ratio. Figures 2 and 3 reflect a discrepancy between PMA and CMA particularly for the deformation range of high thermal load and low mechanical load.

Derivation of the Effective Thermal Expansion Coefficient for Thermally-Assisted Compaction of Granular Beds
Similar in spirit to earlier studies, the multi-scale approach discussed in this work is used to link the particle level information to the continuum level description of the thermally-assisted compaction process. We present a methodology to derive an effective thermal expansion coefficient for confined granular systems. The superscript mn is used to refer to the particle interactions at the contact of individual particles m and n. In this section, the quantities defined at the particle-level description and the continuum-scale description are denoted slightly different. For instance, stress and heat flux are defined as σ mn and q mn in PMA and σ(x mn ) and q(x mn ) in CMA.
Considering the case of an infinite chain of identical particles, we assume that: (i) the temperature difference between consecutive pairs is negligible compared to the change with respect to the reference temperature (i.e., (T m − T n )/(T mn − T re f ) 1 ); and (ii) the particles are locally subject to uniform compaction. Moreover, we consider an average porosity for the resembling systems of particle-scale and continuum scale analysis; therefore, φ s for a given compaction strain ε is calculated as π/6(1 − ε). Based on the above assumptions, the average stress calculated in PMA can be expressed as: whereas the continuum-level approach predicts the particular compression stress as (valid for compressive stresses σ(x mn ) > 0): whereα is the effective thermal expansion coefficient. Next, we enforce stress and heat flux expressions in both descriptions to be equal, that is σ(x mn ) = σ mn and q(x mn ) = q mn . x mn , T(x mn ) and ε(x mn ) are: Finally, the equivalence of Equations (17) and (18) provides the following continuum level effective thermal expansion expression that is dependent on the applied mechanical and thermal load, σ and T − T re f , as well as the bulk thermal expansion property of the solid.
The first order approximation of the above expression is: To quantify the overall effect of this approximation, we implement both equations, Equations (22) and (23), in the continuum-scale solution, and we compare the required stress and heat flux particularly for the previously discussed deformation range. It is found that the two compared CMA solutions differ less than 1% under these conditions.

An Application of the Proposed Effective Thermal Expansion Coefficient
Along with the effective mechanical properties and effective thermal conductivity (listed in Section 2.2), we implement the first order approximation of the proposed thermal expansion coefficient, Equation (23), to solve the described thermally-assisted compaction problem. The analytical solution is given in the following set of equations: where L i is the initial length of the system, and σ is given by: It is also worth noting that in the limiting case, when there is no thermal load, the derived solution for stress, Equation (26), approaches the solution given in conventional continuum mechanics approach, Equation (16).

Results and Discussion
In order to evaluate the effect of the proposed effective thermal expansion coefficient on the continuum-scale approach, we compare the three analytical solutions discussed in this study: (i) PMA; (ii) conventional CMA, where effective mechanical properties and effective thermal conductance are implemented; (iii) an improved continuum mechanics solution, where also the proposed thermal expansion coefficient is included, the Improved-CMA (also given in Equations (24)-(26)). We focus on four critical boundary conditions in detail; high/low thermal load and high/low mechanical deformation ranges.
In Figure 4a,b, the force needed to compress the system up to a compaction strain, ε, of 5% is traced under two different thermal load conditions. When the temperature difference between the two boundary walls is only 40 K, Improved-CMA overlaps with the conventional continuum solution, as expected; whereas under high thermal load such as T w 2 − T w 1 = 1000 K, and particularly at low compaction strain, it is seen in Figure 4b that Improved-CMA has a better estimation of compaction force in terms of converging the PMA solution.
Moreover, we evaluate the effect of thermal load on the compaction force (Figure 5a,b). Figure 5a shows that the proposed expression for the effective thermal expansion coefficient significantly improves the continuum-scale solution in predicting the compaction force. Under high mechanical load, Figure 5b suggests that Improved-CMA shows a trend in compaction force similar to the one predicted by PMA. However, there exists a difference between continuum-scale solutions and the particle-level solution at zero thermal load. The authors claim that this discrepancy stems from the calculation of the effective mechanical properties, which is also discussed in detail in the earlier study of Makse et al. [28].  Figures 6a,b and 7a show a good agreement in terms of the heat transferred between particle-scale and continuum-scale models. Even though Figure 7b indicates that there is a difference between these two modeling approaches under varying compaction strains, it is worth noting that the maximum difference between Improved-CMA and PMA is 6%.
In our previous study, we also focused on the position of individual particles under the thermally-assisted compaction process [39]. For a one-dimensional chain of particles, we noted that the relative difference in estimating the nodal position compared to the particle position between the conventional continuum mechanics approach and particle mechanics approach is up to 40% under high thermal and low mechanical load conditions. Therefore, in the current study, we plot the displacement versus the initial position of each node/particle in Figure 8a,b. It is assumed that the node in contact with the non-moving boundary is placed at x = 0. We conclude that the implementation of the proposed effective thermal expansion coefficient significantly improves the conventional continuum solution in terms of predicting the displacement of the individual particle. (b) Figure 8. Comparison of the displacement profile of each node/particle; (a) T w 2 − T w 1 = 1000 K and ε = 0.005; (b) T w 2 − T w 1 = 1000 K and ε = 0.05.

Conclusions
The present work centers on a multi-scale approach to bridge the gap between discrete and finite-scale solutions of a thermo-mechanically-coupled problem while introducing an effective thermal expansion coefficient. The response of a granular system under thermally-assisted compaction shows a high dependence on the thermal expansion of the particles. A discrete-system solution based on a particle mechanics approach carries out this relationship, and it successfully shows nonlinear effects on thermal strains due to thermal expansion. Despite the fact that effective medium theory enhances the conventional continuum mechanics model to a large extent, there still exists a notable discrepancy between these two approaches. In this study, we address this gap by incorporating a previously-suggested methodology to identify the effective thermal expansion property of granular materials. It is shown that the implementation of the proposed effective thermal expansion coefficient significantly improves the conventional continuum mechanics analysis, thereby resulting in better accuracy in predicting particle-level characteristics of the thermally-assisted compaction problem. The extension of the proposed approach to other multi-physics phenomena appearing in granular systems in multi-dimensional analysis is a worthwhile direction for future research.