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

駐留型UUV錨泊系統(tǒng)運(yùn)動(dòng)建模與分析

2016-08-03 01:30:15張斌宋保維

張斌,宋保維

(西北工業(yè)大學(xué) 航海學(xué)院, 陜西 西安 710072)

?

駐留型UUV錨泊系統(tǒng)運(yùn)動(dòng)建模與分析

張斌,宋保維

(西北工業(yè)大學(xué) 航海學(xué)院, 陜西 西安 710072)

摘要:針對(duì)駐留型水下航行器(UUV)錨泊系統(tǒng)對(duì)水流作用的運(yùn)動(dòng)響應(yīng)問題,依據(jù)歐拉-伯努利梁理論,建立了包括彎矩作用在內(nèi)的錨鏈三維運(yùn)動(dòng)模型,并使用四元數(shù)代替歐拉角來描述錨鏈姿態(tài),以消除某些特殊情況下因錨鏈姿態(tài)大幅度變動(dòng)或個(gè)別歐拉角不確定性導(dǎo)致的運(yùn)動(dòng)方程奇異現(xiàn)象,然后通過適當(dāng)?shù)倪吔鐥l件,將UUV、錨鏈和錨塊的運(yùn)動(dòng)控制方程耦合起來,采用有限差分方法對(duì)系統(tǒng)耦合運(yùn)動(dòng)模型進(jìn)行數(shù)值離散處理通過牛頓-拉夫遜方法迭代求解整個(gè)錨泊系統(tǒng)的運(yùn)動(dòng)響應(yīng)。使用Hopland拖曳試驗(yàn)數(shù)據(jù)對(duì)模型進(jìn)行實(shí)例對(duì)比驗(yàn)證。結(jié)果表明:這種建模方法可以取得良好的準(zhǔn)確性與計(jì)算效率,在此基礎(chǔ)上模擬仿真得到了周期變動(dòng)水流作用下的UUV位置及姿態(tài)的響應(yīng)情況為錨泊系統(tǒng)的正常工作提供理論依據(jù)。

關(guān)鍵詞:水下航行器;錨泊系統(tǒng);動(dòng)態(tài)分析;有限差分方法;姿態(tài)四元數(shù);模型仿真

網(wǎng)絡(luò)出版地址:http://www.cnki.net/kcms/detail/23.1390.u.20160127.1137.030.html

在海洋資源探測(cè)、地形地貌測(cè)量及軍事反潛作戰(zhàn)中,需要UUV在某一特定海域進(jìn)行探測(cè)或監(jiān)測(cè)工作,為保證UUV維持其空間位置并保持良好的姿態(tài),必須時(shí)刻通過推進(jìn)器對(duì)UUV進(jìn)行控制,從而消耗掉大量能源,無法達(dá)到長(zhǎng)時(shí)間監(jiān)測(cè)或探測(cè)的目的。為了解決UUV續(xù)航能力差、工作時(shí)間短的局限性,通過錨鏈將UUV與錨塊連接起來,組成駐留UUV錨泊系統(tǒng),使UUV處于水下一定深度并維持其靜平衡狀態(tài)[1],減少能源消耗。為便于UUV水下探測(cè)或監(jiān)測(cè)工作的開展及其二次啟動(dòng)發(fā)射的實(shí)現(xiàn),對(duì)UUV錨泊狀態(tài)時(shí)的姿態(tài)有一定要求,必須研究系統(tǒng)在海流作用下的動(dòng)態(tài)運(yùn)動(dòng)響應(yīng)。

水下錨鏈、纜索廣泛應(yīng)用于拖曳及錨泊系統(tǒng)中,在非均勻海流作用下其動(dòng)態(tài)運(yùn)動(dòng)解析求解十分困難,只能采用數(shù)值方法得到時(shí)域內(nèi)的動(dòng)態(tài)運(yùn)動(dòng)響應(yīng),有限差分方法則是其中最為普遍采用的方法。文獻(xiàn)[2]對(duì)二維狀態(tài)下水面船舶的錨泊狀態(tài)動(dòng)態(tài)響應(yīng)進(jìn)行了求解分析;文獻(xiàn)[3-4]分別研究了水下錨系導(dǎo)彈發(fā)射系統(tǒng)在海流作用下的三維運(yùn)動(dòng)情況以及通訊纜索對(duì)水下航行器運(yùn)動(dòng)的影響,但都沒有考慮彎矩對(duì)錨鏈姿態(tài)的影響;文獻(xiàn)[5-6]將纜索彎曲剛度考慮在內(nèi),建立了三維狀態(tài)下拖曳系統(tǒng)的耦合運(yùn)動(dòng)模型,并采用有限差分方法求解。

