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

聲子晶體串聯(lián)組合桿振動帶隙研究

2014-09-08 03:09:26舒海生高恩武劉少剛王威遠
振動與沖擊 2014年16期
關鍵詞:有限元振動

舒海生,高恩武,劉少剛,趙 丹,王威遠

(哈爾濱工程大學 機電工程學院702所,哈爾濱 150001)

受天然晶體電子能帶理論啟發(fā),已有對周期復合材料或結構中經(jīng)典波傳播的研究[1]。彈性波在彈性常數(shù)、材料密度、結構尺寸等參數(shù)周期性調制下傳播時會產(chǎn)生聲子帶隙,即一定頻帶范圍內(nèi)彈性波傳播受到抑制或禁止。因此聲子晶體或理論上或減振降噪,對精密機械減振平臺、新型濾波器等工程應用領域均具有重要價值。聲子晶體研究主要集中在帶隙形成機理分析及帶隙計算方法探索等[2-7],實際應用研究[8]較少。研究對象多為單一構型聲子晶體,如一維聲子晶體中梁、桿、軸[2,9]及二維聲子晶體[3,10]中格柵、板等,而對實際應用中組合形式結構件研究較少。對聲聲子晶體平面組合桿情況,如T型桿[10]、角梁等在結合處能實現(xiàn)波形轉換,拓寬帶隙頻率范圍,且已取得二維寬頻較好減振效果。而工程實際中激擾多為三維,需研究三維寬頻減振元件。為此,本文將聲子晶體理念引入串聯(lián)型空間組合桿結構構造聲子晶體串聯(lián)組合桿,綜合利用彎曲、縱向、扭轉帶隙,達到多維寬頻減振效果。具體步驟為:①對聲子晶體串聯(lián)組合桿振動特性進行理論分析,推導總傳遞矩陣;②據(jù)總傳遞矩陣編寫MATLAB程序進行數(shù)值計算,重點研究加載方向及桿件間夾角對聲子晶體串聯(lián)組合桿振動特性影響;③用有限元方法分析振動傳遞特性,證明理論分析的正確性;④用遺傳算法對影響減振性能因素的桿件間夾角進行最優(yōu)化設計;⑤給出結論。

1 聲子晶體串聯(lián)組合桿理論分析

圖1為聲子晶體串聯(lián)組合桿一般結構(以3桿為例,對K(K≥3)桿情況分析過程相同),上(桿Ⅲ)、中(桿Ⅱ)、下(桿Ⅰ)材料分布及周期結構尺寸相同,桿Ⅰ、Ⅱ組成的面稱為面Ⅰ,桿Ⅱ、Ⅲ組成的面為面Ⅱ。為研究方便,規(guī)定桿Ⅰ與Ⅱ、Ⅲ分別位于面Ⅰ、面Ⅱ內(nèi)。桿Ⅰ與Ⅱ,桿Ⅱ與Ⅲ,面Ⅰ與Ⅱ間夾角分別為θ1,θ2,θ3,桿與桿間固結。聲子晶體串聯(lián)組合桿振動可解耦為組合桿面內(nèi)振動及垂直于組合桿平面的面外振動。面內(nèi)振動分面內(nèi)彎曲振動、縱向振動,面外振動分面外彎曲振動、扭轉振動。由于聲子晶體串聯(lián)組合桿為空間立體結構,在桿Ⅰ與Ⅱ、桿Ⅱ與Ⅲ結合處組合桿存在面內(nèi)振動與面外振動相互轉化情況,總體上表現(xiàn)出彎振、縱振、扭振的全耦合特性。

圖1 聲子晶體串聯(lián)組合桿模型

當面內(nèi)振動的彎曲波或縱波在單獨桿Ⅰ、Ⅱ或Ⅲ中傳播時,不同周期間傳遞關系[10]可寫為

(1)

(2)

當面外振動彎曲波或扭轉波在單獨桿Ⅰ、桿Ⅱ或桿Ⅲ中傳播時,不同周期間傳遞關系[10]可寫為

(3)

(4)

為便于敘述,將同周期同材料面內(nèi)振動解系數(shù)列陣與面外振動解系數(shù)列陣組合為

(5)

(6)

據(jù)式(1)~式(6)得:

(7)

(8)

