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

基于改進交叉模型交叉模態法的局部損傷識別方法

2016-01-07 00:57:59占超,李東升,任亮
振動與沖擊 2015年7期

第一作者占超男,碩士生,1990年生

通信作者李東升男,副教授,1972年生

郵箱:dsli@dlut.edu.cn

基于改進交叉模型交叉模態法的局部損傷識別方法

占超, 李東升, 任亮, 李宏男

(大連理工大學海岸及近海工程國家重點實驗室,遼寧大連116024)

摘要:傳統的交叉模型交叉模態(CMCM)法由于其核心矩陣的缺秩使其在全局修正時的解不唯一,為得到唯一解必須人為假定約束。以往利用CMCM法進行損傷識別的研究中,通常將損傷前后質量保持不變作為約束求得各單元剛度的修正量,然后將剛度的下降作為判斷損傷的依據。然而,結構在極端運營條件下質量的變化也是損傷的一種。為此,對傳統的CMCM法進行了改進,提出了以下的改進方法:首先,由基準狀態下的有限元模型和損傷后實測的結構低階模態信息求得核心矩陣,將核心矩陣最小奇異值對應的右奇異向量作為損傷指示向量(DIV);然后,基于損傷的局部性利用聚類分析算法自動識別DIV中的異常元素,將異常元素作為損傷定位的依據;最后,由DIV估計各單元質量和剛度的損傷程度。改進后方法的優勢在于:無需人為主觀地添加約束即可一次性求出結構各單元質量和剛度損傷前后的變化,避免了錯誤或不準確的假設給損傷識別結果帶來的誤差。通過數值實驗對該改進方法的可行性、魯棒性以及損傷敏感性進行了研究,并通過美國洛斯阿拉莫斯國家實驗室(LANL)的一個4自由度板柱結構振動臺實驗對改進方法的有效性做了進一步的驗證。

關鍵詞:損傷識別; 全局修正; 交叉模型交叉模態法; 結構健康監測

收稿日期:2014-01-02修改稿收到日期:2014-04-03

中圖分類號:TP206+.3

文獻標志碼:A

DOI:10.13465/j.cnki.jvs.2015.07.021

Abstract:Since the core matrix in the traditional cross-model cross-mode(CMCM) method is rank-deficient, constraints must be applied artificially for getting a unique solution in a complete model-updating situation. In previous damage identification researches based on the CMCM method, structural masses are often supposed to be unchanged for getting the updated coefficients of stiffness, and the reduction of stiffness is regarded as an indicator of damage. However, elemental mass changes under operational conditions are obviously a kind of damages. For detecting damage more effectively, the traditional CMCM method was improved, and an improved method was proposed as follows: evaluate the core matrix from the finite element model under the baseline condition and the measured modal data after damage and take the right singular vector corresponding to the least singular value of the core matrix as a damage indication vector(DIV); identify the abnormal elements in the DIV by using cluster analysis algorithm based on the assumption that damage appears locally and locate damage based on these abnormal elements; derive the extent of damage for each elemental mass and stiffness from the DIV. The advantage of the improved method is its ability to solve the mass and stiffness changes before and after damage for all elements without artificial constraints, thus the error of damage identification caused by imprecise or wrong constraints vanishes. The improved method’s practicability, robustness, and sensitivity to damage were analysed via a numerical simulation, and its effectiveness is further verified by virtue of the experiments on a 4 degrees of freedom test-bed structure done at the Los Alamos National Laboratory.

Local damage identification approach based on improved cross-model cross-mode method

ZHANChao,LIDong-sheng,RENLiang,LIHong-nan(State Key Laboratory of Coastal and Offshore Engineering, Dalian University of Technology, Dalian 116024,China)

Key words:damage identification; complete-updating; cross-modal cross-mode; structural health monitoring

