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

利用固有頻率異常值分析法檢測螺栓擰緊力

2015-05-25 00:34:10緱百勇陸秋海王波王世英
振動與沖擊 2015年23期
關鍵詞:結構檢測方法

緱百勇,陸秋海,王波,王世英

(1.清華大學航天航空學院,北京100084;2.空軍工程大學航空航天工程學院,西安710038)

利用固有頻率異常值分析法檢測螺栓擰緊力

緱百勇1,2,陸秋海1,王波1,王世英1

(1.清華大學航天航空學院,北京100084;2.空軍工程大學航空航天工程學院,西安710038)

為了定量評估螺栓松動損傷程度,依據統計學中的馬氏平方距離(Squared Mahalanobis Distance,MD)異常值分析法,研究了綜合利用結構前5階固有頻率(稱為固有頻率矢量)檢測螺栓擰緊力的方法。首先測試了不同擰緊力下結構的固有頻率;然后計算了各擰緊力工況與緊固工況固有頻率矢量間的標準歐幾里德距離(Standard Euclidean Distance,SED),利用Matlab擬合了SED與擰緊力之間的關系曲線,據此提出了定量檢測螺栓擰緊力的方法;最后對該方法進行了試驗驗證和應用舉例。結果表明,固有頻率矢量的SED與擰緊力呈指數關系,而且利用此關系曲線能夠比較可靠地評估螺栓松動程度。

螺栓松動;損傷檢測;結構健康監測;固有頻率;異常值分析;標準歐幾里德距離

螺栓連接與鉚接、焊接、膠結等結構連接形式相比較,能夠承受大的載荷而且可重復裝配和拆卸,所以被大量用于連接可拆卸的主承力結構[1]。航空(天)器、建筑工地的吊塔、海洋平臺的支撐架以及鋼架橋等工程結構的主體受力部分基本都是通過螺栓連接的。因為預緊力可以提高螺栓連接的可靠性、防松能力和螺栓的疲勞強度,增強連接的緊密性和剛性,所以很多結構對螺栓預緊力有嚴格的要求。結構服役前,螺栓都按照要求處于擰緊狀態。然而,在長期的服役過程中,由于受到疲勞、振動、沖擊等各種載荷的作用,結構的連接螺栓有可能會逐漸松動(表現為螺栓擰緊力的下降),輕者造成設備故障,重則結構完全破壞,從而導致嚴重的災難。因此需要發展一種無損檢測方法在線評估螺栓連接的實際狀態[2]。

徐超等[3]和Wang等[4]分別就螺栓松動損傷檢測方法進行了回顧和總結。提及的方法概括起來主要分四類:①基于聲彈性效應的方法。該方法利用螺栓的聲彈性特性會隨螺栓擰緊力的變化而發生改變這一現象來檢測螺栓的預緊力[5-7]。其缺點是對信號采集設備的采樣頻率要求很高(約1 Gs/s),不但設備昂貴,而且高頻信號易受噪聲的干擾。②基于壓電傳感器的方法。該方法又可分為基于壓電導納(阻抗)的方法和基于能量的方法。前者早期被Sun等[8]用在桁架結構連接處的松動識別上,王丹生等[9]和高峰等[10]則將其用在螺栓松動損傷的檢測上。該方法通過檢測緊貼在結構上的壓電陶瓷片的電阻抗(或導納)來反映結構阻抗的變化,從而實現對基體結構上螺栓松動與否的識別。Wang等[11]則對基于能量的方法進行了研究。該方法利用了主動傳感技術,將螺栓連接構件一側的壓電片作為激勵源,另一側的壓電片作為拾能器,根據兩個壓電片之間傳遞能量的大小反映螺栓松緊程度。基于壓電傳感器的方法利用的也是高頻信號,除了設備昂貴、易受干擾之外,傳感器與被測結構連接的可靠性和耐久性問題也需要解決。③基于非線性動力學理論的方法。在較大外激勵作用下,螺栓松動造成結構連接處界面分離,產生非線性動態響應,直接利用結構的非線性動態響應信號,基于非線性動力學理論方法提取描述結構損傷本質的特征參數,建立起結構損傷狀態與非線性特征參量之間的關系,進而進行狀態監測和辨識。如Timothy等[12]基于混沌理論研究了螺栓松動狀態的識別問題。這類方法測試原理和測試設備都很復雜,不便于工程應用。④基于線性動力學理論的方法。結構損傷會引起結構動力學參數(如固有頻率、振型、阻尼、曲率模態、頻響函數、功率譜等)的變化,測得這些動力學參數變化量可反推結構的損傷狀態。楊智春等[13-14]、Cawley等[15]和Liang等[16]的研究即可歸為此類方法。

