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

耦合尾噴管堵蓋運動的發射箱內流場研究

2014-06-27 05:41:55于邵禎姜毅周笑飛牛鈺森孫璐璐
兵工學報 2014年11期
關鍵詞:區域實驗

于邵禎,姜毅,周笑飛,牛鈺森,孫璐璐

(1.北京理工大學宇航學院,北京 100081;2.海軍航空工程學院青島校區,山東青島 266041)

耦合尾噴管堵蓋運動的發射箱內流場研究

于邵禎1,姜毅1,周笑飛1,牛鈺森1,孫璐璐2

(1.北京理工大學宇航學院,北京 100081;2.海軍航空工程學院青島校區,山東青島 266041)

利用初始沖擊波超壓完成前后蓋開啟過程的貯運發射箱已得到廣泛應用。為研究含尾噴管堵蓋的沖擊波超壓形成過程及對后易碎蓋的作用效果,應用有限元方法并結合動網格技術建立了導彈點火后堵蓋的運動模型,并通過實驗方法對仿真結果進行了驗證。結合計算結果可清晰地看到尾焰流場的形成過程,并得到了沖擊波超壓在后易碎蓋表面的隨時間變化曲線。研究表明:受堵蓋的影響,沖擊波超壓首先形成并沖擊后易碎蓋,燃氣由堵蓋的邊緣向中心匯聚形成主流,在對后易碎蓋的沖擊時間和作用位置上與沖擊波作用有明顯的不同;后易碎蓋主要受到沖擊波超壓作用實現碎裂變形,在堵蓋運動的投影區域首先達到最大受力,瞬時峰值達5×105Pa.

兵器科學與技術;尾噴管堵蓋;沖擊波;易碎蓋;數值仿真

0 引言

沖擊波開蓋技術在目前應用廣泛并日臻成熟。發射箱內沖擊波的形成是導彈在箱內點火后,高壓燃氣流從噴管噴出與周圍空氣形成最初的壓力界面,并以同射流邊界大致近似的形狀隨著燃氣流向外發展而向前推進,同時伴隨能量不斷增強[1]。沖擊波開蓋技術是利用沖擊波能量實現對設計有一定承壓標準的易碎蓋的開啟[2]。利用沖擊波開蓋技術的貯運發射箱結構簡單,質量輕,操作維護方便,經濟成本較低,因此有廣闊的應用前景。

箱體前后易碎蓋的壓力匹配問題是采用易碎蓋技術亟待解決的問題之一,前后易碎蓋的強度設計標準與沖擊波的強度密切相關。沖擊波在形成過程中受到箱內配件如導軌、擾流器、堵蓋等影響。在諸多影響因素中,尾噴管堵蓋對箱內射流形狀和沖擊波的形成及強度影響較大,而對于某些型號的導彈發動機來說,尾噴管堵蓋安放有點火器并起到射前增壓的作用[3],為導彈必不可少配件,因此有必要考慮有堵蓋影響的發射箱內沖擊波的形成過程。

對于易碎蓋技術的研究,國外起步較早,并且已經有成熟的應用,如美國的“陸麻雀”、“愛國者”“戰斧”,意大利的“阿斯派德”、“信天翁”等,同時對于超音速流產生的沖擊波在能量、成因以及數值模擬方法和實驗的研究上取得了很大的進展[4-6]。國內對于箱內沖擊波的研究主要集中在沖擊波形成后在箱內的傳播過程,分析過程以實驗并輔助工程經驗為主[7-11],也取得了一定的研究成果,對于沖擊波在發射箱內部的作用效果也有了一定的認識,但是目前對于利用沖擊波能量完成開蓋動作的發射箱設計主要以工程經驗為主,對于沖擊波的強度分析只能預估在一定的范圍內,不能進行定量分析,嚴重影響了易碎蓋的設計精度和使用的可靠性。

