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

流化床內松木屑氣化過程的MP-PIC數值模擬

2022-11-16 13:33:44萬章豪楊世亮
石油學報(石油加工) 2022年6期

萬章豪, 楊世亮, 王 華

(昆明理工大學 省部共建復雜有色金屬資源清潔利用國家重點實驗室,云南 昆明 650093)

在化石資源日益枯竭及二氧化碳過度排放等全球能源和環境危機背景下,開發和利用第四大能源載體——生物質是實現“2030碳達峰、2060碳中和”的國家戰略目標的重要步驟[1]。生物質作為來源廣、產量大的可再生能源,現已逐漸替代化石燃料并占據全球可再生能源供應量的14%[2-3]。流化床因具有良好的氣-固混合、傳熱和傳質效果[4],已廣泛應用于生物質氣化等熱化學轉化工藝。為提升生物質氣化效率、優化氣化工藝,國內外學者對氣化操作條件進行了大量實驗研究。然而,由于密相流化床內的高溫及不可視特性,實驗中對其內部氣-固兩相運動規律和微觀多相反應過程的探索難以進行。

隨著數學算法的發展和計算機性能的提升,計算流體力學(Computational fluid dynamics, CFD)可以彌補實驗速度慢、成本高、檢測難等缺點,已成為研究流化床中多物理場(流場、溫度場、組分場)變化規律的強有力手段[5]。目前,用于稠密氣-固兩相流動的數值模擬方法主要分為Euler-Euler法和Euler-Lagrange法。Euler-Euler法將氣相和顆粒相均視為連續介質[6-7],且假定固相具有均一的顆粒尺寸和密度,其可以用較低的計算成本精準預測氣-固介質流動的宏觀特性[8-9],但卻無法實現多尺度顆粒的計算。Euler-Lagrange法將氣相和顆粒分別視為連續相和離散相,其從顆粒尺度描述固相運動。目前,基于Euler-Lagrange的計算方法主要有計算流體方法耦合離散元模型(Computational fluid dynamics-discrete element method, CFD-DEM)和多相質點網格法(Multiphase particle-in-cell, MP-PIC)。其中,CFD-DEM模型可準確追蹤單個顆粒的軌跡,并采用軟球模型對粒子間的相互碰撞進行高精度模擬[10-11]。然而,當前的計算機水平難以求解千萬量級的顆粒數量,所以CFD-DEM模型對于模擬工業尺度甚至實驗室尺度的流化床有較大的局限性[12-13]。相比于CFD-DEM模型,MP-PIC模型因精度相當[14-16]、計算成本更小等特點被用于預測實驗室尺度和工業尺度流化床中顆粒的流動特性及生物質氣化過程中熱質傳遞過程。

近幾年,國內外眾多學者基于MP-PIC模型研究了氣化反應影響因素,優化了氣化工藝,提高了氣化效率。流化床中床料不僅作為熱載體為生物質氣化提供熱量,而且還增加了生物質在床內的停留時間[17],當生物質顆粒進入鼓泡床后迅速氣化,此時生物質熱解反應和氧氣的還原反應同時發生于鼓泡床底部[18]。增加鼓泡床內空氣當量比(ER)、蒸汽生物質比(SBR)、氣化溫度和生物質粒徑寬度均可提升碳轉化率[17,19]。通入蒸汽可以提升雙循環床中顆粒的循環效率和H2的產率,但卻降低了生物質氣化的穩定性[20-21],而溫度的提升僅增加了CO的濃度[22]。

然而,目前對于流化床內生物質氣化過程中復雜氣-固湍動、相間熱質傳遞和氣化反應的強非線性耦合機理的了解遠遠不足,且對流化床氣化反應過程中氣相熱物性組成及其分布規律的研究有限。為了加深對流化床內生物質氣化過程中固相流體動力學特性及氣-固熱質傳遞規律的理解。筆者在Euler-Lagrange框架下,借助OpenFOAM開源軟件并耦合氣體-顆粒相互作用、傳熱和化學反應模塊,建立用于模擬三維流化床氣化爐內生物質氣化過程的MP-PIC模型。基于可靠的數值驗證結果,預測了流化床內顆粒的空間分布規律,探尋固相的流動和傳熱特性,分析不同空氣當量比、蒸汽生物質比、氣化溫度對爐內氣體組分和氣相熱物性的影響,以期能為推動綠色能源可持續發展、提升生物質氣化技術和生物質流化床氣化反應器的設計、放大以及優化提供強有力的理論支撐。

