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

半轉翼懸停和前進飛行升力估算方法

2017-08-31 13:18:42王孝義張玉華董銀萍陳富強邱支振
中國機械工程 2017年15期

王孝義 張玉華 董銀萍 邱 晗 陳富強 邱支振

安徽工業大學機械工程學院,馬鞍山,243002

半轉翼懸停和前進飛行升力估算方法

王孝義 張玉華 董銀萍 邱 晗 陳富強 邱支振

安徽工業大學機械工程學院,馬鞍山,243002

在分析半轉翼運動模型和翼面氣流特點的基礎上,建立了懸停和前進兩種飛行狀態下的半轉翼升力計算模型。根據半轉翼的運動特性,推導出適合半轉翼運動的升力計算解析表達式。結合半轉翼樣機參數,應用導出的理論公式和基于CFD軟件的數值仿真模型,分別計算不同飛行條件下半轉翼的升力,獲得半轉翼懸停和前進兩種飛行狀態下升力變化規律。理論計算與數值仿真所得升力曲線的對比驗證了升力估算解析法的有效性和可行性。研究結果可為半轉翼飛行器的參數設計與升力預估提供理論指導。

半轉翼;懸停飛行;前進飛行;升力計算

0 引言

動物肢體運動形式雖各有不同,但本質上都是“不對稱擺動”[1]。擺動是適應肌肉特點的運動形式,動物利用這種不對稱運動產生了高效而巧妙的運動效果,像鳥類和昆蟲的飛行,機動靈活,非常適合復雜環境下的活動。撲翼飛行器是模仿鳥類或昆蟲飛行方式的飛行器,因其在國防和民用領域的廣闊應用前景而成為各國學者研究的熱點[2-5]。特別是在撲翼氣動力研究方面,人們開展了理論建模[6-7]、不同工況升力試驗[8-12]以及基于CFD理論的數值仿真[13-14]等一系列的研究,取得了較為豐富的成果[15]。然而,撲翼的往復快速擺動會產生較大慣性力[16-17],導致撲翼飛行器在大型化方面遇到了難以逾越的障礙[18]。在此背景下,安徽工業大學提出了一種能避免擺動而只產生不對稱運動的轉動機構——半轉機構[1,19],通過對半轉機構的不斷改進,研制出一種簡約化半轉機構作為類撲翼飛行器驅動機構。這種新型飛行驅動機構具有結構簡單、飛行驅動效率高、適應能力強等優點,可嘗試在大尺寸、高升力上取得突破,從而使大尺寸半轉翼飛行器成為可能。

在半轉翼飛行器的研制中,飛行器的升力預估是最基本的環節。半轉翼飛行器作為一種以轉動替代擺動的新型類撲翼系統,針對其升力特性方面的研究尚未展開。本文通過分析基于簡約化半轉機構的半轉翼運動模型和翼面氣流特點,對半轉翼飛行器懸停和前進兩種飛行狀態下的升力形成進行理論研究,旨在建立適合半轉翼運動的氣動力模型,探索半轉翼飛行升力估算方法。

1 半轉翼運動及翼面氣流

1.1 半轉翼運動模型

半轉翼是一種適應于飛行的簡約化半轉機構的輸出構件,也是半轉翼飛行器產生升力的工作構件。如圖1所示,AB為半轉翼,翼展長度為2a;O1C為曲柄,半徑為R。曲柄與半轉翼的鉸接中心C位于半轉翼的中部;定點O(又稱不動點)處安裝轉動滑塊,使得定點O始終位于半轉翼AB上。

圖1 半轉翼的平面運動模型Fig.1 Planar motion model of HRW

因以上幾何約束作用,半轉翼在曲柄的帶動下做平面運動,即繞定點O的轉動以及隨鉸接中心C點的平移。在定點O建立固定坐標系OXY,X軸通過曲柄轉動中心O1,在動點C建立動坐標系Cxy,x軸方向沿翼的展向。設ω為曲柄的轉動角速度,則半轉翼轉動角速度為ω/2,半轉翼上任意位置D的速度為

