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

介電陶瓷/NiZn鐵氧體互擴散行為的第一性原理研究

2023-04-29 00:44:03張凱郭子康劉振濤畢鵬
西南科技大學學報 2023年4期

張凱 郭子康 劉振濤 畢鵬

摘要:針對介電陶瓷/NiZn鐵氧體異質復合材料的低溫共燒陶瓷體系,建立摻雜結構模型,采用基于密度泛函理論的第一性原理計算并結合CI-NEB方法研究陽離子互擴散機制。結果表明:鐵氧體中的Ni, Zn, Fe離子主要取代介電陶瓷MgTiO3和CaTiO3體系中的Ti位,遷移勢壘1.0~5.5 eV;對于介電陶瓷Mg和Ca傾向于占據NiZnFe4O8中Zn位點,Ti則傾向于取代Fe位,遷移勢壘0.6~1.0 eV。對于該材料體系,典型共燒工藝條件下Ca, Mg擴散進入NiZn鐵氧體距離400~1 000 μm。

關鍵詞:介電陶瓷/鐵氧體共燒體系 摻雜 互擴散 第一性原理

中圖分類號:TQ174.1文獻標志碼:A文章編號:1671-8755(2023)04-0045-09

First-principles Study of Diffusion Behavior between Dielectric Ceramics and NiZn-ferrite

ZHANG Kai1, GUO Zikang1, LIU Zhentao1, BI Peng2

Abstract:? For the low temperature co-fired ceramic system of the dielectric ceramic and NiZn-ferrite, the doping model is established, and the cation mutual diffusion mechanism is studied by using the first-principles calculation based on the DFT and CI-NEB method. The results show that Ni, Zn and Fe in ferrite mainly replace the Ti in the dielectric ceramic MgTiO3 and CaTiO3, and the migration barrier is 1.0-5.5 eV; For dielectric ceramics, Mg and Ca tend to occupy the Zn site in NiZnFe4O8, while Ti tends to replace the Fe site, and the migration barrier is 0.6-1.0 eV. For this material system, under typical co-firing process conditions, Ca and Mg diffuse into NiZn ferrite at a distance 400-1 000 μm.

Keywords: Dielectric ceramic/ferrite co-firing system;? Doping;? Mutual diffusion;? First-principles

伴隨5G時代的到來,電子器件在片式小型化的同時朝著高性能、多功能、高可靠的方向發展。對于電子器件來說,單一種類的材料已很難滿足要求,需開發新型共燒復合材料[1-2],如:產業龍頭Skyworks推出的鐵氧體磁芯與介質基板高溫共燒集成微帶或基片集成波導環行器就是其典型應用。由于介質材料與磁性材料的晶體結構差異較大,界面處的相容性較差,共燒連接時存在著熱膨脹系數難匹配、燒結收縮率難調控、界面反應和擴散不可控等問題,兩相共燒之后由擴散形成的結構復雜界面可提高材料的力學和電磁性能。適當的擴散會在過渡區形成新相,并在界面處形成共價鍵,這將有利于環行器復合襯底中陶瓷與鐵氧體的緊密共燒[3]。但擴散過度或擴散區延伸過長,會導致A, B位鐵氧體引入非磁性離子或快弛豫離子,進而影響電阻率、飽和磁化強度和磁性鐵氧體溫度穩定性[4]。因此,研究界面反應和界面擴散機制對提高器件性能的穩定性和可靠性具有重要的指導意義。

本文擬采用第一性原理以及CI-NEB方法研究MgTiO3,CaTiO3介電陶瓷和NiZn鐵氧體異質復合材料共燒過程中的陽離子互擴散機制,并評估各種典型離子影響的界面長度,為后期成分設計與工藝優化提供理論依據。

1 建模與計算方法

1.1 模型建立

本文中所涉及的3種晶體分別為MgTiO3,CaTiO3以及NiZnFe4O8。對3種晶胞做充分的弛豫,使其達到最穩定狀態,弛豫后的晶體結構如圖1所示,空間群以及晶格常數列于表1,并給出了3種晶體的實驗晶格參數。計算得到的晶格參數與實驗值接近。

