Update 1_find_dpa.ipynb

tasks/task_06_CSG_cell_tally_DPA/1_find_dpa.ipynb

@@ -1,7 +1,6 @@
 {
  "cells": [
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -25,7 +24,6 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -84,11 +82,10 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "This sets up the damage energy tally using the score \"damage-energy\". This is the same as MT number 444. More details on MT numbers and their full description can be found [here](https://t2.lanl.gov/nis/endf/mts.html)."
+    "This sets up the damage energy tally using the MT number 444. A list of MT numbers including their reaction discription can be found [here](https://t2.lanl.gov/nis/endf/mts.html)."
    ]
   },
   {
@@ -101,13 +98,12 @@
     "cell_filter = openmc.CellFilter(vessel_cell)\n",
     "dpa_reaction_tally = openmc.Tally(name='DPA')\n",
     "dpa_reaction_tally.filters = [cell_filter]\n",
-    "dpa_reaction_tally.scores = ['damage-energy']\n",
+    "dpa_reaction_tally.scores = ['444']  # note use of 444 in string format, this is the MT reaction number for damage energy more MT numbers here https://t2.lanl.gov/nis/endf/mts.html\n",
     "dpa_reaction_tally.nuclides = ['Fe54', 'Fe56', 'Fe57', 'Fe58'] # this records the tally for each nuclide in the list\n",
     "my_tallies = openmc.Tallies([dpa_reaction_tally])"
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -126,7 +122,6 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -150,11 +145,18 @@
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "This calculates the total number of displacements for all atoms by summing together the seperate damage-energy for each nuclide. From here, DPA can be found."
+    "The damage-energy ($T_{d}$) represents the the kinetic energy available for creating atomic displacements. This variable can be used to compute the total number of displacements per atom - or DPA. The current international standard [1] for quantifying DPA is based on the Kinchin-Pease model or NRT equation [2], which states that the number of Frenkel pairs produced $N_{d}$ in a material with displacement energy $E_{d}$ for a nuclear deposited energy $T_{d}$ is\n",
+    "\n",
+    "$$N_{d} = \\dfrac{0.8T_{d}}{2E_d},$$\n",
+    "\n",
+    "where the factor 0.8 was introduced by the authors to allow for realistic atomic scattering instead of the hard core approximation.\n",
+    "\n",
+    "[1] ASTM Standard E693-94, 'Standard practice for characterising neutron exposure in iron and low alloy steels in terms of displacements per atom (dpa)', 1994\n",
+    "\n",
+    "[2] https://doi.org/10.1016%2F0029-5493%2875%2990035-7"
    ]
   },
   {
@@ -163,29 +165,33 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "damage_energy_in_ev = df['mean'].sum()\n",
+    "damage_energy_in_ev = df['mean'].sum() #T_D in NRT equation\n",
+    "E_D = 40 # threshold displacement energy of Fe in eV\n",
     "\n",
-    "print('Damage energy deposited per source neutron = ', damage_energy_in_ev, 'eV\\n')\n",
+    "print('Damage energy deposited per source neutron = ', f\"{damage_energy_in_ev:.2f}\", 'eV\\n')\n",
     "\n",
     "print('Two times the threshold energy of 40eV is needed to displace an atom')\n",
-    "displacements_per_source_neutron = damage_energy_in_ev / (2*40)\n",
-    "print('Displacements per source neutron = ', displacements_per_source_neutron, '\\n')\n",
-    "\n",
-    "print('Assuming about 80% remains after 20% recombine to original lattice locations')\n",
-    "displacements_per_source_neutron_with_recombination = displacements_per_source_neutron*0.8\n",
-    "print('Displacements per source neutron after recombination = ', displacements_per_source_neutron_with_recombination, '\\n')\n",
-    "\n",
+    "displacements_per_source_neutron = 0.8*damage_energy_in_ev / (2*E_D)\n",
+    "print('(NRT) Displacements per source neutron = ', f\"{displacements_per_source_neutron:.2f}\", '\\n')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
     "fusion_power = 3e9  # units Watts\n",
     "energy_per_fusion_reaction = 17.6e6  # units eV\n",
     "eV_to_Joules = 1.60218e-19  # multiplication factor to convert eV to Joules\n",
     "number_of_neutrons_per_second = fusion_power / (energy_per_fusion_reaction * eV_to_Joules)\n",
-    "print('Number of neutrons per second', number_of_neutrons_per_second, '\\n')\n",
+    "print('Number of neutrons per second', f\"{number_of_neutrons_per_second:.2e}\", '\\n')\n",
     "\n",
     "number_of_neutrons_per_year = number_of_neutrons_per_second * 60 * 60 * 24 * 365.25\n",
-    "print('Number of neutrons per full power year ', number_of_neutrons_per_year)\n",
+    "print('Number of neutrons per full power year ', f\"{number_of_neutrons_per_year:.2e}\")\n",
     "\n",
-    "displacements_for_all_atoms = number_of_neutrons_per_year * displacements_per_source_neutron_with_recombination\n",
-    "print('Displacements for all atoms in the volume ', displacements_for_all_atoms, '\\n')\n",
+    "displacements_for_all_atoms = number_of_neutrons_per_year * displacements_per_source_neutron\n",
+    "print('(NRT) Displacements for all atoms in the volume ', f\"{displacements_for_all_atoms:.2e}\", '\\n')\n",
     "\n",
     "print('Now the number of atoms in the volume must be found to find displacements per atom (DPA)')\n",
     "\n",
@@ -194,15 +200,40 @@
     "iron_atomic_mass_in_g = 55.845*1.66054E-24  # molar mass multiplier by the atomic mass unit (u)\n",
     "number_of_iron_atoms = volume_of_firstwall_cell * density_of_iron_in_g_per_cm3 / (iron_atomic_mass_in_g)\n",
     "\n",
-    "print('Number of iron atoms in the firstwall ', number_of_iron_atoms)"
+    "print('Number of iron atoms in the firstwall ', f\"{number_of_iron_atoms:.2e}\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Therefore, as the total number of atoms and the total number of displacements is known, NRT-DPA can be found."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "DPA_NRT = displacements_for_all_atoms / number_of_iron_atoms\n",
+    "\n",
+    "print('NRT-DPA =', f\"{DPA_NRT:.2f}\")"
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Therefore, as the total number of atoms and the total number of displacements is known, DPA can be found."
+    "Recently, a new method [3] was proposed in 2018 to solve the limitations of NRT equation. Namely, experiments have shown that \"the number of radiation defects produced in energetic cascades in metals is only ~1/3 the NRT-DPA prediction\". The arc-DPA (athermal recombination corrected) extends the NRT-DPA approach by assuming that, after the so-called 'thermal spike', \"almost all atoms regain positions in the perfect lattice sites [...] only interstitials transported to the cascade outer periphery will result in stable defects\". \n",
+    "\n",
+    "Hence, in order to account for this behaviour, the authors modified the original NRT equation as follows:\n",
+    "\n",
+    "$$N_{d} = \\dfrac{0.8T_{D}}{2E_d} \\cdot \\xi(T_d)= \\dfrac{0.8T_{D}}{2E_d}\\left( \\dfrac{1-c}{{(2E_D/0.8)}^b} {T_D}^b + c\\right),$$\n",
+    "\n",
+    "where $\\xi$ is a surviving defect fraction factor and $b$ and $c$ are constants for a given metal. For Fe, $b=−0.568$ and $c=0.286$.\n",
+    "\n",
+    "[3] DOI: 10.1038/s41467-018-03415-5"
    ]
   },
   {
@@ -211,13 +242,25 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "DPA = displacements_for_all_atoms / number_of_iron_atoms\n",
+    "b = -0.568 # constant for Fe\n",
+    "c = 0.286\n",
+    "\n",
+    "surviving_defect_fraction_factor = (1-c)*(displacements_per_source_neutron)**b+c\n",
+    "displacements_per_source_neutron_arc = displacements_per_source_neutron*surviving_defect_fraction_factor\n",
+    "\n",
+    "print('(arc) Displacements per source neutron = ', f\"{displacements_per_source_neutron_arc:.2f}\", '\\n')\n",
+    "\n",
+    "displacements_for_all_atoms_arc = number_of_neutrons_per_year * displacements_per_source_neutron_arc\n",
+    "print('(arc) Displacements for all atoms in the volume ', f\"{displacements_for_all_atoms_arc:.2e}\", '\\n')\n",
+    "\n",
+    "arc_NRT = displacements_for_all_atoms_arc / number_of_iron_atoms\n",
     "\n",
-    "print('DPA =', DPA)"
+    "print('arc-DPA =', f\"{arc_NRT:.2f}\", '\\n')\n",
+    "\n",
+    "print('arc/NRT =', f\"{arc_NRT/DPA_NRT:.2f}\")"
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -231,12 +274,6 @@
     "- Damage energy deposited can be found with the OpenMC MT 444 tally.\n",
     "- Post tally calculations are sometimes required to convert neutronics numbers into something more useful."
    ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": []
   }
  ],
  "metadata": {
@@ -258,7 +295,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.12"
+   "version": "3.11.2"
   }
  },
  "nbformat": 4,