劉 鵬 杜曉超 洪 鋒 袁顯寶

（三峽大學機械與動力學院 宜昌 443002）

液態(tài)鉛鉍共晶合金（Liquid Lead-Bismuth alloy，LBE）具有良好的中子學性能、抗輻射性能、熱物性和傳熱性能，是第4代反應(yīng)堆以及加速器驅(qū)動次臨界 系 統(tǒng)（Accelerator Driven Sub-critical System，ADS）的散裂靶兼冷卻劑［1-2］。鉛鉍合金冷卻劑跟其他流體一樣，會對其流經(jīng)的冷卻劑通道產(chǎn)生腐蝕沖刷、FAC效應(yīng)等［3-4］，使其結(jié)構(gòu)材料腐蝕脫落形成顆粒；另外，鉛鉍合金中不可避免地會存在一些雜質(zhì)，其中一些熔點比鉛鉍合金高的雜質(zhì)不能熔化，就會形成顆粒隨鉛鉍合金在回路中流動；綜合上面這兩種情況，若這些顆粒沉積在管壁時，不僅影響鉛鉍合金的流動，而且還會影響傳熱，造成管道局部溫度升高，這會造成極大的安全隱患。另外，鉛鉍合金開展實驗困難［5］，模擬不失為一種好的方法。華北電力大學的周濤等［6］對亞微米顆粒在Sierpinski海綿模型中熱泳沉積做了研究；劉亮等［7］對鉛鉍合金中顆粒物的沉積做了研究，并且設(shè)計了一種防鉛鉍合金中顆粒物沉積的管道；楊旭等［8］對矩形通道內(nèi)不溶性腐蝕產(chǎn)物的沉積分布進行了研究；陳娟等［9］對納米顆粒在液態(tài)金屬內(nèi)的熱泳規(guī)律進行了研究等。但是這些顆粒最終還是存在于管路中，考慮到顆粒的密度均小于液態(tài)鉛鉍合金的密度，因而考慮基于密度差，研究顆粒在液態(tài)鉛鉍合金中的上浮規(guī)律，這對于顆粒的除去研究是很有意義的。

1 物理模型

標準k-ε湍流模型是目前工程計算中應(yīng)用最為廣泛的模型，可以滿足大部分仿真需要，所以模擬過程選用標準k-ε湍流模型作為計算模型。k和ε的方程分別為：

其中：pk為湍動能k的生成項，由下式確定：

式中：μt= ρCμk2/ε；Cμ、σk、σε、C1ε和C2ε為湍流模型系數(shù)，根據(jù)Rodi建議的標準值分別?。篊μ=0.09，σk=1.0，σε=1.3，C1ε=1.44，C2ε=1.92［10］；且本文Fluent模擬計算采用此值。

本研究不考慮相間的質(zhì)量和熱量傳遞，連續(xù)相的歐拉控制方程為：

連續(xù)性方程：

式中：α、ρ、u分別表示體積分數(shù)、密度和速度，下標i表示液相。

動量守恒方程：

式中：p、μ、g、F分別表示壓力、動力粘度、重力加速度以及相間作用力。

顆粒在流場中的運動軌跡是其在流場中所受的各種作用力共同作用的結(jié)果，在DPM離散相模型中，流場中顆粒的運動軌跡是通過對拉氏坐標系下的顆粒作用力微分方程進行積分得到的，顆粒的作用力平衡方程為：

式中：FD(u-up)是單位顆粒質(zhì)量所受到的流體曳力；Fx為其他作用力，包括浮力、阻力、壓力梯度力、附加質(zhì)量力、Basset力、Saffman力和Magnus力等；其中豎直方向主要的力是重力、浮力和阻力，在運動初始階段壓力梯度力和附加質(zhì)量力也較為重要；在徑向方向（以下稱x方向）主要的力是Saffman力，在不考慮顆粒自身的轉(zhuǎn)動時Magnus力可以忽略。由于顆粒物的粒徑為微米級，因而不考慮布朗力和熱泳力等。其中：

式中：u為流體速度；up為顆粒速度；μ為流體動力粘度；ρ為流體的密度；ρp為顆粒的密度；dp為顆粒直徑；CD為曳力系數(shù)；Rep為相對雷諾數(shù)，定義為：

其中：曳力系數(shù)CD可由式（9）計算得到［11］：

式中：a1、a2、a3為常數(shù)。具體取值見表1。

表1 不同雷諾數(shù)下a取值Table 1 Value of a under different Reynolds numbers

2 研究對象與方法

