999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

Thermodynamic analysis and modification of Gibbs–Thomson equation for melting point depression of metal nanoparticles

2021-05-18 11:06:48NanhuaWuXiaohuaLuRongAnXiaoyanJi

Nanhua Wu,Xiaohua Lu*,Rong An,Xiaoyan Ji*

1 Division of Energy Science/Energy Engineering,Lule? University of Technology,97187 Lule?,Sweden

2 State Key Laboratory of Materials-Oriented and Chemical Engineering,Nanjing Tech University,Nanjing 211816,China

3 Key Laboratory of Advanced Catalytic Materials and Technology,School of Petrochemical Engineering,Changzhou University,Changzhou 213164,China

4 Herbert Gleiter Institute of Nanoscience,Nanjing University of Science &Technology,Nanjing 210094,China

Keywords:Melting point depression Metal nanoparticle Gibbs–Thomson equation Substrate Interfacial effect

ABSTRACT Abnormal melting point depression of metal nanoparticles often occurs in heterogeneous catalytic reactions,which leads to a reduction in the stability of reactive nanoclusters.To study this abnormal phenomenon,the original and surface-energy modified Gibbs–Thomson equations were analyzed in this work and further modified by considering the effect of the substrate.The results revealed that the original Gibbs–Thomson equation was not suitable for the particles with radii smaller than 10 nm.Moreover,the performance of the surface-energy modified Gibbs–Thomson equation was improved,and the deviation was reduced to(-350-100)K,although further modification of the equation by considering the interfacial effect was necessary for the small particles (r <5 nm).The new model with the interfacial effect improved the model performance with a deviation of approximately -50 to 20 K,where the interfacial effect can be predicted quantitatively from the thermodynamic properties of the metal and substrate.Additionally,the micro-wetting parameter αw can be used to qualitatively study the overall impact of the substrate on the melting point depression.

1.Introduction

Nanoparticles have been widely used in academia and industry,and their novel physical and chemical properties are different from those in the bulk phase[1].Due to the special nano-size effect and specific surface area,nanoparticle-supported metal catalysts play an important role [2–4].Nanoparticle-supported metal catalysts have shown great potential in the production of clean fuels,chemicals,and pharmaceuticals [5,6].However,maintaining the stability of the catalysts is always a key issue for practical applications[3].The agglomeration due to particle surface melting and then particle growth via particle diffusion and coalescence or Ostwald ripening [5,7]was considered as the main reason leading to catalyst deactivation.For instance,in a low-temperature (below 200 K) CO oxidation reaction,supported Au catalysts show high activity;however,the activity declines significantly when the temperature increases to room temperature [8]because of surface melting and the increased mobility of particles at high temperatures [9].In addition,it was suggested that both the particle size and the interfacial contacted environment(i.e.,the support materials) are important for controlling stability and activity of the catalyst.Moreover,the melting point depression of the metal used as the catalyst is an intrinsic factor that controls the catalyst stability.Thus,it is necessary to predict to the melting point of metals using an efficient method,for evaluating the catalyst stability.

The melting points (or melting temperatures) of metals have been studied both experimentally and theoretically [10–16].The results showed that the melting point decreased with the a decrease in particle size,and it decreased significantly when the particle size was reduced to nanoscale.In theory,Gibbs,for the first time,applied thermodynamics to study the effects of geometry on properties such as vapor pressure and melting point,which is the well-known Kevin equation [17–19].In particular,due to the increasingly vital role of nanotechnology and nanoscience,a semi-empirical engineering model with the consideration of sizeeffect,the Gibbs–Thomson equation,has been widely used to calculate the melting point of metals in the fields of bio-materials and nano-materials [20].The Gibbs–Thomson equation was originally proposed by Jackson and McKenna et al.[21],in which the surface energy was included in the total Gibbs energy and considered as the capillary effect due to the curvature of the interfacial surface.The Gibbs–Thomson equation can also be used to calculate the solid and gas solubilities for various particle geometries.

Recently,the Gibbs–Thomson equation was further modified by Lee et al.considering the surface-energy (i.e.,the modified Gibbs–Thomson equation) with a significant improvement for small particles [15].However,similar to the original Gibbs–Thomson equation,the interfacial effect has not been considered.Meanwhile,the interfacial effect on the melting point depression has been studied theoretically [22–28].For example,Lee et al.developed a theoretical method based on CALPHAD to reveal both size and interfacial effects on the chemical potential of nanoparticles[22].A number of other studies used the curvature diameter and contact angle to describe the interfacial effect [23,26].However,none of them pointed out the particle size range in which the interfacial effect is dominated.At the same time,the application of these models is generally complicated because of the difficulty in obtaining the model parameters.

