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

基于固有頻率的呼吸式裂紋梁損傷識別方法

2017-04-07 09:07:46劉文光郭隆清賀紅林
中國機械工程 2017年6期
關鍵詞:裂紋振動深度

劉文光 郭隆清 賀紅林 顏 龍

南昌航空大學航空制造工程學院,南昌,330063

基于固有頻率的呼吸式裂紋梁損傷識別方法

劉文光 郭隆清 賀紅林 顏 龍

南昌航空大學航空制造工程學院,南昌,330063

將呼吸裂紋梁簡化為由扭轉彈簧連接的兩段彈性梁,在假定振動響應隨振幅變化的基礎上推導出呼吸裂紋梁的固有頻率方程;考慮振動過程中呼吸裂紋的開合情況,假定裂紋梁的剛度是振幅的非線性函數,建立了呼吸裂紋梁的多項式剛度模型;結合等高線裂紋識別理論和方法,提出了一種基于固有頻率的呼吸裂紋梁損傷識別方法,算例驗證了方法的可行性與有效性。研究表明,該方法的識別精度取決于實驗固有頻率的精度。

呼吸裂紋;固有頻率;損傷識別;等高線法

0 引言

梁是機械結構或系統中常見的構件之一,引致梁失效的主要形式之一是疲勞裂紋。服役于振動環境之下的梁,諸如直升機的旋翼葉片、飛機的機翼以及發動機的葉片等,均不可避免地會因裂紋的疲勞擴展及其呼吸效應使得梁的振動特性隨振幅變化而改變,在動力學行為上表現為振動非線性。非線性給呼吸式裂紋梁的振動預測帶來非常大的困擾,不少學者在研究裂紋梁動力問題時,通常都假定裂紋在振動過程中始終處于張開狀態[1-2],以避開裂紋呼吸效應給振動分析帶來的麻煩。這種建模的思路與方法,從物理現象來看,顯然與客觀事實不符。

為了推動基于振動理論的結構健康監測技術的發展,很多學者圍繞含裂紋構件的損傷識別振動方法進行了深入研究。CORNWEL等[2]利用板損傷前和損傷后的振型,基于應變能量法對板結構的損傷進行識別和定位,該方法實質就是STUBBS等[3]提出的一維應變能法的延伸。劉文光等[4]基于結構單元損傷前后模態應變能的變化率構建損傷指標,提出了彈性梁損傷識別模態應變能法。YANG等[5]在假定裂紋附近應變能會發生改變的基礎上,研究了張開式裂紋對彈性梁振動特性的影響。SWAMIDAS等[6]結合能量法和斷裂力學建立了張開型裂紋模型,求解了張開裂紋和無損傷梁的固有頻率。吳國榮[7]利用無質量扭轉彈簧等效裂紋截面,建立了含裂紋梁的等效動力學模型,推導出了含張開式裂紋梁的傳遞矩陣和固有頻率方程??紤]呼吸裂紋閉合效應對裂紋梁振動行為的影響,KISA等[8]建立了考慮振動后半周期裂紋閉合產生的接觸剛度模型,但是此模型僅考慮了裂紋全開和全閉狀態,忽略了裂紋開合中的部分開合狀態。于是,CHENG等[9]通過周期函數描述呼吸裂紋梁的剛度模型,建立了含時變剛度的單自由度模型;BOVSUOVSKY等[10]應用雙線性剛度模型描述裂紋的呼吸行為;ABRAHAM等[11]忽略系統阻尼的影響,通過傅里葉多項式級數模擬了呼吸裂紋剛度的變化。

盡管研究者在含裂紋的結構動力學分析方法上做了大量研究,但很少文獻考慮結構剛度因裂紋呼吸行為隨振幅變化而發生變化,更少研究考慮裂紋接觸剛度和阻尼對損傷識別的影響。本文在考慮裂紋梁剛度隨振幅變化的基礎上,研究一種呼吸式裂紋梁的固有頻率計算方法,結合等高線裂紋識別方法與理論,提出一種基于固有頻率的呼吸式裂紋梁損傷識別技術。

1 呼吸式裂紋梁的固有頻率方程

