矿井通风网络的数学建模

用图表示通风网络#

在矿井《矿井通风》课程中，我们已经学习到煤矿井下的巷道相互连接呈网络结构，而风流在巷道中流动，通过通风机提供动力，当自然分风无法满足风量分配要求时，进一步通过各种通风设施进行风量调节。一般将矿井通风系统抽象为单连通的有向图（Directed Graph），即通风网络图，如图1所示。

图1 矿井通风网络图示例

图（Graph）是图论（一个数学分支）的主要研究对象，它是由若干给定的节点（也称顶点）及连接两节点的分支（也称边）所构成的图形，这种图形通常用来描述某些事物之间的某种特定关系。通风网络图中不存在环（两端点相同的分支），但可能存在平行分支（有公共起点并有公共终点的两条分支）。图论已经系统地给出了图的数学表示和相关计算，因此可以据此构建通风网络图的数学模型。假设通风网络图中有 m 个节点，n 条分支，以列向量分别表示图的节点和分支所构成的集合：

其中 vi 表示编号为 i 的节点，ej 表示编号为 j 的分支。为了表述方便，也常用两个下标表示分支，如 eij 表示从节点 vi 到节点 vj 的分支。这里 V 和 E 以加粗、斜体的形式表示，表明他们不是单值标量，而是包含多个分量的向量或矩阵，下面的符号表示同样采用此法，将不再解释。

图中节点和分支的关系可用关联矩阵（Incidence matrix）表示：

对于有向图，上式中 bij 的取值为：

如图2中左侧为一个有向图，右侧为其关联矩阵。可以看出，该矩阵每列的代数和均为 0。

根据图论知识，连通图的关联矩阵的秩 rank(M) = m – 1 ，即关联矩阵 M 的各行是线性相关的。从关联矩阵 M 中去掉任意一行，得到各行线性无关的矩阵，称为基本关联矩阵，记作 B。

将各种通风参数看作是与通风网络图的分支 E 或节点 V 关联的数据，以列向量表示这些参数。定义在分支上的参数有风阻 R(e)、风量 Q(e)、通风阻力 hR(e)、风机风压 hf(e)、自然风压 HN(e) 等，定义在节点上的参数有节点绝对全压 Pt(v)、相对全压 ht(v)、重力位能 EP(v) 等。这些参数的带圆括号的上标 (e) 和 (v) 仅起指示作用，并无任何计算意义。其中上标 (e) 表示该参数是定义在分支上的向量，其元素个数为 n；上标 (v) 表示该参数是定义在节点上的向量，其元素个数为 m。这些向量的分量的下标标注方法与节点、分支的下标标注方法相同，如 Pt i 表示第 i 个节点的全压，Qi 表示第 i 个分支的风量，而 Qij 表示从节点 vi 到节点 vj 的分支的风量。分支中风量 Q(e)、风机风压 hf(e) 和自然风压 HN(e) 具有方向性，在建模时假定他们的方向都与分支方向相同，当他们实际的方向与分支方向相反时，这些参数将取负值。

构建通风系统方程组#

矿井通风系统中空气流动遵循风量守恒定律和能量守恒定律，通风网络建模的主要任务就是用方程以解析的形式描述这些定律。井巷内空气的流动是一个非常复杂的力学过程，为了便于建模和求解，必须对该过程做大幅度的简化，这些包括：将井巷内空气视作不可压缩的牛顿流体，将其流动视作一维定常流，空气的流态为紊流。

风量守恒定律#

风量守恒定律可分别针对节点或割集用不同的方法进行表述，这里针对节点进行表述：单位时间内流入与流出某节点的各分支的空气质量的代数和等于零。

由于已将空气看作是不可压缩的，这里直接用体积流量代替质量流量，得如下方程：

(1)

式中， qij 和 qjk分别为与节点 vj 关联的分支 eij 和 ejk 的风量，m3/s。

上式可表示为矩阵形式：

(2)

式中，0(v) 是元素个数等于图节点个数 m 的零向量。

