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

        ?

        山地超高層建筑風致響應研究

        2011-01-25 00:44:16李正良魏奇科黃漢杰
        振動與沖擊 2011年5期
        關鍵詞:風洞試驗風場山體

        李正良,魏奇科,黃漢杰,孫 毅

        (1.重慶大學 土木工程學院,重慶 400045;2.中國空氣動力研究與發(fā)展中心低速空氣動力研究所,四川 綿陽 621000)

        山地超高層建筑風致響應研究

        李正良1,魏奇科1,黃漢杰2,孫 毅1

        (1.重慶大學 土木工程學院,重慶 400045;2.中國空氣動力研究與發(fā)展中心低速空氣動力研究所,四川 綿陽 621000)

        風在山地地形的干擾下,其幅值和空間分布規(guī)律相對平地均會發(fā)生較大改變,尤其是山頂處風速有明顯增大。如果不考慮山地的影響,仍然用平地邊界層風場進行高層建筑維護結構和整體結構計算將偏于不安全。以前對山地風場的研究多限于平均風的變化,而對于脈動風的湍流特征和頻域特性較少提及,且缺少山地與平地風場中超高層建筑風致響應計算的對比。本文建立了鐘形、高斯形、余弦型幾種軸對稱三維山體模型,通過數值計算得到了不同山體坡度下山頂平均風加速比。計算結果顯示,山頂平均風有較大的增大效應,最大加速比可達1.7,且反函數擬合結果與模擬結果更加吻合。通過與風洞試驗結果的對比可見,平均風的數值模擬有較高準確性,湍流與試驗相比還有一定差別。通過某超高層建筑的風振響應分析表明,山地風的增大效應對超高層建筑整體響應計算不可忽略,位移響應增大比例最大可達20%。

        山地風場;CFD;風洞試驗;超高層建筑;風致響應

        大氣邊界層在山地地形的影響下,其幅值和空間分布規(guī)律均會發(fā)生較大改變。很多學者針對山地風場進行了風洞試驗、數值模擬以及現場實測研究。Taylor[1,2]總結了早期的一些研究成果,并運用混合譜有限差分法(MSFD)及非線性混合譜有限差分法(NLMS-FD)對山地風流動進行了計算,結合風洞試驗結果提出了山地風計算模型;Miller,Davenport[3]進行了一系列二維復雜山體的風洞試驗,并與加拿大、英國規(guī)范和MSFD模型進行了對比,得出了規(guī)范偏保守的結論;Weng[4]基于MSFD和NLMSFD的計算結果,考慮了表面粗糙度和風速增大系數的非線性,提出了新的風場模型;Kim[5,6]對四處典型山地地形進行了數值模擬,選用了標準k-ε模型和RNG k-ε模型兩種湍流模型,通過與現場實測數據的對比,認為數值模擬可以作為復雜山地風場預測可靠的工具;Carpenter[7]進行了1:1 000比例的二維山地風洞試驗,包括坡度較小和較大的正弦截面山體、多個連續(xù)山體以及不規(guī)則山體;T.Takahashi[8]利用熱線風速儀測量了兩種不同表面粗糙度下二維山體近地平均風剖面、湍流度、雷諾應力、湍動能和湍流粘性剖面,分析了不同粗糙表面山體的湍流特性;同時,一些國家和地區(qū)規(guī)范按地形條件提供了地形影響系數以供參考。

        但已有研究均只針對山地平均風,對于山地風的湍流特征和頻域特性研究較少,對山地風作用下超高層建筑的響應研究也較少。

        1 山地平均風加速比

        1.1 加速比模型

        加速效應(speed-up effect)指在山地地形中,某高度平均風速比平地相應高度平均風速有所增加的效應,一般在山頂的近地面最為明顯。通常用一個無量綱的參數:加速比(speed-up ratio)來定量描述加速效應:

        式中:U(z)表示山地地面以上z高度處的風速,U0(z)表示平地地面以上z高度處的風速。

        針對山頂加速比,Jackson和Hunt[9]最先提出了一種計算二維光滑山體加速效應的解析算法,但這種算法僅僅適用于未發(fā)生流動分離的情況;隨后Jackson[10]運用風洞試驗和數值模擬檢驗了這種算法;Hunt[11]又對其提出了修改,使之更加趨于完善。Mason和Sykes[12]在Jackson等人的基礎上,將其算法延伸到單個三維山體的情況。Kaimal和Finnigan[13]通過研究指出,以Jackson和Hunt的方法為主體的算法,在計算未發(fā)生流動分離時的最大加速效應誤差可以小于15%到18%。

        Jackson和Hunt[9]提出的算法中,對于最大加速比的規(guī)定為:

        Taylor和 Lee[14]的在 Jackson算法基礎上提出了“原始算法”,可以計算山頂不同高度處的加速比:

        圖1 山體形狀及加速效應示意圖Fig.1 Schematic diagram of a hill and the speed-up effect

        1.2 一些國家和地區(qū)風荷載規(guī)范關于地形修正系數的建議

        1.2.1 中國規(guī)范[15]

        風壓高度變化系數考慮地形條件的影響系數在山頂處可用以下公式計算:

        tgα是山峰或山坡在迎風面一側的坡度;κ是系數,對山峰取3.2,對山坡取1.4;H是山峰或山坡高度;z是建筑物計算位置離地面的高度(當z>2.5H時,取z=2.5H)。對于山峰和山坡的其他部位,影響系數按η的線性插值確定。

        1.2.2 美國最小荷載規(guī)范(ASCE)[16]

        風速地形增大因子Kzt,其表達式如下:

        表1 美國規(guī)范中地形影響系數各參數的取值Tab.1 Coefficients for Kztin ASCE

        美國規(guī)范總結了 Taylor[1,2]和 Weng[4]的研究成果,綜合考慮了地貌類別,山體形狀,山體坡度等的影響。

        1.2.3 加拿大規(guī)范及ISO標準

        加拿大規(guī)范[17]及 ISO 標準[18]的思想同美國規(guī)范基本相同,僅K1的系數取值不同:

        (1)加拿大規(guī)范不考慮地貌類別差異,山脊、陡坡和小山各地形 K1的系數分別為2.2、1.3 和1.6,比美國規(guī)范各類地貌K1的系數都大;

        (2)ISO標準也不考慮地貌類別差異,山脊、陡坡和小山各地形 K1的系數分別為 2.0、1.8 和1.6。

        1.2.4 日本規(guī)(AIJ)[19]

        地形影響系數可由式(6)確定。

        1.2.5 澳大利亞規(guī)范(AS/NZS)[20]

        ① 當 H/2Lu<0.05 時,Mh=1.0;

        ② 當0.05≤H/2Lu≤0.45 時,

        ③ 當H/2Lu>0.45時,山頂近地面的流動分離區(qū)

        H是小山、山脊或陡坡的高度;Lu是小山、山脊或陡坡頂至半山高的水平距離;x是建筑物上風向或下風向至山頂的水平距離;L1取0.36Lu或0.4H的較大值;L2除了陡坡的下風向取10L1外,其他均取4L1;Z是離地面高度。

        1.2.6 不同規(guī)范的對比分析

        一個三維山體,高150 m,半山高到山頂的距離300 m,坡度h/L1=0.5。根據各規(guī)范計算的山頂平均風增大系數見圖2。

        可見,根據各參考模型計算的山頂位置平均風增大系數相差較大。中國規(guī)范計算值最大,美國規(guī)范次之,日本和澳大利亞規(guī)范接近。并由此說明,山地風場平均風的加速效應仍然沒有一個相對一致的研究結果。

        2 山地風場特性數值模擬

        2.1 幾何建模及網格劃分

        本文采用Fluent建立模型進行數值模擬。計算域在豎直方向取為10 H,模型上游10 H,下游20 H,寬度20 H,模型采用六面體映射網格劃分計算區(qū)域。山體近壁網格劃分較密,取網格豎向尺寸為H/150,遠離壁面網格尺寸為H/2。山體的幾何模型及網格劃分如圖3所示。

        2.2 邊界條件設定和求解設置

        入流面選用速度進口邊界。速度剖面選用中國規(guī)范B類場地,取α=0.16。湍流度剖面和湍流積分尺度參考日本規(guī)范。

        在入流處直接給定湍動能k和湍流耗散率ε:

        出流面采用自由出流邊界,流域頂部和兩側采用對稱邊界,等價于自由滑移壁面邊界。地面采用無滑移的壁面條件,且在地面引入粗糙壁面修正[21]。由于鈍體繞流將出現分離、再附、沖撞、環(huán)繞及旋渦等復雜的流動,所以選用非平衡壁面函數來模擬壁面附近復雜的流動現象。

        壓力-速度求解選用耦合算法求解,同時選用二階迎風格式離散控制方程,以減少假擴散誤差。湍流模型采用RNG k-ε模型,可較好的處理高應變率及流線彎曲程度較大的流動[22]。

        2.3 山體形狀影響

        常用的山體形狀模型包括鐘型、高斯型和余弦型。

        式中:r為到山頂的水平距離,h是山體高度,L是山頂到半山高度處的水平距離。圖4是不同山體形狀加速比計算結果。

        可見,山體形狀對山頂平均風加速比影響不大。因此本次模擬和相關試驗就采用同一種山體形狀——余弦型山體。

        表2 不同山體形狀山頂加速比對比Tab.2 The comparison of speed-up ratio at the hilltop of different hill shapes

        2.4 山頂加速比沿高度的分布

        綜合前面的分析,可以看出加速比沿高度的分布可以歸納為兩種形式:一種是以美國規(guī)范為代表的指數函數形式;一種是以澳大利亞規(guī)范為代表的反函數形式。取較緩和較陡兩種坡度,分別運用公式(14a)和(14b)對山頂平均風加速比數值計算的結果進行擬合,如圖5所示。

        式中ΔSmax反映的山頂表面的最大加速比,而k反映的是加速比的影響高度。

        圖4 數值模擬結果的擬合Fig.4 The fitting of the numerical simulation result

        可見,不管坡度是緩還是陡,用反函數的形式進行擬合的結果與數值模擬吻合更好。因此,本文都采用反函數的形式進行擬合。

        2.5 山體坡度影響

        由澳大利亞規(guī)范可得:

        可以看出澳大利亞規(guī)范中ΔSmax隨山體坡度線性增加,一直到坡度為0.9時趨于0.714,而k與山體坡度呈反函數關系,最后趨于0.4 H。

        對不同坡度下山頂加速比數值模擬數據進行參數ΔSmax和k二次擬合的結果可見圖5、圖6。

        其中坡度β=H/LH。

        由圖6可見,參數k用式(16b)擬合的結果與計算結果非常吻合,但與澳大利亞規(guī)范相比,明顯要小。

        2.6 湍流特征

        湍流的強弱主要由湍動能來衡量。類比于平均風速,山地風場湍動能的改變也可以用湍動能增大系數來表示,如公式(17)所示。

        由圖6可見,與澳大利亞規(guī)范相比,參數ΔSmax數值模擬和擬合結果比規(guī)范要大,同時曲線形狀也有所不同。規(guī)范采用的是折線形模型,而計算結果可以用經過原點的三次多項式如式(16a)擬合,擬合效果很好。

        式中:k(z)表示山地地面以上z高度處的湍動能,k0(z)表示平地地面以上z高度處的湍動能。

        由圖可以看出,當坡度較小時,山頂近地湍動能比來流湍動能略小,分析原因可能是數值模擬大氣邊界層的不平衡原因引起的。對于這個問題,很多學者進行了研究[23],但是有效的保持大氣邊界層水平平衡性,尤其是近地表水平平衡性的方法還沒有得到。

        隨著坡度增大,底部湍動能逐漸增大,但增大并不明顯。而邊界層上部湍動能卻略有減小,分析原因應該為上部邊界層基本沒有受到山地的干擾,因此如果以絕對高度來比較,應該是相等的,但是增大系數是以相對于地表的高度來比較的,而湍流度隨高度越高而約小,因此上部邊界層湍動能增大系數還略有減小??偟膩碚f,山頂湍動能的變化不大,沒有平均風速的增大效應明顯。

        3 山地風場風洞試驗對比

        試驗在中國空氣動力研究與發(fā)展中心低速所1.4 m×1.4 m風洞中進行。幾何模型縮尺比為1∶1 000,風速縮尺比約為1∶2.5,測試儀器采用排管和電子掃描閥。對大氣邊界層進行了細致的調試,評價指標包括平均風速剖面、湍流度剖面、脈動風功率譜、各指標均滿足要求。試驗湍流積分尺度結果約為0.2m,換算到實際尺度約200 m,與大氣邊界層的湍流積分尺度基本一致。

        山體坡度分別為0.5和1.0時,山頂平均風加速比數值模擬和試驗的對比見圖8。

        圖7 湍動能增大系數Fig.7 Increasing coefficient of turbulent kinetic energy

        圖8 平均風加速比數值模擬和風洞試驗對比Fig.8 Comparison of the speed-up ratio

        由圖8可知,不管較緩的山體還是陡峭的山體,數值模擬和風洞試驗的結果都非常接近,數值模擬的結果略大,加速比沿高度的分布基本相同。數值模擬在平均風速風壓的模擬方面已經比較準確。本文擬合的山地平均風加速比公式為對山地平均風取值提供了參考。

        但是從許多學者進行數值模擬的經驗來看,風場湍流模擬方面與試驗相比還有一定差距。以下列出本文數值模擬與風洞試驗湍動能的對比,如圖9。

        由湍動能增大系數對比可以看出,數值模擬和風洞試驗結果基本都分布在-0.2~0.2之間,在邊界層上部比較吻合,但是在近地部分有較大差別。數值模擬比風洞試驗在近地高處湍動能要高,估計與模擬采用的湍流模型有關,兩參數的高雷諾數湍流模型都過大的估計了近地處的湍動能。因此,數值模擬的結果比試驗要高。

        由于數值模擬采用的定常計算,并不能反映風場的頻域特性,因此山地風場的頻域特性只能通過風洞試驗得到。為此將風洞試驗中,山頂位置脈動風速功率譜與前方來流參考脈動風速功率譜進行對比如圖10。

        由圖可見,山頂位置縱向脈動風速功率譜與來流脈動風速基本一致。山頂位置與平地相比,頻域特征基本相同。

        4 山地風場下超高層建筑風振響應分析

        某超高層建筑總高300 m,特征寬度50 m,高寬比為6.0,單位高度質量5×105kg/m。由于超高層建筑各頻率間隔較大,振型間的耦合影響較小,同時高階振型影響較小,因此本文僅計算超高層建筑的一階位移響應。一階振型取為線性振型,一階頻率0.15 Hz,周期約為7 s,阻尼比1.5%,基本風壓取為0.55 kPa。通過抖振頻域分析方法計算建筑的風振響應,以得到山地對超高層建筑風振響應的影響范圍

        由數值模擬和風洞試驗的結果可知,山頂位置風場的湍流特性和頻域特性變化不大,因此本文平地和山頂位置的風速譜均直接選用廣泛應用的Davenport譜,見公式(19);相干函數采用指數形式,如公式(20),取Cy=7。

        則一階廣義力譜可由各層荷載互譜乘以振型的平方在空間上積分得到,如公式(21);

        一階廣義力譜乘以傳遞函數就得到一階廣義位移譜,然后再乘以振型的平方就可以得到各層位移譜,如公式(22);

        對各層位移譜頻域積分可得各層位移根方差值,如公式(23);再乘以峰值系數3.0,疊加上平均靜力響應,即可得到超高層建筑的總響應。

        表3 平地與山地高層1建筑頂部位移響應對比Tab.3 Displacement response of building in flat terrain and in hilly terrain

        平地超高層建筑和坡度為0.5的山頂位置超高層建筑風振響應對比見表3,由表3可見,當超高層建筑位于坡度為0.5的山地的山頂時,超高層建筑的平均位移響應、動力位移響應以及總位移響應都比相應的平地要大。其中平均位移響應增大約28%,動力位移響應增大約20%,總的位移響應增大24%。

        雖然山體坡度越大,近地風場平均風加速比越大,但是影響的高度也越小。因此坡度越大,對超高層建筑的影響并不一定一直增大,會存在一個極值坡度。與平地情況相比,將不同坡度下超高層建筑的頂部位移風振響應放大系數計算如圖11所示。

        由圖11可見,對于算例特定的超高層建筑而言,按美國規(guī)范的公式,超高層建筑的頂部位移響應在坡度為0.2時達到最大,與平地相比頂部位移放大系數達到1.3。按澳大利亞規(guī)范計算頂部位移放大系數隨坡度增大而增大,直到坡度為0.9時趨于平穩(wěn),最大達1.25。按數值模擬結果計算的頂部位移響應在山體坡度為0.7時達到最大,頂部位移放大系數1.2。如果把澳大利亞規(guī)范在坡度大于0.9的部分可以看作是保守的歸并,那么數值模擬的極值坡度正好在美國規(guī)范和澳大利亞規(guī)范之間。

        圖11 不同坡度建筑頂部位移放大系數Fig.11 Displacement increasing coefficient by different slope

        5 結論

        由數值模擬的分析,以及和風洞試驗的對比可以得出以下結論:

        (1)山地對平均風有較大的增大效應,山頂位置最大加速比可達1.7,對維護結構的計算要嚴加考慮。

        (2)山頂平均風最大加速比先增大后減小,而影響高度系數隨坡度增大呈反函數減小。本文擬合的山地平均風加速比公式為山地平均風取值提供了參考。

        (3)數值模擬和風洞試驗的對比可以看出,平均風的模擬有較高準確性,湍流模擬與試驗結果相比還有一定差距。

        (4)通過某超高層建筑的風振響應分析表明,山地風的增大效應對結構整體響應計算不可忽略,位移響應增大比例最大可達20%。

        [1]Taylor P A,Mason P J,Bradley E F.Boundary-layer flow over low hills-A review[J].Boundary-Layer Met,1987,39:107-132.

        [2]Taylor P A.Turbulent boundary-layer flow over low and moderate slope hills[J].Journal of Wind Engineering and Industrial Aerodynamics,1998,74 -76;25 -47.

        [3] Miller C A,Davenport A G.Guidelines for the calculation of wind speed-ups in complex terrain.JournalofWind Engineering and Industrial Aerodynamics 74-76(1998)189-197.

        [4]Weng W,Taylor P A,Walmsley J L.Guidelines for airflow over complex terrain:model developments.Journal of Wind Engineering and Industrial Aerodynamics, 2000, 86:169-186.

        [5] Kim H G,Patela V C,Leeb C M.Numerical simulation of wind flow over hilly terrain.Journal of Wind Engineering and Industrial Aerodynamics,2000,87:45 -60.

        [6] Kim H G,Lee C M,Lim H C,et al.An experimental and numerical study on the flow over two-dimensional hills.Journal of Wind Engineering and Industrial Aerodynamics,1997 ,66:17-33.

        [7]Carpenter P,Locke N .Investigation of wind speeds over multiple two-dimensional hills.Journal of Wind Engineering and Industrial Aerodynamics,1999,83:109 -120.

        [8] Takahashi T,Ohtsu T,Yassin M F,et al.Turbulence characteristics of wind over a hill with a rough surface.Journal of Wind Engineering and Industrial Aerodynamics,2002,90:1697-1706.

        [9]Jackson P S,Hunt J C R.Turbulent wind flow over a low hill[J].Q.J.Roy.Meteorol.Soc,1975,101:929 -955.

        [10] Jackson P S.The influence of local terrain features on the site selection for wind energy generating systems[R].Boundary Layer Wind Tunnel Laboratory Internal Report,University of Western Ontario,BLWT -1 -1979,1979.

        [11] Hunt J C R,Leibovich S,Richards K J.Turbulent shear flow over low hills[J].Q.J.Roy.Meteorol.Soc,1988,114:1435-1470.

        [12] Mason P,Sykes R.Flow over an isolated hill of moderate slope[J].Q.J.Roy.Meteorol.Soc,1979,105:383 -395.

        [13] Kaimal J,Finnigan J.Atmospheric Boundary Layer Flows-Their Structure and Measurement[M].Oxford,Oxford University Press,Oxford,1994.

        [14] Taylor P A,Lee R J.Simple guidelines for estimating wind speed variations due to small scale topographic features[J].Climatol.Bull,1984,18(2):3 -22.

        [15]中華人民共和國國家規(guī)范.建筑結構荷載規(guī)范(GB50009-2006)[S].北京:中國建筑工業(yè)出版社,2006.

        [16] ANSI/ASCE Standard,Minimum Design Loads for Buildings and Other Structures,2002[S].

        [17] National Commission of Building and Fire Code[S],National Building Code of Canada 1995(NBC),National Research Council of Canada,1995.

        [18] International Standard ISO4354.Wind Actions on Structures[s],2004.

        [19] AIJ Recommendations for Loads on Buildings.Architectural Institute of Japan,2004[s].

        [20] AS/NZS 1170.2:2002.Australian/New Zealand Standard,Structural design actions,Part 2:Wind actions[s].

        [21] Blocken B,Stathopoulos T,Carmeliet J.CFD simulation of the atmospheric boundary layer:wall function problems.Atmos.Environ,2007,41(2),238 -252.

        [22] Fluent Inc.Fluent 6.3 User’Guide.Fluent Inc.,Lebanon,2006.

        [23] Richards P J,Hoxey R P.Appropriate boundary conditions for computational wind engineering models using the k-e turbulence model.J. Wind Eng.Ind. Aerodyn,1993,46&47,145-153.

        Wind-induced response of super tall buildings in hilly terrain

        LI Zheng-liang1,WEI Qi-ke1,HUANG Han-jie2,SUN Yi1

        (1.College of Civil Engineering,Chongqing University,Chongqing 400045,China;2.China Aerodynamics Research& Development Center,Mianyang 621000,China)

        Considerable changes of amplitude and spatial distribution of wind velocity take place under interfering of hilly terrain.Wind load and response of super tall buildings at special locations in hilly terrain are significantly different from those in flat terrain,therefore,the current methods for calculating wind-induced response of super tall buildings in flat terrain may be unsafe under condition of hilly terrain.Previous studies were mainly focused on mean wind speed-up effect,and very few studies on turbulent characteristic of pulsive wind are involved with this subject.Besides,there is no comparison between results of wind induced response of super tall buildings in hilly terrain and those in flat terrain.Herein,several kinds of 3D hill model including bell,Gaussian and cosine shaped ones were established with CFD software Fluent.Vertical profiles of mean wind velocity at crest were calculated under different hill slopes.The numerical simulation showed that the great speed-up effect is found at crest and the maximum speed-up ratio reaches 1.7;the results with those using inverse functions;compared to the test data,the simulation of mean velocity is accurate,but the numerical simulation of turbulent has a certain difference.A numerical example of wind-induced response of a super-tall building showed that the effect of hilly terrain can't be ignored,and its displacement response in hilly terrain increases 20%compared with that in flat terrain.

        hilly terrain wind field;CFD;wind tunnel test;super tall buildings;wind-induced response

        TU973+.31

        A

        國家自然科學基金重大研究計劃(90715024)

        2010-02-05 修改稿收到日期:2010-03-01

        李正良 男,教授,博士生導師,1963年生

        猜你喜歡
        風洞試驗風場山體
        基于FLUENT的下?lián)舯┝魅S風場建模
        濟南市山體修復中的植物應用與技術——以濟南市臥虎山山體公園為例
        低風壓架空導線的風洞試驗
        電線電纜(2017年5期)2017-10-18 00:52:03
        “最美風場”的贏利法則
        能源(2017年8期)2017-10-18 00:47:39
        滾轉機動載荷減緩風洞試驗
        側向風場中無人機的飛行研究
        山體別墅設計分析
        人間(2015年21期)2015-03-11 15:23:42
        城市山體保護規(guī)劃的實踐與探索——以武漢市新洲區(qū)山體保護規(guī)劃為例
        遮擋條件下超高層建筑風洞試驗研究
        重慶建筑(2014年12期)2014-07-24 14:00:32
        高速鐵路接觸線覆冰后氣動力特性的風洞試驗研究
        精品国产一级毛片大全| 午夜精品久久久久久久| 中文无码一区二区不卡av| 国产亚洲美女精品久久久| 亚洲成a人网站在线看| 国产一级黄色片一区二区| 强开小婷嫩苞又嫩又紧视频| 日本乱码一区二区三区在线观看| 亚洲av无码一区二区三区网址 | 午夜国产视频一区二区三区| 亚洲av无码乱码国产精品| 亚洲AV无码一区二区三区人| 成人av天堂一区二区| 亚洲av毛片在线网站| 亚洲av福利院在线观看| 柠檬福利第一导航在线| 国产免费av片在线观看播放| 一区二区三区视频免费观看在线 | 免费a级毛片无码a∨中文字幕下载| 一本大道东京热无码| 国产视频嗯啊啊啊| 国产午夜精品视频在线观看| 久久久久成人精品无码中文字幕| 1000部精品久久久久久久久| 久久精品国产亚洲AV高清wy| 人妖一区二区三区视频| 性无码免费一区二区三区在线| 日本午夜福利| 少妇激情一区二区三区| 午夜dv内射一区二区| 女人扒开下面无遮挡| 国产亚洲精品不卡在线| 99久久婷婷国产精品综合| 精品无码无人网站免费视频| 国产人成精品综合欧美成人 | 国产午夜在线观看视频播放| 成人综合亚洲欧美一区h| 国产69精品麻豆久久| 亚洲字幕av一区二区三区四区| 亚洲无毛片| 97超碰中文字幕久久|