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

攻角拉起時前體非對稱渦誘導機翼搖滾運動

2015-03-19 08:24:52徐思文鄧學鎣王延奎
北京航空航天大學學報 2015年11期
關鍵詞:模型

徐思文,鄧學鎣,王延奎

(北京航空航天大學 航空科學與工程學院,北京 100191)

采用細長體機身的飛行器在大攻角下會產生前體非對稱渦,進而誘導出非指令的搖滾運動[1].現代戰機包括 F-18 HARV[2]和 X-31[3]均會產生前體非對稱渦誘導的搖滾運動,尤其是X-31[3]隨攻角會呈現包括滾轉偏離、極限環振蕩和發散等多種運動形態.Brandon和 Nguyen[4]采用細長體和小后掠翼(Λ=26°)的組合體模型首次對前體非對稱渦誘導機翼搖滾進行了實驗研究.隨后采用類似的模型,Ericsson[5-6]、馬寶峰[7]和鄧學鎣等[8]都對這類搖滾運動進行了研究.研究發現,快速建立的極限環振蕩是前體非對稱渦誘導搖滾的典型運動形態,前體非對稱渦左右渦位隨滾轉角的切換被認為是維持極限環機翼搖滾運動所必不可少的條件[4-6,9].另外,人工擾動的研究明顯促進了對于前體非對稱渦誘導搖滾問題的研究:源于模型頭尖部幾何加工公差(自然擾動)的不確定性將使得細長體非對稱渦形態也具有不確定性[10-11],而在細長體頭尖部添加的人工擾動[11-15]的周向位置與前體非對稱渦渦型間具有確定的相關關系[16-17],從而可以解決前體非對稱渦的不確定性問題.在此基礎上,馬寶峰[7]和鄧學鎣等[8]發現在細長體頭尖部添加微顆粒人工擾動將使得自然擾動下原本沒有確定性的自由搖滾運動變得具有確定性,并開展了人工擾動周向位置與前體非對稱渦誘導搖滾運動相關關系的研究.

以上關于前體非對稱渦誘導搖滾問題的研究都是在靜態大攻角下展開的[1-9],然而現代戰機通常是在攻角拉起機動中經歷大攻角的.比如Su-27的“眼鏡蛇”機動[18]就要求在2 s左右的時間內將戰機從小攻角拉起到 90°攻角;X-31“Herbst”機動[19]的第一步也是要求將戰機拉起到過失速的大攻角.因此,本文著重研究攻角拉起時前體非對稱渦誘導非指令搖滾運動及其機理.

1 實驗模型及設備

實驗在北京航空航天大學流體所D4低速風洞展開,風洞來流湍流度為0.08%.開口實驗段尺寸為:2.5m ×1.5m ×1.5m.文中實驗的來流風速固定在Uo=25m/s,相對模型圓柱段直徑(D=90mm)的雷諾數 Re=UoD/ν≈1.6 ×105,ν為來流運動黏性系數.模型通過尾支撐固定在圖1(a)所示的動態拉起機構上,該機構能夠驅動模型以恒定的角速度進行拉起運動,拉起攻角范圍為α =0°~90°,拉起速度范圍為 ω =0.5 ~75(°)/s,對應的拉起減縮頻率范圍為 ψ=ω·L/(360·Uo)=4.0 ×10-5~6.0 ×10-3,其中L=720mm為模型長度.攻角運動由松下伺服電機驅動,連接模型支桿的搖臂和電機軸之間通過傳動比分別為1∶5和1∶2的兩級齒輪進行傳動;攻角由與電機同軸的絕對式編碼器(Koyo-TRD-NA10NW,精度10 bit)進行測量,其測量的攻角精度為0.035°.在攻角動態拉起時,通過引入比例-積分-微分(PID)的控制方法,使得機構對模型攻角的控制誤差不超過0.5°.

實驗模型如圖1(b)所示.模型機身為一個尖拱旋成體,30°后掠的梯形翼作為產生滾轉力矩的前體非對稱渦作用面安裝在距模型尖端x/D=-4~ -5.5的位置,機翼厚度 4 mm.在 x/D=-2.5的位置有24個等間距測壓孔;模型迎風駐點位置定義為測壓孔的0°/360°方位角 θ,如圖1(b)中A-A視圖所示;截面側向力系數Cy由式(1)計算,其中 Cp為測壓孔的壓力系數.在x/D=-4.85截面,兩側機翼各均布6個測壓孔,截面滾轉力矩系數Cl由式(2)計算.拉起運動的轉心位于x/D=-5.5截面.模型頭尖部可以通過安裝在前體的電機進行驅動.整個實驗模型由鋁制成,繞x軸的轉動慣量約為0.007 kg·m2.

