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

湍流特性對矩形高層建筑風荷載影響的研究

2021-05-04 03:27:14陳水福
空氣動力學學報 2021年2期

劉 奕,陳水福

(浙江大學 建筑工程學院,杭州 31005)

0 引 言

風荷載是高層建筑的主要控制荷載,而風的湍流特性是影響風荷載分布特性尤其是脈動與極值特性的重要因素。探討來流湍流特性對建筑風荷載的影響一直是學術界和工程界十分關注的課題。

湍流度和湍流積分尺度是表征風的湍流特性的兩個主要參數,現有研究也主要從這兩個方面出發研究其對建筑風荷載及周圍繞流特性的影響。文獻[1-3]研究了湍流特性對二維矩形棱柱繞流的影響,其中Lee[1]發現湍流度增大使側風面風壓恢復提前、背風面負壓(均指絕對值)減小;Laneville[2]發現湍流度增大使平均剪切層曲率半徑減小,從而使阻力系數最大值對應的深寬比減小;Nakamura和Ohya[3]發現湍流強度和湍流積分尺度均對背風面負壓有較大影響。文獻[4-9]研究了湍流特性對二維平板鈍體周圍分離再附流的影響,其中Hillier和Cherry[4]發現湍流度對平均風壓和脈動風壓均有很大影響,而湍流積分尺度對平均風壓影響不大,對脈動風壓影響很大;Nakamura和Ozono[5]發現分離泡內平均風壓分布在湍流積分尺度比(即湍流積分尺度與模型特征尺寸之比)小于2時變化不大,而超過2時逐漸趨近均勻流下的分布;Li和Melbourne[8-9]發現平均風壓在尺度比大于3.6時受湍流積分尺度影響明顯,脈動風壓和極值負壓隨湍流度和湍流積分尺度增大而增大,湍流積分尺度對極值負壓的影響程度隨湍流度增大而提高。顧明和葉豐[10-11]研究了不同風場下不同體型高層建筑的氣動力幅值和頻域特性,發現隨地面粗糙度增加,脈動升力、阻力系數增大,升力系數譜頻帶變寬、能量減小;蘇萬林[12]研究了湍流對高層建筑表面風壓的影響,發現在湍流強度很大時,背風面和側風面可能出現小正壓;Akon和Kopp[13]通過粒子圖像測速試驗,研究了湍流度和湍流積分尺度對低矮房屋屋頂分離再附流的影響,發現湍流強度增大會使流動更早再附,且影響風壓恢復速率,而湍流積分尺度對再附長度影響不大。總體而言,目前國內外湍流特性研究主要針對深寬比較小的矩形截面鈍體或建筑,或是深寬比無限大的二維平板,但對于介于二者之間的有限深寬比的鈍體研究較少。但是,這類鈍體隨著深寬比的變化,其表面風壓是否會因流動分離再附的改變而表現得顯著不同,其中來流湍流特性的影響程度如何,目前并不明確。

鑒于此,本文以國內廣泛存在的具有較大深寬比范圍的板式住宅高層建筑為背景,對4種不同地貌條件下9種深寬比的矩形截面高層建筑進行了風洞測壓試驗,分析了湍流度與湍流積分尺度對不同深寬比建筑平均、脈動與極值風壓,以及橫風向氣動力譜和周圍繞流特性的影響規律。

1 風洞試驗

1.1 風場模擬

本次試驗在加拿大西安大略大學(University of Western Ontario,UWO)邊界層風洞Ⅱ的高速試驗段中進行。該試驗段長30.0 m,寬3.4 m,高2.1 m,最大風速可達40 m/s。試驗段前端擁有全自動的粗糙元控制系統,配合入口處可選放置的格柵、尖劈和擋板等湍流產生裝置,可以精確模擬各種地貌下的大氣邊界層。

本次試驗共模擬了4種地貌,其中O1、O2地貌與S1、S2地貌使用了兩種不同的粗糙元配置,分別模擬開闊地貌(Open terrain)與郊區地貌(Suburban terrain)風場;O1、S1地貌較O2、S2地貌在風洞入口處多使用了一塊0.38 m高的擋板。根據Holmes和Osonphasop[14],擋板的存在可以產生更大尺度的湍流,實現對湍流強度和湍流積分尺度的同步控制。試驗風速采用Cobra Probe進行測量,采樣頻率為1 250 Hz,采樣時間30 s。來流風場縮尺比為1∶200。平均風速剖面與湍流度剖面根據工程科學數據庫(ESDU)[15-16]建議的對數率進行模擬,試驗值與理論值吻合良好。表1列出了各地貌的地表粗糙長度z0(對應對數率剖面)與粗糙度指數α(對應指數率剖面)。圖1給出了4種地貌風剖面的比較,圖中H為建筑頂部高度,和Iuz分別為高度z處的平均風速和湍流強度。由圖可見,由于粗糙元配置的不同,O類地貌與S類地貌的風剖面差異較大;由于擋板的存在,S1地貌較S2地貌的湍流強度有小幅提高,但O1地貌與O2地貌下的湍流度相差不大。