與拖曳系統(tǒng)數(shù)學(xué)模型相比,錨泊系統(tǒng)的數(shù)學(xué)模型與之類似卻有不同之處。兩者不僅具有截然不同的邊界條件,而且更值得注意的是,錨泊系統(tǒng)與海流之間的相對(duì)速度較低,海流非均勻性造成的錨鏈流體動(dòng)力非線性更加顯著,錨鏈易呈現(xiàn)低應(yīng)力狀態(tài),錨鏈彎曲剛度的影響就變得不可忽視。另外,由于錨泊系統(tǒng)的下端點(diǎn)固定,UUV由于具有較大的運(yùn)動(dòng)慣性,其運(yùn)動(dòng)狀態(tài)的改變往往滯后于錨鏈,這使得某些特殊情況下錨鏈姿態(tài)大幅度變化、個(gè)別歐拉角具有不確定性,會(huì)導(dǎo)致的運(yùn)動(dòng)方程奇異。

本文采用有限差分方法建立包含彎矩在內(nèi)的錨鏈運(yùn)動(dòng)數(shù)學(xué)模型,使用四元數(shù)代替歐拉角,避免數(shù)值求解過程產(chǎn)生奇異,并通過耦合條件把駐留UUV、錨鏈與錨塊系統(tǒng)綜合考慮,建立適合駐留UUV錨泊系統(tǒng)的動(dòng)態(tài)運(yùn)動(dòng)方程,對(duì)錨泊系統(tǒng)不同條件下的動(dòng)態(tài)運(yùn)動(dòng)響應(yīng)進(jìn)行預(yù)報(bào)分析。

1錨泊系統(tǒng)三維運(yùn)動(dòng)數(shù)學(xué)模型

1.1坐標(biāo)系選擇

為便于分析錨泊系統(tǒng)在不規(guī)則海流作用下的動(dòng)態(tài)運(yùn)動(dòng)響應(yīng),需要建立兩個(gè)直角坐標(biāo)系:地面坐標(biāo)系SE(O0,x0,y0,z0)、航行器體坐標(biāo)系SB(O,x,y,z)。如圖1所示,地面坐標(biāo)系原點(diǎn)選取在錨點(diǎn)O0處,O0x0軸位于水平面內(nèi),O0y0軸位于豎直面內(nèi),鉛直向上為正,O0z0軸的指向參照右手系規(guī)則確定;航行器體坐標(biāo)系原點(diǎn)位于UUV浮心O處,Ox軸沿航行器縱軸且向后為正,Oy軸垂直于Ox軸并指向上方,Oz軸垂直于Oxy平面。另外,引入局部坐標(biāo)系ST(t,n,b),其中軸t表示錨鏈切向,方向?yàn)殄^鏈長(zhǎng)度s的收縮方向;n為錨鏈的法向;b為錨鏈的副法線方向且位于水平面內(nèi)。錨塊沉于海底并與錨鏈下端固聯(lián),錨鏈上端在錨系點(diǎn)處與UUV連接,構(gòu)成完整的錨泊系統(tǒng)。

圖1 錨泊系統(tǒng)坐標(biāo)系示意圖Fig.1 Coordinate systems of mooring system

坐標(biāo)系之間可通過相應(yīng)的姿態(tài)角相互轉(zhuǎn)換,歐拉角(ε,γ)為錨鏈微元相對(duì)于地面坐標(biāo)系的姿態(tài)角。其中ε為方位角,即錨鏈偏離X0軸的角度,而γ為抬升角,如圖2所示。它們可以定義為

(1)

(2)

式中:dl為錨鏈微元的長(zhǎng)度,dx0、dy0、dz0為微元在地面坐標(biāo)系下的3個(gè)分量。通過坐標(biāo)系旋轉(zhuǎn),可以獲得局部坐標(biāo)系到地面坐標(biāo)系的轉(zhuǎn)換矩陣:

圖2 錨鏈微元姿態(tài)角Fig.2 Attitude angles of the cable element

1.2姿態(tài)四元數(shù)與歐拉角映射關(guān)系

