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

快凝Pd82Si18合金原子團簇的演化特性及遺傳機制*

2020-02-28 10:58:02高明鄧永和文大東田澤安趙鶴平彭平
物理學報 2020年4期
關鍵詞:結構

高明 鄧永和 文大東 田澤安 趙鶴平 彭平

1) (吉首大學物理與機電工程學院, 吉首 416000)

2) (湖南工程學院理學院, 湘潭 411104)

3) (湖南大學材料科學與工程學院, 長沙 410083)

采用分子動力學(MD)模擬計算, 對Pd82Si18合金快凝過程中基本原子團簇的遺傳特性、演化趨勢和結構穩定性進行了研究.團簇類型指數法(CTIM)分析表明: 非晶固體中Si原子為中心的(10 2/1441 8/1551)雙帽阿基米德反棱柱(BSAP)團簇數目占據優勢.快凝過程中, BSAP結構團簇具有最大的遺傳分數, 并且其他以Si原子為中心的Kasper團簇大多都會向BSAP結構團簇轉變.通過對Si原子為中心的Kasper基本團簇電子性質第一性原理計算發現, 體系中BSAP團簇的結合能最低, 結構穩定性較高, 與分子動力學計算結果一致.

1 引 言

Pd基非晶合金由于在力學、磁學、化學催化和耐腐蝕等方面都具備良好的性能而備受關注.自20世紀80年代Kui等[1]僅在冷速為10 K/S的條件下制備出厘米級Pd?Ni?P棒狀非晶合金以來, 實驗室可制備出最大的Pd基非晶合金樣品的直徑已經達到了75 mm, Pd基非晶合金被認為是具有最好玻璃形成能力(GFA)的合金體系[2].非晶合金的優異性能和GFA與其局域原子結構密切相關, 已成為非晶材料研究者的共識.已有研究發現Cu原子為中心的二十面體團簇的比例可以表征Cu?Zr快凝合金的玻璃轉變[3,4], 甚至可以用來評估非晶合金的GFA[5,6].Wu等[7]還發現過冷液體的二十面體中程序(IMRO)具有較長的壽命, 從而抑制結晶化, 促進金屬玻璃的形成.由于金屬玻璃被稱為“凍結的液體”, 本課題組從金屬玻璃液體和固體局域原子結構的關聯上找到了一條理解GFA的新途徑.鄧永和等[6]和Cheng等[8]采用分子動力學模擬發現: 快凝Cu?Zr合金二十面體團簇的遺傳性是快凝TM?TM合金的固有屬性, 與GFA有著密切的聯系, 二十面體團簇的遺傳分數和遺傳起始溫度都可以用來表征Cu?Zr非晶合金GFA的大小.在過渡金屬?金屬 (TM?M)的 Pd?Si非晶合金系統中, 雙帽阿基米德反棱柱(BSAP)和三帽三棱柱(TTP)的數目較多, 其中尤以BSAP為特征結構[9?13].然而, 針對TM?M合金系列在快凝形成非晶過程中的結構演化、遺傳特性和遺傳機制的研究還很少, 因此本文將針對結構的演化和遺傳特性以及產生結構遺傳的機制進行深入的研究.

由于在共晶點x = 17.2附近的Pd100?xSix合金都具有很強的GFA[14], 特別是Pd82Si18合金, 因此, 本文以Pd82Si18合金為研究對象來認識其特征團簇的結構遺傳特性及電子機制.本文首先采用分子動力學(MD)方法對Pd82Si18合金的快凝過程進行研究, 并采用團簇類型指數法(CTIM)表征和跟蹤其微結構的演化過程, 接著通過第一性原理,計算各團簇的能量分布, 找出能量較低的穩定團簇, 并分析結構的遺傳機制與團簇的演化過程和結構穩定性的關系, 發現結合能與結構的遺傳性有著較好的對應關系.

2 計算條件和方法

