Prediction of Damage Behaviors in Asphalt Materials Using a Micromechanical Finite-Element Model and Image Analysis Qingli Dai1; Martin H. Sadd2; Venkit Parameswaran3; and Arun Shukla4 Abstract: A study of the micromechanical damage behavior of asphalt concrete is presented. Asphalt concrete is composed of aggregates, mastic cement, and air voids, and its load carrying behavior is strongly related to the local microstructural load transfer between aggregate particles. Numerical simulation of this micromechanical behavior was accomplished by using a finite-element model that incorporated the mechanical load-carrying response between aggregates. The finite-element scheme used a network of special frame elements each with a stiffness matrix developed from an approximate elasticity solution of the stress and displacement field in a cementation layer between particle pairs. Continuum damage mechanics was then incorporated within this solution, leading to the construction of a microdamage model capable of predicting typical global inelastic behavior found in asphalt materials. Using image processing and aggregate fitting techniques, simulation models of indirect tension, and compression samples were generated from surface photographic data of actual laboratory specimens. Model simulation results of the overall sample behavior and evolving microfailure/ fracture patterns compared favorably with experimental data collected on these samples. DOI: 10.1061/共ASCE兲0733-9399共2005兲131:7共668兲 CE Database subject headings: Damage; Asphalt concrete; Finite element method; Image analysis.

Introduction Asphalt concrete is a complex heterogeneous material composed of aggregate, mastic cement, additives, and void space. For such materials, the macroload carrying behavior and resulting failure depend on many microphenomena that occur at the aggregate/ mastic level. Important microbehaviors are related to mastic properties including volume percentage, elastic moduli, inelastic/timedependent response, aging hardening, microcracking, and debonding from aggregates. Other microstructural features include aggregate size, shape, texture, and packing geometry. Because of these issues it appears that a micromechanical model is needed to simulate this material. Micromechanics offers the possibility to more accurately predict asphalt failure behavior and to relate such behavior to particular mix parameters such as mastic properties, aggregate gradation, and sample compaction. Recently, many studies have been investigating the micromechanical behavior of particulate, porous, and heterogeneous materials. For example, studies on cemented particulate materials by Dvorkin et al. 共1994兲 and Zhu et al. 共1996兲 provided information 1 Research Associate, College of Engineering, Texas A&M Univ.-Kingsville, Kingsville, TX 78363. 2 Professor, Dept. of Mechanical Engineering and Applied Mechanics, Univ. of Rhode Island, Kingston, RI 02881. 3 Visting Professor, Dept. of Mechanical Engineering and Applied Mechanics, Univ. of Rhode Island, Kingston, RI 02881. 4 Professor, Dept. of Mechancial Engineering and Applied Mechanics, Univ. of Rhode Island, Kingston, RI 02881. Note. Associate Editor: Arif Masud. Discussion open until December 1, 2005. Separate discussions must be submitted for individual papers. To extend the closing date by one month, a written request must be filed with the ASCE Managing Editor. The manuscript for this paper was submitted for review and possible publication on March 24, 2003; approved on December 7, 2004. This paper is part of the Journal of Engineering Mechanics, Vol. 131, No. 7, July 1, 2005. ©ASCE, ISSN 0733-9399/ 2005/7-668–677/$25.00.

668 / JOURNAL OF ENGINEERING MECHANICS © ASCE / JULY 2005

on the normal and tangential load transfer between cemented particles, and would be fundamental in developing a micromechanical theory for load distribution and failure of such materials. Recent applications of such contact-based micromechanical analysis for asphalt concrete behavior have been reported by Chang and Gao 共1997兲, Cheung et al. 共1999兲, and Zhu and Nodes 共2000兲. In a related study, Krishnan and Rao 共2000兲 used mixture theory to explain air void reduction in asphalt materials under load. Numerical modeling of cemented particulate materials has generally used both finite and discrete-element methods. The discrete element method 共DEM兲 analyzes particulate systems by modeling the translational and rotational behaviors of each particle using Newton’s second law. Using appropriate interparticle contact forces, the scheme establishes a time-stepping procedure to determine each of the particle motions. The DEM studies on cemented particulate materials include the work by Rothenburg et al. 共1992兲, Chang and Meegoda 共1993兲, Trent and Margolin 共1994兲, Buttlar and You 共2001兲, Ullidtz 共2001兲, and Sadd et al. 共1992兲, and Sadd and Gao 共1997兲. In regard to finite element modeling 共FEM兲, Sepehr et al. 共1994兲 used an idealized finite element microstructural model to analyze the behavior of an asphalt pavement layer. Soares et al. 共2003兲 used cohesive zone elements to develop micromechanical fracture model of asphalt materials. A common finite element approach to simulate particulate and heterogeneous materials has used an equivalent lattice network system to represent the interparticle load transfer behavior. This type of microstructural modeling has been used previously: Bazant et al. 共1990兲, Mora 共1992兲, Sadd and Gao 共1998兲, and Budhu et al. 共1997兲. Along similar lines, Guddati et al. 共2002兲 recently presented a random truss lattice model to simulate microdamage in asphalt concrete and demonstrated some interesting failure patterns in an indirect tension test geometry. Bahia et al. 共1999兲 have also used finite elements to model the aggregate-mastic response of asphalt materials. Mustoe and Griffiths 共1998兲 developed a FEM, which was

