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

相近韋伯數(shù)條件下激波波后液滴初期變形的影響機制*

2018-05-21 11:01:49易翔宇朱雨建楊基明
爆炸與沖擊 2018年3期
關鍵詞:變形實驗

易翔宇,朱雨建,楊基明

(1.中國科學技術大學近代力學系,安徽 合肥 230027;2.中國航天空氣動力技術研究院,北京 100074)

在激波誘導的高速氣流中,液滴會在氣流氣動力作用下發(fā)生變形和破碎;這是高速兩相流研究中的一個典型問題。同時液滴的氣動破碎廣泛存在于各類工程實踐中,如發(fā)動機液態(tài)燃料噴射霧化中的二次破碎[1-2]和飛行器表面的雨滴侵蝕防護等。因此,對這一現(xiàn)象的研究具有重要學術價值和工程意義。

高韋伯數(shù)條件下,液滴迅速變形并被液霧籠罩,破碎初期未被液霧遮擋的液滴變形圖像是研究液滴破碎機理的重要依據(jù)。另一方面,液滴變形初期的形態(tài)直接決定了完成破碎時液滴的拓撲結構,對破碎效率有著不可忽視的影響。在SIE機制下的液滴破碎實驗中,液滴整體形狀的演化過程并不單一,而是呈現(xiàn)出多種不同的模態(tài)。一般認為,在完全破碎前,液滴主要呈現(xiàn)出一種扁平的橢球形或碟形[8-9];Wierzba等[10]以激光干涉法捕捉液滴變形流場則揭示出更復雜的液滴變形形態(tài):在液滴的背風面出現(xiàn)一個環(huán)形的脊狀突起,并在其頂端有液霧生成。在此后的相關研究中,類似的背風面突起形態(tài)被多次觀測和描述,而隨實驗條件的不同,這些突起形態(tài)也有所差異[5]。到目前為止,對液滴破碎初期變形的研究多集中于低韋伯數(shù)條件[11]或高韋伯數(shù)條件下液滴表面的不穩(wěn)定性現(xiàn)象[7,12],對于高韋伯數(shù)條件下液滴初期變形不同模態(tài)出現(xiàn)的原因以及揭示液滴初期變形對液滴破碎效率影響的研究工作相對較少。

本文中通過實驗,對韋伯數(shù)在2 100~2 700區(qū)間內(nèi)液滴破碎的初期變形進行觀測,獲得在不同實驗條件下具有顯著差異的液滴變形圖像。在已有工作[13]的基礎上,利用數(shù)值模擬結合理論分析,解析液滴周圍流場的發(fā)展過程,探討不同變形模態(tài)的形成原因。

1 研究方法

1.1 實驗設備和方法

實驗在設有電控破膜系統(tǒng)的矩形截面(40 mm×70 mm)水平激波管中進行[13]。通過激光束和光電二極管探測下落的液滴,并經(jīng)同步控制系統(tǒng)觸發(fā)電控破膜裝置產(chǎn)生激波,保證液滴與激波在觀察窗預設區(qū)域內(nèi)相互作用并發(fā)生破碎,從而盡可能地提高實驗照片的時間及空間分辨率。實驗照片像素密度為33~40 pixels/mm,拍攝速率為5×104s-1,單幀曝光時間為0.37~1.00 μs。本實驗液滴破碎的特征時間為0.5~1.0 ms,其中受液霧干擾較小、變形圖像較清晰的初期破碎過程持續(xù)0.15~0.20 ms;因此在實驗總有效時長內(nèi)可獲得25~50幅圖像,其中初期變形圖像7~10幅。

表1 實驗參數(shù)Table 1 Experimental parameters