表1 各地貌下的風場參數Table 1 Parameters for the upstream flow under different terrains

圖1 平均風速剖面與湍流度剖面Fig. 1 Mean wind velocity and turbulence intensity profiles

來流風場的湍流積分尺度按ESDU 74[17]建議的經驗式進行模擬。湍流積分尺度試驗值按下式計算:

式中,Lu為湍流積分長度,Lt為湍流積分時間尺度,ρuu(τ)為脈動風速自相關系數,τ為時差,τ0為自相關系數收斂至0時對應的時差。各地貌下z= 0.762H(靠近駐點)高度處的湍流積分長度試驗值及ESDU 74結果可見表1,其中O1地貌下試驗值與ESDU值較為接近。O1、S1地貌分別較O2、S2地貌具有更大的湍流積分長度,可見擋板的存在能使湍流積分長度顯著增大,增幅在50% ~ 100%左右。

脈動風速譜按ESDU 74[17]建議的Von Karman譜進行模擬。圖2給出了試驗值與理論值的擬合情況,圖中f為頻率,Suu為脈動風速功率譜密度,由圖可見二者吻合良好。圖3給出了4種地貌風速譜的比較,圖中曲線在豎向的相對位置代表湍流度的相對大小,在水平向的相對位置代表湍流積分尺度的相對大小。由圖可見,由于S類地貌的湍流度總體大于O類地貌,故S類地貌曲線整體在O類地貌上方;由于1類地貌湍流積分尺度總體大于2類地貌,1類地貌曲線整體較2類地貌略微偏左。

圖2 O1地貌下脈動風速譜Fig. 2 Wind velocity spectrum under O1 terrain

圖3 各地貌脈動風速譜Fig. 3 Wind velocity spectra under different terrains

1.2 模型及試驗參數

本次模型為一矩形截面高層建筑剛性模型,縮尺比為1∶200。模型共由11段組成,如圖4所示,其中段1~段6的深寬比為0.5(按圖中風向),段7~段11的深寬比為1,通過拼接組合及調整擺放角度可獲得深寬比為0.125~8的模型工況。沿模型高度方向共布置了7個測點層,各層高度分別為0.1H、0.3H、0.5H、0.65H、0.8H、0.9H、0.98H,各層布置方式相同。深寬比為8的工況的模型參數如表2給出,平面測點布置如圖4所示。測壓孔通過測壓管道系統接入至掃描閥上,掃描閥最終通過2048通道全自動測壓掃描系統被接入至PC機上,關于測壓管道與掃描系統的詳細信息可參見文獻[18]。本次試驗對深寬比為1/8、1/5、1/3、1/2、1、2、3、5、8的共計9種深寬比工況進行了4種不同地貌下的同步風壓測試,風洞試驗照片如圖5所示。試驗風壓采樣頻率為400 Hz,采樣時長為90 s,相當于足尺1H以上。風洞內邊界層高度以上的參考風速為14.0 m/s。對于每個深寬比工況,進行了不同風向的測試,本文僅針對正交風向情況進行研究。

圖4 模型平面測壓點布置Fig. 4 Plan layout of pressure taps on the testing model

表2 D /B = 8工況模型參數Table 2 Model parameters for the configuration of D /B = 8

圖5 模型及風洞試驗照片Fig. 5 Photos of the models in the wind tunnel

1.3 數據處理

風壓系數按風工程界慣例,采用建筑頂部高度H作為歸一化高度,定義如下:

式中,pi為測點i的風壓值,p∞為靜壓值,ρ為空氣密度。

為考察建筑某一表面的整體風壓系數CpA或風荷載合力系數CF的大小,將其定義如下:

式中:Cpi為測點i的風壓系數,Ai為測點i的控制面積,A為進行計算的建筑表面的面積。由上式可見,整體風壓系數即為測壓點風壓系數在整個表面上的面積加權平均值。

