Numerical analysis of full-scale tests on rammed aggregate piers

The purpose of this research was to evaluate the uplift capacity of single Rammed Aggregate Piers (RAPs) using numerical analysis and compare the results to the measured values from Full-Scale Tests (FSTs) conducted in the field. Each single RAP, uplift plate, and surrounding matrix soil was analyzed using a 3-dimensional Finite Difference Method (FDM). An important aspect of the RAP numerical modeling is construction sequence, which was simulated with three different techniques: Lift Compaction, Cavity Expansion, and the newly proposed Cavity Expansion with Variable Displacement technique. The constitutive model parameters were estimated from in situ tests and material properties reported in the literature. Two constitutive models were used in the numerical analyses - an elastic-perfectly plastic Mohr-Coulomb model and a Plastic-Hardening model that employs nonlinear stress-strain behavior in the elastic and plastic ranges. The RAP numerical analyses were performed with FLAC3D (Fast Lagrangian Analysis of Continua 3D) models, and the numerical models were validated against the result of five full-scale RAPs uplift tests that were performed in May and June of 1998 underneath the I-15 bridges near South Temple Street in Salt Lake City, Utah. The five tested RAPs had the same diameter of 3 ft (0.91 m) and heights of 3, 6, 9, 12, and 15 ft (0.91, 1.83, 2.74, 3.66, and 4.57 m). The FLAC3D models predicted the limiting uplift capacity of the RAPs measured during the FSTs with a 1% to 40% error depending on the RAP construction technique and the RAP length. The average accuracy of the three iv construction methods is nearly the same for each RAP length. The numerical models showed that the numerically expensive RAP construction process could be replaced with a horizontal stress initialization with insignificant loss of accuracy. The numerical model successfully captured the slight bulging outward of the lower portion of the RAP and formation of a short failure wedge near the surface at the maximum uplift force that were observed in the FSTs. The FLAC3D model with the Plastic-Hardening constitutive model could not predict the load-displacement of the FSTs with the same level of accuracy of the FLAC3D model with Mohr-Coulomb constitutive model. The lower accuracy of the more sophisticated Plastic-Hardening model is likely due to the estimation of the numerous parameters of this constitutive model with average material properties reported in the literature.

