CN104933230A - 考虑大气压影响的矿井采空区温度场仿真方法 - Google Patents

Abstract

本发明涉及一种考虑大气压影响的矿井采空区温度场仿真使用的进风巷道口处的全压值p_in，回风巷道口处的全压值p_out方法，将平面温度场模型方程联立并加入风压边界条件及热力学边界条件方程，在计算流体动力学软件FLUENT中进行数值求解，获得在一定地面大气压下采空区温度场仿真云图，得出大气压增大时，采空区内部气体的绝对压力升高并引起流速降低，进而因多场耦合导致热对流及传导效率低，温度分布的梯度更小的规律，与矿山实际状况吻合，误差小，精度高。

Description

考虑大气压影响的矿井采空区温度场仿真方法

技术领域

本发明涉及煤矿自燃火灾防治技术领域，具体涉及一种考虑大气压变化影响的矿井采空区温度场的仿真方法。

背景技术

煤自然发火是煤矿自燃火灾的主要原因，我国煤矿自燃发火形势十分严峻，全国约有56％的煤矿存在自燃发火问题。如何判断煤自然发火危险区域是煤炭安全生产的重要研究课题，煤自然发火预测预报作为其中的第一道防线，其重要性不言而喻。煤自然发火主要发生在采空区，因此采空区内煤自然发火是煤矿自燃火灾防治的重点，如何准确判定自燃危险区域是采空区自然发火预测预报的主要内容，而确定采空区自燃高温区域范围是煤矿自燃火灾防治措施制定的重要依据。

采空区自燃危险区域判定可分为直接测温法、气体分析法及计算机模拟法。直接测温法是在地面或井下向采空区可能发生自燃的地方打钻,在钻孔中埋设温度元件,根据测温结果来判定自燃危险区域，这种方法可靠性较好，但成本高、效率低；气体分析法具有工作量小、测定时间短、投资少等优点，但其探测误差较大；计算机模拟法是根据流体力学、热力学、传热学及多孔介质理论建立的采空区流-固耦合模型，结合相应的初始条件和边界条件，对采空区内流场、温度场及标志气体浓度场分布规律进行数值模拟，从而实现采空区自燃危险区域划分。计算机模拟法在理论上兼备了直接法及间接法的优点，受到学术界广泛重视。

目前，针对采空区温度场的数值模拟研究普遍只考虑了漏风量、环境温度、氧气浓度等因素，并未考虑采空区瓦斯抽放负压、大气压对计算结果的影响。而实际大气压的变化可达10kPa，瓦斯抽放负压变化也达几百帕，这使得数值模拟结果可靠性下降，对煤矿现场防灭火工作的指导力度减小

针对以上不足，需要提供一种考虑大气压变化影响的矿井采空区温度场的仿真方法。

发明内容

本发明目的在于根据实际气流流动及温度变化情况，进行采空区温度场的仿真，解决现有技术中的采空区温度场模型未考虑大气压的边界条件，使仿真结果不准确，从而影响采空区自燃危险区域判定的可靠性问题。

为了解决现有技术存在的问题，本发明采用的技术方案是提供了一种考虑大气压影响的矿井采空区温度场仿真方法，包括如下步骤：步骤①确定采空区边界，包括第一边界、第二边界和第三边界构成的矩形目标采空区，确定矩形目标采空区长度L值和宽度W值，确定坐标轴x、y，确定矩形目标采空区网格单元值，量取进风巷道口入风口宽度a值，测取气体的密度ρ_g值,确定采空区孔隙率ε值，确定采空区物体动力粘度μ值，确定采空区的渗透系数k值，确定气体扩散系数D_i,确定气体的比热C_pg值，确定煤的密度ρ_c和比热C_pc值，确定气体导热系数λ_g，确定煤的导热系数λ_c，确定采空区遗煤氧化反应产生的热量Q，确定采空区遗煤氧化速率r，确定工作面温度T₀值，确定工作面的热流密度q₀值，确定采空区边界温度T₁、T₂、T₃值，确定采空区边界的对流换热系数h₁、h₂、h₃值，测取x,y方向的气体流速u,v值，特点在于：步骤②建立以下方程获取考虑大气压影响的矿井采空区温度场仿真使用的进风巷道口处的全压值p_in，回风巷道口处的全压值p_out：

p_in＝s_n1·p_n1+s_m1·p_m1

p_out＝s_n2·p_n2+s_m2·p_m2

首先分时段实测多组进风巷道口处的绝对静压值p_n1，回风巷道口处的绝对静压值p_n2，进风巷道口处机械通风压力值p_m1，回风巷道口处机械通风压力值p_m2后，采用最小二乘法对上述两式中的s_n1、s_m1、s_n2、s_m2大气压影响因子进行参数辨识，得到大气压影响因子s_n1、s_m1、s_n2、s_m2值后，在测取当地大气压下同时，测取一组进风巷道口处的绝对静压值p_n1，回风巷道口处的绝对静压值p_n2，进风巷道口处机械通风压力值p_m1，回风巷道口处机械通风压力值p_m2后，代入上述两式中计算出考虑大气压影响的矿井采空区温度场仿真使用的进风巷道口处的全压值p_in，回风巷道口处的全压值p_out。

其中：步骤③

建立平面温度场模型方程：

连续性方程为：

$\frac{\∂u}{\∂x}+\frac{\∂v}{\∂y}=0$

