Weijing Dai

Joerg Reimann

Dorian Hanaor

Claudio Ferrero

Yixiang Gan

Modes of wall induced granular crystallisation in vibrational packing

Keywords: granular materials, packing, vibration, boundary effects, crystallisation

published or not. The documents may come

Introduction

The behaviour of granular matter subjected to shear, vibration, flow and mixing is of tremendous significance in diverse applications involving the handling and processing of materials in particulate form [START_REF] Jaeger | Physics of the Granular State[END_REF][START_REF] De Gennes | Granular matter: a tentative view[END_REF]. The evolution of granular packing is particularly important as system behaviour and effective properties are closely correlated to packing structure [START_REF] Wang | Predictions of effective physical properties of complex multiphase materials[END_REF]. Recently, granular crystallisation under vibration has drawn increasing interest [START_REF] Pouliquen | Crystallization of non-Brownian Spheres under Horizontal Shaking[END_REF][START_REF] Tsai | Internal Granular Dynamics, Shear-Induced Crystallization, and Compaction Steps[END_REF][START_REF] Carvente | Crystallization of Confined Non-Brownian Spheres by Vibrational Annealing[END_REF][START_REF] An | DEM study of crystallization of monosized spheres under mechanical vibrations[END_REF][START_REF] Saadatfar | Pore configuration landscape of granular crystallization[END_REF]. Similar to its counterpart in materials science, crystallisation in granular matter can be characterised by the formation of ordered structures. Analogous to heat in materials science, vibration acts as an energy source to agitate particles to jump around and form crystalline structures. However, as granular matter is an athermal and generally repulsive system, the underlying mechanisms differ from bonding at the atomic level. For decades, vibration has been utilised to excite granular matter for fluidisation, segregation and packing [START_REF] Dong | A novel method based on orientation discretization for discrete element modeling of non-spherical particles[END_REF][START_REF] Wu | DEM simulation of cubical particle packing under mechanical vibration[END_REF][START_REF] Tai | The flow regime during the crystallization state and convection state on a vibrating granular bed[END_REF][START_REF] D'anna | Observing brownian motion in vibrationfluidized granular matter[END_REF][START_REF] Sufian | Pore shapes, volume distribution and orientations in monodisperse granular assemblies[END_REF][START_REF] An | Quasi-universality in the packing of uniform spheres under gravity[END_REF], and is thus utilised in this study to explore how granular crystallisation occurs. Early studies show how internal dynamics of particles determine the macroscopic behaviour of entire granular systems [START_REF] Barker | Transient phenomena, self-diffusion, and orientational effects in vibrated powders[END_REF][START_REF] Zhao | Micro origins for macro behavior in granular media[END_REF][START_REF] Ratnaswamy | Evolution of solids fraction surfaces in tapping: simulation and dynamical systems analysis[END_REF]. Depending on the intensity of vibration, macroscopic phenomena manifesting during vibration are compaction [START_REF] Boutreux | Compaction of granular mixtures: a free volume model[END_REF][START_REF] Richard | Slow relaxation and compaction of granular systems[END_REF][START_REF] Yu | Self-Assembly of Particles for Densest Packing by Mechanical Vibration[END_REF][START_REF] Rosato | Microstructure evolution in density relaxation by tapping[END_REF] and convection [START_REF] Philippe | Granular Medium under Vertical Tapping: Change of Compaction and Convection Dynamics around the Liftoff Threshold[END_REF][START_REF] Lan | Convection related phenomena in granular dynamics simulations of vibrated beds[END_REF][START_REF] Ratnaswamy | Evolution of solids fraction surfaces in tapping: simulation and dynamical systems analysis[END_REF]. However, granular crystallisation lies at the intersection of these two phenomena, as it requires a certain fluidisation to facilitate particle rearrangement, yet results in a compacted state.

Experimentally, granular crystallisation has been studied by vibrating existing packings [START_REF] Carvente | Crystallization of Confined Non-Brownian Spheres by Vibrational Annealing[END_REF][START_REF] Carvente | Self-assembling of dry and cohesive non-Brownian spheres[END_REF] or by adding particles at controlled rates to horizontally or vertically vibrating boxes [START_REF] Pouliquen | Crystallization of non-Brownian Spheres under Horizontal Shaking[END_REF][START_REF] Nahmad-Molinari | Epitaxial Growth of Granular Single Crystals[END_REF]. In the former scenario, a gradual deceleration of agitation, analogous to annealing, plays an important role in the emergent crystallisation, as do pre-set templates and particular container geometries in the latter case [START_REF] Pouliquen | Crystallization of non-Brownian Spheres under Horizontal Shaking[END_REF]. The use of three-dimensional vibration, rather than uniaxial, is found to further enhance crystallisation by disrupting granular arching [START_REF] Yu | Self-Assembly of Particles for Densest Packing by Mechanical Vibration[END_REF]. The behaviour of sinusoidal-vibration driven rearrangement of granular packing is manipulated by filling rate and vibration intensity 𝛤, defined as 𝛤 = 𝐴(2𝜋𝑓) 2 /g , where g is the gravitational acceleration, and 𝐴 and 𝑓 are the vibration amplitude and frequency, respectively. Maximum crystallisation is found to occur for an intermediate value of 𝛤 [START_REF] Panaitescu | Epitaxial growth of ordered and disordered granular sphere packings[END_REF][START_REF] Yu | Self-Assembly of Particles for Densest Packing by Mechanical Vibration[END_REF][START_REF] Nowak | Reversibility and irreversibility in the packing of vibrated granular material[END_REF]. Earlier works have demonstrated the possibility to control granular packing structures through geometry parameters, cohesion and agitation, with the 3 minimisation of internal energy proposed to drive crystallisation. However, a mechanistic understanding of how granular crystallisation occurs and develops under simple and continuous vibration still remains elusive. Insights can be gained from the densification of granular matter, where mechanisms have been discussed in terms of the preservation and reconstruction of contact networks [START_REF] An | Micromechanical Simulation and Analysis of One-Dimensional Vibratory Sphere Packing[END_REF] and pore size distribution [START_REF] Richard | Analysis by x-ray microtomography of a granular packing undergoing compaction[END_REF][START_REF] Boutreux | Compaction of granular mixtures: a free volume model[END_REF].

In order to establish accurate mechanisms, numerical simulations, such as discrete element methods and Monte Carlo methods [START_REF] Dshinde | Shaking-induced crystallization of dense sphere packings[END_REF], have been implemented together with X-ray computed tomography (XCT) to examine structural evolution [START_REF] Rosato | Microstructure evolution in density relaxation by tapping[END_REF] and the appearance of crystalline structures [START_REF] Saadatfar | Pore configuration landscape of granular crystallization[END_REF]. In recent studies, the emergence of crystallisation by vibration was identified for conditions where the random close packing fraction (0.64) was exceeded [START_REF] Hanifpour | Structural and mechanical features of the order-disorder transition in experimental hardsphere packings[END_REF][START_REF] Reimann | Pebble bed structures in the vicinity of concave and convex walls[END_REF]. By examining the internal structure of various packings formed under different vibration protocols, the formation of polytetrahedral patterns and octahedral cavities was shown to impart geometrical frustration [START_REF] Francois | Geometrical Frustration in Amorphous and Partially Crystallized Packings of Spheres[END_REF][START_REF] Saadatfar | Pore configuration landscape of granular crystallization[END_REF]. Additionally, using high-speed cameras, the role of boundaries in the nucleation of crystalline regimes in two-dimensional systems has been studied across multiple scales in terms of pattern formation and densification as well as their relation to grain mobility [START_REF] Lumay | Experimental Study of Granular Compaction Dynamics at Different Scales: Grain Mobility, Hexagonal Domains, and Packing Fraction[END_REF] and energy dissipation [START_REF] Komatsu | Roles of Energy Dissipation in a Liquid-Solid Transition of Out-of-Equilibrium Systems[END_REF]. Three-dimensional systems with planar, convex and concave boundaries have been recently investigated experimentally [START_REF] Reimann | X-ray tomography investigations of mono-sized sphere packing structures in cylindrical containers[END_REF][START_REF] Reimann | Pebble bed structures in the vicinity of concave and convex walls[END_REF], demonstrating the significance of boundary geometries.

