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

球窩/球凸結構下的U型通道蒸汽冷卻性能

2015-03-07 00:34:04李森李亮杜長河
西安交通大學學報 2015年9期
關鍵詞:結構

李森,李亮,杜長河

?

球窩/球凸結構下的U型通道蒸汽冷卻性能

李森,李亮,杜長河

近些年來,可以起到強化傳熱效果而阻力損失較小的球窩/球凸結構[1]日益受到人們的關注,這種結構可以廣泛用于透平葉片、微通道、燃燒器等的冷卻中。Afansayev等指出,在平板上布置球窩時總體傳熱系數增加30%~40%,而阻力損失基本與光滑表面持平[1]。Moon等研究指出,隨通道高度的變化,傳熱增強效果和阻力損失幾乎保持常數,橫橢圓球窩的傳熱要優于圓形球窩[2-3]。Burgess等研究結果顯示,隨直徑高度比例的增加,傳熱系數線性增大,但阻力系數也迅速增大[4-5]。Ligrani等指出,球窩結構傳熱得到強化的原因在于球窩及下游流場內存在強烈的三維渦旋結構[6-7]。Rao等通過實驗研究了球窩、翅片以及二者混合布置時的傳熱及阻力特性,指出球窩/翅片的混合布置有著更好的傳熱特性和更低的阻力系數,其綜合傳熱性能優于單純的球窩結構[8]。Shen等研究發現,抽吸孔可以在不引起額外流動阻力的情況下強化傳熱過程[9]。

已有文獻中,大多數關于球窩/球凸結構傳熱特性的研究都是針對空氣進行的。與空氣相比,蒸汽的導熱率和比熱容更大,黏性系數更小,有著更優異的傳熱、載熱、氣動特性,帶肋通道中的研究表明[10],蒸汽的冷卻效能比空氣約高15%。

本文以蒸汽為工質,在不同雷諾數、來流溫比和出口壓力條件下數值分析了U型通道內球窩/球凸結構的流動與傳熱特性,比較了蒸汽和空氣流經球窩/球凸結構的傳熱和流動特性,以期為透平葉片冷卻通道的設計提供參考。

1 計算模型

為了研究球窩或球凸結構的影響,分別建立了2個幾何模型:模型A中,在U型通道[11]的流入段和流出段分別對稱布置了12個球窩結構,球窩表面直徑D為6.35 mm,U型通道橫截面為2D×2D,長度L為35D,球窩距進、出口的距離T為5D,球窩間距S為2D,球窩的相對深度δ/D為0.2;模型B中,只是將球窩換成球凸,球凸的相對高度δ/D為0.2,其余的結構及尺寸與模型A相同。U型通道幾何模型及網絡如圖1所示。

(a)U型通道結構

(b)球窩/球凸截面特征及局部網格圖1 U型通道幾何模型及局部球窩/球凸網格

計算中球窩/球凸表面采用O型網格,其他區域采用H型網格,并在壁面附近進行了加密處理,以保證計算域離開壁面第一層網格的y+值小于1。經過網格無關性驗證,布置球窩、球凸的U型通道網格數分別確定為360萬和270萬。

2 數值處理方法

2.1 數值方法

首先,將直通道中布置單個球窩時流動傳熱的實驗結果[12]與湍流模型計算結果進行了比較。計算時邊界條件與實驗相同[12]:進口雷諾數為21 000,總溫為293 K,湍流強度為5%,出口靜壓為0.1 MPa,壁面應用無滑移條件,給定恒定壁溫323 K。然后,利用CFX軟件分別采用4種湍流模型求解了三維RANS方程。圖2給出了數值模擬結果與實驗數據[12]比較。可以看到,標準k-ω湍流模型預測的努塞爾數Nu分布與實驗值最為吻合。另外,本文模型中也包括了U型通道180°彎頭對流動換熱的影響。文獻[11]研究指出,采用k-ω湍流模型能夠準確預測U型通道內部的流動換熱過程,因此本文后續計算均采用標準k-ω湍流模型。

X:流動方向圖2 球窩尾跡區努塞爾數的分布

文中以蒸汽作為冷氣工質,蒸汽的熱物理性質由IAPWS-IF97真實蒸汽性質模型計算,流動與傳熱邊界條件:進口湍流度為5%,上、下壁面的溫度保持在573.15 K,其余壁面均為絕熱無滑移邊界。進口速度Ui由當量直徑Dh和雷諾數Re決定,即

