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

環(huán)形通道中湍流減阻技術(shù)及其機理研究

2018-12-26 07:55:34王春光尤軍鋒許桂陽鞏倫昆
關(guān)鍵詞:溝槽

王春光,尤軍鋒,許桂陽,鄧 哲,鞏倫昆

(1. 西安近代化學(xué)研究所,西安,710065;2. 西安航天動力技術(shù)研究所,西安,710025)

0 引 言

有關(guān)減阻的研究可追溯到20世紀30年代但直到20世紀60年代中期,研究工作主要集中在減小表面的粗糙度上。湍流非光滑減阻技術(shù)的研究開始于20世紀70年代,并隨著能源的緊缺越來越多[1]。湍流是自然界和工程技術(shù)中普遍存在的一種流動狀態(tài),降低湍流阻力對節(jié)約能源、提高效率以及消聲減振等方面的研究具有十分重大的學(xué)術(shù)和實際意義。目前關(guān)于減小湍流邊界層摩擦阻力的方法很多,且減阻效果均比較明顯,主要包括以下3類[2]:a)壁面添加其他物質(zhì)法:這類方法主要包括微氣泡減阻法和表面高分子涂層減阻法2種。b)邊界層控制法:對邊界層的控制就是設(shè)法延緩或避免邊界層的轉(zhuǎn)挾,保持更大的層流區(qū)域。具體控制邊界層的方法有彈性壁面法、邊界層加熱法以及壁面吹氣和吸氣法等。c)非光滑表面法:這種方法最常用的就是溝槽表面減阻法,這種方法是在表面加工條紋溝槽的方法來實現(xiàn)減阻目的。近年來人們對各種形狀的溝槽進行了減阻性能的研究,包括矩形、T型、U型、V型以及半圓型等多種形狀,都得到了一定程度上的減阻效果,大量研究結(jié)論已經(jīng)表明,溝槽減阻的確是一種十分有效的減阻方法[3~10]。

本文主要研究了一種環(huán)形管道中湍流流動的結(jié)構(gòu)形式。這種結(jié)構(gòu)流動形式在石油工程、熱交換器、渦輪機械,燃料電池、航空發(fā)動機和各種化工設(shè)備中廣泛存在,對這些設(shè)備的效率產(chǎn)生重大影響。近些年,對于環(huán)狀流的穩(wěn)定性問題和介于層流和湍流狀態(tài)的經(jīng)典問題有諸多研究[11~14]。其中,Moradi[15,16]做了大量關(guān)于環(huán)形管道中具有順流向溝槽的層流流動研究,分析了影響減阻效果的各種因素,包括波數(shù)、管道半徑、溝槽的高度,以及不同形狀的溝槽對減阻效果的影響,并分析了減阻系數(shù)的變化機理,所得結(jié)論對工程應(yīng)用具有重要的參考意義。

本文主要研究環(huán)形通道中順溝槽方向湍流狀態(tài)下減阻系數(shù)的變化規(guī)律,分析其內(nèi)在機理。首先將任意形狀的溝槽形狀表達為一階傅里葉模式,同時建立該模式的單位化無量綱模型用來計算不同狀態(tài)下的壓降損失。通過Fluent流場分析軟件研究不同的半徑、幅值和波數(shù)對減阻系數(shù)的影響。分別分析內(nèi)壁與外壁具有溝槽時對減阻系數(shù)的影響,并分析二者差別的內(nèi)在原因。同時研究了內(nèi)壁、外壁同時具有溝槽時的減阻特性。輸出不同情況下流體速度場隨幾何形狀的變化情況,研究流動分布的變化規(guī)律,分析影響減阻系數(shù)的內(nèi)在機理。

1 計算模型和方法

1.1 計算模型

假設(shè)在一個充分發(fā)展的環(huán)形通道中,流體的流動狀態(tài)為湍流狀態(tài)。流體被同心內(nèi)外光滑表面約束。內(nèi)表面和外表面的函數(shù)表達形式分別為inr和outr,如圖 1所示。

圖1 單位化模型Fig.1 Unit Model

內(nèi)、外表面的平均位置分別在 Rin= R 和Rout= R + 2 H,取H為特征長度,進行單位化處理。

內(nèi)、外表面均為任意形狀的溝槽,利用一階傅里葉模式對其進行近似。內(nèi)、外表面的函數(shù)表達式如式(1)所示:

式中inr,outr表示溝槽的幾何形狀;S代表溝槽的幅值;M為波的數(shù)目,MRα=,其中α為波數(shù),量綱化模型中槽的半高R為1;?為相位角。

流場分析軟件Fluent中的計算模型如圖2所示,為節(jié)省計算資源,建立周期對稱模型。網(wǎng)格由ANSYS ICEM CFD 14.5生成,利用塊技術(shù)生成高質(zhì)量的結(jié)構(gòu)化網(wǎng)格,且沿著壁面方向進行加密,以保證近壁面位置的計算精度,網(wǎng)格總數(shù)為40 000。

圖2 Fluent計算網(wǎng)格Fig.2 Fluent Calculation Grids

對網(wǎng)格的有效性進行校驗,輸出上下壁面 y+的大小, y+均小于1,說明Fluent的自帶函數(shù)可以較精確地計算近壁面的層流底層內(nèi)的問題。同時為了進一步提高計算精度,結(jié)果的收斂殘差設(shè)置為小于當?shù)貧埐?0-7,并且選取二階迎風格式來計算所有的平衡方程。

1.2 計算方法

不可壓黏性流體的RANS平衡方程如下:

本文利用商業(yè)通用流體計算軟件 ANSYS Fluent 14.5對流動過程進行模擬并求解 RANS方程,選取kω- SST作為湍流模型,同時選取Low-Re correction選項。

剪切應(yīng)力輸運模型是標準kω-模型的一種變形形式,結(jié)合了原有kω-模型對近壁面位置的精確計算以及標準kω-模型對遠離壁面位置的混合函數(shù),渦的黏性規(guī)則被適當修改以使湍流剪切應(yīng)力輸運效果精確計算。

Low-Re correction選項可以對湍流的層流邊界層直接計算而不必采取近似的壁面函數(shù),因此 Low-Re correction具有更高的精度。

假設(shè)存在一種理想的牛頓流體,設(shè)其密度1ρ=,動力黏度1μ=。在單位模型中溝槽的特征尺寸取一半槽的高度為1,因此雷諾數(shù)Re可由下列公式計算:

入口邊界條件為質(zhì)量流率入口,由下式給出:

計算入口面積時需要注意,必須保證有溝槽通道的面積與光滑通道面積相等,因此對于光滑通道半徑R1或 R1+ 2需根據(jù)面積相等原理進行計算修正。而對于內(nèi)表面或外表面單獨具有溝槽,或者雙壁面都具有溝槽的情況,通道面積 A的計算公式是不同的。這里只給出最終結(jié)果,如式(7)~(9)所示。

僅有內(nèi)壁具有溝槽情況下的通道面積:

僅有外壁具有溝槽情況下的通道面積:

兩側(cè)都具有溝槽情況下的通道面積:

對壓強及坐標進行單位化,并進行變化如下:

從而可得壓降的計算方法如式(12)所示:

摩擦系數(shù)由式(13)求得:

將式(12)帶入到式(13)中,可以得到摩擦系數(shù)計算公式如下:

式中 d p*/dz*可以直接由Fluent計算輸出。

減阻系數(shù)由式(15)給出,下標“0”代表無溝槽的光滑環(huán)形通道,下標“1”代表壁面有溝槽環(huán)形通道。

式中 下標“0”代表無溝槽的光滑環(huán)形通道;下標“1”代表壁面有溝槽環(huán)形通道。

2 計算模型的校驗

2.1 試驗對比

Moradi[15,16]做了大量關(guān)于環(huán)形管道中具有順流向溝槽的層流流動研究。通過公式推導(dǎo),獲得了環(huán)形管道具有順向溝槽層流流動的離散化方程,結(jié)合 Matlab數(shù)值計算,得到了減阻系數(shù)隨幾何形狀的變化規(guī)律。為了校核本文所建模型的合理性,選取相同狀態(tài)下的流動參數(shù)進行流動仿真分析。計算工況為1R=2、M=5和M=10、S=0.4,與Moradi文章中的1R=1、M=5和M=10、S=0.4相對應(yīng)。計算模型以及Moradi的計算結(jié)果如 圖3所示(本文的計算結(jié)果只是其曲線中的2個點)。

圖3 計算模型驗證Fig.3 The Validation of Calculation Model

續(xù)圖3

本文計算方法的計算結(jié)果分別為 f '5=0.1182和f '10=0.2927,與 Moradi的計算結(jié)果 ( f1/ f0)M=5=0.1179和(f1/ f0)M=10=0.2926非常接近,證明本文的計算模型設(shè)置合理。

