USFDC Home  USF Electronic Theses and Dissertations   RSS 
Material Information
Subjects
Notes
Record Information

Full Text 
PAGE 1 Analysis of Water Seepage Through Earthen St ructures Using the Particulate Approach by Kalyani Jeyisanker A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy Department of Civil and Environmental Engineering College of Engineering University of South Florida Major Professor: Manjriker Gunaratne, Ph.D. Alaa Ashmawy, Ph.D. Mark Luther, Ph.D. Gray Mullins, Ph.D. Muhammad Mustafizur Rahman, Ph.D. Date of Approval: November 03, 2008 Keywords: Navier Stokes Equations, MonteCarlo Simulation, Volumetric Compatibility, Staggered Grid, Critical Localized Condition, Unsaturated, Soil Water Characteristic Curve Copyright 2008, Kalyani Jeyisanker PAGE 2 DEDICATION This work is dedicated to my husband, Jeyisanker and my daughter, Harinie. Jeyisanker has been a tremendous support. He was always there in the hard days throughout the research. Without his support a nd encouragement, I would never achieve this. PAGE 3 ACKNOWLEDGMENTS I would like to thank my advisor, Dr. Manj riker Gunaratne, for his guidance and support throughout the research. Dr. Guna ratne was always there to suggest ideas and to improve my writing skills, and to expand my academ ic knowledge in Geotechnical Engineering. Additionally, I would like to extend my thanks to the members of my examining committee Dr. Alaa Ashmawy, Dr. Gray Mullin s, Dr. Mark Luther and Dr. Muhammad Mustafizur Rahman for all their support a nd courage given to me to complete my research. I would like to thank Dr. Nimal Se neviratne for his advice to complete my studies successfully. My thanks go to my colleagues, Ms. S.Amar asiri, Mr. L.Fuentes, Mr. M.Rajapakse, Mr. J.Kosgolle and Ms. K.Kranthi for their help especially when I am out of office. PAGE 4 i TABLE OF CONTENTS LIST OF TABLES ...............................................................................................................v LIST OF FIGURES ........................................................................................................... vi ABSTRACT ...................................................................................................................... ...x CHAPTER ONE: INTRODUCTION ..................................................................................1 1.1 Background ........................................................................................................1 1.2 Literature Review...............................................................................................2 1.3 Objectives .........................................................................................................4 1.3.1 Pavement Filter Design .......................................................................4 1.3.2 Retention Pond Design .......................................................................6 1.4 Numerical Modeling ..........................................................................................8 1.4.1 Solid Continuum Methods ..................................................................8 1.4.2 Discrete Element Methods ..................................................................8 1.5 Assembling of Particles in the Particle Flow Code ...........................................9 1.5.1 Theoretical Background of PFC .......................................................10 1.5.1.1 Particle Interactions ...........................................................10 1.6 Organization of Dissertation ............................................................................11 CHAPTER TWO: ANALYSIS OF S TEADY STATE WA TER SEEPAGE IN A PAVEMENT SYSTEM USING THE PARTICULATE APPROACH .......12 2.1 Preliminary Studies Using Existing Software .................................................12 2.2 Methodology Followed to Develop a Novel Algorithm ................................ 14 2.2.1 Modeling of Soil Structure (3D Random Packing of Soil Particles) ..........................................................................................17 2.2.1.1 Simulation of Maximum and Minimum Void Ratios Using PSD ..............................................................17 PAGE 5 ii 2.2.1.2 Implementation of the Packing Procedure .........................18 2.2.1.3 Determination of the Natura l Void Ratio Distribution ......20 2.3 Modeling of Steady State Flow .......................................................................23 2.3.1 Need for a Staggered Grid ................................................................25 2.3.2 Numerical Solution Technique .........................................................28 2.3.3 SemiImplicit Method for PressureLinked Equations (SIMPLE) Algorithm ........................................................................31 2.3.3.1 Derivation of the Pressure Correction Formula .................32 2.3.3.2 Boundary Conditions for Pavement System ......................36 2.4 Validation of the Numerical Model .................................................................37 2.5 Numerical Illustration ......................................................................................39 2.5.1 Results of Numerical Illustration ......................................................41 2.6 Conclusions ......................................................................................................43 CHAPTER THREE: ANALYSIS OF TR ANSIENT WATER SEEPAGE IN A PAVEMENT SYSTEM USING TH E PARTICULATE APPROACH ................44 3.1 Modeling of Transient Flow ............................................................................44 3.1.1 Transient Navier Stokes Equations ...................................................46 3.1.2 Volumetric Compatibility of Solid and Fluid Phases .......................47 3.1.3 Numerical Solution Technique .........................................................48 3.1.4 Boundary Conditions ........................................................................52 3.2 Interface Effect in FDM ...................................................................................53 3.3 Validation of the Numerical Models ...............................................................54 3.4 Numerical Illustration ......................................................................................56 3.5 Selection of Time Step .....................................................................................56 3.6 Results of Numerical Illustration .....................................................................57 3.7 Limitations of the Model .................................................................................65 3.8 Conclusions ......................................................................................................66 CHAPTER FOUR: TRANSIENT SEEPAGE MODEL FOR PARTLY SATURATED AND SATURATED SOILS USING THE PARTICULATE APPROACH ..............................................................................68 4.1 Introduction ......................................................................................................68 PAGE 6 iii 4.2 Existing Particulate Models .............................................................................69 4.3 Flow in Partly Saturated Soils .........................................................................70 4.4 Overview of the New Model 73 4.5 Modeling of Soil Structure (3D Ra ndom Packing of Soil Particles) ...............76 4.5.1 Simulation of Maximum and Minimum Void Ratios Using PSD .........................................................................................76 4.5.2 Implementation of the Packing Procedure ........................................77 4.5.3 Determination of the Natura l Void Ratio Distribution .....................79 4.6 Modeling of Fluid Flow in Pa rtially Saturated Soils ...................................... 82 4.6.1 Soil Water Characteristic Curve/Water Moisture Retention Curve ................................................................................83 4.6.2 Determination of the Unsaturated Permeability Function .............. 85 4.7 Development of the Analytical Model .............................................................86 4.7.1 Navier Stokes Equations ...................................................................86 4.7.2 Modification of Drag For ce for Water Flow through Unsaturated Soils ..............................................................................88 4.8 Numerical Solution Technique ........................................................................89 4.9 Numerical Applications of the Model..............................................................97 4.9.1 Modeling Confined Flow thro ugh a Partly Saturated Fine Sand Layer (Case I) ..........................................................................97 4.9.2 Boundary Conditions ..........................................................................98 4.9.3 Results ................................................................................................99 4.10 Conclusions ..................................................................................................102 CHAPTER FIVE: MODELING UNC ONFINED FLOW AROUND A RETENTION POND ...........................................................................................103 5.1 Introduction ....................................................................................................103 5.2 Existing Design of Retention Pond ................................................................104 5.3 Modeling of Soil Structure (3D Ra ndom Packing of Soil Particles) .............106 5.4 Numerical Procedure for Determ ination of the FreeSurface ........................107 5.4.1 Navier Stokes Equations .................................................................109 5.4.2 Boundary Conditions ......................................................................110 PAGE 7 iv 5.5 Determination of FreeSurf ace from Dupuits Theory ..................................112 5.6 Results on FreeSurface .................................................................................112 5.7 Prediction of Recovery Time .........................................................................114 5.8 Conclusions ....................................................................................................118 REFERENCES ................................................................................................................119 ABOUT THE AUTHOR ................................................................................... END PAGE PAGE 8 v LIST OF TABLES Table 1: Size Characteristics of Pavement Layers .............................................................19 Table 2: Comparison of 3D Coefficients of Hydraulic Conductivity Derived from SIMPLE Algorithm .....................................................................................38 Table 3: Model Parameters ................................................................................................40 Table 4: Comparison of 3D Coefficients of Hydraulic Conductivity Derived from the Numeri cal Methodology .....................................................................54 Table 5: Soil Characteristics ..............................................................................................78 Table 6: Parameters Used to Define the SWCC for Sand [23] ..........................................84 PAGE 9 vi LIST OF FIGURES Figure 1: Typical Pavement Structure ...............................................................................5 Figure 2: Pavement Layer Interfaces .................................................................................5 Figure 3: Computational Cycle in PFC ............................................................................10 Figure 4: Illustration of InterParticle Forces ..................................................................10 Figure 5: Analysis of Flow Problem Using Existing Software ........................................13 Figure 6a: Highway Pavement Layers ...............................................................................13 Figure 6b: Comparison of Hydraulic Gr adients Obtained Using Continuum and Discrete Methods (Steady State) ...............................................................14 Figure 7: Flow Chart Illustrating the Analytical Model for Steady State Flow ...............16 Figure 8: Particle Size Distributions fo r Gravel, Coarse Sand and Fine Sand ...............18 Figure 9: Random Packing of Soil Particles ....................................................................20 Figure 10: Probability Density Functi on (PDF) of Maximum and Minimum Void Ratios ..................................................................................................... 22 Figure 11: Cumulative Frequencies of Maximum and Minimum Void Ratios .................23 Figure 12: Discrete Checkerboard Veloci ty Distribution at Each Grid Point: At Each Node, the Top Number is u and the Bottom Number is v .................25 Figure 13: Discrete Checkerboard Pressure Distribution ..................................................26 Figure 14: A Typical Staggered Grid A rrangement (Solid Circles Represent Pressure Nodes and Open Circle s Represent Velocity Nodes) ........................28 Figure 15: Computational Module for the XMomentum Equation ..................................29 Figure 16: Flow Chart for SIMPLE Al gorithm for Steady State Condition ......................32 Figure 17: Designation of Nodal Points on a Grid Used in SIMPLE Algorithm ..............34 Figure 18: Boundary Conditions Incor porated in the Flow through the Pavement System .............................................................................................37 PAGE 10 vii Figure 19: Relationship between Hydr aulic Conductivity Vs Diameter for Uniformly Distributed Soils ...........................................................................38 Figure 20: Plot of Flow Rate thr ough Coarse SandGrav el Interface Vs Iteration Steps with Differ ent Levels of Compaction (For Pressure Differential of 23 kPa) .............................................................41 Figure 21: Plot of Flow Rate throug h Fine SandCoarse Sand Interface Vs Iteration Steps with Different Leve ls of Compaction (For Pressure Differential of 23 kPa) ....................................................................................41 Figure 22: Comparison of Hydraulic Gr adients Obtained Using Continuum and Particulate Approaches (Usi ng SIMPLE Algorithm Steady State) ...........42 Figure 23: Plot of Flow Rate Vs Iteration Step for High Hydraulic Gradient (i = 1.7 icr) ......................................................................................................42 Figure 24: Flow Chart Illustrating th e Analytical Model for Transient and Steady State Flow ...........................................................................................45 Figure 25: Computational Grid for the XMomentum Equation ......................................49 Figure 26: Designation of Nodal Points on a Grid Used in the Iterative Procedure .........50 Figure 27: Boundary Conditions Incorporat ed in the Flow through the Pavement System (Saturated Soil)....................................................................................53 Figure 28: Development of Water Pressu re within the Coarse SandGravel Layers with Time (For Pressure Differential of 23 kPa) .................................57 Figure 29: Development of Water Veloci ty within the Coarse SandGravel Layers with Time (For Pressure Differential of 23 kPa) .................................58 Figure 30a: Plot of Flow Rate through Coarse SandGravel Layers Vs Time Steps with Different Levels of Compaction (For Pressure Differential of 23 kPa) .....................................................................................58 Figure 30b: Plot of Flow Rate through Fine SandCo arse Sand Layers Vs Time Steps with Different Levels of Compaction (F or Pressure Differential of 23 kPa) .....................................................................................59 Figure 31a: Plot of Maximum Hydraulic Gradient within the Coarse SandGravel Layers Vs Time Steps fo r Different Levels of Compaction (For Pressure Differential of 23 kPa) ...............................................................60 Figure 31b: Plot of Differe nce in Maximum Hydraulic Gradients Within the Coarse SandGravel Layers Vs Ti me Steps for Different Levels of Compaction (For Pressure Differential of 23 kPa) ....................................61 PAGE 11 viii Figure 32: Comparison of Hydraulic Gradients Obtained Using Continuum and Particulate Approach es (Steady State) ....................................................62 Figure 33: Porosity Distribution fo r Coarse SandGrave l Interface with NonUniform Compaction (Circl ed Area is Packed with Dr of 5% and the Remaining Area is Packed with Dr of 95%) ...........................63 Figure 34: Velocity Vector Plots: (a) Uniformly Compacted Soil Media (Dr = 95%) (b) Deficiency of Co mpaction in Local Area .............................64 Figure 35: Plot of Flui d Pressure Distribution acro ss the Coarse SandGravel Interface: (a) Uniformly Compacted Coarse SandGravel Interface (Dr = 95%) (b) Deficiency of Compaction in Local Area .............................64 Figure 36: Plot of Localized H ydraulic Gradients for Coarse SandGravel Interface with Weak Region ...............................................................65 Figure 37: Schematic Diagram to Show the Movement Cycle of Water from Atmosphere to the Groundwater Table [20] ..........................................71 Figure 38: Flow Chart Illustrati ng the Analytical Model for Flow through Unsaturated Soils ..............................................................................75 Figure 39: Particle Size Distribution for Fine Sand .........................................................76 Figure 40: Packing Procedure ..........................................................................................79 Figure 41: Cumulative Di stribution for Natural Void Ratios for the Fine Sand with Dr = 50% ........................................................................................81 Figure 42: Initialization of Water Pre ssure in the Saturated and Unsaturated Zones ..............................................................................................................82 Figure 43a: Typical Soil Wate r Characteristic Curve and wm2for a SaturatedUnsaturated Soil [10] ......................................................................................83 Figure 43b: Soil Water Characteristic Cu rve for Specific Soil Types [23] .......................83 Figure 44: Unsaturated Permeability Function for Selected Soil Types [23 ] ..................86 Figure 45: Computational Grid for the XMomentum Equation .....................................90 Figure 46: Soil Water Charact eristic Curve for Fine Sand ..............................................93 Figure 47: Neighborhood Points Used in the Iterative Procedure ..................................94 Figure 48: Illustration of the Flow through the Unsa turated Fine Sand Layer ................97 Figure 49: Boundary Conditions Incorpor ated in the Flow through Saturated Unsaturated Soil Layer ...................................................................................98 Figure 50: The Variation of Water Poro sity with Water Pressure at a Point within the Unsaturated Zone (At Point A in Figure 49) .................................99 PAGE 12 ix Figure 51: The Variation of Degree of Sa turation with Water Pressure (At Point A in Figure 49) .............................................................................................100 Figure 52: Velocity Ve ctor Plots for Flow through Saturated and PartlySaturated Soil Layer ....................................................................................101 Figure 53: Variation of Discharg e Flow Rate with Time Step ......................................102 Figure 54: Typical Retention Pond ...............................................................................103 Figure 55: Particle Size Di stributions for Coarse Sand .................................................107 Figure 56: Determination of Dry a nd Wet Zones around the Retention Pond ..............109 Figure 57: Boundary and Initial Conditi ons Incorporated in the Flow around the Retention Pond .......................................................................................111 Figure 58: Analytical Method to Determine the FreeSurface ......................................112 Figure 59: Comparison of FreeSurface around the Retention Pond .............................113 Figure 60: Location of FreeSurface around the Retention Pond with Time ................114 Figure 61: Saturated Infiltra tion Analysis Using MODFLOW .....................................115 Figure 62: Prediction of Recovery Time from Numerical Method ...............................116 Figure 63: Sample Plot to Determ ine Recovery Time from MODFLOW ....................117 PAGE 13 x ANALYSIS OF WATER SEEPAGE THROUGH EARTHEN STRUCTURES USING THE PARTICULATE APPROACH KALYANI JEYISANKER ABSTRACT A particulate model is developed to analy ze the effects of steady state and transient seepage of water through a randomlypacked co arsegrained soil as an improvement to conventional seepage analysis based on con tinuum models. In the new model the soil skeleton and pore water are volumetrically c oupled. In the first phase of the study, the concept of relative density has been used to define different compaction levels of the soil layers of a completely saturated pavement filter system and observe the seepage response to compaction. First, MonteCarlo simulation is used to randomly pack discrete spherical particles from a specified Particle Size Dist ribution (PSD) to achie ve a desired relative density based on the theoretical minimum and maximum void ratios. Then, a water pressure gradient is applied across one twola yer filter unit to trigger water seepage. The pore water motion is idealized using Navier Stokes (NS) equations which also incorporate drag forces acting between the wa ter and soil particles. The NS equations are discretized using finite differences and app lied to discrete elements in a staggered, structured grid. The model predicted hydrau lic conductivities are va lidated using widely used equations. The critical water veloci ties, hydraulic gradients and flow within the PAGE 14 xi saturated soil layers are identified under both steady st ate and transient conditions. Significantly critical transient conditions seem to develop. In the second phase of the study the model is extended to analyze the confined flow through a partly saturated pavement layer and unconfined flow from a retention pond into the surrounding saturated granular soil medi um. In partly saturated soil, the water porosity changes resulting from water flow is updated using the Soil Water Characteristics Curve (SWCC) of the soil. The results show how complete saturation develops due to water flow following the wate r porosity Vs pressure trend defined by the SWCC. Finally, the model is used to predict the gradual reduction in the water level of a retention pond and the location of the freesurface. The freesurface is determined by differentiating the wet and dry zones based on the Heaviside step function modified NS equations. PAGE 15 1 CHAPTER ONE INTRODUCTION 1.1 Background Durability of earthen structures such as dams, levees, embankments and pavements is determined by one dominant factor; the nature of interaction of soil particles with water flow. Hence accurate analysis of water seepag e through soils is essential to achieve more durable designs of such structures. The major ity of currently availa ble design criteria are formulated based on either the analysis of steady state laminar flow through saturated soil continua or empiricism. However, very often, field observations are also used to refine or calibrate the design criteria. In the conve ntional models, the dyna mic flow of water through soil pores is commonly idealized usi ng the Darcys law. E xperimental studies show that Darcys law could be inaccurate for modeling transient conditions and high fluid velocities that develop under excessive hyd raulic gradients [1]. It is also known that, under wind and tidal impacts as well as rain fall and rapid reservoir drawdown, it is the transient and nonlaminar flow that plays a more crucial role in dete rmining the stability of earthen and hydraulic infrastructure. In orde r to evaluate localized critical zones, one has to replace the conventional method of anal ysis based on a continuum to an alternative approach with a discrete soil skeleton which allows passage of water through its PAGE 16 2 interstices. Moreover, forensic investiga tions of failures often remind the civil engineering community of 1) The vital role of the discontinuous or particulate nature of soil 2) The importance of analyzing the flow through unsaturated soils, and 3) The importance of incorporating critical transient effects that generally precede the eventual steady state flow. 1.2 Literature Review Modeling of seepage through particulate medi a considering soilwater interaction is relatively new to computational geomechanics. Due to its complexity, Fredlund [2] used Richards equation (Eqn 1) and the continuum approach to obtain approximate solutions for slope stability problems. t h m y h y k x h x k y h k x h kw y x y x 2 2 2 2 2 (1) Where, h is the total hydraulic head, kx and ky are the x and y directional hydraulic conductivities and mw 2 is the water storage coefficient equal to the slope of SoilWater Characteristic Curve. As for nonsteady state or transient flow pr oblems, Fredlund [2] used Richards equation (Eqn 1) and the continuum approach to obtain approximate solutions. Where, h is the total hydraulic head, kx and ky are the x and y directional hydraulic conductivities and mw 2 is the water storage coefficient equal to the slope of SoilWater Characteristic Curve [2]. Ng and Shi [3] also used Eqn1 to numerically investigate the PAGE 17 3 stability of unsaturated soil slopes subjected to transient seepage. In Ng and Shis [3] work, a finite element model was used to i nvestigate the influence of various rainfall events and initial ground water conditions on transient seepage. However, slope stability was analyzed without consider ing the localized e ffects of high pressure buildup and high hydraulic gradients within the slope. Modeling of seepage through particulate medi a considering soilwater interaction is relatively new to computational geomechan ics. The discrete element method (DEM) provides an effective tool to model granular soils in particular based on micro mechanical idealizations. El Shamy et al. [4] presented a computational micromechanical model for coupled analysis of pore water flow and deform ation of granular assemblies. El Shamy et al. [4] have used their model and investigat ed the validity of Darcys law in terms of particle sizes and porosities. In addition, El Shamy and Zeghal [5] conducted simulations to investigate the three dimensional response of sandy deposits when subjected to critical and overcritical upward pore fluid flow us ing a coupled hydromechanical model. These simulations provide valuable information on a number of salient microscale mechanisms of granular media liquefaction under quick sand conditions. In addition, Shimizus [6] particlefluid coupling scheme with a mi xed LagrangianEulerian approach which enables simulation of coupling problems with large Reynolds numbers is implemented in PFC 2D and PFC 3D released by Itasca Cons ulting Group, Inc. [7]. The models used by El Shamy and Zeghal [5] and Shimizu [6] are both based on the work by Anderson and Jackson [8] and Tsuji et al. [9]. Anders on and Jackson [8] mode led pore fluid motion using averaged Navier Stokes equations. Tsuji et al. [9] simulated the process of particle PAGE 18 4 mixing of a twodimensional gasfluidized be d using averaged Navier Stokes equations for comparison with experiments. For all the above cited studies, granular assemblies are modeled using the discrete element model de veloped by Cundall and Strack [10] and the averaged Navier Stokes equations are disc retized using a finite volume technique on a staggered grid [11]. The discrete nature of soil makes the required constitutive relationships to be exceedingly complex needing a large number of parameters to be evaluated in or der to model the soil behavior accurately. However, the stateoftheart hi gh performance comp uter facilities would help the designer save time on th e computations. Furthermore, nanoscale experimentation can be performed to establish model parameters such as the coefficients of normal and shear stiffness between the grains. 1.3 Objectives 1.3.1 Pavement Filter Design Design of durable filters is essential for highway pavements since the filters largely determine the success and failure of the drai nage system and the lasting separation of pavement layers. Inadequate compaction or se gregation of filter layers during placement and excessive cyclic traffic loads can lead to undesired soil particle migration and eventual erosion causing the malfunction of the pavement. Figure 1 shows a typical pavement structure made of three layers ; subgrade, subbase and base. The results presented in the paper are only for the subba sebase layer interface (Figure 2a). Similar results can be obtained for the subgradesubbase layer interface as well (Figure 2b). PAGE 19 5 Currently, the conventio nal criteria (Eqn 2) proposed by U.S. Army Corps of Engineers [12] are used for filter design. In these criter ia, the filter refers to coarser layer and the finer layer is defined as the soil. Figure 1: Typical Pavement Structure (a) (b) Figure 2: Pavement Layer Interfaces Clogging criterion: 585 15soil filterD D (2a) Permeability criterion: 515 15soil filterD D (2b) Additional criterion: 2550 50soil filterD D (2c) Base (Example: Gravel) Subbase (Example: Coarse sand) Subgrade (Example: Fine sand) Gravel (Base) Coarse sand (Subbase) 500 mm (20 in.) 750 mm (30 in.) Coarse sand (Subbase) Fine sand (Subgrade) PAGE 20 6 D15, D50 and D85 are the diameters of a soil mixt ure which correspond to cumulative weight percentages of 15%, 50% and 85% resp ectively in a particle size distribution of the mixture as shown in Figure 8. While the above criteria generally enable the designer to select the gradation of different structural layers, the flow characteristics a nd hydraulic gradients within the system are determined using the Darcys law. An obvious inadequacy of the current drainage design is the use of flow parameters and hydraulic conditions that represent only the overall average and steadystate conditions of each la yer. Therefore, the conventional drainage and filter design techniques are unable to incorporate, perhaps more critical, localized, random and transient effects. Hence, the m odel presented here would equip designers with an analytical tool to address the deficiencies of curre nt design techniques. 1.3.2 Retention Pond Design Retention ponds are manmade or natural de pressions into which stormwater runoff is directed for temporary storage with the exp ectation of disposal by infiltration into a shallow groundwater aquifer. They are often created near areas of development and in many instances required with new development of buildings, parking lots, roads, etc by the permitting agencies. Retention ponds are de veloped primarily to serve two functions such as limit flooding and removal of pollutants. These ponds generally comprise a sedimentation forebay and a larger basin sized to hold the water quality volume (WQV). They retain larger storm volumes for 24 to 48 hours, thus protecting the channels (streams, etc.) th at receive the effluent. They also can be PAGE 21 7 designed to retain larger volumes genera ted by 10to 100year rain events. Water treatment is achieved naturally when partic les settle along the flow path between inlet and outlet of the pond, and between storms wh en additional settling occurs. Nutrient removal occurs between storms via plant upta ke. Rain events provide a fresh influx of stormwater runoff, which forces standi ng water out of the system. Maintenance requirements of retention ponds include the pe riodic removal of sediment and vegetation to restore storage capacity. Sediment remova l is performed primarily in the forebay, which can be designed for easy equipment access. The model presented in this dissertation firs t uses a selfdevelope d packing algorithm to randomly pack a three dimensional discrete soil skeleton. Then the model is used to determine the water flow behavior of th e particulate soil medium consisting of volumetrically coupled water co ntinuum and the disc rete soil skeleton. The flow of water through the particulate medium is modeled using the Navier Stokes (NS) equations which are discretized using the fin ite difference method (FDM) [13] The new model is capable of predicting both transient and steady state flow effects. The model is applied to a pavement structure designed based on the U.S. Army Corps of Engineers filter designed criteria [12] to analyze the lo calized transient and steady stat e seepage effects in terms of water velocities and hydraulic gradients at different degree s of compaction. Appropriate boundary conditions have been employed to simu late the conditions resulting from of a sudden surge of ground water just beneath the pavement subbase. It is assumed in the new model that reasonably accurate estimates of the above parameters can be obtained by volumetric coupling of wate r and the soil skeleton. PAGE 22 8 1.4 Numerical Modeling Numerical modeling can be performed to solve geomechanics problems using two different approaches: 1) Solid continuum methods 2) Discrete element methods 1.4.1 Solid Continuum Methods In the solid continuum approach, the entire so il body is first divide d into a number of small elements and the governing equations are mathematically solved for each element. In this regard, the following th ree methods are widely used: 1) Finite Difference method 2) Finite Element method 3) Finite Volume method 1.4.2 Discrete Element Methods Discrete element methods (DEM) comprise a su ite of numerical techniques developed to model granular materials, rock, and other discon tinua at the scale of grains. In most cases, the granular particles modeled as 3D spheres or 2D discs are indivi dually packed in the structure. Then, appropriate in terparticle characteristics such as coefficients of normal and shear stiffness, friction between adjacent particles and friction between particles and other structures are introduced in the analys is. As in the case of continuum methods, the governing equations are solved numerically. Due to the nature of this analysis, DEM are also known as particle modeling methods. Part icle Flow Code (PFC) [7] is one of the discrete element methods. PAGE 23 9 1.5 Assembling of Particles in the Particle Flow Code Using the PFC program, individual soil part icles can be packed to model a given geotechnical structure such as an earthen da m or a pavement layer by closely simulating the transient dynamics of that particulate me dium involved with the construction of that structure. The following laws govern the p acking (construction) mechanism to achieve the final force equilibrium. 1) Law of motion (Newtons second law) 2) Appropriate force displacement (constitutive) laws for normal and shear deformation The computational cycle in PFC2D is a time stepping algorithm that consists of the repeated application of the la w of motion to each particle, a forcedisplacement law to each contact, and constant updating of boundari es. The forcedisplacement law based on the contact constitutive model is repeatedly appl ied to each contact to update the contact forces based on the relative motion between the two entities at the contact. Until the particles reach equilibrium, the above laws will be applied in a loop as shown in Figure 3. PAGE 24 10 Figure 3: Computational Cycle in PFC 1.5.1 Theoretical Background of PFC 1.5.1.1 Particle Interactions Figure 4: Illustration of InterParticle Forces The abovementioned rigid circular particles (or spherical in 3D) interact by way of normal and shear contacts modeled by a simplified massspring system as shown in Figure 3. For each pair of particles, the in teractive force can be written as (3a) s s k N N k Fs n i. S (Shear) N (Normal) PAGE 25 11 Where, kN Coefficient of normal stiffness, ks Coefficient of shear stiffness, N Normal deformation and s Shear deformation. Then, the total force acting on a particle given by, i totalF F (3b) 1.6 Organization of Dissertation Chapter 2 describes the methodology used to randomly pack the soil particles and the mathematical formulation for steady state fl ow through a particulate pavement system. Chapter 3 presents the modification made on the governing equations for analyzing the transient behavior of water through particul ate pavement system using the volumetric compatibility of the two phase media, such as water and solid particles. Chapter 4 describes the mathematical formulation of flow through partlysaturated and saturated soils using Navier Stokes equations. Chapter 5 illustrates the applica tion of the developed models for analyzing flow around retention pond s in order to determ ine the freesurface (Phreatic surface) using Navier Stokes equa tions. The different approaches are proposed to determine the freesurface for flow from retention pond. Results are included in each Chapter. PAGE 26 12 CHAPTER TWO ANALYSIS OF STEADY STATE WATER SEEPAGE IN A PAVEMENT SYSTEM USING THE PARTICULATE APPROACH 2.1 Preliminary Studies Using Existing Software The existing finiteelement software (Seep /W) and discrete element software (PFC2D) were used at the time of preliminary studies in order to visualize the effects of continuum and discrete approaches on fl ow problems. A simple model shown in Figure 5 was used to identify the localized effects within a twolayer pavement system. Under same pressure gradient, the flow rates are compared using continuum and discrete methods. Flow rate obtained using discrete method (0.0025 m3/s/m width) is less than that of continuum method (q = 0.0027 m3/s/m width) due to additional drag forces in the particulate matrix. PAGE 27 13 Figure 5: Analysis of Flow Problem Using Existing Software Figure 6a: Highway Pavement Layers DEM (PFC 2D) FEM (Seep/W) (Contours for pressure head are shown) Subbase Subgrade PAGE 28 14 Figure 6b: Comparison of Hydraulic Gr adients Obtained Using Continuum and Discrete Methods (Steady State) The seepage through twolayer pavement s hown in Figure 6a was modeled using the PFC2D. As illustrated in Figure 6b, there is a si gnificant difference in hydraulic gradients using the continuum and discrete approaches. At the interface, ilocal exceeds the iConventional analysis. However, since the threedimensional porosities are realistic, the authors developed a threedimensional packing al gorithm. The results obtained from the preliminary studies motivated the authors for analyzing the seepage phenomena through particulate soil media. 2.2 Methodology Followed to Develop a Novel Algorithm The comprehensive analytical procedure and the computer code developed for its implementation are illustrated in Figure 7. The analytical proce dure consists of two primary tasks such as random assembly of the particulate medium (granular soil) and solution of the fluid flow governing equatio ns using partial coupling between the two Plot of hydraulic gradient across a twolayer soil 0 20 40 60 80 010203040 Hydraulic gradientDepth (mm) Localized i DEM i across subgradeDarcy i across subbaseDarcy Sub g radesubbase interface PAGE 29 15 media. The flow chart also includes the s ections, equation numbers and figure numbers corresponding to each stage. PAGE 30 16 Figure 7: Flow Chart Illustrating the An alytical Model for Steady State Flow Packing to determine maximum and minimum void ratios for each particle size distribution (Section 2.2.1) Application of MonteCarlo simulation to obtain randomly distributed insitu void ratios at specified relative densities (Section 2.2.1.3) Assignment of random insitu void ratios for grid elements used in the fluid flow analysis Application of NavierStokes equations to each grid element (Section 2.3) Numerical discretization based on a staggered grid (Section 2.3.2 & Figure 15) Solid phase Fluid phase Execution of SIMPLE algorithm (Fig. 16) Analysis of results critical velocities and hydraulic gradients at the steadystate (Section 2.5) Mass imbalance term < tol. Yes Validating the model by determining the coefficients of hydraulic conductivities (Section 2.4) Iterative solution for steadystate No PAGE 31 17 2.2.1 Modeling of Soil Structure (3D Random Packing of Soil Particles) The particle size distribution curves (PSDs) s hown in Figure 8 were selected to satisfy the US. Army Corps of Engineers filter criteria (Eqn 2) for a typical pavement system made of a gravel base, a coarse sand s ubbase and a fine sand subgrade [12]. 2.2.1.1 Simulation of Maximum and Mi nimum Void Ratios Using PSD In this model the soil particles are assumed to be of spherical shape. To determine the maximum and minimum void ratios (emax and emin) corresponding to each soil type in Figure 8, a customized random packing algor ithm was developed. In the case of gravel, 28 mm 28 mm 28 mm cubes were packed using the corresponding PSD in Figure 8. It can be noted that the PSDs in Figure 8 are cumulative probability distributions of particle size for each soil type. Therefore, by using an adequately large array of random numbers from a uniform distri bution between 0 and 100 (the y axis range in Figure 8), one can select the corr esponding array of particle diameter s that conforms to a selected PSD, from the x axis. This technique, known as the MonteCarlo simulation [14] is used to select the array of packing diameters fo r each cube. It is noted that the resulting maximum and minimum void ratio distributions change in each packing trial since the order of diameters used for packing is changed randomly. Similarly, 16 mm 16 mm 16 mm cubes and 3 mm 3 mm 3 mm cubes were used to obtain emax and emin in coarse sand packing and fine sand packing respectivel y (Figure 8). Smaller cubes are used to pack smaller particle sizes to reduce th e computer running time assuming that an PAGE 32 18 adequate number of particles are within the cube to determine representative emax and emin for each type of soil. 0 10 20 30 40 50 60 70 80 90 100 0.010.1110100Grain size/(mm)% Cumulative frequency Coarse sand Gravel Fine sand Figure 8: Particle Size Distributions for Gravel, Coarse Sand and Fine Sand 2.2.1.2 Implementation of the Packing Procedure In order to obtain the looses t state of each type of soil (emax) within the corresponding cubes described above, different sizes of soil particles inscribed in boxes are packed as indicated Figure 9a using a MATLAB code deve loped by the authors. The side length of each box is the same as the diameter of the inscribed soil particle. Based on the minimum particle diameter of the selected PSD, the side lengths of each cube is divided into a finite number of subdivisions (Figure 9a). The n, the packing algorithm tracks the total number of subdivisions occupied by each incoming box, i.e. each packed soil particle, as packing proceeds based on the MonteCarlo simulation corresponding to a given PSD. Finally, the PAGE 33 19 automated algorithm fills the maximum possi ble number of subdivisions until it finds that no space is available within the consider ed cube for further packing of soil particles from the given PSD. On the other hand, in order to obtain the minimum void ratio (densest packing), the maximum number of spheres smaller than th e inscribing sphere in each box is also packed into the unoccupied space at the co rners of each box (Figure 9c and 9d) before that box is placed in the cube (Figure 9b). Fo r each type of soil, 400 such cubes (trials) were packed randomly using the MonteCarlo simulation. Since the maximum and minimum void ratios obtained in each trial w ould be different as explained above, they are considered as random variables. T hus, the probability distributions of emax and emin obtained from those 400 trials are show n in Figure 10. The ranges derived for emax and emin as shown in Table 1 agree with the typical values in [15]. Table 1: Size Characterist ics of Pavement Layers Soil type D15 (mm) D50 (mm) D85 (mm) emin emax Gravel 9.85 13 19 0.59 0.8 0.98 1.18 Coarse sand 1.85 2.5 3.5 0.73 0.93 0.91 1.10 Fine sand 0.16 0.29 0.49 0.54 0.79 0.911.26 PAGE 34 20 Figure 9: Random Packing of Soil Particles 2.2.1.3 Determination of the Natural Void Ratio Distribution The concept of relative density is helpfu l in quantifying the level of compaction of coarsegrained soils. The relative density of a coarsegrained soil at a given compaction level expresses the ratio of th e reduction in the voids at the given compaction level, to the maximum possible decrease in the voids (Eqn 3). The insitu void ratio di stribution (e) corresponding to a relative density (Dr) is obtained using the previously obtained emax and emin (Eqns 3 or 4). % 100min max max e e e e Dr (4) PAGE 35 21 or min max100 100 1 e D e D er r (5) Where, emax is the randomly distributed void ratio of gravel/sand in its loosest state. emin is the randomly distributed void ratio of gravel/s and in its densest stat e. e is the randomly distributed void ratio of gravel/sand in its natura l state in the field. Thus, r e e enD e f e f F e f ,min max (6) Wheree fen, e fe maxand e fe minare Probability Density Functions (pdf) of insitu natural void ratio, maximum void ratio and minimum void ratio respectively. If e is a function of two random variables (emax = e1, emin = e2) and 1 and 2 are the mean values of these random variables, the expected mean value of e can be expressed using the secondorder Taylor se ries approximation as 2 1 2 1 2 1 2 1 2 2 1, 2 1 e e Cov e e e e e Eii (7) Furthermore, the variance of e can be expressed as 2 1 2 2 1 2 1 1,2 1e e Cov e e e e e Vii (8) PAGE 36 22 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.550.650.750.850.951.051.15void ratioProbability density or frequenc y maximum void ratio for gravel minimum void ratio for gravel Figure 10: Probability Density Function (PDF) of Maximum and Minimum Void Ratios Two alternative approaches can be followed to determine the probability distribution of the natural void ratio (e) from the probabil ity distributions of the maximum and minimum void ratios (Eqn 5): 1) Assume an appropriate Probability Density Function (Example: Lognormal) with mean and variance calculated using Eqn 6 and Eqn 7 from the means (1, 2) and variances (V (e1), V (e2)) of emax and emin. 2) Generating the probability distribution of natural void ratios (e) using Eqn 4 and the MonteCarlo simulation technique. In this model, the method (2) listed above is used. When the probability distributions of emax and emin are known for each soil type (Figure 10) the corresponding cumulative distributions of emax and emin can be derived (Figure 11). Then, once more the Monte 1 2 V (e1) V (e2) PAGE 37 23 Carlo simulation (Section 3.1) can be used to pick two arrays of emax and emin values for that soil type. Finally, with the emax and emin arrays, a corresponding array of natural void ratio (e) of any of the above soil t ypes, for a given relative density (Dr) can be obtained from Eqn 4. The spatial distributions of na tural void ratios (porosities) so obtained assumed to be representative of each soil layer will be used in the analysis of flow (Section 4). Figure 11: Cumulative Frequencies of Maximum and Minimum Void Ratios 2.3 Modeling of Steady State Flow In modeling the pavement layers, the three di mensional porosities obtained from particle packing are coupled with two dimensional water flow. This is because flow is constrained in the third direction due to the two dimensional nature of the pavement. Due to the incompatibility of the sizes of the particles of the three layers a nd since the same grid element size is used for entire pavement structure, only two pavement layers are modeled at a time. Navier Stokes e quations [13] are given by: Random selection PAGE 38 24 Mass Conservation (Continuity Equation): 0 y v n x u n (9) Momentum Conservation (Momentum Equations): X direction: 2 2 2 2 2y u n x u n D x p n y v u n x u n t u nx (10) Y direction: y yg n y v n x v n D y p n y v n x v u n t v n 2 2 2 2 2 (11) Where, n porosity at the locati on (x, y) at time t, u, v fluid velocities in the x and y directions respectively, fluid density, p fluid pressure, fluid viscosity, gy gravitational force per unit mass in the y direction. Averaged fluidparticle inte ractions (drag forces) are qu antified using semiempirical relationships provided in Eqn 9 [6]. 2 2 2 2 21 75 1 1 150 u d n n u d n n Dp p x (12) Where, pd averaged particle diameter; Dx x directional average fluidparticle iteration force per unit volume. A similar expression for Dy is used for y directional averaged fluidparticle interactions. PAGE 39 25 2.3.1 Need for a Staggered Grid The central difference form of Eqns 6 8 requi res a staggered grid due to the two issues explained below: 0 2 21 1 , 1 1 y v v x u uj i j i j i j i (13) Figure 12: Discrete Checkerboard Velocity Distribution at Each Grid Point: At Each Node, the Top Number is u and the Bottom Number is v Figure 12 illustrates an arbitrar y zigzag type of distributi on of both the x component and y components of velocity, u and v, respect ively. When these arbitrary numbers are substituted in Eqn 12, the central differen ce form of the continuity equation for incompressible fluid flow is unconditionally sa tisfied. In other word s, the checkerboard velocity distribution shown in Figure 12 produces a trivial solution inapplicable to a real physical flow field. PAGE 40 26 Considering a twodimensional discrete, checkerboard pressure pattern as illustrated in Figure 13, the second order central difference formulation for the pressure gradients which appear in the momentum equa tions can be written as follows: x P P x Pj i j i 2, 1 1 (14) y P P y Pj i j i 21 1 (15) Figure 13: Discrete Checkerb oard Pressure Distribution The checkerboard pressure distribution (Figure 13) gives zero pressure gradients in the momentum equations in the x and y directions written in terms of the central difference scheme (Eqns 13 and 14 respectively). In othe r words, the pressure field would not be PAGE 41 27 incorporated in the discretized Navier Stoke s equations and hence the numerical solution would effectively see only uniform pressure distributions in the x and y directions. In order to address the above issue, an innovative remedy for the checkerboard distribution is to use a staggered mesh wh ere discrete pressures and velocities are expressed only wherever required. A t ypical twodimensional staggered mesh arrangement is shown in Figure 14. In the st aggered grid, the velo city components are computed for the points that lie on the faces of the flow elements. Thus, the xdirectional velocity, u, is computed on the planes or su rfaces that are normal to the xdirection and similarly, the ydirectional velocity, v, is computed on the planes or surfaces that are normal to the ydirection. On the other hand, th e pressures are computed at the center of the flow elements. By introducing a staggered gri d, the mass flow rates across th e flow element faces can be evaluated without any interpolation of the re levant velocity components. Moreover, for a typical flow element, it will be easy to see that the discretized continuity equation would contain the differences of adjacent velocity components thus preventing a wavy velocity field resulting from the con tinuity equation. When the st aggered grid is used, only realistic velocity fields would have the pos sibility of being acceptable to the continuity equation. Consequently, pressure fields in the momentum equations would no longer be felt as uniform pressure fields. PAGE 42 28 Figure 14: A Typical Staggered Grid Arrangement (Solid Circles Represent Pressure Nodes and Open Circles Represent Velocity Nodes) 2.3.2 Numerical Solution Technique Based on the staggered grid arra ngement introduced in Section 4.2, a finitedifference approach is used to discretize the Eqns 9, 10 and 11. The scheme is based on forward difference in time and central difference in spac e. The numerical form of the x directional momentum equation is written in Eqn 10, refe rring to the notation in Figure 15 for the sequential iteration steps of N and N+1. Y X PAGE 43 29 Figure 15: Computational Module for the XMomentum Equation 2 j 2 1 i j i p 2 j i j i j 2 1 i 2 j i p 2 j i 2 j i j i j 1 i j i 2 1 j 2 1 i 1 j 2 1 i 1 j i 2 1 j 2 1 i 1 j 2 1 i 1 j i 2 j, 2 1 i j 1 i 2 j 2 3 i j 1 i N j 2 1 i j i 1 N j 2 1 i j iu d n n 1 75 1 u d n n 1 150 x p p n y 2 v u n v u n x 2 u n u n t u n u n 2 j 2 1 i j i 1 j 2 1 i 1 j i 1 j 2 1 i 1 j i 2 j 2 1 i j i j 2 1 i j 1 i j 2 3 i j 1 iy u n 2 u n u n x u n 2 u n u n (16a) Where, 2 1 j 1 i 2 1 j i 2 1 j 2 1 iv v 2 1 v a int po At (16b) PAGE 44 30 2 1 j 1 i 2 1 j i 2 1 j 2 1 iv v 2 1 v b int po At (16c) Within the staggered grid some velocities n eed to be interpolated as shown in Eqns 16b and 16c. In a concise form, Eqn 16a can be written as t A x p p ) t ( A t A u n u nj i j 1 i N j 2 1 i j i 1 N j 2 1 i j i (16d) Where, x 2 u n u n A2 j 2 1 i j 1 i 2 j 2 3 i j 1 i y 2 v u n v u n2 1 j 2 1 i 1 j 2 1 i 1 j i 2 1 j 2 1 i 1 j 2 1 i 1 j i (16e) j in A (16f) and 2 j 2 1 i j i p 2 j i j i j 2 1 i 2 j i p 2 j i 2 j iu d n n 1 75 1 u d n n 1 150 A 2 j 2 1 i j i 1 j 2 1 i 1 j i 1 j 2 1 i 1 j i 2 j 2 1 i j i j 2 1 i j 1 i j 2 3 i j 1 iy u n 2 u n u n x u n 2 u n u n (16g) The numerical form of the y directional mo mentum equation can be written similarly based on Eqn 11. PAGE 45 31 2.3.3 SemiImplicit Method for PressureLi nked Equations (SIMPLE) Algorithm An iterative process called the pressure correction technique has found widespread application in the numerical solution of the incompressi ble, viscous Navier Stokes equations. This technique is a vehicle by wh ich the velocity and pressure fields are directed towards a solution that satisfies both the discrete cont inuity and momentum equations. This technique is embodied in an algorithm called SemiImplicit Method for PressureLinked Equations (S IMPLE) [11]. The primary idea behind SIMPLE is to create a discrete equation for pressure in terms of the pressure correction, from the discrete continuity equation. Figure 16 shows th e stepbystep procedure for the SIMPLE algorithm. PAGE 46 32 Figure 16: Flow Chart for SIMPLE Algor ithm for Steady State Condition 2.3.3.1 Derivation of the Pressure Correction Formula At the beginning of each new iteration, j ip, is set to* j ip, where j ip is the pressure from the previous iteration. Then, Eqn 16d will be transformed to Eqn 17. t A x p p ) t ( A t A u n u n* j i j 1 i * N j 2 1 i j i 1 N j 2 1 i j i (17) Similarly rearranged y directi onal equation would be Eqn 18. y j i j i 1 j i * N 2 1 j i j i 1 N 2 1 j i j ig n t B y p p t B t B v n v n (18) Iteration for steadystate Assume trial values of (p*) n, (u*)n, (v*) n for nth iteration Solve momentum equations for (n+1)th iteration [(u*) n+1 and (v*) n+1] Solve a pressurecorrection equation for p' using a relaxation technique Calculate pn+1 = (p*) n + p' Correct the velocities u n+1 = (u*) n + u' v n+1 = ( v* ) n+ v' Visualize steadystate results Determine the mass imbalance term (mit) mit < tol. Yes N o PAGE 47 33 Next, Eqns 17 and 18 are subtracted from E qn 16d and its y direc tional counterpart to obtain Eqns 19 and 20 respectively. x p p ) t ( A u n u nj i j 1 i N j 2 1 i j i 1 N j 2 1 i j i (19) y p p t B v n v nj i 1 j i N 2 1 j i j i 1 N 2 1 j i j i (20) Where, the pressure correction , j i j i j ip p p ,* 'A A A 'A A A j in A A A and. B B ,' and B can be defined similarly from the y directional coefficients. Eqn 9 can be rewritten in the discretized fo rm as in Eqn 21 using the central difference scheme with the staggered grids, (21) The semiimplicit terminology refe rs to arbitrary setting of A A' 'and the corresponding y directional coefficients equa l to zero, thus allowing the ultimate pressure correction formula, (Eqn 21), to have p appearing only at four nei ghborhood grid points illustrated in Figure 17 Omission of N j iu, 2 1', N j iv2 1 ,', A A' 'and the correspondin g y directional coefficients in the pressure correction deriva tion does not affect the final solution because 0 y v n v n x u n u n2 1 j i j i 2 1 j i 1 j i j 2 1 i j i j 2 1 i j 1 i PAGE 48 34 the pressure correction and velocity corrections would in fact be zero for a converged solution. Arbitrarily setting N j iu, 2 1' and N j iv2 1 ,' equal to zero in Eqns 19 and 20 produce, x p p ) t ( A u n u n u nj i j 1 i 1 N j 2 1 i j i 1 N j 2 1 i j i 1 N j 2 1 i j i (22) y p p t B v n v n v nj i 1 j i 1 N 2 1 j i j i 1 N 2 1 j i j i 1 N 2 1 j i j i (23) Using Eqns 22 and 23, the discretized continu ity equation (Eqn 21) can be expressed as a pressure correction formula in terms of j ipvalues of only the neighboring four grid points as illustrated in Figure 17 (Eqn 21). Figure 17: Designation of Nodal Points on a Grid Used in SIMPLE Algorithm 0 f p e p d p c p b p aj i 1 j i j i 1 j i j i j 1 i j i j 1 i j i j i j i (24a) Where, 2 2 j i j iy t x t n 2 a (24b) 2 j i j ix t n b (24c) PAGE 49 35 2 j i j ix t n c (24d) 2 j i j iy t n d (24e) 2 j i j iy t n e (24f) and y v n v n x u n u n f1 N j 1 i j i 1 N j 1 i j i 1 N 1 j i j i 1 N 1 j i j i j i (24g) After determining the pressure correction,' j ip from Eqn 24, the pressure and velocity components at every node are updated as follows: j i j i j ip p p (25) j i j i j iu u u (26) j i j i j iv v v (27) Where, x p p n ) t ( A uj i j 1 i j i j 2 1 i (28) and y p p n t B vj i 1 j i j i 2 1 j i (29) The pressure correction Equation (Eqn 24), whic h is of Poisson format in terms of p can be solved by employing a numerical relaxation technique. The term j if,(Eqn 24g) is called the mass imbalance term which must va nish theoretically in the last iteration PAGE 50 36 where the velocity field converges to a field th at satisfies the continuity equation (Eqn 9). In the numerical algorithm developed in this work, the mass imbalance term is used as a stopping criterion to assure that the solution converges to the correct velocity field. The function of the pressure corre ction formula is to set the iterative process in such a direction that, when the velocity distri bution is determined from the momentum equations, it will eventually converge to the correct distribution which satisfies the continuity equation. Because the pressure co rrection method is designed to solve for the steady flow condition via an iterative process, the superscripts N and N+1 used in the equations are the sequential iteration steps, with no significance to any real transient variation. Under the steady st ate conditions, the term of t can be treated as a parameter which has some effect on the speed at which the convergence is achieved. 2.3.3.2 Boundary Conditions for Pavement System The boundary conditions appropriate for the actual water flow through the pavement system are specified in the computer code as indicated in Figure 18. 1) At the inflow boundary, the pressure and velocities are specified. Hence, the pressure correction, p is zero at the inflow boundary. 2) At the outflow boundary, the pressure is specified and veloc ity components are allowed to float. Hence, p is zero at the outfl ow boundary as well. 3) At the vertical walls, the slip condit ions are maintained. Thus, the velocity component, u, normal to the walls is set to zero. PAGE 51 37 Figure 18: Boundary Conditions Incorporated in the Flow through the Pavement System 2.4 Validation of the Numerical Model In order to validate the flow model described in Section 4, the coe fficients of hydraulic conductivity of several unif ormly graded soil types which are not used in the current illustration are determined from the numerical model and compared with widely used empirical relationships proposed by Hazen (Eqn 30) and Chapius (Eqn 31) [16]. The results of this comparison are summarized in Table 2. 2 10D c ) s / cm ( k (30) 7825 0 e 1 e D 4622 2 ) s / cm ( k3 2 10 (31) Where, c a constant that varies from 1.0 to 1.5 and D10 the effective size in mm. PAGE 52 38 Table 2: Comparison of 3D Coefficients of Hydraulic Conductivity Derived from SIMPLE Algorithm Figure 19: Relationship between Hydr aulic Conductivity Vs Diameter for Uniformly Distributed Soils It can be concluded from Table 2 that th e coefficients of hydr aulic conductivities determined from the model agree fairly well with those computed from common empirical relationships. Furthermore, it was de termined that the coefficients of hydraulic conductivity for gravel, coarse sand and fine sand used in the Numerical Illustration Soil type (1) Effective size (D10) mm (2) K model (cm/s) (3) KHazen (cm/s) (4) KChapius (cm/s) (5) a 0.5 0.979871 0.375 0.497428 b 1 3.108229 1.5 1.47206 c 2 6.811264 6 4.532172 d 3 9.374428 9 8.247057 e 4 12.34681 16 12.94549 PAGE 53 39 (Section 7) are 0.097 m/s ,0.043 m/s and 0. 00024 m/s respectively. These values are within the ranges of typical hydraulic conductivities presented in the literature [16] for the corresponding PSDs shown in Figure 8. According to the developed model, the relationship between hydraulic c onductivity Vs square of the te nth percentile diameter for uniformly distributed soils is shown in Figur e 19. Therefore, the empirical equation (Eqn 45) does not seem to be accurate. 2.5 Numerical Illustration The numerical model was applied to the pave ment layer interfaces shown in Figure 2. The material properties and other model para meters shown in Table 1 and 3 were used for this purpose. PAGE 54 40 Table 3: Model Parameters Particles Number of gravel particles considered for determining the porosity distribution 137668 Number of coarse sand particle s considered for determining the porosity distribution 205332 Number of fine sand particle s considered for determining the porosity distribution 855450 Gravel particle size 2 mm24 mm Coarse sand particle size 1 mm5mm Fine sand particle si ze 0.04 mm 0.62 mm Average density of saturated soil 1900 kg/m3 Average compression index for granular soils 0.37 [18] Specific gravity of soil 2.65 Average void ratio of gravel packing (Dr = 50%) 0.86 Average void ratio of coarse sand packing (Dr = 50%) 0.90 Average critical hydraulic gradient 0.88 Water Density 1000 Kg/m3 Dynamic viscosity 103 Pa.s Flow calculations Gravitational acceleration 9.81 m/s2 Number of grid elements 50 by 50 PAGE 55 41 2.5.1 Results of Numerical Illustration Figure 20 and Figure 21 illustrate the progression of flow rate with iterations through different soil layers for different levels of compaction. They clearly reveal the reduction in flow rate due to densification a nd the impact of th e hydraulic conductivity. Coarse sandgravel interface0.E+00 2.E05 4.E05 6.E05 8.E05 1.E04 1.E04 050100150200250300 iteration stepsFlow Rate (m3/s) Dr = 5% Dr=50% Dr = 75% Dr = 95% Figure 20: Plot of Flow Rate through C oarse SandGravel Interface Vs Iteration Steps with Different Levels of Compaction (For Pressure Differential of 23 kPa) 5.0E08 0.0E+00 5.0E08 1.0E07 1.5E07 2.0E07 2.5E07 3.0E07 3.5E07 4.0E07 4.5E07 050100150200250Flow Rate(m3/s)Iteration stepsFine sand/coarse sand interface Dr = 5% Dr = 50% Dr = 75% Dr = 95% Figure 21: Plot of Flow Rate through Fine SandCoarse Sand Interface Vs Iteration Steps with Different Levels of Compaction (For Pressure Differential of 23 kPa) PAGE 56 42 Figure 22 shows that at the steady state the local hydraulic gradient s at the interface are larger than that at other locations. They also clearly exceed the values predicted by conventional analysis (Darcys law). 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.40.50.60.70.80.911.1 Local hydraulic gradient / icrDepth of layer/ m Localized i Particulate i across subgradeD'Arcy i across subbaseD'Arcy Figure 22: Comparison of Hydraulic Gr adients Obtained Using Continuum and Particulate Approaches (Using SI MPLE Algorithm Steady State) It is seen from Figure 23 that when the hydr aulic gradient is 1.7 times the critical hydraulic gradient, the water flow never reaches a steady state. Coarse sand / gravel interface with Dr = 95%0.0E+00 1.0E04 2.0E04 01002003004005006007008009001000iteration stepFlow Rate/ (m3/s) Figure 23: Plot of Flow Rate Vs Iteration Step for High Hydraulic Gradient (i = 1.7 icr) Coarse sand gravel interface PAGE 57 43 2.6 Conclusions The current practice of designing paveme nt filter systems does not consider 1) The interaction between pore water a nd the individual soil particles. 2) The localized effects. In order to determine the realistic limiting hyd raulic gradients that can be applied within pavement layers, a design me thodology based on a particulate approach that incorporates particlesoil interaction needs to be used. In this Chapter, the steady state water seepage in two saturated filter interfaces with varyi ng levels of compaction was analyzed using a soil particulate model. The part iculate effects of soil with different levels of compaction were incorporated convenient ly in the model using a random packing technique, while the flow of water within the particulate assembly was modeled by the Navier Stokes flow equations. Two separate filter interfaces, i.e. coarse sandgravel and fine sandcoarse sand were assembled using partic le size distributions that satis fied the conventional filter design criteria. Then, a pressure differentia l that corresponded to the critical hydraulic gradient was applied across the layer interf ace. In order to verify the model, the coefficients of hydraulic conduc tivity predicted from the mode l are compared with those computed from the empirical relationships wi dely used in drainage design. The steady state flow algorithm is modified in the next Chapter to analyze the transient behavior of seepage through twolayer par ticulate pavement system wh ich predicts the critical conditions for erosion, piping and clogging. PAGE 58 44 CHAPTER THREE ANALYSIS OF TRANSIENT WATER SEEPAGE IN A PAVEMENT SYSTEM USING THE PARTICULATE APPROACH 3.1 Modeling of Transient Flow The comprehensive analytical procedure a nd the computer code developed for its implementation are illustrated in Figure 24. The analytical procedure consists of two primary tasks such as the random assembly of the particulate medium (granular soil) and the solution of the fluid flow governing equations using par tial coupling between the two media. The flow chart also includes the s ections, equation numbers and figure numbers corresponding to each stage. PAGE 59 45 Figure 24: Flow Chart Illustrating the An alytical Model for Transient and Steady State Flow No Yes No Yes No Packing to determine maximum and minimum void ratios for each particle size distribution (Section 2.2.1) Application of MonteCarlo simulation to obtain randomly distributed insitu void ratios at specified relative densities (Section 2.2.1.3) Assignment of random insitu void ratios for grid elements used in the fluid flow analysis Application of NavierStokes equations to each grid element (Section 3.1.1) Numerical discretization based on a staggered grid (Section 3.1.3 & Figure 25) Analysis of results critical velocities and hydraulic gradients during transient and steadystate conditions (Section 3.4) Solid phase Fluid phase Iterative solution for water pressure (p) distribution in the grid (Eqn. 44 & Figure 26) Time marching p < Tol Timebased change in water velocity ( v) < Tol Validating the model by determining the coefficients of hydraulic conductivities (Section 3.3) Yes mit < tol. Determination of porosity change satisfying the volumetric compatibility (Eqn. 43). Determine the mass imbalance term (mit) PAGE 60 46 3.1.1 Transient Navier Stokes Equations In modeling the pavement layers, the three di mensional porosities obtained from particle packing are coupled with two dimensional water flow. This is because flow is constrained in the third direction due to the two dimensi onal nature of the pavement. Navier Stokes equations [13] are given by: Mass Conservation (Continuity Equation): 0 y v n x u n t n (32) Momentum Conservation (Momentum Equations): X direction: 2 2 2 2 2y u n x u n D x p n y v u n x u n t u nx (33) Y direction: y yg n y v n x v n D y p n y v n x v u n t v n 2 2 2 2 2 (34) Where, n porosity at the locati on (x, y) at time t, u, v fluid velocities in the x and y directions respectively, fluid density, p fluid pressure, fluid viscosity, gy gravitational force per unit mass in the y direction. Averaged fluidparticle inte ractions (drag forces) are qu antified using semiempirical relationships provided in Eqn 35 [6]. 2 2 2 2 21 75 1 1 150 u d n n u d n n Dp p x (35) PAGE 61 47 Where, pd averaged particle diameter; Dx x directional average fluidparticle iteration force per unit volume. A similar expression for Dy is used for y directional averaged fluidparticle interactions. 3.1.2 Volumetric Compatibility of Solid and Fluid Phases The effective stress in the soil phase can be expressed using Eqn 36 [16]. u (36) Where, effective vertical stress, total vertical stre ss and u pore water pressure After each time step, when the effective stre ss changes due to the change in pore water pressure, the void ratio must change according to the compressibility characteristics of the soil. The typical loglinear void ratio versus effective stre ss relation [16] that characterizes saturated finegrained soil (clay, silt) is assumed for coarsegrained soil (gravel, sand) as well (Eqn 37). ) log( c oC e e (37a) d C dec4343 0 (37b) Where, Cc is the equivalent compression index of the soil. The corre sponding void ratio change can be obtained usi ng Equations (38a) and (38b). e 1 e n (38a) 2n 1 de dn (38b) PAGE 62 48 3.1.3 Numerical Solution Technique A finitedifference approach [12] is used to discretize the Eqns 32, 33 and 34. The scheme is based on forward difference in ti me and porosity, and central difference in space. A typical twodimensional staggered mesh arrangement [11, 13] was used to discretize the Eqns 32, 33 and 34. On the othe r hand, Patankars original formulation [11] is based on the pressure co rrection method involving a finitevolume approach. The pavement layers were divided into a 50 50 grid for the finite difference analysis. One fluid cell, known herein as a grid elemen t, was selected to be 25 mm 25 mm 1m (2D flow) to be compatible with the c ube sizes (28 mm 28 mm 28 mm or 16 mm 16 mm 16 mm) used for packing (Section 3. 1). 1 m has been selected in the third direction assuming two dimensional water flow Since the Navier Stokes equations are written in terms of porosities, the spatial proba bility distribution of natural void ratios are first converted to porosities (E qn 41a). First, natural porositie s (void ratios) were assigned to each grid element based on the spatial probability distribution of natural void ratios determined from packing (Section 3.3). With in each grid element, the fluidparticle interaction is quantified consid ering the averaged particle di ameter which is defined as the arithmetic mean of the di ameters of all the particles in each grid. The authors modified numerical scheme is described as follows: Eqn 39a expresses the x direc tional momentum equation in the numerical form for the sequential time steps of N and N+1 re ferring to the notation in Figure 25. PAGE 63 49 Figure 25: Computational Grid for the XMomentum Equation 2 j 2 1 i j i p 2 j i j i j 2 1 i 2 j i p 2 j i 2 j i j i j 1 i j i 2 1 j 2 1 i 1 j 2 1 i 1 j i 2 1 j 2 1 i 1 j 2 1 i 1 j i2 j 2 1 i j 2 1 i 2 j 2 3 i j 2 3 i N j 2 1 i j 2 1 i 1 N j 2 1 i j 2 1 iu d n n 1 75 1 u d n n 1 150 x p p n y 2 v u n v u n x 2 u n u n t u n u n 2 j 2 1 i j i 1 j 2 1 i 1 j i 1 j 2 1 i 1 j i 2 j 2 1 i j i j 2 1 i j 1 i j 2 3 i j 1 iy u n 2 u n u n x u n 2 u n u n (39a) Where, At point a, 2 1 j 1 i 2 1 j i 2 1 j 2 1 iv v 2 1 v (39b) PAGE 64 50 At point b, 2 1 j 1 i 2 1 j i 2 1 j 2 1 iv v 2 1 v (39c) Similarly, the y directional momentum equation can be written numerically. Equations (37b) and (38b) can be used to express the change in porosity within the time interval t using the volumetric compatibility between solid and fluid phases. Then, Eqn 32 can also be discretized as follows: t n y v n v n x u n u n1 N 2 1 j i 2 1 j i 1 N 2 1 j i 2 1 j i 1 N j 2 1 i j 2 1 i 1 N j 2 1 i j 2 1 i t 1 p p p n 1 C 4343 0N j i N j i 1 N j i N j i 2 j i c (40) To compute the change in porosity in the tr ansient continuity e quation in an explicit manner, the water pressure difference between the time steps (N1) and N has been used in Eqn 40. By substituting ve locities of the new time step by those of the previous time step, Eqn 40 is reduced to Eqn 41, which is a Poisson equation in terms of water pressure (p). Figure 26: Designation of Nodal Points on a Grid Used in the Iterative Procedure 0, 1 , 1 , 1 , 1 , j i j i j i j i j i j i j i j i j i j i j iS p e p d p c p b p a (41a) Where, PAGE 65 51 N j i j i j i c j i j i j i j i j ip n C t y n n x n n a. 2 2 2 2 1 2 1 2 2 1 2 1 ,1 4343 0 (41b) 2 2 2 1 ,t x n bj i j i (41c) 2 2 2 1 ,t x n cj i j i (41d) 2 2 2 1 ,t y n dj i j i (41e) 2 2 2 1 ,t y n ej i j i (41f) and 1 , 2 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 2 1 2 1 2 2 1 2 1 2 2 1 2 1 2 2 1 2 1 ,1 4343 0 N j i N j i j i j i c N j i j i N j i j i N j i j i N j i j i N j i N j i N j i N j i N j i N j i N j i N j i j ip p n C t y v n v n t x u n u n y t B B y t B B x t A A x t A A S y t g n ny j i j i 2 2 1 2 1 (41g) Where, A A B and B is a collection of velocity de rivatives from the discretized momentum equations in both x and y directions. PAGE 66 52 Similar to the pressure correction method introduced by Patankar [11], the authors modified algorithm solves the Poisson equation (Eqn 41) iteratively in terms of water pressures (p) at the neighboring nodes shown in Figure 26. Within each time step, when the water pressure difference at any node comp uted in each iterati on becomes negligible, the algorithm checks the volum etric compatibility of solid and fluid phases throughout the entire pavement system (Eqn 40). Until the velocities become constant at each node, the time marching is continued while solving for water pressures iter atively at each time step (Figure 24). 3.1.4 Boundary Conditions The boundary conditions appropriate for the actual water flow through the pavement system are specified in the computer code as indicated in Figure 27. 1) At the inflow boundary, the pressure and velocities are specified. 2) At the outflow boundary, the pressure is specified and veloc ity components are allowed to float. 3) At the vertical walls, the slip condit ions are maintained. Thus, the velocity component, u, normal to th e walls is set to zero. PAGE 67 53 Figure 27: Boundary Conditions Incorporated in the Flow through the Pavement System (Saturated Soil) 3.2 Interface Effect in FDM Since the hydraulic properties are different in different layers, Eqn 44 must be solved separately within each layer ex cluding the interface grid points. On the other hand, at the interface, Eqn 44h is used to satisfy the continui ty of flow from one layer to the other. 2 2 1 2 1 2 1 2 1 2 2 1 2 1 2 1 2 1 1 2 1 1 2 1 int 2 1 2 1 j i j i j i j i j i j i j i j i j i j i j i j i erface j i j i j iv k v k v k v k p n p n p n n (41h) Where, y d n n kj ip j i j i j i 2 2 2 1 2 2 1 2 1 ,2 1 ,1 150 and y d n n kj i j ip j i j i 2 1 2 1 ,2 2 2 1 2 1 ,1 75 1 Saturated soil PAGE 68 54 3.3 Validation of the Numerical Models In order to validate the flow mo del described in Section 3.2, the coefficients of hydraulic conductivity of several unif ormly graded soil types which are not used in the current illustration are determined from the numerical model and compared with widely used empirical relationships proposed by Hazen (Eqn 42) and Chapius (Eqn 43) [16]. The results of this comparison are summarized in Table 4. 2 10D c ) s / cm ( k (42) 7825 0 e 1 e D 4622 2 ) s / cm ( k3 2 10 (43) Where, c a constant that varies from 1.0 to 1.5 and D10 the effective size in mm. Table 4: Comparison of 3D Coefficients of Hydraulic Conductivity Derived from the Numerical Methodology Soil type Effective size (D10) mm K from steady flow model using SIMPLE algorithm (Chapter 2) (cm/s) Kfrom current model (cm/s) KHazen (cm/s) KChapius (cm/s) a 0.5 0.979871 1.01984 0.375 0.497428 b 1 3.108229 3.230112 1.5 1.47206 c 2 6.811264 6.766702 6 4.532172 d 3 9.374428 9.603897 9 8.247057 e 4 12.34681 11.62177 16 12.94549 PAGE 69 55 In order to solve the transi ent continuity equation and Navier Stokes equations, the authors have developed the model presented in this study considering the volumetric compatibility between solid and fluid phases. As opposed to Patankars work [11], there is no pressure correction method involved in the authors modification. Instead, an elliptic equation in terms of the water pressu re is iteratively solved until the volumetric imbalance in any grid element becomes ne gligible. The coefficients of hydraulic conductivity determined using th e current model are labeled as K from current model in Table 4. As an alternative, the authors have also analyzed flow thr ough the same pavement system solving timeindependent (steady state) cont inuity equation and Navier Stokes equations using a finitedifference approach similar to the pressure correction method in the SIMPLE (Semi Implicit Method for PressureL inked Equations) algo rithm [11]. For the purpose of comparing the solutions obtained fr om the current model for the special case of steady state flow, the K value was also co mputed from steady state flow analysis in Chapter 2. In Table 4, these values are la beled as K from steady flow model using SIMPLE algorithm. It can be concluded from Table 4 that th e coefficients of hydraulic conductivities determined from the current model agree fair ly well with those computed from common empirical relationships. It also agrees with the K values predicted from a separate numerical model which the authors have de veloped based on the SIMPLE algorithm in Chapter 2 to solve just the steady fluid flow problems. Furthermore, it was determined that the coefficients of hydraulic conductiv ity for gravel and coarse sand used in the PAGE 70 56 Numerical Illustration (Section 3.4) are 0. 097 m/s and 0.043 m/s respectively. These values are within the ranges of typical hydraulic conductivities presented in the literature [16] for the corresponding PSDs shown in Figure 8. 3.4 Numerical Illustration The numerical model was applied to the pave ment layer interface shown in Figure (2a). The material properties and other model para meters shown in Table 1 and 3 were used for this purpose. 3.5 Selection of Time Step Eqn 41a can be rewritten in the matrix form as S p A (44) Where A will be a banded matrix with elements that are formed from ai,j bi,j ci,j di,j and ei,j evaluated at nodes where (pi,j ) are unknown nodal pressures. p is the vector of unknown pressures. Meanwhile, the S vector contains Si,j terms evaluated at all the nodes and the known (fix ed) nodal pressures (pi,j) at boundaries. Extremely small time steps can make the coefficient matrix [A] become illconditioned which could result in numerical instability in the solu tion procedure. On the other hand, if the time step is too large the numerical methods would not provide precise solutions due to truncation errors associated with the difference expression for the time derivative [13]. Hence, an optimum time step of 1e4 seconds is used in this paper. PAGE 71 57 3.6 Results of Numerical Illustration Figures 28 and 29 respectively illustrate how the pore water pressures and water velocities are developed w ithin the layers in real time upon a sudden pressure buildup at the bottom of the coarse sandgravel layer sy stem. The transient effects are clearly seen from the isochrones shown in Figures 28 and 29. It is also seen that it takes 2.5e2 seconds for the steady state flow to establish. Coarse sandgravel layers Dr = 95% iapplied = icr = 0.880 0.25 0.5 0.75 1 1.25 1.5 50000500010000150002000025000Water pressure / PaDepth of pavement/ m t=1e4 sec t = 4e4 sec t = 2.5e3 sec Steady condition t = 2.5e2 sec Figure 28: Development of Water Pressure within the Coarse SandGravel Layers with Time (For Pressure Differential of 23 kPa) PAGE 72 58 0 0.25 0.5 0.75 1 1.25 0.0100.010.020.030.04Depth of pavement/mWater velocity (m/s)Coarse sandgravel layers Dr= 95% iapplied= icr= 0.88 t=1e4 sec t = 2e4 sec t = 4e4 sec Figure 29: Development of Water Velocity within the Coarse SandGravel Layers with Time (For Pressure Differential of 23 kPa) Figure 30 illustrates the progre ssion of flow rate with tim e for different levels of compaction. They clearly reveal the reduction in flow rate due to densification and the impact of compaction on th e hydraulic conductivity. Coarse sandg ravel la y ers iapplied = icr = 0.882.0E05 2.0E05 6.0E05 1.0E04 1.4E04 02004006008001000 Time steps ( t = 1e4 sec)Flow Rate/ Area (m/s) Dr = 95% Dr = 75% Dr = 50% Dr = 5% Figure 30a: Plot of Flow Rate through C oarse SandGravel Layers Vs Time Steps with Different Levels of Compaction (F or Pressure Differen tial of 23 kPa) PAGE 73 59 Figure 30b: Plot of Flow Rate through Fine SandCoarse Sand Layers Vs Time Steps with Different Levels of Compaction (For Pressure Differential of 23 kPa) Determination of the localized maximum hydr aulic gradients and their locations would be helpful in the design of pavements. Figu re 31a was produced to depict the highest value of imax/icr anywhere within the coarse sandgra vel layer system. During initial time steps, the high imax/icr appears close to the bottom boundar y where a high water pressure of 23 kPa is applied suddenly to initiate flow. This is natura l, because in reality, sudden application of high pressure cr eates instability at the point of application. With time, the magnitude of imax/icr reduces and the high imax/icr gradually moves closer to the interface (Figure 31). As illustrated in Figure 31a, wh en the steady state is reached the spatial maximum hydraulic gradient is nearly equal to the average hydraulic gradient applied across the pavement layers i.e. the critic al hydraulic gradient. Moreover, the overall (temporal) maximum hydraulic gr adient that occurs after th e initial condition effects on PAGE 74 60 bottom boundary is seen close to the interface during transient condition, varying slightly with the compaction level. From Figure 31b, it can be seen that th e impact of compaction on the maximum hydraulic gradient is only ma rginal even under tr ansient conditions. 1 6 11 16 21 26 31 020406080100imax/ icrTime stepsCoarse sandgravel layers iapplied= icr = 0.88 Dr = 95% Dr = 75% Dr = 50% Dr = 5% Figure 31a: Plot of Maximum Hydraulic Gr adient within the Coarse SandGravel Layers Vs Time Steps for Different Levels of Compaction (For Pressure Differential of 23 kPa) Max imax 3 Initial condition effects on bottom boundar y PAGE 75 61 Coarse sandgravel layers iapplied = icr = 0.880.2 0.1 0 0.1 0.2 0.3 0.4 0.5 35097144191 Time steps(imax(Dr) imax(Dr = 5%)) / icr Dr = 95% Dr = 5% Dr = 75% Dr =5% Dr = 50% Dr = 5% Figure 31b: Plot of Difference in Maxi mum Hydraulic Gradients Within The Coarse SandGravel Layers Vs Time Step s for Different Levels of Compaction (For Pressure Differential of 23 kPa) Figure 32 shows that at the steady state the ma ximum spatial hydraulic gradients occur at the interface. First, the locali zed hydraulic gradient is comp uted at each grid from the water pressures, velocities and its elevation using Eqn 34. At any depth in the pavement system, the average localized hydraulic gradient is determined by averaging the localized hydraulic gradient of the indivi dual grid elements in each row along the pavement width. Figure 32 shows that the localized hydraulic gradients also clearly exceed the values predicted by the conventional analysis (Darcys law). PAGE 76 62 Gravelcoarse sand layers Dr = 95% iapplied = icr = 0.88 0 0.25 0.5 0.75 1 1.25 0.20.40.60.811.21.4 Averaged localized hydraulic gradient Depth of pavement/m Localized iTransient model(at steady state) Localized iSteady state SIMPLE algorithm i across coarse sand (deterministic Darcy) i across gravel (deterministicDarcy) Figure 32: Comparison of Hydraulic Gr adients Obtained Using Continuum and Particulate Approaches (Steady State) As illustrated in Figure 33, when the soil is not uniformly compacted, some localized porosities would be high. In this study, in order to clearl y visualize the local effects, the coarse sand/gravel interface is assembled with a local region poorly compacted (Dr = 5%) compared to other regions (Dr = 95%). Coarse sand gravel interface PAGE 77 63 Figure 33: Porosity Distribution for C oarse SandGravel Interface with NonUniform Compaction (Circled Area is Packed with Dr of 5% and the Remaining Area is Packed with Dr of 95%) Figure 34 shows the impact of the localized weak zone on the flow pattern. As illustrated in Figure 35, fluid pressure varies linearly within each layer and shows a sudden change at the interface due to differences in grai n size and porosity. On the other hand, a deficiency of compaction make s the fluid pressure distribution nonlinear which affects the localized flow quantities. That is, an abrupt pressure change between two adjacent grid elements leads to undesirable higher locali zed hydraulic gradients as seen in Figure 36. An advantage of this mode l is that by making the size of the grid elements as small as possible, the desired localized effect s can be observed more effectively. Deficiency of compaction, D r = 5% PAGE 78 64 (a) (b) Figure 34: Velocity Vector Plots: (a) Uniformly Compacted Soil Media (Dr = 95%) (b) Deficiency of Compaction in Local Area (a) (b) Figure 35: Plot of Fluid Pressure Dist ribution across the Coarse SandGravel Interface: (a) Uniformly Compacte d Coarse SandGravel Interface ((Dr = 95%) (b) Deficiency of Compaction in Local Area Localized effects Localized effects PAGE 79 65 Figure 36: Plot of Localized Hydraulic Gr adients for Coarse SandGravel Interface with Weak Region 3.7 Limitations of the Model In this simplified particulate model, the fluid phase is partially coupled with the solid phase only in terms of the volumetric compatib ility. Therefore, the intricate dynamics of individual soil particle inte ractions with water, which would not impact the flow characteristics significantly, ar e not taken into account. Al so, the soil particles are assumed to be perfect spheres neglecting the angularity and ro ughness. One realizes that the angularity of soil particles could wi den the ranges of both the maximum and minimum void ratio and hence that of the natu ral void ratio. Because of the probability of having relatively larger changes in void ratio of adjacent grid elements due to angularity, excessive localized hydraulic gradients of larg er magnitude could resu lt. This issue could be addressed partly, by incorporating the compression indices (Eqn 43) modified for subrounded, subangular and angular particles as detailed in [18]. Furthermore, in this Localized effects PAGE 80 66 model, the three dimensional porosities obtai ned from particle packing are coupled with the flow constrained in the th ird direction due to the assume d two dimensional nature of the pavement system. The above 2D assump tion seems to be reasonable because the predicted hydraulic co nductivities do not de viate much from those determined using widely used empirical equations. To compute the change in porosity in the transient continuity equation, the water pressure diffe rence during the previous time step [(N1) N] was used in Eqn 16 to make the numerical scheme explicit. When the time step is adequately small, this limitation would not affect the solution si nce the variation in water pressure within small time duration can be considered linear and hence equal for all equal time steps. 3.8 Conclusions In this study, water seepage in a couple of tw olayer filter interface with different levels of compaction was analyzed using a soil pa rticulate model with Navier Stokes flow principles. The particulate effects of soil with different levels of compaction were conveniently incorporated in the model using a random packing techni que, while the flow of water within the particulate assembly wa s modeled by the Navier Stokes equations. Compaction effects were incorporated in the m odel using the concept of relative density. Two separate twolayer interfaces, i.e. coarse sand/gravel and fine sand/coarse sand were assembled using particle size distributions that satisfied the conventional filter design criteria. Then, a pressure diffe rential that co rresponded to the critical hydraulic gradient was applied across the layer inte rface. The relative density certa inly impacts the flow rate PAGE 81 67 but has an insignificant effect on the locali zed hydraulic gradient s. Furthermore, the localized hydraulic gradients determined usi ng the particulate approach differed widely from their average values across the entire layer high lighting the importance of considering localized effects in formulat ing more applicable design criteria. The preliminary results obtained from this model illustrate the advantage of the particulate approach in predicting the crit ical conditions for erosion, pi ping and clogging, especially in transient flow conditions, well before the filter fails to function. Extended and more extensive research in this direction is exp ected to provide deeper insight into more accurate filter design criteria. PAGE 82 68 CHAPTER FOUR TRANSIENT SEEPAGE MODEL FOR PARTLY SATURATED AND SATURATED SOILS USING TH E PARTICULATE APPROACH 4.1 Introduction Functionality of earthen struct ures such as dams, levees, retention ponds and pavements is determined by one dominant factor; the nature of interaction of soil particles with water flow. Hence accurate analysis of water seepage through soils is essential to achieve more economic designs of such structures. The major ity of currently available design criteria are formulated based on either empiricism or the analysis of steady state laminar flow through saturated soil continua. However, very often, fiel d observations are also used to refine or calibrate the design criteria. In the conventional models, the dynamic flow of water through soil pores is commonly ideali zed using the Darcys law. Experimental studies show that Darcys law would not be applicable to mo del transient conditions and high fluid velocities that deve lop under excessive hydrauli c gradients and unsaturated soils [1]. It is also known that it is the tr ansient flow caused by ab rupt fluctuations in groundwater conditions plays a more crucial role in determining the durability of pavement structures. Furthermore, hydraulic in frastructures like retention ponds almost always operate under transient flow cond itions. Hence, one has to replace the conventional method of analysis based on a steady state continuum to an alternative PAGE 83 69 approach with a discrete soil skeleton with the ability to incorpor ate transient effects which includes presaturation flow that gene rally precede the eventual steady state flow. In this study, the author develops a finite difference model that uses the SoilWater Characteristic Curves concepts [19] and the Navier Stoke equations for analyzing transient fluid flow through partly saturate d and saturated soils. Then, the model is applied to two specific problems: 1) Confined flow occurring in a partly saturated pavement layer during sudden rise of the groundwater table and 2) Transient seepage from retention pond in to surrounding granular soil medium. In the first case, the results show how satura tion achieves following the water porosity Vs water pressure trend defined by SWCC. In th e second case, the model is able to predict the gradual reduction in the water level of the retention pond including the location of the freesurface. 4.2 Existing Particulate Models In the conventional seepage models, the flow of water through saturated soil pores is idealized using the Darcys law based on a c ontinuum approach. In the above approach, a single representative value of hydraulic conductivity determined from field tests or laboratory tests is used to model the entire so il domain. Alternatively, if seepage through soils is analyzed using the pa rticulate approach then one can certainly account for the spatial variation of hydraulic conductivity more accurately by incorporating drag forces acting between soil particles and water in th e Navier Stokes equations. Moreover, when PAGE 84 70 using the particulate approach, the effects of in tricate particle dynami cs on the water flow can also be considered whenever needed. However, modeling of seepage through particulate media by considering soilwater interaction is relatively new to computational geomechanics. The discrete element method (DEM) does provide an effective tool to model granular soils in partic ular based on micro mechanical idealizations. El Shamy et al. [4] presented a computational micromechan ical model for coupled analysis of pore water flow and deformation of saturated granul ar assemblies. Shimizus [6] particlefluid coupling scheme with a mixed LagrangianE ulerian approach which enables simulation of coupling problems with large Reynolds num bers is implemented in PFC 2D and PFC 3D released by Itasca Consulting Group, Inc. [7]. The models used by El Shamy et al. and Shimizu are both based on the work by Anderson and Jackson [8] and Tsuji et al. [9]. Anderson and Jackson [8] modeled pore fluid motion through saturated soil mass using averaged Navier Stokes equations. Tsuji et al. [9] simulated the process of particle mixing of a twodimensional gasfluidized be d using averaged Navier Stokes equations for comparison with experiments. For all the above cited studies, granular assemblies are modeled using the discrete element model de veloped by Cundall and Strack [10] and the averaged Navier Stokes equations are discre tized using a finite volume technique on a staggered grid [11]. 4.3 Flow in Partly Saturated Soils Although partly saturated s ubsurface soils are common in most building sites, the concepts of partly saturated soil mechanics ar e rarely introduced in geotechnical designs. PAGE 85 71 Furthermore, most geotechnical engineering problems are dealt with assuming positive pore water pressures when suction can occur even in saturated soils. The solutions of many practical problems involving partly saturated soils requ ire an understanding of the hydraulics, mechanics and interfacial physics of partially saturate d soils [20]. Slope stability analysis that also incorporates unsaturated soil mechanics would provide more accurate results in most cases. However, water fl ow in partly saturated soils is, in general, complicated and difficult to de scribe quantitatively, since it often entails changes in the state and content of soil water. Figure 37: Schematic Diagr am to Show the Movement Cycle of Water from Atmosphere to the Groundwater Table [20] Figure 37 is the schematic representation of the water movement cycle between the atmosphere to the groundwater table. Below the water table, soil is generally saturated with positive porewater pressures. Immediat ely above the water table is the capillary PAGE 86 72 fringe where the degree of sa turation approaches 100 per cent with negative pore water pressures (suction). In partly saturated soils, th e two phases of air an d water coexist in the interconnected pore channels. In order to an alyze seepage through saturatedunsaturated soils under transient conditions, Fredlund and Morgenstern [21] proposed Eqn 45a using the constitutive equations relating the volumetri c strain in a soil due to total stress and pore water pressure changes. t u u m t u m y h k y x h k xa w w 2 a w 1 y x (45a) z y x3 1 (45b) Where, h total hydraulic head mean normal stress ,ua poreair pressure, uw porewater pressure, (ua uw) matric suction, kx and ky x and y directional hydraulic conductivities, wm1 and wm2 are the corresponding partia l compressibilities due to ( ua) and (uw ua) respectively. Darcys law is also applicable to flow through partly satu rated soils with the use of variable hydraulic conduc tivities with respect to suction as shown late r in Section 3.2. In saturated soils, any minor ch anges in the hydraulic conduc tivity due to changes in porosity are neglected whereas in unsatura ted soils, hydraulic conductivity is significantly affected by combined changes in the void ratio and the degree of saturation (or water content) of the soil. From a geotechnical engineering point of view, it is the flow in the water phase of the unsaturated zone that is of pr actical interest. Therefore, the flow through unsaturated soils can be simplified assuming that the air phase is continuous throug hout the soil matrix PAGE 87 73 with a pressure equal to that of the atmosp here. Using the singlephase flow approach which is accurate enough for many practical purposes, an unsaturated flow model was developed by Lam et al. [22]. If no external loads are applied on the soil mass during water flow ( 0 t ) and if the air phase is continuous and open to the atmosphere (ua = 0), Eqn 45a can be simplified to Eqn 45c. t u m y h k y x h k xw w y x 2 (45c) Ng and Shi [3] also used Eqn 48c to numerical ly investigate the stability of unsaturated soil slopes subjected to transient seepage. In Ng and Shis [3] work, a finite element model was used to investigate the influence of various rainfall events and initial ground water conditions during transi ent seepage. However, slope stability was analyzed without considering the localized zones of high pressure buildup and high hydraulic gradients within the slope. 4.4 Overview of the New Model The model presented in this paper first us es a newlydeveloped packing algorithm to randomly pack a three dimensional discrete soil skeleton resembling a natural soil deposit that exists at a given relative density or part icle size distribution. Then the model is used to determine the water flow behavior in a particulate saturated partly saturated soil medium. The seepage of water through the particulate medium is modeled using the Navier Stokes (NS) equations which are disc retized using the finite difference method (FDM) [13]. The new model is capable of pred icting both transient and steady state flow PAGE 88 74 effects under both saturated and pa rtly saturated conditions. Firs t, the model is applied to simulate flow in a single layer of granular soil resulting from of a sudden surge of groundwater. Second, the model is applied to si mulate flow around a retention pond built in granular soil. The comprehensive anal ytical procedure and the computer code developed for its implementation are illustra ted in Figure 38. The analytical procedure consists of two primary tasks such as the random assembly of th e particulate medium (granular soil) and the solution of the fluid flow governing equations. The flow chart also includes the section, equation and figure numbers corresponding to each stage. PAGE 89 75 Figure 38: Flow Chart Illustrating th e Analytical Model for Flow through Unsaturated Soils No Yes Yes No Packing to determine maximum and minimum void ratios using a particle size distribution for a insitu soil (Section 4.4) Application of MonteCarlo simulation to obtain randomly distributed insitu void ratios corresponding to insitu relative density (Section 4.4.3) Assignment of random insitu void ratios for grid elements used in fluid flow analysis Application of NavierStokes equations to each grid element (Section 4.6) Numerical discretization based on a staggered grid (Section 4.7 & Fig.45) Analysis of results Plots on degree of saturation Plots on water velocities and pressures Plots of hydraulic gradients Solid phase Fluid phase Iterative solution for water pressure (p) distribution in the grid (Eqn 74 & 75 & Fig.47) Time marching p < Tol Timebased change in water velocity ( v) < Tol Assignment of initial negative water pressure due to capillary effects in unsaturated zone (Section 4.5) Determination of initial volumetric water content (water porosity, ) using Eqn (57) Updating the water porosity (nw) using Eqn 71b. Use available coefficients of permeability data (Fig.44) to modify drag force coefficients (Eqn 66) and analyze flow through unsaturated soils (Section 4.2). Determination of instantaneous volumetric water content (Eqn 57) Solve for velocities (Eqns 60, 61) PAGE 90 76 4.5 Modeling of Soil Structure (3D Ra ndom Packing of Soil Particles) Fine sand represented by the Pa rticle Size Distribution (PSD) curve in Figure 39 is used in this study. Figure 39: Particle Size Di stribution for Fine Sand 4.5.1 Simulation of Maximum and Mi nimum Void Ratios Using PSD In this model the soil particles are assumed to be of spherical shape. To determine the maximum and minimum void ratios (emax and emin) corresponding to fine sand in Figure 39, a customized random packing algorithm wa s developed. Using this algorithm, 3 mm 3 mm 3 mm cubes (Figure 40a) were packed using soil particles picked from the PSD in Figure 39. It must be noted that the PSD in Figure 43 is in fact the cumulative probability distribution of the fine sand particles. Therefor e, by using an adequately large array of random numbers from a uniform dist ribution ranging from 0 and 100 (the y axis range in Figure 39), one can select the corr esponding array of particle diameters that PAGE 91 77 conforms to the selected PSD, from the x axis. This technique, known as the MonteCarlo simulation [14] is used to select the arra y of packing diameters for each cube (Figure 40a). It is noted that the resulting maximum and minimum void ratio distributions change in each packing trial since the order of diameters used for packing is changed randomly. An adequate number of particle s must be within the cube to determine the representative emax and emin for the fine sand. 4.5.2 Implementation of the Packing Procedure In order to obtain the loosest state of fine sand (emax) within the cubes described above, different sizes of soil particle s inscribed in boxes are packed as indicated in Figure 40a using a MATLAB code developed by the auth ors. The side length of each box is the same as the diameter of the inscribed soil particle. Based on the minimum particle diameter of fine sand (0.04 mm), the side le ngth of each cube is divided into a finite number of subdivisions (F igure 40a). Then, the packing algorithm tracks the total number of subdivisions occupied by each in coming box, i.e. each pack ed soil particle, as packing proceeds based on the MonteCarlo si mulation corresponding to a given PSD (Figure 39). Finally, the automa ted algorithm fills all subdivisions until it finds that not a single subdivision is available within the considered cube for further packing of soil particles from the given PSD. As packing of each cube approaches completion, the void ratio reaches a constant value of 0.9098 whic h can be identified as the maximum void ratio for simple cubical packing [15]. PAGE 92 78 On the other hand, in order to obtain the minimum void ratio (densest packing), the maximum number of spheres smaller than th e inscribing sphere in each box is also packed into the unoccupied spac e at the corners of each box (Figure 40c and 40d) before that box is placed in the cube (Figure 40b). 50 such cubes (trials) were packed randomly using the MonteCarlo simulation. Since the minimum void ratio obtained in each trial would be different depending on how the uno ccupied space in each box filled, it can be considered as random variable. The range derived for emin as shown in Table 5 agrees with the typical values in [15]. Table 5: Soil Characteristics Soil type Fine sand Sizes of particle (mm) (Figure 39) 0.04 0.62 D15 (mm) (Figure 39) 0.16 D50 (mm) (Figure 39) 0.29 D85 (mm) (Figure 39) 0.49 emax (Packing) 0.9098 emin (Packing) 0.54 0.79 Mean of emin 0.48 Standard deviation of emin 0.0037 PAGE 93 79 Figure 40: Packing Procedure 4.5.3 Determination of the Natural Void Ratio Distribution The property of relative density is helpful in quantifying the level of compaction of coarsegrained soils. The relative density of a coarsegrained soil expresses the ratio of the reduction in the voids in the current stat e, to the maximum possible reduction in the voids (Eqn 2). The insitu void ratio distribution (e) co rresponding to a given relative density (Dr) is obtained using the previously obtained emax and emin (Eqns 46 or 47). min max maxe e e e Dr (46) or PAGE 94 80 r rD e e D e 1max min (47) Where, emax is the void ratio of fine sand in its loosest state (= 0.9098). emin is the randomly distributed void ratio of fine sand in its densest state. e is the randomly distributed void ratio of sand in its natural state in the field. For a relative density of 50%, Eqn 51 reduces to 45 0 e 5 0 emi n (48) The mean ( e) and the standard deviation (Se) of natural void ratios are determined from Eqns 49 and 50 ) 1 ( 9098 0min r rD e D e (49) min e r eS D S (50) Where, mi n e (Eqn 51a) and Se min (Eqn 51b) are the mean and standard deviation of the minimum void ratio (Table 5) obtained from the packing procedure described in Section 4.4. 50 min min50 1i ie e (51a) 50 2 min min min50 1i i ee e S (51b) Knowing the mean and standard deviation of the natural void ratio (Eqns 49 and 50) and assuming a normally distributed probability density function for the natural void ratio (e), its cumulative distribution (Figur e 41) is obtained using Eqn 52. eZ S e e (52a) PAGE 95 81 Cumulative frequency, *Z Z P Z (52b) Where, Z* is the standard normal variate corresponding to a cumulative probability of *Z in the standard normal distribution (Z). The cumulative distribution of natural void ratios corresponding to a Dr of 50% is shown in Figure 41. Figure 41: Cumulative Distribution for Natura l Void Ratios for the Fine Sand with Dr = 50% Then, once more the MonteCarlo simulation can be used to pick an array of natural void ratios for fine sand using Fi gure 39. Since the NS equatio ns are written in terms of porosities (Section 4.6), the spatial probability distribution of natural void ratios is first converted to porosities using Eqn 53. e 1 e n (53) PAGE 96 82 The spatial distributions of natural void rati os (porosities) so obtained assumed to be representative of the fine sand will be used in the seepage analysis model presented in Section 4.6. 4.6 Modeling of Fluid Flow in Partially Saturated Soils When the developed model is applied to a pa rtially saturated gran ular soil layer, the initial water pressures within the saturated and unsaturated zones are assumed to be due to hydrostatic pressure and the capillary suct ion effects respectively, as shown in Figure 42. Figure 42: Initialization of Water Pressure in the Saturated and Unsaturated Zones Saturated soil Unsaturated soil Capillary fringe soil Suction due to capillary Hydrostatic pressures Groundwater table PAGE 97 83 4.6.1 Soil Water Characteristic Curve/Water Moisture Retention Curve The Soil Water Characteristic Curve (SWCC) [1 7] depicts the relationship between soil water content ( ) and soil water pressure potential (suction, ). The soil parameters defining a typical Soil Water Characteristic Cu rve are shown in Figure 43a while the Soil Water Characteristic Curves for different pa rticle size distributions are presented in Figure 43b. From Figure 43b, it is seen that sands and grav els lose water quite readily upon suction induced drainage while loams and clays lose much less water upon drainage. There are a number of empirical eq uations proposed to bestfit the Soil Water Characteristic Curves [19] using the soil parame ters shown in Figure 43a. In this model, Eqn 54a Eqn 54c are used to determine the volumetric water content ( ) for a particular suction ( ). Volumetric water content, V Vw (54a) Where, V is the total volume of a soil sample and Vw is the volume of water in that sample. Figure 43a: Typical Soil Water Characteristic Curve and wm2for a Saturated Unsaturated Soil [10] Figure 43b: Soil Water Characteristic Curve for Specific Soil Types [23] PAGE 98 84 m n sa e C m n a ln , (54b) Where, C( ) is a correction function defined as: 1 10 1 ln 1 ln6 r rC (54c) Where, e = 2.71828; r is the suction corresponding to the residual water content, r. s is the (saturated) volumetric water content at zero pressure and a, m, and n are fitting parameters. In this model, s in each grid element is set to the porosity of that element obtained from initial packing (Section 4.4). In the particulate approach, although the s varies with the spatial distribution of porosity, the average s corresponding to 50% of degree of compaction is used in this study for simplicity. The parameters for Eqn 54b Eqn 54c relevant to the fine sand modeled in the current study are shown in Table 6. Table 6: Parameters Used to Define the SWCC for Sand [23] r 50 kPa Average s 0.45 a 1.948 kPa m 2.708 n 1.084 PAGE 99 85 In Section 4.6, will also used as the water porosity in formulating the continuity equation and Navier Stokes eq uations for partly saturated soils. Furthermore, using fundamental soil mechanics, the degree of saturation can be shown that n S (55) Where, Sdegree of saturation, n total porosity of the soil. 4.6.2 Determination of the Unsaturated Permeability Function The hydraulic conductivity is commonly referred to as the coefficient of permeability in geotechnical engineering. In unsaturated soils, the coefficient of permeability is primarily determined by the poresize distribution and the degree of saturation. The permeability functions for unsaturated soils are used to represent the relationship between the coefficient of permeability and soil suction. Figure 44 illustrates the variation of the coefficient of permeability with suction for different types of soils. There are several empirical equations and statistical models us ed to determine the permeability functions [19]. One such approach uses the Soil Water Characteristic Curve (Figure 43a) to predict the permeability function [19]. PAGE 100 86 Figure 44: Unsaturated Permeability F unction for Selected Soil Types [23] The flow model presented in this paper does not directly involve the use of the coefficient of permeability. Alternatively, in this mode l the authors use an appropriate plot from Figure 44 to modify the drag forces in the Navier Stokes equations to account for the variation of the permeability of uns aturated soils with suction. 4.7 Development of th e Analytical Model 4.7.1 Navier Stokes Equations Since the flow model presented in this pape r is twodimensional, the three dimensional porosities obtained from particle packing ar e coupled with two di mensional water flow. In this study, it is also assumed that the water paths through the pore channels are PAGE 101 87 continuous for partly saturated soils near satu ration. Hence, Navier Stokes equations can be written as follows in a generalized form in terms of the water porosity, (= Sn) expressed by Eqn 56: Water mass Conservation (Continuity Equation): 0 y v n S x u n S t n S (56) Momentum Conservation (Momentum Equations): X direction: 2 2 2 2 2y u n S x u n S D x p n S y v u n S x u n S t u n Sx (57) Y direction: 2 2 2 2 2y v n S x v n S D y p n S y v n S x v u n S t v n Sy yg n S (58) Where, n and S porosity and degree of satu ration at the location (x, y) at time t respectively, u, v fluid velocities in the x and y directions respectively, fluid density (1000 Kg/m3), p fluid pressure, fluid viscosity (103 Pa.s), gy gravitational force per unit mass in the y direction (9.81m/s2). The following conditions are assumed in the analysis: 1) Under saturated conditions: S = 1 and n vari es due to the compressibility of the soil skeleton as described in Section 4.4. 2) Under unsaturated conditions: n is constant and S increases as suction decreases. PAGE 102 88 4.7.2 Modification of Drag Force for Water Flow through Unsaturated Soils The authors developed a modified expression for the drag forc e in unsaturated soils based on the variation of the hydraulic conductivi ty of unsaturated soils. Averaged waterparticle interactions (drag forces) can be qua ntified using the semiempirical relationship provided by the standard Ergun equation [6]. Eqn 59 shows the drag force modified by the authors for application to unsaturated soils. f u d n S n S 1 75 1 u d n S n S 1 150 D2 p 2 2 p 2 2 x (59) Where, Dx averaged waterparticle interaction force per unit volume of soil, pd averaged particle diameter. f is the correcti on factor to account for the space occupied by air that does not introduce a significant drag force. A similar expression for Dy is used for y directional averaged wate rparticle interactions. For fully saturated soils, S = 1.0 and f = 1.0 with no air. Assuming that the flow through soil pores to be primarily laminar flow with low velocities, the nonlinear term co ntributing to the drag force can be neglected resulting in the simplified form of Eqn 6 for unsaturated soils as follows: u d n S n S f Dp US x 2 2 21 150 (60) Similarly, the drag force under saturated condi tions for the same porosity and flow can be written as PAGE 103 89 u d n n Dp S x 2 2 21 150 (61) Drag forces cause resistance to water flow and hence they can be considered as being inversely proportional to the coefficien t of hydraulic condu ctivity (K) as, f K K D DUS S S x US x (62) From Eqns 60 and 61, 2 2 21 1n S n S K K fUS S (63) Where, KS and KUS are the hydraulic conductivities of saturated and unsaturated soils respectively. 4.8 Numerical Solution Technique A finitedifference scheme is used to discre tize the Eqns 56, 57 and 58. This scheme is based on forward difference in time and poro sity, and central difference in space with a twodimensional staggered mesh arrangement [11, 13]. Eqn 64a expresses the x directional momentum equation in the numerical form for the sequential time steps of N and N+1 referring to the notation in Figure 45. PAGE 104 90 Figure 45: Computational Grid for the XMomentum Equation j i p j i j i j i j i j i j i j i j i j i j i j i j i j i j i j i N j i j i N j i j iu d Sn Sn f x p p Sn y v u Sn v u Sn x u Sn u Sn t u Sn u Snj i, 2 1 2 2 2 , 1 2 1 2 1 1 2 1 1 2 1 2 1 1 2 1 1 2 2 1 2 1 2 2 3 2 3 2 1 2 1 1 2 1 2 1,1 150 2 2 2 j 2 1 i j i 1 j 2 1 i 1 j i 1 j 2 1 i 1 j i 2 j 2 1 i j i j 2 1 i j 1 i j 2 3 i j 1 iy u Sn 2 u Sn u Sn x u Sn 2 u Sn u Sn (64a) Where, 2 1 1 2 1 2 1 2 12 1j i j i j iv v v (64b) PAGE 105 91 2 1 1 2 1 2 1 2 12 1j i j i j iv v v (64c) Within the staggered grid, some velocities need to be interpolated as shown in Eqns 64b and 67c. In a concise form, Eqn 64 a can be rewritten as t A x p p t A t A u Sn u SnN i N j i N j i N j i N j i N j i j i N j i j i 2 1 , 1 2 1 2 1 2 1 2 1 1 2 1 2 1) ( (64d) Where, y 2 v u Sn v u Sn x 2 u Sn u Sn A2 1 j 2 1 i 1 j 2 1 i 1 j i 2 1 j 2 1 i 1 j 2 1 i 1 j i 2 j 2 1 i j 2 1 i 2 j 2 3 i j 2 3 i j 2 1 i (64e) j 2 1 i j 2 1 iSn A (64f) and j i p j i j i j iu d Sn Sn f Aj i, 2 1 2 2 2 , 2 1,1 150 2 2 1 1 2 1 1 1 2 1 1 2 2 1 , 2 1 1 2 3 12 2 y u Sn u Sn u Sn x u Sn u Sn u Snj i j i j i j i j i j i j i j i j i j i j i j i (64g) The numerical form of the y directional mo mentum equation can be written similarly based on Eqn 58 as Eqn 65. PAGE 106 92 t g Sn t B y p p t B t B v Sn v Sny 2 1 j i N 2 1 j i N j i N 1 j i N 2 1 j i N 2 1 j i N 2 1 j i 2 1 j i 1 N 2 1 j i 2 1 j i (65) The discretized continuity equation of water is given by Eqn 66. In computing the change in porosity in the transient continuity equati on, the numerical scheme is made explicit by using porosity during the previous time step (N1) in Eqn 66. y v Sn v Sn x u Sn u Sn1 N 2 1 j i 2 1 j i 1 N 2 1 j i 2 1 j i 1 N j 2 1 i j 2 1 i 1 N j 2 1 i j 2 1 i t Sn Sn t SnN j i N j i 1 , (66) Finally, using either Eqn 67a or Eqn 67b, the change in water porosity (Sn ) during a time step, t, can be determined and used in the di scretized continuity equation (Eqn 66). In the saturated zone (p > 0 and S = 1), the change in porosity is determined using Eqn 67a. t p p p n C t nN j i N j i N j i N j i j i c 1 1 4343 0, 1 , 2 (67a) Where, Cc equivalent compressi on index of the soil, total vertical stress in soil. Meanwhile, in the partly satu rated zone (p < 0 and S < 1), the SoilWater Characteristic Curve is used to determine the change in water porosity using Eqn 67b. PAGE 107 93 t p p m t SnN j i N j i j i w 1 , 2 (67b) The slope of the Soil Water Characteristic Curve for fine sand, wm2,at a particular suction is obtained by curvefitting the correspondi ng SWCC assumed for fine sand (Figure 46) in this study. Figure 46: Soil Water Charact eristic Curve for Fine Sand Now, Eqn 64 and Eqn 65 can be written fo r two adjacent control volumes as shown in Eqn 68 and Eqn 69. t A x p p ) t ( A t A u Sn u SnN 2 1 i N j 1 i N j i N j 2 1 i N j 2 1 i N j 2 1 i j 2 1 i 1 N j 2 1 i j 2 1 i (68) PAGE 108 94 y p p t B t B v Sn v SnN 1 j i N j i N 2 1 j i N 2 1 j i N 2 1 j i 2 1 j i 1 N 2 1 j i 2 1 j i t g Sn t By 2 1 j i N 2 1 j i (69) By substituting the Eqns 64, 65, 68, 69, 67a and 67b into Eqn 66 and rearranging the terms results in Eqn 70 and Eqn 71 for sa turated zone and par tly saturated zone respectively. Figure 47: Neighborhood Points Us ed in the Iterative Procedure For the saturated zone: 0, 1 , 1 , 1 , 1 , j i j i j i j i j i j i j i j i j i j i j iT p e p d p c p b p a (70a) Where, N j i j i j i c j i j i j i j i j ip n C t y n n x n n a. 2 2 2 2 1 2 1 2 2 1 2 1 ,1 4343 0 (70b) 2 2 2 1 ,t x n bj i j i (70c) PAGE 109 95 2 2 2 1 ,t x n cj i j i (70d) 2 2 2 1 ,t y n dj i j i (70e) 2 2 2 1 ,t y n ej i j i (70f) and 1 , 2 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 2 1 2 1 2 2 1 2 1 2 2 1 2 1 2 2 1 2 1 ,1 4343 0 N j i N j i j i j i c N j i j i N j i j i N j i j i N j i j i N j i N j i N j i N j i N j i N j i N j i N j i j ip p n C t y v n v n t x u n u n y t B B y t B B x t A A x t A A T y t g n ny j i j i 2 2 1 2 1 (70g) For the partly saturated zone: 0 T p e p d p c p b p a' j i 1 j i j i 1 j i j i j 1 i j i j 1 i j i j i j i (71a) Where, j i w 2 2 2 2 1 j i 2 1 j i 2 j 2 1 i j 2 1 i j im t y Sn Sn x Sn Sn a (71b) PAGE 110 96 2 2 j 2 1 i j it x Sn b (71c) 2 2 j 2 1 i j it x Sn c (71d) 2 2 2 1 j i j it y Sn d (71e) 2 2 2 1 j i j it y Sn e (71f) and y t g Sn Sn p m t y v n S v Sn t x u Sn u Sn y t B B y t B B x t A A x t A A Ty j i j i N j i j i w N j i j i N j i j i N j i j i N j i j i N j i N j i N j i N j i N j i N j i N j i N j ij i 2 2 1 2 1 1 , 2 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 2 1 2 1 2 2 1 2 1 2 2 1 2 1 2 2 1 2 1 ', (71g) PAGE 111 97 4.9 Numerical Application of the Model 4.9.1 Modeling Confined Flow through a Part ly Saturated Fine Sand Layer (Case I) Figure 48: Illustration of the Flow through the Unsaturated Fine Sand Layer The fine sand layer (Figure 48) was divided into a 50 50 gr id for the finite difference analysis. One fluid cell, known herein as a gr id element, was selected to be 2.5 mm 2.5 mm 1m. A unit width has been selected in the third direction assuming two dimensional water flow. First, natural porositie s (void ratios) were assigned to each grid element based on the spatial probability distribution of natural void ratios determined from packing (Section 4.4). Within each grid element, the fluidparticle interaction is quantified considering the averaged particle diameter which is defined as the arithmetic mean of the diameters of all the particles in that grid. 0.125 m PAGE 112 98 4.9.2 Boundary Conditions The boundary conditions appropriate for the wate r flow through the si ngle soil layer are specified in the computer code as indicated in Figure 49. 1) At the inflow boundary, the pressure and velocities are specified. 2) At the outflow boundary, the pressure is specified and veloc ity components are allowed to float. 3) At the vertical walls, the slip condit ions are maintained. Thus, the velocity component, u, normal to th e walls is set to zero. Figure 49: Boundary Conditions Incorporat ed in the Flow through Saturated Unsaturated Soil Layer A PAGE 113 99 4.9.3 Results Upon the application of the elevated water pressure that causes a sudden hydraulic gradient across the pavement layer in Figu re 50, the variation of water porosity (volumetric water content) with water pressure is illustrated in Figure 50. Figure 50 also shows that suction within the soil mass is gra dually converted to a po sitive pressure under the application of elevated pressure at the bottom of the soil layer. The corresponding variations of degree of satu ration are shown in Figure 51. 0.4 0.41 0.42 0.43 0.44 0.45 0.46 15001000500050010001500200025003000Water pressure (Pa)Water porosity Figure 50: The Variation of Water Porosity with Water Pressure at a Point within the Unsaturated Zone (At Point A in Figure 49) PAGE 114 100 90 91 92 93 94 95 96 97 98 99 100 15001000500050010001500200025003000Water pressure (Pa)Degree of saturation Figure 51: The Variation of the Degree of Saturation with Water Pressure (At Point A in Figure 49) Figures 52a 52c show how the partly satura ted zone gradually becomes fully saturated due to the upward water seepage. From Figure 53, it is seen that it takes approximately 50 time steps (5103 seconds) for complete saturation of the layer and 80 time steps (8103 seconds) to reach the steady state fl ow through the fully saturated zone. PAGE 115 101 Figure 52: Velocity Vector Plots for Flow through Saturated and PartlySaturated Soil Layer (d) t = 100 (a) t = 1 (c) t = 50 (b) t = 25 PAGE 116 102 Figure 53: Variation of Discharge Flow Rate with Time Step 4.10 Conclusions The results show that 1) The partly saturated zone having capillary pressures (suction) gradually becomes fully saturated due to seepage caused by an elevated water pressure at the bottom of the fine sand layer. 2) It takes the fine sand layer almost 80 time steps (8103 seconds) to become fully saturated. 3) The water porosity Vs suction plot and th e degree of saturation Vs suction plot show good agreement with the corresponding characteristics plots for fine sand. Based on the above it can be concluded that mo del is applicable to partly saturated soils as well. Time of saturation Time of steady flow PAGE 117 103 CHAPTER FIVE MODELING UNCONFINED FLOW AROUND A RETENTION POND 5.1 Introduction Retention ponds are manmade or natural de pressions into which stormwater runoff is directed for temporary storage [Figure 54] with the expectatio n of disposal by infiltration into a shallow groundwater aquifer. They are often created near areas of development and in many instances required with new developmen t of buildings, parking lots, roads, etc by the permitting agencies. Retention ponds are developed primarily to serve two functions such as limit flooding and the removal of pollutants. Figure 54: Typical Retention Pond PAGE 118 104 5.2 Existing Design of Retention Pond The required size of a retention pond depends on the rate of inflow to the pond and the total quantity of flow as well as the rate of infiltration during a storm event, given the antecedent conditions of the receiving aquifer. In order to avoid the overflow or flooding, the retention volume recovery time should be adequate to completely dissipate the retained stormwater from th e pond after the design storm event. The retention volume recovery time is defined as the time it takes to infiltrate the retention volume. This depends on subsurface soil conditions and the input due to surface runoff. In the conventional models that are used to design retention ponds, the soil skeleton is treated as a continuum and the dynamic flow of water through soil pores is commonly idealized using the Darcys law. MODFLOW [24] is a U.S. Geological Surveys modular finitedifference saturated flow model, which is currently the most widely used numerical model in the U.S. Geological Survey (USGS) for analyzing groundwater flow problems. The above program is used by hydrogeologi sts to simulate the flow of groundwater through aquifers. In MODFLOW, a groundwat er flow equation derived from Darcys law is solved using the finitedifference approximation. The partial differential equation of groundwater flow used in MODFLOW is t h S W z h K z y h K y x h K xs z y x (72) Where, Kx, Ky and Kz are values of hydraulic conductivity along x, y and z coordinate axes respectively; h is the hydraulic head ; W is a volumetric flux per unit volume representing sources and/or sinks of water (W < zero for flow out of the groundwater PAGE 119 105 system and W > zero for flow in) and Ss is the specific storage of the soil. Eqn 72, when combined with boundary and in itial conditions, describes transient threedimensional groundwater flow in a heteroge neous and anisotropic medium. MODRET (Computer MODEL to Design RETENTION Ponds) is an interactive computer program that can be used to calcula te the unsaturated and saturated infiltration losses from stormwater retention ponds into un confined aquifers usi ng a modified version of MODFLOW. MODRET was originally developed in 1990 as a complement to a research and development project for South West Florida Water Management District (SWFWMD). The Infiltration Module of the MODR ET program uses a modified Green and Ampt infiltration equation [25] to calculate the unsaturated infiltration and a modified USGSs model [24] to calculate th e saturated infiltration. The Green and Ampt infiltration equation [25] was originally pr esented as an empirical description for unsaturated infiltration analysis; later, its empirical constants were theoretically investigated [26]. For all a bove analysis, the coefficients of hydraulic conductivities determined from the field tests or laborato ry experiments are used in determining the infiltration losses through retention ponds. But, instead of using the a single value of coefficient of hydraulic conduc tivity for entire flow domain on a continuum basis, the current particulate model uses spatially locali zed coefficients of hydraulic conductivities in terms of drag forces between solid part icles and pore fluid which are instantaneously determined from the average particle diameters and the localized porosities. PAGE 120 106 5.3 Modeling of Soil Structure (3D Ra ndom Packing of Soil Particles) The particle size distribution curve shown in Figure 55 was selected to model the retention pond. The soil partic les are assumed to be of sphe rical shape and to determine the natural void ratio for coarse sand shown in Figure 55, a customized random packing algorithm was developed. In the case of co arse sand, 16 mm 16 mm 16 mm cubes were packed. It must be noted that the PSD in Figure 55 is cumulative probability distribution of particle size fo r coarse sand. Therefore, by using an adequately large array of random numbers from a uniform distributi on between 0 and 100 (t he y axis range in Figure 55), one can select the co rresponding array of particle diameters that conforms to a selected PSD, from the x axis. This t echnique, known as the MonteCarlo simulation [14] is used to select the a rray of packing diameters for each cube. It is noted that the resulting natural void ratio distribution cha nge in each packing trial since the order of diameters used for packing is changed random ly. The detailed packing procedure can be found in [27]. A thousand of such cubes were ra ndomly packed in orde r to model the half of the retention pond defined by mesh of 100 by 150. The asymmetry of the retention pond caused by the nature of random pack ing is neglected in this study. PAGE 121 107 0 10 20 30 40 50 60 70 80 90 100 1 10Grain size/(mm)% Cumulative frequenc Figure 55: Particle Size Dist ributions for Coarse Sand 5.4 Numerical Procedure for Determination of the FreeSurface Determination of the freesurface is required in the analysis of unconfined seepage problems such as flow through earthen dams and retention ponds. One commonly used approach to determine the fr eesurface is the graphical approach based on Dupuits theory [1]. In addition, there are numerical techniques proposed by various researchers in the literature [28 32]. In this study, a modified technique resembling the extended pressure (EP) method first propos ed by Brezis et al. [33] is employed with Navier Stokes equations in order to determine the freesu rface for flow around a retention pond. The dry zone and wet zone around the retent ion pond are defined using the fundamentals of streamlines, i.e., the locus of the path of flow of an indi vidual particle of water and continuity of flow quantity through a cross se ction. Along any flow line, the path of PAGE 122 108 water passing through point P(x1, y1) can be expressed (Figure 56) in mathematical form as 1 1u v dx dy (73) Hence knowing one point, P(x1, y1) on the freesurface, the next point, Q (x2, y2) along the freesurface can be obtained using Eqn 74. 1 2 P 1 1 1 2y y v u x x (74a) y v u x xP 1 1 1 2 (74b) In generalized form, y v u x xj 1 1 j 1 j (74c) Once, the extreme boundary point, Q, is dete rmined along the freesurface, the dry zone is identified using a Heaviside step function in Section 5.5. PAGE 123 109 Figure 56: Determination of Dry and Wet Zones around the Retention Pond 5.4.1 Navier Stokes Equations In the current formulation, water seepage is analyzed using Navier Stokes equations written in terms of water pressures and veloci ties. Hence, the author s use a Heaviside step function with water pressures and velocities in order to identify the dry and wet zones in conjunction with Navier Stokes equations. For this wo rk, the Heavisidestep function (Eqn 75) is expressed in terms of the coordinates of the extreme boundary point on the freesurface at any horizontal level (xj) as follows: 0 if x > xj H (x) = (75) 1 if x < xj To achieve this objective, Navier Stokes equations are written incorporating the Heaviside step function as follows: Freesurface Dry zone Wet zone x y v1 u1 Q(x2, y2) P(x1, y1) Groundwater level Center line (xj,yj) (xj+1,yj+1) PAGE 124 110 Mass Conservation (Continuity Equation): t n y Hv n x Hu n (76) Momentum Conservation (Momentum Equations): X direction: x p Hn HD y Hu n x u nH y v u nH x Hu n t u nHx 2 2 2 2 2 (77) Y direction: y yng y p Hn HD y Hv n x Hv n y Hv n x v Hu n t Hv n 2 2 2 2 2 (78) It is seen that Eqns 76 78 are inactivated in the dry zone (i i*). This can be assured by using the Heavisidestep functi on in Eqns 75. Now, Eqns 76 78 can be applied both in the wet and dry zones. 5.4.2 Boundary Conditions The boundary conditions appropriate for the act ual water flow around the retention pond are specified in the computer co de as indicated in Figure 57. 1) Hydrostatic pressure distribution (h w) is specified along the retention pond walls based on the water level in the pond (h). 2) Along the walls of the retent ion pond, nonslip velocity boundary conditions are specified as shown in Figure 57. PAGE 125 111 3) At the exterior vertical wall boundaries within the twophas e media, the slip conditions are maintained. Thus, the velocity component, u, normal to the walls is set to zero. 4) At the groundwater table, water pressure is set to zero; water velocities allow floating. 5) As an initial condition, the completely filled retention pond is considered. Figure 57: Boundary and Initial Conditions Incorporated in the Flow around the Retention Pond Groundwatertable PAGE 126 112 5.5 Determination of FreeSurface from Dupuits Theory Figure 58: Analytical Method to Determine the FreeSurface Using the configuration of th e retention pond irrespective of soil properties, the freesurface is analytically de termined using Eqn 79. (79a) (79b) 5.6 Results on FreeSurface The freesurfaces around the retention pond obt ained using the numerical procedure and Dupuits theory are compared as shown in Figure 58. There are discrepancies which can be due to the fact that the Dupuits theo ry does not take into account the hydraulic conductivities around the retention pond whereas the numerical procedure incorporates 2 2 L e L B xL y H B L 2 PAGE 127 113 the localized variations of hydraulic c onductivities. Therefore, the freesurface determined from the numerical method is more accurate than that of Dupuits theory. Moreover, it is seen from Figure 59, that th e influence zone covered by the freesurface is getting wider as time increases. Figure 59: Comparison of FreeSurface around the Retention Pond PAGE 128 114 Figure 60: Location of FreeSurface around the Retention Pond with Time 5.7 Prediction of Recovery Time Recovery time is defined as a length of tim e required for the design treatment volume in a pond to subside to the normal level or bo ttom of the pond. Threedimensional Finite Difference Groundwater Flow Model (MODF LOW) developed by U.S. Geological Survey [24] for saturated groundwater flow modeling is used to compare the recovery time predicted using the model developed in this study. A schematic of an unconfined aquifer is shown in Figure 60 with the notations used in the equations. The aquifer is assumed to be homogeneous and isotropic, and it is underlain by an impervious horizontal layer. PAGE 129 115 Figure 61: Saturated Infiltrati on Analysis using MODFLOW Dimensionless parameters, Fx and Fy are defined in terms of 1) Hydraulic parameters of the aquifer 2) Recharge rate and duration, 3) The physical configuration of th e pond and the desired time period. (80a) (80b) (80c) Where, hc Height of water level in the re tention pond above the initial groundwater table, hb Height of the pond bottom a bove the groundwater table, hv Maximum height 2 1 24 t D K W FH x T c yH h F v b Th h H PAGE 130 116 of water level in pond, W Average width of the pond, KH Average horizontal hydraulic conductivity of the aquifer, D Aver age saturated thickness of the aquifer (= H + hc/2), H Initial saturated thickness of aqui fer, f Effective porosity and t Recovery time. A family of dimensionless curves similar to Figure 63 is plotted by USGS for different effective porosity of soil. According to E qn 80, the estimated recovery time for the selected pond is 4.7 days whereas the numerical model predicts that of 4 days from the Figure 62. Determination of recovery time0 0.2 0.4 0.6 0.8 1 050000100000150000 Time (sec)Water level in Pond (m ) Figure 62: Prediction of Recov ery Time from Numerical Model PAGE 131 117 Figure 63: Sample Plot to Determin e Recovery Time from MODFLOW PAGE 132 118 5.8 Conclusions The application of the model de veloped in the first phase of the study has been extended to predict the recovery time for the retent ion pond including the location of the freesurface. The recovery time obtained from the developed model agrees with that of the MODFLOW analysis. Furthermore, the fr eesurface around the rete ntion determined from the numerical model shows good agreem ent with that predic ted by the Dupuits theory. PAGE 133 119 REFERENCES 1. Harr ME. Groundwater and seepage. Dover publications: New York, 1991. 2. Fredlund DG. Second symposium and short course on unsaturated soils and environmental Geotechnics. University of Saskatchewan, Canada, 2003. 3. Ng CWW. and Shi Q. A Numerical Investig ation of the Stability of unsaturated soil slopes subjected to transient seepage. Co mputers and Geotechni cs 1998; 22(1):118. 4. El Shamy U, Zeghal M, Shephard M, Dobry R, Fish J, and Abdoun T. Discrete element simulations of water flow through granular soils. 15th ASCE Engineering Mechanics Conference, 2002. 5. El Shamy U., and Zeghal M. Coupled C ontinuumDiscrete model for saturated granular soils, Journal of Engineering Mechanics, April 2005. 6. Shimizu Y. Fluid coupling in PFC2D and PFC 3D, Manual, Itasca Consulting Group, Numerical Modeling in Micromechanic s via Particle Methods, 2004. 7. Itasca. Particle Flow Code PFC3D, release 3.0. Itasca Consulting Group, Inc., Minneapolis, Minnesota, 2003. 8. Anderson TB and Jackson R. A Fluid M echanical Description of Fluidized Beds, Industrial Engineering Chemical Fundamentals, 1967; 6(4): 527539. 9. Tsuji Y., Kawaguchi T., and Tanaka T., Di screte particles simulation of 2dimensional fluidized bed. Powder T echnology, 1993; 77: 7987. 10. Cundall PA and Strack ODL. A Discrete Numerical Model for Granular Assemblies, Geotechnique 1979; 1: 4765. 11. Patankar S. Numerical Heat Transfer a nd Fluid Flow: Taylor and Francis, London. 1980. PAGE 134 120 12. Engineer Manual No. 111021901, Depart ment of the Army, U.S. Army Corps of Engineers, Washington, DC, April, 1993. 13. Anderson J. Computational Fluid Dynamics the Basics with applications, McGraw Hill, Inc. 1995. 14. Harr ME. Mechanics of Particulate Medi a. A Probabilistic Approach: New York, McGraw Hill, 1977. 15. Bowles JE. Physical and Geotechnical Properties of Soils, McGraw Hill, Inc.1984. 16. Das BM. Principles of Geotechnical E ngineering: Brooks/Cole Thomson learning, sixth edition, 2006. 17. Jeyisanker, K. and Gunaratne, M, Anal ysis of Steady State Water Seepage in a Particulate Pavement System, under review by Transportation Research Board, 2008. 18. Pestana JM, Whittle AJ. and Salvati LA. Evaluation of a constitutive model for clays and sands: Part 1 Sand behavior, Internati onal Journal for Numeri cal and Analytical Methods in Geomechanics 2002; 26:10971121. 19. Fredlund DG, and Xing A., Equations for the SoilWater Characteristic Curve. Canadian Geotechnical Journal, 31; 3:521532. 20. Lu.N, Likos WJ, Unsaturated Soil M echanics, Wiley Publishers, May, 2004. 21. Fredlund DG., and Morgenstern NR., Stre ss State Variables for Unsaturated Soils, Journal of the Geotechnical Engineering Divisi on, American Society of Civil Engineers, 103;GT5: 447466. 22. Lam.L, Fredlund, DG., and Barbour.SL ., Transient Seepage Model for SaturatedUnsaturated Soil Systems: A Geotechni cal Engineering Approach, Canadian Geotechnical Journal, 1987; 24: 565580. 23. Clifton AW, Barbour SL, and Wilson GW. The Emergence of Unsaturated Soil Mechanics, Fredlund Volume, NRC Research Press, Ottawa, 1999. 24. Harbaugh, A.W., and M.G. McDonald. 1996. User's documentation update to the U.S. Geological Survey modul ar finitedifference U.S. Geological Survey OpenFile Report 96485, 56p. 25. Green, W. H., and Ampt, G. A. Studies on soil physics. 1. The flow of air and water through soils. J. Agric. Sci., 1911, 4(1), 124. PAGE 135 121 26. Bouwer, H., Infiltration of Wa ter into Nonuniform soil, Journal of the Irrigation and Drainage Division, ASCE, 1969, No. IR4, 451462. 27. Jeyisanker, K. and Gunaratne, M. Analys is of Water Seepage in a Pavement System Using the Particulate Approach, accept ed for the publication, Computers and Geotechnics, 2008. 28. Taylor RL, Brown CB, Darcy flow solutions with a free surface, Journal of Hydraulic Division, ASCE, 1967; 93: 2533. 29. Bardet J.P, Tobita T, A Practical Method for Solving FreeSurface Seepage Problems. Computers and Ge otechnics, 2002; 29; 451475. 30. Desai CS, Li GC A residual flow proce dure and application for free surface flow in porous media. Advanced Wate r Resources 1983; 6; 2735. 31. Alt HW, Numerical Soluti on of Steady State Porous Fl ow Free Boundary Problems. Numerical Mathematics 1980; 31: 7398. 32. Kikuchi N. An Analysis of the variati onal Inequalities of S eepage Flow by FiniteElement Methods. Quartitative App lied Mathematics 1977; 35: 149163. 33. Brezis H, Kinderlehrer D, Stampacchia G, Sur une nouvelle formulation du probleme de lecoulement a travers une digue, Serie A. paris: C.R.Academie des Sciences; 1978. PAGE 136 ABOUT THE AUTHOR Kalyani Jeyisanker received her B.Sc. (Eng) with honors, in Civil Engineering from University of Peradeniya, Sri Lanka in 2001. She earned her Masters degree in Mechanical Engineering from Texas Tech University in 2004 and enrolled in PhD program in 2005 at University of South Florid a. While pursuing the PhD, Ms. Jeyisanker actively participated in numerous research projects. She also coauthored two journal papers and two conference papers. xml version 1.0 encoding UTF8 standalone no record xmlns http:www.loc.govMARC21slim xmlns:xsi http:www.w3.org2001XMLSchemainstance xsi:schemaLocation http:www.loc.govstandardsmarcxmlschemaMARC21slim.xsd leader nam Ka controlfield tag 001 002007176 003 fts 005 20090617121943.0 006 med 007 cr mnuuuuuu 008 090617s2008 flu s 000 0 eng d datafield ind1 8 ind2 024 subfield code a E14SFE0002783 035 (OCoLC)401711218 040 FHM c FHM 049 FHMM 090 TA170 (Online) 1 100 Jeyisanker, Kalyani. 0 245 Analysis of water seepage through earthen structures using particulate approach h [electronic resource] / by Kalyani Jeyisanker. 260 [Tampa, Fla] : b University of South Florida, 2008. 500 Title from PDF of title page. Document formatted into pages; contains 121 pages. Includes vita. 502 Dissertation (Ph.D.)University of South Florida, 2008. 504 Includes bibliographical references. 516 Text (Electronic dissertation) in PDF format. 520 ABSTRACT: A particulate model is developed to analyze the effects of steady state and transient seepage of water through a randomlypacked coarsegrained soil as an improvement to conventional seepage analysis based on continuum models. In the new model the soil skeleton and pore water are volumetrically coupled. In the first phase of the study, the concept of relative density has been used to define different compaction levels of the soil layers of a completely saturated pavement filter system and observe the seepage response to compaction. First, MonteCarlo simulation is used to randomly pack discrete spherical particles from a specified Particle Size Distribution (PSD) to achieve a desired relative density based on the theoretical minimum and maximum void ratios. Then, a water pressure gradient is applied across one twolayer filter unit to trigger water seepage.The pore water motion is idealized using Navier Stokes (NS) equations which also incorporate drag forces acting between the water and soil particles. The NS equations are discretized using finite differences and applied to discrete elements in a staggered, structured grid. The model predicted hydraulic conductivities are validated using widely used equations. The critical water velocities, hydraulic gradients and flow within the xi saturated soil layers are identified under both steady state and transient conditions. Significantly critical transient conditions seem to develop. In the second phase of the study the model is extended to analyze the confined flow through a partly saturated pavement layer and unconfined flow from a retention pond into the surrounding saturated granular soil medium. In partly saturated soil, the water porosity changes resulting from water flow is updated using the Soil Water Characteristics Curve (SWCC) of the soil.The results show how complete saturation develops due to water flow following the water porosity Vs pressure trend defined by the SWCC. Finally, the model is used to predict the gradual reduction in the water level of a retention pond and the location of the freesurface. The freesurface is determined by differentiating the wet and dry zones based on the Heaviside step function modified NS equations. 538 Mode of access: World Wide Web. System requirements: World Wide Web browser and PDF reader. 590 Advisor: Gunaratne, Manjriker Ph.D. 653 Navier Stokes equations MonteCarlo simulation Volumetric compatibility Staggered grid Critical localized condition 690 Dissertations, Academic z USF x Civil and Environmental Engineering Doctoral. 773 t USF Electronic Theses and Dissertations. 4 856 u http://digital.lib.usf.edu/?e14.2783 