U.S. Department of Transportation
Federal Highway Administration
1200 New Jersey Avenue, SE
Washington, DC 20590
202-366-4000


Skip to content
Facebook iconYouTube iconTwitter iconFlickr iconLinkedInInstagram

Federal Highway Administration Research and Technology
Coordinating, Developing, and Delivering Highway Transportation Innovations

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 6. Enhancement of VEPCD-FEP++ for Pavement Modeling

A complete redesign of FEP++ was one of the major tasks completed in this study. The redesign activity resulted in a well designed and modular code base of around 40,000 lines. All of the old features of FEP++ have been revamped and tested, and several new features have been added.

Currently, FEP++ supports the following:

  1. Two-dimensional (2D) and 3D pavement models analysis.
  2. Elastic, viscoelastic, and VECD material models.
  3. Separate analysis for static and dynamic loads.
  4. Linear, nonlinear, and quasilinear analysis.
  5. Special elements for modeling pavements.
  6. Elastic and viscoelastic materials thermal analysis.
  7. Multipoint constraints.
  8. Graphical user interface (preprocessor).

The following sections discuss the upgraded VECD model, the overall organization of the modules in FEP++, and the newly developed preprocessor.

6.1. Damage in Viscoelastic Materials

For the sake of completeness, the original VECD model, as detailed elsewhere, is presented first in this section.(13) Then, the upgraded VECD model is presented. Last, the implementation of temperature dependency for these materials is examined.

6.1.1. The Original Model

A one-parameter continuum damage model was incorporated into the viscoelastic material module. The new elements account for the evolution of isotropic damage. In this formulation, it is assumed that the stress-strain relationship obeys the following form:

Equation 181. Stress-strain relationship. The stress, sigma, is given by the product of the damage function, C, times the Young’s modulus, E, times the strain, epsilon.       (181)

Where:

C(S) = Effect of damage on the stiffness of the material.

S = Parameter representing damage evolution.

The following evolution law governs the damage parameter:

Equation 182. Damage evolution law. The partial derivative of the damage parameter, S, with time is given by the negative partial derivative of W, with respect to S raised to the power of alpha.       (182)

Where:

W = Total energy in the material.

α = Parameter that depends on the type of loading.

Currently, based on the available experimental data, the following form for the damage function is assumed:

Equation 183. Damage function. C is given by the exponential of two terms. The first term is a constant, a, and the second term is given by S raised to the power of another constant, b.       (183)

A plot of the above function, assuming a = 0.001334 and b = 0.5725, is shown in figure 173. The formulation presented here is not limited to the form of the damage function shown in equation 183. To show the effect of continuum damage, a 10 by 10 patch of 2D plane stress elements has been modeled and is shown in figure 174. The specimen is fixed at the bottom and loaded (in tension) on five nodes at the top. The loading history is sinusoidal as shown in equation 184.

Equation 184. Loading history. The applied load, f, is given by the product of the loading amplitude, gamma, times the sum of 1 plus sine 2 pi times the time, t.       (184)

Where:

γ = The loading amplitude.

The time history of the vertical displacement for point A () is shown in figure 175. The evolution of damage in the specimen and, consequently, a gradual increment in the displacement amplitude can be seen from the time history. This figure shows that damage has an effect on the response of the material. More quantitative tests on the effect of damage will be performed in the future.

Figure 173. Graph. Damage characteristic relationship used in the finite element implementation of one-dimensional VECD model. This figure shows the relationship of damage function, C, with the damage parameter, S. The x axis is damage parameter parenthesis 0 to 150,000 close parenthesis, and the y axis is the damage function in parenthesis 0 to 1 close parenthesis. This graph is a plot of the relationship in equation 3 with the constants a and b being assigned the values 0.001334 and 0.5725, respectively. The damage function, C, is equal to at 1 when S is equal to 0 and decreases exponentially to 0.3 when S is equal to 150,000.

Figure 173. Graph. Damage characteristic relationship used in the finite element implementation of one-dimensional VECD model.

Figure 174. Graph. Layout of the numerical experiment specimen. This figure shows a discretization of the sample used in the numerical experiment. The discretized specimen is a two-dimensional grid with 10 square elements along either axis. The displacement is measured at a node, A, which can be represented as the 6, 8 node in the two-dimensional grid.

Figure 174. Graph. Layout of the numerical experiment specimen.

Figure 175. Graph. Effect of continuum damage evolution on the vertical displacement of the test specimen at point A in the test simulation. This figure shows the result of the numerical experiment on the specimen shown in figure 109. The x axis is time in seconds from parenthesis 0 to 20 close parenthesis, and the y axis is vertical displacement in meters parenthesis 0 to 0.004 close parenthesis. Two data plots are shown in the graph with the dashed line corresponding to the case when the material is modeled as purely viscoelastic and the solid line corresponds to the case when modeled as VECD. Both plots show the generic trend of increasing displacement over time and have the same frequency. However, the graph clearly shows that the amplitude of the oscillations is higher for the case of VECD model, and the effect is particularly significant at longer times.

