• <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分子篩的形貌和粒徑控制
    亚洲国产日韩欧美精品在线观看 | 一边摸一边做爽爽视频免费| 色综合欧美亚洲国产小说| 国产精品二区激情视频| 999久久久精品免费观看国产| 婷婷丁香在线五月| 两性夫妻黄色片| 日日干狠狠操夜夜爽| 99热只有精品国产| 19禁男女啪啪无遮挡网站| 日韩精品青青久久久久久| 黄网站色视频无遮挡免费观看| 国产又黄又爽又无遮挡在线| 黄频高清免费视频| 免费在线观看完整版高清| 精品电影一区二区在线| 超碰成人久久| 50天的宝宝边吃奶边哭怎么回事| 人成视频在线观看免费观看| 成人亚洲精品一区在线观看| 欧美激情 高清一区二区三区| 欧美日本亚洲视频在线播放| 侵犯人妻中文字幕一二三四区| 欧美乱码精品一区二区三区| 婷婷丁香在线五月| 99riav亚洲国产免费| 母亲3免费完整高清在线观看| 亚洲av电影在线进入| 少妇粗大呻吟视频| 成年版毛片免费区| 亚洲国产日韩欧美精品在线观看 | 午夜久久久在线观看| 国产97色在线日韩免费| 可以在线观看的亚洲视频| 丝袜人妻中文字幕| 成人午夜高清在线视频 | 亚洲一卡2卡3卡4卡5卡精品中文| 69av精品久久久久久| 日韩三级视频一区二区三区| 男女下面进入的视频免费午夜 | 中文字幕人妻丝袜一区二区| 午夜福利18| 我的亚洲天堂| 亚洲成国产人片在线观看| 欧美在线一区亚洲| 日韩三级视频一区二区三区| 国产精品久久久人人做人人爽| 午夜两性在线视频| 91麻豆av在线| 国产单亲对白刺激| 午夜精品久久久久久毛片777| 国产单亲对白刺激| 久久久精品国产亚洲av高清涩受| 男人操女人黄网站| 男人的好看免费观看在线视频 | 精华霜和精华液先用哪个| 国产99久久九九免费精品| x7x7x7水蜜桃| 成人亚洲精品一区在线观看| 香蕉av资源在线| 欧美一区二区精品小视频在线| 一级黄色大片毛片| 天天一区二区日本电影三级| 女性被躁到高潮视频| 一区二区三区国产精品乱码| 成年人黄色毛片网站| 嫩草影院精品99| 亚洲成人国产一区在线观看| 在线观看www视频免费| 久久久久免费精品人妻一区二区 | 精品一区二区三区av网在线观看| 天堂影院成人在线观看| 丰满人妻熟妇乱又伦精品不卡| 午夜免费激情av| 国产精品自产拍在线观看55亚洲| 国产不卡一卡二| 国产三级在线视频| 悠悠久久av| 国产视频一区二区在线看| 精品福利观看| 看片在线看免费视频| 美女扒开内裤让男人捅视频| 熟女少妇亚洲综合色aaa.| x7x7x7水蜜桃| 亚洲一区高清亚洲精品| 两人在一起打扑克的视频| 在线观看免费日韩欧美大片| 老熟妇仑乱视频hdxx| 国产熟女xx| 啦啦啦韩国在线观看视频| 欧美日韩亚洲国产一区二区在线观看| 中文字幕人妻熟女乱码| 制服诱惑二区| 欧美色视频一区免费| 久9热在线精品视频| 欧美乱妇无乱码| 又大又爽又粗| 亚洲精品中文字幕一二三四区| 久久精品国产亚洲av高清一级| 亚洲黑人精品在线| 日本 欧美在线| 黄频高清免费视频| av有码第一页| 国产精品av久久久久免费| 日本 欧美在线| 后天国语完整版免费观看| 欧美日韩福利视频一区二区| 亚洲成a人片在线一区二区| 久久草成人影院| 午夜成年电影在线免费观看| 亚洲欧洲精品一区二区精品久久久| 欧美大码av| 国产人伦9x9x在线观看| 夜夜夜夜夜久久久久| 一级片免费观看大全| 色精品久久人妻99蜜桃| 18美女黄网站色大片免费观看| 成人午夜高清在线视频 | 欧美日韩中文字幕国产精品一区二区三区| 18禁国产床啪视频网站| 免费电影在线观看免费观看| 久久中文看片网| 琪琪午夜伦伦电影理论片6080| 丝袜美腿诱惑在线| 精品国产乱码久久久久久男人| 久久中文字幕人妻熟女| 亚洲国产精品成人综合色| 日韩高清综合在线| 1024手机看黄色片| 国产成人av教育| 亚洲第一欧美日韩一区二区三区| 亚洲欧美精品综合久久99| 国产片内射在线| 国产精品电影一区二区三区| 午夜久久久在线观看| 国产成人影院久久av| 高清毛片免费观看视频网站| 精品久久久久久久末码| 亚洲aⅴ乱码一区二区在线播放 | 国产av在哪里看| 日韩成人在线观看一区二区三区| 午夜福利18| 淫秽高清视频在线观看| 午夜亚洲福利在线播放| 亚洲中文av在线| 嫩草影院精品99| a级毛片在线看网站| 在线观看免费日韩欧美大片| 精品一区二区三区四区五区乱码| 制服诱惑二区| 99热6这里只有精品| 精品熟女少妇八av免费久了| 91麻豆av在线| 中文在线观看免费www的网站 | 禁无遮挡网站| 久久久久九九精品影院| 极品教师在线免费播放| 亚洲一码二码三码区别大吗| 欧美成狂野欧美在线观看| 午夜老司机福利片| 香蕉国产在线看| 欧美成人一区二区免费高清观看 | 中文资源天堂在线| 婷婷亚洲欧美| 999久久久国产精品视频| 一a级毛片在线观看| 19禁男女啪啪无遮挡网站| 亚洲av第一区精品v没综合| 男人操女人黄网站| 狂野欧美激情性xxxx| 久久久国产成人免费| av超薄肉色丝袜交足视频| 国产又黄又爽又无遮挡在线| 桃色一区二区三区在线观看| 日本三级黄在线观看| 精品久久久久久久人妻蜜臀av| 99国产精品一区二区三区| 好看av亚洲va欧美ⅴa在| 精品国产乱子伦一区二区三区| 亚洲一区二区三区不卡视频| 久久精品91蜜桃| 啦啦啦免费观看视频1| 在线天堂中文资源库| 99久久无色码亚洲精品果冻| 99精品久久久久人妻精品| 午夜福利成人在线免费观看| 一级毛片精品| 欧美 亚洲 国产 日韩一| 精品一区二区三区视频在线观看免费| 最好的美女福利视频网| 久久久久久久久久黄片| 国产成人精品久久二区二区91| a级毛片在线看网站| 国内精品久久久久久久电影| 最新美女视频免费是黄的| 欧美性长视频在线观看| 伊人久久大香线蕉亚洲五| 久久精品夜夜夜夜夜久久蜜豆 | 999精品在线视频| 一级a爱片免费观看的视频| 久久精品aⅴ一区二区三区四区| 亚洲国产看品久久| 老汉色av国产亚洲站长工具| 久久久国产精品麻豆| 久久天躁狠狠躁夜夜2o2o| 免费在线观看成人毛片| 国产在线精品亚洲第一网站| 国产免费男女视频| 精品少妇一区二区三区视频日本电影| 少妇熟女aⅴ在线视频| 麻豆一二三区av精品| 9191精品国产免费久久| 亚洲人成网站高清观看| 精品国产国语对白av| 亚洲国产高清在线一区二区三 | 日本 av在线| 琪琪午夜伦伦电影理论片6080| 最近最新中文字幕大全免费视频| 日本撒尿小便嘘嘘汇集6| 久久这里只有精品19| 亚洲精品久久国产高清桃花| 久久久久国内视频| 真人一进一出gif抽搐免费| 黄色片一级片一级黄色片| 国产区一区二久久| 精品第一国产精品| 一本综合久久免费| 丝袜在线中文字幕| 久久精品国产亚洲av香蕉五月| 国产午夜精品久久久久久| 欧美午夜高清在线| 亚洲第一青青草原| 国产亚洲欧美在线一区二区| 亚洲av电影在线进入| 国产亚洲欧美在线一区二区| 国产亚洲av嫩草精品影院| 久久这里只有精品19| 观看免费一级毛片| 欧美不卡视频在线免费观看 | 国产1区2区3区精品| 欧美乱码精品一区二区三区| 国产视频一区二区在线看| 午夜a级毛片| 国产精品九九99| 99国产精品一区二区蜜桃av| 嫩草影院精品99| 一级片免费观看大全| 无人区码免费观看不卡| 国产精品久久电影中文字幕| 在线看三级毛片| 国产视频内射| 丝袜人妻中文字幕| 久久性视频一级片| 99久久国产精品久久久| 国产片内射在线| 又黄又粗又硬又大视频| 久久精品人妻少妇| av欧美777| av在线天堂中文字幕| 可以在线观看的亚洲视频| 男女那种视频在线观看| 久久久久久免费高清国产稀缺| 一二三四社区在线视频社区8| 禁无遮挡网站| 欧美成人性av电影在线观看| 久久久精品欧美日韩精品| 大型av网站在线播放| 欧美日韩一级在线毛片| 最近在线观看免费完整版| 成年人黄色毛片网站| 少妇的丰满在线观看| 看黄色毛片网站| 妹子高潮喷水视频| 国产又色又爽无遮挡免费看| 欧美性猛交黑人性爽| 免费高清在线观看日韩| 午夜免费观看网址| 免费电影在线观看免费观看| 日本免费a在线| 麻豆av在线久日| 成年女人毛片免费观看观看9| 久久 成人 亚洲| 一级a爱视频在线免费观看| 91成人精品电影| 亚洲国产精品sss在线观看| 哪里可以看免费的av片| 亚洲黑人精品在线| 欧美日韩亚洲国产一区二区在线观看| 国产成人精品久久二区二区免费| 亚洲精品国产区一区二| 在线免费观看的www视频| 久久精品国产亚洲av高清一级| 亚洲精华国产精华精| 色精品久久人妻99蜜桃| 美女高潮喷水抽搐中文字幕| 黄色视频不卡| 性色av乱码一区二区三区2| 日本一区二区免费在线视频| 淫妇啪啪啪对白视频| 日韩国内少妇激情av| 伊人久久大香线蕉亚洲五| www.999成人在线观看| 一本久久中文字幕| 男人操女人黄网站| 日本撒尿小便嘘嘘汇集6| 色尼玛亚洲综合影院| 啦啦啦韩国在线观看视频| 免费在线观看视频国产中文字幕亚洲| 成人国产综合亚洲| 亚洲国产毛片av蜜桃av| 国产精品爽爽va在线观看网站 | 欧美日韩亚洲综合一区二区三区_| 欧美日韩福利视频一区二区| 婷婷丁香在线五月| 久久久久久久久免费视频了| 天天一区二区日本电影三级| 老汉色∧v一级毛片| 久久精品夜夜夜夜夜久久蜜豆 | 国产伦一二天堂av在线观看| 亚洲熟女毛片儿| 国产精品免费一区二区三区在线| 国产成人精品无人区| 夜夜夜夜夜久久久久| 脱女人内裤的视频| 久久国产精品影院| 亚洲第一青青草原| 国产精品二区激情视频| 国产精品爽爽va在线观看网站 | 欧美在线一区亚洲| 免费看日本二区| 99国产精品99久久久久| 国产精品免费视频内射| 国产一区在线观看成人免费| 91成年电影在线观看| 精品欧美国产一区二区三| 黄片播放在线免费| 岛国在线观看网站| 亚洲人成网站在线播放欧美日韩| 一区二区三区激情视频| 欧美av亚洲av综合av国产av| 两个人免费观看高清视频| 一卡2卡三卡四卡精品乱码亚洲| 激情在线观看视频在线高清| 亚洲久久久国产精品| 高清毛片免费观看视频网站| 欧美大码av| 亚洲五月色婷婷综合| 亚洲第一青青草原| 欧美乱妇无乱码| 久久久国产成人精品二区| 亚洲国产精品久久男人天堂| 搞女人的毛片| 天天躁狠狠躁夜夜躁狠狠躁| 国产亚洲精品久久久久5区| 国产激情久久老熟女| 亚洲成人免费电影在线观看| 亚洲av成人av| e午夜精品久久久久久久| 欧美激情久久久久久爽电影| 美女午夜性视频免费| 国产精品免费一区二区三区在线| 男女下面进入的视频免费午夜 | √禁漫天堂资源中文www| 日韩欧美一区二区三区在线观看| 美女午夜性视频免费| 国产精品九九99| 男女床上黄色一级片免费看| 国产av又大| 人妻久久中文字幕网| 人妻丰满熟妇av一区二区三区| 亚洲av中文字字幕乱码综合 | 免费人成视频x8x8入口观看| 国产视频内射| 精品久久久久久久久久免费视频| 精品不卡国产一区二区三区| 自线自在国产av| av有码第一页| 国产日本99.免费观看| 麻豆一二三区av精品| 日韩欧美 国产精品| 18禁国产床啪视频网站| 色综合亚洲欧美另类图片| 久9热在线精品视频| 一本精品99久久精品77| 黑人巨大精品欧美一区二区mp4| 国产亚洲欧美在线一区二区| 日本一区二区免费在线视频| 脱女人内裤的视频| 此物有八面人人有两片| 亚洲国产欧美一区二区综合| 国产精品野战在线观看| 国产精品久久久久久亚洲av鲁大| 日本a在线网址| 精品国产乱子伦一区二区三区| 好男人电影高清在线观看| 久久亚洲真实| 亚洲 欧美一区二区三区| 美女免费视频网站| 久久久久九九精品影院| 国产成人av激情在线播放| 日本成人三级电影网站| av天堂在线播放| 中国美女看黄片| 亚洲av中文字字幕乱码综合 | 色av中文字幕| 欧美三级亚洲精品| 一a级毛片在线观看| 久久精品成人免费网站| 夜夜夜夜夜久久久久| 国产精品电影一区二区三区| 国产又爽黄色视频| 99久久国产精品久久久| 一本综合久久免费| 美女高潮喷水抽搐中文字幕| 国产视频内射| 亚洲精品av麻豆狂野| 老司机在亚洲福利影院| 少妇被粗大的猛进出69影院| 色精品久久人妻99蜜桃| 18禁黄网站禁片午夜丰满| 午夜久久久在线观看| 久久久久国产精品人妻aⅴ院| 亚洲三区欧美一区| 亚洲av成人不卡在线观看播放网| av电影中文网址| 亚洲中文字幕日韩| 18禁国产床啪视频网站| 两性夫妻黄色片| 观看免费一级毛片| 日韩精品青青久久久久久| 亚洲第一青青草原| 90打野战视频偷拍视频| 成人手机av| 91老司机精品| 亚洲国产精品sss在线观看| 最近在线观看免费完整版| 欧美大码av| 久久精品人妻少妇| 国产亚洲欧美在线一区二区| 国产又黄又爽又无遮挡在线| 免费在线观看影片大全网站| 91成年电影在线观看| 婷婷精品国产亚洲av在线| 国产精品九九99| 亚洲人成网站高清观看| 久久狼人影院| 18禁裸乳无遮挡免费网站照片 | 国产精品98久久久久久宅男小说| 欧美性猛交╳xxx乱大交人| 一级毛片精品| 丰满的人妻完整版| 国产成人欧美在线观看| 久久久久久九九精品二区国产 | 别揉我奶头~嗯~啊~动态视频| 日日夜夜操网爽| 嫩草影视91久久| 亚洲七黄色美女视频| 日韩欧美国产在线观看| 欧美亚洲日本最大视频资源| 丝袜美腿诱惑在线| 国产av一区在线观看免费| 天天躁狠狠躁夜夜躁狠狠躁| 亚洲五月婷婷丁香| 白带黄色成豆腐渣| 久久久精品欧美日韩精品| 久久午夜综合久久蜜桃| 精品久久久久久,| 亚洲国产欧美一区二区综合| 色婷婷久久久亚洲欧美| 91成年电影在线观看| 精品国内亚洲2022精品成人| 日本精品一区二区三区蜜桃| 午夜激情福利司机影院| 精品不卡国产一区二区三区| av在线天堂中文字幕| 精品国产乱子伦一区二区三区| 国产蜜桃级精品一区二区三区| 99久久国产精品久久久| 中文亚洲av片在线观看爽| 丝袜美腿诱惑在线| 国产又黄又爽又无遮挡在线| 男人舔奶头视频| 听说在线观看完整版免费高清| 午夜成年电影在线免费观看| 一区二区日韩欧美中文字幕| 亚洲avbb在线观看| 亚洲国产毛片av蜜桃av| 久久午夜亚洲精品久久| 国产熟女xx| 精品国产一区二区三区四区第35| 热99re8久久精品国产| 欧美成人午夜精品| 亚洲黑人精品在线| 两性午夜刺激爽爽歪歪视频在线观看 | 久久国产精品人妻蜜桃| 婷婷六月久久综合丁香| 99久久精品国产亚洲精品| 久久香蕉精品热| 欧美性猛交╳xxx乱大交人| 中文字幕av电影在线播放| 成人三级黄色视频| 亚洲成人久久爱视频| 国产亚洲精品av在线| 国产三级黄色录像| 淫妇啪啪啪对白视频| 国产精品,欧美在线| 国产又黄又爽又无遮挡在线| 久久国产精品男人的天堂亚洲| 欧美色欧美亚洲另类二区| 悠悠久久av| 嫁个100分男人电影在线观看| 亚洲 欧美 日韩 在线 免费| 丝袜美腿诱惑在线| 中文在线观看免费www的网站 | 色尼玛亚洲综合影院| 草草在线视频免费看| 中文字幕最新亚洲高清| 国产91精品成人一区二区三区| 天堂影院成人在线观看| 99re在线观看精品视频| 在线av久久热| 亚洲精品在线观看二区| 操出白浆在线播放| 欧美一区二区精品小视频在线| 久99久视频精品免费| 国产精品综合久久久久久久免费| 后天国语完整版免费观看| 身体一侧抽搐| 18禁黄网站禁片午夜丰满| 女人高潮潮喷娇喘18禁视频| 精品无人区乱码1区二区| 村上凉子中文字幕在线| 国产亚洲欧美精品永久| 久久国产精品男人的天堂亚洲| 法律面前人人平等表现在哪些方面| 草草在线视频免费看| 久久人妻av系列| 免费av毛片视频| 伦理电影免费视频| 久久天堂一区二区三区四区| 亚洲男人的天堂狠狠| 在线观看一区二区三区| 日韩欧美国产在线观看| 波多野结衣高清作品| 看片在线看免费视频| а√天堂www在线а√下载| 亚洲一区中文字幕在线| 欧美成人一区二区免费高清观看 | 熟女少妇亚洲综合色aaa.| 婷婷六月久久综合丁香| 制服人妻中文乱码| 亚洲熟妇中文字幕五十中出| 在线国产一区二区在线| 性色av乱码一区二区三区2| 国产黄片美女视频| 在线观看午夜福利视频| 国产99白浆流出| 可以在线观看的亚洲视频| 国产精品综合久久久久久久免费| 最新在线观看一区二区三区| 老熟妇仑乱视频hdxx| 18禁黄网站禁片午夜丰满| 国产亚洲精品第一综合不卡| 在线av久久热| 精品国产美女av久久久久小说| 人妻丰满熟妇av一区二区三区| 叶爱在线成人免费视频播放| 国产精品1区2区在线观看.| 一本久久中文字幕| 日韩免费av在线播放| 在线观看免费午夜福利视频| 日本五十路高清| 操出白浆在线播放| 久久久国产欧美日韩av| 免费看日本二区| av欧美777| 制服人妻中文乱码| 琪琪午夜伦伦电影理论片6080| 国产亚洲精品一区二区www| 欧美成人性av电影在线观看| 国内久久婷婷六月综合欲色啪| 日韩欧美一区视频在线观看| 久久人妻av系列| 亚洲av电影在线进入| 日本熟妇午夜| 又大又爽又粗| 999久久久精品免费观看国产| 黄色片一级片一级黄色片| 亚洲国产毛片av蜜桃av| 国产精品 欧美亚洲| 波多野结衣高清作品| 久久久久久久久中文| 最近最新中文字幕大全免费视频| 国语自产精品视频在线第100页| 国产欧美日韩一区二区精品| 日韩成人在线观看一区二区三区| 色尼玛亚洲综合影院| 啦啦啦 在线观看视频| 久久精品国产综合久久久| 两个人视频免费观看高清| 国产黄片美女视频| 一本一本综合久久| 亚洲av成人不卡在线观看播放网| 亚洲第一欧美日韩一区二区三区| 日韩欧美一区视频在线观看| 日日夜夜操网爽| 51午夜福利影视在线观看| 亚洲av五月六月丁香网| 亚洲aⅴ乱码一区二区在线播放 | 国产99白浆流出| 精品无人区乱码1区二区| 啦啦啦免费观看视频1| 国内久久婷婷六月综合欲色啪| 国产一区二区在线av高清观看| 亚洲午夜理论影院|