1 數學模型

1.1 氣-固兩相控制方程

MP-PIC模型中,氣相在歐拉框架下被視為連續相并基于Navier-Stokes方程求解,而氣相湍流采用大渦模擬(Large-eddy simulation,LES)來計算[23]。顆粒則被視為離散相,在拉格朗日框架下基于多相網格質點法求解。氣相的連續性、動量、能量方程如式(1)~式(3)所示。

(1)

(2)

(3)

氣相的組分方程如式(4)所示。

(4)

(5)

式中:μt為氣相湍流黏度,kg/(m·s);筆者選擇Schmidt數SC為0.9。

式(3)中q表達式可列為:

(6)

式中:λg代表氣相的熱傳導率,J/(s·m·K);Tg為氣體溫度,K。

(7)

Sgw=hgwAgw(Tw-Tg)

(8)

式中:hi為氣相組分i的生成熱,J/kg;Agw表示氣體與壁面之間的局部接觸面積,m2;hgw表示氣體與壁面之間的對流熱交換系數,W/(m2·K);Tw為壁面溫度,K。

為平衡計算成本和計算效率,MP-PIC模型將物理屬性相同(如直徑、密度、溫度等)的實際顆粒群打包成為一個計算顆粒[24],并采用顆粒分布函數fs(Particle distribution function,PDF)和固相應力分別表征固相運動和顆粒碰撞[25]。

(9)

式中:fs為顆粒分布函數;us為顆粒速度矢量,m/s;fD為局部質量平均的顆粒分布函數;τD為顆粒碰撞的阻尼時間,s;x表示為x軸方向的位移,m。考慮到空氣阻力、壓力梯度、碰撞和重力,式(9)中的顆粒加速度Γs(m2/s)可表示為:

(10)

(11)

(12)

式中:CD表示曳力系數;εg表示空隙率;us為顆粒速度矢量,m/s;ds為顆粒直徑,m;Res為雷諾數。

1.2 反應動力學模型

生物質氣化過程主要包含干燥、熱解和均相/非均相反應,其化學反應動力學模型基于單速率的Arrhenius方程進行描述,速率常數Ri方程如式(13)所示。

(13)

式中:mi、Ai和Eai分別為組分i的質量(kg)、指前因子和活化能(kJ/mol)。

圖1為顆粒流化和生物質氣化示意圖。生物質由給料口進入流化床氣化爐后,在氣相作用力和顆粒碰撞力的共同作用下做無規則運動,并在爐內高溫氣氛熱傳遞的作用下逐漸升溫。當生物質達到臨界溫度后開始脫水和熱解。其中,生物質熱解階段釋放的揮發分包含碳氫、碳氧化合物(如CH4、C2H4、CO、CO2)及少量H2。空氣從爐底引入后在高溫條件下分別與揮發分中的可燃氣體和生物質碳轉化產生的焦炭發生均相和非均相反應。生物質中水分的蒸發、顆粒熱解和均相與非均相反應幾乎同時發生,氣-固兩相在復雜的化學反應過程中的傳質傳熱行為劇烈且迅速。筆者所采用的生物質原料為松木屑,其工業分析和元素分析結果如表1所示[27]。

圖1 氣化爐床內顆粒流化及生物質氣化示意圖Fig.1 Schematic diagram of particle fluidization and biomass gasification in the gasifier bed

表1 生物質的工業分析和元素分析Table 1 Proximate and ultimate analysis of biomass

計算過程中難以考慮生物質實際氣化過程中的所有化學反應,因此筆者僅考慮包括氧化反應(R3)、Boudouard反應(R4)、蒸汽氣化反應(R5)在內的非均相反應,以及包括CO、H2、CH4、C2H4的氧化反應(R6~R9)、水煤氣置換反應(R10~R11)、蒸汽甲烷重整反應(R12)在內的均相反應。反應方程式[28]及相關速率[28-31]如表2所示。

表2 生物質氣化過程反應方程式及動力學參數Table 2 Reaction equations and kinetic parameters for the biomass gasification process

2 模擬條件及模型驗證

2.1 三維模型和邊界條件