Figure 175. Graph. Effect of continuum damage evolution on the vertical displacement of the test specimen at point A in the test simulation.

6.1.2. Upgraded Constitutive Model

The continuum damage formulation was upgraded to a more rigorous model based on the work of others.(57) However, a one-parameter damage model was chosen instead of the two-parameter model used in the work because the experimental data needed to characterize a two-parameter damage model are not yet available.

The model assumes that a material is isotropic when undamaged and that the growth of damage under loading leads to local transverse isotropy (i.e., the material has a local axis of symmetry oriented along the maximum principal stress direction). The current framework is formulated for the axisymmetric case but can easily be extended to 3D.

Starting with the work potential theory for an elastic material and making use of viscoelastic fracture mechanics and the correspondence principle for viscoelastic materials, a pseudo strain energy density function, WR, can be written in terms of the pseudo strains in the local axis as follows: (8,9)

Equation 185. Pseudo strain energy density function in terms of all pseudo strain values. The pseudo strain energy density function, W superscript R, is given by one half of the sum of four terms. The first term is the product of a constant, A subscript 11, times the volumetric strain, e subscript v superscript R. e subscript v superscript R in turn is given by the sum of the pseudo strains along the local axis, epsilon subscript 11 superscript R plus epsilon subscript 22 superscript R plus epsilon subscript 33 superscript R. The second term is the product of a constant, A subscript 22, times the deviatoric strain, e subscript d superscript R. e subscript d superscript R in turn is given by the difference epsilon subscript 33 superscript R minus one third e subscript v. The third term is the product of 2 times a constant, A subscript 12, times e subscript d superscript R times e subscript v superscript R. The fourth term is the product of a constant, A subscript 44, times the sum of two terms given by shear strain along 1–3 plane, gamma subscript 13 superscript R, raised to the power 2 plus shear strain along 2–3 plane, gamma subscript 23 superscript R, raised to the power 2. The final term is a product of a constant, A subscript 66, times the sum of shear strain along 1–2 plane, gamma subscript 12 superscript R, raised to the power 2 and square of e subscript s superscript R. e subscript s superscript R in turn is given by the difference epsilon subscript 11 superscript R minus epsilon 22 superscript R.       (185)

Where:

ε11R, ε22R, ε33R, γ12R, γ13R, γ23R are the pseudo strains along the local axis.

When the local axis is also a principal axis, the shear strains are zero, and equation 185 becomes the following:

Equation 186. Pseudo strain energy density function in terms of principle pseudo strains. W superscript R is given by one half the sum of four terms. The first term is given by the product A subscript 11 times the square of e subscript v superscript R. The second term is given by A subscript 22 times the square of e subscript d superscript R. The third term is given by the product of 2 times A subscript 12 times e subscript v superscript R times e subscript d superscript R. The fourth term is given by A subscript 66 times the square of e subscript s superscript R.       (186)

In this case,ε11R, ε22R, and ε33R are the principal pseudo strains (which lie in the local axis) with the axis of isotropy oriented along direction 3. These principal pseudo strains are obtained from the pseudo strains along the global axis using standard tensor transformation and knowing the angle between the local axis and the global axis. (For the axisymmetric analysis, the hoop direction is already a principal axis because there are no shear strains along the Gamma Theta-Theta Zeta plane.) The LVE pseudo strains along the global axis,εklR(t), are calculated from strains along the global axis, εkl, using the convolution integral.

Equation 187. Pseudo strain along global axis. Pseudo strain along the global axis, epsilon subscript k l superscript R, is given by the product of two terms. The first term is 1 divided by reference modulus, E subscript R. The second term is a convolution integral with the variable of integration, tau, going from 0 to t. The integrand is given by the product E evaluated at t minus tau, times the partial derivative of actual strain along kl, epsilon subscript kl, with tau.       (187)

Where:

ER = Reference modulus that has the same dimensions as the modulus and is usually taken as 1 E(t).

E(t) = Relaxation modulus for uniaxial loading.

The calculation of the convolution integral can be very expensive in terms of computation time. In practice, therefore, the pseudo strains are calculated using a state variable approach to reduce the computational expense. When the relaxation modulus of the material E(t) is represented using the Prony series of the form, shown in equation 188, an approximation can be obtained to the convolution integral in equation 187.

Equation 188. Prony series for relaxation modulus. The relaxation modulus, E, evaluated at t, is given by the sum of two terms. The first term is the relaxation modulus at t equals infinity, E subscript infinity. The second term is a summation over i and the summand is the product of Prony coefficient, E subscript i, times the exponential of the term, negative t divided by relaxation time, rho subscript i.       (188)

Where:

E = The relaxation modulus at t = ∞.

Ei = The Prony coefficients corresponding to the relaxation times, Ρ subscript i.

The pseudo strains along the global axis can be calculated from the following:

Equation 189. Pseudo strain along global axis. The kl component of the pseudo strain along the global axis, epsilon subscript kl superscript R, is given by 1 over E subscript R times the sum of two terms. The first term is the product E subscript infinity times the mechanical strain evaluated at the n plus 1 time step, epsilon subscript kl parenthesis t subscript n plus 1 close parenthesis. The second term is a summation over j from 1 to M of the summand given by the product E subscript i times the ith internal state variable, epsilon subscript kl superscript i evaluated at t subscript n plus 1. Epsilon subscript kl superscript i evaluated at t subscript n plus 1, is in turn given by the sum of two terms. The first term is given by the product of exponential negative quotient of delta t subscript n plus 1 divided by rho subscript i, times epsilon subscript kl superscript i evaluated at t subscript n. The second term is given by the product of exponential negative quotient of delta t subscript n plus 1 divided by twice rho subscript i, times delta epsilon subscript kl evaluated at t subscript n plus 1.      (189)

Where:

εkli(tn+1) (i = 1..M) are the internal state variables that record the history of the material up to time tn+1, Δεkl(tn+1) = εkl(tn+1) - εkl(tn), and Δt=tn+1-tn. The pseudo strains at time tn+1 are calculated based on the strain at times tn and tn+1 which are available at the end of the finite element solution step at time tn+1.

The factors A11, A22, A-12, and A66 are stiffness terms that can be related to a damage function, C(S) as follows: (48)

Equation 190. Constants. A subscript 11 is given by one ninth the sum of C plus the product of E subscript R times the quotient of twice 1 plus Poisson’s ratio, nu, divided by 1 minus twice nu. A subscript 22 is given by C plus the product of E subscript R times the quotient of 1 minus twice nu divided by twice of 1 plus nu. A subscript 12 is given by one third of the difference C minus E subscript R. A subscript 44 is equal to A subscript 66 which in turn is equal to the quotient of E subscript R divided by twice the sum of 1 plus nu.       (190)

Where:

ν = Poisson's ratio of the material.

C(S) = A stiffness function that depends on damage in the material.

S = A damage parameter used to track the growth of damage in the specimen.

The principal stresses along the local axis can be found from equation 186 using the following:

Equation 191. Principal stress. The principal stresses along the local axis, sigma subscript ii where i goes from 1 to 3, is given by partial derivative of W superscript R with pseudo principal strains, epsilon subscript ii.     (191)

which gives the following:

Equation 192. Principal stresses. Sigma subscript 11 is given by the sum of three terms. The first term is given by the product of e subscript v superscript R and the difference A subscript 11 minus one third of A subscript 12. The second term is given by the product of e subscript d superscript R and the difference A subscript 12 minus one third A subscript 12. The third term is given by the product of negative A subscript 66 times e subscript s superscript R. Sigma subscript 22 is given by the sum of three terms. The first term is given by the product of e subscript v superscript R and the difference A subscript 11 minus one third of A subscript 12. The second term is given by the product of e subscript d superscript R and the difference A subscript 12 minus one third A subscript 12. The third term is given by the product of A subscript 66 times e subscript s superscript R. Sigma subscript 33, is given by the sum of two terms. The first term is given by the product of e subscript v superscript R and the sum A subscript 11 plus two third of A subscript 12. The second term is given by the product of e subscript d superscript R and the sum A subscript 12 plus two third A subscript 12. The third term is given by the product of negative A subscript 66 times e subscript s superscript R.       (192)

The stresses along the global axis are then obtained by standard stress transformation and the orientation of the local axis with respect to the global axis.

6.1.3. Damage Model

The growth of the damage parameter S is modeled by extending the concepts of viscoelastic fracture to microcracking.(8)

Equation 193. Damage evolution law. The partial derivative of the damage parameter, S, with time is given by the negative partial derivative of W superscript R, with respect to S raised to the power of alpha.       (193)

Where:

WR = Pseudo strain energy density function shown in equation 185.

α = Material-dependent parameter.

From equation 186, the quantity WR/S can be calculated as a function of pseudo strains in the local axis as follows:

Equation 194. Derivative of pseudo strain energy density function. The partial derivative of W superscript R with S is given by half the product of two terms. The first term is the partial derivative of C with S. The second term is given by the sum of one ninth the square of e subscript v superscript R plus two third the product of the e subscript v superscript R times e subscript d superscript R plus the square of e subscript d superscript R.       (194)

Where:

C(S) = Damage function that is assumed to be of the form shown in equation 195, based on experimental data.

Equation 195. Damage function. C is given by exponential of the product of two terms. The first term is a constant, negative a. The second term is given by S raised to the power of another constant, b.       (195)

6.1.4. Finite Element Implementation

The finite element solution of the problem requires the material tangent stiffness matrix that is used to assemble the global tangent stiffness matrix used for the solution of the nonlinear system of equations by the Newton-Raphson method. Not to be confused with the damage function, C(S), the material tangent stiffness matrix, [C], is given by the following:

Equation 196. Material tangent stiffness matrix. The material tangent stiffness matrix is given by partial derivative of the stress vector with the strain vector. This is also equal to the product of the partial derivative of the stress vector with pseudo strain vector times the partial derivative of the pseudo strain vector with strain vector.       (196)

Where:

{σ} = {σrr, σ θθ, σzz, σrz} and {ε}={εrr, ε θθ, εzz, εrz}are the stresses and strains for the axisymmetric problem.

[C] = Material tangent stiffness matrix oriented along the global axis. Because the stresses along the global axis are obtained by transforming the stresses in the local axis, it is easier to construct the material tangent stiffness matrix along the local axis and then transform it to the tangent stiffness along the global axis.

The tangent stiffness matrix in the local axis, [CL], is given as follows:

Equation 197. Local material tangent stiffness matrix. The material tangent stiffness matrix along the local axis is given by partial derivative of the stress vector with the pseudo strain vector. The matrix is a symmetric matrix and the non-zero elements in the upper triangular region are as follows: the 11 element is the sum A subscript 11 minus two thirds A subscript 12 plus one ninth A subscript 22 plus A subscript 66. The 12 element is given by the sum A subscript 11 minus two thirds A subscript 12 plus one ninth A subscript 22 minus A subscript 66. The 13 element is given by A subscript 11 plus one third A subscript 12 minus two ninth A subscript 22. The 2, 2 element is given by A subscript 11 minus two thirds A subscript 12 plus one ninth A subscript 22 plus A subscript 66. The 23 element is given by A subscript 11 plus one third A subscript 12 minus two ninth A subscript 22. The element 33 is given by A subscript 11 plus four third A subscript 12 plus four ninth A subscript 22. The element 44 is A subscript 66. The element 55 is A subscript 44 and the element 66 is A subscript 44.       (197)

Where:

L} = {σ11, σ22, σ33, σ12, σ13, σ23} is the stress vector along the local axis.

LR} = {ε11R, ε22R, ε33R, γ12R, γ13R, γ23R} is the pseudo strain vector along the local axis.

The pseudo shear strains are zero along the local axis, but the stiffnesses corresponding to them are non-zero. The matrix, [C], can be obtained by transforming [CL], as shown in the following:

Equation 198. Matrix C. The matrix C is given by the product of the inverse of rotation and permutation matrix, T subscript R, times the matrix C subscript L times the matrix T subscript R.       (198)

Where:

[TR] = The rotation and permutation matrix that changes the order of the vector components (the axis 3 along the local axis is always oriented along the maximum principal pseudo strain direction) and transforms a vector from the local axis to the global axis.

As seen in equation 193 and 194, the damage growth involves a nonlinear differential equation that can be expensive to solve in a large finite element method (FEM) problem. Hence, a semi-implicit method is used to predict the damage parameter in the next time-step, Sn+1, using the damage parameter in the current time-step, Sn, and the pseudo strain vector in the local axis, {εLR}n+1, for the next time-step, as follows:

Equation 199. Time stepping of S. S at time step n plus 1, S superscript n plus 1, is given by S at time step n, S superscript n, plus the product of the time step size, delta t, times the time derivative of S evaluated using the pseudo strain vector at time step n plus 1, parenthesis epsilon subscript L superscript R close parenthesis, superscript n plus 1, and S superscript n.       (199)

This method should provide results similar to those of an exact nonlinear analysis when the time-steps are made small enough.

6.1.5. Verification

The continuum damage material model is verified for a monotonic uniaxial test. The test is part of the uniaxial testing on cylindrical specimens that is conducted for the damage characterization of FHWA mixtures. Verification uses the 5 °C test because there is little or no viscoplastic deformation at this temperature. For the FEM simulation, the cylindrical specimen is modeled as a single quadrilateral (Q4) element that is 150 mm in height and 75 mm in width and is subjected to a uniform axisymmetric loading in the vertical direction. The stress history measured from the test is used as an input for the problem, with the strains calculated by analysis. Horizontal displacements are constrained along the axisymmetric axis of the element, and vertical displacements are constrained along the bottom edge of the element, which results in a uniform state of stress and strain in the element. Predicted strains are then compared against measured strains to verify the FEM implementation of the material model. The damage function is taken as C(S) = exp(-0.001S0.5737), based on the monotonic test results from several different specimens. The Prony series for the material is obtained from test data for the asphalt concrete mixture used in the testing.

The results of the verification are shown in figure 176 through figure 178. As seen from these figures, the prediction of strain is very good, the damage parameter consistently increases with time and the pseudo stiffness consistently decreases with time. These combined findings verify that the model has been implemented correctly.

Figure 176. Graph. Verification of strain prediction for monotonic uniaxial test using the new continuum damage material model. This figure shows the validation of the VECD model against experimental data. The x axis is reduced time in seconds parenthesis 0 to 150 close parenthesis, and the y axis is strain parenthesis 0 to 0.0025 close parenthesis. The experimental or measured data is plotted in a magenta color using a solid line plot, and the predicted or simulation results are plotted in a blue color using the plus sign as the data marker. Both plots start at 0 strain for reduced time of 0 and then monotonically increase with reduced time. The two plots show good agreement.

Figure 176. Graph. Verification of strain prediction for monotonic uniaxial test using the new continuum damage material model.

