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

基于模型預測控制的高速柔性并聯機構振動控制

2014-09-05 02:02:48胡俊峰張憲民徐貴陽
振動與沖擊 2014年1期
關鍵詞:模態振動模型

胡俊峰, 張憲民, 徐貴陽

(1. 江西理工大學 機電工程學院,江西 贛州 341000; 2. 華南理工大學 機械與汽車工程學院, 廣州 510640)

并聯機構具有高速、高精、高承載能力等特點,在許多領域得到應用。為實現并聯機構的高速精密運行,必須考慮其構件變形,對其進行振動控制。柔性并聯機構為分布參數的非線性、時變、多輸入多輸出系統,這給振動控制帶來較大困難。由于壓電陶瓷等智能材料和現代控制理論的發展,采用智能材料作傳感器和作動器的振動主動控制得到了應用[1-6]。

但由于柔性并聯機構的非線性因素、各參數之間的耦合、模型參數的不確定性和時變性等影響,很難精確建立系統的數學模型,這給其振動主動控制帶來困難。模型預測控制具有以下基本特征:模型預測、滾動優化和反饋校正。模型預測控制方法對系統模型要求低且形式靈活,能夠實時補償因系統參數時變和非線性因素等引起的控制誤差,有效克服系統不確定因素的影響,控制效果好且魯棒性強。

本文以一種考慮構件的彈性變形的新型兩自由度平動并聯機構為研究對象,采用有限元法和實驗模態方法建立其離散狀態空間模型,研究應用模型預測控制策略設計魯棒振動控制器。

1 系統動力學模型

應用有限元法,建立含有壓電陶瓷的柔性并聯機構的彈性動力學方程[3-4]

(1)

采用模態技術處理,忽略高階模態的影響,引入變換

U=ψηc

(2)

其中,ψ為振型矩陣,ηc為受控模態的振型坐標,將式(2)代入式(1)可得:

(3)

(4)

其中,ω1,…,ωc為系統前c階固有頻率,ζ1,…,ζc為系統前c階阻尼比,Dac、Das分別為系統作動器、傳感器分布矩陣,它們是與作動器、傳感器位置和模態振型有關,Nc為在模態坐標下的廣義力,稱為模態力。

定義受控狀態變量為

(5)

則式(4)可以寫成狀態空間方程形式

(6)

式中,

(7)

2 模型預測控制器設計

模型預測控制結構如圖1所示。在模型預測控制中,將系統輸出的動態預估問題分為預測模型的輸出預測和基于偏差的預測校正兩部分。式(6)所表示的預測模型只是對柔性并聯機構動態特性的粗略描述,而實際系統中存在非線性、時變性和隨機干擾等因素,因此,預測模型不可能與實際對象完全相符,預測模型的輸出與實際系統輸出之間必然存在偏差。采用這種偏差進行在線校正,使系統構成具有負反饋環節的系統,從而提高預測控制系統的魯棒性,以滿足柔性并聯機構的魯棒振動控制。

圖1 模型預測控制結構

2.1 系統的預測模型

由于模型預測控制是基于離散狀態方程進行預測,并且采用計算機進行控制,能方便存取各時刻的采樣信號,便于計算。首先對系統連續狀態方程(6)進行離散化,離散控制狀態方程可寫成如下形式:

(8)

式中,設采樣周期為Ts,k表示kTs時刻,Ψ(k)、Γ(k)、Θ(k)、Ω(k)分別表示由連續狀態方程轉換為離散狀態方程系數矩陣,Nc(k)為不可測擾動,并且考慮測量噪聲對輸出的影響,則測量輸出ym(k)為

ym(k)=Ω(k)Xc(k)+m(k)

(9)

其中,m(k)為測量噪聲。不可測量干擾Nc(k)的狀態空間模型為

(10)

(11)

X(k+1)=AX(k)+BuVin+Bdnd(k)
y(k)=CX(k)+Ddnd(k)

(12)

為了求解預測控制問題,在狀態空間表達式中必須計算被控變量的預測值。它可由當前狀態的最佳估計值以及假設的未來輸入計算,或等價地可由最近的輸入V(k-1)的值以及假設的未來輸入變化ΔV(k+i|k)來計算。為了表述簡單,假設從時刻k=0開始,根據式(12)可以預測模型的未來輸出值,第i個時間步的輸出可以表達為:

