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

具有解析式位置正解的三平移并聯機構設計與分析

2020-03-09 07:36:08沈惠平曾博雄尤晶晶許正驍楊廷力
農業機械學報 2020年2期

沈惠平 曾博雄 尤晶晶 李 菊 許正驍 楊廷力

(1.常州大學現代機構學研究中心, 常州 213016; 2.南京林業大學機械電子工程學院, 南京 210037)

0 引言

三自由度的三平移(3T)并聯機器人機構驅動元件少、結構緊湊、設計制造及控制成本較低,在工程應用中得到了廣泛應用。CLAVEL[1]于1988年發明Delta 3-DOF平移并聯機構,之后一些學者對Delta的衍生操作手進行了研究[2-4];TSAI等[5-6]提出一種移動副驅動、支鏈含4R平行四邊形機構的三自由度移動并聯機構;文獻[7-8]對3-RRC型三平移并聯機構進行了運動學和工作空間分析;KONG等[9]提出了一種三自由度3-CRR機構,該機構具有良好的運動性能,且沒有明顯的奇異位置;LI等[10-11]提出了3-UPU型三平移機構,并對該機構的瞬時運動學性能進行了分析;于靖軍等[12]基于螺旋理論對三維平動并聯機構構型進行了綜合分析;陸晶等[13]提出了一種3-RRRP(4R)三平移并聯機構,并進行了運動學和工作空間的分析;楊廷力等[14-16]基于單開鏈單元理論對3T0R型并聯機構進行了型綜合,綜合出多種新機型,并對它們進行了分類;ZHAO[17]考慮運動學的各向異性,對一種三平移并聯機構進行尺度綜合及其運動學研究;ZENG等[18-20]設計了一種三平移Tri-pyramid并聯機構,并對其位置方程的正反解、雅可比矩陣、各向同性等運動學特性進行了分析;PRAUSE等[21]對多種三平移并聯機構分別進行了維數綜合、邊界狀況、工作空間等特征的比較,并選出了性能較好的機構;MAHMOOD等[22]提出了一種三自由度3-[P2(US)]機構,并進行了運動學和靈巧度分析。

上述三平移并聯機構一般存在兩大問題:① 機構耦合度不為零,即κ≥1時,一般得不到解析式位置正解,而只能得到數值解,這對誤差分析、尺度綜合、剛度分析及動力學研究帶來較大的不便。② 不具有輸入-輸出運動解耦特性[23-24],這使運動控制及軌跡規劃等較為復雜。這兩者給應用帶來了困難。

當并聯機構具有解析式位置正解時,后續的研究內容易于進行:用輸入量表示全局奇異位形方程,進行操作度的性能評估及結構參數的全域優化;建立運動學誤差模型,并進行影響因素的敏感度分析;推導靈活工作空間的解析表達式;動力學方程的正解求解精度、效率與穩定性評估。具有解析式位置正解的并聯機構的拓撲設計與分析一直是機構學工作者不斷研究的方向之一,但目前設計的具有解析式位置正解的并聯機構拓撲類型較少。

本文根據基于方位特征(Position and orientation characteristics,POC)方程的并聯機構拓撲結構設計理論與方法[15-16],設計一種僅由移動副和轉動副組成且具有解析位置正解、部分運動解耦的新型三平移(3T)并聯機構,并對其位置正逆解、機構奇異性位形、工作空間以及速度與加速度等進行分析。

1 并聯機構設計及拓撲分析

1.1 機構的拓撲設計

本文提出的三平移機構,由動平臺1、靜平臺0以及2條混合支鏈(HSOC)組成,如圖1所示,其拓撲結構如下:

圖1 3T并聯機構簡圖

(1)右側混合支鏈Ⅰ由一個平面兩滑塊四轉動副6桿機構(簡稱2P4R平面機構)回路的中間構件11上,串聯兩個軸線相互平行的轉動副R3與R4,且R4副與動平臺1相連。顯然,混合支鏈Ⅰ上2P4R平面機構的中間桿11的輸出運動為兩平移一轉動(2T1R)。混合支鏈Ⅰ的拓撲結構等效地記為:

(2)左側混合支鏈Ⅱ由移動副P3與2個4R平行四邊形機構串聯而成,從P3副到動平臺1相連的平行四邊形RaiRbiRciRdi(i=1,2),分別記為①、②;其中,P3副與平行四邊形①在同一平面內剛性連接后,與平行四邊形②在其垂直平面內連接。混合支鏈Ⅱ上平行四邊形①的輸出桿S點的輸出運動為兩平移(2T),而平行四邊形②的輸出桿T點的輸出運動為三平移(3T)。混合支鏈Ⅱ的拓撲結構等效地記為:HSOC2{-P3-P(4R)-P(4R)-}。

(3)移動副P1、P2、P3與靜平臺0相連,P1與P2副為共軸線布置,且P1‖P3;機構運動時,2P4R平面機構與平行四邊形①的運動平面始終平行。

當靜平臺上的3個移動副以相同的速度運動時,該機構可實現大范圍的操作移動;而當其取不同的速度時,可實現小范圍內的精確作業。因此,該機構適合于長度方向大尺寸工件的機加工、噴涂、鉚接等操作工藝。

1.2 機構拓撲特性分析

1.2.1機構POC計算

并聯機構POC方程[15]為

(1)

(2)

式中MJi——第i個運動副的POC集

m——運動副數量

Mbi——第i條支鏈末端的POC集

n——支鏈條數

MPa——機構動平臺的POC集

取動平臺1上的任一點為基點O′,則由式(1)可得

由式(2)可確定動平臺POC集為

由此表明:動平臺1上任一點的POC集為三平移零轉動;即機構中的混合支鏈Ⅱ本身已實現三平移的設計要求,是它約束了混合支鏈Ⅰ末端構件的2個轉動輸出,從而使得動平臺僅產生三平移輸出。

1.2.2機構自由度計算

并聯機構的全周DOF公式[15-16]為

(3)

(4)

v=m-n+1

式中F——機構自由度

fi——第i個運動副的自由度

v——獨立回路數

ξLj——第j個獨立回路的獨立位移方程數

Mb(j+1)——第j+1條支鏈末端構件POC集

該機構包含兩個獨立回路,具體為:

顯然,其獨立位移方程數ξL1=3。

②(R3‖R4和混合支鏈Ⅱ)與上述2P4R平面機構組成第2個獨立回路,即

LOOP2{-R3‖R4-P(4R)-P(4R)-P3-}

由式(4)可得

由式(3)可得機構的自由度F為

因此,該機構自由度為3,當取靜平臺0上的移動副P1、P2、P3為驅動副時,動平臺1可實現3個平移的運動輸出。

1.2.3機構耦合度κ計算

由基于序單開鏈(SOC)單元的機構組成原理[15-16]知,一個機構可以分解為若干個最小子運動鏈(Sub- kinematic chain, SKC),每個SKC僅包含一個基本運動鏈BKC,所謂BKC是指DOF為零且其任意一個子運動鏈DOF大于零的最小運動鏈;而一個SKC又可分解為若干個單開鏈,第j個單開鏈SOCj的約束度定義為

(5)

式中mj——第j個SOCj的運動副數量

fi——第i個運動副的自由度

Ij——第j個SOCj的運動副數量

對一個SKC而言,需滿足

因此,其耦合度

(6)

κ揭示了機構回路運動變量之間的關聯、依賴程度;κ越大,機構的耦合性越強,運動學、動力學分析復雜度越高。

1.2.2節已計算出上述兩個回路的獨立位移方程數,分別為ξL1=3,ξL2=5;因此,兩個回路的約束度由式(5)分別計算得

于是,由式(6)得

這樣,該機構只包含一個SKC,且該SKC的κ=1;一般情況下,求解這類機構位置正解時,僅需在約束度為正值(Δj>0)的回路上設定一個虛擬變量;然后,在約束度為負值(Δj<0)的回路上建立一個含這個虛擬變量的位置約束方程,再通過一維搜索法可求出該虛擬變量的真實值,從而求得該機構的位置正解。

但由于此機構特殊的三平移方位特征約束,可直接通過約束度為負值(Δj<0)的回路作用于約束度為正值(Δj>0)的回路的幾何約束(即:桿11的運動始終平行于靜平臺),求出該虛擬變量,從而直接求得機構的解析式位置正解。

2 機構位置分析

2.1 基于SKC-SOC的機構位置正解求解原理

由式(5)可知,每個SKC可分解為一系列約束度為正值、零、負值的回路,因此,機構位置正解的求解,可轉換為該SKC內回路的位置求解,而3種回路的約束特性及其建模方法分別為:

2.2 坐標系的建立及參數標注

該機構的運動學建模如圖2所示,設靜平臺0為長、寬分別為2a、2b的矩形,以靜平臺0的幾何中心為原點,建立笛卡爾靜坐標系OXYZ,X、Y軸分別垂直、平行于A1A2連線,Z軸由右手法則確定;動坐標系O′X′Y′Z′原點位于動平臺1的中心,X′、Y′軸分別重合、垂直于D2F3連線,Z′軸由右手法則確定。

設3個驅動桿2的長度為l1,混合支鏈Ⅰ上連桿9、10的長度為l2,中間桿11、12的長度分別為l3、l4。

混合支鏈Ⅱ上平行四邊形短桿3、6的長度為l5,長桿4、7的長度為l6;平行四邊形之間的連接桿5的長度為l7,連接桿8的長度為l8;動平臺1上D2F3連線的長度為2d。

圖2 3T機構的運動學建模

設B1C1與Y軸正向的夾角γ為虛擬變量;D1D2、D3E3與X軸正向的夾角分別為α、β。

2.3 位置正解分析

該機構的位置正解求解歸結為:已知靜平臺0上3個點Ai(i=1,2,3)移動位置yA1、yA2、yA3,求動平臺1上O′點的坐標(x,y,z)。

(1)約束度為正的第1回路(LOOP1)的求解

LOOP1:A1-B1-C1-C2-B2-A2

在靜坐標系OXYZ中,易知點Ai、Bi(i=1,2,3)的坐標分別為A1=(-b,yA1,0)、A2=(-b,yA2,0)、A3=(b,yA3,0);B1=(-b,yA1,l1)、B2=(-b,yA2,l1)、B3=(b,yA3,l1)。

由1.2.1節可知,由于動平臺1三平移的特殊方位約束,機構運動過程中,2P4R平面機構的中間構件11始終平行于靜平臺0,即C1C2‖A1A2,則有

zC1=zC2

(7)

因此,點C1、C2的坐標分別為

C1=(-b,yA1+l2cosγ,l1+l2sinγ)
C2=(-b,yA1+l2cosγ-l3,l1+l2sinγ)

由幾何約束B2C2=l2,并整理、簡化有

ABcosγ+B2=0

當B=0時,無法求出虛擬變量γ的值,此時,機構發生輸出奇異,機構在運動過程中應當避免這種情況發生。

當B≠0時,可求出虛擬變量γ的值,此時有

(8)

其中A=2l2B=yA1-l3-yA2

這樣,約束度為負的第2回路作用于約束度為正的第1回路上的特殊幾何約束方程式(7),是直接求出虛擬變量γ解析解的關鍵,這是本機構拓撲結構具有的求出位置解析解的一個特點。γ求出后,機構第2回路上其余運動副的位置可容易求解。

(2)約束度為負的第2回路(LOOP2)的求解

LOOP2:D1-D2-F3-E3-D3-C3-B3-A3

由點C1、C2求出點D1、D2的坐標分別為

D1=(-b,yA1+l2cosγ-l3/2,l1+l2sinγ)D2=(-b+l4cosα,yA1+l2cosγ-l3/2,l1+l2sinγ+l4sinα)

同時,可計算得O′點的坐標為

O′=(x,y,z)=(-b+l4cosα+d,yA1+l2cosγ-l3/2,l1+l2sinγ+l4sinα)

(9)

進一步,點F3、E3、D3、C3的坐標用O′點的坐標分別表示為

(10)

由幾何約束B3C3=l6,并令

l4sinα-l6sinβ=t

(11)

則有

(H1+t)2-H2=0

(12)

由于機構運動時2P4R平面機構與平行四邊形①的運動平面始終平行,因此,恒有

yD1=yD3

(13)

l4cosα+2d-l6cosβ=2b

(14)

再從式(11)、(14)中消除β,則有

J1sinα+J2cosα+J3=0

(15)

其中

J1=2l4tJ2=4l4(b-d)

最后,將式(8)、(15)所求得γ、α代入式(9),即可得動平臺1上O′點的坐標(x,y,z)。

由式(8)知

γ=f1(yA1,yA2)

由式(15)知

α=f2(yA1,yA2,yA3)

因此,由式(9)知

即該機構具有部分輸入-輸出運動解耦性,這對動平臺的軌跡規劃和運動控制是有利的。

