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

基于密度泛函理論的KF-NaF-AlF3低溫電解質體系勢參數擬合及分子動力學模擬

2021-01-19 02:35:14王景坤國輝張紅亮李劼
中南大學學報(自然科學版) 2020年12期
關鍵詞:體系結構

王景坤,國輝,張紅亮,3,李劼,3

(1.中南大學冶金與環境學院,湖南長沙,410083;2.中國宏橋集團有限公司,山東鄒平,256200;3.中南大學難冶有色金屬資源高效利用國家工程實驗室,湖南長沙,410083)

電解鋁工藝是一種能耗較高的工藝,生產1 t Al 電解鋁消耗電12 000~14 000 kW·h(電力成本占電解鋁總成本的35%~40%)。100 多年來,電解鋁所用的電解質幾乎沒有變化,一直以NaF-AlF3體系為主[1]。NaF-AlF3體系的初晶溫度較高,一般為910~960 ℃(NaF 的物質的量與AlF3的物質的量之比為2.0~2.7),這極大地增加了電解鋁的生產成本[2]。降低電解質體系的初晶溫度可以有效降低能耗,因此,尋找一種具有較低初晶溫度和穩定物理化學性質的新型電解質是一種較好的解決辦法。KF-AlF3基電解質體系具有比NaF-AlF3體系更低的初晶溫度(低于800 ℃),且對氧化鋁的溶解能力較強,得到了研究者們的廣泛關注[3-5]。但由于電解過程中鈉的積累(電解過程中原料氧化鋁帶來的鈉),KF-AlF3電解質體系會逐漸轉變為KF-NaFAlF3電解質體系,因此,KF-NaF-AlF3熔鹽體系是較適宜的低溫鋁電解質體系。目前,人們對KFNaF-AlF3體系的實驗研究較少,現有研究主要集中在一些基本的物理化學性質如熔體的密度和液相線溫度等[1-3],對于KF-NaF-AlF3體系的離子結構研究較少。開展鋁電解熔鹽實驗有以下難點:1)鋁電解質體系的高溫強腐蝕性對實驗設備的損耗較大;2)熔融氟化物易揮發,難以準確測定其結構和性質;3)熔融氟化物的結構較復雜,實驗研究困難。隨著拉曼光譜[6]、核磁共振[7]以及計算機模擬技術的發展[8-12],鋁電解質高溫熔體的研究取得了許多成果。相比于拉曼光譜以及核磁共振技術,計算機模擬技術尤其是第一性原理和分子動力學模擬技術具有很大的成本優勢[13-18],但直接從第一性原理出發,對熔融鋁電解質體系進行完全量子化學處理,不僅在計算方法上存在一定困難,而且受時間限制,難以全面地進行計算機模擬。近年來,研究者們嘗試結合第一性原理和經典的分子動力學模擬提出第一性原理分子動力學模擬(FPMD)方法,在KCl-LiCl[13],Li2BeF4[14]和LiFNaF-KF[15]等體系得到了很好的應用。FPMD 方法的優勢是無需選擇勢函數及勢參數,而是直接基于第一性原理計算。但這種方法能計算的原子個數較少,一般只有100~200個原子,這與真實體系有較大差距,并且模擬時間較短,僅為10 ps 左右[13]。盡管FPMD方法已經日趨成熟,但對于KFNaF-AlF3熔鹽體系的FPMD 研究還很少。Lü 等[19]采用FPMD方法對KF-NaF-AlF3體系的離子結構和輸運性質進行了模擬,但由于體系的離子數較少,模擬時間較短,無法完全反映KF-NaF-AlF3體系真實的離子結構形式,并且Lü 等[19]的研究體系與常規的低溫電解質體系的組成不同[3],不同組成的電解質體系的性質存在一定差別。在一定的物理模型基礎上,采用分子動力學模擬方法研究鋁電解質熔體的結構和性質是一種較好的手段[20]。開展KF-NaF-AlF3體系的分子動力學模擬最關鍵的就是勢函數以及相關勢參數的選擇,Buckingham 勢作為一種二體勢已經廣泛應用于NaF-AlF3體系模擬[17],但KF-NaF-AlF3體系相應的勢參數較少,因而,采用分子動力學模擬方法存在一定困難,為此,本文作者首先基于DFT理論計算構建KF-NaFAlF3體系的Buckingham 勢參數,然后,基于構建的勢參數采用經典分子動力學模擬研究KF-NaFAlF3體系的離子結構及其變化規律,以便為低溫電解質的實驗研究提供理論支撐。

