陳 猛，劉馬林，唐亞平，劉榮正，劉 兵，邵友林，常家興

（清華大學(xué)核能與新能源技術(shù)研究院；先進(jìn)核能技術(shù)協(xié)同創(chuàng)新中心，北京100084）

顆粒包覆過程，即在顆粒的表面均勻包覆上一層功能材料，具有保護(hù)、緩釋、提高環(huán)境相容性等用途。在工業(yè)生產(chǎn)中，尤其是制藥行業(yè)應(yīng)用廣泛，同時(shí)也是高溫氣冷堆核燃料元件制備過程的重要環(huán)節(jié)。顆粒包覆是顆粒性能改進(jìn)的主要方法之一，而包覆層的均勻性是顆粒包覆過程追求的目標(biāo)。高溫氣冷堆核燃料元件采用的是4層包覆燃料顆粒，因?yàn)楹巳剂系墓逃邪踩卣髟O(shè)計(jì)，要求包覆層滿足一定的技術(shù)指標(biāo)，其中顆粒包覆層的均勻性則是決定高溫氣冷堆安全性的關(guān)鍵因素之一。

實(shí)驗(yàn)研究顆粒均勻包覆工藝是很重要的一個(gè)方面，如果通過數(shù)學(xué)建模和數(shù)值模擬，能分析出操作參數(shù)對(duì)顆粒均勻包覆的影響規(guī)律，可指導(dǎo)實(shí)驗(yàn)研究方向，減少實(shí)驗(yàn)次數(shù)。對(duì)顆粒包覆過程進(jìn)行準(zhǔn)確的數(shù)值模擬是研究均勻包覆機(jī)理的一種有效手段，進(jìn)而可優(yōu)化包覆過程，進(jìn)行包覆工藝的放大，也是提高顆粒包覆過程經(jīng)濟(jì)性的有效途徑。如何進(jìn)行包覆過程的準(zhǔn)確數(shù)值模擬，以往的研究給出了很多顆粒包覆過程數(shù)值模擬方法，但缺少比較系統(tǒng)和深入的對(duì)比研究[1]。

本文中首先分析了顆粒均勻包覆過程的主要特征，指出顆粒包覆過程的研究層次，即單顆粒表面均勻和全顆粒均勻2個(gè)層次，指出目前大量的模型僅局限在顆粒之間包覆層質(zhì)量均勻的模擬研究。然后闡述不同數(shù)學(xué)模型的具體實(shí)施過程，指出其優(yōu)勢(shì)和劣勢(shì)，重點(diǎn)對(duì)CFD-DEM方法包覆過程模擬進(jìn)行深入分析，并基于運(yùn)動(dòng)－吸附－沉積機(jī)理發(fā)展新型包覆過程模擬方法。最后基于現(xiàn)有文獻(xiàn)中的顆粒包覆過程數(shù)值模擬方法對(duì)比分析研究，指出包覆模擬研究的規(guī)律，給出了顆粒包覆過程準(zhǔn)確數(shù)值模擬的發(fā)展趨勢(shì)，即基于單顆粒尺度的表面均勻包覆機(jī)理研究和基于反應(yīng)器尺度的物理場(chǎng)分布研究相耦合的多尺度研究方法。

1 顆粒包覆技術(shù)

目前工業(yè)上成功采用的主要有2種顆粒包覆過程，即流化床過程和轉(zhuǎn)鼓過程，2種過程在核燃料制備過程中均有采用，如圖1所示。高溫氣冷核反應(yīng)堆采用的燃料元件由包覆燃料顆粒和石墨基體組成。其中陶瓷型二氧化鈾（UO2）核燃料顆粒制成之后，需要在其表面均勻包覆上疏松熱解炭層、內(nèi)致密熱解炭層、碳化硅層、外致密熱解炭層（即3層結(jié)構(gòu)各向同性（TRISO）顆粒制備，采用流化床過程）還需要顆粒表面均勻包覆上一層石墨層（即穿衣顆粒制備，采用轉(zhuǎn)鼓過程），才能將之與石墨基體混合壓制成燃料元件。TRISO顆粒包覆層是核燃料顆粒安全性第1道保障，而穿衣層則決定了顆粒在壓制過程的完整性，因此，顆粒包覆過程，即在顆粒的表面均勻包覆上一層功能材料，是高溫氣冷堆燃料元件制備過程的重要環(huán)節(jié)，必須予以高度重視。顆粒包覆裝置如圖1所示，盡管流化床和轉(zhuǎn)鼓2個(gè)過程差別較大，但是將其用于顆粒包覆過程具有內(nèi)在的相似性，很多數(shù)值模擬方法可以在2個(gè)過程通用[2]。

