U.S. patent application number 14/624696 was filed with the patent office on 2015-10-01 for asphalt composition.
The applicant listed for this patent is NUtech Ventures. Invention is credited to David Allen, Yong-Rak Kim, Roberto Soares, Flavio Souza.
Application Number | 20150275440 14/624696 |
Document ID | / |
Family ID | 44066823 |
Filed Date | 2015-10-01 |
United States Patent
Application |
20150275440 |
Kind Code |
A1 |
Allen; David ; et
al. |
October 1, 2015 |
ASPHALT COMPOSITION
Abstract
A computational method is provided for predicting roadway
failure due to degradation of the roadway over time as a function
of the input loads, the roadway geometry, the material properties
of the constitutents in the asphaltic pavement, the shape,
distribution, orientation and volume fractions of the constituents,
and environmental conditions. The unique and new feature of the
method is that it employs several physically based predictive
methodologies simultaneously.
Inventors: |
Allen; David; (Lincoln,
NE) ; Souza; Flavio; (Lincoln, NE) ; Kim;
Yong-Rak; (Lincoln, NE) ; Soares; Roberto;
(Lincoln, NE) |
|
Applicant: |
Name |
City |
State |
Country |
Type |
NUtech Ventures |
Lincoln |
NE |
US |
|
|
Family ID: |
44066823 |
Appl. No.: |
14/624696 |
Filed: |
February 18, 2015 |
Related U.S. Patent Documents
|
|
|
|
|
|
Application
Number |
Filing Date |
Patent Number |
|
|
13512183 |
Sep 24, 2012 |
9091027 |
|
|
PCT/US2009/066071 |
Nov 30, 2009 |
|
|
|
14624696 |
|
|
|
|
Current U.S.
Class: |
703/12 |
Current CPC
Class: |
G06F 2113/26 20200101;
G06F 30/23 20200101; G06F 2119/04 20200101; E01C 1/002
20130101 |
International
Class: |
E01C 1/00 20060101
E01C001/00; G06F 17/50 20060101 G06F017/50 |
Claims
1. (canceled)
2. A computer-based method of predicting roadway failure under
loads applied to an asphalt surface of a roadway by vehicles,
comprising: inputting into a graphical user interface one or more
values relating to vehicular forces applied to an asphalt surface
of a roadway; determining macroscale stresses resulting from the
vehicular forces throughout a volume below the asphalt surface; and
determining microscale fracture process due to the macroscale
stresses; wherein the microscale fracture process is predictive of
roadway failure.
3. The method of claim 2, wherein the volume comprises a plurality
of layers, and determining macroscale stresses comprises:
determining a global scale finite element mesh having a boundary
coincident with a boundary of the roadway; determining boundary
conditions at the boundary of the global scale mesh from parameters
associated with the geometry of the roadway and the forces applied
by vehicles; and determining values of the stresses over the global
scale mesh by solving the Initial Boundary Value Problem (IBVP)
subject to the boundary conditions.
4. The method of claim 3, wherein determining values of the
stresses comprises solving a global scale Initial Boundary Value
Problem (IBVP) for the global mesh, obtaining stresses and
displacement fields as outcomes.
5. The method of claim 4, wherein the global scale IBVP comprises
conservation of linear momentum, conservation of angular momentum,
strain-displacement relationship, and constitutive equations for
the global scale object, in this case, the roadway.
6. The method of claim 5, wherein the constitutive equations depend
upon the current values of strains or the history of strains.
7. The method of claim 5, wherein the equations in the global scale
IBVP are solved using the finite element method.
8. The method of claim 3, wherein the microscale fracture process
is expressed in terms of a set of microscale fracture parameters,
the set of microscale fracture parameters comprising a damage
parameter defined in the cohesive zone ahead of crack tips.
9. The method of claim 8, wherein determining microscale fracture
process comprises determining: a set of length scales having a
hierarchical order, wherein the largest length scale is associated
with the global scale mesh; a plurality of local meshes associated
to integration points in the global mesh, wherein a local mesh has
a boundary and is associated with a length scale different from the
length scale associated with the global mesh; boundary conditions
of the local scale meshes from values of the strains or stresses at
locations on the global scale mesh, in particular, the location of
integration points in the global scale mesh; and values of the
stress field, the displacement field, and the damage parameter over
the local scale meshes by solving local scale IBVPs for each local
mesh using boundary conditions associated with the local mesh;
wherein the stress field on a length scale associated with a local
mesh determines the stresses at locations on the length scale
associated with the global mesh, and the values of the damage
parameter over the local scale mesh describes the fracture process
within the local scale mesh.
10. The method of claim 9, wherein the local scale IBVP comprises
conservation of linear momentum, conservation of angular momentum,
the strain-displacement relationship, constitutive equations, and a
fracture criterion.
11. The method of claim 10, wherein the stresses at locations on
the length scale associated with the global mesh are related to the
stress field over the local scale mesh by homogenization
principles, wherein homogenization principles comprise an averaging
of the stress field on the length scale associated with the local
scale mesh.
12. The method of claim 11, wherein determining microscale fracture
process further comprises determining: a plurality of cohesive
zones ahead of crack tips, wherein the length scale associated with
a cohesive zone is smaller than the length scale associated with a
local scale mesh; boundary conditions on the cohesive zones,
wherein a cohesive zone has a boundary and the boundary conditions
of a cohesive zone are determined from the traction forces observed
at locations on the local scale meshes over the surfaces of the
cohesive zones; and values of the traction forces, the displacement
field, and a scalar damage parameter on the local scale meshes by
solving a microscale IBVP for a cohesive zone using the boundary
conditions of the cohesive zone; wherein the values of the traction
forces and the displacement field on a length scale associated with
a cohesive zone determine the values of the traction forces and the
displacement field on a length scale associated with a local scale
mesh containing the cohesive zones, and the values of the scalar
damage parameter describes the fracture process within the cohesive
zone.
13. The method of claim 12, wherein the microscale IBVP comprises a
traction-displacement relation.
14. The method of claim 13, wherein the values of the traction
forces on the length scale associated with the local scale mesh are
related to the values of the traction forces on a length scale
associated with a cohesive zone by homogenization principles,
wherein the homogenization principles comprise an averaging of the
values of the traction forces on the length scale associated with
the cohesive zone.
15. The method of claim 12, wherein the scalar damage parameter is
a sum of relative areas of undamaged fibrils throughout a cohesive
zone.
16. The method of claim 3, wherein the parameters are chosen from:
aggregate volume fraction, aggregate shape, aggregate angularity,
aggregate stiffness, aggregate toughness, and asphalt binder
viscoelastic properties.
17. The method of claim 2, wherein the roadway surface material
comprises an asphalt concrete material.
18. The method of claim 17, wherein the roadway volume comprises a
hot-mix asphalt (HMA) layer, a base layer, a subbase layer, and a
subgrade layer, where the subgrade layer is the existing soil, the
subbase layer is built on top of the subgrade layer, the base layer
is built on top of the subbase layer and the HMA layer is built on
top of the base layer.
19. The method of claim 18, wherein the HMA layer is modeled as
viscoelastic material.
20. The method of claim 19, wherein the base layer, the subbase
layer and the subgrade layer are modeled as elastic materials.
21. A computer-based method of determining a desired composition of
a roadway material, comprising: providing initial values of
parameters in a first set of parameters representing the
com-position of the roadway material, the first set of parameters
having one or more values associated with it; predicting the
failure of the roadway material as recited in claim 2; determining
a lifetime of the roadway material from the failure prediction;
adjusting the values of the first set of parameters in the set of
parameters to form a second set of parameters; and determining a
lifetime of the roadway material from the second set of parameters
so as to improve the lifetime of the roadway material sufficiently.
Description
1 TECHNICAL FIELD
[0001] This invention is a computational method for predicting the
life of asphaltic roadways as a function of loading distribution in
both time and space, geometry of the roadway, material properties
of the constituents, and environmental conditions.
2 BACKGROUND
[0002] The worldwide annual cost of roadways constructed with
asphalt has been estimated at $500B, thus making asphalt arguably
the most used material on the planet. Despite this fact, no robust
methods are in place in any country that are capable of accurately
predicting the life of a roadway constructed from asphaltic
pavement as a function of loads, geometry of the roadway, material
properties of the pavement, and environmental conditions. Such a
method would have the potential to reduce the cost of paving with
asphalt by perhaps tens of billions of dollars annually with
improved performance and attendant safety.
3 SUMMARY
[0003] The disclosure provides a computational method for
predicting roadway failure due to degradation of the roadway over
time as a function of the input loads, the roadway geometry, the
material properties of the constitutents in the asphaltic pavement,
the shape, distribution, orientation and volume fractions of the
constituents, and environmental conditions. The unique and new
feature of the method is that it employs all of the following
physically based predictive methodologies simultaneously: [0004] 1.
Viscoelastic properties of the asphaltic matrix [0005] 2. Nonlinear
viscoelastic crack extension at both the local (aggregate) and
global (road-way) scales [0006] 3. Temperature dependence of
constituent material properties [0007] 4. A multiscale algorithm
that is capable of predicting the interaction between cracks on the
local and global scales [0008] 5. Ability to account for arbitrary
spatial and time dependent distributions of tire loadings [0009] 6.
Ability to account for any geometry of the roadway [0010] 7. The
multiscale algorithm results in the ability to predict roadway as a
function of asphaltic pavement mix properties, such as aggregate
volume fraction, angularity, size and spatial distribution, and
orientation at the constituent scale
[0011] Clearly, the ability to account for all of the phenomena
described above will result in a method that is capable of being
utilized as an accurate tool for designing improved asphaltic
roadways
4 CURRENT STATE OF THE ART
[0012] Current state-of-the-art methods for predicting performance
behavior of asphalt concrete mixtures and pavement structures are
generally based on two mechanics approaches: continuum damage
mechanics approach and micromechanical approach.
[0013] The continuum damage models have been widely used and
demonstrated great success in predicting structural degradation of
asphalt mixtures and pavement structures due to damage. By
representing the damaged body as a homogeneous continuum, the state
of damage is quantified by damage parameters so-called internal
state variables within the context of the thermodynamics of
irreversible processes. The damage parameters are selected somewhat
arbitrarily and determined usually by fitting damage evolution
characteristics to the available laboratory performance test
results. The continuum damage models thus provide a simplicity and
efficiency to the modeling of macroscopic mechanical behavior of
materials and structures with distributed damage. These models,
therefore, have shown limitations as they do not account for the
explicit microstructure distribution in the mixtures and
hypothesize that all damage in the mixtures can be averaged and
represented by the phenomenological internal state variables. This
infers that the approach is incapable of predicting detailed energy
dissipation mechanisms and damage evolution characteristics based
on the underlying physics of the damage which is mostly controlled
by the interactions of mixture constituents. Furthermore, new tests
must be performed whenever mixture constituents or applied loading
conditions are changed. The time and costs associated with this
technique limit its usability to select materials and control
design variables such that the microstructure distribution is
optimized and damage is minimized.
[0014] In the micromechanics approach, the aforementioned
shortcomings of continuum damage mechanics approaches can be
overcome to a certain extent, since the micromechanics modeling
approach directly accounts for the effect of mixture heterogeneity
by dealing with mixture constituents separately. Furthermore, with
the aid of computational techniques such as the finite element
method and the discrete element method, the approach has
successfully simulated geometric heterogeneity and predicted
stiffness reduction of mixtures due to damage and fracture of
individual constituents. However, most of micromechanical models
currently available are still limited to provide realistic
information about the macroscopic constitutive framework because of
the intrinsic complexity of the asphalt mixture microstructure and
the microscopic mechanisms associated with interactions among the
heterogeneities. Furthermore, because of technical challenges, the
majority of micromechanics models did not take into account either
material inelasticity (e.g. viscoelasticity) or rate-dependent
nonlinear damage growth (including energy dissipation due to
micro/macro cracking) that are commonly observed phenomena in
asphaltic composites. Another big challenge of the micromechanics
modeling is computational effort. Considering thousands of
irregularly shaped, randomly oriented aggregate particles, along
with thousands of potential cracking sites, a highly refined mesh
with a huge degree of freedom is necessary to model global
pavement-scale, which results in an impractical computational
cost.
[0015] In summary, the currently used methods for predicting life
of asphaltic roadways all suffer from one or more of the following
shortcomings: [0016] 1. They are semi-empirical [0017] 2. They make
physical approximations that are not supported by observations
[0018] 3. They do not accurately account for energy dissipation via
viscoelasticity and/or plasticity [0019] 4. They do not account for
energy dissipation due to crack growth at both the global and the
local scales [0020] 5. They do not account for interactions between
global and local cracks [0021] 6. They are not capable of
accounting for input design variables at the local scale (such as
fines, constituent volume fractions, aggregate angularity, particle
orientation, size, and spatial distribution, etc.) [0022] 5. They
cannot account for the effects of variations in tire loading in
both space and time
[0023] While not all of the currently available methodologies
suffer from all of the above short-comings, they all suffer from at
least one or more of them. Therefore, there is a pressing need for
reliable models that can account for the microstructural details
but with only a fair amount of computational effort.
[0024] The multiscale modeling approach disclosed herein takes
advantage of the strengths of both micromechanics and continuum
damage mechanics, that is, the ability of the micromechanical model
to describe mixture behavior using local scale component material
properties and that of the continuum damage model to describe the
global scale stress-strain behavior of the pavement materials. The
currently proposed method thus overcomes all of these challenges,
which renders it a robust tool for modeling asphaltic pavement
life.
5 DESCRIPTION OF THE FIGURES
[0025] FIG. 1 contains a schematic representation of a three scale
Initial Boundary Value Problem (IBVP).
[0026] FIG. 2 contains a schematic representation of a cohesive
zone.
[0027] FIG. 3 contains an illustration of multiple scales of a
roadway.
[0028] FIG. 4 contains a two-dimensional representation (and FE
mesh) of a roadway.
[0029] FIG. 5 contains a method flowchart.
[0030] FIG. 6 contains a two-dimensional finite element mesh of the
a close-up view of a roadway (local scale).
[0031] FIG. 7 contains a multiscale simulation of an asphalt
pavement (FIG. 7A) where excessive microcracking occurred (FIG.
7B), thus indicating failure of the pavement.
6 DESCRIPTION OF THE METHOD
[0032] The method proceeds from the fundamental assumption that
failure of asphaltic roadways occurs due to the accumulation of a
series of interconnected energy dissipating events. These include
the dissipation of energy in the asphaltic mastic due to material
viscoelasticity (such as polymer chain sliding), the accumulation
of local cracks between the mastic and the aggregate, the
development and propagation of global cracks that can be seen with
the naked eye on the surface of the roadway, and chemical changes
due to environment and load accumulation over time. The method uses
a finite element formulation to account for spatial variations in
the state variables (stresses, strains, and displacements), and the
evolution of cracks, both local scale and global scale. The
evolution of cracking is accounted for in the formulation via a
micromechanically based nonlinear viscoelastic cohesive zone model
previously developed by one of the inventors and collaborators.
Through the deployment of this algorithm, the method has the
capability to predict the evolution of the thermomechanical
response of the roadway to the input conditions over time. Using
existing failure criteria, this evolving picture of the roadway is
then utilized to predict when the roadway will fail due to typical
failure mechanisms such as: rutting, reflective thermal cracking,
rutting/surface cracking; and environmental degradation.
[0033] In order to account for the fact that cracks develop on
widely differing length scales in asphaltic pavement, the method
utilizes a totally new and ingenious approach developed by the
inventors called herein multiscaling. While this term is often
utilized by other researchers, the context of this term in the
current application is different than that utilized by other
researchers. This approach essentially consists of utilizing global
finite element meshes wherever global scale cracks are expected,
and local finite element meshes wherever local scale cracks are
predicted. These global and local analyses are two way coupled via
the deployment of mathematically rigorous homogenization theorems
also previously developed by the inventors and coworkers. This
approach is described conceptually in FIG. 1, wherein three
simultaneous length scales are shown.
[0034] The smallest length scale (called the micro-scale) is
associated with the cohesive zones wherein cracks initiate and
grow. A depiction of this length scale is shown in FIG. 2. This
portion of the method has been developed by the inventors (and
coworkers), and is entirely analytical.
[0035] The remaining two length scales in the figure denote the
length scale of the typical aggregate dimension (local scale), and
the roadway dimension (global scale). These two length scales are
modeled concomitantly by coupled finite element codes that simulate
the response of the roadway on the scales shown in FIG. 3, and the
codes have the capability to simulate the evolution of cracks on
both length scales using the cohesive zone model briefly described
above.
[0036] When this algorithm is applied to a typical roadway
simulation in two dimensions, the general representation shown in
FIG. 3, results in the depiction shown in FIG. 4.
[0037] A simplified flowchart of the multiscale finite element code
is shown in FIG. 5, wherein it can be seen that the global scale
and local scale computations are two-way coupled, thus producing an
algorithm in which local scale objects (such as aggregate) affect
the global scale response, and vice versa.
[0038] The resulting model has been cast into a user friendly
format in which cyclic loading with an arbitrarily imposed surface
traction distribution (due to tires or other imposed loadings) can
be applied as an input. Furthermore, the algorithm can be used to
construct a geometric model of the roadway, as shown in FIG. 3.
Fracture properties of the cohesive zones are quite easily
implemented.
[0039] A user can input one or more parameter values relating to
vehicular forces applied to an asphalt surface of a roadway into
one or more graphical user interfaces (GUIs). Such GUIs can be part
of a front end for the integrated software package. A front end
feeds user inputs from the GUI into a back end for the integrated
software package. The purpose of the back end would be to perform
computations of, for instance, roadway lifetime, corresponding to
the roadway model using the parameters.
[0040] FIG. 4 shows a typical global scale mesh constructed with
the code for purposes of performing a two dimensional simulation of
a typical roadway. FIG. 6 shows a close-up view of the local scale
mesh, with mastic and aggregate depicted in different colors. Note
that these local scale meshes are repeated at all locations
throughout the asphaltic layer of the roadway, thus producing
hundreds of local scale meshes embedded within the single global
scale mesh. The resulting algorithm is highly computationally
intensive, and is most efficiently utilized on a massively parallel
mainframe computer due to time constraints resulting from the
time-stepping nature of the algorithm.
[0041] The results of the time-stepping scheme are that all of the
resulting state variables (stresses, strains, displacements, and
evolving cracks) are modeled as continuous functions of time and
space. In particular, the displacement history at every point in
the roadway is modeled over the life of the roadway due to the
cyclic load history imposed by vehicular traffic. When exported to
the proper visualization tools, these output state variables can be
used to plot the time history of the entire roadway as a function
of all of the input design variables. For example, a road designer
could make successive predictions of the life of a roadway
subjected to specific traffic conditions as a function of the
aggregate volume fraction or the average aggregate particle size.
Such capabilities do not currently exist.
[0042] Finally, using currently available failure criteria for
different failure modes, film clips of the evolution of the roadway
are used within a post-processing algorithm to predict when the
roadway will fail. FIG. 7 presents an example of a simulation where
excessive microcracking has occurred thus indicating that the
pavement is close to failure.
7 DETAILED DESCRIPTION OF THE METHOD
[0043] FIG. 1 presents a schematic depiction of a three scale
Initial Boundary Value Problem (IBVP) of a roadway. In this figure,
superscripts are used to denote the scale index whereas subscripts
are used to denote the conventional rectilinear coordinate system.
The problem of interest (roadway) is referred as scale 0 and
considered to be statistically homogeneous. The local scale object
(Representative Volume Element-RVE), referred as scale 1, on the
other hand, is heterogeneous. Various sources of heterogeneity can
be considered such as aggregates, voids, cracks, additives and
others. The cohesive zone ahead of crack tips is herein considered
to be small enough so it can be treated as a third scale (scale 2).
The third scale has however been analytically solved using
homogenization techniques and considering some simplifications in
order to produce a homogenized traction-displacement relationship
for the cohesive zone. V.sup..mu., .differential.V.sub.E.sup..mu.
and .differential.V.sub.I.sup..mu. denote the volume, the external
and the internal boundary surfaces of the object at scale .mu.,
respectively. The length scale associated with the object and with
the cracks at scale .mu. are respectively denoted by l.sup..mu. and
l.sub.e.sup..mu..
7.1 Multiscale Approach
[0044] The method presented herein is described in terms of
infinitesimal kinematics in order to simplify the constitutive
material modeling, but it can be extended to finite deformations
with careful attention to the description of the cohesive zone
deformation and constitutive law.
7.1.1 Micro Scale--Scale 2
[0045] Beginning with the smallest scale considered in the model, a
micromechanical viscoelastic cohesive zone model is used to
describe the fracture of behavior of asphalt materials. The
micromechanical viscoelastic cohesive zone model is postulated to
be represented by a fibrillated zone that is small compared to the
total cohesive zone area, as shown in FIG. 2. This model is
inherently two scale in nature in that it utilizes the solution to
a microscale continuum mechanics problem, together with a
homogenization theorem to produce a cohesive zone model on the next
larger length scale. Furthermore, this model has also been shown to
be consistent with the fact that the cohesive zone requires a
nonstationary critical energy release rate in order for a crack to
propagate.
[0046] The Initial Boundary Value Problem (IBVP) for scale 2 is
thus given by
.sigma. ji , j 2 = 0 in V 2 ( 1 ) .sigma. ij 2 = .sigma. ji 2 in V
2 ( 2 ) ij 2 = 1 2 ( u i , j 2 + u j , i 2 ) in V 2 ( 3 ) .sigma.
ij 2 ( x p 2 , t ) = .intg. - .infin. t C ijkl 2 ( x p 2 , t -
.tau. ) .differential. kl 2 ( x p 2 , .tau. ) .differential. .tau.
.tau. in V 2 ( 4 ) ##EQU00001##
[0047] Adjoining appropriate boundary conditions, the IBVP above is
said to be well-posed. Taking a closer look into the mechanics
ahead of the crack tip, a locally averaged solution of the cohesive
zone IBVP (scale 2) can be obtained analytically, so the following
traction-displacement relationship results for the viscoelastic
cohesive zones on the next larger scale (scale 1):
t i 1 ( t ) = 1 .lamda. 1 [ u i 1 ] .delta. i * 1 [ 1 - .alpha. 1 (
t ) ] { .sigma. i f 1 + .intg. t 0 t E CZ 1 ( t - .tau. )
.differential. .lamda. 1 .differential. .tau. .tau. } on
.differential. V CZ 1 ( 5 ) ##EQU00002##
where t.sub.i.sup.1 is the traction vector acting on the cohesive
zone boundary, [.XI..sub.i.sup.1] is the cohesive zone opening
displacement, .delta.*.sub.i.sup.1 are empirical material length
parameters which typically reflect a length scale intrinsic to the
cohesive zone, .sigma..sub.i.sup.f1 is the required stress level to
initiate damage, E.sub.CZ.sup.1 is the uniaxial relaxation modulus
of a single fibril in the cohesive zone, .alpha..sup.1(t) is the
internal damage parameter reflecting the area fraction of voids
with respect to the cross-sectional area of the idealized cohesive
zone, and .lamda..sup.1 is the Euclidean norm of the cohesive zone
opening displacements:
.lamda. 1 = ( [ u r 1 ] .delta. r * 1 ) 2 + ( [ u n 1 ] .delta. n *
1 ) 2 + ( [ u s 1 ] .delta. s * 1 ) 2 ( 6 ) ##EQU00003##
where n indicates normal direction and r and s indicate tangential
directions.
7.1.2 Local Scale--Scale 1
[0048] The IBVP for the local RVE is well-posed by uniform initial
boundary conditions adjoined to the following set of equations:
.sigma. ji , j 1 = 0 in V 1 ( 7 ) .sigma. ij 1 = .sigma. ji 1 in V
1 ( 8 ) ij 1 = 1 2 ( u i , j 1 + u j , i 1 ) in V 1 ( 9 ) .sigma.
ij 1 ( t ) = .OMEGA. .tau. = - .infin. .tau. = t { kl 1 ( .tau. ) }
in V 1 ( 10 ) G i 1 .gtoreq. G ci 1 .differential. .differential. t
( .differential. V I 1 ) > 0 in V 1 ( 11 ) ##EQU00004##
where equation 11 describes the fracture criterion, G.sub.i.sup.1
is the fracture energy release rate at a particular position in the
local scale, G.sub.ci.sup.1 is the critical energy release rate of
the material at that particular position and the index i refers to
the mode of fracture. This form of the fracture criterion will
herein be naturally replaced by the cohesive zone model given in
Eq. 5. It can be shown that equation 5 is mathematically equivalent
to equation 11. Even though the mathematical proof is not shown
herein, conceptually, as the damage parameter .alpha..sup.1(t)
reaches unity, the traction vector t.sub.i.sup.1(t) becomes zero
valued, meaning that a free surface has been created or,
equivalently, a crack has propagated,
.differential. .differential. t ( .differential. V I 1 ) > 0.
##EQU00005##
[0049] Furthermore, the functional
.OMEGA..sub..tau.=-.infin..sup..tau.=t describing the constitutive
behavior of each particle at the local position x.sub.i.sup.1 is
known a priori for all constituents. In the current study, this
functional can represent linear elastic or linear viscoelastic
behavior depending on the constituent. For linear viscoelastic
materials, functionals of single integral type are considered:
.sigma. ij 1 ( x p 1 , t ) = .intg. - .infin. t C ijkl 1 ( x p 1 ,
t - .tau. ) .differential. kl 1 ( x p 1 , .tau. ) .differential.
.tau. .tau. ( 12 ) ##EQU00006##
where x.sub.p.sup.1 denotes the particle under consideration.
7.1.3 Global Scale--Scale 0
[0050] The mechanical IBVP for the global object (roadway) can then
be well-posed by appropriate initial boundary conditions and a set
of equations: conservation of linear and angular momentum
(equations 13 and 14), infinitesimal strain-displacement relations
(equation 15) and constitutive equations (equation 16).
.sigma..sub.ji,j.sup.0=0 in V.sup.0 (13)
where Einstein's indicial notation has been used,
.sigma..sub.ij.sup.0 is the Cauchy stress tensor. Note that
inertial effects and body forces have been neglected since the
roadway problem will be considered quasi-static.
.sigma..sub.ij.sup.0=.sigma..sub.ji.sup.0 in V.sup.0 (14)
.epsilon..sub.ij.sup.0=1/2(u.sub.i,j.sup.0+u.sub.j,i.sup.0) in
V.sup.0 (15)
where .epsilon..sub.ij.sup.0 is the infinitesimal strain tensor
defined on the global scale and u.sub.i.sup.0 is the displacement
vector at a certain location x.sub.i.sup.0 inside the statistically
homogeneous object at the global length scale and at a certain time
t.
.sigma..sub.ij.sup.0(x.sub.m.sup.0,t)={right arrow over
(.OMEGA.)}.sub..tau.=-.infin..sup..tau.=t{.epsilon..sub.kl.sup.0(x.sub.m.-
sup.0,.tau.)} in V.sup.0 (16)
where {right arrow over (.OMEGA.)}.sub..tau.=-.infin..sup..tau.=t
is a functional describing the constitutive behavior at each
position x.sub.i.sup.0 in the object at the global scale, which may
account for damage accumulation and history dependent effects, such
as viscoelasticity, and is obtained by locally averaging the
constitutive behavior of the RVE corresponding to that particular
global scale position. 7.1.4 Two-Way Coupling between Global and
Local Scales In the most general sense, the local RVE should obey
the same set of governing equations as the global object. However,
complex inhomogeneous boundary conditions might have to be applied
at the RVE boundaries in order to satisfy that condition. Therefore
some simplifications are assumed so that the complexities may be
neglected and one may seek an approximate, but satisfactorily
accurate, solution for the problem. The first assumption is that
the physical length scales are widely separated, allowing the use
of spatially homogeneous boundary conditions (uniform boundary
tractions or linear boundary displacements) and simplifying the
homogenization relationships connecting global and local
scales:
l.sup.0>>l.sup.1>>l.sup.2 (17)
[0051] The most common boundary conditions applied to the
boundaries of the RVE are the linear displacement, constant
tractions, periodic displacement and anti-periodic tractions. Any
choice of imposed boundary conditions should however satisfy the
Hill condition which requires that the total work-rate per unit
reference volume (or the virtual work) at the global (macro) scale
be equal to the total work-rate per unit reference volume on the
local (micro) scale:
.sigma. ij 0 . ij 0 = 1 V 1 .intg. V 1 .sigma. ij 1 . ij 1 V ( 18 )
##EQU00007##
[0052] The second assumption is that the local length scale is much
larger than the length scale associated with cracks at the local
scale 1, thus mitigating the possibility of statistical
inhomogeneity at the local scale but narrowing the range of
problems for which the method can generate accurate
predictions:
l.sup.1>>l.sub.c.sup.1 (19)
[0053] For the cases where the previous assumptions are valid, one
can show that the relations connecting the global and local scales
are given by the following (mean field) homogenization
principles:
_ ij 1 = 1 V 1 .intg. V 1 ij 1 V = ij 0 + .alpha. ij 0 ( 20 ) ij 0
= 1 V 1 .intg. .differential. V E 1 1 2 ( u i 1 n j 1 + u j 1 n i 1
) S ( 21 ) .alpha. ij 0 = 1 V 1 .intg. .differential. V I 1 1 2 ( u
i 1 n j 1 + u j 1 n i 1 ) S ( 22 ) .sigma. ij 0 = 1 V 1 .intg.
.differential. V E 1 .sigma. ki 1 n k 1 x j 1 S = 1 V 1 .intg. V 1
.sigma. ij 1 V = .sigma. _ ij 1 ( 23 ) ##EQU00008##
where n.sub.i.sup.1 is the outward unit normal vector to the
external or internal boundary, .differential.V.sub.E.sup.1 or
.differential.V.sub.I.sup.1, of the RVE.
[0054] From equation 21, one can observe that the strain tensor in
the global scale is equal to the external boundary average of the
displacement vector on the local scale, which differs from the
volume average of the strain tensor by the internal boundary
average of the displacement vector, .alpha..sub.ij.sup.0. Note,
however, that in the case where there is no separation along the
internal boundaries of the RVE, no voids or cracks,
.alpha..sub.ij.sup.0 is zero, so that .alpha..sub.ij.sup.0 can be
regarded as an averaged measure of damage.
[0055] On the other hand, equation 23 allows us to conclude that
the external boundary average of the traction vector,
t.sub.i.sup.1, is equal to the volume average of the stress tensor.
This is true only if the traction vector over the internal boundary
surfaces is zero (free surfaces, cracks) or self-equilibrating
(cohesive zones).
[0056] The use of equations 21 and 23 is termed a mean field theory
of homogenization because the behavior of the global object is
determined only in terms of the mean stress and strain tensors
evaluated at the external boundary of the local RVE. It is
important to note that the assumption that .epsilon..sub.ij.sup.0
is spatially uniform along the external boundary of the RVE breaks
down in regions of high gradients, such as in the vicinity of
cracks or where the size of the RVE is comparable to the global
length scale.
[0057] For quasi-static problems, such as the life prediction of a
roadway, the computation of the homogenized constitutive tensor is
required in order to correctly calculate the tangent stiffness
matrix. In other words, one is required to compute the homogenized
functional {right arrow over
(.OMEGA.)}.sub..tau.=-.infin..sup..SIGMA.=t. If one uses an
incremental solution scheme, however, it is preferable to seek an
incremental constitutive equation of the form:
.DELTA..sigma..sub.ij.sup.0={right arrow over
(.OMEGA.)}.sup.t(.DELTA..epsilon..sub.kl.sup.0) (24)
7.2 Incremental Formulation for Time Dependence
[0058] General purpose Finite Element codes usually employ
incremental algorithms for solving IBVP's. When the IBVP is time
and/or history dependent, incrementalization is necessary to allow
integration of the governing equations, at least in an approximate
manner. This disclosure presents an incremental homogenization
scheme for viscoelastic solids containing cracks suitable for
multiscale computational codes for predicting life of asphaltic
roadways.
[0059] Most approaches proposed in the literature to compute the
homogenized tangent constitutive tensor are based on purely
numerical techniques. For example, some researchers employ the
technique of condensation of the total RVE stiffness. Others use
the approximate numerical differentiation technique, which requires
the solution of four (2D) or six (3D) IBVP's for every solution
step. Basically, numerical constitutive tests are performed by
applying unit deformation in every mode. Note that these works do
not consider cracking or viscoelasticity. The approach to be
presented in this disclosure, on the other hand, is based on
continuum mechanics homogenization principles and allows the
computation of the full homogenized anisotropic tangent
constitutive tensor of viscoelastic solids with evolving cracks by
solving the IBVP only once at each coordinate location and for each
time step, which greatly reduces the required amount of
computational time. This procedure is presented in a scientific
journal article by the first two inventors of this disclosure.
[0060] As will be shown, the introduction of a localization
quantity relating the local displacement field to the deformation
on the external boundary of the RVE is the key feature that allows
the computation of the full homogenized anisotropic tangent
constitutive tensor by solving the IBVP only once. For the problem
at hand, first assume that the stress tensor at both global and
local scales can be written in incremental form:
.sigma. ij 0 ( t + .DELTA. t ) = .sigma. ij 0 ( t ) +
.DELTA..sigma. ij 0 ( 25 ) .sigma. ij 1 ( t + .DELTA. t ) = .sigma.
ij 1 ( t ) + .DELTA..sigma. ij 1 Therefore , ( 26 ) .sigma. ij 0 (
t ) = 1 V 1 .intg. V 1 .sigma. ij 1 ( t ) V ( 27 ) .DELTA..sigma.
ij 0 = 1 V 1 .intg. V 1 .DELTA. .sigma. ij 1 V ( 28 )
##EQU00009##
[0061] A very stable and accurate incrementalization scheme has
been proposed by one of the inventors (and collaborators) which can
model viscoelastic constitutive behavior of single integral type
and for which the relaxation moduli can be represented by a Prony
series:
C ijkl 1 ( t ) = C ijkl .infin. 1 + m = 1 M ijkl 1 C ijkl m 1 - t /
.rho. ijkl m 1 ( no sum on i , j , k , l ) ( 29 ) ##EQU00010##
where C.sub.ijkl.sub..infin..sup.1 are the long term moduli,
C.sub.ijkl.sub.m.sup.1 are the moduli coefficients and
M.sub.ijkl.sup.1 is the number of Prony terms necessary to
accurately describe the material relaxation moduli, and
.rho..sub.ijkl.sub.m.sup.1 are the so-called relaxation times given
by,
.rho..sub.ijkl.sub.m.sup.1=.eta..sub.ijkl.sub.m.sup.1/C.sub.ijkl.sub.m.s-
up.1 (no sum on i, j, k, l) (30)
where .eta..sub.ijkl.sub.m.sup.1 are viscosity terms.
[0062] If one assumes that the strain rate is constant for each
time step, the following incremental constitutive equation can be
obtained:
.DELTA..sigma. ij 1 = C ijkl '1 .DELTA. kl 1 + .DELTA..sigma. ij R
1 where ( 31 ) C ijkl '1 .ident. C ijkl .infin. 1 + 1 .DELTA. t m =
1 M ijkl 1 .eta. ijkl m 1 ( 1 - - .DELTA. t / .rho. ijkl m 1 ) ( no
sum on i , j , k , l ) ( 32 ) .DELTA..sigma. ij R 1 = - k = 1 3 l =
1 3 A ijkl 1 ( 33 ) A ijkl 1 = m = 1 M ijkl 1 ( 1 - - .DELTA. t /
.rho. ijkl m 1 ) S ijkl m 1 ( t ) ( no sum on i , j , k , l ) ( 34
) S ijkl m 1 ( t ) = - .DELTA. t / .rho. ijkl m 1 S ijkl m 1 ( t -
.DELTA. t ) + .eta. ijkl m 1 .DELTA. kl 1 .DELTA. t ( 1 - - .DELTA.
t / .rho. ijkl m 1 ) ( no sum on i , j , k , l ) ( 35 )
##EQU00011##
where t denotes the previous time step and
.DELTA..epsilon..sub.kl.sup.1/.DELTA.t is determined from the
previous time step t. Note that all time dependence in the material
behavior resides in .DELTA..sigma..sub.ij.sup.R1 which is
recursively computed at each time step. Also note that the values
of S.sub.ijkl.sub.m.sup.1 from the previous time step must be
stored so that S.sub.ijkl.sub.m.sup.1 (t) can be determined
recursively. Moreover, S.sub.ijkl.sub.m.sup.1=0 for t.ltoreq.0.
Further details about the incrementalization of viscoelastic
constitutive equations can found in the specialized literature.
[0063] The same can be done for the cohesive zone
traction-displacement relation. Following the same procedure, one
can show that the cohesive zone constitutive equation can be
written in the following incremental form:
t i 1 ( t + .DELTA. t ) = t i 1 ( t ) + .DELTA. t i 1 ( 36 )
.DELTA. t i 1 = k ij 1 .DELTA. [ u j 1 ] + .DELTA. t i R 1 where (
37 ) k ij 1 = [ 1 - .alpha. 1 ( t + .DELTA. t ) ] .delta. i * 1 E
'1 .delta. ij ( no sum on i ) ( 38 ) E '1 .ident. E .infin. 1 + 1
.DELTA. t j = 1 M 1 .eta. j 1 ( 1 - - .DELTA. t / .rho. j 1 ) ( 39
) .DELTA. t i R 1 = [ 1 - .alpha. 1 ( t + .DELTA. t ) ] .delta. i *
1 [ - j = 1 M 1 ( 1 - - .DELTA. t / .rho. j 1 ) s ij 1 ( t ) ] -
.DELTA..alpha. 1 .delta. i * 1 { E .infin. 1 [ u i 1 ( t ) ] + j =
1 M 1 s ij 1 ( t ) } - .DELTA..alpha. 1 .sigma. i f 1 ( no sum on i
) ( 40 ) s ij 1 ( t ) = - .DELTA. t / .rho. j 1 s ij 1 ( t -
.DELTA. t ) + .eta. j 1 .DELTA. [ u i 1 ] .DELTA. t ( 1 - - .DELTA.
t / .rho. j 1 ) ( no sum on j ) ( 41 ) ##EQU00012##
where .delta..sub.ij is the Kronecker delta, M.sup.1 is the number
of Prony terms necessary to accurately describe the uniaxial
relaxation modulus of a single fibril in the cohesive zone, and
.DELTA.[u.sub.i.sup.1] is the increment of the displacement jump
across the cohesive zone boundaries, which rate is assumed to be
constant over each time step.
[0064] Now assume that the constitutive behavior of the global
object can be approximated by an incremental form similar to
equation 31:
.DELTA..sigma..sub.ij.sup.0=C.sub.ijkl.sup.0(t).DELTA..epsilon..sub.kl.s-
up.0+.DELTA..sigma..sub.ij.sup.R0 (42)
where C.sub.ijkl.sup.0(t) is the instantaneous (tangent)
constitutive tensor evaluated at the previous time step t which is
a function of time through its dependence on the amount of damage
accumulated at the local RVE, thus producing a nonlinear behavior
at the global scale. Substituting equations 42 and 31 into 28,
gives:
C ijkl 0 .DELTA. kl 0 + .DELTA..sigma. ij R 0 = 1 V 1 .intg. V 1 (
C ijkl '1 .DELTA. kl 1 + .DELTA. .sigma. ij R 1 ) V ( 43 )
##EQU00013##
where C.sub.ijkl.sup.0(t) is now written as C.sub.ijkl.sup.0 for
the shortness of notation. Now it is necessary to solve the above
equation for C.sub.ijkl.sup.0. The key to the following development
is to relate the field variables of the local scale to the external
boundary average quantities corresponding to field variables of the
global scale. Using the concept of localization tensors, it is
herein assumed that the local displacement field can be related to
the global strain tensor (external boundary average of the local
displacement field) according to:
u.sub.i.sup.1(t)=x.sub.j.sup.1.epsilon..sub.ij.sup.0(t)+.lamda..sub.ijk.-
sup.1(t).epsilon..sub.jk.sup.0(t) in V.sup.1 (44)
where .lamda..sub.ijk.sup.1 is called the localization tensor. The
first term in equation 44 corresponds to a homogeneous displacement
field which would be experienced by the RVE if it was a homogeneous
material. The second term, on the other hand, represents the
contribution to the local displacement field due to the history
dependent constitutive behavior and the existence of local
asperities and internal boundaries, including cracks and cohesive
zones. Note that the displacement vector is chosen to be described
in such a manner because it is generally used as the primary
variable in finite element formulations. Some researchers have have
used a fourth order strain localization tensor which relates
.epsilon..sub.kl.sup.1 to .epsilon..sub.kl.sup.0. This approach may
be sometimes impractical or inconvenient, especially when using the
finite element method to solve the local IBVP.
[0065] Note that the same procedure can be used in the case of
periodic displacement boundary conditions. It can be shown that
equation 44 is actually equivalent to the representation usually
admitted in the asymptotic homogenization theory, say
u.sub.i.sup.1=x.sub.ijk.sup.1.epsilon..sub.jk.sup.0, where
.chi.x.sub.ijk.sup.1 is the so-called influence function. Note
that, in the referred works, x.sub.ijk.sup.1 is not assumed to be a
function of time because neither viscoelasticity nor cracking are
considered.
[0066] Linearizing the local scale displacement field gives:
u i 1 ( t + .DELTA. t ) = u i 1 ( t ) + .differential. u i 1 ( t )
.differential. t .DELTA. t ( 45 ) ##EQU00014##
[0067] From equation 44 and the chain rule of differentiation, one
can show that:
.differential. u i 1 ( t ) .differential. t = x j 1 .differential.
ij 0 ( t ) .differential. t + .lamda. ijk 1 ( t ) .differential. jk
0 ( t ) .differential. t + .differential. .lamda. ijk 1 ( t )
.differential. t jk 0 ( t ) ( 46 ) ##EQU00015##
[0068] It is important to note that, from the incrementalization of
the global and local viscoelastic constitutive equations,
.differential..epsilon..sub.ij.sup.0(t)/.differential.t as well as
.differential.u.sub.i.sup.1(t)/.differential.t have been assumed to
be constant over the time step. However, no assumption will be made
about the behavior of .lamda..sub.ijk.sup.1(t) and
.differential..lamda..sub.ijk.sup.1(t)/.differential.t in order to
maintain the same degree of accuracy, but equation 46 is redefined
as follows:
.differential. u i 1 ( t ) .differential. t = x j 1 .differential.
ij 0 ( t ) .differential. t + .lamda. ijk 1 ( t ) .differential. jk
0 ( t ) .differential. t + .differential. u i R 1 ( t )
.differential. t where ( 47 ) .differential. u i R 1 ( t )
.differential. t .ident. .differential. .lamda. ijk 1 ( t )
.differential. t jk 0 ( t ) ( 48 ) ##EQU00016##
and one can show from the previous assumptions that
.differential.u.sub.i.sup.R1(t)/.differential.t is constant over
the time step, thus giving:
.DELTA.u.sub.i.sup.1=x.sub.j.sup.1.DELTA..epsilon..sub.ij.sup.0+.lamda..-
sub.ijk.sup.1(t).DELTA..epsilon..sub.jk.sup.0+.DELTA.u.sub.i.sup.R1
(49)
[0069] Substituting 49 into the infinitesimal strain-displacement
relationship for scale 1, then into 43 yields:
C ijkl 0 .DELTA. kl 0 + .DELTA..sigma. ij R 0 = 1 V 1 .intg. V 1 {
C ijkl '1 + C ijpq '1 [ 1 2 ( .lamda. pkl , q 1 + .lamda. qkl , p 1
) ] } V .DELTA. kl 0 + 1 V 1 .intg. V 1 ( C ijkl '1 .DELTA. kl R 1
+ .DELTA..sigma. ij R 1 ) V ( 50 ) ##EQU00017##
where symmetry of the local constitutive tensor C.sub.ijkl.sup.1
and the fact that .DELTA..epsilon..sub.kl.sup.0 is constant in
V.sup.1 have been used, .lamda..sub.ijk,l.sup.1(t) is written as
.lamda..sub.ijk,l.sup.1 to simplify equations, and
.DELTA..epsilon..sub.ij.sup.R1=1/2(.DELTA.u.sub.ij.sup.R1+.DELTA.u.sub.j-
,i.sup.R1) (51)
Since .DELTA..epsilon..sub.kl.sup.0 is arbitrary, one obtains:
C ijkl 0 = 1 V 1 .intg. V 1 { C ijkl '1 + C ijpq '1 [ 1 2 ( .lamda.
pkl , q 1 + .lamda. qkl , p 1 ) ] } V ( 52 ) .DELTA..sigma. ij R 0
= 1 V 1 .intg. V 1 ( C ijkl '1 .DELTA. kl R 1 + .DELTA..sigma. ij R
1 ) V ( 53 ) ##EQU00018##
[0070] Note that the global scale tangent constitutive tensor,
C.sub.ijkl.sup.0, is affected by heterogeneity, internal boundaries
and cracks at the local scale through the localization tensor,
.lamda..sub.ijk.sup.1, and affected by viscoelasticity through
C'.sub.ijkl.sup.1 and .lamda..sub.ijk.sup.1. It is also important
to note that the global scale history dependence, expressed by
.DELTA..sigma..sub.ij.sup.R0, also depends on the heterogeneity,
internal boundaries and cracks initiated at the local scale through
.DELTA.u.sub.i.sup.R1, defined by equation 48, which depends on the
rate of change of .lamda..sub.ijk.sup.1(t).
[0071] Also note that the homogenized tangent constitutive tensor,
C.sub.ijkl.sup.0, is only equal to the volume average of the local
scale constitutive tensor if
.lamda..sub.ijk.sup.1=constant=0.sup.1, meaning that the local
scale object is actually homogeneous. Therefore, the rule of
mixtures (volume average of the local scale constitutive tensor),
which does not account for the localization of field variables,
should be used only as a rough estimate which may be very
inaccurate in some cases. Note that, in order to satisfy the
displacement boundary conditions,
u.sub.i.sup.1(t)=x.sub.j.sup.1.epsilon..sub.ij.sup.0(t), all
components of .lamda..sub.ijk.sup.1 have to be zero at least on
some regions of .differential.V.sub.E.sup.1 (or all over
.differential.V.sub.E.sup.1 if no heterogeneity (cracks and
inclusions) is allowed to cross the external boundary of the RVE
and periodic fluctuations on .differential.V.sub.E.sup.1 are null).
Thus, if .lamda..sub.ijk.sup.1 is constant over the volume V.sup.1,
it has to be zero.
[0072] It can be shown that the expression for computing the
homogenized tangent constitutive tensor is the same for both
elastic and viscoelastic media containing cracks. The only
difference is the additional homogenized incremental history
dependent stress tensor that appears in the viscoelastic
formulation.
[0073] It can also be shown that equation 52 is equivalent to the
one obtained from the asymptotic homogenization theory for both
linear elastic and elasto-plastic materials, even though cracking
has not been considered in these works.
7.3 Finite Element Formulation for Spatial Dependence
[0074] The finite element formulation for the global scale IBVP
follows the standard formulation presented in reference books on
the subject, and in this case, an implicit time integration scheme
is used to solve the global scale (roadway) IBVP.
[0075] In this section, the finite element formulation for the
local IBVP is developed. For the local IBVP, one should now seek
the solution for the localization tensor .lamda..sub.ijk.sup.1(t)
and .DELTA.u.sub.i.sup.R1 so that the displacement increments can
be computed using equation 49, and C.sub.ijkl.sup.0 and
.DELTA..sigma..sub.ij.sup.R0 can be computed using equations 52 and
53, respectively.
[0076] Consider the variational form of conservation of linear
momentum at the local scale, neglecting the body force vector:
.intg..sub.V.sub.1.sigma..sub.ji,j.sup.1.delta.u.sub.i.sup.1dV=0
(54)
[0077] From the chain rule of differentiation and conservation of
angular momentum, the above equation results in:
.intg..sub.V.sub.1(.sigma..sub.ji.sup.1.delta.u.sub.i.sup.1).sub.jdV=.in-
tg..sub.V.sub.1.sigma..sub.ij.sup.1.delta..epsilon..sub.ij.sup.1dV
(55)
[0078] Applying the divergence theorem to the above equation at
time t+.DELTA.t and using Cauchy's formula
(t.sub.i.sup.1=.sigma..sub.ji.sup.1n.sub.j.sup.1) results in:
.intg. V 1 .sigma. ij 1 ( t + .DELTA. t ) .delta. ij 1 ( t +
.DELTA. t ) V = .intg. .differential. V E 1 t i 1 ( t + .DELTA. t )
.delta. u i 1 ( t + .DELTA. t ) S + .intg. .differential. V I 1 t i
1 ( t + .DELTA. t ) .delta. u i 1 ( t + .DELTA. t ) S ( 56 )
##EQU00019##
[0079] Besides the incremental equations 26, 36 and 45, consider
that the strain field can also be approximated by the following
incremental form:
.epsilon..sub.ij.sup.1(t+.DELTA.t)=.epsilon..sub.ij.sup.1(t)+.DELTA..eps-
ilon..sub.ij.sup.1 (57)
[0080] Thus, substituting into 56, but noting that the variation of
.epsilon..sub.ij.sup.1(t) and u.sub.i.sup.1(t) are null because
these are known quantities, yields:
.intg. V 1 .sigma. ij 1 ( t ) .delta..DELTA. ij 1 V + .intg. V 1
.DELTA..sigma. ij 1 .delta..DELTA. ij 1 V = .intg. .differential. V
E 1 t i 1 ( t ) .delta..DELTA. u i 1 S + .intg. .differential. V E
1 .DELTA. t i 1 .delta..DELTA. u i 1 S + .intg. .differential. V I
1 t i 1 ( t ) .delta..DELTA. u i 1 S + .intg. .differential. V I 1
.DELTA. t i 1 .delta..DELTA. u i 1 S ( 58 ) ##EQU00020##
[0081] Also noting that,
.intg. .differential. V E 1 .differential. V I 1 t i 1 ( t )
.delta..DELTA. u i 1 S = .intg. V 1 .sigma. ji , j 1 ( t )
.delta..DELTA. u i 1 V + .intg. V 1 .sigma. ji 1 ( t )
.delta..DELTA. ij 1 V = .intg. V 1 .sigma. ij 1 ( t )
.delta..DELTA. ij 1 V ( 59 ) ##EQU00021##
results in:
.intg..sub.V.sub.1.DELTA..sigma..sub.ij.sup.1.delta..DELTA..epsilon..sub-
.ij.sup.1dV=.intg..sub..differential.V.sub.E.sub.1.DELTA.t.sub.i.sup.1.del-
ta..DELTA.u.sub.i.sup.1dS+.intg..sub..differential.V.sub.l.sub.1.DELTA.t.s-
ub.i.sup.1.delta..DELTA.u.sub.i.sup.1dS (60)
[0082] Under the circumstances where the previously discussed
assumptions regarding the global and local length scales are valid,
either displacements or tractions can be applied on the boundary of
the RVE. If displacements are applied, then the first integral on
the right-hand side of the equation above is zero because the
variation of the displacement field on the external boundary is
null. If tractions are applied instead, this integral is also zero
because the tractions are self-equilibrated over the external
boundary of the RVE. Thus,
.intg..sub.V.sub.1.DELTA..sigma..sub.ij.sup.1.delta..DELTA..epsilon..sub-
.ij.sup.1dV=.intg..sub..differential.V.sub.l.sub.1.DELTA.t.sub.i.sup.1.del-
ta..DELTA.u.sub.i.sup.1dS (61)
[0083] Now notice that the virtual work done by tractions over the
internal boundary (right-hand side of equation above) is zero
unless there is separation along that boundary. Therefore, the
integral over the internal boundaries can be reduced to an integral
over cracks and cohesive zones. Since
n.sub..differential.V.sub.CZ.sub.1.sup.1=-n.sub..differential.V.sub.l.sub-
.1.sup.1, and
.DELTA.t.sub.i(n.sub..differential.V.sub.CZ.sub.1.sup.1)=-.DELTA.t.sub.i(-
n.sub..differential.V.sub.l.sub.1.sup.1), the equation above can be
rewritten as:
.intg..sub.V.sub.1.DELTA..sigma..sub.ij.sup.1.delta..DELTA..epsilon..sub-
.ij.sup.1dV+.intg..sub..differential.V.sub.CZ.sub.1.DELTA.t.sub.i.sup.1.de-
lta..DELTA.u.sub.i.sup.1dS=0 (62)
[0084] Substituting 49 into the strain-displacement relationship at
the local scale, then into 31, gives:
.DELTA..sigma..sub.ij.sup.1={C.sub.ijkl.sup.1+C.sub.ijpq.sup.1[1/2(.lamd-
a..sub.pkl,q.sup.1+.lamda..sub.qkl,p.sup.1]}.DELTA..epsilon..sub.kl.sup.0+-
C.sub.ijkl.sup.1.DELTA..epsilon..sub.kl.sup.R1+.DELTA..sigma..sub.ij.sup.R-
1 (63)
Also, from equation 19,
.DELTA.[u.sub.i.sup.1]=[x.sub.j.sup.1.DELTA..epsilon..sub.ij.sup.0+.lamd-
a..sub.ijk.sup.1.DELTA..epsilon..sub.jk.sup.0+.DELTA.u.sub.i.sup.R1].sub.s-
urface
1-[x.sub.j.sup.1.DELTA..epsilon..sub.ij.sup.0+.lamda..sub.ijk.sup.1-
.DELTA..epsilon..sub.jk.sup.0+.DELTA.u.sub.i.sup.R1].sub.surface 2
(64)
or
.DELTA.[u.sub.i.sup.1]=[.lamda..sub.ijk.sup.1|.sub.surface
1-.lamda..sub.ijk.sup.1|.sub.surface
2].DELTA..epsilon..sub.jk.sup.0+[.DELTA.u.sub.i.sup.R1|.sub.surface
1-.DELTA.u.sub.i.sup.R1|.sub.surface 2] (65)
where one should note that the initial coordinates, x.sub.i.sup.1,
of the cohesive zone surfaces coincide.
[0085] The incremental localization field can now be discretized in
space by expressing it as a function of its nodal values:
[.lamda..sup.1].sub.e=[N.sup.-1].sub.e[{tilde over
(.lamda.)}.sup.1].sub.e (66)
where the above equation has been written in matrix form,
[N.sup.1].sub.e is the matrix of nodal shape functions for each
finite element, the subscript e denotes quantities defined for each
finite element, and the over-tilde indicates nodal value. It is
important to note that the shape functions used to interpolate
.lamda..sub.ijk.sup.1 are the same as for the displacement field.
Therefore, one can also write:
{.delta..DELTA.u.sup.1}.sub.e=[N.sub.CZ.sup.1].sub.e{.delta..DELTA.
.sup.1}.sub.e on .differential.V.sub.CZ.sup.1 (67)
{.DELTA.[u.sup.1]}.sub.e=[N.sub.CZ.sup.1].sub.e[{tilde over
(.lamda.)}.sup.1].sub.e
{.DELTA..epsilon..sup.0}.sub.e+[N.sub.CZ.sup.1].sub.e{.DELTA.
.sup.R1}.sub.e (68)
{.delta..DELTA..epsilon..sup.1}.sub.e=[B.sup.1].sub.e{.delta..DELTA.
.sup.1}.sub.e (69)
where [B.sup.1].sub.e is the matrix of derivatives of the shape
functions with respect to the local scale spatial coordinates.
Substituting equations 37, 63-69 into 62 and noting that
.delta..DELTA. .sub.i.sup.1 and .DELTA..epsilon..sub.ij.sup.0 are
arbitrary quantities results in the following two sets of linear
equations:
[{tilde over
(.lamda.)}.sup.1].sub.a=-[K.sup.1].sub.a.sup.-1[G.sup.1].sub.a
(70)
{.DELTA. .sup.R1}.sub.a=[K.sup.1].sub.a.sup.-1{F.sup.R1}.sub.a
(71)
where the subscript a denotes quantities assembled over the entire
finite element mesh.
[ K 1 ] a = .intg. V 1 [ B 1 ] T [ C 1 ] [ B 1 ] V + .intg.
.differential. V CZ 1 [ N CZ 1 ] T [ k CZ 1 ] [ N CZ 1 ] S ( 72 ) [
G 1 ] a = .intg. V 1 [ B 1 ] T [ C 1 ] V ( 73 ) { F R 1 } a = - {
.intg. V 1 [ B 1 ] T { .DELTA..sigma. R 1 } V + .intg.
.differential. V CZ 1 [ N CZ 1 ] T { .DELTA. t R 1 } S } ( 74 )
##EQU00022##
where the subscripts e have been dropped for a simpler
notation.
[0086] Once .lamda..sub.ijk.sup.1(t) and .DELTA.u.sub.i.sup.R1 have
been obtained using equations 70 and 71, respectively, one can
calculate the incremental displacements using equation 49 and
compute the global scale (homogenized) tangent constitutive tensor,
C.sub.ijkl.sup.0, and the global scale history dependent stress
tensor, .DELTA..sigma..sub.ij.sup.R0 from equations 52 and 53,
respectively. Then standard FE procedures can be used to compute
the incremental strain and stress tensors from the incremental
displacements, thus minimizing the code implementation effort.
[0087] Now one can clearly see that the incremental localization
tensor, .lamda..sub.ijk.sup.1, and consequently the homogenized
tangent constitutive tensor, C.sub.ijkl.sup.0, are affected by
cracks and cohesive zones through the stiffness matrix, [K.sup.1].
Also note that the global scale history dependent term,
.DELTA..sigma..sub.ij.sup.R0 also depends on the damage accumulated
at the local scale through .DELTA.u.sub.i.sup.R1. It is important
to note that the stiffness matrix, [K.sup.1], is the same stiffness
matrix usually assembled in common FE algorithms.
* * * * *