• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    液氮滴撞擊壁面相變行為的數(shù)值研究*

    2019-12-24 08:21:54趙可佘陽梓蔣彥龍秦靜張振豪
    物理學(xué)報(bào) 2019年24期
    關(guān)鍵詞:氣膜潤濕液膜

    趙可 佘陽梓 蔣彥龍 秦靜 張振豪

    (南京航空航天大學(xué), 飛行器環(huán)境控制與生命保障工業(yè)和信息化部重點(diǎn)實(shí)驗(yàn)室, 南京 210000)

    采用Level Set-VOF方法建立單液氮滴撞擊壁面的數(shù)值模型, 探索壁面潤濕性(30°—150°)、撞擊速度(0.1和1.6 m/s)及壁面溫度(300—500 K)對液滴撞壁演化過程中相變行為的影響, 并理論推導(dǎo)了氣膜生長數(shù)學(xué)模型.結(jié)果表明:增強(qiáng)壁面潤濕性、提高撞擊速度有利于液滴沿徑向鋪展, 從而增大了換熱面積并降低熱阻, 使換熱性能得到顯著提升; 提高壁面溫度增大了換熱溫差, 熱流密度隨之上升; 三相接觸線處熱阻較小導(dǎo)致邊緣處熱流密度高于中心處, 不同潤濕壁面上熱流分布的差異性因初始速度的增大而縮小, 呈現(xiàn)明顯的速度效應(yīng); 在膜沸騰區(qū), 傳熱過程主要集中在撞擊初期, 氣膜是主要換熱熱阻; 基于質(zhì)量守恒和能量守恒建立氣膜生長數(shù)值模型, 模型預(yù)測結(jié)果與本文模擬結(jié)果和其他研究結(jié)果非常吻合.

    1 引 言

    噴霧冷卻技術(shù)作為一種高效的熱控制技術(shù)引起眾多學(xué)者的廣泛關(guān)注, 其冷卻過程涉及單液滴撞擊壁面、液滴與液滴、液滴與液膜之間的相互作用,屬于復(fù)雜的多相流問題[1,2].

    由于液氮具有冷卻溫度低、成本低和無殘留等優(yōu)點(diǎn)成為低溫控制領(lǐng)域(低溫風(fēng)洞、低溫試驗(yàn)條件的開發(fā)、大功率電子設(shè)備的熱管理、食品速凍以及低溫手術(shù)等)最理想的工質(zhì)之一[3-6].然而液氮沸點(diǎn)較低, 霧化液滴與壁面接觸后瞬間蒸發(fā)并在液滴底部和壁面之間產(chǎn)生一層連續(xù)的氣膜, 發(fā)生膜沸騰, 即出現(xiàn)Leidenfrost現(xiàn)象[7,8], 熱阻較大的氣膜往往會導(dǎo)致?lián)Q熱效果劇烈下降.研究低溫液滴撞擊壁面的鋪展動(dòng)態(tài)及相變行為過程有助于深入理解低溫噴霧冷卻技術(shù)的換熱機(jī)理, 尤其是發(fā)生Leidenfrost效應(yīng)時(shí)的熱傳輸機(jī)理.

    液滴撞壁是跨尺度、多物理場耦合控制的強(qiáng)瞬變過程, 其影響因素除液滴自身的熱物理性質(zhì)外,還包括熱沉面結(jié)構(gòu)或理化性質(zhì)等方面[9].許多學(xué)者利用理論分析、實(shí)驗(yàn)研究及數(shù)值模擬手段開展了大量細(xì)致的研究工作.Karl和 Frohn[10]研究了在Leidenfrost溫度以上液滴撞擊熱壁的力學(xué)行為,根據(jù)實(shí)驗(yàn)結(jié)果建立的關(guān)系式可用于改善兩相流的數(shù)值模擬, 并且從理論方面利用動(dòng)量損失的相關(guān)性, 推導(dǎo)出液滴最大變形的理論近似值, 其測量結(jié)果與其他文獻(xiàn)報(bào)道的結(jié)果非常吻合.Scheller和Bousfield[11]通過實(shí)驗(yàn)研究撞擊條件(初始直徑、初始速度)和壁面條件(壁面種類)的改變對撞擊行為的影響, 得出液滴撞擊壁面擴(kuò)展過程不發(fā)生回縮時(shí)無量綱潤濕長度最大值的經(jīng)驗(yàn)公式.

    在液滴物性方面, 劉海龍等[12]采用碳納米管、石墨烯、納米石墨粉制備三種穩(wěn)定的納米流體, 利用顯微高速數(shù)碼攝像技術(shù)捕捉上述三種納米流體撞擊固體壁面的動(dòng)態(tài)過程, 結(jié)果表明, 三種納米顆粒均能促使基液表現(xiàn)出顯著的剪切變稀特性; 流體的剪切黏度在液滴撞擊壁面的擴(kuò)展過程中起到了重要的作用; 液滴的無量綱高度和鋪展因子的變化與納米流體剪切黏度的變化呈負(fù)相關(guān); 在液滴的最初鋪展階段, 慣性力起主導(dǎo)作用, 隨著液滴撞擊速度的增大液滴的擴(kuò)展范圍和速度亦增大.沈勝強(qiáng)等[13]采用高速攝像儀觀察了水和乙醇液滴撞擊不銹鋼壁面的蒸發(fā)過程, 分析了液滴撞擊加熱壁面后的蒸發(fā)特性參數(shù), 實(shí)驗(yàn)得到水滴的平均熱流密度為 0.014—0.110 W·mm—2.

    在熱沉面結(jié)構(gòu)或理化性質(zhì)等方面, 已有研究表明熱沉面潤濕性對液滴的蒸發(fā)和二次霧化影響較大, ?ikalo等[14]研究發(fā)現(xiàn)已有的動(dòng)態(tài)接觸角經(jīng)驗(yàn)?zāi)P?Hoffman-Voinov-Tanner模型)不能很好地預(yù)測動(dòng)態(tài)接觸角的變化, 特別是在毛細(xì)數(shù)較大的情況下.梁剛濤等[15]實(shí)驗(yàn)研究低速液滴撞擊潤濕球面的現(xiàn)象, 發(fā)現(xiàn)液滴撞擊潤濕球面上發(fā)生反彈、局部反彈和鋪展等現(xiàn)象, 考察了黏度對撞擊過程的影響, 定量討論了液滴鋪展特征參數(shù)隨撞擊速度、球體直徑和黏度的變化規(guī)律.另外, 壁面溫度是影響相變過程的直接因素.Tran等[16]利用高速攝像機(jī)拍攝水滴 (0.5 ≤We≤ 500)撞擊高溫光滑壁面(250 ℃ ≤T≤ 560 ℃)過程呈現(xiàn) 3 個(gè)不同的機(jī)制:接觸沸騰、溫和膜沸騰、霧化膜沸騰.溫和膜沸騰、霧化膜沸騰機(jī)制下的液滴無量綱最大擴(kuò)展因子與We2/5成正比, 隨著We增加TDL也隨之增大;從溫和膜沸騰到霧化膜沸騰的過渡溫度與液膜內(nèi)的沸騰氣泡有關(guān), 并且隨著We的增加過渡溫度降低.

    綜上所述, 國內(nèi)外學(xué)者對液滴撞擊壁面的實(shí)驗(yàn)和數(shù)值研究, 多集中于液滴撞壁動(dòng)力學(xué)特性, 如三相接觸線行為、動(dòng)態(tài)接觸角變化、自由表面流動(dòng)等,而往往忽略液滴撞擊壁面的相變行為, 尤其是缺乏低溫液滴撞壁后的熱傳輸機(jī)制研究.

    本文利用復(fù)合Level Set-VOF相界面追蹤的方法, 建立單個(gè)液氮液滴撞擊不同潤濕壁面的計(jì)算模型, 并進(jìn)行了模型的可靠性驗(yàn)證.進(jìn)一步從微觀角度探究壁面熱流密度分布和氣膜厚度隨壁面潤濕性、初始速度和壁面溫度等因素的變化規(guī)律, 進(jìn)而揭示Leidenfrost效應(yīng)下低溫液滴撞擊壁面時(shí)的運(yùn)動(dòng)形態(tài)和蒸發(fā)換熱耦合機(jī)制.

    2 計(jì)算模型及模型驗(yàn)證

    2.1 復(fù)合Level Set-VOF相界面追蹤方法

    采用復(fù)合Level Set-VOF方法追蹤不可壓縮兩相流相界面, 可有效解決VOF和Level Set兩種氣液兩相流動(dòng)相界面追蹤方法存在的計(jì)算收斂性、穩(wěn)定性及準(zhǔn)確性不高等問題[17].復(fù)合Level Set-VOF相界面追蹤方法計(jì)算過程包括以下幾個(gè)步驟:Level Set函數(shù)和 VOF 函數(shù)初始化、流動(dòng)控制方程求解、Level Set函數(shù)和VOF 相函數(shù)(?函數(shù)和F函數(shù))的對流輸運(yùn)方程求解、相界面重構(gòu)以及函數(shù)再次初始化等.而氣液相界面重構(gòu)和函數(shù)再次初始化使得函數(shù)滿足質(zhì)量守恒特性, 并且能夠更加精確地計(jì)算表面張力.圖1為復(fù)合Level Set-VOF氣液兩相流動(dòng)計(jì)算流程圖, 虛線部分為氣液相界面追蹤和壁面潤濕模型的計(jì)算流程.

    復(fù)合Level Set-VOF法通過?函數(shù)和F函數(shù)共同構(gòu)造相界面, 采用分段線性界面重構(gòu)思想, 在法向方向上移動(dòng)界面使得單元液態(tài)區(qū)域面積比率和F函數(shù)值相匹配, 單元中心到相界面的垂直距離由割線法迭代求解得到.

    界面法向量n計(jì)算方法為

    界面曲率κ(?) 的表達(dá)式為

    復(fù)合Level Set-VOF方法的相界面的追蹤主要通過求解相函數(shù)對流輸運(yùn)方程來實(shí)現(xiàn).

    式中u為速度.?函數(shù)經(jīng)過對流輸運(yùn)方程求解后,將不再保持距離函數(shù)的性質(zhì), 因此必須對?函數(shù)重新初始化, 再次重新初始化主要包括?函數(shù)符號以及?函數(shù)值的確定.

    ?函數(shù)值是單元中心到相界面的最小距離d,而?函數(shù)的符號由F函數(shù)值決定:

    當(dāng)F< 0.5 時(shí), 單元中心在氣態(tài)區(qū)域,?函數(shù)為正; 當(dāng)F> 0.5 時(shí), 單元中心在液態(tài)區(qū)域,?函數(shù)為負(fù); 當(dāng)F= 0.5 時(shí), 表示相界面,?=0 ; s ign(·) 是符號函數(shù).

    圖1 復(fù)合 Level Set-VOF 相界面追蹤方法計(jì)算流程圖Fig.1.Coupled Level Set-VOF phase interface tracking method calculation flow chart.

    2.2 控制方程

    1)連續(xù)性方程和動(dòng)量方程

    式中,ρ為密度,p為壓強(qiáng),τ為黏性應(yīng)力張量,g為重力加速度,Fs為源項(xiàng).

    采用 CSF (continuum surface force)模型對表面張力進(jìn)行求解, 在CSF模型中附加的表面張力通過在動(dòng)量方程中添加源項(xiàng)的方式實(shí)現(xiàn):

    式中,σ為表面張力系數(shù),H(?) 為 Heaviside函數(shù).

    Heaviside函數(shù)表達(dá)式如下:

    式中,h為網(wǎng)格大小.黏性應(yīng)力張量τ=2μ(?)S,應(yīng)變率張量S表達(dá)式為

    不同區(qū)域ρ和μ可由Heaviside函數(shù)計(jì)算:

    2)能量方程

    式中,cp為計(jì)算單元內(nèi)定壓熱容,λ為計(jì)算單元內(nèi)熱導(dǎo)率,T為計(jì)算單元內(nèi)溫度.

    3)壁面黏附

    模型中還考慮了壁面黏附作用(wall adhesion),如圖2所示, 通過給壁面接觸角賦值的形式來定義壁面的潤濕性, 當(dāng)液滴在固體表面上鋪展時(shí), 定義壁面靜態(tài)接觸角為θ, 用于調(diào)整壁面附近單元表面的法向量.nw和τw分別為單位法向量和壁面切向量, 壁面黏附作用被定義在表面向量中:

    圖2 壁面黏附Fig.2.Wall adhesion diagram.

    2.3 蒸發(fā)模型

    液氮沸點(diǎn)較低, 當(dāng)液氮滴撞擊常溫壁面時(shí), 液滴撞壁過程存在相變, 需考慮液滴撞壁后的中相變行為.本文采用應(yīng)用最廣泛的Lee[18]模型, 其具體表達(dá)式如下.

    蒸發(fā)過程,

    冷凝過程

    蒸發(fā)的氣體質(zhì)量作為質(zhì)量源項(xiàng)離散后加載到控制方程中.計(jì)算涉及的物性參數(shù)如密度、熱導(dǎo)率、比熱容等均取兩相數(shù)學(xué)的平均值.當(dāng)蒸發(fā)發(fā)生在交界面時(shí), Hertz Knudsen方程可以表達(dá)為

    式中,M為氣體分子質(zhì)量,R為理想氣體常數(shù).Clapeyron-Clausius建立了壓力和飽和溫度的狀態(tài)關(guān)系:

    式中L為汽化潛熱.當(dāng)溫度和壓力接近飽和狀態(tài)時(shí), (18)式可表達(dá)為

    則Hertz Knudsen方程可表達(dá)為

    為了正確添加源項(xiàng), 將(20)式乘以單位體積內(nèi)的表面積數(shù)Ai=6αgαl/Dsm,Dsm為 Sauter平均直徑, 則蒸發(fā)源項(xiàng)可表達(dá)為

    2.4 模型設(shè)置及網(wǎng)格無關(guān)性驗(yàn)證

    圖3為單液氮液滴撞擊壁面和發(fā)生Leidenfrost蒸發(fā)后的形貌演變圖, 其計(jì)算域?yàn)?3 mm × 6 mm的二維四邊形均勻網(wǎng)格.定義初始直徑為D0, 初始速度為U0, 方向垂直壁面向下, 環(huán)境壓力為0.1 MPa,接觸直徑為Dt, 最大鋪展直徑為Dmax,δt為中心氣膜厚度, 考慮重力g的影響, 環(huán)境氣體為氮?dú)? 頂部及兩側(cè)為壓力出口邊界, 出口邊界壓力梯度為零, 底部為無滑移壁面.

    圖3 初始時(shí)刻液滴撞擊壁面幾何模型及邊界條件Fig.3.Model of droplet impact on wall at initial time.

    實(shí)際噴霧液滴的幾何形狀較為復(fù)雜, 本文僅研究完美球形液滴的撞壁過程, 且初始時(shí)刻液滴底部與壁面相切; 由于液滴與環(huán)境及壁面之間的溫差在 222—422 K, 此時(shí)輻射換熱較弱, 因此不考慮低溫液氮液滴與環(huán)境、壁面之間的熱輻射效應(yīng); 模擬液滴撞擊速度較低, 環(huán)境壓力為常壓, 可將氣液兩相均視為不可壓縮流體, 且忽略環(huán)境氣體與液滴間的剪切力.

    利用軟件Fluent 17.0求解數(shù)值模型, 控制方程采用有限體積法進(jìn)行離散, 壓力與速度耦合采用Coupled方法, 壓力求解采用PRESTO方法,采用QUICK算法對Level Set方程和對流項(xiàng)進(jìn)行求解, 采用 Geometric-Reconstruction Scheme 求解VOF方程, 動(dòng)量和能量方程采用二階迎風(fēng)格式,時(shí)間步長 Δt= 10—6s, Δt內(nèi)迭代次數(shù)為 20.計(jì)算相關(guān)參數(shù)見表1.

    表1 相關(guān)參數(shù)Table 1.Related parameters.

    網(wǎng)格無關(guān)性驗(yàn)證分別采用網(wǎng)格數(shù)200 × 400,250 × 500, 300 × 600 及 350 × 700 的四邊形結(jié)構(gòu)化網(wǎng)格, 具體工況為:撞擊初速度為 0.1 m/s, 液滴直徑為 0.5 mm, 壁面溫度為 300 K, 壁面靜態(tài)接觸角為90°.對比分析四種網(wǎng)格的鋪展因子β隨時(shí)間的變化關(guān)系, 結(jié)果如圖4所示.其中液滴鋪展因子定義為液滴接觸直徑Dt與液滴初始直徑D0的比值:

    從圖4可以看出, 隨著網(wǎng)格數(shù)密度增大, 鋪展因子β逐漸提高, 當(dāng)網(wǎng)格數(shù)增加至300 × 600后,鋪展因子β增加趨勢平緩, 網(wǎng)格數(shù)對鋪展因子影響較小.綜合考慮計(jì)算精度及計(jì)算效率, 本文模擬采用 300 × 600 網(wǎng)格數(shù).

    圖4 網(wǎng)格無關(guān)性驗(yàn)證Fig.4.Verification of grid independence.

    2.5 模型驗(yàn)證

    為了對模型進(jìn)行驗(yàn)證, 本文選取文獻(xiàn)[21]中的兩種潤濕壁面上 (靜態(tài)接觸角θ= 57°, θ = 25°)的液滴撞擊實(shí)驗(yàn)數(shù)據(jù)進(jìn)行模擬計(jì)算(工況:壁面溫度為 500 ℃, 液滴初始直徑為 2 mm, 液滴初始溫度為 20 ℃, 液滴初始速度為 2.05 m/s), 并將計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果進(jìn)行對比(圖5(a)).對比文獻(xiàn)[21]實(shí)驗(yàn)的結(jié)果可以看出, 隨著時(shí)間推移, 鋪展因子逐漸增大, 且在相同撞擊速度和初始直徑的情況下,接觸角較小的壁面上鋪展因子要高于接觸角較大的壁面, 兩種工況下模擬結(jié)果與實(shí)驗(yàn)結(jié)果高度吻合;另外, 為進(jìn)一步驗(yàn)證模型的準(zhǔn)確性, 選取文獻(xiàn)[22]中一組數(shù)據(jù)進(jìn)行液滴撞壁后的形貌演變對比(工況:壁面溫度 384 ℃, 液滴初始直徑為 2 mm, 液滴初始溫度為 20 ℃,We= 2), 對比結(jié)果見圖5(b),可以看出液滴撞擊光滑壁面后先后經(jīng)歷了鋪展、回縮和反彈等過程, 模擬液滴撞壁后的運(yùn)動(dòng)行為與實(shí)驗(yàn)結(jié)果較吻合.

    Chandra和Aziz[23]設(shè)計(jì)了一種低溫液滴發(fā)生裝置, 利用該裝置實(shí)驗(yàn)觀測到直徑約為1.9 mm的液氮撞擊銅面后發(fā)生了Leidenfrost現(xiàn)象, 液滴撞擊壁面的速度由下落高度控制, 實(shí)驗(yàn)發(fā)現(xiàn)液氮滴在距壁面27 mm處自由下落撞擊銅面后在其底部形成氣膜, 經(jīng) 8.2 s 后液氮滴才蒸發(fā)完全, 而本文模擬的最大撞擊時(shí)間尺度tmax=D0/U0= 5 ms, 即液滴完全蒸發(fā)時(shí)間尺度遠(yuǎn)大于液滴撞擊時(shí)間尺度;另外, 實(shí)際液氮噴霧冷卻過程中, 噴嘴距熱沉面距離較短, 液氮滴初始速度較大, 液氮滴在撞擊壁面之前不會因?yàn)榄h(huán)境溫度過高而完全蒸發(fā).通過上述分析可看出本文的數(shù)值模擬方法及建立具有一定初速度和初始直徑的液氮滴撞擊水平壁面的物理模型是可行的.

    圖5 實(shí)驗(yàn)結(jié)果與模擬結(jié)果對比 (a)鋪展因子隨時(shí)間變化; (b)液滴形貌對比Fig.5.Comparison of experimental and simulated results:(a) Spreading factor changes with time; (b) comparison of droplet morphology evolution.

    3 結(jié)果討論與分析

    3.1 壁面潤濕性和初始速度的影響

    液滴撞擊壁面的鋪展行為是潤濕現(xiàn)象的一種.壁面的潤濕性是影響表面沸騰和凝結(jié)性能的重要參數(shù).一般壁面潤濕性可通過靜態(tài)接觸角來表征:親水壁面 (θ< 90°)、疏水壁面 (θ> 90°)和超疏水壁面 (θ> 150°).

    圖6為0.5 mm的液滴以 0.1和1.6 m/s分別撞擊親水壁面 (θ= 30°)、普通壁面 (θ= 90°)和疏水壁面 (θ= 150°)的鋪展行為隨無量綱時(shí)間 (t*=tU0/D0)的變化趨勢.分析圖6(a)可知, 當(dāng)初速度為 0.1 m/s,θ= 30°時(shí), 撞擊初期液滴底部與壁面直接接觸, 接觸區(qū)域發(fā)生劇烈相變, 從而在液滴底部形成空穴(t*= 0.02), 與此同時(shí)液滴動(dòng)能向表面能和黏性耗散轉(zhuǎn)化, 三相接觸線處的合力指向外側(cè), 三相接觸線不斷外延, 固液接觸面積增大, 促使更多液相向氣相轉(zhuǎn)移, 空穴體積逐漸增大; 在t*=0.1時(shí), 液膜的壓迫使空穴內(nèi)部蒸氣沿徑向逃逸,液滴底部初步形成連續(xù)氣膜.親水壁面具有較大的表面能有利于液滴撞壁后液膜邊緣沿徑向向外鋪展, 而液體的表面張力則阻礙這種趨勢的發(fā)生, 導(dǎo)致液膜頸部出現(xiàn)極薄液膜區(qū), 最終伴隨著換熱過程的進(jìn)行, 在頸部發(fā)生斷裂 (t*= 0.16).

    在普通和疏水壁面上, 撞擊前期液滴在壁面上的鋪展?fàn)顟B(tài)以及連續(xù)氣膜的生長與親水壁面較為相似, 但由于接觸角增大, 壁面對液體吸附力減弱,液滴在壁面上的鋪展速度相對減慢, 鋪展面積相對減小, 液滴邊緣難以發(fā)生破碎; 尤其是當(dāng)液滴撞擊疏水壁面時(shí), 液體在其表面難以潤濕, 黏性耗損減少, 使得液膜表面能較大, 從而表現(xiàn)出在t*= 6.4時(shí)液膜開始回縮, 液體向液膜中心聚集.

    液滴在回縮過程中表面能的變化可表示為

    式中,A1和A2表示液滴回縮過程中的潤濕面積,可以看出壁面潤濕性直接影響液滴形貌演變過程中表面能變化的大小, 若潤濕面積固定, 隨著接觸角增大(0°—180°), 液滴收縮需要克服的表面能減小[24], 因此液滴碰撞疏水壁面后更容易發(fā)生回縮,導(dǎo)致液體向中心聚集.

    圖6 初始速度為 (a) 0.1 m/s和 (b) 1.6 m/s的液滴撞擊不同潤濕壁面的形貌演變 (D0 = 0.5 mm)Fig.6.Morphology evolution of droplets impinging on different wetted walls with initial velocity (D0 = 0.5 mm):(a) 0.1 m/s;(b) 1.6 m/s.

    分析圖6(b)可知, 當(dāng)速度為 1.6 m/s時(shí), 液滴在三種潤濕性壁面上快速鋪展, 并且在邊緣處發(fā)生二次小液滴飛濺(t*= 0.16), 原因是液滴速度增大后We隨之增大, 慣性力起主導(dǎo)作用并能夠完全克服液體的表面張力和黏性應(yīng)力, 因此在邊緣處產(chǎn)生破碎的小液滴, 體積較小的二次液滴受到周圍氣體的阻滯和液膜內(nèi)部徑向流動(dòng)的作用致使其斜向上運(yùn)動(dòng); 在t*= 1.1 時(shí), 三種潤濕壁面連續(xù)氣膜均初步形成, 伴隨液滴在三種壁面上的快速鋪展, 液體在壁面上的覆蓋面積逐漸增大; 在t*= 3.2 時(shí), 鋪展成極薄的液膜, 由于表面張力和慣性力的拮抗作用導(dǎo)致頸部液滴不斷剝落, 此時(shí)二次液滴粒徑較大, 周圍氣體的阻滯作用相對減弱, 剝落的小液滴受液膜徑向流動(dòng)作用而不斷外移; 在t*= 6.4時(shí),液滴破碎完全形成鏈狀小液滴群.重要的是, 與低速液滴相比, 液滴在三種潤濕壁面上形貌演變的差異性逐漸減弱, 在壁面的鋪展面積逐漸接近, 且潤濕面積均大于同等條件下的低速液滴, 表現(xiàn)出明顯的速度效應(yīng).

    圖7為0.5 mm的液滴以 0.1和1.6 m/s撞擊不同潤濕壁面的熱流分布, 壁面溫度均為300 K.圖8對應(yīng)速度為0.1 m/s的液滴撞壁面后的溫度分布, 靜態(tài)接觸角θ= 30°.

    結(jié)合圖6(a)液滴的形貌演化和圖7(a)熱流分布圖可以看出, 當(dāng)初始速度為 0.1 m/s 時(shí), 在撞擊初期(t* = 0.01), 液滴底部局部與壁面有短暫接觸, 壁面中心接觸區(qū)域傳熱方式以熱傳導(dǎo)為主, 因此三種潤濕壁面的熱流分布集中在壁面中心, 疏水、普通和親水壁面中心的熱流峰分別為2 × 106,4 × 106和 5 × 106W/m2, 熱流密度隨著潤濕性增強(qiáng)而增大; 隨后 (t*= 0.02—0.1), 液滴撞壁后在壁面上沿徑向逐漸鋪展并伴隨著劇烈相變, 液滴底部空穴內(nèi)部壓力逐漸增大直至將液膜托起, 發(fā)生Leidenfrost效應(yīng), 熱流密度急劇下降, 親水壁面上的最高熱流峰值僅為初期的10%.由圖8(a2)可以看出, 在t* = 0.1時(shí)氣膜的存在顯著增大了熱阻,壁面中心處A點(diǎn)溫度梯度僅為370 K/mm, 遠(yuǎn)小于邊緣處B點(diǎn)的2469 K/mm, 這對于噴霧冷卻是非常不利的.由圖6(a)的形貌圖可以看出, 壁面上鋪展的液膜中心和四周凸起, 壁面各處熱阻不同,反映在圖7(a2)中為熱流分布不均勻.當(dāng)t* = 0.3時(shí),整體熱流峰值顯著下降, 最高峰值不足35 W/m2,分析圖8(a3)發(fā)現(xiàn), 在撞擊后期氣膜充分生長, 厚度顯著增大, 壁面中心處A點(diǎn)與邊緣處B點(diǎn)的溫度梯度較為接近約為222 K/mm, 可見在撞擊后期換熱被嚴(yán)重削弱.另外, 隨著接觸角的增大, 液滴收縮需要克服的表面能變小, 因此液滴碰撞疏水壁面后更容易發(fā)生回縮, 導(dǎo)致液體向中心聚集, 中心處換熱加強(qiáng), 從而呈現(xiàn)出熱流峰寬逐漸減小且熱流向壁面中心轉(zhuǎn)移的趨勢.

    圖8 不同時(shí)刻的溫度分布 (θ = 30°, U0 = 0.1 m/s, Tw = 300 K)Fig.8.Temperature distribution at different times (θ = 30°, U0 = 0.1 m/s, Tw = 300 K).

    分析圖7(a)和圖7(b)發(fā)現(xiàn), 當(dāng)撞擊初速度增至 1.6 m/s時(shí), 壁面熱流密度在不同時(shí)期 (t*=0.01—0.3)均大幅提升; 在撞擊初始階段 (t*=0.01), 液滴在三種潤濕壁面上的形貌演變較為相似, 壁面熱流曲線幾乎重合, 并且熱流集中在液滴與壁面接觸區(qū), 且熱流峰面積明顯大于同期的低速液滴; 當(dāng)t*= 0.1 時(shí), 液滴的動(dòng)能向持續(xù)表面能和黏性耗損轉(zhuǎn)化, 液膜邊緣沿徑向鋪展, 鋪展面積增大, 熱流分布區(qū)域也隨之增大, 此時(shí)液膜邊緣處較薄, 邊緣處換熱得到加強(qiáng), 從而在兩側(cè)邊緣處出現(xiàn)兩個(gè)熱流峰, 但此時(shí)液相主體仍集中在壁面中心區(qū), 中心處最高熱流峰仍維持在 5 × 106W/m2, 而同期低速液滴的最高熱流峰僅為其1/10; 當(dāng)t*=0.3時(shí), 液滴底部已形成連續(xù)氣膜, 此時(shí)親水壁面的熱流密度峰值比初期下降了近4 × 106W/m2,但仍高于此時(shí)低速液滴的最高熱流峰, 而且三相接觸線處熱流峰值(1250 W/m2)明顯高于液膜中心(75 W/m2), 熱流由中心向邊緣轉(zhuǎn)移; 原因主要有兩點(diǎn), 一是在三相接觸線附近壁面溫度與液膜之間的溫差大于液膜中心的溫差, 二是在壁面中心區(qū)域氣膜較厚熱阻較大.

    觀察圖6和圖7的液滴形貌演變和熱流密度分布可以看出, 提高撞擊速度后液滴在不同潤濕壁面上的運(yùn)動(dòng)行為和相變行為差異性逐漸減弱, 壁面潤濕作用呈現(xiàn)出較強(qiáng)的速度效應(yīng).Kim等[25]實(shí)驗(yàn)發(fā)現(xiàn), 隨著壁面靜態(tài)接觸角由 63°降低至 2°, 液滴的蒸發(fā)速率由 0.0208 mg·s—1增大至 0.0483 mg·s—1,可見增強(qiáng)壁面的潤濕性有利于液滴在壁面的蒸發(fā).李大樹等[24,26]同樣發(fā)現(xiàn), 壁面潤濕作用隨初始速度的增加而減弱, 呈現(xiàn)出較強(qiáng)的速度效應(yīng), 并將其歸因于液滴的初始速度增大導(dǎo)致壁面接觸角所產(chǎn)生的毛細(xì)力作用降低所致, 并且還發(fā)現(xiàn)增大液滴初始速度, 熱流密度隨之增大, 這與本文的模擬結(jié)果是一致的.

    圖9給出了兩種初速度的液滴撞擊親水壁面后在 0.05 ms時(shí)的壓力分布圖; 速度為 1.6 m/s的液滴撞壁后在壁面中心A(26000 Pa)和邊緣區(qū)B(2000 Pa)的壓力遠(yuǎn)大于速度為0.1 m/s的液滴在中心A(14000 Pa)和邊緣處B(1000 Pa)的壓力,可見高速液滴沿徑向的壓力梯度較大可以快速沿徑向鋪展, 導(dǎo)致?lián)Q熱面積相對增大, 熱流密度因而隨著撞擊速度的增大而增大.

    總體看來, 壁面潤濕性越強(qiáng)換熱效果越好, 撞擊初速度增大帶來的速度效應(yīng)不僅提高了壁面的熱流密度, 也縮小了不同潤濕壁面上的相變行為的差異性; 另外, 兩種速度的液滴撞擊壁后的熱流密度隨時(shí)間推移均劇烈下降, 分析認(rèn)為Leidenfrost效應(yīng)是誘發(fā)上述現(xiàn)象的主要原因.Leidenfrost效應(yīng)涉及的氣膜和液膜都是動(dòng)態(tài)的, 氣膜通過液體的蒸發(fā)供給, 并且由于液體施加在其上方的壓力而流動(dòng), 其特征是在液滴壓在蒸氣層之上, 形成液滴懸浮現(xiàn)象, 并且由于蒸氣反向壓力和液滴慣性力二者的拮抗作用固定了蒸氣層的厚度, 可見氣膜的厚度依賴于反向壓力和慣性力的相對大小.初始速度的改變直接影響慣性力的大小, 圖10給出了兩種初始速度下無量綱氣膜厚度δt/D0隨無量綱時(shí)間t*的變化曲線, 無量綱氣膜厚度定義為壁面中心氣膜厚度δt與液滴初始直徑D0的比值.

    圖9 不同速度的液滴撞壁后在 0.05 ms 時(shí)刻的壓力分布Fig.9.Pressure distribution at 0.05 ms after droplets impact wall with different velocities.

    圖10 初始速度為 (a) 0.1 m/s和 (b) 1.6 m/s 液滴撞擊不同潤濕壁面的無量綱氣膜厚度隨無量綱時(shí)間的變化Fig.10.Dimensional film thickness with dimensionless time curve of droplets impinging on different wetted walls with initial velocity:(a) 0.1 m/s; (b) 1.6 m/s.

    結(jié)合圖6(a)和圖10(a)可以看出, 在氣膜生長初期, 液滴慣性力相對較大, 氣膜厚度較小, 隨著鋪展過程的進(jìn)行, 液滴垂直向下的慣性力逐漸削弱, 同時(shí)部分液相逐漸向氣相轉(zhuǎn)化, 進(jìn)入底部空穴,氣化反向壓力對慣性力的比值快速增大, 反映在圖10(a)中為t*< 0.06階段氣膜厚度快速增長;當(dāng)t*= 0.06時(shí), 接觸角對氣膜厚度的影響逐漸顯現(xiàn), 隨著接觸角的減小氣膜厚度逐漸增大, 親水壁面上穩(wěn)定的氣膜厚度約為疏水壁面的1.3倍, 這與增強(qiáng)潤濕性促進(jìn)換熱效果是符合的, 觀察此時(shí)液滴的形貌可以發(fā)現(xiàn)親水壁面上液膜鋪展的面積明顯大于普通壁面和疏水壁面, 使得親水壁面上液相蒸發(fā)速率遠(yuǎn)高于普通壁面和疏水壁面; 當(dāng) 0.1 <t*<0.16時(shí), 連續(xù)氣膜已形成, 這意味著氣化反向壓力起主導(dǎo)作用, 液膜被托起, 在疏水和普通壁面上空穴內(nèi)部分氣相沿橫向逃逸, 氣膜厚度增長趨勢減緩;當(dāng)t*> 0.16時(shí), 液滴的動(dòng)能完全轉(zhuǎn)化為液體的表面能和黏性耗損, 此時(shí)表面能達(dá)到最大并驅(qū)使液體向液膜中心聚集, 換熱面積隨之減小, 因此在普通和疏水壁面中心氣膜厚度出現(xiàn)輕微下降; 但是此時(shí)親水壁面上氣膜厚度仍持續(xù)增加, 原因是親水壁面上液相向氣相遷移的質(zhì)量較大, 并且在液膜邊緣處生成破碎的小液滴, 液相主體慣性力下降顯著, 有利于氣化反向壓力將液相主體托起.

    分析圖6(b)和圖10(b)發(fā)現(xiàn), 當(dāng)初始速度增大至 1.6 m/s時(shí), 液滴動(dòng)能大幅增加, 液滴在三種壁面上快速鋪展, 固液接觸面積相對增大, 氣液相變劇烈, 因此反映在圖10(b)中為氣膜厚度短期內(nèi)(t*< 0.75)快速增大; 當(dāng)t*> 0.75 時(shí), 氣膜厚度增長趨勢變得緩慢, 隨時(shí)間推移至t*= 4.5時(shí)附近三者氣膜厚度相當(dāng)(無量綱氣膜厚度約為0.23), 并在隨后的過程中氣膜厚度出現(xiàn)輕微波動(dòng), 這可能是因?yàn)樵谧矒艉笃? 換熱被嚴(yán)重惡化, 液滴在壁面破碎成鏈狀的液滴群, 導(dǎo)致氣膜內(nèi)部各處壓力不均勻.另外, 對比圖10(a)和圖10(b)發(fā)現(xiàn)高速液滴撞壁后氣膜厚度整體上顯著小于低速度的液滴, 原因是撞擊速度增大后液體的慣性力遠(yuǎn)大于蒸氣反向壓力[27].Chandra和Aziz[23]開展了液氮滴撞擊壁面的實(shí)驗(yàn)研究, 基于動(dòng)量守恒和能量守恒建立了氣膜厚度與壁面/液體之間的溫差 ?T=Tw-Tl, 鋪展直徑Dt以及撞擊初速度U0的數(shù)學(xué)模型:

    式中,λv為氣相導(dǎo)熱系數(shù),ηv為氣相黏度,L為液相汽化潛熱.

    由(24)式可以看出氣膜厚度與液滴初速度成反比, 這很好地驗(yàn)證了本文中液滴速度增大氣膜厚度減小的結(jié)果.

    液滴初始速度增大有利于降低氣膜厚度, 并且促使液滴頂部溫度較低的液滴能夠沿徑向快速擴(kuò)展, 最終導(dǎo)致熱流密度大幅增加, 這對換熱行為是非常有利的.對于低溫噴霧冷卻而言, 提高熱沉面的潤濕性和液滴撞擊速度對噴霧冷卻有促進(jìn)作用;在膜沸騰區(qū)換熱過程主要集中在撞擊初期, 隨著時(shí)間推移氣膜的厚度逐漸增大, 換熱性能劇烈惡化,因此保持新鮮液滴以一定撞擊速度持續(xù)沖擊熱沉面是影響換熱效果的關(guān)鍵[28].

    3.2 壁面溫度的影響

    壁面溫度是影響液滴撞壁相變行為的關(guān)鍵參數(shù), 尤其是伴隨Leidenfrost效應(yīng)的低溫液滴.研究發(fā)現(xiàn), 液滴撞擊不同溫度壁面的運(yùn)動(dòng)特性差異較大, 可能出現(xiàn)接觸沸騰、接觸反彈、反彈時(shí)液滴部分破碎、液滴下部蒸氣吹掃導(dǎo)致破裂、從液滴上表面噴射、液滴完全分解以及完全分解時(shí)伴隨的破碎片段的快速移動(dòng)等現(xiàn)象[29].圖11和圖12分別給出了直徑為0.5 mm的液滴撞擊溫度分別為300,400, 500 K的壁面后的熱流密度分布和局部溫度分布圖, 靜態(tài)接觸角為 30°.

    圖11 (a) 0.05 ms和 (b) 1.6 ms時(shí)不同溫度壁面上的熱流密度分布Fig.11.Heat flux distribution on different temperature walls at (a) 0.05 ms and (b) 1.6 ms.

    圖12 (a) 0.05 ms和 (b) 1.6 ms不同溫度壁面上的液滴溫度分布Fig.12.Droplets temperature distribution on different temperature wall at (a) 0.05 ms and (b) 1.6 ms.

    綜合分析圖11(a)和圖12(a)的熱流分布、形貌演變以及壁面溫度分布可以看出, 在0.05 ms時(shí)液滴底部未形成連續(xù)氣膜, 液滴底部局部液相仍與壁面接觸(2.7和3.3 mm), 換熱以熱傳導(dǎo)為主從而表現(xiàn)出較高的熱流峰, 而在 2.85, 3.00 和 3.15 mm處出現(xiàn)空穴, 溫度梯度沿徑向分布不均勻, 導(dǎo)致壁面熱流沿徑向由內(nèi)向外呈現(xiàn)出現(xiàn)多個(gè)熱流峰, 而且提高壁面溫度顯然增大了壁面和液體之間的溫差,熱流密度隨之顯著提高; 由圖11(b)可以看出, 在1.6 ms時(shí), 液滴在三種溫度的壁面上均破碎完全形成鏈狀小液滴群, 大量液相蒸發(fā)成氣相, 并且形成較厚的氣膜, 此時(shí)換熱性能急劇下降, 而且不同溫度壁面上的熱流分布曲線差異性也逐漸縮小; 另外觀察圖12(b)發(fā)現(xiàn), 在破碎小液滴區(qū)仍維持一定的溫度梯度, 而在斷裂區(qū)的溫度梯度則顯著下降, 導(dǎo)致壁面出現(xiàn)多個(gè)熱流峰.吳蘇晨等[30]發(fā)現(xiàn)提高壁面和液滴的溫差, 在撞擊開始階段平均熱流密度出現(xiàn)一個(gè)較大的峰值, 此后熱流密度變化緩慢, 分析認(rèn)為在隨后的階段氣膜托起液滴, 液滴無法直接潤濕壁面, 導(dǎo)致熱流密度下降, 換熱效果變差, 但未觀察到液滴破碎, 并且未詳細(xì)分析氣膜厚度隨時(shí)間的變化規(guī)律.

    圖13給出了1.6 m/s液滴撞擊溫度分別為300, 400和500 K的潤濕壁面上氣膜厚度隨無量綱時(shí)間的變化, 可以看出液滴撞擊不同溫度壁面后氣膜厚度隨著時(shí)間推移呈現(xiàn)相同的變化趨勢, 在撞擊初始階段(t* < 0.75)是氣膜厚度快速增長區(qū),而后氣膜增長趨于平緩; 觀察t* = 0.75時(shí)液滴形態(tài)發(fā)現(xiàn), 在500 K壁面上, 液體和壁面的溫差最大,空穴內(nèi)部壓力隨著溫度的升高和液相蒸發(fā)供給而快速增大, 致使此時(shí)底部空穴內(nèi)部壓力大于液滴慣性力, 足以將液滴托舉, 而此時(shí) 300 和 400 K 壁面上連續(xù)氣膜未完全形成, 因而在500 K的壁面上氣膜較厚; 當(dāng)t* > 1.1 時(shí), 觀察 500 K 壁面上液滴的形貌可以看出, 由于表面張力的作用液膜邊緣處產(chǎn)生了破碎小液滴, 此時(shí)氣膜內(nèi)部沿橫向逃逸的蒸氣對小液滴的阻滯作用和液膜內(nèi)部流體的徑向運(yùn)動(dòng), 致使其斜向上運(yùn)動(dòng), 最終液膜邊緣不斷向外剝落小液滴, 在壁面形成鏈狀小液滴群, 換熱性能急劇惡化, 可見Leidenfrost效應(yīng)的存在導(dǎo)致液滴與壁面之間的換熱行為主要發(fā)生在撞擊初期.

    圖13 不同壁面溫度上無量綱氣膜厚度隨無量綱時(shí)間的變化Fig.13.Curves of dimensionless film thickness with dimensionless time on different wall temperatures.

    3.3 氣膜生長行為

    通過上述液滴形貌演變、壁面熱流密度分布和氣膜厚度隨時(shí)間變化規(guī)律的分析可以看出, 氣膜增長是影響低溫液滴換熱的關(guān)鍵, 主要受液體與壁面?zhèn)鳠峥刂? 類似于過熱液體中氣泡的生長, 而換熱性能與氣膜厚度成反比, 且傳熱過程主要集中在撞擊初期.因此下文著重分析撞擊初期氣膜生長行為.為便于理論分析, 假設(shè)液滴撞壁后底部為扁平狀的球體(圖3), 液滴和氣膜之間存在連續(xù)氣膜,其厚度δt遠(yuǎn)小于液滴初始直徑D0和接觸直徑Dt,接觸直徑Dt與液滴初始直徑D0存在如下關(guān)系[31]:

    式中κ-1為毛細(xì)特征長度.

    ΔT在300—500 K時(shí), 輻射換熱約占對流換熱的5%, 因此不考慮輻射換熱的影響[32], 壁面主要通過氣膜的一維導(dǎo)熱向液滴傳遞熱量, 導(dǎo)熱熱通量與壁面和液體之間的溫差 ΔT、鋪展面積 π (Dt/2)2、蒸氣導(dǎo)熱系數(shù)λv成正比, 與氣膜厚度δt成反比,引入蒸發(fā)潛熱L可以推導(dǎo)出單位時(shí)間內(nèi)液體蒸發(fā)質(zhì)量 dm/dt, 依據(jù)能量守恒:

    在撞擊初期氣膜內(nèi)部的壓力主要由撞擊慣性力決定, 且內(nèi)部壓力僅在液滴發(fā)生形變的撞擊初期才有意義.因此在撞擊初期, 液膜壓迫下方氣膜沿徑向逃逸, 逃逸的氣膜通過液體蒸發(fā)產(chǎn)生的蒸氣供給, 此時(shí)壁面和液相界面可近似看成無滑移條件,單位時(shí)間內(nèi)氣膜質(zhì)量增長率與液滴底部蒸發(fā)速率應(yīng)該保持一致[32]:

    聯(lián)立方程(27)和(28)可以得出

    觀察(29)式可以看出在撞擊初期氣膜厚度與溫差和時(shí)間的1/2次方成正比.Breitenbach等[33]基于能量守恒分析在撞擊前期的氣膜增長行為, 得出液滴撞擊高溫壁面后的氣膜厚度與液滴物性、溫差以及撞擊時(shí)間的關(guān)系:

    式中,K為無量綱系數(shù),ew為壁面蓄熱系數(shù),el為液相蓄熱系數(shù).Breitenbach等[33]還發(fā)現(xiàn)其模型的預(yù)測結(jié)果與Tran等[7]的實(shí)驗(yàn)結(jié)果非常吻合.

    圖14 不同壁面溫度上氣膜厚度隨時(shí)間的變化 (a)擬合結(jié)果; (b) Breitenbach 等[33]的分析結(jié)果Fig.14.Curves of film thickness with time on different wall temperatures:(a) Fitting results; (b) Breitenbach et al.[33]analysis results.

    由(29)和(30)式可以看出Breitenbach等[33]的理論模型與本文分析結(jié)果在表達(dá)形式上保持較高的一致性; 圖14(a)擬合了在撞擊初期300—500 K 壁面上氣膜厚度隨時(shí)間的變化 (0.01 <t* <1.05), 這與圖14(b) Breitenbach 等[33]的變化曲線非常相似, 可見在傳熱控制階段氣膜的厚度不僅與壁面/液體之間溫差正相關(guān), 還與時(shí)間的1/2次方呈正相關(guān); 而在撞擊后期, 主要換熱過程基本結(jié)束,換熱性能惡化嚴(yán)重, 氣膜內(nèi)部壓力的作用減小, 氣膜的內(nèi)部應(yīng)力起主要作用, 并在最后階段控制氣膜的微小波動(dòng), 從而表現(xiàn)出圖10和圖13中撞擊后期階段的現(xiàn)象, 氣膜厚度隨時(shí)間變化曲線較為平緩.

    4 結(jié) 論

    采用復(fù)合Level Set-VOF相界面追蹤方法建立單液氮滴撞擊不同潤濕壁面的數(shù)值模型, 并結(jié)合實(shí)驗(yàn)數(shù)據(jù)驗(yàn)證了模型的可靠性, 為低溫液氮噴霧冷卻提供了數(shù)值研究基礎(chǔ), 得出如下結(jié)論:

    1) 增強(qiáng)壁面潤濕性 (30°—150°)、提高液滴的撞擊速度 (0.1和 1.6 m/s)以及提高壁面溫度(300—500 K)均有利于換熱性能的提升, 且不同潤濕壁面上相變行為的差異性因撞擊速度的增大而逐漸減弱, 表現(xiàn)出明顯的速度效應(yīng), 換熱過程主要集中在撞擊初期, 隨著時(shí)間推移, 換熱性能急劇惡化;

    2) 基于質(zhì)量守恒和能量守恒建立了氣膜生長的數(shù)學(xué)模型, 氣膜厚度與壁面溫度和撞擊時(shí)間成正比, 模型預(yù)測結(jié)果很好地驗(yàn)證了數(shù)值分析結(jié)果, 并且與他人實(shí)驗(yàn)和理論分析結(jié)果非常吻合;

    3) 保持新鮮液滴以一定的速度持續(xù)沖擊熱沉面是影響換熱效果的關(guān)鍵, 在膜沸騰階段, 通過強(qiáng)化壁面條件和提高撞擊初速度來促進(jìn)換熱在撞擊初始階段仍可獲得較為顯著的效果, 這對于指導(dǎo)低溫噴霧冷卻有重要的參考意義;

    實(shí)際低溫液氮噴霧冷卻的過程涉及多種復(fù)雜物理場, 耦合具有三維性, 而本文通過建立二維模型研究球形單液滴撞壁后的相變行為, 存在一定的誤差與局限.今后的相關(guān)研究工作將結(jié)合相關(guān)實(shí)驗(yàn)建立精度更高的三維模型開展相關(guān)撞擊過程的數(shù)值研究.

    猜你喜歡
    氣膜潤濕液膜
    考慮軸彎曲的水潤滑軸承液膜建模方法
    T 型槽柱面氣膜密封穩(wěn)態(tài)性能數(shù)值計(jì)算研究
    高空高速氣流下平板液膜流動(dòng)與破裂規(guī)律
    液膜破裂對PCCS降膜的影響*
    基于低場核磁共振表征的礦物孔隙潤濕規(guī)律
    氣膜孔堵塞對葉片吸力面氣膜冷卻的影響
    靜葉柵上游端壁雙射流氣膜冷卻特性實(shí)驗(yàn)
    躲避霧霾天氣的氣膜館
    乙醇潤濕對2種全酸蝕粘接劑粘接性能的影響
    預(yù)潤濕對管道潤濕性的影響
    99热这里只有精品一区 | 精品久久久久久,| 宅男免费午夜| 一本久久中文字幕| 成人三级做爰电影| 欧美日韩瑟瑟在线播放| 在线观看一区二区三区| 国产成人精品久久二区二区91| 欧美日韩福利视频一区二区| 观看免费一级毛片| 女同久久另类99精品国产91| 国产一区在线观看成人免费| 欧美激情在线99| 国产视频内射| 免费看光身美女| 国产av在哪里看| 亚洲av免费在线观看| 嫩草影院入口| 久久午夜综合久久蜜桃| 亚洲欧美精品综合一区二区三区| 香蕉av资源在线| 亚洲精品中文字幕一二三四区| xxxwww97欧美| 日本熟妇午夜| 日本 av在线| 日日干狠狠操夜夜爽| 国产亚洲欧美98| 99久久99久久久精品蜜桃| 精品电影一区二区在线| 长腿黑丝高跟| 亚洲人成电影免费在线| 国产欧美日韩精品一区二区| 神马国产精品三级电影在线观看| 国产精品自产拍在线观看55亚洲| 久久草成人影院| 亚洲无线观看免费| 18禁黄网站禁片免费观看直播| 久久中文字幕一级| 狂野欧美激情性xxxx| 中文资源天堂在线| 久久久国产欧美日韩av| 国产一区二区三区在线臀色熟女| 亚洲人成伊人成综合网2020| 久久久久久久久中文| 12—13女人毛片做爰片一| 2021天堂中文幕一二区在线观| АⅤ资源中文在线天堂| 综合色av麻豆| 好男人在线观看高清免费视频| 精品熟女少妇八av免费久了| 九色国产91popny在线| 国产精品1区2区在线观看.| 精华霜和精华液先用哪个| 51午夜福利影视在线观看| 欧美色欧美亚洲另类二区| 成人精品一区二区免费| 999久久久精品免费观看国产| 精品国产美女av久久久久小说| 欧美黑人巨大hd| 色综合欧美亚洲国产小说| 天天一区二区日本电影三级| 国产免费av片在线观看野外av| 欧美zozozo另类| 欧美一级毛片孕妇| 久久久久国产精品人妻aⅴ院| 在线视频色国产色| 亚洲aⅴ乱码一区二区在线播放| 久久草成人影院| 不卡av一区二区三区| 久久人人精品亚洲av| 久久久精品欧美日韩精品| 精品熟女少妇八av免费久了| 国产亚洲精品av在线| 午夜福利在线在线| 亚洲无线在线观看| 最新在线观看一区二区三区| 亚洲一区二区三区不卡视频| 午夜福利18| 欧美性猛交╳xxx乱大交人| 非洲黑人性xxxx精品又粗又长| 国产激情偷乱视频一区二区| 日本免费一区二区三区高清不卡| 九九热线精品视视频播放| 亚洲国产色片| 亚洲中文日韩欧美视频| 村上凉子中文字幕在线| 日日干狠狠操夜夜爽| 老司机福利观看| 欧美日韩瑟瑟在线播放| 日韩欧美一区二区三区在线观看| a在线观看视频网站| 欧洲精品卡2卡3卡4卡5卡区| 国产精品1区2区在线观看.| 亚洲片人在线观看| 美女cb高潮喷水在线观看 | 亚洲熟妇中文字幕五十中出| 99久久国产精品久久久| 欧美日韩国产亚洲二区| 精品人妻1区二区| 成人特级av手机在线观看| 狂野欧美白嫩少妇大欣赏| 波多野结衣巨乳人妻| 91字幕亚洲| 国产一区二区激情短视频| 国产真实乱freesex| 免费在线观看日本一区| 亚洲av成人不卡在线观看播放网| 久久久久久大精品| 免费看日本二区| 无限看片的www在线观看| 亚洲成人精品中文字幕电影| а√天堂www在线а√下载| 欧美色欧美亚洲另类二区| 黑人操中国人逼视频| 此物有八面人人有两片| 搡老妇女老女人老熟妇| 日本三级黄在线观看| 一二三四在线观看免费中文在| 婷婷六月久久综合丁香| 精品欧美国产一区二区三| 人妻丰满熟妇av一区二区三区| 麻豆久久精品国产亚洲av| 亚洲熟妇中文字幕五十中出| 一进一出抽搐动态| 国产精品99久久久久久久久| 国产蜜桃级精品一区二区三区| 免费看十八禁软件| 国产av麻豆久久久久久久| 好男人电影高清在线观看| 亚洲av电影在线进入| 老司机午夜十八禁免费视频| 亚洲无线观看免费| 久久婷婷人人爽人人干人人爱| 老鸭窝网址在线观看| 天堂影院成人在线观看| ponron亚洲| 国产欧美日韩精品亚洲av| 麻豆国产av国片精品| 99久久精品一区二区三区| 神马国产精品三级电影在线观看| 色播亚洲综合网| 国产av一区在线观看免费| 亚洲中文字幕一区二区三区有码在线看 | 宅男免费午夜| 国产成人av教育| 国产亚洲精品一区二区www| 国产三级在线视频| 91久久精品国产一区二区成人 | 国产 一区 欧美 日韩| 麻豆久久精品国产亚洲av| 亚洲av电影在线进入| 狂野欧美激情性xxxx| 国产乱人伦免费视频| 国产精品久久久久久精品电影| 成人无遮挡网站| 精品国产乱子伦一区二区三区| 国产精品综合久久久久久久免费| 淫秽高清视频在线观看| 亚洲欧美激情综合另类| 国产精品亚洲一级av第二区| 亚洲精品美女久久久久99蜜臀| 99在线视频只有这里精品首页| 一级毛片精品| 嫩草影视91久久| 精品久久久久久久毛片微露脸| 亚洲人成伊人成综合网2020| 一个人看的www免费观看视频| 99精品欧美一区二区三区四区| 欧美日韩一级在线毛片| 国产精品99久久久久久久久| 精品国产三级普通话版| 一本一本综合久久| 亚洲自偷自拍图片 自拍| 久久精品国产亚洲av香蕉五月| 久久九九热精品免费| 黑人欧美特级aaaaaa片| 神马国产精品三级电影在线观看| 国产精品 国内视频| 99热这里只有是精品50| 午夜免费成人在线视频| 精品久久久久久久人妻蜜臀av| 三级国产精品欧美在线观看 | 真人做人爱边吃奶动态| 午夜a级毛片| 俺也久久电影网| 99热这里只有是精品50| 久久久国产欧美日韩av| www.www免费av| 人妻丰满熟妇av一区二区三区| 此物有八面人人有两片| 99久久无色码亚洲精品果冻| 99国产精品99久久久久| 美女午夜性视频免费| 麻豆av在线久日| 美女被艹到高潮喷水动态| 日本免费一区二区三区高清不卡| 国产成人精品无人区| av黄色大香蕉| 嫩草影视91久久| 90打野战视频偷拍视频| 日本黄色片子视频| 亚洲熟女毛片儿| 19禁男女啪啪无遮挡网站| 夜夜夜夜夜久久久久| 757午夜福利合集在线观看| 国产欧美日韩一区二区三| 别揉我奶头~嗯~啊~动态视频| 午夜日韩欧美国产| 在线视频色国产色| 日本一二三区视频观看| 一个人看视频在线观看www免费 | 美女高潮喷水抽搐中文字幕| 麻豆国产av国片精品| 中文资源天堂在线| 精品国产美女av久久久久小说| 精品一区二区三区视频在线 | 国产精品综合久久久久久久免费| 亚洲精华国产精华精| 日日夜夜操网爽| 欧美av亚洲av综合av国产av| 午夜福利在线观看吧| 国产av麻豆久久久久久久| 啦啦啦免费观看视频1| 黄片大片在线免费观看| 欧美日本视频| 日本 av在线| 99精品久久久久人妻精品| 亚洲一区二区三区不卡视频| 波多野结衣高清作品| 日韩有码中文字幕| 亚洲专区字幕在线| 宅男免费午夜| 国产日本99.免费观看| 欧美精品啪啪一区二区三区| 亚洲18禁久久av| 精品熟女少妇八av免费久了| 黑人操中国人逼视频| 日韩欧美精品v在线| 日本 欧美在线| 久久热在线av| 变态另类成人亚洲欧美熟女| 亚洲成人精品中文字幕电影| 高清在线国产一区| 小蜜桃在线观看免费完整版高清| 欧美日韩精品网址| 国产三级中文精品| 噜噜噜噜噜久久久久久91| 亚洲成人久久爱视频| 精品一区二区三区av网在线观看| www日本在线高清视频| 日本五十路高清| 亚洲片人在线观看| 国产精品综合久久久久久久免费| 18禁国产床啪视频网站| 久久精品91无色码中文字幕| 欧美色欧美亚洲另类二区| 91在线观看av| 日本与韩国留学比较| 一个人免费在线观看的高清视频| 国产精品自产拍在线观看55亚洲| 在线观看日韩欧美| 国产精品一区二区免费欧美| 国产精品亚洲一级av第二区| 久久久国产精品麻豆| 级片在线观看| 少妇丰满av| bbb黄色大片| 欧美日韩综合久久久久久 | 国产黄a三级三级三级人| 99久久综合精品五月天人人| 97碰自拍视频| 欧美色视频一区免费| 亚洲最大成人中文| 欧美一区二区国产精品久久精品| 免费在线观看日本一区| 久久精品国产清高在天天线| 国产69精品久久久久777片 | 欧美3d第一页| 日韩人妻高清精品专区| 久久午夜亚洲精品久久| 国产精品一区二区三区四区免费观看 | 亚洲乱码一区二区免费版| 久久天躁狠狠躁夜夜2o2o| 99久久精品一区二区三区| 男女做爰动态图高潮gif福利片| 一进一出抽搐动态| 亚洲电影在线观看av| 欧美成人性av电影在线观看| 欧美最黄视频在线播放免费| 热99在线观看视频| 又爽又黄无遮挡网站| 国产一区二区三区视频了| 51午夜福利影视在线观看| 怎么达到女性高潮| 18禁国产床啪视频网站| 亚洲国产中文字幕在线视频| 精品久久久久久,| 久9热在线精品视频| 亚洲av五月六月丁香网| 国产一区二区三区在线臀色熟女| 脱女人内裤的视频| 亚洲av美国av| 日本在线视频免费播放| 日本五十路高清| 两性夫妻黄色片| 成年版毛片免费区| 在线观看美女被高潮喷水网站 | 桃红色精品国产亚洲av| av福利片在线观看| cao死你这个sao货| 欧美绝顶高潮抽搐喷水| 国产一区二区三区在线臀色熟女| 麻豆一二三区av精品| 在线免费观看的www视频| 欧美中文综合在线视频| 99re在线观看精品视频| 欧美丝袜亚洲另类 | 成人18禁在线播放| 2021天堂中文幕一二区在线观| 99在线视频只有这里精品首页| 小说图片视频综合网站| 最近最新免费中文字幕在线| 色老头精品视频在线观看| 九色成人免费人妻av| 伦理电影免费视频| 在线观看一区二区三区| 午夜福利成人在线免费观看| 免费看a级黄色片| av女优亚洲男人天堂 | 熟女人妻精品中文字幕| 999久久久精品免费观看国产| 夜夜爽天天搞| 香蕉久久夜色| 亚洲无线在线观看| 天堂av国产一区二区熟女人妻| 美女高潮喷水抽搐中文字幕| 国产精品永久免费网站| 波多野结衣高清无吗| 久久性视频一级片| 亚洲成人免费电影在线观看| 久久精品91无色码中文字幕| 成年版毛片免费区| 亚洲aⅴ乱码一区二区在线播放| 亚洲成av人片在线播放无| 亚洲国产欧洲综合997久久,| 国产精品 欧美亚洲| 美女大奶头视频| 国内精品久久久久久久电影| 国产三级中文精品| 狂野欧美白嫩少妇大欣赏| 两个人的视频大全免费| 久久久久九九精品影院| av福利片在线观看| 欧美激情在线99| 亚洲专区字幕在线| 叶爱在线成人免费视频播放| 国产激情偷乱视频一区二区| 天堂√8在线中文| 级片在线观看| 国产精品av视频在线免费观看| 亚洲精品456在线播放app | 日本黄色片子视频| 人人妻,人人澡人人爽秒播| 国产真实乱freesex| 999久久久精品免费观看国产| 99riav亚洲国产免费| 亚洲色图av天堂| 搡老妇女老女人老熟妇| 欧美绝顶高潮抽搐喷水| 岛国在线免费视频观看| 免费av不卡在线播放| 欧美黑人巨大hd| 国产乱人伦免费视频| 亚洲国产精品成人综合色| 亚洲欧美日韩东京热| 国产精品女同一区二区软件 | 99国产综合亚洲精品| 一级毛片高清免费大全| 亚洲av日韩精品久久久久久密| 久久草成人影院| 欧美日韩综合久久久久久 | 天堂√8在线中文| 成年女人毛片免费观看观看9| www.自偷自拍.com| 美女 人体艺术 gogo| 12—13女人毛片做爰片一| 欧美乱妇无乱码| 亚洲电影在线观看av| 免费观看的影片在线观看| 在线观看免费视频日本深夜| 成人无遮挡网站| 国产精品一区二区精品视频观看| 久久精品91无色码中文字幕| 国产亚洲av高清不卡| 麻豆av在线久日| h日本视频在线播放| 69av精品久久久久久| 天天一区二区日本电影三级| 亚洲av熟女| 91九色精品人成在线观看| 亚洲av成人av| 老司机福利观看| 免费在线观看影片大全网站| 国产精品久久久久久人妻精品电影| 欧美最黄视频在线播放免费| 宅男免费午夜| 91老司机精品| АⅤ资源中文在线天堂| 人人妻人人看人人澡| 国产亚洲欧美在线一区二区| netflix在线观看网站| av在线天堂中文字幕| 国产成+人综合+亚洲专区| 我要搜黄色片| 老熟妇乱子伦视频在线观看| 又大又爽又粗| 国产99白浆流出| 别揉我奶头~嗯~啊~动态视频| 日韩免费av在线播放| 国产成人av教育| 午夜福利在线在线| 成年女人毛片免费观看观看9| 巨乳人妻的诱惑在线观看| 亚洲精品456在线播放app | 亚洲精品456在线播放app | 熟女少妇亚洲综合色aaa.| www日本黄色视频网| 国产麻豆成人av免费视频| 一个人看的www免费观看视频| 久久久国产精品麻豆| 亚洲avbb在线观看| 无遮挡黄片免费观看| 欧美丝袜亚洲另类 | 人妻久久中文字幕网| 午夜福利18| 久久精品国产99精品国产亚洲性色| 亚洲国产精品999在线| 久久精品国产99精品国产亚洲性色| 国产一区二区在线观看日韩 | 亚洲男人的天堂狠狠| 男女做爰动态图高潮gif福利片| 国产乱人视频| 亚洲最大成人中文| 国产欧美日韩精品亚洲av| 婷婷丁香在线五月| ponron亚洲| 日本 欧美在线| 韩国av一区二区三区四区| 1024香蕉在线观看| 成年女人毛片免费观看观看9| 一进一出抽搐gif免费好疼| 亚洲精华国产精华精| 国产午夜精品论理片| 国产欧美日韩精品亚洲av| 成熟少妇高潮喷水视频| 日韩精品青青久久久久久| 男人的好看免费观看在线视频| 不卡av一区二区三区| 国产精品爽爽va在线观看网站| 精品久久久久久久久久久久久| 亚洲国产看品久久| 两人在一起打扑克的视频| 不卡一级毛片| 又大又爽又粗| 国产精品久久久久久久电影 | 免费在线观看亚洲国产| 午夜福利欧美成人| 色尼玛亚洲综合影院| 国产精品香港三级国产av潘金莲| 日本黄色片子视频| 99国产极品粉嫩在线观看| 18禁黄网站禁片午夜丰满| 国产精品久久久久久久电影 | 欧美日韩乱码在线| 国产伦精品一区二区三区四那| 黄片小视频在线播放| 99久久成人亚洲精品观看| 啦啦啦免费观看视频1| 欧美色欧美亚洲另类二区| 成人三级黄色视频| 观看免费一级毛片| 久久久久久大精品| 少妇熟女aⅴ在线视频| 啪啪无遮挡十八禁网站| 熟女电影av网| 成人性生交大片免费视频hd| 村上凉子中文字幕在线| 午夜成年电影在线免费观看| 久久久久久久久久黄片| 成人欧美大片| 99国产综合亚洲精品| 久久久久国产精品人妻aⅴ院| 久久精品91无色码中文字幕| 婷婷精品国产亚洲av| 色精品久久人妻99蜜桃| 手机成人av网站| 国产高清有码在线观看视频| 国产亚洲精品一区二区www| 欧美激情在线99| 国产伦一二天堂av在线观看| 久久久久久大精品| 精品一区二区三区av网在线观看| 午夜a级毛片| 亚洲欧洲精品一区二区精品久久久| АⅤ资源中文在线天堂| 日韩欧美在线乱码| 亚洲av第一区精品v没综合| 亚洲色图av天堂| 给我免费播放毛片高清在线观看| 大型黄色视频在线免费观看| 欧美一级毛片孕妇| 丰满人妻一区二区三区视频av | 午夜福利在线观看吧| 成人永久免费在线观看视频| 观看美女的网站| 精品福利观看| 香蕉国产在线看| 人人妻,人人澡人人爽秒播| 成人三级做爰电影| 岛国视频午夜一区免费看| 欧美高清成人免费视频www| 亚洲av熟女| 久久精品91蜜桃| 99久久国产精品久久久| av国产免费在线观看| 天堂影院成人在线观看| 久久久久久久午夜电影| 岛国在线观看网站| 精品一区二区三区视频在线 | 国产三级黄色录像| 成人高潮视频无遮挡免费网站| 高清毛片免费观看视频网站| 男人舔奶头视频| 人妻久久中文字幕网| 国产高清视频在线观看网站| 国产精品99久久久久久久久| 亚洲五月婷婷丁香| 久久久国产欧美日韩av| 最近最新中文字幕大全电影3| 午夜福利免费观看在线| 少妇熟女aⅴ在线视频| 黄片大片在线免费观看| 午夜福利在线在线| 色尼玛亚洲综合影院| 久久久国产精品麻豆| 成人国产综合亚洲| 后天国语完整版免费观看| 欧美日本亚洲视频在线播放| 国内久久婷婷六月综合欲色啪| 亚洲国产欧洲综合997久久,| 欧美日本视频| 午夜福利在线观看吧| 97碰自拍视频| 午夜成年电影在线免费观看| 麻豆国产97在线/欧美| 一个人看视频在线观看www免费 | 国产综合懂色| 欧美性猛交╳xxx乱大交人| 国产视频一区二区在线看| 国产成人aa在线观看| 丁香欧美五月| 天堂网av新在线| 在线观看美女被高潮喷水网站 | 狂野欧美激情性xxxx| 88av欧美| 欧美绝顶高潮抽搐喷水| 久久久色成人| 97人妻精品一区二区三区麻豆| 中文字幕最新亚洲高清| 露出奶头的视频| 亚洲精品一卡2卡三卡4卡5卡| 久久这里只有精品19| 国产伦精品一区二区三区视频9 | 欧美高清成人免费视频www| 亚洲精品久久国产高清桃花| 精品久久久久久,| 精品国产乱子伦一区二区三区| 男女之事视频高清在线观看| 一区二区三区高清视频在线| 亚洲av电影不卡..在线观看| 亚洲无线在线观看| 日本黄色视频三级网站网址| 成年女人毛片免费观看观看9| 床上黄色一级片| 一二三四在线观看免费中文在| 桃色一区二区三区在线观看| 国产高清视频在线观看网站| 白带黄色成豆腐渣| 国产欧美日韩一区二区精品| 免费看a级黄色片| 欧美中文综合在线视频| 亚洲精品美女久久av网站| 真实男女啪啪啪动态图| 国产亚洲欧美在线一区二区| 国产熟女xx| av视频在线观看入口| 黄色成人免费大全| 嫁个100分男人电影在线观看| 国内精品美女久久久久久| 亚洲人成伊人成综合网2020| 日韩成人在线观看一区二区三区| 欧美成人免费av一区二区三区| 国产精品一及| 国产伦精品一区二区三区四那| 91av网站免费观看| 婷婷精品国产亚洲av| 后天国语完整版免费观看| 午夜免费成人在线视频| 亚洲男人的天堂狠狠| 日本精品一区二区三区蜜桃| 99久久99久久久精品蜜桃| 日韩av在线大香蕉| 人妻夜夜爽99麻豆av| e午夜精品久久久久久久| 搡老岳熟女国产|