Through well-developed frameworks used to describe analogous crystallisation from glasses or gels in the domain of materials science [START_REF] Russo | The microscopic pathway to crystallization in supercooled liquids[END_REF][START_REF] Tanaka | Bond orientational order in liquids: Towards a unified description of water-like anomalies, liquid-liquid transition, glass transition, and crystallization[END_REF], the field of granular crystallisation remains ripe for further exploration, with a view towards enhancing high-value packing-dependent properties [START_REF] Goodrich | Solids between the mechanical extremes of order and disorder[END_REF][START_REF] Hanifpour | Structural and mechanical features of the order-disorder transition in experimental hardsphere packings[END_REF]. The present study explores transient states of crystallisation, using a discrete element method (DEM) to simulate dynamic behaviour of vibrated granular matter. Different container geometries were used to examine the boundary influence, with results compared to XCT data. Simulations show that in the absence of cohesion or attraction, granular matter exhibits a clear tendency to crystallise under vibration. Different modes of crystallisation were identified through the analysis on the packing structure transitions. The favoured propagation direction of crystallisation is examined here as a function of container geometry as is the resulting spatial distribution of crystallites.

Investigation setup and motivation

Discrete element method

Dynamics of granular matter subjected to external vibration is simulated by the open source software LIGGGHTS [START_REF] Kloss | Models, algorithms and validation for opensource DEM and CFD-DEM[END_REF] based on DEM. In this method, the propagation of external agitation is interpreted as the result of inter-particle collisions, causing the motion of granular particles. For an individual collision between particles 𝑖 and 𝑗 , the normal component 𝑭 n 𝑖,𝑗 and tangential component 𝑭 t 𝑖,𝑗 of the inter-particle force are described by the following equations,

𝑭 n 𝑖,𝑗 = 𝑘 n 𝑖,𝑗 𝜹 n 𝑖,𝑗 -𝛾 n 𝑖,𝑗 𝒗 n 𝑖,𝑗 , 𝑭 t 𝑖,𝑗 = 𝑘 t 𝑖,𝑗 𝜹 t 𝑖,𝑗 -𝛾 t 𝑖,𝑗 𝒗 t 𝑖,𝑗 , (1)
where 𝑘 n 𝑖,𝑗 and 𝛾 n 𝑖,𝑗 and (𝑘 t 𝑖,𝑗 and 𝛾 t 𝑖,𝑗 ) are contact stiffness and viscoelastic damping coefficient for normal and (tangential) contact, respectively. These quantities are derived from the Hertz-Mindlin contact theory, and thus depend on the instantaneous contact configuration. Here, 𝜹 n 𝑖,𝑗 is the overlap distance, 𝒗 n 𝑖,𝑗 and (𝒗 t 𝑖,𝑗 ) are the relative velocities in the normal and tangential directions, and 𝜹 t 𝑖,𝑗 is the tangential displacement vector between particles 𝑖 and 𝑗. In addition, 𝑭 t 𝑖,𝑗 is limited by the Coulomb friction limit, |𝑭 t 𝑖,𝑗 | ≤ 𝜇 |𝑭 n 𝑖,𝑗 |, in which 𝜇 is the friction coefficient.

Order characterisation

In the present study, we used void fraction distributions, coordination numbers, contact angle and radial density distributions as well as Voronoi tessellation to characterise the packings [START_REF] An | Micromechanical Simulation and Analysis of One-Dimensional Vibratory Sphere Packing[END_REF][START_REF] Reimann | X-ray tomography investigations of mono-sized sphere packing structures in cylindrical containers[END_REF][START_REF] An | Effect of vibration condition and inter-particle frictions on the packing of uniform spheres[END_REF].

However, we concentrate here on the use of bond orientation order parameters to distinguish crystalline structures and represent the transitions between ordered and disordered states. An advantage of this measure is its insensitivity to particle separation, enabling its use for the transient characterisation of moving granular matter. Together, static and dynamic measures are used to describe the structural changes across order transitions.

Bond orientation order parameters, initially defined by Steinhardt et al [START_REF] Steinhardt | Bond-orientational order in liquids and glasses[END_REF], represent the rotational symmetry of sphere assemblies as,

𝑄 𝑙𝑚 (𝑟 ⃗) ≡ 𝑌 𝑙𝑚 (𝜃(𝑟 ⃗), 𝜑(𝑟 ⃗)), (2) 5
where a bond 𝑟 ⃗ is defined as a vector that points from the centroid of a central particle to one of its neighbour particles, 𝑌 𝑙𝑚 (𝜃, 𝜑) are spherical harmonics, 𝜃(𝑟 ⃗) and 𝜑(𝑟 ⃗) are the polar and azimuthal angles of the bond in a reference spherical coordinates system, and 𝑙 and 𝑚 are integers indicating the order of spherical harmonics with the condition that 𝑙 ≥ 0 and |𝑚| ≤ 𝑙 . By averaging 𝑄 𝑙𝑚 (𝑟 ⃗) over the 𝑛 b 𝑖 closest neighbours of a central particle 𝑖, the following expression is obtained,

𝑞 ̂𝑙𝑚 (𝑖) = 1 𝑛 b 𝑖 ∑ 𝑄 𝑙𝑚 𝑖,𝑘 (𝑘) 𝑛 b 𝑖 𝑘=1 . (3)
In the current study, the number of neighbours 𝑛 b 𝑖 of a central particle 𝑖 is selected as 12 [START_REF] Russo | The microscopic pathway to crystallization in supercooled liquids[END_REF][START_REF] Dshinde | Shaking-induced crystallization of dense sphere packings[END_REF],

the largest coordination number of non-overlapping mono-sized particles, because it gives significantly different feature values between crystalline structures like BCC, HCP, FCC and even icosahedral. Finally, the local bond orientation order for particle 𝑖 is constructed as [START_REF] Kansal | Diversity of order and densities in jammed hardparticle packings[END_REF] 𝑄 𝑙 local (𝑖) ≡ (

4𝜋 2𝑙+1 ∑ |𝑞 ̂𝑙𝑚 (𝑖)| 2 𝑙 𝑚=-𝑙 ) 1/2 , (4)
𝑊 ̂𝑙(𝑖) = ∑ ( 𝑙 𝑙 𝑙 𝑚 1 𝑚 2 𝑚 3 )𝑞 ̂𝑙𝑚 1 (𝑖) 𝑞 ̂𝑙𝑚 2 (𝑖) 𝑞 ̂𝑙𝑚 3 (𝑖) 𝑚 1 ,𝑚 2 ,𝑚 3 [∑ |𝑞 ̂𝑙𝑚 (𝑖)| 2 6 𝑚=-6 ] 3 2 ⁄ , (5)
where the term in the parentheses is Wigner-3j symbol. 𝑙 = 4 and 6 are widely used due to their unambiguous value for regular structures [START_REF] Russo | The microscopic pathway to crystallization in supercooled liquids[END_REF][START_REF] Kansal | Diversity of order and densities in jammed hardparticle packings[END_REF].

Order indices of neighbourhood configuration of a central particle 𝑖 are defined on the basis of a vector 𝒒 6 ⃗⃗⃗⃗⃗(𝑖) = [𝑞 ̂6𝑚 (𝑖)] , with 𝑚 = -6, -5,…0,…5, 6. The cosine similarity of a pair of such vectors of neighbouring particles 𝑖 and 𝑗 is calculated as

CosSimi(𝑖, 𝑗) = 𝐑𝐞 [ 𝒒 6 ⃗⃗⃗⃗⃗⃗(𝑖) | 𝒒 6 ⃗⃗⃗⃗⃗⃗(𝑖)| • 𝒒 6 ⃗⃗⃗⃗⃗(𝑗) | 𝒒 6 ⃗⃗⃗⃗⃗⃗(𝑗)| ] = 𝐑𝐞 [ ∑ 𝑞 ̂6𝑚 (𝑖) 6 𝑚=-6 •𝑞 ̂6𝑚 * (𝑗) | 𝒒 6 ⃗⃗⃗⃗⃗⃗⃗(𝑖)||𝒒 6 ⃗⃗⃗⃗⃗(𝑗)| ]. (6)
CosSimi(𝑖, 𝑗) between particles 𝑖 and 𝑗 is positively correlated with the similarity of their individual neighbourhood configurations, and a pair of connected particles has CosSimi(𝑖, 𝑗) ≥ 0.7 [START_REF] Ten Wolde | Simulation of homogeneous crystal nucleation close to coexistence[END_REF]. Hence, the order of individual particles can be positively represented by the parameter [START_REF] Russo | The microscopic pathway to crystallization in supercooled liquids[END_REF],

