盘管组倾角对浮顶油罐内含蜡原油融化过程的影响研究

王敏，邵倩倩，杨晓帆，李敬法

（1 北京计算科学研究中心，北京100193； 2 广东石油化工学院石油工程学院，广东茂名525000；3 北京师范大学地理科学学部，北京100875； 4 北京石油化工学院机械工程学院，北京102617）

引 言

浮顶油罐是最重要的储油设施。含蜡原油在储存过程中会因为向外界传热而冷却、析蜡和胶凝[1]。为保证罐内含蜡原油能够随时外输，通常需要定期、快速地对浮顶油罐进行加热。而盘管组是最常用的油罐加热融化方式之一。盘管组的布置方式将对罐内含蜡原油融化过程产生较大的影响[2-6]。因此，开展不同盘管组布置方式下，浮顶油罐内含蜡原油融化规律研究，对于科学设计盘管组的布置方式、合理制定浮顶油罐加热融化方案、削减成本等具有重要的意义。

盘管组加热的浮顶油罐内含蜡原油加热融化问题是当前油罐研究的热点之一，而且这方面研究目前已经取得了一定的成果[2-8]。针对油轮燃料舱重油加热问题，Magazinovic[2]提出了一种紧凑型盘管组布置方式。为了考察新方案的加热效果，建立了盘管组加热的燃料舱传热物理数学模型，采用数值计算方法对舱内重油加热温度场、流场以及盘管组对流传热系数进行计算，评估新型盘管组布置方式的加热效果。针对1.0×105m3浮动油罐，Sun 等[3]分别采用实验和数值模拟方法研究了盘管加热条件下，罐内含蜡原油的温升规律。此外，从盘管传入油罐的有效能利用角度出发，对盘管的管径和长度进行优选。在考虑外界气温、太阳辐射、盘管内蒸气压力以及原油物性的基础上，Sun 等[4]建立了盘管组加热的原油储罐传热物理数学模型，研究了盘管组结构、液高以及盘管组热流量对罐内原油温升过程的影响。Zhao 等[5]采用数值计算方法对盘管组加热的1.0×104m3浮顶油罐内原油区域流动与传热过程进行计算，对罐内油温随时间的变化进行分析，并对罐内纵向和横向油温随时间的演化规律进行研究。孙巍[6]建立了大型浮顶油罐内盘管加热的原油温升物理数学模型，研究了不同盘管排布方式、油品液位、加热蒸汽温度、外界环境等条件下，罐内原油温度场与流场的变化规律。陆雅红等[7]利用Fluent 软件对原油储罐蒸气盘管加热温度场进行模拟，以热流量为主要指标，对蒸气盘管的结构进行优化设计。刘凤荣[8]以大庆地区某2 万方浮顶油罐为研究对象，建立了原油储罐管式加热物理数学模型，研究了冬季浮顶油罐加热作用下罐内原油温度场和流场分布规律。Saifed-Din 等[9]分别采用实验和数值模拟方法研究了卧式储罐内水平圆管加热的罐内热水流场和温度场变化规律，分析了圆管数量对圆管传热系数、流场和温度场以及加热效率的影响。考虑到盘管组布置条件下，浮顶油罐内加热融化区域的复杂性，一些研究人员将罐底盘管组简化为一个恒热流的加热平面，以此来研究罐内介质的温升规律[10-12]。针对盘管加热的10×104m3浮顶油罐温升过程，刘佳等[10]将盘管组简化为水平加热面，通过调整水平面的热通量来模拟不同的盘管组热流量，基于该模型研究水平面加热下罐内原油的加热温度场、流场的演化规律，在此基础上，对不同液位、不同蒸汽流量对罐内原油温度场、流场的影响进行分析。以上这些研究主要是针对盘管组加热时罐内介质的温升过程展开，并未考虑介质可能发生的融化相变问题。

