主题053：石油管道输送换热仿真

一、石油管道输送概述

1.1 石油管道输送的重要性

石油管道输送是现代能源运输体系的核心组成部分，承担着将原油从开采地输送至炼油厂、将成品油从炼油厂配送至终端用户的重要任务。全球范围内，管道输送占石油运输总量的70%以上，其安全性、经济性和可靠性直接影响国家能源安全。

管道输送过程中，石油的温度控制是确保输送安全和经济运行的关键因素。原油在低温环境下粘度急剧增加，流动性变差，可能导致管道堵塞；而温度过高则会造成能源浪费和安全隐患。因此，深入理解石油管道输送过程中的换热机理，建立准确的数值仿真模型，对于优化管道设计、制定运行策略具有重要意义。

1.2 石油管道输送的特点

石油管道输送具有以下显著特点：

（1）长距离输送

现代石油管道动辄数百甚至数千公里，如中俄原油管道全长近1000公里，中亚天然气管道更是超过1800公里。长距离输送意味着石油在管道中停留时间长，与外界环境的热交换充分，温降显著。

（2）变工况运行

管道输送量随市场需求波动，输量变化导致流速改变，进而影响对流换热系数和温降速率。此外，季节变化导致地温变化，也显著影响管道的热损失。

（3）多相流动

原油中含有蜡质、胶质、沥青质等成分，在低温下会析出蜡晶，形成固液两相流动。蜡沉积不仅减小管道流通面积，还增加热阻，加剧温降。

（4）复杂地质条件

管道穿越不同地质环境，包括冻土区、沙漠、沼泽、河流等，各段的土壤导热系数、地温差异显著，需要分段考虑热边界条件。

（5）高粘度特性

含蜡原油的粘度对温度极为敏感，在凝点以上粘度可能低至几十mPa·s，而在凝点以下粘度可飙升至数万mPa·s。这种强非线性特性给数值模拟带来挑战。

1.3 管道输送中的关键换热问题

石油管道输送涉及多个关键换热问题：

（1）热油输送的温降计算

热油输送是最常见的输送方式，通过加热站将原油加热至一定温度后输送。准确预测沿程温降，对于确定加热站间距、优化加热温度至关重要。

（2）保温层设计优化

保温层是降低热损失、减少加热能耗的关键。保温层材料选择、厚度优化需要综合考虑投资成本、运行能耗、施工难度等因素。

（3）蜡沉积预测与控制

蜡沉积是含蜡原油输送面临的主要问题。预测蜡沉积速率、确定清管周期、制定防蜡措施，需要建立蜡沉积与温度、流速、管壁剪切应力的耦合模型。

（4）停输再启动安全分析

管道因检修或事故停输后，原油温度逐渐降低，粘度增加。再启动时需要克服高粘度阻力，可能面临启动压力超限的风险。停输时间、再启动压力的计算是安全运行的重要保障。

（5）地温场影响分析

管道埋深、土壤导热系数、地下水流动等因素影响地温场分布，进而影响管道的热损失。特别是在冻土区，管道散热可能导致冻土融化，引发地质灾害。

二、石油管道输送的理论基础

2.1 管道流动基本方程

石油在管道中的流动遵循质量守恒、动量守恒和能量守恒三大基本定律。

（1）连续性方程

对于不可压缩流动，连续性方程简化为：

$$\nabla \cdot \vec{u}=0$$

在一维管道流动中，若管道截面积不变，则流速沿程恒定：

$$u=\frac{Q}{A}=\frac{4Q}{\pi D^2}$$

其中，$Q$为体积流量，$D$为管道内径，$A$为流通截面积。

（2）动量方程

一维管道流动的动量方程（伯努利方程考虑摩阻损失）为：

$$\frac{dp}{dz}=-\rho g\sin \theta -\frac{f\rho u^2}{2D}$$

其中，$p$为压力，$z$为沿程坐标，$\rho$为密度，$g$为重力加速度，$\theta$为管道倾角，$f$为达西摩阻系数。

摩阻系数$f$与雷诺数$Re$和管壁相对粗糙度$\epsilon /D$相关。对于层流（$Re<2300$）：

$$f=\frac{64}{Re}$$

对于湍流，采用Colebrook-White公式或显式近似公式如Haaland公式：

$$\frac{1}{\sqrt{f}}=-1.8{\log }_{10}\left[{\left(\frac{\epsilon /D}{3.7}\right)}^{1.11}+\frac{6.9}{Re}\right]$$

（3）能量方程

管道内流体的能量方程为：

$$\rho c_{p}u\frac{dT}{dz}=-\frac{4q_w}{D}+\frac{u}{J}\frac{dp}{dz}+\frac{f\rho u^3}{2DJ}$$

其中，$T$为流体温度，$c_p$为定压比热容，$q_w$为管壁热流密度（向外为正），$J$为热功当量。

对于输油管道，焦耳-汤姆逊效应（压力降导致的温度变化）和粘性耗散通常较小，能量方程可简化为：

$$\frac{dT}{dz}=-\frac{4q_w}{\rho c_{p}uD}$$

2.2 管道总传热系数

管道的热损失通过多层结构传递：管内对流、管壁导热、防腐层导热、保温层导热、土壤导热。总热阻为各环节热阻之和：

$$R_{total}=R_i+R_{steel}+R_{coating}+R_{insulation}+R_{soil}$$

各环节热阻计算如下：

（1）管内对流热阻

$$R_i=\frac{1}{\pi D_i{h}_i}$$

其中，$h_i$为管内对流换热系数，可用Dittus-Boelter公式计算：

$$Nu=\frac{h_i{D}_i}{k}=0.023R{e}^{0.8}P{r}^{0.3}$$

（2）钢管壁热阻

$$R_{steel}=\frac{\ln (D_o/D_i)}{2\pi k_{steel}}$$