(1)

計算收斂標準為各控制方程的殘差小于10-5。

2.2 數據處理

在等壁溫條件下,流動截面上的局部溫差沿流動方向是不斷變化的,平均的對流傳熱溫差

(2)

式中:Twi、Two、Tmi和Tmo分別為U型通道沿流向微元段的進、出口壁面溫度和進、出口截面流體的平均溫度。定義來流溫比φ=Tin/Tw,其中Tin為冷氣進口總溫,Tw為壁面溫度。努塞爾數

(3)

式中:qw為壁面熱流密度;λ為流體的導熱系數。為了直觀起見,以Nu/Nu0表示球窩/球凸結構引起的強化傳熱效果,其中Nu0為相應流動參數下光滑通道中充分發展流動的努塞爾數[3]

(4)

式中:f0為光滑通道中充分發展流動的阻力系數[3]

f0=(0.79lnRe-1.64)-2

(5)

實際通道中阻力系數

(6)

式中:u為各截面上的平均流速。

3 計算結果與分析

3.1U型通道中球窩/球凸結構的流動特性

首先分析了球窩/球凸結構附近的流動特征,由于各球窩/球凸附近的流動結構類似,所以此處僅給出進口Re為30 000、出口壓力為1 MPa時流入段第6排球窩/球凸結構附近的流動細節。

如圖3a所示,在球窩前緣,由幾何結構形成的非光滑過渡產生較強的逆壓梯度,引起球窩前緣發生流動分離,并在球窩的尾緣再附。如圖3b所示,流動分離使得球窩內形成了2個關于中心流面對稱的反向旋轉渦對,這一低壓渦將球窩兩側的主流卷吸進入球窩,并與邊界層流體混合后在球窩后緣處加速射入主流。在主流的作用下,射流渦量在球窩尾跡區轉向且引起二次流。

對于球凸情況,在球凸前緣流體出現滯止點,而邊界層內的剪切流動以滯止點為界向兩側發生偏轉,從而形成了馬蹄渦結構(見圖3c),該渦在主流攜帶下向下游遷移,并在球凸尾跡區與越過球凸表面且發生分離的流體摻混后再附(見圖3d)。

(a)球窩底部再附流線 (b)球窩兩側卷吸流線

(c)球凸前緣分離流線 (d)球凸尾緣分離流線圖3 球窩/球凸結構表面的三維流線

3.2 雷諾數對蒸汽流動傳熱的影響

為研究Re變化對球窩/球凸通道傳熱[13]與阻力特性的影響,在來流溫比保持0.79不變、出口壓力為1 MPa的情況下分析了Re從10 000變化至50 000時的流動傳熱情況。

(a)流入段球窩 (b)流出段球窩

(c)流入段球凸 (d)流出段球凸圖4 流入段和流出段第6個球窩/球凸表面Nu及速度矢量分布(Re=30 000)

圖4給出了Re=30 000時流入段和流出段第6個球窩/球凸結構表面Nu及速度矢量分布。可以看到,流出段Nu明顯大于流入段。

在球窩前緣附近的分離區Nu較小,尾緣再附區Nu較大,尾緣駐點處Nu達到峰值;在球凸前緣區Nu較大,尾緣分離區Nu較小,前緣駐點處Nu達到峰值。

圖5給出了不同Re下球窩/球凸通道流出段的Nu分布。同一Re下,通道中各個球窩/球凸附近的Nu分布類似,整體趨勢為沿流動方向降低。隨著Re增大,球窩和球凸通道的換熱效果均增強。這是由于Re增大時主流加速,壁面動量邊界層厚度在主流的拖曳作用下減薄,溫度梯度增大,使得傳熱得到了強化。

從圖5a可以看到,每個球窩結構分別在尾緣駐點和下游尾跡區出現了2個Nu的峰值。該結果可以從流動與傳熱相互影響的微觀機理來分析。球窩通道流速增大使得凹面附近曲率對湍流邊界層的放大作用得以增強[14],邊界層與主流的摻混更加強烈,更多的主流低溫流體被帶入邊界層內,并在球窩尾緣形成如圖4a、圖4b所示的高速射流,同時導致溫度邊界層的最薄點出現在球窩尾緣,由此形成了Nu的第1個峰值;射流內不同流層之間具有速度梯度,射流在球窩尾跡內的流動轉折造成了尾跡區的分離渦,而分離渦的存在強化了尾跡區邊界層與主流的熱質交換,同時Re增加又使得分離渦強度增大,由此在球窩尾跡區內形成了Nu的第2個峰值。

