A time domain boundary element method for modeling the quasi-static viscoelastic behavior of asphalt pavements

ARTICLE IN PRESS

Engineering Analysis with Boundary Elements 31 (2007) 226–240 www.elsevier.com/locate/enganabound

A time domain boundary element method for modeling the quasi-static viscoelastic behavior of asphalt pavements Jianlin Wang, Bjorn Birgisson Department of Civil and Coastal Engineering, University of Florida, 365 Weil Hall, P.O. Box 116580, Gainesville, FL 32611, USA Received 6 March 2006; accepted 5 September 2006 Available online 13 November 2006

Abstract A time domain boundary element method (BEM) is presented to model the quasi-static linear viscoelastic behavior of asphalt pavements. In the viscoelastic analysis, the fundamental solution is derived in terms of elemental displacement discontinuities (DDs) and a boundary integral equation is formulated in the time domain. The unknown DDs are assumed to vary quadratically in the spatial domain and to vary linearly in the time domain. The equation is then solved incrementally through the whole time history using an explicit time-marching approach. All the spatial and temporal integrations can be performed analytically, which guarantees the accuracy of the method and the stability of the numerical procedure. Several viscoelastic models such as Boltzmann, Burgers, and power-law models are considered to characterize the time-dependent behavior of linear viscoelastic materials. The numerical method is applied to study the load-induced stress redistribution and its effects on the cracking performance of asphalt pavements. Some benchmark problems are solved to verify the accuracy and efficiency of the approach. Numerical experiments are also carried out to demonstrate application of the method in pavement engineering. r 2006 Elsevier Ltd. All rights reserved. Keywords: Time domain BEM; Viscoelasticity; Displacement discontinuities; Time marching; Viscoelastic models; Asphalt pavements; Stress redistribution

1. Introduction The loading response of asphalt pavements is, in general, time dependent because of the rheological behavior of the hot mix asphalt (HMA) mixtures. Ignoring this fact, linearelastic pavement models always predict the highest tensile stresses at the bottom of the pavement, which suggests that fatigue cracks initiate and propagate upwards, called bottom-up cracking. In many cases, however, the cracks are observed to initiate and propagate downwards, called top-down cracking, which has been accepted as one of the primary deterioration mechanisms of flexible pavements [1–4]. Although the viscoelastic behavior of asphalt mixtures has been well recognized, it has rarely been taken into account in direct simulation of the time-dependent behavior and cracking performance of asphalt pavements. Corresponding author. Tel.: +1 352 392 9537x1480; fax: +1 352 392 3394. E-mail address: [email protected]fl.edu (J. Wang).

0955-7997/$ - see front matter r 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2006.09.007

Due to its versatility and wide spread in commercial packages, the finite element method (FEM) has been the primary tool used for pavement modeling. One well-known disadvantage associated with the FEM is the need for domain discretization, which makes the method relatively expensive to set up and use, especially for simulation of pavement cracking, which typically requires elaborate remeshing in FEM to update the geometry when a crack grows. Compared to the FEM, the boundary element method (BEM) reduces the dimension of the problem by one and provides an attractive alternative to the FEMs for elasticity and viscoelasticity problems. Birgisson et al. [5,6] have recently applied a special boundary element-based method, called the displacement discontinuity (DD) method, to pavement modeling and have proven the superiority of the method over the traditional FEM for the particular problem in consideration. The DD method, originally developed by Crouch [7], solves for stresses and displacements within an elastic solid in terms of DDs (dislocations), which makes it uniquely suitable for crack

ARTICLE IN PRESS J. Wang, B. Birgisson / Engineering Analysis with Boundary Elements 31 (2007) 226–240

problems [8–10]. This method has been refined over years to include new features such as viscoelasticity [11], elastodynamics [12,13], and high-order elements [14] and has been successfully applied in several engineering fields such as mining and oil industry [15–17] as well as pavement engineering [5,6,18,19]. This paper describes the development of a viscoelastic DD method for efficiently modeling the quasi-static responses of linear viscoelastic materials. Particular consideration is placed on the implementation of the method to investigate the time-dependent response of bituminous materials for better understanding of the cracking mechanism in asphalt pavements. There are basically three approaches that can be used in linear viscoelastic analysis: Laplace transformation, Fourier transformation, and direct time integration method [20]. The first two approaches involve solving an associated elastic problem in the transformed domain and inverting the result to obtain the actual solution in the time domain [11,20–22]. The use of transformation offers a simple and straightforward implementation and makes it easy to handle the timedependent boundary conditions. In the Fourier domain analysis, however, it must be performed over a time duration and the computational storages are required at all frequencies [20]. The Laplace transform seems to be more efficient in computational storage, but the Laplace inversion cannot be obtained analytically for nontrivial problems in general and the numerical inversion has been a difficult task for most practical implementations [11]. In the time integration approach, the boundary integral equation is formulated directly in the time domain and solved incrementally through the whole history with a stepby-step time integration scheme [23]. In the present work, the time domain approach is adopted to avoid numerical Laplace inversion. The approach presented in this paper combines the DD method with an explicit time-marching procedure [24] for linear viscoelastic analysis. The time-dependent behavior of the material is characterized by one of the three commonly used viscoelastic models—Boltzmann, Burgers, or powerlaw model. A substructuring approach is employed for solving problems associated with piecewise homogeneous materials or structures, for example, fiber-reinforced composite materials or layered asphalt pavements. Quadratic functions are used to approximate the DDs on the boundaries, which, in the time domain, are approximated by linear functions. At each time step, all the spatial integrations can be performed analytically and a collocation formulation is used to obtain a linear algebraic system for solving the unknown DDs. In the time-marching process, the boundary unknowns for a given time step are found from the solutions obtained at all previous steps. The integrations over each time step can also be performed analytically. Several numerical examples are presented to illustrate the versatility and accuracy of the approach. These examples include some benchmark problems such as a

227

pressured crack or circular hole in an infinite viscoelastic plane that have analytical solutions, as well as the problem of circular holes and elastic inclusions in a viscoelastic plane, which has been solved numerically by Huang et al. [25] using two different approaches and has demonstrated its application in composite materials. The solutions to these problems are used as independent checks for the current approach presented. Finally, focus will be placed on implementation of the approach to investigate the timedependent behavior of viscoelastic asphalt materials. Detailed analyses are performed to study the load-induced stress redistribution and its effects on the cracking performance of asphalt pavements.

2. Fundamentals of linear viscoelasticity 2.1. Constitutive equations The stress–strain relation for a linear viscoelastic material can be expressed in an integral form by the following two equations [26]: Z t qekk ðx; tÞ skk ðx; tÞ ¼ 3 dt, Kðt  tÞ qt 0 Z t qeij ðx; tÞ dt, ð1Þ sij ðx; tÞ ¼ 2 Gðt  tÞ qt 0 where x is the position vector, KðtÞ and GðtÞ are the relaxation functions of bulk modulus and shear modulus, skk and ekk are the volumetric (or dilatational) stress and strain, respectively, and sij and eij are the deviatoric components of the stress and strain tensors: sij ¼ sij  13 dij skk ;

eij ¼ eij  13 dij ekk .

ð2Þ

The constitutive equations for a linear viscoelastic material can alternatively be written in a differential form as [11,26] Pv ðLÞskk ðx; tÞ ¼ Qv ðLÞekk ðx; tÞ, Ps ðLÞsij ðx; tÞ ¼ Qs ðLÞeij ðx; tÞ,

ð3Þ

where L ¼ q=qt and Pv ðLÞ, Qv ðLÞ, Ps ðLÞ, Qs ðLÞ are time differential operators. By the correspondence principle [27], the viscoelastic governing equations can be transformed into a set of corresponding elastic governing equations using the Laplace transform. As a result, the transformed constitutive equations corresponding to (1) and (3) in the Laplace domain are written as skk ðx; sÞ ¼ 3KðsÞekk ðx; sÞ, sij ðx; sÞ ¼ 2GðsÞeij ðx; sÞ

ð4Þ

