ML042600455

ORNL Report, ORNL/NRC/LTR-04-15, Probabilistic Structural Mechanics Analysis of the Degraded Davis Besse RPV Head.

ML042600455
Person / Time
Site:	Davis Besse
Issue date:	09/01/2004
From:	Bass B, Williams P, Yin S Oak Ridge
To:	Matthew Kirk Division of Engineering Technology
References
DE-AC05-00OR22725, DOE 1886-N653-3Y, Job Code Y6533
Download: ML042600455 (139)

Text

ORNL/NRC/LTR-04/15 Contract Program or Heavy-Section Steel Technology (HSST) Program Project

Title:

Subject of this Document: Probabilistic Structural Mechanics Analysis of the Degraded Davis-Besse RPV Head Type of Document: Letter Report Authors: P. T. Williams S. Yin B. R. Bass Date of Document: September 2004 Responsible NRC Individual M. T. Kirk and NRC Office or Division Division of Engineering Technology Office of Nuclear Regulatory Research Prepared for the U. S. Nuclear Regulatory Commission Washington, D.C. 20555-0001 Under Interagency Agreement DOE 1886-N653-3Y NRC JCN NoY6533 OAK RIDGE NATIONAL LABORATORY Oak Ridge, Tennessee 37831-6085 managed and operated by UT-Battelle, LLC for the U. S. DEPARTMENT OF ENERGY under Contract No. DE-AC05-00OR22725

ORNL/NRC/LTR-04/15 Probabilistic Structural Mechanics Analysis of the Degraded Davis-Besse RPV Head P. T. Williams S. Yin B. R. Bass Oak Ridge National Laboratory Oak Ridge, Tennessee Manuscript Completed - May 2004 Date Published - September 2004 Prepared for the U.S. Nuclear Regulatory Commission Office of Nuclear Regulatory Research Under Interagency Agreement DOE 1886-N653-3Y NRC JCN No. Y6533 OAK RIDGE NATIONAL LABORATORY Oak Ridge, Tennessee 37831-6085 managed and operated by UT-Battelle, LLC for the U. S. DEPARTMENT OF ENERGY under Contract No. DE-AC05-00OR22725

CAUTION This document has not been given final patent clearance and is for internal use only. If this document is to be given public release, it must be cleared through the site Technical Information Office, which will see that the proper patent and technical information reviews are completed in accordance with the policies of Oak Ridge National Laboratory and UT-Battelle, LLC.

This report was prepared as an account of work sponsored by an agency of the United States government. Neither the United States government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise, does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States government or any agency thereof.

Probabilistic Structural Mechanics Analysis of the Degraded Davis-Besse RPV Head P. T. Williams, S. Yin, and B. R. Bass Oak Ridge National Laboratory P. O. Box 2008 Oak Ridge, TN, 37831-6085 Abstract The Heavy-Section Steel Technology (HSST) Program at Oak Ridge National Laboratory (ORNL) has performed a probabilistic structural mechanics (PSM) analysis of the damaged Davis-Besse reactor pressure vessel head in support of the United States Nuclear Regulatory Commissions ongoing forensic investigations. This report documents the results of that PSM analysis, including a description of the Davis-Besse wastage-area damage model, the technical basis for the model, and the results of sensitivity studies based on a cladding capacity analysis (CCA) and an Accident Sequence Precursor (ASP) investigation of the wastage cavity. A companion report describes the HSST experimental program carried out at ORNL in parallel with the PSM analysis.

The CCA and ASP studies provide approximate answers to three questions regarding the Davis Besse event:

(1) What applied pressure would have failed the wastage-cavity cladding at the time of discovery (TOD) on 16 February 2002? (CCA)

(2) How much longer could the Davis-Besse RPV have continued in service without failure of the pressure boundary if the wastage cavity had not been discovered on 16 February 2002? (CCA)

(3) Including uncertainties in the As-Found damage state of the wastage cavity, what was the probability of failure one year before the TOD, and how do these uncertainties affect the estimated probability of failure at TOD? (ASP)

The answer to question No. #1 required the construction of a detailed finite-element model (FEM) of the Davis Besse wastage cavity which incorporated the results of extensive laboratory measurements and metallographic examinations of the damaged site after it had been removed from the RPV head. The Davis-Besse cladding material was carefully characterized in terms of both strength (plastic-flow properties) and fracture toughness (ductile-tearing initiation and flaw growth). The fracture-toughness characterization (at a service temperature of 600 °F) was carried out by the HSST Program using pre-cracked Charpy V-Notch specimens taken from the Davis-Besse wastage cavity. All of the characterization studies included investigations of the uncertainties in the property measurements.

The results of the deterministic FEM analysis indicate that, for the most conservative assumptions made in the study regarding flaw size and depth, the estimated failure pressures all exceeded the relief-valve set-point pressure of 2500 psi. This result is in agreement with the forensic finding that exposure of the wastage cavity to the nominal operating pressure of 2165 psi did not produce any evidence of ductile crack initiation at the TOD. Median pressures needed to fail the wastage cavity (by ductile-tearing initiation of the Model Flaw) were estimated to range between 3000 and 5200 psi, representing a 1.4 to 2.4 margin against the operating pressure. A 90 percent confidence interval covering the median estimates was 2710 psi to 6500 psi (or 1.25 to 3.0 margin against the operating pressure).

Question No. #2 was addressed with a PSM analysis that reflected the very limited state-of-knowledge of how the wastage cavity might be expected to evolve over time beyond the known damage state at the TOD. An Expert Elicitation was carried out by the NRC staff to provide estimates for the wastage-cavity growth rate and the rate of stress-corrosion crack development due to exposure of the unbacked cladding iii

to the concentrated boric acid solution inside the wastage cavity. The resulting best-estimate predicts that the cumulative probability of survival decreases to 50% after approximately 230 days of additional operation. When applying the most conservative flaw model, this median time decreased to approximately 150 days of additional operation. In these studies, the consequences of failure were expressed in terms of a range of break sizes leading to a loss-of-coolant accident (LOCA). Additionally, the results from the CCA sensitivity study can be used to provide an approximate 90 % confidence interval covering the best estimate for the median failure of the cladding. These results predict that, at a confidence level of 90%, Davis-Besse could have continued operating from 2 to 22 months before cladding failure would be expected.

Question No. #3 was addressed by a modification to the PSM analysis that incorporated additional uncertainty in damage state of the wastage cavity at the TOD. At the TOD, the best-estimate probability of cladding failure was estimated to be approximately 20%. At 1 year before the TOD, the ASP analysis estimated a low probability of failure of approximately 1%.

iv

CONTENTS Section Page 1 Introduction........................................................................................................................................... 1 1.1 Method of Analysis ..................................................................................................................... 1 1.2 Background of the Davis-Besse Degraded RPV Head Problem ................................................. 2 1.3 Accident Sequence Precursor Program ....................................................................................... 3 1.4 Previous ORNL Studies of the Davis-Besse RPV Head Degradation Problem .......................... 6 1.4.1 Detailed Finite-Element Global and Submodels of Wastage Area and Cavity....................... 6 1.4.2 Results of Wastage-Cavity Growth-Pattern Study [13] .......................................................... 7 1.5 Scope of this Report .................................................................................................................. 12 2 Davis-Besse Wastage-Area Damage Model ....................................................................................... 13 2.1 Detailed 3D FEM Model for Deterministic Analyses ............................................................... 13 2.1.1 Modeled Wastage Cavity Footprint ...................................................................................... 13 2.1.2 Wastage Cavity ..................................................................................................................... 13 2.1.3 Modeled Surface Flaws......................................................................................................... 18 2.2 Simplified Wastage-Area Damage Model for Monte Carlo Simulations.................................. 19 2.2.1 Damage Model Geometry, Loading, and Constraint Conditions .......................................... 21 2.2.2 Damage-State Parameters ..................................................................................................... 21 2.2.3 Failure Mechanisms - Development of a Plastic-Collapse Fragility Curve ......................... 23 2.2.4 Failure Mechanisms - Development of a Ductile-Tearing Model........................................ 42 2.2.5 Assumed Wastage-Area Damage Model Accident Sequence............................................... 67 2.2.6 Results of an Expert Elicitation for Input Distributions........................................................ 79 2.2.7 Best-Estimate, More Conservative, and Less Conservative Sampling Distributions............ 93 2.2.8 Variance Reduction - The Method of Antithetic Variates.................................................... 93 2.2.9 LOCA Screening Rules......................................................................................................... 96 2.2.10 LOCA Break Size Definitions ........................................................................................ 100 2.2.11 Cavity Growth-Shape Scaling Rules .............................................................................. 100 3 Results and Discussion ..................................................................................................................... 101 3.1 Results of Deterministic FEM Analysis .................................................................................. 101 3.2 Cladding Capacity Analysis (CCA) - Sensitivity Study Results ............................................ 104 3.2.1 Convergence of Monte Carlo Simulations .......................................................................... 104 3.2.2 Best-Estimate CCA Results ................................................................................................ 107 3.2.3 Results of the CCA Sensitivity Study ................................................................................. 107 3.3 Accident Sequence Precursor Analysis - Best Estimate and Sensitivity Study Results ......... 111 3.3.1 Case Matrix for ASP Study................................................................................................. 118 3.3.2 Best-Estimate ASP Results ................................................................................................. 118 3.3.3 Results of the ASP Sensitivity Study .................................................................................. 118 4 Summary and Conclusions................................................................................................................ 124 References................................................................................................................................................. 126 Appendix A............................................................................................................................................... 130 Appendix B ............................................................................................................................................... 135 Appendix C ............................................................................................................................................... 270 v

LIST OF TABLES Table Page

1. Plastic-Flow Properties for SS Cladding ................................................................................................ 26

2. Goodness of Fit Statistics for Yield and Ultimate Stress Statistical Distributions.............................. 26

3. Results of Bivariate Sampling Protocols for a Range of Correlation Coefficients, ................................ 40

4. Ductile-Tearing Data Extracted from Table 13 of NUREG/CR-5511 [35]. ........................................... 45 5a. Ductile-Tearing Data Used in Development of JIc Statistical Distributions ......................................... 46 5b. Statistical Distributions Fitted to JIc Data from Three Sources............................................................. 46

6. Summary of Candidate Input Sampling Distributions ............................................................................ 81

7. LOCA Sizes Defined by Effective Break Sizes.................................................................................... 100 8a. DB Submodel Pressures to Initiate Ductile Tearing for the 5th, 50th , and 95th Percentiles with Margins Against Ductile Initiation at the Operating Pressure of 2165 psi ....................................... 102 8b. DB Submodel Pressures to Initiate Ductile Tearing for the 5th, 50th , and 95th Percentiles with Margins Against Ductile Initiation at the Relief-Valve Setpoint Pressure of 2500 psi .................... 102

9. Case Matrix for Cladding Capacity Analysis (CCA)............................................................................ 105

10. Monte Carlo Results - Summary of LOCA Probabilities (N = 50,000) ............................................. 105

11. Times to Failure after TOD................................................................................................................. 108

12. Case Matrix for ASP Sensitivity Study............................................................................................... 119

13. Case Matrix Layout Key ..................................................................................................................... 119

14. Summary of LOCA Probabilities for ASP Sensitivity Study ............................................................. 119 vi

LIST OF FIGURES Figure Page

1. Depictions of Davis-Besse RPV and wastage cavity: (a) Davis-Besse Nuclear Power Station RPV [7]; (b) schematic of a typical nuclear power reactor showing the relationship of the CRDM nozzles to the RPV head and the location of the damaged area [7]; ........................................ 4

1. (continued) (c) sketch of RPV head degradation; and (d) photograph presenting a plan view of the wastage-area footprint showing remaining J-groove weld, irregular topography of exposed unbacked cladding, and irregular sidewalls of the wastage cavity [7]. ................................................. 5

2. Finite-element global model and submodels applied in previous analyses [11,13] of the Davis-Besse head and wastage area................................................................................................................. 8

3. Geometry of the hemispherical RPV head and closure flange used in the global model [11,13]......... 8

4. Geometry of submodel relative to nozzles #3, #11, #15, and #16 in the RPV head. ............................ 9

5. Wastage-cavity submodel [11,13] developed from a ProEngineer solid model and then meshed with MSC Patran................................................................................................................... 9

6. A computational study was carried out in [13] to determine the sensitivity of burst pressures, predicted by finite-element analyses, to the geometry of estimated growth patterns of the wastage cavity footprint: (a) self-similar and ellipsoidal growth patterns and (b) circular burst disk. ..................................................................................................................................................... 10

7. The FEM-predicted burst pressures of irregularly-shaped footprint growth patterns are compared in [13] to burst pressure predictions for circular diaphragms under lateral pressure loading based on an instability theory proposed by Chakrabarty and Alexander [15]. ...................... 11

8. Davis-Besse global model of RPV closure head for the current study includes a simplified wastage-cavity definition and inner cladding layer............................................................................. 14

9. Davis-Besse submodel includes: (a) nozzles #3, #11, #15, #16 (Alloy 600); (b) a refined definition of the wastage-cavity morphology in the base material (A 533B pressure vessel steel), ... 15

9. (continued) (c) wastage cavity footprint; (d) irregular geometry of wastage cavity walls based on geometric measurements from the dental mold taken by BWXT [16]. ......................................... 16

9. (continued) (e) irregular exposed cladding topography based on a sinusoidal waveform, and (f)

J-groove weld near a modeled surface flaw. ....................................................................................... 17

10. Geometry, loading, and constraint conditions for the Davis-Besse RPV head damage model - a 1/4 symmetry burst disk with a centered surface flaw, where h0 = 0.25 in. and L = 2.0 in................... 20

11. Spherical geometry of deformation assumed in Chakrabarty and Alexanders plastic instability theory. ................................................................................................................................................. 22

12. Plots of true stress vs. true strain from the data in Table 1. ............................................................... 27

13. Probability density functions for fitted statistical marginal distributions based on data in Table 1: (a) yield true stress and (b) ultimate true stress............................................................................... 28

14. Statistical distributions approximating the material variability in the yield and ultimate stresses of cladding at 550 °F to 600 °F: (a) yield true stress CDFs with adjustments to match Davis-Besse best-estimate value of 30.9 ksi and (b) ultimate true stress CDFs with adjustments to match Davis-Besse best-estimate value of 72 ksi. .............................................................................. 33

15. Frequency histograms of sampled fragility curve slopes based on (a) perfect correlation between yield and ultimate stresses, (b) bivariate lognormal joint distribution with correlation coefficient of 0.06931, (c) bivariate gamma joint distribution with correlation coefficient of 0.06931, and (d) uncorrelated yield and ultimate stresses and, therefore, statistically independent. The correlation coefficient of 0.06931 was statistically inferred from the data in Table 1........................ 40 vii

LIST OF FIGURES (continued)

Figure Page

16. The intersection of the fragility curve with the y-axis is equal to the initial undeformed cladding thickness, h0 , which is assumed known in the current analysis: (a) a sampling of fragility curves based on their cumulative probability resulting from the bivariate lognormal joint distribution for yield and ultimate stress and (b) damage state partitioning based on an uncertain fragility curve slope with cumulative probability of 97.5%. ............................................................................ 41

17. Net-section plastic collapse pressures estimated by modified C&A (1970) theory compared to PNI values determined from 3D finite-element ABAQUS solutions: (a) a/t = 0.05, (b) a/t = 0.125, (c) a/t = 0.25, (Ductile tearing not simulated for these cases.) ................................. 43 17 (continued) Net-section plastic collapse pressures estimated by modified C&A (1970) theory compared to PNI values determined from 3D finite-element ABAQUS solutions: (d) a/t = 0.5, and (e) combined plot. (Ductile tearing not simulated for these cases.) ............................................. 44

18. Ductile-tearing data for cladding from several sources: (a) JIc data and (b) average tearing modulus data. ...................................................................................................................................... 47

19. Statistical distributions for JIc : (a) probability density of lognormal distribution fitted to the 12 Davis-Besse PCCVN data points and (b) lognormal cumulative distribution function compared to cumulative probabilities of Davis-Besse PCCVN JIc data estimated by the median rank order statistic p = (i-0.3)/(n+0.4). ................................................................................................................. 50

19. (continued) fitted distributions for separate samples of JIc data from 0.5T C(T) clad overlay (NUREG/CR-5511 and 6363) and PCCVN specimens from PVRUF and Davis-Besse cladding materials: (c) probability densities and (d) cumulative distribution functions.................................... 51

20. JR curves (a) fitted from data available in ref. [1] for unirradiated 0.5TCS fracture specimens at 550 °F and (b) comparing NUREG/CR C(T) JR curves with PCCVN JR curves................................ 53

21. Finite-element models used in calculating applied J-integrals produced by pressure loading of burst disk: (a) (a/t = 0.5, 2L/a = 16) (b) (a/t = 0.25, 2L/a = 32), and (c) (a/t = 0.125, 2L/a =

64) ....................................................................................................................................................... 54

21. (continued) Finite-element models used in calculating applied J-integrals produced by pressure loading of burst disk: (d) (a/t = 0.05, 2L/a = 160) (e) (a/t = 0.8, 2L/a = 10) .................................... 55

22. J-integral driving forces - applied pressure as a function of JI and a/t for a 2-inch long flaw........... 56

23. A-tip tearing instability plots with JR vs TR curves from ref. [35] and Japplied vs Tapplied curves from finite-element simulations of burst disks with 2-inch long surface flaws centered in the burst disk.. ........................................................................................................................................... 57

24. A-tip tearing instability plots with JR vs TR curves from ref. [35] and Japplied vs Tapplied curves from finite-element simulations of burst disks with 2-inch long surface flaws centered in the burst disk. Applied J-T curves are for a range of flaw depths ( 0.0125 in. a 0.125 in. ) and varying applied pressures.................................................................................................................... 58

25. C-tip tearing instability plots with JR vs TR curves from ref. [35] and Japplied vs Tapplied curves from finite-element simulations of burst disks for As-Found flaw with initial 2-inch long with A-tip collapse. ..................................................................................................................................... 58

26. Driving-force curves used in C-Tip ductile-tearing instability analyses: (a) models were developed for varying Rcavity and constant 2L = 2a = 2.0 in, tclad = 0.25 in., (b) applied driving force Japplied at C-tip at 2.165 ksi for varying Rcavity, and (c) applied dJ/da at C-tip at 2.165 ksi for varying Rcavity. ...................................................................................................................................... 60

27. Given a JR curve in power-law model form and estimated flow stress, f , the initiation toughness, JIc , and local tearing modulus, TR, are uniquely defined (see ASTM E-1820)................. 61

28. Lognormal distribution characterizing the uncertainty in the JR curve power-law exponent, m. ....... 63 viii

LIST OF FIGURES (continued)

Figure Page

29. Plot of Japplied as a function of Rcavity and aA-tip for a burst disk with a 2 inch centered flaw under 2.165 ksi applied pressure. Japplied solutions from FEM models (using best-estimate stress vs strain data) were fitted to a surface function by TableCurve 3D. .................................................... 66

30. Time line for cladding capacity analysis - accident sequence of wastage-area damage model. ....... 68

31. For the cladding-capacity analysis, the As-Found damage state is assumed known and fixed. For constant cavity and flaw growth rates, the accumulation of damage follows, in the absence of ductile tearing, a linear path towards the fragility curve................................................... 68

32. An example flaw growth history (with ductile tearing) calculated using MathCad. The case conditions are based on the median values from the more conservative (MC) sampling distribution for dR / d , best-estimate distribution for da / d , JIc, and m. ........................................ 72

33. Evolution of the (a) driving force, Japplied, and (b) tearing modulus, Tapplied for the example problem. .............................................................................................................................................. 75

34. Failure-assessment diagram using the damage-state growth path calculated for the example problem. .............................................................................................................................................. 76

35. Schematic descriptions of the model flaws (a) A-tip and C-tip before A-tip ligament collapse and (b) C-tip location after A-tip ligament collapse. .......................................................................... 77

36. Candidate sampling distributions for cavity size at one year before discovery, 0 : (a) probability densities and (b) cumulative probabilities. (ASP analysis only)......................................................... 83

37. Candidate sampling distributions for cavity wastage rate, dR / d : (a) probability densities for two beta distributions, (b) probability densities for two triangular distributions, ............................... 84

37. (continued) Candidate sampling distributions for cavity wastage rate, dR / d : (c) probability densities for two Weibull distributions, and (d) cumulative probabilities. ......................................... 85

38. Candidate sampling distributions for the Case 1 flaw initiation time , flaw-ini (1) :

(a) probability densities and (b) cumulative probabilities. (ASP analysis only)................................. 86

39. Candidate sampling distributions for the Case 2 flaw initiation time , flaw-ini (1) :

(a) probability densities and (b) cumulative probabilities. (ASP analysis only)................................. 87

40. Candidate sampling distributions for flaw growth rate , da / d : (a) probability densities and (b) cumulative probabilities. ............................................................................................................... 88

41. Event trees based on LOCA screening rules representing (a) best-estimate (BE), ............................. 97

41. (continued) Event trees based on LOCA screening rules representing: (b) more conservative (MC),................................................................................................................................................... 98

41. (continued) and (c) less conservative (LC) rules. ............................................................................... 99

42. Load paths (Japplied vs Applied Pressure) for three model surface flaws placed in the Davis-Besse finite-element submodel compared to the cladding materials ductile-tearing initiation fracture toughness with cumulative probabilities from the lognormal model shown in Fig. 19. ................... 102

