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

長實驗時間對撞活塞壓縮風洞原理研究

2015-06-23 09:12:11龍鐵漢徐勝利
實驗流體力學 2015年2期
關鍵詞:實驗

龍鐵漢, 徐勝利

(1. 中國科學技術大學 近代力學系, 合肥 230027; 2. 清華大學 航天航空學院, 北京 100084)

長實驗時間對撞活塞壓縮風洞原理研究

龍鐵漢1, 徐勝利2,*

(1. 中國科學技術大學 近代力學系, 合肥 230027; 2. 清華大學 航天航空學院, 北京 100084)

對撞活塞壓縮風洞可獲得長時間的高壓高焓實驗氣流,是建造長時間運行高超聲速風洞的新思路。采用簡化分析和數值計算分析對撞活塞壓縮風洞的工作原理,以驗證對撞活塞壓縮風洞的可行性。簡化分析考慮定壓比熱隨溫度變化,基于流場均勻性假設,給出了實驗時間以及實驗氣體狀態參數隨時間變化歷程。采用商業軟件Fluent對全尺寸風洞流場運行過程進行數值模擬,模擬結果表明:擠壓式恒壓裝置可使風洞在較長時間內(約25ms)保持壓強近似不變,與簡化分析結果相符,但對撞活塞壓縮會使管道內流場溫度分布不均勻,導致溫度出現波動(相差不超過180K),這是對撞活塞壓縮風洞一個需要解決的問題。

活塞驅動器;風洞;長實驗時間;數值計算

0 引 言

活塞驅動器以近似絕熱壓縮的方式,可簡單有效地得到高總溫、高總壓氣體,對實驗氣體無污染,常被用于在激波風洞和膨脹管、二級輕氣炮中獲得高溫高壓驅動氣體,特別適合為M>8超燃發動機模擬實驗提供和飛行狀態匹配的來流總溫、總壓和氣體成分、Re數等,但存在實驗時間短的問題,典型實驗時間約1ms量級。當不計壁面摩擦,活塞驅動器相當于能量轉換器,將驅動活塞高壓氣體內能轉變為被壓縮氣體內能,基本原理是被壓縮氣體遵守能量守恒、理想氣體狀態方程和多方氣體關系。

活塞驅動器大致分為2類:(1)自由活塞。在壓縮終點不控制活塞回彈,用于在激波風洞和膨脹管、二級輕氣炮中獲得高溫高壓驅動氣體,例如已用于超燃實驗研究的自由活塞激波風洞T4[1-2]、HEG[3-4]、HIEST[5-6]和膨脹管風洞X3[7]。自由活塞在終點回彈,導致實驗氣體狀態參數急劇變化,實驗時間很短。(2)受控活塞。通過活塞壓縮直接得到所需高溫高壓實驗氣體,同時采用止退裝置防止活塞在壓縮終點回彈。早期采用圓柱剪切止退活塞的激波風洞,不能解決慣性力大(幾百到上千kN)和高溫燒蝕密封件問題。TsNIIMASH改進了止退方法,結合MCC(Multicascade Compression)方法,建造了長實驗時間(幾十到幾百ms)U系列風洞[8],其中PGU-11是完成了30msM=10超燃模型發動機自由射流實驗的活塞風洞[9-10]。MCC原理是:在駐室上游設置帶單向閥的多個不同構型腔體,通過活塞壓縮并止退在第1級腔室得到高壓高溫氣體,氣體在壓差驅動下依次流經各級腔體,產生漩渦導致總壓損失,但總溫和熵增大,即解決了被壓縮氣體總壓偏高、總溫偏低的問題,又可獲得長實驗時間較為穩定的實驗氣流。MCC方法非常復雜,需要精確設計腔室結構和腔室之間的閥門,獲得不同總溫總壓實驗氣體需改變腔體結構和閥門。另一種較為簡單的獲得長實驗時間的方法是采用擠壓裝置,此類風洞有ITAM研制的長時間A1和AT-303風洞[11],AT-303風洞最終未采用A1結構,是因為A1未解決單管活塞壓縮斜撞擊力對設備造成的損傷。總的看來,MCC方法和擠壓結構均可將活塞類風洞實驗時間延長至幾十到幾百ms。 本文提出新的對撞活塞壓縮風洞結構,在壓縮終點加裝止退裝置防止活塞回彈,并通過采用擠壓方式延長實驗時間。該結構可減小單活塞前向壓縮對設備產生的沖擊力,空氣經對撞活塞壓縮后再經恒壓活塞擠壓從噴管噴出,產生長時間的超聲速氣流。以模擬高超聲速飛行器發動機燃燒室的對撞活塞壓縮風洞為例,利用簡化理論分析和Fluent商業計算軟件,本文分析了駐室氣體流場參數分布和實驗時間,結合全尺寸風洞全過程流場計算結果,探索這類風洞運行原理,為后續風洞調試和改進提出建議。

