混合顆粒系蒙特卡羅消光模型及反演方法

黃 茜, 蘇格毅, 孫存金, 鄧 飛, 陳 軍, 楊薈楠, 蘇明旭

上海理工大學(xué)能源與動(dòng)力工程學(xué)院, 上海 200093

引 言

消光光譜法[1]原理簡(jiǎn)單、 測(cè)試快捷、 結(jié)果準(zhǔn)確, 廣泛應(yīng)用于懸浮粉塵[2]、 火焰煙塵[3]、 濕蒸汽[4]、 乳劑[5]中顆粒粒徑和濃度分析。 近年來(lái), 很多學(xué)者對(duì)消光法顆粒測(cè)量進(jìn)行深入研究, Zhang等[6]研究顆粒物的折射率與波長(zhǎng)關(guān)系及對(duì)消光法的影響, 指出折射率的變化可導(dǎo)致消光法粒徑測(cè)量偏差, 處理消光譜數(shù)據(jù)時(shí), 應(yīng)考慮顆粒的色散特性; Tuersun等[7]關(guān)注到金納米球溶液中消光光譜對(duì)粒徑的敏感性, 研究了共振粒子的消光光譜法測(cè)量多分散金納米球的粒徑分布和濃度問(wèn)題; Krogs?e等[8]根據(jù)消光光學(xué)微粒計(jì)數(shù)器輸出信號(hào), 測(cè)得油液中顆粒污染物的真實(shí)粒度。 以上研究成功應(yīng)用至單一類型顆粒系的測(cè)量, 不過(guò)在現(xiàn)實(shí)中諸多被測(cè)對(duì)象包含了兩種甚至多種顆粒物, 稱為混合顆粒體系, 此時(shí)既有的Mie散射理論和Lambert-Beer定律構(gòu)建的消光模型將不再適用。

蒙特卡羅方法(Monte Carlo method, MCM)是一種概率統(tǒng)計(jì)方法, 已成功應(yīng)用于顆粒系中光復(fù)散射(多次散射)和顆粒兩相及多相流研究。 Dap等[9]發(fā)展了基于蒙特卡羅原理的消光光譜法顆粒粒徑和濃度的反演程序, 并通過(guò)可見光和紅外光譜實(shí)驗(yàn)數(shù)據(jù)進(jìn)行驗(yàn)證; Lebovka等[10]建立了蒙特卡羅數(shù)值模型, 模擬光束在圓柱體顆粒懸浮液中的傳遞過(guò)程, 研究顆粒等效截面及與Lambert-Beer定律的偏差, 探究碳納米管的光學(xué)性質(zhì); 肖新宇等[11]運(yùn)用蒙特卡羅方法研究了消光法顆粒粒徑測(cè)量中發(fā)散光束的影響, 并實(shí)現(xiàn)了粒徑反演修正。 鑒于蒙特卡羅法模擬具有過(guò)程清晰可追蹤、 接近實(shí)際的優(yōu)勢(shì), 作者建立基于蒙特卡羅原理的消光模型, 預(yù)測(cè)混合顆粒體系消光譜, 并研究顆粒類型、 混合比對(duì)消光特性的影響。 同時(shí)借助全局優(yōu)化算法, 實(shí)現(xiàn)混合顆粒系的顆粒粒徑和混合比同步反演。

1 原 理

1.1 消光光譜法原理

光譜消光法基于Lambert-Beer(LB)定律, 當(dāng)一束光強(qiáng)為I0的平行單色光, 穿過(guò)離散顆粒物兩相體系會(huì)發(fā)生散射和吸收并導(dǎo)致透射光的衰減, 將入射光強(qiáng)I0與透射光強(qiáng)I之比的對(duì)數(shù)形式稱為消光值, 可以描述顆粒物對(duì)某一波長(zhǎng)光吸收/散射的強(qiáng)弱與顆粒物濃度及光程的關(guān)系, 即

ln(I0/I)=πcnLR2kext=cnLCext

(1)

式(1)中,cn為顆粒數(shù)目濃度,R為被測(cè)顆粒半徑,L為光程,kext為消光系數(shù),Cext為消光截面。