equivalent to a particular discrete element approach. They pointed out that the FEM model has an advantage over the discrete element scheme for static problems. Damage mechanics provides a viable framework for the description of distributed material damage including material stiffness degradation, microcrack initiation, growth and coalescence, and damage-induced anisotropy. Continuum damage mechanics is based on the thermodynamics of irreversible processes to characterize elastic-coupled damage behaviors 共Chaboche 1988兲. Simo and Ju 共1987兲 developed strain- and stress-based anisotropic continuum damage models, and Kachanov 共1987兲 proposed a microcrack-related continuum damage model for brittle materials. These models focus on the relation between damage and effective elastic properties. With respect to asphalt materials, Lee and Kim 共1998兲 developed a viscoelastic constitutive model to study the rate-dependent damage growth and damage healing behaviors. Sangpetngam et al. 共2003兲 recently used a displacement discontinuity boundary element approach to simulate the cracking behavior of asphalt mixtures. Papagiannakis 共2002兲 employed imaging techniques to capture the asphalt concrete microstructure and conducted FEM studies for the viscoelastic response. Work by Wang et al. 共1999兲 used imaging techniques for measuring strains within the mastic by comparing deformations before and after the load application. This paper presents a numerical modeling scheme for asphalt concrete based on micromechanical simulation using the finiteelement method. The model first incorporates an equivalent lattice network approach whereby the local interaction between neighboring particles is modeled with a special frame-type finite element. The element stiffness matrix is constructed by considering the normal, tangential and rotation behaviors between cemented particles using an approximate elasticity solution within the interparticle cementation. Inelastic mastic damage behaviors are then simulated by incorporating a continuum damage mechanics theory within the FEM model. This theoretical formulation was then implemented into the commercial ABAQUS finite-element analysis code using user-defined elements. In our previous work Sadd et al. 共2001, 2002, 2004a,b兲 applied these modeling procedures on numerically generated asphalt samples with idealized microstructure. The current work explores the use of real asphalt samples to generate the numerical model used in the finite element simulation. To capture the microstructure of real asphalt materials, numerical imaging models were generated from surface scans of actual asphalt specimens. Using image processing, digital photographs of the sample’s surface microstructure provided numerical information for curve fitting algorithms to determine idealized aggregate dimensions and locations. Imaging models of indirect tension and compression specimens were generated to conduct numerical simulations. Comparisons were made of the overall sample behavior and the internal damage pattern between the imaging model and test specimen.

Fig. 1. Multiphase asphalt materials

binder and fine particles兲 and air voids 共see Fig. 1兲. The load transfer between the aggregates plays a primary role in determining the load carrying capacity and failure of such complex materials. In order to develop a micromechanical model of this behavior, proper simulation of the load transfer between the aggregates must be accomplished. The aggregate material is normally much stiffer than the mastic, and thus aggregates are taken as rigid particles. On the other hand, the mastic cement is a compliant material with elastic, inelastic, and time-dependent behaviors. In order to properly account for the load transfer between aggregates, it is assumed that there is an effective mastic zone between neighboring particles. It is through this zone that the micromechanical load transfer occurs between each aggregate pair. This loading can be reduced to resultant normal and tangential forces and a moment as shown in Fig. 1. In order to model the interparticle load transfer behavior, some simplifying assumptions must be made about allowable aggregate shape and mastic geometry. Studies on aggregate geometry have commonly quantified particle size, shape, angularity and texture. However, for the present modeling only size and shape will be considered. In general, asphalt concrete contains aggregate of very irregular geometry as shown in Fig. 2共a兲. Our approach is to allow variable size and shape using an idealized elliptical aggregate model as represented in Fig. 2共b兲. The finite-element model then uses an equivalent lattice network approach, whereby the interparticle load transfer is simulated by a network of specially created frame-type finite elements connected at particle centers as shown in Fig. 2共c兲. From granular materials research, the material microstructure or fabric can be characterized to some extent by the distribution of branch vectors which are the line segments drawn from adjacent particle mass centers. Note that the proposed finite-element network coincides with the branch vector distribution. Cementation between neighboring particles was generated using two different schemes: freely filled and maximally filled

Micromechanical Finite-Element Model The basics of our micromechanical finite-element model have been developed in previous studies 共Sadd et al. 2001, 2002, 2004a,b兲. In order to present the new applications in the current work, we briefly review some of the basic model developments. Bituminous asphalt concrete can be described as a multiphase material containing aggregate, mastic cement 共including asphalt

Fig. 2. Asphalt modeling concept: 共a兲 typical asphalt material; 共b兲 model asphalt system; and 共c兲 network finite-element model JOURNAL OF ENGINEERING MECHANICS © ASCE / JULY 2005 / 669

Fig. 3. Cementation between two neighboring particles: 共a兲 freely filled cementation; and 共b兲 maximally filled cementation

cementation as shown in Fig. 3. Freely filled cementation is parallel to the branch vector and can be asymmetrically distributed 关Fig. 3共a兲兴. Maximally filled cementation is symmetrical about the branch vector and completely occupies the area between the two ellipse axes perpendicular to the branch vector 关Fig. 3共b兲兴. The cementation geometry parameters 共width w and average thickness ¯h兲 are calculated as shown in Fig. 3. The current network model uses a specially developed, twodimensional frame-type finite element to simulate the interparticle load transfer. These two-noded elements have the usual three degrees of freedom 共two displacements and a rotation兲 at each node and the element equation is written as

