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

汽車風(fēng)窗玻璃外流場特性及其對雨滴運動影響的研究?

2018-11-15 01:47:40姜立標(biāo)丘華川
汽車工程 2018年10期
關(guān)鍵詞:模型

姜立標(biāo),丘華川

(1.吉林大學(xué),汽車仿真與控制國家重點實驗室,長春 130025;2.華南理工大學(xué)機械與汽車工程學(xué)院,廣州 510640)

前言

傳統(tǒng)橡膠雨刮器具有成本低廉、結(jié)構(gòu)較為簡單、功耗低的優(yōu)點[1-2],但性能一般。為滿足人們對汽車駕駛性能越來越高的要求,英國邁凱倫公司提出了使用超聲波系統(tǒng)取代傳統(tǒng)雨刷結(jié)構(gòu)的方案,不僅可降低雨刷骨架風(fēng)阻,還可取消雨刷電機,減輕車質(zhì)量。為達(dá)到流動減阻有利于雨滴的運動,通常采用超疏水材料對風(fēng)窗玻璃表面進(jìn)行處理。文獻(xiàn)[3]中借助超聲波微流體驅(qū)動技術(shù),液滴產(chǎn)生受控運動,利用壓電陶瓷的逆壓電效應(yīng)在疏水玻璃表面產(chǎn)生超聲行波實現(xiàn)對液滴的驅(qū)動,液滴移動速度約為10mm/s。要在汽車風(fēng)窗玻璃上通過超聲波使雨滴產(chǎn)生受控運動,應(yīng)同時考慮超聲波的形式與傳播特點和玻璃表面疏水特性與氣流流速等因素的影響。在壓強變化、氣體流速和重力等多種因素作用下,液滴的液/氣界面會發(fā)生失穩(wěn)現(xiàn)象,并以不同的形態(tài)進(jìn)行演化,因此必須對汽車風(fēng)窗玻璃外流場特性進(jìn)行研究,包括壓力分布特性和附著氣流速度場特性。

隨著計算機技術(shù)和湍流理論的不斷發(fā)展,計算流體力學(xué)數(shù)值模擬方法廣泛應(yīng)用于汽車空氣動力學(xué)的研究。本文中利用CFD數(shù)值模擬方法,在三維建模軟件CATIA中建立風(fēng)窗玻璃模型,然后導(dǎo)入多物理場耦合分析軟件COMSOL中進(jìn)行數(shù)值研究,分析了風(fēng)窗玻璃外表面的壓力分布特性和附著氣流速度場特性,通過控制變量法研究了氣流流入速度和風(fēng)窗玻璃傾角對流場特性的影響。最后研究了風(fēng)窗玻璃外流場特性對雨滴運動的影響機理,指出不同質(zhì)量的雨滴在不同氣流速度下的運動狀態(tài)。

1 基本方程

文獻(xiàn)[4]中采用標(biāo)準(zhǔn)k-ε模型和Realizable k-ε模型對車身外流場進(jìn)行分析后指出,兩種模型的對稱面上壓力分布基本吻合。風(fēng)窗玻璃外流場特性可進(jìn)一步用于雨滴的運動機理分析,這須添加更多的物理場來進(jìn)行模擬,其次標(biāo)準(zhǔn)k-ε模型具有計算量小、阻力收斂速度快、收斂后阻力波動值小等優(yōu)點,而且計算結(jié)果與實驗值能很好地吻合,誤差范圍在3%以內(nèi),采用多物理場耦合分析軟件COMSOL提供的標(biāo)準(zhǔn)k-ε模型[5]。

1.1 標(biāo)準(zhǔn)k-ε模型

汽車風(fēng)窗玻璃表面的外流場一般為定常、等溫和不可壓縮的三維流動,基于動量守恒和質(zhì)量守恒,可得空氣流動的動量方程和連續(xù)性方程[6-7]為

式中:ρ為流體密度;xi和xj分別為流場中沿i和j方向上的坐標(biāo)分量;ui和uj分別為流場中沿i和j方向上的平均相對速度分量;μe為流體有效動力黏度;p為壓強;Si為廣義源項。

為得到方程的解,須用湍流模型來封閉。k-ε雙方程模型中引入兩個附加的傳輸方程:湍流動能k方程和耗散率ε方程。

