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

液滴撞擊移動及旋轉(zhuǎn)表面過程研究綜述

2023-10-24 10:09:10周易金哲巖楊志剛
關(guān)鍵詞:實驗模型研究

周易,金哲巖,,楊志剛

(1.同濟(jì)大學(xué) 航空航天與力學(xué)學(xué)院,上海 200092;2.上海市地面交通工具空氣動力與熱環(huán)境模擬重點實驗室,上海 201804)

液滴撞擊表面現(xiàn)象廣泛存在于人類日常生產(chǎn)生活中.人們對液滴撞擊有不同的要求,例如噴涂打印、噴霧冷卻、農(nóng)藥噴灑等場景中需要液滴盡可能均勻附著在目標(biāo)位置,而在風(fēng)力發(fā)電機(jī)葉片結(jié)冰、飛機(jī)防除冰等場景中需要盡量減少液滴附著在表面.液滴撞擊現(xiàn)象的深入理解有助于優(yōu)化相關(guān)工程技術(shù),解決對應(yīng)的工程問題.隨著高熱流密度設(shè)備的廣泛運用,傳統(tǒng)散熱方式無法滿足散熱需求,液滴噴霧冷卻技術(shù)則可以更好地解決問題[1];噴涂打印技術(shù)可以實現(xiàn)數(shù)字化的零件生產(chǎn)[2];農(nóng)藥噴灑技術(shù)要求液滴需要均勻地鋪展[3].許多危害現(xiàn)象也與液滴撞擊表面的過程有關(guān),如飛行過程中云中的水滴對飛機(jī)表面沖擊附著帶來的表面結(jié)冰等現(xiàn)象會對飛機(jī)飛行造成影響[4-5].當(dāng)水滴撞擊并附著在旋轉(zhuǎn)帽罩時,由于離心力以及振動對冰的影響,結(jié)出的冰可能會脫落并打壞發(fā)動機(jī)[6].旋轉(zhuǎn)葉片受水滴撞擊結(jié)冰后也會破壞葉片的氣動特性及載荷分布,輕則降低葉片的工作效率,重則引發(fā)安全事故[7-9].

水滴沖擊飛行中的飛機(jī)表面往往涉及到水滴與移動表面間的撞擊.飛機(jī)旋轉(zhuǎn)帽罩、風(fēng)力機(jī)旋轉(zhuǎn)葉片等構(gòu)件的沖擊現(xiàn)象以及覆膜工業(yè)常用的旋涂技術(shù)都與水滴、旋轉(zhuǎn)表面的動力學(xué)規(guī)律有關(guān).其中包含了多個關(guān)鍵的科學(xué)問題,例如液滴撞擊移動及旋轉(zhuǎn)表面后在表面的鋪展行為是如何的,發(fā)生飛濺的條件和飛濺行為是如何的,發(fā)生反彈的條件和反彈行為是如何的等等.要解決上述的科學(xué)問題,需要從液滴撞擊靜止表面的研究體系與結(jié)論出發(fā),通過實驗、理論分析和數(shù)值模擬相結(jié)合的方法,建立起移動及旋轉(zhuǎn)表面的液滴沖擊理論.對于解決相關(guān)工程問題和規(guī)避危害現(xiàn)象,開展液滴撞擊表面的相關(guān)研究尤其對移動及旋轉(zhuǎn)表面的針對性研究具有重要的意義.

在1876 年Worthington[10-11]實驗研究液滴撞擊壁面及液膜的現(xiàn)象.早期的研究往往選擇黏度較大的液體以增加現(xiàn)象持續(xù)的時間[12].隨著高速攝影技術(shù)的發(fā)展,研究人員對液滴撞擊現(xiàn)象進(jìn)行更深入的研究.在1955 年Engel[13]使用化學(xué)方法和高速攝像機(jī)追蹤撞擊面內(nèi)液體的流動研究液滴的鋪展.在1985 年施明恒[14]建立適用于潤濕和非潤濕接觸2 種情況的鋪展模型.Sikalo 等[15-18]于2005 年進(jìn)行一系列研究發(fā)現(xiàn)液滴撞擊表面的現(xiàn)象受到液滴尺寸、撞擊角度、速度、黏度、表面張力、壁面粗糙度等變量的影響.在液滴撞擊表面的研究中,關(guān)于靜止表面(包括固體表面和液膜)的研究已經(jīng)十分深入.針對移動表面的研究仍在逐步開展,雖然目前已經(jīng)有了一些優(yōu)秀的研究成果,但是仍未系統(tǒng)、完善地建立相關(guān)的機(jī)理分析和模型體系.

對液滴撞擊移動表面的相關(guān)文獻(xiàn)進(jìn)行收集和整理,沿著固體和液體表面2 條路線,簡述液滴撞擊移動表面的現(xiàn)象、實驗方法及系統(tǒng);對目前研究的結(jié)果和相關(guān)模型的原理進(jìn)行描述;提出一些液滴撞擊移動表面相關(guān)的有待解決的問題.本研究為科研人員開展后續(xù)的液滴撞擊移動表面的相關(guān)研究提供參考.

1 液滴撞擊移動表面的實驗系統(tǒng)與方法

在液滴撞擊靜止表面的研究中,研究人員常用無量綱數(shù)描述撞擊條件或建立模型,因此在移動表面的研究中引入與表面運動相關(guān)的無量綱數(shù).本研究以下標(biāo)n為表面法向(液滴下落速度方向)量,下標(biāo)s為表面運動(液滴下落速度切向方向)量,下標(biāo)f 為流動液膜.研究中常見的無量綱數(shù)有液滴撞擊韋伯?dāng)?shù)Wen、表面移動韋伯?dāng)?shù)Wes、液滴撞擊雷諾數(shù)Ren、液膜流動雷諾數(shù)Ref、液滴撞擊毛細(xì)數(shù) Can、奧內(nèi)佐格數(shù) Oh ,以及特征量無量綱液膜厚度H和無量綱液膜局部流速U,具體計算式為

式中:D0為液滴撞擊前特征尺寸, ρ 為液滴密度,σ為液體表面張力, μ 為黏性系數(shù),w為渠道寬度,q為液膜流量,Vn為液滴撞擊速度,Vs為表面移動速度,hf為液膜厚度,Uf為液膜局部流速.

1.1 液滴撞擊移動固體表面的實驗系統(tǒng)與方法

