Discontinuous enrichment in finite elements with a partition of unity method.
We present an approximate analytical method to evaluate efficiently and accurately the call blocking probabilities in wavelength routing networks with multiple classes of calls. The model is fairly general and allows each source-destination pair to service calls of different classes, with each call occupying one wavelength per link. Our approximate analytical approach involves two steps. The arrival process of calls on some routes is first modified slightly to obtain an approximate multiclass network model. Next, all classes of calls on a particular route are aggregated to give an equivalent single-class model. Thus, path decomposition algorithms for single-class wavelength routing networks may be readilt extended to the multiclass case. This article is a first step towards understanding  the issues arising in wavelength routing networks that serve multiple classes of customers.
Introduction
The modeling of evolving discontinuities with the finite element method is
cumbersome due to the need to update the mesh topology to match the
geometry of the discontinuity. In this paper, we present a technique to model
discontinuities in the finite element framework in a general fashion. The
essential feature is the incorporation of enrichment functions which contain a
presently Assistant Professor of Civil and Environmental Engineering, Duke University
y Research Associate in Mechanical Engineering, Northwestern University
z Walter P. Murphy Professor of Mechanical Engineering, Northwestern University
discontinuous field. In the application of the technique to fracture mechanics,
functions spanning the appropriate near-tip crack field can also be included
to improve accuracy. The enrichment of the finite element approximation in
this manner provides for both the modeling of discontinuities and accurate
moment intensity factors with minimal computational resources.
The concept of incorporating crack fields in a finite element context is
not new, see for example [1]. In addition, there are several well established
techniques for modeling cracks and crack growth such as boundary element
methods, finite elements with continuous remeshing [2], and meshless methods
[3]. Recently, the trend has focused on the development of finite element
methods which model discontinuities independently of element boundaries.
These include the incorporation of a discontinuous mode in an assumed strain
framework [4], and enrichment with near-tip fields for crack growth with
minimal remeshing [5]. The latter method has recently been extended by
enriching with a discontinuous function behind the crack tip [6], [7], such
that no remeshing is necessary.
In this paper, the method of discontinuous enrichment is cast in a general
framework, and we illustrate how both two-dimensional and plate formulations
can be enriched to model cracks and crack growth. The enrichment
of the approximation with discontinuous near-tip fields requires a mapping
technique, and so an alternative near-tip function is developed. The present
method offers several advantages over competing techniques for modeling
crack growth. In contrast to traditional finite element methods, this technique
incorporates the discontinuity of the crack independently of the mesh,
such that the crack can be arbitrarily located within an element. The present
technique has a distinct advantage over boundary element methods as it is
readily applicable to non-linear problems, anisotropic materials, and arbitrary
geometries. The method does not require any remeshing for crack
growth, and as it is an extension of the finite element method, it can exploit
the large body of finite element technology and software. Specific examples
of augmenting well established two-dimensional and plate elements with both
discontinuous and asymptotic near-tip functions are presented.
The present technique exploits the partition of unity property of finite elements
first cited by [8], which allows global enrichment functions to be locally
incorporated into a finite element approximation. A standard approximation
is 'enriched' in a region of interest by the global functions in conjunction with
additional degrees of freedom and the local nodal shape functions. The application
of this idea to capture a specific frequency band in dynamics can
be found in [9]. The utility of the method has found application in solving
the scalar Laplacian problem for domains with re-entrant corners [10]. In
the context of fracture mechanics, the appropriate enrichment functions are
the near-tip asymptotic fields and a discontinuous function to represent the
jump in displacement across the crack line. In contrast to the work of [4], the
enrichment is not through an assumed-strain method, so the displacement
field is continuous along either side of the crack.
This paper is organized as follows. Following this introduction, we review
the construction of an enriched approximation and we develop a discontinuous
near-tip function which does not require a mapping. For specific applications
we consider two-dimensional linear elastic fracture mechanics and the
fracture of Mindlin-Reissner plates in Section 3. Numerical results to verify
the accuracy of the formulation are given in Section 4, with a summary and
some concluding remarks provided in the last section.
2 Construction of a Finite Element
Approximation With Discontinuities
In this section, we present the construction of a finite element approximation
with discontinuous enrichment. Emphasis is placed on modeling cracks, in
which a standard approximation is enriched with both the asymptotic near-tip
functions and a discontinuous 'jump' function. The incorporation of
discontinuous near-tip functions requires a mapping for kinked cracks, and
so an alternative near-tip function is presented. The manner in which nodes
are selected for enrichment and the modifications to the numerical integration
of the weak form are also given.
2.1 General Form
To introduce the concept of discontinuous enrichment, we begin by considering
the
domain\Omega bounded by \Gamma with an internal boundary \Gamma c as shown in
Fig. 1a. We are interested in the construction of a finite element approximation
to the field u
2\Omega which can be discontinuous along \Gamma c .
Consider the uniform mesh of N nodes for the domain shown in Fig. 1b
which does not model the discontinuity. The discrete approximation u h to
the function u takes the form
I
I ((x))u I (1)
where N I is the shape function for node I in terms of the parent coordinates
((x)), and u I is the vector of nodal degrees of freedom. The nodal shape
function N I is non-zero over the support of node I, defined to be union of
the elements connected to the node.
We now pose the question of how to best incorporate the discontinuity in
the field along \Gamma c . The traditional approach is to change the mesh to conform
to the line of discontinuity as shown in Fig. 1c, in which the element edges
align with \Gamma c . While this strategy certainly creates a discontinuity in the
approximation, it is cumbersome if the line \Gamma c evolves in time, or if several
different configurations for \Gamma c are to be considered.
In this paper we propose to model the discontinuity along \Gamma c with extrinsic
enrichment [11], in which the standard approximation (1) is modified as
I ((x))@u I
a Il G l (x)A (2)
where G l (x) are enrichment functions, and a Il are additional nodal degrees
of freedom for node I. In the above, the total number of enriched degrees of
freedom for a node is denoted by nE (I). If the enrichment functions G l are
discontinuous along the boundary \Gamma c , then the finite element mesh does not
need to model the discontinuity. For example, the uniform mesh in Fig. 1d
is capable of modeling a jump in u when the circled nodes are enriched with
functions which are discontinuous across \Gamma c .
The above form of a finite element approximation merits some discussion.
We note that the enrichment functions G l are written in terms of the global
coordinates x, but that they are multiplied by the nodal shape functions N I .
In this fashion the additional enrichment takes on a local character. This
concept of multiplying global functions by the finite element partition of
unity was first suggested in [8]. The change in the form of the approximation
from (1) to (2) is only made locally in the vicinity of a feature of interest,
such as a discontinuity.
We now turn to the precise form of the enrichment functions used to model
discontinuous fields, with the goal of modeling cracks and crack growth.
Three distinct regions are identified for the crack geometry, namely the crack
interior and the two near-tip regions as shown in Fig.2. In the set I of all
nodes in the mesh, we distinguish three different sets which correspond to
each of these regions. The set J is taken to be the set of nodes enriched for
the crack interior, and the sets K 1 and K 2 are those nodes enriched for the
first and second crack tips, respectively. The precise manner in which these
sets are determined from the interaction of the crack and the mesh geometry
is given in Section 2.3.
The enriched approximation takes the form
I
I u I
l
l
where b J and c 1
Kl , c 2
are nodal degrees of freedom corresponding
to the enrichment functions H(x), F 1
l (x) and F 2
l (x), respectively. The
function H(x) is discontinuous across the crack line, and the sets F 1
l (x) and
l (x) consist of those functions which span the near-tip asymptotic fields.
For two-dimensional elasticity, these are given by
fF l (r; ')g 4
r sin( '
r cos( '
r sin( '
r cos( '
where (r; ') are the local polar coordinates for the crack tip [12]. Note that
the first function in (40),
rsin( '), is discontinuous across the crack faces
whereas the last three functions are continuous. The form of the near-tip
functions for plates is similar and is developed in Section 3.
The jump function H(x) is defined as follows. The crack is considered to
be a curve parametrized by the curvilinear coordinate s, as in Fig. 3. The
origin of the curve is taken to coincide with one of the crack tips. Given a
point x in the domain, we denote by x   the closest point on the crack to x.
At x   , we construct the tangential and normal vector to the curve, e s and
with the orientation of e n taken such that e s " e . The function
H(x) is then given by the sign of the scalar product In the
case of a kinked crack, the cone of normals at x   needs to be considered (see
[6]). Roughly speaking, the function H(x) takes the value of 1 'above' the
crack, and \Gamma1 'below' the crack.
2.2 An Alternative Near-Tip Function
The jump function H(x) is in general not capable of representing the discontinuity
in the displacement field along the entire crack geometry. For
example, if the crack tip is not aligned with an element edge, then a near-tip
function must also be used (see [6]). For cracks which are not straight, a
mapping is required to align the near-tip discontinuities with the crack edges.
Due to the local form of the enrichment, the mapping procedure is only necessary
in those elements with nodes enriched with the near-tip functions. In
this section, we review the mapping procedure and present an alternative
near-tip function.
The crack is modeled as a series of straight line segments connecting ver-
tices, with new crack segments added as the crack grows. The discontinuities
in the near-tip fields are aligned with each segment by using a procedure developed
in [12] and [5]. In this procedure, the discontinuity in the near-tip
functions are aligned with the crack by a mapping technique that rotates
each section of the discontinuity onto the crack model.
A key step in technique is the modification of the angle ' in F l (r; ').
Given a point define an angle
terms of the angle of
the segment ' R (see Fig 4) and the sampling point angle ff(x
The coordinates of the sampling point are mapped to coordinates in
the crack tip frame ("x; " y) as shown in Fig. 5:
r cos(
r sin(
where l is the distance between r is as shown in the
figure. The variables r and ' in the enrichment functions are then computed
in terms of the local ("x; " y) coordinates. This procedure is repeated similarly
for each segment of the crack and the sequence of mappings leaves the length
of the crack invariant.
In [5], the entire crack was modeled with the near-tip fields and the above
mapping procedure. The use of the discontinuous function H(x) eliminates
the need for the mapping on the crack interior, so that the above procedure
is only necessary at each crack tip. In the following, we propose the use of a
smooth 'ramp' function in conjunction with the function H(x) to model the
near-tip region.
Consider the following function defined in terms of the crack tip coordinates
x
l c
l c
where the length l c is taken to be the characteristic length of the element
containing the crack tip. The above function and its derivative vanishes at
the crack tip.
When this smooth ramp function is multiplied by the function H(x), i.e.,
~
a near-tip function which is discontinuous across the crack edges and vanishes
in front of the crack tip results. This function in turn does not require any
mapping to align the discontinuity with the crack edges. The function ~
R is
shown in Fig. 6. in the vicinity of a crack tip and two consecutive segments.
It is clear from the figure that the resulting near-tip function is continuous
in the
domain\Omega and discontinuous across the crack line.
The above near-tip function is useful on several levels. In the first in-
stance, in conjunction with the function H(x) for the crack interior it is
perhaps the simplest means to model the entire crack discontinuity. In addi-
tion, for non-linear problems the exact near-tip functions may not be known.
The concept of multiplying a smooth function which vanishes at the crack
tip by the jump function H(x) can also be extended to three-dimensional
problems. In linear elastic fracture mechanics, however, the incorporation
of the asymptotic near-tip fields is still useful to obtain greater accuracy at
the crack tip. To some extent, this advantage can be maintained by simply
replacing the first function in (4) with the function ~
R:
fF l (r; ')g 4
r cos( '
r sin( '
r cos( '
The last three functions need not be mapped into the crack faces, as they
are all continuous in the
domain\Omega\Gamma
2.3 Node Selection for Enrichment
In the preceding development, three distinct regions were identified for en-
richment, corresponding to the nodal sets J , K 1 and K 2 . In this section we
define these sets precisely, and present the methodology by which nodes are
identified for inclusion in each set.
We begin with some preliminary notations. The support of node I is
denoted by ! I , with closure ! I . Essentially, a node's support is the open set
of element domains connected to the node, and the closure is the closed set
which includes the outer boundary. The distinction as it applies to nodal
selection will be discussed shortly. We also denote the location of crack tips
respectively, and by C the geometry of the crack.
With these definitions, the sets J , K 1 and K 2 are defined as follows.
The sets K 1 and K 2 consist of those nodes whose support closure contains
crack tip 1 or 2, respectively. The set J is the set of nodes whose support is
intersected by the crack and do not belong to K 1 or
Note that the set J consists of nodes whose support, as opposed to support
closure, is intersected by the crack. This distinction implies that if the crack
only intersects the boundary of a node's support, the node will not be enriched
with the function H(x). This prevents any nodes from being enriched
with a constant function (either -1 or +1) over their entire support, which
is important in order to avoid creating a linear dependency in the approxi-
mation. We note that for the alternative near-tip function proposed in the
previous section, the support closure must be changed to the open set for
similar reasons.
In practice the above sets are determined as follows. All elements intersected
by the crack are first determined. From this set of elements, we
distinguish three disjoint sets of 'tip elements' (for either tip 1 or 2) and
'interior elements'. The set of tip elements are given by those which contain
either crack tip. The nodes of the 'tip elements' correspond to either set K 1
or K 2 . The nodes of the 'interior elements' in turn correspond to the set J .
Fig. 7 illustrates the nodes that are selected as tip nodes and interior nodes
for the cases of a uniform mesh and an unstructured mesh.
An additional step is taken to remove those nodes from the set J whose
support closure, but not support, is intersected by the crack. For this pur-
pose, subpolygons which align with both the crack and element boundaries
are generated as shown in Fig. 8b. These subpolygons are generated easily
enough by triangulating the polygons formed from the intersection of the
crack and element boundaries.
The generation of the subtriangles allows the computation of the amount
of a node's support `above' and 'below' the crack, which can then be compared
against a tolerance. In two-dimensional analysis, we denote the area
of a nodal support by A ! , which is calculated from the sum of the areas of
each element connected to the node. With the aid of the subtriangles, we
also calculate the nodal support area above A ab
and below A be
! the crack. We
then calculate the ratios
A ab
(13a)
A be
If either of these ratios are below a tolerance (we have used 0.01, or 1%), the
nodes are removed from the set J .
2.4 Numerical Integration of the Weak Form
For elements cut by the crack and enriched with the jump function H(x),
we make a modification to the element quadrature routines for the assembly
of the weak form. As the crack is allowed to be arbitrarily oriented in
an element, standard Gauss quadrature may not adequately integrate the
discontinuous field. For those nodes in the set J , it is important that the
quadrature scheme accurately integrate the contributions to the weak form
on both sides of the discontinuity. If the integration of the discontinuous
enrichment is indistinguishable from that of a constant function, spurious
singular modes can appear in the system of equations. In this section, we
present the modifications made to the numerical integration scheme for elements
cut by a crack.
The discrete weak form is normally constructed with a loop over all ele-
ments, as the domain is approximated by
where m is the number of elements,
and\Omega e is the element subdomain. For
elements cut by a crack, we define the element subdomain to be a union of
a set of subpolygons whose boundaries align with the crack geometry:
s
s denotes the number of subpolygons for the element. The subtriangles
shown in Fig. 8 already generated for the selection of the interior nodes
also work well for integration. It is emphasized that the subpolygons are
only necessary for integration purposes; no additional degrees of freedom are
associated with their construction. In the integration of the weak form, the
element loop is replaced by a loop over the subpolygons for those elements
cut by the crack.
3 Application to Fracture Mechanics
In this section, we review the pertinent equations for linear elastic fracture
mechanics. In this paper, emphasis is placed on plate fracture, although
some two-dimensional plane-strain studies are discussed. A key difference in
the plate formulation is the enrichment of the displacement components with
different sets of near-tip functions. After reviewing the governing equations
for Mindlin-Reissner plates, we examine the form of the asymptotic crack tip
fields. A domain form of the J-integral for plates is derived for the calculation
of the energy release rate and the moment intensity factors. Finally, the
enriched finite element approximation is presented.
3.1 Mindlin-Reissner Plate Formulation
Two main formulations exist to model a plate: the classical theory or Kirchhoff
plate theory and the Mindlin-Reissner plate theory. Allowing three
boundary conditions instead of two for the Kirchhoff theory, the Mindlin
theory gives a more realistic shear and moment distribution around a crack
tip (see [13] and [14]). A summary of the Mindlin theory follows.
3.1.1 Governing Equations
There are several different ways to introduce the Mindlin theory. As we
are also interested in examining problems in two-dimensional elasticity, the
theory is presented here as a degeneration of the three-dimensional elasticity
problem using the principal of virtual work with the appropriate kinematic
assumptions.
Consider a plate of thickness t whose mid-plane lies in the x
plane. The conventions used throughout this paper are shown in Fig. 9.
The main assumptions of the Mindlin theory state that the in plane dis-
placements, u 1 and u 2 vary linearly through the thickness with the section
rotations / 1 and / 2 . In addition, the normal stress oe 33 is assumed to vanish in
the domain. For the sake of simplicity, we make the additional assumptions
that the surface of the plate and any crack faces are traction-free.
In the (e 1 is the unit normal vector to the plate,
the deformation components at a point are given by
where w is the transverse displacement and / 1 and / 2 are the rotations
about the x 2 and x 1 axes, respectively. The above can be expressed in a
more compact form as
The strain is given by2 (ru
with the bending contribution
and a shear contribution
We note that the x 3 related components are zero for both ffl b and ffl s .
The virtual internal work is defined by
Z
d\Omega (21)
where oe is the symmetric stress tensor, and ffiu is an arbitrary virtual displacement
from the current position. After a few manipulations, we obtain
the relation
where the  superscript indicates a reduction of the operator to the in plane
component and oe s is the shear stress vector oe
Making the substitution (22) into (21) and integrating through the thickness
gives the work expression
Z
A
where the moment M and shear Q are defined by
Z t=2
\Gammat=2
Z t=2
\Gammat=2
oe s dx 3 (24)
The virtual external work is composed of the action of the bending and
twisting moments gathered in a couple vector C, and of the shear traction
T . We assume there is no external pressure acting on the plate. The virtual
external work is then given by
Z
Z
Equating the internal and external virtual work, and applying the divergence
theorem yields the equilibrium equations
in\Omega
r
and the traction boundary conditions on \Gamma
where n is the unit outward normal to the boundary.
The constitutive relationships are obtained by energetic equivalence between
the plate and the three-dimensional model. Assuming the plate is
made of an isotropic homogeneous elastic material of Young's modulus E
and of Poisson's ratio , the constitutive relations are given by4 M 11
22
ffl b22
and
Ekt
5=6 is a correction factor which accounts for the parabolic variation
of the shear stresses through the plate thickness.
These are rewritten in a more compact form using the fourth order bending
stiffness tensor D b and the second order shear stiffness tensor
3.1.2 Weak Form
Let the boundary \Gamma be divided into a part \Gamma u on which displacement boundary
conditions are imposed and a part \Gamma t on which loads are applied with
the restrictions
The kinematics constraints are given by a prescribed transverse displacement
w and prescribed rotations / while the loads come from the prescribed couples
C and prescribed shear tractions T . As in Section 2, we also designate
c as an internal boundary across which the displacement field is allowed to
be discontinuous.
be the space of kinematically admissible transverse displacements
and rotations:
where V is a space of sufficiently smooth functions on \Omega\Gamma The details on
this matter when the domain contains an internal boundary or re-entrant
corner may be found in [15] and [16]. We note that the space V allows for
discontinuous functions across the crack line.
The space of test functions is defined similarly as:
The weak form is to find (w; /) 2 V g such that
Z
(D
Z
(D s s(w; /)) \Delta s(ffiw; ffi/)
Z
Z
It can be shown that the above is equivalent to the equilibrium equations (26)
and traction boundary conditions (27). When the space V is discontinuous
along \Gamma c , the traction-free conditions on the crack faces are also satisfied.
In contrast to boundary element techniques, this enables the method to be
easily extended to non-linear problems.
In the finite element method, the space V is approximated with a finite
dimensional space V h ae V. The space V h is typically made discontinuous
across \Gamma c by explicitly meshing the surface, as in Fig. 1c. In the present
method, the approximating space is constructed with discontinuous enrichment

3.2 Plate Fracture Mechanics
Consider the problem of a through crack in a plate as shown in Fig. 10,
where for convenience we adopt a local polar coordinate system centered at
the crack tip. In contrast to the stress intensity factors obtained in classical
linear elasticity, in plate theory the quantities of interest are moment and
shear force intensity factors. The moment intensity factors are denoted by
K I and K II , while the shear force intensity factor is denoted by K III . These
are defined as
The relationship between these factors and the energy release rate G is similar
to the three-dimensional theory:
\Theta
I
II
5Et
III (37)
The form of the asymptotic near-tip displacement fields differs significantly
from the three-dimensional theory. In particular, the transverse displacement
w is only singular when a K III mode is present. The asymptotic
displacement fields in Mindlin-Reissner plate theory are given in [17], and
they are provided here for the sake of completeness.
5h
\Gamma3
(38a)
\Gammasin(
\Gamma2cos(
(38c)
For the purposes of defining the near-tip enrichment functions in the plate
theory, we consider only the terms proportional to
r for the rotations / 1
and / 2 . For the transverse displacement, we consider terms proportional to
both
r and r 3=2 . With these restrictions, the near-tip fields are contained
in the span of the sets
where
fg I (r; ')g j
ae p
oe
(40a)
ff I (r; ')g j
ae p
oe
(40b)
The discrete approximation for the plate which incorporates the above near-tip
functions is presented in Section 3.4.
3.3 Domain Form of the J-integral
Several different domain and path-independent integrals have been developed
for the extraction of mixed mode moment and shear force intensity factors in
plates [17]. These integrals typically consist of a contour integral enclosing
the crack-tip singularity. With finite elements, the numerical evaluation of
these integrals usually involves some kind of smoothing technique, as the
required field quantities are discontinuous at element interfaces.
In this section, we illustrate the use of a weighting function q to recast
these line integrals into their equivalent domain form. The development
presented here closely follows that given for two-dimensional elasticity in
[18]. The domain forms of crack contour integrals are particularly well suited
for use with finite elements, as the same quadrature points used for the
integration of the weak form can be used to calculate the domain integral.
The construction of additional quadrature points or the use of a smoothing
procedure is not required.
Consider the open contour \Gamma surrounding a through crack as shown in
Fig. 11. In the following, we use indicial notation where the Greek indices
(ff; fi) range over the values and a comma denotes a derivative with
respect to the following argument. The contour integral proposed by [17] in
the absence of an externally applied pressure is given by
I
where W is the strain energy density of the plate. This is defined as
We now introduce a weight function q 1 which is defined over the domain
of interest. Consider the simply connected curve
shown in Fig. 11. The function q 1 is defined to be sufficiently smooth in the
area A enclosed by C, and is given on the surfaces by:
ae 1 on \Gamma
on C
We then use this function to rewrite (41) as
I
Z
where we have used \Gamman i on \Gamma, and
on the crack faces. The last integral above vanishes for traction free crack
faces. Applying the divergence theorem to the closed integral, we then obtain
Z
A
Z
which is the equivalent domain form of the J 1 integral proposed by [17].
The measure number J 1 is domain independent and its magnitude is
equivalent to the energy release rate (37). Therefore, under pure mode I
loading, the moment intensity factor K I is given by
r
In more general mixed-mode conditions, the values K I ,K II and K III cannot
be separated so easily. While the comparison of J 1 to analytical values is
adequate to assess the effectiveness and qualities of the enrichment strategy,
crack growth laws are typically expressed in terms of the mixed-mode intensity
factors. In classical linear elasticity, the interaction integral approach
[19] has proven effective to extract mixed mode stress intensity factors. The
application of this method to plate fracture is currently under development
3.4 Enrichment of the MITC4 plate element
When discretizing the plate equations (26), some care must be taken to avoid
shear locking. As the plate becomes very thin (i.e. t ! 0), the following
relationship must be satisfied to keep the strain energy in the plate bounded:
In other words, the shear strain ffl s must vanish as t ! 0. Standard displacement
based elements, such as the four node isoparametric element, have
difficulty satisfying this constraint. The consequence is a structure which
exhibits an overly stiff response, often referred to as shear locking.
To discretize the plate displacements (16), we begin with the MITC4
element. To avoid shear locking, the MITC formulation modifies the approximation
for the section rotations / in the expression for the shear stiffness
(see [21]). In the following, we express this modification using the notation
~
I , where it is understood that only those expressions relating to the shear
components are modified.
The enriched discretization takes the form
I w I
c w
Kl G l (r; ')
(48a)
~
I / I
~
c /
Kl F l (r; ')
(48b)
where N I are the standard bilinear shape functions. In the above, we have
collapsed the sums over each crack tip into one for compactness.
The sets of near-tip functions G l and F l are derived from (40) in the
following fashion. We take G l to be only those functions in g l which are
proportional to r 3=2 . The set F l is taken to be equivalent to f l . In addition to
having four additional degrees of freedom for each displacement component,
this choice for G l and F l satisfies the following relation
such that a linear combination of the near-tip enrichment functions can satisfy
(47). We note that this relationship does not ensure that the enriched
formulation will be completely free of shear locking. However, the numerical
examples presented in the next section indicate that the above formulation
performs well for a wide range of plate thicknesses.
In this section we present several different numerical calculations. We first
examine some problems in two-dimensional elasticity, including a robustness
test and the simulation of crack growth. Then a benchmark and an additional
study are presented for Mindlin-Reissner plates.
4.1 Two-dimensional problems
We begin with a simple example of an edge crack to demonstrate the robustness
of the discretization scheme, and then present results for more complicated
geometries. In all of the following examples, the material is taken to be
isotropic with Young's modulus
and plane strain conditions are assumed. The calculation of the stress intensity
factors is performed with the domain form of the interaction integral,
and the maximum hoop stress law is used to govern crack growth ( see [19]
and [5]).
4.1.1 Robustness tests
Consider the geometry shown in Fig. 12: a plate of width w and height L
with an edge crack of length a, subjected to a far-field stress oe o . We analyze
the influence of the location of the crack with respect to the mesh on the K I
stress intensity factor when the position of the crack is perturbed by ffix in
the X direction and ffiy in the Y direction. The geometry is discretized with
a uniform mesh of 24x48 4-noded quadrilateral elements.
In this study, several different discretizations are obtained depending on
the position of the crack with respect to the mesh. Two cases are shown in
Fig. 13. In this investigation, we wish to examine the performance of the
modified tip function ~
R(x), and the accuracy of the formulation when it is
used in conjunction with the other near-tip functions as in (9).
The exact solution for this problem is given by [22]
a (50)
where C is a finite-geometry correction factor:
a
a
a
a
The numerical results normalized by the exact solution when only the
function ~
R(x) is used to model the near-tip region are given in Table 1.
Depending on the location of the crack tip, the total number of degrees of
freedom varies from 2483 to 2503. The results vary by approximately 4%
over all crack tip locations tested. When the near-tip functions are added,
the accuracy improves as shown in Table. 2. These results are consistent with
those reported in [6] in that the best results are obtained when the crack is
aligned with mesh boundaries. We note that the results are not as accurate
as when the exact asymptotic function
rsin(') is used, in which case the
error is less than 2% (see [6]).
4.1.2 Crack Growth from a Fillet
This example shows the growth of a crack from a fillet in a structural mem-
ber, and serves to illustrate how the present method can be used as an aid to
design against failure. The configuration to be studied is shown in Fig. 14,
with the actual domain modeled as indicated. The setup is taken from experimental
work found in [23]. In this example, we investigate the effect of
the thickness of the lower I-beam on crack growth. Only the limiting cases
for the bottom I-beam of a rigid constraint (very thick beam) and flexible
constraint (very thin beam) are considered. In addition, the welding residual
stresses between the member and the I-beam are neglected.
The structure is loaded with a traction of and the initial
crack length is taken to be a 5mm. The geometry is discretized with 8243
three-node triangular elements. To model a rigid constraint, the displacement
in the vertical direction is fixed along the entire bottom of the domain. A
flexible constraint is idealized by fixing the vertical displacement at both
ends of the bottom of the domain. For both sets of boundary conditions, an
additional degree of freedom is fixed to prevent a rigid body rotation.
For each load case, we simulate crack growth with a step size of \Deltaa =
5mm for a total of 14 steps. Fig. 15 shows the mesh in the vicinity of the
fillet and compares the crack paths for the cases of a thick I-beam (upper
crack) and a thin I-beam (lower crack). It is emphasized that the same mesh
is used throughout the simulation, and that no remeshing is required. As
new crack segments are added, additional enriched degrees of freedom are
generated for each new segment. The results shown are consistent with both
the experimental [23] and previous numerical results [24].
4.2 Plate examples
In this section, we present some examples using the enriched MITC4 plate
formulation developed in Section 3.4. We first examine the accuracy of the
method as a function of plate thicknesses for a benchmark problem, and
then present a more general example. Throughout this section, the material
properties are assumed to be isotropic with Young's modulus of
GPa, and Poisson's ratio
As a benchmark problem we consider a through crack in an infinite plate
subjected to a far-field moment M o . The crack is oriented at an angle fi
with respect to the x 1 axis as shown in Fig. 16. Recently, very accurate
calculations were carried out by [25] for various plate thicknesses for the case
when . In this case, the loading is purely mode I, and the domain form
of the J-integral for plates (45) is used in conjunction with (46) to determine
the moment intensity factor K I . In the finite element model, only one-half
of a square plate is modeled, with symmetry conditions along the x 2 axis.
To approximate the infinite plate, the plate width w is taken to be 10 times
the half crack length a. The crack length for all of the results presented in
this section is taken to be
Fig. 17 shows the normalized K I for four discretizations, two standard and
two enriched. The lower curve corresponds to a non-enriched formulation,
and the values for K I are within 5% of the exact for the entire range of
plate thicknesses t. These values are improved when the mesh is refined for
a total of 2463 degrees of freedom as shown. We observe that the enriched
solution with only 755 degrees of freedom is as accurate as the solution with
2463 degrees of freedom without enrichment. The last curve for the enriched
case with 3087 degrees of freedom exhibits less than 1% error. The enriched
solutions show good correlation with the analytical solution for the full range
of plate thicknesses tested.
As a last example, moment intensity factors are calculated for a finite
plate as a function of crack length for various plate thicknesses. The geometry
of the plate is taken to be the same as the previous example, and the
results are compared to those given in [26]. In this study, the mesh does
not model the crack discontinuity; the jump in the rotations and transverse
displacement is created entirely with enrichment. Table 3 gives the results
for four different plate width to thickness ratios for the case when the plate
is modeled with 1424 MITC4 elements. These results show excellent correlation
for the cases when w=8, in which the maximum
error is 1.2%. For the remaining cases the maximum difference between the
numerical solutions and those given in [26] is 9.4%. We note, however, that
the reference [26] is not as current as [25]. In the latter, the moment intensity
factors are shown to be significantly greater than the classical results as
the thickness t ! 0. The results shown in Table 3 are consistent with these
findings.

Summary

A method of constructing finite element approximations with enrichment
functions was presented which allows for the simulation of evolving discontinuities
in a straightforward fashion. The specific examples of cracks and
crack growth in two-dimensional elasticity and Mindlin-Reissner plate theory
were examined. By incorporating the appropriate asymptotic near-tip
fields, accurate moment and stress intensity factors were obtained for coarse
meshes. A new near-tip function was also developed to remove the need for a
mapping in the case of kinked cracks. The methodology for the construction
of the discrete approximation from the interaction of the crack geometry and
the mesh was provided, and numerical tests served to illustrate the algo-
rithm's robustness. Additional numerical studies for Mindlin-Reissner plates
demonstrated the extent to which stress intensity factors can be calculated
accurately for a wide range of plate thicknesses.
The present method has a lot of potential to extend the finite element
method for the modeling of evolving interfaces and free surfaces. A key feature
of the enrichment in conjunction with numerical integration is the capability
of modeling geometrical features which are independent of the mesh
topology. As was shown in this paper, several different crack configurations
can be considered for a single mesh of a component, simply by changing
the enrichment scheme according to the crack geometry. Future work will
focus on the application of the method to three-dimensional and dynamic
fracture, as well as other areas of mechanics in which moving interfaces are
of importance.

Acknowledgements

The support of the Office of Naval Research and Army Research Office, to
Northwestern University, is gratefully acknowledged. The authors are grateful
for the support provided by the DOE Computational Science Graduate
Fellowship program, to John Dolbow.



--R

A hybrid-element approach to crack problems in plane elasticity
Modeling mixed-mode dynamic crack propagation using finite elements: Theory and applications

Modelling strong discontinuities in solid mechanics via strain softening constitutive equations.
Elastic crack growth in finite elements with minimal remeshing.
A finite element method for crack growth without remeshing.
An extended finite element method with discontinuous enrichment for applied mechanics.

Multiple scale finite element methods.

Meshless methods: An overview and recent developments.
Enriched methods for singular fields
On the bending of an elastic plate containing a crack.
Mechanics of fracture 3: Plates and shells with cracks.

Elliptic Problems in Nonsmooth Domains.
Computation of stress intensity factors for plate bending via a path-independent integral
Crack tip and associated domain integrals from momentum and energy balance.

Modeling fracture in Mindlin- Reissner plates with the extended finite element method
Displacement and stress convergence of our MITC plate bending elements.
Fracture Mechanics.
Morphological aspects of fatigue crack propagation.
The Element-free Galerkin Method for Fatigue and Quasi-Static Fracture
Bending of a thin Reissner plate with a through crack.
Internal and edge cracks in a plate of finite width under bending.
--TR

--CTR
Kenjiro Performance assessment of generalized elements in the finite cover method, Finite Elements in Analysis and Design, v.41 n.2, p.111-132, November 2004
L. B. Tran , H. S. Udaykumar, A particle-level set-based sharp interface cartesian grid method for impact, penetration, and void collapse, Journal of Computational Physics, v.193 n.2, p.469-510, 20 January 2004
H. S. Udaykumar , L. Tran , D. M. Belk , K. J. Vanden, An Eulerian method for computation of multimaterial impact with ENO shock-capturing and sharp interfaces, Journal of Computational Physics, v.186 n.1, p.136-177, 20 March
partition-of-unity;fracture;discontinuous enrichment