Recently,on the basis of molecular simulations,Gubbins and his co-workers [29–32]indicated that when the pore size reduces to the nano-scale (pore diameter less than 15 nm),the interfacial behavior affects the fluid properties;specifically,(1)the molecular interaction between the particles and the contacted wall as well as the particle–wall molecular diameter,are two key factors affecting the melting points of particles confined in pores;and(2)the effect of the contacted wall on interfacial properties relates to the wettability (wetting parameter αw).This work provides a parameter to relate the melting point depression of fluids confined in nano-slit pores with molecular interactions.However,it is still unclear how this parameter can be used to describe the melting point depression of metals when the particle size is on the nanoscale.

This study aimed to investigate the melting point depression of metals using the engineering model,i.e.,the Gibbs–Thomson equation-based models.The performance of the original and modified Gibbs–Thomson equations was analyzed.A further modification by considering the interfacial effect was proposed to calculate the melting point depression for nanosized particles,and the performance was analyzed to illustrate the effect of the modification in different particle-size levels.Finally,the relationship between the wetting parameter αwand the adjustable parameter of the model that was related to the interfacial effect was further discussed.

2.Available Data

The experimental data are important for verifying the model performance or developing a theoretical model.In this work,the size-dependent melting points of metals by experiment were surveyed and are summarized in Table 1.As the effect of the contacted substrate on the melting point of metal was considered in this work,the substrates used for experimental measurements are also listed in Table 1.There are other three sets from molecular simulations(Pt,Ni,Mg),the information of substrates cannot be obtained.

The literature survey showed that for each metal,only one set of data is available,and evaluation of the available data is infeasible.Due to the limited availability of data,all data sets were used for the investigation conducted in this work.

3.Model Analysis

Experimental research has shown that the melting point of metals decreases with a decrease in particle size,and it decreases significantly when the particle size reduces to a certain small scale.Although the original and modified Gibbs–Thomson equations have been proposed to study the melting point of metals,the exact particle size for which a particular theoretical model can be used to provide reliable results remains unclear.Our anlysis was conducted in this regard.

3.1.Original Gibbs-Thomson equation

The total Gibbs free energy of a particle(Gs)is the summation of the Gibbs free energy in the bulk (Gs,bulk) and that of the surface(Gs,interface) [33],i.e.,

When the metal particle is large,the contribution of the surface part is relatively small compared to the bulk contribution.However,for the relatively small particles,the contribution of Gs,interfaceis significant.At equilibrium,the Gibbs free energy in the liquid phase(Gs)is the same as that in the solid phase(Gs),including the contributions of the bulk and interface terms in both cases.A schematic graph illustrating a nanoscale metal particle on a substrate is presented in Scheme 1.

Thus,we can write,

Based on classical thermodynamics and assuming that the fusion enthalpy is constant,ΔGs-l,bulklinked to the melting point of the bulk phase can be expressed,

where ΔHmis the fusion enthalpy in the bulk phase,ΔSmis the entropy in the bulk phase,Tmis the bulk melting point of the particle,T(r) is the melting point of the particle with radius r.

For a spherical particle,the contribution of ΔGs-l,interfacederived from the Kelvin equation [33]is as follows:

where γslis the surface tension for the solid–liquid interface,Vs,mmeans the molar volume of the solid particle,and r is the radius of the spherical solid particle.As a result,for a small isolated spherical particle with a radius r,the melting point depression ΔTmcan be calculated using the Gibbs–Thomson equation proposed by Jackson and McKenna [21]:

In the original Gibbs–Thomson equation,the surface tension γslwas considered to be size-independent,and the effect of the contacted environment was negligible.This is a reasonable assumption when the particle size is not very small.

Table 1 Experimental data for the melting points of nanoscale mental particles on substrates

Scheme 1.The description of supported metal nanoparticles.

3.2.Modified Gibbs–Thomson equation