圖1表明,桿Ⅰ在桿Ⅰ與桿Ⅱ結合處為第n周期,而桿Ⅱ在結合處為第1周期。據(jù)桿Ⅰ、桿Ⅱ結合處位移、彎曲角、扭轉角、彎矩、扭矩、剪力、正應力等協(xié)調條件,有

(9)

(10)

(11)

(12)

(13)

(14)

(15)

(16)

(17)

(18)

(19)

(20)

將式(9)~(20)寫成矩陣形式為

(21)

據(jù)桿Ⅱ、桿Ⅲ結合處位移、彎曲角、扭轉角、彎矩、扭矩、剪力、正應力等協(xié)調條件,有:

(22)

(23)

(24)

(25)

(26)

(27)

(28)

(29)

(30)

(31)

(32)

(33)

式中:字符號含義同上。

將式(22)~式(33)寫成矩陣形式為

(34)

聲子晶體串聯(lián)組合桿(3桿)總傳遞矩陣寫為

一般的聲子晶體串聯(lián)組合桿(K根桿)總體傳遞矩陣可寫為

2 聲子晶體串聯(lián)組合桿數(shù)值分析

2.1 聲子晶體單桿數(shù)值分析

便于對比,利用傳遞矩陣法對3周期聲子晶體單桿進行數(shù)值分析。基本參數(shù)設置為:各子段直徑D=5 mm,長a1=50 mm;鋁的彈性模量E1=7.21×1010Pa,質量密度ρ1=2 799 kg/m3,泊松比ν1=0.345 1;有機玻璃彈性模量E2=0.32×1010Pa,質量密度ρ2=1 062 kg/m3,泊松比ν2=0.333 3。單桿縱向、彎曲及扭轉振動傳遞率計算結果見圖2。由圖2看出,在0~8 kHz范圍內(nèi)存在明顯縱向振動帶隙5 900~8 000 Hz;兩明顯彎曲振動帶隙為1 420~2 680 Hz,3 920~5 720 Hz;扭轉振動帶隙為3 600~8 000 Hz。單桿在三種不同形式下振動帶隙重疊部分較小,只能獲得二維減振效果,故需研究能綜合運用三種振動帶隙的結構,實現(xiàn)三維減振目的。

圖2 聲子晶體單桿彎曲、縱向及扭轉振動傳遞率曲線

2.2 聲子晶體串聯(lián)組合桿振動特性研究

據(jù)傳遞矩陣法推導結果,利用MATLAB數(shù)值計算,由加載方向及桿件間夾角兩方面研究3周期單桿組成的串聯(lián)組合桿(三桿)結構振動特性。基本參數(shù)設置同聲子晶體單桿。

2.2.1 加載角對振動特性影響

工程中激勵負載大多為多維,為更好模擬外部各種激擾對桿件同時作用情況,引入加載向量r,方向余弦為cosα1,cosα2,cosα3。激擾在3坐標平面加載見圖3。

圖3 加載方向

針對加載方向余弦(cosα1,0,cos(90°-α1)),(α1從0°逐漸增到90°,步長30°)情況,數(shù)值計算結果見圖4(a)。由圖4(a)看出,α1逐漸增大時在第一彎曲帶隙范圍(1 420~2 680 Hz)內(nèi)衰減量逐漸增大,縱向帶隙范圍(7 100~8 000 Hz)內(nèi)衰減量逐漸減小,第二彎曲帶隙變化規(guī)律不明顯,但總體能反映組合桿由縱向振動占主導地位逐漸轉變?yōu)閺澢駝诱贾鲗У匚唬粡澢鷰蹲饔迷矫黠@,位于縱向帶隙中彎曲波傳輸共振峰影響(6 600,7 640 Hz附近,圖3)亦越顯著,導致α1趨于90°時峰附近衰減效果明顯削弱; 而α1較小時情況

相反,縱向振動逐漸占據(jù)主導地位,縱向帶隙衰減作用得以更好體現(xiàn),衰減顯著,彎曲帶隙內(nèi)衰減有所減弱;位于彎曲帶隙中縱波傳輸共振峰(5 660 Hz附近,圖3)影響加強,導致此峰值附近減振效果明顯減弱。因此,該結構能獲得較好的二維寬頻減振性能。考慮加載方向余弦(cos(90°-α2), cosα2,0)(α2從0°逐漸增到90°,步長30°)結果見圖4(b)。由圖4(b)看出,加載角逐漸增大時第一彎曲帶隙范圍內(nèi)衰減量逐漸減小,第二彎曲帶隙范圍(1 400~2 400 Hz)內(nèi)情況較復雜,規(guī)律不明顯,縱向帶隙范圍內(nèi)衰減量逐漸增大,總體有組合桿由彎曲振動占主導地位向縱向振動占主導地位轉變趨勢。

