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

等截面一維聲子晶體角梁結(jié)構(gòu)振動特性

2021-07-19 09:56:44史冬巖何東澤
科學(xué)技術(shù)與工程 2021年17期
關(guān)鍵詞:振動結(jié)構(gòu)

張 穎,史冬巖,何東澤

(1.黑龍江八一農(nóng)墾大學(xué)工程學(xué)院,大慶 163319;2.哈爾濱工程大學(xué)機(jī)電工程學(xué)院,哈爾濱 150001)

伴隨著時(shí)代的進(jìn)步和社會的發(fā)展,機(jī)械制造業(yè)以及加工業(yè)對精密儀器的加工越來越重視,尤其對于機(jī)械加工制造而言,如何減振降噪已經(jīng)成為工程機(jī)械設(shè)計(jì)和生產(chǎn)的關(guān)鍵的因素。聲子晶體概念是類比于光子晶體的概念而提出的,通過對比與光子晶體,可以明顯地比較看出,彈性波傳播時(shí),會產(chǎn)生較為特殊的振動帶隙。正因?yàn)槁曌泳w結(jié)構(gòu)獨(dú)特的振動帶隙特性,是一個(gè)較為新穎的領(lǐng)域,其工程應(yīng)用價(jià)值也存在著較大的潛在價(jià)值。目前研究中,對于聲子晶體結(jié)構(gòu)振動特性的掌握,仍不足以滿足實(shí)際應(yīng)用需求,研究相關(guān)結(jié)構(gòu)減振特性具有重要意義。

1 研究現(xiàn)狀分析

聲子晶體結(jié)構(gòu)存在特殊的周期性結(jié)構(gòu),當(dāng)彈性波在其中傳播時(shí),會使其在傳播過程之中產(chǎn)生較為特殊的色散關(guān)系,即能帶關(guān)系[1]。1883年的Floquet較早對周期材料結(jié)構(gòu)進(jìn)行討論,進(jìn)行了關(guān)于一維Mathieus方程的研究。通常將振動衰減范圍內(nèi)的頻率定義為帶隙,進(jìn)而提出了聲子晶體的概念。近年來,基于光學(xué)以及固體物理學(xué)研究領(lǐng)域的擴(kuò)展,光子能帶理論和電子能帶理論對各類交叉學(xué)科領(lǐng)域產(chǎn)生了深遠(yuǎn)的影響,人們提出了一種新的物理概念——聲子晶體,借助于能帶理論中的晶格、倒格子、布里淵區(qū)以及帶隙等相關(guān)研究結(jié)果,對能帶結(jié)構(gòu)以及振動特性進(jìn)行計(jì)算[2]。為了進(jìn)一步對聲子晶體進(jìn)行研究,需要對聲子晶體進(jìn)行分類,通常可以將聲子晶體分為一維、二維以及三維,具體的分類規(guī)則是依據(jù)周期分布排列的方向數(shù)而確定的[3]。在工程實(shí)際的應(yīng)用之中,梁結(jié)構(gòu)較為常用,并且聲子晶體梁結(jié)構(gòu)由于其獨(dú)特性能而具有較為廣闊的應(yīng)用前景;在理論方面,一維聲子晶體結(jié)構(gòu)是一種散射體與基體相互不可區(qū)分的一維結(jié)構(gòu),所以其理論模型應(yīng)具有一定的特殊性,具有重要的研究價(jià)值。在工程應(yīng)用以及裝備生產(chǎn)的過程之中,一維聲子晶體結(jié)構(gòu)具備穩(wěn)定機(jī)械振動的作用,降低其在生產(chǎn)過程以及對工作環(huán)境產(chǎn)生的影響和危害[4]。尤其聲輻射、聲噪會對用戶及操作人員的身心健康產(chǎn)生危害[5]。在國防工業(yè)之中,減振更是值得考慮的問題,例如飛機(jī)、潛艇、航母等高科技裝備的制造和工作中,尤其是在惡劣的環(huán)境條件下都需要具有較好的減振特性,明顯的振動噪聲會造成較多不良后果[6]。