研究人員對液滴撞擊移動固體表面設(shè)計了多種實驗,常見的實驗可以分為4 個系統(tǒng):液滴發(fā)生系統(tǒng)、移動表面系統(tǒng)、圖像捕捉系統(tǒng)和后處理系統(tǒng).液滴發(fā)生系統(tǒng)用于產(chǎn)生液滴,控制液滴各項參數(shù),通常使用毛細(xì)管針頭產(chǎn)生液滴;移動表面系統(tǒng)用于設(shè)置移動的固體表面,控制表面各項參數(shù);圖像捕捉系統(tǒng)一般為高速相機(jī),用于捕捉實驗現(xiàn)象;后處理系統(tǒng)負(fù)責(zé)數(shù)據(jù)處理.其中運動表面系統(tǒng)是最重要的系統(tǒng)之一,也體現(xiàn)了幾種常見的設(shè)計思路,常見的運動表面平臺結(jié)構(gòu)可以分為3 大類.如圖1 所示,分別為轉(zhuǎn)輪式運動表面、平移式運動表面和旋轉(zhuǎn)平臺表面.由于平移式表面實現(xiàn)難度往往較大,大部分實驗設(shè)計者都選擇轉(zhuǎn)輪式和旋轉(zhuǎn)平臺式運動表面模擬平移平面.在使用轉(zhuǎn)輪式運動表面模擬平移平面時,通常需要保證實驗景深內(nèi)表面與水平面間的偏移低至可以忽略[19],否則實驗表面被視為曲面將影響實驗結(jié)果.

圖1 液滴撞擊移動表面的實驗系統(tǒng)Fig.1 Experimental systems of droplet impacting moving surface

在表面選擇上,大部分研究人員使用鋁或不銹鋼作為表面材料.當(dāng)需要研究表面性質(zhì)時,多使用多晶硅[22]作為表面或使用聚四氟乙烯涂層[23]、SiO2涂層[24]等表面涂層進(jìn)行處理.液體選擇上多使用去離子水、乙醇、甘油、硅油及其水混合物.具體實驗材料與部分工況如表1 所示.

表1 液滴撞擊移動固體表面實驗條件Tab.1 Experimental conditions of droplet impacting moving solid surface

經(jīng)過對上述文獻(xiàn)整理發(fā)現(xiàn),現(xiàn)有的實驗對表面移動參量如Wes、Res等的覆蓋范圍較大,而針對液滴撞擊大部分參量如Wen等在較低的數(shù)值范圍內(nèi).

1.2 液滴撞擊移動液膜的實驗系統(tǒng)與方法

相較于液滴撞擊靜止表面的研究,液滴撞擊移動液膜的實驗研究起步較晚.在2011 年Alghoul 等[33]發(fā)表在液滴撞擊移動液膜方面的實驗研究成果.如圖2 所示的液滴撞擊移動液膜的實驗系統(tǒng)可以分為液滴發(fā)生系統(tǒng)、移動液膜系統(tǒng)、圖像捕捉系統(tǒng)和后處理系統(tǒng).液滴發(fā)生系統(tǒng)常使用毛細(xì)管針頭產(chǎn)生液滴,移動液膜系統(tǒng)常為一個橫置或斜置槽,用于提供穩(wěn)定的移動液膜,圖像捕捉和后處理系統(tǒng)則使用高速相機(jī)和匹配的計算機(jī)處理軟件,用于捕捉液滴撞擊現(xiàn)象和圖像處理.

圖2 常見液滴撞擊移動液膜實驗平臺[34]Fig.2 Common experiment platform of droplets impacting moving water film[34]

大部分實驗采用去離子水產(chǎn)生液滴和液膜,在考慮不同液體性質(zhì)時,可選擇使用不同比例的甘油和水的混合物[35],同時實驗所使用的液滴和液膜液體均保持相同.表2 為液滴撞擊移動液膜的部分實驗條件.經(jīng)過對上述文獻(xiàn)整理發(fā)現(xiàn),現(xiàn)有實驗針對表面液滴撞擊參量同樣在較低的數(shù)值范圍內(nèi),此外較多實驗針對的無量綱液膜厚度H較薄,針對H>1 的工況研究相對較少.

表2 液滴撞擊移動液膜實驗條件Tab.2 Experimental conditions of droplet impacting moving liquid film

2 液滴撞擊靜止和移動表面的行為與現(xiàn)象

液滴撞擊移動表面的過程是在液滴撞擊靜止表面的過程中加入一個垂直于液滴撞擊方向(或與液滴撞擊方向成一定角度)的表面運動,因此液滴撞擊移動表面現(xiàn)象與撞擊靜止表面現(xiàn)象具有一定的關(guān)聯(lián),在影響液滴撞擊靜止表面的變量同時,也會影響撞擊移動表面的過程.隨著表面運動形式的復(fù)雜化,液滴撞擊后的現(xiàn)象更復(fù)雜多樣,但是描述液滴撞擊靜止表面現(xiàn)象的一些方法仍用在撞擊移動表面上.

2.1 液滴撞擊靜止固體表面的行為與現(xiàn)象

目前研究人員針對液滴撞擊靜止固體表面的現(xiàn)象和機(jī)理已經(jīng)開展了十分深入的研究和討論.推薦讀者閱讀于2006 年Yarin 等[40]發(fā)表的文獻(xiàn),文獻(xiàn)中對液滴撞擊靜止表面的相關(guān)內(nèi)容進(jìn)行了詳盡的總結(jié).

Rioboo 等[41]于2001 年的研究中提出液滴撞擊靜止硬表面廣泛存在的6 種現(xiàn)象,如圖3 所示,分別為鋪展、噴濺、冠狀飛濺、收縮破裂、部分反彈和完全反彈.Mehdizadeh 等[42]還指出高速撞擊下液滴鋪展過程中有明顯的液指生成,而表面的粗糙度[13]、撞擊參數(shù)[43]、環(huán)境氣壓[44-45]、撞擊角度[46]等因素對液滴撞擊后的現(xiàn)象有明顯的影響.

圖3 液滴撞擊硬表面的典型現(xiàn)象[41]Fig.3 Typical phenomenon of droplet impacting solid surfaces [41]

在液滴撞擊過程中,常被科研人員關(guān)注的行為和現(xiàn)象有鋪展、飛濺和反彈.它們的發(fā)生閾值和具體現(xiàn)象受到表面性質(zhì)、液滴性質(zhì)、撞擊條件、環(huán)境條件等因素的影響.在鋪展過程中,研究人員主要關(guān)注的特征參數(shù)有鋪展因子、鋪展時間等.目前已經(jīng)有大量的實驗研究及模型,推薦讀者參考春江等[47]針對鋪展模型的詳盡整理.