However,with the decreases in particle size,the interfacial properties become increasingly important,and the assumption on γslneeds to be re-considered.Molecular simulations have shown that the surface energy of solid particles increases with a decrease in size,while that of liquid remains unchanged[15];thus,the effect of size on the surface energy(γsl)needs to be considered.

The effect of surface energy on γslhas been considered by Lee et al.[15]

where reis the first nearest neighbor distance for crystalline materials.

Thus,the modified Gibbs–Thomson equation (i.e.,Lee’s model)can be written as,

3.3.Properties

The values of size-independent Tm,bulk,γsl,ΔHmand Vsare necessary to calculate the melting point depression.For the metals listed in Table 1,their corresponding properties of Tm,bulk[34],ΔHm[34],Vs[35]and re[36]obtained from the literature are listed in Table 2.

The size-independent surface tension γslcan be determined experimentally or obtained from the fitting of the experimental melting point depression using Eq.(5)in the valid size range(large size scale).Experimentally,γslcan be measured by different methods,such as dihedral angles (DA),contact angles (CA),and the shape of the grain-boundary-grooves(GBG)[37].Among the three methods,DA has been commonly used to measure the surface tension for most metals,except Mg.The surface tensions measured with different methods are listed in Table 3 and are denoted as γsl,exp.It can be seen that the surface tensions varies with the measurement method used.For example,the surface tensions of Au were measured by the both DA and CA methods,and the corresponding measured results were 0.156 and 0.190 J·m-2,respectively.The surface tensions of Pb measured by the DA and GBG methods were 0.076 and 0.054 J·m-2,respectively.It is worth noting that even when using the same experimental method,the measured surface tensions showed a large variation (~30%) for a same-size metal.For example,for Bi,the surface tensions measured by the commonly used DA method were 0.061,0.074,and 0.082 J·m-2,respectively [37].Therefore,additional experiments need to be carried out to measure the surface tension and thus obtain more reliable results.

γslcan also be obtained from the fitting of the melting point depression using Eq.(5) with the known properties of Tm,bulk,ΔHmand Vslisted in Table 2.Following the method described in reference [38],the parameters γsl(γsl,fit) for different metals were fitted.The results are listed in Table 3.For most metals,except Bi,γsl,fitcan be obtained by fitting with high accuracy,i.e.,R2>0.99.

Table 2 Physical properties of metals

The comparison of γsl,fitwith γsl,expin Table 3 shows that the fitted and experimental values agree with each other for most metals,except Bi.γsl,fitof Bi is 0.196 J·m-2,which is much larger than the experimental values (0.061,0.074,or 0.082 J·m-2).To further study the effect of γslon the performance of the Gibbs–Thomson equation,the melting points of Bi were calculated using γsl,exp(=0.082 J·m-2) and γsl,fit(=0.196 J·m-2),respectively.The results show that the melting points differ significantly when the particle size decreases to 19 nm(i.e.,the calculated melting points at r=19 nm are 537 and 525 K,respectively).The discrepancy increases with the decreases in particle size,for example,at a particle size of 3 nm,the calculated melting points are 490 and 411 K,respectively,with γsl,exp(=0.082 J·m-2) and γsl,fit(=0.196 J·m-2).Because of the lack of experimental data from other sources and the lack of the available data for particles with relatively larger sizes(40–80 nm),it is very difficult to determine the more accurate γslvalue.

3.4.Analysis

To illustrate the performance of the model and clarify the model prediction capacity,the model results were compared with the available experimental data in this section.Analysis was carried out to show the capability of the model with or without modifications.

3.4.1.Original Gibbs-Thomson equation

The original Gibbs–Thomson equation was used to calculate the melting point depression for all metals listed in Table 1 with the physical properties of Tm,bulk,ΔHm,and Vslisted in Table 2.Because the fitted γsl(γsl,fit) listed in Table 3 is reasonable compared to experimental γsl,expfor most metals,γsl,fitwas used in the calculation.The model results were compared with the available data to evaluate the performance of the Gibbs–Thomson equation.The comparison is illustrated in Fig.1a.

As shown in Fig.1a,the original Gibbs–Thomson equation with size-independent properties (i.e.,constant properties listed in Table 2 and γsl,fitin Table 3) is suitable for calculating the melting point depression when the particle size is relatively large.When the particle size r is larger than 10 nm(r-1<0.1 nm-1),the results of the Gibbs–Thomson equation agree well with the data available in the literature,and the deviation is within ±10 K.When the particle size further decreases down to 5 nm(0.1 nm-1<r-1<0.20-nm-1),the results obtained using the original Gibbs–Thomson equation show some discrepancies,and the deviation increases to ±20 K.

