单盘式浮顶油罐内含蜡原油温降规律数值计算研究

王敏，李敬法，张欣雨，宇波

1 油气管道输送安全国家工程实验室/城市油气输配技术北京市重点实验室/中国石油大学(北京)，北京 102249 2 中国石油化工集团国际石油勘探开发公司，北京 100029 3 机械工程学院，北京石油化工学院，北京 102617

单盘式浮顶油罐内含蜡原油温降规律数值计算研究

王敏1，李敬法1，张欣雨2，宇波3*

1 油气管道输送安全国家工程实验室/城市油气输配技术北京市重点实验室/中国石油大学(北京)，北京 102249 2 中国石油化工集团国际石油勘探开发公司，北京 100029 3 机械工程学院，北京石油化工学院，北京 102617

建立了考虑含蜡原油非牛顿性和析蜡相变过程的单盘式浮顶油罐内含蜡原油温降物理数学模型，并发展了一体化耦合求解方法。模型中，采用幂律方程描述罐内含蜡原油的非牛顿特性；采用湍流大涡模拟方法计算罐内含蜡原油的湍流流动过程；采用焓多孔介质理论计算罐内含蜡原油的析蜡相变过程，并跟踪罐内含蜡原油的相变界面。利用文献中的结果对本文模型进行验证。基于验证的程序，对实际单盘式浮顶油罐内含蜡原油温降过程进行计算，研究罐内含蜡原油温降过程的演化规律，并对罐顶、罐底凝油层的增长以及罐顶、罐底和罐壁处传热量的变化规律进行研究。研究结果可为单盘式浮顶油罐的设计和生产方案的制定提供一定的指导。

单盘式浮顶油罐；温降规律；非牛顿性；析蜡相变；大涡模拟

0 背景

中国所产原油大多为含蜡原油。含蜡原油具有易凝、高黏、流变性复杂等特点，在实际生产储存中通常需要加热处理。单盘式浮顶油罐是常用的含蜡原油储罐之一[1]。含蜡原油在罐内储存时会向外界环境放热而被冷却。当油温降低到含蜡原油析蜡点时，含蜡原油开始析蜡；当析蜡分数达到原油总质量的2%～3%时，含蜡原油将会胶凝[2]，严重影响实际生产，甚至造成生产事故。研究单盘式浮顶油罐内含蜡原油的温降过程，准确把握罐内含蜡原油的温降规律，对科学设计浮顶油罐的保温结构、合理制定浮顶油罐的周转方案具有重要意义。

目前，研究油罐内含蜡原油温降规律的方法主要有两种：实验研究[3-5]和数值模拟研究[6-9]，其中数值模拟因为投资小、研究周期短以及研究结果全面等优点成为近年来最主要的研究手段。在数值模拟研究方面，Lin等[6]采用数值计算方法对高径比为1/3～1，Pr数为1～1 000、Ra数为6×106～6×1010的油罐内温降规律进行研究。在此基础上，采用尺度分析方法拟合了无量纲温度与时间、储罐高径比、Pa数和Ra数的关系式。Oliveski等[7]采用实验方法和数值计算方法对Ra数分别为3.1×109和4.3×109的实验油罐内原油温降15 h的原油温度场和速度场进行计算，给出了不同温降时间下，油罐竖直中心线上的油温分布以及罐顶和罐壁处原油的速度矢量场分布。Oliveski[8]采用有限容积法对不同容积的实验罐内原油的温降过程展开研究，并采用量纲分析法回归了油罐的平均Nu数与储罐无量纲高径比、无量纲总传热系数以及Ra数和Pr数的换热关联式。李旺[9]采用大涡模拟方法对10万方双盘式浮顶油罐内原油温降过程进行计算，并分析了浮盘类型、保温层厚度、液位、油品物性、太阳辐射以及风速对罐内原油温降过程的影响。研究发现，浮盘类型、保温层厚度、液位和太阳辐射对油罐温降过程具有较大影响，而风速和油品物性对油罐温降过程影响较小。

虽然国内外学者已经对单盘式浮顶油罐内原油的温降规律开展了一定的研究，但这些研究通常将原油简化为牛顿流体，没有考虑含蜡原油的析蜡相变过程。而且研究中往往只给出油罐内某些特定位置处油温随温降时间的变化。就作者查阅的文献而言，目前还没有见到全面研究罐内含蜡原油温降演化规律以及罐内含蜡原油湍流流动与传热特性的研究成果。基于此，本文建立了考虑含蜡原油非牛顿性和析蜡相变过程的单盘式浮顶油罐内含蜡原油湍流温降物理数学模型，并发展了一体化耦合求解技术。通过对实际单盘式浮顶油罐内原油温降过程进行计算，研究罐内含蜡原油温度场及流场的演化规律，并对凝油层增长以及传热量变化进行分析。