1 理論方法

1.1 勢函數模型

KF-NaF-AlF3體系中每個離子對在經典的二體勢中主要包含3種相互作用勢[21]:體系中不同離子對間的庫侖力引發的靜電相互作用勢Vcharge、體系中離子的電子云重疊引發的斥力相互作用勢Vrepulsion和色散力引發的色散相互作用勢Vdispersion。總相互作用勢Vtotale=Vcharge+Vrepulsion+Vdispersion。色散力引發的色散相互作用勢Vdispersion是體系中離子的正負電荷瞬間不重合所致,這里只考慮色散力中的偶極-偶極相互作用。本文采用Buckingham 勢函數描述體系的斥力相互作用和偶極-偶極相互作用[17]。體系的Buckingham 勢以及靜電相互作用勢的函數形式如下:

式中:rij為2個粒子之間距離;Aij,βij和Cij為用于描述兩粒子之間的相互作用的參數;qi為i 粒子所帶電荷;ε0為真空介電常數。

1.2 勢參數的構建

在硅酸鹽熔體的分子動力學模擬過程中,研究者針對硅酸鹽的基本結構單元即硅氧四面體,對其第一性原理的單點能進行計算,通過其勢能面隨原子間距離的變化導出原子間的相互作用勢,進而用該勢能擬合相應的勢參數[20]。考慮到鋁電解熔鹽體系同樣屬于高溫體系,并且熔體結構復雜多變,這里借鑒硅酸鹽熔體的勢參數擬合方法,對KF-NaF-AlF3體系中不同離子對第一性原理單點能進行計算,得到Al-Al,Al-F,F-F,K-Al,K-F,K-K,K-Na,Na-Al,Na-F 和Na-Na 離子對的勢能曲線,對勢能進行分解得到Buckingham 勢能并進行擬合。

1.2.1 不同離子對單點能的計算

采用廣義梯度近似(GGA)中的Perdew-Burke-Ernzerhof(PBE)相互交聯函數[18]計算第一性原理離子對的單點能。使用Ultrasoft超軟贗勢來處理離子實和價電子的相互作用,其中,電子Na 2s22p63s1,Al 3s23p1和K 3s23p64s1和F 2s22p5被視為價電子。DFT-D2 方法[12]用于處理粒子間的范德華弱相互作用。對計算所用的截斷能、K點網格和模型盒子的收斂性進行多次測試,采用截斷能為850 eV及10×10×1 的K 點網格分布和長×寬×高為(15×10-10)m×(15×10-10)m×(15×10-10)m 的盒子,所有計算均在CASTEP[22]模塊中完成。根據計算得到的不同離子對的單點能,導出勢能隨離子間距離變化的相互作用勢,用于后續的勢能分解。需說明的是,在基于密度泛函理論的第一性原理計算中,離子所帶電荷隨離子間距離變化而變化,而分子模擬中離子所帶電荷是固定不變的,因此,有必要對KFNaF-AlF3體系中不同離子所帶電荷進行研究。

1.2.2 KF-NaF-AlF3體系中離子電荷的確定

