# Dynamic Self-assembly of Non-Brownian Spheres

J. Salazar, J Marc Simon, J Carlos Ruiz-Suárez, Francisco Peñuñuri, Osvaldo Carvente

► To cite this version:

J. Salazar, J Marc Simon, J Carlos Ruiz-Suárez, Francisco Peñuñuri, Osvaldo Carvente. Dynamic Self-assembly of Non-Brownian Spheres. International, 2017, Montpellier, France. hal-02368512

HAL Id: hal-02368512

https://hal.science/hal-02368512

Submitted on 18 Nov 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Dynamic Self-assembly of Non-Brownian Spheres

J. Marcos Salazar<fnref n="1" />,*, J. Marc Simon<fnref n="1" />, J. Carlos Ruiz-Suárez<fnref n="2" />, Francisco Peñuñuri<fnref n="3" /> and Osvaldo Carvente<fnref n="3" />

**Abstract.** Granular self-assembly of confined non-Brownian spheres under gravity is studied by Molecular Dynamics simulations. Starting from a disordered phase, dry or cohesive spheres organize, by vibrational annealing into BCT or FCC structures, respectively. During the self-assembling process, isothermal and isodense points are observed. The existence of such points indicates that both granular temperature and packing fraction undergo an inversion process. Around the isothermal point, a sudden growth of beads having the maximum coordination number takes place. We show by a density fluctuation analysis that a transition form a disordered phase to a crystalline structure may be associated to a first-order transition.

## 1 Introduction

Spontaneous pattern formation by self-assembly is an appealing research subject given its potentiality as a manufacturing process. Self-assembly processes are common through nature and technology. The formation of regular cristalline lattices is a fundamental process in self-assembly and is a method for converting ≈ 100nm structures into photonic materials [1]. In natural systems crystal formation is widely present, e.g., in the geological production of opals and in the biological formation of proteins and virus [2]. Some years ago the methodology of self-assembly has been used for the fabrication of electronic devices and appears as a practical strategy for manufacturing nanostructure. Beyond the technological interest of self assembly it has always fascinated human curiosity for its æsthetic aspects. But, the physical processes governing self-assembly is an open question and the present manuscript address some issues that may be of interest.

Several years ago we serendipitously discovered that a collection of dry spheres subjected to vibrational annealing, forming a body centered tetragonal (BCT) lattice in a commensurable container, self-assembly instead, upon the addition of a small drop of oil, into a face centered cubic (FCC) [3]. We speculated that the observed polymorphism was due to cohesion introduced by capillary forces. In the present article, we perform Molecular Dynamics simulations (MD) with the aim to confirm our earlier speculations. Intriguingly, we found that the interrelating dynamics of the grains, subjected to a vibration that changes with time, gives rise to temperature and density inversions. This creates two peculiar points easily recognizable as isosbestic points. When the system starts

to cross these points density fluctuations decrease and the number of beads having the maximum coordination number (12 for FCC and 10 for BCT) suddenly appear to finally give place to a granular crystallization. Our density fluctuation analysis shows that this transition may be associated to a first order phase transition.

The grains are indifferent to thermal fluctuations so one must provide some sort of stochastic mobility. Such mobility is simply achieved by shaking the box where the grains are contained. Thermodynamic-like macroscopic quantities can be well defined to describe the collective behavior of these constituents. In these systems, the packing fraction, pressure and granular temperature are also well defined [4–14]. However, due to the fact that granular materials are dissipative and subjected to gravity influence, they exhibit a non-uniform temperature and packing fraction when energy is injected to them. Indeed, in a vibrating box filled by spheres one can clearly observe that both quantities change with height [15] and these gradients persist as we change the frequency of the vibration.

## 2 Molecular Dynamics

MD simulations were performed for tracking the packing evolution of dry and cohesive ensembles composed of mono-dispersed steel beads (radius $r = 1.19 \times 10^{-3}$ m and mass $m = 5.54 \times 10^{-5}$ kg) inside a container subjected to a vertical vibration strength. We used ensembles having 666, 1099 and 6120 beads. Given that the results obtained for 1099 and 6120 are similar to those of 666 beads, we focus mainly on the results for 666 beads. The dimensions of the container ensures (within a small tolerance) either a BCT structure ($L_x = L_y = 2r(1 + (N_{cells} = 4)\sqrt{3}) = 18.87 \times 10^{-3}$ m) or a FCC cohesive structure

<footnote n="1">¹Laboratoire Interdisciplinaire Carnot de Bourgogne, UMR 6303 CNRS Université de Bourgogne
Franche-Comté, 9 Avenue Alain Savary, B.P. 47870, 21078 Dijon Cedex, France</footnote>