在盘管组加热的储罐介质融化研究方面，Wang等[13]采用层流模型对小尺寸模拟罐内单根水平圆管恒温加热的含蜡原油融化过程进行计算，研究了罐内含蜡原油的融化规律，分析了含蜡物性对水平圆管热流量及罐内含蜡原油融化液相分数的影响。相变储能也是储罐内材料融化研究的重要方面[14-17]。Ahmadi 等[14]建立了螺旋状加热盘管罐内石蜡融化储能物理数学模型，研究了盘管直径、螺旋管直径、盘管流量、进口温度等对罐内石蜡融化过程的影响。然而，以上研究均针对小尺寸油罐的层流融化展开，研究结果与实际浮顶油罐内含蜡原油湍流融化的结果存在较大差异。

综上所述，目前在盘管组加热的浮顶油罐内含蜡原油融化研究方面存在以下几方面的不足：（1）在实际储罐加热融化研究方面，现有的研究主要关注油罐温升过程，较少涉及罐内含蜡原油的融化相变过程；（2）仅有的一些油罐融化研究基本都是针对小尺寸油罐内层流融化开展的，该过程与实际浮顶油罐内含蜡原油湍流融化过程差异较大；（3）在盘管组结构对油罐加热效果的影响研究方面，现有的研究主要关注盘管组结构对温升速率、罐内原油温度场均匀程度的影响，并未从罐内流场、温度场以及盘管组热流量演化角度来分析盘管组布置方式对油罐加热和融化过程的影响。基于此，针对1000 m3双盘式浮顶油罐，在考虑含蜡原油融化过程中形态和物性变化的基础上，本文建立了盘管组加热的浮顶油罐内含蜡原油湍流融化物理数学模型。考虑到罐内含蜡原油区域的不规则性，引入浸入边界法，在结构化网格上处理浮顶油罐内各区域耦合传热问题，实现浮顶油罐的一体化耦合求解。利用该模型，对不同盘管组布置方式下罐内含蜡原油的融化规律进行研究，并对盘管组热流量及其分配情况进行分析。

1 物理数学模型

1.1 物理模型

浮顶油罐内传热过程包括罐内气-液-固耦合传热（空气层、含蜡原油与盘管组、钢板层、保温层之间的传热）、浮顶油罐与外界大气的强制对流换热以及浮顶油罐与罐底土壤层的导热，传热过程复杂。融化过程中，含蜡原油的形态和物性还将发生显著变化。具体而言，初始时刻罐内油温较低，蜡晶多孔介质形成。此时，含蜡原油表现出非牛顿流体特性。加热时，盘管附近蜡晶多孔介质温度快速升高，蜡晶多孔介质逐渐破坏，固液分散体系形成。此时含蜡原油仍表现出非牛顿流体特性。继续加热时，蜡晶开始融化，纯液态含蜡原油逐渐形成。具体过程见文献[13]。

图1 盘管组加热的浮顶油罐内含蜡原油融化物理模型Fig.1 Physical model of melting of waxy crude oil in a floating roof tank heated by coils

此外，浮顶油罐内含蜡原油流动处于湍流状态（对于1000 m3浮顶油罐，当特征温差为50℃时，Rayleigh 数高达3.0×1013），计算时间过长。考虑到油罐的轴对称特性，本文将三维油罐简化为二维油罐，如图1 所示。此外，罐底部布置有半径为r0盘管，壁面温度恒为T0，间距恒为δr，盘管组中心连线与r轴正方向夹角为θ，如图1(a)所示。

1.2 数学模型

盘管组加热的浮顶油罐内含蜡原油湍流融化和温升控制方程如式（1）～式（3）所示。其中，罐内湍流由大涡模拟方法计算[18-22]。

源项Sui和ST计算式为

其中，Kd为蜡晶多孔介质渗透率，由Kozeny-Carman 公式计算[13,22]。fi和q分别为由浸入边界法引起的分布力密度和分布热密度、ΔH为相变潜热，本文中ΔH=(1-gs)L，其中，L为含蜡原油相变总潜热。