NUMERICAL ANALYSIS OF FULL-SCALE TESTS ON RAMMED AGGREGATE PIERS by Emad Ghodrati A thesis submitted to the faculty of The University of Utah in partial fulfillment of the requirements for the degree of Master of Science Department of Civil and Environmental Engineering The University of Utah December 2019 Copyright © Emad Ghodrati 2019 All Rights Reserved The University of Utah Graduate School STATEMENT OF THESIS APPROVAL The thesis of Emad Ghodrati has been approved by the following supervisory committee members: Evert C. Lawton , Chair 10/24/2019 Date Approved Steven F. Bartlett , Member 10/28/2019 Date Approved Chris P. Pantelides , Member 10/24/2019 Date Approved and by Michael E. Barber the Department/College/School of , Chair/Dean of Civil and Environmental Engineering and by David B. Kieda, Dean of The Graduate School. ABSTRACT The purpose of this research was to evaluate the uplift capacity of single Rammed Aggregate Piers (RAPs) using numerical analysis and compare the results to the measured values from Full-Scale Tests (FSTs) conducted in the field. Each single RAP, uplift plate, and surrounding matrix soil was analyzed using a 3-dimensional Finite Difference Method (FDM). An important aspect of the RAP numerical modeling is construction sequence, which was simulated with three different techniques: Lift Compaction, Cavity Expansion, and the newly proposed Cavity Expansion with Variable Displacement technique. The constitutive model parameters were estimated from in situ tests and material properties reported in the literature. Two constitutive models were used in the numerical analyses – an elastic-perfectly plastic Mohr-Coulomb model and a Plastic-Hardening model that employs nonlinear stress-strain behavior in the elastic and plastic ranges. The RAP numerical analyses were performed with FLAC3D (Fast Lagrangian Analysis of Continua 3D) models, and the numerical models were validated against the result of five full-scale RAPs uplift tests that were performed in May and June of 1998 underneath the I-15 bridges near South Temple Street in Salt Lake City, Utah. The five tested RAPs had the same diameter of 3 ft (0.91 m) and heights of 3, 6, 9, 12, and 15 ft (0.91, 1.83, 2.74, 3.66, and 4.57 m). The FLAC3D models predicted the limiting uplift capacity of the RAPs measured during the FSTs with a 1% to 40% error depending on the RAP construction technique and the RAP length. The average accuracy of the three construction methods is nearly the same for each RAP length. The numerical models showed that the numerically expensive RAP construction process could be replaced with a horizontal stress initialization with insignificant loss of accuracy. The numerical model successfully captured the slight bulging outward of the lower portion of the RAP and formation of a short failure wedge near the surface at the maximum uplift force that were observed in the FSTs. The FLAC3D model with the Plastic-Hardening constitutive model could not predict the load-displacement of the FSTs with the same level of accuracy of the FLAC3D model with Mohr-Coulomb constitutive model. The lower accuracy of the more sophisticated Plastic-Hardening model is likely due to the estimation of the numerous parameters of this constitutive model with average material properties reported in the literature. iv To my loving wife, Leila, and my adorable baby girl, Liana TABLE OF CONTENTS ABSTRACT....................................................................................................................... iii LIST OF TABLES ........................................................................................................... viii LIST OF FIGURES ........................................................................................................... ix ACKNOWLEDGMENTS ............................................................................................... xiv Chapters 1. INTRODUCTION .......................................................................................................... 1 2. LITERATURE REVIEW ............................................................................................... 4 2.1 Introduction .......................................................................................................... 4 2.2 Experimental Tests and Case Histories................................................................ 5 2.3 Numerical Modeling and Analysis .................................................................... 10 3. FULL-SCALE TESTS AT SOUTH TEMPLE 1998 ................................................... 16 3.1 Introduction ........................................................................................................ 16 3.2 Soil Profile ......................................................................................................... 16 3.3 Full-Scale Tests .................................................................................................. 17 4. NUMERICAL MODELING AND ANALYSIS .......................................................... 26 4.1 Introduction ........................................................................................................ 26 4.2 Model Geometry ................................................................................................ 30 4.3 Constitutive Model............................................................................................. 33 4.3.1 Structural Elements ................................................................................ 33 4.3.2 Soil Zones .............................................................................................. 38 4.4 Rammed Aggregate Piers Construction ............................................................. 59 4.4.1 K0-Blade Test Stress Distribution .......................................................... 59 4.4.2 Lift Compaction ..................................................................................... 59 4.4.3 Cavity Expansion ................................................................................... 60 4.4.4 Cavity Expansion with Variable Displacement ..................................... 62 4.5 Uplift Test .......................................................................................................... 71 4.5 Verification Models ........................................................................................... 72 5. RESULTS ..................................................................................................................... 96 5.1 FLAC3D Verification Models ........................................................................... 96 5.2 FLAC3D Full-Scale Tests Models..................................................................... 97 5.2.1 FLAC3D Models with Mohr-Coulomb Constitutive Model ................. 97 5.2.2 FLAC3D Models with Plastic-Hardening Constitutive Model ........... 108 6. CONCLUSION ........................................................................................................... 153 6.1 Concluding Remarks ........................................................................................ 153 6.2 Recommendations ............................................................................................ 155 Appendices A. LIST OF NOTATIONS ............................................................................................. 157 B. FLAC3D CODE ......................................................................................................... 164 REFERENCES ............................................................................................................... 195 vii LIST OF TABLES Tables 3.1. Soil classification and BSTs result (Modified from Hsu 2000) ................................. 19 3.2. Summary of FSTs on RAPs (Hsu 2000; Lawton 2000) ............................................ 19 4.1. Model size and grid statistic ...................................................................................... 73 4.2. The typical range for Poisson ratio, ν (Modifed from Bowles 1996) ........................ 73 4.3. The typical range for γt70 (Modifed from Itasca Consulting Group 2019) ................ 74 5.1. The uplift capacity compression for different RAP construction techniques .......... 113 LIST OF FIGURES Figures 2.1. Construction sequence of compressive RAPs (Modified from Fox and Cowell 1998). a. Drilling cavity, b. Aggregate placement at the bottom, c. High-frequency, low amplitude compaction of the first lift to create the bulb at the bottom, and d. introducing and compaction of aggregate in 12 in. (0.3 m) lift to construct the RAP ................................ 14 2.2. Uplift RAPs element (Modified from Wissmann and FitzPatrick 2016) .................. 15 3.1. Representative CPT sounding near the location of the FSTs (Modified from Lawton 2000) ................................................................................................................................. 20 3.2. Typical BST result (Modifed from Hsu 2000) .......................................................... 21 3.3. K0-Blade test results (Modified from Hsu 2000) ....................................................... 22 3.4. The uplift plate and four threaded bars (Modifed from Lawton 2000)...................... 23 3.5. Uplift setup on a RAP (Modified from Lawton 2000) .............................................. 24 3.6. Typical force-displacement results from the full-scale uplift test (Modifed from Hsu 2000). a. 3 ft (0.91 m) RAP, and b. 15 ft (4.57 m) RAP................................................... 25 4.1. Flowchart for the RAP uplift simulation in FLAC3D ............................................... 75 4.2. Flowchart for the RAP construction in FLAC3D, K0-Blade Stress technique .......... 76 4.3. Flowchart for the RAP construction in FLAC3D, Lift Compaction technique ......... 77 4.4. Flowchart for the RAP construction in FLAC3D, Cavity Expansion technique ....... 78 4.5. Flowchart for the RAP construction in FLAC3D, Cavity Expansion with Variable Displacement technique .................................................................................................... 79 4.6. Schematic RAP model in FLAC3D. (a) Uplift Plate (b) Elevation view of the RAP model in FLAC3D ............................................................................................................ 80 4.7. Zone generation by extruding an unstructured grid (6 ft RAP model). a. Plan view of the unstructured grid, and b. Elevation view of the grid ................................................... 81 4.8. FLAC3D model of the 6 ft RAP ................................................................................ 82 4.9. Uplift plate modeling. a. Triangle liner elements of the uplift plate, and b. Duplicated grid points at the bottom of the RAP ................................................................................ 83 4.10. Uplift plate and threaded rod modeling. a. Connection of the cable element to the uplift plate, and b. Embedded uplift plate and threaded rod (6 ft RAP model) ................ 84 4.11. FLAC3D liner element. a. The coordinate system and 6 DOFs of the liner element, b. Normal stress versus relative normal displacement, c. shear stress versus relative shear displacement, and d. Shear-strength criterion (Modifed from Itasca Consulting Group 2019) ................................................................................................................................. 85 4.12. FLAC3D cable element. a. The coordinate system and 1 DOF of the cable element, b. Cable element elastic, perfectly plastic material behavior, c. Grout material shear force per unit length versus relative shear displacement, and d. Grout material shear-strength criterion (Modifed from Itasca Consulting Group 2019) .................................................. 86 4.13. Native soil saturated unit weight and total vertical stress. a. Saturated unit weight (γsat ), and b. Total vertical stress (σv0 ) ........................................................................... 87 4.14. Native soil elastic saturated unit weight and total vertical stress. a. OCR, and b. Initial effective vertical stress (σ'v0 ) ............................................................................................ 88 4.15. The coefficient of lateral earth pressure at rest (K 0 ). a. Comparison of K 0 estimated by Equations 4.12 and 4.19 with measured values at the field, and b. Selected K 0 to initialize the horizontal stress in the FLAC3D model ...................................................... 89 4.16. Native soil elastic properties for the MC model. a. Poisson ratio (ν), and b. Young’s modulus (Es,OC ) ................................................................................................................ 90 4.17. The bulk (K) and shear (G) modulus of native soil ................................................. 91 4.18. The effective friction angle (ϕ'). a. Comparison of ϕ' estimated by Equations 4.14 and 4.17 with measured BST values at the field, and b. Selected ϕ' for FLAC3D modeling ........................................................................................................................................... 92 4.19. Stress-dependency of the E50 and Eur in the PH model. a. relatively low soil strength c cotϕpref = 0.5, and b. relatively high soil strength c cotϕpref = 3.0 .......................... 93 4.20. Shear modulus degradation curves. a. Normalized shear modulus versus normalized shear strain, and b. Normalized tangent shear modulus curve with recommended range for γt70 ..................................................................................................................................... 94 4.21. The process to find the Jacobian matrix components in cavity expansion with variable displacement technique ..................................................................................................... 95 5.1. Load-displacement curves for the anchor verification problem. a. 3 ft anchor, and b. 9 ft anchor .......................................................................................................................... 114 x 5.2. Coefficient of lateral earth pressure after FLAC3D RAP construction process of 3 ft RAP. a. Lift compaction, and b. Cavity Expansion ........................................................ 115 5.3. Coefficient of lateral earth pressure after FLAC3D RAP construction process of 6 ft RAP. a. Lift compaction, and b. Cavity Expansion ........................................................ 116 5.4. Coefficient of lateral earth pressure after FLAC3D RAP construction process of 9 ft RAP. a. Lift compaction, and b. Cavity Expansion ........................................................ 117 5.5. Coefficient of lateral earth pressure after FLAC3D RAP construction process of 12 ft RAP. a. Lift compaction, and b. Cavity Expansion ........................................................ 118 5.6. Coefficient of lateral earth pressure after FLAC3D RAP construction process of 15 ft RAP. a. Lift compaction, and b. Cavity Expansion ........................................................ 119 5.7. Increase in effective lateral stress in few steps to build the 3 ft RAP by Cavity expansion with variable displacement technique. a. Displacement pattern, b. Effective lateral stress, c. Displacement magnitude, and d. Vector alignment .............................. 120 5.8. Increase in effective lateral stress in a few steps to build the 6 ft RAP by Cavity expansion with variable displacement technique. a. Displacement pattern, b. Effective lateral stress, c. Displacement magnitude, and d. Vector alignment .............................. 121 5.9. Increase in effective lateral stress in a few steps to build the 9 ft RAP by Cavity expansion with variable displacement technique. a. Displacement pattern, b. Effective lateral stress, c. Displacement magnitude, and d. Vector alignment .............................. 122 5.10. Increase in effective lateral stress in a few steps to build the 12 ft RAP by Cavity expansion with variable displacement technique. a. Displacement pattern, b. Effective lateral stress, c. Displacement magnitude, and d. Vector alignment .............................. 123 5.11. Increase in effective lateral stress in a few steps to build the 15 ft RAP by Cavity expansion with variable displacement technique. a. Displacement pattern, b. Effective lateral stress, c. Displacement magnitude, and d. Vector alignment .............................. 124 5.12. Load-displacement curves of 3 ft RAP uplift test .................................................. 125 5.13. Load-displacement curves for 6ft RAP uplift test ................................................. 126 5.14. Load-displacement curves for 9 ft RAP uplift test ................................................ 127 5.15. Load-displacement curves for 12 ft RAP uplift test .............................................. 128 5.16. Load-displacement curves for 15 ft RAP uplift test .............................................. 129 5.17. Maximum shear strain rate during pull-out with the applied uplift load higher than the uplift capacity. a. 3 ft RAP, and b. 9 ft RAP ............................................................. 130 5.18. The uplift capacity compression for different RAP construction techniques ........ 131 xi 5.19. Change in lateral stress and lateral displacement at the matrix soil close to the 3 ft RAP during the FLAC3D uplift test for RAP construction. a., b., and c. change in lateral stress, and d., e., and f. change in interface lateral displacement between the threaded rods (θ = 0°) in front of the threaded rod (θ = 45°).............................................................. 132 5.20. Change in lateral stress and lateral displacement at the matrix soil close to the 6 ft RAP during the FLAC3D uplift test for RAP construction. a., b., and c. change in lateral stress, and d., e., and f. change in interface lateral displacement between the threaded rods (θ = 0°) in front of the threaded rod (θ = 45°).............................................................. 133 5.21. Change in lateral stress and lateral displacement at the matrix soil close to the 9 ft RAP during the FLAC3D uplift test for RAP construction. a., b., and c. change in lateral stress, and d., e., and f. change in interface lateral displacement between the threaded rods (θ = 0°) in front of the threaded rod (θ = 45°).............................................................. 134 5.22. Change in lateral stress and lateral displacement at the matrix soil close to the 12 ft RAP during the FLAC3D uplift test for RAP construction. a., b., and c. change in lateral stress, and d., e., and f. change in interface lateral displacement between the threaded rods (θ = 0°) in front of the threaded rod (θ = 45°).............................................................. 135 5.23. Change in lateral stress and lateral displacement at the matrix soil close to the 15 ft RAP during the FLAC3D uplift test for RAP construction. a., b., and c. change in lateral stress, and d., e., and f. change in interface lateral displacement between the threaded rods (θ = 0°) in front of the threaded rod (θ = 45°).............................................................. 136 5.24. The FLAC3D models' vertical force distribution through the RAP length at the maximum uplift load. a. lift compaction, b. cavity expansion, and c. variable displacement ......................................................................................................................................... 137 5.25. Threaded bar elongation and maximum deformation of the uplift plate (differential vertical displacement of the node at the center of the plate and the node connected to the threaded bar) during the FLAC3D uplift test of 3 ft RAP. a. Lift compaction, b. Cavity expansion, and c. Cavity expansion with variable displacement .................................... 138 5.26. Threaded bar elongation and maximum deformation of the uplift plate (differential vertical displacement of the node at the center of the plate and the node connected to the threaded bar) during the FLAC3D uplift test of 6 ft RAP. a. Lift compaction, b. Cavity expansion, and c. Cavity expansion with variable displacement .................................... 139 5.27. Threaded bar elongation and maximum deformation of the uplift plate (differential vertical displacement of the node at the center of the plate and the node connected to the threaded bar) during the FLAC3D uplift test of 9 ft RAP. a. Lift compaction, b. Cavity expansion, and c. Cavity expansion with variable displacement .................................... 140 5.28. Threaded bar elongation and maximum deformation of the uplift plate (differential vertical displacement of the node at the center of the plate and the node connected to the threaded bar) during the FLAC3D uplift test of 12 ft RAP. a. Lift compaction, b. Cavity expansion, and c. Cavity expansion with variable displacement .................................... 141 xii 5.29. Threaded bar elongation and maximum deformation of the uplift plate (differential vertical displacement of the node at the center of the plate and the node connected to the threaded bar) during the FLAC3D uplift test of 15 ft RAP. a. Lift compaction, b. Cavity expansion, and c. Cavity expansion with variable displacement .................................... 142 5.30. Structural elements force and stress at the maximum uplift load for FLAC3D 3 ft RAP model with variable displacement RAP construction. a. cable force, and b. liner maximum fiber stress ...................................................................................................... 143 5.31. Structural elements force and stress at the maximum uplift load for FLAC3D 6 ft RAP model with variable displacement RAP construction. a. cable force, and b. liner maximum fiber stress ...................................................................................................... 144 5.32. Structural elements force and stress at the maximum uplift load for FLAC3D 9 ft RAP model with variable displacement RAP construction. a. cable force, and b. liner maximum fiber stress ...................................................................................................... 145 5.33. Structural elements force and stress at the maximum uplift load for FLAC3D 12 ft RAP model with variable displacement RAP construction. a. cable force, and b. liner maximum fiber stress ...................................................................................................... 146 5.34. Structural elements force and stress at the maximum uplift load for FLAC3D 15 ft RAP model with variable displacement RAP construction. a. cable force, and b. liner maximum fiber stress ...................................................................................................... 147 5.35. Increase in effective lateral stress in a few steps to build the 9 ft RAP by Cavity expansion with variable displacement technique (PH Model). a. Displacement pattern, b. Effective lateral stress, c. Displacement magnitude, and d. Vector alignment ............... 148 5.36. Load-displacement curves for 9 ft RAP uplift test (PH Model) ............................ 149 5.37. Change in lateral stress and lateral displacement at the matrix soil close to the 9 ft RAP during the FLAC3D uplift test for RAP construction (PH Model). a., b., and c. change in lateral stress, and d., e., and f. change in interface lateral displacement between the threaded rods (θ = 0°) in front of the threaded rod (θ = 45°) ................................. 150 5.38. Threaded bar elongation and maximum deformation of the uplift plate (differential vertical displacement of the node at the center of the plate and the node connected to threaded bar) during the FLAC3D uplift test of 9 ft RAP (PH Model). a. Lift compaction, b. Cavity expansion, and c. Cavity expansion with variable displacement .................... 151 5.39. Structural elements force and stress at the maximum uplift load for FLAC3D 9 ft RAP model with variable displacement RAP construction (PH Model). a. cable force, and b. liner maximum fiber stress.......................................................................................... 152 xiii ACKNOWLEDGMENTS I would like to thank my advisor at the University of Utah, Dr. Lawton, and my mentor at the ITASCA Consulting Group Inc., Mr. Augusto Lucarelli, for their kind support and valuable guidance throughout my education at the University of Utah. Also, I would like to thank the ITASCA Consulting Group Inc. for providing us with the FLAC3D software through the ITASCA Education Partnership (IEP) research program. CHAPTER 1 INTRODUCTION Rammed Aggregate Piers (RAPs) were introduced as an effective and economical soil improvement method in the 1980s. The RAPs system has been successfully used in the settlement control, bearing capacity improvement, liquefaction mitigation, slope stabilization, and uplift enhancement. The unique construction procedure of the RAPs is the main reason for its overall high performance, efficiency, and application as an economical alternative to traditional soil improvement techniques such as stone columns, over-excavation and replacement, deep compaction, chemical stabilization, and pile foundations. A RAP is constructed by drilling a borehole, removing the native soil and building the RAP by placing and compacting small lifts of granular material. The special beveled head compactor exerts high compaction energy to each lift by vibrating in high frequency and pushing downward and outward. As a result of this high energy, highfrequency compaction, the RAP material is densified to a significant degree, and large residual horizontal stresses are generated within the adjacent matrix soil, which further stiffens and strengthens the RAPs system by providing extra confining pressure. Researchers have performed several experimental tests on RAPs to better understand the effect of RAPs construction on the properties of the matrix soil, the residual horizontal and vertical stress in the matrix soil and within the RAPs, compressive bearing 2 capacity and uplift capacity of the RAPs, cavity expansion, and stress path changes due to RAP construction and loading. Based on the results of these experimental tests, a few numerical and simple theoretical methods have been developed to estimate the compressive bearing capacity and uplift capacity of single and groups of RAP. However, the complex nature of the RAPs system and limitation in experimental procedure and budget hinders us from getting the full picture of their behavior. To overcome this difficulty, different numerical models of the RAP have been built and analyzed. Even though the numerical analysis could provide a great insight into the problem, its validity and accuracy are highly dependent on the numerical analysis method, the refinement in the model, restriction of the constitutive models, steps in numerical model for construction consideration, and the accuracy of the test data that provide the material properties. Also, the field and laboratory tests to estimate different soil properties, which are used in numerical modeling, are usually limited in number and the soil property that they represent. As a result, the constitutive modeling is usually limited to simplistic analyses. More sophisticated constitutive modeling is possible, but more assumptions are required to parameterize the whole model. On the other hand, more sophisticated numerical models, with a high degree of refinement, are usually computationally expensive, and they have a higher chance of numerical instability. Therefore, numerical modeling requires a careful balance for model refinement and sophistication to model the physical phenomenon with acceptable accuracy, computational cost and soil test cost required to parameterize the model. In this research, the results of five full-scale RAPs uplift tests that were performed in May and June of 1998 at the South Temple Street in Salt Lake City, Utah, are modeled with FLAC3D (Fast Lagrangian Analysis of Continua 3D). The experimental tests were 3 led by the Department of Civil and Environmental Engineering at the University of Utah. The five tested RAPs had the same diameter of 3 ft (0.91 m) and heights of 3, 6, 9, 12, and 15 ft (0.91, 1.83, 2.74, 3.66, and 4.57 m) (Hsu 2000; Lawton 2000). The results of the Cone Penetration Tests (CPT), Borehole Shear Tests (BST), K0-Blade Tests, and laboratory tests of the matrix soil and the RAPs material reported by Hsu (2000) are used to parameterize the constitutive model of numerical analysis. Moreover, the numerical model was validated against the results of the uplift tests. The single RAP and surrounding matrix soil were modeled in the FLAC3D and analyzed by the Finite Difference Method (FDM). The FDM methods discretize the body into small zones within a grid of nodes and solve the kinematic equations and constitutive model on these nodes and zone. As a result, compared to the Finite Element Method (FEM), FDM can analyze the model much faster because it does not need to construct the stiffness, mass, and damping matrices. In addition, FDM could be used more efficiently for large deformation analysis compare to FEM. In this thesis, Chapter 2 is the literature survey and discusses the RAP concept and reviews the state of the art theoretical, empirical, and numerical methods to estimate the uplift capacity. Also, it provides a review of the numerical modeling and analysis performed by several researchers — Chapter 3 overviews the full-scale test (FST) procedure and result from the South Temple 1998 research project. In Chapter 4, the modeling of the RAP, parametrization, and analysis steps in FLAC3D are presented. The results are discussed in Chapter 5. Finally, concluding remarks and recommendations for future research are summarized in Chapter 6. CHAPTER 2 LITERATURE REVIEW 2.1 Introduction The RAPs were introduced as an effective and economical soil improvement method in the 1980s. Although It was initially developed as an economical alternative to overexcavation and replacement (Lawton et al. 1994; Lawton and Fox 1994) its efficacy and cost-effectiveness have grown its application to settlement control, bearing capacity improvement, liquefaction mitigation, slope stabilization, and uplift enhancement (Kwong et al. 2002; Lawton et al. 1994; Lawton and Fox 1994; Rudolph et al. 2011; White et al. 2001). The special construction process of RAP increases the confining pressure in the matrix soil and the RAP element by displacing the matrix soil outward, which enhances the stiffness and strength of the soil-RAPs composition significantly. The construction of RAPs is well documented in the literature (Lawton et al. 1994; Lawton and Fox 1994). The steps to construct a compressive RAPs element are as follows (see Figure 2.1): 1. A cavity is created in the soil. The cavity diameter is usually between 24 in. to 36 in. (0.61 m to 0.91 m). And, the height of the cavity usually ranges from 7 ft to 30 ft (2 m to 9 m). 2. The soil at the bottom of the cavity is compacted using a special beveled head 5 tamper. A thin lift of well-graded aggregate, not exceeding 12 in. (0.3 m), is placed at the bottom of the cavity. This first lift is compacted using the beveled head tamper with high frequency and low amplitude vibration to exert high compaction energy to this layer. As a result, tamper densifies the RAP material and pushes the particle outward to expand the cavity in diameter and depth and create the “bulb.” The “bulb” formation would increase the horizontal and vertical stress in the adjacent soil. 3. The well-graded aggregate material is placed in lifts of 12 in (0.3 m) or less on top of the previous layer and gets compacted with the beveled head tamper. The process is repeated to build the RAPs element to the desired height. The construction of the uplift RAP follows a similar procedure with an extra step. After compacting the first lift and formation of the “bulb” a horizontal steel plate and attached vertical threaded steel bars are lowered in the cavity and placed on top of the first lift. Then the rest of the RAPs element is built by placing and compacting the 12 in. (0.3 m) Lifts similar to the compression RAPs. The threaded bars are extended through the RAPs element and are embedded in the concrete foundation to form an integrated foundation system (see Figure 2.2). 2.2 Experimental Tests and Case Histories Lawton and Fox (1994) proposed the division of the soil reinforced with RAPs to the upper zone (UZ) and the lower zone (LZ) for settlement calculation. The UZ, which represents the composite zone of RAP and matrix soil, extended from the ground surface to the bottom of the RAPs plus a RAP width to consider the effect densification of the soil 6 beneath the RAP. The settlement of the UZ zone is calculated by using a parallel spring analogy for matrix soil and RAPs assuming they have the same settlement, and the total force is the sum of the force on RAPs and matrix soil. The LZ consisted of the soil beneath the UZ zone. The settlement of the UZ zone was calculated by estimating the transferred stress from the UZ to the LZ, considering the effect of the higher stiffness of the UZ in distributing the force on a larger area and reducing the transferred stress. The authors summarized the application of RAPs to control the settlement for projects. The RAPs system successfully reduced the predicted unreinforced foundation settlement. Also, the measured settlements were in reasonable agreement with the predicted RAPs reinforced settlement calculated using the proposed method by authors. The authors discussed three case histories in more detail. In the first case history, RAPs were proposed as an economical alternative to auger-cast piles for a five-story office building at Columbia South Carolina. Static load tests on 4 ft (1.2 m) deep test RAPs with 2 ft (0.61 m) diameter showed a factor of 2 increase in the subgrade modulus. The installed RAPs with diameters of 2.5 ft and 3 ft (0.76m and 0.91 m) and heights of 5 ft and 6 ft (1.5 m and 1.8 m) reduced the predicted 1.3 in. to 4.0 in (33 mm to 102 mm) settlement of square footing on unreinforced soil to the predicted value of 0.7 in. (18 mm) and the measured value of less than 0.06 in (1.5 mm) for the RAPs reinforced system. The second case history discussed the construction of rectangular prismatic RAPs to reinforce the foundation of a 40 ft (12 m) high milk silo at Atlanta, Georgia. The RAPs are 2 ft (0.61 m) wide and 5.0 ft (1.5 m) deep. The predicted settlement was 1.9 in. to 4.1 in (48 mm to 104 mm) for unreinforced foundation and 0.5 in. (13 mm) for RAPs reinforce the foundation. However, the measured settlement after three months was less than 0.07 in. (1.8 mm). The third project was an industrial manufacturing 7 building at Winterset, Iowa. The 2 ft (0.61 m) and 2.5 ft (0.76 m) deep RAPs with an aspect ratio of two, reinforced the soil beneath the square footings of the one-story building. The result of the two static tests indicated the insensitivity of RAPs behavior to the construction condition and low water table. Also, the result of the K0-blade test and Borehole Shear Tests (BST) suggested that the limiting passive pressure was achieved during RAPs construction at the RAP-soil interface. The predicted bearing capacity increased by a factor of two as a result of RAPs construction and the predicted settlement of 6 in. to 9 in. (150 mm to 230 mm) for footings on unreinforced soil reduced to 0.9 in. (23 mm) for RAPs system, which was close to the measured settlement of less than 0.75 in. (19 mm). Lawton et al. (1994) described two case histories where RAPs were successfully used to control settlement and provide uplift capacity (Lawton et al. 1994). RAPs were constructed beneath both individual and mat foundations in the expansion project of the regional hospital in Atlanta, Georgia, to control the settlement of two towers. The prismatic rectangular RAPs were able to reduce the predicted settlement of 1.1 - 3.9 in. (29 - 98 mm) to less than 0.25 - 0.75 in. (6 to 20 mm) at a different location. The settlement of the mat and footings on the RAPs reinforced soil were calculated by the proposed technique by the authors in their previous publication (Lawton and Fox 1994). The second case history explained the application of RAPs in providing uplift capacity for the Mississippi air national guard hanger at Meridian, Mississippi. The considerable wind load on the airplane hangar necessitates an uplift resisting system. Prismatic rectangular RAPs were used in place of initially designed helical anchors to reduce cost and increase ease in construction. The RAPs performed almost elastically for the measured uplift force-displacement during the test. However, the RAPs were not tested to failure. The authors suggested that the 8 failure happened at the interface of the RAP and matrix soil. Moreover, they reasoned that the high horizontal stress that was built up during the special RAP construction increased the uplift capacity because of the arching of the principal stress at the RAP-soil interface. Also, they stated that the horizontal stress at the interface could decrease during the pullout test. However, they estimated the uplift capacity of the RAP by assuming the development of full passive Rankine horizontal pressure at the RAP-soil interface without any decrease at failure due to the lack of experimental data. Hsu (2000) performed a series of five uplift FSTs were on RAPs with a varying height from 3 to 15 ft (0.91 to 4.57 m). The soil properties were determined by a series of CPTs, BSTs, K0-blade testes, BSTs, and laboratory tests. The lateral stress increased significantly as a result of the RAP construction. And, the lateral stress at the about the top 6 ft (1.83 m) of the RAPs reached the passive Ranking stress. The RAP construction highly disturbed the clayey matrix soil, which resulted in complete loss of cohesion intercept due to the remolding of the soil. The author reported a decrease in the lateral stress of about 144 psf (6.9 kPa) at the maximum uplift load for short piers. Also, the inclinometers near the RAPs showed the formation of a truncated cone near the top portion of the RAP and bulging outward of the bottom portion of the RAP at the maximum uplift load. An analytical model was proposed to estimate the uplift capacity of the RAP by assuming a cylindrical shear failure surface at the perimeter of the RAP. The uplift capacity was estimated by calculating the shear stress along the RAP shaft with Mohr-Coulomb strength parameter of the matrix soil. The full passive lateral stress is assumed for the top 6 ft (1.83 m), and the average value of the coefficient of the lateral earth pressure for at-rest and passive conditions was used to calculate the lateral stress below 6 ft (1.83). To account for 9 the loss of lateral stress at the maximum uplift load, 144 psf (6.9 kPa) was subtracted from the calculated effective lateral stress. The proposed analytical method could predict the uplift capacity with less than 30% error. The author proposed a semiempirical method to estimate the uplift capacity based on the CPT tip resistance and sleeve friction. However, the CPT data were not enough to develop a statistically meaningful correlation. And, it was recommended to use the semiempirical method alongside the analytical method to provide a preliminary estimate of the uplift capacity (Hsu 2000; Lawton 2000). White et al. performed full-scale compression tests on two square footings supported by groups of four RAPs and three isolated RAPs with concrete caps (White et al. 2007). The soil profile consisted of 3.3 ft (1.0 m) of desiccated layer overlaying uniform soft alluvial clay stratum. All the RAPs had a diameter of 2.5 ft (0.76m). The group RAPs were 9.2 ft (2.8 m) and 16.7 ft (5.1 m) deep. The first two isolated RAPs had the same depth and diameter of the RAPs groups, and the third RAP was 9.75 ft (2.97 m) deep. The RAPs were constructed using a well-graded crushed limestone. The Consolidated Drain (CD) triaxial tests on the RAP material showed dilative behavior under deviatoric loading. Effective cohesion intercept (𝑐 ′ ) of 83.5 psf (4 kPa), friction angle (∅′ ) of 47º, dilatancy angle of 9º to 16 º with an average of 12 º, and secant elastic modulus of 1775 ksf (85 MPa) were determined for the RAPs material. The compression tests’ results showed that the difference in the displacement of the top of the RAPs and bottom of the RAPs (measured by telltale) are highest for longer RAPs, which suggested bulging for RAPs with higher length to diameter ratio and tip movement for RAPs with lower length to diameter ratio. The authors reported group effect (ratio of a total load of RAPs group divided by the number of RAPs to the isolated RAP load) of about 1.0 for the two groups of RAPs at 10 lower compression load. However, the authors claimed the group effect of 4.7 was achieved for shorter RAPs ground due to the contribution of the matrix soil in load-carrying capacity and initiation of bearing capacity failure under the reinforced zone. The stress concentration ratio of 4 to 5 was reported for the RAPs. 2.3 Numerical Modeling and Analysis LaMares (2005) used an axisymmetric FLAC model to simulate RAPs and drilled concrete piers for uplift and bearing capacity. The author proposed the simulation of the RAP construction process by creating a cavity and placing the RAP material in single lifts in the borehole. After placing each lift of the RAP, a specific pressure was applied at the top of the lift to compact it statically and increase the vertical and horizontal confining pressure in the matrix soil. After compressing the RAP lift for a certain number of steps in FLAC, the applied pressure was removed, and the next RAP lift was placed on top of the previous one. The process was repeated until the full height of the RAP was constructed in the FLAC model. The uplift load was applied by removing a layer of the soil beneath the RAP element and applying the uplift load at the bottom of the RAP element. LaMeres (2005) used the Mohr-Coulomb constitutive model for the soil and RAP material. The hyperbolic elastic model was evaluated for modeling soil, as well. The FLAC uplift capacity was within 5% of the FST maximum uplift load. The drilled concrete pier underpredicted the ultimate uplift capacity by 32%. The FLAC bearing capacity underpredicted the FST bearing capacity by 10% and 21% for the RAP and drilled concrete piers, respectively. The RAP model with hyperbolic constitutive mode and Mohr-Coulomb constitutive model resulted in similar displacement for uplift and compression tests. The 11 hyperbolic model was unstable near the ultimate uplift capacity (LaMeres 2005). Pham and White (2007) modeled a single RAP and a group of RAPs by an axisymmetric finite element model in PLAXIS (Pham and White 2007) and compared their results with the FSTs results in their previous paper (White et al. 2007). The group of RAPs was modeled by a unit cell which consisted of a single RAP and matrix soil within its tributary area. A concrete cap was modeled on top of the unit cell and the single RAP. The elements at the interface of the soil-cap modeled with 50% strength reduction. The soilhardening constitutive model (Schanz et al. 1999) was selected for matrix soil and RAP elements, which has stress-dependent stiffnesses and Mohr-Coulomb criterion with a cap for yield failure. The model was initialized for the K0 condition after the cavity was modeled. The RAPs construction was modeled by expanding the borehole cavity until the displacement of 5% and 10% of the nominal RAP diameter was achieved at the shaft and bottom of the borehole, respectively. The target values of displacement were based on the field measurement of their FSTs. The stress-controlled compression load was applied gradually on top of the concrete cap in a few steps. The load-displacement of the single and group RAPs for the bearing capacity simulation had an overall good agreement with the test result, especially at the lower load. Using the same process as Lawton and Hsu (Hsu 2000; Lawton 2000), the design bearing capacity predicted by PLAXIS modeling overestimated the measured values by 6% to 31%. The predicted lateral displacement was 65% of the measured values by inclinometers. The predicted stress concentration factor of 2 to 10 by the numerical modeling covered the range of 4 to 5 measured in the FSTs. The stress path analysis of at two depth of 2.46 ft (0.75 m) and 4.92 ft (1.5 m) indicated the soil at the interface with the RAP and group of RAPs unit cell reached Rankine passive stress 12 after the RAPs construction. The finite element results suggested that the vertical stress decreases linearly within the upper zone and within the lower zone. Therefore, vertical stress could be estimated with a bilinear curve to calculate the settlement. Halabian and Shamsabadi (2014) used 2D axisymmetric models to investigate the RAP construction process in RAP behavior. The authors used a Distinct Element Method model for the RAP material coupled with a continuum model for the matrix soil, which was modeled with a Finite Difference Method (FDM) and Mohr-Coulomb constitutive model. They validated their model with the FSTs on RAPs performed by White et al. (2007). The results of their numerical model of the FSTs and parametric studies showed that the RAP length and matrix soil friction angle and cohesion affect the cavity expansion, and a using a constant value for lateral cavity expansion based on a percentage of the RAP diameter could result in errors. The authors suggested a constant lateral displacement of 4.5% to 5% of RAP diameter within the top 4.9 ft (1.5 m) of the RAP, zero lateral displacements within the bottom 6.6 ft (2 m) of the RAP, and a linearly variable lateral displacement within a middle portion of the RAP. The ramming process changed the stress state within the radial distance of 3 to 5 times the RAP diameter in their numerical model, and the increase in the cohesion intercept enlarged the mechanically improved matrix soil zones around the RAP. Demir et al. (2017) performed a series of bearing capacity laboratory tests on RAPs installed in a sandy clay matrix soil inside a test tank and modeled one of the single RAPs with a 2D axisymmetric FEM model. Three single RAP and a group of three RAPs were constructed by pushing a column case downward to create cavity, placing well-graded crushed aggregate at the bottom of the cavity, and compacting the RAP lift. The RAPs had 13 a diameter of 3.9 in. (100 mm) and height of 23.6 in. (600 mm). The result of their preliminary analytical and numerical calculation suggested that the tank was small, and the boundary effect could have affected their results. The three single RAPs were constructed in soils with different undrained shear strength. The bearing capacity laboratory tests showed that the bearing capacity could increase significantly with a slight increase in the undrained shear strength. The single RAPs had a stress concentration ratio of 6 to 8, and the RAP group showed a stress concentration ratio of 2 to 4. The measured lateral stress showed an improved matrix soil zone of 4.5 times the diameter of the RAP. The authors used Plaxis with PH constitutive model for the 2D axisymmetric numerical modeling of one of the RAPs. The RAP construction was simulated by expanding the cavity laterally and vertically to 10% of the RAP diameter. The load-displacement curve of the numerical model and the laboratory test showed a good agreement. And, the numerical model could predict the bulging of the RAP at the top portion of the RAP at the bearing capacity. 14 a. b. c. d. Figure 2.1. Construction sequence of compressive RAPs (Modified from Fox and Cowell 1998). a. Drilling cavity, b. Aggregate placement at the bottom, c. High-frequency, low amplitude compaction of the first lift to create the bulb at the bottom, and d. introducing and compaction of aggregate in 12 in. (0.3 m) lift to construct the RAP 15 Uplift load Figure 2.2. Uplift RAPs element (Modified from Wissmann and FitzPatrick 2016) CHAPTER 3 FULL-SCALE TESTS AT SOUTH TEMPLE 1998 3.1 Introduction The full-scale RAPs uplift testing that is investigated in this thesis was performed under the old southbound I-15 bridge over South Temple in Salt Lake City, Utah, before its demolition during May and June of 1998. The five RAPs have nominal diameters of 3 ft (0.91 m) and heights of 3, 6, 9, 12, and 15 ft (0.91, 1.83, 2.74, 3.66, and 4.57 m)(Hsu 2000; Lawton 2000). A circular A-36 steel plates with a diameter of 34 in. (864 mm) and a thickness of 0.5 in. (12.7 mm) were placed at the bottom of the boreholes before placing and compacting RAPs lifts to provide the uplift capacity. The uplift plates were connected to four No. 7 or No.8 Grade 75 threaded bars. The bars were extended through the height of each RAP and were pulled out by a hydraulic jack during the uplift testing. The matrix soil and RAP material were characterized by performing sieve, hydrometer, and Atterberg limits tests in the laboratory on samples collected in the field and Cone Penetration Tests (CPT), Borehole Shear Tests (BST), K0-Blade tests and Sand cone tests in situ. 3.2 Soil Profile The soil profile at the research site consists of about 19 ft of alluvial deposits over the Bonneville clay. Sedimentation of the sand, silt, and clay of the streams coming from 17 the major canyons formed the alluvial stratum during the Holocene epoch. The Bonneville clay was deposited during the Pleistocene epoch about 25,000 to 12,000 years ago. It consists of clayey silt and silty clay, and some fine sand mainly interbedded in the middle of the Bonneville (Hsu 2000; Lawton 2000). A total of ten CPTs at ten locations and nineteen BSTs at one location were conducted to characterize the strength, stiffness, and stress history of the native soil and RAPs material. A representative CPT and BST results are shown in Figures 3.1 and 3.2, respectively. After adjusting the depth of CPTs to consider the elevation difference at each CPT location, the CPTs data showed good agreement and at all locations, which indicates approximately horizontal soil layers (Lawton 2000). The soil strength parameters are obtained from the BSTs results up to a depth of 15.5 ft (4.72 m) (see Table 3.1). However, FLAC3D simulation should model the soil to a much deeper depth to minimize the boundary effect. Therefore, the strength parameters calculated using CPT data complemented the BST data. The soil elastic stiffness parameters were calculated from the CPT sounding. The calculation of material properties is discussed in detail in Chapter 4. Lawton and Hsu performed K0-Blade tests to calculate the coefficient of lateral earth pressure at 2.5 ft (0.76 m), 5.0 ft (1.52 m) and far away from the center of a 15 ft (4.57 m) RAP element (Hsu 2000; Lawton 2000) (see Figure 3.3). The far filed K0-Blade test represents the native soil’s coefficient of lateral earth pressure at rest (𝐾0 ). 3.3 Full-Scale Tests The full-scale uplift tests were performed on each of the RAPs by pulling out the threaded bars which were connected to the uplift plate. A system of steel beams and two 18 compression RAPs were used to provide the required reaction for the hydraulic jack to pull out the thread bar. The 15 ft (4.57 m) long compression RAPs with 36 in. (0.91 m) the diameter was located at two sides of the uplift test RAP with the center to center span of 8 ft (2.44 m). Figure 3.4 shows the uplift plate and its threaded bars before installation, and Figure 3.5 shows the actual tests setup. The threaded bars were pulled out at approximately a constant rate during the uplift test. The FSTs were continued until the maximum load has been achieved for the 3 ft (0.91 m) long RAP or the top of the soil reached the bottom of the reaction steel beams for the 6, 9, and 15 ft (1.83, 2.74, and 4.57 m) RAPs. During the test on 12 ft (3.66 m) RAP, one of the nuts at the bottom of the bar came loose. Restarting the test after rearranging the loading system to load two opposite bars resulted in excessive bending in one of the reaction beams. The test was stopped at that point without reaching the peak failure uplift load. The force-displacement results for the 3 ft (0.91 m) and 15 ft (4.57 m) long RAPs are shown in Figure 3.6. The design uplift force (Td) is determined by drawing a straight line tangent the initial flat slope portion of the curve and finding its intersection with a second straight line tangent to the steep slope at the end of the curve where a small increase in uplift force results in a significant increase in displacement (Hsu 2000; Lawton 2000). The design and maximum uplift force and the displacement of the bars corresponding to the design uplift force (Td) are summarized in Table 3.2. In addition to the force-displacement, changes in horizontal stress were measured during FSTs by placing horizontal pressure cells at depth 1 ft (0.30 m) and 5 ft (1.53 m) close to the interface of RAPs and matrix soil of 3 ft (0.91 m) and 6 ft (1.83 m) long RAPs respectively. Hsu reported the induced vertical compression stress within the 6 ft (1.83 m) and 9 ft (2.74 m) tall during FSTs (Hsu 2000). 19 Table 3.1. Soil classification and BSTs result (Modified from Hsu 2000) Note: The groundwater table was 4 ft (1.2 m) below the ground surface. USCS = united soil classification system, NP = Nonplastic, LL = liquid limit, PL = plastic limit, and PI = plasticity index Table 3.2. Summary of FSTs on RAPs (Hsu 2000; Lawton 2000) 20 Tip Resistance, qc (ksf) 0 100 200 Friction ratio, fR (%) Sleeve friction, fs (ksf) 0 2 4 0 2 4 6 0 Adjusted depth below ground surface, z (ft) -5 -10 -15 -20 -25 -30 Figure 3.1. Representative CPT sounding near the location of the FSTs (Modified from Lawton 2000) 21 1.5 τ (ksf) BST C1 Z = 3.0 ft CL, Lean Clay ∅′𝑁𝐶 = 38.6° ′ 𝑐𝑁𝐶 =0 1.0 0.5 ∅′𝑂𝐶 = 14.0° ′ 𝑐𝑂𝐶 = 0.298 𝑘𝑠𝑓 0.0 0.0 0.5 1.0 1.5 2.0 σ' (ksf) Figure 3.2. Typical BST result (Modifed from Hsu 2000) 2.5 22 Coefficient of lateral earth pressure, K 0 1 2 3 0 -2 Depth below ground surface, z (ft) -4 -6 -8 -10 K at far field K at r = 2.5 ft K at r = 5.0 ft -12 Figure 3.3. K0-Blade test results (Modified from Hsu 2000) 4 23 Uplift Plate Figure 3.4. The uplift plate and four threaded bars (Modifed from Lawton 2000) 24 Threaded bars RAP Figure 3.5. Uplift setup on a RAP (Modified from Lawton 2000) 25 Displacement of uplift bars, δv (in.) 0.8 0.6 0.4 0.2 0.0 0 Td 10 15 Uplift force, T (kips) 5 20 a. Displacement of uplift bars, δv (in.) 3.5 3.0 2.5 2.0 1.5 1.0 0.5 0.0 0 20 40 60 80 100 120 Uplift force, T (kips) Td 140 160 180 b. Figure 3.6. Typical force-displacement results from the full-scale uplift test (Modifed from Hsu 2000). a. 3 ft (0.91 m) RAP, and b. 15 ft (4.57 m) RAP CHAPTER 4 NUMERICAL MODELING AND ANALYSIS 4.1 Introduction The five FST’s of the South Temple 1998 research project are modeled and analyzed using the FLAC3D 7.0 (Itasca Consulting Group 2019). Only one-quarter of the actual FSTs are modeled due to the symmetry of the problem. The matrix soil and the RAP continua are discretized into hexahedron zones using an unstructured mesh. In the FLAC3D, the structural elements were implemented with FEM elements and nodes, and the soil continuum media were formulated with the FDM zones and grid points. The structural FEM nodes interact with the FDM zones (not grid points) through links that transfer force and velocity between the two entities based on the interpolation location of a node within a zone. The structural element links which connect the structural nodes to zones can be free, rigid, and deformable. In this research, the uplift plate is modeled with two-dimensional liner structural elements embedded beneath the bottom of RAP. Six degrees of freedom (three translational and three rotational) at each node of the triangle liner elements are activated to model the full shell-membrane behavior of the uplift plate. The embedded liner elements interact with zones on both sides with the user-defined interface properties. They have linear stressdisplacement spring behavior in the normal direction to the liner surface with a cut-off 27 tension strength to initialize separation. In the shear direction, they have linear stressdisplacement spring behavior with the Mohr-Coulomb failure criterion. The threaded rod is modeled with one-dimensional cable structural elements. The cable element has one degree of freedom at each node for the axial displacement. In addition to truss-like behavior, the cable elements have a grout interface along their length to interact with soil zones through distributed linear springs with the Mohr-Coulomb failure criterion. The pullout of the threaded bars and development of the shear stress along the perimeter of the bars can be modeled by adjusting the grout adhesion and friction angle of the cable elements to soil-steel adhesion and friction angle. The matrix soil and RAP element are modeled with two sets of constitutive models, Mohr-Coulomb (MC) and Plastic Hardening (PH) models. The MC model, which is a widely used and accepted model in engineering practice, requires only five parameters. The FLAC3D MC model has a linear hooks’ law elasticity with nonhardening MohrCoulomb failure criterion and nonassociative flow rule for plastic strain. Bulk and shear elastic moduli, cohesion intercept and friction angle for yield surface, and dilatation angle for the plastic potential function are used to parameterize the FLAC3D MC model. The PH model in FLAC3D, which is implemented based on Schanz et al. (1999) formulation (Cheng and Detournay 2016; Schanz et al. 1999), has stress-dependent elastic moduli based on the Duncan and Chain model (Duncan et al. 1980). A Mohr-Coulomb yield criterion with a cap constitutes the PH model yield surface. The isotropic hardening of the yield surface is controlled by the evolution of a hardening parameter, which depends on both plastic shear and volumetric strain. The Mohr-Coulomb yield surface of the PH model has a nonassociative flow rule similar to the MC model in FLAC3D, and the cap has an 28 associative flow rule. A recently implemented Small-Strain Stiffness option in the PH model enables the shear modulus degradation due to an increase in shear strain. The PH model requires nine essential parameters and seventeen advanced optional parameters for full parameterization and activation of all of the model features. It can be parameterized with laboratory tests (Cheng and Lucarelli 2016) or in situ tests (Lucarelli and Cheng 2016). In this research, twenty parameters for the PH model were estimated from CPT sounding results and suggested material properties in the literature. The PH model parameterization is discussed in detail in the Constitutive Modeling section of this chapter. The modeling of the construction sequence of the RAPs is of paramount importance because it increases the lateral and vertical confining stresses in the matrix soil and RAP elements, which dominates the load-displacement behavior of the RAP system. RAPs construction is modeled with four different techniques. The first technique is a simple initialization of the horizontal stress within the matrix soil, and RAP element to match the measured values in the field after the RAPs construction. This technique does not consider the vertical stress increase due to RAP construction in modeling. The second technique simulates the compaction of RAP lifts with the application of static compaction stress on top of each lift (LaMeres 2005). The first RAP lift is placed at the bottom of the borehole, and increasing compressive stress is applied on top of the RAP lift to compact it statically and to build up the confining stress. After completion of the first lift compaction, the compressive stress is removed, and the second lift is placed on top of it. The lift placement and compaction are repeated to build the entire RAP. The third technique is based on the cavity expansion method that was proposed by Pham and White (Pham and White 2007) for numerical modeling of RAPs based on their FSTs measurement and observation (White 29 et al. 2007). In this technique, the shaft and the bottom of the excavated borehole are pushed outward and downward simultaneously to expand the cavity and increase the radial stress and vertical stress in the matrix soil. After reaching a predetermined radial and vertical strain at the shaft and bottom of the borehole, the cavity expansion is stopped, all the displacement restraints on the shaft and bottom of the borehole are removed, the RAP material is placed in the expanded borehole, and the model is solved to reach equilibrium. The fourth technique expands the cavity by applying variable velocity to the borehole shaft and bottom. The magnitude of the lateral displacement is variable through the depth of the borehole and throughout the simulation steps. The cavity expansion is done in steps, and a linearization scheme is used to estimate the new lateral displacement magnitude and pattern for each step. By assuming linear behavior for each step, each lift was expanded separately, and the increase in the effective lateral stress at the target locations was acquired in FLAC3D. The Jacobian matrix was assembled, which relates the borehole displacement to lateral earth pressure at each lift. By solving a system of linear equations, the new displacement pattern and magnitude were calculated to achieve the desired increase in effective lateral stress for the next step. This procedure is similar to calculating the stiffness matrix for a linear elastic structure by fixing all Degree of Freedoms (DOFs) and releasing and displacing them one at a time to get the increase in reaction in other DOFs. However, soil behavior is highly nonlinear, and the stress and strain history dependency of its behavior affects the efficiency of this method. After the cavity expansion phase is complete, RAP material is placed in the borehole, and the effective lateral and vertical stresses of the RAP zones in FLAC3D are initialized to match the induced stresses at the borehole boundary. This would reduce the amount of elastic rebound and loss of the 30 confining pressure. Finally, the model is solved to reach equilibrium. The uplift test was performed with a force control scheme by applying an upward incrementally increasing load at the top of the threaded rod and solving the model to reach equilibrium for each load increment. Figures 4.1 to 4.5 show the flowcharts for the numerical modeling of the RAPs uplift test in FLAC3D. 4.2 Model Geometry The RAPs are symmetric. However, the uplift plates consist of a circular plate and four threaded bars (see Figure 3.4). Therefore, a 2D axisymmetric model is not capable of capturing the 3D behavior of the uplift plate during the uplift test. The uplift plate has four planes of symmetry (see Figure 4.6 a). Modeling one-eighth of the RAP, matrix soil, and uplift plate with one-half of the threaded bar circular cross-section is sufficient to model the 3D behavior of the RAP during the uplift tests. However, one-quarter of the uplift RAPs systems were modeled in FLAC3D to simplify the modeling process. The FLAC3D model domain is a quarter of a cylinder with a 10 ft (3.05 ft) radius. The height of the cylinder depends on the length of the RAP. It is the height of the RAP plus four times the diameter of the RAP. As a result, the domain height is 15, 18, 21, 24, and 27 ft (4.57, 5.49, 6.40, 7.32, and 8.23 m) for the RAPs with 3, 6, 9, 12, and 15 ft (0.91, 1.83, 2.74, 3.66, and 4.57 m) height, respectively. The four times the diameter of the RAP of matrix soil beneath the RAP bottom was modeled to minimize the boundary effect during the RAP construction simulation, which induces considerable effective compressive vertical stress beneath the RAP. On the other hand, during the uplift simulation, the tensile force in the threaded bar is transferred to the uplift plate, which pushes the bottom of the RAP upward. The matrix 31 soil beneath the RAP has minimal effect on the load-displacement behavior of the uplift simulation. The height of the matrix soil beneath the RAP bottom is reduced to two times the RAP diameter (reduced from 12 ft (3.66 m) to 6 ft (1.83 m)) before the uplift simulation to decrease the computation time (see Figure 4.6 b). The RAP and matrix soil continua are discretized with an unstructured mesh (see Figure 4.7 a). The RAP and matrix soil zones close to it are discretized with a finer mesh, and the size of the soil zones increases farther away from the RAP. The unstructured mesh divides the RAP into seven zones in x and y direction and ten equidistant zones on the curved section. The length of the first six RAP zones counting from the center of the RAP, in x and y directions adjacent to the boundary, is 2.48 in. (63.0 mm) and the length of the last RAP zone adjacent to the model boundary is 3.12 in (79.4 mm). The matrix soil is divided into eighteen zones in the x and y direction with increasing grid spacing at a ratio of 1.05 and twenty equidistant zones on the curve section. The length of matrix soil in the x and y direction adjacent to boundary varies from 3.62 in (92 mm) to 8.32 in (211 mm). The mesh size of 3 in (76.2 mm) is constant through the height of the RAP in the zdirection. It increases at a ratio of 1.04 from the bottom of the RAP to the bottom of the model boundary (see Figure 4.7 b). The zones at the bottom of the model are 9.12 in. (232 mm) long in the z-direction. The unstructured mesh in the x-y plane is extruded through the z-direction to create the soil zones (see Figure 4.8). Generally, the FLAC3D unstructured mesh is not symmetric. In this research, the FLAC3D mesh parameters were tuned to generate an unstructured mesh that is very close to symmetric mesh. The advantage of the unstructured mesh is the smooth transition from finely meshed zones to coarsely meshed zones and maintaining an aspect ratio close to one for the mesh. 32 In other words, the unstructured mesh generates zones that are square or close to a square in the x-y plane. The accuracy of the numerical model increases as the zones’ aspect ratio tends to one, and the aspect ratio larger than ten should be avoided as a general rule of thumb (Itasca Consulting Group 2019). The maximum aspect ratio for zones within RAP is 1.2, and the maximum aspect ratio of the matrix soil zones at the bottom boundary and curved boundary are 3.7 and 2.8, respectively. Model size and grid statistic are presented in Table 4.1. The uplift plate is modeled with 2D liner structural elements (see Figure 4.9 a). The triangular shape liner elements have 6 DOFs (three translational and three rotational components) to model the full shell-membrane behavior of the uplift plate. The liner elements can interact with soil zones through their interface in both normal and tangent directions. They are defined as an embedded liner in the FLAC3D model to activate the interface on both sides of the uplift plate. At the bottom of the RAP, the RAP and matrix soil zones were separated, and another set of grid points at the same location was defined. The “zone separate” command in FLAC3D separates the zones at the specified location and duplicate the grid points (see Figure 4.9 b). The uplift plate is attached to the bottom of the RAP. The threaded rod is modeled using one-dimensional cable structural elements (see Figure 4.10). The threaded rod extends from the uplift plate to 0.1 ft (0.03 m) above the surface of the RAP. The threaded rod is discretized into 6, 12, 18, 24, and 30 elements for 3, 6, 9, 12, and 15 ft (0.91, 1.83, 2.74, 3.66, and 4.57 m) long RAPs, respectively. The length of each cable element is about 6 in. (152 mm). The threaded rod is connected to the uplift plate though a rigid structural link between the bottom node of the first cable element 33 and the adjacent structural node of the uplift plate. The rigid structural link constraint the displacement of the two nodes in the z-direction. The boundary condition for the model is imposed on both the grid points and structural nodes on the boundary of the model. The soil grid points on the x-side, y-side, and the curve-side of the model are restraint by setting the normal velocity to zero. Zero normal velocity would simulate a roller boundary for these grid points. The grid points at the bottom of the model are restraint in all directions by applying a zero velocity in x, y, and z directions. The grid points at the top of the model are free. The uplift plate structural nodes on the x-side of the boundary are restrained on their translational DOF in the ydirection and rotational DOFs in the x and z directions. Likewise, the translational DOF in the x-direction and rotational DOFs in y and z directions were set to zero for structural nodes on the y-side of the model boundary. The boundary conditions that are mentioned here are presented in the global coordinate system for convenience. However, the structural elements use local coordinate to form mass and stiffness matrices in FLAC3D, and the boundary conditions are specified in the local coordinate system. 4.3 Constitutive Model 4.3.1 Structural Elements The liner elements in FLAC3D have a triangular shape with three nodes. The number of DOFs at the nodes of a liner element could vary depending on the type of the FEM formulation that is chosen for the liner element. For example, a liner element with two translational DOFs (FLAC3D element name: CST Plane-Stress Element) at each node can only capture the membrane type behavior. The selected element type in this research 34 (FLAC3D element name: DKT-CST Hybrid Shell Element) has 6 DOFs at each node, which enables the liner element to capture full shell-membrane behavior (see Figure 4.11). The liner structural elements have a linear elastic constitutive model. In addition to shellmembrane structural behavior, the liner element can interact with the FLAC3D grid on both sides in the normal and tangent direction. It can carry both compressive and tensile stress int the normal direction. The liner element can break free from the grid and come into contact with it again. The interface behavior in the normal direction is formulated in FLAC3D by a linear spring at each node with finite tensile strength. In the tangent direction, the liner element interacts with the FLAC3D grid through a frictional interaction. A spring-slider at each node of the liner node is used to formulate the tangent interface in FLAC3D. Figure 4.11 shows the normal and tangent behavior of the FLAC3D interface nodal spring. The modeled uplift plate has the radius of the 18 in. (457 mm) and thickness of 0.5 in (12.5 mm). The modeled uplift plate radius is 1 in. (25.4 mm) larger than the uplift plate in FSTs. The slightly larger uplift plate eliminates the need to discretize the model to 1 in. zones near the edge of the uplift plate. Also, the slightly larger uplift plate is the same diameter as the RAPs in the FSTs. The selected isotropic material properties for the A36 steel plate are Young’s modulus (𝐸𝑠𝑡 ) equal to 29.0 × 106 psi (200 GPa) and Poisson’s ratio equal to 0.26. The steel unit weight is set to 486.9 pcf (7800 Kg/m3) in the FLAC3D model. FLAC3D manual suggests that stiffness of the interface springs in normal and shear direction be set to ten times the equivalent stiffness of the neighboring zone. The following formula is used in FLAC3D model to calculate and assign the normal and shear coupling spring stiffness (𝑘𝑛 and 𝑘𝑠 in Figures 4.11 c and 4.11 c, respectively): 35 3 (𝐾 + 4 𝐺) 𝑘𝑛 = 𝑘𝑠 = 10 × [ ] ∆𝑧𝑚𝑖𝑛 (4.1) where: 𝑘𝑛 = Normal and shear coupling spring stiffness 𝑘𝑠 = Shear and shear coupling spring stiffness 𝐾 = Bulk modulus 𝐺 = Shear modulus ∆𝑧𝑚𝑖𝑛 = The smallest dimension of an adjoining zone in the normal direction Bowles (1996) suggested 22° for the interface friction angle between “clean gravel, gravel-sand mixture, well-graded rock fil with spalls” and steel sheet piles (Bowles 1996). Potyondy (1961) proposed the ratio of the interface friction angle between soil and steel to the soil internal friction value of 0.76 for dry dense sand and rough steel, and 0.80 for saturated dense sand and rough steel (Potyondy 1961). The threads on the threaded bar may cause the development of passive resistance during the pullout, which should be considered in the selection of the interaction friction angle. The steel strap in Mechanically Stabilized Earth (MSE) wall may have the friction coefficient varying from tan∅ to 1.5 for deeper depth to shallow depth where the ∅ is the friction angle of the soil. The suggested interface friction angle for rod reinforcement for soil nailing is equal to the soil fiction angle (Mitchell and Villet 1987). The uplift plate is embedded in the RAP material. The fiction angle between the 1 uplift plate and the RAP is selected as 𝑡𝑎𝑛−1 (2 𝑡𝑎𝑛∅𝑅𝐴𝑃 ) based on the suggested values for interface friction angle and engineering judgment. It is assumed no adhesion exists between the uplift plate and the granular RAP material. Therefore, the values of the normal 36 coupling spring tensile strength (𝑓𝑡 in Figure 4.11 d), and shear coupling spring cohesion (𝑐𝑙𝑖𝑛𝑒𝑟 in Figure 4.11 d) are set to 0.0 psf. The cable element in FLAC3D is a one-dimensional element with two nodes (see Figure 4.12 a). The cable elements have elastic, perfectly plastic behavior. It can yield in tension and compression (see Figure 4.12 b). The cable element FEM formulation of stiffness and mass matrices are based on truss element FEM formulation. Shear stress can develop along the side of the cable element through the grout that surrounds the perimeter of the cable element. The grout is numerically formulated with spring-slider at nodal point s of the cable element. The cable elastic material properties and unit weight are similar to the uplift plate. Young’s modulus (𝐸𝑠𝑡 ) of 29.0 × 106 psi (200 GPa), Poisson’s ratio of 0.26, and steel unit weight of 486.9 pcf (7800 Kg/m3) were assigned to the FLAC3D cable elements. No. 8 grade 75 threaded bar was modeled for all the RAP models. The modeled threaded bar has 1 in. (25.4 mm) diameter and minimum yield tensile strength of 75 ksi (517 MPa). The calculated cross-sectional area of 0.7854 in2 (507 mm2) and compressive/tensile strength of 58.9 kips (262 kN) were assigned to the cable elements. The grout annulus behavior is defined by grout shear stiffness per unit length (𝑘𝑔 ), grout cohesive strength per unit length (𝑐𝑔 ), grout friction angle (∅𝑔 ), and grout exposed perimeter (𝑝𝑔 ) (see Figure 4.12). FLAC3D manual suggested the following formula to calculate the grout shear stiffness: 𝑘𝑔 ≈ where: 2𝜋𝐺𝑔 2𝑡 10 × 𝑙𝑛 (1 + 𝐷𝑎𝑛𝑛𝑢𝑙𝑢𝑠 ) 𝑐𝑎𝑏𝑙𝑒 (4.2) 37 𝑘𝑔 = Grout shear stiffness per unit length 𝑡𝑎𝑛𝑛𝑢𝑙𝑢𝑠 = Grout annulus thickness 𝐷𝑐𝑎𝑏𝑙𝑒 = Reinforcement diameter 𝐺𝑔 = Grout shear modulus The FSTs’ threaded rods were surrounded with RAP material without any grouting. A virtual grout with annulus thickness equal to one-tenth of the threaded rod diameter is assumed to enable the FLAC3D model to develop frictional shear stress on the perimeter of the threaded rod during the pullout test. The shear modulus of the grout is equal to the RAP material shear modulus. Therefore, Equation 4.2 is simplified to: 𝑘𝑔 ≈ 2𝜋𝐺𝑅𝐴𝑃 = 3.446𝐺𝑅𝐴𝑃 2 × 0.1 × 𝐷𝑐𝑎𝑏𝑙𝑒 10 × 𝑙𝑛 (1 + ) 𝐷𝑐𝑎𝑏𝑙𝑒 (4.3) where: 𝐺𝑅𝐴𝑃 = RAP material shear modulus The grout exposed perimeter is: 𝑝𝑔 = 𝜋 (𝐷𝑐𝑎𝑏𝑙𝑒 + 2𝑡𝑎𝑛𝑛𝑢𝑙𝑢𝑠 ) = 1.2𝜋𝐷𝑐𝑎𝑏𝑙𝑒 (4.4) The grout cohesive strength per unit length is set to 0.0 because the granular RAP material is assumed to be cohesionless. The threaded bar is embedded in RAP material, which is well-graded gravelly sand above the water table and 2 to 4 in. rock below the water table. The well-graded granular material has more points of contact with the threaded rod. However, the larger size RAP material below the water table most likely has less point of contact with the threaded bar. The tangent of the interface friction coefficient (𝑡𝑎𝑛∅𝑔 ) is selected equal to two-third of 2 the tangent of the RAP material friction angle (3 𝑡𝑎𝑛∅𝑅𝐴𝑃 ) above the water table and half 38 1 the tangent of the RAP material friction angle (2 𝑡𝑎𝑛∅𝑅𝐴𝑃 ) below the water table. The soilthreaded bar interaction friction angle is as follows: 2 𝑡𝑎𝑛−1 ( 𝑡𝑎𝑛∅𝑅𝐴𝑃 ) , 3 ∅𝑔 = { 1 𝑡𝑎𝑛−1 ( 𝑡𝑎𝑛∅𝑅𝐴𝑃 ) , 2 𝑧 < 4 𝑓𝑡 (4.5) 𝑥 ≥ 4 𝑓𝑡 where: ∅𝑔 = Grout friction angle ∅𝑅𝐴𝑃 = RAP material friction angle 4.3.2 Soil Zones 4.3.2.1 General Native Soil Properties Many of the soil parameters in this research are estimated based on the selected CPT sounding that was performed in the vicinity of the FSTs area. Seven out of a total of ten CPT soundings were performed within the native soils. By comparing CPTs tip resistance and sleeve friction, Lawton showed that soil strata are approximately horizontal at FSTs site location (Lawton 2000). CPT 3 is selected to characterize the native soil because it was performed close to the FSTs location, and it went down to the depth of 30 ft (9.14 m) below the ground surface. The ground elevation of the CPT 3 is about 1 ft (0.30 m) higher than the CPT 9 ground surface, which is the closest to the FSTs location but only was pushed down to the depth of 15 ft. Figure 3.1 shows the CPT 3 sounding with depth correction. The saturated unit weight of the soil is estimated with CPT correlation by Mayne and Ozer et al. (Mayne 2007; Ozer et al. 2013). The saturated unit weight (𝛾𝑠𝑎𝑡 ) for soil profile from the ground surface to the depth of 19 ft is calculated using the CPT correlation 39 suggested by Mayne as following: 𝛾𝑠𝑎𝑡 = 2.6 𝑙𝑜𝑔(𝑓𝑠 ) + 15 𝐺𝑠 − 26.5 (4.6) where: 𝛾𝑠𝑎𝑡 = Unit weight of saturated soil (kN/m3) 𝑓𝑠 = CPT sleeve friction (KPa) 𝐺𝑠 = Specific gravity Equation 4.6 is unit specific. The same correlation in the empirical unit is as follows: 𝛾𝑠𝑎𝑡 = 6.3659 × [2.6 𝑙𝑜𝑔(47.88 𝑓𝑠 ) + 15 𝐺𝑠 − 26.5] (4.7) where: 𝛾𝑠𝑎𝑡 = Unit weight of saturated soil (pcf) 𝑓𝑠 = CPT sleeve friction (ksf) 𝐺𝑠 = Specific gravity The soil above the groundwater table is considered to be saturated because of the capillarity effect. The total unit weight for the Lake Bonneville clay (soil below 19 ft) is calculated based on the correlation proposed by Ozer et al. as follows: 𝑞𝑡 0.148 (𝑓𝑅 )0.0144 𝛾 = 1.27𝛾𝑤 ( ) 𝑃𝑎 where: 𝛾 = Unit weight of soil (pcf) 𝛾𝑤 = Unit weight of water 𝑞𝑡 = CPT corrected, total tip resistance 𝑃𝑎 = Atmospheric pressure (4.8) 40 𝑓𝑅 = Friction ratio in % The groundwater table was measured at 4 ft below the ground surface. Therefore, the unit weight calculated from Equation 4.8 is the saturated unit weight. The corrected, total tip resistance (𝑞𝑡 ) is calculated by taking into account the unbalanced pore pressure force behind the tip of the cone penetrometer. The 𝑞𝑐 and 𝑞𝑡 values are assumed to be the same in this research because 𝑞𝑐 and 𝑞𝑡 values are usually close, and no conversion coefficient is available. The length of FLAC3D model zones in the z-direction varies from 3 in (76.2 mm) to 9.12 in. (232 mm). Each zone is the continuum representation of soil particles within its domain. Therefore, the average value for soil parameters within the zones’ domain should be used for the saturated unit weight (𝛾𝑠𝑎𝑡 ) and constitutive model parameters. The CPT data are collected every 2 in. (50.8 mm). Therefore, the calculated values for the saturated unit weight (𝛾𝑠𝑎𝑡 ) from Equations 4.7, and 4.8 are averaged every 6 in. (152 mm), which is about the average size of the FLAC3D zone’s length in the z-direction. Taking the moving average of the data smoothed out the saturated unit weight (𝛾𝑠𝑎𝑡 ) slightly. The saturated unit weight (𝛾𝑠𝑎𝑡 ) is imported in FLAC3D as a table. The FLAC3D finds the saturated unit weight (𝛾𝑠𝑎𝑡 ) of each zone by linearly interpolating the imported saturated unit weight (𝛾𝑠𝑎𝑡 ) table for the z-coordinate of the center of each zone, divides it by the gravity acceleration and assigns it to the zone density property. The saturated unit weight (𝛾𝑠𝑎𝑡 ) and calculated initial total vertical stress (𝜎𝑣0 ) are plotted in Figure 4.13. The over consolidation pressure (𝜎𝑝′ ) of the soil profile is estimated using the CPT correlation by Mayne for soil from the ground surface to 19 ft (5.79 m) depth and the CPT correlation by Ozer for Lake Bonneville clay at depth below 19 ft (5.79 m) (Mayne 2007; 41 Ozer 2005). Mayne suggested CPT correlation to estimate the 𝜎𝑝′ is: 𝜎𝑝′ = 0.33(𝑞𝑡 − 𝜎𝑣0 ) (4.9) where: 𝜎𝑝′ = Over consolidation pressure 𝜎𝑣0 = Initial total vertical stress Ozer’s CPT correlation for Lake Bonneville clay is: 𝜎𝑝′ = 0.632𝑃𝑎 ( 𝜎𝑣0 0.564 𝑞𝑡 − 𝜎𝑣0 0.342 ) ( ) 𝑃𝑎 𝑃𝑎 (4.10) The Over Consolidation Ratio (OCR) is: 𝜎𝑝′ 𝑂𝐶𝑅 = ′ 𝜎𝑣0 (4.11) where: 𝑂𝐶𝑅 = Over Consolidation Ratio ′ 𝜎𝑣0 = Initial effective vertical stress The moving average of the OCR is calculated and imported in FLAC3D as a table. ′ Figure 4.14 shows the initial effective vertical stress (𝜎𝑣0 ), and smoothed out OCR. The coefficient of lateral earth pressure at rest (𝐾0 ) is measured by the K0-Blade test (see K at far-field in Figure 3.3). However, the measured values are only available to the depth of less than 12 ft (3.66 m). The 𝐾0 values are needed to initialize the lateral stress in the FLAC3D model. Two CPT correlations were evaluated by comparing the estimated 𝐾0 from correlations and K0-Blade test measured values. The first CPT correlation is applicable for uncemented sands and clays of low to medium sensitivity (Mayne 2007): ′ 𝐾0 = (1 − 𝑠𝑖𝑛∅′ )𝑂𝐶𝑅 𝑠𝑖𝑛∅ ≤ 𝐾𝑝 (4.12) 42 𝐾𝑝 = 1 + 𝑠𝑖𝑛∅′ 1 − 𝑠𝑖𝑛∅′ (4.13) where: 𝐾0 = Coefficient of lateral earth pressure at rest 𝐾𝑝 = Rankine’s coefficient of passive earth pressure ∅′ = Effective friction angle Mayne and Campanella (2005) proposed the following equation to estimate the effective friction angle (∅′ ) for mixed soil type: ∅′ = 29.5° × 𝐵𝑞0.121 (0.256 − 0.336𝐵𝑞 + 𝑙𝑜𝑔𝑄𝑡 ) (4.14) where: 𝐵𝑞 𝑎𝑛𝑑 𝑄𝑡 = Normalized CPT parameters The normalized CPT parameters are calculated using the following formulas: 𝐵𝑞 = 𝑢2 − 𝑢0 𝑞𝑡 − 𝜎𝑣0 (4.15) 𝑄𝑡 = 𝑞𝑡 − 𝜎𝑣0 ′ 𝜎𝑣0 (4.16) where: 𝑢2 = Pore water pressure measured during CPT sounding 𝑢0 = In situ equilibrium pore water pressure Equation 4.14 is only applicable for 0.1 < 𝐵𝑞 < 1, and 20° < ∅′ < 45°. For 𝐵𝑞 < 0.1, Mayne suggested the use of the following correlation, which is proposed for clean sands: ∅′ = 17.6° + 11.0° × 𝑙𝑜𝑔𝑞𝑡1 (4.17) 43 𝑞𝑡1 = 𝑞𝑡 ⁄𝑃𝑎 (4.18) ′ ⁄ 𝑃𝑎 √𝜎𝑣0 where: 𝑞𝑡1 = Normalized CPT tip resistance Mayne (2007) proposed a second 𝐾0 CPT correlation for clean sand as follows: 0.31 𝑞𝑡 0.22 𝑃𝑎 𝐾0 = 0.192 ( ) ( ′ ) 𝑃𝑎 𝜎𝑣0 𝑂𝐶𝑅 0.27 ≤ 𝐾𝑝 (4.19) The 𝐾0 values are estimated using Equation 4.12 by selecting the relevant correlation for ∅′ from Equations 4.14 and 4.17 at each depth, and using Equation 4.12. These 𝐾0 values are plotted against the measured 𝐾0 values in Figure 4.15. The soil type at the top 12 ft (3.66 m) of the soil profile is mostly lean clay. The Equation 4.19, which is proposed for clean sand, has a better prediction of the coefficient of lateral earth pressure at rest (𝐾0 ) compared to Equation 4.12, which is applicable for clay with low to medium sensitivity. Each CPT correlation introduces limitations, uncertainty, and error in the predicted value. And, Equation 4.12 uses two correlation, one to predict the ∅′ (Equations 4.14 and 4.17) and another one to estimate the 𝐾0 (Equation 4.12). This could potentially increase the error in predicted 𝐾0 values. Equation 4.18 is used to compliment the K0-Blade test results and estimate the 𝐾0 where measured values are not available. Figure 4.16 shows the selected 𝐾0 for initializing the horizontal stress in FLAC3D. The 𝐾0 values from CPT correlation are smoothed out (moving average) in this plot. 4.3.2.2 General RAP Material Properties Lawton (2000) assumed total unit weight (𝛾𝑅𝐴𝑃 ) of 135 pcf (21.2 kN/m3) above the ′ groundwater table and an effective unit weight (𝛾𝑅𝐴𝑃 ) of 75 pcf (11.8 kN/m3) below the 44 groundwater table for the RAP material. Assuming 62.4 pcf (9.81 kN/m3) for unit weight of water, the saturated unit weight of the RAP material (𝛾𝑠𝑎𝑡,𝑅𝐴𝑃 ) would be 137.2 pcf (21.6 kN/m3). The value of 135 pcf (21.2 kN/m3) is assumed for both 𝛾𝑅𝐴𝑃 above the groundwater table and 𝛾𝑠𝑎𝑡,𝑅𝐴𝑃 below the groundwater table in the FLAC3D model. 4.3.2.3 Mohr-Coulomb Model The well-known and widely used Mohr-Coulomb model was chosen as the first constitutive model for matrix soil and RAP material. The FLAC3D’s MC constitutive model is parametrized by elastic bulk modulus (𝐾), elastic shear modulus (𝐺), cohesion (𝑐), effective friction angle (∅′ ), and dilatation angle (𝜓). The bulk modulus and shear modulus are calculated from Young’s modulus and Poisson’s ratio. Bowles suggested the following CPT correlation for normally consolidated clayey sand’s Young’s modulus (Bowles 1996): 𝐸𝑠,𝑁𝐶 = (3 𝑡𝑜 6)𝑞𝑐 (4.20) where: 𝐸𝑠,𝑁𝐶 = Young’s modulus for normally consolidated soil 𝑞𝑐 = CPT tip resistance He suggested the following correlation to calculate Young’s modulus of normally consolidated silts, sandy silt, or clayey silt: 𝐸𝑠,𝑁𝐶 = (1 𝑡𝑜 2)𝑞𝑐 (4.21) LaMeres (2005) used the range suggested by Equations 4.9 and 4.10 to estimate the 𝐸𝑠,𝑁𝐶 for his 2D axisymmetric FLAC model of the RAP FSTs that was performed in the winter of 2002 under the northbound bridge of the I-15 at South Temple Street in Salt Lake 45 City, Utah (LaMeres 2005). The CPT correlation that LaMeres used is: 𝐸𝑠,𝑁𝐶 = (1 𝑡𝑜 6)𝑞𝑐 (4.22) The same CPT correlation is adopted in this research because of the two FSTs’ sites are very close, and LaMeres (2005) numerical model results had good agreement with FSTs load-displacement measurement. Young’s modulus should be modified with the following formula to account for the Overconsolidation of the soil (Bowles 1996): 𝐸𝑠,𝑂𝐶 = 𝐸𝑠,𝑁𝐶 √𝑂𝐶𝑅 (4.23) where: 𝐸𝑠,𝑂𝐶 = Young’s modulus for overconsolidated soil 𝑂𝐶𝑅 = Over consolidation ratio The recommended range of values for Poisson’s ratio (𝜈) is shown in Table 4.2. The RAP construction induces a considerable amount of confining pressure within the RAP and matrix soil (Handy and White 2005a; Lawton 2000; Lawton et al. 1994). The high confining pressure and high permeability of the RAP material facilitate and expedite the excess pore water pressure dissipation. In addition to the RAP material drainage, the formation of radial tension cracking and hydraulic fracturing at the matrix soil elastic zone1 opens up a new drainage path during the RAP construction process (Handy and White 2005b). Therefore, high values of 0.45 to 0.5 for Poisson’s ratio (𝜈) are avoided for the FLAC3D matrix soil elastic properties because the RAP construction and uplift test (due 1 The stress state in the matrix soil in contact with RAP and in the vicinity of the RAP can reach a plastic state (plastic zone). The stress state at the shallow depth where the overburden pressure is low may reach a passive stress state (passive zone). And the matrix soil farther away from the RAP remains elastic (elastic zone) (Handy and White 2005a). 46 to a slow rate of loading) are in drained condition for modeling purposes. The drained condition assumption validity may diminish for matrix soil farther way from the RAP. However, the preliminary modeling results and FSTs results suggest that the matrix soil immediately surrounding the RAP has a prominent role in the uplift behavior of the RAP system (Hsu 2000; LaMeres 2005; Lawton 2000). The Poisson’s ratio values are selected based on the recommended values for each soil type and engineering judgment. Figure 4.16 shows the soil profile’s Poisson’s ratio and Young’s modulus modified for over consolidation. The bulk and shear modulus are calculated using the following elasticity formulas: 𝐾= 𝐸𝑠,𝑂𝐶 3(1 − 2𝜈) (4.24) 𝐺= 𝐸𝑠,𝑂𝐶 2(1 + 𝜈) (4.25) The moving average of the bulk and shear modulus are calculated and imported to FLAC3D (see Figure 4.17). The effective friction angle (∅′ ) and cohesion (𝑐 ′ ) of the soil were measured in the filed with a series of BSTs (see Table 3.1 and Figure 3.2). The soil strength parameters are available for normally consolidated (NC) and overconsolidated (OC) conditions. The RAP construction process disturbs the soil immediately surrounding the RAP element and remolds the cohesive matrix soil in the plastic and passive zone (in granular matrix soils, the RAP construction liquefies the soil close to RAP element). The remolding generally would result in complete loss of cohesion intercept. The cohesive soil will regain its strength in time due to the thixotropy process, which is the attraction of ions in the diffuse double layer and edge to face bonding of clay particles (Handy and White 2005a; Hsu 47 2000; Lawton 2000). However, the uplift FSTs were performed shortly after the RAPs construction, and it is assumed that the soil strength gain was negligible in a short time due to thixotropy. Lawton (2000) proposed using the NC soil strength parameters with 𝑐 ′ = 0. The matrix soil farther away from the RAP remains elastic (elastic zone), and no remolding would happen. Setting the cohesion intercept to zero in the FLAC3D model would lead to underestimation of the ultimate capacity. The cohesion intercept can preserve the horizontal confining pressure that is induced during the RAP construction process. The RAP construction process in this research is either using a method to expand the cavity and place RAP inside the borehole after the expansion or using static compaction by applying pressure on top of the RAP material. After reaching the desired level for horizontal confining pressure, the constraints on the borehole are relaxed, and the model is solved to reach equilibrium. Without a cohesion intercept, soil at shallow depth (with low overburden pressure) would lose a considerable amount of the build-up horizontal confining pressure. There are few mechanisms that can preserve the built-up horizontal stress after RAP construction in the field. For example, the vibratory dynamic compaction of the RAP granular material rearranges and interlocks the soil particles, which creates large interparticle stresses. The compacted soil is very stiff and resists the rebound and confining stress relaxation. The numerical modeling of this process requires mesoscale modeling of the soil particle which is beyond FLAC3D’s limitation. Another process that can preserve the confining pressure is the intrusion of the soil particle in the radial cracks and wedge formation in the elastic zone (Handy and White 2005b). The numerical modeling of this process requires fracture modeling which is outside the scope of this thesis 48 and FLAC3D capability. The measured matrix soil cohesion intercept is considered in RAP construction modeling efforts by other researchers as well (LaMeres 2005; Pham and White 2007). Therefore, to preserve the confining pressure at shallow depth, the matrix soil cohesion is linearly decreased from 150 psf (7.2 kPa) to 50 psf (2.4 kPa) from the ground surface to the depth of 6 ft (1.83 m) in the FLAC3D model. The 150 psf (7.2 kPa) is the average measured cohesion with BST. The cohesion intercept is set to 50 psf (2.4 kPa) for the matrix soil below 6ft (1.83 m). The NC friction angle from BSTs and the described cohesion values are used to parameterize the FLAC3D MC constitutive model. The BST data are only available to the depth of 15.5 ft (4.72 m), which is deeper than the longest modeled RAP with 15 ft (4.57 m) length. In other words, the high-quality BST strength parameters are available for matrix soil surrounding the RAP elements. The effective friction angle of soil below 15.5 ft (4.72 m) was estimated by CPT correlations (Equations 4.14 and 4.17). The measured ∅′ from BSTs and estimated ∅′ from CPT correlation are plotted in Figure 4.18 a. The moving average of the friction angle below 15.5 ft is calculated and complemented the BST friction angle data. These data were imported in FLAC3D to assign the effective friction angle for the MC model (see Figure 4.18 b). The dilatation angle (𝜓) is assumed to be zero for cohesive material. Lawton (2000) estimated the effective friction angle of ∅′𝑅𝐴𝑃 = 50° for the RAP element from the ground surface to the depth of 2.5 ft (0.76 m) , and ∅′𝑅𝐴𝑃 = 47° for the RAP element from a depth of 2.5 ft (0.76 m) to 9 ft (2.74 m) base on a CPT sounding ′ within the RAP material. The cohesion intercept for the RAP granular material 𝑐𝑅𝐴𝑃 is considered zero. The effective friction angle for the RAP element is set to 47°, and the 49 cohesion intercept is to 1 psf (0.05 kPa) in the FLAC3D model. The bulk and shear modulus are calculated based on Young’s modulus of 4000 psf (192 MPa), which was used by LaMeres (2005) for his FLAC modeling of RAP systems, and Poisson’s ratio of 0.4 for dense granular material. Using Equations 4.24 and 4.25, the bulk modulus (𝐾𝑅𝐴𝑃 ) of 6,670 ksf (319 MPa), and shear modulus (𝐺𝑅𝐴𝑃 ) of 1,430 ksf (68.5 MPa) were selected and assigned to RAR material in FLAC3D. White et al. (2007a) performed consolidated-drained (CD) triaxial tests on the reconstituted samples of their RAP material. They reported the effective friction angle (∅′𝑅𝐴𝑃 ) of 47° and dilatation angle (𝜓𝑅𝐴𝑃 ) ranged from 9° of 16° for confining stress of 731 psf (35 kPa) to 3,590 psf (172 kPa). Their measured value for the dilatation angle has a good agreement with the recommended dilatation angle estimated by subtracting 30° from the effective friction angle (𝜓 ≈ ∅′ − 30°). White et al. (2007b) used the average value of 12° for the numerical modeling of RAP system. The same value of 𝜓𝑅𝐴𝑃 = 12° is adopted in this research for FLAC3D modeling due to unavailability of triaxial test results for the modeled RAP material. 4.3.2.4 Plastic-Hardening Model 4.3.2.4.1 Theoretical background and FLAC3D formulation This section discusses the FLAC3D PH model formulation (Cheng and Detournay 2016; Itasca Consulting Group 2019). The PH constitutive model has a hyperbolic elastic model (Duncan-Chang model) with a Mohr-Coulomb yield surface enclosed by a cap. The PH model yield surface is not fixed in the principal effective stress and can expand according to its hardening law after either a shear failure (Mohr-Coulomb yield surface) or 50 a volumetric failure (cap yield surface). The PH model was formulated based on effective stress. The stresses are taken positive in tension and negative in compression to be consistent with the FLAC3D sign convention. The principal stresses and effective principal stress are 𝜎1 ≤ 𝜎2 ≤ 𝜎3 and 𝜎1′ ≤ 𝜎2′ ≤ 𝜎3′ , respectively. In the case of a drained triaxial test, the hyperbolic elastic model defines the relationship between the axial strain and deviatoric stress by: 𝜀1 = 𝑞𝑎 𝑞 𝑓𝑜𝑟 𝑞 < 𝑞𝑓 𝐸𝑖 (𝑞𝑎 − 𝑞) (4.26) 𝑞 = 𝜎3 − 𝜎1 = 𝜎3′ − 𝜎1′ 𝑞𝑓 𝑅𝑓 (4.28) 2𝑠𝑖𝑛𝜙(𝑐 𝑐𝑜𝑡𝜙 − 𝜎3′ ) 1 − 𝑠𝑖𝑛𝜙 (4.29) 𝑞𝑎 = 𝑞𝑓 = (4.27) where: 𝜀1 = Axial strain 𝑞 = Deviatoric stress 𝑞𝑎 = Deviatoric stress at failure, (𝜎3 − 𝜎1 )𝑓 𝑞𝑓 = Deviatoric stress from Mohr-Coulomb failure criterion 𝐸𝑖 = Initial tangent modulus 𝑅𝑓 = Failure ratio The 𝑅𝑓 has a value of less than one, and it varies from 0.5 to 0.9 for most soil (Duncan et al. 1980). The recommended value of 𝑅𝑓 for most soil is 0.9 (Itasca Consulting Group 2019; Schanz et al. 1999). The PH model uses the secant modulus at half the 1 deviatoric stress from the Mohr-Coulomb failure criterion (2 𝑞𝑓 ) instead of the initial 51 tangent modulus. The tangent modulus is easier to determine experimentally, and it is related to the initial tangent modulus by: 𝐸𝑖 = 2𝐸50 2 − 𝑅𝑓 (4.30) where: 𝐸50 = Secant modulus at half the deviatoric stress from the Mohr-Coulomb failure criterion The 𝐸50 follows a power law: 𝑚 𝐸50 = 𝑟𝑒𝑓 𝐸50 𝑐 𝑐𝑜𝑡𝜙 − 𝜎3′ ( ) 𝑐 𝑐𝑜𝑡𝜙 + 𝑝𝑟𝑒𝑓 (4.31) where: 𝐸50 = Secant modulus at half the deviatoric stress from the Mohr-Coulomb failure criterion 𝑟𝑒𝑓 𝐸50 = Reference secant modulus at half the deviatoric stress from the Mohr- Coulomb failure criterion corresponding to the reference stress (𝑝𝑟𝑒𝑓 ) 𝑚 = Exponent that determine the amount of stress dependency of 𝐸50 𝑝𝑟𝑒𝑓 = Reference stress The exponent “m” is close to 1 for clay, and it is usually within the range of 0.5 and 1 for sands (Itasca Consulting Group 2019). The reference stress (𝑝𝑟𝑒𝑓 ), the minor stress 50 that corresponds to 𝐸𝑟𝑒𝑓 , is a user input parameter and has a positive value for compressive stress (unlike the stress sign convention). The unloading follows a steeper stress-strain curve. The unloading secant modulus is calculated with a stress-dependent relationship similar to Equation 4.31: 52 𝑚 𝐸𝑢𝑟 = 𝑟𝑒𝑓 𝐸𝑢𝑟 𝑐 𝑐𝑜𝑡𝜙 − 𝜎3′ ( ) 𝑐 𝑐𝑜𝑡𝜙 + 𝑝𝑟𝑒𝑓 (4.32) where: 𝐸𝑢𝑟 = Unloading Young’s modulus 𝑟𝑒𝑓 𝑟𝑒𝑓 𝐸𝑢𝑟 = Unloading Young’s modulus at the reference stress (𝑝 ) Figure 4.19 shows the plots of the 𝐸50 and 𝐸𝑢𝑟 normalized by their reference values versus the normalized minor stress for different values of “m”. The amount of stressdependency is lower for soils with higher relative strength parameters. The PH model yield surface is defined by: 𝑞𝑎 𝐸𝑖 𝐸𝑖 ) 𝑞 − 𝛾𝑝 = 0 𝑓𝑀𝐶 = ( − 𝑞𝑎 − 𝑞 𝐸𝑢𝑟 2 (4.33) where: 𝑓𝑀𝐶 = PH model shear yield surface 𝛾 𝑝 = Internal shear hardening parameter The 𝛾 𝑝 is the internal hardening parameter that controls the evolution of the yield surface, and it is defined in its incremental form by: 𝑝 𝑝 𝑝 ∆𝛾 𝑝 = −(∆𝜀1 − ∆𝜀2 − ∆𝜀3 ) (4.34) where: Δ𝛾 𝑝 = Internal shear hardening parameter increment The ∆𝜀1𝑝 , ∆𝜀1𝑝 , and ∆𝜀1𝑝 are the increment in principal plastic strain. The shear yield surface defined in Equation 4.33 can increase in size by evolving isotopically, but it is ultimately confined by a Mohr-Coulomb failure criterion. The PH model has a nonassociated flow rule consistent with the Mohr-Coulomb yield criterion, and it is parameterized with a dilation angle (𝜓). The shear yield surface is confined by a cap. And, 53 the cap yield function is defined by: 𝑓𝐶𝐴𝑃 𝑞̃ 2 = 2 − (𝑝′ )2 − 𝑝𝑐2 = 0 𝛼 𝑞̃ = −[𝜎1′ + (1 − 𝛿)𝜎2′ − 𝛿𝜎3′ ] 𝜁= 1 + 𝑠𝑖𝑛𝜙 ′ 1 − 𝑠𝑖𝑛𝜙 ′ 1 𝑝′ = − (𝜎1′ + 𝜎2′ + 𝜎3′ ) 3 (4.35) (4.36) (4.37) (4.38) where: 𝑞̃ = A special shear stress measure 𝛼 = A constant derived internally in FLAC3D based on a virtual oedometer test 𝑝′ = Effective mean pressure 𝑝𝑐 = Volumetric hardening parameter 𝜁 = An internal parameter The initial value for the volumetric hardening parameters which denotes the preconsolidation pressure is calculated in FLAC3D by: 𝑝𝑐,𝑖𝑛𝑖 = 𝑂𝐶𝑅√ 2 𝑞̃𝑖𝑛𝑖 ′ )2 + (𝑝𝑖𝑛𝑖 𝛼2 (4.39) where: 𝑝𝑐,𝑖𝑛𝑖 = The initial value of the volumetric hardening parameter 𝑞̃𝑖𝑛𝑖 = The initial value of the special shear stress measure ′ 𝑝𝑖𝑛𝑖 = Initial effective mean pressure The volumetric hardening parameters evolve incrementally with the following relationship: 54 𝑐 𝑐𝑜𝑡𝜙 + 𝑝𝑐 𝑚 𝑣 ) ∆𝛾 Δ𝑝𝑐 = 𝐻𝑐 ( 𝑐 𝑐𝑜𝑡𝜙 + 𝑝𝑟𝑒𝑓 (4.40) ∆𝛾 𝑣 = −(∆𝜀1𝑝 + ∆𝜀2𝑝 + ∆𝜀3𝑝 ) (4.41) where: 𝐻𝑐 = Hardening modulus for cap hardening ∆𝛾 𝑣 = Internal volumetric hardening parameter increment Although the FLAC3D calculates the parameters 𝛼 and 𝐻𝑐 internally, it is the author’s experience that the internally calculated 𝛼 and 𝐻𝑐 may result in a numerically unstable model. Calculating these parameters using a FISH code function in FLAC3D with the suggested formulas by FLAC3D manual resolves the numerical instability issue. The FLAC3D manual suggests calculating the 𝐻𝑐 from the coefficient of lateral earth pressure at rest for normally consolidated soil (𝐾0,𝑁𝐶 ) and tangent oedometer stiffness at the 𝑟𝑒𝑓 reference pressure (𝐸𝑜𝑒𝑑 ) using the following relationship: 𝐻𝑐 = 𝑘𝐾𝑝𝑝 𝐾𝑝𝑝 = 𝐾1 𝐾2 𝐾1 − 𝐾2 (4.42) (4.43) 𝑟𝑒𝑓 𝐸𝑢𝑟 𝐾1 = 3(1 − 2𝜈) 𝐾2 = 1 𝑟𝑒𝑓 𝐸 (1 + 2𝐾0,𝑁𝐶 ) 3 𝑜𝑒𝑑 where: 𝑘 = Correction factor with a typical value of 0 < 𝑘 ≤ 1 𝑟𝑒𝑓 𝐸𝑜𝑒𝑑 = tangent oedometer stiffness at the reference pressure (4.44) (4.45) 55 𝐾0,𝑁𝐶 = Coefficient of lateral earth pressure at rest for normally consolidated soil with a maximum value of 𝜈 ⁄(1 − 𝜈) The 𝐾p , 𝐾1 , and 𝐾2 are intermediary variables. The PH elastic formulation that is discussed so far assumes the strain is very small. However, the soil stiffness is strain-dependent, and the elastic modulus decreases as the shear strain increase. The “small-strain stiffness” option in FLAC3D reduces shear modules nonlinearly with increasing shear strain. The secant shear modulus and tangent shear models are calculated by: 𝐺𝑠𝑒𝑐 = 𝐺0 0.385𝛾 𝑡 1+ 𝑡 𝛾70 (4.46) 𝐺0 𝐺𝑡𝑎𝑛 = 2 (1 + 0.385𝛾 𝑡 ) 𝑡 𝛾70 (4.47) where: 𝛾 𝑡 = Shear strain 𝑡 𝛾70 = The reference shear strain corresponding to the secant shear modulus with the value of 72.2% of the initial shear modulus 𝐺0 = Initial shear modulus (at very small strain) 𝐺𝑠𝑒𝑐 = Secant shear modulus 𝐺𝑡𝑎𝑛 = Tangent shear models The FLAC3D requires two properties for its “small-strain stiffness” option: The 𝑟𝑒𝑓 𝑡 reference shear strain (𝛾70 ) and initial Young’s modulus at reference pressure (𝐸0 ). The stress-dependent initial Young’s modulus is calculated using the following relationship: 56 𝑚 𝐸0 = 𝑟𝑒𝑓 𝐸0 𝑐 𝑐𝑜𝑡𝜙 − 𝜎3′ ( ) 𝑐 𝑐𝑜𝑡𝜙 + 𝑝𝑟𝑒𝑓 (4.48) where: 𝐸0 = Initial Young’s modulus (at very small strain) 𝑟𝑒𝑓 𝐸0 = Initial Young’s modulus (at very small strain) at the reference pressure And, the initial shear modulus is calculated by: 𝐺0 = 𝐸0 2(1 + 𝜈) (4.49) 𝑡 Typical values for the 𝛾70 are listed in Table 4.3 (Itasca Consulting Group 2019). The normalized shear modulus degradation curves (Equations 4.46 and 4.47) are plotted in 𝑡 Figure 4.20 for normalized values of 𝛾 𝑡 /𝛾70 ,and the absolute value of 𝛾 𝑡 with 𝑡 recommended range for 𝛾70 . 4.3.2.4.1 FLAC3D parameterization The parameterization of the PH model is done with a FISH function in FLAC3D that converts the constitutive model of each zone from MC to PH. The original FISH function that was provided by Mr. Augusto Lucarelli, Principal Engineer at ITASCA Consulting Group Inc., had been modified for the FLAC3D model in this research. The PH model assignment involves reading the MC model material properties for each zone, getting the extra parameters from FISH tables or direct input, calculating the PH model input parameters for each zone, and assigning the parameters. The FLAC3D model needs to be solved in the initial state before calling the FISH code to convert the constitutive model. The FLAC3D PH model can be parameterized directly as well. However, converting the MC model to the PH model with FISH coding ensures the parameters 57 comparability PH and MC models. In addition, the relatively new PH constitutive model in FLAC3D lacks the robustness of the well-established MC model. And, using the FISH coding to change the constitutive model of a FLAC3D model in an equilibrium state would reduce the PH model numerical issues significantly. The elastic modulus and the yield surface of the PH model are stress-dependent (see Equations 4.31, 4.32, 4.33, and 4.35). And, the PH model requires the initial principal effective stresses as input. The FISH function reads the principal effective stress for each zone and assigns it as an input parameter for the PH model. Mohr-Coulomb friction angle (∅′ ), cohesion intercept (𝑐 ′ ), dilatation angle (𝜓), and Poisson’s ratio (𝜈) are imported directly from the MC model in the same manner as the effective principal stresses. The reference stiffnesses are estimated by following the procedure described by Lucarelli and Cheng (2016) to parameterize the PH model from in situ tests. The Initial shear modulus at very small strain (𝐺0 ) is calculated by: 𝐺0 = 𝜌𝑉𝑠2 = 𝛾 2 𝑉 𝑔 𝑠 (4.50) where: 𝜌 = Soil density 𝑔 = Earth gravitational acceleration 𝑉𝑠 = Shear wave velocity The soil unit weight is calculated with Equations 4.7 and 4.8 (see Figure 4.13). The shear wave velocity is estimated using Mayne and Rix CPR correlation (Mayne 2007) for clay: 𝑉𝑠 = 1.75(𝑞𝑡 )0.627 (4.51) The initial Young’s modulus at very small strain (𝐸0 ) is calculated by Equation 58 4.49. The 𝐸𝑢𝑟 and 𝐸50 are estimated as a fraction of the 𝐸0 by: 𝐸𝑢𝑟 = 𝐸0 𝑑𝑖𝑣𝑢𝑟 (4.52) 𝐸50 = 𝐸𝑢𝑟 𝑑𝑖𝑣50 (4.53) where: 𝑑𝑖𝑣𝑢𝑟 = A coefficient ranging from 2 to 5 and the most common value of 3 𝑑𝑖𝑣50 = A coefficient ranging from 2 to 5 and the most common value of 3 The recommended value of 3 is chosen for 𝑑𝑖𝑣𝑢𝑟 and 𝑑𝑖𝑣50 coefficients in this research. By choosing the atmospheric pressure as the reference pressure (𝑝𝑟𝑒𝑓 = 𝑃𝑎 ), the 𝑟𝑒𝑓 𝑟𝑒𝑓 𝐸50 and 𝐸𝑢𝑟 are calculated from equations 4.31 and 4.32 and assigned as an input for the PH model. The tangent oedometer stiffness at the reference pressure (𝐸𝑜𝑒𝑑 ) usually ranges 𝑟𝑒𝑓 from 0.8 to 1.2 (with 1 as the most common value) times the 𝐸50 (. And, the 𝐸𝑜𝑒𝑑 is chosen 𝑟𝑒𝑓 equal to the 𝐸50 value in the FLAC3D model. The failure ratio (𝑅𝑓 ) is set the recommended value of 0.9 (see Equation 4.28), and the exponent “m” is set to 0.99 for matrix soil, which mostly consists of clay, and 0.48 for the RAP material based on the Duncan et al. (1980) recommended value for similar material. The OCR calculated from Equation 4.11 (see Figure 4.14) and assigned to the PH model to initialize the cap yield surface. The coefficient of lateral earth pressure at rest for normally consolidated soil (𝐾0,𝑁𝐶 ) is calculated by Jacky’s (1948) method and is limited to the maximum allowable value for the PH model by: 𝜈 )] 1−𝜈 𝐾0,𝑁𝐶 = 𝑚𝑖𝑛 [(1 − 𝑠𝑖𝑛𝜙 ′ ) 𝑎𝑛𝑑 ( (4.54) FLAC3D manual recommends a range of 0.6 to 2 for the parameter 𝛼. The 59 average value of 1.3 is selected for this parameter in the FLAC3D model. The value of 𝑡 the reference shear strain (𝛾70 ) for shear module degradation curve is selected from the recommended range by FLAC3D manual (see Table 4.3). 4.4 Rammed Aggregate Piers Construction 4.4.1 K0-Blade Test Stress Distribution The first technique is a simple initialization of the horizontal stress within the matrix soil, and RAP element to match the measured values in the field after the RAPs construction. The coefficient of the lateral earth pressure measured at 1 ft (0.30 m) from the RAP interface (2.5 ft (0.76 m) from the center of the RAP element) is used to initialize the effective horizontal stress in FLAC3D. This technique does not push the soil to failure, and soil zones stress state in FLAC3D remain elastic after the RAP construction. Also, it does not consider the vertical stress increase due to RAP construction in modeling. This RAP construction technique was implemented to verify the overall FLAC3D results at earlier modeling stages and to debug the FLAC3D code. Only the load-displacement curve of the FLAC3D uplift test with this RAP construction technique is shown in Chapter 5. 4.4.2 Lift Compaction The lift compaction technique simulates the compaction of RAP lifts with the application of static compaction stress on top of each lift (LaMeres 2005). The first RAP lift is placed at the bottom of the borehole, and increasing compressive stress is applied on top of the RAP lift to compact it statically and to build up the confining stress. The applied stress is increased gradually from zero to a predetermined maximum value to minimize the 60 inertial effect in the FLAC3D simulation. After completion of the first lift compaction, the compressive stress is removed, and the second lift is placed on top of it. The lift placement and compaction are repeated to build the entire RAP. The model is solved to reach equilibrium after the final RAP lift placement and compaction. LaMeres (2005) compared the effective lateral stress at 0.5 ft (0.15 m) from the RAP interface to the theoretical effective lateral stress proposed by Lawton (2000) and Hsu (2000) to determine the optimal lift compaction stress. He proposed the lift compaction pressure of 25.5 ksf (1.22 MPa). In this thesis, the coefficient of the lateral earth pressure from FLAC3D RAP construction simulation were compared directly to the K0-Blade test measurement to determine the optimal lift compaction pressure. It was determined that the compaction stress of 7.0 ksf (0.34 MPa) for the 3 ft (0.30 m) RAP and 22.0 ksf (1.05 MPa) for 6 ft (1.83 m), 9ft (2.74 m), 12 ft (3.66 m), and 15 ft (4.57 m) RAP results in a coefficient of lateral stress that better match the filed data. The results are shown and discussed in Chapter 5. 4.4.3 Cavity Expansion The Cavity expansion technique is based on the cavity expansion method that was proposed by Pham and White (2007) for numerical modeling of RAPs based on their FSTs measurement and observation (White et al. 2007). In this technique, the shaft and the bottom of the excavated borehole are pushed outward and downward simultaneously to expand the cavity and increase the radial stress and vertical stress in the matrix soil. After reaching a predetermined radial and vertical strain at the shaft and bottom of the borehole, the cavity expansion is stopped, all the displacement restraints on the shaft and bottom of the borehole are removed, the RAP material is placed in the expanded borehole, and the 61 model is solved to reach equilibrium. Pham and White (2007) expanded the shaft and bottom of the borehole to %5 and 10% of the RAP nominal diameter respectively to increase the lateral stress and simulation of the RAP construction process. In this research, the bottom displacement was set to twice the shaft displacement (the same ratio that Pham and White (2007) observed). The RAP material was placed in the cavity after cavity expansion, and the model was solved to reach equilibrium. The cavity expansion in this technique represents the average lateral expansion of the entire length of the RAP. The matrix soil properties are different at each depth, and the same amount of lift compaction effort for each lift in the field could expand the cavity differently at each depth. To account for soil property variation with depth, the amount of cavity expansion was evaluated separately for each RAP. For each RAP model, the cavity was expanded with a few different values for the lateral displacement, and the coefficient of lateral earth pressures was calculated in FLAC3D and compared with the measured values by the K0-Blade test at 1 ft (0.30 m) away from the RAP interface. The lateral displacement that resulted in the best match between simulated and measured coefficient of lateral stress was selected for RAP construction. The lateral borehole expansion displacement is set to 0.33 in. (8.4 mm) for the 3 ft (0.30 m) RAP, 0.21 in. (5.3 mm) for the 6 ft (1.83 m) RAP, 0.27 in. (6.9 mm) for the 9 ft (2.74 m) RAP, 0.27 in. (6.9 mm) for the 12 ft (3.66 m) RAP, and 0.21 in. (5.3 mm) for the 15 ft (4.57 m) RAP. The results are shown and discussed in Chapter 5. 62 4.4.4 Cavity Expansion with Variable Displacement This technique expands the cavity by applying variable velocity (displacement increment per step) to the borehole shaft and bottom. The magnitude of the lateral displacement is variable through the depth of the borehole and throughout the simulation steps. The cavity expansion is done in steps, and a linearization scheme is used to estimate the new lateral displacement magnitude and pattern for each step. The Python programming language, which is embedded in version 7.0 of FLAC3D, is used to implement a subroutine to simulate the RAP construction by cavity expansion with variable displacement technique. The linear algebra package of the “NumPy” library in Python was essential in implementing this technique. The embedded Python program in FLAC3D provides access to FISH function and FLAC3D command through the “Itasca” library in Python. A linear approximation of a nonlinear function can be estimated by using the firstorder term of the Taylor series around a specific point. For a single variable function of 𝑦 = 𝑓(𝑥), the linear approximation of the dependent variable increment is: Δ𝑦 = 𝑦 − 𝑦0 ≈ 𝑓 ′ (𝑥)|𝑥0 × (𝑥 − 𝑥0 ) where: 𝑥 = A scalar independent variable 𝑥0 = The point of interest of the independent variable 𝑦 = A scalar dependent variable 𝑦0 = The dependent variable value at 𝑥0 Δ𝑦 = The independent variable increment 𝑓 ′ (𝑥)|𝑥0 = The first derivative of 𝑦 with respect to 𝑥 evaluated at 𝑥0 (4.55) 63 The generalized form of Equation 4.55 for a multivariable function of 𝑦⃗ = 𝑓⃗(𝑥⃗) is: Δ𝑦⃗ ≈ [∇𝑓⃗|𝑥⃗ ] . Δ𝑥⃗ 0 (4.56) where: 𝑥⃗ = A vector of independent variables with the size of (𝑟) 𝑥⃗0 = The independent variables vector at points of interest 𝑦⃗ = A vector of dependent variables with the size of (𝑡) 𝑓⃗ = A multivariable function [∇𝑓⃗|𝑥⃗ ] = The Jacobian matrix evaluated at 𝑥⃗0 with the size of (𝑡 × 𝑟) 0 The dot (.) in Equation 4.55 represents the inner product of the Jacobian matrix to the Δ𝑥⃗ vector. The Jacobian matrix is defined by: 𝜕𝑓1 𝜕𝑥1 [ 𝐽 ] = [∇𝑓⃗] = ⋮ 𝜕𝑓𝑡 [𝜕𝑥1 𝜕𝑓1 𝜕𝑥𝑟 ⋱ ⋮ 𝜕𝑓𝑡 ⋯ 𝜕𝑥𝑟 ] ⋯ (4.57) For a specific RAP model in FLAC3D, we can assume a nonlinear relationship between the lateral displacement during the borehole expansion (𝛿⃗) and the effective radial stress in the matrix soil as a close distance from the borehole shaft (𝑃⃗⃗ ). 𝑃⃗⃗ = 𝑓⃗(𝛿⃗) (4.58) where: 𝛿⃗ = A vector consisting of shaft’s lateral displacement 𝑃⃗⃗ = A vector consisting of the effective radial stress in the matrix soil at a close distance from the borehole shaft The lateral displacement that is applied at the borehole shaft is assumed to be 64 constant within the depth of each RAP lift. Therefore, the size of the vector 𝛿⃗ is equal to the number of RAP lifts (𝑛𝑙𝑖𝑓𝑡). The vertical displacement of the bottom of the borehole was constraint to twice the magnitude of the lateral displacement of the first lift. This constraint was applied based on White et al. (2007) observation of the borehole expansion. The effective lateral stress at the radial distance of 1 ft (0.30 m) from the borehole shaft (2.5 ft (0.76 m) from the center of the RAP element) is averaged within the depth of each RAP lift. These averaged values are the components of the vector 𝑃⃗⃗. The vector 𝑃⃗⃗ has the same number of the components as the vector 𝛿⃗ (the number of RAP lifts (𝑛𝑙𝑖𝑓𝑡)). Because of the highly nonlinear behavior of the soil, Equation 4.58 should be linearized in several steps. By replacing Equation 4.58 in Equation 4.56 for each step, we have: Δ𝑃⃗⃗𝑛 ≈ [∇𝑓⃗|⃗𝛿⃗𝑛 ] . Δ𝛿⃗𝑛 = [𝐽𝑛 ] . Δ𝛿⃗𝑛 (4.59) where: 𝛿⃗𝑛 = Lateral displacement vector at step 𝑛 with the size of (𝑛𝑙𝑖𝑓𝑡) 𝑃⃗⃗𝑛 = Effective radial stress vector at step 𝑛 with the size of (𝑛𝑙𝑖𝑓𝑡) Δ𝛿⃗𝑛 = Lateral displacement increment vector at step 𝑛, (𝛿⃗𝑛+1 − 𝛿⃗𝑛 ) Δ𝑃⃗⃗𝑛 = Effective radial stress increment vector at step 𝑛, (𝑃⃗⃗ 𝑛+1 − 𝑃⃗⃗𝑛 ) [𝐽𝑛 ] = Jacobian matrix at step 𝑛 with the size of (𝑛𝑙𝑖𝑓𝑡 × 𝑛𝑙𝑖𝑓𝑡) We can find the components of the Jacobian matrix by performing a series of FLAC3D simulation that pushes the borehole shaft one lift at a time (see Figure 4.21). If the cavity is being expanded at only the ith lift, the lateral displacement vector would be: 65 0 ← 0 ⋮ ⋮ 0 ←𝑖−1 Δ𝛿⃗𝑖𝑛 = 𝐷𝑖𝑛 1 ← 𝑖 0 ←𝑖+1 ⋮ ⋮ {0} ← 𝑛𝑙𝑖𝑓𝑡 (4.60) where: th Δ𝛿⃗𝑖𝑛 = Lateral displacement increment vector with only the i component having a nonzero value, (𝛿⃗𝑖𝑛 − 𝛿⃗𝑛 ) 𝐷𝑖𝑛 = Lateral displacement magnitude of the cavity expansion at the ith lift in step 𝑛 during single lift expansion. By replacing Equation 4.60 in into Equation 4.59, we have: 𝑛 𝐽0,𝑖 ⋮ 𝑛 𝐽𝑖−1,𝑖 𝑛 𝐽𝑖,𝑖 Δ𝑃⃗⃗𝑖𝑛 ≈ [𝐽𝑛 ] . Δ𝛿⃗𝑖𝑛 = 𝐷𝑖𝑛 𝑛 𝐽𝑖+1,𝑖 ⋮ 𝑛 {𝐽𝑛𝑙𝑖𝑓𝑡,𝑖 } (4.61) where: Δ𝑃⃗⃗𝑖𝑛 = Effective lateral stress increment vector after expanding the cavity at the ith lift, (𝑃⃗⃗𝑖𝑛 − 𝑃⃗⃗ 𝑛 ) The right-hand side of Equation 4.61 is the ith column of the Jacobian matrix multiplied by the 𝐷𝑖𝑛 . The Δ𝑃⃗⃗𝑖𝑛 and 𝐷𝑖𝑛 are known, and the ith column of the Jacobian matrix is estimated by: 66 𝑛 𝐽0,𝑖 ⋮ 𝑛 𝐽𝑖−1,𝑖 1 𝑛 𝐽𝑖,𝑖 ≈ 𝑛 Δ𝑃⃗⃗𝑖𝑛 𝐷𝑖 𝑛 𝐽𝑖+1,𝑖 ⋮ 𝑛 𝐽 { 𝑛𝑙𝑖𝑓𝑡,𝑖 } (4.62) By repeating this process for all the RAP lifts, we can assemble the Jacobian matrix. However, the chosen value for the 𝐷𝑖𝑛 could affect the Jacobian matrix due to the nonlinear behavior of the soil. The value of the 𝐷𝑖𝑛 should be chosen close to the value of the ith component of the lateral displacement increment vector at the current step (ith component of Δ𝛿⃗𝑛 ). However, the Δ𝛿⃗𝑛 vector is unknown. Therefore, the 𝐷𝑖𝑛 is estimated by the ith component of the lateral displacement increment vector at the previous step (ith component of Δ𝛿⃗𝑛−1). The displacement pattern to expand the cavity from step 𝑛 to step 𝑛 + 1 is estimated by solving Equation 4.59 for the Δ𝛿⃗𝑛 . However, the 𝑃⃗⃗𝑛+1 vector is unknown. Therefore, 𝑛 the Δ𝑃⃗⃗𝑛 is replaced with Δ𝑃⃗⃗𝑡𝑎𝑟𝑔𝑒𝑡 and Equation 4.59 is rearranged to: 𝑛 𝑛 Δ𝛿⃗𝑡𝑎𝑟𝑔𝑒𝑡 ≈ [𝐽𝑛 ]−1 . Δ𝑃⃗⃗𝑡𝑎𝑟𝑔𝑒𝑡 = [𝐽𝑛 ]−1 . ( 𝑃⃗⃗𝑡𝑎𝑟𝑔𝑒𝑡 − 𝑃⃗⃗𝑛 ) (4.63) where: 𝑃⃗⃗𝑡𝑎𝑟𝑔𝑒𝑡 = Effective radial stress vector at 1 ft (0.30 m) from the borehole shaft (2.5 ft (0.76 m) from the center of the RAP element) measured with K0-Blade test after the RAP construction [𝐽𝑛 ]−1 = The inverse of the Jacobian matrix 𝑛 The Δ𝛿⃗𝑡𝑎𝑟𝑔𝑒𝑡 vector is a linear estimation of the amount of displacement that is needed to increase the effective lateral stress at step 𝑛 to the measured value at the end of 67 𝑛 the RAP construction. However, if we used the calculated Δ𝛿⃗𝑡𝑎𝑟𝑔𝑒𝑡 vector, we would have assumed a linear soil behavior for the entire process. Therefore, only the pattern of the 𝑛 vector Δ𝛿⃗𝑡𝑎𝑟𝑔𝑒𝑡 is used for cavity expansion (not its magnitude). The new displacement pattern is defined by: ⃗⃗⃗⃗⃗⃗⃗ 𝑛 = 𝑝𝑎𝑡 𝑛 Δ𝛿⃗𝑡𝑎𝑟𝑔𝑒𝑡 𝑛 ‖Δ𝛿⃗𝑡𝑎𝑟𝑔𝑒𝑡 ‖ (4.64) where: ⃗⃗⃗⃗⃗⃗⃗ 𝑝𝑎𝑡 𝑛 = The new displacement pattern 𝑛 𝑛 ‖Δ𝛿⃗𝑡𝑎𝑟𝑔𝑒𝑡 ‖ = Magnitude (l2-norm) of the Δ𝛿⃗𝑡𝑎𝑟𝑔𝑒𝑡 vector The Lateral displacement increment vector at step 𝑛 (Δ𝛿⃗𝑛 ) is calculated by multiplying the new displacement pattern with a chosen magnitude as follows: 𝑛 ⃗⃗⃗⃗⃗⃗⃗ 𝑛 Δ𝛿⃗𝑛 = 𝛼𝑠𝑐 × 𝑝𝑎𝑡 (4.65) where: 𝑛 𝛼𝑠𝑐 = The scaling factor for the displacement pattern of step 𝑛 To adjust the scaling factor that determines how far the cavity would expand at each step, we need to have an estimate of the increase in effective lateral stress. A normalized measure of effective lateral stress is used for this purpose. The Normalized Lateral Force (NLF) is defined by: ∑ 𝑃⃗⃗𝑛 − ∑ 𝑃⃗⃗ 0 𝑁𝐿𝐹 = ∑ 𝑃⃗⃗𝑡𝑎𝑟𝑔𝑒𝑡 − ∑ 𝑃⃗⃗0 𝑛 where: ∑ 𝑃⃗⃗𝑛 = The sum of the effective lateral stress at step 𝑛 (4.66) 68 ∑ 𝑃⃗⃗ 0 = The sum of the initial effective lateral stress at step 𝑛 ∑ 𝑃⃗⃗𝑡𝑎𝑟𝑔𝑒𝑡 = The sum of the measured effective lateral stress at the end of RAP construction The sum of the effective lateral stress multiplied by the area of each zone normal to the radial direction gives the sum of the effective lateral force in FLAC3D. In the FLAC3D model, the zones are discretized evenly in the z-direction from the ground surface to the bottom of the RAP element. Also, the lateral stress is constant at each depth. Therefore, the sum of the zones area in the radial direction at each depth (surface area of the side a cylinder with a height of a FLAC3D zone) can be factored out of the force summation. Because of the normalization of the NLF equation, the constant area multiplier cancels out from the Numerator and Denominator of Equation 4.66. Therefore, the 𝑁𝐿𝐹 𝑛 at each step describes the increase of the effective lateral force from the beginning of the RAP construction normalized by the desired increase in the effective lateral force (the increase of the effective lateral force by the RAP construction process). 𝑛 The scaling factor 𝛼𝑠𝑐 is estimated by: 𝑛 𝛼𝑠𝑐 = 𝑔𝑠𝑐 (𝑁𝐿𝐹 𝑛 ) 𝑛−1 × 𝛼𝑠𝑐 𝑁𝐿𝐹 𝑛 − 𝑁𝐿𝐹 𝑛−1 (4.67) where: 𝑛−1 𝛼𝑠𝑐 = The scaling factor at the previous step (step 𝑛 − 1) 𝑔𝑠𝑐 (𝑁𝐿𝐹 𝑛 ) = A function that determines the desired increase in 𝑁𝐿𝐹 for the next step Let’s examine Equations 4.66 and 4.67 with a simple example. Let’s assume that the sum of the initial effective lateral stress (∑ 𝑃⃗⃗0 ) is 20 ksf, the sum of the measured effective lateral stress at the end of RAP construction (∑ 𝑃⃗⃗𝑡𝑎𝑟𝑔𝑒𝑡 ) is 120 ksf, the sum of 69 the effective lateral stress at the previous step (∑ 𝑃⃗⃗ 𝑛−1 ) is 40 ksf, and the sum of the effective lateral stress at the current step (∑ 𝑃⃗⃗𝑛 ) is 65 ksf. And, we want to increase the effective lateral force at each step by 20% of the total increase in the effective lateral force after RAP construction. Therefore, the function 𝑔𝑠𝑐 (𝑁𝐿𝐹 𝑛 ) would have a constant value of 0.2 for all the 𝑁𝐿𝐹 𝑛 values. Using Equation 4.66, 𝑁𝐿𝐹 𝑛−1 = 0.2 and 𝑁𝐿𝐹 𝑛 = 0.45, which means that the cavity expansion at step 𝑛 − 1 increased the effective lateral force 25% of the total increase of the effective lateral force after the RAP construction process. 0.2 𝑛) As a result, Equation 4.46 reduces the step 𝑛 scaling factor (𝛼𝑠𝑐 to 80% (0.25 = 0.8) of 𝑛−1 ). the step 𝑛 − 1 scaling factor (𝛼𝑠𝑐 The soil behavior shows more nonlinearity as the cavity expands. And, the efficiency of the linearization techniques is directly affected by soil nonlinearity. Therefore, the desired increase in the effective lateral force should decrease as the cavity expands. The following linear function is chosen to decrease the desired effective lateral force increment at each step from 25% to 2% of the total increase of the effective lateral force by the RAP construction process: 𝑔𝑠𝑐 (𝑁𝐿𝐹 𝑛 ) = 0.25 − 0.23 × 𝑁𝐿𝐹 𝑛 (4.68) The cavity expansion by variable displacement technique is implemented with the following steps: - 0 ). ⃗⃗⃗⃗⃗⃗⃗ 0 ) and scaling factor (𝛼𝑠𝑐 Choose an initial displacement pattern (𝑝𝑎𝑡 - Expand the cavity in FLAC3D by the Δ𝛿⃗0 (Equation 4.65) and obtain 𝛿⃗1 and 𝑃⃗⃗1 - Expand the cavity one lift at a time and assemble the Jacobian matrix (Equation 4.62) 70 - Calculate the new displacement pattern (Equations 4.63 and 4.64) - Calculate the new scaling factor (Equations 4.66 and 4.67) - Calculate the new lateral displacement increment vector (Δ𝛿⃗1 ) (Equation 4.65) - Repeat the cavity expansion process until the 𝑁𝐿𝐹 close to unity is achieved We can evaluate the effectiveness of this technique by comparing the effective radial stress from FLAC3D at each step with the measured values at the filed. In addition to a direct comparison of the two vectors in a plot, we can use the cosine of the virtual angle between two vectors as a quantitative measure of vector alignment as follows: 𝑉𝑒𝑐𝑡𝑜𝑟 𝐴𝑙𝑖𝑔𝑛𝑚𝑒𝑛𝑡 = 𝑐𝑜𝑠(𝜃𝑥,𝑦 ) = 𝑥⃗. 𝑦⃗ ‖𝑥⃗‖‖𝑦⃗‖ (4.69) where: 𝜃𝑥,𝑦 = The virtual angle between two vectors 𝑥⃗ and 𝑦⃗ For 3-dimensional position vectors 𝑥⃗ and 𝑦⃗, Equation 4.69 calculates the cosine of the actual angle between the two vectors. The vector alignment ranges from -1 to 1. The 1 and 1 values represent a perfect vector alignment (vectors have the same direction), and 0 vector alignment represents no alignments (vectors are perpendicular).The term “virtual” is used here to refer to higher dimension vectors such as the effective radial stress vector at step 𝑛 (𝑃⃗⃗𝑛 ). The vector alignment for the 𝑃⃗⃗𝑛 is defined as: (𝑉𝑒𝑐𝑡𝑜𝑟 𝐴𝑙𝑖𝑔𝑛𝑚𝑒𝑛𝑡)𝑃⃗⃗𝑛 = 𝑃⃗⃗𝑛 . 𝑃⃗⃗𝑡𝑎𝑟𝑔𝑒𝑡 ‖𝑃⃗⃗𝑛 ‖‖𝑃⃗⃗𝑡𝑎𝑟𝑔𝑒𝑡 ‖ (4.70) Equation 4.70 gives a quantitative measure of alignment between the effective radial stress vector at step 𝑛 with the effective radial stress vector measured in the field. Another vector alignment is defined by: 71 (𝑉𝑒𝑐𝑡𝑜𝑟 𝐴𝑙𝑖𝑔𝑛𝑚𝑒𝑛𝑡)𝑙𝑖𝑛𝑛 = (𝑃⃗⃗𝑛+1 − 𝑃⃗⃗𝑛 ). ( 𝑃⃗⃗𝑡𝑎𝑟𝑔𝑒𝑡 − 𝑃⃗⃗𝑛 ) ‖𝑃⃗⃗𝑛+1 − 𝑃⃗⃗𝑛 ‖‖𝑃⃗⃗𝑡𝑎𝑟𝑔𝑒𝑡 − 𝑃⃗⃗𝑛 ‖ (4.71) Equation 4.71 directly measures the effectiveness of the linearization Equation 4.63 by calculating the alignment between the vector that linearization is solving for ( 𝑃⃗⃗𝑡𝑎𝑟𝑔𝑒𝑡 − 𝑃⃗⃗𝑛 ) to get the new displacement pattern and the outcome vector (𝑃⃗⃗ 𝑛+1 − 𝑃⃗⃗𝑛 ) of the cavity expansion with the calculated displacement pattern. 4.5 Uplift Test The simulation of the uplift test is performed with a force control method. An incremental load is applied at the top of the threaded bar, and model is solved to reach equilibrium. The force control approach is similar to the loading procedure of the FSTs. The force control method is not able to capture the possible softening and drop of the applied load. However, the result of the FSTs and the preliminary FLAC3D models did not show any sign of softening during the uplift load. The main advantage of the force control method to the displacement control method in FLAC3D is reaching the equilibrium at the end of each load increment. Usually, modeling the FEM structural elements embedded in FDM soil zones in FLAC3D results in large unbalanced internal forces that requires a higher number of iterations to equalize. The displacement control method applies a velocity (displacement increment per step in static analysis) on top of the thread bar and cycles through each step to increase the displacement of the top of the threaded bar gradually. This process results in large unbalanced internal forces which decrease the accuracy of the numerical simulation. 72 4.5 Verification Models The FLAC3D code is verified against the well-established relationship to estimate the uplift capacity of the vertical anchors. The RAP construction subroutine is disabled for these simulations. Similar to the primary simulation for the RAP FTSs, only one-quarter of the anchor and native soil is modeled. The anchor is modeled with a very stiff uplift plate and cable to minimize the deflection of structural elements. The interface friction angle between the uplift plate and soil is set to the internal friction angle of the soil. And, the interface friction angle between cable and the soil is set to zero to ensure the complete transfer of the uplift load to the uplift plate. Soil is modeled as a homogenous medium in verification models. Soil unit weight is 110 pcf (17.3 KN/m3), and the groundwater table is assumed to be very deep. Three anchors with 3 ft (0.91 m) diameters were modeled and analyzed in cohesionless soil. The first anchor was embedded 3 ft (0.91 m) deep in a granular soil with a friction angle of 35° and Young’s modulus of 400 ksf (19.2 MPa). The second anchor model has the same depth (3 ft) and friction angle (35°) with Young’s modulus of 800 ksf (38.3 MPa). And, the third anchor is 9 ft deep (2.74 m) embedded in granular soil with a friction angle of 30 ° and Young’s modulus of 800 ksf (38.3 MPa). 73 Table 4.1. Model size and grid statistic Model Total zones Model radius Model height Liner elements Cable elements 3 ft RAP 12,636 10 ft 15 ft 90 6 6 ft RAP 16,848 10 ft 18 ft 90 12 9 ft RAP 21,060 10 ft 21 ft 90 18 12 ft RAP 25,272 10 ft 24 ft 90 24 15 ft RAP 29,484 10 ft 27 ft 90 30 Table 4.2. The typical range for Poisson ratio, 𝜈 (Modifed from Bowles 1996) Type of soil Poisson ratio, 𝝂 Saturated Clay 0.4 – 0.5 Unsaturated Clay 0.1 – 0.3 Sandy clay 0.3 – 0.35 Silt 0.3 – 0.35 Sand, and gravelly sand 0.3 – 0.4 74 𝑡 Table 4.3. The typical range for 𝛾70 (Modifed from Itasca Consulting Group 2019) Type of soil 𝜸𝒕𝟕𝟎 Sand: At depth 0 to 20 ft (0 to 6 m) 1.1 × 10−4 At depth 20 to 49 ft (6 to 15 m) 1.6 × 10−4 At depth 49 to 121 ft (15 to 37 m) 2.3 × 10−4 At depth 121 to 249 ft (37 to 76 m) 3.3 × 10−4 At depth 249 to 492 ft (76 to 150 m) 4.5 × 10−4 At depth 492 to 984 ft (150 to 300 m) 7.0 × 10−4 Plastic Index = 30 3.5 × 10−4 Plastic Index = 50 7.7 × 10−4 Plastic Index = 100 1.8 × 10−3 Clay: 75 threaded Figure 4.1. Flowchart for the RAP uplift simulation in FLAC3D 76 Figure 4.2. Flowchart for the RAP construction in FLAC3D, K0-Blade Stress technique 77 Figure 4.3. Flowchart for the RAP construction in FLAC3D, Lift Compaction technique 78 Figure 4.4. Flowchart for the RAP construction in FLAC3D, Cavity Expansion technique 79 Figure 4.5. Flowchart for the RAP construction in FLAC3D, Cavity Expansion with Variable Displacement technique 80 Hole for threaded bar Threaded bar Quarter of the uplift plate is modeled in FLAC3D Uplift plate Matrix soil boundary during the uplift simulation Planes of symmetry (a) Matrix soil boundary before uplift simulation (b) Figure 4.6. Schematic RAP model in FLAC3D. (a) Uplift Plate (b) Elevation view of the RAP model in FLAC3D 81 Matrix soil boundary RAP and Matrix soil interface a. Bottom of the RAP b. Figure 4.7. Zone generation by extruding an unstructured grid (6 ft RAP model). a. Plan view of the unstructured grid, and b. Elevation view of the grid 82 Figure 4.8. FLAC3D model of the 6 ft RAP 83 a. Duplicated grid points at the bottom of the RAP b. Figure 4.9. Uplift plate modeling. a. Triangle liner elements of the uplift plate, and b. Duplicated grid points at the bottom of the RAP 84 Structural link a. b. Figure 4.10. Uplift plate and threaded rod modeling. a. Connection of the cable element to the uplift plate, and b. Embedded uplift plate and threaded rod (6 ft RAP model) 85 a. b. liner liner c. d. Figure 4.11. FLAC3D liner element. a. The coordinate system and 6 DOFs of the liner element, b. Normal stress versus relative normal displacement, c. shear stress versus relative shear displacement, and d. Shear-strength criterion (Modifed from Itasca Consulting Group 2019) 86 a. b. ‘ c. d. Figure 4.12. FLAC3D cable element. a. The coordinate system and 1 DOF of the cable element, b. Cable element elastic, perfectly plastic material behavior, c. Grout material shear force per unit length versus relative shear displacement, and d. Grout material shear-strength criterion (Modifed from Itasca Consulting Group 2019) 87 Saturated unit weight, 𝛾𝑠𝑎𝑡 (pcf) 90 100 110 120 0 Initial total vertical stress, 𝜎v0 (psf) 1000 2000 3000 4000 0 Adjusted depth below ground surface, z (ft) -5 -10 -15 -20 -25 -30 a. b. Figure 4.13. Native soil saturated unit weight and total vertical stress. a. Saturated unit weight (𝛾𝑠𝑎𝑡 ), and b. Total vertical stress (𝜎𝑣0 ) 88 Initial effective vetical stress, 𝜎'𝑣0 (psf) OCR 0 20 40 60 80 0 1000 2000 0 Adjusted depth below ground surface, z (ft) -5 -10 -15 -20 -25 -30 a. b. Figure 4.14. Native soil elastic saturated unit weight and total vertical stress. a. OCR, and ′ b. Initial effective vertical stress (𝜎𝑣0 ) 89 Coefficient of lateral earth pressure at rest, K0 0 1 2 3 4 Coefficient of lateral earth pressure at rest, K0 5 0 1 2 3 0 Adjusted depth below ground surface, z (ft) -5 -10 -15 -20 -25 K 0 from K0-Blade test K by CPT corrolation 0 (Equation 4.12) K 0 by CPT corrolation 9 (Equation 4.18) Selected K 0 for FLAC3D modeling -30 a. b. Figure 4.15. The coefficient of lateral earth pressure at rest (𝐾0 ). a. Comparison of 𝐾0 estimated by Equations 4.12 and 4.19 with measured values at the field, and b. Selected 𝐾0 to initialize the horizontal stress in the FLAC3D model 90 0.3 Poisson's Ratio, ν 0.4 0.5 0 Young's modulus, Es,OC (ksf) 2000 4000 6000 0 CL, Lean Clay CH, Fat Clay CL, Lean Clay -5 CH, Fat Clay with Sand CL, Lean Clay Adjusted depth below ground surface, z (ft) ML, Sandy Silt CL, Lean Clay CL, Lean Clay with Sand -10 CL-ML, Silty Clay ML, Silt with Sand CL-ML, Silty Clay SC, Clayey Sand CL-ML, Silty Clay with Sand -15 CL, Lean Clay -20 Bonneville Clay -25 -30 a. b. Figure 4.16. Native soil elastic properties for the MC model. a. Poisson ratio (𝜈), and b. Young’s modulus (𝐸𝑠,𝑂𝐶 ) 91 Elastic modulus (ksf) 0 2000 4000 6000 0 Adjusted depth below ground surface, z (ft) -5 -10 -15 -20 -25 Bulk modulus, K Shear modulus, G -30 Figure 4.17. The bulk (K) and shear (G) modulus of native soil 92 0 Effective friction angle, 𝜙 ′ (Degrees) 10 20 30 40 50 0 Effective friction angle, 𝜙 ′ (Degrees) 10 20 30 40 50 0 Adjusted depth below ground surface, z (ft) -5 -10 -15 -20 -25 BST 𝜙′ CPT correlation ⬚ ⬚ Selected for FLAC3D modeling -30 a. b. Figure 4.18. The effective friction angle (𝜙 ′ ). a. Comparison of 𝜙 ′ estimated by Equations 4.14 and 4.17 with measured BST values at the field, and b. Selected 𝜙 ′ for FLAC3D modeling 93 Figure 4.19. Stress-dependency of the 𝐸50 and 𝐸𝑢𝑟 in the PH model. a. relatively low soil 𝑐 𝑐𝑜𝑡𝜙 𝑐 𝑐𝑜𝑡𝜙 strength ( 𝑝𝑟𝑒𝑓 = 0.5), and b. relatively high soil strength ( 𝑝𝑟𝑒𝑓 = 3.0) 94 Figure 4.20. Shear modulus degradation curves. a. Normalized shear modulus versus normalized shear strain, and b. Normalized tangent shear modulus curve with 𝑡 recommended range for 𝛾70 95 Borehole shaft displacement at step nth of the cavity expansion with variable displacement Expanding the shaft at only one lift to assemble Jacobian matrix Measured effective lateral stress in the field 𝛿⃗𝑖𝑛 𝑃⃗⃗𝑖𝑛 𝑃⃗⃗𝑡𝑎𝑟𝑔𝑒𝑡 𝑃⃗⃗ 0 𝑃⃗⃗ 𝑛 𝛿⃗𝑛 Borehole Initial effective lateral stress in FLAC3D Effective lateral stress at step nth of the cavity expansion with variable displacement Effective lateral stress after expanding the ith lift at step n Figure 4.21. The process to find the Jacobian matrix components in cavity expansion with variable displacement technique CHAPTER 5 RESULTS 5.1 FLAC3D Verification Models Figure 5.1a shows the load-displacement for the 3 ft (0.91 m) anchor verification simulations. The two FLAC3D models have the same friction angle of 35º, cohesion of 0, anchor diameter of 3 ft (0.91 m), and the anchor embedment depth of 3 ft (0.91 m). The only difference is Young’s modulus in the two models. The first model has Young’s modulus of 400 ksf (19.2 MPa), and the second model has Young’s modulus of 800 ksf (38.3 MPa). The FLAC3D models have the same uplift capacity of 5.6 kips (24.9 kN). This result is consistent with the Mohr-Coulomb constitutive model assumption, which determines the uplift capacity by the strength parameters, and the elastic properties do not change the ultimate uplift capacity. The displacement of the top of the uplift cable for the model with the 400 ksf (38.3 MPa) Young’s modulus is twice the displacement of the model with Young’s modulus of 800 ksf (38.3 MPa). This ratio is almost constant during the pull-out test. Therefore, the stiffer soil resulted in smaller displacement during the pullout test. The theoretical uplift capacity of the anchor with Mayerhof and Adams’ method and Vesic’s method (Das 1987) are plotted in Figure 5.1a for comparison. The theoretical uplift capacity is 6.1 kips (27.1 kN) by Meyerhof and Adams’ methods and 5.3 kips (23.6 kN) by Vesic’s method. The FLAC3D anchor uplift capacity is 8.6% less than the 97 Meyerhof and Adams’ uplift capacity and 6.3% more than the Vesic’s uplift capacity. Figure 5.1b shows the load-displacement for the 9 ft (2.74 m) anchor verification simulations. The load-displacement shows reasonable behavior. It has a higher ultimate uplift capacity, and displacement at failure compare to the 3 ft (0.91 m) anchor. The FLAC3D model uplift capacity is 46.0 kips (204.6 kN). The numerical model uplift capacity is 17.9% higher than the Meyerhof and Adams’ uplift capacity of 39.0 kips (137.5 kN) and 19.5% higher than the Vesic’s uplift capacity of 38.5 kips (171.3 kN). 5.2 FLAC3D Full-Scale Tests Models 5.2.1 FLAC3D Models with Mohr-Coulomb Constitutive Model 5.2.1.1 Induced Effective Lateral Stress Due to RAP Construction Figures 5.2 to 5.6 show the coefficient of lateral earth pressure versus depth after the RAP construction process by lift compaction and cavity expansion techniques for 3 ft (0.91m) to 15 ft (4.57 m) RAPs. The coefficient of the lateral earth pressure (K) is calculated by obtaining the current effective lateral stress and dividing it by the initial effective vertical stress (not the current effective vertical stress). Therefore, the coefficient of the lateral earth pressure (K), in this section, is a normalized representation of the effective lateral stress. The RAP construction by lift compaction (see Figures 5.2a to 5.6a) was performed for different compaction pressure to determine the compaction pressure that results in the best match between the measured and FLAC3D coefficient of lateral earth pressure (hence, the best match for the effective lateral pressure). By compacting the RAP lifts with higher compaction effort, a higher coefficient of the lateral earth pressure (K) is achieved in the 98 matrix soil below about 3 ft (0.91 m). However, higher compaction effort usually decreases the coefficient of the lateral earth pressure (K) above about 3 ft (0.91). The coefficient of the lateral earth pressure (K) reaches a peak at about 3ft (0.91 m) depth. Similar behavior is reported by LaMeres (2005). The RAP construction by cavity expansion (see Figures 5.2b to 5.6b) was performed with different lateral displacement magnitude to determine the amount of cavity expansion that results in the best match between the measured and FLAC3D coefficient of lateral earth pressure (K). As the cavity expands, the coefficient of the lateral earth pressure (K) increases consistently in the matrix soil below about 3 ft (0.91 m). However, this trend was not observed for soil above about 3 ft (0.91 m). The more the cavity expands, the coefficient of the lateral earth pressure (K) decreases for the 3 ft (0.91m) RAP. The cavity expansion with variable displacement technique expands the cavity with different lateral displacement at each lift. The lateral displacement pattern, effective lateral stress, cavity expansion magnitude, and vector alignment (a quality measure for this technique) are plotted in Figures 5.7 to 5.11 for the 3 ft (0.91 m) to 15 ft (4.57 m) RAPs. The initial lateral displacement pattern (see Figures 5.7a to 5.11a) in step 0 is modified in each step to achieve the best match between the FLAC3D and measured value for effective lateral stress (see Figures 5.7b to 5.11b). The displacement pattern changes substantially in the 1st step. However, the change in the displacement pattern is minimal after the 1st step, and the pattern retains its overall shape. The maximum values for the displacement patterns happen at a depth of about 3 ft (0.91 m) to 6 ft (1.83 m). This depth corresponds to the lowest Young’s modulus for the native soil (see Figure 4.16). The linearization algorithm expands the soil within this depth more to achieve the peak in the 99 effective lateral stress that was measured within this depth (see Figures 5.7b to 5.11b). The step by step increase in the effective lateral stress shows a good match to the target value (measured effective lateral stress) in Figures 5.7b to 5.11b. The magnitude of the cavity 𝑛 expansion is controlled by the displacement pattern magnitude at each step (𝛼𝑠𝑐 ) to ensure a more efficient linearization process (see Figures 5.7c to 5.11c). The displacement pattern magnitude decreases with each step for the 3 ft (0.91 m) RAP. However, it increases to reach a peak and then decreases for other RAPs. The vector alignment that was defined as the cosine of the virtual angle between two vectors in Chapter 4 is plotted in Figures 5.7d to 5.11d. The values of 1 and -1 show a perfect alignment, and 0 represents the worst alignment. The red line shows the alignment of the effective lateral stress at the end of each step with the measured values. The FLAC3D lateral effective stress shows an excellent alignment with the measured values for all the RAPs. The blue line represents a quantitative measure of the linearization process efficiency. This criterion shows a good overall efficiency for the linearization process. However, the efficiency of the linearization decreases at higher steps. As the cavity expands, more soil zones reach the failure (plastic stress state) which increases the nonlinearity in the model. And, the increase in non-linear behavior of the soil decreases the linearization process efficiency. 5.2.1.2 Load-Displacement Curves Figures 5.12 to 5.16 show the load-displacement curves for the FSTs and FLAC3D simulation of the uplift tests. The load-displacement curves show the result of the load increments that resulted in equilibrium after solving the model. The last load increment that does not reach the equilibrium state is higher than the FLAC3D uplift capacity. At this 100 applied uplift load, the FLAC3D model stops the iteration after reaching the maximum number of iteration steps (maximum number of iterations was set to 50,000 steps for FLAC3D models). Although the applied uplift load itself is not useful at this loading increment, the maximum shear strain rate (maximum shear strain increment per iteration step) during the iteration process can be used to have an estimate of the shear failure surface at the end of the pull-out simulation. The maximum shear strain rate happens at the matrix soil adjacent to the RAP element for all the RAP models. The maximum shear strain rate for the 3 ft (0.91 m) and 9 ft (2.74 m) FLAC3D RAP models by the cavity expansion with variable displacement (variable displacement for short) RAP construction technique is shown in Figure 5.17. Because of the force control scheme for the pull-out simulation, the actual FLAC3D model uplift load could be up to one load increment higher than the reported maximum uplift load in this section. The load increment was approximately determined to divide the load-displacement curve into about 25 segments to have an acceptable smooth curve. Therefore, the actual FLAC3D uplift force could be about 4% higher than the reported value. Most of the FLAC3D model underestimated the maximum uplift lift of the FSTs. And, using a smaller load increment for the pull-out simulation could improve the accuracy of the FLAC3D for most RAP models. The load-displacement curves of the FLAC3D simulations include the RAP models with the three main RAP construction techniques of lift compaction, cavity expansion, and cavity expansion with variable displacement (variable displacement for short). In addition, the load-displacement curves for the FLAC3D models without any RAP construction consideration and the FLAC3D model with the simple K0-Blade stress distribution are presented for comparison. By comparing the load-displacement curves of the FLAC3D 101 models with a RAP construction technique to the load-displacement curve of the FLAC3D model without any RAP construction consideration, we can see that the simulation of the RAP construction simulation process increases the predicted uplift capacity for all the RAPs. In addition, RAP construction consideration increases the stiffness of the RAP system and decreases the displacement at each load increment for 6 ft (1.83 m) to 15 ft (4.57 m) RAP models. The 3 ft (0.91 m) RAP model is an exception in this regard. The low overburden pressure and at shallow depth reduces the efficiency of the RAP construction process. The higher compaction effort in the lift compaction technique and more lateral expansion of the borehole in the cavity expansion technique resulted in less lateral confining pressure built-up after the RAP construction process (see Figure 5.2). The load-displacement curves for 3 ft (0.91m) RAP in Figure 5.12 show that the FLAC3D models underpredicted the maximum uplift load and the maximum displacement of the FST tests. The three main RAP construction techniques have similar loaddisplacement curves, and their displacement values are within 0.04 in. (1 mm) of each other up to their maximum uplift load. The lift compaction technique maximum uplift load is 40% less than the FST maximum uplift load. The cavity expansion technique maximum uplift load is 38% less than the FST maximum uplift load. And, the variable displacement technique maximum uplift load is 36% less than the FST maximum uplift load. Figure 5.13 shows the load-displacement curves for the 6 ft (1.83 m) RAP. The three main RAP construction techniques overpredicted the maximum uplift load. The FST maximum uplift load is overpredicted by 17%, 1%, and 12% by the lift compaction, cavity expansion, and variable displacement, respectively. The load-displacement curves of the three main RAP construction techniques have a good match with the FST load- 102 displacement, on average, up to about 26 kips (116 kN). However, they underestimated the displacement at maximum uplift load. The displacement of the FLAC3D models with the three main RAP construction techniques are within about 0.05 in (1 mm) of each other up to about 30 kips (133 kN), but their displacement diverges near the maximum uplift load. The load-displacement curves of the 9 ft (2.74 m) RAP (see Figure 5.14) show an acceptable agreement between the maximum uplift capacity of the three main RAP construction techniques and the FST. The uplift capacity is underestimated by 12%, 16%, and 20% with the lift compaction, cavity expansion, and variable displacement RAP construction techniques, respectively. The displacement of these FLAC3D models have a great match with the FST displacement up to about 26 kips (116 kN), but the FST displacement follows a steeper curve after 26 kips (116 kN), and the FLAC3D models underpredicted the displacement at maximum uplift load. The displacements of the three main RAP construction models are within about 0.05 in. (1 mm) of each other up to about 50 kips (222 kN). Figure 5.15 shows the load-displacement for the 12 ft (3.66 m) RAP. The FST was not performed to completion because of some technical difficulties. The FST was stopped before reaching the maximum uplift capacity. The FST load-displacement curve up to the point that the test was performed is plotted in Figure 5.15. Therefore, the FLAC3D uplift capacity is not compared with the FST values. The FLAC3D models for the three main RAP construction again have good agreement with each other, and their displacement is within 0.1 in. (3 mm) of each other up to 65 kips (289 kN). The load-displacement curves of the FLAC3D models underestimated the uplift capacity of the 15 ft (4.57 m) RAP by 27%, 31%, and 40% for the lift compaction, cavity 103 expansion, and variable displacement respectively (see Figure 5.16). The FLAC3D models with the three main RAP construction techniques have a displacement within 0.1 in (3 mm) of each other up to about 95 kips (423 kN). The variable displacement technique’s displacement increased with a steeper slope after 95 kips (423 kN) up to its maximum uplift load. The ratios of the maximum uplift force FLAC3D models to the FST maximum uplift capacity are given in Table 5.1 and Figure 5.18. The FLAC3D models with the main three RAP construction process show similar overall accuracy in predicting the maximum uplift load. The cavity expansion technique has the overall best match with FST maximum uplift load. Interestingly, the simple horizontal stress initialization to the measured K0Blade stress distribution resulted in a similar prediction of the maximum uplift load. This technique could be used as a simple alternative to the more computationally expensive RAP construction process for estimating the maximum uplift load. By comparing the models with and without RAP construction process, we can see that all of the RAP construction techniques improved the predicted maximum uplift load accuracy. 5.2.1.3 Lateral Earth Pressure and Lateral Displacement During the Pull-out Test The changes in lateral stress (∆𝜎ℎ ) and lateral displacement (𝛿ℎ ) at the matrix soil adjacent to the RAP element during the pull-out test simulation are shown in Figures 5.19 to 5.23 for 3 ft (0.91 m) to 15 ft (4.57 m) RAPs models, respectively. The top plots in Figures 5.19 to 5.23 depict the changes in the lateral stress in the lift compaction, cavity expansion, and variable displacement FLAC3D models at different applied uplift loads. The plots at the bottom of Figures 5.19 to 5.23 show the changes in the lateral displacement 104 in the lift compaction, cavity expansion, and variable displacement FLAC3D models at different applied uplift loads. The lateral stress at the matrix soil decreases in almost all of the RAPs models (13 out of 15 models that are investigated in this section) during the pullout test. The 3 ft (0.91 m) RAP with lift compaction technique and the 6 ft (1.83 m) RAP with the variable displacement techniques are an exception to this trend and show an increase in the lateral stress. The average decrease in the lateral stress within the length of each RAP at the maximum uplift load ranges from 26 psf (1.2 kPa) to 783 psf (37.5 kPa). The amount of the decrease in lateral stress is highly dependent on the RAP construction technique in FLAC3D and the RAP length. Hsu (2000) reported a maximum of 158 psf (7.6 kPa) decrease in the lateral stress of the 3 ft (0.91 m) RAP and a maximum of 122 psf (5.9 kPa) decrease in the lateral stress of the 6 ft (1.83 m) RAP. The measured decrease in the lateral stress falls in the predicted range of the lateral stress changes with the FLAC3D models. The FLAC3D predicted values for the average change in the lateral stress increases as the RAP length increases for the FLAC3D models with cavity expansion techniques and variable displacement techniques. However, the FLAC3D predicted values for the average change in the lateral stress increase to a peak value as the RAP length increases and decreases afterward by an increase in the RAP length for the FLAC3D models with lift compaction. The lateral displacement of the RAP element interface with the matrix soil is evaluated at two locations in the FLAC3D model’s x-y plane. The RAP interface with the matrix soil right in front of the threaded rod location (𝜃 = 45°) are shown with dashed lines. And, the RAP interface with the matrix soil between the threaded rods (𝜃 = 0°) are 105 shown with solid lines (see the plots at the bottom of the Figures 5.19 to 5.23). The lateral displacement at the RAP element interface with the matrix soil increases as the applied uplift load increases. The lateral displacement profile at the maximum uplift load shows that the matrix soil is pushed outward at about the top 1 ft (0.30 m), which suggests the formation of a truncated cone near the surface during the pull-out test. All the RAP models show a similar behavior near the surface. At the bottom of the RAP element, matrix soil is pushed outward well. The lateral displacement is highest at the bottom of the RAP near the uplift plate and decreases as the depth decreases. This behavior indicates the bulging of the RAP element at the bottom during the uplift test. The bulging displacement increases as the RAP length increase. Both bulging out of the RAP element at the bottom portion of the RAP element and the formation of the truncated cone near the surface were observed by inclinometers for a 15 ft (4.57 m) RAP by Hsu (2000). The lateral displacement of the RAP element and matrix soil interface near the threaded bar (𝜃 = 45°) and between the threaded bars (𝜃 = 0°) are similar for the top and middle portion of the RAP element. However, the lateral displacement near the uplift plate is a little different at these two locations. At the bottom portion of the RAP element, the matrix soil is pushed out more near the threaded bar compare to the matrix soil between the threaded bars. The matrix soil between the threaded bars at the same depth as the uplift plate is pulled inward at the maximum uplift load as the uplift plate moved upward. The difference between the RAP and matrix soil lateral displacement close to the threaded bar and between the threaded bars decreases farther away from the uplift plate. This 3dimensional behavior usually diminishes about 2 ft (0.61 m) above the uplift plate. 106 5.2.1.4 Vertical Force Distribution in the RAP Element The vertical force distribution through the length of the RAP element at the maximum uplift load is shown in Figure 5.24. The vertical force is normalized by the maximum vertical force at the bottom of the RAP elements that the uplift plate pushed the RAP material upward. The normalized vertical force was plotted versus the normalized RAP length. The vertical force is decreased almost linearly from the bottom of the RAP element to zero at the surface for all the RAP models. 5.2.1.5 Structural Elements The threaded bar elongation and the maximum deflection of the uplift plate during the FLAC3D pull-out test are shown in Figures 5.25 to 5.29. The threaded bar elongation is calculated by calculating the difference between the displacement of the top of the threaded bar and the bottom of the threaded bar. The uplift plate maximum deflection is the difference between the vertical displacement of the node that was restrained to the threaded bar and the node at the center of the uplift plate. The threaded bar elongation and uplift plate maximum deflection were evaluated for the FLAC3D model with the three main RAP construction processes. The bar elongation increases linearly with the increase in the applied uplift force. The threaded bar elongation is very similar for different RAP construction techniques. The maximum threaded bar elongation that corresponds to the maximum uplift load ranges from 0.01 in. (0.3 mm) for the 3 ft (0.91 m) RAP to 0.25 in (6.4 mm) for the 15 ft (4.57 m) RAP. The threaded bar average axial strain ranges from 0.02% to 0.14% for the 3 ft (0.91 m) RAP and the 15 ft (4.57 m) RAP, respectively. The average axial strain of the FLAC3D 107 models is always below the yield axial strain, which is 0.26% for the Grade 75 bars. The uplift plate maximum deflection increases nonlinearly by the increase in the applied uplift force. The liner structural elements have a linear elastic material. However, the amount of the load that is transferred at different locations of the uplift plate is different and depends on the interaction of the highly nonlinear soil zones with the liner structural elements. As a result, the maximum deflection of the uplift plate versus the applied uplift load shows a nonlinear behavior and some dependency on the RAP construction process. The maximum deflection of the uplift plate ranges from 0.03 in. (0.8 mm) for the 3 ft (0.91 m) RAP to 0.42 in. (10.7 mm) for the 15 ft (4.57 m) RAP. The maximum deflection of the uplift plate is much smaller than the measured maximum deflection of 3 in. (76.2 mm) for the 15 ft (4.57 m) RAP at the end of the FST (Hsu 2000). However, the measured maximum deflection of the FST uplift plate measured the permanent plastic deformation of the uplift plate. The FLAC3D liner elements have linear material behavior, which restricts the amount of deformation by allowing the liner fiber stress goes beyond the yield stress. This limitation of the FLAC3D liner elements caused some errors in the modeling of the RAPs with a higher length. Figures 5.30 to 5.34 show the axial force of the cable elements of the threaded bar and the maximum fiber stress for the liner elements of the uplift plate FLAC3D 3 ft (0.91) to 15 ft (4.57 m) RAP models with the variable displacement RAP construction technique. The presented results are for the maximum uplift load of the FLAC3D models. The threaded bar axial force decreases from the surface to the bottom of the RAP because of the load transfer through the frictional shear stress development during the pull-out test at the perimeter of the cable elements. The frictional shear stress reduces the transferred axial 108 force from the top of the threaded bar to the uplift plate by 8% to 16%. By comparing the maximum fiber stress of the liner elements to the yield strength of the uplift plate which is 5.18 × 106 psf (248 MPa), we see that the maximum fiber stress of the uplift plate remains in the elastic range for 3 ft (0.91 m) and 6 ft (1.83 m) RAP. However, the maximum fiber stress of the uplift plate for the 9 ft (2.74 m), 12 ft (3.66 m), and 15 ft (4.57 m) RAPs exceed the steel A36 yield stress. As a result, the uplift plates showed less deflection and introduced some error in the FLAC3D models due to the linear elastic behavior of the liner elements material. 5.2.2 FLAC3D Models with Plastic-Hardening Constitutive Model The use of the PH constitutive model in the FLAC3D simulation of the RAP element for the pull-out test resulted in a less accurate prediction of the maximum uplift load compare to the FLAC3D models that used the MC constitutive model. Generally, the PH constitutive model gives a better estimation of nonlinear soil behavior. And, it has a nonlinear elasticity, plastic hardening, and strain-dependent shear modulus (shear modulus degradation curve). These features can improve the accuracy of any numerical model provided the required laboratory and in situ tests were performed in advance to parameterize the numerous PH model parameters. However, no triaxial tests were performed during the South Temple 1998 research project. And, the PH model was parameterized based on the average suggested values found in the literature for several parameters. As a result, the FLAC3D models with PH constitutive models could not provide the same level of accuracy that was achieved by the MC constitutive model. In this section, the results of the 9 ft (2.74 m) RAP FLAC3D model with the PH constitutive model 109 are discussed. The lateral displacement pattern, effective lateral stress, cavity expansion magnitude, and vector alignment (a quality measure for this technique) are plotted in Figure 5.35 for the 9 ft (2.74 m) RAP model. In the FLAC3D RAP models with the PH constitutive mode, the displacement pattern changes substantially in the 1st step, and it continues to change noticeably during the later steps as well (see Figure 5.35a). The step by step increase in the effective lateral stress shows a fair match with the target values (measured effective lateral stress) in Figure 5.35b. The cavity expansion at the first step resulted in a considerable increase in effective lateral stress. The linearization algorithm decreases the displacement pattern magnitude considerably to control the increase in the effective lateral stress in the next steps (see Figure 5.35c). The vector alignment that was defined as the cosine of the virtual angle between two vectors in Chapter 4 is plotted in Figure 5.35d. The values of 1 and -1 show a perfect alignment and 0 represents the worst alignment. The red line shows the alignment of the effective lateral stress at the end of each step with the measured values. The FLAC3D lateral effective stress shows a good alignment with the measured values (see Figure 5.35d). However, the quantitative measure of the linearization process efficiency (the blue line) shows a poor alignment. The PH model was developed based on nonlinear elasticity and plasticity. And, the FLAC3D model with PH constitutive model resulted in a highly nonlinear behavior for the RAP system compared to the FLAC3D model with MC constitutive model at each step. And the increase in the nonlinear behavior of the soil at each step decreases the efficiency of the linearization process drastically. Therefore, the linearization process with relatively few steps is not suggested for highly nonlinear constitutive models such as the PH model. 110 Figure 5.36 shows the load-displacement curves for the 9 ft (2.74 m) FST and FLAC3D models with PH constitutive model. The FLAC3D curves were plotted for the models with lift compaction, cavity expansion, and variable displacement RAP construction techniques. The FLAC3D models resulted in similar load-displacement curves. The FST’s maximum uplift load is underestimated by 35%, 28%, and 31% with the lift compaction, cavity, expansion, and variable displacement techniques, respectively. The displacements of FLAC3D models with the main three RAP construction techniques are within 0.05 in. (1 mm) of each other. However, the displacement is always less than the FST’s measured values. By comparing the load-displacement curves of the PH model (see Figure 5.36) with the MC model load-displacement curves (see Figure 5.14), we can see that the PH models resulted in a much stiffer RAP system behavior with a lower maximum uplift load. The accuracy of the PH model prediction of the maximum uplift load is about 11% to 23% lower than the MC model accuracy. The average decrease in the lateral stress within the length of each RAP at the maximum uplift load ranges from 252 psf (12 kPa) to 518 psf (24.8 kPa) for the 9 ft (2.74 m) RAP model with the main three RAP construction techniques (see Figure 5.37). The PH model and MC model for FLAC3D models with the cavity expansion and variable displacement RAP construction techniques resulted in a similar average decrease in the horizontal stress during the pull-out simulation (see Figures 5.21 and 5.37). However, the PH model predicted values for the average decrease in lateral stress is about 37% lower than the MC model average decrease of the lateral stress for the lift compaction techniques. The lateral displacement profile of the RAP and matrix soil interface for the PH model (see Figure 5.37) follows similar behavior that was predicted by the MC model (see 111 Figure 5.21). The lateral displacement profile at the maximum uplift load shows that the matrix soil is pushed outward at about the top 1 ft (0.30 m) and the bottom 1 ft (0.30 m) to 2 ft (0.61 m) of the RAP elements. Therefore, the PH model was able to predict the formation of a truncated cone near the surface and bulging out the bottom portion of the RAP during the pull-out test. However, the predicted values for the lateral displacement by the PH model are lower than the MC model predicted lateral displacements by a factor of about 2. The lateral displacement of the RAP element and matrix soil interface near the threaded bar (𝜃 = 45°) and between the threaded bars (𝜃 = 0°) are similar for the top and middle portion of the RAP element. At the bottom portion of the RAP element, the matrix soil is pushed out more near the threaded bar compare to the matrix soil between the threaded bars. The matrix soil between the threaded bars at the same depth as the uplift plate is pulled inward at the maximum uplift load as the uplift plate moved upward. The threaded bar elongation and the maximum deflection of the uplift plate during the FLAC3D pull-out test are shown in Figure 5.38. The bar elongation increases linearly with the increase in the applied uplift force. The threaded bar elongation is very similar for different RAP construction techniques. The maximum threaded bar elongation that corresponds to the maximum uplift load is 0.05 in. (1.3 mm) for the 9 ft (2.74 m) RAP (an average axial strain of 0.04%). The PH model predicted 33% less elongation for the threaded bar compare to the MC model estimation of the threaded bar elongation (see Figure 5.27). The uplift plate maximum deflection increases nonlinearly by the increase in the applied uplift force. The maximum deflection of the uplift plate ranges from 0.05 in. (1.3 mm) for the 3 to 0.08 in. (2.0 mm) for different RAP construction process in FLAC3D. 112 The maximum deflection of the uplift plate by the PH model is about half the maximum deflection of the uplift plate by the MC model (see Figure 5.27). Figure 5.39 shows the axial force of the cable elements of the threaded bar and the maximum fiber stress for the liner elements of the uplift plate FLAC3D 9 ft (2.74 m) RAP models with the variable displacement RAP construction technique. Similar to the MC model results, the threaded bar axial force decreases from the surface to the bottom of the RAP. The frictional shear stress reduces the transferred axial force from the top of the threaded bar to the uplift plate by 16%. This value is similar to the portion of the uplift load that was transferred directly to the RAP material through the frictional shear development in the MC model at the maximum uplift load (see Figure 5.32). The maximum fiber stress of the liner elements is slightly less than the yield stress of the uplift plate material. Unlike the uplift plate in the FLAC3D RAP model with the MC model, the uplift plate stress in the PH model RAP model remained in the elastic range. 113 Table 5.1. The uplift capacity compression for different RAP construction techniques 114 0.10 Displacement of uplift bars, δv (in.) FLAC3D Anchor (H = 3 ft, E = 400 ksf, fi = 35) FLAC3D Anchor (H = 3 ft, E = 800 ksf, fi = 35) 0.08 Meyerhof and Adams Vesic 0.06 0.04 0.02 0.00 0 1 2 3 4 Uplift load, T (kips) a. 5 6 7 Displacement of uplift bars, δv (in.) 1.0 FLAC3D Anchor (H = 9 ft, E = 400 ksf, fi = 30) Meyerhof and Adams Vesic 0.8 0.6 0.4 0.2 0.0 0 10 20 30 Uplift load, T (kips) b. 40 50 Figure 5.1. Load-displacement curves for the anchor verification problem. a. 3 ft anchor, and b. 9 ft anchor 115 a. b. Figure 5.2. Coefficient of lateral earth pressure after FLAC3D RAP construction process of 3 ft RAP. a. Lift compaction, and b. Cavity Expansion 116 a. b. Figure 5.3. Coefficient of lateral earth pressure after FLAC3D RAP construction process of 6 ft RAP. a. Lift compaction, and b. Cavity Expansion 117 a. b. Figure 5.4. Coefficient of lateral earth pressure after FLAC3D RAP construction process of 9 ft RAP. a. Lift compaction, and b. Cavity Expansion 118 a. b. Figure 5.5. Coefficient of lateral earth pressure after FLAC3D RAP construction process of 12 ft RAP. a. Lift compaction, and b. Cavity Expansion 119 a. b. Figure 5.6. Coefficient of lateral earth pressure after FLAC3D RAP construction process of 15 ft RAP. a. Lift compaction, and b. Cavity Expansion 120 Figure 5.7. Increase in effective lateral stress in few steps to build the 3 ft RAP by Cavity expansion with variable displacement technique. a. Displacement pattern, b. Effective lateral stress, c. Displacement magnitude, and d. Vector alignment 121 Figure 5.8. Increase in effective lateral stress in a few steps to build the 6 ft RAP by Cavity expansion with variable displacement technique. a. Displacement pattern, b. Effective lateral stress, c. Displacement magnitude, and d. Vector alignment 122 Figure 5.9. Increase in effective lateral stress in a few steps to build the 9 ft RAP by Cavity expansion with variable displacement technique. a. Displacement pattern, b. Effective lateral stress, c. Displacement magnitude, and d. Vector alignment 123 Figure 5.10. Increase in effective lateral stress in a few steps to build the 12 ft RAP by Cavity expansion with variable displacement technique. a. Displacement pattern, b. Effective lateral stress, c. Displacement magnitude, and d. Vector alignment 124 Figure 5.11. Increase in effective lateral stress in a few steps to build the 15 ft RAP by Cavity expansion with variable displacement technique. a. Displacement pattern, b. Effective lateral stress, c. Displacement magnitude, and d. Vector alignment 125 Figure 5.12. Load-displacement curves of 3 ft RAP uplift test 126 Figure 5.13. Load-displacement curves for 6ft RAP uplift test 127 Figure 5.14. Load-displacement curves for 9 ft RAP uplift test 128 Figure 5.15. Load-displacement curves for 12 ft RAP uplift test 129 Figure 5.16. Load-displacement curves for 15 ft RAP uplift test 130 a. b. Figure 5.17. Maximum shear strain rate during pull-out with the applied uplift load higher than the uplift capacity. a. 3 ft RAP, and b. 9 ft RAP 131 2.0 FST 1.8 No RAP construction K0-Blade stress 1.6 Lift compaction TFLAC3D/TFST 1.4 Cavity displacement Variable displacement 1.2 1.0 0.8 0.6 0.4 0.2 0.0 0 3 6 9 RAP length, Hg (ft) 12 15 Figure 5.18. The uplift capacity compression for different RAP construction techniques 132 Figure 5.19. Change in lateral stress and lateral displacement at the matrix soil close to the 3 ft RAP during the FLAC3D uplift test for RAP construction. a., b., and c. change in lateral stress, and d., e., and f. change in interface lateral displacement between the threaded rods (𝜃 = 0°) in front of the threaded rod (𝜃 = 45°) 133 Figure 5.20. Change in lateral stress and lateral displacement at the matrix soil close to the 6 ft RAP during the FLAC3D uplift test for RAP construction. a., b., and c. change in lateral stress, and d., e., and f. change in interface lateral displacement between the threaded rods (𝜃 = 0°) in front of the threaded rod (𝜃 = 45°) 134 Figure 5.21. Change in lateral stress and lateral displacement at the matrix soil close to the 9 ft RAP during the FLAC3D uplift test for RAP construction. a., b., and c. change in lateral stress, and d., e., and f. change in interface lateral displacement between the threaded rods (𝜃 = 0°) in front of the threaded rod (𝜃 = 45°) 135 Figure 5.22. Change in lateral stress and lateral displacement at the matrix soil close to the 12 ft RAP during the FLAC3D uplift test for RAP construction. a., b., and c. change in lateral stress, and d., e., and f. change in interface lateral displacement between the threaded rods (𝜃 = 0°) in front of the threaded rod (𝜃 = 45°) 136 Figure 5.23. Change in lateral stress and lateral displacement at the matrix soil close to the 15 ft RAP during the FLAC3D uplift test for RAP construction. a., b., and c. change in lateral stress, and d., e., and f. change in interface lateral displacement between the threaded rods (𝜃 = 0°) in front of the threaded rod (𝜃 = 45°) 137 Figure 5.24. The FLAC3D models' vertical force distribution through the RAP length at the maximum uplift load. a. lift compaction, b. cavity expansion, and c. variable displacement 138 Figure 5.25. Threaded bar elongation and maximum deformation of the uplift plate (differential vertical displacement of the node at the center of the plate and the node connected to the threaded bar) during the FLAC3D uplift test of 3 ft RAP. a. Lift compaction, b. Cavity expansion, and c. Cavity expansion with variable displacement 139 Figure 5.26. Threaded bar elongation and maximum deformation of the uplift plate (differential vertical displacement of the node at the center of the plate and the node connected to the threaded bar) during the FLAC3D uplift test of 6 ft RAP. a. Lift compaction, b. Cavity expansion, and c. Cavity expansion with variable displacement 140 Figure 5.27. Threaded bar elongation and maximum deformation of the uplift plate (differential vertical displacement of the node at the center of the plate and the node connected to the threaded bar) during the FLAC3D uplift test of 9 ft RAP. a. Lift compaction, b. Cavity expansion, and c. Cavity expansion with variable displacement 141 Figure 5.28. Threaded bar elongation and maximum deformation of the uplift plate (differential vertical displacement of the node at the center of the plate and the node connected to the threaded bar) during the FLAC3D uplift test of 12 ft RAP. a. Lift compaction, b. Cavity expansion, and c. Cavity expansion with variable displacement 142 Figure 5.29. Threaded bar elongation and maximum deformation of the uplift plate (differential vertical displacement of the node at the center of the plate and the node connected to the threaded bar) during the FLAC3D uplift test of 15 ft RAP. a. Lift compaction, b. Cavity expansion, and c. Cavity expansion with variable displacement 143 a. b. Figure 5.30. Structural elements force and stress at the maximum uplift load for FLAC3D 3 ft RAP model with variable displacement RAP construction. a. cable force, and b. liner maximum fiber stress 144 a. b. Figure 5.31. Structural elements force and stress at the maximum uplift load for FLAC3D 6 ft RAP model with variable displacement RAP construction. a. cable force, and b. liner maximum fiber stress 145 a. b. Figure 5.32. Structural elements force and stress at the maximum uplift load for FLAC3D 9 ft RAP model with variable displacement RAP construction. a. cable force, and b. liner maximum fiber stress 146 a. b. Figure 5.33. Structural elements force and stress at the maximum uplift load for FLAC3D 12 ft RAP model with variable displacement RAP construction. a. cable force, and b. liner maximum fiber stress 147 a. b. Figure 5.34. Structural elements force and stress at the maximum uplift load for FLAC3D 15 ft RAP model with variable displacement RAP construction. a. cable force, and b. liner maximum fiber stress 148 Figure 5.35. Increase in effective lateral stress in a few steps to build the 9 ft RAP by Cavity expansion with variable displacement technique (PH Model). a. Displacement pattern, b. Effective lateral stress, c. Displacement magnitude, and d. Vector alignment 149 Figure 5.36. Load-displacement curves for 9 ft RAP uplift test (PH Model) 150 Figure 5.37. Change in lateral stress and lateral displacement at the matrix soil close to the 9 ft RAP during the FLAC3D uplift test for RAP construction (PH Model). a., b., and c. change in lateral stress, and d., e., and f. change in interface lateral displacement between the threaded rods (𝜃 = 0°) in front of the threaded rod (𝜃 = 45°) 151 Figure 5.38. Threaded bar elongation and maximum deformation of the uplift plate (differential vertical displacement of the node at the center of the plate and the node connected to threaded bar) during the FLAC3D uplift test of 9 ft RAP (PH Model). a. Lift compaction, b. Cavity expansion, and c. Cavity expansion with variable displacement 152 a. b. Figure 5.39. Structural elements force and stress at the maximum uplift load for FLAC3D 9 ft RAP model with variable displacement RAP construction (PH Model). a. cable force, and b. liner maximum fiber stress CHAPTER 6 CONCLUSION 6.1 Concluding Remarks Numerical analyses were conducted using the program Fast Lagrangian Analysis of Continua 3D (FLAC3D) to simulate full-scale tests (FSTs) in uplift that were performed on five Geopier Rammed Aggregate Piers (RAPs). All five RAPs tested had a nominal diameter of 3 ft (0.91 m) but varied in height from 3 to 15 ft (0.91 to 4.57 m) in increments of 3 ft (0.91 m). Two constitutive models were used in the numerical analyses: (1) An elastic-perfectly plastic Mohr-Coulomb model, and (2) a Plastic-Hardening model that employs nonlinear stress-strain behavior in the elastic and plastic ranges. Three methods were used to simulate the process of constructing the RAPs. Two of these methods had been used by previous researchers (lift compaction and cavity expansion), while the third method (cavity expansion with variable displacement) was developed as part of this research. The primary conclusions determined by comparing the results of the FLAC3D analyses with those from the FSTs are as follows: 1- The FLAC3D analyses with MC constitutive model predicted the FST’s uplift capacity of the single RAP during the uplift test with a 1% to 40% error depending on the RAP construction technique and the RAP length. The average accuracy of the lift compaction, cavity expansion, and cavity expansion with 154 variable displacement RAP construction techniques is almost the same for each RAP length. The FLAC3D model for the 6 ft (1.83 m) RAP had the best prediction for the maximum uplift load with 1% to 17% error and FLAC3D model for the 3 ft (0.91 m) RAP has the least accuracy with 36% to 40% error for different RAP construction techniques. The low overburden pressure at the shallow depth decreases the efficiency of the RAP construction techniques. 2- Results from the FLAC3D analyses showed that the computationally expensive RAP construction process could be replaced with a horizontal stress initialization without any considerable loss of accuracy. 3- The proposed linearization process successfully adjusted the lateral displacement at the borehole shaft to induce effective lateral stresses similar to the measured values in the field. 4- The majority of the FLAC3D simulations showed a decrease in the effective lateral stress during the uplift test. However, the lateral stress decrease in the FLAC3D models was higher than the measured values during the FST test for short piers. This could be the reason for the underestimation of the maximum uplift load by the FLAC3D model. 5- The FLAC3D models predicted the formation of a truncated cone near the surface and bulging out of the RAP element near the uplift plate, which had been observed in the field tests. 6- The FLAC3D model with linear elastic liner elements for the uplift plate could not predict the measured plastic deformation in the field for longer RAPs. The 3D displacement of the soil with the elastic uplift plate is limited to about the 155 bottom 2 ft (1.83 m) of the RAP element. However, the large plastic deformation of the uplift plate, which was measured in the field, could increase the 3D zone. 7- The FLAC3D model with the Plastic-Hardening constitutive model could not predict the load-displacement of the FSTs with the same level of accuracy as the FLAC3D model with Mohr-Coulomb constitutive model. The lower accuracy of the more sophisticated Plastic-Hardening model is due to the estimation of the numerous parameters of this constitutive model with average material properties reported in the literature. 6.2 Recommendations The following recommendations are provided regarding potential improvements to the procedures used to model the performance of RAPS numerically: 1- The modeling of the RAP construction procedure can be improved by using numerical modeling methods that are able to model the interparticle interaction (such as Discrete Element Method) for the RAP element and the matrix soil fracture. 2- Advanced constitutive models can improve numerical modeling. However, they usually require more advance laboratory and in situ testing in a greater quantity. Also, planning the type of the numerical model to be used in the numerical analyses prior to conducting the field investigation could improve the constitutive model parameterization. 3- Using a numerical model that is able to model the large plastic deformation of 156 the uplift plate for the longer RAPs could improve the prediction of the 3D soil displacement. A simplified 2D axisymmetric equivalent model could be recommended by having a better understanding of the realistic 3D distribution of the displacement and lateral stress. 4- The cavity expansion at each lift could be estimated in the field by measuring the volume of the RAP material, that is placed in the borehole, and the relative density before and after the lift compaction. Then the estimated lateral displacement could be used in the numerical model to expand the cavity and simulate the RAP construction process. APPENDIX A LIST OF NOTATIONS 158 𝛼 = A constant derived internally in FLAC3D based on a virtual oedometer test 𝑛 𝛼𝑠𝑐 = The scaling factor for the displacement pattern of step 𝑛 𝑛−1 𝛼𝑠𝑐 = The scaling factor at the previous step (step 𝑛 − 1) 𝛾 = Total unit weight of soil 𝛾 𝑝 = Internal shear hardening parameter in FLAC3D PH model 𝛾 𝑡 = Shear strain 𝑡 𝛾70 = The reference shear strain corresponding to the secant shear modulus with the value of 72.2% of the initial shear modulus 𝛾𝑅𝐴𝑃 = Total unit weight of RAP material ′ 𝛾𝑅𝐴𝑃 = Effective unit weight of RAP material 𝛾𝑠𝑎𝑡 = Unit weight of saturated soil 𝛾𝑠𝑎𝑡,𝑅𝐴𝑃 = Saturated unit weight of RAP material 𝛿⃗ = A vector consisting of shaft’s lateral displacement 𝛿⃗𝑛 = Lateral displacement vector at step 𝑛 with the size of (𝑛𝑙𝑖𝑓𝑡) 𝛿ℎ = Lateral displacement 𝛿𝑣 = Vertical displacement Δ𝛾 𝑝 = Internal shear hardening parameter increment in FLAC3D PH model Δ𝛿⃗𝑛 = Lateral displacement increment vector at step 𝑛, (𝛿⃗𝑛+1 − 𝛿⃗𝑛 ) th Δ𝛿⃗𝑖𝑛 = Lateral displacement increment vector with only the i component having a nonzero value, (𝛿⃗𝑖𝑛 − 𝛿⃗𝑛 ) 159 𝑛 ‖Δ𝛿⃗𝑡𝑎𝑟𝑔𝑒𝑡 ‖= 𝑛 Magnitude (l2-norm) of the Δ𝛿⃗𝑡𝑎𝑟𝑔𝑒𝑡 vector Δ𝑃⃗⃗𝑛 = Effective radial stress increment vector at step 𝑛, (𝑃⃗⃗ 𝑛+1 − 𝑃⃗⃗𝑛 ) Δ𝑃⃗⃗𝑖𝑛 = Effective lateral stress increment vector after expanding the cavity at the ith lift, (𝑃⃗⃗𝑖𝑛 − 𝑃⃗⃗𝑛 ) Δ𝑦 = The independent variable increment ∆𝑧𝑚𝑖𝑛 = The smallest dimension of an adjoining zone in the normal direction 𝜀1 = Axial strain 𝜁 = An internal parameter in FLAC3D PH model 𝜃𝑥,𝑦 = The virtual angle between two vectors 𝑥⃗ and 𝑦⃗ 𝜈 = Poisson’s ratio 𝜌 = Soil density 𝜎1 = Major principal stress 𝜎2 = Median principal stress 𝜎3 = Minor principal stress 𝜎1′ = Major principal effective stress 𝜎2′ = Median principal effective stress 𝜎3′ = Minor principal effective stress ′ 𝜎𝑚 = Effective confining stress that acts at the plane perpendicular to the FLAC3D cable element axis 𝜎𝑛 = Interface normal stress of the FLAC3D liner element 𝜎𝑝′ = Over consolidation pressure 𝜎𝑣0 = Initial total vertical stress 160 ′ 𝜎𝑣0 = Initial effective vertical stress ∑ 𝑃⃗⃗0 = The sum of the initial effective lateral stress at step 𝑛 ∑ 𝑃⃗⃗𝑛 = The sum of the effective lateral stress at step 𝑛 ∑ 𝑃⃗⃗𝑡𝑎𝑟𝑔𝑒𝑡 = The sum of the measured effective lateral stress at the end of RAP construction ∅ = Friction angle ∅′ = Effective friction angle ∅𝑔 = Grout friction angle in FLAC3D cable element ∅𝑙𝑖𝑛𝑒𝑟 = Shear coupling spring friction angle of FLAC3D liner element ∅𝑅𝐴𝑃 = RAP material friction angle ∅′𝑅𝐴𝑃 = RAP material effective friction angle 𝜓 = Dilatation angle 𝜓𝑅𝐴𝑃 RAP material dilatation angle 𝐷𝑖𝑛 = Lateral displacement magnitude of the cavity expansion at the ith lift in step 𝑛 during single lift expansion. 𝑐 = Cohesion 𝑐𝑙𝑖𝑛𝑒𝑟 = Shear coupling spring cohesion of FLAC3D liner element 𝑐𝑔 = Grout cohesive strength of FLAC3D cable element ′ 𝑐𝑅𝐴𝑃 = RAP material effective cohesion 𝐷𝑐𝑎𝑏𝑙𝑒 = Reinforcement diameter of FLAC3D cable element 𝐸50 = Secant modulus at half the deviatoric stress from the Mohr-Coulomb failure criterion 161 𝑟𝑒𝑓 𝐸50 = Reference secant modulus at half the deviatoric stress from the Mohr- Coulomb failure criterion corresponding to the reference stress (𝑝𝑟𝑒𝑓 ) 𝐸𝑖 = Initial tangent modulus 𝐸0 = Initial Young’s modulus (at very small strain) 𝑟𝑒𝑓 𝐸0 = Initial Young’s modulus (at very small strain) at the reference pressure 𝑟𝑒𝑓 𝐸𝑜𝑒𝑑 = tangent oedometer stiffness at the reference pressure 𝐸𝑠,𝑁𝐶 = Young’s modulus for normally consolidated soil 𝐸𝑠,𝑂𝐶 = Young’s modulus for overconsolidated soil 𝐸𝑠𝑡 = Uplift plate and threaded bar Young’s modulus 𝐸𝑢𝑟 = Unloading Young’s modulus 𝑟𝑒𝑓 𝑟𝑒𝑓 𝐸𝑢𝑟 = Unloading Young’s modulus at the reference stress (𝑝 ) 𝑓⃗ = A multivariable function 𝑓 ′ (𝑥)|𝑥0 = The first derivative of 𝑦 with respect to 𝑥 evaluated at 𝑥0 [∇𝑓⃗|𝑥⃗ ] = The Jacobian matrix evaluated at 𝑥⃗0 with the size of (𝑡 × 𝑟) 0 𝑓𝑀𝐶 = PH model shear yield surface 𝑓𝑠 = CPT sleeve friction 𝑓𝑡 = Normal coupling spring tensile strength of FLAC3D liner element 𝑔 = Earth gravitational acceleration 𝑔𝑠𝑐 (𝑁𝐿𝐹 𝑛 ) = A function that determines the desired increase in 𝑁𝐿𝐹 for the next step 𝐺 = Shear modulus 𝐺0 = Initial shear modulus (at very small strain) 162 𝐺𝑔 = Grout shear modulus of FLAC3D cable element 𝐺𝑅𝐴𝑃 = RAP material bulk modulus 𝐺𝑠𝑒𝑐 = Secant shear modulus 𝐺𝑡𝑎𝑛 = Tangent shear models [𝐽𝑛 ] = Jacobian matrix at step 𝑛 with the size of (𝑛𝑙𝑖𝑓𝑡 × 𝑛𝑙𝑖𝑓𝑡) [𝐽𝑛 ]−1 = The inverse of the Jacobian matrix 𝑘 = Correction factor with a typical value of 0 < 𝑘 ≤ 1 𝑘𝑔 = Grout shear stiffness per unit length of FLAC3D cable element 𝐾 = Bulk modulus 𝐾0 = Coefficient of lateral earth pressure at rest 𝐾0,𝑁𝐶 = Coefficient of lateral earth pressure at rest for normally consolidated soil with the maximum value of 𝜈 ⁄(1 − 𝜈) 𝐾𝑝 = Rankine’s coefficient of passive earth pressure 𝐾𝑅𝐴𝑃 = RAP material bulk modulus 𝑚 = Exponent that determine the amount of stress dependency of 𝐸50 𝑝′ = Effective mean pressure 𝑝𝑐 = Volumetric hardening parameter 𝑝𝑟𝑒𝑓 = Reference stress 𝑃⃗⃗ = A vector consisting of the effective radial stress in the matrix soil at a close distance from the borehole shaft 𝑃⃗⃗𝑛 = Effective radial stress vector at step 𝑛 with the size of (𝑛𝑙𝑖𝑓𝑡) 𝑃⃗⃗𝑡𝑎𝑟𝑔𝑒𝑡 = Effective radial stress vector at 1 ft (0.30 m) from the borehole shaft 163 (2.5 ft (0.76 m) from the center of the RAP element) measured with K0-Blade test after the RAP construction 𝑞 = Deviatoric stress, 𝜎3 − 𝜎1 𝑞̃ = A special shear stress measure 𝑞𝑎 = Deviatoric stress at failure, (𝜎3 − 𝜎1 )𝑓 𝑞𝑐 = CPT tip resistance 𝑞𝑡 = CPT corrected, total tip resistance 𝑞𝑓 = Deviatoric stress from Mohr-Coulomb failure criterion 𝑝𝑟𝑒𝑓 = Reference stress 𝑅𝑓 = Failure ratio 𝑢0 = In situ equilibrium pore water pressure 𝑢2 = Pore water pressure measured during CPT sounding 𝑉𝑠 = Shear wave velocity 𝑥 = A scalar independent variable 𝑥⃗ = A vector of independent variables with the size of (𝑟) 𝑥0 = The point of interest of the independent variable 𝑥⃗0 = The independent variables vector at points of interest 𝑦 = A scalar dependent variable 𝑦⃗ = A vector of dependent variables with the size of (𝑡) 𝑦0 = The dependent variable value at 𝑥0 APPENDIX B FLAC3D CODE 165 FLAC3D code Used for Uplift Simulation of a single RAP ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ========================================================================== ========================================================================== Simulation of pull-test for uplift RAPs South Temple 1998 research project Emad Ghodrati University of Utah Thesis 2019 Units: Length in feet Force in lb Stress in psf Subroutine name: "Main.dat" Programing language: FLAC3D commands and FISH Description: This is the main subroutine that calls other subroutines ========================================================================== ========================================================================== ; Building the model geometry and initializing horizontal stress call 'Model_Initialize' call 'K0_Stress' ;--------------------------------------------------------------------------------; Pullout simulation with Mohr-Coulomb Model ;--------------------------------------------------------------------------------; Assigning Constitutive model model restore 'Ini_Stress.sav' [global Constit_model = 'MC'] model save 'Ini_Stress' call 'Constitutive_Model' ;---------------------------------------; Borehole Excavation call 'RAP_Excavate' ;---------------------------------------; Anchor pullout test call 'Anchor_Const' call 'Plate_Rod' call 'LoadTest' ;---------------------------------------; Pullout test of RAP build with K_Blade stress call 'RAP_Const_K_Blade' call 'Plate_Rod' call 'LoadTest' ;---------------------------------------; Pullout test of RAP build with Cavity Expansion call 'RAP_Const_Cavity_velocity' call 'Plate_Rod' call 'LoadTest' ;---------------------------------------; Pullout test of RAP build with Lift Compaction call 'RAP_Const_Lift_Pressure' call 'Plate_Rod' call 'LoadTest' ;---------------------------------------; Pullout test of RAP build with Variable Velocity Cavity Expansion call 'RAP_Const_Cavity_variable_velocity.py' call 'Plate_Rod' call 'LoadTest' 166 ;--------------------------------------------------------------------------------; Pullout simulation with Plastic-Hardening Model ;--------------------------------------------------------------------------------; Initializing horizontal stress model restore 'Ini_Stress.sav' [global Constit_model = 'PH'] model save 'Ini_Stress' call 'Constitutive_Model' ;---------------------------------------; Borehole Excavation call 'RAP_Excavate' ;---------------------------------------; Anchor pullout test call 'Anchor_Const' call 'Plate_Rod' call 'LoadTest' ;---------------------------------------; Pullout test of RAP build with K_Blade stress call 'RAP_Const_K_Blade' call 'Plate_Rod' call 'LoadTest' ;---------------------------------------; Pullout test of RAP build with Cavity Expansion call 'RAP_Const_Cavity_velocity' call 'Plate_Rod' call 'LoadTest' ;;---------------------------------------; Pullout test of RAP build with Lift Compaction call 'RAP_Const_Lift_Pressure' call 'Plate_Rod' call 'LoadTest' ;---------------------------------------; Pullout test of RAP build with Variable Velocity Cavity Expansion call 'RAP_Const_Cavity_variable_velocity.py' call 'Plate_Rod' call 'LoadTest' ; ; ; ; ; ; ; ; ; ; ; ========================================================================== ========================================================================== Subroutine name: "Model_Initialize.dat" Programing language: FLAC3D commands and FISH Description: This subroutine set the geometry parameters, calls the “geometry.dat” to build the model, initialize the model properties, set the boundary conditions, creates a double grid points layer at the uplift plate location, and solve the initial model to reach equilibrium ========================================================================== ========================================================================== model new fish automatic-create off model title "Pullout test simulation of 9 ft RAPs" model largestrain off ;--------------------------------------------------------------------------------;Model geometry parameters ;--------------------------------------------------------------------------------;The length unit is ft and the force unit is lb fish define GenSetting ;Uplift RAP geometry parameters global rPl = 18./12 ;Plate radius global zPl = -9.0 ;Plate depth 167 global ThPl = 0.5/12 global dRod = 1./12 global sin45 = math.sin(math.pi/4) global ClearRod = 1.5625/12 global xyRod = (rPl-ClearRod)*sin45 ;Size of the model global rMod = 10. global zbot = zPl-4*(2*rPl) ;Water table global grav = 32.174 global wdens = 62.4/grav global zwt = -4. ;Initial matrix soil properties for global fiSoil = 37.8 global cSoil = 0.0 global KSoil = 1.363e6 global GSoil = 1.022e6 global RoSoil = 109./grav ;RAPs properties global LiftTh = 1. global fiRAP = 47. global cRAP = 1.0 global saiRAP = 12. global KRAP = 6.670e6 global GRAP = 1.430e6 global RoRAP = 135./grav ;Plate thickness ;Rod diameter ;Rod clearance ;x and y coordinate of the the rod ;Radius of the model boundary ;z coordinate of the model boundary ;Gravitational acceleration ;Water density ;GWT depth model initialization ;Friction angle of the soil ;Cohesion of the soil ;Bulk modulus of the soil ;Shear modulus of the soil ;Mass density of the soil ;RAP lift thickness ;Friction angle of the RAP ;Cohesion of the RAP ;Dialation angle of RAP ;Bulk modulus of the RAP ;Shear modulus of the RAP ;Mass density of the RAP end @GenSetting ;--------------------------------------------------------------------------------;Calling the geometry subroutine ;--------------------------------------------------------------------------------call 'geometry' ;--------------------------------------------------------------------------------;Assigning group name ;--------------------------------------------------------------------------------zone group 'MSoil' slot 'Soil' zone group 'RAP' slot 'Soil' range cylinder ... end-1 (0,0,0) end-2 (0,0,@zPl) radius @rPl ;---------------------------------------zone face group 'FC_Xside' slot 'BC' range position-y -0.01 0.01 zone face group 'FC_YSide' slot 'BC' range position-x -0.01 0.01 zone face group 'FC_Bot' slot 'BC' range position-z [zbot-0.01] [zbot+0.01] zone face group 'FC_ModelSide' slot 'BC' range cylinder ... end-1 (0,0,0) end-2 (0,0,@zbot) radius [rMod-0.01] [rMod+0.01] zone face group 'FC_Plate' slot 'Int' internal range cylinder ... end-1 (0,0,[zPl-0.01]) end-2 (0,0,[zPl+0.01]) radius @rPl zone face group 'FC_SideRAP' slot 'Int' internal range cylinder ... end-1 (0,0,0) end-2 (0,0,@zPl) radius [rPl-0.01] [rPl+0.01] ;--------------------------------------------------------------------------------;Create zone separation and re-attach for initialization ;--------------------------------------------------------------------------------zone separate by-face new-side group 'FC_PlateTop' slot 'Int' ... origin (0,0,@zbot) range group 'FC_Plate' slot 'Int' zone attach by-face range cylinder end-1 (0,0,[zPl-0.01]) ... end-2 (0,0,[zPl+0.01]) radius @rPl ;--------------------------------------------------------------------------------; General BCs of the domain ;--------------------------------------------------------------------------------zone face apply velocity-normal 0 range group 'FC_Xside' slot 'BC' zone face apply velocity-normal 0 range group 'FC_YSide' slot 'BC' zone face apply velocity-normal 0 range group 'FC_ModelSide' slot 'BC' 168 zone face apply velocity (0,0,0) range group 'FC_Bot' slot 'BC' ;--------------------------------------------------------------------------------; Soil constitutive model ;--------------------------------------------------------------------------------zone cmodel assign elastic zone property bulk @KSoil shear @GSoil ;--------------------------------------------------------------------------------; Initialize gravity, density, pore water pressure ,and stress state ;--------------------------------------------------------------------------------model gravity @grav zone water density @wdens zone water plane normal (0,0,-1) origin (0,0,@zwt) zone initialize density @RoSoil zone initialize-stress ratio 0.5 ;--------------------------------------------------------------------------------;Getting convergence parameter history and solving the initial model ;--------------------------------------------------------------------------------model update-interval 250 zone ratio local model history mechanical ratio model history mechanical unbalanced-maximum model solve model save 'initial' ; ; ; ; ; ; ; ; ========================================================================== ========================================================================== Subroutine name: "geometry.dat" Programing language: FLAC3D commands and FISH Description: This subroutine builds the model geometry and discretize the model with an unstructured mesh ========================================================================== ========================================================================== extrude point create (0.0,0.0) ;Point 1 extrude point create (0.0,[rPl-2*ClearRod]) ;Point 2 extrude point create ([rPl-2*ClearRod],0.0) ;Point 3 extrude point create (0.0,@rPl) ;Point 4 extrude point create (@rPl,0.0) ;Point 5 extrude point create (0.0,@rMod) ;Point 6 extrude point create (@rMod,0.0) ;Point 7 ; Creating edges extrude edge create by-points 1 2 type polyline ;Edge 1 extrude edge create by-points 1 3 type polyline ;Edge 2 extrude edge create by-points 2 4 type polyline ;Edge 3 extrude edge create by-points 3 5 type polyline ;Edge 4 extrude edge create by-points 4 6 type polyline ;Edge 5 extrude edge create by-points 5 7 type polyline ;Edge 6 extrude edge create by-points 2 3 type arc ;Edge 7 extrude edge create by-points 4 5 type arc ;Edge 8 extrude edge create by-points 6 7 type arc ;Edge 9 ; Adjusting the curve extrude edge id 9 control-point add ([rMod*sin45],[rMod*sin45]) extrude edge id 8 control-point add ([rPl*sin45],[rPl*sin45]) extrude edge id 7 control-point add (@xyRod,@xyRod) ; Mesh parameter extrude mesh type unstructured extrude mesh mode quad extrude mesh target-size 0.1 extrude mesh gradation 0.01 extrude mesh optimization 10 ; Meshing 169 extrude edge size 6 range id-list 1 2 extrude edge size 1 range id-list 3 4 extrude edge size 10 range id-list 7 8 extrude edge size 22 range id-list 9 extrude edge size 18 range id-list 5 6 extrude edge id 5 ratio 1.05 extrude edge id 6 ratio 1.05 extrude set system normal-axis (0,0,-1) extrude block create automatic extrude segment index 1 position [-1*zbot] extrude segment add position [-1*zPl] extrude segment index 1 size [math.round(-1*zPl*12/3)] extrude segment index 2 size 24 extrude segment index 2 ratio 1.04 zone generate from-extruder model save 'geometry' ; ; ; ; ; ; ; ; ; ========================================================================== ========================================================================== Subroutine name: " K0_Stress.dat" Programing language: FLAC3D commands and FISH Description: This subroutine imports the measured coefficient of lateral earth pressure at rest and initializes the effective lateral stress in the model ========================================================================== ========================================================================== model restore 'initial' ;--------------------------------------------------------------------------------zone gridpoint initialize displacement (0,0,0) zone gridpoint initialize velocity (0,0,0) ;--------------------------------------------------------------------------------; K0 table from the K0-blade test complemented with CPT data for deeper depth ;--------------------------------------------------------------------------------table 'TB_K0' delete table 'TB_K0' import 'Input_K0' ;--------------------------------------------------------------------------------; Import the measured K for comparison ;--------------------------------------------------------------------------------table 'TB_K25_target' delete table 'TB_K50_target' delete table 'TB_K25_target' import 'Input_K25-target' table 'TB_K50_target' import 'Input_K50-target' ;--------------------------------------------------------------------------------; Update horizontal stress ;--------------------------------------------------------------------------------fish define HorStress(gpname,sl) local zp = zone.head local zdep local effszz local sxx local syy local K0 ;-----------------------------------loop foreach zp zone.list if zone.model(zp) # 'null' then if zone.group(zp,sl) = gpname then zdep = -1*zone.pos.z(zp) K0 = table('TB_K0',zdep) effszz = zone.stress.zz(zp)+zone.pp(zp) ; Effective vertical stress sxx = K0*effszz-zone.pp(zp) ; Effictive hor stress 170 syy = sxx zone.stress.xx(zp) = sxx zone.stress.yy(zp) = syy end_if end_if end_loop end @HorStress('MSoil','Soil') @HorStress('RAP','Soil') ;---------------------------------------model cycle 100 model solve ratio 1e-5 ;---------------------------------------model save 'Ini_Stress' ; ; ; ; ; ; ; ; ========================================================================== ========================================================================== Subroutine name: " Constitutive_Model.dat" Programing language: FLAC3D commands and FISH Description: This subroutine imports soil properties and assign the selected constitutive model (MC or PH) to soil zones ========================================================================== ========================================================================== model restore 'Ini_Stress.sav' ;--------------------------------------------------------------------------------zone gridpoint initialize displacement (0,0,0) zone gridpoint initialize velocity (0,0,0) ;--------------------------------------------------------------------------------;Mohr-Coulomb Model ;--------------------------------------------------------------------------------zone cmodel assign mohr-coulomb range group 'MSoil' slot 'Soil' zone cmodel assign mohr-coulomb range group 'RAP' slot 'Soil' ;---------------------------------------table 'TB_cSoil' delete table 'TB_cSoil' add 0.0 150 ... 6.0 50 ... 30.0 50 ;---------------------------------------; NC friction angle table 'TB_fiSoil_NC' delete table 'TB_fiSoil_NC' import 'Input_NCfiSoil_modified' ;---------------------------------------; Importing the bluk and shear moduli tables calculated from CPT data table 'TB_UnitWeight' delete table 'TB_KSoil' delete table 'TB_GSoil' delete table 'TB_UnitWeight' import 'Input_UnitWeight' table 'TB_KSoil' import 'Input_Bulk' table 'TB_GSoil' import 'Input_Shear' ;---------------------------------------fish define ApplyMC(gpname,sl) local zp = zone.head local zdep local gammaSz local KSz local GSz local cSz local fiSz ;-----------------------------------loop foreach zp zone.list 171 if zone.model(zp) # 'null' then if zone.group(zp,sl) = gpname then zdep = -1*zone.pos.z(zp) ;-----------------------;Elastic properties gammaSz = table('TB_UnitWeight',zdep) KSz = table('TB_KSoil',zdep) GSz = table('TB_GSoil',zdep) ;-----------------------cSz = table('TB_cSoil',zdep) fiSz = table('TB_fiSoil_NC',zdep) ;-----------------------zone.prop(zp,'density') = gammaSz/grav zone.prop(zp,'bulk') = KSz zone.prop(zp,'shear') = GSz zone.prop(zp,'cohesion') = cSz zone.prop(zp,'fric') = fiSz end_if end_if end_loop end ;--------------------------------------------------------------------------------; Plastic-Hardening Model ;--------------------------------------------------------------------------------table 'TB_GSoil_max' delete table 'TB_OCR' delete table 'TB_GSoil_max' import 'Input_Shear_max' table 'TB_OCR' import 'Input_OCR' ;---------------------------------------table 'TB_GRAP_max' delete table 'TB_GRAP_max' add 00.0 21.378e6 ... 2.50 21.378e6 ... 2.51 17.530e6 ... 30.0 17.530e6 ;---------------------------------------table 'TB_g07' delete table 'TB_g07' add 00.0 2.30e-4 ... 18.0 2.30e-4 ... 30.0 2.55e-4 ;---------------------------------------fish define ApplyPH(gpname,sl,EurEmax,E50Eur,EoedMult,pref,m,Rf,alpha,multK,smstrflag,soilflag) local zp = zone.head local zid local zdep local sig1 local sig2 local sig3 local fri local dil local coh local ten local nu local Gmax local Emax local Emax_ref local g07 local Eur_ref local Eur local E50_ref local E50 172 local Eoed_ref local OCR local coef local cosf local sinf local cotf local num local denum local Hc local Knc local K1 local K2 local Kp ;-------------loop foreach zp zone.list if zone.model(zp) # 'null' then if zone.model(zp) = 'mohr-coulomb' then if zone.group(zp,sl) = gpname then zid = zone.id(zp) zdep = -1*zone.pos.z(zp) if soilflag = 'Granular' then Gmax = table('TB_GRAP_max',zdep) OCR = 1.0 else Gmax = table('TB_GSoil_max',zdep) OCR = table('TB_OCR',zdep) end_if nu = zone.prop(zp,'poisson') Emax = Gmax*(2.*(1.+nu)) fri = zone.prop(zp,'friction') coh = zone.prop(zp,'cohesion') if coh < 0.1 then coh = 0.1 end_if dil = zone.prop(zp,'dilation') ten = zone.prop(zp,'tension') if smstrflag = 1 then if soilflag = 'Granular' then g07 = 1.1e-4 else g07 = table('TB_g07',zdep) end_if end_if sinf = math.sin(fri*math.degrad) cosf = math.cos(fri*math.degrad) sig1 = zone.stress.min(zp)+zone.pp(zp) sig2 = zone.stress.int(zp)+zone.pp(zp) sig3 = zone.stress.max(zp)+zone.pp(zp) cotf = cosf/sinf num = coh*cotf-sig3 denum = coh*cotf+pref coef = (num/denum)^m ;-----------------------------Eur = Emax/EurEmax E50 = Eur/E50Eur Eur_ref = Eur/coef E50_ref = E50/coef Eoed_ref = E50_ref*EoedMult Emax_ref = Emax/coef ;-----------------------------; calculate _HC 173 Knc = math.max(1-sinf,nu/(1.-nu)+1e-4) K1 = Eur_ref/(3.*(1.-2.*nu)) K2 = Eoed_ref*(1.+2.*Knc)/3. if K1 > K2 then Kp = K1*K2/(K1-K2) end_if Hc = multK*Kp ; change to Hardening soil zone.model(zp) = 'plastic-hardening' ;-----------------------------if smstrflag = 1 then zone.prop(zp,'flag-smallstrain') = 'on' zone.prop(zp,'stiffness-0-reference') = Emax_ref zone.prop(zp,'strain-70') = g07 end_if ;-----------------------------zone.prop(zp,'stiffness-50-reference') = E50_ref zone.prop(zp,'stiffness-ur-reference') = Eur_ref zone.prop(zp,'stiffness-oedometer-reference') = Eoed_ref zone.prop(zp,'pressure-reference') = pref zone.prop(zp,'exponent') = m zone.prop(zp,'poisson') = nu zone.prop(zp,'failure-ratio') = Rf zone.prop(zp,'over-consolidation-ratio') = OCR ;---------------------------zone.prop(zp,'friction') = fri zone.prop(zp,'cohesion') = coh zone.prop(zp,'dilation') = dil zone.prop(zp,'tension') = ten ;---------------------------zone.prop(zp,'stress-1-effective') = sig1 zone.prop(zp,'stress-2-effective') = sig2 zone.prop(zp,'stress-3-effective') = sig3 ;---------------------------zone.prop(zp,'constant-alpha') = alpha zone.prop(zp,'stiffness-cap-hardening') = Hc zone.prop(zp,'coefficient-normally-consolidation') = Knc end_if end_if end_if end_loop end ;--------------------------------------------------------------------------------; Applying Constitutive Model ;--------------------------------------------------------------------------------define Apply_Constitutive_Model if Constit_model = 'MC' then command @ApplyMC('MSoil','Soil') @ApplyMC('RAP','Soil') ;---------------------------model cycle 100 model solve ratio 1e-5 end_command else if Constit_model = 'PH' then command @ApplyMC('MSoil','Soil') @ApplyMC('RAP','Soil') ;---------------------------model cycle 100 model solve ratio 1e-5 174 ;---------------------------zone gridpoint initialize displacement (0,0,0) zone gridpoint initialize velocity (0,0,0) ;---------------------------;ApplyPH(gpname,sl,EurEmax,E50Eur,EoedMult,pref,m,Rf,alpha,multK,smstrflag,soilflag) @ApplyPH('MSoil','Soil',3.,3.,1.0,2116.22,0.99,0.9,1.3,0.9,1,'Cohesive') @ApplyPH('RAP','Soil',3.,3.,1.0,2116.22,0.99,0.9,1.3,0.9,1,'Cohesive') ;---------------------------model cycle 100 model solve ratio 1e-5 end_command end_if end ;---------------------------------------@Apply_Constitutive_Model ;---------------------------------------model save [Constit_model+'_Model'] ; ; ; ; ; ; ; ; ; ; ========================================================================== ========================================================================== Subroutine name: " RAP_Excavate.dat" Programing language: FLAC3D commands and FISH Description: This subroutine excavate a borehole by relaxing the zones at the RAP element location and defines FISH functions to calculate the coefficient of the lateral earth pressure and lateral displacement during different modeling stages ========================================================================== ========================================================================== model restore [Constit_model+'_Model'] ;--------------------------------------------------------------------------------zone gridpoint initialize displacement (0,0,0) zone gridpoint initialize velocity (0,0,0) ;--------------------------------------------------------------------------------;Large strain option ;--------------------------------------------------------------------------------model largestrain off ;--------------------------------------------------------------------------------; Saving K and displacement history in tables ;--------------------------------------------------------------------------------;DEfining groups to calculate K zone group 'Z_Khist00' slot 'Out' range cylinder ... end-1 (0,0,0) end-2 (0,0,[zPl-4*rPl]) radius 0.0 0.2 zone group 'Z_Khist15' slot 'Out' range cylinder ... end-1 (0,0,0) end-2 (0,0,[zPl-4*rPl]) radius 1.5 1.75 zone group 'Z_Khist25' slot 'Out' range cylinder ... end-1 (0,0,0) end-2 (0,0,[zPl-4*rPl]) radius 2.38 2.55 zone group 'Z_Khist50' slot 'Out' range cylinder ... end-1 (0,0,0) end-2 (0,0,[zPl-4*rPl]) radius 4.9 5.1 ;---------------------------------------;Define map for the selected zones and calculating the averaging parameter ;Extra data storage -> 1: radius of zone centroid, 2: effective vertical stress, 3: K fish define ZMap local zp = zone.head local zr local zz local effszz global izn = 0 local i00 = 0 local i15 = 0 local i25 = 0 175 local i50 = 0 global div15 global div25 global div50 global Map00 = map global ZoneID00 = map global Map15 = map global ZoneID15 = map global Map25 = map global ZoneID25 = map global Map50 = map global ZoneID50 = map ;-----------------------------------loop foreach zp zone.list ;Counting the number of zones in z-direction if zone.group(zp,'Out') = 'Z_Khist00' then izn = izn + 1 ;Initializing tables to store z-K values zz = math.round(100*zone.pos.z(zp))/100 table('K00',zz) = 0.0 table('K15',zz) = 0.0 table('K25',zz) = 0.0 table('K50',zz) = 0.0 end_if ;creating a map of zones to store z-K at the center of if zone.group(zp,'Out') = 'Z_Khist00' then i00 = i00 + 1 zr = math.sqrt(zone.pos.x(zp)^2+zone.pos.y(zp)^2) effszz = zone.stress.zz(zp)+zone.pp(zp) ;Effective zz = math.round(100*zone.pos.z(zp))/100 zone.extra(zp,1) = zr zone.extra(zp,2) = effszz Map00(zone.id(zp)) = zp ZoneID00(zone.id(zp)) = zz end_if ;creating a map of zones to store z-K at 1.5 ft if zone.group(zp,'Out') = 'Z_Khist15' then i15 = i15 + 1 zr = math.sqrt(zone.pos.x(zp)^2+zone.pos.y(zp)^2) effszz = zone.stress.zz(zp)+zone.pp(zp) ;Effective zz = math.round(100*zone.pos.z(zp))/100 zone.extra(zp,1) = zr zone.extra(zp,2) = effszz Map15(zone.id(zp)) = zp ZoneID15(zone.id(zp)) = zz end_if ;creating a map of zones to store z-K at 2.5 ft if zone.group(zp,'Out') = 'Z_Khist25' then i25 = i25 + 1 zr = math.sqrt(zone.pos.x(zp)^2+zone.pos.y(zp)^2) effszz = zone.stress.zz(zp)+zone.pp(zp) ;Effective zz = math.round(100*zone.pos.z(zp))/100 zone.extra(zp,1) = zr zone.extra(zp,2) = effszz Map25(zone.id(zp)) = zp ZoneID25(zone.id(zp)) = zz end_if ;creating a map of zones to store z-K at 5.0 ft if zone.group(zp,'Out') = 'Z_Khist50' then i50 = i50 + 1 zr = math.sqrt(zone.pos.x(zp)^2+zone.pos.y(zp)^2) RAP vertical stress vertical stress vertical stress 176 effszz = zone.stress.zz(zp)+zone.pp(zp) ;Effective vertical stress zz = math.round(100*zone.pos.z(zp))/100 zone.extra(zp,1) = zr zone.extra(zp,2) = effszz Map50(zone.id(zp)) = zp ZoneID50(zone.id(zp)) = zz end_if end_loop ;number of selested zone at the same depth for averaging at 1.5 ft div15 = i15/izn ;number of selested zone at the same depth for averaging at 2.5 ft div25 = i25/izn ;number of selested zone at the same depth for averaging at 5.0 ft div50 = i50/izn end @ZMap ;---------------------------------------;Calculating radial effective stress and K in radial direction fish define K_hist local zp local sxx local syy local sxy local effszz local effsrr local zx local zy local zr local zz local zK ;-----------------------------------;Reset K00 table loop local ii (1,table.size('K00')) table.y('K00',ii) = 0.0 end_loop ;Populate the K00 table with K values ;-----------------------------------loop foreach zp Map00 sxx = zone.stress.xx(zp) syy = zone.stress.yy(zp) sxy = zone.stress.xy(zp) ;effszz = zone.stress.zz(zp)+zone.pp(zp) ;Effective vertical stress effszz = zone.extra(zp,2) ;Effective vertical stress ;-------------------------------zx = zone.pos.x(zp) zy = zone.pos.y(zp) zr = zone.extra(zp,1) ;Stress trasnformaation effsrr = (sxx*zx^2+syy*zy^2+2*sxy*zx*zy)/zr^2+zone.pp(zp) zK = effsrr/effszz ;curent K value zone.extra(zp,3) = zK ;curent K value ;-------------------------------zz = ZoneID00(zone.id(zp)) table('K00',zz) = zK end_loop ;-----------------------------------;Reset K15 table loop local jj (1,table.size('K15')) table.y('K15',jj) = 0.0 end_loop ;Populate the K15 table with K values 177 ;-----------------------------------loop foreach zp Map15 sxx = zone.stress.xx(zp) syy = zone.stress.yy(zp) sxy = zone.stress.xy(zp) ;effszz = zone.stress.zz(zp)+zone.pp(zp) ;Effective vertical stress effszz = zone.extra(zp,2) ;Effective vertical stress ;-------------------------------zx = zone.pos.x(zp) zy = zone.pos.y(zp) zr = zone.extra(zp,1) ;Stress trasnformaation effsrr = (sxx*zx^2+syy*zy^2+2*sxy*zx*zy)/zr^2+zone.pp(zp) zK = effsrr/effszz ;curent K value zone.extra(zp,3) = zK ;curent K value ;-------------------------------zz = ZoneID15(zone.id(zp)) table('K15',zz) = table('K15',zz) + zK/div15 end_loop ;-------;Reset K25 table loop local kk (1,table.size('K25')) table.y('K25',kk) = 0.0 end_loop ;Populate the K25 table with K values ;-----------------------------------loop foreach zp Map25 sxx = zone.stress.xx(zp) syy = zone.stress.yy(zp) sxy = zone.stress.xy(zp) ;effszz = zone.stress.zz(zp)+zone.pp(zp) ;Effective vertical stress effszz = zone.extra(zp,2) ;Effective vertical stress ;-------------------------------zx = zone.pos.x(zp) zy = zone.pos.y(zp) zr = zone.extra(zp,1) ;Stress trasnformaation effsrr = (sxx*zx^2+syy*zy^2+2*sxy*zx*zy)/zr^2+zone.pp(zp) zK = effsrr/effszz ;curent K value zone.extra(zp,3) = zK ;curent K value ;-------------------------------zz = ZoneID25(zone.id(zp)) table('K25',zz) = table('K25',zz) + zK/div25 end_loop ;-----------------------------------;Reset K50 table loop local ll (1,table.size('K50')) table.y('K50',ll) = 0.0 end_loop ; Populate the K50 table with K values ;-----------------------------------loop foreach zp Map50 sxx = zone.stress.xx(zp) syy = zone.stress.yy(zp) sxy = zone.stress.xy(zp) effszz = zone.extra(zp,2) ;Effective vertical stress ;-------------------------------zx = zone.pos.x(zp) zy = zone.pos.y(zp) zr = zone.extra(zp,1) ; Stress trasnformaation 178 effsrr = (sxx*zx^2+syy*zy^2+2*sxy*zx*zy)/zr^2+zone.pp(zp) zK = effsrr/effszz ;curent K value zone.extra(zp,3) = zK ;curent K value ;-------------------------------zz = ZoneID50(zone.id(zp)) table('K50',zz) = table('K50',zz) + zK/div50 end_loop end ;---------------------------------------@K_hist table 'K00' export ['K00_verification'+Constit_model] truncate csv table 'K25' export ['K25_verification'+Constit_model] truncate csv ;---------------------------------------; Calculating Radial displacement at a specified gridpoint fish define Rad_disp(x,y,z) local gp local gpdispx local gpdispy local gpx local gpy global gpdisp ;-----------------------------------gp = gp.near(x,y,z) gpdispx = gp.disp.x(gp) gpdispy = gp.disp.y(gp) gpx = gp.pos.x(gp) gpy = gp.pos.y(gp) gpdisp = (gpdispx*gpx+gpdispy*gpy)/math.sqrt(gpx^2+gpy^2) end ;---------------------------------------; Calculating Radial displacement at a specified x and y location ; and save data in a table at depth interavl of half of the lift thickness fish define Rad_disp_RAP_table(TBname,x,y,LfTh) local jj = int(2*math.round(-1*zPl/LfTh)+1) local ll local zz local gp local gpdispx local gpdispy local gpx local gpy local gpz local rad_disp ;-----------------------------------loop ll (1,jj) zz = zPl+(ll-1)*LfTh/2 gp = gp.near(x,y,zz) gpdispx = gp.disp.x(gp) gpdispy = gp.disp.y(gp) gpx = gp.pos.x(gp) gpy = gp.pos.y(gp) gpz = gp.pos.z(gp) rad_disp = (gpdispx*gpx+gpdispy*gpy)/math.sqrt(gpx^2+gpy^2) table(TBname,gpz) = rad_disp end_loop end ;Radial displacement at a RAP interface table 'Rad_disp' delete @Rad_disp_RAP_table('Rad_disp',[rPl*sin45],[rPl*sin45],1.0) table 'Rad_disp' export ... ['Disp_r15_th45_Verification_'+Constit_model] truncate csv 179 ;--------------------------------------------------------------------------------; Excavating the borehole ;--------------------------------------------------------------------------------zone relax excavate name 'EX_RAP' step 5000 range group 'RAP' slot 'Soil' model cycle 5000 model save ['RAP_Excavate_'+Constit_model] ; ; ; ; ; ; ; ; ========================================================================== ========================================================================== Subroutine name: " Anchor_Const.dat" Programing language: FLAC3D commands and FISH Description: This subroutine places the RAP material in the excavated borehole without any RAP construction consideration ========================================================================== ========================================================================== model restore ['RAP_Excavate_'+Constit_model] [global RAP_Const = 'Anchor'] ; Construction method ;---------------------------------------zone gridpoint initialize displacement (0,0,0) zone gridpoint initialize velocity (0,0,0) ;--------------------------------------------------------------------------------;Assigning soil properties ;--------------------------------------------------------------------------------zone cmodel assign mohr-coulomb range group 'RAP' slot 'Soil' zone property density @RoRAP bulk @KRAP shear @GRAP cohesion @cRAP ... fric @fiRAP dilation @saiRAP range group 'RAP' slot 'Soil' ;---------------------------------------model cycle 100 model solve ratio 1e-5 ;---------------------------------------fish define Anchor_PH local zp local sig_xx local sig_xy local sig_xz local sig_yy local sig_yz local sig_zz if Constit_model = 'PH' then ApplyPH('RAP','Soil',3.,3.,1.0,2116.22,0.48,0.9,1.3,0.9,1,'Granular') command model cycle 100 model solve ratio 1e-5 end_command end_if loop foreach zp zone.list if zone.model(zp) # 'null' then if zone.group(zp,'Soil') = 'RAP' then ;Extra data storage -> 4: Sigma_xx, 5: Sigma_xy, 6: Sigma_xz, ; 7: Sigma_yy, 8: Sigma_yz, 9: Sigma_zz sig_xx = zone.stress.xx(zp) sig_xy = zone.stress.xy(zp) sig_xz = zone.stress.xz(zp) sig_yy = zone.stress.yy(zp) sig_yz = zone.stress.yz(zp) sig_zz = zone.stress.zz(zp) ;-----------------------zone.extra(zp,4) = sig_xx zone.extra(zp,5) = sig_xy zone.extra(zp,6) = sig_xz 180 zone.extra(zp,7) = sig_yy zone.extra(zp,8) = sig_yz zone.extra(zp,9) = sig_zz end_if end_if end_loop end @Anchor_PH ;---------------------------------------; Saving K @K_hist table 'K00' export ['K00_'+RAP_Const+'_final_'+Constit_model] truncate csv table 'K15' export ['K15_'+RAP_Const+'_final_'+Constit_model] truncate csv table 'K25' export ['K25_'+RAP_Const+'_final_'+Constit_model] truncate csv table 'K50' export ['K50_'+RAP_Const+'_final_'+Constit_model] truncate csv ;--------------------------------------------------------------------------------; Defining function to assign RAP properties during RAPs construction ;--------------------------------------------------------------------------------fish define Place_RAP(gpname,sl) local zp local sig_xx local sig_xy local sig_xz local sig_yy local sig_yz local sig_zz ;-----------------------------------command zone cmodel assign mohr-coulomb range group @gpname slot @sl zone property density @RoRAP bulk @KRAP shear @GRAP cohesion @cRAP ... fric @fiRAP dilation @saiRAP range group @gpname slot @sl end_command ;-----------------------------------loop foreach zp zone.list if zone.group(zp,sl) = gpname then ;---------------------------sig_xx = zone.extra(zp,4) sig_xy = zone.extra(zp,5) sig_xz = zone.extra(zp,6) sig_yy = zone.extra(zp,7) sig_yz = zone.extra(zp,8) sig_zz = zone.extra(zp,9) ;---------------------------zone.stress.xx(zp) = sig_xx zone.stress.xy(zp) = sig_xy zone.stress.xz(zp) = sig_xz zone.stress.yy(zp) = sig_yy zone.stress.yz(zp) = sig_yz zone.stress.zz(zp) = sig_zz end_if end_loop ;-----------------------------------if Constit_model = 'PH' then ApplyPH(gpname,sl,3.,3.,1.0,2116.22,0.48,0.9,1.3,0.9,1,'Granular') end_if end ;---------------------------------------model save ['RAP_const_Anchor_'+Constit_model] 181 ; ; ; ; ; ; ; ; ; ========================================================================== ========================================================================== Subroutine name: " RAP_Const_K_Blade.dat" Programing language: FLAC3D commands and FISH Description: This subroutine simulates the RAP construction by a simple initialization of the effective lateral stress to the measured value ========================================================================== ========================================================================== model restore ['RAP_const_Anchor_'+Constit_model] ;--------------------------------------------------------------------------------zone gridpoint initialize displacement (0,0,0) zone gridpoint initialize velocity (0,0,0) ;--------------------------------------------------------------------------------;Correct horizontal stress ;--------------------------------------------------------------------------------[RAP_Const = 'K_Blade'] fish define HorStress_K_Blade(gpname,sl) local zp = zone.head local zdep local szz local sxx local syy local K0 ;-----------------------------------loop foreach zp zone.list if zone.model(zp) # 'null' then if zone.group(zp,sl) = gpname then zdep = -1*zone.pos.z(zp) K0 = table('TB_K25_target',zdep) szz = zone.stress.zz(zp)+zone.pp(zp) ;Effective vertical stress sxx = K0*szz-zone.pp(zp) ;Correct effective hor. stress syy = sxx zone.stress.xx(zp) = sxx zone.stress.yy(zp) = syy end_if end_if end_loop end @HorStress_K_Blade('MSoil','Soil') @HorStress_K_Blade('RAP','Soil') ;---------------------------------------model cycle 100 model solve ratio 1e-5 ;---------------------------------------model save ['RAP_const_K_Blade_'+Constit_model] ; ; ; ; ; ; ; ; ========================================================================== ========================================================================== Subroutine name: " RAP_Const_Cavity_velocity.dat" Programing language: FLAC3D commands and FISH Description: This subroutine simulates the RAP construction by the cavity expansion technique ========================================================================== ========================================================================== model restore ['RAP_const_Anchor_'+Constit_model] [global RAP_Const = 'Cav_Vel'] ; Construction method ;---------------------------------------zone gridpoint initialize displacement (0,0,0) 182 zone gridpoint initialize velocity (0,0,0) ;--------------------------------------------------------------------------------; Cavity expansion ;--------------------------------------------------------------------------------zone cmodel assign null range group 'RAP' slot 'Soil' table 'ramp' delete table 'ramp' add [global.step] 0.0 ... [global.step+10000] 5.0e-7 ... [global.step+500000] 5.0e-7 ;---------------------------------------zone face apply velocity-normal -1 table 'ramp' range group 'FC_SideRAP' slot 'Int' zone face apply velocity-normal -2 table 'ramp' range group 'FC_Plate' slot 'Int' ;---------------------------------------zone gridpoint group 'GP_Rside' slot 'Int' range cylinder end-1 (0,0,[zPl-0.01]) ... end-2 (0,0,[zPl+0.01]) radius [rPl-0.05] [rPl+0.05] zone gridpoint free velocity range group 'GP_Rside' slot 'Int' ;---------------------------------------zone gridpoint group 'GP_Xside' slot 'Int' range cylinder end-1 (@rPl,0,0.01) ... end-2 (@rPl,0,[zPl-0.01]) radius 0.05 zone gridpoint group 'GP_Yside' slot 'Int' range cylinder end-1 (0,@rPl,0.01) ... end-2 (0,@rPl,[zPl-0.01]) radius 0.05 ;---------------------------------------zone gridpoint fix velocity-y range group 'GP_Xside' slot 'Int' zone gridpoint fix velocity-y range group 'GP_Yside' slot 'Int' ;---------------------------------------table 'Rad_disp45_z50' delete ;---------------------------------------fish define CavityExpanV(DelStep,nRun) local NameFile ;-----------------------------------loop local ii (1,nRun) NameFile = RAP_Const+'_Step_'+string(global.step)+'_'+Constit_model command model cycle @DelStep ;---------------------------------------@K_hist table 'K15' export ['K15_'+NameFile] truncate csv table 'K25' export ['K25_'+NameFile] truncate csv table 'K50' export ['K50_'+NameFile] truncate csv ;---------------------------------------model save [NameFile+'.sav'] ;---------------------------------------;Radial displacement at a RAP interface, z = zPl/2 @Rad_disp([rPl*sin45],[rPl*sin45],[zPl/2) table 'Rad_disp45_z50' insert ([global.step],@gpdisp) end_command end_loop command table 'Rad_disp45_z50' export 'Cav_Vel_Disp_Key' truncate csv end_command end ;Run the cavity expansion for 4 times with 10000 steps in each run @CavityExpanV(10000,5) ;--------------------------------------------------------------------------------;RAP placement ;--------------------------------------------------------------------------------zone face apply-remove velocity-normal range group 'FC_SideRAP' slot 'Int' zone face apply-remove velocity-normal range group 'FC_Plate' slot 'Int' ;--------------------------------------------------------------------------------; Assigning the constitutive model and solving the model ;--------------------------------------------------------------------------------- 183 ;Constitutive Model @Place_RAP('RAP','Soil') ;---------------------------------------model cycle 100 model solve ratio 1e-5 ;--------------------------------------------------------------------------------; Saving K and displacement ;--------------------------------------------------------------------------------@K_hist table 'K15' export ['K15_'+RAP_Const+'_final_'+Constit_model] truncate csv table 'K25' export ['K25_'+RAP_Const+'_final_'+Constit_model] truncate csv table 'K50' export ['K50_'+RAP_Const+'_final_'+Constit_model] truncate csv ;---------------------------------------; Radial displacement at a RAP interface table 'Rad_disp' delete @Rad_disp_RAP_table('Rad_disp',[rPl*sin45],[rPl*sin45],1.0) table 'Rad_disp' export ... ['Disp_r15_th45_'+RAP_Const+'_final_'+Constit_model] truncate csv ;---------------------------------------model save ['RAP_const_Cav_Vel_'+Constit_model] ; ; ; ; ; ; ; ; ========================================================================== ========================================================================== Subroutine name: " RAP_Const_Lift_Pressure.dat" Programing language: FLAC3D commands and FISH Description: This subroutine simulates the RAP construction by the lift compaction technique ========================================================================== ========================================================================== model restore ['RAP_const_Anchor_'+Constit_model] [global RAP_Const = 'Lift_Press'] ; Construction method ;---------------------------------------zone gridpoint initialize displacement (0,0,0) zone gridpoint initialize velocity (0,0,0) ;--------------------------------------------------------------------------------;RAP construction in lifts ;--------------------------------------------------------------------------------zone cmodel assign null range group 'RAP' slot 'Soil' ;Building RAP with sleeve for borehole stability zone face apply velocity-normal 0 range cylinder ... end-1 (0,0,0) end-2 (0,0,@zPl) radius [rPl-0.01] [rPl+0.01] ;---------------------------------------fish defin LiftP(Press,LfTh) local ii = 1 local jj = int(math.round(-1*zPl/LfTh)) local kk = 0 local ll = 0 local zl local NameFile local idx ;-----------------------------------loop ll (1,jj) ;-------------------------------zl = zPl + ll*LfTh ;-------------------------------idx = int((ii-kk)/1) if idx = 1 then NameFile = RAP_Const+'_'+string(Press)+'_Lift_'+... string(ll)+'_'+Constit_model kk = ii 184 end_if command zone face apply-remove velocity-normal range cylinder end-1 ... (0,0,[zl-LfTh]) end-2 (0,0,[zl]) radius [rPl-0.01] [rPl+0.01] ;---------------------------;A table to increase pressure applied to each lift over 10,000 steps table 'ramp_p' delete table 'ramp_p' add [global.step] 0.0 ... [global.step+2000] @Press ... [global.step+7000] @Press ;---------------------------;Placing RAP lift zone group ['RAP_Lift_'+string(ll)] slot 'RAP_Const' range ... group 'RAP' slot 'Soil' position-z [zl-LfTh] @zl [Place_RAP('RAP_Lift_'+string(ll),'RAP_Const')] ;Applying compaction pressure on RAP lift zone face group ['FC_lift_'+string(ll)] slot 'Int' range cylinder ... end-1 (0,0,[zl-0.1]) end-2 (0,0,[zl+0.1]) radius [rPl+0.01] zone face apply stress-normal -1 table 'ramp_p' range ... group ['FC_lift_'+string(ll)] slot 'Int' model cycle 7000 ;Removing compaction pressure on RAP lift zone face apply-remove stress-normal range ... group ['FC_lift_'+string(ll)] slot 'Int' end_command ;----------------if idx = 1 then command @K_hist table 'K15' export ['K15_'+NameFile] truncate csv table 'K25' export ['K25_'+NameFile] truncate csv table 'K50' export ['K50_'+NameFile] truncate csv model save [NameFile+'.sav'] ;----------------;Radial displacement at a RAP interface table 'Rad_disp' delete @Rad_disp_RAP_table('Rad_disp',[rPl*sin45],[rPl*sin45],1.0) table 'Rad_disp' export ['Disp_r15_th45_'+NameFile] truncate csv end_command end_if ii = ii+1 end_loop command model solve ratio 1e-5 ;---------------------------------------@K_hist table 'K15' export ... ['K15_'+RAP_Const+'_'+string(Press)+'_final_'+Constit_model] truncate csv table 'K25' export ... ['K25_'+RAP_Const+'_'+string(Press)+'_final_'+Constit_model] truncate csv table 'K50' export ... ['K50_'+RAP_Const+'_'+string(Press)+'_final_'+Constit_model] truncate csv ;---------------------------------------;Radial displacement at a RAP interface table 'Rad_disp' delete @Rad_disp_RAP_table('Rad_disp',[rPl*sin45],[rPl*sin45],1.0) table 'Rad_disp' export ... ['Disp_r15_th45_'+RAP_Const+'_'+string(Press)+'_final_'+Constit_model] truncate csv ;---------------------------------------model save ['RAP_const_Lift_Press_'+Constit_model] end_command 185 end ;Applying 22e3 pcf pressure on top of each 1 ft lift @LiftP(22.e3,@LiftTh) # # # # # # # # ========================================================================== ========================================================================== Subroutine name: "RAP_Const_Cavity_variable_velocity.py" Programing language: Python Description: This subroutine simulates the RAP construction process by the cavity expansion with variable displacement technique ========================================================================== ========================================================================== import itasca as it import numpy as np it.command("python-reset-state false") #----------------------------------------------------------------------------# Getting the model and stress state data from FLAC3D # and defining loading functions for linearization it.command(""" model restore ['RAP_const_Anchor_'+Constit_model] [global RAP_Const = 'Cav_Variable'] [global nstep = 0] ;---------------------------------------zone gridpoint initialize displacement (0,0,0) zone gridpoint initialize velocity (0,0,0) ;------------------------------------------------------------------------; Removing RAP material and restraining the borehole shaft and bottom zone cmodel assign null range group 'RAP' slot 'Soil' ;Building RAP with sleeve for borehole stability zone face apply velocity-normal 0 range cylinder ... end-1 (0,0,0) end-2 (0,0,@zPl) radius [rPl-0.01] [rPl+0.01] zone face apply velocity-normal 0 range cylinder ... end-1 (0,0,0) end-2 (0,0,@zPl) radius 0 [rPl+0.01] ; Solving the model model solve ratio 1e-5 ;------------------------------------------------------------------------; Getting the initial effective vertical stress in native soil table 'TB_EffSigv0' delete fish define EffSigv0 local zp local zz loop foreach zp Map00 zz = zone.pos.z(zp) if zz > zPl then table('TB_EffSigv0',zz) = zone.extra(zp,2) end_if end_loop end @EffSigv0 table 'TB_EffSigv0' export ['EffSigv0_'+Constit_model] truncate csv ;------------------------------------------------------------------------; Defining RAP lift groups fish define RAPliftgroup loop local ii (1,math.round(-1*zPl)) io.out(ii) command zone face group ['FC_lift_'+string(ii)] slot 'RAPlift' ... internal range cylinder end-1 (0,0,[zPl+ii] end-2 ... (0,0,[zPl+ii-1]) radius [rPl-0.01] [rPl+0.01] end_command 186 end_loop end @RAPliftgroup ;------------------------------------------------------------------------; Defining some gridpoint group zone gridpoint group 'GP_Rside' slot 'Int' range cylinder end-1 ... (0,0,[zPl-0.01]) end-2 (0,0,[zPl+0.01]) radius [rPl-0.05] [rPl+0.05] zone gridpoint group 'GP_Xside' slot 'Int' range cylinder end-1 ... (@rPl,0,0.01) end-2 (@rPl,0,[zPl-0.01]) radius 0.05 zone gridpoint group 'GP_Yside' slot 'Int' range cylinder end-1 ... (0,@rPl,0.01) end-2 (0,@rPl,[zPl-0.01]) radius 0.05 ;------------------------------------------------------------------------; Defining fish function to apply load/displacement at a specific zone face ; Defining function to apply displacement at each lift fish define UnitDisp(UDisp,DispStep,gpname,sl) command ; Defining table to apply and remove velocity gradually table 'RampVel' delete table 'RampVel' add [global.step] 0.0 ... [global.step+DispStep] 36.0e-7 ... [global.step+2*DispStep] 0.0 zone face apply velocity-normal [-1*UDisp] table 'RampVel' ... range group [gpname] slot [sl] end_command ; Applying the vertical velocity to the bottom of the borehole if gpname = 'FC_lift_1' then command zone face apply velocity-normal [-2*UDisp] table 'RampVel' ... range group 'FC_Plate' slot 'Int' end_command end_if ; BC revision command zone gridpoint free velocity range group 'GP_Rside' slot 'Int' zone gridpoint fix velocity-y range group 'GP_Xside' slot 'Int' zone gridpoint fix velocity-y range group 'GP_Yside' slot 'Int' end_command end ; Saving the lateral earth pressure coefficient at r = 2.5 ft @K_hist table 'K25' export 'K25_step' truncate csv ; Saving the model as step 0 model save ... ['RAP_const_Cav_variable_step_'+string(nstep)+'_'+Constit_model] """) #----------------------------------------------------------------------------# Inputs #----------------------------------------------------------------------------# Initial scaling factor for the pattern alpha_0 = 2.5 alpha = alpha_0 alpha_save = [alpha_0] scale = lambda NF: 0.25-NF*0.23 # getting the model data nstep = it.fish.get('nstep') Constit_model = it.fish.get('Constit_model') zPl = it.fish.get('zPl') rPl = it.fish.get('rPl') liftTH = it.fish.get('LiftTh') # x and y coordinate to get the displacement for linearization Disp_loc = rPl*np.sin(np.pi/4) # x=y at 45 degrees 187 # getting the effective vertical stress at target horiz. stress from FLAC3D tables effsigv0 = np.loadtxt("EffSigv0_"+Constit_model+".tab",skiprows = 2,delimiter=",") n_lift = np.round(np.abs(zPl)/liftTH).astype(int) # number of lifts in the RAP n_z_lift = (len(effsigv0[:,0])/n_lift).astype(int) # number of zones per lift # calculating effective target horz. stress at r = 2.5ft K25_target = np.loadtxt("Input_K25-target.tab",skiprows = 2) K25_target[:,0] = -1*K25_target[:,0] K25_target = np.sort(K25_target, axis=0) K25_target_interpolated = np.interp(effsigv0[:,0], K25_target[:,0], K25_target[:,1]) effsigh25_target = effsigv0[:,1]*K25_target_interpolated sum_effsigh25_target = sum(effsigh25_target) # initial displacement pattern for cavity expansion pat = np.arange(n_lift,1,-1) pat = np.append(pat,0.1) pat = pat/np.linalg.norm(pat) # Normalization: unit vector pat_save = pat[:,np.newaxis] print('End of Input!') # getting the K at r = 2.5ft at the end of step 0 K25_step = np.loadtxt("K25_step.tab",skiprows = 2,delimiter=",") K25_step = K25_step[K25_step[:,0]>zPl] # Calculating the difference beween the target lateral stress and # the lateral stress at r = 2.5 ft at the end of step 0 effsigh25 = K25_step[:,1]*effsigv0[:,1] effsigh25_save = effsigh25[:,np.newaxis] Delta_y_target_0 = effsigh25_target-effsigh25 # at zones Delta_y_target_0 = np.mean(Delta_y_target_0.reshape(-1, n_z_lift), axis = 1) # at lifts Delta_y_target_solve = np.copy(Delta_y_target_0) # Initializing the normalized sum of the horizontal stress norm_hor_force = np.array([0]) #---------------------------------------------------------------------------------# Linearization to find the new displacement pattern for cavity expansion #---------------------------------------------------------------------------------for totstep in range(30): print('check point 1') # Initializing the Jacobian matrix print(f'K25_step = {K25_step}') Jac = np.zeros((n_lift,n_lift)) # Running the single lift expansion function and assembling the Jacobian matrix for i, pati in enumerate(pat): print(f'nstep = {nstep}, i = {i}, pati = {pati}') # Getting the displacement before the lift zone push out it.fish.call_function('Rad_disp',(Disp_loc,Disp_loc,zPl+0.5+i*liftTH)) Disp_step = it.fish.get('gpdisp') # setting the applied velocity to zero it.command(""" zone face apply velocity-normal 0 range group 'FC_SideRAP' slot 'Int' zone face apply velocity-normal 0 range group 'FC_Plate' slot 'Int' """) # Applying the ramped velocity to expand the cavity the specified lift it.fish.call_function('UnitDisp', (alpha*pati,1000,'FC_lift_'+str(i+1),'RAPlift')) # Cycling the model to expand the cavity at the lift and solve for equilibrium it.command(""" model cycle 2100 model solve ratio 1e-5 @K_hist table 'K25' export 'K25_temp' truncate csv """) # calculating effective target horz. stress at r = 2.5ft K25_temp = np.loadtxt("K25_temp.tab",skiprows = 2,delimiter=",") K25_temp = K25_temp[K25_temp[:,0]>zPl] 188 # Calculating the increase in lateral stress at r = 2.5 ft Delta_y = (K25_temp[:,1]-K25_step[:,1])*effsigv0[:,1] # at zones Delta_y = np.mean(Delta_y.reshape(-1, n_z_lift), axis = 1) # at lifts print(f'Delta_y = {Delta_y}') # Getting the displacement at the end of the lift zone push out it.fish.call_function('Rad_disp',(Disp_loc,Disp_loc,zPl+0.5+i*liftTH)) Disp_temp = it.fish.get('gpdisp') print(f'Disp_temp = {Disp_temp}, Disp_step = {Disp_step}') # Calculating the increase in borhole radius due to lift expansion Delta_x = Disp_temp - Disp_step print(f'Delta_x = {Delta_x}') # Assembling the Jacobian matrix Jac[:,i] = Delta_y/Delta_x # Restoring the model to the end of the last step it.command(""" model restore ['RAP_const_Cav_variable_step_'+string(nstep)+'_'+Constit_model] """) # Solving a system of linear scalar equations pat = np.linalg.solve(Jac, Delta_y_target_solve) print(f'Jac = {Jac}') print(f'Delta_y_target_solve = {Delta_y_target_solve}') print(f'pat solve = {pat}') pat[pat < 0.0] = 0.0 # set the negative values to 0 pat = pat/np.linalg.norm(pat) # normalize pat[pat < 0.01] = 0.01 # set the minimum value to 0.01 print(f'pat modified = {pat}') pat = pat/np.linalg.norm(pat) # new displacement pattern print(f'new pat = {pat}') pat_save = np.append(pat_save,pat[:,np.newaxis],axis=1) #---------------------------------------------------------------------------------# Expanding the borehole with new lateral displacement pattern #---------------------------------------------------------------------------------# setting the applied velocity to zero it.command(""" zone face apply velocity-normal 0 range group 'FC_SideRAP' slot 'Int' zone face apply velocity-normal 0 range group 'FC_Plate' slot 'Int' """) for i, pati in enumerate(pat): it.fish.call_function('UnitDisp', (alpha*pati,1000,'FC_lift_'+str(i+1),'RAPlift')) nstep += 1 # updating the step number in python it.fish.set('nstep',nstep) # updating the nstep in fish # Solving for the new displacement pattern it.command(""" [io.out(string(nstep))] model cycle 2100 model solve ratio 1e-5 ; Saving the lateral earth pressure coefficient at r = 2.5 ft @K_hist table 'K25' export 'K25_step' truncate csv table 'K15' export 'K15_step' truncate csv ;Radial displacement table 'Rad_disp' delete @Rad_disp_RAP_table('Rad_disp',[rPl*sin45],[rPl*sin45],1.0) table 'Rad_disp' export ... ['Disp_r15_th45_'+RAP_Const+'_step_'+string(nstep)+'_'+Constit_model] ... truncate csv ; Saving the model at the end of the step model save ['RAP_const_Cav_variable_step_'+string(nstep)+'_'+Constit_model] """) #---------------------------------------------------------------------------------- 189 # Getting lateral stress parameters at the end of the current step #---------------------------------------------------------------------------------# getting the K at r = 2.5ft at the end of the current step K25_step = np.loadtxt("K25_step.tab",skiprows = 2,delimiter=",") K25_step = K25_step[K25_step[:,0]>zPl] # Calculating the difference beween the target lateral stress and # the current lateral stress at r = 2.5 ft effsigh25 = K25_step[:,1]*effsigv0[:,1] effsigh25_save = np.append(effsigh25_save,effsigh25[:,np.newaxis],axis=1) Delta_y_target = effsigh25_target-effsigh25 # at zones Delta_y_target = np.mean(Delta_y_target.reshape(-1, n_z_lift), axis = 1) # at lifts Delta_y_target_solve = 1.05*effsigh25_target-effsigh25 # at zones Delta_y_target_solve = np.mean(Delta_y_target_solve.reshape(-1, n_z_lift),axis = 1) # Scale factor norm_hor_force = np.append(norm_hor_force, 1-sum(Delta_y_target)/sum(Delta_y_target_0)) if norm_hor_force[nstep] > norm_hor_force[nstep-1]: alpha = alpha*scale(norm_hor_force[nstep])/( norm_hor_force[nstep]-norm_hor_force[nstep-1]) alpha_save.append(alpha) print(f'alpha = {alpha}') if sum(effsigh25) < 1.05*sum_effsigh25_target: print('The sum of radial stress is higher than the sum of target') print(f'sum(effsigh25) = {sum(effsigh25)}, sum_effsigh25_target = {sum_effsigh25_target}') break elif abs(1-norm_hor_force[nstep])<0.04: print('NF satisfied') print(f'diff_force = {1-norm_hor_force[nstep]}') break # Saving the displacement pattern np.savetxt('Disp_Pattern_'+Constit_model+'.out', pat_save, delimiter=',') np.savetxt('effsigh25_'+Constit_model+'.out', effsigh25_save, delimiter=',') np.savetxt('effsigh25_target_'+Constit_model+'.out', effsigh25_target, delimiter=',') np.savetxt('alpha_'+Constit_model+'.out', np.array(alpha_save), delimiter=',') #------------------------------------------------------------------------------------# RAP placement #------------------------------------------------------------------------------------it.command(""" zone face apply-remove velocity-normal range group 'FC_SideRAP' slot 'Int' zone face apply-remove velocity-normal range group 'FC_Plate' slot 'Int' ; RAP placement zone cmodel assign mohr-coulomb range group 'RAP' slot 'Soil' ; Getting the total vertical stress at the bottom of the borehole zone gridpoint group 'BotRAP' slot 'Out' range group 'FC_Plate' slot 'Int' fish define VertStress local gplist global BH_bot_Vstress gplist = gp.list(gp.group(::gp.list,'Out') == 'BotRAP') BH_bot_Vstress = list.sum(gp.force.unbal.z(::gplist)) / (0.25*math.pi*rPl^2) end @VertStress ; Defining table to linearly change the RAP effective stress from zPl to ground table 'TB_EffSigv' delete table 'TB_EffSigv' add @zPl [-1*BH_bot_Vstress-wdens*grav*math.min(zPl-zwt,0)] ... 0.0 0.0 ; Initializing the RAP stress to the stress of the native soils adjacent to RAP fish define RAPStress(gpname,sl) local zp = zone.head local zdep local effszz 190 local effszz0 local szz local sxx local syy local K15 local KpRAP = (1+math.sin(math.pi*fiRAP/180))/(1-math.sin(math.pi*fiRAP/180)) local KaRAP = (1-math.sin(math.pi*fiRAP/180))/(1+math.sin(math.pi*fiRAP/180)) local KK ;-----------------------------------loop foreach zp zone.list if zone.group(zp,sl) = gpname then zdep = zone.pos.z(zp) K15 = table('K15',zdep) ; admissible stress state KK = math.min(0.98*KpRAP,K15) ;Kp is set as the maximum value for K KK = math.max(1.02*KaRAP,KK) ;Ka is set as the minimum value for K effszz = table('TB_EffSigv',zdep) ; Effective vertical stress effszz0 = table('TB_EffSigv0',zdep) ; Initial effective vertical stress szz = effszz-zone.pp(zp) ; Total vertical stress sxx = KK*effszz0-zone.pp(zp) ; Total radial stress syy = sxx ; Zone stress adjustment zone.stress.xx(zp) = sxx zone.stress.yy(zp) = syy zone.stress.zz(zp) = szz ; Extra variable adjustment for PH constitutive model zone.extra(zp,4) = sxx zone.extra(zp,5) = 0.0 zone.extra(zp,6) = 0.0 zone.extra(zp,7) = syy zone.extra(zp,8) = 0.0 zone.extra(zp,9) = szz end_if end_loop end @RAPStress('RAP','Soil') ; Assigning the constitutive model and solving the model @Place_RAP('RAP','Soil') model cycle 100 model solve ratio 1e-5 ; Saving the lateral earth pressure coefficient @K_hist table 'K15' export ['K15_'+RAP_Const+'_final_'+Constit_model] truncate csv table 'K25' export ['K25_'+RAP_Const+'_final_'+Constit_model] truncate csv table 'K50' export ['K50_'+RAP_Const+'_final_'+Constit_model] truncate csv ;Radial displacement table 'Rad_disp' delete @Rad_disp_RAP_table('Rad_disp',[rPl*sin45],[rPl*sin45],1.0) table 'Rad_disp' export ... ['Disp_r15_th45_'+RAP_Const+'_final_'+Constit_model] truncate csv ; Saving the model at the end of the step model save ['RAP_const_Cav_variable_'+Constit_model] """) ; ========================================================================== ; ========================================================================== ; Subroutine name: " Plate_Rod.dat" ; Programing language: FLAC3D commands and FISH ; Description: This subroutine places the cable and liner structural ; elements and assign their properties ; ========================================================================== 191 ; ========================================================================== model restore ['RAP_const_'+RAP_Const+'_'+Constit_model] ;---------------------------------------zone gridpoint initialize displacement (0,0,0) zone gridpoint initialize velocity (0,0,0) ;--------------------------------------------------------------------------------; Removing the extra soil deeper than 6 ft below the bottom of RAP and resetting BC ;--------------------------------------------------------------------------------zone group 'Extra_Soil' slot 'Soil' range cylinder ... end-1 (0,0,[zPl-4*rPl]) end-2 (0,0,@zbot) radius [rMod+0.1] zone cmodel assign null range group 'Extra_Soil' zone face group 'FC_Bot_mod' slot 'BC' range position-z [zPl-6.15-0.05] [zPl-6.15+0.05] zone face apply velocity (0,0,0) range group 'FC_Bot_mod' slot 'BC' ;Re-activating the zone seperation zone attach delete ; Fish function to calculate the liner and rod properties struct damping combined-local fish define LinerProp ;Liner local DeltaZMin = ClearRod global ksLiner = 10*(KSoil+GSoil*4/3)/DeltaZMin ;Liner interface stiffness global cLiner = 0.0 ;interface adhesion ;Cable global ARod = (math.pi*dRod^2)/4. ;Rod cross section area global kgCable = 2*math.pi*GSoil/(10*math.ln(1.2));Grout stiffness global pgCable = 1.2*math.pi*dRod ;Grout exposed perimeter global segRod = int(2*math.abs(zPl)) ;Rod number of segments ;Steel properties global fiSt_top = math.atan(math.tan(math.pi*fiRAP/180)*2/3)*180/math.pi global fiSt_bot = math.atan(math.tan(math.pi*fiRAP/180)*1/2)*180/math.pi ;... interface friction angle global RoSt = 489.8/grav ;Steel mass density global ESt = 29e6*12^2 ;Steel young modulus global NuSt = 0.26 ;Steel poisson's ratio global YSt = 75e3*12^2 ;Rebar yield strength end @LinerProp ; Create the uplift plate and set its associated properties ;element-type 'dkt-csth' struct liner create by-face element-type 'dkt-csth' group 'Plate' slot 'Str' ... embedded range group 'FC_PlateTop' slot 'Int' struct liner property density @RoSt thick @ThPl isotropic @ESt @NuSt ;Liner interface properties of side 1 and 2 struct liner property ... coupling-stiffness-normal @ksLiner coupling-stiffness-normal-2 @ksLiner ... coupling-stiffness-shear @ksLiner coupling-stiffness-shear-2 @ksLiner ... coupling-cohesion-shear @cLiner coupling-cohesion-shear-2 @cLiner ... coupling-friction-shear @fiSt_bot coupling-friction-shear-2 @fiSt_bot ... coupling-yield-normal 0.0 coupling-yield-normal-2 0.0 struct node group 'ND_PlateLink' slot 'Int' range group 'Plate' slot 'Str' ... sphere center (@xyRod,@xyRod,@zPl) radius 0.05 ; Create a single cable and set its associated properties struct cable create by-line (@xyRod,@xyRod,@zPl) (@xyRod,@xyRod,0.1) ... group 'Rod' slot 'Str' segments @segRod struct cable property density @RoSt cross-sectional-area @ARod young @ESt ... yield-compression [YSt*ARod] yield-tension [YSt*ARod] ... grout-perimeter @pgCable grout-stiffness @kgCable ... grout-cohesion 0.0 grout-friction @fiSt_top range group 'Rod' slot 'Str' struct cable property grout-friction @fiSt_bot range position-z @zwt @zPl ... group 'Rod' slot 'Str' 192 struct node group 'ND_RodLink' slot 'Int' range group 'Rod' slot 'Str' ... sphere center (@xyRod,@xyRod,@zPl) radius 0.01 struct node group 'ND_RodHead' slot 'Int' range group 'Rod' slot 'Str' ... sphere center (@xyRod,@xyRod,0.1) radius 0.01 ; Create a link to connect the liner to cable struct link tolerance-node 0.05 struct link create group 'LK_PlateRod' slot 'Int' side 2 target node ... range group 'ND_RodLink' slot 'Int' struct link attach x rigid range group 'LK_PlateRod' slot 'Int' ; Applying symmetry BC on the liner elements struct node group 'ND_PlateXside' slot 'BC' range ... position-x 0.05 [rPl+0.01] position-y -0.05 0.05 struct node group 'ND_PlateYside' slot 'BC' range ... position-y 0.05 [rPl+0.01] position-x -0.05 0.05 struct node group 'ND_PlateCenter' slot 'BC' range ... position-x -0.05 0.05 position-y -0.05 0.05 model cycle 0 struct node fix system-local range group 'Plate' ;---------------------------------------struct node fix velocity-y range group 'ND_PlateXside' slot 'BC' struct node fix rotation-x rotation-z range group 'ND_PlateXside' slot 'BC' struct node fix velocity-x range group 'ND_PlateYside' slot 'BC' struct node fix rotation-y rotation-z range group 'ND_PlateYside' slot 'BC' struct node fix velocity-x velocity-y range group 'ND_PlateCenter' slot 'BC' struct node fix rotation-x rotation-y rotation-z range group 'ND_PlateCenter' slot 'BC' ;---------------------------------------model solve ratio 1e-5 model save ['Plate_Rod_RAP_const_'+RAP_Const+'_'+Constit_model] ; ; ; ; ; ; ; ========================================================================== ========================================================================== Subroutine name: " LoadTest.dat" Programing language: FLAC3D commands and FISH Description: This subroutine simulates the pull-out test ========================================================================== ========================================================================== model restore ['Plate_Rod_RAP_const_'+RAP_Const+'_'+Constit_model] ;---------------------------------------zone gridpoint initialize displacement (0,0,0) zone gridpoint initialize velocity (0,0,0) structure node initialize displacement (0,0,0) structure node initialize velocity (0,0,0) ; Getting history of force and displacement history interval 50 struct node history name 'RodHead' displacement-z position (@xyRod,@xyRod,0.1) struct node history name 'RodBase' displacement-z position (@xyRod,@xyRod,@zPl) struct node history name 'PlateCenter' displacement-z position (0.0,0.0,@zPl) struct cable history name 'RodForce' force position (@xyRod,@xyRod,0.1) ; Fish function for force-control pull out test and recording data fish define LoadRAP(deltaload) local rh local rb local plc local nd local gp local ii = 1 local jj = 0 local kk = 0 local idx local ndrh 193 local ndrb local ndplc local ndzdrh local ndzdrb local ndzdplc local Frh local Frb local Fplc global TopLoad local NameFile local MaxTopDisp local MaxDeltaStep local OldStep ;-------------ndrh = struct.node.near(xyRod,xyRod,1.0) ;Node at the top of the rod ndrb = struct.node.near(xyRod,xyRod,zPl) ;Node at the bot. of the rod ndplc = struct.node.near(0.0,0.0,zPl) ;Node at the center of plate rh = struct.near(xyRod,xyRod,0.1) ;Elem. at the top of the rod rb = struct.near(xyRod,xyRod,zPl+(0.1-zPl)/segRod/2);Elem. at the bot. of the rod plc = struct.near(0.0,0.0,zPl) ;Elem. at the center of plate Frh = struct.cable.force.nodal(rh,2,3) ;Rod axial force at the top node Frb = -1*struct.cable.force.nodal(rb,1,3);Rod axial force at the bot. node Fplc = struct.shell.force.nodal(plc,3,3) ;Liner normal force at the center ;-------------ndzdrh = math.abs(struct.node.disp.local(ndrh,1)) ;Rod top node disp ndzdrb = math.abs(struct.node.disp.local(ndrb,1)) ;Rod bot. node disp ndzdplc = math.abs(struct.node.disp.local(ndplc,3));Plate center node disp ;x: force, y: displacement table.x('LdCurRH',1) = Frh ;Top of the rod table.y('LdCurRH',1) = ndzdrh table.x('LdCurRB',1) = Frb ;Bot. of the rod table.y('LdCurRB',1) = ndzdrb table.x('LdCurPlC',1) = Fplc ;Center of the plate table.y('LdCurPlC',1) = ndzdplc ;-----------------------------------------------------; Set criteria for loop elimination MaxTopDisp = 7.0/12 MaxDeltaStep = 50000 ;----------------------loop while ndzdrh <= MaxTopDisp TopLoad = float(ii)*deltaload jj = ii+1 OldStep = global.step ;------------------idx = int((ii-kk)/5) if idx = 1 then NameFile = 'Load_'+string(TopLoad)+'_'+RAP_Const+'_'+Constit_model kk = ii end_if command model save ['Load_Final_Stable_Load_'+RAP_Const+'_'+Constit_model] ;--------;K tables @K_hist table 'K15' export ... ['K15_Uplift_Final_Stable_Load_'+RAP_Const+'_'+Constit_model] truncate csv table 'K25' export ... ['K25_Uplift_Final_Stable_Load_'+RAP_Const+'_'+Constit_model] truncate csv table 'K50' export ... ['K50_Uplift_Final_Stable_Load_'+RAP_Const+'_'+Constit_model] truncate csv ;Radial displacement 194 table 'Rad_disp' delete @Rad_disp_RAP_table('Rad_disp',[rPl*sin45],[rPl*sin45],1.0) table 'Rad_disp' export ... ['Disp_r15_th45_Final_Stable_Load_'+RAP_Const+'_'+Constit_model] truncate csv table 'Rad_disp' delete @Rad_disp_RAP_table('Rad_disp',[rPl],0.0,1.0) table 'Rad_disp' export ... ['Disp_r15_th00_Final_Stable_Load_'+RAP_Const+'_'+Constit_model] truncate csv ; Applying load structure node apply force (0.0,0.0,@TopLoad) range ... group 'ND_RodHead' slot 'Int' model cycle 500 model solve ratio 1e-5 cycles-total [OldStep+MaxDeltaStep] structure node apply force (0.0,0.0,0.0) range ... group 'ND_RodHead' slot 'Int' end_command ;----------------ndzdrh = math.abs(struct.node.disp.local(ndrh,1)) ndzdrb = math.abs(struct.node.disp.local(ndrb,1)) ndzdplc = math.abs(struct.node.disp.local(ndplc,3)) Frh = struct.cable.force.nodal(rh,2,3) Frb = -1*struct.cable.force.nodal(rb,1,3) Fplc = struct.shell.force.nodal(plc,3,3) table.x('LdCurRH',jj) = Frh table.y('LdCurRH',jj) = ndzdrh;compare to [TopLoad] table.x('LdCurRB',jj) = Frb table.y('LdCurRB',jj) = ndzdrb table.x('LdCurPlC',jj) = Fplc table.y('LdCurPlC',jj) = ndzdplc ;----------------if idx = 1 then command @K_hist table 'K15' export ['K15_Uplift_'+NameFile] truncate csv table 'K25' export ['K25_Uplift_'+NameFile] truncate csv table 'K50' export ['K50_Uplift_'+NameFile] truncate csv ;Radial displacement table 'Rad_disp' delete @Rad_disp_RAP_table('Rad_disp',[rPl*sin45],[rPl*sin45],1.0) table 'Rad_disp' export ['Disp_r15_th45_'+NameFile] truncate csv table 'Rad_disp' delete @Rad_disp_RAP_table('Rad_disp',[rPl],0.0,1.0) table 'Rad_disp' export ['Disp_r15_th00_'+NameFile] truncate csv ;--------model save [NameFile+'.sav'] end_command end_if ii = ii+1 if global.step >= (OldStep+MaxDeltaStep) then exit end_if end_loop end @LoadRAP(600) ;------------------------------------table 'LdCurRH' export ['Load_disp_rod_head_'+RAP_Const+'_'+Constit_model] truncate csv table 'LdCurRB' export ['Load_disp_rod_base_'+RAP_Const+'_'+Constit_model] truncate csv table 'LdCurPlC' export ['Load_disp_plate_center_'+RAP_Const+'_'+Constit_model] ... truncate csv ;------------------------------------model save ['Load_Final_'+RAP_Const+'_'+Constit_model] REFERENCES Bowles, J. E. (1996). Foundation Analysis and Design, Fifth Edition. The McGraw-Hill Compenies, Inc., New York. Cheng, Z., and Detournay, C. (2016). “Plastic hardening model I : Implementation in FLAC3D.” 4th Itasca Symposium on Applied Numerical Modeling, Lima, Peru, 267– 276. Cheng, Z., and Lucarelli, A. (2016). “Plastic hardening model II : Calibration and validation.” 4th Itasca Symposium on Applied Numerical Modeling, Lima, Peru, 393– 402. Das, B. M. (1987). Theoretical Foundation Engineering. Elsevier Science Publishing Company Inc., New York, NY. Demir, S., Özener, P., and Kirkit, M. (2017). “Experimental and numerical investigations of behavior of rammed aggregate piers.” Geotechnical Testing Journal, ASTM International, 40(3), 411–425. Duncan, J. M., Byrne, P., Wong, K. S., and Mabry, P. (1980). Strength, Stress-Strain and Bulk Modulus Parameters for Finite Element Analysis of Stresses and Movements in Soil Masses. Report No. UCB/GT/80-01, Dept. of Civil Engineering, U.C. Berkeley, Berkeley, California. Fox, N. S., and Cowell, M. J. (1998). Geopier Foundation and Soil Reinforcement Manual. Geopier Foundation Company, Inc, Scottsdale, Arizona. Halabian, A. M., and Shamsabadi, P. J. (2014). “Numerical Modeling of the RAP Construction Process and Its Effects on RAP Behavior.” International Journal of Geomechanics, 15(5), 04014085. Handy, R. L., and White, D. J. (2005a). “Stress Zones Near Displacement Piers: I. Plastic and Liquefied Behavior.” Journal of Geotechnical and Geoenvironmental Engineering, 132(1), 54–62. Handy, R. L., and White, D. J. (2005b). “Stress Zones near Displacement Piers: II. Radial Cracking and Wedging.” Journal of Geotechnical and Geoenvironmental Engineering, 132(1), 63–71. Hsu, C. L. (2000). “Uplift Capacity of Geopier Foundations.” Master of Science thesis, University of Utah. 196 Itasca Consulting Group, I. (2019). “FLAC3D – Fast Lagrangian Analysis of Continua in 3 Dimensions, Version 7.0.” Minneapolis, Minnesota. Kwong, H. K., Lien, B., and Fox, N. S. (2002). “Stabilizing Landslides Using Rammed Aggregate Piers.” Fifth Malaysian Road Conference, Kuala Lumpur, Malaysia, 7. LaMeres, B. J. (2005). “Numerical Analysis to Determine the Ultimate Uplift and Bearing Capacities of Single Piers.” Master of Science thesis, University of Utah. Lawton, E. C. (2000). Performance of Geopier Foundations during Simulated Seismic Tests at South Temple Bridge on Interstate 15, Salt Lake City Utah. Rep. No. UUCVEEN 00-03, Dept. of Civil and Environmental Engineering, Univ. of Utah, Salt Lake City, Utah. Lawton, E. C., and Fox, N. S. (1994). “Settlement of structures supported on marginal or inadequate soils stiffened with short aggregate piers.” Geotechnical Special Publication, 962–974. Lawton, E. C., Fox, N. S., and Handy, R. L. (1994). “Control of Settlement and Uplift of Structures Using Short Aggregate Piers.” Geotechnical Special Publication, 121–132. Lucarelli, A., and Cheng, Z. (2016). “Plastic hardening model III : Design application.” 4th Itasca Symposium on Applied Numerical Modeling, Lima, Peru, 403–414. Mayne, P. W. (2007). Cone Penetration Testing - A Synthesis of Highway Practice. NCHRP Synthesis 368, Transportation Research Board, Washington, DC. Mitchell, J. K., and Villet, W. C. B. (1987). Reinforcement of Earth Slope and Embankment. NCHRP Report 290, Transportation Research Board, Washington D.C. Ozer, A. T. (2005). “Estimation of Consolidation and Drainage Properties for Lake Bonneville Clays.” PhD dissertation, University of Utah. Ozer, A. T., Bartlett, S. F., and Lawton, E. C. (2013). “CPTU and DMT for estimating soil unit weight of Lake Bonneville clay.” 4th International Conference on Geotechnical and Geophysical Site Characterization, Porto de Galinhas, Brazil, 291–296. Pham, H. T. V., and White, D. J. (2007). “Support Mechanisms of Rammed Aggregate Piers. II: Numerical Analyses.” Journal of Geotechnical and Geoenvironmental Engineering, 133(12), 1512–1521. Potyondy, J. G. (1961). “Skin friction between various soils and construction materials.” Geotechnique, 11(4), 339–353. Rudolph, R. W., Serna, B., and Farrell, T. (2011). “Mitigation of Liquefaction Potential Using Rammed Aggregate Piers.” Geo-Frontiers 2011, American Society of Civil Engineers, Reston, VA, 557–566. 197 Schanz, T., Vermeer, P. A., and Bonnier, P. G. (1999). “The Hardening Soil Model: Formulation and Verification.” Beyond 2000 in Computational Geotechnics.Ten Years of PLAXIS International.Proceedings of the International Symposium, Balkema, Amsterdam, Netherlands, 281–296. White, D. J., Pham, H. T. V., and Hoevelkamp, K. K. (2007). “Support Mechanisms of Rammed Aggregate Piers. I: Experimental Results.” Journal of Geotechnical and Geoenvironmental Engineering, 133(12), 1503–1511. White, D., Wissmann, K., and Lawton, E. (2001). “Geopier soil reinforcement for transportation applications.” Geotechnical News, 19(4), 63–68. Wissmann, K. J., and FitzPatrick, B. T. (2016). Geopier Uplift Resistance. Technical Bulletin No. 3, Geopier Foundation Company, Inc., Davidson, North Carolina.