钢管导热系数$k_{steel}\approx 45-50\text{ W/(m⋅K)}$，由于管壁薄，此项热阻通常很小。

（3）防腐层热阻

$$R_{coating}=\frac{\ln (D_{coating}/D_o)}{2\pi k_{coating}}$$

（4）保温层热阻

$$R_{insulation}=\frac{\ln (D_{ins}/D_{coating})}{2\pi k_{ins}}$$

保温层导热系数$k_{ins}$是温度的函数，对于常用保温材料：

• 聚氨酯泡沫：$k_{ins}=0.02-0.03\text{ W/(m⋅K)}$
• 岩棉：$k_{ins}=0.035-0.045\text{ W/(m⋅K)}$
• 硅酸钙：$k_{ins}=0.05-0.06\text{ W/(m⋅K)}$

（5）土壤热阻

土壤热阻是最复杂的一项，与埋深、土壤导热系数、地下水位等因素相关。对于埋地管道，采用以下公式：

$$R_{soil}=\frac{\text{acosh}(2H/D_{ins})}{2\pi k_{soil}}$$

其中，$H$为管道中心埋深，$k_{soil}$为土壤导热系数（通常$1.0-2.5\text{ W/(m⋅K)}$）。

总传热系数

基于管道外表面定义的总传热系数：

$$U_o=\frac{1}{\pi D_{ins}{R}_{total}}$$

或基于单位管长：

$$K=\frac{1}{R_{total}}[\text{W/(m⋅K)}]$$

2.3 热油输送温降计算

热油输送的温降遵循指数衰减规律。考虑稳态运行，能量方程简化为：

$$\rho c_{p}Q\frac{dT}{dz}=-K(T-T_g)$$

其中，$T_g$为管道中心埋深处的地温（视为常数）。

分离变量并积分，得到沿程温度分布：

$$T(z)=T_g+(T_0-T_g)\exp \left(-\frac{Kz}{\rho c_{p}Q}\right)$$

或写成：

$$T(z)=T_g+(T_0-T_g)\exp (-az)$$

其中，$a=K/(\rho c_{p}Q)$为温降系数，表征温降的快慢。

加热站间距计算

若要求管道终点温度不低于最低允许温度$T_{min}$，则加热站最大间距为：

$$L_{max}=\frac{\rho c_{p}Q}{K}\ln \left(\frac{T_0-T_g}{T_{min}-T_g}\right)$$

2.4 原油粘度-温度关系

原油粘度随温度变化显著，常用以下经验公式描述：

（1）Walther方程

$${\log }_{10}[{\log }_{10}(\nu +0.7)]=A-B{\log }_{10}T$$

其中，$\nu$为运动粘度（cSt），$T$为绝对温度（K），$A$、$B$为拟合参数。

（2）Arrhenius型方程

$$\mu ={\mu }_0\exp \left(\frac{E_{\mu }}{RT}\right)$$

其中，$\mu$为动力粘度，$E_{\mu }$为粘流活化能。

（3）Vogel-Fulcher-Tammann（VFT）方程

对于含蜡原油，在凝点附近粘度急剧上升，VFT方程能更好描述：

$$\mu =A\exp \left(\frac{B}{T-T_0}\right)$$

其中，$T_0$为理想玻璃化转变温度，通常低于凝点。

2.5 蜡沉积机理与模型

（1）蜡沉积机理

蜡沉积是分子扩散、剪切弥散和布朗运动共同作用的结果：

• 分子扩散：由于径向温度梯度，溶解的蜡分子从高温区向低温管壁扩散
• 剪切弥散：流速梯度导致蜡晶向管壁迁移
• 布朗运动：蜡晶的随机运动

对于输油管道，分子扩散是主要机理。

（2）蜡沉积速率模型

基于Fick定律，蜡沉积通量为：

$$N_w=-D_{wo}\frac{\partial C_w}{\partial r}{|}_{r=R}$$

其中，$D_{wo}$为蜡分子在油中的扩散系数，$C_w$为蜡浓度。

考虑管壁剪切应力的影响，采用修正模型：

$$\frac{d{m}_w}{dt}=\pi DL\left(k_d\frac{\partial C_w}{\partial r}-k_r{\tau }_w\right)$$

其中，$m_w$为沉积蜡质量，$k_d$为沉积系数，$k_r$为剪切剥离系数，${\tau }_w$为管壁剪切应力。

（3）蜡沉积厚度计算

蜡沉积导致管道内径减小，有效流通面积变化：

$$D_{eff}(t)=D_i-2{\delta }_w(t)$$

其中，${\delta }_w$为蜡沉积层厚度，与沉积量相关：

$${\delta }_w=\frac{m_w}{\pi D_{i}L{\rho }_w(1-{\epsilon }_w)}$$

${\rho }_w$为蜡密度，${\epsilon }_w$为沉积层孔隙率。

2.6 停输再启动分析

（1）停输温降

停输后，管内原油通过管壁向土壤散热，温度逐渐降低。停输温降可用瞬态热传导方程描述：

$$\rho c_p\frac{\partial T}{\partial t}=k\frac{{\partial }^{2}T}{\partial r^2}+\frac{k}{r}\frac{\partial T}{\partial r}$$

对于工程估算，可采用集总参数法：

$$T(t)=T_g+(T_s-T_g)\exp \left(-\frac{t}{{\tau }_{th}}\right)$$

其中，$T_s$为停输时油温，${\tau }_{th}$为热时间常数。

（2）再启动压力

停输后原油温度降低，粘度增加，再启动时需要克服的摩阻：

$$\Delta p=\frac{4{\tau }_{y}L}{D}+\frac{8\mu LQ}{\pi R^4}$$

其中，${\tau }_y$为屈服应力（对于非牛顿流体）。

对于宾汉塑性流体：

$$\Delta p=\frac{16{\tau }_{y}L}{3D}+\frac{32{\mu }_{p}Lu}{D^2}$$

