Sedimentation of a suspension of rods: Monte Carlo simulation of a continuous two-dimensional problem

TL;DR

This study uses kinetic Monte Carlo simulations to explore how rods sediment and self-assemble in a two-dimensional suspension, revealing controllable anisotropic electrical properties.

Contribution

It introduces a continuous off-lattice model for rod sedimentation and analyzes how initial density and sedimentation rate influence self-assembly and anisotropy.

Findings

Different sedimentation-driven self-assembly patterns were observed.

Electrical conductivity anisotropy can be finely tuned by initial conditions.

Spatial distribution and order parameters depend on initial density and sedimentation rate.

Abstract

The sedimentation of a two dimensional suspension containing rods was studied by means of kinetic Monte Carlo (MC) simulations. An off-lattice model with continuous positional and orientational degrees of freedom was considered. The initial state before sedimentation was produced using a model of random sequential adsorption. During such sedimentation, the rods undergo translational and rotational Brownian motions. The MC simulations were run at different initial number densities (the numbers of rods per unit area), $\rho_i$, and sedimentation rates, $u$. For sediment films, the spatial distributions of the rods, the order parameters and the electrical conductivities were examined. Different types of sedimentation-driven self-assembly and anisotropy of the electrical conductivity were revealed inside the sediment films. This anisotropy can be finely regulated by changes in the values of $\rho_i$ and $u$.

Figures

Equations

$$\lambda =\frac{k_{B}T}{mg},$$

$$u=\frac{mg}{f_t},$$

$$D_t=\frac{k_{B}T}{f_t},$$

$$\mathrm{Pe}=\frac{r}{\lambda },$$

$$f_r$$

$$f_{\parallel }$$

$$f_{\perp }$$

$$\Delta r_{\perp }/\Delta r_{\parallel }$$

$$l\Delta \theta /\Delta r_{\parallel }$$

$$\Delta t_B=\frac{{\mathcal{A}}_i}{3}\Delta t_{MC},$$

$$t_B=\frac{\Delta t_{MC}}{3}\sum _{i=1}^{t_{MC}}{\mathcal{A}}_i,$$

$$u=mg/f_{\parallel }.$$

$$u_x/u$$

$$u_y/u$$

$$\mathrm{Pe}={\tau }_B/{\tau }_S=u/{\alpha }^2.$$

$$S_f=\frac{1}{N}\sum _{i=1}^{N}2{\cos }^2{\theta }_i-1,$$

$${\rho }_m^{\ast }\propto u^{\beta },$$

Full text

I Introduction

Gravitational or centrifugal techniques are widely used in practice for the separation of suspended particles Schaflinger (1990), industrial filtration Tiller et al. (1995), the clarification of solvents and the flotation of suspended sewage solids Chen and Scott (1997), the production of colloidal phases Chen et al. (2015), and the preparation of electrically conductive and transparent films Petek et al. (2015). An important characteristic of the deposition process is the sedimentation length, $\lambda$, as determined by Perrin in measurements of the Boltzmann constant, $k_{B}$ Perrin (1913). This value is defined as

[TABLE]

where $k_{B}T$ is the thermal energy, $m$ is the buoyant mass of the particle and $g$ is the acceleration. The buoyant mass of the particle is defined as $m=\Delta\rho V$, where $\Delta\rho$ is the difference between the densities of the particle and the liquid, and $V$ is the volume of the particle. In very dilute systems (ideal gaslike systems), the equilibrium sedimentation profile of particle number density, $\rho$, is barometric and decays exponentially with height, $\rho=\rho_{0}\exp(-y/\lambda)$, where $\rho_{0}$ is the density at the base of sediment, $\rho_{0}=\rho(0)$. This approximation works when the value of $\lambda$ is large compared to the inter-particle correlation length, but small compared to the vertical size of the system.

For an isolated sphere of radius $r$ with no-slip boundary conditions on the surface, the sedimentation rate can be evaluated using Stokes law

[TABLE]

where $f_{t}=6\pi\eta r$ is the Stokes’ friction coefficient related to the translation diffusion coefficient as

[TABLE]

and $\eta$ is the viscosity of surrounding liquid Schuck (2016).

