Skip to contentUnited States Department of Transportation - Federal Highway Administration FHWA Home
Research Home   |   Pavements Home
Report
This report is an archived publication and may contain dated technical, contact, and link information
Publication Number: FHWA-HRT-08-073
Date: September 2009

Development of A Multiaxial Viscoelastoplastic Continuum Damage Model for Asphalt Mixtures

Previous | Table of Contents | Next

Chapter 2. Constitutive Modeling with a Multiaxial Viscoelastoplastic Continuum Damage Model

2.1. Historical Perspective

The foundation for the work presented in this report was laid by researchers at Texas A&M University almost two decades ago.(1) These researchers successfully applied Schapery's nonlinear viscoelastic constitutive theory for materials with distributed damage to describe the behavior of sand asphalt under controlled strain cyclic loading. Later research shows that this theory can also describe the behavior of asphalt concrete under both controlled stress and controlled strain cyclic loading.(2-4) Later research shows that the damage characteristics of a material are independent of the mode of loading and can be determined using simpler tests.(5) Further, verification of the t-TS principle at high levels of damage is an equally significant discovery.(6) These two findings significantly reduce the required testing protocol while simultaneously extending the realm of application for the model. Work by Chehab et al. utilizes Schapery's theory, along with strain decomposition, to account for both viscoelastic and viscoplastic strains. (See references 7 through 11.) The most recent work applies this theory to mixtures tested at the FHWA ALF in Mclean, VA,(12) and successfully demonstrates the application of the modeling principles to both modified and unmodified asphalt concrete mixtures. This work is the basis for the current research.

2.2. Modeling Approach

The MVEPCD model is a theoretical and phenomenological extension of the uniaxial VEPCD model presented in previous work and briefly reviewed here. These modeling approaches share many of the same principles, such as linear viscoelasticity, elastic-viscoelastic correspondence principle, continuum damage, and strain-hardening plasticity. A complete review of these principles in the uniaxial sense is given elsewhere, and a brief review is given in the following sections.(13) For a more rigorous treatment of the subject, the reader is referred to previous work and to the work of Schapery for linear viscoelasticity and continuum damage mechanics. (See references 1 through 10.) For a detailed review of strain-hardening viscoplasticity, the reader is directed to the work of others.(7,10,14)

Linear viscoelastic materials exhibit time- and temperature-dependent properties that make them dependent upon the history of loading, unlike elastic materials. This issue complicates the analysis since many continuum theories have been developed assuming an elastic material. In this report, the elastic-viscoelastic correspondence principle is used so that these elastic theories can be applied to the viscoelastic asphalt concrete mixtures. The correspondence principle uses pseudo strains, a quantity calculated from the actual time-dependent strains, in place of the actual strains. It is a general theory in so much that for purely elastic materials, the pseudo strains may be equal to the physical strains. Many damage theories exist, and in this report Schapery's theory based on energy principles is used. For this theory, consideration is given for both the time-dependent nature of energy availability as well as the time dependence of damage resistance. These arguments are considered within a thermodynamic framework. Finally, strain-hardening plasticity is found to be simple to model and commonly observed viscoplastic phenomenon.

A well-known assumption, founded in the theory of plasticity, is that total strain can be decomposed into elastic and plastic strains. Likewise, total strain can be decomposed into elastic strain, plastic strain, and creep strain, according to the theory of viscoplasticity, to account for the rate-dependent plastic strain of materials, as shown equation 1. In some studies that involve such rate-dependent materials, the rate-independent plastic strain and the rate-dependent creep strain are defined as viscoplastic strain because it is difficult to distinguish the plastic deformation from creep deformation. The theoretical background for the viscoplastic strain concept was first discussed by Perzyna, and since then more complicated models have been developed to explain the behavior of a material due to plasticity-creep interaction.(15)

Equation 1.  Strain decomposition equation. The total strain, epsilon superscript total and subscript ij, is given by the elastic strain, epsilon superscript e and subscript ij, plus the plastic strain, epsilon superscript p and subscript ij, plus the creep strain, epsilon superscript c and subscript ij, which also equals elastic strain, epsilon superscript e and subscript ij, plus viscoplastic strain, epsilon superscript vp and subscript ij.       (1)

Where:

Elastic strain. epsilon superscript total and subscript ij. = Elastic strain.

Plastic strain, epsilon superscript e and subscript ij. = Plastic strain.

Creep strain, epsilon superscript c and subscript ij = Creep strain.

Viscoplastic strain, epsilon superscript vp and subscript ij. = Viscoplastic strain.

From a similar perspective, Schapery suggests that total strain may be separated into several components, such as elastic, viscoelastic, and viscoplastic strains.(10) For the MVEPCD model, elastic, linear viscoelastic, and strains due to microcracking damage are combined into a single term, and plastic and viscoplastic strains are combined into another, as shown in equation 2.

Equation 2. Strain decomposition equation modified. The total strain, epsilon superscript total and subscript ij, is given by the viscoelastic strain, epsilon superscript ve and subscript ij, plus the viscoplastic strain, epsilon superscript vp and subscript ij.       (2)

Where:

Viscoelastic strain, epsilon superscript total and subscript ij. = Elastic plus linear viscoelastic strain due to damage.

Viscoplastic strain, epsilon superscript vp and subscript ij. = Viscoplastic strain.

2.2.1. Linear Viscoelasticity

Linear viscoelastic (LVE) materials exhibit time- and temperature-dependent behavior-the current response is dependent on both the current input and all past input (i.e., the input history). By contrast, the response of an elastic material is only dependent on the current input. Constitutive relationships for LVE materials are typically expressed in the convolution integral form, as follows:

Equation 3. Convolution integral for stress. The stress, sigma, is equal to the convolution integral of the relaxation modulus, E, as a function of time, t, minus the time when loading began, tau, and the derivative of the strain, epsilon, with time, tau, with integration between 0 and the time of interest.       (3)

Equation 4. Convolution integral for strain. The strain, epsilon, is equal to the convolution integral of the creep compliance, D, as a function of time, t, minus the time when loading began, tau, and the derivative of the stress, sigma, with time, tau, with integration between 0 and the time of interest.       (4)

Where:

E(t) and D(t) are the relaxation modulus and creep compliance, respectively. In this report, the analytical forms for these functions are given by the common Prony series formulation shown in equation 5 and equation 6. Note that the characterization of LVE behavior is undertaken by performing temperature and frequency sweep tests which are then processed to obtain the coefficients for these functions.(13)

Equation 5. Prony representation of relaxation modulus. The relaxation modulus, E, as a function of time, t, is equal to the long time modulus, E subscript infinity, plus the summation from i equals 1 to the total number of relaxation times, of individual Prony stiffnesses, E subscript i, times the exponential of negative time divided by the relaxation time, rho subscript i, summed for every Prony term.       (5)

Equation 6. Prony representation of creep compliance. The creep compliance, D, as a function of time, t, is equal to the short time or glassy compliance, D subscript g, plus the summation from j equals 1 to the total number of retardation times, of individual Prony compliances, D subscript j, times 1 minus the exponential of negative time divided by the retardation time, tau subscript j, summed for every Prony term.       (6)

2.2.1.1. Linear Viscoelastic Interconversion

The unit response functions presented in equation 3 and equation 4 are often measured in the frequency domain using the complex modulus test because it is often difficult to obtain measurements in the time domain due to limitations of the testing machine's capacity for testing time. The complex modulus provides the constitutive relationship between the stress and strain amplitudes of a material loaded in a steady-state sinusoidal manner. Then, the storage modulus can be determined from the complex modulus and converted to a time-dependent property, such as E(t) and D(t), according to the theory of linear viscoelasticity. When the storage modulus is expressed in terms of reduced angular frequency,ωR, as shown in equation 7, it can be expressed using the Prony series representation given in equation 8.(16,17)

Equation 7. Calculation of storage modulus. The storage modulus, E superscript apostrophe, as a function of reduced angular frequency, lowercase omega subscript R, is equal to the dynamic modulus, in absolute value signs E superscript star, as a function of reduced angular frequency, lowercase omega subscript R, multiplied by the sine of the phase angle, phi, as a function of reduced angular frequency, lowercase omega subscript R.       (7)

Equation 8. Prony representation of storage modulus. The storage modulus, E superscript apostrophe, as a function of reduced angular frequency, lowercase omega subscript R, is equal to the long time modulus, E subscript infinity, plus the summation of individual Prony stiffness, E subscript i, times the individual relaxation time, rho subscript i, times the reduced angular frequency of interest, lowercase omega subscript R, squared, divided by the following: reduced angular frequency, lowercase omega subscript i, squared times the individual relaxation time, rho subscript i, squared, plus 1.       (8)

Where:

Ε∞ = Elastic modulus.

ωR = Angular frequency (=2πfRΔt).

Δt = Time lag between stress and strain.

Ei = Modulus of ith Maxwell element (fitting coefficient).

ρi = Relaxation time (fitting coefficient).

The coefficients determined from this process are then used with equation 5 to find the relaxation modulus. Using the theory of viscoelasticity, the exact relationship between the creep compliance and relaxation modulus is given by equation 9.

