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

    Kinetic modeling of high-temperature oxidation of pure Mg

    2020-12-18 11:28:32FngzhouXingLijunZhng
    Journal of Magnesium and Alloys 2020年3期

    S M,Fngzhou Xing,N T,Lijun Zhng,?

    a State Key Laboratory of Powder Metallurgy,Central South University,Changsha,Hunan 410083,P.R.China

    b Max-Planck-Institut für Eisenforschung GmbH,Max-Planck-Stra?e 1,40237,Düsseldorf,Germany

    Received 14 June 2019;received in revised form 24 December 2019;accepted 26 December 2019 Available online 2 June 2020

    Abstract A variety of experimental tracer diffusivities of Mg and O in magnesium oxide available in the literature were firs assessed.Atomic mobilities including bulk and short-circuit diffusion of Mg and O were then obtained by means of the CALPHAD(Calculation of Phase Diagram)approach.Afterwards,the diffusion-controlled kinetic model of oxidation in a gas-MgO-Mg environment was developed based on the moving boundary model and Fick's law,coupling with the modifie thermodynamic description of MgO.A mathematical expression for parabolic rate constant kp of the oxide scale was derived for magnesia and correlated with the thermodynamic and diffusion kinetic information.The evaluated kp results were in line with the experimental data.Finally,the oxidation process of pure magnesium at 673 K was model-predicted,and the predicted evolution of the oxide thicknesses agreed very well with the experimental data.It was indicated that the grain boundaries diffusion of magnesium cations predominated the high temperature oxidation process.? 2020 Published by Elsevier B.V.on behalf of Chongqing University.This is an open access article under the CC BY-NC-ND license.(http://creativecommons.org/licenses/by-nc-nd/4.0/)Peer review under responsibility of Chongqing University

    Keywords:MgO;Diffusion;Oxidation;Atomic mobility;CALPHAD;Wagner's model.

    1.Introduction

    Magnesium alloys possess excellent properties including low density,high specifi strength and rigidity,good damping and cutting performance,great degradation behavior and biocompatibility,and thus serve as the key materials for realizing the lightweight of automobiles and biomedical applications[1-4].However,the corrosion behavior limits the practical applications of magnesium alloys[5-7].In order to improve the serviceability of magnesium alloys,it is crucial to study the oxidation of pure magnesium and even magnesium alloys,which will be helpful for comprehensive understanding of their oxidation mechanisms.

    So far,there have been several experimental measurements of the oxidation process of pure magnesium[8-12].However,it is very difficul to describe the oxidation behaviors over the wider temperature range through purely experimental measurements due to their relatively high time/money cost.Under this situation,some underlying modeling/numerical simulation techniques are highly needed to describe the oxidation process.In actual,the oxidation process of metals is a really complex one,involving the chemical reactions between metal and oxygen,as well as diffusion of cations and anions in oxides[13].Considering the fact that the rate of chemical reactions depends on the diffusion of ions,the overall oxidation rate is mainly controlled by the mass transport through oxides.The general Wagner's model[14]provides an understanding of mass transport processes occurring within a growing oxide scale,thereby leading to a capacity to predict the effect on the changes of the oxidation rate along with the temperature.In the Wagner's model[14],the complex process,including the evolution of microstructures due to diffusion of ions,grain growth,defect annihilation,spallation,etc.,was simplified and the parabolic equation was derived by considering only the diffusion-controlled process.It should be noted that the parabolic rate constant in Wagner's model should generally depend on both thermodynamic and kinetic parameters,including oxygen activity and bulk diffusion coefficient Very recently,Xing et al.[15]derived an equation for parabolic oxidation rate constant of Ni accounting for various diffusivities within the Wagner's framework.However,in Ref.[15],the empirical degree of defect ionization and vacancy concentration were introduced to replace the thermodynamic factor,which is not a general treatment for different metals/alloys.Thus,there is a need to correlate the parabolic rate constant with the general thermodynamic and kinetic information.

    Fig.1.(Color online)Calculated phase diagram of binary Mg-O system:(a)over the entire composition range;(b)in the vicinity of MgO compound,according to the thermodynamic descriptions originally from Hallstedt[16]together with the modifie thermodynamic descriptions for MgO in this work.

    Up to now,only one set of thermodynamic descriptions[16]for the Mg-O binary system exist in the literature,while no relevant kinetic descriptions have been reported.The Diffusion-Controlled Transformation(DICTRA)software package,operating under the Calculation of PHAse Diagram(CALPHAD)framework,has been widely employed to establish the atomic mobility databases for different materials,from which various diffusivities of composition and temperature dependence can be predicted together with the reliable thermodynamic databases.During the past several years,DICTRA has been employed to develop the atomic mobility descriptions of the oxides for different metals,including pure Fe[17,18],Cr[19]and Al[20],as well as some complex oxides LaCoO3-δ[21].

    Consequently,the oxidation process of pure Mg at high temperatures is chosen as the target in the present work.The major objectives are:(i)to perform critical review of all the experimental tracer diffusivities of Mg and O in MgO in the literature,and obtain the self-consistent atomic mobility descriptions of O and Mg in MgO;(ii)to derive a mathematical expression for growth rate constantkpof oxide scale by fully considering diffusion and thermodynamics in metal and oxygen ions in the framework of the Wagner's theory;and(iii)to predict the growth process of magnesium oxide during oxidation of pure magnesium,and clarify the oxidation mechanism of pure magnesium.

    2.Literature review on various diffusivities in MgO

    The diffusion behavior of magnesium cations in MgO has been extensively studied during the past years[22-29],and different types of diffusion coefficient of ions in Mg-O system are critically reviewed in the following and concisely summarized in Table 1.

    The self-diffusivities of magnesium cations in MgO single crystals over the temperature range of 1400~1600°C were firs measured in 1957 by Lindner and Parfit[22],who used the isotope exchange technology with the isotope28Mg.In their experiments,the impurities in magnesia were ignored,thus the magnesia was regarded as ideal crystals and the diffusion was intrinsic.They found that the activation energy of magnesium diffusion at low temperatures(<1400°C)was lower than that at higher temperature(>1400°C).After that,several further investigations on diffusion of magnesium in pure MgO single crystals were reported by Harding et al.[23,24]and Wuensch et al.[25].Harding et al.[24]compared their experiment data with those by Lindner and Parfit[22]and found that self-diffusion of magnesium in MgO is sensitive to the impurities and easily affected by short circuit diffusion at low temperatures.Wuensch et al.[25]employed26Mg rather than28Mg to measure the self-diffusion coeffi cients in an argon atmosphere.The experimental data measured by Harding et al.[23,24]are thus different from those by Lindner and Parfit[22].Moreover,the initial materials in the study of Wuensch et al.[25]contained even more impurities than those used by Lindner et al.[22],resulting in the inaccurate experimental self-diffusivities by Wuensch et al.[25].

    In the aspect of short circuit diffusion of magnesium ions in MgO,Sakaguchi et al.[27]utilized26Mg to measure the diffusion coefficient of magnesium in single-crystal MgO,taking both bulk and dislocation diffusion into consideration.Orman et al.[28]measured the grain boundary diffusion coefficient of magnesium in polycrystal MgO at 2273 K.However,their experiments were under the condition of high pressure from 15 to 25 GPa.Recently,the molecular dynamic simulations were carried out by Landuzzi et al.[29]and the grain boundary diffusion coefficient of Mg and O in MgO were also evaluated.

    Table 1Summary of all the experimental diffusivities of Mg and O in MgO available in the literature

    Up to now,there have been many studies on diffusion coefficient of oxygen in MgO available in the literature[28-40].In 1960,Oishi and Kingery[30]firs measured the selfdiffusion coefficient of oxygen anions in single-crystal MgO,and their measured self-diffusion coefficient of oxygen are smaller than those of the magnesium cations by two orders of magnitudes.However,their samples comprise a small amount of impurities.In 1971,the grain boundary diffusion coeffi cients of oxygen in MgO with bicrystals were measured by Mckenzie et al.[31]by means of Fisher's model[41].Their results indicated that grain boundary diffusion would dominate predominantly the mass transport property in polycrystalline MgO at low temperatures.After that,the self-diffusion coefficient of oxygen were determined by Hashimoto et al.[32]in polycrystalline MgO by isotope exchange technology and they concluded that the diffusion of oxygen along the grain boundaries are relatively faster than the bulk diffusion of oxygen.One year later,Shirasaki and Hama[33]reported that the self-diffusivities of oxygen in loosely-sintered and wellsintered polycrystalline MgO,following the same approach by Hashimoto et al.[32].Thereafter,Oishi et al.[34]measured two groups of the oxygen diffusivities in 1983,and the obvious differences between diffusion mechanisms at different temperatures were observed.They claimed that diffusion was impurity-insensitive at high temperatures while structuresensitive at low temperatures,i.e.,that the existence of dislocations and subgrains could increase the diffusion coeffi cients,which is actually in line with the diffusion mechanism of magnesium in MgO.In the same year,Henriksen et al.[35]determined the diffusivities of oxygen in pure,doped and deformed samples at 1673 K by means of the secondary-ion mass spectrometer(SIMS)method,and explained that the rate of oxygen diffusion in deformed materials was about fourfold larger than that in undeformed samples.Meanwhile,Reddy and Cooper[36]measured the diffusion coefficient of oxygen in MgO in the temperature range of 1580 to 1820 K.Then,18O was used as a tracer element in single crystals,deformed single crystals and polycrystals of MgO,and thus bulk,dislocation and grain boundary diffusion coefficient of18O were measured by means of SIMS method by Dolhert[37].In 1994,Yang and Flynn[38]measured the diffusivities of oxygen anions using high-purity magnesium oxides produced by molecular beam epitaxy(MBE),which is smaller than others by about two orders of magnitude.After that,Liberatore and Wuensch[39]measured the diffusion coefficient of oxygen in bulk and grain boundaries by using MgO bicrystals with high purity prepared by chemical vapor transport(CVT)and concluded that the enhancement of diffusivities was more conspicuous along grain boundaries with higher energy.In 2002,Yoo et al.[40]evaluated the prefactor and activation energy of oxygen diffusion coefficient in larger range of temperature and the smaller standard deviations were observed.In 2015,the molecular dynamic simulations were also applied by Landuzzi et al.[29]to evaluate the grain boundary diffusion coefficient of oxygen that deviate largely from the experimental data[32,33,37,39].

    3.Model descriptions

    3.1.Thermodynamic model

    The Mg-O phase diagram was reported by Wriedt[42]in 1987.It was pointed out by Wriedt[42]that MgO is not an ideal stoichiometric compound but with a very narrow homogeneity range.However,the specifi homogeneity range has never been reported so far.In 1993,Hallstedt[16]assessed the thermodynamic description of the Mg-O system by means of CALPHAD approach.Hallstedt[16]firstl took the solubility of oxygen in liquid magnesia into account,while the solubility of oxygen was not considered in the solid magnesia,and MgO was thus treated as an ideal stoichiometric compound.In addition,the site occupations of ions in lattice were regular,so the formula unit of solid magnesia with two sublattices was written as(Mg2+)1(O2?)1,where the firs sublattice denotes the octahedral interstices and the second one was occupied by oxygen anions in a face-centered cubic(fcc)lattice.It should be noted that the vacancies were not considered in the formula unit,and thus the MgO phase was treated as a rigid stoichiometric phase.

    Although the homogeneity range of MgO is too narrow to be measured accurately,it should not be neglected in the study of diffusion kinetics.The similar treatment in our previous studies on binary Co-Si[43]and Al-O systems[20]is also applied to the present Mg-O system,i.e.,that very tiny homogeneity range can be introduced in MgO by modifying its thermodynamic descriptions but keeping the deviation of the invariant reaction temperatures within 2 K.Based on the thermodynamic description by Hallstedt[16],the new sublattice model(Mg2+,Mg3+,Va)1(O2?)1was employed for MgO by introducing the magnesium cations with high valence to maintain the electro-neutrality.As a result,the molar Gibbs energy of MgO phase can be written as:whereis site fraction ofXin theisublattice,0GA:Bis the mole Gibbs energy of the end member of compound,andGexthe excess Gibbs energy.In this work,0GMg3+:O2?=0GMg2+:O2?+70000,0GVa:O2?=0 andGex=0 were set,while the other parameter0GMg2+:O2?was kept as the reported values in Hallstedt's original work[16].

    Based on thermodynamic descriptions originally from Hallstedt[16]together with those for MgO modifie in this work,the phase diagram of binary Mg-O system was calculated,as presented in Fig.1.As can be seen in the enlarged part of the Mg-O phase diagram,a very narrow homogeneity range exists in the MgO compound now.With the derived thermodynamic descriptions for MgO in Eq.(1),the site fractions of vacancy at the stoichiometric composition of MgO under the condition ofP(O2)=0.21 atm was computed and plotted in Fig.2(a).Moreover,the corresponding thermodynamic factor and chemical potential of O in MgO as a function of O compositions at 673 K were also calculated and shown in Fig.2(b).As can be seen in Fig.2,the site fraction of vacancy increases as the temperature increases,while the chemical potential profil indicates that the direction of driving force is consistent with the concentration gradient.Moreover,as the mole fraction of oxygen is close to 0.5,the thermodynamic factor will approach to the infinit.Under this circumvent,the so-called average thermodynamic factor as proposed in[43]is also introduced here.The detailed discussion can be referred to the Appendix.

    3.2.Diffusion kinetic model

    The oxidation rate of magnesium at high temperatures is controlled by diffusion[44],which is driven by the gradient of chemical potential.Based on the °Agren's theory,the flu of the ions in the lattice-fi ed reference frame is given as[17]

    whereμiis the chemical potential,V0the volume per mole of anion sites.MiVais atomic mobility of componentijumping from one site to a neighboring site and has the form:MiVa=M0exp(-Q/RT),in whichM0is the frequency factor,which is the function of vibrational frequency(ν),jump distance(δ)and temperature(T)viaM0=δ2ν/RT,whileQis the activation energy.Here,by assuming that atomic mobilities of magnesium cations in the frst sublattice have the relationship MMg2+Va=MMg3+Va=MMgVa,Eq.(2)can be simplifie as

    According to Ref.[17],one can obtain the expression for the tracer diffusion coefficient D?Mg,which is given by

    Compared with Einstein's equation:D?=RTM,the atomic mobility of magnesium in MgO can be inferred

    Fig.2.(Color online)Calculated(a)site fractions of vacancy at the stoichiometric composition of MgO under the condition of P(O2)=0.21 atm;(b)thermodynamic factors(red solid line)and chemical potentials(blue dash line)of O in MgO as function of O composition at 673 K,based on the presently obtained thermodynamic descriptions for MgO with slight modification of the original parameters from Hallstedt[16].

    As mentioned above,the ions in MgO diffuse according to the vacancy mechanism.Although the second sublattice is fully occupied by oxygen anions in the thermodynamic model above,there are a few sites for vacancies,of which the contribution is too small to be considered in the thermodynamic model.Supposing that the site fractions of species in the second sublattice are yIIVaand yIIO,the flu of the tracer O?alongzaxis in a lattice-fi ed frame of reference can also be expressed in the similar way,

    which yields:

    Then,the atomic mobility of oxygen has the form:

    It should be noted that there are several types of diffusivities besides the bulk diffusivity in normal lattice.As a matter of fact,there are various defects in real crystals,like dislocations and grain boundaries in polycrystals,which can largely enhance the diffusion velocity through reducing the activation energy.Considering the influenc caused by defects,it is necessary to use effective diffusion coefficient to replace the bulk diffusivity in the present work.To quantitatively analyze the diffusion behavior in real polycrystals,the model proposed by Hart[45]and Harrison[46]was utilized to describe the effective diffusion coefficient in this work[47].The same frequency factor with bulk diffusion but a redistributed coefficienFredfor activation energy was introduced in the circuit diffusion compared with the bulk atomic mobility.Then,the atomic mobilities along dislocation and grain boundaries can be expressed as

    where the value of redistributed coefficien in the equation is generally less than 1 for the circuit diffusion.Fredis in relation to the energy of the defects,including Burger's vector and thickness of grain boundaries,where Burgers vector's square represents the energy of a dislocation,and the grain boundary thickness,to some extent,represents the energy of a grain boundary.Besides,the effective diffusion coefficient are connected to the quantity of the defects relating to the key factors,namely,the dislocation density and grain size.Thus,the effective mobility and diffusion coefficient can be expressed as[45,46]

    hereδdenotes the average thickness of grain boundaries,daverage grain size which depends on the temperature and time,ρdislocation density andbBurgers vector.Dbulk,DDIS,DGBandDeffare the tracer diffusion coefficient in bulk,dislocation and grain boundary and the effective tracer diffusion coefficients respectively.

    3.3.Oxidation kinetic model

    The diffusion-controlled oxidation process obeys the parabolic law.According to the Wagner's theory[14]:

    whereXrepresents the oxide thickness(unit:m),kpoxide growth rate constant(unit:m2/s)andttime(unit:s).Based on Eq.(11),the key to evaluate the thickness of oxides as a function of time lies in the reasonable value ofkp,which is determined by both thermodynamics(driving force)and kinetics(diffusion coefficients)In fact,the reliable thermodynamic and kinetic information needed forkpcan be provided by the CALPHAD-type thermodynamic and kinetic descriptions.The derivation of the mathematical expression forkpis described in the following.

    The oxidation rate can be determined by the moving rate of interfaces,including both oxide/metal and oxide/gas interfaces.The inward diffusion of oxygen ions results in an internal migrating interface towards pure Mg,while the outward diffusion of magnesium ions can cause an external migrating interface towards gas.The schematic diagram is shown in Fig.3.Based on the moving boundary model[20],the moving rate of boundary has the form:

    wherevdenotes moving rate of interface,the molar volume ofi(i=α,β)phase and Jik the diffusion flu of componentkiniphase,which can be expressed by Fick's frst law:

    Here,Dkdenotes chemical diffusion coefficien of componentk.Dkcan be related to the tracer diffusion coefficienand thermodynamic factor[16]via the following equation:

    Since thermodynamic factor for a nearly stoichiometric compound is infinit,the average thermodynamic factor define as[43]is introduced and the detail discussion is listed in the Appendix.The combination of Eqs.(12)-(14),and may determine the moving rates of two interfaces of oxide phase,respectively.Therefore,the values ofkpcan be determined by the sum of moving rates of interfaces in both sides and finall simplifie into:

    whereΔμO(píng)denotes the difference of chemical potential of oxygen in magnesium oxide.According to the local equilibrium hypothesis,the identical components in the two phases have the same chemical potentials at interface.Thus,the chemical potentials of oxygen at two interfaces equal to those of gas and Mg(hcp)phase respectively.Moreover,since magnesium oxide owns a very narrow homogeneity range,i.e.,that MgO is nearly a stoichiometric compound,the variation of chemical potential of oxygen can be expressed by:

    4.Results and discussion

    4.1.Bulk diffusion in MgO

    The atomic mobilities for bulk diffusion in MgO were assessed on the basis of various experimental diffusivities in the bulk MgO available in the literature by means of a non-linear least-squares optimization scheme incorporated in the PARROT module of the DICTRA software package.The experimental tracer diffusivities of Mg in MgO in Refs.[22-27]as displayed in Fig.4(a)were used to assess the atomic mobility parameters of Mg.On the contrary,as shown in Fig.4(b),the experimental data of O[30,32-40]have a larger scatter because the purity of samples reported in the literature are not identical,which is one of the relating key factors of diffusion coefficients The impurity contents may change the quantity of vacancies,ions interaction and substructures,thereby,influencin the tracer diffusion coefficients The samples studied in the Refs.[30,32,33]with the highest impurity correspond to the highest diffusion coefficient in Fig.4(b).Conversely,the samples in Ref.[38]with the highest purity possess the lowest diffusion coefficients Moreover,the samples with similar purity have the diffusion coefficient close with each other.Therefore,the experimental data from Ref.[38]were accepted to evaluate the atomic mobility descriptions for O in MgO.The other data in Fig.4(b)are considered as the results of the combined contributions of the bulk and short circuit diffusion,which will be discussed inSection4.3.The evaluated atomic mobility parameters are expressed as Eq.(17):

    Fig.4.(Color online)Calculated bulk diffusion coeff cients of(a)Mg and(b)O in MgO as a function of inverse temperature and oxygen partial pressure,compared with the experimental data[22-27,30,32-40].

    The calculated bulk diffusion coefficient of the magnesium cations and oxygen anions as a function of inverse absolute temperature under different oxygen partial pressures and temperatures are shown in Fig.4,compared with the experimental data[22-27].It should be noted that the partial oxygen pressure determines the vacancy concentration in MgO,which affects the diffusion coefficient according to Eq.(4)and(7).As can be seen in Fig.4(a),the calculated temperaturedependent bulk diffusion coefficient of magnesium increase slightly as the partial oxygen pressure increases and agree very well with most of the experimental data[22,24-27].It should be noted that the experimental data in Ref.[23]deviate from the calculated results because the impurities exist in the tracer solution,as explained in Ref.[24].As can be seem in Fig.4(b),the calculated bulk diffusion coeffi cients of oxygen also accord with the experimental data in Ref.[38].

    4.2.Short-circuit diffusion in MgO

    In practice,most materials are polycrystal,and thus a large amount of defects such as grain boundaries and dislocations coexist.As mentioned in Section 3.2,the jumping rates of atoms along the short-circuit paths such as grain boundaries and dislocations are much higher than that in the lattice.Therefore,the contributions of both grain boundaries and dislocations to diffusion should be taken into account.On the basis of Eq.(9),the mobilities along the short-circuit paths were evaluated by the same method declared above.

    The experimental diffusivities of Mg along the dislocations reported in Ref.[27],shown in Fig.5(a),were employed to assess the atomic mobility parameters.Moreover,the experimental data in Ref.[37],as shown in Fig 5(b),were used to assess the atomic mobility parameters of O along the dislocations,the obtained atomic mobility parameters along the dislocations are presented as follows:

    As shown in Fig.5,both calculated diffusivities of Mg and O conform to the experimental data in Ref.[27]and Ref.[37]respectively.Comparing the mobility along the dislocations with that in the bulk,the redistributed coefficient of activation energy along the dislocations can be calculated

    Besides the dislocation diffusion coefficients the grain boundary diffusion coefficient of Mg and O are also studied,and the results are displayed in Fig.6.For the oxygen anions,there are several groups of experimental data on the grain boundary diffusion coefficient available in the literature,including the data from the experimental measurements[31-33,37,39]and the molecular dynamics simulations[29].The molecular simulation results scatter much larger over several magnitudes than the experimental data,and thus were not used in this work.Moreover,the scattering between different sets of experiments[31-33,37,39]for the grain boundary diffusion coefficient of O is about two magnitudes.It is quite normal because the high-angle grain boundaries have higher energy than the low-angle grain boundaries along which the diffusivities are lower[39].In order to balance the influence from different categories of grain boundary,the average grain boundary diffusion coefficient of O were used by fittin to all the experimental data[31-33,37,39],and the results can be seen in Fig.6(b).

    Due to lack of data on magnesium diffusion along grain boundaries(except for one piece of information frommolecular dynamics simulations[29]which cannot be accepted in the present work),no assessment for magnesium mobility for the grain boundary diffusion was performed here.As discussed above,the redistributed coefficient of Mg and O along dislocations are nearly the same.Thus,it seems reasonable to assume that the grain boundary diffusion coefficient of Mg and O also have the same redistributed coefficient as adopted in this work to predict the grain boundary diffusion coefficient of Mg in MgO.The predicted grain boundary diffusion coefficient of Mg in MgO as a function of inverse temperature are presented in Fig.6(a).

    The finall obtained atomic mobility parameters of Mg and O diffusion along grain boundaries are given as

    Here,Eq.(19)represents the average grain boundary mobilities.As shown in Fig.6(a),the predicted grain boundary diffusion coefficient agree with the molecular dynamics results from Landuzzi et al.[29]at temperatures higher than 2000 K.When the temperature is the lower than 2000 K,the noticeable differences between the two values exist,and become larger and larger as the temperature decreases.As can be seen in Fig.6(b),the agreement between calculated results and the experimental data seems to be good.Moreover,the mobility for diffusion of O along the high-angle grain boundaries is higher by about ten times than the average value,which is however higher than the mobility for diffusion of O along the low-angle grain boundaries by around ten times.The redistributed coefficient of grain boundary diffusion coeffi cients of Mg and O can be evaluated as0.619.

    4.3.Effective diffusion in MgO

    With Eq.(10),the effective diffusion coefficient can be evaluated based on the obtained atomic mobility parameters for diffusion along the bulk and short-circuit paths.

    The factors influencin the diffusion along dislocation include dislocation density and burgers vector.As reported in the literature[27,37,48],the dislocation density is insensitive to the change of temperature.The values of dislocation density in MgO are in the range of 1×1013/m2~1×1015/m2[27,37].In this work,the typical valueρ=1×1014/m2was adopted.Besides,the burgers vector of main dislocationsconsidering that MgO is with the face-centered cubic structure[49-51].Based on the lattice parameters available in the literature,the length of burgers vector can be calculated as 0.297 nm[50].Therefore,the value ofρb2in Eq.(10)is calculated in the order of 10?5.One can obviously see that the contribution of diffusion along the dislocation to the effective diffusion is small and thus may be neglected.

    Fig.7.(Color online)Calculated average grain size of MgO as a function of temperature,compared with the experimental data[29,31-33,37,48,52,53].

    As mentioned inSection 4.1,all the experimental data except the data from Ref.[38]were not used to evaluate the atomic mobility parameters because the purity of samples was not satisfactory.That is because the impurities in those samples may largely enhance the diffusion coefficients To some extent,the content of impurities can be viewed as the dislocation since both of them can serve as the fast path for diffusion.As presented in Fig.4(b),the influence of different dislocation densities on the bulk diffusion are compared with the experimental data with different impurities,and the results show a good agreement with the experimental data.As for the grain boundary diffusion,the influenc factors include the thickness of grain boundary and the average grain size.As reported in the literature,the thickness can be viewed as a constant for simplification and thus the typical valueδ=1 nm[28,29,31,39]was adopted here.Moreover,as the temperature increases,the grains may coarse,resulting in the decrease of interface energy.Therefore,the average grain size should be a function of temperature,i.e.,d=Aexp(BT),whereddenotes the average grain size,whileAandBare the undefine constants.As shown in Fig.7,the temperaturedependent average grain sizes can be obtained by fittin to the experimental data[29,31-33,37,48,52,53]using a least-square method:

    Fig.8.(Color online)Model-predicted effective diffusion coeff cients of(a)Mg and(b)O as a function of inverse temperature,compared with the bulk and short-circuit diffusion coeff cients.

    wheredis the average grain size in the dimension of m,Ttemperature in K.

    Based on the obtained atomic mobility parameters and the values for factorsρ,b,d,δ,the effective diffusion coeffi cients can be evaluated,and displayed in Fig.8.As can be seen in Fig.8,at higher temperatures,the effective diffusion coefficient are close to the bulk diffusion coefficient but show large scattering at lower temperature conversely.It can be explained by the different grain sizes existing at different temperatures.As the temperature decreases,the grains become smaller and there are more grain boundaries to enhance the diffusion of ions.By contrast,at high temperatures,the large grains eliminate the effect of grain boundaries on diffusion.Therefore,the contributions of short-circuit paths to diffusion can be ignored at higher temperatures but turn to be more and more essential as the temperature decreases.For instance,as shown in Fig.8(b),the extrapolated value ofDObulkat 500 K is close to~10?60m2/s,while those ofDOGBandDODISare in the range of 10?40~10?35m2/s.It should be noted that the reliability of the obtained effective diffusion coefficient at lower temperatures needs further confirmatio because it is very difficul to measure the grains with several nanometers accurately.Moreover,by comparing Fig.8(a)with Fig.8(b),the effective diffusion coefficient of O are much lower than those of magnesium cations.Base on Eq.(15),it can be concluded that the high-temperature oxidation process of pure Mg is predominantly controlled by diffusion of the magnesium cations.

    4.4.Oxide growth rate constant kp and prediction of oxidation process of pure Mg

    Fig.9.(Color online)Logarithmic values of growth rate constant of oxide kp as a function of inverse temperature from 500 K to 2000 K,compared with the experimental data[8-12].

    Based on the kinetic model of oxidation inSection 3.3,the growth rate constant of oxideskpcan be evaluated based on Eqs.(15)and(16)by the thermodynamic and kinetic descriptions obtained in this work.In order to truly reveal the oxidation mechanism,different sets ofkpvalues were calculated by considering different types of diffusion coefficient(i.e.,bulk,grain boundary and/or effective diffusion coeffi cients).The calculated results and the experimental growth rate constants are shown in Fig.9,where the experimentalkps were calculated based on Eq.(11)together with the experimental thickness in Refs.[8-12],while the model-predicted values ofkpwere calculated by the bulk diffusion coefficient(black line),grain boundary diffusion coefficient(blue line)and effective diffusion coefficient with different grain boundary angles(red,green,orange lines),respectively.Based on Eq.(11),the oxide thickness as a function of time can be directly predicted by the evaluatedkp.The model-predicted thickness of MgO atP(O2)=0.21 atm and 673K is presented in Fig.10.It can be found that the model-predicted thickness according to the effective diffusivities with the high-angle grain boundaries can reproduce all the experimental data perfectly,which indicates the high-angle grain boundaries exist predominantly in the MgO polycrystals at 673 K.

    Fig.10.(Color online)Model-predicted evolution of thickness of MgO oxide at 673K and P(O2)=0.21 atm using different sets of diffusivities,compared with the experimental data[11,12].

    From Fig.9 and Fig.10,one can conclude(i)that the calculated results based on the effective diffusivities reproduce the experimental data satisfactorily,which illustrates the reliability of the presently evaluated atomic mobility descriptions for Mg and O in MgO;(ii)that the grain boundaries can largely enhance the oxidation rate,especially at low temperatures;and(iii)that the grain boundaries with high angles can enhance the diffusion coefficient more effectively than those with low angles.

    Though the predictedkpvalues and thicknesses of MgO during oxidation are consistent with the experimental data,one should bear in mind that the naturally formed oxide fil on Mg is quite diverse and depends on various factors.For example,the contributions of defects like grain boundary and dislocation were included in Eq.(10),but some empirical or typical parameters of defects were used due to lack of suffcient experimental data on i.e.,the thickness of grain boundary,and the dislocation density.Therefore,in order to perform the quantitative modeling of the real oxidation process in different metals/alloys,the accurate description of microstructure evolution during the oxidation process should be a prior,which will be the topic for next paper.

    5.Conclusions

    ·Atomic mobility descriptions for Mg and O in MgO were firs obtained by means of CALPHAD approach based on various diffusion coefficient critically reviewed from the literature,together with the reported thermodynamic descriptions.Both bulk diffusion and short-circuit diffusion were taken in account.Different diffusivities evaluated according to the obtained atomic mobilities agree well with most of the experimental data.

    ·A mathematical expression for growth rate constantkpof oxide scale by fully considering diffusion and thermodynamics in metal and oxygen ions was proposed in the framework of the Wagner's theory.The growth rate constantkpof magnesium oxide was predicted by the obtained atomic mobility descriptions,and the predictedkps show a good agreement with the experimental data.

    ·The growth of magnesium oxide as a function of time at 673 K was simulated based on the Wagner's model.The simulated results with high-angle grain boundaries reproduced the experiments satisfactorily.Furthermore,the mechanism of magnesium oxidation over a wide range of temperature was discussed.The results indicate that grain boundary diffusion of magnesium cations predominated during the high temperature oxidation process.

    Acknowledgements

    The financia support from the National Key Research and Development Program of China(Grant No.2016YFB0301101),the Hunan Provincial Science and Technology Program of China(Grant No.2017RS3002)-Huxiang Youth Talent Plan,and the Youth Talent Project of Innovation-driven Plan at Central South University(Grant No.2282019SYLB026)is acknowledged.Sa Ma acknowledges the financia support by the Fundamental Research Funds for the Central Universities of Central South University(Grant No.2019zzts485).The authors thank Prof.Bengt Hallstedt at Aachen University,Germany for kindly providing thermodynamic database of Mg-O system.

    Appendix.Derivation of oxide growth rate constant kp

    The modifie formula unit of solid magnesium oxide phase is given with two sublattices as:(Mg,Va)1(O)1,where the valences of ions are neglected.Nevertheless,to ensure the diffusion of oxygen in the second sublattice based on vacancy diffusion mechanism,the formula unit can be rewritten as(Mg,Va)1(O,Va)1.As reported in Ref.[13],the chemical diffusion coefficient of magnesium and oxygen ions in magnesium oxide can be expressed by:

    Here,ni(i=Mg,O)represents the number of speciesiper mole firs or second sublattice sites,which can be regarded as one during calculation since the fraction of vacancy is close to zero.Because the mole fractionxihas a relationship withni(xi=ni/(nMg+nO)),the chemical diffusion coefficien of O in Eq.(A1)can thus be deduced as:

    whereφis thermodynamic factor which can be determined from thermodynamic description.And the chemical diffusivity of Mg has the same form.Using the Gibbs-Duhem relationship,andxMg=xO=0.5,one can express the chemical diffusivity of Mg as:

    In the stoichiometric compound,the average thermodynamic factorφave[41]needs to be introduced to replaceφf(shuō)or evaluation of the diffusion coefficients

    wherex1andx2respectively denote minimum and maximum of mole fraction of oxygen in MgO,μ1andμ2are chemical potentials.According to Eq.(16),the chemical potential differenceΔμO(píng)could be calculated.

    As is mentioned in Section 3.3,the distance between the two interfaces of MgO phase related to the interface moving rate,is the thickness of oxide film Based on Eq.(12),interface moving rates could be deduced by chemical diffusion flux It should be noted that the mole fraction of O in magnesium and the mole fraction of Mg in gas are approximated to be zero,resulting in zero flu of component O and Mg in pure magnesium phase and gas phase,respectively.Moreover,the mole fractions of magnesium and oxygen in MgO are both 0.5.Therefore,combined with the Fick's frst law,the absolute values of interface moving rates could be derived as:

    where v1and v2respectively denote moving rate of metal/oxide interface and oxide/gas interface,Vmthe volume occupied by per mole ions(sum of magnesium and oxygen ions).Vmis considered as independent of composition.It is reasonable to assume that the concentration profile within MgO phase region is linear at any time as magnesium oxide has very narrow homogeneity range.Thus,the growth rate of magnesium oxide could be determined by interface moving rates:

    whereΔxOandΔxMgrespectively represent mole fraction difference of O and Mg,andΔxO=ΔxMg.Combine Eq.(11)with the solution of Eq.(A6),growth rate constant of oxide kpis derived as:

    国产精品一二三区在线看| 国产又色又爽无遮挡免| 成人国产av品久久久| 国产有黄有色有爽视频| 高清黄色对白视频在线免费看| 插逼视频在线观看| av福利片在线| av卡一久久| 美女视频免费永久观看网站| 女人久久www免费人成看片| 日韩av在线免费看完整版不卡| 天天操日日干夜夜撸| 永久免费av网站大全| 欧美精品人与动牲交sv欧美| 亚洲精品自拍成人| 亚洲国产欧美在线一区| 欧美日韩av久久| 国产日韩欧美在线精品| 久久热精品热| 中国国产av一级| 午夜激情av网站| 免费黄频网站在线观看国产| 夫妻性生交免费视频一级片| 精品国产露脸久久av麻豆| 视频中文字幕在线观看| √禁漫天堂资源中文www| 一区二区三区四区激情视频| 搡女人真爽免费视频火全软件| 菩萨蛮人人尽说江南好唐韦庄| av免费在线看不卡| av一本久久久久| 美女福利国产在线| 简卡轻食公司| 少妇的逼水好多| 伊人亚洲综合成人网| 自拍欧美九色日韩亚洲蝌蚪91| 插阴视频在线观看视频| 欧美精品高潮呻吟av久久| 美女脱内裤让男人舔精品视频| 草草在线视频免费看| 中文字幕av电影在线播放| 飞空精品影院首页| 成人免费观看视频高清| 精品久久久久久电影网| 一级毛片我不卡| 亚洲少妇的诱惑av| 精品熟女少妇av免费看| 久久久久久伊人网av| 一本—道久久a久久精品蜜桃钙片| 一级毛片电影观看| 亚洲欧美清纯卡通| 精品一区二区三卡| 欧美日韩成人在线一区二区| 91久久精品国产一区二区成人| 美女视频免费永久观看网站| 久久久久久久久久人人人人人人| 久久人人爽人人爽人人片va| 午夜激情福利司机影院| 精品少妇黑人巨大在线播放| 日韩制服骚丝袜av| 日韩亚洲欧美综合| 天天影视国产精品| freevideosex欧美| 亚洲国产精品一区三区| 蜜桃国产av成人99| 成人亚洲欧美一区二区av| 日本vs欧美在线观看视频| 日韩欧美精品免费久久| h视频一区二区三区| 欧美精品人与动牲交sv欧美| 在线观看美女被高潮喷水网站| 999精品在线视频| 国产在线视频一区二区| 亚洲av男天堂| 五月天丁香电影| 大又大粗又爽又黄少妇毛片口| 两个人的视频大全免费| 人人妻人人添人人爽欧美一区卜| 高清毛片免费看| 久久女婷五月综合色啪小说| 久久久国产一区二区| 亚洲欧美中文字幕日韩二区| 少妇的逼水好多| 免费人妻精品一区二区三区视频| 久久97久久精品| 免费黄网站久久成人精品| 波野结衣二区三区在线| 日韩av免费高清视频| 最近的中文字幕免费完整| 国产精品人妻久久久久久| 成人手机av| 久久久久久久精品精品| 涩涩av久久男人的天堂| 日韩av在线免费看完整版不卡| 久久午夜综合久久蜜桃| 国产爽快片一区二区三区| 久久精品人人爽人人爽视色| 国产亚洲欧美精品永久| 少妇熟女欧美另类| 黄色配什么色好看| 桃花免费在线播放| 国产黄色视频一区二区在线观看| 亚洲精品美女久久av网站| 日韩伦理黄色片| 好男人视频免费观看在线| 久久久国产欧美日韩av| 国产精品 国内视频| 亚洲av电影在线观看一区二区三区| 999精品在线视频| 国产白丝娇喘喷水9色精品| 色婷婷av一区二区三区视频| 最近2019中文字幕mv第一页| a级片在线免费高清观看视频| 亚洲天堂av无毛| 男男h啪啪无遮挡| 极品少妇高潮喷水抽搐| 亚洲精品亚洲一区二区| 丝袜美足系列| 日韩欧美一区视频在线观看| 观看av在线不卡| 亚洲成人手机| 黑人高潮一二区| 一本大道久久a久久精品| 国产色婷婷99| 哪个播放器可以免费观看大片| 美女福利国产在线| 精品视频人人做人人爽| 男人操女人黄网站| 中文字幕精品免费在线观看视频 | 尾随美女入室| 女人精品久久久久毛片| 久久精品国产亚洲网站| 欧美激情极品国产一区二区三区 | 一级毛片我不卡| 久久国产精品男人的天堂亚洲 | 丰满乱子伦码专区| 啦啦啦啦在线视频资源| 国产有黄有色有爽视频| 又粗又硬又长又爽又黄的视频| 好男人视频免费观看在线| 18禁动态无遮挡网站| 中文字幕人妻熟人妻熟丝袜美| 在线观看免费视频网站a站| 国产精品久久久久久久电影| 久久久久精品久久久久真实原创| 久久国产亚洲av麻豆专区| 国产精品久久久久久久久免| 少妇被粗大的猛进出69影院 | 国产欧美日韩一区二区三区在线 | 天美传媒精品一区二区| 99久久精品国产国产毛片| 久久精品国产亚洲网站| 天天影视国产精品| 另类亚洲欧美激情| 亚洲国产日韩一区二区| 中文天堂在线官网| 亚洲av.av天堂| 人人妻人人添人人爽欧美一区卜| 丝袜脚勾引网站| 亚洲少妇的诱惑av| 黄色配什么色好看| 中文欧美无线码| 久久久国产一区二区| 欧美日韩国产mv在线观看视频| 亚洲成人一二三区av| 日韩一区二区视频免费看| 天堂中文最新版在线下载| 九九久久精品国产亚洲av麻豆| 青青草视频在线视频观看| 国产精品99久久99久久久不卡 | 伦理电影大哥的女人| 久久99精品国语久久久| 国产精品免费大片| 午夜福利影视在线免费观看| 内地一区二区视频在线| 亚洲,欧美,日韩| 久久久亚洲精品成人影院| 我的女老师完整版在线观看| 亚洲色图综合在线观看| 婷婷色av中文字幕| 少妇被粗大的猛进出69影院 | 黄色一级大片看看| 亚洲av国产av综合av卡| 又粗又硬又长又爽又黄的视频| 国产在视频线精品| 亚洲国产欧美日韩在线播放| 免费看不卡的av| 国产成人精品无人区| 最近最新中文字幕免费大全7| 美女福利国产在线| 国产无遮挡羞羞视频在线观看| 日韩视频在线欧美| 全区人妻精品视频| 亚洲美女视频黄频| 校园人妻丝袜中文字幕| 男人爽女人下面视频在线观看| 国产乱人偷精品视频| 久久女婷五月综合色啪小说| 欧美日韩视频精品一区| 久久综合国产亚洲精品| 97精品久久久久久久久久精品| 免费黄频网站在线观看国产| 久久国产精品大桥未久av| av播播在线观看一区| 久久久国产一区二区| 亚洲色图综合在线观看| a 毛片基地| 久久精品国产a三级三级三级| 欧美变态另类bdsm刘玥| av.在线天堂| 极品人妻少妇av视频| 丰满迷人的少妇在线观看| 春色校园在线视频观看| 啦啦啦啦在线视频资源| 亚洲av成人精品一区久久| 成人免费观看视频高清| 日产精品乱码卡一卡2卡三| 又大又黄又爽视频免费| 九色成人免费人妻av| av国产精品久久久久影院| 午夜福利影视在线免费观看| 国产精品国产三级国产专区5o| 欧美 日韩 精品 国产| 国产亚洲av片在线观看秒播厂| 久久国产精品男人的天堂亚洲 | 男人添女人高潮全过程视频| 亚洲在久久综合| 亚洲第一区二区三区不卡| 欧美日韩在线观看h| 成人国语在线视频| 在线观看国产h片| 欧美日韩国产mv在线观看视频| 99热这里只有是精品在线观看| 欧美日韩亚洲高清精品| 国产精品不卡视频一区二区| 在线看a的网站| 国产熟女午夜一区二区三区 | 国产精品国产三级国产专区5o| 亚洲精品乱码久久久v下载方式| 久久久精品免费免费高清| 亚洲国产日韩一区二区| 国产免费福利视频在线观看| 免费人妻精品一区二区三区视频| 国产一区二区在线观看av| 制服人妻中文乱码| 最后的刺客免费高清国语| 亚洲五月色婷婷综合| 18禁观看日本| 国产伦理片在线播放av一区| 国产免费一级a男人的天堂| 少妇丰满av| 国产精品一区二区在线不卡| 丝袜美足系列| 亚洲精品国产av蜜桃| 看十八女毛片水多多多| 中国三级夫妇交换| 欧美日韩国产mv在线观看视频| 亚洲图色成人| 久久久久久久久久久免费av| 亚洲精品自拍成人| 国产欧美另类精品又又久久亚洲欧美| 最近中文字幕高清免费大全6| 中文字幕免费在线视频6| 欧美激情国产日韩精品一区| 欧美日韩视频精品一区| 九色成人免费人妻av| 简卡轻食公司| 亚洲国产毛片av蜜桃av| 老司机影院毛片| 久久久久久久久大av| 春色校园在线视频观看| 国产爽快片一区二区三区| 亚洲精华国产精华液的使用体验| freevideosex欧美| 校园人妻丝袜中文字幕| 日韩一区二区视频免费看| videosex国产| 蜜臀久久99精品久久宅男| a级毛片免费高清观看在线播放| 久久精品久久久久久久性| 免费大片黄手机在线观看| 菩萨蛮人人尽说江南好唐韦庄| 在线观看免费高清a一片| 亚洲图色成人| 国产免费现黄频在线看| 免费人妻精品一区二区三区视频| 久久97久久精品| 午夜激情福利司机影院| 熟女av电影| 亚洲熟女精品中文字幕| 99久久综合免费| 老女人水多毛片| 久久人妻熟女aⅴ| 日韩大片免费观看网站| 欧美日韩视频高清一区二区三区二| a 毛片基地| 少妇被粗大的猛进出69影院 | 日韩欧美精品免费久久| 丝袜在线中文字幕| 香蕉精品网在线| 久久av网站| av免费观看日本| 卡戴珊不雅视频在线播放| 内地一区二区视频在线| 国产成人午夜福利电影在线观看| 99re6热这里在线精品视频| 夜夜看夜夜爽夜夜摸| 永久网站在线| 亚洲熟女精品中文字幕| 成年人午夜在线观看视频| 国产精品一区二区在线观看99| 久久综合国产亚洲精品| 日韩成人伦理影院| 久久午夜综合久久蜜桃| 午夜免费男女啪啪视频观看| 男女边吃奶边做爰视频| 亚洲精品久久成人aⅴ小说 | 亚洲人成网站在线播| 大香蕉久久成人网| a级毛色黄片| 国产永久视频网站| 美女国产高潮福利片在线看| 亚洲欧美一区二区三区黑人 | 校园人妻丝袜中文字幕| 一本一本久久a久久精品综合妖精 国产伦在线观看视频一区 | 国产高清不卡午夜福利| 男男h啪啪无遮挡| 久久精品国产亚洲av涩爱| 黑人猛操日本美女一级片| 国产日韩欧美在线精品| xxx大片免费视频| 极品少妇高潮喷水抽搐| 久久久久久久久大av| 99热这里只有精品一区| 99热国产这里只有精品6| 啦啦啦啦在线视频资源| 国产一区二区三区av在线| 亚洲欧美一区二区三区国产| 一级,二级,三级黄色视频| 国产精品熟女久久久久浪| 99视频精品全部免费 在线| 黄色毛片三级朝国网站| 国产在线免费精品| 卡戴珊不雅视频在线播放| 日本黄色日本黄色录像| 国产av国产精品国产| 人人妻人人爽人人添夜夜欢视频| 一级二级三级毛片免费看| h视频一区二区三区| xxx大片免费视频| 免费黄频网站在线观看国产| 一本—道久久a久久精品蜜桃钙片| 亚洲av在线观看美女高潮| 少妇丰满av| 欧美日韩亚洲高清精品| 国产极品粉嫩免费观看在线 | 亚洲精品亚洲一区二区| 亚洲精品成人av观看孕妇| 人成视频在线观看免费观看| 精品人妻熟女av久视频| 日韩人妻高清精品专区| 99热6这里只有精品| 亚洲精品aⅴ在线观看| 大香蕉97超碰在线| 永久网站在线| 亚洲欧美中文字幕日韩二区| 久久免费观看电影| 一级黄片播放器| 久久热精品热| 一边摸一边做爽爽视频免费| 91国产中文字幕| a级毛色黄片| 超色免费av| 人妻系列 视频| 久久免费观看电影| 国产精品国产三级国产专区5o| 亚洲精品自拍成人| 建设人人有责人人尽责人人享有的| 亚洲av中文av极速乱| 亚洲人与动物交配视频| 日产精品乱码卡一卡2卡三| 秋霞在线观看毛片| 免费人妻精品一区二区三区视频| 婷婷色av中文字幕| 欧美日韩视频精品一区| 亚洲精品国产av蜜桃| 99热国产这里只有精品6| 女的被弄到高潮叫床怎么办| 精品国产一区二区久久| 丁香六月天网| 大码成人一级视频| 一级毛片aaaaaa免费看小| 免费久久久久久久精品成人欧美视频 | 我要看黄色一级片免费的| 肉色欧美久久久久久久蜜桃| 在线 av 中文字幕| 亚洲精品视频女| 一边摸一边做爽爽视频免费| 欧美日韩视频精品一区| 99热这里只有精品一区| 女人精品久久久久毛片| 国产色婷婷99| av一本久久久久| 日韩在线高清观看一区二区三区| 国产成人av激情在线播放 | av一本久久久久| 国产不卡av网站在线观看| 91久久精品国产一区二区三区| 日韩在线高清观看一区二区三区| 观看美女的网站| 精品人妻熟女毛片av久久网站| 观看美女的网站| 免费大片黄手机在线观看| 久久 成人 亚洲| 国产免费现黄频在线看| 免费观看的影片在线观看| 寂寞人妻少妇视频99o| a级毛片黄视频| 22中文网久久字幕| 一级黄片播放器| 黄片无遮挡物在线观看| 欧美日韩亚洲高清精品| 国产成人免费观看mmmm| 大又大粗又爽又黄少妇毛片口| 午夜久久久在线观看| 亚洲人成网站在线观看播放| 久久午夜综合久久蜜桃| 国产欧美亚洲国产| 成年美女黄网站色视频大全免费 | 免费高清在线观看视频在线观看| 91精品国产九色| 国产国语露脸激情在线看| 国产乱人偷精品视频| 久久久久视频综合| 久久久精品免费免费高清| 久久ye,这里只有精品| 久久av网站| 中文字幕av电影在线播放| 久久国产精品大桥未久av| 狂野欧美激情性xxxx在线观看| 久久久久久人妻| 91精品国产九色| 国产精品国产av在线观看| 国产午夜精品久久久久久一区二区三区| 大片电影免费在线观看免费| 亚洲精品中文字幕在线视频| 国产成人精品在线电影| 777米奇影视久久| 午夜激情av网站| kizo精华| 91久久精品电影网| 国产欧美另类精品又又久久亚洲欧美| 狂野欧美激情性bbbbbb| 久久久久国产网址| 中国三级夫妇交换| 麻豆乱淫一区二区| 3wmmmm亚洲av在线观看| 丁香六月天网| 男女国产视频网站| 18禁裸乳无遮挡动漫免费视频| 精品人妻一区二区三区麻豆| 国产亚洲欧美精品永久| 热99久久久久精品小说推荐| 日本av手机在线免费观看| 国产免费现黄频在线看| 国产亚洲精品久久久com| av网站免费在线观看视频| 狂野欧美白嫩少妇大欣赏| 亚洲精品,欧美精品| 99热这里只有精品一区| 自拍欧美九色日韩亚洲蝌蚪91| 蜜桃在线观看..| 国产av国产精品国产| 国产熟女欧美一区二区| 午夜激情福利司机影院| 日韩精品有码人妻一区| 大又大粗又爽又黄少妇毛片口| 制服人妻中文乱码| 久久精品人人爽人人爽视色| av女优亚洲男人天堂| 只有这里有精品99| 伊人久久国产一区二区| 免费av中文字幕在线| 亚洲在久久综合| 最后的刺客免费高清国语| 母亲3免费完整高清在线观看 | 我的女老师完整版在线观看| 考比视频在线观看| 日日摸夜夜添夜夜爱| 夜夜看夜夜爽夜夜摸| 亚洲精品第二区| 国产熟女午夜一区二区三区 | av天堂久久9| 女人久久www免费人成看片| 99久国产av精品国产电影| 亚洲av中文av极速乱| 久久久久久人妻| 午夜激情久久久久久久| 高清视频免费观看一区二区| 精品久久久精品久久久| 亚洲精品乱码久久久久久按摩| 久久久久久久大尺度免费视频| 免费不卡的大黄色大毛片视频在线观看| 国产av一区二区精品久久| av国产精品久久久久影院| av.在线天堂| 国产日韩欧美视频二区| 搡女人真爽免费视频火全软件| 日韩电影二区| 日日爽夜夜爽网站| 成人毛片60女人毛片免费| 另类精品久久| 亚州av有码| 尾随美女入室| 婷婷色综合www| 欧美日韩成人在线一区二区| 十八禁高潮呻吟视频| 国国产精品蜜臀av免费| 久久久久久久久久久免费av| 黄色视频在线播放观看不卡| av福利片在线| 一级毛片aaaaaa免费看小| 亚洲精品久久午夜乱码| 亚洲欧美精品自产自拍| 嘟嘟电影网在线观看| 爱豆传媒免费全集在线观看| 久久久久久久久久成人| 欧美 亚洲 国产 日韩一| 黄片无遮挡物在线观看| 国产高清三级在线| 亚洲国产精品一区三区| 国产有黄有色有爽视频| 性色av一级| av免费在线看不卡| 国产伦理片在线播放av一区| 亚洲欧美清纯卡通| 中文欧美无线码| 少妇精品久久久久久久| 黑人欧美特级aaaaaa片| 99国产综合亚洲精品| 三级国产精品片| 99久久中文字幕三级久久日本| 国产精品人妻久久久久久| 99九九在线精品视频| 亚洲一区二区三区欧美精品| 亚洲少妇的诱惑av| 丰满迷人的少妇在线观看| 国产精品久久久久久久电影| 国产69精品久久久久777片| 免费不卡的大黄色大毛片视频在线观看| 日日摸夜夜添夜夜添av毛片| 亚洲综合精品二区| 街头女战士在线观看网站| 精品亚洲成国产av| 人人妻人人添人人爽欧美一区卜| 亚洲av日韩在线播放| 91精品国产九色| 高清av免费在线| av有码第一页| 免费黄色在线免费观看| 桃花免费在线播放| 伦精品一区二区三区| 2021少妇久久久久久久久久久| 免费观看a级毛片全部| 日韩电影二区| 亚洲成人手机| 一区二区三区精品91| 丝袜喷水一区| 黑人欧美特级aaaaaa片| 亚洲av日韩在线播放| 中文字幕免费在线视频6| 欧美少妇被猛烈插入视频| 色网站视频免费| 久久精品久久久久久久性| 大码成人一级视频| 黄色配什么色好看| 欧美xxxx性猛交bbbb| 国产成人91sexporn| 国产成人一区二区在线| 日韩亚洲欧美综合| 久久久久精品久久久久真实原创| 国产成人午夜福利电影在线观看| 亚洲欧美清纯卡通| 欧美日韩精品成人综合77777| 99久久精品一区二区三区| 精品国产国语对白av| 国产精品偷伦视频观看了| 久久国产精品男人的天堂亚洲 | 亚洲av在线观看美女高潮| 丰满迷人的少妇在线观看| 亚洲欧美一区二区三区黑人 | 国产69精品久久久久777片| 亚洲av福利一区| 国产精品.久久久| 秋霞在线观看毛片| 黑人高潮一二区| 五月开心婷婷网| 三级国产精品片| 中文精品一卡2卡3卡4更新| 看十八女毛片水多多多| 麻豆精品久久久久久蜜桃| 精品久久蜜臀av无| 性高湖久久久久久久久免费观看| 一边摸一边做爽爽视频免费| 七月丁香在线播放| 国产成人freesex在线| 在线 av 中文字幕| 中文欧美无线码| 精品人妻一区二区三区麻豆| 日韩强制内射视频| 午夜福利在线观看免费完整高清在| 亚洲成人一二三区av| av免费在线看不卡| 国产精品人妻久久久久久| 王馨瑶露胸无遮挡在线观看|