Figure 177. Graph. Evolution of damage parameter, S, for the monotonic test. This figure shows the evolution of the damage parameter, S, over time in a uniaxial monotonic test. The x axis is time in seconds parenthesis 0 to 150 close parenthesis, and the y axis is damage parameter parenthesis 0 to 180,000 close parenthesis. The graph shows that the damage parameter is 0 at initial time 0 and monotonically grows to around 160,000 at time equal to 140 s.

Figure 177. Graph. Evolution of damage parameter, S, for the monotonic test.

Figure 178. Graph. Plot of function C(S) with time. This figure shows the evolution of the damage function, C, with time in a monotonic uniaxial test. The x axis is time in seconds parenthesis 0 to 150 close parenthesis, and the y axis is damage function or pseudo-stiffness parenthesis 0 to 1.2 close parenthesis. The graph shows that the pseudo-stiffness parameter is 1 at initial time 0 and monotonically falls to around 0.3 at time equal to 140 s.

Figure 178. Graph. Plot of function C(S) with time.

6.1.6. Implementation of Temperature Dependency

The variation of temperature in a pavement has two effects: (1) a change in stiffness of the asphalt concrete and (2) a change of the thermal stresses due to thermal expansion of the material. The thermal stresses are generated in the pavement depending on boundary conditions. These two effects of temperature have been implemented in FEP++.

The change in stiffness in the asphalt concrete due to temperature is taken into account using the concept of reduced time. The reduced time is calculated as follows:

Equation 200. Reduced time. The reduced time, eta, is given by the integral over time, t, between the limits 0 to 1 with the integrand given by the quotient of 1 divided by the time temperature shift factor, a.       (200)

Where:

a(T) = Time-temperature shift factor of the material. The time-temperature shift factor is obtained from the characterization of the relaxation modulus of the material using dynamic modulus tests at different frequencies and temperatures. The reduced time is calculated with respect to a reference time and captures the history of temperature variation on the material. The thermal stresses on the material are incorporated by defining mechanical strains, εm(t), as follows:

Equation 201. Mechanical strain. The mechanical strain vector, epsilon subscript m, is given by the material strain vector, epsilon, minus the product of the coefficient of isotropic thermal expansion vector, alpha, times the difference of current temperature, T, minus reference temperature, T subscript 0.       (201)

Where:

{ε (t)} = Strain in the material.

{α} = Coefficient of isotropic thermal expansion in the material (which is assumed constant for asphalt concrete) and is given by [α,α,α,0,0,0] in three dimensions.

T = Current temperature.

T0 = Reference temperature at which there are no initial stresses.

The constitutive law for a viscoelastic material undergoing damage is given as follows:

Equation 202. Stress. The stress vector, sigma, is given by the product of stiffness matrix, D, times the vector of mechanical pseudo strains, epsilon subscript m superscript R.      (202)

Where:

{σ} = Stress vector.

[D] = Stiffness matrix that is a function of damage.

{ mR} = Vector of mechanical pseudo strain.

The ith component of the pseudo strain vector is given by the following:

Equation 203. Mechanical pseudo strain. The ith component of pseudo strain, epsilon subscript Ri superscript m, is given by an integral with the variable of integration, eta prime, going from 0 to eta, of the integrand given by the product E evaluated at eta minus eta prime, times the derivative of epsilon subscript i superscript m with eta prime.       (203)

Where:

εim = ith component of the mechanical strain vector.

E(ζ) = Relaxation modulus of the material characterized at reduced temperature. ζ is computed in accordance with equation 200 and is a function of the temperature history of the material.

6.2. Redesign of FEP++

The older code base of FEP++ was redesigned during the course of the project, leading to a robust, modular, and standards-compliant version. These enhancements greatly improve the maintainability of the code base and significantly reduce the training necessary for a new user.

The FEP++ code base is divided into two major modules, domain and analysis. The following sections discuss the features and functionalities of these modules.

6.2.1. Domain Module

The domain module represents the core finite element component of FEP++. It encompasses the finite element classes that are involved in the modeling of the computational domain. Figure 179 presents a block diagram and the key components of the module. The following sections briefly discuss the main features of the domain module.

Figure 179. Diagram. Domain module. This figure shows the relationship between the different classes in the domain module of FEP++. The domain module is an aggregation of the following classes: Domain, LoadCase, MultiptConstraint, SinglePtConstraint, Material, ElementSet, and NodeSet. The ElementSet class is an aggregation of Element Class, the NodeSet class is an aggregation of Node class and the LoadCase class is an aggregation of NodalLoad class. The Material class is the parent class for Elastic, Viscoelastic and VECD classes. Also, the Element class is a parent class for PlaneStress and PlaneStrain classes. The illustration also shows that all of LoadCase, MultiPtConstraint, SinglePtContraint, and ElementSet classes are aware of the Utilities class.

Figure 179. Diagram. Domain module.

6.2.1.1. Material Model

FEP++ provides 2D and 3D models of elastic, viscoelastic, and VECD materials for modeling a pavement layer. Detailed discussion of these models and their implementation is provided in section 6.1.