(1)

式中,φ為曲柄的轉角;x為D位置處沿翼展向的位置坐標。

令x=0,得C點的速度

vC=ωR

(2)

由于C點位于半轉翼AB的中點,為防止B點從O點脫出, 必有a>2R。如圖1所示,當半轉翼逆時針轉動到不動點O的左側時,半轉翼逐漸向下運動,相當于昆蟲翅的下拍行程;當半轉翼逆時針轉動到不動點O的右側時,半轉翼逐漸向上運動,相當于昆蟲翅的上揮行程。據此,可將半轉翼運動區間分為下拍區和上揮區[1]。由圖1可知,下拍區和上揮區區間分別為-2Rsin(φ/2)

1.2 半轉翼翼面氣流與迎角

半轉翼飛行器有懸停和前進兩種飛行狀態,兩種飛行狀態下的翼面氣流組成不同。考慮到半轉翼翼片運動的周期性,對曲柄轉動一周的過程進行分析。

當半轉翼飛行器處于懸停飛行狀態時,半轉翼的運動是平面運動,因其在運轉中有明顯的展向伸縮運動,故翼片上會產生沿翼面展向分布的展向氣流,如圖2所示??紤]到除翼片兩頂端的展向氣流外,翼面中間的展向氣流不能越過翼片,因而并不產生繞流作用,也就是說展向氣流因繞流產生的升力甚小,因此,在翼面氣動力元分析計算中可忽略展向氣流的作用,只考慮在翼面弦向截面內的氣流作用。

圖2 懸停時半轉翼翼面氣流Fig.2 Airflow distribution on the surface of HRW in hovering flight

考慮翼面與靜止空氣的相對速度,氣流的速度沿半轉翼弦向的分布是均勻的,而沿展向的分布是變化的。半轉翼上任意位置D的氣流速度v可用式(1)計算,方向與vD相反,氣流與半轉翼的迎角α為vD與翼面切向夾角之銳角,規定下拍區α為正,上揮區α為負。由圖1可得,懸停狀態下氣流迎角

(3)

當半轉翼飛行器以速度v0前進飛行時(圖3),相當于飛行器固定而前方有-Z方向的來流,速度大小為v0,來流在半轉翼上產生弦向氣流。此時翼片上既有展向氣流也有弦向氣流,將半轉翼任意點的展向氣流速度分解為平行于翼面的分量vx和垂直于翼面的分量vy。如前所述,平行于翼面的展向氣流分量對飛行升力的貢獻甚微,可忽略不計(故圖3中只給出弦向氣流示意);前方來流產生的弦向繞流可產生較大升力,此外在弦向截面內還應考慮vy對氣流相對速度的影響,于是可得氣流相對速度

圖3 前進時半轉翼翼面氣流Fig.3 Airflow distribution on the surface of HRW in forward flight

(4)

vy=ω(Rsin(φ/2)+x/2)

前進狀態下氣流迎角用γ表示, 則由圖3中v0和vy關系可得

(5)

由式(5)可見,氣流迎角與作用點在弦向的位置無關,而與來流速度v0、曲柄轉速ω、轉角φ及作用點在展向的位置x有關。令x=0,C點在前進狀態下迎角

(6)

因sin(φ/2)∈[0,1],故C點迎角最大值

γCmax=arctan(ωR/v0)

(7)

2 半轉翼懸停飛行升力估算

2.1 懸停飛行升力元

在靜止空氣中懸停飛行時,氣流相對翼片的速度只由翼片上相應點的速度決定。沿展向的位置不同,氣流的速度和迎角都是變化的,半轉翼上各點處的氣動力也不同。