MD模擬采用 LAMMPS程序[15].首先將32000個原子(2624個Pd原子和5760個Si原子)隨機放入立方盒中, 同時給每個原子編號設置ID,采用Sheng等為Pd?Si合金發展的嵌入原子勢(EAM)[16], 通過將g(r)曲線理論計算值與實驗值作對比論證了該勢函數的準確性.在NPT系綜和周期性邊界條件下, 壓強保持為0, 步長為2 fs.首先在溫度為2100 K下等溫運行1 × 1 06步, 使體系熔化并達到平衡狀態, 接著以1 × 1 011K/s的冷速快凝至300 K.第一性原理計算基于密度泛函理論(DFT)的DMOL3軟件包[17], 在過冷液相區溫度810 K和玻璃態溫度300 K下, 對Si為中心的基本團簇進行幾何優化和總能計算的時候, 有效勢采用相對論修正的Effective Core Potentials(ECP)贗勢, 原子波函數采用帶一個p軌道極化函數的雙數值基組(d?polarization functions)[18],電子交換關聯勢選取GGA (general gradient approximate)近 似 的 Perdw?Burke?Emzerhof(PBE)[19]交換關聯泛函, 能量偏差小于 1.0 × 10—5Ha, 應力改變小于0.002 Ha/?、位移偏差小于0.005 ?, 平面波自洽場的迭代誤差為 1.0 × 10—6Ha,熱展寬(Smearing因子)為0.005 Ha, 在對弛豫前后團簇的能量進行計算時未采用任何對稱性限制.

3 模擬結果與討論

3.1 雙體分布函數和勢能演化

雙體分布函數可以有效地描述液體、晶體和非晶的結構特征[20].由于Pd82Si18合金的熔點約為1088 K[14], 所以圖1給出了Pd82Si18合金從1300 K熔體到300 K固體快凝過程的雙體分布函數變化的曲線圖.從圖1(a)中可以看出Pd82Si18總雙體分布函數g(r)tot曲線并無長程峰, 且第二峰隨溫度的降低而升高, 最終分裂形成兩個次峰, 說明Pd82Si18合金在快凝過程中形成了非晶結構.從圖1(b)中可以進一步看出, Pd82Si18合金g(r)tot曲線的第一峰在形成玻璃態時也發生了劈裂, 結果與前期研究Pd80Si20非晶合金的結果相同[13], 但是對于第一峰劈裂形成的詳細解釋, 文獻[13]并沒有給出.

由于Pd82Si18合金是低溶質體系, 從圖1(c)中可以看出Si原子與Si原子的距離已經超出了以Si原子為中心團簇第一近鄰的距離, 因此以Si原子為中心的短程序團簇的配位原子都是Pd原子, 這一點與溶質?溶質規避準則[11]完全符合.因此在快凝過程中Pd82Si18中Si原子與Si原子部分雙體分布函數g(r)Si?Si對g(r)tot第一峰的貢獻要遠小于Pd原子與Si原子部分雙體分布函數g(r)Pd?Si和Pd原子與Pd原子部分雙體分布函數g(r)Pd?Pd, 即Si與Si原子的規避作用也對圖1(b)劈裂形成起到一定的作用; 通過圖1(d)中g(r)Pd?Si和g(r)Pd?Pd的第一峰放大圖可以看出, 在冷卻過程中兩峰的峰寬都在逐漸縮小, 進而導致了Pd?Si鍵和Pd?Pd鍵的鍵長分布更加集中了, 而兩峰的相交區域卻在不斷減小, 即隨溫度降低, 原子振動范圍的不斷減小, 使得中間范圍(2.5—2.7 ?)鍵長分布減少; 從圖1(d)中 g(r)Pd?Si和 g(r)Pd?Pd第一峰峰值隨溫度的降低而增加可以看出, 以Pd為中心的Pd配位原子和Si配位原子都在不斷增加,以及Si為中心的Pd配位原子也在不斷增加.這些原因最終導致了圖1(b)中Pd82Si18非晶合金g(r)tot曲線第一峰劈裂的形成.

為了得到較為精確的玻璃轉化溫度, 圖2進一步給出了Pd82Si18體系每個原子的勢能隨溫度的變化關系, 從圖2可以看出: 勢能隨溫度的變化曲線并沒有發生突變, 這也進一步說明了在快凝過程中并沒有形成晶體相; 在2100—1000 K以及750—300 K范圍內勢能和溫度呈現出近線性的關系, 而在750—300 K范圍內的斜率卻明顯變小了, 說明在1000 —750 K的溫度區間內體系發生了玻璃化轉變.通過對勢能和溫度曲線的插值和外推得出該體系的Tg≈ 785 K, 該值明顯比實驗值634 K[14]高出許多, 導致這種差異的主要原因是受到計算資源的限制, 在MD模擬中采用的冷速高達1 ×1011K/s, 很難模擬出體系在低冷速下的演化過程.

