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

基于熵權法的含間隙和柔性的機構定量分析方法及應用

2023-10-18 03:47:50潘佳煊錢孟波孫福興
振動與沖擊 2023年19期
關鍵詞:方向影響模型

潘佳煊, 錢孟波, 孫福興, 虞 浪, 陳 強

(浙江農林大學 光機電工程學院,杭州 311300)

并聯機構具有整體結構穩定、精度高和承載能力大的優點,廣泛應用于多個領域。將2-RR&2-PR并聯機構作為農業中谷物清選振動篩的主要機構,相比于傳統振動篩,使篩面獲得空間多維運動,加快物料的散開速度,提高透篩效率,增加使用壽命。

在以往的振動篩系統運動特性研究中,學者們都將整個振動篩系統簡化為一個多剛體系統[1-2]。這種方法忽視了系統中構件的材料彈性屬性,在實際工作過程中相應的構件會產生彈性變形,并且由于機械加工制造和機器在高速情況下工作而造成的摩擦磨損等因素產生間隙,以上都可能會影響機構運動時的穩定性和使用壽命。因此,有必要綜合考慮間隙與構件柔性對機構動態特性的影響,并量化對機構各指標的影響程度,從而提供理論參考。

近年來,學者們在對機構運動特性分析中考慮并引入了間隙或構件柔性。Flores[3]建立了連續接觸碰撞力模型,并對具有間隙的平面四桿機構進行了分析。Bai等[4-5]提出了一種改進的非線性連續接觸力混合模型,以四連桿機構為對象分析了轉動關節間隙的平面機械系統的運動特性。Qian[6]等用帶間隙的三維平移關節進行了實驗,以揭示滑道和滑塊之間的相對運動特性,利用Poincaré映射和Lyapunov指數驗證了帶間隙平移關節中存在混沌行為。Erkaya等[7]提出含間隙四連桿機構的動力學模型,為了控制連桿的柔性改變連桿的截面尺寸,來研究柔性對間隙機構動態特性的影響。孫杰等[8]提出一種含間隙鉸接的航天器剛柔耦合模型建立與控制方法,來減緩間隙對航天器動態特性造成的影響。

信息熵作為一種非常重要的非線性分析方法,被廣泛運用于計算機、機械、數學、金融等領域。基于信息熵提出了很多相關的方法,如交叉熵[9]、功率譜熵[10]、基于EEMD的奇異譜熵[11]等。Qian等[12]運用信息熵理論,提出了一種含間隙機構動力響應非線性分析的新方法,對圓柱間隙關節機構的非線性行為進行定量分析。熵權法是一種基于信息熵理論的目標值分配方法,關于熵權法的研究,林巨廣等[13]采用熵權法對軸承故障進行診斷,取得了較好的故障識別效果。蔣榮超等[14]提出基于熵權法的結構綜合貢獻系數計算方法,以此為評價指標篩選出對懸架性能影響較大的結構參數作為整車性能匹配優化的設計變量。

現有基于信息熵理論對含間隙機構運動特性非線性定量分析研究國內外鮮有人做,本文基于熵權法對含間隙和構件柔性的并聯機構進行非線性定量分析,以2-RR&2-PR并聯機構為研究對象,首先利用ADAMS二次開發功能,通過Fortran語言自主編寫間隙接觸力子程序,加載到ADAMS函數求解庫中,建立了考慮間隙的機構動力學模型,然后開展考慮間隙和構件柔性的機構運動特性仿真,通過仿真數據對比分析在不同條件下對機構運動特性的影響,最后采用Lyapunov指數驗證其運動非線性,在此基礎上提出一種基于熵權法的非線性定量分析新方法,評估間隙和構件柔性對機構非線性的影響程度,對機構進行更準確地描述。

1 含間隙2-RR&2-PR并聯機構動力學模型與分析模型

1.1 并聯機構描述

2-RR&2-PR并聯機構的機構示意圖如圖1所示,該機構由動平臺、靜平臺和兩條支鏈組成。其中兩條支鏈上,A1、A2、A3、A4、A5、A6為6個轉動副,B1、B2為2個移動副。整個機構可視為對稱布置、同步驅動的兩曲柄滑塊機構A1A3A5B1與A2A4A6B2,因其連桿A3A5、A4A6運動規律一致,故可將A3A4A6A5視為同一剛體,即動平臺。通過2個移動副的驅動,實現動平臺的上下振動,從而實現整個機構的運動。機構的模型參數如表1所示,本節先考慮轉動副A1處的間隙,并將構件A3A5和A4A6柔性化處理。

