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

        ?

        基于接觸區(qū)域的球頭刀五軸加工無顫振刀具姿態(tài)研究

        2023-10-12 07:14:38植俊杰張立強(qiáng)許洋洋楊青平曹珍珍
        關(guān)鍵詞:球頭切削力姿態(tài)

        植俊杰,張立強(qiáng)+,許洋洋,楊青平,曹珍珍

        (1.上海工程技術(shù)大學(xué) 機(jī)械與汽車工程學(xué)院,上海 201620;2.成都永峰科技有限公司,四川 成都 610511)

        0 引言

        五軸加工技術(shù)用于生產(chǎn)航空航天、模具、汽車和生物醫(yī)學(xué)工業(yè)中的復(fù)雜零件,不同于三軸加工,五軸加工中刀具姿態(tài)的引入對加工力學(xué)和動力學(xué)有顯著的影響,使得預(yù)測更為復(fù)雜。

        刀具姿態(tài)通常由前傾角和側(cè)傾角描述,OZTURK等[1]指出切削刀具與工件接觸的區(qū)域取決于前傾角與側(cè)傾角、刀具的幾何形狀以及切削深度,準(zhǔn)確高效地獲取加工中刀具工件的接觸區(qū)域是進(jìn)行切削力預(yù)測和穩(wěn)定性分析的前提,具有非常重要的意義。接觸區(qū)域的獲取方法總體劃分為實體建模法、離散法和解析法。ARAS等[2]通過實體建模進(jìn)行布爾運算獲得刀具工件的接觸區(qū)域。MA等[3]基于VERICUT進(jìn)行二次開發(fā),實現(xiàn)了在仿真過程中提取刀具工件的接觸區(qū)域。BOZ等[4]運用實體建模法獲取刀具工件接觸區(qū)域并基于此進(jìn)行了切削力分析。KISWANTO等[5]提出一種基于混合分析和離散的方法來確定五軸加工中刀具工件的接觸區(qū)域。代月幫等[6]運用半解析法搭建了球頭銑刀與工件接觸區(qū)域的邊界投影方程。GHORBANI等[7]提出一種新的刀具工件接觸區(qū)域的解析模型,采用表面主曲率的二次數(shù)學(xué)表示來描述瞬時切削區(qū)域的周圍的工件表面。離散法獲取接觸區(qū)域的效率受限于網(wǎng)格離散的精度,當(dāng)離散精度高時,將會降低計算效率,難以平衡兩者。實體建模法雖然精度高,但是刀位點或刀具姿態(tài)發(fā)生變化時,都需要重新進(jìn)行布爾運算,計算效率低。相對于離散法和實體建模法,解析法采用數(shù)學(xué)模型表示刀具工件接觸區(qū)域的邊界,在保持一定的精度時也具有較高的計算效率。然而,由于沿著刀具軌跡工件表面的曲率連續(xù)變化,解析提取球頭銑刀在五軸銑削加工中的瞬時接觸區(qū)域邊界仍然是一個挑戰(zhàn)。

        此外,對于銑削加工中的穩(wěn)定性,李堯等[8]提出一種基于小波包及Hilbert-Huang變換的數(shù)控銑削顫振診斷技術(shù)。韓振宇等[9]提出一種基于ESPRIT頻譜估計和隱馬爾可夫模型辨識銑削過程中穩(wěn)定性的方法,實現(xiàn)了切削顫振的有效識別。穩(wěn)定性預(yù)測是獲得無顫振工藝參數(shù)的基礎(chǔ),ALTINTAS等[10]提出一種繪制穩(wěn)定性葉瓣圖的解析預(yù)測方法,用以獲得穩(wěn)定的加工參數(shù)。MINIS等[11]利用Nyquist準(zhǔn)則求得了銑削過程的穩(wěn)定域解。隨后,時域法[12]、有限元法[13]和半離散法[14]也被提出用于預(yù)測銑削加工穩(wěn)定性,DING等[15]在對半離散法總結(jié)的基礎(chǔ)上通過離散實際的時間—周期項提出一種全離散法,提高了求解效率。以上方法大部分是通過繪制穩(wěn)定性葉瓣圖,并基于葉瓣圖改變切削深度和主軸轉(zhuǎn)速來實現(xiàn)顫振抑制。若穩(wěn)定性葉瓣圖中某可行轉(zhuǎn)速區(qū)域非常狹窄,則很難通過調(diào)整主軸轉(zhuǎn)速達(dá)到穩(wěn)定性加工要求。SHAMOTO等[16]的研究表明,刀具姿態(tài)的改變會影響多軸加工的穩(wěn)定性,可通過改變刀具姿態(tài)避免多軸銑削加工過程中的振動。隨后,ALAN等[17]提出一種預(yù)測球頭銑刀穩(wěn)定性的迭代方法,該方法利用圓盤單元來確定精確的切入角和切出角,提高了銑削力的計算精度。SUN等[18]計算了加工系統(tǒng)的特征方程,通過搜索算法獲得了每個刀位點處無顫振的刀具姿態(tài)。WANG等[19]提出了構(gòu)造姿態(tài)穩(wěn)定性圖的方法用于顫振預(yù)測,并通過實驗驗證了姿態(tài)穩(wěn)定性圖的有效性,因此如何準(zhǔn)確地構(gòu)造出針對球頭銑刀的姿態(tài)穩(wěn)定性圖意義重大。

        綜上所述,為解決球頭刀銑削過程接觸區(qū)域獲取困難和難以調(diào)整主軸轉(zhuǎn)速達(dá)到穩(wěn)定性銑削的問題,本文建立了考慮刀具姿態(tài)與螺旋角影響的動態(tài)切削力模型,通過NX二次開發(fā)提取刀具與工件接觸區(qū)域的初始數(shù)據(jù),并采用一種坐標(biāo)變換法獲得不同刀具姿態(tài)下的接觸區(qū)域。然后,使用全離散法判斷在不同刀具姿態(tài)下的銑削穩(wěn)定性,構(gòu)建了姿態(tài)穩(wěn)定性圖。最后,通過實驗驗證了該方法的有效性。

        1 球頭銑刀五軸加工動態(tài)切削力模型

        切削力作為加工系統(tǒng)的激勵輸入,是進(jìn)行顫振預(yù)測的基礎(chǔ)。在五軸銑削加工中,球頭銑刀通常用于半精加工或精加工,僅球頭刀端部球體部分參與銑削。刀軸矢量在一定范圍內(nèi)變化,刀具工件接觸體并不會發(fā)生變化,但沿著刀具路徑,每個刀位處的刀軸矢量、進(jìn)給方向和工件間的關(guān)系在不斷變化。

        如圖1a所示,為參數(shù)化表達(dá)刀軸矢量、進(jìn)給方向和工件間的關(guān)系,分別定義了工件坐標(biāo)系、進(jìn)給坐標(biāo)系和刀具坐標(biāo)系。工件坐標(biāo)系為OW-XWYWZW,是CAM加工軟件中最基本的坐標(biāo)系之一,在此坐標(biāo)系下生成刀位點和刀軸矢量。刀具坐標(biāo)系為O-XTYTZT,其ZT沿著刀軸方向。進(jìn)給坐標(biāo)系為O-FCN,F表示進(jìn)給方向,C表示橫向進(jìn)給方向,N為加工表面法向,其通過CAM軟件生成的刀位文件信息確定。根據(jù)圖1b可知,刀具坐標(biāo)系通過進(jìn)給坐標(biāo)系先繞C軸旋轉(zhuǎn)前傾角l,再繞F軸旋轉(zhuǎn)側(cè)傾角t確定,兩者間的旋轉(zhuǎn)矩陣為:

        圖1 五軸加工球頭銑刀幾何坐標(biāo)系

        (1)

        球頭銑刀銑削幾何模型如圖2所示[20],在刀具坐標(biāo)系中進(jìn)行切削力建模。其中R(z)為局部半徑,κ(z)為軸向浸入角,ψ(z)為徑向滯后角,φj(z)為浸入角。由于球頭銑刀球頭部分刀刃為螺旋形,因此將刀沿刀軸方向離散為n層,并假定層上的切削微元是筆直且處于傾斜切削狀態(tài)。

        圖2 球頭銑刀銑削幾何模型

        根據(jù)LEE等[21]提出的切削力模型,單個微元切削產(chǎn)生的切削力可以被分解為切向力dFt,徑向力dFr和軸向力dFa,表示為:

        dFt=KtedS+Ktchuct(φj(z),k(z))db;dFr=KredS+Krchuct(φj(z),k(z))db;

        dFa=KaedS+Kachuct(φj(z),k(z))db。

        (2)

        其中:Kte,Kre,Kae為犁切力系數(shù);Ktc,Krc,Kac為剪切力系數(shù);dS為微元切削刃的長度;db為切削寬度;huct(φj(z),k(z))為瞬時未變形切削厚度,可計算如下:

        huct(φj(z),k(z))=u·d。

        (3)

        其中:

        (4)

        (5)

        (6)

        式中:Ω為主軸轉(zhuǎn)速;β為球頭銑刀螺旋角;N為球頭銑刀的刀刃數(shù)。每個切削微元切削力可以通過以下等式轉(zhuǎn)換到刀具坐標(biāo)系中:

        (7)

        其中:

        (8)

        (9)

        因此,將每個切削微元切削力進(jìn)行累積求和即可得出作用在刀具上的總切削力。

        2 接觸區(qū)域提取與轉(zhuǎn)換

        由式(8)可知,在加工過程中,切削刃只在刀具工件接觸區(qū)域參與切削。因此,獲取球頭刀五軸銑削加工中刀具工件接觸區(qū)域?qū)η邢髁Φ慕O為重要。為了提高刀具-工件接觸區(qū)域的獲取效率,本文提出一種改進(jìn)的實體建模法。首先,基于NX 12.0二次開發(fā)提取仿真加工中初始刀具工件接觸區(qū)域。然后,通過坐標(biāo)變換計算獲取當(dāng)前刀位處所需刀具姿態(tài)下的刀具工件接觸區(qū)域。整體上分為接觸區(qū)域的提取和接觸區(qū)域的轉(zhuǎn)換兩大步驟,其總體流程如圖3所示。

        圖3 接觸區(qū)域獲取流程

        2.1 接觸區(qū)域提取

        刀具工件接觸區(qū)域是基于實體相交法提取沿刀路加工中每一刀位的刀具工件相交體,該方法的基本思想是從原材料模型中去除刀具沿數(shù)控軌跡運動產(chǎn)生的掃掠體積,從而獲得加工至此刀位的工件體。如圖4所示,使用NX軟件獲得當(dāng)前刀位的刀具與工件相交體,便可以提取此刀位刀具工件的接觸區(qū)域。根據(jù)圖4b,球頭刀沿刀軸方向被逐層離散化,通過將圓盤與接觸區(qū)域的邊界求交,獲得點集的坐標(biāo)。由于此時點集均是在工件坐標(biāo)系下獲取,可通過坐標(biāo)轉(zhuǎn)換矩陣TW-T將點集轉(zhuǎn)換到刀具坐標(biāo)系中。TW-T表示在此刀位工件坐標(biāo)系到刀具坐標(biāo)系的轉(zhuǎn)換矩陣,其計算如下:

        圖4 基于NX12.0二次開發(fā)獲取接觸區(qū)域

        (10)

        其中:

        (11)

        (12)

        式中:RW-T為工件坐標(biāo)系到刀具坐標(biāo)系的旋轉(zhuǎn)矩陣,CLW為在工件坐標(biāo)系中刀具球心處的坐標(biāo)值,RF-W為進(jìn)給坐標(biāo)系到工件坐標(biāo)系的旋轉(zhuǎn)矩陣。

        (13)

        為求得沿刀軸方向上每一切削層所對應(yīng)的切入角和切出角,在計算出每一點對應(yīng)的浸入角后,需首先相對于z坐標(biāo)將浸入角進(jìn)行排序。然后,將具有相同z坐標(biāo)的浸入角分為一組。最后,將每一組相對于浸入角按升序排序,并成對依次從各組中獲得切入角和切出角。如圖5所示為使用MATLAB繪制的D12R6球頭銑刀在圖4所示位置處的ztφb圖,其刀具球端離散層數(shù)為60層,切削深度為1 mm,前傾角為10°,側(cè)傾角為15°,橫坐標(biāo)表示浸入角,縱坐標(biāo)表示軸向切深,藍(lán)色點表示切出角,紅色點表示切出角。

        圖5 l=10°,t=15°的ztφb圖

        2.2 接觸區(qū)域轉(zhuǎn)換

        由于在加工中的某一刀位點處,僅球頭刀端部球體部分參與銑削,刀具姿態(tài)在一定范圍的變化并不會使刀具工件接觸體發(fā)生變化,只會影響該接觸區(qū)域在球頭表面的分布,如圖6所示。因此,可以通過初始刀具姿態(tài)中的接觸區(qū)域信息重新構(gòu)建任意刀具姿態(tài)下的刀具工件接觸區(qū)域。

        圖6 不同刀具姿態(tài)接觸區(qū)域分布

        圖7 初始刀具坐標(biāo)系與新刀具坐標(biāo)系

        (14)

        通過相應(yīng)的矩陣逆運算,可求得新的刀具姿態(tài)對應(yīng)的刀軸矢量ta在初始刀具坐標(biāo)系下的表達(dá)式。設(shè)ta=(i,j,k),則α角和β角可通過以下等式計算:

        (15)

        (16)

        因此,在獲得初始接觸區(qū)域信息后即可通過轉(zhuǎn)換矩陣TT-T′求得所需刀具姿態(tài)下的接觸區(qū)域信息。

        3 五軸球頭銑刀動力學(xué)建模

        如圖8所示,在再生顫振機(jī)理的模型中,假定刀具具有兩個自由度,兩個自由度相互正交并且均與刀具的軸線正交,因此可以通過兩個獨立的模態(tài)錘擊實驗獲得沿這兩個方向的模態(tài)參數(shù)。模態(tài)坐標(biāo)系OM-XMYMZM附著在刀具上,其中XM和YM與兩個錘擊方向一致,ZM與刀具軸向方向一致。

        圖8 再生顫振機(jī)理圖

        顫振主要是自激振動,受周期性變化的銑削力影響。對于式(2),在加工過程中,不受瞬時未變形切削厚度huct(φj(z),k(z))影響的靜態(tài)部分的切削力不會對顫振的產(chǎn)生造成影響,因此可以舍去。結(jié)合圖2d可得db=dz/sink(z),則式(7)可以重寫如下:

        (17)

        將每個切削微元切削力進(jìn)行累積求和即可得出作用在刀具上的總動態(tài)切削力:

        (18)

        (19)

        式(19)可以重寫為:

        (20)

        到目前為止,刀具位移和切削力均是在刀具坐標(biāo)系中計算。但是,由于在機(jī)床上進(jìn)行錘擊實驗時,刀具坐標(biāo)系會隨著刀具位置的改變而發(fā)生變化,無法將模態(tài)坐標(biāo)系與每個刀具坐標(biāo)系對齊,需構(gòu)建刀具坐標(biāo)系到模態(tài)坐標(biāo)系的轉(zhuǎn)換矩陣。進(jìn)行實驗加工所使用的機(jī)床為C-100P臥式五軸加工中心,其主運動由X、Y、Z軸的3個直線運動和A、B軸兩個旋轉(zhuǎn)運動組成,如圖9所示。

        圖9 機(jī)床結(jié)構(gòu)示意圖

        為方便起見,當(dāng)每個旋轉(zhuǎn)軸位于機(jī)床上的零初始位置時,將錘擊方向選擇為沿著機(jī)床工作臺的X方向和Y方向,當(dāng)銑削加工采用某種刀具姿態(tài)時,模態(tài)坐標(biāo)也隨著變動。因此,模態(tài)坐標(biāo)系與刀具坐標(biāo)系之間的轉(zhuǎn)換矩陣將由機(jī)床的結(jié)構(gòu)確定。模態(tài)坐標(biāo)系到進(jìn)給坐標(biāo)系的轉(zhuǎn)換可由下式求得:

        TM-F=(TF-W)-1·TTa-W·Tz(θB)·Tx(90°)·
        Ty2(90°)·Tz(θA)·Ty1(90°)。

        (21)

        其中TTa-W為從機(jī)床工作臺到工件的轉(zhuǎn)換矩陣,由在工作臺上設(shè)置工件的方式?jīng)Q定。本文將工件坐標(biāo)系和機(jī)床工作臺坐標(biāo)系設(shè)為對齊,此時式(21)中其他矩陣計算如下:

        (22)

        (23)

        (24)

        (25)

        (26)

        根據(jù)式(1)和式(21),可以通過以下等式得到從刀具坐標(biāo)系到模態(tài)坐標(biāo)系的變換矩陣:

        (27)

        然后,通過在式(20)的兩邊同時乘TT-M,便可將切削力轉(zhuǎn)換至模態(tài)坐標(biāo)系表示為:

        (28)

        由于z方向即軸向的剛度很高,在建立五軸銑削加工動力學(xué)方程時,通常會舍棄z方向的項,因此考慮再生顫振的五軸銑削系統(tǒng)的動力學(xué)方程表達(dá)為:

        (29)

        式(29)可改寫為:

        (30)

        其中M,C,K分別表示通過錘擊實驗得到的模態(tài)質(zhì)量、模態(tài)阻尼和模態(tài)剛度。

        式(30)是一個時滯微分方程,多年來,學(xué)者們采用了多種方法對時滯微分方程進(jìn)行求解,根據(jù)求解結(jié)果來判斷實際加工過程中的穩(wěn)定性,如時域仿真、零階近似、單頻率和多頻率法以及半離散法和有限元法等不斷地被提出,用以進(jìn)行加工穩(wěn)定性判斷。為了達(dá)到更高的計算效率,DING等[15]在半離散法的基礎(chǔ)上提出了全離散法,隨后DAI等[22]、DONG等[23]和ZHI等[24]分別基于黃金搜索法、更新的數(shù)值積分法和隱式Adams法對全離散法進(jìn)行了改進(jìn),進(jìn)一步提升了全離散法的計算效率。全離散法是目前穩(wěn)定性預(yù)測最準(zhǔn)確、最先進(jìn)的方法,因此本文將使用一種高效的全離散法[23]對球頭刀五軸銑削加工穩(wěn)定性進(jìn)行判斷。

        4 實驗驗證

        如圖10所示,進(jìn)行切削力系數(shù)辨識實驗[25]的機(jī)床為立臥轉(zhuǎn)換式五軸加工中心,型號是C-100P;所用的刀具為2齒硬質(zhì)合金球頭銑刀,刀具直徑為12 mm,螺旋角為30°,工件材料為AL7050-T7451。

        圖10 切削力系數(shù)標(biāo)定

        圖10中的測力硬件為Kistler測力儀、電荷放大器和信號收集器,通過使用Manuware和Hrsoft_DW_V2.03兩款商用軟件進(jìn)行切削力系數(shù)標(biāo)定實驗,測得的切削力系數(shù)具體結(jié)果如表1所示。

        表1 切削力系數(shù)

        在模態(tài)實驗中,刀具系統(tǒng)X方向和Y方向上的模態(tài)參數(shù)通過錘擊實驗進(jìn)行獲取。根據(jù)實驗情況,由于力錘激發(fā)的頻率范圍與撞擊的接觸時間成反比,且敲擊物體為硬質(zhì)尖端和短撞擊接觸時間,因此所選用的力錘型號為Kistler 9722A2000。此外,在模態(tài)實驗中,測量的位移非常小,使用位移傳感器很難在較寬的頻率范圍內(nèi)測量振動,加速度傳感器在振動測量中更為常見,因此所選用的加速度傳感器型號為PCB 353B11 S/N 9690。

        使用加速度傳感器的測量設(shè)置如圖11所示,其通過使用特殊的石蠟安裝在刀具上,使用商用測量軟件Cutpro的AMTIF模塊進(jìn)行信號采集,圖中虛線顯示了敲擊的方向,在這種情況下,將測得Y方向上的傳遞函數(shù),單一方向上連續(xù)5次采集成功后,便可使用MODAL模塊進(jìn)行計算獲得該方向的模態(tài)參數(shù)。如圖12所示為根據(jù)5次測量結(jié)果得到的擬合曲線,根據(jù)擬合曲線與實際測量的吻合度可判斷此次為有效測量。

        圖11 刀具系統(tǒng)模態(tài)實驗

        圖12 測量擬合效果圖

        同理,可以改變方向進(jìn)行錘擊實驗獲得X方向上的模態(tài)參數(shù),多次實驗取平均值后的最終計算結(jié)果如表2所示。

        表2 模態(tài)參數(shù)

        在獲取切削力系數(shù)和模態(tài)參數(shù)后,便可使用高效的全離散法構(gòu)建每一刀位處姿態(tài)穩(wěn)定性圖,根據(jù)姿態(tài)穩(wěn)定性圖和刀軸光順輸出優(yōu)化后的刀具路徑。本文構(gòu)建了在固定主軸轉(zhuǎn)速為4 800 r/min,切削深度為1 mm情況下的平面槽銑加工的姿態(tài)穩(wěn)定性圖,如圖13所示。圖中黑色輪廓線λ=1為穩(wěn)定區(qū)域與顫振區(qū)域之間的邊界;當(dāng)λ>1時,加工會產(chǎn)生振紋;當(dāng)λ<1時,加工處于穩(wěn)定性狀態(tài)。由于在平面槽銑加工中,每一刀位處的刀具工件接觸體不變,可根據(jù)圖13得到滿足全局穩(wěn)定加工的刀具姿態(tài)。

        圖13 姿態(tài)穩(wěn)定性圖

        為了驗證生成的姿態(tài)穩(wěn)定性圖的有效性,根據(jù)姿態(tài)穩(wěn)定性圖選取了10組不同的刀具姿態(tài)進(jìn)行了槽銑加工實驗,排除了加工中的幾何干涉和機(jī)床旋轉(zhuǎn)軸限制的影響,具體結(jié)果如表3所示。由表3可知,序號1,5,7,8,9對應(yīng)的刀具姿態(tài)下,其工件加工表面出現(xiàn)明顯的顫振痕跡;序號2,3,4,6,10對應(yīng)的刀具姿態(tài)下,工件處于穩(wěn)定性加工,工件表面無振刀紋。由預(yù)測結(jié)果與實驗結(jié)果之間的比較可知,預(yù)測結(jié)果良好,預(yù)測與實驗結(jié)果之間已經(jīng)達(dá)到了較高的一致性。

        表3 預(yù)測結(jié)果與實驗結(jié)果的對比

        5 結(jié)束語

        本文利用分層離散化思想并結(jié)合NX 12.0的二次開發(fā)獲取球頭銑刀與工件的接觸區(qū)域,求得離散層與刀具工件接觸區(qū)域邊界的交點,根據(jù)幾何關(guān)系方程求出各離散層對應(yīng)的切入角和切出角,用于確定瞬時參與切削的刀刃的實際切削部位。建立了考慮刀具姿態(tài)與螺旋角影響的動態(tài)切削力模型,將接觸區(qū)域的獲取運用到模型中,采用高效的全離散法構(gòu)造出了姿態(tài)穩(wěn)定性圖,可用于解決傳統(tǒng)的穩(wěn)定性葉瓣圖可行轉(zhuǎn)速區(qū)域狹窄的問題。在難以調(diào)整主軸轉(zhuǎn)速情況下,通過調(diào)整刀具姿態(tài)以實現(xiàn)穩(wěn)定性加工。后續(xù)工作中,將開展考慮五軸加工特殊性的刀具切削力系數(shù)標(biāo)定研究,進(jìn)一步提高預(yù)測的準(zhǔn)確性。此外,在穩(wěn)定性約束下,結(jié)合幾何約束、機(jī)床旋轉(zhuǎn)軸約束以及刀軸光順約束,輸出可用于加工復(fù)雜幾何零件的多約束優(yōu)化刀路也值得研究。

        猜你喜歡
        球頭切削力姿態(tài)
        攀爬的姿態(tài)
        帶滾珠調(diào)節(jié)的面—面接觸形式負(fù)荷傳感器球頭副
        30CrMnSiNi2A插銑加工的切削力研究
        錐面球頭銷設(shè)計
        全新一代宋的新姿態(tài)
        汽車觀察(2018年9期)2018-10-23 05:46:40
        圓刀片切削力計算方法
        跑與走的姿態(tài)
        中國自行車(2018年8期)2018-09-26 06:53:44
        城軌用帶串聯(lián)間隙球頭避雷器壽命分析
        電氣化鐵道(2016年2期)2016-05-17 03:42:34
        QP-16型球頭掛環(huán)疲勞壽命試驗研究與分析*
        一種深孔鉆削切削力測量裝置
        真实国产乱啪福利露脸| 91精品国产综合久久久蜜| 吃奶摸下高潮60分钟免费视频| 国产熟女露脸大叫高潮| 在线观看av手机网址| 国产99视频一区二区三区| 一区二区三区人妻少妇| 99久久久无码国产精品6| 国产精品美女白浆喷水| 一本久久伊人热热精品中文| 精品国产一区二区三区三级| 少妇饥渴偷公乱a级无码| 中文字幕一区二区人妻出轨| 日韩精品资源在线观看免费| 午夜影视免费| 国产人妻久久精品二区三区| 青青青伊人色综合久久亚洲综合| 青青草视频在线播放观看| 日韩av午夜在线观看| 亚洲成人小说| 中文字幕一区二区网站| 五月婷婷六月丁香久久综合| 久久99精品九九九久久婷婷| 未满十八勿入av网免费| 国产视频免费一区二区| 狠狠色丁香婷婷久久综合| 久久人人爽人人爽人人av东京热| 亚洲欧洲无码精品ⅤA| 国产在线91精品观看| 国产av麻豆mag剧集| 国产激情久久99久久| 美女视频在线观看一区二区三区| 亚洲熟妇丰满多毛xxxx| 亚洲永久无码动态图| 亚洲一区二区三区国产精品视频| 日韩 亚洲 制服 欧美 综合 | 中文字幕久久精品一二三区| 亚洲AV无码精品一区二区三区l| 亚洲国产女同在线观看| 国内免费高清在线观看| 国产福利午夜波多野结衣|