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

基于頻響主元的衛星結構損傷定位與評估方法

2016-03-16 07:09:40馮彥軍石川千周宇趙發剛周徐斌
航天器工程 2016年4期
關鍵詞:區域結構

馮彥軍 石川千 周宇 趙發剛 周徐斌

(上海衛星工程研究所,上海 201109)

基于頻響主元的衛星結構損傷定位與評估方法

馮彥軍 石川千 周宇 趙發剛 周徐斌

(上海衛星工程研究所,上海 201109)

鑒于目前國內衛星預、復振頻響結果的評估主要依靠設計師的經驗分析,缺乏高效準確的量化手段,文章基于頻響函數(FRF)的主元分析(PCA),提出了兩種針對衛星結構的損傷定位指標,分別是多主元綜合指標和FRF主元殘差均方根;并采用平方預測誤差(SPE)指數對衛星結構的整體損傷程度進行評估。以衛星主承力部件中心承力筒為分析對象,驗證了基于衛星結構頻響主元的損傷定位和評估方法的有效性。文章對基于預、復振對比分析的損傷識別、量化試驗評估標準、提高衛星研制效率有促進作用。

損傷識別;多主元綜合指標;頻響函數主元殘差均方根;平方預測誤差指數;衛星預、復振;頻響函數

1 引言

動力學環境試驗是對衛星等航天產品設計和工藝狀態適用性的考核,是研制過程中不可缺少的環節[1]。在每個滿量級的動力學環境試驗前后,應分別進行預、復振試驗檢驗,通過對比預、復振頻響曲線獲得傳遞特性、共振頻率和各階幅值放大系數等結構特性,識別結構在振動試驗過程中可能發生的結構損傷[2]。大量的振動試驗結果表明:預、復振試驗響應曲線不一致現象非常普遍,典型現象有共振峰漂移、共振峰數量變化和共振峰幅值變化等,究其機理,主要受邊界非線性、材料非線性和結構損傷等因素的影響[3-9]。目前對于預、復振試驗頻響結果不一致現象的分析,主要依靠設計師的工程經驗,檢查手段局限于對試驗件的拆解和復查故障,準確性不高,效率低下,大大制約了衛星的研制進度。因此,基于衛星預、復振頻響數據的損傷識別是極具應用價值并亟待解決的問題。

將損傷識別技術應用于衛星預、復振不一致問題的分析,可以有效地利用試驗數據信息,并找到故障所在。基于結構動力特性的損傷識別是近幾十年國內研究的熱點。目前大多數損傷識別方法是通過頻響函數(Frequency Response Function,FRF)數據進行模態擬合,提取模態參數作為損傷標識量,但存在測量噪聲、模態提取誤差、數據冗余和分析難度大等問題,因此一些學者直接利用FRF中豐富的信息進行損傷識別[10-12],或用主元分析(Principal Component Analysis,PCA)方法對FRF數據進行壓縮、提取和處理[13-16],得到包含損傷信息的特征數據,以此構造損傷識別指標。其中,損傷定位指標的優劣嚴重影響損傷識別的精度:為了避免構造不適當的損傷定位指標,一些學者[14-16]將PCA壓縮的FRF數據直接或間接輸入到神經網絡,進行損傷模式識別,但需要大量有效的實測數據來訓練神經網絡,而且識別精度受測量噪聲和模型誤差的影響很大。楊彥芳等人[17]采用橢圓控制圖和T2控制圖實現對網架結構的損傷定位,但網架損傷形式單一,且不能進行多損傷識別和損傷程度評估。相比于FRF,FRF曲率對損傷更敏感,Maia[18]、姜增國[19]和朱新圓[20]等人利用頻響函數曲率及其變形構造損傷識別指標,對傳力形式相對單一的梁、板和桁架結構進行了損傷識別。但現有的基于FRF的損傷識別理論,僅應用于梁桿件、平板、桁架等傳力形式相對單一的特殊結構,而衛星結構形式多樣,傳力路徑復雜,主承力結構以復合材料殼為主,單純利用以上現有方法不能很好地識別局部損傷位置,更不能對損傷程度進行評估。

