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

變角度水面矢量推進器性能分析

2018-03-13 02:04:28呂建剛張仲志趙正龍劉金華
農業機械學報 2018年2期

呂建剛 張仲志 宋 彬 趙正龍 劉金華

(陸軍工程大學車輛與電氣工程系, 石家莊 050003)

0 引言

水陸兩棲平臺能夠擴展人類的活動方式,并提高水上與陸地農業生產的轉換效率。在現行大方形系數陸面汽車結構的基礎上進行改進,能夠維持陸上性能、節省結構設計成本。然而,已量產的VW166型、Amphicar Model 770型等大方形系數兩棲平臺在水中排水航行[1],受到的興波阻力與速度的6次方成正比,航速受限難以超過15 km/h[2];Aquada、Watercar等兩棲平臺的速度有所提升,但是結構已向船艇發展,方形系數大大減小而偏離了車輛結構[3]。

近年來,普通水面平臺不斷發展,為大方形系數兩棲平臺的改進提供了參照。其中,尾部阻流板設計通過增加托舉力而減小船體濕面積,進而降低了雙體船的阻力[4];水翼能夠增加托舉力并抬升船艇前端,進而改變流場獲得轉矩,實現縱傾,平臺阻力下降了20%以上[5-6];滑行面設計的高速艇滑水航行,由水動力提供托舉力和轉矩,將艇體托出水面并維持一定的縱傾角,大大提高了航速[7]。因此,增加托舉力和轉矩驅動,是水面平臺減阻提速的有效辦法[6-8]。但是將上述研究用于基于汽車結構的兩棲平臺,則改動較大而嚴重影響陸上性能,且增加了一套獨立的噴水或螺旋槳推進系統,結構復雜化。

南美叢林中的蛇怪蜥蜴踩踏水面,能夠獲取前進推進力、平衡體重的托舉力和維持奔跑姿態的轉矩,為推進器設計提供了靈感[9-10]。以蛇怪蜥蜴為仿生對象,徐林森等[11-16]以平板型腳掌、橢圓形軌跡和5 Hz以上高頻率踏水模擬設計,獲取了推進力、托舉力和轉矩三維驅動,初步實現了水面平臺0.5 m/s的低速運動,正在開展結構穩定性研究;而傳統的明輪推進器同樣基于平板型葉片結構,以葉片平面徑向布置為主,通過葉片圓形軌跡的旋轉運動,8~10個葉片、1 r/s左右的低轉速設置,提供垂直于輪軸的推進力,但是推進效率還有待提高[17]。結合以上高頻靈活、三維驅動輸出的平板設計和穩定旋轉運動兩個特性,課題組前期已提出了仿生葉輪推進器,并開展了初步的葉片結構設計和分析[18-19]。

本文進一步針對蛇怪蜥蜴踏水過程中腳掌的轉動變化,提出一種變角度的水面矢量推進器;通過平板與輪輻的夾角變化實現驅動力矢量輸出,結合轉速控制提高驅動調節能力。為獲取推進器的三維驅動輸出性能,建立流體動力學模型,研究運動參數和結構參數對托舉力、推進力和轉矩的影響;進行推進器敞水試驗,結合高速流場分析驗證數值計算模型的正確性。

1 結構仿生設計

基于蛇怪蜥蜴的踏水運動,前期整體結構固連的仿生葉輪推進器初步增加了托舉力和轉矩輸出。然而,在蛇怪蜥蜴踏水時,腳掌平面相對小腿的角度在時刻變化,以此調整托舉力、推進力和轉矩的輸出,滿足托起體重、推動身體和維持前傾奔跑姿態的需求變化[9],其運動姿態如圖1a所示。

圖1 蛇怪蜥蜴踏水運動姿態變化和腳掌結構圖Fig.1 Kinematic phases of a stride and foot structure of basilisk lizard

不同質量的蛇怪蜥蜴腳掌平面(含腳趾)均近似長方形,如圖1b所示。基于腳掌與水流的相互作用,通過姿態角的變化,蛇怪蜥蜴實現了踏水驅動力的矢量輸出,其中后腿關節處肢體角度一個運動周期內的變化如圖2所示。

圖2 腿部關節角度Fig.2 Angular excursions of limb joints