1 簡化分析和結果討論

1.1 工作原理

為達到較長實驗時間,需解決2個問題:(1)破膜前,活塞到達壓縮終點受被壓縮氣體高壓作用回彈,反向加速度高達10000g。(2)破膜后,實驗氣體經噴管流出,駐室內實驗氣體總壓總溫下降。為此,本文提出的對撞活塞壓縮風洞專門設計了止退結構和擠壓式恒壓裝置(見圖1),對撞活塞在直線排列的壓縮管中相向運動,在壓縮終點被止退,恒壓管與壓縮管垂直。

1 壓縮管; 2 恒壓管; 3 過渡段; 4 噴管; 5 膜片; 6 止退裝置; 7 對撞活塞; 8 恒壓活塞

圖1 對撞活塞壓縮風洞結構示意圖

Fig.1 Schematic of shot tunnel driven by opposite compression piston drivers

如圖2所示,風洞運行過程可分如下階段:

(1)t0時刻擊發對撞活塞,對撞活塞在t1時刻運動至壓縮終點并被止退,p(t1)

(2)t2時刻擊發恒壓活塞,開始擠壓駐室氣體,有p(t2)=p(t1)

(3)t3時刻,p(t3)=pbreak,pbreak為破膜壓力。膜片破裂后噴管啟動。t4時刻有p(t4)=pm,忽略膜片破裂時間和噴管啟動時間,有t3≈t4。

(4) 恒壓活塞擠壓駐室氣體升壓,以補償氣體經噴管流出導致的壓降,p(t>t4)=pm。以速度v2(t>t4)=v2(t4)運動的恒壓活塞在t5時刻到達終點(膜片處),有效實驗時間結束。

圖2 對撞活塞壓縮風洞過程中氣體p-t曲線

Fig.2 Time history of pressure in the shot tunnel driven by opposite compression piston drivers

風洞運行過程中,通過2次活塞壓縮得到所需高溫高壓氣體。第1級為對撞活塞壓縮過程(t0→t1),對撞活塞在t1時刻被止退在壓縮終點,p(t1)

1.2 主要參數

對撞活塞壓縮管長l1=15m、內徑d1=300mm。恒壓活塞壓縮管長l2=0.92m、內徑d2=100mm。過渡段內徑100mm。根據圖紙尺寸計算對撞活塞壓縮體積V1和恒壓活塞壓縮體積V2分別為2.12m3和7.23×10-3m3,過渡段體積V3為1.89×10-3m3,被壓縮空氣初始體積2.13m3(Vi=V1+V2+V3)。對撞活塞驅動氣源壓強為pr1,恒壓活塞擠壓氣源壓強為pr2=pm,控制膜片厚度和刻槽深度使pbreak≈pm。

1.3 簡化分析

1.3.1 基本假設和計算參數

簡化理論分析假設氣體量熱完全,定壓比熱Cp為溫度T多項式函數[12],忽略破膜和噴管啟動時間,忽略管壁摩擦和傳熱,活塞速度較小(<60m/s),忽略壓縮過程中活塞上游產生的壓縮波系,認為駐室流場均勻且靜止,活塞壓縮為等熵過程。

針對M=7.5和飛行高度28km,當地大氣靜溫和靜壓為223K和1644Pa[13],飛行器來流速度和動壓為2238m/s和63.86kPa,總壓和總溫為16MPa和2400K。采用直連式實驗模擬燃燒室流場時,經進氣道斜激波壓縮,燃燒室進口M數降低、總溫不變但總壓減小。給定燃燒室進口氣流M數為2.9,靜溫和靜壓為1040K和68kPa,總溫和總壓為2400K和2.38MPa。

1.3.2 控制方程

在活塞壓縮過程中,空氣滿足如下能量守恒定律

(1)

