Numerical analysis of asphalt mixture and comparison with physical Marshall test.

Rimsa, Vytautas ; Kacianauskas, Rimantas ; Sivilevicius, Henrikas 等

Introduction

Asphalt mix is a heterogeneous material composed of asphalt binder, aggregates and air voids. Approx. 85% of the total volume of the mix is the aggregates, asphalt binder constitutes around 10% of the total volume and the rest is air voids (Mitra et al. 2012). Required quantities of mineral fillers and binder are determined using asphalt mixture design methods such as Marshall test (Ozgan et al. 2013), Superpave (Guarin et al. 2013) or Hveem (Brown et al. 2009), which suggest various behaviours of an asphalt mix. Gradation of the inner structure of mineral materials is inconstant due to unavoidable segregation; the asphalt mixture from mineral materials can be designed using stochastic analysis methods (Vislavicius, Sivilevicius 2013; Sivilevicius, Vislavicius 2008). An asphalt mix designed under laboratory conditions and later reproduced in an asphalt mixing plant (AMP) must maintain the same parameters and deviations of the job mix formula (JMF) must not exceed maximum/minimum values (Braziunas et al. 2013; Braziunas, Sivilevicius 2010). Deviations of coarse aggregates, fine aggregates, mineral fillers and binder quantity from a JMF often lead to deterioration of mechanical characteristics of the asphalt mixture, such as decreased durability and life cycle.

Compared to structural analysis, application of numerical methods to asphalt mix modelling is also limited. Numerical application to asphalt mix modelling is basically limited to two discretization concepts, i.e. the Finite Element Method (FEM) and Discrete Element Method (DEM).

The FEM is one of the continuum discretization methods widely used for simulating structures and materials. Regarding asphalt mixtures, the straightforward application of the FEM to an asphalt mix model is rather limited due to the high heterogeneity of material structure. It is difficult to develop physically realistic constitutive model of the macroscopic continuum. Viscoelastic models extensively used in earlier applications are given in Schapery (2000), while recent developments to enhance continuum models of asphalt mixtures are discussed in Sun and Zhu (2013). Review of different simplified 2D and more realistic 3D modelling techniques using coupled viscoelastic, viscoplastic and viscodamage is given by Abu Al-Rub et al. (2012) and references herein.

Generally, macroscopic properties of heterogeneous solids are predefined by the interaction of stiffer particles with a weaker matrix. Series of FEM applications were applied to consider heterogeneity of asphalt mix type materials on meso scale (Dai et al. 2006; Hahn et al. 2010; Chareyre et al. 2011; Fakhari Tehrani et al. 2013; Huang et al. 2011; Ameri et al. 2011; Xia, Pan 2011; Gonzalez et al. 2007). Mixt discrete model presenting asphalt layer by the packing particles described as the system of the FEM was considered by Mo et al. (2008).

The Discrete Element Method (DEM) is another numerical technique that was recently suggested by Cundall and Strack (1979). This approach is elaborated for modelling of discontinuous materials, were individual particles are modelled as deformable interacting bodies. Calculation of contact interaction forces between particles is the most important issue of the DEM technique.

Typically, when DEM approach is used for analysis of granular materials, this method evaluates a unilateral contact; however, bilateral contact dominates in the interaction of particles found in the real asphalt mix type materials. Modelling of asphalt mix with the help of DEM techniques started in the early 1990s with Meegoda and Chang (1993). They created DE model based on granular approach. Here, aggregates were simulated as spheres, while bituminous binder was represent by two sets of viscoelastic elements in normal and tangential directions at each contact. Various interaction techniques were developed within the framework of the DEM methodology. Later, this modelling technique that uses bonded particles was improved and implemented into 2D and 3D particle flow codes PFC2D and PFC3D, respectively.

The micro fabric (MDEM) model in which various material phases (e.g. aggregates, mastic) are modelled with clusters of very small, discrete elements was suggested by Buttlar and You (2001). The advantage of this model is in analysis of particle forms and orientation as well as various contact models. Particle-based approach with empirically evaluated elastic normal interaction and the tangential spring-dashpot interaction was considered in models by Collop et al. (2006) and Liu et al. (2012). The contact model with adhesion between particles for simulating the interaction of bitumen binder with particles was applied by Magnanimo et al. (2012).

Continuum approach and lattice model were used by Kim and Buttlar (2009). Here, interaction stiffness has been derived base on continuum properties. Extensive review of various modelling approaches is given by Liu et al. (2011).

Most of development is relevant to a two-phase composite structure. In reality, contribution of smaller particles is also important. Multi scale analysis including filler-sized particles was presented by Underwood and Kim (2013). The model of the layered particle with elastic layer to simulate the bonding between asphalt mastic and aggregate was considered by Zhu et al. (2011)

It should be noted that asphalt mixtures are characterized by high strain gradients in the bond between particles and the matrix. The strain localization effects between two particles contacting via interface by applying of finite elements was considered in the previous paper issued by the authors (Rimsa et al. 2011). Recently, the semi-analytical model (Kacianauskas et al. 2013) of the normal interaction between two stiff spherical particles contacting via weaker interface has been derived. The model reflects the localisation of strain in the contact zone and is highly sensitive to the difference between the stiffness of the particles and the interface.

Our investigation addresses the interaction of layered particles. This article focuses on evaluating the contribution of the bitumen film layer of a particle. The model is applied to simulation of the Marshall test.

The paper is organized as follows. The modelling approach and basic data of the asphalt mixture as a particulate solid is formulated in Section 1, Discretization approach is described in Section 2. The analytical model and its evaluation are provided in Section 3. Numerical results and influence of various factors are discussed in Sections 4, 5 and 6. Final conclusions are drawn at the end.

