Davids - EverFE theory Manual (2003)

Prévia do material em texto

EverFE Theory Manual 
 
 
 
 
Bill Davids, Ph.D., P.E. 
University of Maine 
Dept. of Civil and Environmental Engineering 
 
February 2003 
 
EverFE version 2.23
 1
1. Introduction 
EverFE (current version 2.23) is a 3D finite-element analysis tool for simulating the response of 
jointed plain concrete pavement systems to axle loads and environmental effects. EverFE couples a highly 
interactive graphical user interface for model development and result visualization written in 
Tcl/Tk/Tix/vTk with finite-element code written in object-oriented C++. 
Some significant features of EverFE include: 
• The ability to model 1, 2, or 3 slab and/or shoulder units longitudinally and/or transversely (up to 
9 slab-shoulder units total in a 3x3 configuration). Transverse tie bars between adjacent slab-
shoulder units can be explicitly modeled. 
• Up to three elastic base layers with a bonded or unbonded base can be specified. Slab-base shear 
transfer can be captured via an elastic-plastic distributed horizontal stiffness between the slabs 
and base. A tensionless or tension-supporting dense liquid foundation underlies the bottom-most 
model layer. 
• Linear or nonlinear aggregate interlock shear transfer can be simulated at transverse joints. 
• Dowels can be precisely located across transverse joints, and dowel looseness modeled. In lieu of 
modeling dowel looseness, a dowel-slab support modulus can be specified to model dowel-slab 
interaction. 
• Dowel misalignment and mislocation can be modeled. 
• A variety of different axle configurations can be easily defined with a minimum amount of input. 
• Linear, bilinear, and trilinear thermal gradients through the slab thickness can be captured. This 
allows the simulation of thermal effects as well as slab shrinkage. 
• EverFE’s extensive post-processing capabilities permit the visualization of stresses, 
displacements, and internal dowel forces and moments. Critical response values at any point in 
the model can be easily retrieved. 
This manual details the finite-element implementation of EverFE, and includes descriptions of critical 
features and references to appropriate literature that augment the material presented here. For specific 
instructions on the use of the software, see the help file EverFE_help.chm, which is integrated with the 
Help menu in the EverFE interface. This manual is organized into the following topics: Finite-Element 
Discretization; Modeling of Dowels and Ties; Aggregate Interlock Modeling; Treatment of Axle and 
Thermal Loads; Finite Element Solution Strategies; and Program Architecture and File Structure. 
2. Basic Finite-Element Discretization 
There are five elements in EverFE’s finite-element library: 20-noded quadratic brick elements are 
used to discretize the slab and elastic base and sub-base layers; 8-noded planar quadratic elements 
incorporate the dense liquid foundation below the bottom-most elastic layer; 16-noded quadratic interface 
elements implement both aggregate interlock joint shear transfer and shear transfer at the slab-base 
interface; and 3-noded embedded flexural elements are coupled with conventional 2-noded shear beam 
elements to model dowels at transverse joints and ties at longitudinal joints. Figure 1 shows a finite-
element mesh of a four-slab model and the corresponding elements. The flexural elements used to model 
dowels and ties are detailed separately in Section 4, and aggregate interlock modeling is covered in 
Section 5. 
 
 
 2
 
 
Figure 1: Basic Finite-Element Discretization 
 
2.1 Model boundary conditions 
The boundary conditions are the minimum required to prevent rigid-body motion, but differ slightly 
depending on whether or not an elastic base layer is explicitly modeled. In the case where a base layer is 
modeled, the slabs are restrained in the horizontal x-y plane by the shear stiffness of the slab-base 
interface as discussed in Section 3, and receive vertical support from contact with the base. Rigid body 
motion of the base and sub-base layers is prevented by restricting the x- and y- displacements of one node 
on the –x face, and restricting the x-displacement of a second node on the –x face. Vertical support of the 
entire system is provided by the dense liquid foundation, which is always incorporated below the bottom-
most layer of the model. 
If the slabs are founded directly on a dense liquid, i.e. no base layer is modeled, each slab is 
restrained against x- and y-direction displacements at one node on its –x face, and against x-direction 
displacement at a second node on its –x face to prevent rigid-body motion of each slab. Again, vertical 
support is provided by the dense liquid foundation. 
2.2 Modeling of the slab, base and sub-base layers 
In all EverFE models, the slab, base and sub-base layers are treated as 3D, linearly elastic, isotropic 
continua. Each layer is discretized with standard 20-noded “serendipity” brick elements. The finite-
element meshes are rectilinear, and the same number of element divisions is used for each slab and the 
base/sub-base layers below the slab in the x-y plane to ensure compatibility at the slab-base interface. 
 Details of the brick element formulation and implementation can be found in finite-element texts 
such as Zienkiewicz and Taylor (1994). To maintain generality, an isoparametric element formulation is 
20-noded brick element
zero
thickness
16-noded interface element
8-noded dense liquid element
x
y
z
20-noded brick element
zero
thickness
16-noded interface element
8-noded dense liquid element
x
y
z
 3