6 𝑆 6 𝑖 = ∑ CosSimi(𝑖, 𝑗) 𝑛 b 𝑖 𝑗 , (7)
ranging between 0 and 12. For each particle, high 𝑆 6 𝑖 means that particles in an ensemble have similar neighbourhood configurations, characterising local crystalline perfection. Additionally, to describe the overall structural order of an assembly of particles, we make use of [START_REF] Goodrich | Solids between the mechanical extremes of order and disorder[END_REF],

𝐹 6 = 1 𝑁 P ∑ 𝑓 6 (𝑖) 𝑁 P 𝑖=1 , (8)
where

𝑓 6 (𝑖) = 1 𝑛 b 𝑖 ∑ 𝛩 𝑛 b 𝑖 𝑗=1 [CosSimi(𝑖, 𝑗) -0.7] , (9)
with the step function 𝛩(⋅). Spanning the range 0 to 1, the term 𝐹 6 is used to characterise the structural evolution during vibration, as an alternative to the conventional packing fraction.

A bottom-up Coarse-graining methodology is applied in this work to convert the discrete data set of the order index 𝑆 6 𝑖 into a higher-scale continuum form to reveal the spatial distribution of order in the structure. In general, the coarse-graining approach used here takes a set of discrete points 𝑷 𝑖 = (𝑥 𝑖 , 𝑦 𝑖 , 𝑧 𝑖 ) and their corresponding scalar data ℎ 𝑖 as input [START_REF] Goldhirsch | Stress, stress asymmetry and couple stress: from discrete particles to continuous fields[END_REF][START_REF] Weinhart | From discrete particles to continuum fields near a boundary[END_REF]. Instead of describing the density at any point 𝑷 = (𝑥, 𝑦, 𝑧) by 𝜌 dis (𝑷) = ∑ ℎ 𝑖 𝛿(𝑷 -𝑷 𝑖 )

𝑖

where 𝛿(∆𝒑) is the Dirac delta function, the coarse-graining approach transforms this discrete field into a continuous one by replacing 𝛿(∆𝒑) with a positive semi-definite function 𝜑(∆𝒑) . This function fulfils the requirement that the integral of the continuous density function

𝜌 con (𝑷) = ∑ ℎ 𝑖 𝜑(𝑷 -𝑷 𝑖 ) 𝑖 , ( 10
)
is equal to the sum of ℎ 𝑖 for a given volume. The exact form of 𝜑(∆𝒑) is not determinative but the width 𝑤 at which 𝜑(∆𝒑) vanishes holds significance [START_REF] Goldhirsch | Stress, stress asymmetry and couple stress: from discrete particles to continuous fields[END_REF]. Here, a three-dimensional Gaussian function is employed,

𝜑(𝑷 -𝑷 𝑖 ) = 𝑒 - 1 2 [( 𝑥-𝑥 𝑖 𝑟 ) 2 +( 𝑦-𝑦 𝑖 𝑟 ) 2 +( 𝑧-𝑧 𝑖 𝑟 ) 2 ] (2𝜋) 3 2
.

(11)

7

This function vanishes at 𝑤 = 3𝑟 where 𝑟 is the radius of particles. This selection of 𝑤 results from the fact that the order index 𝑆 6 𝑖 is based on particle ensembles roughly occupying a spherical region with the radius of 3𝑟.

Dynamic characterisation

Granular-temperature is often used as a term, analogous to thermal energy at the atomic scale, to describe the kinematics of granular matter. While crystallisation and ordering at the atomic scale are studied with respect to thermal conditions, here we examine the local and global ordering of packing with respect to the dynamic status of particles. Granular temperature quantifies the velocity fluctuation of particles in granular matter rather than being a measure of thermal energy [START_REF] Goldhirsch | Introduction to granular temperature[END_REF]. To obtain the granular temperature, granular matter is first discretised into particle ensembles [START_REF] Hsiau | Experimental investigations of granular temperature in vertical vibrated beds[END_REF]. For each ensemble of 𝑛 particles, the average velocity is calculated in 𝑥, 𝑦, and 𝑧 axis as

Finally, the granular temperature is the mean of the three axial granular temperatures GT = 1

3

(GT 𝑥 + GT 𝑦 + GT 𝑧 ). To maintain consistency with the structural order characterisation, the ensembles used to derive the granular temperature are chosen as the 12-neighbour configuration used in the 𝑆 6 calculation. In this way, every particle has its granular temperature as well as its order index 𝑆 6 , and relations between granular temperature and granular crystallisation can be investigated.

Application to experimental observations

The aforementioned order characterisation methods are applied to two previously studied cases of vibrated granular media, both confined by cylindrical boundaries [START_REF] Reimann | X-ray tomography investigations of mono-sized sphere packing structures in cylindrical containers[END_REF]. ExpC (same notation as in [START_REF] Reimann | X-ray tomography investigations of mono-sized sphere packing structures in cylindrical containers[END_REF]) has a slender shape, 30 mm in its cylindrical diameter and 70 mm in its height, and ExpD (same notation as [START_REF] Reimann | X-ray tomography investigations of mono-sized sphere packing structures in cylindrical containers[END_REF]) is relatively flat, 50 mm in diameter and 40 mm in height. The coordinates of their particles were obtained by XCT techniques. One apparent difference between these two 8 granular media is the order of the superficial particles. ExpC shows extensive hexagonal packing on the cylindrical surface and disordered packing on the top surface (Fig. 1a). However, ExpD shows an opposite trend, despite the existence of discontinuous hexagonal packing (Fig. 1b). In order to expose the internal packing of these granular media, the 𝑆 6 𝑖 of individual particles are computed and fed into coarse-graining approach to construct spatial density mappings of 𝑆 6 , as

given by Figure 1c and Figure 1d. Since the order of packing is positively related to 𝑆 6 , the redder region indicates stronger crystallisation. In consistence with the superficial packing of these two granular media, crystallised regions in ExpC and ExpD exhibit radial and axial accumulation, respectively. The redder colour in the mapping of ExpD implies that fewer defects exist in the packing compared with ExpC. This motivates further investigations into the relations between granular crystallisation and geometrical characteristics.

Results

Simulation parameters

To investigate the influence of container boundaries on granular crystallisation, mono-dispersed frictionless spheres were generated in cylindrical containers of different radii and heights. Then,

9
the evolution of the packing structure inside these granular media during continuous sinusoidal vibration was simulated and characterised to study the crystallisation process. The simulations followed a similar scheme. Firstly, 5000 particles of size 𝑑 were randomly dispersed in a finite cylindrical container with a diameter 𝐷 and were allowed to settle under gravity. This state was considered as the initial state before vibration. Sinusoidal vibration function with an amplitude 𝐴 and frequency 𝑓 was introduced by moving the bottom of the cylindrical container in the axial direction along the gravitational direction. At the upper surface of the granular matter, a lid that could move freely in axial direction confined the granular matter in a negligible pressure of a few Pa. This lid is particularly used to flatten the upper surface for structural characterisation, e.g.

packing fraction and Voronoi tessellation, and reproduce experimental conditions [START_REF] Reimann | X-ray tomography investigations of mono-sized sphere packing structures in cylindrical containers[END_REF]. The vibration was applied over 1000 periods (𝑇 = 1/𝑓), followed by sufficient relaxation (1s) to reach static equilibrium. The vibration amplitude was varied to study the influence of the energy input on the crystallisation. Table 1summarises the material and simulation parameters used in the study.

The diameters and the initial heights of the cylindrical volumes were varied to study the geometrical influence on the crystallisation process. 𝐷/𝐻 =30/75 (with 𝐷/𝑑 =13, 𝐻/𝑑 =33) is characteristic for a rather slender container, while the case 𝐷/𝐻=60/19 (with 𝐷/𝑑=26, 𝐻/𝑑=8) represents a rather flat one. The chosen geometrical parameters 𝐷/𝐻=30/75 and 𝐷/𝐻=50/25 are similar to those investigated by XCT [START_REF] Reimann | X-ray tomography investigations of mono-sized sphere packing structures in cylindrical containers[END_REF] for the purpose of reliable comparison between simulation and experiment.

Selective crystallisation by wall effect and epitaxial growth

By employing DEM simulations to track transient states of particles during vibration, the dynamic development of crystalline structures is unveiled, which supplies information difficult to be obtained by XCT. Following the method introduced in section 2.4, 𝑆 6 density mapping is constructed to spatially expose the propagation of the periodic crystalline structures, featured by high 𝑆 6 regimes. In the sequential 𝑆 6 mappings of the granular media experiencing vibration, the dynamic process of crystallisation is reconstructed by the morphology and intensity evolution as shown in Figure 2.