沿展向取氣動力元(圖4),dD是對應迎角α的阻流力元,它們在固定坐標系Y方向上的投影構成阻流升力元;dL是對應迎角α的繞流力元。懸停時繞流只能沿展向發生,可忽略展向氣流的繞流作用,故懸停時只需計算阻流升力。由于懸停飛行時絕大多數時間處于大迎角狀態,故采用平板大攻角氣動力的近似計算公式[20]來估算其阻流升力。需要指出的是,半轉翼在下拍區產生向上的升力,在上揮區產生向下的負“升力”。

圖4 懸停時半轉翼氣動力元Fig.4 Aerodynamic element of HRW in hovering flight

2.2 懸停飛行的升力

懸停飛行時的升力按阻流升力計算,阻流力元計算式為

dD=0.5CD1ρhv2sin2αdx

(8)

式中,ρ為空氣的密度;CD1為懸停飛行時的壓差阻力系數,由實驗測定。

進一步可得下拍區和上揮區的阻流升力元統一表達式:

dFLB=-dDcos(φ/2+α)

(9)

對阻流升力元沿展向在全翼進行積分,得全翼升力FLB(即懸停飛行總升力FL):

(10)

3 半轉翼前進飛行升力估算

由1.2節分析可知,前進飛行狀態下半轉翼飛行升力主要考慮弦向氣流的影響。氣流的速度和迎角沿弦向是固定不變的(不計翼邊界效應),但兩者隨來流速度v0、曲柄轉速ω、曲柄轉角φ以及作用點沿展向的位置x的變化而變化。為區分小迎角和大迎角工況,以半轉翼展向對稱中點C的迎角γC作為判定量。令進入小迎角區的臨界迎角為γ0,則|γC|<γ0為小迎角狀態,而|γC|≥γ0為大迎角狀態。

3.1 半轉翼前進飛行升力元

如圖5所示,沿半轉翼的展向截取弦向氣流氣動力元,氣動力元的前緣為半轉翼的前緣。半轉翼沿展向任意微段dx受到繞流力元dL和阻流力元dD,dF代表dL和dD在垂直翼面方向分力的合力。飛行升力元dFL為dF在鉛垂Y方向的分力。迎角γ規定逆時針為正,順時針為負,由圖5可知下拍區γ為正,上揮區γ為負。

圖5 前進時半轉翼弦向氣動力元Fig.5 Aerodynamic element of chordwise HRW in forward flight

3.2 小迎角條件下弦向繞流升力

小迎角時,繞流升力由繞翼環量確定。圖5中,半轉翼弦向氣流氣動力元上的環量

(11)

任意展向長度為dx的翼元上由弦向氣流產生的繞流力元

dLS=ρvΓdx=πρv2hsinγdx

(12)

進一步可得下拍區和上揮區的繞流升力元統一表達式:

(13)

對繞流升力元在全翼進行積分,可得小迎角條件下(|γC|<γ0)全翼弦向繞流升力

(14)

3.3 大迎角條件下弦向繞流升力

對半轉翼大迎角條件下弦向氣流產生升力的分析,仍采用平板大攻角繞流氣動力的近似計算公式估算繞流升力與阻流升力。

大迎角(|γC|≥γ0)時,任意展向長度為dx的翼元上由弦向氣流產生的繞流力元

dLL=0.5CD2ρv2hsinγcosγdx

(15)

其中,CD2為前進飛行時的壓差阻力系數,一般也通過實驗測定。需要指出的是,前進和懸停飛行條件下壓差阻力系數并不相同。

進一步可得下拍區和上揮區弦向氣流產生的繞流升力元統一表達式:

(16)

對繞流升力元在全翼進行積分,可得大迎角條件下全翼弦向繞流升力

(17)

3.4 半轉翼前進飛行阻流升力

半轉翼前進時的阻流升力主要出現在大迎角范圍內。大迎角時,任意展向長度為dx的翼元上由弦向氣流產生的氣動阻流力元

dD=0.5CD2ρv2hsin2γdx

(18)

進一步可得下拍區和上揮區弦向氣流產生的阻流升力元統一表達式:

(19)

對阻流升力元在全翼進行積分,可得大迎角條件下全翼弦向阻流升力:

(20)

