Add Fuel cycle benchmarking

doc/content/bib/tmap8.bib

@@ -287,3 +287,14 @@ @article{li2010analytical
   year={2010},
   publisher={American Society of Civil Engineers}
 }
+
+@article{meschini2023modeling,
+  title={Modeling and analysis of the tritium fuel cycle for ARC-and STEP-class DT fusion power plants},
+  author={Meschini, Samuele and Ferry, Sara E and Delaporte-Mathurin, R{\'e}mi and Whyte, Dennis G},
+  journal={Nuclear Fusion},
+  volume={63},
+  number={12},
+  pages={126005},
+  year={2023},
+  publisher={IOP Publishing}
+}

doc/content/verification_and_validation/figures/comparison_fuel_cycle_benchmark.py

@@ -0,0 +1 @@
+../../../../test/tests/fuel_cycle_benchmark/comparison_fuel_cycle_benchmark.py

doc/content/verification_and_validation/fuel_cycle_benchmarking.md

@@ -0,0 +1,132 @@
+# Fuel Cycle Benchmarking
+
+## Case and Model Description
+
+This benchmarking case is taken from [!cite](meschini2023modeling) to simulate a fuel cycle model using resident time method. The model only assign a resident time to simulate the tritium flow for each system in fuel energy plant. The model avoid the complex techinical detail and the un finished technical. The model help us to understand the challenges and potential solutions to lowering the requirement for the tritium inventory and accelerate the development of fusion energy.
+
+In the model, we have 11 systems to finish the tritium recycling in fuel cycle. The detail and corresponding ODE for each system are listed in below:
+
+!table id=tritsystems caption=Systems and labels used in this example.
+| System Name | System number | Tritium inventory variable| system equation |
+| --- | --- | --- | --- |
+| Breeding Blanket                       | 1  | `T_01_BB`       | [eqn:t1] |
+| Tritium Extraction System              | 2  | `T_02_TES`      | [eqn:t2] |
+| First Wall                             | 3  | `T_03_FW`       | [eqn:t3] |
+| Divertor                               | 4  | `T_04_DIV`      | [eqn:t4] |
+| Heat Exchanger                         | 5  | `T_05_HX`       | [eqn:t5] |
+| Detritiation System                    | 6  | `T_06_DS`       | [eqn:t6] |
+| Vaccum Pump                            | 7  | `T_07_vacuum`   | [eqn:t7] |
+| Fuel Clean-up                          | 8  | `T_08_FCU`      | [eqn:t8] |
+| Isotope Separation System              | 9  | `T_09_ISS`      | [eqn:t9] |
+| Storage and Management                 | 10 | `T_10_storage`  | [eqn:t10] |
+| Fueling System                         | 11 | -               | -         |
+| Tritium Separation Membrane            | 12 | `T_11_membrane` | [eqn:t11] |
+
+!alert note title=The fuel cycle of Tritium in Fueling System is ignored
+Fueling system only injects fresh fuel in the vacuum chamber and is not modeled in the fuel cycle to simplify the model. Instead, an outflux equal to the tritium fueling rate is added to the equation describing the storage and management system.
+
+<!-- Waiting updates -->
+
+\begin{equation}
+\label{eqn:t1}
+\frac{dI_1}{dt} = TBR \cdot \dot{N}_{T,burn} + \frac{I_3}{\tau_3} + \frac{I_4}{\tau_4} + f_{5-1} \frac{I_5}{\tau_5} - \frac{I_1}{\tau_1}- \frac{I_1\varepsilon_1}{\tau_1} - I_1\lambda ,
+\end{equation}
+
+\begin{equation}
+\label{eqn:t2}
+\frac{dI_2}{dt} = \frac{I_1}{\tau_1} - \frac{I_2}{\tau_2} - \frac{I_2\varepsilon_2}{\tau_2} - I_2\lambda ,
+\end{equation}
+
+\begin{equation}
+\label{eqn:t3}
+\frac{dI_3}{dt} = f_{p-3} \frac{ \dot{N}_{T,burn} }{TBE} + f_{5-3}\frac{I_5}{\tau_5} - \frac{I_3}{\tau_3} - \frac{I_3\varepsilon_3}{\tau_3} - I_3\lambda ,
+\end{equation}
+
+\begin{equation}
+\label{eqn:t4}
+\frac{dI_4}{dt} = f_{p-4}\frac{\dot{N}_{T,burn}}{TBE} + f_{5-4}\frac{I_5}{\tau_5} - \frac{I_4}{\tau_4} - \frac{I_4\varepsilon_4}{\tau_4}  - I_4\lambda ,
+\end{equation}
+
+\begin{equation}
+\label{eqn:t5}
+\frac{dI_5}{dt} = (1 - \eta_2)\frac{I_2}{\tau_2} - \frac{I_5}{\tau_5} -\frac{I_5\varepsilon_5}{\tau_5} - I_5\lambda ,
+\end{equation}
+
+\begin{equation}
+\label{eqn:t6}
+\frac{dI_6}{dt} = f_{5-6}\frac{I_5}{\tau_5} + f_{9-6}\frac{I_9}{\tau_9} - \frac{I_6}{\tau_6} - \frac{I_6\varepsilon_6}{\tau_6} - I_6\lambda ,
+\end{equation}
+
+\begin{equation}
+\label{eqn:t7}
+\frac{dI_7}{dt} = (1- TBE - f_{p-3} - f_{p-4}) \frac{\dot{N}_{T,burn}}{TBE} - \frac{I_7}{\tau_7} - \frac{I_7\varepsilon_7}{\tau_7} - I_7\lambda ,
+\end{equation}
+
+\begin{equation}
+\label{eqn:t8}
+\frac{dI_8}{dt} = (1 - f_{DIR}) \frac{I_7}{\tau_7} - \frac{I_8}{\tau_8} - \frac{I_8\varepsilon_8}{\tau_8} - I_8\lambda ,
+\end{equation}
+
+\begin{equation}
+\label{eqn:t9}
+\frac{dI_9}{dt} =  \frac{I_6}{\tau_6} + \frac{I_8}{\tau_8} - \frac{I_9\varepsilon_9}{\tau_9} - \frac{I_9}{\tau_9} - I_9\lambda ,
+\end{equation}
+
+\begin{equation}
+\label{eqn:t10}
+\frac{dI_{10}}{dt} = (1 - f_{9-6}) \frac{I_9}{\tau_9} + f_{DIR} \frac{I_7}{\tau_7} + \frac{I_{11}}{\tau_{11}} - \frac{\dot{N}_{T,burn}}{TBE} - I_{10}\lambda ,
+\end{equation}
+
+\begin{equation}
+\label{eqn:t11}
+\frac{dI_{11}}{dt} = \eta_2 \frac{I_2}{\tau_2} - \frac{I_{11}}{\tau_{11}} - \frac{I_{11}\varepsilon_{11}}{\tau_{11}} - I_{11}\lambda ,
+\end{equation}
+
+where $TBE$, the tritium burn efficiency, is defined as:
+
+\begin{equation}
+\label{eqn:TBE}
+TBE = \eta_f f_b
+\end{equation}
+
+
+## Case and Model Parameters
+
+We use the ScalarKernels in MOOSE to calculate the ODEs from 11 systems. All the model parameters are listed in [fuel_cycle_benchmark_table2]:
+
+!table id=fuel_cycle_benchmark_table2=Values of material properties.
+| Parameter | Description                          | Value                                                       | Units                 |
+| --------- | ------------------------------------ | ----------------------------------------------------------- | --------------------- |
+| TBR     | Tritium breeding ratio                   | 1.067                                  | -                   |
+| TBE     | Tritium burn efficiency in plasma        | 0.025                                          | -         |
+| AF | Availability factor | 0.75           | -               |
+| $\varepsilon_i$ | Fraction of tritium lost from non-radioactive phenomena in the $i$th component | 1e-4 except for FW and DIV system            | -               |
+| $\eta_f$     | Fueling efficiency        | 0.25                                                        | -                    |
+| $f_b$     | Tritium burn fraction in the plasma        | 0.10                                                        | -                    |
+| $\eta_2$     | Tritium extraction efficiency     | 0.7                                         | -              |
+| $f_{DIR}$       | Direct internal recycling fraction                         | 0.5                                       | -              |
+| $I_i$ | Tritium inventory in the $i$th component     | -                                                           | kg                     |
+| $\dot{N}_{T,burn}$ | Tritium burn rate     | 8.99e-7                                       | -                     |
+| $\lambda$ | Tritium decay rate     | 1.73e-9                                       | s$^{-1}$                     |
+| $t$ | time   | -                                                  | s                     |
+| $\tau_i$ | Tritium residence time in the $i$th component       | 4500 in $\tau_1$, 86400 in $\tau_2$, 1000 in $\tau_3$, $\tau_4$, $\tau_5$, 3600 in $\tau_6$, 600 in $\tau_7$, 585 in $\tau_8$, 22815 in $\tau_9$, 100 in $\tau_11$                                                | s                     |
+
+
+
+## Results
+
+We compare results from TMAP8 with MatLab results from [!cite](meschini2023modeling) in [comparison_fuel_cycle_benchmarking]. The agreements between the four chosen systems are quite good.
+
+!media comparison_fuel_cycle_benchmark.py
+       image_name=fuel_cycle_comparison.png
+       style=width:50%;margin-bottom:2%;margin-left:auto;margin-right:auto
+       id=fuel_cycle_comparison
+       caption=Comparison of TMAP8 calculation with the MatLab data on the fuel cycle model.
+
+## Input files
+
+!style halign=left
+The input file for this case can be found at [/fuel_cycle_benchmark/fuel_cycle.i]. The input file is different from the input file used as test in TMAP8. To limit the computational costs of the test case, the test runs a version of the file with fewer time steps. More information about the changes can be found in the test specification file for this case, namely [/fuel_cycle_benchmark/tests].
+
+
+!bibtex bibliography

doc/content/verification_and_validation/index.md

@@ -41,6 +41,7 @@ TMAP8 also contains [example cases](examples/tmap_index.md), which showcase how
 | Case    | Title                                                                              |
 | ------- | ---------------------------------------------------------------------------------- |
 | ver-1fc | [Conduction in Composite Structure with Constant Surface Temperatures](ver-1fc.md) |
+| fuel cycle | [Tritium fuel cycle in fusion energy plant](fuel_cycle_benchmarking.md) |
 
 
 # List of validation cases

test/tests/fuel_cycle_benchmark/comparison_fuel_cycle_benchmark.py

@@ -0,0 +1,116 @@
+import matplotlib.pyplot as plt
+import numpy as np
+from matplotlib import gridspec
+import pandas as pd
+from scipy import special
+import os
+
+# Changes working directory to script directory (for consistent MooseDocs usage)
+script_folder = os.path.dirname(__file__)
+os.chdir(script_folder)
+
+# ================================= Functions ================================ #
+def interpolation_on_expected_input(date_x, data_y, expected_input):
+    """Get new numerical solution based on the experimental input data points
+
+    Args:
+        expected_input (float): expected input data points
+        date_x (float, ndarray): numerical input data points
+        data_y (float, ndarray): numerical output data points
+
+    Returns:
+        float: updated expected output based on the data points in expected_input
+    """
+    left_limit = np.argwhere((np.diff(date_x < expected_input)))[0][0]
+    right_limit = left_limit + 1
+    return (expected_input - date_x[left_limit]) / (date_x[right_limit] - date_x[left_limit]) * (data_y[right_limit] - data_y[left_limit]) + data_y[left_limit]
+
+
+def read_csv_from_TMAP8(file_name, parameter_names):
+    # Read simulation data
+    if "/TMAP8/doc/" in script_folder:     # if in documentation folder
+        csv_folder = f"../../../../test/tests/fuel_cycle/gold/{file_name}"
+    else:                                  # if in test folder
+        csv_folder = f"./gold/{file_name}"
+    simulation_data = pd.read_csv(csv_folder)
+    simulation_results = []
+    for i in range(len(parameter_names)):
+        simulation_results.append(simulation_data[parameter_names[i]])
+    simulation_results = np.array(simulation_results)
+    return simulation_results
+
+# ================================ parameters ================================ #
+# reserve inventory
+time_unit = 3600 * 24 # time unit - days
+two_year = 3600 * 24 * 365 * 2 # double time
+tritium_burn_rate = 8.99e-7 # kg/s
+TBE = 0.02
+q = 0.25
+t_res = 24 * 3600 # s
+initial_inventory = 1.14 # kg
+AF = 0.7
+reserve_inventory = tritium_burn_rate / TBE * q * t_res * AF
+print(f"Required reserve inventory = {reserve_inventory} kg")
+
+# =========================== TMAP8 data extraction ========================== #
+file_name = "fuel_cycle_out.csv"
+parameter_names = ['time','T_01_BZ','T_02_TES','T_09_ISS','T_10_storage']
+simulation_results = read_csv_from_TMAP8(file_name, parameter_names) # read csv file
+simulation_results[parameter_names.index('time')] = simulation_results[parameter_names.index('time')] / time_unit # update time unit
+
+# inflection point
+inflection_y = np.min(simulation_results[parameter_names.index('T_10_storage')])
+inflection_x = simulation_results[parameter_names.index('time')][np.argmin(simulation_results[parameter_names.index('T_10_storage')])]
+print(f"TMAP8_base: Inflection time = {round(inflection_x,5)} days, and inventory = {round(inflection_y,5)} kg")
+
+# end inventory
+end_inventory = interpolation_on_expected_input(simulation_results[parameter_names.index('time')],
+                                                simulation_results[parameter_names.index('T_10_storage')],
+                                                two_year / time_unit)
+print(f"TMAP8_base: End inventory = {round(end_inventory,5)} ({round(end_inventory / initial_inventory, 3)} I_startup)")
+
+# ======================== Experiment data extraction ======================== #
+file_name = "inventory_paper.csv"
+parameter_names_experiment = ['time [s]','blanket inventory [kg]','TES inventory [kg]','ISS inventory [kg]','storage inventory [kg]']
+experiment_results = read_csv_from_TMAP8(file_name, parameter_names_experiment) # read csv file
+experiment_results[parameter_names_experiment.index('time [s]')] = experiment_results[parameter_names_experiment.index('time [s]')] / time_unit # update time unit
+
+# inflection point
+inflection_y = np.min(experiment_results[parameter_names_experiment.index('storage inventory [kg]')])
+inflection_x = experiment_results[parameter_names_experiment.index('time [s]')][np.argmin(experiment_results[parameter_names_experiment.index('storage inventory [kg]')])]
+print(f"experiment: Inflection time = {round(inflection_x,5)} days, and inventory = {round(inflection_y,5)} kg")
+# end inventory
+end_inventory = interpolation_on_expected_input(experiment_results[parameter_names_experiment.index('time [s]')],
+                                                experiment_results[parameter_names_experiment.index('storage inventory [kg]')],
+                                                two_year / time_unit)
+print(f"experiment: End inventory = {round(end_inventory,5)} ({round(end_inventory / initial_inventory, 3)} I_startup)")
+
+figure_base = 'fuel_cycle_comparison'
+# =================================== Plot =================================== #
+fig = plt.figure(figsize=[6.5, 5.5])
+gs = gridspec.GridSpec(1, 1)
+ax = fig.add_subplot(gs[0])
+
+for i in range(len(parameter_names)-1):
+    if i==0:
+        ax.plot(simulation_results[parameter_names.index('time')], simulation_results[i+1], linestyle='-', label=r"TMAP8", c='tab:grey')
+        ax.plot(experiment_results[parameter_names_experiment.index('time [s]')], experiment_results[i+1], '--', label=r"MatLab", c='tab:blue')
+    else:
+        ax.plot(simulation_results[parameter_names.index('time')], simulation_results[i+1], linestyle='-', c='tab:grey')
+        ax.plot(experiment_results[parameter_names_experiment.index('time [s]')], experiment_results[i+1], '--', c='tab:blue')
+ax.text(10, 5.5e-3, 'BZ',fontweight='bold')
+ax.text(10, 1e-1, 'TES',fontweight='bold')
+ax.text(10, 2.05e-1, 'ISS',fontweight='bold')
+ax.text(10, 9.5e-1, 'storage',fontweight='bold')
+
+ax.set_xlabel(u'time (days)')
+ax.set_ylabel(u"Tritium Inventory (kg)")
+ax.legend(loc="best",ncols=3)
+ax.set_ylim(bottom=0.001,top=1e2)
+ax.set_xlim(left=0.1)
+plt.xscale('log')
+plt.yscale('log')
+plt.grid(visible=True, which='major', color='0.65', linestyle='--', alpha=0.3)
+ax.minorticks_on()
+plt.savefig(f'{figure_base}.png', bbox_inches='tight', dpi=300)
+plt.close(fig)