The Péclet number is defined as the ratio of the Brownian, $\tau_{B}$, and sedimentation, $\tau_{S}$, times of a particle over the same distance (e.g., a distance of one radius), i.e. $\mathrm{Pe}=\tau_{B}/\tau_{S}$ Royall et al. (2007). Using the values $\tau_{B}=r^{2}/D_{t}$ and $\tau_{S}=r/u$, we can obtain $\mathrm{Pe}=ru/D_{t}$. Accounting Eqs. (1), (2), and (3),

[TABLE]

i.e., the Péclet number simply presents the ratio of particle size $r$ and sedimentation length $\lambda$.

Uniformity in particle arrangement is expected when the effects of diffusion are dominant ($\mathrm{Pe}\ll 1$), whereas in the opposite case, when the effects of sedimentation are dominant ($\mathrm{Pe}\gg 1$) a spatial gradient in the distribution of particles is typically observed. For example, at room temperature, $T=298$ K, $\Delta\rho=10^{3}$ kg/m3, and $g=g_{0}=9.8$ m/s2 ($g_{0}$ is the acceleration due to gravity) at $r=0.5$ $\mu$m, we have $\lambda\approx 0.8$ $\mu$m ($\mathrm{Pe}\approx 0.62$), i.e., a uniform arrangement, while, at $r=15$ $\mu$m, we have $\lambda\approx 0.1$ $\mu$m ($\mathrm{Pe}\approx 10$), i.e., this represents a non-uniform arrangement. The experimental investigations performed for regimes of $\mathrm{Pe}\ll 1$ Thies-Weesie et al. (1995) and $\mathrm{Pe}\gg 1$ Segre et al. (2001) revealed different types of particle arrangement. Computer simulation methods in different Péclet regimes have also been applied to study the steady state sedimentation of particles with spherical Royall et al. (2007); Padding and Louis (2004) and non-spherical Wachs (2009); Choi et al. (2013); Karimnejad et al. (2018); Singh et al. (2018) shapes.

There is special interest in the sedimentation effects in suspensions of impenetrable rodlike particles (rods). For example, electrically conductive and transparent films prepared from sediments of carbon nanotubes are of particular interest in the production of electrodes for super-capacitors, thin film transistors, and fuel cells Hu et al. (2010).

For systems of elongated particles, self-assembly and phase transition effects can introduce supplementary complications to the sedimentation regimes. For example, experimental studies have revealed the presence of jamming and orientational ordering in colloidal rod sedimentation Mohraz and Solomon (2005). A density-driven isotropic-nematic (IN) phase transition in three dimensional (3D) homogeneous system of rods with infinite aspect ratio (length-to-diameter ratio $a=\infty$) was theoretically predicted in the 1940s Onsager (1949). The theory predicted coexisting isotropic and nematic phases between particle number densities $\rho_{i}\approx 3.29$ and $\rho_{n}\approx 4.19$, and a transition to the nematic phase with strong ordering at $\rho\geq\rho_{n}$. Spontaneous ordering in concentrated rod suspensions has been experimentally confirmed in many studies Adams et al. (1998); Maeda and Maeda (2003); Alargova et al. (2004).

In two-dimensions (2D), the Onsager theory predicts a continuous IN transition at a critical density of $\rho_{i}=3\pi/2\approx 4.7$ Kayser and Raveché (1978). Monte Carlo (MC) simulations have revealed the IN transition as the number density of rods, $\rho$, increases Frenkel and Eppenga (1985); Bates and Frenkel (2000). The ordered phase becomes absolutely unstable with respect to disclination unbinding at $\rho<\rho_{i}$, but a transition to a quasi-nematic phase at $\rho\geq\rho_{n}\approx 7.25$ occurs. The quasi-nematic phase demonstrates algebraic order (quasi long-range order) and the occurrence of a disclination-unbinding transition of the Kosterlitz–Thouless (KT) type has been suggested. MC simulations of continuous 2D systems of rods revealed a 2D nematic phase of the KT type at high densities for relatively long rods with the high aspect ratio of $a\gtrsim 7$ Bates and Frenkel (2000). For non-equilibrium 2D systems of rods obtained using a random sequential adsorption (RSA) model, further self-assembly is possible owing to deposition-evaporation processes or to the diffusion motion of particles. Several problems related to such types of self-assembly of rods have previously been discussed Ghosh and Dhar (2007); López et al. (2010); Lončarević et al. (2010); Matoz-Fernandez et al. (2012); Kundu and Rajesh (2013); Kundu et al. (2013); Lebovka et al. (2017).