其中，${\mu }_p$为塑性粘度。

三、热油输送系统仿真

3.1 热油输送系统概述

热油输送系统由加热站、输油泵、管道、保温层等组成。加热站将原油加热至设定温度（通常高于凝点10-20°C），输油泵提供输送所需压头，保温层减少热损失。

热油输送的设计目标是：在满足输送压力约束的前提下，最小化加热能耗和泵站能耗。这需要建立准确的温降模型和压降模型，优化加热站间距和加热温度。

3.2 热油输送数学模型

（1）稳态温降模型

对于稳态运行，沿程温度分布为：

$$T(z)=T_g+(T_0-T_g)\exp \left(-\frac{Kz}{\rho c_{p}Q}\right)$$

考虑地温随季节变化：

$$T_g(z,t)=T_{g,avg}+A_g\sin \left(\frac{2\pi t}{365}+\varphi \right)$$

其中，$T_{g,avg}$为年平均地温，$A_g$为地温年变幅，$\varphi$为相位角。

（2）压降计算

沿程压降由摩阻损失和位能变化组成：

$$\Delta p={\int }_0^L\left(\frac{f\rho u^2}{2D}+\rho g\sin \theta \right)dz$$

由于粘度随温度变化，摩阻系数沿程变化，需要分段计算：

$$\Delta p=\sum _{i=1}^N\left(\frac{f_i\rho u^2}{2D}+\rho g\sin {\theta }_i\right)\Delta z_i$$

（3）非稳态模型

考虑输量和地温随时间变化，建立非稳态能量方程：

$$\rho c_{p}A\frac{\partial T}{\partial t}+\rho c_{p}Q\frac{\partial T}{\partial z}=-K(T-T_g)$$

采用数值方法求解，常用迎风格式：

$$\frac{T_i^{n+1}-T_i^n}{\Delta t}+u\frac{T_i^n-T_{i-1}^n}{\Delta z}=-\frac{K}{\rho c_{p}A}(T_i^n-T_g)$$

3.3 Python仿真实现

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle, FancyBboxPatch
import matplotlib.patches as mpatches

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

# ==================== 热油输送基础参数 ====================
# 管道参数
D_i = 0.5          # 管道内径，m
D_o = 0.52         # 管道外径，m
delta_ins = 0.08   # 保温层厚度，m
D_ins = D_o + 2 * delta_ins  # 保温层外径，m
H = 1.5            # 埋深（管中心），m
L = 100000         # 管段长度，m

# 物性参数
rho_oil = 850      # 原油密度，kg/m³
cp_oil = 2100      # 原油比热容，J/(kg·K)
k_steel = 45       # 钢管导热系数，W/(m·K)
k_ins = 0.03       # 保温层导热系数，W/(m·K)
k_soil = 1.5       # 土壤导热系数，W/(m·K)

# 运行参数
Q = 0.5            # 体积流量，m³/s
T0 = 60            # 起点温度，°C
T_g = 5            # 地温，°C

# 计算流速
A = np.pi * D_i**2 / 4
u = Q / A

print("=" * 60)
print("热油输送基础参数")
print("=" * 60)
print(f"管道内径: {D_i*1000:.1f} mm")
print(f"保温层厚度: {delta_ins*1000:.1f} mm")
print(f"管段长度: {L/1000:.1f} km")
print(f"体积流量: {Q*3600:.1f} m³/h")
print(f"平均流速: {u:.2f} m/s")
print(f"起点温度: {T0:.1f} °C")
print(f"地温: {T_g:.1f} °C")

3.4 总传热系数计算

# 计算各环节热阻
# 1. 管内对流热阻（先假设一个对流换热系数）
h_i = 500  # W/(m²·K)，后续根据流动状态计算
R_i = 1 / (np.pi * D_i * h_i)

# 2. 管壁热阻
R_steel = np.log(D_o/D_i) / (2 * np.pi * k_steel)

# 3. 保温层热阻
R_ins = np.log(D_ins/D_o) / (2 * np.pi * k_ins)

# 4. 土壤热阻
import numpy as np
R_soil = np.arccosh(2*H/D_ins) / (2 * np.pi * k_soil)

# 总热阻和总传热系数（基于单位管长）
R_total = R_i + R_steel + R_ins + R_soil
K = 1 / R_total  # W/(m·K)

# 基于管道外表面的总传热系数
U_o = K / (np.pi * D_ins)

print("\n" + "=" * 60)
print("总传热系数计算")
print("=" * 60)
print(f"管内对流热阻: {R_i:.4f} (m·K)/W")
print(f"管壁热阻: {R_steel:.6f} (m·K)/W")
print(f"保温层热阻: {R_ins:.4f} (m·K)/W")
print(f"土壤热阻: {R_soil:.4f} (m·K)/W")
print(f"总热阻: {R_total:.4f} (m·K)/W")
print(f"总传热系数K: {K:.2f} W/(m·K)")
print(f"总传热系数Uo: {U_o:.4f} W/(m²·K)")

3.5 沿程温降计算

# 计算温降系数
a = K / (rho_oil * cp_oil * Q)

# 沿程温度分布
z = np.linspace(0, L, 500)
T = T_g + (T0 - T_g) * np.exp(-a * z)

# 计算加热站间距（假设最低允许温度为凝点+5°C）
T_pour = 30  # 凝点，°C
T_min = T_pour + 5  # 最低允许温度

if T_min > T_g:
    L_max = np.log((T0 - T_g)/(T_min - T_g)) / a
    n_stations = int(np.ceil(L / L_max))
else:
    L_max = L
    n_stations = 1

print("\n" + "=" * 60)
print("沿程温降分析")
print("=" * 60)
print(f"温降系数a: {a*1000:.6f} 1/m")
print(f"终点温度: {T[-1]:.2f} °C")
print(f"最大加热站间距: {L_max/1000:.2f} km")
print(f"所需加热站数量: {n_stations}")