式中:llocal為測壓孔所在截面半展長.

直徑為0.2mm的球形顆粒[20]作為人工擾動被黏貼于模型頭尖部,如圖1(c)所示.人工顆粒擾動的軸向黏貼位置和周向位置γ的定義也如圖1(c)所示.在本文中人工擾動固定于γ=45°,此時前體非對稱渦呈現為左側渦位更低(沿x方向)的左渦型[15,20].

圖1 實驗裝置、模型及人工擾動Fig.1 Experimental set up,model and artificial perturbation

自由搖滾運動時,模型與自由搖滾支桿[7]連接,安裝在自由搖滾支桿尾端的編碼器(12 bit)測得的滾轉角精度為0.088°.在自由搖滾實驗之前模型通過支桿內的電磁閘固定在滾轉角φ=0°處.在壓力測量時,自由搖滾支桿將被替換為強迫搖滾支桿[9].通過安裝在強迫搖滾支桿內的伺服電機,驅動模型復現自由搖滾時間歷程,從而能夠排除測壓管道對搖滾運動的影響,實現拉起搖滾過程中的瞬時同步壓力測量.強迫搖滾時模型滾轉角由安裝于支桿尾部的編碼器測量,滾轉角的測量精度為0.018°;強迫搖滾系統的滾轉角復現誤差不超過1.5°,該數值小于自由搖滾重復性實驗時的滾轉角誤差.壓力測量采用量程為±1 psi(6.895 kPa)、精度為 0.1%FS 的 64 通道壓力掃描閥,其最大采集頻率為330Hz,文中動態測壓時的最大采集頻率為18Hz,靜態測壓時的采集頻率為200Hz.模型表面測壓孔與掃描閥之間通過外徑1.0mm的鋼管和軟管連接,管路的長度不超過80 cm;根據曹博超[21]的研究,100 cm長度的管路對10Hz以下壓力信號的濾波和遲滯作用很弱,而文中所關注的非對稱渦所誘導的側向力及滾轉力矩等壓力信號的頻率小于10Hz,因此本文認為所測量的壓力結果真實可信.文中展示的動態壓力和靜態壓力數據均為多次采樣后的平均結果,采樣次數分別為7次和600次.

2 實驗結果與討論

2.1 攻角拉起時搖滾運動形態

2.1.1 攻角靜態時搖滾運動隨攻角的變化曲線

在之前關于前體非對稱渦誘導搖滾運動隨攻角演化規律[4]的研究中并沒有提及不確定性問題,但實驗結果證實這一問題是實際存在的.圖2(a)給出了相同完全外形的兩個不同頭尖部下,搖滾運動隨攻角的演化規律.圖中的數據點表示搖滾運動穩定后第5 s至第20 s搖滾運動的平衡位置,豎線表示第5 s至第20 s搖滾運動的均方差,位于φ=±90°的數據點表示自0°滾轉角釋放的模型將越過90°或-90°滾轉角呈現為發散的運動形態.從圖2(a)中可以看出在非對稱渦產生的27.5°~70°的攻角范圍[22]內,兩條曲線具有明顯的不同.然而,如果在頭尖部黏貼人工擾動,并將人工擾動固定于某一位置(γ=45°),自由搖滾隨攻角的變化規律就是確定的,如圖2(b)所示.

圖2 外形相同的兩個頭尖部下,搖滾運動平衡位置及均方差隨攻角的變化規律(ω=0(°)/s)Fig.2 Mean roll angle and standard deviation of free-to-roll motion with angle of attack for two nose tips with the same shape(ω =0(°)/s)

由于前體非對稱渦隨攻角的變化[22-25],在27.5°~70°攻角范圍內的搖滾運動形態隨攻角增加而變化,如圖3(a)所示;其中在31°~47.5°為較為劇烈的發散和振蕩發散運動,運動時間歷程分別如圖3(b)、圖3(c)所示,在其余攻角下則呈現為小振幅的滾轉偏離和極限環振蕩.在前體非對稱渦對稱的15°攻角附近,存在著機翼本身流動引起的滾轉偏離運動.

