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

基于特征波導納法和敏感性分析的大型周期支撐結構損傷識別

2017-11-30 05:49:39尹孟林
振動與沖擊 2017年22期
關鍵詞:模態結構模型

尹 濤, 尹孟林

(武漢大學 土木建筑工程學院, 武漢 430072)

基于特征波導納法和敏感性分析的大型周期支撐結構損傷識別

尹 濤, 尹孟林

(武漢大學 土木建筑工程學院, 武漢 430072)

提出了一種基于特征波導納法和敏感性分析的大型周期支撐結構損傷識別方法。通過分析存在單一損傷單元有限周期支撐結構的自由波傳播規律,建立了周期結構無量綱自振頻率變化對基本周期單元剛度變化的敏感性分析矩陣;求解敏感性識別方程組,實現基于測量自振頻率變化的大型周期支撐結構損傷檢測;同時,分析表明所提出無量綱自振頻率的敏感性與結構參數無關,即無需知道原始結構的準確幾何物理參數就能確定敏感性系數。分別以一周期支撐梁與周期支撐法蘭接頭管道為例進行研究,表明該方法僅用量測損傷前后的前幾階自振頻率變化就能較準確地進行周期支撐結構多損傷識別。

特征波導納; 敏感性分析; 周期結構; 損傷識別; 法蘭連接管道

基于振動的損傷檢測方法主要分為基于動力學模型的有模型方法和基于信號分析的無模型方法。其中,無模型方法通常僅能判斷損傷發生的部位,而有模型動力檢測法則可同時實現損傷部位的判定以及損傷程度的評估,基于敏感性分析的損傷檢測方法就是其中重要一類,并受到廣泛關注[1]。結構損傷識別中的敏感性分析方法是基于固有模態參數對結構參數變化敏感度的分析計算,判斷結構是否出現損傷并確定損傷位置及其程度。該方法的關鍵在于選取可測且對結構損傷敏感的參數,報道較多的損傷檢測參數有結構的固有頻率、模態振型以及模態振型曲率等。如Cawley等[2]最早將固有頻率敏感性分析運用于結構損傷識別研究,通過固有頻率敏感性分析,采用一組非完備的固有頻率測量值對簡單結構的損傷位置和損傷程度進行了識別。Biswas等[3]利用振型和振型曲率的差異進行損傷定位,取得了一定的效果;Wahab等[4]進一步分析了利用振型進行結構損傷定位的缺陷,提出了基于振型曲率的敏感性參數。薛松濤等[5]推導了結構損傷前后頻率變化與結構剛度變化之間的二階敏感性公式,利用多層框架結構在多種損傷下的振動試驗數據確定其損傷位置和程度;Yin等[6]結合自振頻率、模態振型及模態應變能的敏感性分析方法與攝動法開展了框架結構的概率損傷識別研究。

上述研究均建立在依據健康結構模型參數獲得的模態參數敏感性系數基礎上,而對于那些無法準確獲知結構幾何物理參數的情況,以上方法就存在著較大局限性。此外,實際工程結構如高層建筑、多跨橋梁、長輸油氣管線以及鐵路軌道系統等都可視作由相同子結構首尾串聯而成的鏈狀周期結構系統。周期結構動力學分析的顯著優勢在于,對其自由振動進行精確的波動分析不要求對結構進行完全的模擬[7],利用周期結構這一特性,可以簡化敏感性分析的過程,提高計算效率。朱宏平等[8-9]運用波傳播理論得到了周期彈簧-質量系統的頻率特征方程,并通過模態參數敏感性分析對多層建筑模型進行損傷識別,但該方法應用范圍僅限于能簡化為彈簧-質量系統的結構。

本文在朱宏平等的研究基礎上,運用波傳播理論分析了具有N個基本周期單元的大型周期支撐結構自由振動問題,提出一種無量綱自振頻率概念,并得到無量綱自振頻率與基本周期單元整體剛度變化率之間的關系。通過建立無量綱自振頻率變化率對單元損傷的敏感性識別方程組,并采用約束優化方法求解該識別方程組,以實現周期支撐結構的損傷識別。此外, 本文周期支撐結構無量綱自振頻率的敏感性不依賴于結構的具體幾何物理參數(如質量密度、彈性模量、結構長度及橫截面形狀、尺寸等),具有顯著的優越性和實用性。分別通過對一周期支撐梁模型(考慮單元整體剛度降低)和周期支撐法蘭管道(僅單元局部剛度降低與附加質量影響)模型的損傷識別數值仿真,對本文方法的正確性和有效性進行驗證。

