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

雷諾數對圓柱氣動力和流場影響的試驗研究

2016-07-05 12:53:35劉慶寬鄭云飛李聰輝馬文勇劉小兵
實驗流體力學 2016年4期
關鍵詞:模型研究

劉慶寬,邵 奇,鄭云飛,李聰輝,馬文勇,劉小兵

雷諾數對圓柱氣動力和流場影響的試驗研究

劉慶寬1,2,*,邵 奇3,鄭云飛3,李聰輝3,馬文勇1,2,劉小兵1,2

(1.石家莊鐵道大學大型結構健康診斷與控制研究所,石家莊 050043;2.河北省大型結構健康診斷與控制重點實驗室,石家莊 050043;3.石家莊鐵道大學土木工程學院,石家莊 050043)

通過剛性模型測壓風洞試驗,研究了圓柱的氣動阻力、氣動升力系數和風壓系數隨雷諾數的變化規律,從流場分布的角度分析了氣動力變化的原因,并研究了雷諾數影響下的流場在圓柱軸向的相關性。結果表明:在亞臨界雷諾數區域,在時間平均上流場沿模型兩側呈對稱分布,雷諾數對平均阻力系數和流場影響較小,平均升力系數基本為零。在臨界雷諾數區域,隨著特定區域大負壓區的出現,流場不再對稱,出現不容忽視的平均升力和脈動升力。在超臨界雷諾數區域,隨著對稱側大負壓區的出現,流場恢復對稱狀態,平均升力基本消失。雷諾數對流場的軸向相關性有顯著的影響。在雷諾數較低時(亞臨界區域),卡門渦在軸向上的尺度相對較大,而隨著雷諾數的提高,該尺度逐漸減小,各斷面流場的相關性降低。

圓柱體;氣動力;流場分布;雷諾數;軸向相關性

0 引 言

隨著斜拉橋設計和施工工藝的發展,其跨徑逐漸增大,斜拉索上的風荷載占全橋風荷載的比例也隨之增大。以蘇通長江大橋為例,在橫向風的作用下,斜拉索上的風荷載對主梁位移和內力的貢獻占整個風荷載的60%~70%[1]。因此,準確掌握斜拉索上的風荷載,對于斜拉橋的設計具有重要意義。同時,由于斜拉索長細比大,柔度大,阻尼小,容易發生風致振動,尤其是風雨振因其振幅大、危害嚴重等,廣受設計和研究人員的重視。

作為斜拉橋主要受力構件的斜拉索為典型的圓形斷面,作為基本的結構斷面形狀,針對其氣動特性研究者進行了較為廣泛研究。V·Karmen、Taylor、Wieselsbeger、Fage、Warsap等人先后進行了試驗研究,Roshko發現在超臨界區邊界層內出現分離氣泡并伴隨湍流再附[2];Bearman發現在亞臨界與超臨界區之間存在單側分離氣泡,結構出現升力[3]。

上述研究表明,圓形斷面的氣動力與雷諾數、來流湍流度及結構表面的粗糙度密切相關。因為大跨徑斜拉橋多建在近海地區,經常受到臺風的侵擾,基本風速較大,同時由于通航的要求,其主梁和斜拉索等構件離水面較高,加上江面開闊(地表粗糙度低),這些因素使得其設計基準風速較高,斜拉索的雷諾數數值也較大,從低風速到設計基準風速,對應的雷諾數經常跨越亞臨界、臨界和超臨界雷諾數。如蘇通長江大橋A1索的索梁端和索塔端靜陣風風速對應的雷諾數分別為4.56×105和5.20×105[4-5]。處在不同雷諾數區的氣動力和流場都有明顯的不同。

因此無論是斜拉索氣動力計算,還是風致振動的機理研究,雷諾數效應引起的特殊的氣動力及流場的變化、流場沿結構軸向相關性的變化等問題是需要考慮的重要因素[6-10],有必要進行專門的研究。

