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

低速旋轉尾翼式彈箭氣動特性數值研究

2018-01-08 02:21:27王學德曹友銓
彈道學報 2017年4期

張 超,王學德,王 勇,曹友銓

(南京理工大學 能源與動力工程學院,江蘇 南京 210094)

低速旋轉尾翼式彈箭氣動特性數值研究

張 超,王學德,王 勇,曹友銓

(南京理工大學 能源與動力工程學院,江蘇 南京 210094)

為了研究低速旋轉對尾翼式彈箭氣動特性的影響,采用三維非定常N-S方程并結合滑移網格技術,在小攻角和全馬赫數下,對某尾翼彈在低轉速狀態下的繞流流場進行了數值模擬。以美國陸-海軍動導數計算標模驗證該文算法的有效性。結果表明該方法有較高的精確度。由不同馬赫數、轉速和滾轉角條件下的計算結果發現:縱向氣動特性(即升力、阻力、俯仰力矩)不隨轉速而變化,平均滾轉力矩系數和轉速為定比例關系,平均馬格努斯力系數隨轉速呈非線性變化,瞬時馬格努斯力系數隨滾轉角呈正弦變化。

彈箭;馬格努斯效應;滑移網格;氣動特性;數值計算

彈箭通常采用繞彈體縱軸旋轉的飛行方式以提高射擊密集度以及飛行穩定性,然而,旋轉的同時也會帶來諸如馬格努斯效應、陀螺效應等問題。相比于無旋彈箭,旋轉彈箭的空氣動力特性非常復雜。因此,研究旋轉彈箭相關氣動特性對彈箭設計,尤其是無控彈箭[1],有著極其重要的影響。

在21世紀以前,人們主要通過實驗和工程估算方法研究自旋彈箭的氣動特性。伴隨著計算機技術性能的大幅度提高以及相關數值計算軟件的日臻完善,加之具有計算周期短、精度高、成本低廉等優勢,這使得利用成熟的計算軟件對彈箭進行數值仿真成為主流趨勢。美國航天宇航局[2]分別采用風洞實驗、工程估算和數值計算3種方式對自旋彈丸的氣動特性進行了專門的研究。Sturek W B與Pechier M[3-4]對尖拱-圓柱型旋轉彈丸同時進行了風洞試驗和數值仿真,并將兩者的結果做了對比,發現兩者吻合程度較高,尤其是在小攻角范圍內。在國內,吳甲生[5]對旋轉彈箭空氣動力學效應進行了系統深入的分析,并指出旋轉彈箭未來的研究重點以及發展方向。雷娟綿[1]對大攻角、高轉速條件下的SOCBT(Secant-Ogive-Cylinder-Boattail)進行了數值模擬,并從流場結構闡述了馬格努斯效應的形成機理。劉周[6]采用不同的湍流模型對SOCBT彈丸進行數值模擬,發現不同的湍流模型對計算精度有較大影響。

從公開的文獻來看,國內外對旋轉彈丸進行了廣泛的研究,在高轉速、大攻角條件下,對不同形狀的無翼彈丸采用不同的數值方法和湍流模型進行了大量研究,并對馬格努斯效應作出了相應的機理分析,但缺少在低轉速且小攻角條件下,對尾翼式彈箭氣動特性隨馬赫數、轉速以及滾轉角進行定量分析與研究。本文采用滑移網格技術,先對ANF標模進行數值方法的有效性驗證,而后對某自旋式尾翼彈進行模擬,從給出的壓力云圖對馬格努斯效應產生原因進行分析與解釋,并且定量地給出了馬格努斯力、升力系數等氣動特性隨轉速、滾轉角以及馬赫數的變化關系,為飛行控制系統設計提供了相關依據。

1 數值模擬方法

1.1 滑移網格技術

滑移網格技術是上世紀90年代興起的一種基于非穩態思想的數值模擬方法,近年來得到了廣泛的應用,尤其是在彈箭數值仿真方面。滑移網格原理如圖1所示,整個流動區域被分為包含123單元和456單元的2個不同的區域,其中一個區域相對另一個區域做相對滑移運動來模擬真實的流動過程。ABC和EFG是各子區域與滑移面重合的邊界面。在數值模擬過程中,123單元和456單元在運動過程中會形成B′、C′、D′、E′、F′、G′節點,并沿著交界面做相對滑移運動。其中子塊2的信息通過面C′D′和面D′E′傳遞到子塊4和子塊5,其他單元之間的信息傳遞與此類似。

