杨 伟， 薛思瀚， 刘秦见， 张树光

（1.辽宁工程技术大学建筑工程学院，阜新 123000；2.辽宁工程技术大学土木与交通学院，阜新 123000）

沉降分布孔隙率多孔介质模型及数值研究

杨 伟1， 薛思瀚1， 刘秦见1， 张树光2

（1.辽宁工程技术大学建筑工程学院，阜新 123000；2.辽宁工程技术大学土木与交通学院，阜新 123000）

研究沉降分布孔隙率多孔介质流动和传热，根据“O”形圈理论和现场测定确定孔隙率系数，建立坐标方向孔隙率分布函数；考虑流体密度变化，并引入Brinkman-Forchheimer的扩展Darcy模型，能量方程采用界面连续条件，建立沉降分布孔隙率多孔介质流动和传热求解模型.采用差分法对模型进行离散化，应用高斯–赛德尔方法迭代求解.数值分析表明：沉降分布孔隙率条件下多孔介质内流体流动速度在壁面附近较大，中心部位较小，壁面附近孔隙率的增大使得低流速区域减小，较高流速区域增大；当孔隙率小值时，温度按线性减小；当孔隙率大值时，温度在高低温壁面附近迅速减小，在中部减小较缓，热量按导热和对流共同传递；孔隙率增大能使平均怒谢尔数增大，对流换热作用增强.

沉降分布孔隙率；Darcy模型；多孔介质；“O”型理论

0 引言

沉降是一种自然现象，多孔介质沉降使孔隙率发生变化，形成二次分布.由于煤矿地下开采扰动原来地质结构，很容易发生煤炭、岩石的沉降，形成孔隙率的二次分布，二次分布的结果是边界部位孔隙率增大，引起平均孔隙率的增大，导致对流换热作用增强.钱鸣高［1］对覆岩采动裂隙进行研究，发现孔隙率二次分布的“O”形圈特征.对于含多孔介质孔隙率均匀分布方腔内自然对流研究较多［2－3］，而关于沉降导致孔隙率二次分布影响的多孔介质内自然对流的研究却较少［4－5］，研究结果表明孔隙率二次分布后多孔介质内的风流发生较大变化［6］；同样，孔隙率发生变化后多孔介质流场、温度场要发生一定的变化.本文在“O”型圈理论基础上，根据实测数据建立孔隙率分布函数模型，考虑密度随温度变化等因素，用Brinkman-Forchheimer的扩展Darcy模型建立孔隙率二次分布后的多孔介质流场、温度场求解模型.

1 物理问题及数学模型

煤矿地下开采后上覆岩层发生向下沉陷，沉陷区域发生二次采动孔隙率分布［7］，孔隙率图像分析如图1.“O”型圈理论认为在沉陷区域中心部分基本被压实，孔隙率维持在一平均水平，四周存在一连通的离层裂隙发育区，孔隙率发生较大变化.

连续性方程考虑了密度随温度变化［5］，ρ＝β－1（∂ρ／∂T）变化，动量方程用Brinkman-Forchheimer的扩展Darcy模型表示由于微孔结构引起的体积力，能量方程采用界面连续条件处理，即流体区域和多孔介质接触壁面温度与热流密度相等.则Ts＝Tf＝T.二维数学模型可表示为：连续性方程动量方程

上述方程中，ρ为密度，ρref为参考密度，kg·m－3；＾分别为x、y方向的流体流动速度，u＝＾，v＝＾，m·s－1；t为时间，s；ε为孔隙率；μ为有效粘度，Pa·s；g为重力加速度，m·s－2；T为温度，K；c为比热容，J·kg·K；k为导热系数，W·m－1·K－1，下标s为固体，下标f为流体；K为多孔介质渗透率，K＝（1－ε）－2，m2，db为颗粒直径，系数C，Ergun［8］给出C＝150，Bear［9］给出C＝180，Kozeny-earmen［10］给出C＝172.8，本文取C＝175［11］；引入瑞利数Ra＝gβ（Th－Tc）H3ν－1a－1，β为热膨胀系数，K－1，β＝－ρ－1∂p／∂T；a为热扩散率，m2·s－1.