本研究將斜拉索簡化為二維圓柱結構,通過風洞試驗研究了氣動力和流場隨雷諾數的變化,并研究了不同雷諾數下流場沿模型軸向分布相關性的問題,為氣動力計算和機理分析提供依據和參考。

1 風洞試驗介紹

試驗在石家莊鐵道大學風工程研究中心的STU-1風洞內進行。該風洞為串聯雙試驗段回/直流邊界層風洞,其低速試驗段轉盤中心處寬4.4m,高3.0m,長24.0m,最大風速大于30.0m/s;高速試驗段寬2.2m,高2.0m,長5.0m,最大風速大于80.0m/s[11]。風洞結構如圖1所示。

圖1 風洞結構示意圖Fig.1 Sketch of wind tunnel

本研究共進行了2種試驗,第1種試驗在低速試驗段中進行,通過表面測壓研究雷諾數對斜拉索氣動力和流場的影響,風場為均勻低湍流度風場(湍流度I≤0.4%),模型安裝及測點布置如圖2所示,模型兩端與風洞頂部和底部固定連接,為了保證二維流動,在模型兩端各安裝Φ900mm的端板。在模型圓周方向布置了一圈180個測壓孔,其中編號A001的測壓孔正對來流方向,試驗工況如表1所示。

第2種試驗在高速試驗段中進行,研究軸向流場的相關性問題,風場為均勻低湍流度風場(湍流度I≤0.2%),模型安裝及測點布置示意圖如圖3所示。

圖2 低速試驗段模型安裝及測點布置圖Fig.2 Sketch of installation and pressure tap arrangement of test model in low speed test section

表1 低速試驗段模型試驗工況Table 1 Test cases in low speed test section

圖3 高速試驗段模型安裝及測點布置圖Fig.3 Sketch of installation and pressure tap arrangement of test model in high speed test section

該工況共布置了A、B、C、D、E共5圈測點,共有8種軸向間距,將測點軸向間距與模型直徑d(d=120mm)的比值定義為無量綱軸向間距s,選取7個無量綱軸向間距s(s=1、2、3、4、5、6、7)進行相關性研究。同樣每圈共布置60個測壓孔,其中編號A01的測壓孔正對來流方向,試驗工況如表2所示。

表2 高速試驗段模型試驗工況Table 2 Test cases in high speed test section

壓力測試采用PSI電子掃描閥進行,壓力信號采用文獻[12]的方法進行了修正。以減小模型的阻塞度和獲得較大雷諾數為原則,選取了低速試驗段和高速試驗段的模型直徑,兩者的阻塞度分別為6.8%和6.0%,通過Two-step Maskell的方法[13-14]對氣動力進行了相應的阻塞度修正。

氣動力的計算采用風軸坐標系,計算簡圖如圖4所示,積分方法如公式(1)和(2)所示。

圖4 模型氣動力計算圖Fig.4 Aerodynamic force calculation diagram of test model

式中:CD為平均阻力系數;CL為平均升力系數;pi為第i個測壓點的動壓;Li為第i個測壓點的弧長;θi為第i個測壓點與束流方向的夾角;Cpi為第i個測壓點的平均風壓系數。

值得說明的是,因為是通過測壓孔測試模型表面法向壓力,并通過積分得到的阻力和升力,空氣與模型表面的摩擦力無法計入其內,所以理論上測試結果應比實際值略小。

2 雷諾數對氣動力的影響

模型的氣動阻力系數和氣動升力系數隨雷諾數的變化情況如圖5所示。

圖5 氣動力系數隨雷諾數的變化曲線Fig.5 Curve of aerodynamic force coefficient with Re

