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

        ?

        斜拉索干索馳振機(jī)理的數(shù)值模擬與試驗(yàn)研究

        2017-06-19 19:35:13李壽英曾慶宇溫曉光陳政清
        振動(dòng)與沖擊 2017年11期
        關(guān)鍵詞:風(fēng)攻角氣動(dòng)力升力

        李壽英, 曾慶宇, 溫曉光, 陳政清

        (湖南大學(xué) 風(fēng)工程與橋梁工程湖南省重點(diǎn)實(shí)驗(yàn)室, 長沙 410082)

        斜拉索干索馳振機(jī)理的數(shù)值模擬與試驗(yàn)研究

        李壽英, 曾慶宇, 溫曉光, 陳政清

        (湖南大學(xué) 風(fēng)工程與橋梁工程湖南省重點(diǎn)實(shí)驗(yàn)室, 長沙 410082)

        采用數(shù)值模擬和風(fēng)洞試驗(yàn)方法,對(duì)斜拉索干索馳振機(jī)理進(jìn)行研究?;贔LUENT商業(yè)軟件平臺(tái),研究斜索和直索的風(fēng)壓系數(shù)、平均氣動(dòng)力系數(shù)、脈動(dòng)氣動(dòng)力系數(shù)的軸向相關(guān)性等;進(jìn)行粒子圖像測(cè)速(PIV)風(fēng)洞試驗(yàn),對(duì)斜索尾流進(jìn)行可視化顯示,研究斜拉索背后的軸向流特性。研究結(jié)果表明:斜索的平均升力系數(shù)在特定的傾角和風(fēng)攻角下會(huì)發(fā)生突降,其馳振力系數(shù)為較大的負(fù)值(-4.74),可能會(huì)發(fā)生馳振;斜索各截面氣動(dòng)力的軸向相關(guān)系數(shù)比直索要小,甚至?xí)霈F(xiàn)負(fù)的相關(guān)系數(shù);在平均氣動(dòng)力系數(shù)突降的風(fēng)攻角下,尾流中的軸向流不是非常明顯,軸向流對(duì)于干索馳振的影響值得進(jìn)一步定量研究。

        斜拉索; 干索馳振; 數(shù)值模擬; 風(fēng)洞試驗(yàn); 機(jī)理

        斜拉索質(zhì)量輕、頻率低、阻尼小,極易在風(fēng)荷載的作用下發(fā)生大幅振動(dòng)。在斜拉索各種風(fēng)致振動(dòng)中,風(fēng)雨激振是危害最大的一種。風(fēng)雨激振是斜拉索在風(fēng)和雨共同作用下發(fā)生的一種大幅、低頻振動(dòng),自20世紀(jì)80年代以來,研究者對(duì)其機(jī)理進(jìn)行了大量的研究[1-3],研究結(jié)果表明:斜拉索表面形成的上水線對(duì)風(fēng)雨激振的起振關(guān)系密切。與此同時(shí),人們也發(fā)現(xiàn),在有風(fēng)無雨的情況下,斜拉索也會(huì)發(fā)生大幅振動(dòng)。對(duì)于這種類型的拉索振動(dòng),主要有兩種解釋:第一種是高風(fēng)速下的渦激共振,主要由日本學(xué)者M(jìn)atsumoto等[4]提出;第二種是干索馳振(Dry galloping),與平均氣動(dòng)力系數(shù)的突降有關(guān)。Chen等[5]進(jìn)行了三維斜索節(jié)段模型測(cè)振試驗(yàn),在試驗(yàn)室中重現(xiàn)了干索弛振現(xiàn)象;Chen等[6]則在對(duì)斜索表面進(jìn)行剛性模型測(cè)壓的基礎(chǔ)上,通過積分得到了斜索的平均氣動(dòng)力系數(shù),研究結(jié)果表明:在特定的傾角和風(fēng)攻角下,平均升力系數(shù)確實(shí)會(huì)出現(xiàn)突降;Matsumoto等[7]認(rèn)為軸向流對(duì)干索馳振有重要影響,且與雷諾數(shù)從亞臨界區(qū)向臨界區(qū)過渡有關(guān);劉慶寬等[8]對(duì)二維拉索進(jìn)行測(cè)力試驗(yàn)發(fā)現(xiàn),當(dāng)雷諾數(shù)進(jìn)入臨界區(qū)域時(shí),會(huì)產(chǎn)生較大的平均升力。

        本文采用FLUENT軟件平臺(tái)的LES模塊,對(duì)直索和斜索進(jìn)行數(shù)值模擬,并在風(fēng)洞實(shí)驗(yàn)室中進(jìn)行了PIV風(fēng)洞試驗(yàn)。研究了直索和斜索表面風(fēng)壓系數(shù)、平均氣動(dòng)力系數(shù)、氣動(dòng)力系數(shù)軸向相關(guān)性、尾流特性等,在此基礎(chǔ)上,對(duì)干索馳振機(jī)理進(jìn)行了一定的解釋。

        1 數(shù)值模型

        斜拉索空間姿態(tài)采用傾角α和風(fēng)攻角β來定義,如圖1所示。拉索直徑D取為150 mm。數(shù)值模擬對(duì)象包括兩種:直索和斜索。對(duì)于直索模型,計(jì)算流域的長、寬、高分別為34D、6D和20D,如圖2(a)所示。其中,拉索中心線與計(jì)算流域上、下表面的距離均為10D,與入口的距離為10D,與出口的距離為24D。對(duì)于斜索模型,拉索上、下端點(diǎn)與計(jì)算流域上、下表面的距離分別為10D,拉索前、后端點(diǎn)與計(jì)算流域的入口、出口的距離分別為10D、24D,如圖2(b)所示。入口采用速度邊界條件;出口采用壓力邊界條件;上、下表面采用對(duì)稱邊界條件;左、右表面采用自由滑移壁面條件;拉索表面采用無滑移壁面條件,如圖3所示。

        圖1 拉索的空間姿態(tài)

        (a) 直索模型

        (b) 斜索模型

        (a) 直索模型

        (b) 斜索模型

        在拉索表面設(shè)置7個(gè)截面的風(fēng)壓監(jiān)測(cè)點(diǎn),每個(gè)截面與拉索軸線方向垂直,相鄰截面的間距為0.1 m,其中截面4為拉索的中間截面。每個(gè)截面上共均勻布置風(fēng)壓監(jiān)測(cè)點(diǎn)36個(gè),相鄰監(jiān)測(cè)點(diǎn)的圓心角為10°,如圖4所示。

        網(wǎng)格劃分對(duì)數(shù)值模擬結(jié)果的影響很大。采用結(jié)構(gòu)化網(wǎng)格進(jìn)行網(wǎng)格劃分,為保證計(jì)算精度,第一層網(wǎng)格的無量綱高度y+保持在0.30以下,直索和斜索模型的網(wǎng)格總數(shù)分別約為44萬和60萬,并進(jìn)行了網(wǎng)格收斂性驗(yàn)證。拉索截面附近的網(wǎng)格,如圖5所示。從圖5中可以看出,在與來流平行的豎截面內(nèi),直索的截面是圓形,而斜索的截面是橢圓形。

        直索模擬工況只有1個(gè);斜索的傾角α為55°,風(fēng)攻角β在0°~60°之間,基本間隔5°,在25°~35°之間加密至1°,斜索工況總數(shù)為21個(gè)。來流采用均勻流,風(fēng)速為30 m/s和50 m/s,對(duì)應(yīng)的雷諾數(shù)分別為3.0×105和5.0×105。

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

        (a) 直索

        (b) 斜索

        (c) 截面測(cè)點(diǎn)

        (a) 直索

        (b) 斜索

        粒子圖像測(cè)速法(Particle Image Velocimetry, PIV)是一種瞬態(tài)、多點(diǎn)、無接觸式的流體力學(xué)測(cè)速方法,能在同一瞬態(tài)記錄下大量空間點(diǎn)上的速度分布信息,并可提供豐富的流場(chǎng)空間結(jié)構(gòu)以及流動(dòng)特性。試驗(yàn)在湖南大學(xué)風(fēng)工程試驗(yàn)研究中心HD-2大氣邊界層風(fēng)洞實(shí)驗(yàn)室的高度試驗(yàn)段進(jìn)行。該試驗(yàn)段長17 m、寬3 m、高2.5 m,最高風(fēng)速可達(dá)58 m/s。拉索模型直徑150 mm,與數(shù)值模擬相同。拉索長度為1 500 mm,其固定裝置及風(fēng)向角定義如圖6所示,試驗(yàn)中主要采集拉索軸線的豎直平面內(nèi)的風(fēng)速數(shù)據(jù)。試驗(yàn)在均勻流場(chǎng)中進(jìn)行,風(fēng)速為30 m/s。拉索傾角固定為55°,風(fēng)攻角與數(shù)值模擬工況相同,風(fēng)攻角通過轉(zhuǎn)動(dòng)轉(zhuǎn)盤來實(shí)現(xiàn)。

        圖6 拉索模型裝置示意圖

        3 結(jié)果及討論

        模型監(jiān)測(cè)截面氣動(dòng)力的軸向相關(guān)系數(shù)γxy定義如下:

        (1)

        3.1 平均風(fēng)壓系數(shù)

        數(shù)值模擬結(jié)果表明:7個(gè)截面的平均和脈動(dòng)風(fēng)壓系數(shù)結(jié)果非常一致,后文僅分析中間截面結(jié)果。圖7給出了直索模型中間截面的平均風(fēng)壓系數(shù)與文獻(xiàn)結(jié)果的比較,包括Achenbach[9]的風(fēng)洞試驗(yàn)結(jié)果(Re=5.0×105)、Catalano等[10]的數(shù)值模擬結(jié)果(Re=1.0×106)。從圖7中可以看出,在氣流分離區(qū)附近,平均風(fēng)壓系數(shù)有一定的差異,但總體上來說,本文的平均風(fēng)壓系數(shù)結(jié)果與文獻(xiàn)結(jié)果吻合較好。圖8給出了風(fēng)攻角β=0°、10°、20°、30°時(shí)斜索模型中間截面的平均風(fēng)壓系數(shù)周向分布規(guī)律。從圖8中可以看出,與直索模型相比,駐點(diǎn)處的平均風(fēng)壓系數(shù)約為0.7左右,小于直索模型的1.0,這主要是因?yàn)閬砹骱屠鞅砻娌淮怪?,一部分風(fēng)速能量分解到與拉索表面相切的方向上;另外,分離點(diǎn)和背風(fēng)面的風(fēng)壓吸力也小于直索模型的值。斜索平均風(fēng)壓系數(shù)的這些特點(diǎn)與文獻(xiàn)[11]的結(jié)果吻合較好。

        圖7 平均風(fēng)壓系數(shù)結(jié)果的驗(yàn)證

        圖8 斜索的平均風(fēng)壓系數(shù)(Re=3.0×105)

        3.2 氣動(dòng)力系數(shù)

        雷諾數(shù)為3×105時(shí),直索模型的平均阻力系數(shù)為0.88,脈動(dòng)升力系數(shù)的根方差值為0.40,這與經(jīng)典結(jié)果比較吻合[12]。與直索模型不同的是,斜索模型的平均氣動(dòng)力力系數(shù)隨著風(fēng)攻角的變化而變化。圖9給出了斜索模型的平均阻力系數(shù)和平均升力系數(shù)隨風(fēng)攻角的變化曲線,包括雷諾數(shù)為3×105和5×105。從圖9中可以看出,平均升力系數(shù)在25°~32°之間出現(xiàn)了較大的突降,這將使得平均升力系數(shù)對(duì)風(fēng)攻角的導(dǎo)數(shù)為負(fù)值。根據(jù)經(jīng)典的Denhartog馳振理論[13],若馳振力系數(shù)小于零,即:

        (2)

        則表示氣流和結(jié)構(gòu)的相互耦合作用會(huì)引起負(fù)的氣動(dòng)阻尼,從而使得結(jié)構(gòu)發(fā)生劇烈振動(dòng)。圖10給出了根據(jù)圖9的平均氣動(dòng)力系數(shù)得到的馳振力系數(shù)。從圖10中可以看出,兩種雷諾數(shù)情況下,馳振力系數(shù)均在β=25°~32°范圍內(nèi)出現(xiàn)較大的負(fù)值,其中,雷諾數(shù)為3×105時(shí),最大負(fù)馳振力系數(shù)達(dá)到-3.85;雷諾數(shù)為5×105時(shí),最大負(fù)馳振力系數(shù)達(dá)到-4.74。因此,從平均氣動(dòng)力系數(shù)的特點(diǎn)來看,斜索確實(shí)存在發(fā)生馳振的可能性。需要說明的是,圖9和圖10的結(jié)果是拉索傾角α=55°情況下的值,拉索傾角變化,平均氣動(dòng)升力突降的風(fēng)攻角范圍及其對(duì)應(yīng)的馳振力系數(shù)大小也會(huì)發(fā)生變化。從圖5(b)中可以看出,在斜風(fēng)作用下,拉索的迎風(fēng)截面為橢圓形,且隨著風(fēng)攻角的變化迎風(fēng)橢圓截面也會(huì)變化,這應(yīng)該是斜索平均升力系數(shù)發(fā)生突降的主要原因。

        3.3 氣動(dòng)力系數(shù)的軸向相關(guān)性

        圖11和圖12分別給出了雷諾數(shù)為3×105時(shí)、直索和斜索各個(gè)截面與中間截面氣動(dòng)力系數(shù)之間的相關(guān)系數(shù),其中拉索截面的氣動(dòng)力系數(shù)通過監(jiān)測(cè)點(diǎn)風(fēng)壓積分得到。從圖11中可以看出,直索各截面氣動(dòng)力系數(shù)之間的相關(guān)性較好,且均呈現(xiàn)為正相關(guān),截面1、2、3與截面4之間阻力系數(shù)的相關(guān)系數(shù)分別為0.27、0.46、0.47,對(duì)應(yīng)的升力系數(shù)的相關(guān)系數(shù)分別為0.55、0.74、0.86,升力系數(shù)的相關(guān)性高于阻力系數(shù)的相關(guān)性。從圖12中可以看出,斜索各截面氣動(dòng)力系數(shù)之間的相關(guān)性要小于直索的值,其中阻力系數(shù)的相關(guān)系數(shù)主要為正相關(guān),但升力系數(shù)的相關(guān)系數(shù)出現(xiàn)了較大的負(fù)值,特別是在風(fēng)攻角β=20°~35°范圍內(nèi),截面3和截面4的升力系數(shù)的負(fù)相關(guān)系數(shù)可達(dá)到-0.60左右,這說明斜索表面的氣動(dòng)力時(shí)空分布比直索更為紊亂。

        (a) Re=3.0×105

        (b) Re=5.0×105

        3.4 尾流特征

        圖10 斜索的馳振力系數(shù)

        圖11 直索各截面氣動(dòng)力系數(shù)的相關(guān)系數(shù)(Re=3.0×105)

        (a) 阻力系數(shù)

        (b) 升力系數(shù)

        圖13給出了通過PIV成像風(fēng)洞試驗(yàn)得到的斜索尾流的流線圖及速度云圖,其中,風(fēng)攻角β=0°、25°、30°、35°,雷諾數(shù)Re=3×105。圖13中同時(shí)還給出了通過CFD數(shù)值模擬得到的流線圖。從圖14中可以看出,在風(fēng)攻角β=0°、25°時(shí),斜索尾流的軸向流分量非常明顯,但隨著風(fēng)攻角的增大(當(dāng)β=30°、35°時(shí)),軸向流分量有減小的趨勢(shì)。日本學(xué)者M(jìn)atsumoto認(rèn)為軸向流是引起干索馳振的一個(gè)重要原因[7],從本文數(shù)值模擬結(jié)果來看,負(fù)馳振力系數(shù)的風(fēng)攻角范圍內(nèi)(β=25°~32°),軸向流并不明顯,軸向流對(duì)于干索馳振的影響需進(jìn)一步的定量研究。值得注意的是,從數(shù)值模擬得到的流線圖可以發(fā)現(xiàn),當(dāng)風(fēng)攻角β=20°~35°時(shí),存在兩個(gè)方向的軸向流,并交織在一起。

        4 結(jié) 論

        對(duì)直索和斜索模型進(jìn)行了CFD數(shù)值模擬,并在風(fēng)洞中進(jìn)行了斜索的PIV成像風(fēng)洞試驗(yàn),得到如下主要

        (a) β=0°

        (c) β=30°

        (d) β=35°

        圖13 斜索尾流的流線圖(Re=3.0×105)

        Fig.13 Streamline in the wake of the 3-d cable (Re=3.0×105)

        結(jié)論:

        (1) 斜索傾角為55°、風(fēng)攻角在25°~32°范圍內(nèi)時(shí),斜索的馳振力系數(shù)出現(xiàn)負(fù)值,最大負(fù)馳振力系數(shù)可達(dá)到-4.74,這表明斜索有發(fā)生干索馳振的可能性,主要原因可能為斜索迎風(fēng)截面變?yōu)闄E圓形。

        (2) 斜索軸向氣動(dòng)力的相關(guān)系數(shù)遠(yuǎn)小于直索的軸向氣動(dòng)力相關(guān)系數(shù),且升力的軸向相關(guān)系數(shù)大于阻力的軸向相關(guān)系數(shù)。當(dāng)風(fēng)攻角為20°~35°范圍內(nèi)時(shí),斜索升力的軸向相關(guān)系數(shù)為較大的負(fù)值。

        (3) 從PIV成像風(fēng)洞試驗(yàn)結(jié)果來看,當(dāng)斜索傾角為55°時(shí),風(fēng)攻角為0°和25°時(shí)的軸向流非常明顯,但當(dāng)風(fēng)攻角為30°和35°時(shí),尾流的軸向流分量減弱,軸向流的作用值得進(jìn)一步定量研究。

        [1] HIKAMI Y, SHIRAISHI N. Rain-wind induced vibrations of cables in cable stayed bridges[J]. Journal of Wind Engineering and Industrial Aerodynamics, 1988, 29: 409-418.

        [2] GU M, DU X Q. Experimental investigation of rain-wind-induced vibration of cables in cable-stayed bridges and its mitigation[J]. Journal of Wind Engineering and Industrial Aerodynamics, 2005, 93: 79-95.

        [3] LI S Y, CHEN Z Q, WU T, et al. Rain-wind-induced in-plane and out-of-plane vibrations of stay cables[J]. Journal of Engineering Mechanics, 2013, 39(12): 1688-1698.

        [4] MATSUMOTO M, SHIRATO H, YAGI T, et al. Field observation system of cable aerodynamics in natural wind[J]. In: Proceedings of the Fourth International Symposium on Cable Dynamics, Montreal, Canada, 2001: 219-225.

        [5] CHENG S, LAROSE G L, SAVAGE M G, et al. Experimental study on the wind-induced vibration of a dry

        inclined cable—Part I: Phenomena[J]. Journal of wind engineering and industrial aerodynamics, 2008, 96: 2231-2253.

        [6] CHENG S, IRWIN P A, TANAKA H. Experimental study on the wind-induced vibration of a dry inclined cable--Part II: Proposed mechanisms[J]. Journal of Wind Engineering and Industrial Aerodynamics, 2008, 96: 2254-2272.

        [7] MATSUMOTO M, YAGI T, HATSUDA H, et al. Dry galloping characteristics and its mechanism of inclined/yawed cables[J]. Journal of Wind Engineering and Industrial Aerodynamics, 2010, 98(6): 317-327.

        [8] 劉慶寬,張峰,馬文勇,等.斜拉索雷諾數(shù)效應(yīng)與風(fēng)致振動(dòng)的試驗(yàn)研究,振動(dòng)與沖擊,2011, 30(12): 114-119.

        LIU Qingkuan, ZHANG Feng, MA Wenyong, 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] ACHENBACH E. Distribution of local pressure and skin friction around a circular cylinder in cross flow up toRe= 5×106[J]. Fluid Mech, 1968, 34(4):625-39.

        [10] CATALANO P, WANG M, IACCARINO G, et al. Numerical simulation of the flow around a circular cylinder at high Reynolds numbers[J]. Int. J. Heat Fluid Flow, 2003, 24: 463-469.

        [11] 李壽英,顧明.斜、直圓柱繞流的CFD模擬[J].空氣動(dòng)力學(xué)學(xué)報(bào),2005, 23(2): 222-227.

        LI Shouying, GU Ming. Numerical simulation for flow around perpendicular and oblique circular cylinders[J]. Acta Aerodynamica Sinica, 2005, 23(2): 222-227.

        [12] SCHLICHTING H, GERSTEN K. Boundary-laryer theory[M]. New York: McGraw Hill, 1979.

        [13] DENHARTOG J P. Mechanical Vibrations[M]. New York: Dover Publications, INC., 1985.

        Numerical simulations and tests for dry galloping mechanism of stay cables

        LI Shouying, ZENG Qingyu, WEN Xiaoguang, CHEN Zhengqing

        (Key Laboratory for Wind and Bridge Engineering of Hunan Province, Hunan University, Changsha 410082, China)

        The dry galloping mechanism of stay cables was studied by means of CFD simulations and wind tunnel tests. Firstly, CFD simulations with 2-D and 3-D cable models were conducted by utilizing the LES method based on FLUENT software to obtain wind pressure coefficients, mean drag and lift coefficients, correlation coefficients of fluctuating aerodynamic forces along the cable axis. Secondly, wind tunnel tests using the particle image velocimetry (PIV) technique were performed to specially investigate characteristics of stay cable wake flow. The results showed that a sudden decrease in mean lift coefficient of a 3-D stay cable is observed, and a maximum minus galloping force coefficient -4.74 is found, so large amplitude galloping vibration of the stay cable is possible to take place; correlation coefficients between aerodynamic forces of a 3-D cable’s cross sections are smaller than those of a 2-D cable, and correlations of drag coefficients are smaller than those of lift coefficients; the axial flow in a 3-D cable’s wake flow is not very obvious; the effects of axial flow on cable galloping vibration need further quantitatively studying.

        stay cable; dry galloping; CFD simulation; wind tunnel tests; mechanism

        國家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃(2015CB057701; 2015CB057702); 國家自然科學(xué)基金(51578234)

        2015-09-14 修改稿收到日期:2016-04-16

        李壽英 男,教授,博士生導(dǎo)師,1977年5月出生

        U448.27

        A

        10.13465/j.cnki.jvs.2017.11.015

        猜你喜歡
        風(fēng)攻角氣動(dòng)力升力
        Scruton數(shù)對(duì)小寬高比H型斷面典型攻角風(fēng)致振動(dòng)的影響
        高速列車車頂–升力翼組合體氣動(dòng)特性
        風(fēng)攻角對(duì)某大跨斜拉橋氣動(dòng)系數(shù)影響研究
        大科技(2022年20期)2022-05-25 01:53:54
        飛行載荷外部氣動(dòng)力的二次規(guī)劃等效映射方法
        無人機(jī)升力測(cè)試裝置設(shè)計(jì)及誤差因素分析
        不同風(fēng)攻角下薄平板斷面顫振機(jī)理研究
        基于自適應(yīng)偽譜法的升力式飛行器火星進(jìn)入段快速軌跡優(yōu)化
        側(cè)風(fēng)對(duì)拍動(dòng)翅氣動(dòng)力的影響
        基于CFD的流線型橋梁斷面阻力系數(shù)測(cè)壓結(jié)果修正研究
        升力式再入飛行器體襟翼姿態(tài)控制方法
        国产日韩欧美网站| 久久午夜羞羞影院免费观看| 又大又粗又爽的少妇免费视频| 中国年轻丰满女人毛茸茸| 无码8090精品久久一区| 宅男视频一区二区三区在线观看| 极品老师腿张开粉嫩小泬| 亚洲人成色777777老人头| 91精品91| 亚洲一区在线二区三区 | 冲田杏梨av天堂一区二区三区| 亚洲最大在线视频一区二区| 少妇人妻精品一区二区三区| 中文幕无线码中文字蜜桃| 日本岛国大片不卡人妻| 成人国产精品三上悠亚久久| 中国人妻与老外黑人| av网站免费线看| 看全色黄大色大片免费久久久| 狼人伊人影院在线观看国产| 人妻 色综合网站| 久久久久成人亚洲综合精品| 国产精品久久中文字幕亚洲| 国产精品女主播福利在线| 国产精品久久久久久久久岛| a级福利毛片| 亚洲免费一区二区av| 人妻精品久久久久中文字幕| 一道久在线无码加勒比| 精品人妻av一区二区三区不卡| 一本久道竹内纱里奈中文字幕| 野花社区视频在线观看| 免费av在线国模| 福利视频在线一区二区三区| 麻豆国产精品va在线观看不卡| 国产999精品久久久久久| 国产成年无码久久久久下载| 亚洲成人av在线蜜桃| 国产中文欧美日韩在线| 无码不卡一区二区三区在线观看| 丝袜美腿亚洲综合一区|