由于結構動力學參數測試方便而且精度高[17-18],所以,成為狀態監測和識別的經典方法。不過,利用結構動力學參數識別螺栓松動損傷的研究中,大都將螺栓分為“松”和“緊”兩種狀態來研究,定量檢測擰緊力的文獻還很缺乏。

本文從統計學中基于馬氏平方距離(the squared Mahalanobis Distance,MD)的異常值分析法出發,提出了一種基于結構前5階固有頻率標準歐幾里德距離(Standard Euclidean Distance,SED)的螺栓松動定量檢測方法,并進行了試驗驗證和應用舉例。

1 基于SED的異常值分析法

在基于數據驅動的損傷檢測方法中,如果將結構完好情況下測得的特征參數定義為正常數據,損傷情況下測得的特征參數定義為異常數據,損傷檢測的本質就是從實測的大量數據中將異常數據分離出來,并盡可能從異常數據中提取更多有關損傷的信息(損傷程度、損傷位置等)。由于受測量誤差等因素的影響,實測數據存在一定的分散性,即使結構在完好情況下,每次測得的特征參數數據也不可能完全一樣,這給異常數據的分離帶來了困難。基于馬氏平方距離的異常值分析法正是解決這一問題有效方法。

如果實測特征參數是一個單變量分布,其均值為μ,標準差為σ,則樣本x到均值的距離有兩種定義方式;①x-μ,稱為絕對距離,②(x-μ)/σ,稱為相對距離。顯然,相對距離是一種標準化的距離,不僅反映了該樣本與均值的差別程度,還可從分布概率的角度反映出該樣本在所有樣本中的“位置”。

拓展開來,對一個均值為μ,協方差矩陣為Σ的多變量分布而言,其樣本x到均值的相對距離可以定義為Z=L-1(x-μ),其中有LLT=Σ。該定義不僅考慮了各變量的分布情況,還考慮了不同變量之間的相互影響。為了去掉計算中的平方根計算,取上述相對距離的平方值,即定義D=ZTZ,則

最后一行即為馬氏平方距離的定義。

要將馬氏平方距用在損傷檢測中,首先要從實測的數據中提取能夠反映損傷狀態的多變量特征參數(本文以結構前5階固有頻率組成的向量為特征參數)。假設實測多變量特征參數由p個變量組成,且有n個觀測值,則可看作一個p維空間的n個測點,也可以表示為p×n維的矩陣。利用上述定義,第i個測點的馬氏平方距離為:

式中:Di為異常值指標;{fi}為實測p維變量{f}的第i次測量值,為可能異常的p維矢量,例如不同時期得到的結構前p階固有頻率矢量;{f}和Σ分別為實測的n個p維變量{f}的平均值矢量和協方差矩陣。T代表轉置運算。Di與某種給定的異常值閾值進行比較,大于閾值的即判斷為異常,而且從文獻[19]可以看出,Di越大,代表偏離正常情況越遠[20]。

如果只考慮各變量的分布而不考慮各變量之間的相互影響,馬氏距離可進一步簡化成為標準歐幾里德距離,其定義為:

其中,diag(Σ)表示協方差矩陣Σ的對角矩陣。標準化歐幾里德距離因為歸一化協方差矩陣為對角陣,計算量比馬氏平方距離要少很多,但卻保留了對異常值很好的檢測能力。

