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

NaX分子篩吸附天然氣中酸性氣的分子模擬

2023-09-06 12:34:56魯家榮劉熠斌陳小博
石油煉制與化工 2023年9期
關鍵詞:擴散系數模型

王 翀,魯家榮,閆 昊,劉熠斌,陳小博

(1.中石化石油工程設計有限公司,山東 東營 257026;2.中國石油大學(華東)重質油國家重點實驗室)

天然氣的主要成分為CH4,一般還存在部分的CO2、H2S等酸性氣。酸性氣的存在不僅會在水蒸氣的作用下腐蝕設備和管路,而且會影響天然氣的品質及后續加工過程。因此,天然氣中酸性氣的脫除是天然氣工業的重要過程。目前酸性氣主要脫除技術為濕式化學脫除,具體做法是采用堿性液體作為吸收劑,使其與天然氣中的H2S和CO2發生化學反應以達到脫除酸性氣的目的[1]。盡管濕式化學法脫除酸性氣的效率較高,但能耗高,且存在溶劑廢液污染。近年來節能環保型酸性氣脫除技術有了長足的發展,如出現了分子篩變壓吸附技術、膜分離技術,其中分子篩變壓吸附技術具有脫除率高、能耗低、無腐蝕、操作簡單等優點。

研究天然氣中不同組分在分子篩中的吸附擴散行為,對于優化分子篩設計、提升酸性氣脫除能力具有重要的指導意義。分子模擬是通過理論計算模擬分子結構與行為的方法,可用于多孔材料內的吸附與擴散研究。Maurin等[2]以不同離子交換的X分子篩為研究對象,采用巨正則蒙特卡羅(GCMC)方法結合根據離子本征特性導出的力場模型進行了模型分子吸附行為的模擬,計算的等溫線和微分等量熱與試驗結果具有良好的一致性。Demir等[3]采用GCMC方法對CH4、CO2體系在類沸石金屬有機骨架(ZMOFs)中的吸附行為進行了研究,發現在二元吸附體系下CH4和CO2在ZMOFs中的吸附位點相同,ZMOFs對CO2的吸附優于對CH4的吸附。Sung Chunyi等[4]通過分子模擬方法研究了多價陽離子交換的Y分子篩對克勞斯反應尾氣中微量H2S的選擇性吸附,并通過密度泛函理論(DFT)方法計算了不同材料對H2S的吸附能,提出使用GaY分子篩進行H2S的吸附是有效的脫硫策略。丁雪等[5]采用GCMC方法研究了干氣中不同組分在ZSM-5分子篩中的吸附行為,并計算了吸附過程的焓變、熵變等熱力學數據。黨宇等[6]則將GCMC方法與分子動力學方法結合,研究了壓力對噻吩、吡咯、呋喃3種雜環化合物在HY分子篩中的吸附與擴散性能,獲得了概率密度、徑向分布函數、吸附能等。

作為一種具有廣泛應用前景的X分子篩,NaX分子篩被廣泛應用于天然氣中酸性氣的脫除,但是對天然氣中不同組分在NaX分子篩中競爭吸附擴散行為的分子模擬,尤其是針對同時含有CO2、H2S的體系的研究則鮮有報道。本課題采用GCMC與分子動力學結合的方法,研究3個溫度(273,293,313 K)下天然氣中CH4,H2S,CO2在NaX分子篩中的吸附與擴散行為,為理解天然氣吸附脫硫、脫碳的微觀行為和吸附材料的理論設計提供指導。

1 計算模型與方法

1.1 分子篩吸附劑模型

X分子篩是一種硅鋁比為1.0~1.5的多孔硅酸鹽材料,其拓撲結構為FAU型,屬于立方晶系,空間群為Fm-3m,晶胞參數a=b=c=2.474 nm,α=β=γ=90°,籠口孔徑為0.74 nm[7]。13X分子篩的晶胞由β籠和六方柱籠構成,β籠與β籠之間通過六方柱籠連接,形成窗口直徑為0.74 nm的十二元環和對應的超籠結構。