對于球凸通道,Re增加使得前緣滯止區增大,邊界層流體擾動加強,流速提高又使得溫度邊界層變薄、傳熱增強。由于凸面曲率對湍流具有壓制作用[14],所以尾緣流動的分離提前、再附位置推后,尾緣分離渦增強,這使得有更多的主流低溫流體被帶入邊界層內,尾緣回流渦低速區空間被壓縮,回流渦強度降低,傳熱效果增強,從而造成圖5b的Nu分布。

另外,從圖5還看到,Re增大時U型通道180°轉折處的Nu逐漸增大。這表明轉折處二次流對傳熱的強化作用不容忽視,其機理如圖6所示。經180°轉折后流體湍動能增加,這也是流出段Nu高于流入段的原因。

(a)球窩U型通道

(b)球凸U型通道圖5 不同Re下球窩/球凸流出段的Nu分布

圖6 180°轉折處的三維流線及湍動能變化

圖7 Re變化時球窩/球凸通道的傳熱與阻力特性

3.3 來流溫比對蒸汽流動傳熱的影響

對進口Re=30 000、出口壓力為1 MPa、φ從0.79到0.93變化時的流動傳熱情況進行了分析,結果如圖8所示。球窩/球凸通道的強化傳熱效果隨來流溫比的增加變化不明顯,而阻力系數逐步降低。

圖8 不同來流溫比下球窩/球凸通道的傳熱與阻力特性

圖9給出了流入段第6、7排球窩/球凸結構表面Nu和極限流線隨來流溫比的變化。在球窩尾緣和球凸前緣分別形成了半月牙形和W形的高Nu區,來流溫比增大時兩種結構表面的Nu都略有增加,這是來流溫比增大導致流速增大引起的,其微觀機理與上節分析類似。隨著來流溫比的增大,水蒸氣的普朗特數Pr由0.991降至0.969,傳熱能力減弱,但是溫比上升又使得水蒸氣的密度減小、流速增加,從而導致壁面的邊界層減薄、傳熱增強,所以綜合作用的結果是傳熱系數略有增加。在流動阻力方面,來流溫比上升使得水蒸氣的黏性系數增大、摩擦阻力增大,速度增大使得湍流的脈動增強、分離推遲、壓差阻力減小,所以綜合作用結果是阻力系數降低。

3.4 出口壓力對蒸汽流動傳熱的影響

(a)球窩結構

圖11給出流入段第6排球窩/球凸結構D/2處與主流方向相垂直的橫截面內的二次流線分布。由于流動結構對稱,圖中只給出1/2圓弧附近的流線分布。圖11a中,球窩內形成對稱的旋轉渦對,結合圖3知,這是球窩前緣發展出的馬蹄渦結構,主要作用是將主流流體卷吸進入球窩內部,但在強化傳熱的同時并不會帶來較大的阻力損失。當出口壓力上升到2 MPa時,原來的渦對結構中又嵌套出現了一對渦,而馬蹄渦結構也被下層流體抬升。渦系結構的改變導致在垂直于流動方向的截面內流體微團之間的相互作用增強,動能耗散為熱能的比例增加,流動阻力顯著增強。另外,出口壓力上升使得流速急劇下降,黏性邊界層增厚,這也是流動阻力上升的一個原因。如圖11c和圖11d所示,球凸結構在不同出口壓力下的流動結構變化與球窩類似,此處不再分析。

(a)球窩,0.5 MPa (b)球窩,2 MPa

(c)球凸,0.5 MPa (d)球凸,2 MPa圖11 不同出口壓力下球窩/球凸附近流線

3.5 綜合傳熱性能評估

為了全面評估球窩/球凸結構的傳熱與阻力特性,定義綜合傳熱增強因子[3]

(7)

它是在推動流體流動的機械功一定的條件下球窩/球凸表面與光滑表面的傳熱性能之比。

圖12為綜合傳熱性能對比。可以看到,球凸通道的綜合傳熱性能要優于球窩通道,這與文獻[9]的結論一致。隨著Re和φ的增加,球窩和球凸通道的綜合傳熱性能提升,但φ的增加使得有效熱通量降低;隨著出口壓力Po的上升,阻力系數對綜合換熱效能的影響逐漸增大,球窩和球凸通道的綜合傳熱性能逐漸降低。綜合考慮傳熱效能和熱通量的影響,選取中等壓力(1.5~2 MPa)和中等來流溫比(0.84~0.88),相應的綜合傳熱性能能夠獲得最為理想的效果。

