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

    高壓下硝基甲烷分子溫度的密度泛函計(jì)算

    2015-03-23 11:58:30董光興葛素紅李德華
    關(guān)鍵詞:聲子硝基甲烷

    董光興, 葛素紅,李德華

    (1.河西學(xué)院物理與機(jī)電工程學(xué)院, 張掖 734000; 2. 中國(guó)科學(xué)院武漢物理與數(shù)學(xué)研究所波譜與原子分子物理國(guó)家重點(diǎn)實(shí)驗(yàn)室, 武漢 430071; 3. 四川師范大學(xué)物理與電子工程學(xué)院, 成都610066 )

    高壓下硝基甲烷分子溫度的密度泛函計(jì)算

    董光興1, 葛素紅1,2,李德華3

    (1.河西學(xué)院物理與機(jī)電工程學(xué)院, 張掖 734000; 2. 中國(guó)科學(xué)院武漢物理與數(shù)學(xué)研究所波譜與原子分子物理國(guó)家重點(diǎn)實(shí)驗(yàn)室, 武漢 430071; 3. 四川師范大學(xué)物理與電子工程學(xué)院, 成都610066 )

    用密度泛函理論在B3LYP/6-311++G(2d,2P)計(jì)算水平上對(duì)硝基甲烷分子進(jìn)行了結(jié)構(gòu)優(yōu)化、 頻率和熱化學(xué)分析. 發(fā)現(xiàn): 在相同溫度條件下改變壓強(qiáng), 分子熵函數(shù)產(chǎn)生了改變, 當(dāng)溫度和壓強(qiáng)條件相同時(shí), 對(duì)于不同物質(zhì)熵函數(shù)的改變是相同的. 以熱力學(xué)理論中麥克斯韋關(guān)系為基礎(chǔ), 通過(guò)計(jì)算等溫過(guò)程中分子的熵函數(shù)對(duì)壓強(qiáng)的變化率, 用數(shù)值擬合方法得到不同壓強(qiáng)條件下分子溫度的表達(dá)式:T=T0+(1-B)[18.3858+0.5392P]V0, 式中T0、 V0分別表示分子系統(tǒng)初態(tài)的溫度和體積, T、 V分別表示系統(tǒng)在末態(tài)的溫度和體積, B是體積的壓縮比. 在選定參數(shù)的情況下該表達(dá)式可以計(jì)算不同壓強(qiáng)條件下CHNO含能材料的分子溫度. 同時(shí), 以硝基甲烷為驗(yàn)證, 選取基本參數(shù)V0和B, 計(jì)算其在C-J條件對(duì)應(yīng)的爆壓14GPa下, 分子溫度為3461K, 對(duì)應(yīng)愛(ài)因斯坦溫度, 相當(dāng)于3228cm-1的能量, 在實(shí)驗(yàn)中該能量足以激發(fā)硝基甲烷分子內(nèi)振動(dòng)能量重新分配過(guò)程, 有可能激發(fā)C-N鍵的紅外振動(dòng)而引起單分子分解反應(yīng)的發(fā)生. 因此, 此表達(dá)式可用于預(yù)測(cè)含能材料撞擊點(diǎn)火過(guò)程單分子分解可能的反應(yīng)通道.

    熱力學(xué); 麥克斯韋關(guān)系; 多聲子遷移模型; 含能材料; 硝基甲烷(NM)

    1 引 言

    大多數(shù)含能材料在常溫下是固相, 也有一部分是液態(tài). 液體在其受力作用時(shí)間很短的情況下表現(xiàn)出固體的脆性, 在一些具體的研究中可以將液體的含能材料作為瞬態(tài)的固體來(lái)處理. 含能材料撞擊點(diǎn)火過(guò)程是指材料在一個(gè)極短的時(shí)間內(nèi)受到外界一定閾值的撞擊從而發(fā)生爆炸的過(guò)程, 因此, 在對(duì)由有機(jī)分子組成的單相的含能材料撞擊點(diǎn)火的研究中, 可以將單質(zhì)的材料看成分子晶體. 對(duì)于它們撞擊點(diǎn)火過(guò)程研究的出發(fā)點(diǎn)主要有兩個(gè)方面, 其一, 從基本物理過(guò)程的角度分析, 外界撞擊在很短的時(shí)間內(nèi)給固體材料施加了強(qiáng)大的壓力, 導(dǎo)致材料內(nèi)部熱力學(xué)狀態(tài)發(fā)生變化, 溫度升高的同時(shí)內(nèi)能增加, 繼而發(fā)生爆炸. 其二, 從含能材料爆炸過(guò)程的化學(xué)性質(zhì)來(lái)看, 它是以凝聚相和氣相化學(xué)反應(yīng)為基礎(chǔ)的多組分和多階段過(guò)程, 是由于材料內(nèi)部溫度升高而導(dǎo)致的分子的分解過(guò)程及后續(xù)的系列化學(xué)反應(yīng)過(guò)程, 需要深入研究熱分解的各基元反應(yīng)歷程, 建立反應(yīng)物、 過(guò)渡態(tài)(包括中間體)和產(chǎn)物之間的聯(lián)系[1]. 含能材料的同一個(gè)撞擊點(diǎn)火過(guò)程同時(shí)包括了這兩方面的物理化學(xué)的變化過(guò)程, 材料內(nèi)部的壓強(qiáng)和溫度有急劇的變化, 點(diǎn)火過(guò)程中有作為分子解體的第一步反應(yīng)的吸熱反應(yīng)和接下來(lái)一系列的放熱反應(yīng), 這就使得研究材料溫度和壓強(qiáng)的閾值條件對(duì)撞擊點(diǎn)火過(guò)程的微觀機(jī)理有非常重要的意義.

    在含能材料領(lǐng)域內(nèi), 基于熱力學(xué)基本理論和物態(tài)方程的建立對(duì)含能金屬材料沖擊壓縮特性的實(shí)驗(yàn)和理論計(jì)算已有很多研究[2-7]. 而對(duì)由CHNO有機(jī)分子組成的含能材料在高壓下的理論和實(shí)驗(yàn)的研究比較有限. Perger W F用密度泛函方法計(jì)算了4Gpa以?xún)?nèi)太安炸藥晶體在高壓下的晶格參數(shù)和內(nèi)坐標(biāo)[8]. 逯來(lái)玉等對(duì)季戊四醇分子晶體在高溫下的結(jié)構(gòu)和振動(dòng)性質(zhì)進(jìn)行計(jì)算研究[9,10]. 他們的研究發(fā)現(xiàn)除了極少數(shù), 如C-C鍵的伸縮振動(dòng)模式是隨溫度的升高有緩慢的升高, 大多數(shù)晶體分子內(nèi)振動(dòng)的頻率隨溫度的升高都有比較平緩的降低, 他們認(rèn)為這說(shuō)明分子內(nèi)的鍵強(qiáng)隨著溫度的升高存在變?nèi)醯内厔?shì). 在分析分子晶體含能材料在撞擊點(diǎn)火過(guò)程中材料分解的微觀機(jī)理時(shí), 這樣的研究結(jié)論有非常重要的意義, 因?yàn)楣腆w材料的解體必須首先發(fā)生在材料分子內(nèi)部, 是由分子的某一個(gè)化學(xué)鍵首先離解之后而導(dǎo)致系列反應(yīng)的發(fā)生的.

    硝基甲烷(NM)是一個(gè)非常典型且分子較小的由CHNO組成的硝基復(fù)合物, 它可以用于制造炸藥, 還可用作火箭的液體燃料. 作為分子結(jié)構(gòu)較為簡(jiǎn)單的含能材料, 對(duì)它的研究能夠作為其他含能大分子體系, 尤其是由20-50個(gè)原子組成的分子系統(tǒng)的基準(zhǔn). 因此, 在正常和極端的熱力學(xué)條件下的對(duì)于硝基甲烷材料的理論和實(shí)驗(yàn)的研究是非常廣泛的. 早在1959年Hubbard和Johnson理論地計(jì)算了NM爆炸過(guò)程的熱力學(xué)參數(shù), 他們的計(jì)算結(jié)果與實(shí)驗(yàn)的結(jié)果吻合很好[11,12]. 2003年Manaa對(duì)硝基甲烷在均勻受壓和軸向受壓進(jìn)行了量子化學(xué)計(jì)算[13], 研究認(rèn)為當(dāng)體積壓縮至原來(lái)的一半時(shí), 壓強(qiáng)上升至50GPa, 最高占據(jù)軌道和最低空軌道的能級(jí)間隔減小了約0.6eV, 但是, 最高占據(jù)軌道和最低空軌道能級(jí)同時(shí)升高, 其差值卻幾乎隨體積的變化單調(diào)下降, 帶隙減小從理論上可以認(rèn)為感度在增加. 2011年Han對(duì)硝基甲烷在高溫過(guò)程進(jìn)行了模擬計(jì)算[14], 認(rèn)為硝基甲烷在300K時(shí)分解起始步驟為分子間質(zhì)子遷移, 生成CH3NOOH和CH2NO2. 在較低的溫度(2500~2000K)第一步反應(yīng)則為異構(gòu)化生成CH3ONO的過(guò)程. 同時(shí)在反應(yīng)剛開(kāi)始極短時(shí)間內(nèi)可伴隨分子內(nèi)質(zhì)子遷移生成CH3NOOH. 而對(duì)于硝基甲烷分子的分解機(jī)理的研究則有很多報(bào)道[15-21]. 對(duì)它的光解研究表明, 氣態(tài)情況下分子沿C-N鍵的均裂主要生成甲基和二氧化氮[18]. 電子碰撞光譜研究報(bào)道了硝基甲烷分子的最低單態(tài)到三態(tài)的激發(fā)躍遷[19]. 有作者在實(shí)驗(yàn)中觀測(cè)分析認(rèn)為在硝基甲烷光解為甲基和二氧化氮的過(guò)程中, 分子態(tài)躍遷起著很大的作用, 他們認(rèn)為合理的反應(yīng)機(jī)理應(yīng)該是分子先激發(fā)為三態(tài)然后再分解為甲基和二氧化氮的[20].

    但是, 從整個(gè)撞擊點(diǎn)火反應(yīng)過(guò)程來(lái)看, 首先是外界的機(jī)械撞擊引起了材料的壓強(qiáng)和溫度等物態(tài)參數(shù)的改變, 其次是材料分子的某個(gè)化學(xué)鍵受到激發(fā)而發(fā)生斷裂, 引起了整個(gè)的系列反應(yīng). 撞擊機(jī)械能引起的材料內(nèi)部溫度和壓強(qiáng)條件的變化是材料分子分解反應(yīng)發(fā)生的先決條件, 這就需要在研究由CHNO有機(jī)分子組成的含能材料撞擊點(diǎn)火機(jī)理時(shí)考慮外界撞擊對(duì)材料內(nèi)部產(chǎn)生的影響, 將材料的熱力學(xué)狀態(tài)變化和材料單分子分解條件相互聯(lián)系. 因此, 本文選取具有代表性的含能材料---硝基甲烷為研究對(duì)象, 以熱力學(xué)理論為基礎(chǔ), 在密度泛函算法基礎(chǔ)上推導(dǎo)了含能材料分子高壓下的溫度壓強(qiáng)關(guān)系式, 并借鑒多聲子遷移模型的基本理論預(yù)測(cè)了硝基甲烷在撞擊誘導(dǎo)點(diǎn)火過(guò)程的基元反應(yīng).

    2 理論模型和計(jì)算方法

    2.1 多聲子遷移模型的基本思想

    在撞擊點(diǎn)火過(guò)程中外界沖擊和撞擊產(chǎn)生的機(jī)械能是被轉(zhuǎn)移到材料的分子反應(yīng)中心的, 這說(shuō)明撞擊點(diǎn)火過(guò)程是由兩個(gè)不同的過(guò)程組成的. 多聲子遷移模型認(rèn)為這說(shuō)明外界撞擊產(chǎn)生的機(jī)械能是經(jīng)過(guò)在固體中激發(fā)聲子, 再經(jīng)聲子將能量傳遞轉(zhuǎn)移至分子鍵, 激發(fā)分子鍵的紅外振動(dòng), 導(dǎo)致分子鍵的斷裂, 激發(fā)隨后而來(lái)的放熱化學(xué)反應(yīng), 導(dǎo)致區(qū)域溫度急速上升, 熱點(diǎn)形成, 材料迅速點(diǎn)火爆炸[22-25]. 晶體材料的聲子振動(dòng)模式頻率較低, 而對(duì)應(yīng)分子紅外振動(dòng)頻率卻分布較廣. 多聲子遷移模型認(rèn)為, 被機(jī)械撞擊激發(fā)的聲子經(jīng)雙聲子吸收過(guò)程, 先激發(fā)分子中較低的紅外振動(dòng)模式, 將這種低頻紅外振動(dòng)模式稱(chēng)為門(mén)模式, 門(mén)模式的頻率一般在(200,700)cm-1. 然后, 門(mén)振動(dòng)模式作為低頻振動(dòng)模式, 可與聲子耦合. 當(dāng)能量持續(xù)從過(guò)熱聲子區(qū)流入到門(mén)振動(dòng), 經(jīng)過(guò)多聲子遷移, 高一級(jí)的分子振動(dòng)模式會(huì)通過(guò)分子內(nèi)振動(dòng)能量的遷移而被門(mén)模式所激發(fā), 最終達(dá)到分子內(nèi)振動(dòng)能量的某一閾值, 然后, 形成了分子紅外振動(dòng)能量在分子內(nèi)的重新分配過(guò)程, 這種振動(dòng)態(tài)的遷移一直持續(xù)到有足夠的能量被定域到導(dǎo)致反應(yīng)的分子反應(yīng)路徑上, 導(dǎo)致鍵的裂解, 分子解體, 第一步吸熱反應(yīng)發(fā)生, 接下來(lái)會(huì)在局部發(fā)生一系列的放熱反應(yīng), 導(dǎo)致溫度升高, 熱點(diǎn)形成, 點(diǎn)火完成[23]. 硝基甲烷分子振動(dòng)能量重新分配過(guò)程的實(shí)驗(yàn)和理論模擬都有研究報(bào)道[26,27]. 實(shí)驗(yàn)中首先激發(fā)了硝基甲烷分子C-H鍵的紅外振動(dòng)(~3000cm-1), 發(fā)現(xiàn)分子的紅外振動(dòng)遷移主要由三個(gè)過(guò)程組成, 先是被激發(fā)的高能量紅外振動(dòng)在幾個(gè)皮秒內(nèi)振動(dòng)冷卻, 激發(fā)了化學(xué)鍵的能量較低的紅外振動(dòng), 范圍大致為(~1600-1400cm-1). 然后能量較低的這些振動(dòng)發(fā)生弛豫過(guò)程, 激發(fā)了化學(xué)鍵的能量更低的紅外振動(dòng)模式, 范圍大致為(~1100-480cm-1). 最后, 被激發(fā)的這些紅外振動(dòng)模式發(fā)生振動(dòng)冷卻過(guò)程, 從而激發(fā)聲子振動(dòng)模式[26].

    因此, 在撞擊點(diǎn)火過(guò)程中, 當(dāng)存在足夠的能量以激發(fā)材料分子能量較高的紅外振動(dòng)模式, 且達(dá)到分子內(nèi)振動(dòng)能量重新分配的閾值, 則會(huì)有這樣階梯式的分子紅外振動(dòng)能量遷移. 當(dāng)作為分子分解的化學(xué)鍵的紅外振動(dòng)被激發(fā), 就會(huì)造成該化學(xué)鍵的斷裂, 單分子分解反應(yīng)發(fā)生, 引發(fā)接下來(lái)的點(diǎn)火爆炸過(guò)程.

    2.2 高壓下分子溫度的預(yù)測(cè)和計(jì)算原理

    到目前為止, 雖然多聲子遷移模型是一個(gè)比較完美的理論, 用它來(lái)解釋含能材料的撞擊點(diǎn)火過(guò)程是合乎邏輯的. 但是, 由于聲子振動(dòng)模在實(shí)驗(yàn)上難以檢測(cè)和驗(yàn)證, 多聲子遷移模型理論被認(rèn)為缺乏直接的實(shí)驗(yàn)依據(jù)[28]. 無(wú)論如何, 含能材料的撞擊點(diǎn)火過(guò)程是將撞擊機(jī)械能傳遞到分子內(nèi)部而導(dǎo)致分子的解體的, 多聲子遷移模型認(rèn)為機(jī)械能轉(zhuǎn)化為聲子振動(dòng)的能量后遷移至分子內(nèi)部的振動(dòng)能, 從能量守恒的角度分析, 也可以認(rèn)為是撞擊機(jī)械能轉(zhuǎn)化為熱能而后遷移至分子內(nèi)部激發(fā)分子內(nèi)部化學(xué)鍵的紅外振動(dòng).

    對(duì)于多原子體系由N個(gè)單原子組成的分子本身, 或者由M個(gè)分子組成的晶體原胞都可以看成是宏觀熱力學(xué)系統(tǒng). 則系統(tǒng)的熱力學(xué)狀態(tài)函數(shù)的改變只取決于狀態(tài)而與過(guò)程無(wú)關(guān), 依據(jù)熱力學(xué)基本理論麥克斯韋關(guān)系式[29]:

    (1)

    等壓過(guò)程中體積隨溫度的變化率等于等溫過(guò)程狀態(tài)函數(shù)熵對(duì)壓強(qiáng)的偏導(dǎo)數(shù)的負(fù)值. 若將撞擊點(diǎn)火過(guò)程中系統(tǒng)溫度和壓強(qiáng)的急劇變化過(guò)程看成是一系列等壓過(guò)程或者等溫過(guò)程的積累, 在撞擊點(diǎn)火過(guò)程中分子系統(tǒng)內(nèi)部的溫度壓強(qiáng)關(guān)系可以通過(guò)計(jì)算熵函數(shù)的改變而獲得:

    (2)

    在一無(wú)限小等壓過(guò)程中, 可將偏導(dǎo)數(shù)表示為物理近似關(guān)系式:

    (3)

    由(3)式可以獲得系統(tǒng)的溫度壓強(qiáng)關(guān)系式:

    (4)

    式中T0、V0分別表示分子系統(tǒng)初態(tài)的溫度和體積,T、V分別表示系統(tǒng)在末態(tài)的溫度和體積, B是體積的壓縮比, 在正常情況下, 要求其小于1, 定義為B=V/V0.

    (5)

    由于第一性原理的密度泛函理論(DFT)方法系統(tǒng)誤差要比HF方法小[31], 且密度泛函B3LYP雜化方法, 不但對(duì)交換能做密度梯度校正, 而且對(duì)相關(guān)能也做了梯度校正, 提高了密度泛函計(jì)算精確度[32]. 本文采用密度泛函理論方法, 在B3LYP/6-311++(2d,2P)計(jì)算水平上, 對(duì)硝基甲烷的結(jié)構(gòu)進(jìn)行了全優(yōu)化, 并在同一水平上計(jì)算了298.15K溫度、 不同壓強(qiáng)下的熱化學(xué)參數(shù). 所有計(jì)算工作用Gaussian09程序包完成. 頻率計(jì)算驗(yàn)證無(wú)虛頻, 而且有六個(gè)幾乎為零的頻率, 說(shuō)明計(jì)算研究是在最小勢(shì)能面上進(jìn)行的.

    表1 硝基甲烷分子的幾何參數(shù)

    鍵長(zhǎng)單位(?);鍵角和二面角單位(o)

    表2 不同壓強(qiáng)下硝基甲烷分子的熵和自由能的計(jì)算值

    3 計(jì)算結(jié)果與討論

    3.1 硝基甲烷單分子溫度的模擬公式

    圖1 等溫過(guò)程分子熵隨壓強(qiáng)變化率與壓強(qiáng)的關(guān)系(1-100)atmFig.1 The molecular entropy change rate with pressure versus pressure of NM in (1-100)atm

    圖2 等溫過(guò)程分子熵隨壓強(qiáng)變化率與壓強(qiáng)的關(guān)系(100-1000)atmFig.2 The molecular entropy change rate with pressure versus pressure of NM in (100-1000)atm

    首先運(yùn)用密度泛函理論中的B3LYP/6-311++G(2d,2P)方法對(duì)基態(tài)的硝基甲烷分子進(jìn)行幾何優(yōu)化, 得到穩(wěn)定的幾何構(gòu)型, 確保計(jì)算是在勢(shì)能面的最低點(diǎn)進(jìn)行. 分子的幾何構(gòu)型參數(shù)列在表1中. 根據(jù)文獻(xiàn), 基態(tài)的硝基甲烷分子穩(wěn)定構(gòu)型參數(shù)是:C-H鍵鍵長(zhǎng)為1.088 ?, C-N鍵鍵長(zhǎng)為1.489 ?, N-O鍵鍵長(zhǎng)為1.224 ?, 鍵角∠NCH為107°, ∠ONO為125.3°[33]. 表1中基態(tài)分子的幾何構(gòu)型參數(shù)與其實(shí)驗(yàn)值相比較, 鍵長(zhǎng)的最大偏差為0.005?, 鍵角的誤差基本是1°, 這樣小的誤差說(shuō)明分子的構(gòu)型優(yōu)化是比較成功的, 在這個(gè)基礎(chǔ)上再進(jìn)行頻率分析計(jì)算就會(huì)得到更可靠的結(jié)果. 然后, 保持溫度在298.15K的條件, 改變壓強(qiáng)參數(shù), 在B3LYP/6-311++G(2d,2P)計(jì)算這個(gè)構(gòu)型下的硝基甲烷分子的熱化學(xué)參數(shù), 計(jì)算發(fā)現(xiàn)分子的熵函數(shù)和自由能(經(jīng)過(guò)零點(diǎn)修正)的數(shù)值隨壓強(qiáng)的改變而改變. 計(jì)算數(shù)據(jù)在表2中列出.

    圖3 等溫過(guò)程分子熵隨壓強(qiáng)變化率與壓強(qiáng)的關(guān)系(1000-8000)atmFig.3 The molecular entropy change rate with pressure versus pressure of NM in (1000-8000)atm

    圖4 等壓過(guò)程分子溫度隨體積的變化率與壓強(qiáng)的關(guān)系圖線Fig.4 The molecular temperature change rate with volume versus pressure of NM

    根據(jù)表2數(shù)據(jù)計(jì)算等溫過(guò)程分子體系的熵函數(shù)隨壓強(qiáng)的變化率并擬合變化率與壓強(qiáng)的關(guān)系圖形, 按壓強(qiáng)分區(qū)繪制如圖1、 圖2和圖3所示. 根據(jù)表2數(shù)據(jù), 依據(jù)公式(2)計(jì)算不同壓強(qiáng)條件下溫度隨體積的變化率, 擬合得到溫度對(duì)體積的變化率與壓強(qiáng)呈線性, 如圖4所示, 其相關(guān)系數(shù)R為0.99927,R2等于0.9985, 這個(gè)數(shù)據(jù)說(shuō)明數(shù)據(jù)的線性相關(guān)性較好.擬合所得方程為:

    (6)

    擬合方程(6)給出了硝基甲烷材料分子在等壓過(guò)程中溫度對(duì)體積的變化率.將(6)式代入(4)式可得最終的擬合推導(dǎo)結(jié)果, 即不同壓強(qiáng)下分子溫度的計(jì)算公式:

    =T0+(1-B)[18.3858+0.5392P]V0

    (7)

    圖5 硝基甲烷分子結(jié)構(gòu)圖Fig.5 Optimized geometrical configurations of NM

    3.2 不同條件下硝基甲烷分子分解通道分析

    常溫下硝基甲烷是液體狀態(tài), 人們沒(méi)有獲得它的撞擊感度實(shí)驗(yàn)值. 一般情況下人們將含能材料爆轟時(shí)所產(chǎn)生的爆轟波陣面的壓力稱(chēng)為爆壓, 也稱(chēng)為C-J壓力.實(shí)驗(yàn)上可以得到常用炸藥材料的爆壓為10-40GPa, 硝基甲烷的爆轟壓力可達(dá)到14GPa[34]. 傳統(tǒng)的固體含能材料熱力學(xué)爆炸理論是有名的Chapman-Jouguet模型[35], 這個(gè)模型認(rèn)為C-J溫度和C-J壓力就是含能材料的爆轟條件. 因此選取14GPa壓強(qiáng)條件作為硝基甲烷的點(diǎn)火爆轟條件.

    圖5所示為硝基甲烷分子結(jié)構(gòu)圖, 考慮電子云分布, 微觀上將單個(gè)分子看成球型結(jié)構(gòu), C-N鍵中點(diǎn)為球心, C-N鍵長(zhǎng)的一半加上N—O鍵和C-H鍵中最長(zhǎng)的數(shù)值可以認(rèn)為是球的半徑, 由此估算得硝基甲烷分子的體積為6π×10-30m3. 由于撞擊點(diǎn)火過(guò)程中含能材料受到猛烈的外界撞擊, 外界撞擊的力遠(yuǎn)大于分子間力的作用, 因此, 對(duì)于結(jié)構(gòu)疏松的分子晶體, 忽略分子間力的作用,將未受撞擊作用時(shí)硝基甲烷分子的體積作為初態(tài)體積, 則有V0=6π×10-30m3, 在撞擊發(fā)生之前材料與環(huán)境達(dá)到熱平衡, 選取分子的初始溫度T0等于環(huán)境溫度298.15K. 同樣,在這種近似條件下, 可將疏松結(jié)構(gòu)的分子晶體在撞擊前后的體積的壓縮比B值選取為0.793[23]. 依據(jù)(7)式可計(jì)算得在壓強(qiáng)達(dá)到14GPa時(shí), 硝基甲烷分子的溫度值為:3461K, 這樣的分子溫度相當(dāng)于0.4eV的能量, 對(duì)應(yīng)愛(ài)因斯坦溫度, 相當(dāng)于3228cm-1的能量.

    表3 硝基甲烷分子的振動(dòng)頻率

    a本文理論計(jì)算的頻率值;b實(shí)驗(yàn)文獻(xiàn)[26]中的頻率值

    硝基甲烷分子C-H鍵的紅外振動(dòng)能量就在3000cm-1左右, 如表3所示是硝基甲烷分子的紅外振動(dòng)光譜理論計(jì)算值和實(shí)驗(yàn)值. 按照Deak等人的實(shí)驗(yàn)結(jié)果, 在激發(fā)了硝基甲烷分子C-H鍵的紅外振動(dòng)(~3000cm-1), 之后, 發(fā)現(xiàn)激發(fā)的高能量紅外振動(dòng)在幾個(gè)皮秒內(nèi)振動(dòng)冷卻,激發(fā)了化學(xué)鍵的能量較低的紅外振動(dòng), 范圍大致為(~1600-1400cm-1), 然后能量較低的這些振動(dòng)發(fā)生弛豫過(guò)程, 激發(fā)了化學(xué)鍵的能量更低的紅外振動(dòng)模式, 范圍大致為(~1100-480cm-1). 最后, 被激發(fā)的這些紅外振動(dòng)模式發(fā)生振動(dòng)冷卻過(guò)程, 從而激發(fā)聲子振動(dòng)模式[26]. 那么, 當(dāng)外界撞擊導(dǎo)致材料分子所受壓強(qiáng)達(dá)到14GPa時(shí), 分子溫度所對(duì)應(yīng)的能量相當(dāng)于3228cm-1, 足以激發(fā)硝基甲烷分子C-H鍵的紅外振動(dòng), 這將導(dǎo)致分子振動(dòng)能量重新分配過(guò)程的發(fā)生, 實(shí)驗(yàn)顯示, 這個(gè)過(guò)程會(huì)激發(fā)化學(xué)鍵的能量較低的紅外振動(dòng), 其范圍分別為(~1600-1400cm-1)和(~1100-480cm-1). 當(dāng)CN鍵的伸縮振動(dòng)模式918(1379)cm-1被激發(fā)后, 硝基甲烷分子就有可能沿著這個(gè)反應(yīng)坐標(biāo)分解產(chǎn)生分解反應(yīng), 導(dǎo)致材料分子解體, 引發(fā)接下來(lái)的一系列鏈?zhǔn)椒磻?yīng), 使點(diǎn)火反應(yīng)發(fā)生:

    CH3NO2→CH3+NO2

    (8)

    因此, 可以認(rèn)為本文的理論計(jì)算公式與硝基甲烷分子振動(dòng)能量重新分配過(guò)程的實(shí)驗(yàn)和理論模擬的研究結(jié)果相一致[26,27]. 可以認(rèn)為14GPa就是硝基甲烷材料C-J條件. 在本文的近似條件下, 選取相同的參數(shù)V0和B, 根據(jù)公式(7)還可以計(jì)算不同壓強(qiáng)條件下的分子溫度和其對(duì)應(yīng)的振動(dòng)能量, 計(jì)算結(jié)果在表4中列出. 從表4數(shù)據(jù)可以看出當(dāng)壓強(qiáng)條件達(dá)到12GPa時(shí), 分子溫度對(duì)應(yīng)的能量數(shù)值已經(jīng)達(dá)到了實(shí)驗(yàn)中需要激發(fā)分子C-H鍵紅外振動(dòng)的能量, 但材料實(shí)驗(yàn)中硝基甲烷的爆壓卻是14GPa. 造成這一偏差的原因在于熱能的傳遞和利用有效率的問(wèn)題, 在被有效利用的過(guò)程中熱能總是存在著散失. 因此, 在表4的計(jì)算中選取數(shù)據(jù)時(shí), 計(jì)算值都要大于實(shí)驗(yàn)給出的數(shù)值, 由此來(lái)推斷在對(duì)應(yīng)壓強(qiáng)下可能被激發(fā)的分子的紅外振動(dòng)模式, 這正是推斷晶體材料單分子分解情況的前提.

    表4 不同壓強(qiáng)下硝基甲烷分子溫度和能量

    就目前的實(shí)驗(yàn)事實(shí)發(fā)現(xiàn), 在激發(fā)了硝基甲烷分子C-H鍵的紅外振動(dòng)之后, 分子中發(fā)生了振動(dòng)能量的重新分配過(guò)程, 使得CN鍵的伸縮振動(dòng)模式或者其他的紅外振動(dòng)模式有可能被激發(fā), 而導(dǎo)致單分子分解反應(yīng)的發(fā)生. 但是, 實(shí)驗(yàn)并沒(méi)有給出這樣一種情形, 那就是直接激發(fā)分子的C-H鍵的紅外振動(dòng)模式或者其他的紅外振動(dòng)模式之后, 會(huì)不會(huì)存在晶體材料單分子分解反應(yīng)的發(fā)生. 因此, 在表4中, 只有當(dāng)壓強(qiáng)達(dá)到或者大于14GPa, 才可以確定地認(rèn)為從本文的理論預(yù)測(cè)和實(shí)驗(yàn)兩方面都能判斷分子的CN鍵的伸縮振動(dòng)模式可能被激發(fā), 這個(gè)條件下會(huì)有分子分解反應(yīng)的發(fā)生, 反應(yīng)路徑沿著C-N鍵. 當(dāng)壓強(qiáng)達(dá)到3GPa時(shí), 分子溫度所提供的能量就足夠激發(fā)CN鍵的伸縮振動(dòng)模式了, 在其他小于14GPa的壓強(qiáng)條件下, 分子溫度所提供的能量也足夠激發(fā)另外的紅外振動(dòng)模式, 分子解體的可能情況在表4中可以對(duì)應(yīng)分析, 但是反應(yīng)發(fā)生的可能性無(wú)從判斷. 對(duì)于這樣的情形, 需要實(shí)驗(yàn)研究給出證明.

    4 結(jié) 論

    撞擊點(diǎn)火過(guò)程是指含能材料受到猛烈的外界撞擊, 從而導(dǎo)致材料的點(diǎn)火反應(yīng)發(fā)生引起爆炸的過(guò)程. 由于外界撞擊的力遠(yuǎn)大于分子間的相互作用, 因此, 對(duì)于結(jié)構(gòu)疏松的分子晶體, 忽略分子間力的作用. 本文基于這樣的近似處理, 以熱力學(xué)基本理論為基礎(chǔ), 將材料的撞擊點(diǎn)火過(guò)程看成系列等壓過(guò)程, 通過(guò)計(jì)算多粒子的分子體系熵函數(shù)的改變推導(dǎo)擬合得到一個(gè)隨壓強(qiáng)而變化的分子溫度函數(shù)關(guān)系式, 進(jìn)而計(jì)算不同壓強(qiáng)條件下材料分子的溫度. 在愛(ài)因斯坦特征溫度的概念下, 將分子溫度換算成能量數(shù)值, 配合現(xiàn)存的實(shí)驗(yàn)結(jié)論, 由此判斷該溫度閾值是否足夠激發(fā)分子紅外振動(dòng)正則模式, 并在分子內(nèi)部形成振動(dòng)能量的重新分配過(guò)程. 本文選取硝基甲烷為研究分子, 計(jì)算并擬合出高壓下分子溫度公式, 選擇材料的初狀態(tài)和分子的基本參數(shù), 計(jì)算得到硝基甲烷爆壓條件14GPa下的分子溫度為3461K, 這樣的分子溫度相當(dāng)于0.4eV的能量, 對(duì)應(yīng)愛(ài)因斯坦溫度, 相當(dāng)于3228cm-1的能量.對(duì)比DeaK等人的實(shí)驗(yàn)研究, 該能量正好可以激發(fā)硝基甲烷分子C-H鍵的紅外振動(dòng)模式而導(dǎo)致分子內(nèi)振動(dòng)能量重新分配過(guò)程, 這樣的分配過(guò)程最終會(huì)激發(fā)C-N鍵的紅外振動(dòng)而引起單分子分解反應(yīng)的發(fā)生. 這說(shuō)明本文擬合的理論公式計(jì)算的結(jié)果與硝基甲烷材料的爆轟條件相符合, 可用于預(yù)測(cè)計(jì)算CHNO含能材料在高壓下的分子溫度, 或者計(jì)算預(yù)測(cè)含能材料產(chǎn)生點(diǎn)火反應(yīng)的壓強(qiáng)條件, 在有關(guān)于分子振動(dòng)能量重新分配過(guò)程的實(shí)驗(yàn)結(jié)果支持的條件下, 本文的計(jì)算結(jié)果也可用于預(yù)測(cè)外界撞擊條件下單分子分解反應(yīng)的可能路徑和通道.

    [1]JuXH,YeCC,XuSY.Overviewonquantumchemicalcomputingandmoleculardynamicsimulationsofenergeticmaterials[J].ChineseJournalofExplosives&Propellant, 2012, 35(2): 1 (in Chinese)[ 居學(xué)海, 葉財(cái)超, 徐司雨.含能材料的量子化學(xué)計(jì)算與分子動(dòng)力學(xué)模擬綜述[J]. 火炸藥學(xué)報(bào), 2012, 35(2): 1]

    [2] Shi A S, Zhang X F, Qiao L,etal. Theoretical calculation on shock compression characteristics of multifunctional energetic structural material[J].ExplosionandShockWaves, 2013, 33(2): 148(in Chinese) [史安順, 張先鋒, 喬良, 等. 多功能含能結(jié)構(gòu)材料沖擊壓縮特性的理論計(jì)算[J]. 爆炸與沖擊, 2013, 33(2): 148]

    [3] Zhang X F, Qiao L, Shi A S. A cold energy mixture theory for the equation of state in solid and porous metal mixtures[J].JournalofAppliedPhysics, 2011, 110(1): 013506.

    [4] Wu Q, Jing F Q. Unified thermodynamic equation-of-state for porous materials in a wide pressure range [J].JournalofAppliedPhysics, 1995, 57(1): 49.

    [5] Jordan J L, Ferranti L, Austin R. Equation of state of aluminum-iron oxide-epoxy composite [J].JournalofAppliedPhysics, 2007, 101(9): 093520.

    [6] Yuan X L, Wei D Q, Cheng Y,etal. Thermodynamics properties of Zr4Al2under high pressure from first-principle calculations[J].J.At.Mol.Phys., 2012, 29(2): 325 [苑曉麗, 魏冬青, 程艷, 等. Zr4Al2高壓下熱動(dòng)力學(xué)性質(zhì)的第一性原理研究[J]. 原子與分子物理學(xué)報(bào), 2012, 29(2): 325]

    [7] Wang J, Guo F, Cheng X L,etal. Reactive molecular dynamics simulations of initial decomposition of TATB under high temperature and high pressure[J].JournalofSichuanUniversity:NaturalScienceEdition, 2013, 50(3): 580(in Chinese) [王君, 郭峰, 程新路, 等. TATB高溫高壓下初始分解反應(yīng)的分子動(dòng)力學(xué)模擬[J].四川大學(xué)學(xué)報(bào): 自然科學(xué)版, 2013, 50(3): 580]

    [8] Perger W F, Vutukuri S, Dreger Z A,etal. First principles vibrational studies of pentaerythritol crystal under hydrostatic pressure[J].Chem.Phys.Lett., 2006, 422: 397.

    [9] Lu L Y, Zhou X L, Ji G F,etal. Structural and vibrational properties underhigh temperature of solid pentaerythritol by molecular dynamics simuations[J].J.At.Mol.Phys., 2012, 29(5): 908 (in Chinese) [逯來(lái)玉, 周小林, 姬廣富, 等. 高溫季戊四醇晶體的結(jié)構(gòu)和振動(dòng)性質(zhì)的分子動(dòng)力學(xué)模擬[J]. 原子與分子物理學(xué)報(bào), 2012, 29(5): 908]

    [10] Lu L Y, Wei D Q, Chen X R,etal. First principle calculations of structures and electronic properties of solid pentaerythritol under pressure[J].Chin.Phys.Lett., 2008, 25: 3368.

    [11] Hubbard H W, Johnson M H. Initiation of detomations [J].J.Appl.Phys., 1959, 30: 765.

    [12] Campbell A W, Davis W C, Travis J R. Shock initiation of detonation in liquid explosives [J].Phys.Fluids., 1961, 4: 498.

    [13] Manaa M R, Fried L E, Reed E J. Explosive chemistry: simulating the chemistry of energetic materials at extreme conditions [J].J.Comput.Aid.Mol.Des, 2003, 10(2): 75.

    [14] Han S P, Van Duin A C T, Goddard Ⅲ W A,etal. Thermal decomposition of condensed-phase nitromethane from molecular dynamics from ReaxFF reactive dynamics [J].J.Phys.Chem. B, 2011, 115: 6534.

    [15] Blais N C, Engelke R, Sheffield S A. Mass spectroscopic study of the chemical reaction zone in detonating liquid nitromethane [J].J.Phys.Chem. A, 1997, 101: 8285.

    [16] Gruzdkov Y A, Gupta Y M. Emission and fluorescence spectroscopy to examine shock-induced decomposition in nitromethane [J].J.Phys.Chem. A, 1998, 102: 8325.

    [17] Miller P J, Piermarini G J, Block S. An FT-IR micro-spectroscopic method for kinetic measurements at high temperatures and high pressures [J].AppliedSpectroscopy, 1984, 38: 680.

    [18] Moss D B, Trentelmen K A, Houston P L. 193 nm photodissociation dynamics of nitromethane [J].J.Chem.Phys., 1992, 96: 237.

    [19] Flicker W M, Mosher O A, Kuppermann A. Variable angle electron-impact excitation of nitromethane [J].J.Chem.Phys., 1980, 72: 2788.

    [20] Honda K, Mikuni H, Takahasi M. Photolysis of nitromethane in gas phase at 313 nm[J].Bull.Chem.Soc.Jpn., 1972, 45: 3534.

    [21] Liu H, Zhao J J, Wei D Q,etal. Structural and vibrational properties of solid nitromethane under high pressure by density functional theory [J].J.Chem.Phys., 2006, 124: 124501.

    [22] Dlott D D. Multi-phonon up-pumping in energetic materials [J].Overviewsofrecentresearchonenergeticmaterials.Adv.Ser.Phys.Chem., 2005, 16: 303.

    [23] Tokmakoff A, Fayer M D. Chemical reaction initiation and hot-spot formation in shocked energetic molecular materials [J].J.Phys.Chem., 1993, 97: 1901.

    [24] Dlott D D, Fayer M D. Shocked molecular solids: vibrational up pumping, defect hot spot formation, and the onset of chemistry [J].J.Phys.Chem., 1990, 92: 3797.

    [25] Dlott D D. Picosecond dynamics behind the shock front [J].LeJournaldephysiqueⅣ, 1995, 5(4): C4-337.

    [26] Deak J C, Iwaki L K, Dlott D D. Vibrational energy redistribution in ployatomic liquids: ultrafast IR-Raman spectroscopy of nitromethane [J].J.Phys.Chem. A, 1999, 25: 971.

    [27] Kabadi V N, Rice B M. Molecular dynamics simulations of normal mode vibrational energy transfer in liquid nitromethane [J]J.Phys.Chem. A, 2004, 108: 532.

    [28] Zeman S. Sensitivities of high energy compounds [J].Struct.Bond., 2007, 12: 195.

    [29] Wang Z C.Thethermodynamicstatisticalphysics[M]. Beijing: Higher Education Press, 2003(in Chinese)[汪志誠(chéng). 熱力學(xué)統(tǒng)計(jì)物理學(xué)[M]. 北京: 高等教育出版社, 2003]

    [30] Huang K.Solidstatephysics[M]. Beijing: Higher Education Press, 1988(in Chinese)[黃昆. 固體物理[M]. 北京: 高等教育出版社, 1988]

    [31] Wiberg K B, Hadad C M, Lepage T J,etal. An analysis of the effect of electron correlation on charge density distribution [J].J.Chem.Phys., 1992, 96: 671.

    [32] Joseph W, Ochterski Ph D.ThermochemistryinGaussian[M]. Gausssian Inc. 13, 2000.

    [33] Huber K P, Herzberg G.Molecularspectraandmolecularstructure:constantsofdiatomicmolecules[M]. London: Van Nostrand Reinhold Co, 1979.

    [34] Hardesty D R. An investigation of the shock initiation of liquid nitromethane [J].CombustionandFlame, 1976, 27: 229.

    [35] Zeldovich Y B, RAizer Y P.Physicsofshockwavesandhightemperaturehydrodynamicphenomena[M]. New York: Academic Press, 1966.

    Density functional theory study on the temperature of nitromethane molecule under pressure

    DONG Guang-Xing1,GE Su-Hong1,2, LI De-Hua3

    (1. Department of Physics, Hexi University, Zhangye 734000, China; 2. State Key Laboratory of Magnetic Resonance and Atomic and Molecular Physics, Wuhan Institute of Physics and Mathematics, Chinese Academy of Sciences, Wuhan 430071, China; 3. College of Physics and Electronic Engineering, Sichuan Normal University, Chengdu 610066, China)

    In the present paper B3LYP method complied with the 6-311++G(2p,2d) basis set was employed to investigate the molecular property of nitromethane. Based on a systematic study of the molecular structure, electronic energy, normal vibration mode and thermal chemical functions in nitromethane via the lowest singlet, it is found that the decrease of entropy with pressure is evident if we keep the temperature constant. Based on Maxwell relations of the thermodynamic we deduced the fitting equation as follow:T=T0+(1-B)[18.3858+0.5392P]V0, which can calculate the temperature of the energetic molecules under different pressures. Further more, the temperature of the nitromethane under 14GPa is 3461K, which equates to 3228cm-1. This energy value could excite the C-H stretching and the intramolecular vibrational redistribution process, the process could generate the decomposition of the molecule. Thus, the fitting equation can expect the decomposed path way of the energetic molecules.

    Thermodynamic; Maxwell relations; Multiponon transfer mode; Energetic Materials; Nitromethane

    103969/j.issn.1000-0364.2015.02.002

    2014-06-25

    國(guó)家自然科學(xué)基金(41171175)

    董光興(1971—),男,甘肅平?jīng)鋈?,副教授,碩士,主要從事光學(xué)、量子力學(xué)等的教學(xué)科研工作.

    葛素紅.E-mail: gesuhong@163.com

    O642

    A

    1000-0364(2015)02-0181-09

    猜你喜歡
    聲子硝基甲烷
    半無(wú)限板類(lèi)聲子晶體帶隙仿真的PWE/NS-FEM方法
    納米表面聲子 首次實(shí)現(xiàn)三維成像
    硝基胍烘干設(shè)備及工藝研究
    化工管理(2021年7期)2021-05-13 00:46:24
    聲子晶體覆蓋層吸聲機(jī)理研究
    液氧甲烷發(fā)動(dòng)機(jī)
    論煤炭運(yùn)輸之甲烷爆炸
    水上消防(2020年1期)2020-07-24 09:26:02
    Gas from human waste
    基于聲子晶體理論的導(dǎo)線防舞方法及數(shù)值驗(yàn)證
    高塔硝基肥,科技下鄉(xiāng)助農(nóng)豐收
    九硝基三聯(lián)苯炸藥的合成及表征
    国产精品98久久久久久宅男小说| 久久亚洲精品不卡| 小蜜桃在线观看免费完整版高清| 小说图片视频综合网站| 在线观看美女被高潮喷水网站 | 很黄的视频免费| 丁香欧美五月| 久久久久久久亚洲中文字幕 | 天堂动漫精品| 国产三级中文精品| 欧美中文综合在线视频| 午夜福利免费观看在线| 国产亚洲精品久久久久久毛片| 人妻久久中文字幕网| 国产精品,欧美在线| 欧美日韩瑟瑟在线播放| 狂野欧美激情性xxxx| 亚洲狠狠婷婷综合久久图片| 午夜福利高清视频| 男女视频在线观看网站免费| 亚洲av一区综合| 18禁黄网站禁片免费观看直播| 热99re8久久精品国产| 国内久久婷婷六月综合欲色啪| 国产一区二区亚洲精品在线观看| 色在线成人网| 亚洲在线观看片| 麻豆国产97在线/欧美| 国产高清视频在线播放一区| 啦啦啦观看免费观看视频高清| 国产精品女同一区二区软件 | 国产高清有码在线观看视频| 日韩欧美精品v在线| 又粗又爽又猛毛片免费看| xxx96com| 黄色片一级片一级黄色片| 少妇的逼好多水| 国产一区二区三区视频了| 狂野欧美激情性xxxx| 88av欧美| 国产精品99久久99久久久不卡| 啦啦啦韩国在线观看视频| 欧美色视频一区免费| 97超视频在线观看视频| 亚洲成人精品中文字幕电影| 少妇熟女aⅴ在线视频| 性欧美人与动物交配| 久久久久久久精品吃奶| 女生性感内裤真人,穿戴方法视频| 亚洲乱码一区二区免费版| 日韩精品中文字幕看吧| 中文字幕久久专区| 少妇人妻精品综合一区二区 | 狂野欧美白嫩少妇大欣赏| 在线免费观看的www视频| 欧美区成人在线视频| 成人欧美大片| 嫁个100分男人电影在线观看| 少妇的逼好多水| 91麻豆精品激情在线观看国产| www日本在线高清视频| 欧美中文日本在线观看视频| 久久6这里有精品| 久久中文看片网| 12—13女人毛片做爰片一| 麻豆一二三区av精品| 亚洲av不卡在线观看| 男插女下体视频免费在线播放| 国产黄a三级三级三级人| 两个人的视频大全免费| 天天一区二区日本电影三级| www.熟女人妻精品国产| 欧美乱妇无乱码| 我要搜黄色片| 亚洲中文字幕日韩| 99久久精品一区二区三区| 精品久久久久久久久久久久久| bbb黄色大片| 嫁个100分男人电影在线观看| 国产av不卡久久| 久久精品国产自在天天线| 亚洲色图av天堂| tocl精华| 亚洲18禁久久av| 神马国产精品三级电影在线观看| 国语自产精品视频在线第100页| 免费一级毛片在线播放高清视频| 成熟少妇高潮喷水视频| 香蕉丝袜av| 宅男免费午夜| 亚洲自拍偷在线| 九九在线视频观看精品| 9191精品国产免费久久| 免费电影在线观看免费观看| 国产成+人综合+亚洲专区| 午夜激情福利司机影院| 精品乱码久久久久久99久播| 久久久久国产精品人妻aⅴ院| 熟女少妇亚洲综合色aaa.| 亚洲黑人精品在线| 12—13女人毛片做爰片一| 非洲黑人性xxxx精品又粗又长| 九九热线精品视视频播放| 成人国产综合亚洲| 熟女人妻精品中文字幕| 国产av不卡久久| 日韩中文字幕欧美一区二区| 99久久成人亚洲精品观看| 日韩欧美精品免费久久 | 3wmmmm亚洲av在线观看| 在线播放无遮挡| 国产单亲对白刺激| 午夜老司机福利剧场| 精品一区二区三区人妻视频| 欧美国产日韩亚洲一区| 亚洲熟妇中文字幕五十中出| 一区二区三区激情视频| 亚洲av五月六月丁香网| 国产色婷婷99| 蜜桃久久精品国产亚洲av| 久久久国产成人免费| 日韩欧美免费精品| 国产精品一区二区三区四区久久| 日本黄大片高清| 在线播放无遮挡| 啦啦啦观看免费观看视频高清| 欧美中文综合在线视频| 国产亚洲欧美98| 91在线精品国自产拍蜜月 | 2021天堂中文幕一二区在线观| 黑人欧美特级aaaaaa片| 日韩成人在线观看一区二区三区| 色精品久久人妻99蜜桃| 国产伦人伦偷精品视频| 亚洲av免费在线观看| 亚洲真实伦在线观看| 国产日本99.免费观看| 精品免费久久久久久久清纯| 琪琪午夜伦伦电影理论片6080| 婷婷六月久久综合丁香| 亚洲精品亚洲一区二区| 激情在线观看视频在线高清| 无人区码免费观看不卡| 18禁美女被吸乳视频| 五月玫瑰六月丁香| 97人妻精品一区二区三区麻豆| 麻豆久久精品国产亚洲av| 国产成人a区在线观看| 99riav亚洲国产免费| 听说在线观看完整版免费高清| 深夜精品福利| 动漫黄色视频在线观看| 高清在线国产一区| tocl精华| 日本精品一区二区三区蜜桃| 搡老岳熟女国产| 不卡一级毛片| 免费一级毛片在线播放高清视频| 亚洲一区高清亚洲精品| 最近在线观看免费完整版| 51午夜福利影视在线观看| 成人特级av手机在线观看| 国产私拍福利视频在线观看| 老汉色∧v一级毛片| 成人精品一区二区免费| 精品人妻一区二区三区麻豆 | 中文字幕人成人乱码亚洲影| 欧美高清成人免费视频www| 操出白浆在线播放| 性色av乱码一区二区三区2| 九色成人免费人妻av| 宅男免费午夜| 99热只有精品国产| 他把我摸到了高潮在线观看| 成人三级黄色视频| 一进一出抽搐gif免费好疼| 全区人妻精品视频| 熟女电影av网| 欧美精品啪啪一区二区三区| 99精品在免费线老司机午夜| 国产69精品久久久久777片| 国产美女午夜福利| 最近最新免费中文字幕在线| 亚洲七黄色美女视频| av中文乱码字幕在线| 91在线观看av| 窝窝影院91人妻| 欧美另类亚洲清纯唯美| 国产亚洲欧美98| 亚洲人成电影免费在线| 国产成人aa在线观看| 一本久久中文字幕| 精华霜和精华液先用哪个| 一a级毛片在线观看| 少妇人妻一区二区三区视频| 国内少妇人妻偷人精品xxx网站| 小说图片视频综合网站| 757午夜福利合集在线观看| 91九色精品人成在线观看| 一级毛片女人18水好多| 国产高清三级在线| 蜜桃亚洲精品一区二区三区| 精品熟女少妇八av免费久了| 少妇的丰满在线观看| 尤物成人国产欧美一区二区三区| 亚洲天堂国产精品一区在线| 午夜精品一区二区三区免费看| 午夜免费男女啪啪视频观看 | 亚洲精华国产精华精| 偷拍熟女少妇极品色| 1024手机看黄色片| 日本三级黄在线观看| 国产精品99久久久久久久久| 1000部很黄的大片| 成人国产一区最新在线观看| 亚洲人成网站在线播放欧美日韩| 国产亚洲精品一区二区www| 天堂av国产一区二区熟女人妻| 国产高清三级在线| 法律面前人人平等表现在哪些方面| 亚洲成a人片在线一区二区| e午夜精品久久久久久久| av黄色大香蕉| 在线看三级毛片| 色精品久久人妻99蜜桃| 亚洲无线在线观看| 香蕉av资源在线| 一级a爱片免费观看的视频| 午夜影院日韩av| av欧美777| 97超视频在线观看视频| 日本精品一区二区三区蜜桃| 亚洲精品国产精品久久久不卡| 有码 亚洲区| 在线观看美女被高潮喷水网站 | 国产免费男女视频| 日日夜夜操网爽| 婷婷丁香在线五月| 国产淫片久久久久久久久 | 国产激情欧美一区二区| 成人鲁丝片一二三区免费| 99riav亚洲国产免费| 亚洲五月天丁香| 亚洲精品国产精品久久久不卡| 哪里可以看免费的av片| 丁香欧美五月| 人妻夜夜爽99麻豆av| 亚洲精品亚洲一区二区| 波多野结衣巨乳人妻| 欧美国产日韩亚洲一区| 国产精品久久电影中文字幕| 欧美日韩福利视频一区二区| 欧美色欧美亚洲另类二区| 国模一区二区三区四区视频| 最近视频中文字幕2019在线8| 丁香欧美五月| 欧美日韩瑟瑟在线播放| 日日摸夜夜添夜夜添小说| 亚洲天堂国产精品一区在线| 一级黄片播放器| 国产单亲对白刺激| 亚洲国产中文字幕在线视频| 法律面前人人平等表现在哪些方面| 99久久成人亚洲精品观看| xxx96com| 小蜜桃在线观看免费完整版高清| 又黄又粗又硬又大视频| 变态另类成人亚洲欧美熟女| 成人av一区二区三区在线看| 国产精品久久久久久亚洲av鲁大| 久久人人精品亚洲av| 色噜噜av男人的天堂激情| 此物有八面人人有两片| 久久这里只有精品中国| 午夜福利欧美成人| 国产午夜福利久久久久久| 国产精品国产高清国产av| 免费在线观看影片大全网站| 午夜福利18| 搡老岳熟女国产| 精品乱码久久久久久99久播| 国产伦人伦偷精品视频| 久99久视频精品免费| 99在线视频只有这里精品首页| 久久久色成人| bbb黄色大片| 在线免费观看不下载黄p国产 | 99国产精品一区二区蜜桃av| 久久精品影院6| 九色成人免费人妻av| 一级作爱视频免费观看| 美女cb高潮喷水在线观看| 男人和女人高潮做爰伦理| 国产一区二区三区在线臀色熟女| 叶爱在线成人免费视频播放| 岛国在线免费视频观看| 久久性视频一级片| 国产精品美女特级片免费视频播放器| 亚洲男人的天堂狠狠| 国产精品美女特级片免费视频播放器| 在线免费观看不下载黄p国产 | 1024手机看黄色片| 久久久色成人| 欧美zozozo另类| 精品久久久久久久末码| 国产成人影院久久av| 久9热在线精品视频| 亚洲人与动物交配视频| 亚洲欧美激情综合另类| 亚洲男人的天堂狠狠| 国产精品一区二区免费欧美| 国产麻豆成人av免费视频| 欧美日韩一级在线毛片| 女人被狂操c到高潮| 丰满人妻熟妇乱又伦精品不卡| 男女那种视频在线观看| 亚洲av成人精品一区久久| 真人做人爱边吃奶动态| 久久人妻av系列| 久久久久久久亚洲中文字幕 | 18禁裸乳无遮挡免费网站照片| 一进一出抽搐动态| 国产单亲对白刺激| 国产av一区在线观看免费| 国内精品久久久久精免费| 一本综合久久免费| 精品一区二区三区人妻视频| 淫妇啪啪啪对白视频| 国产探花在线观看一区二区| 成人欧美大片| 18+在线观看网站| 网址你懂的国产日韩在线| ponron亚洲| 一进一出好大好爽视频| 久久亚洲真实| 国产精品99久久99久久久不卡| 91久久精品电影网| 一进一出抽搐gif免费好疼| 男女床上黄色一级片免费看| 在线观看午夜福利视频| 高清在线国产一区| 亚洲人成网站在线播放欧美日韩| 国产精品久久久人人做人人爽| 久久久色成人| 国产高清三级在线| 一级黄色大片毛片| 久久99热这里只有精品18| 母亲3免费完整高清在线观看| 午夜福利成人在线免费观看| 老熟妇乱子伦视频在线观看| 日韩成人在线观看一区二区三区| 久久久精品欧美日韩精品| 国产美女午夜福利| 人妻丰满熟妇av一区二区三区| 女人高潮潮喷娇喘18禁视频| 黄色视频,在线免费观看| 夜夜夜夜夜久久久久| 久久精品亚洲精品国产色婷小说| 少妇人妻一区二区三区视频| 91在线精品国自产拍蜜月 | 麻豆成人av在线观看| av天堂在线播放| 99久久99久久久精品蜜桃| tocl精华| 国产伦精品一区二区三区视频9 | 我的老师免费观看完整版| 亚洲欧美日韩高清在线视频| 国产精品野战在线观看| 亚洲成人精品中文字幕电影| 毛片女人毛片| 亚洲欧美日韩高清在线视频| 亚洲中文字幕日韩| 久久人人精品亚洲av| 国产欧美日韩精品亚洲av| 国产精华一区二区三区| 亚洲欧美精品综合久久99| 99久久九九国产精品国产免费| 亚洲国产精品sss在线观看| 国产精品久久久久久人妻精品电影| 黄色丝袜av网址大全| 老汉色∧v一级毛片| 日韩av在线大香蕉| 特级一级黄色大片| 色综合站精品国产| 两人在一起打扑克的视频| 国产伦一二天堂av在线观看| 无遮挡黄片免费观看| 舔av片在线| xxxwww97欧美| 久久久久国产精品人妻aⅴ院| 国产高清videossex| 日本黄色视频三级网站网址| a级毛片a级免费在线| www国产在线视频色| 麻豆成人av在线观看| 国产国拍精品亚洲av在线观看 | 黄色丝袜av网址大全| 精品一区二区三区视频在线观看免费| 国产主播在线观看一区二区| 精品乱码久久久久久99久播| 中亚洲国语对白在线视频| 中文字幕久久专区| 国产男靠女视频免费网站| 日韩av在线大香蕉| 午夜精品一区二区三区免费看| ponron亚洲| 制服丝袜大香蕉在线| 成人亚洲精品av一区二区| 又紧又爽又黄一区二区| 亚洲熟妇熟女久久| 桃色一区二区三区在线观看| 欧美激情久久久久久爽电影| 99热只有精品国产| 国产成人av教育| 手机成人av网站| 热99在线观看视频| 久9热在线精品视频| 欧美乱码精品一区二区三区| 欧美在线黄色| 亚洲欧美日韩高清专用| 在线观看舔阴道视频| 少妇丰满av| 久久久久久久精品吃奶| 国产在线精品亚洲第一网站| 精品人妻偷拍中文字幕| 一个人免费在线观看电影| 国产aⅴ精品一区二区三区波| 夜夜看夜夜爽夜夜摸| 叶爱在线成人免费视频播放| 国产成人啪精品午夜网站| 国语自产精品视频在线第100页| 久久99热这里只有精品18| 嫩草影院入口| 日韩大尺度精品在线看网址| 一区二区三区高清视频在线| 国产三级黄色录像| av黄色大香蕉| 大型黄色视频在线免费观看| 在线a可以看的网站| 日日干狠狠操夜夜爽| 日本三级黄在线观看| 成人av一区二区三区在线看| 亚洲av一区综合| 最近最新中文字幕大全电影3| 男人舔奶头视频| 中亚洲国语对白在线视频| 亚洲av中文字字幕乱码综合| 亚洲精品456在线播放app | 久久久久久久久久黄片| 国产美女午夜福利| 老司机在亚洲福利影院| 女生性感内裤真人,穿戴方法视频| 中出人妻视频一区二区| 三级国产精品欧美在线观看| 欧美激情在线99| 窝窝影院91人妻| 午夜a级毛片| 久久九九热精品免费| 一本综合久久免费| 久久精品国产亚洲av香蕉五月| 久9热在线精品视频| 久久精品国产亚洲av涩爱 | 中文字幕久久专区| 国产真人三级小视频在线观看| 给我免费播放毛片高清在线观看| 9191精品国产免费久久| 亚洲欧美日韩无卡精品| 午夜福利在线观看吧| 色综合婷婷激情| 91av网一区二区| 久久午夜亚洲精品久久| 亚洲国产精品合色在线| 国产成人福利小说| 欧美中文综合在线视频| 91字幕亚洲| 亚洲精品久久国产高清桃花| 国产午夜精品论理片| 国产私拍福利视频在线观看| xxxwww97欧美| 男女那种视频在线观看| 神马国产精品三级电影在线观看| 日本 欧美在线| 免费在线观看日本一区| 男女午夜视频在线观看| tocl精华| 内射极品少妇av片p| 噜噜噜噜噜久久久久久91| 真实男女啪啪啪动态图| 国产午夜精品久久久久久一区二区三区 | 国产91精品成人一区二区三区| 亚洲成a人片在线一区二区| 动漫黄色视频在线观看| 91久久精品国产一区二区成人 | 国内精品美女久久久久久| 岛国在线免费视频观看| 窝窝影院91人妻| 99久久精品热视频| 色综合欧美亚洲国产小说| 成人国产一区最新在线观看| 桃红色精品国产亚洲av| www日本在线高清视频| 90打野战视频偷拍视频| 中文字幕人妻熟人妻熟丝袜美 | 久久久久性生活片| 亚洲一区二区三区色噜噜| 亚洲aⅴ乱码一区二区在线播放| 91在线观看av| 哪里可以看免费的av片| 精品人妻一区二区三区麻豆 | 婷婷精品国产亚洲av| 免费av毛片视频| 91麻豆av在线| 无遮挡黄片免费观看| 熟女人妻精品中文字幕| 午夜两性在线视频| 色播亚洲综合网| 国产精品亚洲美女久久久| 天堂动漫精品| 一卡2卡三卡四卡精品乱码亚洲| 午夜福利在线在线| 亚洲欧美日韩高清专用| 又紧又爽又黄一区二区| 精品免费久久久久久久清纯| 窝窝影院91人妻| 欧美黄色片欧美黄色片| 99热这里只有精品一区| 亚洲自拍偷在线| 免费在线观看亚洲国产| 三级国产精品欧美在线观看| 亚洲精品久久国产高清桃花| 老熟妇乱子伦视频在线观看| 老汉色av国产亚洲站长工具| 日本三级黄在线观看| 午夜福利在线在线| 久9热在线精品视频| 日本成人三级电影网站| 搡老妇女老女人老熟妇| 成人精品一区二区免费| 天堂√8在线中文| 精品人妻一区二区三区麻豆 | 最近最新中文字幕大全电影3| 在线观看一区二区三区| 欧美色欧美亚洲另类二区| 悠悠久久av| 给我免费播放毛片高清在线观看| 在线a可以看的网站| 亚洲aⅴ乱码一区二区在线播放| 欧美乱码精品一区二区三区| 成人特级av手机在线观看| 亚洲欧美日韩高清在线视频| 丰满的人妻完整版| 色尼玛亚洲综合影院| 色噜噜av男人的天堂激情| 亚洲第一欧美日韩一区二区三区| 精品人妻偷拍中文字幕| 午夜精品在线福利| 在线a可以看的网站| 天美传媒精品一区二区| 亚洲国产日韩欧美精品在线观看 | 中出人妻视频一区二区| 最好的美女福利视频网| 一进一出抽搐动态| 亚洲 国产 在线| 真人一进一出gif抽搐免费| 亚洲 欧美 日韩 在线 免费| 亚洲精华国产精华精| 国产精品98久久久久久宅男小说| 天美传媒精品一区二区| 在线观看舔阴道视频| 中文字幕熟女人妻在线| 脱女人内裤的视频| 亚洲熟妇熟女久久| 欧美一级a爱片免费观看看| 国产精品免费一区二区三区在线| 内地一区二区视频在线| 亚洲自拍偷在线| 狂野欧美激情性xxxx| av中文乱码字幕在线| xxxwww97欧美| 波多野结衣高清无吗| 国产精品美女特级片免费视频播放器| 久久久久久人人人人人| 午夜精品久久久久久毛片777| 国产极品精品免费视频能看的| 亚洲五月婷婷丁香| 国产亚洲精品久久久久久毛片| 中亚洲国语对白在线视频| 国产精品1区2区在线观看.| 免费人成视频x8x8入口观看| 99国产精品一区二区三区| 国产精品精品国产色婷婷| 好看av亚洲va欧美ⅴa在| 90打野战视频偷拍视频| 成人欧美大片| 国产黄a三级三级三级人| 啦啦啦观看免费观看视频高清| 国产91精品成人一区二区三区| 级片在线观看| 黄色女人牲交| 色av中文字幕| 国产极品精品免费视频能看的| 色精品久久人妻99蜜桃| 老司机在亚洲福利影院| 日韩高清综合在线| 中文亚洲av片在线观看爽| 欧美日韩乱码在线| 熟妇人妻久久中文字幕3abv| 欧美日韩乱码在线| 欧美日韩国产亚洲二区| 深夜精品福利| 国产高清三级在线| 99热只有精品国产| 午夜精品一区二区三区免费看| 亚洲精品粉嫩美女一区| 国产精品野战在线观看|