間歇冷卻結晶過程的CFD-PBM數值模擬

覃 渫 蘭, 張 磊, 都 健

(大連理工大學 化工學院 化工系統工程研究所, 遼寧 大連 116024)

0 引 言

溶液結晶是化學、制藥和食品工業中應用較為廣泛的分離和純化技術之一[1].在結晶過程中,由于對物理化學和機械特性的強烈依賴性,需要控制結晶產品質量[2].粒度分布(crystal size distribution,CSD)是評價最終晶體質量的重要指標之一[3-4].

自20世紀60年代中期以來,粒數衡算模型(population balance model,PBM)已被廣泛用于結晶過程的模擬[1-2].粒數衡算模型研究的是既包括外部坐標(空間)又包括內部坐標(顆粒屬性)的分散系統[1].離散方法是應用較為廣泛的數值方法,該方法首先將整個連續范圍劃分為若干個相鄰的小區間,然后將粒數衡算方程轉換為一定數量的離散方程進行求解[2].盡管該方法計算成本較高,但其優點在于數值魯棒性好[5].

實際結晶系統多為固液兩相流體系,流場分布和結晶動力學等各個過程共同決定最終產品質量,進而影響晶體粒度分布[6-7].計算機模擬是研究結晶過程的一種有效方法,將計算流體動力學(computational fluid dynamics,CFD)和粒數衡算模型相結合,以解釋流場和物質濃度場中強烈的空間不均勻性[2].

關于CFD-PBM耦合模型描述結晶過程的研究已取得了一定成果.2018年,Fu等[8]建立了CFD-PBM耦合模型來模擬對乙酰氨基酚的連續冷卻結晶過程,結合生長動力學和成核動力學,探究混合效應和壁面溫度對顆粒粒度分布的影響.2019年,Mousavi等[9]對間歇式攪拌結晶器中磷酸銨鎂晶體的結晶過程,建立了CFD-PBM耦合模型,僅考慮晶體生長動力學,模擬計算的磷酸銨鎂晶體的粒度分布和實驗值相差不大.2021年,De Souza等[6]通過CFD-PBM耦合模型模擬了間歇式攪拌結晶器中硫酸鋁鉀的冷卻結晶過程,僅考慮晶體生長動力學,探究實際流動條件對晶體生長動力學的影響.

關于CFD-PBM耦合模型模擬結晶過程的研究中,大多對結晶動力學模型進行簡化,關于破碎過程的研究工作有限.然而,結晶過程是一個復雜的動態相平衡過程,晶體的破碎行為時有發生,特別是在攪拌結晶器中,顆粒的破碎會影響離散相的混合和流體流動行為[10].為了準確分析結晶過程,需要更全面地考慮晶體行為[11].

綜上所述,本文針對攪拌結晶器中撲熱息痛-乙醇體系的間歇冷卻結晶過程,采用CFD-PBM耦合模型進行模擬與優化,充分考慮生長、初級成核和二次成核過程,同時采用Luo破碎模型補充結晶動力學模型,以探究晶種加入、破碎過程、攪拌速度和降溫速率對晶體粒度分布的影響.

1 數學模型

1.1 粒數衡算模型

粒數衡算方程(population balance equation,PBE)描述了多相體系中離散相狀態在時間和物理空間上的分布[12-13].基于顆粒特征尺寸L的粒數衡算模型的表達式如下所示:

?[Γeff?n(L,y,t)]+ρl?[G(L)n(L,y,t)]/?L=

(1)

晶體的破碎速率為

g(L)n(L,t)

(2)

式中:g為破碎頻率;β(L,L′)為從顆粒特征尺寸L′到L的破碎顆粒概率密度函數;P為一個顆粒破碎產生子顆粒的數量.

1.2 CFD-PBM耦合模型

