亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        循環(huán)流化床鍋爐爐膛流動特性數值模擬進展

        2021-11-10 03:25:22劉雪敏張肖陽齊國利
        煤炭學報 2021年10期
        關鍵詞:曳力流態(tài)爐膛

        楊 照,劉雪敏,張肖陽,齊國利,董 勇

        (1.山東大學 燃煤污染物減排國家工程實驗室,山東 濟南 250061; 2.中國特種設備檢測研究院,北京 100029)

        爐膛內的流動特性是循環(huán)流化床鍋爐(CFBB)燃燒反應和傳熱的基礎,具有重要的研究價值和意義,并得到了充分的研究[1-4]。整體和局部流動結構的不均勻性是爐膛氣固兩相流動最顯著的特征,根據已有的研究,爐膛整體由底部高固體顆粒體積分數的密相區(qū)和上部低固體顆粒體積分數的稀相區(qū)組成,孔隙率分布呈現為軸向的S型結構和徑向的環(huán)-核結構,而爐膛局部顆粒則處于不斷地聚集和分裂狀態(tài),具有團聚趨勢。

        然而實驗方法存在研究成本高、獲得數據有限等不足,盡管測量工具和研究方法隨著科技進步得到了一定的發(fā)展,但截至目前,爐膛內的許多關鍵參量或是難以通過實驗測得,或是存在一定的局限性。例如爐膛中的壓力尚可通過實驗測得,固體體積分數也可以通過所測的壓降推導,但所得數據僅能反映某一高度或某一時刻的特定物料體積分數。因此,數值模擬方法受到了越來越多研究者的青睞。

        筆者分析了近10 a爐膛流動特性模擬研究的進展,重點分析了爐膛流動特性模擬的關鍵因素——曳力模型,對模擬中多相流模型和曳力模型的耦合進行闡述,分析當前研究存在的不足,探討未來CFBB氣固兩相流動模擬的發(fā)展方向。

        1 爐膛流動特性模擬研究進展

        爐膛內固相的流態(tài)化現象是CFBB最典型的特征,模擬須從氣固多相流模型入手。

        根據對固相的描述,目前已有的多相流模型分為基于歐拉-歐拉(E-E)方法的雙流體模型(Two Fluid Model,TFM)和基于歐拉-拉格朗日(E-L)方法的離散元(Discrete Element Model,DEM)模型以及計算顆粒流體力學(Computational Particle Fluid Dynamics,CPFD)模型。E-E方法基于連續(xù)介質假設(亦被稱為連續(xù)介質模型),將流體與離散的固相作為共同存在且相互滲透的連續(xù)介質;E-L方法則只將氣體相當作連續(xù)相,分別在歐拉坐標系和拉格朗日坐標系下描述氣相和固相。E-L方法中,DEM模型充分計算固相間的力和碰撞,詳細描述顆粒運動,能夠揭示顆粒相的宏觀流動特性和確定固相在系統(tǒng)內的停留時間;CPFD模型基于MP-PIC(Multi-Phase Particle In Cell)[5]方法,引入顆粒碰撞模型和空間梯度來計算固相間的相互作用力和碰撞,并將一部分位置、速度相近且粒徑、密度等物理量相同的顆粒作為整體進行計算[6-7]。

        相較于實驗,模擬可以捕獲更詳細的流動細節(jié),并提供可視化的結果,有助于理解整個系統(tǒng)內的氣固流動行為。目前已有許多研究者通過模擬得到了爐膛的典型流動特性,證明了模擬的可信程度。

        LIU等[8]和ZHANG等[9]在采用TFM模型的CFBB全回路模擬中,捕獲了爐膛內固相下濃上稀的分布特點,得到了爐膛內 “下濃上稀”和徑向“環(huán)-核”的分布特征(圖1);NIKOLOPOULO等[10]在研究中進一步將爐膛分為底部密相區(qū)域、中部飛濺區(qū)域和頂部自由發(fā)展區(qū)域3個區(qū)域(圖2)。而WANG等[11]在采用DEM模型對實驗室規(guī)模的CFBB全回路的模擬中除了捕獲了爐膛固相軸向S型分布(圖3)外,還發(fā)現湍流現象貫穿整個爐膛,團聚主要在爐膛底部壁面附近形成(圖4)。

        圖1 CFB鍋爐內固體體積分數Fig.1 Solid volume fraction distribution in CFBB

        圖2 爐膛內固體體積分數[10]Fig.2 Solid volume fraction in furnace[10]

        圖3 CFB中固體瞬態(tài)運動[11]Fig.3 Transient motion of solid phase in CFB[11]

        圖4 爐膛內固體運動及分布Fig.4 Movement and distribution of solids in furnace

        除了能直觀地反映爐膛內的宏觀流態(tài)特征外,通過模擬還能捕獲顆粒詳細的運動軌跡,用以研究系統(tǒng)內固相的微觀流動特性。

        LUO等[12]和WANG等[13]在采用DEM模型的CFBB模擬中觀察到了強烈的固相返混現象,且主要發(fā)生在爐膛壁面和角落處;FANG等[14]進一步將爐膛內固相的混合過程分為啟動階段、遷移階段、強烈混合階段和穩(wěn)定階段4個階段,并發(fā)現對流混合、剪切混合和擴散混合分別主要發(fā)生在爐膛的下部、中部和頂部,固相軸向擴散在整個區(qū)域中占主導地位。LUO等[12]還根據爐內固相平均速度捕獲了爐膛中心區(qū)域的固相向上運動和壁附近的固相向下流動,并將爐膛底部粒子無序運動的原因歸為固相反混、回料閥回流、流化風和粒子間相互作用。WANG等[15]還發(fā)現固相速度在爐膛中心區(qū)域和4個角處較大,而在環(huán)帶區(qū)域中較小(圖5)。

        圖5 爐膛局部區(qū)域固體速度[15]Fig.5 Solid velocity in the local area of furnace[15]

        模擬方法還能提供爐膛內更為詳細的物理參數,通過研究參數的動態(tài)變化過程,有助于加深對爐內氣固流動行為的理解和認識。

        曾勝庭等[16]在300 MW CFBB的 CPFD 模擬中,根據爐膛內的固相體積分數和顆粒速度分布證實了密相區(qū)后墻壁面附近存在的床料富集現象,并發(fā)現固相體積分數分布隨著爐膛高度的增加而逐漸降低,稀相區(qū)內的床料分布更為均勻,部分床料顆粒在氣流曳力的作用下被攜帶上升,另一部分貼近壁面的顆粒則呈現下降的流動狀態(tài)。李德波等[17]在300 MW循環(huán)流化床鍋爐的CPFD模擬中,根據爐膛速度矢量圖發(fā)現,爐膛密相區(qū)中部氣流速度在一次風和部分沿壁面返回風量的作用下大幅增加,同時壁面處出現數個中心區(qū)域速度極低的旋流,爐膛底部前壁處的顆粒體積分數較其他處更高。張瑞卿等[18]在采用CPFD模型的工業(yè)示范循環(huán)流化床鍋爐爐膛模擬中研究了爐內的物料體積分數,研究表明,固相物料體積分數在爐膛邊壁處存在明顯的上升趨勢,在爐膛下部中心區(qū)出現明顯的體積分數峰值,而體積分數峰值并沒有在爐膛頂端中心區(qū)出現。JIANG等[19]在對具有六旋風分離器布置的大型CFBB的CPFD模擬中研究了爐膛內的氣相速度場,觀察到爐膛出口處的氣流速度矢量方向從垂直轉變?yōu)樗?,爐膛底部氣流具有向上和向下流動的不均勻性,中心處固相的垂直速度高于壁面處的顆粒。另外,實驗研究中用于定義爐內流態(tài)的物理參數通常局限于空隙率、氣固滑移速度和方向,而通過數值模擬研究,有可能利用更詳細的物性參數完善對爐內流態(tài)的定義。TU等[20]采用CPFD多相流模型模擬了帶有6個旋風分離器的CFBB,并根據壓力時間序列等參數確定了爐膛中隨著表觀氣速增加而出現的3種流態(tài),即多氣泡流態(tài)(Multiple Bubble)、鼓泡流態(tài)(Bubble Regime)和湍流流態(tài)(Turbulent Fluidization Regime)。

        總而言之,通過數值模擬方法研究爐膛內的流動特性,研究者們已獲取更詳細的流動細節(jié)和可視化結果,加深了對爐膛內氣固兩相宏觀和微觀流態(tài)特征的認識,有助于完善爐內氣固兩相流態(tài)的定義,并會在今后的模擬研究中不斷深化對流動特性的認知。

        2 爐膛流動特性數值模擬的關鍵模型

        顆粒間的碰撞只是氣固兩相流模型需要處理的問題之一。研究證明,相對于對固相的準確描述,氣固兩相流模擬的關鍵在于系統(tǒng)中氣體-粒子和粒子-粒子相互作用的正確表達[21-23]。曳力作為多相流動中最重要的相間相互作用,決定了爐膛內的氣流對顆粒的夾帶和輸送以及內循環(huán)流動等行為[24-26]。因此,無論是E-E法或E-L法,曳力模型都是決定數值模擬準確性的關鍵。

        作為多相流動中最主要的相間相互作用,針對曳力已建立了大量的模型。目前曳力模型主要可以分為傳統(tǒng)均勻曳力模型、關聯型曳力模型和極值型曳力模型3類[27-29],表1對現有的曳力模型進行了歸納。

        表1 曳力模型的分類Table 1 Classification of drag model

        復雜系統(tǒng)中的介尺度關聯著微尺度和宏尺度,在CFBB系統(tǒng)中,微尺度、宏尺度和介尺度分別指顆粒尺度、設備尺度和顆粒團尺度。

        研究表明,爐膛內的固相趨向保持最小重力勢能,處于不斷地聚集、形成團聚并分裂的狀態(tài);氣相則趨向于選擇阻力最小的路徑,傾向于繞過而非穿透顆粒團[63-64]。2種運動機制之間的相互協(xié)調得到稀、密相共存的穩(wěn)定性條件,并實現了流動優(yōu)勢的折中[65],造成了爐膛內流動的復雜性和多樣性。同時,由于爐膛內的主要流動曳力來自密相團聚體內部及其界面,因此爐膛內的有效氣固動量交換低于均勻懸浮液[57]。

        綜上所述,介尺度結構是導致爐膛氣固流動結構非均勻性分布的根本原因,研發(fā)曳力模型的重點在于考慮介尺度結構和系統(tǒng)中固相的異質空間分布對曳力的影響[66-67]。

        2.1 均勻曳力模型及關聯型曳力模型

        均勻曳力模型的構建基于單顆粒曳力模型和均質系統(tǒng)(固定床或均勻液固流化床)的實驗結果,采用平均化方法對曳力系數進行修正,并應用于CFBB多顆粒體系模擬。但平均化方法簡化和忽略了爐膛多尺度結構的復雜性,會高估氣固相之間的曳力[68-69],更適用于顆粒體積分數較低、氣固流動較為均勻的體系。

        基于建立新曳力模型的困難和均勻曳力模型捕獲氣固兩相基本特征的能力[70],研究者們試圖通過修正均勻曳力模型以獲得能夠用以描述介尺度的次網格模型,即關聯型曳力模型[27]。根據對均勻曳力模型的修改依據可以將關聯型曳力模型分為修改顆粒直徑法、細網格歸納法和修正曳力系數法3類。

        修改顆粒直徑法將顆粒團當成一個整體,增大直徑的顆粒被當作顆粒團以考慮顆粒趨于成團的流動特性。采用修改顆粒直徑法的前提是準確描述顆粒團聚現象,但爐膛內的顆粒團受到氣相和固相的性質、爐膛的幾何形狀和運行條件等因素的影響,導致該方法不具有普適性,實用性和準確性也有待深入研究。

        細網格歸納法從概念上將一個復雜的多尺度閉合模擬任務分解為解決描述微尺度和介尺度的任務,采用傳統(tǒng)曳力模型對爐膛內的微小區(qū)域進行高精度模擬,通過求解區(qū)域內參數的時間平均得到適用于粗網格模擬的曳力模型,但其本質依然是均勻曳力模型。

        修正曳力系數法通過修正均勻曳力模型中的曳力系數,減小均勻曳力模型對于曳力的估計,根據修正方法的不同可以分為直接修正方法和經驗修正系數法2類。直接修正方法即通過給曳力模型添加一個小于1的系數減小曳力,缺點在于曳力系數的修正需要大量試算,且修正值取決于研究者的主觀意愿,缺乏完善可靠的理論依據;經驗修正系數法以O’BRIEN和SYAMLAL[54]提出的O-S模型為代表。O’BRIEN和SYAMLAL基于FCC系統(tǒng)實驗數據提出了顆粒團聚效應的修正公式,并以此修改均勻模型曳力參數,保障了模型的準確性。但由于需要根據模擬對象進行廣泛的研究而缺乏普遍性,O-S模型通常作為檢驗曳力模型的基準。

        均勻曳力模型和O-S曳力模型的曳力函數分布如圖6所示(其中,β為單位體積的曳力系數。定義為β=εgnsFD/ur,其中εg為氣體體積分數,ns為顆粒數密度,FD為顆粒所受曳力,ur為氣固相對速度[61])。圖6(a)匯總了部分均勻曳力模型,由圖6(a)可知,不同均勻曳力模型的曳力曲線幾乎完全重合,曳力隨顆粒體積分數增加而單調增大,沒有體現出非均勻流動中因介尺度的存在而導致的曳力降低。由圖6(b)可以看出, O-S曳力模型預測在顆粒體積分數εs=0.01~0.50內非均勻流的曳力會發(fā)生顯著的降低,且在εs≈0.10時曳力下降幅度最大,爐膛稀、密相兩側的流動趨于均勻,符合實驗現象[71]。

        圖6 均勻曳力模型與O-S曳力模型對比[28,71]Fig.6 Comparison of homogeneous model and O-S drag model under different Gs[28,71]Fig.6 Comparison of homogeneous model and O-S drag model Gs[28,71]

        2.2 極值型曳力模型

        極值型曳力模型以EMMS模型為代表的非均勻曳力模型。EMMS曳力模型基于“多尺度最小能量理論”(Theory of Energy Minimum Multi-scale,EMMS),核心在于對介尺度結構和多結構折衷特性的描述[72],爐膛非均勻區(qū)域被分為“顆粒密相”、“顆粒稀相”和“相互作用相”3個均勻的子系統(tǒng),通過穩(wěn)定性條件(單位質量顆粒的懸浮輸送能耗Nst趨于最小)描述介尺度上兩種主要趨勢的折衷[73]。

        原始EMMS模型由李靜海等提出[57],最初被應用于穩(wěn)態(tài)氣固流動和結構效應的量化[60],僅考慮了重力對曳力的影響。

        而后,楊寧等[58-59]考慮了顆粒的平均加速度和曳力、重力的不平衡,建立了Y-EMMS模型。Y-EMMS將顆粒局部加速度納入考慮,能求解特定工況下的曳力系數,但模型中包含經驗公式,且忽略了氣固兩相流速對曳力的影響。

        隨后,王維[60]進一步引入稠密相、稀相和相間的慣性加速用以區(qū)別稀、密相內部顆粒加速度,建立了EMMS/matrix模型。EMMS/matrix模型采用2步EMMS分析法,首先對循環(huán)流化床的操作參數進行最小能量分析得到顆粒團參數關于氣體體積分數的關系式,然后將關系式與CFD耦合得到非均勻曳力修正系數,綜合考慮了操作參數和局部流動參數對曳力的影響。但該模型密相顆粒的受力表達式并未采用實際的受力分析。

        后續(xù),祁海鷹等[28]基于顆粒團密度εsc修正了顆粒團直徑,通過完善局部顆粒的受力平衡提出了不含經驗系數的QL-EMMS模型和QL-EMMSn模型[61]。在祁海鷹等的基礎上,陳程等[62]進一步引入非均勻因子和基于整體氣固滑移速度的狀態(tài)雷諾數提出了QC-EMMS模型,實現了曳力模型的流動自適應,但該階段模型尚未考慮顆粒尺寸分布。

        EMMS曳力模型的對比如圖7所示。EMMS曳力模型在一定局部體積分數范圍內都預測到了不同程度的曳力下降,Y-EMMS模型預測的曳力系數在空隙率0.7~1.0內會大幅降低, EMMS/matrix模型預測曳力系數在除了爐膛兩端外的大部分空隙率范圍內都會降低,QL-EMMS曳力模型預測曳力在εs=0.01~0.35內顯著下降。

        圖7 不同曳力模型對比[77]Fig.7 Comparison of different drag model[77]

        但EMMS模型曳力函數隨顆粒體積分數的變化趨勢與基于實驗的相對準確可靠的O-S模型相比,或包含沒有物理意義的轉折,或存在不同的變化趨勢,原因主要在于對顆粒團直徑預測的缺陷[74-77]。

        合理描述顆粒團特性的QC-EMMS模型能夠體現非均勻流曳力隨顆粒體積分數呈“凹狀”上升,且在稀、濃兩端趨于均勻態(tài)的非均勻流曳力下降的本質特征(圖8)。

        圖8 QC/QL-EMMS模型與O-S模型的對比[62]Fig.8 Comparison of QC/QL-EMMS model and O-S model[62]

        決定CFBB爐膛內兩相流模擬成敗的關鍵和未來研究的重點仍是曳力模型[78]。采用均勻曳力模型、關聯型曳力模型將導致低估爐膛(提升管)底部稠密區(qū)域中固體的體積分數和高估顆粒循環(huán)流率,無法捕獲顆粒動態(tài)集群的形成和消失等特性[58],且即使在采用高網格密度的前提下仍然不能實現本質上的改變[79]。而EMMS曳力模型捕獲爐膛(提升管)內非均勻結構的能力更強,可以明顯區(qū)分和預測爐膛(提升管)內的稀相、密相區(qū)域及異質結構[20,80]。

        現有的EMMS曳力模型均未考慮顆粒粒徑的寬篩分特性,且目前用以檢驗的O-S模型基于FCC系統(tǒng)實驗數據而非CFBB系統(tǒng),同時研究爐膛內的氣固流動特性時需考慮各部件之間的相互影響。因此,現有的曳力模型還需要進一步改進以適用于CFBB全回路多流態(tài)的模擬。

        3 模型的耦合與選擇

        盡管多相流模型通常被認為僅對模擬結果產生次要影響[81],但出于科學的嚴謹性考慮,在CFBB的模擬中考慮多相流模型和曳力模型的耦合對于模擬結果的影響十分必要。考慮到CFBB屬于CFB的應用之一,而目前針對CFBB的模擬相對較少,因此對近20 a CFB模擬中采納的多相流模型、曳力模型及模擬對象和顆粒類型進行了歸納(表2~4)。

        表2 基于TFM多相流模型的CFB模擬Table 2 CFB Simulation based on TFM model

        固體顆粒依據顆粒粒徑和與流體密度之差可以分為 A,B,C,D 四類。CFB鍋爐系統(tǒng)中,A類顆粒飛灰僅在爐膛中短暫停留, B類顆粒作為外循環(huán)回路的主要成分將直接影響流動情況,故CFBB模擬的顆粒類型將以B類顆粒為主。

        TFM模型和均勻曳力模型或EMMS曳力模型耦合是目前被采納最多的方法。相比于與均勻曳力模型耦合,TFM模型與EMMS的耦合基于EMMS穩(wěn)定性條件約束,能更準確預測流態(tài),是目前準確性最高的模擬方法;模擬對象多為化學鏈燃燒(CLC)裝置、氣化爐和碳酸化器等小型/中試CFB提升管的研究,應用于大型CFB全回路模擬的研究較少;涉及的顆粒類型主要是A類和B類。但TFM模型基于連續(xù)介質假設,導致難以實現粒徑分布[60,105]和缺乏尺度分離[106-107],同時出于計算機性能考慮而簡化了計算[108-110],與真實情況存在一定差異。

        表4 基于CPFD多相流模型的CFB模擬Table 4 CFB Simulation based on CPFD model

        DEM模型在詳細追蹤顆粒運動的同時增加了計算壓力[97],因此多與均勻曳力模型耦合并應用于小型循環(huán)流化床的模擬,涉及的顆粒類型主要為粒徑較大的D類顆粒,不適用于CFB鍋爐的模擬。

        CPFD模型兼?zhèn)涮幚韺捄Y分、多粒徑的顆粒負載情況和追蹤單個顆粒、處理多分散顆粒群的能力,目前多與均勻曳力模型耦合,并被大量應用于B類顆粒和大型CFB系統(tǒng)的模擬,同時由于CPFD模型能夠模擬非流化顆粒的運動而被不少研究者應用于立管、返料閥等設備的研究。相較于完全依賴網格精度的TFM模型和大量追蹤粒子運動的DEM模型,CPFD模型兼?zhèn)涓咝屎蜏蚀_性,更適用于CFB鍋爐的模擬。

        CPFD模型目前主要的不足在于對曳力把握的不準確,但目前CPFD模型多與均勻曳力模型耦合,且因為高估曳力而產生誤差[24,78,111]。而將CPFD模型與EMMS模型耦合則能夠解決CPFD模型高估曳力的問題,能夠獲得更準確的模擬結果[100];同時CPFD模型中考慮了顆粒粒徑的寬篩分特性,得以彌補當前EMMS曳力模型的不足。EMMS曳力模型CPFD模型與EMMS模型耦合的方法在理論上兼顧對氣固兩相分布和氣固、固相之間作用力的準確描述,同時還優(yōu)化了數值計算效率,具有很高的應用前景。

        4 結 論

        (1)數值模擬能提供更詳細的流動細節(jié)和可視化結果,加深了對爐膛內氣固兩相宏觀和微觀流態(tài)特征的認識,根據爐膛內氣固兩相詳細的物理參數能完善對爐膛內流態(tài)的定義。

        (2)曳力模型是影響爐膛流動特性數值模擬的關鍵。曳力模型可分為均勻曳力模型、關聯型曳力模型和極值型曳力模型:均勻曳力模型忽視了介尺度和微尺度不均勻結構的影響,并將高估氣固間的曳力系數;關聯型曳力模型基于均勻曳力模型,無法從本質上解決均勻曳力模型高估曳力的問題;EMMS極值型曳力模型有效解決了介尺度的問題,并在近年得到了不斷的優(yōu)化和改進,是目前最好的曳力模型。但由于爐膛內氣固流動的復雜性及各部件之間的相互影響,目前的曳力模型仍存在不足,有待于建立符合爐膛流態(tài)的且適用于CFBB全回路系統(tǒng)的多流態(tài)曳力模型。

        (3)目前大多采用EMMS模型與TFM模型耦合。然而TFM模型對固相的描述不夠準確,無法捕捉流態(tài)的本質特征;相比較而言,CPFD模型對于多相流的模擬要更加優(yōu)越,發(fā)掘新的流動特性的可能更高。因此將EMMS曳力模型與CPFD多相流模型耦合是一種很有前景的循環(huán)流化床模擬方法,在未來會被更多研究者采用。

        猜你喜歡
        曳力流態(tài)爐膛
        基于鼓泡流化床的新型曳力模型的驗證分析
        預測天然氣斜井臨界攜液流量新方法
        側邊機組故障對泵站前池流態(tài)的影響
        大電機技術(2022年1期)2022-03-16 06:40:24
        二維爐膛氣液兩相對沖流動數值模擬
        層燃型垃圾焚燒鍋爐的爐膛與爐膛溫度簡析
        船用燃油輔鍋爐爐膛爆燃分析
        水上消防(2020年2期)2020-07-24 09:27:06
        改進邊界條件的非恒定流模型在城市河流橡膠壩流態(tài)模擬中的應用
        基于EMMS模型的攪拌釜內氣液兩相流數值模擬
        化工學報(2016年7期)2016-08-06 07:11:55
        動態(tài)流態(tài)冰蓄冷系統(tǒng)在千級凈化廠房的應用
        機電信息(2015年3期)2015-02-27 15:54:46
        基于TM遙感影像的河口流態(tài)信息半定量化研究
        欧美老肥婆牲交videos| 无码专区亚洲avl| 亚洲av网一区天堂福利| 中文字幕亚洲乱码熟女1区2区| 华人在线视频精品在线| 久久久久久久亚洲av无码| 欧美日韩视频在线第一区| 日韩一欧美内射在线观看| 无码视频一区=区| 美腿丝袜视频在线观看| 精品国产亚洲亚洲国产| 精品区2区3区4区产品乱码9| 精品久久久久久中文字幕| 一区二区三区福利在线视频| 热门精品一区二区三区| 亚洲av香蕉一区二区三区av| 国产在线精品一区二区三区| 久久人与动人物a级毛片| 国产精品丝袜在线不卡| 亚洲综合久久久中文字幕| 在线视频一区二区国产| 免费观看a级毛片| 亚洲精品无码久久久久牙蜜区| 含羞草亚洲AV无码久久精品| 亚洲av熟女天堂系列| 亚洲男人的天堂av一区| 真人抽搐一进一出视频| 亚洲精品92内射| 999精品免费视频观看| 精品日韩av专区一区二区| 亚洲国产av一区二区四季| 亚洲国产精品无码中文字| 国产乱子伦视频大全| caoporon国产超碰公开| 亚洲精品456在线播放狼人| 国产自国产自愉自愉免费24区 | 久久精品国产亚洲av麻豆会员| 丁字裤少妇露黑毛| 亚洲精品有码在线观看| 精品国产一区二区三区毛片| 91色老久久偷偷精品蜜臀懂色|