晶胞經充分弛豫后,使用1×1×1的MgTiO3和NiZnFe4O8晶胞以及2×2×2的CaTiO3超胞構建陽離子摻雜介電陶瓷和鐵氧體的晶體模型來模擬其界面擴散過程,并以此模型進一步研究陽離子在晶格中的擴散路徑與勢壘。對于NiZnFe4O8,其體內的陽離子主要來源為介電陶瓷MgTiO3和CaTiO3,即Mg,Ca和Ti;同理,進入介電陶瓷晶體內的陽離子為NiZnFe4O8中Ni,Zn和Fe。

陽離子摻雜主要為取代晶體內固有的陽離子位點O,每種陽離子有多個摻雜位點,如在NiZnFe4O8中有3種取代情況,分別是Ni,Zn所在的四面體位點以及Fe所在的八面體位點。同樣,在兩種介電陶瓷MgTiO3和CaTiO3晶體中都有兩種取代結構。構建的摻雜模型記為X[M],其中X表示取代的陽離子,而M記為被取代的陽離子,如Mg[Fe]則表示Mg取代NiZnFe4O8中的Fe。為了確定陽離子在晶體中最穩定的摻雜結構,使用摻雜形成能來評估摻雜模型的穩定性。摻雜形成能[8-9]的計算公式為:

式中:Ef為摻雜形成能;Ebulk和Edefect分別表示摻雜前晶胞和摻雜后晶胞的總能量;μ為單個原子的化學勢。摻雜缺陷的產生涉及向晶胞中替換原子,式中μY為摻雜原子Y的化學勢,μX為被替換原子X的化學勢。原子化學勢μ的大小與原子所處的化學環境有關,在介電陶瓷和NiZn鐵氧體共燒過程中,體系氧原子占多數,可認為是富氧環境。因此在富氧環境下有:

式中:μM為金屬的化學勢;MaOb為穩定金屬氧化物。計算得到μM和μO的值列于表2。

1.2 計算方法

所有的計算均使用基于密度泛函理論框架下的計算軟件VASP(Vienna ab initio Simulation Package)完成[10-11]。在廣義梯度近似(GGA)下,使用PBE泛函[12]來描述交換關聯能量,并通過Hubbard參數修正(即DFT+U)[13],其中,Ni,Fe,Ti的d軌道U值分別設置為6.3[14],5.3 eV和4.0[15]eV。截斷能設置為520 eV,自洽計算(SCF)迭代的收斂閾值為10-5 eV,而幾何結構的收斂以力為標準,設置為-0.2 eV/nm。使用Monkhorst-Pack方案對布里淵區進行k點采樣,對于CaTiO3和NiZnFe4O8,使用4×4×4取樣方案,對于MgTiO3,使用9×9×2的k點取樣。使用CI-NEB[16]方法,通過插入5個擴散中間態構型來研究摻雜陽離子在不同摻雜體系中的擴散行為。

2 結果與討論

2.1 陽離子取代結構模型

使用VASP對建立的模型進行結構優化,計算各模型對應的摻雜形成能并列于表3,摻雜能的大小表示摻雜的難易程度。從表3可以看出,對于MgTiO3體系,Ni[Ti],Zn[Ti]和Fe[Ti]的形成能都低于Ni[Mg],Zn[Mg]和Fe[Mg]的形成能,表明3種陽離子都更傾向于占據體系中Ti所處位置,CaTiO3體系具有相同的情況;而在NiZnFe4O8體系中,3種取代原子最低形成能的取代為Mg[Zn],Ca[Zn]以及Ti[Fe],進入NiZn鐵氧體的Mg和Ca優先占據Zn所處的位點,而Ti則優先占據Fe所在的位點。值得注意的是,摻雜形成能的值有正有負,負的摻雜能表示摻雜可以自發進行,而正的摻雜能則表示摻雜需要額外的能量。

在研究了各類陽離子分別進入3種體系后的優先占據位點的基礎上,將陽離子摻雜在其最可能的位點建立模型并進行充分弛豫,并將完全弛豫后的晶格常數列于表4。從表4可以看到各摻雜體系的晶格常數相比于未摻雜前雖有變化,但變化較小。摻雜元素的引入會使晶體內部平衡狀態受到破壞,使摻雜元素周圍的原子發生不同程度的畸變[17],以摻雜元素與其近鄰氧鍵長的變化表征其畸變大小,表5列出了摻雜元素與最近鄰O的鍵長。從表5可以看出,對于MgTiO3體系,陽離子引入后的鍵長都有所增加,其中ZnO鍵鍵長變化最大,拉伸達到5.1%; 而對于CaTiO3體系, ZnO鍵鍵長也是處于拉伸狀態,卻只有1.5%,其他兩種陽離子在CaTiO3體系中鍵長處于縮短狀態。在NiZnFe4O8體系中,涉及到Zn和Fe的替換,其中Ti替換Fe后鍵長縮短,而Mg,Ca替換Zn后鍵長拉伸,尤以Ca取代后畸變最為明顯,畸變達到9.1%,是所有摻雜原子中畸變度最大的。