3.5 前進飛行總升力

半轉翼前進飛行總升力FL為弦向繞流升力FLA和弦向阻流升力FLB的合力。在一個運動周期內,FL根據迎角不同而分段計算,即

(21)

其中,FLAS、FLAL、FLB分別用式(14)、式(17)和式(20)計算。

4 半轉翼升力估算方法的數值驗證

為驗證半轉翼升力估算理論模型,根據已研制的半轉翼樣機(圖6)模型參數,采用前述升力公式對不同飛行條件下的升力進行預估計算,并將理論計算結果與相同條件下利用XFLOW軟件數值計算結果進行對比。

圖6 半轉翼樣機(對稱布置雙翼片)Fig.6 HRW prototype with double symmetrical wings

4.1 升力估算實例

半轉翼樣機主要尺寸參數:曲柄半徑R=0.06 m,翼片展向長度2a=0.282 m,弦向長度h=0.2 m??諝饷芏圈?1.225 kg/m3。懸停飛行狀態下壓差阻力系數CD1由實驗測得(其值為4.1)。前進飛行時取γ0=π/9,壓差阻力系數CD2采用Dickinson等[12]通過實驗修正得到的經驗值3.46。設置四種懸停飛行條件(曲柄轉速ω為6π rad/s,10π rad/s,20π rad/s,30π rad/s)和六種前進飛行條件(ω=20π rad/s,v0=2 m/s;ω=30π rad/s,v0=2 m/s;ω=20π rad/s,v0=3 m/s;ω=30π rad/s,v0=3 m/s;ω=20π rad/s,v0=6 m/s;ω=30π rad/s,v0=6 m/s),分別計算上述飛行條件下曲柄運轉一周內不同位置對應的半轉翼懸停和前進飛行升力理論值。另外,根據半轉翼樣機的尺寸參數,在流體動力學仿真軟件XFLOW中建立數值仿真模型,設定相關參數,分別針對以上飛行條件,對半轉翼懸停和前進飛行的升力進行數值計算,獲得升力仿真值。

4.2 理論計算結果與仿真結果的比較分析

圖7所示為半轉翼懸停飛行理論升力和仿真升力曲線比較,同時反映出不同曲柄轉速下的升力變化規律。隨著曲柄轉速的增大,升力增大,且升力增大倍數近似等于轉速增大倍數的平方。圖8所示為半轉翼前進飛行時的理論升力和仿真升力曲線比較,同時也反映出不同前進飛行條件下升力變化規律。來流速度相同時,升力隨曲柄轉速增大而增大;曲柄轉速相同時,來流速度增大,升力變化卻較小,說明平行于翼面的氣流對半轉翼升力影響不大。

(a)ω=6π rad/s

(b)ω=10π rad/s

(c)ω=20π rad/s

(d)ω=30π rad/s圖7 半轉翼懸停飛行理論升力和仿真升力曲線Fig.7 Lift curves of HRW in hovering flight by theoretical calculation and numerical simulation

綜合圖7和圖8可見:①理論計算升力和仿真升力的變化趨勢基本一致,兩者誤差小于20%,說明半轉翼飛行升力的理論估算方法是有效的。②在曲柄行至中間位置,即半轉翼水平下拍瞬間,半轉翼飛行升力達到最大值,情況與實際相符。③懸停飛行時,半轉翼升力只與曲柄轉速有關;前進飛行時,半轉翼升力主要取決于曲柄轉速,來流速度對升力的影響比曲柄轉速的影響小得多。

導致理論計算升力和仿真升力存在誤差的主要原因是理論與仿真的計算條件不同,理論估算中假定翼片在理想的靜止流場中運動,而仿真計算中模仿真實流場考慮了翼片對流場的擾動。因此,圖7中仿真值小于理想狀態的理論值,而且轉速越高,翼片對流場的擾動越強烈,導致理論值與仿真值的偏差越大。另外,由于翼片的運動平面與來流方向垂直,來流會把翼片運動區域內的擾動帶向下游,減少了翼片運動流場中的擾動程度,從而提高了升力數值;而且來流速度越大,帶走擾動的能力越強,所以圖8中當轉速相同時,隨著來流速度提高,仿真值與理論值的誤差變小。

