一种用于模拟湍流中沉积的耦合组分输运、传热与移动边界方法

《International Journal of Heat and Mass Transfer》：A coupled species transport, heat transfer and moving boundary approach for modeling deposition in turbulent flow

　　

**摘要**

表面沉积是管道和传热设备中的主要挑战，因为不断增长的沉积层会同时限制流动、增加热阻、降低工艺效率，并增加能耗和运行成本。大多数现有的连续介质模型仅估算沉积质量，而不将沉积层表示为演化的物理壁面；因此，它们无法一致地捕捉沉积层对近壁流动、界面温度和传热的反馈作用。本研究提出了一种基于物理的数值框架，将颗粒输运和表面附着与流体—沉积层界面的连续演化直接耦合。沉积层被表示为一种移动浸没边界，根据局部沉积质量增长，无需反复重新网格划分。一个次级计算网格采用统一的能量方程，同时求解流体中的对流和不断增长的固体沉积层中的热传导。这使得可以直接计算演化中的沉积层表面温度，而该温度决定了流体属性和温度相关的附着速率。虚拟时间策略进一步使长期沉积层增长模拟在可行的计算成本下得以实现。采用单一固定参数集，该框架在流速、壁温和沥青质浓度变化条件下复现了实验沉积速率，归一化均方根误差分别约为5.1%、1.2%和14.0%。统一的热解算还复现了参考共轭传热结果（在常温和常热通量条件下），最大局部偏差约为8%。随着沉积层增长，模型捕捉到了热阻的逐步增加、流体温度的降低，以及流体和固体区域中温度分布的演化。尽管沉积层表面持续位移，预期的壁面尺度速度关系得以保持，表明移动边界处理方法维持了物理一致的近壁行为。增加流速抑制了沉积，因为附着的减弱超过了向壁面输运的增强；而增加表面温度和颗粒浓度则促进了沉积。与以往基于组分的方法相比，本方法的关键进步在于统一求解了质量传递、表面附着、共轭传热和沉积层引起的几何演化。因此，该框架能够获取演化中的界面温度和近壁量，这对于物理一致地预测沥青质结垢的耦合水力与热学后果是必需的。

**1. 引言**

沉积被定义为颗粒在表面上的沉降和积累[1]，是导致工艺设备性能退化和运行成本增加的主要现象之一。沉积层的生长可能起始于单个颗粒与表面的相互作用，并逐渐改变甚至阻碍流动通道[2]。这一现象的重要性在原油输送和加工系统中尤为突出，因为沥青质沉积改变了有效流动几何形状，并同时影响流体力学、质量传递和传热。Shilpa等人研究了热耦合纳米润滑剂流动中的热泳颗粒沉积，表明沉积可强烈响应于局部热质传递条件[3]。Srilatha等人进一步表明，纳米颗粒聚集可以改变流动和热输运，突出了即使在未显式建模沉积本身的情况下，颗粒群体效应的重要性[4]。这些耦合效应最终表现为压降增加、传热效率降低以及泵送和能源需求增加，使得对沉积层增长的可靠预测对于热流体设备的设计和操作至关重要。

从机理角度来看，净沉积层生长受控于颗粒向壁面的输运、随后的附着以及可能的脱落或剥离。Kern和Seaton提出了早期的渐近结垢模型，其中恒定的沉积贡献与随层发展而增加的脱落项相竞争[5]。Watkinson随后强调了粘附概率，将颗粒到达近壁区域与稳定表面附着的概率区分开来[6]。Han等人将沉积和脱落机制纳入欧拉CFD框架，并考虑了布朗运动、湍流扩散和重力效应[7]。这些研究确立了到达壁面并不意味着永久沉积；最终生长速率是输运、附着和脱落之间竞争的结果，而每种机制都受局部界面和流体力学条件的影响。

温度敏感附着通常通过活化能或Arrhenius型关系来表示。Jamialahmadi等人在絮凝沥青质沉积中通过实验区分了输运和稳定附着，并量化了粘附过程强烈的温度依赖性[8]。Maddahian等人随后采用了一种Arrhenius型局部附着速度，将温度敏感性与壁面剪应力的抑制效应相结合，用于沥青质结垢模型[9]。在更广泛的热流体输运研究中，Reddy等人在磁流体混合纳米流体流动中使用了Arrhenius活化能模型[10]，而Nagapavani等人则在非定常纳米流体系统中结合了Arrhenius活化能和二元化学反应模型[11]。后几项研究并不直接描述沥青质粘附，但它们展示了活化输运闭合在温度敏感质量传递过程中的广泛应用。对于沉积建模，主要限制仍然是有效附着参数具有系统特异性，必须针对所考虑的原油表面体系进行论证。

随着计算能力的进步，沉积建模已从零维关联方法发展到基于场的数值描述，包括组分输运、多流体、界面捕捉和欧拉—拉格朗日公式化方法。Bayat等人使用组分输运框架预测原油结垢和热阻演化，通过极高的粘度表示沉积相[12]；这种方法有效地求解了标量输运，但没有将沉积层重建为独立移动的物理壁面。Guha综述了颗粒输运和沉积机制，说明了为什么连续介质输运描述适用于足够小且扩散主导的颗粒，而较大颗粒则需要显式处理惯性和额外的输运机制[13]。Maddahian等人使用多流体方法复现了沥青质结垢趋势，具有显式相耦合，但以额外的相方程和半经验附着闭合为代价[9]。最近的研究还探索了强耦合热质传递系统的物理信息方法：Sanju等人将物理信息神经网络应用于化学触发的双重扩散对流[14]，Sanju等人随后将此策略扩展到具有交叉扩散和重力调制的三重扩散Jeffrey流体输运[15]。后几项研究并非沉积模型，但它们说明了紧密耦合输运物理和可重复数值闭合日益增长的重视程度。

显式表示沉积层或界面演化的方法提供了额外的几何信息，但也引入了其他限制。Sileri等人使用两流体公式建模壁面沉积和老化，但计算是从现有沉积层开始的，而不是从初始阶段复现生长过程[16]。Zheng等人开发了一种机理结垢模型，其中沉积和脱落过程产生了非线性沉积层生长行为，但该公式仍然具有应用特异性[17]。Khodadadi和Maddahian将欧拉—拉格朗日颗粒描述与浸没边界方法耦合以复现物理结垢层生长；该方法提供了详细的颗粒级信息，但需要显式颗粒追踪，因此在长期沉积模拟中产生更大的计算成本[18]。代表性方法的主要优势和局限性总结于表1。

**表1a** 代表性沉积建模方法及其相对于本目标的局限性比较。

| 代表性研究 | 方法 | 主要贡献 |
|---|---|---|
| Kern和Seaton [5] | 半经验沉积—脱落模型 | 建立了沉积与脱落之间的竞争及渐近结垢行为。 |
| Bayat等人 [12] | 组分输运/高粘度沉积表示 | 以相对较低的计算成本求解标量输运和热阻演化。 |
| Sileri等人 [16] | 两流体壁面沉积与老化 | 在两流体框架中表示沉积老化和壁面积累。 |
| Han等人 [7] | 欧拉沉积—脱落模型 | 包含布朗、湍流、重力、沉积和脱落效应。 |
| Maddahian等人 [9] | 多流体沥青质结垢模型 | 以耦合相输运复现实验沥青质沉积趋势。 |
| Zheng等人 [17] | 机理结垢模型 | 展示了沉积和脱落机制如何产生非线性结垢生长。 |
| Khodadadi和Maddahian [18] | 欧拉—拉格朗日+IBM | 将颗粒沉积转化为物理浸没边界生长。 |
| 本研究 | 组分输运+IBM+次级网格传热 | 局部沉积质量直接驱动物理壁面生长；流动、界面温度和固体热传导反馈至后续沉积。 |

**表1b** 模拟管道内沥青质沉积层生长所使用的原油温度相关性质[8]。

| 性质 | 单位 | 表达式 |
|---|---|---|
| 流体密度 $\rho_f$ | kg·m?3 | $1.230\times10^{-1}T + 820.5$ |
| 流体粘度 $\mu_f$ | mPa·s | $4.5029\times10^{-4}T^2 - 3.6407\times10^{-1}T + 74.38$ |
| 流体比热容 $c_{p,f}$ | J·kg?1·K?1 | $4.14T + 0.83$ |
| 流体导热系数 $k_f$ | W·m?1·K?1 | $1.224\times10^{-4}T + 0.16299$ |

尽管可用方法多样，长期湍流沉积模拟仍存在一个关键缺口：在许多欧拉和基于组分的方法中，沉积物质通过浓度、体积分数或源项来表示，而非作为施加边界条件的物理移动壁面。如图1所示，这一区别很重要，因为壁面剪应力和壁面尺度量（如 $u^+$ 和 $y^+$）取决于流体—固体界面的实际位置。动态网格方法可以表示这种运动，但反复变形或重建可能在长期增长过程中增加计算成本并降低网格质量。浸没边界方法（IBM）在固定背景网格上表示复杂或移动边界，提供了一种替代方案。Goncharuk等人开发了基于SIMPLE的浸没表面流动IBM[19]。Lauber等人展示了由移动薄膜驱动的流动浸没边界模拟[20]。Khodadadi和Maddahian将IBM应用于沉积引起的结垢层生长[18]。在更广泛的OpenFOAM背景下，Bai等人使用大涡模拟复现了复杂溢洪道几何中的速度和湍流动能分布，证明了OpenFOAM适用于强几何相关湍流流动计算[21]。然而，这些进展本身并不能提供本问题所需的封闭反馈环路——即连接组分输运、局部附着、移动沉积边界生长和热反馈的环路。

