Weijing Dai

Dorian Hanaor

Yixiang Gan
email: yixiang.gan@sydney.edu.au.

The effects of packing structure on the effective thermal conductivity of granular media: A grain scale investigation

Keywords: granular media, packing structure, grain-scale, effective thermal conductivity, discrete element method

Structural characteristics are considered to be the dominant factors in determining the effective properties of granular media, particularly in the scope of transport phenomena. Towards improved heat management, thermal transport in granular media requires an improved fundamental understanding. In this study, the effects of packing structure on heat transfer in granular media are evaluated at macro-and grain-scales. At the grain-scale, a gas-solid coupling heat transfer model is adapted into a discrete-element-method to simulate this transport phenomenon. The numerical framework is validated by experimental data obtained using a plane source technique, and the Smoluschowski effect of the gas phase is found to be captured by this extension. By considering packings of spherical SiO2 grains with an interstitial helium phase, vibration induced ordering in granular media is studied, using the simulation methods developed here, to investigate how disorder-to-order transitions of packing structure enhance effective thermal conductivity. Grain-scale thermal transport is shown to be influenced by the local neighbourhood configuration of individual grains. The formation of an ordered packing structure enhances both global and local thermal transport. This study provides a structure approach to explain transport phenomena, which can be applied in properties modification for granular media.

Introduction

Effective thermal management and waste heat dissipation are keys to improving efficiency in diverse industrial applications [START_REF] Peng | Thermal safety of lithium-ion batteries with various cathode materials: A numerical study[END_REF][START_REF] Markl | Multiscale Modeling of Powder Bed-Based Additive Manufacturing[END_REF][START_REF] Gan | Thermal Discrete Element Analysis of EU Solid Breeder Blanket Subjected to Neutron Irradiation[END_REF]. Towards this end, granular media are extensively used and are often studied by researchers in disciplines of thermal science and engineering [START_REF] Richard | Slow relaxation and compaction of granular systems[END_REF][START_REF] Chen | Heat transfer enhancement in dense suspensions of agitated solids. Part I: Theory[END_REF][START_REF] Abyzov | Effective thermal conductivity of disperse materials. I. Compliance of common models with experimental data[END_REF][START_REF] Xie | Heat transfer characteristics of silica aerogel composite materials: Structure reconstruction and numerical modeling[END_REF].

Understanding and managing heat transfer processes in granular media necessitates the accurate evaluation of their effective thermal conductivity. In a broader context, granular media are regarded as multi-phase substances of which the effective thermal conductivity can be described by the Hashin-Shtrikman bound, the Maxwell's equation and the effective medium theory [START_REF] Carson | Thermal conductivity bounds for isotropic, porous materials[END_REF].

However, these methods do not provide sufficient mechanistic insights, motivating further fundamental studies into thermal transport through granular media.

For general binary solid-gas or solid-fluid granular media, major heat transfer pathways include:

(1) heat conduction through solid particles, (2) heat transfer through solid-solid interfaces, [START_REF] Gan | Thermal Discrete Element Analysis of EU Solid Breeder Blanket Subjected to Neutron Irradiation[END_REF] heat transfer through gas/fluid films near particle surfaces, [START_REF] Richard | Slow relaxation and compaction of granular systems[END_REF] heat radiation between solid surfaces in solid-gas cases and from solid surfaces into nearby liquids in solid-fluid cases in high temperature condition (𝑇 > 450K), [START_REF] Chen | Heat transfer enhancement in dense suspensions of agitated solids. Part I: Theory[END_REF] convective heat transfer between solid and mobile gas/fluid, [START_REF] Abyzov | Effective thermal conductivity of disperse materials. I. Compliance of common models with experimental data[END_REF] natural convection in solids and fluids of high Rayleigh number, [START_REF] Xie | Heat transfer characteristics of silica aerogel composite materials: Structure reconstruction and numerical modeling[END_REF] heat transfer by macroscopic fluid flow [START_REF] Yagi | Studies on effective thermal conductivities in packed beds[END_REF], and (8) granular convective heat transfer due to the motion of solid [START_REF] Chen | Heat transfer enhancement in dense suspensions of agitated solids. Part I: Theory[END_REF]. Analytically, representative geometry methods have been widely used [START_REF] Zehner | Thermal conductivity of packings at moderate temperatures[END_REF][START_REF] Hsu | Modified Zehner-Schlunder models for stagnant thermal conductivity of porous media[END_REF], which combine contributions from these aforementioned heat pathways according to the corresponding artificial geometries. In bi-phasic systems, cylinders [START_REF] Zehner | Thermal conductivity of packings at moderate temperatures[END_REF][START_REF] Hsu | Modified Zehner-Schlunder models for stagnant thermal conductivity of porous media[END_REF] and cubic [START_REF] Hsu | A Lumped-Parameter Model for Stagnant Thermal Conductivity of Spatially Periodic Porous Media[END_REF] have been used to determine the relative contributions of each pathway. To account for the dissimilarity between the representative geometries and the meso-scale structure of granular media, including grain size distribution, contact areas, and porosity, adjustable parameters are frequently applied. But such parameters are generally too coarse or artificial to meaningfully describe granular packing structures. Rather than using representative geometries, unit cell methods that derive effective thermal conductivity by integrating grain-scale contact unit cells [START_REF] Batchelor | Thermal and Electrical Conduction Through a Granular Material[END_REF][START_REF] Gusarov | Model of thermal conductivity in powder beds[END_REF][START_REF] Bahrami | A Compact Model for Spherical Rough Contacts[END_REF][START_REF] Slavin | Theoretical model for the thermal conductivity of a packed bed of solid spheroids in the presence of a static gas, with no adjustable parameters except at low pressure and temperature[END_REF][START_REF] Dai | Influence of gas pressure on the effective thermal conductivity of ceramic breeder pebble beds[END_REF] which are commonly formed two adjacent grains have also been developed. One advantage of unit cell methods is that they facilitate the incorporation of grain-scale heat transfer mechanisms along with appropriate abstractions of packing structure, such as coordination number [START_REF] Batchelor | Thermal and Electrical Conduction Through a Granular Material[END_REF][START_REF] Dai | Influence of gas pressure on the effective thermal conductivity of ceramic breeder pebble beds[END_REF]. But similar to the representative geometry methods, it is difficult to define irregular packing structures using quantifiable parameters [START_REF] Van Antwerpen | A review of correlations to model the packing structure and effective thermal conductivity in packed beds of mono-sized spherical particles[END_REF]. With sufficient computational power, discrete element methods become favourable to study realistic irregular packing structures [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF][START_REF] Luding | A discrete model for long time sintering[END_REF][START_REF] Zhu | Discrete particle simulation of particulate systems: A review of major applications and findings[END_REF] for heat transfer purpose. One solution is to generate realistic packing structures to justify and improve the accuracy of model parameters used in the analytical methods [START_REF] Gan | Computer simulation of packing structure in pebble beds[END_REF]. The other way is to directly simulate the heat transfer process by implementing the apt grain-scale thermal interactions [START_REF] Gan | Discrete element modelling of pebble beds: With application to uniaxial compression tests of ceramic breeder pebble beds[END_REF][START_REF] Vargas | Heat conduction in granular materials[END_REF][START_REF] Feng | Discrete thermal element modelling of heat conduction in particle systems: Basic formulations[END_REF][START_REF] Moscardini | Discrete element method for effective thermal conductivity of packed pebbles accounting for the Smoluchowski effect[END_REF], such as the models developed in the unit cell methods [START_REF] Batchelor | Thermal and Electrical Conduction Through a Granular Material[END_REF][START_REF] Cheng | Evaluation of effective thermal conductivity from the structure of a packed bed[END_REF][START_REF] Yun | Three-dimensional random network model for thermal conductivity in particulate materials[END_REF][START_REF] Zhai | Interfacial electro-mechanical behaviour at rough surfaces[END_REF][START_REF] Bahrami | Effective thermal conductivity of rough spherical packed beds[END_REF].