圖1 顆粒包覆裝置示意圖Fig.1 Schematic diagram of particle coating device

2 顆粒包覆層均勻性分析

顆粒包覆層是否均勻是顆粒包覆過程追求的目標(biāo)，包括單顆粒均勻（顆粒內(nèi)均勻）和全顆粒均勻（顆粒間均勻）2個(gè)方面。單顆粒均勻包覆指的是該顆粒上的包覆層性質(zhì)均勻，包括厚度、密度、各向異性度等；全顆粒均勻包覆包括2個(gè)方面，一是每個(gè)顆粒上的包覆層質(zhì)量均勻，二是每個(gè)顆粒上的包覆層性質(zhì)均勻，如圖2所示。

圖2 顆粒包覆層均勻性分析Fig.2 Analysis of particle coating uniformity

目前關(guān)于包覆過程的數(shù)值模擬多是針對(duì)全顆粒包覆層的質(zhì)量均勻而言的，關(guān)于單顆粒性質(zhì)均勻方面的數(shù)值模擬非常困難，研究較少，僅有針對(duì)單顆粒包覆層厚度均勻提出的一種耦合球諧函數(shù)的模擬方法。本文中針對(duì)全顆粒質(zhì)量均勻包覆和單顆粒包覆層厚度均勻包覆提出的模型做詳細(xì)的對(duì)比分析。

3 數(shù)值模擬方法

在顆粒間包覆層質(zhì)量均勻方面，數(shù)值模擬方法包括表面更新模型（surface renewal model）、蒙特卡羅模擬方法（Monte－Carlo method）、分區(qū)群體平衡模型（compartment－based population balance model）、DEM-C-PBM模型（discrete element methodcompartment-population balance model）、CFD-DEM模型（computational fluid dynamics-discrete element method）等。在單顆粒性質(zhì)均勻方面研究較少，僅有耦合球諧函數(shù)的CFD-DEM模型（CFD-DEM&spherical harmonic formulation）等少數(shù)方法。

3.1 全顆粒包覆均勻數(shù)值模擬

3.1.1 表面更新模型

顆粒包覆的一般過程是顆粒在不斷運(yùn)動(dòng)的過程中循環(huán)往復(fù)通過含有待沉積產(chǎn)物的區(qū)域，待沉積產(chǎn)物附著在顆粒表面形成包覆層[3]，基本原理見圖3。

表面更新模型是一個(gè)基礎(chǔ)的數(shù)學(xué)分析模型，起源于連續(xù)循環(huán)體系的數(shù)學(xué)分析[4]，很早被用于分析顆粒包覆工藝參數(shù)的影響，其基本原理認(rèn)為反應(yīng)器內(nèi)分為包覆區(qū)和非包覆區(qū)，顆粒每通過一次包覆區(qū)質(zhì)量增加為xi，則通過n次后總的包覆質(zhì)量為：

圖3 顆粒包覆過程原理示意圖Fig.3 Schematic diagram of particle coating process

顆粒包覆質(zhì)量均勻性，即變化系數(shù)cv，m（coefficient of variation）可以用顆粒包覆質(zhì)量的標(biāo)準(zhǔn)偏差來表示，

式中：NP是顆粒數(shù)量；mc，i是顆粒i的包覆層質(zhì)量；ma_c是所有顆粒的包覆層平均質(zhì)量。

經(jīng)過數(shù)學(xué)推導(dǎo)，可得到顆粒包覆層質(zhì)量的變化系數(shù)和單顆粒循環(huán)時(shí)間以及每次包覆量的關(guān)系

式中：μcpp和σcpp是在一次循環(huán)時(shí)間內(nèi)全體顆粒每次通過包覆區(qū)的包覆層增加值的平均值和標(biāo)準(zhǔn)偏差；μct和σct分別是所有顆粒循環(huán)時(shí)間的平均值和標(biāo)準(zhǔn)偏差；Tcoat是顆?？偟陌矔r(shí)間。