近年來很多學(xué)者針對一維聲子晶體結(jié)構(gòu)的振動特性開展相關(guān)研究。宿星亮等[7]研究了功能梯度材料復(fù)合而成的一維聲子晶體結(jié)構(gòu)的彈性波帶隙特征。左曙光等[8]考察了一維局域共振聲子晶體帶隙受材料黏彈性的影響。張法[9]采用理論推導(dǎo)、仿真分析及實(shí)驗(yàn)驗(yàn)證相結(jié)合的研究方法,分析和比較了聲子晶體組合桿振動特性。舒海生等[6]針對直梁低維振動特性的不足,提出一類聲子晶體角梁結(jié)構(gòu)的理論模型,并開展數(shù)值分析,提供該結(jié)構(gòu)減振的有效方法。朱學(xué)治等[10]基于歐拉梁理論,分析了有周期分布轉(zhuǎn)動振子的聲子晶體梁帶隙變化的一般規(guī)律。何東澤等[11]將一維聲子晶體振動特性的研究拓展到均勻變截面結(jié)構(gòu)當(dāng)中。涂靜等[12]更進(jìn)一步將振動帶隙特性研究拓展到雙層歐拉梁聲子晶體結(jié)構(gòu)中,指出雙層梁結(jié)構(gòu)在減振方面具有特有的優(yōu)勢。

而在以往的研究過程中,大多數(shù)的研究主要集中于單一梁的研究,對組合梁的振動特性的研究較少,等截面聲子晶體角梁可以視為工程應(yīng)用上較為常用的梁單元,是一種常見的組合梁。在工程應(yīng)用之中存在著較多的應(yīng)用,比如擔(dān)架的角梁扶手,房屋屋脊框架的結(jié)構(gòu)組成等均是由角梁單元變化而成[11]。基于此,現(xiàn)對一維等截面聲子晶體、組合桿件——角梁結(jié)構(gòu)的振動特性進(jìn)行研究分析,研究各參數(shù)對聲子晶體角梁振動特性的影響,對發(fā)揮其材料和結(jié)構(gòu)優(yōu)勢、有效控制振動和聲輻射具有重要意義。

2 角梁模型建立與描述

一維等截面聲子晶體角梁結(jié)構(gòu)是由多段等截面一維短粗梁在直角方向上的排列組合而成(圖1),本節(jié)主要研究三周期的一維等截面聲子晶體角梁,聲子晶體角梁的橫截面尺寸與軸向尺寸相當(dāng)。

圖1 一維等截面曲梁示意圖

選取材料A為有機(jī)玻璃,材料B為鋁,其性能參數(shù)如下:鋁:密度為2 799 kg/m3,彈性模量為7.21×1010Pa,泊松比為0.33;有機(jī)玻璃:密度為1 442 kg/m3,彈性模量為0.32×1010Pa,泊松比為0.33。

對于一維聲子晶體角梁結(jié)構(gòu),其物理參數(shù)如下:材料A的長度為l1=60 mm,材料B的長度為l2=60 mm,則周期常數(shù),即晶格常數(shù)N為l1+l2=120 mm;g為組分比,表示結(jié)構(gòu)中兩種材料的比例;材料A和材料B的橫截面均為矩形,橫截面面積相等,即A1=A2,具體參數(shù)如下:截面寬度b1=b2=10 mm,高度h1=h2=10 mm。定義角梁的夾角為θ,則本節(jié)所研究的角梁夾角θ為90o。

如圖2可以表示直角連接處,定義其節(jié)點(diǎn)數(shù)為2,則可以明確地表示出單元12和單元23中對偶坐標(biāo)系的建立情況,根據(jù)各自坐標(biāo)系的建立,結(jié)合節(jié)點(diǎn)處的廣義力平衡條件以及廣義位移協(xié)調(diào)條件進(jìn)行計(jì)算,得出整個(gè)系統(tǒng)的回傳射線矩陣,結(jié)合邊界條件進(jìn)行求解,得出相對應(yīng)的頻率相應(yīng)函數(shù)曲線,研究一維等截面聲子晶體角梁的振動特性。

圖2 坐標(biāo)定義

3 回傳射線矩陣法

回傳射線矩陣法的主要求解思路為:建立整體坐標(biāo)系,將整個(gè)結(jié)構(gòu)劃分為若干的單元,建立一對對偶的局部坐標(biāo)系,根據(jù)相應(yīng)的物理關(guān)系得出對應(yīng)的控制微分方程,進(jìn)行傅里葉變換,將時(shí)域范圍內(nèi)的偏微分方程轉(zhuǎn)化為頻域范圍內(nèi)的常微分方程,求出對應(yīng)的解[13]。將控制微分方程的解寫成波的形式,分為入射波和出射波兩大類;在同一單元之中,因一端的節(jié)點(diǎn)處的出射波,經(jīng)過在該單元之中的傳播過程,在另一單元之中將變?yōu)槿肷洳ǎ紤]二者之間的轉(zhuǎn)化關(guān)系;將各個(gè)單元的相位關(guān)系按照節(jié)點(diǎn)順序并且與散射關(guān)系中的順序相同組成總體相位關(guān)系矩陣;通過聯(lián)立總體散射關(guān)系矩陣以及總體相位關(guān)系矩陣,并且引入轉(zhuǎn)換矩陣,計(jì)算得出整體結(jié)構(gòu)的回傳射線矩陣,進(jìn)而求出相應(yīng)的入射波與出射波的波幅向量,求解對應(yīng)的物理量[14-16]。