圖12 綜合傳熱性能對比

(8)

采用回歸算法分別對球窩/球凸通道的16個數據點進行公式擬合,得到球窩傳熱關聯式Nudim和球凸的傳熱關聯式Nupro分別為

(9)

(10)

其中1×104≤Re≤5×104,0.79≤φ≤0.93,0.5≤ζ≤3。

圖13給出了球窩和球凸Nu關聯式解y與數值解x的相關性分布。結果表明,大多數數據點集圍繞在斜率為1的直線附近,并保持在斜率為0.95和1.05的夾角范圍內。這表明關聯式在±5%范圍內保持與數值解的相關性,即傳熱關聯式曲線與數值計算點吻合良好。利用該傳熱關聯式能夠可靠地預測相應的傳熱系數。

圖13 球窩和球凸關聯式與數值解的相關性分布

3.6 蒸汽與空氣冷卻效果對比

為對比相同幾何及邊界參數下蒸汽與空氣冷卻的傳熱特性,定義傳熱增強系數[10]

(11)

圖14 蒸汽與空氣冷卻性能的對比

圖14為蒸汽與空氣冷卻性能的對比。球凸通道和球窩通道的η分布基本相同。Re增加時,η保持為21.2%左右;φ從0.79增至0.93,η由21.2%下降至19.7%;出口壓力Po從0.5 MPa增至3 MPa,η由21%增大至24.5%。上述變化趨勢主要是由兩種工質Pr的不同變化特性[15]所導致。空氣的Pr對邊界參數的變化不太敏感。對于蒸汽,隨著進口雷諾數變化Pr基本不變,進口溫比增加時Pr減小,出口壓力增加時Pr增大。值得注意的是,出口壓力對η的影響比φ更為強烈。

4 結 論

對球窩/球凸結構下U型通道內的蒸汽傳熱及阻力性能進行了數值研究,主要結論如下。

(1)進口雷諾數、來流溫比和出口壓力對通道內蒸汽冷卻有重要影響。來流溫比不變,進口雷諾數增加,傳熱強化,流動阻力降低。相同雷諾數下,來流溫比增大對傳熱的強化作用不明顯,但流動阻力降低。提高出口壓力可以有效強化傳熱,但同時也使流動阻力大大增加。所以,選取中等壓力(1.5~2 MPa)和中等來流溫比(0.84~0.88)時蒸汽冷卻能夠獲得最為理想的綜合換熱效果。

(2)對球窩/球凸通道內蒸汽傳熱關聯式進行了擬合并與數值計算結果進行了對比,結果表明,擬合的傳熱關聯式與數值計算結果吻合良好,可預測相應通道內的換熱系數。

(3)在本文計算參數下,與空氣相比,蒸汽的傳熱增強系數為19.7%~24.5%,蒸汽的普朗特數受溫度和壓力的影響更明顯,其傳熱增強系數隨來流溫比和出口壓力的變化更為顯著。Re變化對傳熱增強系數基本無影響。

[1] AFANASYEV V N, CHUDNOVSKY Y P, LEONTIEV A I, et al. Turbulent flow friction and heat transfer characteristics for spherical cavities on a flat plate [J]. Experimental Thermal and Fluid Science, 1993, 7(1): 1-8.

[2] MOON H K, O’CONNELL T, GLEZER B. Channel height effect on heat transfer and friction in a dimpled passage [J]. ASME Journal of Engineering for Gas Turbines and Power, 2000, 122(2): 307-313.

[3] MOON S W, LAU S C. Turbulent heat transfer measurements on a wall with concave and cylindrical dimples in a square channel [C]∥ASME Turbo Expo 2002: Power for Land, Sea, and Air. New York, USA: ASME, 2002: 459-467.

[4] BURGESS N K, LIGRANI P M. Effects of dimple depth on channel Nusselt numbers and friction factors [J]. Journal of Heat Transfer, 2005, 127(8): 839-847.

[5] BURGESS N K, OLIVEIRA M M, LIGRANI P M. Nusselt number behavior on deep dimpled surfaces within a channel [J]. Journal of Heat Transfer, 2003, 125(1): 11-18.