從圖5中可以看出,在Re<3.42×105的范圍內,平均阻力系數緩慢減小,平均升力系數基本為0,是亞臨界雷諾數區域,根據前人研究結果,此區域結構表面的邊界層在分離之前仍處于層流狀態。在3.42×105≤Re≤3.73×105的范圍內,平均阻力系數快速下降,同時出現較大的平均升力,并隨雷諾數的增大,阻力系數的變化趨于平緩,平均升力基本消失,是臨界雷諾數區域;在進行風壓分析中發現,結構兩側的風壓不對稱,說明出現了單側分離氣泡;在Re>3.72×105的范圍內,阻力系數和升力系數基本不再隨雷諾數的增大而變化,是超臨界雷諾數區域。

目前對于圓柱形斜拉索的氣動力計算,相關規范[15]只考慮了阻力,沒有考慮升力。該氣動升力對整體風荷載的貢獻,是值得考慮的問題。而目前斜拉索氣動力的研究,針對阻力和脈動升力的較多,對臨界雷諾數區域的平均升力沒有引起足夠的重視[16-18]。

3 雷諾數對風壓分布的影響

由上述內容可知氣動力系數隨著雷諾數的改變而發生較大變化。這是因為雷諾數的改變導致圓柱周圍的流場發生了變化,為了分析模型周圍的流場,分別取亞臨界、臨界和超臨界雷諾數范圍內的3個雷諾數數值,對其表面風壓進行了分析。

圖6為Re=2.05×105(亞臨界)時模型表面平均風壓和脈動風壓的分布情況,本文結果與文獻[19]的試驗值總體趨勢一致,在60°~300°范圍內大于文獻[19]的試驗值,本研究中模型的雷諾數、模型表面的光滑度等因素與文獻[19]之間存在的差異可能是導致差別的原因。由圖可知,0°~180°的分布同180°~360°的分布基本對稱,說明模型兩側的流場在時間平均上呈對稱狀態,平均升力為0。

圖6 亞臨界雷諾數區域模型風壓系數周向分布(Re=2.05×105)Fig.6 Wind pressure coefficient distribution around test model in subcritical regime(Re=2.05×105)

圖7 和8分別是瞬時表面壓力系數和駐點在一個卡門渦脫落周期內的變化情況。隨著卡門渦的交替脫落,駐點、分離點等在模型表面交替擺動。駐點擺動的周期對應卡門渦的脫落周期、脈動升力的周期和2倍脈動阻力的周期。其中圖8描述的是在亞臨界雷諾數區域一個周期T內,模型駐點位置隨時間的擺動。擺動的角度是指駐點附近最大風壓系數值對應位置與擺動范圍的中心點所夾的角度。從圖中可以看出擺動的角度隨時間大致呈正弦曲線有規律地變化。

圖7 亞臨界雷諾數區域模型瞬時風壓周向分布向量圖(Re=2.05×105)Fig.7 Instantaneous wind pressure distribution vector diagram around test model in subcritical regime(Re=2.05×105)

圖8 亞臨界雷諾數區域模型駐點位置擺動圖(Re=2.05×105)Fig.8 Stagnation point position of test model in subcritical regime(Re=2.05×105)

圖9 為Re=3.55×105(臨界)時模型表面平均風壓和脈動風壓的分布情況。由圖9可知,此時模型兩側的流場已經不再對稱,一側的負壓絕對值明顯大于另一側,同時該側的脈動風壓也明顯變大。

圖10是瞬時表面壓力系數的變化情況。此時的流場隨時間基本不再變化,在特定區域有不隨時間變化的較大負壓出現。該不對稱的負壓,是平均升力產生的原因。

圖11和12為Re=3.79×105(超臨界)時模型表面平均風壓和脈動風壓的分布、瞬時表面風壓系數的變化情況。由圖可知,此時模型另一側也出現大絕對值的負壓,同原來出現的負壓接近;同時2個大負壓區的脈動風壓也接近一致。隨著流場恢復到對稱狀態,平均升力基本消失。

需要說明的是,本研究只針對光滑表面模型進行了研究。針對斜拉索,為了防止風雨振采用氣動措施(表面纏繞螺旋線、設置凹坑等)后的情況有待進一步研究。

