UFDC Home  myUFDC Home  Help 



Full Text  
EVALUATION OF THE DISPLACEMENT DISCONTINUITY METHOD WITH TESSELLATIONS TO SIMULATE THE INDIRECT TENSION STRENGTH TEST By CHOTE SORANAKOM A THESIS PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF ENGINEER UNIVERSITY OF FLORIDA 2003 Copyright 2003 by Chote Soranakom ACKNOWLEDGMENTS Many wishes were sent with a very special acknowledgment to my advisor, Dr. Bjorn Birgisson, Assistant Professor of Civil and Coastal Engineering, who supervised the project. His enthusiasm and research creativity in the field of pavement made the project possible and successful. I would also like to thank Dr. Reynaldo Roque for his assistance in providing test data and useful recommendation for my research. Special thanks go to Dr. John A.L. Napier, who let us use his boundary element code and provided technical assistance throughout the research. I also would like to thank Dr. Bhavani V. Sankar for his kindness, teaching and participation as a member of my committee. I would like to take this opportunity to acknowledge the Florida Department of Transportation for providing us financial support. Also, I would like to thank many of my friends, Boonchai Sangpetngam, Daniel Darku, Oscar Tazoe, Chirstos Drakos, whose advice, assistance, support and friendship helped me throughout the project. Finally, a very special acknowledgement is sent to my parents, Chewpore and Wongsa Soranakom, for their love, encouragement and constant advice throughout my academic life away from home. I would also like to thank my brothers and sister in law Charn, Chuan and Bongkoch for their love and emotional support throughout the good and difficult times of doing research. TABLE OF CONTENTS page A C K N O W L E D G M E N T S ...................... .. ..................................................................... iii LIST OF TABLES ............. ................... ........... .. .. ............... .. vii LIST OF FIGURES ......... ......................... ...... ........ ............ ix A B S T R A C T .............................................. ..........................................x iii CHAPTER 1 IN TR OD U CTION ............................................... .. ......................... .. 1.1 B background ......... ...... ................................................................... ........... 1 1.2 Problem Statement ....................... ............ ... ............... 2 1.3 Research Hypothesis............................................... ...... ......... 4 1.4 Objectives ......................... ............ ...................4. 4 1.5 Scope..................................................... . 4 1.6 R research A approach ....... ........................... ............ .. .... .. ...... .. .. 2 LITER A TU R E REV IEW ............................................................. ....................... 6 2.1 O overview ................ ......................................................................... ............. 2.2 C classical F atigue A approach ............................................................. .. ............... 2.3 Continuum D am age Approach ........................................ .......................... 9 2.4 Fracture M echanics Approach .................................... ... ................ ........... 11 2.3.1 Conventional Fracture M echanics ......................................... .................12 2.3.2 Application of Fracture M echanics ....................................... .................16 2.5 Displacement Discontinuity Method with Tessellation Schemes ......................18 2.6 Sum m ary ............... .................................. ..........................2 1 3 DISPLACEMENT DISCONTINUITY METHOD WITH TESSELLATION SCH EM E ............................................................... ..... ...... ........ 22 3.1 Overview ................................ .... ............ ...... ....... ..... ........... 22 3.2 D isplacem ent D iscontinuity M ethod ........................................ .....................22 3.3 N um erical Im plem entation ............................................................................ ...25 3.4 Example of a Crack Growth Model using DIGS...............................................27 4 NUMERICAL EXPERIMENTS OF THE SUPERPAVE INDIRECT TENSION TEST USING THE DISPLACEMENT DISCONTINUITY METHOD WITH A TESSELLATION SCHEM E ........................................................ ............... 31 4.1 O overview ................. .. ............................................................. 31 4.2 Laboratory Experiment of Nova Scotia Mixture...............................................32 4 .3 N um erical Sim ulations ..................... .......... ....... ......... ..............................34 4.4 Determination of Suitable Tessellation Schemes and Crack Growth Rules.........35 4 .4 .1 N um erical E xperim ent..................................................... .....................36 4.4.2 Results from Numerical Simulations............................................. 38 4.4.2.1 C rack G row th R ules............................................ .....................38 4.4.2.2 Evaluation of Delaunay Tessellation Scheme.............................44 4.4.2.3 Evaluation of Voronoi Tessellation Scheme.............................44 4.4.2.4 Evaluation of Voronoi with Internal Fracture Paths Tessellation Schem e ................... ...... ............ .......................................... 45 4.4.3 Summary of Findings on the Determination of a Suitable Crack Growth Rule and Tessellation Schem e ................................. ......................... 46 4.5 Determination of Average Voronoi Aggregate Size to Represent Gradation ......46 4.5.1 N um erical Experim ent.......................................... ........... ............... 46 4.5.2 N um erical R results ............... ............... .. ........................ ............... 48 4.6 Sensitivity Analysis of Material Parameters Defined for Mastics...................52 4.6.1 N um erical Experim ent.......................................... ........... ............... 52 4.6.2 N um erical R results ............................................. ...... ....................... 53 4.7 Sensitivity Analysis of Material Parameters Defined for Internal Fracture Paths ............... .............................. ................................ 61 4.7.1 N um erical Experim ent.......................................... ........... ............... 61 4 .7 .2 N u m erical R esu lts ........................................................... .....................6 1 4.8 Sum m ary of Sensitivity A analysis .................. ...................................... .... 62 4.9 Recommendations for a Successful Numerical Simulation of IDT Strength T e s t .................................................................................................................. 6 7 4.10 Summary ......... ...... ........ .. ............... .........69 5 EVALUATION OF TENSILE STRENGTH AND FRACTURE ENERGY D E N S IT Y .............................................................................7 1 5.1 Overview .................. ........................................... 71 5.2 Numerical Model V.S. Asphalt Mixtures......................... .. ...............73 5.3 Numerical Simulation of Three Mixtures............... ........... .................75 5.4 N um erical Test Results........................ ................... ................. ............... 77 5 .5 C o n clu sio n s.................................................. ................ 8 3 6 SUMMARY, CONCLUSIONS AND RECOMMENDATIONS ...........................84 6 .2 C o n clu sio n s................................................. .................. 8 5 6 .3 R ecom m en nation s.................................................................................... .. 85 APPENDIX USER MANUALS FOR PRE/POST PROCESSOR (DDM)...................87 A 1 Introduction ............................................ 87 A .2 C coordinate System s ...................................................... ............ ............... 88 A.3 Text Files Associated with PRE/POST Processor ............... ....... ............89 A .3.1 Input File (*.IN File) ........................................ ........................... 89 A .3.2 Segm ent File (*._SG ) ........................................ .......................... 94 A .3.3 D ata File (*.D A T )........................................................... ............... 97 A .3.4 O utput F ile (*.O U T ) .................................................................... ..... 106 A .3.5 R request File (*.RE Q ) ..................................... ....................... ............ 108 A .3.6 R report F ile (*.R P T ) .................................... ................... ... ............ 111 A.4 The PRE/POST Processor (DDM) .......................... .... ............................112 A.4.1 Symbols and Terminologies .......................................... ...............112 A .4.2 M enu B ar ................................. ................... ........ .. ......... 113 L IST O F R E FE R E N C E S ............................ ................................................... ............ 136 BIOGRAPHICAL SKETCH ........................................................... ............... 142 LIST OF TABLES Table page 41. Gradation of the N ova Scotia m ixture ............. .................... ............. ............... 33 42. M material properties from IDT strength test ................................... .................34 43. Modeling matrix for the study of suitable tessellation schemes and crack growth ru le s ................................................................................3 7 44. Global linear elasticity parameters for the study of suitable tessellation schemes and crack grow th rules ....................... .. ...................... .... ...... .... ...... ...... 37 45. Local material parameters for the study of suitable tessellation schemes and crack grow th rules ...................... .... ................ ................... .... ....... 37 46. Modeling matrix for the study of size effect: ............................. ............... 47 47. Global linear elasticity parameters for the study of size effect.............................48 48. Local material parameters for the study of size effect ............................................48 49. Modeling matrix for the parametric study of mastic.................... .................52 410. Global linear elasticity parameters for the parametric study of mastic ..................52 411. Local material parameters (calibrated material parameters) for the parametric study of m astic .......................................................................................... ........... 53 412. Varying parameters from the calibrated mastic shown in Table 411 for the parametric study of mastic .......................... ........................... ...............53 413. Modeling matrix for the parametric study of internal fracture paths ..................63 414. Global linear elasticity parameters for the parametric study of internal fracture p a th s ............................................................................. 6 3 415. Local material parameters (calibrated material parameters) for the parametric study of internal fracture paths ..................... ................................................ 63 416. Varying parameters from the calibrated fracture path shown in Table 415 for the parametric study of internal fracture paths.....................................................63 51. Correction factors accounting for bulging effect ............................................. 72 52. Correction factors accounting for stress at center of specimen.............................73 53. Gradation of three mixtures used in simulations.............. .... .................74 54. Material properties of three mixtures obtained from Superpave IDT tests..............75 55. Modeling matrix for three mixtures used:............................. .................76 56. Global linear elasticity parameters for three mixtures used..............................76 57. Local material parameters for three mixtures used ..............................................77 LIST OF FIGURES Figure page 31. Displacement discontinuity element in local coordinate (yz)..............................23 32. Tessellations for granular structures. ............................................ ............... 29 33. Mohr Coulomb type of failure criterion for determining crack mobilization. .........29 41. Indirect tension strength test ............................................................................. 34 42. Predicted number of cracks at each load step........... .... .................39 43. Comparison of predicted and measured horizontal and vertical stressstrain curves, using D elaunay tessellation .............................................. ............... 40 44. Comparison of predicted and measured horizontal and vertical stressstrain curves, using V oronoi tessellation ........................................ ....................... 41 45. Comparison of predicted and measured horizontal and vertical stressstrain curves, using Voronoi with internal fracture path tessellation.............................. 42 46. Crack patterns predicted by two crack growth rules of three tessellations ..............43 47. Comparison of predicted and measured horizontal and vertical stressstrain curves for five aggregate sizes (4.5, 6.0, 8.0, 10.0 and 12.0 mm.), using Voronoi with internal fracture path tessellation ..........................................50 48. Ultimate tensile strength averaged from 3 samples for each aggregate size (4.5, 6.0, 8.0, 10.0 and 12.0 m m ). ........................................ ........................ 51 49. Variance of residual error of tensile stress strain curve for each aggregate size (4.5, 6.0, 8.0, 10.0 and 12.0 m m ). ........................................ ........................ 51 410. Change of stressstrain responses due to the varying of cohesion defined in mastic from 2.00 2.80 mm ....................................... ............... 55 411. Sensitivity of the parameter cohesion defined in mastic to the ultim ate tensile strength............................................. ............................. 56 412. Sensitivity of the parameter tension cutoff defined in mastic to the ultim ate tensile strength........................................................... ............... 56 413. Sensitivity of the parameter friction angle defined in mastic to the ultim ate tensile strength............................................. ............................. 57 414. Sensitivity of the parameter cohesion softening slope defined in mastics to the ultim ate tensile strength........................................................... ............... 57 415. Sensitivity of the parameter tension softening slope (Linear Model) defined in m astic to the ultim ate tensile strength. ......................................... ...............58 416. Sensitivity of the parameter opening crack limit defined in mastic to the ultim ate tensile strength........................................................... ............... 58 417. Sensitivity of the parameter residual cohesion defined in mastic to the ultim ate tensile strength............................................. ............................. 59 418. Sensitivity of the parameter residual friction angle defined in mastic to the ultim ate tensile strength........................................................... ............... 59 419. A comparison of stressstrain response between using a Linear (default) and Nonlinear (still in experiment) tension softening model for the range of opening crack limit DNCR 0.03 0.19 mm .................................... ............... 60 420. Sensitivity of both parameters cohesion and tension cutoff defined in fracture path of aggregate to the ultimate tensile strength. ......................................64 421. Sensitivity of opening crack limit defined in internal fracture path of aggregate to the ultimate tensile strength. ...................................... ............... 64 422. Sensitivity of friction angle defined in internal fracture path of aggregate to the ultim ate tensile strength........................................................... ............... 65 423. Sensitivity of cohesion softening slope defined in internal fracture path of aggregate to the ultimate tensile strength. ...................................... ............... 65 424. Sensitivity of residual cohesion defined in internal fracture path of aggregate to the ultimate tensile strength. ...................................... ............... 66 425. Sensitivity of residual friction angle defined in internal fracture path of aggregate to the ultimate tensile strength. ...................................... ............... 66 51. Comparison of predicted and measured horizontal and vertical stressstrain curves for three mixtures: Nova Scotia, NW39 1C and 195SJN BWP ..................79 52. Simulated deformation differentials for the Nova Scotia, NW39 1C and I95SJN BW P m ixtures. ..... ........................... .........................................80 53. Tensile strength at fracture for three mixtures (Nova Scotia, NW39_1C and 195 SJN B W P) ................................................................... ................. 80 54. Fracture energy density at center of specimen for three mixtures (Nova Scotia, NW39_1C and 195SJNBWP) ......................................................81 55. Predicted crack patterns for the three mixtures (Nova Scotia, NW39_1C and 195SJN _BW P) at three load steps ................................................. ....... ........ 82 A1. Global (YZ) and local (yz) coordinate systems used by DIGS ............................88 A2. Square plate with a circular hole for the input file ................................................93 A3. Square plate with a circular hole for the segment files ........................................95 A 4. M ohr Coulom b type of failure. ............. ........................................ .....................99 A5. Square plate with a circular hole for the data file............................................104 A6. The DDM m enu bar. ........................ ..... ................. ............... 113 A7. Select elem ent group dialog box. ......................................................... .... ......... 116 A 8. M odel setting dialog box .......................................................................... ....... 119 A 9. Tessellation dialog box................................................. .............................. 122 A 10.Export segm ent dialog box ...................................................................... 123 A 11.Line segm ent dialog box. ....................................................................... 124 A 12. Circular segm ent dialog box. ............................................................................ 125 A 13. L inear field point dialog box ................................................................ ........... 126 A 14.M atrix field point dialog box. ............................................................. ............126 A15. Boundary condition dialog box. ........................................ ......................... 128 A16. General param eter dialog box. ........................................ ......................... 129 A17. Cemented material dialog box................................................... 130 A18.Load step dialog box. .............................. ...... .................. 131 A19. Include segment file dialog box. ................................. ............... 132 A 20. A m plification dialog box. ............................................. ............................. 133 A21. Color for DD element dialog box ........................................................... 133 A 22. Probe option dialog box ......... ......... ...................... .................. ............... 134 A23.Request output file dialog box. ........ ........... ..... .......... .. ....... ........ 135 Abstract of Thesis Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Engineer EVALUATION OF THE DISPLACEMENT DISCONTINUITY METHOD WITH TESSELLATIONS TO SIMULATE THE INDIRECT TENSION STRENGTH TEST By Chote Soranakom May, 2003 Chair: Bjom Birgisson Cochair: Reynaldo Roque Major Department: Civil and Coastal Engineering The research focused on the use of the displacement discontinuity method to simulate cracking behavior of the indirect tension strength test of asphalt mixtures. From numerical experiments, it was found that the granular structure of asphalt mixtures could be best represented by Voronoi tessellations with internal fracture paths. Two distinct crack growth rules, SEQUENTIAL and PARALLEL, were evaluated and no difference was found in terms of the predicted crack patterns and stress strain responses. However, the PARALLEL crack growth rule that activates all crack sites when the stress exceeds failure limit is more efficient than the SEQUENTIAL crack growth rule that activates the most overstressed crack site one at a time. The study of the aggregate size effects revealed a small size effect. An increase in average aggregate diameter from 4.5 mm to 12.0 mm increased the ultimate tensile strength of the specimen by 11%. It was also found that the optimal particle size was around 40 to 60% passing of gradation. The sensitivity analysis of material parameters revealed that the mastic properties were a key factor controlling the cracking behavior of mixtures. Finally the method was used to simulate crack propagation of three mixtures and evaluate their tensile strength and fracture energy density. The numerical results showed reasonable predictions of crack patterns and stress strain responses. The predicted tensile strength and fracture energy density were acceptable. CHAPTER 1 INTRODUCTION 1.1 Background It has long been accepted that cracking of hotmix asphalt (HMA) pavements is a major mode of premature failure. A solid understanding of the mechanisms of crack initiation and crack growth is essential to predict pavement performance in the context of thickness design, as well as in the design and optimization of mixtures. Any hope of developing more crack resistant mixes relies on improved understanding of the mechanics of crack initiation and crack growth in HMA mixtures. Unfortunately, the complexity of modeling crack initiation and crack growth has been an obstacle to the incorporation of fracture mechanicsbased approaches in the bituminous pavement area. With an increasing of computer power, numerical analysis of cracking becomes possible. Bazant (1986) provided a good review for cracking models that have been used for granular material such as concrete and rocks. Cracking in material is normally analyzed by either a fracture mechanics approach or smeared crack approach. The former assumes a single line of crack initiates from a preexisting flaw and propagates through material according to a certain crack growth criterion such as maximum energy release rate. In contrast, the latter assumes cracks are smeary distributed over a finite region such that an average tensile strain can adequately represent crack. Even though fracture mechanics and smeared crack approaches can predict cracking behavior, they do not really capture the characteristic of cracks in granular material in which cracks randomly initiate at critical locations and coalesce, forming major crack bands. Explicit fracture modeling using random assemblies of displacement discontinuity boundary elements provides a more realistic approach for analyzing discrete cracks in granular material (Napier, 1990; Napier and Hildyard, 1992; Napier and Pierce, 1995(a b); Malan and Napier, 1995; Kuijpers and Napier, 1996; Napier et al., 1997; Napier and Malan, 1997, Steen et al., 2001, Birgisson et al., 2002(ab) and Soranakom et al., 2003). The method takes advantage of the known solutions of displacement discontinuity in an infinite medium to formulate system of equation; thus, the method is efficient and yields good results for cracktype problems. The change of geometry due to crack propagation in material can easily be handled by placing displacement discontinuity elements along the anticipated crack paths, which are normally assumed at the boundaries of the grain particles and perhaps internal to the grains as well. Another advantage of using the boundary element approach is that the dimension of the problem decreases by one. That means a twodimensional problem needs onedimensional element with a simple failure law; therefore, the computational cost is greatly reduced. For the advantages mentioned above, this method provides an elegant way to simulate and analyze cracks in asphalt mixtures. 1.2 Problem Statement The displacement discontinuity boundary element method has successfully been used in the field of rock mechanics and mining engineering to model the cracking behavior in rock mass and/or excavations underground. Based on this work, it was determined that the method had the potential to be used in the field of pavements to study the effects of cracking on asphalt mixtures. Since no previous work has been performed on the use and implementation of this method to simulate and evaluate the cracking behavior of asphalt mixtures, several questions regarding the appropriate set of numerical components arise, which are summarized below: * What is the best tessellation schemes to simulate the cracking behavior of asphalt mixtures? This includes the type of tessellation scheme (Delaunay triangles, Voronoi polygons or Vornonoi with internal fracture paths) and the optimal particle size used in tessellations to represent aggregate structure of IDT specimen. * Identification of appropriate crack growth rules for using in simulation of cracking behavior of mixtures. This includes an evaluation of appropriate numerical schemes for predicting the crack growth. Currently, two basic schemes have been evaluated, namely the SEQUENTIAL and PARALLEL crack growth rules. In the PARALLEL crack growth rule, all elements that reach failure in any given loading step are allowed to crack. In the SEQUENTIAL crack growth law, only the most overstressed element is allowed to crack, and the resulting redistribution of stresses is performed and the process is repeated until all overstressed elements at any given load step have been allowed to crack. * Evaluation of the current failure law implemented in the boundary element code DIGS. In the displacement discontinuity boundary element method, all nonlinear behavior is lumped into the displacement discontinuity. The rest of the material is considered linear elastic. This allows for the formulation of various local "failure laws" to simulate the fracture behavior of materials. For example, a typical assumption would be that the local failure in a sliding mode at a given displacement discontinuity is governed by a MohrCoulomb failure law, such that once the stress state reaches the failure envelope, sliding is initiated. Similarly, in tension, a typical failure law would be that the material would open in tension once a given local tensile strength criteria is reached. The nature of these failure laws can be very flexible, allowing for peak and residual strength criteria, with degradation laws describing the transition from the peak strength to the residual strength. * It is also important to determine consistent means of obtaining key material parameters for mixtures. This includes the identification of parameters that can be obtained directly from existing laboratory testing procedures, as well as from calibration. In particular, any parameters that need to be obtained from calibration need to be scrutinized in terms of the consistency of the calibration approach used from one mixture to another. * The potential of the method needs to be evaluated in terms of how well it captures key mixture characteristics and trends observed under carefully controlled laboratory testing conditions. In particular, to what extent this method can predict cracking behavior of asphalt mixtures in the laboratory is of interest. This includes determining how well the method can simulate o The pre and postpeak stressstrain behavior under tensile loading conditions o The initiation of fracture o The fracture energy and tensile strength at the point of fracture 1.3 Research Hypothesis The displacement discontinuity boundary element method with an appropriate tessellation scheme can be used to simulate the observed cracking behavior of asphalt mixtures in the laboratory under controlled loading and environmental conditions. Therefore, it is possible to identify an appropriate set of model conditions to describe the fracture behavior of HMA from a micromechanical point of view. 1.4 Objectives The objectives of this research are summarized below: * Identify the appropriate tessellation schemes and optimal particle size to simulate aggregate structure of asphalt mixtures. * Identify the appropriate crack growth rule for simulating indirect tension strength test under monotonic loading conditions. * Evaluate the performance of the current MohrCoulomb type of failure to model cracking behavior of mixture. * Establish numerical calibration procedure to obtain material parameters from a given mixture. * Evaluate the potential of using displacement discontinuity method with tessellation to determine tensile strength and fracture energy density of mixtures 1.5 Scope This study focuses on the implementation of the displacement discontinuity method to the field of pavements. The project will start with a study to identify the appropriate tessellation schemes and crack growth rules. Then, the study will focus on the evaluating the sensitivity of stressstrain predictions to the various input parameters used. These include material parameters for the mastic, material parameters for fracture path inside aggregate, and the determination of the optimal average size of aggregates. The sensitivity of these parameters will be used to establish a standard numerical procedure to simulate indirection strength test. Finally, the numerical model will be used to simulate and evaluate how well the method can predict mixture properties such as tensile strength and fracture energy density. 1.6 Research Approach The research was conducted in five phases: * Learning how to use the boundary element code entitled "Discontinuity Interaction and Growth Simulation" (DIGS). * Writing preand post processors to help input data and interpret results from DIGS * Conducting the following numerical studies: o Compare tessellations schemes (Delaunay, Voronoi and Voronoi with internal fracture paths). o Compare crack growth laws (SEQUENTIAL and PARALLEL). o Evaluate the current failure law, which is a Mohr Coulomb type of failure law by means of evaluating the sensitivity of material parameters used on the resulting predictions of the stressstrain behavior of mixtures. o Evaluate any potential aggregate size effects. * Establishing a numerical procedure to calibrate material properties from a given mix for use in a numerical model to simulate indirect tension strength test. * Using numerical models to evaluate tensile strength and fracture energy of three mixtures: Nova Scotia, NW39 IC and 195SJN BW. CHAPTER 2 LITERATURE REVIEW 2.1 Overview The purpose of this chapter was to conduct a literature review of commonly used approaches in the analysis of cracking and damage in asphalt mixtures. Another popular analytical tool, displacement discontinuity method, which has been used to analyze cracking in field of rock mechanics, was also reviewed. The four approaches that the literature review focuses on are: * Classical fatigue approach * Continuum damage approach * Fracture mechanics approach * Displacement discontinuity approach with tessellation schemes 2.2 Classical Fatigue Approach Fatigue cracking due to repeated traffic loading is one of the major distresses in asphalt concrete pavements. The mechanism of fatigue cracking is caused by the tensile stress at either the top or the bottom of asphalt layer induced by wheel load. With several number of load passes, microcracks will develop and later coalesce to form major crack bands, weakening the pavement structure and finally leading to pavement failure. In order to characterize the fatigue resistance of a mixture for pavement design, it was necessary to describe the behavior of AC mixtures under repeated stress or strain cyclic loading. In early developments, many researchers have proposed fatigue models to predict the number of load passes to failure based on the tensile strain at the bottom of the asphalt layer alone. Later on, the models have been improved by including other related variables such as mixture stiffness and volume of bitumen in the mix. The Illinois Department of Transportation presented a strainbased equation for pavement thickness design based on the studies of Thompson and Cation (1986), Gomez and Thompson (1984), and Thompson (1987). The equation had accounted for mixture composition factors, field calibration and split strength characteristics. The following equation was developed for a particular densegraded mixture used at Illinois DOT N, =5x10 6 3 (21) where Nf is the number of repetitions to failure and ; is the maximum tensile strain. Based on extensive laboratory test data covering a wide range of mixtures, bitumen's, and testing conditions, The Royal Dutch Shell Oil Company (Bonnaure et al., 1980) developed separate fatigue models for constant stress and constant strain beam fatigue tests. For constant stress tests, number of repetitions to failure can be defined as Nf = [0.0252PI 0.00126PI(Vb)+ 0.00673Vb 0.0167]5 8EtS14 (22) where PI is the penetration index, Vb is the percentage of bitumen volume in the mix, Sm is the stiffness modulus of the mix in N/m2 and ;t is the tensile strain, which was defined by S028 002 st =[36.43PI1.82PI(V)+9.71Vb 24.04]x106 Sx 1N (23) j y 106 (5x10 23) In a very similar format, the number of repetitions to failure for constant strain can be defined as Nf = [0.17PI0.0085PI(V)+ 0.0454Vb 0.1 12]' 5Sm'18 (24) and the tensile strain ;t was defined by S )036 02 t, =[36.43PI1.82PI(Vb)+9.71Vb 24.04]x0l 6 Sm 0 (25) Sx 1010 106 (25) The Asphalt Institute (1982) also developed a fatigue equation, which accounted for the asphalt volume and the percent air voids in the mixture. The equation was based on the constant stress criterion, which can be expressed as Nf = 0.00432Cs,3 291 E 0 854 (26) where C is the correction factor defined by C = 10 and M =4.84 b 0.69 By assuming a standard mix with an asphalt volume Vb of 11% and air void volume Va of 5%, which makes M = 0 and C = 1. Also multiplying the equation with a factor of 18.4 to account for the difference between laboratory and field condition, the final fatigue criterion then becomes Nf =0.0796 3 291 E* 0 854 (27) where E* is the dynamic modulus. Even though several fatigue equations have been proposed so far, none of them include all the variables that may affect the fatigue failure on asphalt mixtures. The discrepancy between the number passes to failure predicted from laboratory testing and those observed in the field can be very large due to several environmental variables such as varying temperature and resting period. The studies conducted by Francken (1979) found that the fatigue life of asphalt mixtures increased with longer rest periods. Despite the imperfections of the proposed models mentioned so far, the method still provide practical solutions to approximate the failure life of pavements while waiting for a more rational mechanistic design. 2.3 Continuum Damage Approach The continuum damage approach approximates the progression of fatigue cracking in asphalt concrete as a reduction of stiffness due to damage incurred by repeated loading. Lee and Kim (1998a) proposed a viscoelastic constitutive model for asphalt concrete under cyclic loading, which was based on three mechanisms: * Damage due to micro and macrocracks * Relaxation of stress in viscoelastic material * Healing effect of asphalt concrete Lee and Kim employed the extended elasticviscoelastic correspondence principle (Schapery, 1984) to convert viscoelastic stressstrain relationships to stress (o) pseudo strain (ER) relationships such that undesirable relaxation effects in the mixture were removed. According to the theory, small tensile loading that does not cause damage in material yields a linear relationship of stress pseudo strain S= ERER (28) where ER is the reference modulus, which is an arbitrary constant. Higher tensile loading that causes damage in material will change the stress pseudo strain from a linear relationship to a nonlinear response with a hysteresis loop. The pseudo stiffness SR, which defined as the ratio of a stress value to a pseudo strain value at the peak pseudo strain of each cycle, decreases as repeated loading continues. The material is considered to fail when the stiffness is reduced to 50% its original value. Lee and Kim proposed a viscoelastic constitutive model for asphalt concrete under cyclic load, which was applicable for both stress and strain control R R Ca=I R 1+ F(Sl)+G ER, (29) 0 SL where I is the initial pseudo stiffness, sR =R S, which eR is the pseudo strain and R R ER is ER is the shift of pseudo strain (for stress control mode), L = eL which E is the largest ER during the ER history up to the current time, F(Sn) is the function representing the change in normalized pseudo stiffness, Sn = Sp/Sf, which Sp is the damage parameter and Sf is the parameter corresponding to the 50% of the initial pseudo stiffness, The G E0, is the function accounting for the hysteretic behavior of stress pseudo strain relationship, which is the function of the amplitude of pseudo strain (R ),R and e. Lee and Kim (1998b) extended their model by including the healing effect that prolonged fatigue life. The healing function was added to the existing model, which was briefly presented in a simple form a = I[F + G + H] (210) where I is the initial pseudo stiffness, F is the function for change in pseudo stiffness, G is a function account for the hysteretic behavior of stresspseudo strain relationship and the H is the function account for healing effect, eR =R which eR is the pseudo strain and esf is the shift of pseudo strain (for stress control mode). The application of continuum damage approach was first used to evaluate modified and unmodified mixtures for Inchon airport pavements in Korea, Lee, Daniel and Kim (2000). In their work, they used uniaxial tensile test results with cyclic loading to construct the relationship between the number of cycles to failure and strain amplitude. By using elastic layered analysis of the pavement structure, the tensile strain at the bottom pavement can be calculated and used to rank the fatigue life of each mixture. The best mixture was then recommend for the airport project. The continuum damage approach provides another way to assess the fatigue resistance of mixtures. The method was rational and applicable for both stress and strain control modes. The model was suitable for testing and evaluating mixtures for each project. However, by looking at cracks as continuum damage, the method lacks characteristics associated with the presence of a physical crack such as stress concentration ahead of crack tip and stress redistribution due to crack growth. Therefore the approach may only be suitable for design purposes, but not for analysis of fracture. 2.4 Fracture Mechanics Approach The science of fracture mechanics was first introduced by Griffith (1920), explaining a quantitative connection between fracture stress and flaw size. Griffith invoked the first law of thermodynamics to formulate a fracture theory based on a simple energy balance. According to his theory, the crack will increase in size if sufficient potential energy is greater than the surface energy of the material. Since then, several researchers embraced fracture mechanics concepts to quantify fracture toughness of asphalt mixtures and develop fracture mechanicsbased models to predict crack growth in materials. The fracture mechanics approach used in the field of pavements will be reviewed under two subtopics: conventional fracture mechanics and the application of fracture mechanics. 2.3.1 Conventional Fracture Mechanics Asphalt concrete is known as viscoelasticplastic and temperature dependent material. Under certain conditions such as at low temperature and high strain rate, the material behavior can be approximated with elastic analysis; therefore, early work in the field of pavement fracture mechanics assumed that Linear Elastic Fracture Mechanics (LEFM) is applicable. Material characterization with the fracture mechanics approach normally is conducted in two ways: monotonic loading and cyclic loading. Monotonic loading is used to determine critical stress intensity (KIc), which refers to fracture toughness of the tested mixtures. The cyclic loading is also performed to obtain fatigue test data for fitting material parameters to specific fatigue cracking models. There are several test setups that are possible for obtaining the critical stress intensity factor (KIc) for asphalt concrete. Winnie and Wundt (1958) proposed the following equation for determining KI, (PauN) under three point bending test KI = 0.521o2B (211) where o is an applied critical remote bending stress (Pa), B = (da) (mm), d is a depth of the beam (mm), and a is a crack (notch) depth (mm). Ewalds and Wanhill (1986) used a collocation method to derive the formula for determining fracture toughness Kic under three points bending test PI a KIC 3/2 f() (212) bd3/ d where P is a maximum load on the load deflection curve (N), 1 is a span of the beam (mm), b is a width of the beam (mm), d is a depth of the beam (mm), f(a/d) is a function of crack geometry = A/B, with A=3 l/2 1.99 1il f2.153.93 + 2.7 21 and d d) d d) d) B= 2 1+2 (a I a 3/2 d) d) Under normal circumstances, traffic and thermal loading create a stress intensity factor (KI) less than the critical stress intensity (KIo); therefore, no sudden failure occurs. However, a pavement does crack due to repeated loading. Paris and Erdogan (1963) first discovered the power law relationship for fatigue crack growth in material, which is commonly known as Paris law da= A(AK)n (213) dN where da/dN is the crack growth rate (mm/cycle), AK is the range of stress intensity factor (Mpamm.5) during repeated loading condition, and A and n are material parameters determined from fitting fatigue test data. Since the introduction of Paris's law, there have been several proposed models embedded into the da/dN AK relationship. Weertman (1966) developed an alternative semiempirical equation da CK K4 da C=K 4 (214) dN KIc K Ic max where Kmax is the maximum service stress intensify factor. Forman (1967) proposed the following equation da CAKm N1 (215) dN KIc K1 max where C and m are the material constants Elber (1970) proposed a modified Paris law da =/ CAK," (216) dN where AKeff is the effective stress intensity range, defined as Kmax Kop; in which Kop is the stress intensity factor at which the crack open. Klesnil and Lukas (1972) also proposed a modified Paris law to account for a threshold a C(AKm AK ") (217) dN where the threshold (AKth) is a fitting parameter to be determined from a fatigue test. It is important to note that all Eq.213 through Eq.217 were not developed from basic mechanics considerations. They are empirical relationships containing two or more constants, which can be determined from experimental test data. A more mechanicsbased model was developed by Ramsamooj (1980). He used the nonlinear differential equation governing subcritical growth of a crack embedded in an elasticplastic matrix up to the point of gross instability, derived by Wnuk (1971), to formulate a fatigue crack growth model for a beam on elastic foundation da A (K4 K4) (218) dN 24KI2 2 where ot is the yield stress in flexural tension, KI is stress intensity at service load and Ko is the value of stress intensity at endurance limit. Ramsamooj (1991) then again used the same nonlinear differential equation with an estimated value of plastic zone for asphalt concrete A = 0.125(Ki/ot)2 to formulate a general expression for fatigue behavior of asphalt concrete under any configuration of loading and boundary conditions d K K K 2 K2 K dc i In 1 In 1 L+ ] (219) dN 12a 2 Kf Kf K2 K Zhang (2000 and 2001) used fracture mechanics principles with linear viscoelasticity to develop a fatigue crack growth model for crack propagation in IDT with center hole specimen under repeated loading. The approach utilized the concept of a plastic zone ahead of crack tip to accumulate damage for each cycle of loading. For the haversine 0.1 seconds load pulse, with a 0.9 second rest period; the damage due to dissipated creep strain energy (DE) per cycle is determined as follows DE o Dcy = 0avg sin(10t)s, sin(10to)dt (220) cycle where oavg is an average stress in the zone of interest, Sp is a maximum creep strain rate determined from 100 second creep test. According to the energy threshold concept introduced by Zhang (2000), the process zone ahead of the crack tip will accumulate damage, described by dissipated creep strain energy, for each cycle of loading until a dissipated creep strain energy threshold is reached. If the damage level exceeds the energy limit, macrocracks will propagate for the length of the process zone. On the other hand, if the damage level stays below the threshold, only microcracks develop and they are healable after a resting period. With this energy threshold concept, the model predicted the progression of cracks in a stepwise manner similar to the crack propagation observed in laboratory. 2.3.2 Application of Fracture Mechanics Since the fracture mechanics approach has been well accepted for analyzing crack growth in material, many researchers have adopted the concepts and applied them to the field of pavements. Jacobs (1995 and 1996) used fracture mechanics principles to characterize fracture toughness and construct master curves for various tested mixtures. He conducted uniaxial static testing of double edge notch specimen (50x50x150mm) to obtain maximum tensile strength (om) and fracture energy (F). He also conducted the fatigue tests of the double edge notch specimen to determine material parameters A and n. According to his research, he concluded that theoretical derivations for A and n for viscoelastic materials by Schapery (1973, 1975, 1978) appeared to be valid A 2 2 w(t)"dt and (221) 6,g2 2I 2F f n = 2 1+ for force controlled tests; and 2 n = for displacement controlled tests m where Ii is a result of the integration of stress near the crack tip over a small region ahead of the crack tip known as the failure zone, v is a Poisson ratio, D(t) = Do+D2tm with, Do is an initial creep compliance, D2 is the t = 1 sec. intercept of a line drawn tangent to the doublelog creep compliance {D(t)Do} vs. time plot, t is a current time, m is a slope of doublelog creep compliance vs. time plot, At is the period of loading to complete one cycle of loading, w(t) is the wave shape of the stress intensity factor. Ramsamooj (1993) modified an analytical solution for a thin plate resting on elastic foundation by including 3 cracking conditions: crack through in transverse direction; crack through in longitudinal direction and semi elliptical crack at the bottom of the plate. The formulas can be used to evaluate stress intensity ahead of crack tip for those three types of cracking in pavement. Collop and Cebon (1995) investigated the causes of cracking on the top and bottom layer of an asphalt concrete pavement. According to their findings, the transverse contact traction may cause a short crack (10 mm) on the pavement surface and the repeated bending tensile stress caused crack at the bottom. They then formulated an analytical solution for determining stress intensity factor for a surface crack in a semiinfinite plate subjected to a general remote stress. They also derived a formula based on parametric studies to quantify the stress intensity factor due to traffic and thermal loading for various pavement structures at progressive crack growth in pavement. Myers et al. (2000 and 2001) also investigated the mechanism of surface initiated cracking along the wheel path but they used measured tire contact stresses instead of the traditionally assumed circular load. According to their findings, the surface cracking was caused by transverse shear stress of a radial tire and perhaps combined with thermal stress due to rapid cooling. Then they used a fracture mechanics approach with the finite element software ABAQUS to study the propagation of top down cracking. By assuming cracks propagate perpendicular to the major principal (tensile) stress, they found the cracks grew down vertically in the pavement and then bent 30 degrees toward the wheel load. The predicted crack path was found similar to what has been observed in the field. They also performed parametric studies to quantify the stress intensity factor ahead of the crack tip for various load positions, stiffness ratios and thickness ratios for progressive crack growth in a flexible pavement. Sangpetngam (2003) has adopted a crack growth model developed by Zhang (2000 and 2001) into a displacement discontinuity boundary elementbased pavement fracture simulator. He used crack tip elements to improve accuracy of stress distribution ahead of crack tip and replaced the process zone with yield element. The implementation of crack growth model in the numerical scheme increases the application of fracture mechanics for predicting crack growth in asphalt mixture for various geometries, loading and boundary conditions. The use of fracture mechanics approach to analyze cracking in material is suitable for problems that experience either a single line or a few distinct lines of fracture. The mechanics of stress concentration and stress redistribution ahead of crack tip that play an important role to the failure in material are very well established. However, the accurate modeling of actual distributed crack bands found in granular materials is still missing in conventional fracture mechanics. 2.5 Displacement Discontinuity Method with Tessellation Schemes The displacement discontinuity method developed by Crouch (1983) has been extensively used in the fields of rock mechanics and geological engineering. The method has the potential to be an analytical tool for assessing cracking in granular materials such as asphalt concrete. The clear advantage of using this method to analyze crack type problem is that cracks can be explicitly modeled in the medium and the accuracy of stress ahead of crack tip is very good. The method can be modified to simulate cracking in granular material by creating a numerical model with two types of elements: exterior boundary elements and potential crack elements. The exterior boundary elements are placed along the boundary of a problem to simulate the edge of specimen and potential crack elements are randomly placed inside specimen to simulate predefined crack paths, which normally is assumed to be along the grain boundary or perhaps through grain as well. The displacement discontinuity boundary element method can be used coupled with various tessellation schemes. The use of tessellations to represent granular structure in simulation of fracture process zone is very well accepted. It improves the realism of predicted failure mechanisms at the particle level. Two basic tessellation schemes Delaunay and Voronoi have been used in various fields. Both tessellations can be used to simulate polycrystalline of ductile material such as OstojaStarzewski (1987), Ostoja Starzewski and Wang (1989), Van der Burg and Van der Giessen (1993), Helms, Allen and Hurtado (1999). The tessellation schemes are also applicable to simulate granular structure of brittle rocks such as Napier and Peirce (1995), Napier et al. (1997), Sellers and Napier (1997) and Steen, Vervoort and Napier (2000). The suitable choice of tessellations to represent granular structure depends on the realistic looks of failure pattern and observed responses. Napier and Peirce (1995) invented a new boundary element solution technique, termed the multiplee method" for solving multiple interacting crack problem that involved several thousands boundary elements. They used the new technique to study the different failure mechanisms of a rectangular rock sample under displacement control using two tessellations schemes (Delaunay and Voronoi) for three levels of grain densities. It appeared that Voronoi assemblies are less prone to shed load than the Delauney triangulations. With increasing the density of Voronoi polygons, it seemed to not change this conclusion. Steen, Vervoort and Napier (2001) took a similar approach but with 3 tessellation patterns, Delaunay, Voronoi and Voronoi with internal fracture path to simulate a confined compression test of a rock sample. It found that the use of Voronoi tessellations with internal fracture paths best simulated the formation of shear band in the specimen. The Delaunay tessellation scheme was the second best while the Voronoi tessellation was not able to predict observed shear banding. The experiments also tested the use of two failure criteria; Rankine criterion (Chen and Han 1988) and Coulomb failure criterion to identify the appropriate failure law that allowed the formation of shear bands. The result of the simulation revealed that only Coulomb failure criterion enabled localization of shear band. Birgisson et al. (2002a) first introduced the displacement discontinuity boundary element method to simulate crack growth in Superpave IDT specimens. They used exterior boundary elements to create a 2D plain stress specimen and randomly laid down potential crack elements forming Voronoi tessellation inside the specimen. With an appropriate set of material parameters for local failure at potential crack elements, the numerical prediction can capture stress strain responses and crack pattern. Birgisson et al. (2002b) continued studying the tessellation schemes and crack growth rules for simulating crack growth in asphalt mixture. They found that the Voronoi tessellation with internal fracture path provided realistic simulation of crack growth in asphalt mixture because aggregate is allowed to break down at high load level, which results better fit of stress strain responses. The results of using different two crack growth rules, sequential and parallel, are not significant, but the parallel crack growth rule that activates several cracks at a time is more efficient than the sequential crack growth rule that activates crack one at a time. Soranakom et al. (2003) used the same numerical schemes as Birgisson et al. (2002b) to evaluate mixture properties. They found that the method was capable of evaluating mixture properties with acceptable accuracy. The predicted results matched vertical compressive stress strain curve, horizontal tensile stress strain curve, crack pattern, tensile strength and fracture energy of the mixes. The displacement discontinuity method with tessellation scheme as predefined crack paths has demonstrated its ability to simulate crack propagation in granular material such as asphalt concrete. The use of appropriate tessellation to represent grain structure improved the realism of distributed crack in granular material that macrocracks were accompanied by several microcracks. 2.6 Summary From the literature review presented above, it was concluded that the most appropriate approach for studying the cracking behavior of asphalt mixtures was the displacement discontinuity method with tessellation schemes. The approach provides an elegant way to model aggregate and mastic separately such that the detail of aggregate structure and strength of mastic can be investigated. The approach can be used to identify the key factors that control the cracking mechanism in asphalt mixture. The understanding of microstructure and its cracking behavior will provide the rational recommendation for improving crack resistant mixture. CHAPTER 3 DISPLACEMENT DISCONTINUITY METHOD WITH TESSELLATION SCHEME 3.1 Overview Displacement discontinuity method (DDM) is an indirect boundary element method developed by Crouch (e.g. Crouch and Starfield, 1983). The displacement discontinuity boundary element method is uniquely well suited for dealing with problems involving fracture and crack growth. Even though the method is inherently simple to use and very flexible for many kinds of problems in various fields of engineering, it has received little attention compared to the wellknown finite element and finite different methods. These are some reasons: * The misconception that boundary element method is derived for a particular class of problems, and it cannot be applied to other classes of problems. The method may in fact be easily modified for solving other problems. * Theoretical papers on boundary element method are viewed by many engineers as being somewhat abstruse and sometimes difficult to understand, due to the heavy mathematical detail required. Thus, the complexity of the mathematics involved in deriving the system of equations may have been an obstacle in the proliferation of the boundary element method. The intention of this chapter is to provide a theoretical background of the displacement discontinuity boundary element method to solve problems involving fracture of materials, particularly cracking in the indirect tension specimen. 3.2 Displacement Discontinuity Method A physical crack in material can be viewed as a discontinuity surface. The displacement discontinuity method assumes displacements in a body are continuous everywhere except at a line of discontinuity (Figure 31). When the displacement crosses the line of discontinuity, its value jumps by the amount of the displacement discontinuity Di, in which its components in local axis coordinates yz are D q) = u, (y,O)u,(y,O ) i= y, z; b< yq where z = 0+ is the positive side and z = 0 is the negative side of the discontinuity element. b Cl Q(yq,0) c2 b z vp(yp,Zp) Figure 31. Displacement discontinuity element in local coordinate (yz). A linear discontinuity element has a total length of 2b with two collocation points (ci and c2) at y = + (/2/2)b. The distance r is measured from the field point p(yp,Zp) in the domain to the source point Q(yq,0) at boundary. In twodimensional plane strain problems, if the line of discontinuity has a length of 2b, centered on the yaxis of a local coordinate system yz, and has normal vector components ny and nz along the surfaces, it can be shown (Peirce and Napier 1994) that the contribution of a given element to the total displacement components at point p is Dy(y )n y u ( p) 1 bL M N 1 (P)LY M NY (D(yq)n, +D (yq) dq (32) D, (Yq)q, where Uy, Uz are the local components of the displacement vector, v is Poisson ratio and the symbols in the matrix are given by Ly = (1 v)Y~ + (2 v)YW MA = ZH + (1 v)Yz Ny = VJy' (1 v)Wyz L = (1 v)Yw + O , h = (1 v)fYnt Vp N, =(2 v)Y +(1 v)Y The biharmonic function for plain strain problem is given by (33) 1 2 ogr2) Y=(r r2 logr2) 2 and r = (y y, It can also be demonstrated (Peirce and Napier 1994) element to the total stress tensor components at point p is a (y, q) b T = 4 aJ a (yv, q ) = f . Y TY a (y,, q ) Y Y_ Y 0zz~ypp) JL,yyzz ,~yyz ,yyy ) +z (34) that, the contribution of the Dy(yq)ny D, (y )n + D (y )nY dy D (y,)n (35) For numerical implementation, the displacement discontinuity Di is approximated by a polynomial function, y = ao + alx + a22 + ... + anX". More accurate results can be obtained by using several terms in the approximation, but it will increase computing time. It is sufficient for most practical applications to approximate the displacement discontinuity with a linear function as was done in the Discontinuity Interaction and Growth Simulation (DIGS) program, (Napier, 1990; Napier and Hildyard, 1992; Napier ) and Pierce, 1995a; Malan and Napier, 1995; Kuijpers and Napier, 1996; Napier at el., 1997; Napier and Malan, 1997). The linear variation of discontinuity can be written as D, (yq) =a, + b,yq ai and bi are constants (36) By substituting Eq.36 into Eq.35, and carrying out the mathematical manipulation, the analytical solution for normal stress on the yaxis of the local coordinate system (yz) will have the form: E 1 1 y+b O. (y,O) 8= r('2) 2(a + iP )i yb)+ fz log(Y )2 (37) 8^(1i2) y+b yb yb It is obvious from Eq.37 that the analytical normal stress yzz approaches infinity as y approaches the tip of the element y = + b. However, in boundary element formulations, stress in Eq.37 is evaluated at suitable collocation points y = + c, namely the Gauss Chebyshev points (Crawford and Curran, 1982) y, =cos(2i 1) i = 1, 2 ..., n (38) 2n The stresses at collocation points will have finite values and are solvable using a numerical algorithm. 3.3 Numerical Implementation For numerical implementation, the displacement discontinuity method (DDM) employs the known solutions of the discontinuity surface to formulate a system of governing equations. For a problem with one crack in the infinite elastic body without far field stresses, the general system of governing equations can be written as: = Z(ADj+ + ADI) J S S s (39) ' = (AD+ A DI) Jn ns nn n where o and a, are the shear and normal stress of the element i. A", A" A" and A' are the influence coefficients due to element j on element i, and D' and D are the displacement discontinuity components of element j, which are the unknowns of the system. In simulations of crack interaction problems, a displacement discontinuity element will slide when driving shear stress exceeds shear strength or open up when applied tensile stress exceeds tensile strength. Even if the DDM was initially developed for an open crack, it can easily be extended to include contacting crack surfaces and sliding cracks in softening mode. When two crack surfaces are in contact, the shear and normal stress components o and o, depend on the stiffness (Ks and Kn) and the displacement discontinuity components (D' andD' ). The relationships can be written in matrix forms as follows o = KDg c = KD' (310) By substituting Eq.310 into Eq.39 and rearranging the terms so that the unknowns are on the right hand side, the system of governing equations then becomes 0= X(A"Dl + AtD) KD 0 = (AADj + A'D) KD (311) When an element mobilizes, the crack surface will deform according to softening models, which will be described shortly in the next section. The residual strength of the element i, namely oa andao, can be assumed to decrease as a function of the displacement discontinuities D' and D', which can generally be expressed by Sf (DI, D) a =f n(D, D ) (312) Substituting Eq.312 in Eq.39, and rearranging the unknowns to the right hand side, the system of equations becomes 0 = Y(A Ds + AsDI) f(D, D') 0 = C(AsDJ + A D) f,(DI,D') (313) Finally, a set of algebraic equations that consist of known driving forces, influence coefficients and the unknown displacement discontinuities can be written in a matrix form, {0} = [A]{D}. Since the mobilization of cracks is associated with a softening model in which stresses depend on the unknown discontinuities, an iterative technique needs to be employed in solving the equations. DIGS use an efficient multipole method (Napier and Peircel995a) for iterative solving. Once displacement discontinuities at all boundary elements have been determined, displacements and stresses at any designated points can be computed by using Eq.32 and Eq.35. 3.4 Example of a Crack Growth Model using DIGS A numerical model for a Superpave indirect tension test (IDT) specimen consists of two types of elements: exterior boundary elements and potential crack elements. The exterior boundary elements are regular boundary elements that are employed to construct an edge of a specimen. Potential crack elements on the other hand allow the element to mobilize (slip or open up) when the stress exceeds a failure limit. Therefore, they are randomly placed inside the specimen to form any desired tessellations to represent granular structure. Figure 32 shows examples of using three tessellations, Delaunay, Voronoi and Voronoi with internal fracture paths, to simulate the aggregate structure of an IDT specimen. During the analysis, at each load step, stresses are computed at collocation points inside potential crack elements. The stress state at these locations is then checked against a failure limit to determine the activation of a crack. In DIGS, a nonlinear failure law, shown in Figure 33, is adopted for the cracking criterion. The failure law consists of a linear portion in the compression region (similar to the MohrCoulomb failure law) and a power law curve in the tension region, with a continuous slope at on = 0. The linear portion has the following form r = S tan(Q)c, when on <0 (314) where S is the cohesion, 4 is the friction angle and on is the normal stress (tension positive) across the discontinuity. The power law curve is defined by: r = a(ac a,)b when on > 0 (315) where ct is the tensile strength and the two constant "a" and "b" are the parameters that fit the power law curve to the tension cutoff To at shear stress as = 0 and the linear portion at normal stress on = 0. /r / B) Voronoi and C) Voronoi with internal fracture path. Each tessellation has 503 particles, which yields an equivalent aggregate size of 6.8 mm. '* "' i . iH.<9^ tx TO',^WT T i )<' i ,: ._ *.... ,1; ^ / ... .tM^ '* .w . A B C Figure 32. Tessellations for granular structures. Aggregate structure in indirect tension specimen can be represented by three basic tessellations: A) Delaunay, B) Voronoi and C) Voronoi with internal fracture path. Each tessellation has 503 particles, which yields an equivalent aggregate size of 6.8 mm. Failure envelope Residual failure Stress Co T = a(ct Cn)b I Normal To Stress ct= To Tsoft Dn =0 if Dn>DNCR Figure 33. Mohr Coulomb type of failure criterion for determining crack mobilization. For every load step, a crack search algorithm is performed to detect cracks that may occur in the specimen. There are two crack search algorithms available in DIGS, namely parallel and sequential. The parallel crack search activates all elements where the stress exceeds the failure limit. Alternatively, the sequential crack search activates only the most overstressed element found in system; thus, only the most critical crack element is activated. At the end of each search step, the detected cracks are added to the system and resolved for new stress pattern to determine new cracks in the next search step, until no more cracks found. When potential crack element turn into a cracked state, the cohesion S is assumed to weaken as a linear function of the slip Ds S = C, Cof ID, 1 (316) where Co is the original cohesion intercept and Csoft is the cohesion softening slope. Similarly, the tensile strength ot is also assumed to weaken as a linear function of the opening displacement Dn Ct = To T,, I D,I (317) where To is the tension cutoff and Tsoft is the tension softening slope. When crack slip occurs, the tensile strength is implicitly degraded as the cohesion softens, congruently with the extent of cohesion softening. The tensile strength completely gone when the opening crack (Dn) exceeding the opening crack limit (DNCR) CHAPTER 4 NUMERICAL EXPERIMENTS OF THE SUPERPAVE INDIRECT TENSION TEST USING THE DISPLACEMENT DISCONTINUITY METHOD WITH A TESSELLATION SCHEME 4.1 Overview This chapter provides the results of a series of parametric studies that systematically evaluate the factors and parameters needed to successfully simulate the SUPERPAVE indirect tension (IDT) strength test. The appropriate set of tessellation and crack growth rules are identified, and the use of different grain sizes to represent actual gradation was also conducted to identify the optimal size that would yield the most consistent result for numerical simulations. Finally, a sensitivity study of material parameters defined for mastics and aggregate were performed to study their effects on the loaddeformation and cracking behavior in the IDT strength test. Since the displacement discontinuity method can be used to simulate different kinds of materials such as rock, concrete and asphalt concrete, a typical IDT test result of granite mixture (Nova Scotia NS619) was used to provide the target for modeling asphalt concrete. First, the input parameters (global and micromechanicsbased parameters) for the displacement discontinuity boundary element method were calibrated so that the resulting predicted stressstrain curves matched the experimental curves. Then each component and material parameter was altered from the calibrated fit model to observe the sensitivity of the predicted results (stressstrain curves) and identify the best combination. The scopes of numerical studies were summarized below. * Determination of appropriate tessellation schemes and crack growth rules * Determination of uniform aggregate size to represent gradation * Sensitivity analysis of material parameters defined for mastics * Sensitivity analysis of material parameters defined for internal fracture path (aggregate): After performing all analysis stated above, the appropriate set of tessellation, crack growth rules and material parameters were proposed for numerical simulation of indirect tension strength test. 4.2 Laboratory Experiment of Nova Scotia Mixture A SuperPaveTM IDT strength test at 0C was performed to obtain the experimental results used in this study. The IDT specimen was prepared according to the SuperPaveTM specifications. The aggregate used was granite from Nova Scotia. The gradation shown in Table 41 was used to design a 19mm (0.75 in) nominal maximum aggregate size coarsegraded mix. AC30 asphalt was used with the design content of 5.5% at Ndesign = 109 revolutions in the SuperPaveTM gyratory compactor. The asphalt and the aggregate were mixed at 150 OC (3020F) and subjected to shortterm oven aging for two hours at 1350C (2750F). For the strength testing, specimens were compacted to produce three 4500 g (9.9 lbs), 152.4 mm (6 in) diameter gyratory compacted specimens with an air void content of 7% (+ 0.5%). All aggregate and volumetric properties of the mixture met SuperPaveTM requirements. Table 41. Gradation of the Nova Scotia mixture Sieve size (sieve number) Percent passing 25 mm (1) 100 19 mm (3/4) 97 12.5 mm (1/2) 83 9.5 mm (3/8) 66 4.75 mm (#4) 38 2.36 mm (#8) 23 1.18 mm (#16) 18 600 [tm (#30) 14 300 [tm (#50) 10 150 [tm (#100) 6.5 75 tm (#200) 3.5 The gyratory compacted specimens were sliced to obtain SuperPaveTM IDT test specimens, with an approximate diameter of 152.4 mm (6 in) and 50.8 mm (2 in) thickness. The IDT strength test was performed according to the SHRP Indirect Tension test procedure, presented by Roque et al (1997). A summary result of IDT strength test is shown in Table 42. The applied load, horizontal and vertical deformation measured at each face of a sample can be represented by horizontal and vertical stress strain curves at the center using IDT plane stress formula as follows Gh = 2P/(7nDt) Eh = A/Lgauge o, = 6P/(tDt) Ev = AV/Lgauge (41) where, ch is a horizontal stress at the center of specimen, P is a total load applied, D is a diameter of a specimen, t is a thickness of a specimen, Sh is an average horizontal strain over a horizontal strain gauge, AH is a corrected horizontal deformation at a horizontal strain gauge, Lgauge is a gauge length, o, is a vertical stress at the center of specimen, e, is an average vertical strain over a vertical strain gauge, and AV is a corrected vertical deformation at a vertical strain gauge. Finally, the horizontal and vertical stress strain curves on each face of the sample were averaged and used for comparison to the numerical prediction, which will be described shortly in the next section. Table 42. Material properties from IDT strength test Parameters Value Tensile strength (MPa) 2.05 Fracture energy (MPa) 0.83 Poisson ratio 0.35 Load STop plate Strain gauges 38.1 D q\ I 38 1 D mm mn mDm mm 38.1 3 mm Bottom plate 38.1 (25.4 xT mm) mm A B Figure 41. Indirect tension strength test: A) Experimental setup and B) Displacement discontinuity model. 4.3 Numerical Simulations The experimental setup and the associated numerical model for the Superpave IDT strength test are shown in Figure 41. The specimen size is 152.4 mm (6") in diameter and varying in thickness (T mm). Two end plates 25.4 mm x T mm (1" x T") at the top and the bottom are used to protect the damage to specimen due to high concentrated load. At the center of specimen, two strain gauges 38.1 mm (1.5") are placed in the vertical and horizontal directions to measure vertical and horizontal deformations in the central zone. Figure 41B shows an IDT numerical sample using Voronoi with internal fracture path to represent aggregate structure. Four field points (acting like virtual strain gauges) are placed at the locations of the actual strain gauges to monitor the displacements at each load step. The perimeter of the numerical model is comprised of 152 exterior boundary elements. Of those, 8 elements are placed at the top and forced to move downward to simulate displacement control. Other 8 elements at the bottom are fixed to provide a static condition. The remaining 136 elements along the circumference are specified with zero traction to simulate traction free surface. The displacement controlled loading was performed in 12 equal steps for a total displacement of 0.36 mm. For perfect displacement control, the predicted stress distributions at the load plates are nonuniform and they become significantly fluctuated when crack approaches the end plates. For the comparison to the actual testing results, the applied load was approximated by the average stress at the top and bottom plates multiplied with the plate area. The horizontal and vertical deformations in the center part of the specimen can be obtained by determining the relative movement of the displacements requested at the virtual strain gauges. Finally, the predicted stress strain curves for horizontal and vertical direction were computed and represented in the same way as the experiment by using Eq.41. 4.4 Determination of Suitable Tessellation Schemes and Crack Growth Rules The identification of an appropriate tessellation scheme and crack growth rule was an iterative process since the effects of other relevant components such as particle size and material parameters for the mastic and aggregate were not yet fully known. Several attempts to model IDT strength test Birgisson et at (2002 ab) and Soranakom (2003) had previously been performed in which it was shown that it was possible to fit the stress strain curves and correctly predict the crack pattern observed in the laboratory. However, some issues still needed to be studied further, namely the establishment of a representative optimal Voronoi particle size for a given gradation, and other issues related to the "tuning" of the micromechanical material parameters. In the following, the insight gained from previous work (Birgisson, et al., 2002ab) will be used to identify a suitable tessellation scheme and crack growth rule. 4.4.1 Numerical Experiment Three tessellations Delaunay triangulation, Voronoi polygons and Voronoi polygons with internal fracture paths, shown in Figure 32, were used to model the aggregate structure of the Nova Scotia mixture. Each tessellation has 503 particles filling the circular area of 18240 mm2, which yields an equivalent aggregate size of 6.8 mm or at D50 (aggregate size at 50% passing of the gradation). All three tessellation schemes were used to simulate the Superpave IDT strength test. The predicted stressstrain curves and crack patterns were then compared to the experimental results obtained for the Nova Scotia mixture. For each tessellation scheme, two crack growth rules, parallel and sequential, were compared. The modeling matrix, global and local material parameters used for the numerical experiments are provided in Table 43 through Table 45. It is note that material parameters defining mastic strength and its residual strength for Delaunay tessellation are substantially higher than those for Voronoi and Voronoi with internal fracture path in order to ovoid unstable collapsing and produce the best fit to the experimental stress strain curves. Table 43. Modeling matrix for the study of suitable tessellation schemes and crack growth rules Tessellation schemes Delaunay Voronoi Voronoi with internal fracture path Crack growth rule Parallel Sequential Parallel Sequential Parallel Sequential Equivalent aggregate size (mm) 6.8 6.8 6.8 6.8 6.8 6.8 Number of potential crack element for mastic 720 720 1408 1408 1408 1408 Number of potential crack element for fracture path 2717 2717 Table 44. Global linear elasticity parameters for the study of suitable tessellation schemes and crack growth rules Parameter Value Young modulus (MPa) 7500 Poisson ratio 0.35 Table 45. Local material parameters for the study of suitable tessellation schemes and crack growth rules Tessellation schemes Delaunay Voronoi Voronoi with internal fracture paths Internal fracture Material parameters Mastic Mastic Mastic pahts Tension cutoff, To (MPa) 4.00 2.30 2.30 6.50 Opening crack limit, DNCR (mm) 0.15 0.11 0.11 0.07 Tension softening slop, Tsoft (MPa/mm) 26.67 20.91 20.91 92.86 Cohesion, Co (MPa) 4.00 2.40 2.40 6.50 Residual cohesion, CR (MPa) 1.00 0.15 0.15 0.15 Cohesion softening slop, Csoft (MPa/mm) 0.50 5.00 5.00 1.00 Friction angle, 4o (degree) 45 38 38 40 Residual friction angle, 4R (degree) 44 30 30 32 4.4.2 Results from Numerical Simulations The result of the numerical simulations will be discussed in the next four sections entitled: Crack growth rules, Delaunay tessellations, Voronoi tessellations, Voronoi tessellations with internal fracture path, followed by a summary of findings for the optimal crack growth rule and suitable tessellation scheme. 4.4.2.1 Crack Growth Rules Figure 42 shows the amount of cracking predicted by the two crack growth rules: parallel and sequential. Both predicted approximately the same number of cracks for all load steps up to the load step at the peak of the stressstrain curves (load step 10 for the Delaunay tessellation scheme, load step 9 for Voronoi tessellation scheme and Voronoi with internal fracture path scheme (Figure 43 45). In the post peak region, the sequential crack growth rule predicted fewer cracks than the parallel did. Parallel and sequential crack growth rules also yielded approximately the same stress strain responses up to the peak, but started to deviate slightly in the post peak. The small differences in the postpeak stressstrain response were due to the differences in predicted postpeak crack patterns. Figure 46 shows a comparison of the predicted crack patterns obtained with the three tessellations and the two crack growth rules studied. All the predicted crack patterns showed similar patterns a major vertical opening crack in the middle region of the specimen, accompanied by several smaller size cracks, distributed over a wider band through the center of the specimen. Overall, two crack growth rules yielded approximately the same result for both predicted stress strain responses and crack patterns. However, the sequential crack growth rule that activates the most overstressed element one by one is more computationally intensive than the parallel crack growth rule that activates all overstressed elements at the same time. Thus, in the following, the parallel crack growth rule is used in all remaining simulations. 800 8 A Delaunay (Parallel) 700 T 600 A " Delaunay (Sequential) 500 S400 Voronoi (Parallel) o 400 0 300  E 0 Voronoi (Sequential) z 200 100 i Voronoi with internal fracture path 0 .I.... (Parallel) 0 1 2 3 4 5 6 7 8 9 101112 0 Voronoi with internal fracture path Load Step (Sequential) Figure 42. Predicted number of cracks at each load step. Three tessellations, Delaunay, Voronoi and Voronoi with internal fracture path were used with two crack growth rules, Parallel and Sequential 3.50 3.00 2.50 2.00 1.50 1.00 0.50 0.00 A' 0.0 1.0 10.0 8.0 6.0 4.0 2.0 A Delaunay (Parallel)  A Delaunay (Seqential) Experimental  Delaunay (Parallel)  A Delaunay (Sequential) Experimental 0.0 1.0 2.0 3.0 4.0 5.0 6.0 Vertical Strain (x 1000) B Figure 43. Comparison of predicted and measured horizontal and vertical stressstrain curves, using Delaunay tessellation: A) Horizontal stressstrain response at center of specimen and B) Vertical stressstrain response at center of specimen. 2.0 3.0 4.0 5.0 6.0 Horizontal Strain (x 1000) A 3.50 S3.00 a 2.50 S2.00 C, S1.50 .N 1.00 I 0 I 0.50 0.00  Voronoi (Parallel)  1l Voronoi (Seqential) Experimental 0.0 1.0 2.0 3.0 4.0 5.0 6.0 Horizontal Strain (x 1000) A 10.0  Voronoi 10.0 r 8.0 (Parallel) 6.0 S f1 Voronoi 4.0 (Sequentila) 2.0 Experimental 0.0 1.0 2.0 3.0 4.0 5.0 6.0 Vertical Strain (x 1000) B Figure 44. Comparison of predicted and measured horizontal and vertical stressstrain curves, using Voronoi tessellation: A) Horizontal stressstrain response at center of specimen and B) Vertical stressstrain response at center of specimen. *Voronoi with fracture path (Parallel)  0 Voronoi with fracture path (Sequential) Experimental 1.0 2.0 3.0 4.0 5.0 6.0 Horizontal Strain (x 1000) A Voronoi with fracture path (Parallel)  0 Voronoi with fracture path (Sequential) Experimenta I 0.0 1.0 2.0 3.0 4.0 5.0 6.0 Vertical Strain (x 1000) B Figure 45. Comparison of predicted and measured horizontal and vertical stressstrain curves, using Voronoi with internal fracture path tessellation: A) Horizontal stressstrain response at center of specimen and B) Vertical stressstrain response at center of specimen. 3.50 3.00 2.50 2.00 1.50 1.00 0.50 0.00 10.0 8.0 6.0 4.0 2.0 I Mesh D Mesh D Parallel crack growth rule A Sequential crack growth rule NT, Mesh V Mesh F Parallel crack growth rule B si^ Parallel crack growth rule C Sequential crack growth rule Sequential crack growth rule Figure 46. Crack patterns predicted by two crack growth rules of three tessellations: A) Delaunay, C) Voronoi and C) Voronoi with internal fracture path. /' 4.4.2.2 Evaluation of Delaunay Tessellation Scheme Figure 46A shows a crack pattern predicted by a Delaunay tessellation scheme. The dominant large crack pattern is similar to those predicted using the Voronoi tessellation schemes (Figure 46B and 46C). The narrower band of distributed cracks, as compared to the Voronoi schemes is simply the result of using higher strength and residual strength of the mastics. In Delaunay scheme, cracks appearing on one side of a triangular particle were able to work their way to the other side easily, creating floating particles (particles completely surrounded by cracks), which weaken aggregate structure in carrying compression load. In some cases, these "floating particles" occurred in the area of the virtual strain gauge points, resulting in almost a rigid body motion. This unstable movement could be minimized by increasing the residual strength parameters of the mastic in the model. In order to use the Delaunay tessellation structure to model the IDT strength test, substantially higher mastic strength than Voronoi tessellations and Voronoi tessellations with internal fracture path was needed. Despite using the best possible set of material parameters, the Delaunay can only obtain a rough match to the experimental stress strain responses as shown in Figure 43, the predicted stressstrain curves were not smooth and did not fit very well. The parameters for residual strength were unreasonable high for such a brittle material asphalt concrete. Therefore, it was concluded that the use of Delaunay tessellation scheme to model the aggregate structure for the Nova Scotia asphalt mixture studied was unsuccessful. 4.4.2.3 Evaluation of Voronoi Tessellation Scheme Figure 46B shows the crack pattern at the final load step predicted using the Voronoi tessellation scheme to represent the aggregate structure of the Nova Scotia mixture. Multiple cracks formed a wide band of cracks, within which, the material appeared to form a columnar structure defined by the overall cracking pattern. Those "broken columns" were still able to provide some strength in carrying compressive load from the top to the bottom of the specimen. This mechanism resulted in a simulated compressive softening of the ascending (near peak) portion of the predicted vertical stress strain curve (Figure 44B). At higher load levels, the number of distributed cracks and the size of opening crack increased, resulting in a rapidly decreased strength. Once the major opening cracks occur, the strength of specimen was diminished. By specifying appropriate values of residual strength parameters, the simulations were able to capture the post peak stressstrain response. Finally, the simulations of IDT strength test using Voronoi tessellation schemes revealed a difficulty in tuning the material properties for the mastics to adequately match experimentally observed peak and post peak load deformation responses. 4.4.2.4 Evaluation of Voronoi with Internal Fracture Paths Tessellation Scheme Figure 46C shows the predicted crack pattern at the final load step predicted using a Voronoi tessellation scheme with internal fracture paths to represent aggregate structure. The predicted cracking patterns were similar to that obtained from the Voronoi tessellation scheme, except that some of the aggregates broke down at high load levels. The additional parameters for defining the aggregate can be modeled such that the predicted stress strain responses (Figure 45) can fit to the experimental results better. Finally, the Voronoi tessellations with internal fracture paths seemed to reasonably represent the observed load deformation behavior of the granular structure of the Nova Scotia mixture. 4.4.3 Summary of Findings on the Determination of a Suitable Crack Growth Rule and Tessellation Scheme From the results presented, the use of Voronoi tessellations with internal fracture paths and the parallel crack growth rule resulted in reasonable simulations that corresponded with experimental stress strain curves and fracture patterns. 4.5 Determination of Average Voronoi Aggregate Size to Represent Gradation 4.5.1 Numerical Experiment This study was conducted to identify the optimal aggregate size for the simulation of IDT strength test. A Voronoi tessellation scheme with internal fracture path (medium size 8.0 mm) was first tuned to match the experimental stress strain curves for the Nova Scotia mixture to obtain the material parameters needed to describe the mastic and the internal fracture paths of the aggregates. Once the appropriate material parameters for the asphalt concrete had been determined, a series of numerical IDT strength test simulations was performed using five uniform aggregate sizes 4.5, 6.0, 8.0, 10.0 and 12.0 mm. For more reliable test result, three samples were performed to obtain average value for each aggregate size. The averaged ultimate tensile strength and variance of residual error of the tensile stress strain curves were used to evaluate the size effects and the associated variability. All details of the numerical models and material parameters used for the study of the aggregate size are listed in Table 46 through Table 48. Table 46. Modeling matrix for the study of size effect: Agg. 4.5 mm Agg. 6.0 mm Tessellation Crack growth rule Equivalent aggregate size (mm) Exterior boundary elements Sample No. Potential crack elements for mastic Potential crack elements for fracture path in aggregate Voronoi with internal fracture path Parallel 4.5 152 Voronoi with internal fracture path Parallel 6 152 Agg. 8.0 mm Voronoi with internal fracture path Parallel 8 152 Agg. 10.0 mm Voronoi with internal fracture path Parallel Agg. 12.0 mm Voronoi with internal fracture path Parallel 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 3129 3123 3128 6168 6157 6163 1826 1837 1829 3551 3569 3554 1012 1165 1015 1923 3090 1929 681 685 682 1261 1269 1263 517 505 505 933 910 910 Table 47. Global linear elasticity parameters for the study of size effect Parameter Value Young modulus (MPa) 7500 Poisson ratio 0.35 Table 48. Local material parameters for the study of size effect Potential crack To DNCR Tsoft Co CR Csoft io tR element (MPa) (mm) (MPa/mm) (MPa) (MPa) (Mpa/mm) (degree) (degree) Mastic 2.40 0.11 21.82 2.40 0.20 1.00 38 30 Internal fracture path 6.40 0.10 64.00 6.40 0.30 1.00 38 30 4.5.2 Numerical Results Figure 47 shows the averaged predicted horizontal and vertical stress strain curves for each aggregate size studied. Within the range of uniform aggregate sizes used, a small size effect is noticeable, where a larger aggregate size results in a stiffer response. This could be explained by the fact that the larger the aggregate size the fewer the resulting cracks in the specimen. Fewer cracks yield less compression softening of the vertical compressive stress strain curves and resulting in stiffer responses. Similar trend were obtained for the averaged values of the ultimate tensile strength (Figure 48), where the larger the aggregate size the stronger the specimen. A change of size from 4.5 mm to 12.0 mm increased the predicted strength by about 11%. The small increasing of strength could be explained by the observed cracking behavior and the differences in the formation of major crack bands. Smaller average aggregate sizes tend to yield more numerous crack planes in the specimen, such that the formation of a crack band can happen easier and sooner than for the case of the larger aggregate, where fewer crack planes are observed. Figure 49 shows the variance of the residual error of the tensile stress strain curves plotted against the aggregate size. This variance is a measure of the variability among the three tested samples referenced to the experimental data. The variability of the test results is caused primarily by the initiation of distributed cracks and the formation of major crack bands. Small aggregate size yields more distributed cracks that can form more crack patterns than the large aggregate sizes that yield fewer cracks. However, the formation of a large crack has a more pronounced effect on a specimen than the formation of a smaller crack. In this experiment, it was found that the use of aggregate sizes between 5.5 to 8 mm (or at 42 to 57% of passing) yielded the least variability (below 0.011 MPa2). Therefore, the optimal aggregate size at 40 60% of passing of gradation (D40D60) was recommended for reliable numerical simulations of the Nova Scotia coarse mix asphalt specimen. 3.50 3.00 2.50 2.00 1.50 1.00 0.50 0.00 0.0 1.0 2.0 3.0 4.0 5.0 Horizontal Strain (x 1000) A Agg. 4.5 mm mAgg. 6.0 mm AAgg. 8.0 mm *Agg. 10.0 mm Agg. 12.0 mm  Experimental 0.0 1.0 2.0 3.0 4.0 5.0 6.0 Vertical Strain (x 1000) B Figure 47. Comparison of predicted and measured horizontal and vertical stressstrain curves for five aggregate sizes (4.5, 6.0, 8.0, 10.0 and 12.0 mm.), using Voronoi with internal fracture path tessellation: A) Horizontal stressstrain response at center of specimen and B) Vertical stressstrain response at center of specimen. 10.00 8.00 6.00 4.00 2.00 0.00 *Agg. 4.5 mm  Agg. 6.0 mm  Agg. 8.0 mm *Agg. 10.0 mm xAgg. 12.0 mm  Experimental 7" 4 5 6 7 8 9 10 11 12 Aggregate Size (mm) Figure 48. Ultimate tensile strength averaged from 3 samples for each aggregate size (4.5, 6.0, 8.0, 10.0 and 12.0 mm). 0.020 0.015 0.010 0.005 0 .0 0 0 1 1 1 1 11 1 1 . 4 5 6 7 8 9 10 11 12 Aggregate Size (mm) Figure 49. Variance of residual error of tensile stress strain curve for each aggregate size (4.5, 6.0, 8.0, 10.0 and 12.0 mm). 2.75 2.76 2.81 2.89 2.59 0 4.6 Sensitivity Analysis of Material Parameters Defined for Mastics 4.6.1 Numerical Experiment This part of the study focused on the determination of the sensitivity of the parameters used to describe mastic strength. Voronoi tessellations (without internal fracture paths) that have uniform particle size of 6.8 mm were used with parallel crack growth rule in this experiment. First the numerical model was calibrated to match the stress strain curves of the Nova Scotia mixture to obtain material parameters for mastics. Then each of those calibrated parameters was altered to the higher and lower values and the resulting stressstrain response was observed. In particular, for comparison purposes, the ultimate strength was monitored and used in the comparison. The details of the numerical models and material parameters used for sensitivity analysis of the mastic are listed in Table 49 through Table 412. Table 49. Modeling matrix for the parametric study of mastic Component Description Tessellation Voronoi Crack growth rule Parallel Equivalent aggregate size (mm) 6.8 Exterior boundary elements 152 Potential crack elements for mastic 1408 Potential crack elements for fracture path in aggregate 0 Table 410. Global linear elasticity parameters for the parametric study of mastic Parameter Value Young modulus (MPa) 7500 Poisson ratio 0.35 Table 411. Local material parameters (calibrated material parameters) for the parametric study of mastic Potential crack To DNCR Tsoft Co CR Csoft (o iR element (MPa) (mm) (MPa/mm) (MPa) (MPa) (Mpa/mm) (degree) (degree) Mastic 2.30 0.11 20.91 2.40 0.15 5.00 38 30 Table 412. Varying parameters from the calibrated mastic shown in Table 411 for the parametric study of mastic To DNCR Tsoft Co CR Csoft (o (R (MPa) (mm) Linear Nonlinear (MPa) (MPa) (Mpa/mm) (degree) (degree) 1.90 0.03 76.67 1.00 2.00 0.000 1.00 32 24 2.10 0.07 32.86 2.20 0.075 3.00 35 27 2.30 0.11 20.91 1.00 2.40 0.150 5.00 38 30 2.50 0.15 15.33 2.60 0.225 7.00 41 33 2.70 0.19 12.11 1.00 2.80 0.300 9.00 44 36 4.6.2 Numerical Results Figure 410 shows typical results of horizontal and vertical stress strain responses with variation in the cohesion (Co). The subsequent Figures 411 through 418 show only the variation in ultimate tensile strength as the following parameters are varied: cohesion (Co), tension cutoff (To), friction angle (4o), cohesion softening slope (Csoft), tension softening slope (Tsoft), opening crack limit (DNCR), residual cohesion (CR) and residual friction angle (4R). It is noted that Tsoft and DNCR are redundant parameters. In DIGS, two types of tension softening slope, linear (L) and nonlinear (N), can be selected. The default linear model must set Tsoft equal to To/DNCR. Without doing so, the nonlinear (power law) will be chosen in which the tensile strength decreases in a nonlinear fashion from the selected tension cutoff (To) to zero at a specified opening crack limit (DNCR). All these parameters are used to define Mohr Coulomb type of failure in Figure 32. As expected, the ultimate tensile strength of specimen is strongly affected by the parameters controlling strength of failure envelope Co, To and 4o as shown in Figure 411 through 413). It was noticeable from Figure 411 and 412 that cohesion was slightly more sensitive than the tension cutoff. The parameters that control the rate of softening Csoft and Tsoft also affected the strength of specimen as can be seen in Figure 414 and 415. The opening crack limit (Figure 416) has an influence to the specimen strength in the same way as the tension softening slope. The last two parameters CR and 4R (Figure 417 and 418) controlling residual strength of failure envelope were found to be less sensitive to the specimen strength. Figure 419 shows the results of stress strain curves predicted by two tension softening models, Linear and Nonlinear. There were not significantly different between these two, which may imply that the approximation of rate of tension softening were not critical for simulating the IDT strength test. 5 4Cohesion 2.0 3.0 MPa .0 4Cohesion 2.2 2.5 MPa .0 A Cohesion 2.4 S2.0 MPa MPa 1.5 *Cohesion 2.6 MPa 0 1.0 *Cohesion 2.8 50 MPa I 0.5 .5 Experimental 0 .0 .. . 0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 Horizontal Strain (x 1000) A 10.0 4Cohesion 2.0 MPa M 8.0  Cohesion 2.2 ( MPa 6.0 aCohesion 2.4 MPa C 4Cohesion 2.6 S4.0 MPa *Cohesion 2.8 > 2.0 MPa Experimental 0.0 0.0 1.0 2.0 3.0 4.0 5.0 6.0 Vertical Strain (x 1000) B Figure 410. Change of stressstrain responses due to the varying of cohesion defined in mastic from 2.00 2.80 mm: A) Horizontal stressstrain response at center of specimen and B) Vertical stressstrain response at center of specimen. 3.50 3.00 2.50 2.00 1.50 1.00 0.50 0.00 1.80 2.00 2.20 2.40 2.60 2.80 3.00 Cohesion (MPa) Figure 411. Sensitivity of the parameter cohesion defined in mastic to the ultimate tensile strength. 3.50 3.00 2.50 2.00 1.50 1.00 0.50 0.00 1.80 2.00 2.20 2.40 2.60 2.80 Tension Cut Off (MPa) Figure 412. Sensitivity of the parameter tension cutoff defined in mastic to the ultimate tensile strength. 2. 8 2.)9 2.94 2.71 2.67 2.75 2.76 3.50 3.00 2.50 2.00 1.50 1.00 0.50 ' 0 .00 ........ 30 32 34 36 38 40 42 44 46 Friction Angle (degree) Figure 413. Sensitivity of the parameter friction angle defined in mastic to the ultimate tensile strength. 3.50 3.00 2.50 2.00 1.50 1.00 0.50 0.00 10.0 Cohesion Softening Slope (MPa/mm) Figure 414. Sensitivity of the parameter cohesion softening slope defined in mastics to the ultimate tensile strength. 2.90 2.75 2.75 2.73 2.64 4 2.88 S2.70 2 75 2.88 2.55 3.50 3.00 2.50 2.00 1.50 1.00 0.50 0.00 20.0 40.0 60.0 Tension Softening Slope (MPa/mm) 80.0 Figure 415. Sensitivity of the parameter tension softening slope (Linear Model) defined in mastic to the ultimate tensile strength. 3.50 3.00 2.50 2.00 1.50 1.00 0.50 0.00 0. 00 0.05 0.10 DNCR (mm) 0.15 0.20 Figure 416. Sensitivity of the parameter opening crack limit defined in mastic to the ultimate tensile strength 3.01 A2.76 22.75 235 2. 3.01 970 2.75 2.76 2.35 ow 0..)  3.0 2.74 2.75 2.80 2.85 S2.5 2.0 1.5 1.0  0.5 S 0 .0 . . . 0.00 0.05 0.10 0.15 0.20 0.25 0.30 Residual Conhesion (MPa) Figure 417. Sensitivity of the parameter residual cohesion defined in mastic to the ultimate tensile strength. 3.50 3.00 2.50 2.00 1 50 S1.00 0.50 5 0.00 22 24 26 28 30 32 34 36 38 Residual Friction Angle (degree) Figure 418. Sensitivity of the parameter residual friction angle defined in mastic to the ultimate tensile strength. 2.3 2.74 2.75 2.78 2.77 :  *    3.50 3.50 *DNCR 0.03 Z 3.00 _r (LINEAR) 2.50 DNCR 0.03 2.50 u, (NONLINEAR) S2.00 DNCR 0.11 f___ (LINEAR) 1.50 S E * DNCR 0.11 S1.00 (NONLINEAR) N 0.5 *DNCR 0.19 o 0.50 I (LINEAR) 0.00 ( DNCR 0.19 0.0 1.0 2.0 3.0 4.0 5.0 6.0 (NONLINEAR) Experimental Horizontal Strain (x 1000) A 10.0 DNCR 0.03 8. (LINEAR) 8.0 *0 DNCR 0.03 (NONLINEAR) S6.0 I DNCR 0.11 C (LINEAR) wt 4.0 A DNCR 0.11 E (NONLINEAR) E 2.0 *DNCR 0.19 > (LINEAR) 0.0 __._. __ 0 DNCR0.19 0.0 1.0 2.0 3.0 4.0 5.0 6.0 (NONLINEAR) Experimental Vertical Strain (x 1000) B Figure 419. A comparison of stressstrain response between using a Linear (default) and Nonlinear (still in experiment) tension softening model for the range of opening crack limit DNCR 0.03 0.19 mm: A) Horizontal stressstrain response at center of specimen and B) Vertical stressstrain response at center of specimen. 4.7 Sensitivity Analysis of Material Parameters Defined for Internal Fracture Paths 4.7.1 Numerical Experiment This part of the study focused on the identification of the sensitivity of the predicted ultimate tensile strength to the changes in material strength parameters defined for internal aggregate fracture paths. In the experiment, the same Voronoi tessellation scheme, previously used, with a particle size of 6.8 mm was filled with internal fracture paths to simulate faults in the aggregates. The same parallel crack growth rule was also used throughout, and the numerical model was first calibrated using the experimental stress strain curves from the Nova Scotia mixture. Material parameters for the mastic defined in Table 411 were kept the same from the previous study while calibrating material parameters for internal fracture paths of aggregate. Then each of the internal fracture path strength parameter was changed to higher or lower values from the calibrated parameters to observe the sensitivity of the predicted stressstrain responses. All details of the models and material parameters used for the sensitivity analysis of the internal fracture paths are listed in Table 413 though Table 416. 4.7.2 Numerical Results Figure 420 shows the effects of cohesion (Co) and tension cutoff (To) on the predicted ultimate strength. It can be seen that an increase of cohesion and tension cutoff from 4.5 to 5.5 MPa increases the strength of the specimen. Later increasing of both parameters to 8.5 MPa does not improve the strength. Figure 421 shows the effect of the opening crack limit (DNCR) on the ultimate tensile strength. Too small DNCR 0.01 mm cause brittle failure of specimen and higher DNCR 0.04, 0.07, 0.10, and 0.13 mm improves the strength. Figures 422 through 425 show the effects of other parameters, friction angle (4o), cohesionsoftening slope (Csoft), residual cohesion (CR), and residual friction angle (4R) on the ultimate tensile strength. The ultimate tensile strength was found to be insensitive to these parameters. 4.8 Summary of Sensitivity Analysis Numerical experiments of material parameters for the mastic and internal fracture paths revealed that the tensile strength of an IDT specimen was preliminary controlled by the mastic, which was the weakest link in specimen. The tensile strength of specimen was most sensitive to the parameters controlling strength (cohesion Co, friction angle 4o and tension cutoff To). Parameters that control rate of softening (cohesion softening slope Csoft, tension softening slope Tsoft and opening crack limit DNCR) also affect the strength of specimen. The residual parameters (residual cohesion CR and residual friction angle 4R) have the least effects on the strength of. Even though the last two residual parameters did not have a significant impact on strength, they do provide the smoothness for both ascending and descending curves of stress strain responses. When aggregate has internal fracture paths, the strength of the specimen is slightly weakened by the aggregate breakdown. In studies, low aggregate toughness defined by cohesion and tension cutoff of 4.5 MPa (or 1.9 times that of mastic strength) decreased the specimen strength. However, when aggregate is reasonable tough 6.0 to 8.5 MPa (or 2.5 to 3.5 times that of mastic strength), the toughness does not affect the strength of the specimen. Table 413. Modeling matrix for the parametric study of internal fracture paths Component Description Tessellation Voronoi with internal fracture path Crack growth rule Parallel Equivalent aggregate size (mm) 6.8 Exterior boundary elements 152 Potential crack elements for mastic 1408 Potential crack elements for fracture path in aggregate 2717 Table 414. Global linear elasticity parameters for the parametric study of internal fracture paths Parameter Value Young modulus (MPa) 7500 Poisson ratio 0.35 Table 415. Local material parameters (calibrated material parameters) for the parametric study of internal fracture paths Potential crack To DNCR Tsoft Co CR Csoft io (R element (MPa) (mm) (MPa/mm) (MPa) (MPa) (Mpa/mm) (degree) (degree) Mastic 2.30 0.11 20.91 2.40 0.15 5.00 38 30 Internal fracture path 6.50 0.07 92.86 6.50 0.15 1.00 40 32 Table 416. Varying parameters from the calibrated fracture path shown in Table 415 for the parametric study of internal fracture paths To DNCR CO CR Csoft 0o 4R (MPa) (mm) (MPa) (MPa) (Mpa/mm) (degree) (degree) 4.50 0.01 4.50 0.000 0.10 35.0 27.0 5.50 0.04 5.50 0.075 0.50 37.5 29.5 6.50 0.07 6.50 0.150 1.00 40.0 32.0 7.50 0.10 7.50 0.225 5.00 42.5 34.5 8.50 0.13 8.50 0.300 10.00 45.0 37.0 3.50 3.00 2.50 2.00 1.50 1.00 0.50 0.00 40.0 50.0 60.0 70.0 80.0 90.0 Cohesion and Tension (MPa) Figure 420. Sensitivity of both parameters cohesion and tension cutoff defined in fracture path of aggregate to the ultimate tensile strength. 3.50 a S3.00 S2.50 C 2.00 1.50  1.00 E 0.50 : 0.00 ^ ,n q) 7n 0.0 2.0 4.0 6.0 8.0 10.0 12.0 14.0 DNCR (mm) Figure 421. Sensitivity of opening crack limit defined in internal fracture path of aggregate to the ultimate tensile strength. 242.64 2.70 2.70 2.73 2.48 2.57 2.1 . :L '^ ^^  ^ 3.50 3.00 2.50 2.00 1.50 1.00 0.50 0.00 34.0 36.0 38.0 40.0 42.0 44.0 46.0 Friction Angle (degree) Figure 422. Sensitivity of friction angle defined in internal fracture path of aggregate to the ultimate tensile strength. 3.50 3.00 2.50 2.00 1.50 1.00 0.50 0.00 0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 Cohesion Softening Slope (MPa/mm) Figure 423. Sensitivity of cohesion softening slope defined in internal fracture path of aggregate to the ultimate tensile strength. 2.56 2.68 2. 2.68 2.68 ____ 4* " vt* ~_ ~ 2.0 ______ S.7( 2.70 .O 2.6,  2. 70  17T _ ' 3.50 a 3.00 S2.50 " 2.00 1.50  1.00 " 0.50  0.00 0.00 0.05 0.10 0.15 0.20 0.25 0.30 Residual Cohesion (MPa) Figure 424. Sensitivity of residual cohesion defined in internal fracture path of aggregate to the ultimate tensile strength. S3.50 a 3.00 S2.50 " 2.00 1.50  1.00  0.50 5 0.00 270 2.70 2. 0 2.66 2.68 , ,=  26.0 28.0 30.0 32.0 34.0 36.0 38.0 Residual Friction Angle (degree) Figure 425. Sensitivity of residual friction angle defined in internal fracture path of aggregate to the ultimate tensile strength. .7 6 2.85 2.76 2. 2.68.85 ^^4  4.9 Recommendations for a Successful Numerical Simulation of IDT Strength Test This section summarizes and proposes a standard procedure for successfully model Superpave IDT strength test with the displacement discontinuity boundary element method. * Select "reasonable" IDT test results. The IDT strength test that yields extremely high or low values of Poisson ratio (i.e. less than 0.1 or greater than 0.4) may show a sign of nonisotropic material behavior and should be disregarded. * Select "reasonable" loaddeformation test results from the Superpave IDT test. A large discrepancy of stress strain curves obtained from each face of an IDT sample may reveal unsymmetrical aggregate structure in a specimen or highly nonhomogeneous behavior at the location of strain gauges Numerical analysis starts with creating an IDT numerical model using exterior boundary elements for the edge of specimen and filling the inside with potential crack elements, thus forming Voronoi tessellations or Voronoi tessellations with internal fracture paths. With a suitable set of material parameters for the mastic and fracture paths in the aggregate, the predicted stress strain curves and crack patterns can match the experimental results. The steps for numerical simulation of the IDT strength test can be summarized as follows: 1. The first trial Young modulus (E) for numerical analysis can be approximated by a secant modulus at half ultimate strength Es = 1( va Subsequent trial values could be adjusted slightly higher or lower, such that it yields the best fit to the ascending portion of the vertical compressive and horizontal tensile stress strain curves. 2. The Poison ratio can be obtained directly from IDT strength test. Acceptable Poisson ratios for asphalt concrete should vary between 0.15 and 0.35. Extreme values such as less than 0.1 or greater than 0.4 should not be used. 3. Use Voronoi tessellation with parallel crack growth rule to calibrate material parameters for the mastic. 4. The particle size around D40 to D60 of gradation is recommended for consistent results. Using coarser aggregate size is also allowed to reduce computing time since the size effect is not very strong. 5. Approximate material parameters for the mastic: a. For the first trail, cohesion (Co) and tension cutoff (To) can be set to the ultimate tensile strength of a given sample. On subsequent trials, cohesion and tension cutoff may be adjusted slightly, according to the best fit to the vertical and horizontal stress strain curves. b. Friction angle (4,o) around 38 to 40 degrees and opening crack limit (DNCR) between 0.10 to 0.12 mm have been found to be reasonable values for asphalt mixture. c. A linear tension softening model Tsoft must be defined equal to To /DNCR. d. Other parameters for mastics such as residual friction angle (4R), residual cohesion (CR) and cohesionsoftening slope (Csoft) are chosen such that the predicted post peak responses fit the experimental post peak responses. The suggested values for the first trial are OR 30 degree, CR = 0.10 MPa and Csoft = 1.0 MPa/mm. 6. Replace Voronoi tessellation with its dual Voronoi tessellation with internal fracture paths and calibrate material parameters for fracture paths within the aggregates. 7. Material parameters for internal fracture paths within aggregates (in case the strength of the aggregate is not known, it is assumed that the aggregate is reasonable tough) a. Cohesion (Co) and tensile strength (To) can be set to any high values about 2.5 to 3.5 times of those used for mastic strength. b. Friction angle (4,o) of 40 degree and DNCR between 0.03 0.07 mm have been found to be good values for fracture paths inside aggregate. c. A linear tension softening slope Tsoft must be set to To/DNCR. d. Other parameters for internal fracture paths such as residual friction angle (4R), residual cohesion (CR) and cohesion softening slope (Csoft) are chosen such that the predicted post peak responses fit to the experimental post peak responses. The suggested values for the first trial are OR = 30 degree, CR 0.10 MPa and Csoft = 1.0 MPa/mm. 8. For known soft aggregates the exact values for the material parameters can be used. 9. Adjust material parameters until the predicted stress strain curves fit to the experimental stress strain curves for both the vertical and horizontal directions. 4.10 Summary From the numerical studies conducted in this Chapter, it can be concluded that the appropriate tessellation scheme to represent the aggregate structure in the Superpave IDT strength test was a Voronoi tessellation scheme with internal fracture paths. This tessellation scheme predicted the stress strain responses closest to the experimental results. The internal fracture paths allowed the aggregates to break down at high loads, which realistically simulate the cracking behavior of IDT specimen observed in laboratory. The study to identify the appropriate crack growth rule for the indirect tension strength test resulted in the recommendation of using the parallel crack growth rule. Both parallel and sequential methods predicted approximately the same stress strain responses and crack patterns, but the parallel consumed much less computing time than the sequential did. Furthermore, it is highly possibly that several cracks in a specimen can occur simultaneously for fast loading test. The parametric studies of the aggregate size, material parameters for mastics, and material parameters for internal fracture path revealed that the primary factor that controlled the strength of IDT specimen were the strength of mastics. The secondary factor was the strength of aggregate. The least important factor was the aggregate size. Even though the aggregate size was not really important for the strength of specimen, from a calibration point of view, the larger the aggregate size the more difficult it is to tune the material parameters to fit the experimental results and the smaller the aggregate size, the longer the computing time. For the gradations studied, it was found that an 70 optimal aggregate size around 4060% passing of gradation (D40 D60) yields the most consistent results with a reasonable amount of computational time. CHAPTER 5 EVALUATION OF TENSILE STRENGTH AND FRACTURE ENERGY DENSITY 5.1 Overview Tensile strength at fracture and fracture energy density have identified in previous studies (Sedwick, 1998 and Garcia, 2002) to be related to the fracture resistance of mixtures. The Superpave tensile strength is referred to as the strength at the fracture point (Roque (1997) and Buttlar et al.,(1997)). The strength at this point is normally slightly lower than the tensile strength at ultimate load. Another parameter, the fracture energy density at the center of the IDT specimen, has also been used to evaluate the cracking resistance of mixtures. The fracture energy can be obtained by integrating the area under tensile stressstrain curve up to the point of fracture. During Superpave IDT strength test, load is recorded as well as vertical and horizontal deformations at strain gauges. To account for threedimensional effects, bulging correction factors Cbh and Cbv shown in Table 51 are applied to correct the measured horizontal and vertical deformations to the deformations in flat planes Hcorrected = Hmea,, edCbh (51) corrected measured Cb (52) Table 51. Correction factors accounting for bulging effect D = 4" or 6" Poisson's ratio Diameter to Thickness ratio (t/D) v 0.167 0.333 0.500 0.625 0.750 0.20 0.9816 0.9638 0.9461 0.9358 0.9294 CBX 0.35 0.9751 0.9518 0.9299 0.9179 0.9108 0.45 0.9722 0.9466 0.9234 0.9111 0.9040 0.20 0.9886 0.9748 0.9677 0.9674 0.9688 CBY 0.35 0.9808 0.9588 0.9479 0.9473 0.9493 0.45 0.9759 0.9492 0.9361 0.9358 0.9380 The corrected horizontal and vertical deformations are then divided by the gauge length GL to obtain the average strains, and again the strain center correction factors Ceh = 1.072 and Cev = 0.977 are applied to change the average strains to the strains at the center of specimen reacted H corrected Ch (53) GL S coected corrected Cev (54) GL In the same way as strains, the corresponding stresses in the horizontal and vertical directions evaluated at center of specimen can be computed by twodimensional plane stress formulas and are then corrected with the center stress correction factors Coh and Co shown in Table 52 to convert these stresses in a plane to the stresses on the surfaces of 3D specimen 2P Ch corrected 2 h (55) izDT 6P v corrected o, (56) TilDT Fracture point in specimen can be detected by monitoring the deformation differential (Vcorrected Hcorrected) recorded at each load step and marking the point when it starts to deviate from the smooth curve. Corresponding to the fracture point, the recorded applied load and the measured horizontal and vertical deformations are used to calculate the tensile strength and fracture energy density as previously described. Table 52. Correction factors accounting for stress at center of specimen D = 4"or 6" Poisson's ratio Diameter to Thickness ratio (t/D) v 0.167 0.333 0.500 0.625 0.750 CoxCTR 0.20 0.9471 0.9773 1.0251 1.0696 1.1040 0.35 0.9561 1.0007 1.0871 1.1682 1.2321 0.45 0.9597 1.0087 1.1213 1.2307 1.3171 CoyCTR 0.20 0.9648 0.9754 0.9743 0.9693 0.9611 0.35 0.9732 0.9888 0.9844 0.9710 0.9538 0.45 0.9788 0.9971 0.9864 0.9646 0.9395 5.2 Numerical Model V.S. Asphalt Mixtures The displacement discontinuity boundary element method assumes the material to be linear elastic until fracture. When cracks appear in the specimen, the linear elastic behavior still remains but with a plastic failure and localization in the crack. For this part of the study, three mixtures, Nova Scotia, NW39_1C and I95SJNBWP described in Table 53 were chosen. These mixtures range from having a predominantly ductile failure mode (Nova Scotia mix) to a very brittle failure mode (I95SJN BWP) to evaluate the capabilities of the displacement discontinuity method of simulating very different mixture responses. The material properties obtained from IDT strength tests of three mixtures are listed in Table 54. In this study, a single "representative" sample was selected out of the three samples tested for each mixture, based on the closeness of the measured loaddeformation responses at both sides of the sample. In the case of the Nova Scotia and NW39_1C mixtures, the middle sample was selected, but for the 195SJNBWP the upper bound sample was selected. Once a "representative" sample for each mix was selected, the measured load and deformations were used to calculate vertical and horizontal stress strain curves on the faces of the specimen as described in the previous section. Then the averaged stress straincurves from both faces were used for comparisons to the numerical simulation. Table 53. Gradation of three mixtures used in simulations Sieve size Percent passing Nova Scotia NW39 1C 195SJN BWP 25 mm (1) 100 19 mm (3/4) 97 100 100 12.5 mm (1/2) 83 93 91 9.5 mm (3/8) 66 87 76 4.75 mm (#4) 38 62 53 2.36 mm (#8) 23 46 40 1.18 mm (#16) 18 39 36 600 tm (#30) 14 33 34 300 tm (#50) 10 27 28 150 tm (#100) 6.5 14 19 75 tlm (#200) 3.5 8 5 (1) Nova Scotia was a virgin coarsegraded mix ni ilh 19 mm nominal maximum aggregate size, designed according to SuperPaveT specifications. The mixture has 7. 0% air void and 5.5% design asphalt c "ine n i/th an asphalt binder of grade AC30. (2) NW391C was a core excavated on Section 1 (crackedpavement), from Alachua County, FL, on NW39th Ave., eastbound, section limit 34 St 24 St. The pavement was 14 years old at the time of coring. The mixture was afinegraded mix n i/lh 12.5 mm nominal maximum aggregate size, 7.3% air void, 4.7% effective asphalt content (the recovered unknown binder has Brookfield viscosity 16000 Poises at 60 C). (3) I95SJN BWP was a core excavated between wheel paths, from St. Johns County on Interstate 195, northbound, county milepost limit 26.21834.855. The age of the pavement was 5 years old. The mixture was afinegraded mix n ilh 12.5 mm nominal maximum aggregate size, 4.0% air voids, 4. 7% effective asphalt content (AC30 was believed to be used i1 ith the recovery Brookfield viscosity 36548 Poises at 60 C). Table 54. Material properties of three mixtures obtained from Superpave IDT tests Mixtures Test Failure Tensile Fracture Poisson Do D1 m temp. strain strength energy ratio (micro density (1/Gpa* (C) strain) (MPa) (KJ/m3) (1/GPa) sec) Nova Scotia 0 618 2.05 0.83 0.35 0.0645 0.0271 0.567 NW39 1C 10 786 2.70 1.38 0.30 0.0521 0.0192 0.524 I95SJN BWP 10 590 3.19 1.44 0.22 0.0936 0.0014 0.574 5.3 Numerical Simulation of Three Mixtures Modeling matrix of the displacement discontinuity models for simulating three mixtures, Nova Scotia, NW39 1C and I95SJN BWP are shown in Table 55. Based on the size effect study of Voronoi particles in Chapter 4, Voronoi particles with average diameters of 8.1 mm were used for the Nova Scotia coarse mix, and slightly smaller particles 6.8 mm for NW39_lC fine mix and 6.6 mm for 195SJNBWP finegraded mix. The numbers of exterior boundary elements for constructing the edge of three specimens were 152, of those, 8 elements each were used for simulating fix condition at the bottom and displacement control at the top. The displacement control for all three mixtures was divided into 20 equal load steps for the total vertical displacement of 0.35, 0.30 and 0.19 respectively. The number for potential crack elements, simulating mastics and internal fracture paths inside Voronoi particles are provided in Table 55. Table 56 shows global material parameters (Young modulus and Poisson ratio) for numerical models simulating three mixtures. Poisson ratios were obtained directly from IDT strength test. The secant modulus at 50 percent of the ultimate load was found to be a good value for the Young modulus in the numerical model. Other local material parameters for mastics and internal fracture paths for the three numerical models in Table 57 were determined by numerical tuning. From the sensitivity analysis presented in Chapter 4, it was found that cohesion Co and tension cutoff To of the mastic control the tensile strength of the specimen. Therefore, by tuning these two parameters, the predicted ultimate load in the stressstrain curves can fit to the ultimate load of given stressstrain curves. The residual strength of the mastic was also used to provide additional strength when the material damages. They provided the fit to the stressstrain responses in the post peak region of the stressstrain curves. Other local parameters for the mastic and fracture path in the aggregates are less important; they can be used to adjust the predicted stress strain curves to have the best fit of given stressstrain curves. The complete calibrated local material parameters that make numerical models fit to the test results are listed in Table 57. Once the numerical models fit the test results, the tensile strength at fracture and fracture energy density were evaluated according to SHRP method described earlier. Table 55. Modeling matrix for three mixtures used: Mixtures Specimen Equivalent Number Number Number of Total size Voronoi of exterior of potential vertical (diameter x diameter boundary potential crack deformation thickness) elements crack elements for at load plate elements for mastic internal for 20 load fracture step (mm x mm) (mm) paths (mm) Nova Scotia 150.0 x 50.3 8.1 152 978 1,791 0.350 NW39 1C 150.6 x 28.2 6.8 152 1,409 2,651 0.325 I95SJN BWP 143.5 x 53.6 6.6 152 1,249 2,287 0.190 Table 56. Global linear elasticity parameters for three mixtures used Mixtures Young modulus Poisson ratio (MPa) Nova Scotia 7,250 0.35 NW39 1C 8,500 0.30 I95SJN BWP 13,250 0.22 Table 57. Local material parameters for three mixtures used Mixtures Potential To DNCR Tsoft Co CR Csoft 0o (R crack (MPa/ (Mpa/ elements (MPa) (mm) mm) (MPa) (MPa) mm) (degree) (degree) Mastic 2.30 0.10 1.0 2.50 0.15 1.0 38 30 Nova Scotia Internal fracture path 6.25 0.05 1.0 6.25 0.10 1.0 40 30 Mastic 3.05 0.07 1.0 3.20 0.10 1.0 35 20 NW39 1C Internal fracture path 8.00 0.05 1.0 8.00 0.10 1.0 40 30 Mastic 4.50 0.03 1.0 4.90 0.10 1.0 35 20 195SJN BWP Internal fracture path 9.00 0.03 1.0 9.00 0.10 1.0 40 30 5.4 Numerical Test Results The predicted stress strain curves of three mixtures (Figure 51) agreed well with the experimental results. The predictions matched both horizontal tensile and vertical compressive stress strain curves up to the peak and they partially captured the post peak responses. Figure 52 shows the deformation differential of three mixtures: Nova Scotia, NW39_1C and I95SJNBWP. The fracture points marked in Figure 51A were determined at the peaks of the curves where the deformation differentials start to drop. The corresponding tensile strength at fracture point is shown in Figure 53. The differences between the predicted values and the average values for Nova Scotia, NW39_1C and I95SJN BWP were 5.3%. 5.6% and 12.5% respectively. Figure 54 compares the predicted fracture energy densities to the trim mean average values of the mixes (Table 54). The trim mean average value of a given mix was calculated by excluding the highest and the lowest of the reported values on 6 faces of the 3 samples and average the remaining four. The differences between the predicted values and the trim mean average values for Nova Scotia, NW39_1C and 195SJNBWP were 3.6%, 0.7% and 31.2% respectively. The large discrepancy of the 195SJNBWP was likely primary due to the variability in the experimental test results. The test result of mix 195SJNBWP samples 1 to 3 had tensile strengths of 3.75, 3.38 and 2.45 MPa respectively and fracture energy densities reported on each face of sample 1 to 3 were (1.56, 1.68), (0.70, 1.33), and (1.20, 2.38) KJ/m3 respectively. Overall, the method yields a good approximation for tensile strength at fracture and fracture energy density when the material variability is small and becomes less reliable when the variability is high. Figure 55 shows crack patterns of three mixtures at three critical load steps: A) crack initiation, B) crack at fracture points and C) major opening cracks at final load step. All three samples have the same crack pattern that several micro cracks randomly initiate in the tension zone and coalesce to form a major opening crack, which can be observed in laboratory. Nova Scotia (Experiment) e Nova Scotia __(Experiment) NW39_1C (DDM) . 195SJN_BWP (Experiment) 195SJN BWP (DDM) 0 1 2 3 4 5 6 K Fracture Point Horizontal Strain (x 1000) A Nova Scotia (Experiment) Nova Scotia (DDM) NW39 1C (Experiment) NW39_1C (DDM)  195SJN BWP (Experiment) 195SJN_BWP (DDM) 0 1 2 3 4 5 6 Vertical Strain (x 1000) B Figure 51. Comparison of predicted and measured horizontal and vertical stressstrain curves for three mixtures: Nova Scotia, NW39 IC and 195SJN BWP. A) Horizontal stress strain responses at center of specimen and B) Vertical stress strain responses at center of specimen.  Nova Scotia 15 (DDM) " NW39_1 C (DDM) S10 10 A195SJN_B I WP (DDM) 5 )K Fracture Point 0 1 0 5 10 15 20 Load Step Figure 52. Simulated deformation differentials for the Nova Scotia, NW39_1C and 195SJN BWP mixtures. 4.0 3.59 3.19 2.85 a. 3.0 3.0 ^2.70 5 2.16 S2.2.005 1.0 SIDDM O DDM Nova Scotia NW39 1 C 195SJN BWP Figure 53. Tensile strength at fracture for three mixtures (Nova Scotia, NW39_1C and 195SJN BWP). The average value of each mixture was compared to the tensile strength predicted by displacement discontinuity method (DDM). 0.99 EI Trim Mean Average ODDM Nova Scotia NW39 1 C 195SJN BWP Figure 54. Fracture energy density at center of specimen for three mixtures (Nova Scotia, NW39_1C and 195 SJNBWP). The trim mean average value of each mixture was compared to the fracture energy density predicted by displacement discontinuity method (DDM). 1.44 1.38 1.39 r 0.83 0.80 Nova Scotia Load Sten 7 NW39 1C Load Sten 9 I95SJN BWP Load Sten 14 Load Sten 14 Load Sten 18 Load Sten 20 Load Sten 20 Load Sten 20 Figure 55. Predicted crack patterns for the three mixtures (Nova Scotia, NW39_1C and 195SJN BWP) at three load steps: A) Crack initiation in specimen (first cracks that appear in specimen), B) Crack pattern at fracture point and C) Cracks at the final load step. 83 5.5 Conclusions Explicit fracture modeling with the displacement discontinuity boundary element method successfully simulated the cracking behavior of three asphalt mixtures and evaluated their tensile strength and fracture energy density. It can be concluded, according to three mixtures demonstrated here, that the explicit fracture model with displacement discontinuity elements has the potential to simulate the cracking behavior and reasonably evaluate tensile strength and fracture energy density. CHAPTER 6 SUMMARY, CONCLUSIONS AND RECOMMENDATIONS 6.1 Summary Several studies using the displacement discontinuity boundary element method were conducted to identify the most appropriate set of components for a numerical model that can realistically simulate the Superpave IDT strength test. Parametric studies were also performed to identify the key factors that control the cracking behavior of asphalt mixtures. Finally, the method was applied to evaluate the fracture behavior of mixtures. A summary of the findings obtained from each study is listed below: * The most suitable tessellation scheme to represent aggregate structure was the Voronoi with internal fracture paths. Voronoi tessellations alone can also yield a reasonable results but it was found difficult to calibrate material parameters to fit the experimental stress strain curves, especially at the peak and in the postpeak part of the stressstrain curves. Deluanay tessellation schemes were found unsuitable to represent aggregate structure because the resulting cracked structure causes floating regions, which sometimes lead to instability at strain gauge locations. * The mastic properties are key to controlling the cracking behavior of asphalt mixtures. At the beginning, the initiation and propagation of microcracks in the mastic weakens the aggregate structure and causes compression softening in the stressstrain curve. As the number of active cracks and their sizes increase, the sample stiffness decreases rapidly. The strength is terminated when the random major cracks appear in the mastic. The material parameters defined for the mastic that are directly related to the strength of the specimen are cohesion, tension cutoff and friction angle. The rate of cohesion softening, tension softening and opening crack limit also affected the strength of specimen. The residual cohesion and residual friction angle had little effect to the strength of specimen, but they do control the post peak behavior of stress strain responses. * Toughness of rock represented by the strength of the fracture paths inside aggregates is not sensitive to the specimen strength as long as it is about 23 times larger than the mastic strength. The break down of aggregate at high load levels yields a smooth compression softening curve near the peak of the vertical compression stress strain curve. * In the numerical experiments, using different uniform aggregate size revealed a small size effect on the strength of specimen. Smaller aggregate sizes yield more distributed cracks and cause softer responses than the larger aggregate sizes. * The current version of DIGS that uses MohrCoulomb type of failure can simulate the Superpave IDT strength test and evaluate its mixture properties with acceptable accuracy. 6.2 Conclusions Based on the finding previously discussed, the following conclusion were made: * The displacement discontinuity boundary element method with tessellation scheme can simulate the Superpave IDT strength test. * Voronoi tessellation with internal fracture paths is the best representative of aggregate structure. * The parallel crack growth rule is suitable for simulating the Superpave IDT strength test because it is more economical than the sequential crack growth rule. * The factors that control the strength of a Superpave IDT specimen ranked from the most to the least are strength of mastic, strength of fracture path in aggregates, and aggregate size. 6.3 Recommendations Based on the findings and conclusions of numerical experiments the following recommendations are made: * From the numerical studies, the mastic, which is the weakest link in the specimen, plays the most important role in cracking behavior in the Superpave IDT test. The mastic in a numerical model refers to the apparent strength of the mastic in the presence of air voids. The improvement of mastic properties will directly result in the increase of cracking resistance. * Despite the success of the numerical modeling, the simulations were merely the results of the assumptions made in the models. Whether the assumptions are reasonable may require the development of testing procedures to more directly obtain values for the micromechanical parameters used. * Other failure laws that can simulate viscoplastic and cyclic loading should be developed so that the method will yield more useful applications. 86 * The incorporation of multiregion boundary element codes into DIGS that can simulate layered pavement is recommend for future research. These will provide analytical tools to simulate and observe cracking behavior of specific mixtures in a pavement. 