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

湍流積分尺度對高層建筑風荷載影響的大渦模擬

2018-06-01 02:59:41祝志文鄧燕華
西南交通大學學報 2018年3期
關鍵詞:風速

祝志文, 鄧燕華

(湖南大學土木工程學院, 湖南 長沙 410082)

結構抗風研究中,無論是風洞試驗還是數值模擬,準確模擬大氣邊界層的湍流特性是保證研究結果可信的必要條件.早期的研究主要關心風速剖面和湍流強度剖面,以期得到與大氣邊界層相似的風速和湍流強度.盧占斌等[1]的風洞試驗研究表明,湍流強度主要改變風壓脈動值的大小和能量在各頻譜范圍的分布,湍流積分尺度的增大將使風壓系數平均值的絕對值和均方根值均降低.JR等[2]研究矩形橋梁斷面在風速和湍流強度不變的前提下,湍流積分尺度對風壓系數平均值和均方根值、彎矩譜峰值的影響.Nakamura和Ozono[3]指出,當湍流積分尺度小于2D(D為結構特征長度)時,可以忽略其對湍流流態的影響;當大于2D時,湍流積分尺度增大,將對風壓系數平均值產生顯著影響.Li和Melbourne[4]研究湍流積分尺度在2.8D~7.8D范圍內對結構風效應的影響,隨湍流積分尺度的增加,小尺度湍流脈動運動量降低,壓力系數峰值呈現增大趨勢,當湍流積分尺度長度達到5D時出現降低現象.

風洞試驗中要準確測定湍流積分尺度和保持湍流強度不變的前提下調整其值絕非易事.本文利用大渦模擬研究湍流積分尺度對高層建筑風荷載的影響.大渦模擬入口脈動風場生成的方法主要包括隨機流生成技術(random fow generation, RFG)、離散和合成隨機流生成法(discretizing and synthesizing random flow generation, DSRFG).文獻[5]研究指出,RFG方法通過湍動能和耗散率求湍流積分尺度,DSRFG方法通過調整湍流積分尺度來調整流場特性,具有接近真實流的優點,但計算費用是RFG的fmax倍(fmax為采樣頻率的最大值),且普通計算機難于實現,只能在超算中心或超大集成機群上完成.本文通過調整湍流積分尺度,提高RFG計算精度,從而規避DSRFG方法中的不足.

1 數值模型

1.1 幾何模型和測點布置

高層建筑模型采用CAARC (commonwealth advisory aeronautical research council)建筑,縮尺比為 1∶300,縮尺后的幾何尺寸:長(Dy)152 mm×寬(Dx)102 mm×高(H)610 mm.豎向設8個測點層,各層高度見圖1.

每個測點層設20個測點,測點總數160個.定義來流風速垂直于截面長邊,測點1~5、6~10、11~15和16~20所在的面分別為迎風面、上側面、背風面和下側面,見圖2.

1.2 計算域及網格劃分

順流向(x軸向)入口至模型中心的距離為7.5H,模型中心至出口的距離為15H,上下側面距模型中心均為3.5H,模型頂部到上邊界的距離為6H,計算得到的堵塞率為0.6%,顯著小于風洞試驗堵塞度上限值3%.

1.3 邊界條件設定

邊界條件設置為:入口設為速度,出口為自由流出口,地面和模型表面均為無滑移壁面,其他均為對稱邊界,模型壁面和地面的首層網格高度均為2×10-4Dx.

圖1 測點層分布Fig.1 Arrangement of pressure layers

圖2 測點平面布置Fig.2 Cross-section arrangement of pressure points

1.4 湍流模型

本文基于LES (large eddy simulation)湍流模型開展CFD (computational fluid dynamics)模擬,其控制方程為

(1)

(2)

Cij=pij+Dij+Φij-Eij,

(3)

(4)

pij=-

(5)

(6)

(7)

(8)