顆粒包覆層質(zhì)量分布如圖4所示。

圖4 顆粒包覆層質(zhì)量分布Fig.4 Mass distribution of particle coating layer

圖5 顆粒循環(huán)時(shí)間分布及每次通過包覆區(qū)質(zhì)量增量分布Fig.5 Particle circulation-time distribution and coating-per-pass distribution

μct和σct可以通過加入磁性顆粒以及檢測(cè)線圈、或者其他示蹤方法，例如熒光增強(qiáng)顆粒示蹤技術(shù)（PEPT）來獲得，其測(cè)量相對(duì)比較容易，研究結(jié)論較多一些。μcpp和σcpp則可以通過加入一定顏色的包覆染料來測(cè)得，研究較少，可以給出包覆層不均勻性的主要起源，用于和實(shí)驗(yàn)結(jié)果進(jìn)行分析。典型的研究例子是文獻(xiàn)[5]中利用磁性顆粒示蹤，并用反應(yīng)器外部側(cè)壁不同高度纏繞金屬線圈來檢測(cè)磁性顆粒的運(yùn)動(dòng)軌跡，獲取顆粒循環(huán)時(shí)間分布，得到不同氣速下的顆粒循環(huán)時(shí)間分布，如圖5a所示。由圖5a可以看出，氣速越大，包覆層質(zhì)量分布越集中。利用藍(lán)色染料噴射包覆，用熒光光度計(jì)測(cè)量單次循環(huán)中（氣速300 m3/h），每個(gè)顆粒通過包覆區(qū)的包覆層質(zhì)量增量，得到其分布，如圖5b所示，可以看出包覆層質(zhì)量增量分布非常寬。

表面更新模型可以用于深入分析包覆過程影響因素，即2種決定因素對(duì)顆粒包覆均勻性的影響，認(rèn)為包覆區(qū)內(nèi)的濃度分布、顆粒之間的相互阻擋、顆粒運(yùn)動(dòng)的脈動(dòng)不穩(wěn)定性是決定包覆層增加量的不均勻分布的主要因素，并可以通過該模型的研究給出顆粒均勻包覆的改進(jìn)方向，即取消高濃度包覆區(qū)、減少噴嘴附近的變化以及降低顆粒運(yùn)動(dòng)的脈動(dòng)不穩(wěn)定性等。文獻(xiàn)[6]對(duì)此方法進(jìn)行了發(fā)展，通過以軸向Peclet準(zhǔn)數(shù)為特征的平流運(yùn)動(dòng)和擴(kuò)散運(yùn)動(dòng)的組合來精確地描述顆粒的軸向運(yùn)動(dòng)，開發(fā)了基于更新理論的數(shù)學(xué)框架，并且獲得顆粒間包覆層變化性的表達(dá)式，其解釋了顆粒在涂布機(jī)內(nèi)部的停留時(shí)間的變化。

圖6 利用蒙特卡羅方法計(jì)算得到的顆粒循環(huán)時(shí)間分布圖Fig.6 Particle circulation-time distribution calculated by Monte Carlo method

圖7 包覆時(shí)間對(duì)顆粒包覆層質(zhì)量變化的影響Fig.7 Effect of coating time on mass coating variability

3.1.2 蒙特卡羅模擬方法

蒙特卡羅方法是基于重復(fù)的隨機(jī)性采樣的一種隨機(jī)分析方法。此方法本身是普適性的，可以應(yīng)用在不同領(lǐng)域，在包覆過程中的應(yīng)用主要取決于顆粒隨機(jī)運(yùn)動(dòng)行為和噴射通量的空間分布行為，文獻(xiàn)[7-8]中曾用該方法研究顆粒包覆過程。蒙特卡羅方法模擬顆粒包覆過程的主要原理是：顆粒在包覆區(qū)和非包覆區(qū)之間循環(huán)往復(fù)運(yùn)動(dòng)，當(dāng)顆粒通過包覆區(qū)時(shí)，每個(gè)顆粒表面的顆粒包覆層增量取決于顆粒暴露與噴口的截面積，局部噴射通量以及顆粒與噴射源之間的距離。當(dāng)顆粒的初始位置選定后，顆粒在下一個(gè)時(shí)刻的位置取決于一個(gè)隨機(jī)量決定的值。例如在文獻(xiàn)[12]中對(duì)顆粒在導(dǎo)流板流化床中包覆均勻性進(jìn)行了蒙特卡羅模擬，顆粒速度計(jì)算公式為