本文通過數值仿真計算模擬了尾噴管堵蓋在箱體內受到燃氣流作用的運動過程,并通過監測發射箱后易碎蓋的平均靜壓研究沖擊波作用力的變化過程,同時根據實驗結果對沖擊波作用效果進行驗證,為貯運發射箱易碎蓋設計提供參考。

1 基本控制方程

燃氣射流的控制方程組采用三維非定常方程組[12]。湍流方程采用標準k-ε模型方程(具體方程見文獻[12])。

湍動能k的輸運方程為

式中:μt為湍流粘性系數;Gk為由層流速度梯度產生的湍流動能;Gb為由浮力產生的湍流動能;C1ε、C2ε、C3ε、σk、σε為經驗常數;Cμ為湍流常數。

在動網格計算區域內,在任意控制體V中任意標量的積分形式的控制方程可以表示為

式中:ρ表示流體密度;φ為通量變量;u表示流體速度矢量;ug表示動網格運動速度;Γ表示耗散系數; Sφ表示源項。

2 計算實例

2.1 計算模型

仿真模型模擬發射箱內導彈點火時刻尾噴管堵蓋從初始位置運動到后易碎蓋端面過程。計算模型包括發射箱箱體、前后導軌、導彈彈體、發動機噴管、堵蓋在內的三維模型。模型中根據堵蓋的設計外形將堵蓋簡化處理為圓形薄片,并將原堵蓋質量均勻加載。計算區域如圖2所示。由于堵蓋運動行程短,速度大,同時導彈的發射姿態以傾斜和垂直發射為主,因此將堵蓋的運動軌跡按照直線形式處理。模型中取堵蓋面向發射箱前易碎蓋端面為前端面,反向為后端面。計算模型根據對稱性采用1/2模型計算。

模型網格劃分如圖2所示,采用6面體結構化網格,對于噴管區域采用加密處理,同時設置堵蓋運動區域的邊界面為Interface邊界,在Interface邊界的兩側采用等比例網格尺寸劃分,其網格高度為0.5 mm,比例系數為1.2,如圖2所示。網格總數為40萬。

圖1 仿真計算物理模型Fig.1 The physical model for simulation calculation

圖2 計算網格模型Fig.2 Mesh model

2.2 計算條件

如圖3所示,對實驗所測得的發動機壓力變化曲線進行數值離散處理并擬合出發動機燃燒室內壓強變化方程,假設總溫變化在點火初始狀態至發動機工作穩定時刻為線性變化,將總溫與總壓變化曲線編寫自定義函數輸入到仿真計算中。

計算使用Fluent軟件,應用有限體積法。對彈體、導軌、噴管及發射箱等固壁表面采用標準壁面函數方程處理。堵蓋的運動速度為自定義函數:通過讀取堵蓋兩側的壓強分布并積分轉換成作用力加載在堵蓋表面進行計算。根據實驗測得的破膜壓力為1.5 MPa,因此仿真計算中將燃燒室壓力達到1.5 MPa時刻設置為堵蓋運動初始時刻。

圖3 壓力入口邊界條件Fig.3 Boundary conditions on the pressure inlet

2.3 數值計算驗證

模擬實驗可有效檢驗數值計算結果的精度和可靠性,但是在實驗中不能直接對堵蓋的運動狀態和受到的作用力進行測量。而發射箱后蓋采用易碎蓋設計,同樣不能直接獲取后蓋的受力情況。因此只能通過對其他相關數據進行測量對比進行間接驗證。發動機點火后在發射箱內形成的初始沖擊波超壓是與燃氣流特性最相關的物理量,利用壓力傳感器可直接對超壓值進行測量,如圖1所示的發射箱模型,實驗中在發射箱壁面上設置監測點1、2兩個壓力監測點對發射箱內沖擊波數據采集。

圖4所示為發動機點火后發射箱內沖擊波的傳播過程。白色區域為沖擊波超壓值影響區域,從圖4中可以看出發射箱內初始沖擊波超壓值接近170 000 Pa,并且波前和波后超壓值有明顯的階躍,因此在實驗中能夠對沖擊波流經監測點處的作用時長和峰值進行監測并獲得相關數據。

