水利论文翻译原文及译文课堂基础.doc

Study of flow and transport in fracture network using percolation theoryHaihong Mo a,Mao Bai b,*,Dezhang Lin b,Jean-Claude Roegiers b Department of Civil Engineering, South China University of Technology, Guangzhou, China b Rock Mechanics Institute, The University of Oklahoma, Norman, OK 73019,USAReceived 27 August 1996; received in revised form 3 February 1998; accepted 3 February 1998AbstractA fracture network model based on algebraic topology theory has been developed to study fluid flow and solute transport in fracture dominant media. The discrete fracture network is generated stochastically while the flow and transport in each individual fracture are solved using a continuum approach. Alternative to the traditional formula-tions, the fracture geometries, connectivities and flow characteristics follow a topological structure appropriate to the class of fracture networks or systems. The innovative natures also include the development of a transport mecha-nism in which dispersive-convective solute migration through discrete fracture network can be evaluated. 1998 Elsevier Science Inc. All rights reserved.Keywords: Percolation; Fractures; Topology; Flow; Transport4.2. IntroductionIn the fractured dominant media,fluid flows mainly through fractures. Since the hydraulic be-havior of a rock mass is largely determined by the geometry of the fracture system,the field in-vestigations of flow in fractured rock have clearly demonstrated that modeling the rock as a homogeneous continuum can create serious simplifications. Goodness-of-fit statistics show that the semi-analytic models do not calibrate as well to measured heads as do the continuum models. A porous media modeling approach is restricted because of the heterogeneities,such as fractures, which cannot be explicitly included in the modeling. Thus,in order to analyze the hydraulic be-havior of massive rocks, the necessary step includes a detailed description of the characteristics and geometries of individual discontinuities and of the discontinuity system.The discrete approach involves modeling fluid flow and mass transport in a network intended to quantify hydrogeological factors such as storativity and conductivity in complex domains. Due to the importance and practicality of the subject,numerous discrete fracture network models have been developed and tested 1-7. Using topology as a framework for the fracture system refers that the topological information is explicitly available and that it serves as the organizing factor in the data structure used in the presentation of flow and transport systems. The advantage ofusing this approach is to provide a unified total structure where all topological information is as-sociated in a concise form. In view of topological approaches,Barker 8 provided a proof that the reciprocity principle applies to Darcian flow in a network of heterogeneous branches. Lin and Fairhurst 9 are among the first to apply algebraic topological theory to describe the flow in the network.Dershowitz and Einstein 7 provided a comprehensive discussion to characterize the natural fractures by their shapes,roughness,sizes,locations,spacings,density and orientations. Among these factors, the fracture lengths and orientations may be more critical in the determination of the induced parametric anisotropy. Although the fracture system plays a dominant role in the flow and transport process 10,significant testing efforts still focus on the characteristics of a sin-gle fracture due to experimental constraints 11.Many properties of a macroscopic system are determined by the connectivity of the system el-ements. The special properties of a system,emerging at the onset of macroscopic connectivity within it,are known as percolation phenomena 12. It is a misconception that all fractures in a system contribute to the total flow. Only those fractures,which are connected from initial to ter-minal points of interest,are contributing to the total flow (i.e.,percolating fractures). More rep-resentative work in applying the percolation theory to hydrology can be referred to Jerauld et al. 13,Charlaix et al. 14,Sahimi 15,Silliman and Wright 16,Sahimi and Imdakm 17,Renault 18,Balberg et al. 19,and Kueper and McWhorter 20.Topological concepts may appear rather abstract. However,the steps to accomplish the topo-logical approaches are in fact quite simple. When the network elements may be characterized by linear equations (algebraic or differential) with constant coefficients,a network formulation can be conveniently derived essentially via purely algebraic steps. In this paper,topological aspects of networks composed of fracture branches and sources or sinks are introduced, followed by the algebraic definitions of differential operators involved in the characterization of flow and transport systems. Rules of flow and transport in the fracture network are established to depict critical flow paths while ignoring the nonpercolating dead-end fractures. The tasks are minimized through producing an efficient computer program. Examples are given to illustrate the changing pressure and concentration profiles due to fluid flow and solute transport from a stochastically generated fracture network.4.3. Topological presentation of fracture networkA computer-generated sample of randomly aligned line segments is illustrated in Fig. 1. The line segments may be considered as fractures embedded in a rock with very low matrix permeabil-ity (e.g.,granite). The terms network and one-dimensional complex are synonymous. They both refer to a mathematical structure that consists of the various interconnections of the branches from a fracture network. This structure is constructed in such a form that allows the study of fluid flow in somewhat abstract manner. The topology concept of a network is important for the sys-tematic formulation. One of the essential topological properties of a network is to use simple line segments to obtain the graph of the network.4.3. Oriented graph for fracture flow networkAn oriented graph G(N, E) consists of a set,N (elements of nodes),and a set E of ordered pairs of the forms (n“ n7),while n“ nG N (oriented edge of G(N, E) or branches). nt and n) are the ini-tial and terminal nodes (end points),Fig. 1. Stochastically generated fracture network.The branch is the connection in an orientated path. The simple rule is that when the branch is connected from its initial point to its final point, the mark +1 is given. The connection with the opposite direction is marked -1. Thus each path determines a vector e = (e1, e2,., ek,. .)T with integer coefficients, where the coordinates of the vector e are labeled by the branches of the graph; with ek = +1, for example, the branch is indicated in the positive direction. Similarly, ek = -1 marks a branch with a negative direction.A flow rate distribution of the network can be viewed as a vector q = (qe1, e2,., qe. .)T whose coordinates are labeled by the branches, where q, for example, is the real number related to the flow rate through the branch et while assuming that the one-dimensional fracture network is the complex of branches of a two-dimensional fracture system. In other words, the fractures are located in a 2-D domain where the flow in the individual fracture is in 1-D setting.4.3. Basic topological conceptsA one-dimensional network is a collection consisting of two sets: a set of nodes (zero-dimen-sion) and a set of branches (one-dimension). The chain,denotes a vector space. The vector space, C0, whose components are indexed by the nodes, are called zero-chains. Similarly, C1, with components being indexed as branches, denotes one-chains. A branch is an element of the space C1. Each node i is identified with a vector that has one in the ith position and zero elsewhere. Thus, e1 = (1,0,0,. t)T, e2 = (0,1,0,. .)T, etc. More specifically, a branch vector space S = (51,52,Si,.) can be expressed with components indexed by the branches. It is defined that dim C1 is the number of branches and dim C0 the number of nodes.The boundary map,9,is defined as a linear map from C1 to C0. To define the map d it is suf-ficient to prescribe its values on each of the branches,since the branches form a basis for C1. Each branch has an initial point and a terminal point and a branch is equal to the subtraction between the terminal and initial nodes. Thus,for example,if a branch goes from nt to n,ek = n) n.)1 2U(The above concept can be illustrated by Fig. 2. Considering an oriented graph G(N, E) of frac-ture network,one has(3) (5)(6)Further rules are set as: the elements of boundary operator 9 are (a) set to 0 if the branch j is not associated with node k; (b) set to +1 if fluid in branch j flows toward associated node k; and (c) set to -1 if fluid in branch j flows away from associated node k. Since each branch leaves a single node,each column of matrix 9 contains a single +1 and a single -1,with all other elements being equal to zero. For a connection with four nodes and four branches,this relationship can beexpressed in a matrixformas10+11 99=+11000+110000+1Fig. 2. A network mesh and oriented graph.4.3. Connectivity of a networkRecall that a path joins the node na to the node 叫 if na is the initial point of the first segment of the path and nb is the final point of the last segment of the path. If w is the one-chain with integer coefficients corresponding to this path, then it is clear that dw = nb na. Since we can compute dq by adding the boundaries of all the individual branches, all the intermediate nodes vanish. This is the essence of the topological concept.4.3. Meshes for the networkIf the subspace of C1 consists of one-chains satisfying dZ 1 = 0, Z1 is called the one-dimensional cycle. A closed path (0, n1,., nr) is the one with properties that 0 = nr and the branch (0, n1),. (nr1, nr) are pairwise distinct. A simple closed path is called a mesh M. Clearly, each of the meshes has no boundary, or dM = 0. For example, in the complex shown in Fig. 3, it is possible to define three meshes: M1 = e1 + e4 + e5, M2 = e2 + e3 e5, and M3 = e1 + e2 + e3 + e4. However, only two of these meshes are independent because M3 = M1 + M2. Each mesh deter-mines an element, M, of C1, whose coordinate is +1, -1 or 0, and dM = 0. Every mesh is a cycle, but not every cycle is a mesh.A connected complex containing no meshes is called a tree. In a tree there exists at least one node which is a boundary point of only one branch. Thus, in a tree, the number of nodes is exactly one more than the number of branches: if Nb denotes the number of branches and Np denotes the number of nodes, then in any tree: Np = Nb + 1 (see Fig. 4).4.3. Structure of co-chainsUp to now, only spaces C0 and C1 are introduced. Because the nodal pressure pn may be con-sidered as the quantities of another but related space, thus a pressure vector p = (p1,p2,.), and is assumed to lie within the dual space,of the space C1. This dual space is called space of one-co-chain, and denoted by , which is also the space of linear functions on C1. Similarly, the spaceFig. 3. A closed path.(7)Fig. 4. A family of two meshes.n12of linear functions on C0 can be introduced as the space of zero-cochains. The adjoint boundary operator,which is denoted by d,is a linear map from to ,and is defined as coboundary op-erator. Since e is a branch and 0e = n2 n1,there exists a function,pn,on the nodes,i.e.,a zero- cochain, such that pb = pni pnx. For the network shown in Fig. 2,the coboundary operator d has the following matrix representation:1+10001+10d=+1010100+1Recall Eq. (6),itcanbe seen th4.3. Matrix form of restrict boundary operators and restrict coboundary operator(8)When evaluating fluid flows through a network,a reference nodal pressure must be known. The restricted boundary operator 9* can be obtained by simply deleting components correspond-ing to prescribed nodes from expression as a vector in C0. For the network shown in Fig. 2,one can eliminate one row from the matrix expressed in Eq. (6),since this row is the negative sum ofall the other rowsThus, 9* becomes10+1 19* =+110 00+11 0The mapping d*,induced from a map d,is defined as the restricted coboundary map,whose adjoint is the map 9*. The matrix representing d* is just the transpose of 9*,or,1+1001+1+1011004.4. Modeling flow and transport in the fracture network4.4. Steady state flowUnder creeping-flow conditions the equations governing the flow in the network are linear. As-suming quasi-steady state, the problem for a certain flow situation is reduced to solving a system of linear algebraic equations.The nodes of the network (excluding the boundary ones) are numbered in a convenient but ar-bitrary manner by assigning indices k = 1,2,.,Np,where Np is the number of interior nodes. The branches of the network are given indices,=1,2,Nb,where Nb is the total number of branches. The positive flow direction on each branch can be chosen arbitrarily. The flow rates and pressures for the network are defined as: (a) qk and pbk are the flow rate and pressure across branch k; and (b) qsk and psk are the independent flow rate source and pressure in branch k.If办 is a one-chain describing the flow rate distribution of a network, then9*?b + 9s = 0,(10)where 9* is given in Eq. (8),and 仏 is the node flow-source vector,which has Np elements.Let pn be the vector of node pressure,i.e.,a zero-cochain,the branch pressures pb is designated by the relationshipPb = d*Pn(11)where d* is described by Eq. (9).The sum of pressure drop around each mesh, due to both the mesh flow rates and the sources, is zero as required by Kirchhoffs law. The relationship between branch flow rate 9b,branch pres-sure pb and the branch source pressure 夕s can be represented by9b = G(pb - pj,(12)where G is the branch conductivity matrix. For the chosen oriented graph of the network,G has the dimension Nb x Nb,and is of isotropic nature.Substituting Eqs. (11) and (12) into Eq. (10),yields,(9*Gd*)n = 9*Gs 9s.(13)Each diagonal entry in 9*Gd* is the sum of the conductivities of all the branches connected to a node. The component of (9*Gd*)n represents the net flow rate which would flow out of each node if there were no sources. The components of (9*Gs 9s) are the net flow rates which would flow into each node if all node potentials were equal to zero.The node pressure pn can thus be determined from the expressionA = (9*Gd*)- 1(9*Gs 9s).(14)Boundary conditions at the inlet and outlet of the network can be prescribed. In our calcula- tion,constant pressure boundary conditions are assumed.For convenience and simplicity, adding boundary conditions to Eq. (13), the following formu-lae may be obtained:where g are conductivity of branch ij. With the help of Eq. (15), an efficient computer program can be easily developed.4.4. Unsteady state flowEach branch of the network is bounded by two nodes i and j, and there is a flow along the branch which is characterized by a diffusivity tensor kj. A one-dimensional flow equation is ap-plied to each fracture branch of the network, i.e.,In the fracture network, the drawdown p in each branch obeys the flow equation, with the initial condition Pij(x, 0) = 0 and two boundary conditions described asPij(0, t) = Mt),(19)Pij(Lih t) = Pj(t)-(20)where L” is the length of the fracture being examined.The flux at each node i from branch ij, or Darcy,s law, can be described asq = kij(0)(21)and, similarly,For any node i the flux is described by29 Through substituting Darcys solved asPij (x, s) = Aij exp (ijx) + Bijwherelaw into flow equation, the pressure in Laplace space can be exp (AjX),(24)where Aij and Bij are determined from the above boundary conditions and s is the Laplace trans-form parameter.After the determination of Atand B,the solution isk=i k=ip(s),can be obtained by solving Eq. (28). Then this function can be inverted into real space by means of algorithm described by stehfest 21. The solution for the drawdown at the nodes for the entire network can be completed by summing up all paths from the inlet to the outlet of the net-work.4.4. Solute transportThe developed mechanism can be applied to the analysis of solute transport using conventional dispersion-convection concept. For the radial flow geometry,the following equation can be givenwhere r is the distance from the a source/sink (e.g.,a well),v the average flow velocity,c the mean concentration,and Dr the radial dispersion coefficient.For the fracture network, initial and boundary conditions can be expressed asCi(Xi,0) = 0,(32)Ci(0,i) = Ci(?),(33)i(Li,i) = c(i). (34) The flux at node i for branch ij isq*(i) = SiDr Dfi),(35)where L” is the fracture length,and Siis the fracture cross-section area.Assuming that there is no storage in the nodes, mass balance for any node i is described byEq*(i) = Q*(i).(36)i=Using Laplace transform, the nodal concentration can be expressed asCi(x,s) = Aiexp(ai7x) + Btexp(x),(37)a, p = (0.5Dr)vij 士 (v+ 4Drs) 4.1. Fluid flowThe presented percolation method is evaluated through examining a scenario where fluid flow across a domain exposed in the vertical cross-section, as shown in Fig. 5. Within the rectangular domain, upper and lower boundaries are assumed to be impermeable. The flow direction is pri-marily in the horizontal direction starting from the right-hand side where the constant flux is spec-ified. Fluid exits from the lef