在液滴撞擊表面飛濺的判別模型中,以K為參數(shù)的模型是最常見的一種模型,結(jié)構(gòu)通常比較簡單,并且可以包含黏性力、慣性力、表面張力等因素.研究者們通常以理論分析或經(jīng)驗公式的方式建立K的表達(dá)式,并結(jié)合實驗數(shù)據(jù)得到判別飛濺的臨界值Kc.一些常見的靜止表面(包括固體和液膜)發(fā)生飛濺的閾值如表3,其中列出部分受到較多引用的早期經(jīng)典模型和近年的新模型.不同模型中K的表達(dá)式和臨界值Kc都存在最佳適用范圍,當(dāng)撞擊中的部分參數(shù)例如表面潤濕性、撞擊韋伯?dāng)?shù)、液膜厚度等超出一定范圍,則當(dāng)前模型可能會失效;同樣當(dāng)撞擊的場景不同時,例如撞擊的表面是固體還是液膜、表面是否運動等,也會令現(xiàn)有模型出現(xiàn)誤差.針對運動表面的研究中,同樣需要建立新的K模型來判別飛濺閾值.

表3 部分液滴撞擊靜止表面飛濺閾值模型Tab.3 Partial splash threshold models for droplet impingement on stationary surface

值得注意的是有一種靜止表面較為特殊,此靜止表面與液滴下落方向之間存在夾角,即傾斜的靜止表面.如圖4 所示,液滴撞擊傾斜靜止表面的前期現(xiàn)象與液滴撞擊平移表面具有較高的相似性,當(dāng)2 個系統(tǒng)中液滴對表面保持相同的法向和切向速度,兩者的前期鋪展現(xiàn)象和飛濺現(xiàn)象幾乎等效[28].

圖4 液滴撞擊傾斜與移動固體表面鋪展對比[28]Fig.4 Comparison of liquid drop spreading between solid and moving surfaces [28]

2.2 液滴撞擊移動固體表面的行為與現(xiàn)象

在移動表面上,液滴的鋪展過程與靜止表面最明顯的不同在于鋪展的不對稱性,圖5 列出液滴在移動固體表面上的一些典型幾何不對稱現(xiàn)象.在運動的表面上,液滴一般會在沿表面運動方向的鋪展比垂直運動方向的鋪展尺度更大[54].液滴在移動表面的飛濺過程也受到一定程度的影響,例如在某些工況下,可能會在上下游兩側(cè)表現(xiàn)出不同的飛濺形式[27],液滴鋪展后上下游的飛濺行為會受到表面速度影響[27]、運動表面上液滴能夠發(fā)生飛濺的最低氣壓值顯著小于靜止表面[21]等.液滴的反彈過程同樣與靜止表面有著區(qū)別,例如表面運動速度達(dá)到一定閾值會引起液滴撞擊后反彈[55]、液滴與表面接觸時間縮短[24]等情況.

圖5 液滴撞擊移動表面的不對稱鋪展和飛濺Fig.5 Asymmetric spreading and splashing of liquid droplets impacting moving surfaces

由于液滴撞擊移動表面具有的特殊性,研究人員引入一些定義以描述移動表面上的現(xiàn)象.如圖6 所示,本研究統(tǒng)一定義液滴鋪展的最大寬度a為垂直于表面運動方向上的最大尺寸;液滴鋪展最大長度b為沿表面運動方向上的最大尺寸;液滴鋪展的上下游分別為表面運動方向所指的上游與下游,其分界線為液滴當(dāng)前鋪展最大寬度線.端點為側(cè)視圖中液滴最高點,端點偏移x為液滴端點位置與液滴鋪展幾何最大寬度位置間的偏移,固液接觸時間tc為液滴撞擊移動表面后到完全反彈間固液接觸的時間.在上述定義后,可以把部分實驗結(jié)果匯總[54-57]如表4 所示,初步描述表面運動速度對液滴撞擊后各現(xiàn)象的影響,其中“↑”為變量的增加、“↓”為變量的縮減.

表4 表面運動速度與液滴鋪展、飛濺、端點偏移和接觸時間的關(guān)系[54-57]Tab.4 Relationship between surface velocity and droplet spread, splash, apex offset and contact time[54-57]

圖6 液滴鋪展現(xiàn)象及對應(yīng)參量定義Fig.6 Droplet spreading phenomenon and corresponding parameter definition

在液滴連續(xù)撞擊移動表面的過程中,表面速度越高、液滴連續(xù)撞擊頻率越低,則液膜在表面的覆蓋沉積越平滑[58].當(dāng)需要較高的表面移動速度時,較低的頻率和較高的沖擊能量可以減小波動的生成,從而實現(xiàn)較為平滑的表面液膜覆蓋[59],而鋁液滴等高導(dǎo)熱系數(shù)的液滴,則會發(fā)生逐層凝固.液滴的凝固、鋪展和回縮的耦合作用使得凝固液滴表面出現(xiàn)一系列“L”形隆起[60].

旋轉(zhuǎn)固體表面作為移動表面的一種特殊類型,由于其離心效應(yīng)而存在一些有別于平移表面的特殊現(xiàn)象.根據(jù)液滴撞擊位置的不同,可以分為旋轉(zhuǎn)中心撞擊和非旋轉(zhuǎn)中心撞擊.液滴在旋轉(zhuǎn)中心撞擊后的現(xiàn)象有鋪展邊緣的毛細(xì)脊、液指的發(fā)生[61-62].液滴在非旋轉(zhuǎn)中心撞擊后的典型現(xiàn)象有如圖7 所示的“月牙形”的鋪展形狀、幾何不對稱的液指拉伸、內(nèi)側(cè)的波的形成和向外推進(jìn)等[29-32,63].

圖7 液滴撞擊旋轉(zhuǎn)表面非中心位置的不對稱鋪展、鋪展薄層上波的出現(xiàn)及液指現(xiàn)象Fig.7 Asymmetrical spreading of liquid droplets impacting noncentral position of rotating surface, appearance of waves on spreading layer and fingering phenomenon

Moghtadernejad 等[32]對液滴撞擊旋轉(zhuǎn)表面非中心位置的過程進(jìn)行初步的實驗研究,發(fā)現(xiàn)各項運動參數(shù)對液滴撞擊旋轉(zhuǎn)表面非中心位置的影響,其中較為特殊的有降低旋轉(zhuǎn)表面的角速度會使鋪展薄層上向外推移的波速度降低,提升液滴撞擊點的線速度會促進(jìn)二次液滴的形成等.本文沿用文獻(xiàn)中對各項變量的定義,即Vw為波速度,Dw為波寬度,Nf為液指數(shù), L wθ為周向鋪展尺寸,Lw為徑向鋪展尺寸.在實驗控制變量中, ω 為角速度,Rn為液滴撞擊位置半徑,實驗結(jié)果如表5 所示.

表5 液滴撞擊旋轉(zhuǎn)平臺非中心位置實驗結(jié)果[32]Tab.5 Experimental results of droplet impact on non-central position of rotating platform[32]