實際使用過程中,異常值如果是非先驗已知的,往往進行的是內插異常值檢測,即待檢測的測量矢量參與計算統計平均值,包括均值和方差等;如果有先驗已知的無損傷數據,則可以采用外推方式對可能異常的矢量進行異常值檢測。

2 螺栓松動定量檢測方法及其驗證

下面以單螺栓連接的懸臂梁為例,介紹螺栓松動程度的定量檢測方法。

2.1 試驗裝置

如圖1所示,試驗件由兩個金屬長梁搭接而成,搭接部分通過一個精制不銹鋼螺栓(M10(A2-70))貫穿連接;試驗件右端固定在基礎上,左端自由,成為一個懸臂梁;距試驗件左端10 mm的位置上、下表面分別設置測點和激勵點。

圖1 試驗件及測點設置示意圖Fig.1 Schematic of test structure

用激振器(SINOCERA JZK-5)激勵試件,并用加速度傳感器(PCBM352)在測點拾振,用LMSTest.Lab 12A模態測試儀及附帶軟件(Spectral Testing)測量系統的前5階彎曲固有頻率。測試現場見圖2。

圖2 測試現場照片Fig.2 Specie and test instruments

2.2 試驗過程

(1)將螺栓完全擰緊(擰緊力矩為24 N·m),采用隨機激勵方式,測得擰緊工況下懸臂梁的前5階固有頻率30組。這些數據用于計算初始狀態下固有頻率矢量的平均值矢量{f}和協方差矩陣的對角矩陣diag(Σ)。

(2)調整螺栓的擰緊力矩分別為24 N·m、15 N·m、5 N·m和手緊狀態,測量結構的前5階固有頻率各5組。其中,除了手緊狀態外,其余狀態均用力矩扳手調整。測試工況及對應測量次數見表1。這部分數據用于尋求固有頻率矢量的標準歐幾里德距離與擰緊力的關系。

表1 用于數據擬合的工況和測量次數Tab.1 Bolt conditions and No.of tests for curve fitting

(3)再次調整螺栓的預緊力分別為24 N·m、20 N·m、15 N·m、10 N·m、5 N·m和手緊狀態,測量結構的前5階固有頻率各5組。測試工況及對應測量次數見表2。這部分數據用于分析擰緊力和單一模態固有頻率之間的關系以及后期的試驗驗證。

表2 驗證工況和測量次數Tab.2 Bolt conditions and No.of tests for verification

2.3 SED和螺栓擰緊力的關系

設擰緊工況下測得的30組5維固有頻率矢量分別為{f1},{f2},…,{f30},計算這30組矢量的平均值矢量{?f}和協方差矩陣的對角矩陣diag(Σ),然后根據式(3),分別計算這30組矢量的SED值,結果見圖3。顯然,各樣本的SED值是隨機分布的。這種計算方法屬于內插法。

圖3 參考樣本的SED值Fig.3 SED values of initial state

將表1中24 N·m工況對應的5組測量結果分別代入式(3),得到5個SED值(這5個數據點是外推法計算的結果),將其與圖3的30個SED值繪制在一起,見圖4。從圖4可知,用外推法得到的這5個數據與前面30個內插法得到的數據混在一起時并不能將它們分離出來。這是合乎情理的,因為它們所代表的工況都是24 N·m的擰緊力矩。

圖4 內插法和外推法比較Fig.4 SED values of inclusive and exclusive

再用外推法求得表1中其它各工況對應的SED值,將結果與前期得到的數據一起繪制在圖5中。從圖5中可知:①4種擰緊力矩狀態對應的SED之間區分明顯,沒有混淆現象出現;②同一種擰緊力矩狀態對應的SED相互之間差別很小;③SED隨著擰緊力矩的減小單調地增大。