本研究采用Materials Studio 8.0軟件平臺,分子篩的空間群和晶胞參數如上文所述。為簡化模型,在13X分子篩的建模中采用硅鋁比為1的結構。由于分子篩骨架為負電荷,因此引入Na+進行電荷平衡,平衡后的晶胞組成為:Na96Al96Si96O384。模型建立后,根據文獻[8-9]分別為各種原子進行電荷分配:Na(+0.70),Al(+1.75),Si(+2.35),O(-1.20)。結構優化后的NaX分子篩結構模型如圖1所示。

圖1 NaX分子篩的結構模型

采用Materials Studio 8.0軟件中分子動力學模塊對NaX分子篩結構進行優化,并獲得最低能量態。為了驗證晶體模型的合理性,使用Reflex模塊模擬了該結構的X射線衍射圖譜(圖2)。將模擬結果與數據庫[10]中相應結構的標準圖譜(圖2中FAU)對比,發現兩者具有相同的特征峰,證明結構優化后的晶體類型為FAU結構。

圖2 優化結構的XRD模擬結果與數據庫中FAU結構特征峰的對比

NaX分子篩的結構由帶有負電性的硅氧鋁骨架結構和陽離子構成,其中陽離子的類型以及分布對于其吸附性能具有重要影響[11]。Frising等[12]詳細總結了不同文獻中合成的NaX分子篩的陽離子分布情況。由于各種文獻對Na+位點的定義繁雜并且Na+在分子篩中的分布較為復雜,在總結各種文獻后給出Na+一般分布情況:Na+在骨架中主要分布在如圖3所示的3種位點,分別是Ⅰ(方鈉石籠,另有Ⅰ′位點與該位點十分接近,統稱為Ⅰ位點)、Ⅱ(方鈉石籠和超籠之間的六邊形窗口中心)、Ⅲ(超籠中靠近中心四邊形窗口,另有超籠中靠近周邊四邊形窗口的位點Ⅲ′與之十分接近,文獻中對Ⅲ和Ⅲ′的界定較為模糊,故統稱為Ⅲ位點)。不同文獻所報道的Na+分布情況如表1所示。將結構優化后的NaX分子式模型與表1的數據相對比,發現優化模型具有上述4種陽離子分布位點,位點Ⅰ,Ⅱ,Ⅲ的個數分別為24,30,31個,與文獻中的結果近似,證明模型具有合理性。

表1 不同文獻所報道的NaX分子篩骨架外陽離子分布情況

圖3 NaX分子篩的陽離子分布情況

1.2 吸附質分子模型

選取天然氣中的CH4,CO2,H2S作為吸附質分子,3種分子均采用Materials Studio軟件建模并采用Dmol 3量子力學模塊優化,獲得其ESP電荷(見表2)。3種分子的分子尺寸、四極矩等數據如表3所示。由表3可知,3種分子的動力學直徑均小于NaX分子篩的孔口直徑,因此3種分子均能擴散進入分子篩孔道進行吸附。模擬中將調節不同分子的摩爾比,探究其在分子篩模型上的吸附行為。

表2 天然氣不同組分分子中各原子的計算ESP電荷

表3 不同分子的性質

1.3 模擬參數設置

NaX分子篩對吸附質分子的吸附模擬采用Sorption模塊,抽樣方法選用Metropolis方法,精度選擇Ultra-fine。客體分子和分子篩的相互作用包括靜電作用和范德華作用,其中范德華作用采用Lennard-Jones勢函數描述,其數學表達式為:

(1)

式中:ULJ表示客體分子和吸附劑之間相互作用能,eV;i和j表示不同原子;Rij表示原子間距,nm;Dij和(R0)ij為Lennard-Jones參數;qi和qj表示原子所帶電荷,e。模擬中采用COMPASSⅡ力場,靜電相互作用采用Ewald方法處理,非鍵相互作用的加和采用Atom based算法,范德華作用截斷距離設置為1.25 nm,正好小于晶胞邊長(2.502 8 nm)的一半。計算平衡步數為1×106步,生產步數為1×107步。

吸附質分子的Dmol 3結構優化參數:精度選擇Fine,選用廣義梯度近似泛函(GGA)中的PBE近似,DFT-D色散校正選用TS方法。核處理選用全部電子,基組選用DNF數值基組。

分子在分子篩中的擴散采用Forcite模塊,根據牛頓力學的原理進行分子動力學計算。模擬中選用正則系綜(NVT),采用Velocity-Verlet算法進行積分,模擬步長1 fs,模擬總時長200 ps,其中前50 ps用于平衡體系。擴散系數Ds采用Einstein方程計算:

