Compl.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "e2d3f6c8",
   "metadata": {},
   "outputs": [],
   "source": [
    "import openpyxl\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from tqdm import tqdm\n",
    "from gekko import GEKKO\n",
    "import matplotlib.pyplot as plt \n",
    "from scipy.optimize import fsolve\n",
    "from scipy.interpolate import UnivariateSpline\n",
    "import scipy.integrate\n",
    "import matplotlib.animation as animation\n",
    "from sklearn.metrics import r2_score "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "f243a9ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Test  Cavenati\n",
    "#Parameters for the simulation \n",
    "pressure                              = 3.2                                      # pressure during adsorption step, bar\n",
    "pressure_init                         = 0.1\n",
    "seg                                   = 100\n",
    "T                                     = 303                                   #Solid initial Temperature, K\n",
    "V_tank                                = 6e-5\n",
    "\n",
    "ep                                    = 0.33                                     #Particle porosity \n",
    "ei                                    = 0.45                                     #External porosity\n",
    "et                                    = ep+ei                                     #Total porosity\n",
    "L                                     = 0.83                                    #Length of the bed, m\n",
    "D                                     = 0.021                                #Diameter of the bed, m\n",
    "rho                                   = 715                                     #Apparent(bulk) density, kg/m3\n",
    "R                                     = 8.314                                   #Universal gas constant, J/(mol K)\n",
    "Flow_rate                             = 1.5*6.72701607326312E-07\n",
    "Flow_rate_purge                       = 0.6*6.72701607326312E-07\n",
    "adsorption_time                       = 100                                     #s\n",
    "\n",
    "y1_feed_pressurization                = 1                                    #Methane molar fraction     \n",
    "y2_feed_pressurization                = 0                                   #Ethane molar fraction\n",
    "y3_feed_pressurization                = 0                                   #Propane molar fraction\n",
    "y4_feed_pressurization                = 0                                     #Butane molar fraction\n",
    "y5_feed_pressurization                = 0                                   #CO2 molar fraction\n",
    "y6_feed_pressurization                = 1 - y1_feed_pressurization  - y2_feed_pressurization - y3_feed_pressurization - y4_feed_pressurization  - y5_feed_pressurization                                 #Nitrogen molar fraction\n",
    "\n",
    "                                                                                #Feed Parameters \n",
    "y1_feed_adsorption                    = y1_feed_pressurization                                    #Methane molar fraction     \n",
    "y2_feed_adsorption                    = y2_feed_pressurization                                  #Ethane molar fraction\n",
    "y3_feed_adsorption                    = y3_feed_pressurization                                    #Propane molar fraction\n",
    "y4_feed_adsorption                    = y4_feed_pressurization                                     #Butane molar fraction\n",
    "y5_feed_adsorption                    = y5_feed_pressurization                                   #CO2 molar fraction\n",
    "y6_feed_adsorption                    = y6_feed_pressurization                                   #Nitrogen molar fraction\n",
    "\n",
    "y1_feed_purge                         = 0\n",
    "y2_feed_purge                         = 0\n",
    "y3_feed_purge                         = 0\n",
    "y4_feed_purge                         = 0\n",
    "y5_feed_purge                         = 0\n",
    "y6_feed_purge                         = 1\n",
    "\n",
    "Pressure_purge                        = 0.1                                     #Bar\n",
    "\n",
    "number_cycles = 46                                                                               ### Isotherm parameters ###\n",
    "# Isotherm parameters\n",
    "IP1_CH4 = 1.079e-6 #kmol/kg\n",
    "IP2_CH4 = 2208 #1/bar\n",
    "IP3_CH4 = 0.916 # 1/bar\n",
    "IP4_CH4 = -130.9 # K\n",
    "\n",
    "  \n",
    "IP1_C2H6 = 0.0027 #kmol/kg\n",
    "IP2_C2H6 = 0.0 #1/bar \n",
    "IP3_C2H6 = 2.66E-04 # 1/bar\n",
    "IP4_C2H6 = 2833.77 # K\n",
    "\n",
    "IP1_C3H8 = 0.0062 #kmol/kg\n",
    "IP2_C3H8 = 0.0 #1/bar \n",
    "IP3_C3H8 = 3.75E-04 # 1/bar\n",
    "IP4_C3H8 = 2795.28 # K\n",
    "\n",
    "IP1_C4H10 = 0.0076 #kmol/kg\n",
    "IP2_C4H10 = 0.0 #1/bar \n",
    "IP3_C4H10 = 0.0015 # 1/bar\n",
    "IP4_C4H10 = 2600 # K\n",
    "\n",
    "IP1_CO2 = 1.127e-5 #kmol/kg\n",
    "IP2_CO2 = 1880 #(kmol/kg)/bar \n",
    "IP3_CO2 = 0.145 # 1/bar\n",
    "IP4_CO2 = 777.9 # K\n",
    "\n",
    "IP1_N2 = 0 #kmol/kg\n",
    "IP2_N2 = 0 #(kmol/kg)/bar \n",
    "IP3_N2 = 0 # 1/bar\n",
    "IP4_N2 = 0 # K\n",
    "\n",
    "V_tank1                             = 1e-3                                                                               #Mass transfer coefficients, 1/s\n",
    "MTC                                   = [0,0,0,0,0,0]                     #CH4,C2H6,C3H8,C4H10,CO2,N2\n",
    "\n",
    "T_initial = 303\n",
    "r_p = 5e-3\n",
    "mu_g = 1.185e-2*0.001\n",
    "k_gas = 2.198e-2\n",
    "cp_gas = 34000 #J/kmol K\n",
    "cp_solid = 878 #J/kgK\n",
    "k_solid = 16 # W/mK\n",
    "dH_ch4 = 16*10**6\n",
    "dH_c2h6 = 25*10**6\n",
    "dH_c3h8 = 28*10**6\n",
    "dH_c4h10 = 30*10**6\n",
    "dH_n2 = 10*10**6\n",
    "dH_co2 = 16*10**6\n",
    "\n",
    "exp = 2.7182\n",
    "\n",
    "#Pump char\n",
    "Qmax_pump = 0.000149805471647187\n",
    "b_pump = 768250.064952021\n",
    "n_pump = 0.970650968902519\n",
    "k_pump = 3.2633941869566 \n",
    "\n",
    "B_mu = 0.001\n",
    "K = 0.058"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "e7f0da50",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Staarted cycle 1\n",
      "FINIIIISHED Press\n",
      "Staarted\n",
      "FINIIIISHED ads\n",
      "Staarted\n",
      "Finished DESORP\n",
      "Staarted\n",
      "Finished PURGE\n",
      "Staarted cycle 2\n",
      "FINIIIISHED Press\n",
      "Staarted\n",
      "FINIIIISHED ads\n",
      "Staarted\n",
      "Finished DESORP\n",
      "Staarted\n",
      "Finished PURGE\n",
      "Staarted cycle 3\n",
      "FINIIIISHED Press\n",
      "Staarted\n",
      "FINIIIISHED ads\n",
      "Staarted\n",
      "Finished DESORP\n",
      "Staarted\n",
      "Finished PURGE\n",
      "Staarted cycle 4\n",
      "FINIIIISHED Press\n",
      "Staarted\n",
      "FINIIIISHED ads\n",
      "Staarted\n",
      "Finished DESORP\n",
      "Staarted\n",
      "Finished PURGE\n",
      "Staarted cycle 5\n",
      "FINIIIISHED Press\n",
      "Staarted\n",
      "FINIIIISHED ads\n",
      "Staarted\n",
      "Finished DESORP\n",
      "Staarted\n",
      "Finished PURGE\n",
      "Staarted cycle 6\n",
      "FINIIIISHED Press\n",
      "Staarted\n",
      "FINIIIISHED ads\n",
      "Staarted\n",
      "Finished DESORP\n",
      "Staarted\n",
      "Finished PURGE\n",
      "Staarted cycle 7\n",
      "FINIIIISHED Press\n",
      "Staarted\n",
      "FINIIIISHED ads\n",
      "Staarted\n",
      "Finished DESORP\n",
      "Staarted\n",
      "Finished PURGE\n",
      "Staarted cycle 8\n",
      "FINIIIISHED Press\n",
      "Staarted\n",
      "FINIIIISHED ads\n",
      "Staarted\n",
      "Finished DESORP\n",
      "Staarted\n",
      "Finished PURGE\n",
      "Staarted cycle 9\n",
      "FINIIIISHED Press\n",
      "Staarted\n",
      "FINIIIISHED ads\n",
      "Staarted\n",
      "Finished DESORP\n",
      "Staarted\n",
      "Finished PURGE\n",
      "Staarted cycle 10\n",
      "FINIIIISHED Press\n",
      "Staarted\n",
      "FINIIIISHED ads\n",
      "Staarted\n",
      "Finished DESORP\n",
      "Staarted\n",
      "Finished PURGE\n",
      "Staarted cycle 11\n",
      "FINIIIISHED Press\n",
      "Staarted\n",
      "FINIIIISHED ads\n",
      "Staarted\n",
      "Finished DESORP\n",
      "Staarted\n",
      "Finished PURGE\n",
      "Staarted cycle 12\n",
      "FINIIIISHED Press\n",
      "Staarted\n",
      "FINIIIISHED ads\n",
      "Staarted\n",
      "Finished DESORP\n",
      "Staarted\n",
      "Finished PURGE\n",
      "Staarted cycle 13\n",
      "FINIIIISHED Press\n",
      "Staarted\n",
      "FINIIIISHED ads\n",
      "Staarted\n",
      "Finished DESORP\n",
      "Staarted\n",
      "Finished PURGE\n",
      "Staarted cycle 14\n",
      "FINIIIISHED Press\n",
      "Staarted\n",
      "FINIIIISHED ads\n",
      "Staarted\n",
      "Finished DESORP\n",
      "Staarted\n",
      "Finished PURGE\n",
      "Staarted cycle 15\n",
      "FINIIIISHED Press\n",
      "Staarted\n",
      "FINIIIISHED ads\n",
      "Staarted\n",
      "Finished DESORP\n",
      "Staarted\n",
      "Finished PURGE\n",
      "Staarted cycle 16\n",
      "FINIIIISHED Press\n",
      "Staarted\n",
      "FINIIIISHED ads\n",
      "Staarted\n",
      "Finished DESORP\n",
      "Staarted\n",
      "Finished PURGE\n",
      "Staarted cycle 17\n",
      "FINIIIISHED Press\n",
      "Staarted\n",
      "FINIIIISHED ads\n",
      "Staarted\n",
      "Finished DESORP\n",
      "Staarted\n",
      "Finished PURGE\n",
      "Staarted cycle 18\n",
      "FINIIIISHED Press\n",
      "Staarted\n",
      "FINIIIISHED ads\n",
      "Staarted\n",
      "Finished DESORP\n",
      "Staarted\n",
      "Finished PURGE\n",
      "Staarted cycle 19\n",
      "FINIIIISHED Press\n",
      "Staarted\n",
      "FINIIIISHED ads\n",
      "Staarted\n",
      "Finished DESORP\n",
      "Staarted\n",
      "Finished PURGE\n",
      "Staarted cycle 20\n",
      "FINIIIISHED Press\n",
      "Staarted\n",
      "FINIIIISHED ads\n",
      "Staarted\n",
      "Finished DESORP\n",
      "Staarted\n",
      "Finished PURGE\n",
      "Staarted cycle 21\n",
      "FINIIIISHED Press\n",
      "Staarted\n",
      "FINIIIISHED ads\n",
      "Staarted\n",
      "Finished DESORP\n",
      "Staarted\n",
      "Finished PURGE\n",
      "Staarted cycle 22\n",
      "FINIIIISHED Press\n",
      "Staarted\n",
      "FINIIIISHED ads\n",
      "Staarted\n",
      "Finished DESORP\n",
      "Staarted\n",
      "Finished PURGE\n",
      "Staarted cycle 23\n",
      "FINIIIISHED Press\n",
      "Staarted\n",
      "FINIIIISHED ads\n",
      "Staarted\n",
      "Finished DESORP\n",
      "Staarted\n",
      "Finished PURGE\n",
      "Staarted cycle 24\n",
      "FINIIIISHED Press\n",
      "Staarted\n",
      "FINIIIISHED ads\n",
      "Staarted\n",
      "Finished DESORP\n",
      "Staarted\n",
      "Finished PURGE\n",
      "Staarted cycle 25\n",
      "FINIIIISHED Press\n",
      "Staarted\n",
      "FINIIIISHED ads\n",
      "Staarted\n",
      "Finished DESORP\n",
      "Staarted\n",
      "Finished PURGE\n",
      "Staarted cycle 26\n",
      "FINIIIISHED Press\n",
      "Staarted\n",
      "FINIIIISHED ads\n",
      "Staarted\n",
      "Finished DESORP\n",
      "Staarted\n",
      "Finished PURGE\n",
      "Staarted cycle 27\n",
      "FINIIIISHED Press\n",
      "Staarted\n",
      "FINIIIISHED ads\n",
      "Staarted\n",
      "Finished DESORP\n",
      "Staarted\n",
      "Finished PURGE\n",
      "Staarted cycle 28\n",
      "FINIIIISHED Press\n",
      "Staarted\n",
      "FINIIIISHED ads\n",
      "Staarted\n",
      "Finished DESORP\n",
      "Staarted\n",
      "Finished PURGE\n",
      "Staarted cycle 29\n",
      "FINIIIISHED Press\n",
      "Staarted\n",
      "FINIIIISHED ads\n",
      "Staarted\n",
      "Finished DESORP\n",
      "Staarted\n",
      "Finished PURGE\n",
      "Staarted cycle 30\n",
      "FINIIIISHED Press\n",
      "Staarted\n",
      "FINIIIISHED ads\n",
      "Staarted\n",
      "Finished DESORP\n",
      "Staarted\n",
      "Finished PURGE\n",
      "Staarted cycle 31\n",
      "FINIIIISHED Press\n",
      "Staarted\n",
      "FINIIIISHED ads\n",
      "Staarted\n",
      "Finished DESORP\n",
      "Staarted\n",
      "Finished PURGE\n",
      "Staarted cycle 32\n",
      "FINIIIISHED Press\n",
      "Staarted\n",
      "FINIIIISHED ads\n",
      "Staarted\n",
      "Finished DESORP\n",
      "Staarted\n",
      "Finished PURGE\n",
      "Staarted cycle 33\n",
      "FINIIIISHED Press\n",
      "Staarted\n",
      "FINIIIISHED ads\n",
      "Staarted\n",
      "Finished DESORP\n",
      "Staarted\n",
      "Finished PURGE\n",
      "Staarted cycle 34\n",
      "FINIIIISHED Press\n",
      "Staarted\n",
      "FINIIIISHED ads\n",
      "Staarted\n",
      "Finished DESORP\n",
      "Staarted\n",
      "Finished PURGE\n",
      "Staarted cycle 35\n",
      "FINIIIISHED Press\n",
      "Staarted\n",
      "FINIIIISHED ads\n",
      "Staarted\n",
      "Finished DESORP\n",
      "Staarted\n",
      "Finished PURGE\n",
      "Staarted cycle 36\n",
      "FINIIIISHED Press\n",
      "Staarted\n",
      "FINIIIISHED ads\n",
      "Staarted\n",
      "Finished DESORP\n",
      "Staarted\n",
      "Finished PURGE\n",
      "Staarted cycle 37\n",
      "FINIIIISHED Press\n",
      "Staarted\n",
      "FINIIIISHED ads\n",
      "Staarted\n",
      "Finished DESORP\n",
      "Staarted\n",
      "Finished PURGE\n",
      "Staarted cycle 38\n",
      "FINIIIISHED Press\n",
      "Staarted\n",
      "FINIIIISHED ads\n",
      "Staarted\n",
      "Finished DESORP\n",
      "Staarted\n",
      "Finished PURGE\n",
      "Staarted cycle 39\n",
      "FINIIIISHED Press\n",
      "Staarted\n",
      "FINIIIISHED ads\n",
      "Staarted\n",
      "Finished DESORP\n",
      "Staarted\n",
      "Finished PURGE\n",
      "Staarted cycle 40\n",
      "FINIIIISHED Press\n",
      "Staarted\n",
      "FINIIIISHED ads\n",
      "Staarted\n",
      "Finished DESORP\n",
      "Staarted\n",
      "Finished PURGE\n",
      "Staarted cycle 41\n",
      "FINIIIISHED Press\n",
      "Staarted\n",
      "FINIIIISHED ads\n",
      "Staarted\n",
      "Finished DESORP\n",
      "Staarted\n",
      "Finished PURGE\n",
      "Staarted cycle 42\n",
      "FINIIIISHED Press\n",
      "Staarted\n",
      "FINIIIISHED ads\n",
      "Staarted\n",
      "Finished DESORP\n",
      "Staarted\n",
      "Finished PURGE\n",
      "Staarted cycle 43\n",
      "FINIIIISHED Press\n",
      "Staarted\n",
      "FINIIIISHED ads\n",
      "Staarted\n",
      "Finished DESORP\n",
      "Staarted\n",
      "Finished PURGE\n",
      "Staarted cycle 44\n",
      "FINIIIISHED Press\n",
      "Staarted\n",
      "FINIIIISHED ads\n",
      "Staarted\n",
      "Finished DESORP\n",
      "Staarted\n",
      "Finished PURGE\n",
      "Staarted cycle 45\n",
      "FINIIIISHED Press\n",
      "Staarted\n",
      "FINIIIISHED ads\n",
      "Staarted\n",
      "Finished DESORP\n",
      "Staarted\n",
      "Finished PURGE\n",
      "Staarted cycle 46\n",
      "FINIIIISHED Press\n",
      "Staarted\n",
      "FINIIIISHED ads\n",
      "Staarted\n",
      "Finished DESORP\n",
      "Staarted\n",
      "Finished PURGE\n"
     ]
    }
   ],
   "source": [
    "for cycle_count in range(0,number_cycles):\n",
    "    print(f\"Staarted cycle {cycle_count+1}\")\n",
    "    MTC = [5.7E-5,0.025]\n",
    "    m = GEKKO(remote = False)  # create GEKKO model\n",
    "    tf = 70\n",
    "    nt = int(tf/1) + 1\n",
    "    m.time = np.linspace(0,tf,nt)\n",
    "    if cycle_count == 0:\n",
    "        y1_feed = 0.55\n",
    "        y2_feed = 0.45\n",
    "        y3_feed = 0\n",
    "        L_seg = L/seg\n",
    "        Area = np.pi * D**2/4\n",
    "        P = [m.Var(pressure_init) for i in range(seg)] \n",
    "        y1 =  [m.Var(0) for i in range(seg)] \n",
    "        y2 =  [m.Var(0) for i in range(seg)]\n",
    "        y3 =  [m.Var(1) for i in range(seg)]\n",
    "        v =  [m.Var(0) for i in range(seg)]\n",
    "        q1 = [m.Var(0) for i in range(seg)]\n",
    "        q2 = [m.Var(0) for i in range(seg)]\n",
    "        P1 = m.Var(pressure_init)\n",
    "        u = m.Var(Flow_rate*R*T/(Area*P1.value*10**2))\n",
    "        P2 = m.Var(pressure_init)\n",
    "    else:\n",
    "        y1_feed = 0.55\n",
    "        y2_feed = 0.45\n",
    "        y3_feed = 0\n",
    "        L_seg = L/seg\n",
    "        Area = np.pi * D**2/4\n",
    "        P = [m.Var(pressure_after_purge[-1][i]) for i in range(seg)] \n",
    "        y1 =  [m.Var(met_after_purge[-1][i]) for i in range(seg)] \n",
    "        y2 =  [m.Var(co2_after_purge[-1][i]) for i in range(seg)]\n",
    "        y3 =  [m.Var(n2_after_purge[-1][i]) for i in range(seg)]\n",
    "        v =  [m.Var(0) for i in range(seg)]\n",
    "        q1 = [m.Var(q1_after_purge[-1][i]) for i in range(seg)]\n",
    "        q2 = [m.Var(q2_after_purge[-1][i]) for i in range(seg)]\n",
    "        P1 = m.Var(pressure_after_purge[-1][0])\n",
    "        u = m.Var(Flow_rate*R*T/(Area*pressure_after_purge[-1][0]*10**2))\n",
    "        P2 = m.Var(pressure_after_purge[-1][seg-1])\n",
    "\n",
    "    m.Equation(u == Flow_rate*R*T/(Area*P1*10**2))\n",
    "    m.Equation(P1.dt() ==  (P[0]/(V_tank))*(Flow_rate*10**(-2)*R*T/(P[0])-v[0]*Area))\n",
    "    #Pressure Equation\n",
    "    m.Equation(P2.dt() ==  (P2/(V_tank))*(v[seg-1]*Area))\n",
    "\n",
    "    m.Equation(P[0]*y1[0].dt()+y1[0]*P[0].dt() == - (u/(et))*(y1[0]*(P[0]-P1)/(L_seg)+P[0]*(y1[0]-y1_feed)/(L_seg)) \n",
    "               - (P[0]*y1[0]/(et))*((v[0]-u)/(L_seg)) - (rho/(10**2*et))*R*T*q1[0].dt())\n",
    "\n",
    "    m.Equation([P[i]*y1[i].dt()+y1[i]*P[i].dt() == - (v[i]/(et))*(y1[i]*(P[i]-P[i-1])/(L_seg)+P[i]*(y1[i]-y1[i-1])/(L_seg)) \n",
    "               - (P[i]*y1[i]/(et))*((v[i]-v[i-1])/(L_seg)) - (rho/(10**2*et))*R*T*q1[i].dt() for i in range(1,seg-1)])\n",
    "\n",
    "    m.Equation(P[seg-1]*y1[seg-1].dt()+y1[seg-1]*P[seg-1].dt() == - (v[seg-1]/(et))*(y1[seg-1]*(P[seg-1]-P[seg-2])/(L_seg)\n",
    "                                                                                     +P[seg-1]*(y1[seg-1]-y1[seg-2])/(L_seg)) \n",
    "               - (P[seg-1]*y1[seg-1]/(et))*((v[seg-1]-v[seg-2])/(L_seg)) - (rho/(10**2*et))*R*T*q1[seg-1].dt())\n",
    "\n",
    "\n",
    "    m.Equation(P[0]*y2[0].dt()+y2[0]*P[0].dt() == - (u/(et))*(y2[0]*(P[0]-P1)/(L_seg)+P[0]*(y2[0]-y2_feed)/(L_seg)) \n",
    "               - (P[0]*y2[0]/(et))*((v[0]-u)/(L_seg)) - (rho/(10**2*et))*R*T*q2[0].dt())\n",
    "\n",
    "    m.Equation([P[i]*y2[i].dt()+y2[i]*P[i].dt() == - (v[i]/(et))*(y2[i]*(P[i]-P[i-1])/(L_seg)+P[i]*(y2[i]-y2[i-1])/(L_seg)) \n",
    "               - (P[i]*y2[i]/(et))*((v[i]-v[i-1])/(L_seg)) - (rho/(10**2*et))*R*T*q2[i].dt() for i in range(1,seg-1)])\n",
    "\n",
    "    m.Equation(P[seg-1]*y2[seg-1].dt()+y2[seg-1]*P[seg-1].dt() == - (v[seg-1]/(et))*(y2[seg-1]*(P[seg-1]-P[seg-2])/(L_seg)\n",
    "                +P[seg-1]*(y2[seg-1]-y2[seg-2])/(L_seg)) \n",
    "               - (P[seg-1]**y2[seg-1]/(et))*((v[seg-1]-v[seg-2])/(L_seg)) - (rho/(10**2*et))*R*T*q2[seg-1].dt())\n",
    "\n",
    "    m.Equation(P[0]*y3[0].dt()+y3[0]*P[0].dt() == - (u/(et))*(y3[0]*(P[0]-P1)/(L_seg)+P[0]*(y3[0]-y3_feed)/(L_seg)) \n",
    "               - (P[0]*y3[0]/(et))*((v[0]-u)/(L_seg)))\n",
    "\n",
    "    m.Equation([P[i]*y3[i].dt()+y3[i]*P[i].dt() == - (v[i]/(et))*(y3[i]*(P[i]-P[i-1])/(L_seg)+P[i]*(y3[i]-y3[i-1])/(L_seg)) \n",
    "               - (P[i]*y3[i]/(et))*((v[i]-v[i-1])/(L_seg)) for i in range(1,seg-1)])\n",
    "\n",
    "    m.Equation(P[seg-1]*y3[seg-1].dt()+y3[seg-1]*P[seg-1].dt() == - (v[seg-1]/(et))*(y3[seg-1]*(P[seg-1]-P[seg-2])/(L_seg)\n",
    "                                                                                      +P[seg-1]*y3[seg-1]*(y3[seg-1]-y3[seg-2])/(L_seg)) \n",
    "               - (P[seg-1]*y3[seg-1]/(et))*((v[seg-1]-v[seg-2])/(L_seg)))\n",
    "\n",
    "\n",
    "    m.Equation(v[0] == -(1/K)*(P[1]-P[0])/(L_seg))\n",
    "    m.Equation([v[i] == -(1/K)*(P[i+1]-P[i])/(L_seg) for i in range(1,seg-1)])\n",
    "    m.Equation(v[seg-1] == -(1/K)*(P2 - P[seg-1])/(L_seg))\n",
    "\n",
    "    m.Equation([y1[i] +  y2[i]  == 1 - y3[i] for i in range(0,seg)])\n",
    "\n",
    "    m.Equation([q1[i].dt() == -MTC[0]*(q1[i] - (IP1_CH4*y1[i]*P[i]*exp**(IP2_CH4/T)/(1+y1[i]*P[i]*IP3_CH4*exp**(IP4_CH4/T)))) for i in range(0,seg)])\n",
    "    m.Equation([q2[i].dt() == -MTC[1]*(q2[i] - (IP1_CO2*y2[i]*P[i]*exp**(IP2_CO2/T)/(1+y2[i]*P[i]*IP3_CO2*exp**(IP4_CO2/T)))) for i in range(0,seg)])\n",
    "\n",
    "    # simulation\n",
    "    m.options.IMODE = 7\n",
    "    m.solve(disp = False)\n",
    "    print(\"FINIIIISHED Press\")\n",
    "    \n",
    "    if cycle_count == number_cycles - 1:\n",
    "        \n",
    "        velocity_out = np.array(v[seg-1])\n",
    "        pressure_out = np.array(P[seg-1])\n",
    "        methane_out = np.array(y1[seg-1])\n",
    "        co2_out = np.array(y2[seg-1])\n",
    "        \n",
    "        flow_methane =  velocity_out*Area*methane_out*pressure_out*10**8/(R*T)\n",
    "        flow_co2 = velocity_out*Area*co2_out*pressure_out*10**8/(R*T)\n",
    "        \n",
    "        \n",
    "        old_indices = np.arange(0,len(flow_methane))\n",
    "        new_length = 100\n",
    "        new_indices = np.linspace(0,len(flow_methane)-1,new_length)\n",
    "        spl = UnivariateSpline(old_indices,flow_methane,k=1,s=0)\n",
    "        flow_methane  = np.abs(spl(new_indices))\n",
    "        \n",
    "        old_indices = np.arange(0,len(flow_co2))\n",
    "        new_length = 100\n",
    "        new_indices = np.linspace(0,len(flow_co2)-1,new_length)\n",
    "        spl = UnivariateSpline(old_indices,flow_co2,k=1,s=0)\n",
    "        flow_co2  = np.abs(spl(new_indices))\n",
    "        \n",
    "        old_indices = np.arange(0,len(pressure_out))\n",
    "        new_length = 100\n",
    "        new_indices = np.linspace(0,len(pressure_out)-1,new_length)\n",
    "        spl = UnivariateSpline(old_indices,pressure_out,k=1,s=0)\n",
    "        pressure_out  = np.abs(spl(new_indices))\n",
    "        \n",
    "        old_indices = np.arange(0,len(m.time))\n",
    "        new_length = 100\n",
    "        new_indices = np.linspace(0,len(m.time)-1,new_length)\n",
    "        spl = UnivariateSpline(old_indices,m.time,k=1,s=0)\n",
    "        time  = np.abs(spl(new_indices))\n",
    "        \n",
    "        np_array_rows = time,pressure_out,flow_methane,flow_co2\n",
    "        \n",
    "        with open('results12.csv','a') as csvfile:\n",
    "            np.savetxt(csvfile, np_array_rows , delimiter=',',header='Press',fmt='%s', comments='')\n",
    "        \n",
    "        \n",
    "    met_after_press = np.empty((seg,len(m.time)))\n",
    "    for i in range(seg):\n",
    "        met_after_press[i] = np.array(y1[i].value)\n",
    "    met_after_press= met_after_press.T\n",
    "\n",
    "    co2_after_press = np.empty((seg,len(m.time)))\n",
    "    for i in range(seg):\n",
    "        co2_after_press[i] = np.array(y2[i].value)\n",
    "    co2_after_press= co2_after_press.T\n",
    "\n",
    "\n",
    "    n2_after_press = np.empty((seg,len(m.time)))\n",
    "    for i in range(seg):\n",
    "        n2_after_press[i] = np.array(y3[i].value)\n",
    "    n2_after_press= n2_after_press.T\n",
    "\n",
    "    q1_after_press = np.empty((seg,len(m.time)))\n",
    "    for i in range(seg):\n",
    "        q1_after_press[i] = np.array(q1[i].value)\n",
    "    q1_after_press= q1_after_press.T\n",
    "\n",
    "    q2_after_press = np.empty((seg,len(m.time)))\n",
    "    for i in range(seg):\n",
    "        q2_after_press[i] = np.array(q2[i].value)\n",
    "    q2_after_press= q2_after_press.T\n",
    "    \n",
    "    pressure_after_adsorption = np.empty((seg,len(m.time)))\n",
    "    for i in range(seg):\n",
    "        pressure_after_adsorption[i] = np.array(P[i].value)\n",
    "    pressure_after_adsorption = pressure_after_adsorption.T\n",
    "\n",
    "    print(\"Staarted\")\n",
    "    MTC = [5.7E-5,0.025]\n",
    "    m = GEKKO(remote = False)  # create GEKKO model\n",
    "    tf =  100\n",
    "    nt = int(tf/10) + 1\n",
    "    m.time = np.linspace(0,tf,nt)\n",
    "    y1_feed = 0.55\n",
    "    y2_feed = 0.45\n",
    "    y3_feed = 0\n",
    "    L_seg = L/seg\n",
    "    Area = np.pi * D**2/4\n",
    "    P = [m.Var(pressure) for i in range(seg)] \n",
    "    y1 =  [m.Var(met_after_press[-1][i]) for i in range(seg)] \n",
    "    y2 =  [m.Var(co2_after_press[-1][i]) for i in range(seg)]\n",
    "    y3 =  [m.Var(n2_after_press[-1][i]) for i in range(seg)]\n",
    "    u_initial = Flow_rate*R*T/(Area*pressure*10**2)\n",
    "    v =  [m.Var(0) for i in range(seg)]\n",
    "    q1 = [m.Var(q1_after_press[-1][i]) for i in range(seg)]\n",
    "    q2 = [m.Var(q2_after_press[-1][i]) for i in range(seg)]\n",
    "\n",
    "    # first segment Pressure 1\n",
    "    m.Equation(P[0]*y1[0].dt()+y1[0]*P[0].dt() == - (u_initial/(et))*(y1[0]*(P[0]-pressure)/(L_seg)+P[0]*(y1[0]-y1_feed)/(L_seg)) \n",
    "               - (P[0]*y1[0]/(et))*((v[0]-u_initial)/(L_seg)) - (rho/(10**2*et))*R*T*q1[0].dt())\n",
    "    m.Equation([P[i]*y1[i].dt()+y1[i]*P[i].dt() == - (v[i]/(et))*(y1[i]*(P[i]-P[i-1])/(L_seg)+P[i]*(y1[i]-y1[i-1])/(L_seg)) \n",
    "               - (P[i]*y1[i]/(et))*((v[i]-v[i-1])/(L_seg)) - (rho/(10**2*et))*R*T*q1[i].dt() for i in range(1,seg-1)])\n",
    "    m.Equation(P[seg-1]*y1[seg-1].dt()+y1[seg-1]*P[seg-1].dt() == - (v[seg-1]/(et))*(y1[seg-1]*(P[seg-1]-P[seg-2])/(L_seg)+P[seg-1]*(y1[seg-1]-y1[seg-2])/(L_seg)) \n",
    "               - (P[seg-1]*y1[seg-1]/(et))*((v[seg-1]-v[seg-2])/(L_seg)) - (rho/(10**2*et))*R*T*q1[seg-1].dt())\n",
    "\n",
    "\n",
    "    m.Equation(P[0]*y2[0].dt()+y2[0]*P[0].dt() == - (u_initial/(et))*(y2[0]*(P[0]-pressure)/(L_seg)+P[0]*(y2[0]-y2_feed)/(L_seg)) \n",
    "               - (P[0]*y2[0]/(et))*((v[0]-u_initial)/(L_seg)) - (rho/(10**2*et))*R*T*q2[0].dt())\n",
    "\n",
    "    m.Equation([P[i]*y2[i].dt()+y2[i]*P[i].dt() == - (v[i]/(et))*(y2[i]*(P[i]-P[i-1])/(L_seg)+P[i]*(y2[i]-y2[i-1])/(L_seg)) \n",
    "               - (P[i]*y2[i]/(et))*((v[i]-v[i-1])/(L_seg)) - (rho/(10**2*et))*R*T*q2[i].dt() for i in range(1,seg-1)])\n",
    "\n",
    "    m.Equation(P[seg-1]*y2[seg-1].dt()+y2[seg-1]*P[seg-1].dt() == - (v[seg-1]/(et))*(y2[seg-1]*(P[seg-1]-P[seg-2])/(L_seg)\n",
    "                +P[seg-1]*(y2[seg-1]-y2[seg-2])/(L_seg)) \n",
    "               - (P[seg-1]**y2[seg-1]/(et))*((v[seg-1]-v[seg-2])/(L_seg)) - (rho/(10**2*et))*R*T*q2[seg-1].dt())\n",
    "\n",
    "    m.Equation(P[0]*y3[0].dt()+y3[0]*P[0].dt() == - (u_initial/(et))*(y3[0]*(P[0]-pressure)/(L_seg)+P[0]*(y3[0]-y3_feed)/(L_seg)) \n",
    "               - (P[0]*y3[0]/(et))*((v[0]-u_initial)/(L_seg)))\n",
    "    m.Equation([P[i]*y3[i].dt()+y3[i]*P[i].dt() == - (v[i]/(et))*(y3[i]*(P[i]-P[i-1])/(L_seg)+P[i]*(y3[i]-y3[i-1])/(L_seg)) \n",
    "               - (P[i]*y3[i]/(et))*((v[i]-v[i-1])/(L_seg)) for i in range(1,seg-1)])\n",
    "    m.Equation(P[seg-1]*y3[seg-1].dt()+y3[seg-1]*P[seg-1].dt() == - (v[seg-1]/(et))*(y3[seg-1]*(P[seg-1]-P[seg-2])/(L_seg)\n",
    "                                                                                      +P[seg-1]*y3[seg-1]*(y3[seg-1]-y3[seg-2])/(L_seg)) \n",
    "               - (P[seg-1]*y3[seg-1]/(et))*((v[seg-1]-v[seg-2])/(L_seg)))\n",
    "\n",
    "    m.Equation(v[0] == -(1/K)*(P[0]-pressure)/(L_seg))\n",
    "    m.Equation([v[i] == -(1/K)*(P[i]-P[i-1])/(L_seg) for i in range(1,seg-1)])\n",
    "    m.Equation(v[seg-1] == -(1/K)*(P[seg-1]-P[seg-2])/(L_seg))\n",
    "\n",
    "    m.Equation([y1[i] +  y2[i]  == 1 - y3[i] for i in range(0,seg)])\n",
    "\n",
    "    m.Equation([q1[i].dt() == -MTC[0]*(q1[i] - (IP1_CH4*y1[i]*P[i]*exp**(IP2_CH4/T)/(1+y1[i]*P[i]*IP3_CH4*exp**(IP4_CH4/T)))) for i in range(0,seg)])\n",
    "    m.Equation([q2[i].dt() == -MTC[1]*(q2[i] - (IP1_CO2*y2[i]*P[i]*exp**(IP2_CO2/T)/(1+y2[i]*P[i]*IP3_CO2*exp**(IP4_CO2/T)))) for i in range(0,seg)])\n",
    "\n",
    "    # simulation\n",
    "    m.options.IMODE = 7\n",
    "    m.solve(disp = False)\n",
    "    print(\"FINIIIISHED ads\")\n",
    "\n",
    "    met_after_ads= np.empty((seg,len(m.time)))\n",
    "    for i in range(seg):\n",
    "        met_after_ads[i] = np.array(y1[i].value)\n",
    "    met_after_ads= met_after_ads.T\n",
    "\n",
    "    co2_after_ads = np.empty((seg,len(m.time)))\n",
    "    for i in range(seg):\n",
    "        co2_after_ads[i] = np.array(y2[i].value)\n",
    "    co2_after_ads= co2_after_ads.T\n",
    "\n",
    "\n",
    "    n2_after_ads = np.empty((seg,len(m.time)))\n",
    "    for i in range(seg):\n",
    "        n2_after_ads[i] = np.array(y3[i].value)\n",
    "    n2_after_ads= n2_after_ads.T\n",
    "\n",
    "    q1_after_ads = np.empty((seg,len(m.time)))\n",
    "    for i in range(seg):\n",
    "        q1_after_ads[i] = np.array(q1[i].value)\n",
    "    q1_after_ads= q1_after_ads.T\n",
    "\n",
    "    q2_after_ads = np.empty((seg,len(m.time)))\n",
    "    for i in range(seg):\n",
    "        q2_after_ads[i] = np.array(q2[i].value)\n",
    "    q2_after_ads= q2_after_ads.T\n",
    "    if cycle_count == number_cycles - 1:\n",
    "        velocity_out = np.array(v[seg-1])\n",
    "        pressure_out = np.array(P[seg-1])\n",
    "        methane_out = np.array(y1[seg-1])\n",
    "        co2_out = np.array(y2[seg-1])\n",
    "        flow_methane =  velocity_out*Area*methane_out*pressure_out*10**8/(R*T)\n",
    "        flow_co2 = velocity_out*Area*co2_out*pressure_out*10**8/(R*T)\n",
    "        \n",
    "        \n",
    "        old_indices = np.arange(0,len(flow_methane))\n",
    "        new_length = 100\n",
    "        new_indices = np.linspace(0,len(flow_methane)-1,new_length)\n",
    "        spl = UnivariateSpline(old_indices,flow_methane,k=1,s=0)\n",
    "        flow_methane  = np.abs(spl(new_indices))\n",
    "        \n",
    "        old_indices = np.arange(0,len(flow_co2))\n",
    "        new_length = 100\n",
    "        new_indices = np.linspace(0,len(flow_co2)-1,new_length)\n",
    "        spl = UnivariateSpline(old_indices,flow_co2,k=1,s=0)\n",
    "        flow_co2  = np.abs(spl(new_indices))\n",
    "        \n",
    "        old_indices = np.arange(0,len(pressure_out))\n",
    "        new_length = 100\n",
    "        new_indices = np.linspace(0,len(pressure_out)-1,new_length)\n",
    "        spl = UnivariateSpline(old_indices,pressure_out,k=1,s=0)\n",
    "        pressure_out  = np.abs(spl(new_indices))\n",
    "        \n",
    "        old_indices = np.arange(0,len(m.time))\n",
    "        new_length = 100\n",
    "        new_indices = np.linspace(0,len(m.time)-1,new_length)\n",
    "        spl = UnivariateSpline(old_indices,m.time,k=1,s=0)\n",
    "        time  = np.abs(spl(new_indices))\n",
    "        \n",
    "        np_array_rows = time + 70,pressure_out,flow_methane,flow_co2\n",
    "        \n",
    "        with open('results12.csv','a') as csvfile:\n",
    "            np.savetxt(csvfile, np_array_rows , delimiter=',',header='ADS',fmt='%s', comments='')\n",
    "\n",
    "    print(\"Staarted\")\n",
    "    MTC = [5.7E-5,0.025]\n",
    "    m = GEKKO(remote = False)  # create GEKKO model\n",
    "    tf = 100\n",
    "    nt = int(tf/1) + 1\n",
    "    m.time = np.linspace(0,tf,nt)\n",
    "\n",
    "    L_seg = L/seg\n",
    "    Area = np.pi * D**2/4\n",
    "    P = [m.Var(pressure) for i in range(seg)] \n",
    "    y1 =  [m.Var(met_after_ads[-1][i]) for i in range(seg)] \n",
    "    y2 =  [m.Var(co2_after_ads[-1][i]) for i in range(seg)]\n",
    "    y3 =  [m.Var(n2_after_ads[-1][i]) for i in range(seg)]\n",
    "    v =  [m.Var(0) for i in range(seg)]\n",
    "    q1 = [m.Var(q1_after_ads[-1][i]) for i in range(seg)]\n",
    "    q2 = [m.Var(q2_after_ads[-1][i]) for i in range(seg)]\n",
    "    P1 = m.Var(pressure)\n",
    "    P2 = m.Var(pressure)\n",
    "\n",
    "    Q_pump   =  m.Var(0.8 / 3600)\n",
    "    m.Equation(P2.dt() ==  (P2/(V_tank))*(v[seg-1]*Area))\n",
    "    m.Equation(P1.dt() == (-((P[0])*Q_pump)/(V_tank) - (P[0]*v[0]*Area)/(V_tank)))\n",
    "    m.Equation(Q_pump ==  0.0001974563*P[0]**(0.627027))\n",
    "\n",
    "    # first segment Pressure 1\n",
    "    m.Equation(P[0]*y1[0].dt()+y1[0]*P[0].dt() == (1/(K*et))*((y1[0]*P[0]*(P[1] - 2*P[0] + P1)/(L_seg**2))\n",
    "                                                          + P[1]*((P[1]-P[0])/(L_seg))*((y1[1]-y1[0])/(L_seg))\n",
    "                                                          + y1[1]*((P[1]-P[0])/(L_seg))**2) \n",
    "    - ((1/(et*100))*(R*T)*rho*q1[0].dt()))\n",
    "\n",
    "    # middle segments Pressure 1\n",
    "    m.Equations([P[i]*y1[i].dt()+y1[i]*P[i].dt() ==  \\\n",
    "                (1/(K*et))*((y1[i]*P[i]*(P[i+1]-2*P[i]+P[i-1])/(L_seg**2)) + P[i+1]*((P[i+1]-P[i])/(L_seg))*((y1[i+1]-y1[i])/(L_seg)) \n",
    "                + y1[i+1]*((P[i+1]-P[i])/(L_seg))**2) - (1/(et*100))*(R*T)*rho*q1[i].dt()   for i in range(1,seg-1)])\n",
    "\n",
    "    # last segment Pressure 1\n",
    "    m.Equation(P[seg-1]*y1[seg-1].dt()+y1[seg-1]*P[seg-1].dt() == \\\n",
    "            (1/(K*et))*((y1[seg-1]*P[seg-1]*(-P[seg-1]+P[seg-2])/(L_seg**2))) - ((1/(et*100))*(R*T)*rho*q1[seg-1].dt()))\n",
    "\n",
    "    # first segment Pressure 2\n",
    "    m.Equation(P[0]*y2[0].dt()+y2[0]*P[0].dt() == (1/(K*et))*((y2[0]*P[0]*(P[1] - 2*P[0] + P1)/(L_seg**2)) + P[1]*((P[1]-P[0])/(L_seg))*((y2[1]-y2[0])/(L_seg)) \n",
    "    + y2[1]*((P[1]-P[0])/(L_seg))**2) - ((1/(et*100))*(R*T)*rho*q2[0].dt())) \n",
    "\n",
    "    # middle segments Pressure 2\n",
    "    m.Equations([P[i]*y2[i].dt()+y2[i]*P[i].dt() ==  \\\n",
    "                (1/(K*et))*((y2[i]*P[i]*(P[i+1]-2*P[i]+P[i-1])/(L_seg**2)) +  P[i+1]*((P[i+1]-P[i])/(L_seg))*((y2[i+1]-y2[i])/(L_seg)) \n",
    "                + y2[i+1]*((P[i+1]-P[i])/(L_seg))**2) - ((1/(et*100))*(R*T))*rho*q2[i].dt()   for i in range(1,seg-1)])\n",
    "\n",
    "    # last segment Pressure \n",
    "    m.Equation(P[seg-1]*y2[seg-1].dt()+y2[seg-1]*P[seg-1].dt() == \\\n",
    "            (1/(K*et))*((y2[seg-1]*P[seg-1]*(-P[seg-1]+P[seg-2])/(L_seg**2))) - ((1/(et*100))*(R*T))*rho*q2[seg-1].dt())\n",
    "\n",
    "\n",
    "\n",
    "    # first segment Pressure 3\n",
    "    m.Equation(P[0]*y3[0].dt()+y3[0]*P[0].dt() == (1/(K*et))*((y3[0]*P[0]*(P[1] - 2*P[0] + P1)/(L_seg**2)) \n",
    "    + P[1]*((P[1]-P[0])/(L_seg))*((y3[1]-y3[0])/(L_seg)) + y3[1]*((P[1]-P[0])/(L_seg))**2))\n",
    "\n",
    "    # middle segments Pressure 3\n",
    "    m.Equations([P[i]*y3[i].dt()+y3[i]*P[i].dt() ==  \\\n",
    "                (1/(K*et))*((y3[i]*P[i]*(P[i+1]-2*P[i]+P[i-1])/(L_seg**2)) + P[i+1]*((P[i+1]-P[i])/(L_seg))*((y3[i+1]-y3[i])/(L_seg)) \n",
    "                + y3[i+1]*((P[i+1]-P[i])/(L_seg))**2)  for i in range(1,seg-1)])\n",
    "\n",
    "    # last segment Pressure 3\n",
    "    m.Equation(P[seg-1]*y3[seg-1].dt()+y3[seg-1]*P[seg-1].dt() == \\\n",
    "            (1/(K*et))*((y3[seg-1]*P[seg-1]*(-P[seg-1]+P[seg-2])/(L_seg**2))))\n",
    "\n",
    "    m.Equation(v[0] == -(1/K)*(P[1]-P[0])/(L_seg))\n",
    "    m.Equation([v[i] == -(1/K)*(P[i+1]-P[i])/(L_seg) for i in range(1,seg-1)])\n",
    "    m.Equation(v[seg-1] == -(1/K)*(P2 - P[seg-1])/(L_seg))\n",
    "\n",
    "    m.Equation([y1[i] +  y2[i]  == 1 - y3[i] for i in range(0,seg)])\n",
    "\n",
    "    m.Equation([q1[i].dt() == -MTC[0]*(q1[i] - (IP1_CH4*y1[i]*P[i]*exp**(IP2_CH4/T)/(1+y1[i]*P[i]*IP3_CH4*exp**(IP4_CH4/T)))) for i in range(0,seg)])\n",
    "    m.Equation([q2[i].dt() == -MTC[1]*(q2[i] - (IP1_CO2*y2[i]*P[i]*exp**(IP2_CO2/T)/(1+y2[i]*P[i]*IP3_CO2*exp**(IP4_CO2/T)))) for i in range(0,seg)])\n",
    "\n",
    "\n",
    "    # simulation\n",
    "    m.options.IMODE = 7\n",
    "    m.solve(disp=False)\n",
    "    print(\"Finished DESORP\")\n",
    "    \n",
    "            \n",
    "    met_after_desorp= np.empty((seg,len(m.time)))\n",
    "    for i in range(seg):\n",
    "        met_after_desorp[i] = np.array(y1[i].value)\n",
    "    met_after_desorp= met_after_desorp.T\n",
    "\n",
    "    co2_after_desorp= np.empty((seg,len(m.time)))\n",
    "    for i in range(seg):\n",
    "        co2_after_desorp[i] = np.array(y2[i].value)\n",
    "    co2_after_desorp= co2_after_desorp.T\n",
    "\n",
    "\n",
    "    n2_after_desorp = np.empty((seg,len(m.time)))\n",
    "    for i in range(seg):\n",
    "        n2_after_desorp[i] = np.array(y3[i].value)\n",
    "    n2_after_desorp= n2_after_desorp.T\n",
    "\n",
    "    q1_after_desorp = np.empty((seg,len(m.time)))\n",
    "    for i in range(seg):\n",
    "        q1_after_desorp[i] = np.array(q1[i].value)\n",
    "    q1_after_desorp= q1_after_desorp.T\n",
    "\n",
    "    q2_after_desorp = np.empty((seg,len(m.time)))\n",
    "    for i in range(seg):\n",
    "        q2_after_desorp[i] = np.array(q2[i].value)\n",
    "    q2_after_desorp= q2_after_desorp.T\n",
    "\n",
    "    v_after_desorp = np.empty((seg,len(m.time)))\n",
    "    for i in range(seg):\n",
    "        v_after_desorp[i] = np.array(v[i].value)\n",
    "    v_after_desorp= v_after_desorp.T\n",
    "\n",
    "    pressure_after_desorp = np.empty((seg,len(m.time)))\n",
    "    for i in range(seg):\n",
    "        pressure_after_desorp[i] = np.array(P[i].value)\n",
    "    pressure_after_desorp= pressure_after_desorp.T\n",
    "\n",
    "    v_after_desorp = np.empty((seg,len(m.time)))\n",
    "    for i in range(seg):\n",
    "        v_after_desorp[i] = np.array(v[i].value)\n",
    "    v_after_desorp= v_after_desorp.T\n",
    "    \n",
    "    if cycle_count == number_cycles - 1:\n",
    "        velocity_out = -1*np.array(v[0])\n",
    "        pressure_out = np.array(P[0])\n",
    "        methane_out = np.array(y1[0])\n",
    "        co2_out = np.array(y2[0])\n",
    "        flow_methane =  velocity_out*Area*methane_out*pressure_out*10**8/(R*T)\n",
    "        flow_co2 = velocity_out*Area*co2_out*pressure_out*10**8/(R*T)\n",
    "\n",
    "        \n",
    "        old_indices = np.arange(0,len(flow_methane))\n",
    "        new_length = 100\n",
    "        new_indices = np.linspace(0,len(flow_methane)-1,new_length)\n",
    "        spl = UnivariateSpline(old_indices,flow_methane,k=1,s=0)\n",
    "        flow_methane  = np.abs(spl(new_indices))\n",
    "        \n",
    "        old_indices = np.arange(0,len(flow_co2))\n",
    "        new_length = 100\n",
    "        new_indices = np.linspace(0,len(flow_co2)-1,new_length)\n",
    "        spl = UnivariateSpline(old_indices,flow_co2,k=1,s=0)\n",
    "        flow_co2  = np.abs(spl(new_indices))\n",
    "        \n",
    "        old_indices = np.arange(0,len(pressure_out))\n",
    "        new_length = 100\n",
    "        new_indices = np.linspace(0,len(pressure_out)-1,new_length)\n",
    "        spl = UnivariateSpline(old_indices,pressure_out,k=1,s=0)\n",
    "        pressure_out  = np.abs(spl(new_indices))\n",
    "        \n",
    "        old_indices = np.arange(0,len(m.time))\n",
    "        new_length = 100\n",
    "        new_indices = np.linspace(0,len(m.time)-1,new_length)\n",
    "        spl = UnivariateSpline(old_indices,m.time,k=1,s=0)\n",
    "        time  = np.abs(spl(new_indices))\n",
    "        \n",
    "        np_array_rows = time+170,pressure_out,flow_methane,flow_co2\n",
    "        \n",
    "        with open('results12.csv','a') as csvfile:\n",
    "            np.savetxt(csvfile, np_array_rows , delimiter=',',header='des',fmt='%s', comments='')\n",
    "\n",
    "    print(\"Staarted\")\n",
    "    MTC = [5.7E-5,0.025]\n",
    "    m = GEKKO(remote = False)  # create GEKKO model\n",
    "    tf = 70\n",
    "    nt = int(tf/1) + 1\n",
    "    m.time = np.linspace(0,tf,nt)\n",
    "\n",
    "    y1_feed = 1\n",
    "    y2_feed = 0\n",
    "    y3_feed = 0\n",
    "\n",
    "    L_seg = L/seg\n",
    "    Area = np.pi * D**2/4\n",
    "    P = [m.Var(pressure_after_desorp[-1][i]) for i in range(seg)] \n",
    "    y1 =  [m.Var(met_after_desorp[-1][i]) for i in range(seg)] \n",
    "    y2 =  [m.Var(co2_after_desorp[-1][i]) for i in range(seg)]\n",
    "    y3 =  [m.Var(n2_after_desorp[-1][i]) for i in range(seg)]\n",
    "    v =  [m.Var(v_after_desorp[-1][i]) for i in range(seg)]\n",
    "    q1 = [m.Var(q1_after_desorp[-1][i]) for i in range(seg)]\n",
    "    q2 = [m.Var(q2_after_desorp[-1][i]) for i in range(seg)]\n",
    "    P1 = m.Var(pressure_after_desorp[-1][0])\n",
    "    P2 = m.Var(pressure_after_desorp[-1][-1])\n",
    "\n",
    "    Q_pump   =  m.Var(0.6 / 3600)\n",
    "\n",
    "    m.Equation(Q_pump ==  0.0001974563*P[0]**(0.627027))\n",
    "\n",
    "    m.Equation(P1.dt() == (-((P[0])*Q_pump)/(V_tank) - (P[0]*v[0]*Area)/(V_tank)))\n",
    "\n",
    "    m.Equation(P2.dt() ==  (P[seg-1]/(V_tank))*(Flow_rate_purge*10**(-2)*R*T/(P[seg-1])+v[seg-1]*Area))\n",
    "\n",
    "\n",
    "    #Pressure Equation\n",
    "    # first segment Pressure 1\n",
    "    m.Equation(P[0]*y1[0].dt()+y1[0]*P[0].dt() == (1/(K*et))*((y1[0]*P[0]*(P[1]-2*P[0]+P1)/(L_seg**2))\n",
    "                                                          + (P[1])*((P[1]- P[0])/(L_seg))*((y1[1]-y1[0])/(L_seg)) \n",
    "    + (y1[1])*((P[1]-P[0])/(L_seg))**2) - ((1/(et*100))*(R*T))*rho*q1[0].dt())\n",
    "\n",
    "    # middle segments Pressure 1\n",
    "    m.Equations([P[i]*y1[i].dt()+y1[i]*P[i].dt() ==  \\\n",
    "                (1/(K*et))*((y1[i]*P[i]*(P[i+1]-2*P[i]+P[i-1])/(L_seg**2)) + (P[i+1])*((P[i+1]-P[i])/(L_seg))*((y1[i+1]-y1[i])/(L_seg)) \n",
    "                + (y1[i+1])*((P[i+1]-P[i])/(L_seg))**2) - ((1/(et*100))*(R*T))*rho*q1[i].dt()  for i in range(1,seg-1)])\n",
    "\n",
    "    # last segment Pressure 1\n",
    "    m.Equation(P[seg-1]*y1[seg-1].dt()+y1[seg-1]*P[seg-1].dt() == \\\n",
    "            (1/(K*et))*((y1[seg-1]*P[seg-1]*(P2-2*P[seg-1]+P[seg-2])/(L_seg**2)) \n",
    "            + (P2)*((P2-P[seg-1])/(L_seg))*((y1_feed-y1[seg-1])/(L_seg))\n",
    "                    + (y1_feed)*((P2-P[seg-1])/(L_seg))**2) - ((1/(et*100))*(R*T))*rho*q1[seg-1].dt())\n",
    "\n",
    "    # first segment Pressure 2\n",
    "    m.Equation(P[0]*y2[0].dt()+y2[0]*P[0].dt() == (1/(K*et))*((y2[0]*P[0]*(P[1]-2*P[0]+P1)/(L_seg**2)) \n",
    "    + (P[1])*((P[1]- P[0])/(L_seg))*((y2[1]-y2[0])/(L_seg))\n",
    "                                                          + (y2[1])*((P[1]-P[0])/(L_seg))**2) \n",
    "               - ((1/(et*100))*(R*T))*rho*q2[0].dt())\n",
    "\n",
    "    # middle segments Pressure 2\n",
    "    m.Equations([P[i]*y2[i].dt()+y2[i]*P[i].dt() ==  \\\n",
    "                (1/(K*et))*((y2[i]*P[i]*(P[i+1]-2*P[i]+P[i-1])/(L_seg**2)) \n",
    "                + (P[i+1])*((P[i+1]-P[i])/(L_seg))*((y2[i+1]-y2[i])/(L_seg)) + (y2[i+1])*((P[i+1]-P[i])/(L_seg))**2) \n",
    "                - ((1/(et*100))*(R*T))*rho*q2[i].dt()   for i in range(1,seg-1)])\n",
    "\n",
    "    # last segment Pressure 2\n",
    "    m.Equation(P[seg-1]*y2[seg-1].dt()+y2[seg-1]*P[seg-1].dt() == \\\n",
    "            (1/(K*et))*((y2[seg-1]*P[seg-1]*(P2-2*P[seg-1]+P[seg-2])/(L_seg**2)) \n",
    "            + (P2)*((P2-P[seg-1])/(L_seg))*((y2_feed-y2[seg-1])/(L_seg)) \n",
    "                    + (y2_feed)*((P2-P[seg-1])/(L_seg))**2) - ((1/(et*100))*(R*T))*rho*q2[seg-1].dt())\n",
    "\n",
    "    # first segment Pressure 3\n",
    "    m.Equation(P[0]*y3[0].dt()+y3[0]*P[0].dt() == (1/(K*et))*((y3[0]*P[0]*(P[1]-2*P[0]+P1)/(L_seg**2)) \n",
    "    + (P[1])*((P[1]- P[0])/(L_seg))*((y3[1]-y3[0])/(L_seg)) \n",
    "                                                          + (y3[1])*((P[1]-P[0])/(L_seg))**2) \n",
    "               )\n",
    "\n",
    "    # middle segments Pressure 3\n",
    "    m.Equations([P[i]*y3[i].dt()+y3[i]*P[i].dt() ==  \\\n",
    "                (1/(K*et))*((y3[i]*P[i]*(P[i+1]-2*P[i]+P[i-1])/(L_seg**2)) \n",
    "                        + (P[i+1])*((P[i+1]-P[i])/(L_seg))*((y3[i+1]-y3[i])/(L_seg)) \n",
    "                        + (y3[i+1])*((P[i+1]-P[i])/(L_seg))**2) \n",
    "                 for i in range(1,seg-1)])\n",
    "\n",
    "    # last segment Pressure 3\n",
    "    m.Equation(P[seg-1]*y3[seg-1].dt()+y3[seg-1]*P[seg-1].dt() == \n",
    "            (1/(K*et))*((y3[seg-1]*P[seg-1]*(P2-2*P[seg-1]+P[seg-2])/(L_seg**2))\n",
    "                    + (P2)*((P2-P[seg-1])/(L_seg))*((y3_feed-y3[seg-1])/(L_seg)) + (y3_feed)*((P2-P[seg-1])/(L_seg))**2) \n",
    "            )\n",
    "\n",
    "\n",
    "    m.Equation(v[0] == -(1/K)*(P[1]-P[0])/(L_seg))\n",
    "    m.Equation([v[i] == -(1/K)*(P[i+1]-P[i])/(L_seg) for i in range(1,seg-1)])\n",
    "    m.Equation(v[seg-1] == -(1/K)*(P2 - P[seg-1])/(L_seg))\n",
    "\n",
    "    m.Equation([y1[i] +  y2[i]  == 1 - y3[i] for i in range(0,seg)])\n",
    "\n",
    "    m.Equation([q1[i].dt() == -MTC[0]*(q1[i] - (IP1_CH4*y1[i]*P[i]*exp**(IP2_CH4/T)/(1+y1[i]*P[i]*IP3_CH4*exp**(IP4_CH4/T)))) for i in range(0,seg)])\n",
    "    m.Equation([q2[i].dt() == -MTC[1]*(q2[i] - (IP1_CO2*y2[i]*P[i]*exp**(IP2_CO2/T)/(1+y2[i]*P[i]*IP3_CO2*exp**(IP4_CO2/T)))) for i in range(0,seg)])\n",
    "\n",
    "\n",
    "    # simulation\n",
    "    m.options.IMODE = 7\n",
    "    m.solve(disp = False)\n",
    "\n",
    "    print(\"Finished PURGE\")\n",
    "    \n",
    "    if cycle_count == number_cycles - 1:\n",
    "        velocity_out = -1*np.array(v[0])\n",
    "        pressure_out = np.array(P[0])\n",
    "        methane_out = np.array(y1[0])\n",
    "        co2_out = np.array(y2[0])\n",
    "        flow_methane =  velocity_out*Area*methane_out*pressure_out*10**8/(R*T)\n",
    "        flow_co2 = velocity_out*Area*co2_out*pressure_out*10**8/(R*T)\n",
    "        \n",
    "        \n",
    "        old_indices = np.arange(0,len(flow_methane))\n",
    "        new_length = 100\n",
    "        new_indices = np.linspace(0,len(flow_methane)-1,new_length)\n",
    "        spl = UnivariateSpline(old_indices,flow_methane,k=1,s=0)\n",
    "        flow_methane  = np.abs(spl(new_indices))\n",
    "        \n",
    "        old_indices = np.arange(0,len(flow_co2))\n",
    "        new_length = 100\n",
    "        new_indices = np.linspace(0,len(flow_co2)-1,new_length)\n",
    "        spl = UnivariateSpline(old_indices,flow_co2,k=1,s=0)\n",
    "        flow_co2  = np.abs(spl(new_indices))\n",
    "        \n",
    "        old_indices = np.arange(0,len(pressure_out))\n",
    "        new_length = 100\n",
    "        new_indices = np.linspace(0,len(pressure_out)-1,new_length)\n",
    "        spl = UnivariateSpline(old_indices,pressure_out,k=1,s=0)\n",
    "        pressure_out  = np.abs(spl(new_indices))\n",
    "        \n",
    "        old_indices = np.arange(0,len(m.time))\n",
    "        new_length = 100\n",
    "        new_indices = np.linspace(0,len(m.time)-1,new_length)\n",
    "        spl = UnivariateSpline(old_indices,m.time,k=1,s=0)\n",
    "        time  = np.abs(spl(new_indices))\n",
    "        \n",
    "        np_array_rows = time+270,pressure_out,flow_methane,flow_co2\n",
    "        \n",
    "        with open('results12.csv','a') as csvfile:\n",
    "            np.savetxt(csvfile, np_array_rows , delimiter=',',header='P',fmt='%s', comments='')\n",
    "            \n",
    "    met_after_purge= np.empty((seg,len(m.time)))\n",
    "    for i in range(seg):\n",
    "        met_after_purge[i] = np.array(y1[i].value)\n",
    "    met_after_purge= met_after_purge.T\n",
    "\n",
    "    co2_after_purge= np.empty((seg,len(m.time)))\n",
    "    for i in range(seg):\n",
    "        co2_after_purge[i] = np.array(y2[i].value)\n",
    "    co2_after_purge= co2_after_purge.T\n",
    "\n",
    "\n",
    "    n2_after_purge = np.empty((seg,len(m.time)))\n",
    "    for i in range(seg):\n",
    "        n2_after_purge[i] = np.array(y3[i].value)\n",
    "    n2_after_purge= n2_after_purge.T\n",
    "\n",
    "    q1_after_purge = np.empty((seg,len(m.time)))\n",
    "    for i in range(seg):\n",
    "        q1_after_purge[i] = np.array(q1[i].value)\n",
    "    q1_after_purge= q1_after_purge.T\n",
    "\n",
    "    q2_after_purge = np.empty((seg,len(m.time)))\n",
    "    for i in range(seg):\n",
    "        q2_after_purge[i] = np.array(q2[i].value)\n",
    "    q2_after_purge= q2_after_purge.T\n",
    "\n",
    "    pressure_after_purge = np.empty((seg,len(m.time)))\n",
    "    for i in range(seg):\n",
    "        pressure_after_purge[i] = np.array(P[i].value)\n",
    "    pressure_after_purge= pressure_after_purge.T\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3f2af26b",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}