Colloidal shear thickening fluids using variable shaped functional star particles: A molecular dynamics study

Complex colloidal fluids, depending on particulates’ shapes and packing fractions, may have a wide range of shear thinning and thickening behaviors. A particular interesting way to transition between different types of such behavior is by infusing functional complex particles that can be manufactured using modern techicques such as 3D printing. In this paper, we display 2D molecular dynamics simulations of such fluids with infused star-shaped functional particles, with variable leg length and number of legs, as they are infused in a non-interacting, coarse-grained fluid. We vary the packing fraction (φ) of the system, and for each different system we apply shear with various strain rate that turns the fluid into a jammed state and rise the apparent viscosity of fluid. We demonstrate the dependence of viscosity with the particles’ packing fraction. We show the role of shape and design dependence of the functional particles towards the transition to a shear thickening


Introduction
The apparent viscosity of a fluid is an important component of various applications, such as the ones that include control of mechanical vibrations. In addition, the concept of a diverging viscosity has traditionally been underneath basic theories of glasses, as they emerge in granular and colloidal systems [1]. Nevertheless, in practice, and by using common fluids and macroscopic components, using a variable viscosity concept might not be possible or would be difficult to implement in several cases. For this reason, the magnetorheological (MR) [1,2] and electro-rheological (ER) [3][4][5][6][7][8] fluids emerged, that have the ability to control the fluid viscosity but needs external magnetic and electric fields, respectively. For a continuous control of viscosity, and by using a common fluid [9][10][11] but without applying external force fields, would require additive particulates that could provide such a viscosity control. For that reason, in this paper we focus on simulations of controlling the fluid viscosity [12][13][14][15][16][17][18] through infusing complex-shaped functional particles [19][20][21][22] (like the Ferrous particles infused in MR or ER fluid) based on the concept of jamming [22][23][24][25][26][27][28][29][30] that takes place in randomly packed ensembles and therefore control the damping characteristics [31][32][33] of generic fluids through 3D printed functional particles. The idea of fluid viscosity increase lies behind studies of jamming in particulate matter. Jamming could occur both in frictional [28] or frictionless particles [34] interaction. The Geometric-Asperity(GA) model [22] that has bidisperse disks having asperities centred at disk edge considers the asperity interactions that makes the system mechanically stable and isostatic. In this paper, functional particles consist of spheres that bring asperities throughout their surface and asperity collisions with the colloidal fluid are considered during Stochastic Rotation Dynamics (SRD) runs.
There are several existing approaches to study the rheology for particle suspensions in fluids. Among them, there are dissipative particle dynamics (DPD) [35,36], Brownian dynamics [37], the lattice-Boltzmann method [38,39],Stokesian Dynamics [40][41][42], SRD [18,43,44] etc. In order to take into account temperature equilibration but in an approximate, fast manner, in this paper we utilize SRD that is available in the simulation software LAMMPS [45,46]. The focus in this work involves the calculation of the fluid viscosity as the shape and packing fraction of the infused particulates are added. Selecting particular particle shape can provide particular mechanical responses [47]. The emergent strength of packings made out of these functional particles depends on their shape and they are stiffest/softest for compact and rod-like shape respectively. Star particles in this work are formed through a computer algorithm that adds particulates stochastically towards a macroscopic target geometry ( Figure 1) and the shape we choose is stiffest at centre and softest at outer sides. Strain stiffening happens due to soft packing at the outer side near jamming.
Soft rod-like shapes are often used in nanoparticles coating. Such coated nanoparticles have many applications in biomedical science such as drug delivery, molecular imaging, biosensing etc. Studies have been performed to understand the effect of coated surfaces with various patterns for nanoparticles on their effectiveness in such applications. For example, Li et al. studied four different pattern formed due to coatings using such line shaped ligands around nanoparticles [48]. They found that, among those four patterns anisotropic ligand pattern around the nanoparticle enables easy translocation of coated nanoparticle. The dynamics of such penetration confirms the effect of surface pattern in such coated nanoparticles. However, Bogart et al reported that mixtures of multiple ligand are important for biocompatibility purposes [49]. In addition to such surface designing in nanoparticles, underlying challenges regarding material charracterization as well as the magnetism in nanaoparticles. Such material charracterization is further carried out by Lehn et al. who studied one such case with gold nanoparticles for membrane penetration [50]. Such nanoparticles could provide a basis for colloidal fluids such as the ones investigated in the current work.
The main focus of this work is shear-driven transition of colloidal fluids [25] of functional particles whose leg length and number of leg determines the packing fraction and these can be varied by the algorithm that is used for particle generation. Due to jamming, the system transformed from liquid-like state to a rigid disordered solid-like state. Other shapes like spherical or circular [51,52], ellipsoidal [53][54][55], rods [56][57][58], Polyhedra [59] and U-shaped [25] particles have already been studied considering either their packing or jamming. But the effects of more complex shape particles in jamming remain an interesting open question.
In this paper, we utilize molecular dynamics simulation using LAMMPS [45,46], performed by infusing star shaped particles with various number of legs and leg-lengths into particles that act like background fluids. We perform analysis of shear stress and viscosity with respect to strain rates to find that the liquid like state is transforming into a solid like state with a drastic increase of the apparent fluid viscosity. From the pressure vs packing density and pressure to strain rate vs strain rate we identify approximate jamming packing fractions, where the transition to shear thickening takes place. The higher the packing fraction the higher the total strain and time required to reach a steady state. Most importantly we show the existence of a jump of viscosity after reaching the jamming transition for different functional particle shapes. Finally, the structure of the system is analyzed and the main results are that the particle diffusion is consistent with a shear thickening fluid, and with the longer the leg length the larger the primary peaks in structure's pair correlation function. The rest of the paper is organized as follows. In Section 2 the molecular model and simulation method will be described in detail. In Section 3 we present detailed MD simulation results and discussion is followed by our summary in Section 4.