圖1 含間隙和柔性桿的2-RR&2-PR并聯機構示意圖Fig.1 Schematic diagram of 2-RR&2-PR parallel mechanism with clearance and flexible rod

表1 機構構件參數Tab.1 Mechanism component parameter

1.2 含間隙旋轉副的接觸力模型

機構在運動過程中,元素之間會相對運動,在運動副處往往會因為磨損產生間隙,間隙的存在會導致元素間的接觸和碰撞。本文采用L-N接觸力模型和修正的Coulomb摩擦力模型,來建立含間隙旋轉副的接觸力模型。

1.2.1 間隙矢量模型

建立合理的間隙矢量模型,用來描述鉸接處兩元素相對位置的變化,如圖2所示。

圖2 間隙矢量模型Fig.2 Clearance vector model

(1)

當兩元素發生碰撞時,碰撞深度為

δ=e-c

(2)

式中:e為兩元素的偏心距;c為兩元素的半徑差。將δ作為判斷軸套和軸銷是否發生碰撞的標準,當δ<0時,兩元素未接觸,自由運動;當δ=0時,兩元素開始接觸或開始脫離;當δ>0時,兩元素發生碰撞。

1.2.2 法向接觸力模型

法向接觸力模型表示運動副處元素間接觸時的法向作用力。目前應用最廣泛的接觸力模型是由Lankarani和Nikravesh兩位科學家在Hertz接觸力數學模型基礎之上提出來的L-N模型,該模型考慮了由材料阻尼引起的能量損失[15]。

基于L-N模型的法向接觸力Fn表達式如下

(3)

接觸剛度系數K的大小取決于接觸表面的幾何形狀,對于兩圓柱面的接觸,K表達式為

(4)

(5)

式中:σi、σj為接觸體的材料系數;Ri、Rj為接觸體的接觸半徑;E為材料彈性模量;v為材料泊松比。

阻尼系數D的表達式為

(6)

將式(4)和式(6)代入式(3),即可得到法向接觸力Fn表達式

(7)

1.2.3 切向摩擦力模型

目前Coulomb摩擦力模型是應用最廣泛的摩擦力模型,適用于碰撞中的摩擦行為。然而當摩擦力由靜摩擦向動摩擦過渡時,該模型與實際情況不符,本文采用修正的Coulomb摩擦力模型,其表達式如下[17-18]

(8)

(9)

式中:vt、vs和vd分別為相對滑動速度,靜摩擦速度和動摩擦速度;μs為靜摩擦因數;μd為動摩擦因數。

1.3 熵權法

熵權法就是依據各指標數據所包含的信息量大小來確定指標的權重。設影響機構運動非線性的因素的向量集合為{E1,E2,…,En},每個集合存在m個判定單元,那么關于機構運動非線性影響因素的判斷矩陣為

(10)

式(10)中每一行表示同一個影響因素的不同判斷單元,每一列代表不同影響因素的同一個指標。因為每個判斷單元的含義并不完全相同,所以需要對原始數據進行標準化處理。

采用min-max標準化對數據進行歸一化,使原始數據統一落入[0,1]區間內,目的是保證數據的可靠性,減少由于數據相差太大,導致數值較高的指標影響過大,從而削弱數據較小的指標的作用。

假設給定了k個指標X1,X2,…,Xk,其中Xi={x1,x2,…,xn}。假設對各指標數據標準化后的值為Y1,Y2,…,Yk,那么

(11)

根據標準化之后的數據Y1,Y2,…,Yk可以建立信息熵權的標準判斷矩陣F。

根據信息論中信息熵的定義,一組數據的信息熵為

(12)

其中

(13)

如果pij=0,則定義

(14)

根據式(12)得到的各指標信息熵值E1,E2,…,Ek,可以得到各指標的權重。信息熵權的公式為

(15)

式中:k是信息熵權值的總數;Ei是信息熵。

2 考慮間隙和構件柔性的2-RR&2-PR并聯機構運動特性仿真

為呈現不同條件下的仿真模擬,對2-RR&2-PR并聯機構在考慮間隙及綜合考慮間隙和構件柔性的情況下進行仿真分析。其中,在考慮間隙的條件下又分別對間隙大小及間隙個數對機構的影響進行分析。考慮到實際機構工作情況,設轉速為240 r/min,同時為了提高仿真效率,將仿真模擬時間設為3 s,仿真步數設為3 000,本文所用到的求解器為GSTIFF求解器。