式中:Cij為雷諾應力在平均運動軌跡上的增長率;為湍動能的生成項;pij為雷諾應力通過平均運動的變形率向湍流脈動輸入平均能量,當其大于0時,平均運動向脈動運動輸入能量,湍動能增大,反之,湍動能減?。籇ij為雷諾應力擴散項;Φij為再分配項,在湍流脈動各分量間起調節作用,實現能量的轉移;Eij為湍動能耗散項;表示相應變量的系綜平均值(即在給定的邊界條件下一切可能運動狀態的算術平均);u下標的i、j和k為速度各方向的分量;u′下標的i、j和k為脈動速度各方向的分量;x下標的i、j和k為坐標軸x的i、j和k方向;p′為壓強脈動值;ν是流體運動黏度;

采用Simple算法進行壓力速度耦合,動量方程采用中心差分離散.每個時間步最大迭代20次,經時間無關性檢查,確定時間步長b=0.000 5 s.

1.5 計算參數設置

為后續更好地與文獻試驗結果進行比較,入口風速和湍流強度均選用文獻[7]的B類地貌風洞試驗結果,風速為

U=12.98(z/1.167)0.16.

(9)

湍流強度為

I=0.085(z/1.167)-0.21.

(10)

湍流積分尺度Lu參照文獻[8]并計入縮尺比后得到,即

Lu=1.054z0.5.

(11)

湍動能K和湍動能耗散率ε均參照文獻[19],即

K=1.2(UI)2,

(12)

ε=0.090.75K1.5/(φLu),

(13)

式(9)~(13)中:

z為距地面高度;

φ為湍流積分尺度調整系數.

本文保持風速和湍流強度不變,將φ分別設置為0.4和1.0,得到相應的湍流積分尺度長度與鈍體特征尺度長度之比k1=2.17和k2=5.42.將上述參數通過Fluent用戶自定義函數(user defined function, UDF)導入Fluent軟件.

2 數據處理

平均風壓系數和脈動風壓系數按文獻[7]定義,其它氣動參數按式(14)~(15)定義.

(14)

(15)

式(14)~(15)中:

Cd和Cl分別為平均阻力系數和升力系數;

CMx和CMy分別為順風向和橫風向基底彎矩系數平均值;

UH和Dy分別為模型頂部平均風速和模型迎風面寬度;

FD、FL分別為基底阻力和升力;

3 計算結果比較與分析

3.1 基底氣動力的對比

從表1可以看出,本文LES模擬結果與文獻結果總體上趨于一致,模擬獲得的氣動力均方根值比風洞試驗值小,但平均值與試驗值吻合較好,模擬結果與文獻[12]的數值模擬結果較接近.

表1 模擬結果與文獻結果的比較Tab.1 Comparison of simulation results with other available

3.2 表面風壓系數

圖3給出了2/3H高處的測點風壓系數平均值Cp和均方根Cprms,圖中,橫坐標N為監測點編號,具體布置見圖2.用機構名稱中的關鍵詞代表其機構名稱,分別為:City代表City University, Bristol University of Bristol; Monash代表Monash University; NAE代表National Aeronautical Establishment(a、b分別表示該機構采用的兩種大氣邊界層的模擬方法);NPL代表Nmional Physical Laboratory; TJ-2B代表同濟大學TJ-2建筑風洞B類地貌風場.LES為本文的大渦模擬結果,具體數據及試驗方法參見文獻[7].

從圖3可知,不同文獻給出的風壓系數平均值和均方根均存在一定差別,特別是風壓系數均方根.本文計算的風壓系數平均值和均方根值無論在正壓區還負壓區,均與風洞試驗值存在一致性,表明湍流模型的選擇和計算參數的設置均是可靠的.湍流積分尺度對風壓系數有不同程度的影響(圖3(a)).正壓區風壓系數平均值隨湍流積分尺度增大而增大2%~5%,即湍流積分尺度對正壓區風壓系數平均值的影響很小.在負壓區,隨湍流積分尺度的增大,平均風壓系數絕對值減小,側面平均風壓系數絕對值降低12%~17%,背面降低13%,即湍流積分尺度主要對側面和背面風壓系數的影響顯著,這與文獻[3]的結論一致.與文獻[10]中City、Bristol、Monash、NAE和文獻[7]的風洞試驗結果相比,湍流積分尺度調整系數φ=0.4 時,平均風壓系數與其結果較一致,尤其是側面和背面.