1. Asphalt mixture as a particulate solid: approach and basic data

1.1. Basic definitions and the approach

In reality, asphalt-type mixtures are heterogeneous particulate solids composed from particles of various sizes and shapes, surrounded by the binder matrix. To facilitate development of the computational model, several assumptions regarding the particle shape and size are taken into account. Generally, the role of differently sized aggregate particles is different and different classifications basing on different assumptions exist. Comprehensive analysis of the role of particular ingredients may be found in the review by Liu et al. (2012). As already pointed out, aggregates with particle sizes less than 0.075 mm are termed as mineral fillers. The largest particles are termed as aggregates, where the cut-off size between coarse and fine aggregates is 2.36 mm. On the macroscopic scale, asphalt mixture is modelled as a particulate solid characterised by two fractions, particles and matrix, respectively. Thus, we follow the concept that fine aggregate, mineral filler and asphalt binder are mixed as an integrated component to interact with coarse aggregate particles. The integrated component known as asphalt sand mastic will be termed as matrix.

The volumes and the mechanical properties of each fraction are evaluated based on experimental measurements. The fraction of particles could be described by applying statistical data of aggregates. To simplify the problem, individual particles are assumed to be homogeneous and isotropic spheres. Voids are not considered.

The matrix is also assumed to be homogeneous. This assumption implies that the single or even of the larger amount of mineral filler particles has negligible influences on the behaviour of the solid. In real asphalt mix structure due to polarity of particles and heterogeneity of bitumen structure around the rigid stone, the transition layer between stiff particles and weak bitumen is created. This layer is termed as bitumen film. The homogenisation procedure is illustrated in Figure 1.

[FIGURE 1 OMITTED]

1.2. Mechanical properties

Mechanical properties of elastic particles are attributed to the features of granite from MatWeb (2011). The Young's modulus E is equal to 3-[10.sup.10] Pa while Poisson's ratio v = 0.25.

Viscoelastic properties of the matrix material correspond to those of the asphalt bitumen and were taken from Huurman and Woldekidan (2007). They are sensitive to environmental temperature and to loading rate. Time-dependent properties are characterised by the time dependent elasticity modulus E(t), which is constant throughout the volume. Here, time-independent elasticity constant [E.sub.0] = 46 MPa:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (1)

The values of constants [a.sup.E.sub.m] and [t.sub.m] are given in the

Table 1. The Poisson's ratio is equal to [v.sub.m] = 0.46 and indicates relatively high incompressibility.

To simplify the problem, the values correspond to conditions of the standard Marshall test and are defined at fixed temperature of 60[degrees]C. Binder mechanical parameters were tested at 20[degrees]C, while Marshall tests were made at 60[degrees]C. These differences were assumed based on William Landei and Ferry equation proposing to shift temperature factors for any temperature from the paper by Chailleux et al. (2006):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (2)

where: [T.sub.ref] - reference temperature, [T.sub.ref] = 20[degrees]C; T required temperature, [T.sub.i] = 60[degrees]C; [c.sub.1] - coefficient, [c.sub.1] = 15; [c.sub.2] - coefficient, [c.sub.2] = 160.

Mechanical parameters of the matrix film are very complicated. In this paper, the authors are focusing on two parameters--film thickness and modulus of elasticity - that impact on normal interaction of two particles. Different bitumen film values are given by Koroliov (1986), according to his book where:

[delta] = [BD.sup.m], (3)

where [delta] is bitumen film thickness in [micro]m; B = 11 and m = 0.84 are parameters of the material, for granite type material; D is diameter of the particle in mm.

Bitumen film modulus of elasticity reported by Kandhal and Chakraborty (1996) is from ten to thirty times greater compared with free bitumen. Analysis contribution of two matrix film parameters: stiffness [E.sub.f] and film thickness [delta]. Investigation focused on four film thicknesses [delta]: when 0 um (no film); 4 um; 12 um and 50 um, respectfully. In addition, the matrix film modulus of elasticity contribution for normal contact stiffness was investigated. This approach allows investigating a multi scale model and the influence of the indeterminate parameters and limits of their contribution for a macro scale model. In this paper, the authors investigated the effect of the matrix film on two particles of normal stiffness on the meso scale, and later, validated the results on the macro scale with physical Marshall test.

2. Discretization approach

2.1. Concept of the discrete model

Discretisation approach applied for modelling of particulate solids is schematically illustrated in Figure 2. It specifies conceptual issues related to the development of the discrete model of the highly heterogeneous solid. The representative specimen of the solid, in our case--the cylindrically shaped specimen is shown in Figure 2a.

[FIGURE 2 OMITTED]

The solid is regarded as a system of the finite number n of discrete spherical particles. Each of the particles i (i = 1, 2, n) is characterised by translational degrees of freedom. The particles are considered in the global Cartesian reference frame OXYZ. Thus, displacements of particle i are described by vector [[U.sub.i] = {[U.sub.xi], [U.sub.yi], [U.sub.zi],}.sup.T], while external forces acting in the centre of the particle by vector [f.sub.i] = [{[F.sub.xi], [F.sub.yi], [F.sub.zi],}.sup.T].

A DEM approach is employed for the development of the discrete model of a particulate solid. Therewith, the behaviour of the solid is described by the discrete model, which presents a structural network (Fig. 2d). An irregular triangular grid shapes the geometry of the network. The spherical particles are embedded into the nodes of the grid, while grid lines between the two adjacent points i and j are supposed to be deformable connection elements.