考慮加載方向余弦(0,cos(90°-α3),cosα3)(α3從0°逐漸增到90°,步長30°)時仿真結果見圖4(c)。由圖4(c)看出,面內(nèi)、外加載對整個聲子晶體組合桿減振頻帶范圍影響不同:第一彎曲帶隙1 400~2 120 Hz、第二彎曲帶隙4 360~5 000 Hz范圍內(nèi),隨加載角逐漸增大減振幅值增大,而第一彎曲帶2 120~2 600 Hz、第二彎曲帶隙3 800~4 360 Hz范圍內(nèi),隨加載角逐漸增大減振效果減弱。故在1 400~2 120 Hz、4 360~5 000 Hz頻帶范圍內(nèi)面內(nèi)彎曲振動占主導地位,在2 120~2 600 Hz、3 800~4 360 Hz范圍內(nèi)面外彎曲振動占主導地位。低頻段0~1 000 Hz內(nèi)聲子晶體串聯(lián)組合桿表現(xiàn)出一定減振性能,由于該頻帶不在聲子晶體桿或梁的帶隙內(nèi),因而該減振能力已不能用帶隙解釋。事實上此頻段內(nèi)無論彎曲波或縱波波長均遠大于桿Ⅰ特征尺寸,因此該彎曲波、縱波在桿中的傳播類似于桿或梁的整體運動,即純剛體運動,此時桿Ⅱ、Ⅲ可視為懸臂梁,與桿Ⅰ形成傳統(tǒng)隔振系統(tǒng)中“懸臂梁-振子”體系,產(chǎn)生相似的隔振效果。由于“懸臂梁-振子”隔振體系固有的不足,在該頻段內(nèi)存在較多共振峰。在桿件結合處存在波形轉換,其它頻帶范圍內(nèi)也出現(xiàn)共振峰,可通過增加阻尼等方法進行削減。

(a) 加載方向余弦(cosα1,0 ,cos(90°-α1))聲子晶體串聯(lián)組合桿傳遞率曲線(b) 加載方向余弦(cos(90°-α2),cosα2,0)聲子晶體串聯(lián)組合桿傳遞率曲線(c) 加載方向余弦(0,cos(90°-α3), cosα3)聲子晶體串聯(lián)組合桿傳遞率曲線

總之,無論采用何種角度加載,除局部共振峰不利影響外,整個面內(nèi)、外振動帶隙范圍內(nèi)均能產(chǎn)生較強減振效果,此表明該聲子晶體串聯(lián)組合桿可實現(xiàn)空間內(nèi)三維寬頻減振目的。由于“懸臂梁-振子”效應,低頻范圍內(nèi)隨加載角的變化也有一定減振效果,如加載方向余弦(0,0,1)時,在920~1 200 Hz范圍內(nèi)能形成達-28 dB的減振衰減區(qū)域。

比較圖4發(fā)現(xiàn),加載角從0°~90°變化時,對應的傳遞率曲線族基本介于0°與90°兩條傳遞率曲線之間,因此可將該兩條傳遞率曲線作為包絡線進行重點研究;加載角為0°或90°時衰減量最顯著,故應用中選加載角時應據(jù)振源特點盡量用0°或90°方向角。

2.2.2 桿件夾角對振動傳遞特性影響