極值分析方法對極值估計結果有很大影響。本文采用Gumbel擬合法進行極值估計,其中Gumbel分布參數應用Lieblein BLUE法確定,由文獻[19]可知,該方法具有較高的準確性與穩定性。關于該方法的細節可參見文獻[20-22]。

2 平均風荷載

2.1 迎風面

圖6給出了各地貌下建筑迎風面整體風壓系數的平均值隨深寬比的變化情況。為降低不同地貌風速剖面的差別產生的影響,該值由以測壓點高度作為歸一化高度的局部風壓系數計算獲得,以CpLA表示。由圖可見,各深寬比工況迎風面上的整體風壓系數差別較小,這是由于對于除邊緣附近的迎風面區域,氣動力受建筑本身引起的非定常流的作用很小,基本僅受來流的湍流特性影響。1類地貌下整體風壓系數的平均值較2類地貌要大,可見迎風面上平均風壓系數與湍流積分尺度呈正相關。湍流度對迎風面上平均風壓系數的影響規律不明確。

圖6 迎風面整體風壓系數平均值隨深寬比的變化情況Fig. 6 Variation of the mean area-averaged pressure coefficients with the depth-to-width ratio on the windward surface

2.2 側風面

Ruderich和Fernholz[23]通過對二維鈍體分離再附流下的平均風壓分布的研究,發現在均勻來流下,無論鈍體體型及雷諾數是否改變,在距鈍體前緣的流向距離d被平均再附長度Xr歸一化時,分離再附流下的約化風壓系數沿 流向的分布具有相同的分布曲線。其中約化風壓系數如下式定義:

上述研究基于的流場為均勻來流,即低湍流度來流。Akon和Kopp[14]基于低矮建筑頂部繞流的粒子圖像測速(PIV)試驗數據及現有文獻中針對二維矩形棱柱側邊繞流結果,給出了平均再附點處值 隨湍流度變化的線性擬合公式:

通過上式即可根據來流某高度處的的湍流度獲得平均再附點處的約化風壓系數值 ,再結合風洞試驗數據并根據式(4)獲得平均再附點處的平均風壓系數值,進而獲得平均再附長度。盡管鈍體體型可能對周圍繞流產生一定影響,但鑒于分離再附流及其下風壓分布特性的相似性,本文在估算高層建筑側邊繞流的平均再附長度時采用文獻[14]給出的方法是合理的。

圖7給出了D/B= 5、2、1的建筑在各地貌下側風面上的平均風壓系數分布,這三個工況分別對應再附、部分再附與近乎完全分離3種不同的繞流特性,由圖可以非常直觀地看到地貌對側風面風壓分布的影響。為了更準確地進行考察,圖8同時給出了平均和脈動風壓系數沿水平方向的分布曲線。此外,表3對平均流動再附的情況給出了按Akon和Kopp[14]方法估算的平均再附長度。

圖7 側風面平均風壓系數分布Fig. 7 Mean pressure coefficient distributions on the side surface

圖8 側風面z/H = 0.65處平均與脈動風壓系數分布Fig. 8 Mean and fluctuating pressure coefficient distributions at z/H = 0.65 on the side surface

表3 側風面z/H = 0.65處平均再附長度Table 3 Mean reattachment lengths at z/H = 0.65 on the side surface

對于D/B= 5和D/B= 2的建筑,即平均流動再附和部分再附(O1地貌下)的情況,S類地貌較O類地貌下分離泡內負壓更小(均指絕對值),風壓恢復更為提前。由表3可見,對D/B= 5的情況建筑,S1地貌下平均再附點出現位置較O1地貌下提前約30%。對于D/B= 2的建筑,隨著湍流度的增大,平均流動狀態已從O1地貌下的未再附,轉變為S1地貌下的再附。值得注意的是,在平均再附點處,風壓尚未完全恢復。這些規律與諸多針對二維平板鈍體的研究(如文獻[4])結果相似,但本文明確給出了三維棱柱尤其是側邊繞流部分再附情況的風壓分布。對于D/B= 1的建筑,即平均流動近乎完全分離的情況,與上述兩種情況相似,S類地貌下分離流下負壓更小;O類地貌下整個側風面上均呈現較大負壓,而S類地貌下在靠近尾緣處已經開始出現風壓恢復的趨勢。由以上可見,湍流度增大使分離流下平均負壓減小、風壓恢復提前、平均再附長度減小。另外,對于深寬比介于1到2的情況,即平均流動近乎完全分離及部分再附的情況,湍流度對平均風壓分布影響較大。