Description of the connection element has to be considered on the scale of the particle by considering representative volume (Fig. 2b). This volume presents a circular cylinder composed by two hemispheres connected by weaker interface material of the finite size usually termed as a bond.

It is assumed that normal interaction between centres of particles occurs due to motion of particle i with respect to j. To develop the simulation methodology, the connection element characterised in terms of the local variables has to be developed. The local coordinate frame pointing to the normal direction n is applied for description of the element. The local force [f.sub.n] presents resultant equilibrium of the nodal forces [f.sub.n] = [F.sub.i] x n = -[F.sub.j] x n, while the local displacements un present elongation of the element [u.sub.n] = ([U.sub.y] - [U.sub.t]) x n. The time-dependent constitutive relationship in terms may be derived for the normal force vector [f.sub.n] and displacements [u.sub.n]:

[f.sub.n](t) = [k.sub.n]([u.sub.n], t)[u.sub.n](t). (4)

The Eqn (4) presents a typical DEM model used for description of the normal behaviour of contacting particles, while [k.sub.n] ([u.sub.n], t) is the resultant normal stiffness parameter. In our case, Eqn (4) is more complex. It reflects various properties of contacting particles, matrix film and solid interface. Various material and geometric nonlinearities may be involved. The connection stiffness may be calculated numerically by conducting various deformation tests or applying analytical methods.

Once the connection element is characterised in the frame of the DEM methodology, the macroscopic analysis will follow the conventional path of the DEM.

On the other hand, such an element may be also regarded as synthetic 1D FE with user-defined properties. In the case of the FEM, the connection element may be characterised by the constitutive relationship in terms of the nodal vectors of the force [F.sub.ij] = [{[F.sub.i], [F.sub.j]}.sup.T] and the displacements [U.sup.ij] = [{[U.sub.i], [U.sub.j]}.sup.T]. Finally, it is defined in the conventional form:

[[k.sub.ij]] [U.sub.ij] = [F.sub.ij], (5)

where: [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is the resultant normal stiffness matrix containing interaction stiffness [k.sub.n].

3. Modelling strategy

The two-level strategy is applied for modelling of a particulate solid, while modelling sequence is schematically illustrated in Figure 2. The forward sequence of operations to discretise the particulate solid (Fig. 2a) is treated as the global analysis, which is aimed at evaluating the particle displacements. It comprises:

--Identification and generation of the structural network (Fig. 2d);

--Specification of the particle properties regarding the film layer (Fig. 2c), if necessary;

--Characterisation of the stiffness by Eqn (5) of connection elements (Fig. 2b);

--Evaluation of the time-histories of particle displacements by solving the DEM or FEM techniques. The backward sequence of operations aimed at evaluating the local variables such as strains or stresses by using the global solution may be also involved.

4. Theoretical background

4.1. Description of the connection element

Implementation of the above strategy implies the evaluation of the forces by Eqn (4) occurring during the binary interaction of two particles. It should be noted that evaluation of the element stiffness may be done in different ways and specific properties of the current problem will be taken into account.

The geometry of this volume is a cylinder composed of two interacting spheres with the radii [R.sub.i] and [R.sub.j], in our case [R.sub.i] = [R.sub.j], and the interface in-between is shown in Figure 3. The local Cartesian frame Oxyz positioned at the centre of the sphere i and attached to the connection line [i.sub.j] is introduced for the sake of convenience. Here, axis Oz points to the normal direction, while axes Ox and Oy belong to the tangential plane. The distance between the centres of interacting spheres is denoted by [L.sub.ij] while the separation distance between surfaces is [L.sub.g] (Fig. 3). The centres of spheres are defined by local coordinates [Z.sub.i] = 0 and [Z.sub.j] = [L.sub.ij], respectively.

The thickness of the film is denoted by o. This parameter is exceptionally used for evaluation of the contact stiffness. This value is considerably smaller when compared to the particle radius, [delta] [much greater than to] R; therefore, its influence to particle geometry is neglected.

It is obvious that a real cylinder presents inhomogeneous three-phase continuum body composed by five symmetrically arranged homogeneous sub regions as in (Fig. 3a). If we refer to the central line on Figure 3, five segments are denoted by points No: 1, 2, 6. Naturally, the cylinder is modelled by composition of five sequentially connected springs, stiffness of which [k.sub.12], [k.sub.56] are denoted by numbers of connecting points as shown in Figure 3.

[FIGURE 3 OMITTED]

Sequential model assumes that the total system stiffness is the sum of individual components. This condition written in terms of element stiffness gives required expression of the total stiffness k, obtained from the classical model:

[([k.sub.n]).sup.-1] = [([k.sub.12]).sup.-1] + [([k.sub.23]).sup.-1] + [([k.sub.34]).sup.-1] + [([k.sub.45]).sup.-1] + [([k.sub.56]).sup.-1]. (6)

Evaluation of particular components requires theoretical justification.

4.2. Constitutive model

For detailed description of the governing relationship, deformation behaviour of the representative volume should be evaluated, and proper model should be derived. Regarding the 3D continuum, the internal state variables, i.e. displacements u(x, y, z), strains [epsilon](x, y, z) and stresses [sigma](x, y, z), are functions of an arbitrary position point. They are attributed to the reference configuration at the time instant [t.sub.0].

Now, we can refer to the beam theory by assuming that the cylinder follows the assumptions used for homogeneous axially loaded bar. These assumptions have several consequences. The cross-section remains plane and rigid. Assumption on the rigid section implies zero cross-sectional displacements:

[u.sub.x] (x, y, z) = [u.sub.y] (x, y, z) = 0. (7)

Further we assumed that deformation shape of the imaginary cylinder bar may be described by comprising deformation of separate lines parallel to the axis of the cylinder, or so-called generatrices.

In this case, we are working with the non-plane sections. Defining location of the line by coordinates x and y, its length is L(x, y). The state of the line could be characterised by the longitudinal variables, i.e. displacement component [u.sub.z](x, y, z), strain [[epsilon].sub.zz] (x, y, z) and stress [sigma].sub.zz](x, y, z). The subscript z will be omitted for the sake of simplicity.