6.2.1.2. Elements

FEP++ provides two main types of elements, the quadrilateral elements for 2D analysis and the brick elements for 3D analysis of pavements. Also, special elements developed for modeling infinite domains in pavement analysis have been incorporated into FEP++. The following subsection describes the special elements and their efficiency in pavement analysis.

6.2.1.2.1.    Incorporation of special elements for pavement modeling:

A finite element mesh generation methodology for incorporating special elements that reduce the computational costs of pavement analysis has been developed. An accurate mesh of special elements can be difficult to generate because it depends on the characteristics of the load applied. Development of such a methodology may also turn out to be somewhat computationally expensive if it tries to provide greater accuracy than is necessary for the problem at hand. Hence, there is a need to develop a mesh generation methodology that allows the user to easily generate a finite element mesh that models the problem at hand with specified accuracy.

To this end, work is being carried out to generate the best mesh for a specified accuracy through optimization techniques. Generating the best mesh is an involved process that could render it unfeasible for everyday analysis. Therefore, simplified rules will be developed that allow the user to easily generate a mesh that is as close as possible to the best mesh.

A simple example is given below to demonstrate the computational efficiency achieved. Figure 180 shows a layer of elastic material loaded on a rigid base. This layer is infinite in the two horizontal directions. For a strip load of width "2B," the 3D problem can be reduced to a 2D-plane strain problem. Also, symmetry allows that only one half of the problem needs to be modeled. Only the deflection and stresses under the load (the interior region) are of interest here.

Figure 180. Illustration. Infinite elastic layer on a rigid base. This figure shows a two-dimensional plane strain model of the pavement. The wheel load is applied as a strip load of width “2B” at the center of the pavement. This width of the pavement is called the interior region and is flanked at either side by the exterior region which is an extension of the pavement to infinity. The exterior region models the fact that in reality the pavement has infinite regions surrounding it.

Figure 180. Illustration. Infinite elastic layer on a rigid base.

However, despite the emphasis on the interior region, the exterior region also must be modeled to a sufficient length to simulate the infinite extent of the layer. As the exterior region is extended, the solution in the interior approaches the solution of the physical problem. The objective is to model the physical problem to achieve an accuracy of at least 1 percent. A typical mesh of finite elements used to simulate this physical problem is Mesh-1, shown in figure 181. For Mesh-1, the exterior is four times the extent of the interior and hence needs four times as many finite elements as the interior to achieve an accuracy of 0.86 percent.

If, on the other hand, special elements are used, then Mesh-2 shown in figure 181 is sufficient and provides an accuracy of 0.84 percent. The total number of elements in Mesh-1 is 320, whereas the total number in Mesh-2 is 64. In order to achieve almost the same degree of accuracy as Mesh-1, far fewer special elements are required for Mesh-2. In short, because the computational time is directly proportional to the number of finite elements used to model the problem, Mesh-2 with the special elements is far more computationally efficient when compared to Mesh-1.

The development of a finite element mesh generation methodology that incorporates special elements is now complete. This methodology is being tested on various finite element models to demonstrate its efficacy. Once the results have been verified against theoretical solutions, the mesh generation methodology will be coded into FEP++ to make the discretization both automatic and efficient. (This task will be performed as part of the ongoing cooperative agreement between NCSU and the FHWA.)

Figure 181. Illustration. Finite element mesh required to model the physical problem with 1-percent accuracy without special elements. This figure shows a schematic of the finite element mesh that is needed to model the pavement described figure 180. This mesh is a regular finite element mesh with 8 elements along depth and 40 elements along the width. The interior region is discretized using 5 elements and the exterior region using 35 elements along the width. The bottom and the right edges are fixed and the top surface is free. Symmetry boundary conditions are applied on the left edge. The wheel load of width “B” is applied as nodal loads on the elements along the top edge of the interior region.

Figure 181. Illustration. Finite element mesh required to model the physical problem with 1-percent accuracy without special elements.

Figure 182. Illustration. Finite element mesh required to model the physical problem with 1 percent accuracy with special elements. This figure shows a schematic of the finite element mesh that is needed to model the pavement described in figure 180. This mesh differs from that shown in the previous figure in the way the exterior is discretized. The exterior region in mesh 2 is discretized using only 24 elements. These are special elements designed to model the effect of an infinite exterior and are different from the elements used to model the interior. Thus, the total number of elements in mesh 2 is only 64 and hence is computationally efficient than mesh 1.

Figure 182. Illustration. Finite element mesh required to model the physical problem with 1-percent accuracy with special elements.

6.2.2. Analysis Module

The analysis module handles the generation and solution of the linear/nonlinear system of equations using the domain module. A block diagram of the classes and relationships is shown in figure 183. The following sections discuss the key components of this module.

6.2.2.1. Solution Algorithm

The solution algorithm is the component that acts as a liaison between the domain and analysis modules. This component uses the information in the domain to generate a set of equations that provides the solutions to the problem. Currently, the following static and dynamic analyses of linear, quasilinear, and nonlinear problems are available.