標(biāo)準(zhǔn)k-ε湍流模型由文獻(xiàn)[8]中提出。有效動力黏度μe等于分子動力黏度μ加上渦團(tuán)黏度μT,即

式中:k為湍流動能;ε為湍流耗散率。

根據(jù)文獻(xiàn)[6],設(shè)定 Cμ= 0.09,C1= 1.44,C2=1.92,σε= 1.3,σk= 1.0。

1.2 疏水風(fēng)窗玻璃表面液滴模型

在疏水風(fēng)窗玻璃表面上的雨滴,在重力、慣性力、表面張力、接觸角滯后阻力、氣流作用力和黏性耗散作用力的共同作用下,可能處于靜止、向下滑動和向上滑動的狀態(tài),主要取決于重力、接觸角滯后阻力和氣流作用力的影響[9]。傾斜疏水表面雨滴模型如圖1所示。

圖1 傾斜疏水表面雨滴模型

2 數(shù)值計算模型

本文中研究對象為某國產(chǎn)自主品牌的風(fēng)窗玻璃,根據(jù)實際風(fēng)窗玻璃1∶1建立幾何曲面模型,然后將所建立的模型導(dǎo)入多物理場耦合分析軟件COMSOL中進(jìn)行數(shù)值分析。

2.1 計算域的確定

綜合考慮計算域尺寸和計算域與風(fēng)窗玻璃均為對稱結(jié)構(gòu),為了節(jié)省計算資源,采用對稱建模方法,取其縱對稱面左側(cè)的一半進(jìn)行計算。構(gòu)建尺寸為1200mm×2000mm×2000mm的“二分之一”計算域模型,如圖2所示。定義風(fēng)窗玻璃傾角α為風(fēng)窗玻璃對稱面外表面的上、下兩端點連線與水平面的夾角。

圖2 計算域的三維結(jié)構(gòu)

2.2 網(wǎng)格控制域

風(fēng)窗玻璃氣流模型研究的重點區(qū)域只在于風(fēng)窗玻璃外表面附近的流場,提取風(fēng)窗玻璃的外表面,沿水平方向平移使其向外產(chǎn)生200 mm厚度的區(qū)域,并定義該區(qū)域為網(wǎng)格控制域,如圖3所示。

圖3 網(wǎng)格控制域

2.3 邊界條件

設(shè)置邊界條件如表1所示。

表1 邊界條件設(shè)置

2.4 網(wǎng)格策略

構(gòu)建的網(wǎng)格為滿足實際工程的要求,采用用戶控制網(wǎng)格。基于網(wǎng)格控制域,采用Delaunay方法自由剖分四面體網(wǎng)格,映射生成金字塔和三棱柱等混合網(wǎng)格。風(fēng)窗玻璃表面附近進(jìn)行邊界層加密處理,同時考慮到模型中成角的部位對求解的收斂有影響,且往往容易產(chǎn)生應(yīng)力集中或出現(xiàn)奇異,使用角細(xì)化網(wǎng)格策略對風(fēng)窗玻璃的4個角(包括4條邊)附近的網(wǎng)格進(jìn)行網(wǎng)格加密,模型的網(wǎng)格劃分如圖4所示,計算域的網(wǎng)格模型單元總數(shù)為552 684。

3 仿真結(jié)果與試驗驗證

本研究中求解線性方程組的迭代方法采用廣義極小殘差算法[10-11],具有穩(wěn)定性好和收斂速度快等優(yōu)點,可以利用有限的內(nèi)存空間求解復(fù)雜的流場代數(shù)方程組。

3.1 流場特性分析

為研究風(fēng)窗玻璃表面流場特性,設(shè)定傾角α為45°,氣流速度u0為7.31 m/s,進(jìn)行仿真,結(jié)果如圖5~圖8所示。其中圖5為對稱面及其周圍的速度分布。圖6為風(fēng)窗玻璃外表面壓力分布。圖中箭頭方向表征附著氣流的流向,顏色深淺表征附著氣流的流速。圖7為風(fēng)窗玻璃外表面附著氣流速度場。圖8為風(fēng)窗玻璃對稱面上的壓力曲線。

圖4 模型的網(wǎng)格劃分

圖5 對稱面及其周圍的速度分布

圖6 風(fēng)窗玻璃表面壓力分布

圖7 風(fēng)窗玻璃表面附著氣流速度分布

圖8 風(fēng)窗玻璃對稱面上的壓力曲線