對于脈動風壓系數,由圖8可見,對于D/B= 5和D/B= 2的建筑,在分離泡內,S類地貌下的脈動風壓系數較O類地貌明顯要大,且曲線峰值出現位置較O類地貌明顯靠前。許多針對均勻流下二維棱柱周圍繞流的研究表明,最大脈動風壓出現于Kelvin-Helmholtz渦卷起、配對并分解成湍流之后,一般該位置剛好出現在平均再附點前[4,24]。而由圖8可見,O1地貌下最大脈動風壓出現在0.59Xr處,顯著提前于平均再附點,S1地貌下則更為提前,出現在0.44Xr處,該規律與文獻[4,24]發現相似。此外,在風壓恢復過程中,對于平均流動再附情況,S類地貌的平均風壓系數仍較O類地貌偏小,但二者之間的差距明顯減小,兩類地貌下脈動風壓系數的差距也明顯減小,分布曲線幾乎趨于一致。對于平均流動部分再附情況,兩類地貌下脈動風壓系數在建筑尾緣處均略微增大,O類地貌增幅較S類地貌略微偏大。對于D/B= 1的建筑,流動近乎完全分離,即整個側風面均處于分離泡下,故S類地貌下整個側風面上的脈動風壓系數均較O類地貌要大。由以上可見,湍流度增大會使分離再附流下脈動風壓增大,脈動風壓峰值出現位置更靠近前緣。

相較湍流度的影響,湍流積分尺度對平均風壓分布的影響較小,但仍對平均繞流有一定的效應。總體上2類地貌較1類地貌平均負壓略微偏小。在靠近前緣處,約0.33Xr(O1地貌)內,2類地貌則明顯較1類地貌要小。對于平均流動再附或部分再附的情況,在風壓恢復過程中,二者之間差距明顯減小;對于流動近乎完全分離的情況,在整個側風面上,1類地貌下的平均負壓均較2類地貌要大。對于D/B= 5情況風壓恢復過程及恢復后的平均風壓系數,S2地貌較S1地貌偏大可能是由于S2地貌下的湍流度較小導致的。由以上分析可知,湍流積分尺度減小使總體平均負壓減小、靠近前緣區域平均負壓明顯減小。

同時,由圖8可見,O2、S2地貌下的脈動風壓分布曲線形狀分別與O1、S1地貌相似,但前者脈動風壓系數值總體上均分別小于后者,說明湍流積分尺度減小會使總體脈動風壓減小,但對脈動風壓分布曲線形狀影響不大。

由圖3可見,在湍流度相同的情況下,湍流積分尺度較小的風場的脈動風速譜,通常在高頻區具有更多的能量分布,即含有較少的大尺度湍流成分、而含有更多的小尺度湍流成分。這些分析表明,大尺度湍流使分離流下平均風壓和脈動風壓值增大,而小尺度湍流促使分離流更早再附,且使分離流下的平均風壓和脈動風壓值減小、峰值出現位置更靠近前緣。

需要說明的是,上述分析大多基于側風面z/H=0.65高度處,在該高度處建筑側風面繞流狀態較為接近二維鈍體。而由圖7可見,建筑側風面風壓分布沿高度方向變化明顯,呈現明顯的三維效應。各地貌下側風面最大平均負壓均出現在建筑前緣頂角處。

2.3 背風面

背風面上平均風壓分布較為非常均勻,故本節僅給出了各地貌下背風面整體風壓系數隨深寬比的變化情況,如圖9所示。由圖可見,背風面上的整體風壓基本均在深寬比為0.5至1時達到最不利。S類地貌下背風面上的平均風壓系數均小于O類地貌,對于平均流動再附的情況及在深寬比非常小時,兩類地貌下的差距較小。但對于流動近乎完全分離及部分再附的情況該規律非常明顯,該深寬比范圍內的建筑的湍流度效應較大。1類與2類地貌下的平均風壓系數差距較小,通常1類地貌略微大于2類地貌。由此可見,背風面上的負壓大小與湍流度通常呈負相關,而與湍流積分尺度呈正相關。

圖9 背風面整體風壓系數平均值隨深寬比的變化情況Fig. 9 Variation of the mean area-averaged pressure coefficients with the depth-to-width ratio on the leeward surface

3 極值風荷載