以圖1a所示的含裂紋懸臂梁為對象,假定梁的長度為L,寬度為w,高度為h,裂紋深度為a,裂紋位置與固定端的距離為L0??紤]裂紋引起的柔度,利用扭轉彈簧代表裂紋,含裂紋梁被描述為通過扭轉彈簧連接起來的兩根無損傷梁,其等效模型如圖1b所示。工作時,一旦梁激發自由振動,裂紋會隨梁振幅的改變交替出現張開與閉合的情況,導致梁的剛度隨振幅變化發生相應變化,使梁的動力學行為發生復雜反應。

(a)幾何模型 (b)等效模型圖1 含疲勞裂紋的懸臂梁模型Fig.1 Model of canlitiever beam with fatigue crack

忽略阻尼和裂紋對梁振型的作用,把裂紋位置的局部剛度描述為非線性無質量扭轉彈簧,則其剛度變化取決于裂紋的開合程度(即與振幅有關)。此疲勞裂紋梁的運動方程可以描述為

(1)

式中,E為材料的彈性模量;I為截面關于中性軸的慣性矩;ρ為材料密度;S為截面面積;w1和w2分別為對應裂紋左右邊彈性梁的橫向位移函數;x為縱向坐標;t為時間。

因為裂紋位置的局部剛度取決振幅函數,所以梁的固有頻率及其振型也隨振幅的變化發生改變。梁的振動響應可表示為

(2)

式中,A(t)為梁軸線運動時指定位置的位移;φ1(x,A(t))和φ2(x,A(t))分別為左右邊彈性梁與振幅有關的特征函數。

把振動響應函數(式(2))代入運動方程(式(1)),可得到

(3)

式中,ci(i=1,2,…,8)為與邊界條件有關的待定系數;ω為圓頻率;λ和ω與時間有關。

呼吸裂紋梁兩端的邊界條件為

(4)

梁裂紋x=L0位置處,相應的連續性條件為

(5)

梁裂紋x=L0位置處,還有一個斜率條件為

EIφ1,xx(L0,A(t))=
kA[φ2,x(L0,A(t))-φ1,x(L0,A(t))]

(6)

式中,kA為隨振幅變化的呼吸裂紋梁剛度。

上述4個邊界條件(式(4))、3個連續性條件(式(5))、1個斜率條件(式(6)),構成了行列式特征值問題,為了得到ci的量綱一解,行列式的特征值必須為零,即

|Δ(ω,L0,A(t))|=0

(7)

由此可推出含呼吸式裂紋梁的固有頻率方程:

(8)

其中

F(β)=[cos(βλL)sinh(βλL)-sin(βλL)cosh(βλL)]-
[cos(βλL)sinh(βλL)+sin(βλL)cos(βλL)+
sinh(βλL)cosh(βλL)+sin(βλL)cosh(βλL)]·
sin(λL)sinh(λL)+[cos(βλL)sinh(βλL)-
sin(βλL)cosh(βλL)]cos(λL)cosh(λL)+
[cos(βλL)cosh(βλL)+cosh2(βλL)+
sin(βλL)sinh(βλL)]sin(λL)cosh(λL)+
[sin(βλL)sinh(βλL)-cos(βλL)cosh(βλL)-
cos2(βλL)]cos(λL)sinh(λL)

式中,β為裂紋相對位置,β=L0/L。

由式(8)可知,呼吸式裂紋梁的固有頻率變化取決于剛度kA的變化,而kA的變化規律與梁的振幅衰減程度、裂紋深度有關。因此,在建立呼吸裂紋梁剛度模型的基礎上,通過固有頻率方程可以分析裂紋處于不同位置、不同深度以及不同振幅條件下的懸臂梁固有頻率。

2 隨振幅變化的含裂紋梁剛度模型

為便于建立含裂紋梁的剛度模型,不考慮系統阻尼的影響。如圖1a所示,假定梁的軸線過平衡位置O向上運動,轉角θ+增大,在此半個運動周期內裂紋將完全閉合,此時,梁的剛度與無裂紋梁的剛度完全一樣;反之,如果梁的軸線過平衡位置O向下運動,則轉角θ-增大,隨著θ-值的增大,裂紋從完全閉合狀態逐漸過渡到完全張開狀態。實際上,受系統阻尼的影響,含裂紋梁的自由振動會隨時間逐漸衰減。換言之,梁自由振動過程中,裂紋的運動會受阻尼的抑制,即梁的剛度變化范圍逐漸縮小。

假設采用二次多項式近似描述裂紋梁剛度與振幅之間的關系,其表達式為

kA=a0+a1A+a2A2

