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

Fcc Pt的熱力學性質和熱物理性質的CALPHAD研究

2016-09-05 06:37:52徐志鋒郭鑫魯曉剛
上海金屬 2016年6期
關鍵詞:實驗模型研究

徐志鋒郭 鑫魯曉剛,2

(1.上海大學材料科學與工程學院,上海 200444;2.上海大學材料基因組工程研究院,上海 200444)

Fcc Pt的熱力學性質和熱物理性質的CALPHAD研究

徐志鋒1郭 鑫1魯曉剛1,2

(1.上海大學材料科學與工程學院,上海 200444;2.上海大學材料基因組工程研究院,上海 200444)

應用CALPHAD方法和改進的Debye-Grüneisen模型對具有fcc(面心立方)結構的純Pt的熱力學性質和熱物理性質進行了優(yōu)化研究。該方法利用一套具有實際物理意義的模型參數(shù),可以準確地描述fcc Pt在寬廣溫度(0 K~熔點)和壓力(常壓~極端高壓)范圍內(nèi)的大部分熱力學性質(包括Gibbs自由能、熵和熱容等)和熱物理性質(包括體積、熱膨脹系數(shù)和彈性模量等),揭示了熱力學性質和熱物理性質的內(nèi)在聯(lián)系。

面心立方Pt CALPHAD Debye-Grüneisen模型 熱力學性質 熱物理性質

經(jīng)過數(shù)十年的發(fā)展,CALPHAD方法在材料科學研究和工程應用上受到越來越多的關注,已經(jīng)成為一種成熟、強大的材料熱力學和擴散動力學計算模擬技術。在傳統(tǒng)CALPHAD方法中,Gibbs自由能是以溫度T和壓強P作為變量,而實驗數(shù)據(jù)一般都是在常壓下測得的,因此目前絕大多數(shù)的CALPHAD數(shù)據(jù)只是常壓數(shù)據(jù)。為了描述壓強的貢獻,需要確定壓力、溫度、體積以及其他相關參數(shù)之間的關系,即確定狀態(tài)方程(Equation of State,EOS)。具體地說,就是采用壓力模型,如Morse模型[1]、Murnaghan模型[2]、Birch-Murnaghan模型[3]、Vinet模型[4]等,依據(jù)高壓范圍的實驗數(shù)據(jù)(包括熱力學和相圖數(shù)據(jù)),優(yōu)化擬合壓力模型中的所有參數(shù)。一般來說,壓力模型中的參數(shù)擬合和常壓下熱力學數(shù)據(jù)的確定是分開進行的。

用傳統(tǒng)CALPHAD建模方法在描述包含熵、焓、熱容和化學勢等熱力學性質,以及體積、熱膨脹系數(shù)和彈性模量等熱物理性質時,都是利用數(shù)學多項式獨立進行的,因此不得不各用一套參數(shù)來擬合。這些參數(shù)沒有實際物理意義,各參數(shù)之間也沒有直接聯(lián)系,擬合過程中也互無約束。而事實上材料的各種性質間存在著緊密的內(nèi)在相互約束關系。

而在第一性原理計算等理論研究中,由于計算晶體的聲子譜隨溫度的變化非常困難,故常常采用所謂的準簡諧近似,即在不同的體積下計算聲子譜,將其視為體積的函數(shù),再積分求得晶格振動產(chǎn)生的能量。如此,以溫度T和體積V作為變量的Helmholtz自由能就是理論計算中一個很自然的建模選擇。計算表明,這種將各種物理現(xiàn)象或過程對自由能的貢獻分開描述的建模方法,能夠保持各種性質間內(nèi)在的相互約束關系。

本文運用Helmholtz自由能結合Debye-Grüneisen模型,用較少的具有實際物理意義的參數(shù),統(tǒng)一優(yōu)化擬合了fcc結構的純Pt在常壓及高壓,從絕對零度到純Pt熔點的包括Gibbs自由能和熱容等在內(nèi)的熱力學,以及包括體積、熱膨脹系數(shù)和彈性模量等在內(nèi)的熱物理性質的實驗數(shù)據(jù)。這種新方法得到的數(shù)據(jù)不但更具有物理意義,而且僅使用一套參數(shù)便可計算各類熱力學性質和熱物理性質。