由圖可見,當(dāng)氣流到達(dá)風(fēng)窗玻璃底端,由于風(fēng)窗玻璃的存在,氣流的速度降低。同時由于底端氣流直接流過風(fēng)窗玻璃,形成了一個滯區(qū),該區(qū)域具有正壓力。壓力瞬間增大,氣流沿風(fēng)窗玻璃表面向后流動,氣流速度逐漸增大,因此壓力逐漸減小。到風(fēng)窗玻璃頂端與左側(cè)邊緣附近,因氣流流速較大而出現(xiàn)負(fù)壓力,可以看出,整個風(fēng)窗玻璃的絕大部分區(qū)域承受正壓力。

3.1.1 壓力分布特性

以零壓力線為分界,分為正壓力區(qū)和負(fù)壓力區(qū)(法向向外為正)。正壓力區(qū)的面積明顯大于負(fù)壓力區(qū),負(fù)壓力區(qū)分布在靠近風(fēng)窗玻璃外表面上邊緣的小部分區(qū)域和左右邊緣的極小區(qū)域。存在一個壓力絕對值最大點,位于風(fēng)窗玻璃對稱面與其外表面的交線上,并靠近下邊緣處,該點的壓力為正壓力,且|p|max= 53.6260Pa。

對于正壓力區(qū),其等壓線比較稀疏,表明正壓力區(qū)的壓力變化速率比較小;對于負(fù)壓力區(qū),其等壓線比較密集,表明負(fù)壓力區(qū)的壓力變化速率比較大。

3.1.2 附著氣流速度場特性

存在一個流速最小值ufmin=6.67218×10-2m/s,位于風(fēng)窗玻璃對稱面與其外表面的交線上,并靠近下邊緣處,附著氣流由這點開始沿著玻璃外表面向四周流動,流速逐漸變大。在邊界處流速最大,ufmax=7.93756m/s。

考慮到邊界上的點受流域空間突變的影響,而且每個點的受影響程度不同,因此剔除邊界上的點。在正壓力區(qū),風(fēng)窗玻璃外表面上附著氣流的流速與風(fēng)窗玻璃外表面的壓力呈負(fù)相關(guān),流速越大的地方正壓力越小;在負(fù)壓力區(qū),風(fēng)窗玻璃外表面上附著氣流的流速與風(fēng)窗玻璃外表面的壓力呈正相關(guān),流速越大的地方負(fù)壓力越大。即玻璃外表面上某點的壓力代數(shù)值越大,流速越小。

3.2 實驗裝置

實驗裝置由鋁合金固定支架、軸流式風(fēng)機和風(fēng)窗玻璃等組成。表2為軸流式風(fēng)機的具體參數(shù)。

表2 軸流式風(fēng)機的具體參數(shù)

軸流式風(fēng)機輸出風(fēng)速為

式中:Q為軸流式風(fēng)機的輸出風(fēng)量;d為軸流式風(fēng)機的風(fēng)道直徑。

將表2數(shù)據(jù)代入式(10)得到風(fēng)機的輸出風(fēng)速u=7.31 m/s。

3.3 實驗結(jié)果

風(fēng)窗玻璃外表面附著氣流的流向特性實驗的原理為:使用軸流式風(fēng)機在風(fēng)窗玻璃前方形成平穩(wěn)氣流,以此來模擬風(fēng)窗玻璃的外流場。風(fēng)窗玻璃前方是否為平穩(wěn)氣流,可以借助實驗驗證。具體方法為:在軸流式風(fēng)機出風(fēng)口的安全網(wǎng)上系上長絲帶,觀察其在流場中的形狀表現(xiàn)為水平來近似保證。進(jìn)一步通過搭建實物進(jìn)行實驗,觀察粘連在玻璃外表面上細(xì)絲帶的飄動方向,獲得風(fēng)窗玻璃外表面附著氣流的流向特性。經(jīng)過多次實驗,絲帶的流向規(guī)律基本如圖9所示。