筆者以中國科學院廣州能源研究所建立的小型流化床氣化爐為研究對象[27],氣化爐三維模型尺寸示意圖如圖2所示。反應器總高1400 mm,其中流化區和自由空域的內徑分別為40 mm和60 mm。生物質顆粒和水蒸氣分別從反應器兩側通入,且生物質入口與氣體分布器相距40 mm。為穩定流化狀態、提高傳熱效率,將熱導率較高的30 g硅砂作為床料堆置于氣化爐底部。初始時爐內的N2和床料被加熱至973.15 K。反應開始后,將預熱至338 K的空氣作為流化劑以0.5 m3/h的恒定流量由底部氣體分布器通入反應器中,而被蒸汽發生器預熱至427 K的蒸汽以1.2 kg/h的質量流量從流化區的側面引入。生物質顆粒和硅砂床料的密度分別為556、2300 kg/m3,直徑分別為0.375、0.25 mm。模擬過程中,反應器在三維空間中以均勻的網格離散化,網格數量最終為42757個。時間步長設定為1×10-4s,總模擬時間為20 s。

圖2 氣化爐三維模型示意圖Fig.2 Schematic diagram of the three-dimensionalmodel for gasifier

2.2 模型驗證

為確保數值模型的可行性,筆者對氣化爐出口處的氣體產物進行了定性和定量驗證。數值計算得到的氣體組分的時間平均摩爾分數與實驗結果[27]的對比如圖3所示。在生物質氣化產物中,含量最高的為CO,H2與CO2的含量大致相同,CH4次之,而C2H4的含量最低。數值模擬與實驗的誤差主要源于三維模型的簡化以及部分副反應動力學模型的忽略。但總體而言,模擬結果與實驗結果對比良好,各氣體組分的含量和分布規律基本一致。因此,MP-PIC模型對于氣化爐中生物質氣化過程的模擬具有合理性和可靠性。

x—Mole fraction of gas species圖3 實驗數據與模擬結果的氣體組分摩爾分數對比Fig.3 Comparison of mole fraction of gas speciesbetween simulation and experimental results

3 結果與討論

在氣化過程中,生物質顆粒受空氣、蒸汽等氣化劑作用,在高溫條件下經熱化學反應轉化為可燃氣。適量空氣的引入可以提升爐內焦炭不完全燃燒的反應強度和生物質的氣化強度。ER表示實際參與生物質燃燒的空氣質量與理論生物質完全燃燒需要的空氣質量之比,其表達式如式(14)所示[18]。

(14)

氣化過程中引入一定量的蒸汽會提高生物質氣化效率,這歸因于其促進了蒸汽重整和水煤氣置換反應[32]。然而,通入過量的蒸汽會降低爐內溫度并抑制氣化強度[33]。SBR表示入口處通入的蒸汽質量與干生物質質量之比,其方程式如式(15)所示。

(15)

3.1 顆粒物理性質

在ER為0.22、SBR為2.2、Tw為1123 K時,顆粒在流化床內的軸向運動速率(Usz)及氣-固相對滑移速率(ΔUslip)的三維空間分布如圖4(a)和(b)所示,流化床氣化爐內的顆粒溫度(T)的空間分布如圖4(c)所示。由圖4(a)和(b)可知:密度較大的床料主要集中于反應器下部區域,且具有較大的滑移速率;而較輕的生物質在流化氣體的裹挾下加速上升,滑移速率逐漸減小。此外,蒸汽的引入增強了氣-固的相間動量傳遞,顆粒在距底部空氣分布器71 mm高處具有最大流化速率和滑移速率分別為3.72、0.138 m/s。然而,管徑的擴張降低了氣相的表觀氣速,此時自由空域顆粒的上升速率明顯低于流化區顆粒。由圖4(c)可知,低溫空氣和高度床料的熱交換降低了硅砂的溫度,而可燃物與氧氣在床料層上方劇烈燃燒所釋放的熱量加熱了床內顆粒。因此位于床底部區域的顆粒溫度隨高度持續增加。而吸熱的生物質氣化反應及顆粒-低溫蒸汽間的熱傳遞致使顆粒溫度在蒸汽入口高度有明顯降低,并在自由空域達到溫度的穩定狀態。

圖4 氣化爐內瞬時顆粒分布及物理性質(t=20 s)Fig.4 Instantaneous distribution and physical properties of particles in the gasifier (t=20 s)(a) Usz; (b) ΔUslip; (c) T