1 理論背景

1.1 單擾亂周期支撐結構頻率特征方程建立

圖1為含N個基本周期單元的單耦合周期結構,假定其第j個單元為擾亂(或損傷)單元,且記兩任意邊界分別為C和D。若在結構左端C處施加激勵力Fc,則振動波將從點C出發在該結構中傳播。當自由波向前傳播到達擾亂單元時,將分為兩部分:一部分為向激勵點反射的反射波,另一部分為越過擾亂單元繼續前進的傳遞波,且傳遞波繼續在擾亂單元右側傳播并在右邊界D處產生反射。因此,結構中任意點波動可表示為相應的傳遞波與反射波之和。設兩者在左端C與右端D處的位移分別為XCt、XCr與XDt、XDr,則擾亂單元左右兩端A與B處的位移分別表示為

XA-=XA-t+XA-r=XCte-(j-1)μ+XCre(j-1)μ=
αwtFAt+αwrFAr
XB+=XB+t+XB+r=XDte(N-j)μ+XDre-(N-j)μ=
αwtFBt+αwrFBr

(1)

式中:XA-t與XA-r分別為傳遞波和反射波在節點A處的位移;XB+t與XB+r分別為傳遞波和反射波在節點B處的位移;FAt、FAr與FBt、FBr分別為傳遞波與反射波在點A與點B處的廣義力,μ為自由波傳播常數。αwt和αwr為特征波導納,其表達式為[10]

αwt=αll-αlre-μ
αwr=αll-αlreμ

(2)

式中,αll和αlr分別為健康周期單元的直接導納和間接導納。

圖1 具有單擾亂的單耦合周期系統波傳播示意圖Fig.1 Wave propagation of mono-coupled periodic structure system with single disorder

擾亂單元j兩端點A、點B處的廣義力FA、FB可表示為

FA=FAt+FAr
FB=FBt+FBr

(3)

則擾亂單元兩端力與位移之間可通過單元導納聯系,即

XA=αAAFA-αABFB
XB=αBAFA-αBBFB

(4)

式中,αAA,αBB和αAB,αBA分別為擾亂單元的直接導納和間接導納。

設周期結構右端部D的導納為αD, 則其位移XD可用力FD和αD表示為

XD=XNt+XNr=αDFD

(5)

依據擾亂單元j兩端點A和點B的位移協調與力平衡條件,利用式(1)~式(5)得到周期結構左端點C處反射波和傳遞波位移X0r和X0t之比為

(6)

其中,

(7)

同時,

(8)

進一步可得單耦合周期結構中分段表示的任意節點i位移為

(9)

其中,

g=αwr(αAA-αwt)e-(j-1)μ+αwt(αAA-αwr)e(j-1)μΦ

(10)

另,健康狀態周期結構需滿足以下關系

αAA=αll,αBB=αrr,αAB=αBA=αlr

(11)

(12)

特別地,在式(9)中,若分母為0,則任意節點處的位移Xi為無限大,即激勵頻率與結構自振頻率相等,產生共振。因此,單擾亂單元周期結構的自振頻率可通過式(13)求解

1+Φ=0

(13)

理論上,式(13)中給出的頻率特征方程適用于任何存在單擾亂的單耦合周期結構體系,但文獻[8]中基于該理論的損傷識別方法僅適用于簡單的彈簧-質量系統,應用范圍受到較大影響,本文將該理論推廣應用于解決更具一般性的大型周期支撐結構損傷識別問題。圖2為本文所研究的具有N個基本周期單元的周期支撐結構,其中,將整體結構左端邊界C固定以消除結構縱向對稱性對損傷識別結果的影響,并假定第j單元為損傷單元,損傷程度以單元整體剛度降低來表征。

圖2 左端固定的周期支撐結構Fig.2 Periodically-supported structure with left-end fixed

圖2所示結構右端邊界D的導納為αD=∞,將其代入式(8)中可得

(14)

周期支撐結構的健康周期單元的直接導納與間接導納可以分別表示為

(15)

式中:ks=ω1/2(mb/D)1/4為周期結構波數;D=EI為結構抗彎剛度;E為彈性模量;I為截面慣性矩;ω為自振圓頻率;mb為單位長度質量;L為基本周期單元長度。相應地,損傷單元j的單元導納為

(16)

