1.1 离心式压缩机特性曲线解析

离心式压缩机的核心性能指标通过压头-流量（$H-Q$）特性曲线来表征。在恒定转速下，随着气体流量的增加，理论压头通常呈下降趋势。该曲线反映了叶轮传递给气体的能量与气体处理量之间的非线性关系。在实际工程设计中，需重点关注喘振线（Surge Line）和阻塞线（Choke Line）。喘振线左侧为不稳定工况区，发生气流倒流并伴随剧烈振动；阻塞线右侧则因进气速度过高导致叶片攻角过大，效率急剧下降。

压头系数 $\psi$ 定义为叶轮出口相对速度与入口速度的比值相关项，其表达式可简化为：
$$H=\frac{u_2{V}_{u2}-u_1{V}_{u1}}{g}$$
其中 $u$ 为圆周速度，$V_u$ 为切向分量。为了量化转速变化对性能的影响，需应用相似定律（Affinity Laws）。当几何尺寸不变时，流量、压头和功率与转速的关系如下：
$$Q_2=Q_1\times \frac{n_2}{n_1}$$
$$H_2=H_1\times {\left(\frac{n_2}{n_1}\right)}^2$$
$$P_2=P_1\times {\left(\frac{n_2}{n_1}\right)}^3$$

以下示例代码演示了如何利用相似定律预测不同转速下的性能点，这对于裂解气压缩机的变频控制至关重要：

# -*- coding: utf-8 -*-
"""
离心压缩机相似定律计算模块
用于根据基准工况计算新转速下的流量、压头和功率。
"""
import math

def calculate_compressor_performance(Q1, H1, P1, n1, n2):
    """
    参数:
        Q1: 基准流量 (m^3/h)
        H1: 基准压头 (J/kg)
        P1: 基准轴功率 (kW)
        n1: 基准转速 (r/min)
        n2: 目标转速 (r/min)
    返回:
        Q2, H2, P2: 新工况下的参数
    """
    ratio = n2 / n1
    Q2 = Q1 * ratio
    H2 = H1 * (ratio ** 2)
    P2 = P1 * (ratio ** 3)
    return Q2, H2, P2

# 示例数据：某乙烯裂解气压缩机基准工况
base_Q, base_H, base_P, base_n = 8000.0, 45000.0, 1200.0, 6000.0
target_n = 7200.0

q_new, h_new, p_new = calculate_compressor_performance(
    base_Q, base_H, base_P, base_n, target_n
)

print(f"当前转速：{base_n} r/min -> 目标转速：{target_n} r/min")
print(f"流量变化：{base_Q:.1f} -> {q_new:.1f} m^3/h")
print(f"压头变化：{base_H:.1f} -> {h_new:.1f} J/kg")
print(f"功率变化：{base_P:.1f} -> {p_new:.1f} kW")

# 输出结果示例：流量增加20%，压头按平方律增长

1.2 系统管网阻力特性分析

裂解气压缩机的运行点由压缩机特性曲线与系统管网阻力曲线的交点决定。系统管路阻力主要来源于沿程摩擦阻力和局部阻力（阀门、弯头等）。对于不可压缩流体，单管段压力损失遵循达西-魏斯巴赫公式：
$$\Delta P=f\times \frac{L}{D}\times \frac{\rho V^2}{2}$$
由于流速 $V$ 与体积流量 $Q$ 成正比（$V=Q/A$），因此系统总阻力特性通常表现为二次曲线：
$$\Delta P_{sys}=K\cdot Q^2$$
其中 $K$ 为管路阻抗系数。

在乙烯装置设计中，需精确计算各段管线的 $K$ 值。局部阻力系数 $\zeta$ 与沿程摩擦系数 $\lambda$ 共同决定了总损失。若忽略气体可压缩性引起的密度变化，阻力特性近似为抛物线；但在高压比下，需考虑马赫数对密度的影响修正。

以下代码示例展示了如何根据管道几何参数和流体物性计算系统阻力曲线：

# -*- coding: utf-8 -*-
"""
管网阻力特性计算器
用于评估裂解气输送系统的压力降随流量变化的规律。
"""

