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

催化裂化管式待生劑分配器顆粒分配特性及內部流動特性的CPFD模擬

2020-05-27 01:44:46楊賽飛張永民付明義楊智君
石油學報(石油加工) 2020年3期
關鍵詞:風速分配

楊賽飛, 張永民, 付明義, 楊智君

(1.中國石油大學 重質油國家重點實驗室,北京 102249;2.山東京博石油化工有限公司,山東 濱州 256500)

再生器作為煉油廠流化催化裂化(Fluid catalytic cracking,簡稱FCC)裝置的重要組成部分,其主要作用是利用空氣燒掉裂化反應后附著在催化劑上的焦炭,從而恢復其活性[1]。工業中要求再生后催化劑的炭含量盡可能低,且消耗盡可能少的空氣,以降低裝置能耗。然而,在傳統設計的再生器中,很多裝置沒有設計分配待生劑的結構,待生劑通過待生循環斜管直接輸送進入再生器,導致待生劑在再生器橫截面上分布極不均勻,因而難以實現待生劑上焦炭與主風中氧氣的均勻接觸,不僅造成再生燒焦強度的降低以及裝置整體能耗的增大,更甚者還可能造成稀相尾燃,從而危及裝置安全平穩操作。目前,隨著國內FCC裝置的大型化發展,再生器的直徑不斷擴大,待生劑在再生器內分布不均的問題日益突出。

近年來,為了提高再生器的再生效率,國內外研究者逐漸認識到待生劑分配器的重要作用,開始提出和應用不同的新型待生劑分配器結構。Chen和Patel[2]提出了一種頂部帶有水平分配口的待生劑提升管,待生劑從該提升管底部引入,在提升風的作用下,從頂部水平分配口流入再生器床層頂部,提升管頂部出口結構有利于待生劑從多個水平分配口均勻流出。然而,這種結構需要消耗的提升風量較高,部分提升空氣還有可能直接走短路離開密相床層,容易加劇再生器尾燃的風險,且容許設置的水平分配口有限,待生劑的分配均勻性還有待提高。Myers等[3]提出了一種類似樹枝狀氣體分布器的待生劑分配器結構,待生劑首先進入一個小型流化床,之后再流入帶有主管和支管結構、略向下傾斜的樹枝狀待生劑分配管,分配管上設置有很多催化劑流出孔,分配管內部設置有流化催化劑顆粒的氣體分布管,其設計思想是利用氣體分布管流入流化風的幫助,使催化劑實現較好的流化狀態,從而能像液體一樣從眾多顆粒流出孔中均勻地流入再生器床層頂部。但是,這種氣體分布管設計很難實現所有待生劑分配管內顆粒的流化,一旦部分催化劑沒有很好地流化,則會形成阻礙其他催化劑顆粒繼續流動的障礙,因此很難實現待生劑顆粒整體分配的均勻性。

Benjamin等[4]提出了一種和Chen和Patel[2]類似原理的待生劑分配器,待生劑顆粒首先被引入一根水平輸送管,輸送管一端通入輸送風使得顆粒以稀相氣力輸送的方式流到垂直向下的下行輸送管,下行管底部設置有多個噴嘴,輸送風和待生劑顆粒一起從噴嘴出口噴射進入再生器床層。張振千等[5]提出了一種水平主管和傾斜向下支管的管式待生劑分配器,水平主管中設置有改善催化劑流動性的流化風分配管,在流化風的作用下,待生劑流入傾斜向下的支管中,在重力作用下通過設置在支管上顆粒分配口流入再生器。該結構和Myers等[3]提出的結構類似,在主管中存在局部顆粒失流化阻礙其他顆粒流動的問題,支管中由于沒有流化風,在重力作用下很難實現顆粒沿所有分配口均勻分配。

張永民等[6]根據空氣輸送斜槽原理開發了一種新型槽式待生劑分配器。該分配器是由上部物料輸送槽和下部通風槽構成,顆粒分配口設置在上部槽體兩側,槽體之間用多孔分布板隔開。在流化風作用下,待生劑進入分配器后在上部輸送槽內實現良好的流化狀態,從而可以實現從多個側方分配口均勻流入再生器內。張永民等[7]利用一套可實現催化劑連續循環操作的大型冷模實驗裝置檢驗了依據現有工業設計的管式分配器和船型分配器時的顆粒分配性能,得出管式分配器的性能優于船型分配器。新型槽式待生劑分配器的冷模實驗研究[8]表明,該結構可以在較小流化風量下實現待生劑的均勻分配,分配均勻性顯著優于傳統的船型和管式分配器,且具有較大的操作彈性和較小的流動阻力,不會影響待生輸送管路的顆粒輸送能力。但是,該冷模實驗的缺點是只能測量各出口催化劑顆粒的流量和分配比,無法對分配器內氣-固流動特性進行詳細觀察和準確測量,也無法準確測量各出口氣體流量的大小和分布。