Equation 9. Relationship between creep compliance and relaxation modulus. The convolution integral of relaxation modulus, E, as a function of time, t, minus time when load increment began, tau, and the derivative of creep compliance, D as a function of tau, with time, tau, for all time is equal to 1.      (9)

If the creep compliance, written in terms of the Prony representation (equation 10), is substituted into equation 9 along with equation 5 and simplified, then the result can be expressed as a linear algebraic system, shown as equation 11. The coefficients {D} are solved by any proper numerical method.

Equation 10. Prony representation of creep compliance. The creep compliance, D, as a function of time, t, is equal to the short time or glassy compliance, D subscript g, plus the summation of individual Prony compliances, D subscript j, times 1 minus the exponential of negative time divided by the retardation time, tau subscript j, summed for every Prony term.       (10)

Equation 11. Linear equation to solve for interconversion. The matrix A, multiplied by the vector D, is equal to the matrix B.       (11)

Where:

Equation 12. Definition of matrix A. The matrix A is equal to the summation from j equals 1 to the total number of retardation times, M, of the summation from m equals 1 to the total number of relaxation times, N, of the individual Prony stiffness, E subscript m, times the individual relaxation time, rho subscript m, divided by the difference between the individual relaxation time, rho subscript m, and the individual retardation time, tau subscript j, multiplied by the difference between exponential to the negative time, t, divided by the individual relaxation time, rho subscript m, and exponential to the negative time, t, divided by the individual retardation time, tau subscript j, and the long time modulus, E subscript infinity, times 1 minus the exponential of negative time divided by the individual retardation time, tau subscript j.       (12)

Equation 13. Definition of vector D. The vector D is equal to D subscript j.       (13)

Equation 14. Definition of matrix B. The matrix B is equal to 1 minus 1 divided by the long time equilibrium modulus, E subscript infinity, plus the summation from m equals 1 to the total number of individual relaxation terms, of the individual Prony stiffness values, E subscript m, multiplied by the long time modulus, E subscript infinity, plus the summation of individual Prony stiffness, E subscript m, times the exponential of negative time divided by the relaxation time, rho subscript m, summed for every Prony term.       (14)

Once the coefficients Dj are determined, they are substituted into equation 10 to find the creep compliance.

2.2.1.2. Linear Viscoelastic Characterization Methodology

Numerous research efforts by the authors have led to a methodology for assessing and analyzing LVE properties through the dynamic modulus. Explicit details of the experimental method are given in subsection 3.4.1; however, it should be known that the test applies cyclic sinusoidal loading at several combinations of frequency and temperature either with or without confining pressure. Load and axial deformations measured at four locations separated by 90-degree intervals are recorded for each combination of frequency and temperature. From these measurements, stresses and strains are calculated based on the specimen area and gauge length of the deformation measurements, respectively. The analysis procedure is the same for both confined and unconfined tests and is given in detail elsewhere.(13,18)

Asphalt concrete in the LVE range is known to be thermorheologically simple (TRS) and, as such, the effects of time and temperature can be combined into a joint parameter, reduced time/frequency, through the time-temperature shift factor (aT) using equation 15.

Equation 15. Definition of reduced frequency. The reduced frequency, f subscript R, is equal to the frequency, f, multiplied by the time-temperature shift factor, a subscript T.       (15)

In practical terms, this behavior allows for the horizontal shifting of the processed data to form a single curve, the mastercurve, for describing the constitutive behavior of asphalt concrete. The amount of horizontal shift is known as the t-TS shift factor, and the relationship between this factor and temperature is known as the t-TS shift factor function. When combined with the mastercurve, these two functions allow for the prediction of the LVE behavior over a wide range of conditions. The process is shown schematically in figure 1 through figure 3.

2.2.2. Correspondence Principle

The damage theory used in this research, which was originally developed for elastic materials, is generalized for viscoelastic materials using the elastic-viscoelastic correspondence principle.(8) In short, this principle states that viscoelastic problems can be solved with elastic solutions when physical strains are replaced by pseudo strains.

Equation 16. Definition of pseudo strain. The pseudo strain, epsilon superscript R, is equal to the inverse of the reference modulus, E subscript R, multiplied by the convolution integral of the relaxation modulus, E, as a function of time, t, minus the time when loading began, tau, and the derivative of the strain, epsilon, with time, tau, with integration between 0 and the time of interest.       (16)

Where:

ER is a particular reference modulus, typically taken as one. Using pseudo strain in place of physical strain, the uniaxial constitutive relationship presented in equation 3 can be rewritten as follows:

Equation 17. Constitutive equation for linear viscoelastic materials through the correspondence principle. The stress, sigma, is equal to the reference modulus, E subscript R, multiplied by the pseudo strain, epsilon superscript R.       (17)