Although the well-established heat transfer theories of continuum media have been successfully adapted to describe granular media without missing the characteristics of packing structure [START_REF] Yu | Self-Assembly of Particles for Densest Packing by Mechanical Vibration[END_REF][START_REF] Dong | Role of Interparticle Forces in the Formation of Random Loose Packing[END_REF][START_REF] An | Quasi-universality in the packing of uniform spheres under gravity[END_REF] by the discrete element methods, how to explain the effects of these characteristics is still unclear. By employing the grain-scale packing fraction as one of the variables [START_REF] Cheng | Evaluation of effective thermal conductivity from the structure of a packed bed[END_REF][START_REF] Wu | Effect of scale on the modeling of radiation heat transfer in packed pebble beds[END_REF][START_REF] Feng | An accurate evaluation of geometric view factors for modelling radiative heat transfer in randomly packed beds of equally sized spheres[END_REF] to derive heat transfer models, the importance of grain-scale packing structure is proved.

Nonetheless, effects of other structural characteristics remain inconclusive. It has been shown

that the packing structure is of equivalent importance as packing fraction [START_REF] Tsotsas | Thermal conductivity of packed beds-A review[END_REF], but the question regarding which structural indicators are appropriate for the quantitative characterisation of granular packings in heat transfer remains unanswered. Attempts has been made to investigate the influence of grain-scale stress network on thermal transport [START_REF] Watson | Stress effects on the conductivity of particulate beds[END_REF], but the link between such networks and granular packing structure requires further study [START_REF] Lim | Topology of force networks in granular media under impact[END_REF]. Relationships have been examined between effective mechanical properties and the grain-scale order level of the packing structures [START_REF] Goodrich | Solids between the mechanical extremes of order and disorder[END_REF][START_REF] Hanifpour | Mechanical Characterization of Partially Crystallized Sphere Packings[END_REF]. Whether similar relationships can be established in the heat transfer context

has not yet been discussed. The disorder-to-order transition of granular media, in which regular packing structures are formed, has been recognised as one of the major events in agitated states [START_REF] Reimann | X-ray tomography investigations of mono-sized sphere packing structures in cylindrical containers[END_REF]. It is therefore worth examining in depth the relationship between the heat transfer and the order of packing structure, which can further reflect the effects of external agitation on the effective thermal conductivity.

In the current study, the open source discrete element software LIGGGHTS [START_REF] Kloss | Models, algorithms and validation for opensource DEM and CFD-DEM[END_REF] is applied in a modified approach to investigate packing and heat transfer in a generalised static solid-gas granular media. It is demonstrated that disorder-to-order transitions induced by vibration enhance the effective thermal conductivity of granular media as ordered packing structures lead to higher grain-scale thermal conductivity.

Grain-scale heat transfer model

In this work, a modified Batchelor & O'Brien model is employed to describe inter-granular heat transfer. Contact units comprising two contacting or separated hemispheres form the basis of the grain-scale model. Here, with finite element analysis, a power function is numerically derived to extend the original method, including the correction of overestimation and the addition of the Smoluschowski effect.

Basic solution

