Construction and Building Materials 30 (2012) 783793

Contents lists available at SciVerse ScienceDirect

Construction and Building Materials

journal homepage: www.elsevier.com/locate/conbuildmat

Development of a micromechanical nite element model from computed

tomography images for shear modulus simulation of asphalt mixtures
Erdem Coleri a,, John T. Harvey a,1, Kai Yang b,2, John M. Boone b,3
a
University of California Pavement Research Center, Department of Civil and Environmental Engineering, University of California, One Shields Avenue, Davis, CA 95616, USA
b
Department of Radiology, University of California, Davis Medical Center, 4860 Y Street, Suite 3100 Ellison Building, Sacramento, CA 95817, USA

a r t i c l e i n f o a b s t r a c t

Article history: The main objective of this study was to develop 2D and 3D micromechanical nite element (FE) models
Received 12 August 2011 to predict the shear modulus of two asphalt mixture types, PG 64-28 PM (dense-graded polymer modi-
Received in revised form 7 December 2011 ed mixture) and RHMA-G (gap-graded rubberized mixture). The internal microstructure of the asphalt
Accepted 23 December 2011
mixtures were determined by X-ray computed tomography (CT) imaging. X-ray CT images were pro-
Available online 20 January 2012
cessed to create 2D and 3D meshed asphalt mixture FE model structures. Meshed model structures were
exported to the FE modeling software ABAQUS. Elastic and viscoelastic assumptions were used for mod-
Keywords:
eling aggregate and mastic subdomains, respectively. Shear modulus predictions for 2D and 3D models
Micromechanical
Microstructure
were compared to the measured values to determine the effectiveness of each model type for the simu-
X-ray computed tomography lation of the shear frequency sweep at constant height (FSCH) test. It was concluded that the 3D micro-
Superpave shear tester mechanical FE model is capable of predicting shear modulus at relatively high test temperatures with
Shear modulus high accuracy across a range of loading frequencies. On the other hand, 2D numerical models always
Finite element modeling under predicted the shear modulus values at all simulation temperatures and loading frequencies
Mastic because the reduced microstructure for the 2D models decreased the accuracy of numerical predictions.
Asphalt concrete It is expected that developed 3D models will serve as a valuable tool to understand certain problems in
Rubberized
the current design methods of asphalt concrete pavements.
Polymer modied
2011 Elsevier Ltd. All rights reserved.

1. Introduction etc.), asphaltaggregate interactions and the distribution of air voids