實驗結(jié)果與圖9的仿真結(jié)果對比,可以看到實驗中除了個別的細(xì)絲帶外,大多數(shù)細(xì)絲帶對應(yīng)點的附著氣流方向與 α=45°,u0=7.31m/s下的仿真結(jié)果基本一致。從實驗結(jié)果可見,左下角區(qū)域的細(xì)絲帶基本不動,說明此處的附著氣流流速很小,為附著氣流的速度值較小的區(qū)域,而在此區(qū)域外,可以很清楚地觀察到細(xì)絲帶向四周飄散。細(xì)絲帶體現(xiàn)出來的附著氣流的流向規(guī)律為:風(fēng)窗玻璃左下角存在附著氣流速度的發(fā)散點,氣流流向從該點附近區(qū)域開始發(fā)散,并沿著玻璃表面向四周流動。

圖9 風(fēng)窗玻璃表面絲帶流向規(guī)律

為進(jìn)一步研究風(fēng)窗玻璃表面氣流流速的大小,測試圖9中所示11個測點的氣流速度,測速裝置采用國產(chǎn)希瑪分體式風(fēng)速計AS8336,其分辨率為0.001m/s、準(zhǔn)確度為±3%rdg±0.1。 將風(fēng)速計設(shè)置為最大值擋位,分別多次測量風(fēng)窗玻璃表面11個測點的氣流速度后取其平均值以減少測量誤差,并與對應(yīng)點的仿真模型的氣流速度進(jìn)行對比,如表3所示。

綜合人為測量手段、風(fēng)速計誤差和軸流式風(fēng)機產(chǎn)生風(fēng)速大小等因素,由表3可見,實驗測量的風(fēng)速與仿真模型的風(fēng)速在誤差允許范圍內(nèi)非常接近。因此可以得出,所建立的風(fēng)窗玻璃湍流模型與實際情況接近,所得結(jié)論能夠反映風(fēng)窗玻璃外流場特性。

4 各參數(shù)對流場特性的影響

圖10 不同u0下的風(fēng)窗玻璃外表面壓力分布

圖11 不同u0下的風(fēng)窗玻璃表面附著氣流速度分布

4.1 氣流速度對流場特性的影響

風(fēng)窗玻璃傾角α為45°,在氣流速度u0分別為7.31,11.1和16.7m/s的條件下對氣流模型進(jìn)行數(shù)值仿真。不同氣流速度u0下的風(fēng)窗玻璃外表面壓力分布如圖10所示。不同氣流速度u0下的風(fēng)窗玻璃外表面附著氣流速度分布如圖11所示。

4.1.1 氣流速度對壓力分布特性的影響

由圖10可見,改變u0將影響風(fēng)窗玻璃外表面各點壓力的絕對值,u0越大,風(fēng)窗玻璃外表面各點壓力的絕對值越大。但從整體的角度觀察,無論是增大u0還是減小u0,圖10中零壓力線區(qū)域、等壓線的形狀、壓力絕對值最大點的位置基本上沒有發(fā)生變化,但是壓力梯度增加,說明改變u0不影響風(fēng)窗玻璃外表面壓力的整體分布規(guī)律。

4.1.2 氣流速度對附著氣流速度場特性的影響

改變u0將影響風(fēng)窗玻璃外表面各點附著氣流的流速值,u0越大,各點附著氣流的流速越大,且流速的最大值和最小值均變大,但從整體的角度觀察,附著氣流流向的變化規(guī)律并沒有發(fā)生變化,說明改變u0不影響風(fēng)窗玻璃外表面附著氣流流向的變化規(guī)律。

4.2 風(fēng)窗玻璃傾角對流場特性的影響

氣流流速u0=16.7 m/s,在風(fēng)窗玻璃傾角α分別為35°,45°和55°的條件下對氣流模型進(jìn)行數(shù)值模擬。不同α下的風(fēng)窗玻璃外表面壓力分布特性如圖12所示。不同α下的風(fēng)窗玻璃外表面附著氣流速度場特性如圖13所示。

4.2.1 風(fēng)窗玻璃傾角對壓力分布特性的影響

減小α,可使風(fēng)窗玻璃外表面的零壓力線區(qū)域增加,表明零壓力線將下移。對于正壓力區(qū),隨著α的增加,等壓線變得越來越密集,說明壓力變化速率越來越大。

圖12 不同傾角下的風(fēng)窗玻璃外表面壓力分布

圖14 液滴在傾斜表面上隨氣流速度變化的加速度

4.2.2 風(fēng)窗玻璃傾角對附著氣流速度場特性的影響