1.2.1 融化相变过程 式（6）等号右边第一项代表了含蜡原油融化相变过程中蜡晶与含蜡原油之间的相互作用。-ui和-ui,s分别为含蜡原油和蜡晶的速度向量。本文计算中假设，当油温高于或等于含蜡原油显触点Tt时，ui,s等于ui；而当油温低于含蜡原油显触点时，ui,s为0[13]。本文中采用达西定律描述该束缚作用。此外，含蜡原油融化过程中还会吸收融化潜热，融化潜热表现为能量方程的内热源项ST。

1.2.2 非牛顿流体 当油温高于含蜡原油反常点Ta时，含蜡原油表现出牛顿流体的特性，其黏温曲线可以由实验数据拟合得到。当油温降低到含蜡原油反常点Ta以下时，含蜡原油表现出非牛顿性，本文中采用幂律方程描述含蜡原油的非牛顿性，如式（8）所示。此外，当油温降低到含蜡原油显触点Tt以下一定范围时，含蜡原油的表观黏度仍可以由幂律方程很好地描述[13]。因此，本文采用幂律方程描述含蜡原油的非牛顿特性。

1.2.3 浸入边界法 由图1 可知，罐内含蜡原油区域是不规则的，而空气层、保温层和钢板层都是规则的。如果采用非结构化网格离散含蜡原油区域，势必破坏地上浮顶油罐一体化求解过程，从而降低计算效率。基于此，本文引入浸入边界法，在结构化网格上对不规则的含蜡原油区域进行离散，在一套结构化网格上统一处理浮顶油罐内规则区域（空气层、钢板层、保温层）和非规则区域（含蜡原油区），从而确保了油罐部分的一体化求解过程。浸入边界法的核心思想是：将固体区域简化为固体边界，将其对流体的影响以分布力和分布热的方式作用于流体上。实际处理中，采用拉格朗日节点描述固体边界，采用欧拉节点描述流体区域，二者之间通过狄拉克函数传递数据[23-25]。浸入边界法的实施过程包括两步：预测步和校正步。其中，预测步为动量方程的求解[式（6）、式（7）中fi和q为0]。本节重点介绍校正步的实施过程。式（9）给出了校正步的控制方程

离散控制方程式（9），并将式（10）代入，整理得

其中，A、b和X的表达式分别为

式中，h为网格尺寸，UjB为拉格朗日节点j的速度。求解式（12），获得拉格朗日节点分布力FkB，并修正欧拉节点速度uj，如式（13）所示

如此就完成了速度场的校正。温度场校正与此类似，具体过程可见文献[23-25]。

1.3 边界条件和初始条件

1.3.1 边界条件 如图1 所示，浮顶油罐融化模型包括地上油罐部分和罐底土壤部分，其边界条件包括油罐边界条件和土壤边界条件。对于油罐，罐顶和罐壁为第三类边界条件。油罐中心轴为轴对称边界条件。罐底与土壤层左上部为耦合边界条件。对于土壤，土壤层左右边界分别为轴对称和绝热边界条件。下边界温度恒定为10℃。上边界包括左侧耦合边界和右侧第三类边界。盘管壁面为第一类边界条件。除此之外，本文边界条件中包含了外界气温。计算中取融化初始时刻气温为1℃，融化12 h 时气温为12℃，其余时刻气温由余弦函数插值计算。

1.3.2 初始条件 浮顶油罐内含蜡原油的融化初场是罐内含蜡原油储存温度场、速度场等。考虑到实际生产中不允许罐内含蜡原油完全胶凝，本文以浮顶油罐内含蜡原油储存一定时间后的罐内未完全融化的含蜡原油温度场、速度场等作为融化初场。

2 数值计算方法

基于有限容积法，本文采用Fortran 语言开发了盘管组加热的浮顶油罐内含蜡原油融化和温升数值计算程序，流程如图2 所示。程序中非稳态项采用全隐格式离散，时间步长取0.01 s。扩散项和对流项分别采用中心差分格式和GAMMA 格式离散。采用同位网格SIMPLE 算法耦合压力和速度[26-28]。本文的收敛标准为连续性方程余量小于1.0×10-6。本文中以1000 m3浮顶油罐为例，对罐内含蜡原油湍流融化规律进行研究。其中，1000 m3浮顶油罐的尺寸x1、x2、x3、x4、x5、r1、r2和r3分别为-5.0、0.01、8.01、8.02、8.52、6.0、6.11和11.11 m[25]。