為進一步研究陽離子摻雜對體系電子結構的影響,對摻雜前后的結構做了態密度分析,研究各體系的總態密度TDOS以及各原子部分軌道的分波態密度PDOS。圖2-圖4分別為MgTiO3,CaTiO3以及NiZnFe4O8 3種體系摻雜前后的總態密度(TDOS)和分波態密度(PDOS),其中虛線表示費米能級的位置。對于MgTiO3體系,通過觀察其TDOS發現Fe,Ni和Zn的摻入導致體系的帶隙降低,從其PDOS圖中可以看到帶隙降低的主要原因是摻雜的過渡金屬的3d軌道在費米能級附近有一定的貢獻,且與O的2p軌道還有一定的雜化,表明其與O成鍵。對于同樣摻雜了過渡金屬的CaTiO3體系,其態密度呈現出和MgTiO3相同的狀態,過渡金屬的3d軌道在

費米能級附近有輕微貢獻,導致摻雜后體系態密度的價帶有輕微右移,這也導致了禁帶寬度的降低。NiZnFe4O8體系(圖4)相比于MgTiO3和CaTiO3體系,摻雜后態密度變化要復雜得多,其涉及到兩種類型的摻雜:一種是Mg和Ca對Zn的替換摻雜,對于這種摻雜,其費米能級從原來導帶底的位置移動到了價帶頂,通過PDOS圖可以看出主要是由于O的2p軌道發生移動,Ca的2p軌道的貢獻在遠離費米能級位置,而Mg的2p軌道在價帶頂附近的貢獻很微弱;另一種則是Ti對Fe的替換摻雜,Ti的引入并沒有改變其費米能級位置。3種摻雜離子的引入都降低了Ni的3d軌道在費米能級處的貢獻,表現為提升了摻雜后體系的禁帶寬度。

2.2 陽離子的擴散機制

為了模擬摻雜原子在晶體內的擴散,選定最穩定的摻雜結構構建擴散路徑。對于擴散路徑的構建,以摻雜離子所在位置為初態,以最近鄰同種離子所處位置為末態,但對于NiZnFe4O8的尖晶石結構,陽離子所處的位置只有8a和16b兩種,Ni和Zn同處于體系中的8a位置,所以當Mg和Ca替換取代Zn時,其擴散路徑設定為從Zn所處位置擴散到最近鄰Ni所處的位置。各擴散路徑如圖5所示,圖5(c)的①為Mg和Ca的擴散路徑,②為Ti的擴散路徑。圖5給出了3種體系中摻雜陽離子的擴散勢壘圖。圖5(a)和圖5(b)是Fe,Ni和Zn離子分別在MgTiO3和CaTiO3中的擴散勢壘圖,可以看出,盡管兩種體系摻雜的是同種離子,但他們的擴散勢壘卻大不相同,陽離子在CaTiO3中的擴散勢壘要遠大于在 MgTiO3 中的擴散勢壘,表明3種陽離子在CaTiO3中的擴散行為難以發生,其主要原因可能是由于摻雜陽離子在擴散過程中的擴散通道較小。從CaTiO3的晶體結構可以看出,兩個Ti之間的擴散通道是由2個Ca和4個O構成的,摻雜陽離子的擴散遷移需要克服離子之間的相互作用,故比較困難,而在 MgTiO3 中,相應的離子擴散勢壘也大大降低,其中Fe的擴散勢壘從4.024 eV降到了1.446 eV,Zn則從5.507 eV直接降到了0.918 eV,而Ni的擴散勢壘變化沒有其他兩種離子那么大,只降低了0.716 eV。Mg,Ca和Ti 3種離子在NiZnFe4O8中的擴散行為較為容易發生,勢壘最低的Ca僅需0.329 eV就可以驅動擴散的進行,遠低于其他陽離子的擴散勢壘,Ti在NiZnFe4O8中最不容易發生擴散,其擴散勢壘為1.181 eV,但相較于在MgTiO3和CaTiO3體系中的擴散行為,也很容易發生。

