基于MP-PIC方法的冶金硅氫氯化流化床反應器模擬

吳詩鳴，陳皓寧，宗原，許志美，趙玲

（華東理工大學化學工程聯合國家重點實驗室，上海 200237）

引 言

多晶硅是目前光伏能源以及芯片制造行業中不可或缺的原材料[1-3]。在目前的主流工藝中，如西門子法[4]，通常會副產大量的四氯化硅（STC）[5-6]。冶金硅氫氯化反應通過將四氯化硅與冶金硅以及氫氣反應可以實現四氯化硅的循環利用[7]。其總反應方程式為

反應為可逆弱放熱反應[8]。工藝中反應溫度通常在500～600℃之間。根據Mui[9]的研究，反應由以下兩步可逆反應串聯構成。

其中第一步反應是氣相均相反應，反應速率較慢，是整個反應過程的速率控制步驟。第二步反應是氣固非均相反應，反應速率較快。丁偉杰[10]的研究表明第一步反應在600℃以下受到動力學限制，冶金硅表面的雜質對氣相反應具有催化作用。因此，在工藝溫度下，氫氯化反應主要發生在冶金硅顆粒表面。Sill[11]的研究表明催化劑對于冶金硅氫氯化反應的進行至關重要。相比冶金硅顆粒表面的雜質鐵元素，添加金屬銅催化劑可以明顯加快反應。Becker[12]認為金屬銅對反應的催化過程與冶金硅顆粒表面固有的鐵元素催化過程是并行的，并提出了相應的動力學模型。

冶金硅氫氯化反應通常在流化床反應器中進行。流化床反應器具有良好的傳熱傳質特性，同時操作壓降小，可以實現反應過程的連續化操作[13]。然而，流化床反應器內復雜的流動與反應過程使得反應器放大較為困難。Colomb 等[14]采用兩相模型[15]對冶金硅氫氯化實驗流化床反應器進行模擬，計算結果與Sill 等[11,16]的實驗結果吻合較好。由于兩相模型中對氣泡相以及乳化相的運動進行了簡化處理，其中的流動模型不能真實反映反應器中的流場。因此，由反應器放大引起的流場變化需在模型中進行修正[17]。

相比經典的兩相模型，計算流體力學（CFD）方法可以模擬反應器內流場以及組分分布情況，已經成為指導反應器放大與優化的重要方法[18]。Liu等[19]采用雙流體模型與群體平衡方法耦合計算的方式對硅烷熱裂解流化床反應器進行模擬，探究了影響硅顆粒生長的內在機制以及操作條件對生長速率的影響。此外，基于該模型又進一步探究了動力學模型對于模擬結果的影響[20]。Li 等[21]采用雙流體模型對生產多晶硅顆粒的內循環流化床反應器進行了模擬，探究了操作氣速對流化床在傳熱性能以及抑制壁面沉積等方面的影響，并給出了氣速的最優操作范圍。 本文結合MP-PIC（multiphaseparticle in cell）方法對小試冶金硅氫氯化反應器進行數值模擬。MP-PIC 方法是一種面向工業級流化床模擬的計算方法。該方法最早由Andrews 等[22]提出，由Snider 等[23-25]發展成熟并開發出商業化軟件。MP-PIC 模型中對流體相采用歐拉方法描述；對于顆粒相采用計算顆粒（parcel）代替真實顆粒群，同時采用經驗的應力模型簡化顆粒間的相互作用，從而大幅減小了計算量[26]。

本文在Sill 等[11,16]的實驗基礎上對實驗級流化床反應器進行模擬。根據反應動力學模型的特點，在網格以及計算顆粒包含的真實顆粒數量無關性的基礎上探究了反應器內初始三氯氫硅（TCS）濃度對模擬結果的影響。基于模擬與實驗驗證結果，進—步對反應器內的反應與流動進行了分析。對不同壓力以及進口H2∕STC 摩爾比操作條件下的TCS產率進行了研究。研究結果可以為反應器的放大與優化提供參考。

1 數值方法

1.1 MP-PIC模型

1.1.1 氣相方程 氣相質量守恒方程為