近年來,隨著計算機技術及顆粒流體動力學的快速發展,Snider[9]提出了一種基于歐拉-拉格朗日方法的計算顆粒流體動力學(Computational particle fluid dynamics,簡稱CPFD)新方法。CPFD基于多相質點網格法(Multiphase particle-in-cell method,簡稱MP-PIC)[10-11],對具有相同密度、體積、速度及位置的顆粒進行打包處理,同時采用正應力梯度來簡化計算顆粒間的相互作用,從而使得顆粒數目達數十億的工業級流化床歐拉-拉格朗日模擬成為可能。

筆者在前期實驗研究基礎上,運用CPFD多相流模型對工業裝置中較為普遍應用的一種待生劑分配器型式——管式待生劑分配器進行系統的數值模擬研究。管式待生劑分配器不僅常用于催化裂化裝置中待生劑分配,也常用于外取熱器中冷催化劑在再生器中的分配。本研究的模擬對象是前期冷模實驗中[7]采用的管式待生劑分配器模型,首先模擬輸送風和松動風對分配均勻性的影響規律,利用已有實驗數據檢驗所用CPFD模型的可靠性和準確性,然后進一步考察分配器傾斜角度對顆粒分配均勻性的影響。在模擬過程中,同時還考察了分配器內部流動特性以及出口氣體流量分布特性的變化規律。

1 數學模型

1.1 控制方程

在CPFD數學模型中,氣相被視作連續相,采用歐拉法進行計算,而顆粒相被視作離散相,采用拉格朗日法求解。其控制方程如下:

氣相連續方程:

(1)

氣相動量方程:

(2)

單位體積相間動量交換速度:

(3)

反映顆粒分布函數f隨時間變化的Liouville方程:

(4)

顆粒加速度方程:

(5)

顆粒正應力方程:

(6)

顆粒體積分數方程:

(7)

1.2 曳力模型

在氣-固流化床模擬中,氣-固相間的相互作用主要由曳力模型來描述,因此選擇合適的曳力模型對于模擬結果的準確性具有重要的影響[12-13]。較為常用的曳力模型有Wen-Yu模型、Ergun模型、Turton&Levenspiel模型及Gidaspow模型等。筆者選用了經典的Gidaspow曳力模型[14]來模擬。該模型結合了適用于稀相系統的Wen-Yu曳力模型[15]和適用于密相系統的Ergun曳力模型[11],其施加在顆粒上的曳力可表示為:

(8)

本研究對應的前期實驗[7]中顆粒相最大堆積體積分數為0.63(堆積密度比真實密度),且考慮到本研究體系有稀相到密相的過渡區域,因此選用修正的Gidaspow曳力模型[16-18],即Wen-Yu方程和Ergun方程組成的混合模型。

當顆粒體積分數小于0.75θcp時,采用Wen-Yu方程計算曳力:

(9)

(10)

當顆粒體積分數大于0.85θcp時,采用Ergun方程計算曳力:

(11)

為避免稀相到密相的過渡區域中由于曳力系數不連續引起的數值問題,當顆粒體積分數在0.75θcp和0.85θcp之間時,采用過渡方程計算曳力:

(12)

其中Re為雷諾數,用下式計算:

(13)

rp為顆粒半徑,表達式如下:

(14)

2 模擬設置

2.1 幾何模型與網格劃分

本模擬所采用的幾何模型是張永民等[7]前期采用的一套評價待生劑分配器性能的大型冷模實驗裝置中所用的模型。由于本研究重點關注管式分配器內的氣-固流動特性及顆粒分配均勻性,這里對模擬對象進行了簡化,即只模擬管式分配器部分,如圖1所示,省去了其上方的流化床、下部的顆粒集料倉和提升管等部分。實驗中的管式分配器是一根內徑為90 mm,厚為5 mm的有機玻璃管,其入口端與待生劑斜管的水平管段相連,末端使用蓋板封閉。在水平管段左側同心安裝了1個直徑為26 mm的輸送風管。同時,在分配器兩側分別交錯地開設了10個直徑為26 mm、等間距布置的排料口。此外,為了松動和流化分配器底部的顆粒,在其下方排布了19個直徑為6 mm的短管來通松動風。