# 可视化
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 温度分布
ax1 = axes[0]
ax1.plot(z/1000, T, 'b-', linewidth=2.5, label='Oil Temperature')
ax1.axhline(y=T_g, color='g', linestyle='--', linewidth=1.5, label='Ground Temperature')
ax1.axhline(y=T_pour, color='r', linestyle='--', linewidth=1.5, label='Pour Point')
ax1.axhline(y=T_min, color='orange', linestyle='--', linewidth=1.5, label='Min Allowable Temp')

# 标记加热站位置
for i in range(1, n_stations):
    z_station = i * L_max
    if z_station <= L:
        ax1.axvline(x=z_station/1000, color='gray', linestyle=':', alpha=0.7)
        ax1.plot(z_station/1000, T0, 'ro', markersize=8)

ax1.set_xlabel('Distance (km)', fontsize=12)
ax1.set_ylabel('Temperature (°C)', fontsize=12)
ax1.set_title('Temperature Distribution Along Pipeline', fontsize=13, fontweight='bold')
ax1.legend(loc='upper right')
ax1.grid(True, alpha=0.3)
ax1.set_xlim(0, L/1000)

# 温降速率
ax2 = axes[1]
dTdz = -a * (T0 - T_g) * np.exp(-a * z)
ax2.plot(z/1000, dTdz * 1000, 'r-', linewidth=2.5)
ax2.set_xlabel('Distance (km)', fontsize=12)
ax2.set_ylabel('Temperature Gradient (°C/km)', fontsize=12)
ax2.set_title('Temperature Drop Rate', fontsize=13, fontweight='bold')
ax2.grid(True, alpha=0.3)
ax2.set_xlim(0, L/1000)

plt.tight_layout()
plt.savefig('hot_oil_temperature_distribution.png', dpi=150, bbox_inches='tight')
plt.show()

四、保温层优化设计

4.1 保温层设计的重要性

保温层是热油管道降低热损失、减少加热能耗的关键。保温层设计需要平衡初期投资与长期运行成本：

• 保温层太薄：热损失大，加热能耗高
• 保温层太厚：材料成本高，施工难度大，投资回收期长

优化设计的目标是找到经济保温厚度，使得总成本（投资+运行）最小。

4.2 保温层优化数学模型

（1）热损失计算

单位管长热损失：

$$q=K(T-T_g)$$

年热损失费用：

$$C_{heat}=\frac{q\cdot L\cdot t_{year}\cdot c_{fuel}}{{\eta }_{heater}}$$

其中，$c_{fuel}$为燃料价格，${\eta }_{heater}$为加热炉效率。

（2）保温层投资

保温层投资包括材料费和施工费：

$$C_{ins}=L\cdot \pi D_{ins}\cdot {\delta }_{ins}\cdot c_{ins}$$

其中，$c_{ins}$为单位体积保温材料价格。

（3）总成本与优化

考虑资金时间价值，采用等年值法计算总成本：

$$C_{total}=C_{ins}\cdot \frac{i(1+i{)}^n}{(1+i{)}^n-1}+C_{heat}$$

其中，$i$为折现率，$n$为项目寿命。

优化问题为：

$$\underset{{\delta }_{ins}}{\min }{C}_{total}({\delta }_{ins})$$

4.3 Python优化实现

from scipy.optimize import minimize_scalar

# 经济参数
c_fuel = 0.1  # 燃料价格，$/kWh
c_ins = 500   # 保温材料价格，$/m³
eta_heater = 0.85  # 加热炉效率
i_rate = 0.08  # 折现率
n_year = 20    # 项目寿命年数

# 运行时间（假设全年运行）
t_year = 8760  # hours/year

# 资本回收系数
crf = i_rate * (1 + i_rate)**n_year / ((1 + i_rate)**n_year - 1)

def calculate_total_cost(delta_ins):
    """计算给定保温厚度下的总成本"""
    if delta_ins < 0.01:
        delta_ins = 0.01
    
    D_ins = D_o + 2 * delta_ins
    
    # 计算热阻
    R_ins = np.log(D_ins/D_o) / (2 * np.pi * k_ins)
    R_total = R_i + R_steel + R_ins + R_soil
    K = 1 / R_total
    
    # 计算热损失（取平均温差）
    T_avg = T_g + (T0 - T_g) * (1 - np.exp(-K*L/(rho_oil*cp_oil*Q))) / (K*L/(rho_oil*cp_oil*Q))
    q = K * (T_avg - T_g)  # W/m
    
    # 年加热费用
    Q_heat = q * L * t_year / 1000  # kWh/year
    C_heat = Q_heat * c_fuel / eta_heater
    
    # 保温层投资
    V_ins = np.pi * (D_ins**2 - D_o**2) / 4 * L
    C_ins = V_ins * c_ins
    
    # 总成本（等年值）
    C_total = C_ins * crf + C_heat
    
    return C_total, C_heat, C_ins, q

# 优化保温厚度
delta_range = np.linspace(0.02, 0.15, 100)
costs = [calculate_total_cost(d)[0] for d in delta_range]

# 找到最优厚度
opt_idx = np.argmin(costs)
delta_opt = delta_range[opt_idx]
C_total_opt, C_heat_opt, C_ins_opt, q_opt = calculate_total_cost(delta_opt)

print("\n" + "=" * 60)
print("保温层优化设计")
print("=" * 60)
print(f"最优保温厚度: {delta_opt*1000:.1f} mm")
print(f"年加热费用: ${C_heat_opt/1000:.2f}k")
print(f"保温层投资: ${C_ins_opt/1000:.2f}k")
print(f"年总成本: ${C_total_opt/1000:.2f}k")
print(f"单位热损失: {q_opt:.2f} W/m")