def calculate_pipe_resistance(length, diameter_mm, roughness, flow_rate):
    """
    计算单根管线的压降 (Pa)。
    注意：此处简化了密度和摩擦系数，实际需查取穆迪图或Colebrook方程。
    """
    # 假设参数示例值，实际应动态获取
    density = 1.25  # kg/m^3, 乙烯常温近似
    
    # 流速计算 (Q in m^3/h to m/s)
    area_m2 = math.pi * (diameter_mm / 1000.0) ** 2 / 4.0
    velocity_ms = flow_rate / 3600.0 / area_m2
    
    # 摩擦系数经验估算 (Blasius公式适用层流过渡区)
    friction_factor = 0.079 * (density ** -0.25) 
    
    # 长度 L, 直径 D(m), 速度 V(m/s) -> 压降 Pa
    length_m = length
    diameter_m = diameter_mm / 1000.0
    
    # 简化计算逻辑，重点展示二次方关系
    resistance_coefficient = friction_factor * (length_m / diameter_m) * density / 2.0
    pressure_drop = resistance_coefficient * velocity_ms ** 2
    
    return pressure_drop

# 工况数据输入
pipe_L = 500.0      # 管道长度 m
pipe_D = 400.0      # 管道内径 mm
pipe_Q_list = [1000, 2000, 3000]  # 不同流量点

print("=== 系统管网阻力特性分析 ===")
for q in pipe_Q_list:
    delta_p = calculate_pipe_resistance(pipe_L, pipe_D, 0.045, q)
    print(f"流量 {q} m^3/h -> 压降 {delta_p:.2f} Pa (比例验证 K*Q^2)")

# 验证阻力与流量的平方关系
# 若 Q 翻倍，Delta P 应约为 4 倍

在运行优化阶段，需结合上述特性曲线进行匹配。例如，通过调节压缩机入口导叶（IGV）或转速改变压缩机工作线，使其始终位于高效区中心。同时，利用管网阻力分析结果，可识别瓶颈管线并实施扩径或加设旁路措施。对于乙烯装置这种大型连续流程，微小的参数波动都可能引起连锁反应，因此必须建立高精度的数学模型进行实时监测与预测控制。通过集成压缩机特性数据库与管网水力模型，可实现全厂能量系统的动态平衡优化，降低能耗并确保安全生产。此外，在装置大修期间，重新核算管路 $K$ 值对于验证设计余量至关重要。利用在线诊断系统采集的压力和流量数据，可以反推实际运行工况点偏离理论曲线的程度，从而及时调整操作参数。

【本章完】

2.1 复杂工况下的防喘振策略

在乙烯装置运行中，压缩机防喘振系统需应对原料气组分波动、下游用氢量变化等复杂工况。传统的固定旁路阀调节方式能耗高且响应滞后，现代控制策略采用前馈-反馈复合算法。

2.1.1 动态喘振线预测模型

由于气体组成变化会改变比热容和分子量，导致喘振边界移动。建立实时更新的喘振线方程：
$$P_{surge}=f(Q,T_{in},M_w,Z)$$
其中 $M_w$为平均分子量，Z为压缩因子。通过在线气相色谱数据修正模型参数可实现±2%的预测精度。

2.1.2 IGV主动调节算法

入口导叶角度 $\beta$与流量存在非线性关系：
$$\beta =k_{1}Q+k_2{Q}^{\log (Q)}$$
采用神经网络辨识系数 $k_1,k_2$，结合前馈补偿消除滞后。以下为防喘振控制器核心代码示例：

# -*- coding: utf-8 -*-
"""
复杂工况防喘振控制系统模块
实现动态喘振边界预测与IGV主动调节策略。
"""

import numpy as np
from scipy.optimize import minimize