為了仿真時更加精確地計算出間隙處的接觸力,根據1.2節中建立的接觸力模型,采用Fortran語言編寫程序,將此程序加載到ADAMS函數求解庫中,用來求解間隙處兩元素間作用力,具體流程如圖3所示。將模型導入到ADAMS中,利用ADAMS的二次開發功能,通過ADAMS/Solver模塊編譯生成的動態鏈接庫文件(*dll),然后利用ADAMS/View模塊調用自定義函數,即生成的動態鏈接庫文件,就可將此自定義函數應用到模型中,再進行模型仿真,其中接觸力模型參數取值如表2所示。

圖3 自編子程序的實現Fig.3 Implementation of self-programmed subroutines

表2 接觸力模型參數Tab.2 Contact force model parameter

此程序代碼的核心內容是將間隙處兩元素的運動學特征和旋轉副自身的幾何特征綜合在一起,包含了間隙處兩元素在運動過程中無接觸、點接觸、線接觸、面接觸四種模式的描述。在仿真過程中通過該程序中的SYSARY子函數從模型中提取系統的狀態參數,完成接觸模式實時判別,從而依據不同接觸模式下的解析表達,計算出接觸力。

2.1 仿真數值計算

數值計算結果影響了仿真數據的準確性。為了確保仿真結果的準確性,有必要分析數值分析的誤差,以證明分析結果是收斂的,并且當時間和步長變化時,計算結果不會改變,因此增加了數據分析過程。

圖4顯示了t=0.9 s時軸銷中心點在X方向上的位移。圖4(a)表明,在相同積分階數條件下,結果在10次迭代后趨于穩定,計算效率最高。值得注意的是,當迭代次數小于3時,仿真無法滿足收斂條件,導致仿真失敗。圖4(b)表明,在相同的迭代次數條件下,結果在積分階數達到2階后趨于穩定。為了確保后續數據的準確性和計算效率,6階、10次迭代組合是最佳選擇。

(a) 不同迭代次數

(b) 不同積分階數圖4 軸銷中心點在X方向上的位移Fig.4 Displacement of shaft pin center point in X direction

2.2 間隙對機構的影響

2.2.1 間隙大小對機構的影響

將間隙處大小設置成0.04 mm、0.06 mm、0.08 mm和0.10 mm,利用MATLAB將仿真加速度數據進行短時傅里葉變換處理,繪制成不同間隙大小下加速度幅值的變化曲線如圖5所示。

(a) c=0.04 mm

(b) c=0.06 mm

(c) c=0.08 mm

(d) c=0.10 mm圖5 間隙大小對機構動平臺加速度的影響Fig.5 Influence of clearance size on the acceleration of the mechanism moving platform

從圖5可以看出,引入間隙后,機構動平臺的加速度有劇烈的波動,表現為明顯的非線性現象。隨著間隙值從0.04 mm增大到0.06 mm,含間隙機構加速度的最大幅值從197增加到了379,間隙值從0.06 mm增加到0.10 mm時,機構加速度的最大幅值又降低到了114。由此可以看出,在一定范圍內隨著間隙值的增大,間隙對加速度的影響也有所增大。但當間隙增加到一定值后,隨著間隙值的增大,間隙對加速度的影響會趨于減小。機構在運動時,間隙使得轉動副處產生接觸碰撞力,并使這些接觸碰撞力產生波動,從而導致了加速度的變化。

2.2.2 間隙個數對機構的影響

以含有兩個關節間隙的并聯機構為研究對象,研究間隙個數對機構的影響。通過仿真得到間隙個數對機構動平臺加速度的影響如圖6所示,其中兩個關節間隙的尺寸都為0.1 mm。當機構中含有一個間隙時,機構加速度的幅值最大如圖5(d)所示為114。當機構中含有兩個間隙時,機構加速度的幅值最大為183。

圖6 間隙個數對機構動平臺加速度的影響Fig.6 Influence of the number of clearances on the acceleration of the mechanism moving platform

通過圖6可以看出,不同間隙個數的機構,其加速度產生波動的規律并不相同,且相對于機構只含一個間隙的情況下,當機構含有兩個間隙時,動平臺加速度的最大幅值有所增大,說明雙間隙情況下機構動平臺明顯波動得更加劇烈。產生這種現象的原因是間隙與間隙之間的相互耦合作用,當存在多個間隙時,間隙之間相互影響,使得加速度明顯增大,意味著產生了更大的接觸碰撞力,機構更難趨于穩定。