(a)ω=20π rad/s,v0=2 m/s (b)ω=30π rad/s,v0=2 m/s (c)ω=20π rad/s,v0=3 m/s

(d)ω=30π rad/s,v0=3 m/s (e)ω=20π rad/s,v0=6 m/s (f)ω=30π rad/s,v0=6 m/s圖8 半轉翼前進飛行理論升力和仿真升力曲線Fig.8 Lift curves of HRW in forward flight by theoretical calculation and numerical simulation

5 結論

(1)半轉翼懸停飛行時,翼面氣動力主要為由法向氣流產生的阻流升力,建立了適合半轉翼運動的阻流升力計算模型,導出了升力計算公式,可以估算半轉翼懸停飛行狀態下的升力大小。

(2)半轉翼前進飛行時,升力主要來源于弦向氣流。建立了適合半轉翼運動的弦向流升力計算模型,根據氣流迎角的不同,導出了升力計算公式,可以估算半轉翼前進飛行狀態下的升力大小。

(3)應用導出的表達式計算不同飛行條件下半轉翼的升力,獲得了半轉翼懸停和前進兩種飛行狀態下升力變化規律。懸停飛行時,隨著曲柄轉速的增大,升力近似按與轉速的平方成正比的規律增大;前進飛行時,半轉翼升力主要取決于曲柄轉速,與平行于翼面的前進速度關系不大。

(4)比較計算升力曲線與CFD數值仿真所得的升力曲線,兩者變化趨勢一致,最大誤差小于20%。驗證了升力估算解析法的有效性和可行性,可為分析半轉翼飛行器的參數設計和升力預估提供重要的、便于應用的理論分析工具。

[1] 邱支振. 半轉機構[M]. 合肥:中國科學技術大學出版社,2011. QIU Zhizhen.Half-rotating Mechanism [M].Hefei: Press of University of Science and Technology of China,2011.

[2] YOON S, KANG L H, JO S. Development of Air Vehicle with Active Flapping and Twisting of Wing[J]. Journal of Bionic Engineering,2011,8(1):1-9.

[3] NGUYEN Q V, CHAN W L, DEBIASI M. Hybrid Design and Performance Test of a Hovering Insect-inspired Flapping-wing Micro Aerial Vehicle[J]. Journal of Bionic Engineering, 2016,13(2):235-248.

[4] FARUQUE I A, HUMBERT J S.Wing Motion Transformation to Evaluate Aerodynamic Coupling in Flapping Wing Flight[J]. Journal of Theoretical Biology,2014,363:198-204.

[5] SMITH M. Simulating Moth Wing Aerodynamics towards the Development of Flapping Wing Technology[J]. AIAA Journal,1996,34(7):1348-1355.

[6] GURSUL I, HO C M.Oscillating Foils of High Propulsive Efficiency [J]. Journal of Fluid Mechanics,1998,360(1):41-72.

[7] TUNCER I H, WALZ R, PLATZER M. A Computational Study on the Dynamic Stall of a Flapping Airfoil [C]//16th AIAA Applied Aerodynamics Conference. Albuquerque,1998:219-225.

[8] FREYMUTH P. Thrust Generation by an Airfoil in Hover Modes [J]. Experiments in Fluids,1990,9(1/2):17-24.

[9] BRICH J M, DICKINSO M H. The Influence of Wing-wake Interactions on the Production of Aerodynamic Forces in Flapping Flight [J]. Journal of Experimental Biology, 2003,206(13):2257-2272.

[10] 劉嵐,方宗德,侯宇,等. 微型撲翼飛行器的氣動建模分析與試驗[J].航空動力學報,2005,20(1):22-28. LIU Lan, FANG Zongde, HOU Yu, et al. Aerodynamic Modeling and Analysis of Flapping-wing MAV [J]. Journal of Aerospace Power,2005,20(1):22-28.