With a further decrease in particle size (r-1>0.2 nm-1,or r <5 nm),the deviation in the results for the original Gibbs–Thomson equation increases significantly.The deviations are positive for Pt,Au,Mg,and Ni,which means that the original Gibbs–Thomson equation gives an overestimation.In particular for Pt,the deviation is approximately 700 K for the given smallest size (0.76 nm,i.e.,r-1=1.31 nm-1).For Pb,In,and Al,the originalGibbs–Thomson equation leads to an underestimation.In particular,for Al,the estimated melting point is about 50 K lower than the experimental value when the particle size decreases to 2.05 nm(r-1=0.49 nm-1).Therefore,we can conclude that for the particles with a small size(r <5 nm),the original Gibbs–Thomson equation with size-independent properties listed in Table 2 and γsl,fitin Table 3 is not suitable,and modification is needed.

Table 3 Interfacial tension of metals

3.4.2.Modified Gibbs-Thomson equation

The melting point depression for all metals listed in Table 1 was calculated using Eq.(7),i.e.,the modified Gibbs–Thomson equation or Lee’s model,where the surface-energy was modified.In the calculation,again,the size-independent physical properties listed in Table 2,γsl,fittaken from Table 3.The model results were compared with the data available in the literature.The absolute deviation between the model results and the available data is shown in Fig.1b.The deviation in the results of the modified Gibbs–Thomson equation is much smaller than that illustrated in Fig.1a.For the original Gibbs–Thomson equation,the absolute deviation is(-100-700) K,while it decreases to (-350-100) K for the modified equation.Therefore,the overall model performance is remarkably improved due to the surface-energy modification.

The comparison is further illustrated in Fig.2,in which the vertical axis represents (Tm-T (r)),i.e.,the left-hand side of Eq.(7),with the data of Tmand T(r)from the literature,and the horizontal axis represents the right-hand side of Eq.(7) calculated by substituting the literature values for the properties listed in Tables 2 and 3.For Pt,the values are located on the cross line (situation 1 with slope=1),which implies that the model with the surface-energy modification can be used to calculate the melting point depression for this metal.For other metals,all the values are below the cross line (situation 2 with slope <1),which implies underestimation,i.e.,the abnormal melting point depression compared to the modification with the surface-energy only.

3.5.New modification

In this work,the effect of the contacted environment (i.e.,the interfacial effect,which refers to the substrate in this work)is considered using an adjustable parameter δ,i.e.,

Fig.1.Deviation of the melting points of metals calculated with the original Gibbs–Thomson equation (a) and modified Gibbs–Thomson equation (surface-energy modification)(b) from literature data.■:Pt;▼:Au;▲:Mg;●:Ni;?:Pb;★:In;◆:Bi;?:Al.

Fig.2.Deviation of melting points of metals calculated with the modified Gibbs–Thomson equation(surface-energy modification)from the literature data.■:Pt;▼:Au;▲:Mg;●:Ni;?:Pb;★:In;◆:Bi;?:Al.

Using the available data on the melting point depression of metals together with the properties listed in Table 2 and γsl,fitlisted in Table 3,the parameter δ was obtained from the fitting based on Eq.(8).The fitted parameter δ and the corresponding average relative deviation (ARD) are listed in Table 4.From the values listed in Table 4,we can see that for Pt,δ is~0,which means that the model with the surface-energy modification (i.e.,Lee’s model) can provide reliable results.The δ value for Au is much smaller than those for other metals,which means that the interfacial effect is not the dominant part.For other six metals,a large value of δ implies that the interfacial effect becomes dominant when the particle size is on the nanoscale.

3.5.1.Discussion on the interfacial effect

Although the model with the surface-energy modification showed improved the performance,it was not sufficiently accurate in the small particle range for some metals,as shown in Figs.1b and 2.That is the equation with surface-energy modification cannot be used to calculate the abnormal melting point depression of metals.During the experimental measurements,the particle was supported on a surface of the substrate.When the particle size is sufficiently small,the substrate might have a significant effect on the properties of the molecules on the surface of the solid particles.If the interaction between the molecules of the substrate and the metal is weak,the effect of the substrate can be neglected,and no abnormal melting point depression will be observed.This corresponds to ‘‘situation 1(slope=1)”in Fig.2 for Pt metal.If the interaction between the molecules of the substrate and the metal is strong,the effect of the substrate is crucial and an abnormal melting depression will be observed.This corresponds to ‘‘situation 2 (slope <1)”in Fig.2 for metals other than Pt.