3.2 ER對氣化反應的影響

3.2.1 ER對氣化產物的影響

在SBR為2.2、Tw為1123 K條件下,ER對氣體組分沿反應器高度方向的分布如圖5所示。圖5(a)~(f)分別為C2H4、CH4、CO、H2、CO2、H2O等氣體在不同ER條件下的質量分數分布。生物質氣化會產生大量的揮發分氣體,而蒸汽的注入則顯著影響合成氣含量。因此,揮發氣體的含量在反應器下部沿軸向迅速增大,但在0.4 m高度時受蒸汽稀釋而陡然降低。此外,蒸汽的引入會促進蒸汽氣化反應(R5)和水煤氣置換反應(R10),而各組分氣體產量在蒸汽入口上方有輕微增幅并隨后達到穩定狀態。ER從0.1增加至0.15時,可燃物的氧化強度持續增大,可燃氣體的含量不斷減小,而CO2的含量隨著ER的增大而增大。然而, ER的不斷提升減少了生物質的停留時間,使生物質顆粒的氣化過程難以充分進行。因此,ER從0.15增加至0.25時,揮發分氣體含量逐漸減小。ER對不同區域的蒸汽含量影響各異,具體表現為:流化區域的蒸汽含量在不同操作參數下無明顯變化,而自由空域的蒸汽含量分布隨ER的增大而減小。這說明增大氧含量對氣化過程中的蒸汽產率基本無影響,而表觀氣速的提升稀釋了從爐側通入的蒸汽含量。

w—Mass fraction of gas species; z—Axial height圖5 空氣當量比(ER)對氣體組分質量分數軸向分布的影響Fig.5 Effect of air equivalence ratio (ER) on the axial distribution of mass fraction of the gas species(a) C2H4; (b) CH4; (c) CO; (d) H2; (e) CO2; (f) H2OConditions: SBR=2.2; Tw=1123 K

3.2.2 ER對氣體熱物性的影響

當SBR為2.2、Tw為1123 K時,不同ER條件下氣體熱物性沿反應器高度方向的分布如圖6所示。由圖6(a)可知,反應器底部的溫度在4種ER條件下均隨著高度的增加而升高,且在蒸汽入口高度處急劇降低。這是由于可燃物的氧化放熱加熱了爐內氣體,而低溫蒸汽的引入降低了反應器內的溫度。隨ER的增大,高溫氣氛與低溫空氣的熱傳遞效應增強,流化區底部的氣相溫度逐漸降低。另外,ER在0.1~0.2范圍內,氧化強度隨空氣流量的增加而增大。當ER大于0.2時,爐內溫度反而有所降低,這歸咎于過大的氣相曳力帶走了部分生物質顆粒,使得可燃物的氧化強度減弱,且大流量低溫空氣降低了爐內溫度。

氣化爐內合成氣的密度取決于氣相分子活躍度和氣體成分的共同作用,而分子活躍性與溫度密切相關。由圖6(b)可知,混合氣體密度在流化區下部隨ER持續增大。這是因為高溫氣氛促進了分子的運動活性,氣相密度隨之減小。相較于H2、CH4、H2O(g)等揮發產物,N2的密度較大,而ER為0.25時高含量的N2中和了合成氣導致的低密度分布。因此,結合圖5和圖6(b)可知,ER為0.1和0.25時的氣相密度大于ER為0.15和0.2時的混合氣體密度。

圖6(c)表示不同ER條件下比熱容沿軸向的分布曲線。氣相比熱容沿軸向持續增加,在蒸汽入口處的增大速率最為明顯。氣體成分是混合氣體比熱容的決定性因素,其中具有較小比熱容的CO2是致使混合氣體比熱容較低的主要因素,而混合氣體具有較大比熱容則歸功于蒸汽[34]。因此,ER為0.25和0.1時較低/較高的比熱容取決于該條件下較低/較高的蒸汽含量。

圖6 空氣當量比(ER)對合成氣熱物性軸向分布的影響Fig.6 Effect of air equivalence ratio (ER) on the axial distribution of thermal properties of syngas(a) Temperature; (b) Density; (c) Specific heat capacity (c); (d) Thermal conductivity (λ)Conditions: SBR=2.2; Tw=1123 K