(13)

其中,p為預測時域,i代表預測時域p的第i時間步,h為控制時域。式(13)也可以寫成如下形式

(14)

其中,

(15)

(16)

2.2 狀態估計器設計

由于狀態變量為模態位移和模態速度,它們是不能測量的,所以需要用一個狀態估計器估計這些狀態量。狀態估計器是基于如圖2所示的模型建立的,狀態估計器是為了估計受控對象、輸入干擾對象的狀態,使用卡爾曼濾波器技術設計。狀態估計量可表示為:

(17)

(19)

2.3 滾動優化

模型預測控制可以通過解下面的優化問題得到,設預測控制最優化目標函數為:

minJ(ΔV)

(20)

其中,J為系統性能指標,它可以表示為:

Vtarget(k+i)]TRV[V(k+i|k)-Vtarget(k+i)]}(21)

其中,ΔV為設計變量,它可以表示為ΔV=[ΔV(0),…,ΔV(p-1)]T,Q是ny×ny矩陣,RΔV和RV分別是懲罰是設計變量和控制輸入電壓的系數矩陣,它們是Na×Na矩陣,它們均為半正定矩陣,Na為作動器的個數,矩陣Q用于懲罰被控輸出的預測值y(k+i+1|k)與傳感器期望輸出r(k+i+1)之間的偏差,Vtarget(k+i)是控制輸入電壓設定值。

圖2 用于狀態估計的模型

受控系統存在的約束有控制電壓及其變化率的約束條件和系統輸出量約束條件:

Vjmin(i)≤Vj(k+i|k)≤Vjmax(i)
ΔVjmin(i)≤ΔVj(k+i|k)≤ΔVjmax(i)
yjmin(i)≤yj(k+i|k)≤yjmax(i)

(22)

其中,下標( )j表示向量的第j個分量,(k+i|k)表示基于時刻k的值的時刻k+i的預測值,r(k)是輸出參考值。

將式(21)展開可以表達為如下形式:

(23)

WV=blkdiag(RV,…,RV)
WΔV=blkdiag(RΔV,…,RΔV)
Wy=blkdiag(Q,…,Q)

(24)

式中,blkdiag(·)代表塊對角矩陣,每一個矩陣的塊重復p次。

由于控制輸入信號V(k)與未來輸入變化ΔV(k)的關系可以表示為:

V(k)=V(k-1)+ΔV(k)

(25)

將式(14)代入式(23),則最終優化目標函數可以寫成如下形式:

(26)

式中,Const表示常數,KΔV、Kr、Ku、KVt和KX分別為目標函數的系數矩陣。

同理,約束條件(22)也可寫成如下形式:

(27)

將式(14)和式(25)代入上式,約束條件可以寫成如下形式:

MΔVΔV≤Mlim+MVV(-1)+MXX(0)

(28)

式中,MΔV、Mlim、MV、MX為約束條件系數矩陣,它們可以根據約束條件的上下邊界獲得。

聯合式(26)和式(28)的最優化問題即是一個二次規劃優化問題,利用二次規劃方法求解可以得到ΔV的最優解ΔV*,它可以表示為:

(29)

3 作動器與傳感器位置優化

由于振動響應的控制和測量大多采用離散分布式的作動器/傳感器,它們在機構上的位置和數目優化問題對系統可控性/可觀性有很大影響。合理的作動器布置可以在較小輸入能耗作用下最大限度的抑制振動,否則可能使系統性能惡化,甚至導致系統不穩定。相應地,不合理的傳感器配置,可能導致無法對系統性能做出正確的評定。

作動器/傳感器位置優化的目的是通過合理選取作動器的位置和數目以確定作動力分布矩陣,從而盡可能大地影響被控模態,以及合理的選取傳感器的位置和數目以確定傳感器的輸出系數向量,使得傳感器測量的響應中盡可能多包含各階的模態分量,同時應能保證閉環控制系統的性能穩定。由于現有的位置優化準則過于復雜,采用一種表征作動能量的可控性指標和表征觀測信號能量的可觀性指標,以確定作動器和傳感器的最優位置。

3.1 確定作動器位置

壓電作動器作動力列向量可表示為[10]

fc=BcVin

(30)

(31)

對作動力系數矩陣Bc進行奇異值分解可得

Bc=USaVT

(32)

將式(32)代入式(31)可得:

(33)