液滴撞擊移動固體表面的現(xiàn)象是由液滴撞擊靜止表面的基礎(chǔ)現(xiàn)象疊加表面運動干擾而產(chǎn)生的,具體現(xiàn)象跟表面的運動高度相關(guān).相比于靜止表面的相關(guān)實驗,液滴撞擊移動表面的實驗系統(tǒng)設(shè)計難度更高,其動力學(xué)過程更復(fù)雜,理論化推進(jìn)難度更高,因此液滴撞擊移動固體表面的研究文獻(xiàn)相較之下更少,模型體系的建立也更淺.

2.3 液滴撞擊靜止液膜的現(xiàn)象

液滴撞擊靜止液膜也是一個常被討論的物理學(xué)過程.通常來說,液滴撞擊液膜的主要現(xiàn)象可以分為懸浮、反彈、飛濺以及融合(如圖8),其中飛濺又分為冠狀飛濺和射流飛濺[33,64-66].

圖8 液滴撞擊液膜主要現(xiàn)象[33]Fig.8 Main phenomena of droplet impacting water film[33]

2.4 液滴撞擊移動液膜的現(xiàn)象

受到液膜移動的影響,液滴撞擊移動液膜和靜止液膜產(chǎn)生的水花形態(tài)會有明顯的變化,因此水花尺寸通常按照上下游區(qū)分.圖9 列舉了部分研究中的水花形態(tài)和劃分方法,圖中變量名為各文獻(xiàn)內(nèi)部定義,本研究中的變量名以本研究定義為準(zhǔn).

圖9 液滴撞擊移動液膜水花形態(tài)Fig.9 Spray morphology of droplet impact moving liquid film

液滴速度、表面速度之比以及液體性質(zhì)影響液滴撞擊移動液膜的現(xiàn)象[69].相對于撞擊靜止液膜,液滴撞擊移動液膜時初期產(chǎn)生的空腔、冠狀飛濺和射流呈不對稱狀態(tài)[34].液滴撞擊后形成冠狀飛濺和空腔,空腔會隨著液膜的流動逐漸趨于對稱,空腔縮小消失后產(chǎn)生向下游傾斜的射流,射流的根部會隨著液膜流動逐漸靠近尖端,使得射流與液膜表面的角度逐漸向垂直發(fā)展[33,70].隨著移動液膜流速的增加,液滴反彈的可能性也隨之增大[36].當(dāng)液膜流動速度超過液滴撞擊速度的10%后,撞擊濺起的液層出現(xiàn)彎曲[64].如果液滴撞擊在如圖10 所示的3 種不同形態(tài)的流動液膜區(qū)域[71-73]:在平液膜區(qū)域、毛細(xì)波區(qū)域和波峰區(qū)域上則會有不同的現(xiàn)象發(fā)生.其中平滑液膜最易發(fā)生飛濺,波峰上最難發(fā)生飛濺[38].當(dāng)液滴撞擊液膜的上方存在剪切氣流時,剪切空氣通過液滴變形和限制膜厚,對碰撞前的相互作用過程產(chǎn)生影響[74].

圖10 3 種不同形態(tài)的流動液膜區(qū)域[38]Fig.10 Three different flow liquid film regions[38]

Cheng 等[75-76]分別模擬液滴傾斜撞擊靜止液膜和液滴垂直撞擊移動表面上的液膜,發(fā)現(xiàn)傾斜撞擊靜止液膜和垂直撞擊移動表面上的液膜對液滴飛濺造成一定程度上具有相似性的影響,但是兩者無法完全等效.飛濺區(qū)域氣體的運動對飛濺液體尖端的形態(tài)產(chǎn)生擾動,使得2 個系統(tǒng)中液滴冠狀飛濺尖端的形態(tài)具有區(qū)別.與2.1 節(jié)中的傾斜和移動固體表面對比類似,不論是固體表面還是液膜,傾斜撞擊和移動表面撞擊之間都存在著一定的關(guān)聯(lián),但是也有著一定的差異.

3 液滴撞擊移動及旋轉(zhuǎn)表面的模型

在液滴撞擊固體表面過程的理論分析中,最常使用的手段是能量法.其原理是液滴撞擊過程中機(jī)械能、表面能、粘性耗散三者之和守恒[77].在靜止表面的研究中,研究人員曾利用能量法得到了最大鋪展面積等模型[78-81].

3.1 液滴撞擊移動固體表面的模型

對液滴撞擊移動固體表面的模型研究大體沿著液滴撞擊靜止固體表面的思路進(jìn)行,研究主要針對鋪展、飛濺、反彈3 個方面.在液滴鋪展階段,最早可以注意到的現(xiàn)象便是液滴鋪展的幾何不對稱性.Chen 等[26]在對去離子水撞擊聚四氟乙烯表面的實驗中,基于液滴撞擊移動表面的幾何不對稱現(xiàn)象,對實驗結(jié)果進(jìn)行擬合,同時指出“部分反彈—沉積”臨界態(tài)與“沉積—沉積分裂”臨界態(tài)下的液滴拉伸比例(即b/a)是不隨Wes變化的常數(shù).在實驗中當(dāng)該比值到達(dá)1.10 左右時,液滴撞擊從部分反彈變?yōu)殇佌梗划?dāng)比值到達(dá)1.46 左右時,從鋪展變?yōu)殇佌埂屏熏F(xiàn)象.為了對比靜止表面上液滴的最大鋪展面積,定義一個無量綱面積增長的計算式為

式中:XA為矢量綱面積增長,A為固液實際接觸面積.

實驗中該無量綱量與Wes呈線性關(guān)系,經(jīng)擬合得

為了深入理解液滴在運動表面鋪展的機(jī)理,Almohammadi 等[23]提出一個理論模型,如圖11 所示,其中r(t) 為在同等條件下液滴撞擊靜止表面的鋪展半徑[78-86].

圖11 液滴在平移固體表面鋪展模型[23]Fig.11 Model of droplet spreading on translational solid surface [23]

液滴撞擊移動表面的過程可以被分解為液滴由于撞擊產(chǎn)生的鋪展過程和表面運動的過程,同時提出2 項假設(shè):1)液體薄層鋪展關(guān)系不受到表面運動的影響;2)液體薄層鋪展的中心隨著表面的移動而移動.基于以上理解,將運動表面液體薄層鋪展分解為每個時刻液滴撞擊靜止表面的鋪展和該時刻鋪展中心隨表面上的位置移動.沿用靜止表面上的圓形鋪展模型,可以提出移動表面液滴鋪展邊緣隨時間變化的計算式為