基于該顆粒循環(huán)時(shí)間分布圖，通過實(shí)驗(yàn)測(cè)得的噴霧通量分布、噴霧區(qū)形狀及面積大小、顆粒投影面積可以計(jì)算得到顆粒包覆質(zhì)量均勻性變化系數(shù)與包覆時(shí)間之間的關(guān)系，如圖7所示。由圖可以看出，顆粒包覆質(zhì)量均勻性變化系數(shù)cv，m與包覆時(shí)間函數(shù)1/t0.5成正比關(guān)系。

文獻(xiàn)[9]中發(fā)展了一種蒙特卡洛方法，研究了實(shí)心圓錐和平行圓柱2種包覆區(qū)形狀對(duì)模擬結(jié)果的影響，圖8為2種模型中顆粒投影表面積離散化的示意圖。結(jié)果發(fā)現(xiàn)實(shí)心圓錐模型的模擬效果沒有平行圓柱模型好，并對(duì)噴射速度和導(dǎo)流板低端縫隙大小進(jìn)行了影響因素分析，發(fā)現(xiàn)對(duì)于具有較低氣體速度和間隙高度的情況，包覆層分布更寬。

圖8 不同模型中顆粒投影表面積離散化示意圖Fig.8 Schematic diagrams of discretization of particle projection area in different models

3.1.3 分區(qū)群體平衡模型

表面更新模型和蒙特卡羅方法模擬實(shí)質(zhì)上都是二元分區(qū)模型，研究者將分區(qū)模型和群體平衡模型耦合，通過分析包覆過程的不同步驟，提出了多種類型的兩室或多室群體平衡模型，用于預(yù)測(cè)多分散顆粒的包覆層質(zhì)量分布。文獻(xiàn)[10]中給出了一種較為簡(jiǎn)潔的處理方式，將整個(gè)反應(yīng)器分成3種區(qū)域——有效包覆區(qū)G、有效干燥區(qū)F、非有效區(qū)E，如圖9所示。

基于顆粒在每個(gè)區(qū)內(nèi)混合均勻，顆粒不會(huì)磨損、破碎，顆粒不會(huì)跨過F區(qū)直接在G區(qū)和E區(qū)之間交換等假設(shè)，顆粒在G區(qū)和F區(qū)之間傳遞速率為a，在F區(qū)和E區(qū)之間傳遞速率為b，根據(jù)每個(gè)區(qū)的物質(zhì)守恒可以給出區(qū)域間的顆粒群體平衡模型：

圖9 分區(qū)群體平衡模型示意圖Fig.9 Schematic diagram of com partment-based population balance model

式中：f為F區(qū)群體密度函數(shù)；g為G區(qū)群體密度函數(shù)；e為E區(qū)群體密度函數(shù)；α為噴射的活性區(qū)占的分?jǐn)?shù)；β為包覆過程整個(gè)床占的分?jǐn)?shù)；r為包覆層質(zhì)量增加速率；t代表時(shí)間；w代表包覆層質(zhì)量。

通過上述模型可以研究包覆區(qū)的大小、包覆時(shí)間、單位時(shí)間的包覆增加量對(duì)包覆均勻性的影響規(guī)律。研究認(rèn)為包覆質(zhì)量是包覆時(shí)間的函數(shù)，考慮了顆粒流動(dòng)的非均勻性帶來的包覆變異系數(shù)變化。分區(qū)模型關(guān)鍵參數(shù)包括噴射區(qū)的相對(duì)尺寸、分區(qū)之間的流動(dòng)速率、噴射速率以及噴射時(shí)間等。文獻(xiàn)[10]中得到結(jié)果如圖10所示。

圖10 以a/b和非活性區(qū)域尺寸為變量的包覆層質(zhì)量分布的相對(duì)方差圖Fig.10 Variance of the coating distribution as a function of a/b and size of non-active domain