class AntiSurgeController:
    def __init__(self, surge_map_file):
        """初始化控制器，加载喘振线特征点数据"""
        # 喘振点数据库：(流量m3/h, 压头J/kg)
        self.surge_points = np.array([
            [5000, 42000], [6000, 38000], 
            [7000, 35000], [8000, 31000],
            [9000, 27000]
        ])
        
        # IGV开度与流量关系系数（需现场辨识）
        self.k_IGV = np.array([0.015, -0.003, 0.8])
    
    def predict_surge_boundary(self, Q_current, T_inlet=293):
        """基于插值预测当前工况下的喘振压头"""
        # 双线性插值计算喘振边界
        Q_min, Q_max = self.surge_points[:, 0].min(), self.surge_points[:, 0].max()
        
        if Q_current < Q_min or Q_current > Q_max:
            return np.nan
        
        t1 = (Q_current - Q_min) / (Q_max - Q_min)
        H_surge = np.interp(Q_current, 
                          self.surge_points[:, 0],
                          self.surge_points[:, 1] + 
                          500 * (T_inlet - 293) / 293)
        
        return H_surge
    
    def calculate_IGV_position(self, Q_target, Q_current, margin=0.1):
        """计算目标IGV开度，保持安全裕量"""
        H_available = self.predict_surge_boundary(Q_target)
        if np.isnan(H_available):
            return 0.0
        
        # IGV开度与压头能力关系（简化模型）
        beta_opt = self.k_IGV[0] * Q_target + \
                   self.k_IGV[1] * (Q_target ** 2) + \
                   self.k_IGV[2]
        
        # 应用安全裕量（通常10-15%）
        beta_safe = max(0, min(90, beta_opt - margin))
        
        return beta_safe
    
    def emergency_bypass_control(self, current_H, measured_Q):
        """紧急旁路调节阀计算"""
        # 喘振裕量实时监测
        H_surge = self.predict_surge_boundary(measured_Q)
        if not np.isnan(H_surge):
            margin_ratio = (current_H - H_surge) / H_surge
            
            # 分级控制策略
            if margin_ratio < 0.05:  # 临界状态
                bypass_open = 100.0  # 全开旁路
                return bypass_open, "CRITICAL"
            
            elif margin_ratio < 0.15:  # 警告区
                bypass_open = min(100.0, 
                                max(0.0, 
                                    (margin_ratio - 0.05) / 0.2 * 60.0))
                return bypass_open, "WARNING"
            
            else:  # 正常区
                bypass_open = 0.0
                return bypass_open, "NORMAL"

# =============================================================================
# 防喘振系统联锁逻辑测试示例
# =============================================================================

def test_anti_surge_system():
    """模拟复杂工况下的防喘振响应"""
    
    controller = AntiSurgeController("surge_map.dat")
    
    # 工况1：正常操作点监测
    print("=== 工况1: 稳定运行测试 ===")
    Q_norm, H_operating = 7500.0, 32000.0
    margin = controller.predict_surge_boundary(Q_norm) - H_operating
    print(f"流量:{Q_norm} m³/h | 操作压头:{H_operating:.1f} J/kg")
    print(f"喘振边界预测:{margin:+.1f} J/kg ({margin/H_operating*100:.2f}%)")
    
    # 工况2：下游负荷骤降模拟
    print("\n=== 工况2: 负荷骤降响应 ===")
    Q_drop, H_drop = 4500.0, 39000.0
    margin = controller.predict_surge_boundary(Q_drop) - H_drop
    
    bypass_valve, status = controller.emergency_bypass_control(H_drop, Q_drop)
    print(f"新工况: Q={Q_drop}, H={H_drop}")
    print(f"喘振裕量:{margin:+.1f} J/kg")
    print(f"旁路阀开度:{bypass_valve:.1f}% | 状态:{status}")
    
    # 工况3：IGV调节响应测试
    print("\n=== 工况3: IGV主动调节 ===")
    for target_Q in [6000, 7500, 9000]:
        beta = controller.calculate_IGV_position(target_Q, 7500)
        H_boundary = controller.predict_surge_boundary(target_Q)
        print(f"目标流量:{target_Q:4d} m³/h | IGV开度:{beta:.1f}% "
              f"| 喘振边界:{H_boundary:.0f} J/kg")

# =============================================================================
# 运行结果示例（模拟输出）
# =============================================================================
"""
=== 工况1: 稳定运行测试 ===
流量:7500.0 m³/h | 操作压头:32000.0 J/kg
喘振边界预测:+850.0 J/kg (2.65%)

=== 工况2: 负荷骤降响应 ===
新工况: Q=4500.0, H=39000.0
喘振裕量:-1200.0 J/kg
旁路阀开度:100.0% | 状态:CRITICAL

=== 工况3: IGV主动调节 ===
目标流量:6000 m³/h | IGV开度:45.3% | 喘振边界:38250 J/kg
目标流量:7500 m³/h | IGV开度:62.1% | 喘振边界:32180 J/kg  
目标流量:9000 m³/h | IGV开度:78.5% | 喘振边界:24500 J/kg
"""

if __name__ == "__main__":
    test_anti_surge_system()

2.2 转子动力学振动监测方法

压缩机轴系振动是设备健康状态的关键指标。理解转子动力学基本原理对故障诊断至关重要：转子在临界转速附近运行时，弹性变形与惯性力耦合产生共振；不平衡量、不对中、齿轮损伤等缺陷都会激发特定频率的振动响应。

2.2.1 频谱分析与FFT处理

时域信号需转换为频域特征才能识别故障源。快速傅里叶变换将振动位移 $x(t)$分解为各谐波分量：
$$x(t)=a_0+\sum _{n=1}^N[a_n\cos (n\omega t)+b_n\sin (n\omega t)]$$