其中:m為空氣質量,V為空氣體積,dV為對撞活塞和恒壓活塞壓縮導致的氣體體積變化,T為空氣溫度,dm和dT為空氣經噴管流出引起的質量和溫度變化,R為空氣氣體常數。

氣體體積變化滿足

(2)

其中:v1和v2分別為對撞活塞和恒壓活塞運動速度,N為壓縮管數。恒壓活塞速度滿足

(3)

(4)

其中:m1和m2分別為對撞活塞和恒壓活塞質量。

空氣質量變化滿足

(5)

其中:ρ*、v*和A*為噴管喉部空氣密度和速度、喉部截面積。

1.4 簡化分析結果與討論

1.4.1 實驗時間

當t>t4,駐室內空氣壓強和溫度保持不變,即dP=dT=0,方程(1)將化為

(6)

聯立方程(2)、(4)和(6)得t4時刻恒壓活塞速度為

(7)

t4時刻,恒壓活塞和噴管膜片的距離s(t4)為

(8)

其中:Vm為t4時刻空氣體積。實驗時間Δt為

(9)

pm=2.38MPa,Tm=2400K,由方程(7)得t4時刻v2(t4)為13.8m/s。t0時刻,有Ti=300K,Vi=2.13m3。由等熵關系得初始壓強pi=691Pa,vm=4.94×10-3m3。由方程(8)得s(t4)為0.39m,由方程(9)得實驗時間Δt=28ms。

1.4.2 理論最大實驗時間

圖1表明:恒壓活塞無法擠出過渡段(對撞活塞頂面之間)的空氣,其體積為38.3%Vm,空氣利用率為61.7%。在對撞活塞頭部加裝附加段(見圖3),可將過渡段空氣全部擠出,這對應著最大實驗時間Δtmax。令方程(9)的V3=0,得到Δtmax=45ms。

1.4.3 恒壓活塞質量

當t2

(10)

其中,ΔE為時間間隔(t2,t4)駐室空氣內能變化,有

(11)

其中:s(t4)為0.39m,T(t2)為2010K,T(t4)為2400K。聯立方程(10)和(11)得恒壓活塞質量m2為35kg。

圖3 對撞活塞頭部加裝附加段示意圖

1.4.4 初始溫度對實驗時間影響

由方程(9)知,Δt和(Vm-V3)成正比,Δtmax和Vm成正比。由等熵關系知,給定Tm和Vi,Vm隨Ti增大而增大,Δt也增大。當Ti為400K和500K,Vm增大為1.03×10-2m3和1.84×10-2m3,Δt分別為77ms和152ms,Δtmax分別為94ms和168ms。

1.4.5 狀態曲線

聯立方程(1)、(2)、(3)、(4)、(5)和理想氣體狀態方程,通過數值積分可給出p、T和v2隨t變化曲線(見圖4),p(t)和T(t)均出現平臺。當t>t4,恒壓活塞保持勻速運動,v2(t≥t4)=13.8m/s,調整m2可實現pm和Tm在噴管運行中保持恒定。

2 數值計算結果分析和討論

2.1 計算模型簡述

簡化理論分析基于駐室流場均勻性假設。不考慮壁面摩擦、傳熱和破膜延時,為認識流場非均勻性對壓強和溫度平臺的影響,采用商用計算軟件Fluent分析風洞運行流場演化過程。采用密度隱式求解器和有限體積法求解三維非定常N-S方程,通量分裂格式選用ROE格式,對流項離散格式選用二階迎風格式,選用標準k-ε湍流模型,湍動能k與湍流耗散率ε離散格式選用一階迎風格式,定容比熱采用關于溫度的分段多項式來計算,采用動網格模擬活塞運動。取對稱域一半為計算域(見圖5),初始網格數為352776,m2取為31kg。對撞活塞運動方向為x軸(壓縮管軸線)正向,恒壓活塞運動方向為y軸(噴管和恒壓管軸線)負向。指定點1坐標為(0,0,0),即壓縮管軸線和恒壓管軸線交點,指定點2坐標為(0,-0.295,0),靠近膜片處設置。

圖4 簡化分析得到的p-t、T-t和v2-t曲線

Fig.4 Time histories of pressure, temperature and balance piston velocity in simplified analysis

圖5 計算網格分布

2.2 計算結果分析與討論

