Light-Water Reactor Pressure Vessel Model: 3D Probabilistic Fracture Mechanics Model

Contact: Ben Spencer, benjamin.spencer@inl.gov

*Model was co-developed by Ben Spencer and Will Hoffman

Model link: Grizzly RPV PFM Model

High Level Summary of Model

Probabilistic fracture mechanics (PFM) analyis assesses the probability of fracture at the location of any flaw in the population of flaws introduced in the RPV during the manufacturing process. PFM relies on the results of a thermomechanical model of the RPV that does not consider flaws to provide local stress and temperature conditions for each of the flaws considered.

This problem demonstrates the probabilistic fracture mechanics analyis of an RPV, modeled in 3D, containing a population of flaws. It uses the results of the 3D Thermomechanical Model for this analysis.

Computational Model Description

The model evaluates the conditional probability of fracture initiation (CPI) for a population of flaws in an RPV, using a 3D model to represent the global thermomechanical response. In addition to any 3D effects that may be present in the thermomechanical response (which are minimal in this case because spatially uniform cooling is applied), the model also includes 3D aspects of the RPV configuration, including the layout of the various subregions that make up the RPV and the spatial distribution of fluence.

The layout of the subregions making up the RPV in this example is shown in Figure 1. The actual RPV is larger than the region represented here, but this subregion layout emcompasses the beltline region adjacent to the reactor core, which experiences the highest neutron flux and is typically the focus of RPV integrity analyses.

Figure 1: Layout of plate and weld regions in the RPV modeled here

Files used by this model include:

  • MOOSE input file

  • CSV (comma-separated value) file defining the properties of the subregions that make up the RPV

  • Data files defining the distributions of flaw densities

  • CSV files generated as output from the global model described in 3D Thermomechanical Model, located in that model's directory

This document reviews the basic elements of the input file, listed in full here:

# ==============================================================================
# Probabilistic fracture model for LWR Reactor Pressure Vessel
# Application : Grizzly
# ------------------------------------------------------------------------------
# Idaho Falls, INL, 2024
# Author(s): Ben Spencer, Will Hoffman
# If using or referring to this model, please cite as explained on
# https://mooseframework.inl.gov/virtual_test_bed/citing.html
# ==============================================================================

[Mesh]
  type = GeneratedMesh
  dim = 1
  nx = 1
[]

[Samplers]
  [sample]
    type = RPVFractureSampler
    num_rpvs = 50000
    subregions_file = plates_and_welds.csv
    surface_vflaw_file = surface_open_access.dat
    plate_vflaw_file = plate_open_access.dat
    weld_vflaw_file = weld_open_access.dat
    vessel_geometry = vessel_geom
    length_unit = M
    plate_sigma_cu = 0.0073
    plate_sigma_ni = 0.0244
    plate_sigma_p = 0.0013
    weld_sigma_cu = 0.167
    weld_sigma_ni = 0.0165
    weld_sigma_p = 0.0013
    ni_addition_weld = false
    global_fluence_coefficient_of_variation = 0.118
    local_fluence_coefficient_of_variation = 0.056
    execute_on = initial
  []
[]