(a) 2/3H高度處平均風壓系數

(b) 2/3H高度處風壓系數均方根圖3 湍流積分尺度對風壓系數的影響Fig.3 Influence of turbulence integral lengths on wind pressure coefficients

圖3(b)可知,無論是正壓區還是負壓區,風壓系數均方根值均隨湍流積分尺度的增大而降低,且角部降幅最大,這可能是受流體的分離與附著的影響.迎風面除角部顯著降低外,其他測點風壓系數均方根相差很小,側面和背面最大降幅分別為15%和10%.湍流積分尺度調整系數φ=0.4時,模擬值與試驗結果更一致,即湍流積分尺度的增大將降低風壓脈動性.

綜上分析可知,隨湍流積分尺度的增大,除迎風面風壓系數平均值較風洞試驗結果略大外,其他部位均偏小,且風壓系數均方根值降低,這與文獻[1]的試驗結果一致.分析其原因,從式(8)可知,增大湍流積分尺度,若保持湍流強度不變,湍動能的耗散率將降低.文獻[4]指出,增大湍流積分尺度將減少小尺度的脈動運動量.Richardson[4]認為小尺度的脈動是將湍動能向熱能的轉變,所以湍動能耗散率的降低和小尺度脈動運動減少勢必造成系統能量不平衡,要維持平衡將迫使湍動能增大,即提高耗散率.由式(7)可知,湍動能增大的物理表現是湍流強度的提高,由式(3)中的pij項可知,湍流強度的提高必須增大平均運動的變形率,平均運動的變形率向湍流脈動輸入平均能量,這將導致平均速度減小,以致出現迎風面平均風壓系數隨湍流積分尺度的增大而略降現象(圖3(a)).因湍流積分尺度增大導致湍流強度提高,文獻[13]的研究指出,湍流強度的提高分流剪切流附著性提前出現,渦脫強度降低,空氣夾帶能力減弱,這將導致負壓區的平均風壓系數絕對值和風壓系數均方根均減小(圖3).側面風壓與分離流的附著、來流風速和湍流強度有關,背面風壓受旋渦結構、旋渦分布等有關,受來流風速和湍流強度的影響很小.

3.3 氣動力沿高度的變化

為研究湍流積分尺度對阻力的影響,計算了不同湍流積分尺度的測點層高度z處的平均阻力系數,見圖4.

圖4 層平均阻力系數沿高度變化Fig.4 Mean drag coefficient of the layer at different height

由圖4可知:平均阻力系數沿高度的變化與文獻[14]和文獻[15]風洞試驗結果基本一致,且模擬結果與文獻[14]更接近,沿高度方向風速逐漸增大,平均阻力系數也逐漸增大,至0.85H處,頂部受三維流態的影響,其值出現降低現象[14],相比風洞試驗,數值模擬受三維流態的影響偏小,這可能是兩者的邊界條件不同所致;湍流積分尺度越大,平均阻力系數越小,且沿高度方向變化,降幅5%~10%,其原因是阻力包括迎風面來流風與壁面的碰撞和背面旋渦脫落產生的負壓,湍流積分尺度的增大雖然使得迎風面風速略增,但背面分離流附著性增強,負壓絕對值減小,以致阻力系數總體上表現為降低.

圖5 斯托羅哈數(Sr)沿高度變化Fig.5 Strouhal number at different height