為了進一步尋找SED與擰緊力矩之間的量化關系,將表1各工況對應的20個SED值與擰緊力矩值用Matlab軟件進行曲線擬合,得到如圖6示的一個指數曲線。該曲線直觀地揭示了SED與擰緊力矩之間的量化關系,可以用作螺栓擰緊力矩的預測曲線,即螺栓松動損傷檢測的參考曲線。

圖5 4種工況對應的SED值Fig.5 SED values of table1

圖6 SED與擰緊力矩的關系曲線Fig.6 Curve of SED and tightening forces

2.4 利用SED評估螺栓松動程度的方法

根據以上分析,利用結構前5階固有頻率評估螺栓松動程度的方法歸納如下:

(1)測試擰緊狀態的固有頻率若干組,組成一個固有頻率矢量矩陣,計算該矩陣的平均矢量和協方差矩陣。其中樣本數必須不小于矢量的維數。

(2)測試有限幾種擰緊力工況的固有頻率矢量,用外推法分別計算各矢量到上述平均矢量之間的SED值,擬合各工況所得SED值與擰緊力之間的關系曲線,作為損傷程度檢測的參考曲線。

(3)測試待測狀態下結構的固有頻率矢量,計算其與上述平均矢量之間的SED值,然后將得到的結果與參考曲線進行比較,給出待測狀態的螺栓擰緊力估計值。

其中前兩步屬于基礎性的工作,需要在結構正式服役前進行。

2.5 試驗驗證

將表2中6種工況共30組測量數據依次代入式(3),用外推法求得對應的30個SED值,將其與緊固工況內插法得到的30個SED一起繪制在圖8中。比較圖7與圖5可知,圖7反映出與圖5相同的規律。再將表2中6種工況對應的30個SED值和圖6中的參考曲線繪制在一起(見圖8)。從圖8可知,用于驗證的6組共30個數據與參考曲線吻合良好(最大誤差出現在10 N·m的位置,絕對誤差約為4 N·m),從而證明了該方法的合理性。

圖7 驗證樣本的SED值Fig.7 SED values of table2

圖8 驗證數據點與參考曲線對比圖Fig.8 Curve and verification points

3 應用舉例

為了探討上述方法在多螺栓連接情況下的有效性,將其應用在某一衛星太陽帆板實物模型的螺栓松動定量化檢測過程之中。

3.1 試驗裝置

衛星太陽帆板通過夾具固定在基礎上(見圖9和圖10)。太陽帆板和夾具通過四個M6的螺栓連接,螺栓設計工作狀態的擰緊力矩為12 N·m。

圖9 試驗件及測點示意圖Fig.9 Schematic of test structure

圖10 測試現場照片Fig.10 Specie and test instruments

3.2 測試參考曲線

參考曲線的測試在太陽帆板模型的裝配過程中進行。首先,將四個螺栓用手均勻擰緊。然后,用力矩扳手將各螺栓的擰緊力矩調整到3N·m,調整方法參考文獻[21]中的方法進行,以保證各個螺栓的擰緊力矩相等,測試結構的固有頻率5組。再用相同的方法分別測試4 N·m、8 N·m工況下的固有頻率各5組,以及12 N·m工況下的固有頻率35組。最后利用前述方法計算各工況的SED值并擬合參考曲線(見圖11)。

圖11 SED與擰緊力矩的關系曲線Fig.11 Curve of SED and tightening forces

3.3 檢測效果

選用三種工況測試文中提出的方法對多螺栓結構的檢測效果。

工況1:四個螺栓都松至10 N·m;

工況2:四個螺栓都松至6 N·m;

工況3:兩側螺栓均松至4 N·m,中間兩個螺栓保持6 N·m。

在每種狀態測試結構的固有頻率各5組,利用式(3)計算SED指標,并與曲線對比,得出擰緊力矩評估結果(由于每種狀態有5組數據,所以可得到5個評估結果)。評估結果與實際擰緊力矩的對比見表3。

表3 擰緊力矩評估結果與實際值對比情況(單位:N·m)Tab.3 Comparison between test results and actual values(Unit:N·m)