滑移網格技術作為一種特殊的動網格[7],適用于具有周期性旋轉運動物體的流場計算,而且具有計算精度高、占用內存少等優勢[1]。當模擬彈箭旋轉運動時,只需要一個非旋轉區域和一個包圍彈體的旋轉區域,兩者之間存在一個交界面,交界面的網格節點數不必相同,僅僅在滑移網格交界面進行數值插值,即可保證2個區域之間通量守恒。圖2(a)和圖2(b)給出了網格區域劃分示意圖,圖中d為彈體直徑。

1.2 控制方程

采用基于滑移網格下的三維積分守恒型非定常N-S方程組,由質量守恒方程、動量守恒方程和能量守恒方程組成[1]:

式中:ρ,U,e分別為氣流的密度、速度和單位體積的能量,vm為彈體的運動速度,p為靜壓力,I為單位張量,τ為粘性張量,q為熱通量,S和V分別為控制體的面積分和體積分的積分區域,r為S的外法向單位向量。

1.3 離散方法與湍流模型

采用有限體積法對空間進行離散,對流項二階通量差分分裂(FDS)AUSM+格式,粘性項采用中心差分格式[8]。

Realizablek-ε模型是對RNGk-ε湍流模型做了修正而產生的,能更為準確地反映流體真實的流動過程,尤其是旋轉流動計算中氣流的流動狀態[9]。

1.4 時間步長技術

對自旋彈箭非定常的研究,采用雙時間步長技術,即采用物理時間步長和內循環時間步長。在進行計算時,一般在一個時間步內彈體旋轉0.1°~0.3°比較合適,在這個范圍內既能夠保證流場空間數據的連續性,又能使得滾轉物體有足夠的旋轉角度,這樣可以大大節省收斂所需的時間[10]。

2 算法有效性驗證

2.1 驗證模型

選取ANF標模進行數值算法的有效性驗證,ANF模型經過大量風洞實驗和飛行試驗驗證,被廣泛采用為動導數的計算標模。表1給出計算條件,其中馬赫數Ma、雷諾數Re、總壓p0和總溫T0與風洞實驗條件一致,彈箭旋轉速度Ω=385.8 rad/s,Ω*為對應的無量綱轉速,Ω*=dΩ/c,d為彈體直徑,c為當地聲速,N為外循環中每轉動一周的總計算的迭代次數,Δt為外循環的計算步長,i為每一次外循環計算中內循環的迭代次數。

表1 計算條件

2.2 網格模型

利用ICEMCFD,采用多塊對接網格生成技術對整個計算區域生成六面體網格,計算區域共生成900多萬網格。縱向采用C型網格,橫向采用O型網格。第一層距離物面為1.0×10-5m,網格伸展比為1.1,控制第一層網格到物面的距離和網格伸展比的目的是為了保證y+<1。圖3(a)和圖3(b)分別是彈體頭部和尾部附近網格。

彈箭表面設置為無滑移壁面條件,彈體跟隨旋轉區域以相同的轉速同步運動。非旋轉區域和旋轉區域的交界面設置為滑移邊界條件;非旋轉區域的外邊界設置為壓力遠場條件。

2.3 算法驗證分析

從圖4(a)和圖4(b)可知,俯仰力矩和法向力系數與實驗值以及參考文獻中計算值高度吻合,誤差均在5%以內;從圖4(c)和圖4(d)可知,在攻角α<30°時,馬格努斯力和力矩動導數與實驗值以及參考文獻中的計算值吻合程度較好。當攻角為60°時,本文數值計算結果與實驗值相差為20%,這是由于當攻角過大時,旋轉彈體附近流場變得更為復雜。綜上所述,在小攻角范圍內,本文采用的數值計算方法具有較高的可信度。

3 某低速尾翼彈馬格努斯效應數值計算

3.1 模型、網格和計算狀態

計算模型為某尾翼彈,由錐形部、圓柱段和尾翼組成,其尺寸外形如圖5所示。選取的計算狀態為:來流馬赫數Ma=0.9,1.1,1.4,1.6,2.0,來流溫度T0=288 K。該彈旋轉速度Ω=π,2π,4π rad/s,此時無量綱轉速分別為Ω*=0.003 45,0.006 9,0.013 8,來流攻角α=2°。一個計算周期采用1 440個物理時間步,即Δt=0.000 69 s(3種轉速采用相同的物理時間步和步長)。氣動力和力矩計算的參考點位于彈體質心位置,參考面積為彈體的最大橫截面積,參考長度為彈體直徑d。