# 可视化
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 成本曲线
ax1 = axes[0]
C_heat_list = [calculate_total_cost(d)[1] for d in delta_range]
C_ins_list = [calculate_total_cost(d)[2] * crf for d in delta_range]

ax1.plot(delta_range*1000, np.array(C_heat_list)/1000, 'r-', linewidth=2.5, label='Heating Cost')
ax1.plot(delta_range*1000, np.array(C_ins_list)/1000, 'g-', linewidth=2.5, label='Insulation Cost')
ax1.plot(delta_range*1000, np.array(costs)/1000, 'b-', linewidth=2.5, label='Total Cost')
ax1.axvline(x=delta_opt*1000, color='k', linestyle='--', linewidth=1.5, alpha=0.7)
ax1.set_xlabel('Insulation Thickness (mm)', fontsize=12)
ax1.set_ylabel('Annual Cost ($k)', fontsize=12)
ax1.set_title('Economic Optimization of Insulation', fontsize=13, fontweight='bold')
ax1.legend()
ax1.grid(True, alpha=0.3)

# 热损失随保温厚度变化
ax2 = axes[1]
q_list = [calculate_total_cost(d)[3] for d in delta_range]
ax2.plot(delta_range*1000, q_list, 'b-', linewidth=2.5)
ax2.axvline(x=delta_opt*1000, color='r', linestyle='--', linewidth=1.5, label=f'Optimal: {delta_opt*1000:.1f}mm')
ax2.set_xlabel('Insulation Thickness (mm)', fontsize=12)
ax2.set_ylabel('Heat Loss (W/m)', fontsize=12)
ax2.set_title('Heat Loss vs Insulation Thickness', fontsize=13, fontweight='bold')
ax2.legend()
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('insulation_optimization.png', dpi=150, bbox_inches='tight')
plt.show()

五、蜡沉积预测与控制

5.1 蜡沉积问题概述

含蜡原油在温度降低至蜡析出温度（WAT, Wax Appearance Temperature）以下时，蜡晶开始析出。当管壁温度低于WAT时，蜡分子向管壁扩散并沉积，形成蜡层。

蜡沉积带来以下问题：

• 流通面积减小：蜡层占据管道空间，有效内径减小
• 摩阻增加：粗糙的蜡层表面增加流动阻力
• 热阻增加：蜡层导热系数低，增加热阻，加剧温降
• 清管成本：需要定期清管，增加运行成本

5.2 蜡沉积数学模型

（1）蜡析出温度（WAT）

WAT与原油组成相关，通常通过DSC（差示扫描量热法）测定。对于仿真，可假设WAT为已知参数。

（2）蜡扩散系数

蜡分子在油中的扩散系数可用Wilke-Chang公式估算：

$$D_{wo}=7.4\times 10^{-8}\frac{(\varphi M{)}^{0.5}T}{\mu V^{0.6}}$$

其中，$\varphi$为溶剂缔合因子（油取1），$M$为溶剂分子量，$V$为蜡分子摩尔体积。

（3）蜡沉积速率

基于Fick定律，考虑分子扩散和剪切剥离：

$$\frac{d{\delta }_w}{dt}=\frac{D_{wo}}{{\rho }_w(1-\epsilon )}\frac{C_{bulk}-C_{wall}}{{\delta }_{bl}}-k_s{\tau }_w$$

其中，${\delta }_{bl}$为边界层厚度，$k_s$为剪切剥离系数，${\tau }_w$为管壁剪切应力。

5.3 Python仿真实现

# 蜡沉积参数
WAT = 40  # 蜡析出温度，°C
C_wax_bulk = 0.15  # 原油中蜡含量，质量分数
rho_wax = 900  # 蜡密度，kg/m³
epsilon_wax = 0.3  # 蜡层孔隙率
D_wo = 1e-9  # 蜡扩散系数，m²/s
delta_bl = 0.001  # 边界层厚度，m
k_shear = 1e-8  # 剪切剥离系数

# 初始条件
D_eff = D_i  # 有效内径
wax_thickness = 0  # 蜡层厚度

# 模拟时间
t_days = 30
t_hours = np.linspace(0, t_days*24, 500)

# 存储结果
wax_thickness_history = []
D_eff_history = []

# 假设油温沿程分布（简化）
T_wall = T_g + 5  # 管壁温度，假设略高于地温

for t in t_hours:
    if T_wall < WAT:
        # 计算蜡溶解度（简化线性模型）
        C_wall = C_wax_bulk * max(0, (T_wall - T_g) / (WAT - T_g))
        
        # 计算沉积速率
        deposition_rate = D_wo * (C_wax_bulk - C_wall) / delta_bl / (rho_wax * (1 - epsilon_wax))
        
        # 计算剪切应力
        tau_w = 0.5 * f * rho_oil * u**2
        
        # 计算剥离速率
        removal_rate = k_shear * tau_w
        
        # 净沉积速率
        net_rate = max(0, deposition_rate - removal_rate)
        
        # 更新蜡层厚度
        wax_thickness += net_rate * (t_hours[1] - t_hours[0]) * 3600
    
    # 更新有效内径
    D_eff = D_i - 2 * wax_thickness
    
    wax_thickness_history.append(wax_thickness * 1000)  # mm
    D_eff_history.append(D_eff * 1000)  # mm

print("\n" + "=" * 60)
print("蜡沉积预测")
print("=" * 60)
print(f"运行时间: {t_days} days")
print(f"最终蜡层厚度: {wax_thickness_history[-1]:.2f} mm")
print(f"有效内径减小: {(D_i - D_eff)*1000:.2f} mm")
print(f"流通面积减小: {(1 - (D_eff/D_i)**2)*100:.2f}%")

# 可视化
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 蜡层厚度随时间变化
ax1 = axes[0]
ax1.plot(t_hours/24, wax_thickness_history, 'b-', linewidth=2.5)
ax1.set_xlabel('Time (days)', fontsize=12)
ax1.set_ylabel('Wax Thickness (mm)', fontsize=12)
ax1.set_title('Wax Deposition vs Time', fontsize=13, fontweight='bold')
ax1.grid(True, alpha=0.3)