<footnote n="2">²CINVESTAV Monterrey, Autopista Nueva al Aeropuerto km 9.5, Apodaca, Nuevo León 66600, México</footnote>

<footnote n="3">³Departamento de Ingeniería Física, Universidad Autónoma de Yucatán, 97310, Mérida, Yucatán, México</footnote>

<footnote>*e-mail: jmarcos@u-bourgogne.fr</footnote>
$$ (L_x = L_y = 2r(1 + (N_{cells} = 5)\sqrt{2}) = 19.3 \times 10^{-3}m). $$
The vibrated ensemble is subjected to an iso-gamma cooling, where the dimensionless acceleration, $\Gamma = A(2\pi v)^2/g$, remains constant while the frequency $v$ (amplitude $A$) is increased (decreased). This iso-gamma annealing drives the system from a dilute state to a compact state [3]. It is important to remark that in the absence of annealing there is no possibility to observe self-assembling (Fig.1).

The interacting force acting between two colliding particles $i$ (mass $m_i$) and $j$ (mass $m_j$) is given by $\vec{F}_{ij} = \vec{F}_N + \vec{F}_T$. The former term is given by: $\vec{F}^n = k_n \xi_{ij} \vec{n} - \frac{\gamma_n}{2} m_{eff} \vec{v}_n$, where $\vec{v}_n = |\vec{v}_j - \vec{v}_i|\vec{n}$, $\vec{n} = \frac{\vec{r}_j - \vec{r}_i}{|\vec{r}_j - \vec{r}_i|}$, $k_n$ is a material constant associated to the hardness of the beads and $\xi_{ij} = R_j + R_i - |\vec{r}_j - \vec{r}_i|$, representing the overlapping between particles ($R_i$, $\vec{r}_i$, $\vec{v}_i$ are respectively, the radius, position and velocity of particle $i$). The second term in $\vec{F}^n$ represents the energy dissipation by the collision and is proportional to the relative velocities of the particles. Where $\gamma_n$ represents the damping coefficient, $m_{eff} = m_i m_j / (m_i + m_j)$ the effective mass.

The tangential forces, $\vec{F}^T$, are modeled by using the tangential relative velocities. For two colliding particles the relative velocity of the spheres at the point of contact is given by: $\vec{v}_{rel}^T = \vec{v}_{rel} - \vec{n}(\vec{n} \cdot \vec{v}_{rel})$, with $\vec{v}_{rel} = \vec{v}_i - \vec{v}_j + R_i \vec{n} \times \vec{w}_i + R_j \vec{n} \times \vec{w}_j$ with $\vec{w}_i, \vec{w}_j$ the angular velocities. We use the force model introduced by Cundall [16], where the static friction is described by a spring acting along the tangential direction of the contact plane. The spring elongation is given by: $\zeta(t) = \int_{k}^{t} v_{rel}^{T}(t') dt'$. This equation allows to determine the restoring tangential force, $|k_t \zeta|$, ($k_t$ represents the spring stiffness) which is limited by the Coulomb's rule. Henceforth, the tangential force can be written as: $F^T = -sign(\vec{v}_{rel}^T) * \min(|k_t \zeta|, \mu|\vec{F}^n|)$ where $\mu$ is the friction coefficient of the beads and $\vec{F}^T = F^T \vec{t}$ with $\vec{t} = \vec{n}^\perp$. In our simulations we consider the restitution coefficient independent of the impact velocity. The parameters used in the simulations are: $k_n(N/m) = 81293.298$, $\gamma_n(s^{-1}) = 0.245$, $\gamma_t(N/m) = 0.122$, $\mu = 0.5$, $\epsilon = 0.9$, $m_{eff}(g) = 2.77 \times 10^{-5}$, $\tau_c(sec) = 5$. The capillary forces of wet particles represented in our simulations are given by:

$$ F_{coh}(s) = \frac{rA\pi\sigma[A(\tan(\beta) - \sec(\beta)) + 1]\tan(\beta)}{A \sec(\beta) - 1}, \quad (1) $$

where $A = (1 + \frac{s}{2\mu})$, $s$ the separation between beads, $\sigma$ the liquid surface tension and $\beta$ the meniscus angle formed by the liquid between two interacting beads. The parameters used here are: $\sigma = 0.03(Nt/m)$, $1 \le s \le 60(\mu m)$ and $\beta = 12.5^\circ$ (see Ref.[17] for more details).

## 3 Temperature and density inversion

The simulation box height is divided in bins, $bin_k$, of thickness $z_{k+1} - z_k$ equal to two beads diameter. The system vibrates at a fixed frequency $\nu_j$ which is active for a period of time (five seconds). In these time intervals, the local granular temperature $\delta v_k(bink, v_j)$ [13], and the local packing fraction $\rho_k(bink, v_j)$ are estimated by bin. Every 6 ms the instantaneous state of the system

is sampled for calculating the mean temperature and packing fraction (Fig.1 and Fig.4). We use the Voronoi tessellation to calculate the volume variation, $V_{voro}(i)$, associated to each particle $i$ [18] for all the frequencies $\nu_j$. Then $\rho_k(bink, v_j) = \langle V_i/V_{voro}(i) \rangle \forall i \in bin_k$, where $V_i$ is the volume of the spheres. The brackets represent average over particles and time.

In Fig.1 we show snapshots of the system for several frequencies of a) dry and c) wet beads, respectively. The colorbar is proportional to the granular temperature (red for hot and blue for cold). From the first figure, we can easily notice that at the beginning of the annealing there is a temperature gradient along the direction of gravity (warmer at the bottom and colder at the top). As the annealing continues, the gradient eventually disappears at a given frequency $\nu_c = 34Hz$ (temperature homogeneity represented by uniform white color). Above this frequency, a temperature inversion appears (the bottom becomes colder than the top). The mean local temperature curves for each $bin_k$, as a function of frequency, are given in Fig. 1(c, d). Accordingly,

