多爆破工作面通风风量分配仿真：MATLAB实现多风机与风窗联合调节

干了这些年矿山通风相关的仿真项目，我越来越觉得多爆破工作面的风量分配是最磨人、也最有意思的一块。多个掌子面同时爆破，各条巷道的需风量不一样，再加上多台风机和多组风窗一起调，光靠经验和手算几乎不可能一次搞对，必须交给仿真去算。我这边用MATLAB整理过一套通风风量分配仿真例程，专门处理“多爆破工作面+多风机/风窗调节”的组合场景，运行一次就能看到全网风量分配结果、风机工况点、风窗调节量。今天把这套例程的思路、建模方法、核心代码和调参踩坑记录都摊开聊一遍，希望能给做矿井通风、隧道通风建模的朋友省点时间。

1. 项目概述：多个爆破工作面同时作业，通风为什么难调

1.1 多工作面爆破通风的核心矛盾

爆破工作面通风和普通巷道通风最大的不同在于，爆破瞬间会产生大量炮烟和有害气体，需要在一定时间内把这些污染物稀释到安全浓度以下。单工作面还好办，风量不够就加大风机；一旦变成多工作面同时爆破，问题就来了：每条巷道都要风，但是整个通风网络的总风量、总风压是有限的，各个分支之间存在强烈的耦合关系。

牵一发而动全身，是我对多工作面通风最直观的感受。调大A工作面的风量，可能会让B工作面的风量降下来；改变一台风机的转速，整个网络的压力分布都要重算。如果这时候还有几道风窗需要调节，那纯粹靠手算，算一个简单的三四个分支网络还行，分支数一多，基本就失控了。仿真工具在这里的意义，是替我们把全网的风量分配问题拆解成可迭代、可收敛的数学计算，快速给出每个调节装置的合理动作量。

1.2 这套MATLAB例程能解决什么问题

这套例程定位很明确：面向多爆破工作面的通风网络，同时支持多台风机的工况调节和多道风窗的阻力调节，最终输出满足需风量要求的风量分配方案。具体来说，它可以做三件事。

第一，根据各爆破工作面的炸药消耗量、巷道断面、通风距离、稀释时间等参数，计算每个工作面需要的风量。第二，对整个通风网络做风量分配解算，得到每条巷道的实际风量和风速。第三，当某些分支的风量不满足需风要求时，自动给出风机转速调节量或风窗增阻量，让全网风量重新分配，直到满足所有用风地点的需风量。

我需要说明一下，这套例程不是某个商业软件那种“开箱即用”的成品，它更像是一套可复用的计算框架。你要做的，是根据自己矿山的实际网络结构修改巷道参数表，然后运行解算和调节模块。但好处也很明显：整个计算过程是你自己可控的，每一步逻辑都摆在明面上，后续要扩展传感器数据联动、变频器控制策略，也都方便。

2. 通风风量分配的底层逻辑与建模

2.1 需风量计算：每个爆破工作面到底要多少风

爆破工作面的需风量不是拍脑袋定的，工程上一般是按几种不同要求分别计算，然后取最大值。我例程里重点考虑两个因素：一个是最低排尘风速要求，一个是炮烟稀释要求。

按最低排尘风速计算，公式很简单：

Q_dust = v_min * S;

其中v_min是最低排尘风速，岩巷一般取0.15~0.25 m/s，煤巷和半煤岩巷要高一些；S是巷道断面积。这个尺寸往往不大，但它是底线，风量再紧张也不能低于这条线。

按炮烟稀释时间计算，公式会稍微复杂一点，常见形式是：

Q_smoke = A * b / (C_limit * t);

式中A为一次爆破的炸药消耗量，b为单位炸药产生的有害气体量（一般取40 L/kg左右，可查规范），C_limit为炮烟稀释到的允许浓度，t为要求稀释时间。这个计算值通常比排尘风速算出来的大得多，是多工作面通风的“大头”。我这边还做了一层保护逻辑：当网络解算后某条用风分支的风量低于需风量一定比例时，程序会把它标记为“需调节分支”，并在后续调节模块中优先处理。

注意：这里的公式用于演示参数关系，实际工程取值必须以现行设计规范和设计手册为准。例程的价值在于把这套计算逻辑自动化，而不是替代规范。

2.2 通风网络解算：风量在网络上怎么分配

通风网络本质上是一个有向图，巷道是分支，分岔和汇合点是节点。风量分配遵循两个基本定律：节点风量平衡定律和回路风压平衡定律。前者说流入某个节点的风量等于流出该节点的风量；后者说在一个闭合回路中，所有分支的风压（包括摩擦风阻、局部阻力、风机升压）之和等于零。

用数学语言描述，就是一组非线性方程组。分支风压h满足二次阻力定律：

h = R * Q^2;

