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

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

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

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

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

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

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

0 引言

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

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

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

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

上述方程中，ρ為密度，ρref為參考密度，kg·m－3；＾分別為x、y方向的流體流動(dòng)速度，u＝＾，v＝＾，m·s－1；t為時(shí)間，s；ε為孔隙率；μ為有效粘度，Pa·s；g為重力加速度，m·s－2；T為溫度，K；c為比熱容，J·kg·K；k為導(dǎo)熱系數(shù)，W·m－1·K－1，下標(biāo)s為固體，下標(biāo)f為流體；K為多孔介質(zhì)滲透率，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為熱擴(kuò)散率，m2·s－1.

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

式中，εx為坐標(biāo)方向平均孔隙率；x值正負(fù)選取，x≥H／2為正值，x＜H／2為負(fù)值.

式（1）～（6）構(gòu)成了沉降分布孔隙率的多孔介質(zhì)模型.

2 模型驗(yàn)證與算例分析

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

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

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

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

3 結(jié)論

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

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

3）溫度沿水平長(zhǎng)度方向逐漸減小，最大值在左側(cè)高溫壁面，最小值在右側(cè)低溫壁面；當(dāng)孔隙率小值時(shí)，溫度按線性減小；當(dāng)孔隙率大值時(shí)，溫度在高低溫壁面附近迅速減小，在中部減小較緩，熱量按導(dǎo)熱和對(duì)流共同傳遞；

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

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

［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

國(guó)家自然科學(xué)基金（50804021）資助項(xiàng)目

楊偉（1965－），男，遼寧阜新，副教授，研究方向?yàn)榱鞴恬詈蟼鳠釞C(jī)理，E-mail：lgdyw＠163.com