much greater than that of the cement layer, the particles are assumed to be rigid. Dvorkin has shown that effects of nonuniform cement thickness for each mastic element are generally small, so the analytical solution with average cement thickness ¯h for each mastic element was used. The two-dimensional model geometry 共uniform thickness case兲 is shown in Fig. 4. Note that freely filled scheme allows arbitrary nonsymmetric cementation, and thus the finite element line will not necessarily pass through the center of the mastic material. Thus in general, w = w1 + w2, but w1 ⫽ w2 ⫽ w / 2, and an eccentricity variable may be defined by e = 共w2 − w1兲 / 2. The stresses ␴x, ␴z, and ␶xz within the cementation layer can be calculated for unit normal, tangential, and rotational particle motion cases. These stresses can then be integrated to determine the total load transfer between particles through the cement mastic. For example, the resultant forces and moments on a pair of particles 共1兲 and 共2兲 for the cases of unit normal, tangential and rotational motions on particle 共1兲 are given by

F共1兲 n =

冕

w

兩␴z兩z=h¯dx =

0

冕

w

0

共␭ + 2␮兲 共␭ + 2␮兲w dx = ¯h ¯h

共1兲 F共2兲 n = − Fn

F共1兲 t =

冕

w

兩␶xz兩z=h¯dx = 0

0

共2兲

共1兲 F共2兲 t = Ft

冤

K11 K12 K13 . K22 K23 . . K33 . . . . . . . . .

K14 K24 K34 K44 . .

K15 K25 K35 K45 K55 .

K16 K26 K36 K46 K56 K66

冥冦 冧 冦 冧 U1 V1 ␪1 U2 V2 ␪2

=

Fn1 Ft1 M1 Fn2 Ft2 M2

M 共1兲 =

冕

w

兩␴z兩z=h¯共x − w1兲dx =

0

冕

w

0

共1兲 =

共␭ + 2␮兲 共x − w1兲dx ¯h

共␭ + 2␮兲w e ¯h

M 共2兲 = − M 共1兲

where Ui, Vi, and ␪i = nodal displacements and rotations; and F . . and M . = nodal forces and moments. The usual scheme of using bar and/or beam elements to determine the stiffness terms is not appropriate for the current applications, and therefore these terms were determined using an approximate elasticity solution from Dvorkin et al. 共1994兲 for the stress distribution in a cement layer between two particles. Since the particle material stiffness is

F共1兲 n =

冕

w

兩␴z兩z=h¯dx = 0

0

共1兲 F共2兲 n = Fn

F共1兲 t =

冕

w

兩␶xz兩z=h¯dx =

0

冕

w

0

␮ ␮ dx = w ¯h ¯h 共3兲

共1兲 F共2兲 t = − Ft

M 共1兲 =

冕

w

兩␶xz兩z=h¯dx · r1 =

0

Fig. 4. Cementation between two adjacent particles 670 / JOURNAL OF ENGINEERING MECHANICS © ASCE / JULY 2005

M 共2兲 =

␮w r2 ¯h

␮w r1 ¯h

F共1兲 n =

冕

w

兩␴z兩z=h¯dx =

0

冕

w

0

共␭ + 2␮兲 共␭ + 2␮兲w 共x − w1兲dx = e ¯h ¯h

共1兲 M 共1兲 = M rn + M rt共1兲 =

共1兲 F共2兲 n = − Fn

冕

w

兩␶xz兩z=h¯dx =

0

冕

M 共2兲 = − w

0

␮r1 ␮w dx = r1 ¯h ¯h 共4兲

共1兲 F共2兲 t = − Ft

关K兴 =

冤

Knn 0 · Ktt ·

·

· ·

· ·

·

·

冕

w

兩␶xz兩z=h¯r21dx

0

共␭ + 2␮兲w 2 ␮w 2 共w2 − w1w2 + w21兲 + r ¯ ¯h 1 3h

共␭ + 2␮兲w 2 ␮w 共w2 − w1w2 + w21兲 + r 1r 2 ¯ ¯h 3h

0

− Knne

− Ktt

Kttr2

Knn 2 Knn 2 共w2 − w1w2 + w21兲 − Knne − Kttr1 Kttr1r2 − 共w2 − w1w2 + w21兲 3 3 · 0 Knn Knne · · − Kttr2 Ktt Knn 2 · · · 共w2 − w1w2 + w21兲 Kttr22 + 3

where Knn = 共␭ + 2␮兲w / ¯h; Ktt = ␮w / ¯h; ␭ and ␮ = usual elastic moduli; w and ¯h = the cementation width and average cementation thickness; r1 and r2 = radial dimensions from each aggregate center to the cementation boundary; w1 and w2 = left and right width of cementation; and e = 共w2 − w1兲 / 2. Each mastic element stiffness matrix is significantly different depending on the two-particle layout and size, and mastic geometry. This procedure establishes the elastic stiffness matrix, which is a function of the material microstructure and mastic moduli.

冥

␴ = ␴0共1 − e−b共␧/␧0兲兲 ⇒

⳵␴ = D0e−b共␧/␧0兲 ⳵␧

共6兲

