A simple discrete-element model for numerical studying the dynamic thermal response of granular materials

This paper aims to investigate the influence of periodicity temperature change on the properties of dry granular materials in macroscopic and microscopic. A series of cyclic thermal consolidation tests have been carried out based on the discrete element method (DEM) that incorporate particles’ volumetric thermal expansion coefficient. The simulation of the direct shear test was carried out on the samples after thermal cycling. Results showed that thermally-induced volumetric strain accumulation of the specimen can be calculated by the DEM model, based on the two-dimensional particle flow code (PFC2D) software. The lateral pressure degraded concomitantly thanks to decreases in particles’ horizontal contact during periodic thermal cycling. In addition, the shear dilatancy level decreases during the shearing process with the number of thermal cycles. Both the size and anisotropy of the normal contact force and contact number and the force chain are affected by the temperature cycle. Finally, the results of this paper have a certain reference for the engineering practice, such as thermal piles or others, when granular materials are subjected to thermal cycling.


Introduction
Previous investigations have been shown that even a small perturbation could have a significant effect on the dry granular materials. For instance, vibrations can densify the packing [1], and vibrations frequency has a significant effect on packing density [2]; stiffness, strength, thermal conductivity, and permeability of granular materials varies with compression and shear [3][4][5]; periodic compression can promote the consolidation speed of granular materials [6]. The above studies were based on the external perturbation on the boundaries. Moreover, the internal perturbation, e.g., temperature change and fluidization, also has a significant effect on granular materials.
Temperature change is one of the most common internal perturbations. For example, soil near the energy pile [7,8], high-temperature submarine pipeline [9], and nuclear waste treatment plant [10] are subjected to a periodic temperature change. A large quantity on the influence of temperature change has been reported. Temperature change has a significant effect on friction resistance [11] and contact pressure [12] of granular materials, inducing displacement [13]. and irreversible settling [14] of energy structure. In addition, thermally induced mechanical responses of the granular materials element have also been investigated based on laboratory tests while taking into account accumulation of stress-strain [15], heat transfer behaviour [16], thermal resistance [17] and thermal conductivity [18]. Based on the stress-strain behaviour, phenomenological constitutive models have also been used to describe related shear dilative behaviour [19][20][21][22][23][24] and are commonly adopted in engineering practice for their efficiency in finite element analysis [25][26][27][28]. However, such studies focused mainly on characterising phenomenological macroscopic-scale behaviour, so that the effects of long-term periodic accumulation thermal loading on the force chain and the fabric anisotropy of the granular materials element in the macro-to-micro scale remain an open question.
Thermally related microscopic-scale studies increasingly rely on laboratory tests to highlight the effects of different temperature treatments on the processes that damage bonded granular materials [29][30][31][32][33]. A microscale computational technique that uses the discrete element method (DEM) [34] can also be used to observe the basic macro-to-micro behaviour of granular materials [35][36][37][38], and indeed many thermally related microinvestigations have been carried out based on this DEM method, whether to assess heat conduction's effects on the thermal expansion [39], incorporate the fabric tensor into thermal conductivity [40], estimate stress-induced thermal conductivity tensors [4], highlight the contributions of heating rate and thermal cycle magnitude to force relaxation [41]. Yun and Evans [42] evaluated the thermal conductivity of granular materials through a three-dimensional discrete element model. Few studies, however, have sought to explain the cyclic thermally induced behaviour of volumetric strain accumulation and shear dilatancy by considering the evolution fabric anisotropy of granular materials in the microscale. This paper analyses the thermally induced stress-strain response of granular materials through DEM direct shear tests. Different from the thermally coupling model [42], a simple approach that incorporates particles' volumetric thermal expansion coefficient, without considering the thermal transfer was developed [41]. A series of cyclical thermal consolidation tests and direct shear tests are then performed using the commercial DEM code [43]. Furthermore, cyclic thermally induced densification and strength degradation are interpreted according to the evolution of contact anisotropy during the corresponding loading process. The research results can serve as a reference for numerical study for energy piles engineering.