LB定律適用于單一顆粒系的單分散情況, 在不同波長(zhǎng)條件下運(yùn)用式(1), 可獲得消光譜信息并確定顆粒系粒徑分布。 不過(guò), 對(duì)于混合顆粒體系, 該定律適用性遇到困難, 不同類型顆粒光散射特性不同, 其消光系數(shù)/截面也不一致, 在模型中很難表達(dá)并使得理論消光值無(wú)法計(jì)算。 此時(shí), 引入蒙特卡羅方法, 構(gòu)建新的消光譜預(yù)測(cè)模型。

1.2 蒙特卡羅方法

根據(jù)消光法原理, 蒙特卡羅方法核心思想是將物理上的入射光束按照“光子”概念抽象并做離散化處理, 通過(guò)光子與體系中顆粒的相互作用將光子歷程分為透射、 散射與吸收, 分類統(tǒng)計(jì)各去向的光子數(shù), 進(jìn)而得到其消光值。 根據(jù)其特點(diǎn), 易于引入到顆粒及混合顆粒多相體系的模擬計(jì)算。

如圖1所示, 兩相體系中混合了兩種球形顆粒: 黑色顆粒Ⅰ、 藍(lán)色顆粒Ⅱ。 幾何參數(shù)L為樣品池厚度(光程),dT、dR分別為發(fā)射器與接收器的直徑, 2H為樣品池上下邊界距離,S為發(fā)射器或者接收器與樣品池的距離,l為隨機(jī)散射自由程。 通過(guò)統(tǒng)計(jì)獲取接收器光子數(shù), 即可模擬在一定顆粒粒徑、 濃度和光波長(zhǎng)條件下混合顆粒系中光波行為并計(jì)算消光值。

圖1 光子在混合顆粒體系中傳播過(guò)程示意圖

圖2給出了蒙特卡羅模型計(jì)算流程, 根據(jù)輸入初始參數(shù)顆粒半徑R、 波長(zhǎng)λ、 體積濃度cv, 進(jìn)行消光值計(jì)算。 假設(shè)光束沿著圖1中x軸方向傳播, 平行光入射。 以發(fā)射器中心點(diǎn)為原點(diǎn), 光子的初始出射坐標(biāo)(x0,y0)

(2)

圖2 蒙特卡羅模型算法流程

式(2)中,ε1為[0, 1]區(qū)間內(nèi)服從均勻分布的隨機(jī)數(shù)。

光子剛進(jìn)入樣品池時(shí)坐標(biāo)為(S,y0), 若光子和顆粒發(fā)生碰撞, 則光子首次發(fā)生散射的坐標(biāo)(x1,y1)可以表示為

(3)

(4)

式中,εl是[0, 1]范圍內(nèi)均勻分布隨機(jī)數(shù),τ為濁度, 其在混合顆粒系中, 仍可定義為[12]

τ=cn×Cext

(5)

其中, 對(duì)于混合顆粒體系, 需要判斷顆粒類型

(6)

式(6)中,ε2是在[0, 1]區(qū)間內(nèi)服從均勻分布的隨機(jī)數(shù), 混合比φ指顆粒Ⅱ在整個(gè)混合顆粒體系中所占的數(shù)目比, 例如:φ=0時(shí)全為顆粒Ⅰ,φ=1時(shí)則全為顆粒Ⅱ。 當(dāng)前顆粒的消光截面Cext可由Mie散射理論計(jì)算[1]。 顆粒系數(shù)目濃度cn為

(7)

通常進(jìn)入兩相體系的光子可能被顆粒吸收、 被接收器直接接收(透射)、 未被接收器接收(逃逸)或在樣品池里再次散射, 定義散射截面和消光截面比值為反照率a

a=Csca/Cext

(8)

其中: 消光截面Cext和散射截面Csca由Mie散射理論計(jì)算得出。 根據(jù)設(shè)定條件判斷光子的下一事件

(9)

式(9)中,ε3是[0, 1]區(qū)間服從均勻分布的隨機(jī)數(shù),n為散射次數(shù)。

對(duì)于多次散射, 光子與顆粒碰撞后的空間散射角分布可以由Henyey-Greenstein相函數(shù)[13]確定, 散射角θ0的抽樣表示為

(10)

式(10)中,εθ0是另一個(gè)[0, 1]范圍內(nèi)隨機(jī)數(shù),g是不對(duì)稱因子。 發(fā)生第一次散射時(shí), 光子散射方向θ1=θ0。 重復(fù)式(4)—式(10), 根據(jù)散射自由程與散射角, 第n次散射時(shí)的光子散射方向?yàn)棣萵=θn-1+θ0(n>1), 光子散射后的坐標(biāo)為