CFD是一種利用流體基本控制方程求解流體流動的數值方法[14].CFD-PBM耦合模型結合計算流體力學模型和粒數衡算模型,即將空間尺度的流場分布和時間尺度的顆粒狀態演變相結合,從而得到產品晶體的粒度分布.CFD模擬計算得到的顆粒流動特性和流場分布,用于求解PBM的成核速率、生長速率和破碎速率等,同時求解PBM方程得到顆粒粒度分布和顆粒直徑,進而修正和計算相間作用力和湍動能等,改進CFD模型的預測能力,因此CFD-PBM耦合是雙向的[7].

為了將離散相的PBM與CFD結合,采用Sauter平均直徑表示離散相的顆粒直徑.對于離散方法,Sauter平均直徑定義為

(3)

如圖1所示,在本文研究中,CFD-PBM耦合模型主要考慮并結合以下建模層次[7,15]:

圖1 CFD-PBM耦合模型

(1)計算流體動力學模型,描述顆粒在結晶器中的運動和分布.

(2)粒數衡算模型,描述所有顆粒相關性質在時間和空間上的演化.

(3)結晶動力學模型,描述結晶過程中晶體的行為.

2 結晶動力學模型

2.1 結晶過程

(1)晶體成核過程

成核速率是指單位時間單位體積溶液中新生成的微小晶體粒子數目[16].成核過程通常分為初級成核(primary nucleation)和二次成核(secondary nucleation)[3,17].

初級成核是在無晶種加入時自發成核[3,17].初級成核通常被經典成核理論(classical nucleation theory,CNT)所描述,初級成核速率b1的表達式如下[18]:

(4)

式中:kb1為初級成核常數;V為一個溶質分子的體積;k為玻爾茲曼常數;σ為晶體和溶液的界面能;T為溫度;s為相對過飽和度,s=c/cs,c為撲熱息痛在乙醇中的溶解度,cs為撲熱息痛在乙醇中的飽和溶解度.

二次成核是在擁有晶體時形成晶體[3,17].結晶器中二次成核的主要原因是晶體與攪拌器的接觸和碰撞,碰撞的概率與攪拌器的轉速成正比[18].二次成核速率的經驗方程為[19]

(5)

式中:kb2為二次成核常數;ms為晶體質量;α為常數2.08;β為常數0.713.

(2)晶體生長過程

晶體生長的機理有表面能理論和吸附層理論等[3].由于晶體生長的復雜性,至今仍未建立通用的生長理論[3].在工業結晶中,生長速率常采用冪函數型的經驗方程,表達式如下[3,19]:

G=kgexp(-Ea/RT)(Δc)γ

(6)

式中:kg為生長速率常數;Ea為晶體活化能;Δc為絕對過飽和度,Δc=c-cs;γ為動力學常數;R為摩爾氣體常數.

晶體破碎過程與顆粒性質和擾動情況有關[20].Luo破碎模型是常用的破碎模型,適用于湍流擴散中顆粒破碎的預測.Luo破碎模型是基于各向同性和概率理論建立的,可以預測給定顆粒尺寸的破碎速率.在湍流中,由于渦流在顆粒表面的轟擊,離散相表面產生相對速度的波動,當渦流的能量達到某個閾值時,可以滿足破碎所需要的表面能量,從而產生破碎過程[21].Luo破碎模型進行了一些假設[21]:湍流被認為是各向同性的;只考慮顆粒的二元破碎;破碎體積分數是一個隨機變量;破碎的發生是由渦旋能量等級決定的;只有長度小于或等于顆粒直徑的渦流才能引起粒子振蕩.

Luo破碎模型的破碎速率為

Bbre-Dbre=Ωbre(L,L′)=(0.923 8ε1/3L-2/3ω)×

2.047-1}ξ-11/3}dξ

(7)

式中:λ為渦流大小;ξ為量綱一渦流大小,ξ=λ/L;ε為離散相的局部分數;f為破碎常數;ρp為顆粒密度;ω為常數1.5.

2.2 撲熱息痛冷卻結晶的動力學模型

Li等[19]利用實驗數據回歸了撲熱息痛-乙醇結晶過程的生長和成核動力學.結晶實驗的裝置如圖2所示,其中結晶系統包括一個溫度探頭、攪拌系統、1 L玻璃結晶器和加熱冷卻金屬夾套.攪拌器為4個槳葉攪拌器.