下载：下载高分辨率图像（284KB）下载：下载全尺寸图像

**图1** 浸没边界（IB）方法与类VOF方法在模拟沉积层生长等移动边界问题中的原理对比。

在本研究中，这些局限性通过在开源OpenFOAM框架内将组分输运方程与浸没边界方法耦合来解决。沉积形成沥青质的输运通过非定常对流—扩散方程表示，而稳定表面附着则通过一阶边界条件施加。Ruckenstein提供了布朗颗粒在颗粒表面相互作用力下沉积的早期机理描述[22]，Guha确立了小尺度扩散主导颗粒连续介质描述的更广泛输运基础[13]。在本公式中，局部沉积质量不仅保留为标量或固定场；它被直接转化为浸没边界的法向位移。因此，沉积层几何形状作为局部输运和附着求解的结果而演化，且结果壁面位置反馈至后续流体力学和沉积计算。

能量方程的求解同样重要，因为沉积层增长引入了递增的固体热阻，并持续改变流体—沉积层界面温度。与静止共轭传热问题不同，固体沉积区域随时间增长并改变几何形状。因此，本框架引入一个次级计算网格，在该网格上使用统一的能量方程求解流体对流和通过不断增长的固体层的热传导。所得界面温度被映射回主IBM网格，并用于更新温度相关的流体属性、分子扩散系数和附着系数。这在沉积层生长与后续沉积之间建立了双向热反馈。

第二个数值挑战是湍流流动计算所需的短时间步长与结垢生长所需长物理时间尺度之间的分离。为使长时模拟在计算上可行，采用了一种虚拟时间策略，其中沉积引起的边界位移通过受控放大因子进行推进。该因子通过与实验沉积趋势进行比较来评估，以便在不改变主要物理响应的情况下加速生长时间尺度。

本框架的方法论区别应予以强调。其新颖性并非来自组分输运、浸没边界方法或移动边界技术的单独使用，而是来自它们为沉积引起的壁面演化而进行的直接耦合。在常规组分输运公式中，沉积物质通常通过浓度、相分数、源项或固定区域场来表示，因此不一定被重建为演化的物理壁面。在本方法中，由欧拉组分输运方程和表面附着条件预测的局部沉积质量被直接转化为浸没边界的法向位移。因此，所得沉积层表面根据局部求解的输运和附着场在空间上增长。

与体网格动态贴合方法相比，该边界演化无需反复进行网格变形或重建。此外，不断生长的沉积层通过一个辅助计算网格被处理为热导固体，在该网格上同时求解流体对流与固体导热。所计算的流体-沉积层界面温度随后被传回主网格，并用于依赖于温度的黏附模型。因此，该框架在物种输运、局部沉积、物理边界生长、流体动力学和传热之间建立了一个闭环反馈回路。贡献声明也已修改，以使该区分更加明确：

1. 开发了一种欧拉空间物种输运-浸入边界方法耦合方案，其中局部预测的沉积质量直接驱动物理流体沉积边界的空间变化生长。
2. 开发了一种辅助网格热力学方法，用于求解持续演化沉积层中的流体对流和热传导，并将所得界面温度反馈至黏附模型。
3. 将长期沉积生长表示为物理移动边界，而无需反复进行体网格重建。
4. 实现了虚拟时间策略，以弥合湍流计算与长期结垢生长之间巨大的时间尺度差异。
5. 在演化沉积表面上直接计算近壁流体动力学子量，并分析在不同速度、温度和沥青质浓度条件下竞争性的输运与黏附机制。

2. 控制方程

2.1. 连续相流体动力学

在本研究中，含有分散微米级沥青质颗粒的原油被视为连续相。其流体动力学由雷诺时均连续性方程和动量方程控制，分别如公式（1）和（2）所示。

（1）?(ρf)/?t + ?(ρf uj)/?xj = 0

（2）?(ρf ui)/?t + ?(ρf ui uj)/?xj = ??P/?xi + ?/?xj(μf(?ui/?xj + ?uj/?xi) ? ρf ui′uj′) + ρf gi

在公式（1）和（2）中，ρf 和 μf 分别表示连续相的密度和黏度。变量 u 和 P 分别对应连续相的时均速度和压力。

2.2. 湍流模型

在本研究中，湍流采用 Menter [23] 提出的 k–ω 剪切应力传输（k–ω SST）模型进行建模。该模型结合了 k–ω 公式在近壁区域的优势和 k–ε 公式在外流区域的优势，从而改进了边界层行为的预测能力，同时降低了对外流条件的敏感性 [23]。之所以选择 SST 模型，是因为准确计算近壁速度梯度和壁面剪切应力对于预测本问题中的颗粒黏附和沉积生长至关重要。湍流动能和比耗散率的相应输运方程如公式（3）和（4）所示。

（3）?(k)/?t + uj ?(k)/?xj = ?/?xj[(νf + σk νt) ?k/?xj] + Pt ? β* ω k

（4）?(ω)/?t + uj ?(ω)/?xj = ?/?xj[(νf + σω νt) ?ω/?xj] + γ/νt Pt ? β ω2 + 2(1?F1)(σω2/ω) ?(k)/?xj ?(ω)/?xj

混合函数 F1 作为主要切换函数，负责在 k–ε 和 k–ω 模型之间进行过渡。关于该公式的详细描述及相关模型常数可在 Menter [23] 的原始文献中找到。

2.3. 物种方程

为计算沉积质量，求解器需要一个额外的控制方程，因为该物理量无法直接从连续性方程和动量方程中提取。相反，其求解需要在流场中追踪特定相、物种或组分的质传输。因此，采用物种输运方程来计算沉积质量。该方程提供了一个综合框架，用于描述目标物种浓度在时间和空间上的演化过程。它包含非稳态项，用于描述时间变化；对流项，表示物种被流场速度场输运；以及扩散项，用于模拟物种在流体介质中的分子扩散或湍流扩散。

因此，为了准确捕捉质量输运现象并实现对整个计算域沉积质量的计算，将物种输运方程（公式（5）所示）纳入求解器。将该方程纳入控制方程体系后，可以同时将目标物种的空间和时间分布与流场进行求解。因此，可以以更高的物理保真度和数值精度来分析流体流动、质量传递和沉积机制的耦合作用。

（5）?(ρf Y)/?t + ?(ρf ui Y)/?xi = ?/?xi(ρf Deff ?Y/?xi)

在公式（5）中，变量 Y 代表颗粒质量分数，是物种输运公式中的主要输运标量。该变量描述了悬浮或沉积颗粒质量在流场中的空间和时间分布。参数 Deff 表示有效扩散系数，在颗粒质量传递建模中起着核心作用。该系数综合了不同扩散机制的效应，定义为分子扩散系数与湍流扩散系数之和，如公式（6）所示。在本研究的物种输运模型中，沉积形成颗粒或物种的输运采用同时包含分子扩散和湍流扩散贡献的有效扩散系数来表示。分子扩散系数描述了由随机分子运动引起的扩散，而湍流扩散系数则考虑了由湍流速度脉动引起的附加混合效应。这种处理方式使两种扩散机制都能纳入物种方程，从而改进了整个计算域中颗粒质量输运的预测精度。

（6）Deff = Dm + Dt = Dm + νt/Sc?

在公式（6）中，Sc? 表示湍流施密特数，定义为湍流动量扩散率与湍流质量扩散率之比。对于在物种输运框架内建模的微米级颗粒，Sc? 通常假设为等于1 [13]。在本研究中，这一假设被认为是合理的，因为分散的沥青质颗粒足够小，其输运主要由载流相湍流混合主导。因此，采用 Sc? = 1 作为湍流扩散项的实用闭合，而无需引入额外的经验参数。

此外，Dm 表示颗粒的分子扩散系数，量化了由随机分子运动引起的质量输运。该系数的值取决于颗粒材质和直径，并受载流流体温度和黏度的强烈影响。一般而言，温度升高和流体黏度降低会导致分子扩散率增大。由于 Dm 在确定近壁区域颗粒输运速率中起着关键作用，特别是在边界层内温度和浓度梯度显著的区域，不同工况下壁温变化和流体黏度的影响已被明确考察并纳入模型 [24]。这种处理使颗粒质量传递行为以及壁面附近相关沉积或黏附机制的预测更加准确。

（7）Dm/T? = (T?/T??) × (μ??/μ?)

在公式（7）中，T? 表示壁温，μ? 表示在壁面温度下评估的流体黏度。下标"0"表示这些参数的相应参考值，通常用于归一化目的以及考虑相对于参考状态下温度相关的流体性质变化。公式（5）至（7）共同描述了颗粒向近壁区域的输运，并纳入了流场中非稳态输运、对流、分子扩散和湍流扩散的效应。