圖9 臨界雷諾數區域模型風壓系數周向分布(Re=3.55×105)Fig.9 Wind pressure coefficient distribution around test model in critical regime(Re=3.55×105)

圖10 臨界雷諾數區域模型瞬時風壓周向分布向量圖(Re=3.55×105)Fig.10 Instantaneous wind pressure distribution vector diagram around test model in critical regime(Re=3.55×105)

圖11 超臨界區域模型風壓系數周向分布(Re=3.79×105)Fig.11 Wind pressure coefficient distribution around test model in supercritical regime(Re=3.79×105)

圖12 超臨界雷諾數區域模型瞬時風壓周向分布向量圖(Re=3.79×105)Fig.12 Instantaneous wind pressure distribution vector diagram around test model in supercritical regime(Re=3.79×105)

4 雷諾數影響下流場軸向的相關性

為了研究雷諾數影響下的流場沿模型軸向的相關性,定義了各圈測點之間的無量綱軸向間距s(1≤s≤7,s=L/d,其中L指2排測點間的距離,d指模型的直徑),以脈動升力系數為指標,來考察流場沿軸向的相關性。評價相關性的程度有2個無量綱參數,相關系數和相干函數,一個是時域上的分析,一個是頻域上的分析,對于本研究的相關性分析,采用的是相關系數,如式(3)所示:

式中:rxy為2個各態歷經平穩過程中的脈動信號x和y的互相關系數;N為采樣的樣本數;n為樣本中第n個樣本;Rxx(0)為脈動信號x初始時的自相關系數;Ryy(0)為脈動信號y初始時的自相關系數。

為了定性和定量地說明相關性問題,通常按表3給出的相關系數與相關程度對照表來進行說明。

表3 相關系數與相關程度對照表Table 3 Table of correlation coefficient and degree

通過改變雷諾數(1.0×105≤Re≤3.2×105)和無量綱軸向間距s(1≤s≤7),得到脈動升力系數的相關系數,如圖13所示。

圖13 相關系數同雷諾數和無量綱軸向間距的關系Fig.13 Relation between correlation coefficient,Re,and dimensionless axial spaceing

從圖中可以看出,隨著雷諾數的增大、無量綱間距的增大,相關性相應降低。對于區域1≤s<3和1.0×105≤Re<1.6×105,模型升力系數相關性較強,為顯著相關或高度相關,而對于其它區域相關系數很小,基本為微相關。

由此可推知,在雷諾數較低時(亞臨界區域),結構表面為層流分離,卡門渦在軸向上的尺度相對較大;隨著雷諾數的提高(臨界區域),結構表面的邊界層出現層流和湍流之間相互轉換。這個過程對結構表面粗糙度較為敏感,導致不同位置處的分離出現差異,使其在軸向的分布不均勻,從而導致軸向相關系數逐漸減小。

由于測試設備的量程所限,該相關性研究的測試只做到了雷諾數為3.2×105,高雷諾數下的相關性問題有待進一步研究。

5 結 論

本文通過風洞試驗,研究了雷諾數對圓柱氣動力和流場特性的影響,并考慮了流場在軸向上的相關性問題,得到如下結論:

(1)在亞臨界雷諾數區域,在時間平均上流場沿模型兩側呈對稱分布,雷諾數對平均阻力系數和流場影響不大,平均升力系數基本為0。在臨界雷諾數區域,隨著特定區域大負壓區的出現,流場不再對稱,平均升力出現,并出現較大的脈動升力。在超臨界雷諾數區域,隨著對稱側大負壓區的出現,流場恢復對稱狀態,平均升力消失。

(2)雷諾數對流場的軸向相關性有顯著的影響。在雷諾數較低時(亞臨界區域),卡門渦在軸向上的尺度相對較大,而隨著雷諾數的提高,該尺度逐漸減小,各斷面流場的相關性降低。