As equation 17 shows, for the uniaxial condition, a form corresponding to that of a linear elastic material (Hooke's Law) is taken when strains are replaced by pseudo strains. Additional theoretical details of this concept can be found elsewhere.(8,19,20) In a practical sense, pseudo strains are simply the LVE stress response to a particular strain input. The most important effect of pseudo strains is seen when plotting with stress, as the time effects are removed from the resulting graph. This property allows the direct quantification of damage independent of any time effects. The results of two typical monotonic tests are presented in figure 4 and figure 5 in both stress strain space and stress-pseudo strain space. The behavior during initial loading is shown as an inset in these figures. In stress-strain space, as seen in figure 4, nonlinearity appears in the initial stage of loading, which might suggest that damage commences from the outset. However, the nonlinearity in this zone is related only to the time effects of the material. When these time effects are removed, as seen in figure 5, it is clear that damage does not commence at the outset of loading and does not begin until the stress level reaches approximately 500 kPa.

Figure 1. Graph. Schematic representation of dynamic modulus shifting process with unshifted data. This figure the shows dynamic modulus, in absolute value signs E superscript star, in megapascals on the y axis from parenthesis 10 to 1 million close parenthesis in logarithm space, and frequency in Hertz on the x axis from parenthesis 1 times 10 superscript -8 to 1 times 10 superscript 4 close parenthesis in logarithm space, with measured dynamic modulus values at -10 degrees Celsius, 5 degrees Celsius, 20 degrees Celsius, 40 degrees Celsius, and 54 degrees Celsius and a continuous function called a master curve, which follows a decreasing pattern with decreasing frequency and changes in an S-shaped pattern.

Figure 1. Graph. Schematic representation of dynamic modulus shifting process with unshifted data.

Figure 2. Graph. Schematic representation of dynamic modulus shifting process with shifted data. This figure shows dynamic modulus, in absolute value signs E superscript star, in megapascals on the y axis from parenthesis 10 to 1 million close parenthesis in logarithm space, and frequency in Hertz on the x axis from parenthesis 1 times 10 superscript -8 to 1 times 10 superscript 4 close parenthesis in logarithm space. The time-temperature shift factors have been applied to the measured data to form a continuous curve, which lies along the master curve drawn in part.

Figure 2. Graph. Schematic representation of dynamic modulus shifting process with shifted data.

Figure 3. Graph. Schematic representation of dynamic modulus shifting process with time-temperature shift factor. This figure plots the logarithm of the time-temperature shift factor, a subscript T, on the y axis from 3 to -6, and temperature in degrees Celsius on the x axis from -15 degrees Celsius to 65 degrees Celsius. The time-temperature shift factor function is plotted and is shown to decrease with increasing temperature.

Figure 3. Graph. Schematic representation of dynamic modulus shifting process with time-temperature shift factor.

Figure 4. Graph. Constant crosshead test results in stress-strain space. This figure shows the stress in kPa on the y axis from parenthesis 0 to 4,000 close parenthesis, and strain on the x axis from parenthesis 0 to 7 times 10 superscript -3 close parenthesis. Results from tests at two different rates are shown and the slower rate shows less strength. An inset figure shows the behavior between stress of parenthesis 0 and 1,000 close parenthesis kPa with a line showing the initial response, these two data curves deviate from this line at approximately 250 kPa.

Figure 4. Graph. Constant crosshead test results in stress-strain space.

Figure 5. Graph. Constant crosshead test results in stress-pseudo strain space. This figure shows the stress in kPa on the y axis from parenthesis 0 to 4,000 close parenthesis, and the pseudo strain on the x axis from parenthesis 0 to 250,000 close parenthesis. Results from the same two tests shown in figure 4 are also shown. An inset in this figure shows the behavior up to 1,000 kPa along with a line representing the initial stiffness in this space. The two tests deviate from this line at approximately 500 kPa.

Figure 5. Graph. Constant crosshead test results in stress-pseudo strain space.

2.2.3. Continuum Damage

On the simplest level, continuum damage mechanics considers a damaged body with some stiffness as an undamaged body with reduced stiffness. Continuum damage theories thus attempt to quantify two values: damage and effective stiffness. Further, these theories ignore specific microscale behaviors and instead characterize a material using macroscale observations (i.e., the net effect of microstructural changes on observable properties). The most convenient method to assess the effective stiffness in the macro sense is the instantaneous secant modulus. On the other hand, damage is oftentimes more difficult to quantify and generally relies on macroscale measurements combined with rigorous theoretical considerations. For the model at hand, Schapery's work potential theory, based on thermodynamic principles, is appropriate for the purpose of quantifying damage. Within Schapery's theory, damage is quantified by an internal state variable, S, that accounts for microstructural changes in the material. (See references 8 through 11.)

2.2.4. Viscoelastic Continuum Damage Theory

The uniaxial viscoelastic continuum damage (UVECD) model combines elements from the preceding sections to arrive at the constitutive relationship. From continuum damage, the stiffness reduction is defined by the pseudo secant modulus (pseudo stiffness). This quantity is typically normalized for specimen-to-specimen variability by the factor I and denoted as C.

Equation 18. Definition of pseudo stiffness. The pseudo stiffness, C, is equal to the stress, sigma, divided by the pseudo strain, epsilon superscript R, which is multiplied by the normalization factor, I.       (18)

The relationship between damage, S, and the normalized pseudo secant modulus, C, is known as the damage characteristic relationship and is a material function independent of loading conditions.(19) With these considerations, the nonlinear constitutive relationships used for this research are given by equation 19 for stresses and equation 20 for strains.

Equation 19. Constitutive equation for viscoelastic materials with damage through the correspondence principle. The stress, sigma, is equal to the pseudo stiffness, C, which is a function of damage, S, multiplied by the pseudo strain, epsilon superscript R.       (19)

Equation 20. Constitutive equation for viscoelastic materials with damage for the strain response. The strains relating to viscoelastic damage, epsilon subscript ve, is equal to the reference modulus multiplied by the convolution, integral, of the creep compliance, D, as a function of time, t, minus the time when loading began, tau, and the derivative of the stress, sigma, divided by the pseudo stiffness, C, which is a function of damage, S, with time, tau, with integration between 0 and the time of interest.       (20)

In comparing equation 20 to equation 4, a striking similarity is observed. Equation 4 is the constitutive relationship for linear viscoelasticity with a stress input. The modeling approach uses the given input (σ) to determine the input if no damage has occurred (ΕR) and then utilizes the LVE constitutive relationships to find the response. Equation 20 is solved in this research using the same state variable type of approach used for solving equation 16. These formulations are presented elsewhere.(13)

2.2.4.1. Refinement of Damage Characteristic Relationship

The work potential theory specifies an internal state variable, S, to quantify damage. This internal state variable quantifies any microstructural changes that result in the observed stiffness reduction. For asphalt concrete in tension, this variable is related primarily to the microcracking phenomenon. This report highlights only components new to this research, as significant theoretical work has been done by others. (See references 2, 7, 20, and 21.)

The derivation of the UVECD model begins with an assumption of damage behavior (equation 21) or damage evolution law.

Equation 21. Viscoelastic continuum damage model damage evolution law. The rate of damage growth, dS divided by dt, is equal to the negative rate of change of the dual pseudo strain energy density, W subscript d superscript R, with respect to damage, S, raised to the damage evolution rate, alpha.       (21)

The method used to solve the damage evolution law is a matter of preference; therefore, two different solutions are proposed for solving equation 21. The first transforms the original form of the equation into an integral form, assumes α >> 1, and defines a new parameter,S overhat.(21) Equation 22 presents this method in discrete form.

Equation 22. Calculation of damage by method of Park. The damage, S, is equal to modified damage, S overhat, multiplied by 1 plus the reciprocal of damage evolution rate, alpha, all raised to the power of 1 divided by the sum of 1 and the reciprocal of damage evolution rate, alpha.      (22)

Where:

S overhat is given by equation 22, as follows:

Equation 23. Definition of modified damage for method of Park. Modified damage, S overhat, at time step i plus 1 is equal to the modified damage at time step i plus the negative of one half multiplied by change in pseudo stiffness between time step i and time step i plus 1, uppercase delta times C, multiplied by the pseudo strain, epsilon subscript i superscript R, squared, multiplied by time, t, raised to the reciprocal of the damage evolution rate, alpha.       (23)

The second possible solution shown in equation 24 utilizes the chain rule and makes no assumption regardingα.(3) Both methods have been successfully applied in asphalt concrete research.(7,20,21)

Equation 24. Calculation of damage by method of Lee. The damage at the next time step, S subscript i plus 1, equals damage at the current step, S subscript i, plus negative one half multiplied by the change in pseudo stiffness between time step i and time step i plus 1, uppercase delta times C, multiplied by the pseudo strain, epsilon subscript i superscript R, squared, raised to the power of the damage evolution rate, alpha, divided by damage evolution rate plus 1, alpha plus 1, and multiplied by the change in time between step i and i plus 1, uppercase delta multiplied by time, t, raised to 1 divided by the damage evolution rate, alpha, plus 1.        (24)

To reconcile the approximations of these methods, an iterative refinement technique is incorporated into this research. In short, this method assumes that the change in material integrity is sufficiently small over some discrete time step. The rate of change of the material integrity with respect to damage is determined at a point near the current value of damage, Si + S, where the extrapolation error is minimized.

This refinement process begins with an initial calculation of S by either of the approximate methods, both of which require results from constant crosshead rate tests for the stress-pseudo strain relationship. The initial S values are plotted with the pseudo stiffness values, C, to obtain the damage characteristic curve. This relationship is then fitted to some analytical form, as shown in equation 25.

Equation 25. Analytical relationship between damage and pseudo stiffness. The pseudo stiffness, C, equals exponent of coefficient a multiplied by damage, S, raised to the power of coefficient b.       (25)

Returning to equation 21 and noting that the increments of time are generally small, the rate of change in damage can be written as follows:

Equation 26. Discretization of damage evolution rate. The derivate of damage with respect to time can be approximated by taking the finite difference of damage, uppercase delta of S, divided by the finite difference of time, uppercase delta of time.       (26)

By substituting this expression into equation 21 and rearranging and writing in the discrete form, equation 27 is established:

Equation 27. Damage evolution equation in discrete form. The damage and the next time step, S subscript i plus 1, equals the damage at the current time step, S subscript i, plus the finite difference of time, uppercase delta of time, multiplied by the negative rate of change of the dual pseudo strain energy density function, del mark W subscript i superscript R divided by del mark S, raised to the damage evolution rate, alpha.      (27)

For the uniaxial case, the work function, WR, is given by equation 28.(19)

Equation 28. Definition of dual pseudo strain energy density function for uniaxial conditions. The dual pseudo strain energy density function, W subscript d superscript R, is equal to one half multiplied by pseudo stiffness, C, which is a function of damage, S, multiplied by the pseudo strain, epsilon superscript R.       (28)

By substituting equation 28 into equation 27 and then simplifying, equation 29 is reached.

Equation 29. Discrete form of damage growth equation used for numerical solution. The damage at the next time step, S subscript i plus 1, equals damage at the current step, S subscript i, plus finite difference in time, uppercase delta of t, multiplied by negative one half multiplied by the pseudo strain at the current time step, epsilon subscript i superscript R, squared, multiplied by the change in pseudo stiffness with respect to damage at the current step, del mark C subscript i divided by del mark S, raised to the damage evolution rate, alpha.       (29)

In equation 29, it is assumed that before loading occurs, S and C are zero and 1, respectively. Furthermore, S must be specified and should be significantly less than the change in damage over a time step (typically, 0.1 is used). After calculating the value of damage, Si, and the incremental damage, Si +S, at a given time step, the corresponding values of C are found by equation 25. The difference between these values (C) is then used to calculate damage at the next time step. The process is repeated until all data points are processed.

After completing this first iteration, the new values of S are plotted against the original pseudo stiffness values, and a new analytical relationship is found. The entire process is repeated until the change in successive iterations is small. In this research, eight such iterations were performed, but very little improvement was made after the third or fourth iteration.

Figure 6 and figure 7 present the initial S calculated by both approximate techniques along with results from the refinement process. From these figures, it is seen that the refinement process results in S values that fall between the two approximate methods. In these figures, the seed values for the refinement process are obtained by the chain rule method. However, trials show that regardless of the method used to find the seed values, iterations collapse to the same curve. Details regarding this refinement process can be found elsewhere.(22)

The final refined values of S are plotted with C to obtain the true damage characteristic curve. As noted before, the benefit of this curve is its mode-of-loading independence; the curve represents a fundamental behavior of the material. From equation 16 and equation 19, if this fundamental relationship is known and the relaxation modulus is given, then the stresses can be directly calculated. Conversely, if this fundamental relationship is known and the creep compliance is given, the strain response can be directly calculated for any stress input using equation 20.

Figure 6. Graph. Comparison of refined and approximate damage calculation techniques. This figure shows the damage, S, is plotted on the y axis from parenthesis 0 to 100,000 close parenthesis, and the reduced time is plotted on the x axis from parenthesis 0 to 125 close parenthesis. The graph shows the results of damage, S, refinement starting from the method of Park and the Chain rule method of Lee and Kim. Although both of these approaches have different initial values, both converge to a middle answer upon refinement. This figure was provided by the Transportation Research Board, National Academy of Science from Appendixes to NCHRP Report 547: Sample Performance Tests and Advanced Materials Characterization Models (2005).

Transportation Research Board, National Academy of Science. From Appendixes to National Cooperative Highway Research Program
(NCHRP) Report 547: Sample Performance Tests and Advanced Materials Characterization Models (2005)

Figure 6. Graph. Comparison of refined and approximate damage calculation techniques.(22)

Figure 7. Graph. Comparison of refined and approximate damage characteristic relationship. This figure shows the pseudo stiffness, C, is plotted on the y axis from parenthesis 0 to 1 close parenthesis, and the damage, S, is plotted on the x axis from parenthesis 0 to 60,000 close parenthesis. The relationship is shown for both the method of Park and the Chain rule method of Lee and Kim, as well as the refined values. All methods show an exponential decay pattern with the refined values lying in between the two methods. This figure was provided by the Transportation Research Board, National Academy of Science from Appendixes to NCHRP Report 547: Sample Performance Tests and Advanced Materials Characterization Models (2005).

Transportation Research Board, National Academy of Science. From Appendixes to NCHRP Report 547:
Sample Performance Tests and Advanced Materials Characterization Models
(2005)

Figure 7. Graph. Comparison of refined and approximate damage characteristic relationship.(22)

2.2.4.2. Calculation of Damage

In a typical application, the user wants to predict strain from only stress and time using equation 20. However, from the preceding discussion, the user must know the strain to compute
S and C. To bypass this requirement, it is first assumed that S and C are initially zero and 1, respectively. Then, at some incremental amount of damage (S), C can be computed by the functional relationship found through the refinement process. Finally, using the relationship in equation 19, equation 29 can be rewritten as follows:

Equation 30. Discrete form of damage growth equation used for numerical solution when stress is known. The damage at the next time step, S subscript i plus 1, equals damage at the current step, S subscript i, plus finite difference in time, uppercase delta of t, multiplied by negative one half multiplied by the stress at the current time step, sigma subscript i, divided by the pseudo stiffness at the current time step, C subscript i, squared, multiplied by the change in pseudo stiffness with respect to damage at the current step, del mark C subscript i divided by del mark S, raised to the damage evolution rate, alpha.       (30)

Through equation 30, damage may be incrementally calculated using a technique similar to that used to refine S using only stress and time.

2.2.4.3. Viscoelastic Damage Characterization Methodology

Viscoelastic damage characterization refers to the development of the characteristic damage relationship (i.e., the C versus S relationship). Although such characterization can be performed under any loading condition, the simplest method is the constant crosshead rate test.(20) For MVEPCD characterization, both confined and unconfined constant crosshead rate tests are performed. The only requirement for this test is that viscoelastic damage mechanisms dominate the material behavior, such as when the material is at a low temperature or when it is loaded at a very fast rate. In this study, characterization is performed at 5 °C. This temperature is convenient because it allows for moderate strain rates and does not require consideration of any dynamic effects that might be related to extremely high strain rates. In practice, at least three different rates are tested at 5 °C, and if the resulting C versus S curves collapse, it is known that viscoelastic damage mechanisms dominate. The appropriate strain rates are material-dependent.

For both VEPCD and MVEPCD characterization, the first step is the calculation of the pseudo strains by equation 16 for the constant rate tests. These values are used along with equation 18 and equation 24 to calculate the initial damage, S, values. The relationship between C and these initial S values are fitted to equation 25 and refined through the aforementioned procedure. The true C versus S relationship is obtained by refitting equation 25 to these refined S values. Results from the tests at different rates are finally averaged to obtain the C versus S relationship for modeling. Typically, each test is sampled at different intervals, which complicates this averaging process. To overcome this difficulty, a common list of C values is compiled, and the corresponding S values are interpolated from the respective test. Averaging can then be easily performed because each test has common C values.

2.3. Multiaxial Viscoelastic Continuum Damage Model

2.3.1. Review of Mechanistic Principles

The arguments presented in this section focus on general stress-strain relationships. For notational simplicity, the derivations in this section are given for linear elasticity; however, the arguments remain valid for viscoelastic materials through the elastic-viscoelastic correspondence principle. Much information about this principle has been provided elsewhere; the form given by Schapery is adopted here.(8)

For linear elastic materials, the most general constitutive relationship is given by the following:

Equation 31. General material constitutive relationship. The general stress tensor, sigma subscript ij, is equal to the general stiffness matrix, C subscript ijkl, multiplied by the strain, epsilon subscript kl.    (31)

Cijkl is a fourth-order tensor with 81 constants. Through symmetry arguments, this number is reduced to 21. In matrix form, equation 31 can be viewed as equation 32.

Equation 32. General stiffness matrix. The stiffness matrix, bold C, is equal to 6 by 6 matrix; terms in the matrix reflect symmetry and total matrix. Beginning at the upper left corner terms are labeled C subscript 11, then C subscript 12 and so on. Across the diagonal of the matrix is labeled C subscript 11, C subscript 22, C subscript 33, etc.       (32)

If further assumptions such as orthotropic isotropy, transverse isotropy, and cubical isotropy are assumed, these numbers reduce further. The simplest of these cases is full isotropy, which assumes no directional dependence on the material properties. In this case, two material constants are necessary to describe the constitutive relationship between stress and strain. In general terms, these are called Lame's constants and together define the more commonly recognized fundamental parameters: Poisson's ratio (ν), shear modulus (G), Young's modulus (E), and bulk modulus (K). In terms of compliances, these constants define the shear compliance (J), longitudinal compliance (D), and bulk compliance (B).

For a general isotropic case, it can be shown that C11, C22, and C33 are equal; C12, C13, and C23 are equal; and C44, C55, and C66 equal either half of the shear modulus or just the shear modulus depending on the definition used for shear strain. All other terms become zero. Written in full, assuming the shear strain in equation 31 is given as the total angle form, the stiffness matrix for the general isotropic case becomes equation 33.

Equation 33. Stiffness matrix for general isotropic material. The stiffness matrix, strong C, is equal to Young’s modulus, E, divided by the product of 1 plus Poisson’s ratio, lowercase nu, and 1 minus two times Poisson’s ratio, lowercase nu, multiplied by 6 by 6 matrix, where the upper left and first three diagonal terms of this matrix are equal to 1 minus Poisson’s ratio, lowercase nu, the off diagonal elements in the first 3 by 3 block of the matrix are equal to Poisson’s ratio, lowercase nu, the remaining diagonal elements are equal to 1 minus two times Poisson’s ratio, lowercase nu, divided by two, and all other elements are 0.        (33)

It is seen from equation 31 that the use of equation 33 recovers the common form of Hooke's law.(23) For the problem at hand, however, isotropy is of less concern than transverse isotropy. For this case, it can be shown that the stiffness matrix is given by equation 34. A change in nomenclature from Cijto Zij for general stiffness is made to avoid confusion with later terminology. Also, axis three is assumed to be the axis of symmetry.

Equation 34. Stiffness matrix for transversely isotropic material. The stiffness matrix for transversely isotropic material, strong Z, is equal to 6 by 6 matrix, where the matrix is symmetric about the diagonal with the first two diagonal elements equal—this is the isotropic plane. The first two elements in the third column are also equal but different than the value in the second column first row. The element in the third column and third row is the symmetry axis and is different than the other elements. The fourth and fifth diagonal elements are equal but different than the sixth diagonal element, and all other elements are 0.       (34)

For isotropy, it is seen that the stiffness matrix terms are related to the engineering parameters, Young's modulus and Poisson's ratio. Using similar terminology, transverse isotropy can be considered to have two Young's moduli and three Poisson's ratios. The nomenclature used in this report for these parameters is as follows:

Equation 35. Definition of engineering parameters in transverse isotropy. The stiffness on the isotropy plane, E, is equal to the stiffness on plane one, E subscript 1, and the stiffness on plane two, E subscript 2. The stiffness along the axis of symmetry is given by E subscript 3. The Poisson’s ratio on the isotropy plane is equal to lowercase nu subscript 12. The Poisson’s ratio between axis of symmetry and isotropy plane, lowercase nu subscript 3,132, is equal to the negative of the strain on the isotropy plane, epsilon subscript 11, divided by strain along the axis of symmetry, epsilon subscript 33, and is also equal to the negative of the strain on the other isotropy direction, epsilon subscript 22, divided by strain along the axis of symmetry, epsilon subscript 33. The Poisson’s ratio between the isotropy plane and axis of symmetry, lowercase nu subscript 1,323, is equal to the negative of the strain along the axis of symmetry, epsilon subscript 33, divided by the strain on the transverse plane, either epsilon subscript 11 or epsilon subscript 22. The shear modulus between the transverse plan and the axis of symmetry is given by G subscript 13 or G subscript 23. The shear modulus between isotropic planes is given by G subscript 12.         (35)

Where:

ν1323 is not equal to ν3132; however, for symmetry, ν3132/E3 = ν1323/E.

The relationships between Zij and the values defined in equation 35 are as follows:

View alternative text.          (36)

Or, solving for the engineering parameters, the relationships are as follows:

View alternative text        (37)

If the compliance matrix with generalized compliances, Sij, is considered, then the resulting relationships are as shown in equation 38.

View for alternative text.         (38)

The compliance matrix is considered here for completeness. Also, the constitutive relationships are more compact and easier to present.

2.3.2. Application of Schapery's Work Potential Theory

Schapery's work potential theory model contains parameters that characterize the material integrity with damage growth.(22,11) To clarify the issue of the physical significance of these parameters, a step-by-step process is followed to link the damage functions to the engineering parameters defined in equation 37.

In terms of principal strains, the strain energy density function as defined by the terms found in equation 34 is shown in equation 39:

View for alternative text      (39)

Schapery shows that a more convenient way to express the strain energy density function for transversely isotropic materials is as follows:

Equation 40. Strain energy density function from Schapery’s derivation with transverse isotropy. The energy, W, is equal to one half multiplied by the product of the first stiffness term, A subscript 11, and dilation, e subscript v, squared, plus the second stiffness term, A subscript 22, multiplied by the primary loading direction deviatoric strain, e subscript 3, squared, plus 2 times the third stiffness term, A subscript 12, multiplied by the product of the dilation, e subscript v, and the primary loading direction deviatoric strain, e subscript 3, plus the fourth stiffness term, A subscript 66, multiplied by the strain difference term, e subscript 2, squared.        (40)

Where:

Equation 41. Definition of dilation. Dilation, e subscript v, is equal to the sum of the strains along the three directions, epsilon subscript 1 plus epsilon subscript 2 plus epsilon subscript 3.       (41)

The stiffness values, Aij, in equation 40 can be damage dependent for both the elastic and viscoelastic case. Through an expansion of equation 40, the Aij values correspond to the more recognizable forms in equation 39, thusly:

View for alternative text       (42)

Schapery presents the energy density function as a dual energy density function so that the damage characteristics of a cylindrical body subjected to confining pressure, p, and axial deformation can be more easily characterized. This dual energy density function, as derived by Schapery, is presented in equation 43. To arrive at the relationship presented in equation 43, Schapery assumes that work potential exists and uses the principle of virtual work on a cylindrical specimen subjected to axial loading and confining pressure.(11) Then, he defines a general work potential theory and expands it into a power series, canceling terms to satisfy the initial conditions. The relationships between the Aij terms and the damage-dependent variables, Cij, are presented in equation 44. These relationships must be satisfied for consistency between equation 40 and equation 43. The Cij terms in equation 43 and equation 44 are not the same as those defined in equation 32.

Equation 43. Definition of the dual strain energy density function. The dual energy density function, W subscript d, is equal to one half multiplied by the product of first material integrity parameter, C subscript 11, and strain, epsilon, squared, plus the product of the second material integrity parameter, C subscript 12, strain, epsilon, and pressure, p, plus one half multiplied by the product of the third material integrity parameter, C subscript 22, and pressure, p, squared.       (43)

 

Equation 44. Relationship between first Schapery stiffness term and material integrity parameters. The first Schapery stiffness term, A subscript 11, is equal to 1 divided by 9 multiplied by the first material integrity term, C subscript 11, minus the second material integrity parameter, C subscript 12, minus 3 squared and divided by the third material integrity parameter, C subscript 22. The second Schapery stiffness term, A subscript 12, is equal to the first material integrity parameter, C subscript 11, divided by 3 plus the quotient of the second material integrity parameter, C subscript 12, by the third material integrity parameter, C subscript 22, multiplied by 1 minus the second material integrity parameter, C subscript 12, divided by 3. The third Schapery stiffness term, A subscript 22, is equal to the first material integrity parameter, C subscript 11, minus the quotient of the second material integrity parameter, C subscript 12, squared and the third material integrity parameter, C subscript 22. The fourth Schapery stiffness term, A subscript 66, is equal to the initial material shear modulus, G subscript 0.         (44)

Making substitutions and canceling terms, the engineering parameters in equation 37 are expressed as a function of Cij:

View for alternative text.       (45)

For the compliance matrix, the following can be seen:

View for alternative text.        (46)

It is worthwhile to consider the definitions of stress and strain used in equation 43. Stress is defined as the stress above stress due to pressure (i.e., the deviator stress). According to the definition of strain given in the formulation of the work potential theory presented in the Park dissertation, strain includes all changes in length, even those associated with volume change.(24) In an earlier paper, strain is defined clearly as the change in axial displacement due to load divided by the initial length.(11) To clarify the confusion, virtual work is defined by the following:

Equation 47. Definition of virtual work. The change in virtual work, lowercase delta open parenthesis the word “work” and close parenthesis, is equal to the summation from step i to step N of the product of the generalized forces, Q subscript i, and the change in virtual displacements, dq subscript i.       (47)

In equation 47, Qi are the generalized loads, and qi are the generalized displacements. Whether the strain is defined as the total strain, the evaluated value of equation 47 does not change because the derivative of the strain is independent of these two definitions. Still, the issue does present some consequences in the characterization stage, and so, equation 40 and the relationships in equation 44 will be used to clarify the issue. For a transversely isotropic axial specimen subjected to confining pressure and axial elongation along the axis of symmetry, the following can be shown:

Equation 48. Finding stress using virtual work principles. The partial derivative of work with respect to strain along the axis of symmetry, del W by del epsilon subscript 3, equals stress along the axis of symmetry, sigma subscript 3, equals strain along the symmetry axis, epsilon subscript 3, multiplied by Schapery’s first stiffness term, A subscript 11, plus 4 divided by 3 multiplied by Schapery’s second stiffness term, A subscript 12, plus 4 divided by 9 multiplied by Schapery’s third stiffness term, A subscript 22, plus 2 multiplied by strain along the transverse direction, epsilon subscript 1, multiplied by Schapery’s first stiffness term, A subscript 11, plus 1 divided by 3 multiplied by Schapery’s second stiffness term, A subscript 12, minus 2 divided by 9 multiplied by Schapery’s third stiffness term, A subscript 22.       (48)

Using the relationships in equation 44, equation 48 can be recast with Upsilon defined the same way as ev from equation 40 as the following:

Equation 49. Showing stress as a function of material integrity terms. The stress along the axis of symmetry, sigma subscript 3, equals the strain along the symmetry axis, epsilon subscript 3, multiplied by the first material integrity parameter, C subscript 11, plus the quotient of the second material integrity parameter, C subscript 12, minus 1 and the third material integrity parameter, C subscript 22, multiplied by dilation, lowercase upsilon, minus the second material integrity parameter multiplied by the strain along the symmetry axis, epsilon subscript 3.       (49)

From equation 43, C22 and C12 are related by the following:

Equation 50. Relationship between the second and third material integrity terms. The second material integrity term, C subscript 22, is equal to dilation, lowercase upsilon, minus the product of the second material integrity parameter, C subscript 12, and strain along the symmetry axis, divided by pressure, p.       (50)

If it is assumed that the strain in equation 50 is the total strain,ε3, and that equation 50 is substituted into equation 49 and rearranged, then the following occurs:

Equation 51. Finding the definition of stress and strain. The stress along the axis of symmetry, sigma subscript 3, plus pressure, p, is equal to the product of strain along the symmetry axis, epsilon subscript 3, and the first material integrity term, C subscript 11, plus the product of the second material integrity parameter, C subscript 12, and pressure, p.       (51)

Equation 51 is compared to the similar form obtained from equation 43 to achieve equation 52:

Equation 52. Definition of stress through virtual work and material integrity parameters. The partial derivative of the dual energy density function with respect to strain, del W subscript d divided by del epsilon, equals deviatoric stress, sigma, equals the product of the total strain along the axis of symmetry, epsilon, and the first material integrity term, C subscript 11, plus the product of the second material integrity parameter, C subscript 12, and pressure, p.       (52)

Observing that the sign convention for pressure is opposite that of stress above pressure, a comparison of equation 51 and equation 52 indicates that the stress in equation 43 is the deviator stress, and the strain is the total strain.

2.4. Theory of Viscoplasticity

The hot mix asphalt (HMA) mixture is a pavement material that exhibits both viscoelastic and viscoplastic behavior and therefore shows complicated rate-dependent behavior for repetitive traffic loadings. Rutting, one of the major distress types in HMA pavements, is directly related to the rate-dependent permanent deformation behavior of HMA. To predict the rutting performance of HMA, much effort has been made to develop constitutive models capable of describing the rate-dependent permanent strain development in HMA. Some researchers suggest a viscoplastic model with strain-hardening features.(25) As the simplest model, it is also able to describe monotonic behavior in tension, as shown by researchers.(7,10) Others have presented a viscoplastic model for HMA that incorporates Perzyna's flow rule with Desai's yield function.(26,27) While still others suggest a simplified hierarchical single surface (HISS)-Perzyna model that shows a reasonable viscoplastic strain prediction.(28) However, the difficulty in developing a constitutive model for HMA is that hardening rules based on the behavior of metals or soils are not necessarily appropriate to describe both viscoplastic and viscoelastic behaviors of HMA. That is, models developed for metals and soils describe only elastic recovery during unloading, whereas HMA shows nonlinear strain recovery during unloading due to the viscoelastic property of the material. An important observation supporting this phenomenon was made by Saadeh.(29)

2.4.1. Flow Rule

The general concepts behind the constitutive equations for plastic deformation were proposed by Von Mises based on the theory of elasticity.(30) As such, the strain tensor is related to the stress through an elastic potential function, the complementary strain energy, Ue.

Equation 53. Relationship between strain tensor and complementary strain energy. The elastic strain tensor, epsilon subscripts ij and superscript e, is equal to partial derivative of the complementary strain energy, U subscript e, with respect to stress tensor, sigma subscripts ij.       (53)

Where:

εije= Elastic strain tensor.

Ue = Elastic complementary strain energy.

σij = Stress tensor.

The plasticity theory based on the above flow rule is called the plastic potential theory. When the state of stress reaches the yield criterion, f, plastic strain develops; this mechanism is called plastic flow. To generalize this concept to the plasticity theory, Von Mises proposes that a plastic potential function, g(σij), exists.(30) The plastic strain rate, dεijp, can then be derived from the following flow rule:

Equation 54. Definition of flow rule for plastic strain. The incremental plastic strain tensor, epsilon subscripts ij and superscript p, is equal to positive scalar factor, lamda, multiplied by the partial derivative of potential function, g, with respect to stress tensor, sigma subscript ij.      (54)

Where:

dεijp = Plastic strain rate.

λ = Positive scalar factor.

the partial derivative of potential function, g, with respect to stress tensor, sigma subscript ij. = Gradient of the plastic potential, gij).