1 物理数学模型及定解条件

1.1 物理模型

单盘式浮顶油罐包括：钢板层、保温层和含蜡原油，其传热过程包括液固(含蜡原油与钢板层)，固固(钢板层和保温层)耦合传热以及罐体与外界大气之间强制对流换热和油罐与罐底土壤层之间的导热。因此，单盘式浮顶油罐的传热过程非常复杂，计算难度较大。

对含蜡原油温降过程分析发现，初始时刻含蜡原油温度较高，含蜡原油以纯液态形式存在，且表现出牛顿流体特性；当油温开始降低到含蜡原油析蜡点Tw(含蜡原油中开始有蜡晶析出时的温度)以下时，含蜡原油中开始有蜡晶析出，固液分散体系开始形成。此时含蜡原油仍然表现出牛顿流体特性；当油温降低到含蜡原油反常点Ta(牛顿流体与非牛顿流体的转换温度)以下时，含蜡原油开始表现出非牛顿流体特性，此时含蜡原油仍以固液分散体系形式存在；当油温继续降低到含蜡原油显触点Tt(含蜡原油显现触变性的最高温度)以下时，含蜡原油中析出的蜡晶开始发生交联，形成蜡晶多孔介质结构[2]。表1给出了温降过程中含蜡原油的状态以及流变性的变化过程。

由于实际单盘式浮顶油罐的尺寸较大，罐内含蜡原油传热方式为湍流自然对流(实际1 000方单盘式浮顶油罐半径为6 m，当温差为30 ℃时，罐内含蜡原油Ra数高达2.0×1013)，因此要对这样的三维实际油罐进行数值计算，其计算难度和计算量都将非常大，计算非常耗时。考虑到单盘式浮顶油罐的对称性，本文将三维油罐简化为二维油罐，如图1所示。采用湍流大涡模拟方法对考虑含蜡原油的非牛顿性和析蜡相变过程的单盘式浮顶油罐内含蜡原油湍流温降过程进行计算，研究罐内含蜡原油的温降规律。

表1 温降过程中含蜡原油的状态和流变性变化Table 1 The state and rheological behavior of waxy crude oil in temperature drop process

图1 二维单盘式浮顶油罐物理模型Fig. 1 Sketch map of single-plate fl oating roof oil tank

1.2 数学模型

对浮顶油罐内含蜡原油温降过程分析发现，其本质是一个非稳态自然对流过程，其中涉及相变、固液相作用、非牛顿性、流固耦合以及湍流。结合1.1节建立的单盘式浮顶油罐温降物理模型，建立考虑含蜡原油非牛顿性[2,10-11]和析蜡相变过程[12-14]的单盘式浮顶油罐内含蜡原油湍流温降数学模型。考虑到物理模型中涉及到流固耦合，为了实现一体化耦合求解，将固体(钢板层和保温层)区域设定为黏度为无穷大的流体，其余参数按照实际参数取值。此外，本文采用大涡模拟方法[15]计算罐内含蜡原油的湍流过程。二维圆柱坐标系下单盘式浮顶油罐温降控制方程如式(1)～(4)所示：

其中，ρ、β和cp分别为密度、体积膨胀系数和定压热容。和分别为采用大涡模拟方法过滤后x方向和r方向的速度分量、压力和温度。sux和sur为蜡晶多孔介质对含蜡原油的束缚作用，本文称之为达西源项。sh为含蜡原油析蜡潜热项。含蜡原油的析蜡相变处理方法将在下文中详细介绍。μeff和λeff分别为有效黏度和有效导热系数，本文中将分别由式(5)和式(6)计算：

其中，μt和αt分别为亚格子涡黏系数和亚格子涡扩散系数。亚格子涡黏系数的构造是大涡模拟方法的关键，本文中将采用Smagorinshy-Lilly模型[16]构造，具体表达式如式(7)所示：

