微米颗粒污垢特性的数值模拟

徐志明 张敏霞 张一龙 王景涛

(1. 东北电力大学能源与动力工程学院；2. 华北电力大学能源与机械工程学院)

符号说明

CD——曳力系数；

d——颗粒直径，m；

f——壁面摩擦系数；

Fa——颗粒浮力，N；

Fp——颗粒重力，N；

Fr——颗粒曳力，N；

Fs——颗粒Saffman力，N；

m——质量；

n——猝发次数；

N——颗粒个数；

Nd——在壁面上沉积的颗粒个数；

Re——雷诺数；

Rf——污垢热阻，m2·K/W；

S——广义源项；

t——沉积时间，s；

td——颗粒沉积时间；

v——速度，m/s；

V*——壁面摩擦速度，m/s；

α——常数；

Γ——广义扩散系数；

θ——粒子湍流猝发时间，s；

λf——污垢导热系数，W/(m·K)；

μ——流体运动粘度；

ν——流体动力粘度，Pa·s；

ρ——密度，kg/m3；

τs——壁面剪切应力，N；

φ——广义变量；

下标：

d——沉积；

f——流体(水)；

p——颗粒。

换热器运行过程中形成的污垢直接影响换热器的换热效率。颗粒污垢是指悬浮在流体中的固体颗粒在换热面上的积聚。目前在污垢模拟中对析晶污垢的模型已经很成熟，但是还没有较好的模型来描述颗粒污垢，颗粒污垢沉积在换热器中特别普遍，具有广泛的工业背景，该过程包含复杂的两相流动和传质理论，对其进行深入研究既有重要的理论意义，又有实际除垢、抑垢的应用价值。

国内外在颗粒污垢方面已经做了大量的研究工作。Chamra L M和Webb R L研究了强化管中液侧颗粒污垢的形成机理，首次基于扩散机制和惯性机制提出一个半理论模型，通过实验数据得到粘附概率和污垢抗分散强度，以此预测较大粒径和浓度分布范围内的颗粒污垢沉积机理[1]。Beal S K重点研究了颗粒污垢的形成过程，得到颗粒污垢沉积机理主要有：布朗运动、湍流扩散和重力或惯性力，且颗粒的沉积量与流体介质粘性有关[2]。徐志明和王建国基于湍流猝发理论，提出了一个以易测参数表述的污垢热阻预测模型，其形式类似于基于Kern-Seaton模型得到的结果，但待定系数少、所需参量更易测量[3]。甘云华和徐志明利用这一模型研究了管壳式换热器管内微粒污垢的积聚特性，发现模型中引入的输运速度经验公式在输运机制发生变化时误差较大[4]。

对于颗粒沉积的数值模拟研究中连续相的模拟大多采用直接数值模拟(DNS)[5，6]、大涡模拟(LES)或雷诺时均模拟(RANS)；DNS和LES计算量大，现在研究较少，而RANS除了计算量小且方便外，还能进行大雷诺数的模拟，在现代计算中应用比较广泛，离散相大都采用拉格朗日和欧拉模型。Sun L等数值模拟了血浆蛋白颗粒在层流下的沉积过程，通过计算指出在层流条件下，颗粒的沉积量与速度无关，与颗粒数量、大小、沉积时间和沉积高度有关[7]。曾敏捷等采用RNGκ-ε模型描述气体的流动，并在拉格朗日框架下追踪颗粒运动，建立了一个新的完整的口腔-咽-喉-气管-支气管三维几何模型[8]，结果表明，颗粒沉积率主要取决于气流流速与颗粒惯性。Martos C等提出了一种新的微米颗粒沉积模型，该模型基于粒子在壁面附近湍流运动的轨迹，得出在垂直和水平管道内颗粒沉积速度的V形变化图[9]。Tian L比较了在湍流中不同模型下的微纳米颗粒沉积，详细对比了RSM、κ-ε、CFWN和DRW各模型的差异，从结果中看出RSM模拟结果相比其他模型更接近直接模拟[10]。朱静等采用雷诺应力结合拉格朗日模型模拟碳酸钙颗粒在圆管内的沉积，比较脉动速度模拟结果，证明了三维模拟结果优于二维模拟结果，且在近壁区内十分接近DNS结果[11]。刘洪涛和张力数值模拟了颗粒在气相中的趋壁沉积过程，假设微米级颗粒一旦附着在壁面上就很难使其从壁面脱离，从而得出了在不同弛豫时间下的无量纲沉积速度[12]。上述对于颗粒沉积的模拟大多只考虑了颗粒的沉积过程，没有涉及到颗粒在沉积过程中还存在剥蚀过程。实际过程中沉积与剥蚀同时存在，假定颗粒一旦附着壁面就很难脱离不能够完全说明颗粒沉积的整个过程，特别是在固液两相中，流体相对颗粒的剥蚀是客观存在的且必须考虑。

