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

重質船用燃料油的瀝青質組分在柴油中的相行為模擬研究

2023-09-06 12:34:54師黎明
石油煉制與化工 2023年9期
關鍵詞:體系模型

師黎明,楊 鶴,任 強,李 妍,龍 軍

(中石化石油化工科學研究院有限公司,北京 100083)

重質船用燃料油(重質船燃)是一種用于大型遠洋船舶發動機的重質燃料,具有黏度大、熱值高等特性。重質船燃主要采用重質油與輕質油調合[1-2],然而受國際海事組織(IMO)限硫協議及其他諸多因素影響[3-4],重質船燃的調合原料來源變得更加廣泛,從石油煉制產物到深加工原料[5],甚至包括煤焦油、輪胎油等[6],然而,不同來源調合原料的組成存在明顯差異,常常會影響原料間的相容性以及各組分的分散性,從而加劇重質船燃中瀝青質的聚集與沉積,導致調合的重質船燃產品穩定性不佳[7-9]。

鑒于此,國內外學者嘗試使用不同的原料進行重質船燃調合試驗,優選出不同組分間的相容性較好、產品相對穩定的重質船燃調合方案。薛倩等[10]研究發現,在由煤焦油組分與瀝青調合的重質船燃中摻入C9或C10低凝點高芳香性組分后,產品中的不穩定組分明顯減少;Kass等[11]評估發現以生物質油與高黏度重質燃料(HFO)混合制備船用燃料油是可行的;Vráblík等[12]向由減壓渣油(VR)與常壓瓦斯油(AGO)調合的船用燃料油中加入輕循環油(LCO)或萘與四氫萘的混合物后,其加速老化沉積數值(TSA)均降低;Kondrasheva等[13]使用三角形相圖發現由VR、超低硫柴油(ULSD)和LCO可以調合得到穩定性較好的船用燃料產品。

雖然通過調合試驗探索可以獲得穩定性較好的重質船燃,但以上結果多基于經驗和實際操作獲得,對其他原料不具有普遍的適用性,對于新的調合原料仍需進行復雜的試驗探索。因此,從微觀結構和分子水平認知重質船燃調合組分間作用的規律性,對于重質船燃產品的調合生產具有重要的指導意義。然而,由于重質船燃的組成極其復雜,目前還無法通過試驗方法獲得其真實的組成和組分結構。近年來,人工智能模擬方法被廣泛用于重質油微觀結構認知及其分子間相互作用的研究[14],如張文等[15]通過分子動力學(MD)模擬發現不同結構VR分子間的相互作用不同,從而導致其微觀結構分布不均勻;關冬等[16]基于耗散粒子動力學(DPD)方法,對重質油體系瀝青質微觀聚集態進行模擬,將瀝青質聚集率與膠體不穩定指數(CII)進行關聯,提出了改進的重質油穩定性判據;Park等[17]采用MD方法對溶劑脫瀝青(SDA)過程進行了模擬,發現SDA過程中尺寸較大的聚集體一般是由多個小聚集體聚集形成;Ramírez等[18]對PA3型結構瀝青質分子在不同溶劑中的溶解行為進行了模擬研究,發現其分子中的雜原子能增加分子的極性和平面性,改善分子的聚集傾向。上述研究表明,借助分子模擬工具,可以深入探究重質船燃的微觀組成結構及分子間相互作用情況,從而提升對重質船燃調合原料匹配規律的認知。

基于此,本研究采用MD等模擬方法,對重質船燃組分瀝青質的分子結構、柴油調合原料的分子結構、不同分子間相互作用及組分間的相容性進行模擬計算,以提升對瀝青質分子在柴油中相行為的認知和理解。

1 模擬方法

1.1 重質船燃分子模型

重質船燃分子模型的構建是基于典型調合原料的烴類組成特點進行的。其中,瀝青質分子為塔河原油中瀝青質的平均分子結構[19],其平均分子式為C47H41NOS2,相對分子質量為699;柴油原料主要有LCO和加氫精制柴油(簡稱加氫柴油)。研究表明[20]:LCO中的飽和烴組分主要為鏈烷烴,而其芳香烴組分以多環芳烴為主,因此選擇與柴油主要組分特征相近的正十四烷、烷基萘及菲作為柴油模型分子;加氫柴油是由LCO選擇性加氫得到,因而其飽和烴組分含有較多的環烷烴,其芳香烴組分則主要是單環芳烴,因此選擇烷基十氫萘、烷基苯作為加氫柴油模型分子。重質船燃的各種模型分子結構如圖1所示。其中灰、白、紅、黃、藍色球分別代表碳、氫、氧、硫、氮原子。