The general and common phenomenon appearing in all cases is that the region of high 𝑆 6 (shown in red colour) in the density mappings expands and intensifies with vibration, and eventually reaches a quasi-equilibrium state in which only minor changes of ordering could be observed (row by row in Figure 2). Larger amplitude accelerates the crystallisation rates and produces larger crystalline regimes. In addition, increasing D results in regions of higher 𝑆 6 value, indicating finer crystalline regimes. This phenomenon can be explained by the structural difference in those crystalline regimes which will be addressed in later sections. Within the scope of current study, the granular crystallisation is initiated automatically once the vibration commences. It could be inferred that highly disordered granular matter is vulnerable to dynamic perturbations like vibration, while the crystallised granular matter at the quasi-equilibrium state manifests the stability of crystalline structure.

Regarding the growth of crystalline regimes, a wall effect is identified in all geometries. According to the 𝑆 6 density mappings at the early periods, reddish regions indicating highly ordered structure always preferentially appear at the bottom walls and cylindrical walls (the second column in Figure 2 (a) and the first column in Figure 2 (b)). This preferential crystallisation can be partly explained by the existence of a partially ordered layer of particles adjacent to the walls before vibration, which can be seen from the first column in Figure 2 (a). These ordered layers emerge during the settlement of granular matter because the particles need to rest in positions with as many contacts as possible to support themselves under frictionless conditions. However, those partially ordered layers can also be developed near the walls during vibration, which highlights the robustness of two-dimensional ordered layers in spite of collisions by other particles. After the similar prior 11 crystallisation at walls, the emerging crystallisation becomes dependent on the container boundary configuration, i.e., 𝐷/𝑑 and 𝐻/𝑑 . Depending on the origin of the crystalline regimes, two crystallisation modes can be differentiated in this period, the cylindrical mode and the bottom mode. The cylindrical mode induces crystalline regimes in a radially inward direction, whereas crystalline regimes grow upwards along the axial direction in the bottom mode. Thus, a competing mechanism between these two modes is introduced for the subsequent crystalline growth stage, which is determined by the H/D ratio. As the H/d decreases and 𝐷/𝑑 increases, crystallisation at the bottom wall becomes the favoured mode. Conversely, the cylindrical mode dominates the crystallisation in geometries with small 𝐷/𝑑 and large 𝐻/𝑑. This trend is responsible for the observed radial and axial accumulation of crystallised regions found in previously reported experimental work, as shown in section 2.4.

Evolutional phase diagrams were constructed by comparing the 𝑆 6 density maps of different granular media at the specific vibration duration. Two operational phases exist in the phase diagram of small amplitude (𝐴=0.1𝑑) while for the large amplitude (𝐴=0.2𝑑) three phases are observed. In the early periods, all granular media exhibit a "dual-mode cooperating" phase where crystallisation progresses on both cylindrical and bottom walls. When the granular media continue to be vibrated, either the cylindrical mode or the bottom mode prevails in crystallisation ("singlemode prevailing" phase), although different crystalline regimes grown in the other modes can still be identified. The exception D40 is always in the dual-mode phase, in which disordered regions act as boundaries and partition different crystalline regimes. The packing structures of the boundary regions are instable and vary with vibration because the mismatch between crystalline regimes grown in two different modes makes it difficult to reach a final state.

Most interestingly, a third operational phase is present for D30 and D60 in the final periods of large amplitude vibration. This phase represents the extreme scenarios where one of the crystallising modes is eliminated, named as "sole-mode dictating" phase. This third operational phase was experimentally observed in vibrated particle geometries with large 𝐷/𝑑 and small 𝐻/𝑑 by one of the authors. The cylindrical wall layer was fairly hexagonally packed for the longest part of the total vibration period. However, this layer eventually became significantly unstable while at the top bed surface well defined hexagonal patterns still appeared, as shown in Figure 1 b.

12

Although this phenomenon was observed in ExpD with smaller cylindrical diameter [START_REF] Reimann | X-ray tomography investigations of mono-sized sphere packing structures in cylindrical containers[END_REF], the final stage D50 case in simulation also nearly shows "sole-mode dictating" phase, which is possibly resulted from the much shorter vibration duration in simulation. Figure 3 shows the variation of 𝐹 6 in the cylindrical wall layer, the bottom layer and the core region separately with extended vibration periods. It is clearly shown that both wall layers possess higher 𝐹 6 than the core region at the initial state, proving the existence of partially ordered regimes. Not only for the initial state, has this comparison underlined the significance of the wall effect during the granular crystallisation. The 𝐹 6 evolution of the two wall layers conforms the phase separation in the evolution diagrams. In dual-mode cooperating phase, the contrast between the patterns is small and both parts exhibit a common trend of increase-to-stabilise behaviour. As the contrast extends or a crossover arises, the granular media move into the single-mode prevailing phase. In regard to the sole-mode dictating phase, it is foreshadowed by the continuous decrease of the 𝐹 6 in one of the two parts. The decrease of the 𝐹 6 proves those wall layers become disordered by vibration, wiping out the corresponding mode.

Since very large granular systems (particle numbers ≫10 4 [START_REF] Rietz | Nucleation in Sheared Granular Matter[END_REF][START_REF] Aste | Geometrical structure of disordered sphere packings[END_REF]) are beyond the aims of the current study, we conclude that granular crystallisation in containers of finite size is initiated by the wall effect and progresses in modes that depend on the wall geometry. The wall effect of both cylindrical and bottom modes generates quasi two-dimensionally ordered layers with as many contacts as possible which we will subsequently prove to be hexagonal packings. The following processes are the expansion of the hexagonal packing across the wall plane and the epitaxial growth of an adjacent layer. The repetition of these stepwise processes establishes the granular crystallisation in these confined systems. Because of imperfections and the competition with other crystallised layers, each epitaxial layer is smaller than the preceding one, resulting in conical crystalline regimes. From an energetic perspective, the underlying microscale mechanism is proposed as the selection of positions with high number of contacts to enhance the propagation of kinetic energy. Owing to the maximised contacts, the kinetic energy induced by the vibration on the one hand dissipates quickly from one particle to its surroundings in such structures, leaving the particle in a less perturbed state and making the crystalline regimes robust; on the other hand, this kinetic energy is efficiently transferred throughout the crystalline arrangements, triggering the relocation of particles from the disordered regions into particular crystallising positions nearby.

Such steps build up a positive feedback loop to promote crystallisation. Therefore, larger amplitude activates more particles near ordered layers and enhances the possibility for these particles to lodge in an energetically favourable position. However, too much energy will reverse those steps, cause drifting particles and deteriorate the crystallisation. Self-nucleation [START_REF] Rietz | Nucleation in Sheared Granular Matter[END_REF] is not the prevailing mechanism due to the repulsive nature of granular matter while the merging of crystalline regimes is commonly observed.

Transient state vs. Static state

In the context of crystallisation in various systems, the interplay between densification and order formation has been debated by researchers to determine whether the densification precedes the order formation or the other way around [START_REF] Tanaka | Bond orientational order in liquids: Towards a unified description of water-like anomalies, liquid-liquid transition, glass transition, and crystallization[END_REF], while some reports suggest that these two processes occurred simultaneously [START_REF] Berryman | The early crystal nucleation process in hard spheres shows synchronised ordering and densification[END_REF]. Densification in granular matter, characterised by static packing fraction (𝛾) in shearing and vibrating-relaxing processes, indicates different stages of granular crystallisation [START_REF] Dshinde | Shaking-induced crystallization of dense sphere packings[END_REF][START_REF] Aste | Geometrical structure of disordered sphere packings[END_REF], presenting cooperative development between these two aspects. However, 15 according to the temporal evolution of 𝛾 and 𝐹 6 in the right column in Figure 4, 𝛾 at transient state exhibits fluctuation due to the granular dilation resulted from repulsive interaction, hiding its relation with the order formation. On the other hand, the order formation, given by the increase of 𝐹 6 , persists in spite of the unstable densification. Such decoupling phenomenon stresses that the order formation should be focused when the granular crystallisation is studied in transient state.

16

