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

低旋流數旋進射流流場結構演變的POD分析

2021-09-15 08:11:56何創新劉應征
實驗流體力學 2021年4期
關鍵詞:模態區域結構

付 豪 ,何創新,*,劉應征

1.上海交通大學機械與動力工程學院 葉輪機械研究所,上海 200240

2.上海交通大學 燃氣輪機研究院,上海 200240

0 引 言

旋進射流被廣泛應用于燃燒室中,用于降低NOx排放以及改善火焰穩定性[1-2]。而低旋流(Low swirl)燃燒作為一種新穎的貧預混燃燒技術,與無旋和高旋流燃燒相比,具有更低的NOx排放[3-5]。實際上,旋進與旋流兩者的結合在預混燃燒室中更加普遍。無旋射流經過噴嘴流入大的軸對稱腔體時會發生偏轉并不對稱地再附到壁面,此后,射流會繞腔體軸線在整個腔體中旋進。這種主要的旋進流動結構會在腔體內部循環,對燃燒性能影響很大。毫無疑問,增加了旋流條件的旋進射流,其流場結構會變得更加復雜,并進一步增強燃料摻混。因此,研究低旋流數旋進射流的流場結構對了解摻混機理和強化控制尤為重要。

針對無旋的旋進射流,已經有一些學者對其進行了細致的研究。Nathan等[6]對旋進射流進行了詳細描述,發現流動是雙穩態的(Bistable),在旋進射流(Precessing Jet,PJ)模式和軸對稱射流(Axial Jet,AJ)模式之間間歇性地、隨機地切換。前者流場的特征是一個連續不穩定的再附射流,同時腔體另一側存在一個強回流區;后者流場與普通軸對稱射流相同。由于流動處于中性平衡狀態,任何不對稱和/或湍流波動都會導致再附著點移動,從而發生旋進。當射流開始沿一個方向旋轉時,進入腔室中的流體的不對稱性導致旋轉壓力場產生。而再循環流體在腔體中逆流而上,其旋轉的方向與旋進相反,從而保持凈角動量等于零,同時極大地促進了腔體內部流場的摻混[7]。此后,一些研究聚焦于確定有利于旋進射流的最佳幾何[2,8-9]。當雷諾數Re>1.0× 104、膨脹比(腔體直徑D與噴嘴直徑d的比值)為5時,旋進僅發生在腔體長徑比(腔體長度L與直徑D的比值)2.00~3.50的范圍內,且當長徑比為2.75時,旋進現象最為明顯。然而過去的研究中,對腔體內流場中的拓撲結構研究較少。Wong等[10]利用相位平均激光多普勒測量技術確定了旋進射流腔體的內部流場,發現PJ模式下的內部速度場比AJ模式下衰減更快。Cafiero等[11]通過層析PIV測量并結合POD分析來描述內部流場,發現流場中前三階模態占主導作用,其中能量最大的前兩階模態與射流旋進動有關,而第三階模態則與軸向運動有關。Ceglia等[12]的研究則發現在旋進中卷吸過程會影響大尺度相干結構的瞬時組織。

低旋流燃燒是一種利用低旋流的流場維持預混火焰穩定的燃燒方式。最初是由香港理工大學的Chan等[13]在1992年發現的,他們將切向射流法產生的低旋流場應用于CH4預混燃燒,發現流場中并沒有出現回流區,而預混火焰依然能夠維持穩定,同時高溫煙氣的停留時間縮短,降低了NOx排放。隨后美國勞倫斯伯克利國家實驗室的 Cheng[14]研究了低旋流燃燒的流場特性,結果表明低旋空氣流誘導的離心力作用在主流預混氣體上,導致主流速度逐漸降低,最終在燃燒器出口處形成發散流場(Diverging flow)。這種發散的流場具有幾個特性:1)氣流速度沿軸向逐漸衰減;2)火焰面穩定在當地湍流火焰傳播速度等于當地氣流速度的位置;3)湍流火焰傳播速度與湍流強度呈線性對應關系。

然而,目前對低旋流數旋進射流流場的研究還比較有限。Dellenback等[15]發現低旋流數旋進射流中旋進方向與旋流方向相反,但并未對其流場結構進行分析。Markovich等[16]利用POD分析了受限空間低旋流火焰中的結構,發現主要的相干結構以傾斜的環渦形式出現,但其流場并未出現旋進現象。為此,本文利用粒子圖像測速技術對3個不同低旋流數的旋進射流流場進行全場測量,在前面工作的基礎上[17],著眼于利用本征正交分解(POD)提取流場中含能大尺度結構,并對比分析流場中占主導的動態結構及其演化過程。

