a r X i v :c o n d -m a t /0603667v 1 [c o n d -m a t .s o f t ] 24 M a r 2006Shear viscosity and shear thinning in two-dimensional Yukawa
liquids
Z.Donk´o 1,J.Goree 2,P.Hartmann 1,and K.Kutasi 11Rearch Institute for Solid State Physics and Optics of the Hungarian Academy of Sciences,H-1525Budapest,P.O.Box 49,Hungary 2Department of Physics and Astronomy,The University of Iowa,Iowa City,Iowa 52242,USA (February 6,2008)Abstract A two-dimensional Yukawa liquid is studied using two different nonequi-librium molecular dynamics simulation methods.Shear viscosity values in the limit of small shear rates are reported for a wide range of Coulomb cou-pling parameter and screening length.At high shear rates it is demonstrated that this liquid exhibits shear ,the viscosity ηdiminishes with increasing shear rate.It is expected that two-dimensional dusty plasmas will exhibit this effect.
Typet using REVT E X
Many-particle systems characterized by the Yukawa potential include a variety of physical dusty plasmas,charged colloids,astrophysical objects and high energy density matter.The Yukawa potentialφ(r)∝Q exp(−r/λD)/r models a Coulomb repulsion that is exponentially suppresd wit
h a screening lengthλD.Yukawa systems behave like a liquid when the temperature exceeds a melting point which depends on Q,λD,and particle spacing, e.g.[1,2].
Transport parameters of Yukawa systems–the diffusion coefficient[3],the shear viscosity [4–6]and the thermal conductivity[6,7]–have mainly been calculated for3D systems,but there is now an increasing interest in2D ttings.For example,in dusty plasma experiments, charged microspheres suspended as a monolayer in a gas discharge make a2D Yukawa system.By creating a shearflow in such a particle suspension,the viscosity was measured in recent experiments using2D suspensions[8]and quasi-2D suspensions consisting of a few monolayers of charged microspheres[9].The transport properties of such ultrathin liquids are also of interest as macroscopic analogs of molecularflow in nanoscience applications[10].
Transport coefficients are meaningful if they are part of a valid“constitutive relation”between the gradients of local variables andfluxes.For shear viscosityη,the constitutive relation j y=−η[dv x(y)/dy]relates a momentumflux j y to the velocity gradient dv x(y)/dy, which is also termed the shear rate.In a non-Newtonianfluid,ηmay vary with the velocity gradient,whereas in Newtonianfluids it does not.In particular,ifηdiminishes as shear is incread,thefluid is said to exhibit“shear thinning”.This occurs in simple liquids[11],as well as in complex mixtures such as foams,micelles,slurries,pastes,ge
ls,polymer solutions, and granularflows[12].Recently,experimenters have claimed to obrve shear thinning in dusty plasma liquids[9].The reports motivate our simulations to arch for the prence of shear thinning in2D Yukawa liquids.
Subquent to the experimental measurement of viscosity in a2D dusty plasma[8],a 2D molecular dynamics simulation was ud to obtain the shear viscosity from the Green-Kubo relations[13].In this Letter we will go beyond the results of Ref.[13],which were performed for equilibrium conditions,by using here nonequilibrium simulations to arch for
non-Newtonian behavior under conditions of a high shear rate.We will also compute the viscosity over a wider range ofΓandκ.
Our simulations u a rectangular cell with edge lengths L x and L y and periodic bound-ary conditions.The number of particles is between N=990and7040.The system is characterized by dimensionless parametersΓ=Q2/4πε0ak B T andκ=a/λD,where a=(1/nπ)1/2is the Wigner-Seitz radius,with n being the areal density.Additional pa-rameters include the thermal velocity v0=(2k B T/m)1/2and the2D analog of the plasma frequencyωp=(Q2/2πε0ma3)1/2,the shear rateγ=dv x/dy,and its normalized value
dt =
p i
dt
=F i,(1)
where r=(x,y),p=(p x,p y)are the positions and the momenta of particles,m is their mass, and F i is the force acting on particle i,are integrated by the velocity Verlet algorithm.
Method2simulates a planar Couetteflow,which is established by the Lees-Edwards periodic boundary conditions which result in a homogeneous streamingflowfield in the simulation box: v x =γ(y−L y/2),where denotes a time average.The system is described by the Gaussian thermostated SLLOD equations of motion[16]:
d r i
m +γy iˆx,
d˜p i
|j y|=ηeq dv x(y)/dy=∆p/2t sim L y,(3) where∆p is the total x-directional momentum exchanged between slabs A and B during the simulation time t sim[15].In method2the off-diagonal element of the pressure tensor is measured during the cour of the simulation:
P xy(t)=
N
i=1
mv ix v iy+N
j>i
x ij y ij
dr ij
φ(r ij) ,(4)
where r ij=r i−r j=(x ij,y ij),and the shear viscosity is obtained as
η=lim
t→∞
P xy(t) /γ.(5) In method1,the spatial profiles for temperature and velocity,Fig.1,develop lf-consistently in respon to the perturbation applied by introducing momentum in slabs A and B.We u method1only for small perturbations,so that the velocity profile has a linear gradient and the temperature is isotropic,with T x=T y,where T x,y= (m/N j k B) N j i=1 [v ix,y(t)−
v jy is negligibly small.
Obtaining reliable results forηat smallγrequires a simulation duration of typically ωp t∼104–105for both methods.The required time step is smallest,and the simulations are most costly,at lowΓ.In method1,system size effects are expected to appear when (i)particles traver the simulation box without significant interaction with the others,or (ii)the compressional sound wave transits the box in a shorter time than the decay time t c of the velocity autocorrelation function.For(i),for our most demanding condition(small size N=990and high temperatureΓ=1)a particle moving at the thermal velocity would transit the cell in a timeωp t≈57if it were undeflected by collisions.Wefind that the decay time is short enough,ωp t c∼5–10,even for the smallestΓvalues of interest.Thus we expect no“
ballistic”trajectories across the entire simulation box.For(ii),the sound speed[14]atκ=1is v=dω/dk∼aωp,and the wave’s transit time∆t across a box with length L y isωp∆t∼L y/a=57for our most demanding ca,N=990particles.Thus,we find both criteria fulfilled for a“sufficiently large”system.Method2is known to produce
accurate results even for small number of particles simulated[16].We verified that the results obtained from both methods did not depend significantly on N.
Figures2(a)and(b)illustrate particle trajectories in simulations bad on method2, for conditionsΓ=10,κ=1at a shear rateγ=0.05, respectively.
Our results forηeq as a function ofΓ,for different values ofκare plotted in Fig.3(a).We find a good agreement with the earlier equilibrium MD simulation of Ref.[13].In contrast with most simple liquids,which have a viscosity that varies monotonically with temperature, a prominent feature of the viscosity of the prent system is a atΓ∼=20 forκ=1),which has been noted previously in both OCP(one-component plasma)and Yukawa liquids.The shape of theηeq(Γ)curve can be explained by the prevailing kinetic and potential contributions to the viscosity at low and high values ofΓ,respectively.The near-equilibrium shear viscosity values obtained with method2forκ=1are also displayed Fig.3(a).Wefind an excellent agreement between the results of methods1and2.
Similar to what was obrved in[5]for3D Yukawa liquids,wefind that the near-equilibrium viscosityηeq obeys a scaling law as demonstrated in Fig.3(b),where viscos-ity has been normalized byηE=mnωE a2.The Einstein frequencyωE depends onκ,and we computed it from Eq.(7)of Ref.[14]using pair-correlation functions measured from our simulations.The horizontal axis is a normalized temperature T′=T y/T m=Γm/Γ, where T m andΓm are melting-point values reported in Ref.[2].Using the normalizations, the data fall on the same curve,demonstrating the existence of a scaling law for the0.5≤κ≤2.0range of the screening parameter.We note that for this purpo we foundωE was more significant thanωp.The near-equilibrium viscosity isfit by an empirical form ηeq/ηE=aT′+b/T′+c with coefficients:a=0.0093,b=0.78and c=0.098.
A shear-thinning effect is revealed in Fig.4(a),which shows thatηdiminishes significantly as the shear rate
γ[11]. Wefind that this scaling also occurs for the Yukawa system,as indicated by data that fall