通過對陽離子在各體系中擴散勢壘的比較可以得出,對于NiZn鐵氧體和介電陶瓷共燒體系,陽離子的擴散是存在的,其可能的形式是:Fe,Ni和Zn在進入介電陶瓷后,首先會與介電陶瓷發生摻雜反應,進入陶瓷晶體內部,隨之發生擴散,但由于這3種陽離子在CaTiO3中的擴散勢壘較大,因此其最可能的擴散通道為MgTiO3晶體內部的擴散;同理,Mg,Ca和Ti在進入NiZn鐵氧體晶體內部后會發生擴散遷移。

擴散系數的阿倫尼烏斯表達形式[18]為:

D=D0exp(-Ea/kBT)(4)

式中:D0為擴散常數;Ea為擴散勢壘;kB和T分別為玻爾茲曼常數和絕對溫度,kB=1.38×10-23 J/K。

根據過渡態理論[19],原子在固體中的躍遷率表達式為:

ν=υexp(-Ea/kBT)(5)

式中ν為摻雜原子的振動頻率。由原子躍遷率可以將擴散系數近似表達為[20-21]:

D=l2υexp(-Ea/kBT)(6)

式中l是原子擴散位移,可通過摻雜原子擴散始末位置得出。根據溫特-齊納理論[22-23],振動頻率ν可以近似地表示為:

ν=(2Ea/ml2)1/2(7)

式中m為單一原子的質量。根據式(7)可得出各摻雜原子的振動頻率。

根據式(4)- 式(7)分別計算了MgTiO3,CaTiO3 以及NiZnFe4O8體系摻雜原子的振動頻率以及在1 500 K下的擴散系數,并連同各體系的擴散勢壘、位移以及原子質量的相關數據列于表6。

從表6可以看出不同體系以及同一體系的擴散系數都有1~10個數量級不等的差距,其主要原因在于擴散勢壘的差距,如擴散勢壘最大的CaTiO3中的擴散與勢壘最小的NiZnFe4O8中的擴散,兩者擴散系數差別在10個數量級甚至更多。而其他如擴散位移和原子質量的差距并非主要因素,如Ti和Ca在NiZnFe4O8中的擴散,盡管Ti的擴散距離小于Ca,但由于Ca的擴散勢壘更小,導致其擴散系數比Ti的擴散系數高3個量級。

2.3 擴散距離評估

圖6為介電陶瓷(標記為A物質)和NiZn鐵氧體(標記為B物質)共燒界面擴散示意圖。A物質和B物質的異質界面處在燒結過程中會發生離子的相互擴散,并最后隨燒結過程的完成達到穩定的擴散濃度。在燒結保溫過程中,A物質的Mg2+,Ca2+,Ti4+ 不斷地向B物質擴散,B物質中的Ni2+,Zn2+,Fe3+ 不斷地向A物質擴散。A,B兩種物質形成一對互擴散偶。

在燒結過程中,A,B物質的互擴散偶模型中離子濃度分布會發生變化,即發生不穩定擴散,其擴散通量會隨位置與時間而變化。假設離子的擴散系數D不隨擴散濃度而變化,此時離子的擴散濃度隨時間t和位置x的變化關系滿足Fick第二定律:

基于上述理論,對A/B異質擴散偶模型作出以下假設:

(1)將A和B物質看作截面一致,且成分均勻的擴散偶;

(2)擴散系數D相對于擴散距離x不發生變化;

(3)A和B兩接觸端的質量濃度仍然保持他們原來的數值,不受擴散的影響;

(4)兩側擴散偶初始值厚度相比于擴散層厚度可看作無限長;

(5)A物質向B物質擴散,A物質的起始濃度為1,B物質的原始濃度設為0;

(6)擴散接觸面設置為擴散位置x=0 μm。

因此,擴散模型可以看作是一維無限長或半無限長的擴散偶,經過數學處理可以得到半無限長的擴散方程,如式(9)所示:

基于式(9),假設整個擴散過程為3 h,則在1 500 K下,陽離子的擴散濃度與距離界面處x位置的關系如圖7所示。

