From e177b298070dad59b3687578ef313fb10e09a2ec Mon Sep 17 00:00:00 2001 From: church89 Date: Fri, 20 Jan 2023 13:11:40 +0100 Subject: [PATCH] remove file --- msre.py | 632 -------------------------------------------------------- 1 file changed, 632 deletions(-) delete mode 100644 msre.py diff --git a/msre.py b/msre.py deleted file mode 100644 index bc394ec..0000000 --- a/msre.py +++ /dev/null @@ -1,632 +0,0 @@ -import os -import openmc -import openmc.deplete -import openmc.lib -import numpy as np -from math import log10, sqrt -from collections import OrderedDict -import matplotlib.pyplot as plt -import pandas as pd -from scipy import stats -import math -from scipy.optimize import curve_fit -import numpy as np -import re - -def calc_density(temp_c): - temp_f = temp_c*9/5 + 32 #Convert from Celsius to F - temp_ref = 1181 # F or 638.3 C - rho_ref = 2.3275 # g/cm3 @ 1181 F - exp_c = -1.18e-4 # 1/F fuel salt expansion coefficient - rho_calc = rho_ref * (1 + exp_c*(temp_f-temp_ref)) #adjusted density - return rho_calc - -def define_materials(fuel_temp, u235_load=None): - # CRITICALITY BENCHMARK, definition based on U235 mass fraction of 0.01408 wt% - # fuel salt temperature of 638.3 C and salt density of 2.3275 - if u235_load is None: - rho = calc_density(fuel_temp) - tot_mass = 4340.43 * 1000 # Total fuel salt in the loop charge, fixed [g] - volume = tot_mass/rho # assuming total volume remains constant - - salt = openmc.Material(name="salt", temperature = fuel_temp + 273.15) - salt.add_nuclide('Li6',1.31480070E-05) - salt.add_nuclide('Li7', 0.262960140146177) - salt.add_nuclide('Be9',1.1863E-01) - salt.add_nuclide('Zr90',1.0543E-02) - salt.add_nuclide('Zr91',2.2991E-03) - salt.add_nuclide('Zr92',3.5142E-03) - salt.add_nuclide('Zr94',3.5613E-03) - salt.add_nuclide('Zr96',5.7375E-04) - salt.add_nuclide('Hf174',8.3786E-10) - salt.add_nuclide('Hf176',2.7545E-08) - salt.add_nuclide('Hf177',9.7401E-08) - salt.add_nuclide('Hf178',1.4285E-07) - salt.add_nuclide('Hf179',7.1323E-08) - salt.add_nuclide('Hf180',1.8370E-07) - salt.add_nuclide('U234',1.034276246E-05) - salt.add_nuclide('U235',1.009695816E-03) - salt.add_nuclide('U236',4.227809892E-06) - salt.add_nuclide('U238',2.168267822E-03) - salt.add_nuclide('Fe54',2.8551E-06) - salt.add_nuclide('Fe56',4.4818E-05) - salt.add_nuclide('Fe57',1.0350E-06) - salt.add_nuclide('Fe58',1.3775E-07) - salt.add_nuclide('Cr50',2.1224E-06) - salt.add_nuclide('Cr52',4.0928E-05) - salt.add_nuclide('Cr53',4.6409E-06) - salt.add_nuclide('Cr54',1.1552E-06) - salt.add_nuclide('Ni58',5.8597E-06) - salt.add_nuclide('Ni60',2.2571E-06) - salt.add_nuclide('Ni61',9.8117E-08) - salt.add_nuclide('Ni62',3.1284E-07) - salt.add_nuclide('Ni64',7.9671E-08) - salt.add_nuclide('O16',5.1437E-04) - salt.add_nuclide('O17',1.8927E-07) - salt.add_nuclide('O18',9.6440E-07) - salt.add_nuclide('F19',5.9409E-01) - salt.set_density('g/cm3', rho_calc) - - else: - # REACTIVITY BENCHMARK - rho = calc_density(fuel_temp) - # Salt composition in kg with 65.25 kg U235 at 649 C - salt_comp = {'Li6':507.27*0.00005,'Li7':507.27*0.99995,'Be':293.96, - 'Zr':513.97,'Hf':0.0029,'U234':0.67,'U235':65.25,'U236':0.27, - 'U238':141.91,'Fe':0.75,'Cr':0.13,'Ni':0.14,'O':2.27,'F':3103.22} - tot_mass = sum(salt_comp.values())*1000 - # Calculate total volume assuming fuel salt composition at first load reamins fixed - volume = tot_mass/rho - - # HEU additions mass fraction - heu = {'Li6':4.93*0.00005,'Li7':4.93*0.99995,'U234':0.59, - 'U235':57.7,'U236':0.24,'U238':3.54,'F':33} - - # Add u235 - mass_add = u235_load - salt_comp['U235'] - for comp,mass in salt_comp.items(): - if comp in heu: - salt_comp[comp] += mass_add/heu['U235']*heu[comp] - # Recalculate total mass - tot_mass = sum(salt_comp.values())*1000 - - salt = openmc.Material(name="salt", temperature = fuel_temp + 273.15) - for comp, mass in salt_comp.items(): - splt = re.split(r'\d+', comp) - if len(splt) == 2: - salt.add_nuclide(comp, mass/tot_mass, 'wo') - else: - for nuc,frac in openmc.data.isotopes(comp): - salt.add_nuclide(nuc, mass*frac/tot_mass, 'wo') - # Recalculate density as new total mass and fixed volume - salt.set_density('g/cm3',tot_mass/volume) - #Get U235 mass - salt.volume = volume #Assign volume as Loop total volume in cm3 - - #moderator blocks170 - graphite = openmc.Material(name='graphite',temperature=fuel_temp + 273.15) - graphite.set_density('g/cm3',1.86) - graphite.add_nuclide('C12',1) - graphite.add_s_alpha_beta('c_Graphite') - - #inor-8 - inor = openmc.Material(name='inor-8',temperature=fuel_temp + 273.15) - inor.set_density('g/cm3',8.7745) - inor.add_element('Ni',(66+71)/2,'wo') - inor.add_element('Mo',(15+18)/2,'wo') - inor.add_element('Cr',(6+8)/2,'wo') - inor.add_element('Fe',5,'wo') - inor.add_element('C',(0.04+0.08)/2,'wo') - inor.add_element('Al',0.25,'wo') - inor.add_element('Ti',0.25,'wo') - inor.add_element('S',0.02,'wo') - inor.add_element('Mn',1.0,'wo') - inor.add_element('Si',1.0,'wo') - inor.add_element('Cu',0.35,'wo') - inor.add_element('B',0.010,'wo') - inor.add_element('W',0.5,'wo') - inor.add_element('P',0.015,'wo') - inor.add_element('Co',0.2,'wo') - - #helium - helium = openmc.Material(name='helium') - helium.add_element('He',1.0) - helium.set_density('g/cm3',1.03*(10**-4)) - - #Control rods inconel clad - trace = 0.01 - inconel = openmc.Material(name='inconel', temperature = 65.6 + 273.15) - inconel.add_element('Ni',78.5,percent_type='wo') - inconel.add_element('Cr',14.0,percent_type='wo') - inconel.add_element('Fe',6.5,percent_type='wo') - inconel.add_element('Mn',0.25,percent_type='wo') - inconel.add_element('Si',0.25,percent_type='wo') - inconel.add_element('Cu',0.2,percent_type='wo') - inconel.add_element('Co',0.2,percent_type='wo') - inconel.add_element('Al',0.2,percent_type='wo') - inconel.add_element('Ti',0.2,percent_type='wo') - inconel.add_element('Ta',0.5,percent_type='wo') - inconel.add_element('W',0.5,percent_type='wo') - inconel.add_element('Zn',0.2,percent_type='wo') - inconel.add_element('Zr',0.1,percent_type='wo') - inconel.add_element('C',trace,percent_type='wo') - inconel.add_element('Mo',trace,percent_type='wo') - inconel.add_element('Ag',trace,percent_type='wo') - inconel.add_element('B',trace,percent_type='wo') - inconel.add_element('Ba',trace,percent_type='wo') - inconel.add_element('Be',trace,percent_type='wo') - inconel.add_element('Ca',trace,percent_type='wo') - inconel.add_element('Cd',trace,percent_type='wo') - inconel.add_element('V',trace,percent_type='wo') - inconel.add_element('Sn',trace,percent_type='wo') - inconel.add_element('Mg',trace,percent_type='wo') - inconel.set_density('g/cm3',8.5) - - # SS 316 control rod flexible hose - ss316 = openmc.Material(name='ss316', temperature = 65.6 + 273.15) - ss316.add_element('C',0.026,'wo') - ss316.add_element('Si',0.37,'wo') - ss316.add_element('Mn',0.16,'wo') - ss316.add_element('Cr',16.55,'wo') - ss316.add_element('Cu',0.16,'wo') - ss316.add_element('Ni',10,'wo') - ss316.add_element('P',0.029,'wo') - ss316.add_element('S',0.027,'wo') - ss316.add_element('Mo',2.02,'wo') - ss316.add_element('N',0.036,'wo') - ss316.add_element('Fe',70.622,'wo') - ss316.set_density('g/cm3',7.99) - - #Control rods bushing posion material - Gd2O3 = openmc.Material() - Gd2O3.add_element('Gd',2) - Gd2O3.add_element('O',3) - Gd2O3.set_density('g/cm3',7.41) - Al2O3 = openmc.Material() - Al2O3.add_element('Al',2) - Al2O3.add_element('O',3) - Al2O3.set_density('g/cm3',3.95) - bush = openmc.Material.mix_materials([Gd2O3,Al2O3],[0.7,0.3],'wo') - bush.name='bush' - bush.temperature = 65.6 +273.15 - - #Concrete block - concrete = openmc.Material(name='concrete') - concrete.add_element('H',0.005,'wo') - concrete.add_element('O',0.496,'wo') - concrete.add_element('Si',0.314,'wo') - concrete.add_element('Ca',0.083,'wo') - concrete.add_element('Na',0.017,'wo') - concrete.add_element('Mn',0.002,'wo') - concrete.add_element('Al',0.046,'wo') - concrete.add_element('S',0.001,'wo') - concrete.add_element('K',0.019,'wo') - concrete.add_element('Fe',0.012,'wo') - concrete.set_density('g/cm3',2.35) - - #Thermal shielding as water and SS305 (50-50) - water = openmc.Material() - water.add_element('H',2) - water.add_element('O',1) - water.set_density('g/cm3',0.997) - - #stainless steel 304 - ss304 = openmc.Material() - ss304.add_element('C',0.08,'wo') - ss304.add_element('Mn',2,'wo') - ss304.add_element('P',0.045,'wo') - ss304.add_element('S',0.03,'wo') - ss304.add_element('Si',0.75,'wo') - ss304.add_element('Cr',19,'wo') - ss304.add_element('Ni',10,'wo') - ss304.add_element('N',0.1,'wo') - ss304.add_element('Fe',67.995, 'wo') - ss304.set_density('g/cm3',7.93) - shield = openmc.Material.mix_materials([water,ss304],[0.5,0.5],'vo') - shield.temperature = 32.2 + 273.15 - shield.name='steelwater' - - # "Careytemp 1600" by Philip Carey Manufacturing Compamy (Cincinnati) from http://moltensalt.org/references/static/downloads/pdf/ORNL-TM-0728.pdf - insulation=openmc.Material(name='insulation') - insulation.add_element('Si',1) - insulation.add_element('O',2) - insulation.set_density('g/cm3',0.16) #https://www.osti.gov/servlets/purl/1411211 - - # sand water, not sure about this material - sandwater=openmc.Material(name='sandwater') - sandwater.add_element('Fe',3) - sandwater.add_element('O',4) - sandwater.set_density('g/cm3',6) - - #Vessel anular steel - steel = openmc.Material(name='steel') - steel.add_element('Fe',1) - steel.set_density('g/cm3',7.85) - - mats = openmc.Materials([salt,graphite,inor,helium,inconel,shield,concrete,steel,ss316,sandwater,insulation,bush]) - return mats - -def build(make_tally=True, plot_geom=True, u235_load=None, fuel_temp=638.3, cr1_pos=51, cr2_pos=51, cr3_pos=51): - #Clean-up - os.system("rm *.xml *.h5 *.out") - - # CAD h5m files - core_h5m = 'h5m/msre_reactor_1e-2.h5m' - control_rod1_h5m = 'h5m/msre_control_rod_1e-2.h5m' - - #Materials - mats=define_materials(fuel_temp,u235_load) - - #Geometry - core = openmc.DAGMCUniverse(filename=core_h5m, auto_geom_ids=True, universe_id=1) - control_rod1 = openmc.DAGMCUniverse(filename=control_rod1_h5m, auto_geom_ids=True, universe_id=2) - core_region = core.bounding_region() - cr1_region = control_rod1.bounding_region(boundary_type='transmission', starting_id=20000) - # Extend control rod region to include upwards translations - cr1_region = cr1_region | cr1_region.translate([0,0,150]) - cr2_region = cr1_region.translate([-10.163255,0,0]) - cr3_region = cr1_region.translate([-10.163255,10.163255,0]) - core_cell = openmc.Cell(region=~(cr1_region | cr2_region | cr3_region) & core_region , fill=core) - cr1_cell = openmc.Cell(name='CR1', region=cr1_region, fill=control_rod1) - cr2_cell = openmc.Cell(name='CR2', region=cr2_region, fill=control_rod1) - cr3_cell = openmc.Cell(name='CR3', region=cr3_region, fill=control_rod1) - #translate control rods - start_pos = 19.2 - setattr(cr1_cell, 'translation', [0, 0, start_pos + cr1_pos*2.54]) - setattr(cr2_cell, 'translation', [-10.163255, 0, start_pos + cr2_pos*2.54]) - setattr(cr3_cell, 'translation', [-10.163255, 10.163255, start_pos + cr3_pos*2.54]) - geometry = openmc.Geometry([core_cell,cr1_cell,cr2_cell,cr3_cell]) - - # Settings - settings = openmc.Settings() - settings.temperature = {'method':'interpolation','range':(293.15,923.15)} - settings.batches = 50 - settings.inactive = 20 - settings.particles = 30000 - settings.photon_transport = False - source_area = openmc.stats.Box([-100., -100., 0.],[ 100., 100., 200.],only_fissionable = True) - settings.source = openmc.Source(space=source_area) - - if make_tally: - tally = OrderedDict() - tally['general'] = openmc.Tally(name="General") - if settings.photon_transport: - heating_score = 'heating' - else: - heating_score = 'heating-local' - tally['general'].scores.append(heating_score) - - e_min, e_max = 1e-5, 20e6 - groups = 500 - energies = np.logspace(log10(e_min), log10(e_max), groups + 1) - energy_filter = openmc.EnergyFilter(energies) - particle_filter = openmc.ParticleFilter(['neutron']) - cell_filter = openmc.MaterialFilter([mats[0]]) - # mesh = openmc.RegularMesh() - # mesh.dimension = [1, 1, 1] - # mesh.lower_left = [39.7, -3.1, 39.2] - # mesh.upper_right = [41.7, -2, 195] - # mesh_filter = openmc.MeshFilter(mesh) - tally['flux']= openmc.Tally(name="flux") - tally['flux'].filters = [energy_filter, particle_filter] #, cell_filter]#, mesh_filter] - tally['flux'].scores = ['flux'] - - mesh = openmc.RegularMesh() - mesh.dimension = [500, 500, 1] - mesh.lower_left = [-100, -100, 50] - mesh.upper_right = [100, 100, 200] - mesh_filter = openmc.MeshFilter(mesh) - tally['mesh'] = openmc.Tally(name="Mesh") - tally['mesh'].scores = ['flux','absorption','fission','scatter'] - tally['mesh'].filters = [mesh_filter] - tally['mesh'].filters.append(particle_filter) - - tally['leak'] = openmc.Tally(name='leakage') - mesh = openmc.RegularMesh() - mesh.dimension = [1, 1, 1] - mesh.lower_left = core.bounding_box[0]-10 - mesh.width = (core.bounding_box[1]+10)*2 - meshsurface_filter = openmc.MeshSurfaceFilter(mesh) - tally['leak'].filters = [meshsurface_filter] - tally['leak'].scores = ['current'] - tallies = openmc.Tallies(tally.values()) - model = openmc.model.Model(geometry,mats,settings,tallies) - else: - model = openmc.model.Model(geometry,mats,settings) - - if plot_geom: - colors = {[i for i in mats if i.name=="salt"][0]: 'yellow', - [i for i in mats if i.name=="graphite"][0]: 'black', - [i for i in mats if i.name=="inor-8"][0]: 'grey', - [i for i in mats if i.name=="helium"][0]: 'cyan', - [i for i in mats if i.name=="inconel"][0]: 'grey', - [i for i in mats if i.name=="bush"][0]: 'blue', - [i for i in mats if i.name=="ss316"][0]: 'grey', - [i for i in mats if i.name=="concrete"][0]: 'brown', - [i for i in mats if i.name=="steelwater"][0]: 'red', - [i for i in mats if i.name=="insulation"][0]: 'green', - [i for i in mats if i.name=="sandwater"][0]: 'lightgreen', - [i for i in mats if i.name=="steel"][0]: 'grey'} - - plot_file = openmc.Plots() - plot1 = openmc.Plot() - plot1.width = [150, 150] - plot1.pixels = [2000, 2000] - plot1.origin = [0,0,150] - plot1.basis = 'xy' - plot1.color_by = "material" - plot1.colors = colors - - plot2 = openmc.Plot() - plot2.width = [1000, 2000] - plot2.pixels = [4000,8000] - plot2.origin = [-5,0,150] - plot2.basis = 'yz' - plot2.color_by = "material" - plot2.colors = colors - - plot3 = openmc.Plot() - plot3.width = [1000, 2000] - plot3.pixels = [4000,8000] - plot3.origin = [0,-5,150] - plot3.basis = 'xz' - plot3.color_by = "material" - plot3.colors = colors - - model.plots.append(plot1) - model.plots.append(plot2) - model.plots.append(plot3) - model.plot_geometry() - - return model - -def run(model, mass, power): - results=model.run() - salt_vol = 1.65058e6 #fuel salt core volume [cm3] from OnShape model - #core_vol = 342.76 # core vol fuel cell [cc] - sp = openmc.StatePoint(results) - heating = sp.get_tally(name="General").get_pandas_dataframe()["mean"].sum()*openmc.data.JOULE_PER_EV - fac = power/heating - t = sp.get_tally(name="flux") - energy_filter = t.filters[0] - energies = energy_filter.bins[:, 0] - mean = t.mean.ravel() - uncertainty = t.get_values(value='std_dev').ravel() - - fig, ax = plt.subplots() - ax.plot(energies, mean*fac/salt_vol, drawstyle='steps-post') - ax.set_xlabel('Energy [eV]') - ax.set_ylabel(r'Flux [neutrons/cm$^2$-s]') - ax.set_xscale('log') - #ax.set_yscale('log') - ax.grid(True, which='both') - #plt.title('Neutrons spectrum') - plt.savefig("norm_spectrum",dpi=600) - # - values = sp.get_tally(name="Mesh").get_slice(scores=['flux']).get_pandas_dataframe()["mean"] - values = values.values.reshape(500,500) - fig, ax = plt.subplots() - pos = ax.imshow(values*fac/salt_vol, - aspect='auto', - origin='lower') - cbar = plt.colorbar(pos,ax=ax,label=r'Flux [neutrons/cm$^2$-s]') - plt.savefig("flux",dpi=600) - - values = sp.get_tally(name="Mesh").get_slice(scores=['fission']).get_pandas_dataframe()["mean"] - values = values.values.reshape(500,500) - fig, ax = plt.subplots() - pos = ax.imshow(values*fac/salt_vol, - aspect='auto', - origin='lower') - cbar = plt.colorbar(pos,ax=ax,label=r'Fission [-/cm$^2$-s]') - plt.savefig("fission",dpi=600) - - values = sp.get_tally(name="Mesh").get_slice(scores=['absorption']).get_pandas_dataframe()["mean"] - values = values.values.reshape(500,500) - fig, ax = plt.subplots() - pos = ax.imshow(values*fac/salt_vol, - aspect='auto', - origin='lower') - cbar = plt.colorbar(pos,ax=ax,label=r'Absorption [-/cm$^2$-s]') - plt.savefig("abs",dpi=600) - - values = sp.get_tally(name="Mesh").get_slice(scores=['scatter']).get_pandas_dataframe()["mean"] - values = values.values.reshape(500,500) - fig, ax = plt.subplots() - pos = ax.imshow(values*fac/salt_vol, - aspect='auto', - origin='lower') - cbar = plt.colorbar(pos,ax=ax,label=r'Scatter [-/cm$^2$-s]') - plt.savefig("scatter",dpi=600) - -def depletion(model, mass, power): - vol=mass/2.3275*1000 # total volume of fuel salt [cm3] - model.materials[0].volume=vol - op = openmc.deplete.CoupledOperator(model, normalization_mode = "energy-deposition", chain_file='/home/lorenzo/Documents/ca_depletion_chains/ENDF-B-VIII.0_chain_msr.xml') - msr = openmc.deplete.msr.MsrContinuous(op,model) - msr.set_removal_rate('salt', ['Xe','Kr'], 4.067e-5) - msr.set_removal_rate('salt', ['Se','Nb','Mo','Tc','Ru','Rh','Pd','Ag','Sb','Te'], 8.777e-3) - - # Sets MSR batchwise features - msr_bw_geom = openmc.deplete.msr.MsrBatchwiseGeom(op, model, axis = 2, - cell_id_or_name = 'CR1', - bracket = [-2,5], #cm - bracket_limit = [-19.2, 51*2.54], #cm - nuc_density_limit=1e20, #cm - tol=0.1) #cm - integrator = openmc.deplete.CECMIntegrator(op, [5,5,30,30,30,180,95], msr_continuous=msr, - msr_batchwise=msr_bw_geom, timestep_units='d', power=power) - - integrator.integrate(final_step = False) - -def control_rod_worth(fuel_temp, u235_load): - drho = [] - model = build(make_tally=False, plot_geom=False, fuel_temp = fuel_temp, u235_load=u235_load) - cell = model.geometry.get_cells_by_name('CR1')[0] - rod_x = [49, 45, 41, 37, 33, 29, 25, 21, 17, 13, 9, 5, 2, 0] - for pos in rod_x: - rho = [] - pos *= 2.54 - dx = 4 * 2.54 - for shift in [-dx,dx]: - setattr(cell, 'translation', [0, 0, 19.2+pos+shift]) - res=model.run() - with openmc.StatePoint(res) as sp: - keff=sp.keff.n - pcm = (keff-1)/keff*1e5 - rho.append(pcm) - drho.append(abs(np.diff(rho))/(2*dx)) - print(drho) - plt.figure() - plt.plot(rod_x, drho, marker='x', label='Sim') - drho_exp = [6.7, 10, 13.9, 17.7, 20.8, 23, 24.2, 24.4, 23.3, 20.9, 17.5, 13.2, 9.6, 7.1] - err_exp = [0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.7, 0.7, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2] - plt.errorbar(rod_x, drho_exp, yerr=err_exp, label='Exp') - plt.legend() - plt.xlabel('Withdrawn of control rod n. 1 [inch]') - plt.ylabel('Reactivity worth') - plt.savefig('reac_rod_worth_4inches.png') - -def func (x,a,b): - return a*x +b - -def feedback_isothermal(temperatures,u235_load): - # temperature arguments in Celsius - y_data = [] - y_err = [] - for temp in temperatures: - model = build(make_tally=False, plot_geom=False, fuel_temp = temp, u235_load=u235_load) - res=model.run() - with openmc.StatePoint(res) as sp: - k=sp.keff.n - y_data.append((k-1)/k) - y_err.append(sp.keff.s) - - x_data = np.array(temperatures) - y_data = np.array(y_data) - y_err = np.array(y_err) - popt, pcov = curve_fit(func, x_data, y_data) - - plt.figure() - plt.errorbar(x_data, y_data, yerr=y_err, label='Sim data') - symb = r'$^\circ$C' - #propagated error - err = r'$\pm${:.4f}'.format(sqrt(sum(y_err**2))*1e5) - plt.plot(x_data, func(x_data, *popt), '--', label='{:.3f} [pcm/{}]'.format(popt[0]*1e5,symb)) - plt.legend() - plt.xlabel('Temperature '+symb) - plt.ylabel(r'$\rho\,\pm\sigma$') - plt.title(r'Isothermal temperature coefficient, {u235_load} kg $^{235}U$ in loop') - plt.savefig(f'iso_temperature_feedback_{str(u235_load)}kg.png') - -def rod_bank(temp, positions, loads): - load_results = {} - for load in loads: - y_data = [] - y_err = [] - # Rods position in inches where 0 is fully inserted and 51 fully withdrawn - for pos in positions: - model = build(make_tally=False, plot_geom=False, fuel_temp=temp, u235_load=load, cr1_pos=pos, cr2_pos=pos, cr3_pos=pos) - res=model.run() - with openmc.StatePoint(res) as sp: - y_data.append(sp.keff.n) - y_err.append(sp.keff.s) - load_results[load] = y_data,y_err - df=pd.DataFrame(load_results) - df.to_csv('rod_bank') - -def mass_reactivity(loads): - y_data = [] - for load in loads: - model = build(make_tally=False, plot_geom=False, fuel_temp=648.9, u235_load=load, cr1_pos=51, cr2_pos=51, cr3_pos=51) - res=model.run() - with openmc.StatePoint(res) as sp: - y_data.append(sp.keff.n) - loads = np.array(loads) - y_data = np.array(y_data) - print(y_data) - x_data=(loads-loads[0])/loads[0] - y_data=(y_data-y_data[0])/(y_data*y_data[0]) - popt, pcov = curve_fit(func, x_data, y_data) - plt.plot(x_data,y_data,'o',label='data') - plt.plot(x_data, func(x_data, *popt), '--', label='fit: f(x)={:.3f}x + {:.3f}'.format(popt[0],popt[1])) - plt.legend() - plt.ylabel(r'Change of reactivity, $\Delta k /(k_1k_2)$') - plt.xlabel(r'Fractional change of $^{235}U$ loading in loop, $\Delta m /m$') - plt.savefig('mass_reactivity') - -def triton_adder(mass): - - df = pd.read_csv('ADDER/msre_simple_fuel_mod.csv') - res = openmc.deplete.Results(f'depletion_results.h5') - vol=mass/2.3275*1000 - mats = dict() - mats['fiss'] = ['U235','Pu239','U239','Pu238','H3','Cs137'] - mats['fp'] = ['Ce144','Ce141','I131','Mo99','Xe135','Kr91','Pr145','Sm149'] - - fuel_id = [str(mat.id) for mat in res.export_to_materials(0) if mat.name == 'salt'][0] - t_omc = res.get_atoms(fuel_id,'U235')[0]/3600/24 - t_add = df.loc[np.where(df[df.columns[0]] == 'times')[0][0]][1:].values.astype(float) - - for cat, nucs in mats.items(): - div = round(len(nucs)/2) - fig,ax = plt.subplots(div,2,figsize=(12,12)) - nuc_index=0 - for i in range(div): - for j in range(2): - ax[i,j].plot(t_omc, res.get_atoms(fuel_id,nucs[nuc_index])[1]/vol, marker='x', markersize=9, color='r',label='openmc-msr') - ax[i,j].plot(t_add, df.loc[np.where(df[df.columns[0]] == nucs[nuc_index])[0][0]][1:].values.astype(float), marker='o', markersize=9, color='b',label='mcnp-adder') - ax[i,j].legend() - nuc_index +=1 - - for a,n in zip(ax.flat, nucs): - a.set(xlabel='EFPD [d]', ylabel=f'{n} [atoms/cc]') - plt.tight_layout() - plt.savefig(f'msre_openmc-vs-adder_{cat}', dpi=600) - - tt_add=np.unique(t_add) - args_add = [] - for _t in tt_add: - args_add.append(np.argwhere(t_add == _t)[0][0]) - - means = [] - lim=30 - for nuc in df[df.columns[0]][9:].values: - try: - if res.get_atoms('1',nuc)[1].mean()>1e15: - omc = res.get_atoms('1',nuc)[1][1:]/vol - add = df.loc[np.where(df[df.columns[0]] == nuc)[0][0]][1:].values.astype(float).take(args_add)[1:] - diff = (add-omc)/add *100 - diff = diff.mean() - if not math.isnan(diff): - if abs(diff) < lim: - means.append(abs(diff)) - except: - continue - pd.DataFrame(means).describe() - params = stats.gamma.fit(means) - x = np.linspace(0, round(max(means)), 1000) - pdf = stats.gamma.pdf(x, *params) - plt.figure() - plt.plot(x, pdf, label='Gamma func. data fit') - plt.hist(means,round(max(means))*2,density=True,label='Data') - plt.legend() - plt.xlabel('Relative error [%]', weight='bold') - plt.ylabel('Probability', weight='bold') - plt.title(f'Relative error distribution for {len(means)} most abundant nuclides in fuel salt below {lim}%',fontsize=9) - plt.savefig(f'rel_error_below{lim}%', dpi=600) - -if __name__ == '__main__': - mass = 4590 #tot fuel salt mass [kg] - power = 8e6 #total thermal power [W] - fuel_temp = 638.3#fuel temperature at the initial criticality point reported [C] - #run(build(make_tally=True, plot_geom=True, u235_load=65.25, cr1_pos=51,fuel_temp=648.9), mass, power) - #control_rod_worth(648.9, 65.25) - #rod_bank(648.9, [0,51],[67.94,69.94,71.71]) - #feedback_isothermal([598.9,648.9,698.9],67.86) - #feedback_isothermal([598.9,648.9,698.9],71.71) - #feedback_isothermal([598.9,648.9,698.9],69.85) - #feedback_fuel([598.9])#,648.9,698.9]) - #mass_reactivity([65.25,66,67,68,69,70,71,72]) - depletion(build(make_tally=False, plot_geom=False, u235_load=71.71, cr1_pos=0, fuel_temp=648.9), mass, power) - #triton_adder(mass)