圖4 初始沖擊波傳播過程Fig.4 The propagation of initial shock wave

圖5為沖擊波流經監測點處的超壓變化曲線,需要說明的是實驗過程中前易碎蓋的開啟是利用多組沖擊波的積聚達到一定強度后實現。因此對于初始沖擊波在監測點處的傳播歷程由于反射作用會出現多次峰值,實驗曲線如圖5(a)所示。在監測點1、2中,監測點1的初始峰值出現時間要早于監測點2,而第2個峰值出現時間在監測點2第2個峰值之后,有力的證明了沖擊波的反射結論。而仿真計算只記錄了發射箱首次傳播過程,因此在進行數據對比中只對沖擊波第1個峰值大小和時程進行對比。另外由于時間基準在實驗和數值計算中很難統一,在圖5(b)和圖5(c)的壓強變化曲線中選取峰值點時刻作為基準進行對比分析。由監測點的壓強變化曲線可知實測值與計算值壓強變化在時間跨度和峰值上都比較接近。從表1中數據對比可知,監測點理論計算強度和實測值誤差最大為4.51%.在時間歷程上二者誤差最大為33.30%.而堵蓋與燃氣流耦合作用主要受強度因素的影響,因此數值計算滿足精度要求。

表1 監測點壓強數據對比Tab.1 Comparison of pressure data

3 計算結果分析

3.1 仿真云圖分析

圖6所示為堵蓋運動0.5 ms時刻燃氣流從堵蓋與噴管間隙流出的仿真云圖。圖6(a)所示為對稱面燃氣流靜壓云圖,可以清晰地看到堵蓋周圍的燃氣壓強緊貼壁面處發生變化,燃氣流從噴管與堵蓋間隙流出后由邊界向堵蓋中心擴散。圖6(b)所示為噴管內溫度云圖,燃氣流開始沿管壁向外擴散,在此狀態下,燃氣流溫度分布與壓強分布一致。發射箱易碎蓋上暫未受到燃氣流擾動的影響。

圖5 監測點壓強變化曲線Fig.5 Curves of pressures at monitoring points

圖7所示為堵蓋運動1.5 ms后發射箱內仿真結果。圖7(a)所示為箱內對稱面上靜壓云圖,在堵蓋前端面受燃氣流的作用存在一個高壓強集中區,出現原因在于燃氣流對堵蓋的沖擊作用使堵蓋端面上存在速度滯止區,在相應區域壓強值較高。同時在發射箱后易碎蓋端面的中心區域也存在高壓區域。圖7(b)所示為燃氣流溫度云圖,燃氣流沿著噴管壁向外流動,主流并未作用到易碎蓋端面中心上。圖7(c)所示為后易碎蓋端面壓強分布,后易碎蓋端面的壓強主要集中在堵蓋的投影區域范圍內,達到2.7×105Pa,易碎蓋上整體平均壓強略有升高,最低值達到1.42×105Pa.圖7(d)為易碎蓋端面上溫度分布云圖,溫度分布由中心向外逐漸減弱,但總體溫升較低,證明燃氣流主流暫未接觸易碎蓋,同時在易碎蓋上出現一個高溫的環形區域,溫度達到400 K.

圖6 堵蓋運動0.5 ms后噴管內燃氣流云圖Fig.6 Contours of jet flow after 0.5 ms

由圖7分析結果可知:受堵蓋影響,燃氣流與初始沖擊波作用位置不同,并且沖擊波形成于燃氣流前端,其形狀與燃氣流邊界形狀也有所不同。

