diff options
| -rwxr-xr-x | hp_height_study_driver.sh | 28 | ||||
| -rwxr-xr-x | krusty.py | 106 |
2 files changed, 98 insertions, 36 deletions
diff --git a/hp_height_study_driver.sh b/hp_height_study_driver.sh new file mode 100755 index 0000000..234bbaf --- /dev/null +++ b/hp_height_study_driver.sh @@ -0,0 +1,28 @@ +#!/bin/bash + +mkdir -p results/outputs +touch results/keff_table +echo "\#" $(date) > results/keff_table + + +if [[ $1 == "lcl" ]]; then + prog=openmc-inst01.py krusty.py +else + prog=./krusty.py +fi + +green="\e[0;32m" +reset="\e[0m" + +for pool_height in $(seq 1 45); do + echo -e "$green Starting runs for $pool_height cm pool height... $reset" + for sat_height in $(seq $pool_height 45); do + output_file=results/outputs/p$pool_height\_s$sat_height.out + $prog -p $pool_height -s $sat_height > $output_file + keff=$(cat $output_file | grep "Combined k-effective" | cut -d "=" -f 2 | cut -d "+" -f 1) + pm=$(cat $output_file | grep "Combined k-effective" | cut -d "=" -f 2 | cut -d "-" -f 2) + md5=$(cat $output_file | grep "MD5:" | cut -d ":" -f 2 | sed "s/ //g") + echo K-eff $keff +/- $pm found for p=$pool_height and s=$sat_height \($md5\) + echo $pool_height $sat_height $keff $pm $md5 \($(date)\) >> results/keff_table + done +done @@ -15,10 +15,20 @@ cli_parser = argparse.ArgumentParser( description="Program to investigate the reactivity effect of heat pipe saturation on the KRUSTY reactor." ) -#cli_parser.add_argument("-t","--technique",default="none",choices=["none","vwh","rpt","rrpt"],help="What homogenization technique to use (if any) [none/vwh/rpt/rrpt], default none") -cli_parser.add_argument("-n","--na-height",default="0.8",type=float,help="Radius for RPT or RRPT homogenization [cm], default 0.8 cm. In the case of RRPT, this is the inner radius.") +cli_parser.add_argument("-s","--sat-height",default=5.5,type=float,help="Height of wick that is saturated with Na starting at the bottom of the heat pipe (0-45) [cm], default 25.0") +cli_parser.add_argument("-p","--pool-height",default=5.0,type=float,help="Height of Na pool starting at the bottom of the heat pipe (0-45) [cm], default 5.0") +cli_parser.add_argument("-l","--lclrs",default=False,action="store_true",help="Flag for specifying LCLRS accessible cross-sections, default false.") cli_args = cli_parser.parse_args() +if cli_args.pool_height < 0.0 or cli_args.pool_height > 100.0: + sys.exit(f"Error: Invalid sodium pool height of {cli_args.pool_height} cm. Value must be between 0.0 and 100.0 inclusive.") + +if cli_args.sat_height < 0.0 or cli_args.sat_height > 100.0: + sys.exit(f"Error: Invalid saturated wick height of {cli_args.sat_height} cm. Value must be between 0.0 and 100.0 inclusive.") + +if cli_args.lclrs: + openmc.config["cross_sections"] = "/libs/endfb81_official/cross_sections.xml" + ############ INTPUT-OUTPUT MAPPING ############ path = sys.argv[0] @@ -37,6 +47,7 @@ print("="*60) mat_list = [] mat_temp_global = 900 mat_temp_fuel = 900 +mat_hp_wick_porosity = 0.7 ## 100 - Structural Materials ## mat_haynes230 = openmc.Material(100,name='Structural HAYNES 230') @@ -114,13 +125,37 @@ mat_fuel.temperature = mat_temp_fuel mat_fuel.set_density('atom/b-cm', 4.96130E-02) mat_list.append(mat_fuel) -mat_hp_na = openmc.Material(201,name="Heat Pipe Sodium Coolant") -mat_hp_na.add_element("Na",1) -mat_hp_na.temperature = mat_temp_global -mat_hp_na.set_density("g/cc",0.763) -mat_list.append(mat_hp_na) +mat_hp_na_vap = openmc.Material(201,name="Heat Pipe Core Sodium (Vapour)") +mat_hp_na_vap.add_element("Na",1) +mat_hp_na_vap.temperature = mat_temp_global +mat_hp_na_vap.set_density("kg/m3",1.7e-2) +mat_list.append(mat_hp_na_vap) + +mat_hp_na_liq = openmc.Material(202,name="Heat Pipe Core Sodium (Liquid)") +mat_hp_na_liq.add_element("Na",1) +mat_hp_na_liq.temperature = mat_temp_global +mat_hp_na_liq.set_density("kg/m3",805) +mat_list.append(mat_hp_na_liq) + +mat_hp_wick_unsat = openmc.Material.mix_materials( + material_id=203, + materials=[mat_ss316,mat_hp_na_vap], + fracs=[1-mat_hp_wick_porosity, mat_hp_wick_porosity], + percent_type="vo", + name="Heat Pipe Wick (Unsaturated)" + ) +mat_list.append(mat_hp_wick_unsat) + +mat_hp_wick_sat = openmc.Material.mix_materials( + material_id=204, + materials=[mat_ss316,mat_hp_na_liq], + fracs=[1-mat_hp_wick_porosity, mat_hp_wick_porosity], + percent_type="vo", + name="Heat Pipe Wick (Saturated)" + ) +mat_list.append(mat_hp_wick_sat) -mat_beo = openmc.Material(202,name="BeO Reflector") +mat_beo = openmc.Material(205,name="BeO Reflector") mat_beo.add_nuclide('Be9', 4.99528E+01,'ao') mat_beo.add_nuclide('O16', 4.99042E+01+1.00046E-01,'ao') mat_beo.add_nuclide('O17', 2.00099E-02,'ao') @@ -134,7 +169,7 @@ mat_beo.add_s_alpha_beta("c_Be_in_BeO") mat_beo.add_s_alpha_beta("c_O_in_BeO") mat_list.append(mat_beo) -mat_b4c = openmc.Material(203,name="B4C Control Rods") +mat_b4c = openmc.Material(206,name="B4C Control Rods") mat_b4c.add_nuclide('B10',4,'ao') mat_b4c.add_element('C',1,'ao') mat_b4c.temperature = mat_temp_global @@ -153,8 +188,11 @@ geo_z_core_top = openmc.ZPlane(z0=25) geo_z_abs_top = openmc.ZPlane(z0=35.16,boundary_type="vacuum") geo_z_region = +geo_z_abs_bot & -geo_z_abs_top +geo_z_hp_sat = openmc.ZPlane(z0=geo_z_abs_bot.z0 + cli_args.sat_height) +geo_z_hp_pool = openmc.ZPlane(z0=geo_z_abs_bot.z0 + cli_args.pool_height) + ## Heat Pipes ## -geo_hp_cyl = [openmc.ZCylinder(r=dia/2, x0=0.0, y0=0.0) for dia in [1.092, 1.270]] +geo_hp_cyl = [openmc.ZCylinder(r=dia/2, x0=0.0, y0=0.0) for dia in [0.8920, 1.0920, 1.270]] # Pattern around the core geo_hp_steps = [math.pi*(theta+90)/180 for theta in range(0,360,45)] @@ -165,7 +203,13 @@ geo_hp_inner_regions = [openmc.ZCylinder(r=geo_hp_cyl[0].r, y0=geo_hp_pitch*math.sin(theta) ) for theta in geo_hp_steps] -geo_hp_outer_regions = [openmc.ZCylinder(r=geo_hp_cyl[1].r, +geo_hp_wick_regions = [openmc.ZCylinder(r=geo_hp_cyl[1].r, + x0=geo_hp_pitch*math.cos(theta), + y0=geo_hp_pitch*math.sin(theta) + ) for theta in geo_hp_steps] + + +geo_hp_outer_regions = [openmc.ZCylinder(r=geo_hp_cyl[2].r, x0=geo_hp_pitch*math.cos(theta), y0=geo_hp_pitch*math.sin(theta) ) for theta in geo_hp_steps] @@ -173,9 +217,16 @@ geo_hp_outer_regions = [openmc.ZCylinder(r=geo_hp_cyl[1].r, geo_hp_outer_union = openmc.Union(-reg for reg in geo_hp_outer_regions) # Fill cells for each pipe -geo_hp_inner_cells = [openmc.Cell(name=f"Heat Pipe {num} Inner Cell",region=-reg&geo_z_region,fill=mat_hp_na) for num,reg in enumerate(geo_hp_inner_regions)] +geo_hp_inner_liq_cells = [openmc.Cell(name=f"Heat Pipe {num} Liquid Inner Cell",region=-reg&+geo_z_abs_bot&-geo_z_hp_pool,fill=mat_hp_na_liq) for num,reg in enumerate(geo_hp_inner_regions)] +geo_hp_inner_vap_cells = [openmc.Cell(name=f"Heat Pipe {num} Vapour Inner Cell",region=-reg&+geo_z_hp_pool&-geo_z_abs_top,fill=mat_hp_na_vap) for num,reg in enumerate(geo_hp_inner_regions)] + +geo_hp_wick_unsat_cells = [openmc.Cell(name=f"Heat Pipe {num} Unsaturated Wick Cell",region=-reg_out&+reg_in&+geo_z_hp_sat&-geo_z_abs_top,fill=mat_hp_wick_unsat) + for num,(reg_out,reg_in) in enumerate(zip(geo_hp_wick_regions,geo_hp_inner_regions))] +geo_hp_wick_sat_cells = [openmc.Cell(name=f"Heat Pipe {num} Saturated Wick Cell",region=-reg_out&+reg_in&-geo_z_hp_sat&+geo_z_abs_bot,fill=mat_hp_wick_sat) + for num,(reg_out,reg_in) in enumerate(zip(geo_hp_wick_regions,geo_hp_inner_regions))] + geo_hp_outer_cells = [openmc.Cell(name=f"Heat Pipe {num} Outer Cell",region=-reg_out&+reg_in&geo_z_region,fill=mat_haynes230) - for num,(reg_out,reg_in) in enumerate(zip(geo_hp_outer_regions,geo_hp_inner_regions))] + for num,(reg_out,reg_in) in enumerate(zip(geo_hp_outer_regions,geo_hp_wick_regions))] ## Core Air Cylinder ## geo_air_inner_cyl = openmc.ZCylinder(r=1.9939, x0=0.0, y0=0.0) @@ -185,7 +236,6 @@ geo_air_inner_cell = openmc.Cell(name="Inner Air Cell",region=geo_air_inner_regi ## Core Fuel Block ## geo_fuel_cyl = openmc.ZCylinder(r=5.4991, x0=0.0, y0=0.0) geo_fuel_region = +geo_z_core_bot & -geo_z_core_top & -geo_fuel_cyl & +geo_air_inner_cyl & ~geo_hp_outer_union -#geo_fuel_region = +geo_z_abs_bot & -geo_z_abs_top & -geo_fuel_cyl geo_fuel_cell = openmc.Cell(name="Core Fuel Cell",region=geo_fuel_region,fill=mat_fuel) ## BeO Reflectors ## @@ -201,47 +251,31 @@ geo_refl_upper_cell = openmc.Cell(name="Upper BeO Reflector Cell",region=geo_ref ## Radial Core Surroundings ## geo_rad_cyl = [openmc.ZCylinder(r=dia/2, x0=0.0, y0=0.0) for dia in [12.70, 13.30, 14.11, 38.40]] geo_rad_cyl[-1].boundary_type="vacuum" -geo_rad_cells = [openmc.Cell(name="Intermediate Air Cell",region=+geo_fuel_cyl&-geo_rad_cyl[0]&geo_z_region,fill=mat_air), +geo_rad_cells = [openmc.Cell(name="Intermediate Air Cell",region=+geo_fuel_cyl&-geo_rad_cyl[0]&geo_z_region&~geo_hp_outer_union,fill=mat_air), openmc.Cell(name="Steel Annulus Cell",region=+geo_rad_cyl[0]&-geo_rad_cyl[1]&geo_z_region,fill=mat_ss316), openmc.Cell(name="Outer Air Cell",region=+geo_rad_cyl[1]&-geo_rad_cyl[2]&geo_z_region,fill=mat_air), openmc.Cell(name="Outer BeO Reflector Cell",region=+geo_rad_cyl[2]&-geo_rad_cyl[3]&geo_z_region,fill=mat_beo)] -geo_root_universe = openmc.Universe(name="Root Universe",cells=[geo_fuel_cell,geo_air_inner_cell,geo_refl_upper_cell,geo_refl_lower_cell]+geo_rad_cells+geo_hp_inner_cells+geo_hp_outer_cells) +geo_root_universe = openmc.Universe(name="Root Universe",cells=[geo_fuel_cell,geo_air_inner_cell,geo_refl_upper_cell,geo_refl_lower_cell]+geo_rad_cells+geo_hp_inner_liq_cells+geo_hp_inner_vap_cells+geo_hp_wick_unsat_cells+geo_hp_wick_sat_cells+geo_hp_outer_cells) geo_obj = openmc.Geometry(geo_root_universe) geo_obj.export_to_xml() ## Settings ## set_obj = openmc.Settings() -set_obj.particles = 10000 -set_obj.batches = 1500 -set_obj.inactive = 100 +set_obj.particles = 20000 +set_obj.batches = 4500 +set_obj.inactive = 500 set_obj.temperature = {"method": "interpolation"} set_src = openmc.IndependentSource() set_src.space = openmc.stats.Point((0.0,0.0,30.0)) set_src.space = openmc.stats.CylindricalIndependent( - r = openmc.stats.Uniform(geo_air_inner_cyl.r,geo_fuel_cyl.r-1.0), - #r = openmc.stats.Uniform(0,2), + r = openmc.stats.Uniform(geo_air_inner_cyl.r,geo_fuel_cyl.r), z = openmc.stats.Uniform(geo_z_core_bot.z0,geo_z_core_top.z0), phi = openmc.stats.Uniform(0.0,2*math.pi), origin = (0.0,0.0,0.0) ) -#set_src.space = openmc.stats.Box(lower_left=(-2,-2,-2),upper_right=(2,2,6)) - -## DEBUG SOURCE - -#for i in range(10000): -# sample_r = set_src.space.r.sample() -# sample_z = set_src.space.z.sample() -# sample_phi = set_src.space.phi.sample() -# #print(f"{sample_r = }, {sample_z = }, {sample_phi = }") -# x = sample_r*math.cos(sample_phi[0]) -# y = sample_r*math.sin(sample_phi[0]) -# -# cell = geo_obj.find((x,y,sample_z))[-1] -# print(f"Particle {i}: ({x}, {y}, {sample_z}) -> Cell: {cell}") -# #set_obj.verbosity=10 set_obj.export_to_xml() |