3.1 迎風面

迎風面上,各深寬比建筑迎風面上的極值風壓分布非常相似,故圖10僅給出D/B= 1工況迎風面豎直中心線處極大值風壓系數分布。由圖可見,1類地貌、S類地貌極值風壓系數分別較2類地貌、O類地貌要大。可見,迎風面上的極值正風壓系數是由湍流度和湍流積分尺度綜合決定的,其中極值風壓系數與湍流積分尺度明顯呈正相關關系。

圖10 D/B = 1工況迎風面豎直中心線處極大值風壓系數分布Fig. 10 Maximum pressure coefficient distributions along the windward vertical centerline for D/B = 1

3.2 背風面

背風面上的極值風壓分布較為均勻,湍流特性對極值風壓的影響規律也較為簡單。由圖11給出的豎直中心線上的極小值風壓系數分布曲線容易看到,1類地貌下的極值負壓普遍較2類地貌要大;但S類地貌與O類地貌下的差別未呈現明顯的規律。這可能是由于S類地貌來流湍流脈動較大、但因分離流更早在側邊再附而使尾流有所減弱這兩種因素相互作用導致的。

圖11 背風面豎直中心線上的極小值風壓系數分布Fig. 11 Minimum pressure coefficient distributions along the leeward vertical centerline

3.3 側風面

圖12 側風面z/H = 0.65處極小值風壓系數分布Fig. 12 Minimum pressure coefficient distributions at z/H = 0.65 on the side surface

圖12給出了側風面z/H= 0.65處的極小值風壓系數分布曲線。由圖可見,對于D/B= 5和D/B= 2的建筑,S1類地貌分離泡內的極值負壓較O類地貌明顯要大,尤其是在靠近前緣位置處,這與湍流度對平均風壓分布的影響規律相反。但相似的是,S類地貌下最大極值負壓出現的位置較O類地貌更為提前,但最大極值負壓一般均出現在0.2Xr~0.3Xr處,這與文獻[6]的結果基本相符。對于平均流動完全再附的情況,極值負壓在達到最大值后一般隨到前緣距離的增大而減小;但對于平均流動部分再附或近乎完全分離情況,極值負壓則在尾緣突然顯著增大,不同地貌對尾緣處較大負壓值的影響規律不明顯。由此可見,湍流度與分離再附流下的極值負壓呈正相關關系。

湍流積分尺度對側風面極值風壓分布也有一定影響。由圖12可見,1類地貌下極值風壓系數總體上較2類地貌下要大,說明湍流積分尺度減小會使總體極值負壓減小。對于湍流積分尺度對分離再附流下的最大極值負壓的影響,考慮到湍流積分尺度較大的風場通常在低頻區具有更多的能量分布,即含有較多的大尺度湍流成分,可以發現更大尺度的湍流可能是產生更大的極值負壓峰值的原因。故本文也部分驗證了文獻[6]中提到的“最大極值負壓由分離剪切流在鈍體表面附近卷起而形成的較大旋渦引起”。

圖13 橫風向氣動力譜Fig. 13 Cross-wind aerodynamic force spectra

4 橫風向氣動力譜

圖13給出了深寬比D/B= 0.33、1、2和5的建筑在4種地貌下的橫風向氣動力譜的比較。對于D/B=1工況,盡管S類地貌來流脈動風速譜較O類地貌要高,但對于氣動力譜,O類地貌下譜峰所含能量較S類地貌要大,這可能是由于S類地貌湍流度更大、分離流已開始部分再附于建筑側面導致的。相比之下,D/B= 0.33工況的建筑尾體較短,氣流在各地貌下均完全分離并脫落于尾流中,S類地貌下氣動力譜則整體高于O類地貌。對于D/B= 2工況,各地貌下分離流均開始再附,橫風向氣動力譜則呈現更寬的譜峰。對于D/B= 5工況,各地貌下分離流均近乎完全再附,譜在高頻區的峰消失。總體而言,湍流度對深寬比介于1到2的建筑、即流動近乎完全分離及部分再附的情況影響較大。

圖14給出了4種地貌下建筑斯托羅哈數St隨深寬比的變化曲線及與文獻結果的比較。文獻[25]數據較本文數據明顯要小,且其給出了深寬比直至4的斯托羅哈數結果,而本文對于深寬比大于2的工況(S類地貌下則對于深寬比大于1的工況)的斯托羅哈數因譜峰過寬已很難由譜峰位置確定出來。但是對于深寬比為1的工況,本文數據與文獻[26]非常接近。由圖14可見,整體上各地貌下的St相差不大。

