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

        ?

        基于特征系統(tǒng)實現(xiàn)算法的輸電塔氣動阻尼風洞試驗研究

        2014-09-18 09:55:58段成蔭鄧洪洲
        振動與沖擊 2014年21期
        關鍵詞:角下阻尼比風向

        段成蔭,鄧洪洲

        (同濟大學 建筑工程系,上海 200092)

        輸電塔是一種由角鋼或鋼管通過螺栓連接或焊接組成的格構式結構,相較于較小的結構阻尼,氣動阻尼對輸電塔結構風振響應具有較大的影響。通常的觀點是結構氣動阻尼與風速之間存在較強的依賴性,根據(jù)準定常理論順風向氣動阻尼應與來流平均風速(或折減風速)成正比[1],近年來一些氣動阻尼風洞試驗研究正是基于這一關系同時考慮其他可能的影響因素如風場類型[2]、結構阻尼[3]、結構外形[4-6]等開展了一些工作,但這些研究的對象是高層建筑或者大跨屋蓋這些鈍體外形的結構。格構式輸電塔迎風面密實度較低,桿件之間存在互相影響,氣動阻尼比較離散[7]。雖然順風向氣動阻尼隨風速有增大趨勢[8],但與準定常假定是否相符仍需進一步的驗證,根據(jù)Takeuchi等[7]的實測研究,氣動阻尼并非隨風速增大而成比例地增大。而且由于輸電塔氣動外形與傳統(tǒng)的鈍體相去甚遠,不太可能發(fā)生整體渦振現(xiàn)象,其橫風向氣動阻尼可能與高層建筑出現(xiàn)負值[9]的規(guī)律有所不同。

        特征系統(tǒng)實現(xiàn)算法(Eigensystem Realization Algorithm,ERA)[10]是一種基于狀態(tài)空間的時域模態(tài)參數(shù)識別方法,利用系統(tǒng)脈沖響應函數(shù)或者相關函數(shù)得到系統(tǒng)的最小實現(xiàn)。ERA識別模態(tài)參數(shù)僅需很短的響應歷程,識別速度快,對密頻模態(tài)的識別效果好,得到了很多研究者的重視并且發(fā)展出了諸如ERA/DC(Eigensystem Realization Algorithm using Data Correlation)[11]、FERA(Fast ERA)[12]等改進算法。然而ERA也存在抗噪能力差的缺點,往往會出現(xiàn)大量的虛假模態(tài),模型定階和系統(tǒng)真實模態(tài)的提取目前尚未形成完善的理論,識別過程中仍存在一定主觀因素的影響。

        本文首先對特征系統(tǒng)實現(xiàn)算法中虛假模態(tài)的剔除問題進行了研究,提出模態(tài)參與因子(Modal Participation Factor,MPF)作為識別真實模態(tài)與虛假模態(tài)的一個指標,并給出了估算模態(tài)能量貢獻(Modal Energy Contribution,MEC)的簡化計算公式,基于MPF和MEC采用自動識別程序對一輸電塔氣彈模型環(huán)境激勵下的模態(tài)參數(shù)進行識別,得到了不同風向角下模型兩個主軸方向的氣動阻尼,并研究了其與來流平均風速的關系。

        1 模態(tài)參數(shù)識別與提取

        1.1 特征系統(tǒng)實現(xiàn)算法

        由Juang等人[10]提出的特征系統(tǒng)實現(xiàn)算法(ERA)基于系統(tǒng)的離散時間域狀態(tài)空間模型,利用脈沖響應序列即Markov參數(shù)的特點,構造Hankel矩陣:

        式中,S 為對角陣,U、V 為酉陣,SN、UN、VN為最大的前N 個奇異值所對應的子矩陣,Em=[Im,O,…,O]T,Ep=[Ip,O,…,O]T,m 為外激勵數(shù),p輸出通道數(shù),I為單位矩陣,O為零矩陣。

        通過對A1進行特征值分解得到以共軛對的形式出現(xiàn)的特征值 λ1,λ2,…,λn,…。

        從而模態(tài)的頻率和阻尼比可以通過如下方法計算,Δt為采樣間隔:

        總阻尼減去結構阻尼即為氣動阻尼:

        為了提高ERA識別的抗噪能力,本文采用脈沖響應的相關系數(shù)r(k)來代替脈沖響應v(k)構造Hankel矩陣,能有效抑制噪聲影響[15]:

        式中,Nr為相關數(shù)據(jù)長度。秦仙蓉等[13]指出在10%和30%的噪聲下ERA/DC算法比ERA精度高6~60倍。

        1.2 隨機減量技術——RDT

        ERA算法的關鍵就是根據(jù)系統(tǒng)的脈沖響應序列集成Hankel矩陣,脈沖響應可基于實測信號(位移、速度、加速度等)采用隨機減量技術(Random Decrement Technique,RDT)[14]得到。隨機減量技術利用隨機響應中由系統(tǒng)固有特性決定的確定性成份和由環(huán)境激勵決定的隨機成份二者在統(tǒng)計上的不同性質,通過指定某種合理條件,然后采用平均的方法將隨機成份濾掉從而分離出自由衰減信號。

        如果選取適當?shù)拈撝祒=x0去截取該樣本曲線x(t),交點對應的時刻ti(i=1,2,…),可得到一系列從ti時刻開始的子樣本曲線xi(t),每條子樣本曲線可看成由ti時刻初位移、初速度引起的自由振動和從ti時刻開始隨機激勵引起的強迫振動三部分組成,RDT通過對各條子樣本曲線xi(t)作算數(shù)平均得到自由衰減信號。

        1.3 模態(tài)識別指標

        受各種噪聲影響,環(huán)境激勵下結構響應通過ERA計算出來的模態(tài)成份是十分豐富的,如何提取系統(tǒng)的真實模態(tài)是ERA應用中的一個問題,通常提取較大奇異值所對應的模態(tài),但這種做法仍無法有效剔除虛假模態(tài),特別是無法對結果進行歸類和自動提取。

        文獻[15]認為各階傳遞函數(shù)的全頻域積分可作為衡量各階模態(tài)能量貢獻MEC的指標:

        式中,Ф為A1的特征向量矩陣,fs為采樣頻率,Δt為采樣間隔,"(·)表示求最大奇異值。式(6)雖然可以反映模態(tài)的能量水平,但需要進行全頻域的積分運算,運算量較大,而且如果虛假模態(tài)阻尼比很小,采用式(6)仍難以判斷。

        系統(tǒng)的脈沖響應可分解為各階模態(tài)之和:

        經推導各階信號vj經可寫為:

        可見,v(k)為代表各階模態(tài)的脈沖信號之和,MPF恰為各階信號的初始振幅,根據(jù)MPF的大小可分離各階模態(tài),理論上噪聲模態(tài)MPF應接近于0,因此前若干階MPF對應的模態(tài)即為系統(tǒng)真實模態(tài)。

        各階模態(tài)的能量貢獻式(6)可以用更為簡單的估算公式表示為:

        Δωe,j= πξjωj為傳遞函數(shù)的白噪聲帶寬,由式(9)和式(10)可見,MEC與識別的頻率和阻尼比有很大關系,當阻尼比很小時(通常為噪聲模態(tài))MEC將很大,而MPF則是與頻率和阻尼比無關的振幅,因而可以避免小阻尼虛假模態(tài)共振峰的干擾。

        1.4 模態(tài)參數(shù)提取步驟

        本文對于ERA識別結果的提取過程如下:

        (1)計算各識別結果的MPFj和MECj;

        (2)若MPFj>閾值則保留,否則判定為虛假模態(tài),為不至于發(fā)生遺漏,閾值的選取應使保留模態(tài)數(shù)量大于系統(tǒng)階次;

        (3)根據(jù)MPF大小對保留模態(tài)進行排序,若某兩個結果的MPF比較接近則根據(jù)其MEC排序,即可區(qū)分不同階次模態(tài)參數(shù);

        (4)最后,將不同Hankel陣階次下的提取結果形成穩(wěn)定圖,確定是否為有效識別結果,并取平均值輸出為模態(tài)參數(shù)提取值。

        該程序基本可以由計算機自動完成,能提高識別效率以及減小人工干預的主觀影響。

        表1 模型主要相似比Tab.1 Main scales for the model

        2 風洞試驗概況

        風洞試驗的原型為一雙回路500 kV輸電塔,總高度101 m。氣彈模型的縮尺比和風速比根據(jù)風洞試驗段的大小分別定為λL=1/80和λL=1/3(模型/原型,下同),模型主要相似比見表1。

        如圖1所示,定義水平橫擔的方向為X向,垂直于橫擔方向為Y向,定義風向角為來流方向與Y軸的夾角,加速度傳感器布置在模型下橫擔處X向和Y向。

        圖1 測點布置和坐標定義Fig.1 Model instrumentation and coordinate definition

        圖2 風剖面和紊流度剖面Fig.2 Simulated mean velocity and turbulence intensity profiles

        通過人工激振法得到的模型動力特性如表2所示,表中理論值表示ANSYS模態(tài)分析結果。

        表2 模型動力特性Tab.2 Eigenvalues of the model

        試驗在同濟大學TJ-3低速邊界層風洞進行,試驗風剖面指數(shù)接近于α=0.16,紊流度接近于30 m(模型0.375 m)高紊流度Iu=16%,試驗風剖面和紊流度剖面如圖2所示。

        試驗分別在 0°、15°、30°、45°、60°、75°、90°七個風向角下進行,風速范圍為4~11 m/s,采樣頻率250 Hz,采樣時長69.632 s,樣本容量17 408,控制風速的參考高度ZR位于模型頂部高度1.263 m。

        3 結果分析與討論

        3.1 模態(tài)提取

        以0°風向6.24 m/s風速下Y向加速度識別為例,其功率譜如圖3(a)所示,采用60階Hankel矩陣時各模態(tài)的MEC譜和MPF譜見圖3(b)和(c)。最大MEC(0.034)對應的頻率為97 Hz,甚至大于32 Hz處真實模態(tài)的值(0.022),但在功率譜圖3(a)中97 Hz附近并無顯著共振峰,因此為噪聲模態(tài),結合式(10)不難看出其MEC過大的原因就在于阻尼比識別值4.6×10-6很小。圖3(c)中97 Hz處識別結果的MPF僅為0.011,表明該階信號振動幅度很小,確為虛假模態(tài),32 Hz處MPF(0.33)明顯大于其他值,故判定其代表一階真實模態(tài),至于高階模態(tài)是否存在從圖3(c)的單一結果中無法確定,因此保守地提取MPF最大的前三個結果,并將不同Hankel陣階次下的提取值形成穩(wěn)定圖。

        圖3 奇異值譜(a)、MEC譜(b)和MPF譜(c)Fig.3 Spectra of singular value(a),MEC(b)and MPF (c)

        圖4 剔除前穩(wěn)定圖(a)和按本文方法提取的前三階頻率(b)、阻尼比(c)穩(wěn)定圖Fig.4 Stabilization diagrams of original results(a),the first three identified frequencies(b)and damping ratios(c)

        圖4 為模態(tài)參數(shù)識別結果的穩(wěn)定圖。圖4(a)為不經任何處理的ERA結果,頻率成份十分豐富,難以直接從穩(wěn)定圖中辨別真實模態(tài)。圖4(b)和圖4(c)為按本文識別方法提取最大三個MPF對應的頻率和阻尼比穩(wěn)定圖,一階顯然是真實模態(tài),頻率十分穩(wěn)定,阻尼比偏差稍大但也基本保持常數(shù),二、三階則并不穩(wěn)定,只能準確識別一階物理模態(tài)。原因一方面是對于輸電塔這類高聳結構,其一階彎曲模態(tài)通常占有統(tǒng)治地位,高階模態(tài)的參與程度很小,很可能被噪聲所掩蓋;另一方面試驗時為了避免加速度傳感器和數(shù)據(jù)線對模型造成太大影響,測點布置在塔身中部偏上,該處與二階彎曲振型的節(jié)點接近,難以識別二階模態(tài)??偠灾?,本文提出的判別指標MPF對真實模態(tài)進行提取的方法是行之有效的,提取過程中人工干預很少,大大減小了識別過程中的主觀因素。

        3.2 氣動阻尼結果與討論

        對輸電塔模型環(huán)境激勵下0°~90°風向角、4~11 m/s平均風速下X向和Y向各工況的1階頻率和阻尼比的識別結果如圖5~6圖所示。

        X向和Y向一階頻率如圖5(a)和圖6(a)所示,兩個方向頻率基本相當,各工況處于31.2~33.2 Hz之間,這與表2中采用人工激振法識別的頻率基本一致。頻率與平均風速有較強的相關性,但與風向角沒有明顯的關系,頻率隨著風速的增大而降低,各風向角下均呈現(xiàn)相同趨勢,體現(xiàn)空氣附加質量和氣動剛度對輸電塔振動頻率具有一定影響。

        圖5(b)和圖6(b)為X向和Y向的一階阻尼比識別結果。絕大多數(shù)工況的阻尼比大于表2給出的1.4%,氣動阻力大多以正阻尼的形式存在,這對結構是有利的。但各工況阻尼比的離散性較大,與風速也并非線性關系,在不同風向角下的變化曲線也各不相同。但是圖5(b)的 X 向阻尼比在60°、75°、90°風向角下隨著風速增大呈上升的趨勢,而 0°、15°、30°風向角下基本在一個定值附近波動。圖6(b)的Y向阻尼比則相反,0°、15°、30°風向角下隨風速增大而增大,60°、75°、90°風向角下基本保持定值。注意到對X向90°角是順風向而對Y向0°角是順風向,兩個方向以0°和90°為代表的縱橫向氣動阻尼與風速的關系如圖5(c)和圖6(c)所示??梢妼τ陧橈L向氣動阻尼,輸電塔結構仍與高層建筑隨風速(折減風速)增大而單調遞增[16]具有相同的變化趨勢;但在橫風向,識別結果中并未觀察到明顯的氣動負阻尼,而是基本在0或者較小的水平周圍,這是與高層建筑[9,16]不同的特點,原因是輸電塔由于格構式的結構特點難以發(fā)生整體渦振,橫風向仍以結構阻尼為主。

        圖5 模型識別結果Fig.5 Model identification results

        圖6 模型識別結果Fig.6 Model identification results

        圖7 氣動阻尼與相應風速分量的關系Fig.7 Model’s aerodynamic damping versus corresponding wind speed components

        如前文所述,在不同風向角下氣動阻尼與風速的關系各異,很難用一個統(tǒng)一的公式進行描述,但X和Y向氣動阻尼在相近角度呈現(xiàn)出相似的規(guī)律,表明其與風向角存在某種內在聯(lián)系,一個最直接的考慮就是X、Y兩主軸的氣動阻尼與相應方向的平均風速分量的關系。

        從圖7可以看出,X和Y向氣動阻尼與相應風速分量呈正相關關系,相關系數(shù)分別為0.623 8和0.774 0,線性相關性仍然不強。圖中同時給出了根據(jù)準定常理論[1]的氣動阻尼取值:

        式中,ρ為空氣密度;f為一階模態(tài)頻率;m(z)為z高度處質量;面積A(z)分別采用正側面投影面積;平均風速(z)取試驗風剖面;阻力系數(shù)CD根據(jù)司瑞娟[20]的研究取1.5;振型φ采用理論分析的一階振型系數(shù)。準定常理論大體能反映試驗氣動阻尼的變化趨勢,取值也跟多數(shù)試驗結果相當,但在較高風速時有所低估。Gabbai等[18]提出氣動阻力可能與相對風速的N次方(而非準定常假定的平方)成正比,這個關系雖然物理意義不如準定常假定明確,但從圖7來看,輸電塔結構氣動阻尼的確表現(xiàn)出了一些高階效應,因此這里在準定常理論的基礎上增加平均風速的高階項來表示氣動阻尼的趨勢:

        式中,a,b,c為擬合參數(shù),見表3所示。兩個方向通過式(12)計算的結果亦繪于圖7中,該取值比準定常理論更為準確地表現(xiàn)了氣動阻尼與風速分量的變化趨勢。

        表3 擬合參數(shù)和標準誤差Tab.3 Fitted parameters and standard error

        4 結論

        (1)提出了一個用以提取ERA識別結果中真實模態(tài)的參數(shù)MPF,具有明確的物理意義:從信號的角度講,MPF代表組成脈沖響應的各階調幅信號的初始振幅;從系統(tǒng)的角度來講,MPF代表系統(tǒng)各階模態(tài)的參與程度。

        (2)給出了計算模態(tài)能量貢獻MEC的簡化公式,根據(jù)MPF和MEC可以有效地剔除虛假模態(tài)并自動提取真實模態(tài),降低了識別過程中的主觀因素。

        (3)環(huán)境激勵下輸電塔一階頻率隨風速的增加而略有減小,空氣附加質量和氣動剛度對輸電塔結構具有一定影響。

        (4)輸電塔氣動阻尼大多為正,順風向隨風速的增大呈增大趨勢,橫風向未出現(xiàn)明顯的負阻尼,而是在較小的取值周圍波動。

        (5)角度風下X、Y方向氣動阻尼與相應方向的平均風速分量正相關,但線性關系不強,準定常理論能反映氣動阻尼的大體趨勢,在較高風速時有所低估。在準定常理論的基礎上增加高階項能更準確地表現(xiàn)試驗結果的變化。

        [1]Holmes J D.Along-wind response of lattice towers.Ⅱ:Aerodynamic damping and deflections[J].Engineering Structures,1996,18(7):483-488.

        [2]曹會蘭,全涌,顧明.風場類型對方形超高層建筑順風向氣動阻尼的影響研究[J].振動與沖擊,2012,31(16):22-26.CAO Hui-lan,QUAN Yong,GU Ming.Effect of roughness exposure on aerodynamic damping of square high-rise buildings[J].Journal of Vibration and Shock,2012,31(16):22-26.

        [3]曹會蘭,全涌,顧明,等.獨立矩形截面超高層建筑的順風向氣動阻尼風洞試驗研究[J].振動與沖擊,2012,31(5):122-127.CAO Hui-lan,QUAN Yong,GU Ming,et al.Effect of roughness exposure on aerodynamic damping of square highrise buildings[J].Journal of Vibration and Shock,2012,31(5):122-127.

        [4]曹會蘭,全涌,顧明.一類準方形截面超高層建筑順風向氣動阻尼[J].振動與沖擊,2012,31(22):84-89.CAO Hui-lan, QUAN Yong, GU Ming. Along-wind aerodynamic damping of high-rise buildings with aerodynamically modified square cross-sections[J].Journal of Vibration and Shock,2012,31(22):84-89.

        [5]樓文娟,盧旦,楊毅,等.開孔建筑屋蓋風振響應中的氣動阻尼識別[J].空氣動力學報,2007,25(4):419-424.LOU Wen-juan,LU Dan,YANG Yi,et al.Identification of aerodynamic damping for roof wind-induced response of opening building[J].Acta Aerodynamica Sinica,2007,25(4):419-424.

        [6]潘峰,孫炳楠,樓文娟.基于Hilbert-Huang變換的大跨屋蓋氣動阻尼識別[J].浙江大學學報(工學版),2007,41(1):65-70.PAN Feng,SUN Bing-nan,LOU Wen-juan.Aerodynamic damping identification of long-span roof based on Hilbert-Huang transform [J]. Journal of Zhejiang University(Engineering Science),2007,41(1):65-70.

        [7]Takeuchi M,Maeda J,Ishida N.Aerodynamic damping properties of two transmission towers estimated by combining several identification methods[J]. Journal of Wind Engineering and Industrial Aerodynamics,2010,98(12):872-880.

        [8]任坤,李正良,肖正直,等.環(huán)境激勵下特高壓輸電塔線體系氣動阻尼的識別[J].重慶工學院學報(自然科學版),2009,23(7):64-68.REN Kun, LI Zheng-liang, XIAO Zheng-zhi, et al.Aerodynamic damping identification of UHV transmission line system under ambient excitation[J].Journal of Chongqing University of Technology(Natural Science),2009,23(7):64-68.

        [9]吳海洋,梁樞果,陳政清,等.強風下方截面高層建筑橫風向氣動阻尼比研究[J].工程力學,2010,27(10):96-103.WU Hai-yang,LIANG Shu-guo,CHEN Zheng-qing,et al.Research on the aerodynamic damping ratios of square tall buildings in across-wind direction under strong wind[J].Engineering Mechanics,2010,27(10):96-103.

        [10] Juang J N,Pappa R S.An eigensystem realization algorithm for modal parameter identification and model reduction[J].Journal of Guidance,Control and Dynamics,1985,8(5):620-627.

        [11] Juang JN,Cooper JE.An einginsystem realization algorithm using data correlation(ERA/DC) for model parameter identification[J].Control Theory and Advanced Technology,1988,4(1):5-14.

        [12]劉福強,張令彌.一種改進的特征系統(tǒng)實現(xiàn)算法及在智能結構中的應用[J].振動工程學報,1999,12(3):316-322.LIU Fu-qiang,ZHANG Ling-mi.An improved eigensystem realization algorithm and its application to modal parameter identification of intelligent space trusses[J].Journal of Vibration Engineering,1999,12(3):316-322.

        [13]秦仙蓉,王彤,張令彌.模態(tài)參數(shù)識別的特征系統(tǒng)實現(xiàn)算法:研究與比較[J].航空學報,2001,22(4):340-342.QIN Xian-rong,WANG Tong,ZHANG Ling-mi.Comparison of different variants of the eigensystem realization algorithm in modal parameter identification[J].Acta Aeronautica Et Astronautica Sinica,2001,22(4):340-342.

        [14] Cole H A J.Methods and apparatus for measuring the damping characteristics of a structure: United States,3620069[P].1971-11-16.

        [15]章國穩(wěn),湯寶平,潘飛.特征系統(tǒng)實現(xiàn)算法的虛假模態(tài)剔除方法[J].重慶大學學報,2012,35(3):20-25.ZHANG Guo-wen,TANG Bao-ping,PAN Fei.The method for removing the spurious modes of eigensystem realization algorithm[J].Journal of Chongqing University,2012,35(3):20-25.

        [16]全涌,顧明.方形斷面高層建筑的氣動阻尼研究[J].工程力學,2004,21(1):26-30.QUAN Yong, GU Ming. Wind tunnel test study of aerodynamic damping of super highrise buildings[J].Engineering Mechanics,2004,21(1):26-30.

        [17]司瑞娟.特高壓輸電塔線體系風振響應及抗風設計參數(shù)研究[D].上海:同濟大學,2011.

        [18] Gabbai R D,Simiu E.Aerodynamic damping in the alongwind response of tall buildings[J].Journal of Structural Engineering,2010,136(1):117-119.

        猜你喜歡
        角下阻尼比風向
        建筑物對塔機順風向風力干擾效應研究
        基于細觀結構的原狀黃土動彈性模量和阻尼比試驗研究
        地震研究(2021年1期)2021-04-13 01:05:24
        兇手是A角
        黏滯阻尼器在時程分析下的附加有效阻尼比研究
        振動與沖擊(2019年4期)2019-02-22 02:33:34
        波形分析法求解公路橋梁阻尼比的探討
        上海公路(2018年3期)2018-03-21 05:55:48
        結構構件阻尼比對大跨度懸索橋地震響應的影響
        自然與風Feeling Nature
        行業(yè)統(tǒng)計帶來哪些風向?
        風向
        風能(2015年8期)2015-02-27 10:15:11
        風向
        風能(2015年4期)2015-02-27 10:14:30
        亚洲天堂av路线一免费观看| 美丽人妻被按摩中出中文字幕| 久久精品无码一区二区2020| 久久无人码人妻一区二区三区 | 国产亚洲精品熟女国产成人| 熟妇激情内射com| 国产亚洲欧美日韩综合一区在线观看 | 婷婷午夜天| 亚洲激情人体艺术视频| 免费看黄片视频在线观看| 97久久婷婷五月综合色d啪蜜芽| 日本少妇被黑人xxxxx| 巨臀中文字幕一区二区| 日韩精品视频av在线观看| 精品久久久久久亚洲综合网| 中文字幕av一区中文字幕天堂| 精品少妇一区一区三区| 国产亚洲av夜间福利在线观看| 久久天天躁夜夜躁狠狠| 无码国产激情在线观看| 水蜜桃一二二视频在线观看免费 | av大全亚洲一区二区三区| 无遮挡又爽又刺激的视频| 国产精品原创av片国产日韩| 日本免费精品免费视频| 天堂а在线中文在线新版| 伊伊人成亚洲综合人网7777| 国产在线观看www污污污| 欧洲一区在线观看| 国产爽快片一区二区三区| 东京热久久综合久久88| 日日摸夜夜添夜夜添无码免费视频 | 国产精品一区二区久久不卡| 亚洲熟妇色xxxxx欧美老妇 | 午夜免费视频| 国产精品麻豆最新AV| 99久久免费精品色老| 亚洲精品乱码久久久久久| 国产精品麻豆成人av电影艾秋 | 日本在线免费一区二区三区| 国产性自爱拍偷在在线播放|