- Email: [email protected]
Tensile cracks in polycrystalline ice under transient creep: Part I — Numerical modeling
MECHANICS Of MATERIALS ELSEVIER
Mechanics of Materials 21 (1995) 265-279
Tensile cracks in polycrystalline ice under transient creep: Part I - Numerical modeling S. Nanthikesan a,1, S. Shyam Sunder h,2 a Prudent Engineering Group, Liverpool, IVY 13088, USA b Department of Civil and Environmental Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139, USA
Received 30 August 1993; revised version received 1 June 1995
Abstract The interaction between creep deformations and a stationary or growing crack is a fundamental problem in ice mechanics. Knowledge concerning the physical mechanisms governing this interaction is necessary: (1) to establish the conditions under which linear elastic fracture mechanics can be applied in problems ranging from ice-structure interaction to fracture toughness testing; and (2) to predict the ductile-to-brittle transition in the mechanical behavior of ice and, especially, the stability and growth of cracks subjected to crack-tip blunting by creep deformations. This requires a quantitative estimate of the creep zone surrounding a crack-tip, i.e., the zone within which creep strains are greater than the elastic strains. The prediction of the creep zone in previous ice mechanics studies is based on the theory developed by Riedel and Rice (1980) for tensile cracks in creeping solids. This theory is valid for a stationary crack embedded in an isotropic material obeying an elastic, power-law creep model of deformation and for a suddenly applied uniform far-field tension road that is held constant with time. The deformation of ice at strain-rates ahead of a crack (i.e., 10 -6 to 10 -2 s -1) is dominated, however, by transient (not steady power-law) creep and the loading, in general, is not instantaneous and constant. A numerical model is developed in this paper to investigate the role of transient creep and related physical mechanisms in predicting the size, shape and time evolution of the creep zone surrounding the tip of a static crack in polycrystalline ice. The model is based on the fully consistent tangent formulation derived in closed form (Shyam Sunder et al., 1993) and used in the solution of the physically-based constitutive theory developed by Shyam Sunder and Wu (1989a, b) for the multiaxial behavior of ice undergoing transient creep. The boundary value problem involving incompressible deformations ahead of a stationary, traction-free mode I crack in a semi-infinite medium is modeled and solved by a finite element analysis using the boundary layer approach of Rice (1968). This model is verified by comparing its predictions with (i) the known theoretical solutions for the elastic and HRR asymptotic stress and strain fields in an elastic-plastic material of the Ramberg-Osgood type, and (ii) the creep zone size for an isotropic material obeying the elastic power-law creep model of deformation. Keywords: Fracture mechanics; Ice mechanics; Creep fracture; Fracture process zone; Crack-tip inelasticity;
Polycrystalline ice; Fracture mechanics of ice; Prediction of crack-tip creep zone
i Senior Engineer, Prudent Engineering Group, Liverpool, NY 13088, USA. Formerly Postdoctoral Associate, Department of Civil and Environmental Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139, USA.
2 Senior Research Scientist, Department of Civil and Environmental Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139, USA. Now at the Building and Fire Research Laboratory, National Institute of Standards and Technology, Gaithersburg, MD 20899, USA.
0167-6636/95/$09.50 © 1995 Elsevier Science B.V. All rights reserved SSD1 0167-6636(95)00014-3
266
S. Nanthikesan, S. Shyam Sunder~Mechanics of Materials 21 (1995) 265-279
1. Introduction
The interaction between creep deformations and a stationary or growing crack is a fundamental problem in ice mechanics. Knowledge concerning the physical mechanisms governing this interaction is necessary: (1) to establish the conditions under which linear elastic fracture mechanics (LEFM) can be applied in problems ranging from ice-structure interaction to fracture toughness testing; and (2) to predict the ductile-to-brittle transition in the mechanical behavior of ice and, especially, the stability and growth of cracks subjected to crack-tip blunting by creep deformations. A qualitative discussion regarding the applicability of L E F M to predict ice forces on offshore structures is provided by Palmer et al. (1983). They characterize the modes of deformation and, in particular, the relative role of creep and fracture processes as a function of loading rate and the problem geometry. For impact times of less than about 10-100 seconds, this work suggests that creep deformations are often negligible and only (linear elastic) fracture processes dominate. The need to analyze the creep-crack interaction problem more rigorously has been recognized by experimentalists interested in characterizing the fracture resistance (i.e., "toughness") of polycrystalline ice. The inelastic (creep-dominated) fracture process zone ahead of a crack must be small relative to the key geometric dimensions of a test specimen to ensure small-scale yielding (SSY) conditions and the validity of LEFM. This requirement is explicitly considered by Goodman (1980), Urabe and Yoshitake (1981), Nixon and Schulson (1987) and Timco and Frederking (1986). More recently, Schulson (1990) predicts the strain-rate at which the ductile-to-brittle transition occurs in polycrystalline ice under uniaxial compression. His wing crack model of brittle compressive failure assumes that a crack can propagate only if the crack-tip creep zone is smaller than a critical size that is proportional to the microstructural scale of the aggregate. Otherwise, there is no growth and the crack remains
stable. A similar argument is used by Schulson et al. (1989) to explain the stability of long nonpropagated cracks under uniaxial tension. The model allows for an increase in fracture toughness due to crack-tip blunting; more specifically, it assumes that, as in steel, the toughness is proportional to the square-root of the crack-tip radius and which, in turn, is proportional to the "radius" of the creep zone. The prediction of the creep zone, defined as the zone surrounding the crack-tip within which creep strains are greater than elastic strains, in the above studies is based on the theory developed by Riedel and Rice (1980) for tensile cracks in creeping solids. Riedel and Rice present approximate analytical solutions for the creepdominated near-tip asymptotic stress and strain fields (the RR fields) under small-scale yielding conditions as well as for the boundary of the creep zone. The solutions are valid for a stationary crack embedded in an isotropic material obeying an elastic, power-law creep model of deformation and one that is subjected to a suddenly applied uniform far-field tension load which is held constant with time. These theoretical predictions were verified by Bassani and McClintock (1981) using finite element analysis. Riedel (1981) and Smith (1990) have extended the theory to include the effects of isotropic strain hardening during creep. The deformation of ice at strain-rates of engineering interest (i.e., 10 -6 to 10 -2 s -a) is dominated, however, by transient (not steady powerlaw) creep and the loading, in general, is not instantaneous and constant. The strain-rates in the vicinity of a crack fall within this range (as will be shown subsequently). Consequently, any analysis of the inelastic processes surrounding the crack-tip must include the effects of transient creep which involves isotropic and kinematic hardening. The analysis must also consider the effects of loading history, temperature sensitivity, variation in material strain hardening among different types of ice, and fabric anisotropy. This paper presents a numerical model to study the role of transient creep and related physical mechanisms in predicting the size, shape and
S. Nanthikesan, S. Shyam Sunder~Mechanics of Materials 21 (1995) 265-279
time evolution of the creep zone surrounding a crack-tip in polycrystalline ice. The physicallybased constitutive theory of Shyam Sunder and Wu (1989a,b) for the multiaxial creep of ice undergoing transient and viscous creep is used in conjunction with a finite element analysis to model the creep-crack interaction problem. The solution to the c r e e p - c r a c k interaction problem using this constitutive theory involves the numerical time integration of a coupled set of very stiff and highly nonlinear first-order differential equations. These equations, governing the evolution of the observable and internal state variables, characterize local material behavior at each integration point in a finite element mesh. A closed-form N e w t o n - R a p h s o n (tangent) formulation in conjunction with an implicit form of the o~-method of time integration is developed to solve the local constitutive equations. A consistent constitutive Jacobian matrix, which is used to assemble the global tangent stiffness matrix in a finite element analysis, is also established in closed form. The algorithm is implemented as a user-specified subroutine in the finite element program A B A Q U S and its numerical stability, accuracy and convergence properties have been verified by Shyam Sunder et al. (1993). The problem of a stationary, traction-free, tensile (mode I) crack is modeled using the equivalent boundary layer idealization proposed by Rice (1968). The crack-tip elements and the mesh configuration are chosen so as to model the singular stress and strain fields within the boundary layer accurately. Finally, the ability of the numerical model to predict the elastic and H R R asymptotic stress and strain fields and the creep zone for an isotropic material obeying a steady power-law creep model of deformation is validated against known analytical solutions. The companion paper (Nanthikesan and Shyam Sunder, 1995) presents the results from a comprehensive set of numerical simulations that analyze the role of transient creep, including loading history, temperature sensitivity, material strain hardening and fabric anisotropy, in predicting the size, shape and the time evolution of the creep zone in polycrystalline ice.
267
2. Summary of the constitutive model
A summary of the physically-based constitutive model of creep in polycrystalline ice developed by Shyam Sunder and Wu (1989a, b) is presented in this section. The mathematical formulation of the transient creep deformation mechanism, founded on the theory of internal state variables, is characterized by: (a) a scalar internal state variable B, the drag stress, to describe isotropic hardening; (b) a second-rank tensorial internal state variable R, the elastic back stress tensor, to describe kinematic hardening or deformation-induced anisotropy; and (c) a generalized equivalent stress measure to describe fabric anisotropy of the Hill-type in the solid. In addition to transient creep, the model takes into account the elastic and steady viscous creep deformations through a summation of the strain-rates for each of the deformation mechanisms. Transient creep in polycrystalline ice is anelastic and fully recoverable upon unloading while viscous creep is irrecoverable. Only small deformations are considered since the magnitude of the transient creep strain does not exceed 1%. At larger strains, the material exhibits strain-softening by progressive distributed microcracking or recrystalization processes if premature failure by brittle fracture does not occur. Readers are referred to the previously referenced papers, as well as Shyam Sunder and Wu (1990) and Wu and Shyam Sunder (1992), for the physical basis of the constitutive model and its material parameters. Matrix notation is used in this paper. A superposed dot denotes a time derivative and the subscript " e q " refers to a scalar equivalent value. For small deformation theory, the total strainrate, i, is additively decomposed as follows: = ~e + iv + ~t, (1) where the subscripts "e", " v " and " t " denote the elastic, viscous and transient components of the strain-rate, respectively. The elastic strain-rate is related to the stress-rate, d-, by:
b=D~e=O(~-iv-~,) where D is the elastic stiffness matrix.
(2)
268
S. Nanthikesan, S. Shyam Sunder / Mechanics of Materials 21 (1995) 265-279
The viscous strain-rate vector, iv, follows the classic power law (Glen, 1955). For an incompressible orthotropic material it is given by: • -iv = ~'0/3
N-I 0"eq S,
(3)
S = Go', (4) where ~0 is a scalar reference strain-rate and can be taken as unity with no loss of generality; /3 = a 1 + a2, where a I and a 2 (as well as a3, a 4, a 5 and a 6) are constants which describe the Hilltype fabric anisotropy of the material (these orthotropic parameters are calculated from experimental data, see Shyam Sunder et a1.(1993)), N is the power law exponent; cr is the stress vector {0-xx, 0-yy, o'~z, cr~y, 0-.~, 0-~y}T; S is the deviatoric stress vector; G is the matrix which transforms stresses into their "deviatoric" components and is given in the Appendix; V is the temperature-dependent flow stress in the reference direction (assumed to be the Y-direction in this paper) given by the Arrhenius relationship:
where V0 is a temperature-independent constant, Q is the activation energy, T is the temperature in Kelvin, and R is the universal gas constant; 0-eq is the equivalent stress measure given by: 3 T o'e2q= ~0- GO-,
(6)
where the superscript " T " denotes the transpose operator. The transient strain-rate vector, it, is given by: ,,
(7)
with Sa = G0-~ = G(0- - R). (8) The evolution equations for the back stress vector, R, and the scalar drag stress B are given by: 1 = ~ [ 3 A E H : i t,
(9)
Be /~ = O [ - - ] ~ t
I_O'o,o ]
eq,
'
(tO)
where 0-d is the reduced stress vector, A E is the anelastic modulus of kinematic hardening; E is the Young's modulus in the reference (Y) direction, and A is a non-dimensional material constant (in a uniaxial creep test, the ratio of the maximum transient creep strain to elastic strain is equal to 1//1); H is the strain (and strain-rate) transformation matrix which transforms the strains into their "deviatoric" components and is given in the appendix; B 0 is a non-dimensional material constant which represents the threshold value of the drag stress; O is a function of the drag stress which ensures that upon unloading (when magnitude of transient strain decreases), B tends asymptotically to the threshold drag stress, B~), and is given in Shyam Sunder and Wu (1990); H is a non-dimensional material constant which describes the rate of isotropic hardening; 0-d,eq is the equivalent measure of reduced stress obtained by: 3 O'2d,eq = ~-0"J60"d"
(11)
gt,eq is the equivalent transient strain-rate and is given by: ~2t,eq = ~ i T H e t ,
(12)
here, the strain-vector is expressed as {exx, e-yy, ezz, 7xy, Yxz, y~y}v. Note that the shear components, 7o, are engineering strains. 2.1. Material parameters
The constitutive model consists of the standard elastic moduli, the universal gas constant R, the flow stress V0, the power-law exponent N, and the Arrhenius activation energy Q. The parameters R, Q, V0 and N are associated with viscous creep behavior• The transient creep is characterized by three additional constants which model strain-hardening: (1) the anelastic constant for kinematic hardening, A, (2) the initial or threshold drag stress, B 0, and (3) the rate of isotropic hardening, H. In addition, there are the five independent Hill constants, a i (i = 1, 5), which characterize the fabric anisotropy of the material•
S. Nanthikesan, S. Shyarn Sunder~Mechanics of Materials 21 (1995) 265-279
The formulation of transient strain-rate given by Eqs. (8) through (11) considers the flow law to be the same for both deformation-rate mechanisms, the only difference being the existence of drag and back stresses in the case of transient flow. It also implies the following: (a) the exponent N is the same as that for viscous flow in Eqs. (3) through (7), and (b) the temperature dependence of the transient creep-rate, represented by the parameter V, is given by an Arrhenius law with an activation energy equal to that of viscous creep. These assumptions are justified in Shyam Sunder and Wu (1989).
3. Integration of constitutive equations To solve creep-crack interaction problems using the finite element method requires a numerical integration of the governing nonlinear, ratedependent constitutive equations summarized above. When the deformation is dominated by creep, as is the case for the problem under consideration, the material behavior is incompressible and displays significant strain hardening. The coupled set of first-order differential equations, which models the local constitutive response at each integration point in a finite element mesh, is very stiff and highly nonlinear. This complexity, in turn, poses severe demands on the global finite element formulation with respect to numerical stability, accuracy and rate of convergence. A comprehensive numerical model to solve ice mechanics problems dominated by transient creep deformations using the finite element method has been developed by Shyam Sunder et al. (1993) and is summarized in this section.
3.1. 77me integration scheme In a displacement-based finite element formulation, estimates of the incremental displacement fields (and thus, the incremental strain fields) are known at the beginning of each constitutive iteration within a time increment. The problem of marching the solution forward in time then reduces to computing the increment in stress (and in the internal state variables) corresponding to
269
the given increment in total strain. The constitutive equations are solved using the one-parameter class of integration algorithms known as the amethod of integration, with:
F~i~ = (1 - a ) t F + a'+~'F
(13)
where the weighting parameter c~ may vary between 0 and 1; in general, the symbol F represents the (Cauchy) stress vector, the backstress vector, and the non-dimensional drag stress; F~(i) denotes the weighted value of F between time t and t + A t after the i th iteration. The two superscripts on the left side of symbol refers to time; and the superscript in parentheses on the right side refers to the iteration number within the time step. If the latter superscript is missing, it refers to the converged value. This nomenclature will be used for the rest of this paper.
3.2. Newton-Raphson iteration for the stress increment The constitutive equations in differential form are expressed in terms of finite differences to solve for the stress increments. Eqs. (1) through (12) in conjunction with Eq. (13) will yield the solution for the increment in stress t+~t~r~i+1). Such a scheme has only linear convergence. In this paper, a summary of the algorithm based on the N e w t o n - R a p h s o n iteration scheme is presented. The rate of convergence is quadratic in the neighborhood of the actual stress t+a/~r. Rewriting Eq. (2) as a residual, and noting A denotes the increment from time t to time t + At: f(A,~)
= zXo, - D z X ~ e
= Ao--D(A~-
A e v - A~t).
(14)
The total strain increment AE is prescribed, and Aev, Aet are determined assuming a Ao~. The objective is to find the root of f ( A c r ) = 0 which provides an improved estimate of the stress increment. Expanding f(Ao-) in a Taylor series, considering only the first two (i.e., linear) terms, and rearranging yields the following equation: ,+~,~,i+,, =,+~,~,)_~,i,[CA~,,)_ A4i)],
(15)
270
S. Nanthikesan, S. Shyam Sunder~Mechanics of Materials 21 (1995)265-279
where C is the compliance matrix and /)(0 is the elastic-creep stiffness matrix: /~(i)=
C+a
~aAo_~i + ~ } j
(16)
The derivation of the expressions for the partial derivatives in Eq. (16) is given subsequently. 3.3. Convergence criterion
Convergence is considered to be achieved when the magnitude of change in stress increments between two successive iterations at time t + At, is smaller than a prescribed fraction (e.g. 1/1000) of the magnitude of stress at time t + At, i.e., Ao-(i +1) _ Ao-(i)
.< 0.001.
ture have established the accuracy of the numerical model in simulating the constitutive behavior of polycrystalline ice (Shyam Sunder et al., 1993).
4. Derivation of tangent stiffness matrices
This section summarizes the derivation of the fully consistent closed-form solution for the following: (a) the elastic-creep stiffness matrix used in the local N e w t o n - R a p h s o n iteration for solving the material constitutive equations, and (b) the constitutive Jacobian matrix used in the global iterations. The latter is required in any "implicit" finite element analysis using a Newton-Raphson iteration scheme to solve the weak form of the linear momentum balance equations.
(17) 4.1. Elastic-creep stiffness matrix
3.4. Numerical implementation and evaluation
The above constitutive model has been implemented as a user-specified subroutine, UMAT, in the finite element program ABAQUS (1988). This subroutine is called at each integration point in the finite element model for every global iteration. At the beginning of each global iteration, U M A T requires the input of the total strain increment and the state variables (stress vector, scalar drag stress and vector back stress). At the end of the iteration, U M A T updates and outputs these state variables, the transient strains and the transient strain-rates, as well as the constitutive Jacobian matrix. Only a fixed (i.e., constant) time increment can be specified by the user. Therefore, the automatic time incrementation scheme of A B A Q U S is used in this analysis. In this scheme, the time step is increased by a factor of 1.5 if global convergence is attained in less than half the user-specified number of iterations in two consecutive time increments. On the other hand, if convergence is not achieved within this prescribed number of iterations, the time increment is reduced by a factor of four. Extensive comparisons with experimental data and theoretical solutions published in the litera-
Rewriting Eq. (16) as: I) (i, = [ C q- F(vi) + F, (') ] -1
18)
where 19)
F (i) ~ a aAo.~i) ,
and
aAe{ i) (20)
F[ i) ~ a aAo.~i ) .
The terms F~(i) and Ft(i) contain the contributions of viscous and transient creep, respectively, to the elastic-creep stiffness matrix. The expressions for these terms are obtained using the product rule of differentiation and then the chain rule, and are summarized below: Fv{/) = g0At [A(/)G 9 (NJr- /~2 ~ V
1) [_(0 "~N-3c(i)[c(i)'~T] \~'a,eq} °or I.°a ) ] (21)
and F[ i , = [ M ~ )1 -'F~O,
(22)
S. Nanthikesan, S. Shyarn Sunder/Mechanics of Materials 21 (1995) 265-279
where
M~' = [I + ½fiAEF~i)H].
(23)
The discretized weak form of the equilibrium equations at time t + At can be written as: fu6ffT t+Ato" dT/=
I is the identity matrix and
F~0 = a[ M~'a)] -it ~,,i, + M(j)),
271
£6uT t+atf B d ~
(24)
+ fg&,T t+zXtfS d ~ ,
(28)
in which
M(j) = goAt I A(~IG k
9
(N-
1) t - t 0
+ 132 (B2i v)
~N-3S~O1' S(e.,,J 0 ~T]]
Wa .... (25)
M~,;)=i+(goAt)2I
3
[
N
(mB(i) 1
\---t,eq]
[X(i)'~2q(i) (~,(i) T H ] ' X ~,'~d ] ~d,oe\~d,a)
lift(i) = *"B2
3 coAt o~~
N
(26)
f//v~TJAE d~'F+ fv~ET to" d ~ ,
(29)
where 0Ao"
( AB (o
[,~_(i) "12 7 ~ d ,c~,eq ] Bct
X ),(i)K,(i {~(i) ,~T] .,d ~.d,.t ~.a,.j ].
where 66 is the variation in the strain field corresponding to the variation in the displacements 6u; the superscript "T", denotes the transpose operator, ,+~fB and ,+~tfs are the body forces and surface tractions, respectively; ~'is the volume and Y' is the surface area. Applied concentrated forces are included as a part of surface tractions in Eq. (28). The constitutive Jacobian matrix, J, is required when the first term in Eq. (28) is linearized around 'u. The linearized term is:
(27)
4.2. Constitutive Jacobian matrix The rate of convergence of "implicit" finite element procedures for nonlinear materials is quadratic when a N e w t o n - R a p h s o n type iterative method is used to solve the weak form of the equilibrium equations. This method exhibits a stronger rate of convergence than the initial stress method, the modified Newton iteration and the quasi-Newton scheme (William, 1978; Argyris et al., 1978; Bathe, 1982). However, the N e w t o n Raphson scheme requires the constitutive Jacobian matrix to be computed at every integration point at each global iteration. The computational effort is considerably reduced by using closedform expressions for the constitutive Jacobian matrix (see, e.g., Lush et al., 1989). This section outlines the key steps in deriving such an expression.
j =
(30)
tK = ~._~,Tj,.~ d ~ ' ,
(31)
0Ae The global Newton iteration scheme is obtained from Eqs. (28),(29), as described in detail in Shyam Sunder et al. (1993) and the stiffness matrix, tK, is written as:
where the B matrix relates the strain fields to the nodal strains fi, as given by • =.~'t~. The stiffness matrix for the system is assembled from the contributions of each finite element. The use of Gaussian quadrature to perform this integration necessitates the computation of J at each integration point within the element. The J matrix for the material model used here is derived by differentiating the difference form of Eq. (2) and Eq. (30). This yields: [ J =D
I
0A'v 0Ae
0A't ] 0Ae "
(32)
Using the chain rule of differentiation as in the previous section, the expression for J is obtained as; J = [C + F v + F t ] - '
(33)
S. Nanthikesan, S. Shyam Sunder/Mechanics of Materials 21 (1995) 265-279
272
It can be seen that this is identical elastic-creep stiffness matrix in matrices bv and F t are given by (20), respectively, when the lower has converged.
in form to the Eq. (18). The Eqs. (19) and level iteration
(~ :
UnfformF~r-FieidTensLieS~'~ss
5. Crack-tip model under small-scale creep
5.1. Formulation of the boundary value problem In this paper, the boundary value problem of a static crack under mode I loading is considered. To obtain the response in the vicinity of the crack-tip, an equivalent but simpler problem of a boundary layer (Rice, 1968) is postulated and solved. U n d e r small scale creep (SSC) conditions, the L E F M stresses dominate at distances from the crack-tip that are large compared to the creep zone size but still small compared to the dimensions of the crack, a. In this simplified model, the actual boundary is replaced by one that is circular and centered around the crack-tip with the actual boundary conditions replaced by the L E F M asymptotic field along this circular boundary. The circular boundary is chosen sufficiently far from the crack-tip where L E F M conditions prevail. This model focuses on the immediate vicinity of the crack-tip and hence is an efficient and convenient idealization to model the asymptotic stress and strain fields in the creep-crack interaction problems of interest. For the purposes of the finite element analysis, the simplified configuration is represented by the crack -R<_xl <0, and the circular boundary centered at x 1 = x 2 = 0 with radius R (see Fig. 1). Here, R is chosen large enough (a >> R >> Rcr) so that the deformation at the boundary is unperturbed by the creep zone. The crack is traction free and the L E F M displacement field, corresponding to the applied load on the actual configuration, is applied at the circular boundary. This displacement field is given by: ui = ~
fi(O)
with i=
1, 2,
fl(O) = (3 - 4u - cos O) cos(O/2), f2(O) = (3 - 4 u - cos O) s i n ( O / 2 ) ,
0" Fig. 1. Schematic of creep-crack interaction problem for a tensile mode I crack in a semi-infinite medium modeled using the bounday layer approach of Rice (1968).
where Kj is the stress intensity factor determined from the actual configuration; G is the shear modulus; u is the Poisson's ratio; and 0 is the angle sustained with the x 1 axis in the counter clockwise direction. Symmetry of this mode I configuration allows the analysis of just a halfplane lying either above or below the plane of the crack, with the boundary conditions on x 2 = 0,x 1 > 0 given by u 2 = o-12 0. U p o n application of the load, the instantaneous response of the material is purely elastic and the crack-tip field is given by LEFM. Sustained application of the load, smaller than that required for crack propagation, gives rise to creep deformations that relax the stresses. The relaxation is most significant in the immediate vicinity of the crack-tip where the strain-rate (and stress) field is singular. This creep-dominated zone grows with time and may eventually occupy the entire solid. During growth, the creep zone is bounded by the L E F M field. =
5.2. Finite element mesh design (34)
Singular stress and strain fields exist within the boundary layer. The singularities associated with
S. Nanthikesan, S. Shyam Sunder~Mechanics of Materials 21 (1995) 265-279
/ 1 0 mm
Traction free crack
t--Crack-tip
Fig. 2. Finite element mesh: radially focussed mesh of eightnoded quadrilateral isoparamteric elements with reduced integration (624 elements); crack-tip elements with r - l singularity.
the creep deformation r a t e s (r -N/N+I) are higher than the r -1/2 singularity of the L E F M field. Moreover, the material is incompressible in the creep dominant region. Under plane strain conditions, this may lead to "locking" problems in a displacement-based finite e l e m e n t analysis. Hence, the choice of elements and mesh design is critical. The crack-tip region can be effectively modelled by the "spider web" finite element configuration (Anderson, 1991) which consists of concentric rings of four-sided elements that are focussed toward the crack-tip (Fig. 2). Each element has straight sides with side nodes located at the midpoints between the corner nodes. The innermost rings are degenerated to triangular crack-tip elements as described subsequently. Since the cracktip region contains steep stress and strain gradients, the mesh refinement should be greatest at the crack-tip. The spider web design facilitates a smooth transition from a fine mesh at the cracktip to a coarser mesh remote to the tip. The final mesh used in most of the numerical studies consists of 624 eight-noded isoparametric elements. Fifty-two rings of elements divided into 12 radial segments are focussed at the crack-tip. The element sizes decrease from the boundary towards the crack tip as follows: the outer ten rings of elements have a radial width equal to R/IO each except for the first two which are half this size to provide a transition; a similar pattern is adopted for the next ten rings where R/lO is
273
replaced with R / 1 0 0 ; and this pattern is repeated. The innermost two rings have a radial width of R / 2 0 0 0 0 0 each. This arrangement effectively captures both the crack-tip singularities and the smaller gradients in the field variables away from the tip. The finite element model contains 2001 nodes with 4002 degrees of freedom. For loading rates higher than 5 kPa m 1/2 s-1, the above mesh with R = 10 mm is found to locate the boundary more than 50 x rcn away from the crack-tip• However, for loading rates below 5 kPa m 1/2 s-1, R = 100 mm is required. Ten more rings are added to the previous mesh with the first two rings of radial length 5 mm and the subsequent eight rings of 10 mm length. The total number of elements for this mesh is 744.
5.3. Locking phenomenon The deformation in the creep-dominated zone near the crack-tip is highly incompressible. Nagtegal et al. (1974) have shown that the tangent stiffness finite element solutions exhibit much too stiff a response under constant volume deformation. The cause of the problem was identified in the very formulation of such finite elements and is summarized as follows: From the principle of virtual power a tangentstiffness finite element must satisfy: fsT//~i d s =
Y'.
f
(TijEij dV,
(35)
elements V elements
where T/ is the surface traction rates on the surface area S, d-ij is the stress rate within each element obtained from the prescribed constitutive law in terms of the corresponding current stress (¢rij) and strainrate (gq). Noting that all creep deformations are purely deviatoric, for the idealized case of isotropic homogeneous material, the right-hand side of Eq. (35) can be decomposed into deviatoric and volumetric components according to: f T/~ii d s
=
~] f [~j~j +k(e~k) ] dV, elements JV elements•
2
(36)
274
S. Nanthikesan, S. Shyam Sunder/Mechanics of Materials 21 (1995) 265-279
where ~ij and ~# are the deviatoric stress and strain-rates, respectively; • = E / 3 ( 1 - 2u) is the elastic bulk modulus; and d~k is the volumetric strain-rate. Under creep loading or constant strainrate loading, the steady-state creep condition is reached as t --* oo and the non-negative t e r m sijeij and ~. tend to vanish point wise. From the incremental finite element formulation given by Eq. (36), kkk = 0 point wise in each element. Nagtegal et al. (1974) have shown that this requirement leads to unrealistic kinematic constraints in an assemblage of finite elements. Consequently, attaining the steady-state condition may be delayed or prevented altogether. The result is an overly stiff system, or the so-called "locking" phenomenon. The techniques adopted to eliminate locking in a displacement-based finite element analysis include: (a) restricted mesh topology, (b) limited element shapes, and (c) reduced integration. Nagtegal et al. (1974) show that the zero volume change condition can be met only when the total number of constraints are less than or equal to the total number of degrees of freedom. They find that certain mesh configurations reduce the number of constraints per degree of freedom. For instance, four constant strain triangles arranged to form a quadrilateral was shown to satisfy the incompressibility condition. DeLorenzi and Shih (1977) determined that meshes constructed with rectangular eight-noded isoparametric elements are also capable of representing the zero volume change mode. Quadrilateral shapes for the eightnoded elements were not studied, but elements with curved edges were found to exhibit severe locking. Extending the studies of DeLorenzi and Shih (1977) and Hinton et al. (1978), Dodds (1982) showed that a reduced (2 × 2) Gaussian integration scheme within arbitrarily shaped eight-noded plane strain finite elements eliminates the locking behavior in elastic-plastic materials. Dodds (1982) also observed that reduced integration yields better predictions of limit loads for an identical mesh. Locking can also be eliminated by using mixed variational formulations which retain the
volumetric strain-rate as an independent variable (see, e.g., Bathe, 1982; Nagtegal et al., 1974). The reduced integration method permits simplification of element mesh refinement schemes near a crack-tip by dropping the rectangular shape restrictions. Moreover, element computations decrease by a factor of 9 / 4 compared to mixed variational formulations, and by a factor of 3 compared to the arrangement of four constant strain triangular elements in the form of a quadrilateral. The use of reduced integration on elements with straight edges (i.e., orthogonal isoparametric coordinates) gives rise to the Barlow point property (Barlow, 1989). Because of this, the strains calculated from the interpolation functions display higher accuracy at integration points than anywhere else in the element. This property is especially advantageous for nonlinear constitutive equations of the type considered here. Bassani and McClintock (1981) successfully applied the reduced integration technique for the power-law (incompressible) creep behavior of metals and determined the relaxation of crack-tip stresses. This paper adopts a reduced (2 × 2) Gaussian integration scheme with eight-noded isoparametric elements to model the incompressible transient and viscous creep response of ice. The mesh chosen for the present study performs very well and does not produce spurious modes of deformation which may be encountered when using reduced integration. 5.4. Crack-tip elements
At the crack-tip, the eight-noded quadrilateral elements are degenerated down to triangles by collapsing the three nodes of one side to occupy the same location. The crack-tip elements in elastic problems are often obtained by moving the mid-side nodes to the quarter points and tying the collapsed nodes;, which results in a r -1/= singularity (Barsoum, 1977). Within the creep zone around the crack-tip, the r -N/CN+I) singularity of the R R strain field becomes dominant (Riedel and Rice, 1980). Consequently, elastic singular elements are not appropriate for an elastic-creep analysis. A crack-
275
S. Nanthikesan, S. Shyam Sunder/Mechanics of Materials 21 (1995) 265-279
tip element that is capable of characterizing this higher singularity efficiently is one that is degenerated to a triangle as before, but the crack-tip nodes are untied and the middle nodes on the other three sides remain at mid-point. This elem e n t geometry with mid-point nodes produces a r-~I singularity on the displacement derivatives (Barsoum, 1977; Bassani and McClintock, 1981) which proves to be a good approximation for the weaker singularities of the R R and L E F M fields.
of the elastic and plastic strains, respectively, and is written as:
i1
E = ~- + c~eo o-
(37)
L°0 J
where ~ is a material constant, e 0 is the yield strain, and ~r0 is the yield stress. The asymptotic stress field under small-scale yielding conditions for this material can be written, in polar coordinates (r, 0) with origin at the crack-tip, as:
(J)N+I~ij(ON), 1
6. Validation of numerical model
O'ij = Or0
The ability of the numerical model to accurately predict crack-tip processes is evaluated by comparing the results of numerical simulations with known analytical solutions for the asymptotic fields and the creep zone under a plane strain m o d e I loading configuration. The results should be, and indeed are, independent of the specific choice of the radius of the boundary layer R, when the boundary is chosen sufficiently far from the creep zone. The first numerical example considers a rateindependent plastic material which exhibits classical power-law hardening of the R a m b e r g - O s g o o d type. The uniaxial form of the constitutive law for this material expresses the total strain as the sum
(38)
60.OCOINr
where J denotes the J-integral, which under plane strain SSY conditions is related to the stress intensity factor K I through J = K~(1 - v)e/E; I u is an integration constant (IN= 5.51 for plane strain conditions and N = 3) and 6is(O, N) is an angular function, both of which can be determined from Sih (1983). The above equation represents the well-known H R R field (Hutchinson, 1968; Rice and Rosengren, 1968). This material model, rather than a power-law creeping solid, is used to verify the effectiveness of the r -1 singular tip elements and reduced integration. While the asymptotic fields are selfsimilar (identical) for the two material models, the size of the inelastic process zone is indepen-
(e)
(A) to6
........
I
........
I
......
"1
.......
"T
.......
lOP
........
I
........
I
.......
I
.......
]
.......
J
.......
I
.......
0- 7"
o ,.90 ° I N-3
° l°s
.~ ]
10-1
!
l "~w°l°s
tF'~'''''" to~
1 ~ ~ ~ " ~
~. ~,.
tO"3
to3
io2
........
10"$
!
10.4
........
!
10-3
.......
J
10-2
. . . . .
,d
I0-I
i
,r
.,.,,I
10.
I00
i
5
10"i
i I o''''!
10.4
........
10.3
10.2
J
. . . . .
10-1
i-/a
Fig. 3. Numerical verification of asymptotic fields for elastic power-law hardening plastic material.
1o°
S. Nanthikesan, S. Shyam Sunder/Mechanics of Materials 21 (1995) 265-279
276
io3
' ' ...... I
....... I
........
I
.......
~
.......
0 = 45 ° IO2
I01
"-----2>
I0 o
i 0 "1
...---NUMERICAL = . , ..... I
i0 "2
10-5
i0 --4
. ,,,,.,i
= i rI~ITT I
10-3
1o.2
* ....... I
10-1
.......
100
r/a
normal opening and shear stresses along a 45 ° radial path emanating from the crack-tip is shown in Fig. 4. This is the instantaneous response (at t = 0 ÷) under a suddenly applied stress intensity of 100 KPa m 1/2 that is held constant with time. As expected, the numerical predictions are in excellent agreement with the well-known L E F M field. Further, the numerical estimate of the J-integral (9.43 × 10 .4 J / m 2) is within about 0.6 percent of the theoretical value (9.38× 1 0 . 4 J/m2). Under sustained loading, a creep zone develops around the crack-tip and grows with time.
Fig. 4. N u m e r i c a l v e r i f i c a t i o n of a s y m p t o t i c fields for elastic power-law c r e e p m a t e r i a l at t = 0 +.
(an)
POWER LAW
:o.oz
.KI = lO0 K P L m I ~
dent of time for the plastic material and this is convenient for purposes of evaluation. The material parameters used in the numerical simulation are: E = 9.5 GPa, v = 0.33, d = 1.6 × 10 .6 , ~r0 = 1 MPa, and N = 3. Note that e 0 = ~ro/E. Numerical predictions of the asymptotic fields for equivalent stress and equivalent strain, corresponding to a suddenly applied stress intensity of 100 kPa m 1/2 that is held constant, are shown in Fig. 3 along a 90 ° radial path emanating from the crack tip. For N = 3, the stress exhibits a r - 1 / 4 singularity while the strain exhibits a r - 3 / 4 singularity. The figure shows that the numerical results are in excellent agreement with the theoretical H R R fields. Only at large radial distances approaching the crack length does the stress field begin to follow the prescribed L E F M field. But this corresponds to large-scale plastic deformations in the test material, which are not of direct interest in this study. The second numerical example considers an isotropic material obeying an elastic power-law creep model of deformation. The elastic constants and power-law exponent is as before with the flow stress for freshwater polycrystalline ice at T = - 1 0 ° C obtained from Shyam Sunder and Wu (1989a) as V = 179.8 MPa. The asymptotic response of a typical pair of
i
~
t=lOs
(m)
0.01
-o:oz Traction-Free Crack
(a)
i0"1
pow'[~ LAW CREEP K.I = 100 K P & m 1/2 10.2
T =-10°C
N
M 10.3
10 "1
10 0
101
TIME (s)
(b) Fig. 5. (a) T i m e evolution of c r e e p zone boundary. (b) T i m e versus c r e e p zone size for a e l a s t i c - p o w e r law c r e e p m a t e r i a l u n d e r step loading.
S. Nanthikesan, S. Shyam Sunder~Mechanics of Materials 21 (1995) 265-279
Riedel and Rice (1980) express the boundary of this creep zone rcr(0, t) as follows: 2
~ 1 K 2 [ / ( N + 1) 2 Egot rcr = 2~[ 2 u 8 U + l vN
U--I
Fcr(O),
(39)
where For(O) is an angular function which depends on N and the loading condition, i.e., plane strain versus plane stress. This function, which is independent of time, defines the shape of the creep zone. Figure 5a shows the shape of the creep zone predicted by the numerical model at three time instants: t = 1, 3.3 and 10 seconds. These are in excellent agreement with the functional variation of Fcr(0) plotted in Fig. 1 of Riedel and Rice (1980). The numerical model predicts the size of the creep zone at a time equal to 10 seconds, which is typical of fracture tests and impact problems in ice, to be equal to 1.84 × 10 .2 mm. This corresponds to the maximum value of rc~ which occurs at 0 = 90 °. The Riedel and Rice equation given above yields a size of 1.73 × 10 .2 mm that is within 6.5 percent of the numerical prediction. In evaluating the above equation, it is assumed that a N~ 1 as shown by Riedel and Rice and the maximum value of F¢r = 0.25, which occurs at 0 = 90 ° for N = 3. The evolution of the size of the creep zone with time is plotted in Fig. 5b. The predicted linear variation is in agreement with the Riedel and Rice equation for N = 3.
7. Conclusions A numerical model has been developed to study the interaction between creep deformations and a stationary crack in polycrystalline ice. The problem is analyzed using the physically-based constitutive model of Shyam Sunder and Wu (1989a,b) in conjunction with the finite element method. The solution to the c r e e p - c r a c k interaction problem using this constitutive theory involves the numerical time integration of a coupled set of very stiff and highly nonlinear first-order differential equations. A closed-form N e w t o n - R a p h -
277
son (tangent) formulation in conjunction with an implicit form of the a-method of time integration is developed to solve the local constitutive equations. A consistent constitutive Jacobian matrix, which is used to assemble the global tangent stiffness matrix in a finite element analysis, is also established in closed form. The algorithm is implemented as a user-specified subroutine in the finite element program ABAQUS. The creep-crack interaction problem involving a traction-free semi-infinite crack embedded in a semi-infinite medium is simulated using eightnoded isoparametric elements which are arranged in a cascading, radially focussed mesh. The crack-tip element, formed by collapsing the three nodes on one side of an element, effectively models the singular stress and strain fields ahead of the crack-tip. A reduced Gaussian integration scheme is used to avoid numerical "locking" and increase computational efficiency. The ability of the numerical model to accurately predict crack-tip processes has been verified by comparing numerical predictions of: (a) the asymptotic fields with the well-known L E F M and H R R fields; and (b) the size, shape and time evolution of the creep zone produced by a suddenly applied load in an elastic power-law creep material with the analytical solutions of Riedel and Rice (1980). The companion paper (Nanthikesan and Shyam Sunder, 1995) presents the results from a comprehensive set of numerical simulations that analyze the role of transient creep, including loading history, temperature sensitivity, material strain hardening and fabric anisotropy, in predicting the size, shape and the time evolution of the creep zone in polycrystalline ice.
Acknowledgements The authors would like to thank Professor David M. Parks of MIT for many helpful discussions. The authors acknowledge financial support from the Office of Naval Research (Grant No. N00014-92-J-1208), the Minerals Management Service, and from an industry consortium consist-
278
S. Nanthikesan, S. Shyam Sunder/Mechanics of Materials 21 (1995) 265-279
ing of AMOCO, ARCO, BP America, CHEVRON, CONOCO, EXXON, MOBIL and TEXACO through the MIT Center for Scientific Excellence in Offshore Engineering.
References ABAQUS User's Manual, Version 4.7 (1988), Hibbit, Karlsson and Sorensen Inc., Providence, Rhode Island. Argyris, J.H., L.E. Vaz and K.J. William (1978), Improved solution methods for inelastic rate problems, Comp. Meth. App. Mech. Eng. 16, 231. Anderson, T.L. (1991), Computational Fracture Mechanics, CRC Press, p. 659. Barlow, J. (1989), More on optimal stress points - reduced integration, element distortions and error estimation, Int. J. Num. Meth. Eng. 28, 1487. Barsoum, R.S. (1977), Triangular quarter-point elements as elastic and perfectly-plastic crack tip elements, Int. J. Num. Meth. Eng. I1, 85. Bassani, J.L. and F.A. McClintock (1981), Creep relaxation of stress around a crack tip, Int..L Solids Struc. 17, 479. Bathe, K.J. (1982), Finite Element Procedures in Engineering Analysis, Prentice Hall, New Jersey. DeLorenzi, H.G. and C.F. Shih (1977), On the use of 2D isoparametric elements for calculations in the fully plastic range; Int. J. Fract. 13, 507. Dodds, R.H. (1982), Effects of reduced integration on the 2-D quadratic isoparametric element in plain strain plasticity. Int. J. Fract. Mech. 9, R75. Glen, J.W. (1955), The creep of polycrystalline ice, Proc. Roy. Soc. London A, 228 (1175), 519. Goodman, J.D. (1980), The Fracture Toughness of Ice, Data Report., Hokkaido University, Dept. of Applied Physics, Sapporo, Japan. Hinton, E., E.M. Salonen and N. Bicanic (1978), The Mathematics of Finite Elements and Applications, III, p. 437. Hutchinson, J.W. (1968), Plastic stress and strain fields at a crack tip, J. Mech. Phys. Solids 16 (13), 337. Lush, A.M., G. Weber and L. Anand (1989), An implicit time-integration procedure for a set of internal variable constitutive equations for isotropic elasto-viscoplasticity, J. Plasticity 5, 521. Nagtegal, J.C., D.M. Parks and J.R. Rice (1974), On numerically accurate finite element solutions in the fully plastic range, Comp. Meth. App. Mech. Eng. 4, 153. Nanthikesan, S: and S. Shyam Sunder (1995), Tensile cracks
in polycrystalline ice under transient creep: Part II Numerical Simulations, Mech. Mater 21,281. Nixon, W.A. and E.M. Schulson (1987), A micromechanical view of the fracture toughness of ice, J. Phys. Colloq. (France) 48, C1-313. Palmer, A.C., D.J. Goodman, M.F. Ashby, A.G. Evans, J.W. Hutchinson and A.R.S. Ponter (1983), Fracture and its role in determining ice forces on offshore structures, Annals Glac. 4, 216. Rice, J.A. (1968), A path independent integral and the approximate analysis of strain concentration on notches and cracks, J. App. Mech. Trans. ASME 6, 379. Rice, J.A., and G.F. Rosengren (1968), Plane strain deformation near a crack tip in a power-law hardening material, J. Mech. Phys. Solids 16, 1. Riedel, H. (1981), Creep deformations at crack tips in elasticviscoplastic solids, J. Mech. Phy. Solids 29, 35. Riedel, H. and J.R. Rice (1980), Tensile cracks in creeping solids, Fracture Mechanics: 12th conference, A S T M STP 700, p. 112. Schulson, E.M. (1990), The brittle compressive fracture of ice, Acta Met Mat. 38 (10), 1963. Schulson, E.M., S.G. Hoxie and W.A. Nixon (1989), The tensile strength of cracked ice, Philos. Mug. 59 (2), 303. Shyam Sunder, S. and M.S. Wu (1989a), A differential flow model for potycrystalline ice, Cold Reg. Sci. Technol. 16 (1), 45. Shyam Sunder, S. and M.S. Wu (1989b), A multiaxial differential mode[ for flow in orthotropic polycrystalline ice, Cold Reg. Sci. Technol. 16 (2), 223. Shyam Sunder, S. and M.S. Wu (1990), On the constitutive modeling of transient creep in polycrystalline ice. CoM Reg. Sci. Technol. 18, 267. Shyam Sunder, S., A. Elvin and S. Nanthikesan (1993), Numerical modeling of transient creep in polycrystalline ice, .L Engr. Mech. 119 (10) 2011. Sih, C.F. (1983), Small scale yielding analysis of mixed mode plane strain crack problem, ASTM STP 560, 187. Smith, D.J. (1990), Crack growth in brittle materials exhibiting primary creep. Eng. Fract. Mech. 37 (1), 27. Timco, G.W. and R.M.W. Frederking (1986), The effects of anisotropy and microcracks on the fracture toughness of freshwater ice, Proc. Fifth Int. Sym. Offshore Mechanics and Arctic Engineering, Tokyo, Japan, Vol. IV, p. 341. Urabe, N., and A. Yoshitake (1981), Strain-rate dependent fracture toughness of pure and sea ice, Proc. 6th Int. lAHR Sym. Ice, Quebec City, Canada, July 27-31, p. 551. William, K.J. (1978), Numerical solution of inelastic rate process, Comp. Struct. 8, 511. Wu, M.S. and S. Shyam Sunder (1992), On the constitutive modeling of transient creep in polycrystalline ice: Closure to Discussion. Cold Reg. Sci. Tech. 20 (3), 315.
S. Nanthikesan, S. Shyam Sunder~Mechanics of Materials 21 (1995) 265-279
Appendix G Matrix: a1 3
( a 1 + a3) 3
a3 3
(al +a2) 3 G=
a2 3 ( a 2 + a3) 3 2a4 Sym.
2a5 2a 6
H Matrix: 3 ( a I + a3) a. 2 a2
3ala2a3 a, 2 3(a1+a2) a*2
3ala2a 3 a*2 3alaea 3 a*2 3(a2 + a3) a*2
H=
a2 3
1 2a4 Sym.
1 2a5
1 2a 6 where a * = (a2a 3 + a3a I + a~a2 )2.
279
Mechanics of Materials 21 (1995) 265-279
Tensile cracks in polycrystalline ice under transient creep: Part I - Numerical modeling S. Nanthikesan a,1, S. Shyam Sunder h,2 a Prudent Engineering Group, Liverpool, IVY 13088, USA b Department of Civil and Environmental Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139, USA
Received 30 August 1993; revised version received 1 June 1995
Abstract The interaction between creep deformations and a stationary or growing crack is a fundamental problem in ice mechanics. Knowledge concerning the physical mechanisms governing this interaction is necessary: (1) to establish the conditions under which linear elastic fracture mechanics can be applied in problems ranging from ice-structure interaction to fracture toughness testing; and (2) to predict the ductile-to-brittle transition in the mechanical behavior of ice and, especially, the stability and growth of cracks subjected to crack-tip blunting by creep deformations. This requires a quantitative estimate of the creep zone surrounding a crack-tip, i.e., the zone within which creep strains are greater than the elastic strains. The prediction of the creep zone in previous ice mechanics studies is based on the theory developed by Riedel and Rice (1980) for tensile cracks in creeping solids. This theory is valid for a stationary crack embedded in an isotropic material obeying an elastic, power-law creep model of deformation and for a suddenly applied uniform far-field tension road that is held constant with time. The deformation of ice at strain-rates ahead of a crack (i.e., 10 -6 to 10 -2 s -1) is dominated, however, by transient (not steady power-law) creep and the loading, in general, is not instantaneous and constant. A numerical model is developed in this paper to investigate the role of transient creep and related physical mechanisms in predicting the size, shape and time evolution of the creep zone surrounding the tip of a static crack in polycrystalline ice. The model is based on the fully consistent tangent formulation derived in closed form (Shyam Sunder et al., 1993) and used in the solution of the physically-based constitutive theory developed by Shyam Sunder and Wu (1989a, b) for the multiaxial behavior of ice undergoing transient creep. The boundary value problem involving incompressible deformations ahead of a stationary, traction-free mode I crack in a semi-infinite medium is modeled and solved by a finite element analysis using the boundary layer approach of Rice (1968). This model is verified by comparing its predictions with (i) the known theoretical solutions for the elastic and HRR asymptotic stress and strain fields in an elastic-plastic material of the Ramberg-Osgood type, and (ii) the creep zone size for an isotropic material obeying the elastic power-law creep model of deformation. Keywords: Fracture mechanics; Ice mechanics; Creep fracture; Fracture process zone; Crack-tip inelasticity;
Polycrystalline ice; Fracture mechanics of ice; Prediction of crack-tip creep zone
i Senior Engineer, Prudent Engineering Group, Liverpool, NY 13088, USA. Formerly Postdoctoral Associate, Department of Civil and Environmental Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139, USA.
2 Senior Research Scientist, Department of Civil and Environmental Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139, USA. Now at the Building and Fire Research Laboratory, National Institute of Standards and Technology, Gaithersburg, MD 20899, USA.
0167-6636/95/$09.50 © 1995 Elsevier Science B.V. All rights reserved SSD1 0167-6636(95)00014-3
266
S. Nanthikesan, S. Shyam Sunder~Mechanics of Materials 21 (1995) 265-279
1. Introduction
The interaction between creep deformations and a stationary or growing crack is a fundamental problem in ice mechanics. Knowledge concerning the physical mechanisms governing this interaction is necessary: (1) to establish the conditions under which linear elastic fracture mechanics (LEFM) can be applied in problems ranging from ice-structure interaction to fracture toughness testing; and (2) to predict the ductile-to-brittle transition in the mechanical behavior of ice and, especially, the stability and growth of cracks subjected to crack-tip blunting by creep deformations. A qualitative discussion regarding the applicability of L E F M to predict ice forces on offshore structures is provided by Palmer et al. (1983). They characterize the modes of deformation and, in particular, the relative role of creep and fracture processes as a function of loading rate and the problem geometry. For impact times of less than about 10-100 seconds, this work suggests that creep deformations are often negligible and only (linear elastic) fracture processes dominate. The need to analyze the creep-crack interaction problem more rigorously has been recognized by experimentalists interested in characterizing the fracture resistance (i.e., "toughness") of polycrystalline ice. The inelastic (creep-dominated) fracture process zone ahead of a crack must be small relative to the key geometric dimensions of a test specimen to ensure small-scale yielding (SSY) conditions and the validity of LEFM. This requirement is explicitly considered by Goodman (1980), Urabe and Yoshitake (1981), Nixon and Schulson (1987) and Timco and Frederking (1986). More recently, Schulson (1990) predicts the strain-rate at which the ductile-to-brittle transition occurs in polycrystalline ice under uniaxial compression. His wing crack model of brittle compressive failure assumes that a crack can propagate only if the crack-tip creep zone is smaller than a critical size that is proportional to the microstructural scale of the aggregate. Otherwise, there is no growth and the crack remains
stable. A similar argument is used by Schulson et al. (1989) to explain the stability of long nonpropagated cracks under uniaxial tension. The model allows for an increase in fracture toughness due to crack-tip blunting; more specifically, it assumes that, as in steel, the toughness is proportional to the square-root of the crack-tip radius and which, in turn, is proportional to the "radius" of the creep zone. The prediction of the creep zone, defined as the zone surrounding the crack-tip within which creep strains are greater than elastic strains, in the above studies is based on the theory developed by Riedel and Rice (1980) for tensile cracks in creeping solids. Riedel and Rice present approximate analytical solutions for the creepdominated near-tip asymptotic stress and strain fields (the RR fields) under small-scale yielding conditions as well as for the boundary of the creep zone. The solutions are valid for a stationary crack embedded in an isotropic material obeying an elastic, power-law creep model of deformation and one that is subjected to a suddenly applied uniform far-field tension load which is held constant with time. These theoretical predictions were verified by Bassani and McClintock (1981) using finite element analysis. Riedel (1981) and Smith (1990) have extended the theory to include the effects of isotropic strain hardening during creep. The deformation of ice at strain-rates of engineering interest (i.e., 10 -6 to 10 -2 s -a) is dominated, however, by transient (not steady powerlaw) creep and the loading, in general, is not instantaneous and constant. The strain-rates in the vicinity of a crack fall within this range (as will be shown subsequently). Consequently, any analysis of the inelastic processes surrounding the crack-tip must include the effects of transient creep which involves isotropic and kinematic hardening. The analysis must also consider the effects of loading history, temperature sensitivity, variation in material strain hardening among different types of ice, and fabric anisotropy. This paper presents a numerical model to study the role of transient creep and related physical mechanisms in predicting the size, shape and
S. Nanthikesan, S. Shyam Sunder~Mechanics of Materials 21 (1995) 265-279
time evolution of the creep zone surrounding a crack-tip in polycrystalline ice. The physicallybased constitutive theory of Shyam Sunder and Wu (1989a,b) for the multiaxial creep of ice undergoing transient and viscous creep is used in conjunction with a finite element analysis to model the creep-crack interaction problem. The solution to the c r e e p - c r a c k interaction problem using this constitutive theory involves the numerical time integration of a coupled set of very stiff and highly nonlinear first-order differential equations. These equations, governing the evolution of the observable and internal state variables, characterize local material behavior at each integration point in a finite element mesh. A closed-form N e w t o n - R a p h s o n (tangent) formulation in conjunction with an implicit form of the o~-method of time integration is developed to solve the local constitutive equations. A consistent constitutive Jacobian matrix, which is used to assemble the global tangent stiffness matrix in a finite element analysis, is also established in closed form. The algorithm is implemented as a user-specified subroutine in the finite element program A B A Q U S and its numerical stability, accuracy and convergence properties have been verified by Shyam Sunder et al. (1993). The problem of a stationary, traction-free, tensile (mode I) crack is modeled using the equivalent boundary layer idealization proposed by Rice (1968). The crack-tip elements and the mesh configuration are chosen so as to model the singular stress and strain fields within the boundary layer accurately. Finally, the ability of the numerical model to predict the elastic and H R R asymptotic stress and strain fields and the creep zone for an isotropic material obeying a steady power-law creep model of deformation is validated against known analytical solutions. The companion paper (Nanthikesan and Shyam Sunder, 1995) presents the results from a comprehensive set of numerical simulations that analyze the role of transient creep, including loading history, temperature sensitivity, material strain hardening and fabric anisotropy, in predicting the size, shape and the time evolution of the creep zone in polycrystalline ice.
267
2. Summary of the constitutive model
A summary of the physically-based constitutive model of creep in polycrystalline ice developed by Shyam Sunder and Wu (1989a, b) is presented in this section. The mathematical formulation of the transient creep deformation mechanism, founded on the theory of internal state variables, is characterized by: (a) a scalar internal state variable B, the drag stress, to describe isotropic hardening; (b) a second-rank tensorial internal state variable R, the elastic back stress tensor, to describe kinematic hardening or deformation-induced anisotropy; and (c) a generalized equivalent stress measure to describe fabric anisotropy of the Hill-type in the solid. In addition to transient creep, the model takes into account the elastic and steady viscous creep deformations through a summation of the strain-rates for each of the deformation mechanisms. Transient creep in polycrystalline ice is anelastic and fully recoverable upon unloading while viscous creep is irrecoverable. Only small deformations are considered since the magnitude of the transient creep strain does not exceed 1%. At larger strains, the material exhibits strain-softening by progressive distributed microcracking or recrystalization processes if premature failure by brittle fracture does not occur. Readers are referred to the previously referenced papers, as well as Shyam Sunder and Wu (1990) and Wu and Shyam Sunder (1992), for the physical basis of the constitutive model and its material parameters. Matrix notation is used in this paper. A superposed dot denotes a time derivative and the subscript " e q " refers to a scalar equivalent value. For small deformation theory, the total strainrate, i, is additively decomposed as follows: = ~e + iv + ~t, (1) where the subscripts "e", " v " and " t " denote the elastic, viscous and transient components of the strain-rate, respectively. The elastic strain-rate is related to the stress-rate, d-, by:
b=D~e=O(~-iv-~,) where D is the elastic stiffness matrix.
(2)
268
S. Nanthikesan, S. Shyam Sunder / Mechanics of Materials 21 (1995) 265-279
The viscous strain-rate vector, iv, follows the classic power law (Glen, 1955). For an incompressible orthotropic material it is given by: • -iv = ~'0/3
N-I 0"eq S,
(3)
S = Go', (4) where ~0 is a scalar reference strain-rate and can be taken as unity with no loss of generality; /3 = a 1 + a2, where a I and a 2 (as well as a3, a 4, a 5 and a 6) are constants which describe the Hilltype fabric anisotropy of the material (these orthotropic parameters are calculated from experimental data, see Shyam Sunder et a1.(1993)), N is the power law exponent; cr is the stress vector {0-xx, 0-yy, o'~z, cr~y, 0-.~, 0-~y}T; S is the deviatoric stress vector; G is the matrix which transforms stresses into their "deviatoric" components and is given in the Appendix; V is the temperature-dependent flow stress in the reference direction (assumed to be the Y-direction in this paper) given by the Arrhenius relationship:
where V0 is a temperature-independent constant, Q is the activation energy, T is the temperature in Kelvin, and R is the universal gas constant; 0-eq is the equivalent stress measure given by: 3 T o'e2q= ~0- GO-,
(6)
where the superscript " T " denotes the transpose operator. The transient strain-rate vector, it, is given by: ,,
(7)
with Sa = G0-~ = G(0- - R). (8) The evolution equations for the back stress vector, R, and the scalar drag stress B are given by: 1 = ~ [ 3 A E H : i t,
(9)
Be /~ = O [ - - ] ~ t
I_O'o,o ]
eq,
'
(tO)
where 0-d is the reduced stress vector, A E is the anelastic modulus of kinematic hardening; E is the Young's modulus in the reference (Y) direction, and A is a non-dimensional material constant (in a uniaxial creep test, the ratio of the maximum transient creep strain to elastic strain is equal to 1//1); H is the strain (and strain-rate) transformation matrix which transforms the strains into their "deviatoric" components and is given in the appendix; B 0 is a non-dimensional material constant which represents the threshold value of the drag stress; O is a function of the drag stress which ensures that upon unloading (when magnitude of transient strain decreases), B tends asymptotically to the threshold drag stress, B~), and is given in Shyam Sunder and Wu (1990); H is a non-dimensional material constant which describes the rate of isotropic hardening; 0-d,eq is the equivalent measure of reduced stress obtained by: 3 O'2d,eq = ~-0"J60"d"
(11)
gt,eq is the equivalent transient strain-rate and is given by: ~2t,eq = ~ i T H e t ,
(12)
here, the strain-vector is expressed as {exx, e-yy, ezz, 7xy, Yxz, y~y}v. Note that the shear components, 7o, are engineering strains. 2.1. Material parameters
The constitutive model consists of the standard elastic moduli, the universal gas constant R, the flow stress V0, the power-law exponent N, and the Arrhenius activation energy Q. The parameters R, Q, V0 and N are associated with viscous creep behavior• The transient creep is characterized by three additional constants which model strain-hardening: (1) the anelastic constant for kinematic hardening, A, (2) the initial or threshold drag stress, B 0, and (3) the rate of isotropic hardening, H. In addition, there are the five independent Hill constants, a i (i = 1, 5), which characterize the fabric anisotropy of the material•
S. Nanthikesan, S. Shyarn Sunder~Mechanics of Materials 21 (1995) 265-279
The formulation of transient strain-rate given by Eqs. (8) through (11) considers the flow law to be the same for both deformation-rate mechanisms, the only difference being the existence of drag and back stresses in the case of transient flow. It also implies the following: (a) the exponent N is the same as that for viscous flow in Eqs. (3) through (7), and (b) the temperature dependence of the transient creep-rate, represented by the parameter V, is given by an Arrhenius law with an activation energy equal to that of viscous creep. These assumptions are justified in Shyam Sunder and Wu (1989).
3. Integration of constitutive equations To solve creep-crack interaction problems using the finite element method requires a numerical integration of the governing nonlinear, ratedependent constitutive equations summarized above. When the deformation is dominated by creep, as is the case for the problem under consideration, the material behavior is incompressible and displays significant strain hardening. The coupled set of first-order differential equations, which models the local constitutive response at each integration point in a finite element mesh, is very stiff and highly nonlinear. This complexity, in turn, poses severe demands on the global finite element formulation with respect to numerical stability, accuracy and rate of convergence. A comprehensive numerical model to solve ice mechanics problems dominated by transient creep deformations using the finite element method has been developed by Shyam Sunder et al. (1993) and is summarized in this section.
3.1. 77me integration scheme In a displacement-based finite element formulation, estimates of the incremental displacement fields (and thus, the incremental strain fields) are known at the beginning of each constitutive iteration within a time increment. The problem of marching the solution forward in time then reduces to computing the increment in stress (and in the internal state variables) corresponding to
269
the given increment in total strain. The constitutive equations are solved using the one-parameter class of integration algorithms known as the amethod of integration, with:
F~i~ = (1 - a ) t F + a'+~'F
(13)
where the weighting parameter c~ may vary between 0 and 1; in general, the symbol F represents the (Cauchy) stress vector, the backstress vector, and the non-dimensional drag stress; F~(i) denotes the weighted value of F between time t and t + A t after the i th iteration. The two superscripts on the left side of symbol refers to time; and the superscript in parentheses on the right side refers to the iteration number within the time step. If the latter superscript is missing, it refers to the converged value. This nomenclature will be used for the rest of this paper.
3.2. Newton-Raphson iteration for the stress increment The constitutive equations in differential form are expressed in terms of finite differences to solve for the stress increments. Eqs. (1) through (12) in conjunction with Eq. (13) will yield the solution for the increment in stress t+~t~r~i+1). Such a scheme has only linear convergence. In this paper, a summary of the algorithm based on the N e w t o n - R a p h s o n iteration scheme is presented. The rate of convergence is quadratic in the neighborhood of the actual stress t+a/~r. Rewriting Eq. (2) as a residual, and noting A denotes the increment from time t to time t + At: f(A,~)
= zXo, - D z X ~ e
= Ao--D(A~-
A e v - A~t).
(14)
The total strain increment AE is prescribed, and Aev, Aet are determined assuming a Ao~. The objective is to find the root of f ( A c r ) = 0 which provides an improved estimate of the stress increment. Expanding f(Ao-) in a Taylor series, considering only the first two (i.e., linear) terms, and rearranging yields the following equation: ,+~,~,i+,, =,+~,~,)_~,i,[CA~,,)_ A4i)],
(15)
270
S. Nanthikesan, S. Shyam Sunder~Mechanics of Materials 21 (1995)265-279
where C is the compliance matrix and /)(0 is the elastic-creep stiffness matrix: /~(i)=
C+a
~aAo_~i + ~ } j
(16)
The derivation of the expressions for the partial derivatives in Eq. (16) is given subsequently. 3.3. Convergence criterion
Convergence is considered to be achieved when the magnitude of change in stress increments between two successive iterations at time t + At, is smaller than a prescribed fraction (e.g. 1/1000) of the magnitude of stress at time t + At, i.e., Ao-(i +1) _ Ao-(i)
.< 0.001.
ture have established the accuracy of the numerical model in simulating the constitutive behavior of polycrystalline ice (Shyam Sunder et al., 1993).
4. Derivation of tangent stiffness matrices
This section summarizes the derivation of the fully consistent closed-form solution for the following: (a) the elastic-creep stiffness matrix used in the local N e w t o n - R a p h s o n iteration for solving the material constitutive equations, and (b) the constitutive Jacobian matrix used in the global iterations. The latter is required in any "implicit" finite element analysis using a Newton-Raphson iteration scheme to solve the weak form of the linear momentum balance equations.
(17) 4.1. Elastic-creep stiffness matrix
3.4. Numerical implementation and evaluation
The above constitutive model has been implemented as a user-specified subroutine, UMAT, in the finite element program ABAQUS (1988). This subroutine is called at each integration point in the finite element model for every global iteration. At the beginning of each global iteration, U M A T requires the input of the total strain increment and the state variables (stress vector, scalar drag stress and vector back stress). At the end of the iteration, U M A T updates and outputs these state variables, the transient strains and the transient strain-rates, as well as the constitutive Jacobian matrix. Only a fixed (i.e., constant) time increment can be specified by the user. Therefore, the automatic time incrementation scheme of A B A Q U S is used in this analysis. In this scheme, the time step is increased by a factor of 1.5 if global convergence is attained in less than half the user-specified number of iterations in two consecutive time increments. On the other hand, if convergence is not achieved within this prescribed number of iterations, the time increment is reduced by a factor of four. Extensive comparisons with experimental data and theoretical solutions published in the litera-
Rewriting Eq. (16) as: I) (i, = [ C q- F(vi) + F, (') ] -1
18)
where 19)
F (i) ~ a aAo.~i) ,
and
aAe{ i) (20)
F[ i) ~ a aAo.~i ) .
The terms F~(i) and Ft(i) contain the contributions of viscous and transient creep, respectively, to the elastic-creep stiffness matrix. The expressions for these terms are obtained using the product rule of differentiation and then the chain rule, and are summarized below: Fv{/) = g0At [A(/)G 9 (NJr- /~2 ~ V
1) [_(0 "~N-3c(i)[c(i)'~T] \~'a,eq} °or I.°a ) ] (21)
and F[ i , = [ M ~ )1 -'F~O,
(22)
S. Nanthikesan, S. Shyarn Sunder/Mechanics of Materials 21 (1995) 265-279
where
M~' = [I + ½fiAEF~i)H].
(23)
The discretized weak form of the equilibrium equations at time t + At can be written as: fu6ffT t+Ato" dT/=
I is the identity matrix and
F~0 = a[ M~'a)] -it ~,,i, + M(j)),
271
£6uT t+atf B d ~
(24)
+ fg&,T t+zXtfS d ~ ,
(28)
in which
M(j) = goAt I A(~IG k
9
(N-
1) t - t 0
+ 132 (B2i v)
~N-3S~O1' S(e.,,J 0 ~T]]
Wa .... (25)
M~,;)=i+(goAt)2I
3
[
N
(mB(i) 1
\---t,eq]
[X(i)'~2q(i) (~,(i) T H ] ' X ~,'~d ] ~d,oe\~d,a)
lift(i) = *"B2
3 coAt o~~
N
(26)
f//v~TJAE d~'F+ fv~ET to" d ~ ,
(29)
where 0Ao"
( AB (o
[,~_(i) "12 7 ~ d ,c~,eq ] Bct
X ),(i)K,(i {~(i) ,~T] .,d ~.d,.t ~.a,.j ].
where 66 is the variation in the strain field corresponding to the variation in the displacements 6u; the superscript "T", denotes the transpose operator, ,+~fB and ,+~tfs are the body forces and surface tractions, respectively; ~'is the volume and Y' is the surface area. Applied concentrated forces are included as a part of surface tractions in Eq. (28). The constitutive Jacobian matrix, J, is required when the first term in Eq. (28) is linearized around 'u. The linearized term is:
(27)
4.2. Constitutive Jacobian matrix The rate of convergence of "implicit" finite element procedures for nonlinear materials is quadratic when a N e w t o n - R a p h s o n type iterative method is used to solve the weak form of the equilibrium equations. This method exhibits a stronger rate of convergence than the initial stress method, the modified Newton iteration and the quasi-Newton scheme (William, 1978; Argyris et al., 1978; Bathe, 1982). However, the N e w t o n Raphson scheme requires the constitutive Jacobian matrix to be computed at every integration point at each global iteration. The computational effort is considerably reduced by using closedform expressions for the constitutive Jacobian matrix (see, e.g., Lush et al., 1989). This section outlines the key steps in deriving such an expression.
j =
(30)
tK = ~._~,Tj,.~ d ~ ' ,
(31)
0Ae The global Newton iteration scheme is obtained from Eqs. (28),(29), as described in detail in Shyam Sunder et al. (1993) and the stiffness matrix, tK, is written as:
where the B matrix relates the strain fields to the nodal strains fi, as given by • =.~'t~. The stiffness matrix for the system is assembled from the contributions of each finite element. The use of Gaussian quadrature to perform this integration necessitates the computation of J at each integration point within the element. The J matrix for the material model used here is derived by differentiating the difference form of Eq. (2) and Eq. (30). This yields: [ J =D
I
0A'v 0Ae
0A't ] 0Ae "
(32)
Using the chain rule of differentiation as in the previous section, the expression for J is obtained as; J = [C + F v + F t ] - '
(33)
S. Nanthikesan, S. Shyam Sunder/Mechanics of Materials 21 (1995) 265-279
272
It can be seen that this is identical elastic-creep stiffness matrix in matrices bv and F t are given by (20), respectively, when the lower has converged.
in form to the Eq. (18). The Eqs. (19) and level iteration
(~ :
UnfformF~r-FieidTensLieS~'~ss
5. Crack-tip model under small-scale creep
5.1. Formulation of the boundary value problem In this paper, the boundary value problem of a static crack under mode I loading is considered. To obtain the response in the vicinity of the crack-tip, an equivalent but simpler problem of a boundary layer (Rice, 1968) is postulated and solved. U n d e r small scale creep (SSC) conditions, the L E F M stresses dominate at distances from the crack-tip that are large compared to the creep zone size but still small compared to the dimensions of the crack, a. In this simplified model, the actual boundary is replaced by one that is circular and centered around the crack-tip with the actual boundary conditions replaced by the L E F M asymptotic field along this circular boundary. The circular boundary is chosen sufficiently far from the crack-tip where L E F M conditions prevail. This model focuses on the immediate vicinity of the crack-tip and hence is an efficient and convenient idealization to model the asymptotic stress and strain fields in the creep-crack interaction problems of interest. For the purposes of the finite element analysis, the simplified configuration is represented by the crack -R<_xl <0, and the circular boundary centered at x 1 = x 2 = 0 with radius R (see Fig. 1). Here, R is chosen large enough (a >> R >> Rcr) so that the deformation at the boundary is unperturbed by the creep zone. The crack is traction free and the L E F M displacement field, corresponding to the applied load on the actual configuration, is applied at the circular boundary. This displacement field is given by: ui = ~
fi(O)
with i=
1, 2,
fl(O) = (3 - 4u - cos O) cos(O/2), f2(O) = (3 - 4 u - cos O) s i n ( O / 2 ) ,
0" Fig. 1. Schematic of creep-crack interaction problem for a tensile mode I crack in a semi-infinite medium modeled using the bounday layer approach of Rice (1968).
where Kj is the stress intensity factor determined from the actual configuration; G is the shear modulus; u is the Poisson's ratio; and 0 is the angle sustained with the x 1 axis in the counter clockwise direction. Symmetry of this mode I configuration allows the analysis of just a halfplane lying either above or below the plane of the crack, with the boundary conditions on x 2 = 0,x 1 > 0 given by u 2 = o-12 0. U p o n application of the load, the instantaneous response of the material is purely elastic and the crack-tip field is given by LEFM. Sustained application of the load, smaller than that required for crack propagation, gives rise to creep deformations that relax the stresses. The relaxation is most significant in the immediate vicinity of the crack-tip where the strain-rate (and stress) field is singular. This creep-dominated zone grows with time and may eventually occupy the entire solid. During growth, the creep zone is bounded by the L E F M field. =
5.2. Finite element mesh design (34)
Singular stress and strain fields exist within the boundary layer. The singularities associated with
S. Nanthikesan, S. Shyam Sunder~Mechanics of Materials 21 (1995) 265-279
/ 1 0 mm
Traction free crack
t--Crack-tip
Fig. 2. Finite element mesh: radially focussed mesh of eightnoded quadrilateral isoparamteric elements with reduced integration (624 elements); crack-tip elements with r - l singularity.
the creep deformation r a t e s (r -N/N+I) are higher than the r -1/2 singularity of the L E F M field. Moreover, the material is incompressible in the creep dominant region. Under plane strain conditions, this may lead to "locking" problems in a displacement-based finite e l e m e n t analysis. Hence, the choice of elements and mesh design is critical. The crack-tip region can be effectively modelled by the "spider web" finite element configuration (Anderson, 1991) which consists of concentric rings of four-sided elements that are focussed toward the crack-tip (Fig. 2). Each element has straight sides with side nodes located at the midpoints between the corner nodes. The innermost rings are degenerated to triangular crack-tip elements as described subsequently. Since the cracktip region contains steep stress and strain gradients, the mesh refinement should be greatest at the crack-tip. The spider web design facilitates a smooth transition from a fine mesh at the cracktip to a coarser mesh remote to the tip. The final mesh used in most of the numerical studies consists of 624 eight-noded isoparametric elements. Fifty-two rings of elements divided into 12 radial segments are focussed at the crack-tip. The element sizes decrease from the boundary towards the crack tip as follows: the outer ten rings of elements have a radial width equal to R/IO each except for the first two which are half this size to provide a transition; a similar pattern is adopted for the next ten rings where R/lO is
273
replaced with R / 1 0 0 ; and this pattern is repeated. The innermost two rings have a radial width of R / 2 0 0 0 0 0 each. This arrangement effectively captures both the crack-tip singularities and the smaller gradients in the field variables away from the tip. The finite element model contains 2001 nodes with 4002 degrees of freedom. For loading rates higher than 5 kPa m 1/2 s-1, the above mesh with R = 10 mm is found to locate the boundary more than 50 x rcn away from the crack-tip• However, for loading rates below 5 kPa m 1/2 s-1, R = 100 mm is required. Ten more rings are added to the previous mesh with the first two rings of radial length 5 mm and the subsequent eight rings of 10 mm length. The total number of elements for this mesh is 744.
5.3. Locking phenomenon The deformation in the creep-dominated zone near the crack-tip is highly incompressible. Nagtegal et al. (1974) have shown that the tangent stiffness finite element solutions exhibit much too stiff a response under constant volume deformation. The cause of the problem was identified in the very formulation of such finite elements and is summarized as follows: From the principle of virtual power a tangentstiffness finite element must satisfy: fsT//~i d s =
Y'.
f
(TijEij dV,
(35)
elements V elements
where T/ is the surface traction rates on the surface area S, d-ij is the stress rate within each element obtained from the prescribed constitutive law in terms of the corresponding current stress (¢rij) and strainrate (gq). Noting that all creep deformations are purely deviatoric, for the idealized case of isotropic homogeneous material, the right-hand side of Eq. (35) can be decomposed into deviatoric and volumetric components according to: f T/~ii d s
=
~] f [~j~j +k(e~k) ] dV, elements JV elements•
2
(36)
274
S. Nanthikesan, S. Shyam Sunder/Mechanics of Materials 21 (1995) 265-279
where ~ij and ~# are the deviatoric stress and strain-rates, respectively; • = E / 3 ( 1 - 2u) is the elastic bulk modulus; and d~k is the volumetric strain-rate. Under creep loading or constant strainrate loading, the steady-state creep condition is reached as t --* oo and the non-negative t e r m sijeij and ~. tend to vanish point wise. From the incremental finite element formulation given by Eq. (36), kkk = 0 point wise in each element. Nagtegal et al. (1974) have shown that this requirement leads to unrealistic kinematic constraints in an assemblage of finite elements. Consequently, attaining the steady-state condition may be delayed or prevented altogether. The result is an overly stiff system, or the so-called "locking" phenomenon. The techniques adopted to eliminate locking in a displacement-based finite element analysis include: (a) restricted mesh topology, (b) limited element shapes, and (c) reduced integration. Nagtegal et al. (1974) show that the zero volume change condition can be met only when the total number of constraints are less than or equal to the total number of degrees of freedom. They find that certain mesh configurations reduce the number of constraints per degree of freedom. For instance, four constant strain triangles arranged to form a quadrilateral was shown to satisfy the incompressibility condition. DeLorenzi and Shih (1977) determined that meshes constructed with rectangular eight-noded isoparametric elements are also capable of representing the zero volume change mode. Quadrilateral shapes for the eightnoded elements were not studied, but elements with curved edges were found to exhibit severe locking. Extending the studies of DeLorenzi and Shih (1977) and Hinton et al. (1978), Dodds (1982) showed that a reduced (2 × 2) Gaussian integration scheme within arbitrarily shaped eight-noded plane strain finite elements eliminates the locking behavior in elastic-plastic materials. Dodds (1982) also observed that reduced integration yields better predictions of limit loads for an identical mesh. Locking can also be eliminated by using mixed variational formulations which retain the
volumetric strain-rate as an independent variable (see, e.g., Bathe, 1982; Nagtegal et al., 1974). The reduced integration method permits simplification of element mesh refinement schemes near a crack-tip by dropping the rectangular shape restrictions. Moreover, element computations decrease by a factor of 9 / 4 compared to mixed variational formulations, and by a factor of 3 compared to the arrangement of four constant strain triangular elements in the form of a quadrilateral. The use of reduced integration on elements with straight edges (i.e., orthogonal isoparametric coordinates) gives rise to the Barlow point property (Barlow, 1989). Because of this, the strains calculated from the interpolation functions display higher accuracy at integration points than anywhere else in the element. This property is especially advantageous for nonlinear constitutive equations of the type considered here. Bassani and McClintock (1981) successfully applied the reduced integration technique for the power-law (incompressible) creep behavior of metals and determined the relaxation of crack-tip stresses. This paper adopts a reduced (2 × 2) Gaussian integration scheme with eight-noded isoparametric elements to model the incompressible transient and viscous creep response of ice. The mesh chosen for the present study performs very well and does not produce spurious modes of deformation which may be encountered when using reduced integration. 5.4. Crack-tip elements
At the crack-tip, the eight-noded quadrilateral elements are degenerated down to triangles by collapsing the three nodes of one side to occupy the same location. The crack-tip elements in elastic problems are often obtained by moving the mid-side nodes to the quarter points and tying the collapsed nodes;, which results in a r -1/= singularity (Barsoum, 1977). Within the creep zone around the crack-tip, the r -N/CN+I) singularity of the R R strain field becomes dominant (Riedel and Rice, 1980). Consequently, elastic singular elements are not appropriate for an elastic-creep analysis. A crack-
275
S. Nanthikesan, S. Shyam Sunder/Mechanics of Materials 21 (1995) 265-279
tip element that is capable of characterizing this higher singularity efficiently is one that is degenerated to a triangle as before, but the crack-tip nodes are untied and the middle nodes on the other three sides remain at mid-point. This elem e n t geometry with mid-point nodes produces a r-~I singularity on the displacement derivatives (Barsoum, 1977; Bassani and McClintock, 1981) which proves to be a good approximation for the weaker singularities of the R R and L E F M fields.
of the elastic and plastic strains, respectively, and is written as:
i1
E = ~- + c~eo o-
(37)
L°0 J
where ~ is a material constant, e 0 is the yield strain, and ~r0 is the yield stress. The asymptotic stress field under small-scale yielding conditions for this material can be written, in polar coordinates (r, 0) with origin at the crack-tip, as:
(J)N+I~ij(ON), 1
6. Validation of numerical model
O'ij = Or0
The ability of the numerical model to accurately predict crack-tip processes is evaluated by comparing the results of numerical simulations with known analytical solutions for the asymptotic fields and the creep zone under a plane strain m o d e I loading configuration. The results should be, and indeed are, independent of the specific choice of the radius of the boundary layer R, when the boundary is chosen sufficiently far from the creep zone. The first numerical example considers a rateindependent plastic material which exhibits classical power-law hardening of the R a m b e r g - O s g o o d type. The uniaxial form of the constitutive law for this material expresses the total strain as the sum
(38)
60.OCOINr
where J denotes the J-integral, which under plane strain SSY conditions is related to the stress intensity factor K I through J = K~(1 - v)e/E; I u is an integration constant (IN= 5.51 for plane strain conditions and N = 3) and 6is(O, N) is an angular function, both of which can be determined from Sih (1983). The above equation represents the well-known H R R field (Hutchinson, 1968; Rice and Rosengren, 1968). This material model, rather than a power-law creeping solid, is used to verify the effectiveness of the r -1 singular tip elements and reduced integration. While the asymptotic fields are selfsimilar (identical) for the two material models, the size of the inelastic process zone is indepen-
(e)
(A) to6
........
I
........
I
......
"1
.......
"T
.......
lOP
........
I
........
I
.......
I
.......
]
.......
J
.......
I
.......
0- 7"
o ,.90 ° I N-3
° l°s
.~ ]
10-1
!
l "~w°l°s
tF'~'''''" to~
1 ~ ~ ~ " ~
~. ~,.
tO"3
to3
io2
........
10"$
!
10.4
........
!
10-3
.......
J
10-2
. . . . .
,d
I0-I
i
,r
.,.,,I
10.
I00
i
5
10"i
i I o''''!
10.4
........
10.3
10.2
J
. . . . .
10-1
i-/a
Fig. 3. Numerical verification of asymptotic fields for elastic power-law hardening plastic material.
1o°
S. Nanthikesan, S. Shyam Sunder/Mechanics of Materials 21 (1995) 265-279
276
io3
' ' ...... I
....... I
........
I
.......
~
.......
0 = 45 ° IO2
I01
"-----2>
I0 o
i 0 "1
...---NUMERICAL = . , ..... I
i0 "2
10-5
i0 --4
. ,,,,.,i
= i rI~ITT I
10-3
1o.2
* ....... I
10-1
.......
100
r/a
normal opening and shear stresses along a 45 ° radial path emanating from the crack-tip is shown in Fig. 4. This is the instantaneous response (at t = 0 ÷) under a suddenly applied stress intensity of 100 KPa m 1/2 that is held constant with time. As expected, the numerical predictions are in excellent agreement with the well-known L E F M field. Further, the numerical estimate of the J-integral (9.43 × 10 .4 J / m 2) is within about 0.6 percent of the theoretical value (9.38× 1 0 . 4 J/m2). Under sustained loading, a creep zone develops around the crack-tip and grows with time.
Fig. 4. N u m e r i c a l v e r i f i c a t i o n of a s y m p t o t i c fields for elastic power-law c r e e p m a t e r i a l at t = 0 +.
(an)
POWER LAW
:o.oz
.KI = lO0 K P L m I ~
dent of time for the plastic material and this is convenient for purposes of evaluation. The material parameters used in the numerical simulation are: E = 9.5 GPa, v = 0.33, d = 1.6 × 10 .6 , ~r0 = 1 MPa, and N = 3. Note that e 0 = ~ro/E. Numerical predictions of the asymptotic fields for equivalent stress and equivalent strain, corresponding to a suddenly applied stress intensity of 100 kPa m 1/2 that is held constant, are shown in Fig. 3 along a 90 ° radial path emanating from the crack tip. For N = 3, the stress exhibits a r - 1 / 4 singularity while the strain exhibits a r - 3 / 4 singularity. The figure shows that the numerical results are in excellent agreement with the theoretical H R R fields. Only at large radial distances approaching the crack length does the stress field begin to follow the prescribed L E F M field. But this corresponds to large-scale plastic deformations in the test material, which are not of direct interest in this study. The second numerical example considers an isotropic material obeying an elastic power-law creep model of deformation. The elastic constants and power-law exponent is as before with the flow stress for freshwater polycrystalline ice at T = - 1 0 ° C obtained from Shyam Sunder and Wu (1989a) as V = 179.8 MPa. The asymptotic response of a typical pair of
i
~
t=lOs
(m)
0.01
-o:oz Traction-Free Crack
(a)
i0"1
pow'[~ LAW CREEP K.I = 100 K P & m 1/2 10.2
T =-10°C
N
M 10.3
10 "1
10 0
101
TIME (s)
(b) Fig. 5. (a) T i m e evolution of c r e e p zone boundary. (b) T i m e versus c r e e p zone size for a e l a s t i c - p o w e r law c r e e p m a t e r i a l u n d e r step loading.
S. Nanthikesan, S. Shyam Sunder~Mechanics of Materials 21 (1995) 265-279
Riedel and Rice (1980) express the boundary of this creep zone rcr(0, t) as follows: 2
~ 1 K 2 [ / ( N + 1) 2 Egot rcr = 2~[ 2 u 8 U + l vN
U--I
Fcr(O),
(39)
where For(O) is an angular function which depends on N and the loading condition, i.e., plane strain versus plane stress. This function, which is independent of time, defines the shape of the creep zone. Figure 5a shows the shape of the creep zone predicted by the numerical model at three time instants: t = 1, 3.3 and 10 seconds. These are in excellent agreement with the functional variation of Fcr(0) plotted in Fig. 1 of Riedel and Rice (1980). The numerical model predicts the size of the creep zone at a time equal to 10 seconds, which is typical of fracture tests and impact problems in ice, to be equal to 1.84 × 10 .2 mm. This corresponds to the maximum value of rc~ which occurs at 0 = 90 °. The Riedel and Rice equation given above yields a size of 1.73 × 10 .2 mm that is within 6.5 percent of the numerical prediction. In evaluating the above equation, it is assumed that a N~ 1 as shown by Riedel and Rice and the maximum value of F¢r = 0.25, which occurs at 0 = 90 ° for N = 3. The evolution of the size of the creep zone with time is plotted in Fig. 5b. The predicted linear variation is in agreement with the Riedel and Rice equation for N = 3.
7. Conclusions A numerical model has been developed to study the interaction between creep deformations and a stationary crack in polycrystalline ice. The problem is analyzed using the physically-based constitutive model of Shyam Sunder and Wu (1989a,b) in conjunction with the finite element method. The solution to the c r e e p - c r a c k interaction problem using this constitutive theory involves the numerical time integration of a coupled set of very stiff and highly nonlinear first-order differential equations. A closed-form N e w t o n - R a p h -
277
son (tangent) formulation in conjunction with an implicit form of the a-method of time integration is developed to solve the local constitutive equations. A consistent constitutive Jacobian matrix, which is used to assemble the global tangent stiffness matrix in a finite element analysis, is also established in closed form. The algorithm is implemented as a user-specified subroutine in the finite element program ABAQUS. The creep-crack interaction problem involving a traction-free semi-infinite crack embedded in a semi-infinite medium is simulated using eightnoded isoparametric elements which are arranged in a cascading, radially focussed mesh. The crack-tip element, formed by collapsing the three nodes on one side of an element, effectively models the singular stress and strain fields ahead of the crack-tip. A reduced Gaussian integration scheme is used to avoid numerical "locking" and increase computational efficiency. The ability of the numerical model to accurately predict crack-tip processes has been verified by comparing numerical predictions of: (a) the asymptotic fields with the well-known L E F M and H R R fields; and (b) the size, shape and time evolution of the creep zone produced by a suddenly applied load in an elastic power-law creep material with the analytical solutions of Riedel and Rice (1980). The companion paper (Nanthikesan and Shyam Sunder, 1995) presents the results from a comprehensive set of numerical simulations that analyze the role of transient creep, including loading history, temperature sensitivity, material strain hardening and fabric anisotropy, in predicting the size, shape and the time evolution of the creep zone in polycrystalline ice.
Acknowledgements The authors would like to thank Professor David M. Parks of MIT for many helpful discussions. The authors acknowledge financial support from the Office of Naval Research (Grant No. N00014-92-J-1208), the Minerals Management Service, and from an industry consortium consist-
278
S. Nanthikesan, S. Shyam Sunder/Mechanics of Materials 21 (1995) 265-279
ing of AMOCO, ARCO, BP America, CHEVRON, CONOCO, EXXON, MOBIL and TEXACO through the MIT Center for Scientific Excellence in Offshore Engineering.
References ABAQUS User's Manual, Version 4.7 (1988), Hibbit, Karlsson and Sorensen Inc., Providence, Rhode Island. Argyris, J.H., L.E. Vaz and K.J. William (1978), Improved solution methods for inelastic rate problems, Comp. Meth. App. Mech. Eng. 16, 231. Anderson, T.L. (1991), Computational Fracture Mechanics, CRC Press, p. 659. Barlow, J. (1989), More on optimal stress points - reduced integration, element distortions and error estimation, Int. J. Num. Meth. Eng. 28, 1487. Barsoum, R.S. (1977), Triangular quarter-point elements as elastic and perfectly-plastic crack tip elements, Int. J. Num. Meth. Eng. I1, 85. Bassani, J.L. and F.A. McClintock (1981), Creep relaxation of stress around a crack tip, Int..L Solids Struc. 17, 479. Bathe, K.J. (1982), Finite Element Procedures in Engineering Analysis, Prentice Hall, New Jersey. DeLorenzi, H.G. and C.F. Shih (1977), On the use of 2D isoparametric elements for calculations in the fully plastic range; Int. J. Fract. 13, 507. Dodds, R.H. (1982), Effects of reduced integration on the 2-D quadratic isoparametric element in plain strain plasticity. Int. J. Fract. Mech. 9, R75. Glen, J.W. (1955), The creep of polycrystalline ice, Proc. Roy. Soc. London A, 228 (1175), 519. Goodman, J.D. (1980), The Fracture Toughness of Ice, Data Report., Hokkaido University, Dept. of Applied Physics, Sapporo, Japan. Hinton, E., E.M. Salonen and N. Bicanic (1978), The Mathematics of Finite Elements and Applications, III, p. 437. Hutchinson, J.W. (1968), Plastic stress and strain fields at a crack tip, J. Mech. Phys. Solids 16 (13), 337. Lush, A.M., G. Weber and L. Anand (1989), An implicit time-integration procedure for a set of internal variable constitutive equations for isotropic elasto-viscoplasticity, J. Plasticity 5, 521. Nagtegal, J.C., D.M. Parks and J.R. Rice (1974), On numerically accurate finite element solutions in the fully plastic range, Comp. Meth. App. Mech. Eng. 4, 153. Nanthikesan, S: and S. Shyam Sunder (1995), Tensile cracks
in polycrystalline ice under transient creep: Part II Numerical Simulations, Mech. Mater 21,281. Nixon, W.A. and E.M. Schulson (1987), A micromechanical view of the fracture toughness of ice, J. Phys. Colloq. (France) 48, C1-313. Palmer, A.C., D.J. Goodman, M.F. Ashby, A.G. Evans, J.W. Hutchinson and A.R.S. Ponter (1983), Fracture and its role in determining ice forces on offshore structures, Annals Glac. 4, 216. Rice, J.A. (1968), A path independent integral and the approximate analysis of strain concentration on notches and cracks, J. App. Mech. Trans. ASME 6, 379. Rice, J.A., and G.F. Rosengren (1968), Plane strain deformation near a crack tip in a power-law hardening material, J. Mech. Phys. Solids 16, 1. Riedel, H. (1981), Creep deformations at crack tips in elasticviscoplastic solids, J. Mech. Phy. Solids 29, 35. Riedel, H. and J.R. Rice (1980), Tensile cracks in creeping solids, Fracture Mechanics: 12th conference, A S T M STP 700, p. 112. Schulson, E.M. (1990), The brittle compressive fracture of ice, Acta Met Mat. 38 (10), 1963. Schulson, E.M., S.G. Hoxie and W.A. Nixon (1989), The tensile strength of cracked ice, Philos. Mug. 59 (2), 303. Shyam Sunder, S. and M.S. Wu (1989a), A differential flow model for potycrystalline ice, Cold Reg. Sci. Technol. 16 (1), 45. Shyam Sunder, S. and M.S. Wu (1989b), A multiaxial differential mode[ for flow in orthotropic polycrystalline ice, Cold Reg. Sci. Technol. 16 (2), 223. Shyam Sunder, S. and M.S. Wu (1990), On the constitutive modeling of transient creep in polycrystalline ice. CoM Reg. Sci. Technol. 18, 267. Shyam Sunder, S., A. Elvin and S. Nanthikesan (1993), Numerical modeling of transient creep in polycrystalline ice, .L Engr. Mech. 119 (10) 2011. Sih, C.F. (1983), Small scale yielding analysis of mixed mode plane strain crack problem, ASTM STP 560, 187. Smith, D.J. (1990), Crack growth in brittle materials exhibiting primary creep. Eng. Fract. Mech. 37 (1), 27. Timco, G.W. and R.M.W. Frederking (1986), The effects of anisotropy and microcracks on the fracture toughness of freshwater ice, Proc. Fifth Int. Sym. Offshore Mechanics and Arctic Engineering, Tokyo, Japan, Vol. IV, p. 341. Urabe, N., and A. Yoshitake (1981), Strain-rate dependent fracture toughness of pure and sea ice, Proc. 6th Int. lAHR Sym. Ice, Quebec City, Canada, July 27-31, p. 551. William, K.J. (1978), Numerical solution of inelastic rate process, Comp. Struct. 8, 511. Wu, M.S. and S. Shyam Sunder (1992), On the constitutive modeling of transient creep in polycrystalline ice: Closure to Discussion. Cold Reg. Sci. Tech. 20 (3), 315.
S. Nanthikesan, S. Shyam Sunder~Mechanics of Materials 21 (1995) 265-279
Appendix G Matrix: a1 3
( a 1 + a3) 3
a3 3
(al +a2) 3 G=
a2 3 ( a 2 + a3) 3 2a4 Sym.
2a5 2a 6
H Matrix: 3 ( a I + a3) a. 2 a2
3ala2a3 a, 2 3(a1+a2) a*2
3ala2a 3 a*2 3alaea 3 a*2 3(a2 + a3) a*2
H=
a2 3
1 2a4 Sym.
1 2a5
1 2a 6 where a * = (a2a 3 + a3a I + a~a2 )2.
279