圖1 前期實驗中采用的管式待生劑分配器[7]
Fig.1 Pipe spent catalyst distributor used in the previous experiments[7]

在模擬中,管式分配器的幾何模型如圖2所示,其幾何結構與實驗裝置基本一致。為了便于直觀地觀察各個出料口顆粒的流出情況,在管式分配器的外圍設置了1個方形計算域(2100 mm×160 mm×170 mm)。管式分配器的長度和直徑分別為 1900 mm 和90 mm。由于在網格劃分時分配器兩側的10個圓形排料口難以捕捉,因此將其簡化為10個方形出口(26 mm×26 mm)。Barracuda軟件采用笛卡爾網格切割單元法對幾何模型進行網格劃分。經網格無關性驗證之后,管式分配器的網格劃分如圖3所示,網格總數為44982。

圖2 模擬中管式待生劑分配器幾何模型的主視圖和左視圖Fig.2 Main view and left view of the geometric model of the pipe spent catalyst distributor in the simulation(a) Main view; (b) Left view

圖3 管式待生劑分配器的網格劃分Fig.3 Meshing of the pipe spent catalyst distributor

2.2 邊界條件及模擬參數設置

在圖2中,管式分配器的左端為顆粒入口邊界和輸送風入口邊界。顆粒選用FCC平衡催化劑,密度為1500 kg/m3,最大堆積體積分數為0.63,索特平均直徑(Sauter mean diameter)為60 μm,質量流率為1.2 kg/s。輸送風選用的氣體壓力和溫度分別為101325 Pa和300 K,密度為1.173 kg/m3,風速為1~6 m/s。在分配器的底部均勻設置了19個射流點源來模擬松動風入口,其位置排布與實驗完全一致。采用點源注射的進氣方式可以規避因捕捉幾何小孔所導致的局部網格過密的問題,進而優化了網格質量和提高了計算速率。松動風選用的也是常溫常壓空氣,風速范圍為0~0.55 m/s。分配器外圍方形計算域的底部設置為壓力邊界,只有氣體可以排出。為了研究分配器傾斜角對顆粒分配均勻性的影響,模擬中采用分解重力加速度的方法來實現5個不同的分配器傾斜角度(進料端固定,末端向下旋轉),分別為0°、2°、6°、10°和14°。初始狀態下,整個計算域只有空氣,沒有顆粒,壓力和溫度分別為101325 Pa和300 K。在分配器前后兩側的10個排料口分別設置流量監測面,以此來實時監測由各個排料口流出的顆粒和氣體的流量。模擬的時間是20 s,初始模擬時間步長的設置要根據模擬對象先給定,本研究中初始時間步長均設置為0.002 s,其后Barracuda會根據每步計算所生成的CFL值自動調整計算時間步長至合適的大小。

3 模擬結果與討論

3.1 輸送風和松動風對顆粒分配均勻性的影響

為了定量評價待生劑分配器的顆粒分配性能,張永民等[7]在前期實驗中提出了一個顆粒分配不均勻指數ξ,其表達式如下所示:

(15)

(16)

其中,ci為各個出料口顆粒實際流出量wi與理想流出量W/n之比。ξ反映了分配器實際分配與理想分配之間的差異性。ξ越大,表明分配器內顆粒分配均勻性越差;ξ越接近于零,表明顆粒橫向分配的均勻性越好。

圖4 顆粒分配不均勻指數(ξ)隨輸送風速(ut)和松動風速(ua)的變化Fig.4 Variation of the heterogeneity index (ξ) of particle distribution under different flowrates of transport air (ut) and fluidizing gas (ua)(a) Experimental result; (b) Simulated result