本文对罐内三种盘管组布置方式下（表1），罐内含蜡原油的湍流融化过程进行计算，研究盘管组布置方式对罐内融化过程的影响，其中算例2 为对照算例。考虑到盘管组内加热循环水的顺利排放，工程设计中通常采用倾斜方式布置盘管组。

图2 本文计算流程Fig.2 Flow chart of this research

表1 三种盘管组倾角Table 1 Three inclinations of heating coils at tank bottom

罐内含蜡原油、空气、蜡晶、钢板、土壤和保温层的热物性如表2所示。

考虑融化过程中含蜡原油形态和流变性变化，含蜡原油黏度由式（14）计算。

表观黏度μa由式（8）计算得到。其中，稠度系数K和流动特性指数n通过实验数据拟合得到。本文中K和n的计算式为

表2 含蜡原油、空气、蜡晶、钢板、土壤和保温层热物性参数Table 2 Physical properties of waxy crude oil，air layer，wax crystal，steel layer，soil and insulating layer

析蜡量gs与油温-T的关系由DSC 实验数据拟合得到，如式（17）所示。

本文中以中国西部地区生产的含蜡原油为基础开展研究，该含蜡原油的析蜡点Tw、反常点Ta和显触点Tt分别为27、25和23℃。

3 模型验证

在计算之前，首先对本文模型进行验证。本节模型验证主要包括浮顶油罐传热模型验证、浸入边界法验证以及网格无关解和时间步长无关解的验证三部分内容。

需要说明的是，在模型验证方面，目前文献中还没有实际浮顶油罐内含蜡原油加热融化的实验数据。而在浮顶油罐内含蜡原油温降胶凝规律实验研究方面，目前已经有一些研究成果，且文献中也给出了相关的实验数据[29-30]。考虑到浮顶油罐内含蜡原油的加热融化过程与罐内含蜡原油的温降胶凝过程属于同一个过程的不同发展方向，因此，本节将采用浮顶油罐温降实验数据对本文传热模型进行验证。本研究中还采用了浸入边界法描述罐底盘管组与罐内含蜡原油的耦合传热过程，本节将利用文献中的数据对本文浸入边界法进行验证。

3.1 浮顶油罐传热模型验证

本节将利用文献中的实验数据对本文数值程序的准确性进行验证。验证中，本文选取中国西部地区和中国东部地区某油罐内量油孔（r=38.7 m）处实测纵向油温的平均值作为比较对象，从而确定模型的准确性。其中，油罐半径为40 m，罐内液高分别为17.7 m和19.44 m，其他数据见文献[29-30]。计算中，罐顶空气层、含蜡原油区、罐底土壤层和保温层的网格数（x×r）分别取97×241、160×241、82×302、266×13，总网格数为95504 个。时间步长取5.0 s。图3 比较了实测结果与本文计算结果，从图中可以看出，本文计算结果与实验结果吻合良好，最大误差控制在1℃以内。

3.2 浸入边界法验证

利用本文程序对方腔内圆管加热的空气自然对流进行计算，通过将计算得到的圆管表面稳态传热量与文献结果[23,31]进行对比来验证本文程序的正确性，结果如表3 所示。其中，r和l分别为方腔内圆管的半径和方腔特征尺寸，其他参数见文献[23,31]。从表中可以看出，本文结果与文献结果吻合良好。

表3 方腔内圆管恒壁温加热的稳定热流量Table 3 Steady heat flux of heating coil in a square cavity heated with constant temperature

3.3 网格无关解和时间步长无关解验证

图3 本文计算的量油孔纵向平均油温与实测结果的比较Fig.3 Comparison of simulated and measured longitudinal average oil temperature at gauge hatch