圖3 攻角靜態時,搖滾運動隨攻角的變化規律及發散、振蕩發散時間歷程曲線Fig.3 Free-to-roll motion changing with angle of attack and time histories for divergence as well as oscillating divergence when angle of attack is static

2.1.2 不同拉起速度下搖滾運動隨攻角變化曲線

為研究攻角拉起速度對搖滾運動的影響,本節在0°~85°的攻角范圍內研究了系列拉起速度下的搖滾運動形態.圖4給出了3期實驗所得搖滾運動振幅隨拉起速度的變化曲線,其中振幅是單次拉起搖滾運動實驗中滾轉角絕對值的最大值,圖5給出了不同拉起速度下搖滾運動隨攻角的變化曲線.從圖中可以看出,在3個不同拉起速度分區內呈現出3種不同拉起搖滾運動形態:①當攻角拉起速度位于0.5~5(°)/s范圍內時,搖滾運動隨攻角會呈現出偏離、發散和極限環振蕩等形態,如圖5(a)所示;②當攻角拉起速度位于10~25(°)/s時,搖滾運動隨攻角呈現為發散的運動形態,如圖5(b)所示;③當攻角速度位于32.5~75(°)/s時,搖滾運動呈現為一個周期左右的類正弦振蕩形式,搖滾振幅在30°左右,運動形態相比前兩個拉起速度分區下的運動形態較為緩和,15°攻角附近的滾轉偏離也不再產生,如圖5(c)所示.另外,從圖4的結果也可以看出,當拉起速度位于3個拉起速度分區之間時,搖滾運動形態不具有重復性.

圖4 搖滾運動振幅隨攻角拉起速度的變化曲線Fig.4 Curves of amplitude of free-to-roll motion changing with pitch rate

圖5 不同拉起速度下搖滾運動隨攻角的變化曲線Fig.5 Curves of free-to-roll motion changing with angle of attack for different pitch rates

2.1.3 拉起減縮頻率對拉起搖滾運動的影響

對于拉起速度對搖滾運動影響的研究是否具有普遍性的意義,圖6給出了兩個不同無量綱拉起減縮頻率ψ下的拉起搖滾運動曲線.從圖6中可以看出,盡管來流風速Uo和拉起速度ω都發生了變化,但只要拉起減縮頻率維持不變,搖滾運動形態就大致相同;而當拉起減縮頻率改變后,搖滾運動形態即發生明顯的變化.在相同拉起減縮頻率下,圖6中不同拉起速度(來流風速)下的搖滾曲線在幅值和相位上仍然有一定差異,比如圖6(b)中,拉起速度(來流風速)越大,幅值就越大,除了搖滾運動自身的重復性影響以外,另一個重要的原因為支桿摩擦力:來流風速越大,驅動搖滾的滾轉力矩的數值就越大,那么摩擦力的影響就相對越小,從而搖滾的振幅就相對較大、可以想象,當風速降為0.5m/s甚至更小時,較小的滾轉力矩將不能克服支桿的摩擦力作用,模型就不會再產生拉起搖滾運動.

圖6 不同拉起減縮頻率下的搖滾運動Fig.6 Free-to-roll motions for different reduced frequencies

除了無量綱拉起減縮頻率以外,來流Re數由于決定了前體非對稱渦的分離類型和強度[26]也必然對拉起搖滾運動形態有重要的影響,但由于本文Re數小于2.0×105,前體流動始終處于亞臨界Re區(前體兩側邊界層分離形態為層流),因此Re數效應還沒有顯現出來.

2.2 攻角拉起時搖滾運動的機理分析

2.2.1 搖滾運動隨拉起速度分區效應的機理分析

從第2.1.2節的實驗結果可以看出,當攻角拉起速度增加到較快的第③分區(32.5~75(°)/s)時,搖滾運動的形態由發散變為振幅約30°左右的類正弦振蕩,運動形態突然變得“緩和”,這一顯著變化可從圖7得出原因.