圖1 重質船燃組分模型分子

1.2 量子化學計算方法

使用分子模擬軟件Materials Studio 2017(簡稱MS)中的DMol3量子化學模塊,優化模型分子的幾何結構,并計算分子結構的靜電勢分布特性。分子結構優化過程中選擇基于廣義梯度近似(GGA)的PBE泛函,在大數值基組DNP(雙數值軌道基組與p軌道極化函數)水平上進行全電子計算,自洽場(SCF)的迭代收斂閾值設置為1×10-5Ha,收斂精度為:能量2×10-5Ha,受力0.004 Ha/nm,位移5×10-4nm。對優化結果進行Potentials選項分析,得到分子結構的靜電勢分布。

1.3 蒙特卡洛抽樣計算方法

使用MS中的Blends模塊進行分子結合能的模擬計算:如圖2所示,將兩種不同結構分子(記作pack和seed)的質心進行重疊,之后移動pack分子,當二者的范德華表面不再重疊時,計算該空間構象下二者的結合能;再改變pack分子在空間的取向以及移動方向,重復上述計算過程,最終得到兩種分子在空間里不同二聚體形貌的集合,對集合進行平均得到分子間的平均結合能(Ebs)。計算過程中,集合樣本數取10 000,迭代次數選擇20,計算環境溫度選擇373.15 K,為重質船燃的調合溫度。

圖2 Pack分子(黃色)與Seed分子(藍色)的空間構象

1.4 分子動力學模擬方法

使用MS中的Amorphous Cell建立初始密度為0.4 g/L的無定形晶胞,并對晶胞進行結構優化與退火,選擇退火過程中能量最低的構象。對低能量構象使用MS中Focite模塊Dynamics任務項進行分子動力學模擬,先進行1 000 ps的NPT系綜模擬以使體系達到合理的密度,再進行2 000 ps的NVT系綜模擬使體系中的分子進行充分相互作用。當體系能量及其分項出現小范圍波動時,表示體系已達到平衡。過程中力場選擇COMPASSⅡ,計算精度為fine,靜電作用和范德華作用的非鍵截斷均采用Group based,溫度與壓力分別設置為373.15 K及標準大氣壓,且溫度與壓力的控制函數分別選擇Nose與Berendsen。對充分平衡后的動力學模擬結果使用Cohesive Energy Density(CED)任務項進行溶解度參數的計算,計算式見式(1)[21]。

(1)

式中:Ecoh為瀝青質內聚能,指將瀝青質體系中的所有分子移動到無窮遠處所需要的能量,kJ/mol;V為摩爾體積,cm3/mol;d為內聚能密度,J/cm3;δ為溶解度參數,(J·cm-3)0.5;ΔHv為瀝青質分子氣化熱,kJ/mol;RT表示瀝青質分子聚集體由固體或液體轉化為氣體時所需要的膨脹功,kJ/mol。

使用Energy任務項進行不同組分間平均作用能的計算,計算式見式(2):

Ea-b=Etotal-Ea-Eb

(2)

式中:Ea-b為a、b組分間的平均作用能,kJ/mol;Etotal為a、b組分共存時整個體系的總能量,kJ/mol;Ea與Eb則分別為a、b組分單獨存在時體系的能量,kJ/mol。

對部分動力學模擬結果會進行徑向分布函數(RDF)分析。RDF的意義是反映從一個任意指定的“中心”粒子出發,到半徑為r的球殼層位置上出現其他粒子的概率[22],如圖3所示。

圖3 徑向分布函數圖示

本研究中以某一粒子為觀測中心,以0.1 nm為球半徑差(Δr)向外畫出一層層的同心球殼,然后統計每一殼層里的粒子數密度與平均數密度的比值,即為徑向分布函數,用g(r)表示,如式(3)。

(3)

式中:n(r)表示總粒子數;r為距中心原子半徑,nm;V為體積,nm3;ρ0為理想晶體粒子數密度(ρ0=n(r)/V),δr為球殼層厚度,nm。RDF曲線越高,表明中心粒子周圍越傾向于出現某種粒子。