In order to characterize the nonlinear damage behavior of asphalt concrete, a particular Weibull distribution function is chosen

共7兲

where the material parameters ␧0 and b are related to the softening strain and damage evolution rate, respectively; ␴0 = material strength; and D0 = ␴0b / ␧0 = initial elastic stiffness. Using the damage stiffness definition from relationship 共6兲, the incremental uniaxial damage stiffness Ds and the damage scalar ⍀ become Ds = 共1 − ⍀兲D0 = D0e−b共␧/␧0兲

To simulate the inelastic damage behaviors observed in asphalt materials, continuum damage mechanics was incorporated within the interparticle cementation model. Among the previous work in this area, the approach by Ishikawa et al. 共1986兲 was selected for use in our FEM. The theory was originally developed for concrete materials whereby the internal microcracks within the matrix cement and around the aggregates are modeled as a continuous defect field. Inelastic behavior is thus developed by the growth of damage within the material with increasing loading. A damage tensor 关⍀兴 is defined by considering the reduction of the effective area of load transfer within the continuum. The damage stiffness matrix 关Ds兴 can be expressed in terms of the initial elastic stiffness matrix 关D0兴 using continuum damage principles

共5兲

to describe the evolution of the defect field within the mastic cement. Such forms have been used before and for the uniaxial hardening response, the constitutive relation is taken as

Damage Mechanics Model

关Ds兴 = 共关I兴 − 关⍀兴兲关D0兴

兩␴z兩z=h¯共x − w1兲dx +

Similar calculations can be done for the case with unit normal, tangential, and rotational motions on particle 共2兲. These resultant force and moment calculations determine the various stiffness terms needed in the element Eq. 共1兲, giving the final result

− Knn 0

Knne Kttr1 Kttr21 +

w

0

=

F共1兲 t =

冕

共8兲

−b共␧/␧0兲

. In this damage model, the critical damage where ⍀ = 1 − e scalar ⍀c and critical strength ␴c are expressed as ⍀c = 1 − e−b and ␴c = ␴0共1 − e−b兲

共9兲

Note for this case, the critical value of the damage scalar ⍀c could be less than 1. Once the damage scalar reaches ⍀c, the material will gradually lose its stiffness. After critical strength, the softening behavior is taken as ␴ = ␴0共1 − e−b兲em共1−␧/␧0兲 ⇒

⳵␴ D 0m =− 共1 − e−b兲em共1−␧/␧0兲 ⳵␧ b 共10兲

where m = material parameter related to the softening rate. The corresponding incremental damage softening stiffness Ds and the damage scalar ⍀ become JOURNAL OF ENGINEERING MECHANICS © ASCE / JULY 2005 / 671

⌬Ut = cttw

where cnt,cnc,ctt = tension, compression, and shear softening factors. These criteria correspond to the average critical strength. Since the cementation geometry ¯h and w will in general be different for each particle pair, it is expected that each element will have different softening criteria related to its local microstructure. This damage behavior is a result of the mastic material’s defects that could include evolving microcracking leading to a separation or debonding between aggregate pairs. In order to simulate such total failure, elements were given a mastic failure criterion for tension, compression or shear based on the average failure strength ␴ f ␴ f = c f ␴c

Fig. 5. Uniaxial stress–strain response for damage model

Ds = 共1 − ⍀兲D0 = −

D 0m 共1 − e−b兲em共1−␧/␧0兲 b

where ⍀=1+

m 共1 − e−b兲em共1−␧/␧0兲 b

共11兲

The uniaxial stress–strain response corresponding to this particular constitutive model is shown in Fig. 5 for the case of ␧0 = 0.2, b = 5 and m = 1. This damage modeling scheme was incorporated into the finite-element network model by modifying the microframe element stiffness matrix given in Eq. 共5兲. Using the uniaxial relation 共8兲, the incremental normal and tangential damage stiffness terms for the hardening behavior can be written as 共Knn兲s = Knne−b共⌬un/⌬Un兲 共Ktt兲s = Ktte−b共⌬ut/⌬Ut兲

共12兲

and using Eq. 共11兲 the corresponding incremental normal and tangential damage softening stiffnesses are given by 共Knn兲s = − 共Knnm/b兲共1 − e−b兲em共1−⌬un/⌬Un兲 共Ktt兲s = − 共Kttm/b兲共1 − e−b兲em共1−⌬ut/⌬Ut兲

共13兲

where ⌬un and ⌬ut = normal and tangential accumulated relative displacements and ⌬Un and ⌬Ut = normal and tangential displacement softening criteria. Thus the microframe element incremental damage stiffness matrix 关Ks兴 is constructed from Eq. 共5兲 by replacing Knn and Ktt with 共Knn兲s and 共Ktt兲s. The initiation of mastic softening behavior for tension, compression and shear is governed by softening criteria based on accumulated relative displacements between particle pairs. A simple and convenient scheme to determine the softening criteria is based on using the dimensions of the interparticle mastic geometry in the form ¯ ⌬U共t兲 n = cnth ¯ ⌬U共c兲 n = cnch

共14兲

672 / JOURNAL OF ENGINEERING MECHANICS © ASCE / JULY 2005

共15兲