笔者的研究内容主要包括：分析固体颗粒在液体中的沉积机理，在比较前人计算颗粒沉积方法的前提下，改进得出了适合本课题的颗粒模型，即基于流体流动的前提下，采用雷诺应力模型和拉格朗日方法建立合适的数学模型，同时考虑颗粒的沉积与剥蚀来研究微米颗粒污垢的污垢热阻；通过运用Fluent软件数值模拟微米氧化镁颗粒在圆管内的沉积过程，得到颗粒的沉积率；结合污垢的沉积机理，通过引入经验公式结合颗粒的沉积率计算得到颗粒的剥蚀率，得到微米氧化镁颗粒污垢的污垢热阻，并与实验数据对比，验证模型的可行性，并进一步分析微米颗粒污垢的沉积特性。

1 数学方程

1.1流体相方程

若要实现颗粒污垢的数值模拟，首先要得到管内的流场及温度场等的分布，这就需要求解连续性方程、动量方程、能量方程和湍流方程。这些方程均可以表示成一种简洁的通用形式[13]：

(1)

Ahmadi G和Chen Q发现κ-ε模型忽略了湍流的各向异性，导致颗粒沉积速度的计算结果偏大，虽然RSM考虑了湍流的各向异性，但还是会得出偏大的结果[14]；Gao N等比较了不同RANS模型与拉格朗日方法的结合对计算颗粒沉积的影响[15]，并指出其他模型的缺陷以及RSM、拉格朗日方法[16]在计算成本和计算精度上的适用性。

Tian L发现用标准壁面法计算得到的结果会偏大，采用增强壁面法比标准壁面法算出的结果更接近直接模拟值[10]。

1.2颗粒沉积模型

很多学者研究过颗粒在湍流中的运动[17～21]，一般来说，颗粒受到的作用力主要包括曳力、压力梯度力、重力、浮力、布朗力及Saffman升力等，在各种不同的具体条件下，颗粒运动方程中的某些作用力项可以合理忽略，笔者考虑颗粒所受的作用力包括重力、浮力、曳力和Saffman升力。单个颗粒的运动方程可直接根据牛顿第二定律得出：

(2)

颗粒运动所受到的重力为：

(3)

浮力为：

(4)

曳力为：

(5)

当雷诺数不太大时，CD可以表示为：

(6)

采用Stokes标准曳力公式：

f(Rep)=1

(7)

颗粒在有速度梯度的流场中运动时，即使没有受到旋转也会受到一个附加侧向力，即Saffman力。在小的颗粒雷诺数下可以表示为：

(8)

在计算中，考虑曳力、重力、浮力和Saffman升力而忽略其他作用在固体颗粒上的力，得到以下方程：

(9)

在Fluent中选择Discrete phase model模型和基于拉格朗日法的Discrete Random Walk (DRW)模型，此模型用来预测湍流影响下的颗粒运动。通过壁面捕捉颗粒数可以得到颗粒的沉积率：

(10)

1.3颗粒污垢剥蚀模型

文献[22]指出粒子剥离是粘性底层中湍流猝发的结果，可导出粒子剥离率的表达式。其中粒子湍流猝发时间为：

θ=n(75ν/V*2)

(11)

从而得到粒子剥离率即剥蚀率为：

(12)

对以上公式进行求导，当时间为0时，可得到：

(13)

式(13)表明，初始剥离率仅取决于壁面剪切应力和运动粘度系数。其中α与流动性质和再沉积有关。考虑到猝发中心到外缘其强度逐步衰减，根据Cleaver和Yates的假定α≈0.01。壁面摩擦速度V*由壁面摩擦系数f求得，忽略管壁的粗糙，可得到光滑管摩擦系数公式：

(14)

由此可计算壁面摩擦速度：

(15)

τs=fρv2/2

(16)

根据文献[22]提出的污垢模型描述可得到污垢沉积的净沉积率：

(17)

如果假定污垢的成分和特性沿换热面和沿垢层厚度方向都是均匀分布的，则热阻为：

Rf=mf·(ρfλf)-1

(18)

2 模型建立

2.1物理模型

笔者的物理模型为光滑圆管三维模型。物理模型如图1所示，管长1 500mm，管直径为22mm，流体与颗粒从左端面流入，从右边流出。

图1 物理模型

2.2网格划分和边界条件