and Pv ðsÞskk ðx; sÞ ¼ Qv ðsÞekk ðx; sÞ, Ps ðsÞsij ðx; sÞ ¼ Qs ðsÞeij ðx; sÞ,

ð5Þ

ARTICLE IN PRESS J. Wang, B. Birgisson / Engineering Analysis with Boundary Elements 31 (2007) 226–240

228

respectively, where s is the Laplace transform parameter. Comparison of (4) and (5) yields GðsÞ ¼

Qs ðsÞ ; 2Ps ðsÞ

KðsÞ ¼

Qv ðsÞ : 3Pv ðsÞ

ð6Þ

The viscoelastic Young’s modulus EðsÞ and Poisson’s ratio nðsÞ in the transformed domain can accordingly be expressed as EðsÞ ¼

9KðsÞGðsÞ ; 3KðsÞ þ GðsÞ

nðsÞ ¼

3KðsÞ  2GðsÞ : 6KðsÞ þ 2GðsÞ

ð7Þ

The transformed material constants expressed in (6) and (7) are only used to derive the viscoelastic fundamental solutions. In viscoelastic analysis, for many cases it can be assumed that the material behaves elastically in dilatation [24]. As a result, we can write KðtÞ ¼ KðsÞ ¼ K,

(8)

where K is a constant. In this case, the Poisson’s ratio is time dependent in general and can be determined from (7). In practical implementations, we sometimes assume that the viscoelastic material behaves the same in dilatation as in shear, which results in a constant Poisson’s ratio n. This assumption is made only for simplicity and could reduce the complexity of the mathematical derivations. For either case, we only need to determine the expression for the shear modulus GðsÞ. In this paper, both cases (constant K and constant n) are considered to solve practical problems. 2.2. Linear viscoelastic models Fig. 1 shows the one-dimensional representation of the Boltzmann and Burgers models. The Boltzmann model (Fig. 1(a)) consists of a spring and Kelvin model connected in series, where the spring is used to represent the instantaneous elastic response (eie ) and the Kelvin model to represent the delayed elastic response (ede ). The Burgers model (Fig. 1(b)) is represented by the arrangement of a Kelvin model and a Maxwell model connected in series. Compared to the Boltzmann model, the additional dashpot 

 R1

ie

ie

R1

Maxwell

1 

R2

de Kelvin

 (a)

cp

2

R2

de

 (b)

Fig. 1. One-dimensional representation of (a) Boltzmann model and (b) Burgers model.

in Burgers model is used to represent the creep behavior of the viscoelastic material, and the creep strain is denoted by ecp as illustrated in Fig. 1. In both models, R1 and R2 denote the spring constants, Z, Z1 , and Z2 denote the viscosity coefficients of dashpots. The one-dimensional stress–strain relations for Boltzmann model and Burgers model can be written as (9) and (10), respectively:     d d R1 þ R2 þ Z sðtÞ ¼ R1 R2 þ Z eðtÞ, ð9Þ dt dt     d d2 d d2 1 þ p1 þ p2 2 sðtÞ ¼ q1 þ q2 2 eðtÞ, ð10Þ dt dt dt dt where eðtÞ ¼ eie ðtÞ þ ede ðtÞ in (9) and eðtÞ ¼ eie ðtÞ þ ecp ðtÞ þ ede ðtÞ in (10); and for the Burgers model, Z Z Z ZZ p1 ¼ 1 þ 1 þ 2 ; p2 ¼ 1 2 , R1 R2 R2 R1 R2 Z1 Z2 . ð11Þ q1 ¼ Z1 ; q2 ¼ R2 Extending the above two models into three dimensions, we can obtain the transformed shear modulus GðsÞ as follows by taking into account (5) and (6): R1 R2 þ R1 Zs for Boltzmann model, 2ðR1 þ R2 þ ZsÞ q1 s þ q2 s 2 for Burgers model. GðsÞ ¼ 2ð1 þ p1 s þ p2 s2 Þ GðsÞ ¼

ð12Þ ð13Þ

The Poisson’s ratio n will either be assumed as a constant (n ¼ n) or be obtained from (7), where G is given in (12) or (13) and K is assumed to be a constant. Another model that is commonly used to characterize the viscoelastic behavior of asphalt mixtures is also considered. This model expresses the creep compliance JðtÞ directly as a power function of time as JðtÞ ¼ l0 þ l1 tm ,

(14)

where l0 , l1 , and m are material constants (for asphalt mixtures, mo1). The l0 constant in the power-law model describes the immediate elastic response; and the l1 tm term represents the time-dependent response. The power-law equation for compliance is usually obtained by applying a nonlinear curve-fitting method to the graphs of compliance versus time that are obtained from measurements in laboratory tests [28,29]. 3. Basic equations In a direct boundary integral formulation, the displacements and stresses in a viscoelastic body can be written in a time integral form of the Somigliana’s formula for elastic problems by using the reciprocal theorem and the Stieltjes convolution [24,30]. Other than the direct approach, an indirect formulation of the DD method is also available through the use of an influence function technique. In this version of the method, the solution to the problem of a DD over a line segment (in two dimensions) is used as the influence functions.

ARTICLE IN PRESS J. Wang, B. Birgisson / Engineering Analysis with Boundary Elements 31 (2007) 226–240

As a simple illustration, consider a DD over the line segment apxpa, y ¼ 0: Di ðxÞ ¼ ui ðx; 0þ Þ  ui ðx; 0 Þ

ði ¼ x; yÞ.

(15)

The sign convention of the DD in (15) is shown in Fig. 2. For the case of plane strain, the stress and displacement solution to the problem of a DD ðDx ; Dy Þ defined in (15) can be expressed as [8,14] ux ¼ 2cx  k2 fy  k1 yðfx;x þ cy;x Þ, uy ¼ 2cy þ k2 fx  k1 yðfx;y þ cy;y Þ,

229

As an example, we calculate K3 ðtÞ for Boltzmann-type material that is elastic in dilatation (bulk modulus K is a constant). Substitution of (7) and (12) into (19) yields k3 ðsÞ ¼

6KR1 ðR2 þ ZsÞðR1 þ R2 þ ZsÞ þ R21 ðR2 þ ZsÞ2 , 3KðR1 þ R2 þ ZsÞ2 þ 2R1 ðR2 þ ZsÞðR1 þ R2 þ ZsÞ (21)

where K, R1 , R2 , and Z are model constants as described earlier. The function K3 ðtÞ is finally obtained by applying inverse Laplace transform to (21):

sxx ¼ k3 ½cy;y þ 2fx;y þ yðfx;yy þ cy;yy Þ,

K3 ðtÞ ¼ a1  dðtÞ  a2  expða3 tÞ  a4  expða5 tÞ,

syy ¼ k3 ½cy;y  yðfx;yy þ cy;yy Þ,

where dðtÞ is the Dirac’s delta function and the coefficients are defined as

sxy ¼ k3 ½cx;y  yðfy;yy þ cx;xx Þ,

ð16Þ

where functions fi and ci are defined as Z 1 a xx fi ðx; yÞ ¼ Di ðxÞ dx 4p a ðx  xÞ2 þ y2 Z ði ¼ x; yÞ, 1 a y ci ðx; yÞ ¼ Di ðxÞ dx 4p a ðx  xÞ2 þ y2

a2 (17) a3

1 1  2n 2G ; k2 ¼ ; k3 ¼ . (18) 1n 1n 1n For an arbitrarily oriented DD, a coordinate transformation needs to be performed to change solution (17) from a local coordinate system to a global coordinate system [8]. The time-dependent solution for the problem of a DD ðDx ; Dy Þ in a viscoelastic plane can be obtained from (16) using the correspondence principle. First, the transformed solution in the Laplace domain is obtained by replacing the constants ki (i ¼ 1; 2; 3) in (16) with ki ðsÞ=s, where ki ðsÞ are obtained from (18) by using the transformed material properties GðsÞ and nðsÞ that are derived in the previous section: k1 ¼

1 ; 1  nðsÞ