為了確定顆粒上浮之后的粒徑分布與濃度分布范圍，計算顆粒的起始位置與終點位置在x方向的偏移量，以此偏移量的長度為半徑的圓形范圍即為顆粒最大概率的分布范圍。由于容器是軸對稱的圓柱形，為了更好的理解這一問題以及定量得到計算結(jié)果，同時也為了減少計算量和增加收斂性，所以將三維的模型降維成二維，如圖1所示，參考華北電力大學在液態(tài)鉛鉍合金氧控實驗中所用的裝置［1］，設(shè)置尺寸為直徑L=400 mm，高為H=500 mm；網(wǎng)格劃分采用結(jié)構(gòu)化網(wǎng)格。但是在后面進行顆粒分布的模擬時，為了便于觀察模擬計算結(jié)果，采用三維模型，并且使用結(jié)構(gòu)化網(wǎng)格進行計算，經(jīng)網(wǎng)格無關(guān)化驗證，二維網(wǎng)格數(shù)量為5萬左右，三維網(wǎng)格為60萬左右，截面的網(wǎng)格情況如圖2所示。

表2 液態(tài)LBE物性參數(shù)經(jīng)驗關(guān)系式Table 2 Empirical relationship of liquid LBE physical parameters

由于Fluent自帶的數(shù)據(jù)庫沒有鉛鉍合金，所以需要在Fluent中創(chuàng)建新的材料的物性數(shù)據(jù)庫，根據(jù)Koji等［12］、蘇子威等［13］的研究結(jié)果［14-15］，計算用的鉛鉍合金的物性參數(shù)列于表2中。代數(shù)方程的求解全部采用SⅠMPLE算法，方程的求解差分方式采用標準差分方式，湍流動能k、以及動量和能量等均采用二階迎風差分格式，這種格式使得計算的速度較快，同時又具有較好的精度和收斂性；選擇雙向耦合方式，即顆粒物能影響流場，同時流場的參數(shù)變化又會反過來對顆粒物的運動產(chǎn)生影響，顆粒物選擇為氫化鋰與氧的反應(yīng)產(chǎn)物氧化鋰（Li2O），為了比較全面的研究粒徑對氧化鋰顆粒物運動規(guī)律的影響，粒徑選擇為0.1～1.0 mm［16］，壁面設(shè)置為吸收面。

圖1 幾何模型Fig.1 Geometric model

圖2 截面網(wǎng)格Fig.2 Section grid

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

3.1 粒徑的分布

當鉛鉍合金的溫度為573.15 K時，不同直徑的顆粒在流場相同的位置，以相同的初始條件進行計算得到的結(jié)果如圖3所示，由圖3可知，在同一溫度下，不同直徑的顆粒，在x方向的速度隨時間的變化趨勢相同，首先在第一階段，速度在短時間內(nèi)達到最大值，斜率較大；接著進入第二階段，速度快速下降，達到平穩(wěn)的階段，隨著上浮到容器的頂端，計算結(jié)束，速度變?yōu)?；在第一階段時，浮力大于重力以及流體的阻力和附加質(zhì)量力的總和，所以向上做加速運動，出現(xiàn)x方向的位移是由于x方向存在薩夫曼升力的作用，大小與顆粒運動速度的平方成正比，方向為垂直于顆粒運動方向。隨著顆粒進入第二階段，上浮速度的降低，薩夫曼升力減小，所以x方向的速度也減小。

圖3 顆粒物在x方向的速度隨時間的變化Fig.3 Ⅴelocity of particles in the x direction with time

由圖4結(jié)合圖3可知，圖3的第一階段對應(yīng)于圖4的初始直線段，斜率較大，因而速度較大；直線段后為曲線段，斜率減小，對應(yīng)圖3的速度減小的第二階段，由于在同一溫度下，顆粒在x方向的速度隨時間變化的趨勢相同，所以，位移隨時間變化的趨勢也相同。從圖4也可以得出顆粒在x方向位移的范圍，介于20～35 mm。

圖4 顆粒物在x方向的位移隨時間的變化Fig.4 Displacement of particles in the x direction with time

圖5 為同一顆粒物直徑，在不同的溫度下，x方向的速度隨時間的變化，與圖3的趨勢是一致的；由圖5可知，在不同溫度的鉛鉍合金里，顆粒物在x方向的速度大小不同，相差不超過2 mm，但是趨勢相同，造成這種現(xiàn)象的原因是不同的溫度導致鉛鉍合金的密度和粘度不同，這些參數(shù)又會影響浮力與阻力以及附加質(zhì)量力的大小，從而顆粒物的上浮速度發(fā)生變化，導致薩夫曼升力大小發(fā)生變化。出現(xiàn)如圖5所示的不同溫度下顆粒物的x方向速度的不同。

