基于matlab的通风网络解算程序编制

本文发表在陕西煤炭 2010 年第 29 卷第 6 期，如需引用，请注明参考文献姜诗明,裴绍宇,郄雷敏.基于 MATLAB 的矿井通风网络解算程序编制J. 陕西煤炭,2010,29(6).111 基于 MATLAB 的矿井通风网络解算程序编制姜诗明（西安科技大学能源学院，陕西 西安 710054）Email：QQ：摘要：介绍了矿井通风网络解算的数学模型，根据 Cross 迭代法基于 MATLAB 编制了矿井通风网路解算程序，并结合实例介绍了程序使用方法。表明用 MATLAB 编制通风网络解算程序具有编程简单、代码简洁、使用方便等优点，为通风系统分析及优化提供了工具。关键词：MATLAB；通风网络解算；程序；Program for Mine Ventilation Network Solution Based On MATLABJiang Shiming(School of Energy Engineering,Xian university of science and technology,Xian 710054,China)Abstract: In this paper, a mathematics model of calculating ventilation network is introduced at first. According to Cross iteration method programmed for mine ventilation network solution based on MATLAB, and taken examples to show how to use it. All of this proof that the program has many advantages, such as easy to use, code simple, convenient. It provides a tool to analysis and optimize the ventilation system.Keywords: MATLAB；Ventilation network solution；Program1 引言矿井通风网络解算是在已知矿井通风网络结构、分支风阻、风机特性的情况下，求解所有分支风量的过程。它作为通风安全管理定量分析工具，其作用贯穿矿井通风系统生命周期的始终，矿井的新建、改建、扩建都离不开它。通风网络解算对于模拟井下通风状况，预测网络工况在网络结构、风阻、风机的参数改变时的变化，一直起着很重要的作用。矿井通风网络解算对矿井通风系统进行理论分析或实验研究，辨识通风系统危险源，做出安全可靠性评价，制定安全技术措施，提高矿井通风系统安全可靠性有着重要的意义。然而矿井通风网络解算程序编制却是一项困难的工作，目前通风网络解算软件一般利用 VB 、C+ 、C# 等高级程序语言编制，在整个软件生命周期中普遍存在着软件开发严重依赖于操作系统及编程语言，难以实现跨平台、跨语言共享代码；软件自身形成封闭系统，难以对现有通风软件进行二次开发，功能扩展性差；程序代码冗长难懂，开发及维护困难等问题。MATLAB 被称作第四代计算机语言，它具有语言简洁高效，简单易学，运算符、库函数及工具箱丰富，计算功能强大，绘图方便，扩展能力强大，可移植性好等特点。因此非常适合于编制矿井通风网络解算程序。本文发表在陕西煤炭 2010 年第 29 卷第 6 期，如需引用，请注明参考文献姜诗明,裴绍宇,郄雷敏.基于 MATLAB 的矿井通风网络解算程序编制J. 陕西煤炭,2010,29(6).2 矿井通风网络解算数学模型矿井空气在通风网络中流动遵循节点风量平衡定律、回路风压平衡定律和阻力定律 1,2。对于节点数 ，分支数 的通风网络 ：VmEn,GVE(1)0BQ(2)CH(3)diagiR其中：(4)fn式（1）可改写为：(5)TCQ将式（3） （4） （5）代入（2）得：(6)()()0TTCdiagdiagFRHfn式中： ，为分支风量列向量； 为对风量列向量每个元素取绝对值所12(,TnQq得的列向量； 为以 为主对角元素的对角矩阵； 为余树枝风量列向量；diagQCQ，为分支风压列向量； ，为分支阻力列向量；12(,)TnHh 12(,)TnHRhR，为风机风压列向量； ，为分支位能差,ff nh列向量； 为分支风阻列向量； 为以 为主对角元素的对角矩阵；12()TnRr diag为回路风压代数和列向量。 ，为图 G 的关联矩阵；1,mFff ()ijmnBb为图 G 的基本回路矩阵；()ijnCc,01,j iij ij ievb分 支 空 气 流 入 节 点分 支 与 节 点 不 相 关分 支 空 气 流 出 节 点 1,0,jij jeici分 支 在 回 路 中 且 与 回 路 反 向分 支 不 在 回 路 中 分 支 在 回 路 中 且 与 回 路 同 向3 通风网络解算算法及程序编制式（6）为非线性方程组，通常采用迭代法求解其数值解。设第 k 次迭代后风量近似值为 ，将式（6）用泰勒级数展开并忽略二阶无穷小项得：kCQ本文发表在陕西煤炭 2010 年第 29 卷第 6 期，如需引用，请注明参考文献姜诗明,裴绍宇,郄雷敏.基于 MATLAB 的矿井通风网络解算程序编制J. 陕西煤炭,2010,29(6).(7)1()()0kCkkCCQFFQ(8)1()kCk式中， ，为独立回路风量修正值列向量； 为雅可比矩阵，1kkCCQCF(9)2TTdiagCdiagCFHfRQ若雅可比矩阵为对角占优矩阵，即：(10)1,2,1niijCCjiffnmq?则可略去雅可比矩阵非对角元素，式（8）简化为： (11)(),2,1)/ci kkiCcifQnq可根据（6）式计算， ，式（11）称为 Cross 迭代式，也()kiCfQ1ni iijjcfdHfr称 Scott-Hinsley 迭代式。此算法计算步骤如下：输入通风网络结构及数据。以风阻为权计算最小生成树及余树。计算基本回路矩阵。余树分支风量赋初始值，根据式（5）计算所有分支的初始风量。计算回路风量修正值，修正回路风量，根据式（5）计算所有分支风量。精度检验，即检查是否满足：(12)max,(12,1)Ciqnm为预设精度值。若满足上式则迭代终止，否则转继续迭代计算。输出计算结果。根据以上计算步骤在 MATLAB 中编辑解算 m 文件 solve.m 如下 3,4：clear;clc;%程序从这里开始load in.txt%通风网络数据文件（分支编号，始节点，末节点，风阻 R，风机索引，自然风压（火风压） ）load fan.txt%风机数据文件（ a0,a1,a2,a3,a4,a5）V=union(in(:,2),in(:,3);%节点集合：始节点和末节点的并集Vsize=size(V);m=Vsize(1);%节点数：节点集合的行数insize=size(in);本文发表在陕西煤炭 2010 年第 29 卷第 6 期，如需引用，请注明参考文献姜诗明,裴绍宇,郄雷敏.基于 MATLAB 的矿井通风网络解算程序编制J. 陕西煤炭,2010,29(6).n=insize(1);%分支数：in 矩阵的行数%计算关联矩阵B=zeros(m,n);%生成 m 行 n 列 0 矩阵for i=1:n%对所有分支循环B(in(i,2),i)=1;%in 的第 2 列始节点B(in(i,3),i)=-1;%in 的第 3 列末节点end%计算最小生成树（prim 算法）s=1;%节点集合,加入第 1 个节点tree=;%树枝集合for(k=1:(m-1)%循环 m-1 次，每次加入 1 树枝a,b=find(B(s,:)=0);%寻找 s 的所有关联分支Rmin=10000;%最小风阻赋初值bsize=size(b);for(i=1:bsize(1)%从 1 到关联分支数（b 的行数）if(size(find(s=in(b(i),2)=1)%跳过此分支endif(in(b(i),4)0)%始节点在节点集合中s=s,in(select,3);%加入始节点elseif(size(find(s=in(select,3)0)%末节点在节点集合中s=s,in(select,2);%加入末节点endend%计算余树E=in(:,1);%in 的第 1 列，分支编号E(tree,:)=;%删除树枝所在行cotree=E;%余树%计算基本关联矩阵temp=1;%参考节点B(temp,:)=;%删除第 temp 行%计算基本回路矩阵本文发表在陕西煤炭 2010 年第 29 卷第 6 期，如需引用，请注明参考文献姜诗明,裴绍宇,郄雷敏.基于 MATLAB 的矿井通风网络解算程序编制J. 陕西煤炭,2010,29(6).B11=B(:,cotree);%余树所在列B12=B(:,tree);%树所在列B=B11,B12;C11=eye(n-m+1);%n-m+1 行 n-m+1 列单位矩阵C12=-(B11)*(inv(B12);C=C11,C12;%基本回路矩阵%将 in 矩阵按照余树在前，树枝在后的方法排序in1=in(cotree,:);%余树所在行in2=in(tree,:);%树所在行in=in1;in2;%余树在前，树在后Hf=zeros(n,1);%风机风压列向量（n 行 1 列 0 矩阵）%风量赋初值Q=zeros(n,1);%风量列向量Q(1:n)=1e-6;%开始迭代计算k=1;maxdq=1e-5;%最大回路风量修正值maxf=1e-5;%最大回路不平衡风压for(k=1:1e+3)%控制最大迭代次数maxdq=0;maxf=0;for(i=1:n-m+1)%对所有回路循环f=0.0;%累加前赋初值df=0.0;%累加前赋初值for(j=1:n)%对所有分支循环fanindex=in(j,5);%in 矩阵的第 5 列为风机索引if(fanindex=1)Hf(j)=fan(fanindex,1)+fan(fanindex,2)*Q(j)+fan(fanindex,3)*(Q(j)2)+fan(fanindex,4)*(Q(j)3)+fan(fanindex,5)*(Q(j)4)+fan(fanindex,6)*(Q(j)5);%计算风机风压elseHf(j)=0;endf=f+C(i,j)*(in(j,4)*Q(j)*abs(Q(j)-Hf(j)-in(j,6);%in 矩阵第 4 列为风阻，第 6 列为自然风压(火风压)df=df+(C(i,j)*2*in(j,4)*abs(Q(j);end;dq=-f/df;%计算回路风量修正值for(j=1:n)Q(j)=Q(j)+C(i,j)*dq;%修正回路风量end;本文发表在陕西煤炭 2010 年第 29 卷第 6 期，如需引用，请注明参考文献姜诗明,裴绍宇,郄雷敏.基于 MATLAB 的矿井通风网络解算程序编制J. 陕西煤炭,2010,29(6).f=0.0;%累加前赋初值for(j=1:n)f=f+C(i,j)*(in(j,4)*Q(j)*abs(Q(j)-Hf(j)-in(j,6);%重新计算回路不平衡风压end;if abs(dq)maxdq;maxdq=abs(dq);%计算最大回路风量修正值end;if abs(f)maxf;maxf=abs(f);%计算最大回路不平衡风压end;end;if(maxdq1e-6)endend;%kHR=diag(in(:,4)*diag(abs(Q)*Q;%阻力列向量H=HR-Hf-in(:,6);%总风压列向量%out=in,Hf,HR,H,Q%输出结果%B*Q%检查误差%C*H%输出网络解算数据fprintf( 分支编号 始节点 末节点 分支类型 风量 风阻 风压 阻力 风机风压 分支位能差nn);for(i=1:n)fprintf(%10d%10d%10d%10d%15.4f%15.4f%15.4f%15.4f%15.4f%15.4f,in(i,1),in(i,2),in(i,3),in(i,5),Q(i),in(i,4),H(i),HR(i),Hf(i),in(i,6);fprintf(n);end%程序到这里结束4 实例分析如下图所示通风网络，风阻列向量 ，TR0.37512.0468751.360.5（ ， ， ， ， ， ， ）安装于分支 7 中的风机风压特性曲线为： 。2hfq本文发表在陕西煤炭 2010 年第 29 卷第 6 期，如需引用，请注明参考文献姜诗明,裴绍宇,郄雷敏.基于 MATLAB 的矿井通风网络解算程序编制J. 陕西煤炭,2010,29(6).图 1：通风网络图在 solve.m 文件同一路径下编辑通风网络数据文件 in.txt，如下表前 6 列所示，一行中每个数据用空格或 Tab 键隔开，每一行用回车符隔开，共 7 行 6 列。编辑风机数据文件fan.txt 为 1046.3 5 -0.85 0 0 0，每个数据用 Tab 键隔开，共 1 行 6 列。运行solve.m 输出解算结果如下表所示：表 1：通风网络解算数据表分支编号始节点末节点 风阻风机编号自然风压风机风压 阻力 风压 风量1 1 2 0.375 0 0 0 163.5340 163.5340 20.88282 2 3 0.15 0 0 0 22.2382 22.2382 12.17603 2 4 2.0 0 0 0 151.6166 151.616