為理解方便,上述計算過程描述如圖3所示。

圖3 機構位置正解流程圖

Fig.3 Flow chart of direct position solutions

由此可知,幾何約束方程式(7),以及式(13)、(14)分別是求出本機構第1、2回路位置方程解析式的關鍵。

2.4 位置反解分析

該機構反解求解可歸結為:已知動平臺1上O′點的坐標(x,y,z),求靜平臺0上3個點Ai(i=1,2,3)移動位置yA1、yA2、yA3。

由式(9)可求出α為

(16)

再由式(14)可求出β為

(17)

進一步,求出點C1、C2的坐標分別為

C1=(-b,y+l3/2,z-l4sinα)
C2=(-b,y-l3/2,z-l4sinα)

另外,點C3的坐標已由式(10)給出,因此,由桿長約束條件,建立位置約束方程

(18)

即可求解yAi(i=1,2,3)為

(19)

綜上可知,當動平臺1上O′點的坐標(x,y,z)已知時,靜平臺0上3個點Ai(i=1,2,3)移動位置yA1、yA2、yA3各有兩組解。故逆解數為2×2×2=8,因此,該機構有8種構型。

2.5 位置正反解實例驗算

2.5.1正解算例

設該并聯機構結構參數為:a=350 mm,b=150 mm,d=50 mm,l1=30 mm,l2=280 mm,l3=140 mm,l4=180 mm,l5=90 mm,l6=230 mm。

設平行四邊形①與②之間連接桿5的長度l7=0,四邊形②與動平臺1之間的連接桿8的長度l8=0。此時,機構的CAD三維模型如圖4所示。

圖4 3T機構三維模型

取靜平臺0上3個點Ai(i=1,2,3)的位置為yA1=154.677 4 mm,yA2=-193.670 7 mm,yA3=31.061 1 mm。

由Matlab計算該機構位置正解,如表1所示。

表1 機構的位置正解數值

2.5.2逆解算例

將表1中第3組正解數值代入式(19),可得yAi(i=1,2,3)的8組逆解數值,如表2所示。

表2 機構的位置逆解數值

可見,表2中第3組的逆解數據和正解求解時給定的3個輸入位置yAi(i=1,2,3)一致,從而驗證了其正逆解的正確性。

3 并聯機構奇異位形分析

3.1 奇異位形分析方法

采用Jacobian法分析該機構的奇異位形。對式(16)、(17)的兩邊同時對時間t求一階導數,得

(20)

(21)

再將式(18)的兩邊對時間t求一階導數,并將式(20)、(21)代入求導后的等式,得

(22)

其中

u11=yC1-yB1u22=yC2-yB2
u33=yC3-yB3f11=cotα(zC1-zB1)
f12=yC1-yB1f13=zC1-zB1
f21=cotα(zC2-zB2)f22=yC2-yB2
f23=zC2-zB2f31=cotβ(zC3-zB3)
f32=yC3-yB3f33=zC3-zB3

因此,該機構動平臺末端執行器的輸出速度v1和Ai(i=1,2,3)輸入移動速度v2的關系為

Jpv1=Jqv2

(23)

其中

依據Jp、Jq是否奇異,將機構的奇異位形分為如下3類:①當det(Jq)=0時,機構發生輸入奇異。②當det(Jp)=0時,機構發生輸出奇異。③當det(Jq)=det(Jp)=0時,機構發生綜合奇異。

3.2 奇異位形分析

3.2.1輸入奇異

當機構發生輸入奇異,意味著每條支鏈靠近驅動桿的兩根桿處于折疊在一起或完全展開狀態。這時,動平臺自由度減少。此時det(Jq)=0,方程解的集合W為

W={W1∪W2∪W3}

(24)

其中

W1={yC1-yB1=0}(A1、B1、C1三點共線)
W2={yC2-yB2=0}(A2、B2、C2三點共線)
W3={yC3-yB3=0}(A3、B3、C3三點共線)

滿足W3的三維構型如圖5所示。

圖5 輸入奇異位形示例

3.2.2輸出奇異

當機構發生輸出奇異,意味著每條支鏈靠近動平臺的桿處于折疊在一起或完全展開的狀態,此時的動平臺自由度數增多,即使鎖住輸入,動平臺也可能存在自由度輸出。設

[fi1fi2fi3]=ei(i=1,2,3)