[UserObjects]
  [vessel_geom]
    type = VesselGeometry
    total_thickness = 0.219202
    clad_thickness = 0.004064
    inner_radius = 2.1971
  []
  [flaw_data]
    type = FlawDataFromSampler
    sampler = sample
    vessel_geometry = vessel_geom
    execute_on = initial
  []
  [clad_axial_coefs]
    type = PolynomialCoefficientsFromFile
    polynomial_coefficient_file = ../thermomechanical/rpv_thermomechanical_3d_out_coefs_axial_clad.csv
    coord_system = 3D_CARTESIAN
    symmetry = QUARTER
    flaw_data = flaw_data
    axial_start = -0.2
    axial_end = 4.2
    axial_num_points = 16
    azimuthal_start = 0.0
    azimuthal_end = 90.0
    azimuthal_num_points = 12
  []
  [clad_hoop_coefs]
    type = PolynomialCoefficientsFromFile
    polynomial_coefficient_file = ../thermomechanical/rpv_thermomechanical_3d_out_coefs_hoop_clad.csv
    coord_system = 3D_CARTESIAN
    symmetry = QUARTER
    flaw_data = flaw_data
    axial_start = -0.2
    axial_end = 4.2
    axial_num_points = 16
    azimuthal_start = 0.0
    azimuthal_end = 90.0
    azimuthal_num_points = 12
  []
  [base_axial_coefs]
    type = PolynomialCoefficientsFromFile
    polynomial_coefficient_file = ../thermomechanical/rpv_thermomechanical_3d_out_coefs_axial_base.csv
    coord_system = 3D_CARTESIAN
    symmetry = QUARTER
    flaw_data = flaw_data
    axial_start = -0.2
    axial_end = 4.2
    axial_num_points = 16
    azimuthal_start = 0.0
    azimuthal_end = 90.0
    azimuthal_num_points = 12
  []
  [base_hoop_coefs]
    type = PolynomialCoefficientsFromFile
    polynomial_coefficient_file = ../thermomechanical/rpv_thermomechanical_3d_out_coefs_hoop_base.csv
    coord_system = 3D_CARTESIAN
    symmetry = QUARTER
    flaw_data = flaw_data
    axial_start = -0.2
    axial_end = 4.2
    axial_num_points = 16
    azimuthal_start = 0.0
    azimuthal_end = 90.0
    azimuthal_num_points = 12
  []
  [base_temp_coefs]
    type = PolynomialCoefficientsFromFile
    polynomial_coefficient_file = ../thermomechanical/rpv_thermomechanical_3d_out_coefs_temp_base.csv
    coord_system = 3D_CARTESIAN
    symmetry = QUARTER
    flaw_data = flaw_data
    axial_start = -0.2
    axial_end = 4.2
    axial_num_points = 16
    azimuthal_start = 0.0
    azimuthal_end = 90.0
    azimuthal_num_points = 12
  []
  [temperature]
    type = FieldValueFromCoefficients
    base_coefficient_calculator = base_temp_coefs
    vessel_geometry = vessel_geom
  []
  [fluence]
    type = FluenceAttenuatedFromSurface
    length_unit = M
    flaw_data = flaw_data
  []
  [embrittlement]
    type = EmbrittlementEONY
    irradiation_time = 1.2614e+9
    irradiation_temperature = 547.0
    fluence_calculator = fluence
    plate_type = CE
    weld_type = OTHER
    flaw_data = flaw_data
    version = FAVOR16_EASON_2006
  []
  [ki_calculator]
    type = KIAxisAlignedROM
    base_axial_coefficient_calculator = base_axial_coefs
    base_hoop_coefficient_calculator = base_hoop_coefs
    clad_axial_coefficient_calculator = clad_axial_coefs
    clad_hoop_coefficient_calculator = clad_hoop_coefs
    flaw_data = flaw_data
    vessel_geometry = vessel_geom
    length_unit = M
    axial_surface_sific_method = FAVOR16
    circumferential_surface_sific_method = FAVOR16
    embedded_sific_method = FAVOR16
  []
  [cpi_calculator]
    type = FractureProbability
    ki_calculator = ki_calculator
    temperature_calculator = temperature
    embrittlement_calculator = embrittlement
    ki_unit = PA_SQRT_M
    temperature_unit = K
    flaw_data = flaw_data
  []
[]

[Executioner]
  type = Transient
  start_time = 0.0
  dt = 60.0
  end_time = 2000.0
[]

[Problem]
  solve = false
[]

[VectorPostprocessors]
  active = 'rpv_failprob cpi_running_statistics flaw_failure_data'
  [rpv_failprob]
    type = RPVFailureProbability
    cpi_calculator = cpi_calculator
    flaw_data = flaw_data
    execute_on = timestep_end
    outputs = none
  []
  [cpi_running_statistics]
    type = VectorPostprocessorRunningStatistics
    vector_postprocessor = rpv_failprob
    vector_name = failure_probabilities
    execute_on = timestep_end
    outputs = final
  []
  [flaw_failure_data]
    type = RPVFlawFailureData
    cpi_calculator = cpi_calculator
    vessel_geometry = vessel_geom
    execute_on = timestep_end
    outputs = final
  []
  [samples]
    type = RPVSamplerData
    sampler = sample
    execute_on = initial
    outputs = final
  []