1 實驗方法

實驗所用的射流循環系統如圖1(a)所示。實驗系統固定在一個大的玻璃水箱中,其尺寸長3000 mm、寬550 mm、深700 mm。該系統由潛水泵、變壓器、圓柱形沉降室、旋流器、軸對稱腔室和一些管路組成。流動由安裝在水箱底部的潛水泵驅動,通過變壓器可以調節流速。射流從內徑為40 mm的圓管中流經軸向旋流器,然后射入到軸對稱腔體中。實驗裝置的幾何尺寸如圖1(b)所示。這里使用的旋流器與文獻[18]中使用的相同,內徑d為40 mm。實驗中選取了旋流數較小的3個旋流器(S= 0、 0.26和 0.41)。軸對稱腔體的長度L和內徑D分別為550 mm和200 mm,膨脹比D/d=5, 腔體長徑比L/D=2.75。在實驗中,將玻璃水箱里裝滿自來水以確保水位遠遠超過軸對稱腔體的頂部。通過變壓器調節流速后保持不變,確保不同噴嘴流動的雷諾數相同,根據旋流器內徑d以及給水管中的平均速度U0求 得Re= 4.5×104。

圖1 實驗裝置示意圖和幾何尺寸Fig.1 Schematic diagram of the experimental setup and geometric sizes

利用平面PIV測量了流向(x-y)平面內的流場。實驗中用密度為1.04 g/mm3的空心玻璃珠作為示蹤粒子,粒徑約為10 μm。照明光源為5 W的連續半導體激光器(波長532 nm),測試區片光源厚度約為1 mm。實驗用高像素密度CCD相機(IPX 16M,IMPERX,USA)捕獲粒子的流動,分辨率為4872 pixel×3248 pixel。PIV測量的區域從噴嘴出口延伸到腔體出口附近,由于激光強度的限制,照亮的區域有限,因此分別測量區域Ⅰ和區域Ⅱ的速度場,如圖1(b)所示。它們的流向測量范圍分別為0~7d和5~14d。在實驗中,相機以1 Hz的頻率對每個區域連續記錄了2000幅流動圖像。使用標準的PIV互相關算法獲得瞬態速度場,其中判讀窗口大小為32 pixel×32 pixel,相鄰窗口重疊率為50%,空間矢量分辨率為1.4 mm×1.4 mm。數據處理過程中,還采用了窗口偏移、亞像素識別和失真校正技術,使得兩幅圖像間粒子位移測量誤差小于0.1 pixel。

2 結果與討論

為了確定低旋流數旋進射流流場中的大尺度流動結構,對其速度場進行了本征正交分解。POD的詳細信息(包括其基本原理以及數學程序)可以參考Sirovich[19]的文章。POD分解的目的是通過提取最高能量的特征模態對流場動態特性進行描述。本文采用快照POD方法,分別對不同旋流數下、區域Ⅰ和Ⅱ中的1 000個瞬態場進行POD分解,分解后得到1 000階模態和對應的特征值。由于POD模態是流場能量最大化的一種形式,前幾階POD模態體現了流場中大尺度含能結構的動態特性,因此分析前幾階模態的特性對了解流場的基本結構具有重要意義。

圖2展示了不同旋流數下區域Ⅰ和區域Ⅱ的前十階POD模態的能量分布(即歸一化后的特征值λn,n為階數)。從圖中可以看到,低模態具有較高的能量,且隨著模態階數的增加,特征值衰減很快,前六階模態之后能量占比迅速減少。在區域Ⅰ中,不同旋流數之間的模態能量分布比較明顯。無旋時,第一階模態能量占比約為11.4%,遠超第二階模態能量占比的4.2%。而在旋流條件下,第一、二階模態的能量衰減隨旋流數增加而減弱。在旋流數為0.26和0.41時,第一、二階能量分別從13.0%衰減到7.4%和9.1%衰減到4.7%。這可能是由于旋流的存在,流場變得更加復雜,使得無旋時流場中占主導地位的大尺度結構的主導性下降。區域Ⅱ中的結果也體現了這一特點。對比同一旋流數下區域Ⅰ和Ⅱ的結果,可以發現,在區域Ⅱ中,前兩階模態的能量較區域Ⅰ有所增加,說明非定常特性在下游區域的優勢有所增強。