(9)

式中,A為軸線自由端的振幅;a0、a1、a2為待定系數,取決于設定的邊界條件。

假設裂紋完全張開時,裂紋梁的剛度等于開裂紋梁的剛度,裂紋閉合時,裂紋梁的剛度等于無損傷梁的剛度。所以,裂紋梁的剛度與振幅之間應滿足以下關系:

式中,AO為裂紋全開時軸線自由端的振幅;AC為裂紋全閉時軸線自由端的振幅;kC為裂紋全閉時梁的剛度;kO為裂紋全開時梁的剛度。

利用插值原理和以上條件,可得到含裂紋梁的剛度模型:

(10)

裂紋全開狀態下,含裂紋梁的局部剛度計算表達式為[12]

(11)

J(α)=0.6272α2-0.4533α3+4.5948α4-9.973α5+20.2984α6-33.0351α7+47.1063α8-40.7556α9+19.6α10

式中,μ為材料泊松比;J(α)為量綱一局部柔度;α為相對裂紋深度,α=a/h。

量綱一柔度計算公式對于α∈[0,0.5]情況比較符合。

3 裂紋識別方法與理論

如果定義由裂紋位置和深度確定固有頻率稱為正問題,那么由固有頻率識別裂紋位置和深度則稱為逆問題?;诠逃蓄l率方程,對應于不同的裂紋位置和裂紋深度,可以得到各階次固有頻率,同樣,如果已知梁的固有頻率,也可以由已知的固有頻率與裂紋位置和深度的關系式,求解出裂紋的位置和深度。含呼吸式裂紋梁的固有頻率方程是考慮了裂紋位置、裂紋深度和振動幅值的三元函數,為方便損傷識別,若取呼吸裂紋開合過程中梁固有頻率的平均值作為對應裂紋位置、裂紋深度的固有頻率值,固有頻率就簡化為裂紋位置、深度的二元函數。從求解線性方程組的角度來看,通過固有頻率識別裂紋位置、裂紋深度,只需知道兩階固有頻率就可以求得唯一解。為了避免振動測量中的節點效應并考慮工程測試中高階模態頻率誤差相對較大的實情,實測頻率一般僅取前三階固有頻率值進行識別。

結合裂紋梁的固有頻率方程和剛度模型,基于等高線法可識別出裂紋位置和深度。首先,擬合出以裂紋位置和深度為自變量的固有頻率響應面;然后,將實測固有頻率代入到頻率響應面上,得到若干條頻率等高線并投影到同一平面,等高線理論上必然存在一個交點,這個交點橫縱坐標即對應裂紋位置和裂紋深度。圖2為對應某一裂紋參數的固有頻率等高線圖,三階固有頻率等高線由三段曲線組成,二階模態頻率等高線由兩段曲線組成,一階模態頻率等高線只有一段曲線[13]。

1.第一階固有頻率等高線 2.第二階固有頻率等高線 3.第三階固有頻率等高線圖2 等高線圖Fig.2 Diagram of contour line

工程實際中,通過實測頻率所得的三條等高線可能無法準確地交于一點,此時可以選取三條曲線互相交點所形成的小三角形的形心來確定裂紋位置和深度。

4 算例分析與驗證

采用兩個算例驗證本文所提方法。算例1為利用裂紋梁解析固有頻率進行驗證,算例2為利用文獻測得的實驗固有頻率進行驗證。每個算例均采用6種不同裂紋工況呼吸裂紋梁的前三階固有頻率進行分析。

4.1 算例1

假設含呼吸裂紋懸臂梁的幾何尺寸:梁長度L=300 mm,橫截面高度h=8 mm,寬度w=20 mm。材料的彈性模量E=68.9 GPa,密度ρ=2750 kg/m3。對于已知裂紋相對位置和相對深度,基于固有頻率方程可以構造考慮裂紋位置、深度的固有頻率響應曲面。為了驗證該方法的合理性,設計出6種不同的裂紋工況,并利用固有頻率方程解出梁的前三階固有頻率,見表1。

表1 不同裂紋條件下的前三階固有頻率

把表1中各種工況條件下的前三階固有頻率代入裂紋梁固有頻率響應面,可得到對應的固有頻率等高線,將等高線投影到水平面,得到對應工況裂紋相對位置和深度的關系曲線,如圖3所示。

(a)工況a

(b)工況b

(c)工況c

(d)工況d