然而，颗粒输运至壁面附近并不必然意味着颗粒沉积。要实现沉积，颗粒不仅需要到达近壁区域，还必须成功黏附到表面上。这一区别在物理上具有重要意义，因为许多颗粒在接近壁面后可能返回主流，或者可能由于剪切力从表面脱落。因此，对颗粒向壁面的黏附进行适当建模对于真实描述沉积过程至关重要。

一般而言，颗粒向表面黏附的建模主要有两种方法。第一种方法是对所有颗粒-壁面相互作用机制进行详细表征，包括分子间力、表面粗糙度、静电效应、颗粒形状和尺寸以及局部流体动力学力。尽管该方法在物理上较为全面，但其计算代价非常高 [18]。此外，许多所需输入参数，如颗粒的实际形状、粒径分布和详细的表面特征，在实际应用中往往未知或难以准确测量。

第二种方法为本研究所采用，其基于对物种输运方程的壁面边界条件进行修改。该方法不显式解析复杂的颗粒-表面相互作用，而是通过施加于物种输运方程的有效边界条件来体现颗粒黏附效应。如文献 [9] 所示，颗粒向表面黏附可通过一阶表面反应进行建模，使得壁面处的物种消耗速率等于颗粒沉积速率，如公式（8）所示。该方法在保留问题本质物理基础的同时，显著降低了计算复杂度。

此外，Johansen（2008）报告指出，基于物种的方法对于以扩散控制为主的颗粒输运是适当且高效的，并且在此条件下不需要求解更复杂的多相流方程。该结论为采用物种输运模型结合表面反应边界条件提供了有力的理论支持，从而进一步验证了本研究所采用建模策略的有效性。

（8）?a = ka cs

在公式（8）中，?a 表示单位面积上颗粒黏附到表面的质量速率，是表征颗粒在壁面黏附或沉积强度的直接指标。该参数表示单位时间内单位面积上黏附到表面的颗粒质量，在沉积现象的定量评估中起着关键作用。参数 ka 表示颗粒黏附速率常数，反映了颗粒黏附到表面的趋势，并将表面物理化学性质、颗粒特性和流动条件的影响综合为一个有效参数。

此外，cs 表示近壁区域的颗粒浓度，代表壁面紧邻处存在的物质量。在一阶表面反应模型中，该量作为颗粒黏附的主要驱动力。近壁颗粒浓度的增加直接导致黏附速率的提高，进而使沉积速率增大。cs 的值通过公式（9）计算，该公式规定了如何从物种输运场中提取壁面附近的有效颗粒浓度。

（9）cs = ρf Y(0)

在公式（9）中，Y(0) 表示壁面处的颗粒质量分数，指示由于输运和黏附过程而积累在表面上的颗粒物质量。作为边界量，该参数在描述物种输运场与壁面之间的相互作用中起着关键作用，并为计算近壁区域颗粒质量通量提供了依据。

假设近壁区域建立了质量平衡，则颗粒从主流输运至壁面的速率 ?t 必须等于颗粒黏附到表面的速率 ?a。该平衡假设意味着在稳态条件下，近壁控制体积中不会发生瞬时的颗粒质量积累或耗散，颗粒向壁面的输运直接导致颗粒黏附。在该框架内，颗粒黏附速率 ?a 被视为等于表面颗粒沉积速率 ?d。该等式反映了输运、黏附和沉积过程的连续性，并实现了流场中颗粒质量输运与壁面沉积物质量之间的直接联系。上述考虑在公式（10）中以数学形式表达，明确给出了近壁区域的质量平衡条件。

（10）?d = ?t = ?a

利用经典质量传递原理，颗粒输运速率 ?t 可根据公式（11）计算。

（11）?t = ?ρf Dm ?Y/?y (0) = kt (cb ? cs)

在公式（11）中，cb 表示主流中沥青质颗粒浓度，kt 表示颗粒向近壁区域输运的速率，或更一般地，表示质量传递系数。通过联立公式（8）至（11），可推导出一个适当的颗粒黏附建模边界条件，如公式（12）所示。

（12）?Dm ?Y/?y (0) = ka Y(0)

公式（12）提供了适用于壁面颗粒黏附建模的边界条件。可以证明，当黏附速率常数 ka 取非常大的值时，公式（12）退化为气体流动颗粒沉积的经典边界条件。在此极限情况下，黏附过程足够快，使得向壁面的质量输运成为控制沉积的速率限制步骤。

需要注意的是，k?的量级强烈依赖于壁面条件和流体性质，这些因素的变化会显著影响粒子的附着和沉积速率。以往的研究表明，粒子附着速率k?遵循阿伦尼乌斯型关系，并依赖于壁面温度T?，即升高表面温度可以增强克服附着能垒的概率。此外，壁面剪切应力的增加会降低附着速率，因为剪切力减少了粒子在表面的停留时间和附着力。因此，在式(13)中提出了粒子附着速率的半经验关系，该关系同时考虑了温度和壁面剪切应力的影响。

(13) k? = A(μ?/τw)exp(?E/RT?)

在式(13)中，E表示粒子附着到表面的活化能，通过实验确定。R为通用气体常数，μ?表示在表面温度下评估的流体黏度，τw为壁面剪切应力，A为经验确定的指数前因子。阿伦尼乌斯型形式源自Jamialahmadi等人[8]为沥青质沉积引入的粘附概率概念。在他们的公式中，粒子向壁面的输运和随后的稳定附着被作为不同的过程处理，附着概率的温度依赖性通过阿伦尼乌斯函数表示。他们的实验测量得到了约65.3 kJ·mol?1的表观附着活化能。因此，式(13)中的阿伦尼乌斯项应被理解为成功附着概率随温度变化的有效表征，而非基元化学反应的速率。流体动力学贡献遵循Maddahian等人[9]的公式，其中μ?/τw的比值为时间量纲，代表近壁特征停留时间。增加壁面剪切应力会缩短这一特征时间，从而抑制稳定附着。因此，式(13)为附着步骤提供了一个半经验但具有物理可解释性的闭合关系：沥青质向壁面的输运通过组分方程独立求解，而k?则控制到达近壁区域的粒子的后续稳定附着。

2.4 能量方程

基于式(13)，可以得出结论，表面温度是粒子沉积过程中的一个关键控制参数。因此，为了正确描述这一现象，必须将能量方程与其他控制流方程同时求解，以准确计算流场内部及表面附近区域的温度分布。然而，根据式(13)，沉积层的表面温度在沉积过程中起着尤为关键的作用。需要注意的是，这一表面并不代表固定的几何边界；而是由于连续粒子沉积和壁面上沉积层的逐渐生长，它会随时间动态演化。因此，沉积-流体界面的位置以及该界面上的温度都是时变的，并随沉积过程的变化而变化。

从传热角度来看，壁面上形成的沉积层起到额外的热阻作用，导致从管道外壁到沉积层表面再到主体流体的传热速率降低。随着沉积层厚度的增加，管道壁温度与直接接触流体的沉积层外表面温度之间会产生显著的温差。一方面，管道壁温度通常在数值模拟中作为边界条件或输入参数给定。另一方面，与流体接触的沉积层表面温度在改变流体热物性方面起着决定性作用，并直接影响粒子沉积和附着速率。换言之，随着沉积层随时间增长，控制流体-表面界面的热条件会发生变化，这些变化必须动态地纳入数值模型中。相应地，随着沉积层的逐渐生长，有必要同时求解固体区域中的传热方程。这一要求源于沉积层内的传热主要由热传导支配，而对流效应不参与固体沉积层内部的热量传递。流体区域和固体区域中传热的耦合求解能够更准确地预测温度分布、流体物性变化以及最终的粒子沉积速率。

该问题的主要数值挑战体现在以下几个方面。首先，浸入边界法(IBM)将流体-沉积界面视为内部边界条件，这从根本上改变了能量方程的系数矩阵。在界面处，既不能显式施加经典的狄利克雷边界条件，也不能施加诺伊曼边界条件；相反，表面温度必须作为解的一部分隐式确定。此外，该方法并不显式地在管道外壁上施加规定温度，同时也无法直接在固体区域内求解热传导方程。常规的共轭传热(CHT)问题的数值方法似乎是一种潜在解决方案；然而，在本案例中的主要困难在于，与标准共轭传热问题不同，固体区域（即沉积层）并非静止不动，而是随时间不断生长和改变几何形状。这一特征严重限制了经典共轭传热方法的直接适用性。

为解决这些挑战，本研究开发了一种原创的数值策略。在该方法中，引入了一套辅助计算网格来克服上述限制。该辅助网格有效地对应于未考虑浸入边界的主网格。具体而言，流体体积和速度场从主网格映射到辅助网格，随后在辅助网格上求解能量方程。求解能量方程后，温度场被传回主网格，用于评估每个计算单元中粒子沉积速率及其时间演化。该数值方法的完整细节，包括实现算法以及两个网格之间的数据交换过程，均在附录A中进行了全面记录。

