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

    顆粒-顆粒接觸力的熱力學(xué)模型?

    2018-03-26 22:06:28蔣亦民劉佑
    物理學(xué)報(bào) 2018年4期
    關(guān)鍵詞:勢(shì)能熱力學(xué)力學(xué)

    蔣亦民 劉佑

    1)(中南大學(xué)物理與電子學(xué)院,長(zhǎng)沙 410083)

    2)(蒂賓根大學(xué)理論物理研究所,德國(guó),蒂賓根 72076)

    1 引 言

    顆粒物質(zhì)(諸如沙堆等)與傳統(tǒng)晶體材料(例如半導(dǎo)體硅)的重要不同是多出一個(gè)介觀(或細(xì)觀)顆粒尺度.這使得從電子-原子核之間的二體庫侖靜電力出發(fā),計(jì)算其物理力學(xué)性質(zhì)的所謂“從頭算(ab initio)”思想,雖然原則概念上是存在的,但實(shí)際操作上沒有可能性.也就是說對(duì)硅晶體可以從頭計(jì)算,但對(duì)顆粒物質(zhì)最小尺度的分析計(jì)算只能從細(xì)觀顆粒出發(fā).另外為了計(jì)算數(shù)目龐大的顆粒系統(tǒng),還不能像赫茲那樣將顆粒處理成可變形的彈性固體[1],而是必須簡(jiǎn)化為只有六個(gè)運(yùn)動(dòng)學(xué)(kinematic)變量的、帶轉(zhuǎn)動(dòng)的“質(zhì)點(diǎn)粒子”.如果要求顆粒i的六個(gè)運(yùn)動(dòng)學(xué)量,即質(zhì)心位矢和角位矢{ri(t),φi(t)}遵守牛頓力和力矩方程的話,則需要建立{ri,φi}與其他顆粒j之間的作用力fij關(guān)系模型來閉合牛頓方程組,然后對(duì)其數(shù)值求解,可得到該質(zhì)點(diǎn)系統(tǒng)細(xì)觀結(jié)構(gòu)和力結(jié)構(gòu)的時(shí)間演化(我們將省略與本文內(nèi)容無關(guān)的邊界條件).如果進(jìn)一步選擇合適的物理體積元做這些細(xì)觀信息的統(tǒng)計(jì)平均(諸如粗?;骄?,還可得各種宏觀平均量及其漲落的時(shí)空變化以及各種關(guān)聯(lián)性質(zhì).這套稱作軟球離散元(DEM)模擬方法的一個(gè)優(yōu)勢(shì)是可以描述以長(zhǎng)程關(guān)聯(lián)為特征的彈性力鏈的呈現(xiàn)與消亡現(xiàn)象,又稱固液轉(zhuǎn)變或jamming轉(zhuǎn)變.也就是說軟球DEM適用于顆粒物質(zhì)的所有類固液氣態(tài)及其轉(zhuǎn)變[2].

    軟球DEM面臨的顆粒-顆粒接觸作用力建模問題,在晶體物性的“從頭算”領(lǐng)域里是不存在的.“從頭算”認(rèn)為庫侖靜電力能夠解決所有凝聚態(tài)物理問題[3,4],但對(duì)顆粒物理而言接觸力建模才是它的一個(gè)基礎(chǔ)課題.由于這是個(gè)帶耗散的不可逆非保守力,原則上需要在熱力學(xué)概念基礎(chǔ)上建模,具體地講就是除了接觸力外,還應(yīng)該給出機(jī)械能功率與熱功率的比例.或者給出計(jì)算彈性勢(shì)能的表達(dá)式,因?yàn)榻佑|時(shí)的總能量變化功率可以將作用力乘以運(yùn)動(dòng)學(xué)速度得到:fij(vi?vj)(這里vi≡dri/dt),動(dòng)能的計(jì)算公式又是熟知的,故一旦明確了接觸過程彈性勢(shì)能的計(jì)算方式,問題也就徹底解決了.注意著名的赫茲法向接觸力fn~|ri?rj|3/2[1],是忽略耗散時(shí)的結(jié)果.有耗散時(shí),特別是有塑性耗散時(shí)(這在切向運(yùn)動(dòng)時(shí)必須考慮),接觸力是不能簡(jiǎn)單地表達(dá)為顆粒幾何運(yùn)動(dòng)學(xué)變量的代數(shù)函數(shù)型這樣的捆綁關(guān)系的.這個(gè)復(fù)雜性來自接觸力中有彈性貢獻(xiàn),而彈性是需要用額外的獨(dú)立變量來表征的概念[3,4]:即所謂的平移對(duì)稱破缺后呈現(xiàn)的Goldstone量(彈應(yīng)變或彈簧長(zhǎng)度變化).很多模型都注意到了有塑性時(shí)彈簧長(zhǎng)度變化與幾何運(yùn)動(dòng)學(xué)相對(duì)量{ri?rj,φi?φj}的關(guān)系需要用微積分來描述,特別是切向.例如Luding模型在處理切向力時(shí),就引入了一個(gè)切向彈簧長(zhǎng)度變量,它與運(yùn)動(dòng)學(xué)量的聯(lián)系是利用增量關(guān)系模型和一個(gè)歷史累積變量來定義[5].有些人不用彈簧長(zhǎng)度,直接建立切向力與運(yùn)動(dòng)學(xué)量的增量關(guān)系模型(參見文獻(xiàn)[2])(工程領(lǐng)域常常將某些微積分方程寫成增量關(guān)系形式,這個(gè)習(xí)慣很可能源自Truesdell[6]).長(zhǎng)期以來,接觸力建模研究一直處于這樣的唯象力學(xué)階段,并未去澄清熱力學(xué)量,即熱功率或彈性勢(shì)能的定義.值得提到的是,很多(例如Hertz-Mindlin)三維接觸力模型的可逆部分采用了虎克力乘以重疊量平方根的Boussinesq模型[7,8].這里雖然不涉及耗散,由于Boussinesq力不是保守的,我們?nèi)詿o法知道彈性勢(shì)能.

    由于這些問題,雖然有少數(shù)DEM工作計(jì)算了彈性勢(shì)能與動(dòng)能比值[9,10],一個(gè)嚴(yán)格熱力學(xué)的接觸力模型仍值得建立.考慮到目前一個(gè)嚴(yán)格熱力學(xué)的顆粒物質(zhì)連續(xù)力學(xué)理論——顆粒固體流體動(dòng)力學(xué)(GSH)已經(jīng)存在[11?14],接觸力的熱力學(xué)模型完全可以用類似的方案得到(第二節(jié)).這時(shí)GSH中處理庫侖屈服、塑性、滯迴等復(fù)雜力學(xué)現(xiàn)象的成功方法,可以直接繼承到接觸力上.例如前面提到的Boussinesq力沒有勢(shì)能的問題,完全可以連同庫侖屈服一起,用一個(gè)類似GSH的勢(shì)能模型來解決[15].由于各種能量耗散已經(jīng)在細(xì)觀接觸的層面上明確定義,我們相信基于它的DEM計(jì)算能更好地描述長(zhǎng)時(shí)間接觸力鏈的熱損耗和彈性勢(shì)能及其演化(今后擬開展的工作).本文的第三節(jié)將針對(duì)二體力特有的碰撞恢復(fù)系數(shù)隨碰撞初速度變化問題,做些簡(jiǎn)單的數(shù)值計(jì)算,以表明熱力學(xué)方法能夠計(jì)入這個(gè)實(shí)驗(yàn)現(xiàn)象.目前常用的Hertz-Mindlin和Luding模型的恢復(fù)系數(shù)都是材料常數(shù),這個(gè)變化關(guān)系長(zhǎng)期以來一直未能被描述(參見文獻(xiàn)[2]).鑒于能量和耗散與接觸力的關(guān)系比較復(fù)雜,我們將盡可能地給出詳細(xì)解釋和討論,以方便讀者了解熱力學(xué)能量方法與以往的純力學(xué)唯象方法的區(qū)別(第四節(jié)).

    2 顆粒作用力模型

    如圖1,兩個(gè)質(zhì)量m1,m2,半徑a1,a2的球形顆粒,實(shí)驗(yàn)室坐標(biāo)下的質(zhì)心位矢是r1,r2,相互作用力f.記相對(duì)位矢r=r1?r2,(法向)重疊量δ=d?r,其中r=|r|,d=a1+a2.兩球位置的連線方向,即法向單位矢量n=r/r=(nx,ny).切向單位矢量t=(?ny,nx).記f的法向和切向投影為fn,t,有

    為簡(jiǎn)單起見,本文假設(shè)圖1的所有矢量都位于二維平面,并且忽略轉(zhuǎn)動(dòng)(即局限于考慮沒有轉(zhuǎn)動(dòng)的二維運(yùn)動(dòng),這些簡(jiǎn)化不影響將討論的內(nèi)容和概念).這時(shí)兩球的運(yùn)動(dòng)將由牛頓方程

    圖1 兩個(gè)運(yùn)動(dòng)顆粒的坐標(biāo)和力示意圖Fig.1.Coordinates and forces of two grains.

    所支配.其中的v1,2=dr1,2/dt是速度,f1,2外力是兩球受到的外力(考慮為已知,例如重力或邊界墻力).另外記兩球的相對(duì)速度為v=v1?v2,它的法向投影vn=vn,切向投影vt=vt.兩球的切向相互位移是δt.用點(diǎn)表示時(shí)間導(dǎo)數(shù),例如˙δ=dδ/dt.由(2)和(3)式可得相對(duì)速度的運(yùn)動(dòng)方程:

    其中m=m1m2/(m1+m2)是折合質(zhì)量.

    顯然一旦有了力fn,t=fn,t(r1,r2,v1,v2)的模型公式,上面的牛頓方程組就閉合了(外力是已知的).所謂的軟球DEM模擬,就是在合適的邊界條件下,對(duì)這樣二體相互作用的大量顆粒系統(tǒng)進(jìn)行數(shù)值求解.注意fn,t的建模不是一件容易的事情,因?yàn)槌吮J氐膹椈尚妥饔昧ν?還需要考慮諸如庫侖屈服打滑(slip)或塑性(plasticity)、黏力(dashpot)、滯迴(hysteresis)等一系列的不可逆耗散現(xiàn)象.我們面臨的不是一個(gè)哈密頓保守力學(xué)問題,而是復(fù)雜的動(dòng)-靜摩擦力模型.另外根據(jù)赫茲的分析,法向力fn中的彈簧力部分肯定不是線性的.如引言里提到的,目前常用的DEM力模型的一個(gè)問題是能量的耗散情況不明確.例如含黏力Hertz-Mindlin模型:

    其中的kn,t,ηn,t,μ是材料常數(shù).文獻(xiàn)中并沒有明確給出相應(yīng)的彈簧勢(shì)能w彈的計(jì)算公式.這時(shí)如果期望分析計(jì)算力鏈的彈性勢(shì)能就沒有標(biāo)準(zhǔn)了.也許一個(gè)最好的做法是取

    沒有滑動(dòng)時(shí),它對(duì)δt的導(dǎo)數(shù)的確是ft公式(6)的彈性部分,但出現(xiàn)打滑就不好說了.另外它對(duì)δ的導(dǎo)數(shù)也不是fn公式(5)的彈性部分.總之對(duì)上述Hertz-Mindlin模型,我們覺得是不好定義彈簧勢(shì)能.實(shí)際上Hertz-Mindlin模型可解釋為用Boussinesq力描述法向和切向彈簧力,再加上黏力和庫侖屈服.由于Boussinesq力不是保守力,因而無從知道相應(yīng)的彈性勢(shì)能.為改善以往顆粒-顆粒作用力模型的這個(gè)能量缺陷,下面將模仿GSH的二級(jí)不可逆耗散圖像(見圖2),提出一套基于能量守恒的建模方法.

    一般而言如果要求彈性能w彈只是幾何相對(duì)位置r(t)的代數(shù)函數(shù),兩球之間不會(huì)有滑動(dòng)或塑性.也就是說含塑性變形的運(yùn)動(dòng)學(xué)(kinematic)變量不能直接用來描述彈性,需要額外引入表征接觸處彈性變形或彈簧長(zhǎng)度變化的位移矢量:

    圖2 顆粒-顆粒作用時(shí)的二級(jí)不可逆熱力學(xué)示意圖,圖中的三個(gè)能量之和是守恒的Fig.2.Two-stage-irreversibility of grain-grain interaction.Sum of the three energies is conserved.

    作為w彈的自變量(注意對(duì)本文考慮的二維運(yùn)動(dòng),只討論矢量的兩個(gè)分量即可),即

    彈 性 能w彈的 導(dǎo) 數(shù)πn=??w彈/?un和πt=??w彈/?ut將給出顆粒-顆粒相互作用力fn,t的彈性部分.可以不失一般地將它們寫成:

    這里的上標(biāo)“熱”和“松”,分別表示圖2中左右兩種耗散機(jī)制的力貢獻(xiàn).類似地可將彈性變形un,ut的運(yùn)動(dòng)方程寫成

    其中R松≥0是它的激發(fā)功率,I≥0是它向熱能弛豫的衰減功率.松動(dòng)能的激發(fā)功率是與πn,t的乘積,加上與vn,t的乘積:

    相應(yīng)地,有熱功率

    其中的I是來自緩沖區(qū)w松的貢獻(xiàn).最后,圖2中機(jī)械能的定義是兩個(gè)顆粒的動(dòng)能與彈性能之和:

    當(dāng)沒有外力時(shí),能量守恒要求機(jī)械能w機(jī)械、松動(dòng)能w松和熱能之和保持常數(shù),或它們的時(shí)間導(dǎo)數(shù):

    這可以直接計(jì)算證明.的確微分方程(16)有

    代入(9)和(10)式,有

    將(18)和(13),(15)式相加,即得(17)式.

    與以往的模型對(duì)比,基于方程(7)—(17)理論構(gòu)架下的接觸相互作用除了給出彈性勢(shì)能和熱功率外,在處理彈性變形量與幾何運(yùn)動(dòng)學(xué)量的關(guān)系上,以及機(jī)械能與熱能的耗散關(guān)系上,引入了不那么僵硬的微分方程和能量緩沖區(qū)等措施.這將使得理論具有足夠的靈活性來描述前面提到的各種復(fù)雜的接觸力學(xué)現(xiàn)象(庫侖屈服、打滑、塑性、黏力、滯迴等).

    將取類似文獻(xiàn)[15,16]的彈性勢(shì)能模型

    其中k,ξ,c是材料常數(shù),c是描述濕黏力的常數(shù)[16].注意un和ut分別是法向和切向彈簧長(zhǎng)度變化,量綱是米.彈性能量w彈量綱焦耳.彈簧系數(shù)k量綱為這里考慮了非線性的赫茲3/2冪率彈簧力.勢(shì)能(19)給出的彈簧力是

    這里要求法向彈簧的長(zhǎng)度變化總是小于零:un<0,但剪切彈簧的ut可正可負(fù).后面將看到un<0可以由熱力學(xué)穩(wěn)定來保障.

    松動(dòng)能的衰減可以寫成弛豫時(shí)間模型的形式

    (9),(10)和(11),(12)式中的耗散項(xiàng),可以用標(biāo)準(zhǔn)的Onsager非平衡熱力學(xué)來處理.仿照GSH的做法,有

    是非對(duì)角(off-diagonal)遷移系數(shù).如果把上面的矩陣方程寫成非矩陣形式,有f松=?η松v,即類似地,有和和將它們代入前面的(7)—(17)式,得顆粒-顆粒相互作用力為:

    彈簧變形量un,t的運(yùn)動(dòng)方程為:

    松動(dòng)能的運(yùn)動(dòng)方程為

    另外熱功率為

    將(28)—(32)與牛頓方程(2),(3)聯(lián)立求解,可計(jì)算兩個(gè)顆粒運(yùn)動(dòng)軌跡.

    以上接觸理論的材料參數(shù)是:這里除彈性勢(shì)能中的三個(gè)參數(shù)k,ξ,c必須是常數(shù)外,其他參數(shù)都是與耗散有關(guān)的遷移系數(shù),并且可以是松動(dòng)能w松、重疊量δ等變量的函數(shù).考慮到接觸力的復(fù)雜性,在這里保留了較多的遷移系數(shù),以便理論有足夠的靈活性來應(yīng)對(duì).

    一個(gè)從GSH繼承來的重要性質(zhì),是非線性彈性勢(shì)能(19)有失穩(wěn)現(xiàn)象,其對(duì)應(yīng)的恰好是庫侖屈服[15].勢(shì)能(19)的力學(xué)穩(wěn)定區(qū)域由不等式

    給出,等號(hào)是穩(wěn)定區(qū)域的邊界方程.另外在穩(wěn)定區(qū)域里總是有un<0,故(19)—(21)不會(huì)出現(xiàn)虛數(shù)(因?yàn)閷?shí)際運(yùn)動(dòng)總是發(fā)生在穩(wěn)定區(qū)域內(nèi)).利用(20),(21)式可將(35)式寫成

    失穩(wěn)意味著一旦觸及(35)或(36)式,接觸力的彈性部分將快速松動(dòng).這個(gè)現(xiàn)象可以方便地通過讓在屈服邊界急劇增大來描述.

    除失穩(wěn)外,材料參數(shù)(34)式的進(jìn)一步建模還需要顧及重疊量δ的影響,因?yàn)榻佑|力、熱功率、松動(dòng)能等都是重疊時(shí)δ>0才能被激發(fā)而出現(xiàn)的東西.一旦顆粒分離δ<0,它們的激發(fā)源都會(huì)消失,同時(shí)也迅速衰減為零(假設(shè)沒有濕黏力c=0).下面試用一組簡(jiǎn)單的材料參數(shù)值,以討論碰撞恢復(fù)系數(shù)問題為例,對(duì)此做些具體說明.

    3 碰撞恢復(fù)系數(shù)

    作為一個(gè)簡(jiǎn)單例子,取下面的參數(shù)模型

    其中Θ(x)是Heaviside階梯函數(shù),d=d1=d2和m=m1=m2分別是顆粒直徑和質(zhì)量(考慮兩個(gè)相同顆粒).(37)式中假設(shè)了沒有濕黏力,另外盡量略去了一些遷移系數(shù).模型(38)意味著彈簧變形的激發(fā)源,即(30)和(31)式右邊的第一項(xiàng)只在重疊時(shí)出現(xiàn),分離后它們是零.黏滯系數(shù)(39)式也是這樣.與GSH類似,(40)式假設(shè)切向彈簧變形的弛豫系數(shù)比例于(相當(dāng)于GSH里的顆粒溫度Tg).(41)式在滿足穩(wěn)定條件(35)時(shí)等于零,一旦不滿足穩(wěn)定條件,它將取適當(dāng)?shù)恼?起到盡快松弛切向彈簧變形、讓系統(tǒng)迅速返回穩(wěn)定區(qū)域的目的.注意上述模型的法向耗散與切向耗散有很大差異,后者比前者多考慮了(40)和(41)式兩個(gè)系數(shù).這顯然有其合理性,因?yàn)榍邢蚝纳⒌拇_比法向的復(fù)雜,需要更多的遷移系數(shù)來描寫.另外δ和un,t有長(zhǎng)度量綱,為了讓Heaviside函數(shù)的自變量是量綱1的,用顆粒粒徑d對(duì)它們做了約化.另外發(fā)生顆粒重疊和屈服的數(shù)學(xué)描寫總會(huì)涉及不連續(xù)性,需要用Heaviside函數(shù)來表達(dá).這里是將它們反映在遷移系數(shù)(38),(39)和(41)式中,并且后者與彈性勢(shì)能失穩(wěn)判據(jù)(35)式有關(guān).這個(gè)基于熱力學(xué)概念的描述方法,與直接在力模型中唯象描寫的傳統(tǒng)做法是有差異的.

    圖3是上述參數(shù)值下數(shù)值計(jì)算的兩個(gè)等質(zhì)量和直徑的、沿著x方向相對(duì)運(yùn)動(dòng)的顆粒做對(duì)心碰撞的情況.碰撞前后的速度分別是±v0和±vf(都沿著x方向),并且都是分離的,fn,t,un,t,w松都是零.顆粒重疊作用時(shí),它們的激發(fā)、衰減、以及圖2中三個(gè)能量之間的轉(zhuǎn)換情況等,都能夠明確計(jì)算(見圖3).注意由于熱功率R≥0,圖3(c)的機(jī)械能(黑色曲線)總是隨時(shí)間減小或不變.如果忽略碰撞過程細(xì)節(jié),可以將這個(gè)接觸作用模型看作以e=|vf/v0|為恢復(fù)系數(shù)的非彈性碰撞.

    在分析長(zhǎng)時(shí)間接觸的力鏈演化問題時(shí),除顆粒初始位置和速度外,還需要un,t,w松的初始值.通??梢匀〕跏挤ㄏ驈椈勺冃蝩n等于初始重疊量δ.另外任何接觸力模型(包括Hertz-Mindlin和Luding)都會(huì)出現(xiàn)Heaviside階梯函數(shù)Θ(x). 數(shù)值分析時(shí)可用近似的解析函數(shù)來代替,例如取

    本文的熱力學(xué)接觸力模型的恢復(fù)系數(shù)e=|vf/v0|不是材料常數(shù),而是初速度v0和(斜碰撞(oblique impact)時(shí)的)碰撞長(zhǎng)度的函數(shù).對(duì)心碰撞時(shí),本節(jié)考慮的模型的恢復(fù)系數(shù)隨初速度增加而減小(圖4i).這與實(shí)驗(yàn)定性符合[2],表明熱力學(xué)建模方法有描述復(fù)雜觀測(cè)細(xì)節(jié)的潛力.

    圖3 顆粒-顆粒對(duì)心碰撞過程的時(shí)間演化曲線 (a)位置坐標(biāo);(b)接觸力(c)彈簧勢(shì)能we,動(dòng)能wk和機(jī)械能w;(d)松動(dòng)能wrelaxation;是兩個(gè)顆粒的初始能量(動(dòng)能)Fig.3.Temporal evolution of normal collision of two grains:(a)Relative position;(b)contact force f=(c)potential of springs we,kinematic en-ergy wk,mechanic energy w;(d)relaxation energywrelaxation.The notation is initial(kine-matic)energy of two grains.

    本節(jié)計(jì)算的目的僅限于用一個(gè)簡(jiǎn)單例子對(duì)熱力學(xué)建模方法做些具體解釋.實(shí)際顆粒的參數(shù)情況不一定是(37)—(41)式那樣簡(jiǎn)單,參數(shù)值的標(biāo)定也需要結(jié)合實(shí)驗(yàn)仔細(xì)處理.另外出于方便解釋概念和方法,本文忽略了三維和顆粒轉(zhuǎn)動(dòng)運(yùn)動(dòng)φ(t).但在分析實(shí)際問題時(shí),它們是應(yīng)該計(jì)入的.

    圖4 對(duì)心碰撞恢復(fù)系數(shù)隨初速度變化曲線Fig.4.Variation of normal restitution coefficient with collision velicity.

    4 討 論

    為方便討論本文熱力學(xué)模型與其他接觸力模型的概念差異,我們將彈性能方程(8)給出的力,和方程(28)—(31)寫成下面的矢量形式:

    1)彈性力π總是要求為保守力.也就是說有明確定義的彈性勢(shì)能w彈(u).目前的DEM計(jì)算,雖然有些采用了保守的線性法向和線性切向彈簧,但有些是非保守的(例如Boussinesq力).也就是說對(duì)保守性的要求并不像熱力學(xué)那樣,是強(qiáng)制性的.非保守的彈簧模型會(huì)導(dǎo)致機(jī)械功率和熱功率之間的分配比例不能清楚定義的問題.

    2)彈簧應(yīng)變u并不是純幾何的位矢r=∫vdt,而是由運(yùn)動(dòng)方程(43)描述的獨(dú)立變量.這在目前的DEM計(jì)算中都注意到了(特別是切向分量),但u的運(yùn)動(dòng)方程與(43)式有很大不同.

    3)接觸力方程(44)中的黏力在當(dāng)前DEM計(jì)算中大都有類似的對(duì)應(yīng)(例如(5),(6)中的ηt,n項(xiàng)).但在對(duì)黏滯系數(shù)?η的處理方式上,沒有圖2所示的二級(jí)不可逆概念背景.

    4)本文的接觸力和幾何運(yùn)動(dòng)學(xué)速度(f,v),與彈簧力和彈簧變形(π,u)之間的關(guān)系方程(43),(44)采用了典型的Onsager非平衡熱力學(xué)形式.但這里除了熟知的機(jī)械能和熱能外,還必須考慮一個(gè)中間緩沖能量w松去控制某些遷移過程和系數(shù).這個(gè)稱作“二級(jí)不可逆”的方案是顆粒物質(zhì)熱力學(xué)的基本特征.它不僅對(duì)宏觀連續(xù)力學(xué)層面的GSH有效[11?14],細(xì)觀層面的顆粒接觸熱力學(xué)也必須這樣處理.它能把庫侖屈服、彈性松弛、力與幾何變形之間的滯迴和棘輪(ratcheting)等復(fù)雜力學(xué)現(xiàn)象,直接自然地與遷移過程聯(lián)系起來,使得用熱力學(xué)語言來理解和分析它們成為可能.從物理角度看這顯然是合理的,因?yàn)轭w粒物質(zhì)的這些復(fù)雜力學(xué)現(xiàn)象都源自于耗散發(fā)熱.

    5)方程(43),(44)右邊的第二項(xiàng)可分別理解為彈性弛豫和黏滯.右邊的第一項(xiàng)涉及到以往未曾遇到的非對(duì)角遷移系數(shù)?G或?α=?1??G.如果忽略其他遷移系數(shù),可將?α單獨(dú)看作一個(gè)具有改變兩端力和形變比值效果的“力學(xué)元件”.我們把它形象地稱作“變速箱(gearbox)”,其能量和力學(xué)關(guān)系情況如圖5所示.值得指出的是,如果變速比?G恒定(α=const.)并且u=(1?α)ε的話,可將圖中的彈簧勢(shì)能和力寫成幾何變形ε的函數(shù):w彈=(1?α)2kε2/2和f=?(1?α)2kε,它們滿足熱力學(xué)關(guān)系f=??w彈/?ε.這意味著如果不考慮α的動(dòng)態(tài)效應(yīng),變速箱是不需要的,因?yàn)樗淖饔弥皇钦{(diào)整與其聯(lián)接的彈簧的彈性系數(shù),完全可以將其歸入彈簧元件中.這也許是以往的各種力學(xué)元件都不曾提到變速箱的原因.但如果α是動(dòng)態(tài)的,在考慮具有熱力學(xué)背景的力學(xué)理論時(shí),變速箱與彈簧就必須區(qū)別開來了.變速箱這個(gè)非對(duì)角遷移系數(shù)是我們近年來為理解顆粒物質(zhì)宏觀彈塑轉(zhuǎn)變行為而提出的[17],它的統(tǒng)計(jì)物理和動(dòng)理學(xué)意義還有待深入研究.唯象地看它的力學(xué)效果是:“在不增加熱損失的前提下減弱彈性”.目前可以肯定的是,作為彈應(yīng)變u的驅(qū)動(dòng)源項(xiàng)的系數(shù)(見(43)式),它與彈性的出現(xiàn)與消亡密切相關(guān),而且是動(dòng)態(tài)的.本文為了簡(jiǎn)單,假設(shè)α隨著由重疊量控制的接觸彈性的出現(xiàn)與消亡,在0和1之間跳變,見模型(38).當(dāng)然實(shí)際情況可以復(fù)雜許多,例如還與松動(dòng)能w松有關(guān),在0和1之間變化等.

    6)GSH熱力學(xué)的一個(gè)重要特征是圖2所示的二級(jí)不可逆耗散現(xiàn)象:除了左邊的直接耗散為熱的過程外,機(jī)械能還可以通過右邊的緩沖能量w松進(jìn)行耗散,同時(shí)遷移系數(shù)也可以受其控制(即?G,?λ,?η可以是w松的函數(shù)).需要注意的是,庫侖屈服現(xiàn)象雖然源自于彈性勢(shì)能的靜力失穩(wěn)(見(35),(36)),我們還需要在某些遷移系數(shù)上采取相應(yīng)的措施來保障任何動(dòng)力學(xué)演化都不會(huì)進(jìn)入該非穩(wěn)區(qū)域,例如模型(41).(因?yàn)槿魏斡^測(cè)到的實(shí)際系統(tǒng)都在熱力學(xué)穩(wěn)定區(qū)域內(nèi)運(yùn)動(dòng)著,不會(huì)進(jìn)入非穩(wěn)區(qū)域.這與廣為熟悉的范德瓦耳斯氣體理論類似,它的熱力學(xué)狀態(tài)空間里也有一個(gè)實(shí)際氣體不可能進(jìn)入的非穩(wěn)區(qū)域).

    圖5 變速箱的能量與力學(xué)關(guān)系(注意α可以是動(dòng)態(tài)的,但m,k假設(shè)為常數(shù))Fig.5.Energetics and mechanics of gearbox.Notably α may temporally change,but m,k are assumed constant.

    由于兩個(gè)顆粒之間的很多相互作用行為是可以實(shí)驗(yàn)觀測(cè)的,任何接觸力模型原則上都應(yīng)該用這些細(xì)觀顆粒尺度的測(cè)量結(jié)果進(jìn)行驗(yàn)證和標(biāo)定.但以往的一些DEM工作發(fā)現(xiàn),這個(gè)細(xì)觀-宏觀聯(lián)系并不都是協(xié)調(diào)的,例如用線性彈簧計(jì)算的宏觀行為,卻比用更合理的赫茲彈簧要好(參見文獻(xiàn)[12]的第六章).也許這些矛盾和困難與所采用的接觸力模型對(duì)熱力學(xué)概念的描述不夠全面有關(guān).如果是這樣的話本文的熱力學(xué)建模方法將有助于問題的澄清.

    目前顆粒物理研究方法可分為純幾何、DEM接觸力、和GSH熱力學(xué)三種.幾何法致力于僅從顆粒的運(yùn)動(dòng)學(xué)變量(r,φ)出發(fā),通過深入考察其幾何時(shí)空結(jié)構(gòu)與關(guān)聯(lián),以及幾何構(gòu)型熵(con fi guration entropy)或漲落耗散關(guān)系( fl uctuation-dissipation relation)等,來解釋顆粒物質(zhì)的各種物性(參見文獻(xiàn)[18—21]).也就是說期望構(gòu)建類似開普勒行星運(yùn)動(dòng)理論那樣的顆粒物理.傳統(tǒng)DEM方法的出發(fā)點(diǎn)是{r,φ,f},即比純幾何法增加了接觸力f,并且期望在牛頓方程的基礎(chǔ)上解釋顆粒物質(zhì)物性.由于接觸力f不是可以用哈密頓描寫的保守力,DEM一般都是直接數(shù)值求解大量顆粒的牛頓方程組.如果忽略f中的不可逆,僅保留保守的彈性力,也可以用傳統(tǒng)的平衡態(tài)統(tǒng)計(jì)物理來分析這個(gè)哈密頓顆粒系統(tǒng)[22].另外對(duì)f中的不可逆部分,還可以嘗試用動(dòng)理學(xué)[23]或mode-coupling[24]等非平衡統(tǒng)計(jì)技術(shù)來分析.但隨著顆粒數(shù)密度增加,這些非平衡統(tǒng)計(jì)技術(shù)都面臨困難,目前還無法滿意地處理顆粒固體,特別是集體彈性行為的納入和描述.因此數(shù)值求解牛頓方程仍是當(dāng)前分析接觸作用系統(tǒng)的最好辦法(如引言里提到的,它能描述宏觀彈性的呈現(xiàn)于消亡).注意幾何法和傳統(tǒng)DEM法都可以不具體涉及熱力學(xué)熵,可稱作“非熱(athermal)”類型理論.與之對(duì)比,熱力學(xué)方法則強(qiáng)制性地要求無論是連續(xù)力學(xué)的宏觀尺度,還是本文考慮的顆粒接觸力細(xì)觀尺度,熱力學(xué)能量守恒和熵增加原理都應(yīng)該得到明確的體現(xiàn)和遵守.當(dāng)然三種方法之間有共同相似或互補(bǔ)的內(nèi)容,但肯定還有不能互通的內(nèi)容,因?yàn)榕nD力不可能取代熱力學(xué),也不可能被純幾何學(xué)取代.考慮到耗散發(fā)熱是顆粒物質(zhì)的基本性質(zhì),物理熱力學(xué)應(yīng)該是最為合理和可靠的方法.

    值得提到的是,固體顆粒間的法向接觸力一般沒有塑性打滑等現(xiàn)象,模型(37)中取很可能是普遍合理的.但在Luding模型里法向滯迴是他的特色[5].為方便今后與其對(duì)比,我們?cè)贠nsager矩陣(23)—(26)里保留了這幾個(gè)遷移系數(shù).遷移系數(shù)那些是無價(jià)值(irrelevant)可以去掉的,那些是有價(jià)值(relevant)需要保留的,及其相關(guān)的動(dòng)力學(xué)意義等,目前還不能確定.這些問題需要結(jié)合DEM計(jì)算和實(shí)驗(yàn)情況來逐步澄清.模型(37)—(41)的參數(shù)數(shù)值僅作為參考,沒有普適性.它們都是材料參數(shù),應(yīng)該依據(jù)實(shí)驗(yàn)來標(biāo)定.另外將本文理論向三維推廣時(shí),只需要對(duì)文中的力、速度、位移等矢量補(bǔ)充一個(gè)切向分量即可.對(duì)球形顆粒,由于對(duì)稱性兩個(gè)切向的遷移系數(shù)是一樣的,理論參數(shù)仍由(34)式列出.如果考慮(37)—(41)式的模型,這個(gè)三維理論給出的二維運(yùn)動(dòng)將退化為本文的公式.這時(shí)它給出的對(duì)心碰撞也就是本文的圖3和圖4.顆粒轉(zhuǎn)動(dòng)在一些接觸力模型有所考慮,例如Luding模型[5].但也有不少DEM工作省略了這個(gè)復(fù)雜性,例如文獻(xiàn)[9].本文出于簡(jiǎn)單和方便討論熱力學(xué)概念的原因也省略了顆粒轉(zhuǎn)動(dòng).它的計(jì)入要復(fù)雜許多,建議參照Luding的做法.

    常用接觸力模型與本文建議的熱力學(xué)方法有很多對(duì)應(yīng)的和相似的地方,例如彈簧和黏滯阻尼.算例模型(37)—(41)式中的彈簧系數(shù)k(對(duì)應(yīng)于剛度系數(shù)),也是常用模型中經(jīng)常出現(xiàn)的常數(shù).但熱力學(xué)肯定還有自己的不同內(nèi)容和改善的地方,例如將屈服和打滑放到了遷移系數(shù)中、給出了機(jī)械能和熱功率的具體表達(dá)式、要求伴隨任何遷移耗散過程的熱功率都有“二次正定多項(xiàng)式”的形式(見(32),(33)式)等.注意“二次正定多項(xiàng)式”的要求來自熱力學(xué)穩(wěn)定性,它也會(huì)改善模型的數(shù)值穩(wěn)定性.常用模型一般只有黏性的耗散熱是平方正定的形式,其他復(fù)雜的不可逆現(xiàn)象往往沒有這個(gè)性質(zhì),例如屈服.顯然熱力學(xué)內(nèi)容的最大優(yōu)勢(shì)將體現(xiàn)在研究能量和耗散問題方面,這在最近已經(jīng)開始受到重視(例如Yu小組的工作[25,26]).

    本文建模方法與常用接觸力模型的仔細(xì)對(duì)比,特別是具體改善和優(yōu)越的地方,將在今后的定量對(duì)比工作中具體體現(xiàn)出來.這里可分為物理和工程兩個(gè)方面來開展.物理方面,僅關(guān)心“二體力”這個(gè)基礎(chǔ)層面即可,也就是說只是對(duì)比不同方法和模型的“兩個(gè)顆粒動(dòng)力學(xué)”情況即可,這樣便于討論和澄清與基本概念相關(guān)的優(yōu)劣對(duì)比,同時(shí)避開大規(guī)模數(shù)值編程的麻煩.工程方面,主要考慮大量顆粒的DEM計(jì)算情況的對(duì)比,涉及大規(guī)模數(shù)值編程或商業(yè)軟件.另外“二體力”和“DEM”兩邊都應(yīng)該盡可能多地考慮與實(shí)驗(yàn)的對(duì)比.

    在考慮計(jì)入顆粒轉(zhuǎn)動(dòng)因素時(shí),首先是增加描述顆粒轉(zhuǎn)動(dòng)角的運(yùn)動(dòng)方程(即牛頓力矩方程),在機(jī)械能公式(16)中增加顆粒轉(zhuǎn)動(dòng)動(dòng)能,以及可能的扭力勢(shì)能等.然后是滾動(dòng)力(rolling)和扭轉(zhuǎn)力(torsion)的建模問題.由于熱力學(xué)的影響主要體現(xiàn)在不可逆耗散部分,建模時(shí)涉及的純幾何運(yùn)動(dòng)學(xué)(kinematic)部分,以及滾動(dòng)力和扭轉(zhuǎn)力的可逆部分,基本可直接挪用Luding的相應(yīng)公式[5].但滾動(dòng)力和扭轉(zhuǎn)力的不可逆部分,需要用與本文類似的能量守恒和Onsager遷移系數(shù)等方法來處理.

    5 結(jié) 論

    雖然耗散發(fā)熱是顆粒間接觸運(yùn)動(dòng)時(shí)的重要現(xiàn)象,當(dāng)前DEM采用的接觸力模型仍屬于唯象力學(xué)類型,對(duì)熱力學(xué)要求的機(jī)械功率和熱功率等的定義并不十分清楚.本文的主要結(jié)論是:基于熱力學(xué)的接觸力建模方法可以參照已有的GSH概念來實(shí)現(xiàn).我們相信從唯象力學(xué)發(fā)展到熱力學(xué)模型,不僅僅從物理概念上更加合理,對(duì)諸如恢復(fù)系數(shù)行為等的一些實(shí)驗(yàn)現(xiàn)象的描述,也會(huì)有所改善.當(dāng)然在嘗試將本文內(nèi)容用于DEM數(shù)值計(jì)算之前,還需要先將其推廣到三維和計(jì)入顆粒轉(zhuǎn)動(dòng)運(yùn)動(dòng).這些將是今后需要陸續(xù)開展的工作.

    感謝厚美瑛、孫其誠(chéng)、程曉輝、劉曉星和張國(guó)華的討論.

    [1]Hertz H 1881J.Reine Angew.Math.92 156

    [2]Thornton C 2015Granular Dynamics,Contact Mechanics and Particle System Simulations—A DEM Study(Particle Technology Series,Volume 24)(eBook DOI 10.1007/978-3-319-18711-2)(Switzerland:Springer International Publishing AG)

    [3]Laughlin R B,Pines D 2000Proc.Natl.Acad.Sci.USA97 28

    [4]Laughlin R B 2004A Different Universe(Changsha:Hunan Science and Technology Press)(in Chinese)[王文浩譯2008不同的宇宙(長(zhǎng)沙:湖南科學(xué)技術(shù)出版社)]

    [5]Luding S 2008Granular Matter10 235

    [6]Truesdell C 1972Rational Thermodynamics(Berlin:Springer-Verlag)

    [7]Boussinesq J 1873C.R.Hebd.Seances Acad.Sci.77 1521

    [8]Bonneau L,Andreotti B,Clément E 2007Phys.Rev.E75 016602

    [9]Sun Q C,Jin F,Zhou G D 2013Granular Matter15 119

    [10]Kumar N 2014Ph.D.Dissertation(Enschede:University of Twente.The Netherlands ISBN:978-90-365-3634-9,DOI:10.3990/1.9789036536349)

    [11]Jiang Y M,Liu M 2009Granular Matter11 139

    [12]Sun Q C,Hou M Y,Jin F 2011Physics and Mechanics of Granular Matter(Beijing:Science Press)(in Chinese)[孫其誠(chéng),厚美瑛,金峰 2011顆粒物質(zhì)物理與力學(xué)(北京:科學(xué)出版社)]

    [13]Jiang Y M,Liu M 2014Acta Mech.225 2363

    [14]Jiang Y M,Liu M 2015Eur.Phys.J.E38 15

    [15]Jiang Y M,Liu M 2003Phys.Rev.Lett.91 144301

    [16]Jiang Y M,Liu M 2004Phys.Rev.Lett.93 148001

    [17]Jiang Y M,Liu M 2007Phys.Rev.Lett.99 105501

    [18]Torquato S,Stillinger F H 2010Rev.Mod.Phys.82 2633

    [19]Edwards S F,Moun field C C 1996Physica A226 1

    [20]Edwards S F,Moun field C C 1996Physica A226 12

    [21]Edwards S F,Moun field C C 1996Physica A226 25

    [22]Parisi G,Zamponi F 2010Rev.Mod.Phys.82 789

    [23]Luding S 2009Nonlinearity22 R101

    [24]Hayakawa H,Otsuki M 2008Prog.Theor.Phys.119 381

    [25]Kuang S B,Zou R P,Pan R H,Yu A B 2012Ind.Eng.Chem.Res.51 14289

    [26]Hou Q F,Kuang S B,Yu A B 2017Chem.Engineer.Sci.161 67

    猜你喜歡
    勢(shì)能熱力學(xué)力學(xué)
    力學(xué)
    “動(dòng)能和勢(shì)能”知識(shí)鞏固
    作 品:景觀設(shè)計(jì)
    ——《勢(shì)能》
    文化縱橫(2022年3期)2022-09-07 11:43:18
    “動(dòng)能和勢(shì)能”知識(shí)鞏固
    “動(dòng)能和勢(shì)能”隨堂練
    弟子規(guī)·余力學(xué)文(十)
    快樂語文(2021年11期)2021-07-20 07:41:32
    弟子規(guī)·余力學(xué)文(四)
    快樂語文(2020年30期)2021-01-14 01:05:28
    Fe-C-Mn-Si-Cr的馬氏體開始轉(zhuǎn)變點(diǎn)的熱力學(xué)計(jì)算
    上海金屬(2016年1期)2016-11-23 05:17:24
    活塞的靜力學(xué)與熱力學(xué)仿真分析
    電子制作(2016年19期)2016-08-24 07:49:54
    力學(xué) 等
    亚洲av.av天堂| 秋霞在线观看毛片| 美女高潮的动态| 内地一区二区视频在线| 亚洲国产精品999| 99精国产麻豆久久婷婷| 亚洲精品成人久久久久久| av又黄又爽大尺度在线免费看| 日韩成人av中文字幕在线观看| 久久国内精品自在自线图片| 亚洲真实伦在线观看| 菩萨蛮人人尽说江南好唐韦庄| 丰满乱子伦码专区| 亚洲国产最新在线播放| 久热久热在线精品观看| 97超碰精品成人国产| 久久久久久久久久人人人人人人| 看非洲黑人一级黄片| 2018国产大陆天天弄谢| 69人妻影院| 欧美日韩国产mv在线观看视频 | 又爽又黄无遮挡网站| 天堂网av新在线| 亚洲欧美日韩东京热| 纵有疾风起免费观看全集完整版| 亚洲怡红院男人天堂| 亚洲精品国产成人久久av| 国产免费视频播放在线视频| 免费看av在线观看网站| 亚洲无线观看免费| 亚洲av国产av综合av卡| 国产国拍精品亚洲av在线观看| 一二三四中文在线观看免费高清| 久久久久性生活片| 波多野结衣巨乳人妻| 欧美激情久久久久久爽电影| 少妇被粗大猛烈的视频| av国产久精品久网站免费入址| 女人被狂操c到高潮| 成人综合一区亚洲| 波野结衣二区三区在线| 国产亚洲5aaaaa淫片| 各种免费的搞黄视频| 国产黄片视频在线免费观看| 国产69精品久久久久777片| 亚洲在线观看片| 国产精品一区二区三区四区免费观看| 日韩一本色道免费dvd| 视频区图区小说| 舔av片在线| 亚洲最大成人手机在线| 久久久久久久国产电影| 免费大片黄手机在线观看| 精品人妻偷拍中文字幕| 国产男女内射视频| 成年人午夜在线观看视频| 狠狠精品人妻久久久久久综合| 另类亚洲欧美激情| 中文字幕av成人在线电影| 一级毛片 在线播放| av国产精品久久久久影院| 国产黄色视频一区二区在线观看| 成年女人在线观看亚洲视频 | 亚州av有码| av在线播放精品| 一级黄片播放器| tube8黄色片| 亚洲成色77777| 久久99蜜桃精品久久| 2018国产大陆天天弄谢| 麻豆成人午夜福利视频| 国产成人a区在线观看| 秋霞伦理黄片| 亚洲国产高清在线一区二区三| 美女国产视频在线观看| 在线观看免费高清a一片| 最近2019中文字幕mv第一页| 久久久a久久爽久久v久久| av在线老鸭窝| 三级经典国产精品| 街头女战士在线观看网站| 我的女老师完整版在线观看| 久久韩国三级中文字幕| a级毛色黄片| 一级毛片aaaaaa免费看小| 精品国产一区二区三区久久久樱花 | 中文字幕免费在线视频6| 国产精品不卡视频一区二区| 亚洲av欧美aⅴ国产| 亚洲av男天堂| 亚洲精品久久午夜乱码| 国产69精品久久久久777片| 欧美潮喷喷水| 欧美日韩在线观看h| 久久久精品94久久精品| 狂野欧美激情性bbbbbb| 国产在线男女| 日韩一本色道免费dvd| av又黄又爽大尺度在线免费看| 欧美xxxx黑人xx丫x性爽| 中文字幕亚洲精品专区| 欧美日本视频| 国产真实伦视频高清在线观看| 免费播放大片免费观看视频在线观看| 好男人视频免费观看在线| 亚洲成人一二三区av| 97超碰精品成人国产| 国产永久视频网站| 高清av免费在线| 日韩免费高清中文字幕av| 欧美成人精品欧美一级黄| 国产精品人妻久久久久久| 欧美日韩视频精品一区| 国语对白做爰xxxⅹ性视频网站| 99热这里只有是精品50| 久久精品人妻少妇| 亚洲色图综合在线观看| 极品少妇高潮喷水抽搐| 亚洲激情五月婷婷啪啪| 国精品久久久久久国模美| 亚洲美女搞黄在线观看| 两个人的视频大全免费| 日韩成人av中文字幕在线观看| 内地一区二区视频在线| 亚洲第一区二区三区不卡| 国产男女内射视频| 亚洲成人一二三区av| 网址你懂的国产日韩在线| 国产成人免费观看mmmm| 少妇 在线观看| 国产一区有黄有色的免费视频| 亚洲最大成人中文| freevideosex欧美| 成人国产麻豆网| 免费大片18禁| 男人狂女人下面高潮的视频| 看十八女毛片水多多多| 久久久久久伊人网av| 六月丁香七月| 免费黄频网站在线观看国产| 边亲边吃奶的免费视频| 精品久久国产蜜桃| 国产精品99久久99久久久不卡 | 成人高潮视频无遮挡免费网站| 日韩欧美精品v在线| 听说在线观看完整版免费高清| 少妇 在线观看| 成人午夜精彩视频在线观看| 午夜福利高清视频| 欧美性猛交╳xxx乱大交人| 欧美3d第一页| 看十八女毛片水多多多| 亚洲一级一片aⅴ在线观看| 成人无遮挡网站| 日日啪夜夜爽| 国产精品久久久久久久电影| 国产国拍精品亚洲av在线观看| 成年人午夜在线观看视频| 精品一区二区三卡| 国产亚洲午夜精品一区二区久久 | 插阴视频在线观看视频| 如何舔出高潮| 久久久久国产网址| 国产精品秋霞免费鲁丝片| 我要看日韩黄色一级片| 建设人人有责人人尽责人人享有的 | 涩涩av久久男人的天堂| 色综合色国产| 亚洲国产欧美人成| 久久精品国产鲁丝片午夜精品| av在线蜜桃| 成人漫画全彩无遮挡| 女人被狂操c到高潮| 精品国产一区二区三区久久久樱花 | 久久99热这里只有精品18| 免费观看无遮挡的男女| 成年av动漫网址| 中文资源天堂在线| 亚洲天堂av无毛| av又黄又爽大尺度在线免费看| 亚洲va在线va天堂va国产| 一级爰片在线观看| 日本免费在线观看一区| 午夜日本视频在线| 在线观看国产h片| 国语对白做爰xxxⅹ性视频网站| 欧美成人午夜免费资源| 老司机影院毛片| 国产亚洲5aaaaa淫片| 黄色怎么调成土黄色| 色视频www国产| 97超碰精品成人国产| 国产午夜福利久久久久久| 免费电影在线观看免费观看| 午夜福利视频精品| 国产大屁股一区二区在线视频| 精品国产一区二区三区久久久樱花 | 老司机影院成人| 国产淫语在线视频| 亚洲精品aⅴ在线观看| 免费观看a级毛片全部| 69av精品久久久久久| 婷婷色综合www| 精品人妻一区二区三区麻豆| av天堂中文字幕网| 亚洲伊人久久精品综合| 卡戴珊不雅视频在线播放| 欧美日韩亚洲高清精品| 精品国产一区二区三区久久久樱花 | 欧美成人a在线观看| 欧美高清性xxxxhd video| 亚洲经典国产精华液单| 国产精品成人在线| 欧美性猛交╳xxx乱大交人| 97热精品久久久久久| 国产大屁股一区二区在线视频| 日韩在线高清观看一区二区三区| 亚洲成人中文字幕在线播放| 免费观看a级毛片全部| xxx大片免费视频| 国产人妻一区二区三区在| 人妻 亚洲 视频| 日产精品乱码卡一卡2卡三| 午夜免费鲁丝| 中文字幕免费在线视频6| 久久久久久久久久人人人人人人| 亚洲精品色激情综合| 亚洲国产精品成人久久小说| 女的被弄到高潮叫床怎么办| 久久久久久久久久成人| 亚洲欧美成人综合另类久久久| av卡一久久| av在线app专区| 美女脱内裤让男人舔精品视频| 菩萨蛮人人尽说江南好唐韦庄| 日本一二三区视频观看| 久久久精品免费免费高清| av在线观看视频网站免费| 视频中文字幕在线观看| 国产乱人偷精品视频| 国产午夜福利久久久久久| 欧美另类一区| 亚洲av不卡在线观看| 老司机影院毛片| 在线观看三级黄色| 欧美3d第一页| 国产av码专区亚洲av| 一边亲一边摸免费视频| 精品99又大又爽又粗少妇毛片| 欧美区成人在线视频| 日日撸夜夜添| 少妇熟女欧美另类| 亚洲经典国产精华液单| 男女国产视频网站| 如何舔出高潮| 人妻夜夜爽99麻豆av| 在线观看国产h片| 韩国高清视频一区二区三区| 大陆偷拍与自拍| 日韩制服骚丝袜av| 亚洲国产最新在线播放| 国产成人精品福利久久| 亚洲av男天堂| 在线免费观看不下载黄p国产| a级一级毛片免费在线观看| 成人毛片a级毛片在线播放| 久久久久久国产a免费观看| 亚洲成人一二三区av| 白带黄色成豆腐渣| 国产精品人妻久久久影院| 丝袜脚勾引网站| 蜜桃亚洲精品一区二区三区| 亚洲,欧美,日韩| 免费播放大片免费观看视频在线观看| 国内精品宾馆在线| 欧美日韩综合久久久久久| 久久久久久久精品精品| 极品教师在线视频| 国产免费一级a男人的天堂| 女人十人毛片免费观看3o分钟| 小蜜桃在线观看免费完整版高清| 自拍欧美九色日韩亚洲蝌蚪91 | 精品一区二区三区视频在线| 三级男女做爰猛烈吃奶摸视频| 国产综合精华液| a级一级毛片免费在线观看| 中国三级夫妇交换| 国产亚洲av片在线观看秒播厂| 亚洲成人中文字幕在线播放| 国产精品久久久久久精品古装| 国产高潮美女av| 精品久久久精品久久久| 三级经典国产精品| 亚洲国产精品成人综合色| 中文资源天堂在线| 亚洲欧洲日产国产| av专区在线播放| 18禁裸乳无遮挡动漫免费视频 | 大码成人一级视频| 亚洲国产精品国产精品| 日韩不卡一区二区三区视频在线| 嘟嘟电影网在线观看| 国产极品天堂在线| 97精品久久久久久久久久精品| 日本-黄色视频高清免费观看| 丝袜美腿在线中文| 黄色配什么色好看| av国产精品久久久久影院| 在线播放无遮挡| 色哟哟·www| 观看免费一级毛片| 国产成人免费观看mmmm| 成年av动漫网址| 成人黄色视频免费在线看| 在线免费观看不下载黄p国产| 人妻夜夜爽99麻豆av| 国产在视频线精品| 欧美三级亚洲精品| 美女cb高潮喷水在线观看| 99久久九九国产精品国产免费| 国产精品精品国产色婷婷| 亚洲av中文字字幕乱码综合| 亚洲国产欧美人成| 欧美变态另类bdsm刘玥| 久久久久久国产a免费观看| 精品一区二区三卡| 日韩av免费高清视频| 老女人水多毛片| 视频区图区小说| 极品教师在线视频| 日韩伦理黄色片| 国产精品久久久久久精品古装| 欧美亚洲 丝袜 人妻 在线| 久久久久精品性色| 男的添女的下面高潮视频| 男人和女人高潮做爰伦理| 国产精品久久久久久久电影| 亚洲欧美清纯卡通| 国产午夜福利久久久久久| 国语对白做爰xxxⅹ性视频网站| 久久久久久久久久久免费av| 精品久久久久久久末码| 国产黄片美女视频| 丰满少妇做爰视频| av.在线天堂| 一区二区三区免费毛片| 精品酒店卫生间| 日韩成人av中文字幕在线观看| 国产毛片在线视频| 欧美性感艳星| 不卡视频在线观看欧美| 嫩草影院新地址| 色5月婷婷丁香| 久久影院123| 国产精品人妻久久久久久| 丰满乱子伦码专区| 国产伦理片在线播放av一区| 午夜免费鲁丝| 99久久九九国产精品国产免费| 国产av国产精品国产| 久久鲁丝午夜福利片| 久久久久久久精品精品| 精品酒店卫生间| 国产成人免费无遮挡视频| 精品酒店卫生间| 丝袜美腿在线中文| 国产高清不卡午夜福利| 一本久久精品| 亚洲电影在线观看av| 18禁裸乳无遮挡免费网站照片| 99精国产麻豆久久婷婷| 日本欧美国产在线视频| 1000部很黄的大片| 狂野欧美白嫩少妇大欣赏| 日韩人妻高清精品专区| 国产精品一二三区在线看| 国产精品偷伦视频观看了| 99热这里只有是精品在线观看| 亚洲精品亚洲一区二区| 日本黄色片子视频| 少妇被粗大猛烈的视频| 熟妇人妻不卡中文字幕| 精品一区二区免费观看| 中文在线观看免费www的网站| 七月丁香在线播放| 国产色婷婷99| 亚洲内射少妇av| 国产欧美日韩精品一区二区| 激情 狠狠 欧美| 国产精品不卡视频一区二区| 最近中文字幕高清免费大全6| 国产毛片在线视频| 亚洲aⅴ乱码一区二区在线播放| 国产精品.久久久| 男男h啪啪无遮挡| 亚洲成人av在线免费| 18禁裸乳无遮挡免费网站照片| 肉色欧美久久久久久久蜜桃 | 久久99蜜桃精品久久| 久久精品久久久久久久性| 日日啪夜夜撸| 欧美成人精品欧美一级黄| 日韩精品有码人妻一区| 丝袜喷水一区| 久久久国产一区二区| 亚洲精品一区蜜桃| 亚洲国产av新网站| tube8黄色片| 国产av不卡久久| 一级毛片aaaaaa免费看小| 久久久久久国产a免费观看| 亚洲精华国产精华液的使用体验| 亚洲av一区综合| 国产探花极品一区二区| 一二三四中文在线观看免费高清| 国产av码专区亚洲av| 亚洲av免费高清在线观看| av在线蜜桃| 少妇 在线观看| 黄片无遮挡物在线观看| 亚洲人成网站在线播| 免费av毛片视频| 午夜精品国产一区二区电影 | 热99国产精品久久久久久7| 在线天堂最新版资源| 在线 av 中文字幕| 亚洲真实伦在线观看| 美女视频免费永久观看网站| 亚洲av成人精品一二三区| 老师上课跳d突然被开到最大视频| xxx大片免费视频| av国产免费在线观看| 久久精品国产亚洲av天美| 丰满人妻一区二区三区视频av| 日本一本二区三区精品| 麻豆国产97在线/欧美| 欧美日本视频| 久久韩国三级中文字幕| 午夜福利在线观看免费完整高清在| 国产人妻一区二区三区在| 欧美亚洲 丝袜 人妻 在线| 欧美 日韩 精品 国产| 一区二区av电影网| 男女那种视频在线观看| 国产精品三级大全| 精品少妇黑人巨大在线播放| 能在线免费看毛片的网站| 我的老师免费观看完整版| 成人漫画全彩无遮挡| 伊人久久精品亚洲午夜| 国语对白做爰xxxⅹ性视频网站| 日韩欧美 国产精品| 午夜激情福利司机影院| 99热国产这里只有精品6| av在线app专区| 日韩av免费高清视频| 国产毛片在线视频| 国产免费一区二区三区四区乱码| 91久久精品电影网| 欧美xxxx性猛交bbbb| 国产午夜精品久久久久久一区二区三区| 中文字幕免费在线视频6| 少妇人妻久久综合中文| 秋霞伦理黄片| 免费观看性生交大片5| 蜜臀久久99精品久久宅男| 天堂俺去俺来也www色官网| 别揉我奶头 嗯啊视频| 一本色道久久久久久精品综合| 久久久精品免费免费高清| 午夜免费鲁丝| 亚洲国产高清在线一区二区三| 最近的中文字幕免费完整| 欧美日韩视频精品一区| 观看免费一级毛片| 在线观看国产h片| 看十八女毛片水多多多| 麻豆成人午夜福利视频| 久久国产乱子免费精品| 久久精品国产亚洲av天美| 一级片'在线观看视频| 午夜免费男女啪啪视频观看| 色视频在线一区二区三区| 精品国产一区二区三区久久久樱花 | av卡一久久| av在线蜜桃| 国产69精品久久久久777片| 菩萨蛮人人尽说江南好唐韦庄| 日韩电影二区| av国产精品久久久久影院| 日本免费在线观看一区| 天堂俺去俺来也www色官网| 亚洲怡红院男人天堂| 国产大屁股一区二区在线视频| 高清av免费在线| 一个人看视频在线观看www免费| videossex国产| 成人毛片60女人毛片免费| 三级国产精品片| 九九久久精品国产亚洲av麻豆| 成人一区二区视频在线观看| 精品少妇久久久久久888优播| 一级毛片黄色毛片免费观看视频| 国产精品国产三级国产专区5o| 亚洲av成人精品一二三区| 免费人成在线观看视频色| 久久综合国产亚洲精品| 一级毛片黄色毛片免费观看视频| 国产精品一二三区在线看| 精品久久久久久电影网| 99久久精品一区二区三区| 69人妻影院| 国产免费福利视频在线观看| 免费看日本二区| 日韩大片免费观看网站| 亚洲国产精品999| 最后的刺客免费高清国语| 国产永久视频网站| 99热这里只有精品一区| 97在线视频观看| 久久久久久久久久久免费av| 国产黄色视频一区二区在线观看| 国内揄拍国产精品人妻在线| 日韩欧美一区视频在线观看 | 深爱激情五月婷婷| 精品国产露脸久久av麻豆| 免费在线观看成人毛片| www.色视频.com| 欧美高清成人免费视频www| 一区二区三区免费毛片| 2022亚洲国产成人精品| 涩涩av久久男人的天堂| 偷拍熟女少妇极品色| 国产精品久久久久久久电影| 热re99久久精品国产66热6| 国产精品国产三级国产av玫瑰| 尾随美女入室| 18禁裸乳无遮挡免费网站照片| 国产亚洲5aaaaa淫片| 国产日韩欧美亚洲二区| 日韩强制内射视频| 人体艺术视频欧美日本| 性色av一级| 国产欧美另类精品又又久久亚洲欧美| 搞女人的毛片| 国产精品嫩草影院av在线观看| 亚洲成人中文字幕在线播放| 久久精品熟女亚洲av麻豆精品| 一区二区三区四区激情视频| 国产高潮美女av| 亚洲久久久久久中文字幕| 久久精品综合一区二区三区| 亚洲色图av天堂| 老司机影院成人| 一区二区三区精品91| 国产综合懂色| 视频中文字幕在线观看| 午夜福利网站1000一区二区三区| 极品教师在线视频| 亚洲欧洲日产国产| 在线a可以看的网站| 国产亚洲一区二区精品| 在线观看三级黄色| 色综合色国产| 最近2019中文字幕mv第一页| 免费av毛片视频| 国产伦在线观看视频一区| 国产在视频线精品| 亚洲成人中文字幕在线播放| 国产成人精品一,二区| 国产一区二区亚洲精品在线观看| www.av在线官网国产| 简卡轻食公司| www.色视频.com| 成人午夜精彩视频在线观看| 欧美成人一区二区免费高清观看| 国产极品天堂在线| 国产成人精品一,二区| 免费观看无遮挡的男女| 午夜日本视频在线| 亚洲天堂国产精品一区在线| 亚洲精品久久午夜乱码| 精品人妻偷拍中文字幕| 国产日韩欧美亚洲二区| 日本爱情动作片www.在线观看| 国产免费一级a男人的天堂| 高清毛片免费看| 久久ye,这里只有精品| 插逼视频在线观看| 极品教师在线视频| 精品午夜福利在线看| 国产成人午夜福利电影在线观看| 最近手机中文字幕大全| 五月伊人婷婷丁香| 成人综合一区亚洲| 亚洲国产精品成人综合色| 国产极品天堂在线| 国产色婷婷99| 好男人视频免费观看在线| 亚洲成人精品中文字幕电影| 欧美潮喷喷水| 国语对白做爰xxxⅹ性视频网站| 哪个播放器可以免费观看大片| 精品久久久久久久人妻蜜臀av| 久久久久九九精品影院| 五月伊人婷婷丁香| 久久久欧美国产精品| 美女cb高潮喷水在线观看| 男女下面进入的视频免费午夜| 国产午夜精品一二区理论片| 99久久精品国产国产毛片| 视频区图区小说| 国产成人一区二区在线| 在线 av 中文字幕| 亚洲综合色惰|