Pioneering work on computer simulation of the sedimentation of dilute dispersions of randomly oriented rods was performed in 1959 Vold (1959). It was demonstrated that the volume fraction of rods in the sediment decreases significantly with increase of the aspect ratio, $a$, of the particles. Self-assembly can be important at late stages of the sedimentation when high density layers are formed at the bottom. For example, pronounced density oscillations near the base have been observed in system of rods Biben et al. (1993).

Sedimentation of non-Brownian rods complicated by the long-range nature of multibody hydrodynamic interactions has intensively been studied in theoretical, experimental, and simulation works (for a review, see Guazzelli and Hinch (2011)). The mechanism of instability and cluster formation in suspension of non-spherical sedimenting particles was described in Koch and Shaqfeh (1989). Theoretical calculations predicted that the hydrodynamic interactions could result in coupling between the particles orientations. Therefore, the flow field generated by the particles can lead to a clustering of the particles and enhancement of the sedimentation rate.

The experimental investigation of the sedimentation of a dilute non-Brownian fibers within a viscous fluid (at very low Reynolds number, $\mathrm{Re}\rightarrow 0$) have been performed under well stirred conditions Herzhaft et al. (1996). The aspect ratio of fibers was $a\approx 10$ and effective number density was $\rho=0.09(2/l)^{3}$ (here $l$ is the length of the particle). The well-stirred suspensions were unstable, formation of particle clusters was observed; within a cluster, individual particles tend to align in the direction of gravity. These investigations were continued for the fibers with different aspect ratios ($a=5$–32) and different concentrations of suspensions ($=(0.019-13)\times(2/l)^{3}$) Herzhaft and Guazzelli (1999). The different regimes of sedimentation were identified. The velocity fluctuations were found to be large and anisotropic. Formation of large-scale streamers from clusters of fibers during early stages of the sedimentation was experimentally observed Metzger et al. (2005, 2007). The streamers presented highly correlated regions of downward velocities. The ‘clumps’ of parallel non-Brownian fibers settled more rapidly than individual fibers.

The MC Mackaplow and Shaqfeh (1998), dynamics simulations Butler and Shaqfeh (2002); Saintillan et al. (2005), and theoretical kinetic model Helzel and Tzavaras (2017) also predicted the formation of clusters and migration of fibers into narrow streamers aligned in the direction of gravity at $\mathrm{Re}\rightarrow 0$. The 3D dynamical simulations has been performed to study the sedimentation of large prolate non-Brownian spheroids of aspect ratio $a=5$ at small, but nonzero Reynolds number ($\mathrm{Re}=0.3$) Kuusela et al. (2001, 2003). The results revealed that, in the diluted suspension, the inertial effects tend to align particles horizontally (i.e., perpendicular to gravity), whereas, with increasing concentration, the orientation of particles becomes opposite (i.e., parallel to gravity). These sedimentation data were qualitatively supported by the experimental investigations of sedimentation of non-Brownian fibers at $\mathrm{Re}\rightarrow 1$ Salmela et al. (2007).

A computational hydrodynamic 2D model was applied to simulate the suspensions of rigid thin fibers in a viscous incompressible fluid with a nonzero Reynolds number and in absence of Brownian motion Wang and Layton (2009). The data confirmed the fiber clustering in the small Reynolds number regime; by contrast, the dispersion instead of aggregation was observed in the high Reynolds number regime. The large scale 3D numerical simulations of the gravity induced sedimentation of non-Brownian fibers in a highly viscous fluid revealed existence of the densification phase (where the cluster densifies and grows), and the coarsening phase (where the cluster becomes smaller and less dense) Gustavsson and Tornberg (2009).