3.1 結(jié)構(gòu)描述

假設(shè)一維聲子晶體由K段組成,將每一段結(jié)構(gòu)視為一個(gè)整體,建立整體坐標(biāo)系并進(jìn)行結(jié)構(gòu)定義。則可以將整個(gè)聲子晶體桿劃分為K個(gè)整體結(jié)構(gòu),即,將每一段材料視為一個(gè)整體結(jié)構(gòu),分別建立整體坐標(biāo)系和對偶坐標(biāo)系進(jìn)行求解。通常采用數(shù)字對節(jié)點(diǎn)以及單元進(jìn)行編號,如圖3所示為對聲子晶體桿進(jìn)行節(jié)點(diǎn)分配以及結(jié)構(gòu)定義,其中,K=6。

圖3 結(jié)構(gòu)描述示意圖

3.2 總體坐標(biāo)與單元局部坐標(biāo)

對于一根桿單元,如i、j,分別引入對偶的局部坐標(biāo)系(x,y,z)ij和(x,y,z)ji,分別以i和j為相應(yīng)的坐標(biāo)原點(diǎn),以單元的軸線方向?yàn)閤軸方向。對于(x,y,z)ij坐標(biāo)系而言,以i到j(luò)的方向?yàn)檎较颍鄬?yīng)的,(x,y,z)ji以j到i的方向?yàn)檎较颉ij和xji在xy平面內(nèi)逆時(shí)針旋轉(zhuǎn)90°即可得到y(tǒng)ij和yji,根據(jù)各自坐標(biāo)系中x和y的正方向,結(jié)合右手螺旋法則,得到對應(yīng)的zij和zji。對偶坐標(biāo)的具體形式如圖4所示。從圖4中可以看出,xij和xji方向相反,yij和yji方向相反,zij和zji方向相同。在對偶坐標(biāo)系中,坐標(biāo)x滿足關(guān)系式:xij=lij-xji以軸向位移u為例,可以得到uij=-uji。

圖4 單元的對偶坐標(biāo)系

3.3 坐標(biāo)變換與對偶變換

圖5 局部坐標(biāo)系與整體坐標(biāo)系之間的關(guān)系

(1)

對于任意單元i、j而言,將每一個(gè)物理量的上標(biāo)用對應(yīng)的坐標(biāo)系標(biāo)注,以區(qū)分不同坐標(biāo)系中的物理量。對于每個(gè)單元,單元中兩個(gè)節(jié)點(diǎn)之間滿足對偶變換關(guān)系,各物理量之間的關(guān)系為

Nij=Nji,Qij=Qji,Mij=-Mji,uij=-uji,

vij=-vji,φij=φji

(2)

式(2)中:Nij、Qij、Mij分別為梁的軸向力、剪切力以及彎矩;uij、vij、φij分別為3個(gè)方向的位移分量。

將單元坐標(biāo)系中的廣義力與廣義位移之間的關(guān)系表示為

(3)

將上述關(guān)系矩陣分別定義為廣義力坐標(biāo)變換矩陣,廣義位移坐標(biāo)變換矩陣為

(4)

3.4 行波解及其矩陣形式

現(xiàn)有的梁理論主要是Euler-Bernoulli梁理論和Timoshenko梁理論。其中在實(shí)際應(yīng)用過程中發(fā)現(xiàn),Euler-Bernoulli梁理論的結(jié)果只有當(dāng)細(xì)長梁的低頻情況下才和理論值相接近,當(dāng)在高頻區(qū)段或者梁模型為短粗梁時(shí),其理論計(jì)算結(jié)果與實(shí)際結(jié)果相差甚遠(yuǎn)[15]。隨著時(shí)代的進(jìn)步以及科學(xué)情況的發(fā)展,人們提出了Timoshenko梁理論,該理論不僅僅考慮到轉(zhuǎn)動慣量的影響,并且還考慮到梁的剪切變形。這個(gè)理論極大地改變了梁的動力學(xué)理論,提高了計(jì)算的精度,具有極大的工程應(yīng)用價(jià)值[17]。對于梁的軸向波的控制微分方程,Timoshenko梁理論與經(jīng)典的梁理論是相同的,在笛卡爾坐標(biāo)系中,考慮梁的彎曲變形之中剪切變形和轉(zhuǎn)動慣量的影響,其控制微分方程[5-6]為