設Vineq=VTVin,由于VTV=INa,則

(34)

如果控制輸入電壓Vin給定,向量Vineq的值也可確定,即

(35)

基于σi在式(35)所示控制能量中的作用,引入可控性指標

(36)

由式(36)可知,控制電壓確定后,指標Ωa正比于作動器提供的能量,指標Ω值越大,作動器提供的能量越大,而作動能量又與分布矩陣Da、模態矩陣有關。因此,對于給定的控制輸入電壓,根據Ωa值的大小可以確定作動器的位置,Ωa值最大的位置即為放置作動器的最優位置。

3.2 確定傳感器位置

(37)

同理,對系數矩陣Cc進行奇異值分解,則性能指標Jc可表示為

(38)

(39)

4 實 例

為了驗證所設計的控制器的有效性,對一新型2自由度并聯機構進行振動主動控制研究。考慮到高速機構振動控制的復雜性,對機構進行了簡化,并為了簡化模型,設計如圖3所示的并聯機構[7-8],其簡圖和系統坐標系統如圖4所示。l1和l2分別為主動臂和從動臂的長度,l1=505 mm,l2=430 mm。由于從動臂相對于主動臂為細長桿件,為了研究方便,將機構的主動臂視為剛性構件,從動臂為柔性構件。板條由鋁合金制成,長360 mm,寬60 mm,厚2 每一柔性從動臂用梁單元進行模擬。系統單元劃分如圖5所示,圖中Ni(i=1~21)為結點編號,Ej(j=1~20)為單元編號,從動臂 1、2均劃分為10單元,系統共劃分為20個梁單元,考慮邊界約束,系統共有79個廣義坐標。系統廣義坐標U是在如圖5所示的坐標系統XOY定義的,系統廣義坐標為結點自由度,每個結點具有縱向位移、橫向位移、彈性轉角和曲率4個自由度,具體的有限元模型建立過程參見文獻[9]。

4.1 作動器與傳感器位置

設每個柔性桿粘貼3個壓電作動器,對機構前4階模態進行控制,并且作動器的尺寸已確定。下面確定放置作動器最佳位置。如圖4所示,綜合考慮彈性構件的長度和作動器不能重疊,設作動器位置為下面4種組合,第1種組合為分別作用在單元3、5、7、13、15、17;第2種組合為作用在單元3、5、8、13、15、18;第3種為3、6、8、13、16、18;第4種為4、6、8、14、16、18;作動器在4 種不同位置的性能指標如圖6所示,它們分別為64.775 2、103.209 0、 93.635 4、82.081 9。由此可知,第2種組合的可控性指標最大,所以選擇該組合,即作動器放置在如圖8所示的單元E3、E5、E8、E13、E15、E18上。

對放置傳感器的位置進行20種組合,對應這20種組合的性能指標值如圖7所示,從該圖可以發現第7種和第18種組合,性能指標值較大,這兩種組合放置傳感器的位置單元分別為單元3、5、8、13、15、18和單元5、6、8、15、16、18,恰好其中第7種組合為放置作動器的位置,而第18組合也有兩個位置是放置作動器的位置,所以作動器和傳感器采用同位配置較好。A1~A3表示對應于放置在單元E3、E5和E8的作動器,S1~S3表示放置在在單元E3、E5和E8中點的傳感器。同樣,A4~A6表示為放置在從動臂2上的作動器,對應于放置在如圖所示的單元E13、E15和E18,S4~S6表示為放置在從動臂2上的傳感器,對應于放置在單元E13、E15和E18中點。根據作動器的位置,并且考慮壓電陶瓷片的長度為50 mm,為方便建模,設每一片壓電陶瓷片放置位置為一梁單元,這樣每一片壓電陶瓷片施加在構件的力矩可視為加在單元節點上外力矩,可以設置每根從動臂10個單元的長度分別為65 mm、15 mm、50 mm、40 mm、50 mm、40 mm、40 mm、50 mm、15 mm、65 mm。

圖3 實驗配置

圖4 機構簡圖及其坐標系統

圖5 系統單元劃分

圖6 作動器位于不同位置的性能指標值

圖7 傳感器位于不同位置的性能指標值

圖8 作動器與傳感器位置

4.2 實驗模態測試