材料熱力學和熱物理性質等是材料設計必不可少的基礎數(shù)據(jù)。本研究的模型參數(shù)是通過對大量實驗數(shù)據(jù)進行計算軟件的分析、甄別、優(yōu)化、擬合得到的,在材料設計基礎數(shù)據(jù)庫中只需存儲這些關鍵的模型參數(shù),經(jīng)過計算模擬,就可以獲得任意設定條件下的材料熱力學及熱物理性質,這對材料的理性設計具有重要的指導意義。

1 研究方法

與傳統(tǒng)的CALPHAD Gibbs自由能方法不同,本研究所采用的是基于Helmholtz自由能結合Debye-Grüneisen模型和自由電子Fermi氣模型的方法。某一體系的Helmholtz自由能往往可以表達為各種物理現(xiàn)象或過程的貢獻之和。總的Helmholtz自由能F主要由以下幾個部分構成:0 K下的總能量(即靜態(tài)晶格能)Etot、晶格振動能FD、電子熱激發(fā)能Fel和磁性貢獻的能量Fmag。對fcc Pt來說,磁性貢獻的能量Fmag非常小,可以忽略不計,則有:

本研究中,0 K下的總能量Etot可以用Morse方程來表達[1]:

式中,x是晶格常數(shù),φ是擬合參數(shù),Ec是內(nèi)聚能,Eref是參考能,Eref的取值要保證體系在298.15 K和1 bar條件下的焓值為0。

體系的晶格振動能FD可以用Debye模型或Einstein模型[5]描述。為了準確描述寬廣溫度和壓力下的各種熱力學性質和熱物理性質,本研究采用了改進的Debye模型[6],即模型中的關鍵參數(shù),Debye溫度θD,寫成如式(3)體積V的函數(shù):

式中,h、kB和Na分別代表普朗克常數(shù)、波爾茲曼常數(shù)和阿伏伽德羅常數(shù)。m為體系的質量,ν是泊松比,k()ν是一個與泊松比ν相關的函數(shù):

模型中的另一個重要參數(shù)就是Grüneisen參數(shù),一般來說,有三種近似用來描述Grüneisen參數(shù):當式(3)中的λ=-1、0和+1時,分別對應Slater近似[7]、Dugdale-MacDonald(DM)近似[8]和自由體積理論近似[9]。

另一方面,電子熱容和電子貢獻的能量也可以用Sommerfeld參數(shù)γel來描述:

為了與體積V聯(lián)系起來,將Sommerfeld參數(shù)γel寫成一個與溫度無關、而與體積相關的函數(shù):

2 結果與討論

近幾十年來,對純Pt的各種熱力學性質和熱物理性質的研究比較廣泛,這給本工作對fcc Pt進行綜合研究提供了良好的基礎。本研究優(yōu)化得到fcc Pt的模型參數(shù)值如表1所示。各熱力學性質和熱物理性質的優(yōu)化結果如圖1~圖7所示。

表1 本工作優(yōu)化得到的fcc Pt模型參數(shù)值Table 1 Model parameters for fcc Pt in the presentwork

2.1 熱容

通過文獻查閱收集到了大量的fcc Pt的熱容實驗數(shù)據(jù)。本研究選用了Berg[10]、Furukawa等[11]、Seville[12]、Yokokawa和Takahashi[13]和AIP[14]手冊上的數(shù)據(jù)進行優(yōu)化。Berg[10]實驗測量了純Pt在0~20 K的熱容,并利用這些極低溫度下的實驗數(shù)據(jù)研究了純Pt的Sommerfeld參數(shù)和0 K的Debye溫度。Furukawa等[11]測量了fcc Pt室溫以下的熱容。Yokokawa和Takahashi[13]用量熱法測得了純Pt在80~1 000 K的熱容值。高溫熱容值(600~1 850 K)本研究選用的是Seville[12]測量的結果。本工作計算了有電子貢獻的熱容和沒有電子貢獻的熱容,并在圖1中分別以實線和虛線表示。從圖1可以看出,本工作的計算結果與實驗數(shù)據(jù)吻合較好。

圖1 常壓時fcc Pt不同溫度下的熱容Fig.1 Calculated heat capacity at atmospheric pressure compared with the experimental data for fcc Pt

2.2 低溫電子熱容