式中，Ls为亚格子混合长度，其定义式为Ls=min(κd, CsΔ)[16]。混合长度的计算中考虑了边界效应，对近壁区做了壁面修正；κ为冯卡门常数，本文中取0.41；d为计算位置距最近边界的距离；Cs为模型常数，本文中取0.1；Δ为过滤尺寸。

对于亚格子涡扩散系数，本文中采用式(8)计算。其中，Prt取0.87[16]。

1.2.1 非牛顿性

当油温降低到含蜡原油的反常点以下时，含蜡原油为假塑性流体，其剪应力与剪切率的关系可由幂律方程很好地描述。此外，当油温继续降低到含蜡原油显触点以下一定范围内时，幂律方程仍能很好的描述触变性含蜡原油的剪应力与剪切率的关系[2]，因此选择幂律方程描述含蜡原油的非牛顿性。式(9)给出了含蜡原油的表观黏度μa与剪切率˙之间的关系式。

其中，K为稠度系数，单位为Pa·sn；n为流动特性指数。K和n与油温T的关系由实验数据拟合得到。为罐内含蜡原油的剪切率，由式(10)计算[17]。

1.2.2 析蜡相变过程

当油温降低到含蜡原油析蜡点以下时，含蜡原油开始析蜡。此时，析出的蜡晶漂浮在含蜡原油中，随着含蜡原油一起流动。当油温降低到含蜡原油显触点以下时，析出的蜡晶开始相互交联，形成蜡晶多孔介质结构，并对包裹在其孔隙中的含蜡原油产生束缚。随着油温的继续降低，蜡晶多孔介质孔隙度逐渐减小，束缚作用逐渐增强。本文采用达西定律描述这种束缚作用，达西源项的表达式如式(11)和式(12)所示。

其中，K0为渗透率常数，本文中取10-7[14,18]。gs为析蜡分数，即某一温度下，含蜡原油中析出的蜡晶质量占原油总质量的比例，由DSC实验确定。

含蜡原油在析蜡相变过程中还将释放相变潜热，相变潜热表现为能量方程中的源项sh，如式(14)所示：

式中，ΔH为相变潜热。本文中ΔH=(1-gs)L，L为含蜡原油的析蜡相变总潜热。

1.3 边界条件和初始条件

1.3.1 边界条件

单盘式浮顶油罐罐顶、罐壁以及浮顶油罐右侧土壤表层均为第三类边界条件，其对流换热系数为25.63 W/(m2·℃)[2]。油罐中心轴为对称边界，土壤层右边界为绝热边界，下边界为恒温边界，温度取10 ℃。此外，在油罐边界条件设定中，包含了外界气温，本文以中国某地区11月到12月的实测气温为基础数据来计算。

1.3.2 初始条件

本文中规定初始时刻罐内含蜡原油、空气层、钢板层和保温层的温度均为30 ℃。地表温度为初始时刻的外界气温，为1 ℃。不同深度处的土壤温度由土壤恒温层和地表温度线性插值得到。

2 数值计算方法及模型验证

2.1 数值计算方法

本文采用有限容积法对控制方程进行离散，其中对流项采用GAMMA格式[19]离散，扩散项采用中心差分格式离散，时间项采用一阶全隐格式离散。采用基于非均分同位网格的SIMPLE算法耦合压力和速度，求解器选用强隐方法[20]。为提高计算速度，本文采用多重网格方法求解罐底土壤温度场，采用一体化求解方法对油罐部分(包括罐顶空气层和钢板层、罐底钢板层和罐壁保温层)温度场和流场等进行计算，再与罐底土壤温度场进行耦合。如此反复，直至两部分均达到收敛标准为止。本文中时间步长取5 s。此时收敛标准为：油罐部分质量不平衡量小于10-4，罐底土壤层导热方程余量小于10-6。

另一个关键问题是相变潜热ΔH的更新。在数值计算中，每个时层内相变潜热根据各时层内油温进行更新，相变潜热能否正确更新会直接影响计算过程的稳定性和收敛性。本文采用式(15)更新相变潜热[12]。

其中，α为亚松弛因子，本文中取0.2；m和m-1分别为当前迭代步和上一迭代步；为相变潜热的反函数[12]，即含蜡原油温度。根据前面的分析可知，当油温低于含蜡原油析蜡点时，ΔH=(1.0-gs)L，析蜡分数gs与油温之间的关系可以由DSC实验测得，由此，可以得到的计算式。

根据建立的单盘式浮顶油罐温降物理数学模型，本文采用Fortran语言编写了计算程序，相应的程序设计流程图如图2所示。