(2)

式中:Ds為擴散系數,m2/s;N表示單位晶胞內的分子數;r(t)是吸附質分子在t時刻的質心坐標;r(0)是吸附質分子質心的初始位置坐標;t表示時間,ps;平方項是擴散分子均方根位移(MSD)的系綜平均。

2 結果與討論

2.1 吸附等溫線和吸附熱

為了分析3種組分在NaX上的吸附強度,采用GCMC方法模擬計算了工業天然氣變壓吸附脫除酸性氣條件下,CH4,CO2,H2S在不同溫度(273,293,313 K)下的吸附等溫線數據(見圖4中散點),并采用Freundlich吸附模型對3種分子在3個溫度下的吸附等溫線數據進行擬合,擬合曲線見圖4中線條。

圖4 純組分CH4,CO2,H2S在NaX分子篩中的吸附等溫線

由圖4可知:3種物質的吸附等溫線均為Ⅰ型等溫線,表示3種分子在相應條件下在NaX分子篩上的吸附形式均為單分子層吸附;隨著平衡壓力的增加,3種組分的吸附量均表現出先增加后趨于飽和的規律;在壓力較低(0~2 MPa)時,CH4的吸附量增加相對緩慢,且在7 MPa左右達到吸附飽和;相同壓力范圍內H2S和CO2的吸附量迅速增加,在2 MPa附近即達到飽和。

Freundlich吸附模型適用此模擬的壓力范圍,其吸附量計算式如式(3)所示。

(3)

式中:N為吸附量,mmol/g;Nm為單層飽和吸附量,mmol/g;K為吸附平衡常數,MPa-1,數值上等于吸附速率常數與脫附速率常數之比,表征分子篩對某物質的吸附強弱;1/n為指數,表示濃度對吸附強度的影響,其值一般小于1;P為模擬條件下的壓力,MPa。不同溫度下CH4,CO2,H2S分子在NaX分子篩中吸附的Freundlich模型擬合參數分別如表4~表6所示。

表4 不同溫度下CH4分子在NaX分子篩中吸附的Freundlich模型擬合參數

表5 293 K下3種分子在NaX分子篩孔道中的自擴散系數

表5 不同溫度下CO2分子在NaX分子篩中吸附的Freundlich模型擬合參數

表6 不同溫度下H2S分子在NaX分子篩中吸附的Freundlich模型擬合參數

由表4~表6可以看出:Freundlich吸附模型可以很好地描述3種分子在NaX分子篩上的吸附等溫線,擬合方程的決定系數R2均大于0.99。表4~表6中的Nm,K,1/n分別為根據Freundlich吸附模型擬合得到的數據,擬合的Nm結果與采用GCMC方法模擬得到的單層飽和吸附量十分接近,表明Freundlich吸附模型可以較為合理地描述3種分子在分子篩中的單層吸附情況。

采用GCMC方法模擬得到在溫度273,293,313 K下,CH4的飽和吸附量分別為6.79,6.43,6.05 mmol/g,CO2的飽和吸附量分別為9.13,8.74,8.42 mmol/g,H2S的飽和吸附量分別為9.84,9.45,9.23 mmol/g,可見NaX分子篩對H2S和CO2的吸附量高于對CH4的吸附量。另外,3種分子的飽和吸附量均隨著溫度的降低而升高,表明高壓和低溫有利于CH4,CO2,H2S在NaX分子篩上的吸附。在擬合壓力范圍內,3種分子在NaX分子篩中吸附的平衡常數由大到小的順序為H2S>CO2>CH4,說明NaX對3種分子的吸附強度由大到小的順序為H2S>CO2>CH4,表明NaX分子篩對于天然氣中的酸性氣分子具有較強吸附能力。3種分子在該溫度下的微分吸附熱也具有上述規律。由于CH4分子的偶極矩和四極矩為零,CH4與NaX分子篩孔道的作用力主要為范德華力,而CO2和H2S分子除了上述作用外,還與分子篩骨架外Na+產生較強的靜電相互作用。H2S由于具有較大的偶極矩和較強的四極矩以及極化率,其與骨架外陽離子靜電相互作用的強度高于CO2,因此出現低壓吸附平衡常數和飽和吸附量由大到小的順序H2S>CO2>CH4。Maghsoudi等[21]通過試驗的方法研究了天然氣中CH4,CO2,H2S在CHA分子篩中吸附和擴散行為,結果表明CHA分子篩在CH4的存在下可以選擇性地分離H2S和CO2,其測定的吸附等溫線對飽和吸附量的結論與本研究結論一致。