針對上述問題,本文提出了兩種線性歸一化的無量綱化指標:多主元綜合指標和FRF主元殘差均方根,分別對衛星結構FRF數據進行分析,實現結構損傷定位,并用平方預測誤差(SquarePredictionError,SPE)指數對結構整體損傷程度進行了評估。通過衛星主承力筒模型在不同損傷工況下的數值仿真分析,對上述損傷定位和程度評估方法進行了驗證。

2 基于頻響主元分析的損傷識別理論

2.1 PCA方法

PCA分析步驟為:

(1)對H(ω)按列進行標準化處理,將各變量化為均值為0,方差為1的標準變量。

(2)對H(ω)的協方差矩陣S進行特征值分解:

(1)

2.2 損傷定位理論

首先對損傷識別矩陣H(ω)進行主元分析,得到包含損傷信息的主元Fn×k,然后基于其主元構造出兩種損傷定位指標:多主元綜合指標和FRF主元殘差均方根。

2.2.1 多主元綜合指標

由多元控制理論知,當測點數n足夠大時,由任意測點j前2階主元決定的坐標(yj1,yj2)服從式(2)的橢圓控制域分布。第j點前2階主成分對應的點距離橢圓控制邊界的距離,一定程度上代表了該點數據的異常程度,因此定義偏離橢圓邊界的相對距離為dj1,見式(3)。

(2)

(3)

(4)

(5)

則第j點的總體相對距離為dj=dj1+dj2。各測點的dj構成向量d,經0~1線性歸一化,即得到多主元綜合指標:

(6)

控制限(ControlLimit,CL)的確定。橢圓控制限和T2控制限對應的相對距離dCL=0,將dCL與d一起歸一化即可得到標準化CL,見式(7),式中L1為CL的值。d中大于L1的元素對應的測點被認為是異常點,即位于損傷區域附近的測點。

(7)

2.2.2FRF主元殘差均方根

對損傷識別矩陣H(ω)的前k階主成分矩陣Fn×k,定義第j個測點對應的FRF主元殘差均方根為rj(見式(8)),各測點構成的FRF主元殘差均方根向量為R,經0~1線性歸一化后,可得到標準化的FRF主元殘差均方根向量R,見式(9)。R中大于CL值的元素對應的測點被認為是異常點。

(8)

(9)

CL的確定,可以通過兩種方法確定:一種是通過統計的方法,統計學中常用3σ準則來排除數據中偏差較大的數據,本文中用L2=E(R)+2D(R)確定CL值,E(R)為R的均值,D(R)為R的標準差;另一種方法,考慮到具體實際試驗中由于噪聲、激勵能量不穩定及其它干擾因素,即使對于同樣工況的同一試件,其實測頻響數據的主元殘差均方根一般不為零。因此,以完好結構相同條件下的連續兩次實測頻響數據的主元殘差均方根向量中最大元素為CL。若某工況下某測點對應的FRF主元殘差均方根大于CL,則被認為是異常點。

2.3 損傷程度評估理論

(10)

式中:I為單位矩陣;SPE指數定義為Hd(ω)在RS投影的方差,表征了有損結構實測數據對無損結構主元模型的偏離程度,見式(11),E為SPE指數的值。因此,當結構無損傷(或損傷程度較小)時,Hd(ω)的主要信息落在主元子空間(PrincipleComponentSubspace,PCS),SPE指數很小,主要為噪聲;當結構有損傷(或損傷程度較大)時,Hd(ω)在RS的投影就會增加,SPE指數也增加,主要為結構損傷引起的異常信息。因而可通過SPE指數來評估結構的整體損傷程度。

(11)

3 衛星承力筒的頻響主元分析算例

3.1 算例模型

