楊 偉， 薛思瀚， 劉秦見， 張樹光

（1.遼寧工程技術大學建筑工程學院，阜新 123000；2.遼寧工程技術大學土木與交通學院，阜新 123000）

沉降分布孔隙率多孔介質模型及數(shù)值研究

楊 偉1， 薛思瀚1， 劉秦見1， 張樹光2

（1.遼寧工程技術大學建筑工程學院，阜新 123000；2.遼寧工程技術大學土木與交通學院，阜新 123000）

研究沉降分布孔隙率多孔介質流動和傳熱，根據(jù)“O”形圈理論和現(xiàn)場測定確定孔隙率系數(shù)，建立坐標方向孔隙率分布函數(shù)；考慮流體密度變化，并引入Brinkman-Forchheimer的擴展Darcy模型，能量方程采用界面連續(xù)條件，建立沉降分布孔隙率多孔介質流動和傳熱求解模型.采用差分法對模型進行離散化，應用高斯–賽德爾方法迭代求解.數(shù)值分析表明：沉降分布孔隙率條件下多孔介質內流體流動速度在壁面附近較大，中心部位較小，壁面附近孔隙率的增大使得低流速區(qū)域減小，較高流速區(qū)域增大；當孔隙率小值時，溫度按線性減小；當孔隙率大值時，溫度在高低溫壁面附近迅速減小，在中部減小較緩，熱量按導熱和對流共同傳遞；孔隙率增大能使平均怒謝爾數(shù)增大，對流換熱作用增強.

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

0 引言

沉降是一種自然現(xiàn)象，多孔介質沉降使孔隙率發(fā)生變化，形成二次分布.由于煤礦地下開采擾動原來地質結構，很容易發(fā)生煤炭、巖石的沉降，形成孔隙率的二次分布，二次分布的結果是邊界部位孔隙率增大，引起平均孔隙率的增大，導致對流換熱作用增強.錢鳴高［1］對覆巖采動裂隙進行研究，發(fā)現(xiàn)孔隙率二次分布的“O”形圈特征.對于含多孔介質孔隙率均勻分布方腔內自然對流研究較多［2－3］，而關于沉降導致孔隙率二次分布影響的多孔介質內自然對流的研究卻較少［4－5］，研究結果表明孔隙率二次分布后多孔介質內的風流發(fā)生較大變化［6］；同樣，孔隙率發(fā)生變化后多孔介質流場、溫度場要發(fā)生一定的變化.本文在“O”型圈理論基礎上，根據(jù)實測數(shù)據(jù)建立孔隙率分布函數(shù)模型，考慮密度隨溫度變化等因素，用Brinkman-Forchheimer的擴展Darcy模型建立孔隙率二次分布后的多孔介質流場、溫度場求解模型.

1 物理問題及數(shù)學模型

煤礦地下開采后上覆巖層發(fā)生向下沉陷，沉陷區(qū)域發(fā)生二次采動孔隙率分布［7］，孔隙率圖像分析如圖1.“O”型圈理論認為在沉陷區(qū)域中心部分基本被壓實，孔隙率維持在一平均水平，四周存在一連通的離層裂隙發(fā)育區(qū)，孔隙率發(fā)生較大變化.

連續(xù)性方程考慮了密度隨溫度變化［5］，ρ＝β－1（?ρ／?T）變化，動量方程用Brinkman-Forchheimer的擴展Darcy模型表示由于微孔結構引起的體積力，能量方程采用界面連續(xù)條件處理，即流體區(qū)域和多孔介質接觸壁面溫度與熱流密度相等.則Ts＝Tf＝T.二維數(shù)學模型可表示為：連續(xù)性方程動量方程

上述方程中，ρ為密度，ρ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為導熱系數(shù)，W·m－1·K－1，下標s為固體，下標f為流體；K為多孔介質滲透率，K＝（1－ε）－2，m2，db為顆粒直徑，系數(shù)C，Ergun［8］給出C＝150，Bear［9］給出C＝180，Kozeny-earmen［10］給出C＝172.8，本文取C＝175［11］；引入瑞利數(shù)Ra＝gβ（Th－Tc）H3ν－1a－1，β為熱膨脹系數(shù)，K－1，β＝－ρ－1?p／?T；a為熱擴散率，m2·s－1.