In this work,an empirical parameter δ was used to consider the interfacial effect.Recently,Gubbins and his co-workers[29,30]proposed another empirical parameter αwto study the interfacial influence(i.e.,the contacted wall)on the melting/freezing behavior of fluids in nanoporous media.Subsequently,furtheranalysis was carried out to investigate the relationship between δ and αw.

Table 4 The fitted δ and the corresponding average relative deviation with new modified Gibbs–Thomson equation (Eq.(8))

(1) Theoretical analysis based on Gubbins’ work

Recently,Gubbins and his co-workers[29,30]studied the melting/freezing behavior of fluids in nanoporous media to describe the interfacial influence on fluid behaviors,and they proposed a wetting parameter αwthat is related to the molecular interaction and size,i.e.,

where ρwis the density of wall atoms,Δ is the interlayer spacing in the solid,σfwis the fluid-wall diameter parameter,and εfwand εffare the Lennard-Jones potentials for fluid-wall and fluid–fluid molecules,respectively.

The molecular simulation results for the fluids in slit-shaped pores in the work by Gubbins and his co-workers [32]indicated the following:

(1) A small αwimplies non-ideal wettability and the interfacial influence could be neglected.

(2) A large αwimplies ideal wettability and the interfacial influence plays an important role in estimating the melting points.

In this work,αwwas used to analyze the interfacial effect.Since the analysis requires the properties of the substrate,only five metals were evaluated.According to Eq.(10),the parameters εfw,ρw,Δ,σfw,and εffare necessary to determine αw.Based on our previous work,the corresponding parameters are listed in Table 5 [39].

Using the values of the properties listed in Table 5 in Eq.(10),αwwas calculated for each metal with the corresponding substrate.The calculation results are listed in Table 6 and illustrated in Fig.3.

As listed in Table 6,the corresponding αwfor Au,with an amorphous carbon film as the substrate,is relatively small.For the other four metals,with their corresponding substrates,the αwvalues are relatively large and the effect of the substrates becomes important.From the physical interpretation of αw,Au belongs to the category of non-ideal wettability,and the interfacial influence (interfacial effect) is not crucial.The other metals belong to the case of ideal wettability,and the interfacial effect is a crucial factor affecting the properties.This observation is consistent with the results in Section 3.4.2,and it also implies that the abnormal melting point depression is indeed caused by the interfacial effect.

(2) Relationship between δ and molecular parameter

To further study the relation between δ and the molecular parameters,/(d·εmm)=αw/(ρ·Δ·d) was calculated for the case with abnormal melting point depression.In calculation,d is the atomic diameter of the metals[40],and the corresponding valuesof σms,εms,and εmmwere taken from Table 5.The relationship between δ andis illustrated in Fig.4.A linear equation can be used to describe δ fromThis implies that the interfacial effect depends on the thermodynamic properties of both the metal and substrate;moreover,the interfacial effect can be estimated from the atomic diameter of the metals and the metal–metal and metal-substrate molecular interactions.

Table 5 Effective ε and σ for metal-substrates

Table 6 αw for different pairs of metal-substrate

Fig.3.The effect of αw for different pair of metal-substrate on the size-dependent melting point depression of metals.■:Pt;▼:Au;?:Pb;★:In;◆:Bi;?:Al.

Owing to the limited availability of experimental data with the information of substrates,the model prediction capacity cannot be further verified.New and systematic data obtained experimentally or by molecular simulation are needed to study the physical meaning of δ,which will be the focus of our future work.

3.5.2.Model performance

To further illustrate the performance of the newly modified model in this work,Fig.5 shows the deviation of the melting point calculated using the new model from that available in the literature data.For most metals,the deviation of the new model is within ±20 K throughout the size scale,which implies that the new model can be used to represent the abnormal melting point depression of metals.Compared to the results shown in Fig.1,the deviation in the results for this model was further decreased compared to the Gibbs–Thomson equation with surface-energy modification (Lee’s model) and the original Gibbs–Thomson equation.For the newly modified model,the deviation is (-50-20) K,while the deviations in the original and modified Gibbs–Thomson equations are (up to 700) K and (-350 to 100) K,respectively.

