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

    鐵顆粒非均相著火模型與著火特性研究

    2024-02-05 02:07:22徐義華陳宣任楊志豪孫???/span>任冠龍
    燃燒科學(xué)與技術(shù) 2024年1期
    關(guān)鍵詞:鐵粉氧氣粒徑

    蔡 平,徐義華,胡 坤,陳宣任,楊志豪,孫???,任冠龍

    鐵顆粒非均相著火模型與著火特性研究

    蔡 平,徐義華,胡 坤,陳宣任,楊志豪,孫???,任冠龍

    (南昌航空大學(xué)飛行器工程學(xué)院江西省微小航空發(fā)動(dòng)機(jī)重點(diǎn)實(shí)驗(yàn)室,南昌 330063)

    基于縮小未反應(yīng)核模型,考慮氧化層對氧氣的擴(kuò)散影響,構(gòu)建了空氣中Fe顆粒的非均相著火模型.模型計(jì)算的著火溫度與實(shí)驗(yàn)著火溫度的最大誤差為6.95%,與Mi等的著火模型對比,兩者所計(jì)算的顆粒升溫速率基本一致,著火延遲時(shí)間的數(shù)量級相同.基于該模型研究了環(huán)境壓強(qiáng)、氧氣濃度、初始氧化層厚度和顆粒粒徑對著火溫度、著火延遲的影響.結(jié)果表明:著火溫度、著火延遲時(shí)間隨環(huán)境壓強(qiáng)、氧氣濃度增大而減小,隨初始氧化層厚度增厚而增大;粒徑小于10μm的范圍內(nèi),著火溫度隨粒徑減小而增大,粒徑在10~50μm范圍內(nèi),粒徑對著火溫度影響不大;隨著粒徑增加,著火延遲時(shí)間增加.

    鐵顆粒;燃料;非均相著火模型;著火溫度;著火延遲

    金屬粉末具有較高的能量密度和更高的比能量,作為推進(jìn)劑輔助添加劑被廣泛應(yīng)用于軍事和商業(yè)領(lǐng)域[1-7].根據(jù)金屬燃料的燃燒特性可知,其產(chǎn)物為固態(tài)金屬氧化物,屬于零碳排放能源.金屬粉末燃燒生成的固態(tài)金屬氧化物利用太陽能等清潔能源電解冶煉可循環(huán)利用.另外,金屬燃料比太陽能、風(fēng)能、地?zé)岬瘸跫壞茉锤菀捉鉀Q能源的儲(chǔ)存、運(yùn)輸和交易問題,使初級能源的生產(chǎn)和最終用途在空間和時(shí)間上分離,從而實(shí)現(xiàn)世界范圍內(nèi)的廣泛使用[8-9].

    在眾多可以利用的金屬燃料中,金屬鐵顆粒在燃燒中氧化產(chǎn)物直接附著在金屬鐵核表面,使氧化物顆粒尺寸最終大于初始鐵顆粒.與鎂、鋁等金屬的液滴蒸發(fā)氣相燃燒模式相比,鐵顆粒的燃燒模式,更有利于燃燒產(chǎn)物的回收與利用;加上金屬鐵在地表含量高、分布廣、易開采、成本低的特點(diǎn),使金屬鐵粉在未來有望成為一種極具潛力的零碳、清潔、可回收的綠色能源[6,8-9].Mandilas等[10]利用鐵粉燃燒特性,搭建了金屬燃料內(nèi)燃機(jī),進(jìn)而研究了鐵粉在內(nèi)燃機(jī)內(nèi)的燃燒特性.劉世寧等[11]則設(shè)計(jì)出鐵粉燃燒器并實(shí)現(xiàn)常溫、常壓下鐵粉與空氣的自持燃燒.然而,目前所設(shè)計(jì)的燃燒器因其燃燒效率低、燃燒室堵塞等問題阻礙了鐵粉燃料的使用,因此需要深入理解鐵顆粒復(fù)雜的著火機(jī)理,建立其著火模型,為鐵粉燃燒器的科學(xué)設(shè)計(jì)、鐵能源的高效利用和燃燒產(chǎn)物的有效回收提供一些理論依據(jù).

    當(dāng)前,鐵的物理化學(xué)特性已經(jīng)被熟知.由于鐵的熔點(diǎn)(1809K)和沸點(diǎn)(3135K)較高,常溫下的鐵為固態(tài).空氣中鐵被氧化后主要生成+2價(jià)和+3價(jià)的鐵氧化物,其中FeO和Fe3O4的熔點(diǎn)分別為1650K和1870K,而Fe2O3在1733K左右時(shí)會(huì)被分解為Fe3O4和O2[12].其余氧化物也分別在3537K和2257K時(shí)被分解Fe和O2[13].已有文獻(xiàn)對鐵的著火進(jìn)行了相關(guān)的研究,高文靜等[14]采用普適積分法測量出鐵粉燃燒時(shí)的活化能和指前因子.孫金華等[15]通過實(shí)驗(yàn)發(fā)現(xiàn)了鐵粉燃燒屬于非均相燃燒,鐵顆粒燃燒狀態(tài)為液態(tài).Li等[16]通過實(shí)驗(yàn)證明了這一點(diǎn),他們測量了鐵顆粒燃燒中的顆粒溫度隨時(shí)間的變化,發(fā)現(xiàn)顆粒的燃燒溫度高于鐵及其氧化物的熔點(diǎn).Mandilas等[10]搭建了鐵粉燃燒器并測量了其燃燒的溫度,得到最高溫度在2200~2300K,該溫度高于鐵熔點(diǎn)及其氧化物的熔點(diǎn),低于鐵沸點(diǎn).

    Leshchevich等[17]通過實(shí)驗(yàn)確定了鐵粉著火燃燒時(shí)的活化能,并根據(jù)Arrhenius方程建立了鐵粉的半經(jīng)驗(yàn)著火模型.丁宋毅等[18]則是基于粒子群算法擬合出鐵粉的燃燒速率方程,進(jìn)而建立其著火燃燒模型. Schiemann等[19]認(rèn)為著火反應(yīng)區(qū)位于鐵顆粒表面,反應(yīng)產(chǎn)物(Fe3O4)吸附在顆粒表面,構(gòu)建了鐵顆粒表面反應(yīng)的著火模型,該模型未考慮氧化層對著火的影響.在此基礎(chǔ)上,Mi等[20]引入了著火反應(yīng)速率與氧化層厚度呈反比關(guān)系,并基于能量守恒定律,建立了鐵顆粒著火模型.

    由于空氣中鐵顆粒的著火屬于氣固反應(yīng),其著火模型的建立可以參考較為成熟的冶金宏觀動(dòng)力模型,即縮小未反應(yīng)核模型[21].該模型指出:由于固體反應(yīng)物是無孔的,初始反應(yīng)僅發(fā)生在顆粒表面并形成很薄且疏松的固體產(chǎn)物層,氣體反應(yīng)物通過該產(chǎn)物層擴(kuò)散并與固體反應(yīng)物核反應(yīng).隨著反應(yīng)進(jìn)行,產(chǎn)物層逐漸增厚,內(nèi)部的未反應(yīng)核區(qū)域逐漸縮小,直至消失.

    綜上所述,目前有關(guān)鐵顆粒著火模型的研究主要分兩種,第1種是由實(shí)驗(yàn)結(jié)果構(gòu)成的經(jīng)驗(yàn)公式,第2種是未考慮氧化層對氧氣擴(kuò)散影響的顆粒表面反應(yīng)模型,或引入反應(yīng)速率與氧化層呈反比的表面反應(yīng)模型.然而,在鐵顆粒的著火過程中,氧氣在氧化層內(nèi)的傳輸過程對顆粒著火產(chǎn)生一定的影響.因此,本文將基于縮小未反應(yīng)核模型和謝苗諾夫著火理論,考慮氧化層內(nèi)的氧氣傳輸以及顆粒所處環(huán)境區(qū)域內(nèi)的溫度和氧氣分布,構(gòu)建空氣中Fe顆粒非均相著火模型.在此基礎(chǔ)上,應(yīng)用該模型討論環(huán)境壓強(qiáng)、氧氣濃度、初始氧化層厚度和顆粒粒徑對著火溫度、著火延遲的影響.

    1 Fe顆粒著火模型

    1.1 Fe顆粒著火過程中氧化層結(jié)構(gòu)

    空氣中鐵顆粒的著火過程,本質(zhì)是鐵在空氣中的高溫氧化過程.由文獻(xiàn)[22-24]可知,當(dāng)溫度低于570℃時(shí),F(xiàn)e顆粒的氧化層由Fe3O4和Fe2O3組成的兩層結(jié)構(gòu).如圖1(a)所示,F(xiàn)e3O4層緊挨著單質(zhì)鐵基體,F(xiàn)e2O3層則是外氧化層.因?yàn)闇囟鹊陀?70℃時(shí),F(xiàn)eO不穩(wěn)定,其易分解成Fe和Fe3O4.當(dāng)溫度高于570℃時(shí),F(xiàn)eO可以穩(wěn)定存在,F(xiàn)e顆粒氧化膜由3層氧化物組成.如圖1(b)所示,F(xiàn)eO緊貼單質(zhì)鐵基體,中間層是Fe3O4,最外層是Fe2O3.但文獻(xiàn)[25]表明,在中等溫度(700~1250℃)的條件下,F(xiàn)e2O3厚度占整個(gè)氧化層厚度的1%,說明Fe2O3層對氧氣擴(kuò)散的影響幾乎為零,且生成Fe2O3的反應(yīng)速率幾乎為零,因此在鐵顆粒著火過程中可以忽略外層Fe2O3,將氧化層視為內(nèi)層為FeO和外層為Fe3O4的兩層結(jié)構(gòu),如圖2所示.由于著火過程中生成的氧化層厚度與粒徑的比值很小,忽略Fe2O3層對計(jì)算結(jié)果幾乎不產(chǎn)生影響.此外,這兩層結(jié)構(gòu)都是多孔結(jié)構(gòu),其形成機(jī)理可參見文獻(xiàn)[26-27].基于以上鐵的高溫氧化機(jī)理以及氧化層的結(jié)構(gòu)構(gòu)建鐵顆粒著火模型.

    (a)溫度低于570℃

    (b)溫度高于570℃

    圖1 顆粒的氧化層結(jié)構(gòu)

    Fig.1 Oxide layer structure of particle

    圖2 不考慮Fe2O3層時(shí)顆粒的氧化層結(jié)構(gòu)

    1.2 著火機(jī)理

    由于鐵顆粒著火過程較為復(fù)雜,影響著火因素較多,為了簡化著火模型,對整個(gè)著火過程做出如下假設(shè):

    (1)顆粒為球形且內(nèi)部溫度均勻.因?yàn)轭w粒尺寸極小且內(nèi)部傳熱速度很大.

    (2)著火過程中,不考慮多孔氧化層內(nèi)部的化學(xué)反應(yīng),反應(yīng)區(qū)位于鐵核與氧化層界面上.

    (3)雙層氧化層結(jié)構(gòu)可以近似簡化為單層結(jié)構(gòu).因?yàn)樵贔eO層內(nèi),F(xiàn)e與O的摩爾比接近1,且中等溫度下摩爾比隨溫度變化而在0.84~0.95內(nèi)變動(dòng).此外,在FeO層存在較多Fe3O4聚集區(qū)[28],可以認(rèn)為FeO層主要由Fe和Fe3O4構(gòu)成.需要說明的是,這種近似簡化會(huì)導(dǎo)致計(jì)算的反應(yīng)放熱速率比實(shí)際略大.因?yàn)樯蒄eO的反應(yīng)熱小于生成Fe3O4的反應(yīng)熱,將氧化層簡化后,使原本一部分氧化成FeO的Fe直接氧化成Fe3O4.

    (4)在顆粒升溫過程中,不考慮顆粒受熱膨脹,也不考慮a-Fe到g-Fe再到δ-Fe的相變過程,因?yàn)橄嘧儩摕徇h(yuǎn)小于反應(yīng)放熱,用于計(jì)算鐵的熱力學(xué)性質(zhì)采用α-Fe的熱力學(xué)性質(zhì).

    (5)在整個(gè)著火過程中,顆粒僅與氧氣反應(yīng),其他氣體稱為惰性氣體.顆粒處于相對靜止、干燥的環(huán)境中.

    (6)為1,即ρcD

    基于以上假設(shè),將著火區(qū)域分為內(nèi)部未反應(yīng)鐵核(Ⅰ區(qū))、多孔Fe3O4氧化層(Ⅱ區(qū))和空間環(huán)境(Ⅲ區(qū)),如圖3所示.空氣中的O2先從環(huán)境傳輸至顆粒表面,再通過多孔Fe3O4氧化層擴(kuò)散到鐵-氧化層界面上,并與Fe反應(yīng)生成Fe3O4,如式(1)所示.反應(yīng)產(chǎn)物堆積在氧化層內(nèi)表面,使氧化層朝顆粒內(nèi)部方向逐漸增厚,未反應(yīng)鐵核區(qū)不斷縮?。?/p>

    此外,有兩種熱源使顆粒升溫,一是外部環(huán)境熱量通過導(dǎo)熱和熱輻射方式傳輸至顆粒,二是反應(yīng)熱.

    圖3 鐵顆粒著火模型原理示意

    1.3 著火控制方程

    依據(jù)上述著火機(jī)理,著火控制方程可簡化為徑向一維球坐標(biāo)描述,以鐵顆粒的中心為坐標(biāo)原點(diǎn),建立球坐標(biāo)系,如圖4所示.由圖可知,氧氣通過濃度梯度在區(qū)域Ⅲ和Ⅱ中擴(kuò)散.

    1.3.1 空間區(qū)域Ⅲ的控制方程

    圖4 球坐標(biāo)系下的氧氣濃度分布

    顆粒處于穩(wěn)定狀態(tài),依據(jù)控制體內(nèi)能量增量等于外界傳入的熱量,可獲得空間區(qū)域的能量方程,如式(5)所示.

    1.3.2 氧化層區(qū)域Ⅱ的控制方程

    由于氧化層為多孔結(jié)構(gòu),考慮多孔結(jié)構(gòu)對氧氣擴(kuò)散影響,引入氧化層的孔隙度p和孔隙的曲折度p兩種參量[29].其中孔隙度為孔隙空間占整個(gè)氧化層的體積比值,一般是小于1的正數(shù);而曲折度是大于1的數(shù),因?yàn)檠趸瘜觾?nèi)的孔徑遠(yuǎn)大于反應(yīng)物分子的平均自由程,同時(shí)孔的走向是蜿蜒曲折的,迫使反應(yīng)物在孔內(nèi)傳輸路徑要比直線傳輸長.此外,文獻(xiàn)[30]表明,孔隙度與曲折度呈反比關(guān)系且乘積為1.

    依據(jù)著火機(jī)理,顆粒內(nèi)部的化學(xué)反應(yīng)是不可逆的,且反應(yīng)物Fe為固態(tài),反應(yīng)速率由O2的濃度及其化學(xué)動(dòng)力學(xué)機(jī)理決定.依據(jù)化學(xué)速率由反應(yīng)物濃度控制的質(zhì)量作用定律,采用O2的消耗速率來表示反應(yīng)速率,如式(9)所示.化學(xué)反應(yīng)速率常數(shù)(re1)可用Arrhenius函數(shù)描述[31],即

    在穩(wěn)態(tài)條件下,氧氣擴(kuò)散速度等于反應(yīng)速率,即

    1.3.3 氧化層內(nèi)徑及顆粒溫度的控制方程

    聯(lián)立方程(11)和(12)得到氧化層內(nèi)徑與時(shí)間的變化關(guān)系,如式(13)所示.

    依據(jù)能量守恒關(guān)系,得到顆粒溫度s與時(shí)間關(guān)系式,即

    式(14)中等號右邊第1項(xiàng)為對流換熱,第2項(xiàng)為輻射傳熱,第3項(xiàng)為反應(yīng)熱.

    2 計(jì)算結(jié)果與分析

    2.1 數(shù)值計(jì)算方法

    圖5 數(shù)值計(jì)算流程

    2.2 模型校核

    Arshad等[32]利用Godbert-Greenwald (G-G)加熱爐構(gòu)建了鐵粉的著火溫度測量系統(tǒng).該測量系統(tǒng)用于觀察鐵粉在某一環(huán)境溫度下是否被點(diǎn)燃.如果鐵粉可以在一定的最低環(huán)境溫度下點(diǎn)燃,認(rèn)為該環(huán)境溫度為著火溫度且這種測量法是依據(jù)謝苗諾夫著火理論[31].他們利用該系統(tǒng)測量出不同環(huán)境壓強(qiáng)下粒徑為32~50μm、質(zhì)量濃度為370g/m3的鐵粉的著火溫度,結(jié)果如表1所示.

    表1 不同壓力下鐵粉在空氣中的實(shí)驗(yàn)著火溫度[32]

    Tab.1 Experimental ignition temperature of iron powder in air with different pressures[32]

    對于表1中的所有實(shí)驗(yàn)工況,使用本文所建立的著火模型來計(jì)算和分析它們的著火過程.其中模型計(jì)算采用的顆粒直徑為實(shí)驗(yàn)粒徑的平均值,即=41μm.通過模型計(jì)算可以得到著火過程中顆粒溫度與時(shí)間的關(guān)系.以該表的工況1為例,繪制了不同環(huán)境溫度下,顆粒溫度隨時(shí)間的變化關(guān)系,如圖6所示.當(dāng)環(huán)境溫度T≤867K時(shí),顆粒的最高溫度停留在環(huán)境溫度附近.一旦將環(huán)境溫度提升至868K,那么顆粒溫度在0.34s左右時(shí)突破環(huán)境溫度,且急劇升高.說明此時(shí)的鐵顆粒已經(jīng)著火成功,因此把環(huán)境溫度868K作為鐵顆粒的著火溫度.

    圖6 不同環(huán)境溫度下顆粒溫度隨時(shí)間變化

    其他工況所計(jì)算的著火溫度如表2所示.由該表可知,模型計(jì)算的結(jié)果與實(shí)驗(yàn)結(jié)果誤差均小于7%.另外,計(jì)算著火溫度的變化趨勢與實(shí)驗(yàn)結(jié)果相反,因?yàn)锳rshad等往加熱爐內(nèi)通入了常溫空氣.點(diǎn)燃鐵粉的同時(shí),爐內(nèi)的空氣也會(huì)被加熱.隨著爐內(nèi)壓強(qiáng)的增大,通入的空氣增多,要使鐵粉在3s內(nèi)被點(diǎn)燃,需要提高爐內(nèi)溫度.

    表2 計(jì)算著火溫度與實(shí)驗(yàn)著火溫度對比

    Tab.2 Comparisonbetween calculated ignition tempera-ture and experimental ignition temperature

    為了進(jìn)一步說明著火模型的特點(diǎn),將本模型與Mi等[20]的著火模型進(jìn)行對比.在Mi模型中,Mi認(rèn)為Fe離子在氧化層中的傳輸性能比O離子強(qiáng),氧化層的生長主要受Fe離子向外擴(kuò)散控制,內(nèi)部鐵核縮小的同時(shí),氧化層向外增厚,其構(gòu)建的反應(yīng)速率與氧化層厚度呈反比關(guān)系.設(shè)鐵顆粒半徑s=250nm,初始氧化層厚度為3nm,初始顆粒溫度298.15K,在環(huán)境溫度T=1100K、壓強(qiáng)p=101325Pa、氧氣占比為21%的干燥空氣中的條件下,分別利用本文模型和Mi模型進(jìn)行數(shù)值計(jì)算,得到兩種模型下顆粒溫度隨時(shí)間變化,如圖7所示.由圖可知,在鐵顆粒著火之前,兩者所計(jì)算的顆粒溫度隨時(shí)間變化的趨勢一致,著火延遲時(shí)間分別是5.65μs和8.83μs,兩者著火延時(shí)時(shí)間的數(shù)量級相同.

    圖7 Mi模型和本文模型所計(jì)算的顆粒溫度隨時(shí)間變化

    2.3 著火過程參數(shù)分析

    設(shè)鐵顆粒半徑s=250nm,依據(jù)文獻(xiàn)[10,33]可取其初始氧化層厚度為3nm,則其初始氧化層內(nèi)徑b=247nm,初始顆粒溫度298.15K,在環(huán)境溫度T=1000K、壓強(qiáng)p=101325Pa、氧氣占比為21%的干燥空氣中的條件下,通過數(shù)值計(jì)算得到氧化層內(nèi)徑和顆粒溫度隨時(shí)間變化(見圖8),以及在不同顆粒溫度時(shí)所對應(yīng)的空間溫度分布(見圖9)和空間氧氣摩爾濃度分布(見圖10).其中對于空間距離采用=/s,對空間坐標(biāo)進(jìn)行無量綱化.

    由圖8可知,隨著反應(yīng)時(shí)間增加,顆粒升溫分為3個(gè)階段,即顆粒在對流換熱作用下的急劇升溫階段(時(shí)間小于5μs)、顆粒溫度接近環(huán)境溫度的緩慢升溫階段(時(shí)間在5μs至16μs之間),以及著火后的顆粒快速升溫階段(時(shí)間大于16μs).反應(yīng)產(chǎn)物堆積在氧化層內(nèi)表面,氧化層的厚度可由其內(nèi)徑表示,厚度增大,則內(nèi)徑減?。稍搱D可知,隨著反應(yīng)時(shí)間增加,氧化層內(nèi)徑先緩慢減小,在16μs時(shí),由于顆粒已經(jīng)被點(diǎn)燃,內(nèi)徑開始急劇減小,即氧化層厚度快速增大.

    圖8 氧化層內(nèi)徑和顆粒溫度隨時(shí)間變化

    由圖9和圖10可知,在時(shí)間從0ms增加至7.5μs過程中,由于顆粒溫度始終不超過環(huán)境溫度(1000K),鐵顆粒所處的空間區(qū)域中,隨著空間距離的增加,空間溫度的走向可分為兩部分,即空間溫度急劇增大(距離/s小于4左右),以及其緩慢增大(距離/s大于4).在此過程內(nèi),反應(yīng)速率主要受化學(xué)動(dòng)力學(xué)控制,隨著時(shí)間增加,顆粒溫度升高,反應(yīng)速率有微小增加,使顆粒周圍的氧氣摩爾濃度微量減?。谥疬^程的7.5μs至16.0μs中,由于顆粒溫度高于環(huán)境溫度,空間溫度隨空間距離的增加而逐漸減?。诉^程中,顆粒溫度接近著火溫度,氧氣在空間區(qū)域(Ⅲ區(qū))和氧化層(Ⅱ區(qū))中的擴(kuò)散效應(yīng)比前一過程明顯,隨著時(shí)間增加,顆粒的反應(yīng)速率受化學(xué)動(dòng)力學(xué)控制逐步過渡到受化學(xué)動(dòng)力學(xué)與氧氣擴(kuò)散機(jī)制共同控制,且反應(yīng)速率有較大增加,使顆粒周圍的氧氣摩爾濃度有了較大的減?。浑S著空間距離的增加,氧氣濃度先急劇增大(距離/s小于4左右),再緩慢增大(距離/s大于4).此外,當(dāng)空間距離足夠大時(shí),空間溫度和氧氣濃度都趨于環(huán)境值.

    (a)氧氣摩爾濃度分布

    (b)圖(a)中局部放大

    圖10 不同時(shí)刻下的氧氣摩爾濃度分布

    Fig.10 O2molar concentration distribution under differ-ent time

    為了更加直觀體現(xiàn)鐵顆粒著火過程中的空間參數(shù)變化規(guī)律,圖11給出了不同時(shí)刻的溫度分布和氧氣摩爾濃度分布云圖.圖11(b)中間白色部分為顆粒.由圖11可知,在著火過程中,隨著時(shí)間增加,在離開顆粒方向的空間溫度梯度由正梯度逐漸轉(zhuǎn)變?yōu)樨?fù)梯度,即隨著著火時(shí)間的推移,顆粒溫度由常溫開始迅速增大到接近環(huán)境溫度,再在反應(yīng)生成熱作用下顆粒溫度進(jìn)一步升高,直至著火成功;而氧氣摩爾濃度梯度在離開顆粒方向保持正梯度,且梯度隨時(shí)間的增加而增大.

    2.4 環(huán)境壓強(qiáng)、氧氣濃度、初始氧化層厚度和顆粒粒徑對著火溫度的影響

    2.4.1 環(huán)境壓強(qiáng)的影響

    設(shè)鐵顆粒半徑s=250nm,初始氧化層內(nèi)徑b=247nm,初始顆粒溫度298.15K,在氧氣占比為21%的干燥空氣中的條件下,計(jì)算了環(huán)境壓強(qiáng)為0.05~0.5MPa的著火過程,得到了著火溫度隨環(huán)境壓強(qiáng)的變化,如圖12所示.由圖可知,著火溫度隨環(huán)境壓強(qiáng)增大而降低,因?yàn)榄h(huán)境壓強(qiáng)與濃度呈正比,隨著環(huán)境壓強(qiáng)增大,氧氣濃度增加,增大了著火過程的顆粒反應(yīng)速率,使單位時(shí)間內(nèi)反應(yīng)熱增多,降低了著火溫度.

    (a)溫度

    (b)氧氣摩爾濃度

    圖11 著火過程中空間參數(shù)分布

    Fig.11 Space parameter distribution during ignition

    圖12 著火溫度隨環(huán)境壓強(qiáng)的變化

    Fig.12 Ignition temperature versus atmospheric pressure

    2.4.2 氧氣濃度的影響

    設(shè)鐵顆粒半徑s=250nm,初始氧化層內(nèi)徑b=247nm,初始顆粒溫度298.15K,在壓強(qiáng)p=101325Pa干燥氣氛中的條件下,由于壓強(qiáng)不變,要改為氧氣濃度則只需改變氧氣摩爾分?jǐn)?shù),設(shè)氣氛中氧氣的摩爾分?jǐn)?shù)從0.21均勻變化到0.66,計(jì)算得到的著火溫度隨氧氣摩爾分?jǐn)?shù)的變化如圖13所示.由圖可知,著火溫度隨氧氣濃度增大而減?。?/p>

    2.4.3 初始氧化層厚度的影響

    設(shè)鐵顆粒半徑s=250nm,初始顆粒溫度298.15K,在壓強(qiáng)p=101325Pa、氧氣占比為21%的干燥空氣中的條件下,計(jì)算了初始氧化層厚度為0~45nm的著火過程,得到了著火溫度隨初始氧化層厚度的變化,如圖14所示.由圖可知,著火溫度隨初始氧化層厚度增厚而增大,因?yàn)殡S著氧化層厚度增大,阻礙了氧氣的擴(kuò)散,鐵顆粒反應(yīng)速率降低,需要更高的環(huán)境溫度彌補(bǔ)顆粒升溫速率,使其達(dá)到成功著火溫度.

    圖13 著火溫度隨氧氣摩爾分?jǐn)?shù)的變化

    圖14 著火溫度隨初始氧化層厚度的變化

    2.4.4 顆粒粒徑的影響

    設(shè)鐵顆粒初始氧化層厚度3nm和初始顆粒溫度298.15K,在壓強(qiáng)p=101325Pa、氧氣占比為21%的干燥空氣中的條件下,計(jì)算了顆粒直徑為100nm至50μm的著火過程,得到了著火溫度隨顆粒粒徑的變化,如圖15所示.由圖可知,粒徑在100nm到10μm左右,著火溫度隨粒徑增大而顯著減小,因?yàn)榉磻?yīng)速率遵循線性規(guī)律,導(dǎo)致著火溫度與顆粒尺寸有關(guān)[34].對于等質(zhì)量的納米級顆粒群[35],其著火規(guī)律與單顆粒相反,即由于納米顆粒群效應(yīng)可大大降低其著火溫度,隨著粒徑減小,比表面積增大,顆粒的著火溫度降低.在粒徑為10μm至50μm之間,著火溫度隨粒徑增大而稍有降低,可以認(rèn)為該粒徑范圍內(nèi),粒徑對著火溫度的影響較?。?/p>

    圖15 著火溫度隨顆粒直徑的變化

    2.5 環(huán)境壓強(qiáng)、氧氣濃度、初始氧化層厚度和顆粒粒徑對著火延遲的影響

    基于控制變量法,將研究環(huán)境壓強(qiáng)、氧氣濃度、初始氧化層厚度和顆粒粒徑對著火延遲的影響.其中,著火延遲時(shí)間為顆粒從初始狀態(tài)到顆粒溫度開始急劇上升時(shí)的時(shí)間變化量,即著火時(shí)刻的時(shí)間量為著火延遲時(shí)間.

    2.5.1 環(huán)境壓強(qiáng)的影響

    設(shè)鐵顆粒半徑s=250nm,初始氧化層內(nèi)徑b=247nm,初始顆粒溫度298.15K,在溫度T=1100K、氧氣占比為21%的干燥空氣中的條件下,計(jì)算了環(huán)境壓強(qiáng)為0.05~0.5MPa的著火過程,得到了不同環(huán)境壓強(qiáng)下顆粒溫度隨時(shí)間的變化,如圖16所示.由圖可知,在顆粒著火成功之前,不同環(huán)境壓強(qiáng)下,顆粒溫度隨時(shí)間變化的趨勢基本不變,即環(huán)境壓強(qiáng)對顆粒升溫速率幾乎沒有影響,因?yàn)橹疬^程中顆粒溫度不高,反應(yīng)熱比較小,顆粒升溫所需的能量主要來自環(huán)境對顆粒的加熱.在環(huán)境壓強(qiáng)為0.05MPa、0.10MPa、0.15MPa、0.30MPa和0.50MPa的條件下,所計(jì)算的著火延遲時(shí)間分別為7.1μs、5.7μs、5.1μs、4.5μs和4μs,說明著火延遲時(shí)間隨環(huán)境壓強(qiáng)增大而縮短,這主要因?yàn)轭w粒溫度達(dá)到著火臨界點(diǎn)溫度后,其反應(yīng)速率受化學(xué)動(dòng)力學(xué)與擴(kuò)散機(jī)制共同控制,則環(huán)境壓強(qiáng)越大,氧氣濃度越高,使得其反應(yīng)速率增大,放出的反應(yīng)熱增多,從而減少顆粒達(dá)到著火點(diǎn)的時(shí)間.

    圖16 不同環(huán)境壓強(qiáng)下顆粒溫度隨時(shí)間變化

    2.5.2 氧氣濃度的影響

    設(shè)鐵顆粒半徑s=250nm,初始氧化層內(nèi)徑b=247nm,初始顆粒溫度298.15K,在壓強(qiáng)p=101325Pa、溫度T=1100K干燥氣氛中的條件下,計(jì)算了氧氣摩爾分?jǐn)?shù)為0.21~0.66的著火過程,得到了不同氧氣摩爾分?jǐn)?shù)下,顆粒溫度隨時(shí)間的變化,如圖17所示.由圖可知,在著火成功之前,不同氧氣濃度下,顆粒溫度隨時(shí)間變化的趨勢基本不變,即氧氣濃度對顆粒升溫速率幾乎沒有影響.在氧氣摩爾分?jǐn)?shù)為0.21、0.26、0.31、0.46和0.66的條件下,所計(jì)算的著火延遲時(shí)間分別為5.7μs、5.4μs、5.1μs、4.7μs和4.3μs,說明著火延遲時(shí)間隨氧氣濃度增大而縮短,這與環(huán)境壓強(qiáng)的影響相似,所導(dǎo)致的機(jī)理也相同.

    圖17 不同氧氣摩爾分?jǐn)?shù)下顆粒溫度隨時(shí)間變化

    2.5.3 初始氧化層厚度的影響

    設(shè)鐵顆粒半徑s=250nm,初始顆粒溫度298.15K,在壓強(qiáng)p=101325Pa、溫度T=1100K、氧氣占比為21%的干燥空氣中的條件下,計(jì)算了初始氧化層厚度為0~45nm的著火過程,得到了不同初始氧化層厚度下,顆粒溫度隨時(shí)間變化,如圖18所示.由圖可知,在著火成功之前,不同初始氧化層厚度下,顆粒溫度隨時(shí)間變化的趨勢基本不變,即初始氧化層厚度對顆粒升溫速率幾乎沒有影響.在初始氧化層厚度為0、20nm、35nm、40nm和45nm的條件下,所計(jì)算的著火延遲時(shí)間分別為5.6μs、5.9μs、6.1μs、6.2μs和6.3μs,說明著火延遲時(shí)間隨初始氧化層厚度增厚而增加,因?yàn)轭w粒溫度達(dá)到著火臨界點(diǎn)溫度后,反應(yīng)速率受化學(xué)動(dòng)力學(xué)與擴(kuò)散機(jī)制共同控制,隨氧化層厚度增厚,氧氣擴(kuò)散速率減少,反應(yīng)速率減小,增加了顆粒達(dá)到著火點(diǎn)的時(shí)間.

    2.5.4 顆粒粒徑的影響

    設(shè)初始氧化層厚度3nm、初始顆粒溫度298.15K,在壓強(qiáng)p=101325Pa、溫度T=1100K、氧氣占比為21%的干燥空氣中的條件下,計(jì)算了顆粒直徑為100~3000nm的著火過程,得到了不同粒徑下,顆粒溫度隨時(shí)間變化,如圖19所示.由圖可知,在著火之前,不同顆粒粒徑下,顆粒溫度隨時(shí)間變化的趨勢變化較大,且隨著粒徑增大,顆粒升溫速率減小.這主要因?yàn)榱皆酱?,顆粒存儲(chǔ)能量越多,環(huán)境對顆粒的加熱時(shí)間增加.在顆粒直徑為100nm、500nm、900nm、1500nm和3000nm的條件下,所計(jì)算的著火延遲時(shí)間分別為0.45μs、5.9μs、16μs、40μs和142μs,說明隨著粒徑增大,著火延遲時(shí)間增加.

    (a)顆粒溫度變化

    (b)圖(a)的局部放大

    圖18 不同初始氧化層厚度下顆粒溫度隨時(shí)間變化

    Fig.18 Particle temperature versus time under different initial oxide-layer thickness

    (a)顆粒溫度變化

    (b)圖(a)的局部放大

    圖19 不同顆粒直徑下顆粒溫度隨時(shí)間變化

    Fig.19 Particle temperature versus time under different particle diameter

    3 結(jié) 論

    基于縮小未反應(yīng)核模型,假設(shè)氧化層為單層Fe3O4層,并考慮氧化層對氧氣擴(kuò)散的影響,構(gòu)建了空氣中Fe顆粒的非均相著火模型.應(yīng)用該模型與Mi模型的計(jì)算結(jié)果進(jìn)行對比,分析了鐵顆粒著火過程中顆粒溫度、氧化層內(nèi)徑隨時(shí)間變化關(guān)系及著火過程中不同時(shí)刻的鐵顆粒所在空間的溫度和氧氣摩爾濃度分布.基于該著火模型,研究了環(huán)境壓強(qiáng)、氧氣濃度、初始氧化層厚度和顆粒粒徑對著火溫度、著火延遲的影響,得到如下結(jié)論:

    (1)所建立的著火模型計(jì)算的著火溫度與實(shí)驗(yàn)結(jié)果的最大誤差為6.95%,與Mi模型所計(jì)算的顆粒升溫速率基本一致,兩者著火延遲時(shí)間的數(shù)量級相同.

    (2)在鐵顆粒所處的環(huán)境空間中,隨著時(shí)間增加,在離開顆粒方向的空間溫度梯度由正梯度逐漸轉(zhuǎn)變?yōu)樨?fù)梯度,即隨著時(shí)間的推移,顆粒溫度由常溫開始迅速增大到接近環(huán)境溫度,再在反應(yīng)生成熱作用下顆粒溫度進(jìn)一步升高,直至著火成功;而氧氣摩爾濃度梯度在離開顆粒方向保持正梯度,且梯度隨時(shí)間的增加而增大.

    (3)在顆粒著火成功之前,顆粒升溫速率受環(huán)境壓強(qiáng)、氧氣濃度和初始氧化層厚度的影響較小,但隨顆粒粒徑增大而減?。?/p>

    (4)著火溫度隨環(huán)境壓強(qiáng)和氧氣濃度增大而降低,隨初始氧化層厚度增大而升高;粒徑小于10μm的范圍內(nèi),著火溫度隨粒徑減小而增大,粒徑在10~50μm范圍內(nèi),粒徑對著火溫度影響不大.著火延遲時(shí)間隨環(huán)境壓強(qiáng)、氧氣濃度增大而縮短,但隨初始氧化層厚度、顆粒粒徑增大而延長.

    1——指前因子,s-1;

    C——混合氣體和組分的摩爾濃度,mol/m3;

    和e——擴(kuò)散系數(shù)和多孔氧化層的擴(kuò)散系數(shù),m2/s;

    1——活化能,J/mol;

    ——對流換熱系數(shù),W/(K·m2);

    re1——化學(xué)反應(yīng)速率常數(shù);

    M——混合氣體和組分的相對分子(原子)質(zhì)量,g/mol;

    ——環(huán)境壓強(qiáng),Pa;

    re1——反應(yīng)熱,J/mol;

    ——通用氣體常數(shù),J/(mol·K);

    、e、b和s——半徑、顆粒內(nèi)部Fe核半徑、氧化層內(nèi)徑和顆粒半徑(氧化層外徑),m;

    ——舍伍德數(shù);

    、melt、s和T——溫度、Fe的熔點(diǎn)、顆粒溫度和無窮遠(yuǎn)處的溫度,K;

    ——速度,m/s;

    ——顆粒表面的發(fā)射率;

    p——多孔氧化層的孔隙率;

    ——導(dǎo)熱率,W/(m·K);

    ρ——混合氣體和組分的密度,g/m3;

    ——斯蒂芬-玻爾茲曼常數(shù),W/(m2·K4);

    p——多孔氧化層的曲折因子.

    角標(biāo):

    r——表示球形坐標(biāo)系中某一位置;

    s——表示顆粒表面、氧化層外表面;

    b——表示氧化層內(nèi)表面、反應(yīng)界面;

    ∞ ——表示無窮遠(yuǎn)處.

    [1] 李 悅,胡春波,胡加明,等. 粉末火箭發(fā)動(dòng)機(jī)研究進(jìn)展[J]. 推進(jìn)技術(shù),2018,39(8):1681-1695.

    Li Yue,Hu Chunbo,Hu Jiaming,et al. Progress of powder rocket engine technology[J].,2018,39(8):1681-1695(in Chinese).

    [2] Li Y,Hu C B,Deng Z,et al. Experimental study on multiple-pulse performance characteristics of ammonium perchlorate/aluminum powder rocket motor[J].,2017,133:455-466.

    [3] Szabo J,Miller T,Herr J,et al. Magnesium bipropellant rockets for martian ascent vehicles[C]//47San Diego,USA,2011:5834.

    [4] Huang H T,Zou M S,Guo X Y,et al. Analysis of the solid combustion products of a Mg-based fuel-rich propellant used for water ramjet engines[J].,2012,37(4):407-412.

    [5] Huang H T,Zou M S,Guo X Y,et al. Analysis of the aluminum reaction efficiency in a hydro-reactive fuel propellant used for a water ramjet[J].,2013,49(5):541-547.

    [6] Beach D B,Rondinone A J,Sumpter B G,et al. Solid-state combustion of metallic nanoparticles:New possibilities for an alternative energy carrier[J].,2007,129:29-32.

    [7] 王金云,王孟軍,楊在林,等. 金屬燃料技術(shù)研究進(jìn)展[C]//第五屆空天動(dòng)力聯(lián)合會(huì)議暨中國航天第三專業(yè)信息網(wǎng)第41屆技術(shù)交流會(huì)論文集(第三冊). 南京,2020:206-219.

    Wang Jinyun,Wang Mengjun,Yang Zailin,et al. Re-search progress of metal fuel technology[C]//41st(3)Nanjing,2020:206-219(in Chinese).

    [8] Bergthorson J M,Goroshin S,Soo M J,et al. Direct combustion of recyclable metal fuels for zero-carbon heat and power[J].,2015,160:368-382.

    [9] Bergthorson J M. Recyclable metal fuels for clean and compact zero-carbon power[J].,2018,68:169-196.

    [10] Mandilas C,Karagiannakis G,Konstandopoulos A G,et al. Study of oxidation and combustion characteristics of iron nanoparticles under idealized and enginelike conditions[J].,2016,30(5):4318-4330.

    [11] 劉世寧,胡春波,胡 旭,等. 基于產(chǎn)物分析的鐵粉燃燒器燃燒特性實(shí)驗(yàn)研究[J]. 推進(jìn)技術(shù),2022,43(6):385-396.

    Liu Shining,Hu Chunbo,Hu Xu,et,al. Experimental study on combustion characteristics of iron powder burner based on product analysis[J].,2022,43(6):385-396 (in Chinese).

    [12] 葉大倫,胡建華. 實(shí)用無機(jī)物熱力學(xué)數(shù)據(jù)手冊[M]. 第2版. 北京:冶金工業(yè)出版社,2002.

    Ye Dalun,Hu Jianhua.[M]. 2nd Edition. Beijing:Metallurgical Industry Press,2002(in Chinese).

    [13] Wilson D B,Steinberg T A,Stoltzfus J M. Thermodynamics and kinetics of burning iron[J].,1997,1319:240-257.

    [14] 高文靜,金 晶,曾武勇. 納米鐵粉的燃燒動(dòng)力學(xué)模型研究[J]. 科學(xué)技術(shù)與工程,2013,13(33):9808-9812.

    Gao Wenjing,Jin Jing,Zeng Wuyong. Kinetic model study on combustion of nano iron powders[J].,2013,13(33):9808-9812 (in Chinese).

    [15] 孫金華,盧 平,劉 義. 空氣中懸浮金屬微粒子的燃燒特性[J]. 南京理工大學(xué)學(xué)報(bào),2005,29(5):582-585.

    Sun Jinhua,Lu Ping,Liu Yi. Combustion behavior of metal particles suspended in air[J].,2005,29(5):582-585 (in Chinese).

    [16] Li S,Huang J Q,Weng W B,et al. Ignition and combustion behavior of single micron-sized iron particle in hot gas flow[J].,2022,241:112099.

    [17] Leshchevich V V,Penyazkov O G,F(xiàn)edorov A V,et al. Conditions and delay time of ignition of iron microparticles in oxygen[J].,2012,85(1):148-154.

    [18] 丁宋毅,孫 波,卓長飛. 基于粒子群算法的微納米鐵粉燃燒實(shí)驗(yàn)研究[J]. 科學(xué)技術(shù)與工程,2022,22(7):2709-2716.

    Ding Songyi,Sun Bo,Zhuo Changfei. Experimental research on combustion of micro/nano iron powder based on particle swarm algorithm[J].,2022,22(7):2709-2716(in Chinese).

    [19] Schiemann M,F(xiàn)ischer P,Bergthorson J M. Iron particles as carbon-neutral fuel in spray roasting reactors[C]//8. Dubrovnik,Croatia,2017:487-492.

    [20] Mi X C,F(xiàn)ujinawa A,Bergthorson J M. A quantitative analysis of the ignition characteristics of fine iron particles[J].,2022,240:112011.

    [21] 肖興國,謝蘊(yùn)國. 冶金反應(yīng)工程學(xué)基礎(chǔ)[M]. 北京:冶金工業(yè)出版社,1997.

    Xiao Xingguo,Xie Yunguo.[M]. Beijing:Metallurgy Industry Press,1997(in Chinese).

    [22] Gleeson B,Hadavi S M M,Young D J. Isothermal transformation behavior of thermally-grown wüstite[J].,2000,17(2):311-318.

    [23] Juricic C. On the Mechanisms of Internal Stress Formation in Multi-phase Iron Oxide Scales[D]. Bochum:Ruhr University Bochum,2008.

    [24] Chen R Y,Yeun W Y D. Review of the high-temperature oxidation of iron and carbon steels in air or oxygen[J].,2003,59(5):433-468.

    [25] Pa?dassi J. The kinetics of the air oxidation of iron in the range 700—1250℃[J].,1958,6(3):184-194.

    [26] Neil B,Gerald H M,F(xiàn)rederick S P.[M]. 2nd EditionCambridge:Cambridge University Press,2006.

    [27] 李鐵藩. 金屬高溫氧化和熱腐蝕[M]. 北京:化學(xué)工業(yè)出版社,2003.

    Li Tiefan.[M]. Beijing:Chemical Industry Press,2003(in Chinese).

    [28] Chen R Y,Yuen W Y D. A study of the scale structure of hot-rolled steel strip by simulated coiling and cooling[J].,2000,53(5):539-560.

    [29] 李文超. 冶金與材料物理化學(xué)[M]. 北京:冶金工業(yè)出版社,2001.

    Li Wenchao.[M]. Beijing:Metallurgical Industry Press,2001(in Chinese).

    [30] 葛慶仁. 氣固反應(yīng)動(dòng)力學(xué)[M]. 北京:原子能出版社,1991.

    Ge Qingren.[M]. Beijing:Atomic Energy Press,1991(in Chinese).

    [31] 李法社,王 華. 高等燃燒學(xué)[M]. 北京:科學(xué)出版社,2016.

    Li Fashe,Wang Hua.[M]. Beijing:Science Press,2016(in Chinese).

    [32] Arshad U,Taqvi S A A,Buang A,et al. SVM,ANN,and PSF modelling approaches for prediction of iron dust minimum ignition temperature (MIT)based on the synergistic effect of dispersion pressure and concen-tration[J].,2021,152:375-390.

    [33] Liu Y X,Liu D,Liu G N. Energy conversion and ignition of iron nanoparticles by flash[J].,2017,60(12):1878-1884.

    [34] Khaikin B I,Bloshenko V N,Merzhanov A G. On the ignition of metal particles[J].,1970,6(4):412-422.

    [35] 楊 麗,朱燕群,王智化,等. 微納米金屬鐵粉的燃燒特性試驗(yàn)研究[J]. 浙江大學(xué)學(xué)報(bào)(工學(xué)版),2010,44(8):1562-1566.

    Yang Li,Zhu Yanqun,Wang Zhihua,et al. Experimental study of combustion characteristics of micron and nano iron powders[J].(),2010,44(8):1562-1566 (in Chinese).

    Heterogeneous Ignition Model and Ignition Characteristics of Iron Particle

    Cai Ping,Xu Yihua,Hu Kun,Chen Xuanren,Yang Zhihao,Sun Haijun,Ren Guanlong

    (Jiangxi Key Laboratory of Micro Aeroengine,School of Aircraft Engineering,Nanchang Hangkong University,Nanchang 330063,China)

    The heterogeneous ignition model of iron particles in air was established based on shrinking unreacted-core model,taking into account the influence of oxygen diffusion in the oxide layer on the ignition process. The maximum error between the calculated ignition temperature by this model and the experimental ignition temperature is 6.95%. Compared with Mi's ignition model,the particle heating rates calculated by the two models are basically the same,and their ignition delay time is of the same order of magnitude. The influences of atmospheric pressure,oxygen concentration,initial oxide-layer thickness and particle diameter on ignition temperature and ignition delay were studied based on this model. The results showed that ignition temperature and ignition delay decrease with the increase of atmospheric pressure and oxygen concentration but increase with the increase of initial oxide-layer thickness. The ignition temperature within the range of particle diameter less than 10μm increases with the decrease of particle diameter,and within the range of 10μm to 50μm,it is hardly affected by particle diameter. The ignition delay time increases as the particle diameter increases.

    iron particle;fuel;heterogeneous ignition model;ignition temperature;ignition delay

    TK16

    A

    1006-8740(2024)01-0091-12

    2023-03-01.

    國家自然科學(xué)基金資助項(xiàng)目(51666012;12102161);航空科學(xué)基金資助項(xiàng)目(20200001056001).

    蔡 平(1998— ),男,碩士研究生,caipinghead@163.com.

    徐義華,男,博士,教授,xuyihua@nchu.edu.cn.

    (責(zé)任編輯:隋韶穎)

    猜你喜歡
    鐵粉氧氣粒徑
    火星上成功制造出氧氣啦
    軍事文摘(2023年22期)2023-12-19 06:41:04
    聚焦空氣與氧氣
    SiO2包覆羰基鐵粉及其涂層的耐腐蝕性能
    氧氣的測定與制取
    木屑粒徑對黑木耳栽培的影響試驗(yàn)*
    天然微合金鐵粉中釩、鈦、鉻含量測定的研究
    昆鋼科技(2020年6期)2020-03-29 06:39:42
    基于近場散射的顆粒粒徑分布測量
    地球上的氧氣能用得完嗎?
    Oslo結(jié)晶器晶體粒徑分布特征的CFD模擬
    SAPO-56分子篩的形貌和粒徑控制
    国产精品 国内视频| 美女主播在线视频| 亚洲精品国产一区二区精华液| 美女福利国产在线| 天天躁夜夜躁狠狠躁躁| 伦理电影大哥的女人| 国产黄频视频在线观看| 久久久精品94久久精品| 美女国产高潮福利片在线看| 亚洲国产看品久久| 欧美xxⅹ黑人| 一本—道久久a久久精品蜜桃钙片| 91精品国产国语对白视频| 韩国av在线不卡| 精品人妻在线不人妻| 免费观看av网站的网址| 久久久久视频综合| 中文字幕精品免费在线观看视频| 丁香六月天网| 制服丝袜香蕉在线| 午夜久久久在线观看| 亚洲国产最新在线播放| 一个人免费看片子| 在线看a的网站| e午夜精品久久久久久久| 久久韩国三级中文字幕| 伊人久久国产一区二区| 欧美激情高清一区二区三区 | 国产精品久久久久久精品电影小说| 亚洲国产精品国产精品| 国产成人精品在线电影| 亚洲av男天堂| 日本欧美国产在线视频| 免费黄网站久久成人精品| 亚洲伊人色综图| 99国产精品免费福利视频| 国产欧美亚洲国产| 别揉我奶头~嗯~啊~动态视频 | 2018国产大陆天天弄谢| 国产黄色视频一区二区在线观看| 午夜福利一区二区在线看| 亚洲av欧美aⅴ国产| 国产男人的电影天堂91| 午夜激情av网站| 日韩视频在线欧美| 你懂的网址亚洲精品在线观看| 国产精品偷伦视频观看了| 国产日韩欧美在线精品| 免费黄网站久久成人精品| 自线自在国产av| 男人舔女人的私密视频| 亚洲专区中文字幕在线 | 永久免费av网站大全| 最新的欧美精品一区二区| 黄片无遮挡物在线观看| 香蕉国产在线看| 国产一区二区在线观看av| 午夜老司机福利片| 高清av免费在线| 国产成人a∨麻豆精品| 另类精品久久| 国产精品麻豆人妻色哟哟久久| 免费黄频网站在线观看国产| 黄色视频不卡| av.在线天堂| 电影成人av| 国产精品无大码| 欧美中文综合在线视频| 日韩欧美精品免费久久| 少妇人妻精品综合一区二区| 亚洲成人免费av在线播放| 精品一品国产午夜福利视频| 考比视频在线观看| 男女床上黄色一级片免费看| 黄片小视频在线播放| 亚洲av电影在线进入| 午夜福利,免费看| 最新的欧美精品一区二区| 青春草亚洲视频在线观看| 秋霞在线观看毛片| 亚洲久久久国产精品| 视频在线观看一区二区三区| 国产一区二区在线观看av| 午夜影院在线不卡| 精品亚洲成a人片在线观看| 秋霞伦理黄片| 热99久久久久精品小说推荐| 不卡视频在线观看欧美| 免费女性裸体啪啪无遮挡网站| 久久免费观看电影| 性少妇av在线| 成年人午夜在线观看视频| 亚洲国产看品久久| 亚洲熟女精品中文字幕| 欧美日韩国产mv在线观看视频| 无遮挡黄片免费观看| 男女之事视频高清在线观看 | 久久这里只有精品19| av卡一久久| av线在线观看网站| 天天操日日干夜夜撸| 国产男女内射视频| 97精品久久久久久久久久精品| 精品一品国产午夜福利视频| 91国产中文字幕| 亚洲精品久久午夜乱码| 精品亚洲成a人片在线观看| 中文字幕人妻丝袜制服| 国产片内射在线| 一级黄片播放器| 黄片小视频在线播放| 亚洲色图 男人天堂 中文字幕| xxxhd国产人妻xxx| 国产亚洲av片在线观看秒播厂| 成人亚洲欧美一区二区av| 免费不卡黄色视频| 狂野欧美激情性bbbbbb| av有码第一页| 精品午夜福利在线看| 老司机影院成人| 午夜福利乱码中文字幕| 日本色播在线视频| 婷婷色综合www| 久久精品熟女亚洲av麻豆精品| 国产精品久久久久久人妻精品电影 | 午夜福利视频精品| 国产又爽黄色视频| 日韩大码丰满熟妇| 9191精品国产免费久久| 亚洲人成网站在线观看播放| 国产欧美日韩综合在线一区二区| 五月天丁香电影| 亚洲视频免费观看视频| 精品国产一区二区久久| 国产精品熟女久久久久浪| av视频免费观看在线观看| 如日韩欧美国产精品一区二区三区| av在线播放精品| 国产野战对白在线观看| 免费高清在线观看视频在线观看| 久久精品国产综合久久久| 在线观看一区二区三区激情| 欧美精品一区二区免费开放| 免费在线观看视频国产中文字幕亚洲 | 亚洲自偷自拍图片 自拍| 国产精品香港三级国产av潘金莲 | 亚洲免费av在线视频| 只有这里有精品99| 国产毛片在线视频| 久久女婷五月综合色啪小说| 成人国产麻豆网| 美国免费a级毛片| 少妇精品久久久久久久| 欧美最新免费一区二区三区| 免费人妻精品一区二区三区视频| 久久久久人妻精品一区果冻| 国产精品99久久99久久久不卡 | 老鸭窝网址在线观看| 2018国产大陆天天弄谢| 国产老妇伦熟女老妇高清| 日日啪夜夜爽| 波野结衣二区三区在线| 成人国产麻豆网| 国产1区2区3区精品| 国产精品久久久久久久久免| 国产精品.久久久| 精品国产一区二区久久| 9191精品国产免费久久| 极品少妇高潮喷水抽搐| 制服丝袜香蕉在线| 激情五月婷婷亚洲| 下体分泌物呈黄色| 免费黄网站久久成人精品| 久久久久精品人妻al黑| 欧美日韩成人在线一区二区| 欧美日韩福利视频一区二区| 亚洲国产中文字幕在线视频| 日日爽夜夜爽网站| 久久人人爽人人片av| 精品一区二区三卡| 精品国产一区二区三区久久久樱花| 国产国语露脸激情在线看| 日日摸夜夜添夜夜爱| 免费黄频网站在线观看国产| 亚洲天堂av无毛| 国语对白做爰xxxⅹ性视频网站| 考比视频在线观看| 十分钟在线观看高清视频www| 丝袜脚勾引网站| 亚洲国产日韩一区二区| svipshipincom国产片| 国产精品蜜桃在线观看| 亚洲av国产av综合av卡| 日本欧美国产在线视频| 精品少妇黑人巨大在线播放| 18禁裸乳无遮挡动漫免费视频| 亚洲精品第二区| 国产精品.久久久| 国产激情久久老熟女| 久久久精品区二区三区| 香蕉丝袜av| 成人影院久久| 久久国产精品男人的天堂亚洲| 免费高清在线观看视频在线观看| 欧美人与性动交α欧美软件| 欧美精品亚洲一区二区| 午夜福利在线免费观看网站| 国产精品女同一区二区软件| 成年动漫av网址| 日本色播在线视频| 亚洲精品第二区| 天天添夜夜摸| 日韩成人av中文字幕在线观看| videos熟女内射| 人成视频在线观看免费观看| 亚洲精品av麻豆狂野| 亚洲国产看品久久| 久久精品久久久久久噜噜老黄| 下体分泌物呈黄色| 你懂的网址亚洲精品在线观看| 国产精品女同一区二区软件| 久久精品aⅴ一区二区三区四区| 成人国产麻豆网| 亚洲欧美精品综合一区二区三区| 夫妻性生交免费视频一级片| 观看美女的网站| av在线老鸭窝| 亚洲欧美成人综合另类久久久| 欧美国产精品va在线观看不卡| 男女床上黄色一级片免费看| 激情五月婷婷亚洲| 美女午夜性视频免费| 99香蕉大伊视频| 亚洲一区中文字幕在线| 成人亚洲欧美一区二区av| 青春草亚洲视频在线观看| 国产亚洲精品第一综合不卡| 少妇被粗大猛烈的视频| 激情视频va一区二区三区| 亚洲第一区二区三区不卡| 青青草视频在线视频观看| tube8黄色片| 国产福利在线免费观看视频| 中国国产av一级| 欧美精品高潮呻吟av久久| 满18在线观看网站| 成人18禁高潮啪啪吃奶动态图| 涩涩av久久男人的天堂| 视频区图区小说| 麻豆乱淫一区二区| 人成视频在线观看免费观看| 嫩草影院入口| 蜜桃在线观看..| 中国三级夫妇交换| 人人妻,人人澡人人爽秒播 | 久久99精品国语久久久| 最近2019中文字幕mv第一页| 久久精品国产亚洲av高清一级| xxxhd国产人妻xxx| 精品人妻在线不人妻| 七月丁香在线播放| 久久久久久人妻| 亚洲精品自拍成人| 国产老妇伦熟女老妇高清| 午夜福利免费观看在线| 丁香六月欧美| 亚洲视频免费观看视频| 多毛熟女@视频| 无遮挡黄片免费观看| 少妇被粗大猛烈的视频| 欧美亚洲日本最大视频资源| 这个男人来自地球电影免费观看 | 亚洲人成电影观看| 97人妻天天添夜夜摸| 国产精品av久久久久免费| 欧美日韩视频高清一区二区三区二| 久久精品熟女亚洲av麻豆精品| 亚洲伊人久久精品综合| 999久久久国产精品视频| 免费高清在线观看视频在线观看| 最近的中文字幕免费完整| 狠狠婷婷综合久久久久久88av| 啦啦啦中文免费视频观看日本| 秋霞伦理黄片| 又粗又硬又长又爽又黄的视频| 国产在线一区二区三区精| 久久性视频一级片| 97精品久久久久久久久久精品| 黄频高清免费视频| 亚洲欧美清纯卡通| 丝袜美腿诱惑在线| 国产1区2区3区精品| 国产毛片在线视频| 亚洲国产精品一区三区| 免费看不卡的av| 成年女人毛片免费观看观看9 | 日韩精品免费视频一区二区三区| 日韩中文字幕视频在线看片| 在线观看人妻少妇| 在线观看免费午夜福利视频| 一本大道久久a久久精品| av视频免费观看在线观看| 人人妻人人澡人人看| 日本欧美国产在线视频| 女性生殖器流出的白浆| 国产精品欧美亚洲77777| 女人被躁到高潮嗷嗷叫费观| 久久精品国产亚洲av高清一级| 国产精品女同一区二区软件| 久久久久久久大尺度免费视频| 久久精品人人爽人人爽视色| 亚洲国产欧美在线一区| 一级a爱视频在线免费观看| 在线精品无人区一区二区三| 日韩制服骚丝袜av| 曰老女人黄片| 天天操日日干夜夜撸| 一本大道久久a久久精品| 人人妻人人澡人人爽人人夜夜| 人人澡人人妻人| 考比视频在线观看| 国产亚洲最大av| 亚洲精品久久成人aⅴ小说| 成年女人毛片免费观看观看9 | 国产av一区二区精品久久| 国产精品人妻久久久影院| 新久久久久国产一级毛片| 亚洲欧洲国产日韩| 老熟女久久久| 免费不卡黄色视频| 丁香六月天网| 中国三级夫妇交换| 在线观看一区二区三区激情| 自拍欧美九色日韩亚洲蝌蚪91| 在线观看人妻少妇| 成年女人毛片免费观看观看9 | 国产精品国产av在线观看| 中文精品一卡2卡3卡4更新| 涩涩av久久男人的天堂| 各种免费的搞黄视频| 亚洲伊人色综图| 日韩不卡一区二区三区视频在线| 99久久人妻综合| bbb黄色大片| 三上悠亚av全集在线观看| 91aial.com中文字幕在线观看| 麻豆精品久久久久久蜜桃| 日韩av在线免费看完整版不卡| 国产片特级美女逼逼视频| 男女之事视频高清在线观看 | 亚洲欧美激情在线| 街头女战士在线观看网站| 国产乱来视频区| 日韩一本色道免费dvd| 97人妻天天添夜夜摸| 欧美精品一区二区大全| 97人妻天天添夜夜摸| 人成视频在线观看免费观看| 最近2019中文字幕mv第一页| 亚洲人成电影观看| 婷婷成人精品国产| 国产极品粉嫩免费观看在线| 日韩制服丝袜自拍偷拍| 亚洲国产av新网站| 丁香六月天网| 啦啦啦在线观看免费高清www| tube8黄色片| 爱豆传媒免费全集在线观看| 精品一区二区免费观看| 亚洲国产精品国产精品| 亚洲精品乱久久久久久| 国产免费又黄又爽又色| 美女扒开内裤让男人捅视频| 午夜免费男女啪啪视频观看| 中国国产av一级| 亚洲国产精品999| 成人漫画全彩无遮挡| 久久av网站| 久久午夜综合久久蜜桃| 精品人妻熟女毛片av久久网站| 综合色丁香网| 亚洲中文av在线| 午夜福利视频精品| 女性生殖器流出的白浆| 天天躁日日躁夜夜躁夜夜| 香蕉丝袜av| 亚洲自偷自拍图片 自拍| 精品免费久久久久久久清纯 | 精品一品国产午夜福利视频| 欧美精品av麻豆av| 99国产精品免费福利视频| 看免费av毛片| 老司机影院毛片| 咕卡用的链子| 亚洲国产精品一区二区三区在线| 国产精品久久久久久精品电影小说| 香蕉丝袜av| 日本欧美国产在线视频| 一级,二级,三级黄色视频| 侵犯人妻中文字幕一二三四区| 日韩av在线免费看完整版不卡| 日韩一区二区三区影片| 一区二区三区精品91| 精品少妇一区二区三区视频日本电影 | 久久久久久人妻| 乱人伦中国视频| 男的添女的下面高潮视频| 亚洲婷婷狠狠爱综合网| 人妻一区二区av| 性少妇av在线| 悠悠久久av| 成人毛片60女人毛片免费| 男男h啪啪无遮挡| 99re6热这里在线精品视频| 国产片内射在线| tube8黄色片| 香蕉国产在线看| 咕卡用的链子| 人人妻人人爽人人添夜夜欢视频| 欧美人与性动交α欧美精品济南到| 精品少妇内射三级| 宅男免费午夜| 亚洲精品乱久久久久久| 老汉色∧v一级毛片| 国产福利在线免费观看视频| 毛片一级片免费看久久久久| 尾随美女入室| 一边亲一边摸免费视频| www.自偷自拍.com| 亚洲精品日韩在线中文字幕| 高清欧美精品videossex| 国产一区二区激情短视频 | 日本黄色日本黄色录像| 欧美少妇被猛烈插入视频| 国产野战对白在线观看| 国产精品国产三级专区第一集| 欧美国产精品va在线观看不卡| 精品酒店卫生间| 国产成人a∨麻豆精品| 你懂的网址亚洲精品在线观看| 国产精品久久久久久精品电影小说| 18在线观看网站| 老司机靠b影院| 亚洲美女视频黄频| 美女午夜性视频免费| 国产视频首页在线观看| 欧美激情极品国产一区二区三区| 中文字幕精品免费在线观看视频| 麻豆av在线久日| 老鸭窝网址在线观看| 美女国产高潮福利片在线看| 中文欧美无线码| 三上悠亚av全集在线观看| 最新在线观看一区二区三区 | 欧美日本中文国产一区发布| 在线观看免费高清a一片| 中文字幕人妻丝袜制服| 秋霞在线观看毛片| 视频区图区小说| 国产精品久久久av美女十八| 在线观看国产h片| 亚洲精品美女久久av网站| 99九九在线精品视频| 在线天堂最新版资源| 国产精品久久久久久精品电影小说| 久久99精品国语久久久| 乱人伦中国视频| 大片电影免费在线观看免费| 人人妻人人爽人人添夜夜欢视频| netflix在线观看网站| 纵有疾风起免费观看全集完整版| 亚洲成人一二三区av| 午夜老司机福利片| 久久精品熟女亚洲av麻豆精品| 久久综合国产亚洲精品| 波野结衣二区三区在线| 国产精品99久久99久久久不卡 | 久久久国产欧美日韩av| 欧美另类一区| 伦理电影大哥的女人| 日本欧美国产在线视频| 成人手机av| 十八禁高潮呻吟视频| 又粗又硬又长又爽又黄的视频| 亚洲国产精品999| 王馨瑶露胸无遮挡在线观看| 中文天堂在线官网| 老汉色∧v一级毛片| e午夜精品久久久久久久| 在线观看免费视频网站a站| 日本av手机在线免费观看| 精品国产一区二区三区久久久樱花| 黄色视频在线播放观看不卡| av免费观看日本| 亚洲国产毛片av蜜桃av| 黄片播放在线免费| 最近中文字幕2019免费版| 国产精品免费视频内射| 日韩不卡一区二区三区视频在线| 在线观看免费视频网站a站| 久久青草综合色| 国产熟女欧美一区二区| 亚洲熟女毛片儿| 亚洲欧美成人综合另类久久久| 精品一区二区三卡| 啦啦啦中文免费视频观看日本| av.在线天堂| 黄网站色视频无遮挡免费观看| 又大又黄又爽视频免费| 欧美黑人精品巨大| 在线观看免费高清a一片| 欧美在线一区亚洲| 女人被躁到高潮嗷嗷叫费观| 午夜激情久久久久久久| 看免费成人av毛片| 一边摸一边抽搐一进一出视频| 成人黄色视频免费在线看| 国产亚洲精品第一综合不卡| 亚洲国产精品一区三区| 一级毛片黄色毛片免费观看视频| 国产亚洲一区二区精品| 亚洲精品乱久久久久久| 一本一本久久a久久精品综合妖精| 亚洲婷婷狠狠爱综合网| 国产黄色免费在线视频| 伊人久久国产一区二区| 制服人妻中文乱码| 夫妻性生交免费视频一级片| 国产亚洲av高清不卡| 久久这里只有精品19| 亚洲人成网站在线观看播放| netflix在线观看网站| 亚洲精品日韩在线中文字幕| 精品人妻一区二区三区麻豆| 精品第一国产精品| 狠狠婷婷综合久久久久久88av| 激情视频va一区二区三区| 十分钟在线观看高清视频www| a级片在线免费高清观看视频| 精品少妇内射三级| 婷婷色麻豆天堂久久| 97在线人人人人妻| 久久午夜综合久久蜜桃| 一级毛片黄色毛片免费观看视频| 妹子高潮喷水视频| 啦啦啦在线观看免费高清www| 在线天堂中文资源库| 曰老女人黄片| 少妇人妻 视频| 另类亚洲欧美激情| 王馨瑶露胸无遮挡在线观看| 十八禁人妻一区二区| 丝袜喷水一区| 日韩成人av中文字幕在线观看| 亚洲国产精品一区二区三区在线| 欧美日韩亚洲高清精品| 欧美日韩亚洲国产一区二区在线观看 | av卡一久久| 视频在线观看一区二区三区| 精品亚洲成a人片在线观看| 亚洲情色 制服丝袜| 制服丝袜香蕉在线| 中文字幕制服av| 夫妻午夜视频| 人人妻人人澡人人爽人人夜夜| 一边摸一边做爽爽视频免费| 日日撸夜夜添| 一本一本久久a久久精品综合妖精| 不卡av一区二区三区| 99精品久久久久人妻精品| 日韩欧美一区视频在线观看| 午夜福利视频精品| 国产精品国产三级专区第一集| 波多野结衣av一区二区av| 天美传媒精品一区二区| 国产精品久久久久久人妻精品电影 | 亚洲第一区二区三区不卡| 悠悠久久av| 黑人猛操日本美女一级片| 欧美日韩综合久久久久久| 国产精品一二三区在线看| av又黄又爽大尺度在线免费看| www.av在线官网国产| 免费黄色在线免费观看| 18禁裸乳无遮挡动漫免费视频| 国产一区二区三区综合在线观看| 国产精品99久久99久久久不卡 | 亚洲美女搞黄在线观看| 成人国语在线视频| 久久ye,这里只有精品| 日本欧美视频一区| 久久久精品94久久精品| 欧美黄色片欧美黄色片| 丁香六月欧美| 天堂中文最新版在线下载| 黄网站色视频无遮挡免费观看| 少妇人妻 视频| 秋霞伦理黄片| 日日摸夜夜添夜夜爱| 久久99一区二区三区| 超碰97精品在线观看| 一本—道久久a久久精品蜜桃钙片| 午夜免费男女啪啪视频观看| 在线观看免费高清a一片| 人成视频在线观看免费观看| 欧美 日韩 精品 国产| xxx大片免费视频| 十八禁网站网址无遮挡| 久久精品亚洲熟妇少妇任你| 国产精品成人在线| 97在线人人人人妻| 国产亚洲最大av| 国产高清不卡午夜福利| 少妇人妻精品综合一区二区| 国产精品秋霞免费鲁丝片|