First, kinematics of the deformed cylinder is considered. Since motion of the cylinder is described by the displacements of the rigid end-sections, u(x, y, 0) = 0 and u(x, y, L(x,y)) = [u.sub.n], these values along with Eqn (7) may be interpreted as boundary conditions for displacements. Assuming a simplified linear approximation law for the longitudinal displacement component, it is defined by simple equation:

u(x, y, z) = Z/L(x, y) [u.sub.n]. (8)

Regarding the above Eqn (8) and restricting by the small strain and small displacement relations, the longitudinal strain [[epsilon].sub.zz] = [partial derivative][u.sub.z]/[partial derivative]z presents relative elongation of the line:

[epsilon](x, y, z) = 1/L(x, y)[u.sub.n]. (9)

The constitutive relationship is attributed to the reference configuration, and the instantaneous linearity is assumed. The axial stress in the line [sigma] is proportional to the strain e. To ensure incompressibility of the section, or the transversal incompressibility of the cylinder, the confined state issued and the confined elasticity modulus [??] = E/[1 - [v.sup.2]] will be applied in the final equation:

[sigma](x, y, z) = [??] 1/L(x, y) [u.sub.n]. (10)

The stiffness of the cylindrical element will be found applying the principal of minimal deformation energy approach. Firstly, we define the density of the deformation energy accumulated in arbitrary point of the volume by a scalar product:

w(x, y, z) = 1/2 [sigma](x, y, z)[epsilon](x, y, z). (11)

Average energy of separate line:

[w.sub.s] (x, y) = [[integral].sup.L(x, y).sub.0] w(x, y, z)dz, (12)

is defined as energy of elementary prism of length L(x, y) with elementary section area dA. By inserting stress (Eqn (9)) and strain (Eqn (8)) into Eqn (11) and by integration along the length yields to expression:

[w.sub.s] (x, y) = 1/2 [??] (1/L(x, y)) [u.sup.2.sub.n]. (13)

The stiffness of the cylindrical bar is defined as the second derivative of the total deformation energy [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] with respect to displacement:

[k.sub.n] = [[partial derivative].sup.2] W/[partial derivative][u.sup.2.sub.n]. (14)

Finally, it is obtained by integration of expression (14) over the section area and differentiation:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (15)

Denoting [r.sup.2] = [x.sup.2] - [y.sup.2] and regarding axial symmetry, the integral stiffness could be defined in cylindrical coordinates:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (16)

Formally, the above Eqn (15) or (16) could be applied for evaluation of stiffness of the cylinder with the non-plane section, as well as for sub volumes of the cylinder.

4.3. Stiffness of particular phases

Local stiffness of each of the three phases (Fig. 3) will be evaluated independently on its neighbours. Precise evaluation would be available, if deformed shape on interface surface was known. In other cases, expression (16) is reasonable for practically undeformed surface.

Let us consider the contribution of individual phases. It should be reminded that geometry of the particle is actually defined as that for the layered particle including the film layer. The total radius of the particle including the film layer may be expressed as R = [R.sub.p] + [delta]; where [R.sub.p] is the radius of the particle without the film. Elasticity properties of the particle are defined by the elasticity modulus E and the Poisson's ratio v. The properties of the film and interface bond are identified by subscripts f and b while elasticity constants are denoted by [E.sub.f], [E.sub.b] and [v.sub.f], [v.sub.b], respectively.

The volume of the particle in the cylinder is shaped by the hemi-sphere, which is limited by two surfaces. The background circular plane is fixed at z = 0 and presents the mid-plane of sphere of the radius R. Simultaneously, it is the normal section of the cylinder.

The spherical surface defined by the radius R is described by the equation z (x, y) = [square root of [R.sup.2] - [r.sup.2]], which automatically predefines the length of the arbitrary coaxial line L(x, y) = z(x, y). The maximal length, which is also computational length of the spring, is at the position of the central line, i.e. [L.sub.max] = L(0, 0) = R.

Now, the stiffness of the particle [k.sub.p] = [k.sub.12] with the radius Rp may be evaluated by the general Eqn (16). Substitution of particle data yields:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (17)

Then, the integral defined by Eqn (17) can be obtained explicitly as follows:

[k.sub.12] = 2[pi][??][R.sub.p]. (18)

Geometry of the film layer (Fig. 4) is simplified, and only normal projection is considered for modelling purposes. Consequently, the stiffness of the film [k.sub.f] = [k.sub.23] is defined as stiffness of the cylinder bar having radius [R.sub.p] and the height [delta]. Formally, it is obtained according to the expression (16). Substituting of L(r) = [delta] yields the final expression of the film stiffness:

[k.sub.23] = [??] [pi][R.sup.2.sub.p]/[delta]. (19)

Assuming the sequential location of the particle and the film, the resultant effective stiffness of the layered particle [k.sub.13] = [k.sub.12] x [k.sup.23]/([k.sub.12] + [k.sub.23]) can be expressed explicitly:

[k.sub.13] = 2[pi][[??].sub.p][R.sub.p]/1 + [[??].sub.f] [delta]/[[??].sub.p][R.sub.p]. (20)

This allows the analytical evaluation of the bitumen film contribution.

[FIGURE 4 OMITTED]