孔隙率是指材料中連通孔隙體積和材料總體積之比.沉陷區(qū)多孔介質的孔隙率與空間位置有關，在x軸方向上，壁面兩端的孔隙率較工作面中部大.參照“O”形圈理論的表達方法用平均孔隙率εx和孔隙率修正系數(shù)ε′乘積來表達沉陷分布的實際孔隙率（式（5））.根據(jù)現(xiàn)場實測結果，將測定數(shù)據(jù)進行歸納為平均值和修正系數(shù)，將修正系數(shù)擬合為式（6）和圖2.

式中，εx為坐標方向平均孔隙率；x值正負選取，x≥H／2為正值，x＜H／2為負值.

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

2 模型驗證與算例分析

根據(jù)沉陷特征，取中心線平面建立二維物理模型如圖3所示，研究區(qū)域選定為長寬為H×H的二維平面.左壁面溫度為高溫Th、右壁面溫度為低溫Tc，上下壁面為絕熱邊界條件，多孔介質骨架與內部流體溫度相等，即采用局部熱平衡假設；流體密度隨溫度變化，其他熱物性為常量，在模型壁面上速度均采用無滑移邊界條件.內部充滿流體（水）.

圖4為速度v在y＝0.5 m處的值.從圖4可明顯看出：速度v在兩側大，在－0.1～0.1范圍內很?。紫堵示鶆蚍植紩r在（0.2～0.8）很小范圍內［13］）；當εx＝0.1，瑞利數(shù)Ra＝100時，v在整個范圍內很小；當Ra相同時，εx值大，在影響區(qū)（－0.5～－0.3，0.3～0.5）范圍流速增加值大；當εx相同時，Ra值大，在影響區(qū)范圍流速增加值大；說明εx和Ra同樣對多孔介質內流體速度起到改變作用.

圖5為溫度在y＝0.5 m處的值.從圖5可明顯看出：溫度沿長度方向逐漸減小，最大值在左側高溫壁面為305 K，最小值在右側低溫壁面為295 K；當εx小時，溫度按線性減小，熱量按導熱傳遞；當εx大值，溫度在高低溫壁面附近迅速減小，在中部減小較緩，熱量按導熱和對流共同傳遞.圖5佐證圖4速度計算結果.

圖6為高溫壁面平均怒謝爾數(shù)Nu 隨瑞利數(shù)Ra變化.從圖4可明顯看出：孔隙率和瑞利數(shù)對腔內對流換熱過程的影響是顯著的.隨著εx、Ra逐漸增加，Nu 增加，增加結果在同一數(shù)量級，多孔介質內換熱強度隨對流作用增加而增強.

3 結論

1）根據(jù)“O”形圈理論和現(xiàn)場測定，建立x方向孔隙率分布函數(shù)；考慮流體密度變化，用Brinkman-Forchheimer的擴展Darcy模型建立沉降分布孔隙率多孔介質求解模型；

2）對沉降分布孔隙率多孔介質求解模型進行數(shù)值計算，多孔介質內流體流動速度在壁面附近較大，中心部位較小，壁面附近孔隙率的增大使得低流速區(qū)域減小，較高流速區(qū)域增大；

3）溫度沿水平長度方向逐漸減小，最大值在左側高溫壁面，最小值在右側低溫壁面；當孔隙率小值時，溫度按線性減小；當孔隙率大值時，溫度在高低溫壁面附近迅速減小，在中部減小較緩，熱量按導熱和對流共同傳遞；

4）沉降分布孔隙率條件下，邊界部位孔隙率增大引起平均孔隙率增大，致使平均怒謝爾數(shù)增大，對流換熱作用增強.

5）本文研究煤礦地下開采后上覆巖層發(fā)生向下沉陷，沉陷區(qū)域發(fā)生二次采動孔隙率分布特征所引起的方腔內自然對流換熱影響，對于其他幾何形狀沒有研究.在受限空間里，幾何形狀對自然對流換熱有著較大影響，因此，當情況不同時，應實測孔隙率的變化和幾何形狀再進行計算.

［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，F(xiàn)uxin 123000，Liaoning，China；2.School of Civil Engineering and Transportation，Liaoning Technical University，F(xiàn)uxin 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