(11)

根據(jù)式(9)條件判斷并統(tǒng)計(jì)其最終去向, 進(jìn)行下一個(gè)光子歷程, 直到完成所有光子計(jì)算, 統(tǒng)計(jì)各個(gè)物理過(guò)程的光子數(shù), 得到散射介質(zhì)的消光特性。 消光譜可表示為

ln(I0/I)=ln(Nset/N)

(12)

式(12)中,N為接收透射光子總數(shù),Nset為設(shè)定光子數(shù)。

2 數(shù)值結(jié)果

2.1 計(jì)算參數(shù)

按照前述模型編制計(jì)算程序, 模型尺寸(參見圖1)、 顆粒性質(zhì)、 可變參數(shù)均列于表1中, 介質(zhì)為水(折射率m=1.33), 顆粒Ⅰ選取聚苯乙烯(折射率m1), 顆粒Ⅱ選取高密度玻璃(折射率m2), 忽略顆粒的吸收特性。 根據(jù)前期初步數(shù)值結(jié)果分析, 設(shè)定光子數(shù)為105進(jìn)行計(jì)算。

2.2 光子去向統(tǒng)計(jì)

使用蒙特卡羅模型計(jì)算程序模擬R=0.2 μm的顆粒Ⅰ的光子去向, 將經(jīng)過(guò)一次與多次散射后被接收器接收到的光子分為單散射與復(fù)散射, 定義從前向邊界出射的其他逃逸光子為前向逃逸, 反之為后向。

圖3可以看出, 在給定粒徑下, 隨著波長(zhǎng)增大透射光子數(shù)逐漸增多。 按光散射理論, 波長(zhǎng)增大, 亞微米區(qū)顆粒的無(wú)因次參數(shù)α(α=2πR/λ)減小, 消光系數(shù)減小, 再結(jié)合式(4)計(jì)算出射光子隨機(jī)自由程, 光子準(zhǔn)直透射的概率增大。 在當(dāng)前光學(xué)厚度下, 散射光較微弱, 且以單散射為主, 其隨光波長(zhǎng)增大而減小。 由于接收器位置設(shè)定離樣品池較遠(yuǎn), 接收器尺寸小, 最終進(jìn)入接收器散射光有限(計(jì)算消光時(shí)限定只考慮透射光子), 逃逸發(fā)生的光子以前向?yàn)橹? 數(shù)目隨光波長(zhǎng)增大逐漸減小。

圖3 光子事件統(tǒng)計(jì)

2.3 蒙特卡羅方法驗(yàn)證

為了驗(yàn)證蒙特卡羅方法預(yù)測(cè)消光譜的準(zhǔn)確性, 首先分別考慮顆粒Ⅰ和顆粒Ⅱ的單一顆粒系情況, 采用蒙特卡羅方法仿真計(jì)算, 研究顆粒系的消光譜變化, 并和Lambert-Beer模型進(jìn)行對(duì)比。

圖4給出兩種顆粒各自的消光譜, 顆粒半徑分別為0.2、 0.8、 2.0 μm。 可以看出, 隨著顆粒粒徑的增大, 消光譜有明顯變化, 表明其對(duì)于粒徑是非常敏感, 有利于粒徑的反演; 當(dāng)粒徑超出亞微米區(qū)后(如2.0 μm), 消光系數(shù)隨著無(wú)因次參數(shù)α的增大而呈現(xiàn)振蕩現(xiàn)象; 而相同粒徑下, 由于顆粒Ⅰ和顆粒Ⅱ自身折射率不同, 導(dǎo)致消光譜趨勢(shì)不同。 通過(guò)結(jié)果對(duì)比看出蒙特卡羅方法預(yù)測(cè)與Lambert-Beer模型在數(shù)值上吻合, 相對(duì)誤差均在±2%以內(nèi), 表明蒙特卡羅方法適用于顆粒系的消光譜預(yù)測(cè)。

圖4 蒙特卡羅方法與Lambert-Beer模型的消光譜對(duì)比

單一顆粒系體積濃度cv與消光值存在線性關(guān)系[見式(1)], 其并不影響消光隨光波長(zhǎng)的變化趨勢(shì)即歸一化的消光譜, 對(duì)于混合顆粒系, 同樣通過(guò)數(shù)值模擬進(jìn)行了驗(yàn)證。 后續(xù)計(jì)算中, 僅以cv=2×10-5進(jìn)行計(jì)算, 不考慮cv對(duì)消光譜的影響。