圖7給出了攻角靜態和兩個典型拉起速度下搖滾運動隨攻角的變化曲線,從圖中可以看出發散運動發生的攻角范圍為31°~47.5°,拉起速度位于第②分區10~25(°)/s時的發散運動均發生在該攻角范圍內.對于搖滾來說,其滾轉角時間歷程可表達為

即搖滾的振幅由滾轉力矩及搖滾時間決定,快速拉起時不再發散也由這兩個因素引起.圖8(a)給出了不同拉起速度下機翼背風面滾轉力矩系數隨攻角的變化曲線,隨著拉起速度增加,攻角31°~47.5°范圍內的力矩系數有所減小,但沒有本質的差異.圖8(b)則給出了攻角靜態和攻角拉起時的發散時間歷程曲線,從圖中可以看出攻角靜態和動態時發散運動時間歷程曲線相近,根據實驗測得的所有發散運動時間歷程曲線的統計,模型從0°滾轉角發散到 -90°滾轉角平均需要0.585 s;因此,當拉起速度大于 28.2(°)/s(計算公式為(47.5°-31°)/0.585 s)時,模型即“沖過”發散的31°~47.5°攻角區域直接進入滾轉偏離的攻角區域,使得拉起搖滾呈現為0°滾轉角附近的類正弦運動形態,如圖7所示.這一分析所得的拉起速度28.2(°)/s正位于第②、第③拉起速度分區之間,也證實了這一分析的合理性.因此,隨著拉起速度增加,搖滾時間的減少,是快速拉起時搖滾沖過發散區不再發散的主要原因.在拉起速度較快時,攻角15°附近滾轉偏離的消失也是同樣的原因.

圖7 不同拉起速度時搖滾運動隨攻角變化曲線對比Fig.7 Comparison of curves of free-to-roll motions changing with angle of attack for different pitch rates

圖8 拉起速度對機翼背風側截面滾轉力矩系數和搖滾時間歷程曲線的影響Fig.8 Effects of pitch rate on leeward sectional rolling moment coefficient and on time history of free-to-roll motion

2.2.2 快速拉起類正弦振蕩的機理分析

圖9 攻角拉起類正弦振蕩與單獨攻角拉起(φ =0°)兩種情形下 Cy-α 及Cl-α 曲線對比(ω=40(°)/s)Fig.9 Comparison of Cy-α and Cl-α curves with φ=0°in pitch-up and during sinusoidal-like oscillation in pitch-up(ω =40(°)/s)

以往對于攻角靜態時搖滾運動的機理分析,往往通過流動及滾轉力矩等隨滾轉角的變化特性來展開.但對于攻角快速拉起時的搖滾,由于涉及到滾轉和攻角拉起兩個自由度的運動,其分析也將有所不同.圖9給出了快速拉起機翼搖滾過程中Cy-α及Cl-α曲線與固定0°滾轉角單獨快速拉起(拉起速度與前者相同)過程中Cy-α及Cl-α曲線對比.圖9(a)表明兩種情形下側向力隨攻角的變化曲線基本相同,即主控搖滾的前體非對稱渦在搖滾過程中沒有發生明顯的變化;圖9(b)所示兩種情形下滾轉力矩系數隨攻角的變化曲線則表明:對于背風側的滾轉力矩,滾轉運動對其也沒有本質上的影響.因此,在快速拉起時,前體非對稱渦及其誘導滾轉力矩隨滾轉角的變化對于搖滾的分析相對次要,二者隨攻角的變化更為重要.

另外,從圖9中也可以看出,在攻角超過70°以后,由于前體流動逐漸進入類卡門渦區[22-25],滾轉力矩接近于0,因此快速拉起時70°攻角以后的類正弦振蕩應該是搖滾的一個收斂過程.

3 結論

1)搖滾運動隨攻角的變化規律隨攻角拉起速度呈現出明顯的分區特性:在拉起速度位于0.5~5(°)/s的第①分區時,搖滾運動隨攻角的演化規律與攻角靜態時的結果相似,呈現出包括偏離、極限環振蕩以及發散等運動形態;在拉起速度位于10~25(°)/s的第②分區時,模型連續旋轉;而拉起速度位于32.5~75(°)/s的第③分區時,非指令搖滾運動形態呈現為周期數較少的類正弦形式振蕩;拉起速度位于3個分區之間時,搖滾運動不具有重復性.