在實際實驗中發(fā)現(xiàn)該模型與實驗結(jié)果仍有部分出入,這是由于液滴在撞擊過程中黏性力還會導(dǎo)致端點偏移,同時引入端點偏移量:

式中:C為偏移量修正數(shù).修正后的計算式為

式中:Vd為下游的液鋪展薄層收縮速度.上述模型公式與實驗結(jié)果[23]符合良好.

Buksh 等[20]為了使該模型可以覆蓋高張力和低張力的液體,提出改進(jìn)的模型為

式中:涉及到的3 個變量,a(t) 與b(t) 為模型的幾何參數(shù),取決于液滴撞擊參數(shù)和液滴性質(zhì),并同樣引入液滴端點的偏移量修正C.如圖12 所示,模型可以很好地同時描述低表面張力和高表面張力液體在達(dá)到最大鋪展尺寸前的鋪展過程.

圖12 橢圓形鋪展模型在不同表面運動速度下的表現(xiàn)[20]Fig.12 Effects of elliptical spreading model under different surface velocities [20]

針對液滴的飛濺,Bird 等[27]基于液滴薄層的不穩(wěn)定性提出飛濺閾值的模型.液滴薄層的動能表達(dá)式為

式中:h為液滴鋪展薄層的厚度,Vl為液滴鋪展薄層的延伸速度,L為液滴薄層的尺寸.由于變形引起的表面能增加可以表達(dá)為 σhL,而要引起飛濺,則需要前者遠(yuǎn)大于后者,即

由此給出了飛濺閾值計算式為

式中:K為模型參數(shù),其具體數(shù)值與表面性質(zhì)等因素有關(guān).模型在K=2.5 與K=5 700 的情況下,預(yù)測結(jié)果與實驗結(jié)果重合較好,可以簡化為表面速度為0,即液滴撞擊靜止表面的模式,并與早期實驗[47,87]相符.然而該模型有2 個局限性:1)該模型的實驗驗證僅限于乙醇一種液體,因此無法體現(xiàn)出液體黏度、表面張力等因素對飛濺現(xiàn)象的影響;2)該模型對液滴幾何不對稱性的描述僅有一個維度,即僅考慮了液滴薄層上游端點和下游端點,而無法描述整個液滴鋪展邊緣其他位置的飛濺情況.

為了解決上述局限,Almohammadi 等[19]以不同比例的水和甘油的混合物進(jìn)行實驗,在疏水表面液體薄層脫離表面的部分多于親水表面,且提高液體的黏性會降低飛濺閾值.如果設(shè)定一個如圖13 所示的位置角 ? 描述液滴鋪展邊緣的位置,則液體薄層邊緣從位置角 ?=0 (即上游的最遠(yuǎn)端)到 ?=180 (下游的最遠(yuǎn)端)的飛濺能力大幅遞減.描述不同位置角飛濺的模型為

圖13 液體薄層邊緣位置角 ? [19]Fig.13 Position angle of the liquid lammella edge ? [19]

式中:K1、K2為模型參數(shù),其具體數(shù)值與表面有關(guān).實驗中親水表面的K2= 4.53,疏水表面的K2=6.59.

Xu 等[44]提出液滴的飛濺可以通過改變環(huán)境氣壓而抑制后,研究人員發(fā)現(xiàn)環(huán)境氣體在液滴的飛濺過程中也起到作用.Hao 等[21]實驗研究移動表面環(huán)境氣壓對液滴撞擊飛濺產(chǎn)生的影響,通過降低氣壓可以抑制液滴上游的冠狀飛濺,在移動表面上飛濺的臨界氣壓小于靜止表面.根據(jù)單位長度液體薄層上升力與表面張力比值超過某閾值產(chǎn)生飛濺這一原理[88-89],給出基于升力的液滴撞擊移動表面的飛濺閾值模型為

式中:下標(biāo)ga 為氣體相關(guān)參數(shù),V1為液滴鋪展薄層的速度,Hf為液滴薄層邊緣厚度,Kl和Ku為與氣體雷諾數(shù)、氣體自由程等量相關(guān)的參數(shù),實驗中得到Kβ的平均值為0.137.

針對液滴的反彈,Povarov 等[55]觀察到移動表面邊界層的氣體升力對液滴撞擊表面反彈的作用,升力系數(shù)表達(dá)式為

式中: δ 為邊界層厚度.從氣體升力方面判斷反彈的機(jī)制,在升力系數(shù)達(dá)到一定條件時液滴可能發(fā)生反彈.此外在高表面速度時,液滴在接觸到固體表面前即受到邊界層氣流影響發(fā)生變形,環(huán)境氣體在液滴的反彈過程中起到明顯影響.

詹海洋[90]提出液滴在撞擊超疏水表面反彈的過程中存在滑移,而滑移受到表面剪切力與升力的共同作用,他提出滑移距離d表達(dá)式為

式中:Ks和KL分別為剪切力和升力的權(quán)重.同時發(fā)現(xiàn)移動超疏水表面上液滴的接觸時間相比靜止表面更短,這一現(xiàn)象也在Zhang 等[24]的研究中被觀察到.基于在靜止表面的模型[91-92],Zhang 等[24]提出撞擊移動超疏水表面的時間tM與液滴撞擊靜止表面的接觸時間tS之比可以描述為

3.2 液滴撞擊旋轉(zhuǎn)固體表面的模型

由于引入旋轉(zhuǎn)帶來的慣性力,液滴撞擊旋轉(zhuǎn)固體表面的物理過程相比于平移固體表面更為復(fù)雜.液滴不對稱鋪展與慣性力共同作用,產(chǎn)生一些在平移固體表面觀察不到的現(xiàn)象,如波的向外推進(jìn)等.液滴撞擊旋轉(zhuǎn)表面可以分為旋轉(zhuǎn)中心撞擊和非旋轉(zhuǎn)中心撞擊.在液滴撞擊旋轉(zhuǎn)固體表面中心時,隨著液體薄層的擴(kuò)散,可以觀察到邊緣不穩(wěn)定性以及液指的拉伸[61-62].Sahoo 等[61]基于線性穩(wěn)定性分析得到液指數(shù)預(yù)測計算式為

式中: ε 為液體的初始體積,Rc為代表液發(fā)生時的臨界鋪展半徑,由實驗數(shù)據(jù)決定.由于式(19)適用于完全潤濕液體,因此實驗中得到的數(shù)據(jù)與公式預(yù)測值有一定偏差.

當(dāng)液滴撞擊非旋轉(zhuǎn)中心位置時,會發(fā)生如波向外側(cè)推進(jìn)以及不對稱的液指拉伸等典型現(xiàn)象.Moghtadernejad 等[32]基于實驗對液體徑向鋪展尺寸 Lw 、二次液滴脫離液體薄層的時間td、液指的數(shù)量Nf、液指發(fā)生的時間tf進(jìn)行了擬合,分別得到