The Brownian motion randomizes the orientational ordering of fibers, however, at large concentration, the Onsager type of nematic ordering can also be important due to the entropic reasons Vroege and Lekkerkerker (1992); Philipse (1997). The Onsager theory and its Parsons–Lee extension have been applied for calculation of the concentration profiles of rods with aspect ratios over the range $a=10$–80 in a gravitational field Allen et al. (1999). The calculations revealed nematic ordering at the bottom of the sediment and a strong dependence of the calculated profiles on the value of $a$. The MC simulations were applied to investigate sedimentation and phase equilibria in suspensions of hard spherocylinders with aspect ratio of Savenko and Dijkstra (2004). The investigations revealed that sedimentation led to multi-phase coexistence of isotropic and different liquid crystalline phases. However, previous studies of sedimentation of Brownian rods have not paid much attention to the effects of the particle density and sedimentation rate on the kinetics of sedimentation, on evolution orientational ordering and the electrical conductivity of the sediments.

This paper analyzes the self-assembly of Brownian rods during sedimentation in 2D colloidal suspensions. The MC simulations have been applied. Note that the MC technique is not aimed at generating dynamics. However, in the limit of long times, MC produces the trajectories equivalent to trajectories of a system under Brownian dynamics. Generally the standard MC is more fast (about one order of magnitude) than the Brownian dynamics Saintillan et al. (2005). Moreover, the proper rescaling the MC accounting for the acceptance rate of simultaneous trial displacements and rotations allows a direct comparison between MC and Brownian dynamics simulations Sanz and Marenduzzo (2010); Romano et al. (2011); Patti and Cuetos (2012); Cuetos and Patti (2015). Recently such approach have been applied to study the rodlike particles in the isotropic phase Corbett et al. (2018). The initial state was produced using an RSA model with isotropic orientations of the rods, after which the sedimentation was started and the rods were allowed to undergo both translational and rotational diffusion. The kinetics of the changes of structure and the electrical conductivity in the horizontal $x$ and vertical $y$ directions have been analyzed.

The rest of the paper is constructed as follows. In Sec. II, the technical details of the simulations are described, all necessary quantities are defined, and some test results are given. Section III presents our principal findings. Section IV summarizes the main results.

II Computational model

The initial state before sedimentation was produced using an RSA model Evans (1993). Rods of length $l$ and thickness of $d$ with a large aspect ratio, $a=l/d\gg 1$ were deposited onto a plane randomly (both their positions and orientations were random) and sequentially until the desired initial number density $\rho_{i}$ (i.e., number of rods per unit area) were obtained. Their overlapping with previously placed particles was forbidden. An isotropic initial orientation of the rods was assumed. Since the aspect ratio is large, we ignored the finite width of the rods when looking for the rod intersections, i.e., in this case, the rods were treated as zero-width objects.

The length of the system was $L$ along the horizontal $x$-direction and periodic boundary conditions were applied along this axis. The height of the system along the vertical $y$-direction was $H$. In the present work, all calculations were performed using $L=H=32l$. Zero flux boundary conditions were applied at the upper and lower borders.

The simulation of the sedimentation assumed simultaneous Brownian motion of the rods and their downward movements. The Brownian diffusion of rods was simulated using the MC procedure. At each step, an arbitrary rod was randomly chosen, and its rotational and translational diffusion motions were taken into consideration. The rotational diffusion coefficient was calculated as $D_{r}=k_{B}T/f_{r}$, the translational diffusion coefficients were calculated as $D_{\parallel}=k_{B}T/f_{\parallel}$ and $D_{\perp}=k_{B}T/f_{\perp}$ (compare with Eq. (3) for spherical particles) for the motions along and perpendicular to the direction of the long axis, respectively. Here $f_{r}$, $f_{\parallel}$, and $f_{\perp}$ are the Stokes’ friction coefficients for rotational and translational motions. For long rods, $a\gg 1$, the following formulas can be used Löwen (1994)

[TABLE]

where $\eta$ is the viscosity of the surrounding liquid, and $\gamma_{r}\approx-0.662$, $\gamma_{\parallel}\approx-0.207$, and $\gamma_{\perp}\approx 0.839$ are the end correction coefficients.

The amplitudes of the Brownian motions are proportional to the square root of the corresponding diffusion coefficients. The amplitude of the displacement along the long axis was put as $\Delta r_{\parallel}=\alpha l$, where $\alpha$ was chosen to be small enough ($\alpha=0.05$) in order to obtain satisfactory acceptance of the MC displacement Landau and Binder (2014).

