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

        ?

        基于無(wú)網(wǎng)格法梁結(jié)構(gòu)多裂紋定量識(shí)別

        2016-08-04 07:09:28羅志鋼蔣占四王衍學(xué)向家偉
        振動(dòng)與沖擊 2016年13期
        關(guān)鍵詞:網(wǎng)格法曲率固有頻率

        羅志鋼, 蔣占四, 王衍學(xué), 向家偉

        (1. 桂林電子科技大學(xué) 機(jī)電工程學(xué)院,廣西 桂林 541004; 2.溫州大學(xué) 機(jī)電工程學(xué)院,浙江 溫州 325035)

        基于無(wú)網(wǎng)格法梁結(jié)構(gòu)多裂紋定量識(shí)別

        羅志鋼1, 蔣占四1, 王衍學(xué)1, 向家偉2

        (1. 桂林電子科技大學(xué) 機(jī)電工程學(xué)院,廣西 桂林541004; 2.溫州大學(xué) 機(jī)電工程學(xué)院,浙江 溫州325035)

        研究基于無(wú)網(wǎng)格法的懸臂梁結(jié)構(gòu)多裂紋定量識(shí)別方法。通過(guò)局部加權(quán)殘量法將裂紋以虛擬邊界融入離散系統(tǒng)方程,建立多裂紋梁結(jié)構(gòu)動(dòng)力學(xué)模型,獲得定量識(shí)別正問(wèn)題數(shù)據(jù)庫(kù)。結(jié)合曲率模態(tài)法先找到裂紋位置,再以實(shí)測(cè)頻率作為反問(wèn)題輸入利用粒子群法反演尋優(yōu)從而定量預(yù)測(cè)出裂紋的損傷程度。數(shù)值仿真表明,該方法前處理簡(jiǎn)單,適合用來(lái)模擬非連續(xù)裂紋問(wèn)題,為結(jié)構(gòu)多裂紋識(shí)別提供了新途徑。

        無(wú)網(wǎng)格法;多裂紋;識(shí)別

        結(jié)構(gòu)的裂紋是一種常見(jiàn)的損傷形式,在長(zhǎng)期交變應(yīng)力和沖擊載荷的作用下,容易引起裂紋擴(kuò)展從而導(dǎo)致結(jié)構(gòu)破壞,引起安全事故。因而對(duì)結(jié)構(gòu)是否產(chǎn)生裂紋并對(duì)其位置和損傷程度進(jìn)行及時(shí)診斷和評(píng)估[1]在工程實(shí)踐中具有重要意義。

        近年來(lái)國(guó)內(nèi)外對(duì)含裂紋的結(jié)構(gòu)進(jìn)行了大量分析研究,比如基于小波包與倒頻譜分析的裂紋診斷方法[2], 基于網(wǎng)格式支持向量機(jī)算法的轉(zhuǎn)軸裂紋故障診斷[3]。這些方法在診斷裂紋故障上取得了較好的效果,但準(zhǔn)確性有待提高,尤其多裂紋問(wèn)題存在困難。在用基于模型的方法分析裂紋方面以有限元法最為普遍。陳雪峰等[4-6]將小波有限元引入到裂紋故障診斷。利用小波多尺度多分辨的特性,針對(duì)求解問(wèn)題的精度要求,可采用不同類型和尺度的小波基函數(shù),因而在處理結(jié)構(gòu)裂紋奇異性問(wèn)題具有一定的優(yōu)越性。文獻(xiàn)[5]在懸臂梁裂紋識(shí)別中利用6階Daubechies小波基單元,取得了良好的計(jì)算結(jié)果。但是Daubechies小波無(wú)顯式表達(dá)式,用于構(gòu)造小波單元時(shí)必須計(jì)算聯(lián)系系數(shù),計(jì)算復(fù)雜。

        用有限元法模擬裂紋問(wèn)題時(shí),由于裂紋尖端的應(yīng)力奇異性和裂紋的動(dòng)態(tài)擴(kuò)展往往需要對(duì)網(wǎng)格細(xì)化或網(wǎng)格的重新劃分,特別在多裂紋問(wèn)題求解中,網(wǎng)格的大量重構(gòu)不僅降低了模型的計(jì)算精度,而且增加了計(jì)算成本。

        無(wú)網(wǎng)格法解決了傳統(tǒng)基于網(wǎng)格的數(shù)值方法缺點(diǎn),它不存在類似有限元中單元邊界要始終和結(jié)構(gòu)內(nèi)部不連續(xù)邊界(裂紋)保持一致的問(wèn)題,避免了網(wǎng)格重構(gòu)帶來(lái)的數(shù)值困難。本文通過(guò)無(wú)網(wǎng)格PETROV-GALERKIN法[7]對(duì)多裂紋懸臂梁結(jié)構(gòu)進(jìn)行研究。從加權(quán)殘量法出發(fā)將裂紋以虛擬邊界引入,建立多裂紋結(jié)構(gòu)的正問(wèn)題動(dòng)力學(xué)模型;結(jié)合前人研究的曲率模態(tài)法[8-11]確定裂紋位置的基礎(chǔ)之上,將反問(wèn)題實(shí)測(cè)系統(tǒng)的前三階固有頻率作為輸入,利用粒子群優(yōu)化(PSO)方法反演尋優(yōu)定量預(yù)示該位置的損傷程度。通過(guò)懸臂梁雙裂紋診斷試驗(yàn),驗(yàn)證該方法的有效性。

        1 MLS形函數(shù)

        根據(jù)移動(dòng)最小二乘法,在域Ω中任一點(diǎn)的近似表達(dá)式為:

        (1)

        式中:φi為對(duì)應(yīng)于點(diǎn)x的支持域中n個(gè)節(jié)點(diǎn)的MLS形函數(shù);節(jié)點(diǎn)i的形函數(shù)φi(x)為:

        (2)

        A(x)為加權(quán)力矩矩陣,定義為:

        (3)

        矩陣B(x)為:

        (4)

        2 多裂紋懸臂梁的無(wú)網(wǎng)格分析

        圖1,2分別多裂紋懸臂梁及其無(wú)網(wǎng)格模型。其中L是梁的長(zhǎng),L1,L2…Li;a1,a2…ai分別為相應(yīng)裂紋對(duì)應(yīng)的位置和深度。假定均勻截面的梁的自由振動(dòng)方程為:

        0

        (5)

        式中:u(x,t)是梁的軸向橫位移,ρ是梁的密度,A是梁的截面面積,EI是梁的彎曲剛度,將u(x,t)空間坐標(biāo)進(jìn)行離散化:

        u(x,t)=u(x)eiwt=ueiwt

        (6)

        將式(6)代入等式(5)化簡(jiǎn)為:

        EIu″″=ω2ρAu0

        (7)

        這樣u=u(x)變成只是空間坐標(biāo)的函數(shù),ω是自然頻率。

        圖1 多裂紋懸臂梁模型Fig.1 Multiple cracks cantilever model

        圖2 多裂紋懸臂梁無(wú)網(wǎng)格模型Fig.2 Multiple cracks cantilever meshless model

        邊界條件為:

        u(x)=u0在Γu上,

        (8)

        M=M0在ΓM上,

        V=V0在ΓV上

        (9)

        Γu∩ΓV=?,Γθ∩ΓM=?

        (10)

        其中:M是彎力矩,V是剪切力,Γu,Γθ,ΓM,ΓV,是邊界上的位移,轉(zhuǎn)角,彎力矩,剪切力。在裂紋處,連續(xù)條件為:

        (11)

        (12)

        (13)

        不連續(xù)條件為轉(zhuǎn)角[12]:

        (14)

        式中:xc為裂紋處的坐標(biāo),λ表示裂紋所在部位梁等效為無(wú)質(zhì)量扭轉(zhuǎn)彈簧的影響系數(shù),λ能被表示成裂紋深度a和梁高h(yuǎn)的函數(shù)[13-14]:

        λ=5.346hg(ξ)

        (15)

        其中:

        g(ξ)=1.862 4ξ2-3.95ξ3+16.375ξ4-

        37.226ξ5+76.81ξ6-126.9ξ7+

        172ξ8-143.97ξ9+66.56ξ10

        (16)

        (17)

        采用局部加權(quán)殘量法[7,15]:

        (18)

        (19)

        (20)

        (21)

        由力矩M和剪力V的關(guān)系為:

        M=EIu″;V=-EIu?

        (22)

        將式(22)和邊界條件式(9)代入式(21)得:

        (23)

        (24)

        由式(11)~(14)可知在裂紋處懸臂梁的位移和它的二階,三階導(dǎo)數(shù)是連續(xù)的,而一階導(dǎo)數(shù)不連續(xù),即:

        由式(24)可得:

        (25)

        將形函數(shù)式(1)和權(quán)函數(shù)代入式(25),由此可得到系統(tǒng)模態(tài)方程:

        (K+Kc-ω2M)U=F

        (26)

        式中:K是懸臂梁系統(tǒng)不含裂紋時(shí)的剛度矩陣,Kc是對(duì)應(yīng)裂紋處的等效剛度矩陣,M是系統(tǒng)的質(zhì)量矩陣,F(xiàn)是外界所施加的載荷。因此,系統(tǒng)節(jié)點(diǎn)的剛度矩陣為:

        〗(27)

        裂紋處的等效剛度矩陣為:

        (28)

        系統(tǒng)節(jié)點(diǎn)的質(zhì)量矩陣為:

        (29)

        系統(tǒng)節(jié)點(diǎn)載荷為:

        (30)

        若為任意數(shù)量的裂紋(28)可寫為:

        (31)

        3 多裂紋診斷方法

        基于無(wú)網(wǎng)格法的裂紋識(shí)別流程如下(見(jiàn)圖3):首先利用曲率模態(tài)法確定裂紋的相對(duì)位置(β1,β2,…βi, βn);然后正反問(wèn)題結(jié)合并利用粒子群優(yōu)化(PSO)算法來(lái)估算該位置的損傷程度(α1,α2,…αi,αn)。該方法結(jié)合了曲率模態(tài)法[8-11]在奇異性分析方面的優(yōu)勢(shì),同時(shí)降低了多裂紋參數(shù)反演尋優(yōu)的維度,簡(jiǎn)化了計(jì)算。

        3.1裂紋損傷位置識(shí)別

        首先用實(shí)驗(yàn)?zāi)B(tài)分析(EMA)測(cè)量獲得裂紋梁的第m階模態(tài)Ym,并計(jì)算得到其模態(tài)曲率Vm。通過(guò)中心差分法近似得到該點(diǎn)的模態(tài)曲率[16]如下:

        (32)

        式中:hc為相鄰兩點(diǎn)的距離,j表示空間不同位置的點(diǎn),Ym,j-1,Ym,j,Ym,j+1,是Ym上三個(gè)連續(xù)的點(diǎn)。裂紋位置可以從模態(tài)曲率的外觀上確定,如果有突變或峰值出現(xiàn),則該位置有損傷存在,從而確定損傷位置β1,β2,…βi,βn。

        3.2裂紋損傷程度求解

        從圖3可知,在正問(wèn)題分析中,利用無(wú)網(wǎng)格局部PETROV-GALERKIN法,求得振動(dòng)特征方程得到了裂紋梁自由振動(dòng)的固有頻率ωm(m=1, 2, 3, …), 即對(duì)于給定的裂紋相對(duì)位置β,求解不同α相關(guān)的各階總體振動(dòng)頻率ωm=ω(α1,α2,…αj,αn)可得到裂紋梁定量識(shí)別數(shù)據(jù)庫(kù)。

        在反問(wèn)題分析中,通過(guò)實(shí)測(cè)的ω*反演優(yōu)化求解出裂紋的相對(duì)深度?,F(xiàn)在構(gòu)造如下目標(biāo)函數(shù)(33)。

        (33)

        式中:ωj為數(shù)值計(jì)算的不同裂紋相對(duì)深度組合的數(shù)據(jù)庫(kù),ωj*表示實(shí)測(cè)頻率。利用粒子群優(yōu)化方法(有關(guān)粒子群優(yōu)化方法詳見(jiàn)文獻(xiàn)[18])從數(shù)據(jù)庫(kù)中尋求最優(yōu)解以獲得真實(shí)裂紋的參數(shù)(α1,α2,…αi,αn)。

        曲率模態(tài)差定位懸臂梁裂紋位置(β1,β2,…βi,βn)

        圖3 裂紋診斷流程圖Fig.3 Crack diagnosis flow chart

        4 數(shù)值仿真

        4.1裂紋梁固有頻率計(jì)算

        為了驗(yàn)證該算法的有效性,以文獻(xiàn)[17]中的雙裂紋懸臂梁進(jìn)行驗(yàn)證,參數(shù)如下:彈性模量E=210 GPa,密度ρ=7 860 kg/m3,懸臂梁有效長(zhǎng)度L=0.5 m,橫截面高度H=0.02 m,寬度B=0.012 m。

        用本文無(wú)網(wǎng)格法計(jì)算雙裂紋梁的前三階固有頻率并與文獻(xiàn)[17]中的頻率對(duì)比,列于表1。其中α表示裂紋相對(duì)深度,β表示裂紋相對(duì)位置,ω表示固有頻率,括號(hào)內(nèi)數(shù)據(jù)是來(lái)自文獻(xiàn)[17]中的頻率??梢钥闯鲈摲椒ㄇ蟮霉逃蓄l率誤差小于1%,滿足精度要求。為裂紋識(shí)別打下了基礎(chǔ)。

        4.2裂紋識(shí)別仿真

        利用曲率模態(tài)法先找到雙裂紋懸臂梁的裂紋位置β1,β2,然后以該裂紋位置實(shí)測(cè)頻率ω*作為輸入,利用PSO[18]算法尋優(yōu)求解出與輸入值相差最小的樣本點(diǎn),從而獲得該裂紋位置相應(yīng)的損傷程度α1,α2。

        為了模擬實(shí)測(cè)頻率中含有噪聲信號(hào)的影響,現(xiàn)加入±2%高斯白噪聲[19](WGN)后的固有頻率為:

        (34)

        由表2可知,在不同的裂紋工況下,定量診斷的兩個(gè)裂紋的相對(duì)深度α誤差在8%以內(nèi)。表明此無(wú)網(wǎng)格法在低水平的噪聲環(huán)境下,可以可靠地進(jìn)行裂紋梁的定量診斷。

        圖4 雙裂紋位置模態(tài)曲率形狀圖Fig.4 Two crack locations detection using the mode shape curvature

        圖5 粒子群優(yōu)化(PSO)的收斂進(jìn)程和對(duì)應(yīng)識(shí)別兩個(gè)裂紋的相對(duì)深度Fig.5 Convergence process of PSO and identify the two corresponding relative crack depths

        5 實(shí)驗(yàn)研究

        圖6 實(shí)驗(yàn)裝置圖Fig.6 Experimental setup

        表1 雙裂紋懸臂梁前三階固有頻率

        表2 四組裂紋工況的識(shí)別結(jié)果

        圖7 二階模態(tài)形狀與曲率Fig.7Thesecondmodeshapeandcurvemodeshape圖8 PSO識(shí)別對(duì)應(yīng)位置裂紋相對(duì)深度Fig.8PSOidentifycorrespondingrelativecrackdepths

        6 結(jié) 論

        (1) 提出一種用于梁結(jié)構(gòu)多裂紋定量識(shí)別的無(wú)網(wǎng)格法。通過(guò)局部加權(quán)殘量法將裂紋以虛擬邊界引入,從而將裂紋融入系統(tǒng)振動(dòng)特征方程中,構(gòu)造懸臂梁裂紋參數(shù)定量識(shí)別正問(wèn)題模型以獲得裂紋診斷數(shù)據(jù)庫(kù);并通過(guò)實(shí)測(cè)的裂紋系統(tǒng)固有頻率,結(jié)合曲率模態(tài)法先定位再利用粒子群優(yōu)化(PSO)算法反演數(shù)據(jù)庫(kù)以尋最優(yōu)解從而定量識(shí)別出懸臂梁系統(tǒng)裂紋的相對(duì)深度。

        (2) 數(shù)值仿真研究表明:基于無(wú)網(wǎng)格法對(duì)梁結(jié)構(gòu)多裂紋定量識(shí)別方法切實(shí)可行,且只需要節(jié)點(diǎn)信息,前處理簡(jiǎn)單,因而非常適合用來(lái)模擬非連續(xù)(裂紋)問(wèn)題,為結(jié)構(gòu)多裂紋識(shí)別提供了新途徑。

        [1] 張敬芬,趙德有. 工程結(jié)構(gòu)裂紋損傷振動(dòng)診斷的發(fā)展現(xiàn)狀和展望[J].振動(dòng)與沖擊, 2002, 21(4): 22-26.

        ZHANG Jin-fen, ZHAO De-you. The current status and future development of damage diagnosis of vibration engineering structural crack [J]. Journal of Vibration and Shock, 2002, 21(4): 22-26.

        [2] 羅毅,甄立敬. 基于小波包與倒頻譜分析的風(fēng)電機(jī)組齒輪箱齒輪裂紋診斷方法[J].振動(dòng)與沖擊,2015,34(3):210-214.

        LUO Yi,ZHEN Li-jing. Diagnosis method of turbine gearbox crack based on wavelet packet and cepstrum analysis[J]. Journal of Vibration and Shock, 2015, 34(3): 210-214.

        [3] 袁勝發(fā) 李秀瓊. 基于網(wǎng)格式支持向量機(jī)算法的轉(zhuǎn)軸裂紋故障診斷[J].振動(dòng)與沖擊, 2009, 28(9): 155-158.

        YUAN Sheng-fa LI Xiu-qiong. Shaft crack fault diagnosis based on grid support vector machines[J]. Journal of Vibration and Shock, 2009, 28(9): 155-158.

        [4] 陳雪峰, 向家偉, 董洪波, 等. 基于區(qū)間 B 樣條小波有限元的轉(zhuǎn)子裂紋定量識(shí)別[J]. 機(jī)械工程學(xué)報(bào),2007,43(3):123-127.CHEN Xue-feng, XIANG Jia-wei, DONG Hong-bo, et al. Quantitative identification of rotor cracks based on finite element of B-SPLINE wavelet on the interval[J]. Chinese Journal of Mechanical Engineering, 2007,43(3):123-127.

        [5] 李兵, 陳雪峰, 何正嘉. 基于小波有限元的懸臂梁裂紋遺傳優(yōu)化辨識(shí)[J]. 振動(dòng)與沖擊, 2009, 28(12): 27-29.

        LI Bin, CHEN Xue-feng, HE Zheng-jia. Genetic optimization identification of Cantilever cracks based on finite element of wavelet [J]. Journal of Vibration and Shock, 2009, 28(12): 27-29.

        [6] Xiang J W, Jiang Z S, Wang Y X,et al. Study on damage detection software of beam-like structures[J]. Structural Engineering and Mechanics,2011,39(1):77-91.

        [7] Atluri S N, Zhu T. A new meshless local petrov-galerkin approach in computational mechanics[J]. Computational Mechanics,1998,2:117-127.

        [8] 彭華,游春華,孟勇. 模態(tài)曲率差法對(duì)梁結(jié)構(gòu)的損傷診斷[J]. 工程力學(xué),2006, 23(7): 49-53.

        PEN Hua, YOU Chun-hua, MENG Yong. Damage diagnosis of beam structures by modal curvature difference method [J]. Engineering Mechanics, 2006, 23(7): 49-53.

        [9] 何欽象, 楊智春, 姜峰, 等. 薄板損傷檢測(cè)的高斯曲率模態(tài)差方法[J]. 振動(dòng)與沖擊, 2010, 29(7): 112-115.

        HE Qin-xiang, YANG Zhi-chun, JIANG Feng, et al. Damage detection for an elastic thin plate based on Gauss curvature modal difference [J]. Journal of Vibration and Shock, 2010, 29(7): 112-115.

        [10] Wang Y F,Liang M, Xiang J W. Damage detection method for wind turbine blades based on dynamics analysis and mode shape difference curvature information[J]. Mechanical Systems and Signal Processing,2014, 48(2): 351-367.

        [11] Ravi Prakash Babu K, Durga Prasad G. Crack detection in beams from the differences in curvature mode shapes[J]. Journal of Structural Engineering, 2012, 39(2): 237-242.

        [12] Rizos P F, Aspragatos N, Dimarogonas A D. Identification of crack location and magnitude in a cantilever beam from the vibration modes[J]. Journal of Sound and Vibration,1990,138(3): 381-388.

        [13] Zheng D Y, Fan S C. Natural frequency changes of a non-uniform beam with multiple cracks via modified Fourier series[J]. Journal of Sound and Vibration,2001,242(5):701-717.

        [14] 何永勇,褚福磊,郭丹,等.基于遺傳算法的旋轉(zhuǎn)機(jī)械轉(zhuǎn)子裂紋識(shí)別的研究[J].機(jī)械工程學(xué)報(bào),2001,37(10):69-74.

        HE Yong-yong, CHU Fu-lei, GUO Dan, et al. Research on rotating machinery rotor cracks identification based on genetic algorithm [J]. Chinese Journal of Mechanical Engineering, 2001,37(10): 69-74.

        [15] Fan S C, Sheng N. Meshless formulation using NURBS basis functions for eigenfrequency changes of beam having multiple open-cracks[J]. Journal of Sound and Vibration, 2004, 269(3):781-793.

        [16] Pandey A K, Biswas M, Samman M M.Damage detection from changes in curvature mode shapes [J]. Journal of Sound and Vibration,1991,145(2): 321-332.

        [17] Lee J. Identification of multiple cracks in a beam using natural frequencies[J]. Journal of Sound and Vibration,2009,320:482-490.

        [18] Birge B. PSOt-a particle swarm optimization toolbox for use with Matlab. Proceedings of the 2003 IEEE Swarm Intelligence Symposium (SIS03).

        [19] Suk M,Gillis D. Effect of mechanical design of the suspension on dynamic loading process [J] Microsystem Technologies,2005,11(8/9/10):846-850.

        Multi-crack detection of beam structures based on meshless method

        LUO Zhi-gang1, JIANG Zhan-si1, WANG Yan-xue1, XIANG Jia-wei2

        (1. College of Mechanical and Electronic Engineering, Guilin University of Electronic Technology, Guilin 541004, China;2. College of Mechanical and Electronic Engineering, Wenzhou University, Wenzhou 325035, China)

        Multi-crack of a cantilever beam structure were detected based on the meshless method. With the local weighted residual method, cracks were taken as virtual boundaries to be integrated into a discrete system’s equations to build a multi-crack beam dynamic model and acquire a data base for identifying quantitatively a positive problem. Combined with the curvature modal analysis method, the crack positions were found. The measured frequencies were taken as inputs, the particle swarm optimization method was used for the optimization of an inverse problem to predict quantitatively the damage level of cracks. Namerical simulations showed that the proposed method is simple and suitable for simulating non-continuous crack problems, and it provides a new way for detecting multi-crack of structures.

        meshless method; mult-crack; identification

        10.13465/j.cnki.jvs.2016.13.033

        國(guó)家自然科學(xué)基金(51565008;51575400);廣西制造系統(tǒng)與先進(jìn)制造技術(shù)重點(diǎn)實(shí)驗(yàn)室開(kāi)放課題(14-045-15-009Z);桂林電子科技大學(xué)研究生創(chuàng)新項(xiàng)目(YJCXS201506)

        2015-05-20修改稿收到日期:2015-07-01

        羅志鋼 男,碩士生,1987年生

        蔣占四 男,博士,副教授,1977年生

        E-mail:jiangzhansi@gmail.com

        TU223;O346.1

        A

        猜你喜歡
        網(wǎng)格法曲率固有頻率
        大曲率沉管安裝關(guān)鍵技術(shù)研究
        一類雙曲平均曲率流的對(duì)稱與整體解
        現(xiàn)場(chǎng)測(cè)定大型水輪發(fā)電機(jī)組軸系的固有頻率
        雷擊條件下接地系統(tǒng)的分布參數(shù)
        半正迷向曲率的四維Shrinking Gradient Ricci Solitons
        角接觸球軸承的優(yōu)化設(shè)計(jì)算法
        基于遺傳算法的機(jī)器人路徑規(guī)劃研究
        基于GIS的植物葉片信息測(cè)量研究
        總溫總壓測(cè)頭模態(tài)振型變化規(guī)律研究
        A novel functional electrical stimulation-control system for restoring motor function of post-stroke hemiplegic patients
        无码人妻丰满熟妇啪啪网站 | 久久99精品久久久久婷婷| 色欲欲www成人网站| 国产呦系列呦交| 老熟妇高潮喷了╳╳╳| 国产成人77亚洲精品www| 久久久久久一级毛片免费无遮挡| 日本肥老熟妇在线观看 | 国产精品99久久久精品免费观看| 中文字幕 在线一区二区| 婷婷久久国产综合精品| 国产精品自在线拍国产手机版| 1717国产精品久久| 国产亚洲精品国产福利在线观看| 国产 无码 日韩| 香蕉蜜桃av一区二区三区| 欧美疯狂性受xxxxx喷水| 欧美操逼视频| 美女高潮流白浆视频在线观看 | 国产国语一级免费黄片| 少妇精品亚洲一区二区成人| 国产精品人妻一区二区三区四| 看全色黄大色大片免费久久| 亚洲电影中文字幕| 在线不卡精品免费视频| 亚洲三级视频一区二区三区 | 日本边添边摸边做边爱| 国产日产精品一区二区三区四区的特点 | 免费人成视频xvideos入口| 久久久久亚洲AV无码专| 精品女人一区二区三区| 久久偷看各类wc女厕嘘嘘偷窃| 婷婷综合缴情亚洲| 69搡老女人老妇女老熟妇| 亚洲国产中文字幕在线视频综合| 又粗又黑又大的吊av| 亚洲中文有码字幕青青| 国产高清一区在线观看| 亚洲一区二区免费在线观看视频 | 一区二区三区四区亚洲免费| 99久久超碰中文字幕伊人|