$$ \left. \frac{\partial(\delta v_k(bink, v_j))}{\partial bin_k} \right|_{\nu_c} = 0. \quad (2) $$

When the system crosses the region of the isothermal point established by the Eq.2, the maximum coordination number suddenly starts to increase as seen in Fig.2 where we show, for each annealing frequency, the mean number of particles having the maximum first neighbors (10 for BCT and 12 for FCC). While the BCT lattice can form only with the annealing of the system, the FCC structure requires also cohesion. Fig.3 shows how the mean Voronoi volume decreases during the annealing process and we can readily see that near 100 seconds an intriguing change of regime occurs (detailed below). Above this time the mentioned volume per bead remains unchanged.

<figure><img src="image_2.png" /><figcaption>**Figure 2.** The mean number of particles having the maximum coordination number. The maximum coordination number sharply grows around 34-36 Hz, when the temperature vs frequency profiles cross.</figcaption></figure>

Previous studies have shown how the density and the temperature of a vibrated granular system changes with height [15]. For the studied system presented here we find similar temperature profiles. Concerning the local density $\rho_k(bink, v_j)$ as a function of the container's height, Fig.4 shows distinctive crossing points. These points define a
<figure><img src="image_3.png" /><figcaption>**Figure 1.** a) Snapshots of the dry system at 22, 32, 42, and 62 (Hz). The system crystallizes into a BCT structure; b) Granular temperature versus frequency for each *bin*<sub>*k*</sub>. We clearly observe an inversion beyond the isothermal point. Panels c) and d) are for the cohesive system crystallizing into a FCC. The colorbar indicates that particles with the highest and lowest temperature are in red and blue respectively.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Time evolution of mean Voronoï volume of beads during the annealing process for a BCT and a FCC structure. It is interesting to remark that near ninety seconds the Voronoï volume decay shows a change of regime.</figcaption></figure>

small region of constant packing fraction that does not depend on vibration conditions. We call it isodense point. Then, it is possible to establish that: $\frac{\partial(\rho_k(bin_k, v_j))}{\partial v_j}|_{bin_k} = 0$. For the dry (BCT) and cohesive (FCC) structures $k^* = 9$ or (6), respectively.

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Local packing fraction for some frequencies as a function of height (*bin*<sub>*k*</sub>) using Voronoi tessellation. a) BCT crystal and b) FCC crystal. Isodense points appear near the final heights for each lattices.</figcaption></figure>

## 4 Phase transition

For characterizing the transition obtained from a disorder ensemble of steel beads to a crystalline structure (Fig.1) we use concepts from statistical mechanics. An efficient way to track a transition is to follow the time evolution of the thermal compressibility, $\chi_T$. This quantity in the grand canonical (GC) ensemble can be computed from density fluctuations [19] according to the following expression:

$$ \frac{N}{V} k_B T \chi_T = \left( \frac{\langle \delta N^2 \rangle}{N} \right)_{GC} = \Gamma^{-1}(V) \quad (3) $$