其中Q是分支风量，R是分支风阻。风阻由巷道断面积、周长、长度、摩擦阻力系数决定，例程里把它作为输入参数直接读入：

R = alpha * L * U / S^3;

alpha是摩擦阻力系数，L是巷道长度，U是巷道周长，S是断面积。巷道越长、断面越小，风阻越大，这个分支就越“难走”。

非线性方程组没有解析解，只能靠迭代。我例程里用的是通风网络解算中非常经典的Hardy Cross法，也就是风量迭代法。思路是先给每个闭合回路一个假定的校正风量，计算回路风压不平衡量，然后用这个不平衡量反过来修正那条回路里所有分支的风量，反复迭代直到全网风压平衡。

2.3 多风机/风窗调节为什么非仿真不可

如果不牵涉风机调节和风窗调节，上面的网络解算其实已经够用了。但实际场景里，光算出一个“自然分风”结果远远不够，因为自然分风的结果几乎不可能正好满足每个工作面的需风量。

这时候就要主动干预：要么调风机，要么调风窗。风机的干预手段通常是变频调速，改变风机特性曲线，进而改变它所在分支甚至整个网络的能量供给；风窗的干预手段是增加局部阻力，增大某条分支的风阻，把风量“压”到其他分支去。

问题在于：多台风机和多道风窗的调节是互相影响的。调A工作面的风机，可能让B工作面附近的风窗两端压差变大，原来定的风窗开度可能就不合适了。这种多变量耦合调节，手工计算几乎没法收敛到合理结果。仿真的做法是把风机特性、风窗阻力都纳入网络方程组，通过迭代解算和调节策略，让机器替我们完成这件“全局寻优”的工作。

3. MATLAB例程实现与关键参数配置

3.1 例程整体框架与输入参数表

整套例程的脚本结构不复杂，主要分四个模块。数据输入模块负责读巷道参数、风机参数和需风量参数；网络解算模块负责Hardy Cross迭代；调节计算模块负责判断哪些分支风量不满足要求并给出风机/风窗的调整量；结果输出模块负责画图和导出数据。

先看输入参数。我习惯将所有固定参数放到一个结构体里，便于统一管理：

net.nBranch = 9; % 分支数 net.nNode = 6; % 节点数 net.R = [0.35, 0.28, 0.42, 0.18, 0.31, 0.25, 0.38, 0.22, 0.30]; % 分支风阻 net.Qreq = [0, 0, 4.5, 0, 3.8, 0, 0, 3.2, 0]; % 各分支需风量，0表示无要求 net.fanBranch = [1, 5]; % 风机所在分支号 net.windBranch = [7]; % 风窗所在分支号

这里我特意把分支编号和风阻值都列出来，方便对照。Qreq数组中，非零值表示该分支是需要保证风量的用风分支，比如分支3、5、8分别对应三个爆破工作面。风阻的单位是Ns²/m⁸，实际值取决于巷道尺寸，我这里用的都是简化示例值。

风机参数用一条二次特性曲线描述，工程里常见表达方式是：

fan.h0 = 1200; % 风压特性曲线常数 fan.rq = 35; % 风压损失系数

风机升压H = h0 - rq * Q^2，其中h0近似于风机的最大静压，rq越大表示风机在大风量时静压跌落越厉害。变频调节时，根据相似定律，按转速比缩放h0和rq，这个后面会在调节模块里专门讲。

3.2 核心解算代码：Hardy Cross风量迭代法

网络解算是整套例程的核心。我先把分支、回路的关系整理成两个矩阵：一个是回路关联矩阵loopMap，行代表回路，列代表分支；另一个是回路方向矩阵dirMap，取值为1或-1，表示该分支在回路中的方向。

核心迭代代码如下：

Q = flowInit(net); % 初值：按节点流量平衡给出各分支风量 for iter = 1:200 maxDev = 0; for c = 1:nLoop sumP = 0; sumGrad = 0; for k = 1:nBranchInLoop(c) b = loopMap(c, k); d = dirMap(c, k); R = net.R(b); h = R * Q(b) * abs(Q(b)); if ismember(b, net.fanBranch) h = h - fanH_derived(b, Q(b), fan); % 风机升压为负项 sumGrad = sumGrad + 2 * R * abs(Q(b)) - fanH_deriv(b, Q(b), fan); else sumGrad = sumGrad + 2 * R * abs(Q(b)); end sumP = sumP + d * h; end dQ = -sumP / sumGrad; for k = 1:nBranchInLoop(c) b = loopMap(c, k); Q(b) = Q(b) + d * dQ; end maxDev = max(maxDev, abs(sumP)); end if maxDev < 1e-4 break; end end

