From 29c8a0faf07d6a020f091d7d577e770bf28200b0 Mon Sep 17 00:00:00 2001 From: church89 Date: Thu, 30 Mar 2023 15:47:59 +0200 Subject: [PATCH] general improvements --- msre_depletion_post_processing.py | 70 +++++++++++++++++++------------ msre_power_history.py | 42 ++++++++++++------- 2 files changed, 70 insertions(+), 42 deletions(-) diff --git a/msre_depletion_post_processing.py b/msre_depletion_post_processing.py index 70e7faf..c71f09d 100644 --- a/msre_depletion_post_processing.py +++ b/msre_depletion_post_processing.py @@ -1,17 +1,22 @@ import openmc +import openmc.deplete import os import re import matplotlib.pyplot as plt import seaborn as sns +import datetime +import pandas as pd +import numpy as np +from pathlib import Path regex = re.compile(r'(\d+|\s+)') -path_to_results = 'depletion_results.h5' +path_to_results = '/home/lorenzo/mnt/uranium/msre/depletion_results.h5' save_dir = Path(os.path.realpath(path_to_results)).parent mats = {mat.name: mat.id for mat in openmc.material.Materials.from_xml(save_dir / 'materials.xml') if mat.depletable} results = openmc.deplete.Results(path_to_results) -t, keff = results.get_keff() +_, keff = results.get_keff() n_xe = 0 n_kr = 0 for nuc,_ in openmc.data.isotopes('Xe'): @@ -19,21 +24,29 @@ for nuc,_ in openmc.data.isotopes('Xe'): for nuc,_ in openmc.data.isotopes('Kr'): n_kr += results.get_atoms(str(mats['salt']), nuc)[1] -# Let's convert time from sec to days -t /= (3600 * 24) +df=pd.read_excel('/home/lorenzo/mnt/uranium/msre/MSRE_235_233_Power_History_R24E.xlsx') +# U235 operation run +dt = df['Duration (h)'][:94] +# First operation in MW range(ORNL-4674) +t0=datetime.datetime.strptime('24/01/1966', '%d/%m/%Y') +date_times = [] +for t in dt: + date_times.append(t0) + t0 += datetime.timedelta(hours=t) +date_times.append(t0) -plt.figure() +power=df['Power (MWth)'][:95] +plt.figure(figsize=(18,10)) ax = plt.subplot() -k1, = ax.plot(t, [k[0] for k in keff], '--', c='red', label='keff') ax1 = ax.twinx() -n2, = ax1.plot(t, n_xe, c='green', label='Xe') -n4, = ax1.plot(t, n_kr, c='blue', label='Kr') -ax.set_xlabel('Time[d]') -ax.set_ylabel(r'$k_{eff}$', color='r') -ax.tick_params(axis='y', colors='red') -ax1.set_yscale('log') -ax1.set_ylabel('Nuclides [atoms]') -ax1.legend(handles=[k1, n2, n4]) +ax.errorbar(date_times, [k[0] for k in keff], [k[1] for k in keff], marker='o', color='black', markersize=4, fmt=' ') + +ax1.step(date_times, power, where='post', color='red') +ax.set_ylabel(r'$k_{eff}\pm\sigma$') +ax.set_ylim(0.98,1.01) +ax.tick_params(axis='y') +ax1.set_ylabel('Power [MWth]',color='red') +ax1.tick_params(axis='y', colors='red') plt.savefig(f'{save_dir}/keff', dpi=600) # Microscopic absorption cross section at 0.0253 eV @@ -43,9 +56,8 @@ _, n_xe135 = results.get_atoms(str(mats['salt']), 'Xe135') _, n_u235 = results.get_atoms(str(mats['salt']), 'U235') # Poison fraction pf = (xs_xe135*n_xe135)/(xs_u235*n_u235)*100 -plt.figure() -plt.plot(t, pf) -plt.xlabel('Time [d]') +plt.figure(figsize=(18,10)) +plt.plot(date_times, pf) plt.ylabel('Xe posion fraction [%]') plt.savefig(f'{save_dir}/fission_products', dpi=600) @@ -56,18 +68,23 @@ for nuc,_ in openmc.data.isotopes('U'): for nuc in ['Pu238','Pu239','Pu240','Pu241','Pu242']: inventory[nuc] = results.get_atoms(str(mats['salt']), nuc)[1] / openmc.data.AVOGADRO * openmc.data.atomic_mass(nuc) / 1000 -plt.figure() +plt.figure(figsize=(18,10)) for nuc, mass in inventory.items(): - plt.plot(t, mass, label=nuc) -plt.xlabel('Time [y]') + plt.plot(date_times, mass, label=nuc) plt.ylabel('Mass [g]') plt.yscale('log') plt.ylim(1e-5) plt.legend() plt.savefig(f'{save_dir}/inventory', dpi=600) -# All nuclides present in the fuel at last time-step -nucs = results.export_to_materials(-1)[0].get_nuclides() +lib = openmc.data.DataLibrary.from_xml(os.environ.get("OPENMC_CROSS_SECTION")) +nuclides = set() +for library in lib.libraries: + if library['type'] != 'neutron': + continue + for name in library['materials']: + if name not in nuclides: + nuclides.add(name) # Let's begin by making some useful groupings gaseos = ['H', 'He', 'Ne', 'Ar', 'Kr', 'Xe', 'Rn'] #gaseous fission products @@ -82,7 +99,7 @@ m_a = ['Ac','Th','Pa','Np','Am','Cm','Bk','Cf','Es','Fm','Md','No','Lr'] # Get fissile nuclides in the fuel, based on Ronen's rule for determining fissile isotopes fissile = [] -for nuc in nucs: +for nuc in nuclides: elm = regex.split(nuc)[0] a = round(openmc.data.atomic_mass(nuc)) z = openmc.data.ATOMIC_NUMBER[elm] @@ -99,7 +116,7 @@ for nuc in fissile: nuclides_stack = dict() groups_stack = {'Gaseos':0, 'Noble metals':0, 'Metals':0, 'Halogens':0 , 'Alkali metals':0, 'Alkali earths':0, 'Lanthanides':0, 'MA':0, 'Others':0} -for nuc in nucs: +for nuc in nuclides: if regex.split(nuc)[0] in ['U','Pu']: nuclides_stack[nuc] = results.get_reaction_rate(str(mats['salt']), nuc, 'fission')[1] @@ -161,7 +178,7 @@ groups_stack.update(u_series) groups_stack.update(pu_series) plt.figure(figsize=(15,10)) -plt.stackplot(t, groups_stack.values(), labels=groups_stack.keys(), +plt.stackplot(date_times, groups_stack.values(), labels=groups_stack.keys(), edgecolor="black", linewidth=0.5,colors=colors, alpha=0.8) handles, labels = plt.gca().get_legend_handles_labels() legend = plt.legend([handles[idx] for idx in list(reversed(np.arange(0,len(handles),1)))], @@ -169,8 +186,7 @@ legend = plt.legend([handles[idx] for idx in list(reversed(np.arange(0,len(handl bbox_to_anchor=(1.05,1), loc='upper left', borderaxespad=0, ncol=2, fontsize=15) -plt.xlabel('Time [d]',weight='bold',fontsize=17) -plt.title('Neutrons absorption distribution per neutron absorbed in fissile isotopes', +plt.title('Fuel salt neutrons absorption distribution per neutron absorbed in fissile isotopes', weight='bold', fontsize=17) plt.xticks(fontsize=13) plt.yticks(fontsize=13) diff --git a/msre_power_history.py b/msre_power_history.py index 6202062..e8f9a7b 100644 --- a/msre_power_history.py +++ b/msre_power_history.py @@ -50,7 +50,7 @@ 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', 2.3275) +salt.set_density('g/cm3', 2.32151) salt.volume = 4560/2.32151*1000 #moderator blocks170 @@ -228,12 +228,10 @@ cr3_cell = openmc.Cell(name='CR3', region=cr3_region, fill=control_rod1) #Fix control rods initial positions start_pos = 19.2 #cm, geometrical distance between lower bottom and starting point -cr1_pos = 150 * 2.54 # cm, initial position of control rod1 with respect to start post -cr2_pos = 150 * 2.54 -cr3_pos = 150 * 2.54 -setattr(cr1_cell, 'translation', [0, 0, start_pos + cr1_pos]) -setattr(cr2_cell, 'translation', [-offset, 0, start_pos + cr2_pos]) -setattr(cr3_cell, 'translation', [-offset, offset, start_pos + cr3_pos]) +top_pos = 51 * 2.54 # cm, initial position of control rod1 with respect to start post +setattr(cr1_cell, 'translation', [0, 0, start_pos + top_pos]) +setattr(cr2_cell, 'translation', [-offset, 0, start_pos + top_pos]) +setattr(cr3_cell, 'translation', [-offset, offset, start_pos + top_pos]) geometry = openmc.Geometry([core_cell,cr1_cell,cr2_cell,cr3_cell]) # Settings @@ -268,16 +266,16 @@ plots = openmc.Plots() plot = openmc.Plot() plot.basis = 'xy' plot.width = (150,150) -plot.pixels = (500,500) -plot.origin = (0,0,start_pos + cr1_pos) +plot.pixels = (1500,1500) +plot.origin = (0,0,160) plot.color_by = 'material' plot.colors = colors model.plots.append(plot) plot = openmc.Plot() plot.basis = 'xz' plot.width = (150,300) -plot.pixels = (500,1000) -plot.origin = (0,-5,start_pos + cr1_pos) +plot.pixels = (750,1500) +plot.origin = (0,-5,150) plot.color_by = 'material' plot.colors = colors model.plots.append(plot) @@ -305,12 +303,26 @@ msr.set_removal_rate('salt', df=pd.read_excel('MSRE_235_233_Power_History_R24E.xlsx') #U235 -timesteps = np.cumsum(df['Duration (d)'][:94]).values -power = df['Power (MWth)'][:94].values +timesteps = df['Duration (d)'][:94].values +power = df['Power (MWth)'][:94].values*1000000 +#Add one further timestep at the end of every power run +timesteps_ext = [] +power_ext = [] +for i in range(len(timesteps)): + power_ext.append(power[i]) + if power[i] != 0.0: + # duplicate power value + power_ext.append(power[i]) + # add timestep - 1 sec + timesteps_ext.append(timesteps[i] - 1/3600/24) + # add 1 sec + timesteps_ext.append(1/3600/24) + else: + timesteps_ext.append(timesteps[i]) -integrator = openmc.deplete.CECMIntegrator(op, timesteps=timesteps, +integrator = openmc.deplete.CECMIntegrator(op, timesteps=timesteps_ext, msr_continuous=msr, #msr_batchwise=msr_bw_geom, - timestep_units='d', power=power) + timestep_units='d', power=power_ext) integrator.integrate()