模仿蛇怪蜥蜴腳踝處腳掌繞小腿的轉動行為,以旋轉運動代替往復踏水,本文提出了一種變角度水面矢量推進器。推進器由葉片、伸縮桿、輪輻和輪轂組成,如圖3所示。其中,近似腳掌的長方形葉片與輪輻末端鉸接,伸縮桿連接輪輻和葉片,通過長短伸縮調整葉片與輪輻的夾角θ,如同腳掌與小腿夾角θ1的改變,控制葉片的拍水角度。

圖3 變角度水面矢量推進器Fig.3 Water-surface vector propeller of variable angle1.葉片(長度0.06 m,寬度0.05 m) 2.伸縮桿 3.輪輻(含輪轂,長度0.05 m) 4.輪轂

結合葉片夾角和轉速參數調節,實現仿生推進器托舉力Fl、推進力Ft和轉矩M的高效輸出,對應方向如圖3所示。因推進器傳動難以實現高頻轉動過程中夾角θ的實時快速變化,所以本文更側重角度調整穩定后的持續轉動性能。

為驗證推進器的三維驅動輸出特性,開展了2對推進器的反向轉動試驗,尺度近似比例增大的原理樣機采用0.3 m×0.2 m葉片,轉速均為1.5 r/s。由圖4靜止和轉動狀態對比可看出,樣機被托起了0.2 m,驗證了托舉力的存在;中部0.03 m×0.03 m截面的2.5 mm厚不銹鋼管向內側下方折彎,驗證了轉矩的存在;因2對推進器的推進力相互抵消,沒有水平位移。

圖4 矢量推進器反向轉動試驗Fig.4 Reverse rotation tests of vector propeller

伴隨著氣穴生成與潰散的流場變化,通過腳掌轉動和踏水頻率的調控,蛇怪蜥蜴踏水高效地獲得了三維驅動力,實現了敵害躲避。而要進一步研究矢量推進器的驅動輸出性能,也應結合流場變化,對應不同的參量,開展流體動力學分析。

2 數值計算模型

首先,通過理論分析獲取影響推進器三維驅動輸出的參數變量;然后,建立流體動力學模型,開展詳細的數值計算。

2.1 理論分析

葉片平板旋轉拍擊水面的固液相互作用可看作水流沖擊固定葉片的平板繞流。位于圖5所示位置時,在葉片微元ds上,產生正向力dD、法向力dT和繞轉軸的轉矩dM[20]

dD=0.5CDρU2ds

(1)

dT=0.5CTρU2ds

(2)

dM=MdT+MdD

(3)

式中CD、CT——正向力和法向力系數

ρ——流體密度U——來流速度

ds——微元面積

MdT、MdD——正向力和法向力對輪軸的轉矩

圖5 葉片受力分析Fig.5 Force analysis model of blade

經坐標變化,微元ds受到的作用進一步表示為推進力dFax、托舉力dFay和轉矩dM

dFax=-dTsinα+dDcosα

(4)

dFay=dTcosα+dDsinα

(5)

dM=dFaxksinβ+dFaykcosβ=

dTcos(α+β)+dDsin(α+β)

(6)

式中α——正向力與水平軸的夾角

β——微元與輪軸連線的轉角

k——微元到轉軸的力臂

推進器的三維驅動輸出可由4個葉片上的微元積分求和獲得。由式(4)~(6)可得:三維驅動輸出是U、β的函數,而U、β又分別與推進器轉速、葉片輪輻夾角直接相關,同時結構參數也是重要的影響因素。

2.2 流體動力學模型

建立流體動力學模型,分析運動參數(轉速w和葉片夾角θ)和結構參數(輪輻長度l、輪轂寬度b和輪轂直徑d)對三維驅動輸出的影響,如圖6所示。

圖6 結構參數示意圖Fig.6 Schematic of structure parameters

基于k-εRNG湍流模型,應用PISO算法和流體體積法跟蹤葉片水面轉動的流固耦合作用。在不可壓非定常流動假設的基礎上,任一單元的流體體積分數方程、密度方程、粘度方程、連續方程、動量方程為

αw+αa=1

(7)

ρ=αwρw+αaρa

(8)

μ=αwμw+αaμa

(9)

(10)

(11)

式中αw、αa——水和空氣在單元中所占體積比

ρw、μw、ρa、μa——水、空氣密度和粘度

ui、uj——速度分量均值(i、j=1,2,3)

p——壓力均值

μ——流體粘性系數

Si——廣義源項