k2 ðsÞ ¼

1  2nðsÞ ; 1  nðsÞ

k3 ðsÞ ¼

2GðsÞ . 1  nðsÞ (19)

Then, an inverse Laplace transform is applied to obtain the viscoelastic version of the fundamental solution expressed in (16). In this procedure, there are only three basic time functions that need to be obtained from inverse Laplace transform: K1 ðtÞ ¼ L1 ½k1 ðsÞ;

K2 ðtÞ ¼ L1 ½k2 ðsÞ,

K3 ðtÞ ¼ L1 ½k3 ðsÞ. + Dn

ð20Þ

+ Ds

Fig. 2. Normal and shear displacement discontinuities.

R1 ð6K þ R1 Þ , 3K þ 2R1 R2 ¼ 1, 2Z R1 þ R2 , ¼ Z 27K 2 R21 ¼ , 2Zð3K þ 2R1 Þ2 2R1 R2 þ 3KðR1 þ R2 Þ . ¼ Zð3K þ 2R1 Þ

a1 ¼

and k1 , k2 , and k3 are introduced for notational convenience:

k1 ðsÞ ¼

(22)

a4 a5

The other functions in (20) can be obtained in a similar way. If the Poisson’s ratio is assumed to be a constant, calculation of the functions in (20) will become much simpler. It is worth mentioning that the fundamental solutions obtained from the indirect formulation are essentially the same as those from the direct boundary integral approach. 4. Numerical procedure 4.1. Approximation of unknown functions In either the direct or indirect formulation, the boundaries of the solution domain are first divided into a number of (say, M) boundary elements. The unknown DD over each boundary element can be approximated locally by a polynomial, which is fitted to a specified number of collocation points along the element. In the present analysis, second-order polynomials are used for approximation of unknowns because the quadratic DD elements generally give better overall results than the constant or linear elements and are able to capture the bending effects that are commonly present in flexible pavements. In addition, the use of quadratic elements allows for more accurate determination of the stresses at locations close to the boundaries [31]. As an illustration, Fig. 3 depicts a quadratically varying DD along the line segment apxpa; y ¼ 0. The DD at any point x ðapxpaÞ along the segment can be expressed in terms of the values of the DD components at three collocation points x1 , x2 , and x3 . In general, these

ARTICLE IN PRESS J. Wang, B. Birgisson / Engineering Analysis with Boundary Elements 31 (2007) 226–240

230

foregoing discussions, the displacements at a collocation point i can be written as Z t ½Aijss ðtÞDjs ðt  tÞ þ Aijsn ðtÞDjn ðt  tÞ dt, uis ðtÞ ¼ 0 Z t ½Aijns ðtÞDjs ðt  tÞ þ Aijnn ðtÞDjn ðt  tÞ dt, ð26Þ uin ðtÞ ¼

y

(Di)2

(Di)1

(Di)3

0

uis

x −a

x1

x2

x3

a

Fig. 3. A quadratically varying displacement discontinuity.

points can be arbitrarily chosen within the element except the element end points, where numerical problems may occur due to the interpolation function not being differentiable across elements. Crawford and Curran [31] suggested that the optimum locations of the nodes are at the points, i.e. x1 ¼ pffiffiffi Gauss–Chebyshev integration pffiffiffi  3a=2, x2 ¼ 0, and x3 ¼ 3a=2 for a second-order polynomial. If we set Di ðxk Þ ¼ ðDi Þk ðk ¼ 1; 2; 3Þ, the expression for Di ðxÞ can be written as Di ðxÞ ¼ a1 ðxÞðDi Þ1 þ a2 ðxÞðDi Þ2 þ a3 ðxÞðDi Þ3 ;

i ¼ x; y, (23)

where 1 x 2 x2 a1 ðxÞ ¼  pffiffiffi þ ; 3a 3 a 1 x 2  x 2 a3 ðxÞ ¼ pffiffiffi þ . 3a 3 a

a2 ðxÞ ¼ 1 

4 x2 , 3 a

where and uin are the shear and normal displacements at the collocation point i, Djs and Djn are shear and normal components of the DD at point j, and Aijss , Aijsn , etc. are coefficients that express the influence of point j on point i. As mentioned previously, the repeated index j in the equation implies summation. Similarly, the shear and normal boundary tractions ss and sn at point i can be written as Z t i ½Bijss ðtÞDjs ðt  tÞ þ Bijsn ðtÞDjn ðt  tÞ dt, ss ðtÞ ¼ 0 Z t i ½Bijns ðtÞDjs ðt  tÞ þ Bijnn ðtÞDjn ðt  tÞ dt. ð27Þ sn ðtÞ ¼ 0

Explicit expressions for the influence coefficients Aij ðtÞ and Bij ðtÞ in (26) and (27) can be obtained by using the functions Ki ðtÞ (i ¼ 1; 2; 3) in (20) and the integration results given in Appendix A, and applying suitable coordinate transformations, but the details are omitted here. 4.3. Time-marching scheme

ð24Þ

4.2. Evaluation of spatial integrals By using the approximation in (23), all the spatial integrals involved in (16) can be evaluated analytically. The results can be written explicitly in terms of the following two types of functions:  k Z a xx x Ik ðx; yÞ ¼ dx 2 2 a a ðx  xÞ þ y ðk ¼ 0; 1; 2Þ, (25)  k Z a y x dx Jk ðx; yÞ ¼ 2 2 a a ðx  xÞ þ y and their derivatives, which are summarized in Appendix A. Once the influence of one single element is calculated, the combined effects of all M DDs at a point in the solution domain are found by summing the contributions from the individual elements, provided suitable coordination transformation is applied to each element to account for its position and orientation (see [8] for details). Since it is convenient to define the boundary conditions in the local coordinate system, the displacements and tractions at the collocation points are expressed in terms of the shear and normal DDs at these points. For example, according to the

The next step is to evaluate the time integrations involved in (26) and (27). Following Sim and Kwak [24], a time-marching scheme is used to evaluate the convolution integrals numerically. However, we use a linear time marching instead of a constant marching scheme to achieve better accuracy and efficiency. For the purpose of numerical calculation, the time interval ½0; t is divided into N time steps, denoted as t1 ; t2 ; . . . ; tN ðtN ¼ tÞ. Then, all the time integrals in (26) and (27) can be written in a piecewise fashion as Z t H ij ðtÞDj ðt  tÞ dt 0 Z t1 Z t2 Z tN  þ þ þ ð28Þ ¼ H ij ðtÞDj ðt  tÞ dt, 0

t1

tN1

where Dj ð ¼ s; nÞ denotes the shear or normal DD at j collocation point, and H ij stands for one of the influence coefficients Aij or Bij . The DDs Dj ð ¼ s; nÞ over each time step ½tn1 ; tn  ðn ¼ 1; 2; . . . ; NÞ are approximated by linear functions as ðn1Þ þ g2n ðtÞDðnÞ Dj ðtÞ ¼ g1n ðtÞDj j ; ðn1Þ Dj

tn1 ptptn ,

(29)

DðnÞ j

and are used to denote Dj ðtn1 Þ and where Dj ðtn Þ for simplicity, and the functions g1n and g2n are defined as g1n ðtÞ ¼ ðtn  tÞ=DT n ;

g2n ðtÞ ¼ ðt  tn1 Þ=DT n ,

(30)

ARTICLE IN PRESS J. Wang, B. Birgisson / Engineering Analysis with Boundary Elements 31 (2007) 226–240

where DT n ¼ tn  tn1 is the length for time step n. The linear approximation is illustrated in Fig. 4. It is worth mentioning that DT n ðn ¼ 1; . . . ; NÞ can generally vary from one step to another, but we assume DT n ¼ DT to be a constant to simplify the derivations. Using the approximation in (29), the individual integrals in (28) can be written as Z tn Z tn H ij ðtÞDj ðt  tÞ dt ¼ Dðn1Þ g1n ðtÞH ij ðt  tÞ dt j tn1

þ DðnÞ j

Z

g2n ðtÞH ij ðt  tÞ dt