2.3 綜合考慮間隙和構件柔性對機構的影響

剛體在受力后不會產生任何變形,而實際中任何物體在受到力的作用后,都會產生變形,為了更精確地揭示機構運動狀態,需要考慮構件的柔性。

將圖1中2-RR&2-PR并聯機構的A3A5、A4A6兩構件柔性化處理,通過ANSYS和ADAMS軟件聯合生成構件柔性化模型,將構件三維模型導入到ANSYS中,在ANSYS中創建單元和材料,設置材料密度屬性的參數選擇如表3所示。

表3 材料參數選擇Tab.3 Material parameter selection

通過劃分網格,在構件兩端的孔處設置剛性連接區域,確定關鍵點,建立剛性連接區域,從而生成模態中性文件(MNF文件)。MNF文件包含了ADAMS中柔性體的所有信息,在ADAMS中直接導入MNF文件,從而生成機構的剛柔耦合模型。

在2.2節中間隙尺寸分別為0.04 mm、0.06 mm、0.08 mm和0.10 mm,針對不同間隙尺寸下的模型,分別導入構件的MNF文件,將構件柔性化,建立含間隙和構件柔性耦合機構的仿真模型。

通過仿真得到綜合考慮間隙和構件柔性情況下對機構動平臺加速度的影響圖。為了更好地揭示給機構帶來的影響,選取間隙尺寸為0.06 mm情況下機構動平臺加速度的影響圖如圖7所示。

圖7 綜合考慮間隙和構件柔性對機構動平臺加速度的影響Fig.7 Comprehensive consideration of the influence of clearance and component flexibility on the acceleration of the mechanism moving platform

由圖5(b)可以得到只考慮間隙時,動平臺加速度的最大幅值為379,而綜合考慮間隙和構件柔性時,由圖7可以得到,動平臺加速度的最大幅值為20,相比之下大幅度減小。

圖8為仿真得到的綜合考慮間隙和構件柔性對間隙鉸鏈處接觸力的影響圖,其中間隙尺寸為0.06 mm,從圖8可以看出,只考慮間隙時最大接觸力約為50 N,而綜合考慮間隙和構件柔性的情況下接觸力有所減小,約為17 N,其中各個時間點所顯示的接觸力沒有明顯的規律,代表機構運動時間隙鉸鏈處具有非線性。而接觸力有所減小是因為當在關節間隙處的兩元素發生碰撞時,柔性構件可以起到明顯的緩沖作用,從而對接觸力和加速度偏差進行補償,大幅度緩和了間隙處兩元素的碰撞。同時可以看出間隙和構件柔性對機構運動及性能的影響不能相互抵消,所以對機構進行分析時有必要綜合考慮間隙和構件柔性。

圖8 綜合考慮間隙和構件柔性對間隙鉸鏈處接觸力的影響Fig.8 Comprehensive consideration of the influence of clearance and component flexibility on the contact force of the mechanism moving platform

3 2-RR&2-PR并聯機構定量分析

現有的評估機構非線性行為的定量分析方法主要是Lyapunov指數法,該方法通過Lyapunov指數的正負來判斷系統是否具有混沌特性,但沒有具體分析不同因素對機構非線性的影響程度,需要一種新的定量分析方法,在此基礎上去衡量不同因素對機構非線性的影響程度。本文提出一種基于熵權法的非線性定量分析新方法,通過計算各影響因素的相對比例,根據熵權法的定義計算出各影響因素下的熵權值,以此評估不同條件下間隙和構件柔性對機構非線性的影響程度。

3.1 Lyapunov指數法定量分析

Lyapunov指數法中正指數表示系統具有混沌特性,其中只需求解最大Lyapunov指數即可。選取機構動平臺加速度為計算對象進行混沌分析,在MATLAB中計算得到

(16)

式中:tw為時間窗口;tau為時間延遲;m為嵌入維數;P為平均周期。

將所得數據代入最大Lyapunov指數程序,計算結果如圖9(a)所示。選擇虛線矩形中所示的線性區域,用最小二乘法在該區域擬合直線。擬合直線的斜率為最大Lyapunov指數,如圖9(b)所示。顯然所得最大Lyapunov指數大于零,驗證了該系統具有混沌特性。

(a) 最初結果