where c f = failure factor related to the average failure strength in each behavior, and ␴c indicates average critical strength in the corresponding behavior. The failure criterion for the uniaxial behavior is shown in Fig. 5 with the case of c f = 0.03. For this case, the failed elements still remain in the computation model, but their stiffnesses are very small and they carry very little load. For the compression behavior between particle pairs, the cementation spacing will decrease with load increment. Eventually, softening behavior will be initiated and the total element stiffness will significantly decrease. This will lead to the closing of the ¯ cementation gap when ⌬u共c兲 n = h, thus creating contact between the aggregates. At this point the element normal stiffness must be modified to account for this change of physics. The aggregate-toaggregate stiffness would be significantly higher than the cementation elastic stiffness and currently the model uses a contact stiffness three orders of magnitude larger than the elastic stiffness. This damage modeling scheme was incorporated into the ABAQUS finite-element code using the nonlinear user-defined element 共UEL兲 subroutine. The UEL subroutine determines the compression and tension force in each element and performs the required damage calculations in the normal and tangential directions based on the prescribed softening and failure criteria. Displacement controlled boundary conditions were employed and the modified Riks method was used in order to provide a more stable solution scheme. Also, because aggregate 共nodal兲 displacements became sizeable, the mesh geometry was updated during each load increment. A series of numerical indirect tension simulations were conducted on a model using different values of mastic softening factors and model parameters b and m. In our previous work, Sadd et al. 共2004b兲 included the indirect tension sample simulation response of vertical load versus displacement for a numerical model using three different compression softening factors. Since each case had identical elastic moduli and other model parameters, the initial responses were essentially the same. However, as the compression softening factor increased, less softening behavior was generated and these cases produced a higher maximum load. Similar trends were also found in simulations with different tension or shear softening factors, and the observed sensitivity could be ranked as compression⬎ shear⬎ tension. Parametric variation of model parameters b and m showed only small effects on the hardening damage evolution and softening rates of the overall model behavior.

Fig. 6. Image processing and model generation for indirect tension sample: 共a兲 grayscale image; 共b兲 binary image with segmentation technique; 共c兲 sieved aggregates 共D⬎1 mm兲 and least-squares elliptical fitting; and 共d兲 generated image model

Indirect Tension Image Model Generation and Simulation In order to capture real asphalt concrete microstructure, simulation material models were generated using imaging analysis procedures from photographic data of actual asphalt samples. The MATLAB image toolbox, DIPimage toolbox for MATLAB 共2004兲, and Adobe Photoshop were used for the image analysis. Surface electronic images were digitally processed and sieved to determine aggregate geometry using the ImagePro Code. A leastsquares curve-fitting routine 共EllipFit code兲 was developed within MATLAB and used to numerically fit the best ellipse to each sieved aggregate. Once the idealized elliptical particulate distribution was established, the finite-element modeling procedures could then be implemented to conduct a loading simulation on the scanned sample. The first model generation was conducted on a standard 102 mm diameter indirect tension sample shown in Fig. 6. A digital camera provided an electronic RGB image of the sectioned specimen, and this was converted to grayscale using Adobe Photoshop as shown in Fig. 6共a兲. The image processing steps can be summarized as follows: 1. Each grayscale pixel has a brightness value ranging from 0 共black兲 to 255 共white兲. A grayscale threshold was applied to convert the original images to a binary 共black and white兲, with white signifying the aggregates and black being the combination of mastic, air void and aggregate fines that can not be captured by the grid 关see Fig. 6共b兲兴. The image size 102⫻ 102 mm was digitized to 784⫻ 784 pixels, and this results in a 0.13 mm/ pixel grid. 2. Segmentation techniques were applied to the binary image to ensure that neighboring aggregate pixels did not artificially

connect together 关see Fig. 6共b兲兴. Aggregates were digitally sieved with 50-pixel area size, which lead to a sieving aggregate size less than 1 mm. 3. Each sieved aggregate was then labeled and selected as a separate image. The boundary pixels of each aggregate were extracted, and the coordinates of these pixels were stored in an array for the ellipse fitting. 4. A least-squares, ellipse-fitting algorithm developed by Pilu et al. 共1996兲 was incorporated to determine the “best” ellipse to represent each irregular aggregate geometry. The fitted ellipse was then generated for each sieved aggregate with the EllipFit code as shown in Fig. 6共c兲, and its geometry 共center coordinates, size, and orientation兲 were stored for use in the finite-element model generation and simulation. Based on the ellipse parameters obtained from the image analysis, a computation model was generated using MATLAB as shown in Fig. 6共d兲. Neighboring ellipses were maximally filled with cementation as illustrated in Fig. 3共b兲. It happened in our computation model that the domain of the mastic cement was repeatedly accounted for the stiffness calculation of the neighboring particles. And these would affect the computational accuracy of sample stiffness. The overlapped elements were recorded in our model generation code, and the stiffness calculation of these overlapped elements will be investigated in our future model study. For this computational model, the porosity was approximately 4% and the sieved aggregate percentage was 62.5%. For indirect tension simulation, a special contact boundary condition was imposed on aggregates in contact with the top and bottom loading plates shown as Fig. 6共d兲. The normal contact behavior was simulated by using very stiff elastic finite elements, and a small sliding displacement was allowed between the contact aggregates and the bearing plates to model frictional behavior. The x and y displacements of the bottom plate and x displacement of top plate were fixed. Sample loading was achieved by incrementing the y displacement of the top plate. The model material parameters were chosen as E = 47.5 MPa, ␯ = 0.3, b = 2, and m = 1.2, and softening factors of cnc = 0.2, and cnt = 0.1, and ctt = 0.1. In order to investigate the nature of the microstructural damage processes within the sample, a detailed numerical simulation was conducted on the image model shown in Fig. 6共d兲. The model assumes that there exists a continuous distribution of defects in the mastic material. This defect field is taken to grow with the material deformation, and this results in a softening response of particular elements in the FEM network. As per Eqs. 共12兲 and 共13兲, this softening behavior can affect the compression, tension, and shear behavior of the element. The image model was subjected to incremental loading, and during this process all elements within the model were monitored for softening behavior. The element softening evolution and load-displacement response for this model are illustrated in Fig. 7. The initial onset of sample loading is shown in Fig. 7共a兲 and this also indicates the initial tension and compression behavior within the finite element network. It was found that the largest compressional behavior occurred in elements near the loading plates, while elements across the loading centerline carried the greatest tension. Later loading steps are shown in Figs. 7共b–d兲, and the locations of these steps for each model are illustrated by the square points on the overall load–displacement curve 关Fig. 7共e兲兴. Darkened lines in Figs.7共b–d兲 indicate softening elements thereby illustrating considerable softening behavior with increasing load. It is observed that the evolution of softening damage primarily occurs in the central portions of the model sample where the element loading is the largest. The overall load–displacement response of the model JOURNAL OF ENGINEERING MECHANICS © ASCE / JULY 2005 / 673