推進器輪軸位于液面,流體計算域尺寸為3 m×1 m×1 m。采用向邊界等比放大的方式優化網格,推進器表面網格尺寸為0.005 m,如圖7所示。計算域上側為空氣,下側為水。

3 性能分析

基于流體動力學模型,數值計算了運動參數和結構參數對三維驅動輸出的影響。基礎參數:葉片夾角60°;轉速1 r/s;輪輻長度0.05 m(含輪轂);輪轂寬度0.02 m;輪轂直徑0.04 m;葉片尺寸0.06 m×0.05 m×0.005 m。

3.1 葉片夾角

設定葉片夾角在-90°~90°之間變化(圖6中夾角θ水面上取正值、水面下取負值),其中6個不同角度的三維驅動輸出曲線如圖8所示。從圖8可看出,驅動輸出呈周期性變化,頻率為4 Hz,這是因為4個葉片依次入水產生了周期性的流固耦合作用;驅動幅值隨著角度的變化先增大后減小,推進力、托舉力和轉矩分別在0°、-20°和0°取得最大幅值。

最大幅值出現在0°附近是因為此時葉片與旋轉方向呈約90°,葉片與水流相互作用劇烈,隨著角度增大或減小,都對流固耦合作用形成了緩沖。相互作用強度可由圖9流場壓力反映,左下側葉片處流場壓力0°時大部分超過1 200 Pa,而-50°和50°時經過緩沖大部分在1 000 Pa左右。

圖8 不同角度驅動輸出曲線Fig.8 Force characteristic curves at different angles

圖9 流場壓力Fig.9 Fluid pressure at different angles

以1 s時間內驅動力和轉矩的平均值作為特征值進行推進器驅動輸出性能分析。三維驅動隨角度的變化如圖10所示,隨角度增大托舉力先減小后增大為正值再減小;推進力和轉矩先增大后減小,在0°夾角時取得最大值。

圖10 驅動性能隨葉片夾角變化曲線Fig.10 Curves of driving properties changing with blade angle

由圖10同時可看出,θ<0°時,托舉力為負值,推進器下壓平臺,會增大浸濕面積而增大阻力;0°<θ<40°時,托舉力由負值逐漸向正值增加,但數值較小;θ>60°時,推進力和轉矩迅速下降,數值較小。綜合托舉力、推進力和轉矩輸出變化,得:40°≤θ≤60°時,三維驅動輸出均衡,沒有過小項,所以常態航行時夾角應設置在此區間。

進一步得到了驅動力的矢量輸出角,即托舉力和推進力合力的夾角,如圖11所示。由圖11可看出,隨葉片角度從0°至90°增大,驅動力輸出角從-6.85°連續增加至51.95°,實現了驅動力輸出的矢量變化,可根據航態匹配驅動力需求。

圖13 流場壓力隨輪輻長度變化圖Fig.13 Changing diagrams of fluid pressure with different spoke lengths

圖11 驅動力輸出角隨葉片夾角變化曲線Fig.11 Curve of driving force angle changing with blade angle

3.2 輪輻

輪輻長度l(含輪轂部分)從0.02 m依次增加0.01 m至0.10 m,得到了9組三維驅動輸出特征值(li,Fi),i=0,1,…,8,如圖12所示。由圖12可看出,三維驅動隨著輪輻長度的增加而增加。這是因為,同轉速下,葉片上任一點的速度隨著輪輻長度增加而提高,因而與水流的相互作用加強,如圖13中壓力場所示。由圖13可看出,輪輻長度從0.02 m增加到0.10 m,前側入水葉片壓力由750 Pa提升至1 250 Pa,后側葉片壓力由部分750 Pa提升至全部1 000 Pa。

圖12 驅動性能隨輪輻長度變化曲線Fig.12 Curves of driving properties changing with spoke length

(12)

3.3 轉速

轉速的提高增加了固液作用的頻率,其中3、5、7 r/s的力學特性曲線如圖14所示。可以看出,3 r/s和7 r/s時驅動周期性輸出的頻率均為轉速的4倍,但是5 r/s轉速時輸出的頻率為轉速的2倍。

圖14 不同轉速下的力學特性曲線Fig.14 Force characteristic curves at different rotation velocities