实际工程中通常使用汉宁窗（Hanning Window）减少频谱泄漏。以下为振动信号分析代码示例：

# -*- coding: utf-8 -*-
"""
压缩机转子动力学振动监测模块
实现频谱分析、故障特征提取及趋势预测功能。
"""

import numpy as np
from scipy.signal import welch, find_peaks, hann
import matplotlib.pyplot as plt

class RotorVibrationAnalyzer:
    def __init__(self, sample_rate=1024):
        """初始化振动分析仪
        
        参数说明：
            sample_rate: 采样频率 (Hz)，乙烯压缩机典型值800-4096 Hz
        """
        self.fs = sample_rate
        # 基准转速对应的转频（rpm → Hz）
        base_rpm = 5900.0
        self.rotor_freq = base_rpm / 60.0  # ≈ 98.33 Hz
        
    def fft_analysis(self, vibration_signal):
        """执行快速傅里叶变换并返回频谱特征"""
        
        N = len(vibration_signal)
        window = hann(N)
        signal_windowed = vibration_signal * window
        
        # FFT计算
        freqs, spectrum_power = welch(signal_windowed, 
                                     fs=self.fs, 
                                     nperseg=N//4,
                                     scaling='density')
        
        # 峰值检测 - 识别显著振动成分
        peaks, properties = find_peaks(spectrum_power, height=0.1*max(spectrum_power))
        
        # 提取转频及其倍频成分强度
        rotor_index = np.argmin(np.abs(freqs - self.rotor_freq))
        rotor_amplitude = spectrum_power[rotor_index]
        
        analysis_result = {
            'frequencies': freqs[:10],      # 前10阶频率点 (Hz)
            'amplitudes': spectrum_power[:10],
            'rotor_frequency_Hz': self.rotor_freq,
            'rotor_amplitude': rotor_amplitude,
            'peak_count': len(peaks),
            'spectral_energy': np.sum(spectrum_power)
        }
        
        return analysis_result
    
    def calculate_unbalance_severity(self, spectrum_data):
        """根据频谱数据评估不平衡严重程度"""
        
        # 转频处振幅与允许值比较（ISO 10816标准）
        rotor_freq = self.rotor_freq
        rotor_idx = np.argmin(np.abs(spectrum_data['frequencies'] - rotor_freq))
        
        if spectrum_data['amplitudes'][rotor_idx] < 2.5e-4:  # mm级别阈值
            severity = "NORMAL"
        elif spectrum_data['amplitudes'][rotor_idx] < 7.5e-4:
            severity = "MODERATE"  
        else:
            severity = "SEVERE - 需停机检查"
        
        return {
            'severity': severity,
            'recommendation': self._get_recommendations(severity)
        }
    
    def _get_recommendations(self, level):
        """生成维护建议"""
        recs = {
            "NORMAL": ["继续保持常规监测", "记录基线数据"],
            "MODERATE": ["增加监测频率至每日一次", 
                        "检查联轴器对中情况",
                        "验证润滑油品质"],
            "SEVERE - 需停机检查": [
                "立即安排计划停机",
                "进行激光对中测量",
                "必要时更换动平衡块"
            ]
        }
        return recs[level]

    def trend_prediction(self, historical_data, future_points=24):
        """基于历史数据的振动趋势外推"""
        
        if len(historical_data) < 10:
            return {"error": "数据点不足，无法进行趋势分析"}
        
        # 简单线性回归预测
        x = np.arange(len(historical_data))
        y = np.array(historical_data).flatten()
        
        slope, intercept = np.polyfit(x, y, deg=1)
        max_trend = slope * future_points
        
        if abs(max_trend) < 0.5:  # 阈值可调整
            prediction = "稳定，无显著恶化趋势"
        elif max_trend > 0 and max_trend < 2.0:
            prediction = f"缓慢上升，预计{future_points}小时后超限风险"
        else:
            prediction = "快速恶化，需立即干预"
        
        return {
            'slope': slope,
            'max_increase': max_trend,
            'prediction': prediction
        }

# =============================================================================
# 振动监测综合应用示例
# =============================================================================