The other amplitudes were evaluated using the following equations

[TABLE]

For example, at $a=10^{3}$, we have (see Eqs. (5)–(6)) $\Delta r_{\perp}/\Delta r_{\parallel}=0.76$ and $l\Delta\theta/\Delta r_{\parallel}=2.37\alpha$.

One MC time step ($\Delta t_{MC}=1$) corresponds to the attempted two displacements and one rotation for the all rods in the system. This time increment corresponds to the Brownian dynamics time increment that can be evaluated as Patti and Cuetos (2012):

[TABLE]

where $\mathcal{A}_{i}$ is the acceptance coefficient for the $i$-th MC step.

The total Brownian dynamics time was evaluated as the sum

[TABLE]

where $t_{MC}$ is the MC time.

The sedimentation rate of a rod, $u$, depends upon the angle $\varphi$ between the long axis and direction of the gravity (Fig. 1). It can be estimated from a balance of the net sedimentation force with the opposing translational frictional force Schuck (2016), and the maximum value is reached when $\varphi=0$

[TABLE]

Here $m=\Delta\rho V$ is the buoyant mass of the rod, $V=\pi l^{3}/(4a^{2})$ is the volume of the rod, $\Delta\rho$ is the difference between the densities of the particle and the liquid (compare with Eq. (2) for spherical particles).

In general case, when $\varphi\neq 0$, the sedimentation rates in the horizontal ($x$) and vertical ($y$) directions were different

[TABLE]

and at each step the displacements of the chosen particle were $\Delta x=lu_{x}/u$ and $\Delta y=lu_{y}/u$. Finally, the Péclet number can be evaluated as the ratio of the Brownian, $\tau_{B}$, and sedimentation, $\tau_{S}$, times required for the motion of a rod over the distance of $l$, i.e.,

[TABLE]

At the initial moment before settling the front of the settling layer coincides with the upper border and it moves in the downward direction in the course of the settling. Time counting was started from the value of $t_{MC}=1$, being the initial moment (before settling and diffusion), and the total duration of the simulation was typically $10^{6}$–$10^{7}$ MC time units.

Figure 2 presents an example of the evolution of the patterns at a sedimentation rate of $u=10^{-3}$ and an initial number density of $\rho_{i}=3$. See Supplemental Material at [URL will be inserted by publisher] for an animation of the sedimentation for different values of the initial number density. To calculate the running height of the front of the settling layer, $H_{f}$, the film was divided along the vertical axis into regular columns with widths $L/32$, then the maximum value of the $y$ coordinates of all rods in each column was determined, and the value of $H_{f}$ was defined as the mean value averaged along all the columns. The relative height of the front of the settling layer was defined as $h_{f}=H_{f}/H$. The averaged number density of the particles in the settling layer with a height below $H_{f}$ can be calculated as $\rho_{f}=\rho_{i}/h_{f}$.

The mean order parameter in the settling layer was calculated as

[TABLE]

where $\theta_{i}$ is the angle between the axis of the $i$-th rod and the horizontal axis $x$, and $N$ is the total number of rods. To characterize the orientation of rods at different heights, we utilized $S(y)$, i.e., the order parameter (11) calculated within a narrow band $y,y+\Delta y$, where $\Delta y=l$. The profiles of the normalized number density, $\rho^{*}(y)=\rho(y)/\rho_{i}$, and the order parameter, $S(y)$, along the vertical axis $y$ were evaluated as the averaged values of $\rho/\rho_{i}$ and $S$ in the layers with height of $y$ ($y\leq H_{f}$).

Figure 3 presents an example of the profiles of the normalized number density, $\rho(h)/\rho_{i}$, and order parameter, $S(h)$, along the vertical axis $y$ at different sedimentation times, $t_{B}$. Here, $h=y/H$ is a relative height. These data were evaluated at a sedimentation rate of $u=10^{-3}$ and an initial number density of rods of $\rho_{i}=3$. At the initial moment of time ($t=1$) the distribution of rods in the layer was homogeneous with $\rho^{*}(h)=1$ and $S(h)=0$ inside the layer. In the course of sedimentation a densified layer is formed near the bottom with increased values of $\rho^{*}$ and $S$. The zone above $y>h_{f}$ corresponds to the supernatant without particles in it.