以实际1 000方的单盘式浮顶油罐为例，对其温降规律展开研究。其中1 000方单盘式浮顶油罐的尺寸x1、x2、x3、x4、r1、r2和r3分别为-5.0 m、0.01 m、8.01 m、8.02 m、6.0 m、6.11 m和11.11 m[1]。罐内含蜡原油、空气层、钢板层、土壤层和保温层的物性如表2。

考虑非牛顿性和析蜡相变过程的含蜡原油黏度由式(16)计算。

含蜡原油的表观黏度μa由式(9)计算得到。其中，稠度系数K和流动特性指数n由实验数据拟合得到，本文中的计算式见式(17)和(18)。

对含蜡原油进行DSC实验，拟合得到析蜡分数gs与油温T的关系见式(19)。

根据析蜡潜热ΔH、析蜡分数gs及油温的关系，可以推导出f-1(ΔH)的表达式如式(20)。

其中，含蜡原油的析蜡点Tw、反常点Ta和显触点Tt分别为27 ℃、25 ℃和23 ℃。

图2 程序设计流程图Fig. 2 Flow chart of the calculation program

2.2 大涡模拟方法验证

为了验证大涡模拟方法的准确性，本文在102×102的非均分网格上对直角坐标系方腔内Pr数和Ra数分别为0.71和1.58×109的湍流自然对流过程进行计算，将计算结果与文献中的实验结果[21]进行比较。图3给出了由本文大涡模拟方法计算的方腔水平中心线上无量纲温度Θ和无量纲速度Uy与文献结果的比较，从图中可以看出，本文计算结果与文献中的实验结果吻合良好。

2.3 非牛顿流体模型验证

为了验证非牛顿流体模型的准确性，在102×102的均分网格上，对Pr数为100，Ra数分别为104和105的方腔非牛顿流体稳态自然对流过程进行计算，并将计算得到的热边界平均Nu数与文献中的结果[10,11]进行比较，结果如表3所示。从表中可以看出，本文计算结果与文献结果吻合良好。

2.4 相变模型验证

为了验证相变模型的正确性，本文采用流固耦合方法对半径为31.9 mm，高为59.8 mm，壁面和底部包裹有厚度分别为2.85 mm的聚碳酸酯和5.97 mm聚丙烯酸固体保温材料的圆柱形区域内正二十烷的融化过程进行计算。计算中，正二十烷的初始温度为23 ℃，相变温度为36.4 ℃，顶部绝热，壁面和底部温度分别为45 ℃和32 ℃，物性等其他数据参见Jones等文献[13]。图4(a)给出了融化9601.6 s后区域内正二十烷的融化界面位置以及液态正二十烷中的温度场和速度矢量场分布。图4(b)给出了不同融化时间下，区域内正二十烷融化界面位置与实验结果的对比，从图中可以看出，本文计算结果与实验结果[13]吻合良好。

表2 含蜡原油、空气层、钢板层、土壤层和保温层的物性Table 2 Thermo-physical properties of different materials

3 单盘式浮顶油罐内含蜡原油温降规律研究

3.1 单盘式浮顶油罐温度场演化

图5给出了不同温降时刻，1000 m3单盘式浮顶油罐内含蜡原油温度场和速度矢量场分布。对图5进行分析可知，初始时刻，罐底部含蜡原油温度为30 ℃，土壤表层温度仅为1 ℃，热量快速从罐内含蜡原油传入罐底土壤层，罐底含蜡原油中形成了层状分布的温度场，如图5(a)。而且，在这一阶段，罐顶部含蜡原油将通过灌顶钢板层向外界大气传热，同时罐顶部含蜡原油损失的热量将通过浮升力，以自然对流方式从油罐中部获得补充，湍流自然对流过程开始形成。此时，流动过程仅限于油罐右上部，如图6(a)。随着温降过程的继续进行，罐内含蜡原油湍流流动区域逐渐扩展，温降100 h时，湍流自然对流区域逐渐扩展整个油罐区域，罐内x方向最大流速达到0.03 m/s，如图5(b)，罐底部层状分布的温度场消失，罐内油温近似均匀，且整体发生温降，如图5(b)～图5(c)。当罐内油温降低到含蜡原油显触点(23 ℃)以下时，罐内含蜡原油中析出的蜡晶开始相互交联，形成蜡晶多孔介质结构，并对包裹在其孔隙中的含蜡原油产生束缚，传热方式由固液分散体系自然对流转变为蜡晶多孔介质传热，类似于导热，此时，罐内含蜡原油流动过程明显减弱。温降300 h时，罐内x方向最大流速仅约为4.0×10-4m/s，如图5(d)。而且在罐顶、罐底和罐壁冷量的共同作用下，罐内出现了明显的冷热油界面，如图5(d)～图5(f)。