極低溫度時,在不考慮相轉變和磁性影響的前提下,某一體系的熱容可以視為由兩部分組成:晶格振動的貢獻和電子的貢獻,對應的關系式如式(7):

式(7)中,(γelT+aTb)是電子貢獻部分,AT3是晶格振動的貢獻。由于在極低溫度下參數(shù)a和b的影響非常小,研究時常常可以忽略,因此式(7)就可以近似地表達為:

式(8)兩邊同時除以T,得:

分析某一體系在極低溫度下電子對熱容的貢獻常用的方法是作cp/T和T2的關系圖,如圖2所示。由式(9)可知,cp/T和T2呈一次函數(shù)的關系,圖中直線與Y軸的截距就是0 K時的Sommerfeld參數(shù)γel,直線的斜率就是常數(shù)A。本研究中,電子對熱容的貢獻是通過和c這4個參數(shù)來描述的。這4個參數(shù)的值是在其他參數(shù)都取得合理的數(shù)值的基礎上,利用極低溫度下的實驗數(shù)據(jù)優(yōu)化得到的。本研究選用的極低溫度下的實驗數(shù)據(jù)來源于Berg[10]、Furukawa等[11]和Boerstoel等15]。他們的實驗測量結果基本一致。本研究的計算結果與Berg[10]、Furukawa等[11]和Boerstoel等[12]的實驗數(shù)據(jù)吻合很好。

2.3 線性熱膨脹系數(shù)

圖3是fcc Pt的熱膨脹系數(shù)隨溫度的變化情況。熱膨脹系數(shù)的實驗數(shù)據(jù)很多,本研究采用了AIP[14]、Kirby[16]和TPRC[17]的實驗數(shù)據(jù)。Touloukian等[17]對1973年前的30套熱膨脹系數(shù)的實驗數(shù)據(jù)進行了分析甄別,列出了一套極具參考價值的數(shù)據(jù),即TPRC。AIP[14]全稱American Institute of Physics handbook,也是獲得學術界普遍認可的一套數(shù)據(jù)。如圖3所示,本研究的計算結果與實驗數(shù)據(jù)非常吻合。

圖2 常壓時fcc Pt極低溫度下的熱容Fig.2 Calculated heat capacity of fcc Pt at extremely low temperatures and atmospheric pressure

圖3 常壓時fcc Pt不同溫度的熱膨脹系數(shù)Fig.3 Calculated coefficient of linear thermal expansion coefficients at atmospheric pressure compared with the experimental data for fcc Pt

2.4 摩爾體積

圖4是fcc Pt在不同溫度下的摩爾體積。采用的實驗數(shù)據(jù)來源于Kirby[16]、Lu等[18]、Schroder等[19]、Evans和Fischer[20]、Edwards等[21]和Owen和Yates[22],摩爾體積的實驗數(shù)據(jù)一般是通過X射線技術測量晶格常數(shù),然后由晶格常數(shù)轉化得到。本研究的計算結果與這些實驗數(shù)據(jù)非常一致。

圖4 常壓時fcc Pt不同溫度的摩爾體積Fig.4 Calculated molar volume of fcc Pt at atmospheric pressure compared with the experimental data

2.5 體彈模量

一般來說,體彈模量的數(shù)據(jù)是通過實驗測得體系在不同溫度下的3個二階彈性常數(shù)c11、c12和c44的值,然后通過關系式計算得到。Collard和McLellan[23]實驗測得了fcc Pt在300~1 480 K的體彈模量,Dorgokupets和Oganov[24]的實驗結果在800 K以內(nèi)與Collard和McLellan[23]的數(shù)據(jù)基本一致,在800 K以上Dorgokupets和Oganov[24]的結果略小一些。本研究分別計算了絕熱體彈模量BS和等溫體彈模量BT,如圖5所示。由于缺乏等溫體彈模量BT的實驗數(shù)據(jù),本研究只采用了絕熱體彈模量BS的實驗數(shù)據(jù)。通過比較可以發(fā)現(xiàn),本研究的計算結果與實驗數(shù)據(jù)吻合得比較理想,在0~700 K偏差在1.2%以內(nèi),700 K以上偏差約為2%。對于實驗數(shù)據(jù)相對比較離散的彈性性質,本研究的計算結果與實驗數(shù)據(jù)相當吻合。

2.6 楊氏模量