[6] LIGRANI P M, HARRISON J L, MAHMMOD G I, et al. Flow structure due to dimple depressions on a channel surface [J]. Physics of Fluids, 2001, 13(11): 3442-3451.

[7] MAHMOOD G I, SABBAGH M Z, LIGRANI P M. Heat transfer in a channel with dimples and protrusions on opposite walls [J]. Journal of Thermophysics and Heat Transfer, 2001, 15(3): 275-283.

[8] RAO Y, WAN C, ZANG S. Comparisons of flow friction and heat transfer performance in rectangular channels with pin fin-dimple, pin fin and dimple arrays [C]∥ASME Turbo Expo 2010: Power for Land, Sea, and Air. New York, USA: ASME, 2010: 185-195.

[9] SHEN Z, QU H, ZHANG D, et al. Effect of bleed hole on flow and heat transfer performance of U-shaped channel with dimple structure [J]. International Journal of Heat and Mass Transfer, 2013, 66: 10-22.

[10]SHUI L, GAO J, SHI X, et al. The effect of cooling conditions on convective heat transfer and flow in a steam-cooled ribbed duct [J]. Journal of Mechanical Science and Technology, 2014, 28(1): 331-341.

[11]STEPHENS M A, SHIH T I P, CIVINSKAS K C. Computations of flow and heat transfer in a rotating U-shaped square duct with smooth walls, AIAA 1996-3161 [R]. Reston, VA, USA: AIAA, 1996.

[12]PATRICK W V, TAFTI D K. Computations of flow structure and heat transfer in a dimpled channel at low to moderate Reynolds number [C]∥ASME 2004 Heat Transfer/Fluids Engineering Summer Conference. New York, USA: ASME, 2004: 401-412.

[13]武俊梅, 陶文銓. 縱向渦強化換熱的數值研究及場協同原理分析 [J]. 西安交通大學學報, 2006, 40(7): 757-761. WU Junmei, TAO Wenquan. Numerical analysis to vortex heat transfer enhancement based on field synergy principle [J]. Journal of Xi’an Jiaotong University, 2006, 40(7): 757-761.

[14]DIXON S L. Fluid mechanics and thermodynamics of turbomachinery [M]. Oxford, UK: Butterworth-Heinemann, 2005.

[15]王文, 劉志剛. 氣體普朗特數變化規律的初步研究 [J]. 西安交通大學學報, 1999, 33(1): 77-80. WANG Wen, LIU Zhigang. Variation of Parandtl number of gas in dilute and dense state [J]. Journal of Xi’an Jiaotong University, 1999, 33(1): 77-80.

(編輯 苗凌)

(西安交通大學能源與動力工程學院,710049,西安)

采用標準k-ω湍流模型研究了球窩/球凸結構下U型通道內蒸汽的傳熱與阻力特性,計算、分析了雷諾數為1×104~5×104、來流溫比為0.79~0.93、出口壓力為0.5~3 MPa時通道內的流動與傳熱性能,獲得了擬合的傳熱關聯式,對比了蒸汽與空氣冷卻特性的差異。研究結果表明:雷諾數增加在提升總體傳熱效果的同時能夠有效降低流動阻力系數;來流溫比增大對強化傳熱的作用并不明顯,但可以有效降低流動阻力損失;提高出口壓力可以有效增強傳熱效果,但大大增加了流動阻力;出口壓力為1.5~2 MPa、來流溫比為0.84~0.88時蒸汽冷卻參數能夠獲得較為理想的綜合換熱效果;與空氣相比,蒸汽的熱物性更容易受到溫度和壓力的影響,導致傳熱增強系數的變化較大,此時雷諾數的變化對傳熱增強系數的貢獻并不顯著;擬合的傳熱關聯式可以準確預測通道內蒸汽冷卻換熱系數。

球窩/球凸;U型通道;蒸汽冷卻;換熱系數;數值模擬

Heat Transfer Performance in Steam-Cooled U-Shaped Channel with Concave/Convex Structure

LI Sen,LI Liang,DU Changhe

(School of Energy and Power Engineering, Xi’an Jiaotong University, Xi’an 710049, China)