圖8所示為堵蓋運動3 ms飛出噴管后仿真結果。圖8(a)所示在堵蓋前端面上壓強較低,證明燃氣流在噴管內充分膨脹。而在堵蓋后部靜壓值較高,接近3×105Pa.圖8(b)所示燃氣流溫度云圖,此刻燃氣流已經沖擊到發射箱易碎蓋上并反向向前流動,發射箱后部空間被燃氣流所填充,溫度趨于一致接近1 900 K,而在堵蓋前后端面被燃氣流包圍的包覆層,溫度略高,在2 600 K左右。圖8(c)所示為發射箱后易碎蓋端面上壓強云圖,如圖8(c)所示,在燃氣流作下,在堵蓋投影的中心區域內壓強最大,達到4.67×105Pa,并向外逐漸減小。圖8(d)所示為發射箱后端面上溫度云圖,從圖中可以看出易碎蓋上大部分區域溫度約為2 200 K,局部有小范圍的高溫區域,接近2 400 K.在易碎蓋下半部分燃氣流溫度略低。造成這種現象的原因在于底部的導軌阻礙燃氣流的擴散,在該部位燃氣流產生壅塞現象。

圖7 1.5 ms流場仿真結果Fig.7 Contours of jet flows at 1.5 ms

3.2 仿真數據分析

為具體說明燃氣流在堵蓋運動過程中的作用效果,通過讀取堵蓋的速度值、堵蓋兩側的壓強值說明燃氣流與堵蓋的相互作用,同時讀取易碎蓋端面上壓強和受力,得到相應參數變化曲線,將后易碎蓋按照堵蓋的幾何投影劃分為堵蓋投影區域和其他區域兩部分進行對比分析,具體結果如圖9~圖13所示。

圖9所示為堵蓋速度變化曲線,從圖中可知堵蓋運動3.6 ms后撞擊后易碎蓋,在前1.5 ms運動過程中,燃氣流對堵蓋的作用基本穩定,堵蓋作勻加速運動,1.5 ms后堵蓋運動速度基本趨于平衡。

圖10所示為堵蓋前后端面靜壓變化曲線,從圖中可知堵蓋前端面壓強在堵蓋初始運動時經過小范圍的波動后逐漸降低并在末尾時刻降低到1×105Pa,后端面壓強變化平緩,基本在2×105Pa附近波動,從初始環境壓強逐漸升高并伴有小幅的波動發生,在堵蓋運動2.2 ms時刻,前后端面壓強相等,隨后后端面壓強高于前端面,堵蓋受力反向作減速運動直至3.5 ms時刻撞擊至后易碎蓋。

圖8 3 ms流場仿真結果Fig.8 Contours of jet flows at 3 ms

圖9 堵蓋速度曲線Fig.9 Velocity curve of nozzle closure

圖10 堵蓋前后端面壓強曲線Fig.10 Pressures on the surface of nozzle closure

圖11 后易碎蓋端面平均壓強曲線Fig.11 Average pressures on the surface of back cover

圖12 后易碎蓋端面最大壓強曲線Fig.12 Maximum pressures on the surface of back cover

為說明后蓋上壓強分布規律,圖11和圖12所示壓強曲線中分別給出了后蓋上在堵蓋投影區域和除投影區域外的平均靜壓和最大靜壓變化曲線以及整個后蓋平均壓強曲線,從圖11所示的壓強變化曲線中看出,在堵蓋投影區域出現峰值的1.5 ms時刻,后蓋上其他區域并未有明顯的波動,而整個后蓋端面上峰值出現時刻在2.5 ms,落后于該點,由此說明在投影區域所出現的峰值由沖擊波形成,后蓋整體的壓強峰值由燃氣流造成。同時可以看出在后易碎蓋端面上堵蓋投影區域壓強較大,首先受到破壞。

從圖12的最大壓強變化曲線上分析,在后蓋端面上除投影區域外壓強最大值點出現兩次峰值,并且第1次峰值出現時刻與堵蓋投影區域峰值出現時間一致,說明兩個表面峰值點由同一因素作用所致。而第2次峰值點的出現時刻堵蓋投影區域并未出現明顯波動,因此可以斷定為燃氣流作用,更有力地證明了沖擊波與燃氣流作用的時間先后性以及作用位置的差別。并且從沖擊波峰值曲線上可以判斷初始沖擊波的作用范圍和強度。