(e)工況e

(f)工況f1.第一階固有頻率等高線投影 2.第二階固有頻率等高線投影 3.第三階固有頻率等高線投影圖3 裂紋參數識別1Fig.3 Case 1 of crack parameters identification

從圖3中可以看出,三條等高線非常精準地相交于一點,交點的橫縱坐標即為對應裂紋的相對位置和相對深度,因而可以非常方便地得到裂紋參數。比較識別出的裂紋參數與試件的實際解析解裂紋參數,計算兩者之間的誤差,結果如表2所示。

表2 等高線法懸臂梁裂紋參數識別結果1

從表2中可以看出,裂紋相對深度的識別效果的相對誤差不會超過0.32%,裂紋相對位置的識別效果的相對誤差不會超過0.93%。對于有效長度為300 mm、截面高度為20 mm的懸臂裂紋梁,裂紋深度識別誤差不超過0.064 mm,裂紋位置的識別誤差不超過2.79 mm,精度相當高,這證明了用等高線裂紋識別方法對裂紋參數進行檢測的可行性。

4.2 算例2

實測固有頻率通過振動測試技術可以得到,它不依賴于所選擇的結構動力學模型,可以與動力學分析同步進行。下面的裂紋識別將以實測固有頻率作為輸入進行裂紋參數的識別。

取文獻[14]中含裂紋梁的材料屬性和幾何尺寸進行損傷識別驗證。材料彈性模量E=206 GPa,密度ρ=7750 kg/m3,梁的有效長度L=300 mm,梁的橫截面高度h=20 mm,寬度w=20 mm?;谝阎膸缀纬叽纭⒉牧蠀岛凸逃蓄l率方程可構造梁的固有頻率響應面。采用文獻[14]給出的實驗測得的各種裂紋工況的前三階固有頻率值(表3)可進一步驗證本文方法的合理性。

將表3中各裂紋工況的前三階實測固有頻率代入固有頻率響應面,得到相應的頻率等高線,并將等高線投影到水平面,得到裂紋梁的相對裂紋深度和裂紋位置的關系曲線,如圖4所示。

從圖4中可以看出,對于實測固有頻率,三條等高線并不交于一點,故無法得到三條等高線的交點,只能確定兩兩之間的交點。此時,可以選取三條曲線交點所形成的小三角形的形心來確定裂紋相對位置和相對深度。比較識別出的裂紋參數與試件實際參數,計算兩者之間的誤差,表4給出了裂紋參數識別結果和誤差。

表3 實測6種工況條件下前三階固有頻率

由表4誤差分析可知,實測固有頻率下,除去工況⑥裂紋深度誤為14.0%,其余均在7.0%以內;裂紋相對位置的識別誤差不超過4.4%。對于橫截面高度為20 mm、懸臂梁有效長度為300 mm的裂紋梁來說,裂紋深度識別誤差在1.4 mm以內,裂紋位置的識別誤差在13.2 mm以內。總的來說,裂紋參數識別精度較高,具有較好的工程應用前景。

表4 等高線法懸臂梁裂紋參數識別結果2

4.3 誤差分析與討論

影響裂紋梁損傷識別精度的因素有很多,但是主要是正問題的求解精度和實測固有頻率的測量精度。原因在于:

(1)如果通過呼吸裂紋梁固有頻率方程構造的固有頻率響應曲面精度很差,也就是裂紋模型的建立不可靠,則裂紋識別精度就難以保證。

(2)實測固有頻率的精度影響因素分為內因和外因兩種。內因是裂紋處于某階次或某幾階次模態振型的節點處或在節點附近時,裂紋變化對梁的固有頻率影響很小,幾乎可以忽略,從而導致測量結果無法用于裂紋識別。外因主要是傳感器引起的附加重量的影響、懸臂梁的懸臂邊界難以模擬、模態測量儀器的測試精度等原因,容易造成固有頻率測量誤差。實測固有頻率誤差直接影響投影所得的固有頻率等高線的準確性,進而影響若干等高線的交點位置,即影響了裂紋相對位置和相對深度等參數的識別效果。

(a)工況①

(b)工況②

(c)工況③

(d)工況④

(e)工況⑤

(f)工況⑥1.第一階固有頻率等高線投影 2.第二階固有頻率等高線投影 3.第三階固有頻率等高線投影圖4 裂紋參數識別2Fig.4 Case 2 of crack parameters identification