隨著我國經濟的不斷發展,大量土木工程結構(尤其是高層建筑和大跨橋梁)伴隨著人們的需要如雨后春筍般涌現。這些結構在以后的服役過程中,由于各種因素(如地震、材料銹蝕等)某些構件往往會出現損傷。損傷的存在勢必會降低結構的使用性能,放任損傷的加劇甚至可能引發重大的土木工程結構事故(如橋梁的突然折斷、房屋的驟然倒塌等),從而造成重大的人員傷亡和財產損失。因此,對重大工程結構的性能進行監測和診斷,及時發現結構的損傷,對可能出現的災害進行預測并對結構的安全性能進行正確的評價,具有極其重要的意義,同時也是現代工程的迫切要求和土木工程學科發展的一個重要領域[1]。

結構健康監測和損傷識別近幾十年來在硬件和理論上均有很大進展[2],其中基于模型修正的損傷識別方法受到諸多學者的關注[3]。該類方法的基本思想為根據健康狀態下的實測振動數據對結構的初始有限元模型進行修正,得到基準模型;結構損傷后,由損傷狀態下的實測振動數據對基準模型進行修正,則損傷前后各單元模型參數的變化即反映出損傷的位置和程度。

傳統的模型修正方法[5]主要有:利用模態信息的矩陣直接修正法、利用模態信息的迭代修正法和利用實測頻響函數的修正法等。這些方法或多或少存在各種缺點,比如直接修正法修正后的質量、剛度矩陣無法保證原有的對稱性、稀疏性、和正定性,缺乏物理意義;迭代修正法存在解的收斂性以及迭代過程中的誤差累計問題,并且計算量較大;而頻響函數修正法的修正參數估計為有偏估計。為了克服傳統方法的不足,近年來不斷有學者提出新的模型修正方法。Doebling等[6]提出通過敏感性分析得到單元連接矩陣,然后由最小秩法求出單元剛度參數的擾動量,并使修正后的剛度矩陣保持了原有的單元連接形式;Kaouk等[7]提出將模型修正分為兩步,即先進行誤差定位,再使用最小秩擾動法對結構質量、阻尼、剛度矩陣同時修正;侯吉林[8]提出基于局部模態的約束子結構模態修正法,并對此進行了系統地研究。

Hu等[9]提出了交叉模型交叉模態(CMCM)算法。該方法可以根據測得的有限低階模態信息對單元質量和單元剛度矩陣同時修正,而且無需模態配對,對振型沒有同比例要求,是一種直接的物理屬性修正方法[10]。本文先對傳統的CMCM法推導過程進行簡要介紹;然后指出其中的不足并提出改進方法;接著通過一個數值實驗對改進方法的可行性、魯棒性和損傷敏感性進行研究;最后通過美國洛斯阿拉莫斯國家實驗室(LANL)的一個4自由度板柱結構實驗對改進方法的可行性做進一步驗證。

1交叉模型交叉模態法

記有限元模型的剛度和質量矩陣分別為K和M,則

KΦ=MΦΛ

(1)

式中Φ為有限元模型振型矩陣,Λ為與各振型對應的圓頻率的平方值構成的Ni階對角陣。

類似式(1),對于實際結構有

K*Φ*=M*Φ*Λ*

(2)

式中,K*和M*分別為實際結構的剛度和質量陣,Φ*實測振型矩陣,Λ*為實測圓頻率的平方值構成的Nj階對角陣。

式(2)等號兩邊同時左乘的Φ轉置,即

ΦTK*Φ*=ΦTM*Φ*Λ*

(3)

假設實際結構的剛度和質量矩陣滿足下列方程,

(4)

式中Kn與Mn分別為在整體坐標系下標號為n的單元剛度和單元質量矩陣,αn和βn分別為n號單元剛度矩陣修正系數和n號單元質量矩陣修正系數,NK和NM分別為單元剛度矩陣和單元質量矩陣的個數。

將式(4)代入式(3),可得

(5)

式中

(6)

式(5)可以寫成如下矩陣形式,

Gγ=f

(7)

式中

(8)