根據穩定狀態下(模擬開始15 s后)分配器的10個出料口監測的瞬時顆粒質量流量,可以計算得到各個出料口的平均顆粒質量流量,進而運用式(15)和式(16)得到不同操作條件下的顆粒分配不均勻指數ξ。圖4(a)和4(b)分別給出了實驗和模擬的ξ隨輸送風速和松動風速的變化情況。由圖4(b)可以看出,隨著輸送風速的增大,顆粒分配不均勻指數明顯減小,增加松動風速能有效改善顆粒橫向分配的均勻性,但是進一步增大松動風速后,顆粒分配不均勻指數并無顯著變化。出現這種現象的原因是:隨著輸送風速的增大,顆粒獲得了更大的橫向動能,進而更容易流向分配器的末端。此外,當在分配器底部注入松動風后,顆粒被流化起來,此時在輸送風的作用下,顆粒更容易從分配器末端的出口流出,從而進一步改善了顆粒分配的均勻性。這與前期實驗中得到的結論基本一致,即增加輸送風和松動風均能改善顆粒的橫向分配均勻性,但總體上,松動風的影響不如輸送風顯著。通過對比圖4(a)和圖4(b)可以發現,實驗結果中不均勻指數ξ在輸送風ut>3 m/s時發生陡降,而模擬結果并沒有出現。這主要是因為,模擬中省去了分配器前端的水平管段以及與之相連的下料立管,而水平管段與下料立管內的氣-固流動形態會對分配器內的氣-固流動產生直接的影響,進而影響顆粒分配的均勻性。這在前期實驗[7]中有詳細的解釋,所以從趨勢上看,此模擬結果是可靠的。

3.2 傾斜角度對顆粒分配均勻性的影響

在前面的分析中提到,顆粒沿分配器的橫向動能是影響顆粒分配均勻性的重要因素,橫向動能越大,顆粒越容易向分配器的末端流去。鑒于此,考察了分配器傾斜角度θ對顆粒分配均勻性的影響。圖5為輸送風速1 m/s時,有無松動風對不均勻指數隨分配器傾斜角度變化的影響。與圖5相對應,圖6為輸送風速1 m/s,無松動風時不同傾斜角度下的分配器平均顆粒體積分數云圖。

由圖5可以看出,當無松動風時,隨著傾斜角度的增大,顆粒分配不均勻指數呈現了一個先急后緩的降低趨勢,表明顆粒的分配均勻性越來越好。這可以結合圖6加以解釋,隨著傾斜角度的逐漸增大,顆粒在軸向上的重力分量也不斷增大,從而獲得更大的橫向動能,使得更多囤積在分配器前端的顆粒向末端流去。然而,分配器的傾角也不是越大越好,當傾角過大時,大部分顆粒會從靠近末端的出口排出,不利于顆粒的均勻分配。這也是不均勻指數降低程度逐漸放緩的原因所在。當有松動風時,隨著傾斜角度的增大,不均勻指數呈先下降后上升的變化趨勢,并且在θ=6° 時達到最小,表明顆粒的分配均勻性在該角度時達到最好。造成這種結果的原因與上述無松動風時類似,傾角的增大使得顆粒獲得更大的軸向推動力,進而使得顆粒更容易流向分配器的末端;然而,過大的傾角會導致大部分顆粒涌向分配器的末端,導致顆粒分配均勻性變差。總體來看,分配器傾斜角度對顆粒的分配均勻性具有顯著的影響,當分配器的傾斜角度在一定范圍內時,有松動風的顆粒分配均勻性明顯好于無松動風。

圖5 有無松動風時不均勻指數(ξ)隨分配器傾斜角度(θ)的變化Fig.5 Change of the heterogeneity index (ξ) under different distributor tilt angles (θ)ut=1 m/s

圖6 不同傾斜角度下無松動風的分配器內平均顆粒體積分數云圖Fig.6 Contours of solids fraction in the distributor at different tilt anglesut=1 m/s; ua=0 m/s

3.3 輸送風速和松動風速對各出口顆粒流量的影響

圖7為分配器各出口平均顆粒流量隨輸送風速和松動風速的變化。為了便于區分,將分配器自入口端至末端的10個出料口依次編號為1,2,……10。由圖7可以看出,在任意松動風速和輸送風速下,均有一個出料口顆粒流量最大,并且隨著輸送風速的增大,最大出料口的位置會逐漸向分配器的末端移動,同時該出料口的顆粒流量也逐漸減小。此外,與不加松動風(見圖7(a))相比,加松動風后分配器各出口的顆粒流量更加接近,顆粒分配更加均勻,并且隨著松動風速的繼續增大,顆粒分配的均勻性變化不大。出現上述現象主要是因為,輸送風速的增大使得分配器入口端更多的顆粒被氣體夾帶并朝著分配器的末端移動,從而使得分配器前密后稀的顆粒分布狀況得到改善。與此同時,當輸送風速固定時,增加松動風速可以使得分配器內部尤其是底部區域的顆粒流化起來,增大了顆粒之間的空隙,降低了氣體流動阻力,使得更多的顆粒更容易被輸送風夾帶移動,從而在減少流化死區的同時也使得分配器內的顆粒流動更加順暢,最終導致各出料口的顆粒流量更加接近。然而,由于松動風主要是對顆粒的縱向運動產生影響,因此進一步增大松動風速對顆粒在分配器內的橫向運動影響甚微。