Fig. 8. Fracture patterns in indirect tension test and simulation: 共a兲 damaged indirect tension specimen from test; and 共b兲 model simulation with removal of failed elements

Fig. 7. Image model softening patterns and load versus displacement curve with test data. Model parameters: E = 47.5 MPa, ␯ = 0.3, b = 2, m = 1.2, cnc = 0.2, ctt = 0.1, cnt = 0.1, and c f = 0.05.

is compared with experimental test data in Fig. 7共e兲, and simulation results compare favorably with the data. In order to simulate the fracture and failure behaviors observed in the actual test shown in Fig. 8共a兲, the failure criterion previously mentioned in Eq. 共15兲 has been incorporated into the damage model. All elements within this image model were monitored for the failure behavior based on the failure criterion with c f = 0.05. The simulated fracture pattern at the final loading step was obtained by removing the failed elements as shown in Fig. 8共b兲. Fig. 8共a兲 shows the actual damaged specimen from the laboratory test near the end of loading, and this compares well with the corresponding simulation result shown in Fig. 8共b兲.

Compression Image Model Generation and Simulation Using the techniques described in the previous section, a compression image model was also generated from an actual asphalt concrete sample with grayscale shown in Fig. 9共a兲. Aggregates were again digitally sieved from the sample binary image with a 50-pixel area size 共sieving size of 1 mm兲. The sieved aggregates 674 / JOURNAL OF ENGINEERING MECHANICS © ASCE / JULY 2005

Fig. 9. Compression image model generation: 共a兲 grayscale image; 共b兲 aggregate ellipse fitting; and 共c兲 image model

Fig. 11. Fracture patterns in compression test and simulation: 共a兲 damaged specimen of compression test; and 共b兲 model simulation with removal of failed element

appears that the microframe damage model has good ability to predict the sample global fracture behaviors. Fig. 10. Compression model softening pattern and load versus displacement curve with test data. Model parameters: E = 75 MPa, ␯ = 0.3, b = 1.0, m = 0.4, cnc = 0.04, ctt = 0.01, cnt = 0.02, and c f = 0.02.

were then least-squares fitted with elliptical particles as shown in Fig. 9共b兲. A numerical image model was generated based on these fitted particles with the maximally filled cement scheme as shown in Fig. 9共c兲. For this compression model, the porosity was approximately 1.5% and the aggregate percentage was 64.5%. The geometry dimensions were the same as the asphalt sample with length= 102 mm 共4 in.兲, height= 63 mm 共2.5 in.兲, and thickness = 22 mm 共7 / 8 in.兲. Displacement controlled boundary conditions were also used for the compression simulation. The x and y displacements of the particles on the bottom layer and the x displacements of the particles on the top layer were constrained. The y-displacement loading was incrementally imposed on particles of the top layer. The model material parameters were chosen as E = 75 MPa, ␯ = 0.3, b = 1.0, m = 0.4, and the softening factors cnc = 0.04, ctt = 0.01, cnt = 0.02, and the failure factor c f = 0.02. Compression model simulation and comparisons with test data are shown in Fig. 10. Fig. 10共a兲 shows the initial tension and compression force distribution within the microframe element network. Figs. 10共b–d兲 show the element softening evolution for three different loading steps indicated by the sample vertical displacement ⌬. Darkened lines again indicate softening elements and the line thickness is related to the softening level. The overall load–displacement response of this compression model favorably matches the test data shown in Fig. 10共e兲. As in the previous simulation, the model failure predictions were used for this compression sample. Fig. 11共a兲 shows the sample crack distribution observed in the actual compression test, while Fig. 11共b兲 illustrates the model fracture pattern in the last loading step by artificially removing the failed elements from the network. It again