隨著風(fēng)窗玻璃傾角α的減小,流速的最大值將變小,附著氣流的速度值最小的位置點也隨之下移,同時流速為0~6m/s的區(qū)域和流速為10~18m/s的區(qū)域顯著縮小,而流速為6~10m/s的區(qū)域顯著擴大,說明α減小使附著氣流的整體流速趨于均勻。從圖13中箭頭方向的變化可發(fā)現(xiàn),減小α,對于附著氣流的流向,將有更多的附著氣流向上流動。同時,驗證了壓力梯度較大時,其對應(yīng)的速度梯度也較大。

5 流場特性對雨滴運動的影響分析

文獻(xiàn)[12]中研究了在傾斜角度為35°的疏水表面上不同質(zhì)量的液滴隨氣流速度變化的加速度,圖14為液滴在傾斜表面上隨氣流速度變化的加速度。由圖可見:在氣流速度為7~11m/s時,液滴處于靜止?fàn)顟B(tài);在氣流速度大于11m/s時,液滴向上滑動。

根據(jù)Laplace效應(yīng),彎曲液面會在液體內(nèi)部產(chǎn)生附加壓強,其大小為

式中:σ為表面張力系數(shù);R1和R2分別為液滴界面曲率的兩個主軸半徑。

在疏水的風(fēng)窗玻璃表面上,若液滴靜止時,液體內(nèi)部各個方向的附加壓強相等且相互抵消,液滴保持球形。在重力、表面張力、接觸角滯后阻力和氣流作用力的共同作用下,液滴可能靜止、向下滑動或者向上滑動。當(dāng)其發(fā)生變形時,由于各個位置的曲率不同,附加壓強也隨之變化,因此液滴不能保持平衡,發(fā)生向下滑動或者向上滑動。

進(jìn)行液氣兩相流實驗分析,所使用的玻璃尺寸為200mm×150mm×2mm,玻璃傾角為45°。使用疏水涂料對玻璃表面進(jìn)行處理,采用接觸角測量儀(Dataphysics,JC200C1,中國)測定疏水玻璃表面的接觸角。在室溫條件下,取水滴體積10μL,選擇不同的位點測量5次,計算其平均值,測得液滴接觸角為96.5°,接觸角測量圖如圖15所示。

圖15 接觸角測量圖

圖16 所示為氣流速度為3.0m/s條件下液滴的宏觀運動狀態(tài),可以看到大液滴向下移動,融合小液滴后繼續(xù)向下移動,其他小液滴處于靜止?fàn)顟B(tài)。用工業(yè)高速攝像儀(Photron,德國)記錄液滴的微觀運動狀態(tài)。圖17所示為氣流速度為3.0m/s條件下液滴的微觀運動狀態(tài),不同時刻下,液氣自由界面在氣流作用力、表面張力、接觸角滯后阻力和重力等共同作用下,流氣自由界面不斷發(fā)生變形,雖然液滴接觸線為靜止?fàn)顟B(tài),但由于液氣界面處于失穩(wěn)狀態(tài),影響接觸線的移動。圖18為氣流速度為15m/s時雨滴的運動狀態(tài)[13]。可以看到大雨滴直接向上移動,然而小雨滴保持靜止?fàn)顟B(tài)。基于超聲波微流體驅(qū)動技術(shù),將風(fēng)窗玻璃表面的壓力映射為雨滴的氣流驅(qū)動力,對雨滴運動控制具有重要的參考價值和指導(dǎo)意義。

圖16 氣流速度為3.0m/s條件下液滴的宏觀運動狀態(tài)

圖17 氣流速度為3.0m/s條件下液滴的微觀運動狀態(tài)

圖18 氣流速度為15m/s時雨滴的運動狀態(tài)

6 結(jié)論

本文中基于對稱建模方法,通過計算域的構(gòu)建與網(wǎng)格控制域的生成、采用控制重點區(qū)域劃分混合網(wǎng)格策略,建立了汽車風(fēng)窗玻璃氣流模型,通過對風(fēng)窗玻璃外流場特性的研究,得出以下結(jié)論。

(1)在風(fēng)窗玻璃傾角α=45°,氣流流入速度u0=7.31m/s的條件下,附著氣流存在一個流速最小點,位于風(fēng)窗玻璃對稱面與其外表面的交線上,并靠近下邊緣處,附著氣流由這點開始沿著玻璃外表面向四周流動,流速逐漸變大,在邊界處流速最大。除邊界點外,玻璃外表面上某點的壓力代數(shù)值越大,流速越小。在誤差允許的范圍內(nèi),通過實驗驗證了數(shù)值模擬所得出的風(fēng)窗玻璃外表面附著氣流的流向特性和速度大小與實驗數(shù)據(jù)基本吻合。