ðn ¼ 1; . . . ; NÞ.

ð31Þ

tn1

According to previous discussions, the influence coefficients H ij ðtÞ can be expressed as a linear combination of the time functions K| ðtÞ ð| ¼ 1; 2; 3Þ in (20). Thus, only the following time integrals need to be evaluated: Z tn ðnÞ G‘| ¼ g‘n ðtÞK| ðt  tÞ dt ð‘ ¼ 1; 2; | ¼ 1; 2; 3Þ. (32) tn1

Taking into account the definition of g‘n ðtÞ given in (30), GðnÞ ‘| can be written as  Z tn  Z tn ð1Þ‘ GðnÞ ¼ tK ðt  tÞ dt  t K ðt  tÞ dt | n‘þ1 | ‘| DT tn1 tn1 ‘

ð1Þ ½tNnþ‘1  PðNnþ1Þ  PðNnþ1Þ , 0| 1| DT ðnÞ where PðnÞ 0| and P1| are defined as ¼

ð33Þ

ðNÞ ðNÞ þ Að1Þ Að1Þ ss Ds sn Dn

ðNÞ Að1Þ þ ns Ds

ð34Þ

D(t) D(tn-1)

D(tn)

D(t2) D(t1) D(t0)

t1

t2

tn-1

Fig. 4. A linear time-marching scheme.

¼ uðNÞ n 

N 1 X

þ AðNkþ1Þ DðkÞ sn n Þ,

ðNkþ1Þ ðkÞ ðAðNkþ1Þ DðkÞ Dn Þ, ns s þ Ann

ð35Þ

k¼1

where Ds and Dn are vectors of the shear and normal DDs at the collocation points, the superscripts in the parenthesis denote the time steps, AðkÞ is the displacement influence matrix obtained from integration of the time-dependent coefficients Aij in (26) over time step k, and uðNÞ and uðNÞ s n are vectors of the prescribed shear and normal displacements at the current time. A similar linear algebraic system can be written for traction boundary conditions as ðNÞ ðNÞ þ Bð1Þ Bð1Þ ss Ds sn Dn

¼ sðNÞ  s ðNÞ Bð1Þ ns Ds

þ

N 1 X

ðNkþ1Þ ðkÞ ðBðNkþ1Þ DðkÞ Dn Þ, ss s þ Bsn

k¼1 ð1Þ ðNÞ Bnn Dn N 1 X

ðNkþ1Þ ðkÞ ðBðNkþ1Þ DðkÞ Dn Þ, ns s þ Bnn

ð36Þ

k¼1

The detailed expressions for Pm| ðtÞ are given in Appendix B for all viscoelastic models considered in the paper under the assumption that the Poisson’s ratio does not change with time, i.e. the material behaves in the same manner in shear as in dilatation. This assumption is made in some practical applications for simplicity [25,32]. For the case where the bulk modulus is a constant, i.e. the material is elastic in dilatation, the function Pm| ðtÞ can be found in a similar way, but the expressions are omitted here. After evaluating all the integrals involved in (31), the results are substituted into (28) and then into (26) or (27) to obtain a linear system of algebraic equations, which are

...

N 1 X

ðAðNkþ1Þ DðkÞ ss s k¼1 ðNÞ Að1Þ nn Dn

¼ uðNÞ  s

¼ sðNÞ n 

PðnÞ m| ¼ Pm| ðtn Þ  Pm| ðtn1 Þ and Z Pm| ðtÞ ¼ tm K| ðtÞ dt ðm ¼ 0; 1; | ¼ 1; 2; 3Þ.

0

used to solve for the unknown DDs step by step. If the boundary condition is displacement prescribed, the linear system is given as follows:

tn1

tn

231

tn

time

where BðkÞ is the traction influence matrix obtained from integration of the coefficients Bij in (27), and sðNÞ and snðNÞ s are vectors of the prescribed shear and normal boundary tractions at the current time. In case of mixed boundary conditions, a combination of (35) and (36) will be used. It is observed from (35) and (36) that the unknown DDs at the current time step depend on their values at previous time steps, which are obtained earlier in the time-marching process. 4.4. A substructuring approach for piecewise homogeneous bodies The viscoelastic DD method can be formulated to model layered structures or more general piecewise homogeneous bodies by simply rearranging the formulation described above. As described by Crouch and Starfield [8], the piecewise homogeneous domain is first divided into a number of homogeneous subregions, each of which can be elastic or viscoelastic. A substructuring approach is then employed to satisfy the continuity condition along the interfaces—the displacements on either side of a bonded interface must be equal, and the tractions must be equilibrated. Unlike the regular boundary elements, where boundary values for tractions or displacements are prescribed, the boundary parameters at the interface are unknown but are continuous across the boundary. Imperfect interface that allows for traction or displacement

ARTICLE IN PRESS J. Wang, B. Birgisson / Engineering Analysis with Boundary Elements 31 (2007) 226–240

232

jumps can also be handled in a similar way, but it is beyond the scope of this paper. Once the DDs are solved for the whole solution domain, the stresses and displacements in the individual subregions can be calculated as in the case of a homogeneous body. The substructuring approach is commonly used to solve inclusion problems with application in composite materials and to model multilayered structures such as pavements.

where HðtÞ is the heaviside step function and the coefficients ci (i ¼ 1; . . . ; 5) are defined as

c2 c3 c4

5. Verification through benchmark problems

c5

In order to verify the proposed numerical procedure, several benchmark problems under plane strain conditions are solved using the current approach. The numerical results are compared with the analytical solutions if available, or with other numerical solutions obtained from independent approaches.

(37)

where G and n are the elastic constants, and Du denotes the crack opening displacement, which is equal to the positive normal DD according to the sign convention defined in Fig. 2. Using the correspondence principle and analytical Laplace transform inversion, the viscoelastic solution can be expressed as pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Duðx; tÞ ¼ 4p b2  x2 ½c1  c2 expðc3 tÞ  c4 expðc5 tÞ  HðtÞ,

ð38Þ

y p −b

K ¼ 2  104 ;

R1 ¼ 2  104 ;

R2 ¼ 104 MPa,

Z ¼ 5  106 MPa h.

Duðx; tÞ ¼ Duðx; tÞ  K=ðp  bÞ;

t ¼ t  K=Z.

The numerical solution is obtained by using 21 quadratic DD elements to approximate the crack and using 20 steps in the time-marching process to arrive at t ¼ 2:4 ðt ¼ 600 hÞ. Good agreement between the numerical and analytical solutions is observed from Fig. 6. To study how the numerical accuracy can be improved by using high-order DD elements and different timemarching schemes, the same problem is solved for three different numerical scenarios: (1) constant strength DD (CSDD) elements and constant time-marching scheme; (2) CSDD elements and linear time-marching scheme; (3) high-order quadratic DD elements (HODD) and linear 6 Normalized crack opening displacement

As the first example, consider the problem of a pressurized crack in a viscoelastic plane of the Boltzmann type. For simplicity, the crack is assumed to be aligned along x-axis with the center located at the origin of the coordinate system, as depicted in Fig. 5. The length of the crack is 2b and the crack surfaces are subjected to uniform normal pressure p at t ¼ 0. Suppose that the material is elastic in dilatation and that the material constants K, R1 , R2 , and Z (see Fig. 1(a)) are given. For the case where the plane is elastic, the analytical solution for the crack opening displacement is given as [8] 2ð1  nÞp pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi b 2  x2 , G

The analytical and numerical results for the crack opening displacement are calculated for the following parameters:

Fig. 6 shows the evolution of the normalized crack opening displacement at the center of the crack (x ¼ 0), where the normalized displacement Duðx; tÞ and the normalized time t are defined as

5.1. A pressurized crack in a viscoelastic plane

DuðxÞ ¼

ðR1 þ R2 Þ½2R1 R2 þ 3KðR1 þ R2 Þ , R1 R2 ½R1 R2 þ 6KðR1 þ R2 Þ 1 ¼ , 2R2 R2 , ¼ Z 3R21 , ¼ 2ð6K þ R1 Þ½R1 R2 þ 6KðR1 þ R2 Þ R1 R2 þ 6KðR1 þ R2 Þ . ¼ Zð6K þ R1 Þ