在该框架下，如式(14)所示，对固体和流体区域均采用单一统一的能量方程，而非分别求解两个独立的能量方程。在固体区域中，速度设为零，这消除了能量方程中的对流项，仅保留非稳态项和热扩散项。这种表述方式使得同一控制方程能够一致地表示两个区域中的传热。(14) ?(ρc?T)/?t + ?(ρc?u?T)/?x? = ?/?x?(k_c·?T/?x?) + S

在式(14)中，k_c表示热导率，ρ代表密度，c?为比热容。在本研究中，式(14)用于捕捉沉积层内的热传导以及在流体区域求解完整的能量方程。相应的热物性采用类似于VOF（流体体积法）的混合方法定义，如式(15)–(17)所示[[25], [26], [27], [28]]。(15) ρ = ρf·γ + ρs(1?γ)(16) c? = c?f·γ + c?s(1?γ)(17) 1/k_c = γ/keff + (1?γ)/ks

在式(17)中，有效热导率考虑了分子和湍流热导率的综合贡献，由式(18)获得。需要注意的是，上标(s)和(f)分别指固体相和流体相。(18) keff = kf + kturb = kf + (c?f·μt)/Prt

在式(18)中，kf表示流体的分子热导率，kturb代表湍流热导率。参数Prt为湍流普朗特数，它将湍流动量扩散率与湍流热扩散率联系起来。本研究采用Prt = 0.85的常数值。Kays [29]表明，当湍流佩克莱数足够大时，在近壁湍流流动中Prt趋近于约0.85。这一值也已被用于经验证的数值传热计算中；例如，Lindqvist等人[30]在共轭传热CFD框架中与k–ω SST模型一起使用了Prt = 0.85。使用该常数值也与本问题的热力学和数值特征一致。在所研究的温度范围内，原油的分子普朗特数约为6.5–12.5，而最大近壁y?保持在约0.9左右。因此，近壁区域被直接解析，且在沉积层-流体界面处湍流黏度趋近于零，限制了在界面温度评估位置的湍流热扩散的直接贡献。一致地，Kügele等人[31]针对Pr = 9的流体表明，当近壁区域被充分解析时，Prt的变化对壁面与流体之间的温差仅产生可忽略的影响。然而，Prt并非通用常数。因为湍流热扩散率与Prt成反比，该参数的变化可以改变湍流热传递，从而影响预测的近壁温度。因此，所选值应被理解为针对当前工况的平均工程闭合而非普遍适用的常数。在强物性变化或浮力主导的条件下，如超临界流动，可能需要对Prt进行可变处理，正如Mohseni和Bazargan [32]所展示的那样。

更一般地，控制闭合关系表明主要模型参数具有不同的敏感性路径。增加表观活化能E会通过阿伦尼乌斯项降低附着系数k?，而指数前因子A则线性改变k?。增加Sc?会降低湍流粒子扩散率(νt/Sc?)，从而削弱向壁面的湍流质量输运。类似地，增加Prt会降低湍流热扩散率(νt/Prt)，从而通过预测的近壁温度和沉积层-流体界面温度间接影响沉积。这些关系阐明了模型对各参数的预期方向性响应。在本研究的模拟中，E、A、Sc?和Prt均固定在其选定值；因此，这些观察结果应被理解为解析敏感性路径，而非单独的全局不确定性分析。

3. 浸入边界法

一般来说，浸入边界法(IBM)，如图2所示，能够在无需体适网格或自适应网格生成的情况下求解涉及复杂几何形状和移动边界的问题。在该方法中，控制流方程在背景欧拉网格的流体单元内求解，而固体边界的位置和形状通过浸入边界表示独立定义[33]。因此，固体边界的运动和变形可以在无需重新划分网格的情况下追踪，与动态网格技术相比，这提供了显著的简便性和更低的计算成本。尽管如此，在每次更新固体边界位置后，需要重新评估单元中心变量以维持数值一致性[34]。

图2. 基于表面的浸入边界法中计算单元分类的示意说明[18]。

浸入边界法在处理复杂几何形状的问题方面具有重大优势，因为网格生成大大简化且网格质量对几何细节不那么敏感[35]。然而，其主要挑战在于边界条件的施加，因为流方程在背景网格上求解，且网格节点通常不与固体边界重合[34]。因此，浸入边界的处理会影响数值精度和守恒性。由于这些特征，浸入边界法已成为具有复杂或演化边界的模拟中对体适方法的有吸引力替代方案[[36], [37], [38], [39]]。

如图2所示，在本研究采用的截断单元方法中，计算单元被分为流体单元、固体单元和边界单元。固体单元通过赋予其零体积来从求解矩阵中排除，流体单元以常规方式离散化，边界单元的体积和表面积根据浸入边界的位置进行修正。所有几何修正和物性更新均由浸入边界库处理。在该框架下，不生成新单元，单元体积和表面积通过精确的几何运算计算。

截断单元方法的稳健性依赖于防止浸入边界与背景网格之间产生不完全几何相交。

由于基于STL的几何体常常与背景网格产生不规则的交叉，不精确的边或顶点交点可能导致形成开放的几何单元，从而引起数值不稳定性。因此，浸入边界在生成和处理过程中的几何质量控制对于确保数值解的准确性和鲁棒性起着决定性的作用。(19) $\sum_{c} S_f = 0$

(20) $S_f^{IB} + \sum_{c} \gamma_f S_f = 0$

为避免形成开放的几何单元，规则单元和交叉单元分别必须满足公式(19)和(20)所定义的条件[37]。在这些关系中，$\gamma_f$表示由单元切割过程获得的面积校正因子。修正后的浸入边界表面积$S_f^{IB}$根据公式(21)计算[37]。

(21) $S_f^{IB} = -\sum_{c} \gamma_f S_f$

在建立物种输运和浸入边界方法的控制方程后，本节展示了这两种方法的耦合策略以及沉积层生长过程在开源OpenFOAM框架中的实现。第一步，利用公式(8)在每个时间步和每个计算单元中计算颗粒沉积质量速率。该量表示单个时间步内从流场中脱离并参与近壁区域沉积的颗粒质量[18]。

在计算沉积质量速率之后，调用浸入边界库，并将边界单元中的沉积质量值传递到其中。完成传递后，通过插值邻域计算单元的平均沉积质量来更新浸入边界计算节点的位置，使其方式与公式(20)一致。该过程建立了欧拉物种输运场与浸入边界拉格朗日表示之间的耦合。为此，使用基于距离的搜索算法为每个边界单元构建邻域列表，以识别每个边界单元的四个最近邻边界单元。随后，根据平均沉积质量以及公式(22)移动浸入边界点和面的位置。空间平均的应用确保了沉积表面的变化被平滑引入，从而防止数值粗糙度和局部几何不连续性。

因此，在计算沉积质量后，根据公式(22)移动浸入边界点和面，同时应用缩放因子。换言之，每个时间步中边界的几何位移取为计算出的节点位移乘以该缩放因子，并同时定义一个对应的虚拟时间。引入该因子的目的是管理控制物理过程的特征时间尺度之间的差异[18]。使用缩放因子的原因是，求解具有温度依赖性质、移动边界和连续网格更新的湍流所需的时间步长极小，而沉积层的发生和生长发生在更大的物理时间尺度上，数量级为天。因此，在没有这种缩放方法的情况下，同时进行非常小的时间步长和非常长的物理时间范围的模拟将是计算上不可行的。

另一方面，缩放因子的使用不能歪曲问题的基础物理。因此，其值通过试错过程并结合数值预测与实验数据的比较来确定。结果表明，该因子取值为100时与实验观测值的吻合度最好。将该值增加到100以上会导致数值结果与实验结果之间出现显著偏差，从而降低模拟的物理保真度。然后，每个边界单元的质量根据公式(22)进行校正。参与结垢层生长计算的边界单元如图3所示。

图3. 结垢层生长的示意图。

如图3所示，每个节点（$\delta$）的位移可以使用公式(22)计算：

(22) $\delta = \frac{1}{2}\left(\frac{m_1}{A_1} + \frac{m_2}{A_2}\right) \frac{1}{\rho_p} \cdot (1-\varepsilon)$

在该方程中，$\delta$是垂直于沉积表面的节点位移，$m_1$和$m_2$是与连接到该节点的两个边界单元相关的校正沉积质量，$A_1$和$A_2$是它们对应的表面积。参数$\rho_p$表示沉积相的密度。系数$\varepsilon$表示沉积层孔隙率，使得$1-\varepsilon$表示固体体积分数，并在公式(22)中保留以维持沉积层质量到位移转换的一般形式。在本模拟中，未单独求解孔隙率演化模型。实验AFM和QCM-D研究表明，石油来源和含沥青质的沉积物可以形成几乎连续的、强吸附的、逐渐刚化的表面层[[40], [41], [42]]。因此，在本研究考虑的大陆尺度上，沉积相被近似为致密的固体层，$1-\varepsilon \approx 1$，因此在所有计算中使用$\varepsilon = 0$。该近似并不意味着微观空洞的完全缺失；而是未解析的孔隙空间和局部不均匀性在连续沉积相中被均质化。在此假设下，单位表面积上的校正沉积质量直接利用沉积密度转换为浸入边界的等效法向位移。

**4. 求解方法**