圖2 POD模態前十階特征值Fig.2 Eigenvalues of the first ten POD modes

由于前六階模態具有比較突出的能量,可以表征流場中的大尺度空間結構,因此選取前六階模態進行深入分析。圖3展示了無旋時區域Ⅰ內的前六階模態,云圖表示流向脈動速度(u')大小。含能最高的第一階模態如圖3(a) 所示,圖中有兩個符號相反的速度脈動明顯的區域,它們對稱分布在射流軸線的兩側。腔體上半部分正的速度脈動明顯的區域表明流體從那里直接流向下游;腔室下半部分負的速度脈動明顯的區域則表明該區域是回流。這說明旋進導致流動在測量平面里交替從一側流出、從另一側流入。圖3(b)中所示的第二階模態中,有兩個負的速度脈動明顯的區域,說明區域Ⅰ基本全是回流,這種結構可能與腔體壁面限制引起的大的回流有關。圖3(c)中所示的第三階模態中,區域Ⅰ末端有兩個符號相反的速度脈動明顯的區域,與第一階模態相比,區域面積變小;同時在這兩個區域上游不遠處出現了一個逆時針渦。兩個符號相反的速度脈動明顯的區域同樣出現在第五階模態中,在圖3(e)中,這兩個區域向上游移動;同時上游的逆時針渦消失,在上下兩側出現了兩個小的順時針渦。第三、五階模態可能代表旋進起始位置附近的小尺度結構。第四、六階模態表示的結構比較復雜,很難將它們與特定的物理機制相關聯,這是由于POD模態是按能量區分,某一模態可能是由多個物理現象疊加而得到。

圖3 旋流數 S = 0 時區域Ⅰ內的前六階模態Fig.3 The first six POD modes in zone Ⅰ at swirl number S = 0

為了更好地理解低旋流數旋進射流大尺度結構演變過程,利用前六階空間模態來重構脈動速度場。圖4展示了無旋時區域Ⅰ內4個典型時刻的重構脈動速度場,為了更好地體現演變過程,原始瞬態速度(u)場也同時給出。從圖4(a)可以看到此時旋進已經發生,主流在距離噴嘴出口大約5d處開始向軸線下方偏轉,體現了旋進時主流再附到壁面的過程。進一步通過圖4(e)可以看到,在旋進發生處存在兩個脈動明顯的區域,導致流體從軸線下方流出,從上方流入。這種流動導致該區域兩側同時形成了兩個順時針渦。這個結果與上面的第三階模態空間結構比較相似,印證了之前推測:此結構表征了旋進起始區域的流動結構。此外,從圖4(e)中可以清晰地看到射流剪切層中的相干渦結構向下游發展,直至旋進發生的區域,渦開始大幅度地向軸線兩側移動。圖4(b)展示了旋進更為強烈的一個階段,這一點是通過觀察主流偏轉另一側的回流速度大小來進行判斷的。這一結論在圖4(f)中可以看得更加清晰。此時,速度脈動明顯的區域向上游發展并變得很大,幾乎占據了整個區域Ⅰ的后半部分,說明此時旋進的起始位置開始向上游發展,下游整個區域幾乎都受到了旋進的影響,導致流體從下半部分流出、上半部分流入。此時,剪切層中的相干渦結構幾乎消失不見,說明旋進會使剪切層中的渦結構受到破壞,導致射流處于一個極度不穩定的狀態。圖4(c)和4(d)則體現了主流向另一側再附的過程,說明旋進會導致主流在測量平面上下來回擺動,而不是僅僅再附到一側,這是由于再附點的周向不穩定性造成的[6];而從圖4(g)和4(h)中可以看到,在旋進過程中,剪切層渦不再向兩側移動,而是在軸線附近交替扭曲著向下游發展。

圖4 原始流場(a)-(d)及重構脈動速度場(e)-(h)云圖(旋流數 S = 0,區域Ⅰ)Fig.4 Contour plot of (a)-(d) the original field and (e)-(h) the reconstructed fluctuating velocity field (swirl number S = 0, zone Ⅰ)

圖5是對無旋時區域Ⅱ內的流場進行POD分解后得到的前六階空間模態。圖5(a)中的第一階模態與區域Ⅰ中的第一階模態流場結構十分相似,都有兩個大的符號相反的速度脈動明顯的區域,它們幾乎對稱分布在軸線兩側,表明流動從一側流出、從另一側流入。相同的流動結構表明它們表征同一個物理現象,即旋進現象。在圖5(b)中的第二階模態中,軸線中心下游處存在明顯的正速度波動區域,說明該結構表征遠場中強烈的流向振蕩。對于圖5(c)中的第三階模態,其流動與旋進非常相似:流體在靠近腔室壁面的位置流出,導致環境流體從腔室出口的中間流入,其整體結構是軸對稱的,這可能是由于該模態結構是旋進模態與軸對稱模態疊加的結果。從圖5(d)和5(f)中可以觀察到第四、六階模態中的流動結構十分相似,通常這種成對的POD模態是由對流引起的,并且這兩階模態代表的是相同的結構,僅僅流向上有空間位移。在第四階模態中,x/d=8.0和12.5處中心位置分別出現了逆時針和順時針的旋渦結構;而在第六階模態中,x/d=7.0、10.5和14.0處中心位置分別出現了順時針、逆時針和順時針旋渦結構。因此,這對耦合的POD模態表征了剪切層中的大尺度旋渦結構,而這些旋渦結構與內射流(Inner jet)的大尺度振蕩有關[20]。

圖5 旋流數 S = 0 時區域Ⅱ內的前六階模態Fig.5 The first six POD modes in zone Ⅱ at swirl number S = 0

圖6是對無旋時區域Ⅱ利用前六階POD模態重構的脈動速度場及其相對應的原始瞬態場。圖6(a)~(d)中4個不同時刻的瞬態流場給出了旋進射流在靠近腔體出口區域再附到腔體壁面并在流向平面上來回擺動的過程,這個擺動過程中還伴隨著扭曲,這是因為射流在下游存在著強烈的振蕩。在圖6(e)和6(h)中,一正一負脈動明顯的區域幾乎占據整個腔體,此時腔體內的流體與腔體出口附近的流體形成一個大的循環:腔體內流體再附到一側壁面并沿著壁面流出,導致腔體出口附近的流體從另一側流入腔體。圖6(f)和6(g)則體現了旋進從腔體一側壁面發展到另一側壁面的過程。此過程中脈動明顯的區域變小,且脈動強度較再附到壁面時有所減弱,說明旋進射流再附到壁面時,同周圍流體摻混的效果最好。此外,區域Ⅱ里流動結構相對區域Ⅰ簡單,沒有旋渦結構,這是因為旋進在下游的發展更加強烈,導致剪切層內的旋渦結構被完全破壞。

圖7是對旋流數為0.26時區域Ⅰ內的流場進行POD分解后得到的前六階空間模態。圖7(a)中的第一階模態與無旋時區域Ⅰ的第一階模態相似,說明這兩個旋流數下區域Ⅰ流場內占主導地位的流場結構是相同的,即旋進。但與無旋時相比,圖7(a)中兩個脈動明顯區域的位置更靠近噴嘴出口,說明旋進發生的起始位置向上游移動。而圖7(b)和7(c)中的第二、三階模態不再同無旋時區域Ⅰ對應的模態相似,反而與無旋時區域Ⅱ的第二、三階模態相似,分別表示射流遠場的流向振蕩以及旋進射流兩種模式的疊加。這可能是由于旋流使軸向速度衰減加劇,從而導致此時區域Ⅰ的位置相當于無旋時區域Ⅱ的位置。而在第四至六階模態中,并沒有像無旋時區域Ⅱ內一樣存在一對模態,這也是由于旋流的存在導致軸向速度衰減,使得射流自身的大尺度振蕩減弱。

圖8是對旋流數為0.26時區域Ⅰ利用前六階POD模態重構的脈動速度場及其相對應的原始瞬態場。從圖8(a)~(d)中可以清楚地看到旋流數為0.26時區域Ⅰ內旋進流場的演變過程:射流在平面內來回擺動。與無旋時相比可以看到,旋流數為0.26時,旋進起始位置向上游移動,進而導致射流幾乎在區域Ⅰ的尾端再附到壁面上。此時射流的偏轉角度也有所增加。此外,在圖8(b)中還發現了無旋時沒有的一個流動現象:此時射流的主流是向下側偏轉,而下游的發展并不是同主流一樣偏轉到下側直至再附到腔體壁面,而是向上側發展。從圖8(f)中可以看到,上游剪切層中旋渦結構還未被破壞,而到了x/d>3.0的區域后,剪切層被完全破壞。這說明此時上游未處于旋進狀態,下游卻處于旋進狀態。這也體現了旋進射流在旋進模式和軸對稱模式之間混亂地切換,導致流場變得十分復雜[6]。

圖8 原始流場(a)-(d)及重構脈動速度場(e)-(h)云圖(旋流數 S = 0.26,區域Ⅰ)Fig.8 Contour plot of (a)-(d) the original field and (e)-(h) the reconstructed fluctuating velocity field ( S = 0.26, zone Ⅰ)

圖9是對旋流數為0.26時區域Ⅱ內的流場進行POD分解后得到的前六階空間模態。可以看到所有模態在x/d>8.0的區域幾乎沒有脈動明顯的區域,這可能是旋流引起的軸向速度衰減引起的。第一、二階模態與區域I的結果比較相似,分別代表旋進和流向振蕩。而第三、四階模態中,未發現比較明顯的結構,雖然有脈動明顯的區域,但很難與某一物理現象相結合。在第五階模態中,可以看到兩個正的脈動明顯區域之間存在一個負的明顯脈動,從而產生了兩個方向相反的旋渦,第六階模態稍下游的位置也發現相似的結構,可以認為是同一對流結構。這一對流結構明顯與無旋時區域Ⅱ中射流自身的振蕩不同,是由于旋流產生的一些旋渦結構向下游發展輸運的結果。

圖9 旋流數 S = 0.26 時區域Ⅱ內的前六階模態Fig.9 The first six POD modes in zone Ⅱ at swirl number S = 0.26

圖10是對旋流數為0.26時區域Ⅱ利用前六階POD模態重構的脈動速度場及其相對應的原始瞬態場。瞬態場的演變及相對應的流場結構與無旋時區域Ⅱ的差別不大,這是因為旋進在下游比較強烈,導致射流剪切層被完全破壞,流場中的大尺度結構均與旋進相關,即一正一負脈動明顯的區域,導致流體從腔體一側流出從另一側流入。與無旋時相比,此時的脈動速度強度有所增加,這也說明了旋流的存在會增強旋進的強度。

圖10 原始流場(a)-(d)及重構脈動速度場(e)-(h)云圖(旋流數 S = 0.26,區域Ⅱ)Fig.10 Contour plot of (a)-(d) the original field and (e)-(h) the reconstructed fluctuating velocity field ( S = 0.26, zone Ⅱ)

圖11和12分別是對旋流數為0.41時區域Ⅰ內的流場進行POD分解后得到的前六階空間模態以及用前六階POD模態重構的脈動速度場及其相對應的原始瞬態場。可以看到,當旋流數為0.41時,區域Ⅰ中的前六階空間模態均能在無旋時或者旋流數為0.26時找到相似的模態,這里不做詳細分析。但對比不同旋流數下結構相似的模態可以發現,隨著旋流數的變化,這些結構出現的位置向上游移動,這也是由于旋流引起流向速度衰減導致的。同時從瞬態場可以看出,隨旋流數增加,旋進起始位置向上游移動且旋進的偏轉角度逐漸增大,強度也有所增強。

圖11 旋流數 S = 0.41 時區域Ⅰ內的前六階模態Fig.11 The first six POD modes in zone Ⅰ at swirl number S = 0.41

圖12 原始流場(a)-(d)及重構脈動速度場(e)-(h)云圖(旋流數 S = 0.41,區域Ⅰ)Fig.12 Contour plot of (a)-(d) the original field and (e)-(h) the reconstructed fluctuating velocity field ( S = 0.41, zone Ⅰ)

圖13和14分別是對旋流數為0.41時區域Ⅱ內的流場進行POD分解后得到的前六階空間模態以及用前六階POD模態重構的脈動速度場及其相對應的原始瞬態場。此時旋流數已經大到使區域Ⅱ中脈動明顯的區域位于與區域Ⅰ重疊的部分,而更下游區域的速度及其脈動強度都非常小,沒有明顯的結構,因此對區域Ⅱ的結果也不做詳細討論。

圖13 旋流數 S = 0.41時區域Ⅱ內的前六階模態Fig.13 The first six POD modes in zone Ⅱ at swirl number S = 0.41

圖14 原始流場(a)-(d)及重構脈動速度場(e)-(h)云圖(旋流數 S = 0.41,區域Ⅱ)Fig.14 Contour plot of (a)-(d) the original field and (e)-(h) the reconstructed fluctuating velocity field (S = 0.41, zone Ⅱ)

3 結 論

本文利用PIV技術測量了3個不同旋流數下的旋進射流流場,并結合本征正交分解提取流場中含能大尺度結構,對比分析流場中占主導的動態結構及其演化過程,得到以下結論:

1)通過POD方法成功提取了受限射流中大尺度旋進現象導致的平面內射流體上下振蕩的現象,而第一階POD模態能量的相對大小是旋進現象強弱的標志。