2)攻角拉起的無量綱減縮頻率是主控攻角拉起時搖滾運動形態的重要無量綱參數,在亞臨界Re數下,拉起無量綱減縮頻率不變,拉起搖滾運動的形態就基本維持不變.

3)當拉起速度較大時,模型將沒有足夠的時間完成從0°~90°滾轉角的發散運動,這是搖滾運動形態隨攻角拉起速度具有分區效應的關鍵原因.

4)與固定攻角時自由搖滾不同,流動隨滾轉角的變化對快速拉起時自由搖滾運動并沒有決定性的影響,取而代之的是流動隨攻角的變化,即Cl-α曲線取代Cl-φ曲線,成為快速拉起搖滾分析的關鍵.

References)

[1] Katz J.Wing/vortex interactions and wing rock[J].Progress in Aerospace Sciences,1999,35:727-750.

[2] Quast T,Nelson R C,Fisher D F.A study of high alpha dynamics and flow visualization for a 2.5%model of the F-18 HARV undergoing wing rock[J].AIAA,1991,3267-CP:524-533.

[3] Klein V,Noderer K D.Aerodynamic parameters of the X-31 drop model estimated from flight data at high angles of attack[J].AIAA,1992,4357:174-181.

[4] Brandon JM,Nguyen L T.Experimental study of effects of forebody geometry on high angle of attack stability[J].Journal of Aircraft,1988,25(7):591-597.

[5] Ericsson L E.Wing rock generated by forebody vortices[J].Journal of Aircraft,1989,26(2):110-116.

[6] Ericsson L E.Further analysis of wing rock generated by forebody vortices[J].Journal of Aircraft,1989,26(12):1098-1104.

[7]馬寶峰.前體渦誘導機翼搖滾的實驗研究[D].北京:北京航空航天大學,2007.Ma B F.Experimental investigation of roll oscillation induced by forebody vortex[D].Beijing:Beihang University,2007(in Chinese).

[8] Wang B,Deng X Y,Ma B F,et al.Effect of tip perturbation and wing locations on rolling oscillation induced by forebody vortices[J].Acta Mechanica Sinica,2010,26(5):787-791.

[9]榮臻.前體渦誘導機翼搖滾的流動特性及機理研究[D].北京:北京航空航天大學,2009.Rong Z.An experimental investigation on flow characteristicsand mechanism of wing rock induced by forebody vortex[D].Beijing:Beihang University,2009(in Chinese).

[10] Dexter P C,Hunt B L.The effect of roll angle on the flow over a slender body of revolution at high angles of attack[J].AIAA,1981,0358:1-13.

[11] Zilliac G G,DeganiD,Tobak M.Asymmetric vortices on a slender body of revolution[J].AIAA,1991,29(5):667-675.

[12] Dexter PC.A study of asymmetric flow over slender bodies at high angles of attack in a low turbulence environment[J].AIAA,1984,0505:1-11.

[13] Moskovitz C A,Raleigh N C,Hall R M,et al.Effects of surface perturbations on the asymmetric vortex flow over a slender body[J].AIAA,1988,0483:1-10.

[14] Moskovitz C A,Hall RM,Dejarnette FR.Effects of nose bluntness,roughness and surface perturbations on the asymmetric flow past slender bodies at iarge angles of attack[J].AIAA,1989,2236-CP:720-732.

[15] Chen X R,Deng X Y,Wang Y K,et al.Influence of nose perturbations on behaviors of asymmetric vortices over slender body[J].Acta Mechanica Sinica,2002,18(6):581-593.

[16] Hunt B L.Asymmetric vortex forces and wakes on slender bodies[J].AIAA,1982,1336:1-41.

[17] Lamont P J.Pressures around an inclined ogive cylinder with laminar,transitional,or turbulent separation[J].AIAA,1982,20(11):1492-1499.

[18] Li Y P,Ge F Y,Lan CE.Estimation of aircraft models in cobra maneuver through dynamic inversion[J].AIAA,1994,3494-CP:322-332.

[19] Alcorn C W,Croom M A,Francis M S,et al.The X-31 aircraft:Advances in aircraft agility and performance[J].Progress in Aerospace Sciences,1996,32:377-413.

[20]田偉.細長旋成體及融合體型機身大迎角背渦流動特性研究[D].北京:北京航空航天大學,2010.Tian W.Study on behaviors of leeward vortices over slender body and chined fuselage at high angle of attack[D].Beijing:Beihang University,2010(in Chinese).