式(7)中矩陣G的行數,亦方程數Nm=NiNj;其列數,亦未知數個數Nu=NK+NM。由于矩陣G是CMCM算法的核心部分,因此也被稱為核心矩陣[11-12]。需要指出的是,當實測模態測得足夠準時,即使式(7)中方程數Nm大于未知數個數Nu,由于G的缺秩,式(7)仍會出現解不唯一的情況。要想得到唯一解,必須添加約束。假設約束方程如下,

(9)

若假設結構質量保持不變,則

(10)

式中,I為維數為NM的單位方陣。將式(7)與式(10)合并,即得

(11)

式中

(12)

(13)

2改進的CMCM局部損傷識別法

在第1節介紹的傳統CMCM法中,由于核心矩陣G的缺秩,使其在全局修正時解不唯一。要想得到唯一解,至少需要人為地添加一個約束。錯誤或者不精確的約束都可能會引起傳統CMCM法修正結果出現誤差。以往在對基于CMCM法的損傷識別研究中,往往假設各單元質量保持不變,而僅對各單元剛度進行修正,并以單元剛度的下降作為損傷判斷的依據。然而,結構在極端運營條件下質量的變化(例如暴雪導致房屋頂層質量的顯著增加等)也是損傷的一種。本節對傳統CMCM法的推導過程稍加變動,提出改進的CMCM局部損傷識別法,如下:

假設損傷后剛度、質量陣K*和M*滿足下列關系

(14)

式中ηn和εn分別為n號單元剛度矩陣和n號單元質量矩陣的“損傷比”,與式(4)中的αn和βn具有如下關系

ηn=1+αn,εn=1+βn

(15)

將式(14)代入式(3),可得

(16)

將式(16)寫成矩陣形式

Gv=0

(17)

其中Cn、Dn和G分別同式(6)和式(8)的定義。v定義如下

(18)

向量v的前半段和后半段分別記錄了各單元剛度、質量損傷前后的比值,求得了v即實現了對結構參數損傷前后變化的全面了解。由式(17)可見,只要v有解,核心矩陣G必缺秩,v必有無窮多個解。求解v的問題即化為求解G的零空間的問題。

數學上往往采用奇異值分解的方法求解一個矩陣的零空間,即

G=USVT

(19)

式中S為Nu×Nu的對角矩陣,其對角線上記錄了G的Nu個從大到小排列的奇異值s1,s2,…,sNu;U(Nm×Nu)和V(Nu×Nu)中各列向量分別稱為G的各奇異值對應的左奇異向量和右奇異向量,并滿足

UTUINu,VTV=INu

(20)

則奇異值s1,s2,…,sNu中0奇異值對應的右奇異向量即為零空間的基。

然而,由于儀器的內部噪聲、計算機的舍入誤差、環境噪聲等因素,實際測得的模態中含有噪聲而使得G趨于滿秩,即所有奇異值均大于0。為此,我們取奇異值中最小的值sNu對應的右奇異向量作為解v的基,記為v*(v*滿足v*Tv*=1)。

v與v*間滿足倍數關系,即v=av*,其中a為待定系數。我們知道,損傷是一種典型的局部現象,兩次動力測試中結構大部分的參數應該基本保持不變,只有損傷處參數會有較顯著的變化,即:v中大部分元素都應該近似等于1,對應地v*中大部分元素應處于同一水平;損傷情況下v*中會出現少數異常元素,表現在與其他正常元素不在同一水平線上。v*中的異常元素可以作為損傷定位的依據。由于向量v*中包含結構的損傷信息,我們又稱v*為損傷指示向量(Damage Indication Vector,DIV)。DIV中正常元素和異常元素的劃分可以通過聚類分析算法[13]自動實現。聚類分析能夠從樣本數據出發,計算樣本之間的空間距離,根據樣本彼此間的親疏關系,對樣本進行自動分類。

(21)

(22)

圖1 四自由度鏈式 結構(數值實驗) Fig.1 4-DOF chain-like system (for numerical simulation)

則av*即為我們對解v的估計,(av*-1)×100%即為我們對各單元剛度、質量損傷程度的估計。

3數值實驗

3.1改進方法的可行性