# 有效内径变化
ax2 = axes[1]
ax2.plot(t_hours/24, D_eff_history, 'r-', linewidth=2.5)
ax2.axhline(y=D_i*1000, color='g', linestyle='--', linewidth=1.5, label='Original Diameter')
ax2.set_xlabel('Time (days)', fontsize=12)
ax2.set_ylabel('Effective Diameter (mm)', fontsize=12)
ax2.set_title('Effective Diameter Reduction', fontsize=13, fontweight='bold')
ax2.legend()
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('wax_deposition.png', dpi=150, bbox_inches='tight')
plt.show()

六、停输再启动安全分析

6.1 停输温降分析

管道停输后，管内原油通过管壁向周围土壤散热，温度逐渐降低。停输温降决定了安全停输时间，是制定应急预案的重要依据。

6.2 停输温降数学模型

停输后管内原油温度变化可用集总参数法描述：

$$\rho c_{p}V\frac{dT}{dt}=-KL(T-T_g)$$

解得：

$$T(t)=T_g+(T_s-T_g)\exp \left(-\frac{t}{{\tau }_{th}}\right)$$

其中，热时间常数：

$${\tau }_{th}=\frac{\rho c_{p}A}{K}$$

6.3 再启动压力计算

停输后原油温度降低，粘度增加。再启动时需要克服的摩阻包括：

（1）牛顿流体

$$\Delta p=\frac{128\mu LQ}{\pi D^4}$$

（2）宾汉塑性流体

$$\Delta p=\frac{16{\tau }_{y}L}{3D}+\frac{128{\mu }_{p}LQ}{\pi D^4}$$

其中，${\tau }_y$为屈服应力，${\mu }_p$为塑性粘度。

6.4 Python仿真实现

# 停输分析参数
T_shutdown = 55  # 停输时油温，°C
t_shutdown_max = 72  # 最大分析时间，hours

# 热时间常数
tau_th = rho_oil * cp_oil * A / K

# 停输温降
t_stop = np.linspace(0, t_shutdown_max*3600, 500)  # seconds
T_stop = T_g + (T_shutdown - T_g) * np.exp(-t_stop / tau_th)

# 粘度-温度关系（Arrhenius型）
def viscosity(T):
    """计算给定温度下的粘度，mPa·s"""
    mu_ref = 50  # 参考粘度
    T_ref = 60   # 参考温度
    E_visc = 25000  # 粘流活化能，J/mol
    R_gas = 8.314
    return mu_ref * np.exp(E_visc/R_gas * (1/(T+273.15) - 1/(T_ref+273.15)))

# 屈服应力-温度关系
def yield_stress(T):
    """计算给定温度下的屈服应力，Pa"""
    if T > WAT:
        return 0
    else:
        return 50 * (WAT - T)  # 简化线性模型

# 计算再启动压力
def restart_pressure(T, Q_restart):
    """计算给定温度下的再启动压力，Pa"""
    mu = viscosity(T) / 1000  # 转换为Pa·s
    tau_y = yield_stress(T)
    
    # 宾汉塑性流体压降
    delta_p = 16 * tau_y * L / (3 * D_i) + 128 * mu * L * Q_restart / (np.pi * D_i**4)
    return delta_p

# 计算不同停输时间后的再启动压力
Q_restart = 0.3  # 再启动流量，m³/s
p_restart = [restart_pressure(T, Q_restart)/1e6 for T in T_stop]  # MPa

# 安全压力限值
p_max = 10  # 最大允许压力，MPa

print("\n" + "=" * 60)
print("停输再启动安全分析")
print("=" * 60)
print(f"停输时油温: {T_shutdown:.1f} °C")
print(f"热时间常数: {tau_th/3600:.2f} hours")
print(f"停输72小时后油温: {T_stop[-1]:.2f} °C")
print(f"停输72小时后粘度: {viscosity(T_stop[-1]):.2f} mPa·s")
print(f"停输72小时后启动压力: {p_restart[-1]:.2f} MPa")

# 可视化
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 停输温降
ax1 = axes[0]
ax1.plot(t_stop/3600, T_stop, 'b-', linewidth=2.5)
ax1.axhline(y=WAT, color='r', linestyle='--', linewidth=1.5, label='WAT')
ax1.axhline(y=T_pour, color='orange', linestyle='--', linewidth=1.5, label='Pour Point')
ax1.set_xlabel('Shutdown Time (hours)', fontsize=12)
ax1.set_ylabel('Oil Temperature (°C)', fontsize=12)
ax1.set_title('Temperature Drop During Shutdown', fontsize=13, fontweight='bold')
ax1.legend()
ax1.grid(True, alpha=0.3)

# 再启动压力
ax2 = axes[1]
ax2.plot(t_stop/3600, p_restart, 'r-', linewidth=2.5)
ax2.axhline(y=p_max, color='k', linestyle='--', linewidth=2, label=f'Max Allowable: {p_max} MPa')

# 找到安全停输时间
safe_idx = np.where(np.array(p_restart) < p_max)[0]
if len(safe_idx) > 0:
    t_safe = t_stop[safe_idx[-1]] / 3600
    ax2.axvline(x=t_safe, color='g', linestyle=':', linewidth=2, label=f'Safe Time: {t_safe:.1f}h')

ax2.set_xlabel('Shutdown Time (hours)', fontsize=12)
ax2.set_ylabel('Restart Pressure (MPa)', fontsize=12)
ax2.set_title('Restart Pressure vs Shutdown Time', fontsize=13, fontweight='bold')
ax2.legend()
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('shutdown_restart_analysis.png', dpi=150, bbox_inches='tight')
plt.show()

七、综合仿真案例

