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

考慮幾何非線性的橋梁后顫振極限環特性

2018-09-26 11:31:58吳長青張志田張偉峰陳政清
湖南大學學報·自然科學版 2018年5期

吳長青 張志田 張偉峰 陳政清

摘 要:根據橋梁斷面的試驗顫振導數,采用階躍函數模擬了自激力的時域表達式,并推導了便于時域分析的自激力遞推公式.采用APDL語言編制了顫振時域分析的程序并在ANSYS中實現.時域分析表明,幾何非線性效應對橋梁顫振臨界狀態影響甚微,而對其后顫振性能影響很大.線性理論揭示的后顫振響應是一種典型的發散現象,而計入幾何非線性效應后,后顫振響應最終演變為小振幅的極限環振動(LCO).此外,能量分析表明,線性發散振動的結構儲能不斷增加,而考慮幾何非線性時,結構儲能維持在一個較低水平(LCO狀態).相比線性發散造成的災難性毀滅而言, LCO只會對結構產生累積損傷;鑒于此,還需要綜合考慮材料的強度及疲勞特性等因素,才能進一步評估橋梁結構的安全性與穩定性.

關鍵詞:后顫振;時域分析;極限環;幾何非線性;能量

中圖分類號:U448.25 文獻標志碼:A

Limit Cycle Oscillation of Bridge Post-flutter with

Geometric Nonlinearities Included

WU Changqing, ZHANG Zhitian?, ZHANG Weifeng, CGEN Zhengqing

(Wind Engineering Research Center, Hunan University, Changsha 410082, China)

Abstract: According to experimental sectional flutter derivatives, indicial functions are applied to simulate the self-excited loads of a bridge deck section, and their recursive formulas which are essential in the implementation of FE analysis procedure are given. The procedure of time-domain flutter analysis, which is achieved by APDL language, is performed by the ANSYS software. Numerical results show that geometric nonlinearities have a negligible effect on the flutter threshold, but a significant effect on the post-flutter properties. When geometric nonlinearities are included, the post-flutter ultimately leads to a limit cycle oscillation (LCO) with a very small amplitude compared to a formidable divergence resulted from a linear method. Furthermore, the analysis shows that, in the case of linear analysis, the energy stored in the structure increases continuously as time progresses; however, this energy is limited in a quite low level (LCO state) when geometric nonlinearities are involved. Compared with the traditional divergence and catastrophic collapses, LCO results in cumulative damage that is of much more moderate; in view of this, other factors, such as the material strength and fatigue properties, are indispensable for further evaluation of the security and stability of the bridge structure.

Key words: post-flutter; time-domain analysis; limit cycle; geometric nonlinearities; energy

大跨度橋梁風致失穩現象包括靜力失穩與動力失穩兩類,其中靜力失穩表現為側向屈曲與扭轉發散[1-2],而動力失穩則表現為顫振、渦激共振、馳振及抖振等引起的失穩現象.顫振問題是大跨度橋梁風致振動研究中的熱點問題之一.多年來,抗風設計者圍繞顫振臨界風速這一問題進行了廣泛地研究,建立了一套成熟的顫振頻域分析方法[3-7].然而,這類方法屬于線性理論的范疇,局限于顫振臨界問題的分析,無法準確地分析與評估大跨度橋梁結構的后顫振特性.基于線性理論可知,當風速高于顫振臨界風速時,結構將會出現大變形的發散振動,這顯然是一種非線性振動.然而,對于非線性顫振研究而言,頻域法不再適用.因此,要準確地評估橋梁后顫振性能,必須建立一套非線性顫振的分析方法.近年來,機翼氣動問題的相關研究表明,當同時考慮結構及氣動力非線性或者只考慮其中之一時,機翼的后顫振行為最終表現為極限環振動(LCO)[8-10].然而,至今為止,大跨度橋梁的非線性顫振與后顫振研究的相關報道較少.究其原因主要有以下兩點:第一,抗風設計者目前最關心的問題依然是結構的顫振臨界問題,而這一問題在線性范圍內即可得到很好地解決;第二,建立一套成熟的橋梁非線性氣動彈性理論存在很大困難,目前這方面的研究還處于起步階段.