本小節通過數值實驗對第2節中改進方法的可行性進行說明。如圖1所示一四自由度鏈式結構,底層柱與小車相連,結構整體可以在水平地面上自由滑動。各集中質量僅考慮其在水平方向上的運動。各層質量、小車質量以及各層柱的水平剛度如圖中所示。實驗中,結構的單處損傷通過將第一層的層間剛度k1折減50%實現;多處損傷工況下則將第二層的層間剛度k2折減20%,并將第一層的集中質量增加20%以模擬極端運營條件下質量的增加。數值實驗中損傷狀態下的模態是通過損傷后的有限元模型直接計算得到的。結構基準狀態下的模態和損傷狀態下的“實測”模態分別如表1和表2所示。結構的第一階模態為平動模態,在表中不予顯示。

表1 基準狀態下有限元模型的模態

將表1和表2中的模態數據代入式(17),即可求得損傷指示向量(DIV),進而求得損傷位置和程度。無噪聲情況下,單損傷工況和多損失工況的損傷識別結果如圖2所示。圖中可見,無論是單損傷工況還是多損傷工況下,DIV中異常元素均正確指示出損傷位置。單損傷工況中,η1被標記為異常元素,表明一號單元剛度出現異常;多損傷工況中,η2和ε2均被標記為異常,表明二號單元剛度以及二號單元質量(即圖1中一層質量)出現了異常。兩種工況下,損傷程度識別值均與預設值均很好地吻合,從而說明改進方法可行。

表2 損傷狀態下的結構模態

圖2 損傷識別結果(數值實驗) Fig.2 The results of damage identification(numerical simulation)

3.2改進方法的魯棒性和損傷敏感性

有關小損傷的識別一直是損傷識別中的難點。小損傷一方面對結構整體參數影響不大,難以從結構的動力響應信號中辨別;另一方面又對結構安全造成隱患,可能隨時發展成大損傷進而導致結構構件失效甚至是結構整體的破壞。同時,實際工程中來自方方面面的噪音也對損傷的正確識別造成困難。本小節以多損傷工況(圖1中第二層層間剛度折減y%,第一層集中質量增加y%,y為損傷程度)為例,對改進方法的魯棒性和損傷敏感性進行研究。

圖3 改進方法識別性能灰度圖 Fig.3 The grayscale of the improved method’s performance

我們以實測模態中的噪聲比例x為橫軸(分辨率0.1%),以多損傷工況中的損傷程度y為縱軸(分辨率1%),作出改進方法的識別性能灰度圖,如圖 3所示。圖中各點(共50×20個點)均重復實驗100次。每次實驗中,首先由該點縱坐標對應的損傷程度求得對應的損傷有限元模型和損傷后的模態數據,然后在模擬的振型中加入該點橫坐標對應比例的比例噪聲,接著由含噪聲的模態數據利用改進方法求得DIV,若DIV正確指示出了損傷位置,則記為一次成功的損傷識別。100次實驗中,成功識別次數z越多,則識別性能越好,在圖中對應的顏色也就越深。圖中可以看到,識別性能隨著損傷程度的減小和噪音比例的增大而變差。無噪聲條件下,任意程度的損傷均可準確識別; 1%比例噪聲下,改進方法可以對10%以上的損傷準確定位(z>95);而在5%的比例噪聲下,小損傷(y<5)的識別則完全無法實現。

4實驗室實驗

為進一步證實改進方法的可行性,本節對LANL實驗室的一個4自由度板柱結構振動臺實驗數據[14-15]進行分析,并利用改進方法對預設的損傷進行識別。實驗裝置如圖4所示。三層“樓板”和“基礎”均為鋁制金屬板(30.5 cm×30.5 cm×2.5 cm),彼此通過四根鋁柱使用螺栓連接?!盎A”被限制在軌道上使其僅可以在x方向上滑動。激振器采用(20 Hz~150 Hz)的限帶白噪聲對“基礎”進行激勵。結構上共安裝5個傳感器(通道1用于測量激振器對“基礎”施加的力,通道2、3、4、5分別用于測量“基礎”和一、二、三層“樓板”的加速度響應),各傳感器采樣頻率均為320 Hz,采樣數目為8 192個。實驗中將損傷歸為兩類,一類是線性損傷,另一類是非線性損傷。線性損傷通過在樓板或基礎上添加集中質量塊、替換某層的一根或兩根鋁柱(替換后的鋁柱的橫截面面積為原始的一半)以減小層間剛度來實現;非線性損傷[16]則通過調節頂層懸掛鋁柱和安裝在第二層樓板上的“碰撞器”之間的間距,在動荷載下二者之間不同程度的相互碰撞來模擬。本文僅對其中的線性損傷數據進行分析,如表3[17]所示。結構基準狀態下的各傳感器時程如圖5所示。