[]

[Postprocessors]
  [rpv_failprob_mean]
    type = VectorPostprocessorStatistics
    vector_postprocessor = rpv_failprob
    vector_name = failure_probabilities
    quantity = Mean
    execute_on = timestep_end
    outputs = 'console out'
  []
  [rpv_failprob_std_dev]
    type = VectorPostprocessorStatistics
    vector_postprocessor = rpv_failprob
    vector_name = failure_probabilities
    quantity = StandardDeviation
    execute_on = timestep_end
    outputs = 'console out'
  []
[]

[Outputs]
  [out]
    type = CSV
    execute_on = 'timestep_end'
  []
  [final]
    type = CSV
    execute_on = 'final'
  []
  [initial]
    type = CSV
    execute_on = 'initial'
  []
[]
(lwr/rpv_fracture/probabilistic_fracture/rpv_pfm_3d.i)

Input File

Mesh

This model does not perform a finite element simulation, because that was done in the previous global thermomechanical simulation step. However, all MOOSE models require a mesh, even if it is not used, so a simple single-element mesh is defined here.

[Mesh]
  type = GeneratedMesh
  dim = 1
  nx = 1
[]
(lwr/rpv_fracture/probabilistic_fracture/rpv_pfm_3d.i)

Samplers

This block is used to define the way in which the parameters defining each flaw in the population are defined. The RPVFractureSampler is a specialized Sampler developed specifically for RPV analysis.

A key parameter in this block is num_rpvs, which is set to a relatively low number here to make the model easy to run quickly. Larger numbers will give better-converged solutions, but will require more computational resources.

This sampler uses distributions of the flaw population defined in files using the VFLAW format Simonen et al. (2004). Separate files are used for the surface-breaking flaws, plates, and weld regions. Each of these files, which are text files, contain 1000 separate blocks of data defining distributions of flaws. While the flaw distribution blocks in general differ from each other, in this example these files are taken directly from the printouts in Simonen et al. (2004), which provide only the first distribution. These same distributions are then repeated 1000 times to match the required file formatting. Because they are based on public data, these files are named plate_open_access.dat, surface_open_access.dat, and weld_open_access.dat for the distributions of embedded flaws in plates, surface-breaking flaws, and embedded flaws in welds, respectively.

This block also requires that the vessel geometry and its units be defined, and requires the values for a number of parameters defining how concentrations of alloying elements are sampled. The file defining the properties of subregions shown in Figure 1, which is plates_and_welds.csv is specified here.

[Samplers]
  [sample]
    type = RPVFractureSampler
    num_rpvs = 50000
    subregions_file = plates_and_welds.csv
    surface_vflaw_file = surface_open_access.dat
    plate_vflaw_file = plate_open_access.dat
    weld_vflaw_file = weld_open_access.dat
    vessel_geometry = vessel_geom
    length_unit = M
    plate_sigma_cu = 0.0073
    plate_sigma_ni = 0.0244
    plate_sigma_p = 0.0013
    weld_sigma_cu = 0.167
    weld_sigma_ni = 0.0165
    weld_sigma_p = 0.0013
    ni_addition_weld = false
    global_fluence_coefficient_of_variation = 0.118
    local_fluence_coefficient_of_variation = 0.056
    execute_on = initial
  []
[]
(lwr/rpv_fracture/probabilistic_fracture/rpv_pfm_3d.i)

UserObjects

Grizzly uses a modular structure for its PFM calculations Spencer et al. (2019), and utilizes a set of code objects that are specializations of each of the objects shown in Figure 2. Each of these objects is defined in the UserObjects section of the input. The set of objects used here is typical of those that would be used in a PFM analysis, but in many cases, there are options to substitute in other types of objects to perform various parts of the computation in different ways.

Figure 2: Dependencies of UserObjects used in PFM calculations

