- Email: [email protected]
A computational probabilistic model for creep-damaging solids
> Pergamon
0045-7949(94)E0253-X
C,m,,wte,s & Srrurrurer Vol. 54. No. I. pp. 161-166. 1995 Copyright J 1994 El&r Science Ltd Prmted in Great Bmain. All rights reserved
A COMPUTATIONAL PROBABILISTIC MODEL FOR CREEP-DAMAGING SOLIDS D. G. Harlow and T. J. Delph Dept.
of Mechanical
Engineering
and Mechanics, (Received
Lehigh I7 May
University,
Bethlehem,
PA 18015, U.S.A.
1993)
Abstract-A
computational probabilistic model for creep-damaging solids is presented which takes into account the effects of randomness in creep and creep-rupture properties observed in experimental data. The model is based upon concepts from continuum damage mechanics, using the well-known Kachanov creep damage model as a prototype, and assuming that various material parameters present in the problem may be considered as random variables. The mapping from the random parameter space to the solution space is computed using exact relations from probability theory. As an example, the distributions for the parameters in the model are estimated from experimental data and the probabilistic growth of creep damage in a thick-walled cylinder under internal pressure is computed.
In the present work, we wish to develop a computational model of probabilistic creep damage which is based upon the principles of continuum damage mechanics and which takes into account both the scatter in creep deformation and rupture data. Although the method which we will present here is sufficiently general to be applied to more complicated creep and creep-damage constitutive models, we will choose, for the sake of illustration, to consider a somewhat simplified version of the well-known Kachanov constitutive model for a creep-damaging solid [8], one for which the Monkman-Grant relation is satisfied. This model involves a scalar damage variable n which is initially zero for an undamaged material and which, under a uniaxial stress CT,grows according to the rate equation
INTRODUCTION
It is by now well-established experimentally that the creep deformation and creep rupture behavior of metals exhibits marked probabilistic behavior, e.g. [l-3]. The scatter observed in creep deformation and failure data is usually regarded as being due to material variability, because it is known that modest variations in chemical composition, heat treatment, etc. can lead to substantial variations in mechanical properties. This scatter is of considerable technological importance, because it greatly complicates the task of making accurate deformation and lifetime estimates for high-temperature components. Accordingly, several efforts, e.g. [4-71, have been made to take account of random variations in materials properties in creep analyses. In the most recent of these, Delph and Yukich [7] considered the probabilistic creep of an elastic/powerlaw creeping body, for which, under a uniaxial stress 6, the strain rate is given by
i =;+Ca”.
f&z
(1 -cl)“’
where A is a material constant. Failure occurs at a value of R = 1. To take account of the effects of damage upon the deformation behavior, eqn (1) is modified to read
(1)
+c-
Here, t is the axial strain, E is Young’s modulus and the superposed dot indicates time rate. The constants C and N in eqn (I) were taken to be correlated random variables, with a joint probability density function (jpdf) derived from experimental data. Subject to the restriction that C and N be constant throughout the body and vary randomly only from realization to realization, a finite element-based method was derived which yielded the probability density function (pdf) for the mechanical variables of interest at a given instant of time.
( >. CT
N
1-R
Damage equations quite similar to the Kachanov model have been previously used to construct other statistical and stochastic damage models [9, IO]. In the present case, we take the point of view that the material constants C, N and A appearing in eqns (2) and (3) are correlated random variables whose distributions may be determined from experimental data. Our goal is to obtain the marginal pdfs for the 161
162
D. G. Harlow
mechanical variables of interest under arbitrary, but deterministic, time-dependent loading, at a given instant in time. To do this, we introduce the simplifying assumption that C, N and A are constant throughout the body, but may vary randomly from realization to realization of the system. Thus, we explicitly exclude any spatial stochastic variation. In what follows, we outline a numerical scheme for obtaining the marginal pdf for any of the mechanical variables of interest, based upon a finite element formulation for a solid whose constitutive behavior is given by the multiaxial generalization of eqns (2) and (3) and for which the parameters C, N and A are random variables having a known joint probability density function (jpdf). We then present a numerical example for the probabilistic evolution of damage in a creeping cylinder under internal pressure, in which the jpdf for C, N and A is determined from experimental data.
NUMERICAL
FRAMEWORK
Under a multiaxial stress state, may be generalized to
eqns (2) and (3)
and T. J. Delph
In order to solve these equations numerically, the constitutive equations (4) and (5) are integrated forward in time from the initial conditions, which we will take to be the elastic solution. The nodal point displacements are then obtained from equations (6-8) and the process is repeated iteratively up to the time of interest. Our goal is the following: given the jpdf for the random variables C, N and A, FC,,v,..A(~,n, a). we wish to calculate the marginal pdf for any of the mechanical variables of interest in the problem (here we make use of a standard notation for jpdfs, in which upper-case variables denote the random variables of interest, while lower-case variables denote the corresponding arguments in the jpdf). The mapping method by which we intend to accomplish this is outlined in [7] and is described in greater generality by Harlow and Delph [ Ill. We will briefly outline the method here, in the context of the present problem. Suppose that P is the variable of interest, where P might represent a nodal point displacement, the value of damage or stress at any point in the body, or any other of the variables in the problem. We represent the solution of the governing equations for the quantity P in the form. P = h(t, C, N, A).
(5) where s,, is the deviatoric strain tensor, c? is the von Mises effective stress, E is Young’s modulus and v is Poisson’s ratio. A standard finite element formulation for the deformation of a creeping body whose constitutive behavior is given by eqns (4) and (5) is
K6 =
BrCt’ dV + s I’
HJt dS.
(6)
sS
We assume that h has continuous partial derivatives with respect to each of its arguments and that it is one-to-one with at least one of the three parameters C, N and A. For simplicity, assume that the solution is one-to-one with respect to C. Then we may invert the solution in the form C = g(t, P, N, A). We now augment identification
the probability
N = Y;
Here K is the elastic stiffness matrix, 6 is the vector of nodal point displacements, B is the straindisplacement matrix, C is the matrix of elastic constants, H is the vector of the finite element shape functions, t is the traction vector, and L’ is the creep strain vector. Also, V is the volume of the body under consideration and S, is the portion of the bounding surface upon which tractions are prescribed. In addition to eqn (6), we have the finite element strain-displacement relation
(9)
(10)
space by making the
A =Z.
We wish to make a connection
(11) between
the known
jpdff,.,V,,d(c, n, a) and the jpdff,,,,(t,p, L’,2). To do this, we refer to a well-known theorem in probability theory governing transformations of random variables (‘the mapping theorem’), which states, e.g. [12], that
n =L’,a ==).
t = BS and Hooke’s
Law, which determines CT =
C(t - t”),
where t is the total strain vector vector.
(12)
(7) the creep strain:
Here J is the Jacobean mapping (P, Y, Z)+(C,
(8)
determinant N, A), i.e.
of the inverse
(13)
and (r is the stress The desired marginal
pdf for P at a given time t may
Model
then be obtained
by integration
solids
for creep-damaging
of eqn (I 2):
163
By virtue of eqn (II), all rows of the determinant below the first will have ones on the main diagonal and zeros elsewhere. Thus
stress levels. The data from these tests evidenced a substantial amount of scatter, both in the creep response and in the creep rupture failure times. We wish to use these data to estimate the jpdf for the parameters C, N, A in eqns (2) and (3). In particular, we will focus upon three data sets resulting from constant-load creep tests carried out at 198.9, 218.1 and 239.1 MPa, respectively, at a test temperature of 593°C. Each data set consisted of six replicated creep-rupture tests, with the mean failure times for the three sets being 1280, 325 and 77 hr, respectively. For each tests, the failure time t,, the strain at failure t,, the steady-state (secondary) strain rate i,, and the strain and the time to the end of the secondary creep phase, t2 and t,, were reported. It is these data that we will use to estimate the jpdffc,,,,(c, n, n). First of all, as long as the damage level is relatively low, the steady-state creep rate will be given by the familar power-law creep equation:
1
i,, = CuN.
.&(LP) =
fp.u.z(t,~>~> 2) dy dz.
(14)
ii The inversion required to calculate C = g(t, P, N, A) can be carried out by standard numerical techniques, which we shall outline subsequently. To calculate the Jacobean determinant J, we consider its inverse, K = l/J, where K is the Jacobean determinant of the direct mapping (C, N, A)+(P, X, Y) given by
c”P (16)
K=5=c.
The quantity aP/aC may be determined as the solution of a separate set of governing equations [7, 111 or by the use of finite difference techniques. In the present work, we shall make use of the latter method. In any event, once having determined the Jacobean determinant J and the inverse solution C = g(t, P, N, A), we may obtain the marginal pdf f,(t,p) from eqn (14) by numerical integration; this completes the basic solution procedure. We now pass to a detailed numerical example.
Moreover, under constant stress conditions, the uniaxial eqns (2) and (3) for the Kachanov constitutive model may be integrated to obtain 1
For the purposes of illustration, we will make use of the experimental results of Garofalo et al. [l]. These authors carried out a number of replicated creep and creep rupture tests on a single heat of AISI type 316 stainless steel at various temperatures and
(18)
“= (N + l)AoN and ~+(l-)I~i*+‘~].
Equation
EXAMPLE
(17)
(19)
(19) gives the creep strain
c,=-.
A
For each of the 18 reported values of (t,, t,), C,,, and (Ed, t2), a least-squares fit of eqns (17-20) was ”
90
:7--i
60
0
-I
0
0 0
70
0
60
0”
30 -
a
0” 0
20 0
i
50 0
20
1
IO
1
5c -49
-48
-47
-46
-45
0 0 0 6.9
6.6
7.0
WC) mN
(1)
Figs
:
0
1
2-
0
0
30
0
5-
0 0
40
0
10
as
C
95
a 40
at failure
I and 2. Ln(C) and N on probability
paper.
7.1
D. G. Harlow
164 95
and T. J. Delph
In order to convert eqn (21) to standard form, we introduce the change of variables:
90
normal
z=M/tx,
60
(24)
0 70
0
a
where M is the matrix having the eigenvectors of V as its columns and n is the diagonal matrix having the quantities JE., , J& and ,/A, on its diagonal, where i., , i2 and i, are the eigenvalues of V. We then have
0
60
0
0 0
50
0
40
0 0
30
0
I
0
20
.L @I?x2 ?-x31= (271 __
0 10
)3;2
exp( -f(xT
0
+ if + xi)). (25)
5 -46
-44
-45 WQ
Fig. 3. Ln(A) on probability paper. conducted, in order to estimate values of the parameters C, N and A for each test. For the sake of consistency with the power-law creep model, the primary creep strain was substracted from the reported values of t, and t2 before carrying out the fits. The results are shown in Figs 1-3, plotted on log-normal, normal and log-normal probability paper, respectively. It can be seen from an inspection of Figs l-3 that C and N can be reasonably well-described by lognormal and normal distributions, respectively. The distribution for A is slightly more problematic; however, the assumption of a log-normal distribution for this quantity does not seem to be too unreasonable. Thus we will assume that the quantities C, N and A have a jointly log-normal, normal, log-normal distribution. The parameters for the jpdf can be obtained from the statistics of the data, which yield p,,c = -46.9, qlnc = 0.89, pn = 6.91, qK = 0.061, pInA= -45.0, R”,~= 1.00, rl,c,l,A=0.923, rlnC.,v= - 0.843, TN,,” A = -0.824, with C and A in units of MPa-N hr-‘. Here p denotes the expected value, rl the standard deviation, and Y the correlation coefficient. The jpdf for C, N and A then has the standard form [12] I
f inC,N.lnA(ln c, n, ln a) =
exp( - fzY--‘z),
We now consider the creep and creep-damage behavior of a cylinder under plane strain conditions, having an inner radius of a = 50 mm, an outer radius of b = 1OOmm and subjected to a constant internal pressure of p = 60 MPa. The constants C, N and A in the creep and creep-damage constitutive laws are distributed as described above. The elastic constants are assumed to be deterministic and are taken to be E = 1.49 x IO5 MPa and v = 0.3. A finite element formulation for this problem, similar to that given in takes advantage of the oneeqn (6) which dimensional nature of the problem, was used to solve the field equations. For the sake of simplicity, we focus attention upon the value of the damage variable at the inner radius, Q(t, a) = W. We assume [e.g. eqn (9)] that W is determinable in terms of the three parameters X, , X,, A’, in the form
w = h(t, x,,
x2, X3).
(26)
Likewise, we assume that h has continuous partial derivatives with respect to each of its arguments and that it is one-to-one with respect to at least X, With the identification Y=X2;
z=x,.
(27)
we may invert eqn (26) in the form X,=g(W,
Y,Zt).
(28)
(2743:qq (21)
where (22)
z=
The connection between the known jpdf ,f,(_~,, .x2, x3) and the jpdf fw,v.z(w,y, z) is given by eqn (12) i.e.
fW.Y,Z(~~~~ Y>J3t) = IJK,X,.X~.X,(XI = g(M.3 I’>22t)7 x2
and V is the variance-covariance
matrix
given by
=
Y, Y3 = z)
(29)
Model for creep-damaging solids and the Jacobean
99 -
J is given by
96 -
aw -1
(>
J= ax, The marginal then obtained
- -
95
90.
pdf for the quantity from eqn (29) by
of interest
W is
Computed,
500 hrs.
Computed,
1000 hn.
A
Monte Carlo, 500 hrs.
0
Monte Carlo, loo0
hrs.
.
80
70 2
if;,.(w, t) = T. J. fw.r.z(w, ~2 z, t) dy dz. ss--1. --71
(31)
60 50 40 30.
In the numerical evaluation of eqn (31), the inverse function X, = g( W, Y, 2, t) was evaluated, for fixed values of the arguments, by a numerical shooting scheme in which X, was systematically varied in the solution of the finite element equations until the solution led to the desired value of W. The Jacobean J was then calculated by using a finite difference perturbation about this point. In order to evaluate the double integral in (31) for a fixed value of t and w, a numerical optimization routine was used to determine, for fixed values of Y, the value of 2 corresponding to the maximum value of the integrand in eqn (31). Gaussian integration about this value of Z was then used to carry out the integration with respect to Z, using truncated limits. This was done over a range of Y values, so that, at the limits of the range, the integrand decreased to a sufficiently small value as to have a negligible contribution to the integral. A cubic spline was then fitted to the result and integrated with respect to Y to yield the value off&w, t). This procedure was repeated for different values of W to yield a series of pointwise values for f;,.(w, t). Figure 4 shows the results of this procedure for t = 500, 750 and 1000 hr. As might be expected, the pdf shows a marked tendency to broaden with increasing time, indicating increasing scatter in W with increasing time. We might note that the damage values shown are too
300
...... 250
-
500 hrs.
- 750 hrs. 1OOOhrs.
200
2
150
100
50
0 ”
100
50 WXlO’
Fig. 4. Calculated pdfs for inner radius damage.
150
20 -
d 0
10 -
_1 2.0
2.5
3.0
3.5
4.0
4.5
5.0
In(WX1o-q
Fig. 5. Calculated cdfs for inner radius damage and Monte Carlo results.
small to have an appreciable effect upon the creep deformation behavior of the body. For the parameter distributions used in the present problem, the effect of damage upon the creep deformation becomes evident in the regime in which the strains become too large to invalidate the use of the linearized kinematics used in the finite element code employed in the calculations. The corresponding cdfs for r = 500 and 1000 hr, obtained by integrating the pdfs, are shown in Fig. 5. Also shown in Fig. 6 are the results of 200 Monte Carlo trials for each value of time. As may be seen, the agreement between the empirical and the calculated cdfs is excellent. DISCUSSION
We have presented here a computational probabilistic model for a creep-damaging material which takes into account both the randomness present in the creep deformation and the creep rupture data. We have formulated the description in terms of the simplest of the creep-damage models, the well-known Kachanov model, but the extension to more complex creep-damage models is immediate. Our formulation assumes that the parameters in the model are constant for each realization, but that the parameters vary in a random fashion from realization to realization, with a known jpdf. We have shown how the distributions for these parameters may be estimated from experimental data. This formulation would, for example, describe the situation in which the variations in material result from parameters random heat-to-heat variations or from variations in factors such as mechanical processing or heat treatment. Specifically excluded from the formulation is the situation in which the parameters vary in a spatially stochastic
D. G. Harlow
166
manner within the body. However, as long as the dimensions of the body are comparatively small relative to the characteristic length scale for stochastic variation, or if the stochastic variation is small compared to the variability between different realizations of the system, then the assumption of spatially-constant parameters is likely to provide good description of the random creep/creep-damage response of the system. The formulation we have proposed has several significant advantages. The first of these is that it may make use of unmodified finite element codes, if, as was done here, a finite difference estimate for a Jacobean J is accepted. The second is that it makes use of exact relations in the theory of probability, with the only approximations introduced being those inherent in the numerical technique used to calculate the deformation, stress and damage fields, e.g. finite element techniques. As a result, it is quite accurate, as may be seen by a comparison of the calculated cdfs with the results of Monte Carlo simulations (Fig. 6). We may contrast this with firstand second-order reliability methods, e.g. [l3, 141, which also yield much the same information as the present technique. These methods, however, use several ad hoc approximations whose applicability to time-dependent, highly nonlinear problems such as we have considered here have yet to be demonstrated. The method presented here also has its disadvantages, however. The principal one is that it is computationally demanding, requiring that the finite element equations for the problem be solved repetitively a fairly large number of times. Moreover, as might be expected, the numerical effort rises steeply as the number of random variables in the problem increases. The evaluation of the marginal pdf from the double integral given by eqn (31) also poses some difficulties. In the present case. it was found useful to truncate the limits of integration. using optimization techniques to locate the maximum of the integrand and hence to define the appropriate limits of integration. We have chosen the problem of a creep-damaging cylinder under constant internal pressure to illustrate the method. The jpdf for the three material parameters C, N and A required by the model were estimated from a set of replicated experimental data for type 316 stainless steel. The results of this example show a clear tendency for the distribution of the damage variable at the inner radius of the cylinder to broaden with time, indicating an increasing degree of scatter with time. Unfortunately, the linearized kinematics used in the present work made it impossible to continue the calculations, the linearized kinematics used in the present work made it impossible to continue damage The
the calculations
parameters
use of linearized
essential
restriction.
in the
up to values
neighborhood
kinematics
is not,
of
of the unity.
however,
an
and T. J. Delph The scatter noted in replicated uniaxial creep and creep rupture data made it clear that the prediction of component creep failure is properly viewed as a probabilistic problem. One of the important features of the method presented here is that it provides a rational approach to probabilistic creep failure predictions. In the context of the Kachanov the probability of failure is just model, I - PrjR < I), a quantity which the method outlined here yields directly. For the example presented herein, such a calculation would have required the use of large-deformation kinematics, which we have not considered in the present work, but which would be a straightforward extension. This calculation represents a topic for future research. Ackno~,k~~d~ernc~,lrTJD support of the National INT-8922408.
would like to acknowledge the Science Foundation under grant no.
REFERENCES
I. F. Garofalo, R. W. Whitmore, W. F. Domis and F. von Gemmingen, Creep and creep-rupture relationships in an austenitic stainless steel. Trms. Me/all. Sot. AIME 221, 310-319 (1961). 2. K. F. A. Walles. Rundon? md S~~.s/erm/ic Fuc/or.s in /he SCU//CI. of Creep Data. C.P. No. 935, Aeronautical Research Council. H.M.S.O.. London (1967). 3. 1. Sprung and V. A. Zilberstein. In’ Under.~/onding Voribhili!,~ in Creep and Ruprure Behavior (Edited by M. Prarer and J. D. Parker), MPC Vol. 28. ASME, New York (1988). and W. N. Huang. Effect of 4. F. A. Cozzarelli random material parameters on nonlinear steady creep solutions. In/. J. Solids Struct. 7, 1477 (1971). 5. H. Broberg and R. Westlund, Creep in structures with random material properties. In/. J. Solids Struct. 14, 365 (1978). and R. Westlund, Creep rupture of 6. H. Broberg specimens with random material properties. I///. J. Solid,r Strut/. 14, 959-970 (1978). I. T. J. Delph and J. E. Yukich, A finite element method for the probabilistic creep of solids. Inr. J. Nun/ Me/l/. Eng. 35, 1171~~1182 (1992). Inrroduciion to Cor7tinuun1 Damqy 8. L. M. Kachanov. Mrchonics. Martinus Nijhoff. Dordrecht (I 986). and M. A. G. Silva. Statistical aspects 9. D. Krajcinovic of the continuous damage theory. I/t/. J. Solid.\ Srruc~r. 7, 551-562 (1982). IO. M. H. J. W. Paas, C. W. J. Oomens. P. J. G. Schreurs and J. D. Janssen, The mechanical behaviour of continuous media with stochastic damage. Engng Frrrcf. Mech. 36, 255-m266 (1990). II. D. G. Harlow and T. J. Delph, The numerical solution of random initial value problems. Mu/h. Compul. Sin!. 33, 243-258
12.
1. Elishakoff,
(1991). Prohahilkric
Me/hod.\
in the
Theory, of
John Wiley, New York (1983). 13. T. A. Cruse. 0. H. Burnside, Y.-T. Wu. E. Z. Polch and J. B. Dias. Probabilistic structural analysis methods for select space propulsion system structural components (PSAM). Cornp. S/ruct. 29, 891-901 (1988). and P.-L. Liu, First and second14. A. Der Kiureghian order finite element reliability methods, Computu/ional Srructurm.
Mechcmic.~
of
Probabilistic
cmd
Relic/hi/i/j,
(Edited by W.-K. Liu and T. Belytschko). Elmepress. Base1 (I 989).
Ancr1y.G.~
pp. 281-298.
0045-7949(94)E0253-X
C,m,,wte,s & Srrurrurer Vol. 54. No. I. pp. 161-166. 1995 Copyright J 1994 El&r Science Ltd Prmted in Great Bmain. All rights reserved
A COMPUTATIONAL PROBABILISTIC MODEL FOR CREEP-DAMAGING SOLIDS D. G. Harlow and T. J. Delph Dept.
of Mechanical
Engineering
and Mechanics, (Received
Lehigh I7 May
University,
Bethlehem,
PA 18015, U.S.A.
1993)
Abstract-A
computational probabilistic model for creep-damaging solids is presented which takes into account the effects of randomness in creep and creep-rupture properties observed in experimental data. The model is based upon concepts from continuum damage mechanics, using the well-known Kachanov creep damage model as a prototype, and assuming that various material parameters present in the problem may be considered as random variables. The mapping from the random parameter space to the solution space is computed using exact relations from probability theory. As an example, the distributions for the parameters in the model are estimated from experimental data and the probabilistic growth of creep damage in a thick-walled cylinder under internal pressure is computed.
In the present work, we wish to develop a computational model of probabilistic creep damage which is based upon the principles of continuum damage mechanics and which takes into account both the scatter in creep deformation and rupture data. Although the method which we will present here is sufficiently general to be applied to more complicated creep and creep-damage constitutive models, we will choose, for the sake of illustration, to consider a somewhat simplified version of the well-known Kachanov constitutive model for a creep-damaging solid [8], one for which the Monkman-Grant relation is satisfied. This model involves a scalar damage variable n which is initially zero for an undamaged material and which, under a uniaxial stress CT,grows according to the rate equation
INTRODUCTION
It is by now well-established experimentally that the creep deformation and creep rupture behavior of metals exhibits marked probabilistic behavior, e.g. [l-3]. The scatter observed in creep deformation and failure data is usually regarded as being due to material variability, because it is known that modest variations in chemical composition, heat treatment, etc. can lead to substantial variations in mechanical properties. This scatter is of considerable technological importance, because it greatly complicates the task of making accurate deformation and lifetime estimates for high-temperature components. Accordingly, several efforts, e.g. [4-71, have been made to take account of random variations in materials properties in creep analyses. In the most recent of these, Delph and Yukich [7] considered the probabilistic creep of an elastic/powerlaw creeping body, for which, under a uniaxial stress 6, the strain rate is given by
i =;+Ca”.
f&z
(1 -cl)“’
where A is a material constant. Failure occurs at a value of R = 1. To take account of the effects of damage upon the deformation behavior, eqn (1) is modified to read
(1)
+c-
Here, t is the axial strain, E is Young’s modulus and the superposed dot indicates time rate. The constants C and N in eqn (I) were taken to be correlated random variables, with a joint probability density function (jpdf) derived from experimental data. Subject to the restriction that C and N be constant throughout the body and vary randomly only from realization to realization, a finite element-based method was derived which yielded the probability density function (pdf) for the mechanical variables of interest at a given instant of time.
( >. CT
N
1-R
Damage equations quite similar to the Kachanov model have been previously used to construct other statistical and stochastic damage models [9, IO]. In the present case, we take the point of view that the material constants C, N and A appearing in eqns (2) and (3) are correlated random variables whose distributions may be determined from experimental data. Our goal is to obtain the marginal pdfs for the 161
162
D. G. Harlow
mechanical variables of interest under arbitrary, but deterministic, time-dependent loading, at a given instant in time. To do this, we introduce the simplifying assumption that C, N and A are constant throughout the body, but may vary randomly from realization to realization of the system. Thus, we explicitly exclude any spatial stochastic variation. In what follows, we outline a numerical scheme for obtaining the marginal pdf for any of the mechanical variables of interest, based upon a finite element formulation for a solid whose constitutive behavior is given by the multiaxial generalization of eqns (2) and (3) and for which the parameters C, N and A are random variables having a known joint probability density function (jpdf). We then present a numerical example for the probabilistic evolution of damage in a creeping cylinder under internal pressure, in which the jpdf for C, N and A is determined from experimental data.
NUMERICAL
FRAMEWORK
Under a multiaxial stress state, may be generalized to
eqns (2) and (3)
and T. J. Delph
In order to solve these equations numerically, the constitutive equations (4) and (5) are integrated forward in time from the initial conditions, which we will take to be the elastic solution. The nodal point displacements are then obtained from equations (6-8) and the process is repeated iteratively up to the time of interest. Our goal is the following: given the jpdf for the random variables C, N and A, FC,,v,..A(~,n, a). we wish to calculate the marginal pdf for any of the mechanical variables of interest in the problem (here we make use of a standard notation for jpdfs, in which upper-case variables denote the random variables of interest, while lower-case variables denote the corresponding arguments in the jpdf). The mapping method by which we intend to accomplish this is outlined in [7] and is described in greater generality by Harlow and Delph [ Ill. We will briefly outline the method here, in the context of the present problem. Suppose that P is the variable of interest, where P might represent a nodal point displacement, the value of damage or stress at any point in the body, or any other of the variables in the problem. We represent the solution of the governing equations for the quantity P in the form. P = h(t, C, N, A).
(5) where s,, is the deviatoric strain tensor, c? is the von Mises effective stress, E is Young’s modulus and v is Poisson’s ratio. A standard finite element formulation for the deformation of a creeping body whose constitutive behavior is given by eqns (4) and (5) is
K6 =
BrCt’ dV + s I’
HJt dS.
(6)
sS
We assume that h has continuous partial derivatives with respect to each of its arguments and that it is one-to-one with at least one of the three parameters C, N and A. For simplicity, assume that the solution is one-to-one with respect to C. Then we may invert the solution in the form C = g(t, P, N, A). We now augment identification
the probability
N = Y;
Here K is the elastic stiffness matrix, 6 is the vector of nodal point displacements, B is the straindisplacement matrix, C is the matrix of elastic constants, H is the vector of the finite element shape functions, t is the traction vector, and L’ is the creep strain vector. Also, V is the volume of the body under consideration and S, is the portion of the bounding surface upon which tractions are prescribed. In addition to eqn (6), we have the finite element strain-displacement relation
(9)
(10)
space by making the
A =Z.
We wish to make a connection
(11) between
the known
jpdff,.,V,,d(c, n, a) and the jpdff,,,,(t,p, L’,2). To do this, we refer to a well-known theorem in probability theory governing transformations of random variables (‘the mapping theorem’), which states, e.g. [12], that
n =L’,a ==).
t = BS and Hooke’s
Law, which determines CT =
C(t - t”),
where t is the total strain vector vector.
(12)
(7) the creep strain:
Here J is the Jacobean mapping (P, Y, Z)+(C,
(8)
determinant N, A), i.e.
of the inverse
(13)
and (r is the stress The desired marginal
pdf for P at a given time t may
Model
then be obtained
by integration
solids
for creep-damaging
of eqn (I 2):
163
By virtue of eqn (II), all rows of the determinant below the first will have ones on the main diagonal and zeros elsewhere. Thus
stress levels. The data from these tests evidenced a substantial amount of scatter, both in the creep response and in the creep rupture failure times. We wish to use these data to estimate the jpdf for the parameters C, N, A in eqns (2) and (3). In particular, we will focus upon three data sets resulting from constant-load creep tests carried out at 198.9, 218.1 and 239.1 MPa, respectively, at a test temperature of 593°C. Each data set consisted of six replicated creep-rupture tests, with the mean failure times for the three sets being 1280, 325 and 77 hr, respectively. For each tests, the failure time t,, the strain at failure t,, the steady-state (secondary) strain rate i,, and the strain and the time to the end of the secondary creep phase, t2 and t,, were reported. It is these data that we will use to estimate the jpdffc,,,,(c, n, n). First of all, as long as the damage level is relatively low, the steady-state creep rate will be given by the familar power-law creep equation:
1
i,, = CuN.
.&(LP) =
fp.u.z(t,~>~> 2) dy dz.
(14)
ii The inversion required to calculate C = g(t, P, N, A) can be carried out by standard numerical techniques, which we shall outline subsequently. To calculate the Jacobean determinant J, we consider its inverse, K = l/J, where K is the Jacobean determinant of the direct mapping (C, N, A)+(P, X, Y) given by
c”P (16)
K=5=c.
The quantity aP/aC may be determined as the solution of a separate set of governing equations [7, 111 or by the use of finite difference techniques. In the present work, we shall make use of the latter method. In any event, once having determined the Jacobean determinant J and the inverse solution C = g(t, P, N, A), we may obtain the marginal pdf f,(t,p) from eqn (14) by numerical integration; this completes the basic solution procedure. We now pass to a detailed numerical example.
Moreover, under constant stress conditions, the uniaxial eqns (2) and (3) for the Kachanov constitutive model may be integrated to obtain 1
For the purposes of illustration, we will make use of the experimental results of Garofalo et al. [l]. These authors carried out a number of replicated creep and creep rupture tests on a single heat of AISI type 316 stainless steel at various temperatures and
(18)
“= (N + l)AoN and ~+(l-)I~i*+‘~].
Equation
EXAMPLE
(17)
(19)
(19) gives the creep strain
c,=-.
A
For each of the 18 reported values of (t,, t,), C,,, and (Ed, t2), a least-squares fit of eqns (17-20) was ”
90
:7--i
60
0
-I
0
0 0
70
0
60
0”
30 -
a
0” 0
20 0
i
50 0
20
1
IO
1
5c -49
-48
-47
-46
-45
0 0 0 6.9
6.6
7.0
WC) mN
(1)
Figs
:
0
1
2-
0
0
30
0
5-
0 0
40
0
10
as
C
95
a 40
at failure
I and 2. Ln(C) and N on probability
paper.
7.1
D. G. Harlow
164 95
and T. J. Delph
In order to convert eqn (21) to standard form, we introduce the change of variables:
90
normal
z=M/tx,
60
(24)
0 70
0
a
where M is the matrix having the eigenvectors of V as its columns and n is the diagonal matrix having the quantities JE., , J& and ,/A, on its diagonal, where i., , i2 and i, are the eigenvalues of V. We then have
0
60
0
0 0
50
0
40
0 0
30
0
I
0
20
.L @I?x2 ?-x31= (271 __
0 10
)3;2
exp( -f(xT
0
+ if + xi)). (25)
5 -46
-44
-45 WQ
Fig. 3. Ln(A) on probability paper. conducted, in order to estimate values of the parameters C, N and A for each test. For the sake of consistency with the power-law creep model, the primary creep strain was substracted from the reported values of t, and t2 before carrying out the fits. The results are shown in Figs 1-3, plotted on log-normal, normal and log-normal probability paper, respectively. It can be seen from an inspection of Figs l-3 that C and N can be reasonably well-described by lognormal and normal distributions, respectively. The distribution for A is slightly more problematic; however, the assumption of a log-normal distribution for this quantity does not seem to be too unreasonable. Thus we will assume that the quantities C, N and A have a jointly log-normal, normal, log-normal distribution. The parameters for the jpdf can be obtained from the statistics of the data, which yield p,,c = -46.9, qlnc = 0.89, pn = 6.91, qK = 0.061, pInA= -45.0, R”,~= 1.00, rl,c,l,A=0.923, rlnC.,v= - 0.843, TN,,” A = -0.824, with C and A in units of MPa-N hr-‘. Here p denotes the expected value, rl the standard deviation, and Y the correlation coefficient. The jpdf for C, N and A then has the standard form [12] I
f inC,N.lnA(ln c, n, ln a) =
exp( - fzY--‘z),
We now consider the creep and creep-damage behavior of a cylinder under plane strain conditions, having an inner radius of a = 50 mm, an outer radius of b = 1OOmm and subjected to a constant internal pressure of p = 60 MPa. The constants C, N and A in the creep and creep-damage constitutive laws are distributed as described above. The elastic constants are assumed to be deterministic and are taken to be E = 1.49 x IO5 MPa and v = 0.3. A finite element formulation for this problem, similar to that given in takes advantage of the oneeqn (6) which dimensional nature of the problem, was used to solve the field equations. For the sake of simplicity, we focus attention upon the value of the damage variable at the inner radius, Q(t, a) = W. We assume [e.g. eqn (9)] that W is determinable in terms of the three parameters X, , X,, A’, in the form
w = h(t, x,,
x2, X3).
(26)
Likewise, we assume that h has continuous partial derivatives with respect to each of its arguments and that it is one-to-one with respect to at least X, With the identification Y=X2;
z=x,.
(27)
we may invert eqn (26) in the form X,=g(W,
Y,Zt).
(28)
(2743:qq (21)
where (22)
z=
The connection between the known jpdf ,f,(_~,, .x2, x3) and the jpdf fw,v.z(w,y, z) is given by eqn (12) i.e.
fW.Y,Z(~~~~ Y>J3t) = IJK,X,.X~.X,(XI = g(M.3 I’>22t)7 x2
and V is the variance-covariance
matrix
given by
=
Y, Y3 = z)
(29)
Model for creep-damaging solids and the Jacobean
99 -
J is given by
96 -
aw -1
(>
J= ax, The marginal then obtained
- -
95
90.
pdf for the quantity from eqn (29) by
of interest
W is
Computed,
500 hrs.
Computed,
1000 hn.
A
Monte Carlo, 500 hrs.
0
Monte Carlo, loo0
hrs.
.
80
70 2
if;,.(w, t) = T. J. fw.r.z(w, ~2 z, t) dy dz. ss--1. --71
(31)
60 50 40 30.
In the numerical evaluation of eqn (31), the inverse function X, = g( W, Y, 2, t) was evaluated, for fixed values of the arguments, by a numerical shooting scheme in which X, was systematically varied in the solution of the finite element equations until the solution led to the desired value of W. The Jacobean J was then calculated by using a finite difference perturbation about this point. In order to evaluate the double integral in (31) for a fixed value of t and w, a numerical optimization routine was used to determine, for fixed values of Y, the value of 2 corresponding to the maximum value of the integrand in eqn (31). Gaussian integration about this value of Z was then used to carry out the integration with respect to Z, using truncated limits. This was done over a range of Y values, so that, at the limits of the range, the integrand decreased to a sufficiently small value as to have a negligible contribution to the integral. A cubic spline was then fitted to the result and integrated with respect to Y to yield the value off&w, t). This procedure was repeated for different values of W to yield a series of pointwise values for f;,.(w, t). Figure 4 shows the results of this procedure for t = 500, 750 and 1000 hr. As might be expected, the pdf shows a marked tendency to broaden with increasing time, indicating increasing scatter in W with increasing time. We might note that the damage values shown are too
300
...... 250
-
500 hrs.
- 750 hrs. 1OOOhrs.
200
2
150
100
50
0 ”
100
50 WXlO’
Fig. 4. Calculated pdfs for inner radius damage.
150
20 -
d 0
10 -
_1 2.0
2.5
3.0
3.5
4.0
4.5
5.0
In(WX1o-q
Fig. 5. Calculated cdfs for inner radius damage and Monte Carlo results.
small to have an appreciable effect upon the creep deformation behavior of the body. For the parameter distributions used in the present problem, the effect of damage upon the creep deformation becomes evident in the regime in which the strains become too large to invalidate the use of the linearized kinematics used in the finite element code employed in the calculations. The corresponding cdfs for r = 500 and 1000 hr, obtained by integrating the pdfs, are shown in Fig. 5. Also shown in Fig. 6 are the results of 200 Monte Carlo trials for each value of time. As may be seen, the agreement between the empirical and the calculated cdfs is excellent. DISCUSSION
We have presented here a computational probabilistic model for a creep-damaging material which takes into account both the randomness present in the creep deformation and the creep rupture data. We have formulated the description in terms of the simplest of the creep-damage models, the well-known Kachanov model, but the extension to more complex creep-damage models is immediate. Our formulation assumes that the parameters in the model are constant for each realization, but that the parameters vary in a random fashion from realization to realization, with a known jpdf. We have shown how the distributions for these parameters may be estimated from experimental data. This formulation would, for example, describe the situation in which the variations in material result from parameters random heat-to-heat variations or from variations in factors such as mechanical processing or heat treatment. Specifically excluded from the formulation is the situation in which the parameters vary in a spatially stochastic
D. G. Harlow
166
manner within the body. However, as long as the dimensions of the body are comparatively small relative to the characteristic length scale for stochastic variation, or if the stochastic variation is small compared to the variability between different realizations of the system, then the assumption of spatially-constant parameters is likely to provide good description of the random creep/creep-damage response of the system. The formulation we have proposed has several significant advantages. The first of these is that it may make use of unmodified finite element codes, if, as was done here, a finite difference estimate for a Jacobean J is accepted. The second is that it makes use of exact relations in the theory of probability, with the only approximations introduced being those inherent in the numerical technique used to calculate the deformation, stress and damage fields, e.g. finite element techniques. As a result, it is quite accurate, as may be seen by a comparison of the calculated cdfs with the results of Monte Carlo simulations (Fig. 6). We may contrast this with firstand second-order reliability methods, e.g. [l3, 141, which also yield much the same information as the present technique. These methods, however, use several ad hoc approximations whose applicability to time-dependent, highly nonlinear problems such as we have considered here have yet to be demonstrated. The method presented here also has its disadvantages, however. The principal one is that it is computationally demanding, requiring that the finite element equations for the problem be solved repetitively a fairly large number of times. Moreover, as might be expected, the numerical effort rises steeply as the number of random variables in the problem increases. The evaluation of the marginal pdf from the double integral given by eqn (31) also poses some difficulties. In the present case. it was found useful to truncate the limits of integration. using optimization techniques to locate the maximum of the integrand and hence to define the appropriate limits of integration. We have chosen the problem of a creep-damaging cylinder under constant internal pressure to illustrate the method. The jpdf for the three material parameters C, N and A required by the model were estimated from a set of replicated experimental data for type 316 stainless steel. The results of this example show a clear tendency for the distribution of the damage variable at the inner radius of the cylinder to broaden with time, indicating an increasing degree of scatter with time. Unfortunately, the linearized kinematics used in the present work made it impossible to continue the calculations, the linearized kinematics used in the present work made it impossible to continue damage The
the calculations
parameters
use of linearized
essential
restriction.
in the
up to values
neighborhood
kinematics
is not,
of
of the unity.
however,
an
and T. J. Delph The scatter noted in replicated uniaxial creep and creep rupture data made it clear that the prediction of component creep failure is properly viewed as a probabilistic problem. One of the important features of the method presented here is that it provides a rational approach to probabilistic creep failure predictions. In the context of the Kachanov the probability of failure is just model, I - PrjR < I), a quantity which the method outlined here yields directly. For the example presented herein, such a calculation would have required the use of large-deformation kinematics, which we have not considered in the present work, but which would be a straightforward extension. This calculation represents a topic for future research. Ackno~,k~~d~ernc~,lrTJD support of the National INT-8922408.
would like to acknowledge the Science Foundation under grant no.
REFERENCES
I. F. Garofalo, R. W. Whitmore, W. F. Domis and F. von Gemmingen, Creep and creep-rupture relationships in an austenitic stainless steel. Trms. Me/all. Sot. AIME 221, 310-319 (1961). 2. K. F. A. Walles. Rundon? md S~~.s/erm/ic Fuc/or.s in /he SCU//CI. of Creep Data. C.P. No. 935, Aeronautical Research Council. H.M.S.O.. London (1967). 3. 1. Sprung and V. A. Zilberstein. In’ Under.~/onding Voribhili!,~ in Creep and Ruprure Behavior (Edited by M. Prarer and J. D. Parker), MPC Vol. 28. ASME, New York (1988). and W. N. Huang. Effect of 4. F. A. Cozzarelli random material parameters on nonlinear steady creep solutions. In/. J. Solids Struct. 7, 1477 (1971). 5. H. Broberg and R. Westlund, Creep in structures with random material properties. In/. J. Solids Struct. 14, 365 (1978). and R. Westlund, Creep rupture of 6. H. Broberg specimens with random material properties. I///. J. Solid,r Strut/. 14, 959-970 (1978). I. T. J. Delph and J. E. Yukich, A finite element method for the probabilistic creep of solids. Inr. J. Nun/ Me/l/. Eng. 35, 1171~~1182 (1992). Inrroduciion to Cor7tinuun1 Damqy 8. L. M. Kachanov. Mrchonics. Martinus Nijhoff. Dordrecht (I 986). and M. A. G. Silva. Statistical aspects 9. D. Krajcinovic of the continuous damage theory. I/t/. J. Solid.\ Srruc~r. 7, 551-562 (1982). IO. M. H. J. W. Paas, C. W. J. Oomens. P. J. G. Schreurs and J. D. Janssen, The mechanical behaviour of continuous media with stochastic damage. Engng Frrrcf. Mech. 36, 255-m266 (1990). II. D. G. Harlow and T. J. Delph, The numerical solution of random initial value problems. Mu/h. Compul. Sin!. 33, 243-258
12.
1. Elishakoff,
(1991). Prohahilkric
Me/hod.\
in the
Theory, of
John Wiley, New York (1983). 13. T. A. Cruse. 0. H. Burnside, Y.-T. Wu. E. Z. Polch and J. B. Dias. Probabilistic structural analysis methods for select space propulsion system structural components (PSAM). Cornp. S/ruct. 29, 891-901 (1988). and P.-L. Liu, First and second14. A. Der Kiureghian order finite element reliability methods, Computu/ional Srructurm.
Mechcmic.~
of
Probabilistic
cmd
Relic/hi/i/j,
(Edited by W.-K. Liu and T. Belytschko). Elmepress. Base1 (I 989).
Ancr1y.G.~
pp. 281-298.