2.4 混合顆粒系消光特性

將模型拓展到由顆粒Ⅰ與顆粒Ⅱ組成的混合顆粒系中,R1=0.2 μm不變,R2不同時(shí), 改變兩種顆粒的數(shù)目混合比φ, 進(jìn)行消光譜預(yù)測(cè)。

圖5(a)、 (b)、 (c)給出了不同混合比條件下顆粒系的消光譜。 隨著入射光波長(zhǎng)的增大, 消光譜呈遞減趨勢(shì), 顆粒Ⅱ的粒徑對(duì)消光譜曲線有明顯影響。 當(dāng)R2為0.15 μm, 圖5(a)中消光譜隨波長(zhǎng)的增大逐漸趨緩, 而當(dāng)R2增至0.3 μm, 由于消光系數(shù)差異, 圖5(c)中消光譜隨波長(zhǎng)增大趨勢(shì)更近于線性。 對(duì)于給定波長(zhǎng)和粒徑, 從圖5(d)可以看出: 顆粒系的消光值隨混合比增大而遞增, 而波長(zhǎng)減小時(shí), 消光值由線性趨勢(shì)向非線性趨勢(shì)轉(zhuǎn)變, 且兩種顆粒消光特性差別越大, 非線性趨勢(shì)越明顯, 因?yàn)閮煞N顆粒類型不同, 當(dāng)入射光波長(zhǎng)改變時(shí), 顆粒的消光特性相應(yīng)變化。 上述結(jié)果說(shuō)明, 在一定濃度下混合顆粒系的消光由顆粒類型、 混合比、 顆粒粒徑和光波長(zhǎng)等共同決定, 如混合比對(duì)消光譜的影響與不同類型顆粒的散射特性差異有關(guān), 其受制于粒徑和波長(zhǎng)取值, 二者同時(shí)影響顆粒的消光和散射效應(yīng), 據(jù)此可根據(jù)消光譜同步反演出顆粒粒徑與混合比。

圖5 混合顆粒系消光預(yù)測(cè)

3 反 演

3.1 反演方法

最優(yōu)化算法是消光譜反演求解顆粒粒徑時(shí)最為常見方法[14], 混合顆粒系可借助蒙特卡羅模型獲得理論消光譜, 不過(guò), 由于兩種甚至多種顆粒類型的存在, 增加了混合比等待定參數(shù), 反演難度會(huì)增加。 采用消光譜誤差構(gòu)建目標(biāo)函數(shù)

(13)

式(13)中: ln(I/I0)sim是理論消光譜, ln(I/I0)set為設(shè)定值的消光譜,j為計(jì)算消光譜的波長(zhǎng)數(shù),j=1, 2, …,M。 在設(shè)定反演參數(shù)的上下限范圍內(nèi)尋優(yōu)搜索最佳粒徑與混合比, 使目標(biāo)函數(shù)達(dá)到最小, 從而實(shí)現(xiàn)參數(shù)的反演。

為更好地驗(yàn)證反演效果, 采用三種全局最優(yōu)算法對(duì)所構(gòu)造的目標(biāo)函數(shù)進(jìn)行尋優(yōu), 分別是優(yōu)化遺傳算法(improved genetic algorithm, IGA)[15]、 粒子群算法(particle swarm optimization, PSO)[16]、 改進(jìn)差分進(jìn)化算法(improved differential evolution, IDE)[17]。 其中, IGA具有較好全局收斂性, 經(jīng)優(yōu)化設(shè)定最大迭代次數(shù)為500、 種群尺度為50、 交叉概率為0.8、 變異概率為0.045, 以提高算法穩(wěn)定性。 PSO采用自適應(yīng)慣性權(quán)重, 對(duì)適應(yīng)度值小于平均值的粒子進(jìn)行小波變異, 增強(qiáng)群體多樣性, 并使粒子在解空間的其他區(qū)域進(jìn)行搜索, 防止粒子群算法陷入局部最優(yōu), 設(shè)定種群尺度為500, 迭代次數(shù)在50～60之間。 IDE引入了自適應(yīng)控制參數(shù)因子對(duì)縮放因子和交叉因子進(jìn)行優(yōu)化, 選定交叉因子為0.85, 縮放因子上限和下限分別為0.5和0.3, 種群尺度為50, 并設(shè)定最大迭代次數(shù)為50, 以提高優(yōu)化效率。

