多匝道隧道一维稳态通风 Matlab 数值计算

刘冰清

上海市政工程设计研究总院（集团）有限公司

0 引言

随着城市交通规划的需要，多匝道隧道在实际工程中出现得愈加频繁，如在建的苏州的金鸡湖隧道、上海的北横通道等。隧道通风计算目前以《公路隧道通风设计细则》（以下简称《细则》）[1]为标准，设计人员依据该规范可手工计算直通隧道的通风量。但是对于多匝道隧道，由于分合流节点的存在，使得隧道的风量计算变为网络化求解过程，需多次迭代计算，人工直接计算困难。

杨铠腾、赵红莉[2-3]利用Matlab建立了以质量流描述通风动力计算模型，并研究了隧道内的温度场与质量流之间的影响。但是计算模型中的三通局部阻力采用的是查表计算，软件的自动化程度不高。

本次研究采用Bassett理论公式[4]计算三通局部阻力，利用 Matlab编写通风网络计算软件 MTC，并对一个多匝道隧道火灾工况进行数值计算分析，检验 MTC的计算精度。

1 数值计算模型

1.1 模型简化

1）计算模型以一维稳态流体为对象。

2）隧道由单元、绝对方向、节点、回路、最简回路组5个要素组成。

单元：单元内的几何信息相同。

绝对方向：由用户定义单元内的绝对方向，当单元内的空气质量流量方向与绝对方向一致时为 ﹢，相反为-。

节点：任意节点处的空气质量流量守恒、能量守恒且静压唯一。

回路：从一个入口沿绝对方向经过一系列相连的单元、节点至另一个出口，所包含的单元、节点构成的集合。

最简回路组：由最少个数的回路构成的集合，该集合应包含隧道的所有单元、节点。

3）火灾时，烟气的发生量忽略不计。

4）火源点设置在单元绝对方向下游端，火源沿隧道纵向长度为0 m，空气经过火源点后，吸收全部火源的热释放量。

5）隧道中空气的温升、车辆的发热量、围护结构的传热量应满足能量守恒定律。

6）任意回路的全压应满足伯努利方程。

7）隧道内空气满足理想气体状态方程，隧道内空气的压力等于室外大气静压。

1.2 压差方程

隧道内单元的压差主要由8 部分组成，分别是通风阻力 ΔP r、交通通风力 ΔP t、射流风机升压力 ΔP j、轴流风机升压力ΔP fan、火源点节流阻力ΔP ft、动压差ΔP d、密度压差ΔP f、三通局部阻力ΔP si。

1）通风阻力 ΔP r、交通通风力 ΔP t、射流风机升压力 ΔP j三项压力计算按照《细则》[1]中 7.4.1-1～3、7.3.2、7.3.3、7.3.4-1～2及7.5.3-1式计算。

2）轴流风机升压力ΔP fan、火源点节流阻力 ΔP ft两项压力计算按照赵红莉[3]中 3-4和3-15式计算。

3）动压差ΔP d

4）密度压差ΔPf

式中：ρa为单元平均空气密度，kg/m3；ρt为室外的空气密度，k g/m3；L为单元长度，m ；g为重力加速度，9.8 m/s2。

5）三通局部阻力ΔP si

三通局部阻力与气流流向、断面积、角度、流量比例有关，软件中三通局部阻力计算采用 Bassett M D 等人的研究成果[4]计算。

6）回路净压差方程

回路内的所有单元各压差项之和应满足伯努利方程，回路的净压差应等于0，即满足式（3）。

式中：ΔP为回路的静压差，Pa。

1.3 节点流量方程

任意一个节点的净流量应为0，可按式（4）计算。

式中：Δm为空气净流量，kg/s为与该节点相连的i单元的空气质量流量，kg/s；ai为节点与相连i单元的方向，绝对方向流入节点为+，流出节点为-。

1.4 能量平衡方程

1）单元内的能量平衡方程

单元内的能量应满足能量守恒方程，即流入单元内的热量等于流出单元内的热量，即单元的净能量应为0，见图1所示。该单元净能量可按式（5～6）计算。

图1 单元内能量平衡图

式中：ΔE、ΔE F为净能量，kW；tu为单元绝对方向上游端的空气温度，℃ ；q为单元的总散热量，k W/m；tw为单元的围护结构平均温度，℃ ；C为单元的传热周长，m ；K为单元的传热系数，k W/K·m2；Qf为火源点热释放量，kW；dx为火源点沿隧道纵向长度，0 m。