熱導率是衡量物質傳熱能力的重要表征參數,而流化床生物質氣化反應過程中熱量的傳遞包括床料-流化劑、氧化產物-爐內氣體以及爐內氣體-流化劑等多種途徑。高溫氣氛通過提升氣體分子的運動強度進而加強氣相的熱傳導效應。可燃物的氧化提升了床內溫度,增加了氣體混合物的熱導率。然而,合成氣的熱導率隨著蒸汽的引入急速減小。

3.3 SBR對氣化反應的影響

3.3.1 SBR對氣化產物的影響

在ER為0.20、Tw為1123 K時,探究不同SBR下氣化產物沿軸向的分布規律,結果如圖7所示。蒸汽的引入增加了反應器中H、O元素含量,影響了均相和非均相反應。由于通入的蒸汽隨流化介質向上流動,因此4種SBR條件下的揮發分在低于蒸汽入口高度時均保持一致。隨著高度的繼續增加,各氣體組分含量受注入蒸汽的影響而顯著降低。然而,蒸汽與焦炭的非均相反應消耗了少量蒸汽,此時氣化產物濃度有略微增加,并在一定高度后保持穩定。當SBR由1.4增大至2.6時,有且僅有蒸汽在出口處的質量分數由28.2%增加至46.5%,而C2H4含量基本不變,其余氣體含量均不斷減小。

w—Mass fraction of gas species; z—Axial height圖7 蒸汽/生物質比(SBR)對氣體組分質量分數軸向分布的影響Fig.7 Effect of steam to biomass ratio (SBR) on the axial distribution of mass fraction of the gas species(a) C2H4; (b) CH4; (c) CO; (d) H2; (e) CO2; (f) H2OConditions: ER=0.20; Tw=1123 K

3.3.2 SBR對氣體熱物性的影響

在ER為0.20、Tw=1123 K的條件下,氣體熱物性參數在不同SBR時沿軸向的分布曲線見圖8。由圖8可知,蒸汽注入反應器后隨流化介質作上升運動,因此SBR對反應器下半部分的熱物性基本無影響。蒸汽通入量的增大強化了低溫蒸汽與高溫氣氛之間的熱傳遞效應,爐內溫度隨之逐漸降低。溫度的變化促使混合氣相密度隨著SBR的增加逐漸增大。與之相反,蒸汽流量的提升降低了氣體分子的運動活性,減弱了混合氣相的傳熱能力,因此氣相熱導率在持續增大的低溫蒸汽流量的作用下不斷減小。此外,蒸汽流量對混合氣相比熱容僅有輕微影響。

圖8 蒸汽/生物質比(SBR)對合成氣熱物性軸向分布的影響Fig.8 Effect of steam to biomass ratio (SBR) on the axial distribution of thermal properties of syngas(a) Temperature; (b) Density; (c) Specific heat capacity (c); (d) Thermal conductivity (λ)Conditions: ER=0.20; Tw=1123 K

3.4 溫度對氣化反應的影響

在ER為0.20、SBR為2.2的條件下,合成氣含量在不同氣化溫度時沿軸向的分布規律見圖9。在均相和非均相反應的綜合影響下,氣體濃度保持不變或隨溫度增大/減小的變化趨勢體現了氣化反應的復雜性。氣化溫度的升高促進了H2和CO2的生成。而氣化溫度的變化對CH4、CO和H2O含量僅有輕微的影響,這可能是因為升高溫度提升了氧氣與可燃物的燃燒強度和蒸汽氣化的反應速率。

w—Mass fraction of gas species; z—Axial height圖9 氣化溫度(Tw)對氣體組分質量分數軸向分布的影響Fig.9 Effect of temperature (Tw) on the axial distribution of mass fraction of the gas species(a) C2H4; (b) CH4; (c) CO; (d) H2; (e) CO2; (f) H2OConditions: ER=0.20; SBR=2.2

在ER為0.2、SBR為2.2的條件下,氣化溫度對混合氣體熱物性軸向分布的影響規律見圖10。氣體通入量一定時,氣化溫度的升高會提高整個爐內的溫度和熱導率,但卻減小了氣體的密度。比熱容和熱導率在流化區具有相同的變化趨勢,即在熱質傳遞的過程中沿軸向逐漸增大,且2種變量在不同條件下的變化趨勢與氣化溫度正相關。