The objects used here are summarized below:

  • The VesselGeometry object defines the important dimensions of the RPV.

  • The FlawDataFromSampler object provides the data for the flaws to be evaluated from the RPVFractureSampler that was described above. Other types of FlawData objects can provide flaw information from other sources, such as CSV files.

  • The PolynomialCoefficientsFromFile objects provide the polynomial coefficients that are used for computing the various stress components and temperatures in the clad and base blocks from the CSV files generated in the global thermomechanical analysis.

  • The FieldValueFromCoefficients object computes the temperature at the flaw locations from one of the polynomial coefficient objects described above.

  • The FluenceAttenuatedFromSurface object computes a fluence attenuated from a value defined at the RPV inner surface. Other options are also available for defining the fluence, including from an Exodus file, which permits accurate transfer of data from neutron transport calculations.

  • The EmbrittlementEONY object computes the embrittlement at specific flaw locations using the EONY model.

  • The KIAxisAlignedROM object computes mode-II stress intensity factor using an influence coefficient method.

  • The FractureProbability object computes the fracture probability for individual flaws using the KIK_I, temperature, and embrittlement computed by the other objects described above.

[UserObjects]
  [vessel_geom]
    type = VesselGeometry
    total_thickness = 0.219202
    clad_thickness = 0.004064
    inner_radius = 2.1971
  []
  [flaw_data]
    type = FlawDataFromSampler
    sampler = sample
    vessel_geometry = vessel_geom
    execute_on = initial
  []
  [clad_axial_coefs]
    type = PolynomialCoefficientsFromFile
    polynomial_coefficient_file = ../thermomechanical/rpv_thermomechanical_3d_out_coefs_axial_clad.csv
    coord_system = 3D_CARTESIAN
    symmetry = QUARTER
    flaw_data = flaw_data
    axial_start = -0.2
    axial_end = 4.2
    axial_num_points = 16
    azimuthal_start = 0.0
    azimuthal_end = 90.0
    azimuthal_num_points = 12
  []
  [clad_hoop_coefs]
    type = PolynomialCoefficientsFromFile
    polynomial_coefficient_file = ../thermomechanical/rpv_thermomechanical_3d_out_coefs_hoop_clad.csv
    coord_system = 3D_CARTESIAN
    symmetry = QUARTER
    flaw_data = flaw_data
    axial_start = -0.2
    axial_end = 4.2
    axial_num_points = 16
    azimuthal_start = 0.0
    azimuthal_end = 90.0
    azimuthal_num_points = 12
  []
  [base_axial_coefs]
    type = PolynomialCoefficientsFromFile
    polynomial_coefficient_file = ../thermomechanical/rpv_thermomechanical_3d_out_coefs_axial_base.csv
    coord_system = 3D_CARTESIAN
    symmetry = QUARTER
    flaw_data = flaw_data
    axial_start = -0.2
    axial_end = 4.2
    axial_num_points = 16
    azimuthal_start = 0.0
    azimuthal_end = 90.0
    azimuthal_num_points = 12
  []
  [base_hoop_coefs]
    type = PolynomialCoefficientsFromFile
    polynomial_coefficient_file = ../thermomechanical/rpv_thermomechanical_3d_out_coefs_hoop_base.csv
    coord_system = 3D_CARTESIAN
    symmetry = QUARTER
    flaw_data = flaw_data
    axial_start = -0.2
    axial_end = 4.2
    axial_num_points = 16
    azimuthal_start = 0.0
    azimuthal_end = 90.0
    azimuthal_num_points = 12
  []
  [base_temp_coefs]
    type = PolynomialCoefficientsFromFile
    polynomial_coefficient_file = ../thermomechanical/rpv_thermomechanical_3d_out_coefs_temp_base.csv
    coord_system = 3D_CARTESIAN
    symmetry = QUARTER
    flaw_data = flaw_data
    axial_start = -0.2
    axial_end = 4.2
    axial_num_points = 16
    azimuthal_start = 0.0
    azimuthal_end = 90.0
    azimuthal_num_points = 12
  []
  [temperature]
    type = FieldValueFromCoefficients
    base_coefficient_calculator = base_temp_coefs
    vessel_geometry = vessel_geom
  []
  [fluence]
    type = FluenceAttenuatedFromSurface
    length_unit = M
    flaw_data = flaw_data
  []
  [embrittlement]
    type = EmbrittlementEONY
    irradiation_time = 1.2614e+9
    irradiation_temperature = 547.0
    fluence_calculator = fluence
    plate_type = CE
    weld_type = OTHER
    flaw_data = flaw_data
    version = FAVOR16_EASON_2006
  []
  [ki_calculator]
    type = KIAxisAlignedROM
    base_axial_coefficient_calculator = base_axial_coefs
    base_hoop_coefficient_calculator = base_hoop_coefs
    clad_axial_coefficient_calculator = clad_axial_coefs
    clad_hoop_coefficient_calculator = clad_hoop_coefs
    flaw_data = flaw_data
    vessel_geometry = vessel_geom
    length_unit = M
    axial_surface_sific_method = FAVOR16
    circumferential_surface_sific_method = FAVOR16
    embedded_sific_method = FAVOR16
  []
  [cpi_calculator]
    type = FractureProbability
    ki_calculator = ki_calculator
    temperature_calculator = temperature
    embrittlement_calculator = embrittlement
    ki_unit = PA_SQRT_M
    temperature_unit = K
    flaw_data = flaw_data
  []