為了驗證上述提出的損傷識別方法,采用衛星典型主承力部件承力筒進行數值驗證。如圖1所示,承力筒筒體柱段外徑Ф1114mm,內徑Ф1090mm,高3260mm,主體采用碳纖維/環氧-鋁蜂窩夾層結構復合材料,橫向加強框和縱向桁條均為碳纖維材料。采用Pro/E軟件進行三維實體建模,Patran軟件進行有限元建模。筒體部分采用殼單元,利用三明治夾芯等效理論模擬等效蜂窩夾層復合材料;面板安裝法蘭和縱向桁條采用梁單元。筒體部分單元總數為10 000,周向和縱向各劃分100個網格節點,編號從星箭連接環端至上端加強框依次增大。

圖1 衛星承力筒模型Fig.1 Satellite bearing cylinder model

通常認為結構損傷并不引起質量的損失,而只引起剛度損失。因此,本文采用剛度折減法模擬損傷,即對損傷部位的材料彈性模量折減來模擬局部損傷。設計9個工況,工況1為完好結構,用于和有損結構形成對比;工況2~5用來驗證損傷定位方法;工況3及工況6~9用于驗證損傷程度評估方法,見表1。每種工況下損傷區域均為沿軸向5行、周向25列的125個剛度折減單元的連續區域,如圖2。承力筒從下至上,每周5個測點,共11周,最下面一周測點編號為1~5,編號依次增大,最上面一周測點編號為51~55,在損傷區域附近布置4個測點,一個測點位于損傷區域中心,其它三個測點位于損傷區域邊緣,圖3是承力筒的平面展開圖。各工況對應的損傷區域附近測點編號見表2。

表 1 工況設計Table 1 Design of working condition

注:h表示損傷區域距離筒體下端的距離占筒體高的百分比,損傷程度s為剛度折減率

圖2 損傷區域Fig.2 Damage region

圖3 測點和損傷區域位置(沿承力筒母線展開)Fig.3 Measuring point and damage region location (unfold drawing along cylinder generatrix line)

工況損傷區域附近測點編號21,2,6,73、6、7、8、916,17,21,22431,32,36,37546,47,51,52

3.2 損傷定位驗證

仿真驗證的總體流程為:采用MSC/Nastran軟件對承力筒進行工況1~9下20~2000 Hz的頻響分析,加速度激勵施加在星箭連接環基礎上,利用PCL提取55個測點的FRF數據。利用Matlab軟件編制PCA程序,提取FRF主元和對應的貢獻率,構造出兩種損傷定位指標,驗證算法可行性。

3.2.1 構造損傷識別矩陣

以工況1(完好結構)和不同損傷工況2~5(損傷結構)下各測點的FRF曲線幅值差為變量,以每一測點得到FRF數據為樣本,組成55行(測點數)900列(譜線數)組成損傷識別矩陣H(ω)。

3.2.2 求各階主元及相應的貢獻率

利用工況2~5下對應的結構損傷識別矩陣,求其協方差陣的900個特征值和特征向量,提取前10階主元,對應的貢獻率如圖4和表3。

圖4 前10階主元貢獻率Fig.4 First ten principle components contribution rate

工況2345前10階累積貢獻率0.90400.91060.89660.9292

從表3和圖4可以看出,各工況下的前10階主元累積貢獻率大于89%,明顯大于其它各階成分,說明前10階主元包含了損傷識別矩陣900個變量的絕大部分信息。因而僅利用前10階主元所包含的數據信息進行分析損傷識別是可行的。

3.2.3 兩種損傷定位指標識別結果及分析

利用多主元綜合指標對工況2~5四種損傷狀態進行識別,置信度為95%。識別結果見圖5。

圖5 多主元綜合指標仿真結果Fig.5 Simulation of multi-principal component comprehensive index

歸一化的FRF主元殘差均方根法,表征了各測點間的相對損傷程度。損傷程度從0~1,0代表無損傷,1代表損傷程度最嚴重。根據2σ準則,各工況相應的CL如圖6所示。