當體彈模量B和泊松比ν確定之后,楊氏模量E就可以通過式(10)計算得到:

本工作計算得到的fcc Pt的楊氏模量如圖6所示,楊氏模量的實驗數(shù)據(jù)相對比較少,本研究選用的實驗數(shù)據(jù)來源于Collard和McLellan[23]。從圖中可以看出,本研究的計算結果與實驗數(shù)據(jù)吻合得很好,偏差在1%以內(nèi)。

圖5 常壓時fcc Pt不同溫度的體彈模量Fig.5 Calculated adiabatic(BS)and isothermal(BT)bulk modulus at atmospheric pressure compared with the experimental data for fcc Pt

圖6 常壓時fcc Pt不同溫度的楊氏模量Fig.6 Calculated Young’smodulus at atmospheric pressure compared with the experimental data for fcc Pt

2.7 高壓下的體積比

圖7為fcc Pt在300 K時,體積比V/V0隨壓強P的變化曲線。實驗數(shù)據(jù)來源于Holmes等[25]、McQueen等[26]、Matsui等[27]、Yokoo等[28]、SUN等[29]、Singh和Sabrawat[30]和Altshuler等[31]。實驗采用的測量方法一般都是沖擊波技術。Holmes等[25]利用該技術測量了fcc Pt高壓下的狀態(tài)方程(壓強高達660 GPa)。McQueen等[26]以同樣的方法對多種金屬元素進行了測試,測得壓強范圍為300 GPa的數(shù)據(jù)。Matsui等[27]測量了壓強范圍為290 GPa的體積比。Yokoo等[28]和SUN等[29]測得的最高壓強都是550 GPa。Singh和Sabrawat[30]采用幾個可靠的高壓狀態(tài)方程計算了fcc Pt高壓下的性質。Altshuler等[31]測得了25種金屬的高壓性質,最高壓強為400 GPa。本研究的計算結果與這些數(shù)據(jù)具有很好的一致性。

圖7 300 K時fcc Pt高壓條件下的體積比V/V0Fig.7 Normalized molar volume of fcc Pt at300 K under high pressures

3 結論

本研究采用CALPHAD Helmholtz自由能結合Debye-Grüneisen模型的方法對fcc Pt的熱力學性質和熱物理性質進行了優(yōu)化,得到了一套具有實際物理意義的模型參數(shù),準確地描述了fcc Pt在常壓及高壓,從絕對零度到純Pt熔點的各種熱力學性質和熱物理性質,并保持了各種性質間內(nèi)在的相互約束關系。

致謝:

感謝國家自然科學基金對本工作的資助(No.51271106)。

[1]MORUZZIV L,JANAK JF,SCHWARZ K.Calculated thermalproperties of metals[J].Physical Review B,1988,37(2):790-799.

[2]MURNAGHAN FD.The compressibility ofmedia under extreme pressures[J].Proceedings of the National Academy of Sciences,1944,30(9):244-247.

[3]BIRCH F.Elasticity and constitution of the Earth's interior[J].Journal of Geophysical Research,1952,57(2):227-286.

[4]VINET P,F(xiàn)ERRANTE J,ROSE JH,et al.Compressibility of solids[J].Journal of Geophysical Research:Solid Earth(1978-2012),1987,92(B9):9319-9325.

[5]GRIMVALL G.Thermophysical properties of materials[M].Amsterdam,The Netherlands:Elsevier,1999.

[6]LU X G,CHEN Q.A CALPHAD Helmholtz energy approach to calculate thermodynamic and thermophysical properties of fcc Cu[J].Philosophical Magazine,2009,89(25):2167-2194.

[7]SLATER JC.Introduction to chemical physics[M].Martindell Press,1939.

[8]DUGDALE JS,MACDONALDD K C.The thermalexpansion of solids[J].Physical Review,1953,89(4):832-834.

[9]VASHCHENKO V I A,ZUBAREV V N.Free volume theory used to derive an expression for the gruneisen constant gamma,including an analysis of the degree of approximation[J].Soviet Physics-Solid State,1963,5:653-655.

[10]BERGWT.The low temperature heat capacity of platinum[J].Journal of Physics and Chemistry of Solids,1969,30(1):69-72.