本研究采用的数值求解过程基于PIMPLE算法，该算法结合了PISO和SIMPLE方案，适用于具有强压力-速度耦合的非稳态问题。在该框架内，基本求解器pimpleDyMIBFoam被用作主要计算平台。在本改进的求解器中，局部颗粒沉积速率在每个时间步进行估算，然后利用浸入边界库更新浸入边界几何。然而，原始的pimpleDyMIBFoam求解器本身并不包含计算颗粒沉积的模型；此功能通过耦合的物种输运和附着公式在此引入。此外，该量的评估需要准确估算表面温度，因为颗粒附着速率和流体的所有热物性都直接依赖于温度。

因此，为了考虑热效应，在求解器中引入了第二求解网格用于求解能量方程。在第二网格上求解能量方程后，温度场被传回主网格，用于计算颗粒沉积速率以及更新浸入边界位置。此外，物种输运方程被纳入控制方程组中，其边界条件明确依赖于表面温度。通过这种公式化，沉积过程中固有的热-质量耦合被一致地捕获。

在求解物种输运方程之后，执行速度预测步骤并求解动量方程。随后，压力校正循环（p-loop）被激活，通过连续的压力校正更新速度场。在完成压力校正过程并使速度场收敛后，求解湍流方程，从而结束当前时间步。这一操作序列通过PIMPLE压力-速度耦合程序在每个时间步内执行，确保能量、物种输运、动量、压力和湍流方程的耦合求解，并使用更新的热物性。为了准确模拟沉积层的生长同时保持数值稳定性，采用可变时间步策略，使最大库朗数始终保持在0.25以下。该方法能够有效控制数值误差并防止因边界运动和流体性质变化引起的不稳定性。

为清晰起见，每个时间步内执行的浸入边界更新被明确总结。在第二网格上求解能量方程并将所得温度场映射到主网格后，求解物种方程并在浸入边界单元中评估局部沉积质量。这些值被传递到浸入边界模块中，在那里识别相邻边界单元并使用第3节描述的求和程序对沉积质量进行空间校正。校正后的沉积质量然后根据公式(22)转换为浸入边界节点的法向位移，并将规定的虚拟时间缩放应用于边界演化。随后移动浸入边界点和面，之后重新计算与背景网格的交叉并重新分类受影响的单元。流求解器所需的几何信息随后更新，并在更新后的构型上求解动量、压力校正和湍流方程。该浸入边界更新在每个物理时间步重复。

在本研究中，时间离散化使用二阶后向差分格式，为不稳定性模拟提供足够的精度。对于对流项，使用带有Van Leer限制器的TVD格式以抑制数值振荡并保持解的稳定性。扩散项使用二阶中心差分格式离散化，该方法在避免不必要的矩阵填充的同时保持高精度。对于压力-速度耦合，使用PIMPLE算法，设置nCorrectors = 2和nNonOrthogonalCorrectors = 2。动量方程使用PBICGStab求解器求解，而预测步中的压力方程使用PCG求解器结合GAMG预条件子求解。在所有模拟中，施加严格为$10^{-14}$的收敛准则，使具有强梯度的流动得到精确解析，同时确保数值稳定性和高解保真度。整体求解过程在图4的流程图中标示。

图4. 耦合的物种输运、能量、PIMPLE和浸入边界算法用于沉积层生长的计算流程。右侧分支详细说明了浸入边界几何更新程序。

**5. 结果与讨论**

**5.1. 问题描述**

为验证物种输运方程，研究了热交换器管道中沥青质颗粒的沉积，如图5所示。由于管道几何形状和施加的边界条件在周向是均匀的，该问题假定为在方位角方向上轴对称。因此，仅模拟管子的5度扇区，这保留了基本的沉积物理同时显著降低了计算成本。

图5. 用于在模拟沉积层生长过程中验证物种输运方程的几何形状[9]。

此外，为评估浸入边界方法的准确性，将该扇区的壁面建模为浸入边界。该方法允许首先验证沉积速率，然后可以在后续模拟中一致地模拟同一边界上沉积层的生长。

**5.2. 热物性和边界条件**

本模拟中使用的原油温度相关热物性总结在表1中。用于评估这些性质作为温度函数的相关性来自Jamialahmadi等人[8]报告实验测量的数据。这些关系被纳入数值模型中以准确考虑管道内沥青质沉积层生长模拟过程中流体性质的温度依赖性。

此外，本研究这部分中施加的边界条件列于表2。这些边界条件既用于验证算例，也用于管道内沥青质沉积层生长的数值模拟，确保验证框架与主要计算研究之间的一致性。

表2. 用于验证和管道内沥青质沉积层生长模拟的问题边界条件。

| 边界条件 | 速度(m/s) | 压力 | k | ω | 温度(K) | 物种(kg/m3) |
|---------|----------|------|---|---|---------|------------|
| 入口 | 0.35 | Zero Gradient | 0.0014 | 80.0358 | 0 | 3.5 |
| 出口 | - | Zero | Zero Gradient | Zero Gradient | Zero Gradient | Zero Gradient |
| 绝热壁I | No Slip | Zero Gradient | $1.0 \times 10^{-14}$ | $1.0 \times 10^{-14}$ | Zero Gradient | Zero Gradient |
| 加热壁 | No Slip | Zero Gradient | $1.0 \times 10^{-14}$ | $1.0 \times 10^{-14}$ | 393.0 | Robin |
| 绝热壁II | No Slip | Zero Gradient | $1.0 \times 10^{-14}$ | $1.0 \times 10^{-14}$ | Zero Gradient | Zero Gradient |

由于原油组成的固有多样性，本模型包含与附着和分子扩散闭式相关的三个经验参数：表观活化能$E$、前指数因子$A$和参考分子扩散系数$D_0$。在本研究中，这些参数并未针对当前验证数据独立重新校准。相反，采用Maddahian等人[9]报告的第一原油测试算例中的值$E = 85\ \text{kJ}\cdot\text{mol}^{-1}$、$A = 5.13 \times 10^8\ \text{m}\cdot\text{s}^{-2}$和$D_0 = 2.6 \times 10^{-10}\ \text{m}^2\cdot\text{s}^{-1}$，他们使用了相同的Arrhenius类型附着公式。