能量守恒定律#

能量守恒定律可分别针对分支或回路用不同的方法进行表述，这里针对分支进行表述：通风网络中任一分支的通风阻力，等于该分支两端点间全压差、分支自然风压及分支中通风机风压的代数和。

该定律可表示为：

(3)

式中：

• hRij——分支 eij 的通风阻力，Pa；
• pi, pj——分别为节点 vi 和 vj 的全压，Pa；
• hf ij——分支 eij 的通风机风压，Pa；
• hNij——分支 eij 的自然风压，Pa。

根据通风阻力定律，当井巷内风流状态为紊流时，式(3)中的通风阻力项 hRij 为：

(4)

式中：

• qij——分支 eij 的风量，m3/s；
• rij——分支 eij 的风阻，N·s2/m8。

式(3)中的通风机风压项 hf ij 是由通风机的个体特性曲线决定的，流过通风机不同的风量，将会有不同的通风机风压 hf ij。为了便于计算，一般事先经过回归分析，将 hf ij 表示为风量 qij 的高次多项式：

(5)

式中：

• c0, c1, …, cn——多项式各项的系数；
• sign()——符号函数；当实际的风流方向与 eij 方向相同时，qij 为正，sign(qij) = 1；否则，qij 为负，sign(qij) = -1。

由于式(5)的长度较大，在后续的公式推导中，并不使用其展开式的形式，而仍然使用 hf ij。

本次设计中忽略式(3)中的自然风压项 hNij。

将式(4)代入式(3)，可得：

(6)

上式也可以表示为矩阵形式：

(7)

式中：

• diag(R(e)), diag(|Q(e)|) ——分别为由 R(e) 和 |Q(e)| 中各元素为对角元素，其余非对角元素值为 0 的对角矩阵；
• 0(e)——元素个数等于图分支个数 n 的零向量。

通风网络方程组#

联立风量守恒方程式(2)和能量守恒方程式(7)，即得到描述通风网络的方程组为：

(8)

方程组的定解条件#

实际的矿井通风网络解算多数是已知表征风网结构的关联矩阵 M、各个分支（巷道）的风阻 R(e)、自然风压 HN(e) 和通风机个体特性曲线 Hf(e)，求解各分支风量 Q(e) 和节点风压 Pt(v)（本设计仅对这种求解方式进行分析）。由于式(8)中 $\mathrm{diag}\left( \boldsymbol{R}^{{\left( \mathrm{e} \right)} \right) \mathrm{diag}\left( \left| \boldsymbol{Q}}{\left( \mathrm{e} \right)} \right| \right) \boldsymbol{Q}^{\left( \mathrm{e} \right)}$ 项和 Hf(e) 项都是分支风量 Q(e) 的非线性函数，因此整个方程组是一个非线性方程组，没法用线性代数中学到的线性方程组求解方法对此方程进行求解。

式(8)中的第 1、2 式分别包含 m、n 个方程。由于 rank(M) = m – 1，因此第 1 式中包含m - 1 个线性无关的方程。又由于 diag(M(e)) 这样的对角矩阵必然是满秩的，因此第 2 式包含 n 个线性无关的方程。这样线性无关的方程的总个数为 m + n – 1，而未知数包含 m 个节点全压和 n 个分支风量，共 m + n 个。未知数个数大于线性无关方程的个数，方程组是欠定的，将有无穷多个解。

为了使此方程组仅有有限个数的解，还需要再给定一个节点全压或分支风量值，通常是给定通风网络中大气节点的全压，即井口的大气压 P0。假设大气节点的编号为 r，则式(8)变为：

(9)

式(9)中的最后一式相对于为式(8)添加了边界条件。在式(9)中，线性无关方程组个数和未知数个数都是 m + n，因此该方程组是良态的，其具有有限数量的解。注意，式(9)并非只有一个解，这和抛物线方程常常有两个解的情形类似。

接下来对于通风网络的解算问题，已经变为求解式(9)表示的非线性方程组问题了。