In equation 54, λ is a proportional positive scalar factor that can be determined from the yield criterion. For some materials, the plastic potential function, g, and the yield function, f, can be assumed to be the same. These kinds of materials are considered to follow the associative flow rule of plasticity. In other words, the normality rule for this material is associated with the yield criterion, f. However, for other materials, the plastic potential function, g, and the yield function, f, are different. These materials are considered to follow the nonassociative flow rule of plasticity, and the flow rule of the material is derived from a plastic potential, g. In this case, ∂g/ ∂σ ij is not proportional to ∂f/∂σij, and therefore, the plastic strain direction is not normal to the yield surface, f. Basically, the viscoplastic flow rule takes a similar form to the plastic flow rule except that it takes advantage of the over-stress function as a replacement forλ,as shown in equation 55.

Equation 55. Definition of flow rule for viscoplastic strain rate. The viscoplastic strain rate, epsilon subscript ij and superscript vp, is equal to lamda multiplied by the partial derivative of yield function, f, with respect to stress tensor, sigma subscript ij, is equal to the overstress function, iota, divided by the viscosity, eta subscript 0, multiplied by the partial derivative of yield function, f, with respect to stress tensor, sigma subscript ij, is also equal to the overstress function, iota, multiplied by the fluidity, gamma, multiplied by the partial derivative of yield function, f, with respect to stress tensor, sigma subscript ij.       (55)

