- Email: [email protected]
A nonlocal elastic damage theory: Mesh-insensitivity under strain softening
Compvfers& S~rucfuresVol. 48, No. 3. pp. 412-422, 1993 Printed in Great Britain.
OMS-7949/93s6.M)+ 0.00 0
1993 Pergamott msa Ltd
A NONLOCAL ELASTIC DAMAGE THEORY: MESH-INSENSITIVITY UNDER STRAIN SOFTENING H. MURAKAMI,~D. M. KENDALL?and K. C. VALANIS~ tDe.partment of Applied Mechanics and Engineering Sciences, University of California at San Diego, La Jolla, CA 92093-0411,U.S.A. $ENDOCHRONICS Inc., Vancouver, WA 98665, U.S.A. (Received 14 ApriI 1992)
AIrstrati-A nonlocal damage theory was recently proposed by Valanis (J. Appl. Mech., ASME 58, 311-316, 1991). In the theory the basic conservation equations take a conventional local form and only the damage evolution law becomes nonlocal. In contrast to local theories, the resulting initial value problem of an elastic damaging solid retains its hyperbolicity even under strain softening. Thus, mesh sensitivity of the numerical solution observed in local damage theories is eliminated. The domain of damage localization, controlled by the heterogeneous microstructural interaction, converges to a Snite material region as the mesh size decreases to zero. This is in contrast to local models which predict an infinitely thin fracture domain of vanishing volume. Nonlocal fracture models based on the theory can be easily implemented into existing finite element wave codes by merely adding the constitutive subroutine for damage evolution. Convergence and mesh insensitivity are demonstrated by solving initial value problems using an updated Lagrangian finite element wave code.
1. INTRODUCTION Strain softening behavior has been observed in heterogeneous materials, such as concrete and rocks. Sandler and Wright [l] reported troublesome mesh sensitivity in finite element analyses of initial value problems with rate-independent strain softening materials. This mesh sensitivity is due to deformation trapping[2] which is caused by the loss of hyperbolicity of the local constitutive model. Thus, in order to avoid mesh sensitivity in finite element analyses, a rate-independent, strain softening constitutive model must retain hype&&city. An excellent review on strain softening based upon existing experimental data was compiled by Read and Hegemier [3]. A nonlocal damage theory recently proposed by Valanis [4] ensures that the ‘wave’ equation retains hyperbolicity even under strain softening. The objective of the present paper is to investigate numerically the mesh insensitivity of the model under strain softening, by implementing it into a finite element wave code. Convergence and mesh insensitivity are demonstrated numerically by solving the problems considered by Sandler and Wright [l] and Baiant and Belytschko [5]. In the former, a semi-infinite rod is subjected to a constant end velocity which generates strain softening in the first element. In the latter, a finite rod is subjected to constant end velocities, which generate tensile waves propagating inward. Strain softening occurs at the center element when the two waves meet. Local damage theories tend to localize in a region of vanishing volume, and converge to zero energy dissipation, which is physically
unrealistic [6]. On the contrary, the nonlocal theory localizes into a finite region whose size is defined by a characteristic length of nonlocal process zone that describes the heterogeneous, microstructural interaction. A 6nite dissipation energy is predicted as the mesh size tends asymptotically to zero. 2. SUMMARY
OF THE NONLOCAL THEORY
EIMXIC
DAMAGE
The basic premise of Valanis’ nonlocal elastic damage theory [4] has evolved from his local damage theory [7J. In these theories, damage is represented by a positive semi-definite second order tensor 4,,, called the ‘integrity tensor’, which appears expicitly in the expression of free energy and constitutive equations. In this section the nonlocal elastic damage theory is briefly summarized with respect to a rectangular Cartesian coordinate system. The free energy density $ is given by Iz
where 1 and p and the Lame constants, and crj represents the strain tensor. For undamaged materials, the integrity tensor takes the form +,j = a,, where 6,j is the Kronecker delta. The above form of free energy leads to a stress-strain relation for elastic damaging materials described by
(2)
415
416
H.
MURAKAMI et al.
where u,, is the stress tensor. The damage force tensor Q ,,, which is thermodynamically conjugate to the integrity tensor, is given by
equality, may be found in [4] and [A. In what follows, the theory is specialized to uniaxial stress loading, and is employed to describe strain softening materials in one-dimensional, initial value problems. 3. ONEDIMENSIONAL
With respect to the damage evolution, Valanis [7] observed that the direction of the damage increment -d& is dictated by the principal directions, .T, a = 1-3, of the strain increment dc, defined by the eigen-equation dc..n? 1, J = dc”n”I,
no sum on a,
(4)
where dt” are the principal strain increments. The components of the damage force and the strain in the direction ns are given, respectively, by Qi = n;n;Qij
and
au
no sum on a,
dt” = max 0, ( JS
k(x,, x() de; dx; dx; dx; >
if dt: > 0, c; > 0, Q:: < 0, = 0
otherwise,
(7)
where max(.) implies the maximum of the arguments and imposes the inequality dc” > 0. In the above, the normal distribution function given below is employed
(8)
The above nonlocal damage evolution law reduces to the local form [7] when the nonlocal process zone characterized by l/,/b tends to zero, i.e., when the distribution function becomes the Dirac delta function k(x,,x()=k,6(xi-xj).
d=Ei+&
E=E&’
(lhb)
i=-au ax’
(12)
where Q is the Cauchy stress, E,, and E are the intact and current elastic moduli, respectively, c is the logarithmic strain, and 4 is the integrity parameter which describes the stiffness degradation and lies in the range: 0 < 4 < 1; the undamaged state is represented by 4(O) = 1. In eqn (11) it is noted that due to the lack of spin in this one-dimensional problem, any objective stress rates reduce to the material time derivative of the Cauchy stress. (c) The damage evolution law 4 = -E,,q5&,
(13)
where 5 (> 0) is the damage coordinate and evolves according to
b k(x,, xi) = ko 3 rc exp - {(xi -xi)* J +(x2-x;)*+(xg--xi)*}.
where p is mass density; (b) The constitutive relation
(6)
where dt” is the increment of the damage coordinate in the ny direction. The evolution of the damage coordinate is described by the following nonlocal expression
(10)
;i;;=p",
(5a, b)
The damage evolves linearly with the damage force according to d4” = - Qt d<“,
VALUE PROBLEMS
Consider one-dimensional wave propagation in an elastic-fracturing solid in the axial direction denoted by coordinate x. The axial deformation is described by a map, x = x(X, t), between the material (Lagrangian) coordinate, X, and the spatial coordinate, x, at time I (see Fig. 1). Axial Cauchy stress, strain rate, and velocity are denoted by Q, i and v (= 2), respectively. In what follows, overdot implies material time differentation. The basic equations of the damage theory pertaining to the problem of wave propagation in one dimension follow. (a) The equation of motion
ci = nyn;cij, no sum on a.
INITIAL
(9)
Details of the above theory which is developed within the framework of irreversible thermodynamics by imposing the Clausius-Duhem dissipation in-
[(x, t) = max 0, (s
k(x - x’)i(x’, t) dx’ >
for .5(x, t) > 0, L(x, t) > 0, = 0
otherwise
Undeformed Configuration X
IX
(14) Deformed
-D
x,=x(X,.t)
Fig. 1. A deformation map.
Configuration
417
A nonlocal elastic damage theory in which the integration is performed over the onedimensional region, (x,,, x, ), occupied by the solid. The kernel function, k(x) is given by
k
z k(x) = k,
J$
exp(-bx*),
(15)
where k,, and b are material constants; in particular, the constant l/,/b represents a characteristic length of nonlocal interaction for the evolution of the damage coordinate t. The solution of the above equations necessitates the specification of (a) appropriate boundary data at x = x0 and x, and (b) initial conditions on x and v at time f = 0. It was shown in [4] that the equations of motion expressed by the displacement become unconditionaIly hyperbolic due to the nonlocal damage evolution, and the wave speed is given in terms of the secant modulus E, instead of the tangent modulus (see Appendix A). A corresponding local damage theory may be obtained from the above by setting the kernel function equal to the Dirac delta function k(x) = k,&x).
(16)
For both nonlocal and local damage theories, the energy dissipation W, of the solid, with an initial cross-sectional area A0 , becomes
Qd dX dt,
W,=A,
(17)
ss where the damage force Q is defined as
(18) The initial boundary value problem defined by eqns (lO)-(lS) for the nonlocal damage theory and eqns (10x14) with eqn (16) for the corresponding local theory is solved numerically by employing an explicit, updated Lagrangian finite element method. 4. UNIAXIAL
STRJBS-STRAIN
RELATION
In this section, the parameters of the above elastic damage model are related to quantities which can be measured by uniaxial tension and compression tests. It is noted, however, that the above ideal model was introduced as the simplest constitutive model which exhibits strain softening for the investigation of meshsensitivity in initial value problems. When the specimen length is large enough compared with the nonlocal reference length l/,/b, the local elastic damage theory becomes applicable, in which case eqn (15) can be employed in eqn (14). Therefore, parameters of the damage theory, except 6, can be determined by monotonic uniaxial tension and com-
b
“0
200
loo e
300
400
(PI
Fig. 2. Cyclic tensile stress-strain relation. pression tests. The initial Young’s modulus E,, can be determined by the slope of the stress-strain curve at t = 0. The parameter k0 in eqns (14) and (15) is related to the critical strain t, in tension [4], in terms of the following c = E,,c exp(-fk,E,,c’),
k0 = 1/(2E,t,‘).
(19)
According to eqn (14) damage evolution in compression is suppressed. Therefore, loading response in compression is perfectly elastic for the present model; if necessary for the modeling of structural elements, the damage evolution in compression can be introduced by modifying the damage evolution law, eqn (14). In the sequel, example calculations are presented for the material with E,, = 34.5 GPa, p = 2400 kg/my, and ce = 0.12 x 10-3. The tensile loading, unloading and reloading responses of the elastic damaging material are presented in Fig. 2. The unloading and reloading responses in tension are elastic, but with reduced elastic modulus. The determination of the nonlocal characteristic length l/Jb may require a suitable set of fracture tests. Since the objective of this paper is to investigate finite element mesh-sensitivity in initial value problems, it is assumed that the parameter b is given. 5.
FINITE ELEMENT
IMPLEMENTATION
The governing equations of motion, (IO), and stress-strain relations, eqns (11) and (12) are in local form. Therefore, updated Lagrangian formulation based upon the principle of virtual work holds. Consequently, the new damage model can be easily incorporated into existing finite element codes by simply adding a subroutine for the nonlocal damage evolution, eqns (13) and (14). An explicit, updated Lagrangian code was developed by using two-node elements. Since the implementation of updated Lagrangian codes is well documented [g-11] (see Appendix B), only the implementation of the damage evolution law (14) is briefly explained. For a typical integration point at x, where i and c are both positive, the increment of the damage coordinate A{ under time increment At is
H. MURAKAMI et al.
418
obtained from eqn (14) by using two-point Gauss quadrature over the elements which cover the support of k (x -x’)
-0.0 VI 2
M(x) = C i k(x - x’(tl,))A+‘(~J) c 6 -I
-0.1
_ _
O.Im
_
0.02s
-
0.0125m
0.05 m m
, (20)
$ 1
where the first summation is over the elements on the support, the second summation is over the Gauss integration points, L, is the element length, and 200 Time(p)
q*= -a,=‘.
(21) J5
This scheme, employed by the authors [12], is economical and can be easily generalized for two- and three-dimensional regions. However, the above numerical integration becomes inaccurate under extreme element stretch in the final stage of tensile failure. Therefore, an alternative numerical scheme with added computational time is presented in what follows. In the present numerical investigation, the nonlocal integration was evaluated analytically on the premise that Ae is known only at the integration points. For one point integration, A6 is constant over the element, while for two point integration, de is constant over the half region of the element. Denoting the range (xi, xi) at the pth integration point, eqn (14) yields
At(x)=C
c
i
9
Mx’(~J)ko
p=I
J XIt X 5
exp{ -b(x
-x’)*)
dx’
(22)
x:
which, upon evaluation of the integral, results in the expression
-erf(&(x
- xi)}].
(23)
The second summation in eqns (22) and (23) is over the integration points, and erf{*} denotes the error function whose value may be found numerically by a
c--x 0 Fig. 3(a). Problem 1: tensile loading of a semi-infinite bar.
+
-L
01
c
L
Fig. 3(b). Problem 2: symmetric tensile loading of a bar of
finite length.
Fig. 4. Velocity history at X = 0.4 m for c = 0.32m/sec computed by using different meshes.
rational approximation (see Appendix C) shown, for example, by Abramowitz and Stegun [13]. 6.
NUMERICAL RESULTS
In the numerical studies, attention was focused on the two problems, illustrated in Figs 3(a and b), involving a semi-infinite bar (Problem 1) and a finite bar (Problem 2). These problems were selected because they offer a means of studying the behavior of the damage models when strain softening occurs at a surface (Problem l), or in the interior of a material (Problem 2). In the subsequent finite element (FE) analyses uniform, initial meshes are employed to investigate the mesh sensitivity. In Problem 1, shown in Fig. 3(a), a particle velocity c is suddenly applied at time t = 0 to the end of a semi-infinite, initially quiescent bar, producing a tensile wave that travels to the right in the bar. In the FE approximation, the semi-infinite region was approximated by a finite region of length L = 1.6 m. The boundary conditions imposed in the FE analysis are u(x(0, t),
t) = -cH(t),
u(x(L, t), t) = 0
(24a, b)
where H(t) (= 1 for t > 0 and 0 for t < 0) is Heaviside’s step function. The nonlocal characteristic length l/Jb = 0.1 m was employed in the calculation. A constant time step defined by the Courant number computed by the intact elastic modulus is used. The calculation was terminated when the wave front hit the fictitious boundary X = L. This problem in compression, instead of in tension, was originally proposed by Sandler and Wright [l] as a means of illustrating the difficulties pertaining to loss of hyperbolicity when local strain-softening models are employed. The critical velocity to cause strain softening is approximately 0.326m/sec. First, a velocity of c = 0.32 m/set was chosen to ensure that the material does not enter the strain softening region. The velocity history at X = 0.4 m, computed by employing the local damage theory, is shown in Fig. 4 for different initial mesh size AX = 0.1 m, 0.05 m, 0.025 m, and 0.0125 m, respectively. The results show
419
A nonlocal elastic damage theory that as mesh size decreases the numerical solution converges to a unique ~lution; this was also true for the nonlocal damage theory. Next, the calcufation was repeated with a larger velocity, c = 0.6 m/set so that the material experiences strain softening. Figures S(a and b) show the calculated velocity histories at X = 0.4 m, by employing the local and nonlocal theories, respectively. Reference to Fig. S(a) reveals that the solution based upon the local theory shows no evidence of converging to a physically realistic solution as the mesh size AX is reduced; instead it appears to be approaching a vertical line. The result of the nonlocai theory, shown in Fig. S(b), on the other hand, clearly show convergence to a unique, physically realistic solution, Note that the solution obtained with the local model for AX = 0.1 m is quite cbse to the converged solution for the nonhz-cal modek The corresponding stress histories at the integration point of the right side element adjacent to X = 0.4 m are shown in Figs 6(a and b) for the local and the nonlocal theories, respectively. With respect to convergence, the results illustrate the same trend as the velocity histories. The dissipated energy is shown in Fig. 7 for A0 = I m2. As shown in this figure, the local solution converges to the situation in which there is zero energy dissipation
Fig. 6(b). Stress history at the right integration point adjacent to X = 0.4 m, predicted by the nonlocal theory for c = 0.6 m/w, and computed by using different meshes.
0.0 n
2 -
A
-0.1
.E
0
J
fig. 6(a)* Stress history at the right integration point adjacent to X ~5:0.4m, predicted by the local theory for c = 0.6 m/xc, and computed by using differeat meshes.
-0.2
5.3
0
loo
_
d.Im
_ _
0.0sIll 0.025m
-
0.012cm
200
300
400
TimeIp5)
Fig. 5(a). Velocity history at X = 0.4 m, predicted by the locak theory for c = Wm/sec, and computed by using different meshes.
as the bar edge softens and eventually breaks. This situation is not physically realistic for materials with heterogeneous ~c~s~cture. The nonlocal solution, shown in Figs 5(b) and 6(b), on the other hand, converges to a solution where softening occurs in a finite edge region characterized by the nonlocal characteristic length l/,/b = 0.1 m and finite dissipation energy is recorded as the mesh size tends to zero. Figure 7 atso illustrates, for comparison, the dissipation energy for the case, shown in Fig. 4, with no strain softening.
0.0 t.3
-c -a. .z
49.1
4
9
5.2
a.3
0
100
200
300
400
0 v.00
Fig. S(b). V&city history at X = OS4m, predicted by the
nonlocal theory for c = 0.6 m/set, and computed by using different meshes.
0.05
5.10
0.15
0.20
AX Cm) Fig. 7. Dissimtion energy veraus mesh size (Problem 1).
420
H. MUMKAMer al.
Problem 2, illustrated in Fig. 3(b), consists of the symmetrical tensile loading of a bar of finite length 2L. The initially quiescent bar is suddenly loaded by applying end velocities, - c and c, at time f = 0. The boundary conditions are u(x(-
2-
1-
L,t),t)=-cH(t), v(x(L,t),t)= d(t). (25a,b)
The velocity c is selected so that (a) the maximum stress in the initial waves that propagate from the ends toward the center at x = 0 is less than that necessary to cause strain softening and (b) the maximum stress that results when the two incoming waves collide at x = 0 is large enough to cause strain softening. In the present study, we set c = 0.3 m/set and took the length of the bar 2L to be 1.6 m. This problem was originally studied by Baiant and Belytschko [5], where an exact solution was obtained for the case in which the material behavior was described by a local strain softening model. The calculated velocity histories at a quarter point, X = 0.4 m, by the local theory and the nonlocal theory are shown in Figs 8(a and b), respectively. The corresponding stress histories at the closest integration point to the right of X = 0.4 m are illustrated
-1
0
100
200
300
Time ( p
400
)
Fig. 9(a). Stress history at the right integration point adjacent to A’ = 0.4m, predicted by the local theory for c = 0.3 m/set, and computed by using different meshes.
2
_
x
cn
-I
0
loo
O.lm
_
0.05
_
O.OZSm
-
0.012Sm
200
111
300
4co
Time ( ps )
0.2
Fig. 9(b). Stress history at the right integration point adjacent to X = 0.4 m, predicted by the nonlocal theory for c = 0.3 m/xc, and computed by using different meshes.
-0.6 0
IO0
200
300
400
Time(p)
Fig. 8(a). Velocity history at X = 0.4 m, predicted by the local theory for c = 0.3 m/set, and computed by using different meshes. 0.2
0.0
_ _
oosm
-
00125 111
in Fig. 9(a) for the local theory and in Fig. 9(b) for the nonlocal theory. As the two tensile waves meet the specimen breaks into two pieces, as shown in the stress histories in Fig. 9. Both the local and the nonlocal theories predict that splitting. The local theory predicts vertical stress drop reflecting localized failure in an infinitely thin layer. On the contrary, the nonlocal theory predicts the splitting failure in a finite volume as shown by the slanted stress drop in Fig. 9(b). This scenario is better explained by showing the dissipated energy in Fig. 10. The local theory converges to the finite dissipation energy (reference data)
Ol125m
-0.2
-0.6 0
loo
200 Time (pi
3ocl
400
)
Fig. 8(b). Velocity history at X = 0.4 m, predicted by the nonlocal theory for c = 0.3 mjsec, and computed by using different meshes.
Fig. 10. Dissipation energy versus mesh size (Problem 2).
A nonlocal elastic damage theory
which is predicted when the bar does not experience any strain softening (c = 0.16 m/set). Finally, it is informative to mention the additional computational time required to perform the calculation using the nonlocal damage evolution, eqns (14) and (15), relative to the local damage evolution, eqns (14) and (16). For the fine mesh calculations, AX = 0.025 and 0.0125 m, the nonlocal analyses typically take 30-50 times more CPU time than the local analyses. This order of magnitude difference is due to the error function calculation (C2). Therefore, it is expected that by creating a look-up table for the error function and eliminating the repeated evaluation of the error functions at adjacent elements, the additional computational effort required for the nonlocal damage evolution will reduce to less than a factor of ten. 7. CONCLUSIONS
This study consisted of applying the nonlocal elastic damage theory to uniaxial stress problems and examining the performance of the constitutive model in wave propagation studies, especially under strain softening. From the results of this study, the following conclusions are drawn: 1. The model is computationally feasible for use in conjunction with finite element wave codes. 2. In wave propagation problems involving strain softening, the model leads to (physically plausible) unique solutions. This was dramatically illustrated in the numerical analysis of the problem involving the tensile loading of a semi-infinite bar (Figs 5 and 6); as the mesh size was reduced, the nonlocal solution converged to a physically realistic solution, while the local solution did not. 3. It is recommended that this damage theory be expanded to incorporate plastic damaging process to simulate the mechanical response of rocks, concrete, and ceramics.
421
5. Z. P. Baiant and T. B. Belytschko, Wave propagation in a strain-softening bar: exact solution. ASCE I. Engng
111, 381-398 (1984).
Mech.
6. Z. P. Baiant, T. B. Belytschko and T:-P. Chang, Continuum theory for strain-softening. ASCE I. Engng Mech. 110, 1666-1692 (1984). 7. K. C. Valanis, A theory of damage in brittle materials. Engng Fract. Mech. 36, 403-416 (1990). 8. K. J. Bathe, Finite Element Procedures in Engineering Analysis. Prentice-Hall, Englewood Cliffs, NJ (1982). 9. T. J. R. Hughes, The Finite Element Method, Linear Static and Dynamic Finite Element Analysis. PrenticeHall, Englewood Cliffs, NJ (1987). 10. J. 0. Hallquist, Theoretical manual for DYNA3D. Report UCID-19401, University of California, Lawrence Livermore National Laboratory, Livermore, CA (1982). 11. G. L. Goodreau and J. 0. Hallquist, Recent development in large-scale finite element Lagrangian hydrocode technology. Comput. Meth. Appi. Mech. Engng 33, 725-757
(1982).
12. H. Murakami, D. M. Kendall and K. C. Valanis, Numerical validation of a nonlocal damage theory. In Constitutive Laws for Engineering Materials, Recent Advances and Industrial and Infrastructure Applications
(Edited by C. S. Desai, E. Krempl, G. Frantziskonis and H. Saadatmanesh), pp. 743-746. ASME Press, New York (1991). 13. M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions, pp. 297-329. Dover, New York (1964).
APPENDIX
A: UNCONDITIONAL
HYPERBOLICITY
The type of initial value problem defined by eqns (lO)-(lS) can be examined by investigating that of the corresponding rate problem. Upon monotonic loading, eqns (llb), (13) and (14) yield R = -2EE,c*
k(x - x’)i(x’,
r)dx’,
s
(Al)
where the integration is performed over (x,,, x,). By substituting the above equation into eqn (1 la), a rate constitutive relation is obtained in the form ti = E(i -2E&j
k(x - x’)C(x’, t) dx’ )
(A2)
From eqn (10) the equation of motion becomes, in rate form Acknowledgement-The
authors express their gratitude
ad
to Dr H. E. Read of S-CUBED, a division of Maxwell
ax=pa,
Laboratories, La Jolla, California, for his helpful advice and comments. By substituting obtains
a=d.
eqn (A2) into the above equation,
one
REFERENCES
1. I. S. Sandler and J. P. Wright, Strain softening. In Theoretical Foundation for Large-scale Compuraiions for Nonlinear Materials (Edited by S. Nemat-Nasser,
R. J. Asaro and G. A. Hegemier), pp. 285-296. Martinus Nijhoff, Dordrecht (1984). 2. F. H. Wu and L. B. Freund, Deformation trapping due to thermoplastic instability in one-dimensional wave propagation. J. Mech. Phys. Solirls 32, 119-132 (1984). 3. H. E. Read and G. A. Hegemier, Strain softening of rock, soil and concrete-a review article. Mech. Materials 3, 271-294 (1984). 4. K. C. Valanis, A global damage theory and the hyperbolicity of the wave problem. ASME I. Appl. Mech. 58, 311-316 (1991).
E{;-ZE,,r’Jk*(x-x’)t(x’,i)dx’
-6E,c’;
k(x - x’)((x’,
t) dx’
s +~{G2E0r3~k(x-x’)c(x’,~)dx’}=pd,
, 1 (A4)
where k*(x) = dk(x)/dx. In terms of the velocity field v(x, r), the wave equation (A4) is of the form E&p$+E;+F
(A5)
H. MURAKAMIet cl.
422
where F is a spatial function for au/ax over the material domain D and depends, parametrically, on the strain field 6, where
t sa0
c=
dt.
F<=j;[-&Nt$N2]6dx,
,ax
Given the functional F{ .} in D at I, the character of eqn (A5) in the infinitesimal nonvanishina time interval At, is determined by the higher order derivative terms’ on the left-hand side of eqn (A5). Therefore, eqn (A5) is hyperbolic if E, the secunt modulus, and p, the density, are positive. This is the case even under conditions of material softening. Furthermore the wave speed C,, is given in terms of E and p by eqn (A7)
c,=
The evaluation of the integrals of eqn (Bl) on a typical element (e) yields the internal force vector F, and the diagonalized mass matrix M,
M,=-
P-G
1
2 [0
(B5) 0
1
1’
E
where the superscript T denotes transpose and L=(=x2 - x’) is the length of the element (e). By assembling the nodal coordinate vector X, the nodal velocity vector V, the nodal acceleration vector A ( = V), and the internal force vector F with respect to the global degrees of freedom, one obtains the following equations of motion atatimestepf=t”
P
MA” = I”’ - F”,
J -
and is always real.
APPENDIX B: EXPLICIT, UPDATED LAGRANGIAN FINITE ELEMENT SCHEME
The principle of virtual velocity applied to the initial value problem defined on the current configuration (x0, x,) yields
+ Wx,, t)T(x,, 1)-WA,,
r)W,, lb
@I)
where 6u denotes the virtual velocity, x0 = x (X0, I), xi = x (X,, I), and T is the x component of the traction vector acting on x = x,, or x, . The finite element discretixation of the above weak equation has been explained by many authors, for example, Bathe [8] and Hughes [9]. In what follows, the outline of the procedure is briefly summarized. The one-dimensional domain (x0,x,) is discretized into two-node elements
(W
Vn+1/2=Vn-1/2+f(Atn+Atn+I)An
(B9)
whereAt,(=t,-t,_,)andAr,+,(=r,+,-r,)aretimeincrements: APPEND=
c: ERROR FUNCTION
The error function is defined in eqn (Cl)
0 (&-a9
exp( - t2) dt.
cc11
r=l
where N, is the total number of the elements, and x$, k = 1.2. is the k th nodal coordinates denoted in vector form as X,1 In a typical isoparametric element (e), the coordinate and the velocity are defined on the standard region - 1 C r Q I
as
It can be evaluated numerically by employing the following rational approximation [13]
xexp(-n2)+t(x) r = &,
where af (k = 1,2) are the elemental nodal-velocities, denoted in vector form as V,, and NL(r) (k = 1,2) are shape functions defined on the standard region as follows: N’({) = f(1 - 5)
(B7)
where M is the diagonal, global mass matrix, and P” is the global, external force vector due to the traction T in eqn (Al). For the class of problems treated in this paper, velocity boundary conditions are prescribed. Therefore, the last two terms on the right-hand side of eqn (Al) vanish, since the virtual, boundary velocities are zero, &(x0, t) = du(x,,f)=O; this yields P”=O. The explicit, central difference time integration scheme, implemented, for example, in DYNA2D and DYNA3D [lo, 111, is summarized as follows: A”=M-I@‘“-F”)
(x,7x,)=
(W
N?(T) = f(l + 0.
Wa, b)
k(x)] c 1.5 x 10-7,
where the constants are defined as follows: p = 0.3275911, a, = 0.254829592, n2= -0.284496736,
oS = 1.421413741,
a. = - 1.453152027. or = 1.061405429.
(C3) \ I
OMS-7949/93s6.M)+ 0.00 0
1993 Pergamott msa Ltd
A NONLOCAL ELASTIC DAMAGE THEORY: MESH-INSENSITIVITY UNDER STRAIN SOFTENING H. MURAKAMI,~D. M. KENDALL?and K. C. VALANIS~ tDe.partment of Applied Mechanics and Engineering Sciences, University of California at San Diego, La Jolla, CA 92093-0411,U.S.A. $ENDOCHRONICS Inc., Vancouver, WA 98665, U.S.A. (Received 14 ApriI 1992)
AIrstrati-A nonlocal damage theory was recently proposed by Valanis (J. Appl. Mech., ASME 58, 311-316, 1991). In the theory the basic conservation equations take a conventional local form and only the damage evolution law becomes nonlocal. In contrast to local theories, the resulting initial value problem of an elastic damaging solid retains its hyperbolicity even under strain softening. Thus, mesh sensitivity of the numerical solution observed in local damage theories is eliminated. The domain of damage localization, controlled by the heterogeneous microstructural interaction, converges to a Snite material region as the mesh size decreases to zero. This is in contrast to local models which predict an infinitely thin fracture domain of vanishing volume. Nonlocal fracture models based on the theory can be easily implemented into existing finite element wave codes by merely adding the constitutive subroutine for damage evolution. Convergence and mesh insensitivity are demonstrated by solving initial value problems using an updated Lagrangian finite element wave code.
1. INTRODUCTION Strain softening behavior has been observed in heterogeneous materials, such as concrete and rocks. Sandler and Wright [l] reported troublesome mesh sensitivity in finite element analyses of initial value problems with rate-independent strain softening materials. This mesh sensitivity is due to deformation trapping[2] which is caused by the loss of hyperbolicity of the local constitutive model. Thus, in order to avoid mesh sensitivity in finite element analyses, a rate-independent, strain softening constitutive model must retain hype&&city. An excellent review on strain softening based upon existing experimental data was compiled by Read and Hegemier [3]. A nonlocal damage theory recently proposed by Valanis [4] ensures that the ‘wave’ equation retains hyperbolicity even under strain softening. The objective of the present paper is to investigate numerically the mesh insensitivity of the model under strain softening, by implementing it into a finite element wave code. Convergence and mesh insensitivity are demonstrated numerically by solving the problems considered by Sandler and Wright [l] and Baiant and Belytschko [5]. In the former, a semi-infinite rod is subjected to a constant end velocity which generates strain softening in the first element. In the latter, a finite rod is subjected to constant end velocities, which generate tensile waves propagating inward. Strain softening occurs at the center element when the two waves meet. Local damage theories tend to localize in a region of vanishing volume, and converge to zero energy dissipation, which is physically
unrealistic [6]. On the contrary, the nonlocal theory localizes into a finite region whose size is defined by a characteristic length of nonlocal process zone that describes the heterogeneous, microstructural interaction. A 6nite dissipation energy is predicted as the mesh size tends asymptotically to zero. 2. SUMMARY
OF THE NONLOCAL THEORY
EIMXIC
DAMAGE
The basic premise of Valanis’ nonlocal elastic damage theory [4] has evolved from his local damage theory [7J. In these theories, damage is represented by a positive semi-definite second order tensor 4,,, called the ‘integrity tensor’, which appears expicitly in the expression of free energy and constitutive equations. In this section the nonlocal elastic damage theory is briefly summarized with respect to a rectangular Cartesian coordinate system. The free energy density $ is given by Iz
where 1 and p and the Lame constants, and crj represents the strain tensor. For undamaged materials, the integrity tensor takes the form +,j = a,, where 6,j is the Kronecker delta. The above form of free energy leads to a stress-strain relation for elastic damaging materials described by
(2)
415
416
H.
MURAKAMI et al.
where u,, is the stress tensor. The damage force tensor Q ,,, which is thermodynamically conjugate to the integrity tensor, is given by
equality, may be found in [4] and [A. In what follows, the theory is specialized to uniaxial stress loading, and is employed to describe strain softening materials in one-dimensional, initial value problems. 3. ONEDIMENSIONAL
With respect to the damage evolution, Valanis [7] observed that the direction of the damage increment -d& is dictated by the principal directions, .T, a = 1-3, of the strain increment dc, defined by the eigen-equation dc..n? 1, J = dc”n”I,
no sum on a,
(4)
where dt” are the principal strain increments. The components of the damage force and the strain in the direction ns are given, respectively, by Qi = n;n;Qij
and
au
no sum on a,
dt” = max 0, ( JS
k(x,, x() de; dx; dx; dx; >
if dt: > 0, c; > 0, Q:: < 0, = 0
otherwise,
(7)
where max(.) implies the maximum of the arguments and imposes the inequality dc” > 0. In the above, the normal distribution function given below is employed
(8)
The above nonlocal damage evolution law reduces to the local form [7] when the nonlocal process zone characterized by l/,/b tends to zero, i.e., when the distribution function becomes the Dirac delta function k(x,,x()=k,6(xi-xj).
d=Ei+&
E=E&’
(lhb)
i=-au ax’
(12)
where Q is the Cauchy stress, E,, and E are the intact and current elastic moduli, respectively, c is the logarithmic strain, and 4 is the integrity parameter which describes the stiffness degradation and lies in the range: 0 < 4 < 1; the undamaged state is represented by 4(O) = 1. In eqn (11) it is noted that due to the lack of spin in this one-dimensional problem, any objective stress rates reduce to the material time derivative of the Cauchy stress. (c) The damage evolution law 4 = -E,,q5&,
(13)
where 5 (> 0) is the damage coordinate and evolves according to
b k(x,, xi) = ko 3 rc exp - {(xi -xi)* J +(x2-x;)*+(xg--xi)*}.
where p is mass density; (b) The constitutive relation
(6)
where dt” is the increment of the damage coordinate in the ny direction. The evolution of the damage coordinate is described by the following nonlocal expression
(10)
;i;;=p",
(5a, b)
The damage evolves linearly with the damage force according to d4” = - Qt d<“,
VALUE PROBLEMS
Consider one-dimensional wave propagation in an elastic-fracturing solid in the axial direction denoted by coordinate x. The axial deformation is described by a map, x = x(X, t), between the material (Lagrangian) coordinate, X, and the spatial coordinate, x, at time I (see Fig. 1). Axial Cauchy stress, strain rate, and velocity are denoted by Q, i and v (= 2), respectively. In what follows, overdot implies material time differentation. The basic equations of the damage theory pertaining to the problem of wave propagation in one dimension follow. (a) The equation of motion
ci = nyn;cij, no sum on a.
INITIAL
(9)
Details of the above theory which is developed within the framework of irreversible thermodynamics by imposing the Clausius-Duhem dissipation in-
[(x, t) = max 0, (s
k(x - x’)i(x’, t) dx’ >
for .5(x, t) > 0, L(x, t) > 0, = 0
otherwise
Undeformed Configuration X
IX
(14) Deformed
-D
x,=x(X,.t)
Fig. 1. A deformation map.
Configuration
417
A nonlocal elastic damage theory in which the integration is performed over the onedimensional region, (x,,, x, ), occupied by the solid. The kernel function, k(x) is given by
k
z k(x) = k,
J$
exp(-bx*),
(15)
where k,, and b are material constants; in particular, the constant l/,/b represents a characteristic length of nonlocal interaction for the evolution of the damage coordinate t. The solution of the above equations necessitates the specification of (a) appropriate boundary data at x = x0 and x, and (b) initial conditions on x and v at time f = 0. It was shown in [4] that the equations of motion expressed by the displacement become unconditionaIly hyperbolic due to the nonlocal damage evolution, and the wave speed is given in terms of the secant modulus E, instead of the tangent modulus (see Appendix A). A corresponding local damage theory may be obtained from the above by setting the kernel function equal to the Dirac delta function k(x) = k,&x).
(16)
For both nonlocal and local damage theories, the energy dissipation W, of the solid, with an initial cross-sectional area A0 , becomes
Qd dX dt,
W,=A,
(17)
ss where the damage force Q is defined as
(18) The initial boundary value problem defined by eqns (lO)-(lS) for the nonlocal damage theory and eqns (10x14) with eqn (16) for the corresponding local theory is solved numerically by employing an explicit, updated Lagrangian finite element method. 4. UNIAXIAL
STRJBS-STRAIN
RELATION
In this section, the parameters of the above elastic damage model are related to quantities which can be measured by uniaxial tension and compression tests. It is noted, however, that the above ideal model was introduced as the simplest constitutive model which exhibits strain softening for the investigation of meshsensitivity in initial value problems. When the specimen length is large enough compared with the nonlocal reference length l/,/b, the local elastic damage theory becomes applicable, in which case eqn (15) can be employed in eqn (14). Therefore, parameters of the damage theory, except 6, can be determined by monotonic uniaxial tension and com-
b
“0
200
loo e
300
400
(PI
Fig. 2. Cyclic tensile stress-strain relation. pression tests. The initial Young’s modulus E,, can be determined by the slope of the stress-strain curve at t = 0. The parameter k0 in eqns (14) and (15) is related to the critical strain t, in tension [4], in terms of the following c = E,,c exp(-fk,E,,c’),
k0 = 1/(2E,t,‘).
(19)
According to eqn (14) damage evolution in compression is suppressed. Therefore, loading response in compression is perfectly elastic for the present model; if necessary for the modeling of structural elements, the damage evolution in compression can be introduced by modifying the damage evolution law, eqn (14). In the sequel, example calculations are presented for the material with E,, = 34.5 GPa, p = 2400 kg/my, and ce = 0.12 x 10-3. The tensile loading, unloading and reloading responses of the elastic damaging material are presented in Fig. 2. The unloading and reloading responses in tension are elastic, but with reduced elastic modulus. The determination of the nonlocal characteristic length l/Jb may require a suitable set of fracture tests. Since the objective of this paper is to investigate finite element mesh-sensitivity in initial value problems, it is assumed that the parameter b is given. 5.
FINITE ELEMENT
IMPLEMENTATION
The governing equations of motion, (IO), and stress-strain relations, eqns (11) and (12) are in local form. Therefore, updated Lagrangian formulation based upon the principle of virtual work holds. Consequently, the new damage model can be easily incorporated into existing finite element codes by simply adding a subroutine for the nonlocal damage evolution, eqns (13) and (14). An explicit, updated Lagrangian code was developed by using two-node elements. Since the implementation of updated Lagrangian codes is well documented [g-11] (see Appendix B), only the implementation of the damage evolution law (14) is briefly explained. For a typical integration point at x, where i and c are both positive, the increment of the damage coordinate A{ under time increment At is
H. MURAKAMI et al.
418
obtained from eqn (14) by using two-point Gauss quadrature over the elements which cover the support of k (x -x’)
-0.0 VI 2
M(x) = C i k(x - x’(tl,))A+‘(~J) c 6 -I
-0.1
_ _
O.Im
_
0.02s
-
0.0125m
0.05 m m
, (20)
$ 1
where the first summation is over the elements on the support, the second summation is over the Gauss integration points, L, is the element length, and 200 Time(p)
q*= -a,=‘.
(21) J5
This scheme, employed by the authors [12], is economical and can be easily generalized for two- and three-dimensional regions. However, the above numerical integration becomes inaccurate under extreme element stretch in the final stage of tensile failure. Therefore, an alternative numerical scheme with added computational time is presented in what follows. In the present numerical investigation, the nonlocal integration was evaluated analytically on the premise that Ae is known only at the integration points. For one point integration, A6 is constant over the element, while for two point integration, de is constant over the half region of the element. Denoting the range (xi, xi) at the pth integration point, eqn (14) yields
At(x)=C
c
i
9
Mx’(~J)ko
p=I
J XIt X 5
exp{ -b(x
-x’)*)
dx’
(22)
x:
which, upon evaluation of the integral, results in the expression
-erf(&(x
- xi)}].
(23)
The second summation in eqns (22) and (23) is over the integration points, and erf{*} denotes the error function whose value may be found numerically by a
c--x 0 Fig. 3(a). Problem 1: tensile loading of a semi-infinite bar.
+
-L
01
c
L
Fig. 3(b). Problem 2: symmetric tensile loading of a bar of
finite length.
Fig. 4. Velocity history at X = 0.4 m for c = 0.32m/sec computed by using different meshes.
rational approximation (see Appendix C) shown, for example, by Abramowitz and Stegun [13]. 6.
NUMERICAL RESULTS
In the numerical studies, attention was focused on the two problems, illustrated in Figs 3(a and b), involving a semi-infinite bar (Problem 1) and a finite bar (Problem 2). These problems were selected because they offer a means of studying the behavior of the damage models when strain softening occurs at a surface (Problem l), or in the interior of a material (Problem 2). In the subsequent finite element (FE) analyses uniform, initial meshes are employed to investigate the mesh sensitivity. In Problem 1, shown in Fig. 3(a), a particle velocity c is suddenly applied at time t = 0 to the end of a semi-infinite, initially quiescent bar, producing a tensile wave that travels to the right in the bar. In the FE approximation, the semi-infinite region was approximated by a finite region of length L = 1.6 m. The boundary conditions imposed in the FE analysis are u(x(0, t),
t) = -cH(t),
u(x(L, t), t) = 0
(24a, b)
where H(t) (= 1 for t > 0 and 0 for t < 0) is Heaviside’s step function. The nonlocal characteristic length l/Jb = 0.1 m was employed in the calculation. A constant time step defined by the Courant number computed by the intact elastic modulus is used. The calculation was terminated when the wave front hit the fictitious boundary X = L. This problem in compression, instead of in tension, was originally proposed by Sandler and Wright [l] as a means of illustrating the difficulties pertaining to loss of hyperbolicity when local strain-softening models are employed. The critical velocity to cause strain softening is approximately 0.326m/sec. First, a velocity of c = 0.32 m/set was chosen to ensure that the material does not enter the strain softening region. The velocity history at X = 0.4 m, computed by employing the local damage theory, is shown in Fig. 4 for different initial mesh size AX = 0.1 m, 0.05 m, 0.025 m, and 0.0125 m, respectively. The results show
419
A nonlocal elastic damage theory that as mesh size decreases the numerical solution converges to a unique ~lution; this was also true for the nonlocal damage theory. Next, the calcufation was repeated with a larger velocity, c = 0.6 m/set so that the material experiences strain softening. Figures S(a and b) show the calculated velocity histories at X = 0.4 m, by employing the local and nonlocal theories, respectively. Reference to Fig. S(a) reveals that the solution based upon the local theory shows no evidence of converging to a physically realistic solution as the mesh size AX is reduced; instead it appears to be approaching a vertical line. The result of the nonlocai theory, shown in Fig. S(b), on the other hand, clearly show convergence to a unique, physically realistic solution, Note that the solution obtained with the local model for AX = 0.1 m is quite cbse to the converged solution for the nonhz-cal modek The corresponding stress histories at the integration point of the right side element adjacent to X = 0.4 m are shown in Figs 6(a and b) for the local and the nonlocal theories, respectively. With respect to convergence, the results illustrate the same trend as the velocity histories. The dissipated energy is shown in Fig. 7 for A0 = I m2. As shown in this figure, the local solution converges to the situation in which there is zero energy dissipation
Fig. 6(b). Stress history at the right integration point adjacent to X = 0.4 m, predicted by the nonlocal theory for c = 0.6 m/w, and computed by using different meshes.
0.0 n
2 -
A
-0.1
.E
0
J
fig. 6(a)* Stress history at the right integration point adjacent to X ~5:0.4m, predicted by the local theory for c = 0.6 m/xc, and computed by using differeat meshes.
-0.2
5.3
0
loo
_
d.Im
_ _
0.0sIll 0.025m
-
0.012cm
200
300
400
TimeIp5)
Fig. 5(a). Velocity history at X = 0.4 m, predicted by the locak theory for c = Wm/sec, and computed by using different meshes.
as the bar edge softens and eventually breaks. This situation is not physically realistic for materials with heterogeneous ~c~s~cture. The nonlocal solution, shown in Figs 5(b) and 6(b), on the other hand, converges to a solution where softening occurs in a finite edge region characterized by the nonlocal characteristic length l/,/b = 0.1 m and finite dissipation energy is recorded as the mesh size tends to zero. Figure 7 atso illustrates, for comparison, the dissipation energy for the case, shown in Fig. 4, with no strain softening.
0.0 t.3
-c -a. .z
49.1
4
9
5.2
a.3
0
100
200
300
400
0 v.00
Fig. S(b). V&city history at X = OS4m, predicted by the
nonlocal theory for c = 0.6 m/set, and computed by using different meshes.
0.05
5.10
0.15
0.20
AX Cm) Fig. 7. Dissimtion energy veraus mesh size (Problem 1).
420
H. MUMKAMer al.
Problem 2, illustrated in Fig. 3(b), consists of the symmetrical tensile loading of a bar of finite length 2L. The initially quiescent bar is suddenly loaded by applying end velocities, - c and c, at time f = 0. The boundary conditions are u(x(-
2-
1-
L,t),t)=-cH(t), v(x(L,t),t)= d(t). (25a,b)
The velocity c is selected so that (a) the maximum stress in the initial waves that propagate from the ends toward the center at x = 0 is less than that necessary to cause strain softening and (b) the maximum stress that results when the two incoming waves collide at x = 0 is large enough to cause strain softening. In the present study, we set c = 0.3 m/set and took the length of the bar 2L to be 1.6 m. This problem was originally studied by Baiant and Belytschko [5], where an exact solution was obtained for the case in which the material behavior was described by a local strain softening model. The calculated velocity histories at a quarter point, X = 0.4 m, by the local theory and the nonlocal theory are shown in Figs 8(a and b), respectively. The corresponding stress histories at the closest integration point to the right of X = 0.4 m are illustrated
-1
0
100
200
300
Time ( p
400
)
Fig. 9(a). Stress history at the right integration point adjacent to A’ = 0.4m, predicted by the local theory for c = 0.3 m/set, and computed by using different meshes.
2
_
x
cn
-I
0
loo
O.lm
_
0.05
_
O.OZSm
-
0.012Sm
200
111
300
4co
Time ( ps )
0.2
Fig. 9(b). Stress history at the right integration point adjacent to X = 0.4 m, predicted by the nonlocal theory for c = 0.3 m/xc, and computed by using different meshes.
-0.6 0
IO0
200
300
400
Time(p)
Fig. 8(a). Velocity history at X = 0.4 m, predicted by the local theory for c = 0.3 m/set, and computed by using different meshes. 0.2
0.0
_ _
oosm
-
00125 111
in Fig. 9(a) for the local theory and in Fig. 9(b) for the nonlocal theory. As the two tensile waves meet the specimen breaks into two pieces, as shown in the stress histories in Fig. 9. Both the local and the nonlocal theories predict that splitting. The local theory predicts vertical stress drop reflecting localized failure in an infinitely thin layer. On the contrary, the nonlocal theory predicts the splitting failure in a finite volume as shown by the slanted stress drop in Fig. 9(b). This scenario is better explained by showing the dissipated energy in Fig. 10. The local theory converges to the finite dissipation energy (reference data)
Ol125m
-0.2
-0.6 0
loo
200 Time (pi
3ocl
400
)
Fig. 8(b). Velocity history at X = 0.4 m, predicted by the nonlocal theory for c = 0.3 mjsec, and computed by using different meshes.
Fig. 10. Dissipation energy versus mesh size (Problem 2).
A nonlocal elastic damage theory
which is predicted when the bar does not experience any strain softening (c = 0.16 m/set). Finally, it is informative to mention the additional computational time required to perform the calculation using the nonlocal damage evolution, eqns (14) and (15), relative to the local damage evolution, eqns (14) and (16). For the fine mesh calculations, AX = 0.025 and 0.0125 m, the nonlocal analyses typically take 30-50 times more CPU time than the local analyses. This order of magnitude difference is due to the error function calculation (C2). Therefore, it is expected that by creating a look-up table for the error function and eliminating the repeated evaluation of the error functions at adjacent elements, the additional computational effort required for the nonlocal damage evolution will reduce to less than a factor of ten. 7. CONCLUSIONS
This study consisted of applying the nonlocal elastic damage theory to uniaxial stress problems and examining the performance of the constitutive model in wave propagation studies, especially under strain softening. From the results of this study, the following conclusions are drawn: 1. The model is computationally feasible for use in conjunction with finite element wave codes. 2. In wave propagation problems involving strain softening, the model leads to (physically plausible) unique solutions. This was dramatically illustrated in the numerical analysis of the problem involving the tensile loading of a semi-infinite bar (Figs 5 and 6); as the mesh size was reduced, the nonlocal solution converged to a physically realistic solution, while the local solution did not. 3. It is recommended that this damage theory be expanded to incorporate plastic damaging process to simulate the mechanical response of rocks, concrete, and ceramics.
421
5. Z. P. Baiant and T. B. Belytschko, Wave propagation in a strain-softening bar: exact solution. ASCE I. Engng
111, 381-398 (1984).
Mech.
6. Z. P. Baiant, T. B. Belytschko and T:-P. Chang, Continuum theory for strain-softening. ASCE I. Engng Mech. 110, 1666-1692 (1984). 7. K. C. Valanis, A theory of damage in brittle materials. Engng Fract. Mech. 36, 403-416 (1990). 8. K. J. Bathe, Finite Element Procedures in Engineering Analysis. Prentice-Hall, Englewood Cliffs, NJ (1982). 9. T. J. R. Hughes, The Finite Element Method, Linear Static and Dynamic Finite Element Analysis. PrenticeHall, Englewood Cliffs, NJ (1987). 10. J. 0. Hallquist, Theoretical manual for DYNA3D. Report UCID-19401, University of California, Lawrence Livermore National Laboratory, Livermore, CA (1982). 11. G. L. Goodreau and J. 0. Hallquist, Recent development in large-scale finite element Lagrangian hydrocode technology. Comput. Meth. Appi. Mech. Engng 33, 725-757
(1982).
12. H. Murakami, D. M. Kendall and K. C. Valanis, Numerical validation of a nonlocal damage theory. In Constitutive Laws for Engineering Materials, Recent Advances and Industrial and Infrastructure Applications
(Edited by C. S. Desai, E. Krempl, G. Frantziskonis and H. Saadatmanesh), pp. 743-746. ASME Press, New York (1991). 13. M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions, pp. 297-329. Dover, New York (1964).
APPENDIX
A: UNCONDITIONAL
HYPERBOLICITY
The type of initial value problem defined by eqns (lO)-(lS) can be examined by investigating that of the corresponding rate problem. Upon monotonic loading, eqns (llb), (13) and (14) yield R = -2EE,c*
k(x - x’)i(x’,
r)dx’,
s
(Al)
where the integration is performed over (x,,, x,). By substituting the above equation into eqn (1 la), a rate constitutive relation is obtained in the form ti = E(i -2E&j
k(x - x’)C(x’, t) dx’ )
(A2)
From eqn (10) the equation of motion becomes, in rate form Acknowledgement-The
authors express their gratitude
ad
to Dr H. E. Read of S-CUBED, a division of Maxwell
ax=pa,
Laboratories, La Jolla, California, for his helpful advice and comments. By substituting obtains
a=d.
eqn (A2) into the above equation,
one
REFERENCES
1. I. S. Sandler and J. P. Wright, Strain softening. In Theoretical Foundation for Large-scale Compuraiions for Nonlinear Materials (Edited by S. Nemat-Nasser,
R. J. Asaro and G. A. Hegemier), pp. 285-296. Martinus Nijhoff, Dordrecht (1984). 2. F. H. Wu and L. B. Freund, Deformation trapping due to thermoplastic instability in one-dimensional wave propagation. J. Mech. Phys. Solirls 32, 119-132 (1984). 3. H. E. Read and G. A. Hegemier, Strain softening of rock, soil and concrete-a review article. Mech. Materials 3, 271-294 (1984). 4. K. C. Valanis, A global damage theory and the hyperbolicity of the wave problem. ASME I. Appl. Mech. 58, 311-316 (1991).
E{;-ZE,,r’Jk*(x-x’)t(x’,i)dx’
-6E,c’;
k(x - x’)((x’,
t) dx’
s +~{G2E0r3~k(x-x’)c(x’,~)dx’}=pd,
, 1 (A4)
where k*(x) = dk(x)/dx. In terms of the velocity field v(x, r), the wave equation (A4) is of the form E&p$+E;+F
(A5)
H. MURAKAMIet cl.
422
where F is a spatial function for au/ax over the material domain D and depends, parametrically, on the strain field 6, where
t sa0
c=
dt.
F<=j;[-&Nt$N2]6dx,
,ax
Given the functional F{ .} in D at I, the character of eqn (A5) in the infinitesimal nonvanishina time interval At, is determined by the higher order derivative terms’ on the left-hand side of eqn (A5). Therefore, eqn (A5) is hyperbolic if E, the secunt modulus, and p, the density, are positive. This is the case even under conditions of material softening. Furthermore the wave speed C,, is given in terms of E and p by eqn (A7)
c,=
The evaluation of the integrals of eqn (Bl) on a typical element (e) yields the internal force vector F, and the diagonalized mass matrix M,
M,=-
P-G
1
2 [0
(B5) 0
1
1’
E
where the superscript T denotes transpose and L=(=x2 - x’) is the length of the element (e). By assembling the nodal coordinate vector X, the nodal velocity vector V, the nodal acceleration vector A ( = V), and the internal force vector F with respect to the global degrees of freedom, one obtains the following equations of motion atatimestepf=t”
P
MA” = I”’ - F”,
J -
and is always real.
APPENDIX B: EXPLICIT, UPDATED LAGRANGIAN FINITE ELEMENT SCHEME
The principle of virtual velocity applied to the initial value problem defined on the current configuration (x0, x,) yields
+ Wx,, t)T(x,, 1)-WA,,
r)W,, lb
@I)
where 6u denotes the virtual velocity, x0 = x (X0, I), xi = x (X,, I), and T is the x component of the traction vector acting on x = x,, or x, . The finite element discretixation of the above weak equation has been explained by many authors, for example, Bathe [8] and Hughes [9]. In what follows, the outline of the procedure is briefly summarized. The one-dimensional domain (x0,x,) is discretized into two-node elements
(W
Vn+1/2=Vn-1/2+f(Atn+Atn+I)An
(B9)
whereAt,(=t,-t,_,)andAr,+,(=r,+,-r,)aretimeincrements: APPEND=
c: ERROR FUNCTION
The error function is defined in eqn (Cl)
0 (&-a9
exp( - t2) dt.
cc11
r=l
where N, is the total number of the elements, and x$, k = 1.2. is the k th nodal coordinates denoted in vector form as X,1 In a typical isoparametric element (e), the coordinate and the velocity are defined on the standard region - 1 C r Q I
as
It can be evaluated numerically by employing the following rational approximation [13]
xexp(-n2)+t(x) r = &,
where af (k = 1,2) are the elemental nodal-velocities, denoted in vector form as V,, and NL(r) (k = 1,2) are shape functions defined on the standard region as follows: N’({) = f(1 - 5)
(B7)
where M is the diagonal, global mass matrix, and P” is the global, external force vector due to the traction T in eqn (Al). For the class of problems treated in this paper, velocity boundary conditions are prescribed. Therefore, the last two terms on the right-hand side of eqn (Al) vanish, since the virtual, boundary velocities are zero, &(x0, t) = du(x,,f)=O; this yields P”=O. The explicit, central difference time integration scheme, implemented, for example, in DYNA2D and DYNA3D [lo, 111, is summarized as follows: A”=M-I@‘“-F”)
(x,7x,)=
(W
N?(T) = f(l + 0.
Wa, b)
k(x)] c 1.5 x 10-7,
where the constants are defined as follows: p = 0.3275911, a, = 0.254829592, n2= -0.284496736,
oS = 1.421413741,
a. = - 1.453152027. or = 1.061405429.
(C3) \ I