根據式(2)可知,在對靜電相互作用勢能進行分解時,首先要確定離子所帶電荷。進行第一性原理計算時,離子所帶電荷隨距離變化而不斷變化,而經典的分子模擬過程認為電荷是固定不變的。為了解決這一問題,本文擬采用第一性原理模擬真實的KF-NaF-AlF3熔鹽體系,通過分析KFNaF-AlF3體系中不同離子的平均Mulliken電荷確定最終的分子模擬中離子所帶電荷。用于FPMD 計算的KF-NaF-AlF3體系的化學組成為11 個K 離子、9 個Na 離子、15 個Al 離子和65 個F 離子,共100個離子。采用固定粒子數、體積和溫度的NVT 系綜進行FPMD 模擬,溫度設置為827 ℃,略高于KF-NaF-AlF3熔鹽體系的初晶溫度[23],模擬盒子的密度[3]設置為1.846 g/cm3。熔鹽體系先在NVT系綜下進行6 000步模擬,得到接近真實體系的熔鹽離子結構,然后進行4 000步結構弛豫,模擬時間步長設置為1 fs,總的模擬時間為10 ps,其他計算條件與1.2.1 節中的相同。采用FPMD 方法模擬得到最終的穩定構型用于分析KF-NaF-AlF3體系中離子所帶電荷。

在KF-NaF-AlF3體系中,更關注的是Al-F離子間的相互作用,因為其他離子對的相互作用較弱,難以形成穩定的化合物[24]。依據這一原則,本文著重考察KF-NaF-AlF3體系中氟離子所帶電荷,而其他離子所帶電荷根據電中性原則得到。通過對FPMD 計算得到的KF-NaF-AlF3體系的穩定構型進行Mulliken 電荷分析,得到體系中65 個F 離子所帶電荷,如圖1 所示。總體來看,氟離子電荷在-0.71 e到-0.63 e之間波動(其中,e為電子),平均值約為-0.67 e,因此,本文將KF-NaF-AlF3體系中F 離子的電荷設為-0.67 e,依據電中性原則,Na,K和Al離子的電荷依次設為+0.67 e,+0.67 e和+2.01e。

圖1 KF-NaF-AlF3體系中氟離子的Mulliken電荷Fig.1 Mulliken charge of F ions in KF-NaF-AlF3 system

1.2.3 靜電相互作用勢能分解

KF-NaF-AlF3體系中不同離子對的Buckingham勢能可以根據總勢能分解得到,根據1.1 節的分析,Buckingham勢能的計算式如下:

根據1.2.2節的KF-NaF-AlF3體系中離子電荷的計算結果,可以得到不同離子對的靜電相互作用勢,然后對總勢能進行分解得到不同離子對的Buckingham 勢能。Al-F 配離子的結構對整個KFNaF-AlF3體系的性質具有重要影響[24]。圖2 所示為Al-F離子對的勢能曲線。從圖2可見:Al-F離子對的相互作用勢隨距離的變化在1.5×10-10~2.0×10-10m之間存在明顯的最低點,這正是KF-NaF-AlF3體系中主要的陽離子與陰離子強相互作用的體現,因為Al-F 離子間的極性相差較大,結合能力強,在熔鹽中易形成穩定的配位結構,并且這一能量最低點所對應的距離與熔鹽中Al-F 離子的平均鍵長相對應[18-19],這說明本文關于離子對勢能面的計算過程是合理的。熔鹽中其他離子對的勢能曲線如圖3所示。從圖3可見:除Na-F離子對外,并未看到明顯的勢能曲線的最低點。這說明其他離子對的極性差異較小,結合能力弱。

1.2.4 勢參數擬合

基于最小二乘原理對1.2.3 節得到的Buckingham 勢能曲線進行非線性曲線擬合,所有擬合過程均在Matlab平臺中完成,原理如下:

圖2 Al-F離子對的勢能曲線Fig.2 Potential energy curves of Al-F ion pairs

圖3 KF-NaF-AlF3體系中不同離子對的勢能曲線Fig.3 Potential energy curves of different ion pairs in KF-NaF-AlF3 system

