SOchat

# -*- coding: utf-8 -*-
"""
Created on Tue Feb  5 01:42:06 2019

@author: paritosh
"""

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import math as m
#import thermo.mixture as tm  # Module de thermodynamique

# Constante des gaz parfaits
R = 8.314  # J / (mol K)

# Réacteur multitubulaire
# Catalyseur
poro = 0.4  # Porosité
eta_isu = 1  # efficacité
dp = 0.002  # taille des particules en m
# Paramètres de tubes
L = 10  # longueur des tubes en m (imposée)
d = 0.04  # diamètre de tube en m (à faire varier)
NT = 54166  # nombre de tubes (à faire varier)
S = (m.pi * d**2) / 4  # section d'un tube en m²
# Paramètres du réacteur R1
Stot = S * NT  # section totale de réaction en m²
Ftot = 500 * 1000 / 3600  # débit massique totale kg/s
# Densité du mélange à l'entrée du réacteur R1
rho_0 = 12.39  # kg/m3
Ftot_v = Ftot / rho_0  # Débit volumique total m3/s
uv0 = Ftot_v / Stot  # vitesse en fut vide m/s
n = (141.18 - 2) / 14  # Nombre moyen de carbones

# Masses molaires des différents constituants dont l'alcane moyen MW_Cn
MW_CO, MW_H2, MW_N2, MW_Cn, MW_H2O = 0.028, 0.002, 0.028, 0.1412, 0.018  # masses molaires en kg/mol

# Fonction de résoluton du sytème PH1


def FT(y, z):
    # Initialisation du vecteur
    dCOz = y[0]  # fraction massique
    dH2z = y[1]  # fraction massique
    dCNz = y[3]  # fraction massique
    dH2Oz = y[4]  # fraction massique
    Tmel = y[5]  # Température K
    Pmel = y[6]  # Pression en Pa
    uv = y[7]  # vitesse en fût vide en m/s

    # Calcul des fractions molaires de CO et H2 dans la tranche dz
    Wtot = dCOz / MW_CO + dH2z / MW_H2 + 0.414 / MW_N2 + dCNz / MW_Cn + dH2Oz / MW_H2O
    xCO = (dCOz / MW_CO) / Wtot
    xH2 = (dH2z / MW_H2) / Wtot
    xN2 = (0.414 / MW_N2) / Wtot
    xCn = (dCNz / MW_Cn) / Wtot
    xH2O = (dH2Oz / MW_H2O) / Wtot

    # Masse molaire moyenne
    Mm = xCO * MW_CO + xH2 * MW_H2 + xN2 * MW_N2 + xCn * MW_Cn + xH2O * MW_H2O

    # Calcul de rho_melange
    rho_mix = Pmel * Mm / (R * Tmel)

    # Capacité calorifique thermique du mélange gazeux
    # Utilisation des fonctions du module de calcul
    # Réference : https://thermo.readthedocs.io/en/latest/modules.html
    Cpg_mix = 2000 #MODIFIED #tm.Mixture(['carbon monoxide', 'hydrogen', 'nitrogen', 'decane', 'water'], ws=[dCOz, dH2z, 0.414, dCNz, dH2Oz], T=Tmel, P=Pmel).Cpg

    # Calcul de la vitesse
    uv = rho_0 * uv0 / rho_mix

    # Paramètres de la loi cinétique
    a = 126 * m.exp(-100000 / (R * Tmel))
    b = 3.7 * m.exp(-2 * 10000 / (R * Tmel))

    # Calculs des concentrations de CO et H2 dans dz pour le calcul de la cinétique
    CO = Pmel * xCO / (R * Tmel)
    H2 = Pmel * xH2 / (R * Tmel)

    # Loi cinétique
    rFT = ((a * CO * H2) / ((1 + b * CO)**2))  # molCO/(g_Co.s)
    rr = rFT * 950 * 1000  # molCO/(m3_reacteur.s) car masse volumique du cata en kg/m3_reacteur
    r = 0.3 * rr / (1 - poro)  # molCO/(m3_cata_actif.s) car 30% de phase active assimilée au Co et Vreacteur = (1-poro)V_cata

    # Bilan de matière ------------------------------------------------------------------------------------
    # Variation des fractions massiques en fonction de la longueur du réacteur
    dCOz = - eta_isu * (1 - poro) * MW_CO / (rho_0 * uv) * r
    dH2z = - ((2 * n + 1) / n) * eta_isu * (1 - poro) * MW_H2 / (rho_0 * uv) * r
    dCNz = 1 / n * eta_isu * (1 - poro) * MW_Cn / (rho_0 * uv) * r
    dH2Oz = eta_isu * (1 - poro) * MW_H2O / (rho_0 * uv) * r

    # Bilan thermique -------------------------------------------------------------------------------------
    # Voir détails corrélations "Daniel Schweich, Genie de la réaction chimique avancée" à la fin du rapport
    Tp = 201 + 273  # Temperature à la paroi K
    drH = -165 * 1000  # en J / mol_CO

    em = 1  # Emissivité
    lambf = 0.06  # Conductivité thermique du mélange gazeux moyen
    lambs = 0.4  # cf cours Génie de la réaction chimique avancée

    # Calcul du reynolds
    nu = 1.9377 * 10**(-6)  # viscosité cinématique m2 / s
    mu = rho_mix * nu  # viscosité dynamique kg / (m.s)
    uv = (uv0 * rho_0) / rho_mix  # vitesse du gaz en m / s
    Re = rho_mix * uv * dp / mu  # Nombre de Reynolds

    # new correlation
    # Calcul du lamb_sfer
    lamb2 = 0.227 * 10 ** (-6) * em / (2 - em) * dp * Tmel ** 3
    lamb1 = 0.227 * 10 ** (-6) * Tmel ** 3 * dp / (1 + poro / (2 * (1 - poro)) * (1 - em) / em)
    lamb_sf0 = (1 - m.sqrt(1 - poro)) * (1 + poro * (lamb2 / lambf))
    beta = 0.95  # hypothèse
    lamb_sf0 = poro * (1 + beta * lamb1 / lambf)
    Pe_tr = 8.65 * (1 + 19.4 * (dp / d) ** 2)
    Pr = mu * Cpg_mix / lambf
    lamb_ter = Re * Pr / Pe_tr
    lamb_sfer = (lamb_sf0 + lamb_ter) * lambf

    # Calcul du h_intsf
    phi_w = 2.410 * 10 ** (-3) * (d / dp) ** 1.58
    h_ints = (1 - poro) * lambf / ((lambf / (3 * lambs) + phi_w) * dp)
    h_intf0 = 2 * poro * (lambf / dp)
    h_intt = 0.0835 * Re ** 0.91 * (lambf / dp)
    h_intf = h_intf0 + h_intt
    h_intsf = h_ints + h_intf

    # Calcul de h
    h_int = d / (8 * lamb_sfer) + 1 / h_intsf
    h = 1 / h_int
    # fin new correlations

    # Variation de la température en fonction de la longueur du réacteur
    dT = 4 * h / (d * rho_0 * uv0 * Cpg_mix) * (Tp - Tmel) + (1 - poro) / (rho_0 * uv0 * Cpg_mix) * (-drH) * eta_isu * r

    # Bilan de quantité de mouvement ------------------------------------------------------------------------------
    # Bilan sur la pression
    f = ((1 - poro) / poro ** 3) * (4.2 * ((1 - poro) / Re) ** (1 / 6) + 150 * (1 - poro) / Re)
    dP = -f * (uv0 * rho_0) ** 2 / (rho_mix * dp)

    print("%.4f" % uv)
    print("LOOKIE HERE {}".format(uv))
    print("TEST {}".format(dCNz))

    return [dCOz, dH2z, 0.414, dCNz, dH2Oz, dT, dP, uv]