计算之前，需要确定本研究中采用的网格数和时间步长。本节将对三组网格和三组时间步长下，1000 m3双盘式浮顶油罐内5 根盘管组加热的含蜡原油融化过程（算例2）进行计算，研究网格数和时间步长变化对罐内含蜡原油融化过程的影响。计算之前，首先需要确定油罐融化初场。本节中，融化之前，罐内核心区油温约为24.55℃，罐底部约有0.485 m 厚的凝油层。图4 给出了三组计算网格和三组计算时间步长下，罐内含蜡原油平均温升随时间的变化曲线。从图4(a)可以看出，当网格总数分别为55312 个（含蜡原油区网格为254×178）和37224个（含蜡原油区网格为170×130）时，各时刻罐内平均油温近似相等，且二者与网格总数为28704个（含蜡原油区网格为142×110）的计算结果相差约3℃，因此，本文计算网格取37224 个（含蜡原油区网格为170×130）。同理，当时间步长分别为0.01 s 和0.005 s 时，罐内平均油温相差不大，且二者与时间步长为0.02 s 的计算结果差距较大，如图4(b)所示，因此本文计算中时间步长选取0.01 s。

图4 罐内平均油温随时间的变化曲线Fig.4 Variation of average oil temperature inside tank along with time

通过以上浮顶油罐传热模型和浸入边界法验证可以确定本文模型的准确性。通过网格无关解和时间步长无关解验证可以确认本文空间步长和时间步长选择的合理性。下文中将利用已经验证的模型以及筛选的空间步长和时间步长对三种盘管组布置条件下罐内含蜡原油的湍流融化规律进行研究。

4 结果及分析

本节主要研究不同盘管组布置方式对含蜡原油融化规律的影响，并对盘管组传入浮顶油罐的热量以及其分配情况进行分析。在研究之前，本文对融化初始温度场、速度场、析蜡情况等进行计算。

4.1 融化初场确定

计算之前，首先需要确定融化初场。根据1.3.2节融化初场的设定规则，以化工区储存一定时间的浮顶油罐温度场、速度场等作为融化初场，如图5所示。

图5 浮顶油罐融化初场设定Fig.5 Initial fields for melting process of floating roof tank

分析图5可知，融化开始之前，罐内的流动主要发生在油罐中上部和罐顶空气层，其中罐内含蜡原油和罐顶空气最大流速分别约为0.025 和0.015 m/s。融化开始后，在盘管组内通入温度恒定为80℃且快速流动的热水来加热融化罐内含蜡原油，罐内含蜡原油温度场和流场分布如图5(a)所示。在融化之前，罐内胶凝含蜡原油主要分布在罐底部，最大析蜡量超过9.5%。随着液高的增加，析蜡量逐渐减小。距罐底约0.4 m 范围内，蜡晶多孔介质已经形成，含蜡原油已经胶凝（根据本文选取的含蜡原油，析蜡量超过2.67%时，含蜡原油将会胶凝），如图5(b)所示。

4.2 管组加热的浮顶油罐内含蜡原油融化规律研究

图6～图8 分别给出了三组算例条件下1000 m3浮顶油罐内含蜡原油湍流融化温度场和速度矢量场的演化过程。分析可知，盘管组加热的浮顶油罐内含蜡原油融化过程可以分为两个阶段：融化阶段和温升阶段。