式中，εg為氣相體積分數；ρg為氣體密度；Ugj為氣體速度在j方向的分量；Ng為氣體組分數；Rgn為氣體組分n的消耗或者生成速率。

氣相動量守恒方程為

式中，Pg是氣體壓力；τgij為氣相應力張量；gi為重力加速度；Ipg為氣固相相互作用動量源項；氣相密度ρg由理想氣體方程求解；μg,eff為有效氣體黏度，由式（4）計算得到；μg為氣體分子黏度，由Sutherland方程[27]計算得到；湍流黏度μe采用標準k-eplison 模型[28]進行求解。

氣相組分守恒方程為

式中，Xgn表示氣相組分n的質量分數；Rgn為化學反應源項；Dgn表示氣相組分n的分子擴散系數，由FSG 關聯式[29]計算得到；Dt為湍流擴散系數，由式（6）計算[30]。

式中，Sct為湍流Schmidt數，取0.9。

氣相能量守恒方程為

式中，Tg為氣相溫度；CPg、CPgn分別為氣相混合比熱容以及氣相組分n的比熱容；λg為分子傳熱系數；λt為湍流傳熱系數，由式（9）進行計算[30]，湍流Prandtl 數Prt取0.9；QD為氣相組分焓擴散項；Spg為顆粒對流體的能量源項。式（7）中忽略了流體黏性耗散項，并且由于氫氯化反應速率較慢以及反應放熱較小，床內氣固相間溫差較小，因此也忽略了氣固相間的輻射傳熱。

1.1.2 顆粒相方程 顆粒的質量守恒方程為

式中，Np為真實顆粒中的化學組分數；Wp為顆粒統計權重，即每個計算顆粒中包含的真實顆粒數；Rpn為真實顆粒中組分n的生成或者消耗速率。

顆粒的組分守恒方程為

式中，Ui為計算顆粒在i方向的速度分量；τp為顆粒應力，由式（13）計算得到；Igp是流體對計算顆粒的動量傳遞項。

式中，εcp為固相堆積體積分數；εs為顆粒相體積分數；Ps為經驗壓力常數；β為無量綱指數；α為一非常小的常數，用于防止分母為零[30]。

顆粒空間位置的變化由式（14）得到。

式中，xi代表計算顆粒在坐標軸i方向的坐標。顆粒的能量守恒方程為

式中，Cp、T分別代表計算顆粒的比熱容和溫度；hpn為組分n的比焓；Sgp為流體對顆粒的能量源項。

1.2 相間作用力模型

流體對計算顆粒的動量交換項由式（16）定義。

1.3 相間能量傳遞模型

流體對計算顆粒之間的能量傳遞源項由簡化表達式（17）定義[30]。

式中，γcp為對流傳熱系數；As為計算顆粒k中真實顆粒的表面積。對流傳熱系數由式（18）計算。

式中，Nu為Nusselt數；dp為真實顆粒的直徑；λg為流體熱導率。Nu由Ranz & Marshall 關聯式[32]進行計算。

1.4 化學反應模型

根據丁偉杰[10]的研究，反應主要發生在顆粒附近，因此可將總反應視為氣固反應，反應速率由總包動力學表達式[式（19）]進行計算[9,14]。其中的參數列于表1。模擬中假設反應器內催化劑分布均勻。

表1 反應動力學參數［9，14］Table 1 Parameters of kinetic equation［9，14］

式中，pSiHCl3、pSiCl4、pH2、patm分別為各氣相組分分壓以及標準大氣壓。

1.5 計算方法

基于Sill 等[11,16]的實驗，對實驗流化床反應器進行模擬。鼓泡床狀態下操作，如圖1所示，流化床高度為0.400 m，內徑為0.050 m，初始床層高度為0.165 m。數值計算過程在開源軟件MFIX 中完成。所采用的實驗數據列于表2。模擬的初始條件按實驗條件或者根據分析的需要進行設置。對于邊界條件，在流化床底部進口處設置為質量進口條件，流化床頂部出口設置為壓力出口條件。對于壁面條件，氣相采用壁面無滑移條件，顆粒相與壁面作用采用不完全彈性碰撞。壁面絕熱，沒有質量流入。相應的模擬參數列于表3。