2.2 與DNS結(jié)果對比

為了驗證上述RANS方程ω-k SST模型在模擬湍流流動時的精度,將相同的邊界條件及網(wǎng)格劃分應(yīng)用在直通道情況,與文獻[16]中的直接數(shù)值模擬 DNS計算結(jié)果進行對比。由于DNS計算結(jié)果將減阻效果增加表示為相同條件下質(zhì)量流率的增加,所以需要將RANS的計算結(jié)果輸出為質(zhì)量流率的變化。計算工況選取為幅值S=0.5,波數(shù)α分別從0.25變化到2.0的7個工況。計算結(jié)果對比如表1所示,兩者的變化規(guī)律如圖4所示。

表1 RANS計算結(jié)果與DNS計算結(jié)果對比Tab.1 Comparison of DNS and RANS Results for the Change in Discharge Flow Rate for a Grooved

圖4 DNS與RANS方法所得質(zhì)量流率變化對比Fig.4 Plot of DNS and RANS Results for the Change in Discharge Flow Rate for a Grooved Channel

由表1和圖4可以看出,RANS計算結(jié)果與DNS的計算結(jié)果變化趨勢一致,最大的絕對誤差僅為3.18%。在x軸的上方,即溝槽增加質(zhì)量流率的情況下,RANS的計算結(jié)果略大于DNS結(jié)果;在x軸的下方,即溝槽降低質(zhì)量流率的情況下,RANS的計算結(jié)果略小于DNS結(jié)果。總之,從7個點值上可以看出,RANS結(jié)果中質(zhì)量流率改變的絕對值均大于DNS的結(jié)果。說明在湍流減阻預(yù)測時,應(yīng)該注意RANS方法所得結(jié)果會略大。

DNS計算結(jié)果更接近實際情況,但其主要缺點是需要非常龐大的計算機內(nèi)存與機時耗費,其消耗大致與3Re成正比。因此在研究雷諾數(shù)變化范圍較大的湍流流動時(本文Re選取1000至30000范圍內(nèi)),本文研究的內(nèi)容,不同幾何形狀模型大概上千個,需要的計算資源較多。因此當精度在可以接受的范圍內(nèi)時,選擇RANS方法(上文已經(jīng)比較),可以節(jié)省大量計算資源和時間,且能夠得到合理的變化規(guī)律。

3 幾何參數(shù)對減阻系數(shù)變化規(guī)律的影響

3.1 內(nèi)壁具有溝槽通道減阻系數(shù)變化規(guī)律

首先計算外壁面為無溝槽光滑表面,內(nèi)壁面為溝槽表面的環(huán)形通道減阻系數(shù)。內(nèi)、外形面的表達式如式(16)所示:

因為需要一定的壓強降低才能保證驅(qū)動一定的質(zhì)量流率,因此定義沿通道的壓強降低即為阻力損失。而通道內(nèi)徑1R、溝槽的幅值S、波數(shù)RM/=α的變化都會影響壓強的損失,下面分別針對不同的半徑、幅值、波數(shù)進行計算,研究減阻系數(shù)隨其變化規(guī)律。

Re=2648,幅值S=0.75,α由0.1變化到1.0,M由 1變化到7,計算不同狀態(tài)下的減阻系數(shù) 'f,計算結(jié)果如圖5所示。從圖5中可以看出,波數(shù)M相同的情況下,隨通道半徑 R1的增加,即α逐漸降低情況下,減阻系數(shù)逐漸降低,減阻效果增加, 'f趨近一個極限值。以1=M為例,1.0→α?xí)r,環(huán)形通道的減阻系數(shù)逐漸趨近于-0.070 49。與直溝槽渠道的結(jié)果-0.070 43非常接近,如圖5b所示。其他M數(shù)也可以得到相同結(jié)論,說明當R大到一定程度時,便可以忽略環(huán)形通道曲率的影響,直接借用直溝槽通道的減阻系數(shù)結(jié)果。另外從圖5中也可以看出,半徑R不變情況下。M數(shù)增加不利于減阻效果的降低,即M增加將會導(dǎo)致α增加,會降低減阻效果,這一結(jié)論與直溝槽渠道結(jié)論一致。

圖5 減阻系數(shù)隨幾何形狀的變化規(guī)律Fig.5 Drag Reduction Factor Changing with the Geometry

3.2 外壁具有溝槽通道減阻系數(shù)變化規(guī)律