由圖可以看出：如果包覆時(shí)間足夠長(zhǎng)，包覆層分布接近于正態(tài)分布；非包覆區(qū)的增加會(huì)增大包覆層的不均勻性；a和b相等時(shí)不均勻性最大，但是無論a和b哪個(gè)增大都會(huì)使得包覆分布更加均勻，但是a的增大，即有效干燥區(qū)和有效包覆區(qū)之間的傳遞速率增大，使得包覆分布均勻的效果更加明顯；非活性區(qū)域尺寸的增大通常會(huì)導(dǎo)致涂層分布的變化的增加；若非活性區(qū)域內(nèi)存在停滯，則會(huì)產(chǎn)生具有較大變化的涂層分布；區(qū)域之間轉(zhuǎn)移速率的增加會(huì)使涂層分布變窄；當(dāng)所有區(qū)域之間的顆粒轉(zhuǎn)移速率相等時(shí)，產(chǎn)生最大的涂層變化；與活性干燥區(qū)域和非活性區(qū)域相比，活性干燥區(qū)域和活性噴霧區(qū)域之間的交換速率的增大會(huì)產(chǎn)生更窄的涂層分布。

文獻(xiàn)[11]中還給出一種耦合DEM-C-PBM的分區(qū)模型，用于轉(zhuǎn)鼓中顆粒包覆過程的數(shù)值模擬。如圖11所示為在沒有擋板的水平旋轉(zhuǎn)鼓室中的顆粒流動(dòng)。細(xì)黑線表示在短時(shí)間內(nèi)的顆粒軌跡，箭頭表示包覆機(jī)中的一般顆粒流動(dòng)，噴霧區(qū)由靠近床自由表面的實(shí)心黑線標(biāo)記，紅色陰影區(qū)域表示被動(dòng)床區(qū)域。在該模型中，靠近自由表面并且靠近鼓壁的顆粒以環(huán)形循環(huán)，其中在自由表面處有一小部分暴露于噴霧中；在顆粒床的中心附近還存在顆粒速度較小的區(qū)域；在外部循環(huán)區(qū)域和中心區(qū)域之間發(fā)生顆粒的連續(xù)交換。

圖11 旋轉(zhuǎn)鼓室中顆粒流動(dòng)示意圖Fig.11 Schematic diagram of flow of particles in a rotating drum

基于這種流動(dòng)模式，在循環(huán)回路中的顆粒流動(dòng)，包括噴霧區(qū)和一部分活性床區(qū)，如圖12所示。參數(shù)α和β分別表示噴霧區(qū)域和床區(qū)域尺寸占顆?？倲?shù)的百分比，數(shù)量γ表示循環(huán)回路中床區(qū)占的百分?jǐn)?shù)，循環(huán)回路中的顆粒的數(shù)量流速由N·loop循環(huán)給出，并且在活性區(qū)和被動(dòng)床區(qū)的各個(gè)隔室之間交換顆粒的數(shù)量流量為N·exch/NB。這種耦合PBM的分區(qū)模型中的參數(shù)可由DEM數(shù)值模擬直接獲得，形成基于DEM的參數(shù)評(píng)估分區(qū)PBM模型，可用于包覆過程模擬，是一種多尺度模型。

圖12 旋轉(zhuǎn)鼓室流動(dòng)分區(qū)模型示意圖Fig.12 Schematic diagram of flow of compartmentpopulation particles in a rotating drum

3.1.4 CFD-DEM耦合模型

DEM是一種描述顆粒移動(dòng)和碰撞行為的模擬方法。其可與CFD耦合，模擬顆粒－氣體相互作用引起的各種現(xiàn)象，例如顆粒流態(tài)化行為等，稱為CFD-DEM耦合模型[12-13]。DEM可以深入到單顆粒層次，所以可以準(zhǔn)確地描述每個(gè)顆粒的變化和運(yùn)動(dòng)軌跡，尤其適應(yīng)于包覆過程的數(shù)值模擬，這種算法的核心在于顆粒受力的計(jì)算。

式中：mi為顆粒i的質(zhì)量；vi為顆粒i的速度；ri為顆粒i的位置矢量；t為時(shí)間；p為壓力；Vi為顆粒i的體積；β為相間動(dòng)量傳遞系數(shù)；ε是孔隙率；ug為氣體速度；g為重力加速度；Fc，i為顆粒i受到的接觸力；FA，i為顆粒i受到的顆粒之間的黏合力。

方程（9）右側(cè)的力分別是由于壓力梯度、阻力、重力、接觸力和黏合顆粒之間的相互作用力。