實驗液滴介質(zhì)為純凈水,與室溫空氣之間的表面張力系數(shù)約為70 mN/m。氣流參數(shù)(氣流密度和速度)的改變通過設置不同的激波管高、低壓段壓力來實現(xiàn)。表1給出了本文實驗工況的具體參數(shù),表中雷諾數(shù)Re=ρgugd0/μg,Mas為運動激波馬赫數(shù),ts,0、td,0和κ分別為分離發(fā)展特征時間、液滴變形特征時間及兩者的比值,定義式由后文給出。測試激波馬赫數(shù)為1.37~2.60,液滴直徑為2.6~3.8 mm,基于激波波后氣體屬性計算所得韋伯數(shù)范圍為2 100~2 700,氣動雷諾數(shù)在104量級;由于液滴介質(zhì)相同且尺度相差不大,所有實驗奧內(nèi)佐格數(shù)基本相當,約為2×10-3。參數(shù)的調(diào)整主要體現(xiàn)在氣流密度和速度的改變,其中密度范圍為0.16~1.31 kg/m3,速度范圍為186~641 m/s。圖1中給出了本文實驗的We-Oh參數(shù)范圍。參照Theofanous等[5]關于液滴破碎模式與We-Oh的依賴關系圖,本文實驗條件均位于KH-SIE破碎機制(以界面KH不穩(wěn)定性和剪切誘導夾帶為特征)區(qū)間內(nèi)。相比于Joseph等[8]、Theofanous等[3]的研究工作,本文實驗的We與Oh均處于一個變化較小的區(qū)間內(nèi)。

1.2 數(shù)值模擬方法

為解釋液滴呈現(xiàn)不同變形模態(tài)的原因,本文中將破碎初期,變形尺度相對不大的液滴,近似看作一個同等直徑的剛性球體;通過軸對稱外流流場的數(shù)值模擬,獲得圓球表面氣動壓力分布和摩擦力(切應力)分布的演化過程,用于對液滴初期變形現(xiàn)象進行評估、分析和預測。

外流模擬采用一套成熟的、基于有限體積方法的二維可壓縮數(shù)值模擬程序—VAS2D[14]。計算域如圖2(a)所示。計算總網(wǎng)格數(shù)約為40 000,在球體壁面周向分布的網(wǎng)格數(shù)為300,邊界層網(wǎng)格最小高度為1 μm。圖2(b)給出了典型實驗條件下邊界層分離點附近的速度場。從圖中可知,速度邊界層厚度為30~50 μm,含15~20個網(wǎng)格,足以刻畫邊界層內(nèi)部及分離發(fā)生的流動細節(jié)。

1.3 液滴變形理論

激波掃過液滴后,流動分離流場在液滴下游建立并發(fā)展穩(wěn)定。外流氣動力對液滴初期變形的驅(qū)動機制可分為2類:一是界面上的剪切摩擦誘導出液滴表面(液滴內(nèi)部邊界層)的周向流動,在表面局部形成液體量的堆積或稀疏;二是表面壓力分布的不均衡,對液滴構成局部擠壓和拉伸效應。2種驅(qū)動機制在液滴表面形成的徑向加速度,可分別用下式:

(1)

表示[13]。在本文所涉及的參數(shù)范圍內(nèi),壓力拉伸所誘導生成的加速度比剪切堆積所誘導生成的加速度高約2個數(shù)量級,因此在下文的分析中,對液滴變形的預估忽略剪切堆積的影響,即任意時刻t1液滴表面θ位置的徑向速度和變形量由:

(2)

給出。式(1)~(2)中τ、p分別為t時刻液滴表面θ位置的瞬態(tài)摩擦和壓力,r0為液滴初始半徑。ητ、ηp為修正參數(shù),用以修正因黏性和表面張力導致的液滴變形量的減小。

2 結果與討論

2.1 液滴的初期變形

對液滴變形和破碎現(xiàn)象的描述一般以We和Oh為主要控制參數(shù)。然而,在維持組合控制參數(shù)(即We和Oh)不變的基礎上變動某些原始參數(shù)(改變氣流密度、速度或液滴直徑),其所得實驗現(xiàn)象的細節(jié)實際也不盡相同,對于這種差異當前仍缺少充分的認識[13]。基于這一考慮,本文中進行了一組液滴破碎實驗,一方面將實驗的We和Oh維持在一個變化較小的區(qū)間內(nèi),以保證液滴的破碎遵循相同的破碎機制;另一方面通過調(diào)整來流的密度和速度,獲得不同的液滴變形圖像。

a0=pd/(ρld0)

(3)

則變形特征時間可定義為:

(4)

式中:pd為激波后氣流動壓。因此,在不同工況下,選取t/td,0相同的時刻進行比較,可以保證液滴處在相同的變形階段。另一方面,Pilch等[9]、Thefanous等[6]給出了SIE機制下液滴破碎時間:

(5)

