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

    包含反應(yīng)階數(shù)變化的分形煤焦顆粒燃燒模型的建立與實(shí)驗(yàn)驗(yàn)證

    2016-03-19 07:31:06劉雨廷何榕
    化工學(xué)報(bào) 2016年1期
    關(guān)鍵詞:實(shí)驗(yàn)驗(yàn)證模擬分形

    劉雨廷,何榕

    (清華大學(xué)熱能工程系,熱科學(xué)與動(dòng)力工程教育部重點(diǎn)實(shí)驗(yàn)室,北京 100084)

    ?

    包含反應(yīng)階數(shù)變化的分形煤焦顆粒燃燒模型的建立與實(shí)驗(yàn)驗(yàn)證

    劉雨廷,何榕

    (清華大學(xué)熱能工程系,熱科學(xué)與動(dòng)力工程教育部重點(diǎn)實(shí)驗(yàn)室,北京 100084)

    摘要:為了深入研究煤焦燃燒的機(jī)理并提高對煤焦燃燒過程的預(yù)測精度,建立了一個(gè)綜合的煤焦燃燒模型。該模型考慮了煤焦顆??紫秲?nèi)二次反應(yīng)與擴(kuò)散的耦合作用、煤焦燃燒反應(yīng)階數(shù)的變化和反應(yīng)過程中CO/CO2比例等問題。使用熱天平(TGA)對11種煤焦的燃燒特性進(jìn)行分析,測得各種煤焦的表觀活化能與指前因子,以確定模型中的待定參數(shù)。在沉降爐(DTF)中對這11種煤焦做燃燒實(shí)驗(yàn),用TGA基于灰分守恒測得DTF爐管出口處的煤焦樣品的轉(zhuǎn)化率。運(yùn)用建立的模型模擬這些煤焦的燃燒過程,預(yù)測的轉(zhuǎn)化率與實(shí)驗(yàn)結(jié)果有較好的吻合度,相比傳統(tǒng)的本征動(dòng)力學(xué)模型,該模型預(yù)測的精度有了較大提高,證明了該模型能適用于從褐煤到無煙煤的較廣煤焦范圍。研究還發(fā)現(xiàn),煤焦燃燒的表觀反應(yīng)階數(shù)在燃燒過程中不斷減小并最終趨于穩(wěn)定。

    關(guān)鍵詞:煤焦燃燒;模型;模擬;沉降爐;實(shí)驗(yàn)驗(yàn)證;表觀反應(yīng)階數(shù);分形

    2015-06-15收到初稿,2015-07-07收到修改稿。

    聯(lián)系人:何榕。第一作者:劉雨廷(1988—),男,博士研究生。

    Received date: 2015-06-15.

    引 言

    在我國,煤是一種儲量相對豐富、價(jià)格相對便宜的化石燃料。隨著石油資源的日益枯竭,煤的價(jià)值與地位越發(fā)突顯,如何高效地利用煤炭資源成為我國經(jīng)濟(jì)發(fā)展不可回避的問題。目前,對煤最主要的利用方式是直接燃燒[1]以獲得熱能,以火力發(fā)電為代表。

    煤的燃燒大致可分為熱解與煤焦燃燒兩個(gè)過程。已有的研究[2-5]表明,煤燃燒發(fā)熱量的80%以上來源于煤焦燃燒,且煤焦燃燒所需時(shí)間比熱解要長幾個(gè)量級。因此,煤的燃燒特性主要由煤焦的燃燒特性決定。所以,建立精確的煤焦燃燒模型,對于煤的高效利用,以及燃煤設(shè)備的設(shè)計(jì)與運(yùn)行,都有著極為重要的意義。

    煤焦是典型的多孔介質(zhì),其孔隙結(jié)構(gòu)極為復(fù)雜,而且具有分形特性[6-8]。煤焦的燃燒反應(yīng)不只在煤焦的外表面進(jìn)行,而更多的是通過氣體擴(kuò)散進(jìn)入孔隙內(nèi)部,在孔隙的內(nèi)表面進(jìn)行,煤焦的孔隙結(jié)構(gòu)對煤焦的燃燒特性有著重要的影響[9-15]。因此,精確的煤焦燃燒模型必須考慮煤焦孔隙結(jié)構(gòu)對煤焦燃燒反應(yīng)的影響。

    研究者們已針對煤焦燃燒模型進(jìn)行了大量的研究[1-3,16-19],建立了許多不同的煤焦燃燒模型。其中包含了煤焦孔隙結(jié)構(gòu)影響的模型之中,經(jīng)典的本征動(dòng)力學(xué)模型[3,17-18]應(yīng)用最為廣泛,而He等[2]建立的具有分形孔隙特性的煤焦燃燒模型則是通過理論推導(dǎo)得到的。該兩種模型都認(rèn)為不同煤焦的本征活化能相同,而煤焦燃燒反應(yīng)速率由于煤焦孔隙結(jié)構(gòu)的不同而存在很大的差異。前者[3,17-18]通過引入彎曲因子表征煤焦孔隙結(jié)構(gòu)的復(fù)雜性,彎曲因子越大,煤焦孔隙結(jié)構(gòu)越復(fù)雜。但彎曲因子并無嚴(yán)格的物理定義,且無法通過獨(dú)立的實(shí)驗(yàn)測量或計(jì)算得到。而后者[2]則采用孔隙率、比表面積、分形維數(shù)[8]3個(gè)參數(shù)表征煤焦孔隙結(jié)構(gòu)的復(fù)雜性,并與顆粒粒徑這一參數(shù)共同表征煤焦的結(jié)構(gòu)特性。與前者相比,后者的優(yōu)點(diǎn)在于其計(jì)算精度更高,能更好地反映煤焦的結(jié)構(gòu)參數(shù)對煤焦燃燒特性的影響,且孔隙率、比表面積、分形維數(shù)定義明確,易于獨(dú)立的實(shí)驗(yàn)測量,用于表征煤焦孔隙的復(fù)雜程度較之彎曲因子更為合理。但是,無論本征動(dòng)力學(xué)模型[3,17-18]還是He等[2]建立的煤焦燃燒模型,都假設(shè)煤焦燃燒反應(yīng)的階數(shù)恒定為1,這與實(shí)際情況不符[20-21],難免會給計(jì)算帶來誤差。而且He等[2]所建立的模型中的待定參數(shù)是通過試算法得到的,方法不太嚴(yán)密,而且用于對比的煤焦燃燒實(shí)驗(yàn)數(shù)目也很有限。

    本文在He等[2]所建立的煤焦燃燒模型的基礎(chǔ)上,引入表觀反應(yīng)階數(shù)的變化規(guī)律[21]和燃燒產(chǎn)物中CO/CO2的比例的定量描述[21-22],并用更加嚴(yán)格的方法確定其模型中的待定參數(shù),從而建立了一個(gè)可靠的煤焦燃燒模型。該模型考慮了煤焦孔隙內(nèi)復(fù)雜的擴(kuò)散與反應(yīng)耦合的過程,不僅能預(yù)測瞬時(shí)的煤焦燃燒速率,可以方便地嵌入到一些商用計(jì)算流體力學(xué)軟件(如Fluent)中,還能獲得煤焦燃燒過程中更多的細(xì)節(jié),便于對煤焦燃燒機(jī)理進(jìn)行研究。使用熱天平(TGA)分別對11種煤焦的燃燒特性進(jìn)行研究[23],并測得它們的表觀活化能與指前因子,用來確定模型中的待定參數(shù)。在沉降爐(DTF)中對這11種煤焦進(jìn)行燃燒實(shí)驗(yàn)[23],用TGA測得這些煤焦在DTF爐管出口處的轉(zhuǎn)化率。運(yùn)用本文中所建立的模型,從最基本的流動(dòng)傳熱傳質(zhì)方程出發(fā),對這些煤焦的燃燒過程進(jìn)行模擬并與實(shí)驗(yàn)結(jié)果對比,以驗(yàn)證模型的可靠性。計(jì)算的結(jié)果還與采用本征動(dòng)力學(xué)模型[3,17-18]的計(jì)算結(jié)果進(jìn)行比較。最后,對煤焦燃燒過程中一些關(guān)鍵量的變化規(guī)律進(jìn)行分析。

    1 模型的建立

    本文所建立的分形煤焦燃燒模型主要包括3個(gè)部分:煤焦燃燒瞬時(shí)速率表達(dá)式、煤焦燃燒的表觀反應(yīng)階數(shù)和燃燒產(chǎn)物中的CO/CO2比例。這些內(nèi)容均是比較復(fù)雜的問題,但在此之前已經(jīng)對它們做了大量的研究工作。

    1.1 煤焦燃燒瞬時(shí)速率表達(dá)式

    He等[2]通過理論推導(dǎo)的方法建立了煤焦燃燒速率的表達(dá)式。煤焦顆粒的示意圖如圖1所示,rp為煤焦顆粒的半徑。假設(shè)燃燒過程中煤焦內(nèi)部達(dá)到穩(wěn)定狀態(tài),煤焦顆粒結(jié)構(gòu)各向同性且均勻。

    在r處,單位時(shí)間內(nèi)擴(kuò)散的氧氣量為

    圖1 煤焦顆粒示意圖Fig.1 Schematic diagram of a char particle

    其中,De為氧氣的等效擴(kuò)散系數(shù),c為r處氧氣的濃度。

    在r~ r+ d r微元體內(nèi),單位時(shí)間內(nèi)消耗的煤焦可表示為

    其中,ρa(bǔ)是煤焦的表觀密度,Sg是煤焦的比表面積,ki為本征反應(yīng)常數(shù),其表達(dá)式為

    穩(wěn)定狀態(tài)下,微元體內(nèi)所消耗的氧氣量等于凈擴(kuò)散進(jìn)入該微元體內(nèi)的氧氣量,故有

    化簡得

    已知氧氣濃度的邊界條件為

    其中,cs為煤焦外表面氧氣濃度。

    根據(jù)邊界條件式(6),求解微分方程式(5),可得煤焦內(nèi)部氧氣濃度的分布

    因此,煤焦的總反應(yīng)速率可表述為

    煤焦燃燒的表觀反應(yīng)常數(shù)ka可定義為

    其中,m為煤焦顆粒質(zhì)量。將Q的定義代入可得

    現(xiàn)引入模數(shù)f

    式中,z為比例系數(shù),M0為氧氣分子摩爾質(zhì)量,Tp為煤焦顆粒的溫度,(Tp/ M0)1/2項(xiàng)是為體現(xiàn)溫度對擴(kuò)散系數(shù)的影響。則有

    為方便分析與處理,引入新的模數(shù)c,將煤焦孔隙的所有結(jié)構(gòu)參數(shù)整合到一起,其定義式為

    其中,ρi為煤焦的真密度。因此,煤焦燃燒的瞬時(shí)速率表達(dá)式為

    其中

    式中,Ai、Ei和z為待定參數(shù),由實(shí)驗(yàn)確定。

    He等[2]進(jìn)一步研究表明,模數(shù)c的值越大,煤焦孔隙結(jié)構(gòu)對氣體擴(kuò)散的阻力越小,該種煤焦就越好燒。

    1.2 煤焦燃燒的表觀反應(yīng)階數(shù)

    目前,包括本征動(dòng)力學(xué)模型在內(nèi),絕大多數(shù)煤焦燃燒模型都將煤焦燃燒的表觀反應(yīng)階數(shù)取為1或其他常數(shù),以簡化計(jì)算,但這與實(shí)際情況并不相符[20-21],難免會給計(jì)算帶來誤差。也有很多學(xué)者對煤焦燃燒反應(yīng)的階數(shù)進(jìn)行了大量實(shí)驗(yàn)研究,但他們大多只給出了在特定反應(yīng)條件下、對特定煤焦的反應(yīng)階數(shù)值,并沒有給出一個(gè)通用的規(guī)律來[21]。這也是前面大多數(shù)模型將反應(yīng)階數(shù)取為常數(shù)的原因。

    Liu等[21]采用數(shù)值模擬的方法對這一問題進(jìn)行研究,歸納出了一個(gè)能夠描述不同煤焦、在不同條件下的反應(yīng)階數(shù)的表達(dá)式,并與實(shí)驗(yàn)結(jié)果進(jìn)行對比,取得了較好的效果。要介紹Liu等[21]的結(jié)果,先要對一個(gè)分形孔隙模型進(jìn)行介紹。

    Liang等[26]采用隨機(jī)漫步算法生成了一種分形孔隙模型,它可以用來模擬真實(shí)的煤焦顆粒。圖2是這種孔隙模型的示意圖,可以看出它與真實(shí)的煤焦在外形上很相近。He等[8]基于對煤焦的壓汞實(shí)驗(yàn),定義了一種煤焦的分形維數(shù)

    其中,Sa是煤焦的累積表面積,rpore是煤焦孔隙的半徑。分形維數(shù)表征了煤焦孔隙結(jié)構(gòu)的復(fù)雜程度,對孔內(nèi)的氣體擴(kuò)散影響很大,是描述煤焦孔隙結(jié)構(gòu)的一個(gè)重要參數(shù)。根據(jù)壓汞實(shí)驗(yàn)的原理,Liang等[26]也給出了該種數(shù)值煤焦顆粒分形維數(shù)的計(jì)算方法。

    Liu等[21]采用圖2所示的分形孔隙模型模擬煤焦燃燒。在他們的模擬中,化學(xué)反應(yīng)(包括碳與O2、碳與CO2及CO與O2的反應(yīng))是用簡單碰撞理論描述的,煤焦孔內(nèi)的氣體擴(kuò)散則是用分形孔隙內(nèi)的氣體擴(kuò)散定律[24]來描述。這里并不需要對煤焦的整個(gè)燃燒過程進(jìn)行模擬,而只需要得到在某一個(gè)具體條件下煤焦總的反應(yīng)速率,也就是說,只需要計(jì)算出一個(gè)極短時(shí)間內(nèi)整個(gè)煤焦顆粒中消耗的碳分子總數(shù)。假設(shè)兩個(gè)反應(yīng)溫度相同、濃度相差一個(gè)小量的工況的反應(yīng)階數(shù)近似相等,通過模擬得到這兩個(gè)工況下的反應(yīng)速率值,便可以很容易地求得該條件下的反應(yīng)階數(shù)了。一共計(jì)算了600組工況,得到了300組有效的數(shù)據(jù)。通過對這些數(shù)據(jù)的歸納分析,得到了煤焦燃燒表觀反應(yīng)階數(shù)的表達(dá)式

    其中,χ是模數(shù),其定義與式(16)中的χ相同,代表了孔隙結(jié)構(gòu)對表觀反應(yīng)階數(shù)的影響;Tp是煤焦顆粒的溫度,K;PO2是顆粒表面O2的分壓,atm。

    Liu等[21]進(jìn)一步證明了模數(shù)χ是煤焦顆??紫稊U(kuò)散阻力的度量,χ越大,對應(yīng)的孔隙擴(kuò)散阻力越小。他們還給出了χ的臨界值為1.3171×10-4,當(dāng)煤焦顆粒的χ值大于此臨界值時(shí),孔隙擴(kuò)散對煤焦燃燒速率的影響可以忽略不計(jì)。

    1.3 煤焦燃燒產(chǎn)物中的CO/CO2比例

    燃燒產(chǎn)物中的CO/CO2比例是影響煤焦燃燒的重要因素之一。同樣質(zhì)量的煤焦氧化生成CO2所釋放的熱量是生成CO所釋放熱量的3.5倍[22-23],由于生成不同產(chǎn)物所釋放的熱量存在巨大差異,會給煤焦顆粒溫度的計(jì)算帶來上百度的誤差。

    He等[22]基于煤焦燃燒活化能高斯分布假設(shè),并結(jié)合氣體分子運(yùn)動(dòng)論的相關(guān)內(nèi)容,對O2分子與碳分子碰撞結(jié)果進(jìn)行了分類,使用概率論的相關(guān)方法計(jì)算出生成各種產(chǎn)物及不發(fā)生反應(yīng)的概率,進(jìn)而得到碳與O2反應(yīng)的一次CO/CO2比例

    其中,x1(ε)和x2(ε)是中間函數(shù),具體形式見文獻(xiàn)[22],f(ε)是O2分子平動(dòng)能的麥克斯維分布函數(shù)。之所以稱它是一次CO/CO2比例,是因?yàn)樗⒉话航诡w??紫秲?nèi)的氣體擴(kuò)散與二次反應(yīng)的作用,這用實(shí)驗(yàn)方法是很難得到的。

    Liu等[21]將式(19)運(yùn)用到煤焦燃燒的模擬中,通過使用圖2中的數(shù)值煤焦顆粒模型,引入煤焦顆??紫秲?nèi)的擴(kuò)散與化學(xué)反應(yīng)等因素,得到煤焦顆粒燃燒的表觀CO/CO2比例

    其中,Tp是煤焦顆粒的溫度,K;PO2是顆粒表面O2的分壓,atm。這與Tognotti等[27]的實(shí)驗(yàn)結(jié)果非常接近。

    2 實(shí) 驗(yàn)

    使用熱天平(TGA)分別對11種煤焦的燃燒特性進(jìn)行研究[23],測得它們的表觀活化能與指前因子,用來確定模型中的待定參數(shù)。在沉降爐(DTF)中對這11種煤焦進(jìn)行燃燒實(shí)驗(yàn)[23],用TGA測得這些煤焦在DTF爐管出口處的轉(zhuǎn)化率。運(yùn)用本文中所建立的模型對這些煤焦的燃燒過程進(jìn)行模擬并與實(shí)驗(yàn)結(jié)果對比,以驗(yàn)證模型的可靠性。

    2.1 實(shí)驗(yàn)材料

    表1和表2分別為這11種煤焦對應(yīng)的原煤的工業(yè)分析和元素分析結(jié)果。這些煤中既有高揮發(fā)分的煙煤,也有不太好燒的無煙煤與褐煤。將這些原煤在沉降爐中、1200℃高溫、N2氣氛下進(jìn)行熱解,以制得后面實(shí)驗(yàn)用的煤焦樣品。熱解時(shí)氣體的流量為0.5 m3·h-1,氣體平均流速為0.33 m·s-1,給粉速率為1 g·h-1。所得煤焦樣品中殘余揮發(fā)分含量很低,對后面煤焦燃燒影響很小。

    要對這些煤焦的燃燒過程進(jìn)行模擬,需要測量這些煤焦樣品的平均粒徑與孔隙參數(shù)。對于平均粒徑的測量,使用的是Mastersizer 2000激光粒徑分析儀,它的量程是0.02~2000 μm,測量誤差在1%以內(nèi),重復(fù)測量的誤差小于0.5%。對于其他孔隙參數(shù)的測量,使用的是Autoscan 33壓汞儀,其操作壓力范圍是0~232 MPa,能測量的最小孔徑是3.23 nm,非線性度為滿量程的±0.05%。測得的這11種煤焦的平均粒徑與孔隙參數(shù)列于表3中,其中分形維數(shù)的定義為式(17)。

    表1 原煤的工業(yè)分析結(jié)果(空干基)Table 1 Proximate analyses of parent coals (dry air basis)

    表2 原煤的元素分析結(jié)果(干燥基)Table 2 Ultimate analyses of parent coals (dry basis)

    HV  71.39  4.78  1.56  9.81  0.46

    表3 煤焦樣品的平均粒徑與孔隙參數(shù)Table 3 Average particle diameters and pore structure parameters of chars

    2.2 熱重分析(TGA)實(shí)驗(yàn)

    將煤焦樣品分散均勻地置于熱天平上,樣品上面所通氣體為空氣,熱重采用程序升溫,溫度與時(shí)間呈線性關(guān)系。由于實(shí)驗(yàn)過程中加熱速率很慢,可以保證煤焦顆粒溫度與氣流溫度近似相等。各種煤焦實(shí)驗(yàn)時(shí)溫度與時(shí)間的關(guān)系如表4所示,其中,t為時(shí)間,s;T為氣體和煤焦顆粒的溫度,K。

    表4 TGA實(shí)驗(yàn)升溫曲線函數(shù)Table 4 Heating curve function in TGA experiments

    由于煤焦顆粒燃燒過程中粒徑和孔隙結(jié)構(gòu)是不斷變化的[3],只有得到與煤焦燃燒前的粒徑和孔隙參數(shù)(表3中所列)相對應(yīng)的燃燒動(dòng)力學(xué)參數(shù),才能確定本文所建立模型中的待定參數(shù)(Ai,Ei和)的值。于是,在處理TGA數(shù)據(jù)時(shí),僅取煤焦轉(zhuǎn)化率在0~20%之間的數(shù)據(jù)進(jìn)行分析。在該段數(shù)據(jù)中,煤焦的轉(zhuǎn)化率最多只有20%,用文獻(xiàn)[3,28]中的經(jīng)驗(yàn)公式估算得到此時(shí)的粒徑與孔隙參數(shù)變化相對初始值不超過5%,因此,顆粒的粒徑與孔隙參數(shù)變化均不明顯,可以近似將用此段數(shù)據(jù)得到的表觀動(dòng)力學(xué)參數(shù)與表3中的粒徑和孔隙參數(shù)對應(yīng)起來。

    具體求解煤焦燃燒表觀動(dòng)力學(xué)參數(shù)(表觀活化能Ea和指前因子Aa)的方法見文獻(xiàn)[23,29]。這里只將通過最后迭代求解的結(jié)果列于表5中。

    表5 TGA實(shí)驗(yàn)得到的煤焦燃燒表觀動(dòng)力學(xué)參數(shù)Table 5 Apparent kinetic parameters of char combustion obtained in TGA experiments

    2.3 沉降爐(DTF)燃燒實(shí)驗(yàn)

    實(shí)驗(yàn)所用沉降爐系統(tǒng)的結(jié)構(gòu)如圖3所示。沉降爐使用電加熱,爐管的內(nèi)徑為50 mm,爐管高為750 mm,煤焦顆粒注入爐管的位置距離爐管進(jìn)口200 mm,煤焦取樣點(diǎn)設(shè)置在爐管出口處。

    圖3 沉降爐結(jié)構(gòu)示意圖Fig.3 Schematic diagram of drop tube furnace

    各種煤焦的燃燒實(shí)驗(yàn)條件是相同的,即沉降爐壁溫為1000℃,爐內(nèi)所通氣體為空氣,流量為0.5 m3·h-1,氣體平均流速為0.33 m·s-1,煤焦顆粒的給粉速率為1 g·h-1,煤焦顆粒在沉降爐內(nèi)的平均停留時(shí)間約為1.67 s。煤焦的給粉速率足夠小,以保證燃燒過程中煤焦顆粒之間互不影響。

    沉降爐爐管出口處所取的煤焦樣品的轉(zhuǎn)化率由式(21)計(jì)算

    其中,X為煤焦的轉(zhuǎn)化率,ash0為反應(yīng)前煤焦顆粒的灰分含量,ash為燃燒后煤焦顆粒的灰分含量。ash0和ash的值均由TGA實(shí)驗(yàn)測得,結(jié)果列于表6中。

    表6 煤焦在沉降爐中燃燒的最終轉(zhuǎn)化率Table 6 Final conversions of chars in DTF experiments

    3 分析與討論

    3.1 模型參數(shù)的確定

    1.1節(jié)中經(jīng)理論推導(dǎo)得到了煤焦顆粒燃燒的表觀反應(yīng)常數(shù)ka的表達(dá)式(16),其中包括待定的模型參數(shù)Ai,Ei和z。同時(shí),表觀反應(yīng)常數(shù)ka還可以由式(22)計(jì)算

    其中,Ea為表觀活化能,Aa為對應(yīng)的指前因子。這樣式(16)與式(22)對相同煤焦、在相同溫度下便可建立起一個(gè)等式關(guān)系

    表5中給出了用TGA測得的11種煤焦的表觀活化能Ea和指前因子Aa的值,表3又給出了式(23)等號左邊項(xiàng)計(jì)算所需的煤焦粒徑值與孔隙參數(shù)值,則對11種煤焦分取800 K、1000 K、1200 K、1400K、1600 K和1800 K這6個(gè)溫度,便可得到66個(gè)等式,運(yùn)用最小二乘法擬合,讓這66個(gè)等式的誤差最小,這時(shí)的Ai,Ei和z值即為所求。計(jì)算得到它們的值分別為

    為了更直觀地顯示擬合精度,現(xiàn)將溫度為800 K和1000 K時(shí)的擬合結(jié)果與實(shí)驗(yàn)結(jié)果的對比分別作于圖4和圖5中。圖中實(shí)線為由式(16)結(jié)合式(24)中的參數(shù)值所作,為擬合結(jié)果;圖中的點(diǎn)則為式(22)結(jié)合表5中的表觀動(dòng)力學(xué)參數(shù)計(jì)算得到,為實(shí)驗(yàn)結(jié)果,點(diǎn)與橫坐標(biāo)c的對應(yīng)關(guān)系可以由對應(yīng)煤焦顆粒的孔隙參數(shù)與粒徑計(jì)算得出??梢钥闯鰯M合的效果還是不錯(cuò)的。

    這樣,式(16)、式(18)、式(20)以及式(24)中的參數(shù)值一起,構(gòu)成了本文所構(gòu)建的分形煤焦燃燒模型。該模型中,模數(shù)c是所有孔隙參數(shù)與粒徑的組合,表征了孔隙擴(kuò)散對煤焦燃燒影響[2],c的值越大,孔隙的擴(kuò)散阻力就越小,對煤焦燃燒速率的影響也就越小。當(dāng)c的值大于1.3171×10-4時(shí),可以認(rèn)為煤焦孔隙擴(kuò)散的影響可以忽略不計(jì)[21]。

    圖4 800 K時(shí)擬合結(jié)果與實(shí)驗(yàn)結(jié)果的比較Fig.4 Comparison of data correlation results with experimental results at 800 K temperature

    圖5 1000 K時(shí)擬合結(jié)果與實(shí)驗(yàn)結(jié)果的比較Fig.5 Comparison of data correlation results with experimental results at 1000 K temperature

    3.2 沉降爐燃燒實(shí)驗(yàn)的模擬方法

    運(yùn)用本文所建立的煤焦燃燒模型對2.3節(jié)中介紹的沉降爐燃燒實(shí)驗(yàn)進(jìn)行模擬,并將預(yù)測的煤焦最終轉(zhuǎn)化率與表6中的實(shí)驗(yàn)結(jié)果進(jìn)行比較,以驗(yàn)證模型的可靠性。同時(shí),再采用本征動(dòng)力學(xué)模型[3,17-18]對這些實(shí)驗(yàn)進(jìn)行模擬,也與本文模型模擬的結(jié)果進(jìn)行比較。

    沉降爐中的流動(dòng)是一個(gè)典型的氣固兩相流,可以用Euler-Lagrange方法對此進(jìn)行模擬。作為連續(xù)相的空氣用基本的Navier-Stokes方程進(jìn)行描述,而作為離散相的煤焦顆粒通過跟蹤已確定的流場中顆粒的運(yùn)動(dòng)軌跡來描述。在計(jì)算當(dāng)中,煤焦顆粒與空氣流場不斷交換質(zhì)量、動(dòng)量與能量。由于煤焦顆粒的尺寸較小,顆粒之間的相互作用以及顆粒對流場的影響被忽略。

    計(jì)算時(shí)使用的三維網(wǎng)格如圖6所示,一共有22萬個(gè)節(jié)點(diǎn)??諝鈴臍怏w入口(gas inlet)進(jìn)入爐管,煤焦顆粒從離進(jìn)口200 mm的顆粒源(particle source)注入,忽略了顆粒注射器對流場的影響。模擬中使用的具體方程見文獻(xiàn)[2],這里只對煤焦燃燒模型嵌入部分進(jìn)行介紹。

    圖6 沉降爐的三維網(wǎng)格Fig.6 3D-grid of drop tube furnace

    煤焦顆粒燃燒速率方程為

    其中,ka為煤焦燃燒的表觀反應(yīng)常數(shù),由1.1節(jié)中建立的燃燒速率模型[式(16)]來描述,m為煤焦顆粒的質(zhì)量,Sg為煤焦顆粒的比表面積,CO2為煤焦顆粒表面O2的濃度,n為煤焦燃燒的表觀反應(yīng)階數(shù),由1.2節(jié)中式(18)描述。燃燒產(chǎn)物中CO與CO2的比例由1.3節(jié)中的式(20)來描述。這樣,就把前面建立的分形煤焦燃燒模型嵌入到了對沉降爐實(shí)驗(yàn)的模擬中。

    3.3 模型的驗(yàn)證及計(jì)算結(jié)果分析

    運(yùn)用3.2節(jié)中介紹的方法,分別采用本文建立的分形煤焦燃燒模型和傳統(tǒng)的本征動(dòng)力學(xué)模型對2.3節(jié)中11種煤焦在沉降爐中的燃燒進(jìn)行了模擬。將實(shí)驗(yàn)測量與模型預(yù)測的煤焦最終轉(zhuǎn)化率列于表7中。為了對結(jié)果進(jìn)行更直觀的分析,將表7中的數(shù)據(jù)作于圖7中。

    表7 煤焦最終轉(zhuǎn)化率的實(shí)驗(yàn)結(jié)果和模擬結(jié)果Table 7 Measured and predicted final conversions of all chars

    圖7 沉降爐實(shí)驗(yàn)結(jié)果與模擬結(jié)果的比較Fig.7 Comparison of DTF experimental results with model predicted results

    可以看到,兩種模型計(jì)算結(jié)果的大部分?jǐn)?shù)據(jù)點(diǎn)都接近直線y= x,說明對實(shí)驗(yàn)所用的大部分煤焦,這兩個(gè)燃燒速率模型都有較好精度。仔細(xì)比較可以發(fā)現(xiàn),對于同一種煤焦,本文構(gòu)建模型的數(shù)據(jù)點(diǎn)要比本征動(dòng)力學(xué)模型的數(shù)據(jù)點(diǎn)更靠近直線y= x,這也就說明本文所構(gòu)建模型的精度要更高。

    以Luoyang、Yongcheng、Hebi和Newlands 4種煤焦為例,圖8和圖9分別給出了用本文構(gòu)建的燃燒速率模型模擬得到的煤焦顆粒溫度Tp和轉(zhuǎn)化率X隨時(shí)間變化的曲線。

    圖8 煤焦燃燒過程中顆粒溫度的變化Fig.8 Variation of char particle temperature during combustion

    圖9 燃燒過程中煤焦顆粒的轉(zhuǎn)化率Fig.9 Char particle conversions during combustion

    由圖8可知,煤焦顆粒初始為室溫(27℃)狀態(tài),由于顆粒尺寸很小,沉降爐壁溫與周圍氣流溫度很高,顆粒升溫很快,隨著燃燒的開始,又有大量的熱量釋放,因而顆粒最高溫度會高出壁溫許多,最后隨著燃燒速率的減慢,顆粒的溫度開始緩緩下降。

    由圖9可知,煤焦轉(zhuǎn)化率隨時(shí)間不斷增大,剛開始增長較快,隨后逐漸減慢。不同煤焦的轉(zhuǎn)化率相差很大,這與煤焦本身特性有關(guān),反映在模型中便是煤焦的孔隙結(jié)構(gòu)是否有利于O2向顆粒內(nèi)部擴(kuò)散,這可以由模數(shù)c的值來度量。

    圖10 煤焦燃燒過程中表觀反應(yīng)階數(shù)的變化Fig.10 Variation of apparent reaction order during combustion

    圖10中給出了Luoyang煤焦燃燒過程中表觀反應(yīng)階數(shù)n的變化,其他煤焦情況與此類似。煤焦燃燒之初,顆粒溫度較低,燃燒反應(yīng)處于RegimeⅠ,屬于動(dòng)力學(xué)控制區(qū),燃燒速率與煤焦顆粒外表面O2濃度呈正比,因此n的值近似等于1。隨著反應(yīng)溫度的升高,燃燒反應(yīng)漸漸向Regime Ⅱ或Regime Ⅲ過渡,孔隙擴(kuò)散對煤焦燃燒的影響越來越大,因此,當(dāng)煤焦顆粒外表面O2濃度發(fā)生變化時(shí),煤焦顆??紫秲?nèi)表面的O2濃度變化總會有延遲,這一作用體現(xiàn)到反應(yīng)階數(shù)上便是使其減小。因此,表觀反應(yīng)階數(shù)n在煤焦燃燒過程中是逐漸減小,最終趨于穩(wěn)定的。

    對于各孔隙結(jié)構(gòu)參數(shù)(孔隙率、比表面積、分形維數(shù))和粒徑對煤焦燃燒表觀速率常數(shù)ka的影響,He等[2]已在其論文中分析過了,這里不再贅述。只給出最后結(jié)論:表觀速率常數(shù)ka隨著分形維數(shù)、比表面積、顆粒粒徑的增大而減小,而隨著孔隙率的增大而增大。

    4 結(jié) 論

    本文建立了一個(gè)綜合的煤焦燃燒模型,考慮了顆粒孔隙內(nèi)擴(kuò)散與反應(yīng)的耦合作用,煤焦燃燒過程中表觀反應(yīng)階數(shù)變化和燃燒產(chǎn)物中的CO/CO2比例也有可靠的定量描述。表觀速率常數(shù)ka隨著煤焦顆粒的分形維數(shù)、比表面積、粒徑的增大而減小,而隨著孔隙率的增大而增大。

    用熱重分析法(TGA)測得11種煤階跨度較大的煤焦的表觀活化能與指前因子,用以確定模型中的待定參數(shù)。在沉降爐內(nèi)對這11種煤焦進(jìn)行燃燒實(shí)驗(yàn),測得其最終轉(zhuǎn)化率。分別用本征動(dòng)力學(xué)模型及本文所建立的煤焦燃燒模型對這些煤焦的燃燒過程進(jìn)行模擬,本文所建立的模型預(yù)測的轉(zhuǎn)化率與實(shí)驗(yàn)數(shù)據(jù)符合較好,而且比本征動(dòng)力學(xué)模型精度更高,從而證明該模型可以適用于從褐煤到無煙煤的較廣煤焦范圍。

    研究還發(fā)現(xiàn),燃燒過程中煤焦顆粒的最高溫度會比周圍的環(huán)境溫度更高,這是由反應(yīng)放熱所致;表觀反應(yīng)階數(shù)n由1開始逐漸減小,并最終穩(wěn)定下來,這說明反應(yīng)過程中煤焦孔隙擴(kuò)散作用隨著顆粒溫度的升高而不斷增強(qiáng)。

    符 號 說 明

    References

    [1] 堯志輝, 曠戈, 林誠, 等. 單顆粒煤焦燃燒反應(yīng)動(dòng)力學(xué)研究方法[J].化工學(xué)報(bào), 2009, 60 (6): 1442-1451. YAO Z H, KUANG G, LIN C, et al. Methodological study on combustion reaction kinetics of single coal char particle [J]. CIESC Journal, 2009, 60 (6): 1442-1451.

    [2] HE W, LIU Y, HE R, et al. Combustion rate for char with fractal pore characteristics [J]. Combustion Science and Technology, 2013, 185 (11): 1624-1643.

    [3] SMITH I W. The combustion rates of coal chars: a review [J]. Symposium (International) on Combustion, 1982, 19 (1): 1045-1065.

    [4] 王芳, 曾璽, 王永剛, 等. 微型流化床與熱重測定煤焦非等溫氣化反應(yīng)動(dòng)力學(xué)對比 [J].化工學(xué)報(bào), 2015, 66 (5): 1716-1722. WANG F, ZENG X, WANG Y G, et al. Comparation of non-isothermal coal char gasification in micro fluidized bed and thermogravimetric analyzer [J]. CIESC Journal, 2015, 66 (5): 1716-1722.

    [5] 虞君武, 陳永利, 何榕, 等. FD模型應(yīng)用于煤熱解過程的數(shù)值模擬 [J].化工學(xué)報(bào), 2014, 65 (9): 3592-3598. YU J W, CHEN Y L, HE R, et al. Numerical simulation of coal pyrolysis using FD model [J]. CIESC Journal, 2014, 65 (9): 3592-3598.

    [6] AVNIR D, FARIN D, PFEIFER P. Surface geometric irregularity of particulate materials: the fractal approach [J]. Journal of Colloid and Interface Science, 1985, 103 (1): 112-123.

    [7] HE R, SATO J, CHEN C H. Modeling char combustion with fractal pore effects [J]. Combustion Science and Technology, 2002, 174 (4): 19-37.

    [8] HE R, XU X, CHEN C, et al. Evolution of pore fractal dimensions for burning porous chars [J]. Fuel, 1998, 77 (12): 1291-1295.

    [9] EVERSON R C, NEOMAGUS H W J P, KAITANO R. The random pore model with intraparticle diffusion for the description of combustion of char particles derived from mineral and inertinite rich coal [J]. Fuel, 2011, 90 (7): 2347-2352.

    [10] FEI H, HU S, XIANG J, et al. Study on coal chars combustion under O2/CO2atmosphere with fractal random pore model [J]. Fuel, 2011, 90 (2): 441-448.

    [11] 費(fèi)華, 胡松, 向軍, 等. 隨機(jī)孔模型研究煤焦O2/CO2燃燒動(dòng)力學(xué)特征 [J].化工學(xué)報(bào), 2011, 62 (1): 199-205. FEI H, HU S, XIANG J, et al. Kinetics of coal char combustion with random pore model under O2/CO2atmosphere [J]. CIESC Journal, 2011, 62 (1): 199-205.

    [12] FENG B, BHATIA S K. Variation of the pore structure of coal chars during gasification [J]. Carbon, 2003, 41 (3): 507-523.

    [13] 任軼舟, 王亦飛, 朱龍雛, 等. 高溫煤焦氣化反應(yīng)的Langmuir-Hinshelwood動(dòng)力學(xué)模型 [J].化工學(xué)報(bào), 2014, 65 (10): 3906-3915. REN Y Z, WANG Y F, ZHU L C, et al. Langmuir-Hinshelwood kinetic model of high temperature coal char gasification reaction [J]. CIESC Journal, 2014, 65 (10): 3906-3915.

    [14] SADHUKHAN A K, GUPTA P, SAHA R K. Characterization of porous structure of coal char from a single devolatilized coal particle: coal combustion in a fluidized bed [J]. Fuel Processing Technology, 2009, 90 (5): 692-700.

    [15] SU J L, PERLMUTTER D D. Effect of pore structure on char oxidation kinetics [J]. AIChE Journal, 1985, 31 (6): 973-981.

    [16] BARRANCO R, ROJAS A, BARRAZA J, et al. A new char combustion kinetic model (1): Formulation [J]. Fuel, 2009, 88 (12): 2335-2339.

    [17] FORTSCH D, ESSENHIGH R H, SCHNELL U, et al. On the application of the Thiele/Zeldovich analysis to porous carbon combustion [J]. Energy and Fuels, 2003, 17 (4): 901-906.

    [18] MA L, MITCHELL R. Modeling char oxidation behavior under Zone Ⅱ burning conditions at elevated pressures [J]. Combustion and Flame, 2009, 156 (1): 37-50.

    [19] ROJAS A, BARRAZA J, BARRANCO R, et al. A new char combustion kinetic model (2): Empirical validation [J]. Fuel, 2012, 96: 168-175.

    [20] HURT R H, CALO J M. Semi-global intrinsic kinetics for char combustion modeling [J]. Combustion and Flame, 2001, 125 (3): 1138-1149.

    [21] LIU Y, HE R. Variation of apparent reaction order in char combustion and its effect on a fractal char combustion model [J]. Combustion Science and Technology, 2015, 187 (10): 1638-1660.

    [22] HE W, HE R, ITO T, et al. Numerical investigations of CO/CO2ratio in char combustion [J]. Combustion Science and Technology, 2011, 183 (9): 868-882.

    [23] 何威. 煤焦燃燒模型[D]. 北京: 清華大學(xué), 2011. HE W. Study of char combustion model [D]. Beijing: Tsinghua University, 2011.

    [24] CAO L, HE R. Gas diffusion in fractal porous media [J]. Combustion Science and Technology, 2010, 182 (7): 822-841.

    [25] KNUDSEN M. The law of the molecular flow and viscosity of gases moving through tubes [J]. Annals of Physics, 1909, 28 (3): 75-130.

    [26] LIANG Z, HE R, CHEN Q, et al. Fractal generation of char pores through random walk [J]. Combustion Science and Technology, 2007, 179 (3): 637-661.

    [27] TOGNOTTI L, LONGWELL J P, SAROFIM A F. The products of the high temperature oxidation of a single char particle in an electrodynamic balance [J]. Symposium (International) on Combustion, 1991, 23 (1): 1207-1213.

    [28] BHATIA S K, PERLMUTTER D D. A random pore model for fluid-solid reactions(I): Isothermal, kinetic control [J]. AIChE J., 1980, 26 (3): 379-386.

    [29] HE R, SATO J I, CHEN Q, et al. Thermogravimetric analysis of char combustion [J]. Combustion Science and Technology, 2002, 174 (4): 1-18.

    Foundation item: supported by the National Natural Science Foundation of China (51176096).

    Development and experimental validation of fractal char particle combustion model with effect of various reaction order

    LIU Yuting, HE Rong
    (Key Laboratory of Thermal Science and Power Engineering of Ministry of Education, Department of Thermal Engineering, Tsinghua University, Beijing 100084, China)

    Abstract:A fractal char particle combustion model was developed to study the char combustion mechanism and improve the prediction accuracy of the char combustion characteristics. This model considered the coupling effects of the secondary reactions and gas diffusion within the char pores. The apparent reaction order of char combustion changed with the particle temperature, the O2partial pressure at the particle surface and the particle pore structure parameters instead of keeping constant 1. The primary CO/CO2ratio was obtained by theoretical derivation. The combustion characteristics of eleven different chars were analyzed by a thermogravimetric analyzer (TGA) with their apparent activation energies and pre-exponential factors measured which were used to solve the undetermined model parameters. These chars were also combusted in a drop tube furnace (DTF) and the char samples were collected at the outlet of the DTF. The conversions of these char samples were measured by TGA based on the ash conservation method. The combustion processes of these chars were simulated by this model. The predicted char conversions by this model were matched very well with the measured data and the accuracy was significantly improved compared with the intrinsic model, proving that this model had good accuracy over a wide range of char particles from lignite to anthracite. The study also showed that the apparentbook=340,ebook=348reaction order of char combustion first decreased and then gradually stabilized in the char combustion process.

    Key words:char combustion; model; simulation; drop tube furnace; experimental validation; apparent reaction order; fractals

    Corresponding author:Prof. HE Rong, rhe@mail.tsinghua.edu.cn

    基金項(xiàng)目:國家自然科學(xué)基金項(xiàng)目(51176096)。

    中圖分類號:TQ 534

    文獻(xiàn)標(biāo)志碼:A

    文章編號:0438—1157(2016)01—0339—10

    DOI:10.11949/j.issn.0438-1157.20150934

    猜你喜歡
    實(shí)驗(yàn)驗(yàn)證模擬分形
    感受分形
    分形之美
    分形空間上廣義凸函數(shù)的新Simpson型不等式及應(yīng)用
    淺析如何培養(yǎng)初中生學(xué)習(xí)物理的興趣
    一個(gè)高分子模擬計(jì)算網(wǎng)格的作業(yè)管理
    工業(yè)機(jī)器人模擬仿真技術(shù)在職業(yè)教育中的應(yīng)用淺析
    淺析柔道運(yùn)動(dòng)員的模擬實(shí)戰(zhàn)訓(xùn)練
    虛擬機(jī)局域網(wǎng)組建技術(shù)應(yīng)用初探
    RoboCup中型組機(jī)器人教練機(jī)系統(tǒng)設(shè)計(jì)
    試觸法電表偏轉(zhuǎn)變化規(guī)律的實(shí)驗(yàn)測試與分析
    成人无遮挡网站| 久久久色成人| 99热这里只有是精品在线观看| 麻豆av噜噜一区二区三区| 九九爱精品视频在线观看| 亚洲图色成人| 观看免费一级毛片| 精品少妇黑人巨大在线播放 | 日本-黄色视频高清免费观看| 99在线视频只有这里精品首页| 国产一区二区在线av高清观看| 一a级毛片在线观看| 亚洲第一区二区三区不卡| a级一级毛片免费在线观看| 蜜桃久久精品国产亚洲av| 男人的好看免费观看在线视频| 国产高清视频在线观看网站| 亚洲国产精品sss在线观看| 欧美最黄视频在线播放免费| 久久6这里有精品| 国产欧美日韩精品亚洲av| 你懂的网址亚洲精品在线观看 | 久久久久久久久久成人| 麻豆一二三区av精品| 久久久久久久亚洲中文字幕| 国产成人a区在线观看| 在线看三级毛片| 欧美色欧美亚洲另类二区| 久久久久久久久久久丰满| 久久久精品94久久精品| 国产人妻一区二区三区在| 精品99又大又爽又粗少妇毛片| 联通29元200g的流量卡| 夜夜看夜夜爽夜夜摸| 国产精品福利在线免费观看| 日本精品一区二区三区蜜桃| 国产一区二区在线av高清观看| 女人十人毛片免费观看3o分钟| 99久久精品国产国产毛片| 欧美xxxx性猛交bbbb| 成年av动漫网址| 国产欧美日韩精品一区二区| 国产精品一区二区三区四区免费观看 | 别揉我奶头~嗯~啊~动态视频| 蜜桃亚洲精品一区二区三区| 看黄色毛片网站| 免费av不卡在线播放| 18禁黄网站禁片免费观看直播| 淫秽高清视频在线观看| 成人性生交大片免费视频hd| 午夜精品一区二区三区免费看| 在线观看66精品国产| .国产精品久久| 精品一区二区三区视频在线观看免费| 99久久精品热视频| 校园人妻丝袜中文字幕| 一个人观看的视频www高清免费观看| 久久午夜亚洲精品久久| av在线老鸭窝| 变态另类成人亚洲欧美熟女| 最新中文字幕久久久久| 丰满乱子伦码专区| 2021天堂中文幕一二区在线观| 久久精品国产自在天天线| 一个人看的www免费观看视频| 国产一区二区亚洲精品在线观看| 色5月婷婷丁香| 久久久欧美国产精品| 夜夜夜夜夜久久久久| 久久精品国产亚洲网站| 国产人妻一区二区三区在| 国产 一区 欧美 日韩| a级毛片a级免费在线| 97超视频在线观看视频| 国产男靠女视频免费网站| 国产精品一区二区性色av| 啦啦啦观看免费观看视频高清| 国产精品一区二区免费欧美| 国产又黄又爽又无遮挡在线| 亚洲欧美日韩高清在线视频| 精品少妇黑人巨大在线播放 | 亚洲欧美成人综合另类久久久 | av天堂在线播放| 日韩精品有码人妻一区| 一区福利在线观看| ponron亚洲| 久久久午夜欧美精品| 22中文网久久字幕| 黄色日韩在线| 国产三级中文精品| 在线观看av片永久免费下载| 国产成人精品久久久久久| 色5月婷婷丁香| 国产精品久久电影中文字幕| 给我免费播放毛片高清在线观看| 色综合站精品国产| 男女边吃奶边做爰视频| 有码 亚洲区| 免费av毛片视频| 男女下面进入的视频免费午夜| 91麻豆精品激情在线观看国产| 深夜精品福利| 中国国产av一级| videossex国产| 99久国产av精品国产电影| 天天躁日日操中文字幕| 国产真实乱freesex| 亚洲熟妇熟女久久| 欧美日本亚洲视频在线播放| 嫩草影院精品99| 寂寞人妻少妇视频99o| 欧美在线一区亚洲| 国产精品乱码一区二三区的特点| 久久综合国产亚洲精品| 99久久久亚洲精品蜜臀av| 青春草视频在线免费观看| 人妻少妇偷人精品九色| 中国国产av一级| 三级男女做爰猛烈吃奶摸视频| 欧美成人免费av一区二区三区| 午夜亚洲福利在线播放| 身体一侧抽搐| 淫妇啪啪啪对白视频| 国产成人aa在线观看| 午夜影院日韩av| 精品日产1卡2卡| 精品国产三级普通话版| 国产视频内射| 国产精品福利在线免费观看| 国产亚洲av嫩草精品影院| 特级一级黄色大片| 久久久精品94久久精品| 国产精品久久久久久久久免| 日本黄大片高清| 狠狠狠狠99中文字幕| 日日摸夜夜添夜夜添小说| 少妇丰满av| 1024手机看黄色片| 国产高清不卡午夜福利| 内地一区二区视频在线| 久久久精品欧美日韩精品| 欧美日韩综合久久久久久| 久久久久久久亚洲中文字幕| 国产伦精品一区二区三区视频9| 中出人妻视频一区二区| 久久精品影院6| 一进一出抽搐gif免费好疼| 亚洲在线观看片| 午夜老司机福利剧场| 国产精品无大码| 尾随美女入室| 久久久色成人| 国产黄a三级三级三级人| 国产极品精品免费视频能看的| 亚洲国产精品国产精品| 俄罗斯特黄特色一大片| 国产黄片美女视频| 久久久a久久爽久久v久久| 99riav亚洲国产免费| 亚洲国产欧洲综合997久久,| 黄色日韩在线| 九九爱精品视频在线观看| 国产成人影院久久av| 国产欧美日韩精品亚洲av| 此物有八面人人有两片| 国产一区二区在线观看日韩| 91麻豆精品激情在线观看国产| 日本与韩国留学比较| 高清午夜精品一区二区三区 | 国产一区二区激情短视频| 国产极品精品免费视频能看的| 亚洲三级黄色毛片| 乱系列少妇在线播放| 亚洲18禁久久av| 在线观看免费视频日本深夜| 亚洲av不卡在线观看| 国产精品福利在线免费观看| 亚洲丝袜综合中文字幕| 日日撸夜夜添| 麻豆一二三区av精品| 国产av不卡久久| 国产极品精品免费视频能看的| 亚洲图色成人| 国产精品久久久久久av不卡| 我的老师免费观看完整版| 亚洲图色成人| 能在线免费观看的黄片| 亚洲aⅴ乱码一区二区在线播放| 久久综合国产亚洲精品| 亚洲最大成人av| 狂野欧美激情性xxxx在线观看| 久久人人爽人人爽人人片va| 亚洲av.av天堂| 亚洲综合色惰| 91麻豆精品激情在线观看国产| 久久久成人免费电影| www.色视频.com| 亚洲激情五月婷婷啪啪| 伦理电影大哥的女人| av专区在线播放| 色尼玛亚洲综合影院| 18禁在线无遮挡免费观看视频 | 国产精品久久视频播放| 有码 亚洲区| 国产伦精品一区二区三区视频9| 成人av一区二区三区在线看| 免费av不卡在线播放| 精品人妻一区二区三区麻豆 | 精品人妻视频免费看| 国产精品亚洲一级av第二区| 联通29元200g的流量卡| 午夜福利视频1000在线观看| 乱人视频在线观看| 男女下面进入的视频免费午夜| 精品福利观看| 老司机福利观看| 国产精品嫩草影院av在线观看| 久久久久久久久久黄片| 91午夜精品亚洲一区二区三区| 婷婷亚洲欧美| 日韩成人伦理影院| 精品久久久久久久久久久久久| 中文亚洲av片在线观看爽| 欧美一级a爱片免费观看看| aaaaa片日本免费| 国产淫片久久久久久久久| 日本爱情动作片www.在线观看 | 亚洲无线在线观看| 在线a可以看的网站| 久久草成人影院| 日韩三级伦理在线观看| 欧美高清成人免费视频www| 69av精品久久久久久| 中文字幕av成人在线电影| 又爽又黄a免费视频| 国产精品女同一区二区软件| 色哟哟哟哟哟哟| 国产综合懂色| 狠狠狠狠99中文字幕| 国产一级毛片七仙女欲春2| 夜夜夜夜夜久久久久| 黑人高潮一二区| 少妇熟女aⅴ在线视频| 人妻久久中文字幕网| 国产成人a∨麻豆精品| 国国产精品蜜臀av免费| 国产精品,欧美在线| 国产精品av视频在线免费观看| 日韩一本色道免费dvd| 精品熟女少妇av免费看| 少妇高潮的动态图| 小蜜桃在线观看免费完整版高清| 又黄又爽又刺激的免费视频.| 亚洲人成网站在线播| 成人二区视频| 色av中文字幕| 成人鲁丝片一二三区免费| 久99久视频精品免费| 久久九九热精品免费| 亚洲五月天丁香| 偷拍熟女少妇极品色| 免费观看人在逋| 久久久精品94久久精品| 亚洲熟妇中文字幕五十中出| 精品国产三级普通话版| 五月伊人婷婷丁香| 国产精品亚洲一级av第二区| 天堂av国产一区二区熟女人妻| 乱码一卡2卡4卡精品| 国产在视频线在精品| 99热这里只有精品一区| 性欧美人与动物交配| 插阴视频在线观看视频| 中文字幕熟女人妻在线| 国产一区二区激情短视频| 久久人人爽人人爽人人片va| 看非洲黑人一级黄片| 日本与韩国留学比较| 国产男靠女视频免费网站| 国模一区二区三区四区视频| 91狼人影院| 99在线人妻在线中文字幕| 国产中年淑女户外野战色| 久久久久久久久大av| 国产午夜精品久久久久久一区二区三区 | 免费看美女性在线毛片视频| 在线观看美女被高潮喷水网站| 中文字幕熟女人妻在线| 少妇猛男粗大的猛烈进出视频 | 亚洲av中文字字幕乱码综合| 国产精品人妻久久久久久| 婷婷亚洲欧美| av在线蜜桃| 搡老岳熟女国产| 欧美精品国产亚洲| 亚洲七黄色美女视频| 国产精品嫩草影院av在线观看| 极品教师在线视频| 午夜精品一区二区三区免费看| 国产av麻豆久久久久久久| 欧美激情国产日韩精品一区| 国产精品亚洲一级av第二区| 国产av在哪里看| 亚洲国产精品sss在线观看| 免费电影在线观看免费观看| av在线天堂中文字幕| 欧美成人免费av一区二区三区| 高清毛片免费看| 日日摸夜夜添夜夜添av毛片| 国内精品美女久久久久久| 99热网站在线观看| 校园人妻丝袜中文字幕| 亚洲欧美日韩高清专用| 亚洲美女搞黄在线观看 | 久久久久久伊人网av| 日日摸夜夜添夜夜爱| 午夜精品在线福利| 午夜精品国产一区二区电影 | 国产伦精品一区二区三区四那| 三级经典国产精品| 日本色播在线视频| 97在线视频观看| 美女cb高潮喷水在线观看| 精品午夜福利视频在线观看一区| 九九爱精品视频在线观看| 日韩精品中文字幕看吧| 我的老师免费观看完整版| 五月玫瑰六月丁香| 卡戴珊不雅视频在线播放| 综合色av麻豆| 国内精品宾馆在线| 免费看av在线观看网站| 久久人人爽人人片av| 精品日产1卡2卡| 人妻久久中文字幕网| 国产精品电影一区二区三区| 天天躁夜夜躁狠狠久久av| 国产蜜桃级精品一区二区三区| 女生性感内裤真人,穿戴方法视频| 听说在线观看完整版免费高清| 非洲黑人性xxxx精品又粗又长| 国产综合懂色| 欧美+日韩+精品| 亚洲成人久久性| 色综合亚洲欧美另类图片| 噜噜噜噜噜久久久久久91| 18禁裸乳无遮挡免费网站照片| 成熟少妇高潮喷水视频| av在线亚洲专区| 日本黄色视频三级网站网址| 日本与韩国留学比较| 国产精品爽爽va在线观看网站| 精品久久久久久久末码| 欧美丝袜亚洲另类| 国产在视频线在精品| 日韩,欧美,国产一区二区三区 | 九色成人免费人妻av| 91久久精品国产一区二区成人| 在线免费十八禁| 午夜精品国产一区二区电影 | 午夜免费男女啪啪视频观看 | 色播亚洲综合网| 简卡轻食公司| 久久久久精品国产欧美久久久| 久久久欧美国产精品| 国产亚洲精品久久久久久毛片| 亚洲欧美日韩无卡精品| 亚洲无线在线观看| 久久久久性生活片| 国产成人a区在线观看| 国产精品99久久久久久久久| 亚洲激情五月婷婷啪啪| 国产亚洲91精品色在线| 性插视频无遮挡在线免费观看| 在线看三级毛片| 精品一区二区三区av网在线观看| 日韩精品中文字幕看吧| 国产精品久久久久久亚洲av鲁大| av免费在线看不卡| 狠狠狠狠99中文字幕| 成人性生交大片免费视频hd| 精品久久久久久久久久久久久| 看免费成人av毛片| 国内少妇人妻偷人精品xxx网站| 国产高清视频在线观看网站| 69av精品久久久久久| 国产亚洲欧美98| 国产精品电影一区二区三区| 老熟妇乱子伦视频在线观看| 成人三级黄色视频| 亚洲国产欧洲综合997久久,| 午夜福利视频1000在线观看| 国产又黄又爽又无遮挡在线| 国产激情偷乱视频一区二区| 国产极品精品免费视频能看的| 亚洲欧美成人综合另类久久久 | 韩国av在线不卡| 亚洲av成人精品一区久久| 国产午夜福利久久久久久| 99久久中文字幕三级久久日本| 免费观看的影片在线观看| 成人亚洲精品av一区二区| 成人美女网站在线观看视频| av天堂中文字幕网| 日韩三级伦理在线观看| 最近手机中文字幕大全| 男女那种视频在线观看| 国产精品1区2区在线观看.| 91av网一区二区| 免费搜索国产男女视频| 免费看日本二区| 天天一区二区日本电影三级| 蜜臀久久99精品久久宅男| 一级毛片电影观看 | 大型黄色视频在线免费观看| 久久精品久久久久久噜噜老黄 | 精品人妻一区二区三区麻豆 | 久久欧美精品欧美久久欧美| 十八禁网站免费在线| 亚洲成人久久爱视频| 亚洲天堂国产精品一区在线| 精品不卡国产一区二区三区| 久久这里只有精品中国| 中国美女看黄片| 亚洲av不卡在线观看| 色5月婷婷丁香| 午夜影院日韩av| 色综合站精品国产| 国产麻豆成人av免费视频| 日本一本二区三区精品| .国产精品久久| 亚洲国产欧洲综合997久久,| 精品欧美国产一区二区三| 麻豆一二三区av精品| 男人舔奶头视频| 18禁在线播放成人免费| 亚洲精品乱码久久久v下载方式| 女同久久另类99精品国产91| 国产精品三级大全| 中文字幕精品亚洲无线码一区| 黄色欧美视频在线观看| av在线蜜桃| 成人特级av手机在线观看| 日本黄色片子视频| 婷婷亚洲欧美| 97超视频在线观看视频| 亚洲av电影不卡..在线观看| av在线蜜桃| 免费人成视频x8x8入口观看| 精品久久久久久久久久免费视频| 中国国产av一级| 性色avwww在线观看| 成人综合一区亚洲| 日韩欧美精品v在线| 亚洲va在线va天堂va国产| 成人一区二区视频在线观看| 欧美日韩精品成人综合77777| 禁无遮挡网站| 亚洲性久久影院| 色吧在线观看| 97超级碰碰碰精品色视频在线观看| 网址你懂的国产日韩在线| 床上黄色一级片| 亚洲国产精品久久男人天堂| 久久热精品热| 国产精品久久久久久av不卡| 国产在线男女| 国产蜜桃级精品一区二区三区| 在线免费观看不下载黄p国产| 白带黄色成豆腐渣| 亚洲精品粉嫩美女一区| 亚洲欧美精品综合久久99| 12—13女人毛片做爰片一| 女人十人毛片免费观看3o分钟| 美女被艹到高潮喷水动态| 亚洲美女视频黄频| 99热全是精品| 国产成人精品久久久久久| 69人妻影院| 日韩一区二区视频免费看| 国产伦精品一区二区三区视频9| 久久人人爽人人片av| 亚洲国产高清在线一区二区三| 国产精品三级大全| 日本撒尿小便嘘嘘汇集6| 亚洲人成网站在线播放欧美日韩| 日韩强制内射视频| 亚洲天堂国产精品一区在线| 国产精品福利在线免费观看| 插阴视频在线观看视频| 国产成人一区二区在线| 亚洲自偷自拍三级| 成熟少妇高潮喷水视频| 亚洲av中文av极速乱| 亚洲精品国产av成人精品 | 久久久成人免费电影| 色视频www国产| av在线蜜桃| 亚洲精品国产av成人精品 | 变态另类丝袜制服| 成熟少妇高潮喷水视频| 一进一出抽搐gif免费好疼| 国产精品久久久久久av不卡| 欧美高清成人免费视频www| 日韩欧美三级三区| 成人二区视频| 久久久久久九九精品二区国产| 欧美zozozo另类| 国产精品电影一区二区三区| 免费av毛片视频| 国产单亲对白刺激| 国产精品综合久久久久久久免费| 日韩av在线大香蕉| 日韩在线高清观看一区二区三区| 亚洲综合色惰| 如何舔出高潮| 大型黄色视频在线免费观看| 女人十人毛片免费观看3o分钟| 在线国产一区二区在线| 精品乱码久久久久久99久播| 最新在线观看一区二区三区| 色哟哟哟哟哟哟| 成人无遮挡网站| 久久鲁丝午夜福利片| 久久精品综合一区二区三区| 男女那种视频在线观看| 18禁在线播放成人免费| 日韩一区二区视频免费看| 丰满人妻一区二区三区视频av| 精品日产1卡2卡| 国产男靠女视频免费网站| 亚洲国产精品成人久久小说 | h日本视频在线播放| 丰满乱子伦码专区| 国内精品久久久久精免费| 久久久成人免费电影| 晚上一个人看的免费电影| 久久欧美精品欧美久久欧美| 一本一本综合久久| 在线观看66精品国产| 两个人视频免费观看高清| 搞女人的毛片| 综合色av麻豆| 观看免费一级毛片| 国产三级在线视频| av在线播放精品| 欧美中文日本在线观看视频| 色尼玛亚洲综合影院| 三级男女做爰猛烈吃奶摸视频| 日本 av在线| 中国国产av一级| 亚洲欧美中文字幕日韩二区| 国产毛片a区久久久久| 国产成人福利小说| 久久精品国产亚洲av天美| 国产一区二区三区av在线 | 99riav亚洲国产免费| 国语自产精品视频在线第100页| 国产精品爽爽va在线观看网站| 精品国内亚洲2022精品成人| 美女黄网站色视频| 直男gayav资源| 日本 av在线| 中文字幕免费在线视频6| 国产一区二区三区av在线 | 亚洲va在线va天堂va国产| 欧美不卡视频在线免费观看| 亚洲最大成人av| a级一级毛片免费在线观看| 亚洲精品一区av在线观看| 校园人妻丝袜中文字幕| 中文字幕人妻熟人妻熟丝袜美| 国产精品综合久久久久久久免费| 嫩草影院入口| 18+在线观看网站| 久久久久久久久大av| 成人av在线播放网站| 欧美激情国产日韩精品一区| 日日撸夜夜添| 韩国av在线不卡| 国产白丝娇喘喷水9色精品| 波多野结衣高清作品| 麻豆av噜噜一区二区三区| 白带黄色成豆腐渣| 国产欧美日韩精品亚洲av| 性欧美人与动物交配| 欧美日韩一区二区视频在线观看视频在线 | 丰满人妻一区二区三区视频av| 亚洲人与动物交配视频| 白带黄色成豆腐渣| 色哟哟·www| 亚洲国产精品久久男人天堂| 国产精品综合久久久久久久免费| 欧美日韩乱码在线| av在线老鸭窝| 亚洲成人精品中文字幕电影| 日韩高清综合在线| 赤兔流量卡办理| 国产蜜桃级精品一区二区三区| 少妇人妻精品综合一区二区 | 男人的好看免费观看在线视频| 国产熟女欧美一区二区| 18禁裸乳无遮挡免费网站照片| 欧美日韩综合久久久久久| 男女视频在线观看网站免费| 亚洲精品一卡2卡三卡4卡5卡| 国产精品久久视频播放| 小说图片视频综合网站| 男人的好看免费观看在线视频| 成年女人看的毛片在线观看| 男插女下体视频免费在线播放| 三级男女做爰猛烈吃奶摸视频| 欧美成人一区二区免费高清观看| 亚洲av二区三区四区| 18禁裸乳无遮挡免费网站照片|