Where:

Φ = Overstress function.

η0 = Viscosity.

Figure 8 represents the viscoplastic flow rule using a mechanical analog, which combines a dashpot and a slip element in parallel.(31) In this case, the overstress function in equation 55 is represented by the difference from the applied stress (σ) and the yield stress (G) to the viscosity (η0) in this model.

Figure 8. Illustration. Mechanical analog for the viscoplastic model. This figure show the mechanical analog consists of dashpot with viscosity labeled as eta subscript 0 and a slip element with friction G. The two components are assembled in parallel.        

Figure 8. Illustration. Mechanical analog for the viscoplastic model.

2.4.2. Yield Criterion

In conventional viscoplasticity, the elastic limit of the material can be defined by a surface-in-stress space. Mathematically, the yield surface for general anisotropic materials is expressed as equation 56.

Equation 56. Definition of yield function in stress space. The yield function, f, from open parenthesis sigma subscript ij, close parenthesis is equal to 0.       (56)

For an isotropic material, the orientation of the principal axis is immaterial, and the principal stresses, σ11, σ22, and σ33, are sufficient to describe the state of stress. The principal stresses form the integrity basis; it is common to use I1, J2, and J3 as the integrity basis. Therefore, the yield function becomes equation 57 for the isotropic material.

Equation 57. Definition of yield function in principal stress space. The yield function, f, parenthesis first principle stress, sigma subscript 11, comma second principle stress, sigma subscript 22, comma the third principle stress, sigma subscript 33, close parenthesis, can be represented by f, parenthesis first stress invariant, I subscript 1, comma second deviatoric stress invariant, J subscript 2, comma third deviatoric stress invariant, J subscript 3 close parenthesis.       (57)