网格划分采用Hex/Wedge、Cooper网格，计算区域3个方向x×y×z的网格划分为1500×22×22(图2)。边界条件设定为：入口inlet设置为VELOCITY_INLET型入口，出口outlet为PRESSURE_OUT型出口；管壁壁面wall设置为WALL型壁面。

图2 网格划分

笔者选用的颗粒为20μm氧化镁颗粒，氧化镁颗粒和流体的物性如下：

氧化镁颗粒密度 1 740kg/m3

氧化镁颗粒直径 0.02mm

流体(水)密度 998kg/m3

流体(水)动力粘度 0.805μPa·s

3 模拟结果与分析

3.1网格无关性验证

物理模型网格无关性验证结果如图3所示，设定工况为微米氧化镁颗粒在光滑圆管内沉积，研究入口浓度200mg/L、入口温度30℃、入口速度0.2m/s时不同网格数对沉积率的影响。在边界条件设置中，将四周壁面设为颗粒捕捉类型，当颗粒碰到管壁的时候，管壁就会把颗粒捕捉，相当于颗粒在壁面上沉积。在模拟计算中每隔1h取一个点，由沉积颗粒的个数来计算沉积率，当划分网格数从580 800增大到800 000时，求解精确度有所提升，这说明网格加密对模拟结果有影响。当网格数不小于726 000时，沉积率基本不再变化，即继续加密网格对计算精确度的提高已无明显影响。结合文中的模型和上述考虑，选择网格数为726 000。得到的沉积率是一些离散的点，先增加然后趋于不变，在某一值附近上下波动。

但在实际颗粒污垢形成过程中，不是所有颗粒碰到壁面就是沉积，其中还包括颗粒受到流体剪切力而导致颗粒的剥蚀，所以在模拟计算污垢热阻时，将污垢分为两部分进行计算，首先通过Fluent模拟固液两相流动计算颗粒的沉积率，然后通过经验公式用沉积率来计算剥离率，最后通过净沉积量计算得到污垢热阻值。

图3 网格无关性验证

3.2实验验证

图4、5为模拟数据与文献[23]中的实验数据的对比和误差分布。图4中的工况为：入口温度30℃、水浴温度50℃、流速0.2m/s，氧化镁浓度400mg/L。图5工况为：入口温度30℃、水浴温度50℃、流速0.2m/s，氧化镁浓度200mg/L。将实验值与模拟值每1h取一个点作对比，从图4b、5b可以看出，刚开始第一个点的相对误差较大，颗粒污垢在结垢的开始阶段会存在微小的诱导期，如果忽略初始阶段微小的诱导期，其他时间的相对误差区间均在20%以内，符合规范的要求，由此验证笔者所建立的颗粒污垢模型是可行的。

图4 模拟与实验(工况一)对比和误差分布

3.3流速对污垢热阻的影响

在其他工况不变的情况下，只改变流速，对比流速分别为0.20、0.25、0.30、0.35、0.40m/s时的污垢特性。不同流速下圆管的污垢特性曲线如图6所示。从图6中可以看出，结垢速率随着流速的增加而减小，污垢热阻渐近值随着流速的增加而明显降低。从计算的两方面看，速度越大，单位时间内到达壁面的颗粒数越多，沉积率越大。由于流速的增大，使得壁面摩擦速度增大，从而导致剥蚀率的增加，且剥蚀率的增长幅度比沉积率更大，导致曲线会有交点，使得速度越大，污垢热阻渐近值会越小。

图5 模拟与实验(工况二)对比和误差分布

图6 流速对圆管污垢特性的影响

3.4浓度对污垢热阻的影响

在其他工况不变的情况下，只改变微米氧化镁浓度，研究浓度对圆管污垢特性的影响，模拟的5组浓度分别为100、200、300、400、500mg/L。不同浓度下的污垢热阻曲线如图7所示，可以看出随着溶液中微米氧化镁浓度的增大，污垢热阻渐近值逐渐增大，并且结垢速率也随之明显加快。这是由于，当颗粒浓度增大时，单位容积中会存在更多的微米氧化镁颗粒，使得被输运到管内壁的颗粒随着浓度的增大而增多，颗粒沉积率随之增加。而微米氧化镁浓度的增加对剥蚀的影响相对较小，当颗粒的沉积量与剥蚀量相当时，污垢热阻渐近值不再发生明显变化。

图7 浓度对圆管污垢特性的影响

3.5入口温度对污垢热阻的影响

在其他工况不变的前提下，只改变入口温度来研究温度对圆管颗粒污垢的影响，入口温度分别为20、25、30、35、40℃。不同入口温度对颗粒污垢热阻的影响如图8所示，可以看出，随着循环工质入口温度的升高，达到渐近值的时间变短。从图8中看出温度的变化对污垢热阻有影响，在计算过程中发现，温度变化对颗粒趋壁沉积影响不大，主要是温度的变化使得流体的粘度系数也随之改变，从而改变颗粒的剥蚀率，温度升高，剥蚀率增加，污垢热阻渐近值也明显减小。