研究內(nèi)表面為無溝槽光滑表面,外表面為溝槽表面的環(huán)形通道減阻系數(shù)。內(nèi)、外形面的表達式如式(17)所示:

Re=2648,幅值S=0.75,α由0.1變化到1.0,M選取3個值分別為1,4,7,計算不同狀態(tài)下的減阻系數(shù) 'f,將結(jié)果與內(nèi)壁具有溝槽情況對比,如圖6所示。

圖6 減阻系數(shù)隨幾何形狀的變化規(guī)律(外壁溝槽)Fig.6 Drag Reduction Factor Changing with the Geometry實線—內(nèi)壁具有溝槽情況;虛線—外壁具有溝槽情況

從圖6可以看出,在具有相同半徑和波數(shù)的情況下,外壁具有溝槽情況下的減阻系數(shù)與內(nèi)壁具有溝槽情況相差較小,變化趨勢相同。隨著半徑 R1的增加,差別逐漸減小,且均趨近一個臨界值。在α較小時,外壁具有溝槽情況下的減阻系數(shù)略小于內(nèi)壁具有溝槽情況。分析認為,相同參數(shù)條件下,外壁溝槽情況的波數(shù) α = M /(R1+2)略小于內(nèi)壁溝槽情況下的波數(shù)α= M/R1,因此減阻系數(shù)略小,減阻效果略高。

3.3 內(nèi)、外壁均具有溝槽通道減阻系數(shù)變化規(guī)律

內(nèi)、外形面的表達式如式(18)所示:

對于不同相位角?對減阻系數(shù)的影響進行分析,層流狀態(tài)下,當?=π時可以獲得最大的減阻系數(shù),即“converging-diverging”位置[3,4,16]。同樣,對于湍流狀態(tài),同樣是當?=π時可以得到最大的減阻系數(shù)。后續(xù)的數(shù)值計算,對于內(nèi)外壁面都具有溝槽情況,都選取converging-diverging形式進行計算分析。

Re=2648,幅值S=0.75,α由0.1變化到1.0,M選取3個值分別為1,4,7,計算不同狀態(tài)下的減阻系數(shù) f ',將結(jié)果與內(nèi)、外壁具有溝槽情況對比,如圖 7所示。從圖7可以看出,雙面具有溝槽情況,可以獲得更大的減阻系數(shù),最大可以達到單壁面具有溝槽通道的3.25倍。同樣,同一M數(shù)條件下,隨通道半徑R1的增加,即α逐漸降低情況下,減阻系數(shù)逐漸降低,減阻效果增加, f '趨近一個極限值。以M=1為例,α→ 0.1時,環(huán)形通道的減阻系數(shù)逐漸趨近于-0.2532。查看前期直溝槽渠道的結(jié)果為-0.2 523,如圖7b所示。其他M數(shù)也可以得到相同結(jié)論,說明當半徑大到一定程度時,便可以直接借用直溝槽通道的減阻系數(shù)結(jié)果。

圖7 減阻系數(shù)隨幾何形狀的變化規(guī)律Fig.7 Drag Reduction Factor Changing with the Geometry

3.4 溝槽幅值對減阻系數(shù)的影響

為了研究幅值對減阻系數(shù)的影響,選取波數(shù)α=0.5的情況(M=4、1R=8),計算不同幅值S下的減阻系數(shù),具體如圖8所示。

圖8 幅值對減阻系數(shù)影響Fig. 8 The Effect of Amplitude on Drag Reduction

從圖 8可以看出,隨幅值S的增加,減阻系數(shù)降低,通道的減阻效果增加。對于內(nèi)壁面溝槽情況,f'S=1.5=-0.1054和 f 'S=0.25=-0.0151,幅值1.5可以達到幅值0.25的7倍。對于雙面溝槽情況,f 'S=0.75=-0.1694和 f 'S=0.25=-0.03371,幅值0.75可以達到幅值0.25的5倍。可見,在減阻范圍內(nèi),增加溝槽的幅值,更有利于增加通道的減阻效果。

3.5 不同形狀溝槽對減阻系數(shù)的影響

本文對其他幾種類似形式的溝槽進行了同樣的計算分析,具體的幾種形狀如圖9所示。波數(shù)(M=4)、幅值(S=0.75)、半徑及邊界條件均相同,計算結(jié)果如圖10所示。

圖9 不同形狀溝槽示意Fig. 9 Different Shape Groove on Annuli

圖10 不同形狀溝槽下的結(jié)果Fig.10 The Results of Different Annuli with Different Shape Grooves