def rotor_health_monitoring_demo():
    """演示完整的转子健康监测系统"""
    
    # 模拟采集的振动数据（mm，时间序列）
    np.random.seed(42)
    t = np.linspace(0, 1.0, 2048)      # 1秒采样
    base_vibration = 1e-5 * np.sin(2*np.pi*98*t) + \
                     0.3e-5 * np.random.randn(len(t))
    
    analyzer = RotorVibrationAnalyzer(sample_rate=1024)
    
    print("="*60)
    print("转子动力学振动监测综合报告")
    print("="*60)
    
    # 1. 频谱分析
    print("\n[1] FFT频谱分析结果:")
    results = analyzer.fft_analysis(base_vibration)
    print(f"  - 转频: {results['rotor_frequency_Hz']:.2f} Hz")
    print(f"  - 转频振幅: {results['rotor_amplitude']*1e6:.3f} μm")
    print(f"  - 频谱能量: {results['spectral_energy']/1e9:.4f} mV²/Hz")
    
    # 2. 不平衡评估
    print("\n[2] 不平衡严重程度评估:")
    severity_report = analyzer.calculate_unbalance_severity(results)
    print(f"  - 等级: {severity_report['severity']}")
    for rec in severity_report['recommendation']:
        print(f"    * {rec}")
    
    # 3. 趋势预测
    historical_vibrations = [1.2, 1.25, 1.28, 1.30, 1.32, 
                            1.34, 1.37, 1.40, 1.42, 1.45]
    print("\n[3] 振动趋势预测:")
    trend = analyzer.trend_prediction(historical_vibrations)
    print(f"  - 变化率: {trend['slope']:.6f} μm/点")
    print(f"  - 结论: {trend['prediction']}")
    
    # 4. 多阶频率识别
    print("\n[4] 主要振动频率成分:")
    for i, freq in enumerate(results['frequencies']):
        if abs(freq - analyzer.rotor_freq) < 2:  # 转频附近
            marker = "★"
        else:
            marker = ""
        print(f"  {marker}{freq:.1f} Hz (第{i+1}阶)")

# =============================================================================
# 运行输出示例（模拟）
# =============================================================================
"""
============================================================
转子动力学振动监测综合报告
============================================================

[1] FFT频谱分析结果:
  - 转频: 98.33 Hz
  - 转频振幅: 0.156 μm
  - 频谱能量: 0.4523 mV²/Hz

[2] 不平衡严重程度评估:
  - 等级: NORMAL
    * 继续保持常规监测
    * 记录基线数据

[3] 振动趋势预测:
  - 变化率: 0.008267 μm/点
  - 结论: 稳定，无显著恶化趋势

[4] 主要振动频率成分:
  ★98.3 Hz (第1阶)
  196.5 Hz (第2阶)
  294.8 Hz (第3阶)
"""

if __name__ == "__main__":
    rotor_health_monitoring_demo()

2.3 综合应用集成建议

将上述模块集成至DCS系统可实现：

1. 实时喘振裕量监控 - IGV控制器每50ms计算一次安全边界
2. 振动阈值报警分级 - 采用三级报警机制（黄色预警→橙色警告→红色停机）
3. 预测性维护窗口 - 基于振动趋势提前24-72小时发出检修建议

# DCS集成接口示例
class IntegrationLayer:
    def __init__(self, dcs_interface):
        self.dcs = dcs_interface
        
    def register_monitoring_tags(self):
        """注册需要监测的DCS标签"""
        tags = {
            'COMPRESSOR_FLOW': self.dcs.get_tag('FT-102'),
            'COMPRESSOR_DIFF': self.dcs.get_tag('DT-105'),
            'VIBRATION_X': self.dcs.get_tag('VTX-201'),
            'IGV_POSITION': self.dcs.get_tag('LV-302')
        }
        return tags
    
    def execute_control_loop(self):
        """执行控制循环"""
        # 获取实时数据
        flow = self.tags['COMPRESSOR_FLOW'].value
        diff_pressure = self.tags['COMPRESSOR_DIFF'].value
        
        # 计算喘振裕量并调节IGV
        margin = anti_surge_controller.predict_surge_boundary(flow)
        target_beta = controller.calculate_IGV_position(flow, flow * 1.05)
        
        # 采集振动数据进行分析
        vib_data = get_vibration_samples('VTX-201')
        spectrum = analyzer.fft_analysis(vib_data)
        
        return {
            'margin': margin,
            'IGV_target': target_beta,
            'vibration_level': spectrum['rotor_amplitude']
        }

if __name__ == "__main__":
    # 实际部署时连接真实DCS接口
    pass

通过上述技术体系的实施，乙烯装置可显著提升压缩机运行安全性与能效。建议每3个月更新一次喘振线数据库，每季度进行一次转子动力学校准验证。【本章完】