由圖5可知,湍流積分尺度并沒有改變Sr沿高度變化趨勢.湍流積分尺度越大,Sr沿高度增大相對越小,0.4H高度以下兩者相差很小,0.4H至0.9H其值降低了約30%,0.9H以上其值降低了約20%,可見湍流積分尺度越大,渦脫落頻率越低.

結合圖4和圖5可知:在1/3H以下,湍流積分尺度對測點層平均阻力系數和斯托羅哈數的影響不明顯;當高度大于1/3H時,其影響越來越明顯,這是因為風速沿高度逐漸增大,底部雖然湍流強度大,但風速偏小.可見湍流積分尺度對高層建筑風荷載有明顯的影響,且風速越大,影響越顯著.

3.4 流線分析

為獲得不同湍流積分尺度下的繞流特征,利用CFD的優勢,獲得2/3H高處x-y平面的流線圖;為重點研究細部繞流特性,選取了關鍵部位流線圖進行分析,見圖6.

由圖6可知:

與φ=1相比,當φ=0.4時,迎風面剪切分離流距側面越遠,旋渦的尺度更大,側面和背面均存在大尺度旋渦,表現為較強的渦脫性,尾流寬度大;當φ=1時,雖然旋渦的尺度小,但尾流寬度小,剪切分離流在側面和背面均表現較好的附著性.

(a) 湍流積分尺度調整系數為0.4時(b) 湍流積分尺度調整系數為1時圖6 z=2/3H高處x-y平面流線Fig.6 x-y plan streamlined diagram at z=2/3H height

3.5 頻譜比較與分析