由于機構的阻尼比難以通過有限元方法獲得,需要借助試驗模態測試得到。采用脈沖錘擊法進行實驗模態測試。用力錘對機構進行敲擊,產生一個寬頻帶的激勵,如圖9所示,固定敲擊位置,測量12個不同位置的加速度信號。力錘是由PCB公司生產,型號為Model 086C03,測量范圍為0~500 lbf,靈敏度系數為2.44 mV/N。加速度傳感器采用Kistler 8690C50型壓電式加速度計,它可以測量垂直于測量點的加速度信號。為了消除噪聲干擾,采用多次平均,設每個測點的測量的次數為5次。使用ZonicBook/618E得到激勵點和各測量點的時間歷程數據,利用eZ-Analyst軟件求出各測點的頻響函數。采用ME’scopeVES對這些頻響數據進行曲線擬合,得到擬合后的頻響曲線如圖10所示,得到系統的固有頻率和阻尼比。

圖9 加速度傳感器布置

圖10 頻響擬合曲線

圖11 振動控制試驗原理圖

由試驗得到的前2階固有頻率、阻尼比和有限元計算得出的固有頻率如表1所示,從該表可以看出,有限元計算值和實驗值相對誤差接近于13-19%,則說明采用有限元建立的模型和實際系統還存在一定的差距。造成這一誤差的原因主要有以下幾點:①機構在鉸鏈處的約束復雜;②機構在鉸鏈處存在間隙;③有限元方法采用梁單元,簡化對實際構件模型的建立;④建立模型時沒有考慮螺栓和螺母之類的構件。所以,應用于控制器設計的模型和實際系統存在誤差,也就是說,理論模型不精確。

4.3 實驗驗證

為了實現機構的高速運動,采用伺服電機驅動控制對象,機構的剛體運動軌跡規劃參見文獻[8]。交流伺服電機型號為Panasonic MHMD042P1V,其配套的伺服驅動器型號為MBDDT2210003。機構是在高速的運動過程中執行操作,它的兩個關節的角速度和角加速度都較高,圖12(a)、(b)是以操作時間最短為目標,對機械手的關節空間進行最優運動軌跡規劃得到的兩個主動關節的角速度和角加速度,它們最大值分別接近15 rad/s、555 rad/s2。

圖12 關節角速度和角加速度

電阻應變計為由中航電測公司生產型號為溫度自補償型BE120-3AA(11),該電阻應變計標稱阻值為120 Ω,靈敏度系數為2.17。動態信號采集系統為ZonicBook/618E,8通道動態信號輸入,每通道分辨率為16-bit,最大采樣率1 MHz。功率放大器是由哈爾濱芯明天科技公司生產,型號為X-505.00,放大倍數為15,由于實驗使用了6對壓電陶瓷片,采用6路功率放大器分別驅動。dSPACE實時仿真系統是由德國dSPACE公司開發的一套基于MATLAB/Simulink的控制系統開發及半實物仿真的軟硬件工作平臺。本試驗使用dSPACE DS1103平臺和MATLAB/Simulink搭建振動控制系統。振動控制試驗原理如圖11所示,電阻應變片與動態電阻應變儀通過1/4橋路連接,應變儀將被測點的應變信號轉換成電壓信號傳送至數據采集卡,采集到的數據以實時方式傳送給dSPACE DS1103處理,按照所設計的控制器實時計算所需的控制電壓,實時計算得到的數字量經過數模轉換模塊輸出,由于D/A轉換模塊輸出電壓范圍為-10~+10 V,需經電壓放大器放大后施加給壓電陶瓷片,完成對系統的控制。實時控制系統設計過程為:首先利用MATLAB/Simulink構建控制系統框圖;然后,利用Real-Time Workshop技術將Simulink框圖程序生成實時代碼并下載到dSPACE快速原型機中;最后,使用dSPACE提供的綜合試驗與測試環境軟件ControlDesk實現試驗過程控制和參數在線修改和實時數據采集。為了實現機構高速運動過程中的實時振動控制,采用如圖11所示的工作原理,通過dSPACE同時接受來自伺服電機和應變片傳感器的信號,并根據該信號進行運動控制和振動控制,通過所設計的控制器同時發出剛體運動控制信號和振動控制信號,這樣,機構高速運動控制與振動控制就能實現同步,以滿足振動控制的實時性要求。