(5)

式(5)中:E(x)、G(x)分別為材料的拉伸和剪切彈性模量系數(shù);ρ(x)為材料的密度;A(x)為橫截面面積;I(x)為橫截面二次矩;Q(x,t)為梁的剪力;M(x,t)為梁的彎矩;v(x,t)為梁的撓度位移;φ(x,t)為梁的轉(zhuǎn)角位移。

對式(5)進(jìn)行Fourier變化,可得

(6)

則微分方程的通解可以表示為

(7)

式(7)中:ω為振動頻率;a2(ω)、a3(ω)為入射波解向量;d2(ω)、d3(ω)為入射波解向量;k2、k3為彎曲波的波數(shù);β2、β3為比例系數(shù),表達(dá)式為

(8)

(9)

式(9)中:ξ2、ξ3、μ2、μ3為比例系數(shù),其具體定義為

(10)

對于梁單元的狀態(tài)方程,可以得出其狀態(tài)向量,并且可以得到狀態(tài)方程的解,將狀態(tài)方程的頻域解寫成矩陣形式,可以得到形式為

(11)

式(11)中:AF、DF為力傳播矩陣;AV、DV為位移傳播矩陣;k1為縱波波數(shù),k1=ω/c1。上述物理量的具體矩陣形式可以表示為

(12)

3.5 單元分析與相位關(guān)系

根據(jù)單元內(nèi)對偶變換以及公式推導(dǎo),可以得到關(guān)系式為

(13)

經(jīng)過定義可以得到對應(yīng)的相位關(guān)系為

(14)

式(14)中:Pij為單元i、j的相位矩陣。

(15)

由上文所得節(jié)點(diǎn)相位關(guān)系,將所有節(jié)點(diǎn)的局部相位矩陣按照節(jié)點(diǎn)編號的順序組成總體相位矩陣為

(16)

式(16)中:P為相應(yīng)的總體相位關(guān)系矩陣,是6N×6N階的對角陣。單元總體的入射波向量和出射波向量可以分別表示為

(17)

引入轉(zhuǎn)換矩陣,得出整體散射矩陣為

(18)

式(18)中:U0為轉(zhuǎn)換矩陣,其每一行、每一列有且只有一個(gè)元素為1,可以得到整體散射矩陣:a=PUd。

3.6 節(jié)點(diǎn)分析與色散關(guān)系

根據(jù)廣義力平衡關(guān)系,以節(jié)點(diǎn)i為例,設(shè)節(jié)點(diǎn)i的相關(guān)單元數(shù)為ni,可以寫出在頻域內(nèi)的平衡方程[7]為

(19)

式(19)中:M、K、u、f為相應(yīng)單元的質(zhì)量矩陣、剛度矩陣、位移向量以及外載荷向量。

將式(12)代入式(19)中,得

fi=(-ω2Mi+Ki)ui

(20)

Tij[Vij(0,ω)]=ui,j=1,2,…,nj

(21)

將式(21)代入式(20)中,可得

(22)

di=Siai+si

(23)

式(23)中:Si為局部散射矩陣;si為外載矩陣。

(24)

同時(shí),按照節(jié)點(diǎn)順序得出相應(yīng)的散射矩陣S[8],即

d=Sa+s

(25)

根據(jù)結(jié)構(gòu)受到的周期載荷,如位移u(t)=u(ω)eiωt,作用時(shí)產(chǎn)生的動力響應(yīng)成為穩(wěn)態(tài)響應(yīng),分別定義力向量和位移向量為

(26)

進(jìn)而可得

(27)

式(27)中:AFO和DFO為整體廣義位移的傳播矩陣;AVO和DVO為整體廣義力的傳播矩陣。

4 數(shù)值計(jì)算與比較

4.1 理論計(jì)算與有限元仿真

為了減小與一維聲子晶體模型與實(shí)際工程應(yīng)用的差別,取一維等截面聲子晶體角梁的左端面為激勵(lì)輸入段,取角梁的最右端面為響應(yīng)拾取端,并且取自由邊界條件進(jìn)行理論計(jì)算。于激勵(lì)輸入端輸入振幅為1 mm的軸向位移激擾,取右端的輸出位移響應(yīng),進(jìn)行數(shù)據(jù)處理。

