間歇冷卻結(jié)晶過程的CFD-PBM數(shù)值模擬

覃 渫 蘭, 張 磊, 都 健

(大連理工大學(xué) 化工學(xué)院 化工系統(tǒng)工程研究所, 遼寧 大連 116024)

0 引 言

溶液結(jié)晶是化學(xué)、制藥和食品工業(yè)中應(yīng)用較為廣泛的分離和純化技術(shù)之一[1].在結(jié)晶過程中,由于對物理化學(xué)和機(jī)械特性的強(qiáng)烈依賴性,需要控制結(jié)晶產(chǎn)品質(zhì)量[2].粒度分布(crystal size distribution,CSD)是評價最終晶體質(zhì)量的重要指標(biāo)之一[3-4].

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

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

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

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

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

1 數(shù)學(xué)模型

1.1 粒數(shù)衡算模型

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

?[Γ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的破碎顆粒概率密度函數(shù);P為一個顆粒破碎產(chǎn)生子顆粒的數(shù)量.

1.2 CFD-PBM耦合模型

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

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

(3)

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

圖1 CFD-PBM耦合模型

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

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

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

2 結(jié)晶動力學(xué)模型

2.1 結(jié)晶過程

(1)晶體成核過程

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

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

(4)

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

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

(5)

式中:kb2為二次成核常數(shù);ms為晶體質(zhì)量;α為常數(shù)2.08;β為常數(shù)0.713.

(2)晶體生長過程

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

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

(6)

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

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

Luo破碎模型的破碎速率為

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

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

(7)

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

2.2 撲熱息痛冷卻結(jié)晶的動力學(xué)模型

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

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

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

表1 物性參數(shù)表

(8)

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

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

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

1.765×10-1T-16.17

(9)

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

3 模擬與計算

3.1 三維建模與網(wǎng)格劃分

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

表2 結(jié)晶器的尺寸

圖3 簡化的三維攪拌結(jié)晶器

(a)三維視圖

(10)

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

3.2 數(shù)值模擬設(shè)置

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

4 模擬結(jié)果與分析

4.1 模型驗證

(1)網(wǎng)格獨立性驗證

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

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

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

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

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

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

4.2 攪拌結(jié)晶器的流場分析

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

圖8 400 r/min攪拌速度下的顆粒體積分?jǐn)?shù)分布

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

圖9 不同X坐標(biāo)下顆粒體積分?jǐn)?shù)縱向分布

4.3 晶種的影響

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

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

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

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

(a)加入晶種

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

4.4 破碎過程的影響

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

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

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

4.5 攪拌速度的影響

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

(a)200 r/min

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

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

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

4.6 降溫速率的影響

線性降溫曲線公式為

(11)

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

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

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

5 結(jié) 論

(1)攪拌結(jié)晶器內(nèi)的流場形成上下兩個循環(huán),顆粒體積分?jǐn)?shù)隨著高度增加呈階梯式分布.

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

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

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

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

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