圖2 實驗裝置示意圖[19]

實驗中,初始溶液體積為0.5 L,加入的晶種質量為0.03 g,晶種粒徑為75～106 μm,結晶器的攪拌速度為400 r/min,采用式(8)的線性降溫曲線.表1為涉及的物性參數.

表1 物性參數表

(8)

撲熱息痛在乙醇中的飽和溶解度為

在大一新生入學時，重視宣傳學科競賽的重要性，組織大一新生赴現場觀摩賽事，通過課堂表現和課外接觸發現可造之材，鼓勵有理想有潛力的新生參與到賽事之中，感受賽事的緊張程度，體驗成就感和失落感，為日后成為參賽骨干打下基礎。

cs=7.915×10-7T3-6.439×10-4T2+

1.765×10-1T-16.17

(9)

撲熱息痛結晶過程的生長速率G、初級成核速率b1和二次成核速率b2公式分別如式(6)、(4)和(5)所示.式中動力學參數由Li等[19]的實驗擬合得到,即生長速率常數kg為45.5 m/s,指數γ為1.24,晶體活化能Ea為41.3 kJ/mol,初級成核常數kb1為0.192 1 s-1/kg,晶體和溶液的界面能σ為4.25 mJ/m2,玻爾茲曼常數k為1.380 6×10-23J/K,二次成核常數kb2、α和β分別為1×105s-1/kg、2.08和0.713.

3 模擬與計算

3.1 三維建模與網格劃分

通過Solidworks軟件建立攪拌結晶器的幾何圖形,如圖3所示.結晶器元件的相關尺寸見表2.ANSYS meshing生成的網格如圖4所示.扭曲度(skewness)用于評價網格質量,即趨于理想網格的程度,取值為0～1,數值越大網格質量越差.基于歸一化的正角度,扭曲度被定義為

表2 結晶器的尺寸

圖3 簡化的三維攪拌結晶器

(a)三維視圖

(10)

式中:θmax為網格中最大角度,θmin為網格中最小角度,θe為理想網格角度.本文劃分的網格最大扭曲度為0.80,平均扭曲度為0.23,網格質量良好.

3.2 數值模擬設置

采用ANSYS FLUENT進行CFD-PBM模擬仿真.模型設置:選擇歐拉(Eulerian)多相流,湍流模型選擇RNGk-ε模型,該模型引入了旋轉、曲率相關項,能夠模擬旋流等中等復雜的流動[22].對于攪拌槳旋轉的模擬,采用滑移網格方法(sliding mesh method,SMM).降溫曲線通過用戶定義函數(UDF)在壁面上進行邊界條件設置.PBM采用離散法進行求解,離散成30個區域,成核和生長動力學模型通過用戶定義函數(UDF)實現.速度與壓力的耦合算法為Coupled,殘差設置為1×10-4,在單節點64核內存256 GB超算平臺進行并行計算.

4 模擬結果與分析

4.1 模型驗證

(1)網格獨立性驗證

通過網格獨立性驗證排除網格數量對模擬結果的影響.一共劃分了網格數量為120 877、155 210、173 774和231 714的4套網格,對比不同網格數量的速度分布,如圖5所示.最終173 774個網格被確定為精度和數值成本權衡之間的最佳選擇.

(2)攪拌結晶器速度流場流型驗證

圖6(a)為ANSYS FLUENT軟件的模擬結果,圖6(b)為文獻[23]中攪拌結晶器的計算結果,通過對比發現,兩者的流場分布大致相符.結晶器的流場形成上下兩個循環,使得物料在結晶器內得到充分的接觸與混合.

(3)模擬粒度分布與實驗值對比