[]
(lwr/rpv_fracture/probabilistic_fracture/rpv_pfm_3d.i)

Executioner

The Executioner block simply defines the time stepping in this case, because no equation solution is performed. The start and end times must fall within those for the global thermomechanical model. If the time steps do not match those in the global thermomechanical model, they will be interpolated, but it is generally recommended for the timestepping in the PFM analysis to match that of the thermomechanical solution.

[Executioner]
  type = Transient
  start_time = 0.0
  dt = 60.0
  end_time = 2000.0
[]
(lwr/rpv_fracture/probabilistic_fracture/rpv_pfm_3d.i)

Problem

The solve=false parameter set here is used to instruct MOOSE to not solve an equation system. MOOSE still does, however, run the specified objects even when solve=false is set.

[Problem]
  solve = false
[]
(lwr/rpv_fracture/probabilistic_fracture/rpv_pfm_3d.i)

VectorPostprocessors

The objects defined here are used to compute vectors of data useful for output. These include:

  • RPVFailureProbability, which computes CPI for each sampled RPV based on the CPI for the individual flaws within that RPV.

  • VectorPostprocessorRunningStatistics, which computes the running statistics for the mean value and standard deviation of CPI over the course of the Monte Carlo iterations.

  • RPVFlawFailureData, which gathers key data about each flaw with nonzero CPI, including its location, CPI, and KIK_I and temperature at the point in time when it reached its maximum CPI.

  • RPVSamplerData, which makes the sampled variables for each flaw available for output. Note that this block is not included in the set of active VectorPostprocessors because the resulting output can be very large. This is useful for debugging purposes.

[VectorPostprocessors]
  active = 'rpv_failprob cpi_running_statistics flaw_failure_data'
  [rpv_failprob]
    type = RPVFailureProbability
    cpi_calculator = cpi_calculator
    flaw_data = flaw_data
    execute_on = timestep_end
    outputs = none
  []
  [cpi_running_statistics]
    type = VectorPostprocessorRunningStatistics
    vector_postprocessor = rpv_failprob
    vector_name = failure_probabilities
    execute_on = timestep_end
    outputs = final
  []
  [flaw_failure_data]
    type = RPVFlawFailureData
    cpi_calculator = cpi_calculator
    vessel_geometry = vessel_geom
    execute_on = timestep_end
    outputs = final
  []
  [samples]
    type = RPVSamplerData
    sampler = sample
    execute_on = initial
    outputs = final
  []
[]
(lwr/rpv_fracture/probabilistic_fracture/rpv_pfm_3d.i)

Postprocessors

The objects defined here provide scalar values of key model outcomes at each time step. These include:

