Study of 3-D Numerical Simulation for Gas Transfer in the Goaf of the Coal Mining

Jun. 2007Joural of China University of Mining & TechnologyVol.17 No.2Available online at www.sciencedirect.comSCIENCE*@o)DIRECT.J China Univ Mining & Technol 2007, 17(2): 0152 - 0157.Study of 3-D Numerical Simulation forGas Transfer in the Goaf of the Coal MiningWU Zheng-yan, JIANG Shu-guang, HE Xin-jian, WANG Lan-yun, LIN Bai quanSchool of Mining and Safety Engineering, China University of Mining & Technology, Xuzhou, Jiangsu 221008, ChinaAbstract: In order to simulate field distribution rules, mathematical models for 3-D air flows and gas transfer in thegoaf of the coal mining are established, based on theories of permeability and dynamic dispersion through porous media.A gas dispersion equation in a 3-D ficld is calculated by use of numerical method on a weighted upstream multi-elementbalance. Based on data of an example with a U type ventilation mode, surface charts of air pressure distribution and gasconcentration are drawn by Graphtoo software. Finally, a comparison between actually measured results in the modeltest and the numerical simulation results is made to proves the numerical implementation feasible.Key words: 3D numerical simulation of gas transfer in the goaf; air pressure distribution in the goaf; weighted upstreammulti-element balance numerical simulation methodCLC number: TD 712.5distributions degenerate into 2-D space. But most1 Introductionquestions cannot actually be regarded as 2-D spaceproblems. Therefore, some scholars have made someAt present, studies on gas distribution in the goafefforts in the study of 3-D steady flow seepage andof the coal mining largely rely on sketchy estimationgas concentration distributions', but few examplesby setting some gas measuring points in a goaf or byhave been reported.model test in the laboratory. But these two methodsIn numerical solutions for advection-diffusionare not very useful, which need a lot of time and ef-equations, there are some computational methods,fort. With the development of computer technology,such as the traditional Galerkin finite element methodthe numerical simulation method has become a majorand finite difference method. The Galerkin finiteresearch method for gas transfer.element method is based on the principle of variation,In the aspect of the theoretical analyses, manysuch as principle of minimum potential energy andscholars have gradually come to realize that the es-principle of virtual work. Finite difference methods insence of gas tansfer and gas concentration distribu-the unit network, such as a triangle network, aretion in the goaf is that matter moves and diffuses in abased on an integral conservation law. With advec-porous medium and they can analyze it qualitativelytion as leader, there are bad numerical oscillations, In the early 1960's, Shi-ning Zhou was the firstwhen these methods are applied. S. V. Patankar sug-to apply Darcy's law to the gas flow in the coal seams.gested that there were some falsehood numerical OS-Later air flow in the goaf was regarded as a 2-D,cillations in some numerical methods when the Pecletsteady and incompressible potential energy flow andnumber is comparatively large. Herich and Zienki-a flow field distributions were obtained by analyticalwicz gave a windward numerical solution to elimi-methods. For more complicated geometric boundaries,nate falsehood numerical oscillations. Sun and Yehnumerical examples were obtained by finite element(1983) advanced a weighted upstream method basedmethods. These days empirical and semi empiricalon a multi-element balance to add an upstream weightformulas of gas dispersion and distribution laws arefunction directly to the basis sets, which made thebeing developed5 o0. The above studies are all basedwindw:中国煤化工popular and con-on ideal cases where air flows and gas concentrationvenientReceived 13 June 200; accepted 10 November 2006YHCNMHGProjects 50534090 and 50674090 supported by the National Natural Science Foundation of China and 2005CB221503 by the National Key Basic ResearchDevelopment Program (973 Program)Corresponding author. Tel: +86 516 83885156; E-mnail adress: wzy. 1998@ 163.comWU Zheng-yanetalStudy of 3-D Numerical Simulation for Gas Transfer in the Goaf of the Coal Mining153To simulate more truthfully the laws of air flowsaC。+1a =div(DgradC,)-v . gradC。(2)and gas transfers, we established a 3-D mathematicalmodel of an unstable, compressible and multicom-where C。 is the concentration of the a component,ponent air flow and gas transfer, based on the data ofIa the source or sink of the a component, D the dy-model test in lab by Jiang (1994) and Jiang and Wangnamic dispersion coefficient tensor of flow force, V(1995)911. Our numerical research on the multi-the average flow velocity when air is flowing andcomponent gas transfer model is developed by av=q/n.weighted upstream multi-element balance method.The initial condition is that C。 =Cao(x, y, z).The first class boundary condition is that, if the2 Mathematical Model of Gas Transfer inconcentration in the boundary of the goaf is given,the Goafthen Ca =Can(xy,z,D).The second class boundary condition is that, if theDarcy's classical law is the basic start to studyconcentration dispersion flux of the a componentflows in porous media, and the Fick's diffusion lawin the boundary of the goaf is given， thendeals with the fluid hydrodynamic dispersion.- - DgradCa=f(x, y,Z, [).2.1 Establishment of air flow equation in theThe third class boundary condition is that, if thegoafflux of the a component of air in part of the bounda-Air flow in the goaf obeys approximately Darcy'sries of the goaf is given, then ( DgradC。- vCa )nlinearlaw. Connected with conservation of mass, the= f2(x,y,z,1).Navier- Stokes equation based on the theory of seep-In Eq. (1), the two parameters K and n must to beage is as follows:determined, both of which are functions of space co-pdiv(Kgradp)= nap+W(1)ordinates. K = kpg/u. k is the percolation rate in thegoaf andk = k(x, y, z) is a constant scalar functionwith the assumption of isotropy, which has an affinitywhere φ is the total air head function of the goaf, p airwith the porosity factor n, k = Con'/(1-n)". Cois adensity, K a seepage coefficient of isotropy, K=K(x, y,constant to be determined. The n can be determinedz), n the porosity factor of medium in the goaf, W theby selecting enough observed points to measure airsource qr sink when gas is gushing or exhausting inpressure or air quantity. Thus, unknown parametersthe goaf and P air pressure of the goaf. q is thesuch as n, k and the like are resolved.seepage velocity and q = -Kgradp. After consider-ing a general 3-D pattern and boundary conditions,3 Numerical Implementation of Gas Tran-when gravity is neglected, φ=P, otherwise φ=P+pgzsfer in the Goaf(g is the acceleration of gravity).The initial condition is that P(=0)=Po(x, y, z).Consider the following general 3-D pattern flowThe first class boundary condition is that, if airequation and concentration dispersion equation:pressure in the boundary of the goaf is given, then φ= P1(x, y,z, [).aKoThe second class boundary condition is that, if airdH_“4 axg )= m--W(3)quantity in the boundary of the goaf is given, thenOtaxa-Kgradp = q(x, y, z, 1).The third class boundary condition is that air pres-a DopaCsure is given in parts of the boundary and that the air(4)quantity is given in other parts of the boundary in thedxagoaf.where Kap(a, β=x, y, z) is a permeability coefficient,2.2 Establishment of multicomponent gas contentH the head, Dap a fluid hydrodynamic dispersion co-dispersion equation in the goafefficient, W and M are the source or sink terms whichThe convection and diffuse effect of the multicom-are positive when the air is flowing out and negativeponent gas is mainly considered when air transfer inwhen the air is flowing in. If air flow is steady in thethe goaf is studied. When an infinitesimal cell is se-goaf, then m=1 and s=0. From Eq. (3) and Eq. (4), itlected in the goaf, in δt time, outflow or inflow fromis clear that there is_ an additional convection term中国煤化工，this infinitesimal cell should be equal to the algebraiche comparability ofsum of outlow and inflow of the a component of theMHCNMHGair in the infinitesimal cell, which are substituted intothe two equatons, we discuss the Eq. (4) mainly, be-Fick's law. Then we have the following equation:，cause the results of the Eq. (4) can conveniently beapplied to the Eq. (3). .154Jourmal of China University of Mining & TechnologyVol.17_ No.23.1 3-D weighted upstream multi-element balan-[h+h=1, h+h=1,名+李=1,ce methoda +a32 +a_=1, 0, +ag+0%=1,For calculating the Eq. (4), a new numerical simu-lation method is selected, that is the 3-D weighted0=(h+h+h), 02=2(0+4+%),upstream multi- element balance method. First, theh≥0,∞≥0, (i=1, 2.-. 6).research area (i.e. the goaf being simulated) has to bedissected to a certain amount of triangular prisms, asThe tri-prism element (e) is divided into sixshown in Fig. 1. Every triangular prism has six verti-tri-prism sub-elements (C1 ) to (e6*) by its four planesces, which are PI, P2, P3, P4, Ps and Po, as shown inof PIP4P12P10, P2PsP12P1o, PsP&P12P1o and PrPgPg.Fig. 2. Secondly, another six virtual vertices are set,We assume that the unknown concentration functionwhich are P7, P&, P9, P1o, Pu and P12, as shown in Fig.of each sub-element can approximately be expressed2. Now, suppose the concentration of the six verticesas a linear combination of its six vertices. For exam-are C1, C2, C3, C4, Cs and C6 separately, then theple, the sub-element (e1 ), whose six vertices are P1,concentration of the six virtual vertices can be repre-P2, P1o, Pr, Pg and Pu, then its concentration functionsented the linear combinations of the CI to C6, ascan be represented as follows:shown in Eq. (5). .C(a,y, z, ()=q1 C1+o2 C2+93' Cio+o4° C+φs Cg+P6 C1.(6)where φi° to φ° are the base functions of thesub-element (e1 ). They area*(x,y,z)=(z~ -z* )(x,y),B"(x,y,z)=(z~ -z )q%(x,y), .'°(x,y,z)=(z~ -z )Pw(x,y),Fig. 1 Mesh partitin drawingp"(x,y,z)=2z* Pu(x,y), .q',"(x,y,z)=2z 9y(x,y),o%'(x,y,z)=2z Pom(x,y),P.f_乙一 zPsOz=石-石where Pki, Paj and Pkm are 2-D linear base functions ofthe bottom triangle ijim of the sub-element. Using Eq.TePg(5) to eliminate the virtual nodes of C7, C&, C10 andC1 in Eq. (6), the concentration functions can be ex-pressed as:2PobE-dioiP:0)C(x,y,z,t)=,之cC,()p,"(x,y,z)(7)Fig. 2 Virtual nodes of an elementOther sub-element expressions, (e2 ) to (e6 ) cansimilarly be deduced. In Eq. (7), n" is the linear[C, = AC1+ AC,expression of Pki, Pxj and Pom.Cg=hC2+ sCs,3.2 Solution of the mass conservation equationC, =;C3+ nC%,aC dC ac(5)For the above preparation,Ot’ dx’ dyandCo=aC;+@2C2+@2C,C = σC1o +σ2C12,acof every sub-element can directly be solved.)z(C2=@.Cs+a;Cs+o%CAccording. to Green's formula and the continuityequat中国煤化工on equation of anywhere h, w and σ are all weight factors. They agreeverte_region can be ob-with the following conditions:tainedMHCNM H G'WU Zheng-yan et alStudy of 3-D Numerical Simuation for Gas Transfer in the Goaf of the Coal Mining155The model is meshed to a finite element grid byM)dR=automatic mesh generation software, shown in Fig. 4.There are three levels of fitting vertices, 60 vertidesaCand 48 elements. Vertex 1, 4, 7, 10 and 13 are as-f{DnC n。ds- fdR(8)n dx,'sumed to be essential boundary conditions of the headand other vertices are correspond tonaturalwhere Rp is the proprietary sub-region of vertex P andboundary conditions, i.e., dH= 0. It is furthermore(Sp) the interface of Rp.)nEq. (8) is solved based on the element with vertexassumed that vertex 1, 7, 10, 13, 16, 19, 22, 25, 28,P as the vertex. For example, supposing the vertex P31, 34, 37, 40, 43, 46, 49, 52, 55 and 58 are essentialis the vertex P shown in Fig. 2, the part of the pro-boundary conditions of gas concentration and othersprietary subregion in the element (e) is the polyhe-are correspond to natural boundary conditions, i.e,dron with Pr, Q1, P1o, Q4, P7, Qz, P11 and Qs as verti-aC_(; 0. We finally assumed that the permeabilityces. Given a space limitation, the integral derivationof each part is omitted. By adding the discrete equa-tion of the spatial parameter from the integral resultscoefficients in three directions of x，y， Z areof each part, we obtain the ordinary differential equa-K,=K,=K=5 m/d, the porosity is 0.2, the longitudinaland latitudinal dispersions are respectively 50 m'/dtion related with C and ac. The following system ofand 15 m'/d.equations can be written:_57|[A]C+[B]+F=0Level #3133dt{C= θC+u +(1-0)C,(9)32dC_ Citu-CLevel#2 14必15名药大。where [A] and [B] involve weight factors. AccordingLevel#I洲一科二3之时之教为to the equation system (9), the following equation can71610be obtained:[T]C+o=RFig. 4 Finite element mesh drawingwhere [T]= e[A]+B), R= 8(B1-(1- A)[A]}C,-F .The solution of the mathematical model of gasStSupposing that the content distribution C in time t istransfer is programmed using a weighted upstreamgiven, then the concentration ditribution C+sr in themethod. Before calculations are carried out, thenext time step (+Ot) is obtained according to Eq. (9).model needs the establishment of raw data files, in-cluding the coordinate values of x, y and z, longitudi-For a loop, all the solutions is arrived step in step.nal and latitudinal dispersion, porosity, location of the4 Example of a Numerical Simulationsource or sink of .the head and concentration, locationof essential boundary conditions of the head andTo simulate the gas transfer in the goaf of a coalconcentration, etc. The program flow chart is shownmining, we consider a calculation model as shown inin Fig. 5 where CO0.DAT and H0O.DAT are files ofFig. 3, where the goaf is simulated to a cube withdimensions of 60 m x 60 m x 40 m. A longwallRaw data inputingworking face with a length of 60 m and U type venti-lation mode. The goaf is 60 m long and 40 m high.leset?Read H0.DAT+NWrite CO.DAT↓I Read H00.DAT| Read C00.DATGoafSolve the cquation (1)] Solve the equation (2)40中国煤化工:t: into CO.DAT]-6YHCN MHG-Wokingface o r( EndFig. 3 Calculation modelFig. 5 Program frame156Journal of China University of Mining & TechnologyVol.17 No.2initial values for concentration and pressure andof gas transfer in time intervals of 100 and 200 stepsC0.DAT and H0.DAT are files for correspondingin the region of level #2 and level #3. In Fig. 1, thecalculation values.coordinates x, y and Z respectively mean the strike ofThe curved surface graphs of gas pressure and con-the goaf, the length of the face and values of pressurecentration distribution are drawn by Graphtool soft-or concentration, where the units respectively areware. Figs. 6 and 7 show respectively the distributionx9.8 Pa and %. .5504030>30404030> 30,302010 3000 5050. 504030<403022010 1020 Lovele2z2010n;10Level料3(b) Time interval of 200 stepsFig. 7 Curved surface graphs for concentration distibution5 Conclusionsgas transfer leads to the following conclusions:1) The distribution of the air pressure head is veryThe references [7, 11, 13, 14] reported that the gassteady and don't changes as time goes on. And thereconcentration of the goaf along the air intake side wasare a similar shape for each level, as shown in Fig.6.little and higher at the return airway side. More farCertainly, if the boundary conditions change, thefrom the working face, the gas concentration in theshape「 中国煤化工from the wokiggoaf was more higher. At the vertical direction in theface, trence is. these so-goaf, the higher the height, the higher the gas concen-lutions{YHCNMH (jions of the modeltration was.test in rau dclciling w ICICICILC l] ).From our study by simulation in 3-D of gas pres-2) Fig 7 shows the gas concentration distribution,sure and concentration in the goaf, our results aboutfrom which we can obtained the same conclusions asWU Zheng-yan etalStudy of 3-D Numerical Simuation for Gas Transfer in the Goaf of the Coal Mining157the references [7, 11, 13, 14], that is the gas concen-creases as the time goes on. .tration of the goaf along the air intake side was lttleIt is clear that the results of this numerical simula-and higher at the return airway side. And more fartion exactly match the observed results. The samefrom the working face, the gas concentration in theconclusions are drawn. Therefore, numerical imple-goaf was more higher. At the vertical direction in thementation is feasible to solve the dynamic dispersiongoaf, the higher the height, the higher the gas concen-of gas at a goaf.tration was. Another, the gas concentration will in-References1]Zhou s N, Lin B Q. Gas Existing and Flowing Theory in Coal Seams. Bejing: China Coal Industry Publishing House, 1999.(In Chinese)2] Nie B s, He X Q, Wang E Y. Diffusion mode of methane gas in coal pores. Mining Safety & Environmental Protection, 2000,10, 27(5): 14-16. (In Chinese)[3] Wang J S, Liang B, Fan H B. The analytic solution of gas flow differential equation employing dissociation variable method.Journal of Heilongjiang Institute of Science, 2001, 6, 11(2): 39- 41. (In Chinese)4] Wang J R, Liang D, Du W H, et al. Gas migration law in tunnel airflow. Journal of Fuxin Mining Institute, 1991, 7, 10(3):24 -28. (In Chinese)5] Zhang D M, Liu J Z. Regularity on the distribution of gas moving in gob area of coal mine. The Chinese Journal of Geologi-cal Hazard and Control, 2003, 3, 14(1): 82- -84. (In Chinese)[6] Ding GX, Di Z Q, Ma W X. Flow function equation of 2D nonlinear vadose in mined out area and its finite element methodsolution. Journal of China Coal Sociery, 1993, 4, 18(2): 19 -25. (In Chinese)7] Qi Q J, Feng S H, Bai F C, et al. Mathematical models and finite element solution of stope three dimension steady seepageand gas distribution. Journal of Fuxin Mining Institute, 1990, 8, 9(3): 38- 46. (In Chinese)8] Liu J. Gas density distribution of work field determined by using three dimensional finite element technique. Joumal ofFuxin Mining Instiute, 1990, 8, 9(3): 47. -52. (In Chinese)9] Jiang S G Study on 3-D Modeling and Numerical Simulation of Methane Transfer in Top-Coal Drawing Field [Ph.D. Disser-tation]. Xuzhou: China University of Mining & Technology, 1994. (In Chinese)10] Jiang S G, Zhang R W. Mathematical model and numerical calculation of air flow field in sub-level caving workings. Journalof China Coal Society, 1998, 23(3): 258 -261. (In Chinese)[11] JiangS G Wang X s. Study on 3-D model of methane transfer in top coal drawing field. Journal of China University ofMining & Technology, 1995, 24(4): 85- 91. (In Chinese)[12] Qin Y P, ZhuJ F, Chen Y Q, et al. Simulation test and research on gas migration law in goaf of fully mechanized coal cavingmining face. Coal Science and Technology, 2003, 11, 31(11); 13-16. (In Chinese)[13] Li z Q, Zhou P s. Gas emission preventing in the goaf of fully mechanized face with top coal drawing. Safery in Coal Mines,1993, (1): 8 -14, 28. (In Chinese)[14] Fuxin Mining Institute, Fushun Branch of China Coal Research Institute, Songzao Bureau of Mines. Study on gas contentdistribution law of a goaf. Coal Science and Technology, 1992, (3): 15- 17. (In Chinese)中国煤化工MYHCNMHG