對一維等截面聲子晶體角梁結(jié)構(gòu)進(jìn)行有限元仿真計(jì)算,需要對其進(jìn)行設(shè)置,其具體參數(shù)如下:①材料參數(shù)與周期結(jié)構(gòu)物理參數(shù):與前文結(jié)構(gòu)描述相同;②離散化處理:對整體模型進(jìn)行自由網(wǎng)格劃分的方法,將每一個(gè)離散化單元的尺寸定義為2 mm;③施加載荷:主要研究一維等截面聲子晶體對軸向波的傳輸特性,因此施加軸向簡諧位移載荷,位移的幅值為1 mm,激振頻率范圍為0~30 kHz;④分析模塊設(shè)置:采用諧響應(yīng)分析模塊,因施加的是位移載荷,因此采用完全法進(jìn)行有限元模擬仿真計(jì)算,具體如圖6所示。

圖6 有限元模擬仿真示意圖

為了進(jìn)行回傳射線矩陣法計(jì)算得出的頻率響應(yīng)曲線與有限元分析軟件計(jì)算的曲線對比,由圖7中可以看出,二者重合程度較為良好,可以證明由回傳射線矩陣法計(jì)算的一維等截面聲子晶體角梁結(jié)構(gòu)得出的頻率響應(yīng)函數(shù)曲線的正確性。可以看出,在0~30 kHz的頻率范圍內(nèi),一維等截面聲子晶體角梁存在著較多的振動帶隙,并且在相應(yīng)的帶隙范圍內(nèi)存在表面局域態(tài)現(xiàn)象,在各自的振動帶隙的范圍內(nèi),最大衰減程度較少,對輸入激勵(lì)的傳播抑制程度較低。

圖7 有限元仿真與數(shù)值計(jì)算的對比

4.2 晶格常數(shù)N的影響

為了分析研究晶格常數(shù)N對一維有限尺寸等截面聲子晶體角梁振動特性的影響,將不同晶格常數(shù)N對應(yīng)的頻率響應(yīng)曲線進(jìn)行比較分析研究。取N的范圍為100~140 mm,得出各自晶格常數(shù)N對應(yīng)的頻率響應(yīng)曲線。通過圖8中各條曲線可以觀察出,在0~35 kHz的頻率范圍內(nèi),隨著晶格常數(shù)N的增加,頻率響應(yīng)曲線對應(yīng)的振動帶隙以及帶隙范圍內(nèi)的最大衰減程度變化規(guī)律不是特別明顯,但是可以比較出隨著晶格常數(shù)N的增加,0~35 kHz頻率范圍內(nèi)的帶隙數(shù)量增加。

圖8 不同晶格常數(shù)N對應(yīng)的頻響曲線

4.3 組分比g的影響

為了比較分析研究單一組成材料的長度對一維聲子晶體角梁振動特性的影響,改變材料A,即材料有機(jī)玻璃的長度,將其長度定義為l1,得到不同的組分比,通過組分比的變化,得出相對應(yīng)的頻率響應(yīng)函數(shù)曲線。在圖9中可以看出,在0~35 kHz的頻率范圍內(nèi),隨著組分比u的變化,各曲線對應(yīng)表達(dá)的振動特性不是十分明顯。

圖9 不同組分比u對應(yīng)的頻響曲線

4.4 角梁夾角角度θ的影響

為了研究比較分析一維聲子晶體角梁的夾角θ對振動特性的影響,通過改變夾角θ的大小,得到夾角θ分別為60°、90°和120°所對應(yīng)的頻率響應(yīng)函數(shù)曲線。通過圖10可以看出,隨著一維等截面聲子晶體角梁夾角θ的變化,各個(gè)夾角所對應(yīng)的曲線基本吻合,變化的幅度不是很大。因此,角梁夾角θ對一維等截面聲子晶體角梁的振動特性影響不大。

圖10 不同夾角θ對應(yīng)的頻響曲線

4.5 橫截面形狀的影響

分別將橫截面的形狀參數(shù)設(shè)計(jì)為:圓形橫截面r1=5 mm,矩形橫截面b1=10 mm,h1=8 mm;正方形橫截面b2=h2=10 mm。通過圖11可以看出,三種橫截面所對應(yīng)的頻率相應(yīng)函數(shù)曲線中,各個(gè)橫截面形狀對應(yīng)的頻率響應(yīng)函數(shù)曲線在0~35 kHz的頻率范圍內(nèi),其帶隙寬度,帶隙的起始頻率,截止頻率以及帶隙范圍內(nèi)的最大衰減程度基本上保持不變,因此,橫截面形狀對一維等截面聲子晶體角梁的振動特性的影響較小。