為進一步了解聲子晶體串聯(lián)組合桿減振性能,對聲子晶體桿間不同夾角β1,β2,β3(β1為桿Ⅰ與桿Ⅱ間夾角,β2為桿Ⅱ與桿Ⅲ間夾角,β3為桿Ⅲ、桿Ⅰ、桿Ⅱ組成的面夾角)進行分析。加載方向余弦為(0,0,1)時β1-90°-90°(β1由30°~90°,步長30°)振動傳遞率曲線對比見圖5(a),加載方向余弦為(0,0,1)時90°-β2-90°振動傳遞率曲線對比見圖5(b),加載角為(0,0,1)時90°-90°-β3振動傳遞率曲線對比見圖5(c)。通過比較不難發(fā)現(xiàn),β1對高頻范圍(7 200~8 000 Hz)有較大影響(圖5(a)),隨β1增大衰減量逐漸增大;β2對中低頻范圍(1 420~2 680 Hz)作用明顯,(圖5(b)),β2為90°時減振效果明顯優(yōu)于其它兩種;β3主要影響中高頻范圍(3 940~5 200 Hz),隨β3增大衰減量逐漸減小。對其它加載方向余弦如(1,0,0),(0,1,0)結論類似。由于夾角對減振性能有顯著影響,實際應用中組合桿夾角選擇尤其重要,以組合桿結構為例,隔振頻率在較高的7 200~ 8 000 Hz范圍內(nèi)時,β1應選90°夾角;隔振頻率為中高頻率3 940~5 200 Hz范圍內(nèi)時,β3應設計的小點;隔振頻率在中低頻率1 420~2 680 Hz范圍內(nèi)時,β2則應選90°夾角。總之,聲子晶體組合桿的減振性能為β1,β2,β3綜合作用結果,應針對不同減振需求,對桿件夾角進行優(yōu)化設計。

(a) 桿件夾角為β1-90°-90°聲子晶體串聯(lián)組合桿傳遞率曲線(b) 桿件夾角為90°-β2-90°聲子晶體串聯(lián)組合桿傳遞率曲線(c) 桿件夾角為90°-90°-β3聲子晶體串聯(lián)組合桿傳遞率曲線

3 聲子晶體串聯(lián)組合桿有限元分析

為驗證理論分析的正確性,進行聲子晶體串聯(lián)組合桿有限元仿真。重點研究桿件間夾角為90°情況,基本設置為:① 材料參數(shù):連接塊密度4 000 kg/m-3,彈性模量2E11Pa ,泊松比0.3,其余材料參數(shù)同數(shù)值計算。② 網(wǎng)格劃分:用自由網(wǎng)格劃分,單元類型SOLID186,桿節(jié)點數(shù)17 010,單元數(shù)3 056。③加載:簡諧位移激勵幅值0.001 m。④ 分析設置:用諧響應分析及用Full算法。

加載角方向余弦為(1,0,0)、(0,1,0)、(0,0,1)的有限元仿真及理論分析對比結果見圖6,其中實線為理論分析結果,虛線為有限元仿真結果。由圖6看出,三種加載的理論分析與有限元分析傳遞率曲線走勢基本相同,即:加載方向為(1,0,0)時,在1 400~2 600 Hz、4 000~5 000 Hz及6 100~8 000 Hz范圍內(nèi)理論分析與有限元結果均出現(xiàn)較大振動衰減量(圖6(a));加載方向為(0,1,0)時,在1 400~2 800 Hz及 4 000~5 600 Hz范圍內(nèi)理論分析與有限元仿真?zhèn)鬟f率曲線較吻合(圖6(b));加載方向為(0,0,1)時,在1 420~ 及3 941~5 201 Hz范圍內(nèi)理論分析與有限元分析傳遞率曲線吻合較好(圖6(c))。雖理論分析結果基本符合有限元結果,但仍存在一定誤差,其主要原因為計算模型與仿真模型有誤差,有限元仿真模型桿件間通過連接塊連接,而計算模型則視桿件連接處為質點。聲子晶體組合桿仿真結果已有力驗證理論分析的正確性,且減振性能良好。

(a) 加載方向余弦(1,0,0)聲子晶體串聯(lián)組合桿傳遞率曲線(b) 加載方向余弦(1,0,0)聲子晶體串聯(lián)組合桿傳遞率曲線(c) 加載方向余弦(0,0,1)聲子晶體串聯(lián)組合桿傳遞率曲線

4 桿件夾角最優(yōu)化

由以上分析知,聲子晶體組合桿振動特性受多因素(如β1,β2,β3)綜合影響。聲子晶體研究大多針對正問題,即已知聲子晶體結構、材料求帶隙特性,而工程中更需據(jù)所需振動特性如帶隙位置設計合理的聲子晶體結構,即聲子晶體反問題。在理論分析基礎上將遺傳算法引入聲子晶體組合桿結構設計,對桿件夾角參數(shù)進行優(yōu)化,篩選符合要求的最優(yōu)結構。設所需減振范圍1 500~2 500 Hz,設計變量為β1,β2,β3,建立優(yōu)化模型。目標函數(shù)為