3.2 網格收斂性分析

為了研究網格的收斂性,在保證y+≤1的情況下,對來流條件為Ma=1.1,攻角α=2°,轉速Ω*=0.013 8,分別以網格數M=6×106,7×106,8×106,9×106共4種不同網格數對該彈箭流場進行數值模擬。圖6給出了瞬時馬格努斯力系數Cmg隨滾轉角Φ的變化規律,可知:隨著網格密度的進一步增大,馬格努斯力系數的最大幅值差也不斷變小,6×106個網格與9×106個網格的幅值差異最大為20.2%,7×106個網格與9×106個網格幅值最大相差7.8%,8×106個網格與9×106個網格幅值最大相差2.3%,可見隨著網格的加密,瞬時馬格努斯力系數隨滾轉角的變化趨于收斂。為了進一步驗證網格的收斂性,表2給出Ma=1.1,攻角α=2°,轉速Ω*=0.013 8,不同網格數下氣動特性計算結果。可知:平均升力系數CL、平均阻力系數CD、平均馬格努斯力系數Cmg隨網格數M的增大逐漸收斂。因此,后序計算以900萬網格為準。

MCLCDCmg6×1060.652190.336020.0033017×1060.618140.317040.0031418×1060.614060.304890.0030779×1060.600740.300390.003001

3.3 不同轉速下的流場結構

圖7和圖8分別給出了彈箭圓柱段和尾翼段在不同轉速下的壓力分布云圖,其中觀察方向為彈體尾部向前看。由于攻角和彈體自轉的存在,會產生一個繞彈體縱軸x旋轉的環流(環流方向即為彈體旋轉方向)和一個指向y軸正向的橫流。由圖7(a)可知,當彈體旋轉時,彈體左側橫流方向與環流方向相反,這使得左側氣流速度降低,而右側環流與橫流方向相同,則氣流速度增大。最終使得氣流對彈體左側壓力大于其對彈體右側的壓力,即彈體左右存在壓力差,所以產生了指向彈體右側的側向力,即馬格努斯力。綜合對比圖7(a)~圖7(c)發現,隨著轉速的增大,左側高壓區明顯擴大,而右側低壓區也在逐漸擴大,彈體左右兩側壓差不斷增大,所以馬格努斯力也在增大。通過比較不同轉速下的壓力云圖,可知彈體轉速是影響馬格努斯力大小的一個重要因素。從圖8可知,當彈箭處于低速旋轉狀態時,在迎風區域,由于橫流的影響遠遠大于環流的影響,因而該區域的尾翼兩側區域壓力分布較為對稱,所以產生的馬格努斯力也較小。而背風區由于只受到環流的影響,所以尾翼兩側出現較為明顯的壓力分布不對稱現象,產生較大的馬格努斯力。而且隨著轉速的增大,對背風區的影響也要遠遠大于對迎風區的影響。通過以上分析,并且綜合比較圖7和圖8可知:轉速是馬格努斯效應的影響因素之一,并且彈體處馬格努斯力方向與環流方向一致,而尾翼處馬格努斯力方向與環流方向相反。

3.4 不同轉速下氣動力特性分析

通過數值模擬獲得了尾翼彈在各種飛行條件下的平均阻力系數、平均升力系數、平均滾轉力矩系數等氣動特性數據,如圖9所示。

圖9(a)~圖9(c)是不同轉速下尾翼彈的平均阻力系數CD、平均升力系數CL以及平均俯仰力矩系數Cm隨馬赫數的變化曲線,可知:在全聲速范圍內,縱向氣動特性(阻力、升力、俯仰力矩)隨馬赫數變大呈先增大后減小的趨勢,且縱向氣動特性不隨轉速的變化而變化。