令無量綱自振頻率Ω=ksL=ω1/2(mb/D)1/4L, 則損傷前后該無量綱自振頻率變化率可表示為

(17)

將式(2)、式(14)、式(15)及γ=iμ代入式(12)中,可得健康周期支撐結構第n階模態振型為

(18)

將式(2)、式(6)、式(7)和式(14)~式(16)及γ=iμ代入式(13)中可得如圖2所示周期支撐結構的頻率特征方程為

(19)

其中,

A1=a2b-bd2-ac2-b3-2bc2;

B1=a2c-cd2+abc-c3-2b2c;

C1=a2b-bd2-b3-2bc2;D1=c3+2b2c;

A2=bc2;B2=a2c-cd2-b2c;

C2=abc+b2c;D2=ac2-bc2

式中:δξj=ΔD/D為損傷單元j在損傷前后的剛度變化率。當δξj=0時,式(19)簡化為健康周期支撐結構的頻率特征方程

(20)

其中,

通過該頻率方程能計算得到周期支撐結構未損傷時的各階無量綱頻率值。

1.2 周期支撐結構無量綱頻率的敏感度系數

周期支撐結構波傳播常數與周期結構單元導納之間的關系為[11]

(21)

將式(15)代入式(21)中得

(22)

(23)

表1 健康周期支撐結構的前5階無量綱頻率值

圖3 10單元周期支撐結構各階無量綱頻率敏感度系數Fig.3 Sensitivities of non-dimensional frequencies for the first five modes of periodically-supported structures with 10 cell elements

在多單元損傷情況下,任意階無量綱自振頻率變化可視作單損傷引起頻率變化的線性疊加,若忽略高次項,則損傷引起的無量綱頻率總變化率為

(24)

(25)

式中:p為周期結構測量模態階數; [S]為無量綱頻率的敏感性矩陣。

通常測量獲得的結構自振頻率階數要少于周期支撐結構單元總數(即plt;N),這會導致式(25)中的敏感性識別方程組欠定,造成求解結果不唯一。為此,考慮到損傷不可能引起剛度增加,將本文周期支撐結構的損傷識別問題轉化成如下的非負最小二乘曲線擬合問題

(26)

運用適合該類問題的數值優化算法[12]求解式(26),可得到周期支撐結構損傷識別結果。

2 數值仿真研究

本節通過一個周期支撐梁與周期支撐法蘭管線損傷識別算例對本文方法進行驗證。

2.1 周期支撐梁模型數值仿真

考慮圖2所示的10單元周期支撐梁模型,其中,通過降低某基本周期單元的彈性模量來模擬該單元的整體剛度降低。假定該周期梁模型為矩形截面,寬、高均為0.03 m,基本周期單元(兩相鄰支座之間或支座與固定端之間)的長度為0.5 m,整體結構總長為5 m,密度為2.7×103kg/m3,泊松比為0.3,彈性模量為71 GPa。待識別損傷位置個數為10,分別為各基本周期單元的彈性模量,而損傷識別僅利用前5階模態頻率,以模擬實際應用中僅較少階模態可測情況。

對于該周期支撐梁模型,共考慮4種損傷工況,其中包含單損傷和雙損傷工,各工況具體損傷位置及損傷程度見表2。此外,表3為各損傷工況對應的損傷前后前5階無量綱自振頻率變化率,從表3可知,損傷引起的無量綱頻率改變較小,且不同模態對于損傷敏感程度不一樣,其中第1階和第5階頻率對于損傷較為敏感。本例損傷識別過程中,為模擬自振頻率測量噪聲影響,在損傷前后的真實自振頻率中均考慮了標準差為1%的高斯白噪聲。無測量噪聲及考慮測量噪聲情況下的損傷識別結果,如圖4和圖5所示,從圖可知,本文方法在兩種情況下均能較正確識別損傷位置及相應損傷程度。此外,兩圖的損傷識別結果也同樣反映出,無噪聲情況下損傷程度識別值與預設損傷約有20%的差別,噪聲工況下該差別更明顯。

表2 周期支撐梁模型損傷工況

表3 周期支撐梁模型損傷前后前5階無量綱頻率變化率

圖4 無測量噪聲情況下周期支撐梁的損傷識別結果Fig.4 Identified damage for the periodically-supported beam model without noise

圖5 考慮測量噪聲情況下周期支撐梁的損傷識別結果的誤差棒圖Fig.5 Error bars of identified damage for the periodically-supported beam model with measurement noise