圖13為堵蓋和易碎蓋上的受力曲線,對比壓強變化曲線可知,二者的受力變化曲線與壓強變化基本一致,從易碎蓋受力曲線上可以判斷在1.5 ms沖擊波達到峰值時刻由于其在易碎蓋上作用區域集中并且相對面積較小,因此對于整個箱蓋的受力影響不明顯。圖中箱蓋受力曲線可以為易碎蓋的設計強度標準提供參考值。

圖13 堵蓋和易碎蓋上作用力曲線Fig.13 Forces on nozzle closure and back cover

4 綜合實驗驗證

實驗針對某型導彈點火過程利用高速攝影捕捉到發射箱后易碎蓋的破碎過程,由于發射箱采用易碎蓋設計,并且在導彈點火后易碎蓋所處環境惡劣,無法直接在易碎蓋上布置傳感器,因此將監測點3放置于箱壁內壁面距離后易碎蓋端面頂點0.55 m位置,如圖1所示。實驗中后易碎蓋的設計開蓋壓力為0.06 MPa.由于在仿真過程中假設發射箱易碎蓋達到開啟壓力后仍然封閉,而在實驗過程中,后易碎蓋被燃氣流擊碎。而在易碎蓋破碎前,燃氣流流動已經轉向,會在監測點處出現第1個沖擊波峰值,該峰值大小的仿真與實驗值具有一定的可比性,因此只針對燃氣流流經監測點的壓力峰值對比,仿真實驗數據峰值點絕對靜壓值為150 471.14 Pa,換算為表壓值為0.049 MPa,實驗過程中實測表壓值0.055 MPa,仿真計算誤差為3.14%,因此可以判斷仿真結果存在較高的可信度。仿真數據曲線如圖14所示。另外從圖15的高速攝影圖像中可以清晰看到在發射箱后易碎蓋上首先是在箱蓋中心出現與堵蓋形狀相同的圓形變形區域并破裂,隨后從箱蓋與箱體的結合部位有燃氣溢出,因此實驗過程驗證了仿真結果的正確性。

圖14 監測點壓力變化曲線Fig.14 Curves of static pressure at the monitoring point

圖15 高速攝影實驗結果Fig.15 High-speed photograph of experiment

5 結論

本文通過仿真與實驗相結合的方法,分析了在尾噴管堵蓋影響下的發射箱內燃氣流場的作用以及初始沖擊波的形成特點,獲得了以下具體結論:

1)在有堵蓋的燃氣射流中,燃氣流通過堵蓋與噴管的縫隙流出并向噴管四周擴散,同時與堵蓋后方滯止空氣接觸并形成初始沖擊波陣面向外傳播。

2)對于發射箱蓋的作用沖擊波先于燃氣射流并且主要集中在堵蓋的投影區域上,燃氣流的沖擊作用在易碎蓋上首先形成圍繞中心的環形區域,隨后向四周擴散。

3)發射箱后蓋同時受到沖擊波和燃氣射流的作用,沖擊波的強度較大,作用區域集中,燃氣流的沖擊強度平均值略小,作用在整個后蓋上,是發射箱后蓋開啟或破碎的主要因素。

4)堵蓋在后蓋的投影區域為整個箱蓋上受力最嚴重的區域,首先被破壞,因此對于設計后易碎蓋的破裂強度時應在此位置進行加強處理以避免在該區域首先被破壞的狀態下達不到整體開蓋壓力。

References)

[1] 柯朗R,弗里德里克斯K O.超聲速流與沖擊波[M].李維新,徐華生,管楚洤,譯.北京:科學出版社,1986.

Courant R,Friedrichs K O.Supersonic flow and shock waves[M]. LI Wei-xin,XU Hua-sheng,GUAN Chu-quan,translated.Beijing:Science Press,1986.(in Chinese)

[2] 苗佩云,袁曾鳳.同心發射筒燃氣開蓋技術[J].北京理工大學學報,2004,24(4):283-285.