where $k_B$ is the Boltzmann constant and $\delta N$ the fluctuations associated to the number of particles $N$ inside a volume $V$. Actually, the thermodynamic factor $\Gamma^{-1}$ gives a measure of the deviation from the ideal behavior of the system. Given that the system sizes simulated are not infinite, we compute the temporal evolution of the density fluctuations by considering (for each configuration) spherical volumes randomly located within the simulation box and by discarding the container surfaces. These spheres have a radius equal to $7.0 \times 10^{-3}$ cm. Even when this choice
is rather arbitrary, it has to be taken between twice the particles size and half the simulation box size [20]. We tested different sizes from 2.0×10⁻³ to 9.0×10⁻³ and we verified that the results are similar. Fig.5 shows that by starting from a disordered system Γ⁻¹ starts to increase indicating the re-organization of beads into an ordered structure. Afterwards an abrupt discontinuity of Γ⁻¹ appears which may be associated to a first order transition. The behavior Γ⁻¹ observed for ensembles of beads leading to a BCT or FCC structure are differentiated by the time at which the transition occurs. In Fig.5 we can readily notice that the transitions to a BCT and FCC structures are accompanied by the presence of small peak around 1600 and 1700, respectively, which can be seen as precursors of the transition. Our results suggest that the transition occurs by jumps to well differentiated stable states for reaching the final state. This conjecture is under analysis and will be presented elsewhere.

<figure><img src="image_6.png" /><figcaption>**Figure 5.** The Γ⁻¹ time evolution shows a phase transition for a BCT crystal and a FCC crystal around 1650 and 1800 iterations steps (≈ 100 sec), respectively. Notice that the transition is characterized by a grow of Γ⁻¹ and the presence a small peak before the transition takes place.</figcaption></figure>

## Acknowledgements

This work has been supported by CONACYT, Mexico, under grants 168281 and ANR-France project MI2C 2015-2019.

## References

[1] T. D. Clark, J. Tien, D. C. Duffy, K. E. Paul, and G. M. Whitesides J. Am. Chem. Soc., **123**, 7677-7682, (2001).

[2] G. M. Whitesides and B. Grzybowski, Science **295**, 2418 (2002).

[3] O. Carvente and J. C. Ruiz-Suárez, Phys. Rev. Lett. **95**, 018001 (2005).

[4] J. B. Knight, C. G. Fandrich, C. N. Lau, H. M. Jaeger, and S. R. Nagel, Phys. Rev. E. **51**, 3957 (1995).

[5] E. R. Nowak, J. B. Knight, E. Ben-Naim, H. M. Jaeger and S. R. Nagel, Phys. Rev. E. **57**, 1971 (1998).

[6] J. S. Olafsen and J. S. Urbach, Phys. Rev. E. **60** R2468 (1999).

[7] S. Warr and J. P. Hansen, Europhys. Lett. **36**, 589 (1996).

[8] F. Rouyer and N. Menon, Phys. Rev. Lett. **85**, 3676 (2000).

[9] C. R. K. Windows-Yule and D. J. Parker, Phys. Rev. E. **87**, 022211 (2013).

[10] G. D'Anna, P. Mayor, A. Barrat, V. Loreto and F. Nori, Nature **424**, 909 (2003).

[11] T. Shinbrot, Nature Physics **9**, 263 (2013).

[12] R. Ramírez and. R. Soto, Physica A. **322**, 372 (2003).

[13] I. Goldhirsch, Powder Technology, **182**, 130 (2008).

[14] A. Bodrova, A. K. Dubey, S. Puri and N. Brilliantov, Phys. Rev. Lett. **109** 178001 (2012).

[15] R. D. Wildman, J. M. Huntley, J.-P. Hansen, D. J. Parker, and D. A. Allen, Phys. Rev. E. **62**, 3826 (2000).

[16] P.A. Cundall and O.D.L Strack. *Géotechnique* **29**, No 1, 45-67 (1979).

[17] J. M. Salazar. *Granular Matter*, DOI 10.1007/s10035-014-0497-1 (2014).

[18] C. Briscoe, C. Song, P. Wang and H. A. Makse, Phys. Rev. Lett. **101**, 188001 (2008).

[19] Computer simulation of liquids, M. P. Allen and D. Tildesley, Oxford University Press, New-York (1987).

[20] S. K. Schnell, T. J. H. Vlugt, J.-M. Simon, D. Bedeaux and S. Kjelstrup Chem. Phys. Letters **504**, 199 (2011).

[21] O. Carvente, M. Salazar-Cruz, F. Peñuñuri and J.C. Ruiz-Suarez. Phys. Rev. E. **93**, 020902(R) (2016).