图3 大涡模拟方法验证Fig. 3 Validation of LES method

表3 Pr=100时，本文计算的方腔热边界平均Nu数与文献结果的对比Table 3 Comparison of Nu at the high temperature boundary with the numerical results obtained by this paper, Matin[10]and Turan[11]when Pr=100

图6给出了不同温降时刻，1 000方单盘式浮顶油罐内含蜡原油流函数的分布。值得指出的是，流函数的分布反映了罐内含蜡原油的流动情况。从图中可以看出，罐内流函数演化过程可以分为两个阶段：(1)温降时间小于200 h。在这一阶段，罐内核心区油温高于含蜡原油的显触点，核心区含蜡原油以纯液态或固液分散体系形式存在，流动性较好。而且随外界气温的剧烈波动，罐内流函数处于瞬态变化之中。此外，在外界大气冷量以及罐内湍流自然对流作用下，罐内含蜡原油的流动在整个油罐区域内进行。(2)温降300～600 h。在这一阶段，罐内核心区油温已经降低到含蜡原油显触点以下，蜡晶多孔介质结构已经形成，含蜡原油的流动性明显减弱，表现为流函数值仅约为温降200 h的1%。随着油温继续降低，析蜡量继续增大，蜡晶多孔介质孔隙度继续减小，流动性继续减弱。在外界大气冷量作用下，罐内含蜡原油的流动仅限于油罐右上部。

图4 正二十烷融化过程中不同时刻融化界面位置与实验数据对比Fig. 4 The comparison between calculated values and experimental results

当罐内含蜡原油还未完全胶凝时，单盘式浮顶油罐内油温近似均匀，且整体发生温降，本文称之为油罐核心区。本文通过考察罐内x=6.16 m，r=4.18 m处的油温变化来研究核心区含蜡原油温降情况，通过考察罐内平均油温来研究油罐整体温降情况，计算结果如图7所示。其中，平均油温＜T＞由式(21)计算。式中，为罐内第i、j控制容积的油温，ΔVij为第i、j控制容积的体积。

图7给出了单盘式浮顶油罐核心区油温和平均油温随温降时间的变化曲线。从图中可以看出，核心区温降曲线存在两个转折点：含蜡原油的析蜡点(27 ℃)和显触点(23 ℃)。具体分析可知，当油温降低到含蜡原油析蜡点以下时，含蜡原油中开始有蜡晶析出，并释放析蜡潜热。而析蜡潜热的释放会对罐内含蜡原油温降过程起到缓冲作用。当油温继续降低到含蜡原油显触点以下时，罐内已经析出的蜡晶开始相互交联，形成蜡晶多孔介质结构，造成罐内含蜡原油的传热方式从自然对流转变为蜡晶多孔介质传热，类似于导热，从而使得罐内含蜡原油的传热过程明显变慢。此外，当油温降低到含蜡原油显触点以下时，析蜡潜热的传递，以及油温降低后从核心区获得的热量补偿均变慢，在二者的共同作用下，罐内温降过程出现了阶梯状下降的特点，这与文献中水冷凝的实验结果[22]类似。

此外，从图7中还可以看出，当油温高于含蜡原油显触点时，罐内核心区油温与平均油温近似相等，而当油温降低到显触点以下时，核心区油温与平均油温的差异逐渐显现。这是因为，当罐内油温高于含蜡原油显触点时，罐内含蜡原油传热方式为湍流自然对流，在湍流自然对流过程的搅动下，罐内油温近似均匀，核心区油温与平均油温近似相等。当油温降低到含蜡原油显触点以下时，罐内含蜡原油已经胶凝，从罐顶、罐底和罐壁处传入的冷量开始冷却罐内含蜡原油，罐内出现了明显的温度梯度，如图5(d)～图5(f)，核心区油温将明显高于平均油温。

3.2 单盘式浮顶油罐罐底和罐顶凝油层的增长