（1）0～300 s，罐内含蜡原油融化阶段。在这一阶段，盘管组传入浮顶油罐的热量主要用于融化油罐中上部悬浮的固态蜡晶和罐底部蜡晶多孔介质，并加热罐内含蜡原油。在融化阶段，罐内含蜡原油的湍流自然对流过程依次经历了形成、发展并最终稳定。具体分析可知，当融化时间为0～30 s，热量从盘管组传入罐内含蜡原油，加热融化盘管组附近的含蜡原油和蜡晶。与此同时，被加热的含蜡原油在浮升力作用下向上流动，湍流自然对流开始形成。当融化时间为30～150 s，热含蜡原油到达罐顶，并逐渐扩展到整个油罐区域。在这一阶段，罐内湍流自然对流区域逐渐扩散，强度逐渐增大，最大流速将达到0.7 m/s。此外，这一阶段，在湍流自然对流的加热作用下，不仅油罐中上部悬浮的固态蜡晶逐渐融化，而且罐底靠近罐中心的蜡晶多孔介质在热含蜡原油的冲刷下，逐渐解体并融化。当融化时间为150～300 s 时，罐内湍流自然对流过程逐渐稳定，形成了从上至下顺时针和逆时针的两个主旋涡，如图6(d)、(e)，图7(d)、(e)和图8(d)、(e)所示。而且，在湍流自然对流搅动下，罐内油温近似均匀，且温度整体升高。而且，在这一阶段，罐底靠近罐壁处的蜡晶多孔介质经过热油的冲刷，最终完全融化。至此，罐内含蜡原油融化过程完全结束。

（2）300 s 以后，罐内含蜡原油温升阶段。这一阶段，罐内含蜡原油已经完全融化，湍流自然对流已经扩展到整个油罐区域，且已经形成了上下两个主要的旋涡。此时，从盘管组传入浮顶油罐的热量主要用于加热罐内含蜡原油，并逐渐向罐顶空气层、罐底土壤层以及外界大气传递。而且，随着罐内油温的升高，盘管组与含蜡原油之间的温差逐渐减小，湍流自然对流强度逐渐减弱，表现为罐内含蜡原油最大流速降低到0.35 m/s。

值得指出的是，在浮顶油罐加热融化过程中，热含蜡原油不断冲击罐顶钢板层和罐底钢板层，导致罐顶空气层形成了一系列正逆相间的旋涡，如图6(e)、(f)，图7(e)、(f)和图8(e)、(f)所示。而对于罐底土壤层，由于土壤区域较大、热导率较小，而且加热时间较短（1200 s），传热方式为导热，因此罐底土壤层温度场未出现明显的变化。

图6 算例1条件下浮顶油罐内融化温度场和速度矢量场分布Fig.6 Temperature and velocity fields in tank during melting process of case 1

从图6(b)、图7(b)、图8(b)可以看出，在融化初期，当盘管组倾角θ逐渐增大时，被加热的含蜡原油上浮过程将倾向于盘管组正上方。分析其原因可知，盘管组左右两侧分别为罐中心和罐壁（图1），且罐中心体积明显小于罐壁处，因此，在上升热油的加热下，罐中心附近含蜡原油的温升速率将明显快于罐壁处。当盘管组非水平布置时（θ=-5.71°或5.71°），被盘管组加热的含蜡原油将沿盘管组向上流动。在以上两方面的作用下，当θ较小时，上升热油将向左偏移[图6(b)、图7(b)]，当θ较大时，热油垂直上升[图8(b)]。而且，当热油到达罐顶，并在罐顶堆积时，造成上升热油温差减小，阻力增大，因此在靠近罐壁的盘管上方将形成上升油流[图6(b)、图8(b)]，融化过程加速。

图9给出了罐内平均油温及融化液相分数随时间的变化曲线。分析可知，罐内平均油温和融化液相分数随时间持续增大，且增长率逐渐减小。结合前文的分析可知，在融化初期（0～70 s），算例1 条件下，融化区域主要位于罐中心处，而算例2 和算例3对应的融化区域位于盘管组上方，因此3 个算例表现出相似的融化规律，即三者液相分数和平均油温相差不大。随着融化过程的进行（70 s 以后），θ=5.71°时（算例3），靠近罐壁处的盘管组上方率先出现了上升油流[图8(c)]，罐内含蜡原油融化效率显著提高，且该过程一直持续到融化过程结束。对于温升过程而言，当融化时间超过300 s 后，罐内含蜡原油已经完全融化，湍流自然对流过程已经稳定，此时，算例1 和算例2 对应的罐内含蜡原油流动过程较为相似，尤其是靠近罐底部的逆时针主旋涡[图6(e)、图7(e)]，因此二者表现出相似的温升规律，且均小于算例3条件下的结果。

