UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Estimation of grout distribution in a fractured rock by numerical modeling Rahmani, Helia 2009

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
24-ubc_2009_spring_rahmani_helia.pdf [ 2.41MB ]
Metadata
JSON: 24-1.0063121.json
JSON-LD: 24-1.0063121-ld.json
RDF/XML (Pretty): 24-1.0063121-rdf.xml
RDF/JSON: 24-1.0063121-rdf.json
Turtle: 24-1.0063121-turtle.txt
N-Triples: 24-1.0063121-rdf-ntriples.txt
Original Record: 24-1.0063121-source.json
Full Text
24-1.0063121-fulltext.txt
Citation
24-1.0063121.ris

Full Text

ESTIMATION OF GROUT DISTRIBUTION IN A FRACTURED ROCK BY NUMERICAL MODELING by HELIA RAHMANI B.Sc., University of Tehran, 2004  A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE  in  THE FACULTY OF GRADUATE STUDIES (Civil Engineering)  THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver)  May 2009  © Helia Rahmani, 2009  ABSTRACT Grouting has been used over the past two centuries to increase the strength, decrease the deformation and reduce the permeability of soils or fractured rocks. Due to its significance in engineering and science predicting grout effectiveness in fractured rocks is of interest. There are different approaches to estimate the effectiveness of grouting, one of which is numerical modeling. Numerical models can simulate distribution of grout inside fractures by which the effectiveness of grout can be estimated. Few numerical studies have been carried out to model grout penetration in fractured rocks. Due to complexities of modeling grout and fracture most of these studies have either used simplifying assumptions or been bound to small sizes of fractures, both resulting in unrealistic simulations.  The current work aims to eliminates some of the simplifying assumptions and develop a model that can improve the reliability of the results. In reality grouts are believed to behave as a Bingham fluid, but many models do not consider a full Bingham fluid flow solution due to its complexity. Also, real fractures have rough surfaces with randomly varying apertures. However, some models consider fractures as planes with two parallel sides and a constant aperture. In the current work the Bingham fluid flow equation are solved numerically over a stochastically varying aperture fracture. To simplify the equations and decrease the computational time the current model substitutes two-dimensional elements by one-dimensional pipes with equivalent properties. The model is capable of simulating the time penetration of grout in a mesh of fracture over a rather long period of time. The results of the model can be used to predict the grout penetration for different conditions of fractures or grout.  ii  TABLE OF CONTENTS ABSTRACT.................................................................................................................................... ii TABLE OF CONTENTS............................................................................................................... iii LIST OF TABLES......................................................................................................................... vi LIST OF FIGURES ...................................................................................................................... vii LIST OF SYMBOLS ...................................................................................................................... x  1  2  INTRODUCTION .................................................................................................................. 1 1.1  Background to Grouting ................................................................................................. 1  1.2  Objectives of the Research.............................................................................................. 2  1.3  Overview of the Research............................................................................................... 5  OVERVEIW OF GROUTING ............................................................................................... 6 2.1  Grouting History ............................................................................................................. 6  2.2  Grouting Techniques....................................................................................................... 9  2.2.1  Grouting in Soil....................................................................................................... 9  2.2.2  Grouting in Rock................................................................................................... 15  2.3  Grout Types Based on Rheological Behavior and Material ......................................... 17  2.3.1 2.4  Bingham Grout Materials ..................................................................................... 19  2.4.1  Fracture Characterization.............................................................................................. 25  2.5  Fracture Aperture .................................................................................................. 28 Fluid Flow in Fractures................................................................................................. 31 iii  2.5.1 2.6  Groutability ........................................................................................................... 31 Methods to Estimate the Grout Distribution in the Ground.......................................... 35  2.6.1  Pre and Post Grouting Tests.................................................................................. 35  2.6.2  Modeling ............................................................................................................... 40  3  THEORETICAL MODEL FOR ONE DIMENSIONAL FRACTURE FLOW................... 48  4  NUMERICAL MODELING OF GROUT PENETRATION IN A FRACTURE ................ 54 4.1  Objectives of the Model................................................................................................ 54  4.2  Simplifying Assumptions in Simulation....................................................................... 54  4.2.1  Grout Property Assumptions................................................................................. 55  4.2.2  Fracture Geometry Assumptions .......................................................................... 55  4.2.3  Grout Flow Assumptions ...................................................................................... 56  4.2.4  Boundary Conditions Assumptions ...................................................................... 56  4.2.5  Mesh Generation Assumptions ............................................................................. 57  4.3  5  Algorithm...................................................................................................................... 58  4.3.1  Random Fracture Generator.................................................................................. 58  4.3.2  Solver .................................................................................................................... 66  VERIFICATION CASES ..................................................................................................... 75 5.1  Overview of the Chapter............................................................................................... 75  5.2  Verification Case I: Fracture Generator........................................................................ 75  5.3  Verification Case II: Continuity Equations in the Solver ............................................. 77  5.4  Verification Case III: Head Drop Equations in the Solver and Network Update......... 81  iv  5.5  6  Verification Case IV: Water Case................................................................................. 84  SAMPLE SIMULATIONS AND RESUTLS....................................................................... 90 6.1  Overview of the Chapter............................................................................................... 90  6.2  Input Properties for Simulations ................................................................................... 90  6.3  Analyzing the Results ................................................................................................... 91  6.4  Sample Simulations ...................................................................................................... 94  6.4.1  Overview of the Simulations ................................................................................ 94  6.4.2  Grout Penetration for Different COV Values ....................................................... 95  6.4.3  Grout Penetration for Different Fracture Apertures............................................ 101  6.4.4  Grout Penetration for Different Yield Stresses................................................... 105  7  SUMMARY AND CONCLUSION ................................................................................... 110  8  REFERENCES ................................................................................................................... 112  APPENDIX A: DERIVATION OF THE BINGHAM FLUID EQUATIONS .......................... 118 APPENDIX B: NEWTON’S METHOD FOR THE SOLUTION OF PIPE NETWORK EQUATIONS.............................................................................................................................. 124 APPENDIX C: RANDOM FRACTURE GENERATOR ........................................................ 129 APPENDIX D: PIPE NETWORK SOLVER CODE ............................................................... 133  v  LIST OF TABLES  Table 2.1: Typical range of properties for Bingham grout (Rombough, 2006)............................ 20 Table 5.1: The input properties of the pipe network for verification case II. ............................... 78 Table 5.2: Code results for verification case II............................................................................. 80 Table 5.3: Hand calculation results for verification case II. ......................................................... 81 Table 5.4: The input properties of the pipe network for verification case III............................... 83 Table 5.5: Code results for verification case III. .......................................................................... 83 Table 5.6: Hand calculation results for verification case III......................................................... 84 Table 5.7: The input properties of the pipe network for verification case IV. ............................. 86 Table 5.8: Code results for verification case IV. .......................................................................... 87 Table 5.9: Hand calculation results for verification case IV. ....................................................... 88 Table 5.10: Transmissivity in different time steps from the code result and hand calculations... 89 Table 6.1: Grout material properties used in the simulations ....................................................... 91 Table 6.2: Summary of the simulations ........................................................................................ 95 Table 6.3: Input properties for random fracture generation.......................................................... 96  vi  LIST OF FIGURES  Figure 2.1: Charles Beringy, about 1810 (Grouting Procedure) with a River Seine bridge, Sevres, France (modified from Kutzner, 1996)................................................................................... 6 Figure 2.2: Jet grouting technique (Henn, 1996). ......................................................................... 10 Figure 2.3: Compaction grouting in soils (modified from Warner, 2004).................................... 11 Figure 2.4: Longitudinal view of compaction grouting used to control settlement in conjunction with excavation of soft ground tunnel (Henn, 1996) ............................................................ 12 Figure 2.5: Groutability based on grout type versus soil particle size (Henn, 1996) ................... 13 Figure 2.6: Rheological behavior of typical grouts (Mongilardi and Tornaghi, 1986). ............... 18 Figure 2.7: Basic rheological concepts (modified after Lombardi, 1985).................................... 18 Figure 2.8: Fracture properties determined by fracture void geometry (from Olsson and Barton 2001). .................................................................................................................................... 27 Figure 2.9: Definition of aperture (Hakami and Larsson, 1996) .................................................. 28 Figure 2.10: Illustration of mechanical aperture (modified from Evans et al, 1992) ................... 29 Figure 2.11: Active forces on grout in a slot ................................................................................ 33 Figure 2.12: Pressures for each run of the five pressure water test (after Houlsby, 1990)........... 38 Figure 2.13: Fracture network model by DFN method (Shuttle and Glynn, 2003)...................... 44 Figure 3.1: Schematic display of the inflow and outflow in a node ............................................. 49 Figure 3.2: Forces in a flowing fluid inside a channel (modified from Chhabra and Richardson, 1999). .................................................................................................................................... 51 Figure 4.1: Schematic display of finger model............................................................................. 57 Figure 4.2: General form of normal distribution (Galambos and Kotz, 1978). ............................ 59  vii  Figure 4.3: Sample input for the fracture generator...................................................................... 61 Figure 4.4: Sample output of the 2D elements of the fracture generator...................................... 63 Figure 4.5: Replacing 2D elements with equivalent 1D pipe network. ........................................ 64 Figure 4.6: Averaging properties of adjacent elements. ............................................................... 65 Figure 4.7: 1D equivalent elements .............................................................................................. 65 Figure 4.8: The pipe network in the first time step....................................................................... 67 Figure 4.9: Network update in each time step. ............................................................................. 68 Figure 4.10: Pressure drop in the pipe network. ........................................................................... 69 Figure 4.11: A general form of pipe network. .............................................................................. 71 Figure 5.1: Sample input for the random generator...................................................................... 76 Figure 5.2: Distribution of the aperture values generated for Verification Case I. ...................... 77 Figure 5.3: Pipe network for verification case II with four similar pipes..................................... 78 Figure 5.4: The input file containing grout and grouting properties for verification case II........ 78 Figure 5.5: Pipe network for verification case III......................................................................... 82 Figure 5.6: The input file containing grout and grouting properties for verification case III....... 82 Figure 5.7: Pipe network for verification case IV with four similar pipes. .................................. 85 Figure 5.8: The input file containing grout and grouting properties for verification case IV. ..... 86 Figure 6.1: Cylindrical penetration of grout. ................................................................................ 93 Figure 6.2: Generated aperture distribution (a), (c), (e) and grout penetration contours for 2 minute (120s) intervals (b), (d) and (e) for different Coefficients of Variation. .................. 97 Figure 6.3: Changes of injection pressure with grout penetration for different Coefficients of Variation. .............................................................................................................................. 99 Figure 6.4: Changes of injection pressure over time for different Coefficients of Variation..... 100  viii  Figure 6.5: Penetration length of grout for different Coefficients of Variation.......................... 101 Figure 6.6: Generated aperture distribution (a), (c), (e) and grout penetration contours for 2 minute (120s) intervals (b), (d) and (e) for different fracture apertures. ............................ 102 Figure 6.7: Changes of injection pressure with grout penetration for different apertures.......... 103 Figure 6.8: Changes of injection pressure over time for different fracture apertures................. 104 Figure 6.9: Penetration length of grout for different fracture apertures. .................................... 105 Figure 6.10: Generated aperture distribution (a) and grout penetration contours for 2 minute (120s) intervals for different yield stresses (b) and (c). ...................................................... 106 Figure 6.11: Changes of injection pressure with grout penetration for different yield stresses (cohesions). ......................................................................................................................... 107 Figure 6.12: Changes of injection pressure over time for different yield stresses (cohesions). . 108 Figure 6.13: Penetration length of grout for different yield stresses (cohesions)....................... 109  ix  LIST OF SYMBOLS  e  Aperture  b  Aperture (considered as half of e in the calculations)  c  Cohesion/Yield Stress for the Bingham fluid  ρ  Density  μ  Dynamic viscosity  μB  Dynamic viscosity of Bingham fluid  Q  Flow rate  q  Flow rate per width  g  Gravitational acceleration  Pg  Grout pressure  L  Penetration length of grout  P  Pressure  Pw  Water pressure  τ yz  Shear stress acting on the yz plane in the z direction  T  Transmissivity  w  Width  τ 0B  Yield Stress/Cohesion for the Bingham fluid  x  1 INTRODUCTION 1.1 Background to Grouting Grouting has been used for more than two centuries to improve the engineering properties of soils and rocks. Grouting has three main applications:  Increasing strength: Strength improvement of the soil or rock is one of the most popular applications of grouting. A typical example for this purpose is grout injection under structure foundations to mitigate bearing failure.  Decreasing deformation: Grouting is often performed to reduce soil and rock deformation, hence minimizing damage to adjacent structures. Grouting immediately before tunnel boring in the vicinity of other tunnels or structures is a good example of this application.  Reducing permeability: The other important use of grouting in both soils and rocks, and likely the most widespread use of grouting, is to reduce the permeability of the formation.  Depending on the purpose of grouting and also conditions of soil or rock, the properties of grout such as viscosity, or setting time might be varied. These properties can be adjusted by changing the amount and/or type of the ingredients in grout mixture.  1  In a highly fractured and conductive rock undesirable fluid flow (typically water for geotechnical applications) is a common engineering problem. Typical examples of these unwanted flows include water flow under dams and leakage from oil wells. One method of reducing these flows is to seal the fractures with grout.  How grout distributes itself inside individual rock fractures during injection is critical to grouting: the effectiveness of any grouting project is related to whether the grout forms a continuous sealing layer in the rock. Methods to estimate the grout effectiveness in a medium is non-trivial, and there are a range of approaches used for this purpose, including pre/post grouting tests and modeling.  One of the most meaningful parameters in measuring the effectiveness of grouting is the penetration length of grout during the grouting process. However, the per/post grouting tests are not capable of accurate estimations of this key parameter. Numerical simulations can predict the grout penetration in different directions under different conditions of grouting. Depending on the input parameters of the model, the properties of grout or the media in which grout is injected can be set to required values. Due to their capabilities, numerical models are the best choice in simulating the grout distribution in fractured rocks.  1.2 Objectives of the Research The primary aim of this study is to develop a code in order to determine the time dependent penetration distribution of a Bingham grout injected into a variable aperture 2  fracture. The ability to more accurately predict the grout spread in a fracture network is of great importance: it not only helps the grouting engineer to choose appropriate grout properties but also determines the effectiveness of grouting. The grout flow studies, including numerical and experimental studies, involve gaining realistic understandings of properties of the fluid, geometrical conditions of the fracture and the behavior of the flow inside the fracture. In this study the complexities of these components are combined to achieve a more realistic modeling of the grout flow.  In practice grout penetration is mostly estimated by experimental methods such as pre/post grouting tests. Also, a few numerical models have been developed to study the penetration of grout inside fractured rocks. However, rheological properties of grout and geometric characteristic of fractures, the two factors that significantly derive the time variation and distribution of grout flow inside the fracture, provide a challenge for experimental setups and simplifications required for empirical estimates. Most of the existing numerical methods to our knowledge are not capable of simultaneous consideration of both rheological properties of grout and geometric features of fracture. The few available numerical studies that take both complex fluid properties and fracture stochastic characteristics into account have been restricted to small fracture sizes and short grouting periods due to computational resources limitations. The aim of the current study is to develop a numerical model which realistically represents the rheological properties of cementitious grout, and incorporates stochastic aperture distribution of a fracture.  The model should also be able to simulate a larger number of elements  3  compared to previous studies over a long enough period of time for grout to reach a large area of fracture. To numerically study the evolution of grout, a square mesh of one-dimensional channels is considered representing the fracture surface. Grout is injected from a certain point in the centre of the mesh. The choice of one-dimensional channels was adopted as it saves a considerable amount of computational time, making the model capable of simulating grouting in larger fractures over longer periods of time. The equations of grout flow networks are solved as an initial value problem knowing the injection flow rate at the borehole.  The specific objectives of this study to consider the Bingham behavior of grout and complex geometrical properties of fracture simultaneously require some additional treatments. The behavior of the grouts especially the cementitious grout as a Bingham fluid is not as simple as a Newtonian fluid like water. Therefore, the equations of flow were derived to consider the complex behavior of grout. The variation of aperture over fracture surface is considered to be stochastic using spatially varying field, produced by a random generator.  As grouting continues, grout flow extends to more elements in the network. Solution of the equations gives the penetration length for each element at any time. Having the penetration distribution and pressure of grout at each time step in the fracture, the grout effectiveness can be predicted. Grout spread under different grouting conditions can be studied by changing the input properties of the model. Therefore, the model can be used  4  to find the suitable properties of grout for defined penetration length in any specific fracture condition.  1.3 Overview of the Research The present work includes 6 chapters, organized as follows. Modeling of grout flow in fractured rock requires an insight into three main areas: grout rheology, fracture geometrical properties and the grout flow in the fracture. Chapter 2 contains a review of relevant literature.  After a general background on grouting, an introduction to the  properties of grout and fractures (essential to model the grout fracture penetration) is presented. Prior to development of numerical models, all grouting techniques were developed based on the knowledge gained from experimental and practical studies. Chapter 2 also gives a review of these studies including the grouting methods, materials and controlling experiments. The details on the theory and equations used to model the grout flow is given in Chapter 3. This chapter also presents the derivation of flow equations for a Bingham fluid. Chapter 4 focuses on the numerical modeling of the grouting in a single fracture. This chapter explains the methodology and algorithm adopted to model the grout penetration in a fracture. Verification cases to check the correctness of the numerical results are presented in 5 .  Chapter 6 presents some sample simulations and compares the results of the code.  Finally, the conclusions from the current work are presented in Chapter 7.  5  2 OVERVEIW OF GROUTING 2.1 Grouting History Injecting a cementitious material in the soil for the purpose of improvement is believed to have first been used by the Frenchman Charles Berigny in 1802. Beringy used slurry to fill up the caves in the foundation of a sluice to increase the foundation strength and retrieve its bearing capacity. The settlement damaged foundation which was built on a layer of alluvial soils needed to be sealed to stabilize it. The first known grouting work sketch is from this job, and is shown in Figure 2.1 (Kutzner, 1996).  Figure 2.1: Charles Beringy, about 1810 (Grouting Procedure) with a River Seine bridge, Sevres, France (modified from Kutzner, 1996).  6  Between 1802 and 1809 Berigny used grouting in construction of the harbor basin at Dieppe to decrease the permeability and reduce the inflow. This was the first time that grouting was used for construction purposes and not repair (Verfel, 1989).  Grouting applications became much more popular after the development of hydraulic binders and the invention of Portland cement in 1821 (Kutzner, 1996). Before the invention of the Portland cement Pozzolana cement was mainly used, especially for works with water contact. Pozzolana cement had a much lower strength than Portland, which limited its application (Hall, 1976). Since the development of Portland cement the injection of cementitious grout has been widely used in the construction of many structures. In 1837, Rayal first used grouting to repair a masonry structure. In the years 1856 to 1858 Kinipple experimented with the idea of strengthening masonry by filling the cracks with grout (Verfel, 1989). The first recorded use of grout in underground construction can be traced back to 1864, when Barlow filled the voids left by the tail of a tunnel shield with grout (Henn, 1996).  The beginning of the rock grouting in United States dates back to 1893 when the rock joints at New Croton Dam were systematically sealed by grout injection. In this project there was a design for the arrangement and depth of the boreholes across the entire foundation. In this project grouting was used to seal fractures and decrease the uplift pressure under the dam (Verfel, 1989; Henn, 1996; Littlejohn, 2003).  7  Further developments of Grouting techniques were achieved between years 1900 and 1920 by the improvement of grouting equipment. The grouting techniques introduced in Section 2.2  are basically defined by controlling the injection pressure, and grout  properties like the fineness or viscosity of the grout and filtration of excess water from the grout.  After 1920 a series of inventions had a considerable effect on the development of grouting. One of these significant forward steps was the invention of chemical grout materials by Joosten in 1926. Joosten injected highly concentrated sodium silicate and calcium chloride into gravel and coarse sands to form an impermeable sandstone. After World War II Chemical grouts with very low viscosity were developed to stabilize the soil under the damaged structures. Chemical grouts allow for stabilization of soils down to fine sand with a small content of silt. Another significant advance was made after the Second World War when the wooden pumps where replaced by modern high pressure hydraulic pumps. These new pumps not only enabled the use of higher pressures for grouting, but also allowed the grout pressure and discharge to be controlled much more precisely (Kutzner, 1996).  Much of the grouting development in design, equipment and material are paralleled by advances in the design and building of dam foundations. One example of these advances is the development of arch dams. Since this kind of dams requires stronger foundation and abutments, more effective grouting over a larger area is required. In fact to have a stronger foundation grout should cover larger depths and penetrates more cracks. This  8  requires stronger apparatus with larger grouting pressures and also more penetrable grouting material (for more details on this subject see Jansen, 1988). Currently grouting is used for different purposes by employing different practical techniques.  2.2 Grouting Techniques This section provides a brief introduction of some widely used grouting methods. Selection of the grouting technique depends on the properties of the soil or rock. In soils the preference for the grouting method is defined by factors like soil gradation, ground water level, depth of the structure below the surface, and surface or subsurface access for the grouting equipment. In rocks the same criteria, except the soil gradation, are valid and also the geology of the rock might determine which approach to use.  2.2.1 Grouting in Soil This section describes different types of grouting in soils based on their methodology. The methodologies are varied by changing the grout properties, injection pressure, or the method of applying grout to the soil.  Jet Grouting: Jet grouting, also called replacement grouting, is a relatively new technique building upon an idea proposed in Japan in 1965 (Xanthakos et al., 1994)). It is generally accepted that after early 1980’s this method became economic and practical. Jet grouting  9  developed in response to the need to treat soils ranging from gravels to clays to random fills in areas where environmental controls restricted the available soil improvement options. Right now it is arguably the fastest growing method of ground improvement.  In this method the grout is ‘jetted’ into the ground under high pressure as illustrated in Figure 2.2. The jet simultaneously excavates and mixes the soil with the grout and the nozzle injects high pressure slurry into the soil. To improve the cutting action of the jet compressed air is often added. The whole procedure results in a column of soil mixed with grout. This grout-soil column is designed to have a lower permeability and higher strength than the virgin ground (Henn, 1996).  Figure 2.2: Jet grouting technique (Henn, 1996).  Jet grouting has been used widely in Canada. As discussed by Byrne et al. (1988) this method of grouting has been employed widely, including numerous excavation works in the Montreal region, and in creating a cut off wall under the John Hart Dam, B.C..  10  Compaction Grouting: Compaction Grouting originated in California in the early 1950’s (Gallavresi, 1992). Brown and Warner presented the first report on the actual mechanism of compaction grouting using the data provided from excavations of previously grouted soils (Warner, 2004).  In this method a very stiff, low mobility soil cement mix is injected slowly into the soil at a high pressure level (up to 3.5 MPa). Grout remains in a homogeneous, expanding mass which compacts adjacent soil (Figure 2.3).  Figure 2.3: Compaction grouting in soils (modified from Warner, 2004)  Compaction grouting is performed at discrete locations in soft, loose or disturbed soils, as shown in Figure 2.4. Because of the low slump of the grout it does not enter the pores of the soil, so it forms a small coherent bulb near the injection point (Xanthakos, 1994).  11  Figure 2.4: Longitudinal view of compaction grouting used to control settlement in conjunction with excavation of soft ground tunnel (Henn, 1996)  This technique of grouting is also used to stabilize the soil under residences and light commercial buildings. Since compaction grouting can result in the lifting of the ground surface if not monitored very carefully, this grouting technique might not be suitable for structures that can tolerate only the smallest differential movements (Gallavresi, 1992; Shrof and Shah, 1993).  Permeation Grouting: Permeation grouting is the process of injecting flowable grout into the soil to strengthen and/or reduce the permeability of the soil by filling the existing voids and fissures. This method of grouting is performed at a low injection pressure not to cause major changes in the soil structure (Bruce, 2006).  12  Permeation grouting is particularly sensitive to the groutability of the soil. Groutability, which is explained in detail in Section 2.5.1, is the ability of the rock or soil to take the grout. This property is based on the grout type and the particle size of the soil as approximated in Figure 2.5 (Henn, 1996).  In permeation grouting the aim is to fill the existing gaps with the minimum destruction of the soil structure. If the grout is not suitable for the soil and high pressure is used hydrofracturing might occur, leading to uncontrolled grout movement and the possibility of structure damage. Generally the grading of the soil defines the groutability of the material.  Figure 2.5: Groutability based on grout type versus soil particle size (Henn, 1996)  Permeation grouting technique is regarded as the oldest method of grouting (Warner, 2004). In fact the permeation grouting was invented in the early 19th century when the initial idea of injecting a cementitious slurry into the soil was exploited for the first time.  13  Permeation grouting has a very low effect on the adjacent structures and the surrounding area, hence it has an extremely wide range of application. Some of the categories where this technique is useful are: different soil and rockv type, temporary or permanent works, remedial works, and also for the purpose of strengthening or water proofing (Xanthakos et al., 1994).  Hydro fracture Grouting: As stated in the previous section, the primary aim of permeation grouting is to fill the existing fissures and pores with grout. In this process if the injection pressure increases the existing stresses in the surrounding area, some new fractures would be created and the grout would enter the new fractures. This phenomenon of forming fractures by the pressure of a fluid is known as hydro-fracture. If grouting is deliberately performed in very high pressures (up to 4MPa (Gallavresi, 1992)) to form fractures the grouting method would be called hydro-fracture grouting.  Depending on the properties of the soil the pressure limit to form the fractures would vary. Generally hydrofracture would be initiated when the pressure applied by the grout exceeds the tensile strength of the soil. Since 1957 researchers have tried to develop simple models which can analyze the initiation of hydrofracture. According to Wong and Farmer (1971) among the earliest models developed included those by Hubbert and Willis (1957), Morgenstern and Vaughan (1963), Fairhurst (1964) and Cambefort (1964).  14  Control of the direction and penetration distance of the hydrofracture is difficult. Thus in this method the grout may penetrate far from the injection point and as a result limiting the potential danger to adjacent structures can be very difficult.  2.2.2 Grouting in Rock In this section methods of grouting in rocks are represented. These methods are listed by the purpose of the grouting by which grouting in rocks can be divided into two methods: consolidation grouting and curtain grouting.  Consolidation Grouting: Consolidation grouting is the process of injecting a cement based grout into open joints, separated bedding planes, faulted zones and cavities in the rock mass for the purpose of improving the mechanical properties of the rock. The mechanical properties of interest include the rock permeability for the control of seepage flow, deformation characteristics of the rock mass, and strength of the material. In addition to reducing seepage, by reducing the permeability consolidation grouting can also reduce the uplift potential beneath the concrete due to water flow in the rock pores (Henn, 1996; Kikuchi et al., 1997).  According to Yesilnacar (2003) consolidation grouting around tunnels is performed in the regions determined to have weak rock conditions and leakage problems to provide the required stability of the rock outside the tunnel lining. This technique of grouting affects  15  the properties of the rock mass a minimum of one tunnel diameter beyond the excavation limit (Yesilnacar, 2003; Henn, 1996).  Consolidation grouting is usually performed by drilling and grouting shallow holes in a grid pattern. But if there are some local weak zones there might be some off pattern holes as well (Bruce, 2006).  Curtain Grouting Curtain grouting is the process of forming a water tight barrier by grouting in deep holes which is normally performed under the dams or in an abutment. The curtain seals seams, fissures, fault zones and cavities.  This type of grouting is generally used in structures that transport or store water or any other fluid, such as natural gas or petroleum (Henn, 1966).  In a typical curtain grouting hydraulic tests are carried out in pre-drilled boreholes to estimate the local in situ hydraulic conductivity. These tests may also be used to give an estimation on the amount of required grout in each borehole section. As discussed in the following chapters, the grout take of each fracture is related to the aperture of the fracture. In practice if the grout take is much lower than expected the grout is thinned to increase the grout take. The process of grouting will be continued until refusal (Houlsby, 1990).  16  2.3 Grout Types Based on Rheological Behavior and Material Grouts can be divided into three groups based on their mechanical properties (AFTES, 1991; Xanthakos et al., 1994). These three categories, listed in order of rheological performance (Figure 2.6) as classified by Xanthakos et al. (1994), are:  1. Particulate grouts: Grouts with a Bingham performance, which are suspension or cementitious grouts. Bingham fluid behavior is explained in Section 2.3.1.  2. Colloidal solutions: These grouts are silicate based chemical grouts which have an evolutive Newtonian behavior. As illustrated in Figure 2.6 an evolutive Newtonian fluid has an increasing viscosity over time. The silicate in the grout forms a gel which has properties that can seal fissures and improve the strength properties of the rock mass.  3. Pure solutions: Resin chemical grouts are non-evolutive Newtonian solutions. In these grouts the viscosity remains constant until the grout sets. This category of grout, which is made by solution of organic products in water, has the lowest viscosity of the three types.  17  Figure 2.6: Rheological behavior of typical grouts (Mongilardi and Tornaghi, 1986).  The difference between the rheological behaviors of a Newtonian and a Bingham fluid is also shown in Figure 2.7. This section is mainly focused on cementitious grouts, and the materials used to form cementitious grouts.  Figure 2.7: Basic rheological concepts (modified after Lombardi, 1985)  18  2.3.1 Bingham Grout Materials Bingham grouts are composed of cement, water, fillers (e.g. sand) and some additives. By changing the relative weight of these ingredients in the mix and also using different types of cement, fillers or additives, grouts with different properties can be produced. The most important properties of grout are:  •  Cohesion: Also called yield stress, cohesion is the maximum shear stress that can be applied to grout in the static condition. In other words, to make the grout move an applied shear stress larger than this value is needed (see Figure 2.7).  •  Viscosity: The coefficient of internal grout resistance to flow (see Figure 2.7). Thicker grouts have larger viscosities, while thinner grouts have smaller viscosities and flow more easily.  •  Stability: Bruce (2006) defined this property as the ability of grout suspension to exhibit little or no settlement of particles, bleed (separation of water as a result of grout particles separation), or shrinkage (changes in the volume due to excess use of cement in grout mix).  •  Set time: Also referred to as gel time, the set time is defined as the time after which there is a considerable increase in grout shear strength, or the liquid changes into a plastic consistency. This property can be affected by temperature.  Among the properties introduced here cohesion and viscosity are used in this study to model the flow of grout. Rombough (2006) performed a series of viscometer tests to measure the typical properties of Bingham grout material. Rombough’s tests included 19  three mixes with water:cement ratio of 0.8:1 using Type III, high early strength Portland cement. Two of these tests were performed on mixes with super-plasticizers. The other mix contained no admixture and was used as a baseline to be compared with other mixes (For the detail of tests see Rombough et al., 2006; Rombough 2006; Shuttle et al., 2007). Based on this study he suggested the ranges in Table 2.1 for cohesion and viscosity of a Bingham grout.  Table 2.1: Typical range of properties for Bingham grout (Rombough, 2006) Yield Stress (c) Pa.s 1-10  Viscosity(μ) Pa 0.011-0.016  Portland cement: Portland cement is the most popular cement used in cementitious grouts. This hydraulic cement is composed of hydraulic calcium silicate and hardens by the chemical process of “hydration”. Portland cement normally has a specific gravity of 3.15 gr/cm3 (Weaver, 1991).  ASTM Designation C150 “Standard Specification for Portland Cement” lists eight types of Portland cement:  •  Type I and II: Normal cement. A common cement for general purposes. Type II provides some degree of sulfate resistance.  •  Type III: high early strength and rapid set. Because it is a more finely ground cement, grout made with this type of cement has a much lower viscosity than types I and II. 20  The special properties of this cement make it more penetrable in fine fractures than types I and II. •  Types IA, IIA, IIIA: The same properties as types I, II and III plus air entering.  •  Type IV: low heat of hydration. This type of cement is used when the rise in temperature or heat generated by the hydration process has a negative effect on the surrounding conditions.  •  Type V: Pozzolnic Portland cement which is highly sulfate resistant. Used where the concrete is exposed to aggressive water.  Ultra fine Cement: Ultra fine cement, also called micro fine and super fine cement, is a Portland cement with an average particle size of 4 microns. This type of cement is used where very small openings have to be filled with a strong, durable, non polluting grout (Shimoda and Ohmori, 1982).  Having extremely fine particles (100% finer than 15microns (Weaver, 1991)) the grout made with ultra-fine cement is much less viscous than other types of grout. The viscosity of MC-500 (micro-fine) grout with a water cement ratio of 2:1 is reported to be almost half of viscosity of the same grout with normal Portland cement (Clarke, 1987).  Although ultra-fine cements are more expensive, their higher penetrability may reduce the number of boreholes required, and as a result may reduce the overall cost of grouting (Weaver, 1991). 21  To improve the properties of cementitious grout, such as corrosion resistance, stability, and also to reduce cost, different additives are often added to the grout mixture.  Additives: Additives are materials, other than water, cement and fillers, added to the grout mixture to modify the physical and chemical properties of the grout in a desired way. The type and amount of additives added to a grout depend on the type of admixture and the recommendations provided by the manufacturer. The most popular additives are described in this section.  Pozzolans: Pozzolans are silicates added to cement in order to increase resistance to chemical attacks by low PH and sulfate water. Different types of artificial or natural pozzolans are used in practice. The most popular ones are (Weaver, 1991): •  Blast Furnace Slag.  •  Condensed Silica Fume.  •  Natural Pozzolans (volcanic ash, volcanic glasses …).  Bentonite/Clay: Littlejohn and Bruce (1977) report that neat cement grouts are unstable with w/c ratios above about 0.4. Bentonite, a colloidal clay, is a very popular additive added to cementitious grouts to form a relatively stable suspension and increase grout stability. It  22  can absorb as much as five times its own weight in water (Brady and Clauser, 1986; Henn, 1996). By increasing the viscosity and cohesion of the grout, Bentonite decreases the settlement of cement particles in the mixture. Although Bentonite is not a lubricant, as was originally assumed (Deere and Lombardi, 1985), it limits the travel distance of the cement particles and as a result reduces the sedimentation or bleeding.  Disparents: Disparents, also called as flocculants, are used in grouts to reduce the tendency of the cement particles to agglomerate or flocculate. By reducing the viscosity and cohesion of the grout, this additive increases the ability of the grout to penetrate tiny fractures and openings. Although disparents improve the penetrability of the grout, they retard the set time significantly. A delayed set time might cause problems in projects with limited working area at a single time. It is advantageous in hot weather though, where high temperature shortens the set time (Weaver, 1991; Henn, 1996).  Accelerators: In order to have a faster grout set time, Type III Portland cement can be used, or an accelerator can be added to the grout. The most commonly used accelerators are calcium chloride, sodium chloride, sodium hydroxide and sodium silicate. •  Calcium chloride should be added in solution form as a part of the mixing water. Calcium chloride can increase shrinkage and may corrode reinforcement.  •  Sodium chloride is often used when calcium chloride is not available. This accelerator is much less effective than calcium chloride.  23  •  Sodium hydroxide is typically used when a very rapid set is needed. It is reported that by adding a 50% solution of sodium hydroxide to cement in a proportion of 2% by weight a set time of an hour or less would be achieved. The main disadvantage of sodium hydroxide is the extreme corrosivity, which limits the use of this chemical.  •  Sodium silicate is also a relatively fast acting accelerator. Because of its rapid effect it is normally added to the grout at the point of injection into the grout hole (Weaver, 1991; Henn, 1996).  Gas Producing Agents: These additives are basically finely divided metals such as zinc, aluminum and magnesium. Gas producing agents produce hydrogen gas through chemical reactions with alkalis. This property is used to reduce the effect of volume shrinkage after grout hardening (Weaver, 1991; Henn, 1996).  Fillers (Sand): Fillers are cheap materials added to a grout to reduce the amount of cement required, and hence reduce cost. Sand is the most common filler used in practice because of its availability. Sand also has the advantage of reducing the shrinkage of the ‘in place’ hardened grout. However, sand has some problems which limit its use as a filler. Sand particles are hard to hold in suspension and also because they do not dissolve in the mixture they can cause clogging in the pump and valves. To reduce this tendency more bentonite or fly ash, or a lower water cement ratio, is often used.  24  Sand is also abrasive, which might damage the grout distribution system (Henn, 1996; Weaver, 1991). In some cases clay is used instead of sand as a filler to form clay-cement grout. Clay-cement grouts were first developed in the Soviet Union, but this use has not been broadly documented in English technical journals (Weaver, 1991). Bozovic (1985) mentions that clay-cement grouts are not widely used because clay is not an industrial product. Not being an industrial material, clays have unknown and unpredictable properties which increase the risk of ruining grout properties (Weaver, 1991).  Water: Generally any drinkable water which has no odor or taste can be used in grouts.  2.4 Fracture Characterization The void geometry of fractures in hard rock is a complex three dimensional structure. The characteristics of this 3D geometry are mainly controlled by the chemical properties of the rock mass, stress conditions and geological history of the host rock. These three factors affect the pattern of the fractures by influencing the strength, weaknesses and sliding in the rock mass. The geometrical parameters of the fractures should be established in a way to make a quantitative description of the void space possible. Hakami (1995) introduces eight parameters which are believed to influence the mechanical and hydraulic behavior of the fracture. The primary advantage of these specified parameters is that they can be quantified by direct measuring. The eight  25  parameters introduced by Hakami (1995), and also shown in Figure 2.8, can be introduced as:  •  Aperture: The separation distance between the fracture surfaces. This is arguably the most important parameter in characterizing a fracture geometrically for grout filtration and/or distribution, and is described in detail in the next section.  •  Roughness: The shape of the surface which varies by the aperture distribution over the surface. As an example a fracture with a very high roughness has abrupt changes of aperture over surface.  •  Contact area: The area where the aperture becomes zero (or smaller than a specified threshold value chosen. Over the contact area the fracture surfaces are assumed to be in contact (and therefore negligible or no flow can pass through this area). The threshold value mentioned above depends on the definition for negligible value of flow and might vary for different works.  •  Matedness: How well the surfaces of the fracture are matched. If the matedness in a fracture is high the aperture variations over the surface will be small.  •  Spatial correlation: The rate of aperture change from one point to another.  •  Tortuosity: Changes in the direction of the streamline due to the fracture aperture variations.  •  Channeling: Change in the flow velocity due to the fracture aperture variations along a specific path.  •  Stiffness: Stiffness is a representation of the mechanical strength of the fracture. Stiffness is affected by factors like the stress conditions in the rock mass, mineralogy  26  or geological history of the rock. It can be studied by measuring changes in the aperture at a certain point over time. The aperture changes are normally a result of changes in the normal stress applied to the fracture.  All the above definitions are either directly or indirectly related to the fracture aperture and its variation. Therefore aperture and its variation is arguably the most important property of the fracture, and is a main focus of this research. Although aperture itself is dependent on the stress conditions of the fracture this effect is neglected here.  Figure 2.8: Fracture properties determined by fracture void geometry (from Olsson and Barton 2001).  27  2.4.1 Fracture Aperture The meaning of the aperture can be defined by the sketch shown in Figure 2.9. If the centerline between the fracture surfaces is assumed to be parallel to the x-y plane, as shown in the Figure 2.9, the aperture at each point is equal to the separation distance between the two fracture surfaces in the z direction.  Figure 2.9: Definition of aperture (Hakami and Larsson, 1996)  In order to quantify the opening size of the fractures two different aperture measurements, based on the measuring methodology, are widely utilized. These two apertures are the mechanical and hydraulic apertures.  Mechanical aperture: If two planes are drawn through the average of the highest and the lowest points along each of the surfaces of the fracture, then the perpendicular distance between the two  28  planes would be the mechanical aperture. In Figure 2.10 this aperture is denoted by em. Usually the mechanical aperture is geometrically measured from a two dimensional (2D) joint section (Olsson and Barton, 2001).  Figure 2.10: Illustration of mechanical aperture (modified from Evans et al, 1992)  Hydraulic Aperture: Since direct measure of the fracture apertures inside a rock mass is not always practical another aperture, called the hydraulic aperture, was established. The hydraulic aperture is an equivalent measure of the aperture which assumes the “cubic law” (defined below) is valid. Hydraulic aperture (eh) can be determined both from laboratory fluid flow experiments and borehole pump tests in the field (Olsson and Barton, 2001).  Measuring the hydraulic transmisivity of the fracture (T) in a flow test, the hydraulic aperture is obtained using the modified form of Darcy’s law (Cubic law). Cubic law can be derived assuming a 1D laminar flow of a Newtonian fluid in a plane parallel fracture (For the details of derivation see Fransson,1999 or Rasekh and Shuttle, 2006)  we 3 dP Q= 12 μ dl  Equation 2.1  29  In the above equation Q is the flow rate, w is the width of the flow and viscosity of the fluid.  μ is the dynamic  dP is the pressure drop over the length l. Transmissivity in terms dl  of flow rate and head loss would be:  T=  Equation 2.2  Q dP dy  Assuming a unit width of flow we would have:  e = 3 12μT  Equation 2.3  In real fracture systems the hydraulic aperture is smaller than mechanical aperture. The reason for the discrepancy between the two measurements is that smaller gaps have a disproportional effect on the head loss along a specific flow path; therefore the mechanical ‘average’ aperture is not the aperture controlling flow. There are a number of studies (e.g. Gale et al., 1990; Cook, 1992; Hakami, 1995) on the relation between mechanical and hydraulic aperture. However, hydraulic aperture is more popular because of the simplicity of its measuring methodology. Also hydraulic aperture is a representation of flow and since we are usually interested in flow behavior in the fractures, this type of aperture is used more than mechanical aperture.  30  A fracture can be represented using a frequency distribution of apertures. The frequency distribution gives the probability that points with a certain aperture occur on the fracture surface. Knowing the aperture magnitude at each point over the surface of the fracture the surface texture can be reproduced.  The frequency histogram of the aperture can be approximated by different mathematical probability functions such as normal, log normal or gamma distribution (Hakami and Larsson, 1996).  2.5 Fluid Flow in Fractures 2.5.1 Groutability Bruce (2006) defines groutability as the ability of soil or rock to accept grout. Based on this definition groutability of rock depends on the properties of the fracture and the injected grout.  For a joint to be groutable the grout must fulfill the conditions of penetrability and flowability. Penetrability means the ability of the grout to enter an opening and flowability means the ability of the grout to flow inside the fracture.  31  Penetrability: Penetrability is most important for cement based grouts where the grout particles might be too large to enter the fracture. Partial penetration of a fracture is termed limited penetrability and is associated with filtration or clogging. Some researchers use a limit aperture below which the grout cannot penetrate into the joint (Hansson 1995, Amadei 2000). By performing lab tests Hansson (1995) suggested that grouts with a d95 of less than 3 times of the hydraulic aperture have good penetrability. Eriksson and Stille (2004) introduce two aperture limits called bmin and bcritical. bmin is the lowest limit of aperture, under which no grout can enter the fracture. bcritical is a limit aperture above which there is an infinite penetration, meaning none of the grout particles will be stopped due to clogging. If the aperture is between these two limits there will be partial penetration and some filtration cakes will be formed. According to Eriksson and Stille (2004) these two limits depend on many factors, including time. Since the properties of the grout change over time, the size of apertures that grout can enter also varies over time. In order to be on the safe side, and ensure the aperture is bigger than the critical aperture, one simple solution is to use of micro-fine cement. But Eriksson et al. (1999) showed that there is a strong time dependency in the penetrability behavior of grouts based on a micro-fine cement (d95<12μm). Therefore, the use of micro-fine cement is not always an appropriate solution to poor penetrability.  Flowability: The other condition for groutability of a joint is the flowability of the grout inside the joint. Grout is flowable when it has sufficient pressure to overcome the resistant forces 32  and move forward in the fracture. Based on this definition flowability should be defined by studying the active forces on the injecting grout. One of the simplest solutions on the flowability of the grout is presented by Gustafson and Stille (1996) where the authors use a simple force balance for Bingham fluid flow to get the grout penetration. At refusal the injection pressure is balanced by the shear stress towards the fracture walls. If the fracture is assumed to have parallel walls with an aperture e and the flow is assumed to be one dimensional along the length of the fracture, the penetration, L, can be calculated by the following equation.  L = ( Pg − Pw ).e 2τ 0  Equation 2.4  In this equation Pg and Pw are grout and water pressure and  τ0  is the shear stress on the  walls. These forces are shown in Figure 2.11.  0  Pg  Pw  τ  0  L Figure 2.11: Active forces on grout in a slot  33  e  The main difference between the ideal case shown above and reality is the random changes of the aperture over the fracture surface. To have more realistic results Gustafson and Stille (1996) introduce an effective aperture which is an equivalent for the aperture used in the previous equation. Gustafson and Stille (1996) suggest that the effective aperture is close to the harmonic mean of the aperture values. The harmonic mean is defined as (e.g. Singh, K., 2007):  H=  Equation 2.5  n 1 1 1 + +L+ x1 x 2 xn  =  n  ∑  1 , xi > 0 for all i i =1 x i n  For some probabilistic distributions above equation can be simplified based on the properties of the distribution. As an example for a lognormal distribution the harmonic mean is calculated as:  ehm = Exp ( μ − σ 2 2)  Equation 2.6  Where ehm ,  μ and σ are the harmonic aperture, mean and standard deviation of the  lognormal distribution.  34  2.6 Methods to Estimate the Grout Distribution in the Ground The effectiveness of grouting is related to the distribution of the grout in the ground. This section discusses several methods to estimate the penetration of grouting based on the site investigations.  2.6.1 Pre and Post Grouting Tests Field tests in grouting projects before or after grout injection are quite popular. Performing pre-grouting tests helps the grouting design (i.e. deciding on the properties of grout and borehole arrangements). Results of testing performed after grouting can be compared with the pre-grouting test results to measure the effectiveness of grouting. Permeability tests are essential if a purpose of grouting is seepage control. Pre/Post grouting field tests, also called as water tests, can be for either: •  Exploratory testing  •  Grout hole testing  The aim of each kind of test is described later in this section. The results of water tests are commonly used to calculate Lugeon. This value is related to the crack conditions which can be used in choosing grout properties. In this section first the Lugeon concept is introduced and then the test methods.  35  Lugeon Unit Before the early 20th century dams were built in locations where there was no need to seal the foundation, and also these early dams were huge gravity structures under which the hydraulic gradient was less than one. In the early 20th century significant seepage developed through the subgrade of a few dams, including Janov Dam (1912 to 1914), after filling the basin. Subsequently methods to predict the in situ permeability of rock were developed.  One of the first and most practical evaluation systems was developed in 1933 by Swiss geologist Murice Lugeon, and has been used all over the world since its invention (Verfal, 1989). The idea of the method is to measure the quantity of water that can be forced out of a specific length of a drill hole in a unit of time under a set pressure. The amount of water is defined as a Lugeon unit. One Lugeon is 1 L of water per meter of the hole length per minute at a pressure of 10 bar (1 MPa). Based on the results criteria were proposed for admissible water loss, including those proposed by Lugeon, Jahde and Terzaghi (Verfal, 1989; Warner, 2004).  Exploratory testing This type of water testing is performed to measure the permeability before, during or after grouting. Exploratory tests are normally performed in core drilling holes and their results are used for grout design. In the cases where no coring is required some water testing holes are drilled.  36  Houlsby (1990) introduces this water test method as five steps of applying water with different pressures right after each other without any pause in between: •  Low pressure for ten minutes  •  Moderate pressure for the next ten minute  •  Peak pressure for next ten minutes  •  Back to Moderate pressure for next ten minutes  •  Low pressure for next ten minutes  Then the Lugeon values are calculated for each of the steps. For the Low, Moderate and Peak pressure Houlsby (1990) gives the values in Figure 2.12. As illustrated in this figure the 10 bar pressure given by Lugeon is the maximum peak pressure here. Figure 2.12 gives 10bar as the maximum pressure (used in third increment of the test) for depths larger than 44m. The reason that Lugeon used higher pressure is that the original water test, specified by Lugeon, was in the water wells while the current tests are performed in core drilled holes and 10 bars is a high pressure for this situation (Warner, 2004).  37  Figure 2.12: Pressures for each run of the five pressure water test (after Houlsby, 1990)  38  Grout hole testing This type of water testing is performed in the grouting boreholes right before the start of grouting at each stage to: •  Estimate the likely grout take.  •  Determine effect of grout penetration from previous stages.  •  Check the values assumed in the design process.  •  Check for surface leaks, connections between holes and any other unpredicted happenings during grouting.  One of the methods of performing this test is to inject water to the hole by a specific pressure. 5, 10, 15 minutes after reaching the defined pressure the amount of water entering the hole per unit time is recorded. Houlsby (1990) recommends 1 bar for the pressure in this test.  There are lots of recommendations relating the rock mass permeability from water pressure tests to the need for grouting. For example, Lugeon proposed, based on the original Lugeon test with pressure equal to 10 bars, that if the measured Lugeon is less than the following values there is no need for grouting (Verfal, 1989): •  1 L/min/m  for dams higher that 30m  •  3 L/min/m  for dams lower that 30m  Houlsby (1990) also suggests that based on the proposed exploratory testing methodology given above where water is precious and a very low permeability is desired  39  the closure criteria for the grouting should be in the range of 1 to 3 Lugeons. But he does not give specific values for what he means by low permeability.  Although the results of water tests are quite useful in providing a closure criteria for grouting, these results cannot be used to predict the grout penetration in fractured rock. Lombardi (1985) and Ewert (1997) discuss the reasons that the results of the water test are often not reliable as grout take estimations: •  Grout is a Bingham fluid while water is a Newtonian one.  •  Grout is a suspension fluid with particles that may settle during grouting causing the grout properties to vary with time.  •  Due to the size of grout particles, grout cannot penetrate all of the fissures that water enters easily.  •  Grout has ‘cohesion’ and will arrive at a maximum penetration distance and stop if a constant pressure is applied. But under the same conditions water with no cohesion will flow continuously.  2.6.2 Modeling Modeling is another approach in prediction of the grout distribution in fractures. To model grouting of a fractured rock, two media of grout and fracture should be modeled. Due to uncertainties and complexities in the modeling of both media simulating the grout spread in a fissured media is a challenge. Development of a model that is simple enough and does not require a considerable amount of computational resources is accompanied by making some simplifying assumptions. Each of the available numerical models has 40  some simplifying assumptions which make the calculations easier. As long as these assumptions do not result in unrealistic conditions, a numerical model can be used as a substitute to the experimental and/or empirical methods. Such models are powerful tools in gaining a better understanding of the penetration behavior and investigation of the effect of grout properties. The best known methods of modeling grouting in fractured rocks, as well as their simplifying assumptions, are reviewed and discussed in this section. The discussions in this section will provide grounds for the choice of methods in the present study. These studies consider simple fracture geometry (e.g. constant aperture) with complex grout properties (e.g. Bingham solution) or simple grout properties (e.g. considering only cohesion) with complex geometry (e.g. stochastic aperture changes). Grout Flow Simulation One of the complexities in modeling of grout penetration in a fractured rock is associated with the complicated behavior of grout. All researchers agree on the Non-Newtonian behavior of cementitious grouts (see Section 2.3 for a more detailed review). Most of the previous studies indicate that the grout flow is governed by Bingham fluid equations (e.g. Hassler et al., 1992; Fransson, 1999; Shuttle and Glynn, 2003). Section 2.3 explained that Bingham fluids can be simulated using a dynamic viscosity and cohesion approximation. Gustafson and Stille (1996) calculate grout penetration based on a force equilibrium equation only considering cohesion for Bingham fluid and ignoring viscosity (for more details see Section 2.5.1). Fransson (1999), Hassler et al. (1992) and Eriksson et al. (2000) use the same approach in predicting grout spread in a fracture. However, some researchers (e.g. Hakansson, 1993) believe that the Bingham model with constant 41  properties does not give an accurate description of cement based grout behavior. To improve the simulation Hassler et al. (1992) and Eriksson et al. (2000) used a Bingham model with time varying grout cohesion. Hassler et al. (1992) showed the significant effect of grout hardening on the penetration process for a long period of grouting. Their results show that considering the time dependency of grout cohesion can change penetration length results up to 30% after 2hrs of grouting. However, in shorter injection periods hardening does not affect the results significantly. Therefore, it can be concluded that assuming constant properties for the grout as a Bingham fluid is a reasonable assumption in a short grouting period.  Although the studies discussed above give reasonable results, they are based on the modeling Bingham grout behavior as a cohesion and do not consider viscosity in their simulations. To model Bingham grout completely, both viscosity and cohesion should be considered in the simulations. The complete solution, discussed in Chapter 3, has been adopted in recent works including Shuttle et al. (2007) and Rasekh and Shuttle (2006) to predict grout penetration in a single fracture.  Fracture Simulation Fractured rocks are generally simulated at two different scales. The first scale considers the whole rock mass as a grouting media. The properties of this media are defined by the fractures density (number of fractures in a specific space), pattern (e.g. parallel or perpendicular fractures) and fracture spacing. The second, smaller, scale considers a single fracture in rock, and focuses on the geometrical changes over the fracture surface. 42  The use of each approach in different studies is reviewed. Since this research is modeling the grout spread in a single fracture this chapter is mainly focused on the methods of modeling a single fracture.  Fracture networks in rock masses: National Research Council (1996) classifies mathematical models of fractured geological media into three groups: •  Equivalent continuum models  •  Discrete fracture network models  •  Hybrid models  Equivalent continuum models divide the rock mass as blocks in the scale of interest and represent the properties of each block using equivalent properties (e.g. equivalent permeability or effective porosity). These properties, defined by deterministic or stochastic methods, are used to model fluid flow in the rock mass (NRC, 1996). The main attraction of equivalent continuum models is their simplicity. The complexity of flow patterns in random fracture network are avoided, as instead of going into the details of flow pattern, continuum models predict fluid flow behavior based on the equivalent properties of the rock mass. These properties are the inputs for the models. However, defining the input parameters for equivalent continuum models is difficult. To find a reasonable coefficient for rock mass the correct scale should be chosen to represent the rock mass in terms of uniform properties. The volume of rock needed for averaging, defined by REV (Representative Element Volume), should be larger in the cases with  43  smaller fractured densities. In such cases REV is bigger than the area covered by field measurements (Eriksson, 2002). Oda et al., 1987, recommend using this method for highly fractured rocks.  The other method to simulate fractured rock masses is Discrete Fracture Network (DFN). In this method fractures are represented by polygons in 3D space (Figure 2.13).  Figure 2.13: Fracture network model by DFN method (Shuttle and Glynn, 2003)  In most large fracture systems, the DFN is constructed using a combination of deterministic fractures, conditioned fractures and purely stochastic fractures (Shuttle and Glynn, 2003). A deterministic fracture is a fracture with known properties, a conditioned fracture a fracture with partly known properties, and purely stochastic fractures are the ones with unknown properties and are generated based on statistical functions.  44  Due to the complexity of large DFN models, typically the grout flow component of the model is quite simple. For example, Shuttle and Jefferies in The Practical Handbook of Grouting (Warner, 2004) provide examples where the grout flow in fractures is modeled without solving the flow equations.  The two models discussed here differ in representation of the heterogeneity of the fractured rock. Both of the methods may use deterministic or stochastic properties to simulate a fractured rock mass. In comparison with the equivalent continuum models, DFN models are more realistic because volume averaging approximations is avoided in the fracture network scale. The third type of simulations adopts a combination of the first two approaches by using discrete network models in building continuum approximations (NRC, 1996).  Modeling a Single Fracture: Single fracture models of grout penetration are used because: •  Due to the smaller scale, it is possible to simulate the aperture variation and small scale roughness of the fracture surface. These variations control the grout penetration in a fracture (Eriksson et al., 2000).  •  At smaller scales, more realistic models of grout flow behavior are practically computable, such as the solution of Bingham fluid used by Shuttle et al. (2007), and which improve the accuracy of the grout penetration prediction.  45  Some of the models simulating grout spread in a single fracture are reviewed in this section.  Shuttle et al. (2007) is a recent example of modeling grouting in a single fracture. They assume constant aperture over the surface of the fracture. The model adopts Bingham fluid solution for a cylindrical grout spread. To consider the aperture variations in a fracture, Fransson (1999) and Rasekh and Shuttle (2006) adopt probabilistic functions to generate the aperture distribution over the surface of fracture. Rasekh and Shuttle (2006) use a finite element based software (Comsol) to model grout penetration in the fracture with stochastic properties. This work considers variations of aperture over the surface of fracture as a log normal distribution. By solving Bingham equations for 2D elements the grout penetration in stochastic fractures is modeled. Solving 2D Bingham equations in small elements in fractures with different apertures gives precise predictions of fluid propagation in fractures. However, the size of fracture and number of elements was limited due to calculation capacity problems. Similar to Rasekh and Shuttle (2006), Fransson (1999) also simulated aperture variation with a lognormal distribution, but as described earlier in this chapter, used a simpler ‘cohesion only’ grout model.  Solving the full Bingham equations in 2-D is computationally intensive, which places limitations on the number of elements in the model (e.g. Rasekh and Shuttle, 2006). To reduce the computational effort, some researchers (Hassler et al., 1992; Eriksson et al., 2000) have simplified the elements by substituting them with one dimensional channels (which have a simple analytical solution for the Bingham equations) This simplification  46  builds upon the observation that fracture flow is reported to follow a limited number of paths (Eriksson et al., 2000), so the 1-D flow elements are used to represent these flow paths within a fracture. Using 1D channels to model the grout penetration has two main advantages in calculations from a simplicity point of view:  •  Use of 1D equations: If 2D elements are substituted by narrow channels, 1D equations can be used. Using 1D equations which are simpler than 2D ones, helps increasing the modeling capacity (i.e. being able to model larger fractures).  •  Modeling more than one fracture using the same algorithm: By changing the orientation of 1D channels in 3D space, more than one fracture crossing each other can be modeled.  Hassler et al. (1992) and Eriksson et al. (2000) use force equilibrium equation in the 1D channels to predict the grout penetration. Similar to Fransson (1999) their models only consider cohesion of grout.  47  3 THEORETICAL MODEL FOR ONE DIMENSIONAL FRACTURE FLOW This chapter gives the details of the theory and equations used for the grout flow in a single fracture in the current work.  To solve the fluid flow three main assumptions are made: •  The grout flow is Laminar: Laminar flow has very low velocity and no turbulence. Grout velocity is generally less than 1m/min. Also, grout viscosity is typically around 10mPa.s (almost 10 times the viscosity of water). Such conditions indicate Reynolds number less than 100 for the flow which is in the range of laminar flow.  •  Grout is an incompressible fluid: The amount of fluid entering a controlled volume is equal to the amount leaving that volume. During the injection time that the grout is flowing this is a valid assumption.  •  Grout is a Bingham fluid with constant properties over the grouting time: Cementitious grouts are regarded as Bingham fluids. Before grout set time the properties of grout can be assumed as constant. Therefore, in modeling short periods of grouting as in this work this is a valid assumption.  The aim of solving the grout penetration in a network of channels is to calculate the unknowns in the system which are penetration length and also the grout flow rate in each channel. So the number of unknowns for each network is twice the number of pipes. To solve the system for these unknowns, generally two sets of equations can be used: •  Continuity 48  •  Force equilibrium  Continuity equations: The basis of continuity equations in fluid mechanics is the conservation of mass in the control volume. This means that for an incompressible fluid the mass volume entering the control volume is equal to mass volume leaving the control volume.  In a pipe network if a control volume is assumed at each of the nodes, the amount of fluid entering the control volume should be equal to the outgoing grout.  Q1  Qi  Q2 Q3  Figure 3.1: Schematic display of the inflow and outflow in a node  Figure 3.1 shows a typical node in a pipe network. In this case the continuity equation will be:  Qi = Q1 + Q2 + Q3  Equation 3.1  To solve the whole network these equations should be satisfied for all of the available nodes.  49  Force equilibrium equations: Figure 3.2 shows the forces acting on an element of fluid inside a channel. Writing the force equilibrium equation for such an element gives:  [P − ( P + ΔP)]⋅ y ⋅ w = τ yz ⋅ L ⋅ w  Equation 3.2  ⇒−  ΔP ⋅ y = τ yz L  Where w , L and y are the width, length and height of the element respectively. P and  P + ΔP are the pressures applying at each of the faces of the element. τ yz is the shear stress in the fluid at the height equal to y .  w  τ B  y=b  y  yz  P  C  P+ΔP  y  z x A  D  L  50  Figure 3.2: Forces in a flowing fluid inside a channel (modified from Chhabra and Richardson, 1999).  It was shown in Figure 2.7 that in a Bingham fluid the shear stress is related to the shear strain rate by the viscosity and cohesion or yield stress of the fluid.  The mathematical format of this definition is:  Equation 3.3  In the above equations  and  μB  τ yz = − μ B  du + τ 0B dy  for  τ yz ≥ τ 0B  τ yz = τ 0B (  du = 0) dy  for  τ yz ≤ τ 0B  is the dynamic viscosity,  du is the rate of the shear strain dy  τ 0B is the yield stress.  Substituting  τ yz by its mathematical definition for a Bingham fluid (Equation 3.3) in  Equation 3.2, the fluid velocity can be calculated. Integring velocity over aperture gives the flow rate per unit width, “q”, as:  Equation 3.4  q=−  Where  φ=−  51  1 2ΔP 3 ⎛ 1 3 3 ⎞ b ⎜1 + φ − φ ⎟ 2 ⎠ μ B 3L ⎝ 2  τ 0B ΔP L ⋅ b  b is the half aperture indicated in Figure 3.2.  For an element with a length tending to zero, and width w, the total flow rate is given by:  Q=−  Equation 3.5  2w dP 3 ⎛ 1 3 ⎞ b ⎜1 + φ 3 − φ ⎟ 3μ B dl ⎝ 2 2 ⎠  Equation 3.4 and Equation 3.5 are gives as standard equations for flow of a Bingham fluid through a smooth parallel sided aperture by Chhabra and Richardson (1999). For the case where  τ 0B = 0 the Bingham fluid will be changed into a Newtonian one. In  such a case Equation 3.5 will be identical to Equation 2.1 (i.e. the cubic law for a Newtonian fluid). If Equation 3.4 is rearranged it can be written as a standard cubic equation for  ⎛  φ 3 + ⎜⎜ − 3 − 3  Equation 3.6  Solving this equation for  ⎝  φ.  q.μ B ⎞ ⎟φ + 2 = 0 τ 0B .b 2 ⎟⎠  φ is standard and gives the relation between the pressure  gradient and the flow rate in a pipe.  θ +π ⎛ −α ⎞ ) ⎟ . cos( 3 ⎝ 3 ⎠  φ=2 ⎜  Equation 3.7  52  ⎛ 27 ⎞ cos(θ ) = ⎜ 3 ⎟ ⎝ −α ⎠ ⎛ ⎛ 27 ⎞ ⎞ ⎟ 3 ⎟⎟ α − ⎠⎠ ⎝ ⎝  Where:  θ = cos −1 ⎜ ⎜ ⎜  ⎛ q.μ ⎞ α = ⎜⎜ − 3 − 3 B B2 ⎟⎟ τ .b ⎝  (φ = −  ⎠  0  τ 0B ΔP L ⋅ b  )  The detailed derivations for all of the equations in this section are given in Appendix A. From Equation 3.7 if either pressure drop or flow rate is known the other can be calculated.  53  4 NUMERICAL MODELING OF GROUT PENETRATION IN A FRACTURE 4.1 Objectives of the Model This study’s focus is the estimation of the penetration and distribution of grout inside a rock joint; considering different grout types and fracture properties. Following the aim of the research, a numerical model is required that is able to predict the propagation pattern and penetration length of the grout. This model should also be able to calculate changes of pressure over the fracture surface during grouting. The following sections describe how these objectives have been obtained. The model introduced in this chapter simulates a network of pipes with randomly distributed apertures representing the fracture surface with stochastic properties. Then the model uses the Bingham solution given in Chapter 3 to compute the grout spread in this stochastic fracture. This chapter describes how the model solves the Bingham 1D equation in the stochastic pipe network.  4.2 Simplifying Assumptions in Simulation In order to provide a practical and simple model simplifying assumptions are made in the simulation process. This section describes the assumptions made in the model in simulation of grout properties, fracture geometry, flow behavior and boundary conditions of fracture.  54  4.2.1 Grout Property Assumptions  •  Particulate cementitous grout with a Bingham behavior. Cementitious grout is the most popular grout in rock grouting, and this type of grout is a Bingham fluid, therefore assuming grout is a Bingham fluid is a reasonable assumption.  •  Constant properties of the grout over grouting time: As discussed in Section 2.3, the rhelogical properties of the grout, for example viscosity, typically change over time. In the current model all of the properties of the grout are assumed to be constant over the grouting time. As discussed in Section 2.6.2, this is considered a reasonable assumption for the short grouting times considered where grout hardening does not have a significant effect on grout behavior.  •  All of the apertures are groutable: In the simulations presented here the apertures are chosen in the range of millimeters. Therefore, we can assume that the apertures are all above the critical limit defined in Section 2.5.1 and hence injectable.  4.2.2 Fracture Geometry Assumptions •  Constant aperture over the surface of each element:  Most of the geometrical  properties of a fracture are related to the aperture changes over the surface of the fracture (see Section 2.4). In the current work the fracture is divided into a number of elements each of which has a different aperture. The random variations of these  55  apertures over the surface of the fracture represent the varying texture of the fracture in a simplified way. The aperture is assumed to be constant over the surface of each of these elements.  •  No change in the aperture resulting from the grout pressure: Depending on the strength of the rock and grouting pressure, there might be some changes in the fracture openings during grouting which are neglected here.  4.2.3 Grout Flow Assumptions •  Ignoring the water effect on the penetration surface of the grout: Many grouting projects, for example grouting to seal fractures under dams, are performed below the water table. In such cases fractures are filled with water. In this model the pressure at each point is assumed to be the grout pressure at that point relative to the water pressure. So the zero pressure at the grouting face means that the grout pressure is equal to the water pressure.  4.2.4 Boundary Conditions Assumptions •  Calculations are based on a constant injection flow rate: The input boundary condition for grouting is the constant grout injection flow rate, and the changes of pressure with time are calculated.  56  •  ‘No flow’ outer boundaries: The outer boundaries of the model are assumed to be zero flow. This means any pipe intersecting the boundary is a dead end.  4.2.5 Mesh Generation Assumptions •  Use of 1D Bingham flow solution: The current research extends the work of Fransson (1999) who represented Bingham flow using only cohesion, by considering both the cohesion and viscosity of the grout in the Bingham idealization. The solution adopted here is for a 1D parallel walled channel.  •  Use of a “Finger Model” as the one dimensional pipe network: In this context a Finger Model means a network of pipes with flow in only one direction through each of the pipes. Considering N pipes connected to one node, only one pipe provides the inflow to the node, while N-1 of the pipes have grout outflow from the pipe.  1  QO  QO2  Qi  QO  3  Figure 4.1: Schematic display of finger model  57  4.3 Algorithm The current model is composed of two main codes which are the fracture generator and the solver. These two codes work in a sequence where the output of the first code (generator) is the input for the next code (solver). The algorithm of each of these codes is described in this section. Also, for more detail the reader is referred to Appendix C and Appendix D which contain a hardcopy of each of the codes.  4.3.1 Random Fracture Generator In Section 2.4 parameters that characterize a joint in the rock such as aperture, aperture distribution, roughness and contact area were introduced. Among these parameters aperture and the aperture distribution were identified as the most useful ones to represent a fracture. This section gives the details on the methodology adapted to generate a random fracture by taking aperture and aperture distribution as inputs.  The generator developed for this work is a Visual Basic Application (VBA) code. It is capable of modeling fractures with assigned dimensions, number of elements and aperture distribution.  A normal distribution was chosen because of the following  advantages:  •  Realistic results: This probability function was proven to be one of the best functions to model the real aperture variations over the surface of a fracture (see review of functions considered by Hakami and Larsson, 1996).  58  •  Simplicity in defining the input parameters of the distribution: It will be explained later in this section that knowing the mean and minimum or maximum aperture a normal distribution can be generated. Since these properties are easy to define for a fracture (in comparison with standard deviation), this distribution is more user friendly.  The normal distribution has the general form shown in Figure 4.2 (e.g. Galambos and Kotz, 1978). In this figure σ is the standard deviation of distribution. The probability of a value 3σ less than mean in a normal distribution is 0.1%. Therefore, mean minus 3σ may be regarded as an approximation for the minimum value (mean - 3σ ) in a normal distribution. This approximation is much more valid as the number of readings increases, for a limited data it is only a rough approximation. Since the shape of the normal distribution function is symmetric the same argument is valid for a value 3σ larger than  Probability Density Function  mean, defining the maximum value.  minimum  3σ  3σ  mean( μ )  Maximum  Value of Variable Figure 4.2: General form of normal distribution (Galambos and Kotz, 1978). 59  Also based on the mathematical function for probability density function of normal distribution probability of any value in a normal distribution can be calculated. As an example probability of a value 1σ less than mean is  24  σ  % (Equation 4.1). Therefore, for  a given set of numbers if the probability of each value matches the result from Equation 4.1it means the given set has a normal distribution. ⎛ (x − μ)2 ⎞ 1 ⎟ PDF = exp⎜⎜ − 2σ 2 ⎟⎠ σ 2π ⎝  Equation 4.1  1 ) − exp( ⎛ (( μ − σ ) − μ ) ⎞ 0.24 1 2 ⎟⎟ = = exp⎜⎜ − PDF = 2 σ 2σ σ 2π ⎠ σ 2π ⎝ 2  Considering the mathematical properties of a normal distribution, if the mean and standard deviation of the distribution is known the minimum and maximum apertures can be calculated. Also by knowing mean and minimum or maximum aperture the standard deviation can be found.  Equation 4.2  σ =  mean − min 3  σ=  max − mean 3  Figure 4.3 shows an example of the configuration of an input interface for the generator. For the normal distribution chosen in this example and the values shown the minimum and maximum aperture are: 60  •  Maximum aperture: 1mm + 3(0.27mm)= 1.81mm  •  Minimum aperture: 1mm - 3(0.27mm)= 0.19mm  Figure 4.3: Sample input for the fracture generator.  In addition to the aperture distribution properties, the dimensions (length and width) and element properties (number of elements along each side and aperture of each element) of fractures are also defined as input parameters of the generator (Figure 4.3). Based on this input information, the code generates two outputs. The first one includes the properties of all 2D elements on the fracture surface as requested in the input. This output, which contains the location, dimensions and aperture of each 2D element, gives an insight into the two-dimensional variation of the aperture over a 2D plane. However, the information in this output file cannot directly be transferred to the solver since flow equations are 61  solved for 1D channels. The second output file contains properties of equivalent 1D elements, including their length, aperture and width. As will be described in the following section, the second output is an input for the solver.  Figure 4.4 shows the output for the 2D elements, for the input properties given in Figure 4.3. The output file has been imported in AutoCAD and the figure is generated using this software. Different grey shades indicate the magnitude of the aperture in each of the elements. According to the input data the input mean and standard deviation are 1.0mm and 0.27mm, respectively. Therefore, ideally the minimum and maximum expected values are mean- 3σ =0.2 mm, and mean+ 3σ =1.8 mm. However, statistically these values are obtained only for an infinite number of elements. In the shown example with 150 elements the minimum and maximum values generated (0.29 to 1.78) are reasonably close to the theoretically predicted values. Although a larger number of elements will decrease the difference between requested properties and output values, the more computational time it imposes in solving the equations makes it an improper option regarding the objectives of the modeling.  62  0mm  1mm  2mm  Figure 4.4: Sample output of the 2D elements of the fracture generator.  In order to be capable of solving large systems of fractures, and also keep the model as simple and computationally efficient as possible, 1D elements has been chosen for the network. As explained above the second output of the generator is a 1D pipe network which is an equivalent for the initially generated 2D mesh shown in Figure 4.3. Below the process of converting a randomly generated fracture aperture network of 2D elements to an equivalent network of 1D elements is explained in more detail.  Figure 4.5 shows a typical 2D network and its equivalent 1D network. In this network each of the 1D equivalent pipes is located in between two adjacent elements connecting centers of those two elements (see the zoomed element in Figure 4.5). 63  Figure 4.5: Replacing 2D elements with equivalent 1D pipe network.  Figure 4.6 shows the two adjacent 2D elements from the previous figure along with their geometrical properties. In this figure w, l and e are the width, length and aperture of the element. Equivalent properties of each horizontal/vertical 1D element are found based on the properties of left/bottom 2D element the pipe is passing through. For the given 2D elements properties in Figure 4.6, properties of the equivalent 1D element are found as:  Equation 4.3  weq =  w1 2  leq = l1  and  eeq = e1  In the previous equations w, l and e are the width, length and aperture of the 1D element. The subscript “eq” stands for the equivalent property and physically shows the area that the 1D element represents. The reason half widths are used in the averaging equation is as 64  follows. In the mesh each 2D element is substituted by two one dimensional equivalent elements, one horizontal and one vertical (Figure 4.7). If the full width is considered in averaging, the total area covered by 1D elements will be two times the area of 2D elements. To avoid this, half width of 2D elements is used in computing the properties of 1D elements.  e1  e2  w1  w2  1  2  Figure 4.6: Averaging properties of adjacent elements.  1D Equivalent Elements  Figure 4.7: 1D equivalent elements 65  These 1D elements are used in the solver as channels in which grout flows. Therefore, flow equations are solved only in one dimension, saving a considerable amount of computation. The process of modeling grout in these elements is described in the following section.  4.3.2 Solver The solver, which is the main part of the code, solves the Bingham laminar flow equations for a network of pipes. To solve the Bingham fluid flow in the pipe network the system in which the grout is flowing should be defined and updated during the grouting. Then the governing equations should be solved for the defined network. The code uses the following steps to achieve these tasks:  •  Creating and updating the network at each time step: Creating the initial network at the first time step or adding new pipes to the end of the filled ones and deleting the dead pipes at later time steps.  •  Building up the matrix of equations and solving the governing equations for the whole network.  •  Updating the new conditions in each of the pipes.  Followed each of these tasks will be explained in more detail.  66  Updating the network  After generating the pipe network using the mesh generator (described in Section 4.3.1), the node to be used as the “injection point”, where the grout enters the fracture, is specified by the user. The solver reads the data for all of the pipes provided by the generator. These data include the location, width, length and aperture of each pipe. The initial mesh comprises only the pipes directly connected to the injection point. In the example illustrated in Figure 4.8 this initial network comprises 4 pipes.  Figure 4.8: The pipe network in the first time step.  At the end of each calculation time step, the penetration distance in all of the pipes in the current network is checked. If any of the pipes is full, pipes attached to the end of the full pipe are added to the system. If the pipe is a dead end and no further pipes are attached, the pipe is reported as “dead” and it is deleted from the network. Figure 4.9 shows the network from Figure 4.8 at time step n (n > 2) where additional pipes are added when a pipe becomes filled, and at time n+1 where the upper pipe is filled but no further pipes are attached (i.e. at a dead end). 67  Figure 4.9: Network update in each time step.  This pipe network update is done at the beginning of each time step, so that in the calculations the equations for the new pipes are added to the system of equations and the flow is solved using the new network.  Solving the network  In Chapter 3 a detailed description of the governing equations for the laminar flow of a Bingham fluid was given. The solver builds up a matrix of equations based on these governing equations. The equations of a network include three sets of equations:  •  Continuity equation in each node.  •  Continuity equation in each pipe.  •  Equal head drop in the branches attached to a single node. 68  The system of governing equations which includes nonlinear equations should be solved by an iterative method. The details of the adopted iterative method are presented in Appendix B.  The first two sets of equations are the continuity equations in each pipe as described in Chapter 3. The third set of equations is based on the known head drops in each branch. We know that the pressure at the penetration face (defined as the grout-water/air interface) is equal to zero, or in the cases where grouting is under water table the pressure at this point is equal to the water pressure. So starting from a node like injection node (shown in Figure 4.10 by a black shaded circle) and following the grout path in each of the branches to the interface with air (or water) the head drop should be equal in all of the routes starting from the same node (Figure 4.10).  Figure 4.10: Pressure drop in the pipe network. 69  In Figure 4.10 the pressure at the injection point is defined as P0 and the pressure at the penetration surface in each of the branches is zero. Therefore, the total head drop is:  ΔPtot = P0 − 0  Equation 4.4  Considering two different routes starting from the same node (see Figure 4.10):  ΔPtot = ΔPa + ΔPb  Equation 4.5  ΔPtot = ΔPc + ΔPd + ΔPe  Therefore, the equal head drop in different branches can be written as:  ΔPa + ΔPb = ΔPc + ΔPd + ΔPe  Equation 4.6  Knowing the flow rate in each pipe, ΔP in each of the pipes can be calculated using Equation 3.7.  70  The equations for a general network shown in Figure 4.11 (where the letters show the node labels) is:  Figure 4.11: A general form of pipe network.  1. Continuity equation in each grid:  Qij =  Equation 4.7  m + N −1  ∑Q  k =m  jk  , (N=number of branches at node j)  Where Qij indicated the flow rate in the pipe between grid i and grid j.  2. Continuity equation in each pipe:  ltjk = l ( t − Δt ) jk +  Equation 4.8 If ltjs > l js ( m ≤ s ≤ m + N − 1 ) then  71  ltjs = l js  Q jk A jk  .Δt  Where ltjk is the penetration length at time step t in the pipe connecting nodes j and k. Also l jk is the length of the pipe between nodes j and k.  3. Equal head drop in the branches attached to a single grid:  Equation 4.9  ΔP = f (Q jm ).ltjm = f (Q jm +1 ).ltjm +1 = ... = f (Q jm + n ).ltjm + n  If ltjs > l js ( m ≤ s ≤ m + N − 1 ) then:  ΔP = f (Q js ).l js + f (Q sr ).l tsr = f (Q jk ).l tjk ( m ≤ k ≤ m + N and k ≠ s )  For a Bingham Fluid above parameters are:  f (Q) =  Equation 4.10  φ =2  ΔP 2τ B =− 0 l 2b.φ  −α 1 . cos( (cos −1 ( 27 )+π) −α 3 3 3 α = −3 −  12μ .Q 4wb 2τ 0B  These three sets of equations show that for a network with “n” pipes there will be “2n” unknowns and “2n” equations. Some of these equations are not linear or simple ones. Therefore, this is a system of 2n nonlinear equations with 2n unknowns which has to be solved by iterative methods. The requirement to use an iterative method to solve the system is one reason the solution is complicated. There are a number of iterative methods  72  available to solve such a system of equations. In the current work Newton’s method was used. The details of the Newton method and the solution to this system of non-linear equations are presented in Appendix B. In this iterative method a value is assumed for the flow rate and penetration length in each pipe of the network at the beginning of each time step. These flow rates and penetration lengths are improved in iterations. A good initial value saves the number of iterations and running time. In the current code values from the previous time step are considered as the initial assumptions. In small time intervals the changes of grout flow (i.e. flow rate and penetration length) in each pipe is so small that the previous time step values as the initial input for the iterative solution are close to the new results.  Another difficulty in solving grout penetration, is building head drop equations for the whole system. Equation 4.9 shows that when further branches are added to the end of a pipe these new branches should be considered in the head drop equations for all the nodes connected directly or indirectly to the new pipes. In other words, to write the head drop equations all of the grout routes from each node should be considered. Since these routes change during grouting time by adding or deleting pipes, they should be checked and updated in each time step.  Updating the conditions in each of the pipes for the calculated values  In this stage the pressure, flow rate and penetration length in each of the pipes has been calculated. So the values from the previous time step should be replaced by these new  73  values and also the penetration lengths should be corrected if they are larger than the length of the pipe.  Finally the developed solver gives an output at each time step which includes the flow rate, penetration length and pressure in each of the pipes. The main advantage of this solver is that it can model the grout penetration as a Bingham fluid over a fracture with stochastic properties and reasonable dimensions. As explained in Section 2.6.2, most of the previous models consider either the stochastic fracture properties (e.g. Fransson, 1999) or Bingham fluid solution (Shuttle et al., 2007). The ones which consider a realistic model for both components at the same time (e.g. Rasekh and Shuttle, 2006) have calculation capacity problems which make these models only able to simulate small fractures in short grouting times.  74  5 VERIFICATION CASES 5.1 Overview of the Chapter To verify the correctness of calculations and suitability of the model, four simple verification cases were performed. The first verification case examines the generator results by comparing the properties of generated set with the requested properties in the input. The other cases check the solver outputs. As explained in Section 4.3.2, the solver solves three sets of equations to compute the penetration and pressure distribution in the fracture. Each verification case is intended to verify the solution of at least one of these sets of equations. The solver results are compared against independent hand calculations. The last case is brought to verify the code for a Newtonian fluid like water which is in fact a Bingham fluid with zero cohesion.  5.2 Verification Case I: Fracture Generator To check the generator results the distribution type, mean and standard deviation of the apertures given by the generator should be compared with the properties requested in the input. Figure 5.1 shows the input for the generator, showing the mean aperture of fracture and standard deviation of normal distribution as three and one respectively. The standard deviation for a normal distribution is defined as: Standard deviation = (mean aperture – min aperture)/3. According to this definition for the input properties shown in Figure 5.1 the minimum and maximum apertures are 0 and 6.  75  Figure 5.1: Sample input for the random generator.  Running the generator for the input given in Figure 5.1, 400 elements with random apertures is generated over a square with side length of five.  Figure 5.2 shows the distribution of the random values generated as the apertures over a single fracture for this input. The curve illustrates a normal distribution with an average equal to 3.0. Since the difference between mean aperture and min aperture is 3.0 (Figure 5.2) standard deviation (σ) of this distribution is 1.0. These properties are the same as the ones requested in the input file. Also, the probability of 2, which is (mean - 1σ), is 24/σ% equal to 24%. As discussed in section4.3.1 this is another indication of the correctness of the generated normal distribution.  76  Figure 5.2: Distribution of the aperture values generated for Verification Case I.  5.3 Verification Case II: Continuity Equations in the Solver The first solver verification checks the solution of the continuity equations in the code. In this case a network including four pipes is considered (Figure 5.3). All pipes have identical properties, which are given in Table 5.1. The grouting properties are shown in Figure 5.4. The flow rate given in this figure is the injection flow rate which is applied to Node 1.  77  Figure 5.3: Pipe network for verification case II with four similar pipes.  Figure 5.4: The input file containing grout and grouting properties for verification case II.  Table 5.1: The input properties of the pipe network for verification case II.  Since the pipes have identical properties there will be an equal flow rate in all of them. Knowing the flow rate in each of the pipes, which is the total flow rate from the input divided by four, the penetration length and pressure can be calculated in each time step.  78  The penetration length and pressure are computed using Equation 4.8 and Equation 4.10. As an example for the first time step in pipe 12:  Q12 =  Equation 5.1  Δt =  Equation 5.2  l12 = l 012 +  Equation 5.3  Qinjected 4  m3 0.004 = = 0.001 s 4  5 totaltime = = 1s numberofincrements 5  Q12 0.001 × 1 = 0 .5 m .Δt = 0 + A12 2 × 0.001 × 1  Knowing Q, ΔP can be calculated from Equation 4.10, using two mathematical parameters α and φ : Equation 5.4  α = −3 −  12 × 0.01 × 0.001 12 μ .Q = − − = −33 3 4 × 1 × (0.001) 2 × 1 4 wb 2τ 0B  Equation 5.5  φ =2  −α 1 1 33 . cos( (cos −1 ( 27 . cos( (cos −1 ( 27 3 ) + π ) = 0.0606 3 ) +π) = 2 −α 33 3 3 3 3  2τ 0B ΔP f (Q) = =− ⇒ l 2b.φ Equation 5.6  2τ 0B 1 ΔP = − ⋅l = − × 0.5 = 8.25kPa 2b.φ 0.001× 0.0606  This process of calculation is performed for each time step. Since the pipes are similar the calculated values in all of them are the same. The summary of the results based on these equations is given in Table 5.3 to be compared to the results of the code shown in Table 5.2.  79  The outputs show a close match in the results from both methods. The penetration lengths and pressures are identical. Also the pressure values which are rounded to 2 significant figures in both tables are the same. The pressure column show the pressure values at the beginning of each pipe. Since in this case all of the pipes start from the same node which is Node 1, pressures given in the results are equal in all of the pipes at each time step.  Table 5.2: Code results for verification case II.  80  Table 5.3: Hand calculation results for verification case II.  5.4 Verification Case III: Head Drop Equations in the Solver and Network Update. This case verifies the head drop in the branches attached to a single node. This verification case also checks the network updating when the pipes become full or dead. Figure 5.5 shows the pipe network. Three pipes with different apertures are attached to the injection point at Node 1. Pipe 45 (the pipe connecting Node 4 to Node 5) is attached to the end of pipe 14. This pipe is added to the system when pipe 14 is full.  The grout and network properties are given in Figure 5.6 and Table 5.4. In the first time step grout flows in three pipes (pipes 12, 13, 14) two of which have the same aperture. 81  Knowing that the flow rate is equal in pipes 13 and 14 and also setting the pressure drop equal in pipes 12 and 13 the penetration length and flow rate are calculated in the first two time steps using equations that are the same as for the previous verification case. In the second time step pipe 13 becomes full. Since it does not have any extensions it is removed from the system and the grout flows in the other two pipes from time step 3. In time step 3 the penetration length in pipe 14 becomes larger than the pipe length. So the extension is added to the end of the pipe. Flow rate and penetration length is computed in each of the pipes using the equations given in Section 4.3.2 . The same process as described here is done in the code to compute the results shown in Table 5.5. To verify the penetration lengths and injection pressure hand calculated results are shown in Table 5.6.  Figure 5.5: Pipe network for verification case III.  Figure 5.6: The input file containing grout and grouting properties for verification case III. 82  Table 5.4: The input properties of the pipe network for verification case III.  Similar to the Verification Case I, the results from the code and hand calculations are identical, except for small differences in the computed pipe pressure. To minimize the size of outputs the values are rounded to 2 significant figures when writing into output files. In this case this rounding causes a subtle difference between the code results and hand calculations.  Table 5.5: Code results for verification case III.  83  Table 5.6: Hand calculation results for verification case III.  5.5 Verification Case IV: Water Case If the initial yield stress ( τ 0B ) of a Bingham fluid is set to zero the fluid will behave as a Newtonian fluid and the Bingham flow equations will be simplified into the Newtonian fluid flow equations. To verify the solution for a Newtonian fluid a small value for the yield stress is used with injection into a pipe network identical to that used in Verification Case II (see Figure 5.7). It should be noted that to avoid mathematical divide by zero errors, the yield stress is set at a value close to zero, but not zero.  The results are verified by comparing the theoretical transmissivity for a parallel plate fracture (derived from the cubic law) with the transmissivity computed based on the results of the code.  84  Based on the definition tranmissivity is: T =  Equation 5.7  q (Δh / L)  q is the flow rate per unit width. Δh / L is the hydraulic gradient which is equal to  ΔP . ρgL  The transmissivity based on the cubic law for a Newtonian fluid is calculated as: Tcubic =  Equation 5.8  ρge 3 12 μ  And based on the results from the code is: Tcode =  Equation 5.9  ρgq ΔP / L  Figure 5.7: Pipe network for verification case IV with four similar pipes.  Figure 5.8 shows the input properties for the grouting. The input properties for the network are also in Table 5.7. These data illustrate that all of the input parameters are identical to the verification case II except for the yield stress which is equal to 10-8 Pa.  85  Figure 5.8: The input file containing grout and grouting properties for verification case IV.  Table 5.7: The input properties of the pipe network for verification case IV.  The results from the code again give identical values to the hand calculations, except for the pressure in the pipes which is slightly different due to rounding errors.  86  Table 5.8: Code results for verification case IV.  87  Table 5.9: Hand calculation results for verification case IV.  To calculate transmissivity, density ( ρ ) and gravitational acceleration ( g ) are assumed as 1000kg/m3 and 9.81m/s. The computed transmissivity values in different time steps are shown in Table 5.10. this table also shows the relative error which is computed by the following equation:  Error =  Equation 5.10  TCode − THandCalc × 100 THandCalc  TCode and THandCalc in this equation are transmissivity from code results and hand calculations. Results show a maximum 0.336% error in the transmissivity values.  88  Table 5.10: Transmissivity in different time steps from the code result and hand calculations.  89  6 SAMPLE SIMULATIONS AND RESUTLS 6.1 Overview of the Chapter This chapter presents a few sample simulations by the model. These examples are intended to show some sample results from the code. Input properties chosen for the simulations are reviewed at the beginning of this chapter. Output of the model contains the variations of pressure, injected flow rate and penetration length of grout over grouting time, giving the capability to present the results in many different ways. In this chapter different methods of analyzing the results based on the results of the code are introduced. The outputs of the sample simulations are presented by different figures to compare the results from different aspects. The results are also compared against the results of a constant aperture fracture model.  6.2 Input Properties for Simulations Fracture Geometry and Mesh Properties:  As discussed in Section 4.3.1, the required input parameters for the fracture generator are dimensions of the fracture, aperture distribution over the surface of the fracture, and number of 2D elements. The properties assumed for the simulation here are:  •  Fracture and element dimensions: The dimensions of the fracture are chosen based on the computing capacities of the code. In all of the simulations fractures are 3m* 3m squares which are divided into 400 square elements with constant apertures over the 90  surface of each of the elements. In all of the simulations the location of the borehole is assumed to be at the center of the fracture.  •  Aperture size and distribution: Normal distributions with a mean of 0.5 to 2mm and a Coefficient of Variation (COV) equal to Standard deviation/mean of 0.03 to 0.3 were used to simulate the fracture surface. The range of mean aperture is chosen similar to the values used by Rasekh and Shuttle (2006), Rombough (2006) and Fransson (1999).  Grout Material Properties:  Grout properties chosen for simulations are shown in Table 6.1. To have realistic results these properties are consistent with the range of values obtained by Rombough (2006) from viscometer tests.  Table 6.1: Grout material properties used in the simulations Material Property Range of values  Injection Flow Rate (l/min)  Viscosity (Pa.s)  Yield Stress (Pa)  0.22  0.011  1 - 10  6.3 Analyzing the Results Since the output of the code gives penetration length in each element and also pressure and the injected volume in each time step the results can be compared in different ways. In this section four different methods are presented:  •  Comparing grout spread in the fracture over time. 91  •  Comparing P versus penetration length against that for constant aperture (constant aperture model results).  •  Comparing pressure changes over grouting time.  •  Comparing penetration length versus P.V against that for constant aperture (constant aperture model results).  Grout Penetration Diagram  One of the methods to present the code results is the grout propagation diagram over the surface of the fracture in each time step for each of the cases. These diagrams which show the limits of grout at each time step, illustrate the grout penetration patterns. These diagrams are most useful in cases where the effect of a blocking zone, for example an area with a smaller aperture, is studied.  Penetration Length from the Code Results  To calculate the code ‘penetration length’ a circular penetration is assumed. The volume at each time step is computed by multiplying the injection flow rate by the elapsed time. The penetration length derived from the computed volume would be:  L=  Equation 6.1  92  V πe  e in Equation 6.1 the is the effective aperture. In the literature different suggestion for this value is proposed. As it is explained in 2.5.1 Gustafson and Stille (1996) suggest the harmonic mean as the effective aperture for the random distributions used to simulate the fracture surface. In the current examples average aperture for the whole fracture surface is assumed as the effective aperture. However, in the future works the results of this model can be used to give suggestions on this issue.  Penetration Length Based on Constant Aperture Model Assuming Cylindrical Penetration  In Chapter 3 the one dimensional solution of Bingham fluid is presented. Shuttle et al. (2007) employed this solution in a VBA code to calculate penetration length of grout for a cylindrical penetration shown in Figure 6.1. In this figure Lg, rBH and e are penetration length of grout, borehole radius and aperture. Pg and C are injection pressure and cohesion of grout.  Figure 6.1: Cylindrical penetration of grout.  This model assumes fracture as an infinite plane with two parallel sides and a constant aperture. The code takes the grout properties and a constant aperture as input and calculates the penetration length of grout for different pressures during grouting. 93  Similar properties as given in Table 6.2 are used in the simulations with this model. The only difference is in the fracture properties. This model takes a single value for aperture which is the average aperture here, to simulate the fracture.  Pressure Changes Over Time  The code gives the grout pressure at the beginning of each of the pipes. To compare the changes of grout pressure over time the injection pressure which is the pressure at the borehole node is considered.  Pressure Changes for Different Penetration Lengths  The constant aperture code (Shuttle et al., 2007) gives the penetration length by varying the injection pressure. To compare these results against the current code, injection pressure versus penetration length curves are presented.  6.4 Sample Simulations 6.4.1 Overview of the Simulations This section presents the results of some simulations from the code. These results are meant to show some sample simulations of the model and also the variation of the results for different input parameters. Table 6.2 gives a summary of the parameters used in these simulations. In each of the simulation sets all of the parameters are kept constant except one to focus on the individual influence of each input. All cases were specified to inject 94  grout for 1200 s (20min), however case 4 was stopped before the end of the simulation because the whole fracture was filled with grout.  Table 6.2: Summary of the simulations  The simulations presented are samples to presents outputs of the code – a comprehensive set of analyses was not the intent. Hence although these results might give some ideas about the effect of each input parameter on the grouting results, no general conclusions can be derived based on the presented figures.  6.4.2 Grout Penetration for Different COV Values In this section three different Coefficients of Variation are considered. These three fractures are generated entering different standard deviations, resulting in three different minimum and maximum values of aperture while the average aperture is the same in all three cases. These values are shown in Table 6.1.  95  Table 6.3: Input properties for random fracture generation. COV  mean aperture (mm)  Minimum aperture (mm)  Maximum aperture (mm)  0.03 0.15 0.30  1.0 1.0 1.0  0.92 0.60 0.06  1.11 1.45 1.82  Figure 6.2 shows the distribution of aperture over the surface of fracture and also contours of grout penetration for a two minute interval in each case. These contours indicate the penetration limits of grout at a specific time. As expected changing the coefficient of variation has a considerable effect on the shape of grout spread.  The other observation in Figure 6.2 is the square shape of grout penetration contours. This is also observed in all other cases presented in this chapter. These square shaped contours are obtained based on the penetration length in each 1D element. Although the area inside each contour is close to the injected volume of grout divided by the average aperture over the area, the shape is not a circle. The reason is the use of 1D pipes which are connected to each node in form of a cross representing a 2D element. While drawing the contours the only data available for the penetration of grout is the penetration length in each of these pipes. Therefore, connecting the penetration points to each other forms a square rather than a circle.  96  (a) Aperture distribution over the surface of fracture for COV=0.03.  (b) Grout penetration contours for COV=0.03.  (c) Aperture distribution over the surface of fracture for COV=0.15.  (d) Grout penetration contours for COV=0.15.  (e) Aperture distribution over the surface of fracture for COV=0.30.  (f) Grout penetration contours for COV=0.30.  Figure 6.2: Generated aperture distribution (a), (c), (e) and grout penetration contours for 2 minute (120s) intervals (b), (d) and (e) for different Coefficients of Variation. 97  Figure 6.3 shows changes of pressure at the injection point versus penetration length. These results are also compared with the result of constant aperture model. In the fractures with lower COV the apertures in different elements are much closer to the average value. Therefore, it is expected that in such cases the results are much closer to the constant aperture. But Figure 6.3 shows that the code results indicate a higher injection pressure than the constant aperture model throughout the simulation. Comparing the stochastic aperture results, initially the higher COV is associated with a higher injection pressure, but with increasing injection volume the three results from all COV converge and above an injection radius of approximately 0.4m are not very sensitive to COV. Since decreasing COV does not give closer values to the constant aperture model results, the differences should be attributed to other reasons such as solving the equations for 1D pipes rather than a 2D case. The current model considers summation of penetration lengths times widths of filling pipes in the network as the injected grout volume, rather than considering a circular area considered in a pure radial flow model. This results in deviation of the simulations by the current model from pure radial flow simulations. However, any conclusion on the reason for this difference requires further studies on this subject.  98  Figure 6.3: Changes of injection pressure with grout penetration for different Coefficients of Variation.  Changes of injection pressure over time are illustrated in the following figure. This figure can be used to study the effect of COV changes on the injection pressure at any time during grouting.  99  Figure 6.4: Changes of injection pressure over time for different Coefficients of Variation.  Figure 6.5 displays the changes of penetration length versus P.V for three different COV values. As it can be observed from the results similar to Figure 6.3 the current results show that in all of the cases Stochastic aperture model is giving lower penetration lengths than Shuttle et al. (2007) results for a specific P.V and these curves are not very sensitive to the changes of COV. The same discussion as Figure 6.3 is valid for these observations.  100  Figure 6.5: Penetration length of grout for different Coefficients of Variation.  6.4.3 Grout Penetration for Different Fracture Apertures In this section the results of three simulations with three apertures of 0.5, 1 and 2 mm and coefficient of variation of 0.14 are presented. In order not to have the effect of changes in the aperture distribution over the fracture surface one set of elements with a mean aperture equal to 1mm is generated and each of apertures is multiplied by 0.5 and 2 in the other two cases (Figure 6.6).  Figure 6.6 shows that changes of aperture have considerable effect on penetration length of grout over time. However, as expected the results of the current simulations indicate little effect of mean aperture changes on the grout spread pattern.  101  (a) Aperture distribution over the surface of the fracture for aperture=0.5mm.  (b) Grout penetration contours for aperture=0.5mm.  (c) Aperture distribution over the surface of the fracture for aperture=1.0mm.  (d) Grout penetration contours for aperture=1.0mm.  (e) Aperture distribution over the surface of the fracture for aperture=2.0mm.  (f) Grout penetration contours for aperture=2.0mm.  Figure 6.6: Generated aperture distribution (a), (c), (e) and grout penetration contours for 2 minute (120s) intervals (b), (d) and (e) for different fracture apertures.  102  Figure 6.7 shows the penetration pressure versus penetration length based on the results of simulations and also the curves based on the results of constant aperture model. The curves show that larger apertures have closer results to the constant aperture code results in the current case.  Figure 6.7: Changes of injection pressure with grout penetration for different apertures.  As it is shown in Figure 6.8 pressures in the fracture with smaller aperture is much higher than the pressure in fractures with larger openings. Also, Figure 6.8 shows a sudden increase in the grout pressure around 1100 sec which is a boundary effect corresponding to the time grout is blocked by fracture outer boundaries. This stepwise increase in pressure as the grout reached the boundaries of the fracture is because of applying no flow boundary conditions for the outer boundaries of the fracture. Since this effect will 103  cause getting unrealistic results the dimensions of fracture has to be chosen large enough to avoid being filled with grout. The fractures in the simulations are chosen based on the capacities of the code which is larger than previous studies with stochastically varying aperture (e.g. Rasekh and Shuttle, 2006). To be able to model larger fractures increasing the code capacities is highly recommended for future work. There are a few methods to raise the capacity of the code like making the computing matrices in the code diagonal or reducing the number of equations solved in the code by some simplifying assumptions.  Figure 6.8: Changes of injection pressure over time for different fracture apertures.  Figure 6.9 shows the penetration lengths versus P.V along with the curves based on the results of constant aperture model. Similar to Figure 6.7 these curves become closer to  104  the constant aperture code results when fractures with larger average aperture are modeled.  Figure 6.9: Penetration length of grout for different fracture apertures.  6.4.4 Grout Penetration for Different Yield Stresses In this section the results of three sample simulations with three different yield stresses of grout are presented. In this case except yield stress of the grout, all of the other input parameters including fracture properties are identical in all of the cases. Figure 6.10 shows the grout spread diagram for different yield stresses. Grout yield stress has no significant effect on the penetration length and pattern of the grout over time in these  105  examples. The penetration diagram for the middle yield stress value (τ0=5 Pa) is not displayed due to similarities in the other two figures.  (b) Grout penetration contours for τ0=1Pa.  (a) Aperture distribution over the surface of the fracture.  (c) Grout penetration contours for τ0=10Pa. Figure 6.10: Generated aperture distribution (a) and grout penetration contours for 2 minute (120s) intervals for different yield stresses (b) and (c).  Figure 6.11 shows pressure versus penetration length for each of the cases from both results of current code and constant aperture code. By changing yield stress the pressure-  106  penetration length curves based on simulation results does not become considerably closer to or farther from the constant aperture results.  Figure 6.11: Changes of injection pressure with grout penetration for different yield stresses (cohesions).  Unlike the little effect of yield stress changes on the penetration length of grouting over time, cohesion has a large effect on the pressure changes during grouting in these sample simulations. Figure 6.12 shows that larger yield stresses give larger pressures at any time during grouting.  107  Figure 6.12: Changes of injection pressure over time for different yield stresses (cohesions).  Figure 6.13 shows penetration length versus P.V curves. Contrary to Figure 6.11 these curves illustrate by increasing the yield stress of the grout the constant aperture results and current simulation results give closer values for the penetration length of the grout. This might be due to the decrease in velocity as the yild stresses increases, which makes it closer to an ideal laminar flow.  108  Figure 6.13: Penetration length of grout for different yield stresses (cohesions).  109  7 SUMMARY AND CONCLUSION  Predicting the grout penetration in fractured rocks is of interest because of its significance in engineering and scientific applications. Although popular empirical and experimental methods give approximate solutions for the grout penetration, to date few models have been developed to improve these predictions. The main problems with modeling grouting in fractured rocks are the complexities of grout behavior and fracture geometry.  The model developed in this study provides an improved simulation of grout spread in a single fracture by:  •  Solving the Bingham flow equations over the grouting time to simulate fluid flow. Unlike many previous models which simulate grout only based on cohesion of the fluid, this model simulates grout by taking both viscosity and cohesion into account.  •  Simulating fractures with randomly changing aperture which is more realistic than constant aperture. In reality grout follows certain paths while penetrating a fracture. Since these routes are highly affected by the aperture changes over the surface of fracture, considering fracture stochastic properties gives more realistic results.  These methodologies help to get a more realistic simulation of each of the two phases of grouting. Previously using these two approaches at the same time caused some limitations in the maximum fracture size and grouting time the model can simulate. To 110  solve this problem, 1D equations in stochastic 1D channels are used instead of 2D ones. This simplification improves the computation capacity limits mentioned above by saving a considerable amount of computation. Although this simplification is a great changes from solving 2D equations in 2D elements, since in reality grout flows in thin routes, this assumption is not far from the real physics of the flow.  The model gives the penetration length, pressure and injected grout volume in each of the elements at each time step. This data can be used in different combinations to illustrate the physics of the grouting. A careful post-processing of the numerical data provides a useful tool in comparing the results from different aspects. The current code provides information to study grout penetration under different conditions by comparing:  •  Grout penetration pattern over time.  •  Pressure changes over time at any point in fracture.  •  Pressure changes for different penetration lengths.  •  Penetration length for different pressure x volume.  Considering the features of the present model in simulating fracture and grout and also giving a variety of outputs, this model can be used for any study on grouting in fractures such as sensitivity analysis on both grout and fracture properties.  111  8 REFERENCES  AFTES. (1991). "Recommendations on Grouting for Underground Works." Tunneling and Underground Space Technol., 6(4), 383-461. ASTM C150. “Standard Specification for Portland Cement.”, American Society for Testing and Materials, 2002. Amadei, B. (2000). “A Mathematical Model for Flow of Bingham Material in Fractures.” Proc.4th NARMS Conf.: Pacific Rocks 2000, “Rock around the Rim”, (J. Girard, M. Liebman, Ch. Breeds & T. Doe eds), Balkema, 2000. Bozovic, A. (1985). “Discussion of question 58: Foundation treatment for control of seepage.” Fifteenth Cong. Large Dams, 2(GP.RS.7), 959-1005. Brady, G. S., and Clauser, H. R. (1986). Materials handbook: an encyclopedia for managers, technical professiona ls, purchasing and production managers, technicians, supervisors, and foremen. McGraw Hill, New York, N.Y. Bruce, D. A. (2006). "Glossary of grouting terminology." Geotech News, 24(4), 50-59. Byrne P.M., Imrie A.S., and Marcuson W.F., (1988), Civil Engineering—ASCE, Vol. 58, No. 12, December 1988, pp. 50-53 Cambefort, H. (1964) “Injection des Sols.” Eyrolles, Paris, 1964. Chhabra, R. and Richardson, J. (1999). “Non-Newtonian Flow in the Process Industries: Fundamentals and Engineering Applications” Chapter 3, Butterworth-Heinemann, ISBN: 0750637706. Clarke, W. J. "Microfine Cement Technology, Dec. 6-9, 1987." 1-13. Cook, NGW. (1992). “Jaeger Memorial Dedication Lecture Natural Joints in Rock: Mechanical, Hydraulic and Seismic Behaviour and Properties under Normal Stress.” Int. J. Rock Mech. Min. Sci., Vol. 29, No.3, pp. 198-223. 112  Deere, D. U., and Lombardi, G. (1985). "Grout slurries - thick or thin?" Issues in Dam Grouting. Proceedings of the session held in conjunction with the 1985 ASCE Convention. ASCE, New York, NY, USA, Denver, CO, USA, 156-164. Eriksson, M. (2002). "Prediction of grout spread and sealing effect: a probabilistic approach." PhD Thesis, Department of Civil and Architectural Engineering, Royal Institute of Technology. Eriksson, M., and Stille, H. (2004). "A Method for Measuring and Evaluating the Penetrability of Grouts." Proceedings of the Geo-Institute and Deep Foundations Institute, Specialty Conference on Grouting, New Orleans, Louisiana . Eriksson, M., Stille, H., and Andersson, J. (2000). "Numerical calculations for prediction of grout spread with account for filtration and varying aperture." Tunneling and Underground Space Technology, 15(4), 353-364. Evans, K.F., Kohl, T., Hopkirk, R.J., Rybach, L., 1992. Modelling of energy production from Hot Dry Rocks systems. Final report to the Swiss National Energy Research Fund, Project 359, Institute for Geophysics, ETH Zürich/Polydynamics Ltd., Zürich, April 1992. Ewert, F. K. (1997). "GIN principle revisited." International Water Power & Dam Construction, 49(10), 33-36. Fairhurst, C. (1964), “Measurement of In-situ Rock Stresses with Particular Reference to Hydraulic Fracturing.” Rock Mech. and Eng. Geol.2 130–155. Fransson, Å. (1999). “Grouting predictions based on hydraulic tests of short duration: analytical, numerical and experimental approaches”, Licentiate thesis, Chalmers University of Technology, Dept. of Geology, Publ. A93, Göteborg, Sweden, 1999. Gale, J., MacLeod, R., and LeMessurier, P. (1990). “Site characterization and validation – measurement of flowrate, solute velocities and aperture variation in natural fractures as a function of normal and shear stress”, Stage 3. Stripa Project Technical Report 90-11, SKB, Stockholm. Gallavresi, F. (1992). "Grouting improvement of foundation soils." Grouting, Soil Improvement and Geosynthetics, 1, Geotechnical Special Technical Publication, 30 1-38.  113  Galambos, J., and Kotz, S. (1978). “Characterizations of probability distributions.”, Springer-Verlag New York. Gustafson, G., and Stille, H. (1996). "Prediction of Groutability from Grout Properties and Hydrogeological Data." Tunneling and Underground Space Technology, 11(3), 325332. Hakami, E. (1995). “Aperture Distribution of Rock Fractures.” Doctoral Thesis. Department of Civil and Environmental Engineering, Royal Institute of Technology, Stockholm. Hakansson, U. (1993). "Rheology of fresh cement-based grouts." PhD Thesis, Division of Soil and Rock Mechanics, Royal Institute of Technology, Stockholm. Hakami, E., and Larsson, E. (1996). "Aperture measurements and flow experiments on a single natural fracture." Int. J. Rock Mech. Min. Sci., 33(4), 395-404. Hall, C. (1976). "On the history of Portland cementafter 150 years." J.Chem.Educ., 53(4), 222-223. Hansson, P., (1995). “Filtration Stability of Cement Grout for Injection of Concrete Structures”, IABSE Symposiun, San Fransisco, pp 1199-1204. Hassler L, Hakansson U, Stille H.(1992), “Computer-simulated Flow of grouts in jointed rock.” Tunneling and Underground Space Technology, 7(4),441-446. Henn, R. W. (1996). Practical guide to grouting of underground structures. ASCE Press; Thomas Telford Publications, New York, N.Y.; London. Houlsby, A. C. (1990). Construction and design of cement grouting : a guide to grouting in rock foundations. Wiley, New York. Hubbert, M.K., and D.G. Willis (1957) “Mechanics of Hydraulic Fracturing”, Trans. Amer. Inst. Min. Eng.210, 153–168. Jansen, R. B. (1988). “Advanced Dam Engineering for Design, Construction, and Rehabilitation”, Kluwer Academic Publishers.  114  Kelley, C. T. (1995). “Iterative Methods for Linear and Nonlinear Equations.”, Society for Industrial and Applied Mathematics (SIAM). Kikuchi, , Igari, T., Mito, Y., and Utsuki, S. (1997). "In situ experimental studies on improvement of rock masses by grouting treatment." International Journal of Rock Mechanics and Mining Sciences Geomechanics Abstracts, 34(3), 594. Kutzner, C. (1996). Grouting of rock and soil. A.A. Balkema, Rotterdam ; Brookfield, VT. Littlejohn, S. (2003). "The Development of Practice in Permeation and Compensation Grouting a Historical Review (1802-2002) Part 1 Permeation Grouting." Geotech Spec Publ, 120 50-99. Littlejohn, S., and Bruce, D. (1977). "Rock anchors: state of the art." Foundation Publication Ltd, Brentwood. Lombardi, G. (1985). "The role of cohesion in cement grouting of rock." Proceedings of the 15th Congress on Large Dams, Lausanne, International Commission on Large Dams (ICOLD), 235–261. Mongilardi, E., and Tornaghi, R. (1986). "Construction of Large Underground Openings and Use of Grouts." Proc. Intl. Conf. on Deep Foundations. Beijing, China, September. Morgenstern, N. R., and P. R. Vaughan (1963) “Some Observations on Allowable Grouting Pressures” Symposium on Grouts and Drilling Muds, Butterworths, London, 1963. NRC, National Research Council, (1996). Rock Fracture and Fluid Flow: Contemporary Understanding and Applications. National Academy Press, Washington, US. Olsson, R., and Barton, N. (2001). "An improved model for hydromechanical coupling during shearing of rock joints." Int. J. Rock Mech. Min. Sci., 38(3), 317-329. Oda M., Y.Hatsuyama, Y.Onishi (1987) “Numerical experiments on permeability tensor and its application to jointed granite at Stripa Mine”, Sweden. d. Geophys. Res., 92, 9037-9048.  115  Rasekh A. and Shuttle D.A. (2006). “Assessment of grouting effectiveness from volume pressure response”, First Interim Report to Shell Exploration and Production, Houston, TX. Rombough V. (2006). “Control on penetrability of GIN mixes during fractured rock grouting” Batchelor of Applied Science Thesis, Geological Engineering Department, University of British Columbia. Rombough, V., G. Bonin, and D.A. Shuttle, “Penetrability control of GIN mixes during fractured rock grouting”. Proceeding of the 59th Canadian Geotechnical Society Conference, Vancouver, Canada, 1-4 October 2006. Shimoda, M., and Ohmori, H. (1982). "Ultra fine grouting material." Grouting in Geotechnical Engineering.ASCE, New York, 77–91. Shroff, A. V., and Shah, D. L. (1993). Grouting technology in tunnelling and dam construction. A.A. Balkema, Rotterdam ; Brookfield. Shuttle DA, and Glynn E. (2003) “Grout Curtain Effectiveness in Fractured Rock by the Discrete Feature Network Approach.”, Proceedings of the grout and grout treatment, 2003, Reston, Virginia, ASCE, February 23–28, 2003. Shuttle, D.A., V. Rombough, and G. Bonin (2007), “Impact Of Grout Rheology on GIN”. Proceedings of the mini-symposium on grouting, Geo-Denver 2007: New Peaks in Geotechnics, 18-21 February, Denver, CO. Singh, K. (2007). “Quantitative Social Research Methods.” Sage Publications Inc. Verfel, J. (1989). Rock grouting and diaphragm wall construction. Elsevier; Distribution for the USA and Canada, Elsvier Science Pub. Co., Amsterdam ; New York; New York. Weaver, Ken (Kenneth D.). (1991). Dam foundation grouting. American Society of Civil Engineers, New York, N.Y. Warner, J. (2004). “Practical Handbook of Grouting: Soil, Rock, and Structures”, Wiley. Wong, H. Y., and Farmer, I. W. (1973). "HYDROFRACTURE MECHANISMS IN ROCK DURING PRESSURE GROUTING." Rock Mechanics, 5(1), 21-41. 116  Xanthakos, P. P., Abramson, L. W., and Bruce, D. A. (1994). Ground control and improvement. J. Wiley, New York. Yesilnacar, M. I. (2003). "Grouting applications in the Sanliurfa tunnels of GAP, Turkey." Tunnelling and Underground Space Technology, 18(4), 321-330.  117  APPENDIX A : DERIVATION OF THE BINGHAM FLUID EQUATIONS Figure A.1 shows the forces in an element of flowing fluid in a pipe. Based on this figure the force equilibrium equations can be written as:  [P − ( P + ΔP)]× y × w = τ yz × L × w  Equation A.1  ⇒−  ΔP × y = τ yz L  Where w, L and y are the width, length and height of the element. P and P + ΔP are the pressures applying at each of the faces of the element. τ yz is the shear stress in the fluid at the height equal to y.  Figure A.1: Forces in a flowing fluid inside a channel  The mathematical equations for a Bingham fluid behavior is: Equation A. 2  τ yz = − μ B  du + τ 0B dy  τ yz = τ 0B (  118  du = 0) dy  for τ yz ≥ τ 0B for τ yz ≤ τ 0B  In the above equations μ B is the dynamic viscosity,  du is the rate of the shear strain and dy  τ 0B is the yield stress. Replacing τ yz by its mathematical definition for a Bingham fluid in the Equation A.1 the fluid velocity equation can be derived. Based on these equations the injection rate per unit width would be: For τ yz ≤ τ 0B we have:  −  Equation A. 3  ΔP × y = τ 0B L  ⇒ yP = −  τ 0B ΔP L  If y ≥ y P then:  ΔP du y = −μ B + τ 0B ⇒ L dy du 1 ΔP = ( y + τ 0B ) dy μ B L  −  u=  Equation A. 4  1  μB  If y = b ⇒ u = 0 therefore: u=  1  μB  c=−  (  ΔP 2 b + τ 0B b) + c = 0 ⇒ 2L  1  μB  119  (  ΔP 2 b + τ 0B b) 2L  (  ΔP 2 y + τ 0B y ) + c 2L  u=  Equation A. 5  1  μB  (  ΔP 2 1 ΔP 2 y + τ 0B y ) − ( b + τ 0B b) 2L μ B 2L  If y ≤ y P then: du = 0 ⇒ u = const. dy  Equation A. 6  If y = y P then:  τ 0B  u y≤ y p = u y≥ y p ( y P = −  u y≥ y p  ΔP L  )  2  τ 0B ⎞ ΔP ⎛ τ 0B ⎞ 1 ΔP 2 B⎛ ⎜ ⎟ ⎜ ⎟) − = ( − +τ 0 ⎜− ( b + τ 0B b) ⇒ ⎜ ⎟ ⎟ μ B 2 L ⎝ ΔP L ⎠ ⎝ ΔP L ⎠ μ B 2 L  u y≥ y p =  1  1  (  τ 0B  2  μ B 2ΔP L  u y≥ y p = −  1  (  τ 0B  −  τ 0B  ΔP L  2  μ B 2ΔP L  2  )−  Equation A. 7  1  μB  )− (  1  μB  (  ΔP 2 b + τ 0B b) ⇒ 2L  ΔP 2 b + τ 0B b) = u y ≤ y p 2L  u y≤ y p = −  1  (  τ 0B  2  μ B 2ΔP L  )−  1  μB  (  ΔP 2 b + τ 0B b) 2L  Knowing the velocity at each level over the height of the pipe the flow rate can be calculated. Flow rate will be the integration of the velocity over the height of the channel. Since half of the pipe is considered in all of the equations until now, the total q will be twice the integration of u over b.  120  b ⎡yp ⎤ q = 2 ⎢ ∫ u y ≤ y p dy + ∫ u y ≥ y p dy ⎥ ⎢⎣ 0 ⎥⎦ yp  Equation A. 8  −  τ 0B  2 ⎡ 1 ⎤ τ 0B 1 ΔP 2 B u dy b b = − − + τ ( ) ( ) ⎢ ⎥ dy ⇒ 0 ≤ y y ∫0 p ∫0 ⎢ μ B 2ΔP L μ B 2L ⎥⎦ ⎣ 2 yp ⎡ 1 ⎤ τ 0B τ 0B 1 ΔP 2 B ∫0 u y ≤ y p dy =⎢⎢− μ B ( 2ΔP L ) − μ B ( 2L b + τ 0 b)⎥⎥ × − ΔP L ⇒ ⎣ ⎦  yp  ΔP L  1 τ 0B 2 τ 0B ∫0 u y ≤ y p dy = μ B ( 2(ΔP L )2 ) + μ B ( 2 b + ΔP L b)  yp  b  ∫u  τ 0B  1  b  y≥ y p  yp  ∫  dy = −  τ 0B  1  μB  (  3  2  1 ΔP 2 ΔP 2 ( b + τ 0B b)dy ⇒ y + τ 0B y ) − 2L μ B 2L  ΔP L b  ⎡ 1 ΔP 3 τ 0B 2 ⎤ 1 ΔP 2 B = u dy ∫y y≥ y p ⎢⎣ μ B ( 6 L y + 2 y ) − μ B ( 2 L b + τ 0 b) y ⎥⎦ − p b  b  ∫u  yp  y≥ y p  τ 0B  ⇒  ΔP L  ⎤ ⎡ 1 ΔP 3 τ 1 ΔP 2 ( b + τ 0B b)b ⎥ b + b2 ) − dy = ⎢ ( 2 μ B 2L ⎦ ⎣ μ B 6L B 0  ⎡ 1 ΔP ⎛ τ B ⎞ 3 τ B ⎛ τ B ⎞ 2 ⎛ τ 0B ⎞⎤ 1 ΔP 2 B ⎜⎜ − 0 ⎟ + 0 ⎜ − 0 ⎟ ) − ⎜⎜ − ⎟⎟⎥ ⇒ −⎢ + τ ( ( ) b b 0 μ B 2L 2 ⎜⎝ ΔP L ⎟⎠ ⎢⎣ μ B 6 L ⎝ ΔP L ⎟⎠ ⎝ ΔP L ⎠⎥⎦ 3 3 b ⎞ τ 0B τ 0B 1 ⎛ ΔP 3 τ 0B 2 ⎞ 1 ⎛⎜ ∫y u y≥ y p dy = μ B ⎜⎜⎝ − 3L b − 2 b ⎟⎟⎠ − μ B ⎜ − 6(ΔP L )2 + 2(ΔP L )2 ⎟⎟ p ⎝ ⎠ 2 ⎞ 1 ⎛⎜ τ 0B 2 τ 0B − b + b⎟ ⇒ μ B ⎜⎝ 2 ΔP L ⎟⎠ 3 2 b ⎞ τ 0B τ 0B 1 ⎛⎜ ΔP 3 B 2 ∫y u y≥ y p dy = μ B ⎜ − 3L b − 3(ΔP L )2 − τ 0 b − ΔP L b ⎟⎟ p ⎝ ⎠  121  b ⎡yp ⎤ q = 2 ⎢ ∫ u y ≤ y p dy + ∫ u y ≥ y p dy ⎥ ⇒ ⎢⎣ 0 ⎥⎦ yp 3 2 ⎡ 1 ⎛ ΔP ⎞⎤ τ 0B τ 0B B 2 3 ⎜ ⎟⎥ ⎢ − − − b − b b τ 0 2 ⎜ ΔP L ⎟⎠⎥ 3(ΔP L ) ⎢ μ B ⎝ 3L q = 2⎢ ⎥⇒ B B2 ⎞ ⎢ 1 ⎛ τ 0B 3 ⎥ τ τ ⎜ + 0 b2 + 0 b⎟ ⎢+ ⎥ 2 ΔP L ⎟⎠ 2 ⎢⎣ μ B ⎜⎝ 2(ΔP L ) ⎥⎦  ⎡ 1 q = 2⎢ ⎢⎣ μ B  3 ⎛ ΔP 3 τ 0B τ 0B 2 ⎞⎟⎤ ⎜− − b + b ⎥⇒ 2 ⎜ 3L ⎟⎥ 2 ( ) Δ 6 P L ⎝ ⎠⎦  3 ⎡ 1 ⎛ ΔP τ 0B τ 0B 2 ⎞⎟⎤ 3 ⎜ q = 2⎢ b + − − b ⎥⇒ 2 2 ⎟⎠⎥ 6(ΔP L ) ⎢⎣ μ B ⎜⎝ 3L ⎦ 3 ⎞ τ 0B 1 ⎛⎜ 2ΔP 3 B 2⎟ τ − − q= b + b 0 2 ⎟ μ B ⎜⎝ 3L 3(ΔP L ) ⎠  Equation A. 9  q is the grout injection rate per unit width. To get the grout injection rate, q can be multiplied by w. Assuming φ = −  τ 0B ΔP L × b  the above equation can be rewritten as:  3 3τ 0B ⎞⎟ τ 0B 1 2ΔP 3 ⎛⎜ q=− b 1− + ⇒ μ B 3L ⎜⎝ 2(ΔP L )3 b 3 2(ΔP L )b ⎟⎠  q=− q=  1 2ΔP 3 ⎛ 1 3 3 ⎞ b ⎜1 + φ − φ ⎟ 2 ⎠ μ B 3L ⎝ 2  1 2 τ 0B 3 ⎛ 1 3 3 ⎞ b ⎜1 + φ − φ ⎟ μ B 3 φ .b ⎝ 2 2 ⎠  1 3 3 q.μ B φ = 1+ φ 3 − φ B 2 2 τ 0 .b 2 2 ⎛ 3 3 q.μ B 1 1 + φ 3 + ⎜⎜ − − B 2 2 ⎝ 2 2 τ 0 .b  ⎞ ⎟φ = 0 ⎟ ⎠  122  ⎛  φ 3 + ⎜⎜ − 3 − 3  Equation A. 10  ⎝  q.μ B ⎞ ⎟φ + 2 = 0 τ 0B .b 2 ⎟⎠  Equation A.1 is a cubic polynomial which has three roots, one of which is accepted as our solution.  φ 3 + α .φ + β = 0  α = −3 − 3  q.μ B τ 0B .b 2  β =2  The first root of the above standard equation is:  θ +π ⎛ −α ⎞ ) ⎟ . cos( 3 ⎝ 3 ⎠  φ=2 ⎜  Equation A. 11 Where:  φ =−  τ 0B ΔP L × b  ⎛ 27 ⎞ cos(θ ) = ⎜ 3 ⎟ ⎝ −α ⎠ ⎛ ⎛ 27 ⎞ ⎞ ⎟ 3 ⎟⎟ α − ⎠⎠ ⎝ ⎝  θ = cos −1 ⎜ ⎜ ⎜  ⎛ q.μ ⎞ α = ⎜⎜ − 3 − 3 B B2 ⎟⎟ τ .b ⎝  0  ⎠  123  APPENDIX B :  NEWTON’S  METHOD  FOR  THE  SOLUTION OF PIPE NETWORK EQUATIONS  Newton’s method is one of the simplest methods in solving a set of nonlinear equations. In the mathematical language Newton’s method is:  [x+ ] = [xc ] − [ f ' (xc )]−1 .[ f (xc )]  Equation B. 1  In the above equation: x+ : The root of the equations.  xc : Initial assumption for the root. f ' ( x c ) : Jacobean matrix of the given equations for the initial assumption. f ( x c ) : The value of the equations given in Section 4.3.2 for the initial assumption.  The brackets in Equation B.1 indicate that all of these variables are matrices. The detailed theory, advantages and disadvantages of the method can be found in “Iterative Methods for Linear and Nonlinear Equations.” by C.T. Kelley (1995).  For the current case where we have 2n unknowns for a general network with n pipes the previous equation can be written as:  124  Equation B. 2  ⎡ ∂f1 ⎢ ⎡ x1+ ⎤ ⎡ x1c ⎤ ⎢ ∂x1c ⎢ M ⎥ ⎢ M ⎥ ⎢ M ⎢ ⎥ ⎢ ⎥ ⎢ ∂f ⎢ x n + ⎥ = ⎢ x nc ⎥ − ⎢ n ⎢ ⎥ ⎢ ⎥ ⎢ ∂x1c ⎢ M ⎥ ⎢ M ⎥ ⎢ M ⎢⎣ x 2 n + ⎥⎦ ⎢⎣ x 2 nc ⎥⎦ ⎢ ∂f 2 n ⎢ ∂x ⎣ 1c  L O  ∂f1 ∂x nc  −1  L  ∂f n ∂x nc O L  ∂f1 ⎤ ∂x 2 nc ⎥ ⎡ f1 ( xc ) ⎤ ⎥ ⎥ ⎢ M ⎥ ⎥ ⎥ ⎢ M ⎥ .⎢ f n ( xc ) ⎥ ⎥ ⎢ M ⎥ ⎥ ⎥ ⎢ ⎢ ∂f 2 n ⎥ ⎣ f 2 n ( xc )⎥⎦ ∂x 2 nc ⎥⎦  All of the equations and matrices in this section are for a single grid. In the current code such matrices for all of the grids are assembled to form a matrix for the whole system of pipes. For a grid that “n+1” pipes are attached to it (1 inflow and n outflow), the governing equations are in the following order: 1.  f1 = ∑ Q = 0  2. For a half full pipe: Q f 2...n +1 = lt − l( t − Δt ) − .Δt = 0 (n=number of pipes with outflow) A For the filled pipe with the equation number m ( 2 ≤ m ≤ n + 1 ): f m = lt − l pipem = 0  3. For a half full pipe: f n + 2...2 n = f (Qij ).ltij − f (Qik ).ltik = 0  For the filled pipe with the equation number m+n ( n + 2 ≤ m + n ≤ 2n )  125  f m+n = f (Qij ).l pipem + f (Q jl ).ltjl − f (Qik ).ltik = 0 (“jl” is the extension attached to the full pipe)  The variables also are assumed in the following order: x1...n = Q jk xn +1...2 n = ltjk  For the filled pipe: xm−1 = Qij xn+m = ltij  For the extension: xe = Q jl xn+e = ltjl  The Jacobean matrix for the general case given above will be:  Equation B. 3  ⎡ ∂f1 ⎢ ∂x ⎢ 1 ⎢ M ⎢ ∂f n ⎢ ∂x ⎢ M1 ⎢ ⎢ ∂f 2 n ⎢⎣ ∂x1  126  L O  ∂f1 ∂xn  L  ∂f n ∂xn O L  ∂f1 ⎤ ∂x2 n ⎥ ⎥ ⎥ M ⎥= ⎥ ⎥ ⎥ ∂f 2 n ⎥ ∂x2 n ⎥⎦  127  The differentiation equation (  ∂φ ∂Q  ) in the above equation is derived as follows:  f (Q) =  ΔP 2τ B =− 0 l 2b.φ  φ=2 (  27 −α 1 + π )) ) . cos( (cos −1 ( 3 3 −α 3  α = −3 −  12μ .Q 4b 2τ 0B  2τ 0B  ∂φ ∂Q  ∂f (Q ) = 2b.φ 2 ∂Q  ⎡ ⎛ ⎢ ⎜ cos −1 ( 27 ) + π ∂φ ⎢ − 1 ⎜ −α3 cos⎜ =⎢ 3 3α ∂Q ⎢ ⎜⎜ ⎢ ⎝ ⎣  ⎞ ⎛ ⎟ ⎜ cos −1 ( 27 ) + π −α ⎟ ⎜ −α3 ⎟ + 2 3 . − sin ⎜ 3 ⎟⎟ ⎜⎜ ⎠ ⎝  ∂α 12μ =− 2 B 4b τ 0 ∂Q  128  ⎞ ⎟ 3 27 −1 ⎟1 .− ⎟. 3 . 2 −α5 27 2 ⎟⎟ 1− ( ) ⎠ −α3  ⎤ ⎥ ⎥ ∂α ⎥ ∂Q ⎥ ⎥ ⎦  APPENDIX C : RANDOM FRACTURE GENERATOR Option Explicit Option Base 1 Dim ApertDrawing(20000, 10), LineDrawing(200000, 10), PipeProp(200000, 10), Elemprop(10000, 10)  Sub RandomNetworkGenerator() Dim OutWorksheet, OutWorksheet1 As Range 'Input Variables Dim Length, Width, minaperture, Maxaperture, Meanaperture As Single Dim Nlength, Nwidth, i, j As Integer Dim Distribution As String 'Calculation variables Dim Sigma, x, y, xl(1000, 1000), yl(1000, 1000), Dw, Dl As Single Dim R, Apert(1000, 1000) As Double Dim L, Output(10000, 10), K, M, N As Integer '*****READING FRACTURE DATA (DIMENSIONS OF THE FRACTURE AND THE REQUIRED MESH) , FILE="FRACTUREDATA" Length = Sheet1.Range("c2") Width = Sheet1.Range("c3") Sigma = Sheet1.Range("c5") Meanaperture = Sheet1.Range("c6") Nlength = Sheet1.Range("c8") Nwidth = Sheet1.Range("c9") Distribution = Sheet1.Range("c11") 'Mesh Dimentions Dl = Length / Nlength Dw = Width / Nwidth 'Calculating the fracture properties For i = 1 To Nlength xl(1, i) = 0# Next i For i = 1 To Nwidth yl(i, 1) = 0# Next i x = Dw / 2 K=1 M=1 N=1 For i = 1 To Nwidth y = 0# For j = 1 To Nlength If Distribution = "Uniform" Then R = Rnd Apert(i, j) = Meanaperture * R / 2 End If If Distribution = "Normal" Then R = Application.WorksheetFunction.NormInv(Rnd, Meanaperture, Sigma) Apert(i, j) = R End If If i = 1 Then  129  xl(i + 1, j) = xl(i, j) + Dw / 2 Else xl(i + 1, j) = xl(i, j) + Dw End If If j = 1 Then yl(i, j + 1) = yl(i, j) + Dl / 2 Else yl(i, j + 1) = yl(i, j) + Dl End If If j > 1 Then PipeProp(N, 1) = xl(i, j) PipeProp(N, 2) = yl(i, j) PipeProp(N, 3) = xl(i + 1, j) PipeProp(N, 4) = yl(i, j) PipeProp(N, 5) = xl(i + 1, j) - xl(i, j) If i = 1 Then PipeProp(N, 6) = Apert(i, j - 1) Else PipeProp(N, 6) = Apert(i - 1, j - 1) End If PipeProp(N, 7) = Dl / 2 PipeProp(N, 8) = "0" LineDrawing(M, 1) = "Pline" M=M+1 LineDrawing(M, 1) = PipeProp(N, 1) LineDrawing(M, 2) = "," LineDrawing(M, 3) = PipeProp(N, 2) M=M+1 LineDrawing(M, 1) = PipeProp(N, 3) LineDrawing(M, 2) = "," LineDrawing(M, 3) = PipeProp(N, 4) M=M+1 LineDrawing(M, 1) = "" M=M+1 N=N+1 End If If i > 1 Then PipeProp(N, 1) = xl(i, j) PipeProp(N, 2) = yl(i, j) PipeProp(N, 3) = xl(i, j) PipeProp(N, 4) = yl(i, j + 1) PipeProp(N, 5) = yl(i, j + 1) - yl(i, j) If j = 1 Then PipeProp(N, 6) = Apert(i - 1, j) Else PipeProp(N, 6) = Apert(i - 1, j - 1) End If PipeProp(N, 7) = Dw / 2 PipeProp(N, 8) = "1" LineDrawing(M, 1) = "Pline" M=M+1 LineDrawing(M, 1) = PipeProp(N, 1) LineDrawing(M, 2) = "," LineDrawing(M, 3) = PipeProp(N, 2) M=M+1 LineDrawing(M, 1) = PipeProp(N, 3) LineDrawing(M, 2) = "," LineDrawing(M, 3) = PipeProp(N, 4) M=M+1 LineDrawing(M, 1) = "" M=M+1  130  N=N+1 End If Elemprop(j + (i - 1) * Nlength, 1) = x Elemprop(j + (i - 1) * Nlength, 2) = y Elemprop(j + (i - 1) * Nlength, 3) = x Elemprop(j + (i - 1) * Nlength, 4) = y + Dl Elemprop(j + (i - 1) * Nlength, 5) = Apert(i, j) L = 250 * R / (0.002) ApertDrawing(K, 1) = "Cecolor" K=K+1 ApertDrawing(K, 1) = "RGB:" ApertDrawing(K, 2) = Round(L, 0) ApertDrawing(K, 3) = "," ApertDrawing(K, 4) = Round(L, 0) ApertDrawing(K, 5) = "," ApertDrawing(K, 6) = Round(L, 0) K=K+1 ApertDrawing(K, 1) = "trace" K=K+1 ApertDrawing(K, 1) = Dw K=K+1 ApertDrawing(K, 1) = x ApertDrawing(K, 2) = "," ApertDrawing(K, 3) = y y = y + Dl K=K+1 ApertDrawing(K, 1) = x ApertDrawing(K, 2) = "," ApertDrawing(K, 3) = y K=K+1 ApertDrawing(K, 1) = "" K=K+1 Next j x = x + Dw Next i 'Output Results for the Aperture Drawing Sheet2.Select Set OutWorksheet = Sheet2.Range(Cells(4, 2), Cells(3 + K, 7)) ' and copy the data across OutWorksheet.Value = ApertDrawing 'Output Results for Linear Elements Drawing Sheet3.Select Set OutWorksheet1 = Sheet3.Range(Cells(4, 2), Cells(3 + M, 4)) ' and copy the data across OutWorksheet1.Value = LineDrawing 'Output Results for Linear Elements Properties Sheet4.Select Sheet4.Range("B1") = "NTOTAL" Sheet4.Range("B2") = N - 1 Set OutWorksheet = Sheet4.Range(Cells(4, 2), Cells(3 + N, 9))  131  ' and copy the data across OutWorksheet.Value = PipeProp 'Output Results for 2D Elements Properties Sheet5.Select Set OutWorksheet1 = Sheet5.Range(Cells(4, 2), Cells(3 + Nwidth * Nlength, 6)) ' and copy the data across OutWorksheet1.Value = Elemprop Sheet1.Select End Sub  132  APPENDIX D : PIPE NETWORK SOLVER CODE  PROGRAM PIPENETWORK C C*****THIS PROGRAM CALCULATES THE PENETRATION OF GROUT IN A PIPE NETWORK C WITH C DIFFERENT PIPE PROPERTIES C C C REAL PIPELENGTH(740,740),PIPEAPERTURE(740,740) 1 ,PIPEWIDTH(740,740),Q(740,740,300),ELEMDATA(801,10) REAL LENGTHT(740,740,300) INTEGER PIPENUM,INCREMENTNUM,TIMESTEP, ROW(800), NOUT(800,6) 1 ,NIN(740),NTOTAL DOUBLE PRECISION PI,PHI(740),PHI0(740) 2 ,X0(670),PHI1,PHI2 & ,ALPHA(740),ALPHA1,ALPHA2,Q0 1 ,SIGMA, DT,THAW0, Q00,TIME,MU0,DP0, DP1, PRESSURE(740) 1 ,APERTURE(740),XNODE(740),YNODE(740),XBH,YBH INTEGER ROWNUM, J, K,ACTIVEPIPES C C  EQUATION SOLVER VARIABLES DOUBLEPRECISION FVEC(670),FJAC(670,670),X(670),S(670) DOUBLEPRECISION R,NORM,TOLR,TOL0,RR,EPS,DELTA,DPHIDQ,DALPHADQ INTEGER N CHARACTER FLAG(740,740,300)*4  C C NIN(I): INDEX OF THE PIPE WITH INFLOW ATTACHED TO A NODE C I:NODE NUMBER C NOUT(I,J): INDEX OF THE PIPE WITH OUTFLOW ATTACHED TO A NODE C I:NODE NUMBER, N(I,6):NUMBER OF ATTACHED PIPES C Q(I,J,T): THE FLOW RATE BETWEEN TWO NODES OF I AND J AT THE C TIME INCREMENT OF T C LENGTHT(I,J,T):THE PENETRATION LENGTH BETWEEN TWO NODES OF I AND J AT THE C TIME INCREMENT OF T C C*****READING GROUTING AND FLOW PROPERTIES, FILE="GROUTINGDATA" C OPEN (10,FILE='GROUTINGDATA',STATUS='OLD',ACCESS='SEQUENTIAL', 1FORM='FORMATTED') READ (10,*) XBH,YBH READ (10,*) THAW0 READ (10,*) MU0 READ (10,*) Q00 READ (10,*) TIME, INCREMENTNUM WRITE(*,33)XBH,YBH,"THAW0:",THAW0,"Mu0:",MU0,"Q0:",Q00 1 ,"GROUTINGTIME:",TIME,"INCREMENTNUM:",INCREMENTNUM 33 FORMAT("BOREHOLE LOCATION:",F5.2,",",F5.2,A9,F5.2,A5,G9.3,A5,G9.3 1 ,A16,F8.2,A16,I3)  133  WRITE(*,*) CLOSE (10,STATUS='KEEP') C C*****READING NETWORK DATA, FILE="NETWORKDATA" C EPS=-0.00000 DELTA=10000000  OPEN (5,FILE='NETWORK',STATUS='OLD',ACCESS='SEQUENTIAL', 1FORM='FORMATTED') READ (5,*) READ (5,*)NTOTAL READ (5,*) C C C C C C C C  ELEMDATA(N,1-4): X,Y OF THE STARTING AND ENDING POINT OF THE ELEMENT ELEMDATA(N,5): LENGTH OF THE ELEMENT ELEMDATA(N,6): APERTURE OF THE ELEMENT ELEMDATA(N,7): WIDTH OF THE ELEMENT ELEMDATA(N,8): DIRECTION OF THE ELEMENT (0:HORIZONTAL, 1:VERTICAL) ELEMDATA(N,9-10):START AND END NODE NUMBER OF THE ELEMENT  1 1 10  DO 10 N=1, NTOTAL READ(5,*) ELEMDATA(N,1),ELEMDATA(N,2) ,ELEMDATA(N,3),ELEMDATA(N,4),ELEMDATA(N,5) ,ELEMDATA(N,6),ELEMDATA(N,7),ELEMDATA(N,8) CONTINUE CLOSE (5,STATUS='KEEP') ROW(1)=0  1 1  CALL NETWORKDATA(1,1,XBH,YBH,ELEMDATA,NTOTAL,PIPENUM,PIPELENGTH ,PIPEAPERTURE,PIPEWIDTH,NIN,NOUT,ROW,XNODE,YNODE ,ACTIVEPIPES)  NODENUM=PIPENUM+1 OPEN (20,FILE='RESULTS',STATUS='NEW',ACCESS='SEQUENTIAL', 1FORM='FORMATTED') OPEN (30,FILE='Jacobian',STATUS='NEW',ACCESS='SEQUENTIAL', 1FORM='FORMATTED') WRITE(20,106)"STEP","I","J","XI","YI","XJ","YJ","Q","PENET","PRES" WRITE(*,106) "STEP","I","J","XI","YI","XJ","YJ","Q","PENET","PRES" 106 FORMAT(7A5,3A7) C C*****FINDING ALL THE NODES CONNECTED TO A SPECIFIC NODE C DEFINING THE PIPES WITH INFLOW AND OUTFLOW C PI=4.0*ATAN(1.0) DT=TIME/INCREMENTNUM DO 45 I=1,2*ACTIVEPIPES X(I)=0.0 X0(I)=0.0  134  45  CONTINUE  C C*****CALCULATING PENETRATION LENGTH AND Q FOR EACH PIPE C C C COMPUTING THE ROOTS OF THE CONTROLING EQUATIONS IN EACH TIMESTEP C BY THE USE OF NEWTON ITERATIVE METHOD C C  DEFINING NUMBER OF UNKNOWNS AND PRECISION FACTORS DO 50 TIMESTEP=1,INCREMENTNUM TOLR=0.001*Q00 TOL0=0.001*Q00  620  DO 62 L=1,2*ACTIVEPIPES FVEC(L)=0.0 DO 620 LL=1,2*ACTIVEPIPES FJAC(L,LL)=0.0 CONTINUE X0(L)=X(L)  62 CONTINUE C C ASSEMBLING THE NETWORK MATRIX C ROWNUM=0 DO 60 I=1,NODENUM C C NUMBER OF UNKNOWNS AT EACH GRID C ROWNUM=ROW(I) N=2*NOUT(I,6) IF (N .EQ. 0) GOTO 60 C C SETTING THE INITIAL VALUES FOR X AND OTHER VECTORS DO 65 L=1,NOUT(I,6) IF (TIMESTEP .EQ. 1) THEN X(ROWNUM+L)=0.0 X(ROWNUM+L+NOUT(I,6))=0.0 GOTO 67 ENDIF  1  X(ROWNUM+L)=Q(I, NOUT(I,L), TIMESTEP-1) IF (I .GT. 1 .AND. FLAG(NIN(I), I, TIMESTEP-1) .NE. 'FULL') X(ROWNUM+L+NOUT(I,6))=0.0 IF (FLAG(I, NOUT(I,L), TIMESTEP-1) .EQ. 'FULL') THEN X(ROWNUM+L+NOUT(I,6))=PIPELENGTH(I,NOUT(I,L)) SIGMA=0.0  C CALCULATING THE INITIAL ASSUMPTION FOR THE FLOW RATE IN THE PIPES ACCORDING TO THEIR APPERTURE C NOT TO HAVE ZERO VAKLUES IN THE SYSTEM OF EQUATIONS AT THE BEGINNING OF THE CALCULATIONS DO 6501 INEW=1, NOUT(NOUT(I,L),6) SIGMA=SIGMA+PIPEAPERTURE(NOUT(I,L),NOUT(NOUT(I,L),INEW))  135  6501  1 1 1 1 1 1 6502  67 65  CONTINUE DO 6502 INEW=1, NOUT(NOUT(I,L),6) Q(NOUT(I,L),NOUT(NOUT(I,L),INEW),TIMESTEP-1)= Q(I,NOUT(I,L),TIMESTEP-1) *PIPEAPERTURE(NOUT(I,L),NOUT(NOUT(I,L),INEW)) /SIGMA X(ROW(NOUT(I,L))+INEW) =Q(NOUT(I,L),NOUT(NOUT(I,L),INEW),TIMESTEP-1) X(ROW(NOUT(I,L))+INEW+NOUT(NOUT(I,L),6)) =X(ROW(NOUT(I,L))+INEW)*DT /PIPEAPERTURE(NOUT(I,L),NOUT(NOUT(I,L),INEW)) CONTINUE ENDIF APERTURE(ROWNUM+L)=PIPEAPERTURE(I, NOUT(I,L)) CONTINUE  C C C  DEFINING THE EQUATIONS Q0=Q00 IF (I.GT.1 .AND. TIMESTEP.NE.1) Q0=Q(NIN(I), I, TIMESTEP-1) IF (I.GT.1 .AND. TIMESTEP.EQ.1) Q0=0 IF (I.GT.1 .AND. TIMESTEP.NE.1 .AND. FLAG(NIN(I),I,TIMESTEP-1) 1 .NE. 'FULL') Q0=0 FVEC(ROWNUM+1)=-Q0  1  1 1  1 1  97 1 1  1 1  DO 100 K=1,N/2 FVEC(ROWNUM+1)=FVEC(ROWNUM+1)+X(ROWNUM+K) FVEC(ROWNUM+K+1)=X(ROWNUM+K+N/2)-X0(ROWNUM+K+N/2) -X(ROWNUM+K)/APERTURE(ROWNUM+K)/PIPEWIDTH(I,NOUT(I,K))*DT IF(TIMESTEP .GT. 1) THEN IF(FLAG(I,NOUT(I,K),TIMESTEP-1).EQ.'FULL') THEN FVEC(ROWNUM+K+1)=0.0 X0(ROWNUM+K+N/2)=X(ROWNUM+K+N/2) ENDIF ENDIF ALPHA(ROWNUM+K)=-3-12*MU0*X(ROWNUM+K)/4/PIPEWIDTH(I,NOUT(I,K)) /APERTURE(ROWNUM+K)**2/THAW0 PHI(ROWNUM+K)=2*(-ALPHA(ROWNUM+K)/3)**.5 *COS((ACOS((-27/ALPHA(ROWNUM+K)**3)**.5)+PI)/3) IF (K .EQ. 1) GOTO 100 DP0=2*THAW0/2/APERTURE(ROWNUM+K-1)/PHI(ROWNUM+K-1) *X(ROWNUM+K-1+N/2) DP1=2*THAW0/2/APERTURE(ROWNUM+K)/PHI(ROWNUM+K) *X(ROWNUM+K+N/2) IF(TIMESTEP .EQ. 1) GOTO 99 IF(FLAG(I,NOUT(I,K-1),TIMESTEP-1).EQ.'FULL')THEN IEXTENSION=NOUT(I,K-1) ALPHA1=-3-12*MU0*X(ROW(IEXTENSION)+1)/4 /PIPEWIDTH(IEXTENSION,NOUT(IEXTENSION,1)) /PIPEAPERTURE(IEXTENSION,NOUT(IEXTENSION,1))**2/THAW0 PHI1=2*(-ALPHA1/3)**.5*COS((ACOS((-27/ALPHA1**3)**.5)+PI)/3) DP0=DP0 +2*THAW0/2/PIPEAPERTURE(IEXTENSION,NOUT(IEXTENSION,1)) /PHI1*X(ROW(IEXTENSION)+1+NOUT(IEXTENSION,6)) IF(FLAG(IEXTENSION,NOUT(IEXTENSION,1),TIMESTEP-1).EQ.'FULL')THEN IEXTENSION=NOUT(IEXTENSION,1) GOTO 97 ENDIF ENDIF IF(FLAG(I,NOUT(I,K),TIMESTEP-1).EQ.'FULL')THEN  136  98 1 1  1 1  99 100 60 C C C C  IEXTENSION=NOUT(I,K) ALPHA2=-3-12*MU0*X(ROW(NOUT(I,K))+1)/4 /PIPEWIDTH(NOUT(I,K),NOUT(NOUT(I,K),1)) /PIPEAPERTURE(NOUT(I,K),NOUT(NOUT(I,K),1))**2/THAW0 PHI2=2*(-ALPHA2/3)**.5*COS((ACOS((-27/ALPHA2**3)**.5)+PI)/3) DP1=DP1 +2*THAW0/2/PIPEAPERTURE(NOUT(I,K),NOUT(NOUT(I,K),1)) /PHI2*X(ROW(NOUT(I,K))+1+NOUT(NOUT(I,K),6)) IF(FLAG(IEXTENSION,NOUT(IEXTENSION,1),TIMESTEP-1).EQ.'FULL')THEN IEXTENSION=NOUT(IEXTENSION,1) GOTO 98 ENDIF ENDIF FVEC(ROWNUM+K+N/2)=DP1-DP0 CONTINUE CONTINUE DEFINING THE DIFFERENTIATION OF UNKNOWNS IN THE JACOBIAN MATRIX  ITERATION=0 lgt=0 111 ROWNUM=0 DO 6000 I=1,NODENUM ROWNUM=ROW(I) N=NOUT(I,6)*2 DO 150 K=1,N/2  1 1 1  1  ALPHA(ROWNUM+K)=-3-12*MU0*X(ROWNUM+K)/4/PIPEWIDTH(I,NOUT(I,K)) /APERTURE(ROWNUM+K)**2/THAW0 ALPHA0=-3-12*MU0*X0(ROWNUM+K)/4/PIPEWIDTH(I,NOUT(I,K)) /APERTURE(ROWNUM+K)**2/THAW0 PHI(ROWNUM+K)=2*(-ALPHA(ROWNUM+K)/3)**.5 *COS((ACOS((-27/ALPHA(ROWNUM+K)**3)**.5)+PI)/3) PHI0(ROWNUM+K)=2*(-ALPHA0/3)**.5*COS((ACOS((-27/ALPHA0**3)**.5)+PI)/3) FJAC(ROWNUM+1,ROWNUM+K)=1 FJAC(ROWNUM+K+1,ROWNUM+K)=-DT/APERTURE(ROWNUM+K) /PIPEWIDTH(I,NOUT(I,K)) IF(TIMESTEP .GT. 1) THEN IF(FLAG(I,NOUT(I,K),TIMESTEP-1).EQ.'FULL')THEN FJAC(ROWNUM+K+1,ROWNUM+K)=EPS FJAC(ROW(NOUT(I,K))+1,ROWNUM+K)=-1 ENDIF ENDIF FJAC(ROWNUM+K+1,ROWNUM+K+N/2)=1 IF (K.EQ.N/2) GOTO 1500 IF (X(ROWNUM+K) .EQ. 0) THEN FJAC(ROWNUM+K+N/2+1,ROWNUM+K)=EPS ELSE  1 1 1 1 1 1  DALPHADQ=-12*MU0/4/PIPEWIDTH(I,NOUT(I,K)) /(APERTURE(ROWNUM+K))**2/THAW0 DPHIDQ=-((-3*ALPHA(ROWNUM+K))**(-.5) *COS((ACOS((-27/ALPHA(ROWNUM+K)**3)**.5)+PI)/3) -(-ALPHA(ROWNUM+K)/3)**.5 *SIN((ACOS((-27/ALPHA(ROWNUM+K)**3)**.5)+PI)/3) *(-27/ALPHA(ROWNUM+K)**5)**.5 /(1-ABS(-27/ALPHA(ROWNUM+K)**3))**.5)*DALPHADQ  137  FJAC(ROWNUM+K+N/2+1,ROWNUM+K)=2*THAW0*DPHIDQ 1 /2/APERTURE(ROWNUM+K)/PHI(ROWNUM+K)**2*X(ROWNUM+K+N/2) ENDIF  IF (X(ROWNUM+K+1) .EQ. 0) THEN FJAC(ROWNUM+K+N/2+1,ROWNUM+K+1)=EPS ELSE  1 1 1 1 1 1  DALPHADQ=-12*MU0/4/PIPEWIDTH(I,NOUT(I,K+1)) /(APERTURE(ROWNUM+K+1))**2/THAW0 DPHIDQ=-((-3*ALPHA(ROWNUM+K+1))**(-.5) *COS((ACOS((-27/ALPHA(ROWNUM+K+1)**3)**.5)+PI)/3) -(-ALPHA(ROWNUM+K+1)/3)**.5 *SIN((ACOS((-27/ALPHA(ROWNUM+K+1)**3)**.5)+PI)/3) *(-27/ALPHA(ROWNUM+K+1)**5)**.5 /(1-ABS(-27/ALPHA(ROWNUM+K+1)**3))**.5)*DALPHADQ  FJAC(ROWNUM+K+N/2+1,ROWNUM+K+1)=-2*THAW0*DPHIDQ 1 /2/APERTURE(ROWNUM+K+1)/PHI(ROWNUM+K+1)**2*X(ROWNUM+K+N/2+1) ENDIF FJAC(ROWNUM+K+N/2+1,ROWNUM+K+N/2)=-2*THAW0/2 1 /APERTURE(ROWNUM+K)/PHI(ROWNUM+K) FJAC(ROWNUM+K+N/2+1,ROWNUM+K+N/2+1)=2*THAW0/2 1 /APERTURE(ROWNUM+K+1)/PHI(ROWNUM+K+1)  1500  IF (TIMESTEP .EQ. 1) GOTO  159  IF (NOUT(I,K) .EQ. 0) GOTO 25 IF(FLAG(I,NOUT(I,K),TIMESTEP-1).EQ.'FULL' .AND. K .LT. N/2) THEN IEXTENSION=I ICOLUMN=K+N/2 NEXTENSION=K FJAC(ROW(I)+K+N/2+1,ROW(IEXTENSION)+ICOLUMN)=EPS  1  1 1 1 1 1 1 1 1 1  1 1 1  150  IF(X(ROW(NOUT(IEXTENSION,NEXTENSION))+1) .EQ. 0.0 .OR. PHI(ROW(NOUT(IEXTENSION,NEXTENSION))+1).EQ.0.0) THEN FJAC(ROW(I)+K+N/2+1,ROW(NOUT(IEXTENSION,NEXTENSION))+1)=EPS ELSE DALPHADQ=-12*MU0/4/PIPEWIDTH(NOUT(IEXTENSION,NEXTENSION) ,NOUT(NOUT(IEXTENSION,NEXTENSION),1))/(APERTURE( ROW(NOUT(IEXTENSION,NEXTENSION))+1))**2/THAW0 DPHIDQ=-((-3*ALPHA(ROW(NOUT(IEXTENSION,NEXTENSION))+1))**(-.5) *COS((ACOS((-27/ALPHA(ROW(NOUT(IEXTENSION,NEXTENSION)) +1)**3)**.5)+PI)/3)-(-ALPHA(ROW(NOUT(IEXTENSION ,NEXTENSION))+1)/3)**.5*SIN((ACOS((-27/ALPHA(ROW( NOUT(IEXTENSION,NEXTENSION))+1)**3)**.5)+PI)/3) *(-27/ALPHA(ROW(NOUT(IEXTENSION,NEXTENSION))+1)**5)**.5 /(1-ABS(-27/ALPHA(ROW(NOUT(IEXTENSION,NEXTENSION)) +1)**3))**.5)*DALPHADQ FJAC(ROW(I)+K+N/2+1,ROW(NOUT(IEXTENSION,NEXTENSION))+1) =2*THAW0*DPHIDQ/2/APERTURE(ROW(NOUT (IEXTENSION,NEXTENSION))+1)/PHI(ROW(NOUT(IEXTENSION ,NEXTENSION))+1)**2*X(ROW(NOUT(IEXTENSION,NEXTENSION))  138  1  +1+NOUT(NOUT(IEXTENSION,NEXTENSION),6)) ENDIF  1  1 1 1  1 1  IF(PHI(ROW(NOUT(IEXTENSION,NEXTENSION))+1) .EQ. 0.0) THEN FJAC(ROW(I)+K+N/2+1,ROW(NOUT(IEXTENSION,NEXTENSION)) +1+NOUT(NOUT(IEXTENSION,NEXTENSION),6))=EPS ELSE FJAC(ROW(I)+K+N/2+1,ROW(NOUT(IEXTENSION,NEXTENSION)) +1+NOUT(NOUT(IEXTENSION,NEXTENSION),6))= -2*THAW0/2/APERTURE(ROW(NOUT(IEXTENSION,NEXTENSION))+1) /PHI(ROW(NOUT(IEXTENSION,NEXTENSION))+1) ENDIF IF(FLAG(NOUT(IEXTENSION,NEXTENSION),NOUT(NOUT(IEXTENSION ,NEXTENSION),1),TIMESTEP-1).EQ.'FULL')THEN FJAC(ROW(I)+K+N/2+1,ROW(NOUT(IEXTENSION,NEXTENSION))+1 +NOUT(NOUT(IEXTENSION,NEXTENSION),6))=EPS ICOLUMN=1+NOUT(NOUT(IEXTENSION,NEXTENSION),6) IEXTENSION=NOUT(IEXTENSION,NEXTENSION) NEXTENSION=1 GOTO 159 ENDIF ENDIF  25  IF (NOUT(I,K+1) .EQ. 0) GOTO 150 IF(FLAG(I,NOUT(I,K+1),TIMESTEP-1).EQ.'FULL' .AND. K .LT. N/2) THEN  IEXTENSION=I ICOLUMN=K+N/2+1 NEXTENSION=K+1 149 FJAC(ROW(I)+K+N/2+1,ROW(IEXTENSION)+ICOLUMN)=EPS  1  IF(X(ROW(NOUT(IEXTENSION,NEXTENSION))+1) .EQ. 0.0 .OR. PHI(ROW(NOUT(IEXTENSION,NEXTENSION))+1).EQ.0.0) THEN FJAC(ROW(I)+K+N/2+1,ROW(NOUT(IEXTENSION,NEXTENSION))+1)=EPS ELSE  1 1  DALPHADQ=-12*MU0/4/PIPEWIDTH(NOUT(IEXTENSION,NEXTENSION) ,NOUT(NOUT(IEXTENSION,NEXTENSION),1))/(APERTURE( ROW(NOUT(IEXTENSION,NEXTENSION))+1))**2/THAW0  1 1 1 1 1 1 1  DPHIDQ=-((-3*ALPHA(ROW(NOUT(IEXTENSION,NEXTENSION))+1))**(-.5) *COS((ACOS((-27/ALPHA(ROW(NOUT(IEXTENSION,NEXTENSION)) +1)**3)**.5)+PI)/3)-(-ALPHA(ROW(NOUT(IEXTENSION ,NEXTENSION))+1)/3)**.5*SIN((ACOS((-27/ALPHA(ROW( NOUT(IEXTENSION,NEXTENSION))+1)**3)**.5)+PI)/3) *(-27/ALPHA(ROW(NOUT(IEXTENSION,NEXTENSION))+1)**5)**.5 /(1-ABS(-27/ALPHA(ROW(NOUT(IEXTENSION,NEXTENSION)) +1)**3))**.5)*DALPHADQ  1 1 1 1  FJAC(ROW(I)+K+N/2+1,ROW(NOUT(IEXTENSION,NEXTENSION))+1) =-2*THAW0*DPHIDQ/2/APERTURE(ROW(NOUT (IEXTENSION,NEXTENSION))+1)/PHI(ROW(NOUT(IEXTENSION ,NEXTENSION))+1)**2*X(ROW(NOUT(IEXTENSION,NEXTENSION)) +1+NOUT(NOUT(IEXTENSION,NEXTENSION),6)) ENDIF IF(PHI(ROW(NOUT(IEXTENSION,NEXTENSION))+1) .EQ. 0.0) THEN  139  1  FJAC(ROW(I)+K+N/2+1,ROW(NOUT(IEXTENSION,NEXTENSION)) +1+NOUT(NOUT(IEXTENSION,NEXTENSION),6))=EPS  ELSE FJAC(ROW(I)+K+N/2+1,ROW(NOUT(IEXTENSION,NEXTENSION)) 1 +1+NOUT(NOUT(IEXTENSION,NEXTENSION),6))= 1 2*THAW0/2/APERTURE(ROW(NOUT(IEXTENSION,NEXTENSION))+1) 1 /PHI(ROW(NOUT(IEXTENSION,NEXTENSION))+1) ENDIF IF(FLAG(NOUT(IEXTENSION,NEXTENSION),NOUT(NOUT(IEXTENSION 1 ,NEXTENSION),1),TIMESTEP-1).EQ.'FULL')THEN FJAC(ROW(I)+K+N/2+1,ROW(NOUT(IEXTENSION,NEXTENSION))+1 1 +NOUT(NOUT(IEXTENSION,NEXTENSION),6))=EPS ICOLUMN=1+NOUT(NOUT(IEXTENSION,NEXTENSION),6) IEXTENSION=NOUT(IEXTENSION,NEXTENSION) NEXTENSION=1 GOTO 149 ENDIF ENDIF 150  CONTINUE  6000  CONTINUE  110  RR=0.0 DO 110 II=1,2*ACTIVEPIPES RR=RR+FVEC(II)**2 CONTINUE RR=SQRT(RR) IF (RR .EQ. 0.0) GOTO 435  C  CALCULATING NEW X VALUES CALL LINEAREQUATION(2*ACTIVEPIPES,FJAC,FVEC,S)  220 C C C  DO 220 M=1,2*ACTIVEPIPES X(M)=X(M)-S(M) CONTINUE  CALCULATING THE NEW VARIABLE VALUES ROWNUM=0  1  1  DO 6450 I=1,NODENUM ROWNUM=ROW(I) Q0=Q00 IF (I.GT.1 .AND. TIMESTEP.NE.1) Q0=Q(NIN(I), I, TIMESTEP) IF (I.GT.1 .AND. TIMESTEP.EQ.1) Q0=0 IF (I.GT.1 .AND. TIMESTEP.NE.1 .AND. FLAG(NIN(I),I,TIMESTEP-1) .NE. 'FULL') Q0=0 DO 6510 L=1,NOUT(I,6) IF (X(L+ROWNUM) .GT. Q0) X(L+ROWNUM)=Q0 IF (X(L+ROWNUM) .LT. 0) then X(L+ROWNUM)=PIPEWIDTH(I,NOUT(I,L))*APERTURE(ROWNUM+L)**2/THAW0 /MU0*10e-12 end if Q(I, NOUT(I,L), TIMESTEP)=X(L+ROWNUM) IF (X(L+ROWNUM+NOUT(I,6)) .LT. 0) X(L+ROWNUM+NOUT(I,6))=0 LENGTHT(I, NOUT(I,L), TIMESTEP)=X(ROWNUM+L+NOUT(I,6))  140  6510 6450  CONTINUE CONTINUE ROWNUM=0 DO 60000 I=1,NODENUM  C C C  NUMBER OF UNKNOWNS AT EACH GRID N=NOUT(I,6)*2 ROWNUM=ROW(I) IF (N .EQ. 0) GOTO 60000  C C  DEFINING THE EQUATIONS  1  Q0=Q00 IF (I.GT.1 .AND. TIMESTEP.NE.1) Q0=Q(NIN(I), I, TIMESTEP) IF (I.GT.1 .AND. TIMESTEP.EQ.1) Q0=0 IF (I.GT.1 .AND. TIMESTEP.NE.1 .AND. FLAG(NIN(I),I,TIMESTEP-1) .NE. 'FULL') Q0=0 FVEC(ROWNUM+1)=-Q0  1  DO 1100 K=1,N/2 FVEC(ROWNUM+1)=FVEC(ROWNUM+1)+X(ROWNUM+K) FVEC(ROWNUM+K+1)=X(ROWNUM+K+N/2)-X0(ROWNUM+K+N/2) -X(ROWNUM+K)/APERTURE(ROWNUM+K)/PIPEWIDTH(I,NOUT(I,K))*DT IF(TIMESTEP .GT. 1) THEN IF(FLAG(I,NOUT(I,K),TIMESTEP-1).EQ.'FULL') FVEC(ROWNUM+K+1)=0.0 X(ROWNUM+K+N/2)=X0(ROWNUM+K+N/2) ENDIF ENDIF  1 1  1 1  970 1 1  1 1  THEN  ALPHA(ROWNUM+K)=-3-12*MU0*X(ROWNUM+K)/4/PIPEWIDTH(I,NOUT(I,K)) /APERTURE(ROWNUM+K)**2/THAW0 PHI(ROWNUM+K)=2*(-ALPHA(ROWNUM+K)/3)**.5 *COS((ACOS((-27/ALPHA(ROWNUM+K)**3)**.5)+PI)/3) IF (K .EQ. 1) GOTO 1100 DP0=2*THAW0/2/APERTURE(ROWNUM+K-1)/PHI(ROWNUM+K-1) *X(ROWNUM+K-1+N/2) DP1=2*THAW0/2/APERTURE(ROWNUM+K)/PHI(ROWNUM+K) *X(ROWNUM+K+N/2) IF(TIMESTEP .EQ. 1) GOTO 990 IF(FLAG(I,NOUT(I,K-1),TIMESTEP-1).EQ.'FULL')THEN IEXTENSION=NOUT(I,K-1) ALPHA1=-3-12*MU0*X(ROW(IEXTENSION)+1)/4 /PIPEWIDTH(IEXTENSION,NOUT(IEXTENSION,1)) /PIPEAPERTURE(IEXTENSION,NOUT(IEXTENSION,1))**2/THAW0 PHI1=2*(-ALPHA1/3)**.5*COS((ACOS((-27/ALPHA1**3)**.5)+PI)/3) DP0=DP0 +2*THAW0/2/PIPEAPERTURE(IEXTENSION,NOUT(IEXTENSION,1)) /PHI1*X(ROW(IEXTENSION)+1+NOUT(IEXTENSION,6)) IF(FLAG(IEXTENSION,NOUT(IEXTENSION,1),TIMESTEP-1).EQ.'FULL')THEN IEXTENSION=NOUT(IEXTENSION,1) GOTO 970 ENDIF ENDIF IF(FLAG(I,NOUT(I,K),TIMESTEP-1).EQ.'FULL')THEN IEXTENSION=NOUT(I,K)  141  980 1 1  60000  ALPHA2=-3-12*MU0*X(ROW(IEXTENSION)+1)/4 /PIPEWIDTH(IEXTENSION,NOUT(IEXTENSION,1)) /PIPEAPERTURE(IEXTENSION,NOUT(IEXTENSION,1))**2/THAW0 PHI2=2*(-ALPHA2/3)**.5*COS((ACOS((-27/ALPHA2**3)**.5)+PI)/3) DP1=DP1 +2*THAW0/2/PIPEAPERTURE(IEXTENSION,NOUT(IEXTENSION,1)) /PHI2*X(ROW(IEXTENSION)+1+NOUT(IEXTENSION,6)) IF(FLAG(IEXTENSION,NOUT(IEXTENSION,1),TIMESTEP-1).EQ.'FULL')THEN IEXTENSION=NOUT(IEXTENSION,1) GOTO 980 ENDIF ENDIF FVEC(ROWNUM+K+N/2)=DP1-DP0 CONTINUE PRESSURE(I)=DP0 CONTINUE  330  NORM=0.0 DO 330 M=1,2*ACTIVEPIPES NORM=NORM+FVEC(M)**2 CONTINUE  1 1  990 1100  NORM=SQRT(NORM) ITERATION=ITERATION+1 IF (ITERATION .GT. 40) THEN TOL0=TOL0*10 TOLR=TOLR*10 ITERATION=1 ENDIF IF (NORM.GT.TOLR*RR+TOL0) GOTO 111 C C*****WRITING NEW VALUES AND DEFINING THE CHANGES IN THE PIPE CONDITIONS C 435 ROWNUM=0 DO 645 I=1,NODENUM ROWNUM=ROW(I) Q(I, NOUT(I,L), TIMESTEP)=X(L+ROWNUM) LENGTHT(I, NOUT(I,L), TIMESTEP)=X(ROWNUM+L+NOUT(I,6)) FLAG(I, NOUT(I,L), TIMESTEP)='FILLIN' IF(LENGTHT(I,NOUT(I,L),TIMESTEP).GE.(1-1e-5)*PIPELENGTH(I,NOUT(I,L)) .AND. NOUT(I,6).GT. 0) THEN FLAG(I, NOUT(I,L), TIMESTEP)='FULL' LENGTHT(I, NOUT(I,L), TIMESTEP)=PIPELENGTH(I,NOUT(I,L)) IF (NOUT(NOUT(I,L),6) .EQ. 0)THEN CALL NETWORKDATA(NOUT(I,L),NODENUM,XNODE(NOUT(I,L)) 1 ,YNODE(NOUT(I,L)),ELEMDATA,NTOTAL,PIPENUM,PIPELENGTH 1 ,PIPEAPERTURE,PIPEWIDTH,NIN,NOUT,ROW,XNODE,YNODE,ACTIVEPIPES) NODENUM=PIPENUM+1 DO 1361 K=1,2*NOUT(NOUT(I,L),6) X(ROW(NOUT(I,L))+K)=0 X(ROW(NOUT(I,L))+NOUT(NOUT(I,L),6)+K)=0 1361 CONTINUE ENDIF ENDIF IF(LENGTHT(I,NOUT(I,L),TIMESTEP).GE.(1-1e-5)*PIPELENGTH(I,NOUT(I,L)) 1 .AND. NOUT(NOUT(I,L),6).EQ.0) 1 FLAG(I, NOUT(I,L), TIMESTEP)='DEAD' IF (I .EQ. 1) GOTO 650 1  142  1 1 650  IF(LENGTHT(I, NOUT(I,L), TIMESTEP).EQ. 0 .AND. FLAG(NIN(I), I,TIMESTEP).NE. 'FULL') FLAG(I, NOUT(I,L), TIMESTEP)='EMPTY' CONTINUE  645  CONTINUE  C C*****WRITING THE OUTPUT DATA C  1 1 1 1  DO 6449 I=1,NODENUM DO 6448 L=1,NOUT(I,6) WRITE(20,107)TIMESTEP,I,NOUT(I,L),XNODE(I),YNODE(I) ,XNODE(NOUT(I,L)),YNODE(NOUT(I,L)),Q(I,NOUT(I,L),TIMESTEP) 1 ,LENGTHT(I,NOUT(I,L),TIMESTEP),PRESSURE(I) ,FLAG(I,NOUT(I,L),TIMESTEP) WRITE(*,107) TIMESTEP,I,NOUT(I,L),XNODE(I),YNODE(I) ,XNODE(NOUT(I,L)),YNODE(NOUT(I,L)),Q(I,NOUT(I,L),TIMESTEP) 1 ,LENGTHT(I,NOUT(I,L),TIMESTEP),PRESSURE(I) ,FLAG(I,NOUT(I,L),TIMESTEP)  6448 6449  CONTINUE CONTINUE WRITE (20,*) WRITE(*,*)  107  FORMAT(I3,2I4,4F6.2," ",G9.3," ",G9.3," ",G9.3 ,A15)  C C C 64511 649  OMITING THE DEAD PIPES FROM THE NETWORK DO 6451 I=1,NODENUM L=1 IF (NOUT(I,L) .EQ. 0) GOTO 6451 IF (FLAG(I, NOUT(I,L), TIMESTEP) .EQ. 'DEAD')THEN Q(I, NOUT(I,L), TIMESTEP)=0.0 LENGTHT(I, NOUT(I,L), TIMESTEP)=PIPELENGTH(I,NOUT(I,L))  DO 6500 LL=L,NOUT(I,6) NOUT(I,LL)=NOUT(I,LL+1) 6500 CONTINUE DO 6505 LL=ROW(I)+L,2*ACTIVEPIPES X(LL+NOUT(I,6))=X(LL+NOUT(I,6)+1) X(LL)=X(LL+1) X0(LL+NOUT(I,6))=X0(LL+NOUT(I,6)+1) X0(LL)=X0(LL+1) 6505 CONTINUE NOUT(I,6)=NOUT(I,6)-1 L=L-1 DO 6452 LL=2,NODENUM IF (ROW(LL) .GT. ROW(I).AND. LL .NE. I) ROW(LL)=ROW(LL)-2 6452 CONTINUE IF (NOUT (I,6) .EQ. 0 .AND. I .GT. 1)THEN FLAG(NIN(I),I,TIMESTEP)='DEAD' GOTO 64511 ENDIF IF(I .EQ. 1 .AND. NOUT(I,6) .EQ. 0)THEN WRITE(*,*) "ALL PIPES ARE FULL, STOP GROUTING!!!!!!!" READ (*,*) STOP ENDIF  143  ENDIF  6451  L=L+1 IF (L .LE. NOUT(I,6)) GOTO 649 CONTINUE  C*****FINDING THE NUMBER OF THE ACTIVE PIPES TO HAVE THE NUMBER OF UNKNOWNS ACTIVEPIPES=0 DO 6455 I=1,NODENUM ACTIVEPIPES=ACTIVEPIPES+NOUT(I,6) 6455 CONTINUE 50 CONTINUE read(*,*) STOP END  SUBROUTINE LINEAREQUATION(N,COEF,FVEC,X) C  THIS SUBROUTINE CALCULATES THE ROOTS OF A LINEAR SYSTEM OF EQUATIONS INTEGER N,M(670) DOUBLEPRECISION A(670,670),COEF(670,670),FVEC(670),X(670) REAL C  C C***** A IS THE MATRIX WHICH INCLUDES ALL THE COEFFICIENTS AND F VALUES C***** M IS THE ARRAY WHICH KEEPS TRACK OF THE ROWS WHICH ARE MOVED C DO 10 I=1,N DO 20 J=1,N A(I,J)=COEF(I,J) 20 CONTINUE A(I,N+1)=FVEC(I) 10 CONTINUE  11  121  30  DO 1 I=1,N ISUBSTITUTE=I+1 C=A(I,I) IF (C .EQ. 0) THEN DO 121 MM=I,N+1 C= A(I,MM) A(I,MM)= A(ISUBSTITUTE,MM) A(ISUBSTITUTE,MM)= C CONTINUE M(ISUBSTITUTE)=I M(I)= ISUBSTITUTE ISUBSTITUTE=ISUBSTITUTE+1 GOTO 11 ENDIF DO 30 K=I,N+1 A(I,K)=A(I,K)/C CONTINUE IF (I.LT.2) GOTO 2 DO 40 J=1,I-1 C=A(J,I) DO 50 K=I,N+1  144  50 40 2  70 60 1  80  A(J,K)=-A(I,K)*C+A(J,K) CONTINUE CONTINUE DO 60 J=I+1,N C=A(J,I) DO 70 K=I,N+1 A(J,K)=-A(I,K)*C+A(J,K) CONTINUE CONTINUE CONTINUE DO 80 I=1,N X(I)=A(I,N+1) CONTINUE RETURN END  C C @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ @@@@@@@@@ C SUBROUTINE NETWORKDATA(STARTPOINT,NODENUM,XI,YI,ELEMDATA,NTOTAL 1 ,PIPENUM,PIPELENGTH,PIPEAPERTURE,PIPEWIDTH,NIN,NOUT,ROW 1 ,XNODE,YNODE,ACTIVEPIPES) C C*****THIS SUBROUTINE READS AND UPDATES THE NETWORK DATA FROM THE NPUT DATA C C DOUBLE PRECISION XNODE(740),YNODE(740),XI,YI REAL PIPELENGTH(740,740),PIPEAPERTURE(740,740) 1 ,PIPEWIDTH(740,740),ELEMDATA(801,10) INTEGER NTOTAL,I,J,K,NOUT(800,6),NIN(740),PIPENUM,STARTPOINT 1 ,NODENUM,ROW(800),MAXROW,ACTIVEPIPES  K=1 I=STARTPOINT J=NODENUM+1 XNODE(I)=XI YNODE(I)=YI DO 10 N=1, NTOTAL  1  IF (ELEMDATA(N,1).EQ.XI .AND. ELEMDATA(N,2).EQ.YI .AND. ELEMDATA(N,9).EQ.0)THEN ELEMDATA(N,9)=I ELEMDATA(N,10)=J XNODE(J)=ELEMDATA(N,3) YNODE(J)=ELEMDATA(N,4) PIPELENGTH(I,J)=ELEMDATA(N,5) PIPEAPERTURE(I,J)=ELEMDATA(N,6) PIPEWIDTH(I,J)=ELEMDATA(N,7) NIN(J)=I NOUT(I,K)=J K=K+1 J=J+1  1  END IF IF (ELEMDATA(N,3).EQ.XI .AND. ELEMDATA(N,4).EQ.YI .AND. ELEMDATA(N,9).EQ.0)THEN ELEMDATA(N,9)=I ELEMDATA(N,10)=J  145  XNODE(J)=ELEMDATA(N,1) YNODE(J)=ELEMDATA(N,2) PIPELENGTH(I,J)=ELEMDATA(N,5) PIPEAPERTURE(I,J)=ELEMDATA(N,6) PIPEWIDTH(I,J)=ELEMDATA(N,7) NIN(J)=I NOUT(I,K)=J K=K+1 J=J+1 END IF 10  51  CONTINUE NOUT(I,6)=K-1 PIPENUM=PIPENUM+NOUT(I,6) ACTIVEPIPES=ACTIVEPIPES+NOUT(I,6) MAXROW =0 IMAXROW=1 DO 51 N=1,NTOTAL IF (ROW(N) .GT. MAXROW)THEN MAXROW=ROW(N) IMAXROW=N ENDIF IF (ROW(N).EQ.MAXROW .AND. NOUT(N,6).GE.NOUT(IMAXROW,6))THEN MAXROW=ROW(N) IMAXROW=N ENDIF CONTINUE RETURN END  146  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0063121/manifest

Comment

Related Items