A simple thermal DEM model
The DEM model can be calculated through the alternate application of the force-displacement law and Newton's first law of motion. Particles' contact force consists of normal contact force and tangential contact force, where F i is the contact force of the particle i, F i n , the normal contact force, and F i s , the tangential contact force, which can be expressed as d where K n is the normal stiffness of the contact point, δ i,n the amount of overlap between particles in contact, F i,s ( * ) is the tangential contact surface at the beginning time step, K s the tangential stiffness of the contact point. The greater tangential stiffness, the smaller the shear strain of the particles under a constant loading level. Δδ i,s is the tangential component of the contact displacement increment, obtained by Δδ i , s =V i,s ·Δt, in which V i,s is the tangential component of the contact point velocity and Δt is the time increment of calculation, d is the centroid distance between particles.
K n and K s are calculated by equation (4).
where A is the effective contact area, E * is the effective modulus, L is the distance between the particles or the distance between the particles and the walls, k * is the stiffness ratio, r, R 1 and R 2 are the radius of the particles, as shown in figure 1.
Particles' motion is determined by the Newton-Euler equations, where m i and I i are the mass and moment of inertia of particle i, respectively; v and ω are the speeds of displacement and rotation; F b is the total volume resultant force; N is the number of particles in contact with the particle; R i is the vector of equivalent radius that value depends on the soil particles following the grain size distribution.
As figure 1 shows, based on commercial software PFC 2D 5.0 sphere-based DEM codes, the effect of thermal expansion was incorporated as a temperature-dependent particle radius R, with other parameters (e.g., stiffness and friction coefficient calculated by (4) and (10)) remaining constant.
where, μ′ is the friction coefficient; f μ' is the angle of internal friction. Spherical sand particle expands uniformly follows the equation (11), which has been validated by Zhao et al [41].
where ΔR is the increment in particle radius at the reference temperature increment ΔT of particle and α is the coefficient of linear thermal expansion associated with the physical property of particles, R 0 is the radius of the particle before temperature changed.