这段代码里有几个地方必须提一下。第一，h = R*Q*abs(Q)而不是R*Q^2，这样做是为了让风压的符号跟着风量方向走，避免负风量分支的风压符号搞错。第二，分母sumGrad是回路风压对各分支风量的偏导数之和，本质上是牛顿法里面的雅可比项，如果漏掉风机特性的导数，多风机网络的收敛速度会明显变慢，甚至发散。第三，迭代终止条件看的是回路风压不平衡量的最大值，达到1e-4量级就认为收敛，这个阈值我一般不调，默认够用。

3.3 运行结果与曲线解读

跑完解算模块，程序会输出一组结果：各分支的风量、风速、风压、风机工况点、风窗阻力等。我习惯先看两个东西：全网风量分配表和风机工况曲线。

全网的分配结果一般用一个表格输出：

从这张表可以直观看到，分支3和分支5的风量都有缺口，分支8满足要求。这就是典型的“自然分风无法满足需风要求”场景，需要启动调节模块。分支3和分支5相差不大，可能通过调整风机转速就能补齐；分支5同时包含风机，它的调节会同时影响分支3和分支8，所以不能单独看待某一个分支，必须用全局调节策略。

风机工况曲线我一般用MATLAB里的plot把风机特性曲线画出来，再把解算得到的工况点标注上去。判断风机是否处于高效区，就看工作点能不能落在特性曲线的中段。工作点偏左，说明风机实际风量偏小，可能有风窗阻力太大、风机选型偏大的问题；工作点偏右，说明接近风机能力上限，再想加风就得谨慎，强行超风量运行容易烧电机。

4. 多风机/风窗联合调节策略

4.1 风机变频与风窗增阻的调节原理

多风机/风窗调节的底层逻辑，本质上是改变局部边界条件，让风量在全网重新分配。风机变频的依据是相似定律：转速从n1变到n2，风量近似按一次方正比变化，风压按二次方正比变化，功率按三次方正比变化。反映在风机特性曲线上，就是整个曲线被“压缩”或“拉伸”。

我例程里的变频处理比较直接：

ratio = n_new / n_rated; fan.h0_new = fan.h0 * ratio^2; fan.rq_new = fan.rq;

注意这里h0按转速比的平方缩放，rq不变。严格来说，风机特性曲线的缩放并不完全是这个关系，但对于轴流式风机在中段工作区的工程近似已经足够了。

风窗调节的本质是在分支上增加一个局部阻力，等效于增加该分支的风阻。风窗增阻之后，这条分支的实测风量会下降，多出来的风量会被“挤”到并联的其他分支去。例程中把风窗阻力作为附加风阻叠加到所在分支上：

net.R(b) = net.R(b) + deltaR_window;

4.2 自动调节迭代流程

考虑到多风机和风窗之间的耦合效应，我的调节模块没有直接一步到位，而是采用分步逼近的思路。具体流程是：先解算一次自然分风，找出所有“不满足需风量”的分支；然后根据缺风情况，对风机所在的分支尝试调整转速，对风窗所在的分支尝试调整阻力；每次调整后重新解算网络，再次检查所有分支的需风量；反复循环，直到所有用风分支的风量都达标，或者达到最大迭代次数。

这个流程用文字描述很简单，但代码里有一个细节特别重要：风机的调节范围有限，不是无限加大转速就能解决问题。当某台风机已经调到上限（比如转速比达到1.2），仍然无法满足下游需风量时，程序会转为增加并联风窗的阻力，把其他无关分支的多余风量压过来。

核心伪代码如下：

for adjustIter = 1:20 [Q, net] = solveNetwork(net, fan); shortBranch = find(Q < net.Qreq * 0.99); if isempty(shortBranch) break; end for b = shortBranch if ismember(b, net.fanBranch) fan.ratio(b) = min(1.2, fan.ratio(b) + 0.02); else net.R(findWindNear(b)) = net.R(findWindNear(b)) + 0.05; end end end

我故意把步长取得稍微小一点，每次只调一点点。因为工程上判断风量是否满足，不需要绝对精确到小数点后四位，只要在允许偏差范围内就行。步长太大会让调节结果在目标值附近来回震荡；步长太小又会导致迭代次数过多。0.02的转速步长和0.05的风窗阻力步长，是我测试下来比较稳的组合。

4.3 调节效果验证与评判指标

调理完以后，我通常会再跑一遍网络解算，然后从三个维度评判调节效果。

第一，用风分支的风量合格率。也就是所有需要保证风量的分支里，达到需风量要求的比例。这个指标最直观，我的目标永远是100%。

第二，全网平均风速和最高风速。爆破工作面巷道最低风速必须达标，最高风速又不能超限——风速过高不但吹得人难受，还会扬起巷道积尘，增加粉尘爆炸风险。所以我在例程里加了一个checkVelocity函数，把超过上限的风速标成警示红色。