圖4 LANL 4自由度振動臺試驗裝置圖 Fig.3 Experiment setup for 4-DOF test-bed structure in LANL

表3 實驗分組表 [17]

圖5 基準狀態下各傳感器所測信號時程 Fig.5 Time histories of signals in each sensor under baseline condition

我們將激振器施加的力作為系統輸入,將“基礎”和各層“樓板”的加速度信號作為系統輸出,使用OKID[18](Observer Kalman filter IDentification)算法進行模態識別。由于實測信號中的噪音和算法階數的不同取值,測得的模態中含有虛假模態和計算模態,即非物理模態。模態分析中可以通過穩態圖[19]的方法鑒別物理模態和非物理模態。與非物理模態不同的是,物理模態不會隨算法階數的增加而變化,因而在穩態圖中呈現近乎一條直線。圖6為基準狀態下模態識別的穩態圖。圖中每個實心點對應一個識別的模態(為了使穩態圖更加清晰,我們濾去了阻尼比大于5%的模態);點畫線為力對第三層樓板加速度的頻響函數的幅值。圖中可見,不同算法階數下物理模態的頻率具有較好的一致性,如虛線所示。限于篇幅,損傷狀態下各傳感器的時程和穩態圖不再列出。不同工況下測得的模態具體數值見表 4。結構的有限元模型參數與數值實驗中圖1所示的各參數相同。將表4中基準狀態下的實測模態與表1中的有限元模態相對比,可以發現二者比較吻合,但仍略有差異。為此我們將初始有限元模型和基準狀態下的實測模態代入式(17)求得v*,然后由式(22)對初始有限元模型進行全局修正。有限元模型修正參數及修正效果見表5。表中可見,修正后的有限元模型各階振型與實測振型的相似度(即MAC值)相比修正前均有所提高、第二階頻率誤差減小,但第三階和第四階的頻率誤差修正后反而略有增大。這主要是因為修正前,第三階和第四階的頻率誤差已經很小(小于0.5%),而全局修正是對結構模型的全局優化,因而使其誤差反而略有增大。

圖6 結構基準狀態下模態識別穩態圖 Fig.6 Stability diagram under baseline condition

圖7 損傷識別結果(LANL) Fig.7 The results of damage identification(LANL)

將表 4中各損傷狀態下的實測頻率與基準狀態下的實測頻率相對比,可以發現損傷后結構各階頻率均有不同程度的下降。將損傷狀態下的實測模態和修正后的有限元模型代入式(17)求得DIV,進而對損傷進行定位和定量,識別結果如圖 7所示。圖中可見,各損傷工況下損傷指示向量中的異常元素均正確指示出損傷位置,損傷程度的估計也與表3中給出的預設損傷值很好地吻合(兩者間的誤差不超過5%)。我們還對工況2利用第1節中的傳統CMCM法進行了計算, 假設結構各單元質量不發生變化,結果得出第一、二層的層間剛度大約下降10%,第三層層間剛度略有增大的錯誤結論。改進的CMCM局部損傷識別法中則無須人為主觀地給出約束,避免了錯誤或者不準確的約束給損傷識別結果帶來的誤差。

表4 不同工況下實測模態

表5 基準狀態下有限元模型全局修正

5結論