在獲得3種分子在NaX中的吸附等溫線的基礎上,采用GCMC模擬方法得到了3種分子等量吸附熱隨著壓力變化的圖線。等量吸附熱ΔQ的計算參見文獻[6]。

由于溫度對等量吸附熱的影響很小,故這里僅給出溫度273 K下純物質CH4,H2S,CO2的吸附熱曲線(如圖5所示)。由圖5可知,3種分子的吸附熱由大到小的順序為H2S>CO2>CH4,表明分子篩對CO2和H2S的吸附作用較強,與吸附活性位作用后產生了較為明顯的熱效應,與吸附等溫線的模擬結果相符合;另外,3種分子等量吸附熱均隨著吸附壓力的增加而降低,最終達到穩定,這與很多小分子的吸附呈相似的規律[6]。

圖5 273 K下純組分在NaX分子篩中的吸附等量熱曲線

2.2 吸附位分析

根據物質在不同溫度下的吸附相互作用能量分布曲線可以得到吸附質分子在分子篩中的吸附位點分布,獲得吸附質的最佳吸附位置,探究溫度對吸附位的影響[13]。圖6為CH4,H2S,CO2分子在273,293,313 K下與NaX分子篩的相互作用曲線。由圖6可知,3個溫度下3種分子的吸附峰的位置基本不發生改變,說明一定溫度范圍內溫度對吸附位點的影響較小。CH4吸附能量分布曲線出峰位置在-32 kJ/mol和-20 kJ/mol,且在-20 kJ/mol處為強峰;CO2吸附能量分布曲線出峰位置在-42 kJ/mol和-20 kJ/mol兩處,在-42 kJ/mol處為強峰;H2S吸附能量分布曲線只有1個出峰位置,為-45 kJ/mol。可見3種分子在分子篩內主要有3種吸附位點,其對應的吸附能量分別為-20,-32,-42 kJ/mol左右,分別對應NaX分子篩中孔道中的活性位點Ⅱ,Ⅰ,Ⅲ。CH4主要吸附在Ⅱ位,H2S和CO2主要吸附在Ⅲ位。

圖6 不同分子與分子篩間的相互作用曲線

以溫度293 K為例,根據模擬計算得到了3種分子在該溫度下的吸附概率密度勢能場分布,如圖7所示(紅色代表分子在該處出現的概率高,相互作用強;綠色代表出現的概率低,相互作用弱)。由圖7可以看出,CH4分子主要分布在方鈉石β籠中以及方鈉石籠和超籠之間的交界;CO2分子和H2S分子則主要分布在超籠中,方鈉石籠中的分布量較少。從吸附質與分子篩相互作用上分析,可以看出CH4分子與NaX骨架的相互作用力較小,且相互作用較強的位置位于β籠中;CO2和H2S與骨架的相互作用更強,且主要作用位點在超籠中的Na附近。CO2在NaX中的分布情況與作用位點與Cui Yongkang等[14]的研究結果相一致。

圖7 293 K下不同分子的吸附密度勢能場分布

2.3 混合組分的競爭吸附

為了探究實際混合天然氣體系在分子篩中的吸附行為,首先構建了3種不同CO2含量的混合氣體,其中H2S的摩爾分數固定為3%,CO2的摩爾分數分別為10%,15%,20%。293 K下的吸附等溫線如圖8所示。由圖8可以看出,盡管混合氣體中CH4的含量占主導,但是NaX分子篩對CH4的吸附量遠低于對CO2和H2S的吸附量,體現出NaX分子篩具有較高的酸性氣脫除選擇性。隨著CO2含量的增加,CO2的飽和吸附量逐漸增加并且超過H2S的飽和吸附量。產生這種現象的原因可能是兩者共同競爭骨架的Ⅲ號吸附位,隨著CO2分子數增加,其占據了更多的吸附位,并獲得了更高的吸附選擇性。隨著CO2摩爾分數從10%增加到15%,CO2的飽和吸附量增加了18.2%;而當CO2摩爾分數從15%增加到20%后,其飽和吸附量僅增加了7.9%。說明隨著CO2分子數的增加,吸附選擇性從CO2濃度主導轉變為由H2S、CO2分子與吸附劑的相互作用主導,因此CO2吸附量增加緩慢。