Materials and Methods
In this work, we use two-dimensional periodic simulations in LAMMPS with side length of 30 (LJ units). We choose SRD particles that act as the background solvent for rigid colloidal star particles. The mass of each SRD atom is set at 1.0 and the mass of each star particle atom is set at 100. We use a packing algorithm to generate packings of non-overlapping star particles using small circular core particles each having a radius r=0.2. The overall packing algorithm contains two steps: In the first step, we perform 1000 trials to place core particles randomly in a circular disk to create the core of star shaped particles with radius R=0.65 that rejects particles that are overlapping or are too close. Then, legs are placed equidistantly around the R-circle with an angle (θ) using the formula θ = 2π/Nleg, where Nleg is the number of star-legs (i.e. 3,5,7). The number of core particles in each leg (n) is an integer number determined by the formula n=(Lg/(2*r) + 1) where Lg is the desired leg length. The radial distances are measured by rx = ((R + r) *Cos(iθ) and new particles are added at a distance rx a using the formula that gives the variable leg length (L=0.0 to 6.0 for 3 leg, L=0.0 to 3.0 for 5-leggeds, L=0.0 to 2.5 for 7-leggeds) as showed in Figure 1 (a, b, c). In the second step, 30 star particles of same number of legs and each leg with same leg lengths are created by this exact same way are randomly dispersed in the square shaped simulation box, allowing the overlap of star particles like Figure 1 (d, e, f). The packing fraction (φ) for each different system is calculated by the area covered by the total number of small circular particles over the area of square where the particles are dispersed by the formula φ = π · r 2 · L 2 where N is the total number of small particles and L is the side length of simulation box. SRD particles are coarse-grained solvent constituents consisting of ideal, point particles. SRD does not interact with each other but with different shaped colloidal particles. The fluid like property is generated by collision and rotation properties of these point SRD particles, through particle streaming and particle velocity updates. The fundamental relation between mass of SRD particles (m SRD ), temperature (T), characteristic time step (∆t) and mean free path (λ) between collisions are related through the following relation λ = ∆t √ k B · T/m SRD . During the streaming step the i th particles position at time n+1 is calculated using an Euler scheme, r n+1 where ∆t is the SRD time step that equals to the simulation timestep multiplied by the number of timesteps after which SRD particle velocities reset. The particle velocity is updated using a rotation matrix R[ξ(r n+1 where u(ξ, t) = ∑ k ξ v k /N ξ,t is average velocity of bin ξ, and N ξ,t is the total number of SRD particles at time t [60]. For the simulations, we use a soft potential for star-star as well as star-SRD particles interaction. It is a repulsive potential which can push apart the overlapping atoms and also valid when r goes to 0. The potential is defined as, where A is the pre-factor ( ) and r c is the cut-off (σ). In the simulations, we use the prefactor to ramp up from 0 to 100 towards configuration equilibration, and kept the global cut-off equals 1.12 σ and local cut-off equals 1.0 σ. After the completion of the equilibration run, the pre-factor remains at its final value during shear simulations.
The key aspect of this work is the investigation of shear-induced effects on fluid viscosity. The system is sheared for different strain rates starting from 0.001 to 0.009 and 0.01 to 0.09. The negative pressure tensor in the diagonal gives the normal stresses. Shear stress is then calculated as following where δ ij is the unit tensor being 1 if i=j and 0 if i =j. The viscosity is related to the shear stress and shear strain at point p and time t as following, where the total pressure is P = 0.5 * (P xx + P yy ) and the Viscosity is then calculated from the slope of shear stress vs strain rate as the strain rate approaches zero. For diffusion coefficients, we calculate the mean square displacement (MSD) of a group of atoms. For calculating the diffusion coefficients of molecules (star particles composed of several number of relatively fixed-position atoms) we calculate the centre of mass for each molecule and then consider that the centre of atoms has MSD: where N is the number of total atoms, r n (0) = r 0 and r n (t) is the initial position and the position after time t of atoms respectively. The radial distribution function, g(r), shows the structure of the suspended star particles with different shapes. The g(r) is calculated for different initial conditions of same leg length and number of legs.

Shear thickening due to functional particle infusion
A key effect of the presence of the functional particles in the fluids is in the response to applied shear strain rates. At small strains, the overall stress is close to zero and represents a zero-shear viscosity state. But an increase in the strain rate can non-linearly increase the stress leading to sufficiently high values. Figure 2(a-c) indeed shows that, the slope of stress vs strain is zero near a very low strain rate ( 0.001). However, the gradual increase in the strain rates up to a strain rate of 0.09 shows that the slope has reaches to a comparatively larger value than the slope for low strain rates. This shows that, the shear stress increases more rises compared to an increase in strain rate. In addition to that, we observed the leg lengths are also important in this viscosity increment. For a system with 3-legged star-particles, if the leg length is small ( 0.0 -1.5) the statement of having zero shear stress at nominal stress is true. But for same strain rate, if the leg lengths are larger ( 6.0), we see an observable shear stress even at zero strain rates. This is even more noticeable if we consider the 5-legged and 7-legged star particles with large leg-lengths.
A standard way of investigating the flow stress of such complex fluids is through the fitting to stress-strain-rate relationship. Because of a non-linear increase in shear stress by infusing star-shaped particles in fluid-like SRD, we can determine the gradual rise of yield stresses in our system by using a non-newtonian fluid model [19,61] of the form where A is the slope, α is the exponent and σ y is the yield stress. Figure 2(d-f) shows the yield stresses for 3-legged, 5-legged and 7-legged star particles with different leg lengths that varies the overall packing fractions. When α is less than 1, shear thinning behavior is observed, while for α > 1 shear thickening is expected.
In the model we studied, the fluid is newtonian for small packing fractions, but develops a yield stress as packing fraction increases, as seen in Figure 2(a-c). One can see the rise of solid-like properties in those fluidic systems. At the same time, these plots show the increase towards a shear-thickening fluid for large packing fractions and large number of legs. The same conclusion can be made from Figure 2(d-f) where log-log plot of normalized shear stress are plotted against strain rate. The slope of the plot indicates how much net stress above yield stress for different system that contains a different packing density of infused star particles with respect to different strain rate. We have found from these that the normalized shear stress is positive for a system with higher packing density even at a lower strain rate and the stress sees further increment for an increment of strain rate. To demonstrate the fluid-like to solid-like transition in jammed state, we also determined how the viscosity varies with strain rate for various packing fractions and number of legs. Figure 3(a-c) shows that, with an increase of strain rate the viscosity also increases which agrees to our previous result showed in figure 2(a-c) where we demonstrate that shear stresses increase with strain rates. The increase of viscosity raises the solid-like properties in our simulated system. This increase of viscosity makes our fluid of interest as a characteristic shear-thickening fluid. The shear thickening behaviour is similar for all the 3-legged, 5-legged and 7-legged star particle system as expected. However, the magnitude of viscosity is observed to be different in different systems where particle shapes are different. Characteristically, the fluid becomes deeply shear thickening as the number or/and length of legs increase.

Pressure, viscosity, diffusivity and the onset of jamming
A way to investigate the jamming onset in the fluid arises through the tracking of the steady-state average of total pressure as a function of packing fraction. At strain rate close to zero the total pressure represents the pressure along the yield stress curve [25]. In principle we expect to see P vanish near jamming packing density and then rise to a finite value above jamming. In Figure-4(a-c) we found that the pressure begins to converge and then diverges to finite values for packing fraction (φ) in between 0.5 to 0.6 for systems infused with 3, 5 and 7 -legged star particles. As it is well known [22,29], it is actually difficult to locate the exact jamming packing density value from tracking features in the curve as simulation system size scaling is mandatory, and critical scaling analysis is necessary to determine the jamming transition [17,62,63]. By using characteristic power law fitting forms P = B(φ − φ c ) β when φ > φ c and P = 0 when φ ≤ φ c where B is slope, we estimate β as the power and in this way, we identify the jamming packing fraction φ c to demonstrate the non-monotonic evolution with respect to strain rate. The power is always greater than 1 for all configurations, indicating shear thickening behavior. The jamming packing fraction increases with respect to the increase in strain rate for 5-legged and 7-legged suspensions. Finally, another way to look at quasi-static behaviour and determine the jamming packing fraction is to investigate the pressure vs strain rate behavior, for fixed packing fraction. For φ < φ c , we get P/γ value to a finite value asγ reaches to zero and for φ > φ c we find that P/γ diverges asγ reaches zero. In this case, we find a lower bound for φ c =0.54 for 3-legged (Fig-4a), φ c =0.54 for 5-legged particles (Fig-4b) and φ c =0.6 for 7-legged ones (Fig-4c).  These plotting is having error bar by calculating standard error of mean over 30 different initial configurations for each packing density. From these plots we can clearly see that after a certain total strain value all the system reach at a steady state as we expected earlier. But for a higher packing fraction, the relaxation time as well as total strain both value is much higher than lower packing for each type of star particles, indicating the glassiness features expected for such frustrated fluids. Further, the average viscosity is investigated for each leg length with 30 different initial condition from the average of slope for first 6 and 8 points of shear stress vs strain rate plots. This average result eliminates the initial condition dependence as well as any specific factor that may lead to the jamming state. The viscosity is nearly zero until the packing fraction reaches to 0.5 for systems with 5-legged and 7-legged star particles whereas 0.6 for system with 3-legged star-particles. The similar jump behaviour can be explained in terms of the leg lengths of these particles. When the leg lengths of 5-legged and 7-legged particles reaches to 1.5 whereas the 3-legged particles reaches to 2.5, the viscosity rises significantly due to jamming. The jamming reaches to a critical value with a higher packing fraction as well as leg lengths because of the less compact shape of respective 3-legged particles compared to the 5-legged and 7-legged particles. This sharp increase of viscosity is the key to crossover of this paper that indicates the transition from liquid-like state to solid-like state is caused by the jamming of star particles. It is worth noting that the star particle shape has a direct effect on diffusivity. Diffusivity is calculated from the slope of mean square displacements for small srd atoms and big star particles centre of mass as described in the methodology section. Figure-8(a-c) shows that star particles with longer legs have high packing fraction and more inability to disperse freely without colliding and jamming with other star particles. The lower the arm length the lower the restrictions for free movement that results higher diffusion of big star particles. As the star particles are infused into the fluid therefore fluid atoms also could not move freely if it finds longer stars with longer legs and it causes the lower diffusion of small atoms surrounded by long leg stars rather than small legged stars. The change in suspension structure for different shaped star particles and packing fraction can be visualized by plotting the radial distributive function g(r). The centre of mass for star particles are considered and plotted for 3,5 and 7-legged star particles in Figure-8(d-f). For 3 leg star particles the greater number of distinct peaks indicates highly ordered structure, but it reduces into larger primary peaks as the arm length increases and the system gets more jammed. Similarly, for 5-legged (Figure-7e) and 7-legged ( fig 7f) the primary peaks becomes much larger and more local ordering as the as the arm lengths increases.

Discussion
The infused star particles are the primary element in our work that causes shear thickening in simulated SRD particles with fluid-like properties. We observed few key features, related to jamming, to draw conclusions related to the shear stresses, pressure, viscosity and diffusivity of our system that varies with applied strain rates as well as packing fraction of the infused star particles. The observation of shear stress is important to identify the shifting of our fluid-like system towards shear thickening behavior. Our results of stress with respect to strain rate clearly shows the increase of shear stress with respect to strain rate that leads to an increase of yield stress. We have also observed the linear relationship of stress and strain rate for systems with high packing density. The pressure and jamming analysis in our result section clearly indicates a jammed state caused by the infused star particles at a critical packing density. The solid-like property has shown an increase after this jammed state is achieved in different types of systems i.e. 3-legged, 5-legged and 7-legged star particle system. The key finding of our work is the rise of viscosity with respect to the packing fraction of the infused star shaped functional particles. The viscosity shows a rise after the packing fraction reaches to 0.5 for 5-legged and 7legged star particles whereas 0.6 for 3 leg star particles. This shows the shape dependency towards the rise of viscosity due to jamming transition. The 3-legged star particles are comparatively less complex structure than 5-legged and 7-legged star particles. Therefore, the viscosity rises at a lower packing fraction for a comparatively complex shaped particles.
Finally, diffusivity of particles shows that diffusion is consistent with shear thickening fluid characters. The particle size is analyzed in this study where we found that small fluidic particles are more diffusive compared to big star particles at any packing density. the probability of finding the star particles increased with reduction of diffusivity.

Conclusions
In this study, we simulated, by using MD, fluid-like SRD particles with infused starshaped particles with variable leg numbers as well as leg lengths. For various strain rate applied to the simulated system, we found that the star particles infused in fluid causes a smooth newtonian fluid to a shear-thickening transition. Higher shear strain with larger packing density causes the shear stress to reach a maximum level. The transition from Newtonian to shear-thickening fluid is occurred through jamming and the packing fraction where this transition takes place is the critical packing fraction that ranges from 0.5 to 0.6 for different shaped star particles. We have also found the system reaches to an equilibrium for after straining with a constant strain rate for sufficiently long time by studying the pressure of the system with constant strain rate. The larger the packing density or leg-lengths, the more amount of time as well as straining are required. The diffusivity of big star partilcles as well as small fluid-like SRD particles also show that, jamming causes the reduction in diffusivity which is ordered in less-compact shaped 3-legged particles and disordered in more-compact shaped 5-legged and 7-legged particles. Finally, we found a sharp rise in the viscosity of the fluid due to the infusion of star particles that causes jamming in the system. The sharp rise in viscosity is observed after the packing fraction surpasses to 0.5 for compact shaped star particles. The dependence of viscosity with respect to packing fraction for other kinds of compact shaped particles could also be of great interest, but it is beyond the scope of this work.