若det(Jp)=0,則向量e1、e2、e3有如下兩種情況:

(1)存在兩個向量線性相關

若e1=ke2,即滿足(f11,f12,f13)=k(f21,f22,f23),其三維構型為向量lB1C1、lB2C2在空間內平行,如圖6所示。

圖6 輸出奇異位形示例

若e1=ke3,即滿足(f11,f12,f13)=k(f31,f32,f33),此時有cotα(zC1-zB1)=kcotβ(zC3-zB3),zC1-zB1=k(zC3-zB3)。則有cotα=cotβ。

因本機構設定的桿長尺寸在機構運動過程中,恒有α≠β,即cotα≠cotβ;則e1≠ke3,同理可知e2≠ke3。

(2) 3個向量線性相關

設e2=k1e1+k2e3(k1k2≠0),此時有

(f21,f22,f23)=k1(f11,f12,f13)+k2(f31,f32,f33)

通過Matlab計算表明,該種情況下k1、k2的解無法解出,因此,此種情況不存在。

3.2.3綜合奇異

此時det(Jq)=det(Jp)=0,即輸入奇異和輸出奇異同時發生。在此位形下,機構的驅動關節和末端執行器都存在著瞬時互不影響的非零輸入和輸出,對應的位姿就是第3類奇異,處于該類奇異時,機構將失去自由度,在機構設計階段應予以避免。

上述奇異位置的求得,對樣機調試時的軌跡規劃與運動控制,具有參考價值。

4 機構工作空間分析

并聯機構的可達工作空間,是指在考慮運動副的轉角范圍、桿長不干涉的情況下,末端執行器的工作區域是衡量并聯機構性能的一個重要指標。本文采用極限邊界搜索法分析該并聯機構的工作空間,即首先根據桿長來設定工作空間的搜索范圍;然后,基于位置反解,搜索所有滿足約束條件的點,由這些點組成的三維圖即為該并聯機構的工作空間。

確定空間三維搜索范圍為:-135 mm≤x≤85 mm,-300 mm≤y≤300 mm,180 mm≤z≤480 mm,通過Matlab編程,得到該并聯機構的三維工作空間,根據第①類奇異判別式det(Jq)=0可得到輸入奇異軌跡,由第②類奇異判別式det(Jp)=0可得到輸出奇異軌跡,如圖7所示,其中,綠色部分為無奇異的工作空間,藍色部分為輸入奇異區域,紅色部分為輸出奇異區域,表明該工作空間內部存在較大的無奇異工作空間。

圖7 工作空間及其奇異情況

本文所設定的桿長參數之一是:2a=2l2+l3,在這個設定條件下,若考慮3個驅動副的移動范圍,即隨著3個驅動副的同速度移動,則圖7所示的機構工作空間相應地可沿Y方向延伸,具有與2PPPaR并聯機構所述的倒置型長方形柱的工作空間[25],適合于長度方向較大尺寸工件的機加工、噴涂、鉚接等操作工藝。

圖8為工作空間在XOY、XOZ和YOZ平面上的投影圖。

圖9為工作空間內Z方向上的4個X-Y截面圖,表明:z≥300 mm時,隨著z的增大,X軸方向上的移動距離逐漸減少,Y軸方向上的移動距離逐漸增加。

圖8 工作空間在XOY、XOZ和YOZ平面上的投影

圖9 工作空間內不同的X-Y截面圖

5 速度與加速度分析

5.1 速度公式

當機構處于非奇異位置時,矩陣Jp可逆,由式(23)可得

(25)

式(25)即為動平臺O′點的輸出速度。

5.2 加速度公式

對式(23)的兩邊同時對時間t求一階導數,得

(26)

其中

K=[K1K2K3]T

當機構處于非奇異位置時,矩陣Jp可逆,由式(26)可得

(27)

式(27)即為動平臺O′點的輸出加速度。

5.3 速度、加速度算例驗算

取3個驅動副的輸入移動運動規律分別為:yA1=154.677 4-50sint,yA2=-193.670 7-50sint,yA3=31.061 1+50sint。則其輸入移動速度、加速度的變化規律分別為

將式(25)、(27)導入Matlab軟件編程計算動平臺O′點的速度與加速度,分別得到速度、加速度曲線;同時,將該并聯機構的三維模型,通過SolidWords導入ADAMS軟件中進行仿真,設置仿真時間為10 s,仿真步長為0.1 s,得到101個時間點3個方向的速度與加速度后,將其導入Matlab軟件得速度與加速度如圖10、11所示。