圖6 FRF主元殘差均方根指標損傷識別結果Fig.6 Simulation of FRF principle component residual mean square root

該兩種損傷定位指標識別得到的異常點,即落在CL以外的測點編號,匯總對比后,如表4所示。

表4 兩種損傷定位指標識別結果對比Table 4 Contrast of identification result between the two damage location indexes

對以上分析得到以下幾點結果:

(1)多主元綜合指標識別結果存在錯判(測點11、41),而FRF主元殘差均方根無錯判。4個工況下,實際損傷區域測點為16個,前者識別到13個異常點,有2個錯判,準確率達84.6%;而后者識別到11個異常點,無錯判,準確率達100%。兩種指標基本能準確識別到損傷位置。測點11、41之所以被錯判,可能是該測點距離損傷區域比較近(測點11位于工況2損傷區域軸向正上方,測點41位于工況5損傷區域軸向正下方),軸向傳力特性比較明顯,且遠離節點,響應較大,因此該測點受損傷區域的影響較大,表現為異常點。

(2)多主元綜合指標沒有識別出測點1、16、31、46,存在漏判。原因可能是測點1、16、31、46位于損傷區域中心,其它測點位于損傷區域和完好區域的過渡區,結構特性發生突變,更容易被檢測出來。而利用FRF主元殘差均方根指標,工況3、4下能夠識別到損傷區域中心的測點。因此,綜合參考兩種損傷識別指標,在識別到的異常點附近進行檢查校核,能夠有效防止漏判情況。

(3)FRF殘差均方根方法在工況2漏判情況較其它工況嚴重,漏判率達50%。原因可能是,承力筒約束方式為固支,1、2測點位于底端,響應較小,信噪比太低,不容易被識別出來。故在實際測試中,對于遠離約束端和節點位置的響應較大的測點,更容易被識別出來。

(4)損傷定位可以采用“先軸向、再周向”原則。從圖5和圖6中可以看出,從軸向和周向兩個維度看,越靠近損傷區域,定位指標數值越大。如圖6(b)中,測點6、11、16、21、26、31位于損傷區域軸向正上(下)方,比同一高度其它測點指標更大;而16、17、18、19、20位于損傷區域周向左右,比同一周向其它測點指標更大。位于承力筒上的損傷區域附近或與損傷區域同周向角度、同軸向高度的測點,兩種損傷定位指標較高,更有利于被識別出來。

3.3 損傷程度評估

工況3和6~9代表結構損傷區域不同的材料剛度折減率,用于驗證損傷程度評估算法,求得對應的SPE指數,如表5和圖7。

表5 不同損傷程度下SPE指數Table 5 SPE indexes under different damage degrees

圖7 不同損傷程度下SPE變化趨勢Fig.7 SPE tendency at various damage degree

從圖7中和表5中可以看出,工況1下(無損結構)頻響數據的對應的SPE指數為3 219.1,隨著剛度折減率逐漸增加,結構損傷程度逐漸加重,從10%增加到80%,對應的SPE指數也逐漸增加,從3 340.4逐漸到10 474.0。因此,可以將工況1對應的SPE指數設為閾值,大于該值說明損傷發生,偏離該值越大,說明結構整體損傷程度越嚴重。由此可以看出,SPE指數既可以作為結構整體損傷程度的評估,又可以作為結構損傷開始的預示。

4 結論

本文以航天器振動試驗中預、復振頻響不一致問題為背景,將主元分析和損傷識別理論應用于預、復振頻響數據的分析,提出了兩種基于頻響主元的損傷定位指標,分別是多主元綜合指標和FRF主元殘差均方根指標,并用SPE指數對結構整體損傷程度進行了評估。通過衛星典型主承力部件承力筒進行了數值仿真,驗證了上述損傷識別方法的有效性,得到了以下結論:

(1) 測點利用頻響數據構造損傷識別矩陣,其前10階主元累積貢獻率達到了90%左右,包含了FRF數據矩陣的大部分信息,因此利用前10階主元進行損傷識別是可行的。