(3)對于具有圓形結構的斜拉索,由于其氣動力與風速的平方和氣動力系數成正比,在臨界區,隨著風速的增大阻力系數減小,同時出現較大的升力,因此整個風速范圍內的最大氣動力不一定在最大風速處取得,準確的氣動力需要根據研究結果具體計算。

(4)本文只是針對光滑表面的圓柱進行了分析,對于采用氣動減振措施(如纏繞螺旋線、設置凹坑等)的斜拉索的氣動力和流場特性隨雷諾數的變化有必要進一步研究。同時流場的軸向相關性在高雷諾數下的狀態也有必要進行深入研究。

[1]裴岷山,張喜剛,朱斌,等.斜拉橋的拉索縱橋向風荷載計算方法研究[J].中國工程科學,2009,11(3):26-30.Pei M S,Zhang X G,Zhu B,et al.Study on longitudinal wind load calculation method of cables for cable-stayed bridge[J].Engineering Sciences,2009,11(3):26-30.

[2]Roshko A.Experiments on the flow past a circle cylinder at very high Reynolds numbers[J].Journal of Fluids Mechanics,1961,10:345-356.

[3]Bearman P W.On vortex shedding from a circle cylinder in the critical Reynolds number regime[J].Journal of Fluids Mechanics,1968,37:577-585.

[4]崔冰,曾憲武.南京二橋南汊大橋主橋結構設計[C]//中國土木工程學會橋梁及結構工程學會第十三屆年會論文集,北京:人民交通出版社,1998:264-272.Cui B,Zeng X W.Structure design of the nancha bridge of nanjing yangtse river bridge[C]//Proceedings of the 13th National Conference on Bridge and Structure Engineering,Beijing:China Communications Press,1998:264-272.

[5]張忠義,劉聰.南京長江第二大橋橋位風速觀測及設計風速的計算[J].氣象科學,2000,20(2):200-205.Zhang Z Y,Liu C.Calculation of desighing wind velocity and wind observation at the Nanjing 2nd Yangtse River Bridge[J].Scientia Meteorologica Sinica,2000,20(2):200-205.

[6]劉慶寬,鄭云飛,白雨潤,等.斜拉索風雨振氣動抑振措施的參數優化[J].振動與沖擊,2015,34(8):31-35.Liu Q K,Zheng Y F,Bai Y R,et al.Parametric optimization of aerodynamic and vibration measure for rain-wind induced vibration of cables[J].Journal of Vibration and Shock,2015,34(8):31-35.

[7]劉慶寬,王毅,鄭云飛,等.水線-雷諾數效應-斜拉索振動關系的試驗研究[J].工程力學,2012,29(11):257-265.Liu Q K,Wang Y,Zheng Y F,et al.Experimental study on the relation of water rivulet-Reynolds number effect-cable vibration[J].Engineering Mechanics,2012,29(11):257-265.

[8]劉慶寬,張峰,馬文勇,等.斜拉索雷諾數效應與風致振動的試驗研究[J].振動與沖擊,2011,30(12):114-119.Liu Q K,Zhang F,Ma W Y,et al.Tests for Reynolds number effect and wind-induced vibration of stay cables[J].Journal of Vibration and Shock,2011,30(12):114-119.

[9]Chen S,Irwin P A,Jakobsen J B,et al.Divergent motion of cables exposed to skewed wind[C]//Proceedings of the 5th International Symposium on Cable Dynamics.Santa Margherita Ligure,Italy,2003:271-278.

[10]Matsumoto M,Yagi T,Shima T,et al.The effect of the Reynolds number in aerodynamic cable vibration of cable-stayed bridge[C]//Proceedings of the 19th National Symposium on Wind Engineering.Tokyo,Japan,2006:507-512.

[11]劉慶寬.多功能大氣邊界層風洞的設計與建設[J].實驗流體力學,2011,25(3):66-70.Liu Q K.Aerodynamic and structure design of multifunction boundary-layer wind tunnel[J].Journal of Experiments in Fluid Mechanics,2011,25(3):66-70.