所采用的数值也与先前实验和建模研究[8,[43], [44], [45], [46]]中报道的活化能和沥青质扩散系数范围相一致。本次计算所使用的数值汇总于表3中。表3 沉积层生长模拟中物种输运方程所采用的模型常数。常数及数量如下：活化能（kJmol?1）为85，A（ms?2）为5.13×108，D0（m2s）为2.6×10?10。5.3. 网格无关性研究为评估数值解的网格无关性，考察了壁面上的沉积质量以及沿管道中心线的轴向速度分布。为此，针对所研究的问题，采用三种不同的计算网格——粗网格、中网格和细网格，其规格如表4所示。表4 用于评估管道中沥青质沉积的计算网格特征。网格名称及数量如下：粗网格200×25，中网格400×50，细网格800×100。本研究该部分所采用的计算网格示意图如图5所示。如图5所示，网格在径向方向上呈非均匀分布。管道核心区域采用相对较粗的网格，而在壁面附近则逐步加密网格，以充分解析近壁梯度。图6展示了管道中心线处的轴向速度分布。根据这些结果，可以得出结论：中分辨率网格能够为流场和物种场提供满意的表征，在计算成本与数值精度之间取得了适当的平衡，满足本研究的需要。此外，采用网格收敛指数（GCI）方法对中网格的网格不确定性进行了量化。所有监测变量的GCI值均低于1%，表明解具有网格无关性，证实了所选网格分辨率的充分性。下载：下载高分辨率图片（105KB）下载：下载全尺寸图片。图6 管道中沥青质沉积的网格无关性评估，对应管道中心线处的轴向速度分布。需要指出的是，y+的最大值约为0.9，这对于解析流动的粘性底层是足够的。本问题的一个关键特征在于，除了粘性底层之外，热边界层和浓度边界层也必须得到充分的解析。这一要求产生的原因在于，原油中沥青质颗粒的普朗特数和施密特数的数量级分别约为10和1000。这些较大的数值表明，原油中的热边界层和浓度边界层比粘性底层显著更薄。因此，为了准确表征壁面附近的热传递和质量传递，紧邻壁面的网格单元尺寸必须小于仅用于解析粘性底层所需的尺寸。5.4. 验证为验证整体控制方程组，模拟了含沥青质颗粒的原油在管道中的流动。为此，采用具有温度依赖性质的原油，如表1所列。问题的边界条件根据表2定义。图7展示了三个雷诺数——分别为5000、15000和25000——下壁面处沥青质颗粒的质量分数。可以看到，沉积颗粒量随着雷诺数的降低而增加；这一趋势的物理原因将在下一节中讨论。本研究的结果与Maddahian等人[9]报道的数值结果在三个雷诺数下均吻合良好。下载：下载高分辨率图片（118KB）下载：下载全尺寸图片。图7 本研究与Maddahian等人[9]数值结果在管道表面沥青质物种质量分数方面的对比。5.5. 沉积层的形成与增长利用上述描述的轴对称扇形模型，对图5所示管道几何中沥青质沉积层的增长进行了研究。为便于可视化，计算结果以三维形式展示在图8中。该图展示了在选定的轴向位置和模拟时间下沉积层厚度的时间演化。基于表2中汇总并在图5中说明的热、流体力学和物种边界条件，轴向沉积行为可划分为三个不同区域：在管道入口段，由于所施加的边界条件，未观察到沉积。这一行为表明，在该区域，热、物种和流体力学边界条件的综合效应不足以提供引发沉积过程所需的条件。在管道的加热中段，施加了壁面热流边界条件和一阶反应性物种边界条件（图5和表2），沥青质在壁面处持续被消耗，导致沿该段向下游方向沉积逐渐增强。如图8所示，在模拟时间范围内，沉积层在加热段出口附近的增长速率快于入口附近。然而，在模拟结束时，出口处沉积层厚度的增加仅约为管道半径的1.5%。这一微小的增加表明，在模拟时间范围内，沉积仍处于早期增长阶段，因此由于流通面积减小导致的壁面剪切应力增加尚不足以产生显著的沉积颗粒去除或脱落。在管道出口段，由于施加了零梯度边界条件，颗粒沉积速率逐渐降低，沉积层的增长速度也随之减缓。这一行为反映了传质驱动力的减小以及沉积过程在计算域末端逐渐受到限制。下载：下载高分辨率图片（136KB）下载：下载全尺寸图片。图8 管道不同轴向截面上颗粒沉积的空间分布及沉积层生长情况。5.6. 沉积层增长对连续相的影响以下研究沉积层增长对连续相温度的影响。该部分的边界条件根据表2定义。在热分析中，在管道加热段施加了0.2 mm的初始沉积层厚度。该值作为模拟开始时的初始结垢条件，以便从物理上有意义的非零厚度出发评估固体沉积层中的热传导。图9a展示了加热段出口处的径向温度分布。该分布的外侧部分对应于固体域内的传导传热，在初始阶段且沉积增长尚未开始时，由于建立了稳态流动条件，呈现出线性行为。如图9a所示，沉积层厚度的增加导致热阻增大，从而使得流体温度降低。这一行为产生于沉积层在从管壁到流体的传热路径中起到了额外的热屏障作用。下载：下载高分辨率图片（222KB）下载：下载全尺寸图片。图9 沉积层生长对管道内连续相温度的影响。此外，图9b展示了管道出口平均温度随时间的演化，显示出口温度随时间逐渐下降。在热预热系统中，出口温度的这种降低通常通过增加施加的热通量来补偿，而这又导致燃料和能源消耗的增加。因此，对沉积层热阻变化的详细分析具有重要的实际意义，将在以下讨论中予以阐述。沉积层的热阻是控制热预热器性能的最关键参数之一，它在决定流体与沉积表面之间能量交换速率方面起着决定性作用。因此，图10展示了加热管道入口段、中段和出口段处沉积层热阻与其初始热阻之比随时间的变化情况。下载：下载高分辨率图片（105KB）下载：下载全尺寸图片。图10 加热管道入口段、中段和出口段因沉积层厚度增加而产生的热阻。基于物种输运方程和能量方程所施加的边界条件，沉积在管道下游段更为显著，导致该区域沉积层的增长更大。结果表明，管道沿程的颗粒沉积并不遵循严格线性的空间分布，加热管道中段沉积增长的量级更接近出口处观测到的水平。在模拟结束时，加热管道出口处的沉积层厚度分别比中段和入口段大约36%和86%。鉴于本管道构型中的沉积过程仍处于线性增长阶段，可以清楚地观察到热阻随时间近乎线性增加。这一行为归因于模拟条件下观测到的有限沉积增长。由于沉积层厚度保持较小，由此导致的流通面积减小以及相应的壁面剪切应力增加也较为有限。因此，假设在模拟时间范围内颗粒去除或脱落可以忽略不计。采用物理壁面表示结合浸没边界方法的主要优势在于，它能够从解析的流场中直接计算物理上一致的壁面剪切应力。因此，所有与湍流相关的物理量，包括壁面剪切应力、近壁变量和湍流标度参数，均以自洽的方式计算，而不依赖于简化假设。这一能力使得对沉积层增长、流场变化以及壁面剪切应力变化之间耦合交互作用进行详细分析成为可能。在图11中，展示了加热壁面出口处模拟初始阶段和沉积层显著增长后最终阶段两者的u?和y?分布。对这些分布的检查揭示，尽管沉积层显著增长并伴随壁面的位移，u?–y?分布仍遵循近壁区域预期的物理行为。这一观察表明，本研究所采用的浸没边界方法即使在沉积增长引起移动边界存在的情况下，也能准确地再现边界层物理和壁面附近的湍流标度。下载：下载高分辨率图片（126KB）下载：下载全尺寸图片。图11 加热壁面在沉积增长开始前和沉积增长模拟结束时壁面标度速度分布（u?–y?）。换言之，在整个模拟过程中保持物理上一致的u?和y?分布，证明了所提出的方法不仅成功模拟了沉积层的增长，而且与湍流近壁流动的物理完全兼容。与传统的欧拉方法相比，这一特征代表了本研究方法的一个主要优势。5.7. 入口流速对沥青质沉积速率的影响图12展示了管道加热段内不同入口油速下的沥青质沉积速率，并与参考文献[8]中报道的实验数据进行了对比。与图7中观察到的情况类似，沉积速率沿管道长度并不保持空间均匀，而是表现出显著的轴向变化。因此，图12中报告的数值对应于整个加热段上的轴向平均沉积速率。对于图12中所示的流速效应研究，主体温度、壁面温度和主体沥青质浓度分别固定在Tb=358 K、Ts=393 K、cb=3.5 kg m?3。下载：下载高分辨率图片（159KB）下载：下载全尺寸图片。图12.

