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

    RDX及其衍生物高溫?zé)峤獾姆磻?yīng)分子動力學(xué)模擬

    2017-05-12 06:58:02彭莉娟王靜波李澤榮李象遠(yuǎn)
    物理化學(xué)學(xué)報 2017年4期
    關(guān)鍵詞:含碳單胞衍生物

    彭莉娟 姚 倩 王靜波 李澤榮 朱 權(quán),* 李象遠(yuǎn),*

    RDX及其衍生物高溫?zé)峤獾姆磻?yīng)分子動力學(xué)模擬

    彭莉娟1姚 倩2王靜波3李澤榮2朱 權(quán)3,*李象遠(yuǎn)3,*

    (1四川大學(xué)空天科學(xué)與工程學(xué)院,成都610065;2四川大學(xué)化學(xué)學(xué)院,成都610064;3四川大學(xué)化學(xué)工程學(xué)院,成都610065)

    本文采用ReaxFF反應(yīng)分子動力學(xué)方法,研究了RDX及其衍生物晶體在高溫條件下(2000、2500和3000 K)的熱分解機(jī)理以及主要產(chǎn)物隨時間的變化情況。結(jié)果表明:RDX及其衍生物晶體熱解的第一步反應(yīng)均為N―NO2鍵斷裂生成NO2分子,隨后反應(yīng)由于六元環(huán)上和側(cè)鏈基團(tuán)的不同導(dǎo)致側(cè)鏈脫除或開環(huán)反應(yīng)的順序不同。在這四種單胞體系熱解中,NO2和NO分子為共同的中間產(chǎn)物,形成之后迅速發(fā)生次級反應(yīng)并最終生成N2。各體系熱解終產(chǎn)物一致,均為N2、H2O和CO2,其中N2分子數(shù)最多,大于20個。由于原始分子結(jié)構(gòu)和組成中C/N比、H/O比的不同,各體系熱解后H2O和CO2分子數(shù)目相差較大。不同溫度下,各單胞體系熱分解生成的最大含碳團(tuán)簇中C原子數(shù)均較小。在進(jìn)一步超胞體系的模擬中,RDX和RDX-D2體系生成的含碳團(tuán)簇中C原子數(shù)分別達(dá)到約30和16個,遠(yuǎn)高于單胞模擬,且受溫度影響較大;而RDX-D1和RDX-D3單胞或超胞模擬結(jié)果相近,均未生成含碳團(tuán)簇,僅存在小分子含碳碎片。因此,初始分子的結(jié)構(gòu)和元素比對含碳團(tuán)簇的生成有明顯影響。

    RDX及其衍生物;高溫?zé)岱纸猓籖eaxFF;含碳團(tuán)簇;分子動力學(xué)

    1 引言

    含能材料是一類儲存高能量的特殊材料,在一定外界條件刺激下,自身能在短時間內(nèi)發(fā)生氧化還原反應(yīng),釋放大量能量,在軍用和民用領(lǐng)域廣泛應(yīng)用1。高溫?zé)岱纸膺^程和反應(yīng)機(jī)理的研究,不僅可提高對含能材料反應(yīng)特性的認(rèn)識,而且還可在原子分子水平上進(jìn)行改進(jìn)和設(shè)計,從而研制出更好性能的炸藥和推進(jìn)劑配方。

    黑索金(RDX)是一種綜合性能優(yōu)良的單質(zhì)炸藥和推進(jìn)劑。通常RDX存在α和β兩種晶型,常溫常壓下α-RDX性質(zhì)穩(wěn)定而β-RDX極不穩(wěn)定,因此α-RDX具有更高的實用價值2,本文研究均采用α晶型。文獻(xiàn)中,針對RDX的熱分解過程開展了大量的理論和實驗研究。Goddard等3采用從頭算方法研究了單分子RDX熱分解過程,提出了三種主要初始分解途徑:協(xié)同分解反應(yīng),N―NO2的斷裂反應(yīng)和連續(xù)HONO消去反應(yīng),其中連續(xù)HONO消去反應(yīng)被認(rèn)為是能量最低的反應(yīng)途徑。趙繼軍等4研究了溫度為3000 K下RDX的熱分解行為,表明熱解后的主要產(chǎn)物為N2、H2O和CO2。Schweigert5采用從頭算分子動力學(xué)研究氣相單分子RDX在高溫2000 K下的詳細(xì)裂解反應(yīng)機(jī)理。Irikura6采用密度泛函理論計算比較了RDX單分子和雙分子熱解反應(yīng)途徑,為凝聚相反應(yīng)機(jī)理的研究提供了理論依據(jù)。Behrens和Bulusu7從實驗上指出RDX的主要分解路徑為RDX分子去除HNO和兩個HONO分子生成1,3,5-三嗪(TAZ)。雖然,氣相單分子RDX熱解機(jī)理研究已經(jīng)比較成熟,但凝聚相RDX由于高溫?zé)峤庵蟹肿娱g復(fù)雜的相互作用,反應(yīng)機(jī)理研究仍不夠深入。

    ReaxFF反應(yīng)力場8-10通過鍵長、鍵級和鍵能之間關(guān)系,以及原子間電荷的動態(tài)傳遞來模擬原子間化學(xué)鍵的斷裂和生成,為研究較大體系的物理化學(xué)性質(zhì)提供了有效途徑。ReaxFF反應(yīng)力場已經(jīng)廣泛應(yīng)用于熱分解11-14、燃燒15,16和爆轟17,18等快速反應(yīng)過程。采用這一方法,Strachan等19研究了不同溫度和密度下RDX的熱分解過程,表明產(chǎn)物中CO和CO2分子數(shù)目很大程度上與系統(tǒng)密度有關(guān),而N2和H2O分子的數(shù)目幾乎不受密度的影響。Zhang等20研究了不同溫度和密度下奧克托金(HMX)和三氨基三硝基苯(TATB)晶體的熱分解過程,表明HMX的初始裂解為N―NO2鍵的斷裂反應(yīng)和HONO消去反應(yīng);TATB在低溫下幾乎沒有發(fā)現(xiàn)C―NO2鍵的斷裂,高溫下主要為C―NO2鍵斷裂反應(yīng)和HONO消去反應(yīng)。但在RDX的熱解模擬中,關(guān)于反應(yīng)物分子結(jié)構(gòu)和元素組成對裂解產(chǎn)物和反應(yīng)途徑的影響,一直缺乏系統(tǒng)的對比分析,難以對含能材料的結(jié)構(gòu)設(shè)計提供指導(dǎo)。

    本工作根據(jù)RDX分子,改變六元環(huán)上和側(cè)鏈的取代基團(tuán)設(shè)計衍生物分子,從而改變含能材料中碳、氮含量比,研究其對反應(yīng)過程中能量釋放和含碳團(tuán)簇形成過程的影響。通過RDX取代基團(tuán)的逐步替代設(shè)計了三種化合物,它們的結(jié)構(gòu)分別如圖1所示,這三種衍生物的分子式分別為:C2H5O6N7(RDX-D1)、C3H7O6N7(RDX-D2)、C2H4O8N8(RDXD3)。相對于RDX,RDX-D1采用NH取代RDX主環(huán)上的CH2,改變了主環(huán)結(jié)構(gòu),提高了N元素含量;相對于RDX-D1,RDX-D2用―CH3取代了主環(huán)上NH基團(tuán)中的氫原子,改變了C元素含量和側(cè)鏈結(jié)構(gòu);相對于RDX-D1,RDX-D3用―NO2取代了主環(huán)上NH基團(tuán)中的氫原子,進(jìn)一步提高了N元素含量和改變了側(cè)鏈結(jié)構(gòu)。本工作采用ReaxFFMD方法,在單胞水平上模擬研究了高溫?zé)峤獾某跏挤磻?yīng)過程和反應(yīng)產(chǎn)物隨時間變化,在超胞水平上深入分析了反應(yīng)過程中含碳團(tuán)簇的演化情況,從而探討分子元素組成和分子結(jié)構(gòu)對RDX類化合物高溫?zé)峤庑阅艿挠绊憽?/p>

    圖1 RDX及其衍生物的分子結(jié)構(gòu)Fig.1 Structures of the RDX molecule and its derivatives

    2 計算模型與方法

    2.1 初始構(gòu)型的構(gòu)建

    RDX模擬采用的初始體系為一個單胞,依據(jù)RDX的X射線衍射單晶數(shù)據(jù)21,運用Materials Studio(MS)程序包Visualizer模塊,構(gòu)建單胞體系,如圖2(a)所示,該單胞結(jié)構(gòu)包含8個RDX分子,共計168個原子,初始密度為1.86 g·cm-3。由于RDX-D1與RDX分子具有相似的分子結(jié)構(gòu),故采用和RDX相同的單晶參數(shù),采用MS程序包Visualizer模塊構(gòu)建單胞體系,密度為1.87 g· cm-3。依據(jù)肖鶴鳴等22采用密度泛函理論計算結(jié)果,RDX-D3密度選取為1.92 g·cm-3。RDX-D2與RDX-D3具有相似的分子結(jié)構(gòu),故選用相同的參數(shù)。各體系單胞結(jié)構(gòu)如圖2所示,其中灰色、白色、藍(lán)色和紅色原子分別表示碳原子、氫原子、氮原子和氧原子,模擬均采用周期性邊界條件。

    針對含碳團(tuán)簇生成過程模擬,對RDX及其衍生物分別搭建了(2×2×1)超胞結(jié)構(gòu)。RDX超胞結(jié)構(gòu)共含有32個分子,共計672個原子;衍生物RDX-D1,RDX-D2和RDX-D3超胞結(jié)構(gòu)所含原子個數(shù)分別為640個、736個和704個。

    2.2 計算方法

    分子動力學(xué)模擬采用ReaxFF力場20,選取Velocity Verlet算法,采用ReaxFF-MD方法分別模擬各晶胞體系的高溫?zé)崃呀膺^程。在ReaxFF反應(yīng)力場中,各初始晶胞體系分別在低溫10 K下進(jìn)行幾何優(yōu)化,獲得最小能量結(jié)構(gòu),然后在相同溫度采用等溫等容(NVT)系綜進(jìn)行弛豫,使得模擬體系達(dá)到平衡狀態(tài),弛豫時間為5 ps,時間步長選取為0.1 fs。其后,各體系得到的平衡結(jié)構(gòu)分別在2000、2500、3000 K三個溫度下采取NVT系綜進(jìn)行高溫分解模擬,模擬時間步長選取0.1 fs,總模擬時間為100 ps。由于升溫過程可能會導(dǎo)致化學(xué)鍵的斷裂,因此,本文不考慮升溫過程的模擬。在所有的NVT-MD模擬中,各系統(tǒng)溫度均用Berendsen方法23進(jìn)行調(diào)節(jié),溫度調(diào)整參數(shù)選取100 fs,使溫度在設(shè)定值附近波動。當(dāng)任意原子對的鍵級大于0.3時,判定新化學(xué)鍵形成。

    3 結(jié)果與討論

    3.1 初始反應(yīng)途徑

    在3000 K條件下對四個單胞體系的熱分解進(jìn)行模擬,通過分析原子成鍵情況及各體系在不同時刻下的裂解產(chǎn)物,最終得到高溫?zé)峤獾某跏挤磻?yīng)途徑,如圖3所示。RDX及其三種衍生物熱解過程,第一步反應(yīng)均為N―NO2鍵發(fā)生斷裂生成NO2分子,表明主環(huán)和側(cè)鏈基團(tuán)的變化并不影響側(cè)鏈硝基斷鍵的初始反應(yīng)步驟。RDX分子第一步首先發(fā)生N―NO2鍵的斷裂,從而生成大量的NO2,與Botcher等在實驗中觀察到RDX快速熱解初期生成大量NO2一致24,25。在隨后反應(yīng)中,RDX主要發(fā)生剩余N―NO2鍵斷裂反應(yīng)和母體環(huán)斷裂反應(yīng),這一結(jié)果與Goddard等3采用密度泛函理論計算得到的RDX中N―N鍵斷裂(解離能約83.68 kJ·mol-1)優(yōu)先于協(xié)同開環(huán)反應(yīng)的結(jié)論一致。

    圖2 RDX(a),RDX-D1(b),RDX-D2(c)和RDX-D3(d)單胞結(jié)構(gòu)示意圖Fig.2 Unit cell structure of RDX(a),RDX-D1(b),RDX-D2(c)and RDX-D3(d) color online

    圖3 RDX及其衍生物的初始反應(yīng)途徑Fig.3 Initial reaction path for RDX and its derivatives

    RDX由于分子本身的對稱性,側(cè)鏈第一和第二個NO2基的消去無選擇性,但在RDX-D1中,由于六元環(huán)上的CH2基團(tuán)被NH基團(tuán)取代,在熱解過程中優(yōu)先脫去NH基團(tuán)鄰位氮原子上的兩個側(cè)鏈NO2基,其后,不同于RDX的開環(huán)反應(yīng),RDX-D1發(fā)生HONO脫去反應(yīng),之后再開環(huán)。在RDX-D1基礎(chǔ)上將環(huán)上NH中的氫原子用甲基取代,RDXD2中兩個NO2基脫除也優(yōu)先發(fā)生于N―CH3鄰位氮原子上,其后發(fā)生的是甲基脫除反應(yīng),然后發(fā)生開環(huán)反應(yīng)。因此,與RDX-D1相比,CH3基的引入對鄰位N―NO2的斷鍵影響不大,而對對位NO2的脫除影響較大,使得后續(xù)的反應(yīng)途徑完全不同。將RDX-D2中的甲基用硝基取代,RDX-D3中相鄰的兩個NO2基優(yōu)先脫除,六元環(huán)上形成氮氮雙鍵,其后雙鍵間位上的NO2基發(fā)生脫除,然后發(fā)生開環(huán)反應(yīng)。

    在RDX及其衍生物晶胞熱解過程中主要觀察到三種初始分解機(jī)理:N―NO2鍵的斷裂反應(yīng),HONO消去反應(yīng)和分子主環(huán)的開環(huán)反應(yīng),與Goddard等采用密度泛函理論計算的結(jié)果一致3,22。對比這四種單胞,第一步反應(yīng)均為N―NO2鍵斷裂,表明這一通道是能量最低的反應(yīng)通道,與Imamura等研究的其它硝胺化合物的熱分解機(jī)理一致26-29。然后才出現(xiàn)HONO消去、甲基取代基脫除和開環(huán)反應(yīng),進(jìn)一步生成小分子碎片。在四種體系的高溫?zé)峤庵?,?cè)鏈取代基團(tuán)對N―NO2第一步斷鍵無影響,主要改變了隨后的側(cè)鏈脫除或開環(huán)反應(yīng)的順序。

    3.2 溫度對勢能和總分子碎片數(shù)的影響

    在不同溫度條件下,我們模擬了各單胞體系(經(jīng)過5 ps初始弛豫)勢能隨時間的變化,如圖4所示。在初始熱解時,系統(tǒng)迅速吸收熱量發(fā)生初級熱解反應(yīng),使得系統(tǒng)內(nèi)勢能升高,隨后的次級反應(yīng)生成一系列中間及最終產(chǎn)物,并放出大量的熱,因而,系統(tǒng)勢能發(fā)生明顯下降。溫度越高,系統(tǒng)達(dá)到化學(xué)反應(yīng)動態(tài)平衡所需時間越短,熱解反應(yīng)越完全,生成的產(chǎn)物越多,勢能下降也越多。因而,各體系勢能的改變與其熱分解反應(yīng)過程中的化學(xué)鍵斷裂和生成緊密相關(guān)30。

    從圖4(a)可見,系統(tǒng)溫度越低,初始反應(yīng)越慢,吸熱階段越長;隨著系統(tǒng)溫度的升高,各體系分子初始反應(yīng)加劇,從而導(dǎo)致次級反應(yīng)越早開始,系統(tǒng)達(dá)到化學(xué)反應(yīng)動態(tài)平衡的時間縮短。當(dāng)溫度為3000 K時,各晶胞體系勢能在20 ps后相對穩(wěn)定;而當(dāng)系統(tǒng)處于2000 K的較低溫度時,各體系勢能在100 ps內(nèi)依然緩慢減少,需要更長時間才能達(dá)到穩(wěn)定,故在低溫條件達(dá)到化學(xué)反應(yīng)平衡的時間更長。

    圖5為各單胞體系在不同溫度下熱解過程中總分子碎片數(shù)隨時間的變化情況。從圖5(a)可見,在高溫初始裂解階段,溫度越高RDX熱解釋放能量的速率越快,物種數(shù)增加越快、產(chǎn)生的總分子碎片數(shù)越多,在10 ps左右達(dá)到峰值后,這些快速熱解產(chǎn)生的分子碎片會進(jìn)一步轉(zhuǎn)化,但整體數(shù)量維持在50左右。其他衍生物具有類似的變化規(guī)律,各晶胞體系的總分子碎片數(shù)在前10 ps內(nèi)快速增加,隨后變化相對比較緩慢,且溫度越高,總分子碎片數(shù)增加速度越大。在50 ps后,各晶胞體系的總分子碎片數(shù)在3000 K溫度條件下基本維持恒定,而在較低的溫度2000 K條件下仍緩慢增加。綜上,隨熱解溫度升高,各晶胞體系產(chǎn)生的總分子碎片數(shù)變多,總分子碎片數(shù)增加到最大值的時間縮短,體系釋放的能量也更多。

    3.3 主產(chǎn)物及中間產(chǎn)物分析

    圖4 四個體系的勢能隨時間變化的關(guān)系Fig.4 Time evolution of potential energy for the four materials

    圖5 四個體系的熱解總分子碎片數(shù)隨時間變化關(guān)系Fig.5 Time evolution of total fragments for the four materials

    RDX及其衍生物在3000 K時的主要分解產(chǎn)物N2、H2O和CO2及中間產(chǎn)物NO和NO2分子數(shù)目隨時間的變化如圖6所示。在各體系主產(chǎn)物中,N2是生成數(shù)目最多的產(chǎn)物,主要來自次級反應(yīng)中NO2和NO的轉(zhuǎn)化,其主要機(jī)理為NO2或NO與含碳分子碎片反應(yīng)生成N2。在前20 ps內(nèi),NO2和NO均迅速增大然后減少,此過程中N2數(shù)目一直迅速增加到約20,說明消耗掉的大部分NO2和NO生成了N2。在這四種不同體系中,主產(chǎn)物分子數(shù)目相對較多的還有H2O和CO2,這三種主要產(chǎn)物數(shù)目的大小順序為:N2>H2O>CO2。

    在四種不同的晶胞體系中,主要分解產(chǎn)物H2O和CO2分子數(shù)目相差較大。表1列出了各體系分子的元素組成比以及主產(chǎn)物CO2和H2O分子的最終數(shù)目。從表1可以清晰看出RDX和RDX-D2分子中C/N比大于等于0.5,高于RDX-D1和RDX-D3分子中的C/N比;RDX和RDX-D2體系中碳原子含量相多較高,熱解穩(wěn)定后生成的CO2分子數(shù)均小于10,這說明體系中的碳原子的最終轉(zhuǎn)化途徑并不以生成穩(wěn)定產(chǎn)物CO2為主,而是轉(zhuǎn)變成了其它的含碳分子碎片。與此對應(yīng),RDX和RDX-D2由于其分子構(gòu)成的H/O元素比(≥1)高于另外兩種衍生物,熱解后生成H2O分子數(shù)大于20,H2O的生成消耗了大量的氧,同時也抑制了CO2的生成。說明在初始的分子設(shè)計中C/N比和H/O比,對熱解反應(yīng)路徑和最終產(chǎn)物分布起決定作用。

    仔細(xì)觀察圖6各體系中N2和CO2分子數(shù)目,發(fā)現(xiàn)它們幾乎在同一時間達(dá)到最大,這與N2和CO2的生成途徑有關(guān);反應(yīng)早期的主要步驟是N―NO2的斷裂生成NO2,所以剛開始體系中NO2數(shù)量最多,隨后NO2進(jìn)一步反應(yīng)生成N2,而NO2中的O則與中間物CO反應(yīng)生成CO2,因而N2和CO2達(dá)到最大值的時間一致。同時還表明,高溫下N―N兩原子主要以N2的形式存在,與Strachan等10,19的研究結(jié)果一致。此外,在反應(yīng)過程中還觀察到一些數(shù)量很少基團(tuán)出現(xiàn),如:H、HO、CO、HNO3、N2O和HONO等基團(tuán),這些小分子基團(tuán)相互聚合或反應(yīng)生成其它產(chǎn)物;因其數(shù)量較少,故在本文沒有對其進(jìn)行詳細(xì)討論。

    由于圖6很難清晰展示中間體NO2和NO的變化,而中間體的變化影響著整個熱分解反應(yīng)過程,故研究了前20 ps中間體NO2和NO隨時間的變化,如圖7所示。RDX-D2中的NO2峰值最高,RDX-D3中的NO2峰值最低,RDX和RDX-D2中的NO2峰值(>13)略高于RDX-D1和RDX-D3的NO2峰值(〈10)。從圖7(a)中可以看出,NO2的數(shù)量在前2.5 ps內(nèi)迅速達(dá)到峰值,表明四種物質(zhì)首先發(fā)生斷裂的鍵是N―NO2鍵,這與圖3中研究表明高溫下N―NO2鍵的斷裂是各體系最主要的初始反應(yīng)途徑相一致。隨后NO2的數(shù)量開始減少,這是由于中間產(chǎn)物NO2參與次級反應(yīng)不斷被消耗,生成NO,NO3和HONO等中間產(chǎn)物。

    圖6 四個體系的主要產(chǎn)物分子(N2,H2O,CO2,NO2,NO)數(shù)目隨時間變化關(guān)系Fig.6 Time evolution of key produced molecules(N2,H2O,CO2,NO2,NO)for the four materials

    同時,NO也是很活躍的中間產(chǎn)物,從圖7(b)可以看出,NO和NO2的分布曲線具有相似的演化特征,均是先增大后減小的過程,但各體系NO達(dá)到峰值前的增長速率以及隨后的衰減速率相比NO2較小。在整個反應(yīng)過程中,NO數(shù)量整體上明顯低于NO2,且NO分子數(shù)達(dá)到最大值的時間相對NO2滯后一些,這是因為后續(xù)反應(yīng)消耗掉一部分NO2生成NO。NO的產(chǎn)生主要有三種途徑,(1)N―NO2→N―ONO,然后NO―NO單鍵斷裂生成NO;(2)氫原子轉(zhuǎn)移反應(yīng)生成的HONO,隨后分解為HO和NO分子;(3)由NO2和NO3等相互反應(yīng)生成NO。在20 ps后,NO2和NO的數(shù)量很少,幾乎被完全消耗。

    表1 各體系的C/N和H/O比值及主產(chǎn)物CO2和H2O的生成量Table 1 C/N and H/O ratios of each system and the amount of the main products of CO2and H2O

    3.4 大分子含碳團(tuán)簇的形成

    實驗研究表明,富碳高能材料在爆轟過程中能形成大分子含碳團(tuán)簇30-34。大分子含碳團(tuán)簇的出現(xiàn)是深入了解含能材料性質(zhì)的重要方面,目前采用分子動力學(xué)模擬,已經(jīng)開展含能材料熱解過程含碳團(tuán)簇形成動力學(xué)機(jī)制的研究35。為了深入了解含碳團(tuán)簇形成過程,首先模擬了各單胞體系高溫?zé)峤膺^程最大含碳團(tuán)簇C原子數(shù)隨時間的演變,如圖8(a,b,c)所示。本文中我們將C原子數(shù)大于8定義為含碳團(tuán)簇,隨著溫度升高,RDX和RDX-D2單胞體系熱解過程中最大含碳團(tuán)簇C原子數(shù)有所增加,在3000 K時含碳團(tuán)簇最大碳原子數(shù)可達(dá)19。而在RDX-D1和RDX-D3單胞體系熱解過程中,均未生成含碳團(tuán)簇,只存在小分子含碳碎片(C原子數(shù)小于等于3),且在不同溫度下幾乎保持不變。由于大分子含碳團(tuán)簇的生成受初始模型的影響較大,因此,為了正確考慮其形成過程,我們進(jìn)一步采用超胞體系進(jìn)行模擬,如圖8(d)所示。

    圖7 四個體系熱解過程中NO2(a)和NO(b)的分子數(shù)量變化關(guān)系Fig.7 Quantities of NO2(a)and NO(b)for the thermal pyrolysis of four systems

    圖8 不同溫度下各體系單胞(a-c)和超胞(d)熱解過程中最大含碳團(tuán)簇所含C原子數(shù)隨時間變化關(guān)系Fig.8 Time evolution of C atomic number in the maximum molecule formed in the thermal pyrolysis of unit cell(a-c)and super cell(d)under different temperatures

    對比圖8(c,d)單胞和多胞體系可以看出,RDX和RDX-D2超胞體系中生成的最大含碳團(tuán)簇C原子數(shù)增加迅速,最終C原子數(shù)目分別達(dá)到約30和16個;RDX-D1和RDX-D3單胞和超胞體系的模擬結(jié)果相近,均未生成大分子含碳團(tuán)簇,生成的小分子含碳碎片中C原子數(shù)均小于等于3??紤]初始分子結(jié)構(gòu)中,RDX和RDX-D2分子中原子比C/N≥0.5,而RDX-D1和RDX-D3分子中C/N〈0.5,說明:當(dāng)分子中C/N≥0.5時,生成的最大含碳團(tuán)簇C原子數(shù)與模型大小有關(guān);而對C/N〈0.5的體系,含碳團(tuán)簇的形成過程基本不受模型大小的影響。

    因而,初始設(shè)計的分子中碳、氮含量比和反應(yīng)物分子結(jié)構(gòu)對RDX及其衍生物中最大含碳團(tuán)簇的形成過程影響較大,對最大含碳團(tuán)簇中C原子數(shù)目多少也影響顯著,且C/N比越大越容易形成大分子含碳團(tuán)簇。根據(jù)上述模擬結(jié)果,在熱分解的初級反應(yīng)中,環(huán)上脫落NO2、NO和CH3等自由基后,剩余部分發(fā)生聚合形成大分子含碳團(tuán)簇,這與實驗和理論上研究硝胺類化合物的高溫?zé)峤膺^程是一致的20,36。

    4 結(jié)論

    采用ReaxFF反應(yīng)力場分子動力學(xué)方法,模擬研究了RDX及其衍生物單胞和超胞在不同溫度下的熱分解行為。研究表明,這四種晶體高溫?zé)峤獾牡谝徊骄鶠镹―NO2鍵的斷裂生成NO2,然后發(fā)生取代基的脫去和環(huán)的斷裂反應(yīng)。RDX在斷裂兩個硝基后就發(fā)生開環(huán)反應(yīng),而RDX的衍生物在斷裂兩個硝基后并沒有發(fā)生開環(huán)反應(yīng),而是先發(fā)生取代基的脫去反應(yīng)再發(fā)生開環(huán)反應(yīng),這是由于各體系分子六元環(huán)上基團(tuán)和側(cè)鏈基團(tuán)的改變使得初始反應(yīng)途徑不同。

    通過對各體系熱分解過程主要產(chǎn)物和中間產(chǎn)物進(jìn)行詳細(xì)分析,發(fā)現(xiàn)NO2和NO是裂解過程中十分重要的中間產(chǎn)物。NO2主要來自N―NO2中的N―N鍵的斷裂產(chǎn)生;NO主要來自N―NO2→N―ONO,然后NO―NO單鍵斷裂生成NO,中間產(chǎn)物HONO分解為HO和NO以及次級反應(yīng)生成NO分子。熱分解后最穩(wěn)定的氣相產(chǎn)物為N2、H2O和CO2,其中分子數(shù)最多的是N2(>20),主要來源于NO2和NO參與的次級反應(yīng)中。各分子體系中的C/ N比和H/O比不同,使得主要分解產(chǎn)物H2O和CO2分子數(shù)目相差較大。C/N和H/O比越大的體系,熱解后生成的H2O分子數(shù)越多,H2O的生成消耗了大量的氧,故生成穩(wěn)定產(chǎn)物CO2分子數(shù)越少,生成含碳團(tuán)簇分子量越大。模型分子設(shè)計中C/N比和H/O比,對于最終熱解過程和產(chǎn)物分布起決定作用。

    單胞體系模擬分析表明:溫度升高有利于RDX和RDX-D2單胞體系熱解中最大含碳團(tuán)簇的形成,在3000 K時初步形成大分子含碳團(tuán)簇(C原子數(shù)約為10)。進(jìn)一步超胞體系模擬結(jié)果表明:相對于單胞模擬,RDX和RDX-D2超胞熱解過程中最大含碳團(tuán)簇C原子數(shù)顯著增加(約為30和16),初始模型大小和反應(yīng)溫度對含碳團(tuán)簇的生成影響明顯。而RDX和RDX-D2體系中的C/N比明顯高于RDX-D1和RDX-D3體系,對于后兩者的模擬,單胞和超胞模擬結(jié)果相近,均未生成含碳團(tuán)簇,僅存在小分子含碳碎片(碳原子數(shù)小于等于3)。含碳團(tuán)簇的大量生成,不利于含能材料能量的釋放,故分子設(shè)計中減小C/N比能很好地改善含能材料的熱解性能。

    (1) Shaw,R.W.;Brill,T.B.;Thompson,D.L.Overviews of Recent Research on Energetic Materials;World Scientific:NJ,2005, Vol.16.

    (2) Teipel,U.Energetic Materials:Particle Processing and Characterization;Wiley-VCH:Weinheim,Germany,2005.

    (3) Chakraborty,D.;Muller,R.P.;Dasgupta,S.;Goddard,W.A.J. Phys.Chem.A 2000,104,2261.doi:10.1021/jp9936953

    (5) Schweigert,I.V.J.Phys.Chem.A 2015,119,2747. doi:10.1021/jp510034p

    (6) Irikura,K.K.J.Phys.Chem.A 2013,117,2233.doi:10.1021/ jp310247z

    (7) Behrens,R.,Jr.;Bulusu,S.J.Phys.Chem.1992,96,8877. doi:10.1063/1.462245

    (8) Van Duin,A.C.;Dasgupta,S.;Lorant,F.;Goddard,W.A.J. Phys.Chem.A 2001,105,9396.doi:10.1021/jp004368u

    (9) Cohen,R.;Zeiri,Y.;Wurzberg,E.;Kosloff,R.J.Phys.Chem.A 2007,111,11074.doi:10.1021/jp072121s

    (10) Strachan,A.;van Duin,A.C.;Chakraborty,D.;Dasgupta,S.; Goddard,W.A.,III Phys.Rev.Lett.2003,91,098301. doi:10.1103/PhysRevLett.91.098301

    (11) Zhou,T.T.;Shi,Y.D.;Huang,F.L.Acta Phys.-Chim.Sin. 2012,28(11),2605.[周婷婷,石一丁,黃風(fēng)雷.物理化學(xué)學(xué)報, 2012,28(11),2605.].doi:10.3866/PKU.WHXB201208031

    (12) Zhang,L.;Chen,L.;Wang,C.;Wu,J.Y.Acta Phys.-Chim.Sin. 2013,29(6),1145.[張 力,陳 朗,王 晨,伍俊英.物理化學(xué)學(xué)報,2013,29(6),1145.].doi:10.3866/PKU. WHXB201303221

    (13) Rom,N.;Zybin,S.V.;van Duin,A.C.;Goddard,W.A.,III.; Zeiri,Y.;Katz,G.;Kosloff,R.J.Phys.Chem.A 2011,115, 10181.doi:10.1021/jp202059v

    (14) Liu,H.;Dong,X.;He,Y.H.Acta Phys.-Chim.Sin.2014,30, 232.[劉 海,董 曉,何遠(yuǎn)航.物理化學(xué)學(xué)報,2014,30,232.] doi:10.3866/PKU.WHXB201312101

    (15) Qian,H.J.;van Duin,A.C.;Morokuma,K.;Irle,S.J.Chem. Theory Comput.2011,7,2040.doi:10.1021/ct200197v

    (16) Liu,L.;Bai,C.;Sun,H.;Goddard,W.A.,III.J.Phys.Chem.A 2011,115,4941.doi:10.1021/jp110435p

    (17)Ge,N.N.;Wei,Y.K.;Ji,G.F.;Chen,X.R.;Zhao,F.;Wei,D. Q.J.Phys.Chem.B 2012,116,13696.doi:10.1021/jp309120t

    (18) Zhang,L.Z.;Zybin,S.V.;Van Duin,A.C.T.;Goddard,W.A., III.J.Energ.Mater.2010,28,92.doi:10.1080/ 07370652.2010.504682

    (19) Strachan,A.;Kober,E.M.;van Duin,A.C.;Oxgaard,J.; Goddard,W.A.,III.J.Chem.Phys.2005,122,054502. doi:10.1063/1.1831277

    (20) Zhang,L.;Zybin,S.V.;van Duin,A.C.;Dasgupta,S.; Goddard,W.A.,III;Kober,E.M.J.Phys.Chem.A 2009,113, 10619.doi:10.1021/jp901353a

    (21) Hakey,P.;Ouellette,W.;Zubieta,J.;Korter,T.Acta Crystallographica Section E:Structure Reports Online 2008,64 (8),o1428.doi:10.1107/S1600536808019727

    (22) Qiu,L.;Gong,X.;Ju,X.;Xiao,H.Sci.China Ser.B:Chem. 2008,51,1231.doi:10.1007/s11426-008-0141-1

    (23) Berendsen,H.J.;Postma,J.V.;van Gunsteren,W.F.;DiNola, A.R.H.J.;Haak,J.R.J.Chem.Phys.1984,81,3684. doi:10.1063/1.448118

    (24) Botcher,T.R.;Wight,C.A.J.Phys.Chem.1994,98,5441. doi:10.1021/j100072a009

    (25) Oyumi,Y.;Brill,T.B.Combust.Flame 1985,62,213. doi:10.1016/0010-2180(85)90147-6

    (26)Kohno,Y.;Ueda,K.;Imamura,A.J.Phys.Chem.1996,100, 4701.doi:10.1021/jp9503223

    (27)Prabhakaran,K.V.;Bhide,N.M.;Kurian,E.M.Thermochim. Acta 1995,249,249.doi:10.1016/0040-6031(95)90704-1

    (28) Oxley,J.C.;Kooh,A.B.;Szekeres,R.;Zheng,W.J.Phys. Chem.1994,98,7004.doi:10.1021/j100079a019

    (29) Brill,T.B.;Oyumi,Y.J.Phys.Chem.1986,90,6848. doi:10.1021/j100284a028

    (30) Mironov,E.V.;Petrov,E.A.;Korets,A.Y.Combust.Explos. Shock Waves.2004,40,473.doi:10.1023/B: CESW.0000033571.82326.6a

    (31)Anisichkin,V.F.;Dolgushin,D.S.;Petrov,E.A.Combust. Explos.Shock Waves.1995,31,106.doi:10.1007/BF00755966

    (32) Titov,V.M.;Anisichkin,V.F.;Mal′kov,I.Y.Combust.Explos. Shock Waves.1989,25,372.doi:10.1007/BF00788819

    (33) Baidakova,M.J.Phys.D:Appl.Phys.2007,40,6300. doi:10.1088/0022-3727/40/20/S14

    (34) Mal'Kov,I.Y.;Filatov,L.I.;Titov,V.M.;Litvinov,B.V.; Chuvilin,A.L.;Teslenko,T.S.Combust.Explos.Shock Waves. 1993,29,542.doi:10.1007/BF00782983

    (35) Chevrot,G.;Sollier,A.;Pineau,N.J.Chem.Phys.2012,136, 084506.doi:10.1063/1.3686750

    (36) Hobbs,M.L.;Kaneshige,M.J.;Gilbert,D.W.;Marley,S.K.; Todd,S.N.J.Phys.Chem.A 2009,113,10474.doi:10.1021/ jp906134f

    Pyrolysis of RDX and Its Derivatives via Reactive Molecular Dynamics Simulations

    PENG Li-Juan1YAO Qian2WANG Jing-Bo3LI Ze-Rong2ZHU Quan3,*LI Xiang-Yuan3,*
    (1School of Aeronautics&Astronautics,Sichuan University,Chengdu 610065,P.R.China; 2College of Chemistry,Sichuan University,Chengdu 610064,P.R.China; 3School of Chemical Engineering,Sichuan University,Chengdu 610065,P.R.China)

    The thermal decomposition mechanisms of cyclotrimethylene trinitramine(RDX)and its derivatives are investigated using ReaxFF reactive molecular dynamics simulation at high temperatures(2000,2500,and 3000 K).It is shown that the first pyrolysis step of RDX and its derivatives is the N―NO2homolytic cleavage to produce NO2;however,the subsequent reaction mechanism is completely different due to the different sixmembered rings and side chain groups.In these four model systems,NO2and NO molecules are common intermediates in the thermal decomposition process,which eventually transform to N2in the following reactions. The most stable products are N2,H2O,and CO2in the thermal process,among which N2has the maximum molecular number(more than 20).The numbers of cracked H2O and CO2molecules are very different in the four model systems due to the different C/N and H/O ratios.At different temperatures,the maximum numbers of carbon atoms in the carbon clusters are all small in the unit cell simulations of the four systems.In the further super cell simulation for RDX and RDX-D2 the numbers of carbon atoms reach about 30 and 16,respectively, which are higher than those from the unit cell simulation.In both the unit and super cell simulations for RDX-D1 and RDX-D3,the carbon cluster cannot be formed and there exist only small carbon molecular fragments. Therefore,the initial molecular structure and elemental ratio have a large effect on the production of carbon clusters.

    RDX and its derivatives;Thermal decomposition;ReaxFF;Carbon clusters;Molecular dynamics

    O643

    Xu,J.;Zhao,J.;Sun,L.Mol.Simul.2008,34,961.

    10.1080/ 08927020802162892

    doi:10.3866/PKU.WHXB201701161

    Received:October 22,2016;Revised:January 16,2017;Published online:January 16,2017.

    *Corresponding authors.ZHU Quan,Email:qzhu@scu.edu.cn;Tel:+86-28-85402951.LI Xiang-Yuan,Email:xyli@scu.edu.cn.

    The project was supported by the National Natural Science Foundation of China(91441132)and Program for New Century Excellent Talents in University of Ministry of Education of China(NCET-13-0398).國家自然科學(xué)基金(91441132)和教育部新世紀(jì)優(yōu)秀人才支持計劃(NCET-13-0398)資助項目

    猜你喜歡
    含碳單胞衍生物
    基于NURBS的點陣材料參數(shù)化建模方法
    復(fù)合材料周期結(jié)構(gòu)數(shù)學(xué)均勻化方法的一種新型單胞邊界條件
    考慮界面層影響的三維機(jī)織復(fù)合材料單胞模型研究
    黑龍江省造林樹種含碳率與土壤性質(zhì)研究
    森林工程(2018年4期)2018-08-04 03:23:20
    烴的含氧衍生物知識鏈接
    新型螺雙二氫茚二酚衍生物的合成
    鐵礦含碳球團(tuán)還原后性能的實驗研究
    Xanomeline新型衍生物SBG-PK-014促進(jìn)APPsw的α-剪切
    基于B-H鍵的活化對含B-C、B-Cl、B-P鍵的碳硼烷硼端衍生物的合成與表征
    內(nèi)蒙古某含碳低品位硫化鉛鋅礦石選礦試驗
    金屬礦山(2013年6期)2013-03-11 16:53:57
    人人澡人人妻人| 99riav亚洲国产免费| 1024香蕉在线观看| 亚洲精品在线美女| 精品久久久精品久久久| 91精品国产国语对白视频| www.精华液| 热99re8久久精品国产| 亚洲av第一区精品v没综合| 大香蕉久久网| 国产欧美日韩一区二区三| 午夜免费成人在线视频| 亚洲一区中文字幕在线| 黄色怎么调成土黄色| 成年女人毛片免费观看观看9 | 国产成人精品久久二区二区91| 99精国产麻豆久久婷婷| 久久精品熟女亚洲av麻豆精品| 亚洲国产毛片av蜜桃av| 建设人人有责人人尽责人人享有的| 久久国产精品影院| 亚洲专区国产一区二区| 99国产精品一区二区三区| 色播在线永久视频| 亚洲色图综合在线观看| 少妇猛男粗大的猛烈进出视频| 91老司机精品| 91麻豆av在线| 无人区码免费观看不卡| 最近最新免费中文字幕在线| 国产男女超爽视频在线观看| 国产欧美日韩精品亚洲av| 国产av一区二区精品久久| 亚洲欧美色中文字幕在线| 欧美av亚洲av综合av国产av| 激情视频va一区二区三区| 女性被躁到高潮视频| 色婷婷久久久亚洲欧美| 亚洲一区二区三区不卡视频| 久久影院123| 亚洲一区高清亚洲精品| 久久精品91无色码中文字幕| 亚洲av日韩在线播放| cao死你这个sao货| 午夜福利免费观看在线| 国产午夜精品久久久久久| 黄色女人牲交| 午夜福利视频在线观看免费| 嫩草影视91久久| 一级黄色大片毛片| 我的亚洲天堂| 日韩欧美在线二视频 | 无限看片的www在线观看| 搡老乐熟女国产| 黑人猛操日本美女一级片| 国产男女内射视频| 中文字幕色久视频| 99久久国产精品久久久| 午夜老司机福利片| 国产精品98久久久久久宅男小说| 日韩大码丰满熟妇| 亚洲精品自拍成人| 国产亚洲精品久久久久5区| 香蕉丝袜av| 亚洲熟妇熟女久久| 人妻 亚洲 视频| 精品一品国产午夜福利视频| 精品亚洲成国产av| 999精品在线视频| 亚洲成人免费电影在线观看| 老司机亚洲免费影院| 国产激情欧美一区二区| 国产激情欧美一区二区| 午夜免费鲁丝| av天堂久久9| 日本vs欧美在线观看视频| 国产精品成人在线| 如日韩欧美国产精品一区二区三区| 亚洲欧美一区二区三区黑人| 成人免费观看视频高清| 一进一出好大好爽视频| 亚洲三区欧美一区| 精品福利观看| 在线视频色国产色| 国产精品亚洲av一区麻豆| 岛国在线观看网站| www.精华液| 亚洲人成电影免费在线| 久久99一区二区三区| 人人妻,人人澡人人爽秒播| 啦啦啦视频在线资源免费观看| 精品久久久久久电影网| 国产成人av教育| 亚洲免费av在线视频| 丝袜人妻中文字幕| 大型av网站在线播放| 国产亚洲欧美在线一区二区| 91成年电影在线观看| 午夜福利,免费看| 少妇被粗大的猛进出69影院| 日韩精品免费视频一区二区三区| 巨乳人妻的诱惑在线观看| 久久国产精品男人的天堂亚洲| 丁香欧美五月| 99久久人妻综合| 欧美日韩亚洲高清精品| 国产精品1区2区在线观看. | 亚洲精品成人av观看孕妇| 国产精品免费视频内射| www日本在线高清视频| 欧美精品亚洲一区二区| 变态另类成人亚洲欧美熟女 | 久久久久久久久久久久大奶| 捣出白浆h1v1| 国产蜜桃级精品一区二区三区 | 久久久国产成人精品二区 | 国产精品 欧美亚洲| 国产高清videossex| 成人影院久久| 久久国产精品人妻蜜桃| 99精品在免费线老司机午夜| 看免费av毛片| 国精品久久久久久国模美| 日日爽夜夜爽网站| 国产成+人综合+亚洲专区| 久久久精品区二区三区| 欧美 日韩 精品 国产| 91字幕亚洲| 国产又色又爽无遮挡免费看| 久久中文看片网| 宅男免费午夜| 五月开心婷婷网| 国产亚洲精品久久久久5区| 91av网站免费观看| 18在线观看网站| 亚洲av欧美aⅴ国产| videosex国产| 亚洲男人天堂网一区| 免费在线观看完整版高清| 亚洲精品国产精品久久久不卡| 亚洲av成人av| 国产一区二区激情短视频| 人成视频在线观看免费观看| 一边摸一边抽搐一进一小说 | 国产视频一区二区在线看| 飞空精品影院首页| 国产又爽黄色视频| 亚洲国产精品一区二区三区在线| 一二三四社区在线视频社区8| 国产精品九九99| 精品免费久久久久久久清纯 | 国产精品1区2区在线观看. | 五月开心婷婷网| 国产精品久久视频播放| 大香蕉久久成人网| tocl精华| 亚洲免费av在线视频| 久久青草综合色| 国产免费男女视频| 99国产精品一区二区三区| 成人特级黄色片久久久久久久| 国产在线一区二区三区精| 在线永久观看黄色视频| 啦啦啦 在线观看视频| 69av精品久久久久久| 免费在线观看黄色视频的| 亚洲成国产人片在线观看| 怎么达到女性高潮| 无限看片的www在线观看| 美女视频免费永久观看网站| 国产在线一区二区三区精| 久热这里只有精品99| 一本一本久久a久久精品综合妖精| 99国产精品一区二区蜜桃av | 日韩免费av在线播放| 欧美中文综合在线视频| 久久 成人 亚洲| 九色亚洲精品在线播放| 午夜免费成人在线视频| 日韩欧美三级三区| 欧美在线一区亚洲| 色婷婷av一区二区三区视频| xxxhd国产人妻xxx| 激情在线观看视频在线高清 | 欧美精品一区二区免费开放| 成人国语在线视频| 色尼玛亚洲综合影院| 丁香六月欧美| 欧美 亚洲 国产 日韩一| 麻豆成人av在线观看| 丁香欧美五月| 一个人免费在线观看的高清视频| 男人的好看免费观看在线视频 | 欧美精品高潮呻吟av久久| 人人妻人人爽人人添夜夜欢视频| 国产伦人伦偷精品视频| 日本一区二区免费在线视频| 精品第一国产精品| 水蜜桃什么品种好| 午夜免费成人在线视频| 岛国毛片在线播放| 18禁裸乳无遮挡免费网站照片 | 精品福利永久在线观看| 18禁黄网站禁片午夜丰满| 看免费av毛片| 国产精品亚洲一级av第二区| www.自偷自拍.com| 黄色女人牲交| 性色av乱码一区二区三区2| 后天国语完整版免费观看| 国产激情久久老熟女| 一本大道久久a久久精品| 老熟妇仑乱视频hdxx| 国产一卡二卡三卡精品| 精品免费久久久久久久清纯 | 超碰97精品在线观看| 成年动漫av网址| 大香蕉久久网| 欧美日韩成人在线一区二区| 人妻久久中文字幕网| 王馨瑶露胸无遮挡在线观看| 丰满饥渴人妻一区二区三| 日韩欧美一区视频在线观看| 国产成人精品在线电影| 香蕉国产在线看| 大码成人一级视频| netflix在线观看网站| 老熟妇仑乱视频hdxx| 一区二区日韩欧美中文字幕| 欧美日韩亚洲综合一区二区三区_| 麻豆国产av国片精品| 成人18禁在线播放| 又黄又爽又免费观看的视频| 黑人巨大精品欧美一区二区蜜桃| 中文字幕人妻丝袜一区二区| 国产精品偷伦视频观看了| 国产人伦9x9x在线观看| 一边摸一边抽搐一进一出视频| 日韩成人在线观看一区二区三区| 久久久久久久午夜电影 | 在线观看午夜福利视频| 在线av久久热| 天天躁日日躁夜夜躁夜夜| 亚洲精品在线美女| 欧美黑人精品巨大| cao死你这个sao货| 自拍欧美九色日韩亚洲蝌蚪91| 国产深夜福利视频在线观看| 韩国av一区二区三区四区| 国产精品98久久久久久宅男小说| 国产精品久久久av美女十八| 丝袜美腿诱惑在线| 亚洲人成伊人成综合网2020| 50天的宝宝边吃奶边哭怎么回事| 精品福利观看| 久久精品国产亚洲av香蕉五月 | av电影中文网址| 成人国产一区最新在线观看| 啦啦啦 在线观看视频| tocl精华| 人人妻人人添人人爽欧美一区卜| 亚洲五月色婷婷综合| 国产精品久久久人人做人人爽| 国产成人精品久久二区二区91| 亚洲成a人片在线一区二区| 色播在线永久视频| 久久久久久久国产电影| 国产欧美日韩精品亚洲av| tocl精华| 脱女人内裤的视频| 久久婷婷成人综合色麻豆| 男女床上黄色一级片免费看| 嫁个100分男人电影在线观看| 法律面前人人平等表现在哪些方面| 成年版毛片免费区| 精品一品国产午夜福利视频| 国产黄色免费在线视频| 高清在线国产一区| av免费在线观看网站| 精品国产美女av久久久久小说| 老鸭窝网址在线观看| 多毛熟女@视频| 一级黄色大片毛片| 久久精品国产亚洲av香蕉五月 | 欧美人与性动交α欧美软件| 免费看a级黄色片| 极品教师在线免费播放| e午夜精品久久久久久久| 午夜91福利影院| 波多野结衣av一区二区av| 亚洲欧美一区二区三区黑人| av网站在线播放免费| 国产精品国产高清国产av | 天堂动漫精品| 男人的好看免费观看在线视频 | 久久中文字幕人妻熟女| 最新在线观看一区二区三区| 亚洲自偷自拍图片 自拍| 久久国产乱子伦精品免费另类| 99精国产麻豆久久婷婷| 很黄的视频免费| 91九色精品人成在线观看| 99久久人妻综合| 精品国产一区二区三区久久久樱花| 久久久久视频综合| 国产精品国产高清国产av | 18禁裸乳无遮挡动漫免费视频| 亚洲va日本ⅴa欧美va伊人久久| 国产真人三级小视频在线观看| 亚洲七黄色美女视频| 国产精品免费大片| 又黄又爽又免费观看的视频| 午夜福利乱码中文字幕| 十八禁人妻一区二区| 国产99久久九九免费精品| 亚洲精品国产区一区二| 久久久久久久午夜电影 | 免费在线观看亚洲国产| av国产精品久久久久影院| 好男人电影高清在线观看| 大陆偷拍与自拍| 可以免费在线观看a视频的电影网站| 99国产精品一区二区三区| 一a级毛片在线观看| 久久久久久久久久久久大奶| 高清在线国产一区| 国产亚洲精品久久久久久毛片 | 色尼玛亚洲综合影院| 亚洲精品中文字幕在线视频| 十八禁网站免费在线| 如日韩欧美国产精品一区二区三区| 亚洲成人手机| 在线视频色国产色| 老鸭窝网址在线观看| 91老司机精品| 国产成人欧美| 岛国毛片在线播放| 老司机午夜福利在线观看视频| 亚洲精品一卡2卡三卡4卡5卡| 日本精品一区二区三区蜜桃| 日本一区二区免费在线视频| av不卡在线播放| 亚洲片人在线观看| 国产精品偷伦视频观看了| 少妇猛男粗大的猛烈进出视频| 国产在线一区二区三区精| 天天操日日干夜夜撸| 伦理电影免费视频| 美女视频免费永久观看网站| 亚洲人成77777在线视频| 欧美黄色淫秽网站| 757午夜福利合集在线观看| 日韩精品免费视频一区二区三区| 中国美女看黄片| 高潮久久久久久久久久久不卡| 日韩免费av在线播放| 亚洲av日韩在线播放| 国产aⅴ精品一区二区三区波| 飞空精品影院首页| 1024香蕉在线观看| 如日韩欧美国产精品一区二区三区| 久久青草综合色| 免费看a级黄色片| 这个男人来自地球电影免费观看| 午夜免费鲁丝| 亚洲专区中文字幕在线| 天堂中文最新版在线下载| 成人18禁在线播放| 老熟妇仑乱视频hdxx| 又黄又爽又免费观看的视频| 成年动漫av网址| 欧美日韩成人在线一区二区| 国精品久久久久久国模美| 精品亚洲成a人片在线观看| 大型av网站在线播放| 黄色怎么调成土黄色| 午夜影院日韩av| 欧美亚洲 丝袜 人妻 在线| 亚洲av成人一区二区三| 国产精品久久久人人做人人爽| 欧美日韩一级在线毛片| 国产免费男女视频| 午夜免费观看网址| 国产区一区二久久| 国产亚洲一区二区精品| 日韩视频一区二区在线观看| 大片电影免费在线观看免费| 妹子高潮喷水视频| 叶爱在线成人免费视频播放| 国产人伦9x9x在线观看| 久久久久久人人人人人| 少妇裸体淫交视频免费看高清 | 国产精品综合久久久久久久免费 | 久久精品国产99精品国产亚洲性色 | 又大又爽又粗| 一级毛片女人18水好多| 天天影视国产精品| 亚洲五月婷婷丁香| 亚洲专区字幕在线| 国产精品香港三级国产av潘金莲| 欧美亚洲 丝袜 人妻 在线| 久久精品国产综合久久久| 男女免费视频国产| tocl精华| 丝袜美足系列| 久久亚洲精品不卡| 99香蕉大伊视频| 国产精品乱码一区二三区的特点 | 黄片播放在线免费| 日本五十路高清| 一级,二级,三级黄色视频| 在线国产一区二区在线| 日韩欧美一区二区三区在线观看 | 中文字幕av电影在线播放| 精品视频人人做人人爽| 热99re8久久精品国产| av免费在线观看网站| 久久久国产精品麻豆| av天堂久久9| 午夜免费成人在线视频| 国产成人影院久久av| 天天躁日日躁夜夜躁夜夜| 日韩欧美三级三区| 亚洲欧美激情综合另类| 国产精品免费视频内射| 国产一区二区三区在线臀色熟女 | av一本久久久久| 中文字幕制服av| 中国美女看黄片| 亚洲av成人一区二区三| 日韩欧美一区视频在线观看| 欧美人与性动交α欧美精品济南到| 亚洲va日本ⅴa欧美va伊人久久| 日本黄色日本黄色录像| 色综合欧美亚洲国产小说| 12—13女人毛片做爰片一| 亚洲精品中文字幕一二三四区| 国产激情欧美一区二区| 久久狼人影院| 人人妻人人澡人人看| 麻豆成人av在线观看| 91成年电影在线观看| xxx96com| 搡老岳熟女国产| 18禁黄网站禁片午夜丰满| 国产无遮挡羞羞视频在线观看| 欧美日韩中文字幕国产精品一区二区三区 | 国产区一区二久久| 国产精品98久久久久久宅男小说| 黑丝袜美女国产一区| 真人做人爱边吃奶动态| 大码成人一级视频| 自线自在国产av| 老鸭窝网址在线观看| 在线观看舔阴道视频| 久久久水蜜桃国产精品网| 正在播放国产对白刺激| 国产在线观看jvid| 美女 人体艺术 gogo| 一区二区三区激情视频| 最近最新中文字幕大全免费视频| 在线国产一区二区在线| 亚洲,欧美精品.| www.精华液| 欧美av亚洲av综合av国产av| 国产精品美女特级片免费视频播放器 | 久久精品国产综合久久久| 两个人免费观看高清视频| 久久久久视频综合| 成人国产一区最新在线观看| 欧美亚洲 丝袜 人妻 在线| 9热在线视频观看99| 午夜亚洲福利在线播放| 成人精品一区二区免费| 深夜精品福利| 国产片内射在线| 精品一品国产午夜福利视频| 大码成人一级视频| 国产精品影院久久| 成人黄色视频免费在线看| 黑人巨大精品欧美一区二区mp4| a级毛片在线看网站| 老司机亚洲免费影院| 免费观看人在逋| 欧美乱色亚洲激情| 久久久久国内视频| 日韩欧美国产一区二区入口| 欧美中文综合在线视频| www.999成人在线观看| 天天躁日日躁夜夜躁夜夜| 最近最新免费中文字幕在线| 激情视频va一区二区三区| 一进一出抽搐动态| 久久精品熟女亚洲av麻豆精品| 人人妻人人添人人爽欧美一区卜| 精品国产一区二区三区久久久樱花| 97人妻天天添夜夜摸| 如日韩欧美国产精品一区二区三区| 极品少妇高潮喷水抽搐| 欧美国产精品一级二级三级| 久久久久久久午夜电影 | 国产主播在线观看一区二区| 50天的宝宝边吃奶边哭怎么回事| 中文字幕人妻丝袜制服| 麻豆乱淫一区二区| 亚洲人成电影免费在线| 国产不卡一卡二| 精品少妇久久久久久888优播| 中文字幕另类日韩欧美亚洲嫩草| 女同久久另类99精品国产91| aaaaa片日本免费| 在线永久观看黄色视频| 久久 成人 亚洲| 69av精品久久久久久| 免费在线观看影片大全网站| 自拍欧美九色日韩亚洲蝌蚪91| 亚洲第一青青草原| 精品久久久精品久久久| 一区二区三区精品91| 涩涩av久久男人的天堂| 正在播放国产对白刺激| 大陆偷拍与自拍| 国产精品免费视频内射| 精品一区二区三区四区五区乱码| 另类亚洲欧美激情| netflix在线观看网站| 夜夜爽天天搞| 嫩草影视91久久| 欧美日韩亚洲国产一区二区在线观看 | a在线观看视频网站| 亚洲精品在线观看二区| 国产精品一区二区在线不卡| 极品教师在线免费播放| 国产亚洲一区二区精品| 久久精品国产99精品国产亚洲性色 | 丁香六月欧美| 一级片免费观看大全| 自线自在国产av| 国产精品永久免费网站| 国产精品自产拍在线观看55亚洲 | 视频区图区小说| 久热这里只有精品99| 热re99久久国产66热| 午夜福利在线观看吧| 亚洲专区中文字幕在线| 成人18禁在线播放| 国产男女超爽视频在线观看| 亚洲中文字幕日韩| 最新在线观看一区二区三区| 女人高潮潮喷娇喘18禁视频| xxx96com| 亚洲国产精品一区二区三区在线| 久久久久久人人人人人| svipshipincom国产片| 女人久久www免费人成看片| 国产精品一区二区精品视频观看| 亚洲人成伊人成综合网2020| 亚洲成人国产一区在线观看| 高清欧美精品videossex| 亚洲五月婷婷丁香| 国产在线观看jvid| 亚洲色图综合在线观看| av超薄肉色丝袜交足视频| 乱人伦中国视频| 久久人妻av系列| 叶爱在线成人免费视频播放| 不卡一级毛片| 美女扒开内裤让男人捅视频| 露出奶头的视频| 亚洲中文字幕日韩| 欧美一级毛片孕妇| 女性生殖器流出的白浆| 一边摸一边做爽爽视频免费| 夜夜夜夜夜久久久久| 天堂动漫精品| 精品福利观看| 99热网站在线观看| 亚洲熟妇中文字幕五十中出 | 国产高清激情床上av| 亚洲av成人一区二区三| 国产精华一区二区三区| 精品国产一区二区三区四区第35| 成人精品一区二区免费| 嫁个100分男人电影在线观看| 丁香欧美五月| 99国产精品99久久久久| 男人的好看免费观看在线视频 | 一级片免费观看大全| 国产男靠女视频免费网站| 亚洲欧美日韩高清在线视频| 国产成人精品久久二区二区91| 亚洲全国av大片| 国产aⅴ精品一区二区三区波| 中文字幕另类日韩欧美亚洲嫩草| 欧美 日韩 精品 国产| 国产激情久久老熟女| 成在线人永久免费视频| 国产亚洲一区二区精品| 欧美日韩亚洲国产一区二区在线观看 | 国产国语露脸激情在线看| 国产一卡二卡三卡精品| 国产亚洲精品一区二区www | 动漫黄色视频在线观看| 黄色片一级片一级黄色片| 美女午夜性视频免费| 精品一区二区三区av网在线观看| 老司机影院毛片| 97人妻天天添夜夜摸| 69精品国产乱码久久久| 他把我摸到了高潮在线观看| 欧美 亚洲 国产 日韩一| 午夜精品久久久久久毛片777| 久久精品国产亚洲av香蕉五月 | 大香蕉久久成人网| 国产不卡一卡二|