式中，u,v分别为x,y方向的气体流速，m/s，

运动方程为：

$\frac{\∂\left({\mathrm{\ρ}}_{g}u\right)}{\∂t}+{\mathrm{\ρ}}_g(u\frac{\∂u}{\∂x}+v\frac{\∂u}{\∂y})=-\frac{\∂p}{\∂x}+\mathrm{\μ}(\frac{{\∂}^{2}u}{{\∂x}^2}+\frac{{\∂}^{2}u}{{\∂y}^2})-\frac{\mathrm{\μ}}{\mathrm{k\ϵ}}u$

$\frac{\∂\left({\mathrm{\ρ}}_{g}v\right)}{\∂t}+{\mathrm{\ρ}}_g(u\frac{\∂v}{\∂x}+v\frac{\∂v}{\∂y})=-\frac{\∂p}{\∂y}+\mathrm{\μ}(\frac{{\∂}^{2}v}{{\∂x}^2}+\frac{{\∂}^{2}v}{{\∂y}^2})-\frac{\mathrm{\μ}}{\mathrm{k\ϵ}}v$

式中，ρ_g为气体的密度，kg/m³，ε为采空区孔隙率，无量纲，p为气体静压，kPa，t时间，s，μ为采空区气体动力粘度，N·s/m²，k为采空区的渗透系数，无量纲，

组分方程为：

${\mathrm{\ρ}}_{g}u\frac{{\∂Y}_i}{\∂x}+{\mathrm{\ρ}}_{g}v\frac{{\∂Y}_i}{\∂y}={\mathrm{\ρ}}_g{D}_i\frac{{\∂}^2{Y}_i}{{\∂x}^2}+{\mathrm{\ρ}}_g{D}_i\frac{{\∂}^2{Y}_i}{{\∂y}^2}+S_i$

式中，D_i是气体扩散系数，m²/s，Y_i为气体成分浓度，mol/m³，S_i为耗氧速率，mol/m³·s，以温度T为变量的能量为：

$({\mathrm{\ϵ\ρ}}_g{C}_{\mathrm{pg}}+(1-\mathrm{\ϵ}){\mathrm{\ρ}}_c{C}_{\mathrm{pc}})\frac{\∂T}{\∂t}+{\mathrm{\ρ}}_g{C}_{\mathrm{pg}}(u\frac{\∂T}{\∂x}+v\frac{\∂T}{\∂y})={\mathrm{\λ}}_{\mathrm{eff}}(\frac{{\∂}^{2}T}{{\∂x}^2}+\frac{{\∂}^{2}T}{{\∂y}^2})+\mathrm{rQ}$

式中，C_pg为气体的比热，J/(kg·℃)，ρ_c,C_pc分别为煤的密度和比热，kg/m³，J/(kg·℃)，λ_eff为有效导热系数，W/m·℃，Q为采空区遗煤氧化反应产生的热量，kJ/kg，r为采空区遗煤氧化速率，kg/s,λ_eff的计算公式为:

λ_eff＝ελ_g+(1-ε)λ_c

式中，λ_g为气体导热系数，W/m·℃，λ_c为煤的导热系数，W/m·℃，

设定初始条件：