關(guān)于CFD-DEM耦合模型用于顆粒運(yùn)動(dòng)體系數(shù)值模擬的工作已經(jīng)有很多進(jìn)展，研究人員將CFDDEM耦合模型和包覆區(qū)的分布結(jié)合起來，研究包覆參數(shù)的影響。例如文獻(xiàn)[14]中對(duì)比分析了2種不同的包覆反應(yīng)器構(gòu)造對(duì)顆粒包覆均勻度的影響，如圖13所示。結(jié)果表明，Wurster制粒機(jī)的特征在于具有比較窄的停留時(shí)間分布，而使用頂部噴霧制粒機(jī)則導(dǎo)致較寬的停留時(shí)間分布，這說明Wurster制粒機(jī)可以獲得更加均勻的顆粒包覆分布。該模擬的核心還是將反應(yīng)器區(qū)域先分成幾個(gè)不同功能的區(qū)域，例如有效包覆區(qū)和非有效包覆區(qū)。顆粒在包覆區(qū)內(nèi)的停留時(shí)間分布決定了顆粒包覆層的厚度，因而決定了顆粒包覆層質(zhì)量的分布。

3.2 單顆粒包覆均勻數(shù)值模擬

以上均是對(duì)全體顆粒包覆質(zhì)量均勻進(jìn)行數(shù)值模擬和描述的方法，此類研究較多，而對(duì)單顆粒表面包覆均勻性的數(shù)值模擬研究較少。僅有的文獻(xiàn)[15]中均采用DEM方法進(jìn)行顆粒運(yùn)動(dòng)的數(shù)值模擬，然后耦合蒙特卡羅方法模擬計(jì)算了單顆粒表面沉積均勻性。具體方法是首先將顆粒表面離散化，如圖14所示，顆粒表面包覆層質(zhì)量m為表面積分，

式中：h為包覆層增加厚度；A是包覆層面積；ρ是包覆層密度。

假定顆粒在均勻的噴射場(chǎng)內(nèi)，而顆粒的轉(zhuǎn)動(dòng)是隨機(jī)的，因此與均勻流場(chǎng)碰撞的面是包覆層增加的面，即

圖13 不同包覆機(jī)中顆粒在噴霧區(qū)的停留時(shí)間分布Fig.13 Residence time distribution of particles in spray zone of different coating devices

圖14 耦合球諧函數(shù)的CFD-DEM模型示意圖Fig.14 Schematic diagram of CFD-DEM&spherical harmonic formulation model

式中：n是外延法向量，如圖14b所示；urel是顆粒與噴射場(chǎng)的相對(duì)運(yùn)動(dòng)速度。

計(jì)算出每一個(gè)離散面的包覆層增加厚度，即可以得到單顆粒表面包覆不均勻性的模擬計(jì)算。

文獻(xiàn)[16]中另外一種模擬計(jì)算單顆粒表面沉積均勻性的模型是耦合球諧函數(shù)的CFD-DEM模型。在該模型中，也是首先將顆粒表面離散化，用顆粒表面分布函數(shù)來表示單顆粒表面沉積的均勻性，如圖14所示。圖14a表示將單個(gè)球形顆粒暴露于均勻噴霧，顆粒在涂層試驗(yàn)之間隨機(jī)旋轉(zhuǎn)，顆粒陰影表示顆粒表面的那部分上的膜厚度；圖14b為單顆粒包覆均勻性模擬時(shí)的顆粒表面離散化示意圖，面積為ΔAi的面板位于半徑為R的球形顆粒上，面板的單位法向量為ni，相對(duì)于粒子的固定球體參照系測(cè)量角度分別為φp和θi；圖14c表示噴霧顆粒在 φp、θp方向上與球表面上的點(diǎn)撞擊。該模型采用的是液滴沉積原理，即假定沉積產(chǎn)物在包覆反應(yīng)器內(nèi)以液滴形式存在，當(dāng)液滴與顆粒從一個(gè)角度碰撞時(shí)，液滴的體積會(huì)映射到顆粒表面上，形成包覆層。該映射方式的計(jì)算，即球體表面分布方式的計(jì)算是采用球面諧波函數(shù)的方式進(jìn)行的。例如，如果考慮2個(gè)液滴（體積為V1和V2）沉積到顆粒表面的過程，則映射到顆粒表面的質(zhì)量分別為V1g（θ）和V2g（θ），顆粒表面沉積產(chǎn)物的增加過程可以表示為

