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

        ?

        橋梁斷面兩個渦振鎖定區(qū)間的數值模擬研究

        2016-08-04 06:13:33俐,帥,
        振動與沖擊 2016年11期
        關鍵詞:渦振尾流風洞試驗

        黃 俐, 周 帥, 梁 鵬

        (1. 華南農業(yè)大學 水利與土木工程學院,廣州 510642;2.湖南大學 風工程試驗研究中心,長沙 410082;3. 中交四航局第一工程有限公司,廣州 510310)

        橋梁斷面兩個渦振鎖定區(qū)間的數值模擬研究

        黃俐1, 周帥2, 梁鵬3

        (1. 華南農業(yè)大學 水利與土木工程學院,廣州510642;2.湖南大學 風工程試驗研究中心,長沙410082;3. 中交四航局第一工程有限公司,廣州510310)

        對于某類大寬高比橋梁斷面或者鈍體形式斷面,在基于同一組試驗參數的節(jié)段模型風洞試驗中,能夠實測到相同振型兩個分離的渦振鎖定區(qū)間現象,并且兩個區(qū)間內振動頻率一致,這與一個斷面對應一個Strouhal數的理論不相符合。為了進一步研究這類非常規(guī)振動形式的氣動機理,以一組寬高比為6的矩形斷面為研究對象,基于風洞試驗中實測的兩個獨立的渦振鎖定區(qū)間響應數據,采用流體動力學軟件Fluent開展了相應的數值模擬研究。數值計算獲取了與風洞試驗一致的兩個獨立分離的渦振鎖定區(qū)間風振曲線,并且在區(qū)間跨度以及幅值關系上均吻合良好,然后通過Fluent提取了前后兩個渦振鎖定區(qū)間內的氣動力和尾流漩渦進行了對比研究。研究結果表明,第一個鎖定區(qū)間內的尾流漩渦呈現出經典的卡門渦街形態(tài),第二個區(qū)間內的尾流渦模態(tài)則主要表現為非典型的“魚尾擺動”形態(tài),兩個渦振區(qū)間的尾流形態(tài)完全不同;在兩個獨立的鎖定區(qū)間內,氣動升力與位移響應之間始終存在著相位差,并且均隨著鎖定區(qū)間的發(fā)展而持續(xù)增大,第一個鎖定區(qū)間相位差的跳躍程度明顯大于第二個鎖定區(qū)間。

        橋梁;渦振;鎖定區(qū)間;數值模擬;尾流漩渦

        渦激共振是鈍體結構在來流作用下的一種常見的振動形式,是由于結構尾流中交替脫落的漩渦頻率與結構固有頻率重合而引起的一種流固耦合共振現象,其顯著特征就是共振區(qū)間鎖定和振動限幅[1-3]。針對特定的結構斷面,尾流的漩渦脫落頻率和來流風速之間可以通過Strouhal數進行換算,一種結構斷面對應一個Strouhal數,利用該參數值可以實現對結構渦振鎖定區(qū)間起振臨界風速的估算,這是目前關于渦振研究的一般性結論[4]。

        然而,對于某類大寬高比的橋梁斷面,例如我國的蘇通長江大橋的主梁斷面,在節(jié)段模型風洞試驗中,+5°風攻角的來流狀態(tài)下,隨著試驗風速的增加,實測到了節(jié)段模型系統扭轉自由度兩個分離的渦激共振鎖定區(qū)間,并且兩個鎖定區(qū)間內的振動卓越頻率均與模型固有頻率一致,也即說明蘇通橋主梁斷面在+5°風攻角來流下對應兩個不同的Strouhal數,這與一個截面對應一個Strouhal數的結論是相矛盾的[5-7]。

        對于蘇通橋主梁斷面風洞試驗中實測到的這類兩個鎖定區(qū)間的渦振現象,在大寬高比的矩形斷面的研究中也很容易觀察到。文獻[8]系統研究了寬高比在2.04.0的試驗工況中均能實測到兩個分離的渦振鎖定區(qū)間。

        Matsumoto[9]針對這類渦振現象開展研究,研究結論認為是結構表面運動渦與尾流卡門漩渦相互作用導致了多個渦振鎖定區(qū)間的出現,即,來流在結構斷面前沿分離,一部分在結構表面再覆并以一定的速度(大約為來流速度的60%)沿結構表面運動最終在尾流脫落,另一部分在尾流直接脫落形成卡門漩渦,運動渦的加入改變了原有卡門漩渦固有的脫落頻率,從而改變了結構斷面的Strouhal數,并且由于運動渦的脫落頻率與來流速度相關,因此,該研究認為同一種結構斷面在不同的來流速度下會對應不同的Strouhal數。

        Tamura[10]針對兩自由度尾流振子渦激力模型開展了相關研究,對渦激振動的氣動機理進行了數學的描述,并明確了數學模型中的每個參數均進行了物理意義的解釋。

        為了進一步研究這類振動的氣動機理,本文基于一組寬高比為6的矩形斷面開展研究,通過CFD數值模擬來觀察前后兩個分離的渦振鎖定區(qū)間內氣動力、位移響應以及尾流漩渦形態(tài)之間的對應關系。

        1模型參數與試驗結果

        針對寬高比為6的矩形斷面,Chen等[9]開展了剛性節(jié)段模型風洞試驗研究,模型的基本參數如表1所示。彈性懸掛剛性節(jié)段模型風洞試驗結果如圖1所示,從圖中可以看出兩個豎向自由度分離的渦振鎖定區(qū)間區(qū)間,并且試驗實測到兩個鎖定區(qū)間內的振動卓越頻率均與模型的固有頻率5.90 Hz基本保持一致。根據前后兩個渦振鎖定區(qū)間的起振風速,振動卓越頻率以及橫風向尺寸得出對應的兩個Strouhal數分別約為0.2和0.1,而第二個鎖定區(qū)間對應的Strouhal數0.1與歐洲抗風規(guī)范中的參考值0.06更為接近[12]。在渦振幅值關系上,前一個渦振鎖定區(qū)間的幅值約為后一個區(qū)間的0.5倍。

        圖1 風洞試驗實測風振響應曲線Fig.1 Wind tunnel tests measured wind induced oscillation curves

        參數類型數值長度L/mm1530橫風向尺寸D/mm40順風向尺寸B/mm240試驗風速U/(m·s-1)<6等效質量m/(Kg·m-1)3.25橫風向固有基頻f0/Hz5.90阻尼比ξ00.0058Scruton數Sc121Reynolds數Re6293

        2數值模擬方法

        CFD數值模擬基于流體動力學軟件Fluent展開,首先,采用Gambit進行計算區(qū)域的規(guī)劃、計算網格的劃分以及邊界條件的設置。如圖2(a)所示,矩形斷面的2D模型布置在計算域的中心位置,近壁面為四邊形剛性網格區(qū)域(包括四邊形正交貼體網格和外圍Pave網格,如圖2(b)所示),計算域外圍是完全正交的四邊形結構網格區(qū)域,介于外圍結構網格和近壁面剛性網格之間的是三角形動網格區(qū)域,整個計算域共生成了31 590個計算網格,壁面Yplus在3.5以內。

        采用Newmark-β方法對矩形斷面計算模型的動力特性進行模擬,首先,根據模型的質量、阻尼和剛度參數編寫模型振動的豎向和扭轉兩自由度迭代程序;然后以User Define Function (UDF)的方式在Fluent求解器中進行編譯。

        針對計算模型的Fluent流固耦合分析過程中,采用k-ωSST湍流模型,加載UDF,迭代計算按照以下步驟循環(huán)進行:第一步,對當前計算時間步模型壁面遍歷積分提取豎向和扭轉自由度氣動力;第二步,根據UDF程序計算出模型在當前時間步豎向和扭轉自由度的運動位移和速度;第三步,DEFINE-CG-MOTION宏強迫模型按當前時間步計算得到的結果產生相應位移并促使動網格區(qū)域的網格重劃分,形成下一時間步新的計算域。計算結果是否穩(wěn)定通過觀察迭代計算殘差和響應時程曲線來判斷。

        圖2 計算網格Fig.2 Calculation cells

        3計算結果

        3.1位移響應

        CFD數值計算得到的模型豎向自由度位移幅值響應曲線與試驗實測曲線對比如圖3所示。橫坐標為無量綱風速,f0為模型豎向自由度固有頻率5.90 Hz,從下文的時程曲線和頻譜中可以觀察到渦振鎖定區(qū)間內模型豎向自由度的振動卓越頻率與固有頻率是基本一致的,U為來流風速,D為模型橫風向特征尺寸,D=40 mm;縱坐標為無量綱幅值,即,位移響應幅值與模型橫風向特征尺寸的比值,位移響應幅值采用固定時間長度(55 s)位移時程曲線取根方差乘以根號2的方式確定。

        圖3 風振響應曲線試驗值與CFD計算值對比Fig.3 The comparison of wind induced oscillation curves between experimental and calculation results

        從圖3中可以看出,基于風洞試驗模型參數開展的數值計算同樣觀察到了兩個分離的渦振鎖定區(qū)間現象,并且在鎖定區(qū)間無量綱風速跨度和幅值關系上,計算值與試驗實測值吻合良好。為了能進一步研究前后兩個鎖定區(qū)間內氣動力以及尾流渦等參數的變化特征,特選定數值計算結果中6個特征響應點的計算結果進行分析。如圖3所示,L1、J1、K1對應第一個渦振鎖定區(qū)間的上升點、幅值點和下降點;L2、J2、K2則分別對應第二個渦振鎖定區(qū)間的三個區(qū)間特征點。

        另外,在動力自由度關系上,CFD數值計算模型與風洞試驗剛性節(jié)段模型是一致的,基于風洞試驗模型豎向和扭轉自由度的實際參數,數值計算模型同樣分別設置了豎向和扭轉兩個自由度。數值計算過程中也獲得了模型扭轉自由度的相關響應數據,由于模型在豎向和扭轉兩自由度振動狀態(tài)下,實際振動中豎向自由度響應絕對占優(yōu)。在現有的計算工況中扭轉自由度的位移響應幅值均在0.2°以下,基本可忽略,限于篇幅,不再詳細展開討論。

        第一個渦振鎖定區(qū)間內特征響應點L1、J1、K1豎向自由度位移響應時程和頻譜如圖4~圖6所示。從圖中可以看出,三個特征響應點的位移時程均為等幅狀態(tài),而振動卓越頻率分別為5.82 Hz、5.88 Hz和5.85 Hz,均與模型豎向自由度固有頻率5.90 Hz基本保持一致,可定性為等幅渦振。

        圖7~圖9分別對應第二個渦振鎖定區(qū)間內L2、J2、K2特征響應點。從圖中同樣可觀察到與第一個鎖定渦振區(qū)間一致的響應特征:位移響應時程基本為等幅狀態(tài),卓越頻率分別為5.86 Hz、5.98 Hz、5.98 Hz,與模型固有頻率5.90 Hz基本一致。

        特征響應點J1和 J2分別對應第一個和第二個渦振鎖定區(qū)間的渦振幅值點,無量綱幅值分別為18.0/1 000和30.5/1 000,對應來流風速分別為5.0 m/s和12.2 m/s,振動卓越頻率分別為5.88 Hz和5.98 Hz,并且對應同一模型斷面,橫風向特征尺寸均為40 mm,因此,按照Strouhal定律可以計算得到兩個不同的參數值,分別為0.2和0.08。

        3.2尾流漩渦形態(tài)

        3.2.1第一個鎖定區(qū)間

        第一個渦振鎖定區(qū)間內特征響應點L1、J1和K1對應風速下單個振動周期內尾流漩渦形態(tài)變化如圖10~圖12所示。起振點L1振幅較小,尾流漩渦脫落也不明顯;J1測點對應第一個渦振鎖定區(qū)間的幅值點,尾流呈現出了明顯的卡門漩渦脫落特征,單個振動周期內有兩個漩渦分別從模型的上下表面脫落,形成了豎向自由度的簡諧氣動升力;K1對應著鎖定區(qū)間的下降段,單個振動周期內,尾流漩渦脫落總體上也表現出了卡門渦街特征,但漩渦特征相對于J1測點已經明顯弱化[13]。

        3.2.2第二個鎖定區(qū)間

        第二個渦振鎖定區(qū)間單個振動周期內L2、J2、K2特征響應點的尾流變化分別如圖13~圖15所示。而從圖中可以觀察到,與第一個鎖定區(qū)間出現的卡門漩渦特征明顯不同的是,第二個鎖定區(qū)間內三個特征響應點的尾流特征沒有漩渦脫落,而是表現出“魚尾擺動”的形態(tài),尾流擺動幅度與模型振動位移幅值是正相關的關系。值得注意的是,同一個結構斷面出現兩個渦振鎖定區(qū)間通常發(fā)生在截面寬高比大于4的扁平的鈍體斷面。本文研究的矩形斷面寬高比為6,所觀測到的第一個鎖定區(qū)間為卡門渦街尾流形態(tài),而第二個鎖定區(qū)間內的尾流中卡門渦街不明顯而呈現出“魚尾擺動”,這實際上與文獻中Tamura兩自由度尾流振子模型中描述的渦激振動形成機理較為吻合,前后兩個渦振鎖定區(qū)間的Reynolds數相差將近兩倍,這可能就是尾流形態(tài)發(fā)生改變的原因之一[14-15]。

        圖4 測點L1豎向自由度響應Fig.4DynamicresponsesofverticalDOFofpointL1圖5 測點J1豎向自由度響應Fig.5DynamicresponsesofverticalDOFofpointJ1圖6 測點K1豎向自由度響應Fig.6DynamicresponsesofverticalDOFofpointK1

        圖7 測點L2豎向自由度響應Fig.7DynamicresponsesofverticalDOFofpointL2圖8 測點J2豎向自由度響應Fig.8DynamicresponsesofverticalDOFofpointJ2圖9 測點K2豎向自由度響應Fig.9DynamicresponsesofverticalDOFofpointK2

        圖10 單個振動周期內測點L1尾流漩渦形態(tài)變化規(guī)律Fig.10ThevorticeschangingrulesinoneoscillationcycleofpointL1圖11 單個振動周期內測點J1尾流漩渦形態(tài)變化規(guī)律Fig.11ThevorticeschangingrulesinoneoscillationcycleofpointJ1圖12 單個振動周期內測點K1尾流漩渦形態(tài)變化規(guī)律Fig.12ThevorticeschangingrulesinoneoscillationcycleofpointK1

        3.3氣動力響應

        圖16為數值計算得到的豎向自由度頻比響應曲線,圖中橫坐標為無量綱風速,縱坐標為頻率比值響應,即不同風速下模型尾流漩渦脫落或者擺動頻率與模型豎向自由度固有頻率的比值。從圖中可以觀察到兩個明顯的渦振鎖定區(qū)間,在這兩個區(qū)間內,尾流頻率與模型固有頻率基本一致,在整個區(qū)間內鎖定,而在區(qū)間外,尾流頻率隨風速的增加而呈線性增長趨勢,符合Strouhal定律的特征。

        圖13 單個振動周期內測點L2尾流漩渦形態(tài)變化規(guī)律Fig.13ThevorticeschangingrulesinoneoscillationcycleofpointL2圖14 單個振動周期內測點J2尾流漩渦形態(tài)變化規(guī)律Fig.14ThevorticeschangingrulesinoneoscillationcycleofpointJ2圖15 單個振動周期內測點K2尾流漩渦形態(tài)變化規(guī)律Fig.15ThevorticeschangingrulesinoneoscillationcycleofpointK2

        圖16 渦振鎖定區(qū)間頻率響應曲線Fig.16 Frequency response curves in vortex-induced vibration curves

        圖17為數值計算得到的相位差響應,橫坐標同樣為無量綱來流風速,而縱坐標則為氣動升力與豎向自由度位移響應之間的相位差。從圖中可以看出,在兩個分離的渦振鎖定區(qū)間內,隨著渦振區(qū)間的發(fā)展,氣動升力與位移響應之間均存在著相位差,并且不是恒定的,表現出隨著鎖定區(qū)間的發(fā)展而逐步增大的特征,也即說明在同一個渦振鎖定區(qū)間內,氣動升力與位移響應之間的相位差是在不斷增大的。以第一個渦振鎖定區(qū)間為例,圖18對應L1、J1和K1三個特征響應點的氣動升力時程與位移響應時程的對應關系,隨著渦振區(qū)間的發(fā)展,位移響應幅值先增大后減小,氣動升力幅值也是先增大后減小,而氣動升力與位移響應之間的相位差卻是在持續(xù)增大的。L1和J1測點處相位差基本在30°以內,氣動升力總體上在為模型的振動提供能量,而當渦振區(qū)間發(fā)展到K1點,氣動升力與位移響應之間的相位差已經達到120°,此時,氣動升力實際上是在抑制著模型的振動。

        圖17 渦激氣動力與位移響應相位差Fig.17 The phase differences between vortex shedding forces and displacement

        圖18 第一個渦振鎖定區(qū)間內渦激氣動力與位移響應相位差時程對比Fig.18 The time history of phase differences between vortex shedding forces and displacement of the first lock-in

        4結論

        (1) 采用Fluent流固耦合分析實現了模型風洞試驗中觀察到的兩個分離的渦振鎖定區(qū)間現象的數值模擬,在起振風速點、區(qū)間跨度以及渦振幅值關系上,計算值與試驗實測值均吻合良好;

        (2) 數值計算結果顯示,第一個渦振鎖定區(qū)間內的尾流總體呈現出卡門漩渦的形態(tài)特征,而第二個鎖定區(qū)間內特征響應點的尾流中基本沒有卡門漩渦對的出現,而是表現出“魚尾擺動”的尾流形態(tài),擺動幅度越大,模型位移響應幅值越大,并沒有觀察到明顯的卡門渦與模型表面運動渦相互作用的現象;

        (3) 兩個渦振鎖定區(qū)間內,氣動升力與位移響應之間的相位差始終存在,并且隨著渦振區(qū)間的發(fā)展而不斷增大,在渦振區(qū)間的下降段,該相位差甚至達到120°,此時的氣動力實際上是抑制模型振動的。

        [ 1 ] 許福友,丁威, 姜峰,等. 大跨度橋梁渦激振動研究進展與展望[J].振動與沖擊,2010, 29(10): 40-46.

        XU Fu-you, DING Wei, JIANG Feng, et al. Research development and forward of vortex induced vibration on large span bridges [J]. Journal of Vibration and Shock,2010, 29(10): 40-46.

        [ 2 ] 秦浩, 廖海黎, 李明水. 大跨度雙幅連續(xù)鋼箱梁橋渦激振動特性風洞試驗研究[J].振動與沖擊,2014, 33(14): 206-210.

        QIN Hao, LIAO Hai-li, LI Ming-shui. Vortex induced vibration performance of long span continuous steel twin box beam bridge based on wind tunnel test[J]. Journal of Vibration and Shock,2014, 33(14): 206-210.

        [ 3 ] 張志田, 陳政清. 橋梁節(jié)段與實橋渦激共振幅值的換算關系[J].土木工程學報,2011, 44(7): 77-82.

        ZHANG Zhi-tian,CHEN Zheng-qing. Vibration amplitude relationship between practical bridge and section model[J]. Journal of Civil Engineering,2011, 44(7): 77-82.

        [ 4 ] 陳政清. 橋梁風工程 [M].北京: 人民交通出版社, 2005.

        [ 5 ] 陳艾榮. 蘇通橋主橋高雷諾數下渦激振動及三分力特性試驗研究總報告[R]. WT-200418,上海:同濟大學土木工程防災國家重點實驗室,2004.

        [ 6 ] 陳艾榮,馬如進,王達磊,等 基于性能的蘇通大橋抗風設計[J]. 公路, 2009(5): 139-145.

        CHEN Ai-rong, MA Ru-jin, WANG Da-lei, et al. Wind resistance design based on performance of Sutong bridge[J]. Highway,2009(5): 139-145.

        [ 7 ] 陳艾榮,等. 蘇通大橋風荷載研究[R]. 上海: 同濟大學土木工程防災國家重點實驗室,2004.

        [ 8 ] Garrett J L. Flow-induced vibration of elastically supported rectangular cylinders[D]. Iowa: Iowa State University,2003.

        [ 9 ] Matsumoto M. Vortex shedding of bluff bodies: a review [J]. Journal of Fluids and Structures,1999, 13(1):791-811.

        [10] Tamura Y, Matsui G. Wake-oscillator model of vortex-induced oscillation of circularcylinder[C]//Proceedings of the 5th international conference on wind engineering, Fort Collins, USA, July,1979.

        [11] Chen Zheng-qing, Chen Wen, Hua Xu-gang, et al. Higher modes’ vertical vortex-induced vibrations of suspension bridges—Part 2: aeroelastic model study[C]// University of Nottingham. EACWE VII, Cambridge: University of Nottingham,2013.

        [12] Eurocode: Actions on structures[S]. Brussels: CEN national Members, 2004.

        [13] Sarpkaya T. A critical review of the intrinsic nature of vortex-induced vibrations [J]. Journal of Fluids and Structures, 2004,19(1): 389-447.

        [14] Bearman P W, Gartshore I S, Maull D J, et al.Experiments on flow-induced vibration of a square section cylinder[J].J. fluids struct., 1987,1(1):19-34.

        [15] Simiu E, Scanlan R H. Wind effects on structures-an introduction to wind engineering[M]. ThirdEdition,1996.

        Numerical simulation for two separate vortex-induced vibration lock-in intervals of bridge sections

        HUANG Li1, ZHOU Shuai2, LIANG Peng3

        (1. College of Water Conservancy and Civil Engineering, South China Agricultural University, Guangzhou 510642, China;2. Wind Engineering Research Center, Hunan University, Changsha 410082, China;3. The First Engineering Company of CCCC Fourth Harbor Engineering Co. Ltd., Guangzhou 510500, China)

        For some bridge sections or similar bluff bodies with a large aspect ratio, two separate vortex-induced vibration (VIV) lock-in intervals of the same model shape can be observed in the same section model wind tunnel tests, and the dominant oscillation frequencies of the two lock-in intervals are the same, these do not agree with the traditional Strouhal law. In order to get a further understanding aerodynamic mechanism of these oscillations, a rectangular section with an aspect ratio of 6 was taken as a study object, and based on its 2D section model wind tunnel tests parameters and results, Fluent-based numerical simulations were performed. The simulation results agree well with test ones in both lock-in ranges and amplitude responses. Then, more details about aerodynamic forces and wake vortices were obtained using the post processing of Fluent. The study showed that the wake vortices of the first VIV lock-in interval is the typical Karmen vortex mode, while the second one is full different from the first one, its wake vortices are observed to be “fish tail waving” mode; in the two separate lock-in intervals, there is always a phase difference between aerodynamic lift and vertical displacement response, and the phase difference continuously increases with developing of the lock-in interval; moreover, the jump of the phase difference of the first lock-in interval is bigger than that of the second one.

        bridges; vortex-induced vibration; lock-in interval; numerical simulation; wake vortices

        10.13465/j.cnki.jvs.2016.11.008

        國家自然科學基金(91215302;51278189);廣東省自然科學基金(2015A030310164)

        2015-09-01修改稿收到日期:2015-12-04

        黃俐 女,博士,講師,1984年生

        U446.1

        A

        猜你喜歡
        渦振尾流風洞試驗
        分體式雙箱梁渦振氣動控制措施數值模擬
        結構工程師(2020年4期)2020-11-12 03:00:04
        板桁結合梁渦振性能及抑振措施研究
        大跨度懸索橋渦振風洞試驗與現場實測比較
        結構工程師(2018年3期)2018-07-14 09:18:44
        低風壓架空導線的風洞試驗
        電線電纜(2017年5期)2017-10-18 00:52:03
        飛機尾流的散射特性與探測技術綜述
        雷達學報(2017年6期)2017-03-26 07:53:06
        滾轉機動載荷減緩風洞試驗
        錐形流量計尾流流場分析
        力-狀態(tài)映射法在橋梁斷面渦振研究中的應用
        水面艦船風尾流效應減弱的模擬研究
        遮擋條件下超高層建筑風洞試驗研究
        重慶建筑(2014年12期)2014-07-24 14:00:32
        精品的一区二区三区| 免费人成网ww555kkk在线| 欧美巨大性爽| 久久一日本道色综合久久大香| 一区二区三区av资源网| 久久免费看黄a级毛片| 狠狠做深爱婷婷久久综合一区| 欧美黑人性色黄在线视频| 亚洲二区三区四区太九| 日韩av毛片在线观看| a级毛片100部免费观看| 国产精品美女| 中文字幕亚洲乱码熟女在线| 天天射综合网天天插天天干| 亚洲国产日韩精品一区二区三区 | 五十六十日本老熟妇乱| 精品久久人人爽天天玩人人妻| 蜜桃一区二区三区在线看| 中文字幕东京热一区二区人妻少妇| 精品av熟女一区二区偷窥海滩| 久久国产精品二国产精品| AV在线中出| 羞羞色院99精品全部免| 爆乳熟妇一区二区三区霸乳 | 国产av无毛无遮挡网站| 久久97久久97精品免视看| 精品国产a∨无码一区二区三区| 亚洲一区二区成人在线视频| 亚洲av高清天堂网站在线观看| 亚洲av日韩精品久久久久久久| 手机看片福利日韩| 国产av一区二区内射| 每日更新在线观看av| 黄色视频在线免费观看| 国产AV无码无遮挡毛片| av黄色在线免费观看 | 荡女精品导航| 国产精品国产三级在线专区| 人妻熟妇乱又伦精品视频| 亚洲av无码之日韩精品| 亚洲精品白浆高清久久|