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

    反應性雙射流中標量輸運和化學反應特性

    2022-05-26 02:57:24李巖田阿慧周毅
    化工學報 2022年5期
    關鍵詞:標量反應物流向

    李巖,田阿慧,周毅

    (南京理工大學能源與動力工程學院,江蘇南京 210094)

    引 言

    雙射流和多束射流廣泛存在于工業(yè)領域,如鍋爐燃燒器、排氣煙囪、發(fā)動機燃燒室等。雙射流作為多束射流研究的基礎模型是眾多學者的研究對象。與均勻各向同性湍流相比,雙射流的時空演化機理更為復雜,涉及諸多湍流過程(界面層卷吸、噪聲的產(chǎn)生、射流擬序結構演化以及射流模態(tài)等)[1]。雙射流流動往往伴隨著化學反應,流動過程影響著標量的輸運,而化學反應則決定了標量的產(chǎn)生及耗散。故研究雙射流間相互作用和雙射流間化學反應中的標量特性具有重要意義。

    目前國內(nèi)外學者在雙射流流場統(tǒng)計特性[2-4]和射流間相互作用[4-7]等領域開展了一系列實驗和數(shù)值模擬研究。在實驗方面,Miller 等[6]、Tanaka[7]開展了關于雙射流流場特性的開創(chuàng)性研究,將雙射流流場結構分為收斂區(qū)、混合區(qū)和合并區(qū)三個區(qū)域。此后,學者多分析不同射流口間距[2-5]和不同射流口形狀[4,8]對流場混合點位置、合并點位置和射流寬度的影響以探究雙射流間相互作用。劉鵬遠等[4]利用粒子圖像測速技術(particle image velocity, PIV) 實驗獲得雙射流速度場和湍流強度的分布,研究了高雷諾數(shù)下不同射流口間距對平行射流合并點位置的影響。Tanaka[7]對不同射流口間距的雙射流流場開展了實驗研究,分析了射流口間距對合并點后射流寬度的影響。與此同時,前人也利用數(shù)值模擬方法來研究流場特性,Mondal 等[9]對二維雙射流進行了數(shù)值模擬,指出混合點位置與射流口間距存在線性關系。Li 等[10]在對雙圓柱射流的大渦模擬數(shù)值研究中指出射流口間距和混合點之間存在指數(shù)對應關系。Zhou 等[5]根據(jù)直接數(shù)值模擬結果提出了雙射流中射流流場相互作用尺度(jet-interaction length scale),指出該尺度可以預測雙射流流場特性并根據(jù)該尺度解釋了射流口間距與混合點位置的線性關系。

    除流動空間演化特性相關研究外,工業(yè)中所遇到的射流問題往往涉及兩個或更多不同組分標量的擴散、混合和化學反應,對此部分研究人員進行了相關實驗研究[11-14]和數(shù)值模擬[15-18],旨在預測相關過程中混合速率以及控制化學反應速率。Stapountzis 等[11]采用化學敏感激光誘導熒光技術(laser-induced fluorescence, LIF)研究了反應性雙平面射流,指出標量混合程度取決于射流口間距和射流軸向位置且流向最大平均濃度位置按照冪律增長。Soltys 等[12]使用雙通道平面激光誘導熒光(planar laser-induced fluorescence, PLIF)系統(tǒng)量化了雙圓孔射流從發(fā)射到合并過程中兩個獨立標量的瞬時協(xié)方差空間結構,研究了標量協(xié)方差空間結構對混合過程以及化學反應的影響。Watanabe等[13]基于光吸收光譜法濃度測量系統(tǒng)實驗研究了反應性平面液體射流,結果表明反應物的濃度相關性越大平均反應速率越小。Cai 等[14]采用平面激光誘導熒光法和瑞利散射法實驗研究了湍流同軸射流中的三標量混合,利用標量平均值、方差、相關系數(shù)、聯(lián)合概率密度分布(joint probability distribution function, JPDF)等詳細分析了近場流動中的混合過程,為理解和模擬反應流中的多尺度混合提供了基礎。就湍流化學反應的數(shù)值模擬而言,考慮復雜詳細機理會導致控制方程龐大,故用單步反應來簡化詳細化學反應機理以降低求解難度[19-20]。早在1962 年,Toor[21]進行了化學反應簡化機理的研究,將傳質(zhì)-反應方程組簡化為無反應傳質(zhì)方程并分析了具有等擴散系數(shù)的低濃度物質(zhì)間快速均勻不可逆反應的傳質(zhì)過程。Bilger[22]針對非預混射流火焰發(fā)展了流體混合和化學反應的理論,用一個守恒的標量來充分描述混合,導出了任意組分的瞬時反應速率的表達式。但目前關于反應性雙射流的數(shù)值模擬研究很少,僅有部分關于反應性單射流的相關模擬。DesJardin 等[15]對非預混反應單射流進行了大渦模擬,評估了幾種亞格子模型的優(yōu)缺點。Watanabe 等[16-18]對具有二級等溫化學反應(A+BP)的平面單射流進行了系列數(shù)值模擬,分析了近場湍流/非湍流界面層附近的標量特性。

    據(jù)文獻統(tǒng)計,前人關于反應性雙射流的研究大都基于實驗分析,缺乏對反應性雙射流中射流相互作用和組分相互作用影響機理的認識,亟待開展雙平行射流標量輸運和化學反應特性的數(shù)值模擬研究。本文對分別包含兩種反應物的雙射流之間的二級非平衡基元反應(A+BR)進行數(shù)值模擬,探究射流相互作用尺度在雙平行反應性射流標量場的適用性,研究化學反應方程非穩(wěn)態(tài)項、對流項、擴散項和化學反應源項對標量生成、消耗和輸運的影響,旨在為射流混合及化學反應等領域的優(yōu)化設計提供指導。

    1 數(shù)值方法

    1.1 流場計算

    研究基于OpenFOAM 中pimpleFoam 求解器中結合非迭代瞬態(tài)PISO算法的PIMPLE算法。采用大渦模擬方法,濾波后的不可壓縮流動控制方程組表達如下:

    1.2 標量場計算

    數(shù)值模擬使用類似于對動量方程的處理方式對化學反應控制標量輸運方程進行求解。濾波后的化學反應控制標量輸運方程表達如下:

    式中,k為該二級反應的反應速率常數(shù)。

    1.3 計算模型和參數(shù)

    研究具體考慮分別含有不同反應物的兩股射流之間的化學反應,不考慮射流與環(huán)境流的反應。圖1 為化學反應流場和坐標系的示意圖。反應物A和B 分別預先混合到射流Ⅰ和Ⅱ中。射流Ⅰ、Ⅱ通過寬度為d、射流口間距為Ld的兩個狹縫噴射到環(huán)境流中。反應物A和B隨著雙射流在收斂區(qū)以對流和擴散的方式向周圍環(huán)境及下游輸運,在中心線附近反應物相互接觸產(chǎn)生二級非平衡基元反應A+BR。

    圖1 數(shù)值模擬示意圖Fig.1 Schematic diagram of the simulation

    本文對不同射流口間距的三個雙射流工況進行模擬,以探究射流口間距對流動-化學反應耦合過程的影響。表1為具體的計算參數(shù)及流場網(wǎng)格節(jié)點設置,Run J4、Run J6和Run J8分別對應射流口間距Ld/d=4, 6 和8。其中,雷諾數(shù)基于射流口平均流向速度UJ和射流口寬度d,Red=UJd/ν。UA為環(huán)境流速度,ΓA0和ΓB0分別為反應物A 和反應物B 的初始濃度。Damk?hler 數(shù)(Da)是表示流動時間尺度和化學反應時間尺度之比的無量綱數(shù),本文中其值由Da=k(ΓA0+ΓB0)d/UJ定義。設置計算條件如下:對于流場,射流口設置兩個相等的平均流向速度UJ=UJetⅠ=UJetⅡ。為使雙射流模擬盡快達到湍流狀態(tài),采用白噪聲脈動法(“white noise”fluctuations)[23]在兩個射流口處x,y和z方向添加速度脈動。本文中兩個射流口處三個方向上添加的速度脈動均服從(-0.05, 0.05)均勻分布,其均方根值為Urms=Vrms=Wrms=0.029,該值小于Watanabe 等[18]研究中的射流口速度脈動均方根值最大值0.1,滿足白噪聲脈動法添加速度脈動的原則,即該值既能觸發(fā)層流向湍流轉換,又能保證其對流場統(tǒng)計結果影響不大。為降低環(huán)境流速度對射流運動的影響,研究采用與Zhou等[5]和Stanley 等[24]直接數(shù)值模擬中相同的環(huán)境流速度(UA/UJ=0.1)。流向出口采用對流流出邊界條件,法向邊界采用自由滑移邊界條件,展向邊界采用周期性邊界條件。對于組分濃度場,設置初始濃度比ΓA0/ΓB0=2,產(chǎn)物初始濃度為0。設置施密特數(shù)Sc=0.71。采用OpenFOAM 中時間離散和空間離散格式。其中,時間遞進采用二階隱式backward,流場對流格式采用Gauss linear,標量場對流格式采用Gauss limitedLinear,擴散格式采用Gauss linear corrected,化學反應源項采用自動顯式隱式fvm::SuSp。為節(jié)省計算資源同時精確捕捉雙射流法向流場及標量場特性,在雙射流法向發(fā)展范圍即y方向(-10,10)范圍內(nèi)采用加密網(wǎng)格。

    表1 反應性雙射流計算參數(shù)及節(jié)點設置Table 1 Simulation parameters and node settings of reactive dual jets

    2 數(shù)值方法驗證

    截至目前,國內(nèi)外學術界尚未有關于兩股射流之間化學反應的數(shù)值模擬或?qū)嶒炑芯?,無法直接驗證數(shù)值方法的準確性。單射流和雙射流同屬于自由剪切射流且均存在湍流/非湍流界面,有類似的標量輸運和化學反應特性,本文認為進行反應性單射流算例模擬并將模擬結果與前人研究結果進行對比同樣能夠證明當前數(shù)值方法的準確性。單射流算例中攜帶反應物A 的射流通過寬度為d的射流口射入到包含反應物B 的環(huán)境流中,隨著流場發(fā)展單射流與環(huán)境流發(fā)生化學反應,計算條件設置與雙射流算例相同。

    圖2為單射流算例無量綱平均流向速度的法向分布及其與前人的實驗結果[25-26]和數(shù)值模擬結果[17]對比。其中,符號“〈〉”表示對時間以及z方向平均,UC為射流中心線上平均速度,bU為基于流向速度U的射流半寬,圖中的Zone-averaged 表示達到自相似后對空間平均(z方向平均、x方向平均或y方向折疊平均)。當前流場在下游x/d=20 和x/d=30 處達到自相似,且與前人實驗結果[25-26]和數(shù)值模擬結果[17]一致。圖3 為無量綱速度脈動均方根值的法向分布,由于速度脈動均方根法向分布受初始場影響很大,前人實驗結果[25-27]和數(shù)值模擬結果[24,28-30]中速度脈動均方根法向分布在一定范圍變化,而當前模擬結果在該范圍內(nèi)表現(xiàn)出良好一致性。流向速度和速度脈動均方根的法向分布均與前人結果達到良好吻合,驗證了流場模擬的準確性。

    圖2 單束反應射流無量綱平均流向速度法向分布Fig.2 Vertical distribution of normalized mean streamwise velocity of a reactive jet

    圖3 單束反應射流無量綱速度脈動均方根法向分布Fig.3 Vertical distribution of normalized turbulence rms value of a reactive jet

    為了進一步驗證化學反應模擬的準確性,驗證混合濃度分數(shù)ξ的標量統(tǒng)計,ξ由式(6)給出:

    在相關數(shù)值計算中ξ可以提高化學反應與組分混合的計算效率。反應物A 和B的瞬時濃度可以由混合濃度分數(shù)和產(chǎn)物濃度表達如下:

    圖4(a)為無量綱平均混合濃度分數(shù)ξ的法向分布。ξ以射流中心線上的平均混合濃度分數(shù)〈ξ〉C無量綱,法向坐標y以基于ξ的射流半寬bξ無量綱?;旌蠞舛确謹?shù)在下游達到自相似且與Watanabe 等[18]的結果吻合。圖4(b)為混合濃度分數(shù)脈動均方根的法向分布?;旌蠞舛确謹?shù)脈動均方根與Watanabe等[18]數(shù)值模擬結果略有偏差,這一偏差是由驗證算例所加射流口初始速度脈動強度與Watanabe 等[18]數(shù)值模擬所加速度脈動強度不同造成的。由此可見,當前標量場也與前人結果吻合較好,驗證了化學反應模擬的準確性。

    圖4 單束反應射流無量綱平均混合濃度分數(shù)和混合濃度分數(shù)脈動均方根法向分布Fig.4 Vertical distribution of mean mixture fraction and mixture fraction fluctuation rms of a reactive jet

    3 結果分析與討論

    3.1 流場和標量場可視化

    研究人員根據(jù)混合點和合并點的位置將雙射流流場結構分為收斂區(qū)、混合區(qū)和合并區(qū)三個區(qū)域[6-7],圖5為典型的雙射流流場結構圖。其中,流場中回流出現(xiàn)且流向速度為零的點為準滯止點xqsp,回流消失且流向速度為零的點為混合點xmp,流場中心線流向平均速度最大的點為合并點xcp。收斂區(qū)指混合點之前的區(qū)域,該區(qū)域內(nèi)射流對環(huán)境流的卷吸作用會使兩股射流間形成低壓區(qū)?;旌蠀^(qū)指混合點與合并點之間的區(qū)域。兩股射流激烈碰撞劇烈混合會形成高壓區(qū)。由高壓區(qū)和低壓區(qū)綜合作用會導致在混合點前形成回流區(qū)。合并區(qū)指合并點之后的區(qū)域,合并區(qū)內(nèi)兩股射流完全合并為單射流并呈一定單射流特性[31]。

    圖5 雙射流流場結構圖Fig.5 Structure diagram of dual jet flow field

    雙射流中心線上兩股射流間相互作用以及化學反應最為劇烈,中心線上的統(tǒng)計特性能直觀地反映射流相互作用和組分相互作用對雙射流流動和化學反應耦合過程的影響,因此下文考察雙射流沿中心線上的流場統(tǒng)計量和標量場統(tǒng)計量的流向演化。圖6為三種射流口間距下中心線上平均流向速度流向演化。隨著流場的發(fā)展平均流向速度均有先減小后增大的趨勢,并在上游都有負值出現(xiàn),這說明在三種不同射流口間距的情況下,均會出現(xiàn)回流現(xiàn)象。表2 給出圖6 三種射流口間距下準滯止點位置xqsp、混合點位置xmp和合并點位置xcp,可以看出準滯止點位置xqsp、混合點位置xmp和合并點位置xcp都隨著射流口間距的增大而增大,且呈一定線性規(guī)律。

    圖6 中心線上平均流向速度流向演化Fig.6 Streamwise evolution of the mean streamwise velocity along the center line

    表2 三種射流口間距下準滯止點位置xqsp、混合點位置xmp和合并點位置xcpTable 2 The quasi stagnation point xqsp,locations of the merge point xmp,and the combine point xcp in the three cases of different jet separation distances

    圖7(a)為算例Run J4、Run J6 和Run J8 的中心線上流向速度脈動均方根的流向演化,流向速度脈動均方根在混合區(qū)內(nèi)均有一個峰值,該點的射流相互作用極強。圖7(b)為三種射流口間距下平均壓力的流向演化。整個回流區(qū)均為負壓,前人的研究中指出該負壓是由中心線上的湍流/非湍流界面和被包圍的非湍流區(qū)域之間的卷吸作用引起的[6-7]。同時,在回流區(qū)后存在一個高壓區(qū),其是由兩股射流在混合點碰撞沖擊引起的。

    圖7 中心線上流場統(tǒng)計特性流向演化Fig.7 Streamwise evolution of the flow characteristics along the center line

    在分析三種射流口間距下中心線上流場統(tǒng)計特性后,本文將探究反應物和生成物的標量統(tǒng)計特性。圖8~圖10分別為射流口間距Ld/d=6下各組分同一時刻瞬時濃度云圖和平均濃度云圖。反應物A、B和產(chǎn)物R濃度分別以ГA0、ГB0和ГR0=ГA0ГB0/(ГA0+ГB0)無量綱化,這與Watanabe等的研究[16]中無量綱過程相同。圖中黑色分別代表反應物A 的濃度為ГA/ГA0=1、反應物B的濃度為ГB/ГB0=1以及產(chǎn)物R的濃度為ГR/ГR0=1,兩條黑色虛線代表回流區(qū)兩端位置。從圖8~圖10中可以看出,射流Ⅰ中的反應物A和射流Ⅱ中的反應物B隨著流場發(fā)展在中心線上接觸混合并擴散到彼此射流內(nèi)部,在此過程中發(fā)生化學反應并產(chǎn)生產(chǎn)物R,回流區(qū)內(nèi)產(chǎn)物濃度明顯高于其附近區(qū)域,且回流區(qū)后隨著流場發(fā)展產(chǎn)物的法向分布更為廣泛。

    圖8 反應物A的瞬時濃度和平均濃度Fig.8 Instantaneous concentration and mean concentration of reactant A

    圖10 生成物R的瞬時濃度和平均濃度Fig.10 Instantaneous concentration and mean concentration of product R

    根據(jù)Watanabe 等[17]的單束反應性射流標量場可視化研究,瞬時化學反應速率可以表達為W^R=DaГAГB/(ГA0ГB0),圖11 為關于產(chǎn)物R 的瞬時化學反應速率W^R云圖。圖中白色代表W^R= 0,黑色代表W^R= 0.01,三條黑色豎虛線分別代表準滯止點位置、混合點位置和合并點位置,黑色橫虛線代表流場中心線。由圖11可知,回流區(qū)和混合區(qū)化學反應主要發(fā)生在流場中心線附近,而在混合點后化學反應發(fā)生在圍繞流場中心線法向更廣的范圍,且在流場下游化學反應較弱。

    圖11 瞬時化學反應速率Fig.11 Instantaneous chemical reaction rate

    圖12(a)、(b)為三種射流口間距下中心線上反應物A 和B 無量綱平均濃度ГA/ГA0和ГB/ГB0流向演化。反應物濃度在回流區(qū)內(nèi)急劇升高,在混合點后維持在較高水平,這是因為回流作用將混合點附近雙射流攜帶的反應物輸運到回流區(qū)前部,混合點后隨著流場發(fā)展反應物不斷向中心線上輸運。圖12(c)為三種射流口間距下中心線上產(chǎn)物R的無量綱平均濃度ГR/ГR0流向演化,產(chǎn)物濃度沿中心線在回流區(qū)內(nèi)存在一個峰值。結合圖7(b)、圖8、圖9,回流區(qū)內(nèi)反應物A和B會充分混合并產(chǎn)生產(chǎn)物R,而混合點兩股射流碰撞產(chǎn)生的高壓區(qū)會阻礙產(chǎn)物向下游輸運,從而導致存在上述峰值。此外,隨著射流口間距增大,反應物平均濃度峰值減小且峰值位置后移,產(chǎn)物平均濃度峰值增大且峰值位置后移。

    圖12 中心線上平均反應物濃度和平均產(chǎn)物濃度流向演化Fig.12 Streamwise evolution of the mean concentration of reactant and product along the center line

    圖13(a)、(b)分別為三種射流口間距下反應物A和反應物B的無量綱濃度脈動均方根ГArms/ГA0和ГBrms/ГB0沿中心線的流向演化。反應物濃度脈動均方根在回流區(qū)急劇增加,這與前人關于雙射流湍動能研究中湍動能流向演化趨勢類似[32],該急劇變化是由回流區(qū)卷吸作用的擾動所引起的。圖13(c)為三種射流口間距下產(chǎn)物R無量綱濃度脈動均方根ГRrms/ГR0沿中心線的流向演化。產(chǎn)物濃度脈動均方根在流場上游存在兩個峰值且峰值位置分別與準滯止點和混合點接近。第一個峰值出現(xiàn)的原因是回流區(qū)前反應物A 和B 在準滯止點附近一個反應-不反應界面剛開始接觸并發(fā)生反應,此時產(chǎn)物R 濃度脈動極強。第二個峰值出現(xiàn)的原因是混合點附近反應物A和B隨著兩股射流碰撞劇烈反應。上述演化都是對流、化學反應、擴散過程綜合作用的結果,3.3節(jié)將結合標量輸運方程給出詳細分析。本節(jié)對流場和標量場的分析表明,隨著射流口間距的改變反應物和產(chǎn)物的標量場都呈一定規(guī)律的改變,流動在準滯止點和混合點附近的行為直接影響著化學反應。

    圖13 中心線上反應物濃度脈動均方根和產(chǎn)物濃度脈動均方根流向演化Fig.13 Streamwise evolution of concentration fluctuation rms value of reactant and product along the center line

    3.2 射流相互作用尺度

    3.1 節(jié)分析發(fā)現(xiàn)中心線上湍流流場和標量場統(tǒng)計量的流向演化與射流口間距存在一定規(guī)律。Zhou等[5]基于單射流冪率和湍流/非湍流界面層概率密度函數(shù)首次提出射流相互作用尺度x*(x*=Ld),指出不同射流口間距下湍流特性的流向演化均以x*為尺度,即當以x*對流向距離x無量綱時不同射流口間距下中心線上湍流統(tǒng)計量的流向演化曲線會重合。本節(jié)將射流相互作用尺度引入反應性雙射流,探究反應性雙射流中標量統(tǒng)計量的空間演化是否遵循上述規(guī)律。

    圖14(a)、(b)分別為三種射流口間距下中心線上無量綱反應物A 平均濃度〈ГA〉/〈ГA〉max和反應物B 平均濃度〈ГB〉/〈ГB〉max的流向演化,其中流向距離x以x*無量綱化。在流場中心線上整個流域無量綱反應物平均濃度曲線有相當好的重合,這表明x*可以作為潛在的描述反應物標量統(tǒng)計量的相互作用尺度。為了充分考察x*在反應物標量場的適用性,圖15(a)、(b)給出了三種射流口間距下中心線上無量綱反應物濃度脈動均方根ГArms/ГArms,max和ГBrms/ГBrms,max流向演化情況,其中流向距離x以x*無量綱化。從兩圖可以看出反應物濃度脈動均方根三條曲線在整個流域也均有相當好的重合,這表明x*可以作為描述反應物標量統(tǒng)計量的相互作用尺度。值得一提的是,Run J4 的曲線與Run J6 和Run J8 的曲線略有偏差,本文認為這是由射流口間距過小進而出現(xiàn)了Mondal 等[9]提出的周期性渦脫落現(xiàn)象所導致的,這與射流口間距在6<Ld<8 范圍內(nèi)反應物標量統(tǒng)計量的空間演化規(guī)律并不沖突。

    圖14 反應物A和反應物B平均濃度沿著中心線上x/x*的流向演化Fig.14 Streamwise evolution of the mean concentration of reactant A and reactant B along x/x*on the center line

    圖15 反應物A和反應物B濃度脈動均方根沿著中心線上x/x*的流向演化Fig.15 Streamwise evolution of reactant A and reactant B concentration fluctuation rms value along x/x*on the center line

    圖16(a)、(b)分別為三種射流口間距下中心線上無量綱產(chǎn)物R 平均濃度〈ГR〉/〈ГR〉max和濃度脈動均方根ГRrms/ГRrms,max的流向演化,其中流向距離x以x*無量綱化。圖16(a)中Run J6 和Run J8 的無量綱平均產(chǎn)物濃度曲線在整個流域內(nèi)有良好重合,表明x*可以作為潛在的描述產(chǎn)物標量統(tǒng)計量的相互作用尺度。圖16(b)中Run J6 和Run J8 的無量綱產(chǎn)物濃度脈動均方根曲線在流場上游(x/x*<1.2)有較好重合,這表明x*可以作為描述產(chǎn)物標量統(tǒng)計量在上游演化的相互作用尺度。此外,Run J4 的曲線與Run J6 和Run J8 的曲線在下游偏差較大,這表明周期性渦脫落現(xiàn)象對產(chǎn)物濃度脈動的影響比對反應物濃度脈動的影響更大。

    圖16 產(chǎn)物R平均濃度和濃度脈動均方根沿著中心線上x/x*的流向演化Fig.16 Streamwise evolution of concentration and concentration fluctuation rms value of product R along x/x*on the center line

    表3給出了三種射流口間距下無量綱射流相互作用尺度x*、反應物A 濃度脈動均方根峰值位置xpeak、產(chǎn)物R 濃度脈動均方根峰值位置(xpeak1和xpeak2)和對應的無量綱化比值。其中,比值分別由Run J4算例的對應值(x*1,xpeak1,xpeak1,1和xpeak2,1)無量綱化。可以看出反應物濃度脈動峰值位置xpeak和產(chǎn)物濃度脈動峰值位置(xpeak1和xpeak2)均與射流相互作用尺度x*呈一定的線性關系。此外,本文也以xpeak等尺度對流向距離無量綱化得到了不同射流口間距下的反應物和產(chǎn)物標量統(tǒng)計量的流向演化曲線,其結果與圖14~圖16 中曲線類似,在此不再贅述。上述結果表明xpeak,xpeak1和xpeak2也可以作為描述反應性雙射流中標量統(tǒng)計量的相互作用尺度。

    表3 無量綱射流相互作用尺度x*、反應物濃度脈動均方根峰值位置xpeak和產(chǎn)物濃度脈動均方根峰值位置xpeak 1、xpeak 2Table 3 Normalized wake-interaction scale x*,locations of reactant concentration fluctuation rms peak xpeakandproduct concentration fluctuation rms peak xpeak 1,xpeak 2

    由于本文考慮的是Da=0.1 的弱化學反應過程,反應物標量特性主要受雙射流間相互作用影響,當前雙射流相互作用尺度可以預測整個流域的反應物平均濃度和反應物濃度脈動均方根的流向演化。產(chǎn)物標量特性受射流間相互作用和標量相互作用影響,故產(chǎn)物平均濃度和產(chǎn)物濃度脈動均方根的流向演化較為復雜。盡管如此,當前相互作用尺度也可以預測上游(x/x*<1.2)的產(chǎn)物標量統(tǒng)計量的流向演化。綜上討論,在一定的射流口間距范圍內(nèi)(6<Ld<8),如果使用適當?shù)某叨龋ㄈ绫疚闹械膞*和xpeak等)對流向距離無量綱化,不同射流口間距下射流相互作用中的標量統(tǒng)計特性演化遵循相似的規(guī)律。因此,下文的研究只考慮Ld/d=6一種情況。

    3.3 標量產(chǎn)生、消耗和輸運

    3.3.1 化學反應標量輸運方程 化學反應標量輸運方程是研究標量產(chǎn)生、消耗和輸運的重要方法。反應物A、B 和產(chǎn)物R 的化學反應標量輸運方程分別為:

    圖17(a)、(b)分別為關于反應物A 和B 的化學反應標量輸運方程中各項瞬時值的流向演化。RHS為殘差項表示公式左右兩側代數(shù)差。RHS 在中心線流向上幾乎為0,表明標量場統(tǒng)計結果有良好收斂性。圖中三條虛線分別代表流場的準滯止點位置(Ld/d=6)、混合點位置(Ld/d=8.5)和合并點位置(Ld/d=20)。因為兩束射流所攜帶的反應物還沒有被輸運到流場中心線,化學反應尚未發(fā)生,所以在流場上游回流區(qū)前(x/d<6)標量輸運方程中各項均為0。當x/d>6 時,對流項A和非穩(wěn)態(tài)項At沿著中心線呈正負值交替變化,而擴散項D和化學反應源項S相對于對流項A和非穩(wěn)態(tài)項At可以忽略。這表明對流輸運作用主導反應物濃度變化,而化學反應和擴散作用較弱且對反應物濃度標量的輸運幾乎沒有影響,故式(9)和式(10)可以分別簡化為:

    此外,在圖17(a)、(b)中還可以看出回流區(qū)內(nèi)(6<x/d<8.5)對流項的值比混合區(qū)(8.5<x/d<20)和合并區(qū)(x/d>20)對流項的值小得多。這是因為兩股射流在混合點直接碰撞,碰撞前回流引起的輸運作用較小而碰撞后對流作用較強。

    圖18 為產(chǎn)物R 的化學反應標量輸運控制方程中各項瞬時值的流向演化。當x/d>6 時,對流項AR和非穩(wěn)態(tài)項AtR沿著中心線呈正負值交替變化,化學反應源項SR相對于對流項AR和非穩(wěn)態(tài)項AtR較小,擴散項DR可以忽略。這表明化學反應產(chǎn)生的產(chǎn)物被對流作用輸運,對流作用仍主導反應物濃度變化,擴散作用對產(chǎn)物的輸運幾乎沒有影響,故式(11)可以簡化為:

    由式(12)~式(14)、圖17、圖18 可以看出對流作用和非穩(wěn)態(tài)變化之間存在一定負相關性,引入相關系數(shù)ρ[33]以定量探究對流項和非穩(wěn)態(tài)項的相關性以及控制方程中其他各項彼此間相關性。以對流項A和非穩(wěn)態(tài)項At為例,相關系數(shù)由式(15)給出:

    圖17 反應物標量輸運方程各項沿著中心線的瞬時流向演化Fig.17 Streamwise evolution of the instantaneous terms in reactant scalar transport equation along the center line

    式中,σ為標準差。由于x/d<6時各項均為0,本文研究x/d>6 時反應物和產(chǎn)物標量輸運方程中各項之間相關性。圖19(a)、(b)分別為反應物A 和B 標量輸運方程中各項間相關系數(shù)在中心線上流向演化,ρ(AA,AtA)和ρ(AB,AtB)幾乎保持在-1,ρ(DA,AtA)、ρ(DB,AtB)、ρ(SA,AtA) 和ρ(SB,AtB)幾乎為0,這表明反應物濃度非穩(wěn)態(tài)變化幾乎由對流過程主導,與上文圖18的分析一致。反應物標量輸運方程中其他各項間相關 系 數(shù)[ρ(AA,DA)、ρ(AA,SA)、ρ(DA,SA)、ρ(AB,DB)、ρ(AB,SB)和ρ(DB,SB)]很小,說明擴散過程、化學反應過程和對流過程彼此間影響很小。圖19(c)為產(chǎn)物R標量輸運方程中各項相關系數(shù)在中心線上流向演化。ρ(AR,AtR)幾乎保持在-1,表明產(chǎn)物R濃度非穩(wěn)態(tài)變化也幾乎由對流過程主導。ρ(SR,AtR)相關系數(shù)在混合區(qū)(8.5<x/d<20)的值大于其他區(qū)域,表明混合區(qū)化學反應較為劇烈,對產(chǎn)物R 濃度非穩(wěn)態(tài)變化的影響減弱。此外,ρ(DR,SR)維持在-0.4~-0.15 之間,表明化學反應產(chǎn)生的產(chǎn)物R也由擴散過程輸運。

    圖18 產(chǎn)物標量輸運方程各項沿著中心線的瞬時流向演化Fig.18 Streamwise evolution of the instantaneous terms in product scalar transport equation along the center line

    上文提到圖19(c)中的ρ(AR,AtR)幾乎保持在-1,對流項與非穩(wěn)態(tài)項呈強烈負相關特性,為進一步探究這一負相關性,圖20 給出了非穩(wěn)態(tài)項脈動AtR-〈AtR〉與對流項脈動AR-〈AR〉在x/d=10,20和30處的聯(lián)合概率密度分布(JPDF)。標示的六個等值線的值分別為0.3,0.1,0.03,0.01,0.003 和0.001。AtR-〈AtR〉與AR-〈AR〉在三個位置的JPDF 幾乎均擠壓在y=-x直線上,印證了產(chǎn)物R 濃度非穩(wěn)態(tài)變化也幾乎由對流過程主導。該JPDF在x/d=10處被擠壓效果不如在x/d=20和30處明顯,這是因為x/d=10所在的混合區(qū)化學反應相對較強,此時非穩(wěn)態(tài)項和對流項間相關性較弱。

    圖19 化學反應輸運方程各項間相關系數(shù)沿著中心線的流向演化Fig.19 Streamwise evolution of the correlation coefficient between items of the chemical reaction transport equation along the center line

    圖20 AtR-〈AtR〉和AR-〈AR〉在x/d=10,20和30處的聯(lián)合概率密度分布Fig.20 JPDF of AtR-〈AtR〉and AR-〈AR〉at x/d=10,20,and 30

    對圖19(c)的分析中還指出ρ(DR,SR)維持在-0.4~-0.15 之間,擴散項與化學反應源項有弱負相關性。同樣,為進一步探究這一負相關性,圖21 給出擴散項脈動DR-〈DR〉與化學反應源項脈動SR-〈SR〉在x/d= 10, 20 和30 處的聯(lián)合概率密度分布。標示的五個等值線的值分別為0.1,0.03,0.01,0.003 和0.001。可以看出JPDF 在流向三個位置關于y= -x都有一定的對稱性,這與相關系數(shù)ρ(DR,SR)的流向演化相符。

    圖21 DR-〈DR〉和SR-〈SR〉在x/d=10,20和30處的聯(lián)合概率密度分布Fig.21 JPDF of DR-〈DR〉and SR-〈SR〉at x/d=10,20,and 30

    圖22(a)~(c)為各項平均值的流向演化。對流項〈A〉和化學反應源項〈S〉的值在整個流域幾乎相等,這表明對流運動與化學反應引起的產(chǎn)物R的濃度變化收支平衡。關于產(chǎn)物R的穩(wěn)態(tài)方程可以寫為:

    圖22 化學反應標量輸運方程中各項平均值沿著中心線的流向演化Fig.22 Streamwise evolution of the mean value of each term in chemical reaction scalar transport equation along the center line

    此外,沿著流場中心線有兩個化學反應較強的峰值點?;亓髯饔檬狗磻顰 和B 不能向下游輸運,在混合點前充分混合并反應,回流區(qū)域內(nèi)出現(xiàn)第一個峰值。隨著流場發(fā)展兩股射流逐漸吸引,在混合點合并為一股射流后對流作用不再增強,合并點后出現(xiàn)第二個峰值。

    在當前研究中(Da=0.1 且Sc=0.71)雙射流標量特性受化學反應影響相對較弱。從圖17、圖18瞬態(tài)輸運方程各項的流向演化可以看出,由雙射流攜帶的反應物A 和B 幾乎由對流過程輸運,化學反應引起反應物消耗很小,由化學反應產(chǎn)生的產(chǎn)物R 也幾乎由對流過程輸運。從圖22 輸運方程各項平均值的流向演化可以看出,對流項與化學反應源項相互抵消,混合點前回流作用使化學反應增強。流場上游化學反應最為劇烈的點存在于回流區(qū)內(nèi)回流速度最大的點附近,流場下游化學反應較為劇烈的點出現(xiàn)在合并點后一段距離。

    3.3.2 反應物和產(chǎn)物濃度相關性 反應物和產(chǎn)物各濃度間相關性可以反映出流動與化學反應的耦合關系。圖23 為各組分濃度相關系數(shù)在中心線上的流向演化,反應物A 瞬時濃度和反應物B 瞬時濃度相關系數(shù)ρ(ΓA,ΓB)維持在-0.9,呈現(xiàn)強烈負相關特性。圖24 為反應物輸運方程中對流項AA和AB間相關系數(shù)ρ(AA,AB)的流向演化。自混合區(qū)后對流項AA和對流項AB呈負相關,這導致了反應物A 和反應物B 濃度呈負相關,這一結果與圖17 的結論相符。反應物濃度的負相關特性表明隨著化學反應的發(fā)生混合區(qū)后射流中心線上反應物A的出現(xiàn)會引起反應物B 的減少,同時反應物B 的出現(xiàn)也會引起反應物A 的減少。此外,反應物A 和產(chǎn)物R 濃度之間相關系數(shù)ρ(ΓA,ΓR)維持在0 左右,這是因為組分濃度間的相關性與化學反應有關,而射流中心線上反應物A濃度變化主要由對流作用引起且化學反應對反應物濃度影響較弱,故反應物A 和產(chǎn)物R 濃度之間相關性很弱。反應物B和產(chǎn)物R 濃度之間相關系數(shù)ρ(ΓB,ΓR)在-0.25~0 之間變化,同樣呈一定弱負相關特性。ρ(ΓA,ΓR)和ρ(ΓB,ΓR)的值不同是由反應物B的初始濃度比反應物A的初始濃度小導致的。

    圖23 組分濃度間相關系數(shù)沿著中心線的流向演化Fig.23 Streamwise evolution of the correlation coefficient between the component concentrations along the center line

    圖24 反應物輸運方程中對流項AA和AB間相關系數(shù)ρ(AA,AB)沿著中心線的流向演化Fig.24 Streamwise evolution of the correlation coefficient ρ(AA,AB)between the advection term AA and AB in the transport equation of reactant along the center line

    反應物濃度脈動和產(chǎn)物濃度脈動的聯(lián)合概率密度分布一定程度上可以反映化學反應對各組分濃度的影響。圖25 給出了反應物A 和產(chǎn)物R 濃度脈動之間的JPDF,標示的五個等值線的值分別為0.1, 0.03, 0.01 , 0.003 和0.001。在x/d=10,20 和30處,該JPDF 呈“心形”。此外,反應物B 和產(chǎn)物R 濃度脈動之間JPDF的形狀與其幾乎相同,本文不再給出。從JPDF 可以看出在三個位置處反應物A 濃度和反應物B 濃度脈動較大時,產(chǎn)物R 濃度大概率維持在較低水平。只有反應物A 瞬時濃度和反應物B瞬時濃度都在平均值時,產(chǎn)物R 濃度有較大概率維持在較高水平,有利于化學反應的發(fā)生。

    圖25 ΓA-〈ΓA〉和ΓR-〈ΓR〉在x/d=10,20 和30 處的聯(lián)合概率密度分布Fig.25 JPDF of ΓA-〈ΓA〉and ΓR-〈ΓR〉at x/d=10,20,and 30

    4 結 論

    對具有二級非平衡基元反應的雙射流進行了數(shù)值模擬研究,具體研究了三種射流口間距下在Da= 0.1、Sc= 0.71 時化學反應、對流作用和擴散作用對反應物消耗、產(chǎn)物產(chǎn)生和組分輸運的影響。主要結論如下。

    (1)流動在準滯止點和混合點的行為直接影響著化學反應。流場上游產(chǎn)物濃度均方根脈動有兩個峰值,兩個峰值位置分別在準滯止點xqsp和混合點xmp附近。

    (2)x*可以作為預測反應性雙射流中標量輸運特性的射流相互作用尺度。在一定的射流口間距范圍內(nèi),使用適當?shù)某叨龋ㄈ鐇*和xpeak)對流向距離無量綱化,不同射流口間距下標量統(tǒng)計量演化遵循相似規(guī)律。

    (3)在Da=0.1 且Sc=0.71 時,雙射流標量輸運特性受化學反應影響較弱。由雙射流攜帶的反應物A和B幾乎由對流過程輸運(AtA=-AA,AtB=-AB),由化學反應產(chǎn)生的產(chǎn)物R 也幾乎由對流過程輸運(AtR=-AR+SR)。統(tǒng)計平均下對流運動與化學反應引起的產(chǎn)物R的濃度變化收支平衡(〈AR〉=〈SR〉)。

    (4) 反應物A 濃度脈動和反應物B 濃度脈動具有較強負相關性。反應物A瞬時濃度和反應物B瞬時濃度均在平均值時,產(chǎn)物R 濃度有較大概率維持在較高水平,有利于化學反應的發(fā)生。

    符 號 說 明

    A——對流項

    At——非穩(wěn)態(tài)項

    D——擴散項

    Dα——分子擴散系數(shù),m2/s

    P——壓力,Pa

    S——化學反應源項

    t——時間,s

    Ui——i方向上的速度U、V和W,m/s

    xi——x、y和z方向坐標,m

    Γ——物質(zhì)的量濃度,mol/L

    ν——運動黏度,m2/s

    下角標

    C——雙射流中心線

    i——自由指標

    j——啞指標

    α——反應物或產(chǎn)物組分A、B或R

    猜你喜歡
    標量反應物流向
    小溪啊!流向遠方
    井岡教育(2020年6期)2020-12-14 03:04:42
    一種高效的橢圓曲線密碼標量乘算法及其實現(xiàn)
    一種靈活的橢圓曲線密碼并行化方法
    初中化學中氣體的制取、凈化與干燥
    中學化學(2017年5期)2017-07-07 08:40:47
    十大漲幅、換手、振副、資金流向
    化學反應中的能量變化考點點擊
    中學化學(2016年4期)2016-05-30 16:20:37
    流向逆轉的啟示
    化學平衡移動對反應物轉化率的影響
    中學化學(2014年1期)2014-04-23 08:59:04
    單調(diào)Minkowski泛函與Henig真有效性的標量化
    秋天的流向(組詩)
    鴨綠江(2013年11期)2013-03-11 19:42:04
    国产伦人伦偷精品视频| 99久久99久久久精品蜜桃| 一级黄色大片毛片| 午夜a级毛片| 国产精品爽爽va在线观看网站| 国产精华一区二区三区| 色av中文字幕| 欧美乱色亚洲激情| 亚洲国产精品999在线| 88av欧美| 欧美人与性动交α欧美精品济南到| 天堂影院成人在线观看| 夜夜看夜夜爽夜夜摸| 国产主播在线观看一区二区| 老鸭窝网址在线观看| 成年女人毛片免费观看观看9| 黑人操中国人逼视频| 久久久久久九九精品二区国产 | 男男h啪啪无遮挡| 不卡一级毛片| 老司机午夜十八禁免费视频| 99热这里只有是精品50| 女同久久另类99精品国产91| 岛国视频午夜一区免费看| 一边摸一边抽搐一进一小说| 午夜久久久久精精品| 777久久人妻少妇嫩草av网站| 美女午夜性视频免费| 欧美日韩乱码在线| 国产v大片淫在线免费观看| 在线观看舔阴道视频| www.自偷自拍.com| 一个人免费在线观看的高清视频| 亚洲免费av在线视频| 亚洲av片天天在线观看| 亚洲成人精品中文字幕电影| 国产成人精品无人区| 国产视频内射| 欧美精品啪啪一区二区三区| 国产主播在线观看一区二区| 午夜福利在线观看吧| 99热这里只有精品一区 | 亚洲中文字幕一区二区三区有码在线看 | 97人妻精品一区二区三区麻豆| 亚洲av电影在线进入| 久9热在线精品视频| 国产aⅴ精品一区二区三区波| 国产精品1区2区在线观看.| 色噜噜av男人的天堂激情| 亚洲七黄色美女视频| 国产单亲对白刺激| 国产精品99久久99久久久不卡| 日韩精品免费视频一区二区三区| 99riav亚洲国产免费| 日韩欧美国产在线观看| 全区人妻精品视频| 99精品欧美一区二区三区四区| www.www免费av| 久久久久精品国产欧美久久久| 黄色成人免费大全| 91麻豆精品激情在线观看国产| 亚洲五月天丁香| 国内久久婷婷六月综合欲色啪| 亚洲欧美精品综合一区二区三区| 男女床上黄色一级片免费看| 99久久国产精品久久久| 欧美日本视频| 国产又色又爽无遮挡免费看| 高清毛片免费观看视频网站| 啦啦啦观看免费观看视频高清| 两个人的视频大全免费| 19禁男女啪啪无遮挡网站| 别揉我奶头~嗯~啊~动态视频| 一级毛片女人18水好多| av天堂在线播放| 亚洲九九香蕉| 法律面前人人平等表现在哪些方面| 国产精品国产高清国产av| 五月伊人婷婷丁香| 午夜免费激情av| 久久久久久久久免费视频了| 国产精品影院久久| 欧美成人性av电影在线观看| 精品第一国产精品| 此物有八面人人有两片| 我要搜黄色片| 欧美一级毛片孕妇| 曰老女人黄片| 99精品欧美一区二区三区四区| 午夜精品久久久久久毛片777| 国产精品电影一区二区三区| 午夜福利18| 一个人免费在线观看电影 | 色老头精品视频在线观看| 国产精品av视频在线免费观看| 亚洲五月天丁香| 最新美女视频免费是黄的| 日本黄大片高清| 国产黄a三级三级三级人| 美女 人体艺术 gogo| 亚洲成av人片在线播放无| av天堂在线播放| 亚洲七黄色美女视频| 日本在线视频免费播放| 久久久久久免费高清国产稀缺| 亚洲精品美女久久av网站| 免费在线观看黄色视频的| netflix在线观看网站| 精品久久久久久,| 国产成人精品久久二区二区91| 法律面前人人平等表现在哪些方面| 黑人巨大精品欧美一区二区mp4| 久久久久久久久久黄片| 亚洲免费av在线视频| 久久久久国产精品人妻aⅴ院| avwww免费| 丝袜人妻中文字幕| 丁香六月欧美| 一区福利在线观看| 亚洲色图 男人天堂 中文字幕| 一区二区三区国产精品乱码| 国内揄拍国产精品人妻在线| 最近在线观看免费完整版| 看片在线看免费视频| 亚洲真实伦在线观看| 狠狠狠狠99中文字幕| √禁漫天堂资源中文www| 熟女电影av网| 成在线人永久免费视频| 丰满人妻熟妇乱又伦精品不卡| 久久精品91蜜桃| 亚洲精品粉嫩美女一区| 欧美丝袜亚洲另类 | 色综合亚洲欧美另类图片| 91国产中文字幕| 亚洲精品av麻豆狂野| 久久人人精品亚洲av| 岛国在线免费视频观看| 一级毛片女人18水好多| 搡老岳熟女国产| 亚洲精华国产精华精| 淫妇啪啪啪对白视频| 男人舔女人的私密视频| 最近最新中文字幕大全免费视频| 国内精品久久久久精免费| 国产精品久久久久久亚洲av鲁大| 国产精品 国内视频| 两个人看的免费小视频| 色在线成人网| 在线观看午夜福利视频| 国产亚洲精品一区二区www| av欧美777| 好男人在线观看高清免费视频| 亚洲国产日韩欧美精品在线观看 | 国产野战对白在线观看| 亚洲黑人精品在线| 午夜久久久久精精品| 18禁观看日本| 亚洲精品久久国产高清桃花| 好看av亚洲va欧美ⅴa在| 国产av一区二区精品久久| 制服诱惑二区| 国产精品亚洲美女久久久| 日韩欧美一区二区三区在线观看| 国产成人一区二区三区免费视频网站| 亚洲欧美精品综合一区二区三区| 国产精品自产拍在线观看55亚洲| 正在播放国产对白刺激| 日本免费一区二区三区高清不卡| 老鸭窝网址在线观看| 99精品在免费线老司机午夜| 亚洲国产看品久久| bbb黄色大片| 亚洲国产精品合色在线| 国产精品日韩av在线免费观看| 久久草成人影院| 国产亚洲精品久久久久久毛片| 国产一区二区在线av高清观看| 欧美丝袜亚洲另类 | 88av欧美| 欧美黄色淫秽网站| 悠悠久久av| 久久精品国产亚洲av高清一级| 亚洲av成人精品一区久久| 亚洲精品在线观看二区| 久久九九热精品免费| 欧美高清成人免费视频www| 国产不卡一卡二| 真人一进一出gif抽搐免费| 国产伦一二天堂av在线观看| 欧美午夜高清在线| 国产成人欧美在线观看| 午夜福利高清视频| 久久久久国内视频| 亚洲成av人片在线播放无| 日韩欧美一区二区三区在线观看| 亚洲欧美激情综合另类| 麻豆久久精品国产亚洲av| 亚洲av第一区精品v没综合| 九色国产91popny在线| 日本五十路高清| 久久精品成人免费网站| 麻豆一二三区av精品| 在线观看一区二区三区| 欧美色欧美亚洲另类二区| 麻豆成人av在线观看| 亚洲狠狠婷婷综合久久图片| 正在播放国产对白刺激| 国产一区在线观看成人免费| 麻豆国产av国片精品| 久久精品亚洲精品国产色婷小说| 亚洲熟妇熟女久久| 国产av一区在线观看免费| 嫩草影视91久久| 日韩av在线大香蕉| 中文资源天堂在线| 久久久水蜜桃国产精品网| 婷婷丁香在线五月| www.自偷自拍.com| 欧美色视频一区免费| 天天一区二区日本电影三级| a级毛片a级免费在线| www.熟女人妻精品国产| 国产一区二区在线观看日韩 | 久久精品综合一区二区三区| 日本一本二区三区精品| 婷婷六月久久综合丁香| 久久久精品国产亚洲av高清涩受| 三级国产精品欧美在线观看 | 蜜桃久久精品国产亚洲av| 欧美一区二区精品小视频在线| 啪啪无遮挡十八禁网站| 欧美绝顶高潮抽搐喷水| 欧美一区二区国产精品久久精品 | 美女大奶头视频| 最近在线观看免费完整版| 婷婷亚洲欧美| 欧美日韩亚洲综合一区二区三区_| 最近最新中文字幕大全电影3| 免费看a级黄色片| 午夜福利视频1000在线观看| 精品福利观看| 欧美成人性av电影在线观看| 国产精品亚洲美女久久久| 一二三四社区在线视频社区8| 国产v大片淫在线免费观看| 久久久久久九九精品二区国产 | 大型黄色视频在线免费观看| 精品欧美国产一区二区三| 成人高潮视频无遮挡免费网站| 女生性感内裤真人,穿戴方法视频| 欧美国产日韩亚洲一区| 国产又色又爽无遮挡免费看| 欧美中文日本在线观看视频| 婷婷精品国产亚洲av在线| 嫁个100分男人电影在线观看| 不卡av一区二区三区| 久久香蕉国产精品| 免费在线观看日本一区| 国产精品免费一区二区三区在线| 婷婷丁香在线五月| 成人三级黄色视频| 国产97色在线日韩免费| 亚洲全国av大片| 一边摸一边抽搐一进一小说| 国产乱人伦免费视频| 舔av片在线| 国产成人av激情在线播放| 一区福利在线观看| 人妻久久中文字幕网| 男人舔奶头视频| 亚洲成av人片免费观看| 国产视频一区二区在线看| 欧美日韩乱码在线| 国产亚洲精品一区二区www| 九色成人免费人妻av| 在线播放国产精品三级| www日本黄色视频网| 精品一区二区三区四区五区乱码| 国产精品综合久久久久久久免费| 后天国语完整版免费观看| av福利片在线观看| 一个人免费在线观看电影 | 精品人妻1区二区| 国产欧美日韩一区二区精品| 妹子高潮喷水视频| 久久性视频一级片| 欧美成人午夜精品| 久久天躁狠狠躁夜夜2o2o| 香蕉久久夜色| 最近视频中文字幕2019在线8| 精品久久久久久久久久免费视频| 久久精品国产清高在天天线| 国内久久婷婷六月综合欲色啪| 97人妻精品一区二区三区麻豆| 热99re8久久精品国产| 一区二区三区激情视频| 久久精品91无色码中文字幕| 黄色成人免费大全| 18美女黄网站色大片免费观看| 无限看片的www在线观看| 国产一区二区激情短视频| 亚洲av成人精品一区久久| 99久久精品热视频| 在线播放国产精品三级| av在线播放免费不卡| 成人国语在线视频| 日韩成人在线观看一区二区三区| 成年免费大片在线观看| 人人妻,人人澡人人爽秒播| 国产真实乱freesex| 女同久久另类99精品国产91| www.www免费av| 欧美日韩精品网址| АⅤ资源中文在线天堂| 最新在线观看一区二区三区| 久久久久久久久中文| 午夜老司机福利片| 久久 成人 亚洲| 日韩 欧美 亚洲 中文字幕| 1024视频免费在线观看| 国产激情欧美一区二区| 亚洲自拍偷在线| 高清毛片免费观看视频网站| 中文字幕高清在线视频| 亚洲专区中文字幕在线| 欧美日韩中文字幕国产精品一区二区三区| 亚洲人成伊人成综合网2020| 色老头精品视频在线观看| 亚洲成人国产一区在线观看| 中亚洲国语对白在线视频| 免费在线观看日本一区| 国产又色又爽无遮挡免费看| 成人av在线播放网站| 少妇熟女aⅴ在线视频| 少妇被粗大的猛进出69影院| 十八禁人妻一区二区| 日日爽夜夜爽网站| 国产男靠女视频免费网站| 最近最新中文字幕大全免费视频| 国产激情久久老熟女| 日韩中文字幕欧美一区二区| 国产又色又爽无遮挡免费看| 日本五十路高清| av有码第一页| 一边摸一边抽搐一进一小说| 操出白浆在线播放| 国产片内射在线| 国产野战对白在线观看| cao死你这个sao货| 一级毛片高清免费大全| 少妇被粗大的猛进出69影院| 禁无遮挡网站| netflix在线观看网站| 草草在线视频免费看| 真人一进一出gif抽搐免费| 亚洲 欧美一区二区三区| 嫩草影视91久久| 熟妇人妻久久中文字幕3abv| 久久久久久九九精品二区国产 | 精品国产乱子伦一区二区三区| 亚洲国产看品久久| 亚洲人成伊人成综合网2020| 不卡一级毛片| 老鸭窝网址在线观看| 国产亚洲精品av在线| 国产成年人精品一区二区| 欧美成人免费av一区二区三区| 91字幕亚洲| 男人的好看免费观看在线视频 | 岛国在线免费视频观看| 午夜精品一区二区三区免费看| 色综合亚洲欧美另类图片| 狂野欧美激情性xxxx| 成人特级黄色片久久久久久久| 两个人视频免费观看高清| 一级毛片高清免费大全| 国产91精品成人一区二区三区| 97超级碰碰碰精品色视频在线观看| 久久国产精品影院| 亚洲无线在线观看| www.精华液| 午夜福利18| 久久久久国内视频| 国产伦在线观看视频一区| 变态另类丝袜制服| 国产高清videossex| 搡老岳熟女国产| 在线看三级毛片| 妹子高潮喷水视频| 大型av网站在线播放| 欧美最黄视频在线播放免费| av欧美777| 国产精品亚洲av一区麻豆| 国产精品一区二区精品视频观看| 国产一区二区在线av高清观看| 日本一区二区免费在线视频| 国产精品1区2区在线观看.| 在线免费观看的www视频| 亚洲成a人片在线一区二区| 国产一级毛片七仙女欲春2| 亚洲av熟女| 黄色毛片三级朝国网站| www日本黄色视频网| tocl精华| 女人爽到高潮嗷嗷叫在线视频| 欧美黄色片欧美黄色片| 女人被狂操c到高潮| 我的老师免费观看完整版| 香蕉丝袜av| 91大片在线观看| 在线视频色国产色| 日本在线视频免费播放| 2021天堂中文幕一二区在线观| 欧美色视频一区免费| 白带黄色成豆腐渣| 国产精品av久久久久免费| 亚洲欧美日韩高清在线视频| 午夜a级毛片| 亚洲欧美激情综合另类| 老熟妇仑乱视频hdxx| 日韩高清综合在线| 欧美精品啪啪一区二区三区| 香蕉丝袜av| 国产1区2区3区精品| 国产av不卡久久| 99精品在免费线老司机午夜| 天堂√8在线中文| 99久久精品国产亚洲精品| 此物有八面人人有两片| tocl精华| 久久婷婷成人综合色麻豆| 午夜福利视频1000在线观看| 亚洲午夜理论影院| 亚洲精品粉嫩美女一区| 两个人看的免费小视频| av天堂在线播放| 在线观看免费午夜福利视频| 欧美精品啪啪一区二区三区| 亚洲18禁久久av| 黄色丝袜av网址大全| 久久久国产欧美日韩av| 久久九九热精品免费| 曰老女人黄片| av在线天堂中文字幕| 精品国产乱子伦一区二区三区| 搡老妇女老女人老熟妇| 久久久久久久久中文| 亚洲熟妇中文字幕五十中出| 在线播放国产精品三级| 日韩精品免费视频一区二区三区| 精品午夜福利视频在线观看一区| 无遮挡黄片免费观看| 久久久国产成人精品二区| 天天躁夜夜躁狠狠躁躁| 18禁裸乳无遮挡免费网站照片| 麻豆成人av在线观看| 亚洲精品久久国产高清桃花| 精品不卡国产一区二区三区| 亚洲av第一区精品v没综合| 欧美日韩国产亚洲二区| 国产蜜桃级精品一区二区三区| 日韩大码丰满熟妇| 一二三四在线观看免费中文在| 日本a在线网址| 在线十欧美十亚洲十日本专区| 黄色视频,在线免费观看| 看黄色毛片网站| 国产精品爽爽va在线观看网站| 高清在线国产一区| 免费看a级黄色片| 黄色片一级片一级黄色片| 一进一出抽搐gif免费好疼| 久久九九热精品免费| 精品一区二区三区视频在线观看免费| 婷婷精品国产亚洲av在线| 黄色视频不卡| 国产激情偷乱视频一区二区| 久久久水蜜桃国产精品网| 久久久精品大字幕| 两个人免费观看高清视频| 国产片内射在线| 一进一出好大好爽视频| 日韩欧美在线二视频| 日韩欧美精品v在线| 非洲黑人性xxxx精品又粗又长| 色哟哟哟哟哟哟| 成年人黄色毛片网站| 全区人妻精品视频| 俺也久久电影网| 香蕉久久夜色| 成年版毛片免费区| 床上黄色一级片| 久久午夜综合久久蜜桃| 亚洲片人在线观看| 免费看日本二区| 美女大奶头视频| 可以免费在线观看a视频的电影网站| 午夜免费成人在线视频| 国产麻豆成人av免费视频| 一进一出抽搐gif免费好疼| 最近在线观看免费完整版| ponron亚洲| 在线观看日韩欧美| 欧美一区二区精品小视频在线| 亚洲一卡2卡3卡4卡5卡精品中文| 两个人的视频大全免费| 亚洲欧美日韩高清专用| 成人特级黄色片久久久久久久| 男人的好看免费观看在线视频 | 深夜精品福利| 黑人欧美特级aaaaaa片| 国产成年人精品一区二区| 亚洲欧美日韩高清在线视频| 变态另类丝袜制服| 麻豆一二三区av精品| a级毛片a级免费在线| netflix在线观看网站| 丝袜美腿诱惑在线| 特级一级黄色大片| 亚洲七黄色美女视频| 亚洲av五月六月丁香网| 一边摸一边抽搐一进一小说| 亚洲熟女毛片儿| 日韩大码丰满熟妇| 亚洲中文字幕一区二区三区有码在线看 | 久久久久久久午夜电影| 人人妻人人看人人澡| 久久中文字幕人妻熟女| 欧美 亚洲 国产 日韩一| 一夜夜www| 国产精品国产高清国产av| 一级片免费观看大全| 亚洲人成网站在线播放欧美日韩| 亚洲一卡2卡3卡4卡5卡精品中文| 99国产极品粉嫩在线观看| 老司机福利观看| 久久精品aⅴ一区二区三区四区| 久久热在线av| 免费看美女性在线毛片视频| 精品国产亚洲在线| 可以免费在线观看a视频的电影网站| 99精品欧美一区二区三区四区| 一本大道久久a久久精品| 成年人黄色毛片网站| 国产在线观看jvid| 亚洲 欧美一区二区三区| 人成视频在线观看免费观看| 又粗又爽又猛毛片免费看| 老司机福利观看| 亚洲精品国产精品久久久不卡| 精品福利观看| 亚洲欧美激情综合另类| 日韩欧美三级三区| 不卡av一区二区三区| 一级作爱视频免费观看| 欧美黑人精品巨大| 男男h啪啪无遮挡| 一边摸一边抽搐一进一小说| 亚洲乱码一区二区免费版| 午夜久久久久精精品| 精品国产乱码久久久久久男人| 亚洲一区二区三区色噜噜| 成人亚洲精品av一区二区| 久久精品综合一区二区三区| 女同久久另类99精品国产91| 性欧美人与动物交配| 国产av不卡久久| 国产精品99久久99久久久不卡| 亚洲欧美激情综合另类| 99精品欧美一区二区三区四区| 欧美精品啪啪一区二区三区| 日韩有码中文字幕| 国产精品久久久人人做人人爽| 欧美乱码精品一区二区三区| 青草久久国产| 午夜福利欧美成人| 亚洲 欧美一区二区三区| a级毛片在线看网站| 午夜免费成人在线视频| 淫妇啪啪啪对白视频| 免费人成视频x8x8入口观看| 午夜福利在线在线| 黄色a级毛片大全视频| 人成视频在线观看免费观看| 人人妻,人人澡人人爽秒播| 亚洲人成网站高清观看| 亚洲aⅴ乱码一区二区在线播放 | 亚洲avbb在线观看| 露出奶头的视频| 动漫黄色视频在线观看| 亚洲精品国产一区二区精华液| 亚洲人成77777在线视频| 国产99白浆流出| 国产成人啪精品午夜网站| 叶爱在线成人免费视频播放| 每晚都被弄得嗷嗷叫到高潮| 欧洲精品卡2卡3卡4卡5卡区| 欧美日韩亚洲综合一区二区三区_| 两个人的视频大全免费| www.熟女人妻精品国产| 黄频高清免费视频| 好男人电影高清在线观看| 在线播放国产精品三级| 性色av乱码一区二区三区2| 亚洲全国av大片| 黑人欧美特级aaaaaa片| 一个人观看的视频www高清免费观看 | 一级黄色大片毛片| 两人在一起打扑克的视频| 床上黄色一级片| 亚洲人成网站在线播放欧美日韩| 精品午夜福利视频在线观看一区| 国产精品爽爽va在线观看网站| 欧美日韩黄片免|