不同入口流速下沥青质沉积速率的对比：本框架与Han等[7]的既有数值结果、Seyyedbagheri和Mirzayi[47]的数值结果以及Jamialahmadi等[8]的实验数据之间的比较。结果表明，增加入口油流速会降低总体的沥青质沉积速率。本框架的预测值与实验测量值吻合良好，归一化均方根误差（NRMSE）约为5.1%。为进行同类比较，对Han等[7]和Seyyedbagheri和Mirzayi[47]已发表的数值预测结果进行了数字化处理，并使用与本框架相同的NRMSE指标进行评估，对应的NRMSE值分别约为6.3%和10.8%。因此，本框架相对于Han等的归一化预测误差降低了约19%，相对于Seyyedbagheri和Mirzayi[47]降低了约53%。尽管三种数值方法均再现了实验中观察到的沉积速率随流速增加而降低的趋势，但本框架与实验测量值在整体上的吻合度最为接近。从物理角度看，增加主体流动速率会增强近壁湍流输运，从而提高颗粒质量传递系数kt；然而，随之而来的壁面剪切应力的增加会通过剪切相关的附着关系降低附着系数ka。因此，与稳定表面附着相关的阻力相对于输运阻力有所增加。在本研究所考察的条件下，附着的降低超过了沥青质向壁面输运的增强作用，导致随着入口流速的增加，沉积速率净降低。这一解释与本模型所采用的耦合输运—附着框架一致：颗粒向近壁区域的输运由组分方程控制，而成功的表面附着则取决于局部附着系数。与实验趋势的吻合表明，该模型在研究流速范围内捕捉到了增强近壁输运和剪切抑制附着之间的竞争效应。5.8 壁面温度对沥青质沉积速率的影响图13展示了本研究中壁面温度对沥青质沉积速率的影响，并将数值预测结果与Jamiolahmadi[8]报道的实验数据进行了对比。在该分析中，所有其他操作条件——包括入口流速、主体流体温度和沥青质浓度——均保持不变，以便隔离壁面温度的影响。对于图13中展示的壁面温度效应研究，主体温度、入口流速和主体沥青质浓度分别固定为Tb=358K、V=0.92m·s?1、cb=3.5kg·m?3。结果清楚地表明，提高壁面温度会导致沥青质沉积速率显著增加。本框架的预测值与实验测量值吻合良好，NRMSE约为1.2%。为进行同类比较，对Han等[7]和Seyyedbagheri和Mirzayi[47]已发表的数值预测结果进行了数字化处理，并使用相同的NRMSE指标进行评估，对应的值分别约为11.8%和11.6%。因此，本框架相对于两种先前的数值方法，归一化预测误差降低了约90%。尽管三种模型均再现了实验中观察到的沉积速率随壁面温度增加而升高的趋势，但本框架的预测值在研究的温度范围内与实验测量值的吻合度始终显著更高。图13：不同壁面温度下沥青质沉积速率的对比——本框架与Han等[7]的既有数值结果、Seyyedbagheri和Mirzayi[47]的数值结果以及Jamialahmadi等[8]的实验数据之间的比较。从物理角度看，提高壁面温度同时影响输运和附着两个方面。首先，沥青质的分子扩散系数随温度升高而增大，从而降低近壁质量传递阻力并增强颗粒向表面的输运。其次，根据式(13)，附着系数对沉积物—流体界面温度表现出阿伦尼乌斯型依赖关系。提高表面温度会增加阿伦尼乌斯因子，从而增强颗粒稳定附着的概率和速率。在研究的温度范围内，由于附着项具有指数型温度依赖性，其贡献尤为重要。虽然提高温度也会促进分子向壁面的输运，但附着系数对温度的敏感性更强，导致附着阻力的降低更为显著。因此，输运阻力和附着阻力的共同降低导致了图13中观察到的沉积速率显著增加。沉积速率随壁面温度单调增加还表明，在当前操作条件下，任何可能阻碍颗粒向较热壁面输运的热泳效应并不足以克服增强分子扩散和温度促进附着的共同作用。由于热泳力并未作为独立的力项在本组分方程模型中表示，因此这一观察结果应理解为主导模型机制之间的比较，而非对热泳输运的独立定量评估。总体而言，与实验数据的吻合表明，该耦合输运—附着公式在研究的温度范围内捕捉到了沥青质沉积的壁面温度依赖性的主导特征。5.9 沥青质浓度对沉积速率的影响图14展示了沥青质颗粒浓度对沉积速率的影响，将数值结果与Jamiolahmadi等[8]报道的实验数据进行了对比。实验测量结果表明，增加颗粒浓度会导致沥青质沉积速率线性增加。这一线性趋势表明，在所考虑浓度范围内，沉积过程处于动力学控制区域，即沉积速率与可用于向壁面输运并随后附着到表面的颗粒量直接成正比。对于图14中展示的浓度效应研究，主体温度、入口流速和壁面温度分别固定为Tb=358K、V=1.25m·s?1、Ts=397K。然而，如图14所示，虽然数值预测和实验数据均表现出随颗粒浓度整体线性增加的趋势，但二者之间的偏差在高浓度下逐渐增大。该案例的归一化均方根误差（NRMSE）估计约为14.0%。这种逐渐增大的偏差表明，数值模型无法完全再现高浓度下实验观察到的物理行为。图14：本研究与Jamialahmadi等[8]实验数据之间，不同沥青质浓度下颗粒沉积速率的对比。该偏差的主要来源可归因于数值模拟中采用的单一恒定平均颗粒直径假设。相比之下，实验观察已证明沥青质颗粒具有多分散粒径分布，涵盖一系列颗粒直径而非单一特征尺寸[9]。高浓度下逐渐增大的偏差在很大程度上可归因于本模型中单一恒定颗粒直径的假设。实际上，沥青质颗粒呈多分散粒径分布，而颗粒尺寸影响决定沉积的两个子模型。首先，分子扩散系数和布朗扩散系数与颗粒直径成反比，因此较小颗粒比大颗粒更有效地输运至壁面。其次，附着步骤也具有尺寸依赖性，因为先前研究表明，活化能、黏附概率和附着速率在不同粒径组之间存在差异[9]。因此，用单一平均直径替代完整的粒径分布会消除细颗粒和粗颗粒对壁面输运和表面附着的差异化贡献。这种尺寸相关物理过程的缺失在高浓度下预计更为显著，因为在高浓度下聚集和颗粒尺寸重新分配效应更加强烈，这为数值结果与实验结果之间偏差逐渐增大的现象提供了合理解释。高浓度下逐渐增大的偏差在很大程度上可归因于本模型中单一恒定颗粒直径的假设。实际上，沥青质颗粒呈多分散粒径分布，而颗粒尺寸影响决定沉积的两个子模型。首先，分子扩散系数和布朗扩散系数与颗粒直径成反比，因此较小颗粒比大颗粒更有效地输运至壁面。其次，附着步骤也具有尺寸依赖性，因为先前研究表明，活化能、黏附概率和附着速率在不同粒径组之间存在差异[9]。因此，用单一平均直径替代完整的粒径分布会消除细颗粒和粗颗粒对壁面输运和表面附着的差异化贡献。这种尺寸相关物理过程的缺失在高浓度下预计更为显著，因为在高浓度下聚集和颗粒尺寸重新分配效应更加强烈，这为数值结果与实验结果之间偏差逐渐增大的现象提供了合理解释。因此，在数值模型中采用恒定颗粒直径引入了对输运和附着机制表示的简化。这种简化在高浓度下变得越来越重要，因为不同尺寸的颗粒对整体沉积过程的贡献更大，从而放大了数值预测与实验观测之间的偏差。尽管如此，数值数据和实验数据中均保持了线性趋势，表明所提出的模型能够捕捉到沥青质沉积随颗粒浓度变化的整体定性行为。因此，将粒径分布纳入建模框架预计能进一步提高定量精度，尤其是在高浓度条件下。5.10 局限性与适用范围本框架涉及若干简化假设，界定了其适用范围。沥青质相采用一个代表性颗粒直径表示，因此颗粒尺寸分布、聚集、破碎和尺寸选择性沉积并未被显式解析；这一局限性在高浓度下变得越来越重要。此外，附着关系是一种半经验闭合关系，其参数可能需要针对不同的原油、表面材料或操作范围进行重新评估。沉积物被处理为紧凑的连续层，因此内部孔隙率、演化的粗糙度和详细的表面化学未被解析。最后，轴对称RANS公式无法捕捉瞬态三维湍流或环向不均匀沉积；这些效应需要扩展到更精细的颗粒群和全三维模型中。6 结论本框架利用OpenFOAM开发了耦合数值方法，通过结合组分输运、温度依赖的表面附着模型、基于表面的浸入边界方法以及基于二次网格的热解策略，模拟了湍流管道流动中的沥青质沉积层生长。结果表明，将沉积层表示为物理移动边界对于近壁量（如壁面剪切应力和壁面标度变量）的一致预测以及捕捉沉积生长对流动和传热的反馈作用具有重要意义。与可用实验数据的对比表明，该模型使用先前建立的参数集和独立评估的虚拟时间缩放因子再现了主导沉积趋势。定量上，在流速和壁面温度案例中吻合良好，NRMSE值分别约为5.1%和1.2%，而在高浓度下偏差增大至约14.0%。这些结果表明，在所考察的条件下，颗粒附着是主导的速率控制机制：增加流速尽管增强了近壁输运仍会抑制沉积，而提高壁面温度主要通过降低附着阻力来加速沉积。本框架的定量性能还使用相同的NRMSE指标与先前数值方法进行了比较。对于流速相关案例，本框架的NRMSE约为5.1%，而Han等[7]约为6.3%，Seyyedbagheri和Mirzayi[47]约为10.8%，对应的归一化预测误差降低分别约为19%和53%。对于壁面温度相关案例，对应的NRMSE值分别约为1.2%、11.8%和11.6%，表明相对于两种先前方法，归一化预测误差降低了约90%。重要的是，这一改进的定量吻合是在额外解析沉积引起的边界演化、近壁流体力学反馈以及通过增长沉积层的导热过程中实现的。对于浓度相关案例，未声明等效的同类比较，本框架的NRMSE达到约14.0%。研究发现热场与沉积生长强耦合。

随着沉积层增厚，其热阻增加，从而降低流体温度并改变附着模型所需的沉积表面温度。这种耦合热效应通过基于次级网格的能量方程求解得以捕捉，使得生长中的固体沉积层内的热传导能够与移动界面一起被一致地处理。虚拟时间步策略提供了一种有效的方法，用于弥合湍流模拟所需的短数值时间尺度与沉积生长的长物理时间尺度之间的差距。在适当范围内，放大因子作为数值加速工具使用，不会改变主要的物理趋势；然而，过大的数值会导致与实验行为的明显偏差。从实际应用角度来看，所提出的框架除了提供总体沉积速率信息外，还能同时预测不断演变的沉积几何形状、壁面剪切应力、沉积表面温度以及相关热阻。这一能力使得在单一数值框架内评估操作参数对结垢形成和传热恶化的耦合效应成为可能。因此，该模型可以作为一种机理工程工具，用于比较操作条件、识别与增加结垢倾向相关的条件，以及在成本更高的实验或现场规模研究之前评估潜在的结垢缓解策略。其更广泛的意义在于将沉积动力学与不断演变的沉积所导致的流体动力学和热力学后果联系起来，而不是将结垢仅视为预设的质量传递过程。然而，向不同原油、表面材料和工业几何结构进行定量推广仍需额外的验证和适当的参数重新评估。总体而言，所提出的框架表明，沥青沉积不应简单地被视为质量传递问题，而应被视为具有不断演变几何形状的流体-固体-热耦合过程。当前模型的主要局限在于假设单一恒定粒径，这在较高浓度下变得更加受限。需要注意的是，当前模拟基于轴对称RANS近似处理湍流管道流动。因此，该模型未解析瞬态三维湍流结构、周向非均匀性或二次流动。在非轴对称加热、局部堵塞、强烈粗糙沉积和几何不均匀系统中，这些效应可能变得重要。因此，未来工作应通过群体平衡模型将该框架扩展至多分散颗粒、非牛顿原油流变学、瞬态操作条件以及使用RANS或LES基模拟的全三维工业几何结构。还应纳入更详细的去除机理和表面化学描述，以提高预测精度和工业适用性。

符号说明

所有量均以SI单位报告。基本罗马字母和希腊字母符号列于前；带下标形式、下标、上标、运算符和缩写在后分别定义。单位"1"表示无量纲量。