Simulation program
To validate the shearing behaviour of granular material, the direct shear test was simulated, based on the particle volume expansion method of DEM. The shear box with 90 mm×36 mm was first generated in the PFC 2D software [44]. Then, noncohesive spherical particles with a diameter of 0.7 mm are randomly generated in the shear box under the gravity of 9.81 N kg −1 . It is aimed to build a simulation model with ideal conditions considering uniform diameter for fine silica granules, as shown in the black portion of figure 2(a). The initial specimen reached mechanical equilibrium (α f < = 1e −5 ) by the cyclical iteration of equations (1)- (9). At this time, the particles are considered to be in a relatively stable state, as shown in figures 2(a), (b), with the force chains (red lines) uniformly distributed between the particles and walls. All simulations were performed using the commercial software PFC 5.0. After the particles reached initial equilibrium, the lateral walls remained fixed.
A servo-controlled vertical load σ y of 1 MPa referenced by Duran-Olivencia et al [45] was applied to the top and bottom walls until the assembly reached a relatively stable state (average imbalance coefficient less than 1 e -5 ).
å å å a = where a f is the average imbalance coefficient, N is the total number of particles, C the total number of contacts, C i the number of contacts of the ith particle, F (ij) the contact force of the j th contact of the ith particle, F b (i) the volume force of the i th particle and F c (i) the contact force of the ith contact. The cyclic thermal consolidation tests were carried out by implementing the simple thermal approach equation (11) as part of the DEM calculating cycle. During which particles' temperature is cycled from an initial temperature of T 0 =20°C, with an amplitude of cyclic temperature T cyc =40°C [46] (the temperature cycle range is T 0 − T cyc ∼ T 0 +T cyc ), and the increment of temperature is fixed at ΔT=10°C (see figure 3). In addition, after each temperature change, the average imbalance coefficient α f must be calculated to be less than 1e −5 to ensure the stability of the model. During temperature cycling, the boundary conditions in which the rigid wall in the horizontal direction remains fixed, and the rigid wall in the vertical direction remains to keep a constant vertical loading. The thermally induced particle volumetric expansion can be calculated according to the equation (11).
In total, 4000 temperature cycles were carried out on the consolidated specimen, with constant vertical compression of 1 MPa, to study the strain accumulation behavior of the specimen under long-term cyclic thermal loading. Figure 2 presents the boundary conditions of the direct shear test, in which the servo stress in the vertical direction s y remains 1 MPa during the shear process and the lower shear box applies a constant shear rate of 4.23 mm s −1 in a horizontal direction, and the upper shear box remains fixed. At the same time, no external force is applied to the sides of the shear box. During shearing, the top wall is regularly monitored to calculate the stress-strain of the specimen, with shear stress τ and shear strain γ given by: where N l +N r is the resultant external force acting on the lower shear box, A is the contact area of the shearing interface that is continuously decreasing, ΔA is the shear area, and H 0 the initial height of the specimen. It can be seen from the equation (11) that the linear thermal expansion coefficient greatly affected the volume change of the specimen. Figure 4(a) shows the changes in the void ratio of the specimen in the cyclic thermal consolidation test when the value of coefficient a with 2e −4°C−1 , 1.5e −4°C−1 , 1e −4°C−1 ,1e −5°C−1 , and 1e −6°C−1 , respectively, in which a = 1 e −4°C−1 corresponds to a relatively stable consolidation curve for numerical study than other value of linear thermal expansion coefficient a. Figure 3(b) also shows the simulation results with different initial void ratio e 0 =0.2091, 0.2037, 0.1996 and 0.1943 respectively. It can be found the initial void ratio e 0 also affects the volumetric strain accumulation. Hereby, we take the coefficient of linear thermal expansion a = 1 e −4°C−1 and the initial void ratio e 0 =0.2091 as input parameters, to make the results of volumetric strain accumulation more obvious during numerical simulation. In addition, the DEM calculating paraments selected from references [47][48][49][50], the values of paraments are close to the silica, as shown in table 1.
To validate this numerical method, the simulation of a real thermal coupling test considering heat transfer and specific heat was carried out, with the following equations (based on the User's Guide of PFC [51]): where, Q (p) is the heat energy transferred from the particles to the heat pipe; N p is the number of heat pipes of a single particle; Q v is the intensity of the heat source; m is the mass of the particle; C v is the specific heat capacity; ΔT and η are the temperature difference and thermal resistance between the two ends of the heat pipe per unit length; L (p) is the length of the heat pipe. k and n are the thermal conductivity and porosity of the granular material respectively; V (b) is the volume of the particle; N b is the number of the particle. Thermal cycling consolidation tests with the number of cycles N cyc = 10 were simulated by these two methods using paraments in table 1, the result as shown in figure 5. The numerical method only considering volume expansion can well reproduce the basic stress-strain behavior of simulation coupling the real thermal transfer with paraments thermal conductivity k = 7.7 W/(m°C), and specific heat capacity C v = 880 J (kg°C) −1 .

Simulation results
The simulation results relating to periodic heating and cooling's effects on specimens' void ratio are plotted in figure 6. For each temperature cycle, particle volume continuously swells in the heating phase, decreasing the void ratio of the specimen. Conversely, in the cooling phase, particle volume gradually shrinks, causing the void ratio to increase, as seen in figure 6(a). In which the grey lines result from void ratios for each thermal cycle, with thermal cycles numbering 1, 1000, 2000 and 4000 marked in red to distinguish different stages of thermal cycling. The value of the void ratio at 20°C for each thermal cycle is recorded and plotted against the corresponding numbers of thermal cycles, as shown in figure 6(b). It can be concluded that during the temperature cycle, due to the irreversible rearrangement of the particle bed, the volumetric strain accumulation of the specimen is continuously increased. During temperature cycles, servo stress was applied only in the vertical direction, with the displacement of lateral walls kept constant, so that lateral stress σ x gradually decreased with periodic temperature changes, as shown in figure 7(a). Changes in lateral confining pressure at 1, 1000, 2000 and 4000 thermal cycles are marked accordingly in figure 7, in which for each temperature cycle the thermally induced expansion of particle volume increases confining pressure and the shrinkage of particle volume decreases it. Due to the densification of the specimen during the cyclic thermal consolidation testing, specimens' internal particles were rearranged, which also led to a cumulatively decrease in lateral confining pressure σ x to maintain energy balance.
To study the interface shearing behaviour of the specimen after thermal consolidation, a series of direct shear tests were simulated on the specimens with 1000, 2000, and 4000 temperature cycles. The investigation of the interface properties acting on a macro level was carried out based on the degradations of shearing strength and  dilatancy of specimens after thermal consolidation, as shown in figure 8. In addition, parallel tests that ignored cyclic thermal effects without considering volumetric expansion were also carried out, in which specimens' void ratio was the same as for tests considering thermal cycling. The simulation results are consistent with the laboratory results [52][53][54][55][56]. Figure 8 shows that the shear stiffness and strength of specimens affected by thermal consolidation were lower than those of specimens without temperature cycling, especially for the specimen that underwent the most thermal cycles (N cyc =4000), corresponding to a relatively dense state (e 0 =0.19): Its strength degraded most obviously during the shearing process.
Similarly, thermal cycling significantly affects volumetric strain during the shearing process. As shown in figures 8(d)-(f), the shearing dilatancy of the volumetric strain e v of the specimen(seen black curves), calculated by equation (17), after temperature cycling weakened significantly compared with the shear test that did not consider the temperature effect (seen red curves). A relatively loose specimen (e 0 =0.2035) corresponding to 1000 cycles was essentially in a continuous shearing contraction state, markedly different from the test results for the specimens without consideration of temperature cycles.
where V is the volume of the specimen in the shearing process; V 0 is the initial volume of the specimen.

Anisotropy analysis
Microstructure and micromechanics are known to be the underlying physics controlling granular materials' behaviours and responses. The anisotropy between granular materials is an important index for studying materials' microscopic mechanical properties. Thermal consolidation and direct shear tests can be used to investigate changes in the anisotropy field to explain the degradation of strength and dilatancy by noting contact numbers, the directions of the contact forces, and distributions of force chains between particles.
Orientations of contacts for the granular assemblies can be characterized by the density function E(q), where E(θ) is the proportion of an angular interval q and it also can be approximately expanded by the Fourier series [57], where a is an anisotropic parameter corresponding to the magnitude of anisotropy in contact orientations and θ a represents the main orientation of normal contact. When a=0, the distribution is isotropic. The level of anisotropy for the specimen's particles can be determined by calculating the anisotropic parameter a: First, the Fourier series of E(θ) is integrated: The degree of the anisotropy inside the container can be calculated using the preceding formula.

Contact number
To study the anisotropy of the contact number, in the cyclic thermal consolidation test the distribution diagram of the contact numbers with 0, 1000, 2000 and 4000 cycles were plotted, as shown in figure 9. The number of contacts was found to be increased, with the distribution in each direction becoming inhomogeneous as a result of the cyclic thermally induced consolidation. In this paper, the relationship of the total number of contacts between particles has been analyzed by equation (22), as shown in figure 10(a). The total contact numbers between particles increased from an initial 11181 to 11723, whereas the corresponding a increased from 0.0259 to 0.0743, per equation (22). It was found that in the process of increasing temperature cycles, as the number of particles in contact increased, anisotropy was increased.
Moreover, in the thermal consolidation test, the numbers of contacts in each particle corresponding to the 0, 1000, 2000, 3000 and 4000 thermal cycles were also counted, as shown in figure 10(b). As can be seen, the number of particles with 4 contacts make up the majority. It was found that as the number of temperature cycles increased, the numbers of particles with 5 and 6 contact numbers were increased, with 3 and 4 contact numbers were decreased.

Normal contact force
To study the microscopic mechanism of strength degradation during the cyclic thermal consolidation test, the behaviour of normal contact force between particles was analysed. The normal contact distribution shape equation is similar to the contact direction equation (20), which can be expressed by where f 0 represents the average normal contact force, θ n the main direction of the normal contact force and a n the anisotropy coefficient of the normal contact force. Based on equation (24), the distributions of the normal contact forces corresponding to various temperature cycles (0, 1000, 2000, 4000) are shown in figure 11, with the distributions of the normal contact forces remaining relatively stable in the vertical direction (θ =90°) due to the constant servo stress applied in this direction but decreasing continuously in the horizontal direction (θ =0°). This behaviour can be calculated by equation (24), with the red curve gradually changing from approximately circular (initial isotropy in all directions) to elliptical (shearing anisotropy in all directions).
The trends of changes to f 0 , a n and θ n in the cyclical thermal consolidation simulation are calculated by combining equations (22)- (24), as shown in figure 12. The average normal contact force f 0 of the particle is reduced from an initial 455 kN to 397.8 kN. The anisotropy coefficient a n increases from 0.091 to 0.262 and the  principal direction θ n of the normal contact force stabilizes at about 1.623 (93°). It can be seen that the temperature cycle has a greater influence on the magnitude and distribution of the normal contact force but a smaller influence on the direction of the normal contact force. Figure 13 is the normal contact force distribution curve fitted by the equation (24). When γ =0%, the temperature cycle has the most obvious influence on the distribution of normal contact force in the horizontal direction. In which, the normal contact force distribution in the horizontal direction when the temperature  cycle number is 4000 is concave, which is smaller than the distribution when the temperature cycle number is 0. However, in the direct shear simulation, the difference in the normal contact force distribution caused by the temperature cycle constantly decreases as the shear strain increases. At the same time, during the direct shear simulation, the direction of the maximum normal contact force distribution is gradually deflected from an initial ∼90°to ∼45°and remains stable after the shear strain γ >5%. Accordingly, the direct shear simulation will reduce the influence of temperature cycling on the distribution of particles' normal contact force.
The changing trend of f 0 , a n, and θ n in the direct shear simulation is obtained by the same method as in figure 13, as shown in figure 14. It can be seen from figure 14(a) that owing to an increase in shear stress (see figure 8), particles' average normal contact force increases rapidly when γ <7% and reaches a relatively stable state after γ >7%. At the same time, after γ >7%, the stable value of f 0 for different temperature cycles is very close. It can be seen from figure 13 that in the direct shear simulation, the normal contact distribution of the specimen under different temperature cycles is gradually unified, so that the anisotropy parameters a n is also stable at similar values, as shown in figure 14(b). In figure 14(c), the principal stress direction θ n of the particles' normal contact force has decreased from an initial 1.56 (about π/2) to 0.7 (about 4π/9), tσ in figure 13. In summary, during the direct shear test, temperature cycling has no significant effect on the change of the specimen's normal contact force.

Force chain distributions
To clarify the transmission path of the contact force between the particles during the temperature cycle, the force chain diagram of the specimen is drawn when the thermal cycle is 0, 1000, 2000, and 4000, as shown in figure 15. Before the temperature cycle, the distribution of force chains is homogeneous. Then, with the increase of the number of thermal cycles, the main force chain is assembled along the vertical direction, while the force chain in the horizontal direction gradually weakens.
To understand the behaviour of contact force distribution in the thermal consolidation test, the distribution curve of the contact force under different cycles (0, 1000, 2000, 4000) was fitted according to the equation (25), as shown in figure 16. It can be seen that a contact force of about 0.35 kN accounts for the largest proportion and that this proportion is increased from 0.78% to 0.95% as the number of temperature cycles increases. The average contact force μ calculated by the equation (26) is decreased from 0.803 kN to 0.693 kN, as shown in figure 16. The standard deviation σ calculated by equation (27) corresponds to the degree of stress inhomogeneity, which is decreased from 566.04 to 520.76, where F i is the contact force of the ith contact; C is the contact number; K(x) is the kernel smoothing function, using a Gaussian kernel equation (28); h is the bandwidth.  where μ is the average contact force and σ is the standard deviation. Figure 17 plots the force chain diagram under different shear strains (0.01, 0.025, 0.05, 0.15) . As the shear strain increases, the force chains are gradually assembled from the lower left and upper right. When the shear strain reaches 0.05, a shear band composed of obvious force chains forms from the lower left to the upper right.
To understand the influence of temperature cycling on the specimen's force chain distribution during the direct shear test when γ =0.05, force chain distributions are compared for specimens numbering 0, 1000, 2000, and 4000 temperature cycles, as shown in figure 18. It can be seen that as the number of temperature cycles increases, the shear band formed by the force chain weakens.

Conclusions
Based on the particle volume expansion method, the cyclical thermal consolidation test and direct shear test were simulated using two-dimensional particle flow code (PFC2D) software. The effects of thermal cycling on specimens were analyzed from macro to micro scale, and certain key conclusions were drawn: (1) A simple DEM simulation platform was established, incorporating particles' volumetric thermal expansion coefficient, that can well simulate specimen void ratio and stress degradation behaviour during the cyclic thermal consolidation test for granular materials.
(2) In the cyclic thermal consolidation test, the contact number of particles constantly increases during the thermally induced densification process. Under the boundary conditions in which the boundary in the y-direction has a stable 1-MPa servo stress, the normal contact force increases with an increasing number of temperature cycles, and the corresponding anisotropy level becomes more obvious. However, the principal stress direction of the normal contact force is not affected by the temperature cycles. The distribution of the specimen's force chain gradually changes to a mainly vertical distribution owing to the degradation of lateral stress.
(3) In the thermal direct shear simulation, the temperature cycle reduces not only the specimen's shear strength but also the specimen's dilatancy characteristics during the direct shear test. It is found that the temperature cycle does not affect the inhomogeneity of the specimen after reaching a steady-state in the direct shear test, through analysis of the following parameters: average normal contact force f 0 , the main direction of the normal contact force θ n, and anisotropy coefficients of the normal contact force an. In addition, it can be concluded that the high cycle number has a significant weakening effect on the formation of the shear band, through analysis of the temperature cycle's effects on force chain distribution.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).