2)旋進剛發生時,上游剪切層內的旋渦結構尚未完全破壞,它們會一直向下游發展直至旋進起始點附近后,開始隨著主流一起偏轉,而下游剪切層內的大尺度結構被完全破壞。

3)旋進射流流場中既包含大尺度的旋進現象,也包含普通射流中的流向振蕩以及縱向振蕩,而入口的旋流增強了射流自身振蕩,從而導致流場更加復雜,加速了大尺度結構的破壞。

猜你喜歡
模態區域結構
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
論《日出》的結構
關于四色猜想
分區域
國內多模態教學研究回顧與展望
基于嚴重區域的多PCC點暫降頻次估計
電測與儀表(2015年5期)2015-04-09 11:30:52
創新治理結構促進中小企業持續成長
現代企業(2015年9期)2015-02-28 18:56:50
基于HHT和Prony算法的電力系統低頻振蕩模態識別
由單個模態構造對稱簡支梁的抗彎剛度
計算物理(2014年2期)2014-03-11 17:01:39
主站蜘蛛池模板: 天天综合网色中文字幕| 亚洲精品国产综合99久久夜夜嗨| 在线一级毛片| 精品国产免费观看| 中国成人在线视频| 国产h视频在线观看视频| 人妻精品久久久无码区色视| 亚洲综合第一区| 一级成人a做片免费| 久久久久国产一级毛片高清板| 91成人在线免费视频| 成人亚洲国产| 国产微拍一区二区三区四区| 久久久精品无码一区二区三区| 67194成是人免费无码| 国产精品3p视频| 国产人人射| 欧美第二区| 爱色欧美亚洲综合图区| 青青青视频91在线 | 久久国产香蕉| 精品伊人久久久久7777人| 中文字幕1区2区| 免费人成网站在线高清| 精品人妻一区二区三区蜜桃AⅤ| 喷潮白浆直流在线播放| 亚洲男人天堂2018| 欧美亚洲国产一区| 71pao成人国产永久免费视频 | 野花国产精品入口| 免费在线a视频| 亚洲综合九九| 激情无码字幕综合| 欧美成一级| 国产日韩久久久久无码精品| 91丨九色丨首页在线播放| 91在线日韩在线播放| 中文字幕免费视频| 亚洲精品在线影院| 精品国产福利在线| 亚洲AⅤ无码日韩AV无码网站| 四虎影视无码永久免费观看| 小13箩利洗澡无码视频免费网站| 欧美日韩成人在线观看| 国产正在播放| 欧美日韩专区| 久久精品亚洲热综合一区二区| 亚洲 欧美 偷自乱 图片| 91久久偷偷做嫩草影院电| 色综合久久无码网| a级毛片视频免费观看| 国产伦片中文免费观看| 亚洲码在线中文在线观看| 成年片色大黄全免费网站久久| 国产成人亚洲毛片| jizz在线观看| 日日拍夜夜嗷嗷叫国产| 亚洲男人天堂2018| 手机在线免费不卡一区二| 一本大道视频精品人妻| 久久国产精品麻豆系列| 亚洲色图欧美视频| 国产真实二区一区在线亚洲| 真实国产乱子伦视频| 亚洲综合久久一本伊一区| 五月综合色婷婷| 国产精品美女在线| 亚洲男女在线| 亚洲美女视频一区| 日本不卡视频在线| 国产国拍精品视频免费看 | 欧美一区二区啪啪| 亚洲欧洲日本在线| 毛片久久网站小视频| 91激情视频| 欧美日韩一区二区在线播放| 国产女人水多毛片18| 国产午夜福利片在线观看| 国产精品欧美激情| 亚洲综合第一页| 欧美成人二区| 国产微拍一区|