圖8 293 K下CO2含量不同的混合氣體中各組分分子的吸附等溫線

將CO2的摩爾分數固定為3%,調節CH4和H2S的比例,控制H2S摩爾分數分別為3%,6%,9%,研究293 K下H2S含量對分子吸附行為的影響,獲得的等溫線如圖9所示。由圖9可以看出,分子篩對CH4分子的吸附量仍遠小于對H2S分子和CO2分子的吸附量。隨著H2S含量的增加,H2S和CO2的飽和吸附量均變化較小,且前者的飽和吸附量遠高于后者,體現出NaX分子篩對H2S具有較高的吸附選擇性。

圖9 293 K下H2S含量不同的混合氣體中各組分分子的吸附等溫線

2.4 擴散過程分析

客體分子在吸附劑中的擴散行為是吸附過程中的重要性質。由于CH4,CO2,H2S分子的動力學直徑均小于NaX分子篩的孔道直徑,因此3種分子的擴散性質主要受到分子的偶極矩以及極性的影響。擴散系數是表征擴散能力和計算傳質速率的重要參數。本研究采用分子動力學的方法計算了293 K下CH4,CO2,H2S在NaX分子篩中的擴散軌跡,并進一步得到了3種分子的均方根位移(MSD)隨時間的變化曲線,如圖10所示。

圖10 293 K下3種分子在NaX分子篩孔道中擴散的均方根位移

根據式(2)Einstein方程計算CH4,CO2,H2S在分子篩中的擴散系數,結果如表5所示。由表5可知,3種分子在分子篩中的自擴散系數由大到小的順序為CH4>H2S>CO2。由吸附性質計算可知,由于H2S和CO2與NaX分子篩表面具有較強的相互作用,因此較難在孔道中擴散。CH4與NaX分子篩表面僅有范德華相互作用,因此在孔道中的擴散阻力最小。擴散系數的結論與吸附性質的計算結果具有較好的符合度。

為了進一步分析3種分子與分子篩表面Na+的詳細作用情況,對優化后的結構采用Forcite模塊分析工具進行分析。徑向分布函數分析法是分析體系微觀結構和相互作用的常用方法。徑向分布函數g(r)的數學定義是以某一分子或原子為中心,在半徑為r到r+dr的球殼中找到其他粒子的概率[15],其表達式為:

(4)

式中:δr為球殼厚度,nm;r為客體原子與中心原子之間的距離,nm;n(r)為球殼中的客體粒子數。對293 K下CH4,CO2,H2S分子與分子篩中Na+的徑向分布函數進行了分析,結果如圖11所示。由圖11可以看出:H2S與CO2中心原子與Na+的徑向分布圖在r=0.4 nm處均出現了較強的峰,說明兩種分子的中心均與分子篩骨架外的Na+形成了相互作用,其中H2S的主峰的對應半徑是最小的,這是因為H2S分子結構為折線形,S原子容易和孔道中的Na+接觸并相互作用,空間位阻較小;CH4由于結構為正四面體,C原子空間位阻較大,不易與孔道中的活性位點相互作用,因此在作用范圍內幾乎不存在尖峰。除此之外,由3種分子的端位原子徑向分布函數可以看出,H2S和CH4在r=1 nm的范圍內幾乎不存在峰,說明兩者的端位原子與活性位點沒有較強的相互作用。從CO2的徑向分布函數可以看出,中心原子C和端位原子O都與骨架中的Na+產生相互作用,相對于C—Na的徑向分布,O—Na的徑向分布具有更近的主峰半徑和更高的概率值,這說明相對于C原子,CO2中的O原子與分子篩中的活性位點有更強的相互作用,這也解釋了CO2擴散系數最小的原因。

—C; —H

3 結 論