圖1 Pd82Si18在1300 →300 K快凝過程中體系的雙體分布函數(ΔT = 100 K) (a) g(r)tot; (b) g(r)tot的第一峰放大圖; (c) g(r)tot第二峰放大圖; (d) g(r)Pd?Si和 g(r)Pd?Pd 第一峰的放大圖Fig.1.Pair distribution functions g(r) for rapidly solidified of Pd82Si18 from 1300 to 300 K (ΔT =100 K): (a) The g(r)tot curve;(b) first peak zoom of g(r)tot curve; (c) second peak zoom of g(r)tot curve; (d) first peak zoom of g(r)Pd?Si and g(r)Pd?Pd curve.

圖2 Pd82Si18合金在快凝過程中體系中原子的勢能隨溫度的變化Fig.2.Average atomic potential energy of per atom in the simulated system as a function of temperature T during rapid solidification.

3.2 微結構分析

采用基于H?A鍵型指數[21]的擴展原子團簇類型指數[22,23]法(CTIM)來表征團簇的局域結構,形如Z表示局域短程序團簇中心原子的配位數, nHA表示各種H?A鍵型指數的數目, ijkl是H?A鍵型的類型.如二十面體團簇是由12個1551 H?A鍵型組成, CTIM指數為(12 12/1551), BSAP結構的CTIM指數則可相應地寫為(10 2/1441 8/1551),同樣三椎三棱柱(TTP)的CTIM指數為(9 3/1441 6/1551).圖3進一步給出了BSAP和TTP的結構示意圖.值得指出的是CTIM表征方法并沒有采用近似處理[23], 可以更深刻地認識非晶的微觀結構.MD模擬的快凝Pd82Si18合金體系中, 所有短程序團簇的種類和屬性都可以很方便地導出.低溫時, 盡管得到的團簇種類有100多種, 但是數量超過100的也只有40多種, 而在這40多種團簇中,Kasper團簇占據了主要分數.

如圖4所示, 與研究Pd80Si20體系[13]類似, 在快凝過程中, (12 12/1551)和它的變形結構(12 2/1441 8/1551 2/1661)的數量增加緩慢, 遠少于其他Kasper團簇, 因此, 二十面體并非該體系的特征結構.而隨著溫度降低, 標準Kasper團簇(見圖4(a))中的BSAP結構團簇(10 2/1441 8/1551)和(11 2/1441 8/1551 1/1661)以及變形Kasper團簇 (見圖4(b))中 (10 1/1441 5/1551 1/1541 3/1431)和(14 2/1441 8/1551 4/1661)的數目明顯增加了, 尤其在Tm— Tg的過冷液相區內.相比于其他Kasper團簇, 在小于Tg的溫區內, 隨著溫度降低, BSAP團簇不僅增速最快, 而且在低溫條件下數量最多, 因此, 在Pd82Si18合金的非晶形成過程中, BSAP團簇起到了重要的作用.通過對體系結構的觀察我們還發現配位數為9, 10和11的Kasper團簇都是以Si原子為中心, 且配位原子都是Pd原子.因此, 從圖4(a)和圖4 (b)中可以看出Si原子為中心的Kasper團簇在低溫時占據了體系的主要分數, 所以Si原子為中心的Kasper團簇是Pd82Si18非晶合金體系的基本團簇.

圖3 CTIM指數為(10 2/1441 8/1551)和(9 3/1441 6/1551)的BSAP和TTP的結構示意圖(紅色的球表示Si原子, 灰色球表示Pd原子)Fig.3.Schematic diagram of BSAP and TTP with CTIM index of (10 2/1441 8/1551) and (9 3/1441 6/1551) (Red ball denote Si atom and gray balls denote Pd atoms).

圖4 在快凝過程中Pd82Si18合金基本團簇的數量隨溫度的變化關系 (a)標準Kasper團簇; (b)變形的Kasper團簇Fig.4.The temperature dependence of the number of typic?al basic clusters in Pd82Si18 alloys: (a) Canonical Kasper clusters; (b) distorted Kasper clusters.

3.3 基本團簇的結構遺傳與演化

Si原子為中心的BSAP團簇在Pd82Si18合金的快凝過程中增速最快, 數目最多, 而且以Si原子為中心的Kasper團簇是該體系的基本團簇, 因此對基本Si原子為中心的團簇結構的遺傳和演化過程進行了深入的研究.當體系的溫度從 T1下降到T2(< T1), 基本團簇中心原子的種類和編號保持不變的條件下, 團簇的類型和殼層原子的種類、數目及原子編號均保持不變的演化模式, 為完全遺傳.如果只是團簇的構型和中心原子種類與編號不變, 而部分配位原子的種類或編號發生了變化, 則稱為核遺傳[22].在核遺傳的模式下, 團簇的成分(即化學序)可能會發生改變, 但結構仍保持不變,如圖5所示.