Where:

I1 = σ 11+ σ 2233      (58)

Equation 59. Definition of second deviatoric stress invariant in yield function. The second deviatoric stress invariant, J subscript 2, in index notation is equal to the product of one half, deviatoric stress, s subscript ij, and deviatoric stress, s subscript ij.      (59)

Equation 60. Definition of third deviatoric stress invariant in yield function. The third deviatoric stress invariant, J subscript 3, in index notation is equal to the product of one third, deviatoric stress, s subscript ij, deviatoric stress, s subscript jk, and deviatoric stress, s subscript ki.       (60)

Equation 61. Definition of deviatoric stress. The deviatoric stress, s subscript ij, in index notation is equal to the difference between total stress, sigma subscript ij, and the product of product of one third, hydrostatic stress, sigma subscript kk, and the kroneker delta, delta subscript ij.       (61)

Physically, I1 represents the hydrostatic pressure, and J2 represents the distortional energy in the material; no clear physical meaning is related to J3. Generally, yield criteria can be classified into two subgroups according to their dependence on hydrostatic stress. The isotropic hardening model is the simplest hardening model and is based on the assumption that the yield surface expands isotropically as the plastic strain develops. The typical isotropic hardening model is presented in equation 62 and figure 9. Because the loading surface expands uniformly, it cannot account for the Bauschinger effect observed in various materials, which describes the reduction of compressive yield strength due to a previously applied tensile stress, or vice versa. Therefore, using only the isotropic hardening model frequently limits the characterization of the material behavior when both tension and compression loads are applied.

Equation 62. Concept of isotropic hardening. The yield function, f, parenthesis stress tensor, sigma subscript ij, comma hardening, K, close parenthesis is equal to parenthesis deviatoric stress tensor, s subscript ij, close parenthesis, minus isotropic hardening parameter, kappa.       (62)

Where:

s-ij = Deviatoric stresses.

Κ = Isotropic hardening parameter.

The kinematic hardening model assumes that during plastic deformation the yield surface translates as a rigid body in the stress space and has the same shape and size as the initial yield surface. The kinematic hardening model is represented by equation 63 and figure 10.

Equation 63.Concept of kinematic hardening. The yield function, f, parenthesis stress tensor, sigma subscript ij, comma hardening coefficient a close parenthesis is equal to f, parenthesis deviatoric stress tensor, s subscript ij, minus kinematic hardening parameters, alpha subscript ij, close parenthesis.       (63)

Where:

s-ij = Deviatoric stresses.

αij = Kinematic hardening parameters (i.e., coordinates of the center of the yield surface in the deviatoric stress space).

Equation 64 and equation 65 conceptually represent the classic kinematic hardening rules suggested by Prager, Armstrong-Frederick, and Chaboche, respectively.(32-34) The kinematic hardening rate is a function of the viscoplastic strain rate in these equations, and the hardening rate is always zero when there is no change in viscoplastic strain. Therefore, when a material is subjected to one-directional loading, such as a constant strain-rate test and a repetitive creep and recovery test, the yield stress increases only in the direction of the developed viscoplastic strain and never decreases.

Equation 64. Hardening rule for isotropic hardening. The rate of hardening, alpha overdot, is equal to the plastic potential function, g, parenthesis viscoplastic strain rate, epsilon superscript vp, close parenthesis.       (64)

Equation 65. Hardening rule for kinematic hardening. The rate of hardening, alpha is equal to the plastic potential function, g, parenthesis viscoplastic strain rate, epsilon superscript vp, comma kinematic hardening stress, alpha, close parenthesis.       (65)

Where:

α = Kinematic hardening stress.

 epsilon superscript vp = Viscoplastic strain rate.

Figure 9. Illustration. Isotropic hardening diagram. This figure shows the concept of isotropic hardening. The second principle stress, sigma subscript 2, is plotted on the y axis, and the first principle stress, sigma subscript 1, is plotted on the x axis. An initial yield surface is shown, and a subsequent yield surface is also shown. For isotropic hardening, the subsequent yield surface is shown centered on the initial yield surface but is larger in diameter.

Figure 9. Illustration. Isotropic hardening diagram.

Figure 10. Illustration. Kinematic hardening diagram. This figure shows the concept of kinematic hardening. The second principle stress, sigma subscript 2, is plotted on the y axis, and the first principle stress, sigma subscript 1, is plotted on the x axis. An initial yield surface is shown, and a subsequent yield surface is also shown. For kinematic hardening, the subsequent yield surface is the same size as the initial yield surface but has moved to the right of the initial yield surface.

Figure 10. Illustration. Kinematic hardening diagram.

2.5. Viscoplastic Models

2.5.1. Simple Strain-Hardening Model

A simple strain-hardening model has been suggested and is shown in equation 66, which assumes that viscosity obeys a power law in viscoplasticity. Researchers have shown that the model is applicable to monotonic behavior in tension.(7,10,25)

Equation 66. Simple strain hardening model. The rate of strain growth, epsilon overdot subscript vp, equals a stress function, g, parenthesis then sigma then close parenthesis, divided by viscosity, eta, as a function of viscoplastic strain level, epsilon subscript vp.       (66)

Where:

g(σ) = Stress function.

η = Viscosity.

Assuming that η is a power law in the viscoplastic strain, equation 67 becomes the following:

