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.

主站蜘蛛池模板: 欧美另类第一页| 欧美一区二区三区香蕉视| 噜噜噜久久| 91精品伊人久久大香线蕉| 久久永久精品免费视频| 亚洲中文无码av永久伊人| 亚洲永久精品ww47国产| 夜夜爽免费视频| 国产美女久久久久不卡| 国产精品一区在线麻豆| 成人第一页| AⅤ色综合久久天堂AV色综合| 狠狠ⅴ日韩v欧美v天堂| 国产欧美精品一区二区| 91亚洲视频下载| 国产打屁股免费区网站| 青青青亚洲精品国产| 中文字幕 日韩 欧美| 久久久久人妻一区精品色奶水| 亚洲一区波多野结衣二区三区| 免费看av在线网站网址| 青青草原偷拍视频| 97se综合| 国产精品美人久久久久久AV| 免费中文字幕在在线不卡| 亚洲视频欧美不卡| 久久免费成人| 国产精品xxx| 91九色最新地址| 在线观看无码av免费不卡网站| 天天综合亚洲| 国模私拍一区二区| 宅男噜噜噜66国产在线观看| 欧美自慰一级看片免费| 色婷婷丁香| 国产农村妇女精品一二区| 久久精品aⅴ无码中文字幕 | 青青操国产| AV网站中文| 国产三级a| 久久久久久尹人网香蕉| 干中文字幕| jizz国产在线| 国产成人av大片在线播放| 国产精品内射视频| 色综合a怡红院怡红院首页| 国产成人欧美| 国产一级做美女做受视频| 性做久久久久久久免费看| 午夜啪啪网| 亚洲天堂网视频| 免费观看成人久久网免费观看| 亚洲久悠悠色悠在线播放| 伊人久久大线影院首页| 99久久国产自偷自偷免费一区| 亚洲无码91视频| 污网站免费在线观看| 日韩精品一区二区三区免费在线观看| 久草网视频在线| 无套av在线| 国模沟沟一区二区三区 | 亚洲激情99| 99re在线免费视频| 国产成人久久777777| 高清不卡一区二区三区香蕉| 欧美成人aⅴ| 亚洲日韩第九十九页| 免费A级毛片无码免费视频| 国产亚洲精品自在久久不卡 | 九色综合视频网| 国产呦精品一区二区三区网站| 国产不卡在线看| 伊人精品成人久久综合| 呦女亚洲一区精品| 青草娱乐极品免费视频| 久久精品无码专区免费| 国产农村1级毛片| 欧美人与动牲交a欧美精品| 久久精品人妻中文系列| 全免费a级毛片免费看不卡| 一级片一区| 超清无码熟妇人妻AV在线绿巨人 |