式(4)為最小二乘法進行非線性擬合的基本原理,其中,F(rij)為擬合的勢函數;VBuckingham(rij)為計算出的勢函數;式(5)為Buckingham勢函數,其中參數Aij,βij和Cij是本文需要擬合的勢參數,基于1.1節的分析可知,本文采用Buckingham勢函數來處理KF-NaF-AlF3體系中不同離子對之間的非靜電相互作用勢;式(6)表示體系中的非靜電相互作用勢。擬合得到最終的Buckingham 勢參數如表1所示。

表1 KF-NaF-AlF3體系中不同離子對的Buckingham勢參數Table 1 Buckingham potential parameters of different ion pairs in KF-NaF-AlF3 system

1.3 分子動力學模擬

本文研究的KF-NaF-AlF3熔鹽體系化學組成和模擬盒子的體積如表2 所示,初始結構模型由Packmol 軟件[25]隨機地將粒子投放到指定大小的盒子中得到。在經典分子動力學模擬中,使用Verlet蛙跳算法求解牛頓運動方程,時間步長設置為1 fs。EWALD求和算法[18]作為處理長程相互作用應用最普遍的方法,在這里被用于處理體系中粒子間的庫侖相互作用和偶極相互作用,同時,Buffer寬度設置為0.5×10-10m,能量計算精度為4.184×10-5kJ/mol[19]。短程相互作用的截斷半徑為15×10-10m,離子的電荷設置為+0.67 e(Na),+2.01 e(Al),-0.67 e(F)和+0.67 e(K),這來源于1.2.2 節的DFT計算結果。周期性邊界條件被用于模擬盒子,用于描述1個沒有邊界的無限液態系統,使得計算結果更加可靠。模擬盒子系統首先在NPT[24]系綜下以0.1 MPa壓力升溫到827 ℃,模擬時長為2 000 ps,這意味著模擬過程中保持系統的粒子數(N)、壓力(P)和溫度(T)為常數。之后,使用NVT系綜繼續進行3 000 ps結構馳豫。最后,使用Matlab軟件統計分析結構馳豫3 000 ps 的模擬盒子中的粒子軌跡數據,計算KF-NaF-AlF3熔鹽的離子結構。所有的分子動力學計算均在Forcite模塊中完成[26]。

表2 分子動力學模擬條件Table 2 Molecular dynamics simulation conditions

2 結果與討論

2.1 KF-NaF-AlF3熔鹽體系的局部離子結構

在溫度為827 ℃、壓力為0.1 MPa下模擬得到KF-NaF-AlF3熔鹽體系的穩定構型如圖4 所示。從圖4 可以看出:K 和Na 離子隨機分布在模擬盒子中,并且離子之間的距離較大,這是因為K-K,K-Na 和Na-Na 各離子極性差異較小,且均無共價特性。模擬盒子中局域的離子結構是由四配位[AlF4]-、五配位[AlF5]2-和六配位[AlF6]3-絡合離子集團主導,分別對應扭曲的四面體、三角雙錐體和八面體構型。同時,由于橋氟離子存在,熔鹽中形成了Al-F-Al 復雜離子構型。KF-NaF-AlF3熔鹽體系雖然失去了長程有序狀態,但局域離子結構仍然保持著短程有序,即熔鹽中存在大量的四配位[AlF4]-、五配位[AlF5]2-和六配位[AlF6]3-絡合離子集團。

圖4 KF-NaF-AlF3熔鹽體系的微觀離子結構(KF質量分數為4%)Fig.4 Micro-structure of KF-NaF-AlF3 molten salt system

2.2 徑向分布函數

徑向分布函數(RDF)用于描述液體結構,可采用從分子動力學模擬軌跡提取出的RDF 分析氟化物熔鹽的局域離子結構[15]。徑向分布函數反映了以r處的粒子為中心,半徑Δr范圍內出現另一個粒子的概率。

式中:V為分子動力學模擬盒子的體積;N為粒子數目;nij(r,Δr)為粒子j在指定的Δr處截斷范圍內圍繞中心粒子i的平均數目。