表2 實驗數據Table 2 Experimental data

表3 模擬參數Table 3 Parameters for simulation

圖1 流化床反應器示意圖Fig.1 Sketch of fluidized bed reactor

2 計算結果與分析

2.1 無關性分析

2.1.1 網格以及計算粒子所包含的顆粒數無關性

基于表2 中第2 組實驗條件，在不考慮反應、傳熱和傳質的前提下，對網格以及計算顆粒包含的真實顆粒數量進行無關性檢驗以確定相應的合理值。采用cell-cut 方法[33]對圓柱進行劃分。采用15×115×15、20×155×20、25×190×25三種方式進行劃分，相應的網格大小分別為3.3 mm×3.5 mm×3.3 mm、2.5 mm×2.6 mm×2.5 mm、2.0 mm×2.1 mm×2.0 mm。圖2給出了計算Wp為200 時，采用三種不同網格劃分方式計算得到的時均氣相相分率軸向分布情況。對比可知，三種網格劃分方式得到的計算結果較為接近。基于20×155×20 網格劃分方式，對Wp值的大小進行檢驗。圖3 給出了Wp值分別為100、200、400時，計算得到的時均氣相相分率軸向分布情況。結果表明，Wp值變化對計算結果的影響可以忽略。綜合考慮計算量以及計算精度，選擇采用20×155×20的網格劃分方式，Wp值設置為200。

圖2 采用不同網格劃分方式得到的時均氣相相分率軸向分布Fig.2 Axial distribution of time-averaged gas phase fraction for different grids

圖3 采用不同計算顆粒包含的真實顆粒數量得到的時均氣相相分率軸向分布Fig.3 Axial distribution of time-averaged gas phase fraction for different number of particle per parcel

2.1.2 TCS初始濃度無關性 由于模擬中采用總包反應動力學對反應器進行模擬，反應動力學模型的計算要求必須設置反應器內三氯氫硅（TCS）初始濃度。為了考察TCS 初始濃度對于計算結果的影響，基于表2中1號實驗數據，采用表4所示三組不同的氣相組分濃度作為初始值進行計算。圖4給出了不同TCS初始濃度計算得到的出口TCS摩爾分數隨流動時間變化情況。從圖中可知，雖然初始濃度不同，隨著流動進行，出口TCS摩爾分數趨于一致。圖5 給出了30 s 時反應器內軸向平均溫度分布情況。結果表明，不同TCS 初始濃度計算得到的反應器內軸向溫度分布較為接近。TCS初始濃度無關性檢驗表明，當TCS初始濃度較高時，反應逆向進行。雖然此時反應吸熱，但是由于吸熱量較少，溫度變化較小。這使得初始濃度下的溫度分布不會對之后的反應造成較大影響。反之，當TCS初始濃度較低時，反應正向進行引起的溫度升高量也同樣較小，進而對后續反應的影響同樣較小。

圖4 不同TCS初始濃度條件下出口TCS摩爾分數隨流動時間變化情況Fig.4 Change of outlet TCS molar fraction with flow time under different initial concentration of TCS

圖5 不同TCS初始濃度條件下反應器內軸向平均溫度分布情況Fig.5 Axial average temperature distribution in the reactor under different initial concentration of TCS

表4 TCS初始摩爾分數Table 4 Initial mole fraction of TCS

2.2 模擬結果驗證

基于以上無關性分析，對表2 中的實驗條件進行模擬。TCS產率由式（22）計算。

式 中，n?SiCl4,IN和n?SiHCl3,OUT分 別 為 反 應 器 進 口 處STC 摩爾流率以及出口處TCS 摩爾流率。圖6 給出了TCS 產率的模擬結果與實驗值之間的對比，包括Colomb等[14]對該反應器的模擬結果。模擬值與實驗值的相對誤差由表5 給出。結果表明，與Colomb等[14]的模擬結果相比，本文采用的計算模型同樣具有相當的精確度。其中，第1、2、5、6 組計算得到的TCS 產率與實驗值更為接近；而在第3、4、7、8 組中，計算得到的產率稍高。導致這一現象的原因可能是模擬中采用的催化劑均勻分布假設在較低進口表觀氣速時與實驗情況較為符合，在較高的進口表觀氣速下與實驗情況存在一定差異。但是，這種差異在本文的研究范圍內是可以接受的。