[4]. The mastic phase consists of ne aggregates, sand and other par-
Fundamental properties of composite materials have been ticles (rubber, polymer, etc.) embedded in asphalt binder. The mastic
investigated since the late 19th century [1]. Mechanical properties part appears to be the most challenging phase to model due to its
of composite materials were rst investigated by Albert Einstein in complex properties that depend on loading time and temperature.
1906 [2] in a study, in which he computed the effective viscosity of Dai and You [5] modeled the mastic properties by using the results
a uid containing rigid spherical particles. The technology of com- of creep tests that were conducted with mastic samples. Dai [6]
posite materials started to emerge in the early 1960s with the ad- developed a microstructure-based computational approach for
vent of modern ber composites [3]. stone based materials to predict both phase angle and complex mod-
Asphalt concrete is a composite material which consists of aggre- ulus under cyclic loadings with different frequencies. Papagiannakis
gates, air-voids and mastic. All of the constituents of asphalt con- et al. [7] also developed a 2D micromechanical FE model to predict
crete are large enough to be regarded as continua and the dynamic shear modulus and phase angle for asphalt mixtures tested
aggregate and mastic are bonded together at their interfaces. Perfor- with the Superpave shear tester. Rothenburg et al. [8] also intro-
mance and physical properties of asphalt mixtures are controlled by duced a discrete element model for asphalt mixtures to evaluate
the properties of the aggregate (shape, texture, gradation, modulus, pavement rutting performance. The model evaluated the effects of
etc.), properties of the asphalt binder (complex modulus, cohesion, both aggregatebinder interaction and inter-granular interactions.
Li et al. [9] developed a two-phase micromechanical FE model to pre-
dict asphalt mixture elastic modulus based on measured aggregate
Corresponding author. Tel.: +1 510 305 6246. and mastic properties.
E-mail addresses: [email protected] (E. Coleri), [email protected] Asphalt mixtures exhibit different performance characteristics at
(J.T. Harvey), [email protected] (K. Yang), [email protected] different temperature levels. Asphalt mixtures act as continuum
(J.M. Boone).
1 materials at low temperatures when the stiffness of the asphalt bin-
Tel.: +1 530 754 6409.
2
Tel.: +1 916 734 0314. der and the aggregate are similar [10]. However, at high
3
Tel.: +1 916 734 3158. temperatures, i.e. the critical conditions for rutting deformation

0950-0618/$ - see front matter 2011 Elsevier Ltd. All rights reserved.
doi:10.1016/j.conbuildmat.2011.12.071
784 E. Coleri et al. / Construction and Building Materials 30 (2012) 783793

accumulation, the aggregate can be two orders of magnitude stiffer of the aggregate was preserved during model development by
than the asphalt binder. As a result, as temperature increases or load using high pixel resolution X-ray CT images reconstructed at every
duration decreases the asphalt concrete pavement layers start to act 1 mm. The asphalt mixture was divided into three phases, air-
as particulate structures where the contact points between the voids, mastic and aggregates, by image processing and meshed
aggregates become more important in terms of rutting deformation for FE simulations. The aggregate and mastic phases were meshed
accumulation [11]. Therefore, at high temperatures permanent so that they share nodes on aggregates. The air-void volumes were
deformation accumulation mechanism in asphalt concrete layers left as empty voids in the model because they would not be able to
predicted by using continuum mechanics may not reect the actual resist any kind of stress applied throughout the simulation. Elastic
mechanism under in situ truck trafc. Thus, interactions between and viscoelastic assumptions were used for modeling aggregate
the particles and the effect of binder viscosity on the asphalt mixture and mastic subdomains, respectively. Shear modulus predictions
performance should be evaluated to identify the rutting perfor- for 2D and 3D models were compared to measured values to deter-
mance at high temperatures. mine the effectiveness of each model type for the simulation of the
Kandhal and Mallick [12] showed that aggregate shape, surface FSCH test.
texture, gradation and skeletal structure signicantly affect asphalt This paper focuses on the results for model development and val-
mixture rutting performance at high temperatures. Thus, micro- idation stages. However, with the proposed model, deciencies of
mechanical models that are developed by preserving the actual the Superpave shear tester can be identied and modications can
aggregate geometry and surface texture must be used to simulate be proposed based on the results of further simulations. In the cur-
laboratory tests in order to identify this proposed effect for differ- rent shear test set up, sample height is kept constant by applying an
ent mix types. Zhu and Nodes [13] investigated the effects of axial force at the top of the test sample in order to avoid any change
aggregate angularity on asphalt pavement response and concluded in sample volume. The effects of the missing shear couple [10], the
that increased angularity signicantly increases the stiffness of the missing moment in the opposite direction to avoid the change in
asphalt mixture. Dai and You [4] used 2D microstructures of as- sample height, on measured mixture performance can be investi-
phalt concrete that were obtained by optically scanning smoothly gated by using the FE model developed in this study. In addition,
sawn asphalt mixture specimens to predict their mixture complex the effects of specimen size and geometry on predicted mixture per-
and relaxation modulus at 20 C. Dai and You [4] preserved the formance will be determined to identify possible deciencies in the
actual 2D aggregate shape and angularity for particles larger than current AASHTO T320 [16] testing procedure.
1.18 mm to identify their effects on predicted material properties. The viscoelastic micromechanical nite element model devel-
You et al. [14] developed a 3D discrete element model by combin- oped in this study can also be used to identify the effects of differ-
ing a number of 2D discrete element models to determine the ent test (temperature, type of loading, load level, loading
improvement in predictions caused by the realistic interactions frequency, etc.) and material (air-void content, aggregate size,
in the third dimension. Although the 3D model resulted in better aggregate shape, angularity, etc.) parameters on shear resistance.
predictions for mixture modulus at temperatures lower than 0 C, With the developed FE model, sensitivity studies can be conducted
actual aggregate geometry and texture and its effects on predicted to identify the effects of these properties on measured shear per-
modulus were not simulated. You et al. [15] compared predicted formance. The effects of conning pressure on predicted shear
asphalt mixture stiffness at temperatures 18 C, 6 C and 4 C modulus can also be determined by using the micromechanical
with those predicted by 2D and 3D discrete element models gener- FE model developed in this study. Results of the simulations per-
ated from X-ray CT images. Results of the study indicated that 3D formed with and without conning pressure can be compared with
prediction results are very close to the measured data for all load- the heavy vehicle simulator (HVS) test results to identify the level
ing frequencies and temperatures. It was further concluded that 2D of improvement in predictions with the inclusion of conning
models under predicted the asphalt mixture stiffness in all cases. pressure to the laboratory tests. In addition, the effects of air-void
In the study presented in this paper, the micromechanical FE content and distribution on predicted shear performance can be
method was used for shear modulus simulation for asphalt mix- evaluated by using the proposed model.
tures at high temperatures. The micromechanical FE model was
developed to simulate the conditions for a standard frequency 3. Materials
sweep at constant height (FSCH) test [16]. Predicted shear modulus
Two asphalt mixture types were used in this study, a dense-graded mixture
values at various loading frequencies and temperatures were com- with PG 64-28 polymer modied (PM) binder and a gap-graded mixture with as-
pared to laboratory measured values to evaluate the effectiveness phalt rubber binder, PG 64-16 with crumb rubber modier (RHMA-G). Both mixes
of the developed model. Prediction accuracies of 2D and 3D micro- met Caltrans specications [17,18] for both mix and binder. Nominal maximum
mechanical models were determined separately to identify the aggregate sizes (NMASs) for the PG 64-28 PM and RHMA-G mixtures were
19 mm and 12.5 mm, respectively. Table 1 shows the aggregate gradations, asphalt
applicability of each model type for rutting performance predic-
binder contents and air-void contents of the asphalt mixtures and the mastic used
tion. Because the 3D model was developed from X-ray CT images in this study. NMAS for the mastic samples was selected to be 1.18 mm because
taken at a 1 mm interval distance, geometric characteristics of particles smaller than the 1 mm scanning interval may not be captured as aggre-
the aggregates in the nite element model appeared to closely gates during the X-ray CT imaging process. The air-void content for the mastic
match with the aggregates of the actual asphalt concrete sample. was assumed to be zero because of the high binder content and small aggregate
size.
The FE model developed in this study will be the rst 3D microme- Shear frequency sweep specimens were prepared for asphalt mixture and mas-
chanical model used for evaluating shear resistance of asphalt mix- tic materials for modeling. Field-mixed, laboratory-compacted (FMLC) specimens
tures at high temperatures. This model will be used for rutting were prepared adjacent to a construction site for HVS test sections by using jigs,
performance evaluation of different asphalt mixtures under vari- molds, and a rolling-wheel compactor. Loose mix was taken from the trucks with
a skip loader immediately prior to it being tipped into the paver and then dumped
ous temperatures and loading frequencies.
next to the preparation area. The required volume of material, based on the densi-
ties determined earlier in the laboratory, was weighed and then compacted into
molds at the same temperatures as those recorded on the HVS test sections. Asphalt
2. Objectives and scope mixtures used in this study were sawed from those ingots. Binders and aggregates
for mastic sample preparation were taken from the plant which produced the as-
phalt mixtures. Because NMAS for mastic samples was 1.18 mm and the binder
The primary objective of this study is to develop a microme- contents are very high in the mastic, mechanical compaction was not performed.
chanical FE model from X-ray CT images for shear modulus simu- However, samples were manually compacted in the laboratory with small size
lation of asphalt mixtures at high temperatures. The 3D geometry 7 kg cylinders.
 E. Coleri et al. / Construction and Building Materials 30 (2012) 783793 785

Table 1 STEP 1: Asphalt mixture preparation and imaging.

Aggregate gradations for asphalt mixtures and mastic.

Sieve size (mm) Percentage passing (%) 1. Mix preparation: Prepare FMLC asphalt mixture samples and
Mix Mastic Mix Mastic cut them to nal dimensions.
PG 64-28 PG 64-28 RHMA-G RHMA-G 2. Mastic preparation: Prepare mastic samples for viscoelastic
PM PM performance evaluation.
25.4 100 100 100 100 3. Acquire the 3D internal structure of the asphalt mixture
19 99 100 100 100 samples by X-ray CT imaging at the University of California
12.5 87 100 98 100 Davis Medical Center (UCDMC).
9.5 74 100 84 100
4.75 50 100 34 100
2.36 36 100 21 100 STEP 2: Image analysis and mesh development.
1.18 27 75 14 67
0.6 21 58 10 48 1. Determine mixture air-void content by standard CoreLok
0.3 16 44 7 33
method [20].
0.15 10 28 5 24
0.075 6.2 17.2 3.5 16.7
2. Process the X-ray CT images based on the measured air void
Asphalt binder 5.1 11.6 7.6 16.7 content to create the air void volumes.
contenta (%) 3. Determine aggregate threshold index (the intensity range
Air void content (%) 4.2 0.0 10.0 0.0 for aggregate images) by analyzing the processed images.
a
Design binder content (by dry weight of aggregate). 4. Subtract air void and aggregate domains from the complete
image to obtain mastic domain.
5. Exclude the air void domain from the mix (after exclusion,
volume occupied by air voids will be empty).
Mix preparation
and imaging

Mix preparation

CT imaging at UCDMC

Separate aggregate, mastic and air-

Image analysis and
mesh development

voids by image processing

Develop 2D and 3D image samples

Mesh development for FEM

Model mastic properties based on

FSCH test results
MM FEM development
and validation

Simulate FSCH test by FEM

Calculate G* Run mix FSCH tests and

measure G*

Compare measured and

calculated G*

Fig. 1. The general framework for the development of micromechanical FE model.

Asphalt mixtures and mastic samples for shear testing had rectangular prism
shapes. The width and length of the asphalt mixture samples were 130 mm while
the height was 48 mm. Standard cylindrical specimens with 150 mm diameter were
not used in this study to avoid complications related to the variable length-to-
height ratio along the circular surface area [10]. In addition, it would not be possible
to develop circular 2D numerical shear samples due to the direction of applied
shear load. The width and length for the mastic samples were selected to be
100 mm while the height was 15 mm. Aggregate modulus tests were not performed
since the effect of aggregate stiffness on shear deformation would be negligible due
to the effect of the low modulus values for the mastic at high temperatures. The
aggregate modulus was assumed to be 55 GPa [19]. It should be noted that aggre-
gate modulus would have a signicant effect on the asphalt mixture modulus when
the testing temperatures are low [14].

4. General procedure for model development

The general framework for the development of the microme- Fig. 2. Three dimensional image of asphalt mixture PG 64-28 PM generated from X-
chanical FE model is given in Fig. 1. The steps followed for model ray CT images. (a) Image with all three domains, air void, mastic and aggregate and
development are as follows: (b) air void distribution with in the asphalt mixture.
786 E. Coleri et al. / Construction and Building Materials 30 (2012) 783793

6. Create the high quality tetrahedral mesh for the aggregate 5. Run asphalt mixture laboratory tests at seven loading fre-
and mastic domains while preserving the empty air void quencies and three temperatures and measure G for each
domain. test.
7. Export the meshed structure to the FE modeling program. 6. Compare measured and predicted asphalt mixture G values
for both 2D and 3D models.
STEP 3: Micromechanical FE model development.

1. Run mastic FSCH tests at seven loading frequencies and 5. X-ray CT image acquisition, analysis and mesh development
three temperatures.
2. Develop viscoelastic model for the mastic phase using mas- 5.1. Image acquisition and processing
tic FSCH test results.
3. Develop FE model for the asphalt mixture for the corre- The internal structure of the asphalt mixture was measured by
sponding boundary conditions, material properties, loading X-ray CT imaging at the University of California Davis Medical Cen-
type and temperature. ter (UCDMC). The X-ray CT scanner reconstructs the 3D spatial dis-
4. Predict shear modulus (G) at seven loading frequencies and tribution of the attenuation coefcients of the materials that were
three temperatures that were used for mastic testing (21 scanned. Three dimensional images of the specimen are generated
simulations). by combining a series of 2D images generated by the scanner. The

Fig. 3. Comparisons between 3D virtual image structures and actual shapes for four of the aggregates. (a) Aggregate 1 from two perspectives. (b) Aggregate 2 from two
perspectives. (c) Aggregate 3 from two perspectives. (d) Aggregate 4 from two perspectives.
 E. Coleri et al. / Construction and Building Materials 30 (2012) 783793 787

quality of the 3D image depends on the X-ray CT image resolution.

In this study, horizontal planar images were acquired at 1 mm
intervals while the resolution cell for the other two dimensions
was determined to be 0.24 mm. At each sampling point within
Mastic Air void
the sample volume, each pixel from a CT image measures a value
that is related to the density and atomic number of the material
at that particular point. These measured values are then converted
to gray scale values [21]. The gray scale intensities for the aggre- (a)
gates, air voids and mastic range from 1000 to 3095, but these
values were rescaled to the range from 0 to 255. Because aggre-
gates are denser than mastic and air voids, they occupy the higher
portion of the intensity scale. On the other hand, low density air Aggregate
voids are at the lowest intensity portion while mastic intensity is (b)
between aggregate and air void intensities.
Simpleware software [21] was used in this study to create
masks for each domain based on measured intensity groups. Air (c)
void percentages for each asphalt mixture are shown in Table 1.
The total volume of the asphalt mixture sample was calculated Fig. 4. Two dimensional asphalt mixture (PG 64-28 PM) image before and after
for the complete intensity range, 0255. The upper limit for the mesh generation. (a) Before mesh generation. (b) Enlarged image with all three
air void intensity range was determined by trial and error to match domains. (c) After mesh generation.
the measured air void contents. The range for the PG 64-28 PM mix
was determined to be 0124 for a specimen with a calculated air ical complexity of the object. To ensure high quality triangulations,
void content of 4.3% while the range for RHMA-G mix was 099 i.e. to generate equilateral triangles, surface optimization was also
with a calculated air void content of 10.0%. These air void contents performed before the meshing procedure [22]. Fig. 4 illustrates a
are close to those measured with the CoreLok test (Table 1). The 2D asphalt mixture structure before and after mesh generation. It
volumes for aggregate and mastic were calculated to be 58.2% can be observed that triangle elements in the mastic and aggregate
and 37.5% for the PG 64-28 PM mix, respectively. For the RHMA- domains are perfectly bonded on the boundaries (Fig. 4). Meshed
G mix, aggregate and mastic volumes were calculated to be structures were exported to FE analysis software ABAQUS [23]
50.2% and 39.8%, respectively. for modeling.
After the thresholds for each domain were determined, masks
with different colors were assigned to clearly visualize the bound- 6. Micromechanical FE model development
aries between each domain. Because densities for the aggregates
and mastic were close to each other at some locations that are 6.1. Testing and viscoelastic model for asphalt mastic
close to the aggregate boundaries, aggregates appeared to be con-
nected at those boundaries. This segmentation problem may create FSCH tests were conducted with mastic samples at seven load-
a 3D structure in which all the aggregates are interconnected at the ing frequencies, 0.1 Hz, 0.2 Hz, 0.5 Hz, 1 Hz, 2 Hz, 5 Hz, and 10 Hz,
interface. This unrealistic material structure will result in over pre- and three temperatures, 30 C, 40 C and 50 C under a controlled
diction of stiffness as a result of those strong bonds. In order to shear strain level of 0.2%, with two replicates. It was observed that
avoid that problem, images were rst ltered with the software due to the small NMAS of mastic samples (1.18 mm), differences
and then manually processed to separate mastic from aggregate. between replicate test results appear to be negligible.
The 3D image for the PG 64-28 PM mix obtained after image pro- The temperature dependency of the mastic was dened by
cessing is given in Fig. 2 with the 3D air void distribution. In order using the WilliamsLandelFerry (WLF) equation [24], given as
to present the achieved 3D image accuracy in terms of aggregate follows:
shape, texture and angularity, four aggregates were selected from
C 1 T  T ref
the processed 3D image structure and compared to the aggregates logaT 1
extracted from the actual tested laboratory sample. Fig. 3 shows C 2 T  T ref
the comparisons between 3D virtual image structures and actual where aT is the timetemperature shift factor, C1 and C2 is the
shapes for those four aggregates. It can be observed that the shapes regression coefcients, Tref is the reference temperature, (40 C in
of the 3D aggregate structures closely match with the actual this study), and T is the test temperature.
shapes of the aggregates. For 2D model development, three slices
were randomly selected from the 3D structure for each mix type.
Predictions for the three 2D models were averaged and compared
to 3D model simulation results and laboratory measured shear
modulus values.
g1 g2 g3 gN-1 gN
5.2. Mesh development

The mask created for the air void domain was excluded from Ge 2 3 4 N-1 N
the material structure because it has no stiffness. Before meshing,
both the aggregate and the mastic domains were slightly smoothed
while preserving the element quality and volume. This pre-
smoothing stage allows the software to create more realistic
meshes. High quality (triangles that are as equilateral as possi-
ble) tetrahedral meshes were created for the aggregate and mastic
domains while preserving the empty air void domain. The size of
surface and volume elements were adjusted based on the topolog- Fig. 5. Schematic of the generalized Maxwell model.
788 E. Coleri et al. / Construction and Building Materials 30 (2012) 783793

Table 2 where a, b, c, d is the regression coefcients, and n is the reduced

Prony series parameters for Maxwell model. time.
i RHMA-G PG 64-28 PM Shift factors were calculated by tting the measured shear mod-
gi ki gi ki
ulus to the sigmoidal function (Eq. (2)). One shift factor was calcu-
lated at each test temperature while the shift factor for the
1 4.136E01 6.978E03 6.801E01 6.490E03
2 2.137E01 2.052E02 1.405E01 6.891E03
reference temperature (40 C) was set at zero. The Excel Solver rou-
3 2.041E01 1.129E01 1.284E01 7.703E02 tine was used to optimize the regression coefcients a, b, c, d and
4 1.043E01 1.228E+00 2.789E02 1.101E+00 shift factors for temperatures 30 C and 50 C. Regression coef-
5 2.932E02 2.500E+01 5.270E03 2.000E+00 cients C1 and C2 were calculated by simply tting the WLF equation
to the calculated shift factors.
The generalized Maxwell-type viscoelastic model was used in
this study to simulate the mastic time dependency. The schematic
Shear load direction of the generalized Maxwell model is given in Fig. 5. The model con-
sists of two basic units, a linear elastic spring and a linear viscous
dash-pot. Various combinations of these spring and dashpot units
dene the type of viscoelastic behavior. In this study, one spring

Constant height
and ve Maxwell elements in parallel were used for the mastic
Y model. Measured shear modulus (G) values were separated into
storage modulus [G0 (x)] and loss modulus [G00 (x)] components
based on the measured phase angles (d) using the following equa-
tions [7]:
X BC1: Change in height = 0;
G0 x jG xj cosd 3
BC2: Bottom platen is fixed in Y direction;
BC3: Peak movement at the top platen in +X = 0.096 mm
G00 x jG xj sind 4
Fig. 6. The loading and boundary conditions for the FSCH tests.
where x is the loading frequency in rad/s, G(x) is the complex
shear modulus in MPa, and d is the phase angle between strain
To optimize the regression coefcients C1 and C2, the shear and stress.
modulus data were rst tted to a sigmoid function, in the form Storage modulus and loss modulus can be expressed as follows
of [25]: [26]:

a X
N
xki 2
logGn d 2 G0 x Ge gi 5
1 expb c logn 1 xki 2
i1

Fig. 7. Vertical displacement (U2) distribution of 2D models at 50 C under a loading frequency of 0.2 Hz. (a) Asphalt mixture PG 64-28 PM and (b) asphalt mixture RHMA-G.
 E. Coleri et al. / Construction and Building Materials 30 (2012) 783793 789

X
N
xki 6.2. Prediction of shear modulus by FE modeling
G00 x gi 6
i1 1 xki 2
Micromechanical viscoelastic FE simulation was conducted on
meshed 2D and 3D samples that were exported to FE analysis soft-
where N is the number of Maxwell units, Ge is the equilibrium mod-
ware ABAQUS [23]. Plane strain assumption (exz = eyz = ezz) is used
ulus, gi is the relaxation strengths (spring constants of Maxwell
in the 2D FE analyses. FSCH tests were simulated for all computer
units), and ki is the relaxation times.
generated image samples at the same loading frequencies and tem-
gi peratures that were used for viscoelastic mastic model develop-
ki 7
gi ment. The purpose of the tests was to predict asphalt mixture
shear modulus at all temperature and loading frequency couples.
where gi = dashpot constants of Maxwell units (Fig. 5) Asphalt mixture tests were also conducted in the laboratory with
By tting Eqs. (5) and (6) to measured G0 and G00 data, the the physical parent of the virtual FE model sample to compare
parameters of a discrete relaxation spectrum can be determined. the numerical predictions to laboratory measured shear modulus
A genetic algorithm was used to optimize these model parameters values. In other words, asphalt concrete sample microstructure
(Ge, gi and ki) by minimizing the calculated residual sum of squares measured by X-ray CT imaging is used for FE model development
(RSSs) [27]. In this study, the following tness function was used to and predictions of the developed FE model were compared with
calculate the RSS for optimization: the results of the FSCH tests that were conducted with the original
02 32 2 32 1 asphalt concrete sample that was used for X-ray CT imaging. It
Xm
B4G 0
x j G00
x j C should be noted that for an asphalt mix type, one asphalt concrete
RSS @  15 4  15 A 8 sample was tested at all loading frequencies and temperatures in
j1
b0
G b 00
G
j j
order to avoid the effects of test results variability on comparisons.
FSCH tests are conducted by applying a cyclic shear load at the
where G b 0 and G
b 00 are the measured data at m frequencies. upper loading platen that is glued to asphalt sample while main-
j j
Optimized viscoelastic model parameters used for FE modeling taining the sample height constant. The bottom loading platen
is given in Table 2. To predict the global viscoelastic behavior of the was xed in the test setup that was used for this study. For the
asphalt mixture, timetemperature dependent FE analysis was FE simulation and laboratory tests, a constant sinusoidal shear
conducted by combining the elastic aggregate, viscoelastic mastic strain of 0.2% was applied. Fig. 6 illustrates the loading and bound-
and empty air void subdomains. ary conditions for the FSCH tests. Since the specimen height was

Fig. 8. The distribution of horizontal displacement for the mixture PG 64-28 PM. (a) 2D model structure and (b) 3D model structure.
790 E. Coleri et al. / Construction and Building Materials 30 (2012) 783793

48 mm for all 2D and 3D asphalt mixture samples, cyclic displace- of 0.2 Hz with the shear load applied at the top of the samples in po-
ment-based loading (strain controlled) was applied by xing the sitive direction 1 (left to right). Although the sample height was kept
upper platen movement at 0.096 mm for all cases. Because the constant to avoid the sample volume change during the tests, there
standard FSCH is an unconned test, conning pressures were were signicant vertical movements in the sample microstructure.
not applied around the sample walls [16]. The conning pressure It can be observed that asphalt mixture RHMA-G was exposed to
effect on shear modulus will be determined in a future study. higher levels of vertical displacement than the mixture PG 64-28
For 2D model simulations, three slices were randomly selected PM. In addition, vectors in the positive and negative two directions
from the 3D structure for each mix type. Predictions for the three were more uniformly distributed in mixture RHMA-G while locali-
2D models were averaged and compared to 3D model simulation zation of positive and negative direction vectors was observed for
results and laboratory measured shear modulus values. Because mixture PG 64-28 PM. However, it should be noted that having good
slices were randomly selected from the 3D structure, air void con- agreement between predicted and measured shear modulus values
tents for the 2D samples maybe different from the 3D sample air for the 3D model does not necessarily indicate similar microstruc-
void content. Although more than three 2D samples would be ture movement in all principal directions for the numerical and
needed to get a representative sample of predicted shear modulus, physical experiments. It is possible that horizontal movement can
only three 2D samples were used for each mix type due to time be effectively simulated based on linearized assumptions while
constraints. The average numbers of nodes for the 2D models were these assumptions may avoid the effective simulation of microstruc-
205,992 and 181,713 for the PG64-28 PM and RHMA-G mix mod- ture movement in the other two principal directions. Thus, displace-
els, respectively. The numbers of nodes for the 3D models were ment predictions in the vertical direction should also be validated by
3780,234 and 3573,522 for the PG64-28 PM and RHMA-G mix comparing the predicted vertical movement with the actual move-
models, respectively. The processing times for the 3D models ment in a shear test. The movement in the vertical direction during
(about 22 h) were about 1720 times higher than the processing a shear test can be measured with a high speed camera that can be
times for the 2D FE models. attached to the side of the sample. In this way, the effectiveness of
Fig. 7 illustrates the distribution of vertical displacement (U2) linearized assumptions for simulating vertical movements and
vectors for both mixture types at 50 C under a loading frequency mix dilative behavior in a shear test can be determined.

350
Predicted 50C Predicted 40C Predicted 30C
Measured 50C Measured 40C Measured 30C
300
Shear Modulus - G* (MPa)

250

200

150

100

50

0
0.1 1 10
Loading Frequency (Hz)
(a)
350
Predicted 50C Predicted 40C Predicted 30C
300 Measured 50C Measured 40C Measured 30C
Shear Modulus - G* (MPa)

250

200

150

100

50

0
0.1 1 10
Loading Frequency (Hz)
(b)
Fig. 9. Comparison of laboratory measured and 2D numerical model predicted shear modulus values. (a) Asphalt mixture PG 64-28 PM and (b) asphalt mixture RHMA-G.
 E. Coleri et al. / Construction and Building Materials 30 (2012) 783793 791

Fig. 8 shows the distribution of horizontal displacement for the modulus values. Comparison of laboratory measured and 3D
mixture PG 64-28 PM for both 2D and 3D simulations at 40 C un- numerical model predicted shear modulus values are given in
der a loading frequency of 0.2 Hz. Because the load was applied at Fig. 10. Results of the analysis indicated that the shear modulus
the top of the sample, highest level of horizontal displacement was values predicted by the 3D numerical model are close to the lab-
observed at that part. Similarly, no displacement was observed at oratory measured shear modulus values and the prediction of
the bottom of the sample since it was xed. shear modulus with the developed 3D model is reasonable. It
can be concluded that inability to capture the 3D microstructure
6.3. Comparison of predicted and measured shear modulus with the 2D models decreases the accuracy of numerical
predictions.
Numerical samples prepared for the two asphalt mixture types Fig. 11 illustrates the distributions for laboratory measured and
were tested at seven loading frequencies, 0.1 Hz, 0.2 Hz, 0.5 Hz, 2D3D numerical model predicted shear modulus values. Although
1 Hz, 2 Hz, 5 Hz, and 10 Hz, and three temperatures, 30 C, 40 C the distributions for 3D numerical model predictions are very close
and 50 C under a controlled shear strain level of 0.2%. For 2D to the laboratory measured values for both asphalt mixture types,
and 3D numerical shear test samples, asphalt mixture shear mod- prediction error increases at higher loading frequencies and lower
ulus values were predicted for each temperature and loading fre- temperatures when the mastic is stiffer. It can also be observed
quency couples. For both asphalt mixture types, calculations were from Fig. 10 that predictions with the highest error levels belong
performed by using the shear stressstrain curves of simulation to the lowest temperature (30 C) and highest loading frequencies
results that were calculated from shear load and horizontal dis- (2 Hz, 5 Hz and 10 Hz).
placement. Shear stress is basically the ratio of total force at In order to determine the effect of aggregate modulus on pre-
the loading surface (top surface) to the area of this surface while dicted mix shear modulus, 3D FE simulations were performed with
shear modulus is the ratio of calculated shear stress to the 0.2% an aggregate modulus of 75 GPa at the test temperature of 30 C
shear strain level used for FSCH testing. Fig. 9 shows the compar- for both PG64-28 PM and RHMA-G mixes. Results of the simula-
isons of measured and 2D numerical model predicted shear tions are given in Fig. 12. It can be observed that increasing the
modulus values for the two asphalt mixtures. It can be observed aggregate modulus from 55 GPa to 75 GPa does not cause any sig-
that the 2D numerical models always under predicted the shear nicant effect on predicted shear modulus.

350
Predicted 50C Predicted 40C Predicted 30C
300 Measured 50C Measured 40C Measured 30C
Shear Modulus - G* (MPa)

250

200

150

100

50

0
0.1 1 10
Loading Frequency (Hz)
(a)
350
Predicted 50C Predicted 40C Predicted 30C
Measured 50C Measured 40C Measured 30C
300
Shear Modulus - G* (MPa)

250

200

150

100

50

0
0.1 1 10
Loading Frequency (Hz)
(b)
Fig. 10. Comparison of laboratory measured and 3D numerical model predicted shear modulus values. (a) Asphalt mixture PG 64-28 PM and (b) asphalt mixture RHMA-G.
792 E. Coleri et al. / Construction and Building Materials 30 (2012) 783793

(a) 0.020 (b) 0.020

2D FEM
2D FEM
0.015 0.015

3D FEM Measured
Density

Density
0.010 Measured 0.010

3D FEM
0.005 0.005

0.000 0.000
0 80 160 240 320 400 0 80 160 240 320 400
Shear modulus (MPa) Shear modulus (MPa)

Fig. 11. Comparing distributions for laboratory measured and 2D3D numerical model predicted shear modulus values. (a) Asphalt mixture PG 64-28 PM and (b) asphalt
mixture RHMA-G.

350
Predicted 30C-Eagg=55 GPa
Predicted 30C-Eagg=75 GPa
300 Measured 30C
Shear Modulus - G* (MPa)

250

200

150

100

50

0
0.1 1 10
Loading Frequency (Hz)
(a)
350
Predicted 30C -Eagg=55 GPa
Predicted 30C -Eagg=75 GPa
300 Measured 30C
Shear Modulus - G* (MPa)

250

200

150

100

50

0
0.1 1 10
Loading Frequency (Hz)
(b)
Fig. 12. The effect of aggregate modulus on predicted shear modulus. (a) Asphalt mixture PG 64-28 PM and (b) asphalt mixture RHMA-G.

7. Summary and conclusions developed in this study is the rst 3D micromechanical model used
for evaluating shear resistance of asphalt mixtures at high temper-
In this study, 2D and 3D micromechanical FE models were atures. The microstructures of the FMLC mixture specimens were
developed to predict the shear modulus of two asphalt mixture determined by X-ray CT images taken at UCDMC. The asphalt mix-
types, PG 64-28 PM (dense-graded polymer modied mixture) ture was divided into three phases, air-voids, mastic and aggre-
and RHMA-G (gap-graded rubberized mixture). The FE model gates by image processing. The air-void phase was excluded from
 E. Coleri et al. / Construction and Building Materials 30 (2012) 783793 793

the model because it has zero resistance. X-ray CT images were views of the authors and do not reect the ofcial views or policies
processed to create 2D and 3D meshed asphalt mixture model of the State of California or the Federal Highway Administration.
structures. The aggregate and mastic phases were meshed to share
nodes on aggregates. Meshed model structures were exported to References
FE modeling software ABAQUS. Elastic and viscoelastic assump-
tions were used for modeling the aggregate and mastic subdo- [1] Rayleigh L. On the inuence of obstacles arranged in a rectangular order upon
the properties of the medium. Philos Mag 1892;34:481.
mains, respectively. Shear modulus predictions for 2D and 3D [2] Einstein A. Eine neue bestimmung der molekldimensionen. Ann Phys
models were compared to the measured values to determine the 1906;19:289306 [1911;34:5912. English translation in Investigations on
effectiveness of each model type for the simulation of FSCH test. the theory of Brownian motion. Dover 1956; p. 36-62].
[3] Hashin Z. Analysis of composite materials: a survey. J Appl Mech
It was concluded that 3D micromechanical FE model is capable 1983;50:481505.
of predicting shear modulus at high test temperatures across a [4] Dai QL, You Z. Prediction of creep stiffness of asphalt mixture with
range of loading frequencies. 2D numerical models always under micromechanical nite-element and discrete element models. J Eng Mech
2007;133(2):16373.
predicted the shear modulus values at all simulation temperatures [5] Dai Q, You Z. Micromechanical nite element framework for predicting
and loading frequencies because limited amount of microstructure viscoelastic properties of heterogeneous asphalt mixtures. Mater Struct
and the lack of aggregate interaction in the third dimension for the 2008;41(6):102537.
[6] Dai Q. Prediction of dynamic modulus and phase angle of stone-based
2D models decreased the accuracy of numerical predictions. It was
composites using micromechanical nite element approach. J Mater Civil
also identied that although the sample height was kept constant Eng 2010;22(6):61827.
and sample volume change was avoided during the tests, there [7] Papagiannakis T, Abbas A, Masad E. Micromechanical analysis of the
viscoelastic properties of asphalt concretes. J Transport Res Rec
were signicant aggregate movements in the sample microstruc-
2002;1789:11320.
ture. Asphalt mixture RHMA-G was exposed to higher levels of ver- [8] Rothenburg L, Bogobowicz A, Haas R. Micromechanical modeling of asphalt
tical displacement than the mixture PG 64-28 PM. concrete in connection with pavement rutting problems. In: Proc 7th int conf
on asphalt pavements, vol. 1; 1992. p. 23045.
[9] Li G, Li Y, Metcalf JB, Pang SS. Elastic modulus prediction of asphaltic concrete. J
Mater Civil Eng 1999;11(3):23641.
8. Future study [10] Weissman SL, Harvey JT, Sackman JL, Long F. Selection of laboratory test
specimen dimension for permanent deformation of asphalt concrete
Although the results obtained in this study provide a great deal pavements. J Transport Res Rec 1999;1681:11320.
[11] Weissman SL. The mechanics of permanent deformation in asphaltaggregate
of insight, further work is required to extend the analysis. The vis-
mixtures: a guide to laboratory test selection. Symplectic Engineering Corp.;
coelastic micromechanical nite element model developed in this 1997.
study should be used to identify the effects of different test and [12] Kandhal PS, Mallick RB. Effect of mix gradation on rutting potential of dense-
graded asphalt mixtures. J Transport Res Rec 2001;1767:14651.
material parameters on shear resistance. The effects of conning
[13] Zhu H, Nodes JE. Contact based analysis of asphalt pavement with the effect of
pressure on predicted shear modulus should be determined by aggregate angularity. Mech Mater 2000;32(3):193202.
using this 3D micromechanical FE model. Results of the simula- [14] You Z, Adhikari S, Dai Q. Three-dimensional discrete element models for
tions performed with and without conning pressure can be com- asphalt mixtures. J Eng Mech 2008;134 (12).
[15] You Z, Adhikari S, Kutay ME. Dynamic modulus simulation of the asphalt
pared with the HVS test results to identify the level of concrete using the x-ray computed tomography images. Mater Struct
improvement in predictions with the inclusion of conning pres- 2009;42(5):61730.
sure to the laboratory tests. In addition, the effects of air-void dis- [16] AASHTO. Standard method of test for determining the permanent shear strain
and stiffness of asphalt mixtures using the Superpave shear tester (SST):
tribution on predicted shear performance should be evaluated by designation T320. Washington, D.C.: AASHTO Provisional Standards; 2003.
using the proposed model. [17] California Department of Transportation (Caltrans). Highway design manual
The effect of specimen size on measured shear performance 2006 update chapter 630 exible pavement. Sacramento, CA: Caltrans; 2006.
[18] California Department of Transportation (Caltrans). Use of rubberized asphalt
should also be determined using the developed micromechanical concrete for new construction. Memorandum. Sacramento, CA: Caltrans; 2006.
FE model. Representative volume element requirements for differ- [19] You Z. Development of a micromechanical modeling approach to predict
ent test temperatures and mixture types can be proposed for FSCH asphalt mixture stiffness using discrete element method. Department of Civil
and Environmental Engineering. Ann Arbor, MI University of Illinois at
tests using this model. After performing factorial simulations with
Urbana-Champaign: UMI, a Bell & Howell Information Company; 2003. p. 512.
the developed FE model, a nomograph or a series of regression [20] AASHTO. T-331 bulk specic gravity and density of compacted asphalt
equations correlating shear modulus with the mastic shear modu- mixtures using automatic vacuum sealing method. Washington, D.C.:
AASHTO Provisional Standards; 2009.
lus, aggregate properties, temperature and loading frequency can
[21] Simpleware. ScanIP and +ScanFE software developed by Simpleware Ltd.,
be developed. In this way, mix shear modulus can be predicted Innovation Centre, Rennes Drive. Exeter EX4 4RN, UK; 2010 <http://
without performing any laboratory mix shear tests or computa- www.simpleware.com>.
tionally expensive FE simulations. Further extension of the micro- [22] Simpleware. Converting 3D images into models ScanIP Software + FE
module + CAD module tutorial guide. UK: Simpleware LTD. Innovation
mechanical model developed in this study will serve as a valuable Centre Exeter; 2010.
tool to improve the current design methods of asphalt concrete [23] ABAQUS. Users manual, version 6.5. Hibitt. Pawtucket, R.I.: Karlsson and
pavements. Sorenson; 2004.
[24] Williams ML, Landel RF, Ferry JD. The temperature dependency of relaxation
mechanisms in amorphous polymers and other glass forming liquid. J Am
Chem Soc 1955;77(14):37017.
Acknowledgments [25] AASHTO. AASHTO 2002 pavement design guide. Washington, D.C.: AASHTO
Provisional Standards; 2002.
This paper describes research activities that were requested and [26] Baumgaertel M, Winter HH. Determination of discrete relaxation and
retardation time spectra from dynamic mechanical data. Rheologica Acta
sponsored by the California Department of Transportation (Cal- 1989;28:5119.
trans), Division of Research and Innovation. Caltrans sponsorship [27] Tsai BW, Kannekanti VN, Harvey JT. Application of genetic algorithm in asphalt
is gratefully acknowledged. The contents of this paper reect the pavement design. J Transport Res Rec 2004;1891:11220.