t=507.7ms,對撞活塞速度55m/s,對撞活塞到達壓縮終點并被止退,圖6給出了對應的流場壓力、溫度和速度云圖。圖6表明:駐室內氣體受絕熱壓縮,壓強和溫度升高。壓強分布最均勻,溫度次之,速度分布最不均勻。原因是:對撞活塞壓縮駐室內空氣,恒壓管與過渡段相交處氣流折轉造成了速度分布非均勻性。同時,對撞活塞端面與管道壁面間形成狹縫,狹縫內氣體受強烈擠壓,存在明顯節流效應,其溫度和壓強遠高于駐室內其他部分氣體溫度和壓強。監控點1處壓強為1.01MPa、溫度為1998K,監控點2處壓強1.01MPa、溫度為1975K。

對撞活塞滯止后,圖7給出指定點1和2的p(t)、T(t)曲線。對撞活塞壓縮后的空氣進入過渡段和恒壓段并充分混合,指定點1和2壓強快速升高并達到平衡。盡管指定點溫度均快速上升,但溫度相差近似為常數。原因是:“十字型”管內氣體壓差靠壓力波平衡,特征時間為ms量級。溫差依賴熱傳導達到平衡,傳熱特征時間尺度大于風洞實驗時間。

當t=510ms,恒壓活塞開始向下壓縮空氣。圖8給出了不同時刻壓力、溫度和速度云圖。圖8(a)表明:當噴管膜片未破裂且恒壓活塞向下壓縮過程中,空氣壓力快速達到平衡,但溫度和速度非均勻性較明顯,未改變“十字架”空氣參數分布特征。當膜片處空氣壓強為24.6MPa(略高于pm)。當t=570.3ms,膜片破裂且噴管啟動。圖8(b)和8(c)表明:破膜后,向下運動的恒壓活塞掃過恒壓管和過渡段交叉口,空氣壓力均勻性較好,受恒壓活塞擠壓作用,過渡段與恒壓管相交部分的高溫氣體向噴管方向偏移,溫度非均勻性較明顯。當t=602.1ms,恒壓活塞到達膜片附近的壓縮終點。

圖9給出了不同時刻空氣溫度y軸分布。圖9表明:受恒壓活塞擠壓作用,沿y軸空氣溫度逐步上升,峰值溫度向噴管方向偏移,但溫度曲線形狀幾乎不隨時間改變,并整體向噴管方向偏移。這表明:“十字架”空氣溫度達到均勻分布所需特征時間較長,Tm在風洞運行過程中并不保持恒定,本算例溫度相差約180K。

圖10給出了數值計算和簡化理論分析結果對比。對撞活塞擠壓空氣射流進入“十字架”和膜片破裂時,指定點2 壓力和溫度曲線出現振蕩。對撞活塞壓縮后期存在明顯節流效應,對撞活塞對駐室空氣做功較等熵壓縮過程做功更多,壓力和溫度數值計算值大于簡化理論分析值。為得到恒定壓強曲線,數值計算的恒壓活塞質量(31kg)比簡化分析恒壓活塞質量(35kg)小。CFD計算結果表明:噴管啟動時間約5ms,總壓平臺延續時間約25ms,和簡化分析值28ms接近。過渡段射流進入恒壓段和駐室導致流場出現非均勻性,熱平衡所需時間大于風洞實驗時間,流場溫度分布呈非均勻性,T-t曲線出現凸起峰,峰值為2580K,比簡化分析值高約180K,凸起峰后的溫度平臺延續約10ms。

圖6 流場壓力、溫度和速度云圖(t=507.7ms)

Fig.6 Contours of pressure, temperature and velocity att=507.7ms

圖7 壓強和溫度隨時間變化曲線(指定點1和2)

(a) t=557ms

(b) t=574.1ms

(c) t=584.1ms

圖9 不同時刻溫度沿恒壓管軸線分布

圖10 CFD和簡化理論分析結果對比

Fig.10 Comparison of pressure and temperature time histories by CFD and simplified analysis

3 結 論

(1) 簡化分析和CFD數值模擬結果表明,對撞活塞壓縮風洞通過兩級壓縮得到所需總壓和總溫實驗氣體,止退結構防止對撞活塞回彈,控制恒壓活塞質量,恒壓活塞在恒壓氣源的驅動下將實驗氣體擠壓出噴管的工作模式可使實驗氣體總壓在較長時間內保持恒定,在本文算例中不小于25ms。活塞對撞壓縮風洞具有可行性。