從圖7可以看出,對于MgTiO3,界面處可以觀察到Zn2+ 的明顯擴散,而CaTiO3中離子的擴散系數都很小,因此在界面處無法觀察到明顯擴散,而對于NiZnFe4O8,3種離子都有一定距離的擴散深度,Ca2+ ,Mg2+擴散進入NiZn鐵氧體的距離為400~1 000 μm。過長的過渡擴散區會對復合基板性能產生較大影響。

3 結論

本文采用基于密度泛函理論的第一性原理方法模擬了介電陶瓷/NiZn鐵氧體共燒過程中離子的互擴散行為,計算了陽離子在MgTiO3,CaTiO3以及NiZnFe4O8體系中的互擴散勢壘,發現陽離子在CaTiO3中的擴散勢壘較大,擴散難以進行,而在MgTiO3 和NiZnFe4O8中的擴散勢壘相對較小,擴散容易發生。基于擴散勢壘估算了陽離子在3種體系中的擴散系數,使用半無限長擴散偶模型來評估陽離子的濃度分布,結果表明Mg,Ca離子在NiZn鐵氧體中有明顯的擴散現象,Zn離子在MgTiO3中也有一定的擴散深度。研究結果可為共燒復合材料的設計及共燒過程中擴散的控制提供理論依據。

參考文獻

[1] 李俊, 冷觀武, 彭虎. 旋磁鐵氧體材料的微波燒結及在環行器中的應用研究[J]. 磁性材料及器件, 2005, 36(5): 47-49.

[2] 盧憲俊. 微波鐵氧體材料鐵磁共振線寬和復介電常數測試技術研究[D]. 成都: 電子科技大學, 2013.

[3] WAKINO K. Recent development of dielectric resonator materials and filters in Japan[J]. Ferroelectrics, 1989, 91(1): 69-86.

[4] IKRAM S, JACOB J, ARSHAD M I, et al. Tailoring the structural, magnetic and dielectric properties of Ni-Zn-CdFe2O4 spinel ferrites by the substitution of lanthanum ions[J]. Ceramics International, 2019, 45(3): 3563-3569.

[5] ZHANG M M, LI L X, XIA W S, et al. Structure and properties analysis for MgTiO3 and (Mg0.97M0.03)TiO3 (M=Ni, Zn, Co and Mn) microwave dielectric materials[J]. Journal of Alloys and Compounds, 2012, 537: 76-79.

[6] YANG X F, FU J X, JIN C J, et al. Formation mechanism of CaTiO3 hollow crystals with different microstructures[J]. Journal of the American Chemical Society, 2010, 132(40): 14279-14287.

[7] LI X, LI Q, XIA Z G, et al. Effects on direct synthesis of large scale mono-disperse Ni0.5Zn0.5Fe2O4 nanosized particles[J]. Journal of Alloys and Compounds, 2008, 458(1/2): 558-563.

[8] WANG Q J, WANG J B, ZHONG X L, et al. Magnetism mechanism in ZnO and ZnO doped with nonmagnetic elements X (X=Li, Mg, and Al): A first-principles study[J]. Applied Physics Letters, 2012, 100(13):132407.

[9] CHEN H, LI X C, WAN R D, et al. A DFT study on modification mechanism of (N, S) interstitial co-doped rutile TiO2[J]. Chemical Physics Letters, 2018, 695: 8-18.

[10]KRESSE G, FURTHMüLLER J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set[J]. Physical Review B, Condensed Matter, 1996, 54(16): 11169-11186.

[11]KRESSE G, FURTHMüLLER J. Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set[J]. Computational Materials Science, 1996, 6(1): 15-50.

[12]PERDEW J P, BURKE K, ERNZERHOF M. Generalized gradient approximation made simple[J]. Physical Review Letters, 1996, 77(18): 3865-3868.

[13]DUDAREV S L, BOTTON G A, SAVRASOV S Y, et al. Electron-energy-loss spectra and the structural stability of nickel oxide: An LSDA+U study[J]. Physical Review B, 1998, 57(3): 1505-1509.

[14]WANG L, MAXISCH T, CEDER G. Oxidation energies of transition metal oxides within the GGA+U framework[J]. Physical Review B, 2006, 73(19): 195107.