為描述錨鏈微元的空間姿態(tài),可以采用式(1)~(2)中描述的歐拉角方法或者采用四元數(shù)及其導(dǎo)出形式。歐拉角方法中采用空間旋轉(zhuǎn)造成的人為割裂給姿態(tài)描述及姿態(tài)控制帶來了便利,但其本身也存在許多不足,其過多的三角運(yùn)算影響計(jì)算速度與精度,同時(shí)也不適用于大幅度的姿態(tài)運(yùn)動(dòng)描述,在特定情況下,運(yùn)動(dòng)方程出現(xiàn)奇異現(xiàn)象[7]。為此,本文引入哈密爾頓四元數(shù)方法描述錨鏈微元的空間姿態(tài)。

由文獻(xiàn)[7-8]可知,定點(diǎn)運(yùn)動(dòng)剛體由某一位置到另一位置的任意有限轉(zhuǎn)動(dòng)可以用四元數(shù)轉(zhuǎn)動(dòng)表示。從而局部坐標(biāo)系到地面坐標(biāo)系的轉(zhuǎn)換矩陣為

式中:哈密爾頓四元素q0、q1、q2、q3與姿態(tài)角ε、γ具有如下關(guān)系:

1.3錨鏈運(yùn)動(dòng)學(xué)模型

本文建立的數(shù)學(xué)模型中,假設(shè)錨鏈為連續(xù)的細(xì)長(zhǎng)圓柱狀纜索,材質(zhì)均勻且具有各向同性,在整個(gè)錨鏈長(zhǎng)度上光滑連續(xù)。錨鏈微元的運(yùn)動(dòng)控制方程為[9]

(3)

(4)

式中:F為局部系下的錨鏈微元受力,Q為局部系下的錨鏈微元內(nèi)部力矩,H為錨鏈微元的流體阻力,W為慣性系下單位長(zhǎng)度錨鏈濕重,Ma為包括附加質(zhì)量在內(nèi)的單位長(zhǎng)度錨鏈質(zhì)量矩陣,r為局部系中的錨鏈微元位移向量,J為單位長(zhǎng)度錨鏈的轉(zhuǎn)動(dòng)慣量,ω為錨鏈微元的角速度向量??刂品匠讨械倪\(yùn)算符號(hào)定義如下,表示對(duì)括號(hào)內(nèi)部的變量進(jìn)行等式右側(cè)的計(jì)算:

由于錨鏈沿長(zhǎng)度s方向具有二階光滑連續(xù)性,可以得到

(5)

由坐標(biāo)系之間的轉(zhuǎn)換關(guān)系以及錨鏈微元的轉(zhuǎn)動(dòng)角速度、曲率的物理含義可知:

(6)

(7)

綜合式(3)~(7),把錨鏈微元的運(yùn)動(dòng)控制方程分解到局部系t、n、b方向并寫為矩陣形式[10]:

(8)

式中:

為錨鏈微元的控制變量;M、N均為10×10矩陣,矩陣M中的非零元素如下所示:

M(5,10)=-SγVb,M(6,10)=SγVn-CγVt

矩陣N中的非零元素如下所示:

N(4,1)=1/EA,N(5,9)=1+Ft/EA,

式中:下標(biāo)t、n、b表示局部系中各軸向,Ct、Cn為錨鏈的切向和法向阻力系數(shù),A為錨鏈橫截面積,E為錨鏈的楊氏彈性模量,ρw為海水密度,ρc為錨鏈密度,U為錨鏈微元相對(duì)海流速度。

1.4邊界條件

要求解方程(8),必須給定錨鏈的初始形態(tài)值以及錨鏈上下兩端點(diǎn)處的邊界條件,形成封閉的系統(tǒng)運(yùn)動(dòng)模型。

錨鏈下端點(diǎn)處與錨塊通過鉸接方式相連接,在該點(diǎn)處彎矩為零,且錨點(diǎn)固定不動(dòng),可知

(9)

(10)

錨鏈上端點(diǎn)處與航行器鉸接耦合,且速度與航行器耦合點(diǎn)處速度相同:

(11)

(12)

2錨泊系統(tǒng)運(yùn)動(dòng)模型數(shù)值求解方法

本文選用有限差分方法在時(shí)間和空間上對(duì)錨鏈微元的控制方程(8)~(12)離散化。通過n+1個(gè)節(jié)點(diǎn)將錨鏈劃分為n段長(zhǎng)度為Δs的微元,并將時(shí)間劃分為一系列時(shí)間步長(zhǎng)Δt。應(yīng)用具有二階精度的中心差分格式,將方程(8)在中間節(jié)點(diǎn)j+1/2和時(shí)間節(jié)點(diǎn)i+1/2處展開,可得到