(2) 簡化分析基于流場均勻性假設,可同時得到恒定的總壓和總溫,恒壓活塞質量應為31kg。考慮到流場的變化,CFD數值模擬得到恒定總壓時,恒壓活塞質量為35kg。由于對撞活塞壓縮后期在局部流場產生過熱空氣體射流,導致流場均勻性被破壞,壓強快速達到平衡之后,不同流場區域之間的溫差在風洞運行時間內保持穩定,總溫隨不同溫度的氣體被擠壓出噴管而波動,需研究提高溫度均勻性的方法。

(3) 簡化分析和CFD數值模擬均未考慮管壁與空氣之間熱量交換、壁面摩擦、膜片破裂過程以及破膜壓力較難精確控制對風洞運行狀態的影響,實現對撞活塞壓縮風洞仍是個挑戰。

[1] Morgan R G, Paull A, Stalker R J, et al. Shock tunnel studies scramjet phenomena[R]. NASA CR-181721, 1988.

[2] Stalker R J, Paull A, Mee D J, et al. Scramjets and shock tunnels—the Queensland experience[J]. Progress in Aerospace Sciences, 2005, 41(6): 471-513.

[3]Gardner A D. HyShot scramjet testing in the HEG[D].Queensland: School of Engineering, University of Queensland,2007.[4] Schramm J M, Karl S, Hannemann K, et al. Ground testing of the HyShot II scramjet configuration in HEG[R]. AIAA Paper-2008-2547, 2008.

[5]TANNO H, ITOH K, UEDA S, et al. Scramjet testing in high-enthalpy shock tunnel (HIEST)[C]//Proceedings of Symposium on Shock Waves, Japan, 2002: 12-16.

[6]Itoh K, Ueda S, Tanno H, et al. Hypersonic aerothermodynamic and scramjet research using high enthalpy shock tunnel[J]. Shock Waves, 2002, 12(2): 93-98.

[7]David G, Jorge S P, Morgan R G. High Mach number scramjet test flows in the X3 expansion tube[C]//Proceedings of the 29th International Symposium on Shock Waves, 2013.

[8]Anfimov N A, Kislykh V V. Multi-cascade compression-effective means to obtain high temperature dense gas in piston gas dynamic units (PGU)[C]//Current Topics in Shock Waves 17th International Symposium on Shock Waves and Shock Tubes. AIP Publishing, 1990, 208(1): 588-593.

[9]Orth R C, Kislykh V V. American and Russian hypersonic combustion research experiments in the TSNIIMASH PGU-11 facility[R]. AIAA Paper-95-0686,1996.

[10]Orth R C, Kislykh V V. Data analysis from hypersonic combustion tests in the TSNIIMASH PGU-11 facility[R]. AIAA paper-96-4584, 1996.

[11]Kharitonov A M, Zvegintsev V I, Fomin V M, et al. New generation hypersonic adiabatic compression with pressure multiplier[J]. Progress in Astronautics and Aeronautics, 2002, 198: 585-616.

[12]Rupley F M, Meeks E, Miller J A. CHEMKIN-III: A FORTRAN chemical kinetics package for the analysis of gas-phase chemical and plasma kinetics[M]. Livermore, CA: Sandia National Laboratories, 1996: 42-45.

[13] Noaa U S, Force U S A. US standard atmosphere, 1976[R]. Washington, DC: NOAA-S/T, 1976: 63-63.

(編輯:楊 娟)

Studies on principles of long-running shot tunnel driven by opposite compression piston drivers

Long Tiehan1, Xu Shengli2

(1. Department of Modern Mechanics, University of Science and Technology of China, Hefei 230027, China; 2. School of Aerospace Engineering, Tsinghua University, Beijing 100084, China)

The theoretical basis of the long-running shot tunnel driven by opposite compression pistons is laid. The long-running shot tunnel driven by opposite compression pistons can provide long-lasting high enthalpy and high pressure experimental gas flow, which is a new idea to build long-running hypersonic tunnels. The operation principles of the long-running shot tunnel driven by opposite compression pistons are analyzed by theoretical simplification and computational fluid dynamics to validate its feasibility. In theoretical analysis, based on the uniform flow field assumption, the test time and the time histories of pressure and temperature are presented with the consideration of the specific heat at constant pressure changing with temperature. The operation of the tunnel is simulated numerically by commercial code Fluent. The results indicate that the pressure of the experimental gas can be kept nearly constant with pressure stabilizer by extruding for about 25 microseconds, that is consistent with theoretical analysis results. Because of the compression of the opposite pistons, the temperature of the flow field is non-uniform and the temperature fluctuation amplitude is less than 180K. This is a notable problem of long-running shot tunnel driven by opposite compression pistons that needs to be resolved.