以下進一步研究湍流積分尺度對基底力矩譜的影響,見圖7(圖中:橫坐標為斯托羅哈數;縱坐標中σMz、σMx和σMy分別為基底扭矩均方根、基底順風向彎矩均方根和基底橫風向彎矩均方根,SMz(f)、SMx(f)和SMy(f)分別為基底扭矩功率譜密度、順風向基底彎矩功率譜密度和橫風向基底彎矩功率譜密度.

由圖7(a)可知,隨湍流積分尺度增大,基底扭矩譜的峰值減小,但降幅很小,頻帶明顯變窄,高頻段幅值降低,原因是增大湍流積分尺度,伴隨小尺度脈動運動量減少[4],附著提前出現(見圖6),這不但導致紊流對扭矩的貢獻降低,而且改變了能量在各頻域帶的分布,譜頻帶變窄.

由圖7(b)可知,順風向基底彎矩譜峰值不明顯,譜頻帶較寬,不同湍流積分尺度的順風向彎矩譜發展趨勢基本一致,湍流積分尺度越大,高頻段幅值越低,這是由于湍流積分尺度的增大導致背面旋渦脫落強度降低所致.

由圖7(c)可知,隨湍流積分尺度增大,橫風向基底彎矩譜峰值基本不變,但渦脫頻率降低、高頻段的幅值減小,低頻段幅值略顯升高趨勢,原因是隨湍流積分尺度的增大,大尺度的脈動運動相對增多和小尺度脈動減少[4].

Richardson[4]認為:

大尺度湍流脈動是蓄能池,小尺度湍流脈動將湍動能轉變成熱能進行耗散,即能量主要集中在大尺度的湍流脈動上,所以低頻段出現能譜升高、高頻段幅值降低的現象;其次,橫風向基底彎矩主要來源于漩渦脫落誘發的氣動力,紊流對其也有一定貢獻.隨湍流積分尺度的增大,湍流強度增大,側面風壓系數絕對值減小也有關.

(a) 基底扭矩譜(b) 順風向基底彎矩譜(c) 橫風向基底彎矩譜圖7 不同湍流積分尺度的氣動力譜Fig.7 Aerodynamic forces spectra with different turbulence integral lengths

3.6 相關性

相關系數是兩個隨機變量存在線性依賴的一種測度,湍流積分尺度對氣動力相關性的影響是一個值得重視的問題,本文定義模型高度方向上任意兩層氣動力時程或同一水平面上模型表面任意兩點壓力時程r、s的相關系數為

(16)

式中:r=1,2,…,R;s=1,2,…,n;CoVF(zr,zs)為高度zr和zs風力協方差;σF(zr)和σF(zs)和分別為高度zr和zs的風力根方差;n(或R)為模型測點層數.

3.6.1水平相關性

以第五層測點(即H=406 mm)為例進行風壓水平相關性分析,按照式(11)計算迎風面、上側面和背面各自面內中間測點與其它測點的風壓相關系數,見圖8.圖中:x為距中點的距離,迎風面下部為負,上部為正;側面左邊為負,右邊為正;背面上部為負,下部為正;d為測點所在面的邊長.由圖8可知:兩測點的距離越大,風壓相關系數越小;相同距離測點相關系數值表現為迎風面的最大,背風面最??;相關系數隨湍流積分尺度增大而減小,且測點間距越大,差異越明顯;隨測點距離的減小,不同湍流積分尺度計算得到的相關系數的差值呈現降低趨勢,且不同的面存在較明顯差異.

迎風面湍流積分尺度小時,風壓水平相關曲線表現良好的對稱性(圖8(a)),這表明旋渦分布更均勻.隨湍流積分尺度的增大,風壓水平相關系數最大降幅為4.5%, 可見湍流積分尺度對迎風面風壓相關性的影響很小.側面風壓水平相關曲線表現為湍流積分尺度越小,風壓水平相關曲線對稱性越差,尤其在下風向更突出,甚至出現驟降現象.隨湍流積分尺度的增大,側面上風向風壓水平相關系數的增幅很小,但下風向增幅達到15%~25%(圖8(b)).這是由于湍流積分尺度小時湍流強度也小,側面旋渦脫落性強,風壓脈動性強(圖3),以致風壓相關性差,且由于側面上下風向旋渦尺度和旋渦強度不同(圖6),所以脈動風壓不同步,以致相關性較差,甚至出現驟降.

由圖8(c)可知:不同湍流積分尺度背面風壓相關系數相差5%~10%,這與背面風壓與背面旋渦的尺度、結構及分布的影響有關,受來流風速和湍流強度的影響相對較小.雖然增大湍流積分尺度會導致湍流強度增大,但湍流強度的變化對背面風壓的影響很小.

可見風壓水平相關性受湍流積分尺度的影響較大,由于湍流積分尺度的增大,湍流強度增大,側面分離剪切流在側面提前附著,導致側面和背面的平均風壓系數絕對值和脈動值均降低,風壓的相關系數提高.

3.6.2豎向相關性

以第4層測點為基準,計算不同湍流積分尺度時各測點層升力相關系數(RL)和阻力相關系數(RD),見圖9.

由圖9可知:在0.8H以下,不同湍流積分尺度的測點層升力和阻力相關系數的差值很小,可不考慮其對相關性的影響;隨高度的增加,湍流積分尺度的影響逐漸增大,值得注意的是層升力相關系數的變化很小,層阻力相關系數的增幅頗為明顯,尤其是靠近建筑物頂部,最大增幅約50%,這可能受頂部三維流態的影響所致[14].雖然湍流積分尺度對0.8H高度以上的層阻力相關系數影響明顯,但影響范圍有限,可以認為湍流積分尺度對氣動力的豎向相關性的影響很小,建議不予考慮.

(a) 迎風面測點(b) 側面測點(c) 背面測點圖8 z=2/3H高處各面測點風壓水平相關系數Fig.8 Correlation coefficients of wind pressure in horizontal on respective plane at z=2/3H height

(a) RL(b) RD圖9 不同測點層升力相關系數和阻力相關系數Fig.9 Correlation coefficients of drags and lifts among different layers

4 結 論

本文采用大渦模擬的方法,在入口處保持速度、湍流強度、湍動能不變的前提下,研究了湍流積分尺度對高層建筑風荷載的影響,得到以下結論:

(1) 增大湍流積分尺度會使得各面的風壓系數均方根均降低,迎風面平均風壓系數略增,側面和背風面平均風壓系數絕對值均減小.湍流積分尺度為k1時的模擬結果與試驗值更吻合.

(2) 隨湍流積分尺度的增大,平均運動的變形率向湍流脈動輸入平均能量,導致平均風速略降,湍流強度略增;橫風向漩渦脫落頻率減小;基底扭矩譜和橫風向基底彎矩譜的峰值及高頻段譜幅值均呈降低趨勢.B類地貌湍流積分尺度調整系數為0.4時,計算精度更高,計算結果與文獻實驗值更吻合.

(3) 隨湍流積分尺度的增大,測點層阻力系數平均值和斯托羅哈數均沿高度方向呈降低趨勢,高度1/3H以下時表現不明顯,但1/3H以上時影響逐漸增強.

(4) 湍流積分尺度的增大,會導致湍流強度增大,側面和背面剪切分離流附著性提前出現,旋渦脫落能力減弱,表現為斯托羅哈數值降低,負壓區風壓得以恢復,風壓水平相關性增強,但對豎向風壓相關性的影響很小.

[1] 盧占斌,魏慶鼎. 網格湍流CAARC模型風洞實驗[J]. 空氣動力學學報,2001,19(1): 16-23.

LU Zhanbin, WEI Qingding. An experiment on a CAARC model in grid turbulent flow[J]. Acta Aerodynamica Sinica, 2001,19(1): 16-23.

[2] JR F L H, KAREEM A, SZEWCZYK A A. The effects of turbulence on the pressure distribution around a rectangular prism[J]. Journal of Wind Engineering and Industrial Aerodynamics, 1998, 77-78(98): 381-392.

[3] NAKAMURA Y, OZONO S. The effects of turbulence on a separated and reattaching flow[J]. Journal of Fluid Mechanics, 1987, 178: 477-490.

[4] LI Q S, MELBOURNE W H. The effects of large scale turbulence on pressure fluctuations in separated and reattaching flows[J].Journal of Wind Engineering and Industrial Aerodynamics.1999, 83(1/2/3): 159-169.

[5] Huang S H, Li Q S, Wu J R. A general inflow turbulence generator for large eddy simulation[J]. Journal of Wind Engineering and Industrial Aerodynamics, 2010, 98(10): 600-617.

[6] 張兆順,崔桂香,許曉春. 湍流大渦模擬的理論與應用[M]. 北京:清華大學出版社,2005: 4-42.

[7] 羅盤. 基于標準模型的風洞試驗研究[D]. 上海:同濟大學土木工程學院,2004.

[8] TAMURA Y, OHKUMA T, KAWAI H, et al. Revision of AIJ recommendations for wind loads on buildings[C]∥Structures Congress. Nashville: [s.n.], 2004: 1-10.

[9] LAKEHAL D. Application of thek-εmodel to flow over a building placed in different roughness sublayers[J]. Journal of Wind Engineering and Industrial Aerodynamics, 1998, 73(1): 59-77.

[10] MELBOURNE W H. Comparison of measurements of the CAARC standard tall building model in simulated model wind flows[J]. Journal of Wind Engineering & Industrial Aerodynamics, 1980, 6(1): 73-88.

[11] OBASAJU E D. Measurement of forces and base overturning moments on the CAARC tall building model in a simulated atmospheric boundary layer[J]. Journal of Wind Engineering and Industrial Aerodynamics,1992, 40(2): 103-126.

[12] 楊有根. 高層建筑的抗風實測分析與大跨度屋蓋的風致響應研究[D]. 長沙:湖南大學,2007.

[13] 蘇萬林,李正農. 湍流對超高層建筑風壓幅值特性影響的研究[J]. 地震工程與工程振動,2016,36(3): 118-126.

SU Wanlin, LI Zhengnong. Research of turbulent effects on characteristics of wind pressure amplitude on tall buildings[J]. Earthquake Engineering and Engineering Dynamics, 2016, 36(3): 118-126.

[14] 唐意. 高層建筑彎扭耦合風致振動及靜力等效風荷載研究[D]. 上海:同濟大學土木工程學院,2006.

[15] 陳建蘭. 矩形截面高層建筑風荷載試驗研究[D]. 武漢:華中科技大學,2005.

猜你喜歡
風速
邯鄲市近46年風向風速特征分析
基于Kmeans-VMD-LSTM的短期風速預測
基于最優TS評分和頻率匹配的江蘇近海風速訂正
海洋通報(2020年5期)2021-01-14 09:26:54
基于時間相關性的風速威布爾分布優化方法
陜西黃土高原地區日極大風速的統計推算方法
陜西氣象(2020年2期)2020-06-08 00:54:38
基于GARCH的短時風速預測方法
快速評估風電場50年一遇最大風速的算法
風能(2016年11期)2016-03-04 05:24:00
考慮風切和塔影效應的風力機風速模型
電測與儀表(2015年8期)2015-04-09 11:50:06
GE在中國發布2.3-116低風速智能風機
考慮風速分布與日非平穩性的風速數據預處理方法研究
主站蜘蛛池模板: 麻豆精品在线| jizz在线观看| 在线观看无码av五月花| 国内自拍久第一页| 精品乱码久久久久久久| 国产全黄a一级毛片| 亚洲 欧美 日韩综合一区| 日本欧美精品| 久久精品免费国产大片| 久久人人爽人人爽人人片aV东京热 | 国产麻豆永久视频| 自拍偷拍欧美日韩| 亚洲成人动漫在线观看| 国内视频精品| 成人国产精品视频频| 无码国产伊人| 欧美三级自拍| 色综合国产| 亚洲欧洲日韩久久狠狠爱| 免费欧美一级| 午夜福利免费视频| 国产精品入口麻豆| 中文字幕不卡免费高清视频| 欧美成人精品在线| 91精品aⅴ无码中文字字幕蜜桃 | 亚洲天堂高清| 欧美a在线看| 人妻少妇乱子伦精品无码专区毛片| 亚洲色图欧美激情| 久久国产亚洲欧美日韩精品| 亚洲第一天堂无码专区| 国产精品第一区| 亚洲成人黄色在线| 国产精品亚洲五月天高清| 国产女人在线观看| 少妇精品久久久一区二区三区| 91娇喘视频| 91成人免费观看| 国产成人艳妇AA视频在线| 日韩精品无码一级毛片免费| 亚洲人视频在线观看| 日韩A级毛片一区二区三区| 亚洲AⅤ无码国产精品| 欧美一区日韩一区中文字幕页| 麻豆国产原创视频在线播放 | 国产欧美精品一区二区| 无码中文字幕乱码免费2| 91视频99| 中文字幕欧美日韩高清| 呦女精品网站| 国产成人免费手机在线观看视频 | 欧美亚洲一区二区三区在线| 亚洲最黄视频| 久久99热66这里只有精品一| 九九热精品视频在线| 米奇精品一区二区三区| 亚洲欧美不卡中文字幕| 国产成人精品优优av| 综合色亚洲| 美女毛片在线| 亚洲国产精品无码久久一线| 日本www在线视频| 国产精品爽爽va在线无码观看| 国产屁屁影院| 精品91自产拍在线| 2022精品国偷自产免费观看| 免费国产高清视频| 露脸一二三区国语对白| 亚洲精品福利视频| 婷婷色一二三区波多野衣| 制服丝袜一区| 亚洲人成网站色7777| 91精品专区| 欧美激情二区三区| 2022国产91精品久久久久久| 99一级毛片| 亚洲永久精品ww47国产| 全免费a级毛片免费看不卡| 午夜精品影院| 亚洲 日韩 激情 无码 中出| 99资源在线| 国产18页|