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

        ?

        矩形通道內脈動湍流流動特性實驗研究

        2017-09-25 07:55:46劉宇生譚思超高璞珍
        核安全 2017年2期
        關鍵詞:慣性力湍流脈動

        劉宇生,許 超,*,譚思超,胡 健,高璞珍

        (1.環(huán)境保護部核與輻射安全中心, 北京 100082;2. 哈爾濱工程大學核安全與仿真技術國防重點學科實驗室, 黑龍江 哈爾濱 150001)

        矩形通道內脈動湍流流動特性實驗研究

        劉宇生1,許 超1,*,譚思超2,胡 健1,高璞珍2

        (1.環(huán)境保護部核與輻射安全中心, 北京 100082;2. 哈爾濱工程大學核安全與仿真技術國防重點學科實驗室, 黑龍江 哈爾濱 150001)

        在流量脈動條件下,本文對矩形通道內的湍流流動特性進行了實驗研究,通過理論分析,得到了影響脈動湍流的主要作用力和關鍵的無量綱數,分析了脈動周期、相對振幅等因素對流量與壓降的相位差、壓降-流量曲線、時均摩阻系數的影響,并與穩(wěn)定流動狀態(tài)下的流動特性進行了對比。實驗結果表明:對于低頻脈動湍流,壓降驅動力、摩擦阻力和流體自身慣性力是影響流動特性的主要作用力;脈動湍流中,流量的變化滯后于壓降變化,存在相位差;由于流體的慣性作用,脈動周期越小,流量脈動的幅值越小;速度相對振幅增大并超過臨界值時,時均摩阻系數會顯著增大。

        脈動湍流; 矩形通道; 阻力特性; 熱工水力; 相似準則

        脈動流是流量周期性波動但時均流量不為零的一種流動現象,屬于非穩(wěn)定流動。脈動流在核動力系統(tǒng)廣泛存在:如采用非能動安全設計的核電廠,事故條件下在非能動系統(tǒng)投入運行的過程中,冷卻劑一直處于非穩(wěn)定流動狀態(tài),會出現流量脈動現象[1];受海洋條件的影響,船舶核動力系統(tǒng)內的冷卻劑也會出現周期性的流量脈動[2]。Elsayed A.M[3]、M.A. Habib[4]、賈輝[5]、王暢[6]等人對圓管內的脈動流進行了研究,其研究結果表明:脈動流的流動換熱特性與穩(wěn)定流動狀態(tài)下的流動換熱特性存在顯著差異,脈動流瞬態(tài)流動換熱過程的研究對核動力系統(tǒng)設計和核安全評價具有重要意義。

        近年來,核動力系統(tǒng)日漸向小型化、模塊化方向發(fā)展[7-9],矩形流動通道由于換熱面積大且結構緊湊,已在研究堆中得到了應用[10]。由于幾何結構特殊,研究人員對矩形通道內的流場分布、阻力特性、換熱特性進行了大量研究[11-15],但所得結果都是針對穩(wěn)定流動狀態(tài),對矩形通道內的脈動流動則研究很少,無法為事故條件下核動力系統(tǒng)瞬態(tài)過程的安全分析提供充分的試驗結果,因此有必要開展這方面的研究,為非穩(wěn)定流動狀態(tài)下的熱工水力設計和安全分析提供依據和參考。

        本文使用試驗手段研究了脈動湍流狀態(tài)下矩形通道內的流動阻力特性,通過對比分析,獲得了脈動周期、脈動相對振幅等參數對流動特性的影響規(guī)律。

        1 實驗裝置及驗證

        在本研究中,對矩形通道內的脈動湍流進行了簡化,認為流體的脈動速度隨時間呈正弦波動變化,脈動速度u可寫為:

        (1)

        1.1 實驗裝置

        實驗回路簡圖和試驗段如圖1、圖2所示。試驗過程中,水箱中的流體在循環(huán)水泵作用下依次流經循環(huán)回路中的控制閥組、流量計和試驗段,最終回到水箱完成循環(huán)。回路中的循環(huán)水泵由變頻器控制,通過變頻器控制泵的轉數,使回路流體產生流量脈動。流體的溫度、流量和壓降分別由溫度計、電磁流量計和電容式壓差傳感器測量獲得。

        圖1 實驗回路簡圖Fig.1 Schematic of experimental apparatus

        圖2 實驗段簡圖Fig.2 Schematic of test section

        本實驗中使用的試驗段為內壁光滑的窄矩形通道,通道的截面尺寸為40.38 mm×3.1 mm,試驗段由透明的有機玻璃管加工而成。試驗段布置有取壓孔,考慮到進出口效應的影響,取壓孔布置在了試驗段上流動充分發(fā)展的位置。

        1.2 實驗參數范圍及儀表精度

        實驗流體采用去離子水,實驗環(huán)境條件和參數范圍如下:實驗回路壓力為0.1 MPa,試驗段壓降變化范圍為0-50 kPa,雷諾數Re為0-15000,流體溫度為14.0 ℃-15.0 ℃,脈動周期為5 s-60 s。

        本實驗測量的參數范圍較大,包含了層流和湍流兩種流動狀態(tài),在測量儀表上使用了測量范圍和精度等級不同的電磁流量計和電容式壓差傳感器,見表1。

        表1 儀表測量范圍和精度等級Table 1 Measurement range and accuracy grade of instruments

        1.3 實驗系統(tǒng)驗證

        在進行脈動湍流實驗之前,利用試驗裝置對不同流量下矩形通道內的穩(wěn)定流動過程進行了試驗,得到相應的l-Re關系,并與現有計算值進行了對比,如圖3所示。圖中層流區(qū)的理論值根據凱斯公式[16]計算得到,其公式為:

        λRe=96(1-1.3553α+1.9467α2-1.7012α3+0.9564α4-0.2537α5)

        (2)

        式中:a為矩形通道橫截面寬高比,本試驗段的寬高比為0.77,于是Re=87。

        湍流區(qū)的對比值根據SADATOMI[17]和勃拉休斯[18]公式計算得到。對比實驗值和計算值可知,二者符合良好,試驗系統(tǒng)穩(wěn)定可靠。

        圖3 穩(wěn)定流動狀態(tài)試驗結果分析Fig.3 Analysis of test results under steady flow state

        2 脈動湍流過程理論分析

        將公式(1)定義的脈動速度對時間求導,得到脈動加速度的表達式為:

        (3)

        將兩端測壓截面內的流體作為控制體,假設水平圓管長為l,A為管道面積,τo為單位面積上摩擦力,Ph為控制體周界長,則管道內的微分方程為:

        (4)

        對微分方程積分后,可得控制體的積分守恒方程為:

        (5)

        利用穩(wěn)定流動的參數值對上述方程進行無量綱化,可得到:

        (6)

        其中:

        (7)

        (8)

        τrs為居留時間,表征了流體質點在控制體內流過所需的時間。

        ΠEu為歐拉數,是驅動壓降與流體動能的比值,表征了壓降作用力與流體慣性力的相對大小。此外,歐拉數還可寫為:

        (9)

        方程(9)中,τp表征了驅動壓降對流體作用的時間,即歐拉數表征了居留時間與壓降特征時間的相對大小。

        ΠF為摩擦數,是沿程摩擦與流體動能的比值,表征了摩擦力與流體慣性力的相對大小。此外,摩擦數還可寫為:

        (10)

        方程(10)中,τf為摩擦力對流體作用的時間,即摩擦數表征了居留時間與摩擦力特征時間的相對大小。

        ΠFr為Froude數,是加速度慣性力與流體動能的比值,表征了加速慣性力與流體慣性力的相對大小。此外,Froude數可以表示為Womersley數與雷諾數的比值,體現了瞬態(tài)慣性力與對流慣性力的相對大?。?/p>

        (11)

        Froude數還可寫為:

        (12)

        方程(12)中,τa表征了加速度對流體作用的時間,即Froude數體現了居留時間與加速度特征時間的相對大小。

        上述無量綱值分別反映了驅動壓力、摩擦阻力和加速度慣性力對流體作用時間與居留時間的相對大小,同時反映了三種作用與流體慣性力的相對大小,因此,無量綱數的數值表明了不同作用力的強弱和作用時間的長短。根據方程(8),計算典型試驗工況的無量綱數值,見表2。

        表2 典型試驗工況無量綱數Table 2 Dimensionless numbers for typical test conditions

        由表2可知:對于本試驗研究的低頻脈動工況,Froude數較小,歐拉數與摩擦數大小接近,因此壓降驅動力和摩擦阻力是低頻脈動湍流中的兩個主要作用力,而流體加速慣性力則可以忽略。又因壓降驅動力由泵提供,屬于控制參數,所以反映摩擦阻力的摩擦壓降和摩阻系數是本文的研究重點。

        3 脈動湍流的相位差及壓降-流量曲線

        3.1 流量與壓降的相位差

        實驗系統(tǒng)內的流體流量處于脈動瞬態(tài)時,試驗段內流體的流量和壓降均隨時間呈周期性地變化,如圖4所示??芍?,脈動湍流的流量隨時間的變化滯后于壓降隨時間的變化,流量與壓降間存在明顯的相位差,這與矩形通道脈動層流的研究結論一致[19]。

        圖4 典型試驗現象Fig.4 Typical experimental phenomena

        已有的研究表明,脈動層流中流量與壓降的相位差主要受到脈動周期的影響[19]。圖5中繪制了脈動湍流工況下流量-壓降間相位差隨脈動周期變化的曲線??芍S著脈動周期的減小,相位差在0°-90°之間增加,當脈動周期超過60s時,相位差趨近于0°;當脈動周期小于10s時,相位差趨近于90°。圖5同時表明:不同壓降相對振幅對流量與壓降的相位差的影響規(guī)律不明顯,需要進一步的實驗研究來確定二者之間的關系。

        圖5 不同壓降相對振幅下相位差隨周期的變化Fig.5 Changing curves of phase difference with pulsating period under different relative pressure drop amplitude conditions

        3.2 不同脈動周期下的壓降-流量曲線

        在矩形通道內的湍流脈動流中,試驗段的流量和壓降隨時間變化呈周期性脈動,因此,需要對壓降、流量的實時變化情況進行研究,即研究壓降-流量曲線。對脈動層流的研究表明,脈動過程中的速度幅值主要受壓降相對振幅和脈動周期兩個因素的影響[20]。試驗中保持壓降相對振幅不變,可得到不同脈動周期下的壓降-流量曲線,如圖6所示,圖中箭頭標示了單個周期內流量脈動的變化規(guī)律??芍?,隨著脈動周期的減小,流量與壓降的相位差不斷增大,二者的夾角沿逆時針方向不斷增大,導致環(huán)形曲線從傾斜的橢圓逐漸演變?yōu)樨Q直的圓形。

        圖6中還繪出了穩(wěn)定流動狀態(tài)下的摩阻壓降,通過比較可知,當出現流量脈動時,在每個周期的過程中,流體在加速狀態(tài)下的摩阻壓降要高于穩(wěn)定流動狀態(tài)下的摩阻壓降,而處于相同的減速狀態(tài)時,其摩阻壓降又低于穩(wěn)定流動狀態(tài)摩阻壓降。在流量最小值和最大值附近,流量脈動的摩擦壓降與穩(wěn)定流動的摩擦壓降相等;結合圖4可知,每個脈動周期內,在流量達到最小值和最大值時,流體的加速度為零,其余時間內,流體一直處于加速流動或減速流動的狀態(tài)。

        該現象表明瞬態(tài)過程中的流體阻力特性與穩(wěn)定狀態(tài)下阻力特性存在顯著不同。根據本文理論分析的結果,在穩(wěn)定流動狀態(tài)下,摩擦力與壓降驅動力相等,摩阻壓降體現了流體粘性導致的能量耗散,因流速不發(fā)生改變,流體的慣性力不會體現;而當流速脈動時,瞬態(tài)過程需要考慮壓降驅動力、摩擦力與流體慣性力之間的相互作用,脈動周期減小,流體加速度增大,都會使流體的慣性作用增強,進而導致流體微團響應驅動壓降變化的時間增長,即導致流量-壓降間相位差增加,最終壓降-流量環(huán)形曲線的形狀也會接近圓形。

        從流量變化的幅值來看,脈動周期的減小,脈動幅值也會減小,對于脈動周期為5 s的工況,這種現象最為明顯。這是因為隨著脈動周期變小,瞬態(tài)變化中流體體現出的慣性會加強,進而導致流量對壓降梯度的響應出現延遲和滯后;當該延遲時間接近或小于脈動周期后,流量的慣性作用進一步變大,流量對壓降梯度的響應更加緩慢,最終導致流量能達到的幅值很低,只能在時均值附近小幅脈動。

        圖6 不同周期下的壓降-流量曲線Fig.6 Curves of pressure drop-flow rate with different pulsating periods

        綜合圖5與圖6可知:脈動紊流中,流量-壓降間相位差對壓降-流量曲線的形狀變化具有重要影響,在小于10s的脈動周期下,流量與壓降之間的相位差接近90°,壓降-流量曲線的形狀接近圓形;在大于60s的當脈動周期下,流量-壓降間相位差接近0°,壓降-流量曲線的形狀為扁長形,與穩(wěn)定流動狀態(tài)下壓降-流量曲線接近。

        4 脈動湍流的時均摩阻系數分析

        對圓管內脈動流研究表明,脈動流的摩阻系數可采用周期平均的方法來進行研究,其主要的影響因素有:雷諾數,脈動周期和速度相對振幅[21]。

        圖7為不同脈動周期條件下,時均摩阻系數隨時均雷諾數變化曲線??芍S著時均雷諾數的增加,時均摩阻系數呈緩慢下降趨勢,該變化規(guī)律與穩(wěn)定流動狀態(tài)下l-Re曲線的變化規(guī)律相同。這是因為脈動湍流中時均雷諾數與壁面邊界層厚度緊密相關,時均雷諾數增大,壁面邊界層厚度減薄,壁面邊界層內的速度梯度也會改變,進而使得摩阻系數發(fā)生變化。圖7中還匯出了勃拉休斯公式和SADATOMI公式的計算值,可知脈動湍流時均摩阻系數與穩(wěn)定流動變化趨勢一致,這表明可采用勃拉休斯公式的形式對其進行公式擬合。

        圖7 時均摩阻系數隨時均雷諾數的變化Fig.7 Variation of period-averaged friction coefficient with period-averaged Reynolds number

        圖8為不同脈動周期工況下時均摩阻系數隨速度相對振幅變化的曲線,可知速度相對振幅增加,時均摩阻系數會隨之增大。

        圖8 不同工況下的時均摩阻系數Fig.8 The period-averaged friction coefficient under different cases

        與穩(wěn)定流動的時均摩阻系數比較后可知,不同脈動周期下,均存在一個速度相對振幅臨界值,小于該值的情況下,流量脈動現象對時均摩阻系數的影響較為微弱,時均摩阻系數與穩(wěn)定流動摩阻系數差別較?。淮笥谠撍俣认鄬φ穹岛?,流量脈動現象的影響變得顯著,隨著脈動流速度變化劇烈程度的增加,時均摩阻系數迅速增大,與穩(wěn)定流動狀態(tài)下的摩阻系數存在明顯不同。此外,脈動周期變短,速度相對振幅的臨界值會變小。出現上述現象的原因主要有兩個:速度相對振幅增大,表示流體速度變化的幅值增加,在近壁區(qū)邊界層內,流體微層之間的速度梯度也會隨之增大,微觀上體現為流體微層間的摩擦切應力增加,宏觀上體現為時均摩阻系數增大;脈動周期變短,會加劇流體微團之間的橫向擾動,導致湍流切應力增大,摩阻壓降增加,時均摩阻系數增大。

        結合圖7和圖8的曲線特點,通過雷諾數、無量綱周期數[22]w′和速度相對振幅三個參數擬合得到矩形通道脈動湍流條件下時均摩阻系數的經驗公式,如式(13)所示:

        (13)

        式中,C和t為系數,與無量綱周期有關,lst為對應穩(wěn)定流動狀態(tài)的摩阻系數,可通過時均雷諾數計算得到。式(14)、(15)、(16)分別給出了C、t、lst的計算方法。

        (14)

        (15)

        λst=0.23Re-0.23

        (16)

        圖8為本文所擬合的經驗公式預測值與實驗值對比的結果,可知預測值與實驗值符合良好,二者的偏差在±5%以內,可滿足非穩(wěn)定流動狀態(tài)下熱工水力設計和安全分析預測計算的精度要求。

        圖9 擬合公式計算值與實驗值的對比Fig.9 Performance of correlation for predicting the magnitude of friction coefficient

        5 結論

        通過對矩形通道內脈動湍流的流動特性進行實驗研究,本文得到如下結論:

        (1)矩形通道內的湍流出現流量脈動時,流量隨時間的變化滯后于壓降隨時間的變化,二者之間存在明顯的相位差;當脈動周期減小時,該相位差增大,脈動周期超過60 s時,相位差趨近于0°,周期小于10 s時,相位差趨近于90°。

        (2)影響流量脈動的主要作用力為壓降驅動力、摩擦力和流體慣性力。流體慣性導致的相位差是壓降-流量環(huán)形曲線呈現橢圓形和圓形等不同形態(tài)的主要原因;隨著脈動周期的減小,脈動速度變化的幅值減小。

        (3)速度相對振幅對時均摩阻系數的影響存在臨界值現象,當流量脈動小于速度相對振幅臨界值時,流量脈動與穩(wěn)定流動狀態(tài)的摩阻系數相接近;當大于該臨界值時,時均摩阻系數迅速增加,與穩(wěn)定流動狀態(tài)的摩阻系數明顯不同。時均摩阻系數擬合關系式預測值與實驗結果符合良好。

        [1] 袁添鴻,于雷,王川. 全廠斷電事故下AP1000 非能動余熱排出系統(tǒng)分析[J]. 原子能科學技術,2010, 44(zl):248-252.

        [2] 譚思超,龐鳳閣. 搖擺運動引起的波動與自然循環(huán)密度波型脈動的疊加[J]. 核動力工程, 2005, 26(2):140-143.

        [3] Elsayed A.M. Elshafei, M. Safwat Mohamed, H. Mansour, M. Sakr. Experimental study of heat transfer in pulsating turbulent flow in a pipe[J]. International Journal of Heat and Fluid Flow, 2008, 29: 1029-1038.

        [4] M.A. Habib, A.M. Attya, S.A.M. Said, A.I. Eid, A.Z. Aly. Heat transfer characteristics and Nusselt number correlation of turbulent pulsating pipe air flows[J]. Heat and Mass Transfer, 2004, 40: 307-318.

        [5] 賈輝, 譚思超, 高璞珍,等. 簡諧脈動流中極點摩擦壓降特性實驗研究[J]. 核動力工程, 2011, 32(3):102-105.

        [6] 王暢, 高璞珍, 許超. 層流脈動流動對流換熱數值分析[J]. 哈爾濱工程大學學報, 2011, 32(7):890-894.

        [7] 楊玨, 孫吉良, 楊偉國,等. 多用途小型堆 ACPR100概念設計[J]. 原子能科學技術, 2014, 48(10):1844-1849.

        [8] Carelli M D. The exciting journey of designing an advanced reactor[J]. Nuclear Engineering & Design, 2009, 239(5):880-887.

        [9] Reyes J N, Groome J, Woods B G, et al. Testing of the multi-application small light water reactor (MASLWR) passive safety systems[J]. Nuclear Engineering & Design, 2007, 237(18):1999-2005.[10] 周媛,王玉林. CARR堆芯熱組件自然循環(huán)條件下特性分析[J]. 原子能科學技術, 2015, 49(3): 433-439.

        [11] 蔣潔,郝英立,施明恒. 矩形微通道中流體流動阻力和換熱特性實驗研究[J]. 熱科學與技術, 2006, 5(3): 189-194.

        [12] 鄭慧凡,秦貴棉,范曉偉,等. 微通道內單相流流動特性的實驗研究進展[J]. 節(jié)能技術,2008,26(147):32-36.

        [13] 泰文波,程惠爾,牛祿,等.大高寬比微小寬度矩形通道內的水力特性實驗研究[J]. 熱科學與技術,2002,1(1):38-41.[14] 幸奠川, 閻昌琪, 王暢,等. 矩形通道高寬比對單相層流特性的影響[J]. 力學學報, 2013, 45(3):331-336.

        [15] Wang C, Gao P, Tan S, et al. Experimental study of friction and heat transfer characteristics in narrow rectangular channel [J]. Nuclear Engineering & Design, 2012, 250(3):646-655.

        [16] HARTNRTT J P, KOSTIC M.Heat transfer to Newtonian and non-Newtonian fluids in rectangular ducts[J]. Advances in Heat Transfer, 1989, 19: 247-356.

        [17] SADATOMI Y, SATO Y, SARUWATARI S. Two-phase flow in vertical noncircular channels[J]. Int J Multiphase Flow, 1982, 8: 641-655.

        [18] 孔瓏.工程流體力學[M]. 3版. 北京: 中國電力出版社, 2007: 108-113.

        [19] 劉宇生,高璞珍,譚思超,幸奠川. 矩形通道內脈動層流相位差實驗研究[J].原子能科學技術, 2013, 47(4):564-569.

        [20] 劉宇生,譚思超,高璞珍,袁其斌. 矩形通道內脈動層流流場特性理論研究[J]. 原子能科學技術, 2012, 46(11): 1318-1323.

        [21] 王暢. 流量波動條件下的流動與傳熱特性研究[D]. 哈爾濱:哈爾濱工程大學,2010.

        [22] Melda ?zdin? ?arpinliog lu, Mehmet Yas,ar Gündog du. A critical review on pulsatile pipe flow studies directing towards future research topics[J]. Flow Measurement and Instrumentation, 2001, 12: 163-174.

        ExperimentalStudyonFlowCharacteristicsofPulsatingTurbulentFlowinRectangularChannel

        LIU Yusheng1, XU Chao1,*, TAN Sichao2, HU Jian1,GAO Puzhen2

        (1.Nuclear and Radiation Safety Center, MEP, Beijing 100082, China;2.Key Discipline Laboratory of Nuclear Safety and Simulation Technology, Harbin Engineering University, Harbin 150001, China)

        Characteristics of pulsating turbulent flow in a rectangular channel were investigated experimentally. The dimensionless numbers for acting force in pulsating turbulent flow are obtained by theoretical analysis. The effects of pulsating period and relative amplitude on flow field in rectangular channel were analyzed. Characteristics of pulsating turbulent flow were compared with that of steady flow. The results show that pressure driving force, friction force and inertial force are the main acting forces affecting flow characteristics of low-frequency pulsating flow. Flow rate change lags behind the variation in pressure drop, a phase difference exists. When pulsating period becomes small, the pulsating amplitude of flow rate will decrease as the inertia of fluid. Besides, the period-averaged friction coefficient will increase rapidly when relative velocity amplitude exceeds the critical value.

        pulsating turbulent flow; rectangular channel; resistance characteristic; thermal hydraulics; similarity criteria

        TL334

        :A

        1672- 5360(2017)02- 0056- 07

        2016- 12- 27

        2017- 01- 30

        國家科技重大專項資助項目,項目編號:2015ZX06002007

        劉宇生(1986—),男,河北唐山人,工程師,碩士,核能科學與工程專業(yè),現主要從事反應堆熱工水力領域的研究工作

        猜你喜歡
        慣性力湍流脈動
        新學期,如何“脈動回來”?
        家教世界(2023年25期)2023-10-09 02:11:56
        RBI在超期服役脈動真空滅菌器定檢中的應用
        新型模切機滑塊運動特性及其慣性力特性研究
        從“慣性力”到“洛希極限”
        重氣瞬時泄漏擴散的湍流模型驗證
        地球脈動(第一季)
        “青春期”湍流中的智慧引渡(三)
        “青春期”湍流中的智慧引渡(二)
        弱分層湍流輸運特性的統(tǒng)計分析
        地脈動在大震前的異常變化研究
        地震研究(2014年1期)2014-02-27 09:29:43
        国产高清黄色在线观看91| 亚洲av无码乱码在线观看富二代| 久久精品国产亚洲av网站| 亚洲一本到无码av中文字幕| 国产AV国片精品有毛| 精品一区二区三区中文字幕在线| 成人水蜜桃视频在线观看| 可以免费看亚洲av的网站| 中文字幕日韩人妻不卡一区| 亚洲av无码成人精品区在线观看| 动漫av纯肉无码av在线播放| 日韩伦理av一区二区三区| 国产精品日韩经典中文字幕| 欧美变态另类刺激| 一本大道久久香蕉成人网| 欧美黑人xxxx性高清版| 亚洲综合新区一区二区| 亚洲天堂二区三区三州| 女女女女女裸体处开bbb| 亚洲av无码成人精品区在线观看| 国产极品喷水视频| 久久国产精品国语对白| 国产三级a三级三级| 九九精品国产亚洲av日韩| 亚洲视频高清| 国产一区二区三区再现| 亚洲一区二区在线观看网址| 99久久久国产精品免费蜜臀| 亚洲成人观看| 手机在线免费看av网站| 穿着白丝啪啪的av网站| 亚洲人成77777在线播放网站| 亚洲女人被黑人巨大进入| 国产香蕉一区二区三区| 国产精品国产三级国产专区不| 俺去啦最新地址| 亚洲国产成人精品无码区99| 一区二区三无码| 亚洲女厕偷拍一区二区| 一本精品99久久精品77| 精品人体无码一区二区三区 |