圖9(d)是不同轉速下平均滾轉力矩系數Crm隨馬赫數的變化曲線,并且每個馬赫數下的滾轉力矩均為負值,這是由于該導彈模型不存在尾翼斜置角,導彈無法獲得使其主動滾轉的力矩(即導轉力矩為0)而只存在阻礙其旋轉的滾轉阻尼力矩,故而平均滾轉力矩系數均為負值;隨轉速增加,平均滾轉力矩系數增大,與轉速成定比例關系;在同一轉速下,隨馬赫數增大,平均滾轉力矩系數的絕對值呈先增加后減小的變化趨勢,與常規滾轉力矩系數變化規律一致。

圖9(e)給出了不同轉速下平均馬格努斯力Cmg隨馬赫數的變化曲線,可知:平均馬格努斯力系數隨馬赫數的增加呈先增大后減小的趨勢,各個馬赫數下的平均馬格努斯力系數與轉速為非線性增大關系。圖9(f)是轉速Ω*=0.013 8、馬赫數Ma=1.1的條件下,瞬時馬格努斯力隨滾轉角的變化曲線,可知:在一個旋轉周期內有4個小周期,每個小周期的馬格努斯力系數隨滾轉角近似呈現正弦曲線變化。在滾轉角從0°到90°的過程中,馬格努斯力系數先增大,當轉到Φ=27°左右時達到峰值,然后減小,在Φ=69°左右時達到波谷,然后再次增大。瞬時馬格努斯力隨滾轉角變化會產生4個相同周期,這是由于該彈箭的尾翼布局為“十字型”,彈箭每旋轉90°時,此時彈箭再一次處于“十字型”布局,彈箭周圍的流場性質與前一次“十字型”布局流場性質相同。自旋彈箭在飛行過程中,馬格努斯力可以分解為粘性分量(由于流體粘性的存在而產生的那一部分馬格努斯力)和壓差分量(彈體兩側的壓力差產生的那一部分馬格努斯力)。表3給出了在轉速Ω*=0.013 8、馬赫數為1.1且不同滾轉角下,各個部件的粘性分量Cmg,S(通過物面摩擦應力積分得到的馬格努斯力)和壓差分量Cmg,P(通過物面壓力積分得到的馬格努斯力)對馬格努斯力系數Cmg的貢獻值。可知,粘性分量對馬格努斯力系數的貢獻始終為正值(對馬格努斯力系數的最大貢獻比例不超過5%)。當滾轉角Φ=0°,50°,90°時,由于彈體和尾翼產生的馬格努斯力始終是一個正值,而另一個為負值,所以總的馬格努斯力在0值上下浮動。當滾轉角Φ=27°時,彈體和尾翼產生的馬格努斯力均為正值,這使得總的馬格努斯力最大(波峰)。而當滾轉角Φ=68°時,彈體和尾翼產生的馬格努斯力均為負值,這使得總的馬格努斯力最小(波谷),即在一個小周期內馬格努斯力系數隨滾轉角近似呈現正弦曲線變化。

圖9(g)給出轉速Ω*=0.013 8下不同部件馬格努斯力系數隨馬赫數的變化曲線。可知:彈體、尾翼、全彈的平均馬格努斯力系數均隨馬赫數的增大呈先增大后減小的趨勢。并且馬格努斯力主要是由尾翼產生的,尤其是在Ma=0.9,2時,彈體的馬格努斯力幾乎可以忽略。各部件平均馬格努斯力系數如表4所示。表中Δ表示部件提供平均馬格努斯系數百分比。

對比圖9(a)~圖9(c),可知轉速對平均阻力、升力、俯仰力矩系數沒有影響,而非定常計算所需的時間遠遠大于定常計算時間,所以可以采用定常方法來得到升力、阻力、俯仰力矩系數。對比圖9(a)~圖9(e)發現:平均滾轉力矩系數與平均馬格努斯力系數的數值大小約為平均升力的1%~10%,這是由于該彈箭的尾翼屬于“十字型”布局,并且當彈箭不旋轉時,彈箭攻角平面兩側流場對稱分布,不會產生馬格努斯力和滾轉力矩。而當彈體逆時針低速旋轉時,由于旋轉產生環流的線速度遠遠小于攻角產生的橫流速度,所以旋轉引起的氣動系數屬于小量。

表3 不同滾轉角下的馬格努斯力系數

表4 不同部件馬格努斯力系數

4 結論

本文選取ANF標模進行數值算法的有效性驗證。數值計算結果表明,本文選取的計算方法具有較高的計算精度。在此基礎上,研究了某低速旋轉尾翼彈箭氣動特性隨馬赫數、轉速以及滾轉角的變化規律,得到以下結論:

①通過比較不同馬赫數和轉速下的縱向氣動數據,可知縱向氣動特性不隨轉速變化而變化。隨馬赫數增加,縱向氣動特性均呈現先增加后減小的變化規律,與常規靜態氣動變化規律相同。

②隨馬赫數增加,平均滾轉力矩系數絕對值呈現先增加后減小的變化趨勢,與常規滾轉力矩變化規律一致。隨轉速增加,平均滾轉力矩系數增大,且平均滾轉力矩系數與轉速呈現典型的線性變化規律。

③隨馬赫數增加,平均馬格努斯力系數呈現先增加后減小的變化趨勢;隨轉速增加,平均馬格努斯力系數增大,且平均馬格努斯力系數與轉速為非線性變化關系。在同一轉速和馬赫數下,瞬時馬格努斯力系數與滾轉角為近似正弦曲線變化規律,在一個滾轉周期內,瞬時馬格努斯力有4個自然變化周期,每一個周期為90°,即尾翼從十字型布局變化到下一個十字型布局為一個變化周期。瞬時馬格努斯力隨滾轉角的變化在一個較大正值和較小負值間振動。

[1] 雷娟綿,李田田,黃燦.高速旋轉彈丸馬格努斯效應數值研究[J].兵工學報,2013,34(6):718-725.

LEI Juan-mian,LI Tian-tian,HUANG Can.A numerical investigation of magnus effect for high-speed spinning projectile[J].Acta Armamentarii,2013,34(6):718-725.(in Chinese)

[2] DESPIRITO J.Lateral jet interaction on finned projectile in supersonic flow:AIAA-2012-0413[R].Reston:AIAA,2012.

[3] STUREK W B,KAYSER L D,NIETUBICZ C J.Computations of Magnus effects for a yawed,spinning body of revolution[J].AIAA,1978,7(16):687-692.

[4] PECHIER M,GUILLEN P,CAYZAC R.A combined theoretical-experimental investigation of Magnus effects[C]//16th AIAA Applied Aerodynamics Conference.Barges:AIAA,1998:778-788.

[5] 苗瑞生,吳甲生.旋轉彈空氣動力學[J].力學進展,1987,17(4):479-488.

MIAO Rui-sheng,WU Jia-sheng.Aerodynamics of spinning projectiles[J].Advances in Mechanics,1987,17(4):479-488.(in Chinese)

[6] 劉周,謝立軍,楊云軍,等.彈丸旋轉空氣動力效應非定常數值模擬[J].航空學報,2016,37(5):1 401-1 410.

LIU Zhou,XIE Li-jun,YANG Yun-jun,et al.Unsteady numerical simulation of aerodynamic effect of a spinning projectile[J].Acta Aeronautica et Astronautica Sinica,2016,37(5):1 401-1 410.(in Chinese)

[7] 米百剛,詹浩,朱軍.基于CFD數值仿真技術的飛行器動導數計算[J].空氣動力學學報,2014,32(6):834-839.

MI Bai-gang,ZHAN Hao,ZHU Jun.Calculation of dynamic derivatives for aircraft based on CFD technique[J].Acta Aerodynamica Sinica,2014,32(6):834-839.(in Chinese)

[8] 趙博博,劉榮忠,郭銳,等.扭曲尾翼彈箭的馬格努斯數值研究[J].固體火箭技術,2015,38(4):465-471.

ZHAO Bo-bo,LIU Rong-zhong,GUO Rui.Numerical prediction of the Magnus effect for twist fin swept flight projectile[J].Journal of Solid Rocket Technology,2015,38(4):465-471.(in Chinese)

[9] 吳萍.二維彈道修正彈氣動布局分析[D].南京:南京理工大學,2014.

WU Ping.Aerodynamic configuration analysis of two-dimension trajectory correction projectile[D].Nanjing:Nanjing University of Science & Technology,2014.(in Chinese)

[10] 宋琦,楊樹興,徐勇,等.滾轉狀態下卷弧翼火箭彈氣動特性的數值模擬[J].固體火箭技術,2008,31(6):552-554.

SONG Qi,YANG Shu-xing,XU Yong.Numerical simulation on aerodynamic characteristic of rolling rocket with wrap around fins[J].Journal of Solid Rocket Technology,2008,31(6):552-554.(in Chinese)