[11] 周驥平,朱興龍,周建華,等.仿生撲翼飛行簡化力學模型及其實驗研究[J].中國機械工程,2007,18(6):631-635. ZHOU Jiping, ZHU Xinglong, ZHOU Jianhua, et al. Simple Mechanics Model and Experimental Research on Flapping-wing Based on Bionics [J]. China Mechanical Engineering, 2007,18(6):631-635.

[12] WANG Z J, BIRCH J M, DICKINSO M H. Unsteady Forces and Flows in Low Reynolds Number Hovering Flight: Two-dimensional Computations vs Robotic Wing Experiments [J]. Journal of Experimental Biology,2004,207(Pt3):449-460.

[13] 魏瑞軒,胡明朗,郭慶,等.仿鳥撲翼飛行器動力學建模[J].系統仿真學報,2009,21(15):4811-4815. WEI Ruixuan, HU Minglang, GUO Qing, et al. Dynamics Modeling of Bird-like Flapping Wing Air Vehicle [J]. Journal of System Simulation, 2009,21(15):4811-4815.

[14] RAMAMURTI R, SANDBERG W, LOHNER R. Simulation of Flow about Flapping Airfoils Using a Finite Element Incompressible Flow Solver [J]. AIAA Journal, 2001,39(2):253-258.

[15] 楊文青,宋筆鋒,宋文萍,等.仿生微型撲翼飛行器中的空氣動力學問題研究進展與挑戰[J]. 實驗流體力學,2015,29(3):1-10. YANG Wenqing, SONG Bifeng, SONG Wenping, et al. The Progress and Challenges of Aerodynamics in the Bionic Flapping-wing Micro Air Vehicle[J].Journal of Experiments in Fluid Mechanic,2015,29(3):1-10.

[16] 胡明朗, 魏瑞軒,崔曉峰,等. 仿昆撲翼飛行器的翅膀慣性力分析[J]. 航空動力學報,2008,23(7):1279-1286. HU Minglang, WEI Ruixuan , CUI Xiaofeng, et al. Inertia Force of Flapping Wing in Entomopter Micro Air Vehicle[J]. Journal of Aerospace Power,2008,23(7):1279-1286.

[17] COMBES S A, DANIEL T L. Into Thin Air: Contributions of Aerodynamic and Inertial-elastic Forces to Wing Bending in the Hawkmoth Manduca Sexta[J].The Journal of Experimental Biology,2003,206(17):2999-3006.

[18] DELAURIER J D. The Development and Testing of a Full-scale Ornithopter [J]. Canadian Aeronautics and Space Journal,1999,45(2):72-82.

[19] 邱晗,王孝義.半轉機構——一種運動仿生機構的構成及其基本運動特性[J].機械科學與技術,2011,30 (2):600-604. QIU Han, WANG Xiaoyi. Half-rotating Mechanism: a Biomimetic Mechanism of Animal Motion and Its Basic Motion Characteristics[J]. Mechanical Science and Technology for Aerospace Engineering,2011,30(2):600-604.

[20] 姜海波,曹樹良,程忠慶.平板大攻角繞流升力和阻力系數的計算[J].應用力學學報,2011,28(5):518-520. JANG Haibo, CAO Shuliang, CHENG Zhongqing. Lift and Drag Coefficients of Flow around a Flat Plate at High Attack Angles [J]. Chinese Journal of Applied Mechanics,2011,28(5):518-520.

(編輯 陳 勇)

Lift Estimation of HRWs in Hovering and Forward Flights

WANG Xiaoyi ZHANG Yuhua DONG Yinping QIU Han CHEN Fuqiang QIU Zhizhen

School of Mechanical Engineering, Anhui University of Technology, Ma’anshan, Anhui,243002