(13)

應(yīng)用Newton-Raphson方法迭代求解差分得到的非線性方程組(13),即可得到錨泊系統(tǒng)的動(dòng)態(tài)運(yùn)動(dòng)響應(yīng)情況。

3模型驗(yàn)證與數(shù)值計(jì)算結(jié)果

3.1模型驗(yàn)證

為驗(yàn)證本文建立的包含彎矩作用在內(nèi)的錨鏈數(shù)值模型精確性,選取Hopland拖曳實(shí)驗(yàn)[11]與本文模擬數(shù)值進(jìn)行對(duì)比。試驗(yàn)參數(shù):錨鏈直徑0.033 2m,錨鏈總長(zhǎng)360m,錨鏈密度3 121kg/m3,錨鏈彈性模量77.5GPa。拖曳過程中,拖船速度歷時(shí)60s,由2.5kn線性降至1.0kn。不同時(shí)刻錨鏈姿態(tài)如圖3所示??梢钥吹?,模型數(shù)值模擬得到的錨鏈姿態(tài)動(dòng)態(tài)值與實(shí)驗(yàn)值吻合良好,錨鏈阻力系數(shù)的不確定性造成了模擬值與實(shí)驗(yàn)值存在細(xì)微偏差。

圖3 拖船減速過程錨鏈形態(tài)變化Fig.3 Cable configurations during the towing ship decelerating

3.2數(shù)值計(jì)算結(jié)果分析

基于上述的錨泊系統(tǒng)運(yùn)動(dòng)數(shù)學(xué)模型,仿真UUV與錨鏈在海流作用下的運(yùn)動(dòng)過程。本文涉及的駐留航行器具有回轉(zhuǎn)體外形,錨泊系統(tǒng)的主要參數(shù)如表1所示。UUV運(yùn)行至距海底高度10m處釋放錨鏈,錨塊及錨鏈在重力作用下下落,實(shí)現(xiàn)錨泊駐留功能,初始時(shí)刻錨鏈姿態(tài)呈豎直狀態(tài)。為模擬海流的復(fù)雜緩慢變化,引入如下正弦波模型模擬海流速度J(T為海流運(yùn)動(dòng)周期):

駐留系統(tǒng)的動(dòng)態(tài)響應(yīng)過程如圖4~9所示。

表1 動(dòng)態(tài)運(yùn)動(dòng)仿真過程主要參數(shù)

圖4 UUV受錨鏈拉力隨時(shí)間變化曲線Fig.4 Curve of cable tension acting on the UUV in time domin

圖5 錨鏈姿態(tài)隨時(shí)間變化曲線Fig.5 Shape of mooring cable under the influence of current

圖6 UUV攻角與側(cè)滑角隨時(shí)間變化曲線Fig.6 Curves of UUV attack angle and sideslip angle changing with time

從仿真結(jié)果可以看出,在海流作用下,錨泊系統(tǒng)位置發(fā)生偏移,由于錨塊固聯(lián)于海底,且駐留UUV相對(duì)于錨鏈具有較大的運(yùn)動(dòng)慣性,故而其運(yùn)動(dòng)狀態(tài)變化滯后于錨鏈,如圖4所示,錨鏈中間部分在錨泊運(yùn)動(dòng)初期呈凸起狀,隨著時(shí)間推移,航行器海流方向速度逐漸增加,錨鏈隨之呈現(xiàn)懸鏈狀。

由于海流波的周期性特征,由圖所示,航行器的姿態(tài)偏降、速度、攻角及錨鏈拉力均以海流波動(dòng)周期波動(dòng),并且比海流提前一定的相位角;航行器的俯仰角在海流及輔助推進(jìn)器的作用下,被限制在[5°,-3°]范圍內(nèi),可以保證良好的駐留姿態(tài),便于水下工作的開展。

圖7 UUV深度隨時(shí)間變化曲線Fig.7 Curve of UUV depth changing with time

圖8 UUV俯仰角角隨時(shí)間變化曲線Fig.8 Curve of UUV pitching angle changing with time

圖9 UUV速度隨時(shí)間變化曲線Fig.9 Curves of UUV velocity changing with time

4結(jié)論

