顆粒包覆過程的數值模擬方法比較研究

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

（清華大學核能與新能源技術研究院；先進核能技術協同創新中心，北京100084）

顆粒包覆過程，即在顆粒的表面均勻包覆上一層功能材料，具有保護、緩釋、提高環境相容性等用途。在工業生產中，尤其是制藥行業應用廣泛，同時也是高溫氣冷堆核燃料元件制備過程的重要環節。顆粒包覆是顆粒性能改進的主要方法之一，而包覆層的均勻性是顆粒包覆過程追求的目標。高溫氣冷堆核燃料元件采用的是4層包覆燃料顆粒，因為核燃料的固有安全特征設計，要求包覆層滿足一定的技術指標，其中顆粒包覆層的均勻性則是決定高溫氣冷堆安全性的關鍵因素之一。

實驗研究顆粒均勻包覆工藝是很重要的一個方面，如果通過數學建模和數值模擬，能分析出操作參數對顆粒均勻包覆的影響規律，可指導實驗研究方向，減少實驗次數。對顆粒包覆過程進行準確的數值模擬是研究均勻包覆機理的一種有效手段，進而可優化包覆過程，進行包覆工藝的放大，也是提高顆粒包覆過程經濟性的有效途徑。如何進行包覆過程的準確數值模擬，以往的研究給出了很多顆粒包覆過程數值模擬方法，但缺少比較系統和深入的對比研究[1]。

本文中首先分析了顆粒均勻包覆過程的主要特征，指出顆粒包覆過程的研究層次，即單顆粒表面均勻和全顆粒均勻2個層次，指出目前大量的模型僅局限在顆粒之間包覆層質量均勻的模擬研究。然后闡述不同數學模型的具體實施過程，指出其優勢和劣勢，重點對CFD-DEM方法包覆過程模擬進行深入分析，并基于運動－吸附－沉積機理發展新型包覆過程模擬方法。最后基于現有文獻中的顆粒包覆過程數值模擬方法對比分析研究，指出包覆模擬研究的規律，給出了顆粒包覆過程準確數值模擬的發展趨勢，即基于單顆粒尺度的表面均勻包覆機理研究和基于反應器尺度的物理場分布研究相耦合的多尺度研究方法。

1 顆粒包覆技術

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

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

2 顆粒包覆層均勻性分析

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

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

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

3 數值模擬方法

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

3.1 全顆粒包覆均勻數值模擬

3.1.1 表面更新模型

顆粒包覆的一般過程是顆粒在不斷運動的過程中循環往復通過含有待沉積產物的區域，待沉積產物附著在顆粒表面形成包覆層[3]，基本原理見圖3。

表面更新模型是一個基礎的數學分析模型，起源于連續循環體系的數學分析[4]，很早被用于分析顆粒包覆工藝參數的影響，其基本原理認為反應器內分為包覆區和非包覆區，顆粒每通過一次包覆區質量增加為xi，則通過n次后總的包覆質量為：

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

顆粒包覆質量均勻性，即變化系數cv，m（coefficient of variation）可以用顆粒包覆質量的標準偏差來表示，

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

經過數學推導，可得到顆粒包覆層質量的變化系數和單顆粒循環時間以及每次包覆量的關系

式中：μcpp和σcpp是在一次循環時間內全體顆粒每次通過包覆區的包覆層增加值的平均值和標準偏差；μct和σct分別是所有顆粒循環時間的平均值和標準偏差；Tcoat是顆粒總的包覆時間。

顆粒包覆層質量分布如圖4所示。

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

圖5 顆粒循環時間分布及每次通過包覆區質量增量分布Fig.5 Particle circulation-time distribution and coating-per-pass distribution

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

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

圖6 利用蒙特卡羅方法計算得到的顆粒循環時間分布圖Fig.6 Particle circulation-time distribution calculated by Monte Carlo method

圖7 包覆時間對顆粒包覆層質量變化的影響Fig.7 Effect of coating time on mass coating variability

3.1.2 蒙特卡羅模擬方法

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

基于該顆粒循環時間分布圖，通過實驗測得的噴霧通量分布、噴霧區形狀及面積大小、顆粒投影面積可以計算得到顆粒包覆質量均勻性變化系數與包覆時間之間的關系，如圖7所示。由圖可以看出，顆粒包覆質量均勻性變化系數cv，m與包覆時間函數1/t0.5成正比關系。

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

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

3.1.3 分區群體平衡模型

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

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

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

式中：f為F區群體密度函數；g為G區群體密度函數；e為E區群體密度函數；α為噴射的活性區占的分數；β為包覆過程整個床占的分數；r為包覆層質量增加速率；t代表時間；w代表包覆層質量。

通過上述模型可以研究包覆區的大小、包覆時間、單位時間的包覆增加量對包覆均勻性的影響規律。研究認為包覆質量是包覆時間的函數，考慮了顆粒流動的非均勻性帶來的包覆變異系數變化。分區模型關鍵參數包括噴射區的相對尺寸、分區之間的流動速率、噴射速率以及噴射時間等。文獻[10]中得到結果如圖10所示。

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

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

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

圖11 旋轉鼓室中顆粒流動示意圖Fig.11 Schematic diagram of flow of particles in a rotating drum

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

圖12 旋轉鼓室流動分區模型示意圖Fig.12 Schematic diagram of flow of compartmentpopulation particles in a rotating drum

3.1.4 CFD-DEM耦合模型

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

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

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

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

3.2 單顆粒包覆均勻數值模擬

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

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

假定顆粒在均勻的噴射場內，而顆粒的轉動是隨機的，因此與均勻流場碰撞的面是包覆層增加的面，即

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

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

式中：n是外延法向量，如圖14b所示；urel是顆粒與噴射場的相對運動速度。

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

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

式中：f為顆粒表面沉積產物分布函數；g為正態分布函數。式子左邊可以用球諧函數來表示，即得到顆粒表面沉積產物的分布。

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

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

4 模型比較分析

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

基于以上對顆粒包覆過程數值模擬模型的分析，針對包覆燃料顆粒制備過程使用的流化床-化學氣相沉積技術，我們提出了CFD-DEM-CVD模型，其基本思想是去除包覆反應器分區概念，將多組分氣體化學反應濃度場、顆粒運動-吸附-沉積結合起來，用于分析在顆粒流化狀態下化學氣相沉積包覆反應過程，如圖16所示。

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

在此模型中，顆粒生長速率與單位時間內顆粒掃過的體積以及顆粒周圍的有效氣體濃度成正比[17]，即

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

目前我們已經實現了單相耦合，基于傳統分區概念，將該模型用于分析顆粒在單孔噴動床的包覆區內的停留時間分布等行為，如圖17所示。

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

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

5 結論

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

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

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

4）顆粒包覆過程準確數值模擬的未來發展趨勢是基于單顆粒尺度的表面均勻包覆機理研究和基于反應器尺度的物理場分布研究相耦合的多尺度研究方法。

（

）：

[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，FREIREICH 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.