[12]馬文勇,劉慶寬,劉小兵,等.風洞試驗中測壓管路信號畸變及修正研究[J].實驗流體力學,2013,27(4):71-77.Ma W Y,Liu Q K,Liu X B,et al.Study on correction and distortion effects caused by tubing systems of pressure measurements in wind tunnel[J].Journal of Experiments in Fluid Mechanics,2013,27(4):71-77.

[13]Hackett J E.Tunnel-induced gradients and their effect on drag[J].American Institute of Aeronautics and Astronautics,1996,34(12):2575-2581.

[14]Cooper K R,Mercker E,Wiedemann J.Improved blockage corrections for bluff bodies in closed and open wind tunnels[C]//Wind Engineering into 21stCentury:Proceedings of the 10thInternational Conference on Wind Engineering.Copenhagen,Denmark,1999:1627-1634.

[15]JTG/T D60-01—2004.公路橋梁抗風設計規范[S].北京:人民交通出版社,2004.JTG/T D60-01—2004.Wind-resistant design specification for highway bridges[S].Beijing:China Communications Press,2004.

[16]Poulin S,Larsen A.Drag loading of circular cylinders inclined in the along-wind direction[J].Journal of Wind Eingineering and Industrial Aerodynamics,2007,95(9/10/11):1350-1363.

[17]王衛華,李明水,陳忻.斜拉索的阻力系數研究[J].空氣動力學學報,2005,23(3):389-393.Wang W H,Li M S,Chen X.Investigation on drag coefficients of stay-cables[J].Acta Aerodynamica Sinica,2005,23(3):389-393.

[18]林志興,楊立波,李文勃.斜拉索順橋向風阻系數的試驗研究[J].鄭州大學學報:工學版,2005,26(1):38-41.Lin Z X,Yang L B,Li W B.Experimental study on drag coefficients of stay-cables corresponding to wind direction along the bridge central line[J].Journal of Zhengzhou University:Engineering Science Edition,2005,26(1):38-41.

[19]Farell C,Blessm A J.On critical flow around smooth circular cylinders[J].Journal of Fluid Mechanics,1983,136:375-391.

Experimental study on Reynolds number effect on aerodynamic pressure and forces of cylinder

Liu Qingkuan1,2,*,Shao Qi3,Zheng Yunfei3,Li Conghui3,Ma Wenyong1,2,Liu Xiaobing1,2
(1.Structural Health Monitoring and Control Institute,Shijiazhuang Tiedao University,Shijiazhuang 050043,China;2.Hebei Province Key Lab of Structural Health Monitoring and Control,Shijiazhuang 050043,China;3.School of Civil Engineering,Shijiazhuang Tiedao University,Shijiazhuang 050043,China)

By wind tunnel tests,Reynolds number effect on drag force coefficient,lift force coefficient and wind pressure coefficient of cylinder were measured,the mechanism of the aerodynamic force variation from the point of view of the flow field distribution was analyzed,and the correlation of the flow field in the cylinder axis direction under the Reynolds number effect was studied.Results show that in the subcritical Reynolds number regime,the time averaged flow field around the cylinder model is symmetric,the Reynolds number has little influence on the averaged drag force coefficient and flow field,and the averaged lift force is around zero.In the critical Reynolds number regime,with the appearance of large amplitude negative pressure in certain areas,flow field becomes asymmetric,and the averaged lift force as well as fluctuation lift force appears.In the supercritical Reynolds number regime,with the appearance of large amplitude negative pressure on the two sides,flow field becomes to symmetric again,and the averaged lift force disappears.Reynolds number has obvious effect on the flow field correlation along the cylinder axis.In the subcritical Reynolds number regime,the scale of the Karman vortex in the cylinder axis direction is relatively large,while the scale becomes small with the increase of Reynolds number.

cylinder;aerodynamic force;flow field distribution;Reynolds number;correlation in axis direction

U441+.3

:A