式中：f為顆粒表面沉積產(chǎn)物分布函數(shù)；g為正態(tài)分布函數(shù)。式子左邊可以用球諧函數(shù)來表示，即得到顆粒表面沉積產(chǎn)物的分布。

將氣體、顆粒、液滴用CFD-DEM進(jìn)行模擬，將顆粒和液滴的碰撞、包覆過程用上式來描述，可以得到整個(gè)體系的顆粒包覆過程，最終模擬結(jié)果如圖15所示，可以看出顆粒沉積變化系數(shù)與包覆時(shí)間、顆粒速度之間的關(guān)系。結(jié)果表明，當(dāng)顆粒速度越大、包覆時(shí)間越長(zhǎng)，則顆粒表面沉積變化系數(shù)越小，即表面包覆越均勻；完全隨機(jī)覆蓋的的理想斜率如圖15中虛線所示。

圖15 CFD-DEM模擬結(jié)果示意圖Fig.15 Schematic diagram of CFD-DEM simulation results

4 模型比較分析

從上述模型比較分析來看，表面更新模型、蒙特卡羅模型、分區(qū)群體平衡模型都是全部顆粒描述性的模型，沒有涉及顆粒包覆表面的演化細(xì)節(jié)過程，而CFD-DEM模型是單顆粒層次的模擬模型，可以準(zhǔn)確描述顆粒的包覆層變化，例如可以考慮單個(gè)顆粒包覆質(zhì)量的變化，耦合顆粒表面離散化分布函數(shù)，就可以對(duì)單顆粒包覆均勻性進(jìn)行描述。另外，上述模型都是基于包覆反應(yīng)器存在不同分區(qū)這一假設(shè)基礎(chǔ)。分區(qū)概念是一種比較合理的近似，但仍是一種不精確的描述，實(shí)際上，包覆反應(yīng)器內(nèi)各個(gè)分區(qū)之間不可能有這么清晰的界限，這一點(diǎn)也始終是顆粒包覆模擬準(zhǔn)確性受質(zhì)疑的最大原因[1-2]。

基于以上對(duì)顆粒包覆過程數(shù)值模擬模型的分析，針對(duì)包覆燃料顆粒制備過程使用的流化床-化學(xué)氣相沉積技術(shù)，我們提出了CFD-DEM-CVD模型，其基本思想是去除包覆反應(yīng)器分區(qū)概念，將多組分氣體化學(xué)反應(yīng)濃度場(chǎng)、顆粒運(yùn)動(dòng)-吸附-沉積結(jié)合起來，用于分析在顆粒流化狀態(tài)下化學(xué)氣相沉積包覆反應(yīng)過程，如圖16所示。

圖16 CFD-DEM-CVD多尺度耦合模型示意圖Fig.16 Schematic diagram of CFD-DEM-CVD multiscale coupling model

在此模型中，顆粒生長(zhǎng)速率與單位時(shí)間內(nèi)顆粒掃過的體積以及顆粒周圍的有效氣體濃度成正比[17]，即

式中：R為顆粒半徑；ka為粘附系數(shù)；kr為沉積系數(shù)；粒子速度Up可以從DEM獲得，鄰近粒子的氣體速度Uf和有效氣體濃度Ceff可以從CFD獲得，從而建立了CFD-DEM-CVD多尺度耦合模型概念。

目前我們已經(jīng)實(shí)現(xiàn)了單相耦合，基于傳統(tǒng)分區(qū)概念，將該模型用于分析顆粒在單孔噴動(dòng)床的包覆區(qū)內(nèi)的停留時(shí)間分布等行為，如圖17所示。

從圖中可以清楚地看到包覆顆粒和未包覆顆粒之間的混合行為，雖然整個(gè)區(qū)域分為有效包覆區(qū)和無效包覆區(qū)，但噴嘴床中顆粒快速混合，2 s后包覆過程變成一個(gè)整體，而包覆效率僅由操作和設(shè)備參數(shù)決定。

圖17 不同包覆時(shí)間下顆粒包覆的模擬結(jié)果Fig.17 Simulation results of particle coating at different coating time

5 結(jié)論

1）顆粒包覆均勻性研究包括2個(gè)層次，即全顆粒均勻包覆和單顆粒均勻包覆，目前文獻(xiàn)研究主要集中在前者。