c1 ¼

5

Analytical solution Numerical solution

4

3

2 b

x

Fig. 5. A uniformly pressurized crack in a viscoelastic plane.

0.0

0.4

0.8

1.2 1.6 Normalized time

2.0

2.4

Fig. 6. Evolution of the normalized crack opening displacement Duð0; tÞ.

ARTICLE IN PRESS J. Wang, B. Birgisson / Engineering Analysis with Boundary Elements 31 (2007) 226–240

5.2. One circular hole in a Burgers viscoelastic plane This example considers a circular hole of radius a in an infinite Burgers viscoelastic plane (Fig. 8). The hole is loaded with constant internal pressure p at t ¼ 0. The timedependent radial displacement around the hole is to be calculated. The analytical viscoelastic solution for the radial displacement ur ðr; tÞ can be written in the following form by using the correspondence principle and analytical Laplace inversion: ur ðr; tÞ

     pa2 p1 q2 1 p p q q ¼  2 þ t þ 2  1 þ 22 exp  1 t , r q1 q 1 q1 q2 q1 q1 q2 ð39Þ

Normalized crack opening displacement

6

where p1 , p2 and q1 , q2 are constants as defined in (11). For the following Burgers parameter (see Fig. 1(b)) constants: R1 ¼ 4  103 ;

R2 ¼ 2  104 MPa;

Z1 ¼ 2  106 ,

Z2 ¼ 5  105 MPa h,

ð40Þ

the radial displacement at r ¼ 2a is calculated from both the analytical solution given in (39) and the numerical approach presented. The results for t 2 ½0; Z1 =R1  are plotted in Fig. 9, where the time is normalized by Z1 =R1 and the displacements are normalized by pa=R1 . In the numerical computation, the hole is approximated by 36 quadratic elements and the time-marching process involves 40 equal time steps. Fig. 9 shows a good agreement between the numerical solution and the analytical one.

1.2

Normalized radial displacement

time-marching scheme. In all the scenarios considered, the crack is approximated by 21 DD elements and 20 time steps are used in the marching process. The results are plotted in Fig. 7 for comparison. It is observed from the figure that the numerical accuracy can be significantly improved by using a linear instead of constant timemarching scheme and that the accuracy is further improved by using quadratic DD elements in lieu of constant DD elements.

233

5

1.0

0.8 Analytical solution Numerical solution 0.6

0.4 0.0

0.2

0.4

0.6

0.8

1.0

Normalized time

4

Fig. 9. Normalized radial displacement ur at r ¼ 2a. Analytical solution CSDD, constant time marching

∞ yy

CSDD, linear time marching

3

y

HODD, linear time marching

r4

2 0.0

0.4

0.8

1.2 1.6 Normalized time

2.0

Fig. 7. The results of Duð0; tÞ for different numerical scenarios.

p

O4

2.4

∞ xx

C•

G1, v1 O1 r1

• D r2

O2

G2, v2

x

∞ xx

r3 O3

a r = 2a

 ∞yy

Fig. 8. A uniformly loaded circular hole in a viscoelastic plane.

Fig. 10. An infinite viscoelastic plane with two holes and two elastic inclusions.

ARTICLE IN PRESS J. Wang, B. Birgisson / Engineering Analysis with Boundary Elements 31 (2007) 226–240

234

5.3. A viscoelastic plane with circular holes and elastic inclusions

0.7

Normalized vertical stress

0.6 0.5 0.4

C - Present approach D - Present approach

0.3 0.2 0.1 0.0 0.0

1.0

2.0 Normalized time

3.0

4.0

Fig. 12. Normalized stress syy =s0 at points C and D computed from the present approach.

0.5 Normalized horizontal displacement

The next example is concerned with a system of two circular holes and two circular inclusions within an infinite viscoelastic plane, as shown in Fig. 10. The inclusions are assumed to be elastic, and the matrix is a Boltzmann-type viscoelastic material of constant Poisson’s ratio. This kind of problem has been solved by Huang et al. [25] using both the finite element software—ANSYS—and a boundary integral method they developed. For the purpose of comparison, the geometry of the problem considered in this paper is taken the same as that in [25]. The two inclusions have the same radii r1 ¼ r2 ¼ a, and their centers are located at points O1 ð2a; 0Þ and O2 ð2a; 0Þ, respectively. The two holes have the radii r3 ¼ r4 ¼ 0:8a, and their centers are located at O3 ð0; 2aÞ and O4 ð0; 2aÞ. The viscoelastic plane is subjected to biaxial far-field stresses 1 s1 xx ¼ s0 and syy ¼ s0 =2 at t ¼ 0. According to Huang et al. [25], the shear modulus and Poisson’s ratios of the two elastic inclusions are chosen as G 1 ¼ 1  102 s0 , n1 ¼ 0:25 and G 2 ¼ 5  102 s0 , n2 ¼ 0:2. The Poisson’s ratio and the Boltzmann model parameters of the viscoelastic matrix are taken as n ¼ 0:3, R1 ¼ 1:6  104 s0 , R2 ¼ 4  103 s0 , and Z ¼ 1:3  104 s0 s. The problem is solved within the time interval ½0; 4g, where g is viscosity coefficient defined as g ¼ Z=½R1 ð1 þ nÞ. In finite element analysis using ANSYS, Huang et al. [25] used 3944 eight-node elements and the computation took about 4 h with an IBM SP workstation. It is worth noting that the fine finite element mesh was mainly used for comparison of the accuracy but not the efficiency [25]. In the present approach, 88 quadratic elements are used (with 24 for each inclusion and 20 for each hole) and the time step is taken as Dt ¼ g=10. It takes about 1 min to solve the problem on a personal computer.

0.4 C - ANSYS C - Present approach

0.3

D - ANSYS D - Present approach

0.2

0.1

0.0 0.0

1.0

2.0 Normalized time

3.0

4.0

Fig. 13. Normalized displacement ðux =aÞ  103 at points C and D computed from ANSYS and the present approach.

Normalized horizontal stress

1.2