頻率變化是因為在低轉速下,流場在下一葉片拍擊之前即可恢復,隨著轉速提高至5 r/s左右,固液作用頻率加快,而下一葉片到來時流場來不及恢復,2個葉片作用過程是一個流場恢復周期;如圖15a、15b中前側入水葉片相互作用的流體較多,而圖15c、15d中下一葉片入水相互作用的流體較少。當轉速進一步加快,葉片與水流相互作用的位置向下移動,靜壓力增加加快了流場的恢復速度,輸出的頻率又變回為轉速的4倍。

轉速由1 r/s向7 r/s增加,三維驅動輸出性能如圖16所示。由圖16可看出,托舉力、推進力和轉矩分別由1 r/s時的0.50 N、0.59 N、0.083 N·m,增加至7 r/s時的6.09 N、11.61 N、1.17 N·m;三維驅動輸出隨著轉速的增加而快速增加,驅動輸出存在較大的調節范圍。

圖15 5 r/s轉速下兩相流分布圖Fig.15 Water volume color contours of propeller when w=5 r/s

圖16 驅動性能隨轉速變化曲線Fig.16 Curves of driving properties changing with rotation velocity

3.4 輪轂

寬度和直徑是圓柱形輪轂的2個結構參數。當寬度由0.01 m向0.05 m增加時,托舉力減小、而推進力和轉矩在0.60 N和0.083 N·m左右浮動,如圖17所示。以0.02 m寬和0.04 m直徑參數為基準,推進力和轉矩的數值浮動均小于6%,影響不大。

圖17 驅動性能隨輪轂寬度變化曲線Fig.17 Curves of driving properties changing with width of hub

當直徑由0.02 m增加至0.08 m時,托舉力隨著直徑的增加而減小,如圖18所示。推進力和轉矩小幅浮動,分別較直徑為0.04 m時的數值最大變化了7.8%和5.4%。

圖18 驅動性能隨輪轂直徑變化曲線Fig.18 Curves of properties changing with diameter of hub

4 推進器試驗

為了驗證數值計算的正確性,搭建了推進器試驗系統。試驗系統由推進器、水平導軌與傳感器、豎直導軌與傳感器和傳動系統組成,如圖19所示,結構參數同數值計算模型。

圖19 矢量推進器試驗系統Fig.19 Test system of vector propeller1.導軌 2.推進器 3.水平力傳感器 4.滑塊 5.傳動系統 6.豎直力傳感器

測得2.5 r/s轉速下的托舉力和推進力曲線,與數值計算結果對比,如圖20所示。由圖20可看出,試驗曲線與數值計算結果作用力變化頻率均為10 Hz,是轉速的4倍。托舉力和推進力試驗平均值分別為2.13 N和1.95 N;數值計算托舉力為2.41 N,較試驗值高13.2%;數值計算推進力為1.60 N,較試驗值低17.9%,誤差在合理的范圍內。

圖20 力學特性對比曲線Fig.20 Contrasts of force characteristic

數值計算與試驗曲線中,托舉力的波峰和推進力的波谷位置一致。但是試驗曲線托舉力峰值小、凸起時間長,推進力峰值大但凸起時間短,作用力曲線未完全對應是因為:在試驗中采用電動機驅動,不能保證恒角速度轉動,在一個作用周期內角速度有所起伏。

進一步測試了不同轉速條件下的托舉力和推進力變化,如圖21所示,可得:驅動力輸出隨著轉速的增加而增加,數值計算值與試驗值符合較好。

圖21 驅動力隨轉速的變化曲線Fig.21 Curves of driving forces changing with rotation velocity

最后,由高速攝像機獲取了轉速為5 r/s時的推進器流場,如圖22所示。流場中黑色部分為輪輻和葉片,葉片剛剛經過的區域產生了負體積的氣穴,如圖中三角形所在區域;葉片與流體的相互作用,一次大體積水流相互作用后緊接著一次小體積水流相互作用,2個葉片拍水過程是一個流場恢復周期,相互作用的流體體積如圖22中環形所示,與圖15所示數值計算流場結果一致。

圖22 推進器試驗流場Fig.22 Fluid fields of propeller test

5 結束語

以蛇怪蜥蜴為仿生對象,設計了變角度仿生矢量推進器,通過角度和轉速調節,為水面平臺提供三維驅動矢量輸出,獲得了減阻提速的新思路。進一步建立流體動力學模型,擴展了平板與水流相互作用的數值計算研究,得到了驅動力隨葉片角度調節矢量輸出的角度變化范圍;推導了托舉力、推進力和轉矩與輪輻長度的函數關系;明確了轉速和輪轂結構參數對三維驅動輸出的影響規律。最后通過試驗驗證了數值計算模型的正確性。