式中:round (·) 函數(shù)為指取函數(shù)內(nèi)部數(shù)值相鄰最近的整數(shù).在液滴撞擊旋轉(zhuǎn)表面非中心位置的研究中,分析波推進(jìn)速度等物理現(xiàn)象存在挑戰(zhàn),部分物理過程如圖14 所示,其中作用力分析較為復(fù)雜,需要借助數(shù)值模擬進(jìn)行輔助分析[32].

圖14 波推進(jìn)過程示意圖[32]Fig.14 Schematic diagram of liquid wave propulsion process[32]

3.3 液滴撞擊移動液膜的模型

液滴撞擊移動的液膜同樣會產(chǎn)生與撞擊靜止液膜不同的現(xiàn)象.當(dāng)前的研究主要基于液滴撞擊靜止液膜的成果,結(jié)合實驗提出液滴撞擊移動液膜的模型.部分模型可以簡化到液滴撞擊靜止液膜或固體表面,這些模型在靜止與移動表面上均具有適用性.

關(guān)于液滴反彈的條件如圖15 所示.由于氣層和液膜的厚度相對于切向長度來說很小,可以使用潤滑理論(lubrication theory)近似氣層和液膜中的流動[36].推導(dǎo)出液滴撞擊移動液膜過程中潤滑力[36]為

圖15 液滴下方氣流和液膜流動[36]Fig.15 Air flow and water film flow under droplets[36]

式中:下標(biāo) a 為環(huán)境氣體,Aga為邊界條件決定的參量.Che 等[36]提出潤滑力大于液滴重力在垂直于液膜運動方向的分量是液滴發(fā)生反彈的最小條件.

關(guān)于液滴的冠狀飛濺,在Roisman 等[93]的理論中,如果撞擊流動液膜的液滴在變形階段仍保持圓形,則在擴(kuò)展階段冠狀的底部可以近似為一個生長的圓,隨著液膜流向下游移動.Gao 等[35]發(fā)現(xiàn)上述理論在實驗中有所偏差,提出當(dāng)H×U?1時,才可以認(rèn)為水膜的流動對液滴的鋪展過程影響很小,此時液滴在擴(kuò)展過程中可以保持圓形,液滴撞擊移動液膜飛濺閾值計算式為

當(dāng)液膜厚度變?yōu)? 時,公式簡化為K=WeRe1/2,即液滴撞擊固體表面飛濺閾值公式[47].

冠狀飛濺的過程與液滴的速度、尺寸以及液膜流速有關(guān)[38].Yarin 等[94]通過擬合實驗數(shù)據(jù)得到了一個冠狀水花發(fā)展的經(jīng)驗計算式:

式中:dc為冠狀水花尺寸.Adebayo 等[38]則在此基礎(chǔ)上,基于努塞爾特厚度:

式中:g為重力加速度.得到了冠狀水花發(fā)展計算式為

式中:α 為渠道傾斜角度,Kθ為參數(shù),在實驗中得到的數(shù)值為1.23[38].

4 液滴撞擊移動及旋轉(zhuǎn)表面的數(shù)值模擬

液滴撞擊移動及旋轉(zhuǎn)表面的數(shù)值模擬計算按模擬尺度可以分為微觀尺度模擬、介觀尺度模擬和宏觀尺度模擬.微觀尺度模擬通常使用分子動力學(xué)模擬(molecule dynamics, MD).MD 模擬結(jié)果準(zhǔn)確與否的關(guān)鍵在于對系統(tǒng)內(nèi)的原子之間相互作用勢函數(shù)的選取[95],主要適用于納米級別微小液滴的場景,若計算場景內(nèi)分子數(shù)量過大,則會嚴(yán)重影響計算效率.介觀尺度模擬多使用格子玻爾茲曼方法(lattice Boltzmann method, LBM)[67],其模擬粒子的尺度比分子尺度大的多[96],可以模擬更大尺度的場景仍具有不俗的計算效率.當(dāng)前的LBM 方法在高雷諾數(shù)的多相流模擬和高液氣密度比的模擬中的數(shù)值穩(wěn)定性存在挑戰(zhàn),且該方法不適合穩(wěn)態(tài)計算[67].宏觀尺度模擬多適用界面追蹤法包括流體體積法(volume of fluid, VOF)和Level-Set 方法,可以模擬宏觀尺度下的液滴撞擊場景.以上兩者都是流體界面的追蹤方法相對容易實施,但是準(zhǔn)確性受到對流方程解的數(shù)值耗散的限制[96].由于VOF 方法和Level-Set 方法各有優(yōu)缺點,耦合2 種方法的模擬方法也開始出現(xiàn).針對液滴撞擊液膜的數(shù)值模擬中常用的耦合方法為耦合水平集流體體積法 (coupled level set and volume of fluid, CLSVOF)方法,該方法可以跟蹤不可壓縮兩相流的相界面,有效地解決了VOF 的計算收斂性、穩(wěn)定性和準(zhǔn)確性問題[97-98].

4.1 液滴撞擊水平移動固體表面的數(shù)值模擬

在微觀尺度模擬中,F(xiàn)ang 等[99]在MD 模擬中使用LJ 勢(Lennard–Jones potential)討論了不同的溫度、平面運動速度和射流角度下納米液滴對移動表面的撞擊.LJ 勢是經(jīng)典的連續(xù)對勢,認(rèn)為原子之間的相互作用是兩兩之間的作用[95].當(dāng)溫度上升時,液體的鋪展寬度會變大,液體的平均能量上升且能量波現(xiàn)象出現(xiàn),平板上液體原子的接觸角減小;隨著平板速度的上升,表面液膜的厚度下降,平板運動方向前部的接觸角更大.

2005 年Abascal 等[100]提出一種潛在的水凝聚相通用模型即TIP4P/2005 模型,該模型可以更針對性地描述水分子,并可以表現(xiàn)冰的微觀組成.Ritos 等[101-102]應(yīng)用TIP4P/2005 模型,設(shè)定硅、石墨和一種假想的超疏水材料進(jìn)行液滴撞擊的MD 模擬.硅表面上的模擬結(jié)果表明,接觸線上的分子位移受固體表面相互作用的影響很大,部分受表面運動引起的黏性耗散效應(yīng)的影響;對于石墨表面,前進(jìn)角和后退角與靜止接觸角接近,黏性耗散效應(yīng)可以忽略不計.與宏觀水滴潤濕動力學(xué)相反,納米尺度水滴撞擊過程中具有較大影響的是固液間作用,而黏性作用反而只起到次要影響,接觸角幾乎不受毛細(xì)數(shù)影響.