(1)在溫度273,293,313 K下,CH4,CO2,H2S分子在NaX分子篩中的吸附量隨著壓力的增加而迅速增大,并趨于飽和吸附量。在一定溫度范圍內,溫度的升高降低了3種分子的飽和吸附量。相同溫度下吸附作用強度和對應飽和吸附量由大到小的順序為H2S>CO2>CH4。

(2)在溫度293 K下,混合氣體中CO2含量提高,CO2的吸附量增加;隨著混合氣體中H2S含量提高,3種分子的飽和吸附量基本保持穩定,但H2S吸附量遠高于CO2。NaX分子篩對混合氣體中3種分子的吸附選擇性由大到小的順序為H2S>CO2>CH4,其表現出較強的選擇性吸附脫硫性能。

(3)CH4,CO2,H2S分子在NaX分子篩內的擴散系數由大到小的順序為CH4>H2S>CO2。徑向分布函數分析結果表明,CO2的中心原子和端位原子均與孔道活性位產生較強相互作用,解釋了CO2的強擴散阻力現象。

猜你喜歡
擴散系數模型
一半模型
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
一類具有變擴散系數的非局部反應-擴散方程解的爆破分析
3D打印中的模型分割與打包
基于Sauer-Freise 方法的Co- Mn 體系fcc 相互擴散系數的研究
上海金屬(2015年5期)2015-11-29 01:13:59
FCC Ni-Cu 及Ni-Mn 合金互擴散系數測定
上海金屬(2015年6期)2015-11-29 01:09:09
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
非時齊擴散模型中擴散系數的局部估計
Ni-Te 系統的擴散激活能和擴散系數研究
上海金屬(2013年4期)2013-12-20 07:57:07
主站蜘蛛池模板: 国产AV无码专区亚洲精品网站| 国产成人精品男人的天堂| 九九免费观看全部免费视频| 国产成人精品亚洲77美色| 国产a网站| 97国产在线观看| 亚洲精品国产日韩无码AV永久免费网| 精品无码人妻一区二区| 日韩成人在线视频| 亚洲成人精品在线| 国产在线视频自拍| 亚洲国产精品不卡在线| 久久久久久久久亚洲精品| 午夜电影在线观看国产1区| 老司机精品99在线播放| 国产人在线成免费视频| 伊人成人在线视频| av午夜福利一片免费看| 91青青视频| 国产黄色视频综合| 久久www视频| 特级精品毛片免费观看| 久久国产拍爱| 97se亚洲综合在线韩国专区福利| 免费一级α片在线观看| 日韩在线影院| 亚洲色图欧美激情| 亚洲熟女偷拍| 免费国产无遮挡又黄又爽| 91精品伊人久久大香线蕉| 日本午夜网站| 国产第一色| 视频在线观看一区二区| 国产精品一区在线麻豆| 久久国产精品无码hdav| 国产男人天堂| 国产一区二区三区视频| 国产精品毛片一区视频播| 国产最新无码专区在线| 三上悠亚一区二区| 国产午夜精品一区二区三区软件| 高清无码一本到东京热| 国产爽爽视频| 成人午夜视频网站| 99视频在线看| 精品综合久久久久久97超人| 国产在线第二页| 国产欧美日韩精品第二区| 妇女自拍偷自拍亚洲精品| 亚洲天堂.com| 亚洲男人天堂2020| a在线亚洲男人的天堂试看| 国产亚洲精品97在线观看| 国产精品视频系列专区| 亚洲人成网站在线播放2019| 久久毛片网| 麻豆精品视频在线原创| 亚洲av无码人妻| 欧美国产日韩在线观看| 777国产精品永久免费观看| 精品视频91| AV无码无在线观看免费| 中文字幕日韩久久综合影院| 网友自拍视频精品区| 欧美不卡视频在线| 欧美精品亚洲精品日韩专区va| 91亚洲精选| 综1合AV在线播放| 天堂成人av| 久久久久88色偷偷| 一级高清毛片免费a级高清毛片| 欧美丝袜高跟鞋一区二区| 色婷婷电影网| 国产精品一区在线麻豆| 国产黑丝一区| 亚洲视频在线观看免费视频| 最新国产你懂的在线网址| 新SSS无码手机在线观看| 日韩精品一区二区三区swag| 97视频精品全国免费观看| 国产JIZzJIzz视频全部免费| 97视频免费在线观看|