[11]FURUKAWA G T,REILLY ML,GALLAGHER JS.Critical Analysis of Heat-Capacity Data and Evaluation of Thermodynamic Properties of Ruthenium,Rhodium,Palladium,Iridium,and Platinum from 0 to 300K.A Survey of the Literature Data on Osmium[J].Journal of Physical and Chemical Reference Data,1974,3(1):163-209.

[12]SEVILLE A H.The heat capacity of platinum at high temperatures[J].The Journal of Chemical Thermodynamics,1975,7(4):383-387.

[13]YOKOKAWA H,TAKAHASHIY.Laser-flash calorimetry II.Heat capacity of platinum from 80 to 1000 K and its revised thermodynamic functions[J].The Journal of Chemical Thermodynamics,1979,11(5):411-420.

[14]GRAY D E.American institute of physics handbook[M].New York:McGraw-Hill,1982.

[15]BOERSTOEL BM,ZWART JJ,HANSEN J.The specific heat of palladium,platinum,gold and copper below 30 K[J].Physica,1971,54(3):442-458.

[16]KIRBY R K.Platinum—a thermal expansion referencematerial[J].International Journal of Thermophysics,1991,12(4):679-685.

[17]TOULOUKIAN Y S,KIRBY R K,TAYLOR R E,et al.Thermophysical Properties of Matter-the TPRC Data Series.Volume 12.Thermal Expansion Metallic Elements and Alloys[R].Thermophysical and Electronic Properties Information Analysis Center,1975.

[18]LU X G,SELLEBYM,SUNDMAN B.Theoreticalmodeling of molar volume and thermal expansion[J].Acta materialia,2005,53(8):2259-2272.

[19]SCHR?DER R H,SCHMITZ-PRANGHE N,KOHLHAAS R.Experimental determination of lattice parameters of platinum metals in the temperature range between minus 190 deg and 1709 C(Analysis of lattice spacings and changes in platinum metals as function of temperature and applied pressure)[J].Z.Metallk.,(Stuttgart),1972,63:12-16.

[20]EVANSD L,F(xiàn)ISCHERGR.The Determination of the Thermal Expansion of Platinum By X-Ray Diffraction[C]//Proceedings of the 1971 Thermal Expansion Symposium.AIP Publishing,1972,3(1):97-104.

[21]EDWARDS J W,SPEISER R,JOHNSTON H L.High Temperature Structure and Thermal Expansion of Some Metals as Determined by X-Ray Diffraction Data.I.Platinum,Tantalum,Niobium,and Molybdenum[J].Journal of Applied Physics,1951,22(4):424-428.

[22]OWEN E A,YATES E L.IX.The thermal expansion of the crystal lattices of silver,platinum,and zinc[J].The London,Edinburgh,and Dublin Philosophical Magazine and Journal of Science,1934,17(110):113-131.

[23]COLLARD SM,MCLELLAN R B.High-temperature elastic constants of platinum single crystals[J].Acta Metallurgica et Materialia,1992,40(4):699-702.

[24]DOROGOKUPETS P I,OGANOV A R.Ruby,metals,and MgO as alternative pressure scales:A semiempirical description of shock-wave,ultrasonic,x-ray,and thermochemical data at high temperaturesand pressures[J].Physical Review B,2007,75(2):024115(1-16).

[25]HOLMESN C,MORIARTY JA,GATHERSG R,et al.The equation of state of platinum to 660 GPa(6.6 Mbar)[J].Journal of Applied Physics,1989,66(7):2962-2967.

[26]MCQUEEN R G,MARSH SP.Equation of State for Nineteen Metallic Elements from Shock-Wave Measurements to Two Megabars[J].Journal of Applied Physics,1960,31(7):1253-1269.

[27]MATSUIM,ITO E,KATSURA T,et al.The temperaturepressure-volume equation of state of platinum[J].Journal of Applied Physics,2009,105(1):013505.

[28]YOKOO M,KAWAIN,NAKAMURA K G,et al.Ultrahighpressure scales for gold and platinum at pressuresup to550 GPa[J].Physical Review B,2009,80(10):104114.

[29]SUN T,UMEMOTO K,WU Z,et al.Lattice dynamics and thermal equation of state of platinum[J].Physical Review B,2008,78(2):024304.

[30]SINGH R S,SAHRAWAT D.Grüneisen Parameter and Its Higher Derivatives for Pt,F(xiàn)e,V and Nb Using Equations of State[J].2015,02(02):0199-0205.