3.4 輸送風和松動風對各出口風量的影響

分配器內氣體的運動帶動了顆粒的流動及排出,反過來,顆粒的分布形態也影響著氣體的運動軌跡,進而影響各出口的排風量。圖8為有無松動風時各出口排風量占總排風量的百分數隨輸送風速的變化。由圖8可以看出,無論有無松動風,分配器各出口的排風量占總排風量的百分數基本呈現一個中間低兩邊高的“碗狀”結構,并且隨著輸送風速的增大,分配器中間區域的出口排風量占總排風量的百分數逐漸增大,而靠近分配器末端的出口排風量占總排風量的百分數逐漸減小。當輸送風速固定時,與無松動風(見圖8(a))相比,加松動風(見圖8(b))后,靠近分配器前端的出口排風量占總排風量的百分數略有降低,而靠近末端的出口排風量占總排風量的百分數明顯增大。

圖7 輸送風速(ut)和松動風速(ua)對出口顆粒流量分布的影響Fig.7 Effects of the transport air (ut) and the fluidizing gas (ua) on the particle flux distribution at the 10 outletsua/(m·s-1): (a) 0; (b) 0.11; (c) 0.22; (d) 0.55

圖8 輸送風和松動風對出口氣體流量分布的影響Fig.8 Effects of transport air and fluidizing gas on the air flow distribution at the 10 outletsua/(m·s-1): (a) 0; (b) 0.22

上述現象可以結合管式分配器有無松動風時的平均顆粒體積分數云圖得到很好的解釋,如圖9所示。由圖9可見,無論有無松動風,顆粒均主要集中在分配器的中間區域,導致氣體優先從分配器的前端和末端的出口排出,進而形成了圖8中所示的“碗狀”分布。隨著輸送風速的增大,當無松動風時,更多中間區域的顆粒被夾帶到分配器的末端,使得更多的氣體可以從中間的出口排出;當有松動風時,由于分配器內的顆粒被流化起來,顆粒間空隙率增大,更多的氣體可以穿過顆粒間隙從中間的出口排出,因而兩者均導致了分配器排風量占總排風量的百分數在中間出口增大而在末端出口減小的趨勢。通過對比圖9(a)和圖9(b)可以看出,加松動風后,分配器內的顆粒在縱向上的高度明顯降低,使得輸送風更容易從顆粒上方的無阻力區流向分配器末端的出口,從而導致靠近末端的出口排風量占總排風量的百分數顯著提高。

圖9 管式分配器有無松動風時的平均顆粒體積分數云圖Fig.9 Comparison of solids fraction contours in the pipe distributor with and without fluidizing gas(a) ua=0 m/s; (b) ua=0.22 m/s; ut=1 m/s

4 結 論

(1) 采用CPFD方法模擬得到的顆粒分配不均勻指數隨輸送風速和松動風速的變化趨勢和數值上總體與實驗結果相吻合,從而驗證了CPFD模型模擬待生劑分配器氣-固流動特性的可行性。同時進一步說明了增加輸送風和松動風能改善顆粒分配的均勻性,但松動風的改善效果不如輸送風顯著。

(2) 待生劑分配器傾斜角度對顆粒分配的均勻性具有顯著的影響,增大傾角能有效降低分配不均勻指數,但傾角過大時,顆粒的分配效果反而變差。在本研究中,當傾角θ=6°時,顆粒的分配均勻性達到最佳。鑒于此,可以考慮改變目前工業裝置中常采用的水平管式分布器結構,使其傾斜一定的角度以達到更好的顆粒分布均勻性。

(3) 隨著輸送風速的增大,出料流率占比最大的出口逐漸遠離分配器入口端,并且該占比也逐漸變小。與不加松動風相比,加松動風后分配器各出口的顆粒流量更加接近,顆粒的分配均勻性更好。

(4) 分配器各出口的氣體流量分布呈“碗狀”結構,即兩端高中間低的趨勢,這與分布器內顆粒料層的高度分布有關。輸送風速的增大有助于縮小各出口氣體流量的差距。

符號說明:

英文符號:

ci——某個出料口顆粒實際流出量與理想流出量之比;

Cd——曳力系數;

C1——線性系數;

C2——非線性系數;

dp——顆粒半徑,m;

D——氣-固相間曳力系數;