used and all required element integration is performed numerically using 8-point (2x2x2) Gauss 
quadrature. The initial public release of EverFE (version 1.02, released January 1998) used 27-point 
(3x3x3) Gauss integration; however, subsequent internal studies showed that the higher-order integration 
scheme added to the computational time without significantly improving accuracy. 
2.3 Modeling of the dense liquid foundation 
The dense liquid foundation can either support tension, or be tensionless. It is important to note that 
the tensionless dense liquid does not account for any pre-compression due to dead load, i.e. the total 
vertical deflection including the effect of dead load must be overcome before the dense liquid foundation 
stress and stiffness become zero. 
The 8-noded element illustrated in Figure 1 is used to discretize the dense liquid. This element was 
formulated specifically for this application, and full details of the implementation can be found in Davids 
(1998). The element incorporates standard quadratic shape functions for interpolation of vertical 
displacements within the element (Zienkiewicz and Taylor 1994), ensuring that it displaces compatibly 
with the 20-noded brick element with which it shares nodes. An isoparametric element formulation is 
used, and all necessary element integrations are performed numerically using 9-point (3x3) Guass 
quadrature to ensure accurate results when the tensionless option is selected. 
The only constitutive parameter needed for this element is the distributed stiffness of the dense liquid 
foundation [force/volume]. For the tensionless foundation, if tension occurs at an element integration 
point during the solution process, the stress and stiffness at that point are set to zero during integration of 
the element stiffness matrix and equivalent force vector. For the conventional, tension-supporting dense 
liquid, the stiffness remains constant at all points. 
3. Treatment of the Slab-Base Interface 
Modeling interaction of the slab and base is crucial to predicting pavement response to axle loads 
near joints and to thermal or shrinkage gradients. EverFE allows the consideration of either perfect bond 
between the slab and base, or separation of the slab and base under tension. In both cases, the slab and 
base do not share nodes, and nodal constraints are used to satisfy the required contact conditions(see 
Figure 2). The solution algorithm relies on a perturbed Lagrangian formulation and a nodal constraint-
updating scheme based on the current normal stress between the slab and base. More details on the global 
nonlinear solution strategy are given in Section 7. 
Shear transfer between the slab and base can be important when analyzing pavements subjected to 
thermal and/or shrinkage strains. Rasmussen and Rozycki (2001) overviewed the factors governing slab-
base shear transfer, noting that both friction and interlock between the slab and base play a role. In 
addition, they calibrated a bilinear, elastic-plastic shear transfer model from results of push tests of slabs 
Base elementInterface
elements transfer
shear stress
Slab element
Pairs of nodes vertically constrained
if compression at interface
x or y
z
x or y
kSB
o
o
Interface constitutive relationship
Figure 2: Modeling Separation and Shear Transfer at the Slab-Base Interface
Base elementInterface
elements transfer
shear stress
Slab element
Pairs of nodes vertically constrained
if compression at interface
x or y
z
x or y
kSB
o
o
Interface constitutive relationship
Figure 2: Modeling Separation and Shear Transfer at the Slab-Base Interface
 4
on various base types. One conclusion of the their study was that the effect of slab-base shear transfer 
should be incorporated in 3D analyses of pavement systems. Another study by Zhang and Li (2001) 
focused on developing a one-dimensional analytical model for predicting shrinkage-induced stresses in 
concrete pavements that accounts for slab-base shear transfer. Like the model developed by Rasmussen 
and Rozycki, their model ultimately relied on a bilinear, elastic-plastic shear transfer model. Zhang and Li 
concluded that the type of supporting base – and thus the degree to which it restrains slab shrinkage – 
significantly affects slab stresses. 
To capture slab-base shear transfer, EverFE employs 16-noded, zero-thickness quadratic interface 
elements that are meshed between the slab and base as shown in Figures 1 and 2. The element 
incorporates standard quadratic shape functions for interpolation of displacements (Zienkiewicz and 
Taylor 1994), ensuring that it displaces compatibly with the 20-node brick elements with which it shares 
nodes. The element tracks relative displacements between the slab and base in the vertical (z) and both 
horizontal (x and y) directions. An isoparametric element formulation is used, and all necessary element 
integrations are performed numerically using 9-point (3x3) Guass quadrature. 
The bilinear element constitutive relationship is based on that given by Rasmussen and Rozycki 
(2001) and Zhang and Li (2001). Figure 2 illustrates the constitutive relationship, which is characterized 
by an initial distributed stiffness kSB�>IRUFH�YROXPH@�DQG�VOLS�GLVSODFHPHQW� o. While kSB has the same units 
as the well-known dense liquid foundation modulus, kSB is a distributed stiffness in the x- and y-directions, 
and the shear stresses in the x-y plane at the slab-base interface are caused by relative horizontal 
displacements between the slab and base layer. This constitutive relationship is assumed to apply 
independently in both the x- and y-directions as long as the slab and base remain in contact (i.e. a 
compressive normal stress exists at the slab-base interface). The fact that there will be little or no shear 
transfer when slab-base separation occurs is accommodated by setting the interface stiffness and shear 
VWUHVV�WR�]HUR�ZKHQHYHU�WKH�UHODWLYH�YHUWLFDO�GLVSODFHPHQW� z > 0. Modeling this loss of shear transfer with 
loss of slab-base contact can be important, especially when thermal gradients are simulated (Davids et al. 
2003). 
Note that unlike a frictional model, the shear stress does not depend on the magnitude of the normal 
stress. However, for very large values of kSB, this model approaches Coulomb friction with a very large 
friction coefficient, and for very small values of kSB, it is equivalent to a frictionless interface. An 
advantage of the modeling scheme used by EverFE is that the symmetry of the system stiffness equations 
is maintained, which allows the use of the highly efficient preconditioned conjugate-gradient solver 
discussed in Section 7. Idealizing slab-base interaction with conventional Coulomb friction would destroy 
this symmetry, requiring the use of more complex (and likely less efficient) solution techniques. 
4. Modeling of Dowels and Ties 
EverFE models dowels and transverse tie bars explicitly with 3-
noded, quadratic embedded flexural finite elements (Davids and 
Turkiyyah 1997; Davids 2000). This approach has the advantage of 
allowing the dowels and tie bars to be precisely located irrespective 
of the slab mesh lines as shown in Figure 3. This embedded element 
formulation also permits significant savings in computation time by 
allowing a range of load transfer efficiencies to be simulated without 
requiring a highly refined mesh at the joint. Additionally, the 
embedded dowel formulation allows the rigorous treatment of gaps 
between the dowels and the slabs (dowel looseness). Alternatively 
the dowels can be modeled as beams on an elastic foundation, where 
Winkler foundation springs are sandwiched between the dowels and 
the slabs. These two options for capturing dowel-slab interaction are 
shown conceptually in Figure 4(a). 
embedded dowel
element 
solid
element 
Figure 3: Dowel Element
embedded dowel
element 
solid
element 
Figure 3: Dowel Element
 5
 