如圖6所示, 進(jìn)行顆粒和混合比的單參數(shù)、 雙參數(shù)和三參數(shù)同步反演。 分別為: ①已知混合顆粒系粒徑, 對(duì)混合比進(jìn)行反演; ②已知顆粒系混合比, 對(duì)兩種顆粒粒徑進(jìn)行反演; ③兩種顆粒類型粒徑相同或不同時(shí), 對(duì)粒徑與混合比進(jìn)行同步反演。 分別選取R1=R2=0.2 μm、φ=0.8和R1=0.2 μm、R2=0.3 μm、φ=0.4兩種情形, 設(shè)置粒徑參數(shù)上下限為R∈[0.05 μm, 1 μm]、 混合比φ∈[0, 1], 計(jì)算波長(zhǎng)數(shù)M=9。

圖6 反演算例

3.2 反演結(jié)果分析

圖7(a)為已知顆粒粒徑R1和R2, 求解其混合比, 可以看出: 無(wú)論顆粒粒徑是否相同, 三種優(yōu)化算法反演混合比與設(shè)定值相對(duì)誤差均小于1.5%。 圖7(b)給定顆粒系的混合比φ=0.4, 反演兩種顆粒的粒徑, 為雙參數(shù)反演, 粒徑相對(duì)誤差在3%以內(nèi), 均較為精確。 為驗(yàn)證算法穩(wěn)定性, 以IGA為例, 進(jìn)行三次重復(fù)反演, 其粒徑R1、R2的標(biāo)準(zhǔn)差為8.2×10-4與1.1×10-3。 圖7(c)、 (d)為三參數(shù)同步反演,R1=R2時(shí),φ的反演誤差仍小于2%, 但R1反演誤差可至約8%(PSO算法); 當(dāng)R1≠R2時(shí), 三參數(shù)的反演誤差增大,R1、R2的誤差接近9%。

圖7 混合顆粒系參數(shù)反演結(jié)果

此外, 三種算法的反演時(shí)間不盡相同, 雙參數(shù)反演時(shí), IGA和IDE的一次計(jì)算時(shí)間約需40 min, 而PSO反演耗時(shí)為其數(shù)倍。 總體上, IGA算法效果較好且更加穩(wěn)定, 就本文算例來(lái)看, 所有算例下反演誤差均在10%以內(nèi)。

4 結(jié) 論

研究了兩相顆粒系光散射效應(yīng), 建立了蒙特卡羅原理的混合顆粒體系消光預(yù)測(cè)模型, 獲得聚苯乙烯和高密度玻璃微珠混合顆粒系的消光譜, 據(jù)此討論了混合比、 混合顆粒粒徑對(duì)消光譜的影響。 根據(jù)預(yù)測(cè)消光譜對(duì)顆粒粒徑與混合比同步反演。 結(jié)果表明:

(1)在單一顆粒系中, 蒙特卡羅方法與傳統(tǒng)Lambert-Beer模型的計(jì)算結(jié)果較一致, 相對(duì)誤差均在±2%以內(nèi)。

(2)對(duì)混合顆粒系的消光譜預(yù)測(cè), 不同混合比的消光譜均分布在φ=0與φ=1之間, 隨顆粒Ⅱ粒徑不同, 消光譜趨勢(shì)相應(yīng)變化, 顆粒Ⅱ的占比對(duì)消光譜有明顯影響, 消光值隨著混合比遞增, 當(dāng)波長(zhǎng)減小時(shí)增長(zhǎng)趨勢(shì)由線性向非線性轉(zhuǎn)變, 即不同顆粒類型影響了消光譜的數(shù)值大小與趨勢(shì)。

(3)對(duì)混合顆粒系的消光譜反演, IGA、 PSO、 IDE算法均可實(shí)現(xiàn)單參數(shù)和多參數(shù)同步反演。 其中, 僅對(duì)混合比參數(shù)反演時(shí)誤差最小, 相對(duì)誤差在1.5%以內(nèi), 隨反演參數(shù)增加, 誤差增大, 三參數(shù)反演時(shí), 混合顆粒的粒徑變化對(duì)反演有影響, 但算例中誤差均在10%以內(nèi)。