[VectorPostprocessors]
  active = 'rpv_failprob cpi_running_statistics flaw_failure_data'
  [rpv_failprob]
    type = RPVFailureProbability
    cpi_calculator = cpi_calculator
    flaw_data = flaw_data
    execute_on = timestep_end
    outputs = none
  []
  [cpi_running_statistics]
    type = VectorPostprocessorRunningStatistics
    vector_postprocessor = rpv_failprob
    vector_name = failure_probabilities
    execute_on = timestep_end
    outputs = final
  []
  [flaw_failure_data]
    type = RPVFlawFailureData
    cpi_calculator = cpi_calculator
    vessel_geometry = vessel_geom
    execute_on = timestep_end
    outputs = final
  []
  [samples]
    type = RPVSamplerData
    sampler = sample
    execute_on = initial
    outputs = final
  []
[]
(lwr/rpv_fracture/probabilistic_fracture/rpv_pfm_3d.i)

Outputs

The blocks defined here are used to define outputs that are done at each time step, and at the initial and final times. The only quantities available for output from a PFM analysis are Postprocessors and VectorPostprocessors, both of which are commonly written to CSV files, so it is essential to define appropriate output blocks for this data. Some of the other blocks refer to these outputs with their outputs parameter.

[Outputs]
  [out]
    type = CSV
    execute_on = 'timestep_end'
  []
  [final]
    type = CSV
    execute_on = 'final'
  []
  [initial]
    type = CSV
    execute_on = 'initial'
  []
[]
(lwr/rpv_fracture/probabilistic_fracture/rpv_pfm_3d.i)

Running the model

To run this model using the Grizzly executable, run the following command:

mpiexec -n 4 /path/to/app/grizzly-opt -i rpv_pfm_3d.i

This runs the model using 4 processors. It is recommended that the model be run using as many processors as are available on the machine to accelerate completion of the simulation because these simulations can require significant computation time, especially as the number of RPV realizations is increased.

This model generates the following output files:

  • A CSV file with the Monte Carlo iteration convergence history: rpv_pfm_3d_final_cpi_running_statistics_0035.csv

  • A CSV file with data on the flaws with nonzero CPI: rpv_pfm_3d_final_flaw_failure_data_0035.csv

  • A CSV file with time histories of all postprocessors: rpv_pfm_3d_out.csv

Plotting scripts are provided to process the data provided in these files. The cpi_convergence.py script can be run to generate the convergence history plot shown in Figure 3. It is evident that a larger number of RPV samples would need to be evaluated to obtain a fully converged solution. Depending on the model, this can require on the order of 100,000 RPV samples.

Figure 3: Convergence history of CPI for this example

The scatter_spatial.py script can be run to generate the scatter plots showing the physical locations of all sampled flaws with nonzero CPI shown in Figure 4 and Figure 5.

Figure 4: Spatial distribution of flaws with nonzero CPI plotted in terms of flaw depth vs. azimuthal position

Figure 5: Spatial distribution of flaws with nonzero CPI plotted in terms of flaw axial vs. azimuthal position

It can also be helpful to plot the values of KIK_I vs. relative temperature (difference between the current temperature and reference nil-ductility temperature) at the point when the maximum CPI is reached. The plot generated in Figure 6 was generated using the scatter_temp_ki.py script in this directory.

Figure 6: Scatter plot of the KIK_I at the time when CPI is a maximum plotted vs. the difference between the current temperature and the reference nil-ductility temperature

References

  1. F. A. Simonen, S. R. Doctor, G. J. Schuster, and P. G. Heasler. A generalized procedure for generating flaw-related inputs for the FAVOR code. Technical Report NUREG/CR-6817, PNNL-14268, Pacific Northwest National Laboratory, Richland, WA, March 2004. URL: http://www.nrc.gov/docs/ML0408/ML040830499.pdf.[BibTeX]
  2. B.W. Spencer, W.M. Hoffman, and M.A. Backman. Modular system for probabilistic fracture mechanics analysis of embrittled reactor pressure vessels in the Grizzly code. Nuclear Engineering and Design, 341:25–37, January 2019. doi:10.1016/j.nucengdes.2018.10.015.[BibTeX]