去年我在做一个区域综合能源系统的年度仿真时,第一次把电转气(Power-to-Gas,P2G)模块完整地写进MATLAB程序里。当时领导给我的任务很直接:风电出力富余的时候,别让电白扔了,看看做成氢气或者合成天然气能回收多少能量,经济上划不划算。说实话,一开始我对电转气的认识还停留在“电解水制氢”几个字上,真正上手建模才发现里面全是细节——电解槽电压效率怎么折、甲烷化反应要不要考虑温度、整年8760小时的时序怎么跑,这些问题把我折腾得不轻。这篇文章就把我整个技术拆解、MATLAB建模思路、程序实现和踩坑记录整理出来,希望给做综合能源、可再生能源消纳、储能仿真的朋友一点参考。
1. 电转气技术路径与MATLAB分析的整体思路
1.1 电转气到底转的是什么,能解决什么问题
电转气,英文叫Power-to-Gas,通俗讲就是用多余的电能去电解水,把电转化成氢气的化学能;更进一步,氢气还能和二氧化碳在甲烷化反应器里合成天然气。为什么这事儿值得做?因为电力的存储一直是个大问题,尤其是风电和光伏的波动性很强,电网不能凭空增加一个几十兆瓦的大电池来应对间歇性。氢气或者合成天然气就不一样了,它可以大规模存储,能灌进天然气管道,也能在燃料电池或者燃气轮机里再发电,等于把电力的时间转移问题变成了燃气的储存和运输问题,这就把电网和气网打通了。
在综合能源系统研究里,电转气通常扮演两个角色:一个是削峰填补,在风电、光伏出力大但用电负荷低的时候,把富余电力转化成燃气,降低弃风弃光率;另一个是碳循环的纽带,把工业过程或者电厂捕集到的二氧化碳和氢气反应,生成甲烷,既消纳了可再生能源,又实现了碳的资源化利用。前者侧重氢能或者甲烷的系统调度,后者侧重碳捕集利用链路的耦合,但不管哪种路径,落到工程计算上都需要物质平衡、能量平衡和效率折算。
1.2 为什么用MATLAB而不是别的工具
我在这个项目里选MATLAB,核心原因有三点。第一,电转气涉及的方程不算复杂,但需要处理长时间序列和大量的敏感性分析,MATLAB的矩阵运算天然适合这种批量计算;第二,后续要跟电力系统调度、热力系统模型耦合,MATLAB和Simulink的生态在这个领域非常成熟,随便一个文献里的算法都能找到对应的工具箱或者开源代码做对照;第三,程序分析和结果可视化太方便了,plot、sensitivity分析、热力图都是现成的。
有人可能会说Python也能做,甚至开源框架更多。我不否认,Python在数据处理和机器学习上确实有优势,但电转气这块大量模型来自电力系统学科,很多老一辈研究者给的参考代码都是MATLAB写的,直接拿过来改参数、跑对比,效率反而更高。而且MATLAB的调试环境对工程计算确实友好,变量面板盯着看,哪里不对马上就能定位。所以这个项目最终定为:MATLAB脚本搭建电转气系统稳态模型,结合时序数据做整年运行仿真,再用优化函数做容量配置的初步分析。
1.3 程序分析的总体框架如何组织
我的程序框架没有一上来就堆复杂公式,而是分成了四层:数据层、单元模型层、系统耦合层和结果分析层。数据层负责读入风电出力、负荷、电价、气象等时间序列;单元模型层就是电解槽、甲烷化反应器、储氢罐、燃气轮机的单独数学模型;系统耦合层把这些单元连起来,按照能量流和物质流的方向逐小时递推;结果分析层负责算效率、经济指标,再输出绘图。这样做的好处是每一层都相对独立,改电解槽参数不会影响甲烷化模块,调电价数据也不用手忙脚乱地去翻主程序,后续扩展容量优化、多能互补调度都能直接往上加。
2. 电转气系统建模的数学基础与关键参数
2.1 电解槽模型:从电压到效率的折算
电解水制氢是电转气的第一道工序,也是最核心的耗能环节。工程上常见的电解槽有碱性电解槽、质子交换膜电解槽和固体氧化物电解槽,我项目里用的是质子交换膜电解槽的简化模型,因为它在负载波动时响应快,适合配合波动性风电一起运行。
电解槽的物理本质是电化学过程,一个关键参数是工作电压。理论分解电压是1.23伏,但实际运行时因为极化损失,工作电压通常在1.8到2.2伏之间。电压越高,意味着同样的电流下消耗的电能越多,但产氢速率不一定按比例提升,所以效率会下降。我建立模型时用的是工程经验公式,把工作电压表达成电流密度的函数,大致形式如下。
% 电解槽工作电压-电流密度关系 % j 是电流密度,单位 A/cm2 j = linspace(0, 2, 100); % 从空载到满载 E_rev = 1.23; % 可逆电压 r = 0.24; % 欧姆电阻系数 s = 0.09; % 极化系数 t = 0.02; % 极化指数 E_cell = E_rev + r .* j + s .* log(t .* j + 1); % 单电池电压这是一个非常典型的极化曲线模型,欧姆项是线性增加的,活化极化项是对数形式,整体趋势符合实测数据。有了单电池电压,我就能算电解槽的实际功率了,进而得到电解制氢的效率和产氢量。
这里必须提醒一个初学者容易犯的错误:效率的定义有两条路线,一条基于电压比,就是1.23除以实际工作电压,另一条基于能量,就是产氢的低热值能量除以实际输入电能。计算的时候一定要统一基准,我自己就曾经在一个报告里把两个效率混着用了,结果整年的能量平衡差了将近300兆瓦时,查了半天才发现是效率口径的问题。如果只是为了算系统能量转换,直接用产氢低热值能量除以耗电量最稳妥。
2.2 甲烷化反应器:物质平衡和能量损耗
电解产生的氢气有两个去向,一个是作为氢产品直接出售,另一个是送入甲烷化反应器,和二氧化碳反应合成甲烷。我会先算物质平衡。基础的甲烷化反应式是CO2加四摩尔氢气生成一摩尔甲烷和两摩尔水,这个反应是强放热过程,标准反应热大约是每摩尔甲烷释放165千焦。工程上比较关心的是氢碳比,通常控制在3到4之间,我模型里取3.5左右,留一点氢气余量,既保证反应速率,又不至于浪费太多氢气。
质量计算的事不难,但量纲很容易出错。摩尔的逻辑理清楚了就好办,一摩尔甲烷是16克,需要4摩尔氢气也就是8克,所以从氢气到甲烷的质量转化系数是2倍。把这个折算关系用程序实现就是几行的事。
% 甲烷化物质平衡 % m_H2_meth 进入甲烷化反应器的氢气质量,单位 kg/h % 根据 4H2 -> CH4 的化学计量,质量比 8:16 m_CH4 = m_H2_meth * 2; % 生成甲烷质量 m_CO2 = m_H2_meth * (44 / 8); % 需消耗二氧化碳质量 n_CH4_kmol = m_CH4 / 16; % 甲烷物质的量 kmol/h reaction_heat = n_CH4_kmol * 165; % 放热量,单位 MJ/h在能量计算上,甲烷化过程并不是百分百把氢气能量转进甲烷,能量转化效率通常在80%到85%左右,剩余部分以反应热排放出来。这部分热量如果直接扔掉,系统总效率就会很低;如果通过换热回收去预热或者供热,整体效率能往上提不少。我做仿真时会同时输出反应热流量,方便后续在热力子系统里对冲掉这部分余热。
2.3 系统能量流分析框架
把电解槽和甲烷化反应器串起来后,整个电转气的能量流就清晰了。输入是一批电力,中间经过电解制氢,产出氢气和氧气,氢气一部分直接存储销售,另一部分进入甲烷化反应器,再加上外部进来的二氧化碳,产出甲烷和水,过程中还有热量释放。我在MATLAB里用一个结构体数组来管理这些能量流和物质流,每个小时都记录电输入、氢气产量、甲烷产量、系统总效率和单位成本。
系统总效率的计算公式可以这样拆解:
P_elec = 1000; % 输入电功率 kW eta_elec = 0.70; % 电解效率 P_H2_lhv = P_elec * eta_elec; % 氢气低热值能量 kW % 假设全部进入甲烷化 eta_meth = 0.83; % 甲烷化能量效率 P_CH4_lhv = P_H2_lhv * eta_meth; % 甲烷低热值能量 kW eta_overall = P_CH4_lhv / P_elec; fprintf('系统电转甲烷总效率: %.2f%%\n', eta_overall * 100);这段代码跑出来的结果大约是58%。我拿这个数字去对比文献,基本在合理范围内:电转氢路线效率约60%到75%,电转甲烷路线效率约50%到65%,如果加上余热回收还能再高几个百分点。这也是为什么很多实际工程在推广电转气时主打冷热电联供,单纯看电转天然气的电效率并不惊艳,但把热用起来之后,全能源链的综合利用效率就非常有竞争力了。
3. MATLAB程序实现:从零搭建电转气仿真模型
3.1 整体程序架构与模块划分
我实际交付的程序包含一个主脚本和三个函数文件:电解槽模型函数、甲烷化模型函数、系统能量平衡计算函数。这种模块化组织方式对排查问题特别友好。主脚本负责加载数据、调用函数、存储结果、出图,每个函数只做单一职责,输入输出接口固定。举个例子,电解槽函数只接收电功率和环境温度,返回氢产量、耗电量、效率;甲烷化函数只接收氢气流量和二氧化碳供给,返回甲烷产量和反应热,其他一概不管。
这种思路对初学者来说可能有点过度设计,但等你要把模型扩大、更换参数、写论文做多次算例的时候就会知道,把所有逻辑糊在一个脚本里,后面改起来会让人崩溃。我在项目中后期加入了容量优化模块,就是因为当初分层清晰,不用重写核心功能,只加了一个目标函数和约束边界就行。
3.2 电解槽模块的关键实现
电解槽模块里最重要的两个输出是产氢量和运行效率。我的实现思路是给定输入电功率,先假设电解槽在允许的功率范围内,通过功率除以单电池电压得到电流,再用法拉第定律算产氢量。不过实际项目中,更常见的做法是根据额定功率和效率直接折算,因为极化曲线的细节参数往往拿不到,工程估算用效率折算已经足够。
实测下来,用额定效率折算的方式在整年仿真里误差可控,而且计算速度极快。下面是我实际使用的核心片段。
function [m_H2, P_input, eta] = electrolyzer(P_set, eta_rated, rho_H2) % 电解槽简化模型:按额定效率折算产氢量 % P_set 输入电功率 kW % eta_rated 额定电解效率 % rho_H2 氢气低热值能量 kWh/kg,约33.3 P_input = P_set; P_H2 = P_set * eta_rated; % 氢气的热值功率 m_H2 = P_H2 / rho_H2; % kg/h eta = eta_rated; end注意,这里有个细节:如果实际功率P_set低于电解槽最小运行功率,比如额定容量的10%以下,电解槽是不能稳定运行的。我就吃过这个亏,风电出力低的时候强行让电解槽运行,程序里算出了很低的产氢量,但实际设备早就停机了。后来我在调用处加了个约束,低于下限直接置零,才跟现场数据对得上。
3.3 甲烷化模块与系统耦合实现
甲烷化模块的逻辑相对简单,就是物质平衡加能量转化。我在代码里额外加入了二氧化碳供给量不足时氢气外送的处理逻辑,这样在整年仿真里不会因为某一小时CO2不够就报错中断。
系统耦合层的核心是逐小时循环,循环体里判断每个小时的电力分配:如果风电出力高于电网可接纳上限,富余电力优先给电解槽;如果电解槽功率超限,多余部分再考虑电锅炉或储能。这样能保证能量流完整闭环。完整循环的结构大概是这样。
for h = 1:8760 P_wind = wind_profile(h); P_grid_accept = grid_limit(h); P_surplus = P_wind - P_grid_accept; if P_surplus > 0 P_p2g = min(P_surplus, P_elec_max); % 制氢 [m_H2_total, ~, eta_elec] = electrolyzer(P_p2g, eta_rated, rho_H2); % 按比例分配氢气去向 m_H2_meth = alpha_meth * m_H2_total; m_H2_product = (1 - alpha_meth) * m_H2_total; % 甲烷化 [m_CH4, heat_meth] = methanation(m_H2_meth, CO2_supply(h)); % 记录结果 results(h, :) = [P_p2g, m_H2_total, m_CH4, heat_meth, eta_elec]; else results(h, :) = [0, 0, 0, 0, 0]; end end这段代码我在实际项目里跑整年数据,8760小时循环不到3秒就出结果,非常轻量。但如果你后续要跑容量优化,就不能用这么简单的循环嵌套,把目标函数改成向量化计算或者用时间序列矩阵直接算,速度还能再提升几十倍。
3.4 数据可视化与结果判读
程序分析如果没有图形输出,干看表格是看不出问题的。我先画了三个图:风电出力与电转气功率的时序对比图、电解槽和甲烷化的逐月产能量柱状图、全年系统效率的概率分布直方图。绘图代码不复杂,但有个小技巧我用了很久:先定义颜色和轴属性,再统一用tiledlayout做多子图,这样导出的图片格式统一,论文里直接能放。
figure tiledlayout(2, 1) nexttile plot(t_hours, wind_profile, 'LineWidth', 0.8) hold on plot(t_hours, P_p2g_series, 'LineWidth', 0.8) legend('风电出力', '电转气耗电功率') xlabel('时间 (h)') ylabel('功率 (kW)') nexttile bar(monthly_production); % 逐月产量 xlabel('月份') ylabel('产量 (kg)')从图形结果里最容易发现的问题有两个:一个是大量小时的电转气功率为零,说明风电富余时段和电解槽运行区间匹配得不好,要么电解槽容量偏大,要么运行下限太高;另一个是系统效率直方图出现双峰,说明部分时段电解槽在低效工况运行。看到这类现象,基本就知道下一步调试方向了。
4. 典型工况算例分析与参数敏感性评估
4.1 输入工况与算例场景设置
我用一套典型风电数据做了算例分析,场景是一台容量10兆瓦的电解槽配合一套2兆瓦的甲烷化装置。风电场的装机容量是50兆瓦,年利用小时数约2000小时,电网接纳上限设定为25兆瓦,这样大概有20%的发电量需要在本地消纳或转化。电价按峰谷平时段区分,氢气销售价格按35元每千克估算。整个仿真跑一年,时间分辨率取1小时。
这种场景在现实中很常见,尤其是北方风资源富集的地区,冬季供暖季风电大发,电网调峰压力大,电转气如果建设得当,白天制氢晚上发电,能很好地缓解高峰时段的气源紧张。但算例做出来之前,我一直对效益持保留态度,因为设备投资太高,单靠氢气销售收入回收周期太长,必须结合辅助服务、碳交易和副产品氧气销售才有希望。
4.2 全年产量与效率结果解读
算例跑完后的核心结果我列在下面,这些数字我都做了单位折算,方便大家对照自己的模型。
| 指标 | 数值 |
|---|---|
| 年富余电量 | 约32800 MWh |
| 电解槽实际消耗电量 | 约24600 MWh |
| 年产氢量 | 约590吨 |
| 其中进入甲烷化的氢量 | 约260吨 |
| 年产甲烷量 | 约520吨 |
| 年产甲烷热值当量 | 约7230 MWh |
| 电转甲烷全链效率 | 约29.4% |
电转甲烷全链效率比前面算的58%低了一大截,原因有两层。一是电解槽并不是全年都在满负荷运行,很多时段只有额定功率的30%到50%,实测电解效率从70%掉到55%左右;二是只有一部分富余电量被电解槽承接,其余富余电力依然被放弃了,这部分在总效率分摊里被摊薄。所以算完才知道,整年系统的实际性能不仅取决于单设备效率,更取决于设备容量配置和运行策略的匹配程度。
4.3 敏感性分析:哪个参数对结果影响最大
为了回答“项目先优化哪个方向”这个问题,我做了一组敏感性分析,方法很简单,把某个参数从基准值上下浮动20%,观察系统整体净收益的变化幅度。结果很有意思,影响最大的居然不是电解槽效率,而是电价峰谷差。电价峰谷差决定了电解槽在什么时段运行最划算,峰谷差越大,越值得在谷段多用电制氢、在峰段少用电。
第二重要的是氢气的销售价格,毕竟如果氢气直接对外卖,收益比走完甲烷化链条来得更快。第三才是电解槽效率。甲烷化能量效率对净收益的影响相对最小,因为它只影响一段中间过程,而且上下游的成本和收益占比更大。基于这个结论,我在后续优化里把电价曲线的获取放在最高优先级,电解效率改进反而放在后面。这也是为什么我一直强调,MATLAB仿真不只是算数字,关键是通过结构化分析找出系统的瓶颈。
5. 程序调试中的常见问题与排查技巧实录
5.1 量纲混乱导致能量平衡对不上
这个坑我几乎每个项目都要踩一次。电转气系统里涉及的参数太多,功率用千瓦、热量用兆焦、氢产量用千克、能量密度有低热值高热值之分,一旦没有统一单位,整年的能量平衡报表就会出现莫名其妙的差异。我后来给自己定了一条规矩:所有内部计算的功率统一用千瓦,能量统一用千瓦时,质量流量统一用千克每小时,转换系数全部写在程序开头,不散落到处。程序跑完后,我会单独算一层全年能量平衡,输入总电量,输出氢气能量、甲烷能量、反应热和损耗,几项相加必须在允许误差内。如果差得多,先查单位,再查公式。
有个更容易忽略的细节是低热值和高热值混用。氢气高热值约141.7兆焦每千克,低热值约120兆焦每千克,很多热力学教材用的是高热值,但工程效率计算基本都用低热值。我一开始混用过,结果电解效率算出来比理论极限还高,一看就不对。后来所有能量折算统一用低热值,才恢复正常。
5.2 迭代不收敛与初值选择问题
电转气系统和电网结合后,如果采用牛顿法或内点法做优化调度,经常会遇到不收敛。我的经验是,大部分不收敛问题不是算法本身不行,而是初值给得太离谱。比如电解槽功率初值给到容量上限之外,导致电压电流曲线出现负值,求解器直接就崩了。解决办法是先跑一版不考虑约束的线性调度,把结果作为非线性模型的初值,实测下来收敛率可以提高到九成以上。
还有一个技巧是把硬约束适度放软,比如把电解槽功率上下限改成惩罚项加入目标函数,这样迭代过程中不会因为越界直接被拉回,数值曲线更平滑,优化更容易收敛。在做年度仿真时,这个技巧尤其好用,因为8760小时的序列非常长,个别小时出现数值抖动就会传导到下游。
5.3 MATLAB版本和工具箱带来的坑
这个项目中期我换了电脑,从MATLAB R2021a升到了R2023b,结果原来跑得好好的脚本报了一堆警告,原因是某些数值函数被新版推荐替换,比如strsplit相关逻辑变了,还有struct数组的默认行为略有调整。虽然整体建模逻辑没受影响,但这种版本差异提醒我,交付程序的时候尽量用最基础的原生函数,不要依赖特别小众的新特性,否则合作方电脑上很可能跑不起来。编程时定期用warning把非致命提示先隐藏掉,方便聚焦真正的问题,但交付前一定要把警告重新打开看一遍,确认没有隐患。
另外,如果用到优化工具箱,注意检查求解器名称在不同版本里的兼容性。早期版本叫fmincon,这个就一直没变过,但某些混合整数求解器在不同版本里的参数设置方式有差异,写通用程序时要多做一次环境判断,或者干脆用自带的optimoptions统一设置,避免踩版本坑。
5.4 常见问题排查速查表
我整理了几个高频问题的对应排查方向,直接给到你们。
| 现象 | 可能原因 | 排查方法 |
|---|---|---|
| 年产氢量偏高 | 低热值高热值混用 | 检查能量折算系数是否统一 |
| 全链效率出现负值 | 反应热符号方向搞反 | 检查放热过程的正负号定义 |
| 电解槽全年运行时间过短 | 最小运行限制设置过高 | 降低最小功率约束至额定10%以下 |
| 优化迭代不收敛 | 初值越界 | 先用线性调度结果作初值 |
| 绘图坐标轴数据错乱 | 时间向量步长不均匀 | 检查采样间隔是否设置为小时 |
| 结果矩阵维度不匹配 | 8760与8761混淆 | 检查循环边界,是否多算或少算一小时 |
这个表是我每次做项目交接时都要附上的。因为程序的问题往往不是单一公式出错,而是几个小坑叠在一起,排查方向对了才能快速定位。
6. 扩展方向与个人实操体会
6.1 从稳态仿真走向混合储能联合优化
电转气模型搭完之后,最自然的扩展是和电池储能、蓄热装置做联合优化。我在最新一版程序里加入了锂电池储能模型,富余风电先看电池能不能存,存不下了再给电解槽,这样短时间尺度的功率波动由电池吸收,长时间尺度的能量转移交给氢气或甲烷,两者互补性非常强。MATLAB里实现这个联合优化,核心是在目标函数里加一个储能的充放电约束,还有SOC的状态转移约束,并不复杂,但运行时间会明显增加。如果你也做类似扩展,建议先把时间分辨率从小时改成典型日,先在日尺度上跑通程序,再放回整年序列。
6.2 我的几点实操体会
做了这么久电转气仿真,我越来越觉得,MATLAB在这里的核心价值不是算那些公式,而是帮你快速把杂乱无章的能源系统想法变成一个可量化的分析框架。很多参数从文献里抄来很简单,但放到特定风电曲线、特定电价体系下到底结果如何,只有建出程序跑一遍才知道。数学公式给人的是方向感,程序跑完给人的才是真实的工程认识。
一个小技巧是,每次跑完算例都把关键输出存成MAT数据文件,或者直接导出Excel,方便后面画图和写报告。我早期吃过亏,脚本跑完没保存中间结果,后面要重新分析时只能整个重跑一遍,浪费了大量时间。另外,加注释真的很重要,多花十分钟把每段代码的逻辑写清楚,三个月后再打开程序,你会感谢当时的自己。
最后我想说,电转气这个技术方向还在快速发展,电解槽的效率在逐年提升,甲烷化的催化剂也在改进,虽然当前测算出来的全链经济性还不能完全和传统储能竞争,但在高弃风率、高碳价的市场条件下,它已经开始显现独特价值。做技术分析的人,能通过MATLAB这块“计算试验田”把趋势量化出来,本身就是一件很有成就感的事。希望这篇文章里分享的思路和代码能帮你在自己的项目里少走几条弯路。