43. A scaling rule is applied to the Driving-Force Model (Eq. (44) to bring the Monte Carlo model driving forces into agreement with the results of the detailed FEM wastage cavity submodel (a) scaling required for Model Flaw 1, (b) scaling required for Model Flaw 3, and (c) linear fit developed between Model Flaws 1 and 3. ........................................................................................ 103

44. Convergence of Case 001 as a function of the number of antithetic paired realizations at 1 year after time of discovery. ..................................................................................................................... 106

45. LOCA probability history for: (a) the best-estimate Cladding Capacity Analysis (CCA) case (Case CCA-001) with a further categorization into small-break LOCA (SBLOCA), medium-break LOCA (MBLOCA), and large-break LOCAs (LBLOCA) and (b) varying model flaw depth at TOD, 0.035 in., 0.065 in., and 0.0995 in............................................................................. 108 ix

LIST OF FIGURES (continued)

Figure Page

46. Sensitivity of CCA LOCA probabilities to: (a) model flaw depth at TOD (b) cavity growth/shape rules. ........................................................................................................................... 109

47. Sensitivity of CCA LOCA probabilities to LOCA binning rules (a) total LOCA probability, (b)

SBLOCA probability, (c) MBLOCA probability and (d) LBLOCA probability.............................. 109

48. Sensitivity of CCA LOCA probabilities to modeling of A-tip ductile tearing (a) total LOCA probability, (b) SBLOCA probability, (c) MBLOCA probability and (d) LBLOCA probability..... 110

49. Comparison of beta distributions for effective cavity radius at TOD-1 year, 0 , where the original distribution (based on the results of an Expert Elicitation) has been revised to produce a sampled R1(median) closer to the observed valued used in the CCA: (a) probability densities and (b) cumulative probabilities. ............................................................................................................. 112

50. Comparison of Weibull distributions for the time of flaw initiation (relative to TOD),

flaw-ini(2) (see the time line in Fig. 1), where the original distribution (based on the results of an Expert Elicitation) has been revised to produce a sampled R1(median) closer to the observed valued used in the CCA: (a) probability densities and (b) cumulative probabilities. ................................... 113

51. Best-estimate sampling distribution for cavity wastage rate compared to sampled values:

(a) probability density and (b) cumulative probabilities. .................................................................. 114

52. Best-estimate sampling distribution for flaw growth rate compared to sampled values:

(a) probability density and (b) cumulative probabilities. .................................................................. 115

53. Distribution of uncertain effective cavity radius calculated for the time of discovery (TOD): (a) frequency distribution and (b) cumulative probabilities. .................................................................. 116

54. Distribution of uncertain effective flaw depth calculated for the time of discovery (TOD):

(a) probability density and (b) cumulative probabilities ................................................................... 117

55. Comparison of LOCA probability histories between the CCA study and the ASP study conditions where the deviation is due to uncertainties in cavity wastage rate, cavity size at TOD-1, flaw growth rate, and flaw initiation time. For the CCA study, the damage state at TOD was treated as known with no uncertainty................................................................................................ 120

56. LOCA probability history for the best-estimate ASP case (Case ASP-001) with a further categorization into small-break LOCA (SBLOCA), medium-break LOCA (MBLOCA), and large-break LOCAs (LBLOCA). ...................................................................................................... 121

57. LOCA probability histories for full case matrix: (a) total LOCA probabilities, (b) SBLOCA probabilities,...................................................................................................................................... 122

57. (continued) LOCA probability histories for full case matrix: (c) MBLOCA probabilities, and (d) LBLOCA probabilities. ............................................................................................................... 123 x

1 Introduction 1.1 Method of Analysis In support of both the United States Nuclear Regulatory Commissions (NRCs) investigation into the structural integrity of the damaged Davis-Besse reactor pressure vessel (RPV) head and the NRCs Accident Sequence Precursor (ASP) Program,1 the Heavy-Section Steel Technology (HSST) Program at Oak Ridge National Laboratory (ORNL) has carried out a probabilistic structural mechanics (PSM) analysis of the degraded RPV head [1]. The objective of this report is to present a description of the Davis-Besse wastage-area damage model, the technical basis for the model, and the results of two PSM studies: (1) a cladding capacity analysis (CCA) of the wastage cavity as it existed at the time of discovery (TOD) and (2) an ASP study for a time period ranging from 1 year before the TOD up to the TOD. A companion report [2] describes a parallel experimental program performed by HSST at ORNL.

PSM is an analysis methodology that combines deterministic damage models with probabilistic representations of unknown or uncertain input parameters [3]. Random variations in the plastic-flow and ductile-tearing properties of the cladding, damage accumulation rates (in terms of cavity-wastage and flaw-growth rates), and flaw-initiation times for the damage model are combined to make an assessment of the reliability of the structure. The method of Monte Carlo Simulation [4] was selected to combine these uncertainties by randomly sampling from their prospective distributions to obtain the required inputs for multiple deterministic realizations of the damage state of the degraded RPV head. By repeating the damage-state realizations a large number of times, the resulting failure states (defined as the consequences of the damage states) constitute a random sample from the probability distribution over the output induced by the combined probability distributions sampled over the several input variables [4].

In both the CCA and ASP studies, a probability of failure is estimated from the fraction of the total number of simulations whose results predict some degree of damage that is judged to be equivalent to structural failure. In the CCA study, the As-Found damage state at TOD is assumed known; however, in the ASP study, estimated uncertainties in the As-Found condition of the wastage cavity are incorporated into the analysis.

1 The statutory ASP Program is under the functional responsibility of the NRCs Office of Nuclear Regulatory Research.

1

The consequences or degrees of failure of the damage state can be ranked by differentiating between small-break (SBLOCA), medium-break (MBLOCA), and large-break (LBLOCA) loss-of-coolant accidents. In this report, the damage states are investigated over a time interval that begins one year before the day of discovery (i.e., 16 February 2001) and extends beyond the day of discovery until failure of the cladding is predicted for all postulated conditions. These results address questions raised as part of the statutory requirements of the ASP analysis as well as the question concerning how long the Davis-Besse RPV could have continued in service (i.e., if the wastage cavity had not been discovered on 16 February 2002) before the failure of the vessel.

1.2 Background of the Davis-Besse Degraded RPV Head Problem Pursuant to the licensees commitments to NRC Bulletin 2001-01 [5], the Davis-Besse Nuclear Power Station2 began a refueling outage [6] on February 16, 2002 that included a 100% volumetric inspection of the sixty-nine RPV head penetrations. With emphasis placed on control-rod drive-mechanism (CRDM) nozzles, these inspections identified five nozzles with cracks that needed repair before restarting [7],

including axial indications in three CRDM nozzles (nozzles # 1, 2, and 3, located near the center of the RPV head) that were experiencing pressure-boundary leakage.

The repair procedure included roll-expanding the CRDM nozzle into the surrounding RPV head and then machining along the axis of the nozzle to a point above the indications. On March 6, 2002, the licensee prematurely terminated the machining process on nozzle # 3. During the removal of the machining apparatus, the nozzle was mechanically agitated and subsequently displaced (or tipped) in the downhill direction (away from the top of the RPV head) until its flange contacted the flange of the adjacent CRDM nozzle. To identify the cause of the displacement, the licensee investigated the condition of the RPV head surrounding nozzle # 3. This investigation included removing the CRDM nozzle from the RPV head, removing boric acid deposits from the top of the RPV head, and ultrasonically measuring the thickness of the RPV head in the vicinity of CRDM nozzles # 1, 2, and 3.

Upon completion of the removal of boric-acid crystal deposits on March 7, 2002, the licensee conducted a visual examination of the area and identified a large cavity in the RPV head on the downhill side of CRDM nozzle # 3 [1]. Follow-up characterization by ultrasonic testing (UT) indicated wastage of the low-alloy steel RPV head adjacent to the nozzle. The wastage region was found to extend approximately 2

The Davis-Besse Nuclear Power Station is located in Oak Harbor, Ohio. Operated by First Energy Nuclear Operating Company (FENOC), Davis-Besse achieved initial criticality on 12 August 1977, and it came online on 31 July 1978. Manufactured by Babcock and Wilcox (B&W), the raised-loop pressurized water reactor (PWR) has a licensed thermal-power output of 2772 MW(t) and a capacity of 873 net MW(e), where capacity is defined by the Energy Information Administration (EIA) as the net summer capability as reported in EIA survey form 860.

2

5 inches downhill on the RPV head from the penetration for CRDM nozzle # 3, with a width of approximately 4 to 5 inches at its widest part. In addition to damage described above, follow-up metallographic studies, carried out by BWXT [16] and discussed in Sect. 2.1.3, discovered a group of surface cracks located in the cladding on the side exposed to the corrosive environment.

The investigation of the causative conditions surrounding the degradation of the RPV head at Davis-Besse is continuing [7]. Boric acid or other contaminants could be contributing factors. Other factors contributing to the degradation might include the environment of the RPV head during both operating and shutdown conditions (e.g., wet/dry), the duration for which the RPV head was exposed to boric acid, and the source of the boric acid (e.g., leakage from the CRDM nozzle or from sources above the RPV head such as CRDM flanges).

See Fig. 1a for a photograph of the Davis-Besse RPV, Fig. 1b for a schematic of a typical nuclear power reactor showing the relative location of the head degradation, Fig. 1c for a cutaway sketch of the wastage cavity, and Fig. 1d for a photograph of the wastage area footprint.

1.3 Accident Sequence Precursor Program The NRC established the Accident Sequence Precursor (ASP) Program [8] in 1979 in response to the findings of the Risk Assessment Review Group as reported in NUREG/CR-0400 [9]. The ASP Program has the primary objective of systematically evaluating the operating experience of the U.S. nuclear power plant fleet to identify, document, and rank those operating events that were most significant in terms of the potential for inadequate core cooling and core damage (precursors). From SECY-99-289 [8], the ASP Program also has the following secondary objectives: (1) to categorize the precursors for plant-specific and generic implications, (2) to provide a measure that can be used to trend nuclear plant core damage risk, and (3) to provide a check on the dominant core damage scenarios predicted by probabilistic risk assessment (PRA) analyses.

3

(a)

(b)

Fig. 1. Depictions of Davis-Besse RPV and wastage cavity: (a) Davis-Besse Nuclear Power Station RPV [7]; (b) schematic of a typical nuclear power reactor showing the relationship of the CRDM nozzles to the RPV head and the location of the damaged area [7];

4

(c)

(d)

Fig. 1. (continued) (c) sketch of RPV head degradation; and (d) photograph presenting a plan view of the wastage-area footprint showing remaining J-groove weld, irregular topography of exposed unbacked cladding, and irregular sidewalls of the wastage cavity [7].

5

An ASP is an operational event or a plant condition that is an important element of a postulated accident sequence3 associated with inadequate core cooling, a sequence that would be expected to result in core damage [10]. The ASP Program identifies nuclear power plant events that are considered precursors to accidents with the potential for severe core damage and uses risk assessment methodologies to estimate the significance of the events. The figure of merit for ASP analyses is the conditional core damage probability (CCDP), where the CCDP is the conditional probability of core damage given the failures observed in the event. Events with CCDPs greater than 1.0 x 106 are judged to be accident sequence precursors.

1.4 Previous ORNL Studies of the Davis-Besse RPV Head Degradation Problem Previous ORNL reports have described (1) the development of a stochastic model for the estimation of the probability of failure of the wastage cavity by plastic collapse under internal pressure loading [11,12]

and (2) the technical basis for applying a surrogate damage model (specifically, a circular burst disk) [13]

to simulate the initial damage state, damage accumulation, and failure consequences for the irregularly shaped wastage region in the Davis-Besse RPV head.

1.4.1 Detailed Finite-Element Global and Submodels of Wastage Area and Cavity Detailed finite-element models were constructed and applied in the studies discussed in [11-13]. The results of these analyses were used in the development of a more detailed 3D finite-element wastage-area damage model to be described in Sect. 2. The results of the detailed 3D FEM model in turn provided guidance in the construction of a simplified model for use in Monte Carlo simulations. The following discussion provides a summary description of those prior FEM analyses.

The submodeling capabilities of the ABAQUS finite-element code [14] were employed in the plastic-collapse stochastic model development described in [11, 12] and in the cavity growth studies described in

[13] to focus the available computational resources on the region around the wastage area cavity at CRDM nozzle # 3. Submodeling can be used to investigate a portion of a model with a refined mesh, where the boundary conditions of the submodel are driven by an interpolation of the displacement solution obtained from a relatively coarse global model. The technique is appropriate and accurate when it is necessary to obtain a refined, detailed solution in a local region, and the detailed modeling of that local region has a negligible effect on the global solution. Communication between models proceeds in one direction only, from the global model to the submodel. As shown in Figs. 2a and 3, the initial global 3

An accident sequence is a combination of events leading from an initiating event that challenges safety systems to an undesired consequence. The sequence is ordered, starting with the initiating event, and proceeds through sequential failures leading to the final consequence [3].

6

model consisted of the full RPV head (with all 69 penetrations) and closure flange. No cladding or CRDM nozzles were included in the global model used in the studies described in refs. [11-13].

The initial submodel (see Fig. 2b) consisted of the cladding (SS308), base (A533B), and CRDM nozzles (A600) #3, #11, #15, and #16. The plan views of the RPV head in Fig. 4a and 4b indicate the position and geometry of the submodel with respect to the global model. Figure 5 shows the ProEngineer solid model of the wastage cavity submodel. This solid model was imported into MSC Patran where the finite-element mesh was constructed.

1.4.2 Results of Wastage-Cavity Growth-Pattern Study [13]

The results of a study of the sensitivity of burst-pressure predictions to the modeled shape of the wastage-cavity footprint were reported in [13]. As shown in Fig. 6, three growth patterns were examined: (1) self-similar (relative to the as-found cavity geometry), (2) ellipsoidal, and (3) circular (i.e., a circular burst disk). Detailed finite-element analyses were carried out, using the global/submodel approach described in Sect. 1.4.1, for the postulated self-similar and ellipsoidal growth patterns, and a plastic-instability theory for circular diaphragms under lateral-pressure loading developed by Chakrabarty and Alexander [15] was applied for the circular growth pattern. The results of that growth-pattern study are summarized in Fig. 7, where the 5th , 50th, and 95th % percentile curves for the circular growth pattern model are based on a scaling of the stochastic model for uncertainties in plastic-collapse burst pressure predictions developed in

[11,12]. Drawn from the results summarized in Fig. 7, a fundamental conclusion is that the unbacked area of the cladding is the dominant geometric parameter for burst-pressure predictions of a model of the Davis-Besse wastage cavity when the unbacked stainless steel cladding does not contain any crack-like flaws; the specific shape of the footprint represents a second-order effect. This conclusion allowed the development of the wastage-area damage-state failure model to be based on failure predictions for a circular diaphragm under lateral-pressure loading. Subsequent discovery of surface flaws in the wastage-cavity cladding added another failure mechanism, namely the potential for cladding failure by ductile tearing.

7

(a) (b)

Fig. 2. Finite-element global model and submodels applied in previous analyses [11,13] of the Davis-Besse head and wastage area. The displacements at the vertical side boundaries of the submodel are driven by the global model. Both models were exposed to the same internal pressure loading.

Fig. 3. Geometry of the hemispherical RPV head and closure flange used in the global model

[11,13] (B&W proprietary dimensions have been blacked out).

8

Fig. 4. Geometry of submodel relative to nozzles #3, #11, #15, and #16 in the RPV head.

Fig. 5. Wastage-cavity submodel [11,13] developed from a ProEngineer solid model and then meshed with MSC Patran.

9

(a)

(b)

Fig. 6. A computational study was carried out in [13] to determine the sensitivity of burst pressures, predicted by finite-element analyses, to the geometry of estimated growth patterns of the wastage cavity footprint: (a) self-similar and ellipsoidal growth patterns and (b) circular burst disk.

10

Fig. 7. The FEM-predicted burst pressures of irregularly-shaped footprint growth patterns are compared in [13] to burst pressure predictions for circular diaphragms under lateral pressure loading based on an instability theory proposed by Chakrabarty and Alexander

[15].

11

1.5 Scope of this Report Section 2 presents the development of the Davis-Besse wastage-area damage model, including both a detailed 3-dimensional FEM model and a simplified model for Monte Carlo simulations. Section 3 presents the results of the deterministic simulations using the FEM model as well as best-estimate and sensitivity Monte Carlo simulations in support of both the CCA and the ASP studies. Section 4 provides a summary and conclusions. Appendix A provides a description of the of the Davis-Besse Monte Carlo damage state code along with a flowchart and guidance for its use. Appendix B presents a detailed listing of the CCA and ASP Monte Carlo results in terms of tables of LOCA probabilities as a function of time starting at 1 year before TOD up to TOD. Appendix C presents the final report on the examination of the reactor pressure vessel head degradation at Davis-Besse carried out by BWXT Technologies, Inc. at Lynchburg, VA.

12

2 Davis-Besse Wastage-Area Damage Model 2.1 Detailed 3D FEM Model for Deterministic Analyses A detailed 3-dimensional FEM submodel of the Davis-Besse wastage cavity was developed from the results [16] of visual inspections, dye-penetrant testing, scanning electron microscopy, energy dispersive spectroscopy, metallography, Knoop microhardness testing, and wastage cavity geometric measurements carried out by the Lynchburg Technology Center for Metallurgical Examinations, BWXT Technologies, Inc., and reported in [16]. The submodeling approach described in [11-13] was used for the current study, where a detailed submodel is driven both by direct pressure loading and by the displacement solution obtained from a coarser global model (see Fig. 8).

The wastage cavity submodel shown in Fig. 9 consisted of four CRDM nozzles (specifically #3, #11, #15, and #16)4, the base material with wastage cavity, the J-groove weld for nozzle #3, the backed cladding, the exposed unbacked cladding, and a modeled cladding surface flaw on the wastage cavity side of the unbacked cladding.

2.1.1 Modeled Wastage Cavity Footprint Details of the wastage cavity footprint were obtained from photographs available in [16] (Fig. 9c) (see Appendix C) and then digitized to a locus of 66 points which were used to construct a spline curve that established the basis for the footprint. The irregular topography of the unbacked cladding region exposed to the corrosive effects of the boric acid was also simulated by a sinusoidal waveform applied to the exposed upper surface.

2.1.2 Wastage Cavity During examinations at BWXT Lynchburg, a dental mold [16] of the wastage cavity was prepared to establish a permanent record of the wastage-cavity shape. The photographs available in [16] were applied in the construction of a solid model (see Fig. 9d) using the Solidworks [17] modeling software.

4 For CRDM nozzles #3, #11, #15, and #16, the corresponding core grid locations are G9, F10, E9, and G11, respectively. Nozzles #15 and #16 are unused spare nozzles [1].

13

(a)

(b)

Fig. 8. Davis-Besse global model of RPV closure head for the current study includes a simplified wastage-cavity definition and inner cladding layer.

14

(a)

(b)

Fig. 9. Davis-Besse submodel includes: (a) nozzles #3, #11, #15, #16 (Alloy 600); (b) a refined defi-nition of the wastage-cavity morphology in the base material (A 533B pressure vessel steel),

15

(c)

(d)

Fig. 9. (continued) (c) wastage cavity footprint; (d) irregular geometry of wastage cavity walls based on geometric measurements from the dental mold taken by BWXT [16].

16

(e)

(f)

Fig. 9. (continued) (e) irregular exposed cladding topography based on a sinusoidal waveform, and (f) J-groove weld near a modeled surface flaw.

17

This solid model of the dental mold was then applied in a loft cut operation to remove material from the base section, thus forming the degraded wastage cavity within a Solidworks solid model of the damaged RPV head. The solid model of the RPV heads base material with wastage cavity was then imported into MSC Patran to create the finite-element mesh of the wastage cavity submodel.

2.1.3 Modeled Surface Flaws Extensive metallographic analyses were carried out by BWXT [16] to better define the nature and depth of the cladding flaws. These analyses identified two key features of the cracks observed in the cladding

[18].

(1) The cracking was more prevalent and the cracks were much deeper than was originally thought, based on the initial preliminary investigation reported in the Summer of 2003

[16]. Cracking existed on the cladding interface, as evident in the initial findings. Each cross section examined near the center of the exposed cavity region revealed several distinct cracks.

(2) The crack morphology was intergranular and/or interdendritic in nature. This morphology is characteristic of intergranular stress corrosion cracking (IGSCC) prevalent in relatively high-carbon stainless steels used in boiling-water-reactor (BWR) service environments.

The following summary of the BWXT findings is quoted from [16]:

The exposed Type 308 stainless steel cladding surface area was approximately 16.5 square inches. The average cladding thickness measured by dial calipers was 0.256", with a minimum thickness of 0.202" and maximum thickness of 0.314". Thickness measurements taken on a transverse metallographic mount prepared through the thinnest portion of exposed cladding indicated a minimum cladding thickness of 0.179". This area, which was located in an uncracked region of the cladding adjacent to the J-groove weld, was not accessible during the dial caliper measurements. The cladding mechanical properties and chemical composition appeared uniform across the cladding thickness. The underside (reactor coolant system or RCS side) of the cladding and the unexposed portions of cladding did not contain any cracks or other signs of deterioration.

Shallow intergranular attack (IGA) was observed on all exposed cladding surfaces examined.

Deeper cracks, which extended a maximum 0.057" below the exposed cladding surface, followed a mixed interdendritic/intergranular path. These cracks initiated in the IGA and propagated under a stress-corrosion cracking mechanism along the ferrite stringers under the influence of an applied tensile stress (i.e., system pressure). The minimum observed distance from an exposed cladding crack tip to the RCS side of the cladding was 0.188".

Based on the results of these examinations, three modeled surface flaws (see Fig. 9f) were developed as representative of a best-estimate flaw and two bounding flaws (see Fig. 10 for a generic description of the model flaw geometry).

Flaw Model #1 - aA-tip = 0.065 in., L = 0.66 in.; (best-estimate) characteristic of deep flaws observed in the Davis-Besse wastage cavity cladding Flaw Model #2 - aA-tip = 0.0995 in., L = 0.606 in.; characteristic of the deepest flaw observed in the Davis-Besse wastage cavity cladding Flaw Model #3 - aA-tip = 0.0995 in., L = 2.0 in.; characteristic of the ASME Boiler and Pressure Vessel 18

Code representation of cracks,5 using the proximity rules and conservative flaw characterization procedures established in Sect. XI, Division 1, Subsection IWA-3300

[19], Division 2, Subsection IWB-3600 [20], and Appendix A [21].

2.2 Simplified Wastage-Area Damage Model for Monte Carlo Simulations The damage states of the Davis-Besse degraded RPV head at different times (on the time-line history of the head) are estimated from the results of PSM analyses of a simplified damage model consisting of a circular diaphragm (or burst disk) under lateral-pressure loading (see Fig. 10). The multiple defects, observed in the Davis-Besse wastage cavity after extraction from the RPV head, are simulated in the model by a single long surface flaw ( L = 2 in.) centered in the unloaded side of the diaphragm with an uncertain (i.e., sampled) flaw growth rate. The assumed transferability of the results of the PSM analysis from the circular-diaphragm damage model back to the actual wastage area in the Davis-Besse RPV head is inferred from the results of the study described in [13] and summarized in Sect. 1.4.2.

For failure by plastic collapse of the cavity and/or the remaining ligament below the modeled flaw, similitude between the calculated structural reliability of the damage model and the Davis-Besse wastage cavity is established by matching the effective unbacked area under load and the depth of the model flaw.

Additional scaling (to be discussed) is required to establish similitude of the elastic-plastic crack driving forces (i.e., J-integrals) between the flaw in the burst-disk damage model and the modeled flaw in the detailed Davis-Besse wastage cavity submodel.

5 It is recognized that the ASME Code does not explicitly treat cracks in the cladding in terms of a structural challenge to the RPV; however, the spirit and intent of the Code is applied herein for this special case of unbacked cladding being required to maintain the integrity of the pressure boundary.

19

Fig. 10. Geometry, loading, and constraint conditions for the Davis-Besse RPV head damage model

- a 1/4 symmetry burst disk with a centered surface flaw, where h0 = 0.25 in. and L = 2.0 in.

20

2.2.1 Damage Model Geometry, Loading, and Constraint Conditions Figure 10 presents a 1/4-symmetry representation of the damage model used in this study in which a circular diaphragm with radius, Rcavity, and undeformed thickness, h0, is loaded on one side by a uniform applied pressure, Papplied. A flat-bottomed surface flaw of length, L, is centered in the diaphragm on its outer, unloaded surface. The depth of the crack tip, aA-tip, is constant along the length (L-2 aA-tip) of the flat bottom of the flaw. This bottom section of the crack tip is designated as the A-tip. The two C-tips of the flaw are circular arcs tracing the path of the crack tip from the A-tip up to the unloaded free surface of the diaphragm. The radius, Rcavity, is scaled so that the diaphragms loaded area matches the as-found or postulated unbacked surface area, Awastage, of the irregularly shaped wastage-cavity footprint, and the as-found model flaw depth, aA-tip, is equal to either a flaw depth representative of the observed range of flaws (0.036 to 0.099 inches) in the wastage cavity (as-found condition) or an uncertain flaw depth sampled from statistical distributions (to be described):

Awastage fixed at time of discovery for CCA studies Rcavity = +  ; a Atip =

or sampled for ASP analyses (1)

The outer edge of the diaphragm is constrained by an encastre boundary condition.

2.2.2 Damage-State Parameters The two uncertain geometric parameters in Eq. (1) (see Figs. 10 and 11) are the effective cavity radius, Rcavity, and the initial flaw depth of the A-tip of the centered modeled flaw, aA-tip. Together, (Rcavity, aA-tip) define a damage state point for the wastage-cavity damage model. Given this damage state point and a prescribed applied pressure, Papplied, the failure state of the model can then be determined based on the position of the state point in damage space relative to the models fragility curve, developed from the instability theory due to Chakrabarty and Alexander [15] (see Sect. 2.2.3). This fragility curve is based on an assumed failure mechanism of plastic collapse of the remaining ligament (termed an A-tip collapse) below the centered flaw in the damage model. Additional potential A-tip failure modes consisting of either ductile-tearing initiation or tearing initiation followed by unstable ductile tearing of the A-tip are also included in the analysis. Upon the failure of the A-tip, failure mechanisms for the C-tips and remaining wastage cavity material are then investigated to determine the final failure state (expressed in terms of LOCA size) of the cavity at the time of discovery. Along the locus of damage-state points defined in damage space by the fragility curve, A-tip failure is predicted for the given applied pressure; therefore, the fragility curve partitions the damage space of the wastage cavity into stable and unstable (failed) regions.

21

Fig. 11. Spherical geometry of deformation assumed in Chakrabarty and Alexanders [15] plastic instability theory. Chakrabarty and Alexander assumed that the burst disk was unflawed, i.e., aA-tip = 0.

22

2.2.3 Failure Mechanisms - Development of a Plastic-Collapse Fragility Curve The fragility curve applied in the PSM analysis for the Davis-Besse problem for an unbacked cladding model with a centered flaw has the following functional form:

a A tip ( crit ) = h0 + ( S x Rcavity ) ; 0 Rcavity Rintercept where a A tip ( crit ) = the critical flaw depth at which plastic collapse occurs h0 = initial undeformed thickness of the cladding (2)

S = fragility curve slope Rcavity = effective radius of the wastage cavity h0 Rintercept =  ; for a specified applied pressure, Papplied S

The slope, S (a function of the specified applied pressure, Papplied, and cladding plastic flow properties), of the fragility curve can be developed from an instability criterion due to Chakrabarty and Alexander [15]

with a modification for this analysis to treat the remaining ligament below the flaw.

When the damage state of the wastage-area intersects the fragility curve, the remaining ligament of the A-tip is predicted to fail by plastic collapse. The value of Rintercept places an upper bound on the size of the wastage cavity (irrespective of the presence of the model flaw) for the specified applied pressure. The flaws A-tip may approach the fragility curve by flaw growth due to exposure to the corrosive environment of the wastage cavity, by flaw growth due to stable ductile tearing, or by flaw growth due to unstable ductile tearing. When the damage state of the wastage-area model intersects the fragility curve due to combined growth of the flaw and cavity, the remaining ligament below the flaw is predicted to fail by plastic collapse, and the bottom of the model flaw blows out of the cavity, leaving a 2-inch long rip in the cladding. Additional failure sequences may then follow toward complete failure of the wastage cavity.

The burst pressure for a circular diaphragm under lateral pressure loading can be estimated using the Chakrabarty and Alexander (1970) [15] instability criterion as discussed in [11]. In the PSM damage model, the failure of an unflawed wastage cavity is simulated by failure due to plastic instability of a pressure-loaded circular diaphragm (i.e., burst disk) with a fixed-grip constraint (encastre boundary condition) along the edges. The radius of the diaphragm matches the effective radius, Rcavity, of the wastage cavity based on a specified unbacked wastage area (see Eq. (1)).

23

For a given material defined by its plastic-flow properties (K, n)6 and unflawed geometry (Rcavity, aA-tip= 0),

the failure pressure (i.e., the burst pressure) can be estimated by the following procedure (see Fig. 11 for a description of the undeformed and deformed geometries):7 2(2 n)(1 + 2n)

• Calculate the effective critical radial true strain at = 0° . crit =

11 4n n

• Calculate the corresponding effective critical true stress. crit ( K , n) = K crit

• Calculate the critical radial deformed thickness at = 0° . (

hcrit = h0 exp crit )

crit

• Calculate the polar height at the critical condition. H crit = Rcavity exp 1 2

2 H crit + Rcavity 2

• Calculate the corresponding bulge curvature radius. Rcrit =

2 H crit 2hcrit crit

• Finally, calculate the predicted burst pressure. pburst =

Rcrit For a surface flaw centered in a circular diaphragm, this procedure can be modified by replacing the initial undeformed diaphragm thickness, h0, with the undeformed remaining ligament thickness below the flaw (see Fig. 11) such that

(

hcrit = ( h0 a Atip ( crit ) ) exp crit ) (3) where aA-tip(crit) is the critical flaw depth at the deepest point of the flaw.

The above procedure can be inverted to solve for the critical flaw depth at failure, aA-tip(crit), as a function of the claddings plastic flow properties, (K, n), the initial cladding thickness, h0, applied pressure, Papplied, and cavity radius, Rcavity, to obtain crit ( n )

Papplied exp a A tip ( crit ) ( K , n, h0 , Papplied , Rcavity ) =

2 x Rcavity + h0 (4) crit ( n )

4 crit ( K , n ) exp crit ( n ) exp 2 1 where n

6 A power-law constitutive model is assumed for plastic deformation, where ( K , n) = K .

7 Following the notational conventions of Dowling [34], and are true stress and true strain, respectively, and and are engineering stress and engineering strain, respectively.

24

2(2 n)(1 + 2n) crit (n) =

11 4n n

crit ( K , n) = K crit (n)

The parameter, S, in Eq. (2) is, therefore, from Eq. (4) for a given applied pressure, Papplied :

crit ( n )

Papplied exp S ( K , n l Papplied ) = 2 (5) crit ( n )

4 crit ( K , n ) exp[ crit ( n )] exp 1 2

If we assume that the initial cladding thickness, h0, is known, then for a specified applied pressure the dominant sources of uncertainty in the construction of the fragility curve are the plastic flow properties of the cladding, (K, n). Note that uncertainties in the initial cladding thickness can be treated independently from the uncertainties in the slope of the fragility curve.

2.2.3.1 Stochastic Model for Plastic-Flow Properties - Sampling for yield and ult In Eq. (5), the dominant sources of uncertainty in the slope of the fragility curve are the plastic flow properties of the cladding as characterized by the coefficient, K, and exponent, n, of the power-law representation of the true effective stress vs. true strain curve for a given temperature. Table 1 presents a database of plastic-flow properties for 550 °F temperature 600 °F used in this study, and Fig. 12 compares the stress vs. strain curves developed from power-law fits of the data in Table 1 to the best-estimate constitutive properties for the Davis-Besse cladding material. It was decided to focus the stochastic modeling task on characterizing the uncertainties of the yield and ultimate stresses which, for a fixed temperature, are presumably due either to material variability between heats or spatial variability within a single heat. It will be shown in Sect. 2.2.3.2 how the parameters K and n can be calculated from sampled values of yield and ult .

As discussed in refs. [22, 23], independent samples of yield and ultimate stress indicate that a statistical correlation exists between the two properties. In refs. [22, 23], a bivariate lognormal distribution was selected as the joint distribution since the fitted marginal distributions for their dataset (TP304 stainless steel pipe material) followed lognormal distributions reasonably well. For this study, several marginal distributions were fitted to the set of yield and ultimate true stress data in Table 1 (see Table 2 for the Goodness of Fit statistics of Anderson-Darling [24], Kolmogorov-Smirnov [25], and equi-probable Chi-square [26]; and also see Figs. 13a and 13b) using the ExpertFit© [27] statistical software. Three marginal distributions were selected for further study, and they have the following forms:

25

Table 1. Plastic-Flow Properties for SS Cladding Specimen Test True Stresses Source Material ID Temp. Yield Ultimate K* n*

(°F) (ksi) (ksi) (ksi) (-)

NUREG/CR-3927 CPC-79 550.0 29.79 73.24 102.09 0.218 NUREG/CR-3927 CPC-80 550.0 29.79 74.40 103.86 0.221 NUREG/CR-5511 3-wire Cladding A20A 550.0 31.70 72.96 100.94 0.206 NUREG/CR-5511 3-wire Cladding A20B 550.0 36.60 69.97 93.18 0.168 NUREG/CR-6363 3-wire Cladding A18C 550.0 28.60 75.37 103.56 0.227 NUREG/CR-6363 3-wire Cladding A24A 550.0 31.00 71.54 99.00 0.206 HSST Clad Burst Tests PVRUF Cladding 1-1 600.0 29.81 68.94 95.45 0.206 HSST Clad Burst Tests PVRUF Cladding 1-2 600.0 28.85 68.05 94.41 0.209 HSST Clad Burst Tests PVRUF Cladding 1-3 600.0 28.14 66.44 92.16 0.209 HSST Clad Burst Tests PVRUF Cladding 1-4 600.0 29.30 66.44 92.16 0.199 HSST Clad Burst Tests PVRUF Cladding 2-1 600.0 29.82 70.03 97.14 0.209 HSST Clad Burst Tests PVRUF Cladding 2-2 600.0 29.93 69.56 96.37 0.207 HSST Clad Burst Tests PVRUF Cladding 2-3 600.0 30.97 68.24 93.98 0.197 HSST Clad Burst Tests PVRUF Cladding 2-4 600.0 31.22 67.59 92.88 0.194 BWXT DB Cladding B2C2A1 600.0 30.53 70.53 97.66 0.206 BWXT DB Cladding B2C2A2 600.0 31.36 73.49 101.96 0.209 NSE Handbook SS 308/308L Table II 600.0 31.00 69.65 114.91 0.228

()

n

• eff = K where eff = true effective stress and = true total strain See Appendix C for a discussion of the BWXT [16] investigations.

Table 2. Goodness of Fit Statistics for Yield and Ultimate Stress Statistical Distributions Sampled Marginal Point Estimates for Parameters Goodness of Test Critical Values for Level of Significance ()

Variate Distribution location scale shape Fit Test Statistic 0.25 0.15 0.1 0.05 0.025 0.01 0.005 Anderson-Darling 0.406 0.452 NA 0.608 0.722 0.836 0.99 NA Yield Stress inverted Weibull 0 29.7279 24.7640 Kolmogorov-Smirnov 0.625 NA NA 0.779 0.843 0.907 0.973 NA equi-probable chi 2 2.706 5.385 6.745 7.779 9.488 NA 13.277 NA Anderson-Darling 0.862 0.447 NA 0.600 0.715 0.830 0.984 1.102 Yield Stress lognormal 0 3.4159 0.0567 Kolmogorov-Smirnov 0.754 NA 0.740 0.782 0.854 0.950 0.988 NA equi-probable chi 2 7.412 5.385 6.745 7.779 9.488 NA 13.277 NA Anderson-Darling 0.768 1.248 NA 1.933 2.492 3.070 3.857 4.500 Yield Stress gamma 28.1119 1.6878 1.4118 Kolmogorov-Smirnov 0.922 NA 1.099 1.182 1.311 1.429 1.572 NA equi-probable chi 2 8.000 5.385 6.745 7.779 9.488 NA 13.277 NA Anderson-Darling 0.255 0.452 NA 0.608 0.722 0.836 0.99 NA Ultimate Stress inverted Weibull 0 69.0596 30.7746 Kolmogorov-Smirnov 0.518 NA NA 0.779 0.843 0.907 0.973 NA equi-probable chi 2 3.294 5.385 6.745 7.779 9.488 NA 13.277 NA Anderson-Darling 0.311 0.447 NA 0.600 0.715 0.830 0.984 1.102 Ultimate Stress lognormal 0 4.2532 0.0375 Kolmogorov-Smirnov 0.551 NA 0.740 0.782 0.854 0.950 0.988 NA equi-probable chi 2 2.706 5.385 6.745 7.779 9.488 NA 13.277 NA Anderson-Darling 0.560 1.248 NA 1.933 2.492 3.070 3.857 4.500 Ultimate Stress gamma 66.241 2.830 1.462 Kolmogorov-Smirnov 0.649 NA 1.099 1.182 1.311 1.429 1.572 NA equi-probable chi 2 2.706 5.385 6.745 7.779 9.488 NA 13.277 NA

• NA = critical point not available 26

Fig. 12. Plots of true stress vs. true strain from the data in Table 1.

27

(a)

(b)

Fig. 13. Probability density functions for fitted statistical marginal distributions based on data in Table 1: (a) yield true stress and (b) ultimate true stress.

28

INVERTED WEIBULL MARGINAL DISTRIBUTIONS Inverted Weibull Distribution for Yield True Stress - W -1 (a,b,c) with ayield = 0, byield = 29.7279, and cyield = 24.764 Inverted Weibull Distribution for Ultimate True Stress - W -1 (a,b,c) with ault = 0, bult = 69.056, and cult = 30.7746 where a = location parameter, b = scale parameter, c = shape parameter, and the inverted Weibull probability density and cumulative distribution functions are defined below.

x a c c bc ( x a )(c +1) exp if x > a density: f ( x) = b 0 otherwise x a c exp if x > a CDF: Pr( X x) = F ( x l a, b, c ) = b 0 otherwise Sampling from a 3-parameter Inverted Weibull Distribution: X i W 1 (a, b, c)

A random number is drawn from a uniform distribution on the open interval (0,1) and then transformed to an inverted Weibull variate with the inverted Weibull percentile function.

U i U (0,1) b X (i ) = a +

[ ln(U i )]1/ c 29

LOGNORMAL MARGINAL DISTRIBUTIONS Lognormal Distribution for Yield True Stress - LN ( , µlog , log ) with

= 0, µlog( yield ) = 3.4159, log( yield ) = 0.0567 Lognormal Distribution for Ultimate True Stress - LN ( , µlog , log ) with

= 0, µlog(ult ) = 4.2532, log(ult ) = 0.0375 where µlog = lognormal mean (scale parameter), log = lognormal standard deviation (shape parameter),

and is the location parameter . The lognormal probability density and cumulative distribution functions are defined below.

ln( x ) µ 2 log exp 1

if x >

density: f ( x) = ( x ) 2 2 2 2 log log 0 otherwise ln( x ) µlog if x >

z 1 2 CDF: F ( x) = log where ( z ) = 2 2 d exp 0 otherwise Sampling from a 3-parameter Lognormal Distribution: X i LN ( , µlog , log)

The log-transformed deviate is sampled from a normal distribution with mean equal to the lognormal mean, µ log , and standard deviation equal to the lognormal standard deviation, log . A standard normal unit variate is first sampled. The following rational function [28] represents an accurate approximation of the standard normal percentile function:

1 p for p <

x= 2 1 p for p 1 2

y = 2ln( x) (6) 1 a + a y + a2 y 2 + a3 y 3 + a4 y 4 Z p = sgn p y + 0 1 2 b0 + b1 y + b2 y 2 + b3 y 3 + b4 y 4 30

where 1 if x < 0 sgn( x) =

+1 if x 0 and the coefficients of the rational function are:

a0 = -0.3222324310880000 b0 = 0.0993484626060 a1 = -1.0000000000000000 b1 = 0.5885815704950 a2 = -0.3422420885470000 b2 = 0.5311034623660 a3 = -0.0204231210245000 b3 = 0.1035377528500 a4 = -0.0000453642210148 b4 = 0.0038560700634 The standard normal deviate is scaled to obtain the required quantile, and the log-transformed deviate is then converted into the required random deviate by the exponential function.

U i U (0,1) 1 2ln(U i ) for U i <

2 y=

2ln(1 U ) 1 for U i i

2 1 a + a y + a2 y 2 + a3 y 3 + a4 y 4 Zi = sgn U i y + 0 1 2 + + 2 + b y3 + b y 4 b0 1 b y b2 y 3 4 Yi = µlog + Zi log X i = + exp(Yi )

GAMMA MARGINAL DISTRIBUTIONS Gamma Distribution for Yield True Stress - gamma( , , ) with

= 1.411176, = 1.6878, = 28.112 Gamma Distribution for Ultimate True Stress - gamma( , , ) with

= 1.4623, = 2.83, = 66.24053 where is the shape parameter, is the scale parameter, and is the location parameter . The gamma probability density and cumulative distribution functions are defined below.

31

( x ) 1 ( x )

exp if x >

density: f ( x ) = ( )

a where ( z ) = t z 1 exp ( t ) dt 0

0 otherwise 1 exp ( x ) (

1 x ) / j if x >

CDF: F ( x) = j =0 j ! if is an integer 0 otherwise If is not an integer, there is no closed form for F(x).

Sampling from a 3-parameter Gamma Distribution: X i gamma( , , )

The gamma distributed deviate is sampled from a standard gamma distribution, gamma ( ,1,0) , and then converted into the required random deviate by the following scaling:

Yi gamma( ,1,0)

X i = + Yi The inverse sampling Yi gamma( ,1,0) is accomplished with the FORTRAN subroutine CDFGAM available in the numerical library dcdflib obtained from the public-domain software repository netlib maintained at ORNL (see http://www.netlib.org). This subroutine uses the numerical procedures and coding specified in ref. [29].

As shown in Fig. 14, the above marginal distributions were then shifted to force their medians to pass through the current best-estimate values for the Davis-Besse cladding material, specifically yield ( med ) = 30.9 ksi and ult ( med ) = 72 ksi , derived from the average of the two sample points obtained from BWXT in Table 1.

{ }

T For the sampled vector X k = yield , ult with k = 1,2, ,n (where n = sample size), the 2 x 2 symmetric covariance matrix can be estimated by n

11 12 ( X ik X i )( X jk X j )

(2x2) = k =1 where ij = (7) 21 22 n and the correlation coefficient, 12 , is then estimated by the maximum-likelihood estimator 12 12 = = 21 (8) 11 22 32

(a)

(b)

Fig. 14. Statistical distributions approximating the material variability in the yield and ultimate stresses of cladding at 550 °F to 600 °F: (a) yield true stress CDFs with adjustments to match Davis-Besse best-estimate value of 30.9 ksi and (b) ultimate true stress CDFs with adjustments to match Davis-Besse best-estimate value of 72 ksi.

33

From the data in Table 1, a maximum likelihood estimate of the covariance matrix is 3.51519 0.35491 (2x2) = and 12 = 0.06931 0.35491 7.45912 and an estimate of the covariance matrix for the log-transformed data is 0.0034112 0.0001978 (2x2)(log) = and 12(log) = 0.08765 0.0001978 0.0014929 The statistically-inferred level of correlation between the true yield and true ultimate stresses for cladding is consistent with the data (after conversion to true stresses) presented in [22,23], where TP304 Stainless Steel Piping 4.4013 1.5280 (2x2) = and 12 = 0.07927 1.5280 84.4235 0.0085228 0.0008016 (2x2)(log) = and 12(log) = 0.08293 0.0008016 0.0109623 Three bivariate distributions along with the uncorrelated case were chosen to test the sensitivity of the Monte Carlo LOCA probability estimates to the degree and method of correlation when sampling the yield and ultimate stresses. (Compare the results of Case CCA-001 with Cases CCA-019, CCA-020, and CCA-021 in Table 10 of Sec. 3.2.)

(1) Bivariate Inverted Weibull Distribution - 12 = 1 (perfect correlation)

A single random number is drawn from a uniform distribution on the open interval (0,1) and then transformed to inverted Weibull variates for both yield and ultimate true stresses:

U i U (0,1) byield yield (i ) = a yield + (9)

[ ln(U i )] 1/ c yield bult ult (i ) = ault +

[ ln(U i )]1/ cult (2) Bivariate Lognormal Distribution - 12 = 0.06931; 12(log) = 0.08765 (inferred from data)

The sampling of a multivariate lognormal random vector follows from the procedures discussed in

[30, 31]. The first step is to carry out a Cholesky decomposition (factorization) of the log-transformed correlation matrix. For the log-transformed data in Table 1:

34

2x2(log) = CC T where C is the 2 x 2 lower triangular Cholesky factor c c 0.05666 0 C = 11 12 =

c 21 c 22 0.013151 0.035097

{ }

T For a bivariate lognormal joint distribution, the sampling protocol for the vector X i = yield (i ) , ult (i )

is as follows:

Z1 N (0,1); Z 2 N (0,1) where Z1 and Z 2 are IID Y1 = µ yield (log) + ( c 11x Z1 )

Y2 = µult (log) + ( c21 x Z1 ) + ( c22 x Z 2 ) (10) exp(Y1 ) yield Xi = =

exp(Y2 ) (i ) ult (i )

where IID refers to independent and identically distributed random variates, and Z1 and Z2 are standard-unit-normal variates.

(3) Bivariate Gamma Distribution - 12 = 0.06931 (inferred from data)

{ }

T The sampling protocol of the bivariate gamma random vector X i = yield (i ) , ult (i ) follows from the procedures discussed in ref. [32] with the requirement that the two variates must be positively correlated and min {1 , 2 } / 1 2 . These requirements are met by the yield and ultimate stress sample dataset in Table 1. The protocol relies on a special property of gamma distributions.

(

Y1 gamma 1 12 1 2 ,1,0 )

Y2 gamma ( 2 12 1 2 ,1,0 ) independent of Y1 Y3 gamma ( 12 1 2 ,1,0 ) independent of Y1 and Y2 (11)

+ (Y + Y ) yield Xi = 1 1 1 3 =

2 + 2 (Y2 + Y3 ) (i ) ult (i )

The above technique is known as trivariate reduction, since the three random variates Y1, Y2, and Y3 are

{ }

T reduced to the two final random variates yield (i ) , ult (i ) , see [33].

(4) Univariate Inverted Weibull Distribution - 12 = 0 (uncorrelated)

Two IID random numbers are drawn from a uniform distribution on the open interval (0,1) and then transformed to inverted Weibull variates for both yield and ultimate true stresses:

35

U1(i ) U (0,1); U 2(i ) U (0,1) where U1(i ) and U 2(i ) are IID byield yield (i ) = a yield + 1/ c yield (12) ln(U1(i ) )

bult ult (i ) = ault + 1/ cult ln(U 2(i ) )

2.2.3.2 Stochastic Model for Plastic-Flow Properties - Calculating (K, n) from estimates for yield and ult Dowling [34] provides a derivation of a method for estimating the true strain at the ultimate true stress based on a power-law representation of the plastic flow properties of the material. The following relation-ships between engineering stress and strain ( , ) and true stress and true strain ( , ) for deformation well beyond yielding hold:

= ln (1 + )

(13)

= (1 + )

where the overhead curl denotes a true stress or strain. The true strain is also related to the initial, Ai, and deformed, A, cross-sectional areas (for example in a tensile roundbar specimen) by Ai

= ln A

(14)

A

()

exp = i A

From the definition of the true effective stress, , we have P P A A

= = i = i A Ai A A (15)

A

= = exp ( )

Ai Differentiating the final expression in Eq. (15) with respect to the true strain, , gives d 1 d

= (16) d exp( ) d 36

At the ultimate stress, we apply the following constraint condition d 1 d ult

= ult = 0 (17) d = ult exp( ult ) d ult where the condition in Eq. (17) holds only when d ult

= ult (18) d ult Applying a power-law constitutive model, the constraint condition produces the following result n

= K d n 1

= Kn d

(19) n 1 n Kn ult = K ult

( )( )

n n = ult ult 1 n n = ult Therefore, the exponent of the power-law constitutive model can be viewed as an estimate of the true strain at ultimate.

From Eq. (13), the 0.2% offset engineering yield strain and engineering yield stress are related to the true strain and true stress by yield = ln (1 + yield )

(20) yield = yield (1 + yield )

From the definition of the 0.2% offset engineering yield strain, we have yield yield = 0.002 + (21)

E where E is the modulus of elasticity. From Eq. (20)

( )

exp yield = yield + 1 (22)

Combining Eqs. (21) and (22), we obtain 37

yield

( )

exp yield = 1.002 + (23)

E From Eq. (20), the engineering yield stress is related to the true yield stress and engineering yield strain by yield yield = (24)

(1 + yield )

Combining Eqs. (22) and (24) in Eq. (23) results in yield

( )

exp yield = 1.002 + (25)

(

E exp yield )

( )

Multiplying Eq. (25) by exp yield produces the following quadratic equation yield 2 1.002 = 0 (26)

E Solving Eq. (26) for its positive root gives 4 yield 1.002 + (1.002)2 +

(

exp yield )

2 E

or (27) 1.002 + (1.002) 2 + 4 yield E

yield = ln 2

Given the estimate for ult by Eq. (19), the relation for yield from Eq. (27), and sampled values of yield (i ) and ult (i ) (from the sampling protocols of Eqs. (9), (10), (11) or (12)), the values for the power-law coefficient, K(i), and exponent, n(i), can now be calculated by solving the following nonlinear system of equations for the two unknowns ( K (i ) , n(i ) ) :

38

n(i )

1.002 + 1.002 + 4 yield (i ) / E yield (i ) K (i ) ln 2

( )

=0 2

(28)

( )

n (i )

ult (i ) K (i ) n(i ) =0

( )

With the sampled values for K (i ) , n(i ) , the slope of the fragility curve can then be calculated from Eq. (5) or crit (n(i ) )

Papplied exp 2

S(i ) = (29) crit (n(i ) )

4 crit ( K (i ) , n(i ) ) exp[ crit ( n(i ) )] exp 1 2

and the sampled fragility curve is then, from Eq. (2):

(

a Atip (crit )(i ) = h0 + S(i ) x Rcavity ) (30)

Table 3 shows the results of one uncorrelated and three correlated bivariate samplings from the yield and ultimate stress stochastic models described by Eqs. (9)-(12). Figure 15 presents histograms (N=10,000) of the resulting sampled fragility curve slopes obtained with the three bivariate joint distributions and the one uncorrelated sampling protocol for yield and ultimate true stresses.

The objectives of multivariate sampling are: (1) to retain the characteristics of the individual marginal distributions for each component of the vector and (2) to reproduce the input correlation coefficient in the generated samples. If these two objectives are met, then the multivariate sampling protocol is said to maintain control over the intended multivariate joint probability distribution. The simulated data in Table 3 indicate that the bivariate lognormal protocol maintains better control over the resulting sampled bivariate distribution than the bivariate gamma protocol. The bivariate lognormal joint distribution was, therefore, chosen for the best-estimate case, and the other three options were included as special cases (cf.

Cases CCA-001, CCA-019, CCA-020, and CCA-021 in Table 10 of Sec. 3.2) in the sensitivity-study case matrix.

Figure 16a provides an example of the varying partitioning of the damage-state space produced by sampling from the bivariate lognormal joint distribution for yield and ultimate stress. The cumulative probabilities for the slope of the fragility curve have been estimated from order statistics with a sample size of 10,000. In Fig. 16b, random bivariate sampling produced a fragility curve associated with a specific cumulative probability of 97.5%.

39

{ }

T Table 3. Results of Bivariate Sampling Protocols for a Range of Correlation Coefficients, X i = yield (i ) , ult (i )

Case Sampling Sampled Input Distribution Statistics Simulation Statistics ( N = 10,000 ) Sampled Quantiles Marginal Distributions Joint Distribution No. Protocol Variate Mean Variance 2.50% Median 97.50% Correlation Mean Variance Correlation Covariance Matrix 2.50% Median 97.50%

Yield Stress inverted Weibull 31.259 2.651 28.99 30.96 35.27 31.272 2.611 2.611 4.825 yield 30.97 30.97 1 Fully Correlated bivariate inverted Weibull 1 0.99996 Ultimate Stress inverted Weibull 72.602 9.055 68.36 72.06 80.00 72.626 8.918 4.825 8.918 68.45 72.09 79.88 Yield Stress lognormal 31.010 3.515 27.71 30.96 34.60 30.993 2.986 2.986 0.324 27.71 30.97 34.49 2 Jones and Miller (1966) bivariate lognormal 0.06931 0.06994 Ultimate Stress lognormal 72.110 7.459 66.96 72.06 77.55 72.080 7.188 0.324 7.188 66.99 72.04 77.56 Arnold (1967) Yield Stress gamma 31.493 4.022 29.26 30.96 36.72 31.512 4.043 4.043 0.373 29.29 30.96 36.81 3 bivariate gamma 0.06931 0.05514 trivariate reduction Ultimate Stress gamma 72.906 11.711 69.10 72.06 81.85 72.932 11.301 0.373 11.301 69.08 72.06 81.48 Yield Stress inverted Weibull 31.259 2.651 28.99 30.96 35.27 31.267 2.670 2.670 -0.070 28.98 30.96 35.31 4 Uncorrelated NA 0 -0.01427 Ultimate Stress inverted Weibull 72.602 9.055 68.36 72.06 80.00 72.600 9.000 -0.070 9.000 68.38 72.05 79.90 (a) (b) (c) (d)

Fig. 15. Frequency histograms of sampled fragility curve slopes based on (a) perfect correlation between yield and ultimate stresses, (b) bivariate lognormal joint distribution with correlation coefficient of 0.06931, (c) bivariate gamma joint distribution with correlation coefficient of 0.06931, and (d) uncorrelated yield and ultimate stresses and, therefore, statistically independent. The correlation coefficient of 0.06931 was statistically inferred from the data in Table 1.

40

(a)

(b)

Fig. 16. The fragility curve stochastic model for Papplied = 2.165 ksi produces a slope with its uncertainty described by the distributions shown in Fig. 15. The intersection of the fragility curve with the y-axis is equal to the initial undeformed cladding thickness, h0 , which is assumed known in the current analysis: (a) a sampling of fragility curves based on their cumulative probability resulting from the bivariate lognormal joint distribution for yield and ultimate stress and (b) damage state partitioning based on an uncertain fragility curve slope with cumulative probability of 97.5%.

41

In Fig. 17, the predictions of the modified Chakrabarty and Alexanders instability theory [15] are compared to the pressures at numerical instability, PNI, predicted by 3-dimensional finite-element models of the 2-inch long model flaw centered in 0.25 inch thick burst disks with varying radius. These results indicated good agreement between the modified instability theory and the detailed finite-element models up to a relative flaw depth of a/t = 0.5.

For flaws with a/t 0.5, the modified instability theory becomes slightly nonconservative, i.e., the theory produces burst pressures greater than those estimated by the finite-element models. With reference to the BWXT forensic investigation [16], however, the deepest crack observed in the Davis-Besse wastage had an a/t < 0.5, and most of the cracks are significantly shallower; therefore, the above model is appropriate for this study.

2.2.4 Failure Mechanisms - Development of a Ductile-Tearing Model A stochastic model for ductile-tearing initiation was developed from data in [35, 36]. Table 4 presents ductile tearing data for three-wire series-arc stainless steel weld overlay cladding published in NUREG/CR-5511 [35]. Table 5a presents additional ductile tearing data for stainless steel cladding published in NUREG/CR-6363 [36] and data from pre-cracked Charpy V-notch (PCCVN) specimens made from PVRUF and Davis-Besse cladding materials [2]. The ductile-tearing data in Tables 4 and 5 are plotted as a function of temperature in Fig. 18. The JIc data at 288 °C (550 °F) from NUREG/CR-5511 have been extrapolated to 318.33 °C (605 °F) using a curve fit developed from the data in Table 4. The extrapolated JIc data are also presented in Table 5a.

The ExpertFit© statistical software [27] was used to fit several statistical distributions (see Table 5b) to different groupings of the JIc data in Table 5a. The resulting fit for the JIc data from NUREG/CR-5511

[35] and NUREG/CR-6363 [36] produced a log-logistic distribution with the following form:

42

(a)

(b)

(c)

Fig. 17. Net-section plastic collapse pressures estimated by modified C&A (1970) theory compared to PNI values determined from 3D finite-element ABAQUS solutions: (a) a/t = 0.05, (b) a/t = 0.125, (c) a/t = 0.25, (Ductile tearing not simulated for these cases.)

43

(d)

(e)

Fig. 17 (continued) Net-section plastic collapse pressures estimated by modified C&A (1970) theory compared to PNI values determined from 3D finite-element ABAQUS solutions: (d) a/t = 0.5, and (e) combined plot. (Ductile tearing not simulated for these cases.)

44

Table 4. Ductile-Tearing Data Extracted from Table 13 of NUREG/CR-5511 [35].

Test Tearing Specimen Temperature J Ic Modulus

(°C) (kJ/m2)

Unirradiated Specimens A13G -75 117 64 H2 -75 137 49 a

A15B 20 165 270 A13D 20 134 209 A10G 20 171 176 A10E 120 128 246 H5 120 119 229 H3 120 120 232 a

A13F 120 159 359 H6 200 90 240 H4 200 111 231 A15D 288 77 267 A13C 288 66 170 H1 288 82 192 Irradiated Specimens A15F -75 78 40 A15G -75 56 36 A13A 30 144 177 A15C 50 124 146 A10F 120 94 175 A15A 288 25 191 a

Specimen was not side-grooved, while all other specimens in table were side-grooved 20%.

45

Table 5a. Ductile-Tearing Data Used in Development of JIc Statistical Distributions J Ic Distribution Values J IC Average Specimen Test Temp. or J Q Reference Remarks Oper. Temp. J IC Tearing 2 2 ID (°C) (kJ/m ) (°C) (kJ/m ) Modulus A15D 288 76 NUREG/CR-5511 Three-wire cladding study 318 72.17 267 A13C 288 70 NUREG/CR-5511 Three-wire cladding study 318 66.47 170 H1 288 83 NUREG/CR-5511 Three-wire cladding study 318 78.82 192 H10 288 85 NUREG/CR-6363 Aged 3-wire cladding at 288 °C for 1605 hrs. 318 80.71 AA04 288 93 NUREG/CR-6363 Aged 3-wire cladding at 288 °C for 1605 hrs. 318 88.31 AA02 288 59 NUREG/CR-6363 Aged 3-wire cladding at 288 °C for 1605 hrs. 318 56.02 AA13 288 91 NUREG/CR-6363 Aged 3-wire cladding at 288 °C for 20,000 hrs. 318 86.41 AA15 288 77 NUREG/CR-6363 Aged 3-wire cladding at 288 °C for 20,000 hrs. 318 73.12 H15 288 111 NUREG/CR-6363 Aged 3-wire cladding at 343 °C for 20,000 hrs. 318 105.4 H16 288 110 NUREG/CR-6363 Aged 3-wire cladding at 343 °C for 20,000 hrs. 318 104.45 R3P1 315.6 166.2 ORNL/NRC/LTR-04/13 Pre-cracked Charpy Test Results PVRUF 316 166.2 N/E R3T1 315.6 197.2 ORNL/NRC/LTR-04/13 Pre-cracked Charpy Test Results PVRUF 316 197.2 214.71 R3T6 315.6 161.7 ORNL/NRC/LTR-04/13 Pre-cracked Charpy Test Results PVRUF 316 161.7 172.69 R4P2 315.6 121.6 ORNL/NRC/LTR-04/13 Pre-cracked Charpy Test Results PVRUF 316 121.6 242.13 R4T2 315.6 185.7 ORNL/NRC/LTR-04/13 Pre-cracked Charpy Test Results PVRUF 316 185.7 160.19 R4T5 315.6 184.2 ORNL/NRC/LTR-04/13 Pre-cracked Charpy Test Results PVRUF 316 184.2 N/E R5P2 315.6 143.7 ORNL/NRC/LTR-04/13 Pre-cracked Charpy Test Results PVRUF 316 143.7 130.56 R5T1 315.6 182.2 ORNL/NRC/LTR-04/13 Pre-cracked Charpy Test Results PVRUF 316 182.2 154.6 R5T5 315.6 216.5 ORNL/NRC/LTR-04/13 Pre-cracked Charpy Test Results PVRUF 316 216.5 181.15 R0T1 315.6 139.5 ORNL/NRC/LTR-04/13 Pre-cracked Charpy Test Results PVRUF 316 139.5 92.15 R0T2 315.6 123.7 ORNL/NRC/LTR-04/13 Pre-cracked Charpy Test Results PVRUF 316 123.7 76.79 R0T3 315.6 112.7 ORNL/NRC/LTR-04/13 Pre-cracked Charpy Test Results PVRUF 316 112.7 113.54 R0T5 315.6 158.3 ORNL/NRC/LTR-04/13 Pre-cracked Charpy Test Results PVRUF 316 158.3 104.26 R0T6 315.6 103.8 ORNL/NRC/LTR-04/13 Pre-cracked Charpy Test Results PVRUF 316 103.8 91.28 TS1 315.6 137.9 ORNL/NRC/LTR-04/13 Pre-cracked Charpy Test Results Davis-Besse 316 137.9 110.9 TS2 315.6 180.85 ORNL/NRC/LTR-04/13 Pre-cracked Charpy Test Results Davis-Besse 316 180.85 180.85 TS4 315.6 121 ORNL/NRC/LTR-04/13 Pre-cracked Charpy Test Results Davis-Besse 316 121 103.52 TS5 315.6 132.41 ORNL/NRC/LTR-04/13 Pre-cracked Charpy Test Results Davis-Besse 316 132.41 132.3 LS1 315.6 96.43 ORNL/NRC/LTR-04/13 Pre-cracked Charpy Test Results Davis-Besse 316 96.43 107.73 LS2 315.6 93.21 ORNL/NRC/LTR-04/13 Pre-cracked Charpy Test Results Davis-Besse 316 93.21 103.03 LS3 315.6 100.92 ORNL/NRC/LTR-04/13 Pre-cracked Charpy Test Results Davis-Besse 316 97.76 117.08 LS4 315.6 101.47 ORNL/NRC/LTR-04/13 Pre-cracked Charpy Test Results Davis-Besse 316 100.92 128.21 LS5 315.6 97.958 ORNL/NRC/LTR-04/13 Pre-cracked Charpy Test Results Davis-Besse 316 101.47 125.11 TL1 315.6 108.15 ORNL/NRC/LTR-04/13 Pre-cracked Charpy Test Results Davis-Besse 316 108.15 74.79 TL3 315.6 128.2 ORNL/NRC/LTR-04/13 Pre-cracked Charpy Test Results Davis-Besse 316 128.2 113.7 TL4 315.6 137.57 ORNL/NRC/LTR-04/13 Pre-cracked Charpy Test Results Davis-Besse 316 137.57 115.08 NUREG data extrapolated to 318.33 °C (605 °F).

Table 5b. Statistical Distributions Fitted to JIc Data from Three Sources J Ic Data Marginal Sample Point Estimates for Parameters Goodness of Test Critical Values for Level of Significance ()

Source Distribution Size location scale shape Fit Test Statistic 0.25 0.15 0.1 0.05 0.025 0.01 0.005 NUREG/CR-5511 Anderson-Darling 0.177 0.416 NA 0.549 0.644 0.75 0.884 0.985 NUREG/CR-6363 log-logistic 10 0 79.9486 9.3076 Kolmogorov-Smirnov 0.390 NA NA 0.679 0.73 0.774 0.823 NA (1/2T C(T)) equi-probable chi 2 NA NA NA NA NA NA NA NA Anderson-Darling 0.470 0.436 NA 0.585 0.698 0.810 0.96 1.075 Davis-Besse lognormal 12 0 4.7656 0.1909 Kolmogorov-Smirnov 0.673 NA 0.726 0.767 0.838 0.932 0.969 NA (PCCVN) equi-probable chi 2 2.667 4.108 5.317 6.251 7.815 NA 11.345 NA Anderson-Darling 0.212 0.450 NA 0.605 0.719 0.833 0.985 NA PVRUF Weibull 14 0 170.30 5.41 Kolmogorov-Smirnov 0.461 NA NA 0.760 0.819 0.880 0.944 NA (PCCVN) equi-probable chi 2 1.429 4.108 5.317 6.251 7.815 NA 11.345 NA 1/2T C(T) = pre-cracked compact tension specimens PCCVN = pre-cracked Charpy V-Notch specimen 46

(a)

(b)

Fig. 18. Ductile-tearing data for cladding from several sources: (a) JIc data and (b) average tearing modulus data.

47

Log-logistic Distribution for NUREG C(T) JIc Data - LL( , , ) with

= 9.30759 , = 79.94864 kJ/m 2 , and = 0 .

The log-logistic probability density and cumulative distribution functions are defined below.

1 x

if x >

2 density: f ( x) = x 1 +

0 otherwise 1

if x >

x CDF: F ( x) = 1 +

0 otherwise The corresponding percentile function is 1 p ln J Ic ( p l , , ) = + exp p

for 0 < p < 1 (31)

A lognormal distribution was fitted to the Davis-Besse PCCVN data (N = 12) where:

Lognormal Distribution for Davis-Besse PCCVN JIc Data - LN ( , µlog , log ) with

= 0, µlog(ult ) = 4.76562, log(ult ) = 0.19088 where µlog = lognormal mean (scale parameter), log = lognormal standard deviation (shape parameter),

and is the location parameter . The lognormal probability density and cumulative distribution functions are defined below.

ln( x ) µ 2 log exp 1

if x >

density: f ( x) = ( x ) 2 2 2 2 log log 0 otherwise 48

ln( x ) µlog if x >

z 1 2 CDF: F ( x) = log where ( z ) = 2 2 d exp 0 otherwise A Weibull distribution was fitted to the PVRUF PCCVN data (N = 14) where:

Weibull Distribution for PVRUF PCCVN JIc Data - W ( a, b, c) ) with location parameter, a = 0, scale parameter, b = 170.4 kJ/m2, and shape parameter, c = 5.41 where x a c c b c ( x a )c 1 exp f ( x) = if x > a density: b 0 otherwise x a c 1 exp F ( x) = b if x > a CDF:

0 otherwise Figures 19a and 19b present plots of the probability density and CDF of the lognormal distribution for the Davis-Besse PCCVN data. Figures 19c and 19d compare the fitted distributions for the three JIc samples, 0.5T C(T) data reported in NUREG/CR-5511 [35] and NUREG/CR-6363 [36] (log-logistic), PCCVN data from the PVRUF cladding (Weibull), and PCCVN data from the Davis-Besse cladding (lognormal).

49

(a)

(b)

Fig. 19. Statistical distributions for JIc : (a) probability density of lognormal distribution fitted to the 12 Davis-Besse PCCVN data points and (b) lognormal cumulative distribution function compared to cumulative probabilities of Davis-Besse PCCVN JIc data estimated by the median rank order statistic p = (i-0.3)/(n+0.4).

50

(c)

(d)

Fig. 19. (continued) fitted distributions for separate samples of JIc data from 0.5T C(T) clad overlay (NUREG/CR-5511 and 6363) and PCCVN specimens from PVRUF and Davis-Besse cladding materials: (c) probability densities and (d) cumulative distribution functions.

51

2.2.4.1 Ductile-Tearing Instability Analyses At a test temperature of 288 °C (550 °F), JR curves are available in [35] for three 1/2 C(T) (20% side grooved) unirradiated fracture specimens. One irradiated test is available at this test temperature but was not used in the analysis, since in the Davis-Besse problem the location of the wastage cavity is far away from the beltline of the RPV and, therefore, is only marginally affected by irradiation. The data-fitted JR curves from [35] and from the PCCVN Davis-Besse specimens are presented in Fig. 20. The applied driving-force curves calculated from finite-element models (see Fig. 21) for a range of flaw depths (all with 2.0 inch long surface flaws centered in the As-Found burst disk model) are shown in Fig. 22. The results from Fig. 20 and 22 can be used to carry out a tearing-instability analysis based on the theory described in [37].

2.2.4.1.1 A-Tip Ductile-Tearing Instability Two J-T A-tip instability analyses are shown in Figs. 23 and 24. Tearing-instability theory [37] employs J-T curves to define regions of stable and unstable ductile tearing. The JR-TR material curves are constructed from JR data for a given material at a specified test temperature where the tearing modulus is defined in [37] by E dJ R TR =

02 da J = J (32)

R The applied tearing moduli are calculated numerically from finite-element results in Fig. 22 such that E J Tapplied =

02 a J = J applied (33)

The differentiation required in Eq. (33) was carried out using central, forward, and backward 2nd-order finite-difference ratios derived for unequal partitions. From the theory discussed in [37], ductile tearing of a given flaw is predicted to become unstable when two conditions are met:

J applied J R and Tapplied TR (34) 52

(a)

(b)

Fig. 20. JR curves (a) fitted from data available in ref. [1] for unirradiated 0.5TCS fracture specimens at 550 °F and (b) comparing NUREG/CR C(T) JR curves with PCCVN JR curves.

A region of instability is therefore defined in (J,T) space to the right of an experimentally-derived (JR,TR) curve. A flaw becomes unstable when its (Japplied,Tapplied) locus intersects a materials (JR,TR) curve.

In Fig. 23 two applied pressures are examined, the nominal operating pressure of 2.165 ksi and the relief-valve setpoint pressure of 2.5 ksi. A range of flaw depths were also investigated. At the operating pressure, flaws with depths up to at least 0.0625 in. are predicted to be stable from the data from all three test specimens. At the set-point pressure, the flaw at a depth of 0.0625 in. is predicted to be unstable from the most conservative of the three test specimens, A13C.

53

(a)

(b)

(c)

Fig. 21. Finite-element models used in calculating applied J-integrals produced by pressure loading of burst disk: (a) (a/t = 0.5, 2L/a = 16) (b) (a/t = 0.25, 2L/a = 32), and (c) (a/t = 0.125, 2L/a =

64) 54

(d)

(e)

Fig. 21. (continued) Finite-element models used in calculating applied J-integrals produced by pressure loading of burst disk: (d) (a/t = 0.05, 2L/a = 160) (e) (a/t = 0.8, 2L/a = 10) 55

3.5 Burst Disk Dimensions a = 0.0625 in diameter = 6.0 in 3.0 flaw length = 2.0 in a = flaw depth 2.5 Japplied (in.-kips/in )

2 a = 0.172 in a = 0.05 in 2.0 a = 0.125 in 1.5 a = 0.03125 in 1.0 a = 0.1 in a = 0.025 in 0.5 a = 0.0125 in 0.0 0 1 2 3 4 5 6 7 8 Pressure (ksi)

Fig. 22. J-integral driving forces - applied pressure as a function of JI and a/t for a 2-inch long flaw.

56

Fig. 23. A-tip tearing instability plots with JR vs TR curves from ref. [35] and Japplied vs Tapplied curves from finite-element simulations of burst disks with 2-inch long surface flaws centered in the burst disk. Applied J-T curves are for two applied pressures (the nominal operating pressure of 2.165 ksi and the relief-valve setpoint pressure of 2.5 ksi) with varying flaw depths 0.0125 in. aA-tip 0.125 in.

57

Fig. 24. A-tip tearing instability plots with JR vs TR curves from ref. [35] and Japplied vs Tapplied curves from finite-element simulations of burst disks with 2-inch long surface flaws centered in the burst disk. Applied J-T curves are for a range of flaw depths ( 0.0125 in. a 0.125 in. ) and varying applied pressures.

Fig. 25. C-tip tearing instability plots with JR vs TR curves from ref. [35] and Japplied vs Tapplied curves from finite-element simulations of burst disks for As-Found flaw with initial 2-inch long with A-tip collapse.

58

In Fig. 24, seven flaws (at depths of 0.0125 in. a 0.125 in. ) are exposed to a range of applied pressures. Above approximately 2.6 ksi, the deepest flaw is predicted to fail by plastic collapse at a burst pressure of 3.4 ksi. The other six flaws all fail by plastic collapse at pressures higher than 3.4 ksi. The smaller flaw (a = 0.0125 in.) is seen to remain stable in Fig. 24 over the full range of applied pressures based on all three test specimens.

2.2.4.2 C-Tip Ductile-Tearing Instability Figure 25 presents the results of a tearing-instability analysis of the C-tip of the model flaw after the assumed collapse of the A-tips remaining ligament for the As-Found condition. This analysis indicates that for pressures below approximately 3.25 ksi, the C-tip will not experience unstable ductile tearing.

The wastage-area damage model analysis carries out a C-tip ductile-tearing instability analysis after failure of the A-tip remaining ligament by plastic collapse. A series of finite-element models (see Fig. 26a for an example) were constructed for varying values of Rcavity with a failed A-tip flaw with constant 2L =

2.0 in. Each model was loaded up to an applied pressure of 2.165 ksi, corresponding to the nominal operating pressure of the Davis-Besse RPV. The resulting applied driving forces, Japplied, and driving force gradients, dJ/da, are shown in Figs. 26b and 26c, respectively, for a range of cavity radii.

To apply the tearing-instability test of Eq. (34), estimates of the claddings tearing modulus, TR, must be developed as a function of JR. If the JR curves of Fig. 20 can be characterized by the power-law form J R = C (a m ) (35) then the local material tearing modulus, TR, is from Eq. (32)

E dJ R E ( m 1)

TR = 2 = 2 x m x C x a (36) f f da Given the elastic modulus, E, and estimated flow stress, f , the remaining three variables required by the ductile-tearing model are JIc, C, and m. Applying the definition of JIc in ASTM E-1820 [38], sampled estimates of two of these parameters (JIc and m) allows the calculation of the third parameter, C. In Fig. 27, the ductile-tearing initiation toughness, JIc , is defined in ASTM E-1820 as the intersection of the JR curve with a 0.2 mm offset blunting line given by J (0.2 mm offset) = 2 f ( a a0 ) (37) 59

(a)

(b)

(c)

Fig. 26. Driving-force curves used in C-Tip ductile-tearing instability analyses: (a) models were developed for varying Rcavity and constant 2L = 2a = 2.0 in, tclad = 0.25 in., (b) applied driving force Japplied at C-tip at 2.165 ksi for varying Rcavity, and (c) applied dJ/da at C-tip at 2.165 ksi for varying Rcavity.

60

Fig. 27. Given a JR curve in power-law model form and estimated flow stress, f , the initiation toughness, JIc , and local tearing modulus, TR, are uniquely defined (see ASTM E-1820).

61

where the prescribed offset is a0 = 0.2 mm (0.008 in) . Therefore, with an estimate of JIc and the power-law exponent, m, the power-law coefficient, C, is J Ic J Ic = C a m C =

a m J Ic J Ic = 2 f (a a0 ) a = + a0 (38) 2 f J Ic C= m J Ic

+ a0 2 f The local tearing modulus then follows from Eq. (36). In the wastage-area damage model, the values of JIc and m are sampled from a bivariate lognormal distribution to be described in Sect. 2.2.4.4.

2.2.4.3 Model for JR Curve Power-Law Exponent, m A lognormal distribution was developed to characterize the uncertainty in the power-law exponent, m, (see Eq. (35)) used in constructing a JR vs a curve (see Fig. 28). A mean value for m of 0.68141 was estimated from specimen H1 in the data presented in NUREG/CR-5511 [35]. In Table B15(a) of NUREG/CR-6004 [22], a variance of 0.0159 is given for m for a stainless steel flux weld. Covariance data for correlated sampling of JIc and m are also presented in the Table B15(a):

Table B15(a): Stainless steel flux weld 8.141

= 2.756 x 10 4

8.141  ; µ J Ic =194.65; µm =0.733; and 0.0159 8.141 12 = = 0.3889

( )

2.756 x 104 ( 0.0159 )

log-transform of covariance matrix (39) ij ij (log) = ln 1 +

µi µ j (log) =

0.5466189 -0.05875121

-0.05875121 0.02916359 and

-0.05875121 12(log) = = 0.46532243

( 0.5466189 ) ( 0.02916359 )

A lognormal correlation coefficient of -0.4653224 was calculated from the data in NUREG/CR-6004 and then applied in the construction of a lognormal bivariate distribution for sampling correlated (JIc, m) data pairs as described in the following section.

62

(a)

(b)

Fig. 28. Lognormal distribution characterizing the uncertainty in the JR curve power-law exponent,

m. The lognormal distribution with (a) probability density and (b) CDF has been fitted using a mean value from the data reported in NUREG/CR-5511 for overlay cladding and a variance from data in NUREG/CR-6004 for stainless steel flux weld.

63

2.2.4.4 Bivariate Lognormal Distribution - 12 = -0.388901; 12(log) = -0.465322 (correlation between JIc and m inferred from data in NUREG/CR-6004, Table B.15(a) for stainless steel flux weld)

The sampling of a multivariate lognormal random vector follows from the procedures discussed in

[30, 31].

Estimate Log-Transformed Covariance Between JIc and m The first step requires the estimation of the log-transformed covariance between JIc and m. From Fig. 19, JIc (log) = 0.19088; JIc 2 (log) = 0.036435 and from Fig. 28, m (log) = 0.1835; m (log) = 0.033672 2

Applying the log-transformed correlation coefficient calculated in Sect. 2.2.4.3:

cov( J Ic , m)(log) = 12(log) 2J (log) m 2 (log) = -0.0162981 Ic where (40) 12(log) = 0.46532243 Perform Cholesky Factorization The second step is to carry out a Cholesky decomposition (factorization) of the log-transformed corre-lation matrix.

2J (log) cov( J Ic , m)(log) 0.036435 -0.0162981 (2x2)(log) = Ic =

-0.0162981 0.033672 (41) cov( J , m) m(log) 2 Ic (log) and 2 x 2(log) = CCT where C is the 2 x 2 lower triangular Cholesky factor c11 c12 0.19088 0 C= =

c21 c22 -0.08538 0.16242 64

Bivariate Sampling Protocol for { J Ic ,m}

T Finally, for a bivariate lognormal joint distribution, the protocol for the correlated sample X i ={ J Ic (i ) ,m(i )}

T is as follows:

Given: µ J Ic (log) = 4.76562; µm(log) = 0.4334; c11 = 0.19088, c22 = 0.16242; c21 = 0.08538 Z 1(i ) N (0,1); Z 2(i ) N (0,1) where Z 1(i ) and Z 2(i ) are IID (42)

(

Y 1(i ) = µ J Ic (log) + c 11x Z 1(i ) )

Y 2(i ) = µm(log) + ( c21 x Z 1(i ) ) + ( c22 x Z 2(i ) )

exp(Y 1(i ) ) J Ic (i )

Xi = =

exp(Y 2(i ) ) m(i )

where IID refers to independent and identically distributed random variates, and Z1 and Z2 are standard-unit-normal variates.

2.2.4.5 Driving Force (Japplied) Surface Function of (aA-tip, Rcavity)

Based on the results of a matrix of finite-element solutions for Japplied as a function of varying aA-tip and Rcavity , a continuous surface function was fitted using TableCurve 3Ds Selective Subset Algorithm [39].

The resulting 7-parameter selective-subset surface function for a constant applied pressure of 2165 psi had the following form (see Fig. 29)

( )

ln J applied = C0 +

C1 a Atip

+

C2 a Atip

( )

+ C3 Rcavity + C4 Rcavity ln Rcavity +

(43)

( ) ( )

2 C5 ln Rcavity Rcavity + C6 ln Rcavity with Rcavity [in], aA-tip [in], and Japplied [kJ/m2] for 2.165 ksi applied lateral pressure.

65

Fig. 29. Plot of Japplied as a function of Rcavity and aA-tip for a burst disk with a 2 inch centered flaw under 2.165 ksi applied pressure. Japplied solutions from FEM models (using best-estimate stress vs strain data) were fitted to a surface function by TableCurve 3D.

or

( )

J applied a Atip , Rcavity = exp C0 +

C1 C

(

+ 2 + C3 Rcavity + C4 Rcavity ln Rcavity +

a Atip a Atip

)

(

C5 ln Rcavity ) ( )

Rcavity + C6 ln Rcavity 2

} (44)

(

J applied a Atip , Rcavity ) =

1 C1

+ 2 C

exp C0 +

C1 C

+ 2 + C3 Rcavity +

a Atip a Atip 2 a Atip a Atip a Atip a Atip Rcavity

( ) (

C4 Rcavity ln Rcavity + C5 ln Rcavity ) ( )

Rcavity + C6 ln Rcavity 2

}

66

with the fitted parameters C0 = 838.324387886275 C3 = -841.0282062799881 C1 = -3.944703731724143 C4 = 112.5438397120073 (45)

C2 = 0.6780555298560729 C5 = 755.9922572599784 C6 = -97.09961384574772 and kJ N J applied 2 or m mm a Atip [ mm ] (46)

Rcavity [in.]

2.2.5 Assumed Wastage-Area Damage Model Accident Sequence As previously discussed, an accident sequence is a combination of events leading from an initiating event to an undesired consequence. The sequence is ordered, starting with the initiating event, and proceeds through sequential failures leading to the consequence [3]. The consequences in this study are formulated in terms of loss-of-coolant-accidents (LOCAs) ranked by break size. The time-line for this analysis is shown in Fig. 30, where seven distinct time stations are identified: (1) cavity-ini , the time when the wastage cavity initiates, (2) flaw-ini, the time when the flaw initiates, (3) 0 , one year before the discovery of the wastage cavity, (4) DT-ini , the time of ductile-tearing initiation, (5) 1, the time of discovery (February-March, 2002), (6) 2, six weeks after the time of discovery, and (7) 3, one year after the time of discovery. With the inclusion of the ductile-tearing model into the analysis, an additional time station is introduced into the problem, the time of initiation of ductile tearing, DT-ini.

For the Davis-Besse Wastage-Area Damage Model, the damage state {now defined by the five parameters [Rcavity, aA-tip, Japplied(aA-tip, Rcavity), Tapplied(aA-tip, Rcavity), TR] at different points in time} can be calculated from sampled values for the cavity growth rate, dRcavity / d , flaw growth rate, da Atip / d ,

and ductile-tearing initiation fracture toughness, JIc.

2.2.5.1 Cladding-Capacity Analysis - Damage State Determination (without A-tip ductile tearing)

In the cladding-capacity analysis, the as-found damage state is assumed known to a sufficient degree of certainty that the (R1, a1) is represented by a fixed point (see Fig. 31) where at 1 67

Fig. 30. Time line for cladding capacity analysis - accident sequence of wastage-area damage model.

Fig. 31. For the cladding-capacity analysis, the As-Found damage state is assumed known and fixed. For constant cavity and flaw growth rates, the accumulation of damage follows, in the absence of ductile tearing, a linear path towards the fragility curve.

68

R1 = 2.443 inches; corresponding to an unbacked wastage area of 18.75 in 2 (47) a1 = 0.065 inches; corresponding to current best-estimate of prototypic flaw depth Cavity growth rates, dRcavity / d , and flaw A-tip growth rates, da Atip / d , are then sampled 8 from statistical distributions derived from an Expert Elicitation process to be discussed.

dRcavity D( dR / d )

d (48) da Atip D( da / d )

d The elapsed time since cavity initiation, cavity ini , can then be calculated from the known value for R1 by R1 cavity ini = (49) dRcavity d

The elapsed time since flaw initiation, flawini , is calculated from the time required to grow the flaw at 1 and is constrained to be no greater than the elapsed time since cavity initiation, cavity ini ,

flawini = min a1

, cavity ini (50) da Atip d

( )

The damage state at the time of flaw initiation, R flawini ,0 , is then dRcavity R flawini = R1 flawini x d (51) a flawini 0 The damage states at 2 and 3 are calculated by dRcavity R2 = R1 + 2 x d

(52) da Atip a2 = a1 + 2 x d

8 The notation X D should be read as the random variate, X , is distributed as the statistical distribution D.

69

dRcavity R3 = R1 + 3 x d

(53) da Atip a3 = a1 + 3 x d

where 2 and 3 have been previously defined 6

2 = 0.11538 year 52 (54) 3 1 year 2.2.5.2 Cladding-Capacity Analysis - Damage State Determination (with A-tip ductile tearing)

As discussed in Sect. 2.2.5.1, the as-found damage state in the cladding-capacity analysis is assumed known to a sufficient degree of certainty that the (R1, a1) is represented by a fixed state point as shown in Fig. 32. For this purely illustrative example, we selected the values R1 = 3 inches; a1 = 0.05 inches .

Cavity growth rates, dRcavity / d , and flaw A-tip growth rates, da Atip / d , are then sampled from statistical distributions derived from an Expert Elicitation process.

dRcavity D( dR / d )

d (55) da Atip D( daAtip / d )

d The position of the model flaw tip is assumed to be a function of the contribution due to the elapsed time of exposure to the corrosive environment after cavity initiation, aA-tip(env), and possibly due to a contribution from stable ductile tearing a Atip (env ) + a Atip ( DT )

( )

for J applied a Atip , Rcavity J Ic a Atip = (56) a Atip (env ) for J applied ( a Atip , Rcavity ) < J Ic where

( )

da Atip a Atip (env ) ( ) = flawini (57) d The elapsed time since cavity initiation, cavity ini , is calculated from the known value for R1 by 70

(

cavity ini = 1 cavity ini = ) R1 dRcavity (58) d All elapsed times are measured relative to the time of cavity initiation, cavity-ini = 0.

Determination of the Time of Flaw Initiation, flaw-ini If the damage state at 1 = R1 / ( dRcavity / d ) is such that J applied ( a1 , R1 ) < J Ic (59) then the time of flaw initiation is calculated directly by a1 flawini = 1 (60) da Atip d

However if J applied (a1 , R1 ) J Ic , then the time of flaw initiation, flawini , is calculated as the solution to the following nonlinear equation da Atip a1 (1 flawini ) aDT (a1 , R1 ) = 0 d

where (61) 1 1 applied J ( a , R ) m J Ic m aDT ( a, R ) =

for J applied ( a, R ) J Ic C C 0 J applied ( a, R ) < J Ic where the sampled JR curve has been offset such that, at the point of ductile tearing initiation, the apparent initial flaw extension due to blunting is ignored, and at Japplied = JIc , aDT = 0 .

71

Fig. 32. An example flaw growth history (with ductile tearing) calculated using MathCad. The case conditions are based on the median values from the more conservative (MC) sampling distribution for dR / d , best-estimate distribution for da / d , JIc, and m. The As-Found condition is assumed known (Rcavity , aA-tip) = (3.0, 0.05) for this example case.

The elapsed time since flaw initiation, flawini , (relative to cavity initiation) is constrained to be no greater than the elapsed time since cavity initiation, cavity ini ,

flawini (c ) = min flawini , cavity ini or (62) flawini 0; cavity ini 0 72

Determination of the Time of Ductile-Tearing Initiation, DT-ini The time of ductile-tearing initiation, DT-ini , is calculated by solving the following system of nonlinear equations

( )

J Ic J applied a Atip ( DT ini ) , Rcavity ( DT ini ) = 0 1

(

da Atip J applied a Atip ( DT ini ) , Rcavity ( DT ini ) ) m J Ic m = 0 1

(

a Atip ( DT ini ) DT ini flawini ) d C C for the flaw tip position, aA-tip(DT-ini0, and DT-ini .

Determination of the flaw-tip position, aA-tip, for Times After Ductile-Tearing Initiation, DT-ini After ductile-tearing has been initiated, flaw tip position, aA-tip, can be calculated by solving the following nonlinear equation 1

1

(

da Atip J applied a Atip , Rcavity ( ) ) m J m a Atip

(

flawini ) d C

Ic = 0 C

(64) dRcavity

(

where for a given time, > DT-ini, Rcavity ( ) = cavity ini ) d Figures 32-34 present the results of an example case using the median values of da Atip / d and dRcavity / d from their best-estimate (BE) distributions.

73

Another critical time to be determined is the point at which the flaw-cavity growth path intersects the fragility curve. This time, * , is the solution of the following nonlinear equation set:

1 1

(

da Atip J applied a Atip , Rcavity ( )

) m Ic = 0

( ) J m a*Atip

• flawini d C C (65)

( )

a*Atip h0 S x Rcavity ( * ) = 0 For a cladding initial thickness of h0 = 0.25 in and a fragility curve slope of S = -0.026117, the intersection point on the fragility curve in Fig. 34 is a*Atip = 0.1701 in; R* =3.058 in; and * = 2.162 years .

The Davis-Besse damage model makes the assumption that the uncertain sampled cavity and flaw growth rates are constant over time. From this assumption and in the absence of ductile tearing of the A-tip, the cavity/flaw growth path, as shown in Fig. 31, follows a straight line approaching the fragility curve with increasing time from flawini . If ductile tearing of the A-tip occurs, then the cavity/flaw growth path becomes nonlinear. At the point of intersection of the cavity/flaw growth path with the fragility curve, the remaining ligament under the A-tip portion of the flaw is predicted by the modified Chakrabarty and Alexander [15] instability theory to collapse by plastic tensile instability for an applied pressure of 2.165 ksi (i.e., the nominal operating pressure for the Davis-Besse RPV). The damage states below the fragility curve are assumed to be stable conditions with no LOCA. At the intersection of the cavity/flaw growth path with the fragility curve, the flaw configuration in Fig. 35a is assumed to instantaneously convert to the configuration shown in Fig. 35b, and the damage state projects vertically to the final damage state, designated as ( R*, h0 ) for example in Fig. 31.

74

(a)

(b)

Fig. 33. Evolution of the (a) driving force, Japplied, and (b) tearing modulus, Tapplied for the example problem.

75

Fig. 34. Failure-assessment diagram using the damage-state growth path calculated for the example problem.

76

(a)

(b)

Fig. 35. Schematic descriptions of the model flaws (a) A-tip and C-tip before A-tip ligament collapse and (b) C-tip location after A-tip ligament collapse.

77

The damage state at the intersection of the cavity/flaw growth path (defined by the line connecting the

( )

damage states R flawini ,0 and ( R2 , a2 ) projected ) is designated ( R*, a*) in Fig. 31. This intersection point is calculated by da h0 + Atip x R flawini dRcavity R* = max ,0 da Atip (66)

S dRcavity a* = max h0 + ( S x R *) ,0 da Atip where S is the uncertain slope of the fragility curve sampled from Eq. (29) and is calculated dRcavity from the uncertain sampled cavity and flaw growth rates da Atip da Atip d

= (67) dRcavity dRcavity d

2.2.5.3 ASP Analysis - Damage State Determination For the ASP analysis, the damage state at the time of discovery, 1 , is treated as uncertain. This uncertainty is expressed in terms of a statistical distribution for the elapsed time since flaw initiation, flawini , relative to the time of discovery (see the Time Line in Fig. 30) and an uncertainty in the cavity size at time station 0 .

flawini D( flawini )

(68) flawini (c ) = min flawini , cavity ini R0 D( R0 ) (69) where the damage state at 1 is now calculated by 78

dRcavity R1 = R0 + 0 x d

(70) da Atip a1 = min max flawini (c ) x ,0 , h0 d

and the values for 0 and h0 are fixed at 0 1 year h0 0.25 inch The damage states at time stations flawini , 2 , and 3 then follow from Eqs. (51)-(54) and (56)-(70).

2.2.6 Results of an Informal Expert Elicitation An informal expert elicitation was carried out within the NRC to provide expert judgments that were then used to develop statistical distributions for the following random variates for both the CCA and ASP studies: (1) the cavity wastage rate in terms of the time rate of change of the effective cavity radius, dR / d , and (2) the rate of A-tip flaw growth, da / d , due only to the exposure of the flaw to a corrosive environment, i.e., da / d does not include flaw growth due to ductile tearing. For the ASP study only, two additional random variate distributions are required: (3) the effective cavity radius at 1 year before discovery, R0, and (4) the time of flaw initiation, flawini , relative to the time of discovery. See Sect. 2.2.5.3 for a description of how these two distributions are applied in the ASP study.

For this elicitation, three subject-matter experts were asked to provide a best estimate (assigned a probability of 50%), a low estimate (assigned a probability of 5%), and a high estimate (assigned a probability of either 95% or 99.9%) for the four random variates described above. These expert judgments (see below) were then combined using the following aggregation procedure. In the case of the cavity wastage rate, dR / d , and flaw initiation time, flawini , two bins were constructed. The dR / d (1) bin contained the extremes of the aggregated sample from the minimum of 0.5 in/yr to 7.0 in/yr. The second bin, dR / d (2), reflected a more optimistic view of the wastage rate and was used only as a sensitivity check case (CCA-010) in the CCA study and not used at all in the ASP study. For the flaw initiation time, flawini , the estimates from Subject-Matter Expert No. # 2 were used exclusively in the ASP study. For R0 and da/d, the highest and lowest values established the selected upper and lower bounds.

79

These six cases were then fitted to a range of distributions including triangular, normal, lognormal, beta, Weibull, and uniform. The final aggregate estimates and distribution fittings are presented in Table 6 and Figs. 36-40.

Judgments Obtained from Informal Expert Elicitation Subject Lower Median (BE) Upper Parameter Parameter Matter Bound Value Bound Description Units Expert Associated Percentiles ID No.# 5% 50% 95%/99.9%

Effective Cavity #1 0 1.25 2 Radius at TOD-1 (in.) #2 0 1.125 2.5 R0 #3 0 1.5 2.25 Effecive Cavity #1 1 2 7 Wastage Rate (in./year) #2 0.5 2 6 dR/d #3 0.75 1.5 (-)

Flaw Initiation #1 12 36 48 Time w.r.t. TOD (months) #2 1 6 60 flaw-ini #3 (-) (-) (-)

Effective Flaw #1 0.001 0.01 0.1 Growth Rate (in./month) #2 0.001 0.01 0.1 da/d #3 0.004 0.01 0.04 80

Table 6. Summary of Candidate Input Sampling Distributions Davis-Besse CCA/ASP Studies - Aggregation of Expert Judgments Random Cavity Size Cavity Wastage Rate Flaw Initiation Time Flaw Growth Variates at 0 , R 0 dR /d (1) dR /d (2) flaw-ini (1) flaw-ini (2) Rate, da /d Percentiles (in) (in/yr) (in/yr) (months) (months) (in/month) 5% 0 0.5 0.1 12 1 0.001 50% (median) 1 2 1 36 6 0.01 95% or 99.9% 2.5 7 3 48 60 0.1 Distributions Distribution Parameters, Percentiles, and Moments Triangular lower bound 0.00 -0.6062 -0.4460 2.6509 -5.0994 upper bound 3.1937 8.7849 3.7679 55.0697 61.9364 shape parameter 0.1801 2 0.9689766 36 6 5% 0 0.5 0.1 12 1 50% 1 3.4105 1.3440 32.22 18.9 95% or 99.9% 2.5 7 3 48 60 mean 1.6368 3.3929 1.4303 31.2402 20.9457 variance 0.4037 3.9173 0.7665 117.3208 215.1634 Normal mean 1.2501 29.9999 standard deviation 0.7599 10.94318 5% 0 12 50% 1.2501 30.00 95% 2.5 48 variance 0.5775 119.7533 Lognormal logmean 0.6264 1.4223 -5.3080 logstandard deviation 0.8022 0.8647 0.9726 5% 0.5 1 0.001 50% 1.87082 4.1465 0.01 95% or 99.9% 7 60 0.1 variance 0.6435 0.7477 0.9459 Beta shape parameter 1 1.0614 0.6478 0.8605 2.0712 0.4451 0.4621 shape parameter 2 1.9386 2.3522 2.1395 0.9288 2.5549 2.5379 lower bound 0.0000 0.4611 0.0429 -0.1129 0.9738 0.0009 upper bound 3.1435 10.6305 4.0986 48.9374 67.1177 0.1115 5% 0.0989 0.5 0.1 12 1 0.001 50% 1 2 1 36 6 0.01 95% or 99.9% 2.5 7 3 48 60 0.1 mean 1.1121 2.6571 1.2062 33.7516 10.7866 0.0180 variance 0.0572 4.3773 0.0511 128.5645 138.1917 0.0004 Weibull location parameter 0.0000 0.3646 -0.0612 0.6492 0.0003 scale parameter 1.2579 2.3224 1.3836 7.8524 0.0141 shape parameter 1.5974 1.0451 1.3816 0.9555 0.9864 5% 0.196 0.5 0.1 1 0.001 50% 1 2 1 6 0.01 95% or 99.9% 2.5 7 3 60 0.1 mean 1.0920 2.6463 1.2025 8.663342 0.0144 variance 0.5933 4.7693 0.8574 70.390404 0.0002 Uniform lower bound -0.00422 upper bound 0.10010 5% 0.001 50% 0.0479 99.9% 0.1 mean 0.0479 variance 0.0009 81

Statistical Distributions Used in the Monte Carlo Analysis A standard uniform distribution on the interval U(0,1) is the starting point for all of the transformation methods that draw random variates from nonuniform continuous distributions. A uniform distribution is defined by the following:

Uniform Distribution - U(a,b) 0 ; x<a 1

PDF: fU ( x l a , b ) =  ; a xb b a 0 ; x > b 0  ; x<a x a CDF: Pr( X x) = FU ( x l a, b) =  ; a xb b a 1  ; x>b Moments for a uniform distribution:

a+b Mean µ=

2 (b a )2 Variance 2 =

12 Sampling from a 2-parameter uniform distribution: U i U ( a, b)

Sampling from a standard uniform distribution, U(0,1), is accomplished computationally with a Random Number Generator. The random number generator (available from the ranlib statistical library (function RANF() ) used in this analysis is based on a combination of two multiplicative linear congruential generators [40] with a combined period of 2.3 x 1018 . The sampled uniform deviate can then be scaled to the required range by Pi U (0,1)

(71)

X i = a + ( b a ) Pi 82

(a)

(b)

Fig. 36. Candidate sampling distributions for cavity size at one year before discovery, 0 :

(a) probability densities and (b) cumulative probabilities. (ASP analysis only).

83

(a)

(b)

Fig. 37. Candidate sampling distributions for cavity wastage rate, dR / d : (a) probability densities for two beta distributions, (b) probability densities for two triangular distributions, 84

(c)

(d)

Fig. 37. (continued) Candidate sampling distributions for cavity wastage rate, dR / d :

(c) probability densities for two Weibull distributions, and (d) cumulative probabilities.

85

(a)

(b)

Fig. 38. Candidate sampling distributions for the Case 1 flaw initiation time , flaw-ini (1) :

(a) probability densities and (b) cumulative probabilities. (ASP analysis only).

86

(a)

(b)

Fig. 39. Candidate sampling distributions for the Case 2 flaw initiation time , flaw-ini (1) :

(a) probability densities and (b) cumulative probabilities. (ASP analysis only).

87

(a)

(b)

Fig. 40. Candidate sampling distributions for flaw growth rate , da / d : (a) probability densities and (b) cumulative probabilities.

88

Triangular Distribution - T(a,b,c)

A triangular distribution is commonly applied as a rough model of a random variable when there are essentially no data available.

2( x a)

(b a)(c a) if a < x c 2(b x)

Density: f ( x l a , b, c ) = if c < x < b (b a)(b c) 0 otherwise

( x a)2 if x a (b a )(c a )

(b x) 2 CDF: Pr ( x X ) = 1 if a < x c (b a)(b c) 1 if b x Moments:

a+b+c Mean µ=

3 a 2 + b 2 + c 2 ab ac bc Variance 2 =

18 Sampling from a 3-parameter triangular distribution: X i T ( a , b, c )

A random variate is sampled from a 3-parameter triangular distribution T(a,b,c) by Pi U (0,1)

(c a )

a + Pi (b a)(c a) if Pi (b a)

Xi = for a < c < b b (1 P )(b a )(b c) otherwise i

89

Beta Distribution - (1 , 2 , a, b )

y (1 1) (1 y )( 2 1) xa for 0 < y < 1; y =

f ( x l a, b,1 , 2 ) = B (1 , 2 ) ba 0 otherwise PDF:

B (1 , 2 ) = t ( 1 ) (1 t )( 2 ) dt 1 1 1 0

where there is no closed form relation for the cumulative distribution function of a beta distribution.

Moments:

1 Mean µ = a + (b a )

1 + 2 1 2 Variance 2 = (b a) 2 (1 + 2 )2 (1 + 2 + 1)

Estimates for the two shape parameters (1 , 2 ) can be developed from an estimated range (a, b) and mode c [33] by

( µ a )( 2c a b ) ( b µ )1 (a + b + c) 1 ( a, b, c )  ; 2 ( a , b, c )  ; µ ( a , b, c )

( c µ )( b a ) (µ a) 3 If we define Q1 ( p,1 , 2 ) as the percentile function9 for the standard beta distribution and assign probabilities for estimates for a lower bound, LB (e.g., pLB = 0.05), a median, MV (e.g., pMV = 0.5), and an upper bound, UB (e.g., pUB = 0.95), then the parameter set ( a, b, c, µ , 1 , 2 ) can be determined by solving the following nonlinear system of equations:

a + ( b a ) x Q1 ( pLB ,1 , 2 ) LB = 0 a + ( b a ) x Q ( pMV ,1 , 2 ) MV = 0 1

a + ( b a ) x Q ( pUB ,1 , 2 ) UB = 0 1

1

( µ a )( 2c a b ) = 0 (c µ )

2

( b µ )1 = 0

(µ a)

µ (a + b + c) = 0 3

9 The percentile function is derived by inverting the cumulative distribution function such that the random variate is expressed as a function of the cumulative probability and the parameters of the distribution.

90

Sampling from 4-term beta distribution - X i (1 , 2 , a, b)

The inverse of a standard beta distribution is calculated in this analysis by the dcdflib FORTRAN statistical routine, CDFBET [41] (designated as Q1 ( Pi ,1 , 2 ) ) and then scaled to the required nonstandard distribution by Pi U (0,1)

X i = a + ( b a ) x Q-1 ( Pi , 1 , 2 )

Normal Distribution - N ( µ , )

1 ( x µ )2 PDF: f N ( x l µ , ) = exp  ; < x < +

2 2 2 1

z 2 xµ CDF: Pr( X x) = ( z ) =

2 2 d ; z =  ; < x < +

exp Moments:

Mean µ Variance 2 Sampling from a 2-parameter normal distribution: X i N ( µ , )

In this analysis, the dcdflib FORTRAN subroutine CDFNOR is called to sample, Zp , from a standard unit-normal distribution. The standard normal deviate is then scaled to obtain the required quantile by X p = Z p + µ (72)

Lognormal Distribution - LN ( µlog , log )

0  ; x0 2

PDF: f ( x l µlog , log ) =

1 (

ln x µlog )

x 2 exp 2  ; 0< x<

log 2 log 0  ; x0 CDF: Pr( X x) = ( z ) = 1 z 2 ln x µlog 2 exp 2

d ; z =

, 0< x<

log 91

Moments:

2 log Mean

µ = exp µlog +

2 Variance 2 = ( 1) exp(2µlog ); = exp( log 2 )

Sampling from a 2-parameter Lognormal Distribution: X i LN µlog , log ( )

The log-transformed deviate is sampled from a standard unit normal distribution using dcdflib FORTRAN routine CDFNOR [28]. With the mean equal to the lognormal mean, µ log , and standard deviation equal to the lognormal standard deviation, log , the log-transformed deviate is then converted into the required random deviate by the exponential function.

Yi N ( µlog , log )

(73)

X i = exp(Yi )

Weibull Distribution - W(a,b,c)

(a = location parameter, b = scale parameter, c = shape parameter) 0  ; xa f W ( x a , b, c ) = c

( )

PDF:

c 1 exp y c ; ( y = ( x a ) / b, x > a, b, c > 0) y b

0 ; xa CDF: Pr( X x) = FW ( x l a, b, c) = c 1 exp y  ; ( y = ( x a ) / b, x > a, b, c > 0)

Moments:

1 Mean µ = a + b 1 +

c 2 1 Variance 2 = b2 1 + 2 1 +

c c where ( x) is Eulers gamma function.

92

Sampling from a 3-parameter Weibull Distribution: X i W (a, b, c)

A random number is drawn from a uniform distribution on the open interval (0,1) and then transformed to a Weibull variate with the Weibull percentile function.

U i U (0,1)

(74)

X i = a + b [ ln(1 U i ) ]

1/ c 2.2.7 Best-Estimate, More Conservative, and Less Conservative Sampling Distributions The sampling distributions shown in Table 6 and Figs. 36-40 were grouped into three categories: (1) best-estimate (BE), (2) more conservative (MC), and (3) less conservative (LC) as shown below. The MC and LC sets of distributions were selected such that, when taken together as a set, that they would be expected to produce either higher or lower LOCA probabilities, respectively, than the BE set of distributions.

Sampling Distributions dR/d da/d R0 flaw-ini Categories (in/yr) (in/month) (in) (months)

BE beta-1 Weibull beta Weibull MC triangluar lognormal triangular lognormal LC Weibull uniform normal triangular In the CCA Case Matrix to be discussed in Sect. 3, additional distributions are included in the study to test the sensitivity of the final CCA results to the applied statistical distributions used to characterize the uncertainties in the analysis. In the ASP study, the BE set of distributions was applied exclusively.

2.2.8 Variance Reduction - The Method of Antithetic Variates A variance reduction technique (VRT) called the Method of Antithetic Variates (AV) [33,42] is applied in the Monte Carlo code developed for this analysis. The objective of this method is to provide a more rapid convergence of the estimated LOCA probabilities to their asymptotic values than the standard Monte Carlo approach, thus producing a more computationally-efficient algorithm. The central idea of the AV method, due originally to Hammersly and Morton [43], is to make pairs of trials (2 trials for each Monte Carlo realization) such that a small observation on one of the trials in a pair tends to be offset by a large observation on the other one; i.e., the two observations are negatively correlated. Taking the average of the results of the two trials as the result for the realization should produce an estimate of the mean (or common expectation) that is closer to the true probability due to the induced negative covariance between the paired trials.

93

For example, at a given point in time, k , let there be a total of N realizations with N (1) f ( k ) from the Trial 1 set and N (2)f ( k ) from the Trial 2 set of realizations experiencing a LOCA at some time, k .

The cumulative probability of failure at time k , is therefore f ( k )

N (1) f ( k )

N (2) f ( k , N ) =

p (1) f ( k , N ) =

p (2)

N N (75) p (1) ( k , N ) + p (2) ( k , N )

p f ( k , N ) =

f f at k for N realizations 2

Equation (75) represents an unbiased point estimator of the statistical expectation, µ LOCA ( k ) =

( )

E p f ( k ) , of the cumulative probability of a LOCA occurring at or before the time k . An estimate for the variance of E p f ( ) for k as a function of the number of Monte Carlo realizations, N, is from

[33, 42]

var p f ( k , N )

var p f ( k ) =

N (76) var p (1) f

( k , N ) + var p(2) f ( k , N ) + 2cov p f ( k , N ) , p f ( k , N )

(1) (2) 4N The individual trial streams (1) and (2) satisfy the requirements of a Bernoulli sequence of discrete random variables such that the estimates provided by Eq. (75) are distributed by a discrete binomial distribution [42]. An unbiased point estimator for the variance of p (fi ) ( k , N ) is, therefore, (1 p (i )

( k , N ) ) p(fi ) ( k , N )

var ( p (fi ) ( k , N ))

f N

for i = 1, 2 (77)

If the two trials within a pair were sampled independently, then the covariance between the two trials would be, cov p (1)( )

f ( k , N ) , p f ( k , N ) = 0 ; however, by inducing a negative covariance between (2) f and p f , the overall variance for the estimate of µ LOCA ( k ) should be reduced.

p (1) (2)

The required negative covariance is applied through the use of complementary random numbers in the sampling protocols required to create the input for each Monte Carlo trial/realization. For the first trial (of 94

the paired trials) for the ith realization (1 i N ) , nine random variates are drawn from a standard uniform distribution using the random number generator, ranf(), such that:

U1((1)i ) U ( 0,1) for sampling R0 (not used in CCA runs) i ) U ( 0,1) for sampling the cavity wastage rate, dR / d (1)

U 2(

i ) U ( 0,1)

(1)

U 3( for lognormal bivariate sampling of ductile tearing model parameters ( J Ic , m )

i ) U ( 0,1)

(1)

U 4(

i ) U ( 0,1) for sampling flaw-ini (not used in CCA runs)

(1)

U 5(

i ) U ( 0,1) for sampling the flaw growth rate, da / d (1)

U 6(

(1)

U 7( (

i ) U ( 0,1) for correlated sampling of the plastic-flow properties yield , ult ;)

i ) U ( 0,1) U 9(i ) used only in the trivariate reduction protocol for the (1) (1)

U 8(

(1)

U 9( i) U ( 0,1 ) bivariate gamma distribution These uniform random variates are then converted into their respective sampling distributions using the Inverse Transform Method as discussed in Sect. 2.2.6. For the second trial of the paired trials in the ith realization (1 i N ) , the complements of each uniform random variate are calculated by:

U1((2) (1) i ) = 1 U1(i ) for sampling R0 (not used in CCA runs)

(2)

U 2( (1) i ) = 1 U 2(i ) for sampling the cavity wastage rate, dR / d (2) (1)

U 3( i ) = 1 U 3(i ) for lognormal bivariate sampling of (2)

= (1) ductile tearing model parameters ( J Ic , m )

U 4( i) 1 U 4(i )

(2)

U 5( (1) i ) = 1 U 5(i ) for sampling flaw-ini (not used in CCA runs)

(2)

U 6( (1) i ) = 1 U 6(i ) for sampling the flaw growth rate, da / d (2)

U 7( (1)

( )

i ) = 1 U 7(i ) for correlated sampling of the plastic-flow properties yield , ult ;

(2) (1) (2)

U 8( i ) = 1 U 8(i ) U 9(i ) used only in the trivariate reduction protocol for the (2)

U 9( (1) i ) = 1 U 9(i )

bivariate gamma distribution with subsequent conversion into their respective sampling distributions. As discussed in [33], synchroni-zation of the random number stream is vital to ensure that the necessary negative correlation is correctly induced for each realization.

95

2.2.9 LOCA Screening Rules Three sets of LOCA screening rules are presented in Fig. 41 in terms of Event Trees that lead from the initiating event to a failure consequence. These screening rules are designated: (a) best-estimate (BE),

(b) more conservative (MC), and (c) less conservative (LC). As was the case with the sampling distributions, the MC and LC LOCA screening rules were developed to produce higher and lower total LOCA probabilities, respectively, than produced by the BE screening rules, thus providing estimated bounding probabilities as a measure of the uncertainties about the predicted best-estimate values.

96

(a)

Fig. 41. Event trees based on LOCA screening rules representing (a) best-estimate (BE),

97

(b)

Fig. 41. (continued) Event trees based on LOCA screening rules representing: (b) more conservative (MC),10 10 Bins 2, 4, 6, and 7 are not used in the MC event-tree definition.

98

(c)

Fig. 41. (continued) and (c) less conservative (LC) rules.

99

2.2.10 LOCA Break Size Definitions Effective break sizes were provided by NRC/RES to define a range of Loss of Coolant Accident (LOCA) characterizations. These effective break sizes are presented in Table 7.

Table 7. LOCA Sizes Defined by Effective Break Sizes LOCA Wastage Effective Cavity Effective Cavity SIZE Area Footprint Diameter Footprint Radius (in2) (in) (in)

SBLOCA 0.1104 to 3.14 0.375 to 2 0.1875 to 1 MBLOCA 3.14 to 72 2 to 9.57 1 to 4.79 LBLOCA > 72 > 9.57 > 4.79 2.2.11 Cavity Growth-Shape Scaling Rules Scaling factors are applied to the sampled cavity wastage rates, dRcavity / d , to approximate different assumptions on how the wastage cavitys unbacked area grows over time. The growth-shape scaling factors are as follows:

1 (1) Rcavity = (designated as LC) to approximate the cavity growing as a semi-ellipse with the major 2

axis growing twice as fast as the minor axis 1

(2) Rcavity = (designated as BE) to approximate the cavity growing as a semi-circle 2

(3) Rcavity = 1 (designated as MC) to approximate the cavity growing as a complete circle where Eq. (55) is now modified to be dRcavity Rcavity x D( dR / d ) (78) d 100

3 Results and Discussion 3.1 Results of Deterministic FEM Analysis The three Davis-Besse Flaw Models, embedded independently as surface flaws on the exposed side of the unbacked cladding, were run with the Davis-Besse wastage-cavity ABAQUS FEM submodel, starting from a zero pressure load up to the point of numerical stability. The resulting driving-force load paths (Japplied vs Pressure) are presented in Fig. 42, where JIc percentiles from the Davis-Besse lognormal model (see Fig. 19) have been overlaid on the plot. At an operating pressure of 2165 psi, the cumulative probabilities that the Model-Flaw driving forces will exceed the ductile-tearing initiation toughness of the Davis-Besse cladding are 9.6 x 1013 % , 3.1 x 1011 % , and 5.0 x 106 % for Model Flaws #1, #2, and

1. 3, respectively. At the relief-valve setpoint pressure of 2500 psi, the probabilities of nonexceedance are 2.6 x 107 % , 7.9 x 106 % , and 1.1 x 101 % for Model Flaws #1, #2, and #3, respectively. Table 8 presents the pressures, predicted by the Davis-Besse FEM submodel, that were required to initiate stable ductile tearing in the three Model Flaws at the 5th, 50th, and 95th JIc percentiles. The margins against ductile-tearing initiation of the Model Flaws at the normal operating pressure of 2165 psi are presented in Table 8a, and the margins at the relief-valve setpoint pressure of 2500 psi are given in Table 8b. For Model Flaw #3 the margin against ductile initiation ranges from 1.2 to 1.5 (with a median of 1.4) at the operating pressure and from 1.08 to 1.3 (with a median of 1.2) at the relief-valve setpoint pressure.

In Fig. 43, the load paths under increasing pressure are presented for Model Flaws 1 (Fig. 43a) and 3 (Fig 43b) in terms of flaw driving forces (Jsubmodel) vs applied pressure. These FEM driving forces are also compared to those predicted by Eq. (44) which was developed from 2 in. long flaws centered in a circular burst disk. It is apparent that flaws in the Davis-Besse wastage-cavity are under significantly different constraint conditions than similar flaws in a circular burst disk. To establish an approximate driving-force similitude between the detailed Davis-Besse submodel (Jsubmodel) and the Monte Carlo model (JMC), a driving-force scaling rule was developed as shown in Fig. 43c. Equation (44) now becomes

( )

J applied a Atip , Rcavity = flaw exp C0 +

C1 C

(

+ 2 + C3 Rcavity + C4 Rcavity ln Rcavity +

a Atip a Atip

)

(

C5 ln Rcavity ) ( )

Rcavity + C6 ln Rcavity 2

} (79) where a Atip flaw = 0.64335 0.71197 tclad 101

Table 8a. DB Submodel Pressures to Initiate Ductile Tearing for the 5th, 50th , and 95th Percentiles with Margins Against Ductile Initiation at the Operating Pressure of 2165 psi Model Flaw Depth Flaw Length Pressures Required for J Ic Probabilities and Margins at Operating Pressure Flaw a A-tip L 5% Margin for 50% Margin for 95% Margin for (in) (in) (psi) Initiation (psi) Initiation (psi) Initiation 1 0.0650 0.660 3963 1.83 5186 2.40 6503 3.00 2 0.0995 0.606 3109 1.44 3553 1.64 4375 2.02 3 0.0995 2.000 2709 1.25 2985 1.38 3306 1.53 Table 8b. DB Submodel Pressures to Initiate Ductile Tearing for the 5th, 50th , and 95th Percentiles with Margins Against Ductile Initiation at the Relief-Valve Setpoint Pressure of 2500 psi Model Flaw Depth Flaw Length Pressures Required for J Ic Probabilities and Margins at Safety Valve Setpoint Flaw a A-tip L 5% Margin for 50% Margin for 95% Margin for (in) (in) (psi) Initiation (psi) Initiation (psi) Initiation 1 0.0650 0.660 3963 1.59 5186 2.07 6503 2.60 2 0.0995 0.606 3109 1.24 3553 1.42 4375 1.75 3 0.0995 2.000 2709 1.08 2985 1.19 3306 1.32 Fig. 42. Load paths (Japplied vs Applied Pressure) for three model surface flaws placed in the Davis-Besse finite-element submodel compared to the cladding materials ductile-tearing initiation fracture toughness with cumulative probabilities from the lognormal model shown in Fig. 19.

102

(a)

(b)

(c)

Fig. 43. A scaling rule is applied to the Driving-Force Model (Eq. (44) to bring the Monte Carlo model driving forces into agreement with the results of the detailed FEM wastage cavity submodel (a) scaling required for Model Flaw 1, (b) scaling required for Model Flaw 3, and (c) linear fit developed between Model Flaws 1 and 3.

103

3.2 Cladding Capacity Analysis (CCA) - Sensitivity Study Results For the Cladding Capacity Analysis (CCA) study, which treats the damage state discovered in Davis Besse on February 16, 2002 as both known and certain, a matrix of 21 cases (see Table 9) was constructed to investigate the sensitivity of the best-estimate case (CCA-001) to many of the assumptions that were applied in the construction of the Davis-Besse Monte Carlo damage-state model. The conditions for each of the CCA cases are presented in Table 9 in which the perturbed conditions (relative to Case CCA-001) are highlighted in red for Cases CCA-002 to CCA-021. The results of all 21 cases are presented in Table 10 for three time stations: (1) the time of discovery (TOD), 16 February 2002, 1 , (2) 6 weeks after the time of discovery, 2 , and (3) 1 year after the time of discovery, 3 . The cumulative probabilities of a LOCA occurring are further partitioned into small-break, medium-break, and large-break LOCAS as defined in Table 7. For the best-estimate Case CCA-001, the calculated probability of a LOCA occurring at 1 is 0%, as required by the fact that at the time of discovery the Davis-Besse wastage cavity had not failed in terms of breaching the integrity of the pressure boundary.

3.2.1 Convergence of Monte Carlo Simulations Figure 44 plots the coefficient of variation, COV, for the total cumulative probability of a LOCA (irrespective of break size) occurring at 3 , one year after the time of discovery. The COV of pf is defined as the ratio of its estimated standard deviation to its estimated mean, where the COV is a function of the number of Monte Carlo trials, N.

(1 p f ( k , N )) p f ( k , N )

var p f ( k , N )

COV p f ( k ) = = N (80) p f ( k , N ) p f ( k , N )

where the cumulative probability of failure is estimated by Eq. (75). The convergence history shown in Fig. 44 indicates that the cumulative LOCA probabilities estimated by the Monte Carlo code are approaching an asymptotic solution after 50,000 realizations (or 100,000 antithetic-paired trials).

104

Table 9. Case Matrix for Cladding Capacity Analysis (CCA)

A-Tip A-Tip As-Found Sampling Scale Sampling Distributions Sampling Cavity LOCA Case Ductile Instability Atip Depth of Plastic Driving dR /d da /d Distribution Growth Binning Number Tearing Tested? (in) Flow Props. Forces (in/yr) (in/month) Group Rules Rules CCA-001 Yes Yes 0.0650 bivariate lognormal Yes beta-1 Weibull BE BE BE CCA-002 Yes Yes 0.0650 bivariate lognormal Yes beta-1 Weibull BE MC BE CCA-003 Yes Yes 0.0650 bivariate lognormal Yes beta-1 Weibull BE LC BE CCA-004 Yes Yes 0.0650 bivariate lognormal Yes triangluar lognormal MC MC MC CCA-005 Yes Yes 0.0650 bivariate lognormal Yes Weibull uniform LC LC LC CCA-006 Yes Yes 0.0650 bivariate lognormal Yes triangluar Weibull mixed BE BE CCA-007 Yes Yes 0.0650 bivariate lognormal Yes Weibull Weibull mixed BE BE CCA-008 Yes Yes 0.0650 bivariate lognormal Yes beta-1 lognormal mixed BE BE CCA-009 Yes Yes 0.0650 bivariate lognormal Yes beta-1 uniform mixed BE BE CCA-010 Yes Yes 0.0650 bivariate lognormal Yes beta-2 Weibull mixed BE BE CCA-011 Yes No 0.0650 bivariate lognormal Yes beta-1 Weibull BE BE BE CCA-012 Yes No 0.0650 bivariate lognormal No beta-1 Weibull BE BE BE CCA-013 No No 0.0650 bivariate lognormal No beta-1 Weibull BE BE BE CCA-014 Yes Yes 0.0650 bivariate lognormal No beta-1 Weibull BE BE BE CCA-015 Yes Yes 0.0650 bivariate lognormal Yes beta-1 Weibull BE BE LC CCA-016 Yes Yes 0.0650 bivariate lognormal Yes beta-1 Weibull BE BE MC CCA-017 Yes Yes 0.0350 bivariate lognormal Yes beta-1 Weibull BE BE BE CCA-018 Yes Yes 0.0995 bivariate lognormal Yes beta-1 Weibull BE BE BE CCA-019 Yes Yes 0.0650 inverse Weibull; r = 1 Yes beta-1 Weibull BE BE BE CCA-020 Yes Yes 0.0650 bivariate gamma Yes beta-1 Weibull BE BE BE CCA-021 Yes Yes 0.0650 inverse Weibull; r = 0 Yes beta-1 Weibull BE BE BE Table. 10. Monte Carlo Results - Summary of LOCA Probabilities (N = 50,000)

LOCA Probabilities (%) LOCA Probabilities (%) LOCA Probabilities (%)

Case Time of Discovery 6 weeks after Time of Discovery 1 year after Time of Discovery Number No LOCA SBLOCA MBLOCA LBLOCA No LOCA SBLOCA MBLOCA LBLOCA No LOCA SBLOCA MBLOCA LBLOCA CCA-001 100.000% 0.000% 0.000% 0.000% 99.757% 0.243% 0.000% 0.000% 29.218% 64.906% 5.188% 0.688%

CCA-002 100.000% 0.000% 0.000% 0.000% 99.562% 0.438% 0.000% 0.000% 24.732% 65.879% 8.220% 1.169%

CCA-003 100.000% 0.000% 0.000% 0.000% 99.809% 0.191% 0.000% 0.000% 33.895% 62.950% 2.832% 0.323%

CCA-004 100.000% 0.000% 0.000% 0.000% 99.288% 0.000% 0.712% 0.000% 16.919% 0.000% 79.181% 3.900%

CCA-005 100.000% 0.000% 0.000% 0.000% 87.919% 12.057% 0.024% 0.000% 5.535% 93.852% 0.613% 0.000%

CCA-006 100.000% 0.000% 0.000% 0.000% 99.693% 0.307% 0.000% 0.000% 20.304% 70.594% 8.054% 1.048%

CCA-007 100.000% 0.000% 0.000% 0.000% 99.732% 0.268% 0.000% 0.000% 29.043% 65.682% 4.700% 0.575%

CCA-008 100.000% 0.000% 0.000% 0.000% 99.679% 0.321% 0.000% 0.000% 46.402% 44.567% 7.803% 1.228%

CCA-009 100.000% 0.000% 0.000% 0.000% 86.892% 13.108% 0.000% 0.000% 4.288% 94.664% 0.944% 0.104%

CCA-010 100.000% 0.000% 0.000% 0.000% 99.835% 0.165% 0.000% 0.000% 39.798% 59.518% 0.659% 0.025%

CCA-011 99.999% 0.001% 0.000% 0.000% 99.685% 0.315% 0.000% 0.000% 28.623% 65.634% 5.100% 0.643%

CCA-012 79.667% 20.333% 0.000% 0.000% 20.309% 79.691% 0.000% 0.000% 0.207% 99.793% 0.000% 0.000%

CCA-013 100.000% 0.000% 0.000% 0.000% 99.763% 0.237% 0.000% 0.000% 29.265% 64.750% 5.172% 0.813%

CCA-014 99.826% 0.174% 0.000% 0.000% 84.032% 15.968% 0.000% 0.000% 6.526% 93.474% 0.000% 0.000%

CCA-015 100.000% 0.000% 0.000% 0.000% 99.757% 0.243% 0.000% 0.000% 29.218% 64.649% 6.133% 0.000%

CCA-016 100.000% 0.000% 0.000% 0.000% 99.757% 0.000% 0.243% 0.000% 29.218% 0.000% 69.830% 0.952%

CCA-017 100.000% 0.000% 0.000% 0.000% 99.964% 0.036% 0.000% 0.000% 40.558% 46.940% 3.533% 8.969%

CCA-018 100.000% 0.000% 0.000% 0.000% 97.923% 2.077% 0.000% 0.000% 17.059% 82.867% 0.074% 0.000%

CCA-019 100.000% 0.000% 0.000% 0.000% 99.753% 0.247% 0.000% 0.000% 29.416% 64.615% 5.165% 0.804%

CCA-020 100.000% 0.000% 0.000% 0.000% 99.759% 0.241% 0.000% 0.000% 29.507% 64.455% 5.172% 0.866%

CCA-021 100.000% 0.000% 0.000% 0.000% 99.760% 0.240% 0.000% 0.000% 29.403% 64.624% 5.189% 0.784%

105

Fig. 44. Convergence of Case 001 as a function of the number of antithetic paired realizations at 1 year after time of discovery.

106

3.2.2 Best-Estimate CCA Results The results of Case CCA-001 present best-estimate cumulative probabilities of a LOCA (see Fig. 45a) occurring over a time period from 1 year before the time of discovery (TOD) up to an elapsed time of 500 days after TOD. These probabilities are intended to answer the question of how long could the Davis-Besse RPV have continued in service if the wastage cavity had not been discovered on 16 February 2002.

From Table 10, the best-estimate case indicates that 6 weeks after 16 February 2002 the estimated cumulative probability of an SBLOCA occurring was approximately 0.24 %. The probabilities of an MBLOCA or LBLOCA occurring for this elapsed time are 0%. After 1 year (on 16 February 2003), the cumulative probability of a failure of the pressure boundary was approximately 70.8 %. These LOCA probabilities have been further binned into 64.9 % for an SBLOCA, 5.2 % for an MBLOCA, and 0.7 %

for an LBLOCA. The SBLOCA probabilities at this time are exclusively Bin 3 (see Fig. 41a) type failures in which, after the A-tip has failed, both the flaws C-tip and the full cavity itself remains stable.

The results from Cases CCA-017 and CCA-018 can be used to provide an approximate 90 % confidence interval covering the best estimate for the median failure of the cladding, as predicted by the results of Case CCA-001. The times to failure after TOD for these cases are presented in Table 11 for 5 %, 50 %,

and 95 % cumulative probabilities, and the failure histories are plotted in Fig. 45b. These results predict that, at a confidence level of 90%, Davis-Besse could have continued operating from 2 to 22 months before cladding failure would be expected.

3.2.3 Results of the CCA Sensitivity Study The sensitivity of the best-estimate Monte Carlo results to variations in the input conditions and statistical distributions are presented in Figs. 46-48.

Figure 46a examines the sensitivity of the overall LOCA probability to the assumed flaw depth at the time of discovery. The sensitivity of the LOCA probability to the assumed cavity growth/shape rule is presented in Fig. 46b. The more conservative cavity growth/shape rule (the cavity grows equally in all directions, i.e., an expanding circular footprint) produces the highest cumulative probabilities of failure.

107

Table 11. Times to Failure after TOD Model Flaw Time to Failure After TOD at Cumulative Probability Depth at TOD 5% 50% 95%

(in) (days)/(months) (days)/(months) (days)/(months) 0.035 103/3.4 304/10 1318/43.3 0.065 79/2.6 232/7.6 982/32.3 0.1 52/1.7 150/4.9 678/22.3 (a)

(b)

Fig. 45. LOCA probability history for: (a) the best-estimate Cladding Capacity Analysis (CCA) case (Case CCA-001) with a further categorization into small-break LOCA (SBLOCA),

medium-break LOCA (MBLOCA), and large-break LOCAs (LBLOCA) and (b) varying model flaw depth at TOD, 0.035 in., 0.065 in., and 0.0995 in.

108

(a) (b)

Fig. 46. Sensitivity of CCA LOCA probabilities to: (a) model flaw depth at TOD (b) cavity growth/shape rules.

(a) (b)

(c) (d)

Fig. 47. Sensitivity of CCA LOCA probabilities to LOCA binning rules (a) total LOCA probability, (b) SBLOCA probability, (c) MBLOCA probability and (d) LBLOCA probability.

109

(a) (b)

(c) (d)

Fig. 48. Sensitivity of CCA LOCA probabilities to modeling of A-tip ductile tearing (a) total LOCA probability, (b) SBLOCA probability, (c) MBLOCA probability and (d) LBLOCA probability.

110

Figure 47a demonstrates, as required, that the LOCA Binning Rules (see Figs. 41a, 41b, and 41c) influence only the partitioning or binning of the LOCA cumulative probabilities into their size categories and not the overall LOCA probability. The sensitivity of the resulting SBLOCAs, MBLOCAs, and LBLOCAs to the binning rules are shown in Figs. 47b, 47c, and 47d, respectively comparing Case CCA-001 to Cases CCA-015 and CCA-016.

The influence of the treatment of A-tip ductile tearing is presented in Figs. 48 comparing case CCA-001 to cases CCA-011 and CCA-013. In Fig. 48b, the inclusion of A-tip ductile tearing in the analysis produces slightly higher probabilities of SBLOCA relative to applying plastic collapse of the remaining ligament. It is important to note that the flaw continues to grow over time, irrespective of any ductile fracture, because of the flaws exposure to the corrosive environment of the wastage cavity; therefore, the increased probabilities of failure over time are due the combined effects of both cavity growth and flaw growth. Ductile tearing of the A-tip tends to cause the cavity to fail slightly earlier in time relative to no ductile tearing, thus increasing the probability of SBLOCAs and reducing the probabilities of MBLOCAs and LBLOCAs.

Additional sensitivities are investigated in the Case Matrix, such as the influence of the statistical distributions assumed for the cavity and flaw growth rates (see Cases CCA-006 to CCA-010) and the influence of the combined Sampling Distribution groups, the cavity growth/shape rules, and the LOCA Binning Rules (see Cases CCA-004 and CCA-005). Cases CCA-001 and CCA-011 indicate that using A-tip ductile-tearing initiation compared to A-tip ductile-tearing instability as the A-tip failure criterion produces only slightly different LOCA probabilities. The much more significant influence of correctly scaling the Monte Carlo model driving forces to establish similitude with the detailed Davis-Besse submodel is demonstrated by comparing Case CCA-001 to Case CCA-012.

3.3 Accident Sequence Precursor Analysis - Best Estimate and Sensitivity Study Results As discussed in Sect. 2.2.5.3, for ASP analyses the damage state at TOD is treated as uncertain, and two additional distributions are sampled in order to provide an estimate for the effective cavity radius, R1 , and model flaw depth, a1 , at TOD. These two distributions are the effective cavity radius at 1 year before TOD, R0 , and the elapsed time since flaw initiation, flawini . Distributions were originally developed based on the results of the Expert Elicitation discussed in Sect. 2.2.6. These distributions were then modified slightly to produce median values for the damage state at TOD that match the assumed fixed values applied in the CCA study. Figures 49 and 50 compare the original derived distributions for R0 and flawini , respectively, to the modified distributions. These distributions are then combined (see Eq. (70))

111

with the sampled values of the cavity wastage rate, dR / d (see Fig. 51), and flaw growth rate, da / d (see Fig. 52), to produce distributions for R1 (Fig. 53) and a1 (see Fig. 54).

(a)

(b)

Fig. 49. Comparison of beta distributions for effective cavity radius at TOD-1 year, 0 , where the original distribution (based on the results of an Expert Elicitation) has been revised to produce a sampled R1(median) closer to the observed valued used in the CCA: (a) probability densities and (b) cumulative probabilities.

112

(a)

(b)

Fig. 50. Comparison of Weibull distributions for the time of flaw initiation (relative to TOD),

flaw-ini(2) (see the time line in Fig. 1), where the original distribution (based on the results of an Expert Elicitation) has been revised to produce a sampled R1(median) closer to the observed valued used in the CCA: (a) probability densities and (b) cumulative probabilities.

113

(a)

(b)

Fig. 51. Best-estimate sampling distribution for cavity wastage rate compared to sampled values:

(a) probability density and (b) cumulative probabilities.

114

(a)

(b)

Fig. 52. Best-estimate sampling distribution for flaw growth rate compared to sampled values:

(a) probability density and (b) cumulative probabilities.

115

(a)

(b)

Fig. 53. Distribution of uncertain effective cavity radius calculated for the time of discovery (TOD):

(a) frequency distribution and (b) cumulative probabilities.

116

(a)

(b)

Fig. 54. Distribution of uncertain effective flaw depth calculated for the time of discovery (TOD):

(a) probability density and (b) cumulative probabilities 117

3.3.1 Case Matrix for ASP Study Table 12 presents the 9 cases developed for the ASP study, where Case ASP-001 represents the best-estimate case. Table 13 gives a layout key that indicates the combination of LOCA screening rules and cavity grow/shape rules used to develop the matrix. Sensitivities to sampling distributions were not addressed in the ASP study; however, they were investigated in the CCA study.

3.3.2 Best-Estimate ASP Results The time histories for the probability of cladding failure (i.e., the probability of a LOCA) are compared in Fig. 55 for the best-estimate cases of the CCA (CCA-001) and the ASP (ASP-001) studies. At TOD, the probability of a LOCA is 0 % for the CCA conditions and 20.1 % for the ASP condition. This difference is due to the combined uncertainties in the cavity wastage rate, the cavity size at TOD-1 (R0), the flaw growth rate, and the time of flaw initiation which affect the ASP results but do not enter into the CCA results. Recall that the damage state at TOD is assumed known (i.e., with no uncertainty) for the CCA study. A breakdown of the LOCA categories in the best-estimate ASP case is shown in Fig. 56, where the dominant SBLOCA probability is due to a combination of Bin 2 and Bin 3 failures using the BE LOCA screening rules.

3.3.3 Results of the ASP Sensitivity Study Figures 57 and Table 14 present the results of the ASP sensitivity study. As demonstrated by the time histories shown in Fig. 57a, the LC and MC LOCA screening rules and the cavity growth/shape rules were developed to provide estimated bounding curves for the BE total LOCA history. When the individual LOCA categories are compared (as in Figs. 57b for SBLOCA and 57c for MBLOCA),

different cases may be required to set the bounding conditions.

118

Table 12. Case Matrix for ASP Sensitivity Study A-Tip A-Tip Sampling Scale Cavity LOCA Case Ductile Instability of Plastic Driving dR/d da/d R0 flaw-ini Sampling Growth Binning Number Tearing Tested? Flow Props. Forces (in/yr) (in/month) (in) (months) Distributions Rules Rules ASP-001 Yes Yes bivariate lognormal Yes beta-1 Weibull beta Weibull BE BE BE ASP-002 Yes Yes bivariate lognormal Yes beta-1 Weibull beta Weibull BE BE LC ASP-003 Yes Yes bivariate lognormal Yes beta-1 Weibull beta Weibull BE BE MC ASP-004 Yes Yes bivariate lognormal Yes beta-1 Weibull beta Weibull BE LC BE ASP-005 Yes Yes bivariate lognormal Yes beta-1 Weibull beta Weibull BE LC LC ASP-006 Yes Yes bivariate lognormal Yes beta-1 Weibull beta Weibull BE LC MC ASP-007 Yes Yes bivariate lognormal Yes beta-1 Weibull beta Weibull BE MC BE ASP-008 Yes Yes bivariate lognormal Yes beta-1 Weibull beta Weibull BE MC LC ASP-009 Yes Yes bivariate lognormal Yes beta-1 Weibull beta Weibull BE MC MC Table 13. Case Matrix Layout Key Cavity LOCA Screening Rules Growth LC BE MC Rule LC ASP-005 ASP-004 ASP-006 BE ASP-002 ASP-001 ASP-003 MC ASP-008 ASP-007 ASP-009 LC - less conservative than best estimate BE - best estimate MC - more conservative than best estimate Table 14. Summary of LOCA Probabilities for ASP Sensitivity Study 1 year before TOD Time of Discovery (TOD) 1 year after TOD Case LOCA Probabilities LOCA Probabilities LOCA Probabilities Number No LOCA LOCA SBLOCA MBLOCA LBLOCA No LOCA LOCA SBLOCA MBLOCA LBLOCA No LOCA LOCA SBLOCA MBLOCA LBLOCA ASP-001 98.6% 1.4% 1.39% 0% 0% 79.9% 20.1% 16.92% 0.50% 2.65% 35.8% 64.2% 44.05% 2.08% 18.06%

ASP-002 99.1% 0.9% 0.21% 0.73% 0% 82.9% 17.1% 7.45% 9.62% 0.05% 41.2% 58.8% 22.49% 32.88% 3.40%

ASP-003 98.6% 1.4% 0.44% 0.94% 0% 80.0% 20.0% 2.93% 14.23% 2.89% 35.8% 64.2% 5.40% 39.44% 19.32%

ASP-004 98.2% 1.8% 1.79% 0.00% 0% 82.2% 17.8% 16.95% 0.24% 0.63% 42.1% 57.9% 47.48% 1.69% 8.78%

ASP-005 98.9% 1.1% 0.25% 0.90% 0% 86.1% 13.9% 6.19% 7.73% 0.00% 50.0% 50.0% 21.69% 27.30% 0.99%

ASP-006 98.2% 1.8% 0.63% 1.14% 0% 82.2% 17.8% 3.84% 13.20% 0.73% 42.1% 57.9% 7.91% 40.41% 9.57%

ASP-007 99.0% 1.0% 1.03% 0% 0% 75.2% 24.8% 15.87% 0.75% 8.13% 29.2% 70.8% 38.94% 2.28% 29.57%

ASP-008 99.3% 0.7% 0.15% 0.58% 0% 77.3% 22.7% 7.96% 13.55% 1.24% 32.3% 67.7% 21.77% 37.83% 8.06%

ASP-009 99.0% 1.0% 0.29% 0.73% 0% 75.2% 24.8% 2.01% 14.18% 8.57% 29.2% 70.8% 3.14% 36.55% 31.11%

119

Fig. 55. Comparison of LOCA probability histories between the CCA study and the ASP study conditions where the deviation is due to uncertainties in cavity wastage rate, cavity size at TOD-1, flaw growth rate, and flaw initiation time. For the CCA study, the damage state at TOD was treated as known with no uncertainty.

120

Fig. 56. LOCA probability history for the best-estimate ASP case (Case ASP-001) with a further categorization into small-break LOCA (SBLOCA), medium-break LOCA (MBLOCA), and large-break LOCAs (LBLOCA).

121

(a)

(b)

Fig. 57. LOCA probability histories for full case matrix: (a) total LOCA probabilities, (b) SBLOCA probabilities, 122

(c)

(d)

Fig. 57. (continued) LOCA probability histories for full case matrix: (c) MBLOCA probabilities, and (d) LBLOCA probabilities.

123

4 Summary and Conclusions This report has presented the results of a PSM analysis of the degraded Davis-Besse RPV head, including a description of the Davis-Besse wastage-area damage model, the technical basis for the model, and the results of a cladding capacity analysis (CCA) and an accident sequence precursor (ASP) analysis of the wastage cavity.

The objectives of CCA and ASP studies were to provide approximate answers to three questions regarding the Davis-Besse event:

(1) What applied pressure would have failed the wastage-cavity cladding at the time of discovery (TOD), 16 February 2002 (CCA)?

(2) How much longer could the Davis-Besse RPV have continued in service without failure of the pressure boundary if the wastage cavity had not been discovered on 16 February 2002 (CCA)?

(3) Including uncertainties in the As-Found damage state of the wastage cavity, what was the probability of failure one year before TOD, and how do these uncertainties affect the estimated probability of failure at TOD (ASP)?

The answer to question No. #1 required the construction of a detailed FEM model of the Davis-Besse wastage cavity which incorporated the results of extensive laboratory measurements and metallographic examinations of the damaged site after it had been removed from the RPV head. The Davis-Besse cladding material was carefully characterized in terms of both strength (plastic-flow properties) and fracture toughness (ductile-tearing initiation and flaw growth). The fracture-toughness characterization (at a service temperature of 600 °F) was carried out by the HSST Program using pre-cracked Charpy V-Notch specimens taken from the Davis-Besse wastage cavity. All of the characterization studies included investigations of the uncertainties in the property measurements. Experimental studies with clad burst disk tests were also carried out by the HSST program at ORNL in parallel with the analytical effort to aid in validating the failure models applied in this analysis.

The results of the deterministic FEM analysis indicate that, for the most conservative assumptions regarding flaw size and depth, the estimated failure pressures all exceed the relief-valve set-point pressure of 2500 psi. This result is in agreement with the forensic finding that exposure of the wastage cavity to the nominal operating pressure of 2165 psi did not produce any evidence of crack initiation. Median pressures needed to fail the wastage cavity (by ductile-tearing initiation of the Model Flaw) were estimated to range between 3000 and 5200 psi, representing a 1.4 to 2.4 margin against the operating pressure. A 90 percent 124

confidence interval covering the median estimates was 2710 psi to 6500 psi (or 1.25 to 3.0 margin against the operating pressure).

Question No. #2 was addressed with a PSM analysis that reflected the very limited state-of-knowledge of how the wastage cavity might be expected to evolve over time, beyond the known damage state at the time of discovery. An Expert Elicitation was carried out by the NRC staff to provide estimates for the wastage-cavity growth rate and the rate of stress-corrosion crack development due to exposure of the unbacked cladding to the concentrated boric acid solution inside the wastage cavity. The resulting best-estimate (Case CCA-001) predicts that the cumulative probability of survival decreases to 50% after approximately 230 days of additional operation. When applying the most conservative flaw model, this median time decreased to approximately 150 days of additional operation (Case CCA-018).

Additionally, the results from Cases CCA-017 and CCA-018 can be used to provide an approximate 90 %

confidence interval covering the best estimate for the median failure of the cladding, as predicted by the results of Case CCA-001. The times to failure after TOD for these cases are presented in Table 11 for 5 %, 50 %, and 95 % cumulative probabilities, and the failure histories are plotted in Fig. 45b. These results predict that, at a confidence level of 90%, Davis-Besse could have continued operating from 2 to 22 months before cladding failure would be expected.

Question No. #3 was addressed in the ASP study wherein additional uncertainty in the damage state at TOD was applied in the analysis. For the TOD, the best-estimate probability of cladding failure increased from 0% for the CCA to approximately 20% for the ASP analysis. At 1 year before TOD, the ASP analysis estimated a low probability of failure of approximately 1%.

125

References

1. S. A. Loehlein, Root Cause Analysis Report, Significant Degradation of Reactor Pressure Vessel Head, CR 2002-0891, Davis-Besse Power Station, April 15, 2002.

2. B. R. Bass, W. J. McAfee, P. T. Williams, S. Yin, R. K. Nanstad, and M. Sokolov, Experimental Program for Investigating the Influence of Cladding Defects on Burst Pressure, ORNL/NRC/LTR-04/13, Oak Ridge National Laboratory, Oak Ridge, TN, May 2004.

3. Risked-Based Inspection - Development of Guidelines, NUREG/GR - 005/CRTD - Vol. 20 - 1, prepared by the Research Task Force on Risk-Based Inspection Guidelines, American Society of Mechanical Engineers, prepared for the U.S. Nuclear Regulatory Commission, 1992.

4. M. G. Morgan and M. Henrion, Uncertainty - A Guide to Dealing with Uncertainty in Quantitative Risk and Policy Analysis, Cambridge University Press, Cambridge, UK, 1990.

5. NRC Bulletin 2001-01, Circumferential Cracking of Reactor Pressure Vessel Head Penetration Nozzles, August 3, 2001.

6. Recent Experience with Degradation of Reactor Pressure Vessel Head, NRC Information Notice 2002-11, United States Nuclear Regulatory Commission, Office of Nuclear Reactor Regulation, Washington, DC, March 12, 2002.

7. W. H. Cullen, Jr., Overview of Reactor Vessel Head Degradation, U. S. Nuclear Regulatory Com-mission internet web site, http://www.nrc.gov/reactors/operating/ops-experience/vessel-head-degradation.html , 2004.

8. William D. Travers, Status Report on Accident Sequence Precursor Program and Related Initiatives, NRC Report SECY-99-289, U.S. Nuclear Regulatory Commission, December 20, 1999.

9. U. S. Nuclear Regulatory Commission, Risk Assessment Review Group Report, NUREG/CR-0400, Washington, D. C., September 1978.

10. R. J. Belles, J. W. Cletcher, D. A. Copinger, B. W. Dolan, J. W. Minarick, and P. D. OReilly, 1994 Accident Sequence Precursor Program Results, Nuclear Safety, 37(1), (1996) 73-83.

11. P. T. Williams and B. R. Bass, Stochastic Failure Model for the Davis-Besse RPV Head, ORNL/NRC/LTR-02/10, Oak Ridge National Laboratory, Oak Ridge, TN, September 2004.

12. P. T. Williams and B. R. Bass, Structural Assessment of a Corrosion-Degraded Reactor Pressure Vessel Head, Transactions of the 17th International Conference on Structural Mechanics in Reactor Technology (SMiRT 17), Paper No. D03-5, held in Prague, Czech Republic, on August 17-22, 2003.

13. P. T. Williams and B. R. Bass, Analysis of the Davis-Besse RPV Head Wastage Area and Cavity, ORNL/NRC/LTR-02/09, Oak Ridge National Laboratory, Oak Ridge, TN, September 2004.

14. ABAQUS/Standard Users Guide, Volume I, v6.2-4, Hibbitt, Karlsson, and Sorensen, Inc., Pawtucket, RI, 2001.

15. J. Chakrabarty and J. M. Alexander, Hydrostatic Bulging of Circular Diaphragms, J. Strain Anal.

5(3), (1970) 155-161.

16. J. W. Hyres, Final Report: Examination of the Reactor Vessel (RV) Head Degradation at Davis-Besse, Report No. 1140-025-02-24, Lynchburg Technology Center for Metallurgical Examinations, BWXT Technologies, Inc., BWXT Services, Inc., Nuclear and Environmental Operations, Lynchburg, VA, June 2003.

17. Solidworks 2001 Plus, Solidworks Corporation, Concord, Massachusetts, 2002.

126

18. Personal communication with Dr. M. T. Kirk, Office of Nuclear Regulatory Research, United States Nuclear Regulatory Commission, Rockville, MD, February, 2004.

19. American Society of Mechanical Engineers Boiler and Pressure Vessel Code,Section XI, Rules for Inservice Inspection of Nuclear Power Plant Components, Division 1, Subsection IWA-3300 - Flaw Characterization, American Society of Mechanical Engineers, New York, 1998.

20. American Society of Mechanical Engineers Boiler and Pressure Vessel Code,Section XI, Rules for Inservice Inspection of Nuclear Power Plant Components, Division 2, Subsection IWB-3600 -

Analytical Evaluation of Flaws, American Society of Mechanical Engineers, New York, 1998.

21. American Society of Mechanical Engineers Boiler and Pressure Vessel Code,Section XI, Rules for Inservice Inspection of Nuclear Power Plant Components, Appendix A, Analysis of Flaws, American Society of Mechanical Engineers, New York, 1998.

22. S. Rahman, N. Ghadiali, D. Paul, and G. Wilkowski, Probabilistic Pipe Fracture Evaluation for Leak-Rate-Detection Applications, NUREG/CR-6004, U. S. Nuclear Regulatory Commission, Washington, D.C, 1995.

23. S. Rahman, Probabilistic Elastic-Plastic Fracture Analysis of Circumferentially Cracked Pipes with Finite-Length Surface Flaws, Nuclear Engineering Design 195, (2000) 239-260.

24. T. W. Anderson and D. A. Darling, A Test of Goodness of Fit, J. Am. Statist. Assoc. 49, (1954) 765-769.

25. M. A. Stephens, EDF Statistics for Goodness of Fit and Some Comparisons, J. Am. Statist. Assoc.

69, (1974) 730-737.

26. K. Pearson, On a Criterion that a Given System of Deviations from the Probable in the Case of a Correlated System of Variables is Such That it can be Reasonably Supposed to Have Arisen in Random Sampling, Phil. Mag. 50(5), (1900) 157-175.

27. A. M. Law, ExpertFit© Users Guide, Averill M. Law & Associates, Tucson, Arizona, May 2002.

28. Kennedy and Gentle, Statistical Computing, Marcel Dekker, NY, (1980) 95.

29. A. R. DiDinato and A. H. Morris, Computation of the Incomplete Gamma Function Ratios and Their Inverse, ACM Trans. Math. Softw. 12, (1986) 377-393.

30. R. M. Jones and K. S. Miller, On the Multivariate Lognormal Distribution, J. Indust. Math. 16, (1966) 63-76.

31. M. E. Johnson and J. S. Ramberg, Transformations of the Multivariate Normal Distribution with Applications to Simulation, Technical Report LA-UR-77-2595, Los Alamos Scientific Laboratories, Los Alamos, New Mexico, 1978.

32. B. C. Arnold, A Note on Multivariate Distributions with Specified Marginals, J. Am. Statist. Assoc.

62, (1967) 1460-1461.

33. A. M. Law and W. D. Kelton, Simulation and Modeling Analysis, 3rd ed., McGraw-Hill Book Co.,

New York, N.Y., 2000.

34. N. E. Dowling, Mechanical Behavior of Materials: Engineering Methods for Deformation, Fracture, and Fatigue, 2nd ed., Prentice Hall, Upper Saddle River, N.J., (1999).

127

35. F. M. Haggag, W. R. Corwin, and R. K. Nanstad, Irradiation Effects on Strength and Toughness of Three-Wire Series-Arc Stainless Steel Weld Overlay Cladding, NUREG/CR-5511 (ORNL/TM-11439), Oak Ridge National Laboratory, February 1990.

36. F. M. Haggag and R. K. Nanstad, Effects of Thermal Aging and Neutron Irradiation on the Mechanical Properties of Three-Wire Stainless Steel Weld Overlay Cladding, NUREG/CR-6363 (ORNL/TM-13047), Oak Ridge National Laboratory, May 1997.

37. A. Saxena, Nonlinear Fracture Mechanics for Engineers, CRC Press LLC, Boca Raton, FL, 1998, pp.

175-179.

38. Standard Test Method for Measurement of Fracture Toughness, ASTM E 1820-1, Annual Book of ASTM Standards 2002, Section Three, Metals Test Methods and Analytical Procedures, Volume 03.01, American Society for Testing and Materials, West Conshohocken, PA (2002).

39. TableCurve 3D Users Manual, SPSS Inc., Chicago, IL, 1997.

40. P. LEcuyer and S. Cote, Implementing a Random Number Package with Splitting Facilities. ACM Transactions on Mathematical Software 17, (1991)98-111.

41. A. R. DiDinato and A. H. Morris, Algorithm 708: Significant Digit Computation of the Incomplete Beta Function Ratios, ACM Trans. Math. Softw. 18, (1993) 360-373.

42. A. Haldar and S. Mahadevan, Probability, Reliability, and Statistical Methods in Engineering Design, John Wiley & Sons, New York, 2000.

43. J. M. Hammersley and K. W. Morton, A New Monte Carlo Technique: Antithetic Variates, Proc.

Camb. Phil. Soc. 52, (1956) 449-475.

128