根據上述定義, 跟蹤分析了以BSAP為主的6種基本Si心團簇的遺傳性.CTIM指數分別是標準構型的(9 3/1441 6/1551), (10 2/1441 8/1551),(11 2/1441 8/1551 1/1661), 以及相應的變形結構(9 1/1441 4/1551 4/1431), (10 1/1441 5/1551 1/1541 3/1431)和(11 1/1441 6/1551 2/1541 2/1431).

圖6所示為非晶合金Pd82Si18從810 到300 K的遺傳分數, 這里,= P和C, 分別表示完全遺傳和核遺傳), NT1表示溫度為T1(> T2)時某類團簇的總數,表示該類團簇以第i種模式從 T1遺傳到 T2的數目, 遺傳分數 f =fp+fc.跟蹤分析了過冷液相區(Tm— Tg)到玻璃態中基本Si心團簇的遺傳過程, 從圖6可以看到: 810—300 K的快凝過程中, BSAP的遺傳分數最大, 遺傳性能最好, 其次是(11 2/1441 8/1551 1/1661),遺傳性能最差的是(11 2/1441 8/1551 1/1661)的變形結構(11 1/1441 6/1551 2/1541 2/1431).

圖5 BSAP基本團簇遺傳示意圖 (a)完全遺傳; (b)核遺傳Fig.5.Basic cluster heredity schematic map of BSAP: (a) Perfect heredity; (b) core heredity.

為了考察810 到300 K快凝過程中這些團簇的相對穩定性, 對演化分數進行了詳細的統計.這里, 定義演化分數的計算方式為fj→i表示由第j類團簇演化為第i類團簇的演化分數,表示第j類團簇從高溫T1(810 K)到低溫T2(300 K)保持核原子的ID不變, 演化為第i類團簇的數目,表示 T1(810 K)溫度時, 第j類團簇的總數.計算結果如表1所列.通過觀察可以發現, 表1的演化分數與圖4相應的團簇在快凝過程中的增加趨勢表現出很好的一致性, 在表1的演化分數的統計中, 除了自身遺傳之外, 其他類型的團簇向BSAP演化的分數最多, 其次是(11 2/1441 8/1551 1/1661), 再者是(10 1/1441 5/1551 1/1541 3/1421), 而向 (9 3/1441 6/1551)團簇的演化相對較少, 向(11 1/1441 6/1551 2/1541 2/1431)團簇的演化分數最少.聯系圖6的遺傳性能,BSAP最好, 其次是(11 2/1441 8/1551 1/1661),(9 3/1441 6/1551)在結構遺傳方面相對較差, 最差的是(11 1/1441 6/1551 2/1541 2/1431).因此,我們得出結論: BSAP的結構遺傳及其他團簇向BSAP團簇的演化對Pd82Si18合金的GFA起重要作用.

圖6 非晶合金Pd82Si18從810 K到300 K的遺傳分數Fig.6.The heredity fractions in amorphous alloy Pd82Si18 from 810 K to 300 K.

3.4 基本團簇遺傳的電子機制分析

為了深入理解不同團簇遺傳產生差異的原因及其對GFA的影響機制, 進一步研究以Si原子為中心的Kasper團簇的能態和電子結構.結合能是度量結構穩定性的重要物理量, 采用第一性原理的方法, 分別對810和300 K的5100個以Si原子為中心Pd原子為殼層原子的Kasper團簇的結合能進行計算, 并做出統計, 結合能的定義為[24]

其中 n1和 n2分別表示每個Kasper團簇中Pd原子和Si原子的數目, E (Pdn1Sin2) 表示單個Kasper團簇的總能, EPd和 ESi分別表示單個Pd原子和Si原子的總能.