Fig.4.Linear fitting for further prediction according to δ and (·εms )/(d·εmm ).

Fig.5.Deviation of the melting point calculated with new model from the literature data and the comparison with the original Gibbs–Thomson equation and modified Gibbs–Thomson equation(surface-energy modification).■:Pt;▼:Au;▲:Mg;●:Ni;?:Pb;★:In;◆:Bi;?:Al.

4.Conclusions

In this work,the melting point depression of metals was studied using the original and modified Gibbs-Thomson equations.The results showed that the deviation of original Gibbs–Thomson equation occurred when the particle radius was less than 10 nm and the deviation was up to 700 K.The surface-energy modified Gibbs–Thomson equation can be used to predict the melting point depression for metals without interfacial effects.Further modification and investigation of the Gibbs–Thomson equations by considering both the surface-energy and interfacial effects revealed that(1) the modification with the interfacial effect was necessary for particles smaller than 5 nm in radius and (2) the interfacial effect can be predicted quantitatively from the thermodynamic properties of the metal and the substrate.The analysis of αwshows that it can be used to qualitatively interpret the interfacial effect,and a large αwvalue refers to a strong interfacial impact.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

Financial supports from Key Project(21838004),Joint Research Fund for Overseas Chinese,Hong Kong,Macao Young Scientists of National Natural Science Foundation (21729601) of China are acknowledged.X.Y.would like to thank the Swedish Research Council and the Kempe Foundation for financial support.

主站蜘蛛池模板: 午夜视频免费一区二区在线看| 亚洲欧洲日产国产无码AV| 99在线国产| 丁香六月激情综合| 日韩A∨精品日韩精品无码| 国产91色在线| 91无码国产视频| 国产精品永久在线| 91久久夜色精品国产网站| 一本视频精品中文字幕| 少妇被粗大的猛烈进出免费视频| 女人18一级毛片免费观看| 手机在线国产精品| 成年免费在线观看| 日本亚洲欧美在线| 五月天婷婷网亚洲综合在线| 亚洲一区二区成人| 国模私拍一区二区| 亚洲无码在线午夜电影| 国产一级视频在线观看网站| 99久久精品免费视频| 亚洲欧美h| 性色一区| 成人国产精品一级毛片天堂| 国产女人综合久久精品视| 日韩色图在线观看| 亚洲欧美不卡| 女人爽到高潮免费视频大全| 色成人亚洲| 国产欧美日韩91| 欧美日韩精品综合在线一区| 国产区免费精品视频| 亚洲第一视频免费在线| 国产丝袜啪啪| 人人91人人澡人人妻人人爽| 欧美日韩亚洲综合在线观看 | 91久久国产热精品免费| 国产电话自拍伊人| 亚洲综合色婷婷| 色国产视频| 国产剧情一区二区| 日韩欧美国产中文| 日韩欧美国产三级| 亚洲中文字幕av无码区| 爆乳熟妇一区二区三区| 青青操国产| 中国成人在线视频| 国产最新无码专区在线| 国产欧美又粗又猛又爽老| 免费99精品国产自在现线| 一级成人a做片免费| 青青青国产视频| 国内视频精品| 亚洲欧美在线综合一区二区三区 | 日韩色图在线观看| 精品自窥自偷在线看| 国产美女免费| 国产日本一线在线观看免费| 国产精品免费p区| 国产综合在线观看视频| 国产va免费精品| 欧美激情综合一区二区| 国产成人无码Av在线播放无广告| 嫩草国产在线| 人妻少妇乱子伦精品无码专区毛片| 黄网站欧美内射| 久久久91人妻无码精品蜜桃HD| 亚洲人成色77777在线观看| 亚洲一区二区日韩欧美gif| 国产亚洲欧美日韩在线一区二区三区| 国产99久久亚洲综合精品西瓜tv| 天天色天天操综合网| 手机在线国产精品| 欧美国产菊爆免费观看| 国产精品亚洲综合久久小说| 国产精品香蕉在线| 国产人碰人摸人爱免费视频| 国产精品无码久久久久久| 91美女视频在线| 伊人AV天堂| 五月丁香在线视频| 亚洲高清资源|