Conclusions A micromechanical model has been used to simulate the twodimensional damage behavior of asphalt concrete. The aggregate– mastic microstructure was simulated with an equivalent finite element network that represented the load-carrying behavior between adjacent aggregates in the multiphase material. These network elements were specially developed from an elasticity solution for cemented particulates. Incorporating a damage mechanics approach with this solution allowed the development of a softening model capable of predicting typical global inelastic behavior found in asphalt materials. Simulations were conducted on image models with realistic aggregate microstructure. These image models were numerically generated from surface photographic data of actual asphalt samples. This was accomplished through image processing and numerical aggregate fitting. Two specific models were developed from indirect tension and compression samples. Comparisons of the overall load–deformation behavior and the internal damage/ failure response were made between the test data and simulation results. The overall load–deformation model results shown in Figs. 7共e兲 and 10共e兲 compared favorably with the experimental data, and it appeared that the damage model was able to correctly predict the inelastic softening response typically found in such tests. A further investigation was made into the evolution of the internal microdamage within the indirect tension and compression models. During incremental loading of the sample, all elements within the model were monitored for softening and failure behaviors, and the evolution of this damage was recorded. Similar photographic data from the actual tests were also collected on the surface damage/failure behaviors of each specimen. Comparisons of the model damage evolution with the experimental photographic data showed reasonable qualitative comparisons, and the JOURNAL OF ENGINEERING MECHANICS © ASCE / JULY 2005 / 675

model was able to generally predict the primary failure/fracture patterns observed. For example, within the indirect tension sample shown in Fig. 8, the model was able to predict the irregular failure/fracture pattern caused by the large aggregates on the top of the sample. The model also appeared to correctly predict the shear or mode II failure in the compression specimen as shown in Fig. 11. The current two-dimensional model is limited to the assumption of uniform behavior through the thickness of the sample. Clearly this assumption is not accurate and a three-dimensional extension of the model is underway.

Acknowledgments The writers would like to acknowledge support from the Transportation Center at the University of Rhode Island under Grant Nos. 01-64 and 02-86. Additional support was also provided from the Cardi Construction Corporation.

References Bahia, H., Zhai, H., Bonnetti, K., and Kose, S. 共1999兲. “Nonlinear viscoelastic and fatigue properties of asphalt binders.” J. Assoc. Asphalt Paving Tech., 68, 1–34. Bazant, Z. P., Tabbara, M. R., Kazemi, Y., and Pijaudier-Cabot, G. 共1990兲. “Random particle simulation of damage and fracture in particulate or fiber-reinforced composites.” Damage Mechanics in Engineering Materials, Trans. ASME, AMD 共Applied Mechanics Division兲 109, 41–55. Budhu, M., Ramakrishnan, S., and Frantziskonis, G. 共1997兲. “Modeling of granular materials: A numerical model using lattices.” Mechanics of deformation and flow of particulate materials, C. S. Chang, A. Misra, R. Y. Liang, and M. Babic, eds., Proc., McNu Conf., Trans. ASCE, Northwestern Univ., Evanston, Ill. Buttlar, W. G., and You, Z. 共2001兲. “Discrete element modeling of asphalt concrete: microfabric approach.” Transportation Research Record, 1757, Transportation Research Board, Washington, D.C., 111–118. Chaboche, J. L. 共1988兲. “Continuum damage mechanics, Parts I and II.” Int. J. Appl. Mech., 55, 59–72. Chang, C. S., and Gao, J. 共1997兲. “Rheological modeling of randomly packed granules with viso-elastic binders of Maxwell Type.” Comput. Geotech., 21, 41–63. Chang, G. K., and Meegoda, N. J. 共1993兲. “Simulation of the behavior of asphalt concrete using discrete element method.” Proc., 2nd Int. Conf. on Discrete Element Methods, MIT Press, Cambridge, Mass., 437– 448. Cheung, C. Y., Cocks, A. C. F., and Cebon, D. 共1999兲. “Isolated contact model of an idealized asphalt mix.” Int. J. Mech. Sci., 41, 767–792. DIPimage—A Scientist Image Toolbox for MATLAB, Version 1.5.0, June 29, 共2004兲. 具http://www.ph.tn.tudelft.nl/DIPlib典 Date accessed, 2004. Dvorkin, J., Nur, A., and Yin, H. 共1994兲. “Effective properties of cemented granular materials.” Mech. Mater., 18, 351–366. Guddati, M. N., Feng, Z., and Kim, Y. R. 共2002兲. “Towards a micromechanics-based procedure to characterize fatigue performance of asphalt concrete.” Proc., 81st Transportation Research Board Annual Meeting, Washington, D.C. Ishikawa, M., Yoshikawa, H., and Tanabe, T. 共1986兲. “The constitutive model in terms of damage tensor.” Proc. of Finite Element Analysis of Reinforced Concrete Structures, ASCE, Tokyo, 93–103. Kachanov, M. 共1987兲. “On modeling of anisotropic damage in elasticbrittle materials—A brief review.” Damage mechanics in composite, A. S. D. Wang and G. K. Haritos, eds., ASME, New York, 99–105. Krishnan, J. M., and Rao, C. L. 共2000兲. “Mechanics of air voids reduction 676 / JOURNAL OF ENGINEERING MECHANICS © ASCE / JULY 2005