7.1 案例描述

某热油管道长200km，管径DN500，设计输量2000m³/d。原油凝点32°C，蜡析出温度42°C。管道采用聚氨酯泡沫保温，厚度80mm。要求：

1. 计算沿程温降，确定加热站间距
2. 优化保温层厚度
3. 预测蜡沉积量，确定清管周期
4. 分析停输安全时间

7.2 Python综合仿真代码

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

print("=" * 70)
print("主题053：石油管道输送换热仿真 - 综合案例分析")
print("=" * 70)

# ==================== 基础参数设置 ====================
# 管道参数
D_i = 0.5          # 管道内径，m
D_o = 0.52         # 管道外径，m
delta_ins = 0.08   # 保温层厚度，m
D_ins = D_o + 2 * delta_ins  # 保温层外径，m
H = 1.5            # 埋深（管中心），m
L_total = 200000   # 管道总长度，m

# 物性参数
rho_oil = 860      # 原油密度，kg/m³
cp_oil = 2200      # 原油比热容，J/(kg·K)
k_steel = 45       # 钢管导热系数，W/(m·K)
k_ins = 0.025      # 保温层导热系数，W/(m·K)
k_soil = 1.2       # 土壤导热系数，W/(m·K)

# 运行参数
Q_design = 2000 / 86400  # 设计输量，m³/s
T_inlet = 65       # 首站出站温度，°C
T_g = 8            # 地温，°C

# 原油特性
T_pour = 32        # 凝点，°C
WAT = 42           # 蜡析出温度，°C

print("\n【基础参数】")
print(f"管道长度: {L_total/1000:.0f} km")
print(f"管道内径: {D_i*1000:.0f} mm")
print(f"保温厚度: {delta_ins*1000:.0f} mm")
print(f"设计输量: {Q_design*86400:.0f} m³/d")
print(f"出站温度: {T_inlet:.0f} °C")
print(f"地温: {T_g:.0f} °C")
print(f"原油凝点: {T_pour:.0f} °C")
print(f"蜡析出温度: {WAT:.0f} °C")

# ==================== 计算总传热系数 ====================
A = np.pi * D_i**2 / 4
u = Q_design / A
Re = rho_oil * u * D_i / (50e-3)  # 假设粘度50mPa·s
Pr = 500

# 对流换热系数
Nu = 0.023 * Re**0.8 * Pr**0.3
h_i = Nu * 0.15 / D_i  # 假设油导热系数0.15 W/(m·K)

# 各环节热阻
R_i = 1 / (np.pi * D_i * h_i)
R_steel = np.log(D_o/D_i) / (2 * np.pi * k_steel)
R_ins = np.log(D_ins/D_o) / (2 * np.pi * k_ins)
R_soil = np.arccosh(2*H/D_ins) / (2 * np.pi * k_soil)

R_total = R_i + R_steel + R_ins + R_soil
K = 1 / R_total

print("\n【总传热系数计算】")
print(f"管内对流热阻: {R_i:.4f} (m·K)/W")
print(f"管壁热阻: {R_steel:.6f} (m·K)/W")
print(f"保温层热阻: {R_ins:.4f} (m·K)/W")
print(f"土壤热阻: {R_soil:.4f} (m·K)/W")
print(f"总传热系数K: {K:.2f} W/(m·K)")

# ==================== 沿程温降分析 ====================
a = K / (rho_oil * cp_oil * Q_design)
z = np.linspace(0, L_total, 1000)
T = T_g + (T_inlet - T_g) * np.exp(-a * z)

# 计算加热站间距
T_min = T_pour + 3  # 最低允许温度
L_station = np.log((T_inlet - T_g)/(T_min - T_g)) / a
n_stations = int(np.ceil(L_total / L_station))

print("\n【沿程温降分析】")
print(f"温降系数: {a*1000:.6f} 1/m")
print(f"终点温度: {T[-1]:.1f} °C")
print(f"最大加热站间距: {L_station/1000:.1f} km")
print(f"所需加热站数量: {n_stations}")

# ==================== 保温层优化 ====================
c_fuel = 0.08      # 燃料价格，$/kWh
c_ins = 400        # 保温材料价格，$/m³
eta_heater = 0.82
crf = 0.08 * (1.08)**20 / ((1.08)**20 - 1)

def total_cost(delta):
    if delta < 0.01:
        delta = 0.01
    D_ins_opt = D_o + 2 * delta
    R_ins_opt = np.log(D_ins_opt/D_o) / (2 * np.pi * k_ins)
    R_total_opt = R_i + R_steel + R_ins_opt + R_soil
    K_opt = 1 / R_total_opt
    
    T_avg = T_g + (T_inlet - T_g) * (1 - np.exp(-K_opt*L_total/(rho_oil*cp_oil*Q_design))) / (K_opt*L_total/(rho_oil*cp_oil*Q_design))
    q_loss = K_opt * (T_avg - T_g)
    
    C_heat = q_loss * L_total * 8760 / 1000 * c_fuel / eta_heater
    C_ins = np.pi * (D_ins_opt**2 - D_o**2) / 4 * L_total * c_ins
    
    return C_ins * crf + C_heat, C_heat, C_ins, q_loss

delta_range = np.linspace(0.03, 0.12, 50)
costs = [total_cost(d)[0] for d in delta_range]
delta_opt = delta_range[np.argmin(costs)]

print("\n【保温层优化】")
print(f"当前保温厚度: {delta_ins*1000:.0f} mm")
print(f"最优保温厚度: {delta_opt*1000:.1f} mm")

# ==================== 可视化 ====================
fig = plt.figure(figsize=(16, 12))
gs = GridSpec(3, 2, figure=fig, hspace=0.3, wspace=0.25)