当油温降低到含蜡原油显触点以下时，罐内含蜡原油中析出的蜡晶开始相互交联，逐渐形成蜡晶多孔介质结构。该结构一旦形成就会严重束缚含蜡原油，造成含蜡原油流动性显著变差。因此，含蜡原油的显触点可看成凝油层开始形成的临界温度。

图5 不同温降时刻，罐内含蜡原油温度场和速度矢量场分布Fig. 5 Oil temperature and velocity vector distributions in the tank at different moments

图8给出了单盘式浮顶油罐内r=5.52 m处罐顶和罐底凝油层随温降时间的增长。从图中可以看出，随着温降过程的进行，罐顶和罐底凝油层均表现出指数增长的特点。当罐内含蜡原油完全胶凝后，罐顶和罐底凝油层厚度达到最大，均为油罐液高。具体分析可知，初始时刻，罐底部将形成层状的油温分布，如图5(a)，凝油层出现。随着罐内湍流自然对流的形成，罐底凝油层将受到来自核心区含蜡原油的加热和冲刷，凝油层缓慢融化，温降45 h后，罐底凝油层消失。随着温降过程的继续进行，罐内油温逐渐降低，罐底和罐顶凝油层又相继出现，且快速增长。其中，罐底是罐内油温最低的区域，且传热方式以导热为主，因此罐底部凝油层出现相对较早，且增长缓慢。而罐顶部湍流自然对流过程比较剧烈，含蜡原油损失的热量能快速从油罐核心区得到补充。因此，罐顶部凝油层出现相对较晚，且随着罐内油温的降低，凝油层将快速增长，如图8所示。

3.3 单盘式浮顶油罐罐顶、罐底和罐壁处的传热特性

3.1节和3.2节分别研究了单盘式浮顶油罐内含蜡原油的温降规律和凝油层增长规律，本节主要研究单盘式浮顶油罐内含蜡原油与罐顶、罐底和罐壁钢板层之间的传热特性，其中，传热量由式(22)计算。

方程中，i代表计算位置编号；N代表计算边界网格总数；is代表传热面积；代表法向温度梯度。

图9给出了单盘式浮顶油罐罐顶、罐底和罐壁处传热量随温降时间的变化曲线。从图中可以看出，初始时刻，单盘式浮顶油罐内含蜡原油与外界环境之间的热交换主要发生在罐底部，其中，最大传热量达到约5.5 kW。随着罐底温度梯度的形成，罐底传热量快速减小。当温降时间超过220 h后，罐底部出现了凝油层，传热方式由自然对流转变为蜡晶多孔介质传热，类似于导热，传热量小于0.2 kW。对于罐顶部，当温降时间小于220 h时，罐顶部以湍流自然对流形式向外传热，且随着湍流强度的增大，传热量缓慢增大，最大传热量超过2.0 kW。随着温降过程的进行，罐顶部开始出现凝油层，传热方式转变为蜡晶多孔介质传热，传热量逐渐减小，最终稳定在0.8 kW左右。此外，由于罐壁处包裹有0.1 m厚的保温层，在整个温降过程中，罐壁处的传热量相对稳定，不超过0.3 kW。

图6 不同温降时刻，罐内含蜡原油流场分布Fig. 6 Stream function distribution in the tank at different moments

图7 罐内核心区油温和平均油温变化曲线Fig. 7 The core region temperature and average temperature in the tank

综上所述，对于单盘式浮顶油罐，在罐壁保温层作用下，罐壁处传热过程相对稳定，传热量较小。相对而言，罐顶和罐底传热量较大，且变化过程比较复杂。在短期储存时，罐顶和罐底部是主要的传热区域；而在长期储存时，罐顶是最主要的传热区域。因此，在单盘式浮顶油罐保温设计中，要重点关注罐顶和罐底部。

图8 罐顶和罐底凝油层的增长Fig. 8 Growth of gel oil at the bottom and top of the tank

4 结论

基于大涡模拟方法，本文建立了考虑含蜡原油非牛顿性和析蜡相变过程的单盘式浮顶油罐内含蜡原油湍流温降物理数学模型，并发展了一体化耦合求解技术。采用该模型，本文对实际1 000方单盘式浮顶油罐内含蜡原油温降过程进行计算，研究罐内含蜡原油温度场和流场的演化、凝油层增长以及传热规律，得到以下结论：