圖10 氣化溫度(Tw)對合成氣熱物性軸向分布的影響Fig.10 Effect of temperature (Tw) on the axial distribution of thermal properties of syngas(a) Temperature; (b) Density; (c) Specific heat capacity (c); (d) Heat conductivity (λ)Conditions: ER=0.20; SBR=2.2

4 結 論

采用MP-PIC模型對實驗室尺度流化床內松木屑氣化反應進行了三維數值模擬,基于良好的計算與實驗結果的對比,研究了氣化過程中顆粒在爐內的分布規律,探究不同操作參數對氣化產物、氣相溫度、密度、比熱容和熱導率的影響。基于數值模擬結果得到以下結論:

(1)固體物料的質量差異使得松木屑和床料出現分離,較重的的床料堆積在床底。由于氧化反應,密相流化區內顆粒溫度隨高度增加,生物質氣化過程的吸熱反應使自由空域顆粒的溫度低于密相區。

(2)增大ER提升了可燃物的氧化強度,但過大的氣速縮短生物質顆粒停留時間而使其氣化過程難以充分進行,且N2的通入稀釋了合成氣的含量。故可燃氣體的含量隨著ER的增大而減小,而CO2含量、爐內溫度、氣相熱導率隨著ER的增大先增加后減小。

(3)增大低溫水蒸氣的流量,自由空域的氣相溫度和熱導率逐漸減小,而氣相密度持續增大。但SBR對比熱容的影響可忽略不計。

(4)氣化溫度的增加顯著提升了H2和CO2的含量,但減小了C2H4的含量。混合氣體的比熱容和熱導率與氣化溫度呈正相關,而密度與氣化溫度呈負相關。

主站蜘蛛池模板: 在线国产资源| 最新亚洲人成无码网站欣赏网| 99免费视频观看| 亚洲人成在线免费观看| 国产91视频免费观看| 69av免费视频| 欧美成人看片一区二区三区| 午夜色综合| 午夜综合网| 国产精品乱偷免费视频| 伊人国产无码高清视频| 直接黄91麻豆网站| 欧美日韩一区二区三区在线视频| 澳门av无码| 亚洲日韩精品综合在线一区二区| 女人av社区男人的天堂| 国产呦视频免费视频在线观看 | 天天色综网| 超碰91免费人妻| AV无码无在线观看免费| 国产美女精品在线| 欧美a√在线| 99re热精品视频国产免费| 久久黄色毛片| 97在线视频免费观看| 九色91在线视频| 日本免费福利视频| 国产美女视频黄a视频全免费网站| 好吊妞欧美视频免费| 亚洲人成在线免费观看| 国产91透明丝袜美腿在线| 成人毛片免费在线观看| 为你提供最新久久精品久久综合| 日韩人妻精品一区| 中国国产高清免费AV片| 91激情视频| 久久精品免费看一| 亚洲精品无码在线播放网站| 国产原创自拍不卡第一页| 国产一区二区影院| 欧美成人区| av尤物免费在线观看| 91啦中文字幕| 亚洲精品在线91| 久久精品中文字幕少妇| 玩两个丰满老熟女久久网| 国产Av无码精品色午夜| 亚洲小视频网站| 国产成人精品一区二区三区| 大香伊人久久| 免费一级全黄少妇性色生活片| 午夜小视频在线| 91精品久久久久久无码人妻| 又污又黄又无遮挡网站| 日韩小视频在线观看| 97国产在线播放| 中文字幕日韩视频欧美一区| 72种姿势欧美久久久大黄蕉| 国产精品第一区在线观看| 日韩小视频在线播放| 成人综合网址| 国产高清又黄又嫩的免费视频网站| 亚洲天堂网2014| 亚洲精品成人福利在线电影| 91蜜芽尤物福利在线观看| 国产超薄肉色丝袜网站| 亚洲无码精品在线播放| 成人精品区| 最新日韩AV网址在线观看| 亚洲欧美另类专区| 在线观看亚洲成人| 91在线免费公开视频| 内射人妻无码色AV天堂| 亚洲一区二区精品无码久久久| 亚洲69视频| 91九色国产在线| 久久国产精品波多野结衣| 国产精品第页| 亚洲欧洲AV一区二区三区| 免费a级毛片18以上观看精品| 亚洲成A人V欧美综合| 日本免费精品|