橋梁后顫振研究中所涉及的非線性主要包括結構的非線性與氣動力非線性;其中結構的非線性包括幾何非線性與材料非線性,而氣動力非線性則體現在自激力與運動的振幅及頻率相關[11]. Tang等[10]、Attar與Dowell等[8]、Shams等[12]在僅考慮結構的非線性效應的情況下,分別揭示了不同類型機翼的顫振LCO現象.在氣彈問題的研究中,通常引入半經驗模型來描述機翼斷面的氣動力非線性特性,如Tran和Petot[13]提出了著名的ONERA動力失速模型,它被廣泛應用于研究直升機槳葉的動態失速問題.隨后,很多學者采用ONERA模型并考慮結構的非線性特性研究了一些機翼結構的非線性氣彈問題[14-17].Leishman等[18]和Larsen等[19]也分別提出了用于研究動態失速問題的半經驗模型.然而,上述提及的半經驗模型均是基于機翼斷面的氣動力學理論建立起來的,很難直接用于橋梁鈍體斷面的氣彈問題研究.假如橋梁的后顫振行為也表現為LCO,那么橋梁破壞與否最終由材料的強度、疲勞特性以及風荷載持續時間等內外因素共同決定.

顫振時域法由于方便考慮結構及氣動力非線性效應的影響,而被廣泛用于非線性氣彈問題的研究中.顫振導數是描述橋梁結構氣動自激力的參數,一般可通過節段模型風洞試驗或CFD數值模擬得到[20].對于時域法而言,應先基于顫振導數擬合出與SCANLAN自激力模型等效的自激力時程[21-23].本文基于某大橋的試驗顫振導數,采用MATLAB遺傳優化算法擬合得到了用于自激力時域表達的階躍函數,并基于ANSYS軟件實現了大跨度橋梁的后顫振時域分析,評估了結構的幾何非線性對橋梁顫振及后顫振性能的影響,并且從能量的角度出發,探究了結構的儲能與顫振穩定性的關系.氣動力非線性對結構后顫振性能的影響將在今后的研究中探討.

1 氣動自激力描述

1.1 SCANLAN自激力模型

鈍體橋梁的非定常自激力一般采用SCANLAN提出的時頻混合模型來描述,模型中引入8個顫振導數,將主梁單位長度的自激升力與升力矩近似描述為運動狀態的線型函數[11,24],其表達式如下:

分別為豎向位移與扭轉位移對時間t的導數.為了簡化,式(1)與(2)中暫未考慮側向運動的貢獻及氣動彈性阻力的影響.

1.2 自激力時域模型

1.2.1 階躍函數及其參數擬合

SCANLAN自激力模型不能直接用于時域分析,因此在進行顫振分析之前需要將SCANLAN自激力時域化,其方式主要有階躍函數與有理函數兩類方法.由于階躍函數氣動自激力模型的平均(或極限)特性具有明顯的物理意義,其收斂的平均值即代表結構真實的平均響應,而有理函數的方法不具備這樣的優勢[25],因此本文采用階躍函數的方法.階躍函數在經典機翼理論中常被用來描述風速或風攻角突然變化而引起的氣動升力瞬態演變過程,表達為:

. (3)

式中:

由表3可知,數值分析與風洞試驗得到的顫振臨界風速與顫振頻率相差較小,即表明本文的時域分析方法是可靠的.圖5為不同風速作用下主梁中點的顫振響應時程曲線.觀察這些曲線,可得到以下兩點結論:1)不管幾何非線性考慮與否,求解的顫振臨界風速及顫振頻率相差很小,這表明對于顫振臨界問題的求解可忽略幾何非線性效應的影響;2)當風速超過顫振臨界值時,線性分析揭示的后顫振響應不斷增大,主梁發生了明顯的扭轉發散現象,而非線性分析對應的后顫振響應則表現為等幅的LCO,且其振幅較小;此外,隨著風速的進一步增大,LCO的振幅有所增加,并且結構進入LCO狀態所需的時間也縮短了,如圖5(f)與(h)所示.綜上可知,結構的幾何非線性效應抑制了橋梁結構的后顫振扭轉響應,將發散振動最終演變為具有穩定幅值的LCO.對于90 m/s的風速(如圖5(h)所示),考慮幾何非線性時主梁中點的扭轉LCO幅值恒為2.5°,而線性情形對應的扭轉響應呈指數增長,其幅值在800 s時就高達15°.由此可見,LCO遠不足以導致大跨度橋梁結構發生失穩破壞.圖6給出了主梁中點的顫振扭轉振動相平面圖,由圖可知,未考慮幾何非線性效應時,扭轉振動的軌線不封閉成環,軌線不斷往外擴張,這是一種典型的不穩定的發散現象;而考慮幾何非線性效應時,扭轉振動的軌線最終處于一個穩定的封閉軌道上,即表明結構達到了穩定的LCO狀態.圖7所示為顫振臨界狀態與后顫振LCO狀態下主梁中點的彎扭時程曲線,由圖可知,橋梁結構的顫振表現為耦合彎扭振動,且振動頻率及振動步調完全一致.