在不精確的理論模型的基礎上應用魯棒模型預測控制設計振動控制器。在設計控制器中,模態力視為白噪聲信號,其幅值同樣分別設為20.5、20.4。量測噪聲設為幅值為6×10-6的白噪聲。為了與實時控制采樣周期相同,設模型預測控制器的采樣周期Ts=0.001 s,預測時域p=8,控制時域h=2,加權矩陣為Q=diag(105,105,105,105,105,105),Vtarget(k)=0,RV=0,RΔV=diag(0.01,0.01,0.01,0.01,0.01,0.01),為了減小機構振動響應,應變片傳感器期望輸出r(k)=0,約束條件為控制輸入電壓的允許最小和最大值,即Vmin(k)=-150V,Vmax(k)=150V,控制輸入電壓在一個時間步的最小和最大變化率設為ΔVmin(k)=-300V/Ts,Vmax(k)=300V/Ts。

表1 系統固有頻率和阻尼比

圖13 在無控制和有控制兩種情況下傳感器S1~S6的輸出應變

圖13(a)、(b)、(c)、(d)、(e)、(f)分別為在無控制和有控制兩種情況下傳感器S1~S6的輸出應變,分別比較6個傳感器的輸出應變變化,可以得知,在模型預測控制器作用下,傳感器的輸出應變以較快的速度衰減,也就是說,并聯機構的彈性振動響應得到較大的抑制,這說明控制器的有效性和魯棒性。圖14(a)、(b)、(c)、(d)、(e)、(f)分別為在模型預測控制器作用下六個作動器的輸入電壓。由于模型預測控制器是基于離散狀態方程模型設計的,控制器的輸出電壓為離散值,控制器的輸出電壓如圖14所示,從該圖可以看出,控制輸入電壓均滿足約束條件,同樣說明在允許的控制輸入電壓下,所設計的控制器能獲得較好的性能。控制電壓在大部分時間步均達到了約束條件的邊界,說明控制器需要消耗較多的能量才能達到預期的控制效果。

圖14 六個作動器的輸入電壓

5 結 論

考慮高速柔性并聯機構的非線性和不確定性,采用模型預測控制理論設計控制器抑制其彈性振動響應。根據系統不精確的理論模型推導系統的預測模型。將不確定性外部擾動和量測噪聲視為白噪聲,采用Kalman濾波估計器估計系統狀態變量,以控制電壓及其變化率為約束條件,將系統性能指標和約束條件化為一個標準二次規劃優化問題,通過求解這一優化問題來得到最優控制輸出。采用一種表征作動能量的可控性指標和表征觀測信號能量的可觀性指標,以確定作動器和傳感器的最優位置。試驗驗證了模型預測控制策略對柔性并聯機構振動響應控制的有效性和魯棒性,取得較為滿意的控制效果。

參 考 文 獻

[1]宋軼民,余躍慶,張 策,等. 柔性機器人動力學分析與振動控制研究綜述 [J]. 機械設計,2003,20(4):1-5.

SONG Yi-min, YU Yao-qing, ZHANG Ce, et al. Summary on the study of dynamics analysis and vibration control of flexible robots [J]. Machine Design, 2003, 20(4):1-5.

[2]Sun D, Mills J K, Shan J, et al. A PZT actuator control of a single-link flexible manipulator based on linear velocity feedback and actuator placement [J]. Applied Acoustics, 2004, 14(4):381-401.

[3]胡俊峰,張憲民. 一種新型兩自由度柔性并聯機械手的主動振動控制 [J]. 中國機械工程, 2010, 21(17):2017-2024.

HU Jun-feng, ZHANG Xian-min. Active vibration control and its simulation of a novel 2-DoF flexible parallel manipulator [J].China Mechanical Engineering, 2010, 21(17):2017-2024.

[4]胡俊峰,張憲民, 朱大昌,等. 柔性并聯機器人動力學建模 [J]. 農業機械學報, 2011, 42(11):208-213.

HU Jun-feng, ZHANG Xian-min, ZHU Da-chang et al. Dynamic modeling of flexible parallel robot [J]. Transactions of the Chinese Society for Agricultural, 2011, 42(11):208-213.

[5]Gasparetto A, Zanotto V. Vibration reduction in a flexible-link mechanism through synthesis of an optimal controller [J]. Meccanica, 2006, 41(6):611-622.

[6]Iorga L, Baruh H, Ursu I. A Review of H∞robust control of piezoelectric smart structures [J]. Appl. Mech. Rev., 2008, 61(4):1-16.