1 薛軍. 水陸兩棲車發展史[J]. 上海汽車,2001(7):46.

XUE Jun. Development history of amphibious vehicle[J].Shanghai Auto,2001(7):46. (in Chinese)

2 楊楚泉. 水陸兩棲車輛原理與設計[M]. 北京:國防工業出版社,2003.

3 辛志坡,王偉. 高速水陸兩棲車輛技術發展[J]. 專用汽車,2007(6):18-21.

XIN Zhipo, WANG Wei. Technology development of high-speed amphibious vehicle[J]. Special Purpose Vehicle, 2007(6):18-21. (in Chinese)

4 鄧銳, 黃德波, 周廣利, 等. 帶有阻流板的雙體船水動力性能數值研究[J]. 華中科技大學學報:自然科學版,2011,39(4):97-100.

DENG Rui, HUANG Debo, ZHOU Guangli, et al. Numerical research on hydrodynamics performances of catamarans with interceptors[J]. Journal of Huazhong University of Science and Technology:Natural Science Edition,2011,39(4):97-100. (in Chinese)

5 唐建飛, 楊帥. 高速艇水翼減阻方案及翼滑艇阻力估算方法[J]. 船舶力學,2015,19(1-2):69-77.

TANG Jianfei, YANG Shuai. Drag reduction regime by hydrofoil and resistance estimation method for gliding-hydrofoil craft[J]. Journal of Ship Mechanics, 2015,19(1-2):69-77. (in Chinese)

6 鐘本基, 倪其軍, 劉玉川. 動力增升降阻技術在快速船上的應用[J]. 中國造船,2008,49(1):11-18.

ZHONG Benji, NI Qijun, LIU Yuchuan. Technology on dynamic enhancing lift & dropping resistance and its using in high speed craft[J]. Shipbuilding of China, 2008,49(1):11-18. (in Chinese)

7 孫華偉. 滑行面形狀對滑行艇阻力與航態影響數值分析[D]. 哈爾濱:哈爾濱工程大學,2012.

SUN Huawei. Numerical analysis of planning-hull surface shape on resistance and sail attitude[D]. Harbin:Harbin Engineering University,2012. (in Chinese)

8 劇冬梅, 項昌樂, 周鵬飛, 等. 縱傾角對輕型輪式兩棲車輛的阻力特性影響研究[J]. 兵工學報,2015,36(1):19-26.

JU Dongmei, XIANG Changle, ZHOU Pengfei, et al. Analysis of the effect of trim angle on the resistance characteristics for wheeled amphibious vehicle[J]. Acta Armamentatii,2015,36(1):19-26. (in Chinese)

9 HSIEH S T, LAUDER G V. Running on water: three-dimensional force generation by basilisk lizards[J]. PNAS,2004, 101(48):16784-16788.

10 GLASHEEN J, MCMAHON T. Size-dependence of water-running ability in basilisk lizard[J]. The Journal of Experimental Biology,1996, 99(12):2611-2618.

11 徐林森, 魏鮮明, 曹凱. 仿生雙足水上行走機器人優化設計及控制方法[J]. 機械工程學報,2014,50(15):12-18.

XU Linsen, WEI Xianming, CAO Kai. Optimization design and control method of biped robot running on water[J]. Journal of Mechanical Engineering, 2014,50(15):12-18. (in Chinese)

12 XU Linsen, MEI Tao, WEI Xianming, et al. A bio-inspired biped water running robot incorporating the Watt-I planar linkage mechanism[J]. Journal of Bionic Engineering,2013,10(4):415-422.

13 PARK H, FLOYD S, SITTI M. Roll and pitch motion analysis of a biologically inspired quadruped water runner robot[J]. The International Journal of Robotics Research,2010,29(10):1281-1297.

14 FLOYD S, ADILAK S, RAMIREZ S, et al. Performance of different foot designs for a water running robot[C]∥IEEE International Conference on Robotics and Automation, 2008:19-23.

15 KIM H, JEONG K, SEO T. Analysis and experiment on the steering control of a water-running robot using hydrodynamic forces[J]. Journal of Bionic Engineering,2017,14(1):35-46.

