• <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免费在线| 首页视频小说图片口味搜索| 一区二区av电影网| 久久精品人人爽人人爽视色| 精品一品国产午夜福利视频| 日韩免费av在线播放| 欧美亚洲日本最大视频资源| 久久国产精品影院| 99国产精品99久久久久| 狂野欧美激情性xxxx| 91国产中文字幕| 18禁黄网站禁片午夜丰满| 免费观看a级毛片全部| 久久性视频一级片| 亚洲精华国产精华精| 国产一区二区三区在线臀色熟女 | 欧美老熟妇乱子伦牲交| 亚洲精品在线观看二区| 伊人久久大香线蕉亚洲五| 一区二区三区激情视频| 欧美日韩av久久| av不卡在线播放| 啦啦啦在线免费观看视频4| 午夜免费鲁丝| 伦理电影免费视频| 少妇的丰满在线观看| 亚洲人成电影免费在线| 成人手机av| 日本a在线网址| 亚洲久久久国产精品| 一个人免费看片子| 成人av一区二区三区在线看| 中文字幕精品免费在线观看视频| 脱女人内裤的视频| 一级毛片女人18水好多| 亚洲第一欧美日韩一区二区三区 | 日本撒尿小便嘘嘘汇集6| 久久精品成人免费网站| 日韩视频一区二区在线观看| 精品人妻1区二区| 两性午夜刺激爽爽歪歪视频在线观看 | 亚洲专区字幕在线| 国产成人免费无遮挡视频| svipshipincom国产片| 成人国语在线视频| 波多野结衣一区麻豆| 亚洲欧美一区二区三区黑人| 日韩熟女老妇一区二区性免费视频| 国产av一区二区精品久久| 亚洲伊人色综图| 嫩草影视91久久| 久久久久久久大尺度免费视频| 在线十欧美十亚洲十日本专区| 午夜免费鲁丝| 午夜免费鲁丝| 脱女人内裤的视频| 成人永久免费在线观看视频 | 久久精品熟女亚洲av麻豆精品| 亚洲精品在线观看二区| 自线自在国产av| 黄色 视频免费看| 国产亚洲一区二区精品| 久久久久网色| 在线十欧美十亚洲十日本专区| 99九九在线精品视频| 亚洲色图综合在线观看| √禁漫天堂资源中文www| 久久久国产一区二区| av网站在线播放免费| 美女扒开内裤让男人捅视频| 欧美大码av| 黄片播放在线免费| 国产一区二区 视频在线| √禁漫天堂资源中文www| 黄色丝袜av网址大全| 国产人伦9x9x在线观看| 窝窝影院91人妻| 50天的宝宝边吃奶边哭怎么回事| 精品一区二区三区av网在线观看 | 91九色精品人成在线观看| 国产成人免费无遮挡视频| 99国产精品一区二区三区| 99国产精品一区二区三区| 亚洲色图av天堂| 777久久人妻少妇嫩草av网站| 国产精品亚洲一级av第二区| 巨乳人妻的诱惑在线观看| 亚洲美女黄片视频| 国产在线观看jvid| 99re6热这里在线精品视频| 欧美变态另类bdsm刘玥| 国产在线观看jvid| 久久久精品免费免费高清| 亚洲美女黄片视频| 日韩 欧美 亚洲 中文字幕| 波多野结衣av一区二区av| 亚洲全国av大片| 婷婷成人精品国产| 久久久久久免费高清国产稀缺| 亚洲午夜精品一区,二区,三区| av电影中文网址| 午夜日韩欧美国产| 免费一级毛片在线播放高清视频 | 法律面前人人平等表现在哪些方面| 大码成人一级视频| 亚洲免费av在线视频| 精品国产一区二区久久| 国产黄频视频在线观看| 日韩人妻精品一区2区三区| 最近最新中文字幕大全免费视频| 狠狠婷婷综合久久久久久88av| 久久精品91无色码中文字幕| 国产精品亚洲av一区麻豆| 淫妇啪啪啪对白视频| 久久午夜综合久久蜜桃| 欧美乱码精品一区二区三区| 精品少妇内射三级| 高清av免费在线| 制服人妻中文乱码| 色播在线永久视频| 一本久久精品| 亚洲色图av天堂| 伊人久久大香线蕉亚洲五| 巨乳人妻的诱惑在线观看| 欧美av亚洲av综合av国产av| av在线播放免费不卡| 99re在线观看精品视频| 国产成人系列免费观看| 美女午夜性视频免费| 久久狼人影院| 久久天堂一区二区三区四区| 久久中文看片网| 日韩熟女老妇一区二区性免费视频| 久久久国产精品麻豆| 国产免费现黄频在线看| 成人亚洲精品一区在线观看| 亚洲人成77777在线视频| 性高湖久久久久久久久免费观看| 亚洲精品国产精品久久久不卡| 国产精品九九99| cao死你这个sao货| 麻豆av在线久日| 在线观看免费日韩欧美大片| 国产色视频综合| 操出白浆在线播放| 一进一出抽搐动态| e午夜精品久久久久久久| 亚洲第一欧美日韩一区二区三区 | 久久人妻av系列| 青青草视频在线视频观看| 在线观看免费高清a一片| 中文字幕高清在线视频| 国产精品一区二区精品视频观看| 国产亚洲一区二区精品| 人人妻,人人澡人人爽秒播| 十八禁高潮呻吟视频| 亚洲成人免费av在线播放| 巨乳人妻的诱惑在线观看| 一区福利在线观看| 1024香蕉在线观看| 丝袜美腿诱惑在线| av超薄肉色丝袜交足视频| 最新的欧美精品一区二区| 女警被强在线播放| 可以免费在线观看a视频的电影网站| 欧美大码av| 咕卡用的链子| 少妇 在线观看| 男人舔女人的私密视频| 宅男免费午夜| 免费女性裸体啪啪无遮挡网站| 汤姆久久久久久久影院中文字幕| 日韩中文字幕视频在线看片| 天堂8中文在线网| 精品久久久精品久久久| 亚洲专区国产一区二区| 不卡一级毛片| 免费观看a级毛片全部| 9色porny在线观看| 国产日韩一区二区三区精品不卡| 精品高清国产在线一区| 一区二区三区乱码不卡18| 国产精品免费视频内射| 亚洲国产成人一精品久久久| 丁香欧美五月| 亚洲,欧美精品.| 俄罗斯特黄特色一大片| 亚洲精品自拍成人| 国产精品秋霞免费鲁丝片| 自拍欧美九色日韩亚洲蝌蚪91| 国产成人av激情在线播放| 欧美日韩国产mv在线观看视频| 黄色视频在线播放观看不卡| www.自偷自拍.com| av视频免费观看在线观看| 性高湖久久久久久久久免费观看| 99riav亚洲国产免费| 亚洲第一欧美日韩一区二区三区 | 精品少妇内射三级| 黄网站色视频无遮挡免费观看| 丁香六月天网| 精品少妇黑人巨大在线播放| 久久99一区二区三区| 精品久久久久久电影网| 午夜91福利影院| 亚洲国产成人一精品久久久| 精品视频人人做人人爽| 天天添夜夜摸| 超色免费av| 一本—道久久a久久精品蜜桃钙片| 黄片播放在线免费| 岛国在线观看网站| 亚洲美女黄片视频| 成年人免费黄色播放视频| 亚洲精品中文字幕一二三四区 | 女同久久另类99精品国产91| 波多野结衣av一区二区av| 老汉色av国产亚洲站长工具| 高清黄色对白视频在线免费看| 女警被强在线播放| 久久久久久久久免费视频了| 激情在线观看视频在线高清 | av免费在线观看网站| 欧美亚洲日本最大视频资源| 午夜福利欧美成人| 亚洲精品国产精品久久久不卡| 母亲3免费完整高清在线观看| 精品国产乱码久久久久久小说| 黑人猛操日本美女一级片| 两性午夜刺激爽爽歪歪视频在线观看 | 午夜成年电影在线免费观看| 在线观看免费高清a一片| 九色亚洲精品在线播放| 亚洲全国av大片| 99国产精品免费福利视频| 亚洲一码二码三码区别大吗| 亚洲精品乱久久久久久| 高潮久久久久久久久久久不卡| 又黄又粗又硬又大视频| 超碰97精品在线观看| 日韩视频在线欧美| 黄色视频不卡| 中亚洲国语对白在线视频| 亚洲伊人久久精品综合| 黄色怎么调成土黄色| a级毛片黄视频| 亚洲少妇的诱惑av| 他把我摸到了高潮在线观看 | 怎么达到女性高潮| 欧美激情久久久久久爽电影 | 热99久久久久精品小说推荐| 男女下面插进去视频免费观看| 可以免费在线观看a视频的电影网站| 午夜福利影视在线免费观看| 97在线人人人人妻| 欧美日韩亚洲综合一区二区三区_| 日本一区二区免费在线视频| 午夜免费鲁丝| 黄色 视频免费看| 亚洲精品美女久久av网站| 亚洲精品自拍成人| 欧美成人免费av一区二区三区 | 国产欧美日韩综合在线一区二区| 天堂8中文在线网| 亚洲中文字幕日韩| 欧美变态另类bdsm刘玥| 啪啪无遮挡十八禁网站| 老司机福利观看| 超碰成人久久| 在线观看一区二区三区激情| 性色av乱码一区二区三区2| 热99re8久久精品国产| 女人爽到高潮嗷嗷叫在线视频| 国产精品影院久久| 国产精品九九99| 欧美人与性动交α欧美精品济南到| 9色porny在线观看| 精品少妇一区二区三区视频日本电影| 久久狼人影院| 十八禁人妻一区二区| 最新在线观看一区二区三区| 搡老熟女国产l中国老女人| 亚洲成人免费av在线播放| 国产伦人伦偷精品视频| 久久久国产一区二区| 九色亚洲精品在线播放| 91麻豆精品激情在线观看国产 | 757午夜福利合集在线观看| 国产成人系列免费观看| 亚洲av日韩在线播放| 高清av免费在线| 精品一区二区三区视频在线观看免费 | 日韩成人在线观看一区二区三区| 女性生殖器流出的白浆| 99热国产这里只有精品6| 精品高清国产在线一区| 日本五十路高清| 免费观看av网站的网址| 精品久久蜜臀av无| av不卡在线播放| videos熟女内射| 丝袜美腿诱惑在线| 欧美成人免费av一区二区三区 | 成年女人毛片免费观看观看9 | 老司机午夜福利在线观看视频 | 精品人妻1区二区| 肉色欧美久久久久久久蜜桃| 少妇的丰满在线观看| 国产精品一区二区精品视频观看| 人妻一区二区av| 汤姆久久久久久久影院中文字幕| 欧美大码av| 超碰97精品在线观看| 国产精品二区激情视频| 在线观看免费日韩欧美大片| 淫妇啪啪啪对白视频| 国产高清视频在线播放一区| 国产男女超爽视频在线观看| 高清av免费在线| 日本黄色视频三级网站网址 | 国产伦理片在线播放av一区| 性少妇av在线| 精品久久久精品久久久| 亚洲美女黄片视频| 亚洲七黄色美女视频| 国产激情久久老熟女| 亚洲欧美日韩另类电影网站| 亚洲,欧美精品.| 亚洲天堂av无毛| 日本av手机在线免费观看| 欧美亚洲日本最大视频资源| 国产不卡一卡二| 午夜福利在线观看吧| 热re99久久国产66热| 无遮挡黄片免费观看| 亚洲五月婷婷丁香| 怎么达到女性高潮| 日韩制服丝袜自拍偷拍| 俄罗斯特黄特色一大片| 女人精品久久久久毛片| 国产精品1区2区在线观看. | 国产99久久九九免费精品| 黄色 视频免费看| av有码第一页| 日韩大片免费观看网站| 久久精品亚洲精品国产色婷小说| 亚洲国产毛片av蜜桃av| 午夜福利在线观看吧| 国产成人免费无遮挡视频| 99国产极品粉嫩在线观看| 亚洲人成77777在线视频| 国产麻豆69| av电影中文网址| 久久国产精品男人的天堂亚洲| 国产成人精品久久二区二区91| 国产高清videossex| 丰满饥渴人妻一区二区三| 精品少妇黑人巨大在线播放| 亚洲天堂av无毛| 在线观看免费高清a一片| 97在线人人人人妻| 性色av乱码一区二区三区2| 久久久国产一区二区| 一级a爱视频在线免费观看| 1024香蕉在线观看| tocl精华| 丝袜美足系列| 国产成人系列免费观看| av网站免费在线观看视频| av一本久久久久| 婷婷成人精品国产| 国内毛片毛片毛片毛片毛片| 国产精品久久电影中文字幕 | 婷婷丁香在线五月| 国产无遮挡羞羞视频在线观看| 国产午夜精品久久久久久| 亚洲国产av影院在线观看| 9热在线视频观看99| 亚洲,欧美精品.| 一进一出抽搐动态| 成年人午夜在线观看视频| 精品少妇黑人巨大在线播放| 久久香蕉激情| 欧美 亚洲 国产 日韩一| 人妻 亚洲 视频| av福利片在线| 午夜福利影视在线免费观看| 91麻豆精品激情在线观看国产 | 亚洲欧美日韩另类电影网站| 十八禁网站免费在线| 十分钟在线观看高清视频www| 黄网站色视频无遮挡免费观看| 亚洲男人天堂网一区| 极品人妻少妇av视频| 高清欧美精品videossex| 亚洲成人免费av在线播放| 最新的欧美精品一区二区| 中文字幕av电影在线播放| 国产成人系列免费观看| 又大又爽又粗| 国产成人精品久久二区二区免费| 嫩草影视91久久| 18禁国产床啪视频网站| 男人操女人黄网站| 精品国产国语对白av| 欧美变态另类bdsm刘玥| 纵有疾风起免费观看全集完整版| 最新美女视频免费是黄的| 国产精品久久久久久精品古装| av在线播放免费不卡| 国产av又大| netflix在线观看网站| 欧美在线黄色| 久久精品成人免费网站| 日本av手机在线免费观看| 午夜视频精品福利| 免费观看人在逋| 久久久久久久久免费视频了| 欧美av亚洲av综合av国产av| 久久午夜亚洲精品久久| 久久精品国产综合久久久| 久久久精品免费免费高清| 黄色片一级片一级黄色片| 五月天丁香电影| 1024香蕉在线观看| 国产精品亚洲av一区麻豆| 亚洲五月色婷婷综合| 肉色欧美久久久久久久蜜桃| 天堂8中文在线网| 国产精品麻豆人妻色哟哟久久| 久久久久国内视频| 黄色毛片三级朝国网站| 国产成人系列免费观看| 日韩大码丰满熟妇| 国产成+人综合+亚洲专区| 一区二区三区激情视频| 51午夜福利影视在线观看| 国产精品一区二区免费欧美| 久久久国产欧美日韩av| 18在线观看网站| 91老司机精品| www.自偷自拍.com| 电影成人av| 国产免费视频播放在线视频| 国产精品久久久av美女十八| 国产色视频综合| 一本大道久久a久久精品| 美女高潮喷水抽搐中文字幕| 18禁裸乳无遮挡动漫免费视频| 亚洲精品在线美女| 亚洲精品中文字幕一二三四区 | 好男人电影高清在线观看| 亚洲精品自拍成人| 午夜福利影视在线免费观看| www.999成人在线观看| 国产精品成人在线| 色播在线永久视频| 中文字幕人妻丝袜一区二区| 亚洲国产av影院在线观看| 另类精品久久| 精品免费久久久久久久清纯 | 在线av久久热| 日本黄色视频三级网站网址 | 日韩精品免费视频一区二区三区| 色精品久久人妻99蜜桃| 午夜91福利影院| 国产xxxxx性猛交| 久久久国产成人免费| 一边摸一边做爽爽视频免费| 免费高清在线观看日韩| 久久ye,这里只有精品| 免费看a级黄色片| 欧美成狂野欧美在线观看| 老司机亚洲免费影院| 美女高潮到喷水免费观看| 国产精品欧美亚洲77777| av天堂在线播放| 一区二区三区精品91| 一区二区三区乱码不卡18| 一本一本久久a久久精品综合妖精| 国产精品久久电影中文字幕 | 国产av国产精品国产| 成年动漫av网址| 99久久精品国产亚洲精品| 精品国产亚洲在线| 狠狠狠狠99中文字幕| 一级毛片精品| 国产精品久久久人人做人人爽| 久久久久久亚洲精品国产蜜桃av| 亚洲精品国产区一区二| 交换朋友夫妻互换小说| 精品少妇内射三级| 久久中文看片网| av又黄又爽大尺度在线免费看| 少妇粗大呻吟视频| 久久天躁狠狠躁夜夜2o2o| 午夜91福利影院| 99国产精品免费福利视频| 亚洲avbb在线观看| 欧美在线黄色| 操美女的视频在线观看| 老汉色∧v一级毛片| 久久亚洲真实| 一本综合久久免费| 岛国毛片在线播放| 正在播放国产对白刺激| 搡老岳熟女国产| 国产国语露脸激情在线看| 在线观看www视频免费| 亚洲一区二区三区欧美精品| 欧美 亚洲 国产 日韩一| 国产亚洲欧美在线一区二区| av天堂久久9| 日韩精品免费视频一区二区三区| 性少妇av在线| 成人av一区二区三区在线看| 亚洲一区二区三区欧美精品| 国产精品熟女久久久久浪| 午夜福利免费观看在线| 精品少妇久久久久久888优播| 成年动漫av网址| 欧美一级毛片孕妇| 新久久久久国产一级毛片| 中文字幕人妻熟女乱码| 国产精品免费视频内射| 亚洲一码二码三码区别大吗| 在线观看舔阴道视频| 久久久水蜜桃国产精品网| 丰满迷人的少妇在线观看| 午夜福利乱码中文字幕| 91麻豆av在线| 少妇精品久久久久久久| 精品熟女少妇八av免费久了| 天堂俺去俺来也www色官网| 精品福利永久在线观看| 又黄又粗又硬又大视频| 国产精品一区二区在线观看99| 亚洲精品久久午夜乱码| 精品乱码久久久久久99久播| www.熟女人妻精品国产| 日韩成人在线观看一区二区三区| bbb黄色大片| 亚洲国产成人一精品久久久| 国产色视频综合| 国产亚洲av高清不卡| 啪啪无遮挡十八禁网站| 免费日韩欧美在线观看| 大陆偷拍与自拍| 午夜日韩欧美国产| 五月开心婷婷网| 一级毛片精品| 国产一区二区三区综合在线观看| 久久人妻av系列| 黄片小视频在线播放| av天堂久久9| 一进一出好大好爽视频| √禁漫天堂资源中文www| 男女无遮挡免费网站观看| 久久精品国产99精品国产亚洲性色 | 18禁黄网站禁片午夜丰满| 亚洲精品国产区一区二| 国产色视频综合| 热re99久久精品国产66热6| 日本av免费视频播放| 妹子高潮喷水视频| 国产欧美日韩一区二区精品| 欧美+亚洲+日韩+国产| 精品国产一区二区久久| 大型黄色视频在线免费观看| 两人在一起打扑克的视频| 少妇精品久久久久久久| av又黄又爽大尺度在线免费看| 国产xxxxx性猛交| 欧美国产精品一级二级三级| 久久精品aⅴ一区二区三区四区| 黑人猛操日本美女一级片| kizo精华| 欧美精品av麻豆av| 欧美激情 高清一区二区三区| 黄片播放在线免费| 成人免费观看视频高清| 丰满少妇做爰视频| 国产伦人伦偷精品视频| 免费日韩欧美在线观看| 欧美精品人与动牲交sv欧美| av免费在线观看网站| 久久性视频一级片| 日韩中文字幕欧美一区二区| 成年人黄色毛片网站| 免费在线观看影片大全网站| 国产国语露脸激情在线看| 免费看a级黄色片| 一本色道久久久久久精品综合| 搡老岳熟女国产| 免费av中文字幕在线| 一进一出好大好爽视频| 精品少妇久久久久久888优播| 人人妻人人澡人人爽人人夜夜| 欧美 亚洲 国产 日韩一| 久久精品91无色码中文字幕| 国产视频一区二区在线看| 两性夫妻黄色片| 女人被躁到高潮嗷嗷叫费观| 精品国产乱码久久久久久小说| 丝瓜视频免费看黄片| 成人国语在线视频| 亚洲精品久久成人aⅴ小说| 最近最新中文字幕大全电影3 |