Within the vibration duration in this study two stages can be distinguished according to the evolution of 𝐹 6 . The first stage is characterised by the monotonic and rapid increase of 𝐹 6 , which corresponds to particles in disorder moving into an order positions. The second stage can be treated as a quasi-equilibrium state, featured by the mild change of 𝐹 6 . Particle arrangement of elevated 𝐹 6 leads to symmetrical collisions between particles and uniform dissipation of energy in the granular media, maintaining the systems' stability. In other words, this phenomenon suggests that the order formation cannot fully be accomplished by single vibration mode and wall effect must be carefully examined when producing single-crystal granular matter [START_REF] Carvente | Crystallization of Confined Non-Brownian Spheres by Vibrational Annealing[END_REF][START_REF] Yu | Self-Assembly of Particles for Densest Packing by Mechanical Vibration[END_REF]. These are the reasons that the current simulation cannot reach fully crystallised state (𝛾 ≥ 0.7) when compared with other simulation with periodic boundary condition as well as additional lateral perturbation [START_REF] Dshinde | Shaking-induced crystallization of dense sphere packings[END_REF][START_REF] An | DEM study of crystallization of monosized spheres under mechanical vibrations[END_REF].

The static final states of the granular media are shown in the left columns of Figure 4. The distributions of the Voronoi cell packing fractions are presented with the legends indicating overall 𝛾 and 𝐹 6 . Both 𝛾 and 𝐹 6 increase in the final relaxed states after vibration. The larger amplitude leads to a bigger increment in 𝛾 and 𝐹 6 except for the case of D60. This exception can be attributed to the relatively strong fluidisation because of the least gravitational potential to overcome, i.e., the least 𝐻. These results are in agreement with previously reported DEM studies [START_REF] An | Effect of vibration condition and inter-particle frictions on the packing of uniform spheres[END_REF] and Monte Carlo simulation [START_REF] Dshinde | Shaking-induced crystallization of dense sphere packings[END_REF] as well as experiments [START_REF] Nowak | Density fluctuations in vibrated granular materials[END_REF][START_REF] Reimann | X-ray tomography investigations of mono-sized sphere packing structures in cylindrical containers[END_REF] in which the packing fraction is maximised at an intermediate vibration intensity, which can also be concluded the maximum 𝐹 6 can be introduced by an intermediate vibration intensity. In order to reveal the change in structure between the initial and final states, the Voronoi cell packing fraction was calculated [START_REF] Rycroft | Analysis of granular flow in a pebble-bed nuclear reactor[END_REF]. The distributions are of bell-shape for all the initial states and are characterised by several peaks at the final states.

With increasing 𝐷 and decreasing 𝐻, a peak at ≈ 0.74 develops, which proves the emergence of a highly crystallised structure. It should be noted that the packing fraction shown in Figure 4 is calculated according to the maximum occupied volume of the corresponding granular matter.

Packing fraction 𝛾 in excluding the surface particles of the corresponding granular matter has been given, which is consistent with the proved concept that granular crystallisation occurs with packing fraction larger than RCP (𝛾 ≈ 0.64).

Granular temperature and crystallisation

Granular crystallisation can be regarded as the development of structural order. Analogous to structural order, as the uniformity of the velocities of particles in an ensemble increases, we consider it as the development of dynamic order. In this way, the granular temperature is able to quantitatively describe this dynamic order. In the temporal evolution of the 𝑆 6 distributions, a general trend is evident in which a singular peak gradually intensifies during vibration, shown in Fig. 4(a). Therefore, we make use of the 𝑆 6 distributions to examine how the granular temperature evolves under different conditions. By using the particle connection criteria based on the magnitude of 𝑆 6 , the particles are categorised into two groups, the solid-like and liquid-like. Within the solid-like group, the particles are categorised again into two sub-groups, the highly crystallised and moderately ordered, by setting a threshold for highly crystallised particles, i.e., 𝑆 6 >10; because no liquid-like particle possesses 𝑆 6 which exceeds 10. After such categorisation, the average granular temperature can be extracted for three groups. As discussed in the previous sections, the particles in the boundary regions are supposed to be more active than those in the crystallised regimes. Thus, the dynamic order should be suppressed in the boundary regions of low order. A general observation is that the highly crystallised particles have the lowest average granular temperature but the liquid-like particles have the highest average granular temperature. The distinction of D60 vibrated via large oscillation shown in Fig. 4 demonstrates that the dynamic order is suppressed in the over-fluidised state. Thus, we argue that the granular temperature has a transition from divergence to convergence depending on the geometry and the vibration intensity, and further research is required to fully demonstrate it.

This phenomenon is consistent with the discussion on the relation between the propagation of kinetic energy and the order formation in the crystallisation presented in the previous sections.

Along with 𝑆 6 rise, the order formation and the particles crystallisation, the local structure of individual particles in granular media becomes periodic and symmetric. As a result, the particle collisions inside such crystalline structures turn out to be counterpoised, causing the structures to attain the ability to suppress velocities deviating from the vibration axis and preventing the particles from over-acceleration. In conclusion, the structural order is the foundation for the dynamic order. The granular crystallisation can stabilise the granular media subjected to agitation by resetting the elevation of the granular temperature.

Structural characterisation and evolution

It was demonstrated in the previous sections that the wall effect and the dual crystallisation modes are the major events during vibration. However, 𝑆 6 can capture the degree of order of the particles but cannot distinguish the structural types of crystalline regimes. Thus, the local bond orientation parameters 𝑄 6 local , 𝑊 4 local and 𝑊 6 local are utilised to differentiate the structure of those regimes.

Firstly, 𝑊 6 local is used to examine whether the cubic-based structures exists based on its sign (+/-) 19 contrast, positive for BCC and negative for HCP and FCC. Nearly all particles in these granular media have negative 𝑊 6 local , denoting that the majority of particles sit in hexagon-based structures (not shown explicitly in the figures). This result is consistent with the crystallised structure in granular matter obtained in other research [START_REF] Dshinde | Shaking-induced crystallization of dense sphere packings[END_REF][START_REF] An | DEM study of crystallization of monosized spheres under mechanical vibrations[END_REF][START_REF] Hanifpour | Structural and mechanical features of the order-disorder transition in experimental hardsphere packings[END_REF], implying that BCC is not the preferential structure in granular crystallisation. Next, the (𝑄 6 local , 𝑊 4 local ) coordinates, exhibiting specific coordinates for HCP (0.485, 0.134) and FCC (0.575, -0.159), are employed to characterise the particular structures in the granular media.