2 結果與討論

2.1 模型分子的靜電勢分布

對不同的重質船燃模型分子進行靜電勢分布模擬計算,過靜電勢分布推測分子的各種性質以及分子間的相互作用。各種模型分子的靜電勢模擬計算結果如圖4所示。

圖4 模型分子的靜電勢分布

計算結果顯示,不同來源的重質船燃模型分子,因其結構差異所呈現的靜電勢分布各不相同。瀝青質分子因具有大芳環結構和多雜原子的結構特征,使得π電子云主要集中在芳環上,瀝青質分子間因π電子云重疊而產生較強的π-π相互作用,在空間上會呈現出面對面式的堆疊形貌;而其芳環周圍的H原子則帶部分正電荷,可以與π電子云形成π-H作用使得分子間呈T型堆疊;瀝青質分子雜原子芳環結構中的N原子以及芳環側鏈上的羥基O原子帶有負電荷,與羥基O原子相連的H原子帶部分正電荷,使得瀝青質分子間易形成O—H、N—H等類型的氫鍵作用。較小芳香烴模型分子的芳環結構中也因p軌道的肩并肩形成了大π鍵,π電子云在該結構中富集,且芳環數越多富集區域越大。芳香烴分子結構中的長烷基側鏈以及飽和烴結構分子電荷分布較為均勻,分子間主要是范德華相互作用。對比不同結構分子間的靜電勢分布可以看到,在結構上芳香烴分子與瀝青質分子較為相似,分子間容易形成π-π相互作用,而飽和烴結構分子與瀝青質分子存在明顯的結構差異,分子間可能主要形成范德華相互作用。

2.2 模型分子間的結合能

通過結合能的模擬計算可量化比較不同結構模型分子間的相互作用強弱,因此計算了瀝青質分子自身及其與其他不同結構模型分子間的結合能,結果如圖5所示。由圖5可知,瀝青質分子間具有很強的相互作用,最左側黑色柱狀圖數值較高,在373.15 K下的分子間結合能為118.18 kJ/mol,因而其表現出很強的自聚傾向。不同結構的模型分子與瀝青質分子的結合能存在明顯的差異,其從高到低的順序為菲>烷基萘>烷基苯>烷基十氫萘>正十四烷。由此可見:柴油原料中多芳環結構的菲、烷基萘分子與瀝青質分子的結合能較高,表明芳環間的π-π相互作用能較高,這是因為π電子云重疊面積越大,π電子的離域范圍越大,致使體系能量降低、穩定性提高;而含長烷基側鏈的烷基苯分子及飽和烴分子與瀝青質分子的結合能較低,可能的原因是長烷基側鏈結構和飽和烴長鏈不利于接近瀝青質分子,分子間形成的范德華作用相對較弱。

圖5 重質船燃模型分子間的結合能

2.3 模型分子的溶解度參數

雖然分子結合能能夠對單個分子間的相互作用強弱進行量化描述,但不能反應不同分子間的相行為。溶解度參數(δ)是衡量分子間相容性的重要參數,是對“相似相溶”原理的量化表征,因而通過溶解度參數的模擬計算可以預測重質船燃中不同結構分子間的相容性,能在一定程度上反映分子在溶液中的相行為。分子的溶解度參數越接近,表明分子間的相容性越好,分子在溶液體系中越不易分相。

對于由50個相同模型分子構成的無定形晶胞體系,在373.15 K下進行動力學模擬,并在體系平衡后進行內聚能密度計算,將計算結果取平方根后即得到該種物質的溶解度參數。為驗證上述溶解度參數計算方法的可靠性,采用相同的模擬方法對3種已知小分子進行密度(ρ)、溶解度參數的模擬計算,結果如表1所示。

表1 3種小分子的密度、溶解度參數模擬計算結果[21]

由表1可知,由上述模擬得到的計算值與試驗值[21]相近,表明該模擬計算方法具有良好的可靠性。進一步計算瀝青質分子及柴油原料分子的溶解度參數,結果如表2所示。

表2 模型分子密度、溶解度參數模擬計算結果

