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

    分子動力學(xué)模擬剪切速度對納米間隙中角鯊?fù)榻缑婊频挠绊?/h1>
    2016-07-04 03:43:57潘伶高誠輝福州大學(xué)機(jī)械工程及自動化學(xué)院福建福州350108
    化工學(xué)報 2016年4期
    關(guān)鍵詞:聚合物

    潘伶,高誠輝(福州大學(xué)機(jī)械工程及自動化學(xué)院,福建 福州 350108)

    ?

    分子動力學(xué)模擬剪切速度對納米間隙中角鯊?fù)榻缑婊频挠绊?/p>

    潘伶,高誠輝
    (福州大學(xué)機(jī)械工程及自動化學(xué)院,福建 福州 350108)

    摘要:采用聚合物一致性力場(PCFF),分別在7種剪切速度V和3種油膜厚度h下對納米間隙中潤滑劑角鯊?fù)檫M(jìn)行分子動力學(xué)模擬,分析固液界面的密度、分子和流速的分布,探究納米薄膜潤滑的潤滑機(jī)理和剪切速度對界面滑移的影響。結(jié)果表明,納米間隙中潤滑劑存在分層現(xiàn)象,各層間距相近,并非越遠(yuǎn)離固體壁面層間距越大,層間距約為角鯊?fù)榉肿訂蝹€C C鍵距離的3~4倍;隨著油膜厚度的減小,納米間隙中潤滑劑層狀分布越明顯,固化層密度越大;當(dāng)油膜厚度為3.44 nm時,固液界面滑移現(xiàn)象明顯,滑移長度b值隨著V先增大后減小,當(dāng)V 為22.8 m·s-1時,b達(dá)到最大值4.35 nm;根據(jù)模擬和計算結(jié)果,給出滑移長度與剪切速度的關(guān)系公式。

    關(guān)鍵詞:分子模擬;納米結(jié)構(gòu);聚合物;界面滑移;薄膜潤滑;角鯊?fù)?/p>

    2015-06-01收到初稿,2015-09-10收到修改稿。

    聯(lián)系人:高誠輝。第一作者:潘伶(1969—),女,博士,副教授。

    Received date: 2015-06-01.

    Foundation item: supported by the National Natural Science Foundation of China (51175085),the Tribology Science Found of State Key Laboratory of Tribology (SKLTKF13A09) and the Natural Science Foundation of Fujian Province (2016J01226).

    引 言

    無滑移邊界條件是經(jīng)典流體力學(xué)和潤滑理論中的重要假設(shè),廣泛應(yīng)用于流體的理論分析、實驗和實際工程。然而,近年來的研究發(fā)現(xiàn)界面滑移在許多情況下有可能發(fā)生,特別是存在于微納米尺度下的流體。界面滑移的相關(guān)研究方法主要有理論數(shù)值分析、實驗和分子動力學(xué)(MD)模擬。理論數(shù)值分析的物理模型主要有線性滑移長度模型[1]、極限剪應(yīng)力模型[2-3]和非線性滑移長度模型[4-5]。實驗方法常用的有原子力顯微鏡(AFM)[6]和表面力儀(SFA)[7]等。AFM方法是通過探針測量擠壓流中流體阻力和膜厚的關(guān)系,再根據(jù)理論公式分析計算測量結(jié)果,而界面滑移的理論和數(shù)值分析還不完善。

    目前,納米薄膜潤滑的研究仍處于初級階段,其潤滑機(jī)理、數(shù)學(xué)模型、微觀分子分布和油膜厚度等的理論分析、模擬計算和實驗檢測都與工程實際有較大差距。對薄膜潤滑進(jìn)行深入研究并開發(fā)相應(yīng)的實際應(yīng)用技術(shù)有著重要的意義。隨著計算機(jī)計算能力、并行技術(shù)以及相關(guān)勢函數(shù)的發(fā)展,分子動力學(xué)模擬方法應(yīng)用日益廣泛[8-11]。采用MD模擬納米間隙潤滑劑,可以動態(tài)適時地顯示全油膜中分子和密度的分布規(guī)律,克服了實物實驗無法測量高速加載時潤滑薄膜的密度和原子運動的缺點。另外,實驗的方法不能檢測新型的尚未研制出來的潤滑劑,而MD方法則可以靈活準(zhǔn)確地建立各種分子結(jié)構(gòu)模型。

    雖然在過去的20年里,納米流體的界面滑移和薄膜潤滑機(jī)理是研究的熱點,但以往的研究,無論是實驗、MD模擬還是數(shù)值分析的結(jié)論多數(shù)是建立在較為簡化的模型和推理上[12-20]。以往的MD研究多數(shù)采用簡單的單粒子硬球模型[5,21-23],簡單的LJ勢函數(shù)描述分子間的相互作用[5,21-22,24],或通過改變某一參數(shù)來模擬不同的濕潤性、固液相互作用勢和硬度等[22-23]對界面滑移的影響,而很少涉及液體的具體分子組成、結(jié)構(gòu)和更為準(zhǔn)確的分子力場。

    本文采用聚合物一致性力場(PCFF)模擬研究納米間隙內(nèi)角鯊?fù)樵诠腆w壁面不同的剪切速度和油膜厚度下的分子、密度及流速分布,探究納米間隙中固液界面的界面滑移及潤滑機(jī)理。

    1 模型的建立

    圖1是納米間隙潤滑劑剪切模型,納米級角鯊?fù)闈櫥な芟抻趦善叫械墓腆w壁面。圖中球形表示原子,橙色球表示鋁Al原子,藍(lán)色球表示碳(C)原子,綠色球表示氫(H)原子。固體壁面采用面心立方fcc晶體鋁Al(001)晶面,因為Al原子的Lennard-Jones勢函數(shù)中的能量參數(shù)e=3.32320 kcal·mol-1(1 cal=4.1868 J),距離參數(shù)s=2.9964?(1 ?=0.1 nm),其能量參數(shù)比Pt和Fe等原子小,固液之間的相互作用勢較小,使得常壓下界面滑移現(xiàn)象比較明顯,而且Al也是常用的金屬材料。在保證計算誤差小于3%的情況下,為減少計算機(jī)集群的計算量,取Al晶體表面尺寸為6.48 nm×4.05 nm,晶格常數(shù)為0.405 nm。上、下固體壁面施加同樣大小的反向速度V,V分別取7種介于3~30 m·s-1不同的速度。常壓下油膜厚度h分別取3.44、5.13 和7.02 nm。

    圖1 不同剪切速度下納米間隙潤滑劑模型Fig.1 Model of confined thin film at different shear velocities

    當(dāng)h=3.44 nm時,模型中總共有13408個原子,其中有3840個Al原子、3120個C原子和6448個H原子,即潤滑劑中共有104個角鯊?fù)榉肿?;?dāng)h=5.13 nm時,模型中總共有18100個原子,潤滑劑中共有155個角鯊?fù)榉肿?;?dāng)h=7.02 nm時,模型中總共有23344個原子,潤滑劑中共有212個角鯊?fù)榉肿印?/p>

    常用于美容產(chǎn)品和醫(yī)藥保健品的角鯊?fù)?,還可以作為電接觸或機(jī)械接觸潤滑油的基礎(chǔ)油,無毒且生物降解性高,是良好的環(huán)境友好潤滑劑。角鯊?fù)榈姆肿邮綖镃30H62,分子結(jié)構(gòu)和模型如圖2所示。

    2 模擬方法

    采用大規(guī)模原子/分子并行模擬器(large-scale atomic/molecular massively parallel simulator,LAMMPS[25])編程模擬納米間隙潤滑劑的剪切流動,分子力場采用聚合物一致性力場(polymerconsistent force field,PCFF),其力場參數(shù)是由大量的實驗數(shù)據(jù)和精確的量子計算結(jié)果擬合得出,因此計算量大,模擬也更準(zhǔn)確,適用于計算聚合物、有機(jī)物和生物分子體系,并含有一些金屬原子的力場參數(shù)。 PCFF勢函數(shù)包含非鍵結(jié)作用項、鍵伸縮項、鍵角彎曲項、二面角扭曲項、離平面振動項和多種交叉作用項,其勢函數(shù)和各參數(shù)的選取詳見文獻(xiàn)[25-28]。

    模型中油膜厚度方向即z方向采用自由邊界條件,其他兩個方向都設(shè)置成周期性邊界條件。由于界面滑移受上、下固體壁面的變形影響較大,故將固體壁面設(shè)置成彈性可變形,壁厚為1.21 nm。非鍵結(jié)作用項中的LJ勢作用項和庫侖作用項的截斷半徑都取12.5 ?,因為當(dāng)油膜厚度為5.31 nm時,其計算結(jié)果與截斷半徑取16 ?時非常接近,誤差小于0.5%,但在美國SUN公司的Blade60000刀片服務(wù)器上單節(jié)點的計算時間可減少近1/3。分子動力學(xué)模擬時間步長取1 fs。采用共軛梯度法,通過迭代調(diào)整原子的位置以減小模擬體系的能量,確定合理的初始位形,其中能量偏差取10-6kcal·mol-1,力偏差取10-8kcal·(mol·?)-1。先在微正則系綜(NVE)下對模擬體系進(jìn)行0.5 ns的弛豫以達(dá)到293 K恒溫條件下的平衡。然后,在上、下固體壁面分別施加7種x向的方向相反的速度V,使納米間隙中潤滑劑受到1.8 ns的剪切作用。

    圖2 角鯊?fù)榉肿覨ig.2 Squalane molecule

    圖3 h=7.02 nm、V=15 m·s-1時吸附層分子Fig.3 Distributions of squalane molecules adjacent to the wall at h=7.02 nm and V=15 m·s-1

    根據(jù)模擬結(jié)果,采用簡單實用的線性滑移長度模型計算滑移長度。線性滑移長度模型[1]認(rèn)為固液交界處流體與固體的相對速度即滑移速度vs與流體在壁面的剪切速率g呈正比[29]

    式中,b為滑移長度,表示固液界面滑移的程度;vx是流體沿x方向的流動速度。

    3 結(jié)果與討論

    3.1潤滑劑分子和密度的分布

    圖3是h=7.02 nm、V=15 m·s-1時,上固體壁面和角鯊?fù)槲侥し肿釉诓煌羟兴矔rt的分布。由圖可見,開始剪切時即t=0時,有較多的潤滑劑分子接近垂直于固體表面分布,且跨越2~3層以上的潤滑劑聚集層。隨著剪切時間的增長,這些接近垂直分布的潤滑劑分子出現(xiàn)彎曲和扭轉(zhuǎn),并與固體表面形成較小的角度或長分子鏈中的一部分與固體表面平行。

    然而,以往的研究多數(shù)認(rèn)為流體分子在近壁面處主要是平行于壁面,距離壁面越遠(yuǎn)的分子越垂直于壁面[12,30-31]。以往的研究多采用Legendre多項式表示的方向參數(shù)S( z )

    式中,z是粒子質(zhì)心的z坐標(biāo),θ是鍵方向與z軸的夾角。當(dāng)鍵平行固體壁面時,q=90°,S( z )=-0.5;當(dāng)鍵垂直于固體壁面時,q=0°,S( z )= 1。鍵的方向是取兩相鄰成鍵原子之間的鍵的方向[12,30],或取幾個連續(xù)成鍵原子中某兩個原子中心的連線方向[31],或以分子中功能基團(tuán)定位,通過統(tǒng)計各區(qū)域的平均S( z )值來判斷分子的方向。方向參數(shù)S( z )只能表明分子中部分鍵的方向,適于研究短鏈分子的方向,而長鏈分子在運動過程中普遍出現(xiàn)多次彎曲和扭轉(zhuǎn)等,故S( z )不宜用于判斷長鏈分子的方向。

    當(dāng)剪切時間t=1.8 ns時,將模擬體系沿z方向的總尺寸L(圖1)分為200等分,以各原子中心所在位置為各原子的質(zhì)心,計算出各等分中潤滑劑的平均密度,并畫出密度分布曲線,如圖4(a)、(c)、(e)所示。由圖可見,3種油膜厚度下,由于受到固壁的約束,在固液交界處都有一層密度很大的潤滑劑層,即出現(xiàn)固化現(xiàn)象。油膜厚度越小固化層密度越大,固化層對界面滑移的影響也越大。

    圖4 沿油膜厚度方向角鯊?fù)橘|(zhì)量密度和部分分子的分布Fig.4 Density and molecular distributions of lubricants in the z direction

    當(dāng)h=3.44 nm時,納米間隙內(nèi)潤滑劑明顯為層狀分布并分為7層,如圖4(a)所示,圖中標(biāo)出了各層之間的相對距離,可以看出各層之間距離沒有明顯的規(guī)律,而以往其他學(xué)者的研究常認(rèn)為離壁面越遠(yuǎn)層間距越大[12-14]。實際上,在文獻(xiàn)[32-33]中,采用原子力顯微鏡(AFM)測量探針在趨近潤滑劑中的基片時所受的阻力,仔細(xì)觀察文中的圖就可以發(fā)現(xiàn)力曲線的層間距是不規(guī)則的,相同的現(xiàn)象在采用MD模擬的文獻(xiàn)[34]的圖中也可以觀察到。當(dāng)h=5.13 nm時,上、下固液交界處附近分別有明顯的3~4層潤滑劑,膜厚中部波動較小,如圖4(c)所示,同樣并非離壁面越遠(yuǎn)層間距越大;隨著油膜厚度的增加,潤滑劑分層現(xiàn)象越來越不明顯。當(dāng)h=7.02 nm時,只在上、下固液交界處附近有較明顯的3~4層潤滑劑,膜厚中部波動更小且無規(guī)則。另外,3種油膜厚度下,固體壁面附近的潤滑劑層間距基本相等,約為0.507 nm。角鯊?fù)榉肿訂蝹€鍵的距離約為0.154 nm,各潤滑劑聚集層中心間距離約為角鯊?fù)榉肿訂蝹€鍵距離的3~4倍。圖4(b)、(d)、(f)是用不同的顏色表示隨機(jī)抽取的少量分子,以清楚地顯示潤滑劑分子在納米間隙中的分布,由圖可見,同一分子常??缭蕉鄬?,而且出現(xiàn)彎曲和扭轉(zhuǎn)。

    以往的多數(shù)MD研究采用單粒子硬球模型,結(jié)果顯示在垂直于油膜厚度方向(即xy平面),潤滑劑分子在固液交界處呈有序化排列,x向密度分布呈周期性變化[13,19]。然而,研究結(jié)果顯示,雖然納米薄膜潤滑中,z方向即沿油膜厚度方向潤滑劑分層現(xiàn)象明顯,但x方向沒有分層的現(xiàn)象。圖5是固液交界處xy平面上部分角鯊?fù)榈姆肿臃植肌D6是角鯊?fù)樵趚方向密度分布曲線,可見曲線波動小,密度分布均勻,不存在周期性變化規(guī)律。

    圖5 固液交界面角鯊?fù)椴糠址肿拥姆植糉ig.5 Molecular distribution of squalane at solid- liquid interface

    圖6 x向角鯊?fù)橘|(zhì)量密度分布Fig.6 Density profile of squalane in x direction

    3.2界面滑移

    當(dāng)h=3.44 nm時,將整個模擬體系沿z方向分成34等分,圖7是當(dāng)V=10 m·s-1時,整個模擬體系沿x方向的速度分布。每個數(shù)據(jù)點的橫坐標(biāo)表示每等分中心的z坐標(biāo)在整個體系中的相對位置;縱坐標(biāo)表示每等分中原子x方向的平均速度。圖7中上、下兩條水平數(shù)據(jù)點分別表示上、下兩固體壁面的速度;中間的數(shù)據(jù)點是角鯊?fù)榱黧w的速度分布。由圖可見,固液交界處有明顯的相對速度即出現(xiàn)界面滑移,潤滑劑流速vx呈非線性分布。固體壁面的速度V與流體在固液交界處的速度vx之差即為滑移速度vs,圖中流體速度分布曲線在固液交界處的變化率除以L即為剪切速率g。

    圖7 V=10 m·s-1時x向速度分布Fig.7 Velocity in x direction at V=10 m·s-1

    圖8是h=3.44 nm、V=10 m·s-1時流體在3個方向的分段平均流速vx、vy和vz分布。與vx相比,vy和vz變化很小,特別是vz變化更小,說明在宏觀上流體屬于層流。圖8中的插圖是放大的vz分布圖,可以看出,在靠近固液交界處vz幾乎為0,而油膜中心附近vz稍大。

    圖8 V=10 m·s-1時角鯊?fù)?個方向的速度分布Fig.8 Velocity in three directions at V=10 m·s-1

    在固體壁面上分別施加7種不同的速度V,并取流體上半部分(z/L≥0.5)進(jìn)行分析,如圖9所示。由圖可見,當(dāng)V=3~30 m·s-1時,固液交界處的流體速度vx和速度變化率隨著V的增大而增大;當(dāng)V=25 m·s-1后,固液交界處的速度變化率明顯增大;同時,隨著V的增大,流體速度vx呈非線性分布越明顯。

    圖9 不同V時納米間隙上半部分流體的速度分布Fig.9 Velocities of upper fluids at different shear velocities

    圖10(a)、(b)中的離散數(shù)據(jù)點分別是根據(jù)圖9的數(shù)據(jù)計算出的固液界面的相對速度vs和剪切速率g,圖中曲線分別是對離散點進(jìn)行線性和3次多項式擬合后的曲線。圖10(c)中的離散數(shù)據(jù)點是根據(jù)圖10(a)、(b)中的離散數(shù)據(jù)點,由線性滑移長度模型的式(1)計算得出的滑移長度b,圖中曲線是二次多項式擬合曲線,其表達(dá)式為

    式(3)方便地表達(dá)了滑移長度b與剪切速度V的關(guān)系。由圖10可以看出,隨著V的增大,vs值基本呈線性增大,剪切速率g呈非線性增大且在V>20 m·s-1后增大更明顯,從而使得滑移長度b值隨著V先增大后減小。當(dāng)V=22.8 m·s-1時,b達(dá)最大值4.35 nm。

    4 結(jié) 論

    在PCFF力場下,采用MD方法模擬不同剪切速度下納米間隙角鯊?fù)榈募羟辛鲃?。分析了流體分子、密度和流速的分布,計算出固液界面的滑移長度,得出以下結(jié)論。

    (1)固液交界處液體分子形態(tài)各異,長鏈分子中的部分平行于固體壁面,部分垂直于固壁,還有部分與固壁呈一定角度。方向參數(shù)S(z)不宜作為長鏈分子方向的判定依據(jù)。

    (2)各種油膜厚度下,在固體壁面附近,都有一層密度很大的潤滑劑固化層。隨著油膜厚度的減小,納米間隙中潤滑劑層狀分布越明顯,固化層密度越大,固化層對界面滑移的影響也越大。層間距基本相等,并非離壁面越遠(yuǎn)層間距越大。層間距約為角鯊?fù)榉肿訂蝹€C C鍵距離的3~4倍。

    圖10 各參數(shù)隨V的變化Fig.10 Factors as functions of shear velocity

    (3)由于Al和squalane原子之間相互作用勢較小時,在常壓下就出現(xiàn)界面滑移現(xiàn)象。當(dāng)油膜厚度為3.44 nm時,隨著剪切速度V的增大,流速vx明顯呈非線性分布,固液交界處的相對速度vs基本呈線性增大,剪切速率g呈非線性增大且在V>20 m·s-1后增大更明顯,從而使得滑移長度b值隨著V先增大后減小。根據(jù)模擬和計算結(jié)果,給出了滑移長度與剪切速度的關(guān)系公式。

    本文的研究克服了實物實驗無法測量納米間隙中流體的密度及分子運動規(guī)律等的缺點,為采用分子動力學(xué)方法定量地、動態(tài)適時地預(yù)測納米間隙流體的密度和分子分布,進(jìn)一步探究納米薄膜潤滑機(jī)理、研究應(yīng)用界面滑移和開發(fā)新型潤滑劑提供了可靠的依據(jù)。

    符號說明

    b ——滑移長度,nm

    h ——油膜厚度,nm

    L ——模擬體系沿z向的總尺寸,nm

    t ——剪切時間,ns

    V ——剪切速度,m·s-1

    vs——滑移速度,m·s-1vx,vy,vz——分別為流體在x、y和z方向的分速度,

    m·s-1

    g——剪切速率,ns-1

    ε ——能量參數(shù),kcal·mol-1

    ρ ——密度,kg·m-3

    σ ——距離參數(shù),?

    References

    [1] OLDSTEIN S. Modern Developments in Fluid Dynamics [M]. Vol Ⅱ. Oxford: Clarendon Press,1957: 676-680.

    [2] BEAVERS G S,JOSEPH D D. Boundary conditions at a naturally permeable wall [J]. Journal of Fluid Mechanics,1967,30 (1): 197-207. DOI: 10.1017/S0022112067001375

    [3] MA G J,WU C W,ZHOU P. Wall slip and hydrodynamics of two-dimensional journal bearing [J]. Tribology International,2007,40: 1056-1066. DOI: 10.1016/j.triboint.2006.10.003.

    [4] THOMPSON P A,TROIAN S M. A general boundary condition for liquid flow at solid surfaces [J]. Nature,1997,389: 360-362. DOI: 10.1038/38686.

    [5]RIEZJEV N V,TROIAN S M. Molecular origin and dynamic behavior of slip in sheared polymer films [J]. Physical Review Letters,2004,92 (1): 018302. DOI: 10.1103/PhysRevLett.92.018302.

    [6] JING D,BHUSHAN B. Boundary slip of superoleophilic,oleophobic,and superoleophobic surfaces immersed in deionized water,hexadecane,and ethylene glycol [J]. Langmuir,2013,29 (47): 14691-14700. DOI: 10.1021/la4030876.

    [7] ESPINOSA-MARZAL R M,ARCIFA A,ROSSI A,et al. Microslips to “avalanches” in confined,molecular layers of ionic liquids [J]. The Journal of Physical Chemistry Letters,2013,5 (1): 179-184. DOI: 10.1021/jz402451v.

    [8] 陳其樂,孔憲,盧滇楠,等. 外壁電荷性質(zhì)對雙壁碳納米管中水分子運動行為的影響 [J]. 化工學(xué)報,2014,65 (1): 319-327. DOI: 10.3969/j.issn.0438-1157.2014.01.042. CHEN Q L,KONG X,LU D N,et al. Molecular simulation of outer surface charge on water transport through double-wall carbon nanotube [J]. CIESC Journal,2014,65 (1): 319-327. DOI: 10.3969/j.issn.0438-1157.2014.01.042.

    [9] 張程賓,趙沐雯,陳永平,等. 流體密度對納通道內(nèi)流動滑移的影響 [J]. 化工學(xué)報,2012,63 (S1): 12-16. DOI: 10.3969/j.issn.0438-1157.2012.zl.003. ZHANG C B,ZHAO M W,CHEN Y P,et al. Effects of fluid density on velocity slip in nanochannels [J]. CIESC Journal,2012,63 (S1): 12-16. DOI: 10.3969/j.issn.0438-1157.2012.zl.003.

    [10] 潘伶,高誠輝. 納米間隙潤滑劑季戊四醇四酯的壓縮性能分子動力學(xué)模擬 [J]. 機(jī)械工程學(xué)報,2015,51 (5): 76-82. DOI: 10.3901/JME.2015.05.076. PAN L,GAO C H. Molecular dynamics simulation on the compressibility of pentaerythritol tetra in nanogap [J]. Journal of Mechanical Engineering,2015,51 (5): 76-82. DOI: 10.3901/JME.2015.05.076.

    [11] PAN L,GAO C H. Confined fluid density of a pentaerythritol tetraheptanoate lubricant investigated using molecular dynamics simulation [J]. Fluid Phase Equilibria,2015,385: 212-218. DOI: 10.1016/j.fluid.2014.11.014.

    [12] TSIGE M,PATNAIK S S. An all-atom simulation study of the ordering of liquid squalane near a solid surface [J]. Chemical Physics Letters,2008,457 (4/5/6): 357-361. DOI: 10.1016/j.cplett.2008.04.026.

    [13] KARNIADAKIS G,BESKOK A,ALURU N. Microflows and Nanoflows: Fundamentals and Simulation [M]. Berlin Heidelberg,New York: Springer,2004: 1-23.

    [14] PERTSIN A J,GRUNZE M. Long-ranged solvation forces in a fluid with short-ranged interactions [J]. Journal of Chemical Physics,2003,118 (17): 8004-8009. DOI: 10.1063/1.1564051.

    [15] LOI S,SUN G,F(xiàn)RANZ V,et al. Rupture of molecular thin films observed in atomic force microscopy (II): Experiment [J]. Physical Review E,2002,66 (3): 031602. DOI: 10.1103/PhysRevE.66.031602.

    [16] FRANZ V,BUTT H-J. Confined liquids:? solvation forces in liquid alcohols between solid surfaces [J]. The Journal of Physical Chemistry B,2002,106 (7): 1703-1708. DOI: 10.1021/jp012541w.

    [17] 陳天星. 利用原子力顯微鏡對近壁面受限液體性質(zhì)的研究 [D].北京: 清華大學(xué),2011. CHEN T X. Study on properties of the confined liquids at solid surface with AFM[D]. Beijing: Tsinghua University,2011.

    [18] SIVEBAEK I M,SAMOILOV V N,PERSSON B N J. Squeezing molecularly thin alkane lubrication films: layering transitions and wear [J]. Tribology Letters,2004,16 (3): 195-200. DOI: 10.1023/B:TRIL.0000009730.31175.82.

    [19] 胡元中,王慧,郭炎. 超薄油膜潤滑的分子動力學(xué)模擬 [J]. 摩擦學(xué)學(xué)報,1995,15 (2): 138-144. DOI: 10.16078/j.tribology.1995.02.007. HU Y Z,WANG H,GUO Y. Molecular dynamics simulation of ultra thin film lubrication (I): Rigid molecule model [J]. Tribology,1995,15 (2): 138-144. DOI: 10.16078/j.tribology.1995.02.007.

    [20] 王慧,胡元中,郭炎. 超薄潤滑膜界面滑移現(xiàn)象的分子動力學(xué)研究 [J]. 清華大學(xué)學(xué)報 (自然科學(xué)版),2000,40 (4): 107-110. WANG H,HU Y Z,GUO Y. Molecular dynamics study of interfacial slip behavior of ultrathin lubricating films [J]. J. Tsinghua Univ. (Sci. & Tech.),2000,40 (4): 107-110.

    [21] NAGAYAMA G,CHENG P. Effects of interface wettability on microscale flow by molecular dynamic simulation [J]. International Journal of Heat and Mass Transfer,2004,47: 501-513. DOI: 10.1016/j.ijheatmasstransfer.2003.07.013.

    [22] ASPROULIS N,DRIKAKIS D. Boundary slip dependency on surface stiffness [J]. Physical Review E,2010,81 (6): 061503. DOI: 10.1103/PhysRevE.81.061503.

    [23] ZHANG H,ZHANG Z,YE H. Molecular dynamics-based prediction of boundary slip of fluids in nanochannels [J]. Microfluidics and Nanofluidics,2012,12 (1-4): 107-115. DOI: 10.1007/s10404-011-0853-y.

    [24] NOORIAN H,TOGHRAIE D,AZIMIAN A R. Molecular dynamics simulation of poiseuille flow in a rough nano channel with checker surface roughnesses geometry [J]. Heat and Mass Transfer,2014,50 (1): 105-113. DOI: 10.1007/s00231-013-1232-x.

    [25] PLIMPTON S. LAMMPS molecular dynamics simulator[EB/OL]. [2015-5-15]. http://lammps.sandia.gov.

    [26] PLIMPTON S. Fast parallel algorithms for short-range molecular dynamics [J]. Journal of Computational Physics,1995,117 (1): 1-19. DOI: 10.1006/jcph.1995.1039.

    [27] SUN H. Ab initio calculations and force field development for computer simulation of polysilanes [J]. Macromolecules,1995,28 (3): 701-712. DOI: 10.1021/ma00107a006.

    [28] SUN H. COMPASS:? an ab initio force-field optimized for condensed-phase applications overview with details on alkane and benzene compounds [J]. Journal of Physical Chemistry B,1998,102 (38): 7338-7364. DOI: 10.1021/jp980939v.

    [29] NAVIER C L M H. Memoire sur les du movement des fluids [J]. Mem l'Acad Roy Sci l’Inst France,1823,6: 389-440.

    [30] VADAKKEPATT A,DONG Y,LICHTER S,et al. Effect of molecular structure on liquid slip [J]. Physical Review E,2011,84 (6): 066311. DOI: 10.1103/PhysRevE.84.066311.

    [31] MENDONCA A C F,F(xiàn)OMIN Y D,MALFREYT P,et al. Novel ionic lubricants for amorphous carbon surfaces: molecular modeling of the structure and friction [J]. Soft Matter,2013,9 (44): 10606-10616. DOI: 10.1039/C3SM51689J.

    [32] LIM R,O’SHEA S J. Solvation forces in branched molecular liquids [J]. Physical Review Letters,2002,88 (24): 246101. DOI: 10.1103/PhysRevLett.88.246101.

    [33] LIM R,LI S F Y,O’SHEA S J. Solvation forces using sample-modulation atomic force microscopy [J]. Langmuir,2002,18 (16): 6116-6124. DOI: 10.1021/la011789+.

    [34] ZHENG X,ZHU H,KOSASIH B,et al. A molecular dynamics simulation of boundary lubrication: the effect of n-alkanes chain length and normal load [J]. Wear,2013,301 (1/2): 62-69. DOI: 10.1016/j.wear.2013.01.052.

    Molecular dynamics simulation of boundary slip in nanogap: effect of shear velocity

    PAN Ling,GAO Chenghui
    (School of Mechanical Engineering and Automation,F(xiàn)uzhou University,F(xiàn)uzhou 350108,F(xiàn)ujian,China)

    Abstract:Molecular dynamics (MD) simulations using the polymer consistent force field (PCFF) were adopted to investigate the density,molecular and velocity distributions of lubricant squalane in nanogap at 293 K,three different film thicknesses and a wide range of shear velocities. The lubrication mechanism and boundary slip were analyzed. The results showed that the lubricant atoms tended to form layers parallel to the confining wall. The distances between the layers of lubricant atoms were irregular rather than broadening far away from the walls and were about three to four times the length of C C bond in the squalane. The tendency of lubricant atoms to form layers and the density of solid-like layer increased with decreasing film thickness. It was clearly to find the boundary slip at the solid-liquid interface from the velocity profile. The slip lengths increased with increasing velocity of substrates at the beginning,and then decreased. When the film thickness was 3.44 nm,the maximum slip length was 4.35 nm at the substrate velocity of 22.8 m·s-1. According to the simulations,the relationship between the slip length and the shear velocity was given.

    Key words:molecular simulation; nanostructure; polymers; boundary slip; thin film lubrication; squalane

    DOI:10.11949/j.issn.0438-1157.20150755

    中圖分類號:TH 117.2

    文獻(xiàn)標(biāo)志碼:A

    文章編號:0438—1157(2016)04—1440—08

    基金項目:國家自然科學(xué)基金項目(51175085);清華大學(xué)摩擦學(xué)國家重點實驗室開放基金項目(SKLTKF13A09);福建省自然科學(xué)基金項目(2016J01226)。

    Corresponding author:Prof. GAO Chenghui,gch@fzu.edu.cn

    猜你喜歡
    聚合物
    注聚污水處理技術(shù)調(diào)研報告
    聚合物改善水泥石滲透性評價方法研究
    價值工程(2017年11期)2017-04-18 17:58:11
    高溫高鹽油藏用化學(xué)驅(qū)油劑的研究
    電化學(xué)法血糖儀電極的改性專利技術(shù)綜述
    環(huán)保型抗高溫抗鹽聚合物降粘劑的合成及評價
    科技傳播(2016年10期)2016-07-15 23:11:56
    采油生產(chǎn)過程中聚合物分散技術(shù)的應(yīng)用分析
    淺談聚合物驅(qū)油區(qū)塊轉(zhuǎn)油站加藥措施
    聚合物驅(qū)站場施工安全管理措施分析
    復(fù)合型高分子導(dǎo)電材料淺析
    磷酸鎂水泥的聚合物改性研究

    精品亚洲乱码少妇综合久久| 亚洲精品久久久久久婷婷小说| 欧美xxⅹ黑人| 咕卡用的链子| 男男h啪啪无遮挡| 超碰成人久久| 999久久久精品免费观看国产| 男女午夜视频在线观看| 久久久水蜜桃国产精品网| 两个人视频免费观看高清| 午夜福利在线观看吧| 18禁黄网站禁片午夜丰满| 一进一出抽搐动态| 日日摸夜夜添夜夜添小说| 97人妻天天添夜夜摸| 91字幕亚洲| 久久婷婷人人爽人人干人人爱 | 欧美日韩乱码在线| 在线观看免费午夜福利视频| 黄色a级毛片大全视频| 777久久人妻少妇嫩草av网站| 日韩精品青青久久久久久| 亚洲成人精品中文字幕电影| 精品久久久久久久毛片微露脸| 久久久久国内视频| 校园春色视频在线观看| 国产免费男女视频| 少妇 在线观看| 色在线成人网| 亚洲人成77777在线视频| 亚洲精品久久国产高清桃花| 午夜久久久在线观看| 国产亚洲精品综合一区在线观看 | 亚洲精品国产一区二区精华液| 亚洲精品av麻豆狂野| 色婷婷久久久亚洲欧美| 国产精品香港三级国产av潘金莲| 搡老妇女老女人老熟妇| 操美女的视频在线观看| 亚洲av片天天在线观看| ponron亚洲| АⅤ资源中文在线天堂| 午夜精品在线福利| 亚洲精品国产色婷婷电影| 香蕉久久夜色| 日韩免费av在线播放| 少妇被粗大的猛进出69影院| 免费看美女性在线毛片视频| 国产极品粉嫩免费观看在线| 亚洲色图综合在线观看| 老熟妇仑乱视频hdxx| 欧美午夜高清在线| 18美女黄网站色大片免费观看| 校园春色视频在线观看| 美国免费a级毛片| 在线观看www视频免费| 国产成人系列免费观看| 日本免费a在线| 国产熟女午夜一区二区三区| 精品久久蜜臀av无| 一边摸一边抽搐一进一出视频| 精品电影一区二区在线| 欧美日韩中文字幕国产精品一区二区三区 | 电影成人av| 可以免费在线观看a视频的电影网站| 中文字幕久久专区| 天天躁夜夜躁狠狠躁躁| 90打野战视频偷拍视频| aaaaa片日本免费| 欧美大码av| 手机成人av网站| 一级a爱片免费观看的视频| 嫩草影院精品99| 男人操女人黄网站| 看免费av毛片| 欧美成人午夜精品| av视频免费观看在线观看| cao死你这个sao货| 亚洲一区高清亚洲精品| 精品久久久精品久久久| 亚洲自拍偷在线| 欧美激情极品国产一区二区三区| 免费在线观看黄色视频的| 午夜免费鲁丝| 精品国产超薄肉色丝袜足j| 国产高清激情床上av| 久久中文字幕一级| 久久久国产成人精品二区| 波多野结衣巨乳人妻| 国产精品久久久久久亚洲av鲁大| 久久人妻福利社区极品人妻图片| 国内精品久久久久久久电影| 韩国av一区二区三区四区| 好男人电影高清在线观看| 久久天堂一区二区三区四区| 久久精品91无色码中文字幕| 制服人妻中文乱码| 日本黄色视频三级网站网址| АⅤ资源中文在线天堂| 视频在线观看一区二区三区| 亚洲情色 制服丝袜| 国产精品98久久久久久宅男小说| 免费少妇av软件| 久久中文字幕人妻熟女| 久久久精品国产亚洲av高清涩受| 亚洲 欧美一区二区三区| 亚洲欧美日韩另类电影网站| 国产精品久久久av美女十八| 亚洲aⅴ乱码一区二区在线播放 | 高清毛片免费观看视频网站| 国产aⅴ精品一区二区三区波| 午夜激情av网站| 99re在线观看精品视频| 久久久久久久精品吃奶| 亚洲免费av在线视频| 热re99久久国产66热| 成人av一区二区三区在线看| 国产三级黄色录像| 国产麻豆69| 国产1区2区3区精品| 久久婷婷成人综合色麻豆| 1024香蕉在线观看| 免费看十八禁软件| 久久 成人 亚洲| 亚洲国产欧美日韩在线播放| 久久久久久大精品| 欧美黄色淫秽网站| 熟妇人妻久久中文字幕3abv| 久热这里只有精品99| 欧美在线一区亚洲| 日韩欧美国产在线观看| 国产成人欧美| 成人亚洲精品一区在线观看| 国产成人一区二区三区免费视频网站| 美女国产高潮福利片在线看| 国产精品久久久久久精品电影 | 三级毛片av免费| 国产成人免费无遮挡视频| 国产欧美日韩一区二区三区在线| 99在线人妻在线中文字幕| 免费不卡黄色视频| 精品久久久精品久久久| 一级a爱片免费观看的视频| 男女之事视频高清在线观看| 国产精品亚洲av一区麻豆| netflix在线观看网站| 久久国产精品男人的天堂亚洲| 国产伦人伦偷精品视频| 中文字幕色久视频| 97人妻天天添夜夜摸| 日韩欧美一区二区三区在线观看| 亚洲五月婷婷丁香| 色综合婷婷激情| 99国产极品粉嫩在线观看| av在线天堂中文字幕| 欧美黄色淫秽网站| 色老头精品视频在线观看| 乱人伦中国视频| 国产亚洲精品久久久久5区| www国产在线视频色| 老司机深夜福利视频在线观看| 侵犯人妻中文字幕一二三四区| 97人妻天天添夜夜摸| 欧美一级毛片孕妇| 午夜久久久久精精品| 久久影院123| 在线国产一区二区在线| 每晚都被弄得嗷嗷叫到高潮| 一a级毛片在线观看| 国产午夜精品久久久久久| 国产成人精品在线电影| 日韩av在线大香蕉| 精品午夜福利视频在线观看一区| 国产亚洲精品久久久久久毛片| 欧美激情 高清一区二区三区| 男人操女人黄网站| 脱女人内裤的视频| 男女下面进入的视频免费午夜 | 日韩 欧美 亚洲 中文字幕| 99国产精品一区二区蜜桃av| 久久影院123| 国产成人欧美| 黄频高清免费视频| 中文字幕最新亚洲高清| 亚洲电影在线观看av| 国产成人精品在线电影| 极品人妻少妇av视频| 亚洲专区字幕在线| 久热这里只有精品99| 久久久久九九精品影院| 啦啦啦韩国在线观看视频| 999久久久国产精品视频| 国产av精品麻豆| 大码成人一级视频| 成年版毛片免费区| 国产精品爽爽va在线观看网站 | 国产av一区二区精品久久| 这个男人来自地球电影免费观看| 国产1区2区3区精品| 日韩欧美国产一区二区入口| 一本久久中文字幕| 亚洲国产欧美网| 亚洲国产精品sss在线观看| 日日干狠狠操夜夜爽| 国产av一区在线观看免费| 国产亚洲精品久久久久久毛片| 亚洲人成77777在线视频| 精品午夜福利视频在线观看一区| 亚洲电影在线观看av| 美女高潮到喷水免费观看| 久久九九热精品免费| 高潮久久久久久久久久久不卡| 麻豆久久精品国产亚洲av| 欧美日韩中文字幕国产精品一区二区三区 | 精品午夜福利视频在线观看一区| 久久国产精品影院| 亚洲熟妇熟女久久| 十分钟在线观看高清视频www| 777久久人妻少妇嫩草av网站| 99热只有精品国产| 女人被狂操c到高潮| 精品久久久久久久人妻蜜臀av | 久久九九热精品免费| 亚洲 欧美一区二区三区| 亚洲美女黄片视频| 亚洲色图 男人天堂 中文字幕| 免费一级毛片在线播放高清视频 | 黄色丝袜av网址大全| 国产99白浆流出| 一进一出好大好爽视频| 亚洲一区高清亚洲精品| 亚洲欧美激情在线| av在线天堂中文字幕| 国产精品 欧美亚洲| ponron亚洲| 亚洲欧美一区二区三区黑人| 亚洲情色 制服丝袜| 亚洲熟女毛片儿| 精品国产乱码久久久久久男人| svipshipincom国产片| 国产一区二区三区视频了| 黑丝袜美女国产一区| 在线av久久热| 18禁观看日本| 欧美在线一区亚洲| 精品久久久久久久人妻蜜臀av | 久久香蕉国产精品| 久久热在线av| 国产野战对白在线观看| 99riav亚洲国产免费| 男男h啪啪无遮挡| 91字幕亚洲| 成在线人永久免费视频| 丰满人妻熟妇乱又伦精品不卡| 丝袜在线中文字幕| 大型av网站在线播放| 成人18禁在线播放| 国产精品一区二区精品视频观看| 变态另类丝袜制服| 校园春色视频在线观看| 男人舔女人的私密视频| 女性生殖器流出的白浆| e午夜精品久久久久久久| 国产乱人伦免费视频| 美女大奶头视频| 男女之事视频高清在线观看| 亚洲欧美日韩高清在线视频| 国产真人三级小视频在线观看| 日韩大码丰满熟妇| 亚洲av成人一区二区三| 国产三级黄色录像| 国产精品一区二区在线不卡| 午夜影院日韩av| 一夜夜www| 99国产精品99久久久久| 少妇被粗大的猛进出69影院| 在线av久久热| 亚洲电影在线观看av| 免费人成视频x8x8入口观看| 精品国产乱子伦一区二区三区| 最好的美女福利视频网| 欧美色视频一区免费| 午夜日韩欧美国产| 国产亚洲欧美在线一区二区| 国产成人啪精品午夜网站| 999久久久国产精品视频| 久久久久久久久久久久大奶| 国产精品电影一区二区三区| 亚洲片人在线观看| 国内精品久久久久久久电影| 亚洲国产欧美日韩在线播放| 久久久久国内视频| 久久精品国产综合久久久| 巨乳人妻的诱惑在线观看| 99久久国产精品久久久| 国产精品98久久久久久宅男小说| av在线天堂中文字幕| 精品欧美一区二区三区在线| 国产精品久久久久久亚洲av鲁大| 午夜福利视频1000在线观看 | 久久这里只有精品19| 黑人巨大精品欧美一区二区mp4| 狠狠狠狠99中文字幕| 男男h啪啪无遮挡| 在线观看免费视频网站a站| 午夜精品久久久久久毛片777| 变态另类成人亚洲欧美熟女 | 黑人巨大精品欧美一区二区mp4| 亚洲av日韩精品久久久久久密| 自拍欧美九色日韩亚洲蝌蚪91| 久久欧美精品欧美久久欧美| tocl精华| 九色亚洲精品在线播放| 一个人免费在线观看的高清视频| 黄色毛片三级朝国网站| 午夜福利一区二区在线看| 久久这里只有精品19| 亚洲成人免费电影在线观看| 亚洲av熟女| 窝窝影院91人妻| avwww免费| 国产欧美日韩综合在线一区二区| 国产99白浆流出| cao死你这个sao货| 国产成+人综合+亚洲专区| 少妇被粗大的猛进出69影院| 亚洲第一电影网av| 国产在线观看jvid| 亚洲人成伊人成综合网2020| 在线观看66精品国产| 最近最新中文字幕大全免费视频| 日日摸夜夜添夜夜添小说| 少妇被粗大的猛进出69影院| 成熟少妇高潮喷水视频| www.自偷自拍.com| 69精品国产乱码久久久| 一二三四在线观看免费中文在| 叶爱在线成人免费视频播放| 国产免费男女视频| 满18在线观看网站| 亚洲一区中文字幕在线| 电影成人av| 亚洲午夜精品一区,二区,三区| 亚洲七黄色美女视频| 欧美绝顶高潮抽搐喷水| av视频在线观看入口| 欧美乱色亚洲激情| 男女之事视频高清在线观看| 伦理电影免费视频| 18禁国产床啪视频网站| 成人18禁高潮啪啪吃奶动态图| 精品国产美女av久久久久小说| 色尼玛亚洲综合影院| 美女午夜性视频免费| 亚洲第一欧美日韩一区二区三区| 一级a爱片免费观看的视频| 国产成人一区二区三区免费视频网站| 日本五十路高清| 三级毛片av免费| 大型黄色视频在线免费观看| 国产麻豆69| 午夜福利高清视频| 99香蕉大伊视频| 9热在线视频观看99| 日韩三级视频一区二区三区| 亚洲第一电影网av| 欧美黑人精品巨大| 国产精品秋霞免费鲁丝片| 男人舔女人的私密视频| 麻豆国产av国片精品| 亚洲欧美精品综合一区二区三区| 日本a在线网址| 亚洲精品av麻豆狂野| 欧美黄色片欧美黄色片| 久久久久精品国产欧美久久久| 日本欧美视频一区| 欧美日韩乱码在线| 制服人妻中文乱码| 亚洲欧美日韩另类电影网站| 午夜影院日韩av| 90打野战视频偷拍视频| 亚洲欧美激情综合另类| 精品福利观看| 亚洲无线在线观看| 亚洲欧洲精品一区二区精品久久久| 欧美+亚洲+日韩+国产| 国产精品 国内视频| 亚洲中文日韩欧美视频| 久热这里只有精品99| 老司机午夜十八禁免费视频| 亚洲av五月六月丁香网| АⅤ资源中文在线天堂| 成人亚洲精品av一区二区| 纯流量卡能插随身wifi吗| 国产精品久久久久久人妻精品电影| 黑人欧美特级aaaaaa片| 久久欧美精品欧美久久欧美| 亚洲av日韩精品久久久久久密| 成年女人毛片免费观看观看9| 九色亚洲精品在线播放| 十八禁网站免费在线| 欧美中文综合在线视频| 免费在线观看影片大全网站| 夜夜爽天天搞| 国产熟女xx| 男女床上黄色一级片免费看| 熟女少妇亚洲综合色aaa.| 嫩草影视91久久| 嫁个100分男人电影在线观看| 美国免费a级毛片| 91精品三级在线观看| 少妇的丰满在线观看| 亚洲午夜理论影院| 久久中文看片网| 亚洲avbb在线观看| 国产精品影院久久| 免费搜索国产男女视频| 91在线观看av| 亚洲精品美女久久久久99蜜臀| 亚洲国产中文字幕在线视频| 国产精品久久久久久亚洲av鲁大| 久久久久九九精品影院| 又黄又爽又免费观看的视频| 琪琪午夜伦伦电影理论片6080| 很黄的视频免费| 女性被躁到高潮视频| 1024香蕉在线观看| 国产欧美日韩一区二区精品| 脱女人内裤的视频| 亚洲精品久久成人aⅴ小说| 久久这里只有精品19| 亚洲一区二区三区色噜噜| 最近最新中文字幕大全免费视频| 88av欧美| 亚洲精品国产一区二区精华液| 日本五十路高清| 麻豆国产av国片精品| 精品少妇一区二区三区视频日本电影| 欧美日韩亚洲综合一区二区三区_| 精品日产1卡2卡| 精品国产超薄肉色丝袜足j| 国产伦人伦偷精品视频| www.精华液| 国产99白浆流出| 在线观看午夜福利视频| 热re99久久国产66热| АⅤ资源中文在线天堂| 国产精品野战在线观看| 欧美黑人欧美精品刺激| 免费久久久久久久精品成人欧美视频| 美女 人体艺术 gogo| 久久久久久免费高清国产稀缺| 国产精品久久视频播放| 每晚都被弄得嗷嗷叫到高潮| 国产欧美日韩综合在线一区二区| 韩国精品一区二区三区| av福利片在线| 亚洲专区中文字幕在线| 亚洲人成伊人成综合网2020| 91成人精品电影| 校园春色视频在线观看| 日韩欧美一区视频在线观看| 级片在线观看| 免费少妇av软件| 中文字幕最新亚洲高清| 国产成年人精品一区二区| 久久人人爽av亚洲精品天堂| 亚洲av日韩精品久久久久久密| 亚洲中文日韩欧美视频| 亚洲九九香蕉| 国产亚洲精品一区二区www| 久久人妻福利社区极品人妻图片| 精品一区二区三区av网在线观看| 国产一区二区三区综合在线观看| 中文字幕av电影在线播放| 免费观看人在逋| 色精品久久人妻99蜜桃| netflix在线观看网站| 别揉我奶头~嗯~啊~动态视频| 咕卡用的链子| 无人区码免费观看不卡| 亚洲av美国av| 久久久久国产一级毛片高清牌| 在线观看66精品国产| 国产精品av久久久久免费| 欧美+亚洲+日韩+国产| 色综合欧美亚洲国产小说| 色在线成人网| 国产精品一区二区免费欧美| 色播亚洲综合网| 亚洲自拍偷在线| 久久国产精品人妻蜜桃| 国产av在哪里看| 亚洲一卡2卡3卡4卡5卡精品中文| 欧美在线黄色| 又黄又爽又免费观看的视频| 少妇的丰满在线观看| 欧美亚洲日本最大视频资源| 久久久国产成人精品二区| 国产不卡一卡二| 在线国产一区二区在线| 如日韩欧美国产精品一区二区三区| 亚洲一区二区三区色噜噜| 午夜精品在线福利| 村上凉子中文字幕在线| 咕卡用的链子| 日韩高清综合在线| 97碰自拍视频| 又黄又爽又免费观看的视频| 欧美日韩亚洲国产一区二区在线观看| 国产欧美日韩精品亚洲av| 午夜福利高清视频| 国产免费av片在线观看野外av| 久久天堂一区二区三区四区| 国产精品av久久久久免费| 亚洲中文字幕日韩| 黑人巨大精品欧美一区二区蜜桃| 18禁观看日本| 国产成+人综合+亚洲专区| 欧美av亚洲av综合av国产av| 韩国精品一区二区三区| 脱女人内裤的视频| 少妇的丰满在线观看| 亚洲七黄色美女视频| 丁香欧美五月| 成人三级黄色视频| 日韩中文字幕欧美一区二区| 亚洲欧洲精品一区二区精品久久久| 18禁美女被吸乳视频| 美女高潮到喷水免费观看| www.999成人在线观看| 日韩欧美免费精品| 电影成人av| 亚洲免费av在线视频| 免费在线观看日本一区| 国产精品99久久99久久久不卡| 国产一区二区三区视频了| 欧美av亚洲av综合av国产av| 美女免费视频网站| АⅤ资源中文在线天堂| 亚洲欧美一区二区三区黑人| 香蕉国产在线看| 国产单亲对白刺激| 欧美乱色亚洲激情| 黄色成人免费大全| 亚洲国产欧美日韩在线播放| 悠悠久久av| 村上凉子中文字幕在线| 国产不卡一卡二| 国产午夜精品久久久久久| 亚洲第一青青草原| 黄片小视频在线播放| 精品人妻1区二区| 亚洲精品久久成人aⅴ小说| 亚洲国产精品999在线| 午夜福利欧美成人| av在线播放免费不卡| 久久天堂一区二区三区四区| 欧美日韩中文字幕国产精品一区二区三区 | 曰老女人黄片| 精品不卡国产一区二区三区| 一区二区三区高清视频在线| 在线十欧美十亚洲十日本专区| www.自偷自拍.com| av中文乱码字幕在线| 久久久久久大精品| 国内毛片毛片毛片毛片毛片| 99riav亚洲国产免费| 99精品久久久久人妻精品| 国产亚洲av嫩草精品影院| 亚洲一区高清亚洲精品| 中出人妻视频一区二区| 桃色一区二区三区在线观看| 久久人人爽av亚洲精品天堂| 久久久久精品国产欧美久久久| 最新美女视频免费是黄的| 精品久久久久久成人av| 亚洲男人的天堂狠狠| 99久久国产精品久久久| 国产成人影院久久av| 90打野战视频偷拍视频| 18禁国产床啪视频网站| 一个人免费在线观看的高清视频| 99精品欧美一区二区三区四区| 亚洲美女黄片视频| 97人妻精品一区二区三区麻豆 | 精品熟女少妇八av免费久了| 免费在线观看黄色视频的| 法律面前人人平等表现在哪些方面| 男人操女人黄网站| 免费在线观看黄色视频的| 一级毛片高清免费大全| 国产亚洲精品av在线| 动漫黄色视频在线观看| 亚洲精品一区av在线观看| 少妇粗大呻吟视频| 在线观看一区二区三区| 丁香欧美五月| 欧美老熟妇乱子伦牲交| 黑丝袜美女国产一区| 女生性感内裤真人,穿戴方法视频| 国产免费av片在线观看野外av| 久久精品aⅴ一区二区三区四区| 日本五十路高清| 免费高清视频大片| 制服丝袜大香蕉在线| 国产成人啪精品午夜网站| 国语自产精品视频在线第100页| 男女午夜视频在线观看| 欧美成人免费av一区二区三区| 亚洲一区高清亚洲精品| 免费不卡黄色视频| www国产在线视频色| 成人18禁在线播放| 757午夜福利合集在线观看|