To characterize the sediment film, the electrical conductivity at the end of simulation was calculated. The discretization scheme of the problem used a supporting mesh with a mesh size $L/m$ ($m\geq 256$) Lebovka et al. (2018a). The cells of the supporting mesh covered by rods were assumed to be conducting while the others were assumed to be insulating. In fact, the use of a supporting mesh for the calculation of electrical conductivity is equivalent to rasterization of a structure of infinitely thin rods and the substitution of them by rods having a finite aspect ratio of the order of $a_{m}\approx m/L$. The system with infinitely thin rods corresponds to that with very large values of $m/L$. However, the applied schema of discretization generates a “zoo of lattice animals”, i.e., a set of polyominoes of different shapes and sizes, especially, for small values of $m$. This set is not complete because discretization of a rod cannot produce polyominoes of all possible shapes Lebovka et al. (2018a). Therefore, we can only use this method for estimation of the behavior of electrical conductivity. Similar approaches have previously been applied for the estimation of electrical conductivities of aligned rods Tarasevich et al. (2018a, b) and dried films of rods Lebovka et al. (2018b).

In our work, the Frank—Lobb algorithm was applied to evaluate the electrical conductivity Frank and Lobb (1988). Note, that when the electrical conductivity of a host matrix is insignificant and the sticks are allowed to intersect, the utilization of Kirchhoff network approach can be also applied Tarasevich et al. (2019). We put $\sigma_{i}=1$, and $\sigma_{c}=10^{6}$ in arbitrary units for the conductivities of the insulating and conducting sites, respectively. The bottom layers with different thicknesses, $\delta$, were selected, the two conducting buses were applied to the opposite borders, and the electrical conductivity was calculated between these buses in the horizontal, $\sigma_{x}$ and vertical $\sigma_{y}$ directions (see Lebovka et al. (2016, 2017) for details).

For each given value of $\rho_{i}$ or $u$, the computer experiments were repeated up to $100$ independent runs. The error bars in the figures correspond to the standard deviations of the means. When not shown explicitly, they are of the order of the marker size.

III Results and Discussion

Figure 4 compares sedimentation patterns at different sedimentation rates, $u$, for a fixed value of the initial number density of $\rho_{i}=3$ with relatively long time of sedimentation, $t_{B}=2\times 10^{5}$. Similar patterns were also observed at other values of $\rho_{i}$ in the interval between 1 and 7.

Comparison of the patterns allows the following preliminary conclusion to be drawn. For relatively high sedimentation rates ($u=10^{-2}$, $\mathrm{Pe}=4$), the formation of a compact bottom layer is observed. Its thickness increases with decreasing values of $u$ and at the smallest sedimentation rate ($u=10^{-4}$, $\mathrm{Pe}=4\times 10^{-2}$) no compact bottom layer can be observed owing to the presence of a sedimentation-diffusion equilibrium (SDE).

Figure 5 presents the mean normalized number density $\rho^{*}$ and order parameter $S$ in the sediment layers versus the sedimentation time $t_{B}$ for the particular value of initial number density of rods of $\rho_{i}=3$. For large sedimentation rates, $u=10^{-2}$, $\mathrm{Pe}=4$, the sedimentation includes two steps with the formation of non-equilibrium porous film (a piled state with a porous structure and poorly oriented rods) at the initial stage ($t_{B}\approx 10^{2}-10^{3}$) followed by equilibration of the porous structure and the formation of compact highly oriented films ($S\approx 1$) subsequently ($t_{B}>10^{4}$).

At smaller sedimentation rates ($u<2.5\times 10^{-3}$), only one stage of sedimentation was observed. The time of film formation $t_{f}$ increased with decrease of $u$. Moreover, the saturation level of the mean normalized number density $\rho^{*}_{m}=\rho_{m}/\rho_{i}$ and order parameter $S_{m}$ for longer sedimentation times ($t_{B}\geq 5\times 10^{5}$) decreased with decreasing values of $u$. At some critical level below $u\approx 10^{-4}$ the value of $\rho^{*}_{m}$ tend to 1. That corresponded to the establishment of a SDE state (see Fig. 4, $u\approx 10^{-4}$). Note that in the SDE state, the mean order parameter was non-zero, $S\approx 0.5$. It reflected the formation of dense oriented layers near the bottom of the sediment.