在810 K時(圖7(a)), Pd82Si18合金快凝至過冷液體的時候, BSAP和TTP總體的結合能分布趨勢大致相同, 難以區分出哪種團簇的結構穩定性最好, 但是相比于其他四種基本Si為中心的團簇,BSAP的結構穩定性更好.從圖7(b)中可以看到,300 K的Pd82Si18玻璃合金中, BSAP團簇的結合能分布在能量較低區間的團簇所占比例最高, 結構穩定性相對最好, 而TTP團簇的結合能分布的峰位處于相對較高的能態, 結構穩定性不如BSAP.圖7(c)進一步給出了不同類型Kasper團簇的平均結合能, 不難看出, 從 810到300 K的快凝過程中, 各類Kasper團簇的平均結合能都在降低, 穩定性都在增強.從810 到 300 K降溫過程中,Pd82Si18體系一直在不停地從高能態向低能態弛豫, 團簇結構也隨著一起弛豫.但無論是810 還是300 K, BSAP團簇始終呈現較低的平均結合能,300 K時表現得最低, 而(11 1/1441 6/1551 2/1541 2/1431)結構團簇的結合能分布和平均結合能一直較高, 很好地解釋BSAP團簇遺傳性能力最強, (11 1/1441 6/1551 2/1541 2/1431)結構團簇的遺傳能力最弱, 以及在810到300 K快凝過程中, 其他Kasper團簇大多數向BSAP演化, 而幾乎不向(11 1/1441 6/1551 2/1541 2/1431)結構團簇演化的原因.

表1 Pd82Si18合金從810 到300 K的幾種基本Si原子為中心的團簇的演化分數Table 1.The evolution fractions of several basic Si?centered clusters in amorphous alloy Pd82Si18 from 810 to 300 K.

圖7 非晶合金Pd82Si18在810 K和300 K的幾種基本Si為中心的團簇的結合能隨團簇的分布 (a) 810 K基本Si為中心的團簇的結合能分布; (b) 300 K基本Si為中心的團簇的結合能分布; (c) 810 與300 K基本Si為中心的團簇的平均結合能分布Fig.7.Binding energies of several basic Si?centered clusters of amorphous alloy Pd82Si18 at 810 and 300 K depend on the distribu?tion of clusters: (a) Binding energy distribution of basic Si?centered clusters at 800 K; (B) binding energy distribution of basic Si?centered clusters at 300 K; (c) distribution of average binding energy of basic Si?centered clusters at 800 and 300 K.

為了消除系統中團簇結構受溫度的影響, 本文又分別通過EAM勢和第一性原理計算了Pd82Si18快凝合金中以Si原子為中心Pd原子為殼層的配位數分別為9, 10和11的Kasper團簇優化后的結合能, 如圖8所示, Pd10Si的結合能最小, Pd9Si最大, 進一步證實了Pd82Si18快凝合金中10配位BSAP團簇確實相對其他Kasper團簇在理論上具有最低的結合能, 具有最高的遺傳分數, 其他團簇向其演化的分數最高.而11配位的(11 2/1441 8/1551 1/1661)結構團簇的結合能比BSAP略高,在Tg以下的快凝合金中數目僅次于BSAP.由于9配位的TTP團簇的結合能較高, 其在快凝合金中的數目也相對較少, 遺傳分數也較低.圖8(a)通過EAM計算所得的結合能趨勢與圖8(b)圖中第一性原理計算結果的基本一致.為了進一步說明TTP和BSAP團簇結合能的差異, 圖9進一步畫出了Pd9Si和Pd10Si優化后結構的局域電荷密度分布.由圖9可知: 殼層Pd與Pd原子之間以金屬鍵結合為主, 電荷密度為0.3處凸起的等高線部分存在著微弱的共價相互作用, 而Pd與Si原子之間則是以離子鍵的相互作用為主, 僅在電荷密度為0.5處凸起的部分存在著微弱的共價相互作用.可以看出Pd10Si團簇的殼層Pd原子與中心的Si原子間的電荷密度重疊比Pd9Si要高.圖10進一步給出了PdnSi (n = 9, 10, 11)團簇的電子態密度(DOS)圖.從圖10(b)中可以看出Pd9Si在費米能級上的電子數最少, 化學穩定性最好, 這與Cheng等[25]的結論非常一致.可以看出Kasper團簇的結合能對其遺傳性起決定作用[26,27].

圖8 基本Si為中心的團簇優化后結構的結合能隨團簇的分布 (a) EAM計算; (b)第一性原理計算Fig.8.Binding energies of several optimized basic Si?centered clusters depend on the distribution of clusters: (a) EAM calculations;(b) first?principle calculations.

圖9 局域電荷密度分布圖 (a)Si原子為中心的Pd10Si團簇的局域電荷密度; (b)Si原子為中心的Pd9Si團簇的局域電子密度(圖中白色和紅色的字體表示切面上的原子)Fig.9.Pattern of local charge density distribution: (a) Local charge density of Si?centered Pd10Si cluster; (b) local charge density of Si?centered Pd9Si cluster(White and red fonts in the figure represents atoms on the tangent plane).