由表2可知,不同結構分子的溶解度參數各不相同,其中瀝青質分子的溶解度參數最大,為22.618(J·cm-3)0.5,與之較為接近的是LCO原料中的菲和烷基萘分子,而正十四烷與加氫柴油原料中的烷基十氫萘、烷基苯分子的溶解度參數最為接近,為16~18(J·cm-3)0.5。同時,正十四烷、烷基十氫萘以及烷基苯分子的溶解度參數與瀝青質分子相差較大,表明這3種分子與瀝青質分子的相容性較差。結合上述分子結構中的靜電勢分布及單分子間的結合能分析,菲與烷基萘在結構上與瀝青質分子較為相似,其分子間因π電子云重疊能形成較強的π-π相互作用,因而其在混合時更容易互溶。飽和烴分子與瀝青質分子存在明顯的結構差異,其分子間僅有較弱的范德華作用,二者混合時更傾向于分相;而烷基苯因其分子結構中較長烷基側鏈的影響,其芳環結構不易接近瀝青質分子,二者混合也存在相分離的趨勢。

2.4 二元模擬溶液

通過單分子性質的模擬計算,對重質船燃分子結構特征、分子間作用強弱及互溶性有了深入的了解。然而,溶解度參數僅可用于對不同結構重質船燃模型分子間互溶性的預測,并不能真實反映分子在溶液中的相行為。同時,上述研究還發現,瀝青質分子間具有較強的相互作用,可能會在體系中聚集,因此在重質船燃調合過程中瀝青質組分在整個體系中的相行為更受關注。基于此,構建了以瀝青質分子為溶質分子、其他某種模型分子為溶劑分子的二元模擬溶液模型,進一步在373.15 K、標準大氣壓下探究瀝青質與不同柴油組分的互溶性,考察瀝青質分子在柴油溶液中的分子分布。

2.4.1瀝青質模型

在構建二元溶液模型前,為模擬計算不同大小瀝青質分子模型體系的動力學特性,分別構建了分子數為10,20,30,40,50的5個瀝青質分子一元體系模型,并分別對其結構進行優化并退火。經動力學模擬體系平衡后,對5個瀝青質分子模型體系分別進行密度、溶解度參數、內聚能的計算,結果如表3所示。

表3 不同大小瀝青質分子一元模型體系性質的計算結果

各瀝青質分子模型體系平衡時,瀝青質分子的分子分布如圖6所示(隱藏了晶格)。從截取的平衡軌跡圖像可知:瀝青質分子間具有較強的聚集傾向,模型中形成了尺寸較大的多聚體結構,聚集體主要以π-π堆疊形貌為主;而且,隨著瀝青質分子數量增多,π-π堆疊尺寸逐步增大,瀝青質聚集體的內聚能也逐漸增大,這表明內聚能可以反映瀝青質的團聚程度。

圖6 瀝青質分子一元模型體系的分子分布

2.4.2二元溶液模型的動力學模擬

將瀝青質分子分別與正十四、烷基十氫、烷基苯、菲的分子混合,構建二元模擬體系,并調整溶劑的質量分數(下同)分別為30%,50%,70%,考察分子間作用強弱與溶劑含量對瀝青質分子在溶液中相行為的影響。截取12個溶液模型動力學模擬平衡時的軌跡圖像,分別如圖7所示。

圖7 含瀝青質分子與柴油組分分子的二元溶液體系動力學模擬平衡時的軌跡圖像

由圖7可知,對于瀝青質-飽和烴溶液模型,飽和烴分子正十四烷、烷基十氫萘與瀝青質分子不相溶,模擬體系中瀝青質分子發生聚集形成了尺寸較大的聚集體;即使增加飽和烴的量,體系中的瀝青質分子依舊發生明顯的聚集,并且有自成一相的趨勢,表明瀝青質分子在飽和烴中的分散性較差。而在瀝青質-菲溶液模型中,當菲含量較低時,模型體相中存在瀝青質分子的多聚體,但其聚集尺寸相較于瀝青質-飽和烴體系較小;增加體系中菲的量,體相中聚集體的尺寸明顯減小,僅局部存在瀝青質分子的二聚體或三聚體,其他瀝青質分子以單分子的形態分散于體相中。對于瀝青質-烷基苯溶液模型,當長鏈烷基苯溶劑的質量分數升至70%時,體相中仍存在尺寸較大的瀝青質聚集體,說明瀝青質分子在長鏈烷基苯中的分散性比在瀝青質-菲溶液體系差。