One significant enhancement over the earlier capabilities of FEP++ is the ability to perform quasilinear analysis. This improvement is particularly beneficial in the context of pavement analysis because the VECD material model is quasilinear. By using a quasilinear solver instead of a nonlinear solver, the analysis time is reduced by almost half.

Figure 183. Diagram. Analysis module. The figure shows the class relationships in the analysis module in FEP++. The Analysis class is shown to be an aggregation of AnalysisUtilities and SolutionAlgo classes. The SolutionAlgo is the parent class for LinearStatic, LinearDynamic, QuasiLinearStatic, QuasiLinearDynamic, NonLinearStatic, and NonLinearDynamic classes. Also, SolutionAlgo class is an aggregation of Stepper and Iterator classes. The Stepper is the parent class for Static and Newmark classes, and the Iterator is the parent class for Linear, NewtonRaphson, and ModifedNewton classes. The illustration also shows that the Iterator and Stepper classes are aware of each other, and both SolutionAlgo and Solver classes are aware of the Mapping class.

Figure 183. Diagram. Analysis module.

6.2.2.2. Solver

The solver is the workhorse of the module. It is responsible for solving the set of equations generated by the solution algorithm. The solver is also an important component in terms of the capabilities of FEP++ in solving large computational problems because the solver determines the memory requirements and performance output. Any improvements to this module would significantly improve the overall capabilities of FEP++.

Currently, a class of direct methods called Gaussian elimination is used to solve the equations. The objective here is to enhance the component by providing an option to interface with external linear algebra packages to solve the system of equations. The benefit of such a feature would be the ability to use mature, well tested, and highly efficient libraries available in the public domain for performing the analysis. This capability would also reduce further interdependence between the domain and the analysis modules, thereby enhancing the flexibility and elegance of the code.

6.3. PreProcessor

The preprocessor is a graphic user interface (GUI) front end developed for the efficient analysis of pavements using FEP++. The motivation behind the tool is to provide the user with a simplified and intuitive interface to FEP++ for pavement analysis. The tool can be used either to run an analysis directly or to generate the input files for FEP++. It also has visualization capabilities to view the mesh discretization, which makes it easy for the user to verify the input data for the analysis. Figure 184 shows a screenshot of the main window in the preprocessor.

Figure 184. Screen capture. Main window of the preprocessor. This figure shows the main screen of the preprocessor for FEP++. It is a typical user interface on a windows platform and has a standard menu and toolbar. The body of the window is divided into three regions. The left region is called Control Panel; the middle region is called Canvas Panel; and the right region is called VTK Control Panel.

Figure 184. Screen capture. Main window of the preprocessor.

As labeled in figure 184, the main window is divided into three functional units, the control panel, the canvas panel, and the popular open-source library called Visualization Toolkit (VTK) control panel, described below.

The left-hand unit is called the Control Panel. This panel contains the controls to invoke the data-entry windows and other operations related to the analysis. Each control is either red or green, indicating the current state of the analysis. For controls corresponding to data entry operations, green indicates that no user action is required. Green, therefore, could either signal that the data entry has already been completed or that it is preset with a default set of parameters. Red indicates that some action is required from the user. The data input is considered complete only when all the controls in the input parameters block (figure 185) are green.

Figure 185. Screen capture. Control panel. This figure shows the Control Panel in the main screen of the preprocessor. The panel shows a hierarchy of controls that can be clicked to enter data corresponding to the control. The top level hierarchy has three levels: Input Parameters, Output Settings, and Analysis Results. Input Parameters has six controls: General Information, Material Properties, Layer Properties, Load Properties, Mesh Properties, and Analysis Options. Output Settings has four controls: Output Options, View Summary, Generate Input File, and Run Analysis.

Figure 185. Screen capture. Control panel.

The middle unit is called the Canvas Panel. This panel is used for the visualization of the mesh discretization and the solid model of the pavement. The underlying graphical engine used to render the graphics is VTK. This toolkit was chosen for its stable and mature set of high-level graphical libraries. Also, it is popular and widely used in the open-source and scientific community. Figure 186 shows a sample mesh discretization in the canvas panel.

The rightmost unit provides the controls for the visualization operations supported in the preprocessor, e.g., zoom, pan, etc., and is called the VTK control panel. The following operations are currently supported:

  • Projections along a plane, e.g., along the x-y plane.
  • Image operations, such as zoom, pan, and rotate.
  • Display nodes and cell numbers on the mesh.

The results of these operations are displayed in the view window. Figure 187 shows the zoomed version of a 2D mesh with node numbers displayed. The visual feedback from these windows is useful to verify the data inputs to the analysis. The different stages involved in performing a pavement analysis and the ways that the various analytical pieces fit together are outlined in figure 188.

The input stage corresponds to the data-entry windows that receive the input data from a user (i.e., it shows the list of data that needs to be entered by the user). The analysis stage shows the list of options available to the user once the input stage is complete. The postprocessing stage again shows the options after the analysis is completed successfully. The above mentioned stages are discussed in detail in the following sections.