(1)当罐内含蜡原油还没完全胶凝时，在湍流自然对流作用下，罐内油温近似均匀，且整体发生温降。当罐内含蜡原油胶凝后，罐内传热过程为蜡晶多孔介质传热，类似于导热。此时，从罐顶、罐底和罐壁处传入油罐的冷量开始向油罐核心区推进，油罐中出现了明显的冷热油界面。

(2)就凝油层增长而言，对于本文算例，温降初期，罐底部出现了凝油层。然而，随着罐内湍流过程的形成，凝油层逐渐被加热和冲刷，最终完全融化。到温降后期，随着罐内油温的进一步降低，罐顶和罐底再次出现凝油层，且凝油层表现出指数的、波动式增长的特点。

(3)通过对罐内传热过程的研究发现，在罐壁保温层作用下，罐壁处传热过程相对稳定，传热量较小，而罐顶和罐底传热量相对较大，且变化过程更加复杂。就短期储存而言，罐顶和罐底是主要的传热区域；而就长期储存而言，罐顶是最主要的传热区域。因此，在单盘式浮顶油罐保温设计中，要重点关注罐顶和罐底部。

[1] 郭光臣, 董文兰, 张志廉. 油库设计与管理[M]. 东营： 中国石油大学出版社, 2006. [GUO G C, DONG W L, ZHANG Z L. The design and management of tank farm[M]. Dongying： China University of Petroleum Press, 2006.]

[2] 杨筱蘅. 输油管道设计与管理[M]. 东营： 中国石油大学出版社, 2006. [YANG X H. Oil pipeline design and management[M]. Dongying： China University of Petroleum Press, 2006.]

[3] 于达. 大型浮顶油罐测温系统的研发[J]. 油气储运, 2005, 24(8)： 41-43. [YU D. The development of the oil temperature measured system in large fl oating roof oil tank[J]. Oil and Gas Storage and Transportation, 2005, 24(8)： 41-43.]

[4] 赵志明. 大庆北油库浮顶储油罐非稳态传热问题的数值计算[D]. 大庆： 大庆石油学院, 2009. [ZHAO Z M. The numerical computation of transient heat transfer problems of the fl oating-roof tank in Daqing N. Tank farm[D]. Daqing： Daqing Petroleum Institute, 2009.]

[5] 李超, 刘人玮, 李旺等. 大型浮顶储罐原油温度场实验测试研究[J]. 工程热物理学报, 2013, 34(12)： 2 332-2 334. [LI C, LIU R W, LI W, et al. Crude oil temperature measurement in a large-scale fl oating roof tank[J]. Journal of Engineering Thermophysics, 2013, 34(12)： 2 332-2 334.]

[6] LIN W X, ARMFIELD S W. Long-term behavior of cooling fluid in a vertical cylinder[J]. International Journal of Heat and Mass Transfer, 2005, 48： 53-66.

[7] OLIVESKI R D C. Correlation for the cooling process of vertical storage tanks under natural convection for high Prandtl number[J]. International Journal of Heat and Mass Transfer, 2013, 57(1)： 292-298.

[8] OLIVESKI R D C, MACAGNAN M H, COPETTI J B, et al. Natural convection in a tank of oil： Experimental validation of a numerical code with prescribed boundary condition[J]. Experimental Thermal and Fluid Science, 2005, 29(6)： 671-680.

[9] 李旺. 大型浮顶油罐温度场数值模拟方法及规律研究[D]. 北京： 中国石油大学(北京), 2013. [LI W. Study on numerical simulation methods and regularities for temperature fi eld of large fl oating roof oil tank[D]. Beijing： China University of Petroleum (Beijing), 2013.]

[10] MATIN M H, POP I, KHANCHEZAR S. Natural convection of power-law fl uid between two-square eccentric duct annuli[J]. Journal of Non-Newtonian Fluid, 2013, 197： 11-23.

[11] TURAN O, SACHDEVA A, CHAKRABORTY N, et al. Laminar natural convection of power-law fl uids in a square enclosure with differentially heated side walls subjected to constant temperatures[J]. Journal of Non-Newtonian Fluid, 2011, 166： 1 049-1 063.

[12] VOLLER V R, PRAKASH C. A fi xed grid numerical modelling methodology for convection-diffusion mushy region phase-change problems[J]. International Journal of Heat and Mass Transfer. 1987, 30(8)： 1 709-1 719.

[13] JONES B J, SUN D, KRISHNAN S, et al. Experimental and numerical study of melting in a cylinder[J]. International Journal of Heat and Mass Transfer. 2006, 49(15-16)： 2 724-2 738.