圖5 顆粒物在x方向的速度隨時間的變化Fig.5 Ⅴelocity of particles in the x direction with time

圖6 為同一顆粒物直徑，在不同的溫度下，x方向的位移隨時間的變化，與圖4趨勢一致；不同的溫度會導致鉛鉍合金的密度、粘度不同，這些參數(shù)進而影響浮力與阻力以及附加質(zhì)量力的大小，導致x方向的薩夫曼升力的大小發(fā)生變化；由圖6可知，隨著溫度的升高，x方向的速度降低，速度的降低導致在x方向的位移也減小，于是出現(xiàn)圖6所示的情形。

圖6 顆粒物在x方向位移隨時間的變化Fig.6 Displacement of particles in the x direction with time

圖7 是粒徑與x方向最大位移的關(guān)系，在同一溫度下，x方向的最大位移隨顆粒物直徑的增大而減??；造成這一現(xiàn)象的原因從物理模型分析，隨著顆粒直徑的增大，浮力的增加量要大于重力的增加量和阻力的增加量，由前面的分析可知，顆粒的加速時間非常短暫，所以附加力的作用時間較短，這就導致顆粒在y方向的合力較大，因而上浮速度更快，在x方向的薩夫曼升力與阻力都與速度的平方成正比，升力系數(shù)對升力的大小影響能力有限，所以x方向的速度變化很小，這一點與圖3所示的一致，因而在上浮時間減小的情況下，x方向的位移也減小。另一方面，溫度對x方向最大位移的影響是通過影響流體的密度、粘度等物性來間接影響顆粒物的分布。通過圖7分析也可知，直徑大的顆粒物分布的靠近軸線，直徑小的顆粒物則較遠離軸線。

圖7 顆粒物x方向最大位移隨顆粒物直徑的變化Fig.7 The maximum displacement of the particles in the x direction varies with the diameter of the particles

3.2 時間的分布

圖8 顯示顆粒的上浮時間隨顆粒物的直徑的變化關(guān)系，隨著顆粒物直徑的增大，上浮時間減小，這與對圖7的分析是一致的，是由于直徑的改變導致各種力的變化，進而影響所作用的對象-顆粒，可以發(fā)現(xiàn)上浮時間的差距還是很大的；但是流體的溫度對上浮時間的影響很小，尤其是對小直徑的顆粒而言，幾乎沒有影響，當直徑大于0.5 mm時，出現(xiàn)了一定的的差別，但差距也沒有超過0.5 s。

圖8 上浮時間隨顆粒物直徑的變化Fig.8 The floating time varies with the diameter of the particles

3.3 濃度的分布

圖9 為三維模型的情況，在初始時刻容器內(nèi)的顆粒物的分布俯視圖，由于液態(tài)鉛鉍合金不透明，實驗測定顆粒物分布有困難，因而顆粒物的坐標暫時采用均勻分布隨機生成數(shù)法由Matlab隨機生成。

圖9 初始顆粒分布俯視圖Fig.9 Top view of the initial particle distribution

圖10 顯示了顆粒物濃度的分布，可以發(fā)現(xiàn)顆粒近似成對稱分布；同時，圍繞中心軸Position=0 mm，以s=90 mm為半徑的范圍內(nèi)，顆粒物的濃度要比其他位置平均高2～3倍。

圖10 顆粒物濃度分布Fig.10 Particle concentration distribution

4 結(jié)語

1）在同一溫度下，直徑大的顆粒物上浮時間短，直徑小的顆粒物上浮時間長，且溫度對上浮時間的影響較小，可以忽略；顆粒物的直徑與其在x方向的最大位移成負相關(guān)，即直徑大的顆粒物在x方向的位移小，直徑小的顆粒物在x方向的位移大。

2）在同一顆粒物直徑下，溫度與顆粒物在x方向的最大位移成負相關(guān)，即溫度高時x方向的位移小，溫度低時x方向的位移大。

3）直徑大的顆粒物分布靠近圓柱軸心位置，直徑小的顆粒物分布相較直徑大的顆粒物分布遠離軸心位置。

4）總的來說，圍繞軸心半徑為90 mm左右的圓內(nèi)，顆粒物的濃度較其位置要高2～3倍，因此，在對鉛鉍合金進行凈化時應(yīng)著重關(guān)注這一范圍，凈化的效率與效果較好。