• <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)頻率的影響
    最近手机中文字幕大全| 五月开心婷婷网| 久久婷婷青草| 亚洲美女搞黄在线观看| 色吧在线观看| 免费观看a级毛片全部| 国产精品偷伦视频观看了| 在线看a的网站| 大香蕉久久网| 亚洲国产欧美在线一区| 国产精品一区二区在线观看99| 久久久久精品性色| 如何舔出高潮| 最近的中文字幕免费完整| 欧美 亚洲 国产 日韩一| 9色porny在线观看| 啦啦啦在线观看免费高清www| www.色视频.com| 日本vs欧美在线观看视频| 最新中文字幕久久久久| 五月伊人婷婷丁香| 亚洲色图 男人天堂 中文字幕 | 免费人妻精品一区二区三区视频| 美女大奶头黄色视频| 国产av精品麻豆| 人人妻人人澡人人爽人人夜夜| 国产成人一区二区在线| 日韩成人伦理影院| 人妻 亚洲 视频| 亚洲欧美成人精品一区二区| 在线观看免费视频网站a站| 在线观看人妻少妇| 国产精品麻豆人妻色哟哟久久| 日韩制服丝袜自拍偷拍| 18禁动态无遮挡网站| 天堂俺去俺来也www色官网| 人妻少妇偷人精品九色| 咕卡用的链子| 国产成人a∨麻豆精品| 久久午夜综合久久蜜桃| 久久午夜综合久久蜜桃| 91久久精品国产一区二区三区| √禁漫天堂资源中文www| 少妇的丰满在线观看| 七月丁香在线播放| 九九在线视频观看精品| 精品99又大又爽又粗少妇毛片| 最近2019中文字幕mv第一页| 久久国内精品自在自线图片| 天美传媒精品一区二区| 国产熟女午夜一区二区三区| 超色免费av| 国产极品天堂在线| 男女免费视频国产| 寂寞人妻少妇视频99o| 中文字幕人妻丝袜制服| 国产精品久久久久久久电影| 大陆偷拍与自拍| 最近中文字幕2019免费版| 免费大片黄手机在线观看| 久久久国产精品麻豆| 夫妻性生交免费视频一级片| 99国产综合亚洲精品| 日本黄大片高清| 90打野战视频偷拍视频| 国产又色又爽无遮挡免| 精品人妻在线不人妻| 午夜免费观看性视频| 日韩欧美一区视频在线观看| 久久这里有精品视频免费| 久久精品aⅴ一区二区三区四区 | 美女国产高潮福利片在线看| 国产一区有黄有色的免费视频| tube8黄色片| 乱人伦中国视频| 一个人免费看片子| 午夜福利视频在线观看免费| 涩涩av久久男人的天堂| 老女人水多毛片| 天天影视国产精品| 巨乳人妻的诱惑在线观看| 久久毛片免费看一区二区三区| 成人午夜精彩视频在线观看| 两个人免费观看高清视频| 免费看不卡的av| 女性生殖器流出的白浆| 寂寞人妻少妇视频99o| a级片在线免费高清观看视频| 国产成人午夜福利电影在线观看| 国产精品国产av在线观看| 日韩成人伦理影院| 成年人免费黄色播放视频| 九九在线视频观看精品| 国产精品久久久久久久电影| 亚洲精品aⅴ在线观看| 亚洲一码二码三码区别大吗| 午夜日本视频在线| 国产极品粉嫩免费观看在线| 国产精品国产三级国产av玫瑰| 欧美另类一区| 亚洲欧洲精品一区二区精品久久久 | 亚洲精品一区蜜桃| 极品少妇高潮喷水抽搐| 精品一区二区三卡| 欧美激情 高清一区二区三区| 蜜桃在线观看..| 18禁国产床啪视频网站| 免费看光身美女| 久久久国产精品麻豆| 国产精品久久久久久久电影| 久久久久久久久久久久大奶| 精品人妻偷拍中文字幕| 国产黄色免费在线视频| 成人黄色视频免费在线看| 人人妻人人澡人人看| 各种免费的搞黄视频| 精品一区在线观看国产| 在线观看国产h片| 激情五月婷婷亚洲| 久久久久人妻精品一区果冻| 中国国产av一级| 欧美日韩视频精品一区| 99热国产这里只有精品6| 自线自在国产av| av国产久精品久网站免费入址| 中国国产av一级| 伦理电影免费视频| 日韩中文字幕视频在线看片| 黄色配什么色好看| 97超碰精品成人国产| 热re99久久精品国产66热6| 国产爽快片一区二区三区| 男女边吃奶边做爰视频| 草草在线视频免费看| 欧美日韩精品成人综合77777| 国产永久视频网站| 日韩一本色道免费dvd| 成人影院久久| 熟妇人妻不卡中文字幕| 亚洲av日韩在线播放| 日本猛色少妇xxxxx猛交久久| 少妇熟女欧美另类| 亚洲精品美女久久av网站| 婷婷色综合大香蕉| 2021少妇久久久久久久久久久| 乱码一卡2卡4卡精品| 男女免费视频国产| 久久久久人妻精品一区果冻| 国产乱人偷精品视频| 熟女电影av网| 亚洲精品国产色婷婷电影| 久久久精品区二区三区| 最近中文字幕2019免费版| 自拍欧美九色日韩亚洲蝌蚪91| 国产精品99久久99久久久不卡 | 中文字幕人妻丝袜制服| 中文字幕最新亚洲高清| 久久精品国产a三级三级三级| 国产精品一区二区在线观看99| 色吧在线观看| 国产精品国产三级国产专区5o| 国产欧美亚洲国产| 又粗又硬又长又爽又黄的视频| 国产精品久久久久久精品电影小说| 日韩成人伦理影院| 亚洲av.av天堂| 久久久久网色| 99re6热这里在线精品视频| 精品久久久久久电影网| 最近的中文字幕免费完整| 色5月婷婷丁香| 69精品国产乱码久久久| 老司机亚洲免费影院| 少妇的逼水好多| 制服诱惑二区| 国产1区2区3区精品| 纯流量卡能插随身wifi吗| 亚洲第一av免费看| 考比视频在线观看| 欧美日韩综合久久久久久| 日韩在线高清观看一区二区三区| 免费人成在线观看视频色| 国产精品.久久久| 性高湖久久久久久久久免费观看| 中文天堂在线官网| 只有这里有精品99| 久久久久久久久久人人人人人人| 亚洲精品美女久久av网站| 国产精品.久久久| 成人毛片a级毛片在线播放| 国产av精品麻豆| 22中文网久久字幕| 国产成人免费无遮挡视频| 在线免费观看不下载黄p国产| 日韩一区二区视频免费看| 国产永久视频网站| 亚洲欧洲精品一区二区精品久久久 | 99九九在线精品视频| 一级爰片在线观看| 日本午夜av视频| 婷婷色综合www| av视频免费观看在线观看| 日本午夜av视频| 99久国产av精品国产电影| 97超碰精品成人国产| 国产精品秋霞免费鲁丝片| 亚洲欧美中文字幕日韩二区| 女性生殖器流出的白浆| 久久热在线av| 欧美97在线视频| 如日韩欧美国产精品一区二区三区| 最后的刺客免费高清国语| 考比视频在线观看| 成人国语在线视频| 亚洲国产成人一精品久久久| 汤姆久久久久久久影院中文字幕| 又黄又粗又硬又大视频| 91精品三级在线观看| 最近最新中文字幕免费大全7| 一级黄片播放器| 黑丝袜美女国产一区| 久久久精品94久久精品| 国产乱人偷精品视频| 日产精品乱码卡一卡2卡三| 好男人视频免费观看在线| 人妻 亚洲 视频| 久久久久精品久久久久真实原创| 熟女人妻精品中文字幕| 99精国产麻豆久久婷婷| 看免费av毛片| 午夜免费男女啪啪视频观看| 日韩三级伦理在线观看| 午夜精品国产一区二区电影| 欧美亚洲 丝袜 人妻 在线| 99九九在线精品视频| 欧美激情 高清一区二区三区| 女人精品久久久久毛片| 熟妇人妻不卡中文字幕| 国产 一区精品| 欧美国产精品va在线观看不卡| 久久鲁丝午夜福利片| 777米奇影视久久| 国产成人欧美| freevideosex欧美| 中文乱码字字幕精品一区二区三区| 捣出白浆h1v1| 搡老乐熟女国产| 国产乱人偷精品视频| 久久99蜜桃精品久久| 五月天丁香电影| 国产在线免费精品| 亚洲情色 制服丝袜| 国产一区有黄有色的免费视频| 欧美日韩精品成人综合77777| 水蜜桃什么品种好| 亚洲一码二码三码区别大吗| 日韩欧美精品免费久久| 黄色一级大片看看| 黑人高潮一二区| 老司机影院成人| 亚洲高清免费不卡视频| 久久久精品94久久精品| 国产又爽黄色视频| 国产一区二区在线观看日韩| 少妇人妻精品综合一区二区| 欧美日韩综合久久久久久| 人人澡人人妻人| 秋霞伦理黄片| 中文精品一卡2卡3卡4更新| 亚洲精品中文字幕在线视频| 免费人妻精品一区二区三区视频| 午夜激情av网站| 久久国产精品大桥未久av| 色视频在线一区二区三区| 久久这里有精品视频免费| 极品人妻少妇av视频| 国产亚洲一区二区精品| 一级,二级,三级黄色视频| 午夜老司机福利剧场| 在线天堂中文资源库| 免费av中文字幕在线| tube8黄色片| 久久久久精品久久久久真实原创| 欧美成人精品欧美一级黄| 天天躁夜夜躁狠狠久久av| 亚洲欧美清纯卡通| 国产高清不卡午夜福利| 久久国产亚洲av麻豆专区| 18+在线观看网站| 日韩中文字幕视频在线看片| 丝袜喷水一区| 久久久久久伊人网av| 国产一级毛片在线| 久久99蜜桃精品久久| 国产精品一区www在线观看| 中文字幕另类日韩欧美亚洲嫩草| 精品少妇内射三级| 九草在线视频观看| 亚洲欧美日韩卡通动漫| 国产探花极品一区二区| 美女xxoo啪啪120秒动态图| 久久精品熟女亚洲av麻豆精品| 精品久久久久久电影网| 亚洲精品一区蜜桃| 人人妻人人澡人人看| 久久久精品区二区三区| 国产精品成人在线| 国产精品国产三级国产专区5o| 亚洲精品aⅴ在线观看| 午夜91福利影院| 亚洲久久久国产精品| 国产精品三级大全| 日韩电影二区| 精品午夜福利在线看| 免费黄网站久久成人精品| 国产精品成人在线| 超碰97精品在线观看| 久久免费观看电影| a 毛片基地| av.在线天堂| 久久精品国产综合久久久 | 久久午夜综合久久蜜桃| 亚洲精品一二三| 高清av免费在线| 国产一区有黄有色的免费视频| 久久这里只有精品19| 国产精品久久久久久久久免| 欧美日韩视频精品一区| 欧美成人精品欧美一级黄| 中文欧美无线码| 久久人人爽人人片av| 纵有疾风起免费观看全集完整版| 女性被躁到高潮视频| 九九在线视频观看精品| 亚洲成人手机| 国产成人精品在线电影| 91精品伊人久久大香线蕉| 亚洲精品一二三| www.色视频.com| 一边摸一边做爽爽视频免费| 人人澡人人妻人| 亚洲伊人色综图| 精品人妻偷拍中文字幕| 中国国产av一级| 男女午夜视频在线观看 | 久久韩国三级中文字幕| 欧美成人精品欧美一级黄| 久久人妻熟女aⅴ| 黑人猛操日本美女一级片| 女的被弄到高潮叫床怎么办| 久久久国产精品麻豆| 一级片'在线观看视频| 黄色视频在线播放观看不卡| 熟妇人妻不卡中文字幕| 午夜免费观看性视频| 高清黄色对白视频在线免费看| 美女主播在线视频| 国产精品熟女久久久久浪| 国产av精品麻豆| 欧美日韩国产mv在线观看视频| 伦理电影大哥的女人| 久久青草综合色| 国产在线视频一区二区| 国精品久久久久久国模美| 99re6热这里在线精品视频| 国产成人av激情在线播放| 欧美人与善性xxx| 丰满迷人的少妇在线观看| 熟女av电影| 欧美性感艳星| 精品少妇黑人巨大在线播放| 啦啦啦中文免费视频观看日本| 99国产精品免费福利视频| 成人18禁高潮啪啪吃奶动态图| 热99久久久久精品小说推荐| 在线观看国产h片| 国产精品人妻久久久久久| www.熟女人妻精品国产 | 黄色怎么调成土黄色| 国产毛片在线视频| 黄色怎么调成土黄色| 亚洲精品美女久久av网站| 国精品久久久久久国模美| 女人久久www免费人成看片| 日日爽夜夜爽网站| 国产免费现黄频在线看| 一级毛片我不卡| 99精国产麻豆久久婷婷| 女性被躁到高潮视频| 交换朋友夫妻互换小说| 国产乱来视频区| 免费高清在线观看视频在线观看| 久久精品久久久久久久性| 国产精品偷伦视频观看了| 亚洲精品中文字幕在线视频| 搡女人真爽免费视频火全软件| 高清毛片免费看| 久久这里只有精品19| 91成人精品电影| 精品午夜福利在线看| 久久精品熟女亚洲av麻豆精品| 黑人猛操日本美女一级片| 97超碰精品成人国产| 色视频在线一区二区三区| 久久久久国产网址| 久久午夜福利片| 中文字幕人妻丝袜制服| 国产精品99久久99久久久不卡 | 制服丝袜香蕉在线| 久久久久国产精品人妻一区二区| 一区二区三区精品91| 另类精品久久| 日韩精品有码人妻一区| 精品亚洲成a人片在线观看| av又黄又爽大尺度在线免费看| 日本欧美视频一区| 国产精品欧美亚洲77777| 亚洲人成网站在线观看播放| 久久综合国产亚洲精品| 亚洲经典国产精华液单| 春色校园在线视频观看| 母亲3免费完整高清在线观看| 一边摸一边做爽爽视频免费| 中文字幕人妻熟女乱码| 十八禁高潮呻吟视频| 国产麻豆69| 国产高清视频在线播放一区| 国产精品成人在线| 黑人欧美特级aaaaaa片| 久久久国产成人精品二区 | 国产亚洲av高清不卡| 黄频高清免费视频| 精品国内亚洲2022精品成人 | 久久九九热精品免费| 曰老女人黄片| 国产精品一区二区在线观看99| 国产aⅴ精品一区二区三区波| 欧美精品高潮呻吟av久久| 极品人妻少妇av视频| 每晚都被弄得嗷嗷叫到高潮| 国产精品成人在线| 国产视频一区二区在线看| 国产精品久久视频播放| 成年人午夜在线观看视频| 免费观看人在逋| 午夜福利影视在线免费观看| 中文字幕人妻丝袜制服| 亚洲专区字幕在线| 黄色视频,在线免费观看| 天天躁日日躁夜夜躁夜夜| 久久中文字幕一级| 超碰成人久久| 男女之事视频高清在线观看| 久久精品国产综合久久久| 免费高清在线观看日韩| 99久久99久久久精品蜜桃| 久久国产精品大桥未久av| 久久国产亚洲av麻豆专区| 成人免费观看视频高清| 午夜免费鲁丝| 午夜福利,免费看| 男女高潮啪啪啪动态图| 日韩欧美一区视频在线观看| 亚洲avbb在线观看| 色在线成人网| 两人在一起打扑克的视频| 亚洲专区中文字幕在线| 99热国产这里只有精品6| 嫩草影视91久久| av中文乱码字幕在线| 久久性视频一级片| 十八禁高潮呻吟视频| 国产激情久久老熟女| a级片在线免费高清观看视频| 看免费av毛片| 国产伦人伦偷精品视频| 女同久久另类99精品国产91| 三级毛片av免费| 国产真人三级小视频在线观看| 黄色丝袜av网址大全| 岛国毛片在线播放| 老司机深夜福利视频在线观看| 一边摸一边抽搐一进一小说 | 免费不卡黄色视频| 午夜福利乱码中文字幕| 国产成人av激情在线播放| 亚洲性夜色夜夜综合| 国产精品国产高清国产av | 色在线成人网| 久久午夜综合久久蜜桃| 天天添夜夜摸| 女同久久另类99精品国产91| 国产麻豆69| 午夜福利乱码中文字幕| 日日摸夜夜添夜夜添小说| 欧美激情高清一区二区三区| 精品少妇一区二区三区视频日本电影| 人人妻人人澡人人看| 日韩精品免费视频一区二区三区| 91老司机精品| 黄频高清免费视频| 热re99久久精品国产66热6| 亚洲成a人片在线一区二区| 久久精品aⅴ一区二区三区四区| 看免费av毛片| 满18在线观看网站| 欧美日韩黄片免| 亚洲一区高清亚洲精品| 免费在线观看完整版高清| 咕卡用的链子| 桃红色精品国产亚洲av| 99国产极品粉嫩在线观看| 欧美中文综合在线视频| 国产亚洲欧美在线一区二区| 丝袜美足系列| 亚洲av电影在线进入| 精品国产国语对白av| 日韩欧美免费精品| 岛国在线观看网站| 久久久久久久精品吃奶| 成年版毛片免费区| 18禁黄网站禁片午夜丰满| 在线观看免费视频网站a站| 国产精品电影一区二区三区 | 天堂俺去俺来也www色官网| av片东京热男人的天堂| 欧美午夜高清在线| 99久久99久久久精品蜜桃| 中国美女看黄片| 日韩欧美国产一区二区入口| 久久久久精品国产欧美久久久| 国产精品.久久久| 97人妻天天添夜夜摸| 18禁黄网站禁片午夜丰满| 一区二区日韩欧美中文字幕| 亚洲中文日韩欧美视频| 啪啪无遮挡十八禁网站| 日韩制服丝袜自拍偷拍| 色尼玛亚洲综合影院| 露出奶头的视频| 亚洲精品一卡2卡三卡4卡5卡| 日本wwww免费看| 日韩欧美一区视频在线观看| 亚洲成人手机| 欧美亚洲 丝袜 人妻 在线| 91成人精品电影| 黄片播放在线免费| 免费观看人在逋| 欧美亚洲日本最大视频资源| av网站在线播放免费| 日本黄色日本黄色录像| 亚洲国产精品合色在线| bbb黄色大片| 国产精品国产av在线观看| 韩国av一区二区三区四区| 国产精华一区二区三区| 午夜福利视频在线观看免费| 国产在线精品亚洲第一网站| 狠狠狠狠99中文字幕| 日韩免费av在线播放| 国产精品永久免费网站| 国产精品一区二区在线观看99| 女同久久另类99精品国产91| 十分钟在线观看高清视频www| 久久中文看片网| 欧美性长视频在线观看| 免费在线观看黄色视频的| 国产精品久久久人人做人人爽| 在线观看舔阴道视频| 成人手机av| 午夜福利欧美成人| 国产成人免费无遮挡视频| av网站免费在线观看视频| 国产在线一区二区三区精| 黑人欧美特级aaaaaa片| 久久99一区二区三区| 一进一出抽搐gif免费好疼 | 黄网站色视频无遮挡免费观看| 国产成人欧美在线观看 | 中文字幕精品免费在线观看视频| 国产熟女午夜一区二区三区| 久久精品国产亚洲av香蕉五月 | 最新的欧美精品一区二区| 在线视频色国产色| 看免费av毛片| 最近最新免费中文字幕在线| 亚洲成人免费av在线播放| 欧美最黄视频在线播放免费 | 热99re8久久精品国产| 另类亚洲欧美激情| 亚洲自偷自拍图片 自拍| 操出白浆在线播放| 国产91精品成人一区二区三区| 久久精品亚洲精品国产色婷小说| 欧美人与性动交α欧美精品济南到| 久久精品熟女亚洲av麻豆精品| netflix在线观看网站| 如日韩欧美国产精品一区二区三区| avwww免费| 天天躁日日躁夜夜躁夜夜| 麻豆国产av国片精品| 在线永久观看黄色视频| av不卡在线播放| 午夜老司机福利片| 欧美黑人欧美精品刺激| 美女国产高潮福利片在线看| 精品国产一区二区久久| 99精品欧美一区二区三区四区| 97人妻天天添夜夜摸| 国产精品亚洲av一区麻豆| 国产一区二区三区综合在线观看| 97人妻天天添夜夜摸| 首页视频小说图片口味搜索| 五月开心婷婷网| 久久久精品区二区三区|