從圖10中可以看出,單側(cè)溝槽減阻能力最弱,圖中單側(cè)三角形溝槽曲線幾乎接近X=0的軸線,隨半徑增加,減阻能力增加并不明顯;圓形等深溝槽與正弦曲線溝槽的減阻能力相當,差別較小;三角形等深溝槽的減阻能力較正弦曲線溝槽小近1倍;矩形溝槽的減阻系數(shù)變化幅度較大,當半徑較小時(約小于8),并無減阻能力,隨半徑增大,其減阻能力迅速增加,變化斜率大于圓形等深溝槽和正弦曲線溝槽。從方便制造加工角度講,矩形溝槽更容易實現(xiàn),但是選擇這種形式時,通道的半徑應(yīng)該足夠大(即保證波數(shù)α較小)。

4 減阻機理討論

為了分析幾何形狀影響減阻效果的內(nèi)在機理,將不同情況下,通道內(nèi)部的流體流速輸出,如圖11所示,波數(shù)M=4,通道內(nèi)徑1R=8,幅值S=0.75。

從圖11中可以看出,相比于光滑通道,由于溝槽的存在,流體的流動分布發(fā)生變化,溝槽中心速度較大,流體大部分質(zhì)量流集中在溝槽中心位置,壁面位置的流速減小。對于湍流流動,近壁面位置有薄薄的附面層,剪切應(yīng)力與速度變化率密切相關(guān),可以通過式(19)求出。附面層的速度梯度變化,必然導(dǎo)致壁面剪應(yīng)力變化。根據(jù)平衡原理,剪應(yīng)力變化將導(dǎo)致壓降變化,從而影響減阻系數(shù)變化。

式中 τ為壁面黏性層的剪切應(yīng)力;μ為黏性系數(shù);du/dy為沿y方向(與流體速度方向垂直)的速度梯度。

圖11 通道的速度分布Fig.11 The Velocity Distribution of Annuli

續(xù)圖11

如圖11b、11c所示,對于相同情況下,內(nèi)、外壁面分別具有溝槽時的速度分布相似。同樣根據(jù)式(19)可知,壁面的剪應(yīng)力分布相似,從而壓降相似導(dǎo)致減阻系數(shù)相差不大。對于內(nèi)外壁面均有溝槽,且位置為“分離-聚合”的形式時,溝槽中心的速度明顯大于一側(cè)具有溝槽的情況,如圖 11d所示,中心速度增加,流體質(zhì)量流集中在中心位置,壁面附近流過的流體質(zhì)量減小,黏性摩擦降低,剪應(yīng)力降低將導(dǎo)致壓降降低,減阻系數(shù)降低,溝槽的減阻效果增加。

5 結(jié) 論

a)對于環(huán)形通道,適當?shù)脑黾訙喜劭梢詫崿F(xiàn)減阻效果。相同波數(shù)情況下,隨通道半徑R1的增加,即α逐漸降低情況下,減阻系數(shù)逐漸降低,減阻效果增加, 'f趨近一個極限值。當R1大到一定程度時,便可以忽略環(huán)形通道曲率的影響,直接借用直溝槽通道的減阻系數(shù)結(jié)果。

b)具有相同半徑 R1和波數(shù)M的情況下,外壁具有溝槽情況下的減阻系數(shù)與內(nèi)壁具有溝槽情況相差較小,變化趨勢相同。隨著半徑R1的增加,差別逐漸減小,且均趨近一個臨界值。在α較小時,外壁具有溝槽情況下的減阻系數(shù)略小于內(nèi)壁具有溝槽情況。分析認為,相同參數(shù)條件下,外壁溝槽情況的波數(shù)α= M /(R1+2)略小于內(nèi)壁溝槽情況下的波數(shù) α = M/R1,因此減阻系數(shù)略小,減阻效果略高。

c)雙面具有溝槽情況,“分離-聚合”形式可以獲得更大的減阻系數(shù),最大可以達到單壁面具有溝槽通道的3.25倍。同樣,同一M數(shù)條件下,隨通道半徑R的增加,即α逐漸降低情況下,減阻系數(shù)逐漸降低,減阻效果增加, f '趨近一個極限值。當R大到一定程度時,便可以直接借用直溝槽通道的減阻系數(shù)結(jié)果。