式中:DB(β1,β2,β3)為組合桿減振幅值。由于遺傳算法要求目標函數(shù)必須為正值,故將每組減振幅值增加100 dB作為補償,取1 500~2 500 Hz范圍內(nèi)均值作為目標函數(shù)。

約束取π/6≤β1,β2,β3≤π/2。遺傳算法具體參數(shù)設置為每一代個體數(shù)NIND =40,最大遺傳代數(shù)MAXGEN=200,代溝GGAP= 0.9。采用二進制編碼,編碼位數(shù)20。運算采用隨機遍歷抽樣算子,交叉運算采用兩點交叉算子,變異運算采用離散變異算子。通過遺傳算法求解獲得目標函數(shù)收斂曲線見圖7。由圖7看出,經(jīng)200代優(yōu)化獲得結果為β1=30°,β2=89.63°,β3=83.48°。在中低頻范圍1 500~2 500 Hz內(nèi),夾角β2對振動特性影響較大,β2=90°時減振效果較好。通過遺傳算法優(yōu)化結果β2=89.63°,與結論基本一致。將幾種未優(yōu)選的聲子組合桿與優(yōu)化后聲子晶體組合桿傳遞率曲線對比見圖8。由圖8看出,優(yōu)化后聲子晶體組合桿減振性能整體優(yōu)于其它情況,尤其在1 500~1 700 Hz范圍內(nèi)減振效果加明顯。不同組合桿形式在1 500~2 500 Hz范圍內(nèi)平均減振幅值即綜合減振性能見圖9,幅值越小表示綜合減振性能越好。由圖9看出,優(yōu)選的組合桿形式綜合減振性能最好。由于影響聲子晶體組合桿減振性能因素較多,如桿徑比(子段A、B直徑比)、子桿長度比、彈性模量等,通過遺傳算法優(yōu)化多參數(shù)目標函數(shù),選出減振性能最優(yōu)的聲子晶體組合桿形式為研究重點。

圖7 目標函數(shù)收斂曲線

圖8 聲子晶體組合桿傳遞率對比

圖9 聲子晶體組合桿目標函數(shù)值對比

5 結 論

用傳遞矩陣法計算聲子晶體串聯(lián)組合桿振動帶隙,并用數(shù)值計算與有限元仿真對振動傳遞特性進行計算及分析,結論如下:

(1) 聲子晶體串聯(lián)組合桿具有三維寬頻減振能力,為良好的寬頻多維減振元件。其機理在于聲子晶體串聯(lián)組合桿可實現(xiàn)波型轉換,使彎曲、縱向、扭轉帶隙同時發(fā)揮作用,達到多維寬頻減振效果。

(2) 加載角對聲子晶體串聯(lián)組合桿減振性能有一定影響。應據(jù)振源頻率特點盡量選用0°或90°加載角。

(3) 桿件間夾角對減振性能影響明顯。實際應用中應綜合考慮各桿件間夾角,進行夾角最優(yōu)化求解,確定最佳角度以滿足不同減振需求。

[1] 溫熙森,溫激鴻,郁殿龍,等。聲子晶體[M].北京:國防工業(yè)出版社,2009.

[2] 宋卓斐,王自東,王艷林,等。一維桿狀聲子晶體的帶隙特性[J].振動與沖擊,2010,29(2):145-148.

SONG Zhuo-fei, WANG Zi-dong, WANG Yan-lin, et al. The characteristic of stop-band of a kind PC rod [J]. Journal of Vibration and Shock,2010,29(2):145-148.

[3] Goffaux C, Sanchez-Dehesa J. Two-dimensional phonon crystals studied using a variational method: application to lattices of locally resonant materials[J]. Physical Review B,2003, 67(14):144301-144301.

[4] Psarobas I E, Sigalas M M. Elastic band gaps in a fcc lattice of mercury spheres in aluminum[J]. Physical Review B,2002, 66(5):052302-052302.

[5] Liu Z Y, Chan C T, Sheng P. Analytic model of phononic crystals with local resonances[J]. Physical Review B,2005, 71(1):014103-014103.

[6] Lai Yun, Zhang Zhao-Qing. Large band gaps in elastic phononic crystals with air inclusions[J]. Applied Physics Letters, 2003,83(19):3900-3902.