從表3可知,工況1時,最大誤差為-0.78 N·m,如果利用5次測量的平均值,誤差為-0.19 N·m;工況2時最大誤差為0.96 N·m,如果利用5次測量的平均值,誤差為0.65 N·m;曲線雖不具備檢測工況3情況的能力,但工況3的平均在參考曲線上對應的預緊力矩值(4.91 N·m)介于6 N·m和4 N·m之間,體現了一種“等效損傷”程度。

4 結論

本文提出了一種利用結構前5階固有頻率的SED進行螺栓擰緊力定量評估的方法。通過上述研究,可以得出以下結論:

(1)結構前5階固有頻率的SED值與螺栓擰緊力矩呈指數關系;

(2)文中提出的方法不僅適用于單螺栓連接結構,對多螺栓連接結構也有一定的檢測能力。

此外,就相關問題討論如下:

(1)通過對太陽帆板模型的檢測結果可以發現,同一種狀態,用多次測量得到的固有頻率分別進行識別,其結果有一定的分散性,所以可取平均值作為最終的評估結果。

(2)該方法目前對各螺栓擰緊力矩相同情況的定量評估效果較好;另外,通過對太陽帆板工況3情況的檢測說明,對不同預緊力情況也有一定的“等效損傷”檢測能力;但對于更為復雜的情況和更高的檢測要求則需要進一步的研究,與其他局部檢測方法相結合是一個探索的方向。

(3)固有頻率反映的是結構的總體特性,而螺栓松動是結構的局部性能變化,要從固有頻率變化檢測螺栓的松動程度需滿一定的條件——局部性能變化影響到了結構的總體特性。從參考曲線可以看出,在螺栓松動程度不大時,對結構總體特性影響較小,所以曲線很平緩,此時該方法的檢測精度一般比較低;當螺栓松動到一定的程度,就會對結構總體特性有較大的影響,對應于曲線中相對陡峭的部分,此時該方法可以得到較準確的結果。

(4)如果不需要量化評估,只是簡單的趨勢監測,則無需參考曲線,以任何一個狀態作為起始狀態,用式(3)計算之后各時間節點測得的特征參數矢量的SED值,并進行基于時間的趨勢分析即可。這樣能夠拓展該方法的應用范圍。

[1]中國航空研究院.復合材料結構設計手冊[M].北京:航空工業出版社,2001.

[2]向志海,黃俊濤.螺栓松緊程度的受控敲擊檢測方法[J].試驗力學,2012,27(5):545-551.

XIANG Zhi-hai,HUANG Jun-tao.A controlled tap detection method for bolt tightness[J].Journal of Experimental Mechanics,2012,27(5):545-551.

[3]徐超,周幫友,劉信恩,等.機械螺栓連接狀態監測和辨識方法研究進展[J].強度與環境,2009,36(2):29-36.

XU Chao,ZHOU Bang-you,LIU Xin-en,et al.A review of vibration-based condition monitoring and identification for mechanical bolted joints[J].Structure&Environment Engineering,2009,36(2):29-36.

[4]Wang Tao,Song Gang-bing,Liu Shao-peng,etal.Review of bolted connection monitoring[J].International Journal of Distributed Sensor Networks,2013,871213(8pages).

[5]Yasui H,Kawashima K.Acoustoelastic measurement of bolt axial load with hypothetical velocity ratio method[J].Transactions of the Japan Society of Mechanical Engineers A,2000,66(642):390-396.

[6]Chaki S,Corneloup G,Lillamand I,et al.Combination of longitudinal and transverse ultrasonic waves for in situ control of the tightening of bolts[J].Journal of Pressure Vessel Technology,2007,129(3):383-390.

[7]Joshi SG,Pathare R G.Ultrasonic instrument formeasuring bolt stress[J].Ultrasonics,1984,22(6):270-274.

[8]Sun F P,Chaudhry Z,Liang C.Truss structure integrity identification using PZT sensor-actuator[J].Journal of Intelligent Material Systems and Structures,1995,6:134-139.