式(4)和(5)中所定義的液滴變形與液滴破碎特征時間具有相同的形式,因此相同的變形階段也對應相同的破碎階段。在本文中,統(tǒng)一采用t/td,0~0.2時刻進行不同工況下液滴變形的比較。一方面,該時刻的液滴已發(fā)生較明顯的變形,易于對不同模態(tài)加以區(qū)分;另一方面,液滴赤道附近產(chǎn)生的液霧尚未對液滴背風面的發(fā)展圖像產(chǎn)生明顯的遮蔽,整個背風面的形態(tài)清晰可見。

圖3給出了5種實驗條件下,t/td,0~0.2時刻的液滴變形照片。基于背風面變形發(fā)展模態(tài)的區(qū)別,液滴的變形可明顯地分為2類:(a)實驗I中,液滴背風面存在孤立的環(huán)形突起,并向外發(fā)展至較大尺度。在后續(xù)發(fā)展中,該突起頂部會因氣流剪切形成液霧。(b)在氣流密度較低的幾次實驗(II~V)中,液滴背風面存在3個環(huán)形突起,且高度較平均。在變形的中期,受表面張力約束,這些突起均難以進一步向外擴展,因此對液霧的形成幾乎沒有貢獻。

利用1.2節(jié)中的數(shù)值模擬和1.3節(jié)中的理論方法,可以對破碎初期的液滴變形進行模擬和估算。圖4給出了實驗I和II中,液滴實際變形與數(shù)值模擬的對比,計算中修正參數(shù)ηp=0.05。2種來流條件下,液滴背風面上突起的數(shù)量發(fā)展幅度均可以得到刻畫。為探索液滴出現(xiàn)不同形態(tài)的原因,本文下一節(jié)對液滴背風面流場的形成和壓力特征進行總結和分析。

2.2 環(huán)形突起形態(tài)差異的原因分析

從2.1節(jié)可知,在相近的We條件下,液滴背風面環(huán)形突起的發(fā)展會呈現(xiàn)出不同的模態(tài),導致破碎前液滴整體的形狀具有較大差別。這種突起發(fā)展模態(tài)的差異主要受密度的影響:發(fā)展程度較高的孤立突起,多出現(xiàn)于氣流密度較高的實驗中;氣流密度較低時,突起發(fā)展幅度較平均。由式(2)可知,利用等直徑剛性球體壓力數(shù)據(jù)可以預測液滴變形。因此,本節(jié)中通過歸納球體表面分離形成過程中的壓力演化,對突起發(fā)展呈現(xiàn)的多種模態(tài)進行解釋。

2.2.1液滴背風面分離的發(fā)展過程

圖6給出了激波繞射階段中,圓球附近流場的發(fā)展過程。圖中利用x*=x/d0,y*=y/d0進行長度參數(shù)的歸一化。入射激波在球體表面發(fā)生常規(guī)反射,并在接近圓球赤道附近時轉變?yōu)轳R赫反射。隨后彎曲的馬赫桿掃過圓球背風面,在后駐點處相交并向迎風面移動,強度減弱直至完全耗散。這一階段持續(xù)時間較短,通常小于20 μs。本文條件下,激波繞射階段液滴表面無明顯變形發(fā)生。

激波掃過圓球表面后,回流區(qū)內(nèi)的復雜流場逐漸發(fā)展形成并趨于穩(wěn)定,穩(wěn)定后的流場具有3個層次的渦結構。圖7給出了圓球附近分離形成和發(fā)展的過程。馬赫桿相交于背風面駐點后向上游移動,在邊界層內(nèi)誘發(fā)形成強烈的逆壓梯度,并形成一個主渦(渦A)。在圓球背風面的流動中,渦A逐漸占據(jù)主導地位,其中心附近壓力最低。此時圓球表面存在2個明顯的低壓區(qū)域,如圖8虛線所示。隨后渦A內(nèi)部形成二次渦(渦B),并向外發(fā)展,將主渦A分割為A1和A2兩段,見圖7(c)。此后流場結構趨于穩(wěn)定,渦系結構沒有顯著變化。穩(wěn)定階段圓球表面的低壓點有4個,分別位于迎風面分離點附近,以及渦A1、B與A2的內(nèi)部,如圖8實線所示。A1和B內(nèi)部可能形成尺度更小的三次渦,但對圓球表面壓力分布影響不大。上述背風面流場結構,與文獻[15]中的α和sub-α模式吻合。