圖6 計算得到的TCS產率與文獻模擬值[14]以及實驗值對比Fig.6 Comparison of simulated TCS yield with simulated results[14]and experimental results

表5 模擬值與實驗值的相對誤差Table 5 Relative error between simulated results and experimental results

2.3 反應與流動特征分析

基于表2 中第1、2、4、6 組實驗條件模擬結果對反應過程中流動與反應特征進行分析。圖7 和圖8分別給出了相應的軸向平均表觀氣速以及壓力分布情況。從圖中可見，雖然反應會使得氣相摩爾流率減小，但是軸向平均表觀氣速分布仍較為均勻。這是由于反應轉化率較低。按照實驗轉化率計算，氣體體積減小的比例僅為2%左右，而由于床層壓降引起的氣相壓力減小所導致的氣相體積膨脹可以彌補由反應減小的氣相體積。因此，反應對軸向平均表觀氣速的影響可以忽略，軸向平均表觀氣速可以視為不變。

圖7 30 s時反應器內軸向平均表觀氣速分布Fig.7 Axial average apparent gas velocity distribution in the reactor at 30 s

圖8 30 s時反應器內軸向平均壓力分布Fig.8 Average axial pressure distribution in the reactor at 30 s

已知反應受平衡限制，引入Qp∕Kp值作為判斷反應遠離平衡程度的指標，對不同操作條件下反應器內反應特征進行分析。圖9 和圖10 分別給出了不同條件下，時均Qp∕Kp值分布情況以及對應的氣相相分率分布情況。對比圖9（a）、（b）可知，在較低的進口表觀氣速下，氣相組分沿徑向分布較為均勻，主要表現為平推流。同時，從圖10中可以發現，在較低的表觀氣速下，氣固相分布較為均勻[圖10（a）]。進口表觀氣速增加對床內流動的影響主要表現在顆粒床層高度略增加，而對床內氣固相間空間分布的均勻性影響較小。相比床層高度增加，進口表觀氣速的增加更加明顯，由此導致反應氣體在床內的停留時間縮短。同時，由速率常數k的量綱可知，反應物的反應級數大于1，由較高的表觀氣速引起的氣相組分返混會使得TCS的產率進一步降低。

圖9 流化床反應器中時均Qp∕Kp值分布情況Fig.9 Distribution of time-averaged Qp∕Kp value in the simulated fluidized bed reactor

圖10 模擬流化床反應器中時均氣相相分率分布情況Fig.10 Distribution of time-averaged gas volume fraction in simulated fluidized bed reactor

對比圖9（b）～（d）可知，升高反應溫度或者增加催化劑量均可以加快反應速率從而提升轉化率。從圖9（c）、（d）中可以發現，反應在氣相反應物進入床層一段距離后已經接近平衡，這使得在床層上部，反應幾乎不再進行。因此，由表觀氣速引起的短路以及返混的影響可以忽略。這里值得說明的是，與添加催化劑不同，升高反應溫度會提高反應平衡常數，但是由于溫度對反應速率的影響更大，因此升高溫度會使得反應過程更接近平衡。由圖10（b）～（d）可以發現，由于反應轉化率受到限制，在反應消耗以及壓差共同作用下，反應器內表觀氣速相差較少，因而增加反應速率對于床內的相分率分布的均勻性影響很小。

綜合以上討論，對于實驗級的冶金硅氫氯化反應器而言，化學平衡與氣固流動是反應產率的共同控制因素。當反應器內反應離化學平衡較遠時，反應轉化率由流動控制，在本文的研究中，流動因素主要表現為氣相在床層中的停留時間；而當反應器內的反應接近平衡時，反應轉化率由化學平衡控制。在反應器設計與優化過程中應該綜合考慮這兩種因素的影響。

2.4 操作條件分析

基于以上分析，在Sill 等[11,16]對冶金硅流化床反應器的實驗研究基礎上，進一步對操作壓力以及進口H2∕STC摩爾比對TCS產率的影響進行探究。