Now, the stiffness of the interface bond [k.sub.b] = [k.sub.34] may be evaluated in the same manner by the general Eqn (16). Computational lengths of co-axial lines of the interface are expressed as follows:

[L.sub.int] (r) = [L.sub.g] + (R - z (r)). (21)

Here, the spherical surface is defined by the equation z (r) = [square root of [R.sup.2] - [r.sup.2]].

The stiffness [k.sub.13] of the particle together with the matrix film is obtained assuming that the springs representing the particle and the matrix film are connected sequentially:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (22)

Integral in Eqn (22) in our case is obtained numerically.

The total normal stiffness defined according to Eqn (6) is now simplified and may be expressed as follows:

[k.sub.total] = [k.sub.13] x [k.sub.34] x [k.sub.46]/[k.sub.13] x [k.sub.34] + [k.sub.34] x [k.sub.46] + [k.sub.13] x [k.sub.46], (23)

where particular stiffness of the particle [k.sub.13] is defined by Eqn (20) while the interface stiffness [k.sub.34] is defined by Eqn (22).

4.4. Analysis of the normal stiffness of two particles for various matrix mechanical parameters

Contribution of the film is investigated numerically by considering variations of the film thickness and the elasticity modulus variation are given in Figures 5 and 6.

As Figure 5 suggests, investigation showed that the increase in the matrix thickness from 12 um to 50 um or by 24% gives up to 56% increase in the total normal contact stiffness. Thinner matrix film does not significantly change the total contact stiffness ratio.

[FIGURE 5 OMITTED]

[FIGURE 6 OMITTED]

It could be observed that the increase in the film thickness or modulus of elasticity of the film reduces the strain gradient in the contact zone between a particle and the matrix. Influence of fill effects is given in Table 2.

Variation of relative matrix stiffness is shown in Figure 6. The main conclusion regarding the film thickness and modulus of elasticity investigation suggests that the increase in bitumen film thickness increases the contact stiffness, while the increase of elasticity modulus is important up to the ratio 100[E.sub.23] = [E.sub.45] and later additional increase in film elasticity modulus has no influence on the contact stiffness of the system.

5. Marshall test

5.1. Set-up

One of the widely used tests for characterization of asphalt mix materials is Marshall test. Marshall asphalt mix design method was originally developed by Bruce Marshall of the Mississippi Highway Department around 1939. The advantage of this test consists in complex characterization of the mechanical stiffness of the entire heterogenic structure. Basic scheme is shown in Figure 5. Figure 6 provides the typical test result force-displacement graphic. Due to high stiffness nonlinearity and temperature sensitivity of the matrix, this test performs in strong controlled load speed and temperature.

According to the Marshall test design method for determination of the optimal mix gradation, Marshall test specimens have to be created and testes for determination of physical and mechanical parameters. Creation and testing of Marshall test specimens requires the use of expensive specialized laboratory instruments and highly qualified staff.

A promising alternative approach uses numerical modelling for asphalt mix stiffness analysis taking into account inner structure granulometry, binder quantity and mechanical properties, and can replace laboratory tests or decrease the number of required tests. Concentrated information about Marshall test is given in European Standards LST EN 12697-34 (2004), White (1985), Roberts et al. (1996). The objective of Marshall test is a simple apparatus for design and control of asphalt paving mixtures. This method is widely used because:

1. It was designed to analyse the entire specimen.

2. It is fast testing with minimum cost, compact and portable.

3. Density of testing specimen could reasonably respond to road densities.

Marshall test is performed for a specially designed asphalt cylinder specimen with dimensions of 102 mm in diameter and 64 mm in height (Fig. 7). The Marshal test curve can be divided in to three stages (Fig. 8): I--adapting the specimen to the form of the metal plates; II--linear loading stage; and III--nonlinear stage with the maximum force result and high plasticity deformations. Most important parameters are maximum force reaction or "Marshall stability" and displacement value or "flow" at the maximum force point.

[FIGURE 7 OMITTED]

[FIGURE 8 OMITTED]

The field of interests of this paper is II linear stage of the Marshall curve. New elaborated models (discussed later) will be validated by tangential angles.

5.2. DE model for Marshall test analysis

To simulate the Marshall test specimen, discrete model shown in Figure 2e was developed.

The discrete model was created in the way that interaction of two contacting particles is represented by 1D element normal stiffness.

The specimen was created by C++ program code by spreading out the particles randomly in space according to Delaunay's triangulation methodology. The programme code was written in a way that virtual particles could not overlap; creation of 1D elements stops once the required density is reached (in this case 77%). This means that no volume correction coefficients are required for full compliance of physical and digital inner volume structures.

FEM specimen represents 1012 nodes and 4812 elements.

Creation of DE model for digital Marshall test analysis is one of the critical steps. Detail representation of inner specimen structure increases model accuracy and drastically increases the CPU time.

[FIGURE 9 OMITTED]

In terms of the granulometric composition of AC 16 PD 70/100 asphalt concrete, the percentage distribution by mass of particles in diameter D < 1 mm amounts to ~34% (Fig. 9). This paper focuses on asphalt concrete type AC 16 PD 70/100, which is widely used for road structures. One of the most important parameters that characterize asphalt concrete is granular particle percentage distribution by size given in Figure 9. Particle size varies form 0.075 mm to 16 mm. In this paper, the inner model structure was created using the average particle size approach, which is discussed in Mo et al. (2008) paper, where particles less than 1 mm in size are not taken in to account; for the average particle size calculation, the average particle size is:

[bar.d] = [n.summation over of (i=1)] [w.sub.i][[bar.d].sub.i]/[n.summation over of (i=1)][w.sub.i]. (24)

