crp/MSRE
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions
MSRE/dynamic_model/jitcdde_implementation/simulation_study.ipynb

914 lines
717 KiB
Text
Raw Permalink Normal View History

include jitcdde implementation of dynamic model
2023-11-10 16:55:30 +01:00
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Imports"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"from jitcdde import jitcdde, y, t\n",
"from parameters import *\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import pandas as pd"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Define System"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# delays\n",
"taus = [tau_c, tau_c_hx, tau_hx_c, tau_hx_r, tau_l, tau_r_hx]\n",
"\n",
"# dT/dt for radiator nodes\n",
"T_out_rc = (W_rp/mn_rp)*(y(11,t-tau_hx_r)-y(0)) + (hA_rpn/mcp_rpn)*(y(1)-y(0)) # T_out_rc: y(0)\n",
"T_out_air = -((W_rs/mn_rs)+(hA_rsn/mcp_rsn))*y(1) + (hA_rsn/mcp_rsn)*y(0) + (W_rs/mn_rs)*Trs_in # T_out_air: y(1)\n",
"\n",
"# dT/dt for heat exchanger nodes\n",
"T_hf1 = -((W_p/mn_p)+(hA_pn/mcp_pn))*y(2) + (hA_pn/mcp_pn)*y(6) + (W_p/mn_p)*y(21,t-tau_c_hx) # T_hf1: y(2)\n",
"T_hf2 = (W_p/mn_p)*(y(2)-y(3)) + (hA_pn/mcp_pn)*(y(6)-y(2)) # T_hf2: y(3)\n",
"T_hf3 = -((W_p/mn_p)+(hA_pn/mcp_pn))*y(4) + (hA_pn/mcp_pn)*y(7) + (W_p/mn_p)*y(3) # T_hf3: y(4)\n",
"T_hf4 = (W_p/mn_p)*(y(4)-y(5)) + (hA_pn/mcp_pn)*(y(7)-y(4)) # T_hf4: y(5)\n",
"T_ht1 = (2*hA_pn/mcp_tn)*(y(2)-y(6)) + (2*hA_sn/mcp_tn)*(y(10)-y(6)) # T_ht1: y(6)\n",
"T_ht2 = (2*hA_pn/mcp_tn)*(y(4)-y(7)) + (2*hA_sn/mcp_tn)*(y(8)-y(7)) # T_ht2: y(7)\n",
"T_hc1 = -((W_s/mn_s)+(hA_sn/mcp_sn))*y(8) + (hA_sn/mcp_sn)*y(7) + (W_s/mn_s)*y(0,t-tau_r_hx) # T_hc1: y(8)\n",
"T_hc2 = (W_s/mn_s)*(y(8)-y(9)) + (hA_sn/mcp_sn)*(y(7)-y(8)) # T_hc2: y(9)\n",
"T_hc3 = -((W_s/mn_s)+(hA_sn/mcp_sn))*y(10) + (hA_sn/mcp_sn)*y(6) + (W_s/mn_s)*y(9) # T_hc3: y(10)\n",
"T_hc4 = (W_s/mn_s)*(y(10)-y(11)) + (hA_sn/mcp_sn)*(y(6)-y(10)) # T_hc4: y(11)\n",
"\n",
"# dn/dt\n",
"n = (y(22)-beta_t)*y(12)/Lam+lam[0]*y(13)+lam[1]*y(14)+lam[2]*y(15)+lam[3]*y(16)+lam[4]*y(17)+lam[5]*y(18) # n (no source insertion): y(12)\n",
"\n",
"# dC_i/dt (precursor concentrations)\n",
"C1 = y(12)*beta[0]/Lam-lam[0]*y(13)-y(13)/tau_c+y(13,t-tau_l)*np.exp(-lam[0]*tau_l)/tau_c # C1: y(13)\n",
"C2 = y(12)*beta[1]/Lam-lam[1]*y(14)-y(14)/tau_c+y(14,t-tau_l)*np.exp(-lam[1]*tau_l)/tau_c # C2: y(14)\n",
"C3 = y(12)*beta[2]/Lam-lam[2]*y(15)-y(15)/tau_c+y(15,t-tau_l)*np.exp(-lam[2]*tau_l)/tau_c # C3: y(15)\n",
"C4 = y(12)*beta[3]/Lam-lam[3]*y(16)-y(16)/tau_c+y(16,t-tau_l)*np.exp(-lam[3]*tau_l)/tau_c # C4: y(16)\n",
"C5 = y(12)*beta[4]/Lam-lam[4]*y(17)-y(17)/tau_c+y(17,t-tau_l)*np.exp(-lam[4]*tau_l)/tau_c # C5: y(17)\n",
"C6 = y(12)*beta[5]/Lam-lam[5]*y(18)-y(18)/tau_c+y(18,t-tau_l)*np.exp(-lam[5]*tau_l)/tau_c # C6: y(18)\n",
"\n",
"# dT/dt core nodes\n",
"T_cg = (hA_fg/mcp_g1)*(y(20)-y(19)) + k_g*P*y(12)/mcp_g1 # T_cg: y(19)\n",
"T_cf1 = W_f/mn_f*(y(5,t-tau_hx_c)-y(20)) + (k_f1*P*y(12)/mcp_f1) + (hA_fg*k_1*(y(19)-y(20))/mcp_f1) # T_cf1: y(20)\n",
"T_cf2 = W_f/mn_f*(y(20)-y(21)) + (k_f2*P*y(12)/mcp_f2) + (hA_fg*k_2*(y(19)-y(20))/mcp_f2) # T_cf2: y(21)\n",
"\n",
"# rho y(22)\n",
"rho = (a_f/2)*(T_cf1 + T_cf2) + (a_g)*(T_cg)\n",
"\n",
"# initial reactivity \n",
"rho_initial = 0.000"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Solve with JiTCDDE"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Generating, compiling, and loading C code.\n",
"Using default integration parameters.\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/mnt/envs/thesis_env/lib/python3.9/site-packages/jitcdde/_jitcdde.py:795: UserWarning: You did not explicitly handle initial discontinuities. Proceed only if you know what you are doing. This is only fine if you somehow chose your initial past such that the derivative of the last anchor complies with the DDE. In this case, you can set the attribute `initial_discontinuities_handled` to `True` to suppress this warning. See https://jitcdde.rtfd.io/#discontinuities for details.\n",
" warn(\"You did not explicitly handle initial discontinuities. Proceed only if you know what you are doing. This is only fine if you somehow chose your initial past such that the derivative of the last anchor complies with the DDE. In this case, you can set the attribute `initial_discontinuities_handled` to `True` to suppress this warning. See https://jitcdde.rtfd.io/#discontinuities for details.\")\n",
"/mnt/envs/thesis_env/lib/python3.9/site-packages/jitcdde/_jitcdde.py:795: UserWarning: You did not explicitly handle initial discontinuities. Proceed only if you know what you are doing. This is only fine if you somehow chose your initial past such that the derivative of the last anchor complies with the DDE. In this case, you can set the attribute `initial_discontinuities_handled` to `True` to suppress this warning. See https://jitcdde.rtfd.io/#discontinuities for details.\n",
" warn(\"You did not explicitly handle initial discontinuities. Proceed only if you know what you are doing. This is only fine if you somehow chose your initial past such that the derivative of the last anchor complies with the DDE. In this case, you can set the attribute `initial_discontinuities_handled` to `True` to suppress this warning. See https://jitcdde.rtfd.io/#discontinuities for details.\")\n",
"/mnt/envs/thesis_env/lib/python3.9/site-packages/jitcdde/_jitcdde.py:792: UserWarning: The target time is smaller than the current time. No integration step will happen. The returned state will be extrapolated from the interpolating Hermite polynomial for the last integration step. You may see this because you try to integrate backwards in time, in which case you did something wrong. You may see this just because your sampling step is small, in which case there is no need to worry.\n",
" warn(\"The target time is smaller than the current time. No integration step will happen. The returned state will be extrapolated from the interpolating Hermite polynomial for the last integration step. You may see this because you try to integrate backwards in time, in which case you did something wrong. You may see this just because your sampling step is small, in which case there is no need to worry.\")\n"
]
}
],
"source": [
"# instantiate jitcdde object\n",
"DDE = jitcdde([T_out_rc,T_out_air,T_hf1,T_hf2,T_hf3,T_hf4,T_ht1,T_ht2,T_hc1,\n",
" T_hc2,T_hc3,T_hc4,n,C1,C2,C3,C4,C5,C6,T_cg,T_cf1,T_cf2,rho])\n",
"\n",
"# set initial conditions\n",
"DDE.constant_past([T0_rp, T0_rs, T0_p1,T0_p2, T0_p3, T0_p4, T0_t1, T0_t2, T0_s1, T0_s2, \n",
" T0_s3, T0_s4, n_frac0, C0[0], C0[1], C0[2], C0[3], C0[4], C0[5], \n",
" T0_g1, T0_f1, T0_f2,rho_initial])\n",
"\n",
"# jitcdde solver parameters \n",
"t0 = 0.0\n",
"tf = 1000.00\n",
"T = np.arange(t0,tf,0.01)\n",
"\n",
"sol_jit = []\n",
"for t_x in T:\n",
" sol_jit.append(DDE.integrate(t_x))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Vary temeprature feedback coefficients"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Generating, compiling, and loading C code.\n",
"Using default integration parameters.\n",
"Generating, compiling, and loading C code.\n",
"Using default integration parameters.\n",
"Generating, compiling, and loading C code.\n",
"Using default integration parameters.\n",
"Generating, compiling, and loading C code.\n",
"Using default integration parameters.\n"
]
}
],
"source": [
"# solution vectors\n",
"sols = []\n",
"sols = [sol_jit]\n",
"\n",
"# feedbacks\n",
"factors_f = [1/2, 1/4, 1/8, 0]\n",
"factors_g = [1/2, 1/4, 1/8, 0]\n",
"\n",
"for i in range(len(factors_f)):\n",
"\n",
" # reinstantiate jitcdde objects\n",
" import sys\n",
" sys.modules.pop('jitcdde')\n",
" from jitcdde import jitcdde, y, t\n",
"\n",
" # rho y(22)\n",
" a_f_i = factors_f[i]\n",
" a_g_i = factors_g[i]\n",
" rho = a_f_i*(a_f/2)*(T_cf1 + T_cf2) + a_g_i*(a_g)*(T_cg)\n",
"\n",
" # instantiate jitcdde object\n",
" DDE = jitcdde([T_out_rc,T_out_air,T_hf1,T_hf2,T_hf3,T_hf4,T_ht1,T_ht2,T_hc1,\n",
" T_hc2,T_hc3,T_hc4,n,C1,C2,C3,C4,C5,C6,T_cg,T_cf1,T_cf2,rho])\n",
"\n",
" # set initial conditions\n",
" DDE.constant_past([T0_rp, T0_rs, T0_p1,T0_p2, T0_p3, T0_p4, T0_t1, T0_t2, T0_s1, T0_s2, \n",
" T0_s3, T0_s4, n_frac0, C0[0], C0[1], C0[2], C0[3], C0[4], C0[5], \n",
" T0_g1, T0_f1, T0_f2,rho_initial])\n",
" \n",
" sol_i = []\n",
" for t_x in T:\n",
" sol_i.append(DDE.integrate(t_x))\n",
" sols.append(sol_i)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Text(0.5, 1.0, '$\\\\rho$')"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAABcIAAAPxCAYAAADZsN7tAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOzdeZxN9f/A8dddZ1+ZMTOyjCWGDEJiCFmjVEhEDZVpk+IbUtnKkkoo35TqK5U2hZ98vyprKZLIVvbdmI3Zt7ue3x937jV3FmaYO3dmvJ+P7uPe8znb+1zTPfe87+e8PypFURSEEEIIIYQQQgghhBBCiBpK7e4AhBBCCCGEEEIIIYQQQghXkkS4EEIIIYQQQgghhBBCiBpNEuFCCCGEEEIIIYQQQgghajRJhAshhBBCCCGEEEIIIYSo0SQRLoQQQgghhBBCCCGEEKJGk0S4EEIIIYQQQgghhBBCiBpNEuFCCCGEEEIIIYQQQgghajRJhAshhBBCCCGEEEIIIYSo0SQRLoQQQgghhBBCCCGEEKJGk0S4uCGpVCpmzJjh7jCKGTVqFA0bNnR3GKIa+eOPP9Dr9Zw5c+aa1v/hhx/w9fUlJSWlgiMTQghnco6D06dPo1Kp+OSTT9wdiqhG+vfvz5gxY655/dtvv51JkyZVYERCiJpq69atqFQqvv32W3eHUuHke4goL7nWrpkkES7c6pNPPkGlUpX4ePHFF90am/1iVaVS8d133xWbP2PGDFQqFRcvXnRDdFdWOParPU6fPu3ucN1mzpw5rFmzxt1hXJeXX36Z4cOH06BBg2LzVq9ezV133UXt2rXR6/VEREQwdOhQNm/e7FimX79+NGnShLlz51Zm2EIINzpx4gRPPPEEjRo1wtPTE39/f2JiYli0aBF5eXnuDu+6bN++nRkzZpCenl6m5UeNGoVKpSI6OhpFUYrNV6lUjB07toKjrBj22K/2GDVqlLtDdZvy/j1URb/99hs//fQTkydPLjYvKSmJF154gebNm+Pt7Y2Pjw/t2rVj1qxZTsc8efJk/v3vf5OYmFiJkQshKktZr/u2bt3q7lAd5Fq75pNrbbnWrqq07g5ACIBXX32VyMhIp7ZbbrnFTdEU9+qrrzJo0CBUKpW7QymTkJAQPvvsM6e2+fPnc/78eRYsWFBs2RvVnDlzGDJkCPfdd5+7Q7kme/fuZePGjWzfvt2pXVEUHn30UT755BPatm3LhAkTCAsLIyEhgdWrV9OzZ09+++03OnfuDMATTzzBCy+8wMyZM/Hz83PHoQghKsl///tfHnjgATw8PHjkkUe45ZZbMBqN/Prrr0ycOJG///6bpUuXujvMa7Z9+3ZmzpzJqFGjCAwMLPN6Bw4cYNWqVQwePNh1wVWwJ554gl69ejmmT506xbRp04iLi6Nr166O9saNG7sjvCrhWv8eqpI333yTnj170qRJE6f2Xbt20b9/f7Kzsxk5ciTt2rUD4M8//+T111/nl19+4aeffgLg3nvvxd/fn/fee49XX3210o9BCOFaRa/7Pv30UzZs2FCsPSoqqjLDKjO51q6Z5FpbrrWrKkmEiyrhrrvuon379u4Oo0Rt2rRh7969rF69mkGDBrk7nDLx8fFh5MiRTm1fffUVaWlpxdprCqvVitFoxNPT84aJY9myZdSvX5/bb7/dqX3+/Pl88sknPP/887z99ttOXypffvllPvvsM7Tayx//gwcP5tlnn2XlypU8+uijLo9bCOEep06dYtiwYTRo0IDNmzcTHh7umPfMM89w/Phx/vvf/7oxQvfw8vKiXr161e5CvFOnTnTq1Mkx/eeffzJt2jQ6depUY8/1OTk5+Pj4uDuMSosjOTmZ//73v7z//vtO7enp6dx///1oNBr++usvmjdv7jR/9uzZfPjhh45ptVrNkCFD+PTTT5k5c2a1+RsXQpRN0c/833//nQ0bNlSLc4Fca1cPcq19mVxrV39SGkVUeaXV827YsGGx233T09N5/vnnqVevHh4eHjRp0oR58+ZhtVqvef/Dhg3j5ptv5tVXXy3xtumiVq5cSbt27fDy8qJ27dqMHDmS+Pj4YsutWbOGW265BU9PT2655RZWr15d4vasVisLFy6kZcuWeHp6UqdOHZ544gnS0tKu+ZjsDAYD06dPp0mTJnh4eFCvXj0mTZqEwWBwWs5+a/jKlStp0aIFXl5edOrUiQMHDgDwwQcf0KRJEzw9PenevXuxW8C6d+/OLbfcwu7du+ncuTNeXl5ERkYWu7C7lphWrFhBy5Yt8fDw4IcffgDgrbfeonPnztSqVQsvLy/atWtXrM6dSqUiJyeH5cuXF7t9vLT6cfZb9MoaR3x8PI8++ih16tTBw8ODli1b8p///KfYdt99911atmyJt7c3QUFBtG/fni+++KLYckWtWbOGO++80ymmvLw85s6dS/PmzXnrrbdKvNh9+OGHue222xzToaGhREdH83//939X3acQovp64403yM7O5uOPP3ZKgts1adKE5557zjFtNpt57bXXaNy4MR4eHjRs2JCXXnqp2OcxwHvvvef4DIyIiOCZZ54pUzmKsnxew+XPWvu50/6Zav+8Bdtn9MSJEwGIjIws823JarWaV155hf3795d6Li4sOTmZxx57jDp16uDp6Unr1q1Zvnx5seXS09MZNWoUAQEBBAYGEhsbW+p7cvjwYYYMGUJwcDCenp60b9+etWvXXjWWsti5cyf9+vUjICAAb29vunXrxm+//ea0jP38dvToUUaOHElAQAAhISFMnToVRVE4d+6co1dxWFgY8+fPd1rfXlP266+/5qWXXiIsLAwfHx8GDhzIuXPnriumf/75h4ceeoigoCC6dOkCwP79+xk1apSjvE9YWBiPPvooly5dclq/tL+HK9VqL/q980pxAHz++eeO733BwcEMGzas2DEfO3aMwYMHExYWhqenJzfddBPDhg0jIyPjCv9ytjs4zGazU89/sH3vio+P5+233y6WBAeoU6cOr7zyilNb7969OXPmDHv37r3iPoUQNVNJ185gu07r3r17sXaLxVJhn+dXItfacq1tJ9faorJIj3BRJWRkZBSr/1W7du1ybSM3N5du3boRHx/PE088Qf369dm+fTtTpkwhISGBhQsXXlNsGo2GV155hUceeeSqv1R/8sknjB49mg4dOjB37lySkpJYtGgRv/32G3/99ZfjttyffvqJwYMH06JFC+bOnculS5cYPXo0N910U7FtPvHEE47tjhs3jlOnTrF48WL++usvfvvtN3Q63TUdl9VqZeDAgfz666/ExcURFRXFgQMHWLBgAUePHi1Wz2vbtm2sXbuWZ555BoC5c+dy9913M2nSJN577z2efvpp0tLSeOONN3j00UedamMBpKWl0b9/f4YOHcrw4cP55ptveOqpp9Dr9Y5fRssb0+bNm/nmm28YO3YstWvXdpxQFy1axMCBAxkxYgRGo5GvvvqKBx54gHXr1jFgwADAdgvh448/zm233UZcXBxw7bePlxRHUlISt99+u+PkHRISwvr163nsscfIzMzk+eefB+DDDz9k3LhxDBkyhOeee478/Hz279/Pzp07eeihh0rdZ3x8PGfPnuXWW291av/1119JTU3l+eefR6PRlPkY2rVrV+1ruAkhruz777+nUaNGjls1r+bxxx9n+fLlDBkyhH/961/s3LmTuXPncujQIacLyhkzZjBz5kx69erFU089xZEjR1iyZAm7du266nmqLJ/Xdr/++iurVq3i6aefxs/Pj3feeYfBgwdz9uxZatWqxaBBgzh69ChffvklCxYscHyPKMttyQ899BCvvfYar776Kvfff3+pPWbz8vLo3r07x48fZ+zYsURGRrJy5UpGjRpFenq644cERVG49957+fXXX3nyySeJiopi9erVxMbGFtvm33//TUxMDHXr1uXFF1/Ex8eHb775hvvuu4/vvvuO+++//6rxl2bz5s3cddddtGvXjunTp6NWq1m2bBl33nkn27Ztc7pQA3jwwQeJiori9ddf57///S+zZs0iODiYDz74gDvvvJN58+axYsUKXnj
"text/plain": [
"<Figure size 1800x1200 with 6 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig,axs = plt.subplots(2,3,figsize=(18,12))\n",
"\n",
"# fuel temps\n",
"#axs[0,0].set_xlim([0,20])\n",
"axs[0,0].plot(T,[s[20] for s in sols[0]],label=\"core 1\") \n",
"axs[0,0].plot(T,[s[21] for s in sols[0]],label=\"core 2\") \n",
"axs[0,0].plot(T,[s[2] for s in sols[0]],label=\"hx 1\") \n",
"axs[0,0].plot(T,[s[3] for s in sols[0]],label=\"hx 2\")\n",
"axs[0,0].plot(T,[s[4] for s in sols[0]],label=\"hx 3\")\n",
"axs[0,0].plot(T,[s[5] for s in sols[0]],label=\"hx 4\") \n",
"axs[0,0].legend()\n",
"axs[0,0].set_title(\"Fuel Node Temperatures (C)\")\n",
"axs[0,0].tick_params(\n",
" axis='x', # changes apply to the x-axis\n",
" which='both', # both major and minor ticks are affected\n",
" bottom=False, # ticks along the bottom edge are off\n",
" top=False, # ticks along the top edge are off\n",
" labelbottom=False) # labels along the bottom edge are off\n",
"\n",
"# coolant temps\n",
"#axs[0,1].set_xlim([0,20])\n",
"axs[0,1].plot(T,[s[8] for s in sols[0]],label=\"hx 1\") \n",
"axs[0,1].plot(T,[s[9] for s in sols[0]],label=\"hx 2\") \n",
"axs[0,1].plot(T,[s[10] for s in sols[0]],label=\"hx 3\") \n",
"axs[0,1].plot(T,[s[11] for s in sols[0]],label=\"hx 4\")\n",
"axs[0,1].plot(T,[s[0] for s in sols[0]],label=\"r 1\")\n",
"axs[0,1].legend()\n",
"axs[0,1].set_title(\"Coolant Node Temperatures (C)\")\n",
"axs[0,1].tick_params(\n",
" axis='x', # changes apply to the x-axis\n",
" which='both', # both major and minor ticks are affected\n",
" bottom=False, # ticks along the bottom edge are off\n",
" top=False, # ticks along the top edge are off\n",
" labelbottom=False) # labels along the bottom edge are off\n",
"\n",
"# tube node temps\n",
"#axs[0,2].set_xlim([0,20])\n",
"axs[0,2].plot(T,[s[6] for s in sols[0]],label=\"hx 1\") \n",
"axs[0,2].plot(T,[s[7] for s in sols[0]],label=\"hx 2\") \n",
"axs[0,2].legend()\n",
"axs[0,2].set_title(\"Tube Node Temperatures (C)\")\n",
"axs[0,2].tick_params(\n",
" axis='x', # changes apply to the x-axis\n",
" which='both', # both major and minor ticks are affected\n",
" bottom=False, # ticks along the bottom edge are off\n",
" top=False, # ticks along the top edge are off\n",
" labelbottom=False) # labels along the bottom edge are off\n",
"\n",
"# precursor concentrations\n",
"#axs[1,2].set_xlim([0,20])\n",
"axs[1,2].plot(T,[s[13] for s in sols[0]],label=\"C1\") \n",
"axs[1,2].plot(T,[s[14] for s in sols[0]],label=\"C2\") \n",
"axs[1,2].plot(T,[s[15] for s in sols[0]],label=\"C3\") \n",
"axs[1,2].plot(T,[s[16] for s in sols[0]],label=\"C4\")\n",
"axs[1,2].plot(T,[s[17] for s in sols[0]],label=\"C5\")\n",
"axs[1,2].plot(T,[s[18] for s in sols[0]],label=\"C6\")\n",
"axs[1,2].legend()\n",
"axs[1,2].set_xlabel(\"t (s)\")\n",
"axs[1,2].set_yscale(\"log\")\n",
"axs[1,2].set_ylabel(r\"concentration (1/cm$^3$)\")\n",
"axs[1,2].legend(loc=\"right\")\n",
"axs[1,2].set_title(\"Precursor Concentrations\")\n",
"\n",
"# multiplication factor temp\n",
"#axs[1,0].set_xlim([0,20])\n",
"axs[1,0].plot(T,[s[12] for s in sols[0]],label=\"n\") \n",
"axs[1,0].set_xlabel(\"t (s)\")\n",
"axs[1,0].set_title(r\"$n$\")\n",
"axs[1,0].set_ylabel(r\"$\\frac{n}{n_0}$\")\n",
"\n",
"# reactivity\n",
"#axs[1,1].set_xlim([0,20])\n",
"axs[1,1].plot(T,[s[22] for s in sols[0]],label=\"n\") \n",
"axs[1,1].set_xlabel(\"t (s)\")\n",
"axs[1,1].set_title(r\"$\\rho$\")\n",
"\n",
"#plt.figure(figsize=(8,6))\n",
"#plt.plot(T,[s[12] for s in sols[0]],label=r\"$(\\alpha_f,\\alpha_g) = $ (-8.71,-6.66) pcm\")\n",
"#plt.plot(T,[s[12] for s in sols[1]],label=r\"$(\\alpha_f,\\alpha_g) = $ (-4.36,-3.33) pcm\")\n",
"#plt.plot(T,[s[12] for s in sols[2]],label=r\"$(\\alpha_f,\\alpha_g) = $ (-2.18,-1.67) pcm\")\n",
"#plt.plot(T,[s[12] for s in sols[3]],label=r\"$(\\alpha_f,\\alpha_g) = $ (-1.09,-0.83) pcm\")\n",
"#plt.plot(T,[s[12] for s in sols[4]],label=r\"$(\\alpha_f,\\alpha_g) = $ (0.0,0.0) pcm\")\n",
"#plt.xlabel(\"t (s)\")\n",
"#plt.ylabel(r\"$n/n_0$\",rotation=0)\n",
"#plt.legend()\n",
"#plt.ylim([0,20])\n",
"#plt.xlim([0,500])\n",
"#plt.title(r\"Reactor Response vs $\\alpha$\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Feedback"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Text(0.5, 1.0, '$\\\\rho$')"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAABcIAAAPxCAYAAADZsN7tAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOzdd7wcVf3/8ddsb3d3by/JTXJzEzoYDIoICApSFVQQUZAAUkSUn4UiXwsEFAREQSyA36+Iig0QvohfQKogHSHUJCQkubm93+1tZs7vj9nde/eW1Jvcks/z8VhmdmZ29uzksmfmvWfO0ZRSCiGEEEIIIYQQQgghhBBilrJNdQGEEEIIIYQQQgghhBBCiB1JgnAhhBBCCCGEEEIIIYQQs5oE4UIIIYQQQgghhBBCCCFmNQnChRBCCCGEEEIIIYQQQsxqEoQLIYQQQgghhBBCCCGEmNUkCBdCCCGEEEIIIYQQQggxq0kQLoQQQgghhBBCCCGEEGJWkyBcCCGEEEIIIYQQQgghxKwmQbgQQgghhBBCCCGEEEKIWU2CcLFL0jSNK6+8cqqLMcaZZ57JggULproYYgZ56aWXcLlctLS0bNPrH374YQKBAL29vZNcMiGEKCV1HGzYsAFN0/jtb3871UURM8hxxx3Hueeeu82v/9CHPsSll146iSUSQsxWTz31FJqmcc8990x1USadnIeIrSXX2rOTBOFiSv32t79F07RxH9/+9rentGyFi1VN07j33nvHrL/yyivRNI2+vr4pKN2mjSz75h4bNmyY6uJOmWuuuYb7779/qouxXb7zne/w+c9/nvnz549Zd99993HsscdSVVWFy+WioaGBU045hSeeeKK4zTHHHMOiRYu49tprd2axhRBT6L333uP8889n4cKFeDwegsEgBx98MDfffDOpVGqqi7ddnnvuOa688kqGhoa2aPszzzwTTdPYb7/9UEqNWa9pGl/96lcnuZSTo1D2zT3OPPPMqS7qlNnav4fp6Nlnn+Wf//wnl1122Zh13d3dXHzxxeyxxx74fD78fj9Lly7lBz/4Qclnvuyyy/jFL35BV1fXTiy5EGJn2dLrvqeeemqqi1ok19qzn1xry7X2dOWY6gIIAXDVVVfR1NRUsmyfffaZotKMddVVV/GZz3wGTdOmuihbpLq6mt///vcly2688Uba2tr46U9/OmbbXdU111zDySefzKc+9ampLso2WbFiBY899hjPPfdcyXKlFGeffTa//e1v2X///fnmN79JXV0dnZ2d3HfffRxxxBE8++yzfPjDHwbg/PPP5+KLL2b58uWUlZVNxUcRQuwk//jHP/jsZz+L2+3mjDPOYJ999iGbzfLvf/+bSy65hLfffpvbb799qou5zZ577jmWL1/OmWeeSTgc3uLXvfnmm/ztb3/jpJNO2nGFm2Tnn38+Rx55ZPH5+vXr+f73v895553HoYceWlze3Nw8FcWbFrb172E6ueGGGzjiiCNYtGhRyfKXX36Z4447jng8zumnn87SpUsBeOWVV/jRj37E008/zT//+U8ATjzxRILBIL/85S+56qqrdvpnEELsWKOv+373u9/x6KOPjlm+55577sxibTG51p6d5FpbrrWnKwnCxbRw7LHHcsABB0x1Mca1ZMkSVqxYwX333cdnPvOZqS7OFvH7/Zx++ukly/785z8zODg4ZvlsYZom2WwWj8ezy5TjjjvuYN68eXzoQx8qWX7jjTfy29/+lq9//ev85Cc/KTmp/M53vsPvf/97HI7hr/+TTjqJr33ta9x9992cffbZO7zcQoipsX79ek499VTmz5/PE088QX19fXHdhRdeyNq1a/nHP/4xhSWcGl6vl8bGxhl3IX7QQQdx0EEHFZ+/8sorfP/73+eggw6atXV9IpHA7/dPdTF2Wjl6enr4xz/+wa233lqyfGhoiE9/+tPY7XZee+019thjj5L1P/zhD/n1r39dfG6z2Tj55JP53e9+x/Lly2fM37gQYsuM/s5/4YUXePTRR2dEXSDX2jODXGsPk2vtmU+6RhHT3kT9eS9YsGDM7b5DQ0N8/etfp7GxEbfbzaJFi7juuuswTXOb3//UU09lt91246qrrhr3tunR7r77bpYuXYrX66WqqorTTz+d9vb2Mdvdf//97LPPPng8HvbZZx/uu+++cfdnmiY33XQTe++9Nx6Ph9raWs4//3wGBwe3+TMVZDIZrrjiChYtWoTb7aaxsZFLL72UTCZTsl3h1vC7776bvfbaC6/Xy0EHHcSbb74JwG233caiRYvweDwcfvjhY24BO/zww9lnn334z3/+w4c//GG8Xi9NTU1jLuy2pUx33XUXe++9N263m4cffhiAH//4x3z4wx+msrISr9fL0qVLx/Rzp2kaiUSCO++8c8zt4xP1H1e4RW9Ly9He3s7ZZ59NbW0tbrebvffem9/85jdj9nvLLbew99574/P5KC8v54ADDuCPf/zjmO1Gu//++/nYxz5WUqZUKsW1117LHnvswY9//ONxL3a/+MUv8sEPfrD4vKamhv3224///d//3ex7CiFmruuvv554PM7//M//lITgBYsWLeL//b//V3yu6zpXX301zc3NuN1uFixYwH/913+N+T4G+OUvf1n8DmxoaODCCy/cou4otuT7Goa/awt1Z+E7tfB9C9Z39CWXXAJAU1PTFt+WbLPZ+O53v8sbb7wxYV08Uk9PD1/60peora3F4/Hwvve9jzvvvHPMdkNDQ5x55pmEQiHC4TDLli2b8JisWrWKk08+mYqKCjweDwcccAAPPPDAZsuyJV588UWOOeYYQqEQPp+Pww47jGeffbZkm0L99u6773L66acTCoWorq7me9/7HkopWltbi62K6+rquPHGG0teX+hT9i9/+Qv/9V//RV1dHX6/nxNOOIHW1tbtKtM777zDF77wBcrLyznkkEMAeOONNzjzzDOL3fvU1dVx9tln09/fX/L6if4eNtVX++jzzk2VA+APf/hD8byvoqKCU089dcxnXrNmDSeddBJ1dXV4PB7mzp3LqaeeSiQS2cS/nHUHh67rJS3/wTrvam9v5yc/+cmYEBygtraW7373uyXLPv7xj9PS0sKKFSs2+Z5CiNlpvGtnsK7TDj/88DHLDcOYtO/zTZFrbbnWLpBrbbGzSItwMS1EIpEx/X9VVVVt1T6SySSHHXYY7e3tnH/++cybN4/nnnuOyy+/nM7OTm666aZtKpvdbue73/0uZ5xxxmZ/qf7tb3/LWWedxQc+8AGuvfZauru7ufnmm3n22Wd57bXXirfl/vOf/+Skk05ir7324tprr6W/v5+zzjqLuXPnjtnn+eefX9zvRRddxPr16/n5z3/Oa6+9xrPPPovT6dymz2WaJieccAL//ve/Oe+889hzzz158803+elPf8q77747pj+vZ555hgceeIALL7wQgGuvvZZPfOITXHrppfzyl7/kK1/5CoODg1x//fWcffbZJX1jAQwODnLcccdxyimn8PnPf56//vWvXHDBBbhcruIvo1tbpieeeIK//vWvfPWrX6WqqqpYod58882ccMIJnHbaaWSzWf785z/z2c9+lgcffJDjjz8esG4hPOecc/jgBz/IeeedB2z77ePjlaO7u5sPfehDxcq7urqahx56iC996UtEo1G+/vWvA/DrX/+aiy66iJNPPpn/9//+H+l0mjfeeIMXX3yRL3zhCxO+Z3t7Oxs3buT9739/yfJ///vfDAwM8PWvfx273b7Fn2Hp0qUzvg83IcSm/f3vf2fhwoXFWzU355xzzuHOO+/k5JNP5lvf+hYvvvgi1157LStXriy5oLzyyitZvnw5Rx55JBdccAGrV6/mV7/6FS+//PJm66kt+b4u+Pe//83f/vY3vvKVr1BWVsbPfvYzTjrpJDZu3EhlZSWf+cxnePfdd/nTn/7ET3/60+J5xJbclvyFL3yBq6++mquuuopPf/rTE7aYTaVSHH744axdu5avfvWrNDU1cffdd3PmmWcyNDRU/CFBKcWJJ57Iv//9b7785S+z5557ct9997Fs2bIx+3z77bc5+OC
"text/plain": [
"<Figure size 1800x1200 with 6 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig,axs = plt.subplots(2,3,figsize=(18,12))\n",
"\n",
"t_stop = 125\n",
"\n",
"# fuel temps\n",
"#axs[0,0].set_xlim([0, 500])\n",
"#axs[0,0].set_ylim([600, 1000])\n",
"axs[0,0].set_xlim([0,t_stop])\n",
"axs[0,0].plot(T,[s[20] for s in sols[0]],label=r\"core 1: $(\\alpha_f,\\alpha_g) = $ (-8.71,-6.66) pcm\")\n",
"axs[0,0].plot(T,[s[20] for s in sols[1]],label=r\"core 1: $(\\alpha_f,\\alpha_g) = $ (-4.36,-3.33) pcm\")\n",
"axs[0,0].plot(T,[s[20] for s in sols[2]],label=r\"core 1: $(\\alpha_f,\\alpha_g) = $ (-2.18,-1.67) pcm\")\n",
"axs[0,0].plot(T,[s[20] for s in sols[3]],label=r\"core 1: $(\\alpha_f,\\alpha_g) = $ (-1.09,-0.83) pcm\")\n",
"axs[0,0].plot(T,[s[20] for s in sols[4]],label=r\"core 1: $(\\alpha_f,\\alpha_g) = $ (-0.00,-0.00) pcm\")\n",
"axs[0,0].legend()\n",
"axs[0,0].set_title(\"Fuel Node Temperatures (C)\")\n",
"axs[0,0].tick_params(\n",
" axis='x', # changes apply to the x-axis\n",
" which='both', # both major and minor ticks are affected\n",
" bottom=False, # ticks along the bottom edge are off\n",
" top=False, # ticks along the top edge are off\n",
" labelbottom=False) # labels along the bottom edge are off\n",
"\n",
"# coolant temps\n",
"#axs[0,1].set_xlim([0, 1000])\n",
"#axs[0,1].set_ylim([500, 1000])\n",
"axs[0,1].set_xlim([0,t_stop])\n",
"axs[0,1].plot(T,[s[8] for s in sols[0]],label=r\"hx 1: $(\\alpha_f,\\alpha_g) = $ (-8.71,-6.66) pcm\")\n",
"axs[0,1].plot(T,[s[8] for s in sols[1]],label=r\"hx 1: $(\\alpha_f,\\alpha_g) = $ (-4.36,-3.33) pcm\")\n",
"axs[0,1].plot(T,[s[8] for s in sols[2]],label=r\"hx 1: $(\\alpha_f,\\alpha_g) = $ (-2.18,-1.67) pcm\")\n",
"axs[0,1].plot(T,[s[8] for s in sols[3]],label=r\"hx 1: $(\\alpha_f,\\alpha_g) = $ (-1.09,-0.83) pcm\")\n",
"axs[0,1].plot(T,[s[8] for s in sols[4]],label=r\"hx 1: $(\\alpha_f,\\alpha_g) = $ (-0.00,-0.00) pcm\")\n",
"axs[0,1].legend()\n",
"axs[0,1].set_title(\"Coolant Node Temperatures (C)\")\n",
"axs[0,1].tick_params(\n",
" axis='x', # changes apply to the x-axis\n",
" which='both', # both major and minor ticks are affected\n",
" bottom=False, # ticks along the bottom edge are off\n",
" top=False, # ticks along the top edge are off\n",
" labelbottom=False) # labels along the bottom edge are off\n",
"\n",
"# tube node temps\n",
"#axs[0,2].set_xlim([0, 500])\n",
"#axs[0,2].set_ylim([500, 1000])\n",
"axs[0,2].set_xlim([0,t_stop])\n",
"axs[0,2].plot(T,[s[8] for s in sols[0]],label=r\"t1: $(\\alpha_f,\\alpha_g) = $ (-8.71,-6.66) pcm\")\n",
"axs[0,2].plot(T,[s[8] for s in sols[1]],label=r\"t1: $(\\alpha_f,\\alpha_g) = $ (-4.36,-3.33) pcm\")\n",
"axs[0,2].plot(T,[s[8] for s in sols[2]],label=r\"t1: $(\\alpha_f,\\alpha_g) = $ (-2.18,-1.67) pcm\")\n",
"axs[0,2].plot(T,[s[8] for s in sols[3]],label=r\"t1: $(\\alpha_f,\\alpha_g) = $ (-1.09,-0.83) pcm\")\n",
"axs[0,2].plot(T,[s[8] for s in sols[4]],label=r\"t1: $(\\alpha_f,\\alpha_g) = $ (-0.00,-0.00) pcm\")\n",
"\n",
"# tube node temps\n",
"axs[0,2].plot(T,[s[9] for s in sols[0]],label=r\"t2: $(\\alpha_f,\\alpha_g) = $ (-8.71,-6.66) pcm\")\n",
"axs[0,2].plot(T,[s[9] for s in sols[1]],label=r\"t2: $(\\alpha_f,\\alpha_g) = $ (-4.36,-3.33) pcm\")\n",
"axs[0,2].plot(T,[s[9] for s in sols[2]],label=r\"t2: $(\\alpha_f,\\alpha_g) = $ (-2.18,-1.67) pcm\")\n",
"axs[0,2].plot(T,[s[9] for s in sols[3]],label=r\"t2: $(\\alpha_f,\\alpha_g) = $ (-1.09,-0.83) pcm\")\n",
"axs[0,2].plot(T,[s[9] for s in sols[4]],label=r\"t2: $(\\alpha_f,\\alpha_g) = $ (-0.00,-0.00) pcm\")\n",
"\n",
"axs[0,2].legend()\n",
"axs[0,2].set_title(\"Tube Node Temperatures (C)\")\n",
"axs[0,2].tick_params(\n",
" axis='x', # changes apply to the x-axis\n",
" which='both', # both major and minor ticks are affected\n",
" bottom=False, # ticks along the bottom edge are off\n",
" top=False, # ticks along the top edge are off\n",
" labelbottom=False) # labels along the bottom edge are off\n",
"\n",
"# precursor concentrations\n",
"#axs[1,2].set_xlim([0, 500])\n",
"#axs[1,2].set_ylim([1e1,1e3])\n",
"axs[1,2].set_xlim([0,t_stop])\n",
"axs[1,2].plot(T,[s[13] for s in sols[0]],label=r\"C1: $(\\alpha_f,\\alpha_g) = $ (-8.71,-6.66) pcm\")\n",
"axs[1,2].plot(T,[s[13] for s in sols[1]],label=r\"C1: $(\\alpha_f,\\alpha_g) = $ (-4.36,-3.33) pcm\")\n",
"axs[1,2].plot(T,[s[13] for s in sols[2]],label=r\"C1: $(\\alpha_f,\\alpha_g) = $ (-2.18,-1.67) pcm\")\n",
"axs[1,2].plot(T,[s[13] for s in sols[3]],label=r\"C1: $(\\alpha_f,\\alpha_g) = $ (-1.09,-0.83) pcm\")\n",
"axs[1,2].plot(T,[s[13] for s in sols[4]],label=r\"C1: $(\\alpha_f,\\alpha_g) = $ (-0.00,-0.00) pcm\")\n",
"axs[1,2].legend()\n",
"axs[1,2].set_xlabel(\"t (s)\")\n",
"axs[1,2].set_yscale(\"log\")\n",
"axs[1,2].legend(loc=\"right\")\n",
"axs[1,2].set_title(\"Precursor Concentrations\")\n",
"axs[1,2].set_ylabel(r\"concentration (1/cm$^3$)\")\n",
"\n",
"# multiplication factor \n",
"#axs[1,0].set_ylim([0,20])\n",
"#axs[1,0].set_xlim([0,500])\n",
"axs[1,0].set_xlim([0,t_stop])\n",
"axs[1,0].plot(T,[s[12] for s in sols[0]],label=r\"n: $(\\alpha_f,\\alpha_g) = $ (-8.71,-6.66) pcm\")\n",
"axs[1,0].plot(T,[s[12] for s in sols[1]],label=r\"n: $(\\alpha_f,\\alpha_g) = $ (-8.71,-6.66) pcm\")\n",
"axs[1,0].plot(T,[s[12] for s in sols[2]],label=r\"n: $(\\alpha_f,\\alpha_g) = $ (-8.71,-6.66) pcm\")\n",
"axs[1,0].plot(T,[s[12] for s in sols[3]],label=r\"n: $(\\alpha_f,\\alpha_g) = $ (-8.71,-6.66) pcm\")\n",
"axs[1,0].plot(T,[s[12] for s in sols[4]],label=r\"n: $(\\alpha_f,\\alpha_g) = $ (-8.71,-6.66) pcm\")\n",
"axs[1,0].set_xlabel(\"t (s)\")\n",
"axs[1,0].set_title(r\"$n$\")\n",
"axs[1,0].set_ylabel(r\"$\\frac{n}{n_0}$\")\n",
"\n",
"# reactivity\n",
"#axs[1,0].set_ylim([0,0])\n",
"#axs[1,0].set_xlim([0,500])\n",
"axs[1,1].set_xlim([0,t_stop])\n",
"axs[1,1].plot(T,[s[22] for s in sols[0]],label=r\"$\\rho$: $(\\alpha_f,\\alpha_g) = $ (-8.71,-6.66) pcm\")\n",
"axs[1,1].plot(T,[s[22] for s in sols[1]],label=r\"$\\rho$: $(\\alpha_f,\\alpha_g) = $ (-8.71,-6.66) pcm\")\n",
"axs[1,1].plot(T,[s[22] for s in sols[2]],label=r\"$\\rho$: $(\\alpha_f,\\alpha_g) = $ (-8.71,-6.66) pcm\")\n",
"axs[1,1].plot(T,[s[22] for s in sols[3]],label=r\"$\\rho$: $(\\alpha_f,\\alpha_g) = $ (-8.71,-6.66) pcm\")\n",
"axs[1,1].plot(T,[s[22] for s in sols[4]],label=r\"$\\rho$: $(\\alpha_f,\\alpha_g) = $ (-8.71,-6.66) pcm\")\n",
"axs[1,1].set_xlabel(\"t (s)\")\n",
"axs[1,1].set_title(r\"$\\rho$\")\n",
"\n",
"#plt.figure(figsize=(8,6))\n",
"#plt.plot(T,[s[12] for s in sols[0]],label=r\"$(\\alpha_f,\\alpha_g) = $ (-8.71,-6.66) pcm\")\n",
"#plt.plot(T,[s[12] for s in sols[1]],label=r\"$(\\alpha_f,\\alpha_g) = $ (-4.36,-3.33) pcm\")\n",
"#plt.plot(T,[s[12] for s in sols[2]],label=r\"$(\\alpha_f,\\alpha_g) = $ (-2.18,-1.67) pcm\")\n",
"#plt.plot(T,[s[12] for s in sols[3]],label=r\"$(\\alpha_f,\\alpha_g) = $ (-1.09,-0.83) pcm\")\n",
"#plt.plot(T,[s[12] for s in sols[4]],label=r\"$(\\alpha_f,\\alpha_g) = $ (0.0,0.0) pcm\")\n",
"#plt.xlabel(\"t (s)\")\n",
"#plt.ylabel(r\"$n/n_0$\",rotation=0)\n",
"#plt.legend()\n",
"#plt.ylim([0,20])\n",
"#plt.xlim([0,500])\n",
"#plt.title(r\"Reactor Response vs $\\alpha$\")"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Generating, compiling, and loading C code.\n",
"Using default integration parameters.\n",
"Generating, compiling, and loading C code.\n",
"Using default integration parameters.\n",
"Generating, compiling, and loading C code.\n",
"Using default integration parameters.\n",
"Generating, compiling, and loading C code.\n",
"Using default integration parameters.\n",
"Generating, compiling, and loading C code.\n",
"Using default integration parameters.\n"
]
}
],
"source": [
"# redefine system and get new solution set \n",
"\n",
"# delays\n",
"taus = [tau_c, tau_c_hx, tau_hx_c, tau_hx_r, tau_l, tau_r_hx]\n",
"\n",
"# adjust delay term \n",
"#tau_hx_c *= 0.5\n",
"tau_c_hx *= 0.5\n",
"#tau_r_hx *= 0.5\n",
"\n",
"# dT/dt for radiator nodes\n",
"T_out_rc = (W_rp/mn_rp)*(y(11,t-tau_hx_r)-y(0)) + (hA_rpn/mcp_rpn)*(y(1)-y(0)) # T_out_rc: y(0)\n",
"T_out_air = -((W_rs/mn_rs)+(hA_rsn/mcp_rsn))*y(1) + (hA_rsn/mcp_rsn)*y(0) + (W_rs/mn_rs)*Trs_in # T_out_air: y(1)\n",
"\n",
"# dT/dt for heat exchanger nodes\n",
"T_hf1 = -((W_p/mn_p)+(hA_pn/mcp_pn))*y(2) + (hA_pn/mcp_pn)*y(6) + (W_p/mn_p)*y(21,t-tau_c_hx) # T_hf1: y(2)\n",
"T_hf2 = (W_p/mn_p)*(y(2)-y(3)) + (hA_pn/mcp_pn)*(y(6)-y(2)) # T_hf2: y(3)\n",
"T_hf3 = -((W_p/mn_p)+(hA_pn/mcp_pn))*y(4) + (hA_pn/mcp_pn)*y(7) + (W_p/mn_p)*y(3) # T_hf3: y(4)\n",
"T_hf4 = (W_p/mn_p)*(y(4)-y(5)) + (hA_pn/mcp_pn)*(y(7)-y(4)) # T_hf4: y(5)\n",
"T_ht1 = (2*hA_pn/mcp_tn)*(y(2)-y(6)) + (2*hA_sn/mcp_tn)*(y(10)-y(6)) # T_ht1: y(6)\n",
"T_ht2 = (2*hA_pn/mcp_tn)*(y(4)-y(7)) + (2*hA_sn/mcp_tn)*(y(8)-y(7)) # T_ht2: y(7)\n",
"T_hc1 = -((W_s/mn_s)+(hA_sn/mcp_sn))*y(8) + (hA_sn/mcp_sn)*y(7) + (W_s/mn_s)*y(0,t-tau_r_hx) # T_hc1: y(8)\n",
"T_hc2 = (W_s/mn_s)*(y(8)-y(9)) + (hA_sn/mcp_sn)*(y(7)-y(8)) # T_hc2: y(9)\n",
"T_hc3 = -((W_s/mn_s)+(hA_sn/mcp_sn))*y(10) + (hA_sn/mcp_sn)*y(6) + (W_s/mn_s)*y(9) # T_hc3: y(10)\n",
"T_hc4 = (W_s/mn_s)*(y(10)-y(11)) + (hA_sn/mcp_sn)*(y(6)-y(10)) # T_hc4: y(11)\n",
"\n",
"# dn/dt\n",
"n = (y(22)-beta_t)*y(12)/Lam+lam[0]*y(13)+lam[1]*y(14)+lam[2]*y(15)+lam[3]*y(16)+lam[4]*y(17)+lam[5]*y(18) # n (no source insertion): y(12)\n",
"\n",
"# dC_i/dt (precursor concentrations)\n",
"C1 = y(12)*beta[0]/Lam-lam[0]*y(13)-y(13)/tau_c+y(13,t-tau_l)*np.exp(-lam[0]*tau_l)/tau_c # C1: y(13)\n",
"C2 = y(12)*beta[1]/Lam-lam[1]*y(14)-y(14)/tau_c+y(14,t-tau_l)*np.exp(-lam[1]*tau_l)/tau_c # C2: y(14)\n",
"C3 = y(12)*beta[2]/Lam-lam[2]*y(15)-y(15)/tau_c+y(15,t-tau_l)*np.exp(-lam[2]*tau_l)/tau_c # C3: y(15)\n",
"C4 = y(12)*beta[3]/Lam-lam[3]*y(16)-y(16)/tau_c+y(16,t-tau_l)*np.exp(-lam[3]*tau_l)/tau_c # C4: y(16)\n",
"C5 = y(12)*beta[4]/Lam-lam[4]*y(17)-y(17)/tau_c+y(17,t-tau_l)*np.exp(-lam[4]*tau_l)/tau_c # C5: y(17)\n",
"C6 = y(12)*beta[5]/Lam-lam[5]*y(18)-y(18)/tau_c+y(18,t-tau_l)*np.exp(-lam[5]*tau_l)/tau_c # C6: y(18)\n",
"\n",
"# dT/dt core nodes\n",
"T_cg = (hA_fg/mcp_g1)*(y(20)-y(19)) + k_g*P*y(12)/mcp_g1 # T_cg: y(19)\n",
"T_cf1 = W_f/mn_f*(y(5,t-tau_hx_c)-y(20)) + (k_f1*P*y(12)/mcp_f1) + (hA_fg*k_1*(y(19)-y(20))/mcp_f1) # T_cf1: y(20)\n",
"T_cf2 = W_f/mn_f*(y(20)-y(21)) + (k_f2*P*y(12)/mcp_f2) + (hA_fg*k_2*(y(19)-y(20))/mcp_f2) # T_cf2: y(21)\n",
"\n",
"# rho y(22)\n",
"rho = (a_f/2)*(T_cf1 + T_cf2) + (a_g)*(T_cg)\n",
"\n",
"# initial reactivity \n",
"rho_initial = 0.000\n",
"\n",
"# instantiate jitcdde object\n",
"DDE = jitcdde([T_out_rc,T_out_air,T_hf1,T_hf2,T_hf3,T_hf4,T_ht1,T_ht2,T_hc1,\n",
" T_hc2,T_hc3,T_hc4,n,C1,C2,C3,C4,C5,C6,T_cg,T_cf1,T_cf2,rho])\n",
"\n",
"# set initial conditions\n",
"DDE.constant_past([T0_rp, T0_rs, T0_p1,T0_p2, T0_p3, T0_p4, T0_t1, T0_t2, T0_s1, T0_s2, \n",
" T0_s3, T0_s4, n_frac0, C0[0], C0[1], C0[2], C0[3], C0[4], C0[5], \n",
" T0_g1, T0_f1, T0_f2,rho_initial])\n",
"\n",
"# jitcdde solver parameters \n",
"t0 = 0.0\n",
"tf = 1000.00\n",
"T = np.arange(t0,tf,0.01)\n",
"\n",
"sol_jit = []\n",
"for t_x in T:\n",
" sol_jit.append(DDE.integrate(t_x))\n",
"\n",
"sols.append(sol_jit)\n",
"\n",
"# feedbacks\n",
"factors_f = [1/2, 1/4, 1/8, 0]\n",
"factors_g = [1/2, 1/4, 1/8, 0]\n",
"\n",
"for i in range(len(factors_f)):\n",
"\n",
" # reinstantiate jitcdde objects\n",
" import sys\n",
" sys.modules.pop('jitcdde')\n",
" from jitcdde import jitcdde, y, t\n",
"\n",
" # rho y(22)\n",
" a_f_i = factors_f[i]\n",
" a_g_i = factors_g[i]\n",
" rho = a_f_i*(a_f/2)*(T_cf1 + T_cf2) + a_g_i*(a_g)*(T_cg)\n",
"\n",
" # instantiate jitcdde object\n",
" DDE = jitcdde([T_out_rc,T_out_air,T_hf1,T_hf2,T_hf3,T_hf4,T_ht1,T_ht2,T_hc1,\n",
" T_hc2,T_hc3,T_hc4,n,C1,C2,C3,C4,C5,C6,T_cg,T_cf1,T_cf2,rho])\n",
"\n",
" # set initial conditions\n",
" DDE.constant_past([T0_rp, T0_rs, T0_p1,T0_p2, T0_p3, T0_p4, T0_t1, T0_t2, T0_s1, T0_s2, \n",
" T0_s3, T0_s4, n_frac0, C0[0], C0[1], C0[2], C0[3], C0[4], C0[5], \n",
" T0_g1, T0_f1, T0_f2,rho_initial])\n",
" \n",
" sol_i = []\n",
" for t_x in T:\n",
" sol_i.append(DDE.integrate(t_x))\n",
" sols.append(sol_i)\n"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.lines.Line2D at 0xffff3c777e50>"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAABcsAAAPxCAYAAAAlpfUnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOzdd3wWRf7A8c+Wp+VJTwhJILTQmyggIigoqNg71hPs5+l5nmc57zwVez/bWX939oqKZzkbooI0RUEQpfeEFNKTp+7u/P54kgdCQgkEgvJ989rX7s7uMzu7PHlmZ3Z2RlNKKYQQQgghhBBCCCGEEEKI/Zje1gkQQgghhBBCCCGEEEIIIdqaVJYLIYQQQgghhBBCCCGE2O9JZbkQQgghhBBCCCGEEEKI/Z5UlgshhBBCCCGEEEIIIYTY70lluRBCCCGEEEIIIYQQQoj9nlSWCyGEEEIIIYQQQgghhNjvSWW5EEIIIYQQQgghhBBCiP2eVJYLIYQQQgghhBBCCCGE2O9JZbkQQgghhBBCCCGEEEKI/Z5Ulov9kqZp3HbbbW2djCYmTpxIly5d2joZ4lfk22+/xe12s3bt2l36/CeffEJiYiKlpaWtnDIhhGhM8jhYs2YNmqbxwgsvtHVSxK/Icccdx6WXXrrLnz/kkEO44YYbWjFFQojfqq+++gpN03j77bfbOimtTu5DREtJWXv/JZXlok298MILaJrW7PTXv/61TdPWUKDVNI133nmnyfbbbrsNTdPYtGlTG6Ru+7ZM+46mNWvWtHVy28zdd9/Ne++919bJ2C1///vfOeecc+jcuXOTbVOmTOHYY48lMzMTt9tNbm4u48ePZ9q0afF9xo0bR/fu3bnnnnv2ZrKFEG1o5cqVXH755XTr1g2v10tycjIjRozg0UcfJRgMtnXydsusWbO47bbbqKys3Kn9J06ciKZpDBw4EKVUk+2apnHVVVe1cipbR0PadzRNnDixrZPaZlr6fdgXzZw5k88++4wbb7yxybbi4mKuu+46evfuTUJCAn6/n8GDB3PnnXc2Oucbb7yRf/3rXxQVFe3FlAsh9padLfd99dVXbZ3UOClr//ZJWVvK2r9mZlsnQAiA22+/na5duzYK69+/fxulpqnbb7+d0047DU3T2jopO6Vdu3a8/PLLjcIeeughNmzYwD//+c8m++6v7r77bs444wxOOeWUtk7KLlmwYAFTp05l1qxZjcKVUlx00UW88MILHHjggVx77bVkZ2ezceNGpkyZwpgxY5g5cyaHHnooAJdffjnXXXcdkyZNIikpqS1ORQixl3z00UeceeaZeDweLrjgAvr3708kEuGbb77h+uuvZ/HixTz77LNtncxdNmvWLCZNmsTEiRNJTU3d6c8tWrSId999l9NPP33PJa6VXX755YwdOza+vnr1am655RYuu+wyDjvssHh4fn5+WyRvn7Cr34d9yQMPPMCYMWPo3r17o/DvvvuO4447jtraWs4//3wGDx4MwLx587j33nuZPn06n332GQAnn3wyycnJPPnkk9x+++17/RyEEHvW1uW+l156ic8//7xJeJ8+ffZmsnaalLV/m6SsLWXtXzOpLBf7hGOPPZYhQ4a0dTKaNWjQIBYsWMCUKVM47bTT2jo5O8Xv93P++ec3CnvjjTeoqKhoEv5b4TgOkUgEr9e736Tj+eefp1OnThxyyCGNwh966CFeeOEFrrnmGh5++OFGN55///vfefnllzHNzT//p59+On/84x+ZPHkyF1100R5PtxCibaxevZqzzz6bzp07M23aNHJycuLbrrzySlasWMFHH33UhilsGz6fj7y8vF9dYX348OEMHz48vj5v3jxuueUWhg8f/pvN6+vq6vD7/W2djL2WjpKSEj766COefvrpRuGVlZWceuqpGIbB/Pnz6d27d6Ptd911F88991x8Xdd1zjjjDF566SUmTZr0q/mOCyF2zta/+XPmzOHzzz//VeQFUtb+dZCy9mZS1t4/SDcsYp+3rf7Fu3Tp0uTV4srKSq655hry8vLweDx0796d++67D8dxdvn4Z599Nj179uT2229v9hXtrU2ePJnBgwfj8/nIzMzk/PPPp6CgoMl+7733Hv3798fr9dK/f3+mTJnSbHyO4/DII4/Qr18/vF4v7du35/LLL6eiomKXz6lBOBzm1ltvpXv37ng8HvLy8rjhhhsIh8ON9mt4DX3y5Mn07dsXn8/H8OHDWbRoEQDPPPMM3bt3x+v1Mnr06Cavm40ePZr+/fvz/fffc+ihh+Lz+ejatWuTwt+upOnVV1+lX79+eDwePvnkEwAefPBBDj30UDIyMvD5fAwePLhJv3uaplFXV8eLL77Y5FX1bfVn1/A64M6mo6CggIsuuoj27dvj8Xjo168f//nPf5rE+/jjj9OvXz8SEhJIS0tjyJAhvPbaa03229p7773HkUce2ShNwWCQe+65h969e/Pggw82WyD+3e9+x8EHHxxfz8rKYuDAgfz3v//d4TGFEL9e999/P7W1tfz73/9uVFHeoHv37vzpT3+Kr1uWxR133EF+fj4ej4cuXbrwt7/9rcnvMcCTTz4Z/w3Mzc3lyiuv3KmuL3bm9xo2/9Y25J0Nv6kNv7cQ+42+/vrrAejatetOvwKt6zo333wzCxcu3GZevKWSkhIuvvhi2rdvj9fr5YADDuDFF19ssl9lZSUTJ04kJSWF1NRUJkyYsM1rsmTJEs444wzS09Pxer0MGTKE999/f4dp2Rlz585l3LhxpKSkkJCQwKhRo5g5c2ajfRryt2XLlnH++eeTkpJCu3bt+Mc//oFSivXr18dbJ2dnZ/PQQw81+nxDH7dvvvkmf/vb38jOzsbv93PSSSexfv363UrTzz//zLnnnktaWhojR44EYOHChUycODHelVB2djYXXXQRZWVljT6/re/D9vqO3/q+c3vpAHjllVfi933p6emcffbZTc55+fLlnH766WRnZ+P1eunYsSNnn302VVVV2/mfi70JYllWozcIIHbfVVBQwMMPP9ykohygffv23HzzzY3CjjrqKNauXcuCBQu2e0whxG9Tc2VniJXTRo8e3STctu1W+z3fHilrS1m7gZS1xb5EWpaLfUJVVVWT/sgyMzNbFEcgEGDUqFEUFBRw+eWX06lTJ2bNmsVNN93Exo0beeSRR3YpbYZhcPPNN3PBBRfs8In3Cy+8wIUXXsjQoUO55557KC4u5tFHH2XmzJnMnz8//grwZ599xumnn07fvn255557KCsr48ILL6Rjx45N4rz88svj8V599dWsXr2aJ554gvnz5zNz5kxcLtcunZfjOJx00kl88803XHbZZfTp04dFixbxz3/+k2XLljXpX2zGjBm8//77XHnllQDcc889nHDCCdxwww08+eST/OEPf6CiooL777+fiy66qFFfXQAVFRUcd9xxjB8/nnPOOYe33nqLK664ArfbHX/C2tI0TZs2jbfeeourrrqKzMzMeKb76KOPctJJJ3HeeecRiUR44403OPPMM/nwww85/vjjgdjripdccgkHH3wwl112GbDrr6o3l47i4mIOOeSQeAbfrl07Pv74Yy6++GKqq6u55pprAHjuuee4+uqrOeOMM/jTn/5EKBRi4cKFzJ07l3PPPXebxywoKGDdunUcdNBBjcK/+eYbysvLueaaazAMY6fPYfDgwb/6PuWEENv3wQcf0K1bt/hroTtyySWX8OKLL3LGGWfwl7/8hblz53LPPffwyy+/NCp03nbbbUyaNImxY8dyxRVXsHTpUp566im+++67HeZTO/N73eCbb77h3Xff5Q9/+ANJSUk89thjnH766axbt46MjAxOO+00li1bxuuvv84///nP+H3EzrwCfe6553LHHXdw++23c+qpp26z5W0wGGT06NGsWLGCq666iq5duzJ58mQmTpxIZWVl/GGDUoqTTz6Zb775ht///vf06dOHKVOmMGHChCZxLl68mBEjRtChQwf++te/4vf7eeuttzjllFN45513OPXUU3eY/m2ZNm0axx57LIMHD+bWW29F13Wef/55jjzySGbMmNG
"text/plain": [
"<Figure size 1800x1200 with 6 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"ename": "",
"evalue": "",
"output_type": "error",
"traceback": [
"\u001b[1;31mThe Kernel crashed while executing code in the the current cell or a previous cell. Please review the code in the cell(s) to identify a possible cause of the failure. Click <a href='https://aka.ms/vscodeJupyterKernelCrash'>here</a> for more info. View Jupyter <a href='command:jupyter.viewOutput'>log</a> for further details."
]
}
],
"source": [
"fig,axs = plt.subplots(2,3,figsize=(18,12))\n",
"\n",
"t_stop = 30\n",
"vert = 10.555\n",
"\n",
"# fuel temps\n",
"#axs[0,0].set_xlim([0, 500])\n",
"#axs[0,0].set_ylim([600, 1000])\n",
"axs[0,0].set_xlim([0,t_stop])\n",
"axs[0,0].plot(T,[s[20] for s in sols[0]],label=r\"core 1: $(\\alpha_f,\\alpha_g) = $ (-8.71,-6.66) pcm\",color='b')\n",
"axs[0,0].plot(T,[s[20] for s in sols[1]],label=r\"core 1: $(\\alpha_f,\\alpha_g) = $ (-4.36,-3.33) pcm\",color='g')\n",
"axs[0,0].plot(T,[s[20] for s in sols[2]],label=r\"core 1: $(\\alpha_f,\\alpha_g) = $ (-2.18,-1.67) pcm\",color='r')\n",
"axs[0,0].plot(T,[s[20] for s in sols[3]],label=r\"core 1: $(\\alpha_f,\\alpha_g) = $ (-1.09,-0.83) pcm\",color='y')\n",
"axs[0,0].plot(T,[s[20] for s in sols[4]],label=r\"core 1: $(\\alpha_f,\\alpha_g) = $ (-0.00,-0.00) pcm\",color='m')\n",
"axs[0,0].legend()\n",
"axs[0,0].set_title(\"Fuel Node Temperatures (C)\")\n",
"axs[0,0].tick_params(\n",
" axis='x', # changes apply to the x-axis\n",
" which='both', # both major and minor ticks are affected\n",
" bottom=False, # ticks along the bottom edge are off\n",
" top=False, # ticks along the top edge are off\n",
" labelbottom=False) # labels along the bottom edge are off\n",
"\n",
"# coolant temps\n",
"#axs[0,1].set_xlim([0, 1000])\n",
"#axs[0,1].set_ylim([500, 1000])\n",
"axs[0,1].set_xlim([0,t_stop])\n",
"axs[0,1].plot(T,[s[8] for s in sols[0]],label=r\"hx 1: $(\\alpha_f,\\alpha_g) = $ (-8.71,-6.66) pcm\",color='b')\n",
"axs[0,1].plot(T,[s[8] for s in sols[1]],label=r\"hx 1: $(\\alpha_f,\\alpha_g) = $ (-4.36,-3.33) pcm\",color='g')\n",
"axs[0,1].plot(T,[s[8] for s in sols[2]],label=r\"hx 1: $(\\alpha_f,\\alpha_g) = $ (-2.18,-1.67) pcm\",color='r')\n",
"axs[0,1].plot(T,[s[8] for s in sols[3]],label=r\"hx 1: $(\\alpha_f,\\alpha_g) = $ (-1.09,-0.83) pcm\",color='y')\n",
"axs[0,1].plot(T,[s[8] for s in sols[4]],label=r\"hx 1: $(\\alpha_f,\\alpha_g) = $ (-0.00,-0.00) pcm\",color='m')\n",
"axs[0,1].legend()\n",
"axs[0,1].set_title(\"Coolant Node Temperatures (C)\")\n",
"axs[0,1].tick_params(\n",
" axis='x', # changes apply to the x-axis\n",
" which='both', # both major and minor ticks are affected\n",
" bottom=False, # ticks along the bottom edge are off\n",
" top=False, # ticks along the top edge are off\n",
" labelbottom=False) # labels along the bottom edge are off\n",
"\n",
"# tube node temps\n",
"#axs[0,2].set_xlim([0, 500])\n",
"#axs[0,2].set_ylim([500, 1000])\n",
"axs[0,2].set_xlim([0,t_stop])\n",
"axs[0,2].plot(T,[s[8] for s in sols[0]],label=r\"t1: $(\\alpha_f,\\alpha_g) = $ (-8.71,-6.66) pcm\",color='b')\n",
"axs[0,2].plot(T,[s[8] for s in sols[1]],label=r\"t1: $(\\alpha_f,\\alpha_g) = $ (-4.36,-3.33) pcm\",color='g')\n",
"axs[0,2].plot(T,[s[8] for s in sols[2]],label=r\"t1: $(\\alpha_f,\\alpha_g) = $ (-2.18,-1.67) pcm\",color='r')\n",
"axs[0,2].plot(T,[s[8] for s in sols[3]],label=r\"t1: $(\\alpha_f,\\alpha_g) = $ (-1.09,-0.83) pcm\",color='y')\n",
"axs[0,2].plot(T,[s[8] for s in sols[4]],label=r\"t1: $(\\alpha_f,\\alpha_g) = $ (-0.00,-0.00) pcm\",color='m')\n",
"\n",
"# tube node temps\n",
"axs[0,2].plot(T,[s[9] for s in sols[0]],label=r\"t2: $(\\alpha_f,\\alpha_g) = $ (-8.71,-6.66) pcm\",color='b')\n",
"axs[0,2].plot(T,[s[9] for s in sols[1]],label=r\"t2: $(\\alpha_f,\\alpha_g) = $ (-4.36,-3.33) pcm\",color='g')\n",
"axs[0,2].plot(T,[s[9] for s in sols[2]],label=r\"t2: $(\\alpha_f,\\alpha_g) = $ (-2.18,-1.67) pcm\",color='r')\n",
"axs[0,2].plot(T,[s[9] for s in sols[3]],label=r\"t2: $(\\alpha_f,\\alpha_g) = $ (-1.09,-0.83) pcm\",color='y')\n",
"axs[0,2].plot(T,[s[9] for s in sols[4]],label=r\"t2: $(\\alpha_f,\\alpha_g) = $ (-0.00,-0.00) pcm\",color='m')\n",
"\n",
"axs[0,2].legend()\n",
"axs[0,2].set_title(\"Tube Node Temperatures (C)\")\n",
"axs[0,2].tick_params(\n",
" axis='x', # changes apply to the x-axis\n",
" which='both', # both major and minor ticks are affected\n",
" bottom=False, # ticks along the bottom edge are off\n",
" top=False, # ticks along the top edge are off\n",
" labelbottom=False) # labels along the bottom edge are off\n",
"\n",
"# precursor concentrations\n",
"#axs[1,2].set_xlim([0, 500])\n",
"#axs[1,2].set_ylim([1e1,1e3])\n",
"axs[1,2].set_xlim([0,t_stop])\n",
"axs[1,2].plot(T,[s[13] for s in sols[0]],label=r\"C1: $(\\alpha_f,\\alpha_g) = $ (-8.71,-6.66) pcm\",color='b')\n",
"axs[1,2].plot(T,[s[13] for s in sols[1]],label=r\"C1: $(\\alpha_f,\\alpha_g) = $ (-4.36,-3.33) pcm\",color='g')\n",
"axs[1,2].plot(T,[s[13] for s in sols[2]],label=r\"C1: $(\\alpha_f,\\alpha_g) = $ (-2.18,-1.67) pcm\",color='r')\n",
"axs[1,2].plot(T,[s[13] for s in sols[3]],label=r\"C1: $(\\alpha_f,\\alpha_g) = $ (-1.09,-0.83) pcm\",color='y')\n",
"axs[1,2].plot(T,[s[13] for s in sols[4]],label=r\"C1: $(\\alpha_f,\\alpha_g) = $ (-0.00,-0.00) pcm\",color='m')\n",
"axs[1,2].legend()\n",
"axs[1,2].set_xlabel(\"t (s)\")\n",
"axs[1,2].set_yscale(\"log\")\n",
"axs[1,2].legend(loc=\"right\")\n",
"axs[1,2].set_title(\"Precursor Concentrations\")\n",
"axs[1,2].set_ylabel(r\"concentration (1/cm$^3$)\")\n",
"\n",
"# multiplication factor \n",
"#axs[1,0].set_ylim([0,20])\n",
"#axs[1,0].set_xlim([0,500])\n",
"axs[1,0].set_xlim([0,t_stop])\n",
"axs[1,0].plot(T,[s[12] for s in sols[0]],label=r\"n: $(\\alpha_f,\\alpha_g) = $ (-8.71,-6.66) pcm\",color='b')\n",
"axs[1,0].plot(T,[s[12] for s in sols[1]],label=r\"n: $(\\alpha_f,\\alpha_g) = $ (-8.71,-6.66) pcm\",color='g')\n",
"axs[1,0].plot(T,[s[12] for s in sols[2]],label=r\"n: $(\\alpha_f,\\alpha_g) = $ (-8.71,-6.66) pcm\",color='r')\n",
"axs[1,0].plot(T,[s[12] for s in sols[3]],label=r\"n: $(\\alpha_f,\\alpha_g) = $ (-8.71,-6.66) pcm\",color='y')\n",
"axs[1,0].plot(T,[s[12] for s in sols[4]],label=r\"n: $(\\alpha_f,\\alpha_g) = $ (-8.71,-6.66) pcm\",color='m')\n",
"axs[1,0].set_xlabel(\"t (s)\")\n",
"axs[1,0].set_title(r\"$n$\")\n",
"axs[1,0].set_ylabel(r\"$\\frac{n}{n_0}$\")\n",
"axs[1,0].legend()\n",
"\n",
"# reactivity\n",
"#axs[1,0].set_ylim([0,0])\n",
"#axs[1,0].set_xlim([0,500])\n",
"axs[1,1].set_xlim([0,t_stop])\n",
"axs[1,1].plot(T,[s[22] for s in sols[0]],label=r\"$\\rho$: $(\\alpha_f,\\alpha_g) = $ (-8.71,-6.66) pcm\",color='b')\n",
"axs[1,1].plot(T,[s[22] for s in sols[1]],label=r\"$\\rho$: $(\\alpha_f,\\alpha_g) = $ (-8.71,-6.66) pcm\",color='g')\n",
"axs[1,1].plot(T,[s[22] for s in sols[2]],label=r\"$\\rho$: $(\\alpha_f,\\alpha_g) = $ (-8.71,-6.66) pcm\",color='r')\n",
"axs[1,1].plot(T,[s[22] for s in sols[3]],label=r\"$\\rho$: $(\\alpha_f,\\alpha_g) = $ (-8.71,-6.66) pcm\",color='y')\n",
"axs[1,1].plot(T,[s[22] for s in sols[4]],label=r\"$\\rho$: $(\\alpha_f,\\alpha_g) = $ (-8.71,-6.66) pcm\",color='m')\n",
"axs[1,1].set_xlabel(\"t (s)\")\n",
"axs[1,1].set_title(r\"$\\rho$\")\n",
"axs[1,1].legend()\n",
"\n",
"#plt.figure(figsize=(8,6))\n",
"#plt.plot(T,[s[12] for s in sols[0]],label=r\"$(\\alpha_f,\\alpha_g) = $ (-8.71,-6.66) pcm\")\n",
"#plt.plot(T,[s[12] for s in sols[1]],label=r\"$(\\alpha_f,\\alpha_g) = $ (-4.36,-3.33) pcm\")\n",
"#plt.plot(T,[s[12] for s in sols[2]],label=r\"$(\\alpha_f,\\alpha_g) = $ (-2.18,-1.67) pcm\")\n",
"#plt.plot(T,[s[12] for s in sols[3]],label=r\"$(\\alpha_f,\\alpha_g) = $ (-1.09,-0.83) pcm\")\n",
"#plt.plot(T,[s[12] for s in sols[4]],label=r\"$(\\alpha_f,\\alpha_g) = $ (0.0,0.0) pcm\")\n",
"#plt.xlabel(\"t (s)\")\n",
"#plt.ylabel(r\"$n/n_0$\",rotation=0)\n",
"#plt.legend()\n",
"#plt.ylim([0,20])\n",
"#plt.xlim([0,500])\n",
"#plt.title(r\"Reactor Response vs $\\alpha$\")\n",
"\n",
"\n",
"#######################################################################################################\n",
"\n",
"# adjusted delay term \n",
"\n",
"# fuel temps\n",
"#axs[0,0].set_xlim([0, 500])\n",
"#axs[0,0].set_ylim([600, 1000])\n",
"axs[0,0].set_xlim([0,t_stop])\n",
"axs[0,0].plot(T,[s[20] for s in sols[5]],linestyle='-.',color='b')\n",
"axs[0,0].plot(T,[s[20] for s in sols[6]],linestyle='-.',color='g')\n",
"axs[0,0].plot(T,[s[20] for s in sols[7]],linestyle='-.',color='r')\n",
"axs[0,0].plot(T,[s[20] for s in sols[8]],linestyle='-.',color='y')\n",
"axs[0,0].plot(T,[s[20] for s in sols[9]],linestyle='-.',color='m')\n",
"axs[0,0].legend()\n",
"axs[0,0].set_title(\"Fuel Node Temperatures (C)\")\n",
"axs[0,0].tick_params(\n",
" axis='x', # changes apply to the x-axis\n",
" which='both', # both major and minor ticks are affected\n",
" bottom=False, # ticks along the bottom edge are off\n",
" top=False, # ticks along the top edge are off\n",
" labelbottom=False) # labels along the bottom edge are off\n",
"axs[0,0].axvline(x=vert, color='black', linestyle=':')\n",
"# coolant temps\n",
"#axs[0,1].set_xlim([0, 1000])\n",
"#axs[0,1].set_ylim([500, 1000])\n",
"axs[0,1].set_xlim([0,t_stop])\n",
"axs[0,1].plot(T,[s[8] for s in sols[5]],linestyle='-.',color='b')\n",
"axs[0,1].plot(T,[s[8] for s in sols[6]],linestyle='-.',color='g')\n",
"axs[0,1].plot(T,[s[8] for s in sols[7]],linestyle='-.',color='r')\n",
"axs[0,1].plot(T,[s[8] for s in sols[8]],linestyle='-.',color='y')\n",
"axs[0,1].plot(T,[s[8] for s in sols[9]],linestyle='-.',color='m')\n",
"axs[0,1].legend()\n",
"axs[0,1].set_title(\"Coolant Node Temperatures (C)\")\n",
"axs[0,1].tick_params(\n",
" axis='x', # changes apply to the x-axis\n",
" which='both', # both major and minor ticks are affected\n",
" bottom=False, # ticks along the bottom edge are off\n",
" top=False, # ticks along the top edge are off\n",
" labelbottom=False) # labels along the bottom edge are off\n",
"\n",
"# tube node temps\n",
"#axs[0,2].set_xlim([0, 500])\n",
"#axs[0,2].set_ylim([500, 1000])\n",
"axs[0,2].set_xlim([0,t_stop])\n",
"axs[0,2].plot(T,[s[8] for s in sols[5]],linestyle='-.',color='b')\n",
"axs[0,2].plot(T,[s[8] for s in sols[6]],linestyle='-.',color='g')\n",
"axs[0,2].plot(T,[s[8] for s in sols[7]],linestyle='-.',color='r')\n",
"axs[0,2].plot(T,[s[8] for s in sols[8]],linestyle='-.',color='y')\n",
"axs[0,2].plot(T,[s[8] for s in sols[9]],linestyle='-.',color='m')\n",
"\n",
"# tube node temps\n",
"axs[0,2].plot(T,[s[9] for s in sols[5]],linestyle='-.',color='b')\n",
"axs[0,2].plot(T,[s[9] for s in sols[6]],linestyle='-.',color='g')\n",
"axs[0,2].plot(T,[s[9] for s in sols[7]],linestyle='-.',color='r')\n",
"axs[0,2].plot(T,[s[9] for s in sols[8]],linestyle='-.',color='y')\n",
"axs[0,2].plot(T,[s[9] for s in sols[9]],linestyle='-.',color='m')\n",
"\n",
"axs[0,2].legend()\n",
"axs[0,2].set_title(\"Tube Node Temperatures (C)\")\n",
"axs[0,2].tick_params(\n",
" axis='x', # changes apply to the x-axis\n",
" which='both', # both major and minor ticks are affected\n",
" bottom=False, # ticks along the bottom edge are off\n",
" top=False, # ticks along the top edge are off\n",
" labelbottom=False) # labels along the bottom edge are off\n",
"\n",
"# precursor concentrations\n",
"#axs[1,2].set_xlim([0, 500])\n",
"#axs[1,2].set_ylim([1e1,1e3])\n",
"axs[1,2].set_xlim([0,t_stop])\n",
"axs[1,2].plot(T,[s[13] for s in sols[5]],linestyle='-.',color='b')\n",
"axs[1,2].plot(T,[s[13] for s in sols[6]],linestyle='-.',color='g')\n",
"axs[1,2].plot(T,[s[13] for s in sols[7]],linestyle='-.',color='r')\n",
"axs[1,2].plot(T,[s[13] for s in sols[8]],linestyle='-.',color='y')\n",
"axs[1,2].plot(T,[s[13] for s in sols[9]],linestyle='-.',color='m')\n",
"axs[1,2].legend()\n",
"axs[1,2].set_xlabel(\"t (s)\")\n",
"axs[1,2].set_yscale(\"log\")\n",
"axs[1,2].legend(loc=\"right\")\n",
"axs[1,2].set_title(\"Precursor Concentrations\")\n",
"axs[1,2].set_ylabel(r\"concentration (1/cm$^3$)\")\n",
"\n",
"# multiplication factor \n",
"#axs[1,0].set_ylim([0,20])\n",
"#axs[1,0].set_xlim([0,500])\n",
"axs[1,0].set_xlim([0,t_stop])\n",
"axs[1,0].plot(T,[s[12] for s in sols[5]],linestyle='-.',color='b')\n",
"axs[1,0].plot(T,[s[12] for s in sols[6]],linestyle='-.',color='g')\n",
"axs[1,0].plot(T,[s[12] for s in sols[7]],linestyle='-.',color='r')\n",
"axs[1,0].plot(T,[s[12] for s in sols[8]],linestyle='-.',color='y')\n",
"axs[1,0].plot(T,[s[12] for s in sols[9]],linestyle='-.',color='m')\n",
"axs[1,0].set_xlabel(\"t (s)\")\n",
"axs[1,0].set_title(r\"$n$\")\n",
"axs[1,0].set_ylabel(r\"$\\frac{n}{n_0}$\")\n",
"axs[1,0].axvline(x=vert, color='black', linestyle=':')\n",
"\n",
"# reactivity\n",
"#axs[1,0].set_ylim([0,0])\n",
"#axs[1,0].set_xlim([0,500])\n",
"axs[1,1].set_xlim([0,t_stop])\n",
"axs[1,1].plot(T,[s[22] for s in sols[5]],linestyle='-.',color='b')\n",
"axs[1,1].plot(T,[s[22] for s in sols[6]],linestyle='-.',color='g')\n",
"axs[1,1].plot(T,[s[22] for s in sols[7]],linestyle='-.',color='r')\n",
"axs[1,1].plot(T,[s[22] for s in sols[8]],linestyle='-.',color='y')\n",
"axs[1,1].plot(T,[s[22] for s in sols[9]],linestyle='-.',color='m')\n",
"axs[1,1].set_xlabel(\"t (s)\")\n",
"axs[1,1].set_title(r\"$\\rho$\")\n",
"axs[1,1].axvline(x=vert, color='black', linestyle=':')\n",
"\n",
"#plt.figure(figsize=(8,6))\n",
"#plt.plot(T,[s[12] for s in sols[0]],label=r\"$(\\alpha_f,\\alpha_g) = $ (-8.71,-6.66) pcm\")\n",
"#plt.plot(T,[s[12] for s in sols[1]],label=r\"$(\\alpha_f,\\alpha_g) = $ (-4.36,-3.33) pcm\")\n",
"#plt.plot(T,[s[12] for s in sols[2]],label=r\"$(\\alpha_f,\\alpha_g) = $ (-2.18,-1.67) pcm\")\n",
"#plt.plot(T,[s[12] for s in sols[3]],label=r\"$(\\alpha_f,\\alpha_g) = $ (-1.09,-0.83) pcm\")\n",
"#plt.plot(T,[s[12] for s in sols[4]],label=r\"$(\\alpha_f,\\alpha_g) = $ (0.0,0.0) pcm\")\n",
"#plt.xlabel(\"t (s)\")\n",
"#plt.ylabel(r\"$n/n_0$\",rotation=0)\n",
"#plt.legend()\n",
"#plt.ylim([0,20])\n",
"#plt.xlim([0,500])\n",
"#plt.title(r\"Reactor Response vs $\\alpha$\")"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"10"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"len(sols)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "thesis_env",
"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.9.2"
},
"orig_nbformat": 4
},
"nbformat": 4,
"nbformat_minor": 2
}
Reference in a new issue Copy permalink