• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    Hydrodynamic Stress Tensor in Inhomogeneous Colloidal Suspensions: an Irving-Kirkwood Extension?

    2019-07-25 02:02:20ZongLiSun孫宗利YanShuangKang康艷霜andYanMeiKang康艷梅
    Communications in Theoretical Physics 2019年7期

    Zong-Li Sun (孫宗利), Yan-Shuang Kang (康艷霜), and Yan-Mei Kang (康艷梅)

    1Science and Technology College,North China Electric Power University,Baoding 071051,China

    2College of Science,Hebei Agricultural University,Baoding 071001,China

    3College of Chemistry and Environmental Science,Hebei University,Baoding 071001,China

    4University of International Relations,Beijing 100091,China

    Abstract Based on statistical mechanics for classical fluids,general expressions for hydrodynamic stress in inhomogeneous colloidal suspension are derived on a molecular level.The result is exactly an extension of the Iving-Kirkwood stress for atom fluids to colloidal suspensions where dynamic correlation emerges.It is found that besides the inter-particle distance,the obtained hydrodynamic stress depends closely on the velocity of the colloidal particles in the suspension,which is responsible for the appearance of the solvent-mediated hydrodynamic force.Compared to Brady’s stresslets for the bulk stress,our results are applicable to inhomogeneous suspension,where the inhomogeneity and anisotropy of the dynamic correlation should be taken into account.In the near-field regime where the packing fraction of colloidal particles is high,our results can reduce to those of Brady.Therefore,our results are applicable to the suspensions with low,moderate,or even high packing fraction of colloidal particles.

    Key words: hydrodynamic stress,colloidal suspension,diffusion,hydrodynamic force

    1 Introduction

    In either equilibrium or non-equilibrium of the colloidal suspensions,stress tensor is of supreme importance in the analysis of its static or dynamic properties,including the interfacial tension,transport coefficients,viscoelasticity,and so on.[1]In their pioneering work,Irving and Kirkwood (IK) derived the stress tensor in atom fluids on a microscopic level.[2]Though the uniqueness of its representation is still in debate,[3]it has been extensively employed to describe both static and dynamic aspects of the fluids.[4?8]In contrast,investigation of the stress in colloidal suspensions is rather primitive.Actually,the stress in colloidal suspension is quite different from that in atom fluids,because of the more complex correlation and multi temporal scales in the suspension.This also makes the investigation on stress in suspension quite tedious.[9]Nevertheless,due to its close relation to the viscoelastic behavior,stress tensor in colloidal suspension has being a fascinating topic in the studies on its rheological and transport behaviors.

    Besides the direct inter-particle correlation,the hydrodynamic interaction (HI) is indispensable in the determination of the stress tensor in the dynamics of suspension.However,it has no counterpart in the equilibrium stress.Intrinsically,the hydrodynamic force on the colloidal particle is nothing but an additional force arising from the flow of solvent,which is caused by the motion of the surrounding particles.[9?10]In fact,a huge body of researches on the stress tensor in the bulk colloidal suspensions have been conducted since the early work of Batcheloret al.,[11]who introduced firstly the effect of HI on the bulk stress of suspension.Moreover,in the study on the creeping flow on the Smoluchowski time scale,Felderhof has expressed the HI contribution of stress tensor by using of a mobility matrix.[12]In addition,in the study of the osmotic pressure of a colloidal dispersion,Brady has constructed the HI contribution of the bulk stress by a set of hydrodynamic resistance tensor.[13?14]Actually,all of these derivations are designed for the bulk suspension.In the confined geometry,however,confinement can dramatically modulate both structural and viscoelastic properties,due either to the interface or to the finite-size effects in the systems.[15?17]Therefore,in investigations on mass,momentum or energy transport,the inhomogeneity and anisotropy of the stress in the colloidal suspensions should be taken into account.

    Actually,in the past decades,dynamics of confined colloidal suspension,especially in the over-damped limit,has attracted considerable attention.[18?19]By solving Langevin’s equations of motion in phase space or the Smoluchowski equation in configurational space,one can achieve the dynamic evolution of the colloidal particles in the suspension.[20?22]Recently,an even more promising approach,say dynamic density functional theory(DDFT),[23?27]has been proposed to describe the dynamics of colloidal particles.In most DDFT formalism,the hydrodynamic force is usually described on a two-particle level,such as that proposed by Rotne and Prager.[28]More importantly,the DDFT contains most of the dynamic information,though some approximation is required in the construction of the DDFT.Therefore,one can resort to DDFT for the structure and motion information of the dynamic suspension,which is indispensable in calculation of the hydrodynamic force and stress.

    Though hydrodynamic correlation has achieved satisfactory consideration in DDFT,its contribution to the stress tensor in colloidal suspension is still less satisfactory,especially in the inhomogeneous suspension.Actually,the original IK stress tensor can be extended to the fluid mixture.[29]However,when the size deference between components is large enough,such simple extension fails to capture all the correlation information in the mixture,because the contribution from hydrodynamic force is not included.In fact,in dense colloidal suspensions,hydrodynamic force may become dominant,and its contribution to the stress becomes important and non-negligible.Therefore,how to take the contribution from hydrodynamic force into account is crucial in the investigation of the colloidal suspensions.In view of this,the hydrodynamic stress in inhomogeneous colloidal suspension will be calculated in this work,and how it relies on the motion of the particles will be discussed.

    This paper is organized as follows.In Sec.2,the IK stress tensor in multi-component fluid mixture is calculated based on the equation of momentum balance.In this section,the IK results are generalized to the colloidal suspension,and the hydrodynamic stress is defined after the determination of the solvent-mediated hydrodynamic forces.In Sec.3,the velocity of the colloidal particles is determined with the aid of the DDFT,and the hydrodynamic stresses in diffusive suspension,with and without shear flow,are calculated.A comparison of our results with those of Brady’s stresslets has also been performed in this section.A discussion on the obtained hydrodynamic stresses is given in Sec.4.

    2 Stress Tensor in Multi-Components Mixtures

    Here,we consider a system of colloidal suspension,which is composed of the solvent molecular and colloidal particles.The former of the massmsand radiusas,while the latter of massmcand radiusac.In their dynamic process,they take instantaneous streaming velocitiesus(r,t) anduc(r,t),respectively.Firstly,a volumeV,with boundaryS,is arbitrarily chosen from the suspension.Then,the time rate of change of the total momentumP(t) inVcan be expressed as:

    withρν(r,t) the instantaneous mass density.Actually,there are usually two contributions to the change of the momentum in the volume.One arises from the momentum flow acrossS,while the other from the action of some force on it.The former can be written as:

    while the latter is given by:

    with Π the stress tensor in the suspension.By utilizing of the Gauss’s integral theorem,it is ready to express the momentum continuity equation for the binary mixture as:

    It should be noted that the above equation is applicable for an arbitrary multi-components systems,and that the stress tensor Π(r,t) should include all the correlations between the components,either static or dynamic correlations.

    Actually,Eq.(4) is nothing but the definition of the stress tensor in the suspension,just as its application in the system of simple atom fluid.For the dynamic colloidal suspension under consideration,however,it may become more complex due to the inclusion of dynamic correlations.Even so,it is believed that the stress tensor is still feasible so long as the correlations in the suspension are known.Therefore,in the following,we proceed to express the stress tensor through calculating the time rate of change of the momentum

    To take into account the contribution of stress tensor from the thermal motion,it is convenient to separate the laboratory velocity of thei-th particle of componentν,into a thermal partvν,iand a streaming partuν(ri,t),that is:

    withrν,ithe position of the particle.Due to the fact that the thermal velocity does not contribute to the momentum current,the time rate of change of the momentum can be rewritten as:

    which can further lead to the following:[30]

    Comparing with Eq.(4),it is obtained that

    Here,fν,idenotes the force exerting on thei-th particle of componentν,arising from its interaction with other particles.Accordingly,the second term on the right hand side of Eq.(8) can be rewritten as:

    Further,using the fact

    and combining Eqs.(8) and (9),it leads to:

    withIthe unit tensor.Here,the first term on the right hand side is the kinetic contribution,while the second term denotes the contribution from the intra- and inter-component interactions.In Eq.(10),the stress contribution from the forcefνν′,ijcan be given as:

    withrνν′,ij=rν′j?rνi.Note that the stress contribution in Eq.(11) has been written in the form of summation,which may facilitate its application in the computer simulation.Actually,for the convenience in the theoretical study,one can usually rewrite it,by performing ensemble average,as:

    Some detail about the derivation of Eq.(12)can be found in Appendix A.

    Evidently,by substituting the direct inter-particle interaction forcefνν′,ij≡??φ(rνν′,ij),Eq.(12) is exactly the celebrated IK formula for the potential contribution to the stress.Nevertheless,fνν′,ijshould never be limited to the direct inter-particle forces.It may also arise from the hydrodynamic interaction,because the hydrodynamic force can also impose an additional impact on the change of momentum in the suspension.In view of this,the IKlike formula should also be applicable for the contribution from the hydrodynamic force,which plays the same role as the direct inter-particle force does.

    In the colloidal suspension,significant difference between the mass or size of solvent molecular and those of colloidal particle is responsible for the multi temporal scales in the system.Accordingly,on the Brownian time scale,the solvent flow is so fast that the solvent degreeof-freedom can usually be integrated out in the consideration of the motion of colloidal particles.[9?10]This is the physical interpretation of the hydrodynamic force between colloidal particles.In view of this,the general expression for the stress in colloidal suspension can be written in the following form:

    Here,the kinetic contribution

    is known exactly.Further,by regarding the suspension as a mixture,the contribution from the direct inter-particle interactions can be decomposed into four parts,that is,each of which can be directly obtained by Eq.(12).

    It should be noted that Πdir(r,t) in dynamic suspension should have the same origin as those in the equilibrium phase,because the direct inter-particle forces are independent on the dynamics of the suspension.In contrast,the contribution from HI between colloidal particles,say Πhyd(r,t),is peculiar in the dynamics,and has no counterpart in the equilibrium stress.Therefore,in the following sections,we focus only on the hydrodynamic stress Πhyd(r,t).Note also that in the following,the subscriptsα,βwill be used to distinguish the colloidal particles dispersed in suspension.

    3 Hydrodynamic Stress in Colloidal Suspension

    3.1 Velocity of Colloidal Particles in the Framework of DDFT

    As mentioned above,significant difference in the temporal scales makes it possible to ignore the dynamical details of the solvent,because the solvent can be regarded as continuum on a Brownian time scale.Besides,the momentum of the particles can also be relaxed to the thermal equilibrium with the heat bath of the solvent.Consequently,the dynamics of the colloidal particles can be described only in the coordinate space by the Smoluchowski equation:[9]

    with the total number of colloidal particleNcand the probability distribution function?c(rNc,t).In addition,the microscopic diffusion tensor for the colloidal particle can be constructed as:[31]

    withδαβthe Kronecker delta.The first two terms on the right hand side of Eq.(15) describe the self diffusion,while the remaining describes the collective diffusion.More specifically,D(0)=D0I,withD0=kBT/6πμacand the shear viscosityμ,is the self diffusion coefficient in the dilute limit case,in which the hydrodynamic interaction is absent.represents the effect of its own motion that reflected by the surrounding particles.D(2)(rα,rβ)=represents the effect arising from the motion of the other particles.With the method of reflection,one can express:

    with?=ac/rαβand=rαβ/rαβ.

    In Eq.(14),the total potential energy Φ(rNc,t) in Eq.(14) is given,up to a pairwise level,as:

    with the external potentialVext(rα,t)and the pair potentialφ(rα,rβ).

    On the Brownian time scale,the hydrodynamic force is balanced with the resultant of Brownian and direct interparticle forces.Accordingly,for the colloidal particleα,it reads

    Here,=??αΦ(rαβ) is the direct force,while=is the Brownian force.[9]The total hydrodynamic forceexerting on the particleαrelates to the motion of particles through:

    Note thatuαis the diffusion velocity of the particleαwith respect to the fixed reference system.This provides us a feasible route to determine the motion of the particles.Combining Eqs.(17) and (18),the diffusion velocity can be given by:

    By integrating both sides withand using the first two members of the Yvon-Born-Green hierarchy,[32]one obtains the ensemble average of the diffusion velocity,for the particle locating atrα,as:

    in whichg(2)c(rαβ) is the radial distribution function.Note that in the overdamped limit,the non-equilibrium Helmholtz free energy functional,sayF[nc],can be written in the same form as that of the corresponding equilibrium system,that is,

    whereFexc[nc]is the excess free energy functional contributed from the interaction between particles.[32?33]In Eq.(19),it has been indicated that the determination ofis a complex many-body problem,because the dynamics of each particle relates simultaneously to all the surrounding particles.Fortunately,dynamic density functional theory,which is originally designed to calculate the dynamic evolution of the particle density,provides an effective way to determine the diffusion velocityudif(rα),just as shown in Eq.(20).

    3.2 The Hydrodynamic Stress in Diffusive Suspension

    As is known,colloidal suspension is a multi-scaled fluid mixture.On the Brownian time scale,however,the solvent should be regarded as continuum,because the Brownian time scale is long enough for the solvent to relax their momentum and position coordinates to the thermal equilibrium.However,the particle in a viscous solvent can feel the motion of the surrounding particles due to the solvent-mediated hydrodynamic force between them.In contrast to the direct force,the hydrodynamic force is more complex in that it depends not only on the interparticle distance,but also on the velocity of particles in the suspension.Therefore,in order to determine its contribution to the stress,one should firstly find out the explicit dependence of hydrodynamic force on inter-particle distance and particle velocity.

    Fig.1 (Color online)A sketch for the decomposition of ,which is exerted on the particle α by the surrounding colloidal particles.In the figure,the blue arrows represent the velocity of particles,while the red arrows stand for different contributions of the solvent-mediated hydrodynamic force.In addition,the pink arrow represents the reflected velocity field by the particle γ.

    In the resistance problem of two translating sphere,the total hydrodynamic forceon particleαcan be obtained by using of the method of reflections.[34]For instance,the result for,up to the fourth reflection,can be given by:

    withf0α=?γuα.Note thatγ=6πμacdenotes the friction coefficient of an isolated sphere withμthe shear viscosity.This indicates that for the two-body resistance problem,the hydrodynamic force on particleαhas three contributions,which arise respectively from the friction force of the isolated sphere,the reflected field of particleαby particleβ,and the field of particleβ.It is expected that this type of decomposition forcan be generalized to the many-body resistance problem.In view of this,one can decomposeinto three different types of contribution,that is

    A sketch for the decomposition ofis given in Fig.1.Specifically,in Eq.(23),denotes the hydrodynamic force of an isolated particle,which can be given as:

    with Υ(0)=γI.As stated above,is exactly the friction force in a very dilute suspensions,where the hydrodynamic force is absent.

    In contrast,represents the contribution from theα-th particle’s movement,which is reflected from the particleγand back to the particleα.It is given as:with the resistance tensor Υ(1)(rαγ).Obviously,depends not only on the distance between particlesαandγ,but also on the relative velocity between them.In addition,denotes also the hydrodynamic force exerted on the particleαby the solvent,but stems from the motion of the particleβ.This contribution can be given as:

    with another resistance tensor Υ(2)(rαβ).Actually,the above-mentioned Υ(0)and Υ(1)are the diagonal elements of the resistance matrix Υ,while Υ(2)is the off-diagonal elements of it.Actually,simple expressions for Υ(1)and Υ(2)can be obtained by comparing the results in Eqs.(24)–(26) with those given in Eq.(22),that is,

    Note that in the above,two types of the second-order tensor,sayDand Υ,are used to describe the dynamics of colloidal particles in the suspension.Actually,they describe,respectively,the mobility and resistance problems of the colloidal in the suspension.In general,they are related with each other through Υ?1=D/kBT,which indicates thatDis the inverse of the entire matrix Υ.

    According to the decomposition ofin Eq.(23),the hydrodynamic stress Πhydcan also be decomposed into:

    Substituting Eqs.(25)and(26)into Eq.(12),one obtains the last two terms,respectively,as:

    Note that in the above equations,andhave been regarded as pairwise forces,which are mediated by the solvent.

    From above,it is obvious thatanddepend closely on the velocity of other particles.In contrast,the friction forcerelates only to the motion of its own.Therefore,the friction force can be treated as a thermodynamic force,whose contribution to the stress tensor can be given by:

    Inserting the expression forand performing ensemble average,can be rewritten as:

    Some details about the derivation of Eq.(32)can be found in Appendix B.Collecting the results in Eqs.(29),(30),and (32),one obtains the hydrodynamic stress contribution from the interactionas:

    Note that Eq.(33) gives the stress contribution from the hydrodynamic force in a diffusive suspension,where the solvent has been considered as a continuum.Further,in the dilute limit,there is no hydrodynamic interaction.In this situation,the last two terms in Eq.(33) disappear,and only the first term remains.

    In practice,at any moment of a dynamic process,the involved pair density of the fluid can be further expressed as:n(2)c(rα,rβ)=nc(rα)nc(rβ)g(rαβ),withg(rαβ) the pair distribution function.In the framework of DDFT,this is also a requirement to close the non-linear equation for the one-particle distribution.[18,24]Actually,g(rαβ)in non-equilibrium fluid can usually be approximated by that of an equilibrium fluid with the same one body density profile,which has been proven to be effective in investigations on non-equilibrium properties of diverse fluids.[18?19,23?25,27,29]

    3.3 Hydrodynamic Stress in Suspension with a Shear Flow

    In the above,we have focused on the hydrodynamic stress tensor in a diffusive colloidal suspension,particularly on that arising from the inter-particle hydrodynamic force.When the suspension is subjected to a shear flow,the involved hydrodynamic interaction may become even more complex because of the non-uniform flow velocity,sayE ·r,in the solvent.Here,E=(1/2)[?v+(?v)T]is the symmetrical part of the rate of strain tensor of the solvent,andvdenotes the velocity field of the solvent.In addition,the flow around each particle can be distorted by the presence of the other particles,which leads to the flow disturbance.Accordingly,the hydrodynamic force on the particleαcan be expressed as

    withCαthe disturbance tensor of rank three,anduβthe velocity of colloidal particleβin the shear flow.Here,the termCα:Estands for the involved hydrodynamic force arising from the disturbance in the solvent.

    On the Brownian time scale,the above hydrodynamic force on a given particle is balanced by other forces exerting on it,that is

    Accordingly,the relative velocity of the particleαwith respect to the shear flow,defined as=uα?E ·rα,can be given by:

    with

    In Eq.(36),it is obvious thathas two contributions,which arise respectively from the sum of inter-particle and Brownian interactions,and the disturbance in the solvent.Further,one can calculate their ensemble averages to obtain local values ofand,that is,Accordingly,the local relative velocity of particle atrαcan be determined by:

    Actually,urel-1(rα)is exactly the result given in Eq.(20),but the involvednc(r) should be the particle density for the suspension with a shear flow.

    Similar to the case of diffusion without shear flow,the hydrodynamic force of particleαin suspension with a shear flow can also be decomposed into three contributions,which arise respectively from the relative motion of particleαwith respect to solvent,the reflected motion of particleαby the surrounding particles,and the motion of other particles.That is,

    On the other hand,bothurel-1(rα) andurel-2(rα) contribute to the stress in the suspension.Consequently,one can denote their contributions to the stress as Πsf-1and Πsf-2,respectively.After a similar calculation as performing in the previous section,one can obtain Πsf-1(r) as:

    Comparing with the case without shear flow,it is obvious that this equation can be rewritten as

    which suggests that Πsf-1(r) is nothing but the diffusion term of the hydrodynamic stress.Following a similar treatment,the stress contribution fromurel-2(r) can be obtained as:

    In contrast to Πsf-1(r),Πsf-2(r) is a disturbance term in the hydrodynamic stress.The results in Πsf-2(r) indicate that the presence of shear flow has complicated the expression of the hydrodynamic stress,through the effect of the flow disturbance embedded inurel-2(r).

    However,this is not the whole.In the above calculation for the hydrodynamic stress,we have taken only into account the relative velocity of particle with respect to the flowing solvent.In fact,even in the dilute limit,the shear flow itself should contribute to the stress in the suspension,which can be given,on the one-particle level,as:

    Here,RSEis a rank-four tensor,and its magnitude is known as|RSE|=(20/3)πμa3c.[10]By summing Πsf-1(r),Πsf-2(r),and ΠE(r),one can finally express the hydrodynamic stress in a shear flow,as:

    3.4 Comparison with Brady’s Stresslets in a Shear Flow

    In suspension mechanics,the concept of stresslet has been a useful tool for understanding the rheology of suspension and its mixture.It is the symmetric portion of the first moment of the force distribution in suspension.Though the stresslet is unnecessary in the equations of motion for a particle,it is an important bridge between the force distribution and the motion of the suspension.By introducing a set of hydrodynamic resistance tensor,Brady had proposed the expressions for the bulk stress in a suspension with shear flow:[13?14]

    with

    The hydrodynamic stressletis a moment of the fluid traction force taken over the surface of a particle.Therefore,it represents the average of the continuum stress inside that particle.The inter-particle stressletarises from the solvent-mediated hydrodynamic forces,while the Brownian stresslet 〈SB〉 stems from the thermal fluctuations of the solvent.Actually,Brady’s result has obtained wide usage in the investigation of viscoelasticity in suspensions.[35?38]Therefore,in the following,we will compare our results with Eq.(45) to show that our derivation includes all the contributions that presented in Eq.(45).

    Actually,within the mean field scheme,one can simplify(r) by takingξ=0 in Eq.(30),that is,

    withrβ=rα+rαβ.Keeping in mind the fact that Υαβ=RFU,and thatrαβ→in the near-field limit,(rα) can be rewritten as:

    Performing similar calculation for(rα) and(rα),and collecting the results for them,one can obtain

    which is consistent with the results given inandNote that Eqs.(17) and (23) have been used in derivation of Eq.(50).Further,it is obvious that Πsf-2(r)in Eq.(42) denotes the contribution from the flow disturbance.According to the above analysis,one can rewrite it as:

    which is also consistent with the first term shown in,with the fact thatCα=RFE.[9]

    From the above analysis,it is obvious that within the mean field scheme,the results shown in Eqs.(41)–(44)can reduce to those of Brady for the bulk stress,though they seem quiet different.Evidently,our results can be used in a wider scope in that they can be easily applied to the suspension in a confined geometry,which is beyond the scope of Brady’s result.Accordingly,it is expected that our results can help to deepen our understanding about the hydromechanical and rheological properties of the colloidal suspension,especially in the cases where the hydrodynamic interaction is dominant.

    3.5 Contribution from Hard Sphere Interaction:on a Molecular Level

    As is known,an additional Brownian force is necessary in the construction of the Smoluchowski equation,because an inhomogeneous density in ideal suspension can always evolve spontaneously to a Boltzmann distribution.Such a Brownian force is usually expressed as

    Actually,by regarding the colloidal suspension as a mixture,some information may be missed if only the Brownian force of colloidal particle is taken into account.Assume that the chosen volumeVconsists ofNccolloidal particle andNssolvent molecular,whose probability density functions are respectively?c(rNc,t) and?s(rNc,t).The total probability distribution function can be written as

    with the total particle numberN=Nc+Ns.In order to include the hard sphere interaction of the solvent molecular,we introduce an additional force similar asfB

    Accordingly,by performing an ensemble average with?(rN,t),its contribution to the stress tensor can be expressed as

    Note that,here,the particleimay be either colloidal particle or solvent molecular.Further calculation can be performed for ΠA

    Note that the second integral term gives exactly the numberN.In addition,because of the holes created by the impermeable colloidal particles or solvent molecular,some new surface,saySh,emerges in the suspension.To include this effect,one can transform the volume integral to a surface integral by using the Gauss’s theorem.Actually,the probability density of each particle vanishes at the boundaryS.Accordingly,one obtains

    withdij=ai+ajdenoting the distance between particlesiandjat contact.Here,we have transformed the surface integral back to the volume integral.=(rj?ri)/dijis the unit vector at the boundary point where the particlesiandjmake contact.By using the fact that

    Equation (58) can be rewritten as

    Introducing the center-of-mass coordinaterc=(rj+ri)/2 and relative coordinaterij=rj?ri,it can be further rewritten as

    In this equation,?kBTδ(|rij|?ai?aj)is the interparticle force arising from a hard sphere potential.Note that in the above equation,the indexi,jcan be arbitrary.That is,{i,j}={c,c},{s,s},{c,s}or{s,c},each of which corresponds to one part of the stress tensor.This coincides with our results for Πdir,which has been given in sec.2.Therefore,instead of stresslet,the contribution from the short-ranged inter-particle and particle-solvent interactions can be taken into account through the hard sphere interaction potential.

    In fact,from the above analysis,it is convenient to justify whether any ignorance or repetition of some parts of the stress has taken place in the calculation,especially for the contributions from short-ranged interactions.Obviously,Eq.(60) has taken into account all the possible contributions from the hard sphere interaction,even that from the colloidal-solvent hard sphere interaction,which is usually neglected due to its infinitesimal magnitude in comparison with other contributions.By comparison between Eq.(60) and that in Brady’s results,[13?14]it is obvious that the{s,s}contribution has been included in the solvent pressurep,while the{c,c}contribution has taken into account by the Brownian force.Therefore,introduction of the additional forcefAis equivalent to the consideration of the inter- and intra-components hard sphere interaction in the expression of Πdir.

    4 Concluding Remarks

    In this work,the hydrodynamic stress tensors in diffusive colloidal suspension,with and without shear flow,have been derived on a molecular level.It has been found that besides the contributions from the thermal motion and direct interaction,contributions from the hydrodynamic force can also be expressed as functional of the dynamic densities.Accordingly,if the dynamic structure is known,the stress tensor in the colloidal suspension can be easily calculated.Actually,our main results for the hydrodynamic stress are neither the same as the original IK formula for the atom fluids,nor the same as Brady’s formula for the bulk stress in suspension.In the following,the characteristics of the derived hydrodynamic stress will be discussed in detail.

    For atom fluids,Schofieldet al.had claimed the nonuniqueness of the microscopic stress tensor in inhomogeneous fluids.[3]In general case,the definition of stress tensor should satisfy the equation of momentum balance,which gives Eq.(8).In the colloidal suspension,if one focuses only the hydrodynamic stress,it reads:

    Obviously,one can arbitrarily introduce a zero-divergence tensor to the hydrodynamic stress,without disturbing the above equation of momentum balance.Therefore,the hydrodynamic stress shares the same characteristics with that of the atom fluids.In addition,it is believed that in an inhomogeneous suspension,some of its viscoelastic properties,which are calculated based on the stress tensor,should suffer from the non-uniqueness of the microscopic stress tensor.The relevant research on this issue is in progress.

    Nevertheless,there are still differences between the stress tensor in colloidal suspension and that in a general mixture.Firstly,the hydrodynamic stress is peculiar in the dynamic colloidal suspension,and has no counterpart in the case of general fluid mixture.In addition,the hydrodynamic stress becomes significant only when the size asymmetry in the mixture is large enough,otherwise the hydrodynamic stress is negligible and the total stress reduces to that in the IK case.Secondly,the underlying causes for these two types of stress are different.Specifically,the stress in general fluid mixture arises only from the direct inter-particle forces,which depend on the inter-particle distance and the strength of interaction.In contrast,the stress in colloidal suspension is more complex in that besides the contributions from the direct inter-particle forces,it has additional contributions,say the hydrodynamic stresses,which stem from the solventmediated hydrodynamic forces.Thirdly,the difference between them lies also in the computation aspect.For the stress in a general mixture,the necessary inputs are no more than the pair density and the pair potential.For the stress in colloidal suspension,however,more information,including viscosity of solvent and the velocity of particles,are needed to calculate the hydrodynamic forces between particles.It is necessary to note that the velocity of particles can be obtained with the aid of DDFT.

    In addition,in the derivation of the hydrodynamic stress for suspension,we have used the general resistance tensors for the involved hydrodynamic force.Actually,the expressions for resistance tensor in far- and near-field regimes can be determined by method of reflections and lubrication theory,respectively.Actually,our derivation is performed in an arbitrary condition,and thus the results should be applicable for both far-field and near-field regimes.In this work,we have shown that in the nearfield situation,our results reduce exactly to the stresslets proposed by Brady.This indicates that these stressletbased results are more suitable for describing the dense suspensions.In view of this,our derivation is more advantageous because it is so powerful that it can be applied to dilute,moderate,and dense suspensions.Therefore,it is expected that our derivation in this work can provide new insights for the researches in the fields of viscoelasticity and rheology.

    Appendix A: Derivation of Eq.(12)

    In Eq.(12) in the text,the contribution of the stress from the interaction force has been obtained in a form or summation.However,the ensemble average is necessary in some theoretical calculation,especially based on density functional theory.Therefore,in this appendix,we will transform Eq.(12) to an ensemble average and find its equivalence to Irving-Kirkwood’s result.For the conciseness in notation,we focus only on the case of single component fluid,whose results can be immediately extended to the fluid mixture.In such a single component situation,Eq.(12) can be rewritten as:

    Introducing an auxiliary variabletand an integralEq.(A1)can be rewritten in the form of ensemble average as:

    Note that in the last step,the definition of pair density function has been used,that is,

    Transforming the auxiliary variabletback torij,the potential contribution of the stress can be written in an ensemble average as:

    which is the result obtained by Irving-Kirkwood.It is obvious that Eq.(A3) can be extended immediately to the case of multi-components fluid mixtures.

    Appendix B: Derivation of Eq.(32)

    In the diffusion of colloidal suspension,two types of diffusion processes are involved: self diffusion and collective diffusion.The former concerns only the dynamics of the a single particle,while the latter relates to the motion of many particles simultaneously,as shown in Eq.(15).As stated in the text,the hydrodynamic force on each particle can be decomposed into three parts:,,.Among these three parts,andare on a two-particle level,whileis a force on a single-particle level.Therefore,its average contribution to stress,can be written as

    两性午夜刺激爽爽歪歪视频在线观看 | 91精品国产国语对白视频| 久久人人爽av亚洲精品天堂| 亚洲欧美激情综合另类| 一级片免费观看大全| 精品国产乱码久久久久久男人| 欧美日韩乱码在线| videosex国产| 男女床上黄色一级片免费看| 久久精品国产亚洲av香蕉五月| 1024香蕉在线观看| 99精品欧美一区二区三区四区| 男女床上黄色一级片免费看| 亚洲av日韩精品久久久久久密| 99热只有精品国产| 一二三四社区在线视频社区8| 午夜福利在线观看吧| 婷婷丁香在线五月| 国产精品亚洲美女久久久| 999久久久国产精品视频| 亚洲久久久国产精品| 这个男人来自地球电影免费观看| 国产高清videossex| 少妇熟女aⅴ在线视频| 最好的美女福利视频网| 黄频高清免费视频| 久久久国产成人精品二区| 亚洲午夜理论影院| 韩国精品一区二区三区| 女性生殖器流出的白浆| 女警被强在线播放| 波多野结衣巨乳人妻| 久久久久久免费高清国产稀缺| 一进一出好大好爽视频| 禁无遮挡网站| 国产激情欧美一区二区| 日韩精品青青久久久久久| 国产黄a三级三级三级人| 亚洲成av人片免费观看| 黄色毛片三级朝国网站| 精品人妻1区二区| 亚洲av日韩精品久久久久久密| 别揉我奶头~嗯~啊~动态视频| 九色国产91popny在线| 麻豆成人av在线观看| 国产精品香港三级国产av潘金莲| 久久久久久久久久久久大奶| 久久人妻av系列| 亚洲精品一卡2卡三卡4卡5卡| 午夜福利欧美成人| 亚洲精品久久国产高清桃花| 日本三级黄在线观看| 欧美不卡视频在线免费观看 | 欧美+亚洲+日韩+国产| 国产亚洲欧美在线一区二区| 国产成人一区二区三区免费视频网站| 中文字幕av电影在线播放| 成人av一区二区三区在线看| 久久久精品欧美日韩精品| 在线观看免费午夜福利视频| 可以在线观看毛片的网站| 国产欧美日韩一区二区三区在线| 久久中文字幕一级| 多毛熟女@视频| 亚洲色图av天堂| 日韩精品免费视频一区二区三区| 又大又爽又粗| 亚洲 欧美一区二区三区| 国产高清视频在线播放一区| 久久久水蜜桃国产精品网| 99在线视频只有这里精品首页| 午夜免费鲁丝| 韩国精品一区二区三区| 日日干狠狠操夜夜爽| 色av中文字幕| 亚洲成国产人片在线观看| 高清黄色对白视频在线免费看| 性色av乱码一区二区三区2| 麻豆av在线久日| 国产色视频综合| 极品人妻少妇av视频| 熟女少妇亚洲综合色aaa.| 91九色精品人成在线观看| 又黄又粗又硬又大视频| 久久国产精品影院| 亚洲全国av大片| 99re在线观看精品视频| 欧美黑人欧美精品刺激| 精品国产乱码久久久久久男人| 满18在线观看网站| 国产亚洲av嫩草精品影院| 在线av久久热| 国产在线精品亚洲第一网站| 搡老熟女国产l中国老女人| 亚洲五月色婷婷综合| 日本欧美视频一区| 免费高清在线观看日韩| 久久天躁狠狠躁夜夜2o2o| 欧美丝袜亚洲另类 | 亚洲全国av大片| 美女午夜性视频免费| 亚洲熟妇中文字幕五十中出| 丝袜美腿诱惑在线| 亚洲五月婷婷丁香| tocl精华| 国产97色在线日韩免费| 亚洲国产中文字幕在线视频| 亚洲情色 制服丝袜| 脱女人内裤的视频| 天天一区二区日本电影三级 | 中国美女看黄片| bbb黄色大片| 性欧美人与动物交配| 级片在线观看| 国产精品亚洲美女久久久| 国产午夜福利久久久久久| 国产三级在线视频| 久久婷婷人人爽人人干人人爱 | 老熟妇仑乱视频hdxx| 琪琪午夜伦伦电影理论片6080| 精品一区二区三区四区五区乱码| 国产精品99久久99久久久不卡| 1024香蕉在线观看| 桃色一区二区三区在线观看| 国产精品免费视频内射| 国产一区二区激情短视频| 免费在线观看亚洲国产| 欧美激情久久久久久爽电影 | 性欧美人与动物交配| 搡老熟女国产l中国老女人| 国产一区在线观看成人免费| 香蕉国产在线看| 一级毛片高清免费大全| 久热这里只有精品99| 欧美乱妇无乱码| 可以在线观看毛片的网站| 亚洲中文字幕一区二区三区有码在线看 | 国产熟女xx| 欧美成人午夜精品| 午夜福利一区二区在线看| 国内毛片毛片毛片毛片毛片| 好男人在线观看高清免费视频 | 精品久久久久久久人妻蜜臀av | 欧美激情高清一区二区三区| 无人区码免费观看不卡| 看免费av毛片| 女同久久另类99精品国产91| 日本 av在线| 男女之事视频高清在线观看| 日韩一卡2卡3卡4卡2021年| 国产成人免费无遮挡视频| 每晚都被弄得嗷嗷叫到高潮| 精品乱码久久久久久99久播| 午夜福利在线观看吧| 午夜久久久在线观看| 在线十欧美十亚洲十日本专区| 成人三级黄色视频| 亚洲性夜色夜夜综合| 午夜老司机福利片| 亚洲 欧美一区二区三区| 动漫黄色视频在线观看| 精品国产美女av久久久久小说| 性色av乱码一区二区三区2| 亚洲一区高清亚洲精品| 9191精品国产免费久久| 韩国精品一区二区三区| 精品一品国产午夜福利视频| 欧美激情高清一区二区三区| e午夜精品久久久久久久| 色播在线永久视频| 91成人精品电影| 久久人妻熟女aⅴ| 无人区码免费观看不卡| 看黄色毛片网站| 女同久久另类99精品国产91| 精品人妻在线不人妻| 满18在线观看网站| 成人三级做爰电影| netflix在线观看网站| 亚洲专区中文字幕在线| 亚洲第一青青草原| 黄频高清免费视频| 91麻豆av在线| 两性午夜刺激爽爽歪歪视频在线观看 | 欧美在线黄色| 日本免费a在线| 波多野结衣一区麻豆| 精品久久久久久成人av| 亚洲国产欧美网| 曰老女人黄片| 日韩成人在线观看一区二区三区| 亚洲熟妇熟女久久| 男人的好看免费观看在线视频 | 午夜福利,免费看| 欧美激情久久久久久爽电影 | 色在线成人网| 日本精品一区二区三区蜜桃| 亚洲九九香蕉| 久久精品成人免费网站| 色av中文字幕| 亚洲va日本ⅴa欧美va伊人久久| 欧美成人午夜精品| 国产av一区二区精品久久| 久久中文字幕一级| 精品国产一区二区久久| 国产99久久九九免费精品| 少妇 在线观看| 亚洲第一av免费看| 久久久精品欧美日韩精品| 欧美日韩一级在线毛片| 精品第一国产精品| 一级黄色大片毛片| 99在线人妻在线中文字幕| 成熟少妇高潮喷水视频| 一进一出抽搐gif免费好疼| cao死你这个sao货| 香蕉久久夜色| 久久久久久国产a免费观看| 午夜成年电影在线免费观看| 91成人精品电影| 禁无遮挡网站| 两性午夜刺激爽爽歪歪视频在线观看 | 亚洲精品一卡2卡三卡4卡5卡| 亚洲无线在线观看| 啦啦啦观看免费观看视频高清 | 精品一区二区三区av网在线观看| avwww免费| 国产精品美女特级片免费视频播放器 | 国产一区二区三区在线臀色熟女| 国内精品久久久久精免费| 美女午夜性视频免费| 女警被强在线播放| 麻豆成人av在线观看| 亚洲九九香蕉| 日韩大尺度精品在线看网址 | 精品久久久久久久久久免费视频| avwww免费| 精品一区二区三区四区五区乱码| 一区二区三区高清视频在线| 男人的好看免费观看在线视频 | 9191精品国产免费久久| 免费高清视频大片| 久久午夜综合久久蜜桃| 国产不卡一卡二| 久99久视频精品免费| 久久国产精品影院| 日韩av在线大香蕉| 久久精品亚洲熟妇少妇任你| 亚洲性夜色夜夜综合| 我的亚洲天堂| 欧美午夜高清在线| 亚洲av电影不卡..在线观看| 午夜福利18| 久久久久国内视频| 手机成人av网站| 国产一区二区三区在线臀色熟女| 成人欧美大片| 99精品欧美一区二区三区四区| 欧美日韩亚洲国产一区二区在线观看| 久久久国产欧美日韩av| 99国产精品一区二区蜜桃av| 亚洲国产日韩欧美精品在线观看 | 亚洲一卡2卡3卡4卡5卡精品中文| 亚洲欧美激情在线| 一级,二级,三级黄色视频| 久久国产亚洲av麻豆专区| 一进一出好大好爽视频| 熟女少妇亚洲综合色aaa.| av天堂在线播放| 香蕉国产在线看| 免费无遮挡裸体视频| 久久久精品欧美日韩精品| 首页视频小说图片口味搜索| av免费在线观看网站| 黑人巨大精品欧美一区二区mp4| 成人国产综合亚洲| av视频免费观看在线观看| 黄片大片在线免费观看| 精品国内亚洲2022精品成人| 日本五十路高清| 一级,二级,三级黄色视频| 亚洲色图综合在线观看| 午夜影院日韩av| 国产精品爽爽va在线观看网站 | 国内毛片毛片毛片毛片毛片| 国产日韩一区二区三区精品不卡| 国产欧美日韩一区二区三区在线| 一级片免费观看大全| 在线国产一区二区在线| 婷婷六月久久综合丁香| 成人手机av| 国产成人精品久久二区二区免费| ponron亚洲| 免费女性裸体啪啪无遮挡网站| 91麻豆av在线| 欧美激情久久久久久爽电影 | av福利片在线| 91精品国产国语对白视频| 妹子高潮喷水视频| 国产欧美日韩精品亚洲av| 免费在线观看黄色视频的| 精品久久久久久久人妻蜜臀av | 精品乱码久久久久久99久播| 中文字幕人妻熟女乱码| 中文字幕精品免费在线观看视频| 亚洲第一青青草原| 欧美一级毛片孕妇| 亚洲成人久久性| 深夜精品福利| 中文字幕色久视频| 身体一侧抽搐| 国产不卡一卡二| 美女 人体艺术 gogo| 女同久久另类99精品国产91| 成人免费观看视频高清| 久久久久国产精品人妻aⅴ院| 18禁观看日本| 久久精品国产亚洲av高清一级| 最近最新中文字幕大全电影3 | 精品久久久精品久久久| 婷婷精品国产亚洲av在线| 久久久久久久精品吃奶| 成熟少妇高潮喷水视频| 激情在线观看视频在线高清| 欧美激情 高清一区二区三区| 给我免费播放毛片高清在线观看| 国产激情欧美一区二区| 亚洲精品在线观看二区| 日韩视频一区二区在线观看| 精品久久蜜臀av无| 国产精品免费视频内射| 国产高清激情床上av| 欧美成狂野欧美在线观看| 午夜日韩欧美国产| 久久人人精品亚洲av| cao死你这个sao货| 成人18禁高潮啪啪吃奶动态图| 99久久久亚洲精品蜜臀av| 日日干狠狠操夜夜爽| www.精华液| 亚洲欧美日韩另类电影网站| 高清毛片免费观看视频网站| 亚洲国产高清在线一区二区三 | 这个男人来自地球电影免费观看| 精品一区二区三区四区五区乱码| 此物有八面人人有两片| 青草久久国产| 又黄又粗又硬又大视频| 久久婷婷成人综合色麻豆| 成年女人毛片免费观看观看9| 日日爽夜夜爽网站| 在线观看www视频免费| 岛国视频午夜一区免费看| 久久久精品国产亚洲av高清涩受| 国产视频一区二区在线看| 日韩三级视频一区二区三区| 非洲黑人性xxxx精品又粗又长| 超碰成人久久| 男女下面进入的视频免费午夜 | 999久久久国产精品视频| 久久久国产精品麻豆| 首页视频小说图片口味搜索| 黑人巨大精品欧美一区二区mp4| 女性生殖器流出的白浆| 少妇的丰满在线观看| 国产视频一区二区在线看| 在线天堂中文资源库| 久久草成人影院| 亚洲 欧美 日韩 在线 免费| 91成年电影在线观看| 国产精品久久电影中文字幕| 国产精品久久久av美女十八| 亚洲第一欧美日韩一区二区三区| 精品人妻1区二区| 99久久精品国产亚洲精品| 久久欧美精品欧美久久欧美| 很黄的视频免费| 伊人久久大香线蕉亚洲五| 国产欧美日韩综合在线一区二区| 久久欧美精品欧美久久欧美| 首页视频小说图片口味搜索| 日韩精品免费视频一区二区三区| √禁漫天堂资源中文www| 日本 av在线| 精品午夜福利视频在线观看一区| 久久精品国产亚洲av香蕉五月| 国内久久婷婷六月综合欲色啪| 美女免费视频网站| 精品日产1卡2卡| 亚洲欧美日韩无卡精品| 亚洲精品国产一区二区精华液| 婷婷六月久久综合丁香| 两个人看的免费小视频| 国产午夜精品久久久久久| 可以在线观看毛片的网站| 男女床上黄色一级片免费看| 真人做人爱边吃奶动态| 日韩欧美免费精品| 久久久久国产精品人妻aⅴ院| 啪啪无遮挡十八禁网站| 成人特级黄色片久久久久久久| 一区二区三区精品91| 亚洲精品国产区一区二| 久久久精品国产亚洲av高清涩受| 天天躁狠狠躁夜夜躁狠狠躁| 熟女少妇亚洲综合色aaa.| av视频免费观看在线观看| 中文字幕最新亚洲高清| 亚洲av成人av| 国产成人精品久久二区二区免费| 51午夜福利影视在线观看| 欧美日韩瑟瑟在线播放| 日韩高清综合在线| 一区二区三区高清视频在线| 亚洲免费av在线视频| 天天躁狠狠躁夜夜躁狠狠躁| 国产精品 欧美亚洲| 欧美老熟妇乱子伦牲交| 中文字幕av电影在线播放| 日韩视频一区二区在线观看| 日韩精品中文字幕看吧| 丝袜人妻中文字幕| 纯流量卡能插随身wifi吗| 非洲黑人性xxxx精品又粗又长| 午夜福利18| 校园春色视频在线观看| 久久中文看片网| 国产99久久九九免费精品| 久久久久久大精品| 中出人妻视频一区二区| 在线观看www视频免费| 99国产精品一区二区蜜桃av| 成年版毛片免费区| 久久精品成人免费网站| 免费在线观看视频国产中文字幕亚洲| 日本vs欧美在线观看视频| 久热爱精品视频在线9| 一级片免费观看大全| 一边摸一边做爽爽视频免费| av福利片在线| 国产野战对白在线观看| 亚洲男人的天堂狠狠| 宅男免费午夜| 久久青草综合色| or卡值多少钱| 美女国产高潮福利片在线看| 这个男人来自地球电影免费观看| 久久久久久人人人人人| 国产亚洲欧美98| avwww免费| 欧美久久黑人一区二区| 国产一区二区三区在线臀色熟女| 成人国产综合亚洲| 午夜精品在线福利| 亚洲精品一区av在线观看| 欧美另类亚洲清纯唯美| 日本精品一区二区三区蜜桃| 咕卡用的链子| 亚洲精华国产精华精| 欧美亚洲日本最大视频资源| 97超级碰碰碰精品色视频在线观看| 97人妻精品一区二区三区麻豆 | 18禁黄网站禁片午夜丰满| 无遮挡黄片免费观看| 久久伊人香网站| 欧美黑人欧美精品刺激| 国产熟女xx| av超薄肉色丝袜交足视频| 欧美亚洲日本最大视频资源| 女人被狂操c到高潮| 日韩欧美一区二区三区在线观看| 国产99久久九九免费精品| 亚洲九九香蕉| 侵犯人妻中文字幕一二三四区| 午夜精品久久久久久毛片777| 欧美日本亚洲视频在线播放| 亚洲精品国产色婷婷电影| 久久久水蜜桃国产精品网| 亚洲专区字幕在线| 久久久久九九精品影院| 亚洲精品在线观看二区| 日本在线视频免费播放| 亚洲中文av在线| 老司机深夜福利视频在线观看| 久久精品成人免费网站| 97碰自拍视频| 国产成人啪精品午夜网站| 亚洲成av片中文字幕在线观看| 色综合婷婷激情| 亚洲黑人精品在线| 国产精品综合久久久久久久免费 | 嫩草影视91久久| 亚洲免费av在线视频| 他把我摸到了高潮在线观看| 热re99久久国产66热| 久久久久久久精品吃奶| 每晚都被弄得嗷嗷叫到高潮| 99久久综合精品五月天人人| 国产精品一区二区三区四区久久 | 精品少妇一区二区三区视频日本电影| 免费在线观看黄色视频的| 欧美在线黄色| 欧美最黄视频在线播放免费| 一夜夜www| 精品国产亚洲在线| 好男人在线观看高清免费视频 | 18禁裸乳无遮挡免费网站照片 | 午夜久久久久精精品| 身体一侧抽搐| 99在线视频只有这里精品首页| 激情视频va一区二区三区| 成年人黄色毛片网站| 在线观看www视频免费| 如日韩欧美国产精品一区二区三区| 亚洲精品中文字幕在线视频| 别揉我奶头~嗯~啊~动态视频| 好男人电影高清在线观看| 两人在一起打扑克的视频| 精品一区二区三区视频在线观看免费| 国产伦人伦偷精品视频| 在线观看舔阴道视频| 51午夜福利影视在线观看| av免费在线观看网站| a在线观看视频网站| 少妇 在线观看| 18禁裸乳无遮挡免费网站照片 | 美女国产高潮福利片在线看| 1024视频免费在线观看| 亚洲熟女毛片儿| 久久性视频一级片| 欧美亚洲日本最大视频资源| bbb黄色大片| 亚洲成人精品中文字幕电影| 亚洲av电影不卡..在线观看| 国产精品久久久人人做人人爽| 两个人看的免费小视频| 国产三级黄色录像| 久久 成人 亚洲| 国产三级黄色录像| 国产亚洲精品一区二区www| 久久精品国产亚洲av高清一级| 亚洲av五月六月丁香网| av超薄肉色丝袜交足视频| 50天的宝宝边吃奶边哭怎么回事| 亚洲精品国产精品久久久不卡| 色综合婷婷激情| 99国产精品一区二区蜜桃av| 日韩国内少妇激情av| 亚洲一码二码三码区别大吗| 国产精品1区2区在线观看.| 亚洲五月天丁香| 老熟妇乱子伦视频在线观看| 色播亚洲综合网| 国产亚洲精品久久久久久毛片| 欧美乱妇无乱码| 自线自在国产av| 很黄的视频免费| 一区二区三区高清视频在线| 可以在线观看的亚洲视频| 亚洲精品在线观看二区| 精品久久蜜臀av无| 国产欧美日韩一区二区三| 看免费av毛片| 午夜福利免费观看在线| 丁香六月欧美| 俄罗斯特黄特色一大片| 亚洲午夜精品一区,二区,三区| 高清毛片免费观看视频网站| 午夜视频精品福利| e午夜精品久久久久久久| 中文字幕人成人乱码亚洲影| 日韩欧美国产在线观看| 长腿黑丝高跟| 91九色精品人成在线观看| 69av精品久久久久久| 成人特级黄色片久久久久久久| 国产欧美日韩综合在线一区二区| 久久国产乱子伦精品免费另类| 亚洲精品国产一区二区精华液| 在线免费观看的www视频| 亚洲色图综合在线观看| 极品教师在线免费播放| 成年版毛片免费区| 欧美成狂野欧美在线观看| av天堂在线播放| 国产精品久久电影中文字幕| 99国产精品免费福利视频| 亚洲国产欧美网| 免费人成视频x8x8入口观看| 妹子高潮喷水视频| 国产成人av教育| 亚洲国产看品久久| 精品第一国产精品| 国产亚洲欧美在线一区二区| 99国产精品99久久久久| 日韩欧美免费精品| 精品福利观看| 一区在线观看完整版| 在线视频色国产色| 亚洲av电影不卡..在线观看| 亚洲成人精品中文字幕电影| 99久久综合精品五月天人人| 18禁国产床啪视频网站| 国产精品野战在线观看| 狠狠狠狠99中文字幕| 亚洲成人免费电影在线观看| 国产成人av教育| avwww免费| 免费在线观看日本一区| av在线天堂中文字幕| 亚洲专区国产一区二区| 国产欧美日韩综合在线一区二区| 午夜福利欧美成人| av天堂久久9| 国产欧美日韩精品亚洲av|