基于MATLAB的通风网络解算程序编制5.doc

本文发表在陕西煤炭2010年第29卷第6期，如需引用，请注明参考文献 姜诗明,裴绍宇,郄雷敏.基于MATLAB的矿井通风网络解算程序编制[J].陕西煤炭,2010,296. 基于MATLAB的矿井通风网络解算程序编制 姜诗明 （西安科技大学能源学院，陕西 西安 710054） E QQ 摘要介绍了矿井通风网络解算的数学模型，根据Cross迭代法基于MATLAB编制了矿井通风网路解算程序，并结合实例介绍了程序使用方法。表明用MATLAB编制通风网络解算程序具有编程简单、代码简洁、使用方便等优点，为通风系统分析及优化提供了工具。 关键词MATLAB；通风网络解算；程序； Program for Mine Ventilation Network Solution Based On MATLAB Jiang Shiming School of Energy Engineering,Xi’an university of science and technology,Xi’an 710054,China Abstract In this paper, a mathematics model of calculating ventilation network is introduced at first. According to Cross iteration programmed for mineventilation 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；Program 1引言 矿井通风网络解算是在已知矿井通风网络结构、分支风阻、风机特性的情况下，求解所有分支风量的过程。它作为通风安全管理定量分析工具，其作用贯穿矿井通风系统生命周期的始终，矿井的新建、改建、扩建都离不开它。通风网络解算对于模拟井下通风状况，预测网络工况在网络结构、风阻、风机的参数改变时的变化，一直起着很重要的作用。矿井通风网络解算对矿井通风系统进行理论分析或实验研究，辨识通风系统危险源，做出安全可靠性评价，制定安全技术措施，提高矿井通风系统安全可靠性有着重要的意义。 然而矿井通风网络解算程序编制却是一项困难的工作，目前通风网络解算软件一般利用VB 、C、C等高级程序语言编制，在整个软件生命周期中普遍存在着软件开发严重依赖于操作系统及编程语言，难以实现跨平台、跨语言共享代码；软件自身形成封闭系统，难以对现有通风软件进行二次开发，功能扩展性差；程序代码冗长难懂，开发及维护困难等问题。 MATLAB被称作第四代计算机语言，它具有语言简洁高效，简单易学，运算符、库函数及工具箱丰富，计算功能强大，绘图方便，扩展能力强大，可移植性好等特点。因此非常适合于编制矿井通风网络解算程序。 2矿井通风网络解算数学模型 矿井空气在通风网络中流动遵循节点风量平衡定律、回路风压平衡定律和阻力定律[1,2]。对于节点数，分支数的通风网络 1 2 3 其中 4 式（1）可改写为 5 将式（3）（4）（5）代入（2）得 6 式中，为分支风量列向量；为对风量列向量每个元素取绝对值所得的列向量； 为以为主对角元素的对角矩阵；为余树枝风量列向量；，为分支风压列向量；，为分支阻力列向量；，为风机风压列向量；，为分支位能差列向量；为分支风阻列向量；为以为主对角元素的对角矩阵；为回路风压代数和列向量。，为图G的关联矩阵；为图G的基本回路矩阵； 3通风网络解算算法及程序编制 式（6）为非线性方程组，通常采用迭代法求解其数值解。设第k次迭代后风量近似值为，将式（6）用泰勒级数展开并忽略二阶无穷小项得 7 8 式中，，为独立回路风量修正值列向量；为雅可比矩阵， 9 若雅可比矩阵为对角占优矩阵，即 10 则可略去雅可比矩阵非对角元素，式（8）简化为 11 可根据（6）式计算，，式（11）称为Cross迭代式，也称Scott-Hinsley迭代式。此算法计算步骤如下 ①输入通风网络结构及数据。 ②以风阻为权计算最小生成树及余树。 ③计算基本回路矩阵。 ④余树分支风量赋初始值，根据式（5）计算所有分支的初始风量。 ⑤计算回路风量修正值，修正回路风量，根据式（5）计算所有分支风量。 ⑥精度检验，即检查是否满足 12 为预设精度值。若满足上式则迭代终止，否则转⑤继续迭代计算。 ⑦输出计算结果。 根据以上计算步骤在MATLAB中编辑解算m文件solve.m如下[3,4] clear;clc;程序从这里开始 load in.txt通风网络数据文件（分支编号，始节点，末节点，风阻R，风机索引，自然风压（火风压）） load fan.txt风机数据文件（a0,a1,a2,a3,a4,a5） Vunionin,2,in,3;节点集合始节点和末节点的并集 VsizesizeV; mVsize1;节点数节点集合的行数 insizesizein; ninsize1;分支数in矩阵的行数 计算关联矩阵 Bzerosm,n;生成m行n列0矩阵 for i1n对所有分支循环 Bini,2,i1;in的第2列始节点 Bini,3,i-1;in的第3列末节点 end 计算最小生成树（prim算法） s[1];节点集合,加入第1个节点 tree[];树枝集合 fork1m-1循环m-1次，每次加入1树枝 [a,b]findBs,0;寻找s的所有关联分支 Rmin10000;最小风阻赋初值 bsizesizeb; fori1bsize1从1到关联分支数（b的行数） ifsizefindsinbi,21跳过此分支 end ifinbi,40始节点在节点集合中 s[s,inselect,3];加入始节点 elseifsizefindsinselect,30末节点在节点集合中 s[s,inselect,2];加入末节点 end end 计算余树 Ein,1;in的第1列，分支编号 E[tree],[];删除树枝所在行 cotreeE;余树 计算基本关联矩阵 temp1;参考节点 Btemp,[];删除第temp行 计算基本回路矩阵 B11B,cotree;余树所在列 B12B,tree;树所在列 B[B11,B12]; C11eyen-m1;n-m1行n-m1列单位矩阵 C12-B11*invB12; C[C11,C12];基本回路矩阵 将in矩阵按照余树在前，树枝在后的方法排序 in1in[cotree],;余树所在行 in2in[tree],;树所在行 in[in1;in2];余树在前，树在后 Hfzerosn,1;风机风压列向量（n行1列0矩阵） 风量赋初值 Qzerosn,1;风量列向量 Q1n1e-6; 开始迭代计算 k1; maxdq1e-5;最大回路风量修正值 maxf1e-5;最大回路不平衡风压 fork11e3控制最大迭代次数 maxdq0; maxf0; fori1n-m1对所有回路循环 f0.0;累加前赋初值 df0.0;累加前赋初值 forj1n对所有分支循环 fanindexinj,5;in矩阵的第5列为风机索引 iffanindex1 Hfjfanfanindex,1fanfanindex,2*Qjfanfanindex,3*Qj2fanfanindex,4*Qj3fanfanindex,5*Qj4fanfanindex,6*Qj5;计算风机风压 else Hfj0; end ffCi,j*inj,4*Qj*absQj-Hfj-inj,6;in矩阵第4列为风阻，第6列为自然风压火风压 dfdfCi,j*2*inj,4*absQj; end; dq-f/df;计算回路风量修正值 forj1n QjQjCi,j*dq;修正回路风量 end; f0.0;累加前赋初值 forj1n ffCi,j*inj,4*Qj*absQj-Hfj-inj,6;重新计算回路不平衡风压 end; if absdqmaxdq; maxdqabsdq;计算最大回路风量修正值 end; if absfmaxf; maxfabsf;计算最大回路不平衡风压 end; end; ifmaxdq1e-6 end end; k HRdiagin,4*diagabsQ*Q;阻力列向量 HHR-Hf-in,6;总风压列向量 out[in,Hf,HR,H,Q]输出结果 B*Q检查误差 C*H 输出网络解算数据 fprintf 分支编号 始节点 末节点 分支类型 风量 风阻 风压 阻力 风机风压 分支位能差\n\n; fori1n fprintf10d10d10d10d15.4f15.4f15.4f15.4f15.4f15.4f,ini,1,ini,2,ini,3,ini,5,Qi,ini,4,Hi,HRi,Hfi,ini,6; fprintf\n; end程序到这里结束 4实例分析 如下图所示通风网络，风阻列向量，安装于分支7中的风机风压特性曲线为 。 图1通风网络图 在solve.m文件同一路径下编辑通风网络数据文件in.txt，如下表前6列所示，一行中每个数据用空格或Tab键隔开，每一行用回车符隔开，共7行6列。编辑风机数据文件fan.txt为1046.35-0.85000，每个数据用Tab键隔开，共1行6列。运行solve.m输出解算结果如下表所示 表1通风网络解算数据表 分支编号 始节点 末节点 风阻 风机编号 自然风压 风机风压 阻力 风压 风量 1 1 2 0.375 0 0 0 163.5340 163.5340 20.8828 2 2 3 0.15 0 0 0 22.2382 22.2382 12.1760 3 2 4 2.0 0 0 0 151.6166 151.6166 8.7068 4 3 5 4.6875 0 0 0 376.2194 376.2194 8.9588 5 3 4 12.5 0 0 0 129.3784 129.3784 3.2172 6 4 5 1.7361 0 0 0 246.8410 246.8410 11.9240 7 5 1 0.5 1 0 780.0369 218.0453 -561.9916 20.8828 5结语 Cross迭代法是目前通风网路解算软件普遍采用的方法，运用Cross法基于MATLAB编制矿井通风网络解算程序具有编程简单、代码简洁易于理解、维护容易、易于扩展等特点。本文提供源程序在MATLAB R2008a中运行通过，具有通用性，可用于解算包含自然风压的通风网络，可供参考。 参考文献 [1]张国枢.通风安全学 修订版 [M].徐州中国矿业大学出版社,2007. [2]刘剑.等.流体网络理论[M].北京煤炭工业出版社,2002. [3]徐金明.MATLAB实用教程[M].北京清华大学出版社;北京交通大学出版社,2005. [4]邓薇. MATLAB函数速查手册[M]. 北京人民邮电出版社,2008. 源代码可在如下地址下载 http//