where: n is the total number of particles in the specimen; d is the average diameter of a particle, mm; [d.sub.i]--i-particle size, mm; [w.sub.i]--percentage distribution by mass; [[bar.d].sub.i] = ([d.sub.i] + [d.sub.i-1])/2.

[FIGURE 10 OMITTED]

According to Eqn (27), the average particle size D = 7.1 mm in AC 16 PD 70/100 asphalt concrete. Theoretical matrix thickness:

h = r[(.sup.k-1/3] - 1), (25)

where: r - average particle radius; r = 3.55 mm; k - volume fraction of particles in a mixture, k = 0.77; for asphalt concrete AC16 PD 70/100, the average matrix thickness is h = 0.32 mm.

The effect of particles which are smaller than 1 mm, effect evaluated by increase matrix stiffness[Jl] (Fig. 10). According to given asphalt AC16 PD 70/100 data, the quantity of particles smaller than 1 mm is 34% by mass or 24% by volume, while the stiffness scale factor is k = 2.

6. Results and discussion

Deformation behaviour of the Marshall specimen during is illustrated in Figure 11. Basic results of Marshall test show that vertical displacements dominate during the specimen loading and maximum values are exceeded at the zones that interact with steel plates. Non uniform displacements distribution can be explained by nonuniform inner structure of the DE model.

[FIGURE 11 OMITTED]

Influence of the matrix film on force reaction is given in Figure 12. Maximum differences between physical digital tests are at the initial loading phase at 0.8 seconds exceeding ~9%. In this physical test, the most important parameter is the maximum force reaction at deformations 3.8 mm [F.sub.AC16PD70/100] = 9.5kN, while [F.sub.DE] = 10.9 kN; at this point, the difference between physical and DE models is ~13%; this indicates an acceptable error.

Influence of the matrix film increases with the increase of the film thickness. The force reaction ratio between no film and model with 4 um is only ~0.29%. The maximum influence on force reaction is achieved when the matrix film thickness exceeds 50 um; then, the force reaction increases by 6.5% compared with the initial model, and the difference between normal contact stiffness of two particles without matrix film and with bitumen film is ~5.7%. This shows that normal stiffness of two particles can be directly recalculated to the stiffness of the entire Marshall specimen, for given particles size, matrix stiffness.

[FIGURE 12 OMITTED]

[FIGURE 13 OMITTED]

For better understanding of the behaviour of micro parameters, additional investigations should be made on the matrix film thickness and the contribution of the matrix stiffness on the macro model (Fig. 13, Tables 3 and 4).

Numerical model of Marshall test shows that the increase of matrix stiffness up to 100[E.sub.f] = [E.sub.b] has a positive influence on the macro model, and additional increase of the film modulus of elasticity does not have any influence on the system of the entire system.

Conclusions

The semi-analytical model for evaluation of the binary normal contact of spherical particles interacting via weaker multi-layered interface was proposed and its applicability for modelling of the asphalt mixtures was investigated. It was shown that the contribution of the film layer could be evaluated by the proposed model. On the other hand, the binary model serves as the base for global discretisation of the structure of particulate composite. Adequacy of the discrete model to the asphalt mixture was verified with experimentally obtained results for typical asphalt Marshall test. Calculations show that the developed discrete approach gives very close secant stiffness compared to real physical test, as the difference exceeds only 0.9% at viscoelastic stage. Full comparison to the model requires evaluation of plastic deformations of the matrix.

References

Abu Al-Rub, R. K.; Darabi, M. K.; Huang, C. W; Masad, E. A.; Little, A.; Dallas, N. 2012. Comparing finite element and constitutive modelling techniques for predicting rutting of asphalt pavements, International Journal of Pavement Engineering 13(4): 322-338. http://dx.doi/org/10.1080/10298436.2011.566613

Ameri, M.; Mansourian, A.; Heidary Khavas, M.; Aliha, M. R. M.; Ayatollahi, M. R. 2011. Cracked asphalt pavement under traffic loading--a 3D finite element analysis, Engineering Fracture Mechanics 78(8): 1817-1826. http://dx.doi/org/10.1016/j.engfracmech.2010.12.013

Braziunas, J.; Sivilevicius, H. 2010. The bitumen batching system's modernization and its effective analysis on the asphalt mixing plant, Transport 25(3): 325-335. http://dx.doi/org/10.1061/(ASCE)MT1943-5533.0000646

Braziunas, J.; Sivilevi?ius, H.; Virbickas, R. 2013. Dependences of SMA mixture and its bituminous binder properties on bitumen batching system, mixing time and temperature on asphalt mixing plant, Journal of Civil Engineering and Management 19(6): 862-872. http://dx.doi/org/10.3846/13923730.2013.843587

Brown, E. R.; Kandhal, P.; Roberts, F.; Kim, Y. R.; Lee, D. Y. 2009. Hot mix asphalt materials, mixture design and construction. Landham, Maryland: NAPA Research and Education Foundation. 720 p.

Buttlar, W.; You, Z. 2001. Discrete element modeling of asphalt concrete: microfabric approach, Journal of the Transportation Research Board 1757: 111-118. http://dx.doi/org/10.3141/1757-13

Chailleux, E.; Ramond, G.; Roche, C.; Chantal, R. 2006. A mathematical-based master-curve construction method applied to complex modulus of bituminous materials, Road Materials and Pavement Design, EATA 2006, 75-92. http://dx.doi/org/10.1080/14680629.2006.9690059

Chareyre, B.; Cortis, A.; Catalano, E.; Barth, E. 2011. Pore-scale modeling of viscous flow and induced forces in den se sphere packings, Transport in Porous Media 92(2): 473-493. http://dx.doi/org/10.1007/s11242-011-9915-6