# Résolution ------------------------------------------------------------------------------------------------------
C0 = [0.516, 0.070, 0.414, 0, 0, 498, 35 * 10**5, uv0]
z = np.linspace(0, L, 100)
y = odeint(FT, C0, z)

# Graphiques -------------------------------------------------------------------------------------------------------
CO = y[:, 0]
H2 = y[:, 1]
N2 = y[:, 2]
CN = y[:, 3]
H2O = y[:, 4]

T = y[:, 5]
P = y[:, 6]
uv = y[:, 7]
uv = np.diff(uv) #MODIFIED

b = len(z)
print('rho_0 = ', " %.2f" % rho_0, 'kg/m3')
print('uv0 = ', "%.2f" % uv0, 'm/s')
print('Les fractions massiques finales sont : wCO_s =' " %.3f" % y[b - 1, 0], ', wH2_s =' " %.3f" % y[b - 1, 1], ', wCN_s =' " %.3f" % y[b - 1, 3], ', wH2O_s =' " %.3f" % y[b - 1, 4])
print('La conversion en CO est :', " %.2f" % (100 * (1 - y[b - 1, 0] / y[0, 0])), '%')
print('La température finale est ', " %.1f" % y[b - 1, 5], 'K')
print('La pression finale est ', " %.0f" % y[b - 1, 6], 'Pa')
print('DeltaT_max = ', " %.0f" % abs(max(T) - T[0]), 'K')


# Graphiques
# Il faut décommenter les graphiques pour les voir, par défaut seul les fractions massiques sont visibles

# --------------------------------------------------------------------
fig, ax = plt.subplots()  # Ne pas mettre cette ligne en commentaire

# ---------------------------------------------------------------------
# Graphique vitesse en fût vide
ax.plot(z[:-1], uv, 'r:', label='uv') #MODIFIED
legend = ax.legend(loc='upper right', shadow=True, fontsize='x-large')
legend.get_frame().set_facecolor('#00FFCC')
plt.xlabel('Longueur du réacteur en m')
plt.ylabel('Vitesse en fût vide en m/s')
plt.title('Vitesse en fût vide en fonction de la longueur du réacteur R1')
plt.show()