(2)改變u0只對表面壓力和附著氣流流速的數(shù)值產(chǎn)生影響,不影響玻璃表面壓力分布規(guī)律和附著氣流的流向規(guī)律。隨著α的減小,流速的最大值變小,附著氣流的發(fā)散點下移,附著氣流的整體流速趨于平穩(wěn)且更多的附著氣流向上流動。同時,驗證了壓力梯度較大時,其速度梯度也較大。

(3)通過液氣兩相流實驗分析得到,當(dāng)氣流速度大于雨滴靜止的臨界速度時,速度越大,對雨滴的影響更為顯著。液氣自由界面在氣流作用力、表面張力、接觸角滯后阻力和重力等共同作用下,液氣自由界面不斷發(fā)生變形,液滴接觸線為靜止?fàn)顟B(tài),液氣界面處于失穩(wěn)狀態(tài),影響接觸線的移動。研究風(fēng)窗玻璃流場特性對雨滴運動狀態(tài)的影響,基于超聲波微流體驅(qū)動技術(shù),將風(fēng)窗玻璃表面的壓力映射為雨滴的氣流驅(qū)動力,對雨滴運動控制具有重要的參考價值和指導(dǎo)意義。

猜你喜歡
模型
一半模型
一種去中心化的域名服務(wù)本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數(shù)模型及應(yīng)用
p150Glued在帕金森病模型中的表達(dá)及分布
函數(shù)模型及應(yīng)用
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 在线精品亚洲一区二区古装| 成人av专区精品无码国产| 91久久精品日日躁夜夜躁欧美| 亚洲an第二区国产精品| 久久黄色影院| 亚欧美国产综合| 日本黄色不卡视频| 婷婷久久综合九色综合88| 中文字幕日韩久久综合影院| 午夜一区二区三区| 国产原创演绎剧情有字幕的| 久久网欧美| 99久久99这里只有免费的精品| 好吊妞欧美视频免费| 999国内精品视频免费| 国产欧美日韩在线一区| 免费一看一级毛片| 综合久久五月天| 国产免费自拍视频| 国产精品yjizz视频网一二区| 夜精品a一区二区三区| 欧美色99| 永久天堂网Av| AV不卡在线永久免费观看| 性网站在线观看| 99热国产在线精品99| 欧美19综合中文字幕| 无码国内精品人妻少妇蜜桃视频| 午夜精品区| 国产精品自拍合集| 国产成人精品在线| 欧洲日本亚洲中文字幕| 日韩人妻无码制服丝袜视频| 欧美亚洲激情| 国产精品无码一区二区桃花视频| 亚洲国产av无码综合原创国产| 日韩一区精品视频一区二区| 国产香蕉97碰碰视频VA碰碰看| 欧美黄网在线| 久久semm亚洲国产| 国产在线视频福利资源站| 九九香蕉视频| 午夜少妇精品视频小电影| 精品福利视频导航| 国产精品人人做人人爽人人添| 1024国产在线| 多人乱p欧美在线观看| 国产91视频观看| 丁香六月激情婷婷| 伊人蕉久影院| 久久国产精品夜色| 亚洲精品国产精品乱码不卞| 国产一二三区视频| 伊人色综合久久天天| 国产在线精彩视频二区| 九九九久久国产精品| 亚洲综合经典在线一区二区| 狠狠操夜夜爽| 伊在人亚洲香蕉精品播放| 亚洲av综合网| a毛片在线免费观看| 欧美成人免费午夜全| 综合色区亚洲熟妇在线| 日韩在线第三页| 亚洲乱亚洲乱妇24p| 色悠久久久久久久综合网伊人| 久久国产精品电影| 亚洲日韩久久综合中文字幕| 国产欧美中文字幕| 91口爆吞精国产对白第三集| 国产午夜精品一区二区三| 少妇精品久久久一区二区三区| 露脸一二三区国语对白| 日本欧美成人免费| 久久免费视频6| 亚洲—日韩aV在线| 国产一区亚洲一区| 婷五月综合| 亚洲三级a| 午夜三级在线| AV在线天堂进入| 欧美区一区|