Figure 186. Screen capture. Mesh discretization. This figure shows the Canvas panel of the main screen and shows the basic visualization capabilities of the preprocessor. The canvas panel has the visualization of a three-dimensional finite element mesh.

Figure 186. Screen capture. Mesh discretization.

Figure 187. Screen capture. Zoom operation on a mesh. This figure is another screen capture of the canvas panel and shows the zoom feature in the preprocessor. The canvas panel has a visualization of a zoomed region of the mesh with the node numbers displayed.

Figure 187. Screen capture. Zoom operation on a mesh.

Figure 188. Diagram. Stages involved in an FEP++ analysis. This figure shows the different stages involved in performing a finite element analysis using FEP++. The analysis is split into three main stages: Input Stage, Analysis Stage, and Post Processing Stage. The figure also shows the sequence in which these stages are completed, with the Input Stage being completed first, then the Analysis Stage, and last the Post Processing Stage. Each of these stages has different operations to be performed. The Input Stage involves data entry operations for the following data: General Information, Material Properties, Layer Properties, Load Properties, Mesh Parameters, and Analysis Options. The Analysis Stage involves the following operations: Data entry for Output Options, View Summary, Generate Input File, and Run Analysis. The Post Processing Stage involves viewing the results of the analysis.

Figure 188. Diagram. Stages involved in an FEP++ analysis.

6.3.1. Input Stage

As discussed in the previous section, this stage corresponds to the data entry by the user. Each main input data entry for FEP++ corresponds to a data entry window. The user manual, with a detailed explanation of each of the data entry windows, is provided in the appendix of this report. Figure 189 shows a data entry window for entering the mesh discretization parameters for a pavement.

The data entry windows have been designed in a modular way, and effective error validation schemes have been incorporated to ensure consistent input data. Each data entry window has a data class associated with it that enforces the validation and consistency checks. There are two levels of error validations that are carried out at each data entry window, namely syntactic and semantic errors. The syntactic errors refer to invalid data entries, and they are checked by the data entry windows themselves, whereas the semantic errors refer to consistency check failures imposed by the data classes.

Once the data entry is complete for a particular block, it is validated against the checks provided in the data classes. If the validation is successful, the next data entry window can be accessed. If it fails, appropriate error messages are shown, and the user is asked to modify the incorrect data. Figure 190 shows a typical error dialog that is displayed when a consistency check fails.

Figure 189. Screen capture. Sample data entry window. This figure shows a sample data entry window that is invoked when the controls in the Control Panel are clicked. The data entry window shown accepts the finite element mesh discretization data.

Figure 189. Screen capture. Sample data entry window.

This sequential data entry process may seem to be somewhat limiting, but it greatly simplifies the maintenance of consistency in the input data. Also, by strictly enforcing these checks in the software, there is less scope for input data corruption.

6.3.2. Analysis Stage

Once the data entry has been completed, the analysis section in the control panel turns green, signaling that the actions in the analysis stage can now be performed. The following sections discuss the actions that can be performed at this stage.

6.3.2.1. Setting the Output Options

Output options can be used to control the results that are generated by FEP++. For example, FEP++ can be configured to output data that correspond to a particular region in the domain. These options are essentially used to control the amount of data output by FEP++.

6.3.2.2. Viewing a Summary of the Input Data

The preprocessor can be used to create a user-readable report of the input data entered for an analysis. This report can be used to verify the final state of all the input data before beginning a simulation.

Figure 190. Screen capture. Error dialog for a semantic error. This figure shows a sample data entry window performing data validation. It shows that an error popup window is launched when invalid data is entered in the data entry window.

Figure 190. Screen capture. Error dialog for a semantic error.

6.3.2.3. Generating an FEP++ Input File

The preprocessor can be used to generate input files for FEP++ for different configurations that can then be saved to a remote location. This feature is helpful if the user wishes to perform multiple analyses in a script without any user intervention.

6.3.2.4. Running an Analysis

Running an analysis is the final step in pavement analysis. The analysis can be started only after the run analysis option in the control panel turns green. The preprocessor uses the user input to create an input file to FEP++ and invokes it with the generated input file. The preprocessor creates a separate instance and displays the progress of the analysis in a window, as shown in figure 191 .

Figure 191. Screen capture. Analysis run-time window. This figure shows a the run-time window. It shows a DOS window with run time information showing the current status of the analysis.

Figure 191. Screen capture. Analysis run-time window.

6.3.3. Postprocessing Stage

After an analysis is completed, the run-time window shown in figure 191 exits, and the results option on the control panel turns green, indicating that result files have been generated that can be viewed with a postprocessor. On choosing view results, the default postprocessor is invoked to view the results.

6.4. PostProcessor

The postprocessor is an external tool for visualizing the results of the pavement analysis. Currently, Tecplot® serves as the default postprocessor, but efforts to provide support for other open-source data visualization tools are underway. It has been decided that the VTK format will be used for the output file format and that MayaVi, a stable open-source python-based viewer that is popular in the finite element community, will be used as a suitable viewer. The support for external tools will be completed in the next project.

Previous | Table of Contents | Next