d)相比于光滑通道,溝槽將導(dǎo)致流體的流動分布變化,溝槽中心速度較大,流體大部分質(zhì)量流集中在溝槽中心位置,壁面位置流過的流體質(zhì)量減小,黏性摩擦降低。根據(jù)平衡原理,黏性摩擦剪應(yīng)力降低將導(dǎo)致壓降降低,減阻系數(shù)降低,通道的減阻效果增加。

e)對于其他形式的減阻溝槽,單側(cè)溝槽減阻能力最弱。圓形等深溝槽與正弦曲線溝槽的減阻能力相當。三角形等深溝槽的減阻能力較正弦曲線溝槽小近1倍。矩形溝槽的減阻系數(shù)變化幅度較大,隨半徑增大,其減阻能力迅速增加。選擇矩形溝槽時,通道的半徑應(yīng)該足夠大(即保證波數(shù)α較小)。

猜你喜歡
溝槽
柔性表皮與微溝槽耦合作用減阻效果數(shù)值模擬
基于數(shù)值模擬的2種條紋溝槽減阻特性對比分析
開槽施工鋼筋混凝土管道的臨界溝槽寬度
一種具有多形式鋼片結(jié)構(gòu)的四季胎
溝槽管旋壓拉拔一體式成形裝置的設(shè)計*
機械制造(2021年8期)2021-08-23 10:12:44
基于曲軸溝槽圓角滾壓系統(tǒng)的刀具壽命提升
一種低噪音的全路況輪胎
二層溝槽織構(gòu)對機床導(dǎo)軌表面潤滑特性的影響
簡析電子線路板微溝槽脈沖鍍銅填充工藝
電子制作(2018年14期)2018-08-21 01:38:30
溝槽爆破參數(shù)優(yōu)化及成本分析
主站蜘蛛池模板: 色综合久久88| 为你提供最新久久精品久久综合| 欧美不卡视频一区发布| m男亚洲一区中文字幕| 久久精品这里只有国产中文精品| 婷婷亚洲视频| 午夜视频免费试看| 色综合a怡红院怡红院首页| 久久久噜噜噜久久中文字幕色伊伊| 亚洲二区视频| 国产成人免费高清AⅤ| 亚洲国产欧美国产综合久久| 中文字幕乱码中文乱码51精品| 91精品aⅴ无码中文字字幕蜜桃| 男女猛烈无遮挡午夜视频| 成人免费黄色小视频| 日韩精品专区免费无码aⅴ| 亚洲精品欧美日本中文字幕| 成人免费视频一区二区三区| 91亚洲精品第一| 国产精品页| 久久久无码人妻精品无码| 亚洲国产av无码综合原创国产| 精品人妻无码区在线视频| 综合天天色| 欧美国产日韩在线播放| 久久综合干| 亚洲国产亚综合在线区| 久久精品66| 乱人伦中文视频在线观看免费| 亚洲精品桃花岛av在线| 中日韩一区二区三区中文免费视频| 成人福利一区二区视频在线| 国产精品午夜电影| 免费99精品国产自在现线| 国产资源免费观看| 谁有在线观看日韩亚洲最新视频| 国产激情影院| 亚洲欧美成人在线视频| 亚洲天堂网视频| 麻豆国产在线观看一区二区| 日韩国产欧美精品在线| 国产精品冒白浆免费视频| 欧美亚洲国产视频| 久久精品人妻中文视频| 亚洲男人的天堂久久香蕉| 99久视频| 欧美伊人色综合久久天天| 91精品专区国产盗摄| 日韩精品高清自在线| a毛片基地免费大全| 久久视精品| 黄色网在线| 久久天天躁狠狠躁夜夜2020一| 亚洲精品无码日韩国产不卡| 丰满人妻中出白浆| 国产精品太粉嫩高中在线观看| 國產尤物AV尤物在線觀看| 久久一本日韩精品中文字幕屁孩| 国产91丝袜在线播放动漫| 999精品视频在线| 亚洲无码精彩视频在线观看| 亚洲美女久久| 欧美a级完整在线观看| 99久久精品免费看国产电影| 99久久精品视香蕉蕉| 国产成人乱无码视频| 欧美三级日韩三级| 欧美一区精品| 亚洲三级影院| 找国产毛片看| 色哟哟国产精品一区二区| 四虎永久免费地址| 亚洲 欧美 中文 AⅤ在线视频| 伊人蕉久影院| 亚洲美女一区| 免费jizz在线播放| 国产一级毛片yw| 青青青国产视频| 国产在线拍偷自揄拍精品| 国产福利免费观看| 18禁高潮出水呻吟娇喘蜜芽|