圖10 動平臺速度曲線

圖11 動平臺加速度曲線

對比圖10和圖11可以發現:① 運用Matlab對式(25)、(27)進行編程得到的速度與加速度曲線與運用ADAMS仿真得到的曲線完全吻合,從而驗證了所推導的速度與加速度公式的正確性。② 在給定本文所設定的驅動副的輸入運動規律情況下,該機構動平臺的速度、加速度曲線變化連續平穩、無突變峰值,表明該機構具有較好的運動平穩性。

6 結論

(1)設計的三平移(3T)并聯機構具有4個優點:僅由移動副和轉動副組成,制造簡單、安裝方便;具有解析式位置正解,這對誤差分析、尺度綜合、剛度分析及動力學研究帶來較大的方便;具有部分輸入-輸出運動解耦性,這對機構的軌跡規劃及運動控制十分有利;具有大的操作工作空間,適合于長度方向較大尺寸工件的機加工、噴涂、鉚接等操作工藝。

(2)該機構通過第1回路輸出構件始終保持水平位置這一特殊的幾何約束條件(該條件由約束度為負的第2回路提供),直接求得整個機構位置正解解析解,這是本機構拓撲結構不同于其他機構的優點所在。

(3)基于位置反解,得到了該機構三類奇異位形發生的條件及其位置,給出了機構的工作空間大小及其奇異區域。

(4)機構速度與加速度仿真曲線變化較平穩、連續,具有較好的動態特性。

主站蜘蛛池模板: 国产福利在线观看精品| 真实国产精品vr专区| 免费在线观看av| 亚洲大尺码专区影院| 高清亚洲欧美在线看| 亚洲性网站| 不卡国产视频第一页| 色吊丝av中文字幕| 成人小视频在线观看免费| 九九九精品视频| 最近最新中文字幕在线第一页| 一级香蕉视频在线观看| 亚洲视频二| 自拍偷拍欧美日韩| 亚洲国产天堂久久综合226114| 波多野结衣亚洲一区| 国产喷水视频| 伊人久久精品亚洲午夜| 欧美日韩一区二区三区在线视频| 欧美午夜一区| 黄片在线永久| 国产精品手机在线观看你懂的| 亚洲精品欧美日本中文字幕| 美女无遮挡免费视频网站| 亚洲第一香蕉视频| 亚洲免费黄色网| 国产高颜值露脸在线观看| 亚洲综合色婷婷| 久久国产高潮流白浆免费观看| 亚洲欧美一级一级a| 欧美日韩国产综合视频在线观看 | 亚国产欧美在线人成| 国产尤物jk自慰制服喷水| 精品夜恋影院亚洲欧洲| 在线精品视频成人网| 亚洲午夜福利精品无码| 99免费视频观看| 91精品免费高清在线| 欧美成一级| 色婷婷亚洲综合五月| 国产美女无遮挡免费视频| 亚洲有无码中文网| 制服丝袜亚洲| 免费高清a毛片| 无码一区18禁| 亚洲男人天堂久久| 欧美亚洲欧美| 亚洲男人天堂久久| 亚洲视频四区| 青青热久免费精品视频6| 亚洲swag精品自拍一区| 国产精品久久久久久久久kt| 国产精品林美惠子在线观看| 不卡色老大久久综合网| 精品视频91| 亚洲综合香蕉| 亚洲日本www| 欧美在线视频a| 国产产在线精品亚洲aavv| 国产熟睡乱子伦视频网站| 久久熟女AV| 一级在线毛片| 国产女人18水真多毛片18精品| 亚洲黄色视频在线观看一区| 国产三区二区| 欧美激情视频二区| 成人在线观看一区| 欧美中文字幕在线视频| 女同国产精品一区二区| 波多野结衣久久高清免费| 91久久国产综合精品| 色哟哟国产精品| 国产精品亚洲一区二区三区z| 精品国产网| 欧美日韩久久综合| 亚洲a免费| 国产成人午夜福利免费无码r| 91年精品国产福利线观看久久| 亚洲熟妇AV日韩熟妇在线| 九九久久精品国产av片囯产区| 操美女免费网站| 国产精品成人观看视频国产|