KF-NaF-AlF3體系中離子對的RDF 如圖5 所示。從圖5 可見:當r>8×10-10m 時,g(r)逐漸趨近于1,這與熔體近程有序、遠程無序的結構形式相吻合;Al-F離子對的RDF在1.75×10-10m左右存在1個高而尖的峰值,這意味著Al-F離子間的相互作用較強,在熔鹽中易形成較復雜的絡合離子,并且這一峰值與Al-F 的平均鍵長相對應[18];Al-Al 離子對的RDF 可以反映熔鹽中離子的聚合程度,在3.69×10-10m 和6.33×10-10m 處分別出現了Al-Al 離子對的第一峰值和第二峰值,這與文獻[21]報道的結果相一致。Al-Al離子對的第一峰值半徑約為Al-F離子對第一峰值半徑的2倍,這說明熔鹽中形成了Al-F-Al 構型的復雜絡合離子,即氟離子作為橋離子連接2個鋁離子。一般徑向分布函數的第一峰值代表著某一離子對的平均鍵長,計算結果如表3所示。

圖5 KF-NaF-AlF3體系中離子對的徑向分布函數(RDF)(KF質量分數為4%)Fig.5 RDF of ion pairs in KF-NaF-AlF3 system(4%KF)

2.3 鍵角分布

KF-NaF-AlF3熔鹽中F-Al-F 鍵角分布如圖6 所示。理想的[AlF6]3-八面體擁有8 個90°和3 個180°的F-Al-F鍵角;理想的[AlF5]2-三角雙錐體擁有6個90°,3個120°以及1個180°的F-Al-F鍵角;而理想的[AlF4]-正四面體則擁有6個109.5°的F-Al-F鍵角。從圖6可以看出:F-Al-F的鍵角分布曲線的第一峰值和第二峰值分別位于91°和116°,分別對應著1個八面體或三角雙錐體的鋁氟絡合離子集團結構,同時,位于173°左右的第三峰是三角雙錐體的特征。本文在108°左右發現了第四峰,對應著正四面體離子集團結構,但其峰值較小,說明KF-NaF-AlF3熔鹽中存在著少量的四配位[AlF4]-絡合離子集團結構。同時,這些峰的位置都略微偏離其在理想構型中標準鍵角的位置,表明這些絡合離子集團結構在熔融條件下是扭曲的。通過對上述鍵角分布情況進行分析,發現KF-NaF-AlF3熔鹽中絡合離子集團主要以六配位[AlF6]3-和五配位[AlF5]2-為主,同時存在少量的四配位[AlF4]-絡合離子。

表3 KF-NaF-AlF3體系中各離子對的平均鍵長Table 3 Average bond lengths of ion pairs in KF-NaF-AlF3 system10-10 m

圖6 KF-NaF-AlF3熔鹽體系中F-Al-F鍵角分布(KF質量分數為4%)Fig.6 F-Al-F bond angle distribution in KF-NaF-AlF3 system(4%KF)

2.4 F離子類型分布

事實上,KF-NaF-AlF3熔鹽中Al-F配離子的類型并不是簡單的[AlF6]3-,AlF5]2-和[AlF4]-,因為熔鹽中的氟離子可能會作為橋離子連接2 個Al 離子形成Al-F-Al 構型,這大大增加了絡合離子的復雜程度。為了解釋這一現象,本文統計分析了KFNaF-AlF3熔鹽體系中氟離子類型分布規律,結果如圖7所示。F離子類型為:橋氟Fb,氟離子以Al-FAl 形式連接2個鋁離子;終端氟Ft,氟離子與1個近鄰鋁離子連接;自由氟Ff,氟離子沒有與近鄰鋁離子相互作用。F 離子類型決定了KF-NaF-AlF3熔鹽離子結構的聚合程度,并對熔鹽的傳輸特性產生至關重要的影響。從圖7 可以看出:隨著KF 濃度增大,熔鹽中橋氟離子的濃度逐漸增大,終端氟離子的濃度相應降低,而自由氟離子的濃度維持在較低值。這說明隨著KF 濃度增大,KF-NaFAlF3熔鹽體系的離子結構逐漸復雜,體系中離子的擴散性能也會相應降低。

