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

    納米通道粗糙內(nèi)壁對流體流動(dòng)行為的影響*

    2019-05-17 06:42:52梅濤陳占秀楊歷王坤苗瑞燦
    物理學(xué)報(bào) 2019年9期
    關(guān)鍵詞:潤濕性固液勢能

    梅濤 陳占秀 楊歷 王坤 苗瑞燦

    (河北工業(yè)大學(xué)能源與環(huán)境工程學(xué)院,天津 300401)

    納米流動(dòng)系統(tǒng)具有高效、經(jīng)濟(jì)等優(yōu)勢,在眾多領(lǐng)域具有廣泛的應(yīng)用前景.因該類系統(tǒng)具有極高的表面積體積比,致使界面滑移效應(yīng)對流動(dòng)具有顯著影響.本文采用分子動(dòng)力學(xué)方法以兩無限大平行非對稱壁面組成的Poiseuille流動(dòng)為對象,分析了壁面粗糙度與潤濕性變化對通道內(nèi)流體流動(dòng)的影響.對于不同結(jié)構(gòu)類型的壁面,需要通過水動(dòng)力位置來確定固液界面位置,準(zhǔn)確計(jì)算固液界面位置有助于更好地分析界面滑移效應(yīng).研究結(jié)果表明,上下壁面不對稱會(huì)引起通道內(nèi)流場參數(shù)分布的不對稱,壁面粗糙度及潤濕性的變化會(huì)影響近壁面附近流體原子的流動(dòng)特性,由于壁面凹槽的存在,粗糙壁面附近的數(shù)密度分布低于光滑壁面一側(cè).壁面粗糙度及潤濕性的變化會(huì)影響固液界面位置,肋高變化及壁面潤濕性對通道中速度分布影響較大,界面滑移速度及滑移長度隨肋高和潤濕性的增大而減小;肋間距變化對通道內(nèi)流體流動(dòng)影響較小,界面滑移速度和滑移長度基本保持恒定.

    1 引 言

    對納米流體的深入研究是科技微型化導(dǎo)致的必然趨勢,納米流體具有優(yōu)于常規(guī)介質(zhì)的特性,包括運(yùn)輸、傳熱、吸附等性能,對納米流體的深入研究將更有利于科技微型化的發(fā)展.納米量級流體在流動(dòng)時(shí)存在許多微尺度效應(yīng),包括分子力效應(yīng)、低雷諾數(shù)效應(yīng)以及表面張力效應(yīng)等.相比于宏觀尺度,微尺度條件下對流體流動(dòng)規(guī)律的研究受實(shí)驗(yàn)條件和測量精度的限制,使得數(shù)值模擬成為了微尺度流動(dòng)的主要研究方法.而以微觀原子之間的作用模型為基礎(chǔ)的分子動(dòng)力學(xué),因不引入常規(guī)假設(shè),能以原子級精度描述流體在近鄰壁面處的流動(dòng)特性,是目前研究微尺度流動(dòng)問題的主要數(shù)值模擬方法.

    通常宏觀尺度下流體在通道內(nèi)的流動(dòng)特性都是基于無滑移邊界條件進(jìn)行研究,因?yàn)楸诿娈a(chǎn)生的表面效應(yīng)對近壁處流體流動(dòng)特性影響較小.然而,隨著流動(dòng)通道尺寸的減小,比表面積急劇增加,微尺度下固體壁面對近壁處流體的影響直接關(guān)系到通道內(nèi)整體的流動(dòng)特性分布[1-3].Sun和Ebner[4]研究了固液界面間作用強(qiáng)度大小對流體流動(dòng)特性的影響,發(fā)現(xiàn)當(dāng)固液原子間作用力較強(qiáng)時(shí),流體在壁面處幾乎無滑移,而當(dāng)作用力較弱時(shí),流體在壁面處會(huì)產(chǎn)生較大的邊界滑移.Voronov等[5]利用分子動(dòng)力學(xué)對Couette流動(dòng)進(jìn)行模擬,結(jié)果顯示滑移長度與接觸角相關(guān),對于疏水性表面,接觸角會(huì)增大,滑移長度也會(huì)隨之增大,但也會(huì)出現(xiàn)隨著接觸角增大滑移長度減小的情況.胡海豹等[6]利用分子動(dòng)力學(xué)對超疏水納米通道內(nèi)流體流動(dòng)特性進(jìn)行研究,結(jié)果表明當(dāng)接觸角大于150°時(shí),滑移速度和滑移長度出現(xiàn)隨接觸角增大而減小的反常現(xiàn)象,并進(jìn)一步證明了改變固液原子間的勢能參數(shù)表征的潤濕性不能用來準(zhǔn)確表示超疏水壁面對流體的影響.Nagayama和Cheng[7]對納米通道內(nèi)的Poiseuille流動(dòng)進(jìn)行分析,發(fā)現(xiàn)改變壁面與流體間的勢能參數(shù)以及添加在流體部分的驅(qū)動(dòng)力均會(huì)影響固液界面間流體原子的運(yùn)動(dòng)特性,固液間相互作用越強(qiáng),邊界滑移速度越小,而添加的驅(qū)動(dòng)力越大,滑移速度也越大.Barisik和Beskok[8]以及Shi等[9]利用智能壁分子動(dòng)力學(xué)方法模擬研究了納米通道內(nèi)氣體原子的流動(dòng)特性,結(jié)果表明近壁區(qū)域內(nèi)氣體的速度、密度和壓力的變化趨勢僅由壁面力場決定,不受通道高度以及密度的影響.張冉等[10]利用分子動(dòng)力學(xué)分析了納米通道內(nèi)近壁區(qū)域流體的流動(dòng)特性,同樣發(fā)現(xiàn)近壁區(qū)域的氣體流動(dòng)特性僅由壁面力場決定,壁面對氣體原子的勢能作用越強(qiáng),氣體在近壁區(qū)域的密度越大,直至形成吸附層.Voronov等[5]利用分子動(dòng)力學(xué)方法對Couette流動(dòng)中流體的滑移現(xiàn)象進(jìn)行了研究,發(fā)現(xiàn)固液間勢能強(qiáng)度以及平衡距離的減小均傾向于增大接觸角,而減小勢能強(qiáng)度會(huì)增加滑移長度,減小平衡距離會(huì)減少滑移長度,說明接觸角與滑移長度間的關(guān)系并不唯一.Cieplak[11]同樣利用分子動(dòng)力學(xué)方法對Couette流動(dòng)進(jìn)行研究,主要探究了固液間的作用強(qiáng)度以及不同流體介質(zhì)對流體滑移的影響,結(jié)果表明,滑移長度與固液間的作用有直接關(guān)系,而與流體介質(zhì)無關(guān).

    可見,目前科研工作者大多集中于進(jìn)行固液間相互作用等因素對通道內(nèi)流體流動(dòng)影響規(guī)律的研究,且固體表面多是光滑壁面.實(shí)際上,任何固體表面都不可能是絕對光滑的,當(dāng)宏觀尺度轉(zhuǎn)為納微尺度,流動(dòng)通道尺度急劇減少,比表面積也隨之急劇增加,固體表面的粗糙程度對流體流動(dòng)的作用也相應(yīng)凸顯.因此,粗糙壁面對流動(dòng)的影響機(jī)理研究已逐漸成為當(dāng)前納微尺度傳熱傳質(zhì)研究中的重點(diǎn)[12].張程賓等[13]利用分子動(dòng)力學(xué)方法對含粗糙壁面納米通道內(nèi)的流體流動(dòng)進(jìn)行了研究,發(fā)現(xiàn)粗糙壁面會(huì)限制近壁區(qū)流體原子的運(yùn)動(dòng),導(dǎo)致流體流動(dòng)速度及滑移速度降低.Rahmatipour等[14]利用分子動(dòng)力學(xué)對含粗糙壁面的納米通道內(nèi)的肋高變化進(jìn)行了研究,結(jié)果顯示,相比于光滑壁面,隨著肋高的增大,壁面附近流體的密度波動(dòng)范圍逐漸增大,但波動(dòng)的峰值均低于光滑壁面.Toghraie等[15]對含納米顆粒的粗糙通道內(nèi)的流體流動(dòng)特性進(jìn)行了研究,通過設(shè)定肋高與肋間距的比值來獲取不同壁面粗糙度,結(jié)果表明通道內(nèi)流體的溫度及速度分布并不隨著粗糙度的增大而增大,而是由凹槽內(nèi)限制的流體原子數(shù)量決定,限制的流體原子數(shù)量數(shù)量越多,流體流動(dòng)速度越小,反之流動(dòng)速度越大.但由該文獻(xiàn)的物理模型可知,在研究粗糙度變化時(shí),粗糙壁面上肋的尺寸及數(shù)量也發(fā)生了變化,故不能很好地詮釋壁面肋間距及肋高變化對近壁區(qū)流體的影響.

    綜上所述,目前研究壁面的粗糙度變化對流動(dòng)的影響還不夠詳細(xì),而針對粗糙壁面潤濕性變化的研究更是較少.為深入揭示納微尺度下粗糙壁面對流體流動(dòng)的影響機(jī)理,本文擬構(gòu)建含粗糙壁面的納米通道內(nèi)流體流動(dòng)的分子動(dòng)力學(xué)模型,并與光滑壁面進(jìn)行對比,分析粗糙壁面肋高及肋間距對流體流動(dòng)特性的影響.在此基礎(chǔ)上,討論壁面潤濕性對壁面附近流體原子的影響,并揭示粗糙壁面與光滑壁面潤濕性變化的異同.

    2 物理模型及模擬設(shè)置

    納米通道流體流動(dòng)的分子動(dòng)力學(xué)模型如圖1(a)所示.模擬通道尺寸為Lx×Ly×Lz= 18.20σ× 19.68σ× 7.87σ,納米通道上壁為光滑壁面,下壁為粗糙壁面,壁面由銅原子構(gòu)成,流體單原子氬均勻分布于上下壁面之間,且均按照面心立方結(jié)構(gòu)進(jìn)行排列.圖1(b)為納米結(jié)構(gòu)示意圖,壁面長度為Lx= 18.20σ,壁面厚度D= 2.95σ,粗糙壁面肋間距為a,肋高為h,肋長度為b,通道整體高度為H,光滑壁面與切面間為平直流道,高度為H′.計(jì)算統(tǒng)計(jì)時(shí),對通道高度H進(jìn)行分層,其中統(tǒng)計(jì)流場中數(shù)密度分布時(shí)分為70層,對速度分布進(jìn)行統(tǒng)計(jì)時(shí)分為44層.整個(gè)模擬系統(tǒng)在x,z方向上設(shè)置周期性邊界條件,在y方向上設(shè)置固定邊界條件.

    流體原子間以及固液原子間的勢能作用均采用Lennard-Jones (12-6)勢能模型,其表達(dá)式為

    式中rij為原子間的距離,ε為能量參數(shù),σ為尺寸參數(shù),rc為截?cái)喟霃?流體氬原子間的能量參數(shù)ε=1.67×10-21J,尺寸參數(shù)σ= 0.3405 nm.壁面固體原子勢能參數(shù)εs=40ε,尺寸參數(shù)σs=0.69σ,截?cái)喟霃饺c=2.5σ.

    固液原子間的作用力大小決定了壁面潤濕性,流體在壁面鋪展得越充分,固液界面間的接觸角越小,可認(rèn)為潤濕性越好.而對于光滑壁面而言,接觸角的理論公式為

    圖1 (a)模擬系統(tǒng)圖;(b)納米結(jié)構(gòu)示意圖Fig.1.(a) Simulation system;(b) schematic of nanostructure.

    式中θ為光滑壁面接觸角,εsl為固液原子間的能量參數(shù).而對于粗糙壁面,(2)式將不再適用,針對本文所研究的物理模型,粗糙壁面潤濕性可表示為

    式中r表示粗糙壁面的粗糙程度;a,b及h為圖1(b)中所示參數(shù).令,c為固液原子間的勢能系數(shù).本文模擬工況與對應(yīng)粗糙度下的接觸角如表1所列.

    對于壁面原子與流體原子間的尺寸參數(shù)σsl,由Lorentz-Berthelot混合法進(jìn)行計(jì)算:

    平板間的Poiseuille流動(dòng)是通過在x方向上施加外部驅(qū)動(dòng)力來驅(qū)動(dòng)流體流動(dòng)的,并利用拋物線求解Navier-Stokes方程.Poiseuille流動(dòng)被定義為在通道高度H(y方向上由0到H)內(nèi)的流動(dòng),Navier-Stokes方程可簡化為[16]

    式中μ為剪切黏度,ux為x方向上的流動(dòng)速度,ρ為流體密度.對(6)式連續(xù)進(jìn)行兩次積分,將中心對稱條件和邊界滑移條件ux|y=0=us或ux|y=H=us依次代入,其中us為邊界滑移速度,Poiseuille流動(dòng)的速度場可表示為

    固液界面間的滑移規(guī)律可以由滑移長度ls表示,表達(dá)式為[17-20]

    式中 (?ux/?y)|y|=H為固液界面處流體的速度梯度.對于粗糙壁面滑移長度的確定,需要根據(jù)水動(dòng)力位置來研究[21].圖2(a)和圖2(b)所示分別為Couette流動(dòng)和Poiseuille流動(dòng)示意圖,圖2(a)中yC為外推線性的速度直線達(dá)到壁面速度時(shí)對應(yīng)的位置,圖2(b)中yP為外推拋物線的速度曲線達(dá)到壁面速度的對應(yīng)位置.通過文獻(xiàn)[21]可知,粗糙壁面滑移長度可表示為,而水動(dòng)力位置可表示為yH=yC+ls.

    表1 不同模擬工況下對應(yīng)的粗糙度與接觸角Table 1.Corresponding roughness and contact angle under different simulation conditions.

    圖2 模型結(jié)構(gòu)示意圖 (a) Couette流動(dòng);(b) Poiseuille流動(dòng)Fig.2.Schematic of nanostructure:(a) Couette flow;(b) Poiseuille flow.

    本文模擬過程中,改變壁面粗糙度會(huì)影響通道內(nèi)流體密度與壓力的變化,但變化幅度較小,通道中流體原子的密度基本保持在(1.23 ± 0.01) g/cm3范圍內(nèi),造成的計(jì)算誤差可以忽略.粗糙壁面肋的長度b= 2.7σ,肋的數(shù)量為3.采用Velocity-Verlet算法對流體原子運(yùn)動(dòng)方程進(jìn)行求解,時(shí)間步長為1 fs.首先利用正則系綜對初始模型進(jìn)行弛豫,利用Nose-Hoover恒溫算法將系統(tǒng)溫度恒定在T=0.827ε/kB,流體原子的速度遵循高斯分布,進(jìn)行50萬步使系統(tǒng)達(dá)到穩(wěn)定.弛豫后對模型入口處厚度為2.0σ的流體區(qū)域施加沿x正方向的水平驅(qū)動(dòng)力F= 0.05ε/σ,采用正則系綜對系統(tǒng)進(jìn)行溫度控制,共運(yùn)行400萬步,前200萬步使系統(tǒng)達(dá)到穩(wěn)定,后200萬步對流場中所需參數(shù)進(jìn)行計(jì)算統(tǒng)計(jì).分析模型水動(dòng)力位置時(shí)需要計(jì)算Couette流動(dòng)的速度場分布,設(shè)置上下壁面在水平位置上以相同的速度反向運(yùn)動(dòng),速度大小為.模擬采用LAMMPS程序編寫[22].

    3 模擬方法驗(yàn)證

    為了驗(yàn)證模型及參數(shù)設(shè)置的正確性,本文基于兩平行光滑平板間的Poiseuille流動(dòng),對不同壁面潤濕性下流體在通道內(nèi)的密度分布進(jìn)行驗(yàn)證.具體參數(shù)設(shè)置如下:流體采用單原子氬,其中ε=1.67×10-21J,σ= 0.3405 nm;固液原子間的距離參數(shù)σsl=0.2872 nm,無量綱化后,取固液原子間的勢能系數(shù)c= 2.0,1.0,0.6,0.2;上下壁面設(shè)置為固體剛性壁面,截?cái)喟霃絩c= 2.5σ.

    圖3所示為不同固液原子間勢能系數(shù)c對應(yīng)的微通道內(nèi)流體的無量綱密度分布.由圖3可知,模擬數(shù)據(jù)與文獻(xiàn)[23]所給數(shù)據(jù)基本相符,在近壁面區(qū)域內(nèi),由于表面效應(yīng)的存在使得流體密度分布均出現(xiàn)了有序振蕩現(xiàn)象,流體密度不均勻;而在通道中心的主流區(qū),流體受壁面的影響較小而趨近于平緩.因此,可認(rèn)為本文建立理論模型、選用算法及編寫的程序準(zhǔn)確可靠.

    4 結(jié)果與分析

    4.1 壁面粗糙度變化對流動(dòng)特性的影響

    本節(jié)以固液原子間勢能系數(shù)c= 0.75,對應(yīng)的光滑接觸角θ= 60°時(shí)的工況為例.通過改變粗糙壁面上肋高及肋間距研究壁面粗糙度對通道內(nèi)流體數(shù)密度及速度分布的影響.

    4.1.1 數(shù)密度分布

    為研究肋高變化對近壁區(qū)流體數(shù)密度的影響,固定肋間距a= 3.6σ不變,分別取肋高h(yuǎn)=0.5σ,1.0σ,1.5σ,2.0σ,研究微通道內(nèi)流體流動(dòng)數(shù)密度分布曲線規(guī)律.如圖4所示,橫坐標(biāo)為沿y方向的高度,縱坐標(biāo)為無量綱數(shù)密度ρ?=ρσ3.由圖4可知,由于壁面效應(yīng)使得近壁區(qū)流體原子分布不均勻,近壁區(qū)域流體數(shù)密度分布出現(xiàn)明顯的振蕩現(xiàn)象,而通道主流區(qū)域流體受壁面影響較小,數(shù)密度分布基本保持恒定.由于通道壁面形狀的不對稱,導(dǎo)致流場中數(shù)密度分布的不對稱,由圖4(a)和圖4(b)所示,光滑壁面附近流體的數(shù)密度波動(dòng)幅度要大于粗糙壁面,且呈現(xiàn)逐漸衰減趨勢.通過改變肋高h(yuǎn)來改變壁面粗糙度,實(shí)現(xiàn)不同粗糙表面構(gòu)造.圖4(a)結(jié)果顯示,當(dāng)h較小時(shí),近壁區(qū)流體數(shù)密度波動(dòng)呈現(xiàn)逐漸衰減趨勢,但隨著肋高的增大,近壁區(qū)流體數(shù)密度分布出現(xiàn)一次回升現(xiàn)象,這是由于凹槽內(nèi)的壁面與切面處的壁面對流體均有影響,導(dǎo)致數(shù)密度分布呈現(xiàn)衰弱、增強(qiáng)、再衰弱的趨勢.為研究肋間距a變化對通道內(nèi)流體數(shù)密度分布的影響,固定肋高h(yuǎn)=2.0σ不變,分別取長度a=2.7σ,3.6σ,4.5σ,5.4σ,結(jié)果如圖5所示.由圖5(a)和圖5(b)可知,納米通道內(nèi)流體在不同肋間距下數(shù)密度的振蕩周期一致,肋間距a變化對近壁區(qū)數(shù)密度影響較小.因此,可認(rèn)為肋間距a的變化基本不影響壁面粗糙度對流體數(shù)密度分布.

    圖3 不同勢能系數(shù)c下流體沿z方向的密度分布 (a)c = 2.0;(b)c = 1.0;(c)c = 0.6;(d)c = 0.2Fig.3.Density profiles in thez-direction with different energy coefficientc:(a)c = 2.0;(b)c = 1.0;(c)c = 0.6;(d)c = 0.2.

    圖4 肋高h(yuǎn)對壁面附近流體數(shù)密度分布的影響 (a)粗糙壁面;(b)光滑壁面Fig.4.Effect of rib heighthon the distribution of fluid number density near wall surface:(a) Rough wall surface;(b) smooth wall surface.

    4.1.2 速度分布

    圖5 肋間距a對壁面附近流體數(shù)密度分布的影響 (a)粗糙壁面;(b)光滑壁面Fig.5.Effect of rib spacingaon the distribution of fluid number density near wall surface:(a) Rough wall surface;(b) smooth wall surface.

    圖6 不同肋高h(yuǎn)下流體沿y方向的速度分布 (a) Couette流動(dòng);(b) Poiseuille流動(dòng)Fig.6.Velocity profiles in they-direction with different rib heighth:(a) Couette flow;(b) Poiseuille flow.

    壁面粗糙度變化引起的通道內(nèi)流體數(shù)密度的變化,也會(huì)導(dǎo)致通道內(nèi)流體流速分布的變化.為確定模型水動(dòng)力位置,分別計(jì)算Couette流動(dòng)及Poiseuille流動(dòng)時(shí)通道內(nèi)的速度分布,更深一步地分析邊界滑移機(jī)理.圖6(a)和圖6(b)為不同肋高對兩種流動(dòng)狀態(tài)下速度分布的影響,橫坐標(biāo)為沿y方向上的通道高度分布,縱坐標(biāo)為無量綱速度u?=u/(ε/m)1/2.圖6結(jié)果顯示,微通道內(nèi)流體速度分布呈現(xiàn)不對稱性,壁面粗糙度的存在使近壁區(qū)流體剪切產(chǎn)生了額外黏性耗散,其次粗糙壁面凹槽內(nèi)流體原子的運(yùn)動(dòng)受到限制,難以從凹槽內(nèi)逃脫進(jìn)布差異越明顯.通過兩種流動(dòng)狀態(tài)下的速度分布可推算出固液邊界位置及滑移長度,本文主要分析Poiseuille流動(dòng)時(shí)的滑移機(jī)理,通過施加不同外部驅(qū)動(dòng)力F= 0.03ε/σ,F= 0.05ε/σ,F= 0.07ε/σ來獲取滑移長度計(jì)算的標(biāo)準(zhǔn)差η,并取不同驅(qū)動(dòng)力下滑移長度的平均值作為分析結(jié)果.圖7所示為不同肋高h(yuǎn)下滑移長度標(biāo)準(zhǔn)差分布.圖8所示為不同肋高h(yuǎn)對Poiseuille流動(dòng)中通道內(nèi)滑移速度及滑移長度的影響.分析結(jié)果顯示,光滑壁面一側(cè)的滑移速度及滑移長度要高于粗糙壁面一側(cè),隨著肋高h(yuǎn)的增大,粗糙壁面一側(cè)滑移速度及滑移長度逐漸減小,與Schmatko等[24]研究粗糙高度對界面處滑移速度的影響結(jié)論一致.圖9為肋間距a變化對兩種流動(dòng)狀態(tài)下速度分布的影響.由圖9可知,兩種流動(dòng)狀態(tài)下肋間距a的變化對流場中速度分布影響較小,這是由于肋間距a變化不會(huì)產(chǎn)生額外的入主流區(qū)域,導(dǎo)致Couette流動(dòng)中壁面更容易帶動(dòng)流體運(yùn)動(dòng),速度大小分布高于光滑壁面一側(cè),而在Poiseuille流動(dòng)中則剛好相反,流體的運(yùn)動(dòng)受到壁面的阻礙,速度大小分布低于光滑壁面一側(cè).隨著肋高h(yuǎn)的增大,凹槽內(nèi)限制的流體原子數(shù)量增多,產(chǎn)生的黏性耗散越大,不同結(jié)構(gòu)壁面附近的速度分黏性耗散,且凹槽限制流體原子的數(shù)量不會(huì)產(chǎn)生太大差異.圖10為不同肋間距a下滑移長度標(biāo)準(zhǔn)差分布.圖11(a)和圖11(b)為肋間距a變化對Poiseuille流動(dòng)時(shí)通道內(nèi)滑移效應(yīng)的影響,可以看出肋間距a變化對滑移速度及滑移長度的影響較小.通過分析壁面粗糙度對速度分布的影響,發(fā)現(xiàn)肋高h(yuǎn)的變化會(huì)影響水動(dòng)力位置的變化,對邊界滑移影響較大;而肋間距a變化對通道內(nèi)速度分布影響較小,水動(dòng)力位置變化不大,故對邊界滑移影響較小.

    圖7 不同肋高h(yuǎn)下滑移長度標(biāo)準(zhǔn)差分布Fig.7.Standard deviation distribution of slip length with different rib heighth.

    圖8 (a)肋高h(yuǎn)對滑移長度的影響;(b)肋高h(yuǎn)對滑移速度的影響Fig.8.(a) Effect of rib heighthon the slip length;(b) effect of rib heighthon the slip velocity.

    圖9 不同肋間距a下流體沿y方向的速度分布 (a) Couette流動(dòng);(b) Poiseuille流動(dòng)Fig.9.Velocity profiles in they-direction with different rib spacinga:(a) Couette flow;(b) Poiseuille flow.

    圖10 不同肋間距a下滑移長度標(biāo)準(zhǔn)差分布Fig.10.Standard deviation distribution of slip length with different rib spacinga.

    4.2 壁面潤濕性變化對流動(dòng)特性的影響

    4.2.1 壁面潤濕性變化對數(shù)密度的影響

    壁面潤濕性決定了固體壁面與流體間的相互作用,不僅會(huì)影響固液界面處動(dòng)量的傳遞,也會(huì)改變近壁區(qū)流體原子的分布狀態(tài)[25].為研究壁面潤濕性對微通道內(nèi)流動(dòng)的影響,固定凹槽高度h=2.0σ,長度a= 3.6σ,對勢能系數(shù)c= 1.0,0.75,0.5,0.25下流體在通道內(nèi)的流動(dòng)特性分別進(jìn)行研究.圖12為不同潤濕性下對應(yīng)的通道內(nèi)流體數(shù)密度分布,由圖12可知,無論是粗糙壁面還是光滑壁面,壁面與流體間的作用力越強(qiáng),壁面潤濕性越好,壁面附近吸附的流體原子數(shù)量越多,近壁區(qū)流體的數(shù)密度均隨著勢能系數(shù)的增大而增大.對于粗糙壁面,當(dāng)固液界面間的勢能系數(shù)c≤ 0.5時(shí),隨著固液間勢能系數(shù)c的減小,凹槽內(nèi)流體原子的數(shù)密度有所下降,但下降幅度小于光滑壁面,且凹槽內(nèi)流體數(shù)密度要大于切面附近.而當(dāng)固液間的勢能系數(shù)c= 0.25時(shí),對應(yīng)的接觸角最大,凹槽內(nèi)流體原子的數(shù)密度明顯減少,如圖12(a)所示,凹槽內(nèi)的數(shù)密度波動(dòng)峰值要低于切面附近,說明此時(shí)凹槽對流體原子幾乎處于排斥狀態(tài).

    圖11 (a)肋間距a對滑移長度的影響;(b)肋間距a對滑移速度的影響Fig.11.(a) Effect of rib spacingaon the slip length;(b) effect of rib spacingaon the slip velocity.

    圖12 勢能系數(shù)c對壁面附近流體數(shù)密度分布的影響 (a)粗糙壁面;(b)光滑壁面Fig.12.Effect of energy coefficientcon the distribution of fluid number density near wall surface:(a) Rough wall surface;(b) smooth wall surface.

    4.2.2 壁面潤濕性變化對速度分布的影響

    圖13 不同勢能系數(shù)c下流體沿y方向的速度分布 (a) Couette流動(dòng);(b) Poiseuille流動(dòng)Fig.13.Velocity profiles in they-direction with different energy coefficientc:(a) Couette flow;(b) Poiseuille flow.

    圖13為勢能系數(shù)c變化對兩種流動(dòng)狀態(tài)下通道內(nèi)速度分布的影響.結(jié)果顯示,隨著勢能系數(shù)c的增大,固液原子間的作用力逐漸增強(qiáng),兩種流動(dòng)狀態(tài)下速度分布呈現(xiàn)相反的變化趨勢.如圖13(a)所示,Couette流動(dòng)中通道內(nèi)的速度分布隨著勢能系數(shù)c的增大而增大,而圖13(b)顯示Poiseuille流動(dòng)中通道內(nèi)的速度分布隨勢能系數(shù)c的增大而減小.另外值得一提的是,對于Couette流動(dòng)來說,粗糙壁面附近流體速度變化幅度要大于光滑壁面一側(cè),而Poiseuille流動(dòng)則剛好相反.圖14為不同勢能系數(shù)c下滑移長度標(biāo)準(zhǔn)差分布.圖15(a)和圖15(b)分別為Poiseuille流動(dòng)中通道內(nèi)的滑移速度及滑移長度分布,可以發(fā)現(xiàn),隨著勢能系數(shù)c的增大,無論是光滑壁面還是粗糙壁面,滑移速度和滑移長度分布均逐漸降低.通過分析可知,改變勢能系數(shù)c會(huì)影響通道內(nèi)的速度分布,使水動(dòng)力位置發(fā)生變化,并對邊界滑移影響較大.

    圖14 不同勢能系數(shù)c下滑移長度標(biāo)準(zhǔn)差分布Fig.14.Standard deviation distribution of slip length with different energy coefficientc.

    圖15 (a)勢能系數(shù)c對滑移長度的影響;(b)勢能系數(shù)c對滑移速度的影響Fig.15.(a) Effect of energy coefficientcon the slip length;(b) effect of energy coefficientcon the slip velocity.

    5 結(jié) 論

    本文采用分子動(dòng)力學(xué)方法研究了含粗糙壁面納米通道內(nèi)流體的流動(dòng)特性,探討了不同壁面粗糙度以及壁面潤濕性對通道內(nèi)流體的數(shù)密度和速度場分布的影響,所得結(jié)論如下.

    1)相比于光滑壁面,粗糙壁面附近流體的數(shù)密度分布較低,隨著肋高的增大,密度波動(dòng)呈現(xiàn)一次回升現(xiàn)象;改變肋間距對近壁區(qū)流體影響較小,數(shù)密度波動(dòng)趨勢基本一致;無論是光滑壁面還是粗糙壁面,增大壁面潤濕性均會(huì)使近壁區(qū)數(shù)密度波動(dòng)增大.

    2)通過分析Couette流動(dòng)和Poiseuille流動(dòng)時(shí)通道內(nèi)的速度場分布確定模型的固液邊界位置.分析結(jié)果表示,改變肋高及壁面潤濕性會(huì)明顯影響通道內(nèi)的速度分布,使固液邊界位置發(fā)生偏離,而肋間距變化對固液邊界位置影響較小.

    3)以Poiseuille流動(dòng)為主,分析了壁面粗糙度及潤濕性對界面滑移的影響.結(jié)果表明,粗糙壁面一側(cè)滑移速度及滑移長度均小于光滑壁面一側(cè),且隨著肋高及壁面潤濕性的增大,滑移速度和滑移長度逐漸減小;肋間距變化對界面滑移影響較小,滑移速度及滑移長度分布基本恒定.

    猜你喜歡
    潤濕性固液勢能
    “動(dòng)能和勢能”知識(shí)鞏固
    作 品:景觀設(shè)計(jì)
    ——《勢能》
    文化縱橫(2022年3期)2022-09-07 11:43:18
    “動(dòng)能和勢能”知識(shí)鞏固
    我國新一代首款固液捆綁運(yùn)載火箭長征六號甲成功首飛
    上海航天(2022年2期)2022-04-28 11:58:46
    “動(dòng)能和勢能”隨堂練
    分子動(dòng)力學(xué)模擬研究方解石表面潤濕性反轉(zhuǎn)機(jī)理
    等離子體對老化義齒基托樹脂表面潤濕性和粘接性的影響
    預(yù)潤濕對管道潤濕性的影響
    固液結(jié)合復(fù)合酶在保育豬日糧上的應(yīng)用研究
    廣東飼料(2016年1期)2016-12-01 03:43:00
    固液分離旋流器壁面磨損的數(shù)值模擬
    午夜激情av网站| 蜜桃在线观看..| 在线十欧美十亚洲十日本专区| aaaaa片日本免费| 叶爱在线成人免费视频播放| 欧美av亚洲av综合av国产av| 九色亚洲精品在线播放| 国产日韩欧美在线精品| 真人做人爱边吃奶动态| 中文字幕人妻熟女乱码| 国产精品二区激情视频| 国产无遮挡羞羞视频在线观看| 亚洲欧美日韩高清在线视频 | 欧美大码av| 啪啪无遮挡十八禁网站| 黄片小视频在线播放| 久久婷婷成人综合色麻豆| 免费在线观看日本一区| 中文字幕人妻熟女乱码| 精品免费久久久久久久清纯 | 国产无遮挡羞羞视频在线观看| 丝袜美足系列| 国产男女内射视频| 久久午夜亚洲精品久久| 欧美老熟妇乱子伦牲交| 欧美精品高潮呻吟av久久| 好男人电影高清在线观看| 久久精品亚洲熟妇少妇任你| 国产成人av激情在线播放| 欧美日韩亚洲国产一区二区在线观看 | 国产麻豆69| 一级,二级,三级黄色视频| 老司机影院毛片| 在线观看人妻少妇| 另类精品久久| 亚洲一区中文字幕在线| 国产99久久九九免费精品| 欧美成人免费av一区二区三区 | 成年动漫av网址| 国产成人系列免费观看| 成人国产一区最新在线观看| 一本大道久久a久久精品| 在线播放国产精品三级| 无人区码免费观看不卡 | 法律面前人人平等表现在哪些方面| 在线永久观看黄色视频| 丰满迷人的少妇在线观看| 女警被强在线播放| 国产黄色免费在线视频| 高清视频免费观看一区二区| 亚洲精品久久成人aⅴ小说| tube8黄色片| 日韩中文字幕欧美一区二区| 国产精品熟女久久久久浪| 最近最新中文字幕大全电影3 | 久久ye,这里只有精品| 99精品在免费线老司机午夜| 午夜免费鲁丝| 青青草视频在线视频观看| 久久国产精品人妻蜜桃| 国产高清videossex| 看免费av毛片| 亚洲国产欧美网| 国产精品一区二区精品视频观看| 亚洲va日本ⅴa欧美va伊人久久| 美女福利国产在线| 青青草视频在线视频观看| 中文字幕精品免费在线观看视频| 无限看片的www在线观看| 色综合婷婷激情| av天堂久久9| 国产区一区二久久| 亚洲av日韩在线播放| 国产高清激情床上av| 伊人久久大香线蕉亚洲五| 欧美人与性动交α欧美软件| 免费在线观看黄色视频的| 国产色视频综合| 丝袜美腿诱惑在线| 久久精品国产亚洲av香蕉五月 | 国产视频一区二区在线看| 国产精品1区2区在线观看. | 久久久久久久精品吃奶| 精品国产一区二区三区久久久樱花| 亚洲中文av在线| 国产在线视频一区二区| 国产欧美日韩一区二区三| 在线观看一区二区三区激情| 少妇粗大呻吟视频| 精品一区二区三区视频在线观看免费 | √禁漫天堂资源中文www| 亚洲伊人久久精品综合| 一进一出抽搐动态| 黄片大片在线免费观看| 一夜夜www| 蜜桃在线观看..| 精品视频人人做人人爽| 飞空精品影院首页| 曰老女人黄片| 老熟妇乱子伦视频在线观看| 热99国产精品久久久久久7| 无限看片的www在线观看| 香蕉久久夜色| 精品高清国产在线一区| 久久99热这里只频精品6学生| 日日夜夜操网爽| netflix在线观看网站| 精品国产乱子伦一区二区三区| 最新美女视频免费是黄的| 正在播放国产对白刺激| 久久久国产一区二区| 色婷婷av一区二区三区视频| 国产成人av教育| 99香蕉大伊视频| 桃花免费在线播放| 中文字幕色久视频| 久久这里只有精品19| 在线看a的网站| 免费人妻精品一区二区三区视频| 亚洲男人天堂网一区| 日韩欧美国产一区二区入口| 12—13女人毛片做爰片一| 亚洲av成人不卡在线观看播放网| 丰满饥渴人妻一区二区三| 国产无遮挡羞羞视频在线观看| 亚洲伊人色综图| 精品人妻在线不人妻| 日韩视频在线欧美| 亚洲第一av免费看| 亚洲av美国av| 性高湖久久久久久久久免费观看| 国产精品av久久久久免费| 国产有黄有色有爽视频| 亚洲精品一二三| 日韩欧美一区视频在线观看| 自线自在国产av| 欧美日韩黄片免| 日日爽夜夜爽网站| 亚洲午夜理论影院| 婷婷丁香在线五月| 亚洲精品自拍成人| 亚洲七黄色美女视频| 亚洲av日韩精品久久久久久密| 久久久水蜜桃国产精品网| 岛国毛片在线播放| 成人亚洲精品一区在线观看| 亚洲av第一区精品v没综合| 国产福利在线免费观看视频| 中亚洲国语对白在线视频| 亚洲国产欧美一区二区综合| a级毛片在线看网站| 亚洲一码二码三码区别大吗| 在线看a的网站| 日本a在线网址| 久久久久久久精品吃奶| 国产福利在线免费观看视频| 高清黄色对白视频在线免费看| 成人亚洲精品一区在线观看| 亚洲av日韩精品久久久久久密| 国产亚洲欧美在线一区二区| 欧美黑人欧美精品刺激| 国产午夜精品久久久久久| 精品国产一区二区三区久久久樱花| 啦啦啦 在线观看视频| 亚洲av欧美aⅴ国产| 桃红色精品国产亚洲av| 18禁黄网站禁片午夜丰满| 免费女性裸体啪啪无遮挡网站| 日本精品一区二区三区蜜桃| 免费在线观看完整版高清| 精品国产一区二区三区久久久樱花| 国产不卡一卡二| 丁香六月天网| 丁香六月天网| 免费日韩欧美在线观看| 亚洲av片天天在线观看| 免费不卡黄色视频| 国产极品粉嫩免费观看在线| 国产精品秋霞免费鲁丝片| 男人操女人黄网站| 亚洲精品自拍成人| 欧美亚洲 丝袜 人妻 在线| 久久亚洲真实| 一区在线观看完整版| av片东京热男人的天堂| 久久久久精品国产欧美久久久| 亚洲色图综合在线观看| 99国产精品一区二区三区| 又紧又爽又黄一区二区| 12—13女人毛片做爰片一| 久久精品aⅴ一区二区三区四区| 国产精品.久久久| 中文字幕另类日韩欧美亚洲嫩草| 亚洲av日韩在线播放| 欧美日韩成人在线一区二区| 嫁个100分男人电影在线观看| 91国产中文字幕| 亚洲国产av影院在线观看| 自拍欧美九色日韩亚洲蝌蚪91| 一级毛片电影观看| 久久久久精品国产欧美久久久| 脱女人内裤的视频| 日韩免费高清中文字幕av| 成年版毛片免费区| 亚洲精品国产区一区二| 97在线人人人人妻| 国产不卡av网站在线观看| 叶爱在线成人免费视频播放| 91国产中文字幕| 日韩大片免费观看网站| 夫妻午夜视频| 免费日韩欧美在线观看| 欧美精品亚洲一区二区| 午夜两性在线视频| 亚洲伊人色综图| 国产91精品成人一区二区三区 | 亚洲精品乱久久久久久| 久久热在线av| 精品国内亚洲2022精品成人 | 日韩中文字幕视频在线看片| 男女高潮啪啪啪动态图| 欧美日韩国产mv在线观看视频| 欧美日韩av久久| 天天躁夜夜躁狠狠躁躁| 国产精品久久久人人做人人爽| 国产精品九九99| 一边摸一边抽搐一进一小说 | 99精品久久久久人妻精品| 首页视频小说图片口味搜索| 2018国产大陆天天弄谢| 成人精品一区二区免费| 岛国在线观看网站| 一区二区三区激情视频| 国产野战对白在线观看| kizo精华| 高清黄色对白视频在线免费看| 人人妻人人澡人人看| 中文字幕最新亚洲高清| 亚洲自偷自拍图片 自拍| 丝瓜视频免费看黄片| 深夜精品福利| 亚洲国产精品一区二区三区在线| 狠狠狠狠99中文字幕| 欧美成狂野欧美在线观看| 国产成人免费观看mmmm| 久久 成人 亚洲| 亚洲欧美日韩另类电影网站| 亚洲欧美激情在线| 99国产极品粉嫩在线观看| videosex国产| 日本av免费视频播放| 欧美另类亚洲清纯唯美| 免费在线观看日本一区| 在线十欧美十亚洲十日本专区| 久久久精品区二区三区| 亚洲第一av免费看| 人妻 亚洲 视频| 国产一区二区激情短视频| 免费观看av网站的网址| 亚洲av国产av综合av卡| 后天国语完整版免费观看| 成人国产av品久久久| 老司机亚洲免费影院| 国产高清视频在线播放一区| avwww免费| 午夜福利免费观看在线| 嫁个100分男人电影在线观看| 男女高潮啪啪啪动态图| 麻豆乱淫一区二区| 1024香蕉在线观看| 久久人妻福利社区极品人妻图片| 天堂8中文在线网| 一区二区三区乱码不卡18| 亚洲国产欧美日韩在线播放| 无遮挡黄片免费观看| 一本一本久久a久久精品综合妖精| 午夜福利,免费看| 国产一区二区三区综合在线观看| 精品福利永久在线观看| 亚洲国产av影院在线观看| 午夜免费鲁丝| 精品国产超薄肉色丝袜足j| 日韩欧美一区视频在线观看| 久久久久久久精品吃奶| 亚洲精品乱久久久久久| 男男h啪啪无遮挡| 丰满饥渴人妻一区二区三| 久热这里只有精品99| 少妇裸体淫交视频免费看高清 | 视频区欧美日本亚洲| 两性夫妻黄色片| 黄色丝袜av网址大全| av不卡在线播放| 乱人伦中国视频| a级片在线免费高清观看视频| 精品福利永久在线观看| 国产男女内射视频| 免费日韩欧美在线观看| 自线自在国产av| 国产成人精品在线电影| 黄色片一级片一级黄色片| 久久青草综合色| 在线观看免费视频日本深夜| 亚洲第一av免费看| 日韩欧美免费精品| 一区二区av电影网| av国产精品久久久久影院| 女同久久另类99精品国产91| 亚洲精品久久成人aⅴ小说| 狠狠婷婷综合久久久久久88av| 中文字幕人妻熟女乱码| 欧美大码av| 下体分泌物呈黄色| 91成人精品电影| 久久国产精品大桥未久av| 久久久国产成人免费| 久久久久精品国产欧美久久久| 777久久人妻少妇嫩草av网站| 国产xxxxx性猛交| cao死你这个sao货| 人人妻人人添人人爽欧美一区卜| 建设人人有责人人尽责人人享有的| 国产精品国产av在线观看| 中亚洲国语对白在线视频| 高清在线国产一区| 国产精品1区2区在线观看. | 中文字幕精品免费在线观看视频| 亚洲专区中文字幕在线| 99在线人妻在线中文字幕 | 1024视频免费在线观看| 亚洲第一欧美日韩一区二区三区 | 成人av一区二区三区在线看| 宅男免费午夜| 黄色片一级片一级黄色片| 一本久久精品| 丝瓜视频免费看黄片| 亚洲av日韩精品久久久久久密| 久久久国产成人免费| 啦啦啦在线免费观看视频4| 久久婷婷成人综合色麻豆| 日本黄色日本黄色录像| 后天国语完整版免费观看| 久久影院123| 激情在线观看视频在线高清 | www.自偷自拍.com| 捣出白浆h1v1| 后天国语完整版免费观看| 50天的宝宝边吃奶边哭怎么回事| 日韩中文字幕欧美一区二区| 国产精品九九99| 中文字幕人妻熟女乱码| 怎么达到女性高潮| 亚洲专区中文字幕在线| 国产精品一区二区在线观看99| 中文字幕色久视频| 亚洲国产欧美网| 国产日韩欧美亚洲二区| 一夜夜www| 国产片内射在线| 国产在线一区二区三区精| 99精国产麻豆久久婷婷| 黄色片一级片一级黄色片| 亚洲国产欧美网| 成年动漫av网址| 中国美女看黄片| 欧美精品av麻豆av| 黑人巨大精品欧美一区二区mp4| 热99久久久久精品小说推荐| 香蕉国产在线看| 怎么达到女性高潮| 国产深夜福利视频在线观看| 天天添夜夜摸| 美女扒开内裤让男人捅视频| 午夜激情av网站| 精品少妇黑人巨大在线播放| 性少妇av在线| 色老头精品视频在线观看| 精品人妻在线不人妻| 久久精品成人免费网站| 国产在线精品亚洲第一网站| 国产无遮挡羞羞视频在线观看| 黄色视频不卡| 99riav亚洲国产免费| 色视频在线一区二区三区| 热re99久久精品国产66热6| 精品少妇黑人巨大在线播放| 亚洲人成电影免费在线| 99九九在线精品视频| 久久精品国产a三级三级三级| 高清黄色对白视频在线免费看| 2018国产大陆天天弄谢| 老熟妇仑乱视频hdxx| 久久精品国产99精品国产亚洲性色 | 亚洲精品一卡2卡三卡4卡5卡| 久久天堂一区二区三区四区| 成人国产av品久久久| 国产精品国产av在线观看| 两性夫妻黄色片| 香蕉久久夜色| 午夜福利一区二区在线看| 99热网站在线观看| 精品福利观看| 免费在线观看完整版高清| 操出白浆在线播放| av电影中文网址| 久久精品aⅴ一区二区三区四区| 亚洲熟女精品中文字幕| 两性午夜刺激爽爽歪歪视频在线观看 | 又黄又粗又硬又大视频| 久久精品国产亚洲av香蕉五月 | av福利片在线| 又黄又粗又硬又大视频| 久久中文看片网| 桃花免费在线播放| 黄色视频在线播放观看不卡| 国产av一区二区精品久久| 别揉我奶头~嗯~啊~动态视频| 一本大道久久a久久精品| 国产91精品成人一区二区三区 | 久久婷婷成人综合色麻豆| 露出奶头的视频| www.精华液| 亚洲精品在线美女| 女人久久www免费人成看片| 一本色道久久久久久精品综合| 亚洲欧美精品综合一区二区三区| 精品人妻熟女毛片av久久网站| 国产成人av教育| 两性午夜刺激爽爽歪歪视频在线观看 | 亚洲av片天天在线观看| 啦啦啦免费观看视频1| 免费人妻精品一区二区三区视频| 精品国产国语对白av| 80岁老熟妇乱子伦牲交| 性色av乱码一区二区三区2| 成人18禁在线播放| 黄网站色视频无遮挡免费观看| 中文字幕人妻丝袜一区二区| 婷婷丁香在线五月| 99久久人妻综合| 久久久久国内视频| 啦啦啦中文免费视频观看日本| tube8黄色片| 在线观看免费日韩欧美大片| 少妇猛男粗大的猛烈进出视频| 亚洲第一青青草原| 午夜福利视频精品| 国产淫语在线视频| 日韩人妻精品一区2区三区| 人人澡人人妻人| 久久久国产欧美日韩av| 母亲3免费完整高清在线观看| 99精品在免费线老司机午夜| 80岁老熟妇乱子伦牲交| 精品久久蜜臀av无| 国产精品久久久久久精品古装| 丝袜美腿诱惑在线| 国产av国产精品国产| 免费在线观看完整版高清| 正在播放国产对白刺激| 欧美精品高潮呻吟av久久| 欧美人与性动交α欧美精品济南到| 水蜜桃什么品种好| 99国产精品一区二区蜜桃av | 亚洲国产中文字幕在线视频| 精品国产亚洲在线| av有码第一页| 最近最新中文字幕大全电影3 | 狠狠狠狠99中文字幕| 国产99久久九九免费精品| 免费黄频网站在线观看国产| 又大又爽又粗| 黑人操中国人逼视频| 另类亚洲欧美激情| 欧美黄色淫秽网站| 精品国产超薄肉色丝袜足j| 91成年电影在线观看| 久久久国产一区二区| 999久久久精品免费观看国产| 日韩欧美免费精品| 久久国产精品人妻蜜桃| 涩涩av久久男人的天堂| 十八禁人妻一区二区| 午夜91福利影院| 午夜福利,免费看| 最新在线观看一区二区三区| 成人特级黄色片久久久久久久 | 亚洲av国产av综合av卡| 91九色精品人成在线观看| 国产精品国产av在线观看| 精品久久久久久电影网| 欧美精品亚洲一区二区| 俄罗斯特黄特色一大片| 天堂8中文在线网| 久久狼人影院| 制服诱惑二区| 欧美精品亚洲一区二区| 国产男女超爽视频在线观看| 免费av中文字幕在线| 国产欧美日韩一区二区三区在线| 后天国语完整版免费观看| 啦啦啦视频在线资源免费观看| 另类亚洲欧美激情| 一进一出抽搐动态| 亚洲色图av天堂| 亚洲熟妇熟女久久| 2018国产大陆天天弄谢| 亚洲全国av大片| 桃红色精品国产亚洲av| 精品国产国语对白av| 热99久久久久精品小说推荐| 老司机午夜十八禁免费视频| 在线观看66精品国产| 国产av精品麻豆| 脱女人内裤的视频| 男女无遮挡免费网站观看| 日本撒尿小便嘘嘘汇集6| 女同久久另类99精品国产91| 国产精品av久久久久免费| 亚洲情色 制服丝袜| av在线播放免费不卡| 国产精品美女特级片免费视频播放器 | 777米奇影视久久| 大型黄色视频在线免费观看| 后天国语完整版免费观看| 国产精品国产高清国产av | 一区二区日韩欧美中文字幕| 亚洲天堂av无毛| 在线十欧美十亚洲十日本专区| 国产亚洲欧美在线一区二区| 看免费av毛片| 亚洲成人国产一区在线观看| 国产男女内射视频| 一个人免费看片子| 精品熟女少妇八av免费久了| 在线亚洲精品国产二区图片欧美| 国产欧美日韩一区二区精品| 久久久久久人人人人人| 操出白浆在线播放| 飞空精品影院首页| 在线观看免费日韩欧美大片| 精品少妇一区二区三区视频日本电影| 欧美在线一区亚洲| 热99国产精品久久久久久7| 精品国产一区二区三区久久久樱花| 高潮久久久久久久久久久不卡| 他把我摸到了高潮在线观看 | 少妇粗大呻吟视频| 国产黄频视频在线观看| 水蜜桃什么品种好| 一夜夜www| 午夜福利影视在线免费观看| 性高湖久久久久久久久免费观看| 国产成人av激情在线播放| 考比视频在线观看| 色综合欧美亚洲国产小说| 精品亚洲成a人片在线观看| 99香蕉大伊视频| 每晚都被弄得嗷嗷叫到高潮| 在线观看免费日韩欧美大片| 亚洲欧美精品综合一区二区三区| 亚洲七黄色美女视频| 午夜福利影视在线免费观看| 五月开心婷婷网| 脱女人内裤的视频| 黄片小视频在线播放| 国产片内射在线| 另类亚洲欧美激情| 国产成人啪精品午夜网站| 国产在线观看jvid| 久久久久久亚洲精品国产蜜桃av| 久久人人97超碰香蕉20202| 亚洲熟女精品中文字幕| 欧美日韩亚洲高清精品| 国产成人一区二区三区免费视频网站| 国产在线一区二区三区精| 国产单亲对白刺激| 80岁老熟妇乱子伦牲交| 欧美精品av麻豆av| 久久久精品免费免费高清| 大码成人一级视频| 天天躁狠狠躁夜夜躁狠狠躁| 色婷婷av一区二区三区视频| 国产深夜福利视频在线观看| 国产欧美日韩综合在线一区二区| 悠悠久久av| 免费观看a级毛片全部| 精品国产亚洲在线| 在线观看免费日韩欧美大片| 午夜两性在线视频| 亚洲五月婷婷丁香| 久久婷婷成人综合色麻豆| 亚洲一卡2卡3卡4卡5卡精品中文| 狠狠精品人妻久久久久久综合| 亚洲人成伊人成综合网2020| 波多野结衣av一区二区av| 欧美一级毛片孕妇| 无遮挡黄片免费观看| 精品午夜福利视频在线观看一区 | 1024香蕉在线观看| 国产精品一区二区免费欧美| 欧美老熟妇乱子伦牲交| 50天的宝宝边吃奶边哭怎么回事| 日韩三级视频一区二区三区| 色综合婷婷激情| 亚洲九九香蕉| 国产精品熟女久久久久浪| 亚洲av欧美aⅴ国产| 又紧又爽又黄一区二区| 免费人妻精品一区二区三区视频| 亚洲精品在线美女| 一进一出抽搐动态| 新久久久久国产一级毛片| 大码成人一级视频| 国产欧美日韩一区二区三|