Equation 67. Simple strain hardening model with viscosity function. The rate of strain growth, epsilon overdot subscript vp, equals a stress function, g, parenthesis then sigma then close parenthesis, divided by coefficient A multiplied by viscoplastic strain level, epsilon subscript vp, raised to the power of coefficient p.    (67)

Where:

A and p are model coefficients. Rearranging and integrating yield the following:

Equation 68. Simple strain hardening model with viscosity function rearranged. The viscoplastic strain level, epsilon subscript vp, raised to the power of coefficient p and multiplied by the change in viscoplastic strain, depsilon subscript vp, equals a stress function, g, parenthesis then sigma then close parenthesis, divided by coefficient A and multiplied by discrete time, dt.       (68)

Equation 69. Simple strain hardening model with viscosity function integrated. The viscoplastic strain level, epsilon subscript vp, raised to the power of coefficient p plus 1 equals the integral of a stress function, g, parenthesis then sigma then close parenthesis, from time equal 0 to the time of interest, t, multiplied by the ratio of coefficient p plus 1 divided by coefficient A.       (69)

Raising both sides of equation 69 to the 1/(p + 1) power yields the following:

Equation 70. Simple strain hardening model with viscosity function integrated and simplified. The viscoplastic strain level, epsilon subscript vp, equals the integral of a stress function, g, parenthesis then sigma then close parenthesis, from time equal 0 to the time of interest, t, multiplied by the ratio of coefficient p plus 1 divided by coefficient A, raised to the power of 1 divided by 1 plus coefficient p.       (70)

Letting g(σ) = Bσ1q and coupling coefficients A and B into coefficient Y, equation 70 becomes the following:

Equation 71. Simple strain hardening model with viscosity function integrated and simplified with stress function. The viscoplastic strain level, epsilon subscript vp, equals the integral of stress, sigma, raised to the power of coefficient q from time equal 0 to the time of interest, t, multiplied by the ratio of coefficient p plus 1 divided by coefficient Y, raised to the power of 1 divided by 1 plus coefficient p.       (71)

In the current work, the coefficients, p, q, and Y, are pressure-dependent quantities.

Typically, viscoplastic models are characterized using creep and recovery tests. These tests allow relatively easy separation of the viscoplastic and viscoelastic components, as shown in figure 11. However, it is difficult (if not impossible in some machines) to maintain zero load during the recovery period of the creep and recovery test in tension. Therefore, in tension, viscoplastic characterization uses constant rate tests in which the VECD model is used to first predict the viscoelastic strains. These viscoelastic strains are then subtracted from the measured strains to provide the viscoplastic strains that are needed for curve fitting to equation 71.

The advantage of this model is that it is easy to implement and does not consume much computational time. This model's ability to predict the HMA behavior under complex loading histories in tension is reported by Underwood et al.(12) However, this model's major weakness is its one-dimensional nature, which is not sufficient to describe the behavior of HMA in pavements. This deficiency is particularly troublesome in compression where the confinement is known to play a major role in permanent deformation behavior of HMA. In the following subsections, more general viscoplastic models are presented.

Figure 11. Illustration. Strain decomposition from creep and recovery testing. This figure shows the strain pulse for a creep and recovery test. The axial strain is plotted on the y axis, and time is plotted on the x axis schematically. During the creep portion of loading, strain increases instantly at first, labeled plastic and elastic strains. Then, there is a period of delayed increase, labeled viscoelastic and viscoplastic, to show the time dependent nature of these deformations. Upon recovery, an instant strain reduction is observed and is labeled as elastic strain. The differences between the projected strain increase and the recovered strains are labeled as viscoelastic strains, and the final irrecoverable strains labeled as plastic and viscoplastic.

Figure 11. Illustration. Strain decomposition from creep and recovery testing.

2.5.2. HISS-Perzyna Model

The HISS plasticity model provides a general formulation for the elastoplastic characterization of material behavior. This model, which is potentially able to incorporate isotropic and kinematic hardening and associated and nonassociated plasticity characterizations, can be used to represent a material response based on the continuum plasticity theory. Therefore, the HISS model allows the selection of a more appropriate derivative model for a given material in a specific engineering application. Various well-known plasticity models, such as the Von Mises, Mohr-Coulomb, Drucker-Prager, continuous yielding critical-state, and capped models, can be derived as special cases of the HISS model.(30,32,35)

2.5.2.1. HISS Model Implemented by Delft University of Technology

Equation 72 represents the HISS criterion.(35) In the criterion, n and α determine the shape of the yield stress in deviatoric-hydrostatic stress space, and γ represents ultimate yield stress. R represents the cohesion of the material and determines the position of the yield stress with respect to the hydrostatic stress axis. Because the yield stress of the HISS criterion varies depending on the first stress invariant, I1, the model exhibits a spindle shape of its yield surface when the shape of the yield surface is assumed to be circular (β =0) in the deviatoric stress space, as shown in figure 12.

Equation 72. HISS model implemented by Delft University of Technology. The yield function, lowercase f, is equal to the second deviatoric stress invariant, J subscript 2, divided by atmospheric pressure, P subscript a, squared, multiplied by square bracket softening parameter, gamma, multiplied by parenthesis first stress invariant, I subscript 1, minus the tensile strength of material when deviatoric stress is 0, R, divided by atmospheric pressure, P subscript a, close parenthesis term squared minus hardening parameter, alpha, multiplied by parenthesis the first stress invariant, I subscript 1, minus the tensile strength of material when deviatoric stress is 0, R, divided by atmospheric pressure, P subscript a, close parenthesis raised to the power of n end square bracket, divided by the square root of parenthesis 1 minus beta multiplied by 3 multiplied by the square root of 3 multiplied by the third stress invariant, J subscript 3, divided by 2 divided by the second stress invariant, J subscript 2, raised to power of 3 divided by 2, close parenthesis.       (72)

Where:

γ = Softening parameter.

α = Hardening parameter.

R = Tensile strength of material when deviatoric stress is zero.

n = Parameter determining shape of yield stress.

β = Parameter determining shape of yield stress in deviatoric stress space.

Pa = Atmosphere pressure.

Figure 12. Illustration. Typical yield surface of HISS model. This figure shows one typical yield stress of HISS model, which is represented in three-dimensional, principal stress space. The first principle stress is shown on the x axis; the second principle stress is shown on the y axis; and the third principle stress is shown on the z axis. The figure shows a closed spindle shape protruding into the three dimensional space that the HISS model describes.

Figure 12. Illustration. Typical yield surface of HISS model.

Researchers using this model suggest an HISS criterion whose parameters are the strain rate-dependent functions for a given HMA mixture.(36) For the characterization, a series of constant strain-rate tests in tension and in compression are performed at several strain rates and temperatures. Then, predictions are made for indirect tension (IDT) specimens subjected to constant strain-rate loading. Because this model requires constant strain-rate testing using several different strain rates and temperatures for characterization, it is important to determine the appropriate range of the strain rates and temperatures to minimize the amount of experimental effort. Although it appears that the model can successfully explain the viscoplastic behavior when subjected to monotonic loading, numerical issues related to parameter determination and prediction still remain. In addition, because the model is characterized using constant strain-rate tests, it cannot fully explain the behavior of HMA under discontinuous loading, such as repetitive creep and recovery loading.

2.5.2.2. HISS Model Implemented by the University of Maryland

A viscoplastic model based on the simplified HISS model and using Perzyna's flow rule has also been suggested.(28) Repetitive creep and recovery tests are used for both calibration and prediction, and numerical optimization techniques are adapted for calibrations, unlike Erkens' model.(36) As shown in equation 73, γ, which represents both the ultimate yield stress and the softening of the material, is considered a constant, and R is a function of the viscoplastic strain. In Erkens' study, however, γ and R are functions of the viscoplastic strain rate and temperature.(36)

Equation 73. Definition of yield stress function. The yield stress function, F, is equal to the second stress invariant, J subscript 2, minus square bracket constant softening function, gamma, multiplied by parenthesis first stress invariant, I subscript 1, plus the tensile strength of material when deviatoric stress is 0, R, as a function of the viscoplastic strain trajectory, xi, close parenthesis raised to power of 2 plus the hardening parameter, alpha as a function of the viscoplastic strain trajectory, xi, multiplied by parenthesis first stress invariant, I subscript 1, plus the tensile strength of material when deviatoric stress is 0, R, as a function of the viscoplastic strain trajectory, xi, close parenthesis raise to power of n, square bracket.      (73)

Where:

ζ = Viscoplastic strain trajectory.

γ = Softening parameter (constant).

α = Hardening parameter (function of viscoplastic strain).

R = Tensile strength of material when the deviatoric stress is zero.

The main contribution of this research is to apply the t-TS principle to a conventional viscoplastic constitutive model and confirm the validity of the superposition principle.

2.5.3. Unified Model

In theory of viscoplasticity, the term unified constitutive model refers to models that describe the rate-dependent viscoplastic strain for steel or polymer. However, in other disciplines, the unified constitutive model is used to represent not only viscoplastic strain but also viscoelastic strain. Because this type of model takes a more flexible form of hardening equations than other viscoplastic models, it is worth reviewing the approaches used in the unified model.

2.5.3.1. Linear Kinematic Hardening Model