2.2 周期支撐法蘭管線模型數值仿真

為進一步驗證本文周期支撐結構敏感性損傷識別方法對于更一般實際工程結構的適用性,具有10個周期單元的周期支撐法蘭連接管線模型,其左端C設為固定邊界,如圖6(a)所示。其中,將相鄰兩支座之間的法蘭管道視作一個基本周期單元,如圖6(b)所示。且規定左固定端與相鄰支座之間部分為1號周期單元。

假定法蘭盤之間通過螺栓連接,為簡化建模分析,將兩對接法蘭之間的螺栓連接簡化為薄層單元,即各基本周期單元的中間部分變為含薄層單元的法蘭接頭,見圖6(b)。在管道結構健康狀態下,各法蘭接頭中間薄層單元的彈性模量相等。通過降低某薄層單元的彈性模量,也即降低該法蘭薄層所在基本周期單元的局部剛度,以模擬實際管線工程中螺栓法蘭連接因螺栓松動而導致的局部剛度降低。

(a) 周期支撐法蘭管道模型示意

(b) 基本周期單元圖6 周期支撐法蘭管道簡化模型與其基本周期單元Fig.6 Schematic model for periodically-supported flanged pipeline and its cell element

本周期支撐法蘭管道橫截面為圓形,管道內外徑分別為0.081 m與0.089 m,法蘭外徑為0.18 m,法蘭厚度為0.02 m。相鄰法蘭接頭之間的管道長度為1 m,而整體模型最左端與最右端的管道長度均為0.5 m,如圖6(a)所示。假定每對法蘭配置6根M16×80 mm螺栓,兩對接法蘭盤面中間薄層部分厚度假定為0.004 6 m,彈性模量與質量密度依據文獻[13]方法確定。管道與法蘭材料相同,其彈性模量與質量密度均分別為68.7 GPa和2.8×103kg/m3,泊松比為0.33。圖7表示依據上述幾何材料參數與邊界條件建立的10周期的法蘭管道實體有限元模型。

針對該法蘭接頭管線模型,共考慮了4種損傷工況,如表4所示。其中,前兩種工況為單損傷,后兩種工況考慮雙損傷,且通過將各對應損傷法蘭接頭中間薄層單元彈性模量降低50%以模擬螺栓松動引起的管線局部法蘭接頭損傷。表5為通過圖7所示實體有限元模型計算得到的周期支撐法蘭接頭管線各工況下損傷前后前5階無量綱自振頻率變化率,且后續損傷識別過程僅利用該前5階無量綱頻率改變,以模擬實際工程應用中僅較少階模態可測的情況。

表4 周期支撐法蘭管道模型損傷工況

表5 周期法蘭管線模型傷前后前5階無量綱頻率變化率

圖7 含薄層單元的法蘭管道有限元模型Fig.7 FE model for periodically-supported flanged pipeline

應指出,本例所考慮的法蘭接頭損傷與前例中的周期支撐梁結構有所不同。前例中周期支撐梁考慮理論化的基本周期單元整體剛度降低,而本例考慮更為實際的周期單元內的局部剛度降低,且同時考慮法蘭接頭附加質量的影響。因而,前例中的識別損傷程度能與事先指定的真實值定量比較,而本法蘭管線算例中的識別損傷程度僅能相對表征。同樣地,損傷前后真實自振頻率均考慮了1%的高斯白噪聲影響,采用本文方法得到各損傷工況無噪聲與有噪聲情況下對應的識別結果,如圖8和圖9所示。從兩情況識別結果中可知,對于單元局部損傷且考慮法蘭接頭附加質量的周期支撐法蘭管線,各損傷工況下均能較準確地識別出損傷位置及相對損傷程度,表明本文方法具有較強的適用性。

圖8 無測量噪聲情況下周期支撐法蘭管道損傷識別結果Fig.8 Identified damage for the periodically-supported flanged pipeline model without noise

圖9 考慮測量噪聲情況下周期支撐法蘭管道損傷識別結果的誤差棒圖Fig.9 Error bars of identified damage for the periodically-supported flanged pipeline model with measurement noise

3 結 論