图7 算例2条件下浮顶油罐内融化温度场和速度矢量场分布Fig.7 Temperature and velocity fields in tank during melting process of case 2

综上所述，在本文计算条件下，与算例1和算例2相比，算例3表现出较高的融化效率和较快的温升速率。

4.3 盘管组热流量及其分配情况分析

盘管组传入浮顶油罐的热量是驱动罐内含蜡原油融化和温升的动力，本节对盘管组的热流量及其分配情况进行研究。式（18）给出了盘管组总热流量的计算式。

式中，qtot、qht、qtr和qpc分别为盘管组总热流量，从罐顶、罐底和罐壁钢板层传出的热流量，单位时间内含蜡原油因为温升和融化而吸收的热量，MW。盘管组总传热量Qtot是从融化开始到第n个时间步内盘管组传入油罐的热量，由式（19）计算，单位为MJ。

图10 给出了盘管组的总热流量及其分配。图10(a)为算例2 条件下，罐内盘管组热流量及其分配情况。从图中可以看出，时层热流量（qtot）经历了先快速减小，然后小幅度反弹，最后缓慢波动式减小的过程，这与浮顶油罐内层流融化类似[13]。分析其原因，融化初始时刻，盘管组与周围含蜡原油的温差高达55.8℃，热量以导热方式快速传入罐内，最大热流量达到100 MW。随着传热过程的进行，盘管组与周围含蜡原油的温差快速减小，造成热流量快速减小到约30 MW。随着融化过程的进行，罐内湍流自然对流过程逐渐形成，热量以对流的方式传递到整个油罐区域，热流量出现了波动式的小幅反弹，最大热流量约40 MW。随着罐内油温的升高，盘管组与含蜡原油的温差逐渐减小，自然对流强度逐渐减小，热流量缓慢波动式减小。就热流量分配而言，当融化时间在0～200 s 时，热流量主要用于加热和融化罐内含蜡原油。当融化时间超过200 s后，热流量基本全部用于加热含蜡原油。此外，从含蜡原油传入罐顶空气层、罐底土壤层和罐壁保温层的平均热流量约400 W。

图8 算例3条件下浮顶油罐内融化温度场和速度矢量场分布Fig.8 Temperature and velocity fields in tank during melting process of case 3

图9 浮顶油罐内含蜡原油融化特性Fig.9 Melting characteristic of waxy crude oil in floating roof tank

图10 盘管组热流量随时间的变化Fig.10 Variations of heat flux of coils along with time

如前文所述，盘管组布置方式将对罐内含蜡原油融化过程产生较大的影响。融化热流量表现出“山峰状”变化的特点，即先快速增大，然后波动式减小。这是因为自然对流形成后，热油上浮到达罐顶之前，融化区域快速扩大，热流量快速增大。当热油上浮到达罐顶之后，自然对流在罐内扩展，罐内流动过程变得非常复杂[图6(d)、图8(d)]，此时罐内中上部油温快速升高，最终导致热流量波动式减小。罐内湍流自然对流的发展对时层热流量和总热量具有决定性的影响。结合前文分析可知，由于融化初期和温升阶段罐内含蜡原油的流动过程类似，因此在这两个阶段时层热流量变化过程类似，而在融化中期，不同盘管组布置方式下，罐内含蜡原油自然对流过程非常复杂，时层热流量出现了一些不同。就总热量而言，由于融化过程中，罐内流动过程存在差异，导致不同盘管组布置方式下，盘管组总热量出现了差异，并在整个温升过程中维持不变（总热量差异约1000 MJ）。此外，从图10 中还可以看出，各时层盘管组的热流量处于剧烈波动之中，这是因为融化和温升过程中，罐内含蜡原油始终处于强湍流自然对流过程的缘故。