孔隙率是指材料中连通孔隙体积和材料总体积之比.沉陷区多孔介质的孔隙率与空间位置有关，在x轴方向上，壁面两端的孔隙率较工作面中部大.参照“O”形圈理论的表达方法用平均孔隙率εx和孔隙率修正系数ε′乘积来表达沉陷分布的实际孔隙率（式（5））.根据现场实测结果，将测定数据进行归纳为平均值和修正系数，将修正系数拟合为式（6）和图2.

式中，εx为坐标方向平均孔隙率；x值正负选取，x≥H／2为正值，x＜H／2为负值.

式（1）～（6）构成了沉降分布孔隙率的多孔介质模型.

2 模型验证与算例分析

根据沉陷特征，取中心线平面建立二维物理模型如图3所示，研究区域选定为长宽为H×H的二维平面.左壁面温度为高温Th、右壁面温度为低温Tc，上下壁面为绝热边界条件，多孔介质骨架与内部流体温度相等，即采用局部热平衡假设；流体密度随温度变化，其他热物性为常量，在模型壁面上速度均采用无滑移边界条件.内部充满流体（水）.

图4为速度v在y＝0.5 m处的值.从图4可明显看出：速度v在两侧大，在－0.1～0.1范围内很小（孔隙率均匀分布时在（0.2～0.8）很小范围内［13］）；当εx＝0.1，瑞利数Ra＝100时，v在整个范围内很小；当Ra相同时，εx值大，在影响区（－0.5～－0.3，0.3～0.5）范围流速增加值大；当εx相同时，Ra值大，在影响区范围流速增加值大；说明εx和Ra同样对多孔介质内流体速度起到改变作用.

图5为温度在y＝0.5 m处的值.从图5可明显看出：温度沿长度方向逐渐减小，最大值在左侧高温壁面为305 K，最小值在右侧低温壁面为295 K；当εx小时，温度按线性减小，热量按导热传递；当εx大值，温度在高低温壁面附近迅速减小，在中部减小较缓，热量按导热和对流共同传递.图5佐证图4速度计算结果.

图6为高温壁面平均怒谢尔数Nu 随瑞利数Ra变化.从图4可明显看出：孔隙率和瑞利数对腔内对流换热过程的影响是显著的.随着εx、Ra逐渐增加，Nu 增加，增加结果在同一数量级，多孔介质内换热强度随对流作用增加而增强.

3 结论

1）根据“O”形圈理论和现场测定，建立x方向孔隙率分布函数；考虑流体密度变化，用Brinkman-Forchheimer的扩展Darcy模型建立沉降分布孔隙率多孔介质求解模型；

2）对沉降分布孔隙率多孔介质求解模型进行数值计算，多孔介质内流体流动速度在壁面附近较大，中心部位较小，壁面附近孔隙率的增大使得低流速区域减小，较高流速区域增大；

3）温度沿水平长度方向逐渐减小，最大值在左侧高温壁面，最小值在右侧低温壁面；当孔隙率小值时，温度按线性减小；当孔隙率大值时，温度在高低温壁面附近迅速减小，在中部减小较缓，热量按导热和对流共同传递；

4）沉降分布孔隙率条件下，边界部位孔隙率增大引起平均孔隙率增大，致使平均怒谢尔数增大，对流换热作用增强.

5）本文研究煤矿地下开采后上覆岩层发生向下沉陷，沉陷区域发生二次采动孔隙率分布特征所引起的方腔内自然对流换热影响，对于其他几何形状没有研究.在受限空间里，几何形状对自然对流换热有着较大影响，因此，当情况不同时，应实测孔隙率的变化和几何形状再进行计算.