Similar dependencies, were observed for other studied values of $\rho_{i}$. The obtained data revealed some impact of the initial number density of rods $\rho_{i}$ on the sedimentation behavior. Figure 6 compares the sedimentation patterns for a long sedimentation time ($t_{B}=2\times 10^{5}$, near saturation regime) for sedimentation rates ($u=10^{-3}$ and $u=10^{-4}$) and different initial number densities ($\rho_{i}=1$ and $\rho_{i}=7$). At the relatively high sedimentation rate of $u=10^{-3}$ ($\mathrm{Pe}=0.4$), compact bottom layers were formed for both the $\rho_{i}=1$ and $\rho_{i}=7$ concentrations. However, at the smaller sedimentation rate of $u=10^{-4}$ ($\mathrm{Pe}=0.04$), the SDE states for the $\rho_{i}=1$ and $\rho_{i}=9$ concentrations were rather different. For the smaller initial number density, $\rho_{i}=1$, the upper part of the system was fairly dilute, while at $\rho_{i}=9$ the rods occupied the whole volume and formed locally oriented structures near the bottom layer.

The initial number density $\rho_{i}^{*}$ also affects the kinetics of the mean values of the normalized number density $\rho^{*}=\rho/\rho_{i}$ and the order parameter $S$ in the sediment layers (Fig. 7). The presented data correspond to the near SDE regime ($\mathrm{Pe}=0.1$) in presence of the equilibrium between the sedimentation and diffusion processes. The data on $\rho^{*}(t_{B})$ revealed that the saturated normalized number density $\rho^{*}$ decreased and the order parameter $S$ increased with increasing initial number density $\rho_{i}$.

Figure 8 summarizes the observed dependencies of the saturated values of the normalized mean number density $\rho^{*}_{m}(u)$ and the order parameter $S_{m}(u)$ in sediment layers at different initial number densities, $\rho_{i}$. It is interesting that, across a wide range of sedimentation rates the $\rho^{*}_{m}(u)$ dependencies can be well fitted with high values of the coefficients of determination ($R^{2}>0.991$) by power relation

[TABLE]

with the exponent of $\beta=0.987\pm 0.019$. Therefore, all these dependencies were almost linear over the investigated range of sedimentation rates.

However, structures with a high mean order parameter, $S\approx 1$ were only observed at large values of $u$, $u\geq 2.5\times 10^{-3}$ ($\mathrm{Pe}\geq 1$) (Fig. 8). From these dependencies we can define the two zones of sedimentation rates that approximately correspond to the manifestation of the sedimentation-diffusion equilibrium state (zone I, $\mathrm{Pe}<1$) and formation of compact sediment films with high order parameter, $S\approx 1$, (zone II, $\mathrm{Pe}>1$). At very high sedimentation rates ($\mathrm{Pe}>10$) formation of non-equilibrium piled states, i.e., porous films with poorly oriented rods were also observed (data are not presented).