D1——顆粒體積分數小于0.75θcp時,采用Wen-Yu方程計算曳力,N;

D2——顆粒體積分數大于0.85θcp時,采用Ergun方程計算曳力,N;

Dp——顆粒相間曳力系數;

f——顆粒概率分布函數;

i——出料口編號;

mp——顆粒質量,kg;

n——出料口個數;

P——氣相壓力,Pa;

PS——常數,Pa;

rp——顆粒半徑,m;

Re——雷諾數;

t——時間,s;

ua——松動風大小,m/s;

ut——輸送風大小,m/s;

Vp——顆粒體積,m3;

W——顆粒流出總量,kg;

wi——某個出料口顆粒實際流出量,kg;

希臘符號:

β——常數,被設置為默認值3;

ε——一個數量級為10-7的常數;

τg——非靜力應力,N/m2;

τp——顆粒正應力,N/m2;

θ——分配器傾斜角度,°;

θcp——顆粒最大體積分數;

θg——氣體體積分數;

θp——顆粒體積分數;

ρg——氣體密度,kg/m3;

ρp——顆粒密度,kg/m3;

μg——氣體動力黏度,Pa/s;

ξ——顆粒分配不均勻指數;

下標:

p——顆粒相;

g——氣體相。

猜你喜歡
風速分配
基于可行方向法的水下機器人推力分配
基于Kmeans-VMD-LSTM的短期風速預測
基于最優TS評分和頻率匹配的江蘇近海風速訂正
海洋通報(2020年5期)2021-01-14 09:26:54
應答器THR和TFFR分配及SIL等級探討
遺產的分配
一種分配十分不均的財富
績效考核分配的實踐與思考
基于GARCH的短時風速預測方法
考慮風切和塔影效應的風力機風速模型
電測與儀表(2015年8期)2015-04-09 11:50:06
GE在中國發布2.3-116低風速智能風機
主站蜘蛛池模板: 国产一级毛片在线| 欧美国产在线一区| 亚洲精品高清视频| 国产精品尤物铁牛tv| 亚洲综合二区| 高清视频一区| 国产原创第一页在线观看| 综1合AV在线播放| 视频二区国产精品职场同事| 色噜噜狠狠狠综合曰曰曰| 啪啪国产视频| 亚洲 成人国产| 亚洲一区二区三区中文字幕5566| 国产成人精品在线1区| 国内精品九九久久久精品| 欧美天堂在线| 国产在线精品99一区不卡| 夜夜拍夜夜爽| 日韩精品成人网页视频在线 | 九九视频在线免费观看| 国产成人高清精品免费| 成年A级毛片| 91久久天天躁狠狠躁夜夜| 乱人伦99久久| 亚洲综合专区| 国产真实乱子伦精品视手机观看 | 欧美成人免费| 无码aaa视频| 黄网站欧美内射| 午夜日b视频| 99久久精品免费看国产免费软件 | 91精品啪在线观看国产| 免费观看国产小粉嫩喷水| 欧美福利在线| 精品人妻无码区在线视频| 婷婷六月综合| 97青草最新免费精品视频| 亚洲最新地址| 日本三区视频| 国产主播喷水| 久久性视频| 亚洲91精品视频| 亚洲精品视频免费| 亚洲国产精品久久久久秋霞影院| 欧美一级黄色影院| www成人国产在线观看网站| 最近最新中文字幕在线第一页 | 亚洲v日韩v欧美在线观看| 亚洲成人网在线观看| 直接黄91麻豆网站| 国产三级成人| 久久综合九九亚洲一区| 国产欧美日韩精品综合在线| 在线观看欧美精品二区| 国产成人精品日本亚洲77美色| 18禁高潮出水呻吟娇喘蜜芽| 区国产精品搜索视频| 四虎永久在线精品影院| 欧美成人区| 亚洲欧美一区在线| 国产美女精品人人做人人爽| 亚洲天堂网在线播放| 乱系列中文字幕在线视频| 成人韩免费网站| 久久视精品| 亚洲an第二区国产精品| 在线不卡免费视频| 国产乱人免费视频| 欧美午夜在线观看| 国产精品白浆在线播放| 99久久精品国产自免费| 五月婷婷亚洲综合| 欧美国产综合色视频| 91青青在线视频| 亚洲欧美精品在线| 香蕉伊思人视频| 亚洲欧美一区二区三区麻豆| 国产中文一区二区苍井空| 日韩免费毛片视频| 在线观看欧美国产| 色香蕉影院| 色婷婷亚洲综合五月|