Chen, J. S.; Kuo, P. H.; Lin, P. S.; Huang, C. C.; Lin, K. Y. 2007. Experimental and theoretical characterization of the engineering behavior of bitumen mixed with mineral filler, Materials and Structures 41 (6): 1015-1024. http://dx.doi/org/10.1617/s11527-007-9302-5

Collop, A. C.; McDowell, G. R.; Lee, Y. W. 2006. Modelling dilation in an idealised asphalt mixture using discrete element modelling, Granular Matter 8(3-4): 175-184. http://dx.doi/org/10.1007/s10035-006-0013-3

Cundall, P. A.; Strack, O. 1979. A discrete numerical model for granular assemblies, Geotechnique 29: 47-65.

Dai, Q.; Sadd, M. H.; You, Z. 2006. A micromechanical finite element model for linear and damage-coupled viscoelastic behaviour of asphalt mixture, International Journal for Numerical and Analytical Methods in Geomechanics 30(11): 1135-1158. http://dx.doi/org/10.1002/nag.520

Fakhari Tehrani, F.; Absi, J.; Allou, F.; Petit, C. 2013. Heterogeneous numerical modeling of asphalt concrete through use of a biphasic approach: porous matrix/inclusions, Computational Materials Science 69: 186-196. http://dx.doi/org/10.1016/j.commatsci.2012.11.041

Gonzalez, J. M.; Miquel Canet, J.; Oller, S.; Miro, R. 2007. A viscoplastic constitutive model with strain rate variables for asphalt mixtures--numerical simulation, Computational Materials Science 38(4): 543-560. http://dx.doi/org/10.1016/j.commatsci.2006.03.013

Guarin, A.; Roque, R.; Kim, S.; Sirin, O. 2013. Disruption factor of asphalt mixtures, International Journal of Pavement Engineering 14(5-6): 472-485. http://dx.doi.org/10.1080/10298436.2012.727992.

Hahn, M.; Wallmersperger, T.; Kroplin, B. 2010. Discrete element representation of continua: proof of concept and determination of the material parameters, Computational Materials Science 50(2): 391-I02. http://dx.doi/org/10.1016/j.commatsci.2010.08.031

Huang, C.-W.; Abu Al-Rab, R. K.; Masad, E. A.; Little, D. N. 2011. Three-dimensional simulations of asphalt pavement permanent deformation using a nonlinear viscoelastic and viscoplastic model, Journal of Materials in Civil Engineering 23: 56-68. http://dx.doi.org/10.1061/(ASCE)MT.1943-5533.0000022

Huurman, M.; Woldekidan, M. F. 2007. LOT, Mortar response; measurements, test interpretation and determination of model parameters. Report 7-07-170-3. December 2007. 167 p.

Ka?ianauskas, R.; Rojek, J.; Pilkavi?ius, S.; Rimsa, V. 2013. Interaction of particles via solid interface, in Proc. of the 3rd International Conference on Particle-Based Methods. Fundamentals and Applications (Particles 2013), 18-20 September, 2013, Stuttgart, Germany, 1-10.

Kandhal, P. S.; Chakraborty, S. 1996. Effect of asphalt film thickness on short and long term aging of asphalt paving. Paper presented at the 75th Annual Meeting of the Transportation Research Board, 7-11 January, 1996, Washington, DC. 16 p.

Koroliov, U. 1986. Ways to save bitumen in road construction. Moscow: Transport. 138 p. (in Russian).

Kim, H.; Buttlar, W. G. 2009. Discrete fracture modeling of asphalt concrete, International Journal of Solids and Structures 46(13): 2593-2604. http://dx.doi.org/10.1016/j.ijsolstr.2009.02.006.

LST EN 12697-34: 2004. Bituminous mixtures--Test methods for hot mix asphalt--Part 34: Marshall Test. Vilnius: Lithuanian Standards Board, 2004. 15 p.

Liu, Y.; You, Z.; Dai, Q.; Mills-Beale, J. 2011. Review of advances in understanding impacts of mix composition characteristics on asphalt concrete (AC) mechanics, International Journal of Pavement Engineering 12(4): 385-405. http://dx.doi/org/10.1080/10298436.2011.575142

Liu, Y.; You, Z.; Zhao, Y. 2012. Three-dimensional discrete element modeling of asphalt concrete: size effects of elements, Construction and Building Materials 37(12): 775782. http://dx.doi.org/10.1016/j.conbuild,mat.2012.08.007

Magnanimo, V.; Huerne, H.; Luding, S. 2012. Modeling of asphalt durability and self-healing with discrete particles method, in Proc. of the 5th Eurasphalt & Eurobitumen Congress, 13-15 June, 2012, Istanbul, Turkey. 10 p.

MatWeb. 2011 [online], [cited 04 Feb. 2011]. Available from Internet: http://www.matweb.com

Meegoda, N.; Chang, K. 1993. A novel approach to develop a performance based test for rutting of asphalt concrete, in Proc. of Conference "Airport pavement innovations. Theory to Practice", 8-10 September, 1993, Vicksburg, Mississippi, USA, 126-140.

Mitra, K.; Das, A.; Basu, S. 2012. Mechanical behavior of asphalt mix: an experimental and numerical study, Construction and Building Materials 27(1): 545-552. http://dx.doi.org/10.1016/j.conbuildmat.2011.07.009

Mo, L. T.; Huurman, M.; Wu, S. P.; Molenaar, A. A. A. 2008. 2D and 3D meso-scale finite element models for ravelling analysis of porous asphalt concrete, Finite Elements in Analysis and Design 44(4): 186-196.