MIAO Pei-yun,YUAN Zeng-feng.Techniques for the automatic cover opening in concentric canister launcher[J].Transactions of Beijing Institute of Technology,2004,24(4):283-285.(in Chinese)

[3] JI H拉夫洛夫.固體火箭發動機結構[M].關正西,譯.北京:中國宇航出版社,2006.

JI H Lavrov.Structure of solid rocket engine[M].Guan Zhengxi,translated.Beijing:China Astronautic Publishing House, 2006.(in Chinese)

[4] Arnab C,Abdellah H.Numerical simulation of transient supersonic nozzle flows[C]∥50th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition.Nashville,TN:AIAA,2012.

[5] Steven S,Jacob D,James V.Experimental performance analysis of a toroidal aerospike nozzle integrated with a N2O/HTPB hybrid rocket motor[C]∥46th AIAA/ASME/SAE/ASEE Joint Propulsion Conference&Exhibit.Nashville,TN:AIAA,2010.

[6] Binu P,Manhar D.Surface pressure fluctuations due to an impinging supersonic under expanded jet[C]∥48th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition.Nashville,TN:AIAA.2010.

[7] 傅德彬,姜毅.某導彈易碎蓋的開啟過程[J].固體火箭技術, 2007,30(4):275-277.

FU De-bin,JIANG Yi.Opening process of friable lid of one missile[J].Journal of Solid Rocket Technology,2007,30(4):275-277.(in Chinese)

[8] 李學民.導彈發射箱(筒)內燃氣流特性分析—熱態試驗研究[J].宇航學報,1997,18(4):67-74.

LI Xue-min.Characteristic analysis of missile combustion-gas in a launching-container(launching-tube)—hot experiment study[J]. Journal of Astronautics,1997,18(4):67-74.(in Chinese)

[9] 劉琦,傅德彬,姜毅.貯運發射箱內燃氣射流的非定常沖擊波流場數值模擬[J].彈箭與制導學報,2005,25(S4):382-384.

LIU Qi,FU De-bin,JIANG Yi.Unsteady simulation of shock wave in launcher[J].Journal of Projectiles,Rockets and Guidance,2005,25(S4):382-384.(in Chinese)

[10] 徐強,李開明,張福祥,等.封閉式導彈發射箱內燃氣流特性實驗研究[J].彈道學報,1995,7(2):52-56.

XU Qiang,LI Kai-ming,ZHANG Fu-xiang,et al.Experimental study of the internal exhaust flow field in closed rectangular missile launch tube[J].Journal of Ballistics,1995,7(2):52-56. (in Chinese)

[11] 徐強,李軍.燃氣射流起始沖擊波形成機理的實驗研究[J].推進技術,2000,21(3):16-18.

XU Qiang,LI Jun.Experimental study on mechanism of initial shock wave in jet flow[J].Journal of Propulsion Technology, 2000,21(3):16-18.(in Chinese)

[12] 趙承慶,姜毅.氣體射流動力學[M].北京:北京理工大學出版社,1998.

ZHAO Cheng-qing,JIANG Yi.The kinetics of gas jet[M]. Beijing:Beijing Institute of Technology Press,1998.(in Chinese)

Research on Distribution of Flow Field in Launching Canister with the Effect of Nozzle Closure

YU Shao-zhen1,JIANG Yi1,ZHOU Xiao-fei1,NIU Yu-sen1,SUN Lu-lu2
(1.School of Aerospace Engineering,Beijing Institute of Technology,Beijing 100081,China;
2.Naval Aeronautical Engineering Academy Qingdao Branch,Qingdao 266041,Shandong,China)

The overpressure of shock wave is widely used in opening the launching canister.In order to study the formation of shock wave and its effect on the fragile lid with a nozzle closure,a simulation model is established by using a numerical method,in which a dynamic mesh technique is used to update the meshes of flow field.The simulation results are verified with the experimental results.The formation process of jet flow is shown clearly and the changing curves of the overpressure on the back lid with time are obtained from the analysis results.The research results show that the overpressure of shock wave forms on the nozzle closure and then impact the fragile lid.The gas converges from the edge of nozzle closure to the center to form a mainstream.The effects of the mainstream and the shock wave on the back fragile lid are obviously different in impact time and position.The back fragile lid is brittlely deformed under the action of overpressure of shock wave.In the projection area of the nozzle closure,the instantaneous peak first reaches to the maximum value of 5×105Pa.