Heinz 等[103]表明與嵌入原子模型(embedded atom model, EAM)和密度泛函方法相比,12-6 LJ勢更加適用于大型模擬系統(tǒng)(106個原子), 并證明了此勢在銅上的可靠性.圖16 為Zhang 等[104]應(yīng)用12-6LJ 勢模擬納米液滴撞擊平移銅表面.

圖16 分子動力學(xué)模擬液滴撞擊移動固體表面[104]Fig.16 Molecular dynamics simulation of droplet impact on moving solid surface [104]

在相對滑動階段,水分子不僅沿平移表面移動,而且圍繞納米液滴的質(zhì)心旋轉(zhuǎn).在較高的表面平移速度下,納米液滴沿垂直于相對滑動方向擴(kuò)散2 次.納米液滴在表面運動方向的速度變化方程為

式中:Vx為液滴沿表面方向的運動速度.Wen的范圍為 7.41 ~6 6.67.

在宏觀尺度模擬中,Hou 等[105]通過VOF 方法獲得氣液界面,并利用ANSYS 框架下的凝固/融化算法對過冷液滴撞擊冷表面的凝固過程進(jìn)行模擬.該方法設(shè)定為過冷液滴接觸冷表面的瞬間轉(zhuǎn)化固液混合物,過冷度由調(diào)節(jié)混合物的參數(shù)表達(dá)[106-107].對橫向移動表面上液滴的形態(tài)演變進(jìn)行模擬,在與Zhang 等[24]實驗相似的條件下,得到了相似的結(jié)果.圖17 為對液滴撞擊移動冷表面的碰撞—結(jié)冰過程進(jìn)行的三維數(shù)值模擬,結(jié)果表明:液滴的初始凍結(jié)位置優(yōu)先出現(xiàn)在液滴的上游側(cè),液滴下游一側(cè)的凍結(jié)滯后于上游一側(cè),且滯后效應(yīng)受到冷表面溫度的強(qiáng)烈影響,同樣的凍結(jié)滯后現(xiàn)象也在孫志成[57]液滴撞擊移動冷表面的實驗中出現(xiàn).

圖17 VOF 方法模擬液滴撞擊移動固體表面的模型和初始條件[105]Fig.17 Model and initial conditions for simulating droplet impact on moving solid surface by VOF method[105]

Raman 等[108]則基于相位場不可壓縮格子玻爾茲曼方法進(jìn)行模擬,該模型使用2 個粒子分布函數(shù)來表達(dá)不可壓縮Navier-Stokes 方程和宏觀界面捕獲相場方程,該模型分別包含應(yīng)力和勢形式的動量方程分子間作用力項以及序參量的相場模型[108].模擬發(fā)現(xiàn)液滴在回縮狀態(tài)尾端會出現(xiàn)尾狀結(jié)構(gòu),該特點在Zhan 等[109]的液滴撞擊移動超疏水表面實驗中也可以被觀察到.Raman 的數(shù)值模擬進(jìn)一步揭示移動超疏水壁傳遞的動量會導(dǎo)致氣流從中心區(qū)域轉(zhuǎn)向下游方向,從而在液滴撞擊后的擴(kuò)散和收縮階段導(dǎo)致對稱性破壞的規(guī)律.

4.2 液滴撞擊移動液膜的數(shù)值模擬

Liang 等[110-111]基于CLSVOF 方法建立液滴撞擊移動液膜的三維模型,并對單液滴、雙液滴先后和同時撞擊等情況進(jìn)行數(shù)值模擬.流動液膜帶來的剪切力使界面和熱傳導(dǎo)的對稱性被打破,下游的飛濺被抑制,上游的飛濺被促進(jìn).Zhao 等[68]則發(fā)現(xiàn)液膜厚度的增加導(dǎo)致流入冠壁的液體量以及冠壁與液膜之間的角度增加,此外,非對稱冠的遷移距離隨著液膜的流動而增大.液滴在動態(tài)液膜上的擴(kuò)散速度等于上下游擴(kuò)散速度的平均值,其直徑與時間的1/2 次成正比.將實驗結(jié)果[35,112]與模擬結(jié)果相對照,如圖18 所示,證實模擬的準(zhǔn)確性,其中K值使用Gao 等[35]的計算方法.Wang等[113]學(xué)者則使用了二維模型以提升CLSVOF 方法的計算效率,同時模擬了雙液滴撞擊移動液膜的熱傳導(dǎo)過程,并發(fā)現(xiàn)非對稱分布特性導(dǎo)致壁面熱流不均勻.

圖18 液滴撞擊流動液膜的三維數(shù)值模擬與實驗對比[68]Fig.18 Comparison between three-dimensional numerical simulationand experimental of droplet impingement on flowing liquid film[68]

Xie 等[114]采用自適應(yīng)非結(jié)構(gòu)網(wǎng)格的控制體積有限元方法(novel control volume finite element method)研究環(huán)形兩相流中的三維液滴沉積過程.該數(shù)值框架包括用于界面捕捉的VOF 方法和用于自適應(yīng)非結(jié)構(gòu)網(wǎng)格上表面張力的力平衡連續(xù)介質(zhì)表面力模型.

在介觀模擬方法中,He 等[115]于1999 年改進(jìn)LBM 的相場模型.隨后Lee 等[116]于2005 年提出模擬高密度比不可壓縮兩相流的格子玻爾茲曼方程的穩(wěn)定離散化方法,以解決非理想氣體格子玻爾茲曼方程引發(fā)的數(shù)值不穩(wěn)定性.Raman 等[117]應(yīng)用該方法研究帶有液膜的移動壁對雙液滴撞擊行為的影響.在較大的液滴間距下,中心射流的形成會延遲;射流高度隨著液膜厚度的增加而先增加后減小;壁面運動速度增強(qiáng)上游射流,而抑制了中心射流和下游射流.Liu 等[118]基于LBM 方法并使用Cahn-Hilliard 式[119]捕捉液滴、液膜表面及周圍液體,發(fā)現(xiàn)控制液膜的厚度可以將液滴撞擊飛濺射流的高度控制在一個范圍內(nèi).He 等[67]采用具有可調(diào)表面張力項的格子玻爾茲曼方法(pseudopotential lattice Boltzmann method)研究液滴對運動液膜的沖擊,發(fā)現(xiàn)液膜的運動導(dǎo)致上下游冠的不對稱發(fā)展,上游和下游冠的不穩(wěn)定性隨著Ren和Wen的增加而增加.

5 結(jié)論與展望

