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

    不凝性氣體對納米通道內(nèi)水分子流動傳熱影響的分子動力學(xué)模擬*

    2024-05-13 07:41:28邢赫威陳占秀楊歷蘇瑤李源華呼和倉
    物理學(xué)報(bào) 2024年9期

    邢赫威 陳占秀 楊歷 蘇瑤 李源華 呼和倉

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

    隨著電子元件高性能化和小型化的發(fā)展,納米通道內(nèi)工質(zhì)的流動傳熱問題受到了更多的關(guān)注.本文采用分子動力學(xué)模擬方法,模擬了300,325,350 K 的納米通道中流體的流動傳熱情況,工質(zhì)為水,水中不凝性氣體用氬氣代替.結(jié)果表明: 流動過程中,氬原子形成高勢能團(tuán)簇,隨著溫度升高,流體勢能上升,團(tuán)簇逐漸減小或消失;少量氣體原子能夠促進(jìn)流動,而較多氬氣會導(dǎo)致通道中心區(qū)域形成較大氣體團(tuán)簇而阻礙流動,同時(shí),被加熱的工質(zhì)能顯著減小流動阻力系數(shù);近壁面區(qū)域流體溫度高于中心區(qū)域,團(tuán)簇內(nèi)部原子活動更加劇烈,平均分子動能更大,溫度更高;水的氫鍵結(jié)構(gòu)可以促進(jìn)納米通道內(nèi)的傳熱,氬原子會影響氫鍵數(shù)量,高溫會破壞水分子形成的氫鍵網(wǎng)絡(luò),使努塞爾數(shù)下降.本研究分析了不凝性氣體影響下微通道內(nèi)水分子流動傳熱的機(jī)理,為電子設(shè)備的強(qiáng)化傳熱提供了理論指導(dǎo).

    1 引言

    隨著微納米科學(xué)技術(shù)的不斷發(fā)展,電子器件的性能需求日益增加[1].微納米處理技術(shù)的進(jìn)步使計(jì)算機(jī)芯片的運(yùn)算速度大幅增大,但芯片產(chǎn)生的熱量也劇烈升高[2].為了滿足高性能芯片的冷卻需求,微納米通道中流體的散熱問題受到越來越多研究學(xué)者的重視.然而在微尺度下,物質(zhì)運(yùn)輸與能量傳遞都被限制在一個(gè)狹小的空間內(nèi),在分子力效應(yīng)、表面張力效應(yīng)等多種微尺度效應(yīng)的作用下,實(shí)驗(yàn)的測量精度受到很大的限制,宏觀尺度下流體流動規(guī)律的原理不再適用[3,4].為了更好地對納米通道進(jìn)行分析,研究者把目光轉(zhuǎn)向分子動力學(xué)方法.分子動力學(xué)方法從微觀原子運(yùn)動的角度出發(fā),運(yùn)用計(jì)算機(jī)模擬技術(shù)準(zhǔn)確地描述分子之間的作用力,并且不受時(shí)間與空間的限制,對揭示微觀層面機(jī)制提供了新的視角.因此,分子動力學(xué)方法成為當(dāng)前研究微納米通道中流動傳熱問題的主要方式[5].

    在納米尺度下,壁面效應(yīng)能夠?qū)α鞯纼?nèi)工質(zhì)的流動產(chǎn)生較大影響,很多學(xué)者對納米通道內(nèi)流體的流動情況和輸運(yùn)特性進(jìn)行了深入研究.Shi等[6]采用分子動力學(xué)方法研究了表面潤濕性與數(shù)密度之間的關(guān)系,結(jié)果表明隨著表面潤濕性的提高,密度振蕩的幅度更大,密度峰和密度谷之間的間距縮短,密度振蕩行為更加明顯.Voronov等[7]研究了流體的Couette 流動,發(fā)現(xiàn)當(dāng)固體壁面與流體之間的接觸角超過某一特定值時(shí),接觸角的增大會導(dǎo)致滑移長度的減小.Zhang 和Chen[8]研究了三維Couette 流動系統(tǒng)中的界面速度滑移,結(jié)果表明,無論液體在橫向還是縱向流動,粗糙表面都會引起額外的能量損失,并導(dǎo)致納米通道中界面速度的降低.Priezjev[9]研究了表面能和剪切速率對簡單流體流動中滑移長度的影響,發(fā)現(xiàn)在弱壁面-流體相互作用下,滑移長度隨剪切速率呈非線性關(guān)系.Semiromi 和Azimian[10]設(shè)置了不同驅(qū)動力作用下的流動模型,發(fā)現(xiàn)平均軸向速度在通道中心線附近具有對稱分布,并且外部驅(qū)動力的增大可以增大流體的最大和平均速度值,此外,還發(fā)現(xiàn)滑移長度和滑移速度是驅(qū)動力的函數(shù).

    納米通道內(nèi)流體的傳熱問題也引起了諸多學(xué)者的關(guān)注.Markvoort等[11]使用分子動力學(xué)研究了潤濕性對界面熱輸運(yùn)性質(zhì)的影響,結(jié)果表明,潤濕性的提高有利于原子之間的熱傳遞.Zhou等[12]構(gòu)建了銅原子壁面與氬原子流體間的傳熱模型,結(jié)果表明,即使對于超親水表面,固-液相互作用強(qiáng)度的增大仍然能夠顯著增強(qiáng)核沸騰傳熱性能.Bao等[13]研究了光滑熱基底約束的薄液膜中滑移長度隨溫度和速率的變化規(guī)律,發(fā)現(xiàn)對于親水壁面和疏水壁面,流體的溫度都會隨著剪切速率的上升而升高.Ge等[14]模擬熱流體流過冷壁面的納米通道,發(fā)現(xiàn)努塞爾數(shù)隨著流體與壁面之間作用力的增強(qiáng)而升高,表明對流換熱可以通過增加表面潤濕性增強(qiáng).Yao等[15]建立了恒定壁面溫度下納米通道內(nèi)的對流換熱模型,研究了固-液相互作用力對流動和傳熱特性的影響,發(fā)現(xiàn)相互作用能夠在一定范圍內(nèi)有效提高整體的傳熱性能,如果相互作用力繼續(xù)增大,流動阻力的顯著增大將削弱納米通道整體傳熱性能的提高.

    目前,對于納米通道中的流動傳熱問題,大多數(shù)研究選取了簡單原子作為流動工質(zhì).水作為實(shí)際工程應(yīng)用廣泛的流體,其流動與傳熱特性也受到了學(xué)者們的關(guān)注.Alkhwaji等[16]采用分子動力學(xué)方法對水分子的熱力學(xué)性質(zhì)進(jìn)行了研究,并與實(shí)驗(yàn)結(jié)果進(jìn)行了比較,結(jié)果表明,模擬水分子的導(dǎo)熱系數(shù)、黏度、密度、比熱容、比熱壓、普朗特?cái)?shù)和體積熱膨脹系數(shù)均于實(shí)驗(yàn)結(jié)果相當(dāng)接近,證明了分子動力學(xué)模擬的有效性.Hong等[17]通過模擬發(fā)現(xiàn)液態(tài)水分子中存在具有一定間距的納米氣泡團(tuán).Xi等[18]研究了固體壁面與水分子間的傳熱機(jī)理,結(jié)果表明,增大相互作用強(qiáng)度有利于水分子形成氫鍵網(wǎng)絡(luò),從而增強(qiáng)界面?zhèn)鳠?

    綜上所述,目前研究水分子作為工質(zhì)的納米通道流動傳熱問題還不夠詳細(xì),而考慮納米氣泡對水分子熱力學(xué)與動力學(xué)性質(zhì)影響的研究更是較少.本文采用分子動力學(xué)模擬方法,構(gòu)建了水分子摻入不凝性氣體情況下的流動傳熱模型,選用氬這一單原子分子作為氣體原子,以達(dá)到簡化計(jì)算與提高計(jì)算精度的目的.通過設(shè)置5 種不同氬原子含量的流體,在300,325,350 K 的熱源溫度下,對流體工質(zhì)的流動與傳熱性能進(jìn)行了深入研究,并通過分析氣體團(tuán)簇的變化,揭示其影響機(jī)理.

    2 仿真模型和模擬設(shè)置

    2.1 物理模型和模擬方法

    納米通道內(nèi)流體流動的模型如圖1 所示.仿真模型的尺寸為15 ? × 110 ? × 80 ? (x×y×z).納米通道壁面由銅原子構(gòu)成,按照面心立方結(jié)構(gòu)(FCC)進(jìn)行排列,晶格常數(shù)為3.62 ?.初始流動模型如圖1 所示,Z方向上最外層的兩層銅原子(綠色原子層)是固定的,厚度為3.62 ?,用于防止模型在運(yùn)算期間發(fā)生遷移或變形;臨近最外層的兩層銅原子(橙色原子層)為加熱層,厚度為3.62 ?,被用作熱恒溫器,由Langevin 恒溫器控制;內(nèi)層原子(黃色原子層)與流體接觸,起到導(dǎo)熱作用;加熱層與導(dǎo)熱層的銅原子受彈簧力約束,以便銅原子在模擬過程中交換能量時(shí)能夠維持在其原始位置振動,彈簧力大小為-Kr,其中K是彈簧常數(shù),大小為0.426 eV/nm,r是原子實(shí)際位置于原始位置之間的距離.流體區(qū)域由水分子和氬原子組成,其中紅色原子為氧原子,白色原子為氫原子,藍(lán)色原子為氬原子,均勻分布于間距為H=50 ?的固體壁面之間,流道內(nèi)的原子數(shù)量如表1 所列.流體在外力F的作用下沿Y方向流動,為提高模型的真實(shí)性和系統(tǒng)的穩(wěn)定性,在X,Y,Z方向上采用周期性邊界條件.

    表1 流道內(nèi)的粒子數(shù)量Table 1. Number of particles in the channel.

    圖1 模擬系統(tǒng)圖Fig.1.Simulation system.

    在本模擬中,水分子模型采用TIP4P 模型,該模型在化學(xué)鍵H—O—H 夾角的對角線上距離氧原子0.0125 nm 處引入一個(gè)虛擬原子,并將原本氧原子上的負(fù)電荷置于該虛擬原子上(該虛擬原子只含電量,不含質(zhì)量),以改善水分子周圍的靜電分布,正電荷位置保持不變.水分子氫氧鍵長為0.09572 nm,鍵角為 104.52°,O 原子和 H 原子所帶電荷分別為1.04ec和0.52ec,ec為單個(gè)電子所帶的電荷.

    本文采用12-6 Lennard-Jones (LJ) 勢能模型應(yīng)用于流體原子間及固-液原子間[19],其表達(dá)式為

    式中,U為原子間的相互作用力,rij為原子間的距離,εij為原子間的勢能參數(shù),σij為原子間的尺寸參數(shù),rc為截止距離.各原子之間的相互作用勢能參數(shù)及作用距離如表2 所列[20].

    表2 各原子之間的相互作用勢能參數(shù)及尺寸參數(shù)Table 2. Interaction potential energy parameters and size parameters of each atom.

    對于其他原子間的勢能參數(shù)εij與尺寸參數(shù)σij,采用Lorentz-Berthelot 混合方法進(jìn)行計(jì)算[21]:

    仿真系統(tǒng)的時(shí)間步長設(shè)置為1 fs,采用 Velocity-Verlet 算法對流體原子運(yùn)動方程進(jìn)行求解.為了使模擬計(jì)算的結(jié)果更加準(zhǔn)確,在初始時(shí)刻,為每個(gè)分子設(shè)定隨機(jī)的速度及位置.將整個(gè)系統(tǒng)在NVT 系綜下進(jìn)行弛豫,系統(tǒng)中的所有原子在280 K的溫度下隨機(jī)運(yùn)動,計(jì)算50 萬步(0.5 ns)以實(shí)現(xiàn)熱力學(xué)平衡.然后撤消整個(gè)系統(tǒng)的NVT 系綜,并將其設(shè)置為NVE 系綜,加熱層原子設(shè)置NVT 系綜,對流體區(qū)域的水分子和氬原子在Y方向上施加一個(gè)外力F=0.0005 kcal/(mol·?),繼續(xù)運(yùn)行50 萬步(0.5 ns),然后重新開始計(jì)時(shí)并運(yùn)行1000 萬步(10 ns),并對整個(gè)系統(tǒng)進(jìn)行區(qū)域劃分,統(tǒng)計(jì)數(shù)據(jù).本文中分子動力學(xué)模擬均采用LAMMPS 開源代碼實(shí)現(xiàn),系統(tǒng)可視化使用OVITO 可視化軟件完成[22].

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

    為了驗(yàn)證模型及參數(shù)設(shè)置的正確性,對納米通道中水分子的流動傳熱特性進(jìn)行驗(yàn)證.具體是采用了Xi等[18]的工作,構(gòu)建了一個(gè)33 ?×124 ?×24 ?的水分子流動傳熱模型,如圖2(a)所示.固體基底原子采用面心立方結(jié)構(gòu)進(jìn)行排列,晶格常數(shù)設(shè)置為4.086 ?,相互作用參數(shù)設(shè)置為εs=4.5608 kcal/mol和σs=0.2633 nm,熱基底與冷基底的間距D=52 ?.系統(tǒng)在NVT 系綜下達(dá)到平衡后,將熱基底的溫度設(shè)置為353 K,冷基底溫度設(shè)置為293 K,并進(jìn)行4 ns 的數(shù)據(jù)收集.模擬的水分子溫度分布如圖2(b)所示,將結(jié)果與Xi等[18]提供的研究數(shù)據(jù)進(jìn)行對比,基于冷熱源的60 K 溫差尺度下,最大相對誤差約為1.32%,因此,可認(rèn)為本文構(gòu)建的模型與算法準(zhǔn)確可靠.

    圖2 水分子流動傳熱模擬驗(yàn)證(a) 模型示意圖;(b) 溫度分布Fig.2.Validation of water molecule flow heat transfer simulation: (a) Model Schematic;(b) temperature distribution.

    3 結(jié)果與分析

    3.1 不凝性氣體對勢能分布的影響

    為了深入探究納米通道內(nèi)不凝性氣體對流動傳熱的影響,首先計(jì)算了流體的勢能分布.勢能可以直觀顯示出原子在納米通道中的排布情況,同時(shí)反映原子之間的相互作用能.流體在300,325,350 K 熱源溫度作用下的勢能分布如圖3 所示,Case 1 中不含不凝性氣體,通道內(nèi)水分子的勢能分布相對均勻,而對于Case 2,Case 3,Case 4,Case 5,由于不凝性氣體的存在,流體在流動的過程中,氬原子會形成團(tuán)簇,并集中分布在納米通道的中心區(qū)域,其表現(xiàn)為一個(gè)明顯的高勢能區(qū)域,Case 2 中氬原子含量最多,形成的團(tuán)簇更大,高勢能區(qū)域也更大.固-液界面處始終存在著低勢能區(qū)域(深藍(lán)色區(qū)域),使更多的流體原子被壁面捕獲.在相同的熱源溫度下,Case 1 的低勢能區(qū)相較于其余4 種情況更加顯著,這意味著氬原子的存在削弱了固體壁面與水分子之間的相互作用力.當(dāng)T0=300 K 時(shí),流體的勢能較低,氬原子形成的團(tuán)簇更大;而隨著熱源溫度的上升,流體的勢能升高,氬原子團(tuán)簇逐漸變小,并且熱源溫度達(dá)到350 K 時(shí),Case 3,Case 4,Case 5 中的氬原子團(tuán)簇基本消失,表明熱源溫度的升高可以顯著增加原子運(yùn)動的劇烈程度,使納米通道中原子整體分布趨向均勻.

    圖3 流體的勢能分布Fig.3.Potential energy distribution of a fluid.

    數(shù)密度反映了通道內(nèi)水分子和氬原子的分布情況,對流體的勢能分布情況作出進(jìn)一步的解釋,分別統(tǒng)計(jì)了水分子和氬原子在納米通道內(nèi)的數(shù)密度分布,如圖4 所示.圖中,橫坐標(biāo)為納米通道的Z方向,單位為?,縱坐標(biāo)為水分子數(shù)密度,單位為?-3.水分子和氬原子的數(shù)密度在近壁面處均出現(xiàn)了明顯的振蕩現(xiàn)象,這意味著壁面附近流體的分層現(xiàn)象更加顯著;壁面捕獲了更多的水分子,從而形成了圖3 中的低勢能區(qū).圖4(a),(b)展示了300 K 熱源溫度下水分子和氬原子的分布情況,水分子在壁面附近存在高密度層,而氬原子更多地集中在了中心區(qū)域.Case 1 的情況下由于納米流道中不存在氬原子,其中心區(qū)域數(shù)密度基本保持恒定;Case 2中水分子的數(shù)量最少,同時(shí)水的數(shù)密度最低;Case 3,Case 4,Case 5 水分子數(shù)密度在中心區(qū)域基本相同,而對于壁面處的高密度層,數(shù)密度峰值分別為0.0738,0.0743,0.0751,這意味著氬原子數(shù)量的增加可以削弱壁面對水分子的吸引力.圖4(c),(e)分別對應(yīng)熱源溫度T0=325,350 K 的情況.通過對比圖4(a),(c),(e)可知,隨著熱源溫度的上升,流體的分層現(xiàn)象被削弱,出現(xiàn)這種現(xiàn)象的原因是高溫使納米通道內(nèi)分子運(yùn)動更加劇烈,導(dǎo)致壁面吸附水分子的能力變?nèi)?水分子有更大的可能性脫離壁面的束縛.圖4(b),(d),(f)展現(xiàn)了不同溫度下氬原子數(shù)密度的變化規(guī)律,隨著熱源溫度的升高,氬原子在壁面附近的分層現(xiàn)象也隨之減弱,整體的數(shù)密度分布更加均勻.

    圖4 流體的數(shù)密度分布(a) T0=300 K,H2O;(b) T0=300 K,Ar;(c) T0=325 K,H2O;(d) T0=325 K,Ar;(e) T0=350 K,H2O;(f) T0=350 K,ArFig.4.Number density distribution of fluids: (a) T0=300 K,H2O;(b) T0=300 K,Ar;(c) T0=325 K,H2O;(d) T0=325 K,Ar;(e) T0=350 K,H2O;(f) T0=350 K,Ar.

    3.2 不凝性氣體對流體速度的影響

    納米通道內(nèi)流體的速度分布是研究流體流動特性的重要參數(shù).當(dāng)實(shí)際流體流過固體壁面時(shí),沿著Y方向的速度可以通過的Navier-Stokes 方程來確定,其表達(dá)式為

    其中,ρ(y,z) 是流體的質(zhì)量密度,μ是流體的剪切黏度,uf,Y是沿Y方向流體的流動速度,是沿Y方向的壓力梯度.在本模擬中,流速在Y方向上保持不變,因此,(5)式可以簡化為[23]

    圖5 為5 種類流體在300,325,350 K 的熱源溫度下的速度分布.橫坐標(biāo)為Z方向,單位為?,縱坐標(biāo)為流體在Y方向的平均流速,單位為m/s,流體速度分布結(jié)果采用二次多項(xiàng)式函數(shù)進(jìn)行擬合.納米通道內(nèi)流體的速度分布呈泊肅葉流動,中心的主流區(qū)域流速較高,近壁面區(qū)域由于受到固體壁面的約束導(dǎo)致流速較低.從圖5(a)可以看出,Case 5 的情況下,流體的平均速度最高,超過了純水情況下的Case 1,而氬原子含量最高的Case 2 平均流速最低.這意味著少量氣體原子的存在可以削弱固體壁面對水分子高密度層的吸引力,從而起到促進(jìn)流動的作用;而更多氣體原子的出現(xiàn)導(dǎo)致中心區(qū)域形成更大的團(tuán)簇,使納米通道內(nèi)流體的流動受到一定的阻礙.對比圖5(a),(b),(c)可得,隨著熱源溫度T0從300 K 升高至350 K,流體流動速度明顯上升,這是因?yàn)闇囟鹊纳邥率狗肿舆\(yùn)動更加劇烈,同時(shí),固體壁面對流體的吸附能力也被削弱.圖5(c)中,在T0=350 K 的情況下,Case 1 的速度有超過Case 5 的趨勢,這表明了高溫對固-液相互作用的削弱能力相較于不凝性氣體而言更強(qiáng),少量不凝性氣體為流動帶來的促進(jìn)作用在高溫下已經(jīng)很小.

    圖5 流體的速度分布(a) T0=300 K;(b) T0=325 K;(c) T0=350 KFig.5.Velocity distribution of fluids: (a) T0=300 K;(b) T0=325 K;(c) T0=350 K.

    圖6 統(tǒng)計(jì)了納米通道內(nèi)流體邊界區(qū)域的速度滑移.從圖6 可以看出,流體的速度滑移始終為正滑移,且更高的熱源溫度會產(chǎn)生更大的速度滑移.Case 3,Case 4,Case 5 的速度滑移在T0=300,325,350 K 的情況下均大于Case 1,這種現(xiàn)象說明了不凝性氣體對水分子高密度層有一定的吸引作用,從而減小了邊界處的流動阻力;而Case 2 中更多的氣體原子使得氬原子團(tuán)簇結(jié)構(gòu)對流體流動產(chǎn)生了阻礙作用,導(dǎo)致速度滑移更低.通過對比Case 2,Case 3,Case 4,Case 5 可以發(fā)現(xiàn),在水中加入少量的氣體原子可以削弱通道壁面的約束力,使得邊界區(qū)域的速度滑移更大,而更多氣體原子的加入會阻礙微通道內(nèi)的流動.此外,在高溫的情況下,納米通道內(nèi)的流體流動的速度滑移顯著增大.

    圖6 流體的速度滑移Fig.6.Velocity slip of fluids.

    計(jì)算了存在不凝性氣體情況下流體的的流動阻力系數(shù).流動阻力系數(shù)是影響流動的重要因素,從流動阻力的角度分析流動的影響是有說服力的.阻力系數(shù)f可定義為

    其中,ρ是流體的質(zhì)量密度,Δp表示沿著通道流動方向長度為LY的壓降,um是通道橫截面處的平均流體速度,Dh是水力直徑.當(dāng)水分子流過納米通道時(shí),由于在Y方向上對N個(gè)水分子中的每一個(gè)施加外力fY,使得流體被驅(qū)動.因此,fY的大小對應(yīng)于Δp·A,在分子動力學(xué)模擬中可以推導(dǎo)出以下公式進(jìn)行計(jì)算[24]:

    根據(jù)(5)式,計(jì)算了上述各種情況下的阻力系數(shù),如圖7 所示.流動阻力系數(shù)f與平均流體速度的平方呈反比關(guān)系,這意味著阻力系數(shù)也會受熱源溫度的影響.當(dāng)熱源溫度相同時(shí),氬原子含量最高的Case 2 的阻力系數(shù)始終最大,而含有最少氬原子的Case 5 阻力系數(shù)最小,Case 3,Case 4 與純水情況下的Case 1 阻力系數(shù)相差不大.隨著熱源溫度T0從300 K 升高至350 K,5 種情況下的阻力系數(shù)均顯著減小,這意味著高溫能夠有效降低微通道內(nèi)的流動阻力.值得一提的是,當(dāng)T0=350 K時(shí),Case 1 的阻力系數(shù)減小至50.7,低于Case 4 情況下的65.5,且與Case 5 情況下的45.2 已經(jīng)相差很小.這表明在高溫的情況下,少量不凝性氣體原子的存在對微通道內(nèi)流體流動的促進(jìn)作用受到削弱,而熱源溫度對流體阻力系數(shù)的影響則更加顯著.

    圖7 流體的阻力系數(shù)Fig.7.Drag coefficient of fluids.

    3.3 不凝性氣體對流體傳熱的影響

    為了探究不凝性氣體對納米通道內(nèi)傳熱特性的影響,計(jì)算了5 種流體在熱源溫度T0分別為300,325,350 K 的情況下的溫度分布.在分子動力學(xué)模擬中,流體的溫度基于動力學(xué)輸出,并通過如下公式得到[25]:

    式中,vi為單個(gè)水分子的速度,為宏觀流體的平均速度,KB為玻爾茲曼常數(shù),為納米通道中水分子的數(shù)量,m為水分子的質(zhì)量.

    納米通道內(nèi)流體的溫度分布如圖8 所示.圖8(a)展示了300 K 熱源溫度下流體的溫度云圖.可以觀察到,近壁面區(qū)域流體的溫度略高于中心區(qū)域;Case 1 情況下流體的平均溫度最高,溫度分布也最為均勻;對于Case 2,Case 3,Case 4,Case 5 這4 種含有氬原子流體的情況,流體的平均溫度比純水更低,但在納米團(tuán)簇處會形成一個(gè)高溫區(qū)域,這是因?yàn)閳F(tuán)簇內(nèi)部氬原子活動更加劇烈,平均分子動能更大.T0=325 K 的情況如圖8(b)所示,隨著熱源溫度的升高,氬原子形成的團(tuán)簇更小,雖然流體的溫度隨之升高,但平均溫度與熱源溫度之間的溫差變大,意味著流體的傳熱能力受到影響.圖8(c)展示了350 K 熱源溫度下流體溫度分布情況,可以明顯地看出,納米團(tuán)簇的減小或消失使得溫度的分布更加均勻,且流體與壁面之間的溫差更大,這一現(xiàn)象表明高溫使得流體原子運(yùn)動更加劇烈,同時(shí)傳熱能力產(chǎn)生惡化.

    圖8 流體的溫度分布(a) T0=300 K;(b) T0=325 K;(c) T0=350 KFig.8.Temperature distribution of fluids: (a) T0=300 K;(b) T0=325 K;(c) T0=350 K.

    水的氫鍵結(jié)構(gòu)可以影響納米通道內(nèi)的傳熱作用[26],圖9 展示了氫鍵結(jié)構(gòu)的示意圖.使用VMD軟件對氫鍵數(shù)量進(jìn)行了統(tǒng)計(jì),氫鍵的判定依據(jù)是兩個(gè)水分子的兩個(gè)氧原子之間的距離小于3.5 ?,并且連接兩個(gè)氧原子的矢量與氧和氫原子的矢量之間的角度小于35°.從熱運(yùn)動的角度看,氫鍵網(wǎng)絡(luò)的形成可以更好地將一側(cè)水分子的動能傳遞到另一側(cè),有利于流體的熱傳遞.

    圖9 氫鍵示意圖Fig.9.Hydrogen bonding schematic.

    圖10 展示了300 K 熱源溫度下納米通道內(nèi)的氫鍵數(shù)量.在模擬開始時(shí),Case 1 的水分子數(shù)量最多,Case 2 水分子最少,而Case 3,Case 4,Case 5水分子數(shù)量相同,因此氫鍵數(shù)量也呈現(xiàn)出相同的規(guī)律.從圖10(a)可以看出: 隨著模擬過程的進(jìn)行,納米通道內(nèi)流體的氫鍵數(shù)量逐漸收斂,流體的流動在3 ns 時(shí)趨于穩(wěn)定;相較于Case 3 和Case 4,Case 5 收斂的數(shù)值最小,這是因?yàn)镃ase 5 內(nèi)總原子含量更少,流動趨于穩(wěn)定后,流體原子的平均間距更大,這不利于氫鍵的形成;另外,統(tǒng)計(jì)了Case 3在不同熱源溫度的情況下氫鍵數(shù)量的變化,如圖10(a)插圖所示,當(dāng)熱源溫度上升時(shí),流體的氫鍵數(shù)量明顯減少,該現(xiàn)象表明高溫會破壞水分子形成的氫鍵網(wǎng)絡(luò),這意味著熱源溫度T0升高時(shí),雖然水分子吸收的總熱量會增加,但它的傳熱能力也會受到一定的限制.圖10(b) 統(tǒng)計(jì)了300 K 溫度下流動穩(wěn)定后平均單個(gè)水分子形成的氫鍵數(shù)量,Case 1中平均單個(gè)水分子形成的氫鍵數(shù)量最多,Case 2形成氫鍵最少,Case 3 與Case 4 略高于Case 5,這一現(xiàn)象表明納米通道中水分子數(shù)量的增加能夠提高平均單個(gè)水分子形成的氫鍵數(shù)量,而氬原子對氫鍵的影響幅度較小.

    圖10 流體的氫鍵數(shù)量(a)氫鍵數(shù)量;(b)平均單個(gè)水分子形成的氫鍵數(shù)量Fig.10.Number of hydrogen bonds in fluids: (a) Number of hydrogen bonds;(b) the average number of hydrogen bonds formed by a single water molecule.

    為了探究納米通道內(nèi)的綜合流動傳熱能力,計(jì)算了努塞爾數(shù)Nu,Nu是衡量對流換熱強(qiáng)度的無量綱數(shù).它表示流體底層的熱傳導(dǎo)阻力與對流熱傳遞阻力之比.計(jì)算公式如下[27]:

    其中hc是對流傳熱系數(shù),Dh是特征尺寸,λ是流體的熱導(dǎo)率.對流換熱系數(shù)hc可由局部平均溫度Tave計(jì)算:

    其中,TW是壁面溫度,Tave是沿Z方向的局部平均溫度,是使用二次多項(xiàng)式函數(shù)擬合固體-液體界面處流體溫度分布而獲得的流體溫度梯度.結(jié)合(9)式和(10)式可以得到:

    式中,H是納米通道的高度.

    圖11 顯示了5 種情況下流體的努塞爾數(shù)(Nu).可以觀察到,Nu受氬原子含量和熱源溫度影響,熱源溫度的升高會導(dǎo)致Nu小幅下降.Case 1 的情況下,更多的水分子形成了氫鍵網(wǎng)絡(luò),有利于強(qiáng)化流體的流動傳熱,因此Nu最大;Case 2 水分子相對較少致使Nu偏低;Case 3 相較于Case 4,Case 5,納米通道內(nèi)原子排列更加緊密,傳熱效果更好.值得一提的是,在T0=350 K 時(shí),Case 2的Nu達(dá)到了3.51,超過了Case 5 的3.48,這意味著高溫破壞了水分子的氫鍵網(wǎng)絡(luò),在這種情況下,Case 2 有更多的流體原子參與了傳熱,因此其努塞爾數(shù)更大.總的來說,氬原子的出現(xiàn)會對納米通道內(nèi)流體的流動傳熱產(chǎn)生一定的影響,但從整體上來看,Nu的變化幅度較小;高溫的情況下,水分子形成的氫鍵網(wǎng)絡(luò)受到影響,升高熱源溫度雖然可以使流體吸收的總熱量增大,但會小幅降低流體的傳熱效率.

    圖11 流體的努塞爾數(shù)Fig.11.Nussle number of the fluid.

    4 結(jié)論

    本文采用分子動力學(xué)模擬的方法,研究了含有不凝性氣體的水分子在納米通道中的流動傳熱特性,探討了在300,325,350 K 熱源溫度作用下氬原子對流體密度、速度和溫度的影響,具體結(jié)論如下:

    1) 流體在流動的過程中,不凝性氣體氬原子易形成團(tuán)簇居于流道的中心區(qū)域,并且具有較高的勢能;固-液界面處存在低勢能區(qū)域,并且氬原子的存在使低勢能區(qū)域減小;熱源溫度的升高可以顯著增加原子運(yùn)動的劇烈程度,使納米通道中流體平均勢能升高;水分子和氬原子的數(shù)密度在近壁面處的分層現(xiàn)象變得顯著,水分子在壁面處會出現(xiàn)高密度層,氬原子集中分布在中心區(qū)域;隨著熱源溫度的上升,水分子數(shù)密度峰值降低,壁面吸附水分子的能力變?nèi)?氬原子分層現(xiàn)象也隨之減弱,原子整體分布趨向均勻.

    2) 少量氣體原子的存在可以促進(jìn)流體的流動,而更多氣體原子的出現(xiàn)會導(dǎo)致中心區(qū)域形成團(tuán)簇,從而對流動產(chǎn)生阻礙作用;高溫同樣會使得納米通道內(nèi)的流體流動速度與速度滑移更大;阻力系數(shù)受氬原子含量和熱源溫度影響,氬原子含量最高的情況下阻力系數(shù)始終最大,而含有少量氬原子對流體的影響較小,隨著熱源溫度升高,阻力系數(shù)顯著減小.

    3) 近壁面區(qū)域流體的溫度略高于中心區(qū)域,但納米團(tuán)簇內(nèi)部氬原子活動更加劇烈,平均分子動能更大,溫度較高;隨著溫度的上升,納米通道內(nèi)流體的勢能升高,氬原子形成的團(tuán)簇逐漸變小或消失.水的氫鍵結(jié)構(gòu)可以影響納米通道內(nèi)的傳熱,隨著模擬進(jìn)行,納米通道內(nèi)流體的氫鍵數(shù)量逐漸收斂,流動趨于穩(wěn)定;氬原子會對氫鍵數(shù)量和流動傳熱產(chǎn)生影響,但從整體上來看,努塞爾數(shù)的變化幅度較小;當(dāng)熱源溫度上升時(shí),流體的氫鍵數(shù)量明顯減少,這意味著高溫會破壞水分子形成的氫鍵網(wǎng)絡(luò),流體的努塞爾數(shù)下降,從而導(dǎo)致流體的傳熱效率受到負(fù)面影響.

    国产伦理片在线播放av一区| 欧美区成人在线视频| 日韩亚洲欧美综合| 亚洲,一卡二卡三卡| 日韩强制内射视频| 啦啦啦在线观看免费高清www| 草草在线视频免费看| 这个男人来自地球电影免费观看 | 曰老女人黄片| 亚洲色图综合在线观看| 亚洲情色 制服丝袜| 亚洲欧洲精品一区二区精品久久久 | 亚洲国产欧美在线一区| 七月丁香在线播放| 婷婷色综合大香蕉| 少妇被粗大猛烈的视频| 久久国产乱子免费精品| 亚洲精品乱码久久久久久按摩| 久久久久久久久久久久大奶| 亚洲精品aⅴ在线观看| 美女大奶头黄色视频| 国产成人精品无人区| 最近中文字幕2019免费版| 欧美最新免费一区二区三区| 色94色欧美一区二区| 伦理电影大哥的女人| 国产一区二区三区综合在线观看 | 亚洲真实伦在线观看| 永久免费av网站大全| 欧美精品人与动牲交sv欧美| 高清午夜精品一区二区三区| 免费av不卡在线播放| 国产av精品麻豆| 99精国产麻豆久久婷婷| 大香蕉久久网| 欧美日韩综合久久久久久| 精品酒店卫生间| 国产成人精品福利久久| 亚洲中文av在线| 久久久久国产网址| 欧美变态另类bdsm刘玥| 赤兔流量卡办理| 亚洲国产精品999| 日日摸夜夜添夜夜添av毛片| 国产精品不卡视频一区二区| 少妇的逼水好多| 国产在线男女| 97超碰精品成人国产| 尾随美女入室| 日韩,欧美,国产一区二区三区| 欧美激情极品国产一区二区三区 | 啦啦啦中文免费视频观看日本| 日本av手机在线免费观看| 日本av手机在线免费观看| 国产成人免费无遮挡视频| 精品少妇内射三级| 2021少妇久久久久久久久久久| 下体分泌物呈黄色| 女性被躁到高潮视频| 婷婷色综合www| 我要看黄色一级片免费的| 亚洲精品日本国产第一区| 丝袜脚勾引网站| 一级黄片播放器| 亚洲国产成人一精品久久久| 丝瓜视频免费看黄片| 日本91视频免费播放| 久久毛片免费看一区二区三区| 国产成人一区二区在线| 麻豆成人av视频| 国产一区二区三区综合在线观看 | 亚洲一级一片aⅴ在线观看| 久久久久久伊人网av| 久久久久精品久久久久真实原创| 人人妻人人澡人人看| 最近最新中文字幕免费大全7| 2018国产大陆天天弄谢| 国产日韩一区二区三区精品不卡 | 99热网站在线观看| 精品午夜福利在线看| 国产精品三级大全| 一级毛片久久久久久久久女| 久久 成人 亚洲| 精品久久久精品久久久| 国产片特级美女逼逼视频| 国产精品一区www在线观看| 99热这里只有是精品在线观看| 曰老女人黄片| 水蜜桃什么品种好| 亚洲精品日本国产第一区| 日韩一区二区三区影片| 国产一级毛片在线| 久久久久久久大尺度免费视频| 久久精品国产鲁丝片午夜精品| 亚洲三级黄色毛片| 日本黄色日本黄色录像| av福利片在线观看| 日韩欧美 国产精品| 日韩av不卡免费在线播放| 91aial.com中文字幕在线观看| 七月丁香在线播放| 免费久久久久久久精品成人欧美视频 | 日本午夜av视频| 国产亚洲91精品色在线| 在线观看国产h片| 色视频在线一区二区三区| 大片免费播放器 马上看| 久久久久网色| 欧美3d第一页| 又爽又黄a免费视频| 新久久久久国产一级毛片| 亚洲精品一二三| av有码第一页| 精品人妻偷拍中文字幕| 亚洲欧美精品专区久久| 久久久久久伊人网av| 噜噜噜噜噜久久久久久91| 亚洲高清免费不卡视频| 国产黄色视频一区二区在线观看| 亚洲av不卡在线观看| av在线老鸭窝| 久久av网站| 国产黄片美女视频| 成人毛片a级毛片在线播放| 国产真实伦视频高清在线观看| 色94色欧美一区二区| 男女免费视频国产| 69精品国产乱码久久久| 18禁裸乳无遮挡动漫免费视频| 亚洲av电影在线观看一区二区三区| 99久久综合免费| 51国产日韩欧美| 久久久久久久久久久丰满| 亚洲精品一二三| 精品国产露脸久久av麻豆| 国产精品成人在线| 我的女老师完整版在线观看| 一级毛片久久久久久久久女| 日本爱情动作片www.在线观看| 啦啦啦视频在线资源免费观看| 久久国内精品自在自线图片| 日韩 亚洲 欧美在线| 亚洲国产欧美日韩在线播放 | 人人妻人人澡人人看| 各种免费的搞黄视频| 国产综合精华液| 一级毛片我不卡| 欧美xxxx性猛交bbbb| 国产乱人偷精品视频| 日本欧美视频一区| 国产伦理片在线播放av一区| 一级毛片黄色毛片免费观看视频| 免费高清在线观看视频在线观看| 国产精品国产三级国产av玫瑰| 久久综合国产亚洲精品| 日韩免费高清中文字幕av| 嘟嘟电影网在线观看| 欧美一级a爱片免费观看看| 看免费成人av毛片| 欧美精品一区二区大全| 少妇丰满av| 9色porny在线观看| 中文乱码字字幕精品一区二区三区| 精品亚洲乱码少妇综合久久| 久久久久久久久久成人| 国产女主播在线喷水免费视频网站| 成人国产麻豆网| 久久ye,这里只有精品| 2021少妇久久久久久久久久久| 亚洲国产精品国产精品| 人妻一区二区av| 天天操日日干夜夜撸| 草草在线视频免费看| 久久99热这里只频精品6学生| 亚洲国产色片| 特大巨黑吊av在线直播| 欧美高清成人免费视频www| 99精国产麻豆久久婷婷| 男的添女的下面高潮视频| 99热网站在线观看| 日日摸夜夜添夜夜爱| 久久精品夜色国产| 国产亚洲av片在线观看秒播厂| 久久亚洲国产成人精品v| 国产精品.久久久| 国产一区二区三区av在线| 男女边吃奶边做爰视频| 国产亚洲精品久久久com| 在线亚洲精品国产二区图片欧美 | 在线播放无遮挡| 永久免费av网站大全| 九九久久精品国产亚洲av麻豆| 久久99热这里只频精品6学生| a 毛片基地| 国产亚洲av片在线观看秒播厂| 色94色欧美一区二区| 成人黄色视频免费在线看| 91久久精品国产一区二区三区| 国产免费一区二区三区四区乱码| 晚上一个人看的免费电影| h视频一区二区三区| 高清午夜精品一区二区三区| 亚洲av成人精品一二三区| 视频中文字幕在线观看| 久久久久久久精品精品| 亚洲国产av新网站| 99国产精品免费福利视频| 亚洲av欧美aⅴ国产| 五月玫瑰六月丁香| 人妻一区二区av| 久久久久久久久久成人| 永久免费av网站大全| 视频区图区小说| 在线观看美女被高潮喷水网站| a 毛片基地| 精品视频人人做人人爽| av天堂中文字幕网| av福利片在线| 99九九线精品视频在线观看视频| 熟女人妻精品中文字幕| 国产免费又黄又爽又色| 99精国产麻豆久久婷婷| 精品久久久久久久久亚洲| 欧美少妇被猛烈插入视频| 男人爽女人下面视频在线观看| 中文在线观看免费www的网站| 久久久国产精品麻豆| 精华霜和精华液先用哪个| 在线观看美女被高潮喷水网站| 日韩 亚洲 欧美在线| 夜夜爽夜夜爽视频| 熟女电影av网| 草草在线视频免费看| av在线播放精品| 免费在线观看成人毛片| 美女中出高潮动态图| av天堂中文字幕网| 免费观看a级毛片全部| 国产欧美日韩精品一区二区| 伦理电影大哥的女人| 观看av在线不卡| 9色porny在线观看| 久久精品久久久久久久性| 日韩成人伦理影院| 亚洲精品亚洲一区二区| 在线观看三级黄色| 男男h啪啪无遮挡| 亚洲国产成人一精品久久久| 一个人看视频在线观看www免费| 在线观看免费日韩欧美大片 | 亚洲欧洲日产国产| 午夜激情福利司机影院| 成人综合一区亚洲| 欧美激情极品国产一区二区三区 | 色视频www国产| 国产在线男女| 国内少妇人妻偷人精品xxx网站| 最近最新中文字幕免费大全7| 毛片一级片免费看久久久久| 国产69精品久久久久777片| 午夜免费观看性视频| 欧美精品人与动牲交sv欧美| 少妇的逼好多水| 日本欧美国产在线视频| 久久久国产欧美日韩av| 肉色欧美久久久久久久蜜桃| 色婷婷av一区二区三区视频| 99九九在线精品视频 | 国产欧美日韩综合在线一区二区 | 在线亚洲精品国产二区图片欧美 | 女性生殖器流出的白浆| 久久精品久久久久久噜噜老黄| 免费观看无遮挡的男女| 高清毛片免费看| 国产精品.久久久| 国产午夜精品久久久久久一区二区三区| a级毛片免费高清观看在线播放| 中文乱码字字幕精品一区二区三区| 一级爰片在线观看| 国产av码专区亚洲av| 免费不卡的大黄色大毛片视频在线观看| 亚洲精品国产av成人精品| 91成人精品电影| 久久久久久久久久久免费av| 视频区图区小说| 最后的刺客免费高清国语| 日本91视频免费播放| 大话2 男鬼变身卡| 久久久久视频综合| 国产片特级美女逼逼视频| 中文欧美无线码| 久久鲁丝午夜福利片| 乱码一卡2卡4卡精品| 一本大道久久a久久精品| 午夜福利视频精品| 日日爽夜夜爽网站| 91久久精品国产一区二区成人| 黄色怎么调成土黄色| 伦理电影免费视频| 国产男人的电影天堂91| 日韩大片免费观看网站| 桃花免费在线播放| 亚洲国产精品专区欧美| 26uuu在线亚洲综合色| 街头女战士在线观看网站| 国产精品成人在线| 国产伦精品一区二区三区视频9| 久久精品久久久久久噜噜老黄| 欧美日韩在线观看h| 亚洲精品,欧美精品| 午夜av观看不卡| 成人国产麻豆网| 七月丁香在线播放| 亚洲av男天堂| 午夜激情久久久久久久| 最新中文字幕久久久久| 国产极品天堂在线| 亚洲在久久综合| 国产又色又爽无遮挡免| 久久久久久人妻| 中文字幕人妻熟人妻熟丝袜美| 国产成人精品一,二区| 老女人水多毛片| 免费看不卡的av| 日韩,欧美,国产一区二区三区| 蜜桃在线观看..| 久久久久国产网址| 99精国产麻豆久久婷婷| 一区二区三区精品91| 亚洲性久久影院| 亚洲,欧美,日韩| 欧美日韩综合久久久久久| 久久国产精品大桥未久av | 成人毛片60女人毛片免费| 国产男女超爽视频在线观看| 伦精品一区二区三区| 日韩视频在线欧美| h视频一区二区三区| 秋霞伦理黄片| 欧美97在线视频| 十八禁高潮呻吟视频 | 不卡视频在线观看欧美| 久久女婷五月综合色啪小说| 色94色欧美一区二区| 在线看a的网站| 亚洲av欧美aⅴ国产| 日韩欧美 国产精品| 22中文网久久字幕| 欧美日韩在线观看h| 91aial.com中文字幕在线观看| 在线观看免费日韩欧美大片 | 伦精品一区二区三区| 成人亚洲欧美一区二区av| 亚洲欧美精品专区久久| 色哟哟·www| 久久国内精品自在自线图片| 99久久人妻综合| 精品酒店卫生间| 久久国产乱子免费精品| 日韩成人av中文字幕在线观看| 三级国产精品片| 国产伦理片在线播放av一区| 在线观看av片永久免费下载| 又粗又硬又长又爽又黄的视频| 日韩亚洲欧美综合| 国产在线一区二区三区精| 国产成人一区二区在线| 久久久久国产网址| 国产精品偷伦视频观看了| 最近手机中文字幕大全| 少妇裸体淫交视频免费看高清| 亚洲欧洲日产国产| av卡一久久| 人妻少妇偷人精品九色| 国产中年淑女户外野战色| 久久久久久久大尺度免费视频| 亚洲av成人精品一区久久| 久久av网站| 日韩av不卡免费在线播放| 最黄视频免费看| 噜噜噜噜噜久久久久久91| 美女cb高潮喷水在线观看| 国产无遮挡羞羞视频在线观看| 伊人久久国产一区二区| 国产又色又爽无遮挡免| 黄色配什么色好看| 永久免费av网站大全| 日本色播在线视频| 久久ye,这里只有精品| a 毛片基地| 国产69精品久久久久777片| 亚洲精品一区蜜桃| av专区在线播放| 欧美区成人在线视频| 黄片无遮挡物在线观看| 多毛熟女@视频| av女优亚洲男人天堂| 老熟女久久久| 成人午夜精彩视频在线观看| www.av在线官网国产| 亚洲精品视频女| 中国三级夫妇交换| 最后的刺客免费高清国语| 亚洲av不卡在线观看| 91久久精品国产一区二区成人| 久久 成人 亚洲| 麻豆成人午夜福利视频| 亚洲精品456在线播放app| 免费少妇av软件| 国产成人精品福利久久| 欧美精品人与动牲交sv欧美| 精品国产一区二区久久| 久久影院123| 国产日韩欧美视频二区| 国产午夜精品一二区理论片| 精品国产露脸久久av麻豆| 国产黄频视频在线观看| 中文欧美无线码| 我要看黄色一级片免费的| 日韩欧美 国产精品| 97超视频在线观看视频| 亚洲色图综合在线观看| 亚洲精品中文字幕在线视频 | 亚洲av中文av极速乱| 国产亚洲欧美精品永久| 日本欧美国产在线视频| 国产亚洲av片在线观看秒播厂| 亚洲国产日韩一区二区| av播播在线观看一区| 一级片'在线观看视频| 国产精品熟女久久久久浪| av黄色大香蕉| 国产欧美亚洲国产| 少妇人妻一区二区三区视频| av天堂久久9| 国产精品一二三区在线看| 久久毛片免费看一区二区三区| 天天操日日干夜夜撸| 亚洲欧洲日产国产| 国产淫片久久久久久久久| 国产精品麻豆人妻色哟哟久久| 欧美三级亚洲精品| 国产国拍精品亚洲av在线观看| 亚洲国产色片| 日韩一区二区视频免费看| 亚洲精品乱码久久久久久按摩| 伦理电影免费视频| 国产成人aa在线观看| 80岁老熟妇乱子伦牲交| 久久99一区二区三区| 精品久久国产蜜桃| 伦理电影免费视频| 久久精品国产自在天天线| 丝瓜视频免费看黄片| 国产免费又黄又爽又色| 91精品伊人久久大香线蕉| 大香蕉97超碰在线| 涩涩av久久男人的天堂| 岛国毛片在线播放| 草草在线视频免费看| 亚洲av免费高清在线观看| 99久久精品一区二区三区| 国产亚洲5aaaaa淫片| 大片免费播放器 马上看| 成人综合一区亚洲| 中文天堂在线官网| 一级毛片久久久久久久久女| 免费不卡的大黄色大毛片视频在线观看| 男女国产视频网站| 久久精品久久久久久久性| 最近中文字幕高清免费大全6| 精品一区在线观看国产| 国产乱来视频区| 免费看光身美女| 国产欧美日韩综合在线一区二区 | 狂野欧美激情性bbbbbb| 高清视频免费观看一区二区| 日韩精品有码人妻一区| 久久国内精品自在自线图片| 中文资源天堂在线| 国产伦精品一区二区三区视频9| 春色校园在线视频观看| 99久久精品一区二区三区| 欧美bdsm另类| 中文字幕精品免费在线观看视频 | 国产亚洲91精品色在线| 国产无遮挡羞羞视频在线观看| 亚洲av二区三区四区| 一个人看视频在线观看www免费| 成人黄色视频免费在线看| 最近中文字幕高清免费大全6| 涩涩av久久男人的天堂| 交换朋友夫妻互换小说| 久久精品久久精品一区二区三区| 熟妇人妻不卡中文字幕| 色网站视频免费| 欧美另类一区| 少妇熟女欧美另类| av视频免费观看在线观看| 精品一区在线观看国产| 久久精品国产a三级三级三级| 免费人妻精品一区二区三区视频| 国产精品久久久久久av不卡| 日本欧美视频一区| 美女主播在线视频| 亚洲精品第二区| 国产免费福利视频在线观看| 日韩成人av中文字幕在线观看| 在线观看一区二区三区激情| 亚洲真实伦在线观看| 老司机影院成人| 插逼视频在线观看| 最近中文字幕高清免费大全6| 国产欧美日韩精品一区二区| 久久久久精品久久久久真实原创| 亚洲精品第二区| 国产成人精品无人区| 国产成人免费无遮挡视频| 国产伦精品一区二区三区四那| 亚洲综合色惰| 免费看av在线观看网站| 高清黄色对白视频在线免费看 | 精品亚洲成国产av| 在线观看免费日韩欧美大片 | 欧美 亚洲 国产 日韩一| 精品人妻熟女毛片av久久网站| 男女边吃奶边做爰视频| av在线老鸭窝| 日本色播在线视频| 日韩一区二区三区影片| 久久久久久久久久人人人人人人| 成人影院久久| 国产日韩欧美在线精品| 极品人妻少妇av视频| 天堂俺去俺来也www色官网| 日韩av在线免费看完整版不卡| 丰满少妇做爰视频| 成年人免费黄色播放视频 | 欧美 亚洲 国产 日韩一| 国产一区二区三区综合在线观看 | 日韩一区二区视频免费看| 一本色道久久久久久精品综合| 亚洲国产欧美日韩在线播放 | 免费av中文字幕在线| 国产乱人偷精品视频| 国产精品无大码| 99热全是精品| 精品人妻偷拍中文字幕| 女性被躁到高潮视频| 免费观看在线日韩| 成人无遮挡网站| 一个人看视频在线观看www免费| 欧美另类一区| 精品视频人人做人人爽| 国产欧美另类精品又又久久亚洲欧美| 99热这里只有精品一区| 中文欧美无线码| 久久久久国产精品人妻一区二区| 天堂中文最新版在线下载| 午夜福利视频精品| 欧美+日韩+精品| 中文精品一卡2卡3卡4更新| 中文字幕人妻熟人妻熟丝袜美| 久久国产乱子免费精品| 青春草国产在线视频| 一级,二级,三级黄色视频| 色婷婷久久久亚洲欧美| 九九在线视频观看精品| 亚洲色图综合在线观看| 日韩在线高清观看一区二区三区| 久久人人爽人人爽人人片va| 精品久久久噜噜| 桃花免费在线播放| 妹子高潮喷水视频| 91午夜精品亚洲一区二区三区| av女优亚洲男人天堂| 人妻夜夜爽99麻豆av| 亚洲四区av| 三级经典国产精品| 综合色丁香网| 久久av网站| 成人影院久久| 两个人免费观看高清视频 | 国产无遮挡羞羞视频在线观看| 日韩av不卡免费在线播放| 精品人妻熟女av久视频| 亚洲第一区二区三区不卡| 啦啦啦中文免费视频观看日本| 亚洲人成网站在线观看播放| 日韩强制内射视频| 蜜桃在线观看..| 久久久久人妻精品一区果冻| 人人妻人人看人人澡| 成人漫画全彩无遮挡| av视频免费观看在线观看| 成人特级av手机在线观看| 人妻 亚洲 视频| av一本久久久久| 高清视频免费观看一区二区| 久久97久久精品| 国产欧美亚洲国产| 免费观看无遮挡的男女| 亚洲经典国产精华液单| 国产高清有码在线观看视频| 亚洲国产日韩一区二区| 一二三四中文在线观看免费高清| 成人免费观看视频高清| 日韩欧美精品免费久久| 国产av精品麻豆| 夜夜骑夜夜射夜夜干| 亚洲精品aⅴ在线观看| 国产精品久久久久久精品电影小说| 寂寞人妻少妇视频99o| 最近手机中文字幕大全| 制服丝袜香蕉在线|