单元内的空气密度可按（7～9）计算。

式中：Δρ为净密度，kg/m3；P为室外的大气压力，Pa。

2）节点的能量平衡方程

任意一个节点与其相连的单元能量应满足能量守恒方程，该节点的净能量应为0，可按式（10）计算。

式中：ti为与该节点相连的i单元端的空气温度，℃ 。

2 MTC计算流程

多匝道隧道风量的计算是一个关于风量的多元三次方程组，其中的风量、温度互相耦合，且方程组的系数随着风量的取值范围不同而变化。对一个由CN个单元，一组最简回路组（由 LN 个回路构成的集合）构成的多匝道隧道而言，其方程组表达式可按（11～12）描述。

针对方程组 F 和 G 的特点，M TC 软件以 Matlab为编写平台，求解流程见图2。

图2 求解流程图

3 MTC实例计算

3.1 隧道模型

隧道模型采用苏州元和塘西隧道南线作为 MTC的计算实例，隧道示意图见图3所示。图中矩形内数字表示单元编号，其中 1～20、101～102、201～202、301～302 为隧道单元编号，901～907 为排风竖井编号；箭头表示绝对方向；圆形 S 表示绝对方向为分流的节点；圆形 M 表示绝对方向为合流的节点；圆形 I 表示绝对方向气流入口的节点；圆形O 表示绝对方向气流出口的节点。仅连接两个单元的节点不做表示（如连接1、2单元的节点）。

图3 隧道示意图

主线行车方向为从单元 1 至 20，由 3 条匝道组成，其中单元101～102、201～202为出口匝道，301～302为进口匝道，并在主线出口处设置7处排风竖井（单元901～907），风井内设置轴流风机。

计算工况选择一个火灾排烟工况作为实例计算：

1）在隧道单元13绝对方向的下游末端发生火灾，火灾热释放率20 MW。

2）单元 7～13、301～302 中的车辆按照 165 pcu/km每条车道计算，其余隧道无车。

3）采用纵向分段排烟策略，新风就近从单元 301匝道入口补入，烟气从7个排风竖井排出。

4）开启主线隧道单元 13 中的 1 组 3 台射流风机，开启匝道单元 301中的2组 4台射流风机，开启 7个排风竖井的轴流风机。

1）车类型参数

表1 车类型表

2）几何信息参数

单元所属区段编号、长度、坡度、车道、面积、湿周、净高见表2所示。

表2 几何参数

3）气象参数

室外通风温度 3.7 ℃，大气压力 102410 Pa，室外空气密度1.29 kg/m3，隧道土壤温度 16.1 ℃，室外风速0 m/s。

4）风机参数

隧道内悬挂φ710型射流风机，出口风速v j=36.8 m/s，出口风量Qj=14.6 m3/ s，位置摩阻损失折减系数η=0.7。主线3台一组，匝道2台一组。

轴流风机的性能曲线4个状态点见表3。

表3 轴流风机状态点表

5）阻力系数

各单元阻力系数见表4所示。

表4 单元阻力系数

3.3 计算结果

各单元计算风量、温度项、密度项、压力项的 MTC计算结果见表5～6所示。

表5 MTC计算结果一

表6 MTC计算结果二

F、G方程一共迭代计算了21次，第21次计算的ΔF误差为 1.5282×10-9、ΔG误差为 2.5153×10-6，第 21次与第20次风量误差的绝对值为6.5335×10-4，总误差Δ为 6.5587×10-4，小于设定计算误差 0.001，符合计算精度要求。

4 结论

综合分合流局部阻力的研究编制了多匝道隧道一维稳态通风数值计算软件MTC，通过对多匝道隧道的计算分析，得到以下结论。

1）M TC最终的计算收敛误差为 6.5587×10-4，计算精度满足工程计算要求。

2）火灾工况时，M TC以质量流建立的计算模型可以更加全面的反应隧道内温度场对空气体积流量的影响，方便设计人员比较隧道的断面风速是否满足临街风速的要求。

3）分合流单元在火灾工况下的局部阻力占本单元全部阻力最高可达86%，表明局部阻力对隧道风量的影响较大。

4）研究成果表明，M TC计算软件能够有效的解决多匝道隧道的风量计算，且其中的局部阻力系数可按照理论公式在迭代中自行计算，不再需要依据经验或相关手册设定阻力系数。