{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "9e47daa3",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import cantera as ct\n",
    "import pde\n",
    "\n",
    "import ipywidgets as widgets\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "%matplotlib widget"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "c49236bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import optimize as opt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b1f191f",
   "metadata": {},
   "source": [
    "# 이론공연비 계산"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "60f9c29b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def stoich_AF (fuel_compostion, mole_fraction=True):\n",
    "    \"Stoichiometric Air/Fuel ratio (volume)\"\n",
    "    \n",
    "    gas = ct.Solution('gri30.xml')\n",
    "\n",
    "    if mole_fraction:\n",
    "        gas.TPX = 600 + 273.15, ct.one_atm, fuel_compostion\n",
    "    else:\n",
    "        gas.TPY = 600 + 273.15, ct.one_atm, fuel_compostion\n",
    "    \n",
    "    fuel_comp = gas.X\n",
    "\n",
    "    # 연료의 원자 1 몰 완전 연소에 필요한 산소 원자 몰 수\n",
    "    coef_O = gas.elemental_mole_fraction('C') * 2 + gas.elemental_mole_fraction('H') / 2 - (gas.elemental_mole_fraction('O'))\n",
    "\n",
    "    # 연료의 원자 1 몰 완전 연소에 필요한 공기 원자 몰 수 = 산소 원자 몰 수 X (1 + 3.762)\n",
    "    coef_N = coef_O * 3.762\n",
    "\n",
    "    # 공기 1몰의 원자 몰수 = 2 ; O2 + N2 만으로 구성되었으면\n",
    "    mean_natoms_air = 2\n",
    "\n",
    "    # 연료 기체 1몰의 원자 몰수 (평균 분자당 원자 수)\n",
    "    mean_natoms_fuel = sum([(gas.X[i] * sum([v for v in gas.species(i).composition.values()])) for i in range(gas.n_species)])\n",
    "\n",
    "    # 연료 기체 1몰에 필요한 공기 몰수 ()\n",
    "    air_fuel_ratio_st = mean_natoms_fuel * (coef_N + coef_O) / mean_natoms_air\n",
    "\n",
    "    return (air_fuel_ratio_st)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "c633e758",
   "metadata": {},
   "outputs": [],
   "source": [
    "coke_oven_fuel = \"H2:6.42, O2:0.39, N2:47.28, CH4:1.79, CO:24.25, CO2:19.72, C2H4:0.13, C2H6:0.04\"\n",
    "\n",
    "air = ct.Solution('gri30.xml')\n",
    "air.X = \"O2:1, N2:3.762\"\n",
    "air_comp = air.X\n",
    "\n",
    "fuel = ct.Solution('gri30.xml')\n",
    "fuel.X = coke_oven_fuel\n",
    "fuel_comp = fuel.X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "a20ce147",
   "metadata": {},
   "outputs": [],
   "source": [
    "def exhaust_stoichiometry (phi = 1, fuel=fuel_comp, return_unburned=False):\n",
    "\n",
    "    mix = ct.Solution('gri30.xml')\n",
    "\n",
    "    # mix.TPX = 25+273.15, ct.one_atm, phi * fuel_comp + st_AF * air_comp\n",
    "    \n",
    "    mix.set_equivalence_ratio(phi, fuel=fuel_comp, oxidizer=\"O2:0.21, N2:0.79\")\n",
    "    \n",
    "    element_X = {ename: mix.elemental_mole_fraction(ename) for ename in mix.element_names}\n",
    "    \n",
    "    exhaust = ct.Solution('gri30.xml')\n",
    "    exhaust.X = {\n",
    "        \"CO2\" : element_X['C'],\n",
    "        \"H2O\" : element_X['H']/2,\n",
    "        \"O2\" : (element_X['O'] - 2*element_X['C'] - element_X['H']/2)/2,\n",
    "        \"N2\" : element_X['N']/2,\n",
    "    }\n",
    "    \n",
    "    if return_unburned:\n",
    "        return mix.mole_fraction_dict(threshold=-1), exhaust.mole_fraction_dict(threshold=-1)\n",
    "    else:\n",
    "        return exhaust.mole_fraction_dict(threshold=-1)\n",
    "\n",
    "def exhaust_equilibrium (phi = 1):\n",
    "\n",
    "    mix = ct.Solution('gri30.xml')\n",
    "\n",
    "    mix.TPX = 25+273.15, ct.one_atm, phi * fuel_comp + st_AF * air_comp\n",
    "    \n",
    "    mix.equilibrate(\"HP\")\n",
    "    \n",
    "    return mix.mole_fraction_dict(threshold=-1)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75feeff2",
   "metadata": {},
   "source": [
    "# A/F 계산기"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "bf8c7fe0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6a0200388bb241a4bc4c75fb070af441",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(Textarea(value='H2:6.42, O2:0.39, N2:47.28, CH4:1.79, CO:24.25, CO2:19.72, C2H4:0.13, C2…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "@widgets.interact(\n",
    "    composition=widgets.Textarea(\n",
    "    value='H2:6.42, O2:0.39, N2:47.28, CH4:1.79, CO:24.25, CO2:19.72, C2H4:0.13, C2H6:0.04',\n",
    "    placeholder='\"species_name_0:X0, species_name_1:X1, species_name_2:X2, ... \"',\n",
    "    description='''Fuel Composition:''',\n",
    "    disabled=False\n",
    "),\n",
    "    O2_in_exhaust=widgets.FloatSlider(\n",
    "    value=0,\n",
    "    min=0,\n",
    "    max=21.0,\n",
    "    step=0.1,\n",
    "    description='excess O2(%)',\n",
    "    readout=True,\n",
    "    readout_format='.1f',\n",
    ")\n",
    ")\n",
    "def calculate_AF (composition, O2_in_exhaust):\n",
    "    gas = ct.Solution('gri30.xml')\n",
    "    try:\n",
    "        gas.TPX = 300, ct.one_atm, composition\n",
    "    except ct.CanteraError:\n",
    "        gas.TPX = 300, ct.one_atm, \"H2:1\"\n",
    "        \n",
    "    print(\"Fuel Gas\")\n",
    "    gas()\n",
    "    \n",
    "    AF_st = stoich_AF(composition)\n",
    "    \n",
    "    print(f'stoichiometric A/F = {AF_st}')\n",
    "    \n",
    "    f_found = opt.root_scalar(lambda x: exhaust_stoichiometry(x, fuel=composition)[\"O2\"] - O2_in_exhaust/100, \n",
    "                          bracket=[1e-1, 1])\n",
    "    \n",
    "    # print(f_found)\n",
    "    \n",
    "    phi = f_found.root\n",
    "    print(f_found)\n",
    "    print(f'phi for {O2_in_exhaust}% O2 in exhaust gas = {phi}')\n",
    "    print(f'A/F(vol) for {O2_in_exhaust}% O2 in exhaust gas = {AF_st / phi}')\n",
    "    # mean molecular weight needed to calc mass A/F\n",
    "\n",
    "    gas.TPX = 300, ct.one_atm, exhaust_stoichiometry(phi)\n",
    "    gas()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c6c962f",
   "metadata": {},
   "source": [
    "# 발열량 계산"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "6c0407fe",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "      converged: True\n",
      "           flag: 'converged'\n",
      " function_calls: 9\n",
      "     iterations: 8\n",
      "           root: 0.6547411452355204\n"
     ]
    }
   ],
   "source": [
    "f_found = opt.root_scalar(lambda x: exhaust_stoichiometry(x)[\"O2\"] - 4.5/100, bracket=[1e-1, 1])\n",
    "print(f_found)\n",
    "phi_O2_045 = f_found.root"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "2f732468",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "Xu, Xb = exhaust_stoichiometry(phi=phi_O2_045, return_unburned=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7c5063d",
   "metadata": {},
   "source": [
    "## 배가스 산소 농도 4.5% 연료+공기 혼합 기체 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "53b2c1b9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "  gri30:\n",
      "\n",
      "       temperature   298.15 K\n",
      "          pressure   1.0133e+05 Pa\n",
      "           density   1.1869 kg/m^3\n",
      "  mean mol. weight   29.039 kg/kmol\n",
      "   phase of matter   gas\n",
      "\n",
      "                          1 kg             1 kmol     \n",
      "                     ---------------   ---------------\n",
      "          enthalpy       -1.5255e+06         -4.43e+07  J\n",
      "   internal energy       -1.6109e+06       -4.6779e+07  J\n",
      "           entropy            6999.5        2.0326e+05  J/K\n",
      "    Gibbs function       -3.6124e+06        -1.049e+08  J\n",
      " heat capacity c_p            1027.5             29837  J/K\n",
      " heat capacity c_v            741.17             21523  J/K\n",
      "\n",
      "                      mass frac. Y      mole frac. X     chem. pot. / RT\n",
      "                     ---------------   ---------------   ---------------\n",
      "                H2         0.0018679          0.026906           -19.333\n",
      "                O2            0.1362           0.12361           -26.764\n",
      "               CH4         0.0041445         0.0075019           -57.401\n",
      "                CO           0.09803           0.10163           -70.646\n",
      "               CO2           0.12525          0.082647           -186.95\n",
      "              C2H4        0.00052635        0.00054483           -12.715\n",
      "              C2H6        0.00017359        0.00016764           -70.088\n",
      "                N2            0.6338             0.657           -23.453\n",
      "     [  +45 minor]                 0                 0  \n",
      "\n"
     ]
    }
   ],
   "source": [
    "gas_u = ct.Solution('gri30.xml')\n",
    "gas_u.TPX = 25 + 273.15, ct.one_atm, Xu\n",
    "gas_u()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d54eebc",
   "metadata": {},
   "source": [
    "##  배가스 산소 농도 4.5% 완전 연소 배가스"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "a1bc3d15",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "  gri30:\n",
      "\n",
      "       temperature   298.15 K\n",
      "          pressure   1.0132e+05 Pa\n",
      "           density   1.2684 kg/m^3\n",
      "  mean mol. weight   31.031 kg/kmol\n",
      "   phase of matter   gas\n",
      "\n",
      "                          1 kg             1 kmol     \n",
      "                     ---------------   ---------------\n",
      "          enthalpy       -2.9804e+06       -9.2483e+07  J\n",
      "   internal energy       -3.0603e+06       -9.4962e+07  J\n",
      "           entropy            6565.1        2.0372e+05  J/K\n",
      "    Gibbs function       -4.9377e+06       -1.5322e+08  J\n",
      " heat capacity c_p            997.71             30960  J/K\n",
      " heat capacity c_v            729.77             22645  J/K\n",
      "\n",
      "                      mass frac. Y      mole frac. X     chem. pot. / RT\n",
      "                     ---------------   ---------------   ---------------\n",
      "                O2          0.046403             0.045           -27.775\n",
      "               H2O          0.026988          0.046486           -123.33\n",
      "               CO2           0.29281           0.20646           -186.03\n",
      "                N2            0.6338           0.70206           -23.387\n",
      "     [  +49 minor]                 0                 0  \n",
      "\n"
     ]
    }
   ],
   "source": [
    "gas_b = ct.Solution('gri30.xml')\n",
    "gas_b.TPX = 25 + 273.15, ct.one_atm, Xb\n",
    "gas_b()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e537f99",
   "metadata": {},
   "source": [
    "##  배가스 산소 농도 4.5% 혼합 기체 kg 당 저위발열량 (J/kg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "067ae187",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'1.454846654525596 MJ/kg'"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "heating_value_J_kg = gas_u.enthalpy_mass - gas_b.enthalpy_mass\n",
    "f'{heating_value_J_kg*1e-6} MJ/kg'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "f8713062",
   "metadata": {},
   "outputs": [],
   "source": [
    "heat_to_battery = 80 # GJ / Rev"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "42f68b81",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'3636.3636363636365 MJ / hr to a combustion chamber'"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "heat_to_chamber = heat_to_battery * 1000 * 3 /  66 # MJ / hr\n",
    "f'{heat_to_chamber} MJ / hr to a combustion chamber'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30be4bfc",
   "metadata": {},
   "source": [
    "## 연소실당 질량 유량 kg/hr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "fd85f6dc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'2499.482419712063 (kg/hr) = 3636.3636363636365 (MJ/hr) / 1.454846654525596 (MJ/kg)'"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mfr_to_chamber = heat_to_chamber * 1e6 / heating_value_J_kg # kg / hr\n",
    "f'{mfr_to_chamber} (kg/hr) = {heat_to_chamber} (MJ/hr) / {heating_value_J_kg  * 1e-6} (MJ/kg)'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed0ec84b",
   "metadata": {},
   "source": [
    "# 단열 화염 온도 (평형 계산)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "fa7ef7a2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "  gri30:\n",
      "\n",
      "       temperature   1990.3 K\n",
      "          pressure   1.0133e+05 Pa\n",
      "           density   0.18983 kg/m^3\n",
      "  mean mol. weight   31.002 kg/kmol\n",
      "   phase of matter   gas\n",
      "\n",
      "                          1 kg             1 kmol     \n",
      "                     ---------------   ---------------\n",
      "          enthalpy        -8.933e+05       -2.7694e+07  J\n",
      "   internal energy       -1.4271e+06       -4.4242e+07  J\n",
      "           entropy            8783.1         2.723e+05  J/K\n",
      "    Gibbs function       -1.8374e+07       -5.6964e+08  J\n",
      " heat capacity c_p            1347.2             41767  J/K\n",
      " heat capacity c_v              1079             33452  J/K\n",
      "\n",
      "                      mass frac. Y      mole frac. X     chem. pot. / RT\n",
      "                     ---------------   ---------------   ---------------\n",
      "                 O        6.6906e-05        0.00012965           -15.939\n",
      "                O2          0.045092          0.043688           -31.879\n",
      "                OH        0.00048958        0.00089246           -30.541\n",
      "               H2O          0.026691          0.045932           -45.143\n",
      "                CO         0.0010791         0.0011944           -41.097\n",
      "               CO2           0.29111           0.20507           -57.036\n",
      "                NO         0.0031141         0.0032175           -29.575\n",
      "                N2           0.63235            0.6998           -27.272\n",
      "     [  +45 minor]        8.9238e-06        7.3926e-05  \n",
      "\n"
     ]
    }
   ],
   "source": [
    "gas_chamber = ct.Solution('gri30.xml')\n",
    "gas_chamber.TPX = 600 + 273.15, ct.one_atm, Xu\n",
    "gas_chamber.equilibrate(\"HP\")\n",
    "gas_chamber(threshold=1e-4)\n",
    "\n",
    "eq_state = gas_chamber.TPX\n",
    "\n",
    "Tad, P0, X0 = eq_state"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "1d6c95a7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Tad = 1717.118741092177 `C\n"
     ]
    }
   ],
   "source": [
    "print(f'Tad = {eq_state[0] - 273.15} `C')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "e437f3b4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NO: 3217.4895678195617 ppm\n"
     ]
    }
   ],
   "source": [
    "print(f'NO: {gas_chamber.mole_fraction_dict()[\"NO\"] * 1e6} ppm')"
   ]
  }
 ],
 "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.9.12"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "384px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}