本文基于特征波導納和敏感性分析法得到有限長周期支撐結構無量綱自振頻率敏感性矩陣的通用表達式,進而建立了一種有效的周期支撐結構敏感性動力損傷識別方法。研究表明,充分利用周期支撐結構的固有周期特性,獲得的無量綱自振頻率敏感性矩陣僅與周期單元總數和結構損傷位置有關,而與結構其它幾何物理參數如密度、彈性模量、長度、橫截面尺寸以及橫截面形狀等無關,極大地簡化了復雜周期支撐結構特征值敏感性分析過程,避免了對于結構準確模型參數的依賴,顯著提高了本文損傷識別方法的實用性。對一周期支撐梁與周期支撐法蘭管線模型的多種損傷工況識別表明,本文方法均能較準確識別出損傷位置與損傷程度。

應指出,本文損傷識別方法基于周期支撐結構中基本周期單元整體剛度下降基礎上推導得來,因而,對于周期支撐梁算例中基本周期單元的整體剛度下降,本文方法能給出較準確的損傷程度估計,而對于基本周期單元的局部剛度降低,本文方法也能較準確地給出各位置的相對損傷程度,表明本文方法具有較好的適用性和擴展性。此外,也應注意到,盡管對于損傷位置判定較準確,但本文基于一階固有頻率敏感性分析的損傷識別方法對于損傷程度的準確估計仍顯不足,因此,進一步優化損傷敏感性參數與運用高階敏感性分析技術有助于提高本文方法對于損傷程度估計的準確性。

[ 1 ] ZHAO J, DEWOLF J T. Sensitivity study for vibrational parameters used in damage detection[J]. Journal of Structural Engineering, 1999, 125(4): 410-416.

[ 2 ] CAWLEY P, ADAMS R D. The location of defects in structures from measurements of natural frequencies[J]. Journal of Strain Analysis for Engineering Design, 1979, 14(2): 49-57.

[ 3 ] BISWAS M, PANDEY A K, SAMMAN M M. Diagnostic experimental spectral/modal analysis of a highway bridge[J]. International Journal of Analytical amp; Experimental Modal Analysis, 1990, 5(1): 377.

[ 4 ] WAHAB M M A, ROECK G D. Damage detection in bridges using modal curvatures: application to a real damage scenario[J]. Journal of Sound and Vibration, 1999, 226(2): 217-235.

[ 5 ] 薛松濤, 錢宇音, 陳镕, 等. 采用二階頻率靈敏度的損傷識別和試驗[J]. 同濟大學學報(自然科學版), 2003, 31(3): 263-267.

XUE Songtao, QIAN Yuyin, CHEN Rong, et al. Damage identification and experiments of frame structure based on second-order rrequency sensitivity[J]. Journal of Tongji University (Natural Science), 2003, 31(3): 263-267.

[ 6 ] YIN Tao, ZHU Hongping, YU Ling. Noise analysis for sensitivity-based structural damage detection[J]. Applied Mathematics and Mechanics, 2007, 28(6): 741-750.

[ 7 ] 朱宏平, 何波. 基于敏感性分析的周期結構損傷檢測[J]. 工程力學, 2003, 20(3): 108-114.

ZHU Hongping, HE Bo. Damage detection in periodic structures based on a sensitivity method[J]. Engineering Mechanics, 2003, 20(3): 108-114.

[ 8 ] ZHU H, WU M. The characteristic receptance method for damage detection in large mono-coupled periodic structures[J]. Journal of Sound and Vibration, 2002, 251(2): 241-259.

[ 9 ] ZHU H P, XU Y L. Damage detection of mono-coupled periodic structures based on sensitivity analysis of modal parameters[J]. Journal of Sound and Vibration, 2005, 285(1/2): 365-390.

[10] MEAD D J, BANSAL A S. Mono-coupled periodic systems with a single disorder: free wave propagation[J]. Journal of Sound and Vibration, 1978, 61(4): 481-496.

[11] MEAD D J. Wave propagation and natural modes in periodic systems: I. Mono-coupled systems[J].Journal of Sound and Vibration, 1975, 40(1): 1-18.

[12] LAWSON C L, HANSON R J. Solving least squares problems[M]. Englewood Cliffs: Prentice-hall, 1974.

[13] 姜東,吳邵慶,史勤豐,等.基于各向同性本構關系薄層單元的螺栓連接參數識別[J]. 振動與沖擊, 2014, 33(22): 35-40.

JIANG Dong, WU Shaoqing, SHI Qinfeng, et al. Parameter identification of bolted-joint using thin-layer element with isotropic constitutive relationship[J]. Journal of Vibration and Shock, 2014, 33(22): 35-40.

Damagedetectioninlargeperiodically-supportedstructuresbasedonthecharacteristicreceptancemethodandthesensitivity-basedapproach

