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

本文发表在陕西煤炭 2010 年第 29 卷第 6 期，如需引用，请注明参考文献 姜诗明,裴绍宇,郄雷敏.基于 MATLAB 的矿井通风网络解算程序编制[J].陕西煤炭,2010,296. 基于基于 MATLABMATLABMATLABMATLAB 的矿井通风网络解算程序编制的矿井通风网络解算程序编制 姜诗明 （西安科技大学能源学院，陕西 西安 710054） Emailbestjsm QQbestjsm 摘要摘要介绍了矿井通风网络解算的数学模型，根据 Cross 迭代法基于 MATLAB 编制了矿井 通风网路解算程序，并结合实例介绍了程序使用方法。表明用 MATLAB 编制通风网络解算 程序具有编程简单、代码简洁、使用方便等优点，为通风系统分析及优化提供了工具。 关键词关键词MATLAB；通风网络解算；程序； P P P Programrogramrogramrogram forforforfor MineMineMineMine VentilationVentilationVentilationVentilation NetworkNetworkNetworkNetwork SolutionSolutionSolutionSolution BasedBasedBasedBased O O O On n n n MATLABMATLABMATLABMATLAB Jiang Shiming School of Energy Engineering,Xi’an university of science and technology,Xi’an 710054,China AbstractAbstractAbstractAbstract In this paper, a mathematics model of calculating ventilation network is introduced at first. According to Cross iteration 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. KeywordsKeywordsKeywordsKeywords MATLAB；Ventilation network solution；Program 1 1 1 1 引言引言 矿井通风网络解算是在已知矿井通风网络结构、分支风阻、风机特性的情况下，求解所 有分支风量的过程。 它作为通风安全管理定量分析工具， 其作用贯穿矿井通风系统生命周期 的始终，矿井的新建、改建、扩建都离不开它。通风网络解算对于模拟井下通风状况，预测 网络工况在网络结构、风阻、风机的参数改变时的变化，一直起着很重要的作用。矿井通风 网络解算对矿井通风系统进行理论分析或实验研究， 辨识通风系统危险源， 做出安全可靠性 评价，制定安全技术措施，提高矿井通风系统安全可靠性有着重要的意义。 然而矿井通风网络解算程序编制却是一项困难的工作， 目前通风网络解算软件一般利用 VB 、C、C等高级程序语言编制，在整个软件生命周期中普遍存在着软件开发严重依赖 于操作系统及编程语言，难以实现跨平台、跨语言共享代码；软件自身形成封闭系统，难以 对现有通风软件进行二次开发， 功能扩展性差； 程序代码冗长难懂， 开发及维护困难等问题。 MATLAB 被称作第四代计算机语言，它具有语言简洁高效，简单易学，运算符、库函 数及工具箱丰富，计算功能强大，绘图方便，扩展能力强大，可移植性好等特点。因此非常 适合于编制矿井通风网络解算程序。 2 2 2 2 矿井通风网络解算数学模型矿井通风网络解算数学模型 本文发表在陕西煤炭 2010 年第 29 卷第 6 期，如需引用，请注明参考文献 姜诗明,裴绍宇,郄雷敏.基于 MATLAB 的矿井通风网络解算程序编制[J].陕西煤炭,2010,296. 矿井空气在通风网络中流动遵循节点风量平衡定律、回路风压平衡定律和阻力定律[1,2]。 对于节点数Vm，分支数En的通风网络,GV E 0BQ1111 0CH2222 diag diag HRRQQ3333 其中 HHRHfHn−−4444 式（1）可改写为 T C QC Q5555 将式（3） （4） （5）代入（2）得 0 TT CdiagCC diag F QC RC QC QHfHn−−6666 式中 12 ,,T n Qq qq⋯，为分支风量列向量；Q为对风量列向量每个元素取绝对值所得 的列向量； diag Q为以Q为主对角元素的对角矩阵； C Q为余树枝风量列向量； 12 ,,,T n Hh hh⋯，为分支风压列向量； 12 ,,,T n HRhR hRhR⋯，为分支阻力列向量； 12 ,,,T n Hfhf hfhf⋯，为风机风压列向量； 12 ,,,T n Hnhn hnhn⋯，为分支位能差列 向量； 12 , ,,T n Rr rr⋯为分支风阻列向量； diag R为以R为主对角元素的对角矩阵； 121 ,,, n m Ffff − ⋯为回路风压代数和列向量。 ijm n Bb ，为图 G 的关联矩阵； 1 ijn mn Cc − 为图 G 的基本回路矩阵； 1, 0, 1, ji ijji ji ev bev ev −⎧ ⎪ ⎨ ⎪ ⎩ 分支 空气流入节点 分支 与节点 不相关 分支 空气流出节点 1, 0, 1, j ijj j ei cei ei −⎧ ⎪ ⎨ ⎪ ⎩ 分支 在回路中且与回路反向 分支 不在回路中 分支 在回路中且与回路同向 3 3 3 3 通风网络解算算法及程序编制通风网络解算算法及程序编制 式（6）为非线性方程组，通常采用迭代法求解其数值解。设第 k 次迭代后风量近似值 为 k C Q，将式（6）用泰勒级数展开并忽略二阶无穷小项得 本文发表在陕西煤炭 2010 年第 29 卷第 6 期，如需引用，请注明参考文献 姜诗明,裴绍宇,郄雷敏.基于 MATLAB 的矿井通风网络解算程序编制[J].陕西煤炭,2010,296. 1 0 k CC kkk CCC C QQ F F QF QQ Q ∂ ∆ ∂ 7777 1 k CC kk CC C QQ F QF Q Q − ⎛⎞∂ ∆ −⎜ ⎟ ∂ ⎝⎠ 8888 式中， 1kkk CCC QQQ ∆−，为独立回路风量修正值列向量； C F Q ∂ ∂ 为雅可比矩阵， 2 TT diagC diag CC FdHf CRC QCC QdQ ∂ − ∂ 9999 若雅可比矩阵为对角占优矩阵，即 1 ,1,2,,1 n ii j CiCj j i ff inm qq ≠ ∂∂ − ∂∂ ∑ ≫⋯10101010 则可略去雅可比矩阵非对角元素，式（8）简化为 ,1,2,,1 / ci k k iC ici f Q qinm fq ∆ −− ∂∂ ⋯11111111 k iC f Q可根据（6）式计算， 1 2 n ii ijjj j cii fdHf c r q qdq ∂ − ∂ ∑ ，式（11）称为 Cross 迭代式，也 称 Scott-Hinsley 迭代式。此算法计算步骤如下 ①输入通风网络结构及数据。 ②以风阻为权计算最小生成树及余树。 ③计算基本回路矩阵。 ④余树分支风量赋初始值，根据式（5）计算所有分支的初始风量。 ⑤计算回路风量修正值，修正回路风量，根据式（5）计算所有分支风量。 ⑥精度检验，即检查是否满足 max,1,2,,1 Ci qinmε∆1跳过此分支 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 行 计算基本回路矩阵 本文发表在陕西煤炭 2010 年第 29 卷第 6 期，如需引用，请注明参考文献 姜诗明,裴绍宇,郄雷敏.基于 MATLAB 的矿井通风网络解算程序编制[J].陕西煤炭,2010,296. 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*Qj3 fanfanindex,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; 本文发表在陕西煤炭 2010 年第 29 卷第 6 期，如需引用，请注明参考文献 姜诗明,裴绍宇,郄雷敏.基于 MATLAB 的矿井通风网络解算程序编制[J].陕西煤炭,2010,296. 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,in i,5,Qi,ini,4,Hi,HRi,Hfi,ini,6; fprintf\n; end程序到这里结束 4 4 4 4 实例分析实例分析 如下图所示通风网络，风阻列向量 T R0.375 0.15 2.0 4.6875 12.5 1.7361 0.5 （，， ，，，， ）， 安装于分支 7 中的风机风压特性曲线为 2 1046.350.85hfqq−。 本文发表在陕西煤炭 2010 年第 29 卷第 6 期，如需引用，请注明参考文献 姜诗明,裴绍宇,郄雷敏.基于 MATLAB 的矿井通风网络解算程序编制[J].陕西煤炭,2010,296. 图 1通风网络图 在 solve.m 文件同一路径下编辑通风网络数据文件 in.txt， 如下表前 6 列所示， 一行中每 个数据用空格或 Tab 键隔开，每一行用回车符隔开，共 7 行 6 列。编辑风机数据文件 fan.txt 为 1046.35-0.85000，每个数据用 Tab 键隔开，共 1 行 6 列。运行 solve.m 输出解算结果如下表所示 表 1通风网络解算数据表 分支 编号 始节 点 末节 点 风阻 风机 编号 自然 风压 风机风 压 阻力风压风量 1120.375000163.5340163.534020.8828 2230.1500022.238222.238212.1760 3242.0000151.6166151.61668.7068 4354.6875000376.2194376.21948.9588 53412.5000129.3784129.37843.2172 6451.7361000246.8410246.841011.9240 7510.510780.0369218.0453-561.991620.8828 5 5 5 5 结语结语 Cross 迭代法是目前通风网路解算软件普遍采用的方法， 运用 Cross 法基于 MATLAB 编 制矿井通风网络解算程序具有编程简单、代码简洁易于理解、维护容易、易于扩展等特点。 本文提供源程序在 MATLAB R2008a 中运行通过， 具有通用性， 可用于解算包含自然风压的 通风网络，可供参考。 参考文献参考文献 [1]张国枢.通风安全学 修订版 [M].徐州中国矿业大学出版社,2007. [2]刘剑.等.流体网络理论[M].北京煤炭工业出版社,2002. [3]徐金明.MATLAB 实用教程[M].北京清华大学出版社;北京交通大学出版社,2005. [4]邓薇. MATLAB 函数速查手册[M]. 北京人民邮电出版社,2008. 本文发表在陕西煤炭 2010 年第 29 卷第 6 期，如需引用，请注明参考文献 姜诗明,裴绍宇,郄雷敏.基于 MATLAB 的矿井通风网络解算程序编制[J].陕西煤炭,2010,296. 源代码可在如下地址下载