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

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

引 言

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

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

1 原 理

1.1 消光光譜法原理

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

ln(I0/I)=πcnLR2kext=cnLCext

(1)

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

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

1.2 蒙特卡羅方法

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

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

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

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

(2)

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

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

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

(3)

(4)

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

τ=cn×Cext

(5)

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

(6)

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

(7)

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

a=Csca/Cext

(8)

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

(9)

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

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

(10)

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

(11)

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

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

(12)

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

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

2.1 計算參數(shù)

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

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

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

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

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

2.3 蒙特卡羅方法驗證

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

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

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

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

2.4 混合顆粒系消光特性

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

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

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

3 反 演

3.1 反演方法

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

(13)

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

為更好地驗證反演效果, 采用三種全局最優(yōu)算法對所構(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)重, 對適應(yīng)度值小于平均值的粒子進(jìn)行小波變異, 增強群體多樣性, 并使粒子在解空間的其他區(qū)域進(jìn)行搜索, 防止粒子群算法陷入局部最優(yōu), 設(shè)定種群尺度為500, 迭代次數(shù)在50～60之間。 IDE引入了自適應(yīng)控制參數(shù)因子對縮放因子和交叉因子進(jìn)行優(yōu)化, 選定交叉因子為0.85, 縮放因子上限和下限分別為0.5和0.3, 種群尺度為50, 并設(shè)定最大迭代次數(shù)為50, 以提高優(yōu)化效率。

如圖6所示, 進(jìn)行顆粒和混合比的單參數(shù)、 雙參數(shù)和三參數(shù)同步反演。 分別為: ①已知混合顆粒系粒徑, 對混合比進(jìn)行反演; ②已知顆粒系混合比, 對兩種顆粒粒徑進(jìn)行反演; ③兩種顆粒類型粒徑相同或不同時, 對粒徑與混合比進(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], 計算波長數(shù)M=9。

圖6 反演算例

3.2 反演結(jié)果分析

圖7(a)為已知顆粒粒徑R1和R2, 求解其混合比, 可以看出: 無論顆粒粒徑是否相同, 三種優(yōu)化算法反演混合比與設(shè)定值相對誤差均小于1.5%。 圖7(b)給定顆粒系的混合比φ=0.4, 反演兩種顆粒的粒徑, 為雙參數(shù)反演, 粒徑相對誤差在3%以內(nèi), 均較為精確。 為驗證算法穩(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時,φ的反演誤差仍小于2%, 但R1反演誤差可至約8%(PSO算法); 當(dāng)R1≠R2時, 三參數(shù)的反演誤差增大,R1、R2的誤差接近9%。

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

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

4 結(jié) 論

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

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

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

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