［1］Qian Minggao，Xu Jialin.Study on the‘O-shape'circle distribution characteristics of mining-induced fractures in the overlaying strata［J］.Journal of China Coal Society，1998，23（5）：466－469.

［2］Zeng Min，Wang Gang，Xie Gongnan，et al.Natural convection in a wavy enclosure filled with porous media［J］.Chinese Journal of Computational Physics，2008，25（4）：445－449.

［3］NieldD A，Bejan A.Convection in porous media［M］.4th Ed.New York：Springer，2013.

［4］Rees D A S.The onset of Darcy-Brinkman convection in a porous layer：An asymptotic analysis［J］.International Journal of Heat and Mass Transfer，2002，45：2213－2220.

［5］PrasadV，Kulack F，Keyhani M，et al.Natural convection in porous media［J］.Journal of Fluid Mechanics，1985，150：89－119.

［6］Kim M C.Analysis of onset of buoyancy-driven convection in a fluid layer saturated in anisotropic porous media by the relaxed energy method［J］.Journal of Chemical Engineering，2013，30：1207－1212.

［7］Liang Yuntao，Zhang Tengfei，Wang Shugang，et al.Heterogeneous model of porosity in gobs and its airflow field distribution ［J］.Journal of China Coal Society，2009，34（9）：1203－1207.

［8］Ergun S.Fluid flow through packed columns［M］.Chemical Engineering Progress，1952：89－94.

［9］Bearj.Dynamics of fluids in Porous media［M］.New York：American Elsevier Publishing Company，inc，1972.

［10］Combamous M A，Bories A.Hydrothermal convection in saturated Porous media［M］.Adv Ydrosci，1975：231－301.

［11］Yang Wei，Yang Linlin，Sun Yue.Numerical simulation of air temperature of laneway in Dalong Mine，Liaoning province［J］. The Chinese Journal of Geological Hazard and Control，2010，21（2）：84－88.

［12］Yang Wei，Yang Linlin，An Jing，et al.Numerical simulation and analysis of mine ventilation temperature field of integrated mining area in high temperature mine［J］.Coal Geology＆Exploration，2011，39（5）：55－58.

［13］Capone F，Gentile M，Hill A A.Convection problems in anisotropic porous media with nonhomogeneous porosity and thermal diffusivity［J］.Acta Applicandae Mathematicae，2012，122：85－91.

Model of Porous Media with Settlement Distributed Porosity and Numerical Study

YANG Wei1，XUE Sihan1，LIU Qinjian1，ZHANG Shuguang2
（1.School of Architecture Engineering，Liaoning Technical University，Fuxin 123000，Liaoning，China；2.School of Civil Engineering and Transportation，Liaoning Technical University，Fuxin 123000，Liaoning，China）

Flow and heat transfer in porous media with settlement distributed porosity is studied.Based on‘O-shape'ring theory and in situ measurement，a distribution function of porosity in x direction is established.Considering change of fluid density，Brinkman-Forchheimer-extended Darcy model is introduced.Energy equation is solved with continuity boundary conditions.A model of flow and heat transfer in porous media with settlement distributed porosity is established.The model is discretized with difference method. Gauss-Seidel iteration method is used.Numerical analysis shows that：Fluid flow velocity in porous media with of settlement distributed porosity is greater near wall surface，and it is small at center.Increase of porosity near wall surface reduces low velocity zone and increases high velocity zone.As porosity is small，temperature decreases linearly.As porosity is large，temperature near high and low temperature wall surface reduces rapidly.The decrease becomes slower at center.Heat is transferred through conduction and convection together.Increase of porosity increases average Nusselt number and enhances heat convection.

settlement distributed porosity；Darcy model；porous media；“O-shape”ring theory

date： 2013－06－26；Revised date： 2013－10－17

TK124 O159

A

1001-246X（2014）03-0331-04

2013－06－26；

2013－10－17

国家自然科学基金（50804021）资助项目

杨伟（1965－），男，辽宁阜新，副教授，研究方向为流固耦合传热机理，E-mail：lgdyw＠163.com