對應圖8中圓球表面的壓力分布,圖9給出了與壓力分布曲線一一對應的液滴表面法向加速度曲線。從圖中可以看出,壓力較低的區(qū)域與液面法向加速度大體一致。在分離過程中,加速度具有2個相對孤立的峰值;在分離穩(wěn)定后,加速度有4個較明顯的峰值,對應實驗中尺度差異相對較小的4個突起。

2.2.2分離發(fā)展特征時間

從2.2.1節(jié)可知,在分離發(fā)展的不同階段,外部流場在液滴表面施加的壓力具有2種典型的分布形態(tài),并誘導出不同的法向加速度分布。2種壓力分布形態(tài)分別對應分離產(chǎn)生過程(主渦A主導)和分離完全穩(wěn)定后(主渦A1、A2和二次渦B共同作用)。2個階段的作用時間決定了2種壓力分布對于液滴變形影響的主次關系,因此,本節(jié)中對分離產(chǎn)生和穩(wěn)定的速度進行討論。

分離形成過程中,二次渦B的出現(xiàn)是一個可量化的時間節(jié)點。從激波經(jīng)過到二次渦B出現(xiàn)的時間間隔,可以認為正比于分離發(fā)展的特征時間。圖10給出了不同來流條件下,從激波經(jīng)過到二次渦B出現(xiàn)的間隔時間。從圖中可以看出,利用:

ts,0=d0/ug

(6)

進行歸一化,則次渦B出現(xiàn)于t/ts,0≈0.92。

將式(4)與(6)進行對比可知,分離生成的特征時間ts,0與液滴變形的特征時間td,0之比為密度比的平方根,因此定義時間比:

(7)

這一量綱一數(shù)標志著分離過程對于突起生成的貢獻程度。圖3列出的5個實驗工況各自的ts,0、td,0和κ值均在表1中給出。顯然,由于實驗中液體密度恒定,氣流密度較大時κ也較大,圖8中虛線給出的2個低壓區(qū)會在液滴表面持續(xù)相對較長的時間,因此在50°處的突起會發(fā)展至相對更大的幅度。

2.2.3環(huán)形突起發(fā)展的歸一化

量綱一化后分離發(fā)展過程主要受雷諾數(shù)和氣流馬赫數(shù)影響,并具有高度的相似性。利用這一流動特征,液滴外部流場演化可以被拆分為t/ts,0<4的非定常過程和t/ts,0>4的近似定常過程,并分別進行參數(shù)歸一化。液滴的整體變形,可以看作2階段影響的線性疊加。量綱一加速度疊加時的疊加系數(shù),即2階段分別持續(xù)時間的不同,將導致液滴形態(tài)上的區(qū)別。利用這一方法,κ對液滴變形的影響可以得到直接的體現(xiàn)。通過這一階段劃分方法,并將整體穩(wěn)定階段的流動假設為定常,使式(2)給出的液滴在t/tb=0.2內(nèi)發(fā)生的變形可以由:

(8)

分解為2部分:Δrunsteady和Δrsteady分別為非定常過程所誘導產(chǎn)生的變形和穩(wěn)定后壁面壓力誘導的變形。式中a0為式(3)中定義的突起發(fā)展的特征加速度,3項量綱一加速度:

分別對應非定常過程中的變形、非定常過程產(chǎn)生初速度在穩(wěn)定階段形成的變形及穩(wěn)定階段加速度產(chǎn)生的變形。圖11給出了典型條件下(工況I),液滴表面量綱一加速度的分布情況。3項量綱一加速度分布均會受氣流馬赫數(shù)和雷諾數(shù)的影響,此處不進行深入討論;但在本文參數(shù)范圍內(nèi),加速度極值的個數(shù)和相對大小在不同氣流馬赫數(shù)Mag和Re條件下具有相同的規(guī)律。圖12給出了實驗I與II中,在t/td,0=0.2時刻,非定常和定常過程分別產(chǎn)生的變形量。從該圖可以看出,氣流密度的降低(即κ的減小)會削弱非定常組分的影響,使流場穩(wěn)定后的壓力所誘導的變形占據(jù)主導地位。在We=2 370,Re=4.2×104,Mag=0.483的條件下,通過改變液滴直徑和氣流密度對κ進行調(diào)節(jié),液滴的變形由圖13給出。在背風面50°附近,環(huán)形突起的尺度受κ影響較小;而其余位置突起高度與κ呈明顯的正相關趨勢。因此,在κ較低的條件下,突起尺度較平均;κ較高時,50°附近區(qū)間內(nèi)會形成尺度較大的孤立突起。這一趨勢與實驗現(xiàn)象一致。