[14] OMARI K E, KOUSKSOU T, GUER Y L. Impact of shape of container on natural convection and melting inside enclosures used for passive cooling of electronic devices[J]. Applied Thermal Engineering. 2011, 31(14-15)： 3 022-3 035.

[15] 张兆顺, 崔贵香, 徐春晓. 湍流大涡数值模拟的理论与应用[M]. 北京： 清华大学出版社, 2008. [ZHANG Z S, CUI G X, XU C X. Theory and application of Large Eddy Simulation[M]. Beijing： Tsinghua University Press, 2008.]

[16] ANSYS FLUENT Theory Guide[M]. ANSYS Inc., 2011.

[17] 袁世伟. 幂律非牛顿流体流动的数值计算与实验研究[D]. 上海： 华东理工大学, 2014. [YUAN S W. Numerical simulation andexperimental study of power-law fl uid[D]. Shanghai： East China University of Science and Technology, 2014.]

[18] WANG M, YU G J, ZHANG X Y, et al. Numerical investigation of melting of the waxy crude oil in an oil tank[J]. Applied Thermal Engineering. 2017, 115： 81-90.

[19] JASAK H, WELLER H G, GOSMAN A D. High resolution NVD differencing scheme for arbitrarily unstructured meshes[J]. International Journal for Numerical Methods in Fluids. 1998, 31(2)： 431-449.

[20] 陶文铨. 计算传热学的近代进展[M]. 北京： 科学出版社, 2001. [TAO W Q. Recent advances in computational heat transfer[M]. Beijing： Science Press, 2001.]

[21] AMPOFO F, KARAYIANNIS T G. Experimental benchmark data for turbulent natural convection in an air filled square cavity[J]. International Journal of Heat and Mass Transfer. 2003, 46(19)： 3 551-3 572.

[22] KOWALCZYK W, HARTMANN C, DELGADO A. Modeling and numerical simulation of convection driven high pressure induced phase change[J]. International Journal of Heat and Mass Transfer. 2004, 47(5)： 1 079-1 089.

Numerical study on the temperature drop characteristics of waxy crude oil in a single-plate fl oating roof oil tank

WANG Min1, LI Jingfa1, ZHANG Xinyu2, YU Bo3
1 National Engineering Laboratory for Pipeline Safety/Beijing Key Laboratory of Urban Oil and Gas Distribution Technology/ China University of Petroleum-Beijing, Beijing 102249, China 2 Sinopec International Petroleum Exploration and Production Corporation, Beijing, 100029 China 3 School of Mechanical Engineering, Beijing Institute of Petrochemical Technology, Beijing 102617, China

Physical and mathematical models for the temperature drop process of the waxy crude oil in the single-plate fl oating roof oil tank are established with the consideration of non-Newtonian behavior and wax precipitation, and a coupled integrative numerical procedure is developed. In this research, the non-Newtonian behavior is described by the power-law equation. The turbulent fl ow is calculated by Large Eddy Simulation (LES) method. Wax precipitation of the waxy crude oil in the tank is simulated by the enthalpy-porous media theory, in which the gel oil interface is tracked. The numerical program is validated by the results obtained from the literatures. Based on the verif i ed program, temperature drop process in a single-plate fl oating roof oil tank is calculated, and the evolution characteristics of the waxy crude oil in the tank are studied. Furthermore, gel oil growth on the tank top and tank bottom and heat fl ux on tank top, tank bottom and tank wall are also investigated in this research. The research results provide guidance for the scientif i c designing of the single-plate fl oating roof oil tank and making schedule of waxy crude oil turnover reasonably.

single-plate fl oating roof oil tank; temperature drop characteristics; non-newtonian behavior; phase change; LES

10.3969/j.issn.2096-1693.2017.02.025

(编辑 马桂霞)

王敏, 李敬法, 张欣雨, 宇波. 单盘式浮顶油罐内含蜡原油温降规律数值计算研究. 石油科学通报, 2017, 02： 267-278

WANG Min, LI Jingfa, ZHANG Xinyu, YU Bo. Numerical study on the temperature drop characteristics of waxy crude oil in a single-plate fl oating roof oil tank. Petroleum Science Bulletin, 2017, 02： 267-278. doi： 10.3969/j.issn.2096-1693.2017.02.025

*通信作者, yubobox@vip.163.com

2017-03-15

国家自然科学基金(No. 51325603)资助