[15]SOLOVYEV I V, DEDERICHS P H, ANISIMOV V I. Corrected atomic limit in the local-density approximation and the electronic structure of d impurities in Rb[J]. Physical Review B, 1994, 50(23): 16861-16871.

[16]HENKELMAN G, UBERUAGA B P, JNSSON H. A climbing image nudged elastic band method for finding saddle points and minimum energy paths[J]. The Journal of Chemical Physics, 2000, 113(22): 9901-9904.

[17]李旭, 王建川, 杜勇. Ni,Co,Al摻雜尖晶石LiMn2O4的第一性原理計算[J]. 粉末冶金材料科學與工程, 2021, 26(5): 387-395.

[18]ARRHENIUS S. On the reaction velocity of the inversion of cane sugar by acids [M]∥Selected readings in chemical kinetics. Pergamon, 1967: 31-35.

[19]VINEYARD G H. Frequency factors and isotope effects in solid state rate processes[J]. Journal of Physics and Chemistry of Solids, 1957, 3(1/2): 121-127.

[20]GOMER R. Diffusion of adsorbates on metal surfaces[J]. Reports on Progress in Physics, 1990, 53(7): 917.

[21]KUTNER R. Chemical diffusion in the lattice gas of non-interacting particles[J]. Physics Letters A, 1981, 81(4): 239-240.

[22]WERT C, ZENER C. Interstitial atomic diffusion coefficients[J]. Physical Review, 1949, 76(8): 1169-1175.

[23]WERT C A. Diffusion coefficient of C in α-iron[J]. Physical Review, 1950, 79(4): 601-605.

收稿日期:2023-02-27;修回日期:2023-05-16

作者簡介:第一作者,張凱(1997— ),男,碩士研究生; 通信作者,畢鵬(1985— ),博士,講師,研究方向為計算材料學,E-mail:bipeng010@swust.edu.cn

主站蜘蛛池模板: 一本大道香蕉中文日本不卡高清二区| 国产99视频在线| 欧美成人午夜视频| 69av免费视频| 免费福利视频网站| 在线亚洲精品自拍| 亚洲无码高清一区二区| 亚洲精品波多野结衣| 国产黄网永久免费| 成人看片欧美一区二区| 福利一区三区| 国产真实乱子伦精品视手机观看| 欧美成人精品一级在线观看| 国产精品香蕉在线观看不卡| 青青操视频在线| 国产区在线观看视频| 国产精品手机视频| 亚洲成人手机在线| 亚洲成av人无码综合在线观看| 天天色综网| 国产乱人伦AV在线A| 被公侵犯人妻少妇一区二区三区| 91 九色视频丝袜| 国产一二视频| 一区二区三区成人| 国产精品99r8在线观看| 国产幂在线无码精品| 9啪在线视频| 亚洲 成人国产| 亚洲午夜国产精品无卡| 亚洲无线国产观看| 亚洲人免费视频| 亚洲第一精品福利| 亚洲三级电影在线播放| 福利在线不卡一区| 国产成年无码AⅤ片在线| 国产无码高清视频不卡| 亚洲最黄视频| 国产网站免费| 精品欧美一区二区三区在线| swag国产精品| 国产三级视频网站| 久久久久久久久久国产精品| 热伊人99re久久精品最新地| 无码啪啪精品天堂浪潮av | 五月天综合网亚洲综合天堂网| 亚洲欧美色中文字幕| 欧美日韩国产在线人成app| 青青青视频免费一区二区| 日韩在线播放欧美字幕| 日本高清有码人妻| 久久永久免费人妻精品| 国产精品九九视频| 波多野结衣一区二区三区AV| 日韩欧美国产区| 欧美精品伊人久久| 国产三级成人| 经典三级久久| 中文字幕欧美成人免费| 在线色综合| 中文字幕欧美成人免费| 国产丝袜91| 国产乱子伦视频三区| 国内精品一区二区在线观看| 国产噜噜在线视频观看| 香蕉视频国产精品人| 亚洲色图在线观看| 国产亚卅精品无码| 青草国产在线视频| 在线观看的黄网| av在线手机播放| 亚洲精品天堂自在久久77| 久热这里只有精品6| 女人毛片a级大学毛片免费| 潮喷在线无码白浆| www.精品国产| 天堂成人av| 国产乱人伦偷精品视频AAA| 国产亚洲精| 国产精品开放后亚洲| 日韩二区三区| 九色91在线视频|