2.5 復雜熔鹽離子結構的形成

圖7 KF-NaF-AlF3熔鹽體系中F離子類型分布Fig.7 F ion type distribution in KF-NaF-AlF3 system

通過前面分析可知,KF-NaF-AlF3熔鹽體系中存在常規的[AlF6]3-,[AlF5]2-和[AlF4]-,這些簡單的絡合離子構型可能會通過橋氟離子形成體積更大更復雜的陰離子構型。為了更直觀地表示KF-NaFAlF3熔鹽體系中復雜離子結構的形成過程,本文統計了熔鹽中主要陰離子Al-F 絡合離子的結構,通過定義相應的截斷距離(Al-F 徑向分布函數的第一谷值半徑)[16],可以計算Al離子周圍F離子的數目,如圖8所示。從圖8可知:熔鹽中陰離子構型主要是F-,[AlF4]-,[AlF5]2-和[AlF6]3-;同時,由于橋氟離子的存在,熔鹽中形成了更復雜的絡合離子,如[Al2F8]2-, [Al3F10]-, [Al4F14]2-, [Al4F15]3-和[Al4F16]4-等。

圖8 KF-NaF-AlF3熔鹽體系中Al-F絡合離子的結構(KF質量分數為4%)Fig.8 Structure of Al-F complex ions in KF-NaF-AlF3 system(4%KF)

通過統計分析3 000 ps的粒子軌跡數據,得到KF-NaF-AlF3熔鹽體系中不同Al-F絡合離子的百分比,所有的統計結果均基于1.3節中的分子動力學模擬過程得到,如圖9 所示,其中Al-F-Al 表示以橋氟離子連接的Al-F 絡合離子構型。從圖9 可見:隨著KF濃度增大,熔鹽中Al-F離子傾向于以更高配位如[AlF6]3-存在,而[AlF5]2-和[AlF4]-則逐漸降低,且[AlF4]-一直以一定比例降低,這與F-Al-F鍵角分布計算結果一致。此外,Al-F-Al 構型所占百分比逐漸增大,這表明隨著KF濃度增大,熔鹽的離子結構也逐漸復雜,更復雜的Al-F 配離子構型會降低Al 離子和F 離子的擴散性能,進而降低整個熔鹽體系的輸運性能。

圖9 KF-NaF-AlF3熔鹽體系中Al-F絡合離子類型分布規律Fig.9 Distribution of Al-F complex ion types in KF-NaFAlF3 system

3 結論

1)基于DFT理論構建了KF-NaF-AlF3低溫電解質體系的勢參數,并采用經典分子動力學模擬對KF-NaF-AlF3體系的微觀離子結構進行了研究。

2)KF-NaF-AlF3體系中主要的絡合離子構型為[AlF4]-,[AlF5]2-和[AlF6]3-,分別對應著扭曲的四面體、三角雙錐體和八面體構型。KF-NaF-AlF3熔鹽體系雖然失去了長程有序狀態,但局域離子結構仍然保持著短程有序。Al-F離子對的RDF在1.75×10-10m左右存在1個高而尖的峰值,這意味著Al-F離子間的相互作用較強,在熔鹽中易形成較復雜的絡合離子,并且這一峰值與Al-F 的平均鍵長相對應。

3)F-Al-F 的鍵角分布曲線的第一峰值和第二峰值分別位于91°和116°,分別對應著1 個八面體或三角雙錐體的鋁氟絡合離子集團結構,同時,位于173°左右的第三峰反映了三角雙錐體的特征。在108°左右出現第四峰,對應著正四面體離子集團結構,但其峰值較小,這說明KF-NaF-AlF3熔鹽中存在著少量的四配位[AlF4]-絡合離子集團結構。