不带下标或上标的罗马字母符号

| 符号 | 定义 | SI单位 |
|------|------|--------|
| A | 附着速率关系中的Arrhenius前指数参数 | m s?2 |
| c | 颗粒质量浓度 | kg m?3 |
| D | 颗粒质量扩散系数 | m2 s?1 |
| E | 颗粒附着活化能 | J mol?1 |
| F | SST混合函数（用作F?） | - |
| g | 重力加速度 | m s?2 |
| k | 湍流动能 | m2 s?2 |
| k_c | 热导率 | W m?1 K?1 |
| m | 沉积质量 | kg |
| n | 沿外表面法线方向测量的坐标 | m |
| Nu | 努塞尔数 | - |
| P | 平均静压 | Pa |
| Pr | 普朗特数 | - |
| q | 传热速率 | W |
| R | 通用气体常数（8.314 J mol?1 K?1） | J mol?1 K?1 |
| Re | 雷诺数 | - |
| S | 体积热源项 | W m?3 |
| Sc | 施密特数 | - |
| t | 时间 | s |
| T | 温度 | K |
| u | 平均速度分量；公式(A1)中的参考平均速度 | m s?1 |
| U | 在初级和次级网格之间传递的速度矢量 | m s?1 |
| V | 主体或入口流速 | m s?1 |
| x | 笛卡尔空间坐标 | m |
| Y | 颗粒质量分数 | - |

不带下标或上标的希腊字母符号

| 符号 | 定义 | SI单位 |
|------|------|--------|
| α | 在初级和次级网格之间传递的固体体积分数 | - |
| β | SST特定耗散率模型系数 | - |
| γ | SST生产项系数；当下标为f时，浸没边界面面积修正因子 | - |
| δ | 浸没边界节点的法向位移 | m |
| ε | 沉积层孔隙率（此处紧凑层使用ε = 0） | - |
| μ | 动力粘度 | Pa s |
| ν | 运动粘度 | m2 s?1 |
| ρ | 密度 | kg m?3 |
| σ | SST湍流扩散系数 | - |
| τ | 剪切应力 | Pa |
| φ | 零法向梯度出口条件中使用的通用传输场 | 取决于场 |
| ω | 特定耗散率 | s?1 |

带下标和上标的变量

| 符号 | 定义 | SI单位 |
|------|------|--------|
| A?, A? | 与连接到浸没边界节点的两个边界单元相关的表面积 | m2 |
| c_b | 主体沥青颗粒浓度 | kg m?3 |
| c_s | 沉积表面的颗粒浓度 | kg m?3 |
| c_p | 混合的常数压力比热容 | J kg?1 K?1 |
| c_p,f | 流体的常数压力比热容 | J kg?1 K?1 |
| c_p,s | 固体沉积的常数压力比热容 | J kg?1 K?1 |
| D_eff | 有效颗粒扩散系数（分子加湍流） | m2 s?1 |
| D_m | 分子颗粒扩散系数 | m2 s?1 |
| D_t | 湍流颗粒扩散系数 | m2 s?1 |
| D?, D_m0 | 参考条件下的分子扩散系数 | m2 s?1 |
| F? | 初级SST混合函数 | - |
| k_a | 颗粒附着速率常数（附着速度） | m s?1 |
| k_eff | 有效流体热导率（分子加湍流） | W m?1 K?1 |
| k_f | 流体的分子热导率 | W m?1 K?1 |
| k_s | 固体沉积的热导率 | W m?1 K?1 |
| k_t | 颗粒传输或质量传递系数 | m s?1 |
| k_turb | 湍流热导率 | W m?1 K?1 |
| k_s/k_f | 固体与流体的热导率比 | - |
| m?, m? | 与连接到节点的两个边界单元相关的修正沉积质量 | kg |
| ?_a | 颗粒附着质量通量 | kg m?2 s?1 |
| ?_d | 颗粒沉积质量通量 | kg m?2 s?1 |
| ?_t | 向壁面的颗粒传输质量通量 | kg m?2 s?1 |
| P_t | 单位质量的湍流动能产生率 | m2 s?3 |
| Pr_t | 湍流普朗特数 | - |
| q″ | 施加的壁面热通量 | W m?2 |
| R_th/R_th,0 | 以初始值为归一化的沉积层热阻 | - |
| Sc_t | 湍流施密特数 | - |
| S_f | 计算单元的外向面面积矢量 | m2 |
| S_IB | 相交单元的修正浸没边界表面积矢量 | m2 |
| T_b | 主体流体温度 | K |
| T_s | 沉积表面（壁面）温度 | K |
| T_s0 | 参考表面温度 | K |
| u_i, u_j | 雷诺平均速度的笛卡尔分量 | m s?1 |
| u_i', u_j' | 瞬时速度波动 | m s?1 |
| u_inlet | 局部充分发展入口速度 | m s?1 |
| u? | 以壁面单位归一化的速度 | - |
| x_i, x_j | 笛卡尔空间坐标 | m |
| y? | 以粘性壁面单位归一化的壁面法向距离 | - |
| Y??? | 在壁面（y = 0）处评估的颗粒质量分数 | - |
| γ_f | 由单元切割操作产生的开放面积或面面积修正因子 | - |
| μ_f | 流体的动力粘度 | Pa s |
| μ_s | 在表面温度下评估的流体动力粘度 | Pa s |
| μ_s0 | 参考表面温度下的流体动力粘度 | Pa s |
| μ_t | 湍流动力粘度 | Pa s |
| ν_f | 流体的运动粘度 | m2 s?1 |
| ν_t | 湍流运动粘度 | m2 s?1 |
| ρ_f | 流体密度 | kg m?3 |
| ρ_s | 固体沉积密度 | kg m?3 |
| ρ_p | 沉积颗粒相的密度 | kg m?3 |
| σ_k | k方程的SST湍流扩散系数 | - |
| σ_ω | ω方程的SST湍流扩散系数 | - |
| σ_ω2 | 外分支SST交叉扩散系数 | - |
| τ_w | 壁面剪切应力 | Pa |
| β* | k方程中SST破坏项系数 | - |
| ?φ/?n | 通用传输场的法向导数 | 场单位 m?1 |

下标

| 下标 | 含义 |
|------|------|
| 0 | 参考或初始条件 |
| 1, 2 | 连接到浸没边界节点的两个边界单元 |
| a | 附着 |
| b | 主体流体 |
| c | 计算单元；也指对该单元所有面的求和 |
| d | 沉积 |
| eff | 有效（分子和湍流贡献的组合） |
| f | 流体相；对于S和γ，指计算单元的面 |
| i, j | 笛卡尔分量索引；重复索引表示求和 |
| IB | 浸没边界 |
| inlet | 入口边界 |
| k | 湍流动能方程 |
| m | 分子 |
| p | 常数压力属性或沉积颗粒相（取决于符号） |
| s | 表面或固体相（取决于符号） |
| t | 湍流贡献或颗粒传输（取决于符号） |
| w | 壁面 |
| ω | 特定耗散率方程 |
| ω2 | SST模型的自由流分支 |

上标和修饰符

| 标记 | 含义 |
|------|------|
| ′ | 瞬时波动分量 |
| 上划线（overbar） | 时间平均 |
| · | 随时间的变化率；在?中，表示单位面积的质量率 |
| ? | 以粘性壁面单位归一化的量 |
| * | 标准SST模型系数修饰符 |
| f | 流体相属性 |
| s | 固体相属性 |
| ″ | 单位面积量；q″表示热通量 |

运算符和符号

| 符号 | 含义 |
|------|------|
| ? | 偏导数 |
| Σ_c | 对计算单元c的所有面进行求和 |
| exp | 指数函数 |
| ??? | 在壁面（y = 0）处评估的值 |
| 粗体符号 | 矢量量 |

缩写

| 缩写 | 定义 |
|------|------|
| CFD | 计算流体力学 |
| CHT | 共轭传热 |
| GAMG | 几何代数多重网格 |
| GCI | 网格收敛指数 |
| IBM | 浸没边界法 |
| LES | 大涡模拟 |
| NRMSE | 归一化均方根误差 |
| PBICGStab | 预条件双共轭梯度稳定求解器 |
| PCG | 预条件共轭梯度求解器 |
| PIMPLE | 组合PISO-SIMPLE压力-速度耦合算法 |
| PISO | 压力隐式算子分裂法 |
| RANS | 雷诺平均Navier-Stokes |
| SIMPLE | 压力关联方程的半隐式方法 |
| SST | 剪切应力传输 |
| STL | 立体光刻表面格式 |
| TVD | 全变差缩减 |
| VOF | 流体体积 |

CRediT作者贡献声明

Sajad Khodadadi：原创草稿撰写、可视化、验证、软件、调查、概念化。Reza Maddahian：审稿与编辑、监督、资源、项目管理、数据管理。Xizhong Chen：审稿与编辑、监督、资源、项目管理、资金获取。

Sajad Khodadadi | Reza Maddahian | Xizhong Chen
上海交通大学化学化工学院化工系协同化学生物合成重点实验室，上海，200240，中国