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

        ?

        典型山地地形豎向風(fēng)速分布特征

        2018-11-05 08:04:24樓文娟梁洪超李正昊章李剛
        關(guān)鍵詞:山脊來流山脈

        樓文娟, 梁洪超, 李正昊, 章李剛, 卞 榮

        (1. 浙江大學(xué) 結(jié)構(gòu)工程研究所, 浙江 杭州 310058;2. 國網(wǎng)浙江省電力公司經(jīng)濟(jì)技術(shù)研究院, 浙江 杭州 310000)

        0 引 言

        豎向風(fēng)速對敞開的水平向結(jié)構(gòu)和輸電線路風(fēng)偏十分不利,然而在目前的工程設(shè)計(jì)中,除大跨徑橋梁顫振穩(wěn)定性分析中考慮一定的風(fēng)攻角外,一般的土木工程結(jié)構(gòu)均忽略了上升氣流的不利影響。在國外,對山地地形風(fēng)場的研究自20世紀(jì)80年代便已經(jīng)開始。Jackson 與 Hunt[1]首先對二維平滑低矮山地的風(fēng)速分布進(jìn)行了理論分析。由于其研究模型坡度平緩,山體風(fēng)場中豎向風(fēng)速分量較小,該部分豎向風(fēng)速僅僅作為水平向風(fēng)速的擾動(dòng)進(jìn)行了考慮。Jackson與Hunt這一研究基本確立了以風(fēng)速加速比為主要參數(shù)的山地風(fēng)場的研究方法。Taylor等[2]對上述理論計(jì)算公式進(jìn)行修正,使當(dāng)時(shí)的理論在緩坡條件下計(jì)算不同高度的風(fēng)速加速比時(shí)更簡便。Weng等[3]在計(jì)算過程中引入了粗糙度影響,使理論計(jì)算更加成熟。截止90年代,國外已經(jīng)積累了相當(dāng)?shù)纳降貙?shí)測。Taylor等[4]對其進(jìn)行整理,并與理論計(jì)算進(jìn)行對比,認(rèn)為加速比的理論計(jì)算結(jié)果在近地面處與實(shí)測結(jié)果差距較大。因此山地風(fēng)場風(fēng)速分布仍需要進(jìn)一步的研究。在數(shù)值模擬與試驗(yàn)方面,Bitsuamlak等[5]整理并對比了二維陡坡、緩坡的數(shù)值模擬與風(fēng)洞試驗(yàn),認(rèn)為數(shù)值模擬結(jié)果在迎風(fēng)坡與試驗(yàn)?zāi)軌蜉^好匹配。Cao Shuyang 等[6-7]與Yassin等[8]采用二維CFD研究了山體粗糙度,山體坡度與來流湍流對山地風(fēng)場加速效應(yīng)的影響,其中Yassin等對比了二維CFD模擬與實(shí)驗(yàn)結(jié)果,認(rèn)為其吻合程度較好,用CFD方法進(jìn)行山地風(fēng)場研究具有可行性。Takeshi等[9]利用自主研發(fā)的三維湍流風(fēng)速測試裝置對典型陡坡山體進(jìn)行了實(shí)驗(yàn)研究,其研究在山體迎風(fēng)坡觀測到了相當(dāng)于來流相同高度28.5%水平風(fēng)速大小的豎向風(fēng)速。然而其對豎向風(fēng)速并未進(jìn)行進(jìn)一步詳細(xì)研究。在國內(nèi),李正良等[10-13]對余弦形狀山體分別進(jìn)行了二維、三維數(shù)值模擬和試驗(yàn)研究。研究考察了山體高度、山體坡度與來流地貌等因素對山地風(fēng)場加速比的影響,并對地形修正系數(shù)沿山坡分布情況進(jìn)行了擬合。此外姚旦[14-15]對余弦形單山、雙山也進(jìn)行了細(xì)致的試驗(yàn)與數(shù)值模擬研究,研究提供了較為詳細(xì)的數(shù)值模擬方法,且模擬結(jié)果與試驗(yàn)結(jié)果吻合程度較好。

        上述研究工作加深了對山地地形風(fēng)場分布的認(rèn)識,提供了較為成熟的風(fēng)洞試驗(yàn)與數(shù)值模擬方法。然而,上述研究存在以下兩點(diǎn)不足:其一,目前的大部分研究著重于順風(fēng)向風(fēng)速分布,風(fēng)洞試驗(yàn)中亦只測定順風(fēng)向風(fēng)速,而對于豎向風(fēng)速分布及其隨著山高、坡度等形狀參數(shù)變化的研究處于空白狀態(tài)。實(shí)際工程中豎向風(fēng)速的存在會(huì)影響山地地形中工程結(jié)構(gòu)的使用與安全。例如在山區(qū)輸電線路中,上升氣流會(huì)將導(dǎo)線上托,減小導(dǎo)線張力,從而導(dǎo)線產(chǎn)生更大的風(fēng)偏角而發(fā)生閃絡(luò)事故。故應(yīng)當(dāng)對豎向風(fēng)速分布及其隨山體形狀參數(shù)變化進(jìn)行研究;其二,目前研究的山體大多為無山脈長度的三維軸對稱山體(“饅頭山”)或山脈長度無窮大的二維山體,并未研究山脈長度介于兩者之間的情況。雖然國外研究了具有長短軸的橢圓山體模型[16],但該研究實(shí)質(zhì)仍是無山脈長度的山包。由于實(shí)際山體往往具有一定的山脈長度,而且山脈長度可能對山地風(fēng)場的分布情況產(chǎn)生顯著影響[17-18],故應(yīng)當(dāng)對不同山脈長度下的山地風(fēng)場進(jìn)行研究。

        在此研究背景下,本文參考以往研究,對具有一定山脈長度的陡坡山地模型進(jìn)行了風(fēng)洞試驗(yàn),在該風(fēng)場中測試得到了豎向風(fēng)速分布。而后對照試驗(yàn)結(jié)果利用求解雷諾時(shí)均N-S方程(Reynolds-Averaged Navier-Stokes,RANS)模擬得到不同形狀參數(shù)下山地地形平均風(fēng)場信息,用以研究典型山地地形豎向風(fēng)速的大小,分布及其隨山高、坡度、山脈長度的變化,為實(shí)際山地中的工程建設(shè)提供參考。

        1 風(fēng)洞試驗(yàn)

        1.1 試驗(yàn)概況

        本次試驗(yàn)在浙江大學(xué)ZD-1風(fēng)洞中進(jìn)行。本文研究的山體形狀為目前研究最廣泛,且與實(shí)際山體較為接近的余弦形。其橫截面定義如式(1)所示。其中H為山體高度,D為山體底面圓直徑。

        試驗(yàn)地形山底直徑300 m,山高100 m,具有300 m的山脈長度。模型表面覆蓋綠色纖維,用以增加山體表面粗糙度,以減小雷諾數(shù)效應(yīng)對試驗(yàn)近地面結(jié)果的影響。模型圖片如圖1所示。試驗(yàn)幾何縮尺比為1∶500,來流為相同比例B類地貌來流。來流平均風(fēng)剖面如圖2(a)所示,試驗(yàn)剖面與規(guī)范平均風(fēng)剖面能夠較好吻合,試驗(yàn)原型10 m高度處風(fēng)速為8.10 m/s。試驗(yàn)測點(diǎn)三維風(fēng)速測定采用澳大利亞TFI(Turbulence Flow Instruments)公司生產(chǎn)的眼鏡蛇湍流風(fēng)速測量儀(Series 100 Cobra Probe)。眼鏡蛇三維脈動(dòng)風(fēng)速探頭是4孔壓力探頭,可用于測量三維風(fēng)速和靜態(tài)壓力。試驗(yàn)對來流風(fēng)速三維脈動(dòng)情況進(jìn)行了測試,測得三維湍流度如圖2(b~d)所示。測試所得順風(fēng)向湍流度略大于B類規(guī)范規(guī)定的湍流度,與ESDU規(guī)定的三維湍流度吻合較好。試驗(yàn)測點(diǎn)布置和風(fēng)向角定義如圖3所示,試驗(yàn)工況來流為90°風(fēng)向角。試驗(yàn)測點(diǎn)布置在山脊與迎風(fēng)坡處,測點(diǎn)分為實(shí)心測點(diǎn)與空心測點(diǎn),實(shí)心測點(diǎn)測定了離地10~400 m高度的三維風(fēng)速,空心測點(diǎn)只測定離地10 m高度處的三維風(fēng)速。

        圖1 試驗(yàn)?zāi)P驼掌現(xiàn)ig.1 Photo of experiment model

        (a) 來流平均風(fēng)剖面 (b) 來流順風(fēng)向湍流

        (c) 來流橫風(fēng)向湍流 (d) 來流垂直風(fēng)向湍流

        圖2試驗(yàn)來流條件
        Fig.2Inletconditionofexperiment

        圖3 風(fēng)向角定義和測點(diǎn)布置Fig.3 Definition of wind angle and layout of measurement points

        1.2 試驗(yàn)結(jié)果

        圖4展示了迎風(fēng)坡離地10 m高度處豎向風(fēng)速分布情況,其中正值豎向風(fēng)速為上升氣流,負(fù)值豎向風(fēng)速為下降氣流。由試驗(yàn)結(jié)果可知,上升氣流在迎風(fēng)坡上普遍存在,其中上升氣流風(fēng)速出現(xiàn)臨近山脊的2/3山體高度以上區(qū)域。

        試驗(yàn)各測點(diǎn)中,A1位置處為山脊位置,該處豎向風(fēng)速在2.5 m/s左右,為來流風(fēng)速(8.10 m/s)的31%。A2位置處為迎風(fēng)坡,豎向風(fēng)速達(dá)約5.5 m/s,占來流水平風(fēng)速68%,由此說明了山地地形中豎向風(fēng)速占到來流較大比重,應(yīng)當(dāng)在工程實(shí)踐中予以重視。

        2 CFD數(shù)值模擬與驗(yàn)證

        2.1 建模概況

        本文采用CFD模擬進(jìn)行大量工況的研究。為了計(jì)算結(jié)果更準(zhǔn)確,收斂速度更快,模擬時(shí)統(tǒng)一流場尺寸設(shè)定,流場橫向尺寸為7(L+D),縱向尺寸為11(L+D)+D,其中L為山脈長度,D為山體底面直徑,流場高度為6倍的山體高度。流域的阻塞比在3%以下,如此較好的保證了流場的充分發(fā)展。網(wǎng)格劃分時(shí),山體表面網(wǎng)格尺寸為5 m,底面網(wǎng)格尺寸擴(kuò)展率為1.05,最大網(wǎng)格尺寸為30 m。豎向首層網(wǎng)格尺寸為1 m,豎向擴(kuò)展率為1.1,豎向最大網(wǎng)格尺寸為20 m。網(wǎng)格劃分采用cooper方法。入口采用速度入口,出口采用自由出口,地面與山體采用壁面,側(cè)面與頂面采用對稱面。地面粗糙高度為0.5 m,山體表面粗糙高度為1 m。

        模擬采用入口條件為規(guī)范B類地貌來流,湍流模型選取Realizablek-ε模型,采用非平衡壁面函數(shù)(Non-Equilibrium Wall Function)。風(fēng)剖面與湍動(dòng)能剖面如式(2)~式(5)定義[19]:

        k(z)=0.5[U(z)×I(z)]2(3)

        Lu=0.2(y/30)0.5(5)

        其中,Z0和U0分別表示標(biāo)準(zhǔn)參考高度和標(biāo)準(zhǔn)參考風(fēng)速,本文取10 m和10 m/s,剖面指數(shù)α對B類地貌取為0.15。I(z)表示z高度下的湍流度,按照日本規(guī)范[20]的定義計(jì)算。Cμ取值為0.09,K取值為0.42,Lu為湍流積分尺度,同樣以日本規(guī)范定義。

        圖5為CFD入口處與距入口1000 m位置風(fēng)剖面情況的對比,有對比可知入口平均風(fēng)剖面具有良好的自保持性。在此條件下入口湍流度剖面的自保持性不佳,但由于文中僅針對平均風(fēng)場進(jìn)行研究,應(yīng)優(yōu)先考慮平均風(fēng)剖面的自保持性,且湍流剖面對平均風(fēng)的計(jì)算結(jié)果影響較小,因此可暫時(shí)降低對入口湍流剖面的要求[17]。

        圖5 CFD入口條件及其自保持性Fig.5 CFD inlet condition and its self-maintenance

        2.2 CFD方法的驗(yàn)證

        本文進(jìn)行了風(fēng)洞試驗(yàn)工況的數(shù)值模擬,并對比了兩者加速比與豎向風(fēng)速結(jié)果,加速比按式(6)定義:

        其中,Ux為水平風(fēng)速,U0為相同離地高度來流水平風(fēng)速。

        CFD結(jié)果對比情況如圖6所示,試驗(yàn)與CFD的水平風(fēng)速加速比在山頂與山腳處吻合較好,在山腰流場復(fù)雜區(qū)域,兩者結(jié)果有些差異,但總體趨勢吻合。此外,試驗(yàn)與數(shù)值模擬結(jié)果中豎向風(fēng)速分量在各個(gè)區(qū)域吻合程度普遍較好。

        為了進(jìn)一步驗(yàn)證本文CFD方法,對參考文獻(xiàn)[9]研究的40 m高,底部直徑200 m的三維軸對稱山體進(jìn)行建模計(jì)算。如圖7所示,CFD加速比結(jié)果與Takeshi的研究結(jié)果同樣具有很好的吻合程度。因此通過本文的CFD方法整理總結(jié)豎向風(fēng)速在典型山地中的分布具有可行性。

        (a) 水平風(fēng)速加速比

        (b) 豎向風(fēng)速

        (a) 山腳位置

        (b) 迎風(fēng)坡位置

        (c) 山頂位置

        3 CFD研究典型地形豎向風(fēng)速分布

        3.1 軸對稱山體

        本文對無山脈長度的“饅頭山”進(jìn)行詳細(xì)研究。研究不同山體高度,坡度下豎向風(fēng)速的分布情況。工況設(shè)置如表1所示。

        表1 軸對稱山體模擬工況設(shè)置Table 1 Case of numerical simulation with axisymmetric hill

        3.1.1 山體高度對豎向風(fēng)速分布的影響

        圖8(a)所示為豎向風(fēng)速比Rz沿山體迎風(fēng)坡中軸線的分布。其中x含義與圖3相同,在不具有山脈長度的情況下,-2x/D為0和1時(shí)分別位于山頂與迎風(fēng)山腳。在本文中,豎向風(fēng)速比Rz按式(7)定義:

        其中,Uz、U0分別為豎向風(fēng)速和相同離地高度來流水平風(fēng)速。

        如圖8(a)所示,不同山體高度下,迎風(fēng)坡最大豎向風(fēng)速比均出現(xiàn)在山坡2/3山高以上的位置,且出現(xiàn)位置隨山體高度增加略為上移。最大豎向風(fēng)速比數(shù)值隨山體高度增加而增加,從山高50 m至山高400 m,迎風(fēng)坡最大豎向風(fēng)速比增加為188.12%。-2x/D小于0.8的區(qū)域,即迎風(fēng)坡1/3山高以上區(qū)域,離地10 m高度的豎向風(fēng)速均占來流10%以上。因此在迎風(fēng)坡1/3山高以上區(qū)域應(yīng)適當(dāng)考慮豎向風(fēng)速。

        圖8(b)所示為不同山體高度,山頂處豎向風(fēng)速比隨離地高度的變化情況。如圖所示離地10 m處的豎向風(fēng)速比隨山體高度先增加后減小,山體高度為100 m時(shí)山頂離地10 m高豎向風(fēng)速比存在最大值。從山頂豎向風(fēng)速比剖面來看,山體高度在100 m以下時(shí),最大豎向風(fēng)速比出現(xiàn)在離地10 m高度處,豎向風(fēng)速比隨離地高度增加而減小。當(dāng)山體高度大于200 m時(shí),豎向風(fēng)速比隨離地高度先增加后減小,最大豎向風(fēng)速比出現(xiàn)的離地高度也隨山體高度增加而遞增,在山高400 m時(shí),最大豎向風(fēng)速比出現(xiàn)在離地75 m左右。

        (a) 迎風(fēng)坡豎向風(fēng)速比分布

        (b) 山頂豎向風(fēng)速比剖面

        3.1.2 山體坡度對豎向風(fēng)速分布的影響

        如圖9(a)所示在山體高度一定時(shí),不同山體底部直徑下豎向風(fēng)速比沿迎風(fēng)坡中軸線的分布情況。如圖所示,在山體高度一定時(shí),迎風(fēng)坡最大豎向風(fēng)速比隨著底部直徑增大而減小,底部直徑8H時(shí)最大豎向風(fēng)速比為底部直徑2H時(shí)的55.64%。最大豎向風(fēng)速比出現(xiàn)的位置同樣保持在2/3山高以上的迎風(fēng)坡處,且隨著底部直徑增大而逐漸下移。在山體底部直徑變化時(shí),-2x/D大于0.8的區(qū)域,即迎風(fēng)坡1/3山高以上區(qū)域,離地10 m高度的豎向風(fēng)速比同樣均占來流10%以上。

        如圖9(b)所示,相同山體高度時(shí),山頂最大豎向風(fēng)速比隨著山體底部直徑增加而遞減。值得注意的是在底部直徑小于5H時(shí),山頂?shù)呢Q向風(fēng)速均占來流風(fēng)速5%以下,此時(shí)山頂位置可不考慮豎向風(fēng)速。

        3.2 不同山脈長度下的豎向風(fēng)速分布

        本節(jié)主要研究山脈長度對豎向風(fēng)速分布的影響。表2展示了含山脈山體的模擬工況。該處H與D分別固定為100 m 和300 m 。其中風(fēng)向角定義如圖3所示。

        (a) 迎風(fēng)坡豎向風(fēng)速比分布

        (b) 山頂豎向風(fēng)速比剖面

        Wind angleRidge length L90°0H, 0.5H, 1.0H, 2.0H, 3.0H, 4.0H, 5.0H, 8.0H0°0H, 1.0H, 3.0H, 5.0H, 8.0H

        3.2.1 來流風(fēng)向垂直山脈(90°風(fēng)向角)

        如圖10所示,迎風(fēng)坡坡最大豎向風(fēng)速比與山脊中部豎向風(fēng)速比在離地10 m 處均隨著山脈長度遞增,且遞增趨勢逐漸放緩。這是由于在該風(fēng)向角下,來流正對山脈,山體中部繞流隨著山脈長度增加逐漸接近二維情況,越山效應(yīng)更為顯著。山坡最大豎向風(fēng)速比與山脊中部豎向風(fēng)速比從“饅頭山”至山脈長度達(dá)8H情況,增幅分別達(dá)15%和87%。當(dāng)山脈長度達(dá)到8H時(shí),山脊中部豎向風(fēng)速比可達(dá)0.29。

        圖10 90°風(fēng)向角10 m高度豎向風(fēng)速比Fig.10 Vertical wind velocity ratio at the height of 10 m under 90° wind angle

        如圖11所示為離地10 m高度豎向風(fēng)速比沿山脊的分布情況。其中x含義如圖3所示,L表示山脈長度。x/L為-0.5與0.5時(shí),分別為山脊左、右側(cè)端部。如圖所示豎向風(fēng)速比在山脊中間較大,在左右側(cè)較小。其原因是山脊中部越山效應(yīng)更加明顯。相同山脈長度下,豎向風(fēng)速比沿山脊的變化主要發(fā)生在左右端部1/10L范圍內(nèi),在山脊中間區(qū)域變化相對平緩。隨著山脈長度增加,山脊端部的豎向風(fēng)速比變化不明顯,而山脊中部豎向風(fēng)速比明顯增加。

        圖11 90°風(fēng)向角10 m高度豎向風(fēng)速比沿山脊變化Fig.11 Vertical wind velocity ratio changes along the ridge at the height of 10 m under 90° wind angle

        3.2.2 來流風(fēng)向順山脈(0°風(fēng)向角)

        圖12所示為來流順山脈情況下,迎風(fēng)坡與山脈最大豎向風(fēng)速比隨山脈長度變化情況。如圖所示,迎風(fēng)坡最大豎向風(fēng)速比隨山脈長度增加顯著減少,并在山脈長度達(dá)到5H后保持不變。山脈長度為5H時(shí)迎風(fēng)坡最大豎向風(fēng)速比僅為“饅頭山”情況的58%。山頂最大豎向風(fēng)速比隨著山脈長度的變化幾乎不發(fā)生改變。

        圖12 0°風(fēng)向角10 m高度豎向風(fēng)速比Fig.12 Vertical wind velocity ratio at the height of 10 m under 0° wind angle

        在該風(fēng)向角下,離地10 m高度處的豎向風(fēng)速比沿山脊分布如圖13所示,x/L為-0.5與0.5時(shí),分別代表迎風(fēng)坡山頂與背風(fēng)坡山頂。來流到達(dá)山體,在迎風(fēng)坡山頂為上升氣流,在山脊中間區(qū)域幾乎為0,當(dāng)接近背風(fēng)坡時(shí),轉(zhuǎn)變?yōu)榱讼陆禋饬?。?dāng)山脈長度大于3H,山脊上豎向風(fēng)速比變化主要集中在靠近迎風(fēng)坡和背風(fēng)坡的1/10L區(qū)域內(nèi)。

        圖13 0°風(fēng)向角10 m高度豎向風(fēng)速比沿山脊變化 Fig.13 Vertical wind velocity ratio changes along the ridge at the height of 10 m under 0° wind angle

        4 結(jié) 論

        本文首次嘗試研究了山地地形中的豎向風(fēng)速分布隨山體高度、坡度和山脈長度的變化。研究發(fā)現(xiàn)豎向風(fēng)速在不同山體高度、坡度和山脈長度的山地地形中普遍存在。由試驗(yàn)結(jié)果可知,在陡坡山地地形中豎向風(fēng)速可達(dá)到來流風(fēng)速60%以上,應(yīng)當(dāng)在工程實(shí)踐中予以重視。本文填補(bǔ)了工程中無法參照已有研究考慮豎向風(fēng)速的空白。

        本文研究中考慮了山脈長度對山地豎向風(fēng)速分布影響,發(fā)現(xiàn)山脈長度顯著影響山地地形中的豎向風(fēng)速分布,且豎向風(fēng)速沿山脊分布主要在接近兩側(cè)山頂?shù)?/10山脈長度范圍內(nèi)發(fā)生變化。在對山地地形風(fēng)場的研究中考慮山脈長度的影響是必要的。

        本文研究僅展示了豎向風(fēng)速在山地地形中的分布情況,而對其分布的擬合及如何有效在工程設(shè)計(jì)中考慮豎向風(fēng)速仍有進(jìn)一步研究的需要。

        猜你喜歡
        山脊來流山脈
        它,就在那里
        兩種典型來流條件下風(fēng)力機(jī)尾跡特性的數(shù)值研究
        能源工程(2022年2期)2022-05-23 13:51:48
        Saving the life of a wolf
        不同來流條件對溢洪道過流能力的影響
        山脊新能源
        人,山脈和海洋
        中國三峽(2017年1期)2017-06-09 11:09:41
        “弋”字山上的斜鉤精靈
        基于約束連接方向的最速上升法提取山脊線
        彈發(fā)匹配驗(yàn)證試驗(yàn)系統(tǒng)來流快速啟動(dòng)技術(shù)研究
        冥王星上有一條山脈
        国模91九色精品二三四| 久久精品国产亚洲AV高清特级| 欧美亚洲尤物久久综合精品| a级三级三级三级在线视频| 人妖av手机在线观看| 无码中文字幕免费一区二区三区| 亚洲人成无码www久久久| 国产欧美曰韩一区二区三区| 亚洲国产线茬精品成av| 人妻中出精品久久久一区二| 美女窝人体色www网站| 亚洲中文字幕第一第二页 | 天堂草原电视剧在线观看图片高清| 夜夜春精品视频| 色人阁第四色视频合集网| 男人天堂亚洲一区二区| 色综合久久中文字幕综合网| 欧美中日韩免费观看网站| 中文无码成人免费视频在线观看| 天天摸天天做天天爽天天舒服| 国产精品伦理久久一区| av天堂午夜精品一区| 国产做a爱片久久毛片a片| 中文字幕在线日韩| 天堂av在线播放观看| 91色老久久偷偷精品蜜臀懂色| 女人被爽到高潮视频免费国产 | 日日摸日日碰人妻无码| 欧美性猛交xxxx三人| 97久久超碰国产精品2021| 26uuu欧美日本在线播放| 国产精品女同av在线观看| 亚洲妇熟xxxx妇色黄| 无码中文字幕免费一区二区三区| 国产精品开放小视频| 风流少妇一区二区三区| 国内自拍速发福利免费在线观看| 欧美黑人性暴力猛交喷水| 色诱久久av| 久久婷婷夜色精品国产| 天堂网av一区二区三区在线观看|