The original Batchelor-O'Brien (BOB) model only considers granular media with large ratio 𝛼 of solid thermal conductivity 𝑘 s to gas thermal conductivity 𝑘 g , i.e., 𝑘 s /𝑘 g ≫ 1 [START_REF] Batchelor | Thermal and Electrical Conduction Through a Granular Material[END_REF]. So the temperature inside individual grains can be assumed as constant. The heat transfer equations determined by this condition have been numerically solved according to the separation distance ℎ or the radius of contact surface 𝑟 cont , as depicted in Figure 1. When two grains are separated with ℎ ≥ 0, heat conductance 𝐻 cont of one unit is

𝐻 cont = O 𝜋 𝑘 g 𝑟 #$ ln 𝛼 2 𝜆 < 0.01 𝜋 𝑘 g 𝑟 #$ ln V1 + X * 2
𝑅 cont = 1 𝐻 cont ⁄ .
However, in realistic scenarios the condition of 𝑘 s /𝑘 g ≫ 1 is rarely satisfied, e.g., glass vs. air, so thermal contact resistors tend to be underestimated. To make the model applicable in more practical scenarios, the temperature variation inside individual grains cannot be neglected. In the modified BOB model, two cylindrical thermal resistors 𝑅 cyl. # &$ of identical 𝑘 s are serially connected to the thermal contact resistor [START_REF] Yun | Three-dimensional random network model for thermal conductivity in particulate materials[END_REF], compensating for the non-inclusion of the thermal resistance of solid hemispheres in each contact unit. Finally, the total thermal resistance of one contact unit is formulated as 𝑅 total = 𝑅 cont + 𝑅 cyl. # + 𝑅 cyl. $ . In this way, the thermal resistance of cylinders is calculated according to their geometry as

𝑅 cyl. # or $ = X Y orZ l u s X eff 2 . 2-3
where the 𝑟 eff is the radius of the cylinders and is further defined as 𝑟 eff = 𝜒 min c𝑟 # , 𝑟 $ d for a more general poly-dispersed scenario. Therefore, the determination of 𝜒 becomes the essential step towards ensuring the accuracy of the modified (BOB) model. In reported methods, 𝜒 is empirically fitted according to experimental data [START_REF] Moscardini | Discrete element method for effective thermal conductivity of packed pebbles accounting for the Smoluchowski effect[END_REF][START_REF] Yun | Three-dimensional random network model for thermal conductivity in particulate materials[END_REF][START_REF] Peeketi | Effective thermal conductivity of a compacted pebble bed in a stagnant gaseous environment: An analytical approach together with DEM[END_REF]. Here, a finite element method is introduced to numerically calculate 𝜒. In addition, a flat wall can be approximated as a sphere with infinitely large radius for the grains pairing with flat wall boundaries, so 𝑟 #$ is replaced by

𝑟 #+ = lim X x →z (2 𝑟 # 𝑟 + 𝑟 # + 𝑟 + ⁄ ) = 2𝑟 # .
The other principles are still applied.

Determination of grain-scale parameters

As discussed in the previous section, 𝜒 is the only parameter to be determined for intergranular contacts. In the literature, 𝜒 = 0.5 is used for 𝛼 = 120 [START_REF] Yun | Three-dimensional random network model for thermal conductivity in particulate materials[END_REF] and 𝜒 = 0.71 is used for 𝛼 = 10 [26],

both giving good agreement with experimental results. But it is prohibitively time consuming to collect experimental data across the entire range of 𝛼 values, not to mention the associated measurement uncertainty. Therefore, an estimation of 𝜒 is useful for prediction purposes. Here, finite element analysis (using ABAQUS software) is employed [START_REF] Dai | Influence of gas pressure on the effective thermal conductivity of ceramic breeder pebble beds[END_REF]45] to evaluate 𝜒. In finite element simulations, contact units are assumed to be cylindrical, so axisymmetric modelling is used to reduce computational cost.

Two representative quadrants (solid) of 𝑘 s = 100 W/mK with the interstitial part (gas) of 𝑘 g = 𝑘 s 𝛼 ⁄ form separated and contacting contact units as shown in the left and right of Figure 2, respectively. The radii of the quadrants are set to be 𝑟 c = 1 m and the separation distance ℎ and the contact radius 𝑟 cont are varied to generate different geometries for simulations. Since material properties have no size-dependence, the exact dimensions can be later normalized. In the contacting scenario, the contacting surface is created by truncating these quadrants. The difference of the remaining height between the Hertzian quadrants and the current one is negligible (less than 0.1%) if the contact units are not severely compressed (indentation smaller than 0.2%). A constant temperature difference ∆𝑇 of 1K is imposed between the top and bottom surfaces of each contact unit while the axial side boundaries are set to adiabatic. The interfaces between solid and gas as well as between solid and solid are tied so no temperature difference exists. 𝛼 is varied from 10 6 to 4 because the original BOB model gives a contradictory trend when 𝛼 is larger than 4. 𝑟 cont /𝑟 c is varied from 0.0477 to 0.002 and ℎ/𝑟 c is varied from 0.000002 to 0.2.

Examples of the simulated temperature fields are shown in Figure 2. The total heat flow in the simulation 𝐻 FEM passing through the top surface is extracted at steady-state and compared with the heat flux calculated based on the original BOB model 𝐻 BOB . In this way the 𝜒 is evaluated by In this work, the finite element analysis is considered to be more accurate than the original BOB model due to its capability to simulate temperature variation within grains. Both heat flows, 𝐻 BOB and 𝐻 FEM , are further nondimensionlised by a factor of (∆𝑇𝑘 s 𝑟 # ) in the following discussion.

𝜒 = }~∆ • € FEM - ∆• € BOB • l u s X Y 2 -1 . 2-4
According to (a) and (b) in Figure 3, the original BOB model shows fairly good accuracy for 𝛼 ≥ 10000, but the overestimation turns out to be large when 𝛼 ≤ 1000. Therefore, the evaluation of 𝜒 is mainly conducted in the 𝛼 range from 1000 to 4 which is common practice, as shown in (c) and (d) in Figure 3. For the contacting units, 𝜒 becomes stable against 𝑟 cont as 𝛼 ≤ 200, implying the compression on the contact units has limited influence on 𝜒 in this condition.