交叉模型交叉模態法是一種很好的模型修正方法,具有諸多優點,然而由于其方法中核心矩陣的缺秩,使得全局修正時必須人為添加約束。本文對傳統CMCM法進行了改進,提出了以下的改進方法:首先,由基準狀態下的有限元模型和損傷后實測的結構低階模態信息求得核心矩陣,將核心矩陣最小奇異值對應的右奇異向量作為損傷指示向量(DIV);然后,基于損傷的局部性利用聚類分析算法自動識別DIV中的異常元素,將異常元素作為損傷定位的依據;最后,由DIV估計各單元質量和剛度的損傷程度。改進后方法的優勢在于:無需人為主觀地添加約束即可一次性求出結構各單元質量和剛度損傷前后的變化,避免了錯誤或不準確的假設給損傷識別結果帶來的誤差。數值實驗和美國LANL實驗室實驗均證實該改進方法切實可行,可以對損傷進行正確的定位并對損傷程度進行準確的估計(誤差不超過5%)。數值實驗還表明該改進方法可以在實測振型含1%的比例噪聲下對損傷程度大于10%的多損傷工況進行準確的損傷定位。關于如何進一步提高改進方法的魯棒性和對小損傷的敏感性,還有待進一步的研究工作。

參考文獻

[1]李宏男,李東升. 土木工程結構安全性評估、健康監測及診斷述評[J]. 地震工程與工程振動,2002,22 (3): 82-90.

LI Hong-nan, LI Dong-sheng. Safety assessment, health monitoring and damage diagnosis for structures in civil engineering[J]. Earthquake Engineering and Engineering Vibration, 2002,22(3):82-90.

[2]宗周紅, 任偉新, 阮毅. 土木工程結構損傷診斷研究進展[J]. 土木工程學報, 2003, 36(5): 105-110.

ZONG Zhou-hong, REN Wei-xin, RUAN Yi. Recent advances in research on damage diagnosis for civil engineering structures[J]. China Civil Engineering Journal, 2003, 36(5):105-110.

[3]Mottershead J E, Friswell M I. Model updating in structural dynamics: a survey[J]. Journal of Sound and Vibration, 1993, 167 (2): 347-375.

[4]李輝, 丁樺. 結構動力模型修正方法研究進展[J]. 力學進展, 2005,35(2): 170-180.

LI Hui, DING Hua. Progress in model updating for structural dynamics[J]. Advances in Mechanics, 2005, 35(2): 170-180.

[5]Friswell M, Mottershead J E. Finite element model updating in structural dynamics[M]. Springer, 1995.

[6]Doebling S W. Minimum-rank optimal update of elemental stiffness parameters for structural damage identification[J]. AIAA journal, 1996, 34 (12): 2615-2621.

[7]Kaouk M, Zimmerman D C, Simmermacher T W. Assessment of damage affecting all structural properties using experimental modal parameters[J]. Journal of Vibration and Acoustics, 2000, 122 (4): 456-463.

[8]侯吉林. 約束子結構模型修正方法[D]. 哈爾濱:哈爾濱工業大學, 2010.

[9]Hu S L J, Li H, Wang S Q. Cross-model cross-mode method for model updating[J]. Mechanical Systems and Signal Processing, 2007, 21 (4): 1690-1703.

[10]賈輝,李東升. 基于交叉模型交叉模態法的損傷識別研究[J]. 土木工程與管理學報,2011,28 (3): 133-137.

JIA Hui, LI Dong-sheng. Damage identification inversitigation based on cross-model cross-mode method[J]. Journal of Civil Engineering and Management, 2011, 28(3):133-137.

[11]王俊榮. 海洋平臺結構物損傷檢測與模型修正方法研究[D]. 青島:中國海洋大學, 2009.

[12]王俊榮, 張敏, 李華軍. 一種基于不完備模態信息的有限元模型修正方法[J]. 工程力學, 2010, 27 (8): 60-65.

WANG Jun-long, ZHANG Min, LI Hua-jun. A finite element model updating approach based on incomplete modal data[J]. Engineering Mechanics, 2010, 27(8): 60-65.