2）包覆過程數(shù)值模擬方法分為表面更新模型、蒙特卡羅模擬方法、分區(qū)群體平衡模型、DEM-C-PBM模型、CFD-DEM模型，其中CFD-DEM模型是單顆粒層次的數(shù)值模擬方法。

3）基于顆粒包覆數(shù)值模擬方法的比較分析，針對(duì)顆粒包覆過程的多場(chǎng)耦合特征，提出CFD-DEM-CVD多尺度模型，可以將宏觀流體尺度、介觀顆粒尺度、微觀材料沉積尺度耦合起來，進(jìn)行包覆過程數(shù)值模擬。

4）顆粒包覆過程準(zhǔn)確數(shù)值模擬的未來發(fā)展趨勢(shì)是基于單顆粒尺度的表面均勻包覆機(jī)理研究和基于反應(yīng)器尺度的物理場(chǎng)分布研究相耦合的多尺度研究方法。

（

）：

[1]TOSCHKOFF G，KHINAST J G.Mathematical modeling of the coating process[J].International Journal of Pharmaceutics，2013，457（2）:407-422.

[2]SAHNI E，CHAUDHURI B.Experimental and modeling approaches in characterizing coating uniformity in a pan coater:a literature review[J].Pharmaceutical Development and Technology，2012，17（2）:134-147.

[3]DENISA C，HEMATIA M，CHULIAB D，et al.A model of surface renewal with application to the coating of pharmaceutical tablets in rotary drums[J].Powder Technology，2003，130（1）:174-180.

[4]MANN U，RUBINOVITCH M，CROSBY E J.Characterization and analysisofcontinuousrecyclesystems[J].AICHEJournal，1979，27（5）:873-882.

[5]SHELUKAR S，HO J，ZEGA J，et al.Identification and characterization of factors controlling tablet coating uniformity in a Wurstercoatingprocess[J].PowderTechnology，2000，110（1）:29-36.

[6]KUMAR R，WASSGREN C.Inter-particle coating variability in a continuous coater[J].Chemical Engineering Science，2014，117:1-7.

[7]PANDEY P，KATAKDAUNDE M，TURTON R.Modeling weight variability in a pan coating process using Monte Carlo simulations[J].AAPS Pharm Sci Tech，2006，7（4）:E1-E10.

[8]NAKAMURA H，ABE E，YAMADA N.Coating mass distributions of seed particles in a tumbling fluidized bed coater[J].Powder Technology，1998，99（1）:140-146.

[9]KUSHAARI K，PANDEY P，SONG Y，et al.Monte Carlo simulations to determine coating uniformity in a Wurster fluidized bed coating process[J].Powder Technology，2006，166（2）:81-90.

[10]MARONGA S J，WNUKOWSKI P.Modelling of the three-domain fluidized-bed particulate coating process[J].Chemical Engineering Science，1997，52（17）:2915-2925.

[11]KUMAR R，F(xiàn)REIREICH B，WASSGREN C.DEM-compartmentpopulation balance model for particle coating in a horizontal rotating drum[J].Chemical Engineering Science，2015，125:144-157.

[12]CUNDALL P A，STRACK O D L.A discrete numerical model for granular assemblies[J].Geotechnique，1979，29（1）:47-65.

[13]ZHU H P，ZHOU Z Y，YANG R Y，et al.Discrete particle simulation of particulate systems:theoretical developments[J].Chemical Engineering Science，2007，62（13）:3378-3396.

[14]FRIES L，ANTONYUK S，HEINRICH S，et al.DEM-CFD modeling of a fluidizedbed spraygranulator[J].ChemicalEngineering Science，2011，66（11）:2340-2355.

[15]FREIREICH B，WASSGREN C.Intra-particle coating variability:analysis and Monte-Carlo simulations[J].Chemical Engineering Science，2010，65（3）:1117-1124.

[16]HILTON J E，YING D Y，CLEARY P W.Modelling spray coating using a combined CFD-DEM and spherical harmonic formulation[J].Chemical Engineering Science，2013，99（99）:141-160.

[17]LIU M L，LIU R Z，WEN Y Y，et al.Investigation of the fluidized bed-chemical vapor deposition（FB-CVD）process using CFD-DEM method[C]//2014 proceeding of the HTR.Weihai，China:2014.