[9]王丹生,朱宏平,魯晶晶,等.基于壓電導納的鋼框架螺栓松動檢測試驗研究[J].振動與沖擊,2007,26(10):157-160.

WANG Dan-sheng,ZHU Hong-ping,LU Jing-jing,et al.Experimental study on detecting loosened bolts of a steel frame based on piezoelectric admittance[J].Journal of Vibration and Shock,2007,26(10):157-160.

[10]高峰,王德俊,江鐘偉,等.壓電阻抗技術用于螺栓松緊健康診斷[J].中國機械工程,2001,12(9):1048-1049.

GAO Feng,WANG De-jun,JIANG Zhong-wei,et al.Research on piezoelectric impedance technology for bolt health monitoring[J].China Mechanical Engineering,2001,12(9):1048-1049.

[11]Wang Tao,Song Gang-bing,Wang Zhi-gang.Proof-ofconcept study of monitoring bolt connection status using a piezoelectric based active sensing method[J].Smart Materials and Structures,2013,22:1-5

[12]Fasel TR,Todd M D,Park G.Active chaotic excitation for bolted jointmonitoring[C]//Proceeding of SPIE,Sensors and Smart Structures Technologies for Civil,Mechanical,and Aerospace System.Conference,San Diego,California,USA,2006,6174:61741S1-61741S8.

[13]楊智春,張慕宇,王樂.飛機加筋壁板緊固件松脫損傷檢測的自相關分析方法[J].振動與沖擊,2011,30(11):13-16.

YANG Zhi-chun,ZHANG Mu-yu,WANG Le.Damage detection of fastener loosening of aircraft stiffened panel by auto correlation analysis[J].Journal of Vibration and Shock,2011,30(11):13-16.

[14]楊智春,于哲峰.飛機壁板緊固件松脫損傷檢測的實驗研究[J].應用力學學報,2008,25(1):99-102.

YANG Zhi-chun,YU Zhe-feng.Experimental research on detection in fastener-loosing of aircraft panel by correlation function[J].Chiness Journal of Applied Mechanics,2008,25(1):99-102.

[15]Cawley P,Adams R D.The location of defects in structures from measurements of natural frequencies[J].Journal of Strain Analysis,1979,14(2):49-57.

[16]Liang R Y,Member,ASCE,et al.Theoretical study of crack-induced eigenfrequency changes on beam structures[J].Journal of Engineering Mechanics,1992,118(2):384-396.

[17]高維成,劉偉,鄒經湘.基于結構振動參數變化的損傷探測方法綜述[J].振動與沖擊,2004,23(4):1-7.

GAO Wei-cheng,LIU Wei,ZOU Jing-xiang.A review of damage detecting methods based on changes of structural vibrationparameters[J].Journal of Vibration and Shock,2004,23(4):1-7.

[18]Adams R D,Cawley P,Pye C J.A vibration technique for non-destructively assessing the integrity of structures[J].Journal of Mechanical Engineering Science,1978,20(2): 93-100.

[19]Worden K,Manson G.Damage detection using outlier analysis[J].Journal of Sound and Vibration,2000,229(3): 647-667.

[20]Nguyen T,Chan T H T,Thambiratnam D P.Controlled Monte Carlo data generation for statistical damage identification employing Mahalanobis squared distance[J].Structural Health Monitoring,2014,13(4):461-472.

[21]He K,Zhu W D.Detecting loosening of bolted connections in a pipeline using changes in natural frequencies[J].Journal of Vibration and Acoustics,2014,136:034503-1-034503-7.

[22]Ritdumrongkul S,Abe M,Fujino Y,et al.Quantitative health monitoring of bolted joints using a piezo ceramic actuator-sensor[J].Smart Materials and Structures,2004,13:20-29.

Bolt tightening force detection using outlier analysis of structural natural frequencies