圖11 不同橫截面形狀對應(yīng)的頻響曲線

5 結(jié)論

基于Timoshenko梁理論,以及回傳射線矩陣法對等截面一維聲子晶體角梁結(jié)構(gòu)對于軸向波的振動特性進(jìn)行研究,對角梁形式的一維聲子晶體角梁結(jié)構(gòu)的振動特性進(jìn)行計(jì)算,并選取有限元分析法作為參考,證明了算法的正確性。通過相應(yīng)的頻率響應(yīng)函數(shù)曲線之間的不同,比較分析了各參數(shù)對聲子晶體角梁振動特性的影響及相應(yīng)規(guī)律。

猜你喜歡
振動結(jié)構(gòu)
振動的思考
噴水推進(jìn)高速艇尾部振動響應(yīng)分析
《形而上學(xué)》△卷的結(jié)構(gòu)和位置
This “Singing Highway”plays music
論結(jié)構(gòu)
中華詩詞(2019年7期)2019-11-25 01:43:04
新型平衡塊結(jié)構(gòu)的應(yīng)用
模具制造(2019年3期)2019-06-06 02:10:54
振動攪拌 震動創(chuàng)新
中國公路(2017年18期)2018-01-23 03:00:38
中立型Emden-Fowler微分方程的振動性
論《日出》的結(jié)構(gòu)
創(chuàng)新治理結(jié)構(gòu)促進(jìn)中小企業(yè)持續(xù)成長
主站蜘蛛池模板: 午夜色综合| 久久精品亚洲热综合一区二区| 欧美一区二区自偷自拍视频| 亚洲精品无码久久毛片波多野吉| 成年A级毛片| 亚洲国产看片基地久久1024| 四虎成人精品在永久免费| 欧美成a人片在线观看| 91欧美在线| 无码日韩视频| 小说 亚洲 无码 精品| 九九九精品成人免费视频7| 亚洲三级电影在线播放 | 亚洲人成色77777在线观看| 在线综合亚洲欧美网站| 欧美日韩导航| 色婷婷在线播放| 久久中文电影| 国产视频一区二区在线观看| 欧美va亚洲va香蕉在线| 国产在线观看一区二区三区| 免费看久久精品99| 欧美另类一区| 午夜日本永久乱码免费播放片| 亚洲91精品视频| 久久不卡国产精品无码| 97国产在线播放| 欧美色伊人| 国产在线观看一区精品| 六月婷婷综合| 国产精品成人免费综合| 久久一色本道亚洲| 成人蜜桃网| 国产毛片不卡| 欧美一区日韩一区中文字幕页| 自拍偷拍欧美日韩| 一区二区三区精品视频在线观看| 日本精品影院| 久久久久久久97| 国产欧美视频在线观看| 日韩成人在线一区二区| 91在线免费公开视频| 狂欢视频在线观看不卡| 中文字幕欧美日韩高清| 色精品视频| 国产人在线成免费视频| 亚瑟天堂久久一区二区影院| 国产va在线观看免费| 精品第一国产综合精品Aⅴ| 国产精品美女免费视频大全 | 国产xx在线观看| 国产美女一级毛片| 毛片免费在线视频| 亚洲男人天堂2018| 少妇被粗大的猛烈进出免费视频| 国产成人乱无码视频| 亚洲狼网站狼狼鲁亚洲下载| 日韩小视频在线观看| 91免费国产在线观看尤物| 国产福利影院在线观看| 日本三级精品| 欧美亚洲网| 亚洲欧美人成电影在线观看| 国产日韩欧美黄色片免费观看| 制服丝袜一区二区三区在线| 亚洲综合18p| 视频国产精品丝袜第一页| 日本欧美成人免费| 国产91透明丝袜美腿在线| 毛片免费试看| 欧美一级夜夜爽www| 亚洲日韩精品伊甸| 久久a毛片| 丰满的熟女一区二区三区l| 国产精品视频白浆免费视频| 狠狠ⅴ日韩v欧美v天堂| 久热中文字幕在线| 欧美日韩久久综合| 久久久久青草大香线综合精品 | 天天综合网亚洲网站| 欧美日韩亚洲综合在线观看| 男女猛烈无遮挡午夜视频|