The lift models of HRW in hovering and forward flights were proposed based on analyses of the motion models of HRW and flow characteristics on the surfaces of the wing. According to the motion characteristics of HRW, the analytical expressions of lift were further derived. The lifts under different flight conditions were calculated using derived formula and numerical simulation based on CFD software respectively with same HRW prototype parameters, by which the changing rule of HRW lift in hovering and forward flights could be also got. The comparisons of lift curves between the theoretical calculation and the numerical simulation demonstrated that the analytical method to estimate lift of HRW is effective and feasible. The research results mentioned above may provide important theoretical guidance for the parameter designs and lift prediction of the HRW air vehicles.

half-rotating wing(HRW); hovering flight; forward flight; lift estimation

2016-10-10

國家自然科學基金資助項目(51375014)

V212

10.3969/j.issn.1004-132X.2017.15.005

王孝義,男,1970年生。安徽工業大學機械工程學院教授、博士。主要研究方向為仿生機械、數字化設計與制造。E-mail:wangxy@ahut.edu.cn。張玉華,男,1961年生。安徽工業大學機械工程學院教授。董銀萍,女,1989年生。安徽工業大學機械工程學院碩士研究生。邱 晗,男,1976年生。安徽工業大學機械工程學院實驗師。陳富強,男,1960年生。安徽工業大學機械工程學院副教授。邱支振,男,1946年生。安徽工業大學機械工程學院教授。

主站蜘蛛池模板: 最新国产高清在线| 色视频国产| 91九色国产在线| 无码福利日韩神码福利片| 国产成人精品视频一区二区电影 | 国产99欧美精品久久精品久久| 黄色网站不卡无码| 999精品色在线观看| 91在线一9|永久视频在线| 无码久看视频| 日本欧美午夜| 亚洲国产清纯| 日韩av在线直播| 91娇喘视频| 99ri国产在线| 99无码中文字幕视频| 欧美啪啪一区| 亚洲国产第一区二区香蕉| 播五月综合| 人妻无码中文字幕第一区| 久热精品免费| www.亚洲国产| 久久久久青草大香线综合精品| 激情六月丁香婷婷| 国产欧美视频在线观看| 午夜视频日本| 亚洲综合经典在线一区二区| 国产精品久久久久婷婷五月| 97人妻精品专区久久久久| 久久天天躁狠狠躁夜夜2020一| 青青网在线国产| 91毛片网| 亚洲A∨无码精品午夜在线观看| 色婷婷狠狠干| 免费又黄又爽又猛大片午夜| swag国产精品| 亚洲成av人无码综合在线观看| 暴力调教一区二区三区| 国产又爽又黄无遮挡免费观看| 国产剧情一区二区| 白浆免费视频国产精品视频 | 伊人天堂网| 色一情一乱一伦一区二区三区小说| 国内丰满少妇猛烈精品播| 亚洲美女视频一区| 一本大道香蕉久中文在线播放| 五月综合色婷婷| 伊人精品成人久久综合| 一级爆乳无码av| 日韩精品欧美国产在线| 全部免费毛片免费播放| 亚洲人成网址| 麻豆精品在线视频| 欧洲欧美人成免费全部视频| 亚洲欧美天堂网| 中文字幕天无码久久精品视频免费 | 国产美女免费| 亚洲免费成人网| 日韩无码视频播放| 欧美区日韩区| 欧美精品一区在线看| 中文字幕首页系列人妻| 国产丝袜无码精品| 青草91视频免费观看| 亚洲国产亚洲综合在线尤物| 亚洲天堂首页| 国产麻豆另类AV| 国产尤物在线播放| 欧洲成人在线观看| 影音先锋丝袜制服| 亚洲国产成人自拍| 久久国产精品波多野结衣| 国产午夜人做人免费视频中文 | 97成人在线视频| 亚洲无码精彩视频在线观看| 久久久精品久久久久三级| 亚洲成人黄色网址| 久久精品最新免费国产成人| 国产一区二区三区精品久久呦| 在线观看欧美国产| 国产成人精品2021欧美日韩| 国产午夜不卡|