1)與試驗(yàn)數(shù)據(jù)對(duì)比后發(fā)現(xiàn),本文所建立的錨鏈三維運(yùn)動(dòng)數(shù)學(xué)模型不僅能有效避免數(shù)值求解過程中產(chǎn)生的奇異現(xiàn)象,同時(shí)還具有較高的計(jì)算精度,能準(zhǔn)確地反映出錨鏈等拖曳系統(tǒng)與錨泊系統(tǒng)中常用的水下纜索張力傳遞及其自身形態(tài)變化情況。

2)駐留UUV的水下錨泊姿態(tài)與運(yùn)動(dòng)響應(yīng)是整個(gè)錨泊系統(tǒng)設(shè)計(jì)的關(guān)鍵內(nèi)容,對(duì)系統(tǒng)的運(yùn)動(dòng)響應(yīng)進(jìn)行仿真分析后發(fā)現(xiàn),航行器俯仰角這一關(guān)鍵參數(shù)在海流作用下逐漸增大,在采用輔助推進(jìn)器進(jìn)行姿態(tài)調(diào)節(jié)后能明顯將其限制在安全范圍內(nèi),避免發(fā)生系統(tǒng)走錨或航行器觸底現(xiàn)象。同時(shí),這些系統(tǒng)的參數(shù)變化規(guī)律可以為航行器姿態(tài)控制和穩(wěn)定性研究提供理論依據(jù),為實(shí)現(xiàn)錨泊系統(tǒng)正常工作及駐留UUV二次啟動(dòng)提供有意義的參考。

參考文獻(xiàn):

[1]朱信堯, 宋保維, 單志雄, 等. 海底定點(diǎn)停駐無人水下航行器流體動(dòng)力特性分析[J]. 上海交通大學(xué)學(xué)報(bào), 2012, 46(4):573-578.

ZHU Xinyao, SONG Baowei, SHAN Zhixiong, et al. Hydrodynamic characteristics analysis of UUV parking on the seabed[J]. Journal of Shanghai Tiao Tong University, 2012, 46(4): 573-578.

[2]HUO Cunfeng, YAO Baoheng, FU Bin, et al. Investigation on transient dynamic behaviors of low-tension undersea cables[J]. Journal of Shanghai Tiao Tong University: Science, 2011, 16(1): 34-39.

[3]邵成, 艾艷輝, 代軍. 水下錨系導(dǎo)彈發(fā)射系統(tǒng)運(yùn)動(dòng)研究[J]. 兵工學(xué)報(bào), 2011, 32(9): 1154-1158.

SHAO Cheng, AI Yanhui, DAI Jun. Study on motion of underwater towed missile launch system[J]. Acta armamentarii, 2011, 32(9): 1154-1158.

[4]FENG Z, ALLEN R. Evaluation of the effects of the communication cable on the dynamics of an underwater flight vehicle[J]. Ocean engineering, 2004, 31(8/9): 1019-1035.

[5]PARK H I, JUNG D H, KOTERAYAMA W. A numerical and experimental study on dynamics of a towed low tension cable[J]. Applied ocean research, 2003, 25(5): 289-299.

[6]BURGESS J J. Bending stiffness in a simulation of undersea cable deployment[J]. International journal of offshore and polar engineering, 1993, 3(3): 197-204.

[7]REDOUANE D, ADIL S, HICHAM M. Euler and quaternion parameterization in VTOL UAV dynamics with test model efficiency[J]. International journal of applied information systems, 2015, 9(8): 25-28.

[8]KATSUKI S, SEBE N. Rotation matrix optimization with quaternion[C]//Proceedings of the 10th Asian Control Conference. Kota Kinabalu: 2015: 1-6.

[9]BUCKHAM B J. Dynamics modelling of low-tension tethers for submerged remotely operated vehicles[D]. Victoria: University of Victoria, 2003: 60-64.

[10]SRIVASTAVA V K. Analyzing parabolic profile path for underwater towed-cable[J]. Journal of marine science and application, 2014, 13(2): 185-192.

[11]VAZ M A, PATEL M H. Transient behaviour of towed marine cables in two dimensions[J]. Applied ocean research, 1995, 17(3): 143-153.

收稿日期:2015-01-13.

基金項(xiàng)目:國(guó)家自然科學(xué)基金項(xiàng)目(51179159).

作者簡(jiǎn)介:張斌(1989-),男,博士研究生; 宋保維(1963-),男,教授,博士生導(dǎo)師. 通信作者:張斌, E-mail: dyjzhangbin@163.com.

