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

    自由分子區(qū)內(nèi)納米顆粒的熱泳力計(jì)算*

    2021-03-11 02:39:52崔杰蘇俊杰王軍夏國(guó)棟李志剛
    物理學(xué)報(bào) 2021年5期
    關(guān)鍵詞:剛體勢(shì)能動(dòng)力學(xué)

    崔杰 蘇俊杰 王軍? 夏國(guó)棟 李志剛

    1) (北京工業(yè)大學(xué)能源與動(dòng)力工程學(xué)院, 傳熱強(qiáng)化與過(guò)程節(jié)能教育部重點(diǎn)實(shí)驗(yàn)室, 傳熱與能源利用北京市重點(diǎn)實(shí)驗(yàn)室, 北京 100124)

    2) (香港科技大學(xué)機(jī)械及航空航天工程學(xué)系, 香港)

    基于非平衡態(tài)分子動(dòng)力學(xué)模擬方法, 研究了自由分子區(qū)內(nèi)納米顆粒的熱泳特性.理論研究表明, 納米顆粒與周圍氣體分子之間的非剛體碰撞效應(yīng)會(huì)明顯地改變其熱泳特性, 經(jīng)典的Waldmann 熱泳理論并不適用,但尚未有定量的直接驗(yàn)證.模擬計(jì)算結(jié)果表明: 對(duì)于納米顆粒而言, 當(dāng)氣?固相互作用勢(shì)能較弱或氣體溫度較高時(shí), 氣體分子與納米顆粒之間的非剛體碰撞效應(yīng)可以忽略, Waldmann 熱泳理論與分子動(dòng)力學(xué)模擬結(jié)果吻合較好; 當(dāng)氣?固相互作用勢(shì)能較強(qiáng)或氣體溫度較低時(shí), 非剛體碰撞效應(yīng)較為明顯, Waldmann 熱泳理論與模擬結(jié)果存在較大誤差.基于分子動(dòng)力學(xué)模擬結(jié)果, 對(duì)納米顆粒的等效粒徑進(jìn)行了修正, 并考慮了氣體分子與納米顆粒之間的非剛體碰撞效應(yīng), 理論計(jì)算結(jié)果與分子動(dòng)力學(xué)模擬結(jié)果吻合較好.

    1 引 言

    微細(xì)顆粒在溫度不均勻的氣體中運(yùn)動(dòng)時(shí), 會(huì)受到周圍氣體分子的差異性碰撞, 從而受到一個(gè)沿溫度梯度相反方向上力的作用, 即為熱泳力.顆粒在熱泳力作用下的運(yùn)動(dòng), 稱為熱泳[1?4].通常情況下,來(lái)自高溫區(qū)域的氣體分子與顆粒之間的碰撞更為劇烈, 這便導(dǎo)致顆粒受到的熱泳力方向與溫度梯度方向相反, 指向冷端.在某些情況下, 也可能出現(xiàn)熱泳力與溫度梯度方向相同的情況, 即反向熱泳[2].對(duì)熱泳現(xiàn)象的研究在薄膜制備、微納米制造、環(huán)境科學(xué)、氣溶膠等領(lǐng)域中都有重要的意義[5?10].近年來(lái), 針對(duì)顆粒輸運(yùn)現(xiàn)象的理論與實(shí)驗(yàn)研究越來(lái)越多, 但是, 相關(guān)現(xiàn)象的內(nèi)在機(jī)理還有待深入研究和分析[11?18].

    顆粒的輸運(yùn)特性與其Knudsen 數(shù)Kn (Kn =λ/ L )密切相關(guān).其中, λ 為氣體分子的平均自由程, L 表示顆粒的特征長(zhǎng)度.對(duì)于球形顆粒來(lái)說(shuō),特征長(zhǎng)度為顆粒半徑R, 并使用符號(hào)KnR表示.根據(jù)Knudsen 數(shù)的大小, 可以將顆粒的動(dòng)力學(xué)行為粗略地劃分為連續(xù)介質(zhì)區(qū)(continuum regime)、過(guò)渡區(qū)(transition regime)和自由分子區(qū)(free mole?cule regime).在連續(xù)介質(zhì)區(qū)(KnR? 1), 流體分子的自由程很小, 分子間相互作用的頻率很高, 此時(shí), 可以把流體看作連續(xù)介質(zhì)處理, 使用經(jīng)典的Navier?Stokes 方程來(lái)對(duì)顆粒的輸運(yùn)特性進(jìn)行分析計(jì)算.Epstein[19]基于此方法研究了連續(xù)介質(zhì)區(qū)的熱泳現(xiàn)象, 但求解過(guò)程中采用的邊界條件并不合適.Brock[20]采用完備的滑移邊界條件得到了連續(xù)介質(zhì)區(qū)熱泳力的一階近似解.其他學(xué)者們也嘗試開(kāi)展了高階近似解的計(jì)算[21].然而, 此方法也受到了一定的質(zhì)疑[22?24].在自由分子區(qū)(KnR? 1), 由于氣體分子的平均自由程遠(yuǎn)大于顆粒的特征長(zhǎng)度, 被顆粒反射后的分子需飛行很長(zhǎng)的距離才會(huì)與其他分子相碰, 因此可以忽略入射氣體分子與被顆粒反射氣體分子間的相互作用.此時(shí), 可假設(shè)氣體分子的速度分布函數(shù)不會(huì)受到懸浮顆粒的影響.1956 年,Waldmann[25]基于氣體分子與顆粒之間的剛體碰撞模型, 得到了自由分子區(qū)中顆粒在單原子氣體中所受熱泳力的計(jì)算公式:

    式中, mg為氣體分子的質(zhì)量, kB為Boltzmann 常數(shù), T 為溫度, κ 為氣體熱導(dǎo)率, R 為納米顆粒半徑.? T 為氣體環(huán)境的溫度梯度.越來(lái)越多的理論與實(shí)驗(yàn)研究表明, 隨著KnR增大, 顆粒所受的熱泳力值逐漸接近Waldmann 公式的計(jì)算結(jié)果[2], 因此, (1)式也被稱為熱泳力的自由分子極限.在連續(xù)介質(zhì)區(qū)(KnR~ 1), 氣體分子的離散特性剛開(kāi)始體現(xiàn), 氣體分子的速度分布受顆粒運(yùn)動(dòng)的影響較大, 但又不夠達(dá)到可以看作連續(xù)介質(zhì)的程度, 此時(shí),求解Boltzmann 輸運(yùn)方程是十分困難的.有學(xué)者曾嘗試將自由分子區(qū)的理論擴(kuò)展到過(guò)渡區(qū), 但結(jié)果并不理想[26].Talbot 等[27]通過(guò)改變幾個(gè)Brock 公式中的參數(shù), 使Brock 的連續(xù)介質(zhì)區(qū)理論公式在KnR? 1 時(shí)收斂于(1)式.

    通常情況下, 氣體分子的平均自由程為100 nm量級(jí).對(duì)于粒徑較小的納米顆粒而言, 其特征尺寸遠(yuǎn)小于氣體分子的平均自由程, 因此納米顆粒在氣體中的動(dòng)力學(xué)行為一般處于自由分子區(qū).Waldmann公式是基于氣體分子與顆粒之間的剛體碰撞模型得到的, 但是, 隨著顆粒尺寸減小至納米尺度, 顆粒與氣體分子之間的非剛體碰撞效應(yīng)將成為影響顆粒與氣體分子之間動(dòng)量傳遞的重要影響因素之一[16].文獻(xiàn)[17]中, Li 和Wang 基于非剛體碰撞模型, 得到了自由分子區(qū)內(nèi)納米顆粒所受熱泳力的理論計(jì)算公式:

    式中, mp為顆粒的質(zhì)量; mr= mgmp/(mg+ mp),為氣體分子與顆粒的約化質(zhì)量; Ω(1,1)*與Ω(1,2)*為無(wú)量綱碰撞積分[11].對(duì)于剛體碰撞來(lái)說(shuō), Ω(1,1)*=Ω(1,2)*= 1, (2)式退化為(1)式.換言之, (1)式為(2)式在剛體碰撞假設(shè)下的特殊情形.但是, (2)式的準(zhǔn)確性尚未在實(shí)驗(yàn)或者模擬中得到定量的驗(yàn)證.一方面受限于納米尺度下測(cè)量技術(shù)的不足, 另一方面顆粒的熱泳運(yùn)動(dòng)一般較弱, 很容易被其布朗運(yùn)動(dòng)所掩蓋, 因此關(guān)于納米顆粒熱泳現(xiàn)象的實(shí)驗(yàn)研究往往比較困難[28].分子動(dòng)力學(xué)模擬方法可以直觀地研究并分析顆粒的輸運(yùn)及受力特性[28?33], 在模擬中, 可以在納米顆粒上人為地施加簡(jiǎn)諧勢(shì)來(lái)對(duì)納米顆粒進(jìn)行約束, 從而明顯降低了顆粒布朗運(yùn)動(dòng)對(duì)其熱泳力計(jì)算的影響.

    對(duì)于實(shí)際納米顆粒而言, 其形狀多為非球形(例如碳納米管等).非球形顆粒在流場(chǎng)中的取向會(huì)對(duì)其受力產(chǎn)生重要的影響, 顆粒的瞬時(shí)熱泳力的大小和方向都與其形狀和取向有關(guān), 較為復(fù)雜[34].在自由分子區(qū), 氣體分子的平均自由程遠(yuǎn)大于顆粒粒徑, 因此, 氣體分子與顆粒的隨機(jī)碰撞會(huì)使得顆粒高速旋轉(zhuǎn).顆粒的宏觀受力體現(xiàn)為氣體分子與顆粒之間微觀動(dòng)量傳遞量在長(zhǎng)時(shí)間內(nèi)的統(tǒng)計(jì)平均結(jié)果.如果顆粒沒(méi)有受到較強(qiáng)外勢(shì)場(chǎng)的作用, 顆粒取向分布較為均勻, 可以采用球形顆粒模型近似研究非球形顆粒的熱泳特性.

    本文基于非平衡態(tài)分子動(dòng)力學(xué)模擬方法, 研究了自由分子區(qū)內(nèi)球形納米顆粒的熱泳現(xiàn)象, 計(jì)算得到了納米顆粒所受熱泳力的大小和方向.計(jì)算結(jié)果表明: 當(dāng)氣?固相互作用勢(shì)能較強(qiáng)或氣體溫度較低時(shí), 納米顆粒熱泳力的計(jì)算需要考慮氣體分子與納米顆粒之間的非剛體碰撞效應(yīng); 進(jìn)行顆粒粒徑修正后, 基于Li?Wang 公式(2)的計(jì)算結(jié)果與分子動(dòng)力模擬結(jié)果吻合較好.

    2 分子動(dòng)力學(xué)建模

    2.1 分子動(dòng)力學(xué)模擬系統(tǒng)

    本文采用分子動(dòng)力學(xué)模擬方法計(jì)算納米顆粒在自由分子區(qū)內(nèi)所受熱泳力.分子動(dòng)力學(xué)模擬系統(tǒng)如圖1 所示, 直徑和質(zhì)量分別為DP和mP的納米顆粒被浸入到氣體分子數(shù)為N 的模擬區(qū)域中, 選擇固體分子與氣體分子質(zhì)量相同, 均為m.在系統(tǒng)y 和z 方向采用周期性邊界條件.在x 方向的邊界上采用鏡面反射邊界條件, 并在鄰近邊界的兩端分別建立溫度控制區(qū)域(高低溫?zé)嵩?, 其對(duì)應(yīng)的溫度分別為T(mén)h和Tl.溫度控制區(qū)域在x 方向上的尺寸為lT= 0.1lx, 此處, lx為x 方向的系統(tǒng)尺寸.

    圖1 分子動(dòng)力學(xué)模型圖Fig.1.Snapshot for the MD simulation system.

    采用經(jīng)典的Lennard?Jones(L?J)勢(shì)函數(shù)來(lái)模擬氣?氣、氣?固、固?固分子間的相互作用:

    式中, 參數(shù)ε 和σ 分別為L(zhǎng)?J 勢(shì)參數(shù), rij為原子i 和原子j 之間的距離.系統(tǒng)中的原子(固體原子、氣體原子)參數(shù)均采用Ar 原子的參數(shù)[35]: σ =3.405 ?, ε/kB= 114 K.截?cái)喟霃絩c= 2.5σ.研究氣?固結(jié)合強(qiáng)度對(duì)顆粒所受熱泳的影響發(fā)現(xiàn), 氣?固結(jié)合強(qiáng)度εGP可以在一定范圍內(nèi)變化,

    其中kij為可在一定范圍內(nèi)變化的參數(shù).此外, 在固體顆粒內(nèi)部相鄰原子之間額外引入finite exten?sible nonlinear elastic(FENE)勢(shì)函數(shù)以保持納米顆粒為準(zhǔn)球形, 其勢(shì)能表達(dá)式為[30]

    式中, 勢(shì)能參數(shù)設(shè)定為R0= 1.5σ; kFENE= 30ε/σ2,為FENE 勢(shì)能的彈性系數(shù).kFENE的大小決定了納米顆粒的導(dǎo)熱特性, 文獻(xiàn)[30]對(duì)熱泳力隨kFENE變化情況進(jìn)行了研究分析, 結(jié)果表明納米顆粒的熱導(dǎo)率對(duì)其所受熱泳力影響不大.球形納米顆粒的直徑使用下式進(jìn)行計(jì)算[30]:

    式中, xc為納米顆粒的質(zhì)心, NS為組成納米顆粒的固體原子數(shù)量, 研究過(guò)程中選擇納米顆粒的固體原子數(shù)量至少為40, 納米顆粒的質(zhì)量mP= NSm.Galliero 和Volz[30]的研究表明, 當(dāng)mP/m > 10 時(shí),顆粒的熱泳特性將不再受到顆粒質(zhì)量的影響.相比于氣體分子, 納米顆粒的質(zhì)量與尺寸已足夠大, 因此可忽略顆粒質(zhì)量對(duì)熱泳的影響[36,37].

    由于納米顆粒的布朗運(yùn)動(dòng), 顆粒會(huì)在系統(tǒng)中隨機(jī)行走, 并受到一個(gè)瞬時(shí)的曳力(一般情況下曳力的大小比熱泳力高2—3 個(gè)數(shù)量級(jí)), 為熱泳力的計(jì)算帶來(lái)較大的困難.因此, 本文在顆粒質(zhì)心與系統(tǒng)中心引入額外簡(jiǎn)諧回復(fù)勢(shì)能Uhar:

    式中, rc和ro分別為顆粒質(zhì)心與系統(tǒng)中心坐標(biāo),khar為簡(jiǎn)諧回復(fù)勢(shì)能的彈性系數(shù).作用在納米顆粒上的簡(jiǎn)諧回復(fù)力將均勻施加在納米顆粒的每一個(gè)原子上.簡(jiǎn)諧回復(fù)勢(shì)在系統(tǒng)三個(gè)方向上給顆粒所施加的回復(fù)力分別為F = (Fx, Fy, Fz).當(dāng)系統(tǒng)達(dá)到穩(wěn)定后, 可得顆粒所受到的熱泳力為

    在熱泳力的計(jì)算中引入簡(jiǎn)諧勢(shì), 可以一定程度限制顆粒的隨機(jī)運(yùn)動(dòng), 消除顆粒運(yùn)動(dòng)過(guò)程中所受曳力對(duì)其熱泳力的影響.文獻(xiàn)[32]研究了熱泳力隨彈性系數(shù)khar的變化情況, 結(jié)果表明, 只要彈性系數(shù)khar大于0.1ε/σ2, 其大小對(duì)統(tǒng)計(jì)結(jié)果影響就可以忽略不計(jì).本文取khar= 30ε/σ2.

    本文中, 除特別說(shuō)明, 所有參數(shù)將采用無(wú)量綱形式(用上標(biāo)“*”表示), 采用下列各物理量進(jìn)行無(wú)量綱化: 長(zhǎng)度x0= σ, 溫度T0= ε/kB, 時(shí)間t0=(mσ2/ε)0.5, 力F0= ε/σ.例如, 溫度的無(wú)量綱形式為T(mén)*= T/T0.

    根據(jù)氣體分子運(yùn)動(dòng)論, 氣體分子的平均自由程為

    式中, n 為氣體的分子數(shù)密度, d 為氣體分子的有效硬球直徑.根據(jù)參考文獻(xiàn)[38], 溫度T*= 2.63時(shí)氬原子氣體分子的有效硬球直徑約為0.94σ.為了驗(yàn)證模擬系統(tǒng)中氣體的狀態(tài), 從模擬系統(tǒng)中選擇一組氣體的狀態(tài)參數(shù)與氣體的狀態(tài)方程進(jìn)行對(duì)比,

    式中, P*為NVT 系綜的壓力, Bi(T*)是維里系數(shù).取ΔT*=0.92.P*隨n 的變化情況如圖2 所示.維里系數(shù)的取值分別為= 0.3734,= 0.0421.

    圖2 壓力P*隨數(shù)密度n 變化圖Fig.2.The pressure P* versus the gas molecular number density n.

    分子動(dòng)力學(xué)模擬中, 使用Velocity?Verlet 算法對(duì)顆粒的運(yùn)動(dòng)方程進(jìn)行積分, 為了保證系統(tǒng)能量的穩(wěn)定, 取時(shí)間步長(zhǎng)Δt*= 10—3(Δt = 2 fs).初始時(shí)刻, 氣體分子隨機(jī)分布于模擬系統(tǒng)中, 初始速度為溫度T*下的麥克斯韋分布.系統(tǒng)首先在NVT 系綜下進(jìn)行弛豫, 在高低溫?zé)嵩礈囟鹊乃阈g(shù)平均溫度下運(yùn)行直到t*== 103, 達(dá)到穩(wěn)定.然后, 采用速度校正法[39]將高低溫?zé)嵩〉臏囟确謩e控制為和, 直到t*== 1.5 × 104, 從而得到一個(gè)恒定的溫度差ΔT*, 系統(tǒng)的采樣時(shí)間將在t*=到時(shí)間段內(nèi)進(jìn)行(= 8 × 105).本文所給出的任何物理量都是在此時(shí)間范圍進(jìn)行的長(zhǎng)時(shí)間統(tǒng)計(jì)平均得到的.為了提高結(jié)果的準(zhǔn)確性, 在相同的熱力學(xué)宏觀狀態(tài)下, 取不同初始速度分布, 至少8 次獨(dú)立模擬結(jié)果的平均值作為最終結(jié)果.

    3 有限空間效應(yīng)

    Waldmann 熱泳力公式的一個(gè)前提為與顆粒相互作用的氣體分子是不受空間尺寸束縛的(空間體積無(wú)窮大).在這種情況下, 氣體熱物性不會(huì)受到空間邊界的影響.然而, 在模擬研究中, 氣體分子的活動(dòng)范圍總是被具有一定空間尺寸的邊界所束縛, 在自由分子區(qū)中, 氣體分子的平均自由程可能會(huì)與空間尺寸的量級(jí)相當(dāng), 此時(shí)必須考慮有限的體積效應(yīng)對(duì)氣體熱物性的影響[40].取氣體環(huán)境所處的空間尺寸L 作為特征長(zhǎng)度, 此時(shí)努森數(shù)可用KnL表示.在圖1 所示的分子動(dòng)力學(xué)模型中, L 為x 方向高低溫?zé)嵩〉拈g距, 即L = 0.8lx.1972 年,Phillips[41]借助Chapman?Enskog 分布函數(shù)對(duì)Bo?ltzmann 輸運(yùn)方程進(jìn)行了求解, 得到了有限空間內(nèi)顆粒的熱泳力表達(dá)式:

    式中, αmt和αmn分別為顆粒表面切向和法向的動(dòng)量協(xié)調(diào)系數(shù), α1和α2分別為高低溫壁面的溫度協(xié)調(diào)系數(shù).在本文的分子動(dòng)力學(xué)模擬系統(tǒng)中, 可以認(rèn)為氣體分子是完全協(xié)調(diào)的, 即= α2≈ 1 以及αmt= αmn≈ 1.

    分 別 采 用 Phillips (11)式 和 Waldmann(1)式進(jìn)行熱泳力理論計(jì)算, 并與分子動(dòng)力學(xué)模擬計(jì)算結(jié)果進(jìn)行對(duì)比分析.Phillips 公式的計(jì)算結(jié)果表示為FT,Phillips, 分子動(dòng)力學(xué)模擬結(jié)果為FT,MD.Waldmann 公式將分別使用氣體的宏觀熱導(dǎo)率κ 和有效熱導(dǎo)率κ'計(jì)算, 計(jì)算結(jié)果分別用FT,B和FT,E表示.圖3(a)與圖3(b)分別表示在kij= 5.26 (強(qiáng)氣?固相互作用下)和kij= 1.0 (弱氣?固相互作用下)情況下FT/FT,B隨l 的變化圖.由圖3 可以看出, 模擬系統(tǒng)中“壁面”間距是影響顆粒所受熱泳力的重要因素之一, 這說(shuō)明有限空間效應(yīng)確實(shí)存在.FT,E與FT,Phillips吻合較好, 表明Phillips 對(duì)Wald?mann 熱泳力計(jì)算式的修正是可以通過(guò)對(duì)氣體熱導(dǎo)率的修正來(lái)實(shí)現(xiàn)的.

    表1 分子動(dòng)力學(xué)模擬系統(tǒng)的幾何特征參數(shù)Table 1.Geometric and characteristic parameters of the simulation systems.

    圖3 FT, MD/FT, B, FT, E/FT, B 和FT, Phillips/FT, B 隨模擬空間尺寸的變化 (a) kij = 5.26; (b) kij =1Fig.3.Influence of the size of simulation box on FT, MD/FT, B,FT, E/FT, B and FT, Phillips/FT, B: (a) kij =5.26, (b) kij =1.

    對(duì)于強(qiáng)氣?固相互作用(圖3(a)), 分子動(dòng)力學(xué)結(jié)果FT,MD與理論結(jié)果FT,E和FT,Phillips存在較大誤差, 這種偏差是由于氣體分子與納米顆粒之間的非剛體碰撞所引起的.對(duì)于弱氣?固相互作用(圖3(b)), 氣體分子與納米顆粒之間的非剛體碰撞效應(yīng)較弱, 此時(shí), 基于剛體碰撞模型假設(shè)的Wald?mann 理論仍然適用于納米顆粒, FT,E和FT,Phillips都與分子動(dòng)力學(xué)結(jié)果FT,MD吻合較好.圖4 為熱泳力隨氣體有效熱導(dǎo)率κ' 的變化曲線.與前文結(jié)果類似, 對(duì)于弱氣固耦合的情況(kij= 1.0), FT,E和與分子動(dòng)力學(xué)結(jié)果FT,MD吻合較好; 對(duì)于強(qiáng)氣固耦合(kij= 5.26), 分子動(dòng)力學(xué)模擬得到的熱泳力則明顯高于基于有效熱導(dǎo)率的理論計(jì)算結(jié)果.

    圖4 熱泳力FT,MD 和FT,E 隨氣體有效熱導(dǎo)率κ' 的變化Fig.4.Influence of the effective thermal conductivity of the media gas on the thermophoretic forces.

    4 結(jié)果與討論

    4.1 系統(tǒng)溫度(氣-固結(jié)合強(qiáng)度)影響

    由前文所述可知, 將氣體的有效熱導(dǎo)率引入到Waldmann 公式中, 可得到有限空間中氣體的Waldmann 自由分子極限.但是, 對(duì)于納米顆粒來(lái)說(shuō), 特別是當(dāng)氣?固相互作用較強(qiáng)時(shí), Waldmann 熱泳力計(jì)算式的誤差較大.基于分子動(dòng)力學(xué)模擬結(jié)果,熱泳力隨系統(tǒng)平均溫度氣?固結(jié)合強(qiáng)度kij和溫度梯度 ? T 的變化情況如圖5 所示(圖5(a)中的插圖為在相同系統(tǒng)參數(shù)下本文分子動(dòng)力學(xué)結(jié)果與文獻(xiàn)[28]中結(jié)果的對(duì)比).由圖5 可以看出, 對(duì)于納米顆粒來(lái)說(shuō), Waldmann 熱泳力計(jì)算式的適用性受系統(tǒng)平均溫度與氣?固結(jié)合強(qiáng)度kij的影響, 當(dāng)系統(tǒng)溫度較低或氣?固結(jié)合強(qiáng)度較大時(shí), Waldmann公式給出的熱泳力計(jì)算結(jié)果誤差較大; 隨著系統(tǒng)溫度的升高或氣?固結(jié)合強(qiáng)度的減小, 納米顆粒所受熱泳力的分子動(dòng)力學(xué)計(jì)算結(jié)果逐漸收斂于Wald?mann 公式的預(yù)測(cè)值.這是由氣體分子與納米顆粒間非剛體碰撞效應(yīng)的強(qiáng)弱受溫度或氣?固結(jié)合強(qiáng)度影響所導(dǎo)致的.當(dāng)結(jié)合強(qiáng)度較大或溫度較低時(shí), 勢(shì)能在氣體分子與納米顆粒進(jìn)行動(dòng)量交換的運(yùn)動(dòng)軌跡中占主導(dǎo)地位; 而當(dāng)結(jié)合強(qiáng)度較小或溫度較高時(shí), 氣?固相互作用勢(shì)能的影響十分微弱, 剛體碰撞模型近似成立.

    圖5 不同參數(shù)下熱泳力的分子動(dòng)力學(xué)結(jié)果FT,MD 與Waldmann 公式結(jié)果FT,E 比較圖 (a) 環(huán)境溫度T*; (b) 氣?固結(jié)合強(qiáng)度kij; (c) 溫度梯度?T?Fig.5.The variation of thermophoretic force FT between present MD simulation result FT, MD and Waldmann equa?tion result FT, E under different parameters: (a) The environ?ment temperature ; (b) the intensity of gas?particle in?teraction kij; (c) temperature gradient ? T?.

    4.2 顆粒尺寸影響

    顆粒尺寸的大小不僅影響著顆粒在高低溫兩側(cè)的動(dòng)量交換, 還影響著顆粒與氣體分子之間的碰撞模型, 是影響顆粒所受熱泳力的重要因素之一.為了便于比較, 定義無(wú)量綱熱泳力:

    式中, FT為熱泳力的理論計(jì)算結(jié)果, FT= FT,E為基于有效熱導(dǎo)率的Waldmann 理論計(jì)算結(jié)果, FT=FT,LW為考慮非剛體碰撞效應(yīng)之后(2)式的計(jì)算結(jié)果, 對(duì)應(yīng)的無(wú)量綱熱泳力分別表示為和圖6 展示無(wú)量綱熱泳力和隨顆粒直徑的變化圖.模擬過(guò)程中系統(tǒng)的參數(shù)為:= 55,ΔT*= 0.92, n*=0.0054.顆粒粒徑變化范圍為 DP?≈ 2.5 到≈9.98, 且模擬系統(tǒng)處于自由分子區(qū).由圖6(a)可以看出, 隨著顆粒尺寸逐漸接近分子尺寸(= 1),無(wú)量綱熱泳力單調(diào)增加, 而對(duì)于顆粒尺寸較大的納米顆粒, 顆粒所受熱泳力的分子動(dòng)力學(xué)結(jié)果則逐漸收斂于Waldmann 公式的預(yù)測(cè)值(由于計(jì)算能力的限制, 并沒(méi)有給出更大尺寸顆粒的計(jì)算結(jié)果).對(duì)于納米顆粒而言, 由于納米顆粒與氣體分子之間的非剛體碰撞效應(yīng), 熱泳力結(jié)果并不能收斂于經(jīng)典的Waldmann 理論值.(2)式對(duì)Waldmann理論進(jìn)行了修正, 修正后的無(wú)量綱熱泳力隨顆粒尺寸變化結(jié)果如圖6(b)所示.相對(duì)于Wald?mann 理論的計(jì)算結(jié)果, 無(wú)量綱熱泳力計(jì)算誤差明星降低, 這驗(yàn)證了(2)式的準(zhǔn)確性.但是, 不同氣?固結(jié)合強(qiáng)度下的計(jì)算結(jié)果還存在一定的誤差,這是由于當(dāng)氣?固結(jié)合強(qiáng)度較大時(shí), 氣體分子的動(dòng)能不足以克服氣?固結(jié)合勢(shì)能的束縛, 氣體分子會(huì)被吸附于納米顆粒表面, 從而引起納米顆粒尺寸的增加, 造成理論與模擬結(jié)果的偏差.

    圖6 無(wú)量綱熱泳力 隨顆粒直徑 變化圖 (a)剛體碰撞模型; (b) 非剛體碰撞模型FT = FT, LW1; (c) 非剛體碰撞模型FT = FT, LW2(考慮氣體分子吸附引起的顆粒尺寸變化)Fig.6.The dimensionless thermophoretic force as a function of for different gas?particle interaction intensity:(a) The rigid body collision model FT = FT, E; (b) the non?rigid body collision model FT = FT, LW1; (c) FT = FT, LW2(wherein the vary of particle size is considered).

    圖7 對(duì)顆粒表面氣?固相互作用勢(shì)能以及氣體分子動(dòng)能分布進(jìn)行了統(tǒng)計(jì), 同時(shí)計(jì)算了顆粒表面氣體分子密度分布函數(shù)g.圖7 中, 黑色與藍(lán)色虛線分別為顆粒表面第一層與第二層原子的位置, r'=.由圖7 可以看出, 當(dāng)氣?固結(jié)合強(qiáng)度較小時(shí)(圖7(a)), 勢(shì)井深度較淺, 氣體分子所攜帶的動(dòng)能足以克服勢(shì)能的束縛.此時(shí), 納米顆粒表面的氣體分子數(shù)密度與環(huán)境中的氣體分子數(shù)密度相當(dāng), 并沒(méi)有發(fā)生吸附現(xiàn)象.當(dāng)氣?固結(jié)合強(qiáng)度較大時(shí)(圖7(b)), 顆粒表面第一層內(nèi)的勢(shì)井深度較深,而氣體分子的動(dòng)能相對(duì)太低, 使得該層氣體分子無(wú)法擺脫勢(shì)能的束縛, 從而被吸附在納米顆粒表面,這就造成了納米顆粒表面的第一層內(nèi)的氣體分子數(shù)密度遠(yuǎn)遠(yuǎn)高于環(huán)境中的氣體分子數(shù)密度.通過(guò)對(duì)納米顆粒的粒徑進(jìn)行修正, 可以得到修正后的等效粒徑, 將等效粒徑代入 (2)式, 無(wú)量綱后的結(jié)果表示為如圖6(c)為無(wú)量綱熱泳力隨顆粒尺寸的變化結(jié)果.結(jié)果表明, 計(jì)算結(jié)果與分子動(dòng)力學(xué)模擬結(jié)果吻合較好.

    圖7 納米顆粒表面氣體原子的勢(shì)能Ep, 動(dòng)能Ek 及密度函數(shù)g 分布圖Fig.7.The potential energy Ep, kinetic energy Ek and dens?ity function g on the surface of nanoparticles.

    5 結(jié) 論

    本文研究了自由分子區(qū)內(nèi)納米顆粒的熱泳特性, 基于非平衡態(tài)分子動(dòng)力學(xué)模擬方法計(jì)算了納米顆粒所受熱泳力.考慮有限空間效應(yīng)對(duì)氣體熱物性的影響, 引入氣體的有效熱導(dǎo)率計(jì)算得到了較為準(zhǔn)確的Waldmann 自由分子極限, 并與分子動(dòng)力學(xué)模擬結(jié)果進(jìn)行了對(duì)比.結(jié)果表明, 對(duì)于納米顆粒來(lái)說(shuō), 當(dāng)氣?固相互作用勢(shì)能較弱或氣體溫度較高時(shí),Waldmann 熱泳力公式與分子動(dòng)力學(xué)模擬結(jié)果吻合較好; 當(dāng)氣?固相互作用勢(shì)能較強(qiáng)或氣體溫度較低時(shí), 由于氣體分子與納米顆粒之間的非剛體碰撞效應(yīng)較為明顯, 基于剛體碰撞模型假設(shè)的Wald?mann 計(jì)算式與模擬結(jié)果存在較大誤差.基于Li?Wang 計(jì)算式, 并考慮由于納米顆粒對(duì)氣體分子的吸附效應(yīng)引起顆粒等效尺寸的變化, 計(jì)算結(jié)果與分子動(dòng)力學(xué)模擬結(jié)果較為吻合.

    猜你喜歡
    剛體勢(shì)能動(dòng)力學(xué)
    “動(dòng)能和勢(shì)能”知識(shí)鞏固
    作 品:景觀設(shè)計(jì)
    ——《勢(shì)能》
    文化縱橫(2022年3期)2022-09-07 11:43:18
    《空氣動(dòng)力學(xué)學(xué)報(bào)》征稿簡(jiǎn)則
    “動(dòng)能和勢(shì)能”知識(shí)鞏固
    “動(dòng)能和勢(shì)能”隨堂練
    差值法巧求剛體轉(zhuǎn)動(dòng)慣量
    車載冷發(fā)射系統(tǒng)多剛體動(dòng)力學(xué)快速仿真研究
    基于隨機(jī)-動(dòng)力學(xué)模型的非均勻推移質(zhì)擴(kuò)散
    剛體定點(diǎn)轉(zhuǎn)動(dòng)的瞬軸、極面動(dòng)態(tài)演示教具
    TNAE的合成和熱分解動(dòng)力學(xué)
    久久精品国产99精品国产亚洲性色| 国产精品日韩av在线免费观看| 亚洲真实伦在线观看| 身体一侧抽搐| 精品熟女少妇av免费看| 日日摸夜夜添夜夜添小说| 99热6这里只有精品| 久久久久国内视频| 国产精品爽爽va在线观看网站| 久久精品国产亚洲av香蕉五月| 日韩一本色道免费dvd| 国产精品野战在线观看| 一级av片app| 日韩强制内射视频| 久久久久久久久大av| 欧美日本亚洲视频在线播放| 亚州av有码| 色综合色国产| 看黄色毛片网站| 久久精品夜色国产| 久久天躁狠狠躁夜夜2o2o| 精品免费久久久久久久清纯| 日韩中字成人| 亚洲一区二区三区色噜噜| 久久久久久久久中文| 人妻丰满熟妇av一区二区三区| 国产高清三级在线| 又爽又黄无遮挡网站| 好男人在线观看高清免费视频| 岛国在线免费视频观看| 欧美日韩在线观看h| 欧美最新免费一区二区三区| 女人被狂操c到高潮| 最近手机中文字幕大全| 久久韩国三级中文字幕| 3wmmmm亚洲av在线观看| 亚洲乱码一区二区免费版| 高清午夜精品一区二区三区 | 少妇猛男粗大的猛烈进出视频 | 最近中文字幕高清免费大全6| 午夜福利在线在线| 亚洲欧美精品自产自拍| 伊人久久精品亚洲午夜| 国产精品人妻久久久久久| 在线播放国产精品三级| 美女cb高潮喷水在线观看| 校园人妻丝袜中文字幕| 少妇的逼好多水| 色吧在线观看| 国产精品av视频在线免费观看| 高清午夜精品一区二区三区 | 精品无人区乱码1区二区| 69人妻影院| 亚洲性夜色夜夜综合| 欧美三级亚洲精品| 国产精品1区2区在线观看.| 伊人久久精品亚洲午夜| 成人一区二区视频在线观看| 51国产日韩欧美| 日韩亚洲欧美综合| 插阴视频在线观看视频| 国产成人a∨麻豆精品| 国产精品三级大全| 日韩三级伦理在线观看| 深夜a级毛片| av在线观看视频网站免费| 亚洲欧美精品自产自拍| 老熟妇乱子伦视频在线观看| 国产伦精品一区二区三区四那| 99热这里只有精品一区| 亚洲一区高清亚洲精品| 人妻少妇偷人精品九色| 97在线视频观看| 女生性感内裤真人,穿戴方法视频| 热99在线观看视频| 亚洲精品日韩av片在线观看| 中文字幕av成人在线电影| av天堂在线播放| 日产精品乱码卡一卡2卡三| 国产精品一区二区性色av| av.在线天堂| 大香蕉久久网| 白带黄色成豆腐渣| 日本熟妇午夜| a级一级毛片免费在线观看| 久久精品国产鲁丝片午夜精品| 亚洲美女搞黄在线观看 | 色5月婷婷丁香| 亚洲无线观看免费| 在线播放无遮挡| 国产激情偷乱视频一区二区| 色av中文字幕| 韩国av在线不卡| 最近2019中文字幕mv第一页| 一个人观看的视频www高清免费观看| 观看美女的网站| a级毛片a级免费在线| 国产亚洲91精品色在线| 中文亚洲av片在线观看爽| 波多野结衣巨乳人妻| 大又大粗又爽又黄少妇毛片口| 久久草成人影院| 欧美高清性xxxxhd video| 99久久成人亚洲精品观看| 精品久久久久久久久久久久久| 在线播放国产精品三级| 精品久久久久久成人av| 亚洲图色成人| 国内久久婷婷六月综合欲色啪| 亚洲第一电影网av| 欧美中文日本在线观看视频| 日韩高清综合在线| 狠狠狠狠99中文字幕| 观看美女的网站| 俄罗斯特黄特色一大片| 长腿黑丝高跟| 亚洲av成人av| АⅤ资源中文在线天堂| 久久这里只有精品中国| 久久99热这里只有精品18| 五月伊人婷婷丁香| www日本黄色视频网| 亚洲精品一区av在线观看| 日本免费一区二区三区高清不卡| 午夜福利高清视频| 国产高潮美女av| 成年女人看的毛片在线观看| 精品福利观看| 日韩成人伦理影院| 亚洲精华国产精华液的使用体验 | 午夜福利在线观看吧| 免费大片18禁| 亚州av有码| 精品久久久久久久久av| 校园人妻丝袜中文字幕| 亚洲国产精品成人综合色| 国产av不卡久久| 国产精品女同一区二区软件| 久久久精品94久久精品| 欧美日韩乱码在线| 欧美激情在线99| 99久久无色码亚洲精品果冻| 亚洲欧美日韩卡通动漫| 成年女人永久免费观看视频| 精品久久国产蜜桃| 国产黄a三级三级三级人| 国产一区二区激情短视频| 亚洲无线在线观看| 精品免费久久久久久久清纯| 男女之事视频高清在线观看| 久久韩国三级中文字幕| 国产精华一区二区三区| 在线看三级毛片| 欧洲精品卡2卡3卡4卡5卡区| 天堂网av新在线| av在线亚洲专区| 婷婷精品国产亚洲av| 人妻夜夜爽99麻豆av| 亚洲成人精品中文字幕电影| 一区福利在线观看| 99热只有精品国产| 国产精品女同一区二区软件| 91久久精品国产一区二区成人| av女优亚洲男人天堂| 蜜桃久久精品国产亚洲av| 亚洲在线观看片| 久久亚洲国产成人精品v| 中文亚洲av片在线观看爽| 色在线成人网| www.色视频.com| 97超视频在线观看视频| 久久久精品94久久精品| 亚洲精品影视一区二区三区av| 国产欧美日韩一区二区精品| 美女高潮的动态| 直男gayav资源| 精品国内亚洲2022精品成人| 看黄色毛片网站| 成人高潮视频无遮挡免费网站| 欧美最新免费一区二区三区| 久久精品人妻少妇| 久久亚洲精品不卡| 亚洲熟妇中文字幕五十中出| 69人妻影院| 日本欧美国产在线视频| 一级av片app| 国产色爽女视频免费观看| av在线亚洲专区| 午夜福利在线观看免费完整高清在 | 丰满的人妻完整版| 久久久精品欧美日韩精品| 欧美xxxx黑人xx丫x性爽| 色视频www国产| 国产精品伦人一区二区| 国产久久久一区二区三区| 成人美女网站在线观看视频| 少妇高潮的动态图| 成人亚洲欧美一区二区av| 自拍偷自拍亚洲精品老妇| 亚洲国产日韩欧美精品在线观看| a级一级毛片免费在线观看| 一个人看视频在线观看www免费| 黄色日韩在线| 深夜a级毛片| 少妇的逼好多水| 欧美日韩乱码在线| 99久久精品国产国产毛片| 1024手机看黄色片| 男人狂女人下面高潮的视频| 成人综合一区亚洲| 不卡视频在线观看欧美| 寂寞人妻少妇视频99o| 亚洲自拍偷在线| 亚洲av.av天堂| 亚洲国产欧美人成| 日本a在线网址| 国产成人aa在线观看| 欧美极品一区二区三区四区| 精品久久久久久久人妻蜜臀av| 免费无遮挡裸体视频| av中文乱码字幕在线| 黑人高潮一二区| 亚洲性夜色夜夜综合| 精品久久久久久久末码| 国产黄色视频一区二区在线观看 | 久久久久久久久中文| 亚洲人成网站在线播放欧美日韩| 色综合亚洲欧美另类图片| а√天堂www在线а√下载| av女优亚洲男人天堂| 99久国产av精品| 丝袜喷水一区| 色5月婷婷丁香| 插逼视频在线观看| 97碰自拍视频| 亚洲18禁久久av| 深爱激情五月婷婷| 俺也久久电影网| 免费看日本二区| 在线国产一区二区在线| 少妇高潮的动态图| 国内少妇人妻偷人精品xxx网站| 丰满人妻一区二区三区视频av| 欧美一区二区国产精品久久精品| 岛国在线免费视频观看| 三级毛片av免费| 最近最新中文字幕大全电影3| 国产伦在线观看视频一区| 日本黄大片高清| 一进一出好大好爽视频| 嫩草影院精品99| 亚洲人与动物交配视频| 国产高清三级在线| 亚洲国产精品合色在线| 国内少妇人妻偷人精品xxx网站| 3wmmmm亚洲av在线观看| 淫秽高清视频在线观看| 午夜免费激情av| 91精品国产九色| 欧美日本视频| 日产精品乱码卡一卡2卡三| 免费看美女性在线毛片视频| 国产淫片久久久久久久久| 日韩欧美国产在线观看| 日韩在线高清观看一区二区三区| 欧美区成人在线视频| 一a级毛片在线观看| 舔av片在线| 最近2019中文字幕mv第一页| 免费人成在线观看视频色| 欧美最新免费一区二区三区| 免费在线观看影片大全网站| 18禁黄网站禁片免费观看直播| 成年女人永久免费观看视频| 成人无遮挡网站| 亚洲成a人片在线一区二区| 国产欧美日韩精品一区二区| 国产老妇女一区| 亚洲欧美成人综合另类久久久 | 日韩 亚洲 欧美在线| 久久久久精品国产欧美久久久| 国产精品乱码一区二三区的特点| 日韩欧美 国产精品| 免费av毛片视频| 综合色丁香网| 少妇熟女欧美另类| 国产 一区 欧美 日韩| 在线观看免费视频日本深夜| 亚洲精品一区av在线观看| 97超视频在线观看视频| 丝袜美腿在线中文| 干丝袜人妻中文字幕| 国产欧美日韩精品一区二区| 久久久久久九九精品二区国产| 免费不卡的大黄色大毛片视频在线观看 | 亚洲三级黄色毛片| 少妇裸体淫交视频免费看高清| 精品无人区乱码1区二区| 真人做人爱边吃奶动态| 国产av不卡久久| 亚洲aⅴ乱码一区二区在线播放| 一级黄片播放器| 色5月婷婷丁香| av国产免费在线观看| 国产精品女同一区二区软件| 免费看日本二区| 亚洲,欧美,日韩| 色在线成人网| 久久久久久伊人网av| av在线播放精品| 国产精品伦人一区二区| 亚洲人成网站在线观看播放| 中文字幕av成人在线电影| 天堂动漫精品| 国产麻豆成人av免费视频| 国产欧美日韩精品一区二区| 日本黄大片高清| 99热网站在线观看| 久久韩国三级中文字幕| 中文字幕免费在线视频6| 少妇熟女aⅴ在线视频| 少妇人妻一区二区三区视频| 成人午夜高清在线视频| 亚洲国产精品成人久久小说 | 成年女人永久免费观看视频| 亚洲av免费高清在线观看| 国产aⅴ精品一区二区三区波| 欧美激情久久久久久爽电影| 人人妻,人人澡人人爽秒播| 国产精品三级大全| 亚洲欧美日韩高清专用| 日韩大尺度精品在线看网址| av福利片在线观看| 亚洲av美国av| 99久久精品国产国产毛片| 色av中文字幕| 又黄又爽又免费观看的视频| 黄片wwwwww| 69av精品久久久久久| 亚洲美女搞黄在线观看 | 九九久久精品国产亚洲av麻豆| 校园人妻丝袜中文字幕| 在线看三级毛片| 色5月婷婷丁香| 日本色播在线视频| 亚洲aⅴ乱码一区二区在线播放| 精品午夜福利视频在线观看一区| 国产精品乱码一区二三区的特点| 国产高清激情床上av| 亚洲aⅴ乱码一区二区在线播放| 一卡2卡三卡四卡精品乱码亚洲| 亚洲国产精品合色在线| 午夜福利18| 免费人成视频x8x8入口观看| 久久精品国产亚洲av天美| 白带黄色成豆腐渣| 午夜福利18| 综合色丁香网| 免费看日本二区| 亚洲经典国产精华液单| 亚洲美女视频黄频| 日本一二三区视频观看| 中文字幕人妻熟人妻熟丝袜美| 久久久久国产精品人妻aⅴ院| 日韩中字成人| 一区二区三区四区激情视频 | 亚洲内射少妇av| 久久久色成人| 小蜜桃在线观看免费完整版高清| 12—13女人毛片做爰片一| 国产成人福利小说| 免费在线观看成人毛片| 成人一区二区视频在线观看| 久久热精品热| 午夜亚洲福利在线播放| 欧美成人免费av一区二区三区| 伦理电影大哥的女人| 在线观看一区二区三区| av视频在线观看入口| 欧洲精品卡2卡3卡4卡5卡区| 免费无遮挡裸体视频| 91麻豆精品激情在线观看国产| 欧美xxxx性猛交bbbb| 国产伦一二天堂av在线观看| 欧美日韩在线观看h| 免费大片18禁| 精品99又大又爽又粗少妇毛片| .国产精品久久| 久久婷婷人人爽人人干人人爱| 久久午夜亚洲精品久久| 国产精品野战在线观看| 日本与韩国留学比较| 我要搜黄色片| 精品99又大又爽又粗少妇毛片| 女同久久另类99精品国产91| 香蕉av资源在线| 国产麻豆成人av免费视频| 少妇熟女aⅴ在线视频| 久久久成人免费电影| 老司机影院成人| 国产在线男女| 亚洲成人精品中文字幕电影| 亚洲精品久久国产高清桃花| 国产aⅴ精品一区二区三区波| 人妻夜夜爽99麻豆av| 露出奶头的视频| 成熟少妇高潮喷水视频| 亚洲av二区三区四区| 99热网站在线观看| 成人性生交大片免费视频hd| 日韩欧美在线乱码| 欧美bdsm另类| 国产成人影院久久av| 大又大粗又爽又黄少妇毛片口| 波多野结衣高清无吗| 在线观看美女被高潮喷水网站| 成人亚洲精品av一区二区| 欧美最新免费一区二区三区| 一进一出好大好爽视频| 在现免费观看毛片| 日韩在线高清观看一区二区三区| 男人舔女人下体高潮全视频| videossex国产| 成人欧美大片| 亚洲乱码一区二区免费版| 一个人免费在线观看电影| 欧洲精品卡2卡3卡4卡5卡区| 两个人的视频大全免费| 亚洲国产欧洲综合997久久,| 中出人妻视频一区二区| 国产熟女欧美一区二区| 99在线人妻在线中文字幕| 中文字幕人妻熟人妻熟丝袜美| 精品不卡国产一区二区三区| 国产国拍精品亚洲av在线观看| 久久精品国产亚洲网站| 亚洲最大成人中文| 97超视频在线观看视频| 美女高潮的动态| 内射极品少妇av片p| 免费电影在线观看免费观看| 麻豆av噜噜一区二区三区| 最近2019中文字幕mv第一页| 久久久国产成人免费| videossex国产| 一级a爱片免费观看的视频| 麻豆国产av国片精品| 久久久a久久爽久久v久久| 最近最新中文字幕大全电影3| 欧美三级亚洲精品| 亚洲成人久久性| 久久精品人妻少妇| 亚洲一级一片aⅴ在线观看| 久久综合国产亚洲精品| 欧美极品一区二区三区四区| 伦精品一区二区三区| 黄片wwwwww| 午夜爱爱视频在线播放| 国产探花在线观看一区二区| 日本-黄色视频高清免费观看| 观看美女的网站| 干丝袜人妻中文字幕| 国产精品乱码一区二三区的特点| 国产欧美日韩精品亚洲av| 99在线视频只有这里精品首页| 亚洲三级黄色毛片| 91av网一区二区| 日日干狠狠操夜夜爽| 国产探花在线观看一区二区| 露出奶头的视频| 可以在线观看毛片的网站| 舔av片在线| 人妻丰满熟妇av一区二区三区| 国产一区二区激情短视频| 麻豆一二三区av精品| 十八禁网站免费在线| 一级毛片aaaaaa免费看小| 在线国产一区二区在线| 亚洲精品日韩在线中文字幕 | 99国产极品粉嫩在线观看| 天堂动漫精品| 欧美另类亚洲清纯唯美| 在线观看免费视频日本深夜| 两性午夜刺激爽爽歪歪视频在线观看| 亚洲第一区二区三区不卡| 狂野欧美激情性xxxx在线观看| 中文字幕人妻熟人妻熟丝袜美| 女同久久另类99精品国产91| 精品午夜福利在线看| 极品教师在线视频| 色综合色国产| 美女免费视频网站| 丝袜美腿在线中文| 亚洲天堂国产精品一区在线| 欧美日韩综合久久久久久| 中文亚洲av片在线观看爽| 亚洲无线在线观看| 99热6这里只有精品| 亚洲中文日韩欧美视频| 国产成人精品久久久久久| 一区福利在线观看| 成人漫画全彩无遮挡| 亚洲图色成人| av天堂中文字幕网| 亚洲熟妇熟女久久| 国产精品一二三区在线看| 日日摸夜夜添夜夜添小说| 欧美不卡视频在线免费观看| 日本色播在线视频| 校园人妻丝袜中文字幕| 中文在线观看免费www的网站| 一区二区三区高清视频在线| 搡老妇女老女人老熟妇| 亚洲在线自拍视频| 国产精品久久久久久精品电影| 99久国产av精品| 亚洲人成网站在线播放欧美日韩| 国产三级在线视频| 两个人视频免费观看高清| 午夜精品国产一区二区电影 | 国产成年人精品一区二区| 午夜激情福利司机影院| 国内精品久久久久精免费| 亚洲精品亚洲一区二区| 午夜福利18| 最近中文字幕高清免费大全6| aaaaa片日本免费| 久久鲁丝午夜福利片| 日韩中字成人| 天堂√8在线中文| 简卡轻食公司| 成年av动漫网址| 伦理电影大哥的女人| 三级国产精品欧美在线观看| 精品人妻视频免费看| 精品久久久久久久久av| 美女高潮的动态| 少妇高潮的动态图| 国产高清激情床上av| 人妻少妇偷人精品九色| 黄色日韩在线| 淫秽高清视频在线观看| 如何舔出高潮| 成年女人看的毛片在线观看| 乱人视频在线观看| 91精品国产九色| 成人特级黄色片久久久久久久| 免费大片18禁| 成年av动漫网址| 中文亚洲av片在线观看爽| 国产一区亚洲一区在线观看| 欧美中文日本在线观看视频| 在线天堂最新版资源| 超碰av人人做人人爽久久| 日韩一区二区视频免费看| 亚洲最大成人av| 女的被弄到高潮叫床怎么办| 十八禁网站免费在线| 看片在线看免费视频| 毛片一级片免费看久久久久| 少妇高潮的动态图| 成人特级av手机在线观看| 亚洲,欧美,日韩| a级毛片a级免费在线| 内射极品少妇av片p| 女人十人毛片免费观看3o分钟| 国产极品精品免费视频能看的| 国模一区二区三区四区视频| av天堂中文字幕网| 男插女下体视频免费在线播放| 亚洲第一区二区三区不卡| 天美传媒精品一区二区| 亚洲丝袜综合中文字幕| 成年女人看的毛片在线观看| 男女视频在线观看网站免费| 亚洲丝袜综合中文字幕| 99九九线精品视频在线观看视频| 中文字幕人妻熟人妻熟丝袜美| 99久久久亚洲精品蜜臀av| 午夜激情福利司机影院| 香蕉av资源在线| 99精品在免费线老司机午夜| 亚洲七黄色美女视频| 黄色日韩在线| 麻豆精品久久久久久蜜桃| 国产不卡一卡二| 中文在线观看免费www的网站| 久久中文看片网| 成年女人毛片免费观看观看9| 久久亚洲精品不卡| 亚洲久久久久久中文字幕| 亚洲精品粉嫩美女一区| 中文在线观看免费www的网站| 有码 亚洲区| 日韩精品中文字幕看吧| 少妇人妻精品综合一区二区 | 亚洲精品日韩av片在线观看| 一本精品99久久精品77| 我要搜黄色片| 亚洲成人久久性| 99久久成人亚洲精品观看| 成年女人看的毛片在线观看| 校园春色视频在线观看| 露出奶头的视频| 不卡一级毛片| 中文字幕免费在线视频6| a级毛片免费高清观看在线播放| 国产成人一区二区在线| 免费大片18禁| 日韩,欧美,国产一区二区三区 | 伦精品一区二区三区| 18禁在线无遮挡免费观看视频 | 亚洲精品影视一区二区三区av| 国产精品99久久久久久久久| 亚洲av一区综合|