(b) 線性擬合結果圖9 最大Lyapunov指數的估計Fig.9 Estimation of the maximum Lyapunov exponent

3.2 基于熵權法的考慮間隙的機構定量分析

將前文2.2節中的仿真參數設為組別1~5,如表4所示。圖10為不同間隙條件下機構動平臺在X、Y、Z方向上加速度的仿真數據圖,利用MATLAB編寫熵權法程序,將仿真數據代入程序中,計算在不同間隙大小及間隙個數情況下的熵權值。

表4 機構仿真參數Tab.4 Mechanism simulation parameters

(a) c=0.04 mm

(b) c=0.06 mm

(c) c=0.08 mm

(d) c=0.10 mm

(e) c=0.10 mm圖10 不同間隙條件下X、Y、Z方向加速度變化圖Fig.10 The change of acceleration in X, Y and Z directions with different clearance conditions

如表5所示,為計算得到的不同間隙條件下機構動平臺X、Y、Z方向加速度熵權值。繪制不同間隙大小下(組別1~4)的熵權值變化圖如圖11所示,其中圖11(a)為折線圖,為了更直觀呈現間隙對加速度影響程度,繪制成雷達圖如圖11(b)所示。繪制不同間隙個數下(組別4~5)的熵權值變化雷達圖如圖12所示。

表5 不同間隙條件下的加速度熵權Tab.5 Entropy weight of acceleration with different conditions

(a) 熵權折線圖

(b) 熵權雷達圖圖11 不同間隙大小下的熵權變化圖Fig.11 Entropy weight change map for different clearance sizes

圖12 不同間隙個數下的熵權值變化雷達圖Fig.12 Radar plot of entropy weight change with different number of clearances

圖11(a)顯示了加速度熵權的變化趨勢,圖11(b)顯示了不同方向上加速度影響程度的變化。如圖11(a)所示,隨著間隙尺寸的變大,加速度在Z方向上的熵權幾乎不發生變化,一直趨于0,在X方向和Y方向上的熵權都會發生變化,其中加速度X方向上的熵權值遠大于Y方向的熵權值,代表X方向的加速度對間隙尺寸的變化比較敏感,對加速度影響權重最高。同時從第1組到第4組可以看出隨著間隙值的增大,X方向加速度熵權值先增大后減小,代表間隙對機構X方向加速度的影響先增大后減小,相應的Y方向加速度熵權值先減小后增大,代表間隙對機構Y方向加速度的影響先減小后增大。

圖12為不同間隙個數下的熵權值變化雷達圖,顯示了不同間隙個數下加速度影響程度的變化。從圖12可以得出隨著間隙個數的增加,X方向加速度熵權值增大,對加速度影響權重提高,代表對機構X方向加速度的影響逐漸增大,約為Y方向上的4倍。

3.3 基于熵權法的考慮間隙和構件柔性的機構定量分析

將前文2.3節中的仿真參數設為組別1和2,如表6所示,其中第一組的仿真參數與3.2節中的第二組一致,作為對照組。圖13為考慮間隙和構件柔性條件下機構動平臺在X、Y、Z方向上加速度的仿真數據圖,采用相同的方法計算在含間隙機構中引入構件柔性情況下的熵權值。

表6 機構仿真參數Tab.6 Mechanism simulation parameters

圖13 考慮間隙和構件柔性條件下X、Y、Z方向加速度變化圖Fig.13 The change of acceleration in X, Y and Z directions with considering clearance and component flexibility

如表7所示,為在含間隙機構中引入構件柔性情況下計算所得的熵權值,其中第一組為對照組,熵權值與3.2節中的第二組一致。將計算的結果繪制成雷達圖如圖14所示。

表7 考慮間隙和構件柔性條件下的加速度熵權Tab.7 Entropy weight of acceleration under consideration of clearance and component flexibility

圖14 考慮間隙和構件柔性下的熵權值變化雷達圖Fig.14 Radar plot of entropy weight change under consideration of clearance and component flexibility

圖14顯示了在含間隙機構中引入構件柔性情況下加速度影響程度的變化。從圖14可以看出,在考慮構件柔性的情況下X方向和Y方向上的熵權都會發生變化,其中加速度X方向上的熵權值遠大于Y方向的熵權值,代表X方向的加速度在考慮構件柔性時比較敏感,對加速度影響權重最高,且X方向加速度熵權值有所增大,代表對機構X方向加速度的影響程度增大,Y方向加速度熵權值有所減小,代表對機構Y方向加速度的影響程度減小。同時綜合表5和表7的數據可以發現,在任何條件下,動平臺X方向加速度熵權值最大,所占權重最高,而Z方向加速度熵權幾乎不發生變化,這也體現出本文所做的各組仿真之間的關聯性。