5 結論

(1)把呼吸式裂紋梁簡化為由扭轉彈簧連接的兩段彈性梁,在假定振動響應隨振幅變化的基礎上推導出含呼吸裂紋梁的固有頻率方程。

(2)考慮呼吸裂紋的開合情況,假定裂紋梁的剛度是振幅的非線性函數,建立了呼吸裂紋梁的多項式剛度模型。

(3)結合等高線裂紋識別理論和方法,提出了一種基于固有頻率的呼吸裂紋梁損傷識別新方法。通過算例驗證了方法的可行性與有效性,方法的識別精度取決于實驗固有頻率的精度。

[1] CADDEMI S, CALIO I. Exact Closed-form Solution for the Vibration Modes of the Euler-bernoulli Beam with Multiple Open Cracks[J]. Journal of Sound and Vibration, 2009, 327(1): 473- 489.

[2] CORNWELL P, DOEBLING S W, FARRAR C R. Application of the Strain Energy Damage Detection Method to Plate-like Structures[J]. Journal of Sound and Vibration, 1999, 224(2):359-374.

[3] STUBBS N, KIM J T, FARRAR C R. Field Verification of a Non-destruction Damage Localization and Sensitivity Estimator Algorithm[C]// Proceeding of the 13th International Modal Analysis Conference. Nashville, Tennessee, 1995:210-218.

[4] 劉文光,李俊,嚴鋮,等.彈性梁損傷識別模態應變能法研究[J].中國機械工程,2014,25(12):1651-1654. LIU Weguang, LI Jun, YAN Chen, et al. Modal Strain Energy Method for Damage Detection of an Elastic Beam[J]. China Mechanical Engineering, 2014, 25(12):1651-1654.

[5] YANG X F, SWAMIDAS A S J, SESHADRY R. Crack Identification in Vibration Beams Using the Energy Method[J]. Journal of Sound and Vibration, 2001, 224(2):339-357.

[6] SWAMIDAS A S J, YANG X F, SESHADRY R. Identification of Cracking in Beam Structures Using Timoshenko and Euler Formulation[J].Journal of Engineering Mechanics, 2004, 130 (11) : 1297-1308.

[7] 吳國榮.含裂紋梁自由振動分析[J].船舶力學,2007,11(5):798-803. WU Guorong. Free Vibration Analysis of Crackedbeam[J].MechanicsofShip, 2007, 11(5):798-803.

[8]KSIAM,BRANDONJ.TheEffectsofClosureofCracksontheDynamicsofaCrackedCantileverbeam[J].JournalofSoundandVibration, 2000, 238(1):1-18.

[9]CHENGSM,WUXJ,WALLACEW,etal.VibrationalResponseofaBeamwithaBreathingCrack[J].JournalofSoundandVibration, 1999, 225(1):201-208.

[10]BOVSUNOVSKYAP,SURACEC.ConsiderationsRegardingSuper-harmonicVibrationofaCrackedBeamandtheVibrationintheDampingCausedbythePresenceoftheCrack[J].JournalofSoundandVibration, 2005, 288:865-886.

[11]ABRAHAMONL,BRANDONJA.ModelingoftheOpeningandClosureofaCrack[J].JournalofVibrationandAcoustics, 1995,117:370-377.

[12]ZHANGW,TESTARB.ClosureEffectonFatigueCrackDetection[J].EngineeringMechanics, 1999, 125(2):1125-1132.

[13]GMOAS.CrackDetectioninBeamsUsingChangesinFrequenciesandAmplitudesofFrequencyResponseFunctions[J].JournalofSoundandVibration, 2003, 265(1): 1-22.

[14]LEEYS,CHUNGMJ.AStudyonCrackDetectionUsingEigen-frequencyTestData[J].ComputersandStructures, 2000, 77(3): 327-342.

(編輯 蘇衛國)

A Damage Identification Method of a Breathing Cracked Beam by Natural Frequency

LIU Wenguang GUO Longqing HE Honglin YAN Long

School of Aeronautic Manufacturing Engineering,Nanchang Hangkong University,Nanchang,330063