Liu Q K,Shao Q,Zheng Y F,et al.Experimental study on Reynolds number effect on aerodynamic pressure and forces of cylinder.Journal of Experiments in Fluid Mechanics,2016,30(4):7-13.劉慶寬,邵奇,鄭云飛,等.雷諾數對圓柱氣動力和流場影響的試驗研究.實驗流體力學,2016,30(4):7-13.

(編輯:楊 娟)

1672-9897(2016)04-0007-07

10.11729/syltlx20150112

2015-08-28;

2016-04-28

國家自然科學基金項目(51378323,51108280,51308359);河北省杰出青年科學基金(E2014210138)

*通信作者E-mail:lqk@stdu.edu.cn

劉慶寬(1971-),男,河北保定人,博士,教授,博士生導師。研究方向:橋梁與結構的風荷載、風致振動與控制。通信地址:石家莊鐵道大學風工程研究中心(050043)。E-mail:lqk@stdu.edu.cn

猜你喜歡
模型研究
一半模型
FMS與YBT相關性的實證研究
2020年國內翻譯研究述評
遼代千人邑研究述論
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
視錯覺在平面設計中的應用與研究
科技傳播(2019年22期)2020-01-14 03:06:54
EMA伺服控制系統研究
新版C-NCAP側面碰撞假人損傷研究
3D打印中的模型分割與打包
主站蜘蛛池模板: 亚洲天堂网在线视频| 亚洲男人的天堂网| 在线国产91| 国产精品免费久久久久影院无码| 国内视频精品| 久久一级电影| 国产微拍精品| 99精品久久精品| 麻豆精品在线播放| 日韩高清中文字幕| 欧美午夜网| 欧美日韩综合网| 在线观看国产精品一区| 福利片91| 免费Aⅴ片在线观看蜜芽Tⅴ | 伊人久久久久久久| 日韩美毛片| 国产69精品久久久久妇女| 欧美日韩导航| 成人在线天堂| 色综合中文| 国产91全国探花系列在线播放| 国产一级妓女av网站| 国产在线高清一级毛片| 精品久久蜜桃| 99ri精品视频在线观看播放| 久久综合伊人77777| 原味小视频在线www国产| 国产永久无码观看在线| 精品视频一区在线观看| 国产1区2区在线观看| 国产午夜看片| 国产91九色在线播放| 亚洲精品午夜无码电影网| 亚洲午夜天堂| 露脸国产精品自产在线播| 午夜高清国产拍精品| 三级国产在线观看| 日韩精品久久久久久久电影蜜臀| 热99re99首页精品亚洲五月天| 免费一极毛片| 精品国产香蕉在线播出| 91在线免费公开视频| 久久国产精品国产自线拍| 91人人妻人人做人人爽男同| 国产精品视频免费网站| 91毛片网| 国产又粗又猛又爽视频| 成人免费视频一区二区三区| 欧美va亚洲va香蕉在线| 99久久精品国产自免费| 国产拍揄自揄精品视频网站| 亚洲一区二区三区国产精品| 日本草草视频在线观看| 久久精品欧美一区二区| 婷婷六月在线| 国产福利一区视频| 色吊丝av中文字幕| 99视频在线免费| 国产精品分类视频分类一区| 国产尤物jk自慰制服喷水| a毛片免费观看| 精品国产中文一级毛片在线看 | 精品国产aⅴ一区二区三区| 狂欢视频在线观看不卡| 欧美在线综合视频| 中文字幕在线日韩91| 999在线免费视频| 国产在线视频导航| 午夜精品区| 91午夜福利在线观看精品| 中文字幕人妻av一区二区| 久久精品国产电影| 色亚洲成人| 一级毛片免费的| 最新无码专区超级碰碰碰| 国产精品色婷婷在线观看| 欧洲日本亚洲中文字幕| 中文字幕波多野不卡一区| 亚洲成A人V欧美综合| 伊人成人在线| 国产制服丝袜无码视频|