(2) 兩種損傷識別指標均能基本準確識別出損傷區域的大致位置。多主元綜合指標可以識別到損傷區域邊緣及其附近遠離節點、響應較大的測點,但無法識別到損傷區域內部的測點;FRF主元殘差均方根指標可以準確識別到損傷區域邊緣,并能夠識別到部分損傷區域內部的測點。實際識別時,應該綜合考慮兩種指標,減小錯判和漏判的可能,提高損傷定位效率。

(3)位于承力筒上的損傷區域附近或與損傷區域同周向角度、同軸向高度的測點,兩種損傷定位指標較高,更利于被識別出來。因此實際識別時,可以采取“先軸向、后周向”的損傷定位原則,先定位出損傷區域所在的軸向高度,再定位出損傷區域的周向角度。

(4) SPE指數可以表征結構的整體損傷程度。隨著損傷程度的加深,SPE指數逐漸增大。以無損結構得到的SPE指數為CL,既可以預示結構發生損傷,也可以評估結構的整體損傷程度。

本文提出的方法對衛星部件(承力筒)的損傷定位和評估有效,是否適用于更復雜的衛星系統還有待進一步驗證。該方法對基于預、復振對比分析的損傷識別、量化試驗評估標準及提高衛星研制效率有促進作用。

References)

[1] 柯受全,金恂叔.衛星環境工程和模擬試驗[M].中國宇航出版社,1996

Ke Shouquan,Jin Xunshu. Satellite environmental engineering and simulation test[M]. Beijing: China Astronautics Press,1996 (in Chinese)

[2]向樹紅. 航天器力學環境試驗技術[M]. 北京: 中國科學技術出版社,2008: 18-58

Xiang Shuhong. Spacecraft mechanics environment test technology[M]. Beijing: China Science and Technology Press,2008:18-58 (in Chinese)

[3]Nobukatsu O kuizumi,M C Natori. Nonlinear vibrations of a satellite truss structure with gaps[C]//45thAIAA/ASME/ASCE/AHS/ASC Structures,Structural Dynamics & Materials Confer. Washington D.C.: AIAA,2004

[4]衛洪濤,孔憲仁,王本利,等. 非線性連接結構對一個典型衛星頻率漂移的影響[J]. 航天器環境工程,2012,29(3): 297-302

Wei Hongtao,Kong Xianren,Wang Benli,et al. Research about the influence of nonlinear joints on a typical satellite frequency shift[J]. Spacecraft Environment Engineering,2012,29(3): 297-302 (in Chinese)

[5]陳昌亞,宋漢文,王德禹,等.衛星振動頻率漂移現象與非線性參數識別[J]. 上海交通大學學,2005,39(7): 1197-1200

Chen Changya,Kong Xianren,Wang Deyu,et al. The natural frequency shift of satellite vibration test and parameter identification of nonlinear in satellite structure[J]. Journal of Shanghai Jiaotong University,2005,39(7): 1197-1200 (in Chinese)

[6]吉桂秀. 含多分層復合材料層合板和加筋板的自由振動特性研究[D]. 大連: 大連理工大學,2007

Ji Guixiu. Analysis of free vibration characteristics for multi-delaminated composite plates and stiffened plates[D].Dalian: Dalian University of Technology,2007 (in Chinese)

[7]Anderson J R,Van Noord J L,Soulas G C. Environmental testing of the NEXT PM1 ion engine[R]. Washington D.C.: AIAA,2007

[8]Dervan J,Robertson B,Staab L,et al. Development testing and subsequent failure investigation of a spring strut mechanism[C]// Proceedings of the 42nd Aerospace Mechanisms Symposium. Greenbelt: NASA Goddard Space Flight Center,2014: 1-14

[9]王東升,周桐,張志旭. 淺析振動試驗中結構非線性引發的現象[J]. 航天器環境工程,2007,24(5): 314-317

