JOURNAL OF THEORETICAL AND APPLIED MECHANICS 44, 3, pp. 585-602, Warsaw 2006 STATISTICAL DAMAGE MECHANICS – CONSTITUTIVE RELATIONS Antonio Rinaldi Dusan Krajcinovic Mechanical and Aerospace Engineering, Arizona State University, Tempe, USA e-mail: antnio.rinaldi@gmail.com Sreten Mastilovic Center for Multidisciplinary Studies, University of Belgrade, Belgrade, Serbia and Montenegro e-mail: misko3210@yahoo.com The statistical damage model presented by the authors in the previous pa- per of this series is used to formulate analytical constitutive relations for the hardening and softening phases of two-dimensional lattices. A proper defi- nition of the damage parameter for the softening is introduced. The results confirm that the analytical model can be used for the study of the softening phase and failure. This research offers a seminal basis forDamageTolerance Principles technology standards of the commercial airplane industry. Key words: failure, statistical damagemechanics, damage tolerance, damage parameter, constitutive relations, scaling 1. Introduction: damage tolerance principles Amaterial microstructure containingmany randomly distributedmicrocracks is initially statistically homogeneous, butbecomesheterogeneousdue topropa- gation of microcracks and their clustering into a macrocrack close to failure. Krajcinovic and Rinaldi (2005a) discuss the homogeneous-to-heterogeneous phase transition using the analytical tools of statistical mechanics, thermody- namics and fractal geometry.The thresholdof failuredependson the structural transformations of the material on the microscale and the structure size. The limit states design is driven by crack growth, but also by the structural re- arrangements in thematerial microstructure in the heterogeneous phase. The 586 A. Rinaldi et al. base of structural design and maintenance principles in Boeing Commercial Airplane Group, U.S.A., is the ”Damage Tolerance Principles” (Goranson, 1993), which is focused on two structural design objectives: 1. DamageTolerance: ability of structure to sustain anticipated loads in the presence of fatigue, corrosion or accidental damage until such damage is detected through inspections or malfunctions and repaired; 2. Durability: ability of the structure to sustain degradation from such so- urces as fatigue, accidental damage and environmental deterioration to the extent that they can be controlled by economically acceptable ma- intenance and inspection programs. Damage tolerance comprises three elements of importance for achieving the desired level of safety. The first element is the determination of the residu- al strength or the maximum allowable damage (including multiple secondary cracks) that the structure can sustain under regulatory fail-safe load condi- tions.The second element is the crack growthdefinedas the interval of damage progression from lengths with negligible probability of failure to an allowable size determined by the residual strength. Finally, a damage detection strate- gy (Inspection Program) must be adopted. The sequence of inspections in a fleet of airplanes requiresmethods and intervals selected to achieve timely the damage detection. Fig.1, adopted from Goranson (1993), shows typical experimental data from both full-scale crack growth testing (600 tests on two different wing- panels of width 200mmand 2300mm) and Linear Elastic FractureMechanics (LEFM).The crack length is normalized to theLEFMlimit Ly andthe streng- ths arenormalized to themaximumstrengthof pristineundamagedpanel.The maximum allowable damage (i.e. the minimum normalized strength) and the corresponding maximum defect length are assigned on such a graph in com- pliance with a ”fail-safe” strategy. Nevertheless, the choice of the tolerable damage is largely based upon experience and the Probability Of Detection (POD) from visual inspections. Full-scale tests are necessary to assess the ef- fect of the structural size on the damage tolerance, but Goranson recognizes the ”impossibility” of conducting full-size fracture and fatigue tests to obtain the data in Fig.1 for all the components. ”The emphasis on residual strength verification has gradually shifted in recent years from wing structures to fu- selage pressure shells” because ”the extended use of jet transport structures raised concerns about multiple site damage in fuselage structures” (Goran- son, 1993), which lead to expensive full-scale tests (of the order of millions of dollars) with large pressure test fixtures. Unfortunately, LEFM is not al- ways applied to multiple-site cracking and diffuse damage, which is nowadays Statistical damage mechanics – constitutive relations 587 still managed in a purely empirical manner. In conclusion, there is an urge to formulate more reliable multiscale analytical models that account for the structural size effect and which can be used for data extrapolation. The ca- pability of estimating the POD a priori from suchmodels would likely have a deep impact on the fail-safe strategy. Fig. 1. Sample data of crack growth experiments on wing panels used at Boeing Co. (Goranson, 1993); Ly – boundary between LEFM and transition behavior Goranson’s ”Damage Tolerance Principles” expresses the concern about the ”lack of interest” in the scientific community to characterize and quantify the POD. This paper shows how the ideas discussed in the first part of this series (Krajcinovic andRinaldi, 2005a), although in their infancy stage, could potentially address this concernandthe otherneeds of damage tolerant design. The time-independent damage model (Kachanov, 1958) σ = E o (1 −D)ε, relating themacrostress σ and themacrostrain ε, is obtained for the uniaxial tensile response of two-dimensional disordered spring networks through the damage parameter D and the corresponding scaling relations in Krajcinovic and Rinaldi (2005a). E o is the elastic modulus of pristine state and the bar above the symbols indicates a macroscale quantity. 2. Lattice simulations and numerical data The goal of this paper is to elucidate the universal behavior of the damage pa- rameters and to infer constitutive relations of the quasi-brittle materials. The 588 A. Rinaldi et al. disordered discrete model used for that purpose, although admittedly simple, captures what we believe to be the key features of the problem at hand. The approach taken in this paper is following footsteps of the extensivework on the central-force lattices modeling the brittle fracture of disorderedmaterials (for example,Curtin andScher, 1990;Hansen et al., 1989,Krajcinovic andBasista, 1991). The model does not aim at describing a specific material but rather a class of quasi-brittle materials, whose primary microstructural failure mode is rupture of interfaces of inferior strength (compared to that of the bulk). Examples include the monolithic polycrystalline ceramics with inferior grain- boundary strength, ceramic composites with inferior particle-matrix interface strength, polyphasematerials with weak phase interfaces, as well as concrete, clastic rocks, etc. The common feature of these materials, which is essential for adequacy of the proposedmicrostructural modeling, is that the interfacial (as an example, grain boundary1) fracture energy is significantly less than the cleavage energy of themicrotextural constituents (correspondingly, individual grains). These issues are discussed in detail in Davidge (1979), Lawn (1993), and are beyond the scope of the present work. The microtexture of the materials of this kind in two-dimensional space could be represented by a randomVoro- noi froth with the dual Delaunay lattice (Krajcinovic and Rinaldi, 2005a,b). A Voronoi polygon represents a grain of ceramic, a concrete aggregate or a granule of clastic rock whereas a bond in the Delaunay lattices is represen- tative of corresponding interface cohesion2. Damage evolution, which reflects accumulation ofmicro-events of degradation, is a stochastic process dependent on the disorder of the microstructure. The lattice is geometrically disordered since the equilibrium link lengths are normally distributed within the ran- ge αlλ ¬ λ ¬ (2−αl)λ with αl = 0.1 (if αl = 1 all grains are perfect 1Grain boundary is the most common example of the weak interface in brittle materials (Lawn, 1993). 2There are issues with representation of materials by spring-networkmodels that have been discussed in detail in the past (see, for example, Ashurst andHoover, 1976; Curtin and Scher, 1990; Jagota and Bennison, 1994). Nonetheless, we believe that these issues are not of crucial importance herein, bearing in mind the motivation of the work and the considered class of materials. In short, the rationale for the latter is that although adequacy of the spring-network model is not so obvious as for, say, granular materials, its application does not represent mere discretization of a continuum but rather reflects the discrete and disordered nature of the materials whose failure is governed by an ”ensemble”, a web, of weak interfaces. As Jagota and Berrensen (1994) conclude: ”This may not be an issue when dealing with the large disorder and when the goal is to study universal scaling relationships...” Statistical damage mechanics – constitutive relations 589 hexagons). Damage is introduced in the network by rupturing of the links, which represent interfacial microcrack formation. The links are nonlinear in compression and linear in tension, with random tensile strength and identi- cal stiffness k. If the critical tensile strain εcr is reached, permanent rupture occurs and the spring turns into a contact element. The tensile stiffness be- comes zero and the link cannot any longer carry tensile forces. Broken links remain active in compression if load reversal occurs in the course of defor- mation to account for crack closure. The values εcr are randomly sampled from a uniform distribution starting at zero. This lattice model considers interfacial (hereinafter referred to, more narrowly, as intergranular) micro- cracks only, which is a reasonable approximation for many ceramics (David- ge, 1979; Lawn, 1993). Since the resolution length of the model is equal to the grain facet, the rupture, i.e. the growth of a microcrack from the ini- tial length to the length of the grain boundary facet, is assumed to be in- stantaneous. The local fluctuations of the energy barriers (quenched within the material) and the local fluctuations of stress are essential features of the problem at hand, since localization is impossible in the absence of disorder (Anderson, 1958). Quasi-static displacement-controlled uniaxial tensile tests are simulated on different lattice sizes. Themolecular dynamics solver based upon the Ver- let’s algorithm (Allen and Tidesley, 1994; Krajcinovic and Rinaldi, 2005b; Mastilovic and Krajcinovic, 1999) was adopted. Each simulation is carried on incrementally up to the threshold of failure by applying small displa- cement steps and by computing the equilibrium configuration at each step. The damage process is tracked during the deformation by recording the num- ber of broken bonds n. The macroscopic data scatter of the F vs. u and n vs. u curves (within-size variability) indicates that D(ε,L) is a random va- riable at any given ε in the softening phase (Fig.2c,d). The average F vs. u and n vs. u curves from the 10 replicates per size N = {24,48,96,192} were considered for the scaling in Krajcinovic and Rinaldi (2005a), with N being the number of grains per lattice side. The original dataset is he- re sensibly expanded to enhance the accuracy, robustness and precision of the regression analysis. Intermediate lattice sizes N = {72,120} are ad- ded and more than 10 runs are collected for smaller lattices. The maxi- mum lattice size is also enlarged. Five extra simulations for N = 288 are analyzed but such runs are limited to the hardening phase only for the sake of computational feasibility. The simulation data are summarized in Table 1, for a total of more than 200 simulations as opposed to the original 40. 590 A. Rinaldi et al. Table 1.Number of runs per lattice size fromMD simulations (old Kraj- cinovic and Rinaldi (2005a) and new) Size N 24 48 72 96 120 192 288 Replicates 100 34 30 25 20 13 5 Fig. 2. 12×12 Voronoi/Delaunay graph of the microstructure in tensile test configuration (a); location of broken bonds during a simulation; localization and reduction of damage rate are evident after the transition (b); F vs. u (c) and n vs. u (d) curves for the 10 replicates of N =192 An example of the diffused uncorrelated damage and of the macrocrack from damage localization is displayed in Fig.3. Only the evolution of such macrocrack canprovide the”CrackGrowth”definedas”the interval ofdamage progression from length below which there is negligible POD to an allowable size determined by residual strength requirements” (Goranson, 1993). Statistical damage mechanics – constitutive relations 591 Fig. 3. Damage distribution withmacrocrack at the threshold of failure for N =192 3. Statistical (damage) mechanics The limits of traditional continuummodels of damage are discussed at length in literature, e.g. in Krajcinovic and Rinaldi (2005b). Models based upon Eshelby’s solution (Mura, 1982), representative volume element (Beran, 1968; Kestin, 1992; Nemat-Nasser and Hori, 1993) and fabric tensor (Krajcinovic, 1996), are suitable only as long as the microstructure is statistically homoge- nous. Only in this case, the expectation value of damage parameter D (the bar above the symbol indicates amacroscopic quantity) is equal to the volume average 〈D〉 of individual microcracks. Materials are homogeneous when mi- crocracks nucleate in absence of cooperative phenomena, and heterogeneous when clusters of microcracks form. The threshold of fracture is reached when one cluster reaches a correlation length equal to the lattice (specimen) size. Our simulation data indicate that the scaling techniques apply to damage mechanics (Krajcinovic and Rinaldi, 2005a). This statistical model of the di- sorderedmicrostructure provides analytical expressions of damage parameter. During the hardening phase, i.e. the part of the equation σ=E o [1−D(L)]ε corresponding to ∂σ­ 0, a simple analytical formula for the damage parame- ter was obtained from the scaled data D(ε,L)= aε+ b ε2 Lα (3.1) The values {α,a,b} = {−0.035,275,−14862} were deduced from simulation data. Instead, the analytical model for the damage parameter during the 592 A. Rinaldi et al. softening phase, i.e. the part of the equation σ=E o [1−D(L)]ε corresponding to ∂σ¬ 0, was D=LαDp+a1L αεs+ b1L αLz [ 1− exp ( − c1εs Lz )] (3.2) where εs = (ε − εpeak)/L α, Dp = 0.5 and {a1,b1,c1,z} = = {15.80,2.2,100,−0.52}3 from simulation data (ref. Fig.7 inKrajcinovic and Rinaldi (2005a)). 4. Damage parameter and constitutive relations The purpose of this section is the deduction of the analytical constitutive relation σ(ε,L)=E o [1−D(ε,L)]ε (4.1) from the correct application of the scaling relations (3.1) and (3.2). Thedama- ge parameter in the hardening phase was defined in Krajcinovic and Rinaldi (2005a) as D= n(ε) µ L2 (4.2) where n(ε) is the number of broken bonds and µ is a constant deduced from simulations. 4.1. Damage parameter for hardening: rationale Definition (4.1) is inspired by three observations. First, the strict simila- rity between the Parallel Bar System (PBS) and the lattice during damage nucleation (Krajcinovic, 1996; Krajcinovic and Rinaldi, 2005b) suggested the following form of damage parameter D= n 2Np (4.3) where Np is the number of broken links at the peak. Second, the log-log plot of the average Np vs. L in Fig.3b obtained from simulations demonstrates that Np for the lattice is a fractal following the power law Np(L)∼=n(εp,L)= e −1.3L∗1.94 (4.4) 3As noted in Krajcinovic and Rinaldi (2005a), only one of parameters a1, b1, and c1 is independent. Statistical damage mechanics – constitutive relations 593 The fractal exponent is 1.94 with 95% confidence interval [1.92,1.96]. From Eqs. (4.2)-(4.4) the value µ=2e−1.3 is estimated,where −1.3 is the (rounded) intercept of the regression line in Fig.4b. Such value of µ differs from the one reported in (Krajcinovic and Rinaldi, 2005a), previously set to µ = 2e−1.6. This discrepancy is due to the different datasets used for the statistical analy- sis. Analogously, the previous estimate of the fractal exponent was 1.96. Since the ”expanded” dataset is used here for the regression, the new estimates sho- uld be more reliable and accurate. In this respect, the mean square error and the confidence intervals on the regression coefficients (slope and intercept) are tighter. Nevertheless, most of our treatment will refer to the original dataset only and, hence,wewill retain the old value µ=2e−1.6 in order to legitimately use Eqs.(3.1) and (3.2). Fig. 4. Fractal behavior of n(ε,L) in the middle of hardening phase (exponent =2) and at the peak (exponent =1.94). 95% confidence intervals (dashed line) are very tight Finally, we replace L1.94 with L2 in Eq. (4.2) because n is a fat fractal with exponent 2 throughout the hardening phase before the transition/peak, as demonstrated by the n vs. L plot in Fig.4a (taken at an arbitrary ε away from the peak). Furthermore, the n vs. ε curves of different lattice sizes collapsed into a single curve D vs. ε during damage nucleation only for the exponent 2 (Fig.4b in Krajcinovic and Rinaldi, 2005a). This fact, that the relevant damage parameter in the hardening phase is proportional to the density of broken bonds, was already recognized by Hansen et al. (1989). The coefficients of determination R2 and R2adj reported by the statistical software MINITAB 14 c© are rounded to 100% in both graphs in Fig.4, which proves the high significance/quality of the regression analysis (R2 =1 indica- tes perfect linear correlation ([17]; Montgomery et al., 2001)). The regression on the seven sizes N = {24,48,72,96,120,192,288} provides a 2-decimal di- gits precision for the fractal exponents, which is satisfactory. The size range spans over one order of magnitude. 594 A. Rinaldi et al. 4.2. Damage parameter for softening With D defined as in (4.2), the D− ε curves were scaled in Krajcinovic and Rinaldi (2005a) to obtain Eqs. (3.1) and (3.2). While Eq. (3.1) is a su- itable definition of damage parameter for the hardening phase, Eq. (3.2) is not a valid damage parameter for the softening phase. The lattice is a spring- network and the corresponding secant stiffness E(L) = E o [1−D(ε,L)] is positive-definite (Krajcinovic, 1996). Hence, the ”proper damage parameter” mustalwaysbe0¬D¬ 1andshouldmonotonically approachD=1at failure (Efailure = 0), either discontinuously or continuously depending on whether the failure transition is of thefirst order or second order.These two constraints were not accounted for in the scaling of the softening data. Therefore, while expression (3.1) is perfectly legitimate, definition (3.2) is still incomplete. The number of broken links n(ε) is an indicator of the total damage. The parameter D supplies just the normalization, like in Eq. (4.3), and po- ssibly renders this information scale invariant. One fundamental idea in our framework is that hardening and softening phases are independent sequential processes. The microcracks nucleation is driven by a fat fractal with fractal exponent 2 for the most part and by a proper fractal with exponent 1.94 at the peak of stress-strain response. Such driving set changes at the transition, i.e. an appropriate denominator (normalization number) should replace µL2 in (4.2). The following form of damage parameter is assumed in the softening phase D(ε,L,n)= n(εpeak,L) µL2 + ∆n(ε,L) X(L) (4.5) where ∆n = n(ε,L)−n(εpeak,L) is the number of broken bonds in the so- ftening phase. The analytical expression for ∆n follows from (3.2) and (4.2) as ∆n(ε,L)=µL2(D−LαDp)=µL 2 { a1L αεs+b1L αLz [ 1−exp ( − c1εs Lz )]} (4.6) which is essentially the information conveyed by Eq. (3.2). Eq.(4.6) is equiva- lent tomagnifying the scaled data in the softening phase by a factor µL2. The constitutive relations for the hardening and softening phases can be rewritten as σ=    E o [ 1− n(ε,L) µL2 ] ε for ε¬ εpeak E o [( 1− Np µL2 ) − ∆n(ε,L) X(L) ] ε for ε> εpeak (4.7) Then, X(L) in thedenominator of (4.5) and (4.7)2 iswhatneeds tobedefined. The 3D plot of microcracks location vs. simulation time in Fig.2b highlights Statistical damage mechanics – constitutive relations 595 the contrast between the uncorrelated and uniformly distributed damage nuc- leation and the highly correlated and localized damage propagation. From this viewpoint, a fractal exponent close to one makes intuitive sense for the invariant set driving the propagation, as much as a fat fractal suits the phe- nomenology of the nucleation. A deductive reasoning based upon statistical methods is used to reach the data driven choice for X(L). By assuming that (4.7)2 is the right form of the model, one can write σ̂i(εi,∆ni,β) =E o [( 1− Np µL2 ) −β∆ni ] εi (4.8) where β(L)= 1/X(L), Q is the numberof numerical observations used for the regression and σ̂i is the estimate of stress value at εi and ∆ni from themodel for i=1,2, . . . ,Q. Since themodel is linear in the unknownparameter β, the ”ordinary least squares”method can be used to compute aminimumunbiased estimator β̂ of β (Montgomery et al., 2001). Themean square error function for this model is Err(L)= Q∑ i=1 (σi− σ̂i) 2 = Q∑ i=1 ( σi−E o [( 1− Np µL2 ) − β̂∆ni ] εi )2 (4.9) which is the sum of the squared deviations between the simulation data and the predicted value over Q. The error function is minimized with respect to β̂ by setting ∂Err/∂β̂ =0. The calculations lead to β̂(L)= ∑Q i=1 [( 1− Np µL2 ) E o εi−σi ] ∆niεi ∑Q i=1E o ε2i∆n 2 i (4.10) whichdependson the lattice sizeL. Finally, the formula for the estimate X̂(L) is X̂(L)= β̂−1(L)= ∑Q i=1E o ε2i∆n 2 i ∑Q i=1 [( 1− Np µL2 ) E o εi−σi ] ∆niεi (4.11) The results of Eq. (4.11) are marked in Fig.5 by circles for the sizes N = {24,48,72,96,120,192}. The linear fit through the six X̂(L) valuesmat- ches the results very well, as indicated by a coefficient of determination close to unity (R2 =0.9971; R2adj =0.997). The equation of the regression line is X(L)=µ2L+µ3 =65.07L−1038 (4.12) 596 A. Rinaldi et al. Fig. 5. Comparison of X(L) estimates from (4.11) for the 6 average ∆n (circles) and for 6 random replicates per each size N = {24,48,72,96,120,192} (asterisks) with 95%confidence intervals (60.38,69.76) and (−1589,−585.6) for the slope and the intercept respectively. This linear model supplies X(L) in Eqs. (4.5) and (4.7)2. Such a result seems to confirm the above supposition from phenomenolo- gy about the linear dependence of X on L. Nevertheless, more data points are required to increase the precision and tighten the relatively wide confi- dence intervals. For very large lattices, when the intercept becomes negligible (1038 � 65.07L), the relations (4.5) and (4.7)2 are nicely approximated for ε> εpeak by σ=E o [( 1− Np µL2 ) − ∆n(ε,L) 65L ] ε (4.13) D(ε,L,n)= n(εpeak,L) µL2 + ∆n(ε,L) µ2L = Np µL2 + ∆n(ε,L) 65L Relation (4.13)2 implies a discontinuity in the derivative of the damage pa- rameter D(ε,L,n) due to the abrupt change of the denominator from L2 to L at the hardening-softening transition. While n(ε) is a C1 continuous func- tion at εpeak, the damage is C 0 continuous and the damage rate ∂n(ε)/∂ε discontinuous, in agreement with Fig.2b and Fig.2d. The choice of the Q data points to be used in the evaluation of (4.11) is crucial for the accuracy of X̂(L) anddependsonwhatpart of the response is of interest. This analysis aimed at capturing the steepest descent of the softening regime after the transition and only the corresponding points (about 25% of softening dataset) were selected tominimize themodel error in that part. The performance of the analytical expressions (4.7) are compared in Fig.6a and Statistical damage mechanics – constitutive relations 597 Fig.6b against the original average simulation data for N =48 and N =96, respectively. Fig. 6. Overall model fitting for N =48 and N =96 with X fromOLS, Eq. (4.11); n(ε,L) and ∆n(ε,L) are estimated via scaling relations. (c) Improved fit for N =96 when numerical data are used directly for n(ε,L) and ∆n(ε,L) instead of scaling relations in the constitutive relations (4.7) 4.3. Interpretation and limits of relations for softening The model closely estimates the numerical data for most of the design space.Concernsmight arise about the accuracy of relations (4.5) and (4.7)2 for the softening phase. The performances of themodel are determined essentially by two aspects: 1. The selected model form (4.5) for the damage parameter D 2. The accuracy of the estimates n(ε,L) and ∆n(ε,L) from scaling. At first, by assuming that the model form (4.5) is correct, we focus on the second issue. The scaling relations for n(ε,L) and ∆n(ε,L) from Eqs. (3.1) and (4.6) refer to the average data, obtained by first averaging the original 598 A. Rinaldi et al. n− ε curves of individual replicates for each size and then, after scaling, by fitting a regression model to the scaled data. This twofold ”filtering” process smeared out the irregularities characteristic of each replicate (compare Fig.2d with Fig.7 inKrajcinovic andRinaldi (2005a). Thus, while on one sidewe ob- tain smooth analytical relations capable of estimating the averagemicrocracks number for any lattice size, on the other side the characteristic ”details” of each individual replicate are lost. Fig.6c shows the model obtained fromEqs. (4.6) for N =96 when average n and ∆n from simulations are used in place of the scaling relations. The new estimates ”shadow” the simulation data of the softening phase better than before. Evidently a better knowledge of n and ∆n leads to a more accurate constitutive model. For these simulations, the difference between the estimates of n and ∆n from scaling relations and from simulations, never exceeded 10% for any lattice size. However, the predictive power does not change significantly. To assess the other issue about the appropriateness of the selected model form, Eqs. (4.5) and (4.7) are tested against individual lattice simulations. Six individual replicates for N = {24,48,96,72,120,192} are randomly pic- ked. This time the values of n, ∆n and Np are obtained directly from the simulation. As the model is satisfactory for the hardening phase, the previo- us estimate of µ is kept for all cases. Instead, the parameter X(L) in (4.5) is re-estimated for each replicate via Eq. (4.11). These new values of X̂(L) are marked in Fig.5 by asterisks and fall reasonably close to the line (4.12) for the average values. The corresponding six constitutive models are shown in Fig.7. The agreement between models (full line) and numerical data (dot- ted) is remarkable. Complex patterns and discontinuities arewell-captured for most of the softening phase. Themodels are in general able to reproduce the response also after the change of curvature following the initial steep descen- dent from the peak. These results are achieved by re-estimating merely the sole parameter X(L). The conclusions are the following: • The choice of the model form in Eq. (4.5) seems appropriate • A simple model can reproduce a variety of softening responses • The greatest accuracy is achieved by calibrating the constitutive rela- tions on an individual simulation; the same constitutive relations of the mean response can be used by simply re-estimating X(L) from the si- mulation data • A greater accuracy for the individual replicate is determined by the usage of un-filtered simulation data for n and ∆n in place of the average estimates. Statistical damage mechanics – constitutive relations 599 Fig. 7. Responses of customizedmodel (4.5) and (4.7)2 (full line) vs. simulation data (dotted line) for six individual random replicates for N = {24,48,72,96,120,192} The findings corroborate and clarify the proposed constitutive model. The last conclusion pointed out that the continuous approximation of n(ε,L) and ∆n(∆ε,L) from scaling relations overlooks large discontinuities in the n− ε data registered in connectionwith large avalanches, both at the onset of dama- ge localization and afterwards. The lattice response of the softening phase has a saw-toothed appearance on the macro-scale due to large avalanches, whose signatures are sudden vertical drops in the stress level alternated to linear ela- 600 A. Rinaldi et al. stic ”climbs”. Such vertical drops observable for each replicate (Fig.7) totally disappear in the average curves inFig.??a, leaving space to a ”fictitious”more gentle descendent. An accurate model matching the data from an individual test is achievable when the specific details of ∆n(∆ε,L) are accounted for on a case-by-case basis. Conversely, the scaling relations (3.1) and (3.2) represent an average behavior and offer the general ”trend” for all lattice sizes. The results outline the potential of a statistical damage mechanics model, which go well beyond the limits of traditional continuum approach. 5. Conclusion and summary The statistical mechanics theory and the application of the theory have been presented in thepreviouspaper (Krajcinovic andRinaldi, 2005a).Thispaper is committed to constitutivemodeling, failure anddamage tolerance. The results offered in this paper display the effectiveness of the statistical damage model in the study of the softening regime and failure. The numerical data produced with this approach could not be generated by classical continuum damage mechanics. The scale-invariant damage curves proposed in Krajcinovic and Rinaldi (2005a) are re-examined and finalized herein to obtain the analytical consti- tutive relations (4.7) for the entire damage process. The damage parameter is properly defined not only in the hardening phase but also in the softening phase to attain this result. In the hardening phase the process is driven by a fractal set of exponent two, while in the softening phase the fractal dimension tends towards one. Additional numerical data improved the precision of the statistical analysis. The present paper not only reorganizes and reinforces so- me of the ideas originated in the inspiring work of Hansen et al. (1989), but also addresses successfully some of the open questions by adopting a rather different methodology. This study outlines the importance of the statistical methods, such as ordinary least squares, maximum likelihood and hypothesis testing, for the selection of the model parameters. These techniques are the basis for data-driven reasoning and decision making in damage tolerance. Acknowledgments This research is sponsored by theMathematical, Information andComputational ScienceDivision,OfficeofAdvancedScientificComputingResearch,U.S.Department of Energy under contract numberDE-AC05-00OR22725withUT-Battelle, LLC. The support of theMinistryof Science andEnvironmentalProtectionofRepublic of Serbia Statistical damage mechanics – constitutive relations 601 (under contract numbers 1793 and 1865) to S.Mastilovic is gratefully acknowledged. The authors address a special thank toDr S. Simunovic atORNL and prof. Y.C. Lai at ASU for valuable support. References 1. Allen M.P., Tildesley D.J., 1994, Computer Simulation of Liquids, Cla- rendon Press, Oxford, UK 2. Anderson P.W., 1958, Absence of diffusion in certain random lattices,Phys. Rev., 109, 1492-1505 3. Ashurst W.T., Hoover W.G., 1976, Microscopic fracture studies in the two-dimensional triangular lattice, Phys. Rev. B, 14, 4, 1465-1473 4. Beran M., 1968, Statistical Continuum Theories, Wiley, NewYork 5. Curtin W.A., Scher H., 1990, Brittle fracture in disordered materials, J. Mater. Res., 5, 3, 535-553 6. Davidge R.W., 1979, Mechanical Behavior of Ceramics, Cambridge Univ. Press, Cambridge, UK 7. Goranson U.G., 1993, Damage tolerance – facts and fiction, 14th Plantema Memorial Lecture Presented at the 17th Symposium of the International on Aeronautical Fatigue, Stockholm, Sweden 8. Hansen A., Roux S., Herrmann H.J., 1989, Rupture of central-force latti- ces, Phys. Fracture, 50, 517-522 9. Jagota A., Bennison S.J., 1994, Spring-network and finite element models for elasticity and fracture,Proceedings of a Workshop on Breakdown and Non- linearity in Soft Condensed Matter, K.K. Bardhan, B.K. Chakrabarti, A. Han- sen (edit.), Springer-Verlag Lecture Notes in Physics, Berlin, Heidelberg, New York 10. KachanovL.M., 1958,On the time to failure under creep conditions, Izv. AN SSSR, Ot. Tekhn. Nauk, 8, 26-31 11. Kestin L., 1992, Local-equilibrium formalism applied to mechanics of solids, Int. J. Solids Structures, 29, 14/15, 1827-1836 12. Krajcinovic D., 1996, Damage Mechanics, North-Holland Series in Applied Mathematics andMechanics, Vol. 41, Elsevier, Amsterdam 13. KrajcinovicD.,BasistaM., 1991,Ruptureof central-force lattices revisited, J. de Physique I, 1, 241-225 602 A. Rinaldi et al. 14. Krajcinovic D., Rinaldi A., 2005a, Statistical damagemechanics – 1. The- ory, J. of Applied Mechanics, 72, 76-85 15. KrajcinovicD.,RinaldiA., 2005b,Thermodynamics and statistical physics of damage processes in quasi-ductile solids,Mech. Mater., 37, 299-315 16. LawnB., 1993,Fracture of Brittle Solids, CambridgeUniv. Press, Cambridge, UK 17. MINITAB 14 c©, on-line manual 18. Mastilovic S., KrajcinovicD., 1999, Statistical models of brittle deforma- tion. Part II: Computer simulations, Int. J. Plasticity, 15, 427-456 19. MontgomeryD.C., PeckE.A.,ViningG.G., 2001, Introduction to Linear Regression Analysis, Wiley 20. Mura T., 1982,Micromechanics of Defects in Solids, Martinus Nijhoff Publi- shers, The Hague 21. Nemat-Nasser S., Hori M., 1993, Micromechanics: Overall Properties of Heterogeneous Materials, North-Holland Series in Applied Mathematics and Mechanics, Vol. 37, Elsevier, Amsterdam Statystyczna mechanika uszkodzenia - związki konstytutywne Streszczenie Statystycznymodel uszkodzenia, przedstawiony przez autoróww innej pracy, za- stosowano tu do sformułowania analitycznych związków konstytutywnych w zakresie wzmocnienia i osłabienia dla dwuwymiarowychmodeli sieciowych.Wprowadzono od- powiednią definicję parametru uszkodzenia w zakresie osłabienia. Uzyskane wyniki potwierdzają, że wprowadzony model analityczny można stosować do badania osła- bienia i zniszczenia materiałów. Praca stanowić może podstawę do ustalenia norm technicznych uszkodzenia stosowanychw przemyśle lotniczym Manuscript received October 24, 2005; accepted for print March 27, 2006