[7] Martinsson P G, Movchan A B.Vibrations of lattice structures and phononic band gaps[J]. Q.J. Mech. Appl. Math., 2003,56:45-64.

[8] 郁殿龍,溫激鴻,劉耀宗.充液周期管路的軸向振動帶隙特性[J].機械工程學報,2009, 45(9):36-41.

YU Dian-long, WEN Ji-hong, LIU Yao-zong. Axial vibration property of periodic pipe system conveying fluid[J]. Journal of Mechanical Engineering,2009,45(9): 36-41.

[9] Wang Gang, Wen Xi-sen, Wen Ji-hong,et al. Quasi one-dimensional periodic structure with locally resonant band gap[J].ASME Journal of Applied Mechanics, 2006,73(1): 167-169.

[10] Shu Hai-sheng,Liu Shao-gang,Zhao Dan,et al. A combined sonic crystal rod of Bragg type with ability of vibration reduction in broad frequency band[J]. Noise Control Engineering Journal,2012,60(1): 42-61.

猜你喜歡
有限元振動
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
噴水推進高速艇尾部振動響應分析
This “Singing Highway”plays music
新型有機玻璃在站臺門的應用及有限元分析
基于有限元的深孔鏜削仿真及分析
基于有限元模型對踝模擬扭傷機制的探討
振動攪拌 震動創(chuàng)新
中國公路(2017年18期)2018-01-23 03:00:38
中立型Emden-Fowler微分方程的振動性
磨削淬硬殘余應力的有限元分析
UF6振動激發(fā)態(tài)分子的振動-振動馳豫
計算物理(2014年2期)2014-03-11 17:01:44
主站蜘蛛池模板: 久久国产拍爱| 国产高颜值露脸在线观看| 久久国产拍爱| 日韩a在线观看免费观看| 国产一区二区三区在线观看免费| 日韩二区三区| 国产精品人莉莉成在线播放| 亚洲中久无码永久在线观看软件| 欧美日韩精品一区二区视频| 久久亚洲高清国产| 免费jizz在线播放| 亚洲最猛黑人xxxx黑人猛交| 欧美午夜视频在线| 亚洲日韩精品欧美中文字幕| 波多野结衣无码视频在线观看| 欧美精品一区在线看| 欧美伊人色综合久久天天| 久久99国产乱子伦精品免| 91亚洲免费| 青青草原偷拍视频| 国产中文在线亚洲精品官网| 亚洲人成网线在线播放va| 日韩欧美中文字幕一本| 午夜久久影院| 日韩小视频在线观看| 国产视频大全| 国产综合亚洲欧洲区精品无码| 亚洲欧美日韩高清综合678| 国产综合另类小说色区色噜噜| 亚洲最大看欧美片网站地址| 特级aaaaaaaaa毛片免费视频 | 日韩第九页| 韩日无码在线不卡| 伊人丁香五月天久久综合 | 精品久久国产综合精麻豆| 国产超薄肉色丝袜网站| 久久中文字幕2021精品| 久久精品午夜视频| 很黄的网站在线观看| 国产伦精品一区二区三区视频优播 | 真实国产乱子伦视频| 99久久精品国产自免费| 久久婷婷六月| 国产情侣一区| 国产一区二区三区日韩精品| 日本免费一级视频| 日本欧美成人免费| 国产精品乱偷免费视频| 成人精品免费视频| AⅤ色综合久久天堂AV色综合| 2021国产精品自产拍在线| 日本国产一区在线观看| 在线无码av一区二区三区| 欧美激情第一区| 性网站在线观看| 国产乱子伦视频三区| 波多野结衣中文字幕一区二区| 国产精品美女免费视频大全| 亚洲经典在线中文字幕| 色精品视频| 乱系列中文字幕在线视频| 国产无码高清视频不卡| 黄色网站在线观看无码| 久久黄色小视频| 亚洲一级毛片在线播放| 国产精品太粉嫩高中在线观看 | 欧美另类一区| 国产一区在线视频观看| av一区二区无码在线| 91在线无码精品秘九色APP | 国产毛片高清一级国语 | 九色视频一区| 欧美亚洲欧美| 欧美精品在线免费| 99精品在线看| 99在线国产| 免费无码网站| 欧美a网站| 国产爽妇精品| 精品视频一区二区观看| 97视频在线观看免费视频| 91视频国产高清|