Wang Dongsheng,Zhou Tong,Zhang Zhixu. A simple phenomenon analysis for structure nonlinearity in vibration tests[J]. Spacecraft Environment Engineering,2007,24(5): 314-317 (in Chinese)

[10] Carden E P. Vibration based condition monitoring: a review[J]. Structural Health Monitoring,2004,3(4):355-377

[11]Lee U,Shin J. A frequency response function-based structural damage identification method[J]. Computers & Structures,2002,80(2):117-13

[12]Maia N M M,Silva J M M,Almas E A M,et al. Damage detection in structures: from mode shape to frequency response functionmethods[J]. Mechanical Systems & Signal Processing,2003,17(3):489-498

[13]Zang C,Imregun M. Structural damage detection using artificial neural networks and measured FRF data reduced via principal component projection[J]. Journal of Sound and Vibration,2001,242(5): 813-827

[14]Ni Y Q,Zhou X T,Ko J M. Experimental investigation of seismic damage identification using PCA-compressed frequency response functions and neural networks[J]. Journal of Sound & Vibration,2006,290(1-2):242-263

[15]Dackermann U,Li J,Samali B,et al. Structural damage identification utilising PCA-compressed frequency response functions and neural network ensembles[C]//Proceedings of 20th Australasian Conference on the Mechanics of Structures and Materials. Toowoomba: Faculty of Engineering and Surveying,2008:803-809

[16]Nguyen V V,Dackermann U,Li J,et al. Damage identification of a concrete arch beam based on frequency response functions and artificial neural networks[J]. Electronic Journal of Structural Engineering,2015,14(1):75-84

[17]楊彥芳,宋玉普,紀衛紅. 基于實測頻響函數主成分的在役網架損傷識別方法[J]. 振動與沖擊,2007,26(9): 128-132

Yang Yanfang,Song Yupu,Ji Weihong. A damage identification method for applied truss structure based on measured frequency response functions principle component[J]. Journal of Vibration and Shock,2007,26(9): 128-132 (in Chinese)

[18]Maia N M M,Silva J M M,Sampaio R P C. Localization of damage using curvature of the frequency res-ponse functions[J]. Proceedings of SPIE,1997,3089: 942-946

[19]姜增國,張楨. 基于應變頻響函曲率的結構損傷識別[J]. 建筑科學與工程學報,2009,26(4): 40-43

Jiang Zengguo,Zhang Zhen. Structural damage identification based on SFRF curvature[J]. Journal of Architecture and Civil Engineering,2009,26(4): 40-43 (in Chinese)

[20]朱新圓,阿肯江·托呼提. 鋼桁架結構損傷識別的頻響函數曲率法[J]. 實驗室研究與探索,2014,33(6): 40-44

Zhu Xinyuan,Akenjiang·Tuohuti. A damage identification study of steel truss structure based on frequency response function[J]. Research and Exploration in Laboratory,2014,33(6): 40-44 (in Chinese)

[21]Johnson R A,Wichern D W.實用多元統計分析[M]. 北京:清華大學出版社,2001

Johnson R A,Wichern D W. Practical multivariate statistical analysis[M]. Beijing: Tsinghua University Press,2001 (in Chinese)

(編輯:張小琳)

Methods of Damage Location and Evaluation for Satellite Structure Based on Frequency Response Principle Component

FENG Yanjun SHI Chuanqian ZHOU Yu ZHAO Fagang ZHOU Xubin

(Shanghai Institute of Satellite Engineering,Shanghai 201109,China)

Considering that current evaluation methods for satellites’ pre and post vibration frequency response mainly depend on designers’ experience,lacking efficient quantization methods,this paper first puts forward two damage location indexes for satellite structure based on principle component analysis (PCA) of frequency response function (FRF). The two indexes are mulit-principal component comprehensive index and FRF principle component residual mean square root. Then,square prediction error (SPE) index is adopted to evaluate the whole damage degree for satellite structure. Finally,by taking satellite bearing cylinder as analysis objective, the efficiency of the damage location and evaluation methods based on PCA is validated. This paper brings positive effect on damage identification based on pre and post vibration analysis,quantization of experiment evaluation standard and improvement of satellite manufacture.