2.4.1 壓力對TCS 產率的影響 圖11 中給出了不同操作壓力下TCS產率的變化情況。從圖中可以發現，在實驗級流化床反應器中，壓力對于產率的影響較小。從動力學表達式中可以發現，壓力雖然可以改變平衡常數大小，但是不改變Qp∕Kp值。在相同的條件下，反應速率的增大倍數為壓力增加倍數的1.086 次方。對于較低表觀氣速操作而言，反應已經接近平衡，增加操作壓力，對反應產率的提高不明顯。對于較高表觀氣速操作而言，隨著壓力增加，TCS產率逐漸接近較低表觀氣速下的產率。

圖11 不同進口表觀氣速下操作壓力對TCS產率的影響Fig.11 Effect of operation pressure on TCS yield under different inlet apparent gas velocity

2.4.2 進口H2∕STC 摩爾比對TCS 產率的影響 圖12 中給出了TCS 產率隨H2∕STC 摩爾比變化情況。由Le Chatelier 原理可知，對于可逆反應而言，當反應平衡時，提高一種反應物的濃度可以增加另一種反應物的轉化率。因此，對于受化學平衡限制的冶金硅氫氯化反應器而言，增加進口H2∕STC 摩爾比，可提高TCS產率。

此外，從圖12 中可以觀察到，隨進口H2∕STC 摩爾比變化，進口氣速對于TCS產率的影響發生變化。在進口H2∕STC 摩爾比較低時（H2∕STC 摩爾比＜2），反應可快速接近平衡，因此，流動因素對于TCS產率的影響不大。但隨著進口H2∕STC 摩爾比的增加，反應逐漸遠離平衡，反應轉化率逐漸轉為流動控制，停留時間成為主要的控制因素。因此，在相同的進口表觀氣速下，TCS 產率增加速率減緩。需要指出的是，雖然提高進口H2∕STC 摩爾比有利于TCS 產率的增加，但卻降低了H2的利用率，需要進行全局考慮。

圖12 不同進口表觀氣速下進口H2∕STC摩爾比對TCS產率的影響Fig.12 Effect of inlet H2∕STC molar ratio on TCS yield at different inlet apparent gas velocity

3 結 論

本文基于MP-PIC 方法采用總包反應動力學模型對實驗級冶金硅氫氯化反應器進行模擬。對TCS初始濃度進行無關性檢驗，結果表明由于反應弱放熱且轉化率較低，TCS 初始濃度對于反應器模擬結果幾乎沒有影響。根據實驗條件模擬得到TCS產率與實驗值吻合良好。對反應器內流動與反應特征以及不同操作條件下的計算結果進行分析得到以下結論。

（1）雖然在冶金硅氫氯化反應中氣相物質的量會減少，但是由于反應器內反應物轉化率較低，由床層壓降引起的氣相壓力降低可以彌補由氣相分子數減少引起的氣相體積變化。

（2）在反應器內，反應過程同時受到流動因素與化學平衡因素控制。當反應尚未接近平衡時，流動因素對TCS 產率起控制作用；當反應已經接近平衡時，流動因素對反應產率不起作用，轉而由化學平衡控制。

（3）提高操作壓力雖然可以加快反應速率，但是對于TCS 產率的影響不明顯。通過增加進口H2∕STC 摩爾比可以使得反應遠離平衡進而提高TCS 產率，但是會降低H2的利用率。

符 號 說 明

a,b,c——反應動力學參數

EA——反應活化能，kJ∕mol

H——床高，m

Kp——化學平衡常數

k——速率常數，mol∕(kg·s·Pa1.086)

mSi——冶金硅顆粒加載量，g

P——氣相總壓，Pa

Qp——反應熵

R——氣體常數，J∕(mol·K)

r——流化床半徑，m

rp——反應速率，mol∕(kg·s)

T——反應溫度，K

u0——進口表觀氣速，cm∕s

WCu——催化劑質量占冶金硅顆粒總質量百分比，%

Wp——每個計算顆粒中包含的真實顆粒數

YTCS,EXP——三氯氫硅產率，%