第三，风机工况点的合理性。风机转速调高以后，风量上升，但功率是三次方关系上升，电耗增长非常明显。所以我会在输出结果里附一份能耗估算。比如一台风机从额定转速80%调到100%，理论上功率系数从0.512变成1，接近翻倍。如果某个工作面只是偶尔爆破需要大风量，长期大风量运行在经济上不一定划算，这时候就要考虑是不是用局部通风机接力更合理。这些问题仿真解决不了决策，但至少能给我们提供决策依据。

5. 常见问题与排查技巧实录

5.1 迭代不收敛或结果震荡

我做这套例程时踩过最大的坑，是风量迭代在几条并联分支之间来回“拉锯”，怎么都收敛不到设定精度。排查下来，大概率是初值给得太随意。Hardy Cross法的初值必须满足节点风量平衡，如果初值流一直在某个回路里明显失衡，迭代就容易绕圈。

处理办法有两种。一是用生成树法确定回路，然后给每棵树枝一个合理的初值风量，保证所有节点自动满足流量平衡。二是迭代后期放松收敛精度要求，从1e-4放宽到1e-3。工程上，风量偏差0.1 m³/s已经很小，过高的精度反而会让程序因为浮点误差一直“纠结”。

另外，风机特性的导数项漏写也会导致震荡。很多初版实现会把风机当成一个恒压源，也就是只减一个固定风压值，完全不考虑风压随风量变化的斜率。遇到这种写法，多风机网络特别容易发散，因为恒压源的数值雅可比是零，整个迭代修正量被严重放大。务必把-fanH_deriv这一项加进去。

5.2 负风量与风流反向问题

巷道网络解算过程中出现负风量，不一定代表程序错了，可能是该分支在给定条件下确实存在风流反向。爆破后，炮烟会沿回风系统排出，但如果某条并联巷道的通风动力不足，就可能出现工作面迎头风流停滞甚至反流，这是安全上绝对不允许的。

例程里我会把负风量直接视为风险信号。处理方式有两种：一是调整与该分支并联的调节设施，增加其他分支的阻力，把风量“逼”回来；二是增加该分支的风机动力。但这里有个经验：负风量分支如果其风阻极大，单纯靠增大并联风窗阻力来逼风，效果非常有限，这时候加装机或者改造风路是更合适的方案，仿真可以帮你判断这种改造带来的全网影响。

5.3 风窗调节越调越乱

还有一次我调试时遇到一个非常诡异的现象：本来是给分支7加风窗阻力，想把分支7的风量压一部分到分支3去，结果分支3的风量不升反降。查到最后发现，问题出在风窗所在的回路方向搞错了——分支7和分支3在同一个回路里的方向不一致，风窗增阻改变了回路风压，反而把风量从分支3“抽走”了。

这提醒我，在编写自动调节模块时，一定要先明确各分支之间的串并联关系，尤其是回路方向矩阵。建议每次调节前把网络拓扑和回路方向打印出来人工核对一遍，或者用颜色标注在图上检查。自动化程序再方便，第一步的拓扑正确性仍然要靠人工兜底。

6. 例程扩展与实用经验总结

6.1 从固定工况到动态调风

我现在的例程还支持一种扩展模式：把爆破时间序列加进去。因为多个工作面不一定同时爆破，错峰爆破完全可以减少同时需风量。例程里可以设置每个工作面的爆破时刻和炮烟稀释需求窗口：

schedule(1).blastTime = 10; % 第10分钟爆破 schedule(1).duration = 30; % 需在30分钟内完成稀释

这样程序可以根据时间窗口自动生成“需风量随时间变化”的动态目标曲线，再据此给出各个时段的风机调节建议。实际矿山上这样做省电效果特别明显，因为不需要让风机全天都维持最大爆破稀释风量，只需要在爆破后的那段时间内开足马力。

6.2 我的几点实操心得

这套例程我前后迭代了好几个版本，有些教训是教科书里不会写的。第一，参数文件的组织方式一定要规范，每个巷道名称、节点编号、风阻来源注释清楚，方便三个月后自己回来看还能看懂。用结构体嵌套比一堆零散的全局变量好维护得多。第二，每次修改网络拓扑或调节策略后，先跑一个已知结果的简单网络做回归验证，别直接上复杂网络。我曾经在一次重构后所有逻辑看起来都对，但结果莫名其妙，最后发现是回路方向矩阵某个元素从1写成了-1。第三，风速和风量的单位一定要统一。这个看似低级，但多人协作时真的很容易混，m³/s和m³/min差很多，一旦某处换算忘了，整个结果全废。

如果你正准备做类似的通风仿真，我建议先把单一工作面的简单网络跑通，确定Hardy Cross解算、风机特性、风窗阻力这三个核心模块都可靠，再逐步增加工作面数量和多调节设施的组合。把基础打牢，后面迁移到复杂网络只是增删参数的问题，不会牵动算法框架的改动。