The simplest unified constitutive model is the linear kinematic hardening model shown in equation 74 and equation 75.(37) The inelastic strain rate has a linear relationship with the overstress, σ - α. The back stress rate is a function of the back stress in the previous time step and the strain rate in the current time step. Because the model is designed to explain inelastic strain, which is the summation of viscoelastic strain, plastic strain, and viscoplastic strain, the inelastic strain rate can be negative depending on the stress history. The back stress rate,Definition of hardening function alpha, is another representation of the Maxwell model of mechanical analog.(50)

Equation 74. Definition of inelastic strain in the simplest unified model. The inelastic strain, epsilon overdot superscript in, is equal to material constant, C, multiplied by parenthesis stress, sigma, minus back stress, alpha, close parenthesis.       (74)

Equation 75. Definition of hardening function alpha. The hardening function, alpha overdot, is equal to coefficient A multiplied by the inelastic strain, epsilon superscript in, minus coefficient B multiplied by back stress, alpha.       (75)

Where:

 epsilon overdot superscript in = Inelastic strain.

α = Back stress.

A,B,C = Material constants.

2.5.3.2. Chaboche Model

In general, the viscoplastic strain of a material obeys a power law, and the hardening of the material can be represented by kinematic and isotropic hardening, as shown in equation 76.

Equation 76. Definition of viscoplastic strain rate. The viscoplastic strain rate, epsilon overdot subscript ij superscript vp, is equal to the product of the strain rate magnitude function, iota parenthesis stress, sigma, comma back stress, alpha, comma isotropic hardening parameter, k, close parenthesis, multiplied by the partial derivative of the stress function, F, with respect to stress tensor, sigma subscript ij. It is also equal to Macauley bracket stress function, F, parenthesis stress, sigma, minus back stress, alpha, close parenthesis minus isotropic hardening parameter, k, divided by viscosity parameter, D, end Macauley bracket raised to the power of n, multiplied by the partial derivative of the stress function, F, with respect to stress tensor, sigma subscript ij.      (76)

Where:

Φ = Magnitude of strain rate.

α = Kinematic hardening function (back stress function).

k = Isotropic hardening function.

The Chaboche model is a viscoplastic model consisting of the above flow rule and hardening function represented by a summation form of back stress.(34) By using decomposed back stress, the model is capable of describing nonlinear hardening with enhanced accuracy for a wider range of viscoplastic strain. Equation 77 shows the back stress function in the Chaboche model; it becomes the Armstrong-Frederick model when n = 1.

Equation 77. Definition of kinematic hardening parameter. The kinematic hardening parameter, alpha, is equal to the summation of the individual back stress functions, alpha subscript i, as a function of parenthesis viscoplastic strain, epsilon subscript vp, close parenthesis, from 1 to n.       (77)

Where:

αi--(εvp) = ith back stress, which is a function of the viscoplastic strain.

2.5.3.3. Krempl and Ho Models

As an advanced and recent form of the linear kinematic hardening model, some researchers have proposed a viscoplastic model based on the overstress concept.(38,39) The constitutive model also includes a description of the time-dependent recoverable strain of the materials. It begins with the assumption that the coefficients of viscoelasticity are nonlinear functions. Equation 78 is the differential equation used to derive the model.

Equation 78. Definition of partial differential equation for Krempl model. The positive, bounded and even function, m, parenthesis stress, sigma, minus strain function, g parenthesis strain, epsilon, close parenthesis, close parenthesis, multiplied by the strain rate, epsilon overdot, plus strain function, g parenthesis strain, epsilon, close parenthesis is equal to stress, sigma, plus positive, bounded and even function k, parenthesis stress, sigma minus strain function g parenthesis strain, epsilon, close parenthesis close parenthesis multiplied by the stress rate, sigma overdot.       (78)

Where:

m(σ-g(ε)) = Positive, bounded and even function.

k(σ-g(ε)) = Positive, bounded and even function.

g(ε) = Odd function of strain.

σ-g(ε) = Overstress.

When the functions m - g(ε)) and k - g(ε)) become constant and when g(ε) is linear in ε, equation 78 reduces to the differential equation of a standard linear solid model that represents viscoelastic behavior. The equation can then be expanded to the regular convolution integrals of linear viscoelasticity. Because equation 79 holds true for slow loading, the relationship between the strain rate and stress can be expressed as equation 80, the simplest form of the viscoplastic model.(40)

Equation 79. Definition of partial differential equation for Krempl model. The material elastic modulus, E, is equal to positive, bounded and even function, m, parenthesis stress, sigma, minus strain function, g, parenthesis strain, epsilon, close parenthesis, close parenthesis, divided by bounded and even function, k, parenthesis stress, sigma minus strain function, g, parenthesis strain, epsilon, close parenthesis, close parenthesis.     (79)

Where:

E = Material elastic modulus.

Equation 80. Definition of partial differential equation for Krempl model. The inelastic strain rate, epsilon overdot, is equal to the stress, sigma, divided by the elastic modulus, E, plus stress, sigma, minus strain function g, parenthesis strain, epsilon, close parenthesis, divided by the material elastic modulus, E, multiplied by bounded and even function, k, parenthesis stress, sigma minus strain function, g, parenthesis strain, epsilon, close parenthesis, close parenthesis.       (80)

More complicated model forms are available for describing hardening, dynamic and static softening, and relaxing, as well as negative and positive rate sensitivities. However, as shown in equation 80, Krempl's constitutive models have neither loading or unloading conditions nor the concept of yield stress because they are derived from the concept of general viscoelasticity. Therefore, it is not appropriate to describe only the viscoplastic behavior of HMA, even though the model gives an idea of the rate dependency of viscoplastic properties under unloading conditions. This potential shortcoming is overcome by introducing Macauley's bracket into the second term of equation 80.(41) Because the model is based on equation 80, the back stress function is related to the kinematic stress function, as shown in equation 83.

Equation 81. Definition of inelastic strain rate. The inelastic strain rate, epsilon overdot superscript in, is equal to coefficient B multiplied by Macauley bracket, parenthesis absolute value of stress, sigma, minus kinematic stress function, H, close parenthesis, minus isotropic hardening function, R, divided by drag stress, D, end Macauley bracket, raised to power of m, multiplied by parenthesis stress, sigma, minus back stress, G, close parenthesis, divided by absolute value of parenthesis stress, sigma, minus back stress, G, close parenthesis.      (81)

Where:

D = Drag stress.

R = Isotropic hardening function.

H = Kinematic stress function Equation 82. Definition of kinematic stress function rate. The kinematic stress hardening rate, H overdot, equals the product of the material’s elastic modulus, E, and the inelastic strain rate, epsilon overdot superscript in.    (82)

G = Back stress.

Equation 83. Definition of resistance to inelastic strain. The resistance to inelastic strain, G overdot, is equal to coefficient psi multiplied by square bracket strain rate, epsilon overdot, minus parenthesis back stress function, G, minus kinematic stress function, H, close parenthesis, divided by isotropic hardening function, R, multiplied by the absolute value of strain rate, epsilon overdot superscript in, end square bracket, plus kinematic stress function rate, H overdot.       (83)

Where:

B, m, and ψ = Material constants characterized from experimental results.

In contrast to Krempl's model, this model does not allow a change in the material state during unloading because the state of the material is a function only of the viscoplastic strain rate.

2.6. Viscoelastoplastic Continuum Damage Model

In the previous sections, the viscoelastic and viscoplastic models are developed separately. It has been shown by using a thermodynamic formulation that the total strain is the sum of the viscoelastic and viscoplastic components even in the non-small strain region.(10) Thus, the viscoelastic and viscoplastic models can now be integrated to form the VEPCD model in which the predicted viscoelastic and viscoplastic responses are combined to obtain the total response for a given stress history.

The resulting equation from combining equation 20 and equation 71 predicts the total strain history for a general loading history at a reference temperature.

Equation 84. Multiaxial viscoelastoplastic continuum damage model. The total strain, epsilon subscript T, is equal to the reference modulus multiplied by the convolution integral of the creep compliance, D, as a function of reduced time, lowercase xi, minus the reduced time when loading began, lowercase xi superscript apostrophe, and the derivative of the stress, sigma, minus pressure multiplied by the second material integrity term, C subscript 12, which is a function of damage, S, divided by the first material integrity term, C subscript 11, with reduced time, lowercase xi superscript apostrophe, with integration between 0 and the time of interest; plus the integral of stress, sigma, raised to the power of coefficient, q, from reduced time equal 0 to the time of interest, lowercase xi, multiplied by the ratio of coefficient p plus 1 divided by coefficient Y, raised to the power of 1 divided by 1 plus coefficient p.       (84)

Equation 84 is presented in terms of reduced time, (ζ), due to the verification of the t-TS principle with growing damage. Stated simply, this principle allows the use of time-temperature shift factors determined from LVE characterization to combine the effects of time and temperature at higher levels of damage. This simplification significantly reduces the required testing protocols while simultaneously expanding the realm of application for the model. Section 4.5 and subsection 5.1.3 present the methodology and experimental work used to verify this principle for both tension and compression modes of loading, respectively.

Previous | Table of Contents | Next

ResearchFHWA
FHWA
United States Department of Transportation - Federal Highway Administration