The structure of the sediment at the boundary between zones I and II can significantly depend upon the initial number density of the rods $\rho_{i}$. Figure 9 presents the enlarged portions of the sedimentation patterns at $t_{B}=3\times 10^{5}$ for different initial number densities $\rho_{i}$ for two sedimentation rates $u=2.5\times 10^{-3}$ (at the boundary between zones I and II, $\mathrm{Pe}=1$ and $u=10^{-2}$ ($\mathrm{Pe}=4$). An increase in $\rho_{i}$ resulted in an increase of the height of the sediments at both sedimentation rates. However, for large initial densities $\rho_{i}>5$, significant stack-like porous structures could be clearly observed for $u=2.5\times 10^{-3}$ but these were absent for $u=10^{-2}$.

Such differences in sediment film structures can affect their physical properties. We performed tests of electrical conductivity on compact sediment films with high order parameters at sedimentation rates ($u=2.5\times 10^{-3}$ and $u=10^{-2}$). Figure 10 shows the electrical conductivity in the horizontal $x$ and vertical $y$ directions versus the mesh aspect ratio of rods, $a_{m}$, defined by the size of the supporting mesh. The data are presented at $t_{B}=3\times 10^{5}$ with examples for fixed initial number density $\rho=5$ and sedimentation rate $u=2.5\times 10^{-3}$, and different thicknesses of the bottom layers $\delta$ (Fig. 10a) and fixed $\delta$ ($\delta=0.128$) and $u=10^{-2}$, and different initial number densities over the interval $\rho\in[1,7]$ (Fig. 10b).

The similar tendencies in the $\sigma(a_{m})$ dependencies were observed in both the $x$ and $y$ directions. At relatively small values of $a_{m}$ $(<100)$, the electrical conductivity was fairly high above the level of the mean geometrical conductivity defined as $\sigma_{g}=\sqrt{\sigma_{i}\sigma_{c}}=10^{3}$. We can treat a system with conductivity $\sigma>\sigma_{g}$ as conducting while a system with conductivity $\sigma<\sigma_{g}$ can be considered to be insulating Lebovka et al. (2018a); Tarasevich et al. (2018a). The values of $\sigma_{x}$ and $\sigma_{y}$ continuously decreased with increasing $a_{m}$, while loss of percolation at electrical conductivities smaller than the value of $\sigma_{g}$ can occur at large values of $a_{m}$. This simply reflect the damage of the connectivity inside the mesh clusters with thinning of the mesh structure (increasing $a_{m}$).

At fixed $a_{m}$ the increase of $\delta$ resulted in decrease of $\sigma_{y}$ and increase of $\sigma_{x}$, i.e., anisotropy of electrical conductivity became larger (Fig. 10a). It can reflect the changes in connectivity for different thicknesses of the bottom layers.

Figure 11 shows examples of the $\sigma(\rho_{i})$ dependencies in the horizontal and vertical directions at fixed $a_{m}=1000$, $\delta=0.128$ and two sedimentation rates ($u=2.5\times 10^{-3}$ and $u=10^{-2}$). In the horizontal direction $x$, the electrical conductivity $\sigma_{x}$ exceeded the values of $\sigma_{g}$ even at small initial densities $\rho_{i}>1$. In the vertical direction $y$, the values of $\sigma_{y}$ were approaching $\sigma_{g}$ at $\rho_{i}=7$. Moreover, the increase in the sedimentation rate resulted in noticeable increase of both the $\sigma_{x}$ and $\sigma_{y}$. This evidently reflected the differences in the sediment film structure displayed in Fig. 9.

IV Conclusion

A continuous 2D model of the sedimentation in suspension of rods was studied using MC simulation. During the sedimentation, the rods underwent translational and rotational Brownian motions. Significant sedimentation-driven self-assembly of the rods was observed. Different steps were observed, with the formation of non-equilibrium porous films during the initial stage followed by equilibration and the formation of compact highly oriented films over longer sedimentation time. Two zones of sedimentation rate were observed approximately corresponding to the manifestation of the sedimentation-diffusion equilibrium state (zone I, $\mathrm{Pe}<1$); the formation of compact sediment films with high order parameter, $S\approx 1$, (zone II, $\mathrm{Pe}>1$). Within zones I and II the densities of the sediment films followed the power relation $\rho^{*}_{m}\propto u^{\beta}$ with the exponent of $\beta=0.987\pm 0.019$.

This study provides the scientific community with information regarding the structure of sediment films of conducting rods, as well as how the properties of these films depend on the sedimentation conditions. The simulation of the electrical conductivity in the horizontal $x$ and vertical $y$ directions evidenced that the electrical conductivity of the resulting sediment films displays high anisotropy with $\sigma_{x}>\sigma_{y}$. The degree of anisotropy can be regulated by changing the rate of sedimentation, $u$, and the initial number density of the rods in suspension, $\rho_{i}$.

Acknowledgements

We acknowledge the funding from the National Academy of Sciences of Ukraine, Projects No. 0117U004046 and 43/19-H (N.I.L., N.V.V.), and the Ministry of Science and Higher Education of the Russian Federation, Project No. 3.959.2017/4.6 (Yu.Yu.T.). The simulations for this work were partially performed using the computing cluster of the Taras Shevchenko National University of Kiev111http://cluster.univ.kiev.ua/eng/.