因結構變形引起結構剛度改變的一類問題都屬于幾何非線性問題.幾何非線性通常分為大應變、大位移和應力剛化三類,其中大應變和大位移中均包括了應力剛化.應力剛化效應是指在應變位移關系中考慮了位移的二階非線性以及前一應力狀態對當前應力狀態的影響.考慮應力剛化效應時,程序會自動計算應力剛化矩陣并將其添加到結構的剛度矩陣中.對于桁架、梁和殼體單元,在大變形分析中應使用應力剛化,否則得不到精確解.基于ANSYS平臺的顫振時域分析,是通過設置NLGEOM,ON命令來考慮結構的幾何非線性的,當NLGEOM,ON命令打開時,應力剛化效應也默認打開.ANSYS中結構的切線剛度矩陣通常包括二部分,其表達式為:

K=K0+Kσ. (24)

式中:K0為小位移剛度矩陣;Kσ為應力剛度矩陣,即幾何剛度矩陣,它在每次迭代中都是變化的,因而它是非線性的.對于線性情形,則關閉大變形開關(即設置NLGEOM,OFF命令).

如果結構的總體剛度由于幾何非線性的影響而增加,剛度的增加將會限制結構的振動,反過來制約了自激力的發展,可能正是由于剛度的非線性特性引起了橋梁結構的后顫振LCO現象.Dowell在平板顫振研究中認為拉伸應變能引起的非線性膜力是觸發平板結構LCO的原因,即隨著振動幅值的增加,拉伸應變剛度的增加限制了振動的幅值[28].此外還有一些研究[29-30]也指出,結構系統的剛度非線性特性可能是引起顫振極限環振動的原因之一.然而,幾何非線性觸發橋梁結構后顫振LCO的內在機理目前仍不明確,作者暫時難以給出一種合理的方法來驗證上述猜想.

需要指出的是,由于本文的顫振分析沒有考慮除自激力之外的其他氣動荷載(如平均風荷載、抖振力等)的作用,因此當風速小于或等于顫振臨界風速時,結構的響應依然較小,此時幾何非線性對顫振臨界狀態的影響并不明顯.然而,當風速超過顫振臨界值時,結構的變形逐漸增大,幾何非線性的影響則主要體現在橋梁的后顫振響應中.

對于顫振系統而言,振動能量主要通過結構阻尼與氣動阻尼耗散.氣流輸入結構的總能量減去系統阻尼耗散的能量后剩余的能量為結構儲存的能量(包括動能與彈性應變能兩部分).結構的儲能越大,結構的穩定性越差.因此,可通過考察顫振過程中結構儲能的變化情況來評估橋梁結構的顫振穩定性.

以扭轉自激振動(即式(31)表達的能量方程)為例,考察大跨度橋梁在扭轉振動中的結構儲能的變化情況.圖8給出了橋梁主跨跨中附近單位長度主梁在振動過程中的結構儲能的時程曲線.當風速低于顫振臨界值(U=60 m/s)時,不論是線性情形還是非線性情形,結構儲能均隨著時間的推移逐漸衰減至零,如圖8(a)與(b)所示;這表明振動系統的阻尼耗散了結構振動的能量,結構的響應衰減至零.當風速達到臨界值時,結構儲能曲線呈現出小振幅的極限環形式,且幅值保持不變,如圖8(c)與(d)所示;這表明阻尼耗散的能量與結構從氣流中獲得的能量最終達到了動態平衡狀態,即在完整的振動周期內,氣流對系統的儲能貢獻為零,最終使得結構呈現出一種穩定的LCO狀態.當風速大于臨界值(U=90 m/s)時,線性情形對應的結構儲能隨著時間的推移而呈指數式增長,如圖8(e)所示,這意味著氣動自激力對結構總體上做正功,氣流不斷向結構輸入能量,結構的響應不斷增大,這對應于一種顫振發散現象;然而,非線性情形對應的結構儲能隨著時間的推移最終呈現出極限環的變化形式,如圖8(f)所示,這與顫振臨界狀態下的情形類似,即結構的儲能最終也處于一種動態平衡狀態.

