Analytical and DEM studies of thermal stress in granular materials

The ability to predict thermal-induced stresses in granular materials is of practical importance across a range of disciplines ranging from process engineering to geotechnical engineering. This study presents an analytical formula to predict thermal-induced stress increments in mono-disperse granular materials subject to an initial isotropic stress state. A complementary series of DEM simulations were carried out to explore the applicability of the proposed analytical formula. The comparative analysis showed that the proposed expression can accurately predict stress changes in packings where there are negligible particle displacements as a consequence of the thermal loading (e.g. regular packings and medium/dense random packings); however large errors were observed in loose samples with a random packing.


Background
Granular materials are encountered in various industrial and engineering fields. Granular materials are often subjected to heating / temperature change due to various internal and external factors such as chemical/nuclear reactions, proximity to geothermal structures such as energy piles, use of packed-bed thermal energy storage, heat application in industrial processes, and daily/seasonal temperature variations. Materials, in general, experience thermal expansion when their temperature increases. Under constrained conditions, this thermal expansion can induce significant increments in stress, which may have serious consequences (e.g. buckling of rails during hot summer periods). The application of thermal cycles to granular media are also of interest since cyclic changes in stresses induced by such loadings lead to volume change as demonstrated in previous studies [1][2][3][4]. Nonetheless, there has been only limited consideration of thermal stresses in granular materials while thermal stresses in continuum materials have been well discussed in the scientific literature.
The behaviour of granular materials is difficult to predict as they are formed of discrete particles that move relative to each other. The Discrete Element Method [5], which simulates movement of granular particles based on Newton's law of motion and contact mechanics, has emerged as a valuable tool to investigate the behaviour of granular materials. One drawback of DEM is its high computational cost. Furthermore, any user of DEM must also have access to DEM software and the skills required to run high-quality simulations. An alternative is to use Effective Medium Theory (EMT). Using EMT, [6] have introduced statistical approximations to estimate bulk properties of granular materials * Corresponding author: t.morimoto@imperial.ac.uk analytically. While EMT has known limitations, an analytical formula which estimates thermal stresses in granular materials is useful to enable overall understanding and to verify numerical models, including DEM models.
Here, thermally-induced stresses in granular materials were analytically estimated developing a new expression based on EMT. DEM simulations considering different packings of monodisperse particles were carried out. By comparing the responses predicted by the analytical and numerical models, the advantages and limitations of the analytical approach are critically assessed.

Analytical estimate of stress
Rothenburg and Bathurst [7] amongst others have predicted the isotropic stress ( ) in a mono-disperse, isotropic granular material from the contact density ( ), the mean distance between the centres of contacting particles ( ̅ ), and the mean normal contact force ( ̅ ) as If the granular material consists of mono-disperse particles with radius, , then ̅ ≅ 2 (2) where is the number of particles considered, and ̅ , and are the mean coordination number, the volume and the void ratio of the sample, respectively. The void ratio is related to the solid fraction ( ) as = (1 − ) ⁄ . Using equations (2) and (3), equation (1) can be rewritten as:

Increment in contact force
According to the simplified Hertzian contact model, at temperature the normal contact force ( ) between two identical smooth particles with radius , shear modulus, , and Poisson's ratio , is: is the distance between the centres of the two particles and is the contact deformation (overlap). If the temperature of the two particles increases to + Δ while the centres of the two particles remain fixed in space, and increase due to the thermal expansion of the particles: is the particles' linear thermal expansion coefficient. Consequently, the contact force between particles increases to:

Analytical estimation of stress increase
To derive an analytical estimation for the isotropic stress induced in a granular medium by heating, two assumptions are made: (1) the change in ̅ due to a temperature change from to + Δ is negligible, and (2) the increment in ̅ due to the temperature change follows equation (9). Then, the isotropic stress after the increment in temperature is estimated to be: where ̅̅̅ is the mean value of in the material, which can be estimated from equations (4) and (5) to be: A material variable, , and a packing variable, , are defined as: Thus an analytical equation which expresses +Δ , as a function of the material properties ( and ) and initial conditions ( and ), is obtained:

DEM simulation approach
DEM simulations were conducted to examine the validity of equation (14). A modified version of the open-source code LAMMPS [8] was used. Two regular packings, a Face Centred Cubic (FCC) packing and a Simple Cubic (SC) packing with different values, as well as various random packings with different and values were considered in a parametric study. The simulation parameters are summarised in Table 1; the material properties of particles are those typical of glass beads that are often used in fundamental laboratory studies.

Sample generation
Each cubic sample used was enclosed by periodic boundaries. For the two regular packings, the particle positions were analytically generated so that the particles were initially very close, but not overlapping.
To obtain isotropic and homogeneous samples in the case of the random packings, the particle positions were incrementally randomly generated in a cubic domain whose side length is 20 mm. Particle overlap was not allowed and particles were added one by one until a void ratio of 2.0 was reached. There were 5093 spheres in the simulations of the random packings.