doi:10.11990/jheu.201501019

中圖分類號(hào):TJ63

文獻(xiàn)標(biāo)志碼:A

文章編號(hào):1006-7043(2016)04-0498-05

Dynamic modeling and simulation of mooring system for an unmanned underwater vehicle

ZHANG Bin, SONG Baowei

(School of Marine Science and Technology, Northwestern Polytechnical University, Xi′an 710072, China)

Abstract:To determine the kinematic performance of an unmanned underwater vehicle (UUV) attached to a mooring line and lurking on the seabed, a three-dimensional cable mathematical model that considers the effects of bending moments was established based on the Euler-Bernoulli beam theory. In addition, a quaternion-based cable attitude model was adopted as a substitute for the traditional Euler-angle form to eliminate the singular behavior under some special circumstances, namely, a drastic change in cable attitude or the existence of some specific uncertain Euler angles. The governing equations of the UUV, cable, and anchor were integrated by using appropriate boundary conditions to obtain the translational and rotational motion equations of the mooring system. Thereafter, the mathematical model of the mooring system was discretized by using the finite difference method, and the Newton-Raphson iterative method was employed to solve the difference equations. Application data from Hopland's towing experiment were extracted to validate the mathematical model. The results show that the modeling algorithm is accurate and efficient. Then, the UUV's position deviation and attitude change were simulated in periodically changing current to provide a theoretical principle for maintaining its normal working state.

Keywords:unmanned underwater vehicle; mooring system; dynamic analysis; finite-difference method; attitude quaternion; model simulation

網(wǎng)絡(luò)出版日期:2016-01-27.

主站蜘蛛池模板: 国产真实二区一区在线亚洲| 天堂成人av| 91在线视频福利| 国产人妖视频一区在线观看| 国产熟睡乱子伦视频网站| 欧美特级AAAAAA视频免费观看| 无码精品国产VA在线观看DVD| AV不卡国产在线观看| 国产一级视频在线观看网站| yjizz视频最新网站在线| 黄色三级毛片网站| 丁香综合在线| 激情综合网激情综合| 一区二区三区精品视频在线观看| 亚洲啪啪网| 在线视频亚洲色图| 国产精品3p视频| 日韩无码黄色网站| 国产凹凸一区在线观看视频| 国产在线视频自拍| 91精品国产一区自在线拍| 狼友av永久网站免费观看| 国产美女精品人人做人人爽| 欧美视频在线第一页| 狠狠亚洲五月天| 热re99久久精品国99热| 亚洲精品在线91| 欧美日本在线| 久久综合成人| 中文字幕 欧美日韩| 91系列在线观看| 毛片网站观看| 亚洲国产精品无码久久一线| 国产精品污污在线观看网站| 国产白浆视频| 欧美一级在线播放| 永久成人无码激情视频免费| 91久久国产成人免费观看| 国产国产人成免费视频77777| 国模私拍一区二区| 香蕉国产精品视频| www中文字幕在线观看| 精品福利网| 九九免费观看全部免费视频| 国产人碰人摸人爱免费视频| 色噜噜狠狠狠综合曰曰曰| 国产91视频观看| 亚洲美女操| 国产精品熟女亚洲AV麻豆| 亚洲欧美在线精品一区二区| 亚洲精品无码久久久久苍井空| 伊人久久婷婷五月综合97色| 国产欧美精品一区二区| 亚洲v日韩v欧美在线观看| 中国国产一级毛片| 久久婷婷综合色一区二区| 伊人久久大香线蕉影院| 中文字幕第1页在线播| 国产一区二区人大臿蕉香蕉| 国产精品高清国产三级囯产AV| 97国产精品视频人人做人人爱| 国产成人午夜福利免费无码r| a级毛片在线免费| 青青草欧美| 免费国产黄线在线观看| 日本人妻一区二区三区不卡影院 | 日韩成人免费网站| 免费国产不卡午夜福在线观看| 国产在线自乱拍播放| 日韩欧美中文亚洲高清在线| 成人国产免费| 亚洲av成人无码网站在线观看| 日韩av无码DVD| 精品国产自在在线在线观看| 91免费国产高清观看| 精品黑人一区二区三区| 精品国产免费观看| 四虎成人免费毛片| 精品国产成人高清在线| 国产成人精品男人的天堂| 国产精鲁鲁网在线视频| 欧美yw精品日本国产精品|