Although the model underestimates the heat flow and gives negative 𝜒 in the separated scenario as 𝛼 ≥ 200, it occurs at relatively large ℎ/𝑟 c and results in small heat flow. Thereby, this negative part is neglected, and 𝜒 of corresponding 𝛼 in both scenarios is comparable.

The representative 𝜒 for each 𝛼 is calculated as the mean over all 𝑟 cont /𝑟 c and ℎ/𝑟 c cases for contacting and separated spheres, respectively. It can be seen from (e) in Figure 3 that 𝜒 exhibits similar values in both scenarios. Thus, the average 𝜒 of two scenarios is considered, and finally a fitting power function is obtained to estimate 𝜒 against 𝛼 as, The heat flow calculated by the modified BOB model with the proposed correlation presents reasonable agreement with the finite element simulation results, but large 𝛼 reduces its accuracy.

𝜒 = 1.3121 𝛼 ƒ0.
The 𝜒 calculated from this function for 𝛼 = 120 is 0.528, close to 0.5 as used in the literature.

While 𝛼 = 10 case gives a larger discrepancy, 0.847 versus 0.71, probably due to the polydispersity in that research [START_REF] Moscardini | Discrete element method for effective thermal conductivity of packed pebbles accounting for the Smoluchowski effect[END_REF]. Nevertheless, this function is useful for predictive purposes and the modified (BOB) model thus has good applicability, especially for conditions of 𝛼 < 1000. For 𝛼 > 10000, the original BOB model is preferable.

Inclusion of the Smoluschowski effect

The Smoluschowski effect relates to a reduced thermal conductivity through gases over scales smaller than the mean free path of gas molecules, as the result of boundary effects at a solid-gas interfaces. This effect becomes significant in granular materials where a gaseous phase sits in gaps between surfaces separated by a short distance, [START_REF] Hsu | Kinetic Theory of Parallel Plate Heat Transfer in a Polyatomic Gas[END_REF]. This effect brings about a dependence of thermal conductivity on gas pressure and grain size [START_REF] Dai | Influence of gas pressure on the effective thermal conductivity of ceramic breeder pebble beds[END_REF][START_REF] Moscardini | Discrete element method for effective thermal conductivity of packed pebbles accounting for the Smoluchowski effect[END_REF], and contributes to the deviation of 𝜒 as discussed in the previous section. In order to improve the accuracy of the modified model, the Smoluschowski effect is implemented using the method described below [START_REF] Moscardini | Discrete element method for effective thermal conductivity of packed pebbles accounting for the Smoluchowski effect[END_REF].

The influence of the Smoluschowski effect on the thermal conductivity of gases is described by

𝑘 g 3 = u g 1"2 … Kn , 2-6
where 𝛾 quantifies the efficiency of energy transfer between gas molecules and the solid surface, and is given in terms of the molar mass ratio of the solid over gas with 𝑚 † = 𝑚 7 /𝑚 g ,