In Figure 6, two series of the probability density distributions of (𝑄 or FCC structures, calculated from the neighbour configuration lacked half space due to the boundary conditions. In accordance with the previous discussion, this path reveals the tendency towards hexagonal units in the wall layer during vibration. More importantly, its precedence over the other paths proves the priority of the wall effect. When the vibration continues, HCP and FCC become the main structures emerging during crystallisation. As stated before, the dominating crystallisation mode in the granular media of D30 and D60 is the cylindrical mode and the bottom mode, respectively. According to the density contrast in the HCP and FCC paths between D30 and D60, the cylindrical mode shows a preferential selection of the HCP structure but the bottom mode has little bias on the two paths. The histograms according to probability density function of 𝑄 6 local of the granular matter at the corresponding snapshots are presented in the rightmost column in Figure 6 to give more clearly quantitative sense. Differently from most studies implementing periodic boundary conditions [START_REF] Dshinde | Shaking-induced crystallization of dense sphere packings[END_REF][START_REF] An | DEM study of crystallization of monosized spheres under mechanical vibrations[END_REF][START_REF] Hanifpour | Structural and mechanical features of the order-disorder transition in experimental hardsphere packings[END_REF] or excluding near wall regions [START_REF] Panaitescu | Nucleation and Crystal Growth in Sheared Granular Sphere Packings[END_REF], the priority of the growth of wall layer is identified here. Furthermore, the dominance of FCC over HCP is not seen within the vibration cycles in this work as result of wall effects.

21

First, the formation of HCP structures due to the cylindrical mode contributes most in D30 and D40 cases. The reason for the incapability of producing FCC structure in the cylindrical mode of large wall curvature is elusive. The distinctive crystalline regime in the cylindrical mode is characterised by a deviation from the perfect HCP coordinates introduced by the structural distortion. This cylindrical curvature effect, discussed also in [START_REF] Reimann | X-ray tomography investigations of mono-sized sphere packing structures in cylindrical containers[END_REF], changes the structure greatly from the perfect HCP structure, causing a decrease of the Voronoi cell packing fraction and the coordination number. Nevertheless, as demonstrated in Figure 7, the detailed examination of the 12-neighbour configuration of the distorted particles reveals that the geometrical symmetry is partly maintained in a HCP fashion. HCP segments (blue particles in Figure 7) and rupture particles (yellow particles in

22

One possible angle to further explain this phenomenon can be formed at the stronger mechanical stability of FCC than HCP [START_REF] Heitkam | Packing Spheres Tightly: Influence of Mechanical Stability on Close-Packed Sphere Structures[END_REF], implying in higher resistance to deform in FCC. On the other hand, lower resistance should bring about potential adequateness in HCP to arrange particles in a distorted field like the cylindrical wall without a complete loss of the structural characteristics.

Thus, this stronger deformability may contribute to the prevalence of HCP in the cylindrical mode.

Further, as 𝐷 increases, the distortion due to the wall curvature becomes smaller, resulting in more FCC particle appearing in the cylindrical wall region as shown in Figure 9. Second, although the difference in mechanical stability helps to explain the HCP dominance in the cylindrical mode, the lack of a bias between HCP and FCC in the bottom mode should be discussed because of the contradictory against the argument above [START_REF] Heitkam | Packing Spheres Tightly: Influence of Mechanical Stability on Close-Packed Sphere Structures[END_REF]. According the current result, the mechanical stability plays an important role in distortion effect, while in the distortion free condition the packing history should be considered [START_REF] Heitkam | Packing Spheres Tightly: Influence of Mechanical Stability on Close-Packed Sphere Structures[END_REF], i.e. the vibration protocols used in current work. To elaborate this point, different vibration modes were applied. Using amplitude 𝐴=0.1𝑑

and frequency 𝑓 =50Hz as a reference mode, increasing amplitude ( 𝐴 =0.2 𝑑 ) and frequency (𝑓=70Hz) populates FCC particles, except for D60 due to the aforementioned strong convection by large amplitude. This implies that high energy input benefits the growth of FCC in the bottom mode. Besides, in simulated colloidal systems [START_REF] Fujine | Effect of container shape and walls on solidification of Brownian particles in a narrow system[END_REF][START_REF] Arai | Surface-assisted single-crystal formation of charged colloids[END_REF], the less prevalence of HCP rather than FCC can also be identified when crystallisation initiated from the walls. The identical two-dimension hexagonal packed layers (111) at the walls are considered as the main reason. The energetic favour 23 of FCC may not be hold in these relatively two-dimensional scenarios, which requires further research to clarify the mechanical as well as energetic stability these near wall regions. Apart from the wall effects, this less preference of FCC can partially be attributed to the global packing fraction of the granular matter (𝛾 in ≤0.67) simulated in this work, which is smaller than 24 that packing factor (>0.68) showing stronger preference of FCC [START_REF] Dshinde | Shaking-induced crystallization of dense sphere packings[END_REF]. In addition, recalling the corresponding 𝛾 in of individual pairs of granular matter in Figure 3, higher 𝛾 in also leads more appearance of FCC particles generated from the bottom mode, as shown in Figure 9.

Influence of friction and comparison with experiments

The effect of the boundary geometry on the granular crystallisation has been demonstrated in frictionless media. It is necessary to clarify the influence of friction on crystallisation and compare it to available experimental data. Figure 10 shows for D30 the influence of friction on the evolution of the overall structural index 𝐹 6 during vibration. The initial state was kept identical for both the frictionless and frictional scenarios. The most obvious contrast between these two systems is the efficiency and degree of the granular crystallisation. The frictional energy losses between particles diminish the intensity of the vibration. Therefore, the crystallisation rate is lower in the frictional media and the required vibration time is longer to reach the final crystallised state. The simulated vibration duration (4000 periods for the longest vibration) is significantly shorter than the experimental study (40000 periods for the longest vibration), but the result still matches quantitatively. Using the structural index obtained from the experimental results (𝐹 6 = 0.57808) as a benchmark, orange dash line in Figure 10, the corresponding transient states in the simulation are extracted to reconstruct their temporary morphologies. Compared to the frictionless granular media, the frictional ones exhibit similar patterns of the 𝑆 6 spatial distributions, despite of the radial expansion, the crystalline regimes varying only slightly. Since the 𝑆 6 spatial distributions of different granular media match at the equivalent 𝐹 6 , it is reasonable to argue that the structural evolution in the frictionless granular media represents a fairly complete crystallisation process for a given geometry. Hence, this morphological resemblance that demonstrates that friction scarcely influences the mechanisms of the granular crystallisation and merely hinders the crystallisation growth. Moreover, we have identified the inherent partially ordered regions existing at walls in various initial states conforming to recent XCT data [START_REF] Reimann | X-ray tomography investigations of mono-sized sphere packing structures in cylindrical containers[END_REF]. Nonetheless, such variation in the initial structures is of negligible influence for the wall effects and order development during vibration.

Based on this argument, the experiments can be considered as an intermediate state in the evolution of the frictionless granular media. Due to the frictional force, the granular medium in the experiment is unable to reach the final state of the frictionless simulations, but maintains the same granular crystallisation mechanisms. The distribution obtained in the experiment shows crystalline patterns in the radial direction comparable with the theoretical predictions, while the results about the merging crystalline regime at the bottom of the container show a major difference between theory and experiment. This discrepancy is probably caused by the purely vertical vibration that ineffectively breaks the force chains sustaining the vertical perturbation in the simulations as well as the prolonged vibration in the experiments.

Conclusion

Granular crystallisation has been investigated from the entire system scale down to the individual particle scale in confined granular matter. The results clearly show that vibration naturally brings about a disorder-to-order transition, proving that crystallised structures are maintained during the transient evolution.

The crystallisation process and the role of the wall effect are explained by coarse-graining approach. Internal nucleation growth is restrained due to the purely repulsive interactions, and crystallisation from walls is preferential. The wall effect produces two-dimensional hexagonal packing as a growth template but the following growth of crystalline arrays divides into a cylindrical mode and a bottom mode. In the cylindrical mode, the crystallised structure can be considered as a distorted HCP structure, while in the bottom mode, a mixture of finer HCP and FCC structure is identified. Depending on the geometry, 𝐷/𝑑 and 𝐻/𝐷, competition between these two crystallisation modes falls into three phases during vibration. By increasing 𝐷, the bottom mode crystallisation gradually dominates, leading the crystalline regime to penetrating throughout the entire granular media in the axial direction. In the other case, i.e., when 𝐷 decreases, the crystalline regions growing in the radial direction towards the axis are promoted from the cylindrical mode. Increasing the amplitude of vibration enhances the efficiency of the crystallisation, which can raise the competition level and leads to a sole crystallisation phase, which is commonly seen in experiments. Through the particle scale characterisation we conclude that particle are driven to lodge themselves in structures with as many contacts as possible, because such structures provide sufficient collisions to dissipate kinetic energy and maintain the stability of the granular packing. The relation between granular crystallisation and granular temperature is further explored and it is seen that granular crystallisation leads not only to structural but also to dynamic order.

In this study, granular crystallisation induced by vibration is proved to follow basic processes resulting in a predictable final structure. These results suggest that mechanical, thermal, electrical and other structure related effective properties can be modified by controlling vibration, motivating the continued study of granular crystallisation. Further research can be conducted to explore the

29

physical mechanisms of inducing granular crystallisation. Meanwhile, this study shows the crystallisation can be connected to the statistic description of granular matter. Thus, it would be interesting to seek more precise correlations between granular crystallisation and granular statistics.

4

 4

  Then the velocity fluctuation along 𝑥, 𝑦, and 𝑧 axes are

Figure 1

 1 Figure 1 Reconstruction (left column) and 𝑆 6 density mappings (right column) of two vibrated granular media. ExpCa and c, ExpDb and d.

Figure 2

 2 (a) 𝐴 = 0.1𝑑 (b) 𝐴 = 0.2𝑑 Evolutional phase diagram of the crystallisation in granular media of different height-todiameter ratio. Each 𝑆 6 density mapping consists of central plane slices of three axes, the colour indicated by coarse grained 𝑆 6 suggests disorder in violet direction and crystallisation in reddish direction. In the phase diagrams three phases are marked by background filling, (1) dual-modes cooperating; (2) single-mode prevailing; and (3) sole-mode dictating which are approximately determined by the competition between cylindrical and bottom modes.

13

 13

Figure 3

 3 Figure 3 Evolution of 𝐹 6 for particles groups separated by position. CWfirst layer near the cylindrical wall, BWfirst layer near the bottom wall and Corethe bulk particles.

Figure 4

 4 Figure 4 Overall evolution of each granular medium subjected to vibrations of different amplitude. The histograms (colouring with transparency) in the left column show distributions of the Voronoi cell packing fraction of the initial state and two final relaxed states after vibration (𝐴=0.1𝑑 and 𝐴=0.2𝑑) with legends giving the corresponding overall packing fraction 𝛾 and structural index 𝐹 6 . The corresponding right column plots demonstrate the time variation of the overall at transient states during vibration.

Figure 5

 5 (a) A similar trend in the 𝑆 6 evolution is observed in all cases. The 𝑆 6 accumulates between 10 and 12 in the final state, while the peak right shifts as D increases. (b) Typical granular temperature evolution (D50) in the small amplitude vibration scenario. (c) Typical granular temperature evolution (D40) in the large amplitude vibration scenario. (d) Granular temperature evolution in the relatively strong fluidisation case (D60). (b) and (c) have the same legend shown in (d) with HC -highly crystallised, LL -liquid like, MO -moderately ordered, AVE -averaged.

Figure 6

 6 Figure 6 Smoothed probability density histograms of the crystallised structures appearing in the granular media during vibration for D30 and D60. The (𝑊 4 local , 𝑄 6 local ) coordinates are used to characterise the structure types. The intersects of pairs of dashed lines in orange and green are the coordinates of the FCC and HCP, respectively. The rightmost column gives the smoothed histograms of 𝑄 6 local distributions in the granular matter.

Figure 7 )

 7 are distinguished by the 𝑊 4 local of rupture particles being close to 0. A layer shift occurs in the neighbour configurations of the segment particles but incomplete layers are the typical feature of the rupture particles. Besides, this distortion causes the weaker ordering of the crystalline regimes near the cylindrical walls when compared with the regimes near the bottom walls shown in Figure 2.

Figure 7

 7 Figure 7 Rupture in the HCP structure near the cylindrical wall in D30. (a) Blue particles are distorted HCP particles (0.465 ≤ 𝑄 6 local ≤ 0.505, 𝑊 4 local ≥ 0.08) while yellow particles are rupture particles ( 0.465 ≤ 𝑄 6 local ≤ 0.505 , |𝑊 4 local | ≤ 0.02 ) with size scaled by 0.5 for visibility. (b) Typical rupture section with Particle 2 being the rupture particle. P1, P2, P3 are the 12-neighbour configurations of Particle 1, 2, 3, respectively.

Figure 8

 8 Figure 8 Smoothed histograms of 𝑄 6 local distributions of D40 (left), D50 (middle) and D60 (right) granular matter subjected to different vibration mode. 𝐴=0.1𝑑 indicates 0.1𝑑 amplitude and 50Hz frequency; 𝐴=0.2𝑑 indicates 0.2𝑑 amplitude and 50Hz frequency; and 𝑓=70Hz indicates 0.1𝑑 amplitude and 70Hz frequency.

Figure 9

 9 Figure 9 The corresponding packing structures of individual particles dyed according to the (𝑊 4 local , 𝑄 6 local ) coordinates in red -HCP, blue -FCC, greensurface hexagon and yellow -others. The diameters of the particles are rescaled for visualisation purposes. Left column -𝐴 = 0.1𝑑 and right column -𝐴 = 0.2𝑑.

Figure 10 11 .

 1011 Figure 10 Top -Evolution of the structural index 𝐹 6 of the frictionless and frictional granular media. The top two dashed lines serve as the extension for the final states in the simulation and the middle one labelled with ExpC represents the final state of the experiment performed with a vibration intensity 𝛤 = 2 [36]. Bottom -𝑆 6 spatial distributions of the labelled states in the simulated evolution and the experimental result. Friction coefficient (𝜇) and amplitude (𝐴) values for the simulations are displayed in the legend

Figure 11

 11 Figure 11 Top -𝑆 6 distributions of the final state in the ExpC in [36] and the transient states in the simulations. Particles dyed according to the (𝑊 4 local , 𝑄 6 local ) coordinates are displayed as insets in the 𝑆 6 distribution in red -HCP, blue -FCC and greensurface hexagon. Bottom -The corresponding (𝑊 4 local , 𝑄 6 local ) coordinates distributions. Friction (𝜇), amplitude (𝐴) and duration (𝑡) parameters for the simulations are displayed in the legend.

Figure 12

 12 Figure 12 Top -𝑆 6 distributions of the final state in the experiment ExpD in [36] and the simulations along with the 𝑆 6 spatial distributions and dyed FCC/HCP particles of ExpD in the insets. Bottom -The corresponding (𝑊 4 local , 𝑄 6 local ) coordinates distributions and the comparison of the smoothed histograms of 𝑄 6 local distributions. Friction (𝜇), amplitude (𝐴) and duration (𝑡) values for the simulations are displayed in the legend.

  Dai, W., Reimann, J., Hanaor, D., Ferrero, C., & Gan, Y. (2019). Modes of wall induced granular crystallisation in vibrational packing. Granular Matter, 21(2), 26.

Table 1

 1 Parameters for DEM simulations

	Young's modulus, 𝐸	63 (GPa)
	Poisson ratio, 𝜈	0.2 (-)
	Density, 𝜌	2230 (kg/m3)
	Friction coefficient, 𝜇	0.0, 0.2, 0.5 (-)
	Coefficient of restitution	0.6 (-)
	Diameter of sphere, 𝑑	2.3 (mm)
	Container Diameter/Height, 𝐷/𝐻 30/75, 40/60, 50/25, 60/19 (mm/mm)
		denoted as D30, D40, D50, D60
	Vibration amplitude, 𝐴	0.23 (0.1𝑑), 0.46 (0.2𝑑) (mm)
	Vibration frequency, 𝑓	50 (Hz)
	Gravitational acceleration, g	9.8 (m/s2)
	Vibration intensity, 𝛤	2.3, 4.6 (-)

  6 local , 𝑊 4 local ) for D30 and D60 demonstrate the typical development pathways of the crystalline structures. The similar broad distribution of the pairs of (𝑄 6 local , 𝑊 4 local ) indicates the disordered nature of those granular media at the initial state. Three structuring paths are clearly identified in the granular crystallisation. Once the vibration starts, the path leading to a non-typical structure coordinates, neither HCP nor FCC, appears first. Such coordinates represents the hexagonally packed surface particles in finite HCP