圖10 優化后基本Si原子為中心的團簇的態密度(DOS)圖 (a) Pd9Si, Pd10Si與Pd11Si團簇的DOS圖; (b) 圖(a)中費米能級附近的放大圖Fig.10.The density of states (DOS) diagrams of optimized basic Si?centered clusters: (a) The DOS of Pd9Si、 Pd10Si and Pd11Si clusters; (b) zoom of the Fermi level in (a) diagram.

4 結 論

采用MD模擬對Pd82Si18合金快凝過程進行了研究和采用第一性原理計算對特征Kasper團簇的能態和電子結構進行分析, 得到以下主要結論:

1)由于Si原子的濃度低, Si原子與Si原子的排斥作用強, 在 Pd82Si18合金 Pd?Si鍵與 Pd?Pd鍵的鍵長分布中重疊區域隨著凝固不斷減小,導致了系統總雙體分布函數g(r)tot第一峰出現顯著劈裂, 表明Pd82Si18合金具有很好的非晶特性.

2)快凝過程中Si原子為中心的Kasper團簇對Pd82Si18合金非晶的形成起重要作用, 其中BSAP團簇在玻璃固體中數目最多, 在Tm—Tg過冷液相區增速最快, 并且從過冷液體遺傳至非晶固體中的比例最高, 是Pd82Si18快凝合金的特征團簇結構.

3)以Si原子為中心的Kasper團簇高的結構遺傳性對應低的結合能, Kasper團簇的結合能對其遺傳性起決定作用.

猜你喜歡
結構
DNA結構的發現
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
新型平衡塊結構的應用
模具制造(2019年3期)2019-06-06 02:10:54
循環結構謹防“死循環”
論《日出》的結構
縱向結構
縱向結構
我國社會結構的重建
人間(2015年21期)2015-03-11 15:23:21
創新治理結構促進中小企業持續成長
現代企業(2015年9期)2015-02-28 18:56:50
主站蜘蛛池模板: 91精品国产自产在线老师啪l| 国产综合网站| 中日韩一区二区三区中文免费视频| 久久精品娱乐亚洲领先| 亚洲欧美在线精品一区二区| 日韩国产精品无码一区二区三区| 成人精品午夜福利在线播放| 免费高清自慰一区二区三区| 国产精品成| 经典三级久久| 全部免费毛片免费播放 | 天天色天天综合| 亚洲精品无码久久毛片波多野吉| 国产亚洲精品97在线观看| 日韩在线永久免费播放| 国产内射一区亚洲| 污视频日本| 亚洲性影院| 男女男精品视频| 97在线碰| 日本成人一区| 欧美精品亚洲精品日韩专区| 亚洲精品无码成人片在线观看| 日本一本在线视频| 国产肉感大码AV无码| 天堂av综合网| 欧美成人a∨视频免费观看 | 欧美视频在线播放观看免费福利资源 | 成人毛片免费在线观看| 久久精品66| 国产精品亚洲专区一区| 亚洲精品桃花岛av在线| 国产激情无码一区二区免费| 高清无码一本到东京热| 欧美国产精品拍自| 亚洲综合在线网| 亚洲第七页| 在线视频97| 久久精品人妻中文系列| 98超碰在线观看| 久热中文字幕在线| 亚洲免费三区| 国产成人欧美| 国产一区二区精品高清在线观看| 亚洲最新在线| 538国产在线| 女人爽到高潮免费视频大全| 日韩福利视频导航| 欧美成人国产| 伊人久久大线影院首页| 9cao视频精品| 成人在线欧美| 日韩免费毛片| 国产成人综合久久精品尤物| 538精品在线观看| 日本福利视频网站| 亚洲日韩高清在线亚洲专区| 国产99久久亚洲综合精品西瓜tv| 久久国产拍爱| 成人国产免费| A级全黄试看30分钟小视频| 午夜福利视频一区| 中文字幕一区二区人妻电影| 国产农村精品一级毛片视频| 亚洲动漫h| 无码福利日韩神码福利片| 免费看久久精品99| 乱人伦视频中文字幕在线| 欧美成在线视频| 国产幂在线无码精品| 噜噜噜久久| 狠狠色狠狠色综合久久第一次| 伊人网址在线| 天堂va亚洲va欧美va国产| 国产在线日本| www亚洲精品| 国产欧美另类| 麻豆AV网站免费进入| 91av国产在线| 亚洲第一在线播放| 国产欧美精品一区二区| AV不卡在线永久免费观看|