piston driver;wind tunnel;long test time;numerical simulation

1672-9897(2015)02-0090-08

10.11729/syltlx20140043

2014-04-11;

2014-10-10

國家自然科學基金面上項目資助(11172293)

LongTH,XuSL.Studiesonprinciplesoflong-runningShotTunnelDrivenbyOppositeCompressionPistonDrivers.JournalofExperimentsinFluidMechanics, 2015, 29(2): 90-96,102. 龍鐵漢, 徐勝利. 長實驗時間對撞活塞壓縮風洞原理研究. 實驗流體力學, 2015, 29(2): 90-96,102.

V211.753

A

龍鐵漢(1989-),男,湖南衡陽人,碩士研究生。研究方向:(高)超聲速設備及其實驗技術。通信地址:中國科學技術大學近代力學系(230027)。E-mail: longtiehan@sina.com

*通信作者 E-mail: slxu@mail.tsinghua.edu.cn

猜你喜歡
實驗
我做了一項小實驗
記住“三個字”,寫好小實驗
我做了一項小實驗
我做了一項小實驗
記一次有趣的實驗
有趣的實驗
小主人報(2022年4期)2022-08-09 08:52:06
微型實驗里看“燃燒”
做個怪怪長實驗
NO與NO2相互轉化實驗的改進
實踐十號上的19項實驗
太空探索(2016年5期)2016-07-12 15:17:55
主站蜘蛛池模板: 久久精品这里只有精99品| 九色综合伊人久久富二代| 亚洲成人福利网站| 国内精自视频品线一二区| 国产欧美成人不卡视频| 午夜无码一区二区三区| AV无码无在线观看免费| 亚洲AV永久无码精品古装片| 国产福利微拍精品一区二区| 一区二区在线视频免费观看| 91青青在线视频| 亚洲人成影院在线观看| 美女一区二区在线观看| 国产毛片高清一级国语| 中文字幕无码av专区久久| 亚洲精品无码抽插日韩| 亚洲国产成熟视频在线多多 | www.精品视频| 国产91九色在线播放| 欧美成人日韩| 四虎永久在线精品影院| 国产正在播放| 日韩成人午夜| 成人午夜视频网站| 国产亚洲视频免费播放| 日韩最新中文字幕| 国产福利观看| 精品久久国产综合精麻豆| 99久久精彩视频| 久草中文网| 久久综合亚洲色一区二区三区| 无码网站免费观看| 欧美综合中文字幕久久| 在线观看无码av五月花| 爽爽影院十八禁在线观看| 伊人成人在线| 色婷婷综合在线| 亚洲日韩Av中文字幕无码| 2021国产精品自拍| 亚洲欧美成人网| 亚洲人成网站日本片| 色综合五月| 久久精品欧美一区二区| 国产精品午夜福利麻豆| 无码AV日韩一二三区| 精品久久久久久成人AV| 丰满人妻久久中文字幕| 欧美成人午夜影院| 亚洲无码视频一区二区三区| 97视频免费在线观看| 天天躁日日躁狠狠躁中文字幕| 国产在线一区视频| 91香蕉国产亚洲一二三区| 欧美日本激情| 天天爽免费视频| 热这里只有精品国产热门精品| 无码内射中文字幕岛国片| 久无码久无码av无码| 国产精品一线天| 国产国产人免费视频成18| 亚洲一级色| 国产香蕉国产精品偷在线观看| 天天躁狠狠躁| 久草性视频| 久久毛片网| 美女无遮挡被啪啪到高潮免费| 最新精品国偷自产在线| 五月婷婷亚洲综合| 97精品久久久大香线焦| 国产精品密蕾丝视频| 美女高潮全身流白浆福利区| 在线观看热码亚洲av每日更新| 亚洲日本中文字幕乱码中文 | 欧美亚洲网| 九色视频在线免费观看| www.精品国产| 国产无遮挡猛进猛出免费软件| 精品一区二区三区波多野结衣| 手机精品福利在线观看| 久久国产高清视频| 999国内精品久久免费视频| 日本欧美中文字幕精品亚洲|