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

    午夜爱爱视频在线播放| 99久久精品一区二区三区| 波多野结衣巨乳人妻| av福利片在线观看| 日韩制服骚丝袜av| 18禁裸乳无遮挡免费网站照片| 久久久久性生活片| 人妻 亚洲 视频| 熟女电影av网| 亚洲av男天堂| 欧美xxxx黑人xx丫x性爽| 嘟嘟电影网在线观看| 夜夜爽夜夜爽视频| 身体一侧抽搐| 久久久久久伊人网av| 色视频在线一区二区三区| 亚洲欧美日韩另类电影网站 | 久久久色成人| 日韩视频在线欧美| 身体一侧抽搐| 男女国产视频网站| 国产黄色视频一区二区在线观看| 国产精品麻豆人妻色哟哟久久| www.av在线官网国产| 国产真实伦视频高清在线观看| 欧美另类一区| 国产中年淑女户外野战色| 99热这里只有是精品50| av在线蜜桃| 久久久成人免费电影| 校园人妻丝袜中文字幕| 亚洲精品国产av蜜桃| 久久久久久久久久成人| 精品久久久久久久久亚洲| 青春草亚洲视频在线观看| 国产片特级美女逼逼视频| 亚洲人成网站高清观看| 亚洲第一区二区三区不卡| av免费在线看不卡| 国内精品宾馆在线| 免费人成在线观看视频色| 最近中文字幕2019免费版| 少妇丰满av| 国产精品女同一区二区软件| 国产69精品久久久久777片| 亚洲精品一二三| 亚洲最大成人中文| 亚洲精品久久久久久婷婷小说| 国产 一区精品| 国产精品国产三级国产av玫瑰| 精品久久久久久久末码| 99热这里只有是精品50| 在线看a的网站| 亚洲精品乱久久久久久| 国内精品宾馆在线| 国产又色又爽无遮挡免| 精品人妻视频免费看| 青春草亚洲视频在线观看| 久久精品国产亚洲网站| 国产成人一区二区在线| 在线播放无遮挡| freevideosex欧美| 在线观看国产h片| 97精品久久久久久久久久精品| 久久久a久久爽久久v久久| 国产av码专区亚洲av| 菩萨蛮人人尽说江南好唐韦庄| 男男h啪啪无遮挡| 免费黄网站久久成人精品| 国产成人精品福利久久| 一区二区av电影网| 婷婷色综合www| 国产综合懂色| 纵有疾风起免费观看全集完整版| 国国产精品蜜臀av免费| 婷婷色av中文字幕| www.av在线官网国产| 国内揄拍国产精品人妻在线| 亚洲精品国产av成人精品| 乱码一卡2卡4卡精品| 三级男女做爰猛烈吃奶摸视频| 国产精品秋霞免费鲁丝片| 国产免费一区二区三区四区乱码| 亚洲三级黄色毛片| 久久久亚洲精品成人影院| 久久久久久国产a免费观看| 一级黄片播放器| 日韩在线高清观看一区二区三区| 婷婷色综合大香蕉| 高清欧美精品videossex| 又黄又爽又刺激的免费视频.| 美女脱内裤让男人舔精品视频| 亚洲综合色惰| 久久久久九九精品影院| 亚洲人成网站高清观看| 国内少妇人妻偷人精品xxx网站| 国产国拍精品亚洲av在线观看| 久久国内精品自在自线图片| 国产午夜精品一二区理论片| 日产精品乱码卡一卡2卡三| 搡老乐熟女国产| 国产视频内射| 白带黄色成豆腐渣| 色吧在线观看| 国产精品99久久99久久久不卡 | 97超视频在线观看视频| 国产精品99久久99久久久不卡 | av在线天堂中文字幕| 啦啦啦啦在线视频资源| 日本熟妇午夜| 97超视频在线观看视频| 国产视频内射| 大陆偷拍与自拍| 国产精品久久久久久精品电影| 日韩av在线免费看完整版不卡| 欧美日韩亚洲高清精品| 日韩视频在线欧美| 免费观看av网站的网址| 精品国产三级普通话版| 在线观看免费高清a一片| 麻豆精品久久久久久蜜桃| 国产亚洲av片在线观看秒播厂| 伦精品一区二区三区| 天天躁夜夜躁狠狠久久av| 777米奇影视久久| 国产在视频线精品| 亚洲久久久久久中文字幕| 久久久久久久久久成人| 少妇人妻一区二区三区视频| 男女边吃奶边做爰视频| 一区二区三区乱码不卡18| 国产精品蜜桃在线观看| 精品视频人人做人人爽| 免费av不卡在线播放| 久久99精品国语久久久| 少妇的逼好多水| 精品少妇久久久久久888优播| av卡一久久| 久久久久久久大尺度免费视频| 久久人人爽人人片av| 日韩一区二区三区影片| 亚洲性久久影院| 国产精品一区二区性色av| 亚洲国产av新网站| 精品久久久久久久末码| 欧美丝袜亚洲另类| 极品少妇高潮喷水抽搐| 不卡视频在线观看欧美| 国产亚洲5aaaaa淫片| 亚洲国产欧美在线一区| 18禁在线播放成人免费| 激情 狠狠 欧美| 一级片'在线观看视频| 久久午夜福利片| 一级黄片播放器| 久久久久久久亚洲中文字幕| 男人添女人高潮全过程视频| 亚洲色图av天堂| 国产高潮美女av| 中文欧美无线码| 亚洲自偷自拍三级| 国产av国产精品国产| 亚洲成人久久爱视频| 在线观看人妻少妇| 特级一级黄色大片| 高清日韩中文字幕在线| 久久精品夜色国产| 精品酒店卫生间| 五月玫瑰六月丁香| 午夜精品国产一区二区电影 | 国产精品国产三级国产av玫瑰| 国产黄色视频一区二区在线观看| .国产精品久久| 午夜激情久久久久久久| 亚洲国产成人一精品久久久| 久久热精品热| 一个人观看的视频www高清免费观看| 欧美一级a爱片免费观看看| 2022亚洲国产成人精品| 大香蕉97超碰在线| 伊人久久精品亚洲午夜| 男女边吃奶边做爰视频| 99久久人妻综合| 三级国产精品欧美在线观看| 亚洲,一卡二卡三卡| 日本一本二区三区精品| 黄片wwwwww| 高清欧美精品videossex| 日韩成人伦理影院| 亚州av有码| 99热全是精品| 天堂中文最新版在线下载 | 免费少妇av软件| 日产精品乱码卡一卡2卡三| 美女内射精品一级片tv| 18+在线观看网站| 久久久久精品久久久久真实原创| 亚洲精品日本国产第一区| 国产毛片在线视频| 777米奇影视久久| 一级av片app| 美女脱内裤让男人舔精品视频| 91久久精品国产一区二区三区| 亚洲国产成人一精品久久久| 婷婷色麻豆天堂久久| 在线观看免费高清a一片| 亚洲丝袜综合中文字幕| 最新中文字幕久久久久| 久久人人爽人人片av| 亚洲成色77777| 国产成人福利小说| 成年av动漫网址| 日韩一本色道免费dvd| 国产精品久久久久久av不卡| 在线播放无遮挡| 国产成人91sexporn| 狂野欧美白嫩少妇大欣赏| 亚洲av二区三区四区| 亚洲欧洲日产国产| 丝袜喷水一区| 日韩不卡一区二区三区视频在线| 三级国产精品片| 老司机影院毛片| 激情 狠狠 欧美| 国产亚洲av片在线观看秒播厂| 2021天堂中文幕一二区在线观| 97热精品久久久久久| 国产精品不卡视频一区二区| 五月玫瑰六月丁香| 高清在线视频一区二区三区| 国产午夜精品久久久久久一区二区三区| 国模一区二区三区四区视频| 国产欧美亚洲国产| 一边亲一边摸免费视频| 国产乱人偷精品视频| 亚洲av国产av综合av卡| 国产精品熟女久久久久浪| 国产探花极品一区二区| 欧美少妇被猛烈插入视频| 亚洲不卡免费看| 亚洲激情五月婷婷啪啪| 毛片女人毛片| 国产黄片美女视频| 激情 狠狠 欧美| 真实男女啪啪啪动态图| 国产精品秋霞免费鲁丝片| 人妻系列 视频| 九九在线视频观看精品| 国产精品女同一区二区软件| 午夜免费观看性视频| 精品国产三级普通话版| 国产v大片淫在线免费观看| 亚洲,欧美,日韩| 亚洲av.av天堂| 成人鲁丝片一二三区免费| 午夜福利在线观看免费完整高清在| 日本av手机在线免费观看| 国产欧美日韩精品一区二区| 亚洲精品乱久久久久久| 高清在线视频一区二区三区| 最新中文字幕久久久久| 天堂俺去俺来也www色官网| 免费观看在线日韩| 丰满少妇做爰视频| 国产白丝娇喘喷水9色精品| 丰满乱子伦码专区| 中文字幕久久专区| 最近最新中文字幕免费大全7| 热re99久久精品国产66热6| 日本黄大片高清| 国精品久久久久久国模美| 波多野结衣巨乳人妻| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 欧美区成人在线视频| 日本av手机在线免费观看| 久久久精品免费免费高清| 国产精品一及| 亚洲av在线观看美女高潮| 自拍偷自拍亚洲精品老妇| 国产亚洲av嫩草精品影院| 伦理电影大哥的女人| 国产精品精品国产色婷婷| 中文字幕久久专区| 大码成人一级视频| 国产一区有黄有色的免费视频| 免费观看a级毛片全部| 日韩 亚洲 欧美在线| 精品国产露脸久久av麻豆| 视频中文字幕在线观看| 午夜免费男女啪啪视频观看| 国产精品三级大全| 免费人成在线观看视频色| 99热网站在线观看| 欧美高清性xxxxhd video| 日本与韩国留学比较| 五月天丁香电影| 人人妻人人澡人人爽人人夜夜| 国产精品女同一区二区软件| 尾随美女入室| 国产视频首页在线观看| 国产伦精品一区二区三区四那| www.色视频.com| 欧美日韩亚洲高清精品| av免费观看日本| 人妻系列 视频| 亚洲av欧美aⅴ国产| 99久久精品国产国产毛片| 五月开心婷婷网| 久久人人爽av亚洲精品天堂 | 久久精品久久久久久久性| 边亲边吃奶的免费视频| 日韩成人av中文字幕在线观看| 在线观看人妻少妇| 国产精品成人在线| 精品一区二区三区视频在线| 一二三四中文在线观看免费高清| 中文资源天堂在线| 黑人高潮一二区| 久久精品国产亚洲网站| 国产有黄有色有爽视频| 性插视频无遮挡在线免费观看| 99久久精品一区二区三区| 亚洲综合色惰| 亚洲不卡免费看| 91aial.com中文字幕在线观看| 男插女下体视频免费在线播放| 国产伦理片在线播放av一区| 国产亚洲91精品色在线| tube8黄色片| 美女主播在线视频| 国产69精品久久久久777片| 交换朋友夫妻互换小说| 大香蕉97超碰在线| 夫妻性生交免费视频一级片| 少妇猛男粗大的猛烈进出视频 | 国产成人免费无遮挡视频| 国产又色又爽无遮挡免| 国产一区二区亚洲精品在线观看| 免费少妇av软件| 99久久精品热视频| 成人亚洲精品av一区二区| 亚洲精品久久久久久婷婷小说| 亚洲四区av| 99久久九九国产精品国产免费| 亚洲av一区综合| 免费人成在线观看视频色| 我的女老师完整版在线观看| 国产欧美日韩一区二区三区在线 | 人妻系列 视频| 三级经典国产精品| 有码 亚洲区| 国产乱人偷精品视频| 亚洲av中文av极速乱| 日本av手机在线免费观看| 亚洲丝袜综合中文字幕| 久久久国产一区二区| 亚洲av成人精品一区久久| 少妇人妻久久综合中文| 高清毛片免费看| 日本午夜av视频| 日本黄色片子视频| 免费高清在线观看视频在线观看| 亚洲欧美日韩另类电影网站 | 汤姆久久久久久久影院中文字幕| 久久久久国产精品人妻一区二区| 国产永久视频网站| 国产人妻一区二区三区在| 免费观看在线日韩| 毛片女人毛片| 寂寞人妻少妇视频99o| 一级毛片我不卡| 亚洲经典国产精华液单| 欧美成人a在线观看| 日韩,欧美,国产一区二区三区| 听说在线观看完整版免费高清| 成人午夜精彩视频在线观看| 欧美另类一区| 人人妻人人看人人澡| 国产精品久久久久久精品电影小说 | 男插女下体视频免费在线播放| 黄色日韩在线| 一区二区三区免费毛片| 久久女婷五月综合色啪小说 | 亚洲人成网站在线观看播放| 国产男女内射视频| 午夜免费观看性视频| 中文精品一卡2卡3卡4更新| 免费黄频网站在线观看国产| 免费av毛片视频| 蜜桃久久精品国产亚洲av| 欧美日韩精品成人综合77777| 久久精品久久久久久久性| 香蕉精品网在线| 秋霞伦理黄片| 久久久久久久久久人人人人人人| 日韩精品有码人妻一区| 欧美日韩视频高清一区二区三区二| 亚洲欧美日韩卡通动漫| 国产91av在线免费观看| 欧美老熟妇乱子伦牲交| 大话2 男鬼变身卡| 欧美日韩一区二区视频在线观看视频在线 | 日韩 亚洲 欧美在线| 免费黄网站久久成人精品| 日韩视频在线欧美| av.在线天堂| 人人妻人人爽人人添夜夜欢视频 | 国产精品一及| 国产男人的电影天堂91| 国产精品久久久久久av不卡| 国产成人aa在线观看| 蜜桃久久精品国产亚洲av| 久久久久久久久久久免费av| 街头女战士在线观看网站| 精品一区在线观看国产| 成人鲁丝片一二三区免费| 亚洲精品久久久久久婷婷小说| 建设人人有责人人尽责人人享有的 | 国产伦理片在线播放av一区| 日本爱情动作片www.在线观看| 亚洲,欧美,日韩| 欧美日韩视频高清一区二区三区二| 搞女人的毛片| 日日摸夜夜添夜夜爱| 大香蕉久久网| 国产成人a∨麻豆精品| 亚洲色图av天堂| 一个人观看的视频www高清免费观看| 精品酒店卫生间| 女人十人毛片免费观看3o分钟| 黄色配什么色好看| 天天躁日日操中文字幕| 欧美区成人在线视频| 一区二区三区乱码不卡18| 国产有黄有色有爽视频| 午夜爱爱视频在线播放| 国产黄频视频在线观看| 极品教师在线视频| 亚洲av男天堂| 夜夜爽夜夜爽视频| 久久久精品94久久精品| 十八禁网站网址无遮挡 | 色哟哟·www| 少妇的逼水好多| 久久久精品94久久精品| 久久精品国产亚洲av涩爱| 精品少妇久久久久久888优播| 人妻制服诱惑在线中文字幕| 一级毛片我不卡| 国内精品美女久久久久久| 亚洲国产精品成人久久小说| 亚洲精品,欧美精品| 婷婷色av中文字幕| 日韩制服骚丝袜av| 亚洲av成人精品一区久久| 麻豆国产97在线/欧美| 国产精品不卡视频一区二区| 日本av手机在线免费观看| 国产白丝娇喘喷水9色精品| 久久久久九九精品影院| 亚洲精品aⅴ在线观看| 久久精品夜色国产| 91aial.com中文字幕在线观看| 国产91av在线免费观看| 国产一区亚洲一区在线观看| 久久久精品免费免费高清| 欧美 日韩 精品 国产| 五月开心婷婷网| 国产成人精品婷婷| 成人鲁丝片一二三区免费| 亚洲,欧美,日韩| 国产av不卡久久| 亚洲人成网站在线播| 欧美日韩综合久久久久久| 亚洲av福利一区| 禁无遮挡网站| 欧美日本视频| 国产亚洲5aaaaa淫片| 噜噜噜噜噜久久久久久91| 99热这里只有精品一区| 天堂中文最新版在线下载 | 丝袜脚勾引网站| 色网站视频免费| 国国产精品蜜臀av免费| 国产乱人偷精品视频| 国产午夜精品久久久久久一区二区三区| 18禁在线无遮挡免费观看视频| 伊人久久精品亚洲午夜| 中国国产av一级| 人人妻人人看人人澡| 国国产精品蜜臀av免费| 男女无遮挡免费网站观看| 中国三级夫妇交换| 人人妻人人看人人澡| 女人十人毛片免费观看3o分钟| 你懂的网址亚洲精品在线观看| 亚洲国产精品成人久久小说| 成人黄色视频免费在线看| 别揉我奶头 嗯啊视频| 爱豆传媒免费全集在线观看| 免费av观看视频| 一区二区三区四区激情视频| 欧美xxⅹ黑人| 成人鲁丝片一二三区免费| 亚洲伊人久久精品综合| 成人欧美大片| 欧美激情在线99| 全区人妻精品视频| 久久精品夜色国产| 日韩不卡一区二区三区视频在线| 精品视频人人做人人爽| 久久99热这里只频精品6学生| 三级男女做爰猛烈吃奶摸视频| 午夜爱爱视频在线播放| 3wmmmm亚洲av在线观看| 男女那种视频在线观看| av一本久久久久| 亚洲欧美日韩东京热| 一级毛片我不卡| 亚洲一区二区三区欧美精品 | 777米奇影视久久| 亚洲国产精品国产精品| 色哟哟·www| 免费大片黄手机在线观看| 亚洲va在线va天堂va国产| 久久热精品热| 免费黄网站久久成人精品| 高清av免费在线| 欧美人与善性xxx| 国产在视频线精品| 婷婷色麻豆天堂久久| 丝瓜视频免费看黄片| 国产亚洲午夜精品一区二区久久 | 99久国产av精品国产电影| 神马国产精品三级电影在线观看| 熟妇人妻不卡中文字幕| 高清av免费在线| 日韩,欧美,国产一区二区三区| 日本欧美国产在线视频| 国产精品偷伦视频观看了| 久久精品国产a三级三级三级| 国国产精品蜜臀av免费| 国产白丝娇喘喷水9色精品| 国产人妻一区二区三区在| 看非洲黑人一级黄片| 国产精品麻豆人妻色哟哟久久| 亚洲最大成人av| 国产黄片美女视频| 丝袜脚勾引网站| 男女下面进入的视频免费午夜| 国产成人精品一,二区| 久久人人爽av亚洲精品天堂 | 一本色道久久久久久精品综合| 午夜视频国产福利| 国产免费一区二区三区四区乱码| 神马国产精品三级电影在线观看| av黄色大香蕉| 亚洲最大成人av| 人体艺术视频欧美日本| 免费黄色在线免费观看| 草草在线视频免费看| 又黄又爽又刺激的免费视频.| 亚洲av电影在线观看一区二区三区 | 久久影院123| 国产精品一区二区在线观看99| 丝袜脚勾引网站| 日韩不卡一区二区三区视频在线| 亚洲人成网站在线观看播放| 97热精品久久久久久| 亚洲av中文av极速乱| 新久久久久国产一级毛片| av一本久久久久| 国产老妇女一区| 久久精品久久久久久久性| 国产一区亚洲一区在线观看| 亚洲欧美日韩东京热| 国产人妻一区二区三区在| 亚洲欧洲国产日韩| 国产伦精品一区二区三区四那| 国产成人免费无遮挡视频| 91午夜精品亚洲一区二区三区| 内地一区二区视频在线| 国产一区二区亚洲精品在线观看| 夫妻性生交免费视频一级片| 建设人人有责人人尽责人人享有的 | 日本与韩国留学比较| 久久久久久久亚洲中文字幕| 少妇裸体淫交视频免费看高清| 精品国产一区二区三区久久久樱花 | 精品久久久久久久末码| 亚洲激情五月婷婷啪啪| 直男gayav资源| 国产精品爽爽va在线观看网站| 插逼视频在线观看| 欧美丝袜亚洲另类| 亚洲一级一片aⅴ在线观看| 久久亚洲国产成人精品v| 插阴视频在线观看视频| av播播在线观看一区| 国产高清三级在线| 超碰av人人做人人爽久久| 两个人的视频大全免费| 亚洲精品中文字幕在线视频 | 国产午夜福利久久久久久| 99热全是精品| 国产欧美另类精品又又久久亚洲欧美| 欧美激情国产日韩精品一区| 亚洲av免费高清在线观看| av国产免费在线观看| 久久久久久久久久久丰满| 精品午夜福利在线看| 国产精品不卡视频一区二区| 永久免费av网站大全| 国产成人福利小说| 久久亚洲国产成人精品v|