The stresses and displacements are calculated at two arbitrarily selected points—C (1:139a, 0:730a) and D (1:887a, 0:900aÞ, as shown in Fig. 10. Point C is located inside the matrix and point D is inside the inclusion. The normalized stresses sxx =s0 and syy =s0 at points C and D are plotted in Figs. 11 and 12; the normalized displacements ðux =aÞ  103 and ðuy =aÞ  104 at the two points are plotted in Figs. 13 and 14. In Figs. 11 and 13, the results for the horizontal stresses and displacements at points C and D are also compared with the ANSYS results provided in [25]. Good agreement can be observed between the results obtained using two independent approaches.

1.0 C- ANSYS C - Present approach D - ANSYS D - Present approach

0.8

0.6 0.4 0.2

0.0 0.0

1.0

2.0 Normalized time

3.0

4.0

6. Numerical implementation for asphalt pavement modeling

Fig. 11. Normalized stress sxx =s0 at points C and D computed from ANSYS and the present approach.

A typical pavement structure consists of multiple layers, such as the surface layer, the base/subbase layers, and the

ARTICLE IN PRESS J. Wang, B. Birgisson / Engineering Analysis with Boundary Elements 31 (2007) 226–240

Normalized vertical displacement

0.7

235

Table 1 Material properties for asphalt concrete at 10 1C

0.6

Resilient modulus (ksi)

0.5

Tensile strength (psia)

C - Present approach

0.4

D - Present approach

1149

0.3

a

0.2

234

Creep compliance: JðtÞ ¼ l0 þ l1 tm

l0 ðpsi1 Þ

l1 ðpsi1 Þ

m

3:325  106

3:801  107

0.796

1 psi ¼ 6:89476 kPa.

Table 2 Structural information of the asphalt pavement shown in Fig. 15

0.1 Layers

0.0 0.0

1.0

2.0 Normalized time

3.0

4.0

Fig. 14. Normalized displacement ðuy =aÞ  104 at points C and D computed from the present approach.

AC Base Subgrade a

Elastic modulus (ksi)

E 2 ¼ 50 E 3 ¼ 20

Poisson’s ratio

Thickness (ina)

n1 ¼ 0:41 n2 ¼ 0:35 n3 ¼ 0:40

h1 ¼ 8 h2 ¼ 12 h3 ¼ 200

Width (in)

2w þ wt ¼ 390

1 in ¼ 25:4 mm.

y wt

w

w

pt AC Base

E2, v2

Subgrade

E3, v3

x T B

h1 h2 Point T: (0, 0) Point B: (0, -h1)

h3

Fig. 15. A three-layer pavement structure.

subgrade layer. For asphalt pavements, the surface layer is usually comprised of hot mix asphalt (HMA) mixture (or asphalt concrete (AC)), which is a viscoelastic material exhibiting time-dependent loading response. The rheological behavior of AC is usually represented by the Burgers model or the power-law model as described in Section 2.2 because these models can capture the instantaneous elastic response, viscoelastic response, as well as the creep behavior. The rest layers below the AC layer are usually comprised of elastic materials whose viscosity can be neglected. Due to the rheological behavior of the asphalt mixture, the AC layer will be subjected to strain relaxation if the pavement is continuously loaded, and unloading may cause significant stress redistribution. The numerical method presented in this paper is implemented to study the load-induced responses of asphalt pavements. As depicted in Fig. 15, a typical three-layer asphalt pavement is considered in this study. The viscoelastic behavior of the AC is assumed to be governed by the

power-law model, with the properties at 10 1C1 summarized in Table 1. The structural information of the pavement including the elastic properties of the base and subgrade layers is given in Table 2. The interfaces of two adjacent layers are assumed to be perfectly bonded. A twodimensional BEM model under plane strain condition is applied (Fig. 15). In this model, the pavement structure is discretized into 100 quadratic elements, including 48 elements for the asphalt layer, 36 for the base layer, and 16 for the subgrade layer. A number of smaller-sized elements are placed in the top layer especially around the tire location to cover the critical areas where there is bending action. It is worth mentioning that the mesh shown in Fig. 15 does not represent the actual mesh but an illustration of the nonuniform mesh, which is employed merely for reducing the computation time. The tire load is modeled as a strip (with width wt ¼ 10 in) of uniformly distributed normal pressure, which is applied to the central elements on the surface boundary. To simulate the pavement response under different loading conditions, three different types of loading shown in Fig. 16 are considered:

(a) long-time static load: pt ðtÞ ¼ p0  HðtÞ; (b) load at t ¼ 0 and unload at t ¼ T d : pt ðtÞ ¼ p0  ½HðtÞ  Hðt  T d Þ; P (c) continual repeated loads: pt ðtÞ ¼ N k¼0 p0  ½Hðt kDtÞ  Hðt  kDt  T d Þ, where DT ¼ T d þ T r is the length of one loading cycle, and T d and T r are the loading and resting periods, respectively. 1 The following analyses were conducted using the mixture properties at 10 1C, which is believed to be representative of the intermediate in-service temperatures associated with conditions that are conducive to top-down cracking [4].

ARTICLE IN PRESS J. Wang, B. Birgisson / Engineering Analysis with Boundary Elements 31 (2007) 226–240

236

40 po

po

ΔT

ΔT

t

Tu

(b)

Horizontal stress (psi)

pt(t)

pt(t)

t

(a)

Point T Point B

20

ΔT

0 -20 -40

pt(t)

po Td

-60

Tr

-80 0

200

t

600

800

1000

Fig. 17. Evolution of the horizontal stress (sxx ) at points T and B in the AC layer with time.

Fig. 16. Three different loading configurations.

For all different loading conditions, the load magnitude p0 is set equal to the tire pressure, which is taken as 100 psi in this study. Although the loading types shown in Fig. 16 are simplified models, they are quite representative for simulation of the real rheological behavior of pavement structures. As the first case, consider the pavement subjected to long-time static loading as shown in Fig. 16(a). The evolution of the horizontal stress sxx in the AC layer, also called the bending stress (with tension defined as positive), is investigated in detail. The results for the bending stress presented below have been adjusted by multiplication of a factor so as to give reasonable approximations for the real three-dimensional bending response. This factor, called the bending stress ratio (BSR), can be estimated using a semiempirical equation proposed by Myers et al. [33]: logðBSRÞ ¼ 0:29655ðE 1 =E 2 Þ0:29531 ðh1 =h2 Þ0:95659 ,

400

Time (second)

(41)

where E 2 , h1 , and h2 are defined in Table 2 and E 1 is taken as the resilient modulus of the AC given in Table 1. Fig. 17 shows the evolution of the horizontal stress sxx at two critical points Tð0; 0Þ and Bð0; h1 Þ at the top and bottom of the AC layer (see Fig. 15). It can be seen clearly from Fig. 17 that the stress state is continuously changing as the pavement is loaded. Due to the effect of strain relaxation, both the tensile and compressive stresses at the bottom and top of the AC layer decrease with increasing loading time. In particular, the tension at the bottom of the AC layer turns into compression at longer times. That is to say, the bending stress is diminished eventually at the bottom of the AC layer. Figs. 18 and 19 plot the stress profiles at the surface and bottom of the AC layer at different times ðt ¼ 0; 20; 100; 1000 sÞ, respectively. These figures show how the stresses on the top and bottom boundaries of the AC layer evolve with time. The stress redistribution due to the rheology of the asphalt mixture is clearly demonstrated through Figs. 18 and 19. The bending

20 Horizontal stress at AC surface (psi)

(c)

0

t=0 t = 20 t = 100 t = 1000

-20

-40

-60

-80 0

20

40

40

80

100

Horizontal distance from the center of the load (in) Fig. 18. Evolution of the horizontal stress ðsxx Þ profile at the surface of AC layer.

effect is usually noticeable only in the vicinity of the load area. As the loading time increases, the magnitude of the bending stresses decreases and the influence area becomes smaller. For the case of continuous static loading, bending stresses approach to zero eventually, and only stresses within the zone of influence of the load remain, i.e. only compressive stresses are left at both the surface and bottom of the AC layer. Outside the zone of influence, the stresses are almost negligible. In the second case, a step load depicted in Fig. 16(b) is applied to study the loading and unloading behavior of the asphalt pavement. The pavement is loaded for 100 s (T u ¼ 100 s) and then unloaded. Fig. 20 shows the evolution of the horizontal stresses sxx at the two critical points T and B during the loading and unloading processes. Bending stress reversals are observed at the top and the bottom of the AC layer upon unloading. That is,

ARTICLE IN PRESS J. Wang, B. Birgisson / Engineering Analysis with Boundary Elements 31 (2007) 226–240

Horizontal stress at AC bottom (psi)

50 40 t=0 t = 20 t = 100 t = 1000

30 20 10 0 -10 -20 0

20

40

60

80

100

Horizontal distance from the center of the load (in)

237

these figures that the pavement response is different every time when the pavement is reloaded. This is because the induced bending stresses (tension at the top and compression at the bottom of the AC layer) upon unloading will not be completely gone after a short resting period. As a result, residual tensile and compressive stresses will accumulate at the top and the bottom of the AC layer, respectively. Correspondingly, the dissipated creep strain energy (DCSE) will also accumulate at both the surface and the bottom of the AC layer over each loading cycle. However, the accumulated DCSE per cycle increases over time at the surface while decreases at the bottom. Depending upon the DCSE limit of the asphalt mixture and the rate of DCSE accumulation at the surface and the bottom of the AC layer, the pavement may experience bottom-up or top-down cracking [4].

Fig. 19. Evolution of the horizontal stress (sxx ) profile at the bottom of AC layer.

20 Horizontal stress at point T (psi)

40

Horizontal stress (psi)

20 0 -20

Point T Point B

-40

0 -20 -40 -60 -80

-60

0

-80 0

200

400

600

800

20

Time (second) Fig. 20. Evolution of the horizontal stress (sxx ) at points T and B in the AC layer.

40

60

80

100

Time (second)

1000

Fig. 21. Evolution of the horizontal stress ðsxx Þ point T under cyclic loading (T d ¼ 1, T r ¼ 4 s).

the sxx at point T jumps from high compression to high tension, while at point B the sxx jumps to high compression. This might be due to the rebound of the AC and base layers upon unloading, which induces significant stress redistribution and causes tension at the top and compression at the bottom of the AC layer. In addition, the stresses do not go to zero immediately after unloading because it takes time for the asphalt mixture to recover due to the effect of delayed elasticity. If not reloaded, both the induced tension and compression due to unloading will diminish gradually after a long resting period. As the last case, cyclic loading depicted in Fig. 16(c) is applied to simulate the real traffic loads. Suppose that each loading cycle lasts 5 s, consisting of 1 s load period and 4 s rest period (i.e. T d ¼ 1, T r ¼ 4 s). For this case, the results for the horizontal stresses at points T and B are plotted in Figs. 21 and 22 for the first 20 cycles. It is observed from

Horizontal stress at point B (psi)

50 40 30 20 10 0 -10 -20 0

20

40

60

80

100

Time (second) Fig. 22. Evolution of the horizontal stress (sxx ) point B under cyclic loading (T d ¼ 1, T r ¼ 4 s).

ARTICLE IN PRESS 238

J. Wang, B. Birgisson / Engineering Analysis with Boundary Elements 31 (2007) 226–240

7. Concluding remarks This paper presents a numerical technique to model the time-dependent behavior of linear viscoelastic materials. The approach is based on formulation of the DD method in the time domain and employment of an explicit timemarching scheme. In the spatial domain, the boundary unknowns are approximated by quadratic functions in a piecewise manner; in the time domain, the unknown parameters are assumed to vary linearly over each time step. Several viscoelastic models are adopted to characterize the rheological behavior of linear viscoelastic materials. In most cases, the spatial and time integrals can be evaluated analytically, which gives highly accurate results and fast convergence of the numerical scheme. Benchmark examples are given to demonstrate the efficiency, reliability, and versatility of the present approach. The numerical approach is also implemented to model asphalt pavement and to study the rheological behavior of AC and its effect on the pavement response. The analysis results show that the stress distribution in a typical asphalt pavement is continuously changing as the pavement is loaded, because the viscoelastic asphalt materials relax stresses. Generally, both tensile stresses at the top and compressive stresses at the bottom of the AC layer reduce with increased loading time or number of repeated loads. Continuous repeated loads may lead to accumulation of residual tensile stresses at the surface of the pavement and residual compressive stresses at the bottom of the AC layer. The rate of the DCSE accumulation at the surface and bottom of the AC layer will determine whether bottom-up or top-down cracking is more likely to occur. The research findings imply that the rheological behavior of asphalt mixtures results in load-induced stress redistributions that may dominate the failure mode of a pavement structure. Refinements and extensions are planned in several aspects to further improve the computational efficiency of the present approach. The first natural refinement is to employ a direct boundary integral formulation that automatically takes into account the continuity conditions at the interfaces so as to avoid the explicit substructuring formulation. Another extension of the method is to develop a convolution-free formulation so that only the results from the previous time step are required to compute the response at the current time step (see, for example, [25]). Work on these topics is in progress. Appendix A. Basic spatial integrals The detailed expressions for Ik , Jk ðk ¼ 0; 1; 2Þ in (25) and their derivatives are summarized below [14,32]. For simplicity, the following notations are defined: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 r1 ¼ ðx  aÞ þ y ; r2 ¼ ðx þ aÞ2 þ y2 , y y y1 ¼ arctan ; y2 ¼ arctan . ðA:1Þ xa xþa

(1) I0 and its derivatives:   r2 I0 ¼ ln , r1 xþa xa I0;x ¼ 2  2 , r r 2 1 1 1 I0;y ¼ y 2  2 , r2 r1 I0;yy ¼

r22  2y2 r21  2y2  . r42 r41

(2) I1 and its derivatives:   yðy1  y2 Þ x r2 þ ln I1 ¼  2, a a r1   1 r2 xa xþa  2  2 , I1;x ¼ ln a r1 r r  1  2 y1  y2 1 1 y 2þ 2 , I1;y ¼ a r1 r2       x 1 1 xa 2 xþa 2   2 . I1;yy ¼  2 a r21 r22 r21 r22

ðA:2Þ

ðA:3Þ

(3) I2 and its derivatives:   2xyðy1  y2 Þ y2  x2 r1 2x þ ln  , a2 a a2 r2     2yðy1  y2 Þ 2x r2 xþa xa 4 I2;x ¼ þ ln  þ  , 2 2 2 2 a a a r1 r r    2  1 2xðy1  y2 Þ 2y r1 1 1 I2;y ¼ þ 2 ln þy 2 2 , a2 a r2 r2 r1

I2 ¼

I2;yy ¼

    2 r1 1 2x  a 2x þ a ln þ þ a2 a r2 r21 r22  2  2 xþa xa 2 . þ2 2 r2 r21

ðA:4Þ

(4) J0 and its derivatives: J 0 ¼ y1  y2 , J0;y ¼ I0;x ,   xþa xa J0;yy ¼ 2y  . r42 r41

ðA:5Þ

(5) J1 and its derivatives:   xðy1  y2 Þ y r1 þ ln J1 ¼ , a a r2 J1;y ¼ I1;x ,     y 1 1 xa xþa J1;yy ¼  þ 4  2y . a r21 r22 r41 r2

ðA:6Þ

ARTICLE IN PRESS J. Wang, B. Birgisson / Engineering Analysis with Boundary Elements 31 (2007) 226–240

(6) J2 and its derivatives: J2 ¼

  ðx2  y2 Þðy1  y2 Þ 2xy r1 2y þ ln þ , 2 2 a a a r2

J2;yy ¼

2ðy2  y1 Þ 2y 1 1 þ þ a2 a r21 r22   xa xþa  2y  . r41 r42

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi p21  4p2

¼ a0  a1 , 2p2 p q a0 þ 2p2 q1  p1 q2 c1 ¼ 2 2 , p2 p q a0  2p2 q1 þ p1 q2 c2 ¼ 2 2 , p2

a2 ¼

J2;y ¼ I2;x , 

p1 þ



ðA:7Þ

239

ðB:5Þ

where p1 , p2 and q1 , q2 are combinations of material constants as defined in (11). As a result, the functions for K3 , P03 , and P13 can be expressed as k1 q ðc1 a2 ea2 t  c2 a1 ea1 t Þ þ k1 2 dðtÞ, 2p2 a0 p2 k1 P03 ðtÞ ¼ ½c1 ea2 t þ c2 ea1 t  þ k1 R1 HðtÞ, 2p2 a0   k 1 c2 c1 a1 t a2 t ð1 þ a1 tÞe  ð1  a2 tÞe . P13 ðtÞ ¼ 2p2 a0 a1 a2 K3 ðtÞ ¼

Appendix B. Basic time integrals The detailed expressions for Pm| ðtÞ in (34) are summarized below for three different viscoelastic models—Boltzmann, Burgers, and power-law model. Under the assumption of constant Poisson’s ratios, the functions K1 ðtÞ and K2 ðtÞ defined in (20) can be written in terms of Dirac’s delta functions for all the models considered: K1 ðtÞ ¼ k1 dðtÞ;

K2 ðtÞ ¼ k2 dðtÞ,

(B.1)

where k1 and k2 are constants defined in (18). As a result, Pm1 and Pm2 ðm ¼ 0; 1Þ are obtained as P0j ðtÞ ¼ kj HðtÞ;

P1j ðtÞ ¼ 0 ðj ¼ 1; 2Þ,

(B.2)

where HðtÞ is heaviside step function defined as HðtÞ ¼ 1 for t40 and HðtÞ ¼ 0 for tp0. The expressions of K3 and Pm3 ðm ¼ 0; 1Þ are given below for different models: (1) Boltzmann model: The model constants R1 , R2 , and Z (see Fig. 1(a)) are given, and the Poisson’s ratio n is a constant. The expressions for K3 , P03 , and P13 are given below:   R1 c0 t K3 ðtÞ ¼ k1 R1 dðtÞ  e , Z P03 ðtÞ ¼ k1 R1 ½c1 ec0 t þ HðtÞ, P13 ðtÞ ¼ k1 Zc21 ð1 þ c0 tÞec0 t ,

ðB:3Þ

where k1 is given in (18) and the constants c0 and c1 are defined as c0 ¼

R1 þ R2 ; Z

c1 ¼

R1 . R1 þ R2

(B.4)

(2) Burgers model: For the Burgers model, the following parameters are defined for notational convenience: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi p21  4p2 a0 ¼ , p2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi p1 þ p21  4p2 a1 ¼ , 2p2

ðB:6Þ

(3) Power-law model: The creep compliance for powerlaw model is defined in (14) as a power function of time. Correspondingly, the functions K3 , P03 , and P13 can be written as k 1 l1 m , t1m ðl0 þ l1 tm Þ2 k1 P03 ðtÞ ¼ , l0 þ l1 tm K3 ðtÞ ¼

Z P13 ðtÞ ¼ tP03 ðtÞ 

P03 ðtÞ dt   k1 t k1 t 1 1 l1 m ; 1 þ ;  ¼  O t , l0 þ l1 tm l0 m m l0

ðB:7Þ

where the material constants l0 , l1 , and m are the same as in (14); the hypergeometric function O is defined as Z 1 GðbÞ Oða; b; tÞ ¼ ð1  xÞb2 ð1  xtÞa dx, Gðb  1Þ 0 R1 where GðzÞ ¼ 0 tz1 et dt. Due to the complexity of the function Oða; b; tÞ, itR is suggested to obtain P13 ðtÞ numerically by evaluating P03 ðtÞ dt with numerical quadratures, for example, the Simpson’s rule. References [1] Jacobs MM. Crack growth in asphaltic mixes. PhD dissertation, Delft Institute of Technology, The Netherlands, 1995. [2] Myers L, Roque R, Ruth BE. Mechanisms of surface-initiated longitudinal wheel path cracks in high-type bituminous pavements. J Assoc Asphalt Paving Technol 1998;67:401–32. [3] Wang LB, Myers LA, Mohammad LN, Fu YR. Micromechanics study on top-down cracking. Transp Res Rec 2003;1853:121–33. [4] Roque R, Birgission B, Drakos C, Dietrich B. Development and field evaluation of energy-based criteria for top-down cracking performance of hot mix asphalt. J Assoc Asphalt Paving Technol 2005;74. [5] Birgisson B, Sangpetngam B, Roque R. Prediction of the viscoelastic response and crack growth in asphalt mixtures using the boundary element method. Transp Res Rec 2002;1789:129–35.

ARTICLE IN PRESS 240

J. Wang, B. Birgisson / Engineering Analysis with Boundary Elements 31 (2007) 226–240

[6] Birgisson B, Soranakom C, Napier JAL, Roque R. Microstructure and fracture in asphalt mixtures using a boundary element approach. J Mater Civ Eng 2004;16:116–21. [7] Crouch SL. Solution of plane elasticity problems by the displacement discontinuity method. Int J Numer Methods Eng 1976;10:301–43. [8] Crouch SL, Starfield AM. Boundary element methods in solid mechanics. London: George Allen & Unwin; 1983. [9] Mack MG. The displacement discontinuity method. In: Manolis GD, Davies TG, editors. Boundary element techniques in geomechanics. Southampton: Computational Mechanics Publications; 1993. p. 63–99. [10] Bobet A, Einstein HH. Numerical modeling of fracture coalescence in a model rock material. Int J Fract 1999;92:221–52. [11] Wang Y-J, Crouch SL. Boundary element methods for viscoelastic media. In: Proceedings of the 23rd U.S. Symposium on Rock Mechanics, Berkeley, CA, 1982. p. 704–11. [12] Siebrits E, Crouch SL. Two-dimensional elastodynamic displacement discontinuity method. Int J Numer Methods Eng 1994;37:3229–50. [13] Siebrits E, Birgisson B, Peirce AP, Crouch SL. On the numerical stability of time domain boundary element methods. Int J Blasting Fragmentation 1997;1:305–16. [14] Shou K-J, Crouch SL. A higher order displacement discontinuity method for analysis of crack problems. Int J Rock Mech Min Sci 1995;32:49–55. [15] Napier JAL. Modeling of fracturing near deep level gold mine excavations using a displacement discontinuity approach. In: Rossmanith HP, editor. Mechanics of jointed and faulted rock. Rotterdam: Balkema; 1990. p. 709–15. [16] Napier JAL, Hildyard MW. Simulation of fracture growth around openings in highly stressed, brittle rock. J S Afr Inst Min Metall 1991;91:145–57. [17] Siebrits E, Gu H, Desroches J. An improved pseudo-3D hydraulic fracturing simulator for multiple layered materials. In: Computer methods and advances in geomechanics. Rotterdam: Balkema; 2001. p. 1341–5. [18] Sangpetngam B, Birgisson B, Roque R. Development of an efficient crack growth simulator based on hot mix asphalt fracture mechanics. Transp Res Rec 2003;1832:105–12. [19] Sangpetngam B, Development and evaluation of a viscoelastic boundary element method to predict asphalt pavement cracking. PhD thesis, University of Florida, 2003.

[20] Xu Q, Rahman MS, Tayebali AA. Finite element analysis of layered viscoelastic system under vertical circular loading. In: Pande G, Pietruszczak S, editors. Numerical methods in geomechanics. London: Taylor & Francis; 2004. [21] Wapner PG, Forsman WC. Fourier transform method in linear viscoelastic analysis: the vibrating viscoelastic reed. Trans Soc Rheol 1971;15(4):603–26. [22] Rizzo FJ, Shippy DJ. An application of the correspondence principle of linear viscoelasticity theory. SIAM J Appl Math 1971;21: 321–30. [23] Carini A, De Donato O. Fundamental solutions for linear viscoelastic continua. Int J Solids Struct 1992;29:2989–3009. [24] Sim WJ, Kwak BM. Linear viscoelastic analysis in time domain by boundary element method. Comput Struct 1988;29:531–9. [25] Huang Y, Crouch SL, Mogilevskaya SG. Direct boundary integral procedure for a Boltzmann viscoelastic plane with circular holes and elastic inclusions. Comput Mech 2005;37. [26] Findley WN, Lai JS, Onaran K. Creep and relaxation of nonlinear viscoelastic materials: with an introduction to linear viscoelasticity. New York: Dover; 1989. [27] Tschoegl NW. The phenomenological theory of linear viscoelastic behavior: an introduction. New York: Springer; 1989. [28] Roque R, Buttlar WG. The development of a measurement and analysis system to accurately determine asphalt concrete properties using the indirect tensile mode. J Assoc Asphalt Paving Technol 1992;61:304–32. [29] Roque R, Zhang Z, Sankar B. Determination of crack growth rate parameters of asphalt mixtures using the superpave indirect tension test (IDT). J Assoc Asphalt Paving Technol 1999;68:404–33. [30] Steketee JA. Some geophysical applications of the elasticity theory of dislocations. Cana J Phys 1958;36:1168–98. [31] Crawford AM, Curran JH. Higher-order functional variation displacement discontinuity elements. Int J Rock Mecha Min Sci Geomech Abstr 1982;19:143–8. [32] Wang J, Birgisson B, Roque R, Sangpetngam B. A viscoelastic displacement discontinuity method for analysis of pavements with cracks, Road Mater Pavement Des 2005, in press. [33] Myers L, Roque R, Birgisson B. Use of two-dimensional finite element analysis to represent bending response of asphalt pavement structures. Int J Pavement Eng 2001;2:201–14.