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

        ?

        預(yù)測(cè)銑削穩(wěn)定性的Hamming線性多步法

        2018-11-30 01:05:02智紅英閆獻(xiàn)國曹啟超
        振動(dòng)與沖擊 2018年22期
        關(guān)鍵詞:模態(tài)方法模型

        智紅英 , 閆獻(xiàn)國, 杜 娟, 曹啟超

        (1.太原科技大學(xué) 應(yīng)用科學(xué)學(xué)院, 太原 030024; 2.太原科技大學(xué) 機(jī)械工程學(xué)院, 太原 030024)

        在高速精密銑削加工過程中,再生顫振是一種常見的自激振動(dòng)現(xiàn)象[1-2],是引起加工過程失穩(wěn)的主要因素[3]。因此,為了防止銑削加工過程中再生顫振的發(fā)生,加工前需要對(duì)銑削穩(wěn)定性進(jìn)行預(yù)測(cè),以便在加工時(shí)獲得較高的表面質(zhì)量,從而提高零件的精度。此外,考慮到再生顫振的銑削動(dòng)力學(xué)模型可以近似地表示為時(shí)滯微分方程[4-5],所以如何求解時(shí)滯微分方程才能精確且快速地獲得穩(wěn)定性葉瓣圖是預(yù)測(cè)銑削穩(wěn)定性的關(guān)鍵。

        迄今為止,預(yù)測(cè)銑削穩(wěn)定性的方法有很多種,主要分為實(shí)驗(yàn)方法和數(shù)值方法,而數(shù)值方法主要有頻域法和時(shí)域法兩種。Budak等[6-7]提出了一種近似零階的半解析方法的頻域法求解銑削動(dòng)力學(xué)方程,但這種方法無法滿足小徑向切深切削的計(jì)算精度。為了解決這個(gè)問題,Merdol等[8]提出了一種可以滿足小徑向切深的多頻域法。隨后,Insperger等[9-11]用時(shí)域法預(yù)測(cè)顫振穩(wěn)定性,即半離散法(SDM)。該方法對(duì)銑削時(shí)滯微分方程中延遲項(xiàng)進(jìn)行離散,其計(jì)算精度取決于離散步長。

        基于這種方法,龍新華等人考慮了刀具與工件之間接觸效應(yīng)的損失[12]和切削厚度變化相關(guān)的時(shí)間延遲效應(yīng)[13],李中偉等[14]提出了基于Magnus-Gaussian 截?cái)嗟你娤飨到y(tǒng)穩(wěn)定性的改進(jìn)的半離散法,Magnus-Gaussian截?cái)喾ū攘汶A半離散法的收斂速度更快、計(jì)算時(shí)間更短。Butcher等[15-16]提出了Chebyshev多項(xiàng)式法和基于Chebyshev配點(diǎn)法來預(yù)測(cè)具有線性時(shí)滯微分方程銑削過程穩(wěn)定性。為了提高計(jì)算精度和計(jì)算效率,Ding等[17-18]提出了全離散方法(FDM),包括一階全離散法(1st-FDM)和二階全離散法(2nd-FDM),這種方法是通過線性插值對(duì)狀態(tài)項(xiàng)和延遲項(xiàng)進(jìn)行離散化處理。隨后,Guo等[19]在此基礎(chǔ)上通過牛頓插值提出了三階全離散法(3rd-FDM),得出了具有比一階、二階全離散法更高的計(jì)算精度和更快的計(jì)算效率。Li等[20]提出了一種完全離散算法,這種方法是針對(duì)銑削穩(wěn)定性預(yù)測(cè)的一種半解析算法。與半離散算法(SDM)和全離散算法(FDM)相比,此算法對(duì)時(shí)滯微分方程的各個(gè)部分都進(jìn)行了離散化處理,簡(jiǎn)化了離散化后迭代方程的復(fù)雜度。為了精確分離刀具自由振動(dòng)時(shí)段和刀具-工件接觸過程中的強(qiáng)迫振動(dòng)時(shí)段,Bayly等[21]提出了時(shí)域有限元分析方法,并具有較高的計(jì)算精度。Ding等[22]提出了一種數(shù)值積分法,這種方法是基于牛頓-柯特斯公式和高斯公式,通過直接離散積分方程中的積分項(xiàng),構(gòu)造出系統(tǒng)在單周期上的Floquet轉(zhuǎn)移矩陣用于預(yù)測(cè)銑削穩(wěn)定性。Ding等[23]隨后又將有關(guān)數(shù)值積分法擴(kuò)展為譜方法。在數(shù)值積分方法思想的指導(dǎo)下,Niu等[24]提出了一種廣義的龍格庫塔方法預(yù)測(cè)具有再生顫振的銑削穩(wěn)定性。Li等[25]在經(jīng)典的四階龍格庫塔的基礎(chǔ)上,基于完全離散法,提出了一種完全離散化的龍格庫塔方法預(yù)測(cè)銑削穩(wěn)定性。Zhang等[26]基于辛普森公式提出了一種全新的方法預(yù)測(cè)銑削穩(wěn)定性。Qin等[27]提出了一種基于亞當(dāng)斯-莫爾頓的線性多步法預(yù)測(cè)銑削穩(wěn)定性。Zhang等[28]根據(jù)泰勒公式,運(yùn)用有限差分法和外推法,提出了一種數(shù)值微分法預(yù)測(cè)銑削穩(wěn)定性。

        隨著研究的不斷深入,越來越多的方法被用于預(yù)測(cè)銑削過程的穩(wěn)定性。針對(duì)銑削加工動(dòng)力學(xué)方程,提出了一種Hamming線性多步法預(yù)測(cè)銑削穩(wěn)定性。

        2 動(dòng)力學(xué)模型

        通過考慮再生顫振的動(dòng)力學(xué)模型預(yù)測(cè)銑削穩(wěn)定性,其動(dòng)力學(xué)模型可由以下時(shí)滯微分方程描述:

        (1)

        式中:M,C和K分別表示刀具的模態(tài)質(zhì)量、模態(tài)阻尼和模態(tài)剛度矩陣;q(t)為刀具模態(tài)坐標(biāo),Kc(t)為周期系數(shù)矩陣,ap為軸向切削深度,T為時(shí)滯量且等于刀齒切削周期,T=60/(NΩ),N為刀具齒數(shù),Ω為主軸轉(zhuǎn)速,單位為r/min。

        (2)

        式(2)中,A0表示常數(shù)矩陣;A(t)表示周期為T的考慮再生效應(yīng)的系數(shù)矩陣,且A(t)=A(t+T)。其中

        (3)

        單自由度和兩自由度系統(tǒng)主要區(qū)別是狀態(tài)向量x(t)、常數(shù)矩陣A0和周期矩陣A(t)的不同,下面將主要介紹這兩種動(dòng)力學(xué)模型。

        1.1 單自由度動(dòng)力學(xué)模型

        具有x方向的單自由度銑削加工動(dòng)力學(xué)模型可由下列方程表示:

        (4)

        式中:mt,ζ,ωn分別表示刀具的模態(tài)質(zhì)量、刀具的阻尼比、自然圓頻率;ap為軸向切削深度;T為時(shí)滯量且等于刀齒切削周期,即T=60/(NΩ)。h(t)由下列方程表示:

        (5)

        (6)

        式中:φst和φex分別表示刀具的切入角和切出角。對(duì)于順銑,φst=arccos(2a/D-1),φex=π; 逆銑時(shí),φst=0,φex=arccos(1-2a/D), 其中a/D表示為徑向切深與刀具直徑之比。

        (7)

        式(7)中

        (8)

        1.2 雙自由度動(dòng)力學(xué)模型

        在x和y方向上的雙自由度銑削動(dòng)力學(xué)模型可由下列方程表示:

        (9)

        式(9)中周期系數(shù)矩陣Kc(t)可表示為:

        (10)

        其中:

        (11)

        (12)

        (13)

        (14)

        (15)

        其中:

        (16)

        (17)

        2 Hamming線性多步法(HAMM)

        假設(shè)初始時(shí)刻為t0,通過空間狀態(tài)變換理論[17-18],式(1)的解可由如下方程表示:

        (18)

        刀齒切削周期T可以劃分為自由振動(dòng)時(shí)間間隔和tf強(qiáng)迫振動(dòng)時(shí)間間隔T-tf。當(dāng)?shù)毒咛幱谧杂烧駝?dòng)時(shí)刻t時(shí),即t∈[t0,t0+tf],狀態(tài)值有如下關(guān)系:

        x(t)=eA0(t-t0)x(t0)

        (19)

        加工時(shí)刀具處于強(qiáng)迫振動(dòng)時(shí)刻t時(shí),即t∈[t0+tf,T], 將切削時(shí)間T-tf平均分成m個(gè)時(shí)間間隔,則每個(gè)時(shí)間間隔可表示為h=(T-tf)/m,相應(yīng)的離散點(diǎn)表示為

        ti=t0+tf+(i-1)h,i=1, 2, …,m+1

        (20)

        當(dāng)t∈[ti,ti+1]時(shí),方程(18)可以轉(zhuǎn)化為如下表達(dá)式:

        (21)

        t=t1時(shí),代入(19)式可得狀態(tài)量x(t1)和時(shí)滯量x(tm+1-T)之間的關(guān)系如下式表示:

        x(t1)=x(t0+tf)=eA0·tfx(t0)=eA0·tfx(tm+1-T)

        (22)

        其余離散點(diǎn)x(ti)(i=2, …,m+1)可以通過線性多步法進(jìn)行求解,為了簡(jiǎn)化表達(dá)式,將A(ti)表示為Ai,x(ti)表示為xi,x(ti-T)表示為xi-T。

        x2和x3利用Adams線性多步法[27]可以表示為:

        (23)

        (24)

        對(duì)式(23)、(24)分別化簡(jiǎn)可得:

        (25)

        (26)

        對(duì)于xi(i=4, 5, …,m+1),用線性多步法中的Hamming公式來進(jìn)行求解,則可表示為:

        (27)

        式(27)整理可得:

        (28)

        聯(lián)立式(22), (25), (26), (28)可得:

        (29)

        其中:

        (30)

        (31)

        通過Hamming方法(HAMM)求得系統(tǒng)的狀態(tài)傳遞矩陣Φ為:

        Φ=P-1Q

        (32)

        根據(jù)Floquet理論,當(dāng)系統(tǒng)傳遞矩陣Φ特征值都小于1時(shí),系統(tǒng)處于穩(wěn)定狀態(tài),否則,該系統(tǒng)則會(huì)進(jìn)入顫振狀態(tài)。

        3 數(shù)值仿真結(jié)果分析

        3.1 收斂性分析

        局部離散誤差的變化可以反映所提出方法收斂速度的快慢程度。通過理論上對(duì)局部離散誤差的分析,1st-SDM[10]和2nd-FDM[18]的局部離散誤差均為O(h3),HAMM的局部離散誤差為O(h5)。利用Matlab軟件編程,在單自由度銑削動(dòng)力學(xué)模型下對(duì)1st-SDM、2nd-FDM和HAMM三種方法的收斂性進(jìn)行分析。切削參數(shù)為:徑向比a/D=1,切向力系數(shù)Kt=6×108N/m2,法向力系數(shù)Kn=2×108N/m2, 主軸轉(zhuǎn)速Ω=10 000 r/min和5 000 r/min,軸向切削深度ap==1.5 mm、1 mm、0.3 mm。模態(tài)參數(shù)采用與文獻(xiàn)[10]一樣的參數(shù):逆銑,刀具齒數(shù)N=2,模態(tài)質(zhì)量mt=0.039 93 kg,固有圓頻率wn=922×2π rad/s,固有阻尼ζ=0.011。

        在1st-SDM、2nd-FDM和HAMM中,分別計(jì)算一個(gè)刀齒周期下不同離散數(shù)m所對(duì)應(yīng)的狀態(tài)傳遞矩陣臨界特征值的模|μ|,局部離散誤差||μ|-|μ0||,而精確值|μ0|是m=500時(shí)采用1st-SDM得到的特征值的模。不同轉(zhuǎn)速和不同軸向切削深度下,1st-SDM、2nd-FDM、HAMM三種方法所對(duì)應(yīng)的局部離散誤差如圖1所示。

        由圖1可知,在不同的主軸轉(zhuǎn)速和軸向切削深度下,隨著離散數(shù)m的不斷增加,三種方法的局部離散誤差越來越接近于0。在離散數(shù)相同的條件下,HAMM的局部離散誤差要遠(yuǎn)遠(yuǎn)小于1st-SDM和2nd-FDM。

        比如主軸轉(zhuǎn)速Ω=10 000 r/min,|μ0|=1.336 4 r/min,m=50時(shí),1st-SDM的局部離散誤差為0.011 72,2nd-FDM的局部離散誤差為0.005 16,HAMM的局部截?cái)嗾`差為0.001 28,1st-SDM的局部離散誤差減少89%,比2nd-FDM的局部離散誤差減少75%。

        此外還可以明顯地看出,HAMM的收斂速度均快于其它兩種方法,主軸轉(zhuǎn)速越高,離散數(shù)少的HAMM能達(dá)到離散數(shù)較多的1st-SDM和2nd-FDM的局部離散誤差,和理論分析結(jié)論一致。因此,在相同的局部離散誤差下,HAMM計(jì)算所用的時(shí)間遠(yuǎn)遠(yuǎn)少于1st-SDM和2nd-FDM。

        圖1 1st-SDM、2nd-FDM、HAMM三種方法收斂速度的比較Fig.1 Convergence rate comparisons of1st-SDM, 2nd-FDM and HAMM

        3.2 穩(wěn)定性葉瓣圖

        3.2.1 單自由度銑削模型

        在單自由度銑削模型中,為了比較1st-SDM、2nd-FDM和HAMM這三種計(jì)算方法的準(zhǔn)確性,將模態(tài)參數(shù)設(shè)置與2.1節(jié)相同。徑向比a/D=1, 0.1, 0.05,主軸轉(zhuǎn)速為5 000≤Ω≤10 000 r/min,并設(shè)置200個(gè)等間距的轉(zhuǎn)速;軸向切削深度0≤w≤10 mm,并設(shè)置100個(gè)等間距切削深度,則銑削穩(wěn)定性葉瓣圖便構(gòu)成了200×100的網(wǎng)格。設(shè)離散數(shù)m為30、40,通過2.1節(jié)可知,HAMM計(jì)算速度快,具有較高階的誤差,故將HAMM離散數(shù)m=200時(shí)所得的穩(wěn)定性葉瓣圖作為理想?yún)⒖记€。在不同離散數(shù)值和不同徑深下,三種方法所得的穩(wěn)定性葉瓣圖與理想?yún)⒖记€的對(duì)比如圖2~4所示。

        由圖2~4可知,在單自由度銑削模型下,在同種方法中,隨著離散數(shù)m的增加,三種方法的計(jì)算精度都有所提高,都更接近理想?yún)⒖记€,當(dāng)離散數(shù)m相同時(shí),HAMM的計(jì)算精度要優(yōu)于1st-SDM和2nd-FDM的計(jì)算精度。

        圖2 單自由度銑削模型a/D=1下1st-SDM、2nd-FDM、HAMM三種方法的穩(wěn)定性葉瓣圖Fig.2 Stability lobe diagrams of1st-SDM, 2nd-FDM and HAMM for a single degree-of-freedom milling model with a/D=1

        圖3 單自由度銑削模型a/D=0.1下1st-SDM、2nd-FDM、HAMM三種方法的穩(wěn)定性葉瓣圖Fig.3 Stability lobe diagrams of1st-SDM, 2nd-FDM and HAMM for a single degree-of-freedom milling model with a/D=0.1

        圖4 單自由度銑削模型a/D=0.05下1st-SDM、2nd-FDM、HAMM三種方法的穩(wěn)定性葉瓣圖Fig.4 Stability lobe diagrams of1st-SDM, 2nd-FDM and HAMM for a single degree-of-freedom milling model with a/D=0.05

        在圖2中,當(dāng)m=30時(shí),1st-SDM的計(jì)算時(shí)間為69.36 s,2nd-FDM的計(jì)算時(shí)間為42.9 s,HAMM的計(jì)算時(shí)間為29.5 s,比1st-SDM所用時(shí)間縮短57.6%,比2nd-FDM所用時(shí)間縮短31.2%。當(dāng)m=40時(shí),1st-SDM的計(jì)算時(shí)間為117.1 s,2nd-FDM的計(jì)算時(shí)間為80 s,HAMM的計(jì)算時(shí)間為46.4 s,比1st-SDM所用時(shí)間縮短60.3%,比2nd-FDM所用時(shí)間縮短42%。

        在圖3,圖4中可以看出,在離散數(shù)m相同時(shí),當(dāng)徑向比a/D比較小時(shí),HAMM方法計(jì)算精度遠(yuǎn)遠(yuǎn)高于1st-SDM和2nd-FDM。故HAMM是一種較優(yōu)的預(yù)測(cè)銑削穩(wěn)定性的計(jì)算方法。

        3.2.2 雙自由度銑削模型

        雙自由度銑削過程的系統(tǒng)參數(shù)與單自由度的參數(shù)相同。設(shè)刀具浸入比a/D=0.1,離散數(shù)m=10和20,主軸轉(zhuǎn)速Ω從5 000 r/min到20 000 r/min,軸向切削深度ap從0~10 mm,穩(wěn)定性葉瓣圖的網(wǎng)格設(shè)置為200×100,將m=200時(shí)由HAMM進(jìn)行數(shù)值計(jì)算得出的穩(wěn)定性葉瓣圖作為理想?yún)⒖记€,三種方法的穩(wěn)定性葉瓣圖與理想?yún)⒖记€的對(duì)比圖,如圖5所示。

        圖5 雙自由度銑削模型下1st-SDM、2nd-FDM、HAMM三種方法的穩(wěn)定性葉瓣圖Fig.5 Stability lobe diagrams of1st-SDM, 2nd-FDM and HAMM for a double degree-of-freedom milling model

        由圖5可以看出,在同種方法中,隨著離散數(shù)的增加,三種方法的計(jì)算精度都有所提高;當(dāng)離散數(shù)m相同時(shí),HAMM的計(jì)算精度遠(yuǎn)遠(yuǎn)高于1st-SDM和2nd-FDM。舉個(gè)例子來說,當(dāng)m=10時(shí),1st-SDM的穩(wěn)定性葉瓣圖和理想?yún)⒖记€間產(chǎn)生了明顯的偏差,2nd-FDM的穩(wěn)定性葉瓣圖雖然降低了和參考曲線間的偏差,但還是存在明顯的不同之處,而HAMM的穩(wěn)定性葉瓣圖幾乎和理想?yún)⒖记€一致。

        當(dāng)m=10時(shí),1st-SDM的計(jì)算時(shí)間為32.9 s,2nd-FDM的計(jì)算時(shí)間為11.34 s,HAMM的計(jì)算時(shí)間為10.51 s。當(dāng)m=20時(shí),1st-SDM的計(jì)算時(shí)間為81.64 s,2nd-FDM的計(jì)算時(shí)間為38.54 s,HAMM的計(jì)算時(shí)間為32.16 s,因此,在離散數(shù)相同時(shí),HAMM的計(jì)算效率和2nd-FDM相當(dāng),遠(yuǎn)遠(yuǎn)高于1st-SDM。故HAMM是一種具有較高的計(jì)算精度和較快的計(jì)算效率的預(yù)測(cè)銑削穩(wěn)定性的方法。

        4 實(shí)驗(yàn)驗(yàn)證

        針對(duì)提出的用于預(yù)測(cè)銑削穩(wěn)定性的HAMM算法,本節(jié)通過實(shí)驗(yàn)來進(jìn)行驗(yàn)證。本次驗(yàn)證實(shí)驗(yàn)所用機(jī)床為五軸立銑加工中心DMU60monoBLOCK,所用銑刀為整體式硬質(zhì)合金4刃刀具,直徑為12 mm,刀長為75 mm。被加工工件材料為航空鋁合金AL7075,實(shí)驗(yàn)采用無冷卻液的干式切削。試驗(yàn)現(xiàn)場(chǎng),如圖6所示。

        圖6 試驗(yàn)現(xiàn)場(chǎng)Fig.6 Testing site

        在切削鋁合金工件的過程中,切削參數(shù)為:徑向比a/D=1,切向力系數(shù)Kt=6×108N/m2,法向力系數(shù)Kn=2×108N/m2,主軸轉(zhuǎn)速為3 000 r/min,3 500 r/min,4 000 r/min,4 500 r/min和5 000 r/min,軸向切削深度為0.4 mm,0.8 mm,1.2 mm,1.6 mm,2.0 mm。模態(tài)參數(shù)由模態(tài)錘擊試驗(yàn)獲得:逆銑,刀具齒數(shù)N=4,模態(tài)質(zhì)量mt=0.1 kg,固有圓頻率ωn=9 420 rad/s,固有阻尼比ζ-0.029 81。

        選用不同的主軸轉(zhuǎn)速和軸向切削深度進(jìn)行切削,利用激光測(cè)振儀搜集振動(dòng)信號(hào),并利用統(tǒng)計(jì)方法進(jìn)行分析得穩(wěn)定點(diǎn)和不穩(wěn)定點(diǎn)。基于HAMM法預(yù)測(cè)銑削穩(wěn)定性的有效性,結(jié)合實(shí)驗(yàn)參數(shù),計(jì)算機(jī)仿真結(jié)果和試驗(yàn)結(jié)果,如圖7所示。

        圖7 HAMM法的計(jì)算機(jī)仿真結(jié)果和實(shí)驗(yàn)結(jié)果Fig.7 Computer simulation and experimental results of HAMM method

        圖8 穩(wěn)定點(diǎn)B的時(shí)域圖Fig.8 Time-domain diagram of stability point B

        穩(wěn)定點(diǎn)B的時(shí)域圖如圖8。

        圖9 不穩(wěn)定點(diǎn)A的時(shí)域圖Fig.9 Time-domain diagram of unstability point A

        不穩(wěn)定點(diǎn)A的時(shí)域圖如9。

        由圖7可以看出,HAMM方法的計(jì)算機(jī)仿真結(jié)果和實(shí)驗(yàn)結(jié)果是吻合的,進(jìn)一步說明HAMM方法是一種有效的預(yù)測(cè)銑削穩(wěn)定性的方法。

        5 結(jié) 論

        針對(duì)銑削動(dòng)力學(xué)模型,提出了HAMM預(yù)測(cè)銑削穩(wěn)定性。通過以上研究,得到如下一些結(jié)論:

        通過對(duì)局部離散誤差的分析可知,HAMM的收斂速度比1st-SDM和2nd-FDM快得多。主軸轉(zhuǎn)速越高,離散數(shù)少的HAMM能達(dá)到離散數(shù)較多的1st-SDM和2nd-FDM的局部離散誤差,和理論分析結(jié)論一致。

        無論是在單自由度還是雙自由度銑削模型中,從所給的三種方法的穩(wěn)定性葉瓣圖中可以得出,隨著離散數(shù)的增加,三種方法的計(jì)算精度都有所提高。當(dāng)離散數(shù)相同時(shí),HAMM的計(jì)算精度高于1st-SDM和2nd-FDM。特別是在雙自由度銑削模型下,HAMM的計(jì)算精度遠(yuǎn)遠(yuǎn)高于1st-SDM和2nd-FDM。在離散數(shù)相同的條件下,HAMM的計(jì)算效率和2nd-FDM相當(dāng),遠(yuǎn)遠(yuǎn)高于1st-SDM。

        HAMM是一種較優(yōu)的預(yù)測(cè)銑削穩(wěn)定性的方法,它具有較高的預(yù)測(cè)精度和較快的計(jì)算效率。實(shí)驗(yàn)驗(yàn)證表明,HAMM方法用來預(yù)測(cè)銑削穩(wěn)定性是有效的。

        DOI:10. 1007/s00170-016-8708-z.

        猜你喜歡
        模態(tài)方法模型
        一半模型
        重要模型『一線三等角』
        重尾非線性自回歸模型自加權(quán)M-估計(jì)的漸近分布
        可能是方法不對(duì)
        3D打印中的模型分割與打包
        用對(duì)方法才能瘦
        Coco薇(2016年2期)2016-03-22 02:42:52
        國內(nèi)多模態(tài)教學(xué)研究回顧與展望
        四大方法 教你不再“坐以待病”!
        Coco薇(2015年1期)2015-08-13 02:47:34
        捕魚
        基于HHT和Prony算法的電力系統(tǒng)低頻振蕩模態(tài)識(shí)別
        国产对白国语对白| 加勒比一本大道大香蕉| 国产av普通话对白国语| 丝袜人妻中文字幕首页| 最新国产女主播在线观看| 亚洲av综合色区无码一区| 日本阿v片在线播放免费| 午夜精品一区二区三区的区别 | 老司机在线免费视频亚洲| 日韩中文字幕素人水野一区| 亚洲国产精品成人天堂| 精品久久久久香蕉网| 亚洲av无码av制服丝袜在线 | 日本中文字幕一区二区高清在线 | 日韩在线一区二区三区免费视频| 少妇的丰满3中文字幕| 精品视频专区| 亚洲av区一区二区三区| 国产一区二区三区天堂| 一女被多男玩喷潮视频| 国产在线一区二区三区av| 波多野结衣视频网址| 亚洲色AV天天天天天天 | 亚洲中文字幕一区av| 蜜桃18禁成人午夜免费网站| 亚洲欧美中文字幕5发布| 国产精品久久久久久久免费看| 青青视频一区| 中文字幕乱码亚洲美女精品一区 | 中文字字幕在线精品乱码| 亚洲AV日韩AV永久无码电影| 最新欧美一级视频| 亚洲熟少妇一区二区三区| 包皮上有一点一点白色的| 看久久久久久a级毛片| 亚洲av无码片vr一区二区三区| 国产精成人品| 免费一级黄色大片久久久| 免费国产自拍视频在线观看| 中文字幕漂亮人妻在线| 成人免费直播|