4)由于橋氟離子的存在,熔鹽中形成了更復雜的絡合離子,如[Al2F8]2-,[Al3F10]-,[Al4F14]2-,[Al4F15]3-和[Al4F16]4-等。隨著KF 濃度增大,熔鹽中Al-F 離子傾向于以更高配位如[AlF6]3-存在,而[AlF5]2-和[AlF4]-的物質的量之比逐漸降低,且[AlF4]-一直以一定比例降低。更復雜的Al-F配離子構型會降低Al 離子和F 離子的擴散性能,進而降低整個熔鹽體系的輸運性能。

猜你喜歡
體系結構
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
構建體系,舉一反三
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
探索自由貿易賬戶體系創新應用
中國外匯(2019年17期)2019-11-16 09:31:14
新型平衡塊結構的應用
模具制造(2019年3期)2019-06-06 02:10:54
論《日出》的結構
創新治理結構促進中小企業持續成長
現代企業(2015年9期)2015-02-28 18:56:50
如何建立長期有效的培訓體系
現代企業(2015年1期)2015-02-28 18:43:18
“曲線運動”知識體系和方法指導
基于BIM的結構出圖
主站蜘蛛池模板: 99在线观看国产| 永久免费精品视频| 国产午夜人做人免费视频| 欧美成人精品高清在线下载| 国产亚洲视频在线观看| 免费激情网址| 亚洲欧美自拍中文| 国产精品流白浆在线观看| 美女扒开下面流白浆在线试听| 亚洲欧美日韩中文字幕一区二区三区| 亚洲AV永久无码精品古装片| 97亚洲色综久久精品| 被公侵犯人妻少妇一区二区三区 | 国产三级韩国三级理| 一本综合久久| 找国产毛片看| 久久91精品牛牛| 日韩精品一区二区三区swag| 色婷婷亚洲综合五月| 日韩精品无码一级毛片免费| 国产a网站| 国产激情无码一区二区免费| 午夜a级毛片| 国产欧美一区二区三区视频在线观看| 久久国产成人精品国产成人亚洲 | 波多野结衣在线se| 国产永久在线视频| 国产福利2021最新在线观看| 久草性视频| 精品国产污污免费网站| 国产精品自在在线午夜| 丰满人妻中出白浆| 在线欧美日韩| 国产毛片久久国产| 欧美午夜在线观看| 国产精品所毛片视频| 人与鲁专区| 巨熟乳波霸若妻中文观看免费| 成人精品在线观看| 免费观看国产小粉嫩喷水| 日本免费精品| 亚洲欧美日韩中文字幕一区二区三区| 亚洲国产日韩在线成人蜜芽| 激情网址在线观看| 久久91精品牛牛| 特级aaaaaaaaa毛片免费视频| 国产欧美日韩专区发布| 国产精品美女网站| 亚洲视频在线网| 毛片在线区| 亚洲精品成人7777在线观看| 国产丝袜第一页| 91视频精品| 久久99国产精品成人欧美| 视频一区亚洲| 国产女人综合久久精品视| 欧美成人第一页| 国产一区二区三区精品欧美日韩| 日韩无码视频播放| 日韩精品毛片人妻AV不卡| www.国产福利| 99热这里只有成人精品国产| 色有码无码视频| 五月婷婷导航| 国产成人无码AV在线播放动漫 | 丁香五月激情图片| 少妇被粗大的猛烈进出免费视频| 成人综合久久综合| 久久久久久久久久国产精品| 久久综合色播五月男人的天堂| 日韩在线影院| 免费在线观看av| 国产成人毛片| 亚洲伊人天堂| 亚洲男人在线| 毛片网站在线看| 精品夜恋影院亚洲欧洲| 国产综合亚洲欧洲区精品无码| 亚洲天堂视频在线观看免费| 欧洲高清无码在线| 男女性色大片免费网站| 欧美在线导航|