$\left\{\begin{array}{c}u=v=0,p=p_0,T=T_0\\ 0\≤x\≤L,0\≤y\≤W,t=0\end{array}\right.$

式中，p₀为采空区初始气压，kPa，T₀为初始温度，℃，u,v分别为坐标系中x,y方向的气体流速，m/s，L为采空区长度，m，W为采空区宽度，m，t为时间，s，

平面温度场模型风压边界条件方程：

工作面：p|_{x＝0,0≤y≤a}＝p_in，p|_{x＝0,W-a≤y≤W}＝p_out

p|_{x＝0,a≤y≤W-a}＝p_in-(p_in-p_out)·(y-a)/(W-2a)

第一边界： $\frac{\∂p}{\∂y}{|}_{y=\mathrm{0,0}\≤x\≤L}=0$

第二边界： $\frac{\∂p}{\∂x}{|}_{x=L,0\≤y\≤W}=0$

第三边界： $\frac{\∂p}{\∂y}{|}_{y=W,0\≤x\≤L}=0$

式中，p为气体静压，kPa，p_in为进风巷道口处的全压值，kPa，p_out为回风巷道口处的全压值，kPa，L为采空区长度，m，W为采空区宽度，m，a为进风巷道口入风口宽度，m，

平面温度场模型热力学边界条件方程：

工作面：T|_{x＝0,0≤y≤W}＝T₀， $-\frac{\∂T}{\∂x}{|}_{x=\mathrm{0,0}\≤y\≤W}=q_0$

第一边界：T|_{y＝0,0≤x≤L}＝T₁(x)， $-\frac{\∂T}{\∂y}{|}_{y=\mathrm{0,0}\≤x\≤L}=h_1(T_1-T_0)$

第二边界：T|_{x＝L,0≤y≤W}＝T₂(y)， $-\frac{\∂T}{\∂x}{|}_{x=L,0\≤y\≤W}=h_2(T_2-T_0)$

第三边界：T|_{y＝W,0≤x≤L}＝T₃(x)， $-\frac{\∂T}{\∂y}{|}_{y=W,0\≤x\≤L}=h_3(T_3-T_0)$

式中，T₀为工作面温度，℃，q₀为工作面的热流密度，W/m²，T₁、T₂、T₃分别为采空区第一边界1、第二边界2、第三边界3的边界温度，℃，h₁、h₂、h₃分别为采空区第一边界1、第二边界2、第三边界3的边界的对流换热系数，W/m²·℃，L为采空区长度，m，W为采空区宽度，m，

将以上平面温度场模型方程联立并加入风压边界条件及热力学边界条件方程，在计算流体动力学软件FLUENT中进行数值求解，获得在一定地面大气压下采空区温度场仿真云图。

其中：采空区边界温度T₁、T₂、T₃，通过U型设置于边界上的光纤温度传感器测量获得，每隔100m采集一个温度值。

本发明的有益效果在于：本发明考虑大气压影响的矿井采空区温度场仿真使用的进风巷道口处的全压值p_in，回风巷道口处的全压值p_out方法，将平面温度场模型方程联立并加入风压边界条件及热力学边界条件方程，在计算流体动力学软件FLUENT中进行数值求解，获得在一定地面大气压下采空区温度场仿真云图，得出大气压增大时，采空区内部气体的绝对压力升高并引起流速降低，进而因多场耦合导致热对流及传导效率低，温度分布的梯度更小的规律，与矿山实际状况吻合，误差小，精度高。

附图说明

下面结合附图和具体实施方式对本发明做进一步说明。

图1是本发明待仿真的采空区平面示意图；

图2是本发明目标采空区模型的网格划分示意图；

图3是本发明第一实施例地面大气压为99.225kPa时采空区温度场仿真云图；

图4是本发明第一实施例地面大气压为98.426kPa时采空区温度场仿真云图。

图中：1.第一边界，2.第二边界，3.第三边界，4.采空区，5.工作面，6.进风巷道口，7.回风巷道口。

具体实施方式

实施例一

参见图1、图2、图3，一种考虑大气压影响的矿井采空区温度场仿真方法，包括如下步骤：步骤①确定采空区边界，包括第一边界1、第二边界2和第三边界3构成的矩形目标采空区，确定矩形目标采空区长度L＝500m值和宽度W＝300m值，确定坐标轴x(m)、y(m)，通过点、线、面的顺序形成整个计算区域，建立采空区的计算域尺寸为500m*300m，并确定各离散单元的界面长度，即网格单元为0.5m，离散采用结构化网格，通过Gambit中的map网格划分方法生成了如图2所示的网格，网格划分成功后再由Gambit输出mesh文件，以便下一步用于Fluent软件求解器的识别和计算，量取进风巷道口6入风口宽度a＝5m，测取气体的密度ρ_g＝1.205kg/m³,确定采空区孔隙率ε＝0.23，确定采空区气体动力粘度μ＝1.8×10^-5N·s/m²，确定采空区的渗透系数k＝0.35，确定气体扩散系数D_i＝1.5×10^-5m²/s,确定气体的比热C_pg＝1005J/(kg·℃)，确定煤的密度ρ_c＝1300kg/m³和比热C_pc＝1000J/(kg·℃)，确定气体导热系数λ_g＝0.026W/m·℃，确定煤的导热系数λ_c＝0.20W/m·℃，确定采空区遗煤氧化反应产生的热量Q＝29.5kJ/kg，确定采空区遗煤氧化速率r＝9.94×10^-6kg/s，确定工作面温度T₀＝19.75℃，确定工作面的热流密度q₀＝9.5W/m²，确定采空区边界温度T₁、T₂、T₃(℃)，参见下表：

确定采空区第一边界1、第二边界2、第三边界3的对流换热系数h₁＝h₂＝h₃＝100W/m²·℃，由进风巷道口处温度传感器测取x方向的气体流速u＝2.25m/s，由回风巷道口处温度传感器10测取y方向的气体流速v＝1.85m/s，特点在于：步骤②建立以下方程获取考虑大气压影响的矿井采空区温度场仿真使用的进风巷道口处的全压值p_in(kPa)，回风巷道口处的全压值p_out(kPa)：

p_in＝s_n1·p_n1+s_m1·p_m1

p_out＝s_n2·p_n2+s_m2·p_m2

首先分时段实测42组进风巷道口6处的绝对静压值p_n1(kPa)，即进风巷道口6处大气压力值，回风巷道口7处的绝对静压值p_n2(kPa)，即回巷道口7处大气压力值，进风巷道口6处机械通风压力值p_m1(kPa)，即为不考虑自然风压影响下回巷道口6处的矿井主要通风机克服工作面通风阻力值，回风巷道口7处机械通风压力值p_m2(kPa)，即为不考虑自然风压影响下回巷道口7处的矿井主要通风机克服工作面通风阻力值，参见下表：

根据所测以上数据，采用最小二乘法对上述两式中的s_n1、s_m1、s_n2、s_m2大气压影响因子进行参数辨识，得到大气压影响因子s_n1＝0.9591、s_m1＝2.1798、s_n2＝0.9485、s_m2＝2.8041值后，在测取当地地面大气压为99.225kPa同时，测取进风巷道口6处的绝对静压值p_n1＝101.85kPa，回风巷道口7处的绝对静压值p_n2＝101.53kPa，进风巷道口6处机械通风压力值p_m1＝0.66kPa，回风巷道口7处机械通风压力值p_m2＝0.79kPa后，代入上述两式中计算出考虑大气压影响的矿井采空区温度场仿真使用的进风巷道口处的全压值p_in＝99.11kPa，回风巷道口处的全压值p_out＝98.50kPa。

其中：步骤③

建立平面温度场模型方程：

连续性方程为：

$\frac{\∂u}{\∂x}+\frac{\∂v}{\∂y}=0$

式中，u,v分别为x,y方向的气体流速，m/s，

运动方程为：

$\frac{\∂\left({\mathrm{\ρ}}_{g}u\right)}{\∂t}+{\mathrm{\ρ}}_g(u\frac{\∂u}{\∂x}+v\frac{\∂u}{\∂y})=-\frac{\∂p}{\∂x}+\mathrm{\μ}(\frac{{\∂}^{2}u}{{\∂x}^2}+\frac{{\∂}^{2}u}{{\∂y}^2})-\frac{\mathrm{\μ}}{\mathrm{k\ϵ}}u$

$\frac{\∂\left({\mathrm{\ρ}}_{g}v\right)}{\∂t}+{\mathrm{\ρ}}_g(u\frac{\∂v}{\∂x}+v\frac{\∂v}{\∂y})=-\frac{\∂p}{\∂y}+\mathrm{\μ}(\frac{{\∂}^{2}v}{{\∂x}^2}+\frac{{\∂}^{2}v}{{\∂y}^2})-\frac{\mathrm{\μ}}{\mathrm{k\ϵ}}v$

式中，ρ_g为气体的密度，kg/m³，ε为采空区孔隙率，无量纲，p为气体静压，kPa，t时间，s，μ为采空区气体动力粘度，N·s/m²，k为采空区的渗透系数，无量纲，

组分方程为：

${\mathrm{\ρ}}_{g}u\frac{{\∂Y}_i}{\∂x}+{\mathrm{\ρ}}_{g}v\frac{{\∂Y}_i}{\∂y}={\mathrm{\ρ}}_g{D}_i\frac{{\∂}^2{Y}_i}{{\∂x}^2}+{\mathrm{\ρ}}_g{D}_i\frac{{\∂}^2{Y}_i}{{\∂y}^2}+S_i$

式中，D_i是气体扩散系数，m²/s，Y_i为气体成分浓度，mol/m³，S_i为耗氧速率，mol/m³·s，以温度T为变量的能量为：

$({\mathrm{\ϵ\ρ}}_g{C}_{\mathrm{pg}}+(1-\mathrm{\ϵ}){\mathrm{\ρ}}_c{C}_{\mathrm{pc}})\frac{\∂T}{\∂t}+{\mathrm{\ρ}}_g{C}_{\mathrm{pg}}(u\frac{\∂T}{\∂x}+v\frac{\∂T}{\∂y})={\mathrm{\λ}}_{\mathrm{eff}}(\frac{{\∂}^{2}T}{{\∂x}^2}+\frac{{\∂}^{2}T}{{\∂y}^2})+\mathrm{rQ}$

式中，C_pg为气体的比热，J/(kg·℃)，ρ_c,C_pc分别为煤的密度和比热，kg/m³，J/(kg·℃)，λ_eff为有效导热系数，W/m·℃，Q为采空区遗煤氧化反应产生的热量，kJ/kg，r为采空区遗煤氧化速率，kg/s,λ_eff的计算公式为:

λ_eff＝ελ_g+(1-ε)λ_c

式中，λ_g为气体导热系数，W/m·℃，λ_c为煤的导热系数，W/m·℃，

设定初始条件：

$\left\{\begin{array}{c}u=v=0,p=p_0,T=T_0\\ 0\≤x\≤L,0\≤y\≤W,t=0\end{array}\right.$

式中，p₀为采空区初始气压，kPa，T₀为初始温度，℃，u,v分别为坐标系中x,y方向的气体流速，m/s，L为采空区长度，m，W为采空区宽度，m，t为时间，s，

平面温度场模型风压边界条件方程：

工作面：p|_{x＝0,0≤y≤a}＝p_in，p|_{x＝0,W-a≤y≤W}＝p_out

p|_{x＝0,a≤y≤W-a}＝p_in-(p_in-p_out)·(y-a)/(W-2a)

第一边界： $\frac{\∂p}{\∂y}{|}_{y=\mathrm{0,0}\≤x\≤L}=0$

第二边界： $\frac{\∂p}{\∂x}{|}_{x=L,0\≤y\≤W}=0$

第三边界： $\frac{\∂p}{\∂y}{|}_{y=W,0\≤x\≤L}=0$

式中，p为气体静压，kPa，p_in为进风巷道口处的全压值，kPa，p_out为回风巷道口处的全压值，kPa，L为采空区长度，m，W为采空区宽度，m，a为进风巷道口入风口宽度，m，

平面温度场模型热力学边界条件方程：

工作面：T|_{x＝0,0≤y≤W}＝T₀， $-\frac{\∂T}{\∂x}{|}_{x=\mathrm{0,0}\≤y\≤W}=q_0$

第一边界：T|_{y＝0,0≤x≤L}＝T₁(x)， $-\frac{\∂T}{\∂y}{|}_{y=\mathrm{0,0}\≤x\≤L}=h_1(T_1-T_0)$

第二边界：T|_{x＝L,0≤y≤W}＝T₂(y)， $-\frac{\∂T}{\∂x}{|}_{x=L,0\≤y\≤W}=h_2(T_2-T_0)$

第三边界：T|_{y＝W,0≤x≤L}＝T₃(x)， $-\frac{\∂T}{\∂y}{|}_{y=W,0\≤x\≤L}=h_3(T_3-T_0)$

式中，T₀为工作面温度，℃，q₀为工作面的热流密度，W/m²，T₁、T₂、T₃分别为采空区第一边界1、第二边界2、第三边界3的边界温度，℃，h₁、h₂、h₃分别为采空区第一边界1、第二边界2、第三边界3的边界的对流换热系数，W/m²·℃，L为采空区长度，m，W为采空区宽度，m，

将以上平面温度场模型方程联立并加入风压边界条件及热力学边界条件方程，在计算流体动力学软件FLUENT中进行数值求解：过程为，首先通过File输入.msh文件Check按钮对整个网格的质量进行测试保证计算的收敛，接着在Fluent中选中Energy和porous选项来对能量守恒方程和物质输送守恒方程求解，然后通过viscous按钮选择考虑了湍流的SST模型，工作面边界条件通过boundary condition界面中的presure-inlet/outlet实现，第一、二、三边界通过界面中的symmetry来实现，物理特性参数则通过materials选项输入各对应参数，接着设置Fluent中Initialize按钮来对计算区域各计算变量初始化，设置Monitors按钮对计算过程进行监视控制，最后进行后处理，通过Display-contours按钮输出采空区温度云图。获得地面大气压为99.225kPa时采空区温度场仿真云图，如图3。

其中：采空区边界温度T₁、T₂、T₃(℃)，通过U型设置于边界上的光纤温度传感器测量获得，每隔100m采集一个温度值。

实施例二

参见图1、图2、图4，一种考虑大气压影响的矿井采空区温度场仿真方法，包括如下步骤：步骤①确定采空区边界，包括第一边界1、第二边界2和第三边界3构成的矩形目标采空区，确定矩形目标采空区长度L＝500m值和宽度W＝300m值，确定坐标轴x(m)、y(m)，通过点、线、面的顺序形成整个计算区域，建立采空区的计算域尺寸为500m*300m，并确定各离散单元的界面长度，即网格单元为0.5m，离散采用结构化网格，通过Gambit中的map网格划分方法生成了如图2所示的网格，网格划分成功后再由Gambit输出mesh文件，以便下一步用于Fluent软件求解器的识别和计算，量取进风巷道口6入风口宽度a＝5m，测取气体的密度ρ_g＝1.205kg/m³,确定采空区孔隙率ε＝0.23，确定采空区气体动力粘度μ＝1.8×10^-5N·s/m²，确定采空区的渗透系数k＝0.35，确定气体扩散系数D_i＝1.5×10^-5m²/s,确定气体的比热C_pg＝1005J/(kg·℃)，确定煤的密度ρ_c＝1300kg/m³和比热C_pc＝1000J/(kg·℃)，确定气体导热系数λ_g＝0.026W/m·℃，确定煤的导热系数λ_c＝0.20W/m·℃，确定采空区遗煤氧化反应产生的热量Q＝29.5kJ/kg，确定采空区遗煤氧化速率r＝9.94×10^-6kg/s，确定工作面温度T₀＝20.61℃，确定工作面的热流密度q₀＝9.5W/m²，确定采空区边界温度T₁、T₂、T₃(℃)，参见下表：

确定采空区第一边界1、第二边界2、第三边界3的对流换热系数h₁＝h₂＝h₃＝100W/m²·℃，由进风巷道口处温度传感器测取x方向的气体流速u＝2.17m/s，由回风巷道口处温度传感器10测取y方向的气体流速v＝1.79m/s，特点在于：步骤②建立以下方程获取考虑大气压影响的矿井采空区温度场仿真使用的进风巷道口处的全压值p_in(kPa)，回风巷道口处的全压值p_out(kPa)：

p_in＝s_n1·p_n1+s_m1·p_m1

p_out＝s_n2·p_n2+s_m2·p_m2

首先分时段实测42组进风巷道口6处的绝对静压值p_n1(kPa)，即进风巷道口6处大气压力值，回风巷道口7处的绝对静压值p_n2(kPa)，即回巷道口7处大气压力值，进风巷道口6处机械通风压力值p_m1(kPa)，即为不考虑自然风压影响下回巷道口6处的矿井主要通风机克服工作面通风阻力值，回风巷道口7处机械通风压力值p_m2(kPa)，即为不考虑自然风压影响下回巷道口7处的矿井主要通风机克服工作面通风阻力值，参见下表：

根据所测以上数据，采用最小二乘法对上述两式中的s_n1、s_m1、s_n2、s_m2大气压影响因子进行参数辨识，得到大气压影响因子s_n1＝0.9591、s_m1＝2.1798、s_n2＝0.9485、s_m2＝2.8041值后，在测取当地地面大气压为98.426kPa同时，测取进风巷道口6处的绝对静压值p_n1＝100.96kPa，回风巷道口7处的绝对静压值p_n2＝100.64kPa，进风巷道口6处机械通风压力值p_m1＝0.61kPa，回风巷道口7处机械通风压力值p_m2＝0.70kPa后，代入上述两式中计算出考虑大气压影响的矿井采空区温度场仿真使用的进风巷道口处的全压值p_in＝98.16kPa，回风巷道口处的全压值p_out＝97.42kPa。

其中：步骤③

建立平面温度场模型方程：

连续性方程为：

$\frac{\∂u}{\∂x}+\frac{\∂v}{\∂y}=0$

式中，u,v分别为x,y方向的气体流速，m/s，

运动方程为：

$\frac{\∂\left({\mathrm{\ρ}}_{g}u\right)}{\∂t}+{\mathrm{\ρ}}_g(u\frac{\∂u}{\∂x}+v\frac{\∂u}{\∂y})=-\frac{\∂p}{\∂x}+\mathrm{\μ}(\frac{{\∂}^{2}u}{{\∂x}^2}+\frac{{\∂}^{2}u}{{\∂y}^2})-\frac{\mathrm{\μ}}{\mathrm{k\ϵ}}u$

$\frac{\∂\left({\mathrm{\ρ}}_{g}v\right)}{\∂t}+{\mathrm{\ρ}}_g(u\frac{\∂v}{\∂x}+v\frac{\∂v}{\∂y})=-\frac{\∂p}{\∂y}+\mathrm{\μ}(\frac{{\∂}^{2}v}{{\∂x}^2}+\frac{{\∂}^{2}v}{{\∂y}^2})-\frac{\mathrm{\μ}}{\mathrm{k\ϵ}}v$

式中，ρ_g为气体的密度，kg/m³，ε为采空区孔隙率，无量纲，p为气体静压，kPa，t时间，s，μ为采空区气体动力粘度，N·s/m²，k为采空区的渗透系数，无量纲，

组分方程为：

${\mathrm{\ρ}}_{g}u\frac{{\∂Y}_i}{\∂x}+{\mathrm{\ρ}}_{g}v\frac{{\∂Y}_i}{\∂y}={\mathrm{\ρ}}_g{D}_i\frac{{\∂}^2{Y}_i}{{\∂x}^2}+{\mathrm{\ρ}}_g{D}_i\frac{{\∂}^2{Y}_i}{{\∂y}^2}+S_i$

式中，D_i是气体扩散系数，m²/s，Y_i为气体成分浓度，mol/m³，S_i为耗氧速率，mol/m³·s，以温度T为变量的能量为：

$({\mathrm{\ϵ\ρ}}_g{C}_{\mathrm{pg}}+(1-\mathrm{\ϵ}){\mathrm{\ρ}}_c{C}_{\mathrm{pc}})\frac{\∂T}{\∂t}+{\mathrm{\ρ}}_g{C}_{\mathrm{pg}}(u\frac{\∂T}{\∂x}+v\frac{\∂T}{\∂y})={\mathrm{\λ}}_{\mathrm{eff}}(\frac{{\∂}^{2}T}{{\∂x}^2}+\frac{{\∂}^{2}T}{{\∂y}^2})+\mathrm{rQ}$

式中，C_pg为气体的比热，J/(kg·℃)，ρ_c,C_pc分别为煤的密度和比热，kg/m³，J/(kg·℃)，λ_eff为有效导热系数，W/m·℃，Q为采空区遗煤氧化反应产生的热量，kJ/kg，r为采空区遗煤氧化速率，kg/s,λ_eff的计算公式为:

λ_eff＝ελ_g+(1-ε)λ_c

式中，λ_g为气体导热系数，W/m·℃，λ_c为煤的导热系数，W/m·℃，

设定初始条件：

$\left\{\begin{array}{c}u=v=0,p=p_0,T=T_0\\ 0\≤x\≤L,0\≤y\≤W,t=0\end{array}\right.$

式中，p₀为采空区初始气压，kPa，T₀为初始温度，℃，u,v分别为坐标系中x,y方向的气体流速，m/s，L为采空区长度，m，W为采空区宽度，m，t为时间，s，

平面温度场模型风压边界条件方程：

工作面：p|_{x＝0,0≤y≤a}＝p_in，p|_{x＝0,W-a≤y≤W}＝p_out

p|_{x＝0,a≤y≤W-a}＝p_in-(p_in-p_out)·(y-a)/(W-2a)

第一边界： $\frac{\∂p}{\∂y}{|}_{y=\mathrm{0,0}\≤x\≤L}=0$

第二边界： $\frac{\∂p}{\∂x}{|}_{x=L,0\≤y\≤W}=0$

第三边界： $\frac{\∂p}{\∂y}{|}_{y=W,0\≤x\≤L}=0$

式中，p为气体静压，kPa，p_in为进风巷道口处的全压值，kPa，p_out为回风巷道口处的全压值，kPa，L为采空区长度，m，W为采空区宽度，m，a为进风巷道口入风口宽度，m，

平面温度场模型热力学边界条件方程：

工作面：T|_{x＝0,0≤y≤W}＝T₀， $-\frac{\∂T}{\∂x}{|}_{x=\mathrm{0,0}\≤y\≤W}=q_0$

第一边界：T|_{y＝0,0≤x≤L}＝T₁(x)， $-\frac{\∂T}{\∂y}{|}_{y=\mathrm{0,0}\≤x\≤L}=h_1(T_1-T_0)$

第二边界：T|_{x＝L,0≤y≤W}＝T₂(y)， $-\frac{\∂T}{\∂x}{|}_{x=L,0\≤y\≤W}=h_2(T_2-T_0)$

第三边界：T|_{y＝W,0≤x≤L}＝T₃(x)， $-\frac{\∂T}{\∂y}{|}_{y=W,0\≤x\≤L}=h_3(T_3-T_0)$

式中，T₀为工作面温度，℃，q₀为工作面的热流密度，W/m²，T₁、T₂、T₃分别为采空区第一边界1、第二边界2、第三边界3的边界温度，℃，h₁、h₂、h₃分别为采空区第一边界1、第二边界2、第三边界3的边界的对流换热系数，W/m²·℃，L为采空区长度，m，W为采空区宽度，m，

将以上平面温度场模型方程联立并加入风压边界条件及热力学边界条件方程，在计算流体动力学软件FLUENT中进行数值求解：过程为，首先通过File输入.msh文件Check按钮对整个网格的质量进行测试保证计算的收敛，接着在Fluent中选中Energy和porous选项来对能量守恒方程和物质输送守恒方程求解，然后通过viscous按钮选择考虑了湍流的SST模型，工作面边界条件通过boundary condition界面中的presure-inlet/outlet实现，第一、二、三边界通过界面中的symmetry来实现，物理特性参数则通过materials选项输入各对应参数，接着设置Fluent中Initialize按钮来对计算区域各计算变量初始化，设置Monitors按钮对计算过程进行监视控制，最后进行后处理，通过Display-contours按钮输出采空区温度云图。获得地面大气压为98.426kPa时采空区温度场仿真云图，如图4。

其中：采空区边界温度T₁、T₂、T₃(℃)，通过U型设置于边界上的光纤温度传感器测量获得，每隔100m采集一个温度值。

从图中可见：大气压增大时，采空区4内部气体的绝对压力升高并引起流速降低，进而因多场耦合导致热对流及传导效率低，温度分布的梯度更小。

Claims (4)

1.一种考虑大气压影响的矿井采空区温度场仿真方法，包括如下步骤：步骤①确定采空区边界，包括第一边界、第二边界和第三边界构成的矩形目标采空区，确定矩形目标采空区长度L值和宽度W值，确定坐标轴x、y，确定矩形目标采空区网格单元值，量取进风巷道口入风口宽度a值，测取气体的密度ρ_g值,确定采空区孔隙率ε值，确定采空区物体动力粘度μ值，确定采空区的渗透系数k值，确定气体扩散系数D_i,确定气体的比热C_pg值，确定煤的密度ρ_c和比热C_pc值，确定气体导热系数λ_g，确定煤的导热系数λ_c，确定采空区遗煤氧化反应产生的热量Q，确定采空区遗煤氧化速率r，确定工作面温度T₀值，确定工作面的热流密度q₀值，确定采空区边界温度T₁、T₂、T₃值，确定采空区边界的对流换热系数h₁、h₂、h₃值，测取x,y方向的气体流速u,v值，特征在于：步骤②建立以下方程获取考虑大气压影响的矿井采空区温度场仿真使用的进风巷道口处的全压值p_in，回风巷道口处的全压值p_out：

p_in＝s_n1·p_n1+s_m1·p_m1

p_out＝s_n2·p_n2+s_m2·p_m2

首先分时段实测多组进风巷道口处的绝对静压值p_n1，回风巷道口处的绝对静压值p_n2，进风巷道口处机械通风压力值p_m1，回风巷道口处机械通风压力值p_m2后，采用最小二乘法对上述两式中的s_n1、s_m1、s_n2、s_m2大气压影响因子进行参数辨识，得到大气压影响因子s_n1、s_m1、s_n2、s_m2值后，在测取当地大气压下同时，测取一组进风巷道口处的绝对静压值p_n1，回风巷道口处的绝对静压值p_n2，进风巷道口处机械通风压力值p_m1，回风巷道口处机械通风压力值p_m2后，代入上述两式中计算出考虑大气压影响的矿井采空区温度场仿真使用的进风巷道口处的全压值p_in，回风巷道口处的全压值p_out。

2.根据权利要求1所述的一种考虑大气压影响的矿井采空区温度场仿真方法，其特征在于：步骤③

建立平面温度场模型方程：

连续性方程为：

$\frac{\∂u}{\∂x}+\frac{\∂v}{\∂y}=0$

式中，u,v分别为x,y方向的气体流速，m/s，

运动方程为：

$\frac{\∂\left({\mathrm{\ρ}}_{g}u\right)}{\∂t}+{\mathrm{\ρ}}_g(u\frac{\∂u}{\∂x}+v\frac{\∂u}{\∂y})=-\frac{\∂p}{\∂x}+\mathrm{\μ}(\frac{{\∂}^{2}u}{\∂x^2}+\frac{{\∂}^{2}u}{\∂y^2})-\frac{\mathrm{\μ}}{\mathrm{k\ϵ}}u$

$\frac{\∂\left({\mathrm{\ρ}}_{g}v\right)}{\∂t}+{\mathrm{\ρ}}_g(u\frac{\∂v}{\∂x}+v\frac{\∂v}{\∂y})=-\frac{\∂p}{\∂y}+\mathrm{\μ}(\frac{{\∂}^{2}v}{\∂x^2}+\frac{{\∂}^{2}v}{\∂y^2})-\frac{\mathrm{\μ}}{\mathrm{k\ϵ}}v$

式中，ρ_g为气体的密度，kg/m³，ε为采空区孔隙率，无量纲，p为气体静压，kPa，t时间，s，μ为采空区气体动力粘度，N·s/m²，k为采空区的渗透系数，无量纲，

组分方程为：

${\mathrm{\ρ}}_{g}u\frac{\∂Y_i}{\∂x}+{\mathrm{\ρ}}_{g}v\frac{\∂Y_i}{\∂y}={\mathrm{\ρ}}_g{D}_i\frac{{\∂}^2{Y}_i}{\∂x^2}+{\mathrm{\ρ}}_g{D}_i\frac{{\∂}^2{Y}_i}{\∂y^2}+S_i$

式中，D_i是气体扩散系数，m²/s，Y_i为气体成分浓度，mol/m³，S_i为耗氧速率，mol/m³·s，

以温度T为变量的能量为：

$(\mathrm{\ϵ}{\mathrm{\ρ}}_g{C}_{\mathrm{pg}}+(1-\mathrm{\ϵ}){\mathrm{\ρ}}_c{C}_{\mathrm{pc}})\frac{\∂T}{\∂t}+{\mathrm{\ρ}}_g{C}_{\mathrm{pg}}(u\frac{\∂T}{\∂x}+v\frac{\∂T}{\∂y})={\mathrm{\λ}}_{\mathrm{eff}}(\frac{{\∂}^{2}T}{\∂x^2}+\frac{{\∂}^{2}T}{\∂y^2})+\mathrm{rQ}$

式中，C_pg为气体的比热，J/(kg·℃)，ρ_c,C_pc分别为煤的密度和比热，kg/m³，J/(kg·℃)，λ_eff为有效导热系数，W/m·℃，Q为采空区遗煤氧化反应产生的热量，kJ/kg，r为采空区遗煤氧化速率，kg/s,λ_eff的计算公式为:

λ_eff＝ελ_g+(1-ε)λ_c

式中，λ_g为气体导热系数，W/m·℃，λ_c为煤的导热系数，W/m·℃，

设定初始条件：

$\left\{\begin{array}{c}u=v=0,p=p_0,T=T_0\\ 0\≤x\≤L,0\≤y\≤W,t=0\end{array}\right.$

式中，p₀为采空区初始气压，kPa，T₀为初始温度，℃，u,v分别为坐标系中x,y方向的气体流速，m/s，L为采空区长度，m，W为采空区宽度，m，t为时间，s，

平面温度场模型风压边界条件方程：

工作面：p|_{x＝0,0≤y≤a}＝p_in，p|_{x＝0,W-a≤y≤W}＝p_out

p|_{x＝0,a≤y≤W-a}＝p_in-(p_in-p_out)·(y-a)/(W-2a)

第一边界： $\frac{\∂p}{\∂y}{|}_{y=\mathrm{0,0}\≤x\≤L}=0$

第二边界： $\frac{\∂p}{\∂x}{|}_{x=L,0\≤y\≤W}=0$

第三边界： $\frac{\∂p}{\∂y}{|}_{y=W,0\≤x\≤L}=0$

式中，p为气体静压，kPa，p_in为进风巷道口处的全压值，kPa，p_out为回风巷道口处的全压值，kPa，L为采空区长度，m，W为采空区宽度，m，a为进风巷道口入风口宽度，m，

平面温度场模型热力学边界条件方程：

工作面： ${T|}_{x=\mathrm{0,0}\≤y\≤W}=T_0,-\frac{\∂T}{\∂x}{|}_{x=\mathrm{0,0}\≤y\≤W}=q_0$

第一边界： ${T|}_{y=\mathrm{0,0}\≤x\≤L}=T_1\left(x\right),-\frac{\∂T}{\∂y}{|}_{y=\mathrm{0,0}\≤x\≤L}=h_1(T_1-T_0)$

第二边界： ${T|}_{x=L,0\≤y\≤W}=T_2\left(y\right),-\frac{\∂T}{\∂x}{|}_{x=L,0\≤y\≤W}=h_2(T_2-T_0)$

第三边界： ${T|}_{y=W,0\≤x\≤L}=T_3\left(x\right),-\frac{\∂T}{\∂y}{|}_{y=W,0\≤x\≤L}=h_3(T_3-T_0)$

式中，T₀为工作面温度，℃，q₀为工作面的热流密度，W/m²，T₁、T₂、T₃分别为采空区第一边界1、第二边界2、第三边界3的边界温度，℃，h₁、h₂、h₃分别为采空区第一边界1、第二边界2、第三边界3的边界的对流换热系数，W/m²·℃，L为采空区长度，m，W为采空区宽度，m，

将以上平面温度场模型方程联立并加入风压边界条件及热力学边界条件方程，在计算流体动力学软件FLUENT中进行数值求解，获得在一定地面大气压下采空区温度场仿真云图。

3.根据权利要求1所述的一种考虑大气压影响的矿井采空区温度场仿真方法，其特征在于：采空区边界温度T₁、T₂、T₃，通过U型设置于边界上的光纤温度传感器测量获得，每隔100m采集一个温度值。

4.根据权利要求2所述的一种考虑大气压影响的矿井采空区温度场仿真方法，其特征在于：采空区边界温度T₁、T₂、T₃，通过U型设置于边界上的光纤温度传感器测量获得，每隔100m采集一个温度值。