(a)線性; U=60 m/s (b)幾何非線性; U=60 m/s

(c)線性; U=64.6 m/s (d)幾何非線性; U=64.8 m/s

(e)線性; U=90 m/s (f)幾何非線性; U=90 m/s

圖8 結構儲能時程曲線

Fig.8 Histories of the energy stored in the structure

3 結論

根據試驗顫振導數,采用MATLAB遺傳優化算法擬合得到了用于自激力時域表達的階躍函數,并基于ANSYS平臺實現了大跨度橋梁結構的顫振時域分析,探究了結構的幾何非線性對大跨度橋梁顫振臨界狀態及后顫振行為的影響.通過對某大跨度橋梁進行顫振分析,得到如下結論:

1)幾何非線性效應考慮與否,得到的顫振臨界風速及顫振頻率幾乎是一致的.這表明橋梁的顫振臨界狀態在線性范圍內即可得到準確解,幾何非線性的影響可以忽略不計.

2)幾何非線性效應對橋梁后顫振行為的影響顯著,線性方法揭示的后顫振響應為振幅無限增大的發散現象,而考慮幾何非線性時的后顫振行為則表現為較小且穩定幅值的LCO.此外,橋梁后顫振LCO下的豎向與扭轉振動頻率相等且振動步調完全一致.

3)從結構儲能的時程曲線可知,考慮幾何非線性效應時,后顫振對應的結構儲能曲線最終呈現出幅值較小的極限環形式;這表明系統阻尼耗散的能量與結構從氣流中獲得的能量最終達到了一個動態平衡狀態,即在完整的周期內,氣流對系統的儲能貢獻為零,最終使得結構呈現出一種穩定的LCO狀態;當不考慮幾何非線性時,結構儲能曲線呈指數式增大,即隨著振動周期的推移,氣流不斷地向結構輸入能量,結構響應不斷增大,這對應于顫振發散狀態.

4)較小振幅的穩態LCO并沒有傳統線性方法下的扭轉發散可怕,因為它并不能在短時間內致使結構發生災難性破壞;鑒于此,要進一步評估橋梁結構的穩定性與安全性,應綜合考慮風荷載的變化特性與材料的強度、疲勞特性以及構件薄弱環節等結構強健性因素.

需要說明的是,本文的顫振分析暫未考慮氣動力非線性的影響,這是由于準確描述橋梁斷面的真實氣動力非線性特性目前依然存在較大困難.然而,有關機翼及平板斷面的后顫振研究表明,氣動力非線性可顯著地影響結構的后顫振性能.因此,綜合考慮結構幾何與氣動力非線性的橋梁后顫振特性研究將是一項具有挑戰且意義重大的課題.

參考文獻

[1] 吳長青, 張志田, 陳政清. 懸索橋靜風扭轉發散的影響因素研究[J]. 湖南大學學報(自然科學版),2016, 43(3): 15-22.

WU C Q, ZHANG Z T, CHEN Z Q. Research of influencing factors on aerostatic torsional divergence of suspension bridges[J]. Journal of Hunan University(Natural Sciences), 2016, 43(3): 15-22. (In Chinese)

[2] 吳長青, 張志田. 懸索橋的靜風扭轉發散有限元精細化分析[J]. 湖南大學學報(自然科學版), 2016, 43(9): 88-97.

WU C Q, ZHANG Z T. Refined analysis of finite element for torsional divergence of suspension bridges[J]. Journal of Hunan University (Natural Sciences), 2016, 43(9): 79-88. (In Chinese)

[3] AGAR T J A. Aerodynamic flutter analysis of suspension bridges by a model technique[J]. Engineering Structures, 1989, 11(2): 75-82.