crystallisation in vibrational packing. Granular Matter, 21(2), 26.


Loïc Joud
email: loic.joud@danielson-eng.fr

Alan Keromnes
email: alan.keromnes@u-bourgogne.fr

Daniela Chrenko
email: daniella.chrenko@utbm.fr

Luis Le Moyne
email: luis.le-moyne@u-bourgogne.fr

Luis Le

Rui Da

Silva Lourenco

Prédiction de cycle et stratégie de gestion énergétique basée sur une méthode d'optimisation globale

Keywords: stratégie de gestion énergétique, prédiction, méthode d'optimisation globale, style de conduite energy management strategy, prediction, global optimisation method, driving behaviour Résumé

La nécessité de prédire le cycle de conduite d'un automobiliste afin d'améliorer sa gestion énergétique et ainsi réduire la consommation de carburant a été démontré (Debert, Quigley). Dès lors, il serait aisé d'appliquer une méthode d'optimisation globale, telle que la Programmation Dynamique (DP en anglais), pour obtenir la stratégie optimale. Cette méthode est déjà utilisée dans de nombreuses applications (Sciarretta), mais au quotidien nombreux sont les facteurs environnants qui peuvent perturber les conditions de circulation (météo, travaux, accidents…. Les premiers résultats ont permis de mettre en évidence la répétabilité calendaire (et horaire) des trajets. L'objectif de cette étude est aussi la prédiction d'un cycle représentatif du style de conduite de l'automobiliste. Le cycle moyen n'étant pas envisageable puisqu'il tend à supprimer toute la dynamique de conduite, plusieurs méthodes de reconstruction d'un cycle ont été envisagées. Il a été choisi de retenir parmi la base de données un seul cycle pour chaque trajet différent qui soit représentatif de toutes les occurrences sur ce trajet, en utilisant une méthode de corrélation croisée (figure 1). A chaque nouveau trajet, le cycle représentatif est redéfini, ce qui permet d'effacer de la mémoire les cycles de vitesse perturbés par des conditions de trafic exceptionnelles.

Pour étudier l'impact qu'a la non adaptabilité en temps réel d'une stratégie de gestion énergétique basée sur une méthode d'optimisation globale (DP), celle-ci est préalablement appliquée au cycle le plus représentatif choisi comme décrit précédemment. La DP fournit la solution optimale de répartition énergétique pour un cycle connu. Cette répartition est alors appliquée aux autres cycles de conduite représentant le même trajet, et comparé à la DP de chaque cycle. On s'aperçoit alors qu'il est nécessaire de mettre en place une stratégie de gestion énergétique en temps réel qui s'appuie sur les résultats préliminaires d'une DP, mais qui est capable de s'adapter à l'évolution des conditions de circulation.

  ) et impacter le cycle de conduite. Il est donc nécessaire d'élaborer une stratégie adaptative en temps réel. Les travaux qui suivent tendent à démontrer et à quantifier l'intérêt d'une telle stratégie. Dans cet article, le véhicule considéré est la Noao, véhicule construit par le pôle de performance de Nevers Magny-Cours (PPNMC). Le groupe moto-propulseur de ce véhicule est un hybride série. La différence de consommation entre la stratégie présente dans la voiture et de type basée sur des règles (RB en anglais), et la DP est présentée sur la figure 1. Le calcul est effectué sur un cycle représentant 4 tours du circuit de Magny-Cours.