YIN Tao, YIN Menglin

(School of Civil Engineering, Wuhan University, Wuhan 430072, China)

This paper developed a novel damage detection method for large periodically-supported structures based on the characteristic receptance method and the sensitivity analysis technique. By analyzing the free vibration of a finite periodically-supported structure with a single disorder based on the wave propagation method, the sensitivity matrix of the non-dimensional natural frequencies with respect to the change in element stiffness was obtained. And then, the damage scenarios in large periodically-supported structures were identified with the damage induced changes of natural frequencies by solving a set of underdetermined equations based on the sensitivity matrix. Furthermore, it was found that the sensitivities of the non-dimensional natural frequencies are independent of the structural physical parameters and thus any prior information of the original structures is never required. The proposed method is demonstrated by the numerical case studies conducted for both a periodically-supported beam and a periodically-supported flanged pipeline with various damage scenarios by utilizing only the frequency measurements for the first few modes before and after damage.

characteristic receptance; sensitivity analysis; periodic structures; damage detection; flanged pipeline

國家自然科學基金(51778506)

2016-05-10 修改稿收到日期: 2016-08-03

尹濤 男,博士,副教授,1979年生

O327

A

10.13465/j.cnki.jvs.2017.22.015

猜你喜歡
模態結構模型
一半模型
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
論《日出》的結構
3D打印中的模型分割與打包
國內多模態教學研究回顧與展望
創新治理結構促進中小企業持續成長
現代企業(2015年9期)2015-02-28 18:56:50
基于HHT和Prony算法的電力系統低頻振蕩模態識別
主站蜘蛛池模板: 亚洲有无码中文网| 91精品啪在线观看国产91九色| 国产无码制服丝袜| 欧美视频在线第一页| 国产永久在线视频| 污网站免费在线观看| 久久久久久久久18禁秘| 99视频只有精品| 波多野结衣爽到高潮漏水大喷| 在线观看无码av五月花| 天堂亚洲网| 91高清在线视频| 国产成人禁片在线观看| 97国产精品视频自在拍| 精品乱码久久久久久久| 天堂岛国av无码免费无禁网站 | 男女性午夜福利网站| 在线精品欧美日韩| 亚洲男人的天堂在线观看| 亚洲色精品国产一区二区三区| 国产主播喷水| 成人无码一区二区三区视频在线观看 | 国产网站一区二区三区| 综合久久五月天| 久久综合亚洲鲁鲁九月天| 狠狠色丁香婷婷| 精品丝袜美腿国产一区| 日本高清免费一本在线观看| 成年女人a毛片免费视频| 扒开粉嫩的小缝隙喷白浆视频| 国产a v无码专区亚洲av| 国产精品亚洲一区二区三区z| 99热这里都是国产精品| av在线无码浏览| 国产精品浪潮Av| 无码日韩人妻精品久久蜜桃| 国产视频一二三区| 亚洲一区二区约美女探花| 国产精彩视频在线观看| 亚洲精品不卡午夜精品| 亚洲国产高清精品线久久| 亚洲成A人V欧美综合| jijzzizz老师出水喷水喷出| 国产精品免费电影| 国产无码网站在线观看| 亚洲男人的天堂久久精品| 日韩一区精品视频一区二区| 婷婷综合亚洲| 在线观看免费黄色网址| 巨熟乳波霸若妻中文观看免费| 免费无码AV片在线观看国产| 毛片最新网址| 青青草国产在线视频| 动漫精品中文字幕无码| 天天视频在线91频| 久久96热在精品国产高清| 波多野结衣亚洲一区| 99中文字幕亚洲一区二区| 国产69精品久久久久孕妇大杂乱| 亚洲一区二区三区在线视频| 曰AV在线无码| 亚洲一区精品视频在线| 国产欧美另类| 91精品国产麻豆国产自产在线 | 丝袜国产一区| 国产伦精品一区二区三区视频优播| 亚洲天堂网站在线| 欧美日韩久久综合| 亚洲欧美精品一中文字幕| 久久99久久无码毛片一区二区 | 国产无码高清视频不卡| 香港一级毛片免费看| 国产黄色视频综合| 黄色一级视频欧美| 九九线精品视频在线观看| 欧美影院久久| 国产成人综合网| 国产91精品久久| 在线免费看黄的网站| 欧美不卡视频在线| 欧美a级在线| 免费观看无遮挡www的小视频|