[11] BHAGWANDIN V A.Numerical prediction of roll damping and magnus dynamic derivatives for finned projectiles at angle of attack[C]//30th AIAA Applied Aerodynamics Conference.New Orleans:AIAA,2012:1-17.

ResearchonAerodynamicCharacteristicsofMissileWithFinsUnderLow-speedRotation

ZHANG Chao,WANG Xue-de,WANG Yong,CAO You-quan

(School of Energy and Power Engineering,Nanjing University of Science and Technology,Nanjing 210094,China)

In order to research the effect of low-speed rotation on aerodynamic characteristics of missile with fins,the flow-field around the missile with fins under small attack-angle and total Mach number was simulated by applying 3D unsteady N-S equations and slide mesh technology.The Army-Navy Basic Finner Missile was used for verifying the validity of the numerical simulation method.The results show that this method has high accuracy.The results under different Mach numbers,rotation speed and roll angle show that the longitudinal aerodynamic characteristics(lift,drag,pitching moment)is independent of spin.The average rolling-moment-coefficient is contantly proportionate to rotation speed,and the average Magnus-force-coefficient exhibits a nonlinear change with the rotation speed.

missile;Magnus effect;sliding mesh;aerodynamic characteristics;numerical prediction

TJ303.4

A

1004-499X(2017)04-0040-08

2017-05-02

江蘇省青藍工程資助項目;國家留學基金委出國留學資助項目

張超(1991- )男,碩士研究生,研究方向為飛行器設計。E-mail:15996298561@163.com。

王學德(1977- ),男,副教授,研究方向為航空/航天飛行器外流場的數值模擬算法。E-mail:wangxuede2000@njust.edu.cn。

主站蜘蛛池模板: 伊人久久久久久久久久| 国产成人精品男人的天堂下载| 亚洲无码精彩视频在线观看| 亚洲综合二区| 另类综合视频| 久久影院一区二区h| 青青草原国产av福利网站| 孕妇高潮太爽了在线观看免费| 国产精品永久免费嫩草研究院| 国产一级在线观看www色| 在线观看亚洲人成网站| 欧美精品二区| 色综合久久久久8天国| 亚洲国产看片基地久久1024| 白丝美女办公室高潮喷水视频| 欧美在线一二区| 国产精品男人的天堂| 久久午夜夜伦鲁鲁片无码免费| 国产成人一区二区| 亚洲国产欧洲精品路线久久| 久久综合丝袜长腿丝袜| 成人一区在线| 香蕉伊思人视频| 亚洲国产日韩一区| 国产成人在线无码免费视频| 91九色最新地址| 亚洲人妖在线| 视频在线观看一区二区| 激情无码字幕综合| 色综合激情网| 亚洲国产精品人久久电影| 亚洲欧美人成人让影院| 欧美丝袜高跟鞋一区二区| www.99在线观看| 国产精品片在线观看手机版| 亚洲成人精品| 国产噜噜噜| 亚洲女同欧美在线| 久久久精品无码一区二区三区| 欧美乱妇高清无乱码免费| 久久久久国产一区二区| 国产精品极品美女自在线| 国产成人91精品| 欧美成人综合视频| 久久中文字幕2021精品| 五月婷婷激情四射| 精品福利国产| 免费一级无码在线网站| 亚洲国产高清精品线久久| 国产96在线 | 国产精品色婷婷在线观看| 欧美精品xx| 亚洲一级毛片在线播放| 欧美日一级片| 久久黄色免费电影| 亚洲不卡网| 91在线一9|永久视频在线| 欧美色香蕉| 亚洲三级色| 丁香婷婷激情网| yy6080理论大片一级久久| 国产欧美日韩专区发布| 久久99国产精品成人欧美| 国产大片黄在线观看| 少妇极品熟妇人妻专区视频| 国产日韩欧美视频| 成人午夜视频网站| 无码精油按摩潮喷在线播放| 久久77777| 精品一区二区无码av| 全部无卡免费的毛片在线看| 欧类av怡春院| 国产一区二区人大臿蕉香蕉| 国产成人免费高清AⅤ| 欧美a在线| 最新无码专区超级碰碰碰| 欧美日韩亚洲综合在线观看 | 久久国产精品波多野结衣| 无码啪啪精品天堂浪潮av| 国产xx在线观看| a级免费视频| 国产福利微拍精品一区二区|