結合模型分子的模擬計算結果可知:飽和烴分子與瀝青質分子間的相互作用弱于含芳環結構烴分子與瀝青質分子間的相互作用,導致瀝青質分子不能很好地被飽和烴溶劑分子膠溶分散;芳香烴分子與瀝青質分子雖然具有良好的相互作用,但當芳香烴量較少時,其無法很好地膠溶分散瀝青質分子,同時含有長烷基側鏈的芳香烴分子對瀝青質分子的膠溶分散效果也不佳。這說明,當分子間相互作用較弱時,二元溶液體系中兩種組分不相溶且有明顯的相分離趨勢;當分子間相互作用較強時,若溶劑分子足夠多,就能很好地膠溶分散極性瀝青質分子。

對溶劑質量分數均為70%的瀝青質-正十四烷、瀝青質-菲、瀝青質-烷基十氫萘、瀝青質-烷基苯二元溶液模型進行RDF分析,對比瀝青質分子周圍出現溶劑分子概率的大小,結果如圖8所示。

圖8 4種二元溶液模型體系的RDF對比

由圖8可知,在相同半徑下,瀝青質-菲體系的RDF最大,瀝青質-烷基苯體系次之,瀝青質-飽和烴體系最小,表明芳環結構分子更容易接近瀝青質分子,而且芳環數越多、烷基側鏈越短,則在瀝青質分子周圍出現的概率越大。此外,瀝青質-菲曲線最大RDF峰值對應的半徑為0.329 nm,依據相關研究[23-24],該半徑在含芳環組分π-π相互作用范圍內。圖9為不同模型分子間堆疊形式和相互作用示意。如圖9所示:菲分子與瀝青質分子形成了面對面式的π-π堆疊相互作用;烷基苯分子也與瀝青質分子形成了π-π堆疊相互作用,但其長烷基側鏈不利于其與瀝青質分子接近,導致其對應的RDF峰值比菲小;飽和烴分子僅能與瀝青質分子的烷基側鏈發生較弱的范德華作用,在模擬溶液體系中與瀝青質分子距離較遠。

圖9 模型分子間的相互作用

RDF分析可以對溶液中的分子分布進行定性分析,而對模擬溶液中分子分布的量化分析則需計算不同組分之間作用能。組分間的作用能既包含單分子間的作用能,也要結合分子數量的影響,從而反映出模擬溶液中組分間的相溶性情況。二元模擬溶液體系中瀝青質與溶劑組分間的作用能計算結果如表4所示。

表4 不同模擬溶液中組分間的作用能模擬計算結果 kJ/mol

由表4可知:當溶劑含量較低時,組分間作用能在較低水平;隨著溶劑含量升高,不同二元體系的組分間作用能均有增長。相比于瀝青質-菲溶液體系,飽和烴分子溶液中的組分間作用能增幅較小,即使飽和烴質量分數達到70%,組分間作用能仍維持在較低水平;瀝青質-烷基苯溶液體系中的組分間作用能也相對較小,可能是因為烷基苯中長烷基側鏈不利于其與瀝青質分子間發生相互作用。隨著瀝青質-菲溶液體系中菲含量的增加,其組分間作用能顯著升高。

將以上二元溶液體系組分間作用能與瀝青質一元體系內373.15 K下的內聚能(-2 765.116 kJ/mol)進行比較,可以發現:在瀝青質-飽和烴、瀝青質-烷基苯二元溶液體系中,即使溶劑質量分數達到70%,其組分間作用能仍比瀝青質的內聚能小;而在瀝青質-菲二元溶液體系中,當溶劑質量分數達到70%時,組分間作用能大于瀝青質的內聚能。由于組分間作用能與瀝青質內聚能的相對大小和體系中分子分布有較好的對應關系,由此可知飽和烴、烷基苯難以膠溶分散瀝青質,而質量分數70%的菲作為溶劑能夠很好地膠溶分散瀝青質。

3 結 論

(1)不同來源的重質船燃分子,因其分子結構不同致使其電勢分布各不相同,從而導致不同分子間的相互作用也各不相同。柴油調合原料中的烷基苯、烷基萘、菲等芳香烴分子結構與重質船燃中瀝青質組分分子的結構較為相似,因而其與瀝青質分子間能形成較強的π-π相互作用,分子間結合能較高、溶解度參數相近;飽和烴分子與瀝青質分子結構差異較大,其分子間僅有較弱的范德華力作用,二者的溶解度參數相差較大。