[21]曹博超.風洞動態測壓技術及其在機翼搖滾運動研究中的應用[D].北京:北京航空航天大學,2008.Cao BC.Dynamic pressure measurement technology in the wind tunnel and its application in the investigation of wing-rock phenomenon[D].Beijing:Beihang University,2008(in Chinese).

[22]王剛,鄧學鎣,王延奎,等.亞臨界雷諾數細長體繞流流態隨迎角的變化和分區[J].流體力學實驗與測量,2003,17(2):19-26.Wang G,Deng X Y,Wang Y K,et al.Zonal study of flow patterns around an ogive-cylinder at subcritical reynolds numbers[J].Experiments and Measurements in Fluid Mechanics,2003,17(2):19-26(in Chinese).

[23] Lamont P J,Hunt B L.Pressure and force distributions on a sharp-nosed circular cylinder at large angles of inclination to a uniform subsonic stream[J].J Fluid Mech,1976,76(3):519-559.

[24] Keener E R,Chapman L C,Talegbani J.Side forces on a tangent ogive forebody with a fineness ratio of3.5 at high angles of attack and mach numbers from 0.1 to 0.7[J].NASA,1977,TM,X-3437:1-109.

[25] Ericsson L E,Reding J P.Aerodynamic effects of asymmetric vortex shedding from slender bodies[J].AIAA,1985,1797:222-256.

[26] Lamont P J.Pressures around an inclined ogive cylinder with laminar,transitional,or turbulent[J].AIAA,1980,1556R:1-10.

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 91亚瑟视频| 国产乱子伦一区二区=| 国产一级片网址| 麻豆国产精品一二三在线观看| 露脸国产精品自产在线播| 亚洲综合在线最大成人| 精品一区二区无码av| 67194成是人免费无码| 伊人精品视频免费在线| 精品国产www| 中文字幕不卡免费高清视频| 国产成在线观看免费视频 | 欧美日韩国产在线播放| 波多野结衣国产精品| 欧美精品一区二区三区中文字幕| 欧美爱爱网| 欧美午夜理伦三级在线观看 | 精品无码一区二区三区在线视频 | 欧美国产日韩在线播放| 欧美国产日韩一区二区三区精品影视| 国产网站在线看| 幺女国产一级毛片| a网站在线观看| 欧美激情福利| yy6080理论大片一级久久| 国产手机在线ΑⅤ片无码观看| 亚洲无码在线午夜电影| 在线看国产精品| 久久96热在精品国产高清| 欧美 国产 人人视频| 国产微拍一区二区三区四区| 欧美日韩国产成人在线观看| 丰满人妻久久中文字幕| av一区二区三区在线观看| 亚洲无码高清视频在线观看 | 色播五月婷婷| 2021最新国产精品网站| 亚洲h视频在线| 国产精品视频猛进猛出| 国模沟沟一区二区三区| 综合久久久久久久综合网| 国产精品高清国产三级囯产AV| 国产人在线成免费视频| 婷婷亚洲综合五月天在线| 亚洲无线视频| 国产精品一线天| 久久综合成人| 国产精品所毛片视频| 日本一区二区不卡视频| 久久亚洲国产视频| 日韩精品少妇无码受不了| 国产成人在线小视频| 黄色网页在线播放| 97青草最新免费精品视频| 成年人视频一区二区| 国产精品亚洲综合久久小说| 91网站国产| 大学生久久香蕉国产线观看 | 国产成人免费视频精品一区二区| 欧美高清国产| 婷婷午夜天| 伊伊人成亚洲综合人网7777| 欧美人与性动交a欧美精品| 少妇高潮惨叫久久久久久| 毛片免费观看视频| 国产高清精品在线91| 亚洲精品自产拍在线观看APP| 国产人免费人成免费视频| 久久a毛片| 麻豆国产原创视频在线播放 | 91一级片| 欧美亚洲欧美| 中文字幕 91| 欧美性精品| 婷婷丁香色| 手机永久AV在线播放| 精品伊人久久久香线蕉 | 国产精品自在在线午夜| 97成人在线观看| 久久6免费视频| 亚洲无码37.| 视频一本大道香蕉久在线播放 |