[4] NAMINI A, ALBRECHT P, BOSCH H. Finite element-based flutter analysis of cable-suspended bridge[J]. Journal of Structural Engineering, 1992, 118(6): 1509-1526.

[5] CHEN Z Q, AGAR T J A. Finite element-based flutter analysis of cable-suspended bridge[J]. Journal of Structural Engineering, 1994, 120(3): 1044-1046.

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

[7] KATSUCHI H, JONES N P, SCANLAN R H, et al. Multimode coupled flutter and buffeting analysis of the Akashi-Kaikyo bridge[J]. Journal of Structural Engineering, 1999, 125(1): 60-70.

[8] ATTAR P J, DOWELL E H, TANG D M. A theoretical and experimental investigation of effects of a steady angle of attack on the nonlinear flutter of a delta wing plate model[J]. Journal of Fluids and Structures, 2003, 17(2): 243-259.

[9] TANG D M, DOWELL E H. Effects of geometric structural nonlinearity on flutter and limit cycle oscillation of high-aspect-ratio wings[J]. Journal of Fluids and Structures, 2004, 19(3): 291-306.

[10] TANG D M, DOWELL E H, HENRY J K. Limit cycle oscillations of delta wing models in low subsonic flow[J]. Journal of AIAA, 1999, 37(11):1355-1362.

[11] SCANLAN R H. Problematics in formulation of wind-force models for bridge decks[J]. Journal of Engineering Mechanics, 1993, 119(7): 1353-1375.

[12] SHAMS S H, LAHIDJANI M H S, HADDADPOUR H. Nonlinear aeroelastic response of slender wings based on Wagner function[J]. Journal of Thin-Walled Structures, 2008, 46(11): 1192-1203.

[13] TRAN C T, PETOT D. Semi-empirical model for the dynamic stall of airfoils in view of the application to the calculation of responses of a helicopter blade in forward flight[J]. Vertica, 1981, 5(1): 35-53.

[14] TANG D M, DOWELL E H. Nonlinear aeroelasticity in rotorcraft[J]. Journal of Mathematical and Computer Modelling, 1993, 18(3/4): 157-184.

[15] TANG D M, DOWELL E H. Effects of geometric structural nonlinearity on flutter and limit cycle oscillations of high-aspect-ratio wings[J]. Journal of Fluids and Structures, 2004, 19(3): 291-306.

[16] SARKAR S, BIJL H. Nonlinear aeroelastic behavior of an oscillating airfoil during stall-induced vibration[J]. Journal of Fluids and Structures, 2008, 24(6): 757-777.

[17] STANFORD B, BERAN P. Direct flutter and limit cycle computations of highly flexible wings for efficient analysis and optimization[J]. Journal of Fluids and Structures, 2013, 36: 111-123.

[18] LEISHMAN J G, BEDDOES T S. A semi-empirical model for dynamic stall[J]. Journal of the American Helicopter Society, 1986, 34(3): 3-17.

[19] LARSEN J W, NIELSEN S R K, KRENK S. Dynamic stall model for wind turbine airfoils[J]. Journal of Fluids and Structures, 2007, 23(7): 959-982.

[20] 祝志文,夏昌.基于兩種湍流模型的橋梁顫振導數識別研究及比較 [J].湖南大學學報:自然科學版,2010,37(11):6-11.

ZHU Z W, XIA C. Comparative study of two turbulence models based on the identification of flutter derivatives of bridges[J]. Journal of Hunan University(Natural Sciences), 2010, 37(11): 6-11. (In Chinese)

[21] FARSANI H Y, VALENTINE D T, ARENA A, et al. Indicial functions in the aeroelasticity of bridge decks[J]. Journal of Fluids and Structures, 2014, 48: 203-215.

[22] ZHANG Z T, CHEN Z Q, CAI Y Y, et al. Indicial functions for bridge aero-elastic forces and time-domain flutter analysis[J]. Journal of Bridge Engineering, 2011, 16(4): 546-557.

[23] MIRANDA S D, PATRUNO L, UBERTINI F, et al. Indicial functions and flutter derivatives: a generalized approach to the motion-related wind loads[J]. Journal of Fluids and Structures, 2013, 42: 466-487.

[24] SCANLAN R H, TOMKO J J. Airfoil and bridge deck flutter derivatives[J]. Journal of the Engineering Mechanics Division, 1971, 97(6): 1717-1737.