(2)采用MD方法對瀝青質-溶劑二元溶液體系模擬計算發現:菲與瀝青質組分間相互作用較強,當其在二元體系中質量分數達70%時,可使瀝青質膠溶分散且不易相互聚集;含長烷基側鏈的烷基苯因側鏈的影響,導致其與瀝青質組分間作用能大幅減小,不利于對瀝青質分子的膠溶分散;鏈烷烴或環烷烴均與瀝青質分子作用較弱,即使含量很高,也無法膠溶分散瀝青質分子。

(3)373.15 K下瀝青質的內聚能為-2 765.116 kJ/mol。為使瀝青質膠溶分散,不僅溶劑與瀝青質組分間的相互作用要強,溶劑含量也必須足夠高。原因在于:組分間作用能隨著溶劑含量的增加而增大,當組分間作用能大于瀝青質內聚能時,才能實現對瀝青質的膠溶分散,且二者差值越大,瀝青質-溶劑的膠溶分散體系越穩定。

猜你喜歡
體系模型
一半模型
構建體系,舉一反三
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
探索自由貿易賬戶體系創新應用
中國外匯(2019年17期)2019-11-16 09:31:14
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
如何建立長期有效的培訓體系
現代企業(2015年1期)2015-02-28 18:43:18
“曲線運動”知識體系和方法指導
“三位一體”德育教育體系評說
中國火炬(2010年7期)2010-07-25 10:26:09
主站蜘蛛池模板: 国产香蕉在线视频| 天天干天天色综合网| 毛片最新网址| 亚洲欧美日本国产综合在线| 中文字幕在线日本| 中文字幕无码av专区久久| 99精品一区二区免费视频| 亚洲人成在线精品| 久久91精品牛牛| 日日拍夜夜操| 亚洲欧洲国产成人综合不卡| 国产精品自拍露脸视频| 久久综合九色综合97婷婷| 综合人妻久久一区二区精品 | 99视频在线看| 波多野结衣无码视频在线观看| 一级片免费网站| 亚洲无码视频一区二区三区 | 91精品国产丝袜| 乱系列中文字幕在线视频| 丁香五月激情图片| 国产精品亚欧美一区二区三区| 中文字幕不卡免费高清视频| 国产麻豆另类AV| 91无码视频在线观看| 日韩无码视频专区| jijzzizz老师出水喷水喷出| 91 九色视频丝袜| 2021亚洲精品不卡a| 自偷自拍三级全三级视频| 国产精品白浆在线播放| 欧美日本视频在线观看| 欧美色香蕉| 日韩成人在线网站| 国产永久无码观看在线| 国产欧美网站| 国产成人精品男人的天堂| 麻豆国产原创视频在线播放| 91丝袜美腿高跟国产极品老师| 丁香婷婷在线视频| 无码人中文字幕| 亚洲国产精品国自产拍A| 日本精品视频一区二区| 欧美成人国产| 91精品啪在线观看国产60岁| 伦伦影院精品一区| 老司国产精品视频91| 波多野结衣中文字幕一区| 欧美日韩免费在线视频| 57pao国产成视频免费播放| 香蕉网久久| 免费在线一区| 亚洲五月激情网| 国产亚洲欧美另类一区二区| 久久无码av三级| 国产福利小视频在线播放观看| 1024你懂的国产精品| 精品少妇人妻一区二区| 欧美中文字幕第一页线路一| 国产成人1024精品下载| 中文字幕在线日本| 91在线精品免费免费播放| V一区无码内射国产| 伊人天堂网| 久久永久精品免费视频| 精品福利视频导航| 伊人色婷婷| 日本在线视频免费| 最新国语自产精品视频在| 亚洲精品自拍区在线观看| 激情综合图区| 亚洲午夜国产精品无卡| 亚洲精品午夜天堂网页| 国产女人18水真多毛片18精品 | 久久美女精品| 免费国产一级 片内射老| 成人免费视频一区| 欧美色亚洲| 亚洲人成网站在线播放2019| 亚洲va视频| 99草精品视频| 99久久精品久久久久久婷婷|