[7]張憲民, 袁劍鋒. 一種二維平動兩自由度平面并聯的機器人機構:中國, CN1903521 [P]. 2007-01-31.

[8]Hu J F, Zhang X M, Zhan J Q. Trajectory planning of a novel 2-DoF high-speed planar parallel manipulator [A]. In:1st International Conference on Intelligent Robotics and Applications, ICIRA 2008 [C]. Wuhan, China, 2008:199-207.

[9]胡俊峰,張憲民. 兩自由度高速并聯機械手的彈性動力學分析 [J]. 華南理工大學學報(自然科學版),2009, 37(11):123-128.

HU Jun-feng, ZHANG Xian-min. Elastodynamic analysis of a novel 2-DOF high-speed parallel manipulator [J]. Journal of South China University of Technology(Natural Science Edition), 2009, 37(11):123-128.

[10]Guney M, Eskinat E. Optimal actuator and sensor placement in flexible structures using closed-loop criteria [J]. Journal of Sound and Vibration, 2008, 312(1-2):210-233.

猜你喜歡
模態振動模型
一半模型
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
重要模型『一線三等角』
振動與頻率
天天愛科學(2020年6期)2020-09-10 07:22:44
重尾非線性自回歸模型自加權M-估計的漸近分布
中立型Emden-Fowler微分方程的振動性
3D打印中的模型分割與打包
國內多模態教學研究回顧與展望
基于HHT和Prony算法的電力系統低頻振蕩模態識別
UF6振動激發態分子的振動-振動馳豫
計算物理(2014年2期)2014-03-11 17:01:44
主站蜘蛛池模板: 亚洲三级a| 午夜福利免费视频| 99热这里只有成人精品国产| 国产日本一线在线观看免费| 香蕉在线视频网站| 亚洲第一视频免费在线| 久久国产精品娇妻素人| 亚洲人成影视在线观看| 国产成人久久777777| 国产高颜值露脸在线观看| 成人综合久久综合| 欧亚日韩Av| 免费观看三级毛片| 国产电话自拍伊人| 国产成人免费| 久久精品无码一区二区日韩免费| 国产精品福利尤物youwu| 91精品综合| 亚洲国产精品久久久久秋霞影院| 91黄视频在线观看| 中文字幕欧美日韩| 欧美日韩中文国产| 日本亚洲最大的色成网站www| 亚洲视频色图| 99r在线精品视频在线播放| 最新亚洲人成网站在线观看| 亚洲日韩精品伊甸| 99这里精品| jizz在线免费播放| 高清不卡毛片| AV在线天堂进入| 69综合网| 天堂网亚洲综合在线| 久草美女视频| 成AV人片一区二区三区久久| 尤物在线观看乱码| 亚洲久悠悠色悠在线播放| av在线5g无码天天| 国产精选小视频在线观看| 国产综合精品一区二区| 国产精品嫩草影院av | 欧美激情第一欧美在线| 美女高潮全身流白浆福利区| 无码乱人伦一区二区亚洲一| 丁香婷婷激情网| 国产在线精品美女观看| 国产亚洲精| 91丨九色丨首页在线播放| 亚洲va在线观看| 波多野结衣亚洲一区| 亚洲 欧美 中文 AⅤ在线视频| 国产精品部在线观看| 91在线日韩在线播放| 免费Aⅴ片在线观看蜜芽Tⅴ| 亚洲第一黄色网| 亚洲一区二区在线无码| 午夜高清国产拍精品| 久久久精品无码一区二区三区| 青青青视频91在线 | 久久久久人妻一区精品| 久久狠狠色噜噜狠狠狠狠97视色| 国产欧美在线观看一区| 国产免费久久精品99re丫丫一| 久草中文网| 国产亚洲欧美另类一区二区| 操国产美女| 国产一级毛片网站| 第一页亚洲| 久久久久久久久久国产精品| 国产日韩欧美成人| 亚洲狼网站狼狼鲁亚洲下载| 无码视频国产精品一区二区| 99热免费在线| 欧美性精品不卡在线观看| 国产嫩草在线观看| 伊人无码视屏| 19国产精品麻豆免费观看| 日本在线视频免费| 老司机久久精品视频| 无码久看视频| 国产剧情国内精品原创| 伊人精品成人久久综合|