在與Li等[19]相同結晶條件下,模擬結果與實驗結果的對比如圖7所示.圖7(a)、(b)分別為體積密度函數V(L,y,t)和數量密度函數n(L,y,t)下的粒度分布.體積密度(volume density)函數是指單位物理空間(體積)、單位顆粒長度的顆粒體積.數量密度(number density)函數是指單位物理空間(體積)、單位顆粒長度的顆粒數.由圖可知,粒度分布的模擬結果與實驗結果偏差較小,且模擬得到的粒度分布與實驗趨勢基本一致,因此模擬結果與實驗結果有較好的一致性.模擬值和實驗值的偏差可能是結晶動力學模型的簡化,以及CFD-PBM耦合模型所采用的數值方法帶來的不確定性造成的.由于體積密度函數下的粒度分布更能體現顆粒分布的特性,后續研究中粒度分布都將采用體積密度函數下的粒度分布.

(a)體積密度函數下的粒度分布

4.2 攪拌結晶器的流場分析

圖8為400 r/min攪拌速度下的顆粒體積分數φ分布.攪拌槳作為上下部分的分界線,各自形成循環流股.因為攪拌槳位置靠下,下部循環流股的流速比上部要快,所以下方顆粒被帶起,使得攪拌槳下方的攪拌相對更均勻.

圖8 400 r/min攪拌速度下的顆粒體積分數分布

圖9為結晶器中不同X坐標下顆粒體積分數的縱向分布曲線圖,其中X為距離攪拌槳中心的水平距離.由圖可知,顆粒體積分數隨著高度增加呈階梯式分布.除了距離攪拌槳中心較近的區域外,其他區域均存在兩個階梯.從結晶器底部離開,顆粒體積分數會出現驟降,然后進入一個平臺,即顆粒分布均勻區域,這是由于距離攪拌槳較近,湍流強度較大,混合效果較好.

圖9 不同X坐標下顆粒體積分數縱向分布

4.3 晶種的影響

圖10為加入晶種和未加入晶種的模擬結果,加入的晶種質量為0.03 g,晶種粒徑為75～106 μm.由圖可知,加入晶種的粒度分布峰值比未加入晶種的高,即產生的顆粒數更多.同時加入晶種獲得的顆粒平均粒徑256.20 μm大于未加入晶種的平均粒徑254.95 μm.由于加入晶種數量較少,因此兩者偏差較小.為了獲得顆粒數多、平均粒徑大的粒度分布,可以適當加入晶種.

圖10 加入晶種和未加入晶種的粒度分布對比

圖11為加入晶種和未加入晶種條件下,生長和成核過程產生的源項Sφ隨時間的變化曲線,即縱坐標體現了晶體生長速率和成核速率的大小.由圖可知,在剛開始結晶過程中,加入晶種的源項高于未加入晶種的源項,這是由于自發成核的過飽和度高于顆粒生長所需的過飽和度.在5 000 s之前,源項曲線呈現快速增長趨勢,是顆粒生長和成核過程的關鍵階段,該階段中加入晶種的源項始終大于未加入晶種的源項,即加入晶種的生長速率和成核速率更大,產生的顆粒數更多、平均粒徑更大.圖12為在3 000 s時顆粒體積分數分布圖,其中圖12(c)為加入晶種和未加入晶種的顆粒體積分數差值的分布.圖中大部分區域的數值為正,說明結晶器內加入晶種的顆粒體積分數大于未加入晶種的.所以在結晶器內,加入晶種條件下,顆粒數更多,顆粒的體積分數更大.

圖11 加入晶種和未加入晶種條件下結晶動力學源項隨時間的變化曲線

(a)加入晶種

因此,在結晶過程中,加入晶種可以避免新晶體成核造成的高過飽和度,從而獲得具有顆粒數多、平均粒徑大的粒度分布晶體.

4.4 破碎過程的影響

Luo破碎模型是常用的破碎模型,適用于湍流中顆粒破碎的預測.圖13為在考慮破碎和不考慮破碎條件下顆粒的粒度分布以及與實驗結果的對比.由圖可知,考慮破碎的粒度分布峰值更大,顆粒數更多,更接近于實驗結果.因此,補充Luo破碎模型可以減少由于結晶動力學模型簡化帶來的偏差.