𝛾 = 2 (1" ‡ ˆ)2 ƒ2.4 ‡ 2.4 ‡ ˆ , 2-7
and the Knudsen number Kn is determined by the feature size 𝐿 f of given space and the mean free path (MFP) of gas molecules as,

Kn = MFP ‰ = u B • √2 l ‹ g 2 OE g ‰ , 2-8
where the gas phase is characeterised by pressure 𝑃 g , thermodynamic temperature 𝑇, the kinetic molecular diameter 𝑑 g and the Boltzmann constant 𝑘 B .

The quantities required to determine 𝑘 g 3 are gas related except for the feature size of space 𝐿 f .

The geometry of contact units is used to obtain 𝐿 f as following. In the contacting scenario,

𝐿 f = 𝑟 # • sin -1 Ž ˆmin ˆY •ƒsin -1 Ž ˆcont ˆY •ƒ ˆmin ˆY " ˆcont ˆY sin -1 Ž ˆˆmin ˆY •ƒsin -1 Ž ˆcont ˆY • • + 𝑟 $ ' sin -1 ' ˆmin ˆZ "ƒsin -1 ' ˆcont ˆZ "ƒ ˆmin ˆZ " ˆcont ˆZ sin -1 ' ˆmin ˆZ "ƒsin -1 ' ˆcont ˆZ "
".

2-9

and for separated particles,

𝐿 f = ℎ + 𝑟 # • sin -1 Ž ˆmin ˆY •ƒ ˆmin ˆY sin -1 Ž ˆmin ˆY • • + 𝑟 $ ' sin -1 ' ˆmin ˆZ "ƒ ˆmin ˆZ sin -1 ' ˆmin ˆZ "
".

2-10

𝑟 min is set as the minimum between 𝑟 # and 𝑟 $ in one unit for a more general purpose to include poly-dispersed granular media.

For the contacting grain-wall pairs,

𝐿 f = 𝑟 # • • 2 ƒsin -1 Ž ˆcont ˆY •ƒ1" ˆcont ˆY • 2 ƒsin -1 Ž ˆcont ˆY •
•.

2-11

and the separated grain-wall pairs,

𝐿 f = 𝑟 # Ž • 2 ƒ1 • 2
•.

2-12

Eventually, 𝛼 = 𝑘 g 3 /𝑘 s is derived and used for the modified BOB model.

Simulation protocol

The grain-scale heat transfer model is coded in C++ within LIGGGHTS framework and is provided in supporting information. For each simulation, a granular medium is generated by the free-falling method in a given volume with specified boundary conditions. Following particle settling under gravitational force onto the flat bottom surface of the container, a flat upper surface is added and is used to apply mechanical compression. The upper surface is set as a heat source while the bottom wall is treated as a heat sink with a constant temperature gap ∆𝑇 gap of 1K between these two boundaries. When a target loading is achieved, the granular medium is frozen, the walls are locked, and adiabatic boundary conditions are imposed in lateral directions.

The granular medium is then heated until the simulation reaches steady states where the fractional difference between inflow heat (top) and outflow heat (bottom) is smaller than 10 -6 .

Finally, the effective thermal conductivity 𝑘 eff is expressed as

𝑘 eff = € inflow ‰ axis -cross-section ×∆• grap . 2-13
where 𝐻 inflow indicates the inflow heat, and 𝐿 axis and 𝐴 cross-section are respectively the height and the cross-sectional area of the granular medium. Alternatively, 𝑘 eff can be calculated also as

𝑘 eff = € inflow -cross-section × 1 ∆• gradient . 2-14
with ∆𝑇 gradient representing the temperature gradient inside the granular medium. In general, by binning the grains according to their height, a linear function can be fitted based on the mean temperature of each bin against the representative height of each bin. Thus, ∆𝑇 gradient can be calculated as the slope of that linear function. Compared with the former expression, using ∆𝑇 gradient can exclude the wall-grain contacts and gives the more inherent 𝑘 eff of the granular medium. However, the former one is considered more appropriate for experimental validation as practical measurement always includes wall-grain contacts.

Experimental validation and discussion

The HotDisk system built using a transient plane source technique [START_REF] Gustavsson | Thermal conductivity, thermal diffusivity, and specific heat of thin samples from transient measurements with hot disk sensors[END_REF][START_REF] He | Rapid thermal conductivity measurement with a hot disk sensor[END_REF] is employed to experimentally measure the effective thermal conductivity of granular media. Two kinds of grains, 2 mm glass (𝑘 s = 1.2 W/mK) beads and 2 mm 304 stainless steel (𝑘 s = 16.2 W/mK) beads, are used separately to form solid-gas media in air under ambient conditions of atmospheric pressure and room temperature (𝑘 g = 0.026 W/mK). The Smoluschowski effect is described by the correlations shown in Table 1. A load-bearing sample container is designed to hold these granular media. In this way, the mechanical influence on effective thermal conductivity of granular media can be investigated.

Measurement results, shown in Figure 4, reveal that the thermal transport of the two granular systems responds differently to mechanical loading. The granular medium consisting of glass beads is insensitive to mechanical loading while the 𝑘 eff of the packed steel grains increases gradually with load. This difference can be understood in light of the fact that steel has a thermal conductivity more than ten times higher than that of glass. Larger load creates larger intergranular contact surfaces, which results in greater heat conduction pathways, with these having greater significance in the more thermally conductive steel. Discrete element simulation results are quantitatively similar, demonstrating that the applicability of the present heat transfer model.

Table 1 Gas properties used in this work Gas type Properties Helium (He) 𝑘 g / (W/mK): 3.366 × 10 -3 𝑇 0.668 (𝑇 in K).

𝑑 g / (m):

2.15×10 -10 . 𝑚 g / (g/mol): 4.

Air

𝑘 g / (W/mK): -10 -11 𝑇 3 -4×10 8 𝑇 2 + 8×10 5 𝑇 + 0.0241 (𝑇 in ℃).

𝑑 g / (m): 3.66×10 -10 . 𝑚 g / (g/mol): 28.96. In addition, the role of the Smoluschowski effect and poly-dispersed particle size distributions are also examined with reference to published experimental results. A poly-dispersed granular medium of lithium orthosilicate Li4SiO4 [START_REF] Pupeschi | Effective thermal conductivity of advanced ceramic breeder pebble beds[END_REF] (𝑘 s = 2.57 W/mK at 20 ̊C and 𝑘 s = 2.15 W/mK at 500 ̊C) with helium filling (𝑘 g = 0.15 W/mK at 20 ̊C and 𝑘 g = 0.26 W/mK at 500 ̊C) is fed to the discrete element simulation. 𝑘 eff of the simulated media shows results closely matching experimental data including gas pressure dependence. It should be pointed out that this heat transfer simulation slightly overestimates 𝑘 eff when compared with experimental data [START_REF] Abou-Sena | Effective Thermal Conductivity of Lithium Ceramic Pebble Beds for Fusion Blankets: A Review[END_REF] despite neglecting thermal radiation.

Effective thermal conductivity vs. packing structure

The role of material properties and grain-grain interactions in overall thermal transport have been examined through the implementation of the discrete element simulation reported here, and have also been widely reported in the literature. Therefore, the following sections focus on establishing relationships between packing order and effective thermal conductivity.

Packing structure transition by vibration

It has been demonstrated that mechanical vibration can effectively change the packing structure of granular media [START_REF] Barker | Transient phenomena, self-diffusion, and orientational effects in vibrated powders[END_REF][START_REF] An | Effect of vibration condition and inter-particle frictions on the packing of uniform spheres[END_REF]. For granular media settling freely under gravity, random disordered packing structures form with a few clusters of hexagonal-close-packed (HCP) and face-centredcubic (FCC) packings. By vibrational agitation of granular media disorder-to-order transitions are provoked. Grains gradually rearrange themselves into structured clusters to achieve dynamic efficiency [START_REF] Reimann | X-ray tomography investigations of mono-sized sphere packing structures in cylindrical containers[END_REF][START_REF] Dai | Modes of wall induced granular crystallisation in vibrational packing[END_REF]. Not only can the packing fraction be varied, but also the order level of the packing structure significantly increases during vibration. The vibration protocol used in this study is as follows: Granular media with a constant number of grains are generated by the freefalling method within finite cylinders of varying diameters under gravitational force using LIGGGHTS. These packings are then subjected to continuous mechanical vibration by sine wave agitation of the bottom surface. As a result, grains are excited and relocate into new positions, forming structured clusters that grow and finally reach quasi-equilibrium states [START_REF] Dai | Modes of wall induced granular crystallisation in vibrational packing[END_REF]. Three types of granular media are simulated and their geometrical dimensions are summarised in Table 2.

For each geometry, five samples are extracted at different vibration durations for subsequent heat transfer simulation. To quantify this disorder-to-order transition, the structural index 𝑆 6 is utilised to measure the rotational symmetry of the neighbour configuration of individual grains. In such measurement, the neighbour configuration of each grain is first represented by a set of vectors pointing from a given grains to its neighbours, i.e., 12 nearest grains in current study, and then the cosine similarity between the vector set of a particular grain and those vector sets of its neighbours is used to evaluate the resemblance between the neighbour configurations of these grains. Thus, the order level of grain-scale packing structure is positively described by the aforementioned resemblance, with values ranging from 0 to 12 [START_REF] Goodrich | Solids between the mechanical extremes of order and disorder[END_REF][START_REF] Dai | Modes of wall induced granular crystallisation in vibrational packing[END_REF][START_REF] Kumar | Bond-orientational analysis of hard-disk and hard-sphere structures[END_REF][START_REF] Russo | The microscopic pathway to crystallization in supercooled liquids[END_REF]. Voronoi tessellation [START_REF] Rycroft | Analysis of granular flow in a pebblebed nuclear reactor[END_REF] is employed to obtain the effective volume 𝑉 V and packing fraction 𝜑 V of individual grains. In this way, the mean values of 𝑆 6 and 𝜑 V , excluding surface grains, can be recognised as the global characteristics of a granular medium. As shown in Figure 5, the mean values of both of 𝑆 6 and 𝜑 V monotonically increase with vibration, which implies a cooperative process between densification and ordering in granular media. The reason for the sudden increase of 𝜑 V in D50 can be explained by the geometrical effect in the ordering process [START_REF] Dai | Modes of wall induced granular crystallisation in vibrational packing[END_REF]. The large cylinder diameter promotes the growth of perfect HCP and FCC clusters, while the small one causes distortion in the ordered clusters with imperfection. By considering individual grain positions, the spatial resolution of 𝑆 6 can be further revealed.

Due to the symmetry of the cylindrical containers, two-dimensional (2D) radial sections are studied and sub-domains are created by mapping a 2D grid with grain diameter resolution onto these sections. By sweeping the cylindrical spaces with the 2D sections, grains are mapped into the relevant sub-domains according to their radial distance 𝑥 radial and height 𝑧, and then the mean 𝑆 6 of each sub-domain can be obtained to yield intensity maps of 𝑆 6 . Spatial patterns of these maps are found to be nearly axisymmetric. A typical evolution of 𝑆 6 distribution is shown in Figure 6 presenting a spatially inhomogeneous order development.

Enhanced effective thermal conductivity by vibration

To examine the role of packing structure in thermal conductivity we consider a medium consisting of glass (SiO2) grains filled with a gas phase of helium (He). By varying the helium pressure in the simulation, 𝛼 can be varied over a large range for a given solid phase. According to Figure 7 This global observation indicates that neither packing fraction nor structure order can determine 𝑘 eff alone. The collective processes of densification and ordering are the reason for the enhancement of 𝑘 eff introduced by vibration. But, by decreasing 𝛼, the prominence in such processes shifts from packing fraction to structure order.

Grain-scale structure and heat transfer

It has been shown that 𝑘 eff of granular media is positively correlated to both 𝜑 V and 𝑆 6 at the macroscale, but mechanisms at finer scales are insufficiently understood. As demonstrated in Figure 6, the structural transition brings about substantial and inhomogeneous grain-scale ordering. Whether this improvement enhances heat transfer in granular media and if so, how this enhancement works remain elusive. In order to address this interplay between grain-scale structure and thermal transport, analysis of heat transfer at the grain-scale is performed.

Analytically, a heat flux density vector 𝑯 density across the surface of an individual grain in a granular medium can be presented as [START_REF] Batchelor | Thermal and Electrical Conduction Through a Granular Material[END_REF],

𝑯 density = 𝑛 ∑ 𝒙 # 𝐻 # # . 4-1
where 𝐻 # is the heat flow passing through each contact, 𝒙 # is the vector pointing from the centroid of the particular grain toward the corresponding contact point and 𝑛 is the grain number density, which can also be calculated at the grain-scale as the reciprocal of the Voronoi volume [START_REF] Rycroft | Analysis of granular flow in a pebblebed nuclear reactor[END_REF], 1/𝑉 V of the corresponding grain. Further, it has also been proved that the temperature gradient vector of each grain is equal to the global temperature gradient vector ∆𝑻 gradient of the entire granular medium [START_REF] Batchelor | Thermal and Electrical Conduction Through a Granular Material[END_REF]. Based on this consideration, an expression is validated for each grain,

𝑯 density = 𝒌 * ∆𝑻 gradient = 𝑛 ∑ 𝒙 # 𝐻 # # . 4-2
with 𝒌 * being the second rank thermal conductivity tensor for each grain. Therefore, 𝒌 * of individual grains can be derived if the corresponding neighbour configuration and the heat flows through these neighbouring grains are known. Such information can be easily extracted from the discrete element simulation, which facilitates this grain-scale heat transfer analysis. Besides, the 20 sum of heat flows of each grain is 0 at steady-state, ∑ 𝐻 # # = 0, so by separately summing inflow or outflow heat for each grain, the scalar quantity of total heat flow 𝐻 flow through one grain can also be calculated. In such measurements, two aspects of heat transfer at the grain-scale can be investigated, the magnitude of the heat flow and the directional thermal conductivity. Thus, the fitted ∆𝑇 gradient in Eqn. Finally, it is necessary to point out that it is challenging to untangle the effects of packing ordering and densification in granular media. Nonetheless, an increase in local ordering can be clearly attributed to two results. One is the enhancement of grain-scale thermal conductivity along with the increase of packing fraction. The second is the reduction of anisotropy, improving the alignment of grain thermal conductivity with the principle heat transfer direction. Further, the effects of ordering would be more pronounced in systems where not all grains are in contact e.g. cemented granular media and grain-filler composites. To check this point, the packing structures of those granular media are reused to form artificial media, which can possibly be considered as grain-filled composites. In these composites, the grains are fixed at the identical positions with the order level maintained (disorder -initial state and well-order -longest vibrated state, corresponding to Figure 5), but the radii of the grains are varied to adjust the packing fraction Firstly, in spite of the dissimilar varying tendencies of 𝑘 eff against 𝜑 in Figure 11, a comparable shape of S can be identified. Although the S shape is flattened in the condition of 𝛼 = 10, a transition step is still observed. This invariant transition reflects the formation and loss of contacts between grains in corresponding to the change of packing fraction, i.e., the grain radius.

which is the reason for the sharp transition in the condition of 𝛼 = 1000. Regarding to the effect of the structural order, it works differently in according to 𝛼. In the large 𝛼condition, the ordered structure has relative uniform separation distance between grains, making the plunge of 𝑘 eff become sharper than the disorder structure. But higher coordination number of the ordered structure improves the possibility of contact formation, resulting in higher 𝑘 eff with the packing fraction increasing. The mergence as the reduction of packing fraction in this scenario can be considered as the loose dispersion limit of the Maxwell model [START_REF] Carson | Thermal conductivity bounds for isotropic, porous materials[END_REF] where the structural order is least influential. However, in the condition of small 𝛼, the preference between packing order and packing fraction further depends on the detailed packing structure that is influenced by the geometry. Due to the presence of the compensation by small 𝛼, the slender medium (D30) of ordered structure generally exhibits higher 𝑘 eff in the present packing fraction range. On the contrary, the flat medium (D50) presents similar change of the preference to the one of 𝛼 = 1000. Such difference has been attributed to the poorly ordered side region, which can also be recognised as the segregation of disordered region. This segregation introduces disorder of larger dimension beyond the grain-scale, so out of the scope of current study.

Conclusion

In with mono-dispersed and poly-dispersed granular media, while deviation becomes larger when grain size distributions are unknown. As thermal radiation and convection have been ignored, the applicability of this framework is optimal for granular media with stagnant gas phases at room temperature.

Granular media comprising glass and helium with different packing structures are applied in this heat simulation framework to study how the order characteristics of packing structure affect the effective thermal conductivity of these media. By investigating thermal properties and order characteristics at the grain-scale, it was demonstrated that the structural disorder-to-order transition enhances the overall thermal conductivity 𝑘 eff . Locally, the increase of thermal conductivity in the principal direction, 𝑘 FF , and the decrease of anisotropy 𝜃 F , are further identified as outcomes of the structural transitions, which result from grain-scale ordering processes. The heat flow through individual grains, as indicated by 𝐻 flow , is found to correlate strongly with the global structural variation but does not exhibit a clear relationship with local structure. The two grain-scale effects are further influenced by the 𝑘 s /𝑘 g ratio. The findings here provide a potential method to manage heat transfer in granular media through manipulation of packing structure and provide a conceptual framework for the broader analysis of transport phenomena in granular media and heterogeneous media in general.

Figure 1

 1 Figure 1 Separated (left) and contacting (right) pairs of grains.

Figure 2

 2 Figure 2 Separated (left) and contacting (right) geometries coloured with steady state temperature field in the finite element simulation.

19

 19

  contacting (d) separated (e) 𝜒 versus 𝛼 (f) uncertainty test Figure 3 (a) and (b) compare the heat flow between the finite element simulation (filled circles) and the original Batchelor & O'Brien model (solid lines), coloured according to 𝛼. (c) and (d) plot 𝜒 for each 𝛼 and 𝑟 cont /𝑟 c pair as well as 𝛼 and ℎ/𝑟 c pair with joining lines, respectively. (e) gives the correlation between 𝜒 and 𝛼 and (f) compares the computed heat flow between the modified model and the simulated result.

Figure 4

 4 Figure 4 Top -Experimental configuration. Bottom -Comparison between simulation results and the experimental measurement in this work (left) and previously published studies (right).

Figure 5

 5 Figure 5 𝜑 V (-left) and 𝑆 6 (-right) increase as the duration of vibration is extended, noticing that the horizontal axis is arbitrary.

Figure 6

 6 Figure 6 Evolution of 𝑆 6 spatial patterns in D40 with vibration duration (from left to right).

Figure 7

 7 Figure 7 (a) The variation of 𝛼 with helium gas pressure in SiO2-He granular media consisting of 2.3-mm-diameter SiO2 grains at room temperature. (b) and (c) show the correlations of (𝑘 eff , 𝜑 V ) and (𝑘 eff , 𝑆 6 ) in the minimum 𝛼 (≈ 1000), respectively. (d) and (e) show similar correlations of the maximum 𝛼 (< 10).

( 2 -

 2 14) is used and resulted in ∆𝑻 gradient = 0, 0, ∆𝑇 gradient ¡ due to the imposed temperature boundary condition. Because of the uni-directional nature of the temperature gradient, only the 𝑧 component of thermal conductivity tensor 𝒌 * is needed, giving 𝑘 F¢ , 𝑘 F£ , 𝑘 FF ¡. The principle component 𝑘 FF is of most interest, and the anisotropy 𝜃 F = ArcCos(𝑘 FF /|𝒌 * | F ) of this 𝑧 component of 𝒌 * is examined here.

Figure 8

 8 Figure 8 Spatial evolution of 𝑘 FF /𝑘 s (top row), 𝜃 (middle) and 𝐻 flow (bottom) of the D40 granular media.

Figure 9 -Figure 10 (

 910 Figure 9 -(𝑆 6 , 𝑘 FF ), middle -(𝑆 6 , 𝜃 F ) and right -(𝜑 V , 𝑘 FF ) plots of D30, D40 and D50 (from top to bottom respectively), created to according to the mapping operation, in the scenario of 𝛼 = 100.

  alone. As shown in Figure 11, two scenarios of 𝛼 = 1000 (left) and 𝛼 = 10 (right) are presented, where the Smoluschowski effect is turned off to avoid the size dependency due to the change of grain radius.

Figure 11

 11 Figure 11 The variation of 𝑘 eff against the change of packing fraction in artificial media with disorder and well-order packing structures in 𝛼 = 1000 (left) and 𝛼 = 10 (right) scenarios.

Table 2

 2 Geometry parameters of granular media

	Diameter of grain	2.3 (mm)
	Number of grains	5000
	Cylinder Diameter/Height	30/75, 40/60, 50/25 (mm/mm)
		denoted as D30, D40, D50

  this work, a bi-phasic heat transfer model based on the Batchelor & O'Brien solutions was successfully implemented in the open source software LIGGGHTS. To overcome the inaccurate prediction of the original model when the assumption 𝑘 s /𝑘 g ≫ 1 is unsatisfied, a compensation method with only one variable 𝜒 is applied. With the help of finite element analysis of the heat transfer process of a series of contact geometries and 𝑘 s /𝑘 g ratios, an empirical fitting function is proposed to determine 𝜒 according to 𝑘 s /𝑘 g for individual contacts. For further accuracy, the Smoluschowski effect is included in the modified model. Finally, this heat transfer simulation framework is validated by experimental measurements employing a transient plane source technique as well as comparison with literature data. The simulated results show good agreement

Dai, W., Hanaor, D., & Gan, Y. (2019). The effects of packing structure on the effective thermal conductivity of granular media: A grain scale investigation. International Journal of Thermal Sciences, 142, 266-279.

A similar analysis to that presented in Figure 6 is applied to examine the spatial distributions of the thermal conductivity tensor 𝑘 FF , its anisotropy 𝜃 F and the grain level heat flow (𝐻 flow ). The spatial evolution of the ratio 𝑘 FF /𝑘 s and 𝜃 F for a D40 sample with 𝛼 ≈ 10 is shown sequentially in the top and middle rows of Figure 8, respectively. The spatial distribution of 𝑘 FF /𝑘 s correlates positively to that of 𝑆 6 in Figure 6, while the patterns of 𝜃 F exhibit an inverse trend. In contrast, the spatial distribution of non-dimensional 𝐻 flow (bottom row in Figure 8) does not show any observable trend that can be correlated to structural patterns, and exhibits only a globally increasing trend along with the overall 𝑆 6 , as vibration proceeds. This suggests that 𝑘 FF and 𝜃 F are closely related to the grain-scale structure but 𝐻 flow is less influenced.

Comparing the evolution of the spatial distributions of 𝑘 FF /𝑘 s and 𝜃 F with the evolution of 𝑆 6 , it can be seen that, at the grain-scale, ordering increases thermal conductivity and decreases anisotropy, resulting in similar inhomogeneous patterns. In fact, due to the cooperativity of the densification and ordering, the spatial patterns of 𝜑 V are nearly identical to those of 𝑆 6 . Thus, it is hard to differentiate the influences of these two processes, but ordering is considered to be precedence of densification in granular media due to the rigidity of grains. However, it is clear that the decrease of anisotropy, 𝜃 F , can be related to a greater symmetry of the neighbourhood configurations of individual grains, which is specifically measured by 𝑆 6 . Further, the bottom ordered region in the patterns of Figure 6 exhibits a higher 𝑘 FF /𝑘 eff and lower 𝜃 F in Figure 8 than the corresponding side ordered region. This contrast can be attributed to the different completeness of ordering in these two regions. As mentioned before, due to the curved cylinder wall, distortion occurs in the side region [START_REF] Dai | Modes of wall induced granular crystallisation in vibrational packing[END_REF], making the ordered packing structures deviate from perfect FCC and HCP and become looser as well as less symmetric. In summary, two grain-scale mechanisms to enhance the macroscopic 𝑘 eff are identified regarding to the ordered packing structure.

By transforming the quantitative results from the mapping operations of Figure 5 and Figure 8 into the coordinate plots shown in Figure 9, parameter correlations are further investigated. The Pearson correlation coefficient is employed to quantitatively evaluate the correlations of (𝑆 6 , 𝑘 FF ) and (𝑆 6 , 𝜃 F ) in the left and middle of Figure 9, respectively. It confirms the positive correlation between 𝑆 6 and 𝑘 FF together with the negative correlation between 𝑆 6 and 𝜃 F .