[13]HE Zeng-you, XU Xiao-fei, DENG Sheng-chun. Discovering cluster-based local outliers[J]. Pattern Recognition Letters, 2003, 24 (9-10): 1641-50.

[14]Figueiredo E, Flynn E. Three-story building structure to detect nonlinear effects.http://institute.lanl.gov/ei/damage_id/data/4DOF_Mat_Format.zip. 2009.

[15]Figueiredo E, Park G, Figueiras J, et al. Structural health monitoring algorithm comparisons using standard data sets[R]. LA-14393, Los Alamos, NM (United States): Los Alamos National Laboratory, 2009.

[16]Figueiredo E, Todd M D, Farrar C R, et al. Autoregressive modeling with state-space embedding vectors for damage detection under operational variability[J]. International Journal of Engineering Science, 2010, 48 (10): 822-34.

[17]Miguel R, Hernandez-Garcia, Masri S F, et al. An experimental investigation of change detection in uncertain chain-like systems[J]. Journal of Sound and Vibration, 2010, 329 (12): 2395-2409.

[18]Juang J N. Applied system identification[M]. New Jersey: PTR Prentice Hall Englewood Cliffs Inc, 1994. 175-228.

[19]沃德·海倫,斯蒂芬·拉門茲,波爾·薩斯. 模態分析理論與實驗(白化同,郭繼忠譯)[M]. 北京:北京理工大學出版社,2001.

主站蜘蛛池模板: 米奇精品一区二区三区| 免费高清毛片| 思思热精品在线8| 日本高清成本人视频一区| 日本伊人色综合网| 亚洲av无码专区久久蜜芽| 免费观看国产小粉嫩喷水| 中文字幕啪啪| 亚洲区一区| 国产精品亚洲一区二区三区在线观看| 国产亚洲精品无码专| 中文字幕欧美日韩| 国产爽爽视频| 国产成人乱码一区二区三区在线| 97国产在线观看| 成人在线观看一区| 国产美女无遮挡免费视频| 在线免费a视频| 日本中文字幕久久网站| 久久久噜噜噜久久中文字幕色伊伊| 在线无码九区| 精品国产网| 亚洲欧洲日韩久久狠狠爱 | 26uuu国产精品视频| 日本人又色又爽的视频| 高清无码手机在线观看| 欧美成人午夜视频免看| 国内99精品激情视频精品| 国产亚洲男人的天堂在线观看| 国产杨幂丝袜av在线播放| yy6080理论大片一级久久| 呦女亚洲一区精品| 久久精品丝袜高跟鞋| 亚洲视频无码| 91免费国产在线观看尤物| 亚洲Av综合日韩精品久久久| 国产一级视频在线观看网站| 色男人的天堂久久综合| 精品无码人妻一区二区| 色播五月婷婷| 一区二区三区国产| 久久精品人人做人人| 中文字幕2区| 国产爽妇精品| 欧美特黄一级大黄录像| 欧美伊人色综合久久天天| 亚洲精品亚洲人成在线| 国产91视频免费| 欧美特黄一级大黄录像| h网站在线播放| 国产三级精品三级在线观看| 91九色最新地址| 国产精品第页| 欧美日韩精品在线播放| 色九九视频| 99久久免费精品特色大片| 国内黄色精品| 亚洲成人77777| 日韩区欧美国产区在线观看| 72种姿势欧美久久久大黄蕉| 欧美日韩国产在线人| 国产一二三区在线| 国产偷倩视频| 成人福利在线免费观看| 在线免费a视频| 久久五月天国产自| 免费高清a毛片| 波多野结衣无码中文字幕在线观看一区二区 | 99色亚洲国产精品11p| 久久综合AV免费观看| 精品福利国产| 欧美成人a∨视频免费观看| 国产精选小视频在线观看| 亚洲色图欧美在线| 亚洲无码电影| 亚洲天堂首页| 免费看a毛片| 国产精品久久久久久久久久久久| 97国产精品视频自在拍| 日韩a在线观看免费观看| 国产亚洲精| AV熟女乱|