圖14為結晶動力學源項Sφ隨時間的變化曲線,其中兩者差值在一定程度上體現了Luo破碎模型破碎速率的大小.由圖可知,考慮破碎的結晶過程粒度分布峰值更大,主要原因是前期考慮破碎的結晶動力學源項明顯高于不考慮破碎的,即由于補充Luo破碎模型,顆粒屬性受到了破碎速率的影響,顆粒數更多.結晶后期,破碎速率的影響逐漸減小,兩者源項相差較小.因此采用合適的破碎模型,考慮更為全面的結晶動力學模型,能夠獲得提高最終顆粒粒度分布預測的能力.

圖14 不同破碎條件下結晶動力學源項隨時間的變化曲線

4.5 攪拌速度的影響

圖15為不同攪拌速度下顆粒體積分數的分布.由圖可知,隨著攪拌速度的增大,傳遞給流體的機械能越來越大,結晶器內的固相體積分數分布趨于均勻.在攪拌速度為300、400 r/min時,結晶器中產品顆粒的分布在上下兩部分有明顯差異,攪拌槳作為上下部分的分界線,各自形成循環流股.在攪拌速度為500 r/min時,結晶器內的顆粒體積分數分布基本均勻.

(a)200 r/min

圖16為不同攪拌速度下粒度分布的對比.由圖可知,隨著攪拌速度增大,粒度分布峰值增加,顆粒數增加,其中300、400 r/min的峰值相差不大.由圖17可知,隨著攪拌速度Nr增大,顆粒的最終平均粒徑減小.主要是因為隨著攪拌速度的增大,結晶器中的流體逐漸混合均勻,有利于結晶器內的傳熱過程,進而使得溶液的平均溫度下降,過飽和度略有增加,導致成核速率增大;同時攪拌速度較大時,湍流耗散率較高,進一步加強了二次成核,產生更多細小晶體,所以顆粒數增加,平均晶體尺寸減小.因此攪拌程度對于粒度分布、平均粒徑等性質有較大的影響,為了得到理想的顆粒產品,必須對結晶工藝的攪拌速度進行合理的設計和控制.

圖16 不同攪拌速度下粒度分布的對比

圖17 平均粒徑隨攪拌速度的變化曲線

4.6 降溫速率的影響

線性降溫曲線公式為

(11)

采用不同降溫速率k對結晶過程進行模擬,如圖18所示.當降溫速率k增大時,粒度分布峰值增加,即顆粒數增加,并且粒度分布變窄,其中降溫速率k為-1/150 ℃/s和-1/200 ℃/s的峰值相近.由圖19可知,隨著降溫速率k的增大,顆粒的平均粒徑逐漸減小,雖然降溫速率k從-1/150 ℃/s到-1/200 ℃/s的平均粒徑略有增大,但是相差很小.因此適當采用較為平緩的線性降溫曲線,能夠獲得更窄的粒度分布以及更多的顆粒數,但是產品最終的平均粒徑會更小.所以在結晶工藝中控制降溫速率是很有必要的.

圖18 不同降溫速率k下的粒度分布

圖19 平均粒徑隨降溫速率的變化曲線

5 結 論

(1)攪拌結晶器內的流場形成上下兩個循環,顆粒體積分數隨著高度增加呈階梯式分布.

(2)在結晶過程中,適當加入晶種可以避免新晶體成核造成的高過飽和度,從而獲得顆粒數多、平均粒徑大的晶體.

(3)采用合適的破碎模型,考慮更為全面的結晶動力學模型,能夠提高最終晶體的粒度分布預測能力.

(4)隨著攪拌速度增大,產品顆粒數增加,顆粒的平均粒徑減小.

(5)線性降溫速率減小,產品顆粒數減少,顆粒的平均粒徑增大.

整體而言,本文CFD-PBM耦合模型的模擬結果對結晶工藝的設計具有一定的指導意義.雖然Luo破碎模型滿足適用條件,但是和真實模型相比仍存在一定的誤差,后續將結合實驗對破碎模型做進一步探究.