Ozgan, E.; Serin, S.; Kap, T. 2013. Multi-faceted investigation into the effects of hot-mix asphalt parameters on Marshall Stability, Construction and Building Materials 40: 419425. http://dx.doi.org/10.1016/j.conbulidmat 2012.11.002

Rimsa, V.; Kacianauskas, R.; Sivilevicius, H. 2011. Finite element simulation of the normal interaction of particles in visco-elastic solid, Mechanics and Control 30(4): 245253.

Roberts, F. L.; Kandhal, P. S.; Brown, E. R.; Lee, D. Y.; Kennedy, T. W. 1996. Hot mix asphalt materials, mixture design, and construction. Asphalt Pavement Association's Research and Education Foundation, Lanham, MD. 490 p.

Schapery, R. 2000. Nonlinear viscoelastic solids, International Journal of Solids and Structures 37(1-2): 359-366.

Shashidhar, N.; Shenoy, A. 2002. On using micromechanical models to describe dynamic mechanical behavior of asphalt mastics, Mechanics of Materials 34(10): 657-669. http://dx.doi.org/10.1016/S0167-6636(02)00166-7

Sivilevicius, H.; Vislavicius, K. 2008. Stochastic simulation of the influence of variation of mineral material grading and lose weight on the homogeneity of hot-mix asphalt, Construction and Building Materials 22(9): 2007-2014. http://dx.doi.org/10.1016/j.conbildmat 2007.07.001.

Sun, L.; Zhu, Y. 2013. A serial two-stage viscoelastic-viscoplastic constitutive model with thermodynamical consistency for characterizing time-dependent deformation behavior of asphalt concrete mixtures, Construction and Building Materials 40: 584-595. http://dx.doi.org/10.1016/j.conbuildmat.2012.10.004

Underwood, B. S.; Kim, Y. R. 2013. Microstructural investigation of asphalt concrete for performing multiscale experimental studies, International Journal of Pavement Engineering 14(5-6): 498-516. http://dx.doi.org/10.1080/10298436.2012.746689.

Vislavi?ius, K.; Sivilevi?ius, H. 2013. Effect of reclaimed asphalt pavement gradation variation on the homogeneity of recycled hot-mix asphalt, Archives of Civil and Mechanical Engineering 13(3): 345-353. http://dx.doi.org/10.1016/j.acme.2013.03.0003.

White, I. 1985. Marshall procedures for design and quality control of asphalt mixtures, in Proc. of Asphalt Paving Technology, 1985, Vol. 5: 265-284.

Xia, K.; Pan, T. 2011. Understanding vibratory asphalt compaction by numerical simulation, International Journal of Pavement Research & Technology 4: 185-193. http://dx.doi.org/10.1061/(ASCE)EM.1943-7889.0000730

Zhu, X.; Yang, Z.; Guo, X.; Chen, W. 2011. Modulus prediction of asphalt concrete with imperfect bonding between aggregate-asphalt mastic, Composites Part B: Engineering 42(6): 1404-1411. http://dx.doi.org/10.1016/j.compositesb.2011.05.023

Vytautas RIMSA (a), Rimantas KACIANAUSKAS (b), Henrikas SIVILEVICIUS (c)

(a) Department of Strength of Material, Vilnius Gediminas Technical University, Sauletekio al. 11, 10223 Vilnius, Lithuania

(b) Institute of Mechanics, Vilnius Gediminas Technical University, J. Basanaviciaus g. 28, 10223 Vilnius, Lithuania

(c) Department of Transport Technological Equipment, Vilnius Gediminas Technical University, Plytines g. 27, 10223 Vilnius, Lithuania Received 02 Feb 2014; accepted 18 Mar 2014

Corresponding author: Vytautas Rimsa

E-mail: vrimsa@vgtu.lt

Vytautas RIMSA. PhD student from 2010 at Department of Strength of Materials, Vilniaus Gediminas Technical University. Field of research modeling of particulate solid composite.

Rimantas KA?IANAUSKAS. Doctor, Associate Professor, he is the Head of the Institute of Mechanical science at Vilnius Gediminas Technical University. Research interests: computational mechanics, finite element method, discrete element method, structural dynamics, mechanics of materials, fracture mechanics, and coupled problems.

Henrikas SIVILEVICIUS. Doctor, Professor at Transport Technological Equipment Department, Vilniaus Gediminas Technical University. Member of Lithuanian Scientific Society (LSS), chairman at LSS department "Technika". Author of more than 150 scientific articles. Field of research: modeling of hot mix asphalt (HMA) production quality control methods improvement, HMA composition optimization and control statistical methods also asphalt pavement recycling technologies development.

Table 1. Time-dependent properties of viscoelastic

Component M   Factor                 Time [t.sub.m]
              [[alpha].sup.E.sub.m]

1             0.905                  0.113
2             0.082                  2.960

Table 2. Total normal stiffness of the system of two pa

Thickness [delta]          Total stiffness [K.sub.tot] in N/m
of the film in ([mecto]m

0 (no film)                7.15 x 1[0.sup.4]
4                          7.18 x [10.sup.4]
12                         7.25 x [10.sup.4]
50                         7.59 x [10.sup.4]

Table 3. Stiffness for given models

[[alpha].sub.1]   71.92[degrees]
[[alpha].sub.4]   73.09[degrees]
[[alpha].sub.M]   73.79[degrees]

Table 4. Force reaction of numerical Marshall test

[E.sub.f];   0 (no film)   4.[10.sup.6]   4.[10sup.7]   4.[10.sup.8]
  pa
F kN         10.99         11.11          11.73         11.75