• <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

    久久免费观看电影| 亚洲,一卡二卡三卡| www.熟女人妻精品国产| 精品久久蜜臀av无| 久久综合国产亚洲精品| 天堂8中文在线网| 丰满少妇做爰视频| 又黄又粗又硬又大视频| netflix在线观看网站| 高清不卡的av网站| 老司机在亚洲福利影院| 久久久久久久久免费视频了| 成年女人毛片免费观看观看9 | 亚洲成人国产一区在线观看 | 欧美日韩国产mv在线观看视频| 我要看黄色一级片免费的| a 毛片基地| 蜜桃国产av成人99| 男人舔女人的私密视频| 亚洲伊人色综图| 亚洲精品久久成人aⅴ小说| 欧美激情极品国产一区二区三区| 久久中文字幕一级| 日韩av在线免费看完整版不卡| bbb黄色大片| 国产亚洲午夜精品一区二区久久| 免费在线观看日本一区| 日韩免费高清中文字幕av| 久久久久国产精品人妻一区二区| 三上悠亚av全集在线观看| 日韩电影二区| www.999成人在线观看| 国产欧美亚洲国产| 51午夜福利影视在线观看| 丝袜人妻中文字幕| 亚洲精品乱久久久久久| 高潮久久久久久久久久久不卡| 少妇猛男粗大的猛烈进出视频| 亚洲成人手机| 欧美日韩精品网址| 欧美 日韩 精品 国产| 精品国产一区二区三区四区第35| 欧美变态另类bdsm刘玥| 五月开心婷婷网| 天天操日日干夜夜撸| 日本色播在线视频| 一级毛片 在线播放| 欧美国产精品一级二级三级| 中文字幕高清在线视频| 美女国产高潮福利片在线看| 老司机靠b影院| 丁香六月天网| 一本色道久久久久久精品综合| 成人18禁高潮啪啪吃奶动态图| 午夜91福利影院| 男人爽女人下面视频在线观看| 秋霞在线观看毛片| 亚洲欧美一区二区三区国产| 欧美 日韩 精品 国产| 青春草视频在线免费观看| 看十八女毛片水多多多| 日韩av不卡免费在线播放| 国产极品粉嫩免费观看在线| 又粗又硬又长又爽又黄的视频| 欧美久久黑人一区二区| 一区二区三区四区激情视频| 久久精品亚洲熟妇少妇任你| 晚上一个人看的免费电影| 欧美精品高潮呻吟av久久| 国产一卡二卡三卡精品| 久久国产精品大桥未久av| 老鸭窝网址在线观看| 麻豆国产av国片精品| 在线av久久热| 激情五月婷婷亚洲| 午夜福利视频精品| 老司机在亚洲福利影院| 亚洲 国产 在线| av视频免费观看在线观看| 色94色欧美一区二区| 母亲3免费完整高清在线观看| 免费观看a级毛片全部| 女人精品久久久久毛片| 天堂8中文在线网| 免费在线观看视频国产中文字幕亚洲 | 久久鲁丝午夜福利片| 久久精品国产亚洲av高清一级| 久久性视频一级片| 又粗又硬又长又爽又黄的视频| av国产精品久久久久影院| 中文乱码字字幕精品一区二区三区| 欧美日韩黄片免| 日韩欧美一区视频在线观看| 久久久国产精品麻豆| 丝袜人妻中文字幕| 各种免费的搞黄视频| 久久99热这里只频精品6学生| 在线 av 中文字幕| 亚洲专区中文字幕在线| 精品一区二区三卡| 成人免费观看视频高清| 日韩大码丰满熟妇| 少妇的丰满在线观看| 久久人人爽av亚洲精品天堂| 老司机深夜福利视频在线观看 | 天天躁夜夜躁狠狠躁躁| 午夜福利一区二区在线看| 黑丝袜美女国产一区| 中国美女看黄片| 十分钟在线观看高清视频www| 久久青草综合色| 国产视频一区二区在线看| 天天影视国产精品| 19禁男女啪啪无遮挡网站| 精品福利永久在线观看| 国产麻豆69| 在线观看国产h片| 久久久久网色| 中文字幕亚洲精品专区| 久热这里只有精品99| 午夜福利在线免费观看网站| 成人黄色视频免费在线看| av天堂久久9| 免费在线观看视频国产中文字幕亚洲 | 91九色精品人成在线观看| 啦啦啦啦在线视频资源| 80岁老熟妇乱子伦牲交| 色网站视频免费| 亚洲五月婷婷丁香| 久久午夜综合久久蜜桃| 欧美日韩av久久| 天天躁狠狠躁夜夜躁狠狠躁| 19禁男女啪啪无遮挡网站| 午夜av观看不卡| 国产黄色免费在线视频| avwww免费| 一本色道久久久久久精品综合| a级毛片黄视频| 久久精品亚洲av国产电影网| 欧美激情 高清一区二区三区| 99国产精品一区二区蜜桃av | 脱女人内裤的视频| 久久国产精品影院| 男女无遮挡免费网站观看| 免费观看人在逋| 欧美日韩福利视频一区二区| 人人澡人人妻人| 水蜜桃什么品种好| 久久国产精品男人的天堂亚洲| 看十八女毛片水多多多| 亚洲美女黄色视频免费看| 天堂8中文在线网| 久久久精品区二区三区| 深夜精品福利| 亚洲欧美日韩另类电影网站| 亚洲少妇的诱惑av| 国产精品.久久久| 视频在线观看一区二区三区| 午夜激情久久久久久久| 精品国产国语对白av| 欧美成人精品欧美一级黄| 久久鲁丝午夜福利片| 国产国语露脸激情在线看| 精品少妇久久久久久888优播| 狂野欧美激情性bbbbbb| 国产精品久久久av美女十八| 亚洲九九香蕉| 女警被强在线播放| 久久影院123| 老司机在亚洲福利影院| 免费少妇av软件| 大话2 男鬼变身卡| 两人在一起打扑克的视频| 精品高清国产在线一区| 日本五十路高清| 免费观看人在逋| 制服人妻中文乱码| 国产一区有黄有色的免费视频| 国产成人影院久久av| 国产成人a∨麻豆精品| 欧美在线黄色| 午夜老司机福利片| 日韩中文字幕欧美一区二区 | a 毛片基地| 国产欧美日韩精品亚洲av| 欧美黄色淫秽网站| 久久久久久人人人人人| 国产免费福利视频在线观看| 丁香六月欧美| 欧美老熟妇乱子伦牲交| av网站在线播放免费| 一级黄片播放器| 国产精品一区二区在线观看99| 免费久久久久久久精品成人欧美视频| 亚洲精品一区蜜桃| 亚洲色图 男人天堂 中文字幕| 国产在视频线精品| 精品人妻熟女毛片av久久网站| 欧美黄色淫秽网站| 三上悠亚av全集在线观看| 国产女主播在线喷水免费视频网站| 人人澡人人妻人| 国产成人系列免费观看| 成人手机av| 又黄又粗又硬又大视频| 亚洲,欧美,日韩| 成人国产av品久久久| 欧美少妇被猛烈插入视频| 久久久久久人人人人人| 只有这里有精品99| a级毛片黄视频| 80岁老熟妇乱子伦牲交| 亚洲美女黄色视频免费看| 精品少妇久久久久久888优播| 日韩中文字幕欧美一区二区 | 人人妻人人爽人人添夜夜欢视频| 亚洲av日韩在线播放| 欧美日韩福利视频一区二区| 啦啦啦视频在线资源免费观看| 丝袜在线中文字幕| 女人久久www免费人成看片| 国产成人91sexporn| 一本综合久久免费| 男女无遮挡免费网站观看| 久久九九热精品免费| 99久久99久久久精品蜜桃| 精品视频人人做人人爽| 91国产中文字幕| 精品卡一卡二卡四卡免费| 黄色 视频免费看| 精品国产一区二区三区久久久樱花| 首页视频小说图片口味搜索 | 亚洲精品自拍成人| 欧美精品av麻豆av| 少妇粗大呻吟视频| 国产av精品麻豆| 黄色片一级片一级黄色片| 女性被躁到高潮视频| 一本一本久久a久久精品综合妖精| bbb黄色大片| 女性被躁到高潮视频| 国产xxxxx性猛交| 18在线观看网站| 天天躁日日躁夜夜躁夜夜| 精品人妻在线不人妻| 亚洲av欧美aⅴ国产| 免费av中文字幕在线| 亚洲精品在线美女| 欧美激情极品国产一区二区三区| 亚洲国产av影院在线观看| 男男h啪啪无遮挡| 新久久久久国产一级毛片| 夫妻性生交免费视频一级片| 操美女的视频在线观看| 男人爽女人下面视频在线观看| 国产精品麻豆人妻色哟哟久久| 视频区欧美日本亚洲| 99热网站在线观看| 免费日韩欧美在线观看| 蜜桃在线观看..| 成人黄色视频免费在线看| 亚洲国产最新在线播放| 在线精品无人区一区二区三| 成人影院久久| 亚洲国产看品久久| av又黄又爽大尺度在线免费看| 天天躁夜夜躁狠狠躁躁| 观看av在线不卡| 一级黄色大片毛片| 日韩人妻精品一区2区三区| 99九九在线精品视频| 国产精品二区激情视频| 中国国产av一级| 日韩视频在线欧美| 蜜桃国产av成人99| 亚洲欧美一区二区三区黑人| 在线观看国产h片| 欧美97在线视频| 少妇人妻久久综合中文| 黄色 视频免费看| 麻豆av在线久日| 一本久久精品| 人体艺术视频欧美日本| 国产免费一区二区三区四区乱码| 黄色视频不卡| 亚洲精品久久午夜乱码| 欧美激情 高清一区二区三区| 自拍欧美九色日韩亚洲蝌蚪91| 久久人人97超碰香蕉20202| 中文欧美无线码| 嫩草影视91久久| 中文欧美无线码| 一二三四在线观看免费中文在| 中文字幕另类日韩欧美亚洲嫩草| av电影中文网址| 国产欧美亚洲国产| 日日夜夜操网爽| 久久毛片免费看一区二区三区| 一本色道久久久久久精品综合| 伊人久久大香线蕉亚洲五| 亚洲色图 男人天堂 中文字幕| 少妇 在线观看| 男女免费视频国产| 女人爽到高潮嗷嗷叫在线视频| 欧美激情极品国产一区二区三区| 国精品久久久久久国模美| 久热这里只有精品99| 精品久久蜜臀av无| 国产在线观看jvid| 一本—道久久a久久精品蜜桃钙片| 亚洲 国产 在线| 成人黄色视频免费在线看| 秋霞在线观看毛片| 美女大奶头黄色视频| 国产精品久久久久成人av| 久久精品人人爽人人爽视色| 亚洲色图综合在线观看| www.av在线官网国产| 青春草亚洲视频在线观看| 亚洲国产看品久久| 国产精品免费大片| 国产精品香港三级国产av潘金莲 | 考比视频在线观看| 激情五月婷婷亚洲| 国产91精品成人一区二区三区 | 国产精品亚洲av一区麻豆| 乱人伦中国视频| 国产极品粉嫩免费观看在线| 男人舔女人的私密视频| 一区二区av电影网| 午夜免费成人在线视频| 亚洲免费av在线视频| 亚洲av在线观看美女高潮| 免费在线观看黄色视频的| 丝袜在线中文字幕| 操出白浆在线播放| 每晚都被弄得嗷嗷叫到高潮| 亚洲人成电影免费在线| 嫩草影视91久久| 欧美日韩福利视频一区二区| 伊人久久大香线蕉亚洲五| 精品国产乱码久久久久久小说| 婷婷成人精品国产| av福利片在线| 亚洲图色成人| 老司机午夜十八禁免费视频| 91成人精品电影| 亚洲av男天堂| 亚洲欧美色中文字幕在线| 亚洲精品乱久久久久久| 国产又爽黄色视频| 国产男女内射视频| 国产精品国产av在线观看| 多毛熟女@视频| 午夜免费观看性视频| 一边亲一边摸免费视频| 久久亚洲精品不卡| 国产片内射在线| 久久毛片免费看一区二区三区| 国产精品 欧美亚洲| av欧美777| 国产精品免费大片| 黑人巨大精品欧美一区二区蜜桃| 18在线观看网站| 最近手机中文字幕大全| 91精品三级在线观看| 亚洲成国产人片在线观看| 侵犯人妻中文字幕一二三四区| 欧美黑人欧美精品刺激| 亚洲成国产人片在线观看| 蜜桃国产av成人99| 亚洲欧美一区二区三区黑人| 我的亚洲天堂| 天天影视国产精品| 亚洲人成电影免费在线| 久久国产精品影院| 久久中文字幕一级| 亚洲精品久久午夜乱码| 国产精品香港三级国产av潘金莲 | 香蕉丝袜av| 中文字幕av电影在线播放| a级片在线免费高清观看视频| 国产人伦9x9x在线观看| 国产成人一区二区在线| 久久 成人 亚洲| 成人亚洲精品一区在线观看| 亚洲国产最新在线播放| 亚洲av日韩在线播放| 国产精品久久久久久精品电影小说| 精品福利永久在线观看| 国产精品国产三级国产专区5o| 国产黄频视频在线观看| 少妇粗大呻吟视频| 中文字幕人妻熟女乱码| 热re99久久精品国产66热6| 大香蕉久久成人网| 天堂俺去俺来也www色官网| 嫁个100分男人电影在线观看 | 中文字幕人妻熟女乱码| 97在线人人人人妻| 国产一区亚洲一区在线观看| 人妻 亚洲 视频| 久久久精品94久久精品| 人人妻人人爽人人添夜夜欢视频| 国产成人免费无遮挡视频| 夫妻性生交免费视频一级片| 欧美国产精品va在线观看不卡| 99国产综合亚洲精品| 欧美 亚洲 国产 日韩一| 欧美xxⅹ黑人| 一级黄色大片毛片| 日本午夜av视频| 欧美+亚洲+日韩+国产| 99热网站在线观看| avwww免费| 国产人伦9x9x在线观看| 久久鲁丝午夜福利片| 精品少妇内射三级| 少妇裸体淫交视频免费看高清 | 久久久久久久国产电影| 操美女的视频在线观看| 久久精品熟女亚洲av麻豆精品| 男女高潮啪啪啪动态图| 国产一区二区三区综合在线观看| 99久久99久久久精品蜜桃| 久久人妻熟女aⅴ| 亚洲欧美中文字幕日韩二区| 青草久久国产| 三上悠亚av全集在线观看| 亚洲伊人久久精品综合| 国产淫语在线视频| av天堂久久9| 亚洲国产日韩一区二区| 一本久久精品| 1024视频免费在线观看| 国产高清视频在线播放一区 | 久久精品国产综合久久久| 欧美xxⅹ黑人| 国产日韩欧美视频二区| 亚洲天堂av无毛| 久久久久久久国产电影| 中文字幕色久视频| 少妇被粗大的猛进出69影院| 高潮久久久久久久久久久不卡| 国产免费现黄频在线看| 在现免费观看毛片| 久久精品亚洲av国产电影网| av线在线观看网站| 日本猛色少妇xxxxx猛交久久| 欧美激情 高清一区二区三区| 国产av国产精品国产| 色网站视频免费| 亚洲一区中文字幕在线| 成年人黄色毛片网站| 男人舔女人的私密视频| 又粗又硬又长又爽又黄的视频| 国产一卡二卡三卡精品| 91精品伊人久久大香线蕉| 久久人妻熟女aⅴ| 久久国产精品影院| 又大又黄又爽视频免费| 亚洲欧美色中文字幕在线| 精品一区在线观看国产| 国产亚洲欧美在线一区二区| 老司机影院成人| 欧美日韩成人在线一区二区| 国产av精品麻豆| 亚洲国产精品一区三区| 国产成人精品在线电影| 精品福利观看| 中文字幕人妻熟女乱码| 亚洲第一av免费看| 国产高清国产精品国产三级| 女性生殖器流出的白浆| 亚洲国产精品成人久久小说| 亚洲人成电影观看| 亚洲伊人色综图| 美女扒开内裤让男人捅视频| 十八禁网站网址无遮挡| 精品久久久精品久久久| 久久精品成人免费网站| 天天躁夜夜躁狠狠躁躁| 高清av免费在线| 亚洲成av片中文字幕在线观看| 99久久综合免费| 交换朋友夫妻互换小说| 蜜桃国产av成人99| 成人国语在线视频| 黄片小视频在线播放| 大片免费播放器 马上看| cao死你这个sao货| 女人被躁到高潮嗷嗷叫费观| 大码成人一级视频| 深夜精品福利| 夜夜骑夜夜射夜夜干| 视频在线观看一区二区三区| 国产爽快片一区二区三区| 色综合欧美亚洲国产小说| 高清视频免费观看一区二区| 国语对白做爰xxxⅹ性视频网站| 久久天躁狠狠躁夜夜2o2o | 中文精品一卡2卡3卡4更新| 蜜桃在线观看..| 亚洲成人手机| 日韩av在线免费看完整版不卡| 午夜精品国产一区二区电影| 少妇的丰满在线观看| 欧美激情 高清一区二区三区| 精品卡一卡二卡四卡免费| 国产精品 国内视频| 在线亚洲精品国产二区图片欧美| 久久久久国产一级毛片高清牌| 黄片小视频在线播放| 欧美日韩国产mv在线观看视频| 日韩av在线免费看完整版不卡| 成年av动漫网址| 国产激情久久老熟女| 亚洲久久久国产精品| 国产在视频线精品| 国产色视频综合| 亚洲国产成人一精品久久久| 熟女少妇亚洲综合色aaa.| 久久鲁丝午夜福利片| 老司机亚洲免费影院| 波多野结衣一区麻豆| 日本欧美视频一区| 黄色a级毛片大全视频| 免费观看av网站的网址| 啦啦啦啦在线视频资源| 久久久久国产一级毛片高清牌| 中文字幕人妻丝袜一区二区| 久久精品亚洲熟妇少妇任你| 久9热在线精品视频| 天天躁狠狠躁夜夜躁狠狠躁| 精品一区二区三区av网在线观看 | 日韩视频在线欧美| av在线播放精品| 大型av网站在线播放| 欧美日韩精品网址| 婷婷色麻豆天堂久久| 性色av一级| 男人操女人黄网站| 欧美日韩av久久| 亚洲精品久久久久久婷婷小说| 亚洲精品第二区| 欧美日韩成人在线一区二区| 国产精品一二三区在线看| 婷婷色综合大香蕉| 自线自在国产av| 日韩制服骚丝袜av| 久久久国产欧美日韩av| 久久久久国产一级毛片高清牌| 亚洲精品久久成人aⅴ小说| 久久久久精品国产欧美久久久 | 91字幕亚洲| 亚洲精品乱久久久久久| 男人添女人高潮全过程视频| 男女高潮啪啪啪动态图| 中文字幕人妻丝袜一区二区| 久久国产亚洲av麻豆专区| 香蕉丝袜av| 亚洲欧美一区二区三区国产| 青草久久国产| 国产爽快片一区二区三区| 十分钟在线观看高清视频www| 美女午夜性视频免费| 国产高清videossex| 国产欧美日韩精品亚洲av| 午夜日韩欧美国产| 成年美女黄网站色视频大全免费| 水蜜桃什么品种好| 男人操女人黄网站| 日本五十路高清| 美女扒开内裤让男人捅视频| 18禁国产床啪视频网站| 日日夜夜操网爽| 最近中文字幕2019免费版| 亚洲精品日本国产第一区| 最新的欧美精品一区二区| 久久久久久人人人人人| 久久久久久亚洲精品国产蜜桃av| 国产精品香港三级国产av潘金莲 | 国产不卡av网站在线观看| 日日摸夜夜添夜夜爱| 亚洲欧美精品自产自拍| 热re99久久精品国产66热6| e午夜精品久久久久久久| 久久精品aⅴ一区二区三区四区| 成人18禁高潮啪啪吃奶动态图| 中文字幕人妻丝袜制服| 欧美 日韩 精品 国产| 男女高潮啪啪啪动态图| 中文字幕人妻丝袜制服| 久久久久久久大尺度免费视频| 激情五月婷婷亚洲| 黄频高清免费视频| 成年av动漫网址| 国产一区二区 视频在线| 满18在线观看网站| 久久人人97超碰香蕉20202| 久久国产精品人妻蜜桃| 久久99一区二区三区| 乱人伦中国视频| 高潮久久久久久久久久久不卡| 天天躁夜夜躁狠狠久久av| 亚洲专区中文字幕在线| 精品第一国产精品| 国产精品久久久久久人妻精品电影 | 精品一区二区三区av网在线观看 | 亚洲精品一二三| 国产老妇伦熟女老妇高清| 日韩一卡2卡3卡4卡2021年| 精品人妻熟女毛片av久久网站| 一二三四在线观看免费中文在| 一本色道久久久久久精品综合|