16 KIM H, SITTI M, SEO T. Tail-assisted mobility and stability enhancement in yaw and pitch motions of a water-running robot[J]. IEEE/ASME Transactions on Mechatronics,2017,22(3):1-11.

17 黃渙青. 內河小型清潔船明輪推進研究[D]. 鎮江:江蘇科技大學,2013.

HUANG Huanqing. Study on paddle wheel propulsion of the water clean boat[D]. Zhenjiang:Jiangsu University of Science and Technology,2013. (in Chinese)

18 白向華,呂建剛,郭劭琰,等.基于蛇怪蜥蜴踏水機理的仿生葉輪負載性能分析[J].華中科技大學學報:自然科學版,2015,43(2):76-80.

BAI Xianghua, Lü Jian’gang, GUO Shaoyan, et al. Performance analysis of bionic impeller under different loads based on mechanism of basilisk lizard treading water[J]. Journal of Huazhong University of Science and Technology:Natural Science Edition, 2015, 43(2):76-80. (in Chinese)

19 張仲志, 呂建剛, 宋彬, 等. 水面矢量推進器仿生設計與試驗[J/OL] 農業機械學報,2016,47(10):405-412.http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?file_no=20161053&flag=1.DOI:10.6041/j.issn.1000-1298.2016.10.053.

ZHANG Zhongzhi, Lü Jian’gang, SONG Bin, et al. Bionic design and experiment of water-surface vector propeller[J/OL]. Transactions of the Chinese Society for Agricultural Machinery,2016,47(10):405-412. (in Chinese)

20 CHRISTINA G, MEYER N, MARTIN B. Simulation of an underwater hexapod robot[J]. Ocean Engineering, 2009,36(1):39-47.

主站蜘蛛池模板: 免费一级无码在线网站| 老司国产精品视频91| 黄色一及毛片| 青青草a国产免费观看| 无码精品福利一区二区三区| 58av国产精品| 四虎成人精品| 日韩欧美中文字幕在线韩免费 | 欧美激情视频一区| 日韩小视频在线观看| 不卡网亚洲无码| 国产精品久久久精品三级| 美女无遮挡免费网站| 2022国产91精品久久久久久| AV不卡国产在线观看| 欧美日本中文| 曰AV在线无码| 国产成人精品高清不卡在线| 四虎永久在线精品国产免费| 国产一区二区三区免费| 亚洲精品视频在线观看视频| 小蝌蚪亚洲精品国产| 久热中文字幕在线| 色婷婷成人网| 一级毛片不卡片免费观看| 97无码免费人妻超级碰碰碰| 亚洲男人的天堂久久香蕉| 男人天堂伊人网| 亚洲综合久久成人AV| аv天堂最新中文在线| 青草视频在线观看国产| 国产经典免费播放视频| 国产精品免费福利久久播放| 一区二区欧美日韩高清免费| 中文字幕乱码中文乱码51精品| 欧美日韩国产系列在线观看| 综合色区亚洲熟妇在线| 日本一本正道综合久久dvd| 国产在线一二三区| 伊人欧美在线| 成人午夜免费视频| 青草精品视频| 野花国产精品入口| 香蕉国产精品视频| www中文字幕在线观看| 亚洲成人播放| 在线一级毛片| 国产美女自慰在线观看| 精品久久综合1区2区3区激情| 青青草久久伊人| 免费观看国产小粉嫩喷水 | 99ri国产在线| 国产成人一区| 久久久久久久蜜桃| 91精品国产自产91精品资源| m男亚洲一区中文字幕| 国产69精品久久| 国产三级成人| 午夜啪啪网| 婷婷色中文网| 一本一道波多野结衣一区二区| 国产无套粉嫩白浆| 国产欧美日韩91| 人妻无码中文字幕第一区| 在线色综合| 国产成人午夜福利免费无码r| 国产日本视频91| 成人久久18免费网站| 日韩精品无码免费专网站| 一区二区自拍| 国产成人综合日韩精品无码首页| 中文字幕伦视频| 精品国产免费人成在线观看| 黄色网站在线观看无码| 国产日韩精品欧美一区喷| 国产91在线|日本| 永久在线播放| 国产成人三级| 国产人人乐人人爱| 69国产精品视频免费| 暴力调教一区二区三区| 国产精品hd在线播放|