Isotropic compression
Each packing was compressed until the isotropic stress ( ) reached 1 kPa with a maximum strain rate of 1%/s. Viscous damping with a damping coefficient ( ) of 0.10 was activated to reduce the kinetic energy of the particles while maintaining = 1 kPa. Then, was gradually increased with a maximum strain rate of 0.2%/s. Three values, i.e. 10kPa, 100kPa, and 1000 kPa were considered.
For the random packings, three inter-particle friction coefficient ( ) values were applied during isotropic compression to achieve three samples with different values. The packings were achieved using = 0.01, 0.10, 0.34 and denoted 'Dense', 'Medium', and 'Loose', respectively. The value was set to be 0.0 during compression of the FCC and SC packings.

Radius expansion
After the isotropic compression procedure, the positions of the boundaries were fixed, and the particle radii were increased to simulate thermal expansion of the particles following [3]. Three temperature increments (Δ ) of 1, 10, and 50 K were considered and for each Δ value, the expanded particle radii were calculated using value of 9.0 × 10 −6 , which is a typical value for glass beads.
In the randomly packed samples two scenarios were considered after the radii were expanded. In the first case (denoted 'Fix'), the positions of the particles were fixed, and the isotropic stress was measured immediately after the radius increase procedure. In the second case (denoted 'Release'), particle movements were allowed with = 0.35 and viscous damping ( = 0.10), and the isotropic stress was measured after the stress value converged to a stable value. Only the Fix case was considered for the regular packings, as all the particles will be in a state of equilibrium even after the radius expansion procedure.

Discussion
As discussed in Section 2, several assumptions were introduced to obtain equation (14): i. The granular material is mono-disperse. ii.
The granular material is isotropic. iii.
The material properties of the particles remain constant under a change in temperature. iv.
The overlaps between the particles ( ̅ ) are neglected when calculating ̅ in equation (2). v.
The change in ̅ due to a change in temperature is negligible. vi.
The change in ̅ due to a temperature change follows equation (9). vii.
The positions of the particles are fixed.
Assumptions (i), (ii), and (iii) hold only for limited cases and in principle this approach can be used to derive a more general expression. However, it is worth exploring the validity of the assumptions (iv) to (vii) before implementing this challenging generalisation. This section compares the analytical predictions with the DEM data for the samples considered here so that the errors induced by these assumptions can be examined.

Regular packings
In a regular packing under an isotropic stress condition, where ̅ remains constant, the contact forces are always completely homogeneous, the particles are always in a state of equilibrium, and assumptions (v) to (vii) listed above apply. Only assumption (iv) might lead to a difference between analytical and DEM +Δ values.
The and values used in equations (13) and (14) were determined analytically. Referring to Tables 2 and 3, nearly perfect agreement was obtained between the analytical and DEM +Δ values, confirming the validity of equation (14) for regular packings. Errors due to assumption (iv) are ≤ 0.1% for the conditions considered in this study. These errors will increase when 'soft' granular materials with a low are considered.

Random packings
This section considers the validity of equation (14) for random packings, where, in contrast to the regular packings, assumptions (v) to (vii) do not hold true. As outlined in Section 2, two simulation cases, namely Fix and Release, were considered. The and ̅̅̅ values of the random packings considered in this study are summarised in Table 4. Analytical +Δ values for the random packings were calculated using these values and compared with DEM +Δ values (for both the Fix and Release cases) in Fig. 1 for the middle stress level considered (100 kPa). Both the analytical model and DEM data provide the same useful insight that for a given Δ , Δ Δ increases with reducing .

Fix case
In the Fix case, where particle movement is prohibited, assumption (vii) is enforced in the DEM model. The differences between the analytical-and DEM-generated data are given in Table 5. A positive value indicates an analytical result that is higher than the corresponding DEM result. In the case of a small temperature increment (e.g. Δ = 1 ), where the increase in ̅ is small, the differences are < 6.3%, and they can be attributed to ignoring the inhomogeneity in the contact forces (assumption (vi)). As Δ increased, the differences decreased because the errors due to assumption (vi) are counteracted by ignoring the increase in ̅ (assumption (v)) which leads to underestimation of +Δ .

Release case
In general, if the radii of particles in a random packing are increased, the particles do not remain in static equilibrium. This is true even if the changes in radii take place simultaneously, due to the non-linearity of the contact stiffness and the generation of new contacts. These particle movements lead to a redistribution of the contact forces and reduction in the overall stress. The differences between the predictions of Δ Δ for analytical and DEM data are shown in Table 6. In all cases considered, the analytical approach overestimated +Δ because the stress reduction was neglected (assumption (vii)). The differences increased with increasing Δ and reducing sample density; both of these factors increase the likelihood of particle movements occurring. While the errors for dense and medium packings were within 20%, equation (14) is clearly not reliable for the loosest sample considered.

Conclusions
This study has presented an analytical expression to predict thermally-induced stress increments in isotropic mono-disperse granular materials. Data generated using the new expression were compared with complementary data generated from DEM simulations. The analytical expression achieved very good predictions for regular packings and good predictions for random packings where particle movements are fixed. However, the analytical expression was shown to be ineffective in the case of loose random packings where particle movements are allowed. This result suggests that ignoring the redistribution of contact forces and the associated reduction in stress in the analytical formulation induced a more critical over-simplification when compared with ignoring the inhomogeneity in the contact forces or neglecting the increase in coordination number due to the temperature increment. The proposed analytical formula is useful for verifying DEM simulations using regular packings; however, the applicability for random packings is limited without further modifications.