From c76a9cd4817f3ad144859b6b380fd73d96861719 Mon Sep 17 00:00:00 2001 From: nali88 Date: Mon, 26 Aug 2024 22:56:42 +0200 Subject: [PATCH] LEU-MET-THERM-006 Cases 1-30 Model of all 30 cases of the LEU-MET-THERM-006 benchmark using OpenMC's Python API. The parameters that differentiate between the configurations of LMT6 are written as variables. An automation file is included to run all 30 cases and save the results. --- icsbep/leu-met-therm-006/case_enum.py | 352 ++++++++++++++++ icsbep/leu-met-therm-006/model.py | 559 ++++++++++++++++++++++++++ icsbep/leu-met-therm-006/run_cases.py | 59 +++ 3 files changed, 970 insertions(+) create mode 100644 icsbep/leu-met-therm-006/case_enum.py create mode 100644 icsbep/leu-met-therm-006/model.py create mode 100644 icsbep/leu-met-therm-006/run_cases.py diff --git a/icsbep/leu-met-therm-006/case_enum.py b/icsbep/leu-met-therm-006/case_enum.py new file mode 100644 index 00000000..5c32b898 --- /dev/null +++ b/icsbep/leu-met-therm-006/case_enum.py @@ -0,0 +1,352 @@ +import copy + +############################################# +# Benchmark: LEU-MET-THERM-006 +# Cases: 1-30 +# Model by: Nour Ali +# Date: 26/08/2024 +############################################# + + +w3 = 'w3' +f1 = 'f1' +f2 = 'f2' + +class info: + """ + c_no=case number\n + w_h=critical water height\n + pitch=assembly pitch\n + beams_number_h=number of beams in horizontal direction\n + beams_number_v=number of beams in vertical direction\n + diag: number of diagonal moves from the corner pin to the center pin\n + vert: number of vertical moves from the corner pin to the center pin\n + lattice: lattice type, rectangular or hexagonal\n + additional_plate: True if there is an additional plate in the assembly\n + assembly_dict: dictionary of all possible assembly configurations""" + + START_IDX = 1 + + #the indices are constant throughout for case numbers + #for example, the parameters for case 1 are all at index 0 + c_no = [x for x in range(1,31)] + w_h = [59.55,57.97,52.91,64.32,53.18,61.77,60.82,59.41,54.65,55.50,48.41,63.56, + 56.44,62.35,61.21,56.47,66.54,63.51,63.94,57.07,61.56,60.26,59.60,57.66, + 63.52,60.95,59.25,58.90,65.65,61.20] + pitch = [10.2,10.2,10.2,10.7,10.7,11.2,11.2,11.2,11.2,11.2,11.2,11.7,11.7,12.2, + 12.2,12.2,12.7,12.7,9.7,9.7,10.2,10.2,10.2,10.2,10.7,10.7,10.7,10.7, + 11.2,11.2] + beams_number_h = [7,7,8,6,7,6,6,6,6,6,7,6,6,7,7,7,8,8, + 6,7,6,6,6,6,6,6,6,6,7,7] + beams_number_v= [7]*30 + + diag=[3,3,3,2,3,3,3,3,2,2,3,3,2,3,3,3,3,3] + rec_d=[19]*12 + diag.extend(rec_d) + + vert=[1,2,2,2,1,1,1,1,2,2,1,1,2,1,1,1,1,2] + rec_v=[20]*12 + vert.extend(rec_v) + + lattice=['hex']*18 + rect=['rect']*12 + lattice.extend(rect) + + additional_plate=[True]*30 + additional_plate[6:9]=[False]*3 + additional_plate[10]=False + + assembly_dict={} + + #case 1 + ring1_4 = [w3] + [f1]*3 + [w3]*5 + [f1] + [w3]*5 + [f1] + [w3]*5 +[f1]*3 + ring1_3= [f1] *9 + [w3] + [f1]*8 + ring1_2 = [f1]*12 + ring1_1 = [f1]*6 + ring1_0=[f1] + + hex_array1=[ring1_4, ring1_3, ring1_2, ring1_1, ring1_0] + assembly_dict[1]=copy.deepcopy(hex_array1) + + #case2 + ring2_5=[w3]*30 + ring2_5[12]=ring2_5[18]=f1 + ring2_4=[w3]*24 + ring2_4[3]=ring2_4[9]=ring2_4[10]=ring2_4[11]=ring2_4[-3]=f1 + ring2_4[-9]=ring2_4[-10]=ring2_4[-11]=f1 + ring2_3= [w3] + [f1]*17 + ring2_2 = [f1]*12 + ring2_1 = [f1]*6 + ring2_0=[f1] + hex_array2=[ring2_5, ring2_4, ring2_3, ring2_2, ring2_1, ring2_0] + assembly_dict[2]=copy.deepcopy(hex_array2) + + # Case 3,18 + ring3_5=[w3]*30 + ring3_5[12]=ring3_5[18]=f1 + ring3_5[11]=f2 + ring3_4=[w3]*24 + ring3_4[3]=ring3_4[9]=ring3_4[10]=ring3_4[11]=f1 + ring3_4[-3]=ring3_4[-9]=ring3_4[-10]=ring3_4[-11]=f1 + ring3_4[4:9]=[f2]*5 + ring3_3= [w3] + [f1]*17 + ring3_2 = [f1]*12 + ring3_1 = [f1]*6 + ring3_0=[f1] + hex_array3=[ring3_5, ring3_4, ring3_3, ring3_2, ring3_1, ring3_0] + assembly_dict[3]=copy.deepcopy(hex_array3) + assembly_dict[18]=copy.deepcopy(hex_array3) + + # Case 4,9,10,13 + ring4_4 = [w3]*24 + ring4_4[3]=ring4_4[9]=ring4_4[10]=ring4_4[14]=f1 + ring4_3= [w3]*18 + ring4_3[1:12]=[f1]*11 + ring4_3[16:18]=[f1]*2 + ring4_2 = [f1]*12 + ring4_1 = [f1]*6 + ring4_0=[f1] + hex_array4=[ring4_4, ring4_3, ring4_2, ring4_1, ring4_0] + assembly_dict[4]=copy.deepcopy(hex_array4) + assembly_dict[9]=copy.deepcopy(hex_array4) + assembly_dict[10]=copy.deepcopy(hex_array4) + assembly_dict[13]=copy.deepcopy(hex_array4) + + # Case 5,11,16 + ring5_4 = [w3]*24 + ring5_4[2]=ring5_4[3]=ring5_4[-2]=ring5_4[-3]=ring5_4[9]=ring5_4[15]=f1 + ring5_3=[f1]*18 + ring5_3[9]=w3 + ring5_2 = [f1]*12 + ring5_1 = [f1]*6 + ring5_0=[f1] + hex_array5=[ring5_4, ring5_3, ring5_2, ring5_1, ring5_0] + assembly_dict[5]=copy.deepcopy(hex_array5) + assembly_dict[11]=copy.deepcopy(hex_array5) + assembly_dict[16]=copy.deepcopy(hex_array5) + + # Case 6,7,12 + ring6_4 = [w3]*24 + ring6_4[15]=f1 + ring6_3=[f1]*18 + ring6_3[0]=ring6_3[9]=w3 + ring6_3[3:7]=[w3]*4 + ring6_2 = [f1]*12 + ring6_1 = [f1]*6 + ring6_0=[f1] + hex_array6=[ring6_4, ring6_3, ring6_2, ring6_1, ring6_0] + assembly_dict[6]=copy.deepcopy(hex_array6) + assembly_dict[7]=copy.deepcopy(hex_array6) + assembly_dict[12]=copy.deepcopy(hex_array6) + + # Case 8 + ring8_4=[w3]*24 + ring8_4[15]=ring8_4[21]=f1 + ring8_3=[f1]*18 + ring8_3[0]=ring8_3[9]=w3 + ring8_3[3:7]=[w3]*4 + ring8_2 = [f1]*12 + ring8_1 = [f1]*6 + ring8_0=[f1] + hex_array8=[ring8_4, ring8_3, ring8_2, ring8_1, ring8_0] + assembly_dict[8]=copy.deepcopy(hex_array8) + + # Case 14 + ring14_4 = [w3]*24 + ring14_4[3]=ring14_4[9]=ring14_4[15]=f1 + ring14_3=[f1]*18 + ring14_3[0]=ring14_3[9]=w3 + ring14_2 = [f1]*12 + ring14_1 = [f1]*6 + ring14_0=[f1] + hex_array14=[ring14_4, ring14_3, ring14_2, ring14_1, ring14_0] + assembly_dict[14]=copy.deepcopy(hex_array14) + + # Case 15 + ring15_4 = [w3]*24 + ring15_4[3]=ring15_4[-3]=ring15_4[9]=ring15_4[-9]=f1 + ring15_3=[f1]*18 + ring15_3[0]=ring15_3[9]=w3 + ring15_2 = [f1]*12 + ring15_1 = [f1]*6 + ring15_0=[f1] + hex_array15=[ring15_4, ring15_3, ring15_2, ring15_1, ring15_0] + assembly_dict[15]=copy.deepcopy(hex_array15) + + #Case 17 + ring17_5=[w3]*30 + ring17_5[3]=ring17_5[-3]=f1 + ring17_4=[f1]*24 + ring17_4[4:9]=[f2]*5 + ring17_4[0]=w3 + ring17_4[10:15]=[w3]*5 + ring17_4[16:21]=[w3]*5 + ring17_3=[f1]*18 + ring17_3[9]=w3 + ring17_2 = [f1]*12 + ring17_1 = [f1]*6 + ring17_0=[f1] + hex_array17=[ring17_5, ring17_4, ring17_3, ring17_2, ring17_1, ring17_0] + assembly_dict[17]=copy.deepcopy(hex_array17) + + # Case 19, 24, 28 + rect_array19=[ + [w3,w3,w3,w3,w3,w3,w3,w3], + [w3,f1,f1,f1,f1,f1,f1,w3], + [w3,f1,f1,f1,f1,f1,f1,w3], + [w3,f1,f1,f1,f1,f1,f1,w3], + [w3,f1,f1,f1,f1,f1,f1,w3], + [w3,f1,f1,f1,f1,f1,f1,w3], + [w3,f1,f1,f1,f1,f1,f1,w3], + [w3,f1,f1,f1,f1,f1,f1,w3], + [w3,w3,w3,w3,w3,w3,w3,w3]] + + assembly_dict[19]=copy.deepcopy(rect_array19) + assembly_dict[24]=copy.deepcopy(rect_array19) + assembly_dict[28]=copy.deepcopy(rect_array19) + + # Case 20 + rect_array20=[ + [w3,w3,w3,w3,w3,w3,w3,w3,w3], + [w3,f1,f1,f1,f1,f1,f1,w3,w3], + [w3,f1,f1,f1,f1,f1,f1,w3,w3], + [w3,f1,f1,f1,f1,f1,f1,f1,w3], + [w3,f1,f1,f1,f1,f1,f1,f1,w3], + [w3,f1,f1,f1,f1,f1,f1,f1,w3], + [w3,f1,f1,f1,f1,f1,f1,f1,w3], + [w3,f1,f1,f1,f1,f1,f1,w3,w3], + [w3,w3,w3,w3,w3,w3,w3,w3,w3]] + + assembly_dict[20]=copy.deepcopy(rect_array20) + + # Case 21 + rect_array21=[ + [w3,w3,w3,w3,w3,w3,w3,w3], + [w3,w3,f2,f1,f1,w3,w3,w3], + [w3,f1,f1,f1,f1,f1,f1,w3], + [w3,f1,f1,f1,f1,f1,f1,w3], + [w3,f1,f1,f1,f1,f1,f1,w3], + [w3,f1,f1,f1,f1,f1,f1,w3], + [w3,f1,f1,f1,f1,f1,f1,w3], + [w3,f1,f1,f1,f1,f1,f1,w3], + [w3,w3,w3,w3,w3,w3,w3,w3]] + + assembly_dict[21]=copy.deepcopy(rect_array21) + + #Case 22 + rect_array22=[ + [w3,w3,w3,w3,w3,w3,w3,w3], + [w3,w3,f1,f1,f1,f1,w3,w3], + [w3,f1,f1,f1,f1,f1,f1,w3], + [w3,f1,f1,f1,f1,f1,f1,w3], + [w3,f1,f1,f1,f1,f1,f1,w3], + [w3,f1,f1,f1,f1,f1,f1,w3], + [w3,f1,f1,f1,f1,f1,f1,w3], + [w3,f1,f1,f1,f1,f1,w3,w3], + [w3,w3,w3,w3,w3,w3,w3,w3]] + + assembly_dict[22]=copy.deepcopy(rect_array22) + + #Case23 + rect_array23=[ + [w3,w3,w3,w3,w3,w3,w3,w3], + [w3,f1,f1,f1,f1,f1,w3,w3], + [w3,f1,f1,f1,f1,f1,f1,w3], + [w3,f1,f1,f1,f1,f1,f1,w3], + [w3,f1,f1,f1,f1,f1,f1,w3], + [w3,f1,f1,f1,f1,f1,f1,w3], + [w3,f1,f1,f1,f1,f1,f1,w3], + [w3,f1,f1,f1,f1,f1,w3,w3], + [w3,w3,w3,w3,w3,w3,w3,w3]] + + assembly_dict[23]=copy.deepcopy(rect_array23) + + #Case 25 + rect_array25=[ + [w3,w3,w3,w3,w3,w3,w3,w3], + [w3,w3,f1,f1,f1,w3,w3,w3], + [w3,f1,f1,f1,f1,f1,f1,w3], + [w3,f1,f1,f1,f1,f1,f1,w3], + [w3,f1,f1,f1,f1,f1,f1,w3], + [w3,f1,f1,f1,f1,f1,f1,w3], + [w3,f1,f1,f1,f1,f1,f1,w3], + [w3,f1,f1,f1,f1,f1,f1,w3], + [w3,w3,w3,w3,w3,w3,w3,w3]] + + assembly_dict[25]=copy.deepcopy(rect_array25) + + # Case 26 + rect_array26=[ + [w3,w3,w3,w3,w3,w3,w3,w3], + [w3,w3,f1,f1,f1,f1,w3,w3], + [w3,f1,f1,f1,f1,f1,f1,w3], + [w3,f1,f1,f1,f1,f1,f1,w3], + [w3,f1,f1,f1,f1,f1,f1,w3], + [w3,f1,f1,f1,f1,f1,f1,w3], + [w3,f1,f1,f1,f1,f1,f1,w3], + [w3,f1,f1,f1,f1,f1,f1,w3], + [w3,w3,w3,w3,w3,w3,w3,w3]] + + assembly_dict[26]=copy.deepcopy(rect_array26) + + #case27 + rect_array27=[ + [w3,w3,w3,w3,w3,w3,w3,w3], + [w3,f1,f1,f1,f1,f1,f1,w3], + [w3,f1,f1,f1,f1,f1,f1,w3], + [w3,f1,f1,f1,f1,f1,f1,w3], + [w3,f1,f1,f1,f1,f1,f2,w3], + [w3,f1,f1,f1,f1,f1,f1,w3], + [w3,f1,f1,f1,f1,f1,f1,w3], + [w3,f1,f1,f1,f1,f1,f1,w3], + [w3,w3,w3,w3,w3,w3,w3,w3]] + + assembly_dict[27]=copy.deepcopy(rect_array27) + + #case 29 + rect_array29=[ + [w3,w3,w3,w3,w3,w3,w3,w3,w3], + [w3,f1,f1,f1,f1,f1,f1,f1,w3], + [w3,f1,f1,f1,f1,f1,f1,w3,w3], + [w3,f1,f1,f1,f1,f1,f1,f1,w3], + [w3,f1,f1,f1,f1,f1,f1,f1,w3], + [w3,f1,f1,f1,f1,f1,f1,w3,w3], + [w3,f1,f1,f1,f1,f1,f1,w3,w3], + [w3,f1,f1,f1,f1,f1,f1,f1,w3], + [w3,w3,w3,w3,w3,w3,w3,w3,w3]] + + assembly_dict[29]=copy.deepcopy(rect_array29) + + #case 30 + rect_array30=[ + [w3,w3,w3,w3,w3,w3,w3,w3,w3], + [w3,f1,f1,f1,f1,f1,f1,f1,w3], + [w3,f1,f1,f1,f1,f1,f1,w3,w3], + [w3,f1,f1,f1,f1,f1,f1,f1,w3], + [w3,f1,f1,f1,f1,f1,f1,f1,w3], + [w3,f1,f1,f1,f1,f1,f1,f2,w3], + [w3,f1,f1,f1,f1,f1,f1,w3,w3], + [w3,f1,f1,f1,f1,f1,f1,f1,w3], + [w3,w3,w3,w3,w3,w3,w3,w3,w3]] + + assembly_dict[30]=copy.deepcopy(rect_array30) + + +class case_enum: + """ + This class assigns the parameters of each case to its corresponding case number. + """ + def __init__(self, case_num) -> None: + idx = case_num - info.START_IDX + self.c_no = info.c_no[idx] + self.w_h = info.w_h[idx] + self.pitch=info.pitch[idx] + self.beams_number_h=info.beams_number_h[idx] + self.beams_number_v=info.beams_number_v[idx] + self.diag=info.diag[idx] + self.vert=info.vert[idx] + self.lattice=info.lattice[idx] + self.additional_plate=info.additional_plate[idx] + self.assembly_lat=copy.deepcopy(info.assembly_dict[case_num]) + + \ No newline at end of file diff --git a/icsbep/leu-met-therm-006/model.py b/icsbep/leu-met-therm-006/model.py new file mode 100644 index 00000000..d16e77b6 --- /dev/null +++ b/icsbep/leu-met-therm-006/model.py @@ -0,0 +1,559 @@ +import openmc +import numpy as np +import math +import os +from case_enum import case_enum + +############################################# +# Benchmark: LEU-MET-THERM-006 +# Cases: 1-30 +# Model by: Nour Ali +# Date: 26/08/2024 +############################################# + + +class model: + + #left corner of the assembly + xcorner=-40 + + def __init__(self, case_info:case_enum, cross_section:str, w_unc=0, temperature=296, + dense_unc= 0, relative_path= './case', batch=300, inactive=50, shannon=False, + Tally=True, run=True) -> None: + + """ + case_info: case number given as case_enum object, EX: case_enum(1). \n + cross_section: path to the cross section library used, EX: '/path/to/file/endfb-viii.0-hdf5/cross_sections.xml'.\n + w_unc: uncertainty in critical water height.\n + temperature: temperature of the model in kelvin.\n + dense_unc: uncertainty in fuel density.\n + relative_path: relative path to the directory where the xml files will be saved and case will be run. \n + batch: number of batches.\n + inactive: number of inactive batches.\n + shannon: Bool, True to calculate the shannon entropy.\n + Tally: Bool, True to calculate the tally.\n + run: Bool, True to run the case.\n + """ + + self.case_info = case_info + self.temp=temperature + self.dense_unc=dense_unc + self.r_path=relative_path + self.cross_section=cross_section + self.batch=batch + self.inactive=inactive + self.w_unc=w_unc + + self.setup() + + if run: + if Tally: + self.Get_Tally() + if shannon: + self.shannon_entropy() + self.export_xml() + + self.run_model() + + def setup(self): + openmc.config['cross_sections']= self.cross_section + self._init_global() + self.material() + self.fuel_cell() + self.U_shaped_plates() + self.beams() + self.additional_plate() + self.assembly() + self.outer_box() + self.root_geo() + self.source_def() + self.export_xml() + self.setup_done=True + + def _init_global(self): + + # case Variables + if self.case_info.additional_plate: + self.w_height= openmc.ZPlane(0.4+self.case_info.w_h+self.w_unc) + else: + self.w_height= openmc.ZPlane(self.case_info.w_h + self.w_unc) + + self.dx=self.case_info.pitch*math.cos(math.radians(30)) + self.dy=abs(self.case_info.pitch*math.sin(math.radians(30))) + + if self.case_info.additional_plate: + self.z04=openmc.ZPlane(0.4) + self.z59=openmc.ZPlane(59 + 0.4) + self.z595=openmc.ZPlane(59 + 0.5 +0.4) + self.z505=openmc.ZPlane(50 + 0.5 +0.4) + self.z50=openmc.ZPlane(50.4) + else: + self.z04=openmc.ZPlane(0) + self.z59=openmc.ZPlane(59) + self.z595=openmc.ZPlane(59 + 0.5) + self.z505=openmc.ZPlane(50 + 0.5) + self.z50=openmc.ZPlane(50) + + self.direc=self.r_path+str(self.case_info.c_no) + self.abspath = os.path.abspath(self.direc) + + if not os.path.exists(self.abspath): + os.mkdir(self.abspath) + + + ## Materials + def fuel_mat(self): + fuel=openmc.Material(temperature=self.temp) + fuel.set_density('g/cc', 18.465 + self.dense_unc) + fuel.add_nuclide('C0',9.2396E-04) #different cross section libraries accept C differently + # fuel.add_nuclide('C12',9.137E-04) + # fuel.add_nuclide('C13',1.024E-05) + fuel.add_nuclide('Al27',2.8722E-04) + + Fe_ao=[5.47E-05,1.2633E-06,1.6812E-07] + for i in range(0,3): + fuel.add_nuclide('Fe'+str(56+i), Fe_ao[i]) + + fuel.add_nuclide('Fe54', 3.4846E-06) + + U_ao=[7.0181E-06,7.5547E-04] + for i in range(0,2): + fuel.add_nuclide('U'+str(234+i), U_ao[i]) + + fuel.add_nuclide('U238', 4.5866E-02) + self.fuel = fuel + + def water_mat(self): + self.water=openmc.Material(temperature=self.temp) + self.water.set_density('g/cc', 0.99754) + self.water.add_nuclide('H1',6.6682e-2) + self.water.add_nuclide('H2',7.6693e-6) + self.water.add_nuclide('O16',3.3265e-2) + self.water.add_nuclide('O17', 1.2671e-5) + self.water.add_s_alpha_beta('c_H_in_H2O') + + + def air_mat(self): + + air=openmc.Material(temperature=self.temp) + air.set_density('g/cc',0.0012) + air.add_nuclide('O16',1.0568E-05) + air.add_nuclide('O17',4.0258E-09) + air.add_nuclide('N14',3.9349E-05) + air.add_nuclide('N15',1.4375E-07 ) + self.air=air + + def steel_mat(self): + steel=openmc.Material(temperature=self.temp) + steel.set_density('g/cc', 7.9017) + steel.add_nuclide('C0',5.5547e-4) + # steel.add_nuclide('C12',5.493E-04) + # steel.add_nuclide('C13',6.155E-06) + steel.add_nuclide('Fe54', 4.9735e-3) + fe_sao=[7.8073e-2,1.8031e-3,2.3995e-4] + for i in range(0,3): + steel.add_nuclide('Fe'+str(56+i), fe_sao[i]) + + self.steel=steel + + def steelw_mat(self): + steelw=openmc.Material(temperature=self.temp) + steelw.set_density('g/cc', 6.4697) + steelw.add_nuclide('C0', 1.0438E-04) + # steelw.add_nuclide('C12',1.0322E-04) + # steelw.add_nuclide('C13',1.1565E-06) + nucl=['H1', 'H2', 'O16', 'O17', 'Si28', 'Si29', 'Si30', 'P31', 'S32', 'S33', 'S34', + 'S36', 'Cr50', 'Cr52','Cr53','Cr54', 'Mn55', 'Fe54', 'Fe56', 'Fe57', 'Fe58', + 'Ni58','Ni60', 'Ni61','Ni62', 'Ni64'] + aden=[1.0338E-02,1.1890E-06,5.1569E-03,1.9644E-06,1.3723E-03,6.9714E-05,4.6010E-05, + 5.3972E-05,3.7137E-05,2.9322E-07,1.6616E-06,3.9096E-09,6.2863E-04,1.2123E-02, + 1.3746E-03,3.4217E-04,1.5213E-03,3.0137E-03,4.7309E-02,1.0926E-03,1.4540E-04, + 4.8473E-04,1.8672E-04,8.1165E-06,2.5879E-05,6.5906E-06] + for i in range(0,len(nucl)): + steelw.add_nuclide(nucl[i],aden[i]) + + steelw.add_s_alpha_beta('c_H_in_H2O') + self.steelw=steelw + + def material(self): + self.fuel_mat() + self.water_mat() + self.air_mat() + self.steel_mat() + self.steelw_mat() + self.materials_file = openmc.Materials([self.fuel, self.water, + self.air, self.steel,self.steelw]) + self.colors={self.fuel: 'red', self.water:'blue', + self.air: 'green', self.steel:'grey', self.steelw:'lime'} + + + ## Geometry + def fuel_cell(self): + + ## Surfaces + ## Fuel + ## Type 1 + inner_r1=openmc.ZCylinder(r=3.75) + outer_r1=openmc.ZCylinder(r=4.85) + + #Type2 + inner_r2=openmc.ZCylinder(r=3.85) + outer_r2=openmc.ZCylinder(r=4.75) + + #Disks + disk_r=openmc.ZCylinder(r=4.85) + + ## Type 1 + fuel1c=openmc.Cell(name='f1c', fill=self.fuel, + region=-outer_r1 & +inner_r1 & +self.z04 & -self.z59) + disk1c=openmc.Cell(fill=self.steel, + region=-disk_r & +self.z59 & -self.z595) + #watrr inside fuel pin + wa1=openmc.Cell(name='wain1', fill=self.water, + region= -inner_r1 & +self.z04 & -self.w_height & -self.z59) + #air inside fuel pin + fair=openmc.Cell(fill=self.air, region=-inner_r1 & +self.w_height & -self.z59) + + pinwa=openmc.Cell(fill=self.water, + region= +outer_r1 & +self.z04 & -self.w_height & -self.z595) + pinair=openmc.Cell(fill=self.air, region= +outer_r1 & +self.w_height & -self.z595) + + self.f1=openmc.Universe(name='f1',cells=[fuel1c, disk1c, pinwa, pinair, fair, wa1]) + + ## Type 2 + fuel2c=openmc.Cell(fill=self.fuel, + region=-outer_r2 & +inner_r2 & +self.z04 & -self.z50) + disk2c=openmc.Cell(fill=self.steel, + region=-outer_r2 & +self.z50 & -self.z505) + #water inside fuel pin + wa2=openmc.Cell(fill=self.water, + region= -inner_r2 & +self.z04 & -self.w_height & -self.z50) + wa23=openmc.Cell(fill=self.water, + region= -outer_r2 & +self.z505 & -self.w_height & -self.z595) + #air inside fuel pin + fair2=openmc.Cell(fill=self.air, region=-inner_r2 & +self.w_height &-self.z50) + fair3=openmc.Cell(fill=self.air, + region=-outer_r2 & +self.w_height & +self.z505 &-self.z595) + + pinwa2=openmc.Cell(fill=self.water, + region= +outer_r2 & +self.z04 & -self.w_height & -self.z595) + pinair2=openmc.Cell(fill=self.air, + region= +outer_r2 & +self.w_height &-self.z595) + self.f2=openmc.Universe(name='f2',cells=[fuel2c,disk2c,pinwa2,pinair2, fair2, fair3, wa2, wa23]) + + ## Water pin + pinw=openmc.Cell(fill=self.water, + region= +self.z04 & -self.w_height & -self.z595) + pina=openmc.Cell(fill=self.air, region= +self.w_height & -self.z595) + self.w3=openmc.Universe(name='w3',cells=[pinw,pina]) + + def U_shaped_plates(self): + + ## U Plates + + top_plate=openmc.model.RectangularParallelepiped(xmin=0.3, xmax=9.7, + ymin=0, ymax=100, + zmin=-0.2, zmax=0) + + left_pl=openmc.model.RectangularParallelepiped(xmin=0.3,xmax=0.5, + ymin=0, ymax=100, + zmin=-4, zmax=0) + + right_pl=openmc.model.RectangularParallelepiped(xmin=9.5, xmax=9.7, + ymin=0, ymax=100, + zmin=-4, zmax=0) + + whole_plate= -top_plate | -left_pl | -right_pl + + + plate_bound= openmc.model.RectangularParallelepiped(xmin=0, xmax=(9.7+(0.6/2)), + ymin=0, ymax=100, + zmin=-4, zmax=0) + + plate_water = +top_plate | +left_pl | +right_pl + + uplate=openmc.Cell(fill=self.steel, region=whole_plate) + uplate_w=openmc.Cell(fill= self.water, region= plate_water & -plate_bound) + + self.uplate0=openmc.Universe(cells=[uplate, uplate_w]) + self.uplate0c=openmc.Cell(fill=self.uplate0, region=-plate_bound) + + uplate_region_dict={} + + for i in range(0,5): + uplate_region_dict['uplate'+str(i+1)]=openmc.model.RectangularParallelepiped( + xmin=(i+1)*(9.4+0.6), xmax=(9.7+(0.6/2))+((i+1)*(9.4+0.6)), + ymin=0, ymax=100, + zmin=-4, zmax=0) + uplate_region_dict['uplaten'+str(i+1)]=openmc.model.RectangularParallelepiped( + xmin=-((i+1)*(9.4+0.6)), xmax=-((i+1)*(9.4+0.6))+(9.4+0.6), + ymin=0, ymax=100, + zmin=-4, zmax=0) + + + self.uplate_dict={} + for i in range(0,5): + self.uplate_dict['uplate'+str(i+1)]=openmc.Cell(fill=self.uplate0, region=-uplate_region_dict['uplate'+str(i+1)]) + self.uplate_dict['uplate'+str(i+1)].translation=[(i+1)*(9.4+0.6),0,0] + self.uplate_dict['uplaten'+str(i+1)]=openmc.Cell(fill=self.uplate0,region=-uplate_region_dict['uplaten'+str(i+1)]) + self.uplate_dict['uplaten'+str(i+1)].translation=[-(i+1)*(9.4+0.6),0,0] + + + uplate_region_dict['uplaten'+str(5+1)]=openmc.model.RectangularParallelepiped( + xmin=-((5+1)*(9.4+0.6)), xmax=-((5+1)*(9.4+0.6))+(9.4+0.6), + ymin=0, ymax=100, + zmin=-4, zmax=0) + + self.uplaten6=openmc.Cell(fill=self.uplate0, region=-uplate_region_dict['uplaten'+str(5+1)]) + self.uplaten6.translation=[-60,0,0] + + def beams(self): + + if self.case_info.additional_plate: + zmin_outer=60.9 + zmax_outer=64.9 + zmin_inner=61.1 + zmax_inner=64.7 + else: + zmin_outer=60.5 + zmax_outer=64.5 + zmin_inner=60.7 + zmax_inner=64.3 + + + ## Beams + beam_outer = openmc.model.RectangularParallelepiped(xmin= model.xcorner-2, xmax=model.xcorner+2, + ymin= 0, ymax=100, + zmin=zmin_outer, zmax=zmax_outer) + beam_inner=openmc.model.RectangularParallelepiped(xmin=model.xcorner-2+0.2, xmax=model.xcorner+2-0.2, + ymin=0, ymax=100, + zmin= zmin_inner, zmax=zmax_inner) + whole_beam= +beam_inner & -beam_outer + beam_cell=openmc.model.RectangularParallelepiped(xmin= model.xcorner-(self.case_info.pitch/2), + xmax=model.xcorner+(self.case_info.pitch/2), + ymin=0, ymax=100, + zmin=zmin_outer, zmax=zmax_outer) + + beam=openmc.Cell(fill=self.steel, region=whole_beam) + beam_cell_wa=openmc.Cell(fill=self.water, region= -beam_cell & +beam_outer & -self.w_height) + beam_wa_in=openmc.Cell(fill=self.water, region= -beam_inner & -self.w_height) + + beam_cell_air=openmc.Cell(fill=self.air, region= -beam_cell & +beam_outer & +self.w_height) + beam_air_in=openmc.Cell(fill=self.air, region= -beam_inner & +self.w_height) + + self.b1=openmc.Universe(cells=[beam, beam_cell_wa, beam_cell_air, beam_air_in, beam_wa_in]) + + beam_region_dict={} + self.beam_dict={} + + if self.case_info.lattice == 'rect': + for i in range (0,self.case_info.beams_number_h): + beam_region_dict['beam'+str(i)]=openmc.model.RectangularParallelepiped( + xmin=model.xcorner-(self.case_info.pitch/2) + ((i)*self.case_info.pitch), + xmax = model.xcorner+(self.case_info.pitch/2) + ((i)*self.case_info.pitch), + ymin=0, ymax=100, + zmin=zmin_outer, zmax=zmax_outer + ) + self.beam_dict['beam'+str(i)]=openmc.Cell(fill=self.b1, region=-beam_region_dict['beam'+str(i)]) + self.beam_dict['beam'+str(i)].translation=[(((i)*self.case_info.pitch)),0,0] + + if self.case_info.lattice == 'hex': + for i in range (0,self.case_info.beams_number_h): + beam_region_dict['beam'+str(i)]=openmc.model.RectangularParallelepiped( + xmin=model.xcorner-(self.case_info.pitch/2) + ((i)*self.dx), xmax = model.xcorner+(self.case_info.pitch/2) + ((i)*self.dx), + ymin=0, ymax=100, + zmin=zmin_outer, zmax=zmax_outer + ) + self.beam_dict['beam'+str(i)]=openmc.Cell(fill=self.b1, region=-beam_region_dict['beam'+str(i)]) + self.beam_dict['beam'+str(i)].translation=[(((i)*self.dx)),0,0] + + + def additional_plate(self): + + adplate=openmc.model.RectangularParallelepiped(xmin=-59.7, xmax=59.7, + ymin=0, ymax=100, + zmin=0, zmax=0.4) + if self.case_info.additional_plate: + self.addplate=openmc.Cell(fill=self.steelw, region=-adplate) + + + def assembly(self): + + assembly_lat=self.case_info.assembly_lat + + translation = {'f1':self.f1, 'w3':self.w3, 'f2':self.f2} + for i in range(0,len(assembly_lat)): + for j in range(0, len(assembly_lat[i])): + assembly_lat[i][j] = translation[assembly_lat[i][j]] + + + if self.case_info.lattice == 'rect': + rect_lat=openmc.RectLattice(name='Case30') + + rect_lat.pitch = (self.case_info.pitch, self.case_info.pitch) + rect_lat.universes= assembly_lat + + + rect_lat.lower_left=(model.xcorner-(self.case_info.pitch/2)- + self.case_info.pitch,10-(self.case_info.pitch/2)-self.case_info.pitch) + + + self.assembly_length_h=(self.case_info.pitch*(self.case_info.beams_number_h+1)) + self.assembly_length_v=(self.case_info.pitch*(self.case_info.beams_number_v+1)) + + assembly_bound=openmc.model.rectangular_prism(width=self.assembly_length_h, + height=self.assembly_length_v, + origin=(model.xcorner+(3*self.case_info.pitch),(3*self.case_info.pitch)+10), + boundary_type='transmission') + + self.assembly_cell=openmc.Cell(fill=rect_lat, region= assembly_bound & +self.z04 & -self.z595) + + if self.case_info.lattice == 'hex': + dx1=self.case_info.diag*self.dx + dy1=self.case_info.diag*self.dy + (self.case_info.vert*self.case_info.pitch) + + c_x=model.xcorner + dx1 + c_y=10 + dy1 + + hex_lat = openmc.HexLattice() + hex_lat.center = (c_x, c_y) + hex_lat.pitch = (self.case_info.pitch,) + hex_lat.outer = self.w3 + hex_lat.universes = assembly_lat + + self.dx_bound=2*dx1+3*self.case_info.pitch + self.dy_bound=2*dy1+4*self.case_info.pitch + + lat_bound=openmc.model.rectangular_prism(width=self.dx_bound, height=self.dy_bound, + origin=(c_x, c_y), boundary_type='transmission') + + + self.assembly_cell=openmc.Cell(fill=hex_lat, region=lat_bound & +self.z04 & -self.z595) + + def outer_box(self): + + self.external=openmc.model.RectangularParallelepiped(xmin= -65, xmax=65, + ymin= -20, ymax= 120, + zmin=-10, zmax= 130) + outw=openmc.Cell(fill= self.water, region= -self.external & -self.w_height) + outa=openmc.Cell(fill=self.air, region= -self.external & +self.w_height) + outbox=openmc.Universe(cells=[outw,outa]) + + ##Zeroth Universe + box_cell=openmc.Cell(fill=outbox, + region= -self.external & ~self.assembly_cell.region + & ~self.uplate0c.region + & ~self.uplaten6.region) + + if self.case_info.additional_plate: + box_cell.region &= ~self.addplate.region + + for i in range(0,5): + box_cell.region &= ~self.uplate_dict['uplate'+str(i+1)].region + box_cell.region &= ~self.uplate_dict['uplaten'+str(i+1)].region + + for i in range(0, self.case_info.beams_number_h): + box_cell.region &= ~self.beam_dict['beam' +str(i)].region + + U1=openmc.Universe(cells=[self.assembly_cell, box_cell, self.uplate0c, self.uplaten6]) + + if self.case_info.additional_plate: + U1.add_cell(self.addplate) + + for i in range(0,5): + U1.add_cell(self.uplate_dict['uplate'+str(i+1)]) + U1.add_cell(self.uplate_dict['uplaten'+str(i+1)]) + + for i in range(0,self.case_info.beams_number_h): + U1.add_cell(self.beam_dict['beam'+str(i)]) + self.U1_cell=openmc.Cell(fill=U1, region= -self.external) + + def root_geo(self): + self.outer_box() + + void=openmc.model.RectangularParallelepiped(xmin= -500, xmax=500, + ymin= -500, ymax= 500, + zmin=-500, zmax= 500, boundary_type='reflective') + + + outer_void=openmc.Cell(region=+self.external & -void) + U0=openmc.Universe(cells= [outer_void, self.U1_cell]) + + self.geometry = openmc.Geometry(U0) + + def source_def(self): + + if self.case_info.lattice == 'hex': + box_src = openmc.stats.Box([model.xcorner-self.case_info.pitch, + 10-self.case_info.pitch, 0], + [model.xcorner-self.case_info.pitch+self.dx_bound, + 10-self.case_info.pitch+self.dy_bound, 60], + only_fissionable=True) + src=openmc.IndependentSource(space=box_src) + self.settings = openmc.Settings() + self.settings.source = src + self.settings.batches = self.batch + self.settings.inactive = self.inactive + self.settings.particles = 100000 + + if self.case_info.lattice == 'rect': + box_src=openmc.stats.Box([model.xcorner-(self.case_info.pitch/2)-self.case_info.pitch, + 10-(self.case_info.pitch/2)-self.case_info.pitch,0.4], + [model.xcorner-(self.case_info.pitch/2)-self.case_info.pitch + self.assembly_length_h, + 10-(self.case_info.pitch/2)-self.case_info.pitch + self.assembly_length_v,59], + only_fissionable=True) + src=openmc.IndependentSource(box_src) + self.settings = openmc.Settings() + self.settings.source = src + self.settings.batches = 300 + self.settings.inactive = 50 + self.settings.particles = 100000 + + + def shannon_entropy(self): + if self.case_info.lattice == 'rect': + entropy_mesh = openmc.RegularMesh() + entropy_mesh.lower_left = (model.xcorner-(self.case_info.pitch/2)-self.case_info.pitch, + 10-(self.case_info.pitch/2)-self.case_info.pitch,0.4) + entropy_mesh.upper_right = (model.xcorner-(self.case_info.pitch/2)-self.case_info.pitch + self.assembly_length_h, + 10-(self.case_info.pitch/2)-self.case_info.pitch + self.assembly_length_v, + 60) + entropy_mesh.dimension = (70, 70, 2) + self.settings.entropy_mesh = entropy_mesh + + if self.case_info.lattice == 'hex': + entropy_mesh = openmc.RegularMesh() + entropy_mesh.lower_left = (model.xcorner-self.case_info.pitch,10-self.case_info.pitch,0) + entropy_mesh.upper_right = (model.xcorner-self.case_info.pitch+self.dx_bound, + 10-self.case_info.pitch+self.dy_bound, + 61) + entropy_mesh.dimension = (70, 70, 2) + self.settings.entropy_mesh = entropy_mesh + + + def Get_Tally(self): + + tally=openmc.Tally() + tally.scores = ['fission'] + energy_filter = openmc.EnergyFilter([0.0, 0.025, 0.625, 20e6]) + tally.filters = [energy_filter] + self.tallies = openmc.Tallies([tally]) + self.tallies.export_to_xml(self.direc+'/tallies.xml') + + def export_xml(self): + + self.materials_file.export_to_xml(self.direc+'/materials.xml') + self.geometry.export_to_xml(self.direc+'/geometry.xml') + self.settings.export_to_xml(self.direc+'/settings.xml') + + def run_model(self): + if self.setup_done: + openmc.run(cwd=self.abspath, path_input=self.abspath) + + def statepoint(self): + state=openmc.StatePoint(self.direc + f'/statepoint.{self.settings.batches}.h5') + return state + + + + diff --git a/icsbep/leu-met-therm-006/run_cases.py b/icsbep/leu-met-therm-006/run_cases.py new file mode 100644 index 00000000..1e347a3b --- /dev/null +++ b/icsbep/leu-met-therm-006/run_cases.py @@ -0,0 +1,59 @@ +from case_enum import case_enum +from model import model +import os +import json + +############################################# +# Benchmark: LEU-MET-THERM-006 +# Cases: 1-30 +# Model by: Nour Ali +# Date: 26/08/2024 +############################################# + + +def verify(cross_section:str, path='../verify/'): + """ + This function is used to run the 30 cases and save the k-effective in a txt and json files.\n + This function will first check if a case has been run before, + if the case has been run, it will not run it again\n + path: relative path to save case directory. """ + + + case={} + + case_path=path+'case' + + #check if k_eff.json exists + if os.path.exists(path+'k_eff.json'): + with open(path+'k_eff.json') as k: + k_eff=json.load(k) + else: + k_eff={} + + + for i in range(0,30): + + #check if k_eff['case'+str(i+1)] exists and has a vlue + if 'case'+str(i+1) in k_eff.keys(): + if k_eff['case'+str(i+1)] != None: + continue + + #run cases + case[i+1]=model(case_enum(i+1), relative_path=case_path, cross_section=cross_section) + + #save k_eff + sp=case[i+1].statepoint() + k_eff['case'+str(i+1)]=str(sp.k_combined) + + with open(path+'k_eff.txt','w') as f: + for key, value in k_eff.items(): + f.write(f'{key}: {value}\n') + + #save k_eff to json file + with open(path+"k_eff.json", "w") as outfile: + json.dump(k_eff, outfile) + + + + +