Heat transfer performance and flow resistance characteristics of a steam-cooled U-shaped channel with concave/convex structure were numerically investigated by standardk-ωturbulent model. The heat transfer correlations were evaluated and analyzed for Reynolds number from 10 000 to 50 000, inlet temperature ratio from 0.79 to 0.93, and outlet pressure from 0.5 to 3 MPa. The different cooling performances of steam and air were discussed. It reveals that the overall heat transfer performance rises while the friction factor falls with the increasing Reynolds number; the increasing inlet temperature ratio effectively enables to reduce the flow resistance and enhance heat transfer performance obviously; the heightening outlet pressure is beneficial to heat transfer enhancement, but it also leads to greatly increased flow resistance. Thus, a combination of moderate steam pressure from 1.5 MPa to 2 MPa, and medium inlet temperature ratio from 0.84 to 0.88 is suggested to achieve an ideal heat transfer performance. The thermophysical properties of steam are more prone to be affected by temperature and pressure compared with those of air. As a result, the heat transfer performance evaluation criterion for steam varies obviously while the Reynolds number exerts slight influence.

concave/convex structure; U-shaped channel; steam cooling; heat transfer coefficient; numerical investigation

2014-12-02. 作者簡介:李森(1991—),男,碩士生;李亮(通信作者),男,博士,副教授。

時間:2015-06-17

http:∥www.cnki.net/kcms/detail/61.1069.T.20150617.0902.002.html

10.7652/xjtuxb201509012

TK474.7

A

0253-987X(2015)09-0063-07

猜你喜歡
結構
DNA結構的發現
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
新型平衡塊結構的應用
模具制造(2019年3期)2019-06-06 02:10:54
循環結構謹防“死循環”
論《日出》的結構
縱向結構
縱向結構
我國社會結構的重建
人間(2015年21期)2015-03-11 15:23:21
創新治理結構促進中小企業持續成長
現代企業(2015年9期)2015-02-28 18:56:50
主站蜘蛛池模板: 一本大道在线一本久道| 色婷婷成人| 国产精品夜夜嗨视频免费视频 | 国产欧美精品一区aⅴ影院| 91在线一9|永久视频在线| 特级欧美视频aaaaaa| 日本五区在线不卡精品| 国产自在自线午夜精品视频| 伦精品一区二区三区视频| 国产不卡在线看| 亚洲精品自拍区在线观看| 国产哺乳奶水91在线播放| 亚洲精品自拍区在线观看| 亚洲人成日本在线观看| 丁香六月激情婷婷| 久久永久视频| a级毛片免费网站| 激情無極限的亚洲一区免费| 久久9966精品国产免费| 夜精品a一区二区三区| 精品国产毛片| 久久人人97超碰人人澡爱香蕉| 日韩毛片免费观看| 偷拍久久网| 亚洲人在线| 三上悠亚精品二区在线观看| a级毛片在线免费| 亚洲欧州色色免费AV| 欧美人人干| 国产AV毛片| 欧美视频在线播放观看免费福利资源| 国产在线自乱拍播放| 国产欧美日韩18| 毛片视频网址| 免费国产无遮挡又黄又爽| 日韩a在线观看免费观看| 2020亚洲精品无码| 欧美天天干| 毛片基地视频| 亚洲乱码精品久久久久..| 精品久久高清| 波多野结衣视频一区二区 | 欧美精品影院| www.精品视频| 成人国产精品一级毛片天堂| 久久久久国产一区二区| 成人午夜网址| 国产精品欧美激情| 亚洲h视频在线| 亚洲最新在线| 亚洲AⅤ永久无码精品毛片| 福利在线不卡一区| 亚洲欧洲自拍拍偷午夜色无码| 婷婷亚洲最大| 日韩在线欧美在线| 国产中文在线亚洲精品官网| 免费av一区二区三区在线| 成人日韩视频| 亚洲第一香蕉视频| 国产在线自乱拍播放| 无码免费视频| 日韩国产综合精选| 黄色网页在线播放| 久久亚洲日本不卡一区二区| 毛片视频网址| 亚洲制服丝袜第一页| 国产尤物jk自慰制服喷水| 欧美在线视频不卡第一页| 理论片一区| 国产精品蜜芽在线观看| 国产丝袜丝视频在线观看| 9966国产精品视频| 亚州AV秘 一区二区三区| 久久亚洲AⅤ无码精品午夜麻豆| 国产日韩欧美视频| 日本国产在线| 成人免费网站久久久| 日本亚洲国产一区二区三区| 伊人五月丁香综合AⅤ| 好久久免费视频高清| 中文字幕首页系列人妻| 亚洲Va中文字幕久久一区 |