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

        ?

        顫振臨界風速計算值與試驗值的一致性

        2018-06-01 02:59:41李志國廖海黎
        西南交通大學學報 2018年3期
        關(guān)鍵詞:氣動力風洞試驗攻角

        伍 波, 王 騎, 李志國, 廖海黎

        (1. 西南交通大學土木工程學院, 四川 成都 610031; 2. 西南交通大學風工程四川省重點實驗室, 四川 成都 610031)

        隨著橋梁建造技術(shù)的進步和橋梁跨度的顯著增加,橋梁的顫振顯得越來越重要,而如何準確計算橋梁的顫振臨界風速又是大跨度橋梁乃至超大跨度橋梁設(shè)計時最需要解決的問題.至今為止,國內(nèi)外已經(jīng)對橋梁顫振計算方法做了相當多的研究,也取得了豐碩成果.現(xiàn)在比較常用的顫振分析方法,主要包括:復(fù)特征值法(complex eigenvalue analysis, CEV)[1]、Masumoto提出的分步分析方法[2]、時域分析方法[3]、雙模態(tài)耦合閉合解法[4].這些方法涉及的自激氣動力,皆源于Scanlan等[5]在1971年提出的以顫振導數(shù)表示的氣動力模型,由此派生出來的顫振導數(shù)識別方法[6-8]也在理論上具有相當高的精度.

        盡管顫振導數(shù)的識別方法和顫振分析方法的準確性已在薄平板算例及機翼算例中得到了很好的驗證,理論本身的正確性是毋庸置疑的.但對于橋梁顫振計算來說,即使在小振幅情況下,理論計算結(jié)果與風洞試驗測試值之間總會存在或多或少的差異,尤其在斷面較鈍或攻角較大的情況下,這種差異更加明顯.Argentini等[9]針對某大跨度懸索橋進行全橋氣彈模型試驗及多模態(tài)復(fù)特征值顫振分析,結(jié)果顯示即使在考慮了全模態(tài)(full-set)的情況下,顫振臨界風速計算結(jié)果與試驗結(jié)果仍有一定差別.丁泉順[10]對不同攻角下江陰大橋的顫振分析結(jié)果顯示,小攻角下,顫振風速理論計算結(jié)果與試驗結(jié)果保持一致,而在3°攻角下,計算值與試驗值的差異明顯增大.Ge等[11]針對H型斷面及箱形斷面進行二維顫振分析,結(jié)果顯示H型斷面的計算誤差明顯大于箱形斷面.這樣的差異性曾被認為是氣動力的高次諧波效應(yīng)引起.但王騎等[12-13]在扁平箱梁和矩形斷面的強迫振動風洞試驗研究中發(fā)現(xiàn),在小振幅條件下,兩種斷面均沒有觀察到顯著的高次諧波分量,即便是在中等振幅條件下(扭轉(zhuǎn)角達到10°),其自激氣動力中的高次諧波分量比重仍然不大,由此可知,造成此差異性的原因還另有他解.

        為了詳細研究橋梁顫振臨界風速計算值和試驗值之間的差異性,并發(fā)現(xiàn)其中可能的原因,本文以扁平箱梁為研究對象.首先利用強迫振動風洞試驗獲得斷面顫振導數(shù),并基于不同的動力參數(shù)組合計算獲得了不同工況下顫振風速.其次采用常規(guī)的彈簧懸掛節(jié)段模型風洞試驗技術(shù),選取與計算相同的動力參數(shù)組合,測試并獲得同一模型在不同風攻角下的不同顫振特性,包括臨界風速、振幅比和相位角,并與計算結(jié)果進行對比.結(jié)果表明,0°風攻角下計算值和試驗值一致性較好,而3°和5°風攻角下計算值和試驗值有較大的差異性.最后基于計算的和試驗的顫振因子對比以及顫振因子與顫振相位角成正變的特性,推斷出在攻角較大時,耦合相位角可能會對顫振導數(shù)產(chǎn)生不可忽略的影響,從而導致顫振風速計算值和試驗值產(chǎn)生差異.

        1 風洞試驗參數(shù)與計算方法

        1.1 強迫振動試驗及顫振導數(shù)識別

        本文以某扁平箱梁斷面為對象,制作了縮尺比為 1∶80 的節(jié)段模型,模型長度為1.1 m,寬度為0.4 m,梁高為0.041 m,高寬比為9.7.模型斷面如圖1所示,強迫振動風洞試驗如圖2 所示.

        圖1 試驗?zāi)P蛿嗝鍲ig.1 Cross-Section of model

        圖2 強迫振動風洞試驗Fig.2 Forced vibration testing for section model

        利用強迫振動試驗裝置,通過記錄模型在不同風速及不同風攻角下的豎向及扭轉(zhuǎn)單自由度運動的氣動自激力和位移,識別模型在不同工況下的顫振導數(shù).試驗中采取的振動頻率為2.5 Hz,豎向試驗振幅為10 mm,扭轉(zhuǎn)試驗振幅為2°.豎向和扭轉(zhuǎn)運動引起橋梁斷面自激氣動力分別為

        (1)

        (2)

        試驗時的折算風速V=U/(fb),其中,f為強迫振動頻率,Hz,折算風速范圍從4~20.

        1.2 顫振計算方法

        本文選擇雙模態(tài)耦合閉合解法[4]進行顫振臨界風速計算.該方法不僅能夠簡便和精確地計算橋梁的顫振臨界風速,也可求解出模態(tài)頻率、阻尼及相位等參數(shù)隨著折減風速的變化.同時,由于該方法將顫振風速表示為動力部分和氣動部分的乘積,使得研究中可以單獨量化這兩部分的影響,有利于分析導致顫振風速計算值和試驗值出現(xiàn)差異的原因.該計算方法具體表達式為

        (3)

        (4)

        (5)

        (6)

        式中:Uc為是顫振臨界風速;ωs1和ωs2分別為豎向運動和扭轉(zhuǎn)運動的圓頻率;m為主梁單位長度等效質(zhì)量;r為等效質(zhì)量的慣性半徑;γ為表征顫振導數(shù)影響的系數(shù),本文稱為顫振因子,反應(yīng)了氣動力變化對顫振的影響;ξs2為扭轉(zhuǎn)運動的結(jié)構(gòu)阻尼比,%;υ=ρb4/I,I為單位長度等效質(zhì)量慣性矩;D為表征彎扭振型相似度的因子,對于節(jié)段模型試驗D=1.

        根據(jù)耦合顫振閉合解簡化計算方法[4],最終的顫振風速則是式(3)的風速值曲線和式(7)的風速值曲線的交點.

        (7)

        此外,還可以通過風洞試驗獲得的顫振風速和試驗系統(tǒng)動力參數(shù)反算出顫振因子,本文稱為試驗顫振因子,如式(8).

        (8)

        式中:Ucs為顫振臨界風速試驗值.

        利用上述顫振計算方法及測量得到的顫振導數(shù),即可計算不同動力參數(shù)組合下的顫振臨界風速,主要的計算工況如表1所示.利用強迫振動試驗裝置得到的決定顫振臨界風速的4個關(guān)鍵顫振導數(shù)如圖3所示.

        表1 顫振性能計算及測試工況Tab.1 Calculating and testing cases for flutter

        注:fh、ft分別為系統(tǒng)豎向頻率及扭轉(zhuǎn)頻率,Hz;M、Im分別為系統(tǒng)總質(zhì)量及質(zhì)量慣性矩,kg、kg·m2;

        ξs1為系統(tǒng)豎向阻尼比,%;η為系統(tǒng)的彎扭頻率比.

        1.3 自由振動試驗?zāi)P图肮r

        制作了縮尺比為 1∶50 的節(jié)段模型進行彈簧懸掛節(jié)段模型風洞試驗.試驗在XNJD-1第二試驗段的均勻流場中進行(試驗段截面尺寸為2.4 m×2.0 m,最大試驗風速為45 m/s),采用可模擬二自由度耦合運動(豎向運動及扭轉(zhuǎn)運動)的彈簧懸掛系統(tǒng)進行不同風攻角下的節(jié)段模型顫振性能測試.

        (a) H3(b) A?1(c) A?2(d) A?3圖3 試驗斷面不同攻角下4個關(guān)鍵顫振導數(shù)的測試結(jié)果Fig.3 Flutter derivatives of the girder under different angle of attacks

        為了較詳細地反映顫振風速計算值和試驗之間的差異,試驗中通過改變配重的質(zhì)量及質(zhì)量慣性矩來改變試驗系統(tǒng)的Scruton數(shù)[14]并滿足與計算工況相同的動力參數(shù)組合.

        2 計算結(jié)果和試驗結(jié)果的對比

        從表2中可以看出:對于0°攻角,顫振臨界風速試驗值與計算值之間具有良好的吻合度;對于工況1~6,不同動力參數(shù)下顫振折算風速最大計算誤差為4.9%,其試驗值為8.35,計算值為7.96,對應(yīng)的為工況3,最小計算誤差為0.12%,計算值為8.63,試驗值為8.62,對應(yīng)的為工況1;0°攻角下的計算結(jié)果表明,利用顫振導數(shù)計算的顫振臨界風速和顫振頻率,無論在何種動力參數(shù)條件下,均能夠和風洞試驗值高度吻合.

        表2 顫振臨界風速試驗值與計算值的對比Tab.2 Comparison of calculated results and tested results of critical flutter speed

        隨著攻角增大,計算結(jié)果與實驗結(jié)果的吻合度逐漸降低.從圖3可以看出:在3°攻角,工況5下顫振折算風速計算值和試驗值分別為6.80和7.37,最大差異達到7.7%;工況1下對應(yīng)的計算和試驗差異為6.2%;在5°攻角下,兩者間的差異進一步增大;工況6下的計算折減風速為6.05,試驗折減風速為6.75,誤差達到了10.4%,對于5°攻角下的其他測試工況,計算值和試驗值的差異性也很明顯,如工況4的差異為9.3%,工況5的差異為9.1%,工況2的差異為6.66%,工況1的差異為5.4%.

        由此可以得出,在有攻角的情況下,除了顫振臨界風速本身較攻角0°下的值有顯著降低外,臨界風速的計算結(jié)果與試驗結(jié)果間的差異也顯著增加,尤其在5°攻角下,該差異達到最大.需要說明的是,在有風攻角的條件下,這樣的差異在其他橋梁斷面的風洞試驗數(shù)據(jù)中可以找到[9-11].

        3 顫振風速差異的原因分析

        3.1 差異離散性的量化

        由于在不同風攻角下相同試驗工況的動力參數(shù)是相同的,因此可以推知顫振臨界風速的差異性不是由動力參數(shù)的改變引起的,而是由氣動力的變化引起的.基于此,采用式(3),可以將動力參數(shù)的影響與氣動參數(shù)的影響分離開.根據(jù)現(xiàn)有顫振理論,顫振導數(shù)只與斷面的氣動外形有關(guān),因此只要在相同風攻角下,無論結(jié)構(gòu)的動力參數(shù)如何改變,顫振導數(shù)均保持不變,那么顫振因子的值會保持穩(wěn)定,同時試驗顫振因子和理論顫振因子的值應(yīng)保持一致.基于此推論,通過分析顫振因子的值是否保持一致以及其理論值和試驗之間的差異,就可推測出引起顫振臨界風速差異的原因.

        采用實測的顫振導數(shù),結(jié)合式(4)可以獲得理論上的γ,由此可量化氣動力對顫振的影響;利用每個工況所對應(yīng)的系統(tǒng)動力參數(shù),結(jié)合節(jié)段模型試驗所獲得的顫振臨界風速,采用式(8)可反算出各個工況對應(yīng)的試驗顫振因子γs.為了分析顫振因子數(shù)值的穩(wěn)定性,再定義參數(shù)C1和C2用以表征顫振因子在不同工況下的相對變化,其表達式如下所示:

        (9)

        計算和試驗條件下的顫振因子和對應(yīng)的C1和C2值如表3所示.從表3可以看出:對于顫振導數(shù)計算獲得的理論顫振因子,無論在哪個攻角下,都比較集中和穩(wěn)定,彼此間差異值C1基本都在2%以下;對于試驗風速反算的顫振因子γs以及差異值C2,除了在攻角0°下顫振導數(shù)的值穩(wěn)定外(C2<5%),在攻角3°和5°下,C2都超過了5%.如在0°攻角時,顫振因子分布在[0.407,0.426]之間,C2最大為4.5%;在攻角3°條件下,顫振因子分布在[0.372,0.409]之間,C2增大到8.8%,折算風速區(qū)間也縮短為[7,8];攻角5°時,顫振因子分布在[0.336,0.381],C2上升到了11.8%,對應(yīng)的折算風速區(qū)間縮短為[6,7].圖4則直觀給出了顫振因子的計算值和測試值的相對變化.圖4中所示理論曲線為利用式(4)和顫振導數(shù)計算得到的不同折算風速下的顫振因子曲線,空心點為相應(yīng)工況的顫振臨界狀態(tài)對應(yīng)的顫振因子計算值γs;圖中實心點為相應(yīng)工況在試驗顫振臨界風速處利用式(8)反算得到的顫振因子γs。

        在攻角0°下,試驗顫振因子與計算的顫振因子保持一致,各個顫振因子也均被包圍在理論顫振因子隨折算風速變化的曲線中.但在攻角3°和5°下,試驗顫振因子明顯偏離了理論曲線,在顫振折算風速區(qū)間內(nèi)的離散性顯著高于理論曲線的離散性.

        表3 顫振因子計算值與試驗反算值對比Tab.3 Comparison of calculated results and tested results

        由于系統(tǒng)的動力參數(shù)不會對顫振因子產(chǎn)生顯著影響,由此還可以推知,在攻角3°和5°條件下,發(fā)生顫振時作用在模型上的實際氣動力明顯受到了某個因素的影響,導致顫振導數(shù)受到了影響.顫振過程中,能夠影響顫振風速的有振幅、結(jié)構(gòu)阻尼比和運動相位角3個因素.

        3.2 振幅的影響

        圖5列舉出了工況1對應(yīng)的顫振臨界狀態(tài)下扭轉(zhuǎn)運動和豎向運動的時程曲線.

        (a) 風攻角0°

        (b) 風攻角3°

        (c) 風攻角5°

        圖4 顫振因子試驗值與計算值對比Fig.4 Comparison of RC of flutter factors

        從圖5中可以看出,豎向運動和扭轉(zhuǎn)運動的位移均很小,扭轉(zhuǎn)振幅均在1°以內(nèi),豎向振幅在2 mm以內(nèi).其他未列舉工況的振幅也在此范圍內(nèi).根據(jù)Noda等[14]對于寬高比為13的矩形斷面顫振導數(shù)的試驗研究結(jié)論,在3°以內(nèi)的扭轉(zhuǎn)運動下和6.5 mm的豎向運動下,測試斷面的顫振導數(shù)值保持不變.由此可以得出,試驗中的顫振振幅不會影響顫振導數(shù)的取值.再討論高次諧波的影響:根據(jù)王騎等[12]的強迫振動風洞試驗結(jié)論,扁平箱梁在8°扭轉(zhuǎn)振幅以內(nèi)和20 mm豎向振幅以內(nèi),氣動力中沒有明顯的高次諧波分量存在,由此排除了高次諧波分量對顫振的影響.

        圖5 顫振響應(yīng)時程圖Fig.5 Time-history response of flutter

        3.3 結(jié)構(gòu)阻尼比的影響

        表4 不同結(jié)構(gòu)阻尼比下顫振臨界風速計算結(jié)果Tab.4 Critical velocities under different damping ratios m/s

        結(jié)構(gòu)阻尼比從0.2%變化到0.5%,顫振臨界風速的最大相對變化僅為3.3%,對應(yīng)的顫振因子最大變化值也為 3.3%.由此可知,結(jié)構(gòu)阻尼比的顯著變化不會引起顫振風速的顯著變化,使其離散性顯著增加,因此阻尼比不是扁平箱梁顫振臨界風速試驗值與計算值產(chǎn)生較大差異的原因.

        表5 不同結(jié)構(gòu)阻尼比下顫振因子計算結(jié)果Tab.5 Flutter factors under different damping ratios

        3.4 相位角的影響

        最后討論顫振時相位角的影響.圖6給出了在不同的攻角下,顫振因子及相位差隨風攻角的變化曲線.從圖6中可以看出:在攻角0°時,顫振因子幾乎不隨相位差的變化而變化,其所對應(yīng)的相位差集中在區(qū)間[-0.55,-0.35],跨度為0.2;在攻角3°時,相位差集中在區(qū)間[-0.5,-0.25],跨度為0.25;在攻角5°時,相位差集中在區(qū)間[-0.45,0.03],跨度為0.5.因此可以得到這樣的結(jié)論:隨著攻角增大,發(fā)生顫振時的相位差變化區(qū)間增大,且在攻角3°和5°條件下,相位差越大,顫振因子的值越大,即顫振因子與相位差成正變關(guān)系.

        圖6 顫振因子及相位差關(guān)系圖Fig.6 Comparison of flutter factor and phase difference

        基于以上分析可以得到:顫振時模型斷面所受到的自激氣動力受到了相位差影響而發(fā)生了變化,并由此引起了顫振導數(shù)的變化,即顫振時模型對應(yīng)的實際顫振導數(shù)不同于單自由度振動下測試的顫振導數(shù).事實上,早在1993年,顫振運動相位角對自激氣動力的影響就被Matsumoto等[15]在風洞試驗中發(fā)現(xiàn)和證實.2015年,Lee等[16]通過機翼的強迫振動試驗證實了相位角對氣動力的顯著影響.2016年,Liao等[17]和Wang等[13]分別基于數(shù)值計算和強迫風洞試驗,發(fā)現(xiàn)了同頻耦合運動中的扁平箱梁和矩形斷面的自激氣動力是隨著耦合相位角的變化而變化的,且不等于兩個單自由度下氣動力的疊加.

        因此可以推測,在攻角3°和5°下,由于顫振時相位角的變化導致的顫振導數(shù)變化可能是計算顫振風速和試驗顫振風速產(chǎn)生差異的原因.

        4 結(jié)論和建議

        通過對比扁平箱梁在不同攻角下顫振臨界風速的理論計算分析與風洞試驗驗證,發(fā)現(xiàn)了大攻角下兩者存在較大的差異,得到如下主要結(jié)論:

        (1) 在攻角0°下,顫振臨界風速計算值及試驗值吻合,對于不同的測試工況,兩者差異均在5%以內(nèi),但隨著攻角增大,計算值與試驗值間的誤差也增大,在攻角5°下最大達到10.8%.

        (2) 在攻角3°和5°下,通過對理論顫振因子和試驗顫振因子差異的分析,可以推斷出耦合顫振條件下,相位角的變化可能引起了顫振導數(shù)的變化,從而導致計算顫振風速和試驗值產(chǎn)生差異.

        基于本文所取得的各項研究數(shù)據(jù)和結(jié)論,提出如下建議:有必要進一步開展大攻角條件下彎扭耦合顫振相位差對扁平箱梁斷面顫振導數(shù)帶來的影響,以修正大攻角下的顫振臨界風速.

        [1] WILDE K, FUJINO Y, MASUKAWA J. Time domain modeling of bridge deck[J]. Structural Engineering Earthquake Engineering, 2010, 13(2): 93-104

        [2] MASUMOTO M. Aerodynamic damping of prisms[J]. Journal of Wind Engineering and Industrial Aerodynamics, 1996, 59(2/3): 159-175.

        [3] CHEN X, MASUMOTO M, KAREEM A. Time domain flutter and buffeting response analysis of bridges[J]. Journal of Engineering Mechanics, 2000, 126(1): 7-16.

        [4] CHEN X. Improved understanding of bimodal coupled bridge flutter based on closed-form solution[J]. Journal of Engineering Mechanics, 2007, 133(1): 22-31.

        [5] SCANLAN R, TOMKO J. Airfoil and bridge deck flutter derivatives[J]. Journal of Engineering Mechanics, 1971, 97: 1717-1737.

        [6] SARKAR P P, JONES N P, SCANLAn R H. Identification of aeroelastic parameters of flexible bridges[J]. Journal of Engineering Mechanics, 1994, 120(8): 1718-1742.

        [7] MATSUMOTO M, NIIHARA Y, KOBAYASHI Y, et al. Flutter mechanism and stabilization of bluff bodies[C]∥Proceedings of the Ninth ICWE. New Delhi: Wiley Eastern Limited, 1995: 827-838.

        [8] CHEN X, KAREEM A. Efficacy of the implied approximation in the identification of flutter derivatives[J]. Journal of Engineering Mechanics, 2004, 130(12): 2070-2074.

        [9] ARGENTINI T, DIANA G, ROCCHI D, et al. A case-study of double multi-modal bridge flutter: Experimental result and numerical analysis[J]. Journal of Wind Engineering & Industrial Aerodynamics, 2016, 151: 25-36.

        [10] 丁泉順. 大跨度橋梁耦合顫抖振響應(yīng)的精細化分析[D]. 上海:同濟大學,2001.

        [11] GE Y J, XIANG H F. Computational models and methods for aerodynamic flutter of long-span bridges[J]. China Civil Engineering Journal, 2008, 96(10/11): 1912-1924.

        [12] 王騎,廖海黎,李明水,等. 大振幅下薄翼和流線型箱梁的氣動遲滯研究[J]. 實驗流體力學,2013,27(1): 32-37.

        WANG Qi, LIAO Haili, LI Mingshui, et al. Aerodynamic hysteresis of thin airfoil and streamline box girder under large amplitude oscillation[J]. Journal of Experiment in Fluid Mechanics, 2013, 27(1): 32-37.

        [13] WANG Qi, LIAO Haili, LI Mingshui, et al. Coupling and nonlinearity and span-wise correlation in aerodynamic force of a rectangular cylinder[C]∥The 8th International Colloquium on Bluff Body Aerodynamics and Applications. Boston: [s.n.], 2016: 1-8.

        [14] NODA M, UTSUNOMIYA H, NAGAO F, et al. Effects of oscillation amplitude on aerodynamic derivatives[J]. Journal of Wind Engineering and Industrial Aerodynamics, 2003, 91(1): 101-11.

        [15] MASUMOTO M, SHIRAISHI N, SHIRATO H, et al. Aerodynamic derivatives of coupled/hybrid flutter of fundamental structural sections[J]. Journal of Wind Engineering and Industrial Aerodynamics, 1993, 49(1/2/3): 575-584.

        [16] LEE T, SU Y Y. Surface pressures developed on an airfoil undergoing heaving and pitching motion[J]. Journal of Fluids Engineering, 2015, 137(5): 051105.

        [17] LIAO H, WAN J W, WANG Q, et al. Numerical study on nonlinear and motion coupling effects on self-excited force of a bridge deck[C]∥The 8th International Colloquium on Bluff Body Aerodynamics and Applications. Boston: [s.n.], 2016: 1-8.

        猜你喜歡
        氣動力風洞試驗攻角
        飛行載荷外部氣動力的二次規(guī)劃等效映射方法
        風標式攻角傳感器在超聲速飛行運載火箭中的應(yīng)用研究
        大攻角狀態(tài)壓氣機分離流及葉片動力響應(yīng)特性
        低風壓架空導線的風洞試驗
        電線電纜(2017年5期)2017-10-18 00:52:03
        側(cè)風對拍動翅氣動力的影響
        滾轉(zhuǎn)機動載荷減緩風洞試驗
        附加攻角效應(yīng)對顫振穩(wěn)定性能影響
        振動與沖擊(2015年2期)2015-05-16 05:37:34
        民用飛機攻角傳感器安裝定位研究
        遮擋條件下超高層建筑風洞試驗研究
        重慶建筑(2014年12期)2014-07-24 14:00:32
        高速鐵路接觸線覆冰后氣動力特性的風洞試驗研究
        久久久久亚洲av综合波多野结衣| 日韩精品大片在线观看| 日韩精品中文字幕无码一区 | 中文字幕乱码亚洲在线| 麻豆md0077饥渴少妇| 国外精品视频在线观看免费| 丝袜美女污污免费观看的网站| 亚洲熟妇中文字幕日产无码| 日本女同av在线播放| 国产大屁股视频免费区| 国产精品人妻一码二码尿失禁| 免费啪啪视频一区| 亚洲A∨日韩Av最新在线| 丝袜美腿在线观看视频| 国产精品亚洲片在线观看不卡| 欧美狠狠入鲁的视频777色| 亚洲自偷自偷偷色无码中文| 不卡国产视频| 91中文在线九色视频| 久久综合伊人77777麻豆| aaa级久久久精品无码片| 精品久久久久久无码不卡 | 国产少妇高潮在线视频| 四虎国产成人永久精品免费| 麻豆亚洲av永久无码精品久久| 怡春院欧美一区二区三区免费| 中文字幕久久久久久精| 大屁股流白浆一区二区| 精品国产av一区二区三区四区| 人人人妻人人澡人人爽欧美一区| 麻豆av一区二区三区| 国内精品久久久久久久久久影院| 9久久精品视香蕉蕉| 日本顶级片一区二区三区| 亚洲色一区二区三区四区| 55夜色66夜色国产精品视频| 99热国产在线| 青青草免费观看视频免费| 国产亚洲一区二区在线观看| 欧美黑人疯狂性受xxxxx喷水| 精品国产18禁久久久久久久|