GOU Bai-yong1,2,LU Qiu-hai1,WANG Bo1,WANG Shi-ying1
(1.School of Aerospace,Tsinghua University,Beijing 100084,China; 2.Aeronautic and Astronautic Engineering College,Air Force Engineering University,Xi'an 710038,China)

In order to evaluate the severity of bolt looseness quantitatively,a study was made to detect bolt tightening force using structural first 5 natural frequencies and the squared mahalanobis distance(MD)outlier analysis.Firstly,the structural natural frequencies under different tightening force condition weremeasured,and then,the standard euclidean distances(SEDs)between the natural frequency vector under fully-tightening condition and those under several not-fully-tightening conditionswere calculated.The fitting curve for SEDs and tightening forces wasmade using Matlab,and amethod was proposed to detect tightening force.Test verification and application examples were given at last.The results showed that there is an exponential relationship between SEDs and tightening forces,this relationship can be used to evaluate the severity of bolt looseness reliably.

bolt looseness;damage detection;structural health monitoring(SHM);natural frequency;outlier analysis;standard euclidean distance

TB123;TH113.1

A

10.13465/j.cnki.jvs.2015.23.014

國家自然科學基金資助項目(11090334,11372154)

2014-09-22修改稿收到日期:2014-12-03

緱百勇男,博士生,講師,1976年生

陸秋海男,博士,教授,1970年生

猜你喜歡
結構檢測方法
“不等式”檢測題
“一元一次不等式”檢測題
“一元一次不等式組”檢測題
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
論《日出》的結構
小波變換在PCB缺陷檢測中的應用
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
捕魚
主站蜘蛛池模板: 久久久黄色片| 特级做a爰片毛片免费69| www中文字幕在线观看| 国产欧美高清| 欧美自慰一级看片免费| 欧美黄网站免费观看| 无码内射在线| 99re这里只有国产中文精品国产精品| 国产91小视频| 国产精品人莉莉成在线播放| 毛片久久久| 亚洲中久无码永久在线观看软件| 激情成人综合网| 91热爆在线| 欧洲欧美人成免费全部视频| 2048国产精品原创综合在线| 国产欧美视频在线观看| 免费中文字幕一级毛片| 欧美成人h精品网站| 国产欧美视频综合二区| 毛片大全免费观看| 91综合色区亚洲熟妇p| 日韩午夜伦| 亚洲最黄视频| 国产区在线观看视频| 日韩福利在线视频| 亚洲视频一区在线| 玖玖精品视频在线观看| 久久久成年黄色视频| 亚洲成人高清无码| 国产成人综合在线观看| 欧美人与牲动交a欧美精品| 亚洲综合久久一本伊一区| 亚洲精品国产综合99| 亚洲无码高清免费视频亚洲| 五月综合色婷婷| 亚洲精品国产成人7777| 国产女人在线| 亚洲国产日韩视频观看| 国产精品久久久精品三级| 国产精品自拍露脸视频| 亚洲国产天堂久久综合226114| 91国内视频在线观看| 99r在线精品视频在线播放| 亚洲第七页| 国产一区二区三区免费| 日本影院一区| 国产精品久久自在自线观看| 国产白浆一区二区三区视频在线| 欧洲欧美人成免费全部视频| 四虎亚洲国产成人久久精品| 日韩资源站| 99在线观看国产| 国产黄视频网站| 99九九成人免费视频精品| 久久无码av一区二区三区| a在线亚洲男人的天堂试看| 久久天天躁狠狠躁夜夜躁| 91年精品国产福利线观看久久| 国产乱人视频免费观看| 国产精品久久久久久久久kt| 国产精选小视频在线观看| 成人日韩视频| 天堂岛国av无码免费无禁网站| 思思热精品在线8| 第一页亚洲| 尤物成AV人片在线观看| 欧美爱爱网| 国产高清不卡| 凹凸精品免费精品视频| 日韩欧美91| 好吊妞欧美视频免费| 尤物国产在线| 一区二区影院| 一本大道香蕉中文日本不卡高清二区| 91视频区| 亚洲国模精品一区| www.亚洲天堂| 国产丝袜无码一区二区视频| 国产精品性| 18禁黄无遮挡网站| 亚洲第一区欧美国产综合|