Firstly, a cracked cantilever beam was simplified into two elastic beams connected by a torsion spring. It was supposed that the vibration response varied as the vibration amplitude of the beam. And the natural frequency equation of a breathing cracked beam was derived. Then, the open and close conditions of the breathing crack during vibrations were considered. A polynomial stiffness model of the breathing cracked beam was built up by the postulation that the stiffness of the cracked beam was one of nonlinear functions of the vibration amplitudes. Thereafter, a damage identification method of a breathing cracked beam by natural frequency was proposed combined with the natural frequency equation and the stiffness model. The feasibility and effectiveness of this method was verified by case studies in the end. Results indicate that the identification precision of this method depends on the testing precision of the experimental natural frequency.

breathing crack; natural frequency; damage identification; contour method

2015-12-24

國家自然科學基金資助項目(51565039)

V215.4DOI:10.3969/j.issn.1004-132X.2017.06.011

劉文光,男,1978年生。南昌航空大學航空制造工程學院副教授、博士。主要研究方向為飛行器結構動力學及疲勞壽命預測。E-mail:liuwg14@163.com。郭隆清,男,1991年生。南昌航空大學航空制造工程學院碩士研究生。賀紅林,男,1967年生。南昌航空大學航空制造工程學院教授、博士。顏 龍,男,1992年生。南昌航空大學航空制造工程學院碩士研究生。

猜你喜歡
裂紋振動深度
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
裂紋長度對焊接接頭裂紋擴展驅動力的影響
深度理解一元一次方程
振動與頻率
天天愛科學(2020年6期)2020-09-10 07:22:44
Epidermal growth factor receptor rs17337023 polymorphism in hypertensive gestational diabetic women: A pilot study
深度觀察
深度觀察
深度觀察
中立型Emden-Fowler微分方程的振動性
預裂紋混凝土拉壓疲勞荷載下裂紋擴展速率
主站蜘蛛池模板: 亚洲人成网站色7799在线播放| 亚洲男人天堂2018| 国产高清在线丝袜精品一区| 一本大道视频精品人妻 | 99精品一区二区免费视频| 18禁影院亚洲专区| 亚洲中文字幕无码mv| 亚洲欧美国产高清va在线播放| 亚洲天堂网2014| 日本手机在线视频| 四虎成人免费毛片| 国产成熟女人性满足视频| 亚洲欧美日韩视频一区| 中文字幕有乳无码| 国产亚洲成AⅤ人片在线观看| 中文字幕2区| 国产女主播一区| 欧美全免费aaaaaa特黄在线| 亚洲国产精品一区二区高清无码久久| 91久久偷偷做嫩草影院免费看| 狠狠操夜夜爽| 偷拍久久网| 青青热久免费精品视频6| 特级毛片8级毛片免费观看| 国产精品大白天新婚身材| 无码日韩精品91超碰| 97se亚洲| 蜜桃视频一区二区三区| 99热这里只有精品免费国产| 男女男精品视频| 国产美女一级毛片| 欧美人与牲动交a欧美精品| 亚洲国产成人久久精品软件| 国产制服丝袜无码视频| 理论片一区| 麻豆国产精品| 亚洲一区二区三区国产精品 | 国产精品综合色区在线观看| 国产欧美日韩在线在线不卡视频| 久久综合九色综合97婷婷| 欧美一级色视频| 97久久精品人人做人人爽| 国产www网站| 99精品视频在线观看免费播放| 在线观看国产小视频| 无码不卡的中文字幕视频| 精品久久国产综合精麻豆| 国产成年女人特黄特色毛片免 | 久久毛片基地| 夜夜操天天摸| 欧美翘臀一区二区三区| AV无码一区二区三区四区| 一本色道久久88亚洲综合| 亚洲天堂自拍| 婷婷亚洲最大| 日韩精品资源| 久久福利片| 中文字幕人妻无码系列第三区| 亚欧成人无码AV在线播放| 91蝌蚪视频在线观看| 噜噜噜综合亚洲| 国产亚洲日韩av在线| 四虎成人在线视频| 精品久久久久久成人AV| 亚洲综合二区| 免费a在线观看播放| 欧美成人手机在线视频| 无码粉嫩虎白一线天在线观看| 好吊妞欧美视频免费| 久久婷婷六月| 91精品啪在线观看国产| a级毛片毛片免费观看久潮| 一本大道香蕉久中文在线播放| 九九精品在线观看| 午夜国产精品视频| 欧美日韩91| 欧美一区中文字幕| 亚洲女同欧美在线| 日韩免费毛片视频| 中文字幕久久亚洲一区| 高h视频在线| 在线免费无码视频|