液滴撞擊移動及旋轉(zhuǎn)表面(包括移動液膜)主要的特點是不對稱性.與撞擊靜止表面相比,液滴撞擊移動及旋轉(zhuǎn)表面的鋪展、飛濺、反彈等行為均具有幾何不對稱性.撞擊過程中的形態(tài)與表面的移動形式和速度相關(guān),反彈、飛濺等撞擊現(xiàn)象發(fā)生的閾值也受到表面移動形式及速度的影響.

在當(dāng)前的研究中,部分參數(shù)如撞擊參數(shù)、無量綱液膜厚度等的覆蓋范圍相對有限,同時旋轉(zhuǎn)表面等復(fù)雜情況缺乏數(shù)值模擬的輔助分析和更詳盡的實驗.這些問題值得相關(guān)學(xué)者進(jìn)一步關(guān)注和討論,同時也是液滴撞擊移動表面相關(guān)研究的重要發(fā)展方向之一.

當(dāng)前研究人員對于液滴撞擊移動及旋轉(zhuǎn)表面過程已經(jīng)開展了廣泛的研究,但是仍然存在一些可以進(jìn)一步探索的問題.

1)使用旋轉(zhuǎn)固體平面的研究已有許多,但是對于旋轉(zhuǎn)平面和平移平面間的區(qū)分仍然較少,旋轉(zhuǎn)帶來的離心力對液滴撞擊鋪展、飛濺等過程的影響機(jī)理尚不清晰.針對液滴撞擊旋轉(zhuǎn)固體表面開展更多研究,包括更大的撞擊速度范圍和表面旋轉(zhuǎn)速度范圍、不同的撞擊角度、不同液體和表面屬性等,為了建立更精確的鋪展及飛濺模型奠定基礎(chǔ).

2)研究液滴撞擊移動固體表面的模型大多從液滴撞擊靜止表面的研究基礎(chǔ)開始,而旋轉(zhuǎn)表面的部分特殊現(xiàn)象如波的推進(jìn)在靜止表面并無類似現(xiàn)象和模型可供參考,因此對液滴撞擊旋轉(zhuǎn)表面物理過程的研究需要更多理論和數(shù)值模擬數(shù)據(jù)支持.

3)針對液滴撞擊移動液膜的研究中,大部分研究設(shè)置的液滴和液膜所用的液體均為相同液體.液滴撞擊與自身不同液體組成的移動液膜過程可以作為進(jìn)一步研究的方向,液膜物性特征在撞擊中的影響.

4)當(dāng)前研究的液滴撞擊速度和液滴尺寸范圍仍較為局限,較高的撞擊韋伯?dāng)?shù)和較小尺寸液滴的相關(guān)研究仍較少.在很多撞擊移動表面和液膜的工程問題中,實際均涉及到微小液滴和高速液滴撞擊等問題.之后可以進(jìn)一步開展相關(guān)研究,完善對應(yīng)領(lǐng)域的理論體系.

5)遵循表面運動形式逐漸復(fù)雜這一推進(jìn)方向,在液滴撞擊固體表面和液膜表面2 個研究方向中,均可以進(jìn)一步探索復(fù)雜的表面運動形式帶給液滴撞擊現(xiàn)象的變化.學(xué)者們可以進(jìn)一步討論表面運動形式對液滴撞擊的影響,為了液滴撞擊移動表面理論體系提供支持,同時研究運動的表面對調(diào)控液滴撞擊后行為的可能性,例如設(shè)計表面的運動特性以減少或控制固液接觸時間、利用反彈和飛濺的不對稱性控制液滴撞擊后液體的運動方向、利用表面的運動特性調(diào)節(jié)液滴的飛濺閾值等,為后續(xù)工程應(yīng)用提供基礎(chǔ).

猜你喜歡
實驗模型研究
一半模型
記一次有趣的實驗
FMS與YBT相關(guān)性的實證研究
遼代千人邑研究述論
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計的漸近分布
視錯覺在平面設(shè)計中的應(yīng)用與研究
科技傳播(2019年22期)2020-01-14 03:06:54
做個怪怪長實驗
EMA伺服控制系統(tǒng)研究
3D打印中的模型分割與打包
主站蜘蛛池模板: 尤物视频一区| 欧美中文字幕在线二区| 免费av一区二区三区在线| 久久国产精品电影| 国产玖玖视频| A级毛片高清免费视频就| 2021国产精品自拍| 亚洲成人黄色在线| 99草精品视频| 色亚洲激情综合精品无码视频 | 福利一区在线| 国产欧美日韩在线一区| 欧美一区二区三区国产精品| 国产精品lululu在线观看| 亚洲精品桃花岛av在线| 久久久久国产一级毛片高清板| 成人免费一区二区三区| 在线观看国产黄色| 四虎成人免费毛片| 国产美女无遮挡免费视频| 欧美色视频网站| 国产高清在线观看91精品| 一区二区午夜| 欧美一级片在线| 网久久综合| 国产激情无码一区二区三区免费| 日韩精品无码免费专网站| 凹凸国产分类在线观看| 伊人久久久久久久| 极品性荡少妇一区二区色欲| 日本一本正道综合久久dvd| 2022国产无码在线| 久久精品中文字幕免费| 亚洲欧美激情小说另类| 久久国产精品波多野结衣| 九九这里只有精品视频| 国产真实乱人视频| 国产又色又爽又黄| 美女扒开下面流白浆在线试听 | 亚洲熟女中文字幕男人总站| 亚洲欧洲日产国码无码av喷潮| 国产地址二永久伊甸园| www亚洲精品| 99精品高清在线播放| 国产1区2区在线观看| 成年人国产视频| 99久久精品免费观看国产| 中文国产成人精品久久| 亚洲中文无码av永久伊人| 超碰aⅴ人人做人人爽欧美 | 国产91av在线| 国产精品视频导航| 香蕉在线视频网站| 青青青草国产| 国产一区二区免费播放| www.亚洲国产| 免费在线视频a| a免费毛片在线播放| 久久不卡精品| 成年人免费国产视频| 成人在线第一页| 国产精品免费福利久久播放| 91久久精品国产| 欧美成人精品在线| 乱系列中文字幕在线视频| 久久天天躁狠狠躁夜夜躁| 中文字幕 欧美日韩| 狼友视频一区二区三区| 亚洲一级无毛片无码在线免费视频 | 欧美日韩精品一区二区在线线| 一级一级特黄女人精品毛片| 在线亚洲天堂| 国产精品人人做人人爽人人添| 亚洲无码精品在线播放| 国产在线91在线电影| a天堂视频| 久久亚洲高清国产| 精品人妻AV区| 亚洲福利片无码最新在线播放| 国产人成在线视频| 午夜福利免费视频| 91久久天天躁狠狠躁夜夜|