of asphalt concrete using mixture theory.” Int. J. Eng. Sci., 38, 1331– 1354. Lee, H. J., and Kim, Y. R. 共1998兲. “Viscoelastic continuum damage model of asphalt concrete with healing.” J. Eng. Mech., 124共11兲, 1224–1232. Mora, P. 共1992兲. “A lattice solid model for rock rheology and tectonics.” The Seismic Simulation Project Tech. Rep. No. 4, 3-28, Institut de Physique du Globe, Paris. Mustoe, G. G. W., and Griffiths, D. V. 共1998兲. “An equivalent model using discrete element method 共DEM兲.” Proc., 12th ASCE Engineering Mechanics Conf., La Jolla, Calif. Papagiannakis, A. T., Abbas, A., and Masad., E. 共2002兲. “Micromechanical analysis of the viscoelastic properties of asphalt concretes.” Proc., 81st Transportation Research Board Annual Meeting, Washington, D.C. Pilu, M., Fitzgibbon, A., and Fisher, R. 共1996兲. “Ellipse-specific direct least-square fitting.” Proc., IEEE Int. Conf. on Image Processing, Lausanne, Switzerland. Rothenburg, L., Bogobowicz, A., and Haas, R. 共1992兲. “Micromechanical modeling of asphalt concrete in connection with pavement rutting problems.” Proc., 7th International Conf. on Asphalt Pavements, Nottingham Univ., U.K., Vol. 1, 230–245. Sadd, M. H., and Dai, Q. L. 共2001兲. “Effect of microstructure on the static and dynamic behavior of recycled asphalt material.” Year 1 Final Rep. No. 536108, Univ. of Rhode Island Transportation Center, Kingston, R.I. Sadd, M. H., Dai, Q. L., Parameswaran, V., and Shukla, A. 共2002兲. “Microstructural simulation of asphalt materials: modeling and experimental verification.” Proc., 15th ASCE Engineering Mechanics Conf., Columbia Univ., New York. Sadd, M. H., Dai, Q. L., Parameswaran, V., and Shukla, A. 共2004a兲. “Microstructural simulation of asphalt materials: Modeling and experimental studies.” J. Mater. Civ. Eng., 16共2兲, 107–115. Sadd, M. H., Dai, Q. L., Parameswaran, V., and Shukla, A. 共2004b兲. “Simulation of asphalt materials using a finite element micromechanical model and damage mechanics.” Transportation Research Record, 1832, Transportation Research Board, Washington, D.C., 86–94. Sadd, M. H., and Gao, J. Y. 共1997兲. “The effect of particle damage on wave propagation in granular materials.” Mechanics of deformation and flow of particulate materials, C. S. Chang, A. Misra, R. Y. Liang, and M. Babic, eds., Proc., McNu Conf., Trans. ASCE, Northwestern Univ., Evanston, Ill. Sadd, M. H., and Gao, J. Y. 共1998兲. “Contact micromechanics modeling of the acoustic behavior of cemented particulate marine sediments.” Proc., 12th ASCE Engineering Mechanics Conf., La Jolla, Calif. Sadd, M. H., Qiu, L., Boardman, W., and Shukla, A. 共1992兲. “Modeling wave propagation in granular materials using elastic networks.” Int. J. Rock Mech. Min. Sci 29, 161–170. Sangpetngam, B., Birgisson, B., and Roque, R. 共2003兲. “Development of an efficient hot mix asphalt fracture mechanics-based crack growth simulator.” Proc., 82nd Transportation Research Board Annual Meeting, Washington, D.C. Sepehr, K., Harvey, O. J., Yue, Z. Q., and El Husswin, H. M. 共1994兲. “Finite element modeling of asphalt concrete microstructure.” Proc., 3rd Int. Conf. Computer-Aided Assessment and Control Localized Damage, Udine, Italy. Simo, J. C., and Ju, J. M. 共1987兲. “Strain- and stress-based continuum damage models I. Formulation.” Int. J. Solids Struct., 23共7兲, 821– 840. Soares, J. B., Colares de Freitas, F. A., and Allen, D. H. 共2003兲. “Crack modeling of asphaltic mixtures considering heterogeneity of the material.” Proc., 82nd, Transportation Research Board Meeting, Washington, D.C. Trent, B. C., and Margolin, L. G. 共1994兲. “Modeling fracture in cemented granular materials.” Fracture Mechanics Applied to Geotechnical Engineering, Proc., ASCE National Convention, ASCE, Atlanta. Ullidtz, P. 共2001兲. “A study of failure in cohesive particulate media using

the discrete element method.” Proc., 80th Transportation Research Board Meeting, Washington, D.C. Wang, L. B., Frost, J. D., and Lai, J. S. 共1999兲. “Non-invasive measurement of permanent strain field resulting from rutting in asphalt concrete.” Transportation Research Record 1687, Transportation Research Board, Washington, D.C., 85–94.

Zhu, H., Chang, C. S., and Rish, J. W. 共1996兲. “Normal and tangential compliance for conforming binder contact I: Elastic binder.” Int. J. Solids Struct., 33, 4337–4349. Zhu, H., and Nodes, J. E. 共2000兲. “Contact based analysis of asphalt pavement with the effect of aggregate angularity.” Mech. Mater., 32, 193–202.

JOURNAL OF ENGINEERING MECHANICS © ASCE / JULY 2005 / 677