[31]ALTSHULER L V,BRUSNIKIN S E,KUZ'MENKOV E A.Isotherms and Grüneisen functions for 25 metals[J].Journal of Applied Mechaics and Technical Physics,1987,28(1):129-141.

收修改稿日期:2016-01-20

Study on Thermodynam ic and Thermophysical Properties of fcc Pt through the CALPHAD Approach

Xu Zhifeng1Guo Xin1LuXiaogang1,2
(1.School of Materials Science and Engineering,Shanghai University,Shanghai200444,China;2.Materials Genome Institute,Shanghai University,Shanghai200444,China)

Assessment of thermodynamic and thermophysical properties of fcc Pt has been performed bymeans of CALPHAD on the basis of amodified Debye-Grüneisenmodel.This approach,even though employed only a few parameters thathave physicalmeaning,could reproduce notonly the experimental data of thermodynamic and thermophysical properties at ambient pressure,but also those at extremely high temperatures and pressures with remarkable accuracy.The internal relations between thermodynamic and thermophysical properties have been demonstrated in this research.

fcc Pt,CALPHAD,Debye-Grüneisen model,thermodynamic property,thermophysical property

國家自然科學基金(No.51271106)

徐志鋒,男,主要從事計算材料學研究,Email:17717375203@163.com

猜你喜歡
實驗模型研究
一半模型
記一次有趣的實驗
FMS與YBT相關性的實證研究
遼代千人邑研究述論
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
視錯覺在平面設計中的應用與研究
科技傳播(2019年22期)2020-01-14 03:06:54
做個怪怪長實驗
EMA伺服控制系統(tǒng)研究
3D打印中的模型分割與打包
主站蜘蛛池模板: 国产伦片中文免费观看| 亚洲国产无码有码| 丝袜久久剧情精品国产| 国产嫖妓91东北老熟女久久一| 日韩精品免费一线在线观看 | 国内熟女少妇一线天| 亚洲日韩在线满18点击进入| 91亚洲国产视频| 一级黄色网站在线免费看| 国产美女无遮挡免费视频网站| 国产在线自乱拍播放| 欧美成人精品欧美一级乱黄| 成人国产一区二区三区| 91精品情国产情侣高潮对白蜜| 欧洲高清无码在线| 国产精品人莉莉成在线播放| 国产精品蜜芽在线观看| 国产99精品久久| 成人小视频在线观看免费| 香蕉久久国产精品免| 亚洲欧洲综合| 伊人久久婷婷| 欧美日韩国产成人在线观看| 激情乱人伦| 亚洲成人黄色在线观看| 在线观看精品自拍视频| 综合久久五月天| 亚洲AV无码一区二区三区牲色| 中文字幕免费播放| 国产成人综合欧美精品久久| 中文字幕第1页在线播| 欧美www在线观看| 欧美日韩亚洲国产| 欧美一级片在线| 久久五月天综合| 亚洲天堂久久新| 国产自在线拍| 色综合婷婷| 久久久久久久久18禁秘| 99re视频在线| 天堂在线亚洲| 国产在线八区| 一级毛片免费的| 亚洲成人手机在线| 538精品在线观看| 欧美精品v欧洲精品| 欧亚日韩Av| 色视频国产| 伊人久久综在合线亚洲91| 国产精品xxx| 日本一本在线视频| 十八禁美女裸体网站| 无码福利日韩神码福利片| 国产又粗又爽视频| 亚洲国产亚综合在线区| 亚洲国内精品自在自线官| 亚洲国产成熟视频在线多多| 国产一区二区三区视频| a级毛片网| 极品私人尤物在线精品首页| 国产精品无码翘臀在线看纯欲| 在线观看亚洲天堂| 亚洲成人动漫在线观看| 91网红精品在线观看| 九九九国产| 欧美第二区| 国产精品久久久精品三级| 91无码网站| 91av国产在线| 国产aⅴ无码专区亚洲av综合网 | 无码在线激情片| 免费国产一级 片内射老| 久久国产成人精品国产成人亚洲 | 亚洲首页国产精品丝袜| 日韩欧美在线观看| 看看一级毛片| 99人体免费视频| 国产麻豆91网在线看| 青青极品在线| 欧美国产在线看| 亚洲欧美综合另类图片小说区| 青青极品在线|