综上所述，在算例3对应的盘管组布置方式下，相同时刻盘管组总热量较大，传热效率较高。因此，为了提高浮顶油罐融化和加热效率，宜选用算例3盘管组布置方式。

5 结 论

利用在盘管组中通入热水来加热融化罐内含蜡原油是最常用的油罐加热方式之一，而盘管组布置方式对于罐内含蜡原油加热融化效率等具有重要的影响。基于此，本研究在考虑含蜡原油融化过程中形态和物性变化的基础上，建立了盘管组加热的浮顶油罐内含蜡原油融化物理数学模型。为提高计算效率，本文采用浸入边界法在结构化网格上对罐内规则区域和非规则区域进行统一处理，实现了油罐融化过程的一体化耦合求解。以实际1000 m3双盘式浮顶油罐为例，研究了不同盘管组布置方式下，罐内含蜡原油的融化和温升规律，分析了盘管组布置方式对浮顶油罐融化温升过程的影响，得到以下结论。

（1）盘管组加热的浮顶油罐内含蜡原油融化过程分为两个阶段：融化阶段和温升阶段。在融化阶段，罐内湍流自然对流依次经历形成、发展和稳定过程。而且，在形成期和发展期，热油在罐内上浮、扩展，流动过程相对复杂且混乱。此时，融化主要发生在油罐中上部和罐底靠罐中心处。在稳定期，罐内含蜡原油流动过程趋于稳定，盘管组传入的热量主要用于加热罐底靠近罐壁处的蜡晶多孔介质。在温升阶段，罐内形成了从上至下顺时针和逆时针的两个稳定的旋涡。而且，在这两个旋涡的搅动下，罐内油温近似均匀，且温度整体上升。

（2）时层热流量表现出先快速减小，然后小幅度反弹，最后缓慢波动式减小的特点。融化吸热流量随时间经历“山峰状”的变化过程，最大融化吸热流量达到约10 MW。三种盘管组布置方式下，盘管组总热量随时间的变化过程相似。在融化后期，算例3的总热量比算例1高约1000 MJ。

（3）本文算例中，算例3 表现出较快的温升速率、较大的传热量以及较高的融化效率。

符 号 说 明

cp——比定压热容，J/(kg·℃)

d——计算位置到最近边界的距离，m

——拉格朗日节点的分布力向量，kg/(m3·s2)

f——欧拉节点的分布力向量，kg/(m2·s2)

fi——由浸入边界法引起的分布力密度，kg/(m2·s2)

gi——重力加速度向量，m/s2

gs——析蜡分数

ΔH——相变潜热，J/kg

K——稠度系数

Kd——蜡晶多孔介质渗透率，m2

L——融化总潜热，J/kg

Ls——亚格子混合长度，m

n——流动特性指数

Prt——湍流Prandtl数

q——由浸入边界法引起的分布热密度，kg·℃/(m2·s)

qht,qpc,qtr,qtot——分别为向外界传出的热流量、融化吸热量、温升吸热量、盘管组总热流量，MW

qj——大涡模拟引起的热通量，kg·℃/(m2·s)

Ra——Rayleigh数

rp1,rp2,rp3,rp4,rp5——从左向右盘管中心径向坐标，m

r0——盘管半径，m

δr——相邻盘管间距，m

ST——能量方程源项，J/(m3·s)

Sui——动量方程源项，kg/(m2·s2)

Δsi——盘管表面表面第i段弧长，m

Ta,Tt,Tw——分别为含蜡原油反常点、显触点和析蜡点，℃

T0——盘管壁面温度，℃

t——时间，s

xp1,xp2,xp3,xp4,xp5——从左向右盘管中心轴向坐标，m

β——体积膨胀系数，℃-1

Δ——过滤尺寸，m

θ——盘管组所在平面与r轴正方向夹角，(°)

κ——冯卡门常数，1

λ——热导率，W/(m·℃)

μ,μt,μa——分别为动力黏度、亚格子黏度和表观黏度，Pa·s

ρ——密度，kg/m3

τij——剪应力，kg/(m·s2)