[25] 張志田, 陳政清,李春光.橋梁氣動自激力時域表達式的瞬態與極限特性[J].工程力學,2011, 28(2): 75-85.

ZHANG Z T, CHEN Z Q, LI C G. Limiting and transient characteristics of time-domain expressions for bridge self-excited aerodynamic forces[J]. Engineering Mechanics, 2011, 28(2): 75-85. (In Chinese)

[26] CARACOGLIA L, JONES N P. A methodology for the experimental extraction of indicial function for streamlined and bluff deck sections[J]. Journal of Wind Engineering and Industrial Aerodynamics, 2003, 91(5): 609-636.

[27] COSTA C, BORRI C. Application of indicial functions in bridge deck aeroelasticity[J]. Journal of Wind Engineering and Industrial Aerodynamics, 2006, 94(11): 859-881.

[28] DOWELL E H, HOPKINS M A. Limited amplitude panel flutter with a temperature differential//Proceedings of the 35th Structures, Structural Dynamics and Materials Conference. Washington, DC: AIAA, 1994: 1343-1355.

[29] ONEIL T, GILLIATT H C, STRGANAC T W. Investigations of aeroelastic response for a system with continuous structural nonlinearities[C]//Proceedings of the 37th Structures, Structural Dynamics and Materials Conference. Salt Lake City, UT, USA: AIAA, 1996: 1-9.

[30] ONEIL T. STRGANAC T W. Aeroelastic response of a rigid wing supported by nonlinear springs[J]. Journal of Aircraft, 1998, 35(4): 616–622.

主站蜘蛛池模板: 伊人久久福利中文字幕| 黄色一及毛片| 久久天天躁狠狠躁夜夜2020一| 91精品久久久久久无码人妻| 天天色综网| 无码视频国产精品一区二区| 欧美在线视频a| 国产99热| 国产精品无码久久久久久| 国产精品福利在线观看无码卡| 国产日韩欧美中文| 五月婷婷伊人网| 亚洲成在线观看 | 亚洲精品va| 国产免费a级片| 国产一线在线| 亚洲AV无码久久精品色欲| 国产视频欧美| 亚洲AV无码一区二区三区牲色| 一区二区理伦视频| 日韩在线永久免费播放| 一级毛片在线播放免费观看| 伊人激情综合网| 狼友av永久网站免费观看| 中文字幕调教一区二区视频| 人妻中文字幕无码久久一区| 国产美女精品一区二区| 久久视精品| 久久这里只精品国产99热8| 日韩午夜片| 国产香蕉97碰碰视频VA碰碰看| 免费三A级毛片视频| 亚洲视频一区在线| 99er精品视频| 狠狠干综合| 亚洲二三区| 国产精品主播| 青青国产成人免费精品视频| 午夜国产小视频| 亚洲成人网在线播放| vvvv98国产成人综合青青| 99在线观看免费视频| 午夜老司机永久免费看片| 国产精品无码久久久久久| 欧美丝袜高跟鞋一区二区| 精品国产Av电影无码久久久| 一本大道香蕉久中文在线播放 | 麻豆精品在线播放| 综合人妻久久一区二区精品 | 国产精品私拍在线爆乳| 亚洲黄色片免费看| 国产91丝袜| 亚洲精品无码在线播放网站| 日本精品视频一区二区| 就去吻亚洲精品国产欧美| 亚洲二区视频| 麻豆精选在线| 亚洲性视频网站| 波多野结衣在线se| 911亚洲精品| 97视频在线观看免费视频| 青青草原国产| 五月天天天色| 亚洲国产精品VA在线看黑人| 超清无码一区二区三区| 中文字幕自拍偷拍| 国产综合精品一区二区| 波多野结衣AV无码久久一区| 亚洲日本中文字幕乱码中文| 日韩一级二级三级| 青青草原国产免费av观看| 亚洲国产成人在线| 亚洲大尺码专区影院| 成人精品在线观看| 亚洲中文精品久久久久久不卡| 91青青草视频| 国产大片喷水在线在线视频 | 久久久久国色AV免费观看性色| 久久婷婷五月综合色一区二区| 国产十八禁在线观看免费| 亚洲Av综合日韩精品久久久| 国产成人精品亚洲77美色|