圖14 斯托羅哈數St隨深寬比變化情況Fig. 14 Strouhal number variation with the depth-to-width ratio

5 結 論

本文對4種不同風場下9種深寬比的矩形截面高層建筑各立面的平均、脈動、極值風壓與橫風向氣動力譜進行了考察,獲得了如下結論:

1)隨湍流度增大,分離流更早再附于側風面上,其中郊區地貌下較開闊地貌再附提前約30%。湍流度增大使分離流下平均負壓(均指絕對值)減小、極值負壓增大,最不利脈動風壓與極值負壓出現位置較再附點更為提前。隨湍流度增大,背風面上平均負壓減小,但極值負壓變化不明顯。對于方形建筑,旋渦脫落頻率處的橫風向氣動力譜能量隨湍流度增大有所降低。

2)隨湍流積分尺度減小,迎風面上、背風面的平均風壓與極值風壓絕對值均減小,側風面分離流下平均負壓、脈動風壓與極值負壓均減小。湍流積分尺度對分離再附流形態及分離流下的風壓分布形狀影響不大。

3)不同深寬比建筑分離再附流下的平均、脈動與極值風壓分布形狀相似,但背風面平均風壓值以及橫風向氣動力譜形狀差別較大。總體而言,湍流特性對深寬比介于1到2的建筑,也即流動近乎完全分離及部分再附的情況影響較大。

4)大尺度湍流使分離流下平均風壓和脈動風壓值增大,而小尺度湍流促使分離流更早再附,且使分離流下的平均風壓和脈動風壓值減小,峰值出現位置更靠近前緣。更大尺度的湍流可能是產生更大極值負壓峰值的原因。

主站蜘蛛池模板: 日本免费高清一区| 午夜一区二区三区| 国产高清色视频免费看的网址| 在线观看欧美国产| 无遮挡国产高潮视频免费观看| 波多野结衣在线一区二区| 国产成人久久综合777777麻豆| 国产精品久久久精品三级| 国产亚洲精品91| 丁香综合在线| 毛片久久网站小视频| 国产最新无码专区在线| 国产黄在线观看| 亚洲欧美在线综合一区二区三区 | 色爽网免费视频| 爱爱影院18禁免费| 青青青亚洲精品国产| 伊人精品视频免费在线| 亚洲成人精品久久| 国产理论最新国产精品视频| 免费 国产 无码久久久| 日韩欧美国产成人| 久久综合激情网| 色婷婷视频在线| 国产经典免费播放视频| 色亚洲成人| 中文字幕佐山爱一区二区免费| 午夜日b视频| 成人在线视频一区| 色屁屁一区二区三区视频国产| 久久综合丝袜长腿丝袜| 国产人成网线在线播放va| 无码一区二区波多野结衣播放搜索| 黄色片中文字幕| 999精品在线视频| 欧美日一级片| 四虎精品黑人视频| 青青国产成人免费精品视频| 亚洲国产成人精品一二区| 国产乱子伦精品视频| 国产欧美专区在线观看| 国产一区二区福利| 真人免费一级毛片一区二区| 免费在线色| 日韩经典精品无码一区二区| 久久影院一区二区h| 91精品免费高清在线| 国产精品自拍露脸视频 | 国产办公室秘书无码精品| 无码在线激情片| 国产一级毛片yw| 54pao国产成人免费视频| 日本欧美视频在线观看| 日本高清视频在线www色| 日本精品视频一区二区| 久久国产香蕉| 日韩精品久久久久久久电影蜜臀| av午夜福利一片免费看| 日韩不卡免费视频| 国产精品大白天新婚身材| 手机精品视频在线观看免费| 一级毛片免费观看久| 波多野结衣一区二区三区四区视频 | 5555国产在线观看| 美女被躁出白浆视频播放| 特级做a爰片毛片免费69| 成人国产免费| 免费三A级毛片视频| 亚洲综合香蕉| 国产97公开成人免费视频| 日本精品αv中文字幕| 伊人久综合| V一区无码内射国产| 视频在线观看一区二区| 尤物精品视频一区二区三区| 国产丝袜一区二区三区视频免下载| 熟妇丰满人妻| 国产成人8x视频一区二区| www成人国产在线观看网站| 亚洲国产成熟视频在线多多 | 不卡视频国产| 曰韩免费无码AV一区二区|