3 結 論

(1)通過實驗對相近We條件下的液滴變形和破碎過程進行觀測,實驗結果表明,盡管We相近,不同實驗條件下液滴的變形仍呈現(xiàn)出不同的模態(tài),其主要區(qū)別在于表面隆起液環(huán)位置和發(fā)展速度的不同。

(2)從分離的生成到穩(wěn)定,液滴表面的低壓區(qū)域由2個向4個演化。對應的液滴表面加速度峰值也由2個變?yōu)?個。這一過程可以用ts,0進行時間歸一化。

(3)液滴呈現(xiàn)出不同的變形模態(tài),其主要原因是分離形成過程在液滴變形中的參與程度不同。在來流密度較小的條件下,分離過程較快結束,穩(wěn)定的分離主導了液滴的變形;在來流密度較大的條件下,分離過程對液環(huán)形成影響顯著。分離與變形特征時間之比κ標志著分離發(fā)展過程對液環(huán)形成的貢獻程度。

(4)液環(huán)形成可以看作分離發(fā)展階段變形、初速度誘導后期變形與穩(wěn)定階段變形3個組分之和。3個組分均可以進行較好的參數(shù)歸一化,并具有鮮明的單峰值/多峰值特征。以κ作為各部分變形的權重,將歸一化后的參數(shù)進行線性疊加,即可預測不同來流條件下液滴的變形情況。這進一步說明了氣液密度比對液滴變形的影響機理。

參考文獻:

[1] 費立森.煤油在冷態(tài)超聲速氣流中噴射和霧化現(xiàn)象的初步研究[D].合肥:中國科學技術大學,2007:1-12.

[2] 萬云霞,黃勇,朱英.液體圓柱射流破碎過程的實驗[J].航空動力學報,2008,23(2):208-214.

WAN Yunxia, HUANG Yong, ZHU Ying. Experiment on the breakup process of free round liquid jet[J]. Journal of Aerospace Power, 2008,23(2):208-214.

[3] THEOFANOUS T G, LI G J, DINH T N. Aerobreakup in rarefied supersonic gas flows[J]. Journal of Fluids Engineering, 2004,126(4):516-527.

[4] THEOFANOUS T G, LI G J. On the physics of aerobreakup[J]. Physics of Fluids, 2008,20(5):052103.

[5] THEOFANOUS T G, MITKIN V V, NG C L, et al. The physics of aerobreakup: II[J]. Physics of Fluids, 2012,24(2):022104.

[6] THEOFANOUS T G. Aerobreakup of Newtonian and viscoelastic liquids[J]. Annual Review of Fluid Mechanics, 2011,43:661-690.

[7] CHANG C H, DENG X, THEOFANOUS T G. Direct numerical simulation of interfacial instabilities: A consistent, conservative, all-speed, sharp-interface method[J]. Journal of Computational Physics, 2013,242:946-990.

[8] JOSEPH D D, BELANGER J, BEAVERS G S. Breakup of a liquid drop suddenly exposed to a high-speed airstream[J]. International Journal of Multiphase Flow, 1999,25(6):1263-1303.

[9] PILCH M, ERDMAN C A. Use of breakup time data and velocity history data to predict the maximum size of stable fragments for acceleration-induced breakup of a liquid drop[J]. International Journal of Multiphase Flow, 1987,13(6):741-757.

[10] WIERZBA A, TAKAYAMA K. Experimental investigation of the aerodynamic breakup of liquid drops[J]. AIAA Journal, 1988,26(11):1329-1335.

[11] 金仁瀚,劉勇,朱冬清,等.初始直徑對單液滴破碎特性影響的試驗[J].航空動力學報,2015,30(10):2401-2409.

JIN Renhan, LIU Yong, ZHU Dongqing, et al. Experiment on impact of initial diameter on breakup characteristic of single droplet[J]. Journal of Aerospace Power, 2015,30(10):2401-2409.

