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

    不同環(huán)境下噴霧燃燒碳煙演化歷程的數(shù)值模擬

    2022-08-25 02:03:12王晨辰鄭尊清堯命發(fā)
    燃燒科學(xué)與技術(shù) 2022年4期
    關(guān)鍵詞:當(dāng)量環(huán)境溫度體積

    王晨辰,鄭尊清,王?滸,堯命發(fā),章?嚴(yán)

    不同環(huán)境下噴霧燃燒碳煙演化歷程的數(shù)值模擬

    王晨辰,鄭尊清,王?滸,堯命發(fā),章?嚴(yán)

    (天津大學(xué)內(nèi)燃機(jī)燃燒學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室)

    針對(duì)柴油發(fā)動(dòng)機(jī)在不同負(fù)荷下碳煙排放特征的演變和發(fā)展,依據(jù)ECN噴霧定容燃燒實(shí)驗(yàn)中碳煙的測(cè)量數(shù)據(jù),基于OpenFOAM開(kāi)源軟件噴霧燃燒求解器,耦合半經(jīng)驗(yàn)現(xiàn)象學(xué)碳煙預(yù)測(cè)模型.通過(guò)-關(guān)系圖呈現(xiàn)環(huán)境因素改變對(duì)碳煙形成及氧化各階段的影響,分析各子過(guò)程的分布區(qū)域和演化規(guī)律,進(jìn)而提出環(huán)境條件對(duì)碳煙生成及演化歷程的影響機(jī)制.研究結(jié)果表明:現(xiàn)象學(xué)模型本身能夠較好地反映碳煙的分布情況.碳煙生長(zhǎng)主要發(fā)生在當(dāng)量比>3區(qū)域;而碳煙發(fā)生氧化的區(qū)域,也是OH主要分布的區(qū)間,大多位于當(dāng)量比<2的位置.環(huán)境溫度和密度的增加初期促進(jìn)了碳煙的成核,加劇了顆粒的聚并速率,使得顆粒數(shù)顯著下降;而環(huán)境氧體積分?jǐn)?shù)的增加抑制了碳煙顆粒數(shù)的增長(zhǎng)和聚并速率,使碳煙整體尺寸趨于小尺寸、高密度的分布情況.

    碳煙;定容燃燒;OpenFOAM;數(shù)值模擬

    隨著環(huán)境問(wèn)題日益受到關(guān)注,有害物的排放已然成為制約內(nèi)燃機(jī)發(fā)展的重要因素之一.碳煙顆粒(PM)是內(nèi)燃機(jī)有害排放的重要組成部分,其形成原因主要是由于燃料的不完全燃燒.內(nèi)燃機(jī)排放的碳煙顆粒尺寸通常為納米級(jí),雖然在總排放物中質(zhì)量占比不高,但數(shù)量很大,呈現(xiàn)典型的彌散性,容易長(zhǎng)期懸浮于空氣中,導(dǎo)致霧霾的產(chǎn)生.過(guò)去的排放法規(guī)主要對(duì)碳煙排放的質(zhì)量進(jìn)行限制,在歐Ⅴ及之后排放法規(guī)中,除對(duì)內(nèi)燃機(jī)顆粒排放的質(zhì)量上有了嚴(yán)格的限制外,更是要求對(duì)碳煙排放的數(shù)密度進(jìn)行控制.近期有關(guān)內(nèi)燃機(jī)排放碳煙顆粒物的研究中[1],顆粒物粒徑和數(shù)密度逐漸成為關(guān)注的重點(diǎn),也必將是今后內(nèi)燃機(jī)排放的主要研究方向之一.

    碳煙的形成包含了復(fù)雜的物理和化學(xué)過(guò)程,為深入分析碳煙生成機(jī)理進(jìn)而為碳煙排放控制提供理論指導(dǎo),需要構(gòu)建碳煙生成過(guò)程的理論模型.在碳煙排放預(yù)測(cè)模型方面,國(guó)內(nèi)外學(xué)者開(kāi)展了大量研究工作,發(fā)展了經(jīng)驗(yàn)?zāi)P?、半?jīng)驗(yàn)?zāi)P鸵约霸敿?xì)模型等不同適用范圍的碳煙預(yù)測(cè)模型[2].早期最為成熟并被廣泛應(yīng)用的碳煙模型為Hiroyasu等[3]提出的兩步法模型,通過(guò)Arrhenius公式,依據(jù)燃料的質(zhì)量分布得到碳煙的生成速率,依據(jù)碳煙的質(zhì)量分布計(jì)算碳煙的氧化速率.Wang等[4-7]在研究PAHs的生成過(guò)程時(shí),提出了脫氫加乙炔(HACA)的反應(yīng)機(jī)理,并在后續(xù)碳煙模型的表面生長(zhǎng)過(guò)程得到了應(yīng)用和發(fā)展.Tao等[8]首次提出以乙炔為前驅(qū)物的多步現(xiàn)象學(xué)模型,并在此基礎(chǔ)上發(fā)展了包含燃料熱分解、前驅(qū)物形成、碳煙初始核形成等在內(nèi)的9步法碳煙模型[9].Jia等[10]提出了一個(gè)改進(jìn)的以乙炔為前驅(qū)物的多步法現(xiàn)象學(xué)碳煙模型,在乙烯對(duì)沖火焰中進(jìn)行了實(shí)驗(yàn)驗(yàn)證,并在柴油均質(zhì)壓燃(HCCI)模式下的碳煙排放研究中得到應(yīng)用.Leung等[11]為模擬分析柴油HCCI模式碳煙的形成和氧化過(guò)程,將簡(jiǎn)化的正庚烷機(jī)理與碳煙現(xiàn)象學(xué)模型進(jìn)行耦合,并用新的模型與乙烯火焰對(duì)沖實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比,很好地預(yù)測(cè)了乙烯火焰中碳煙的產(chǎn)量和分布情況. Vishwanathan等[12]基于Leung發(fā)展的半經(jīng)驗(yàn)碳煙模型,改用A4(pyrene)作為碳煙成核的前驅(qū)物,準(zhǔn)確預(yù)測(cè)了碳煙隨環(huán)境溫度和壓力變化的發(fā)展趨勢(shì).

    目前,基于半經(jīng)驗(yàn)現(xiàn)象學(xué)模型的研究工作取得很大進(jìn)展,可以在一定程度上分析碳煙生成氧化過(guò)程的影響機(jī)制和變化規(guī)律.本文基于Vishwanathan等改進(jìn)的碳煙現(xiàn)象學(xué)模型,應(yīng)用本課題組發(fā)展的TRF/PAH簡(jiǎn)化化學(xué)反應(yīng)動(dòng)力學(xué)機(jī)理[13],補(bǔ)充了從苯到更大PAH芘(A4)的簡(jiǎn)化機(jī)理描述.基于以上兩種模型開(kāi)展正庚烷燃料定容燃燒噴霧實(shí)驗(yàn)?zāi)M研究,分析初始環(huán)境條件變化引起的碳煙生成氧化等過(guò)程的發(fā)展趨勢(shì),為后續(xù)開(kāi)展發(fā)動(dòng)機(jī)排放的數(shù)值模擬工作做準(zhǔn)備.

    1?碳煙模型理論

    以Tao等[8]為代表發(fā)展的半經(jīng)驗(yàn)?zāi)P屯ㄟ^(guò)求解質(zhì)量分?jǐn)?shù)和顆粒數(shù)密度兩個(gè)輸運(yùn)方程,描述包括氣相PAH成核、表面生長(zhǎng)、顆粒間聚并以及氧化在內(nèi)的碳煙發(fā)展過(guò)程.在Vishwanathan模型中,碳煙的形成過(guò)程被描述為有關(guān)碳煙質(zhì)量分?jǐn)?shù)soot和顆粒數(shù)密度soot的兩個(gè)輸運(yùn)方程,這里將soot輸運(yùn)方程兩邊所有項(xiàng)同除以一個(gè)歸一化因子norm,得到的兩個(gè)輸運(yùn)方程如式(1)、(2)所示:

    式中:t為湍流施密特?cái)?shù),本文設(shè)為0.7;代表單位體積混合物內(nèi)碳煙的質(zhì)量;代表載流體的動(dòng)力學(xué)黏度;nuc為單位質(zhì)量混合物中碳煙的數(shù)密度;norm為歸一化因子.為簡(jiǎn)化計(jì)算使用,本模型設(shè)為1015;而表示碳煙的熱泳系數(shù),采用文獻(xiàn)[11]提供的公式計(jì)算得到:

    其中參數(shù)為0.9.式(1)、(2)的最后一項(xiàng)為源項(xiàng),用來(lái)描述碳煙的顆粒動(dòng)力學(xué)過(guò)程:

    顆粒的動(dòng)力學(xué)描述包含成核、表面生長(zhǎng)、聚并以及OH和O2的氧化,其速率分別由式(5)中的1、2、3、4和5表示,c(s)為C元素的摩爾質(zhì)量,nuci為成核前驅(qū)物的摩爾質(zhì)量.大量研究表明[8-12],采用A4作為碳煙成核的前驅(qū)物,用芘的二聚化描述碳煙的成核過(guò)程,可以實(shí)現(xiàn)對(duì)碳煙發(fā)展趨勢(shì)的準(zhǔn)確捕捉.因此這里nuci采用A4的摩爾質(zhì)量.表面生長(zhǎng)遵守Wang等[7]的HACA規(guī)則.顆粒的凝并過(guò)程采用Smoluchowski方程描述,假設(shè)顆粒尺寸是均勻分布的,且凝并過(guò)程與顆粒的體積分?jǐn)?shù)v和顆粒的數(shù)量有關(guān),因此顆粒數(shù)由于凝并減少的速率如下:

    其中co表示碰撞速率:

    假設(shè)碳煙顆粒是球形,則為2,B為玻爾茲曼常數(shù):

    碳煙受O2的氧化作用利用Nagle和Strickland-Constable(NSC)[14]模型描述,模型將碳煙顆粒表面的活性點(diǎn)位分為兩類(lèi),一類(lèi)是活性較強(qiáng)的a點(diǎn),占百分比為;一類(lèi)是活性較弱的b點(diǎn),占比為1-.表面氧化速率為:

    其中、a、b和z等變量分別通過(guò)以下方程求得:

    其中O2為氧氣的分壓.

    碳煙受OH的氧化通過(guò)Neoh等[15]的經(jīng)驗(yàn)公式描述:

    2?模型驗(yàn)證

    2.1?噴霧模型的驗(yàn)證

    數(shù)值模擬研究基于OpenFOAM開(kāi)源軟件平臺(tái),OpenFOAM是一款基于C++語(yǔ)言開(kāi)發(fā)和編譯的通用CFD平臺(tái),具有開(kāi)源、擴(kuò)展性強(qiáng)和接口豐富等優(yōu)點(diǎn),可以修改和制定所需的模型和求解器,適應(yīng)不同的需求.

    模擬參考的實(shí)驗(yàn)設(shè)置和數(shù)據(jù)來(lái)源于美國(guó)圣地亞國(guó)家實(shí)驗(yàn)室提供的正庚烷單孔定容噴霧實(shí)驗(yàn)Spray-H[17].本文中,容彈被簡(jiǎn)化為高 108mm,直徑80mm的圓柱體,噴油器設(shè)置在容彈頂部中心處,噴孔直徑0.1mm,噴油錐角為10°.網(wǎng)格數(shù)約28萬(wàn),網(wǎng)格在靠近軸心處進(jìn)行局部加密,最小的網(wǎng)格尺寸為0.25mm,經(jīng)網(wǎng)格無(wú)關(guān)性驗(yàn)證保證精度,定容燃燒彈模型如圖1所示.

    圖1?定容燃燒彈幾何模型

    本文使用OpenFOAM平臺(tái)提供的噴霧燃燒求解器sprayFoam,耦合經(jīng)驗(yàn)、半經(jīng)驗(yàn)碳煙模型,選正庚烷為燃料,其余噴霧、燃燒等子模型如表1所列.

    表1?主要的CFD模型

    Tab.1?Main sub-models in CFD simulation

    模擬燃燒前需要對(duì)噴霧基準(zhǔn)算例進(jìn)行標(biāo)定,標(biāo)定工況為初始?jí)毫?.33MPa、環(huán)境溫度1000K、環(huán)境氧體積分?jǐn)?shù)為0的無(wú)反應(yīng)工況,其余條件見(jiàn)表2.

    表2?定容燃燒彈環(huán)境參數(shù)

    Tab.2 Environmental parameters of constant volume combustion vessel

    根據(jù)選定的模型和邊界條件,通過(guò)數(shù)值模擬得到Spray-H基準(zhǔn)工況下噴霧的液相貫穿距和氣相貫穿距,并與實(shí)驗(yàn)值進(jìn)行對(duì)比.如圖2所示.模擬的結(jié)果與實(shí)驗(yàn)測(cè)量值基本吻合.

    從模擬得到的數(shù)據(jù)中提取沿軸線的混合分?jǐn)?shù),以及位于軸線坐標(biāo)20mm處、沿徑向混合分?jǐn)?shù)曲線,將兩者與實(shí)驗(yàn)值進(jìn)行對(duì)比,如圖3(a)、(b)所示,圖中虛線表示允許的誤差范圍,結(jié)果表明誤差基本在5%范圍內(nèi).噴霧過(guò)程的模擬結(jié)果與實(shí)驗(yàn)值基本相符,后續(xù)模擬研究均基于此噴霧設(shè)定開(kāi)展.

    圖3 無(wú)反應(yīng)條件下實(shí)驗(yàn)與模擬得到的穩(wěn)定混合后燃料摩爾分?jǐn)?shù)分布

    2.2?碳煙模型對(duì)比和評(píng)估

    模型的選取對(duì)正確預(yù)測(cè)碳煙的發(fā)展過(guò)程至關(guān)重要.這里對(duì)幾種主流碳煙模型的模擬結(jié)果進(jìn)行評(píng)估,如圖4所示.根據(jù)實(shí)驗(yàn)數(shù)據(jù),選取環(huán)境氧體積分?jǐn)?shù)為15%、初始溫度為1000K、壓力為4.27MPa的條件為基準(zhǔn),取軸線上的碳煙體積分?jǐn)?shù)分布曲線,即圖4中黑色曲線,與相同邊界條件下不同模型模擬得到的碳煙體積分?jǐn)?shù)分布曲線進(jìn)行對(duì)比.

    圖4?幾種碳煙模型預(yù)測(cè)結(jié)果對(duì)比

    兩步法模型由于依據(jù)燃料的質(zhì)量計(jì)算碳煙的生成速率,忽略了成核等碳煙形成的過(guò)程,因此預(yù)測(cè)得到的碳煙分布和數(shù)值和實(shí)驗(yàn)相比,有較大差異.碳煙分布位置偏上游,分布形狀并不相似;半經(jīng)驗(yàn)多步現(xiàn)象學(xué)模型預(yù)測(cè)的碳煙分布與實(shí)驗(yàn)值更接近.Tao的九步法模型由于假設(shè)乙炔為碳煙成核前驅(qū)物的來(lái)源,預(yù)測(cè)得到的碳煙峰值偏上游;Vishwanathan等[12]采用A4作為碳煙的成核前驅(qū)物,預(yù)測(cè)的碳煙分布形狀和峰值位置與實(shí)驗(yàn)結(jié)果最接近,但整體輪廓顯得更寬.

    結(jié)合內(nèi)燃機(jī)缸內(nèi)熱力學(xué)情況,考慮以不同初始環(huán)境發(fā)動(dòng)機(jī)缸內(nèi)壓力、氧體積分?jǐn)?shù)以及環(huán)境溫度條件為研究背景,以正庚烷為柴油替代燃料,研究碳煙生成氧化過(guò)程演化趨勢(shì).?dāng)M選取表3所列發(fā)動(dòng)機(jī)環(huán)境進(jìn)行模擬研究,將不同初始環(huán)境氧體積分?jǐn)?shù)及環(huán)境壓力條件下預(yù)測(cè)得到的碳煙空間分布情況與實(shí)驗(yàn)值對(duì)比,結(jié)果如圖5所示.

    表3?初始環(huán)境參數(shù)

    Tab.3?Initial environmental parameters

    圖5中(a)、(b)和(c)分別對(duì)應(yīng)環(huán)境溫度1000K、環(huán)境密度為14.8kg/m3條件下,初始氧體積分?jǐn)?shù)為10%、15%和21%三種工況點(diǎn)中模擬和測(cè)量得到的碳煙空間分布;圖5(d)為相同環(huán)境溫度下,初始氧體積分?jǐn)?shù)為15%、密度30kg/m3條件下模擬和測(cè)量得到的碳煙空間分布.可以看出,Vishwanathan模型能夠準(zhǔn)確描述Spray-H碳煙隨環(huán)境氧體積分?jǐn)?shù)、密度改變情況下,碳煙的體積分?jǐn)?shù)、峰值位置和空間分布的演化趨勢(shì).

    相比實(shí)驗(yàn)值,模擬得到的分布云圖略顯狹長(zhǎng),導(dǎo)致這一現(xiàn)象的原因,可能是由于氧化模型預(yù)測(cè)的O2和OH對(duì)碳煙的氧化速率較實(shí)驗(yàn)值要低.低氧體積分?jǐn)?shù)下預(yù)測(cè)得到的碳煙輪廓較實(shí)驗(yàn)值與高氧體積分?jǐn)?shù)情況誤差較大,可能原因是模型預(yù)測(cè)得到的氧化速率相比實(shí)驗(yàn)情況,對(duì)氧體積分?jǐn)?shù)的敏感性更高,低氧體積分?jǐn)?shù)下氧化速率存在一定的誤差.相比于環(huán)境密度為14.8kg/m3的工況,環(huán)境密度為30kg/m3時(shí)預(yù)測(cè)得到的碳煙體積分?jǐn)?shù)峰值與實(shí)驗(yàn)值相比偏低,且峰值位置更偏下游,這與Vishwanathan等的研究結(jié)論一致.

    3?初始環(huán)境變化對(duì)碳煙分布的影響規(guī)律

    3.1?碳煙生成氧化各過(guò)程的分布特性

    碳煙的生成與當(dāng)?shù)貕毫?、溫度和混合氣成分密切相關(guān).在Vishwanathan模型中,由于碳煙的成核速率依賴(lài)于A4的體積分?jǐn)?shù),表面生長(zhǎng)速率依賴(lài)于C2H2的體積分?jǐn)?shù),而碳煙的氧化與當(dāng)?shù)豋2和OH的體積分?jǐn)?shù)密切相關(guān).由于這些關(guān)鍵氣相成分各自分布在不同的當(dāng)量比()和溫度()范圍,因此碳煙生成各過(guò)程發(fā)生的主要分布區(qū)域,其熱力學(xué)狀態(tài)與混合氣體積分?jǐn)?shù)并不一致.

    在環(huán)境密度為14.8kg/m3、溫度1000K、氧體積分?jǐn)?shù)15%的初始條件下,將各網(wǎng)格中A4、C2H2和OH的質(zhì)量分?jǐn)?shù),以其所在網(wǎng)格對(duì)應(yīng)的和值為坐標(biāo),繪制成如圖6(a)、(b)和(c)所示的氣泡分布圖,其中氣泡的尺寸與相應(yīng)組分的質(zhì)量分?jǐn)?shù)大小成正比.采用同樣方法可得到碳煙顆粒數(shù)密度,以及成核、表面生長(zhǎng)和氧化等過(guò)程的質(zhì)量變化速率在-圖中的分布.為便于分析,僅在圖7中顯示顆粒數(shù)密度在-圖中的尺寸分布,而成核、表面生長(zhǎng)等過(guò)程質(zhì)量變化量的主要分布區(qū)域(依據(jù)相應(yīng)氣泡尺寸)分別用綠色、藍(lán)色和紅色的輪廓線圈出.

    對(duì)圖6中的-分布情況進(jìn)行對(duì)比,由圖6(a)、(b)和(c)可知,模型中碳煙成核的前驅(qū)物A4和表面生長(zhǎng)的主要成分C2H2均分布于>2的區(qū)域,與此相反,OH主要分布于<3的高溫區(qū).C2H2在1400K<<2000K區(qū)間的分布相對(duì)均勻,而A4的分布更集中些.這些分布特征與圖7由上述輪廓線圈出的成核、表面生長(zhǎng)以及氧化等質(zhì)量變化速率的主要分布區(qū)域相一致.

    從圖7中還可以發(fā)現(xiàn),O2對(duì)碳煙的氧化不僅存在于低當(dāng)量比、高溫區(qū)域,同時(shí)也有少量存在于高當(dāng)量比區(qū)域,這與Dec[18]提出的概念模型稍有出入.實(shí)際上,Chishty等[19]的模擬結(jié)果均存在這部分由O2氧化的區(qū)域,其位于碳煙分布輪廓的上游處.

    圖7?碳煙生成氧化各個(gè)階段在Φ-T圖中的主要分布

    3.2?環(huán)境溫度對(duì)碳煙生成歷程的影響

    環(huán)境溫度對(duì)碳煙的生成和氧化歷程均有影響.將初始環(huán)境密度為14.8kg/m3,氧體積分?jǐn)?shù)21%條件下,環(huán)境溫度分別為800K、900K、1000K以及1200K時(shí),Spray-H的碳煙分布繪制如圖8所示.

    圖8?不同環(huán)境溫度下碳煙分布區(qū)域的演化

    分析圖8可以發(fā)現(xiàn),隨著環(huán)境溫度從800K增加到1200K,當(dāng)量比的峰值從3的位置依次增加到10以上.而在當(dāng)量比<2的區(qū)域,即對(duì)應(yīng)圖7中OH氧化速率主要分布的區(qū)域,在初始環(huán)境溫度增加之后,整體分布區(qū)域當(dāng)量比并沒(méi)有發(fā)生明顯變化,但溫度發(fā)生了變化:圖8中0=800K對(duì)應(yīng)的散點(diǎn),該區(qū)域(<2的區(qū)域)整體溫度約2500K,到0=1200K對(duì)應(yīng)的散點(diǎn),該區(qū)域整體溫度峰值約2700K.

    將環(huán)境氧體積分?jǐn)?shù)為21%、環(huán)境密度為14.8kg/m3時(shí),初始溫度分別為1000K和1200K條件下,碳煙顆粒的增長(zhǎng)速率繪制成圖9(a)、(b)中的藍(lán)色曲線,而將碳煙顆粒由粒子間碰撞聚并導(dǎo)致數(shù)量減少的速率繪制成圖9中的紅色曲線.

    圖9 環(huán)境溫度不同時(shí)顆??倲?shù)變化速率與顆粒間聚并速率

    由圖9(a)可知,環(huán)境溫度為1000K時(shí),碳煙在噴射始點(diǎn)后約0.7ms產(chǎn)生,稍有延遲后聚并發(fā)生.起初碳煙增長(zhǎng)速率很快,達(dá)到每秒約7×1015后,隨著聚并效應(yīng)的增強(qiáng),顆粒凈增加速率逐步衰減.如圖9(b)所示,當(dāng)環(huán)境溫度=1200K時(shí),碳煙成核發(fā)生在噴射始點(diǎn)后約0.5ms之前,顆粒聚并速率急劇增加到最大值,顆粒數(shù)增長(zhǎng)速率峰值增加且提前到達(dá),隨后急劇衰減直到趨于穩(wěn)定.初始溫度的提高對(duì)碳煙的成核速率和聚并速率都起促進(jìn)作用,且對(duì)聚并的促進(jìn)作用明顯高于對(duì)成核的促進(jìn)作用.

    3.3?初始氧體積分?jǐn)?shù)對(duì)碳煙生成歷程的影響

    由于內(nèi)燃機(jī)通常需要改變EGR率實(shí)現(xiàn)對(duì)燃燒和排放的控制,其中一個(gè)重要作用來(lái)自于EGR對(duì)氧體積分?jǐn)?shù)的改變.鑒于此,分析了環(huán)境氧體積分?jǐn)?shù)為12%、15%和21%條件下,碳煙顆粒隨氧體積分?jǐn)?shù)增加的整體發(fā)展趨勢(shì).

    在-圖中繪制了環(huán)境溫度為1000K、環(huán)境密度為14.8kg/m3條件下,初始氧體積分?jǐn)?shù)為15%時(shí)碳煙顆粒數(shù)密度的氣泡分布圖.并用前述的輪廓線圈出該條件下各碳煙生長(zhǎng)氧化過(guò)程對(duì)應(yīng)的質(zhì)量變化速率的主要分布區(qū)域,如圖10(a)所示.為了對(duì)比,選取同樣環(huán)境溫度、密度下,初始氧體積分?jǐn)?shù)為12%、15%和21%時(shí)碳煙分布情況,分別繪制成圖10(b)中的紅色、黑色以及藍(lán)色散點(diǎn)圖.

    圖10 碳煙生成氧化各個(gè)階段分布及隨氧體積分?jǐn)?shù)改變的演化趨勢(shì)

    由圖10(b)可知,當(dāng)環(huán)境氧體積分?jǐn)?shù)從12%增加到21%時(shí),整個(gè)碳煙的分布曲線呈現(xiàn)出向高溫區(qū)和高當(dāng)量比區(qū)域平移的趨勢(shì).對(duì)應(yīng)于圖10(a)中的碳煙成核、表面生長(zhǎng)過(guò)程的質(zhì)量變化速率主要分布區(qū)域(即>3的部分),整體溫度升高,當(dāng)量比增大,峰值從約=6增加到=7附近;而在碳煙發(fā)生OH氧化的主要區(qū)域(即<2的部分),當(dāng)量比的增加并不明顯,相比之下溫度升高顯著:隨著氧體積分?jǐn)?shù)從12%增大到21%,該區(qū)域內(nèi)溫度整體上升了約500K.

    分析環(huán)境溫度為1000K、環(huán)境密度為14.8kg/m3條件下,初始氧體積分?jǐn)?shù)分別為15%和21%時(shí)碳煙顆粒的增長(zhǎng)速率曲線(圖11(a)、(b)中的藍(lán)色曲線),以及由于粒子間碰撞聚并過(guò)程,導(dǎo)致碳煙顆粒數(shù)量減少的速率曲線(圖11(a)、(b)中的紅色曲線).

    圖11(a)中,即初始氧體積分?jǐn)?shù)為15%時(shí),碳煙顆粒數(shù)增長(zhǎng)的速率峰值約在3ms處,與圖11(b)中初始氧體積分?jǐn)?shù)為21%條件下的速率峰值時(shí)刻1.6ms相比有所延遲.對(duì)比可知,隨著環(huán)境氧體積分?jǐn)?shù)的增加,顆粒數(shù)增長(zhǎng)的速率峰值降低,聚并速率也降低.

    3.4?環(huán)境密度對(duì)碳煙生成歷程的影響

    環(huán)境壓力(密度)對(duì)燃燒和碳煙具有重要影響,發(fā)動(dòng)機(jī)增壓技術(shù)的應(yīng)用可以對(duì)進(jìn)氣壓力進(jìn)行有效調(diào)節(jié),因此有必要研究因環(huán)境壓力(密度)變化引起的碳煙發(fā)展變化趨勢(shì).考慮在環(huán)境溫度為1000K、O2體積分?jǐn)?shù)為15%條件下,將初始密度分別為14.8kg/m3和30kg/m3工況對(duì)應(yīng)的碳煙分布情況在-圖中用散點(diǎn)描述,如圖12所示.

    圖12?不同初始環(huán)境密度下碳煙分布區(qū)域的演化

    結(jié)合圖10(a)和圖12分析,隨著環(huán)境密度從14.8kg/m3增加到30kg/m3,在>2處,即對(duì)應(yīng)圖10(a)中碳煙成核、表面生長(zhǎng)區(qū)域,當(dāng)量比增大,其峰值從=6.3處增加到=8附近.當(dāng)量比的增加促進(jìn)了碳煙的生成.在<2區(qū)間,即對(duì)應(yīng)圖10(a)中OH氧化碳煙的主要分布區(qū)域,環(huán)境密度的增加并不改變?cè)搮^(qū)域的當(dāng)量比.

    考慮因環(huán)境密度改變引起的碳煙顆粒在數(shù)量上的演化趨勢(shì),圖13(a)、(b)中分別繪制了環(huán)境密度14.8kg/m3和30kg/m3條件下,碳煙總顆粒數(shù)增加速率曲線(藍(lán)線)與由于顆粒間聚并效應(yīng)引起的顆粒數(shù)減少速率曲線(紅線).對(duì)比不同密度下的顆粒數(shù)量變化速率,發(fā)現(xiàn)隨著環(huán)境密度的增大,碳煙開(kāi)始生成時(shí)間從噴油始點(diǎn)后約1ms提前到約0.7ms,顆粒增長(zhǎng)速率峰值從每秒約2.4×1016提高到4.6×1016,聚并速率則增加了一個(gè)數(shù)量級(jí).由于聚并速率的急劇增加,碳煙顆粒數(shù)增長(zhǎng)速率迅速達(dá)到峰值并趨于衰減.環(huán)境密度的增加在初期促進(jìn)了顆粒的成核過(guò)程,而在中期提高了碳煙顆粒的聚并速率.

    圖13 環(huán)境密度不同時(shí)顆??倲?shù)變化速率與顆粒間聚并速率

    4?結(jié)?論

    本文基于OpenFOAM開(kāi)源軟件平臺(tái)提供的噴霧燃燒求解器耦合半經(jīng)驗(yàn)多步現(xiàn)象學(xué)模型,研究了碳煙生成氧化各個(gè)階段在-關(guān)系圖中的分布特性,以及初始環(huán)境條件(環(huán)境溫度、氧氣體積分?jǐn)?shù)和環(huán)境壓力)對(duì)碳煙生成和演化歷程的影響.主要結(jié)論如下:

    (1) 依據(jù)半經(jīng)驗(yàn)?zāi)P停紵熡捎谄涑珊恕⒈砻嫔L(zhǎng)以及氧化過(guò)程產(chǎn)生的質(zhì)量變化速率,受關(guān)鍵參與組分的影響,分布在不同的當(dāng)量比-溫度區(qū)間.

    (2) 環(huán)境溫度的增大提高了碳煙生成區(qū)域(主要分布在>3區(qū)間)的整體當(dāng)量比,同時(shí)增加了碳煙受OH氧化的區(qū)域(主要分布在<2區(qū)間)的溫度.

    (3) 環(huán)境密度的增加使得碳煙生成區(qū)域(主要分布在>3區(qū)間)的整體當(dāng)量比提高,溫度不變.而碳煙受OH氧化的區(qū)域(主要分布在<2區(qū)間)基本不受影響.

    (4) 環(huán)境溫度和密度的增加對(duì)碳煙的成核以及聚并過(guò)程起促進(jìn)作用,且對(duì)聚并過(guò)程的影響遠(yuǎn)比成核作用更顯著.

    (5) 隨環(huán)境氧體積分?jǐn)?shù)的增加,總體碳煙顆粒數(shù)增長(zhǎng)速率和聚并速率降低.由于碳煙生長(zhǎng)區(qū)域整體當(dāng)量比增加,前驅(qū)物的體積分?jǐn)?shù)相應(yīng)增大,顆粒數(shù)密度峰值增大.

    [1] 陳夢(mèng)園,莊?遠(yuǎn),邱?亮,等. 柴油機(jī)碳煙顆粒三維形貌特征的研究[J]. 燃燒科學(xué)與技術(shù),2020,26(5):444-451.

    Chen Mengyuan,Zhuang Yuan,Qiu Liang,et al. 3D morphological characteristics of soot particles from diesel engine[J].,2020,26(5):444-451(in Chinese).

    [2] Kennedy I M. Models of soot formation and oxidation[J].,1997,23(2):95-132.

    [3] Hiroyasu H,Kadota T. Models for combustion and formation of nitric oxide and soot in direct injection diesel engines[J].,1976,6(1):327-335.

    [4] Frenklach M,Wang H. Detailed modeling of soot particle nucleation and growth[J].(),1991,23(1):1559-1566.

    [5] Frenklach M,Wang H. Detailed mechanism and modeling of soot particle formation[M]//. Springer,1994:165-192.

    [6] Appel J,Bockhorn H,F(xiàn)renklach M. Kinetic modeling of soot formation with detailed chemistry and physics:Laminar premixed flames of C2hydrocarbons[J].,2000,121(1/2):122-136.

    [7] Wang H,F(xiàn)renklach M. A detailed kinetic modeling study of aromatics formation in laminar premixed acetylene and ethylene flames[J].,1997,110(1/2):173-221.

    [8] Tao F,Golovitchev V I,Chomiak J. A phenomenologi-cal model for the prediction of soot formation in diesel spray combustion[J].,2004,136(3):270-282.

    [9] Tao Feng,Reitz Rolf D,F(xiàn)oster David E,et al. Nine-step phenomenological diesel soot model validated over a wide range of engine conditions[J].,2009,48(6):1223-1234.

    [10] Jia M,Peng Z J,Xie M Z. Numerical investigation of soot reduction potentials with diesel homogeneous charge compression ignition combustion by an improved phenomenological soot model[J].():,2009,223(3):395-412.

    [11] Leung K M,Lindstedt R P,Jones W P. A simplified reaction mechanism for soot formation in nonpremixed flames[J].,1991,87(3/4):289-305.

    [12] Vishwanathan Gokul,Reitz Rolf D. Development of a practical soot modeling approach and its application to low-temperature diesel combustion[J].,2010,182(8):1050-1082.

    [13] Wang Hu,Deneys Reitz Rolf,Yao Mingfa,et al. Development of an n-heptane-n-butanol-PAH mechanism and its application for combustion and soot prediction[J].

    ,2013,160(3):504-519.

    [14] Nagle J,Strickland-Constable R F. Oxidation of carbon between 1000-2000℃[C]//. Oxford,1962:154-164.

    [15] Neoh K G,Howard J B,Sarofim A F. Effect of oxidation on the physical structure of soot[J].(),1985,20(1):951-957.

    [16] Yoshihara Y,Kazakov A,Wang H,et al. Reduced mechanism of soot formation—Application to natural gas-fueled diesel combustion[J].(),1994,25(1):941-948.

    [17] Sandia National Laboratories. Engine Combustion Network(ECN):Experimental Diagnostics and Experimental Data[EB/OL].https://ecn.sandia.gov/diesel-spray-combustion/,2022.

    [18] Dec John E. A conceptual model of DI diesel combustion based on laser-sheet imaging[J].,1997:10.4271/970873.

    [19] Chishty Muhammad A,Bolla Michele,Hawkes Evatt R,et al. Soot formation modelling for-dodecane sprays using the transported PDF model[J].,2018,192:101-119.

    Numerical Simulation of Soot Evolution in Spray Combustion Under Different Ambient Environments

    Wang Chenchen,Zheng Zunqing,Wang Hu,Yao Mingfa,Zhang Yan

    (State Key Laboratory of Engines,Tianjin University,Tianjin 300072,China)

    Aimed at the soot formation and evolution under different diesel engine operation conditions,this study is carried out based on the measured data of soot generated in constant volume spray combustion by the engine combustion network(ECN). The combustion solver based on an open source CFD platform (OpenFOAM)is employed and the classical phenomenological multi-step soot model is coupled in the numerical simulation. The effect of ambient conditions on each soot sub-formation and oxidation process is revealed through the-diagram,and the distribution region and evolution of soot formation process are analyzed,with the aim of proposing the influencing mechanism of initial operation conditions on soot formation and oxidation. The results show that the phenomenological soot model can well predict the distribution of soot. The surface growth region of soot is mainly located in the area with the equivalence ratio greater than 3.0. The oxidation region is mainly distributed in the area with the equivalence ratio smaller than 2.0,which is also the main distribution area of OH species. The increase in ambient temperature and density promotes both nucleation and the coagulation rate and reduces the particle number density. The increase in ambient oxygen volume fraction suppress both nucleation and the coagulation rate,thus resulting in smaller size and higher density distribution of soot particles.

    soot;constant volume combustion;OpenFOAM;numerical simulation

    TK42

    A

    1006-8740(2022)04-0423-10

    10.11715/rskxjs.R202206009

    2022-01-07.

    國(guó)家自然科學(xué)基金資助項(xiàng)目(51976134;51876140).

    王晨辰(1990—??),男,博士研究生,1017201074@tju.edu.cn.

    鄭尊清,男,博士,教授,zhengzunqing@tju.edu.cn.

    (責(zé)任編輯:梁?霞)

    猜你喜歡
    當(dāng)量環(huán)境溫度體積
    多法并舉測(cè)量固體體積
    Review of a new bone tumor therapy strategy based on bifunctional biomaterials
    Bone Research(2021年2期)2021-09-11 06:02:56
    聚焦立體幾何中的體積問(wèn)題
    小體積帶來(lái)超高便攜性 Teufel Cinebar One
    誰(shuí)的體積大
    雷克薩斯CT200h車(chē)環(huán)境溫度顯示異常
    黃河之聲(2016年24期)2016-02-03 09:01:52
    超壓測(cè)試方法對(duì)炸藥TNT當(dāng)量計(jì)算結(jié)果的影響
    環(huán)空附加當(dāng)量循環(huán)密度的計(jì)算方法
    斷塊油氣田(2014年5期)2014-03-11 15:33:50
    環(huán)境溫度對(duì)連續(xù)剛構(gòu)橋模態(tài)頻率的影響
    国产一区亚洲一区在线观看| 国产视频首页在线观看| 国产精品三级大全| 我的女老师完整版在线观看| 极品人妻少妇av视频| 草草在线视频免费看| 成人手机av| 日本av免费视频播放| 男人操女人黄网站| 这个男人来自地球电影免费观看 | 久久午夜福利片| 人妻人人澡人人爽人人| 下体分泌物呈黄色| 丝袜在线中文字幕| 国产伦理片在线播放av一区| 国产精品一国产av| 国产综合精华液| 在线免费观看不下载黄p国产| 在线观看一区二区三区激情| 久久99一区二区三区| 日日摸夜夜添夜夜爱| 免费人妻精品一区二区三区视频| 激情五月婷婷亚洲| 最近2019中文字幕mv第一页| 久久毛片免费看一区二区三区| 国产成人精品久久久久久| 亚洲av.av天堂| 91精品国产国语对白视频| 激情视频va一区二区三区| 少妇被粗大的猛进出69影院 | 国产1区2区3区精品| 欧美丝袜亚洲另类| 亚洲精品美女久久av网站| 欧美亚洲 丝袜 人妻 在线| 各种免费的搞黄视频| 精品卡一卡二卡四卡免费| 亚洲国产成人一精品久久久| 成年人免费黄色播放视频| 欧美 日韩 精品 国产| 一级爰片在线观看| 99视频精品全部免费 在线| 草草在线视频免费看| 卡戴珊不雅视频在线播放| 亚洲精品中文字幕在线视频| 岛国毛片在线播放| 最新中文字幕久久久久| 亚洲欧美日韩卡通动漫| 亚洲av综合色区一区| 午夜福利,免费看| 汤姆久久久久久久影院中文字幕| 最近中文字幕高清免费大全6| 蜜桃国产av成人99| 亚洲,欧美,日韩| 久久99精品国语久久久| 七月丁香在线播放| 国产精品 国内视频| 亚洲精品色激情综合| 亚洲国产最新在线播放| 国产xxxxx性猛交| 午夜福利在线观看免费完整高清在| 色吧在线观看| 人妻少妇偷人精品九色| 久久亚洲国产成人精品v| 欧美日韩成人在线一区二区| 最近中文字幕高清免费大全6| 飞空精品影院首页| 大香蕉久久网| www.色视频.com| av国产久精品久网站免费入址| 久久精品国产a三级三级三级| 国产成人精品在线电影| 亚洲综合色惰| 热99国产精品久久久久久7| 国产色爽女视频免费观看| 亚洲第一区二区三区不卡| 免费不卡的大黄色大毛片视频在线观看| 一级毛片我不卡| 又黄又爽又刺激的免费视频.| 中文字幕另类日韩欧美亚洲嫩草| 九色成人免费人妻av| 国产精品久久久久成人av| 国产亚洲av片在线观看秒播厂| 精品国产露脸久久av麻豆| 七月丁香在线播放| 欧美xxxx性猛交bbbb| av视频免费观看在线观看| 亚洲在久久综合| 18禁国产床啪视频网站| 国产精品成人在线| 亚洲精品av麻豆狂野| 飞空精品影院首页| 中文字幕人妻丝袜制服| 夫妻午夜视频| 秋霞在线观看毛片| 另类亚洲欧美激情| 最黄视频免费看| 亚洲色图综合在线观看| 在线观看人妻少妇| 在线看a的网站| 亚洲综合色网址| 老司机影院成人| 日韩电影二区| 如日韩欧美国产精品一区二区三区| 综合色丁香网| 国产精品麻豆人妻色哟哟久久| 亚洲欧美一区二区三区黑人 | av不卡在线播放| 国产又爽黄色视频| 精品福利永久在线观看| 啦啦啦中文免费视频观看日本| 国产69精品久久久久777片| 在线精品无人区一区二区三| 国产精品秋霞免费鲁丝片| 午夜老司机福利剧场| 国产精品久久久av美女十八| 亚洲精品久久久久久婷婷小说| 欧美性感艳星| 国产一级毛片在线| 日本-黄色视频高清免费观看| 久久久国产一区二区| 最近2019中文字幕mv第一页| 欧美人与善性xxx| 婷婷色综合大香蕉| 国语对白做爰xxxⅹ性视频网站| 一级爰片在线观看| 99精国产麻豆久久婷婷| 亚洲欧美一区二区三区黑人 | 国产 精品1| 18禁国产床啪视频网站| 草草在线视频免费看| 欧美 亚洲 国产 日韩一| 男人添女人高潮全过程视频| 亚洲美女黄色视频免费看| 久久青草综合色| 久久久久久久国产电影| 日本黄大片高清| 国产片内射在线| 久久久精品区二区三区| 久久久久久久久久人人人人人人| 男女高潮啪啪啪动态图| 欧美成人精品欧美一级黄| 久久久久久久久久久免费av| 777米奇影视久久| 国产成人aa在线观看| 日韩在线高清观看一区二区三区| 日韩不卡一区二区三区视频在线| 国产老妇伦熟女老妇高清| 久久毛片免费看一区二区三区| 国产精品一国产av| 在现免费观看毛片| 成年人免费黄色播放视频| 日日撸夜夜添| 在线免费观看不下载黄p国产| 亚洲精品美女久久av网站| 久久久久久伊人网av| 国产精品一区www在线观看| 人人妻人人爽人人添夜夜欢视频| 日韩三级伦理在线观看| 国产69精品久久久久777片| 90打野战视频偷拍视频| √禁漫天堂资源中文www| 啦啦啦啦在线视频资源| 日本-黄色视频高清免费观看| 男女高潮啪啪啪动态图| 精品人妻一区二区三区麻豆| 精品人妻熟女毛片av久久网站| 人体艺术视频欧美日本| 久久热在线av| 夜夜爽夜夜爽视频| 观看av在线不卡| 婷婷色综合www| 亚洲精品美女久久久久99蜜臀 | 涩涩av久久男人的天堂| 精品人妻偷拍中文字幕| 中文字幕av电影在线播放| 亚洲伊人色综图| 国产成人精品无人区| 久久午夜综合久久蜜桃| 丁香六月天网| 亚洲美女黄色视频免费看| 欧美日韩国产mv在线观看视频| 成年动漫av网址| 国产精品成人在线| 成年人午夜在线观看视频| 日韩成人伦理影院| 国产精品人妻久久久影院| 在线免费观看不下载黄p国产| 免费看光身美女| av免费在线看不卡| 999精品在线视频| 中文字幕另类日韩欧美亚洲嫩草| 啦啦啦中文免费视频观看日本| 亚洲国产毛片av蜜桃av| 欧美人与性动交α欧美精品济南到 | 国产男女内射视频| 纯流量卡能插随身wifi吗| 免费观看av网站的网址| 亚洲,一卡二卡三卡| tube8黄色片| 国产黄频视频在线观看| 激情五月婷婷亚洲| 美女福利国产在线| 亚洲国产看品久久| 精品一区二区三区四区五区乱码 | 女人久久www免费人成看片| 在线观看三级黄色| 乱人伦中国视频| 一级,二级,三级黄色视频| 黄色毛片三级朝国网站| 国产精品国产av在线观看| 又大又黄又爽视频免费| 欧美精品亚洲一区二区| 在线天堂中文资源库| 国产男人的电影天堂91| av线在线观看网站| 欧美亚洲 丝袜 人妻 在线| 精品久久久精品久久久| 成年人免费黄色播放视频| 中文精品一卡2卡3卡4更新| 丰满饥渴人妻一区二区三| www.av在线官网国产| 亚洲av福利一区| 国产成人aa在线观看| 夜夜爽夜夜爽视频| 精品国产乱码久久久久久小说| 少妇猛男粗大的猛烈进出视频| 飞空精品影院首页| 欧美日韩成人在线一区二区| 热99国产精品久久久久久7| 国产在线视频一区二区| 国产男女内射视频| 又黄又爽又刺激的免费视频.| 午夜av观看不卡| 国语对白做爰xxxⅹ性视频网站| 国产毛片在线视频| 一级a做视频免费观看| av线在线观看网站| 国产一区二区在线观看日韩| a级毛色黄片| 精品熟女少妇av免费看| 亚洲精品久久午夜乱码| 一区二区日韩欧美中文字幕 | 国产无遮挡羞羞视频在线观看| 精品熟女少妇av免费看| 国产免费一区二区三区四区乱码| 久久午夜福利片| 成人午夜精彩视频在线观看| 成人国产av品久久久| 亚洲,一卡二卡三卡| 成年女人在线观看亚洲视频| 国产一区二区在线观看av| 久久99热6这里只有精品| 熟女电影av网| 精品第一国产精品| 日韩精品免费视频一区二区三区 | 免费看av在线观看网站| 1024视频免费在线观看| 五月天丁香电影| 看十八女毛片水多多多| 欧美精品av麻豆av| 日韩欧美一区视频在线观看| 精品国产一区二区久久| 一级黄片播放器| 99视频精品全部免费 在线| 日韩成人av中文字幕在线观看| 青春草视频在线免费观看| 日韩不卡一区二区三区视频在线| 国产精品久久久久久久久免| 蜜臀久久99精品久久宅男| 精品久久久精品久久久| 亚洲av电影在线进入| 中国国产av一级| 青春草视频在线免费观看| 国产麻豆69| 黄网站色视频无遮挡免费观看| 男女免费视频国产| 毛片一级片免费看久久久久| 久久久久精品人妻al黑| 在线观看国产h片| 又黄又爽又刺激的免费视频.| 欧美bdsm另类| 香蕉精品网在线| 久久久久久久久久久久大奶| 青青草视频在线视频观看| av不卡在线播放| 少妇高潮的动态图| xxxhd国产人妻xxx| 午夜视频国产福利| 久久久国产一区二区| 97在线人人人人妻| 日韩精品有码人妻一区| 黄色 视频免费看| 午夜日本视频在线| 欧美精品亚洲一区二区| 交换朋友夫妻互换小说| 日本欧美国产在线视频| 国产有黄有色有爽视频| videossex国产| 欧美日本中文国产一区发布| 青青草视频在线视频观看| 成人午夜精彩视频在线观看| 亚洲欧洲日产国产| 精品国产一区二区三区四区第35| 国产一区亚洲一区在线观看| 国产精品麻豆人妻色哟哟久久| 免费播放大片免费观看视频在线观看| 午夜精品国产一区二区电影| 一区在线观看完整版| 国产日韩欧美视频二区| 亚洲精品自拍成人| 精品一区二区三区视频在线| 欧美xxⅹ黑人| www日本在线高清视频| 夫妻性生交免费视频一级片| 在线观看免费高清a一片| 尾随美女入室| 在线观看国产h片| 国产av一区二区精品久久| 777米奇影视久久| 乱码一卡2卡4卡精品| 亚洲美女视频黄频| av线在线观看网站| 精品一区二区免费观看| av线在线观看网站| 成人午夜精彩视频在线观看| 国产亚洲午夜精品一区二区久久| av在线老鸭窝| 草草在线视频免费看| 中文字幕亚洲精品专区| 婷婷色综合大香蕉| 女的被弄到高潮叫床怎么办| 一区在线观看完整版| 国产亚洲一区二区精品| 91aial.com中文字幕在线观看| 最近2019中文字幕mv第一页| 国产黄色免费在线视频| 青春草国产在线视频| 日韩一本色道免费dvd| 国产老妇伦熟女老妇高清| 免费女性裸体啪啪无遮挡网站| 欧美日韩国产mv在线观看视频| 国产白丝娇喘喷水9色精品| 免费观看性生交大片5| 深夜精品福利| 青青草视频在线视频观看| 国产伦理片在线播放av一区| 综合色丁香网| 九草在线视频观看| 久久97久久精品| 十八禁网站网址无遮挡| 日本av手机在线免费观看| 中国美白少妇内射xxxbb| 一级毛片电影观看| 亚洲精品一二三| 一级毛片电影观看| 熟女av电影| 久久午夜综合久久蜜桃| 亚洲,欧美,日韩| 亚洲情色 制服丝袜| 中国美白少妇内射xxxbb| 国产精品一区二区在线不卡| 亚洲精华国产精华液的使用体验| 9色porny在线观看| 国产一区二区三区综合在线观看 | 国产精品一区二区在线观看99| 久久精品久久精品一区二区三区| www.av在线官网国产| 日日爽夜夜爽网站| 久久女婷五月综合色啪小说| 亚洲综合精品二区| av.在线天堂| 大香蕉97超碰在线| 男女免费视频国产| 九色亚洲精品在线播放| 1024视频免费在线观看| 丰满饥渴人妻一区二区三| 午夜免费男女啪啪视频观看| 亚洲国产精品专区欧美| 亚洲图色成人| 一级毛片 在线播放| 日本91视频免费播放| 99热网站在线观看| 亚洲国产毛片av蜜桃av| 国产精品免费大片| 精品一品国产午夜福利视频| 免费av中文字幕在线| 国产成人午夜福利电影在线观看| 青青草视频在线视频观看| 91精品三级在线观看| 天堂8中文在线网| 久久久久精品人妻al黑| 嫩草影院入口| 日本-黄色视频高清免费观看| 99久国产av精品国产电影| 中文字幕另类日韩欧美亚洲嫩草| 男人爽女人下面视频在线观看| 国产免费福利视频在线观看| 大香蕉97超碰在线| 久久女婷五月综合色啪小说| 久久97久久精品| 看十八女毛片水多多多| 亚洲色图综合在线观看| 国产伦理片在线播放av一区| 久久热在线av| 国产一区二区激情短视频 | 内地一区二区视频在线| 少妇人妻久久综合中文| 亚洲国产毛片av蜜桃av| 多毛熟女@视频| 国产精品国产av在线观看| 成人18禁高潮啪啪吃奶动态图| 欧美+日韩+精品| 亚洲欧美日韩卡通动漫| 熟女电影av网| 熟女人妻精品中文字幕| 91国产中文字幕| 欧美 日韩 精品 国产| 少妇人妻 视频| 亚洲精品一区蜜桃| 五月玫瑰六月丁香| 80岁老熟妇乱子伦牲交| 日本av手机在线免费观看| av有码第一页| 高清欧美精品videossex| 亚洲精品aⅴ在线观看| 国产成人aa在线观看| 大香蕉久久成人网| 国产精品蜜桃在线观看| 菩萨蛮人人尽说江南好唐韦庄| 97在线人人人人妻| 午夜激情久久久久久久| 熟女av电影| 久久99热这里只频精品6学生| 国产白丝娇喘喷水9色精品| 午夜老司机福利剧场| 大话2 男鬼变身卡| 人妻系列 视频| 婷婷色综合大香蕉| 久久久国产欧美日韩av| 国产男女内射视频| 在线精品无人区一区二区三| 9191精品国产免费久久| 国产免费福利视频在线观看| 国产日韩欧美视频二区| 久久久久精品久久久久真实原创| 国产精品久久久av美女十八| 国产亚洲最大av| 高清视频免费观看一区二区| 人妻少妇偷人精品九色| 久久99精品国语久久久| 内地一区二区视频在线| 亚洲综合色网址| 免费看不卡的av| 久久精品久久精品一区二区三区| 精品国产露脸久久av麻豆| 欧美日韩av久久| 久久99蜜桃精品久久| 成年av动漫网址| 爱豆传媒免费全集在线观看| 国产精品人妻久久久久久| 综合色丁香网| 熟女电影av网| 99热6这里只有精品| 观看av在线不卡| 精品人妻一区二区三区麻豆| 夫妻午夜视频| www.av在线官网国产| 日韩免费高清中文字幕av| 一级片免费观看大全| 另类精品久久| 亚洲综合色惰| av视频免费观看在线观看| 亚洲精品久久午夜乱码| 日本免费在线观看一区| 亚洲人与动物交配视频| 亚洲国产看品久久| 亚洲美女黄色视频免费看| 高清毛片免费看| 国产精品人妻久久久影院| 国产成人午夜福利电影在线观看| 亚洲精品av麻豆狂野| 免费观看av网站的网址| 99热这里只有是精品在线观看| 国产国语露脸激情在线看| 国产成人精品一,二区| 欧美另类一区| 亚洲中文av在线| 国产亚洲精品久久久com| 婷婷色麻豆天堂久久| 久久久久久人人人人人| 各种免费的搞黄视频| 97超碰精品成人国产| av在线播放精品| 99视频精品全部免费 在线| 最近2019中文字幕mv第一页| 久久久久精品久久久久真实原创| 美女中出高潮动态图| 欧美日韩视频高清一区二区三区二| 久久久久久久久久久免费av| 人妻少妇偷人精品九色| 免费在线观看黄色视频的| 国产色婷婷99| 色视频在线一区二区三区| 国产熟女欧美一区二区| 亚洲五月色婷婷综合| 中文字幕人妻熟女乱码| 亚洲三级黄色毛片| 天美传媒精品一区二区| 日本-黄色视频高清免费观看| 麻豆精品久久久久久蜜桃| 亚洲欧美成人综合另类久久久| 在线观看一区二区三区激情| 欧美日韩视频高清一区二区三区二| 欧美日韩综合久久久久久| 午夜激情av网站| 中文字幕人妻熟女乱码| 青春草视频在线免费观看| 人人妻人人添人人爽欧美一区卜| 久久人妻熟女aⅴ| 建设人人有责人人尽责人人享有的| 少妇人妻久久综合中文| kizo精华| 一区二区三区乱码不卡18| 亚洲激情五月婷婷啪啪| 久久精品国产a三级三级三级| 国国产精品蜜臀av免费| 五月开心婷婷网| 曰老女人黄片| 国产高清国产精品国产三级| 亚洲成人av在线免费| 下体分泌物呈黄色| 亚洲欧洲国产日韩| 在线观看免费视频网站a站| 亚洲欧美日韩卡通动漫| 久久久久久久久久久久大奶| 久久人人97超碰香蕉20202| 午夜激情久久久久久久| 欧美激情极品国产一区二区三区 | 亚洲av国产av综合av卡| 蜜桃国产av成人99| 又粗又硬又长又爽又黄的视频| 国产一区二区三区综合在线观看 | 九九爱精品视频在线观看| 日韩av免费高清视频| 久久精品国产亚洲av涩爱| 国产精品嫩草影院av在线观看| 欧美精品一区二区免费开放| av在线老鸭窝| 亚洲精品久久午夜乱码| 亚洲精品aⅴ在线观看| 国产欧美日韩一区二区三区在线| 久久人人爽人人爽人人片va| 国产一区二区激情短视频 | 两个人免费观看高清视频| 男男h啪啪无遮挡| 视频在线观看一区二区三区| 国产乱人偷精品视频| 大码成人一级视频| 国产欧美亚洲国产| 久久久久精品性色| 人人妻人人澡人人看| 在线观看免费视频网站a站| 毛片一级片免费看久久久久| 最新的欧美精品一区二区| 乱码一卡2卡4卡精品| 国产成人av激情在线播放| 男女免费视频国产| 男人舔女人的私密视频| 国产精品欧美亚洲77777| 国产高清不卡午夜福利| 国产亚洲一区二区精品| 777米奇影视久久| 亚洲中文av在线| 亚洲五月色婷婷综合| 中国三级夫妇交换| 日韩,欧美,国产一区二区三区| 五月天丁香电影| 日韩不卡一区二区三区视频在线| 亚洲精华国产精华液的使用体验| 欧美人与性动交α欧美软件 | 国产又爽黄色视频| 黄色毛片三级朝国网站| 午夜日本视频在线| 久久久a久久爽久久v久久| 免费观看性生交大片5| 精品人妻偷拍中文字幕| 极品少妇高潮喷水抽搐| 午夜免费鲁丝| 欧美精品一区二区大全| 亚洲美女搞黄在线观看| 久久ye,这里只有精品| 久久久久久久久久人人人人人人| 另类精品久久| 午夜久久久在线观看| 大陆偷拍与自拍| 亚洲精品中文字幕在线视频| 国产亚洲午夜精品一区二区久久| 黑丝袜美女国产一区| 亚洲精品中文字幕在线视频| 国产精品三级大全| 插逼视频在线观看| 美女视频免费永久观看网站| 久久久国产一区二区| 国产男女超爽视频在线观看| 青春草国产在线视频| 久久人人爽人人爽人人片va| 亚洲人成77777在线视频| 韩国高清视频一区二区三区| 高清av免费在线| 国产av精品麻豆| 婷婷色av中文字幕| 色哟哟·www| 国产精品久久久久久久电影| 中文字幕制服av| 亚洲欧洲精品一区二区精品久久久 |