ordnance science and technology;nozzle closure;shock wave;fragile lid;numerical simulation

TJ760

A

1000-1093(2014)11-1805-08

10.3969/j.issn.1000-1093.2014.11.011

2013-12-13

于邵禎(1985—),男,博士研究生。E-mail:bitysz@bit.edu.cn;

姜毅(1965—),男,教授,博士生導師。E-mail:jy2818@163.com

猜你喜歡
區域實驗
記一次有趣的實驗
微型實驗里看“燃燒”
永久基本農田集中區域“禁廢”
今日農業(2021年9期)2021-11-26 07:41:24
分割區域
做個怪怪長實驗
NO與NO2相互轉化實驗的改進
實踐十號上的19項實驗
太空探索(2016年5期)2016-07-12 15:17:55
關于四色猜想
分區域
基于嚴重區域的多PCC點暫降頻次估計
電測與儀表(2015年5期)2015-04-09 11:30:52
主站蜘蛛池模板: 国产黄色免费看| 成年A级毛片| 国产激爽大片在线播放| 久久国产免费观看| 精品久久国产综合精麻豆| 亚洲婷婷在线视频| 欧美亚洲国产精品第一页| 久久这里只有精品8| 国产特级毛片aaaaaa| 毛片网站免费在线观看| 国产又大又粗又猛又爽的视频| 免费播放毛片| 亚洲国产精品无码久久一线| 丁香婷婷综合激情| 中文字幕第4页| 亚洲无码91视频| 日韩精品一区二区深田咏美| 少妇被粗大的猛烈进出免费视频| 午夜福利网址| 人妻一区二区三区无码精品一区| 无码福利日韩神码福利片| 国产熟女一级毛片| 国产成人一区免费观看| 91在线一9|永久视频在线| 国模粉嫩小泬视频在线观看 | 国产欧美在线| 欧美h在线观看| 国产人碰人摸人爱免费视频| av在线人妻熟妇| 99久久精品免费视频| 熟女日韩精品2区| 99资源在线| 亚洲精品欧美日韩在线| 成人久久精品一区二区三区| 国产噜噜噜| 日本一区二区三区精品国产| a毛片在线播放| 国产亚洲欧美在线中文bt天堂| 这里只有精品在线播放| 亚洲人成网7777777国产| 国产人成在线视频| 在线国产综合一区二区三区| 一级全免费视频播放| 色综合久久无码网| 性做久久久久久久免费看| 免费jjzz在在线播放国产| 国产呦视频免费视频在线观看| 国产免费怡红院视频| 国产免费久久精品99re丫丫一| 操国产美女| 亚洲精品777| 国产综合欧美| 真实国产精品vr专区| 5555国产在线观看| 一级毛片免费的| 五月六月伊人狠狠丁香网| 国产高清在线丝袜精品一区| 青青草a国产免费观看| 日日碰狠狠添天天爽| 伊人大杳蕉中文无码| 精品久久蜜桃| 一级做a爰片久久免费| 曰韩人妻一区二区三区| 欧美一道本| 成人日韩欧美| a毛片免费在线观看| 伊人网址在线| 少妇人妻无码首页| 国产乱子伦精品视频| 婷婷色婷婷| 99视频免费观看| 亚洲高清在线播放| 精品视频福利| 国产无吗一区二区三区在线欢| 精品久久人人爽人人玩人人妻| 色爽网免费视频| 久久久精品无码一区二区三区| 在线观看国产精品第一区免费| а∨天堂一区中文字幕| 亚洲伦理一区二区| 成人福利在线视频免费观看| 欧美区一区|