[12] 王超,吳宇,施紅輝,等.液滴在激波沖擊下的破裂過程[J].爆炸與沖擊,2016,36(1):129-134.

WANG Chao, WU Yu, SHI Honghui, et al. Breakup process of a droplet under the impact of a shock wave[J]. Explosion and Shock Waves, 2016,36(1):129-134.

[13] 易翔宇,朱雨建,楊基明.激波誘導高速氣流中液滴的初期變形[J].爆炸與沖擊,2017,37(5):853-862.

YI Xiangyu, ZHU Yujian, YANG Jiming. Early-stage deformation of liquid drop in shock induced high-speed flow[J]. Explosion and Shock Waves, 2017,37(5):853-862.

[14] SUN M, SAITO T, TAKAYAMA K, et al. Unsteady drag on a sphere by shock wave loading[J]. Shock Waves, 2005,14(1/2):3-9.

[15] KALITA J C, SEN S. Unsteady separation leading to secondary and tertiary vortex dynamics: The sub-α- and sub-β-phenomena[J]. Journal of Fluid Mechanics, 2013,730:19-51.

猜你喜歡
變形實驗
記一次有趣的實驗
微型實驗里看“燃燒”
談詩的變形
中華詩詞(2020年1期)2020-09-21 09:24:52
做個怪怪長實驗
“我”的變形計
變形巧算
例談拼圖與整式變形
會變形的餅
NO與NO2相互轉化實驗的改進
實踐十號上的19項實驗
太空探索(2016年5期)2016-07-12 15:17:55
主站蜘蛛池模板: 国产日韩欧美视频| 欧美日韩免费| 日本色综合网| 好吊色妇女免费视频免费| 国产在线精彩视频论坛| 伊人色天堂| 色综合婷婷| 乱系列中文字幕在线视频| 91久久精品日日躁夜夜躁欧美| 日本一区二区不卡视频| 最新无码专区超级碰碰碰| 国产在线小视频| 亚洲日韩在线满18点击进入| 九九热视频精品在线| 亚洲人成网站日本片| 色婷婷在线播放| 欧美成人影院亚洲综合图| 欧美a√在线| 香蕉国产精品视频| 青青草91视频| 亚洲精品无码av中文字幕| 久久久久人妻精品一区三寸蜜桃| 日韩免费毛片视频| 亚洲高清在线播放| 成年午夜精品久久精品| a色毛片免费视频| 久久77777| 成人在线观看一区| 在线日本国产成人免费的| 亚洲综合狠狠| 伊人久久婷婷五月综合97色| 精品国产www| 色婷婷成人网| 久久久久国产一区二区| 直接黄91麻豆网站| 欧美区一区| 欧美一级大片在线观看| 啪啪永久免费av| 国产在线视频欧美亚综合| 亚洲欧美日韩成人在线| 亚洲综合色吧| 亚洲国产精品美女| 一级毛片在线免费看| 色精品视频| 亚洲国产无码有码| 国产99久久亚洲综合精品西瓜tv| 99精品这里只有精品高清视频| 在线国产三级| 亚洲国产精品VA在线看黑人| 国产精品嫩草影院视频| 真人高潮娇喘嗯啊在线观看| 欧美日韩精品在线播放| 黑人巨大精品欧美一区二区区| 成人小视频网| 国产亚洲高清视频| 国产成人1024精品| 欧美激情视频一区二区三区免费| 国产成人啪视频一区二区三区 | 日本精品一在线观看视频| 狠狠干欧美| 91在线视频福利| 农村乱人伦一区二区| 欧美精品亚洲日韩a| 久久久波多野结衣av一区二区| 欧日韩在线不卡视频| 日韩二区三区无| 漂亮人妻被中出中文字幕久久| 三级国产在线观看| 中国美女**毛片录像在线| 日a本亚洲中文在线观看| 久久人午夜亚洲精品无码区| 午夜精品久久久久久久2023| 99精品热视频这里只有精品7| 国产激爽爽爽大片在线观看| 91啦中文字幕| 一级毛片a女人刺激视频免费| 日韩精品免费一线在线观看| 91久久精品日日躁夜夜躁欧美| 国产麻豆精品手机在线观看| 白丝美女办公室高潮喷水视频| 99久久精品国产麻豆婷婷| 操国产美女|