4 結 論

(1) 間隙對機構動平臺加速度有明顯的影響,間隙過大或間隙過小都會引起副元素劇烈的碰撞,隨著間隙個數的增多,會大大影響機構的穩定性,柔性構件對間隙的碰撞有明顯的緩沖作用。

(2) 熵權法可以定量分析出不同條件下各指標的權重,基于熵權法研究了不同條件下對機構非線性的影響程度,結果表明在任何條件下,動平臺X方向加速度熵權值最大,說明對非線性影響權重最高,為含間隙和構件柔性的機構非線性分析提供一種新的定量分析方法。

(3) 含間隙和構件柔性的并聯機構是一種復雜的非線性系統,很難進行精準的描述。通過熵權法定量分析不同因素對機構非線性影響程度的趨勢,使得對2-RR&2-PR并聯機構運動的描述更為準確,為含間隙機構非線性定量分析提供了理論參考。

猜你喜歡
方向影響模型
一半模型
是什么影響了滑動摩擦力的大小
2022年組稿方向
計算機應用(2022年2期)2022-03-01 12:33:42
2021年組稿方向
計算機應用(2021年4期)2021-04-20 14:06:36
哪些顧慮影響擔當?
當代陜西(2021年2期)2021-03-29 07:41:24
2021年組稿方向
計算機應用(2021年1期)2021-01-21 03:22:38
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
擴鏈劑聯用對PETG擴鏈反應與流變性能的影響
中國塑料(2016年3期)2016-06-15 20:30:00
主站蜘蛛池模板: 国产精品.com| 国产在线视频自拍| 成人精品亚洲| 亚洲成a人片在线观看88| 国产精品亚欧美一区二区三区| 国产精品美人久久久久久AV| 久久精品日日躁夜夜躁欧美| 国产青青操| 国产xx在线观看| 男女精品视频| 精品无码日韩国产不卡av| 精品国产成人av免费| 成年av福利永久免费观看| 麻豆AV网站免费进入| 一级爆乳无码av| 中文字幕丝袜一区二区| 美美女高清毛片视频免费观看| 国产精品一区在线麻豆| 波多野结衣无码中文字幕在线观看一区二区| 98超碰在线观看| 97视频免费在线观看| 91一级片| 一级毛片在线免费视频| 久久无码免费束人妻| 国产凹凸一区在线观看视频| 国产第一页亚洲| 久久激情影院| 欧美一级高清视频在线播放| 日本a级免费| 97视频在线观看免费视频| 久久久黄色片| 国产91色在线| 四虎精品免费久久| 91精品情国产情侣高潮对白蜜| 99久久精品国产自免费| 国产精品无码久久久久AV| 伊在人亞洲香蕉精品區| 亚洲91精品视频| 老司机久久精品视频| 有专无码视频| 狠狠干欧美| 亚洲天堂网在线观看视频| 欧美、日韩、国产综合一区| 久久国产精品嫖妓| 亚洲美女视频一区| 亚洲AV无码久久精品色欲| 亚洲精品视频在线观看视频| 精品一区国产精品| 国产91av在线| 国产成人精品日本亚洲77美色| 精品久久久久久中文字幕女 | 精品国产自在现线看久久| 国产精品视频999| 乱系列中文字幕在线视频| 婷婷六月色| 黄网站欧美内射| 亚洲一级色| 2020国产免费久久精品99| 亚洲色成人www在线观看| 亚洲精品无码高潮喷水A| 国产精品午夜福利麻豆| 精品久久蜜桃| 亚洲首页在线观看| 亚洲人成影院在线观看| 亚洲一区二区成人| 亚洲中文字幕在线一区播放| 精品福利视频导航| 亚洲中文字幕在线一区播放| 亚洲天堂视频在线播放| 国产人人干| 国产免费自拍视频| 综合久久五月天| www.精品视频| 国产在线无码av完整版在线观看| 四虎国产精品永久一区| 亚洲二区视频| 爆操波多野结衣| 国产一区二区人大臿蕉香蕉| 韩国自拍偷自拍亚洲精品| 午夜精品福利影院| 精品国产自在现线看久久| 米奇精品一区二区三区|