4.1 Embedded Dowel Element Formulation 
The finite-element formulation of the embedded element relies on expressing the nodal displacements of 
the dowel as functions of the nodal displacements of the solid element it is embedded within. The three-noded 
quadratic beam element (that originally has 18 degrees-of-freedom in its element displacement vector, Ud) 
takes on the degrees-of-freedom of the solid element it is embedded within, Ue. The element displacement 
vector of the embedded bending element, Ude, becomes: 
)1(





=
e
d
de
U
U
U 
As shown in Davids and Turkiyyah (1997), Ud can be recovered through the following matrix 
transformation: 
)2(ded TUU = 
The matrix T contains shape functions of the embedding element and constrains the dowel to the 
embedding element. It follows that the tangent stiffness matrix of the embedded element Kde, needed for 
the nonlinear solver, can be determined from the original dowel element stiffness, K, as follows: 
)3(KTTK T=de 
If dowel looseness is explicitly modeled, a nodal contact approach is employed where T encapsulates 
all the necessary information regarding constraints. The advantage of this approach is that the contact 
nonlinearity resulting from the rigorous simulation of gaps between the dowels and slabs is treated like a 
x∆x
y
Plan View
α
Original
position
Misaligned
position
x
z
Elevation View
z∆
Original
position
Misaligned
position
β
(b) Dowel Misalignment
q
r
qs
Slab
Dowel-slab springs
C.L. joint
(a) Dowel-Slab Interaction
Gap between
dowel and slab
Figure 4: Dowel Modeling
x∆x
y
Plan View
α
Original
position
Misaligned
position
x
z
Elevation View
z∆
Original
position
Misaligned
position
β
(b) Dowel Misalignment
q
r
qs
Slab
Dowel-slab springs
C.L. joint
(a) Dowel-Slab Interaction
Gap between
dowel and slab
x∆x
y
Plan View
α
Original
position
Misaligned
position
x
z
Elevation View
z∆
Original
position
Misaligned
position
β
(b) Dowel Misalignment
q
r
qs
Slab
Dowel-slab springs
C.L. joint
(a) Dowel-Slab Interaction
Gap between
dowel and slab
Figure 4: Dowel Modeling
 6
material nonlinearity through the simple stiffness matrix transformation of Equation 3. Internally, the 
contact conditions are checked and updated at each iteration of the nonlinear solver. 
With the foregoing formulation, the embedded bending element has also been extended to permit the 
inclusion of generalbond-slip relations between the dowel and surrounding slab (Davids 2000). If the 
incremental vector of relative displacements between the slab and dowel is denoted as d , and it is 
assumed that the corresponding incremental force vector df can be computed as: 
)4(∆= dd Df 
In Equation 4, D is a 3x3 constitutive matrix. It has been shown that the stiffness matrix of the embedded 
dowel element with general bond-slip relations, Kdt, can then be computed as: 
)5(
1
1∫−+= ηhd
dedt DBBKK T 
In Equation 5, B is a matrix operator containing shape functions of both the embedding solid element and 
the dowel element, and h is the length of the dowel element; integration is performed with respect to the 
dowel element local coordinate, . 
Physically, this formulation of the embedded element with a general bond-slip relationship between 
the dowel and slab is analogous to the classic beam on elastic foundation, but differs in that forces in the 
three coordinate directions can exist between the dowel and the slab. The magnitude of the forces depends 
on the relative displacement between the dowel and the slab and the constitutive relations of the 
dowel/slab interface incorporated in the matrix, D. 
4.2 Implementation in EverFE 
In EverFE, the 3-noded quadratic embedded dowel elements are used to model the portions of the 
dowels embedded in the slabs. To ensure accurate results, 12 embedded elements are used to model the 
embedded portion of the dowels on each side of each transverse joint, and a 2-noded shear beam is used 
to model the portion of the dowel spanning the joint. 
When dowel looseness is modeled, 10 of the 12 embedded dowel elements are used over the portion 
of the dowel where there is a gap between the dowel and slab to ensure a sufficient number of potential 
points of contact between the dowels and slabs. If the dowel is treated as a beam on a dense liquid 
foundation, the 12 elements are evenly distributed along the embedded portion of the dowel. The diagonal 
terms of D corresponding to the model y- and z- coordinates, D(2,2) and D(3,3), are the dowel-slab 
support modulus specified in EverFE. The dowel-slab support modulus is computed as the product Kd, 
where K is the commonly used modulus of dowel support [force/volume], and d is the dowel diameter. 
The paper by Dei Poli et al. (1992) discusses the basis for the development of K from the properties of the 
slab concrete and dowel; additionally, the EverFE help manual includes the results of a parametric study 
showing the effect of the parameter Kd on load transfer efficiency for a simple two-slab model. D(1,1) 
applies in the x-direction, and is the dowel-slab restraint modulus that controls the degree of bond 
between the dowels and slabs. A large value of the dowel-slab restraint modulus implies a high degree of 
bond between the dowels and slabs, and a value of 0 implies no bond. For an example illustrating the 
effect of this parameter, see Davids et al. (2003). 
Transverse ties are modeled in the same manner as the dowels, although only the dense liquid support 
option is available, and fewer elements are used to discretize the ties since their shears and moments are 
of secondary interest. 
4.3 Simulation of Dowel Misalignment/Mislocation 
EverFE also allows the simulation of dowel misalignment and/or mislocation through the specification of 
IRXU�SDUDPHWHUV�� x�� z�� �� ��WKDW�VKLIW�DQ�LQGLYLGXDO�GRZHO�DORQJ�WKH�x- and z-axes and define its angular 
misalignment in the horizontal and vertical planes (see Figure 4(b)). The dowel support and restraint moduli 
coincide with the local dowel coordinate axes (q,r,s), which are rotated from the global (x,y,z) axes by the 
 7
DQJOHV� �DQG� ��7KH�PHVKLQJ�DOJRULWKP�SUHFLVHO\�ORFDWHV�LQGLYLGXDO�IOH[XUDO�HOHPHQWV�ZLWKLQ�WKH�VROLG�
elements by first solving for the intersection of each dowel with solid element faces, and then subdividing 
each dowel into at least 12 individual quadratic embedded flexural elements on each side of the joint face as 
discussed previously. 
5. Aggregate Interlock Modeling 
Aggregate interlock shear transfer is assumed to occur across the entire width of each transverse joint 
in the finite-element model. Both linear and nonlinear options are available for modeling aggregate 
interlock. With the linear option, the shear stress developed between the joint faces is proportional to the 
relative vertical movement at the joint, and the shear stress is independent of the joint opening. The 
nonlinear option includes both the nonlinearity in the shear stress-relative vertical displacement relation 
as well as the nonlinear variation in shear stress transfer with changes in joint opening. The basis for and 
implementation of both of these options is detailed in this section. 
In both cases, EverFE employs a 16-noded, zero-thickness quadratic interface element that is meshed 
between the joint faces as shown in Figure 1. This is the same element detailed in Section 3 of this 
manual that is used to capture shear transfer at the slab-base interfaces. As discussed in Section 3, an 
isoparametric element formulation is used, and all necessary element integrations are performed 
numerically using 9-point (3x3) Guass quadrature. 
5.1 Linear Aggregate Interlock Load Transfer 
The linear option is the simplest approach for modeling aggregate interlock load transfer at 
longitudinal joints. In this case, only a joint stiffness is specified to control the degree of aggregate 
interlock load transfer. The units of this parameter are force/volume, and it is analogous to a dense liquid 
k in that it can be interpreted as a spring stiffness per unit area. 
The specified joint stiffness is constant over the entire area of the joint, and does not vary with 
relative vertical displacement or joint opening. If the joint stiffness is set to zero (the default value), there 
will be no aggregate interlock load transfer, and a very large value will result in high load transfer 
efficiency. The joint stiffness applies only in the vertical (z) direction, and y-direction relative joint 
movement is unrestrained. 
It is worthwhile noting that a number of prior studies have used linear springs to model aggregate 
interlock load transfer across pavement joints (Ioannides and Korovesis 1990; Kuo et al. 1995; Brill et al. 
1997). Further, dowel load transfer has also been idealized with linear springs spanning transverse joints 
(Ionnides and Korovesis 1992). To examine the accuracy and usefulness of this approach, consider the 
following example consisting of a simple two-slab (one row, two column) model with a 250 mm thick 
slab (E = 28000 MPa, = 0.20, density = 0) founded directly on a dense liquid with k = 0.03 MPa/mm. 
Two cases have been analyzed: the first considers only linear aggregate interlock load transfer, with joint 
stiffnesses ranging from 0 to 10 MPa/mm; the second considers no aggregate interlock load transfer, but 
has 11-32 mm dowels spaced at 300 mm on center with dowel-slab support moduli ranging from 0 to 
100,000 MPa. In both cases, the slab is subjected to an 80 kN axle located at the joint face. Each slab is 
meshed with 12x12x2 elements. Figure 5 shows a picture of the doweled model. 
After running multiple models with varying aggregate interlock joint stiffnesses and dowel-slab 
support moduli, joint load transfer efficiencies (LTEs) were computed and peak slab tensile stresses under 
each wheel were tabulated. Figure 6 shows the variation in peak tensile stress with LTE for both joint 
stiffness and dowel-slab support modulus. Note how the doweled model consistently predicts a lower slab 
stress for a given load transfer efficiency, except at the extreme values of 100% and 0% LTE. This can be 
attributed to the fact a dowel falls directly below each wheel (see Figure 5), providing a concentrated 
source of joint load transfer, whereas the constant aggregate interlock joint stiffness is evenly distributed 
across the joint width. Whenthere is perfect or no load transfer, the localized nature of the dowel support 
has no effect on slab stresses. 
 8
 
Figure 5: EverFE Screen Shot of Aggregate Interlock Example 
Figure 6: Variations in Slab Stress with Load Transfer Efficiency 
 
0 10 20 30 40 50 60 70 80 90 100 
0.9 
1 
1.1 
1.2 
1.3 
1.4 
1.5 
Evenly Spaced Dowels 
Linear Aggregate Interlock 
Load Transfer Efficiency (%) 
P
ea
k 
S
la
b 
S
tr
es
s 
(M
P
a)
 
 9
5.2 Nonlinear Aggregate Interlock Load Transfer 
The nonlinear aggregate interlock load transfer option allows the consideration of both the effect of 
relative vertical joint displacement and joint opening on aggregate interlock load transfer effectiveness. 
EverFE relies on a two-phase aggregate interlock model developed by Walraven (1981, 1994) to generate 
nonlinear aggregate interlock crack constitutive relations. The crack is assumed to follow the aggregate 
particle boundaries, and the aggregate particles bearing on the cement paste are taken to be at the point of 
slip. Walraven’s model also assumes that the aggregate particles are graded according to a Fuller 
distribution, and the maximum particle diameter, Dmax, and the aggregate volume fraction, pk, are model 
parameters. 
Given an ultimate strength of the cement paste, pu, a coefficient of friction between the paste and 
DJJUHJDWH�RI� ��DQG�FRPSXWLQJ�WKH�SURMHFWHG�FRQWDFW�DUHDV�EHWZHHQ�WKH�DJJUHJDWH�DQG�SDVWH�XVLQJ�WKH�
deformed geometry, the forces required for equilibrium of a single aggregate particle 
diameter/embedment combination can be computed. Using the probability of occurrence for a particular 
embedment/diameter combination derived by Walraven (1981), the likely forces on all particles are then 
summed to give the total forces acting on a crack plane for a given relative slip displacement and joint 
opening. 
Typical crack shear stress-displacement relations predicted by the model are shown in Figure 7, 
where each curve corresponds to a specific joint opening. Although only 3 curves are shown, EverFE 
internally generates 40 curves over a range of joint openings between 0 and 20 mm to give a very 
complete definition of the shear stress-displacement relations. The majority of these curves apply for joint 
openings between 0 and 2 mm, where the most rapid changes in load transfer effectiveness with joint 
opening occur. As with the linear model, shear stress is transferred only in the vertical direction. 
 
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 
0 
1 
2 
3 
4 
5 
6 
7 
8 
9 
Relative Vertical Joint Displacement (mm) 
S
he
ar
 S
tr
es
s 
(M
P
a)
 
0.02 mm Joint Opening 
0.12 mm Joint Opening 
0.98 mm Joint Opening 
 
Figure 7: Typical Nonlinear Aggregate Interlock Shear Constitutive Relations 
 
 10
The implementation of Walraven’s model in EverFE is consistent with the general formulation of a 
materially nonlinear finite-element analysis. The necessary tangent moduli are computed numerically 
from the constitutive relations, and stresses are numerically interpolated from the constitutive relations for 
a given joint opening and relative vertical displacement. It is important to note that the finite-element 
implementation accounts for the variation in joint opening with vertical position on the joint face that 
develops under loading of the pavement system. One limitation of the model, however, is that it does not 
account for the sawcut at the top of a typical contraction joint where no aggregate interlock load transfer 
would normally occur. 
The parameters necessary to define the nonlinear aggregate interlock model in EverFE are the 
maximum paste strength, pu, the paste-DJJUHJDWH�FRHIILFLHQW�RI�IULFWLRQ� ��WKH�DJJUHJDWH�YROXPH�IUDFWLRQ��
pk, and the maximum aggregate diameter, Dmax. Walraven (1994) suggested the following relationship 
between the compressive strength of a 150 mm concrete cube, fcc��DQG� pu, where both are in MPa: 
)6(0.8 ccpu f=σ 
The compressive strength of a 150 mm concrete cube can be assumed to be approximately 1.25 times the 
strength of a standard concrete cylinder,’cf (Wang and Salmon 1985). In addition, when weak aggregate 
is used that is prone to fracture, Walraven suggested proportioning pu downward by a fracture index, Cf < 
1.0. Recent research on the topic of aggregate interlock joint shear transfer (Jensen and Hansen 2003; 
Wattar et al. 2001) suggests that the basic concepts underlying Walraven’s model are sound. However, its 
accuracy in predicting pavement joint load transfer may vary with specific characteristics of aggregate 
shape and type, as well as the degree of damage at the joint. In addition, a study by Davids and Mahoney 
(1999) has shown good qualitative and reasonable quantitative agreement between existing experimental 
data and predictions of aggregate interlock load transfer efficiency using Walraven’s model. 
6. Treatment of Axle Loads and Thermal Effects 
EverFE allows the consideration of simultaneous axle loads and prestrains due to thermal or 
shrinkage effects. This section documents the methods by which these loads are included in each EverFE 
finite-element model. 
6.1 Axle Loads 
The complex axle configurations available in EverFE are simply collections of single rectangular 
wheel loads, and each wheel is treated identically by the finite-element code. A wheel load is defined by 
the (x,y) location of its geometric center, the length L and width W of the tire contact area, and the 
magnitude of the wheel load P. The load is assumed to produce a constant pressure over the wheel contact 
area. 
The critical issue regarding the application of the wheel loads in the finite-element model is 
determining the set of nodal forces that are equivalent to the uniformly distributed pressure generated by 
the wheel. This is challenging since the wheel load contact area is not restricted to coincide with an 
element face, and in fact can partially load several element faces. EverFE handles this by dividing each 
wheel contact area into smaller rectangular sub-areas by using a grid having nx x ny divisions along each 
edge. The ith sub-area of the wheel defined by the grid thus has an area of LW/(nxny) and sees a total force 
of pi = P/(nxny). 
The equivalent nodal force vector due to each pi is then computed by first determining the solid 
element that it contacts using the same fast geometric search procedures needed for the finite-element 
solver that is discussed in Section 7. The work-equivalent set of nodal forces is then computed as the 
product of pi and the vector of element shape functions evaluated at the point of application of pi. The 
sum of all work-equivalent nodal force vectors is the total nodal force vector applied by the entire wheel. 
This procedure is consistent with the virtual work and energy principles that form the basis of the finite-
 11
element method, using a rectangular rule to numerically integrate the tire contact pressure over its area of 
application. 
Clearly, the critical parameters in this calculation are nx and ny. To examine their effect on solution 
accuracy, consider the following example where a single slab founded on a dense liquid is subjected to a 
single 40 kN wheel load at its edge with L = W = 200 mm. The slab is 4400 mm long, 3600 mm wide, and 
254 mm thick with E = 27,60��03D�DQG� � �������7KH�VODE�LV�PHVKHG�ZLWK����[����HOHPHQWV�LQ�SODQ�DQG���
elements through the slab thickness, and is shown in Figure 8. The model was solved using values of nx = 
ny ranging between 1 and 20. When nx = ny = 1, the wheel load is treated as a single point load. Table 1 
gives the x-direction stresses at the top and bottom of the slab for different values of nx = ny, and shows 
clear convergence by nx = ny = 10. Since the equivalent nodal force calculation presented here is not 
computationally expensive, EverFE uses a fixed value of nx = ny = 20 for all simulations to ensure 
accuracy for a wide range of wheel load sizes and levelsof mesh refinement. It should be noted that other 
researchers have developed different methods of handling the problem of wheel load contact stresses 
(Hjelmstat et al. 1997); however, the method presented here is conceptually simple, easily implemented, 
and accurate. 
Difficulties can arise when one or more of the sub-area loads falls outside the slab boundary, and 
cannot be located in an element. EverFE handles this by moving the point of application of pi around a 
circle with an ever-increasing radius until pi falls within a solid element, and the calculation then proceeds 
as detailed. This may degrade the accuracy of the method, however, and wheel loads having portions of 
their contact areas outside slab boundaries should be avoided. 
 
 
 
 
Figure 8: Plan View of Mesh for Wheel Load Refinement Study 
 
 12
Table 1: Effect of nx and ny on Slab Stresses 
 
nx, ny 
Top of Slab 
Stress (MPa) 
Bottom of Slab 
Stress (MPa) 
1, 1 -2.13 1.79 
2, 2 -1.78 1.67 
3, 3 -1.81 1.68 
4, 4 -1.78 1.67 
5, 5 -1.79 1.67 
10, 10 -1.77 1.67 
15, 15 -1.77 1.67 
20, 20 -1.77 1.67 
 
6.2 Thermal and Shrinkage Effects 
Thermal and shrinkage effects are treated as general element prestrains in a manner consistent with 
fundamental finite-element theory (Zienkiewicz and Taylor 1994). EverFE allows the specification of 
linear, bi-linear, or tri-linear temperature changes throughout the slab thickness, and the prestrain is 
computed as the product of the specified temperature change and the user-specified coefficient of thermal 
expansion. Specification of shrinkage strains can be accomplished by converting the desired shrinkage 
strain into an equivalent temperature change using the coefficient of thermal expansion. The element 
prestrains are converted to an equivalent nodal force vector by the usual element integration, and are 
subtracted from the total strain during the calculation of internal stresses. 
An important detail that must be emphasized is that the 20-noded brick elements used by EverFE are 
capable of capturing an essentially linear variation in strain over their volume. Hence, when a bi-linear 
thermal gradient is specified, the finite-element mesh must have an even number of elements through its 
thickness to ensure accurate prediction of stresses and displacements. Similarly, if a tri-linear thermal 
gradient is specified, there must be 3,6,9, etc. elements through the slab thickness. These mesh restrictions 
are automatically enforced by the EverFE. 
7. Finite Element Solution Strategies 
EverFE was originally developed for use on desktop computers, and as a result uses solution 
strategies designed to minimize computational time and memory usage in its complex 3D finite element 
analyses. At the core of the solver for both linear and nonlinear problems is the solution of a system of 
symmetric, positive-definite linear equations. To accomplish this, EverFE relies on an iterative multigrid-
preconditioned conjugate gradient (MG-PCG) solver developed specifically for use on 3D finite element 
models that incorporate the multiple element types and the phenomena detailed in this manual. This 
section overviews EverFE’s global nonlinear solution strategy, and provides information on the core MG-
PCG solver. For more details, see Davids and Turkiyyah (1999). 
7.1 Global Nonlinear Solution Strategy 
EverFE finite-element models can be either linear or nonlinear. Nonlinearity is due to either material 
properties (use of a tensionless dense liquid foundation, dowel looseness, the nonlinear aggregate 
interlock model, or non-zero slab-base shear transfer when the slab-base interface is unbonded) or contact 
conditions arising from an unbonded slab-base interface. As noted in Section 4, dowel looseness is 
actually a contact nonlinearity that is more conveniently treated as a material nonlinearity. 
The nonlinear solution strategy used by EverFE is essentially a full Newton method with the updating 
of contact constraints performed at each Newton iteration. The procedure is given in Algorithm 1, where 
Kk is the tangent system stiffness matrix at the kth iteration; dUk is the update to the system displacement 
vector, U; F is the vector of applied forces; and r is the residual vector of unbalanced forces. A solution is 
 13
achieved when r is sufficiently close to zero. If the model is linear, Kk is constant, and only a single 
iteration is required. 
r = F 
while ||r||/||F|| > 10-04 
update contact constraints 
update Kk 
solve KkdUk = r 
Uk = Uk-1 + dUk 
update r 
end 
Algorithm 1: Nonlinear Solution Strategy 
 
It is worth noting that for models with a base layer, the solution of KkdUk = r for dUk involves a sub-
iteration using a technique called Uzawa’s method to satisfy the contact constraints that is not detailed 
here. This is necessary to avoid ill conditioning of Kk and maintain the efficiency of the MG-PCG solver. 
For nonlinear problems, the contribution of the 20-noded solid elements used to discretize the linearly 
elastic slabs and base to Kk is not updated. This saves significant computational time, since much of the 
work involved in forming Kk arises from the numerical integration of these large element stiffness 
matrices. 
 
7.2 Overview of Multigrid Solver 
When solving a large, 3D finite element problem – either linear or nonlinear – solution of the system 
stiffness equation 
)7(rKU = 
requires the bulk of the computational resources. For a nonlinear problem, K and U in Equation 7 
correspond to Kk and dUk in Algorithm 1. In most codes, direct solution methods such as LU factorization 
are employed to solve Equation 7, where K is factored into upper and lower triangular matrices, and the 
solution vector U is computed through forward and back substitution. While factorization is 
straightforward and relatively insensitive to poor conditioning of K, the amount of work required to factor 
K grows at least quadratically with the number of unknowns for realistic problems (even when sparse 
matrices are utilized). An additional barrier to the use of LU factorization is the additional memory 
required to store the matrix factors, which is significantly more than that required to store K itself for 
sparse matrices. 
To overcome this problem, EverFE employs a highly efficient, multigrid-preconditioned conjugate 
gradient solver. The conjugate gradient method is a widely-used iterative technique for solving systems 
such as Equation 7 that are characterized by a symmetric, positive-definite coefficient matrix, and the 
basic algorithm is easily implemented. However, it is well known that the efficiency of the conjugate 
gradient method is highly dependent on the efficiency of the preconditioner (Saad 1996). 
EverFE relies on the multigrid method to precondition the conjugate gradient method. Multigrid 
methods are themselves iterative techniques for the solution of discretizations of partial differential 
equations that rely on multiple discretizations of the same domain (Brandt 1977). Denoting the current 
error in the solution vector by e, the vector of residual nodal forces as r = F – KU, and the exact 
(unknown) solution as U*, we can write: 
 14
)9(
)8(*
rKe
UUe
=
−=
 
A small number of Gauss-Seidel iterations are performed for the finer meshes to remove high-
frequency error components, and the low frequency error components are approximated via a direct – and 
inexpensive – solution on the coarsest mesh (Brandt 1977). Figure 9 presents the algorithm and a 
conceptual overview for a two-mesh sequence. EverFE relies on a V-cycle multigrid scheme where r is 
sequentially restricted from the finest to the coarsest mesh with symmetric Gauss-Seidel smoothing 
performed at each step. Following the solution on the coarsest grid, the approximated error vector is 
sequentially interpolated and smoothed from the coarsest to the finest mesh. EverFE’s use of a single 
multigrid V-cycle to precondition a conjugate gradient iteration takesadvantage of the symmetric positive 
definiteness of the system stiffness equations. 
One of the primary difficulties in implementing a multigrid scheme for un-nested mesh sequences of 
spatially inhomogeneous domains such as those found in layered foundations is defining appropriate 
restriction and interpolation operators. Restriction can be viewed as computing a force vector on a coarse 
mesh that is statically equivalent to the known force vector on a finer mesh. This process is often 
expressed in matrix form as follows, where rc denotes a coarse mesh residual force vector, rf denotes the 
known fine mesh force vector, and R is the restriction operator: 
)10(fc Rrr = 
Similarly, interpolation is defined as the process of approximating the fine mesh error in the displacement 
vector, ef, from a known coarse mesh error, ec using the interpolation operator, T: 
)11(cf Tee = 
 
The multigrid implementation used by EverFE relies on the element shape functions to define R and 
T, which has been shown to be advantageous (Davids and Turkiyyah 1999; Fish et al. 1996). This allows 
the restriction and interpolation operations to be performed on a node-by-node basis, provided that for 
each fine mesh node, the coarse mesh element it lies within and the corresponding coarse mesh element 
coordinates are known. Efficiently establishing this information is critical, and is achieved using 
geometric search procedures detailed by Davids and Turkiyyah (1999). This general approach also 
permits the easy meshing of solid and bending elements within the same model, since all calculations are 
cf Tee =fc Rrr =
( )
( )f
cf
c
fc
ffff
f
ff
smooth
smooth
subroutine
e
Tee
rKe
Rrr
eKrr
e
er
=
=
=
−=
−1
),(MultiGrid
restrict
smooth
interpolate
Figure 9: Multigrid Concept
cf Tee =fc Rrr =
( )
( )f
cf
c
fc
ffff
f
ff
smooth
smooth
subroutine
e
Tee
rKe
Rrr
eKrr
e
er
=
=
=
−=
−1
),(MultiGrid
restrict
smooth
interpolate
Figure 9: Multigrid Concept
 15
performed at the element level. This is critical, since the dowels are explicitly modeled as flexural 
elements as discussed in Section 4. 
 
7.2 Demonstration of Solver Efficiency 
To illustrate the efficiency of EverFE’s solver, consider the solution of a model of a single 4600 mm 
long by 3600 mm wide by 250 mm thick slab resting directly on a dense liquid and subjected to a single 
axle load. Since the model is linear, the solution of the system stiffness equations is performed only once. 
Models with increasing numbers of elements and a constant maximum element aspect ratio of 4.6 were 
generated and solved using both a sparse direct solver employing LU factorization and using EverFE’s 
MG-PCG solver. The solutions were generated on a Dell Optiplex with a 2.8MHz Pentium IV processor 
and an 800 MHz front side bus. All solutions were achieved without exceeding the machine’s core RAM 
of 1 GB. It must be noted that accounting for the symmetry of the system stiffness matrix could have 
increased the efficiency of the LU factorization; however, as shown in Table 2, the results are sill 
dramatic. 
Table 2: Solver Efficiency Comparison 
 
MG-PCG LU Factorization 
Mesh 
(nx x ny x nz) 
Degrees-
of-
freedom 
Time 
(seconds) 
Memory 
(MB) 
Time 
(seconds) 
Memory 
(MB) 
4 x 4 x 1 465 1 10 1 - 
8 x 8 x 2 2,511 2 15 1 - 
12 x 12 x 3 7,293 3 28 6 75 
16 x 16 x 4 15,963 8 53 28 230 
20 x 20 x 5 29,673 15 91 123 607 
24 x 24 x 6 49,575 23 150 * * 
28 x 28 x 7 76,821 39 228 * * 
*Solution could not be achieved without exceeding core RAM of 1 GB 
The results clearly illustrate the relative efficiency of EverFE’s MG-PCG solver, especially for 
medium and large-sized problems. The times reported are totals, including time required to read all input 
data and write output stresses and displacements. Also of significance is nearly linear increase in solution 
time and memory usage with the increase in number of unknowns (degrees-of-freedom) when the MG-
PCG solver is used. It must be noted that this linear increase in computational time cannot be expected for 
all models, especially nonlinear models with either dowel looseness or slab-base contact conditions. 
Finally, while the larger meshes used here are not typical for an EverFE model in that they have a large 
number of elements through the slab thickness, their overall size is not unusually large. For example, the 
sample project “nine_slab_base” that is installed with EverFE has 55,398 degrees-of-freedom, which is 
larger than the single-slab mesh with 24 x 24 x 6 elements. 
One disadvantage of the MG-PCG solver used by EverFE – and a disadvantage with all iterative 
solvers – is its sensitivity to ill conditioning of K, which can arise from both the use of high values for 
material stiffnesses and large element aspect ratios. For example, if a large value is specified for the 
stiffness of the slab-base interface (see Section 3) or the dowel support moduli (see Section 4) the 
efficiency of the solver may suffer. Keeping maximum element aspect ratios less than the suggested limit 
of 5.0 easily controls the sensitivity to large element aspect ratios. If elements with aspect ratios greater 
than 5.0 are used in less critical regions of the model, such as shoulders, solution time may increase. 
8. Program Architecture and File Structure 
EverFE consists of four separate programs: the finite-element solver, which is written in ANSI 
standard, object-oriented C++; the meshing software, which is also written in C++; the C++ program that 
 16
generates nonlinear aggregate interlock constitutive relationships; and the user interface, which is written 
in the scripting language Tcl/Tk. 
Because of computational requirements, it was necessary to develop the finite-element, meshing, and 
aggregate interlock model definition code using a low-level compiled language such as Fortran or C/C++. 
Ultimately, C++ was chosen for two reasons: its flexibility (and thus the ease with which it allows the 
implementation of the complex solution algorithms and specialized elements detailed in this manual), and 
its effective dynamic memory allocation capabilities, which are crucial when solving large problems on 
desktop computers. The C++ code is extensive, and details of its architecture are not presented here; for 
more information, see Davids (1998). 
The user interface is written in Tcl/Tk version 8.3, a freely available scripting language that runs on 
multiple Unix, Windows and Macintosh platforms. Many of the higher-level interface widgets (dialog 
boxes, checkboxes, menus, etc.) used by EverFE were taken from the additional Tcl/Tk package Tix 
version 8.2, which is also freely available. Visualization of stresses, displacements, and dowel shears and 
moments is accomplished with the Visualization Toolkit (vtk), freely available visualization software that 
can be used with both Tcl/Tk and C++ (in EverFE it is called directly from Tcl/Tk). All of the EverFE 
source code written in Tcl/Tk is installed with EverFE and used in uncompiled form, and these files must 
never be modified. 
The remainder of this section provides details on directory-file structure and the interaction of 
EverFE’s constituent programs. 
8.1 Directory-File Structure and Program Interaction 
EverFE, like most software, is installed in a user-specified directory (the default location is 
C:\Program Files\EverFE2.23). Figure 10 shows EverFE’s top-level directory structure. 
 
All of the Tcl/Tk scripts are located in the scripts subdirectory. The FE-solver subdirectory 
contains the files driver.exe (the executable that generates the finite-element input files) new_fe.exe 
(the main finite-element executable), and agg_int.exe (a separate executable program that generates 
and saves the nonlinear aggregate interlock constitutive model information used by the finite-element 
code). 
Basic model data is stored in the data subdirectory, which for each project contains a single file witha .prj extension and a subdirectory that is necessary to store the project definition an analysis results. 
The name of both the .prj file and the subdirectory corresponding to an EverFE project are identical to 
the user-specified name the project is saved under. Within each project subdirectory, EverFE saves a file 
called model_params.dat, an ASCII text file containing the information necessary to define a project 
(model geometry, material properties, dowel information, meshing parameters, loading information, etc.). 
Figure 10: EverFE Directory Structure
User-specified installation directory
Stores nonlinear aggregate interlock models
Stores all project definitions and results
Contains interactive help file and theory manual
Contains finite-element executables
Stores Tcl/Tk scripts
Pre-compiled Tcl/Tk/Tix/vtk libraries
Figure 10: EverFE Directory Structure
User-specified installation directory
Stores nonlinear aggregate interlock models
Stores all project definitions and results
Contains interactive help file and theory manual
Contains finite-element executables
Stores Tcl/Tk scripts
Pre-compiled Tcl/Tk/Tix/vtk libraries
 17
The Tcl/Tk scripts read this file when a project is opened, and it also serves as the sole input source for 
the program driver.exe, which generates the input files defining the finite-element mesh needed by 
new_fe.exe. Note that the multigrid-preconditioned conjugate gradient solver detailed in Section 7 
utilizes three finite-element meshes with decreasing levels of refinement, and thus driver.exe actually 
generates three separate input files each time an analysis is executed. To save disk space, these files are 
deleted after the analysis is completed. When an EverFE analysis executes successfully, additional files 
are written to the project subdirectory that contain the model stresses, displacements and dowel results. 
These output files are read directly by Tcl/Tk source code when the user enters any component of the 
visualization panel; they are stored until the model is saved without re-analysis, at which point they are 
deleted to ensure that a project never has a finite-element model definition that is inconsistent with the 
saved output. 
The agg_int directory contains input and output data for each nonlinear aggregate interlock model 
that is generated and saved. There are two subdirectories within agg_int: one for models with metric 
units and one for models with English units. The documentation directory contains two files: 
EverFE_Help.chm, the compiled interactive help file developed with Microsoft’s HtmlHelp utility that 
is accessed directly from the EverFE2.23 Help menu, and the file theory_manual.pdf. Finally, the 
tcl_bins directory contains all of the binary code and libraries necessary to run Tcl/Tk/Tix/vtk. 
When pre-defined projects are run in batch mode, the Tcl/Tk code simply puts the programs 
driver.exe and new_fe.exe in a loop, sequentially analyzing each project. 
9. References 
Brandt, A. (1977). “Multi-Level Adaptive Solutions to Boundary-Value Problems.” Mathematics of 
Computation, 31(138):333 – 390. 
Brill, D.R., Hayhoe, G.F. and Lee, X. (1997). “Three-Dimensional Finite-Element Modeling of Rigid 
Pavement Structures.” Aircraft Pavement Technology: In the Midst of Change, ASCE, pp. 151-165. 
Davids, W. and Turkiyyah, G. (1997). Development of Embedded Bending Member to Model Dowel Action. 
Journal of Structural Engineering, ASCE, 123(10):1312 – 1320. 
Davids, W.G. (1998). Modeling of Rigid Pavements: Joint Shear Transfer Mechanisms and Finite Element 
Solution Strategies. PhD Dissertation, University of Washington, Seattle, WA. 
Davids, W.G. and Mahoney, J. (1999). “Experimental Verification of Rigid Pavement Joint Load Transfer 
Modeling with EverFE.” Transportation Research Record 1684, TRB, National Research Council, 
Washington, D.C., pp. 81-89. 
Davids, W.G. and Turkiyyah, G.M. (1999). “Multigrid Preconditioner for Unstructured Nonlinear 3D FE 
Models.” Journal of Engineering Mechanics, ASCE, 125(2):186-196. 
Davids, W.G. (2000). “Effect of Dowel Looseness on Response of Jointed Concrete Pavements.” Journal of 
Transportation Engineering, ASCE, 126(1):50-57. 
Davids, W.G., Wang, Z.M., Turkiyyah, G., Mahoney, J. and Bush, D. (2003). “3D Finite Element Analysis of 
Jointed Plain Concrete Pavement with EverFE2.2.” Transportation Research Record, TRB, National 
Research Council, Washington, D.C. (in press). 
Dei Poli, S., Di Prisco and Gambarova, P.G. (1992). “Shear Response, Deformations, and Subgrade Stiffness 
of a Dowel Bar Embedded in Concrete.” ACI Structural Journal, 89(6):665-675. 
Fish, J., Pan, L., Belsky, V. and Gomaa, S. (1996) “Unstructured Multigrid Method for Shells.” International 
Journal for Numerical Methods in Engineering, 39:1181 – 1197. 
Hjelmstat, K.D., Kim, J. and Zuo, K.H. (1997). “Finite Element Procedures for Three-Dimensional Pavement 
Analysis.” Aircraft/Pavement Technology: In the Midst of Change, pp. 125-137, ASCE. 
 18
Ioannides, A.M. and Korovesis, G.T. (1990). “Aggregate Interlock: A Pure Shear Load Transfer Mechanism.” 
Transportation Research Record 1286, TRB, National Research Council, Washington, D.C., 2001, pp. 14 
– 24. 
Ioannides, A.M. and Korovesis, G.T. (1992). “Analysis and Design of Doweled Slab-on-Grade Pavement 
Systems.” Journal of Transportation Engineering, 118(6):745-768. 
Jensen, E.A. and Hansen, W. (2003). “A New Model for Predicting Aggregate Interlock Shear Transfer in 
Jointed Concrete Pavements.” Proceedings of EM2003-The 16th ASCE Engineering Mechanics 
Conference, Seattle, WA, July 16-18, (CD-ROM). 
Kuo, C., Hall, K. and Darter, M.. “Three-Dimensional Finite Element Model for Analysis of Concrete 
Pavement Support.” Transportation Research Record 1505, TRB, National Research Council, 
Washington, D.C., 1996, pp. 119 – 127. 
Rasmussen, R.O. and Rozycki, D.K. (2001). “Characterization and Modeling of Axial Slab-Support 
Restraint.” Transportation Research Record 1778, TRB, National Research Council, Washington, D.C., 
pp. 26 – 32. 
Saad, Y. (1996). Iterative Methods for Sparse Linear Systems. PWS Publishing Co., Boston, MA. 
Walraven, J.C. (1981). “Fundamental Analysis of Aggregate Interlock.” Journal of the Structural Division, 
ASCE, 107(ST11):2245 – 2270. 
Walraven, J.C. (1994). “Rough Cracks Subjected to Earthquake Loading.” Journal of Structural Engineering, 
120(5):1510 – 1524. 
Wang, C-K. and Salmon , C.G. (1985). Reinforced Concrete Design (4th Ed.). Harper and Row, New York. 
Wattar, S.W., Hawkins, N.M. and Barenberg, E.J. (2001). “Aggregate Interlock Behavior of Large Crack 
Width Concrete Joints in PCC Airport Pavements.” Technical Report, Department of Civil and 
Environmental Engineering, University of Illinois at Urbana-Champaign. 
Zhang, J. and Li, V.C. (2001). “Influence of Supporting Base Characteristics on Shrinkage-Induced 
Stresses in Concrete Pavements.” Journal of Transportation Engineering, ASCE, 127(6);455 – 642. 
Zienkiewicz, O.C. and Taylor, R.L. (1994). The Finite Element Method, Volume 1 (4th Ed.). McGraw Hill 
Book Company, London.