Figure 1 :

 1 Figure 1: Evolution de la consommation et de l'état de charge de la batterie sur un même cycle pour le RB et la DP Dans le cadre d'une étude expérimentale basée sur le volontariat, des données de conduite ont été relevées chez plusieurs participants, en respectant la méthode des quotas (âge, sexe). La population

Figure 2 :

 2 Figure 2 : Exemple de trajets identiques d'un volontaire de la campagne d'essai, et choix du trajet le plus représentatif


Amel Benmouna
email: amel.benmouna@utbm.fr

Mohamed Becherif

Jian Chen

H Eulerchen

Frédéric Gustin

Daniel Depernet

Interconnection and damping assignment passivity based control for cell and battery vehicle : simulation and experimentation

Keywords: Energy management, hybrid system, fuel cell, battery, control based on passivity

come   L'archive ouverte pluridisciplinaire

I. INTRODUCTION

The automobile transport is actually the sector responsible of the emission of greenhouse gases that are mainly responsible at their turn of climate warming. The hybrid electric vehicle represents a good alternative candidate to replace the internal combustion engine [START_REF] Zeng | A Parallel Hybrid Electric Vehicle Energy Management Strategy Using Stochastic Model Predictive Control With Road Grade Preview[END_REF]. This system offers zero emission technology [START_REF] Benmouna | Fault diagnosis methods for Proton Exchange Membrane Fuel Cell system[END_REF]. However, the control of the energy distribution between the embedded sources is an important task [START_REF] Xia | Power management strategy of hybrid electric vehicles based on quadratic performance index[END_REF] [4] [START_REF] Kraa | Experimental validation of a dual loop control of two phases interleaved boost converter for fuel cell applications[END_REF]. This paper focalises on the control based on passivity approach of the hybrid system composed of FC as primary source and battery as the second source. This paper is organized as follows: the section 2 presents the description of the Interconnection and Damping Assignment-Passivity-Based Control (IDA-PBC) approach. The section 3 explains the hybrid structure and the modeling of hybrid system where the state space equations are given. In this section, the passivity control principle based on the IDA-PBC is also developed. In the section 4, the proposed scenario is given considering the role of different used sources. In the section 5, the obtained results from the simulation is presented and discussed. The experimental part is given in the section 6, where the test bench is presented and experimental results are shown and argued. The section 7 is dedicated to the conclusion.

II. PRINCIPAL OF IDA-PBC APPROACH STUDY

Equation of nonlinear system affine in control is given by:

(1)

Where is the state space vector, is the output vector, f(x), g(x) and h(x) are locally Lipschitz functions and is the control input vector. A Port Control Hamiltonian (PCH) form of the equation system can be written by:

(2) J(x) is an n×n skew symmetric matrix, R is a n×n positive semidefinite symmetric matrix and ∇H(x) is the gradient vector of the energy function H(x) of the system .

(3) is a full-rank left annihilator of g(x), that is . PCH system, with H(x) non-negative, is passive system. A recent and very interesting approach to solve these problems is the IDA-PBC approach which consists on identifying the natural energy function of the system called H(x), and then rewrite the nonlinear system versus the gradient of the energy function [18] [START_REF] Becherif | Passivity-Based Control of Hybrid Sources[END_REF].

ISSH2 ' 2017

*Corresponding author email: amel.benmouna@utbm.fr After this step, the passivity command called IDA (Interconnection and Damping Assignment) is used. Consider the system (1), assumes that matrices , are the desired matrices, and function : that verify the partial differential equation. [START_REF] Kraa | Energy Management of Fuel Cell/ Supercapacitor Hybrid Source Based on linear and sliding mode control[END_REF] is such that [START_REF] Kraa | Experimental validation of a dual loop control of two phases interleaved boost converter for fuel cell applications[END_REF] Where is a desired equilibrium to be stabilized, taking the control u as:

(6)
Then the closed-loop dynamic takes the form:

(7)
is then a stable equilibrium. It will be asymptotically stable if in addition, is an isolated minimum of (x) and the largest invariant set under the closed-loop dynamics (7) contained in

III. HYBRID FC SYSTEM STRUCTURE

In the following, the model is written using the new state space , which is defined as the error between the state and its equilibrium value . The desired closed loop energy function is given by: (12)

The error dynamic equations are expressed using the gradient of the desired energy as follows:

The desired equilibrium are noted as following: The control objective is to control the battery current to its reference. In this work, a power source is used to emulate the FC role.

V. SIMULATION RESULT

VI. EXPERIMENTAL RESULT

The test bench is developed at State Key Laboratory of Industrial Control Technology, Zhejiang University (Hangzhou, China). The experimental is illustrated in the Fig. 3 and the objective is to valid the correct operation of the chosen control strategy, study the role of used sources and the distribution control of the power between different elements in real time.

Fig. 3. Powers of different HEV source

The used test bench in this study has been built by authors (see Fig. 3). This test bench makes it possible to emulate energy exchanges within a vehicle on a reduced scale. Comparing to the classical regulators, the used strategy is suitable for nonlinear systems, i.e. no needs to linearize around operating point. And the stability is over all real set. The test bench is composed of:

 Power source emulating a FC role.

 Battery stack consists of 6 cells of 3.2 V connected in series (36Ah each).  An electrical charge emulating the power demand.  Two DC / DC converters for the power source and the battery.  Ammeters to measure the current  Voltmeters to measure the voltage.  LabVIEW software with DAQ MAX Acquisition system.

Fig. 4. Experimental results of powers using IDA-PBC technique.

VII. CONCLUSION

A studied hybrid system is composed of a power source emulating the FC as primary source and battery. A control strategy along with the energy management is developed using the passivity control. For this purpose, this paper introduces the system mathematical modeling. The IDA-PBC is used to design a smart nonlinear energy management control. The experimental and simulation results show the correct and the efficacy of the proposed control. This is clearly demonstrated in the obtained curve through the simulation and practical results. The obtained results reflect exactly the proposed scenario. The other results will be given in an extended journal version.

Fig. 1 :

 1 Fig.1: hybrid system structure

  are the control of the FC boost converter, and the battery control of the buck-boost converter, respectively. The natural energy function is:(9)With matrix Q= diag(LFC, CDC, LB, LL). If, the controls law at the equilibrium are

*

  Corresponding author email: amel.benmouna@utbm.fr is a skew symmetric matrix defining the interconnection between the variables and is symmetric positive semi definite matrix defining the damping of the system. With control laws which are chosen as:(15)Where r is a positive design parameter.IV. PROPOSED SCENARIOThe major objectives to achieve are: -Voltage control of the load and therefore the voltage of DC bus through the command of the two DC-DC converters.-Ensure that the battery provides the power required during the transient phases.-Supply the load by the battery when the FC is not able, for example in the system start up or under limitation. -The energy management strategy is based on the control of the battery current to a desired battery current reference. -This reference can be imposed regarding different criteria (battery SOC, FC state of heath, hydrogen level, start-up, load transient power,…).

Fig. 2 .

 2 Fig. 2. Powers of different HEV source.

ACKNOWLEDGMENTS

The authors gratefully acknowledge the College of Control Science and Engineering, Zhejiang University (Hangzhou, China), where the practical work has been carried.