{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "ff3dc3d3-f1ba-4f76-8625-8832519245bd",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import sympy as sm\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7eb9455d-669d-40aa-9c71-2146373f1ff3",
   "metadata": {},
   "source": [
    "Homework 5 - NUCE 2100\n",
    "\n",
    "**Dane Sabo**\n",
    "\n",
    "*October 1st, 2024*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1bf2efb4-39b4-4430-b544-4bd712b079ff",
   "metadata": {},
   "source": [
    "# Problem 1:\n",
    "\n",
    "The neutron flux in a bare spherical reactor of radius 50 cm is given by \n",
    "\n",
    "$$ \\phi(r) = 5x10^{13} \\frac{\\sin(0.0628r)}{r} \\frac{\\text{neutrons}}{\\text{cm}^2 \\times \\text{sec}} $$\n",
    "\n",
    "where r is measured from teh center of the reactor. The diffusion coefficient for the system is 0.80 cm.\n",
    "\n",
    "## What is the maximum value of flux in the reactor?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "d918e366-30f5-492e-90e8-b257de968b1b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7b9d11f45cd0>]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 31368 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "neutron_flux = lambda r: 5e13 * np.sin(0.0628*r)/r #neutrons/(cm^2-s)\n",
    "\n",
    "r_range = np.linspace(1e-5, 50, 100)\n",
    "flux_range = neutron_flux(r_range)\n",
    "\n",
    "plt.figure(); plt.title(\"Neutron Flux as a Function of Radius $r$\")\n",
    "plt.xlabel('Radius $r$ [cm]'); plt.ylabel('Neutron Flux $\\\\phi(r)$ $\\\\frac{\\\\text{neutrons}}{\\\\text{cm}^2 \\\\times \\\\text{sec}}$')\n",
    "plt.grid('both')\n",
    "plt.plot(r_range, flux_range, '.g')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "17b06d8d-fc16-47f0-b5c0-8bbf49198a17",
   "metadata": {},
   "source": [
    "Flux technically isn't defined at the center of the reactor, but neutron flux is greatest as r approaches 0."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "080f3f81-dc35-46fd-8cf8-dda1749e5fa9",
   "metadata": {},
   "source": [
    "## Calculate the neutron current density  as a function of position in the reactor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "263c4561-9a00-44e1-8732-a442979aed5b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle J = - \\frac{0.05024 a \\cos{\\left(0.0628 r \\right)}}{r} + \\frac{0.8 a \\sin{\\left(0.0628 r \\right)}}{r^{2}}$"
      ],
      "text/plain": [
       "Eq(J, -0.05024*a*cos(0.0628*r)/r + 0.8*a*sin(0.0628*r)/r**2)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'neutrons/cm^2/s'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "J, r, a = sm.symbols('J, r, a')\n",
    "D = 0.80 #cm\n",
    "\n",
    "psi = a*sm.sin(0.0628*r)/r\n",
    "\n",
    "#Phi only depends on r.\n",
    "#Therefore, we can take the gradient with only the partial w.r.t. r.\n",
    "diffusion_approx = sm.Eq(J, -D*psi.diff(r))\n",
    "display(diffusion_approx, 'neutrons/cm^2/s')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b5a31a40-dd14-40b0-ad0d-af8a4de4baeb",
   "metadata": {},
   "source": [
    "Where $a = 5\\times 10^{13}$ and J is in the $\\hat e_r$ direction."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e2203b1-19bf-4348-a2fb-c59410abde46",
   "metadata": {},
   "source": [
    "## How many neutrons escape from the reactor per second?\n",
    "The neutrons escaping from the reactor per second is equal to the neutron current density multiplied by the area:\n",
    "\n",
    "$$ \\int_A \\left( J(\\vec r, t) \\cdot n \\right) dA $$\n",
    "\n",
    "We can see that our neutron current density is always normal to the surface, and our surface is at a constant value of r. This lets us simplify our expression:\n",
    "\n",
    "$$ \\int_A \\left( J(\\vec r, t) \\cdot n \\right) dA \\rightarrow \n",
    "A J(r,t)|_{r = 50 \\text{[cm]}}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "8f13df8c-04c4-4833-b12e-05d24d0dd2d0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "There are 1.579e+15 neutrons leaking per second from the reactor.\n"
     ]
    }
   ],
   "source": [
    "a = 5e13\n",
    "J = lambda r: -0.05024*a*np.cos(0.0628*r)/50 + 0.8*a*np.sin(0.0628*r)/r**2\n",
    "A = lambda r: 4*np.pi*r**2\n",
    "\n",
    "print(f'There are {J(50)*A(50):.3e} neutrons leaking per second from the reactor.')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "45b572f6-2cfd-43d3-a2ee-6effafd2cbc2",
   "metadata": {},
   "source": [
    "# Problem 2\n",
    "A bare slab of thickness 2a contains uniformly distributed sources emitting Q neutrons/cm^3/sec. Given the expression for the scalar flux below verify the equation of continuity by computing per unit are of the slab the total number of neutrons\n",
    "\n",
    "$$ \\phi(x) = \\frac{Q}{\\Sigma_a} \\left( 1 - \\frac{\\cosh(x/L)}{\\cosh(a/L) } \\right) $$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "efef4bcb-e47b-4d6a-8a84-54d6c934b894",
   "metadata": {},
   "outputs": [],
   "source": [
    "Q, sigma, x, a, L = sm.symbols('Q, Sigma_a, x, a, L')\n",
    "\n",
    "phi = Q/sigma*(1-sm.cosh(x/L)/sm.cosh(a/L))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "906a04a7-def0-4eba-a8ea-d9859fa5b813",
   "metadata": {},
   "source": [
    "## Produced within the whole slab"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2ebdb98-99cf-4413-97fd-e1bcc8855322",
   "metadata": {},
   "source": [
    "$$ \n",
    "\\begin{align}\n",
    "\\int_V s(\\vec r, t) dV &= \\int_V Q dV \\\\\n",
    "&= \\int_A 2aQ dA\\\\\n",
    "\\text{but assuming per unit area...} \\\\\n",
    "&= \\frac{d}{dA} \\int_A 2aQ dA\\\\\n",
    "\\int_V s(\\vec r, t) dV &=  2aQ \\text{ per second per unit area}\\\\\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4687ccaf-dbce-4474-91a4-7fb1412298b0",
   "metadata": {},
   "source": [
    "## Absorbed per second within the slab"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "777059b4-e89d-4062-bd1e-9781067682e3",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align}\n",
    "\\int_v \\Sigma_a(\\vec r, t) \\phi(\\vec r, t) dV &= \n",
    "\\int_v \\Sigma_a(\\vec r, t) \\frac{Q}{\\Sigma_a} \\left( 1 - \\frac{\\cosh(x/L)}{\\cosh(a/L)} \\right) dV\\\\\n",
    "&= \\int_{A} \\int_x Q \\left( 1 - \\frac{\\cosh(x/L)}{\\cosh(a/L)} \\right) dxdA\\\\\n",
    "\\text{and doing things per unit area...}\\\\\n",
    "\\frac{d}{dA} \\left(\\int_v \\Sigma_a(\\vec r, t) \\phi(\\vec r, t) dV \\right) &= \\int_x Q \\left( 1 - \\frac{\\cosh(x/L)}{\\cosh(a/L)} \\right) dx\\\\\n",
    "\\end{align}\n",
    "$$\n",
    "Now we can do this integration for the thickness x using SymPy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "71824f97-740c-4ca8-a857-cffd085d3aa6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle Q \\left(- \\frac{L \\sinh{\\left(\\frac{2 a}{L} \\right)}}{\\cosh{\\left(\\frac{a}{L} \\right)}} + 2 a\\right)$"
      ],
      "text/plain": [
       "Q*(-L*sinh(2*a/L)/cosh(a/L) + 2*a)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "absorption_per_second = sm.integrate(Q*(1-sm.cosh(x/L)/sm.cosh(a/L)), (x, 0, 2*a))\n",
    "absorption_per_second"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8a5504a-08b5-4743-bde5-650a325f5c67",
   "metadata": {},
   "source": [
    "Which is our final expression for the absorption per second per unit area!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "577476c8-0bb8-46b9-896e-4681bb3bb05a",
   "metadata": {},
   "source": [
    "## Escaping per second within the slab\n",
    "$$\n",
    "\\begin{align}\n",
    "J &= -D \\nabla \\phi(\\vec r, t) \\\\\n",
    "&\\text{but since we're dealing with a large slab, all current is in the thickness (x) direction...} \\\\\n",
    "&= -D \\frac{d}{dx} \\left( \\phi(x) \\right)\\\\\n",
    "&= -\\Sigma_a L^2 \\frac{d}{dx} \\left(\\frac{Q}{\\Sigma_a} \\left( 1 - \\frac{\\cosh(x/L)}{\\cosh(a/L)}\\right)\\right)\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Now we know that the leakage rate is the neutron current out of each face. We can see by inspection that we won't have leakage out of the sides that aren't in the thickness direction. This is because we're assuming an infinitely large slab.\n",
    "\n",
    "$$\\begin{align}\n",
    "\\int_a \\left(J(x,t) \\cdot n \\right) dA &= \\int_a \\left(-\\Sigma_a L^2 \\frac{d}{dx} \\left(\\frac{Q}{\\Sigma_a} \\left( 1 - \\frac{\\cosh(x/L)}{\\cosh(a/L)}\\right)\\right) \\right)dA\\\\\n",
    "& \\text{but once again, unit area...} \\\\\n",
    "\\frac{d}{dA}\\left(\\int_a \\left(J(x,t) \\cdot n \\right) dA\\right) &= \\frac{d}{dA}\\left(\\int_a \\left(-\\Sigma_a L^2 \\frac{d}{dx} \\left(\\frac{Q}{\\Sigma_a} \\left( 1 - \\frac{\\cosh(x/L)}{\\cosh(a/L)}\\right)\\right) \\right)dA \\right)\\\\\n",
    "J(x,t) &=-\\Sigma_a L^2 \\frac{d}{dx} \\left(\\frac{Q}{\\Sigma_a} \\left( 1 - \\frac{\\cosh(x/L)}{\\cosh(a/L)}\\right)\\right) \\\\\n",
    "\\end{align}$$\n",
    "And once again just doing the manipulation in SymPy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "7ab93d22-cde3-4408-9ade-b9aeeab96ceb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{L Q \\sinh{\\left(\\frac{x}{L} \\right)}}{\\cosh{\\left(\\frac{a}{L} \\right)}}$"
      ],
      "text/plain": [
       "L*Q*sinh(x/L)/cosh(a/L)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "J = -sigma*L**2*sm.diff(Q/sigma*(1-sm.cosh(x/L)/sm.cosh(a/L)),x)\n",
    "J"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9dd7ea6e-dceb-44ec-bfd5-d183a2ee3f52",
   "metadata": {},
   "source": [
    "But! There is leakage out of two faces: leakage at x = 0, and x = 2a. The total leakage is the addition of these two."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "f1a07559-c152-4185-923e-d3885ff63dcc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 2 L Q \\sinh{\\left(\\frac{a}{L} \\right)}$"
      ],
      "text/plain": [
       "2*L*Q*sinh(a/L)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "total_leakage = J.subs(x,0) + J.subs(x,2*a)\n",
    "sm.simplify(total_leakage)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9b8927d-ce50-4848-9772-d7e674f7c97a",
   "metadata": {},
   "source": [
    "which is the total leakage of neutrons per unit area, per second."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ea81b9e6-fe62-4bad-ab6b-126d5d5b960e",
   "metadata": {},
   "source": [
    "# Problem 3\n",
    "Consider a critical bare slab reactor with the following properties: "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c0f0535-d5e6-4f8c-b12f-e8499c073d9f",
   "metadata": {},
   "source": [
    "## Calculate $k_{eff}$ for this reactor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "cb4ff251-c19a-4f89-bcb5-b11f4969322e",
   "metadata": {},
   "outputs": [],
   "source": [
    "a = 40 #cm\n",
    "Sigma_a = 0.066 #1/cm\n",
    "D = 0.9\n",
    "L = (D/Sigma_a)**(1/2)\n",
    "nu = 2.6\n",
    "nuSigma_f = 0.070 #1/cm\n",
    "e = 200 #MeV / fission\n",
    "L = (D/Sigma_a)**(1/2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d6fd9637-b0bc-43ec-a7ea-0013321db089",
   "metadata": {},
   "outputs": [],
   "source": [
    "B_g = (3.1415/a)\n",
    "k_eff = nuSigma_f/(D*B_g**2 + Sigma_a)\n",
    "k_eff"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "93c7fd3e-918d-44d7-a792-ab51ea9495ab",
   "metadata": {},
   "source": [
    "## Calculate the width, a, that yields k_eff = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "97ec39bf-c1e9-440c-a32e-cd7573bf524a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[-47.1210000000000, 47.1210000000000]"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = sm.symbols('a')\n",
    "B_g = (3.1414/a)\n",
    "k_eq = sm.Eq(1, nuSigma_f/(D*B_g**2 + Sigma_a))\n",
    "sm.solve(k_eq)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a4c2ab9-174b-4af6-8f24-c417b5326a72",
   "metadata": {},
   "source": [
    "The thickness that yields $k_{eff} = 1$ is 47.121 [cm]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58847762-ba38-4764-9a0e-f89ea2ad3a4d",
   "metadata": {},
   "source": [
    "## If the reactor is operating at 100 kW, calculate the flux distribution in the reactor, $\\Phi(x)$ using the width found in b."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "214e6142-4dee-4d85-b71f-1b0f20d1fa3c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\Phi{\\left(x \\right)} = 1.22995582046794 \\cdot 10^{15} \\pi \\cos{\\left(0.0212219604847096 \\pi x \\right)}$"
      ],
      "text/plain": [
       "Eq(Phi(x), 1.22995582046794e+15*pi*cos(0.0212219604847096*pi*x))"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Phi = sm.Function('Phi')(x)\n",
    "Sigma_f = nuSigma_f/nu\n",
    "E = e*1.60218e-13 #J\n",
    "P = 100e3 #W\n",
    "a = 47.121 #cm\n",
    "phi_eq = sm.Eq(Phi, P * sm.pi/2/a/E/Sigma_f * sm.cos(sm.pi*x/a))\n",
    "phi_eq"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0fb10451-32b4-4cc7-b74d-3d1aefa4edb2",
   "metadata": {},
   "source": [
    "Where x is in [cm]."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "92b90596-1991-41f4-84f6-b46ad711a5f3",
   "metadata": {},
   "source": [
    "## Calculate the power density as a function of x (Px)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "47a3f8c6-5fde-405c-8747-8f5591387a07",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\operatorname{Power Density}{\\left(x \\right)} = 1.22995582046794 \\cdot 10^{20} \\pi \\cos{\\left(0.0212219604847096 \\pi x \\right)}$"
      ],
      "text/plain": [
       "Eq(Power Density(x), 1.22995582046794e+20*pi*cos(0.0212219604847096*pi*x))"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Px = sm.Function('Power Density')(x)\n",
    "Power_Density_eq = sm.Eq(Px,phi_eq.rhs*P)\n",
    "Power_Density_eq"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0df4cd19-0dc2-4f0a-8490-3c9c33f273d8",
   "metadata": {},
   "source": [
    "# Problem 4\n",
    "Jezebel is a bare, fast, spherically shaped critical reactor constructed of pure $^{239}Pu$ metal ($\\rho = 15.4$ g/cm^3)\n",
    "\n",
    "## Calculate the critical radius and critical mass of the reactor using the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61536d51-6bdd-40d2-a371-f58e3e9090d2",
   "metadata": {},
   "outputs": [],
   "source": [
    "v = 2.98\n",
    "sigma_f = 1.85\n",
    "sigma_a = 2.11\n",
    "D = 1.26"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