# 1. 沿程温度分布
ax1 = fig.add_subplot(gs[0, :])
ax1.plot(z/1000, T, 'b-', linewidth=2.5, label='Oil Temperature')
ax1.axhline(y=T_g, color='g', linestyle='--', linewidth=1.5, label='Ground Temperature')
ax1.axhline(y=T_pour, color='r', linestyle='--', linewidth=1.5, label='Pour Point')
ax1.axhline(y=T_min, color='orange', linestyle='--', linewidth=1.5, label='Min Allowable')

for i in range(n_stations):
    z_s = i * L_station
    if z_s <= L_total:
        ax1.axvline(x=z_s/1000, color='gray', linestyle=':', alpha=0.5)
        ax1.plot(z_s/1000, T_inlet, 'ro', markersize=6)

ax1.set_xlabel('Distance (km)', fontsize=11)
ax1.set_ylabel('Temperature (°C)', fontsize=11)
ax1.set_title('Temperature Distribution Along Pipeline', fontsize=12, fontweight='bold')
ax1.legend(loc='upper right')
ax1.grid(True, alpha=0.3)

# 2. 保温层优化
ax2 = fig.add_subplot(gs[1, 0])
C_heat_list = [total_cost(d)[1]/1000 for d in delta_range]
C_ins_list = [total_cost(d)[2] * crf / 1000 for d in delta_range]
ax2.plot(delta_range*1000, C_heat_list, 'r-', linewidth=2, label='Heating Cost')
ax2.plot(delta_range*1000, C_ins_list, 'g-', linewidth=2, label='Insulation Cost')
ax2.plot(delta_range*1000, np.array(costs)/1000, 'b-', linewidth=2.5, label='Total Cost')
ax2.axvline(x=delta_opt*1000, color='k', linestyle='--', alpha=0.7)
ax2.set_xlabel('Insulation Thickness (mm)', fontsize=11)
ax2.set_ylabel('Annual Cost ($k)', fontsize=11)
ax2.set_title('Insulation Optimization', fontsize=12, fontweight='bold')
ax2.legend()
ax2.grid(True, alpha=0.3)

# 3. 热损失分布
ax3 = fig.add_subplot(gs[1, 1])
q_list = [total_cost(d)[3] for d in delta_range]
ax3.plot(delta_range*1000, q_list, 'b-', linewidth=2.5)
ax3.axvline(x=delta_opt*1000, color='r', linestyle='--', linewidth=1.5, label=f'Optimal: {delta_opt*1000:.1f}mm')
ax3.set_xlabel('Insulation Thickness (mm)', fontsize=11)
ax3.set_ylabel('Heat Loss (W/m)', fontsize=11)
ax3.set_title('Heat Loss vs Insulation', fontsize=12, fontweight='bold')
ax3.legend()
ax3.grid(True, alpha=0.3)

# 4. 停输温降
ax4 = fig.add_subplot(gs[2, 0])
tau_th = rho_oil * cp_oil * A / K
t_stop = np.linspace(0, 72*3600, 200)
T_stop = T_g + (T_inlet - T_g) * np.exp(-t_stop / tau_th)
ax4.plot(t_stop/3600, T_stop, 'b-', linewidth=2.5)
ax4.axhline(y=WAT, color='r', linestyle='--', linewidth=1.5, label='WAT')
ax4.axhline(y=T_pour, color='orange', linestyle='--', linewidth=1.5, label='Pour Point')
ax4.set_xlabel('Shutdown Time (hours)', fontsize=11)
ax4.set_ylabel('Temperature (°C)', fontsize=11)
ax4.set_title('Temperature Drop During Shutdown', fontsize=12, fontweight='bold')
ax4.legend()
ax4.grid(True, alpha=0.3)

# 5. 参数敏感性
ax5 = fig.add_subplot(gs[2, 1])
params = ['Flow Rate', 'Inlet Temp', 'Ground Temp', 'Insulation']
base_values = [Q_design*86400, T_inlet, T_g, delta_ins*1000]
sens_range = np.linspace(0.8, 1.2, 20)

for i, (param, base) in enumerate(zip(params, base_values)):
    T_ends = []
    for factor in sens_range:
        if i == 0:  # 流量变化
            Q_sens = base * factor / 86400
            a_sens = K / (rho_oil * cp_oil * Q_sens)
        elif i == 1:  # 进站温度
            T_in_sens = base * factor
            a_sens = a
            T_end = T_g + (T_in_sens - T_g) * np.exp(-a_sens * L_total)
            T_ends.append(T_end)
            continue
        elif i == 2:  # 地温
            T_g_sens = base * factor
            T_end = T_g_sens + (T_inlet - T_g_sens) * np.exp(-a * L_total)
            T_ends.append(T_end)
            continue
        else:  # 保温厚度
            delta_sens = base * factor / 1000
            D_ins_sens = D_o + 2 * delta_sens
            R_ins_sens = np.log(D_ins_sens/D_o) / (2 * np.pi * k_ins)
            R_total_sens = R_i + R_steel + R_ins_sens + R_soil
            K_sens = 1 / R_total_sens
            a_sens = K_sens / (rho_oil * cp_oil * Q_design)
        
        T_end = T_g + (T_inlet - T_g) * np.exp(-a_sens * L_total)
        T_ends.append(T_end)
    
    ax5.plot(sens_range*100, T_ends, linewidth=2, label=param)

ax5.axvline(x=100, color='k', linestyle='--', alpha=0.5)
ax5.set_xlabel('Parameter Value (% of Base)', fontsize=11)
ax5.set_ylabel('Outlet Temperature (°C)', fontsize=11)
ax5.set_title('Sensitivity Analysis', fontsize=12, fontweight='bold')
ax5.legend()
ax5.grid(True, alpha=0.3)

plt.savefig('pipeline_comprehensive_analysis.png', dpi=150, bbox_inches='tight')
print("\n✓ 综合分析图已保存")
plt.show()

print("\n" + "=" * 70)
print("仿真完成！")
print("=" * 70)