damage identification; multi-principle component comprehensive index; FRF principle component residual mean square root; SPE index; satellite pre and post vibration; frequency res-ponse function

2016-06-27;

2016-07-12

國家重大航天工程

馮彥軍,男,碩士研究生,從事衛星結構設計和振動測試工作。Email:feng_yan_jun@126.com。

V416.2;V414.6

A

10.3969/j.issn.1673-8748.2016.04.018

猜你喜歡
區域結構
永久基本農田集中區域“禁廢”
今日農業(2021年9期)2021-11-26 07:41:24
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
分割區域
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
新型平衡塊結構的應用
模具制造(2019年3期)2019-06-06 02:10:54
論《日出》的結構
關于四色猜想
分區域
基于嚴重區域的多PCC點暫降頻次估計
電測與儀表(2015年5期)2015-04-09 11:30:52
創新治理結構促進中小企業持續成長
現代企業(2015年9期)2015-02-28 18:56:50
主站蜘蛛池模板: 国产剧情一区二区| 特级aaaaaaaaa毛片免费视频| 亚洲欧美日韩天堂| 国产精品免费入口视频| 午夜影院a级片| 日韩毛片在线播放| 亚洲一区第一页| 一区二区偷拍美女撒尿视频| 女同久久精品国产99国| 99这里精品| 国产精品无码一二三视频| 四虎成人免费毛片| 97在线免费| 国产高潮视频在线观看| 婷婷丁香色| 又爽又大又黄a级毛片在线视频| 看国产一级毛片| 88国产经典欧美一区二区三区| 最新国产精品第1页| 国产成人一二三| 中文字幕佐山爱一区二区免费| 久久综合久久鬼| 国产丝袜啪啪| 国产97视频在线| 99re视频在线| 91精品国产自产在线观看| 欧美人与牲动交a欧美精品| 欧美国产视频| 成人噜噜噜视频在线观看| 刘亦菲一区二区在线观看| 国产成人在线无码免费视频| 波多野结衣一区二区三视频 | 永久免费无码成人网站| 毛片a级毛片免费观看免下载| 亚洲精品男人天堂| 国产精品网拍在线| 51国产偷自视频区视频手机观看| 色综合久久久久8天国| 综合色在线| 亚洲人成网站在线观看播放不卡| 美女无遮挡拍拍拍免费视频| 国产在线观看精品| Jizz国产色系免费| 国产乱人免费视频| 午夜久久影院| 亚洲av色吊丝无码| 国产精品无码久久久久久| 亚洲爱婷婷色69堂| 另类综合视频| 国产精品无码制服丝袜| 国产成人精品第一区二区| 久久久久亚洲AV成人网站软件| 午夜精品久久久久久久99热下载| 免费观看亚洲人成网站| 中文字幕亚洲乱码熟女1区2区| 98超碰在线观看| 国产成人夜色91| 97久久超碰极品视觉盛宴| 亚洲AⅤ无码国产精品| 国产高清国内精品福利| 亚洲视频a| 国产午夜无码专区喷水| 啪啪啪亚洲无码| 亚洲天堂精品在线观看| 久久精品国产亚洲AV忘忧草18| 亚洲天堂免费观看| 亚洲码一区二区三区| 久久无码高潮喷水| 一本色道久久88| 国产精品久久久久久久伊一| 国产一区二区视频在线| 国产成人区在线观看视频| Jizz国产色系免费| 亚洲视屏在线观看| 免费看av在线网站网址| 日本高清成本人视频一区| 日韩欧美中文字幕一本| 日本伊人色综合网| 欧美日韩亚洲国产主播第一区| 麻豆国产在线不卡一区二区| 情侣午夜国产在线一区无码| 国产h视频免费观看|