4 结论

4.1将雷诺应力和拉格朗日方法结合，引入颗粒污垢的沉积率和颗粒污垢的剥蚀率可以计算颗粒污垢热阻，通过与实验数据的对比，发现模拟数据具有较高的精度，证明此方法是可行的。

图8 入口温度对圆管污垢特性的影响

4.2随着流速的增加，圆管颗粒污垢热阻渐近值有所降低；随着颗粒浓度的增加，结垢速率加快，污垢热阻渐近值明显增大；随着入口温度的升高，污垢热阻渐近值逐渐降低。

[1] Chamra L M，Webb R L.Modeling Liquid-side Particulate Fouling in Enhance Tubes[J].International Journal of Heat and Mass Transfer，1994，37(4)：571～579.

[2] Beal S K.Depsition of Particles in Turbulent Flow on Channel or Pipe Walls[J]. Nuclear Science and Engineering，1970，40(1)：1～11.

[3] 徐志明，王建国.微粒污垢剥蚀机制研究[J].工程热物理学报，1998，19(5)：611～615.

[4] 甘云华，徐志明.微粒污垢的特性分析[J].热科学与技术，2002，1(2)：109～112.

[5] Kuerten G M.Can Turbophoresis be Predicted by Large Eddy Simulation[J].Physics of Fluids，2005，17(1)：11.

[6] Zhang J，Li A.CFD Simulation of Particle Deposition in a Horizontal Turbulent Duct Flow[J].Chemical Engineering Research and Design，2008，86(1)：95～106.

[7] Sun L，Lin J，Bao F.Numerical Simulation on the Deposition of Nanoparticles under Laminar Conditions[J].Journal of Hydrodynamics Ser. B，2006，18(6)：676～680.

[8] 曾敏捷，胡桂林，樊建人.微颗粒在人体上呼吸道中运动沉积的数值模拟[J].浙江大学报(工学版)，2006，40(7)：1164～1167.

[9] Martos C，Coto B，Pea J L，et al.Effect of Precipitation Procedure and Detection Technique on Particle Size Distribution of CaCO3[J].Journal of Crystal Growth，2010，312(19)：2756～2763.

[10] Tian L.Particle Deposition in Turbulent Duct Flows——Comparisons of Different Model Predictions[J].Aerosol Science，2007，38(4)：377～397.

[11] 朱静，张红，陶汉中.圆管内水中碳酸钙颗粒趋壁沉积数值模拟[J].化学工程，2013，41(6)：58～62.

[12] 刘洪涛，张力.微细颗粒壁面沉积的数值研究[J].工程热物理学报，2010，31(3)：431～434.

[13] 陶文铨.数值传热学[M].西安：西安交通大学出版社，2001：203～408.

[14] Ahmadi G，Chen Q.Deposition of Particles in a Turbulent Pipe Flow[J].Journal of Aerosol Science，1997，28(5)： 789～796.

[15] Gao N，Niu J，He Q，et al.Using RANS Turbulence Models and Lagrangian Approach to Predict Particle Deposition in Turbulent Channel Flows[J].Building and Environment，2012，48(2)：206～214.

[16] 汪靓，崔小朝，蔡明，等.拉格朗日-雷诺应力模型在后台阶颗粒流中的运用[J].太原科技大学学报，2011，32(5)：406～409.

[17] 朱辉.粉尘颗粒在纤维过滤介质表面沉积的计算机模拟[D].上海：东华大学，2009.

[18] Dehbi A.A CFD Model for Particle Dispersion in Turbulent Boundary Layer Flows[J].Nuclear Engineering and Design，2008，23(3)：707～715.

[19] Izadi M.Experimental and Numerical Investigation of Fouling in Heat Exchangers[D].New York：Clarkson University，2011.

[20] Wood N B.A Simple Method for the Calculation of Turbulent Deposition to Smooth and Rough Surfaces[J].Journal of Aerosol Science，1981，12(3)：275～290．

[21] Brahim F，Augustin W，Bohnet M.Numerical Simulation of the Fouling Process[J].International Journal of Thermal Sciences，2003，42(3)：323～334.

[22] 杨善让，徐志明，孙灵芳.换热设备污垢与对策[M].北京：科学出版社，2004：149～154.

[23] 王景涛.交叉缩放椭圆管颗粒污垢特性的实验研究[D].吉林：东北电力大学，2014.