Pseudo-homogenes Modell
[1] BUSSCHE, KM Vanden; FROMENT, G. F. A Steady-State Kinetic Model for Methanol Synthesis and the Water Gas Shift Reaction on a Commercial Cu/ZnO/Al2O3Catalyst. Journal of Catalysis, 1996, 161. Jg., Nr. 1, S. 1-10.
[2] ASKGAARD, T. S., et al. A kinetic model of methanol synthesis. Journal of Catalysis, 1995, 156. Jg., Nr. 2, S. 229-242.
[3] LØVIK, Ingvild. Modelling, estimation and optimization of the methanol synthesis with catalyst deactivation. 2001.
[4] Froment, Gilbert F. ; Bischoff, Kenneth B. ; Wilde, Juray De: Chemical Reactor Analysis and Design. New York: Wiley, 2010.
[5] Bird, R. Byron ; Stewart, Warren E. ; Lightfoot, Edwin N.: Transport Phenomena. New York: John Wiley & Sons, 2007.
[6] e.V., VDI: VDI-Wärmeatlas. Wiesbaden: Springer Berlin Heidelberg, 2013.
[7] Gmehling, Jürgen ; Kolbe, Bärbel ; Kleiber, Michael ; Rarey, Jürgen: Chemical Thermodynamics for Process Simulation. New York: John Wiley & Sons, 2012.
[8] Ott, J. , Gronemann, V. , Pontzen, F. , Fiedler, E. , Grossmann, G. , Kersebohm, D. B., Weiss, G. and Witte, C. (2012). Methanol. In Ullmann's Encyclopedia of Industrial Chemistry, (Ed.).
[9] Baerns, Manfred ; Behr, Arno ; Brehm, Axel ; Gmehling, Jürgen ; Hinrichsen, Kai-Olaf ; Hofmann, Hanns ; Onken, Ulfert ; Palkovits, Regina ; Renken, Albert: Technische Chemie. New York: John Wiley & Sons, 2014.
[10] Fogler, H. Scott: Essentials of Chemical Reaction Engineering. Amsterdam: Pearson Education, 2011.
[11] Smith, Joseph Mauk ; Ness, Hendrick C. Van ; Abbott, Michael M.: Introduction to chemical engineering thermodynamics. New York: McGraw-Hill, 2005.
$\begin{array}{lcccc} \text{1. Hydrierung von CO2:}& CO_2 &+ 3 H_2 &\rightleftharpoons CH_3OH & +H_2O \\ \text{2. Hydrierung von CO:}& CO &+ 2 H_2 &\rightleftharpoons CH_3OH & \\ \text{3. Reverse-Wassergas-Shift-Reaktion:}& CO_2 &+ H_2 &\rightleftharpoons H_2O & +CO \\ \end{array}$
$\begin{array}{ll} log_{10}(K_1) = \frac{3066K}{T}-10,592\\ log_{10}(K_2) = \frac{5139K}{T}-12,621\\ log_{10}(1/K_3) = \frac{-2073K}{T}+2,029\\ \end{array}$
Zunächst nach Nomenklatur aus Fogler CRT
Stoffmengenbilanzen
Akkumulation und Reaktion (Quelle-Senke):
$\frac{d \dot n_{i}}{d m_{Kat}}=\sum\limits_j^{\mathcal R}{\nu_{ji} r_j}$
$ \frac{d \dot n_{MeOH}}{d m_{Kat}}=+r_{MeOH}\\ \frac{d \dot n_{CO_2}}{d m_{Kat}}=-r_{MeOH}-r_{RWGS}\\ \frac{d \dot n_{CO}}{d m_{Kat}}=+r_{RWGS}\\ \frac{d \dot n_{H_2O}}{d m_{Kat}}=+r_{MeOH}+r_{RWGS}\\ \frac{d \dot n_{H_2}}{d m_{Kat}}=-3r_{MeOH}-r_{RWGS} $
Energiebilanz
$\sum\limits_j^{\mathcal K}\dot n_j C_{p j} \frac{d T}{d m_{Kat}} = \sum\limits_i^{\mathcal R}(-\Delta H_i) r_i-\underbrace{4\frac{U}{d_t}\cdot \frac{1}{\rho_B}(T-T_r)}_{\begin{array}{l}U a= U (A/V) = U \cdot\frac{\pi d_t L_R}{{(\pi/4) d_t^2 L_R}}\\ \frac{Ua}{\rho_B}[=]\frac{W}{kgKat K}\end{array}}$
Kontinuität
$\rho \dot V = \rho_0 \dot V_0$
Druckverlust
Ergun Widerstandsbeiwertgleichung
$\frac{d P}{d m_{Kat}} = -\frac{\alpha}{2}\frac{P_0}{P/P_0}\frac{T}{T_0}\frac{\dot n_T}{\dot n_{T0}}$
$\alpha=\frac{2 \beta_0}{A_c(1-\phi)\rho_c\cdot P_0}$
$\beta_0=\frac{G}{\rho_0 g_c D_p}\cdot\left(\frac{1-\phi}{\phi^3}\right)\left[\frac{150(1-\phi)\mu}{D_p}+1,75G\right]$
Kinetik-Ausdrücke, aus Van den Bussche-Froment, auf $Cu/ZnO/Al_2O_3$ Katalysator
$p_i$ hier eigentlich: $p_i/p^{\circ}$ , $p^{\circ}=1bar$. Es heißt, $r_{MeOH}$ hat die Einheiten $r_{MeOH}[=]k_{5a}'=k_{5a}c_t^2=\frac{mol}{kg_{Kat}\cdot s}$, und die Reaktionsquotienten haben keine Einheiten.
$r_{MeOH}=\frac{k_{5a}' K_2' K_3 K_4 K_{H_2} p_{CO_2}p_{H_2}\left(1-\frac{1}{K_1^*}\cdot\frac{p_{H_2O}\cdot p_{MeOH}}{p_{H_2}^3\cdot p_{CO_2}}\right)}{\left(1+\left(\frac{K_{H_2O}}{K_8 K_9 K_{H_2}}\right)\frac{p_{H_2O}}{p_{H_2}}+\sqrt{K_{H_2}\cdot p_{H2}}+K_{H_2O}\cdot p_{H_2O}\right)^3}$
$r_{RWGS}=\frac{k_1' p_{CO_2}\left(1-K_3^*\cdot\frac{p_{H_2O}\cdot p_{CO}}{p_{H_2}\cdot p_{CO_2}}\right)}{\left(1+\left(\frac{K_{H_2O}}{K_8 K_9 K_{H_2}}\right)\frac{p_{H_2O}}{p_{H_2}}+\sqrt{K_{H_2}\cdot p_{H2}}+K_{H_2O}\cdot p_{H_2O}\right)}$
Konstitutive Gleichungen
$p_i = y_i P$
$y_i = \frac{\dot n_i}{\dot n_T}$
$\dot n_T = c_T\dot V = c_T \cdot u_s \cdot A_c$
$c_T = \frac{P}{R T}$
Zwangsbedingungen
$\dot n_T = \sum\limits_i^{\mathcal K}{\dot n_i}$
| Symbol | Beschreibung | Einheiten [4] | Einheiteh (SI) |
|---|---|---|---|
| $M_m$ | Mittlere Molmasse $M_m=\sum\limits_i^{\mathcal K}{y_i M_i}$ | $\frac{g}{mol}$ | $\frac{g}{mol}$ |
| $c_T$ | Molare Konzentration $c_T = \frac{\dot n}{\dot V}$ | $\frac{mol}{m^3}$ | $\frac{mol}{m^3}$ |
| $m_{Kat}$ | Masse Katalysator | $kg_{Kat}$ | $kg_{Kat}$ |
| $\phi$ | Hohlraum-Verhältnis (Void fraction) | $\frac{m^3_{Hohl}}{m^3_{Schüttung}}$ | $\frac{m^3_{Hohl}}{m^3_{Schüttung}}$ |
| $\rho_c$ | Katalysator Feststoffdichte | $\frac{kg Kat}{m^3_{Feststoff}}$ | $\frac{kg Kat}{m^3_{Feststoff}}$ |
| $\rho_B$ | Katalysator Schüttdichte $\rho_B=(1-\phi)\rho_c$ | $\frac{kg Kat}{m^3_{Schüttung}}$ | $\frac{kg Kat}{m^3_{Schüttung}}$ |
| $\rho_g$ | Gas-Dichte $\rho_g = C_T\cdot M_m$ | $\frac{kg}{m^3_{Gas}}$ | $\frac{kg}{m^3_{Gas}}$ |
| $d_t$ | Innerer Rohrdurchmesser (auch $D_R$) | m | cm |
| $R_r$ | Innerer Rohrradius | m | cm |
| $d_p$ | Partikeldurchmesser (auch $D_P$) | m | mm |
| $A_c$ | Querschnittsfläche $A_c=A_{Quer}=(\pi/4) d_t^2$ | $m^2$ | $m^2$ |
| $f$ | Widerstandsbeiwert: Friction factor (auchn $\xi$) | Dimensionslos | Dimensionslos |
| $\dot m$ | Massenstrom | $\frac{kg}{s}$ | $\frac{kg}{h}$ |
| $\dot n$ | Stoffmengenstrom | $\frac{mol}{s}$ | $\frac{kmol}{h}$ |
| $\dot V$ | Volumenstrom | $\frac{m^3}{h}$ | $\frac{m^3}{s}$ |
| $u_g$ | Lineare Gasgeschwindigkeit [3]: $u_g=\frac{\dot V}{A_{Quer}}$ | $\frac{m_{Schüttung}}{h}$ | $\frac{m_{Schüttung}}{s}$ |
| $u_s$ | Oberflächenspezifische Geschwindigkeit [4]: $u_s=\frac{\dot V}{A_{Quer}}$ | $\frac{m_{Gas}^3}{m_{Schüttung}^2\cdot s}$ | $\frac{m_{Gas}^3}{m_{Schüttung}^2\cdot s}$ |
| $u$ | Oberflächenspezifische Geschwindigkeit [10] $u_s=\frac{\dot V}{A_{Quer}}$ | $\frac{ft}{h}$ | $\frac{m}{s}$ |
| $G$ | Oberflächenspezifischer Massenstrom $G=\frac{\dot m}{A_c}=\rho_g u_s$ | $\frac{kg}{m^2 s}$ | $\frac{kg}{m^2 s}$ |
| $M'$ | Massenstromdichte (Massen-Flux) $M'=G$ | $\frac{kg}{h \cdot m^2}$ | $\frac{kg}{s \cdot m^2}$ |
| $\dot n'$ | Stoffmengenstromdichte (Stoffmengen-Flux) $\dot n'=\frac{\dot n}{A_c}$ | $\frac{kmol}{h \cdot m^2}$ | $\frac{kmol}{h \cdot m^2}$ |
| $C_{p}$ | Stoffmengenbezogene Wärmekapazität | $\frac{kcal}{kmol K}$ | $\frac{J}{mol K}$ |
| $C_{p_g}$ | Massenbezogene Wärmekapazität $C_{p_g}=\frac{\sum\limits_j^{\mathcal K} y_j\cdot C_{p j}}{\sum\limits_j^{\mathcal K }y_j \cdot M_j} $ | $\frac{kcal}{kg K}$ | $\frac{J}{kg K}$ |
| $M'C_{p_g}$ | Konvektive Wärmestromdichte | $\frac{kcal}{h \cdot m^2 \cdot K}$ | $\frac{J}{s\cdot m \cdot K}$ |
| $\dot n C_{p}$ | Konvektive Wärmestromdichte | $\frac{kcal}{h \cdot m^2 \cdot K}$ | $\frac{J}{s\cdot m \cdot K}$ |
| $Re$ | Reynoldszahl $Re=\frac{D_p u_g \rho_g}{\mu}=\frac{G \cdot D_P}{\mu}$ | (Dimensionslos) | (Dimensionslos) |
$m_{Kat}=\underbrace{(1-\phi)A_c\cdot z}_{\text{Feststoff-Volumen}}\cdot \underbrace{\rho_c}_{\text{Stoffdichte Kat.}} = (1-\phi) \cdot \rho_c \cdot z \cdot A_c =\rho_B \cdot z \cdot A_c $
Querschnittsfläche
$A_c = (\pi/4)d_t^2$
Aus den konstitutiven Gleichungen (Molanteil)
$\frac{d \dot n_i}{d m_{Kat}}=\frac{d y_i\dot n_T}{d m_{Kat}}=y_i\frac{d \dot n_T}{d m_{Kat}}+\dot n_T\frac{d y_i}{d m_{Kat}}$
$\dot n_T = c_T\dot V = c_T\cdot u_g \cdot A_c = \frac{G \cdot A_c}{M_m}$
Reaktorposition entdimensionieren, indem man durch Reaktorlänge teilt:
$z=z_{absolut}/L_R$
$\Rightarrow\begin{array}{ll} \frac{c_T\cdot u_g \cdot A_c}{\rho_B \cdot A_c \cdot L_R}\frac{d y_i}{d (z/L_r)}&=\frac{d \dot n_i}{d m_{Kat}}-y_i\frac{d \dot n_T}{d m_{Kat}}\\ &=\sum\limits_j^{\mathcal R}{\nu_{ji} r_j}-y_i\sum\limits_j^{\mathcal R}{{\sum\limits_i^{\mathcal K}{\nu_{ij}r_j}}} \end{array}$
$\begin{array}{lrl} \frac{d P}{d m_{Kat}}&=\frac{1}{\rho_B \cdot A_c \cdot L_R}\frac{d P}{d (z/L_r)}&=-\frac{\alpha}{2}\frac{P_0}{P/P_0}\frac{T}{T_0}\frac{\dot n_T}{\dot n_{T0}}=-\frac{2\beta_0}{2 A_c \underbrace{(1-\phi)\rho_c}_{\rho_B}\cdot P_0}\frac{P_0}{P/P_0}\frac{T}{T_0}\frac{\dot n_T}{\dot n_{T0}}\\ &&=-\frac{\beta_0}{A_c \rho_B} \frac{c_{T 0}}{c_T}\cdot\frac{\dot n_T}{\dot n_{T0}}=-\frac{\beta_0}{A_c \rho_B} \frac{\dot V}{\dot V_0}=-\frac{\beta_0}{A_c \rho_B} \frac{\rho_0 \dot V_0}{\rho \dot V_0}\\ \Rightarrow&\frac{d P}{d(z/L_R)}&=-L_R\cdot \frac{G}{\rho g_c D_p}\cdot\left(\frac{1-\phi}{\phi^3}\right)\left[\frac{150(1-\phi)\mu}{D_p}+1,75G\right] \end{array}$
$\rho_g = c_T M_m = c_T \cdot \sum\limits_i^{\mathcal K}{y_i M_i}$
$\sum\limits_j^{\mathcal K}\dot n_j C_{p j} \frac{d T}{d m_{Kat}}=\frac{\sum\limits_j^{\mathcal K}y_j\cdot \dot n_T C_{p j} }{\rho_B A_c}\cdot \frac{d T}{d z}=\frac{\frac{G\cdot A_c}{\sum\limits_i^{\mathcal K }y_i M_i}\cdot\sum\limits_j^{\mathcal K} y_j\cdot C_{p j} }{\rho_B \cdot A_c}\cdot \frac{d T}{d z}=\frac{G \cdot C_{p_g} }{\rho_B}\cdot \frac{d T}{d z}$
Zusammenfassung der Erhaltungsgleichungen als Funktion der Länge, Molanteile, Temperatur und Druck:
$\begin{array}{ll} \frac{d y_i}{d (z/L_R)}&=L_R \cdot \frac{\rho_B\left(\sum\limits_i^{\mathcal K}{y_i M_i}\right)}{G}\cdot\left(\sum\limits_j^{\mathcal R}{\nu_{ji} r_j}-y_i\sum\limits_j^{\mathcal R}{{\sum\limits_i^{\mathcal K}{\nu_{ij}r_j}}}\right)\\ \frac{d P}{d (z/L_R)}&= -L_R\cdot \frac{G}{\left(c_T \cdot \sum\limits_i^{\mathcal K}{y_i M_i}\right) g_c D_p}\cdot\left(\frac{1-\phi}{\phi^3}\right)\left[\frac{150(1-\phi)\mu}{D_p}+1,75 G\right]\\ \frac{dT}{d (z/L_R)}&=L_R\cdot\frac{1}{G \cdot C_{p_g}}\cdot\left(\frac{2 U}{R_r}(T_r-T)+\rho_B\sum\limits_i^{\mathcal R}(-\Delta H_i) r_i\right) \end{array}$
Alternative Form (Form nach [3]):
$\begin{array}{ll} \frac{d y_i}{d (z/L_R)}&=L_R \cdot \frac{\rho_B}{c_T \cdot u_g}\cdot\left(\sum\limits_j^{\mathcal R}{\nu_{ji} r_j}-y_i\sum\limits_j^{\mathcal R}{{\sum\limits_i^{\mathcal K}{\nu_{ij}r_j}}}\right)\\ \frac{d P}{d (z/L_R)}&= -L_R\cdot\frac{10^{-5}bar}{Pa}\cdot \frac{u_g^2 \left(c_T \cdot \sum\limits_i^{\mathcal K}{y_i M_i}\right)}{ D_p}\cdot\left(\frac{1-\phi}{\phi^3}\right)\left[1,75+150\frac{(1-\phi)}{Re}\right]\\ \frac{dT}{d (z/L_R)}&=L_R\cdot\frac{1}{G \cdot C_{p_g}}\cdot\left(\frac{2 U}{R_r}(T_r-T)+\rho_B\sum\limits_i^{\mathcal R}(-\Delta H_i) r_i\right) \end{array}$
1. Anzahl an Übertragungseinheiten (VDI-WA C1-2 [6])
$NTU_i \equiv \frac{k A}{\dot M_i\cdot Cp_{mi}}$
$NTU = L_R \cdot \frac{1}{G\cdot C_{p_g}}\cdot \frac{2 U}{R_R}$
...denn
$\begin{array}{ll} G&=\frac{\dot m}{A_{Quer}}\\ U\cdot a&=U\cdot \frac{A_{Übertragung}}{V_R}=U\cdot \frac{2\cdot\pi R_R L_R}{\pi R_R^2 L_R}=\frac{2}{R_R}\cdot U\\ A_{Quer} &= \pi\cdot R_R^2\\ \Rightarrow NTU=& L_R \cdot \frac{1}{G\cdot C_{p_g}}\cdot \frac{2 U}{R_R}=L_R \cdot \frac{\pi\cdot R_R^2}{\dot m\cdot C_{p_g}}\cdot \frac{2 U}{R_R}\\ &=\frac{(2\cdot\pi\cdot R_R\cdot L_R) U}{\dot m\cdot C_{p_g}}=\frac{U A_{Wärmeübtragung}}{\dot m\cdot C_{p_g}} \quad \checkmark \end{array}$
Systeme, worin gleiche NTU's überleitet werden, weisen ein gleiches Verhalten auf, und die NTU's lassen sich hauptsächlich durch 5 Größen beschreiben:
$NTU=\frac{(2\cdot\pi\cdot R_R\cdot L_R) U}{\dot m\cdot C_{p_g}}$
Hier ist der Massenstrom auf den einzelnen Rohr bezogen. D. h. Gesammtmassenstrom $\dot m_{ges}=n_T\cdot \dot m$. In diesem Sinn ist die Anzahl an Rohren $n_T$ auch ein Parameter.
2. Damköhler-Zahl (III) [9]
$Da_{III}=k \tau\left(\frac{c\cdot \Delta H_R}{C_{p_g} \cdot \rho \cdot T_0}\right)=\frac{Wärmeaus der Reaktion}{Wärmekonvektion}$
3. Reaktionswärme / Wärmeübertragung
Beim isothermen System (oder bei Temperatur-Maxima) reduziert sich die Wärmebilanz, als Funktion der NTU-Anzahl, zu:
$0=NTU\cdot(Tr-T)+NTU\cdot\frac{\rho_b \sum\limits_i^{\mathcal R}(-\Delta H_i) r_i}{U a}$
Das daraus resultierende Verhältnis wird als $Q_{RW}$ bezeichnet, es nimmt bei Temperatur-Maxima und bei isothermen Prozessen den Wert 1 an.
$Q_{RW} = \frac{\rho_b \sum\limits_i^{\mathcal R}(-\Delta H_i) r_i}{U a (T-T_r)}=\frac{\rho_b \sum\limits_i^{\mathcal R}(-\Delta H_i) r_i}{\frac{2\cdot U}{R_R} (T-T_r)}$
In [1]:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import lines
from scipy.integrate import odeint
%matplotlib inline
# Journal of Catalysis, 1996, 161. Jg., Nr. 1, S. 1-10.
# Bezugdaten: 'Løvik, Ingvild. "Modelling, estimation and optimization
# of the methanol synthesis with catalyst deactivation." (2001).
# Katalysator
rho_b = 1775 # kg Kat/m^3 Feststoff
phi = 0.5 # m^3 Gas/m^3 Feststoff
m_kat = 34.8 / 1000. # kg Kat
d_p = 0.0005 # m Feststoff
# Reaktor
d_t = 0.016 # m Rohrdurchmesser
l_r = 0.15 # m Rohrlänge
# Betriebsbedingungen
t0 = 493.2 # K
p0 = 50. # bar
m_dot = 2.8 * 1e-5 # kg/s
n_t = 1
# Zulaufbedingungen
namen = ['CO', 'H2O', 'MeOH', 'H2', 'CO2', 'N2']
y_i0 = np.array([
4, 0, 0, 82, 3, 11
], dtype=float) / 100.
mm = np.array([
28, 18, 32, 2, 44, 28
], dtype=float) # g/mol
# Wärmetausch
t_r = 25 + 273.15 # K
u = 0 # Adiabat
ntu = 0
# Stoechiometrische Koeffizienten
nuij = np.zeros([len(namen), 3])
# Hydrierung von CO2
nuij[[
namen.index('CO2'),
namen.index('H2'),
namen.index('MeOH'),
namen.index('H2O'),
namen.index('CO'),
namen.index('N2'),
], 0] = np.array([-1, -3, +1, +1, 0, 0], dtype=float)
# Hydrierung von CO
nuij[[
namen.index('CO2'),
namen.index('H2'),
namen.index('MeOH'),
namen.index('H2O'),
namen.index('CO'),
namen.index('N2'),
], 1] = np.array([0, -2, +1, 0, -1, 0], dtype=float)
# RWGS Reverse-Wassergasshiftreaktion (muss gleich sein als die Vorwärtsreaktion,
# wofür die Kinetik verfügbar ist)
nuij[[
namen.index('CO2'),
namen.index('H2'),
namen.index('MeOH'),
namen.index('H2O'),
namen.index('CO'),
namen.index('N2'),
], 2] = - np.array([+1, +1, 0, -1, -1, 0], dtype=float)
# Quelle: Chemical thermodynamics for process simulation
h_298 = np.array([
-110540, -241820, -201160, 0, -393500, 0
], dtype=float) # J/mol
delta_h_r_298 = nuij.T.dot(h_298) # J/mol
# Parameter der Gleichung aus dem VDI-Wärmeatlas
cp_konstanten = np.zeros([len(namen), 7])
cp_konstanten_text = [
'407,9796 3,5028 2,8524 –2,3018 32,9055 –100,1815 106,1141',
'706,3032 5,1703 –6,0865 –6,6011 36,2723 –63,0965 46,2085',
'846,6321 5,7309 –4,8842 –12,8501 78,9997 –127,3725 82,7107',
'392,8422 2,4906 –3,6262 –1,9624 35,6197 –81,3691 62,6668',
'514,5073 3,4923 –0,9306 –6,0861 54,1586 –97,5157 70,9687',
'432,2027 3,5160 2,8021 –4,1924 42,0153 –114,2500 111,1019',
]
for i in range(len(namen)):
cp_konstanten[i, :] = np.array(
cp_konstanten_text[i].replace(
'–', '-').replace(',', '.').split(' '),
dtype=float)
# Lennard-Jones Parameter Prausnitz Anhang B
# epsilon / kB
l_j_epsilon_d_k = np.array([
91.7, 809.1, 481.8, 59.7, 195.2, 71.4
]) # K
# sigma
l_j_sigma = np.array([
3.690, 2.641, 3.626, 2.827, 3.941, 3.798
]) # Angstrom
def cp_ig_durch_r(t):
a, b, c, d, e, f, g = cp_konstanten.T
gamma_var = t / (a + t)
return b + (c - b) * gamma_var ** 2 * (
1 + (gamma_var - 1) * (
d + e * gamma_var + f * gamma_var ** 2 + g * gamma_var ** 3
)) # dimensionslos
def int_cp_durch_r_dt_minus_const(t):
a, b, c, d, e, f, g = cp_konstanten.T
return b * t + (c - b) * (
t - (d + e + f + g + 2) * a * np.log(a + t) +
-(2 * d + 3 * e + 4 * f + 5 * g + 1) * a**2 / (a + t) +
+(1 * d + 3 * e + 6 * f + 10 * g) * a**3 / 2 / (a + t)**2 +
-(1 * e + 4 * f + 10 * g) * a**4 / 3 / (a + t)**3 +
+(1 * f + 5 * g) * a**5 / 4 / (a + t)**4 +
- g * a**6 / 5 / (a + t)**5
)
def int_cp_durch_rt_dt_minus_const(t):
a, b, c, d, e, f, g = cp_konstanten.T
return b * np.log(t) + (c - b) * (
np.log(a + t) + (1 + d + e + f + g) * a / (a + t) +
-(d / 2 + e + 3 * f / 2 + 2 * g) * a**2 / (a + t)**2 +
+(e / 3 + f + 2 * g) * a**3 / (a + t)**3 +
-(f / 4 + g) * a**4 / (a + t)**4 +
+(g / 5) * a**5 / (a + t)**5
)
def mcph(x, t0_ref, t):
return sum(x * (
int_cp_durch_r_dt_minus_const(t) -
int_cp_durch_r_dt_minus_const(t0_ref))
) / sum(x) / (t - t0_ref)
def delta_h_r(t):
cp_m = (
int_cp_durch_r_dt_minus_const(t) -
int_cp_durch_r_dt_minus_const(298.15)
) / (t - 298.15) * 8.3145 # J/mol/K * K = J/mol
return delta_h_r_298 + nuij.T.dot(cp_m) * (t - 298.15)
def mu(t, y_i):
# T* = k T / epsilon
t_st = t / l_j_epsilon_d_k # J/K
# Stoßintegral (Bird Tabelle E.2)
omega_mu = 1.16145 / t_st**0.14874 + \
0.52487 / np.exp(0.77320 * t_st) + \
2.16178 / np.exp(2.43787 * t_st)
konst_1 = 5 / 16. * np.sqrt(
8.3145 * 1000 * 100 ** 2 / np.pi
) * 10**16 / 6.022e23 # g/cm/s
mu_i = konst_1 * np.sqrt(mm * t) / (
l_j_sigma**2 * omega_mu
) * 100 / 1000
# g/cm/s * 100cm/m * 1kg/1000g = kg/m/s = Pa s
phi_ab = 1 / np.sqrt(8) * (
1 + np.outer(mm, 1 / mm)) ** (-1 / 2.) * (
1 + np.outer(mu_i, 1 / mu_i) ** (1 / 2.) *
np.outer(1 / mm, mm) ** (1 / 4.)
) ** 2
mu_mix = sum(y_i * mu_i / phi_ab.dot(y_i))
return mu_mix # Pa s
# T-abhängige Parameter, dem Artikel nach
def k_t(t):
r = 8.3145 # Pa m^3/mol/K
# Geschwindigkeitskonstanten der 3 Reaktionen
k_1 = 10 ** (3066 / t - 10.592)
k_2 = 10 ** (5139 / t - 12.621)
k_3 = 10 ** (2073 / t - 2.029)
# Angepasste Parameter des kinetischen Modells
# A(i) exp(B(i)/RT)
a = np.array([
0.499, 6.62e-11, 3453.38, 1.07, 1.22e10
], dtype=float)
b = np.array([
17197, 124119, 0, 36696, -94765
], dtype=float)
k_h2 = a[0] ** 2 * np.exp(2 * b[0] / (r * t))
k_h2o = a[1] * np.exp(b[1] / (r * t))
k_h2o_d_k_8_k_9_k_h2 = a[2] * np.exp(b[2] / (r * t))
k5a_k_2_k_3_k_4_k_h2 = a[3] * np.exp(b[3] / (r * t))
k1_strich = a[4] * np.exp(b[4] / (r * t))
return np.array([
k_1, k_2, k_3,
k_h2, k_h2o, k_h2o_d_k_8_k_9_k_h2,
k5a_k_2_k_3_k_4_k_h2, k1_strich
])
def r_i(t, p_i):
p_co2 = p_i[namen.index('CO2')]
p_co = p_i[namen.index('CO')]
p_h2 = p_i[namen.index('H2')]
p_meoh = p_i[namen.index('MeOH')]
p_h2o = p_i[namen.index('H2O')]
[k_1, _, k_3,
k_h2, k_h2o, k_h2o_d_k_8_k_9_k_h2,
k5a_k_2_k_3_k_4_k_h2, k1_strich] = k_t(t)
r_meoh = k5a_k_2_k_3_k_4_k_h2 * p_co2 * p_h2 * (
1 - 1 / k_1 * p_h2o * p_meoh / (
p_h2 ** 3 * p_co2
)
) / (
1 + k_h2o_d_k_8_k_9_k_h2 * p_h2o / p_h2 +
np.sqrt(k_h2 * p_h2) + k_h2o * p_h2o
) ** 3
r_rwgs = k1_strich * p_co2 * (
1 - k_3 * p_h2o * p_co / (p_co2 * p_h2)
) / (
1 + k_h2o_d_k_8_k_9_k_h2 * p_h2o / p_h2 +
np.sqrt(k_h2 * p_h2) + k_h2o * p_h2o
) ** 1
return np.array([r_meoh, 0., r_rwgs])
def df_dt(y, _):
# FIXME: Normalize y_i on each iteration.
# Inerts should not change at all, but there is a 0.00056% increase already
# on the first time step. Test:
# g*60**2*n_t*(np.pi / 4 * d_t**2)/sum(y_i*mm/1000.)*y_i[6]*mm[6]/1000.
y_i = y[:-2]
p = y[-2]
t = y[-1]
mm_m = sum(y_i * mm) * 1 / 1000. # kg/mol
cp_m = sum(y_i * cp_ig_durch_r(t) * 8.3145) # J/mol/K
cp_g = cp_m / mm_m # J/kg/K
c_t = p / (8.3145 * 1e-5 * t) # bar * mol/bar/m^3/K*K = mol/m^3
p_i = y_i * p # bar
delta_h_r_t = delta_h_r(t)
mu_t_y = mu(t, y_i)
r_j = r_i(t, p_i) # mol/kg Kat/s
dyi_dz = l_r * rho_b * mm_m / g * (
nuij.dot(r_j) - y_i * sum(nuij.dot(r_j))
) # m * kgKat/m^3 * kg/mol * m^2 s/kg * mol/kgKat/s = dimlos
dp_dz = -l_r * 1 / 10**5 * g / (
c_t * mm_m * d_p
) * (1 - phi) / phi**3 * (
150 * (1 - phi) * mu_t_y / d_p + 1.75 * g
) # Pa * 1bar/(1e5 Pa) = bar
dt_dz = l_r / (
g * cp_g
) * (
2 * u / (d_t / 2) * (t_r - t) + rho_b * -delta_h_r_t.dot(r_j)
)
result = np.empty_like(y)
result[:-2] = dyi_dz
result[-2] = dp_dz
result[-1] = dt_dz
return result
# Berechnung der Parameter
g = m_dot / (np.pi / 4 * d_t**2) / n_t # kg/m^2/s
mm_m_0 = sum(y_i0 * mm) * 1 / 1000. # kg/mol
cp_m_0 = sum(y_i0 * cp_ig_durch_r(t0) * 8.3145) # J/mol/K
cp_g_0 = cp_m_0 / mm_m_0 # J/kg/K
sn = (y_i0[namen.index('H2')] - y_i0[namen.index('CO2')]) / (
y_i0[namen.index('CO2')] + y_i0[namen.index('CO')]
)
y_0 = np.empty([len(y_i0) + 1 + 1])
y_0[:-2] = y_i0
y_0[-2] = p0
y_0[-1] = t0
z_d_l_r = np.linspace(0, 1, 200)
soln = odeint(df_dt, y_0, z_d_l_r)
y_i_soln = soln[:, :len(y_i0)]
p_soln = soln[:, -2]
t_soln = soln[:, -1]
vars_1 = [
[r'\rho_b', rho_b, r'\frac{kg_{Kat}}{m^3_{Schüttung}}'],
['\phi', phi, r'\frac{m^3_{Gas}}{m^3_{Schüttung}}'],
['D_p', d_p, 'm'],
['D_t', d_t, 'm'],
['L_R', l_r, 'm'],
]
vars_2 = [
['n_T', n_t, ''],
['U', u, r'\frac{W}{m^2\cdot K}'],
['\dot m', m_dot, 'kg/s'],
['C_{p_g}', cp_g_0 / 1000., r'\frac{kJ}{kg\cdot K}'],
['NTU', ntu, ''],
]
vars_3 = [
['T_0', t0 - 273.15, '°C'],
['P_0', p0, 'bar'],
# ['T_r', t_r-273.15, '°C_{Kühlmittel}'],
# ['P_{Sät}', p_sat, 'bar_{Kühlmittel}'],
# ['\Delta H_{Sät}', h_sat_l, r'\frac{kJ}{kg_{Kühlmittel}}'],
['SN', sn, '']
]
text_1 = '\n'.join(['$' + ' = '.join([line[0], '{:g}'.format(line[1]) +
' ' + line[2]]) + '$'
for line in vars_1])
text_2 = '\n'.join(['$' + ' = '.join([line[0], '{:g}'.format(line[1]) +
' ' + line[2]]) + '$'
for line in vars_2])
text_3 = '\n'.join(['$' + ' = '.join([line[0], '{:g}'.format(line[1]) +
' ' + line[2]]) + '$'
for line in vars_3])
t_r = np.linspace(100, 1000, 200)
fig = plt.figure(1, figsize=(8, 8))
fig.suptitle('Adiabates System (Journal of Catalysis, ' +
'1996, 161. Jg., Nr. 1, S. 1-10. )')
fig.text(0.05, 0.945, text_1, va='top', fontsize=10)
fig.text(0.33, 0.935, text_2, va='top', fontsize=10)
fig.text(0.66, 0.930, text_3, va='top', fontsize=10)
ax = plt.subplot2grid([2, 3], [0, 0])
ax.plot(t_r, np.array(
[mu(t, y_i0) for t in t_r]) / 1e-5, label='Berechnet')
ax.plot(t_r, np.array(
[67.2e-7 + 0.21875e-7 * t for t in t_r]) / 1e-5, label='Bezugdaten')
ax.set_ylabel(r'$\frac{\mu}{(Pa s) \cdot 1e-5}$')
plt.setp(ax.get_xticklabels(), visible=False)
ax.legend(fontsize='medium')
ax2 = plt.subplot2grid([2, 3], [1, 0], sharex=ax)
delta_h_r_t_0 = np.array([delta_h_r(t) for t in t_r])
ax2.plot(t_r, delta_h_r_t_0[:, 0] / 1000.,
label='$Hyd. CO_2 - berechnet$')
ax2.plot(t_r, delta_h_r_t_0[:, 2] / 1000.,
label='$RWGS - berechnet$')
ax2.plot(t_r, -(57980 + 35 * (t_r - 498.15)) / 1000.,
label='$Hyd. CO_2 - Bezugsdaten$')
ax2.plot(t_r, -(-39892 + 8 * (t_r - 498.15)) / 1000.,
label='$RWGS - Bezugsdaten$')
ax2.set_xlabel('T/K')
ax2.set_ylabel(r'$\frac{\Delta_R H(T)}{(J/mol) \cdot 10^3 }$')
ax2.legend(fontsize='small', loc='center left')
ax3 = plt.subplot2grid([2, 3], [1, 1], colspan=2)
ax3.set_ylabel('Konzentration / Mol%')
ax3.set_xlabel('Reduzierte Position, $z/L_R$')
for item in ['CO', 'H2O', 'MeOH', 'CO2']:
marker = np.random.choice(list(lines.lineMarkers.keys()))
index = namen.index(item)
ax3.plot(z_d_l_r, y_i_soln[:, index] * 100., label=item,
marker=marker)
ax3.legend(loc=1)
ax4 = plt.subplot2grid([2, 3], [0, 1])
ax4.set_ylabel('Temperatur / K')
ax4.set_xlabel('Reduzierte Position, $z/L_R$')
ax4.plot(z_d_l_r, t_soln, label='T / K')
ax5 = plt.subplot2grid([2, 3], [0, 2], colspan=2)
ax5.set_ylabel('Druck / bar')
ax5.set_xlabel('Reduzierte Position, $z/L_R$')
ax5.plot(z_d_l_r, p_soln, label='p / bar')
plt.tight_layout(rect=[0, 0, 0.95, 0.8])
plt.show()
In [2]:
# Chem. Eng. Technol. 2011, 34, No. 5, 817–822
# Katalysator
rho_b = 1190 # kg Kat/m^3 Feststoff
phi = 0.285 # m^3 Gas/m^3 Feststoff
m_kat = 1190 * (1 - 0.285) * np.pi / 4 * (
0.04)**2 * 7 # kg Kat (pro Rohr)
d_p = 0.0054 # m Feststoff
# Reaktor
n_t = 1620 # Rohre
d_t = 0.04 # m Rohrdurchmesser
l_r = 7. # m Rohrlänge
# Betriebsbedingungen
t0 = 225 + 273.15 # K
p0 = 69.7 # bar
m_dot = 57282.8 / 60**2 # kg/s
# Wärmetauschparameter
u = 118.44 # W/m^2/K
# Kühlmitteleigenschaften (VDI-WA)
t_r = (240 - 230) / (33.467 - 27.968) * (
29 - 33.467) + 240 + 273.15 # K
p_sat = 29 # bar
h_sat_l = (1037.5 - 990.21) / (33.467 - 27.968) * (
29 - 33.467) + 1037.5 # kJ/kg
h_sat_v = (2803.1 - 2803.0) / (33.467 - 27.968) * (
29 - 33.467) + 2803.1 # kJ/kg
delta_h_sat = (h_sat_v - h_sat_l)
# Zulaufbedingungen
namen = ['CO', 'CO2', 'H2', 'H2O', 'MeOH',
'CH4', 'N2', 'EthOH', 'PrOH', 'METHF']
m_dot_i = np.array([
10727.9, 23684.2, 9586.5,
108.8, 756.7, 4333.1,
8072.0, 0.6, 0.0,
13.0
], dtype=float) / 60**2 # kg/s
m_dot = sum(m_dot_i)
mm = np.array([
28.01, 44.01, 2.016,
18.015, 32.042, 16.043,
28.014, 46.069, 60.096,
60.053
], dtype=float) # g/mol
y_i0 = m_dot_i / mm / sum(m_dot_i / mm)
# Stoechiometrische Koeffizienten
nuij = np.zeros([len(namen), 3])
# Hydrierung von CO2
nuij[[
namen.index('CO2'),
namen.index('H2'),
namen.index('MeOH'),
namen.index('H2O'),
namen.index('CO'),
namen.index('N2'),
], 0] = np.array([-1, -3, +1, +1, 0, 0], dtype=float)
# Hydrierung von CO
nuij[[
namen.index('CO2'),
namen.index('H2'),
namen.index('MeOH'),
namen.index('H2O'),
namen.index('CO'),
namen.index('N2'),
], 1] = np.array([0, -2, +1, 0, -1, 0], dtype=float)
# RWGS Reverse-Wassergasshiftreaktion (muss gleich sein als die Vorwärtsreaktion,
# wofür die Kinetik verfügbar ist)
nuij[[
namen.index('CO2'),
namen.index('H2'),
namen.index('MeOH'),
namen.index('H2O'),
namen.index('CO'),
namen.index('N2'),
], 2] = - np.array([+1, +1, 0, -1, -1, 0], dtype=float)
# Parameter der Gleichung aus dem VDI-Wärmeatlas
cp_konstanten = np.zeros([len(namen), 7])
cp_konstanten_text = [
'407,9796 3,5028 2,8524 –2,3018 32,9055 –100,1815 106,1141',
'514,5073 3,4923 –0,9306 –6,0861 54,1586 –97,5157 70,9687',
'392,8422 2,4906 –3,6262 –1,9624 35,6197 –81,3691 62,6668',
'706,3032 5,1703 –6,0865 –6,6011 36,2723 –63,0965 46,2085',
'846,6321 5,7309 –4,8842 –12,8501 78,9997 –127,3725 82,7107',
'1530,8043 4,2038 –16,6150 –3,5668 43,0563 –86,5507 65,5986',
'432,2027 3,5160 2,8021 –4,1924 42,0153 –114,2500 111,1019',
'1165,8648 4,7021 9,7786 –1,1769 –135,7676 322,7596 –247,7349',
'506,0032 12,1539 0,0041 –36,1389 175,9466 –276,1927 171,3886',
'650,0705 5,0360 –0,4487 –13,5453 126,9502 –219,5695 151,4586'
]
for i in range(len(namen)):
cp_konstanten[i, :] = np.array(
cp_konstanten_text[i].replace(
'–', '-').replace(',', '.').split(' '),
dtype=float)
# Lennard-Jones Parameter Prausnitz Anhang B
# epsilon / kB
l_j_epsilon_d_k = np.array([
91.7, 195.2, 59.7,
809.1, 481.8, 148.6,
71.4, 362.6, 576.7,
469.8
]) # K
# sigma
l_j_sigma = np.array([
3.69, 3.941, 2.827,
2.641, 3.626, 3.758,
3.798, 4.53, 4.549,
4.936
]) # Angstrom
# T* = k T / epsilon
# Quelle: The properties of Gases and Liquids
h_298 = np.array([
-110.53, -393.51, 0,
-241.81, -200.94, -74.52,
0, -234.95, -255.20,
-352.40
], dtype=float) * 1000. # J/mol
delta_h_r_298 = nuij.T.dot(h_298) # J/mol
# Berechnung der Parameter
g = m_dot / (np.pi / 4 * d_t**2) / n_t # kg/m^2/s
mm_m_0 = sum(y_i0 * mm) * 1 / 1000. # kg/mol
cp_m_0 = sum(y_i0 * cp_ig_durch_r(t0) * 8.3145) # J/mol/K
cp_g_0 = cp_m_0 / mm_m_0 # J/kg/K
# Anzahl an Übertragungseinheiten (NTU)
ntu = l_r * 1 / (g * cp_g_0) * 2 * u / (d_t / 2)
# m * m^2 s/kg * kg K /J * J/s/m^2/K *1/m = [dimensionslose Einheiten]
# Stöchiometrische Zahl
sn = (y_i0[namen.index('H2')] - y_i0[namen.index('CO2')]) / (
y_i0[namen.index('CO2')] + y_i0[namen.index('CO')]
)
y_0 = np.empty([len(y_i0) + 1 + 1])
y_0[:-2] = y_i0
y_0[-2] = p0
y_0[-1] = t0
z_d_l_r = np.linspace(0, 1, 100)
d_z_dimlos = 1 / (len(z_d_l_r) - 1) # dimlos
dlr = d_z_dimlos * l_r # m
soln = odeint(df_dt, y_0, z_d_l_r)
y_i_soln = soln[:, :len(y_i0)]
p_soln = soln[:, -2]
t_soln = soln[:, -1]
n_i_soln = np.zeros_like(y_i_soln)
m_i_soln = np.zeros_like(y_i_soln)
n_soln = np.zeros_like(z_d_l_r)
mm_m_soln = np.zeros_like(z_d_l_r)
m_soln = np.zeros_like(z_d_l_r)
v_soln = np.zeros_like(z_d_l_r)
m_km_soln = np.zeros_like(z_d_l_r)
ums_soln = np.zeros_like(z_d_l_r)
for i in range(len(z_d_l_r)):
mm_m_soln[i] = sum(y_i_soln[i] * mm * 1 / 1000.) # kg/mol
n_soln[i] = g * n_t * 60**2 * (np.pi / 4 * d_t**2) / mm_m_soln[i]
# kg/s/m^2 * 60^2s/h * m^2 / kg*mol = mol/h
m_soln[i] = g * n_t * 60**2 * (np.pi / 4 * d_t**2) # kg/h
n_i_soln[i] = n_soln[i] * y_i_soln[i] # mol/h
m_i_soln[i] = n_soln[i] * y_i_soln[i] * (mm * 1 / 1000.)
# mol/h * g/mol * 1kg/1000g
v_soln[i] = n_soln[i] * 8.3145 * 1e-5 * t_soln[i] / p_soln[i]
ums_soln[i] = (n_i_soln[0][namen.index('CO')] -
n_i_soln[i][namen.index('CO')]
) / n_i_soln[0][namen.index('CO')]
m_km_soln[i] = u * (2 / (d_t / 2)) * (t_soln[i] - t_r) * (
np.pi / 4 * d_t**2) / delta_h_sat * n_t * 60**2 / 1000.
# J/s/K/m^2 * 1/m * K * m^2 * kg/kJ * 60^2s/h * 1kJ/(1000J) = kg/h/m
# Energie-Analyse
t_m_1 = 1 / 2 * (36696 / 8.3145 - np.sqrt(36696 /
8.3145 * (36696 / 8.3145 - 4 * t_r)))
t_m_2 = 1 / 2 * (-94765 / 8.3145 - np.sqrt(-94765 /
8.3145 * (-94765 / 8.3145 - 4 * t_r)))
vars_1 = [
[r'\rho_b', rho_b, r'\frac{kg_{Kat}}{m^3_{Schüttung}}'],
['\phi', phi, r'\frac{m^3_{Gas}}{m^3_{Schüttung}}'],
['D_p', d_p, 'm'],
['D_t', d_t, 'm'],
['L_R', l_r, 'm'],
]
vars_2 = [
['n_T', n_t, ''],
['U', u, r'\frac{W}{m^2\cdot K}'],
['\dot m', m_dot, 'kg/s'],
['C_{p_g}', cp_g_0 / 1000., r'\frac{kJ}{kg\cdot K}'],
['NTU', ntu, ''],
]
vars_3 = [
['T_0', t0 - 273.15, '°C'],
['P_0', p0, 'bar'],
['T_r', t_r - 273.15, '°C_{Kühlmittel}'],
['P_{Sät}', p_sat, 'bar_{Kühlmittel}'],
['\Delta H_{Sät}', delta_h_sat, r'\frac{kJ}{kg_{Kühlmittel}}'],
['SN', sn, '']
]
text_1 = '\n'.join(['$' + ' = '.join([line[0], '{:g}'.format(line[1]) +
' ' + line[2]]) + '$'
for line in vars_1])
text_2 = '\n'.join(['$' + ' = '.join([line[0], '{:g}'.format(line[1]) +
' ' + line[2]]) + '$'
for line in vars_2])
text_3 = '\n'.join(['$' + ' = '.join([line[0], '{:g}'.format(line[1]) +
' ' + line[2]]) + '$'
for line in vars_3])
fig = plt.figure(2, figsize=(8, 8))
fig.suptitle('Nicht adiabates System' +
'(Chem. Eng. Technol. 2011, 34, No. 5, 817–822)')
fig.text(0.05, 0.945, text_1, va='top', fontsize=10)
fig.text(0.33, 0.935, text_2, va='top', fontsize=10)
fig.text(0.66, 0.930, text_3, va='top', fontsize=10)
ax = plt.subplot2grid([2, 3], [0, 0])
ax.plot(z_d_l_r, v_soln, label='$\dot V$')
ax.set_ylabel(r'$\frac{\dot V}{m^3/h}$')
ax.set_xlabel('Reduzierte Position, $z/L_R$')
ax2 = plt.subplot2grid([2, 3], [1, 0])
ax2.plot(z_d_l_r, m_km_soln)
ax2.fill_between(z_d_l_r, 0, m_km_soln, color='orange')
ax2.text(0.3, 1 / 2. * (m_km_soln[0] + m_km_soln[-1]),
'{:g}'.format(sum(m_km_soln * dlr)) + 'kg/h \n')
ax2.set_ylabel(r'$\frac{\dot m_{Kuehlmittel}}{\frac{kg}{h\cdot m}}$')
ax2.set_xlabel('Reduzierte Position, $z/L_R$')
ax3 = plt.subplot2grid([2, 3], [1, 1], colspan=2)
ax3.set_ylabel('Massenstrom / (kg/h)')
ax3.set_xlabel('Reduzierte Position, $z/L_R$')
for item in ['CO', 'H2O', 'MeOH', 'CO2']:
marker = np.random.choice(list(lines.lineMarkers.keys()))
index = namen.index(item)
ax3.plot(z_d_l_r, m_i_soln[:, index], label=item,
marker=marker)
ax3.legend(loc=1)
ax4 = plt.subplot2grid([2, 3], [0, 1])
ax4_1 = ax4.twinx()
ax4_1.set_ylabel('Umsatz (CO)')
ax4.set_ylabel('Temperatur / °C')
ax4.set_xlabel('Reduzierte Position, $z/L_R$')
ax4.plot(z_d_l_r, t_soln - 273.15, label='T / °C')
ax4_1.plot(z_d_l_r, ums_soln, label='Umsatz (CO)', ls='--', color='gray')
ax5 = plt.subplot2grid([2, 3], [0, 2], colspan=2)
ax5.set_ylabel('Druck / bar')
ax5.set_xlabel('Reduzierte Position, $z/L_R$')
ax5.plot(z_d_l_r, p_soln, label='p / bar')
plt.tight_layout(rect=[0, 0, 0.95, 0.75])
print('')
print('NTU= ' + '{:g}'.format(ntu))
print('\n'.join([
namen[i] + ': ' + '{:g}'.format(x) + ' kg/h'
for i, x in enumerate(m_i_soln[-1])
]))
print('T=' + str(t_soln[-1] - 273.15) + '°C')
print('P=' + str(p_soln[-1]) + 'bar')
print('V0=' + str(v_soln[0]) + 'm^3/h')
print('V=' + str(v_soln[-1]) + 'm^3/h')
print('Kühlmittel: Gesättigtes $H_2O(l)$' +
' bei ' + '{:g}'.format(p_sat) + ' bar' +
'\n' + 'Verdampfungsenthalpie: ' +
'{:g}'.format(delta_h_sat) +
'kJ/kg' + '\n' + 'Kühlmittelmassenstrom: ' +
'{:g}'.format(sum(m_km_soln * dlr)) + 'kg/h')
print('Partikeldurchmesser für DeltaP=' +
'{:g}'.format(p_soln[0] - p_soln[-1]) + ' bar: ' +
'{:g}'.format(d_p) + ' m'
)
plt.show()
In [3]:
# Lösung der PAT-Übung 4
# Katalysator
rho_b = 1190 # kg Kat/m^3 Feststoff
phi = 0.3 # m^3 Gas/m^3 Feststoff
m_kat = 1190 * (1 - 0.3) * np.pi / 4 * (
0.03)**2 * 7 # kg Kat (pro Rohr)
# Partikeldurchmesser wählen, damit
# Delta P gesamt=-3bar, denn dies ist der Parameter
d_p = 0.0054 # m Feststoff
# Reaktor
n_t = 4800 # Rohre
d_t = 0.03 # m Rohrdurchmesser
l_r = 6. # m Rohrlänge
# Betriebsbedingungen
t0 = 220 + 273.15 # K
p0 = 50 # bar
# Wärmetauschparameter
u = 105 # W/m^2/K
# Kühlmitteleigenschaften (VDI-WA)
t_r = (240 - 230) / (33.467 - 27.968) * (
30 - 33.467) + 240 + 273.15 # K
p_sat = 30 # bar
h_sat_l = (1037.5 - 990.21) / (33.467 - 27.968) * (
30 - 33.467) + 1037.5 # kJ/kg
h_sat_v = (2803.1 - 2803.0) / (33.467 - 27.968) * (
30 - 33.467) + 2803.1 # kJ/kg
delta_h_sat = (h_sat_v - h_sat_l)
# Zulaufbedingungen
namen = ['CO', 'CO2', 'H2', 'H2O', 'MeOH',
'CH4', 'N2', 'EthOH', 'PrOH', 'METHF']
mm = np.array([
28.01, 44.01, 2.016,
18.015, 32.042, 16.043,
28.014, 46.069, 60.096,
60.053
], dtype=float) # g/mol
n_dot_i = np.array([
1386.42, 1871.83, 13085.1,
402.271, 112.237, 0,
1249.91, 0, 0, 0
]) * 1000 / 60**2 # mol/s
m_dot_i = n_dot_i * mm / 1000. # kg/s
m_dot = sum(m_dot_i)
y_i0 = m_dot_i / mm / sum(m_dot_i / mm)
# Berechnung der Parameter
g = m_dot / (np.pi / 4 * d_t**2) / n_t # kg/m^2/s
mm_m_0 = sum(y_i0 * mm) * 1 / 1000. # kg/mol
cp_m_0 = sum(y_i0 * cp_ig_durch_r(t0) * 8.3145) # J/mol/K
cp_g_0 = cp_m_0 / mm_m_0 # J/kg/K
# Anzahl an Übertragungseinheiten (NTU)
ntu = l_r * 1 / (g * cp_g_0) * 2 * u / (d_t / 2)
# m * m^2 s/kg * kg K /J * J/s/m^2/K *1/m = [dimensionslose Einheiten]
# Stöchiometrische Zahl
sn = (y_i0[namen.index('H2')] - y_i0[namen.index('CO2')]) / (
y_i0[namen.index('CO2')] + y_i0[namen.index('CO')]
)
# Restliche Koeffizienten und stoffliche Parameter
# bleiben ungeändert, wie oben.
y_0 = np.empty([len(y_i0) + 1 + 1])
y_0[:-2] = y_i0
y_0[-2] = p0
y_0[-1] = t0
z_d_l_r = np.linspace(0, 1, 100)
dlr = 1 / (len(z_d_l_r) - 1) * l_r # m
def opt_dp():
# Partikeldurchmesser nach Parameter Delta_P=3bar optimisieren.
# Es wird direkt Fixpunkt-Iteration angewendet, nach der Form der Ergun Gl.
# D_p{n+1} = D_p{n} * int(1/D_p{n}*f(D_p{n}) dz) / -3bar
global d_p
for i in range(5):
soln_dp = odeint(df_dt, y_0, z_d_l_r)
dp_dz = (soln_dp[-1, -2] - soln_dp[0, -2]).item()
d_p = dp_dz * d_p / -3.0
opt_dp()
soln = odeint(df_dt, y_0, z_d_l_r)
y_i_soln = soln[:, :len(y_i0)]
p_soln = soln[:, -2]
t_soln = soln[:, -1]
n_i_soln = np.zeros_like(y_i_soln)
m_i_soln = np.zeros_like(y_i_soln)
n_soln = np.zeros_like(z_d_l_r)
mm_m_soln = np.zeros_like(z_d_l_r)
m_soln = np.zeros_like(z_d_l_r)
v_soln = np.zeros_like(z_d_l_r)
m_km_soln = np.zeros_like(z_d_l_r)
ums_soln = np.zeros_like(z_d_l_r)
for i in range(len(z_d_l_r)):
mm_m_soln[i] = sum(y_i_soln[i] * mm * 1 / 1000.) # kg/mol
n_soln[i] = g * n_t * 60**2 * (np.pi / 4 * d_t**2) / mm_m_soln[i]
# kg/s/m^2 * 60^2s/h * m^2 / kg*mol = mol/h
m_soln[i] = g * n_t * 60**2 * (np.pi / 4 * d_t**2) # kg/h
n_i_soln[i] = n_soln[i] * y_i_soln[i] # mol/h
m_i_soln[i] = n_soln[i] * y_i_soln[i] * (mm * 1 / 1000.)
# mol/h * g/mol * 1kg/1000g
v_soln[i] = n_soln[i] * 8.3145 * 1e-5 * t_soln[i] / p_soln[i]
ums_soln[i] = (n_i_soln[0][namen.index('CO')] -
n_i_soln[i][namen.index('CO')]
) / n_i_soln[0][namen.index('CO')]
m_km_soln[i] = u * (2 / (d_t / 2)) * (t_soln[i] - t_r) * (
np.pi / 4 * d_t**2) / delta_h_sat * n_t * 60**2 / 1000.
# J/s/K/m^2 * 1/m * K * m^2 * kg/kJ * 60^2s/h * 1kJ/(1000J) = kg/h/m
# Energie-Analyse
t_m_1 = 1 / 2 * (36696 / 8.3145 - np.sqrt(36696 /
8.3145 * (36696 / 8.3145 - 4 * t_r)))
t_m_2 = 1 / 2 * (-94765 / 8.3145 - np.sqrt(-94765 /
8.3145 * (-94765 / 8.3145 - 4 * t_r)))
vars_1 = [
[r'\rho_b', rho_b, r'\frac{kg_{Kat}}{m^3_{Schüttung}}'],
['\phi', phi, r'\frac{m^3_{Gas}}{m^3_{Schüttung}}'],
['D_p', d_p, 'm'],
['D_t', d_t, 'm'],
['L_R', l_r, 'm'],
]
vars_2 = [
['n_T', n_t, ''],
['U', u, r'\frac{W}{m^2\cdot K}'],
['\dot m', m_dot, 'kg/s'],
['C_{p_g}', cp_g_0 / 1000., r'\frac{kJ}{kg\cdot K}'],
['NTU', ntu, ''],
]
vars_3 = [
['T_0', t0 - 273.15, '°C'],
['P_0', p0, 'bar'],
['T_r', t_r - 273.15, '°C_{Kühlmittel}'],
['P_{Sät}', p_sat, 'bar_{Kühlmittel}'],
['\Delta H_{Sät}', delta_h_sat, r'\frac{kJ}{kg_{Kühlmittel}}'],
['SN', sn, '']
]
text_1 = '\n'.join(['$' + ' = '.join([line[0], '{:g}'.format(line[1]) +
' ' + line[2]]) + '$'
for line in vars_1])
text_2 = '\n'.join(['$' + ' = '.join([line[0], '{:g}'.format(line[1]) +
' ' + line[2]]) + '$'
for line in vars_2])
text_3 = '\n'.join(['$' + ' = '.join([line[0], '{:g}'.format(line[1]) +
' ' + line[2]]) + '$'
for line in vars_3])
fig = plt.figure(5, figsize=(8, 8))
fig.suptitle('Lösung der PAT-Übung 4')
fig.text(0.05, 0.945, text_1, va='top', fontsize=10)
fig.text(0.33, 0.935, text_2, va='top', fontsize=10)
fig.text(0.66, 0.930, text_3, va='top', fontsize=10)
ax = plt.subplot2grid([2, 3], [0, 0])
ax.plot(z_d_l_r, v_soln, label='$\dot V$')
ax.set_ylabel(r'$\frac{\dot V}{m^3/h}$')
ax.set_xlabel('Reduzierte Position, $z/L_R$')
ax2 = plt.subplot2grid([2, 3], [1, 0])
ax2.plot(z_d_l_r, m_km_soln)
ax2.fill_between(z_d_l_r, 0, m_km_soln, color='orange')
ax2.text(0.3, 1 / 2. * (m_km_soln[0] + m_km_soln[-1]),
'{:g}'.format(sum(m_km_soln * dlr)) + 'kg/h')
ax2.set_ylabel(r'$\frac{\dot m_{Kuehlmittel}}{\frac{kg}{h\cdot m}}$')
ax2.set_xlabel('Reduzierte Position, $z/L_R$')
ax3 = plt.subplot2grid([2, 3], [1, 1], colspan=2)
ax3.set_ylabel('Massenstrom / (kg/h)')
ax3.set_xlabel('Reduzierte Position, $z/L_R$')
for item in ['CO', 'H2O', 'MeOH', 'CO2']:
marker = np.random.choice(list(lines.lineMarkers.keys()))
index = namen.index(item)
ax3.plot(z_d_l_r, m_i_soln[:, index], label=item,
marker=marker)
ax3.legend(loc=1)
ax4 = plt.subplot2grid([2, 3], [0, 1])
ax4_1 = ax4.twinx()
ax4_1.set_ylabel('Umsatz (CO)')
ax4.set_ylabel('Temperatur / °C')
ax4.set_xlabel('Reduzierte Position, $z/L_R$')
ax4.plot(z_d_l_r, t_soln - 273.15, label='T / °C')
ax4_1.plot(z_d_l_r, ums_soln, label='Umsatz (CO)', ls='--', color='gray')
ax5 = plt.subplot2grid([2, 3], [0, 2], colspan=2)
ax5.set_ylabel('Druck / bar')
ax5.set_xlabel('Reduzierte Position, $z/L_R$')
ax5.plot(z_d_l_r, p_soln, label='p / bar')
plt.tight_layout(rect=[0, 0, 0.95, 0.75])
print('')
print('NTU= ' + '{:g}'.format(ntu))
print('\n'.join([
namen[i] + ': ' + '{:g}'.format(x) + ' kg/h'
for i, x in enumerate(m_i_soln[-1])
]))
print('T=' + str(t_soln[-1] - 273.15) + '°C')
print('P=' + str(p_soln[-1]) + 'bar')
print('V0=' + str(v_soln[0]) + 'm^3/h')
print('V=' + str(v_soln[-1]) + 'm^3/h')
print('Kühlmittel: Gesättigtes $H_2O(l)$' +
' bei ' + '{:g}'.format(p_sat) + ' bar' +
'\n' + 'Verdampfungsenthalpie: ' +
'{:g}'.format(delta_h_sat) +
'kJ/kg' + '\n' + 'Kühlmittelmassenstrom: ' +
'{:g}'.format(sum(m_km_soln * dlr)) + 'kg/h')
print('Partikeldurchmesser für DeltaP=' +
'{:g}'.format(p_soln[0] - p_soln[-1]) + ' bar: ' +
'{:g}'.format(d_p) + ' m'
)
In [4]:
import textwrap
# Lösung der ChP-Übung : Zunächst ein falsch ausgelegtes System.
# Katalysator
rho_b = 1190 # kg Kat/m^3 Feststoff
phi = 0.3 # m^3 Gas/m^3 Feststoff
m_kat = 1190 * (1 - 0.3) * np.pi / 4 * (
0.03)**2 * 7 # kg Kat (pro Rohr)
# Partikeldurchmesser wählen, damit
# Delta P gesamt=-3bar, denn dies ist der Parameter
d_p = 0.0054 / 0.0054 * 0.16232576224693065 # m Feststoff
# Reaktor
n_t = 4800 # Rohre
d_t = 0.03 # m Rohrdurchmesser
# n_t = 10000 # Rohre
# d_t = 0.04 # m Rohrdurchmesser
l_r = 6. # m Rohrlänge
# Betriebsbedingungen
t0 = 220 + 273.15 # K
p0 = 50 # bar
# Wärmetauschparameter
u = 105 # W/m^2/K
# Kühlmitteleigenschaften (VDI-WA)
t_r = (240 - 230) / (33.467 - 27.968) * (
29 - 33.467) + 240 + 273.15 # K
p_sat = 29 # bar
h_sat_l = (1037.5 - 990.21) / (33.467 - 27.968) * (
29 - 33.467) + 1037.5 # kJ/kg
h_sat_v = (2803.1 - 2803.0) / (33.467 - 27.968) * (
29 - 33.467) + 2803.1 # kJ/kg
delta_h_sat = (h_sat_v - h_sat_l)
# Zulaufbedingungen
namen = ['CO', 'CO2', 'H2', 'H2O', 'MeOH',
'CH4', 'N2', 'EthOH', 'PrOH', 'METHF']
mm = np.array([
28.01, 44.01, 2.016,
18.015, 32.042, 16.043,
28.014, 46.069, 60.096,
60.053
], dtype=float) # g/mol
n_dot_i = np.array([
3394.174, 3767.258+1013.2, 65334.34,
230.9895, 845.4556, 0,
0, 0, 0,
0
]) * 1000 / 60**2 # / 3.09 * 0.460126 # mol/s
m_dot_i = n_dot_i * mm / 1000. # kg/s
m_dot = sum(m_dot_i)
y_i0 = m_dot_i / mm / sum(m_dot_i / mm)
# Berechnung der Parameter
g = m_dot / (np.pi / 4 * d_t**2) / n_t # kg/m^2/s
mm_m_0 = sum(y_i0 * mm) * 1 / 1000. # kg/mol
cp_m_0 = sum(y_i0 * cp_ig_durch_r(t0) * 8.3145) # J/mol/K
cp_g_0 = cp_m_0 / mm_m_0 # J/kg/K
# Anzahl an Übertragungseinheiten (NTU)
ntu = l_r * 1 / (g * cp_g_0) * 2 * u / (d_t / 2)
# m * m^2 s/kg * kg K /J * J/s/m^2/K *1/m = [dimensionslose Einheiten]
# Stöchiometrische Zahl
sn = (y_i0[namen.index('H2')] - y_i0[namen.index('CO2')]) / (
y_i0[namen.index('CO2')] + y_i0[namen.index('CO')]
)
ratio_h2_co2 = y_i0[namen.index('H2')] / y_i0[namen.index('CO2')]
verhaeltnis_co_co2 = y_i0[namen.index('CO')] / y_i0[namen.index('CO2')]
# Restliche Koeffizienten und stoffliche Parameter
# bleiben ungeändert, wie oben.
y_0 = np.empty([len(y_i0) + 1 + 1])
y_0[:-2] = y_i0
y_0[-2] = p0
y_0[-1] = t0
z_d_l_r = np.linspace(0, 1, 100)
dlr = 1 / (len(z_d_l_r) - 1) * l_r # m
def opt_dp():
# Partikeldurchmesser nach Parameter Delta_P=3bar optimisieren.
# Es wird direkt Fixpunkt-Iteration angewendet, nach der Form der Ergun Gl.
# D_p{n+1} = D_p{n} * int(1/D_p{n}*f(D_p{n}) dz) / -3bar
global d_p
for i in range(5):
soln_dp = odeint(df_dt, y_0, z_d_l_r)
dp_dz = (soln_dp[-1, -2] - soln_dp[0, -2]).item()
d_p = dp_dz * d_p / -3.0
opt_dp()
soln = odeint(df_dt, y_0, z_d_l_r)
y_i_soln = soln[:, :len(y_i0)]
p_soln = soln[:, -2]
t_soln = soln[:, -1]
n_i_soln = np.zeros_like(y_i_soln)
m_i_soln = np.zeros_like(y_i_soln)
n_soln = np.zeros_like(z_d_l_r)
mm_m_soln = np.zeros_like(z_d_l_r)
m_soln = np.zeros_like(z_d_l_r)
v_soln = np.zeros_like(z_d_l_r)
m_km_soln = np.zeros_like(z_d_l_r)
ums_soln = np.zeros_like(z_d_l_r)
for i in range(len(z_d_l_r)):
mm_m_soln[i] = sum(y_i_soln[i] * mm * 1 / 1000.) # kg/mol
n_soln[i] = g * n_t * 60**2 * (np.pi / 4 * d_t**2) / mm_m_soln[i]
# kg/s/m^2 * 60^2s/h * m^2 / kg*mol = mol/h
m_soln[i] = g * n_t * 60**2 * (np.pi / 4 * d_t**2) # kg/h
n_i_soln[i] = n_soln[i] * y_i_soln[i] # mol/h
m_i_soln[i] = n_soln[i] * y_i_soln[i] * (mm * 1 / 1000.)
# mol/h * g/mol * 1kg/1000g
v_soln[i] = n_soln[i] * 8.3145 * 1e-5 * t_soln[i] / p_soln[i]
ums_soln[i] = (n_i_soln[0][namen.index('CO')] -
n_i_soln[i][namen.index('CO')]
) / n_i_soln[0][namen.index('CO')]
m_km_soln[i] = u * (2 / (d_t / 2)) * (t_soln[i] - t_r) * (
np.pi / 4 * d_t**2) / delta_h_sat * n_t * 60**2 / 1000.
# J/s/K/m^2 * 1/m * K * m^2 * kg/kJ * 60^2s/h * 1kJ/(1000J) = kg/h/m
# Energie-Analyse
t_m_1 = 1 / 2 * (36696 / 8.3145 - np.sqrt(36696 /
8.3145 * (36696 / 8.3145 - 4 * t_r)))
t_m_2 = 1 / 2 * (-94765 / 8.3145 - np.sqrt(-94765 /
8.3145 * (-94765 / 8.3145 - 4 * t_r)))
vars_1 = [
[r'\rho_b', rho_b, r'\frac{kg_{Kat}}{m^3_{Schüttung}}'],
['\phi', phi, r'\frac{m^3_{Gas}}{m^3_{Schüttung}}'],
['D_p', d_p, 'm'],
['D_t', d_t, 'm'],
['L_R', l_r, 'm'],
]
vars_2 = [
['n_T', n_t, ''],
['U', u, r'\frac{W}{m^2\cdot K}'],
['\dot m', m_dot, 'kg/s'],
['C_{p_g}', cp_g_0 / 1000., r'\frac{kJ}{kg\cdot K}'],
['NTU', ntu, ''],
]
vars_3 = [
['T_0', t0 - 273.15, '°C'],
['P_0', p0, 'bar'],
['T_r', t_r - 273.15, '°C_{Kühlmittel}'],
['P_{Sät}', p_sat, 'bar_{Kühlmittel}'],
['\Delta H_{Sät}', delta_h_sat, r'\frac{kJ}{kg_{Kühlmittel}}'],
['SN', sn, '']
]
text_1 = '\n'.join(['$' + ' = '.join([line[0], '{:g}'.format(line[1]) +
' ' + line[2]]) + '$'
for line in vars_1])
text_2 = '\n'.join(['$' + ' = '.join([line[0], '{:g}'.format(line[1]) +
' ' + line[2]]) + '$'
for line in vars_2])
text_3 = '\n'.join(['$' + ' = '.join([line[0], '{:g}'.format(line[1]) +
' ' + line[2]]) + '$'
for line in vars_3])
fig = plt.figure(6, figsize=(8, 8))
fig.suptitle('Lösung der Zusammensetzung \n' +
'{:d}'.format(round(ratio_h2_co2.item())) +
':1:' +
'{:d}'.format(round(verhaeltnis_co_co2.item())) +
'(H2:CO2:CO)')
fig.text(0, 0.65, textwrap.fill('Für den Massenstrom von' +
r' $\dot m = ' + '{:g}'.format(m_dot * n_t * 60**2) +
r'\frac{kg}{h}$ falsch dimensioniert. Um NTU=3,09' +
' einzustellen sind folgende Variablen zu ändern:' +
'i) höherer Radius, ii) höhere Länge, iii) höherer ' +
'Wärmedurchgangskoeffizient U, iv) niedrigerer Massenstrom (oder' +
' dementsprechend mehr Rohre), '
'iv) niedrigere Wärmekapazität Cpg (über die Zusammensetzung). ' +
'Letzteres kann eigentlich weggelassen werden, da ' +
'H2 schon den höchsten Molanteil beträgt.', 110),
wrap=True,
fontdict=dict((('size', 10), ('family', 'monospace'),
('ha', 'left'), ('va', 'bottom')))
)
fig.text(0.05, 0.945, text_1, va='top', fontsize=10)
fig.text(0.33, 0.935, text_2, va='top', fontsize=10)
fig.text(0.66, 0.930, text_3, va='top', fontsize=10)
ax = plt.subplot2grid([2, 3], [0, 0])
ax.plot(z_d_l_r, v_soln, label='$\dot V$')
ax.set_ylabel(r'$\frac{\dot V}{m^3/h}$')
ax.set_xlabel('Reduzierte Position, $z/L_R$')
ax2 = plt.subplot2grid([2, 3], [1, 0])
ax2.plot(z_d_l_r, m_km_soln)
ax2.fill_between(z_d_l_r, 0, m_km_soln, color='orange')
ax2.text(0.3, 1 / 2. * (m_km_soln[0] + m_km_soln[-1]),
'{:g}'.format(sum(m_km_soln * dlr)) + 'kg/h')
ax2.set_ylabel(r'$\frac{\dot m_{Kuehlmittel}}{kg/h}$')
ax2.set_xlabel('Reduzierte Position, $z/L_R$')
ax3 = plt.subplot2grid([2, 3], [1, 1], colspan=2)
ax3.set_ylabel('Massenstrom / (kg/h)')
ax3.set_xlabel('Reduzierte Position, $z/L_R$')
for item in ['CO', 'H2O', 'MeOH', 'CO2']:
marker = np.random.choice(list(lines.lineMarkers.keys()))
index = namen.index(item)
ax3.plot(z_d_l_r, m_i_soln[:, index], label=item,
marker=marker)
ax3.legend(loc=1)
ax4 = plt.subplot2grid([2, 3], [0, 1])
ax4_1 = ax4.twinx()
ax4_1.set_ylabel('Umsatz (CO)')
ax4.set_ylabel('Temperatur / °C')
ax4.set_xlabel('Reduzierte Position, $z/L_R$')
ax4.plot(z_d_l_r, t_soln - 273.15, label='T / °C')
ax4_1.plot(z_d_l_r, ums_soln, label='Umsatz (CO)', ls='--', color='gray')
ax5 = plt.subplot2grid([2, 3], [0, 2], colspan=2)
ax5.set_ylabel('Druck / bar')
ax5.set_xlabel('Reduzierte Position, $z/L_R$')
ax5.plot(z_d_l_r, p_soln, label='p / bar')
plt.tight_layout(rect=[0, 0, 0.95, 0.65])
In [5]:
# Lösung der ChP-Übung. Massenstrom angepasst.
# Katalysator
rho_b = 1190 # kg Kat/m^3 Feststoff
phi = 0.3 # m^3 Gas/m^3 Feststoff
m_kat = 1190 * (1 - 0.3) * np.pi / 4 * (
0.03)**2 * 7 # kg Kat (pro Rohr)
# Partikeldurchmesser wählen, damit
# Delta P gesamt=-3bar, denn dies ist der Parameter
d_p = 0.0054 / 0.0054 * 0.16232576224693065 # m Feststoff
d_p = 0.0037 # m
# Reaktor
n_t = 4800 * 1.82 / 0.460126 # Rohre
n_t = 10921 # Rohre
d_t = 0.03 # m Rohrdurchmesser
d_t = 0.04 # m Rohrdurchmesser
l_r = 6. # m Rohrlänge
# Betriebsbedingungen
t0 = 220 + 273.15 # K
p0 = 50 # bar
# Wärmetauschparameter
u = 105 # W/m^2/K
u = 118.44 # W/m^2/K
# Kühlmitteleigenschaften (VDI-WA)
t_r = (240 - 230) / (33.467 - 27.968) * (
29 - 33.467) + 240 + 273.15 # K
p_sat = 29 # bar
h_sat_l = (1037.5 - 990.21) / (33.467 - 27.968) * (
29 - 33.467) + 1037.5 # kJ/kg
h_sat_v = (2803.1 - 2803.0) / (33.467 - 27.968) * (
29 - 33.467) + 2803.1 # kJ/kg
delta_h_sat = (h_sat_v - h_sat_l)
# Zulaufbedingungen
namen = ['CO', 'CO2', 'H2', 'H2O', 'MeOH',
'CH4', 'N2', 'EthOH', 'PrOH', 'METHF']
mm = np.array([
28.01, 44.01, 2.016,
18.015, 32.042, 16.043,
28.014, 46.069, 60.096,
60.053
], dtype=float) # g/mol
# Berechnete Daten
n_dot_i = np.array([
2462.31148689,
6223.84452606,
28373.5423591,
357.063600003,
564.402901433, 0, 0, 0, 0,
0
]) / 60**2 * 1000 # mol/s
# Daten aus Aspen
n_dot_i = np.array([
2206.10385252181, 5488.76285475286,
22222.0112180475, 252.168597341323,
252.168597341323, 0, 0, 0, 0,
0
]) / 60**2 * 1000 # mol/s
m_dot_i = n_dot_i * mm / 1000. # kg/s
m_dot = sum(m_dot_i)
y_i0 = m_dot_i / mm / sum(m_dot_i / mm)
# Berechnung der Parameter
g = m_dot / (np.pi / 4 * d_t**2) / n_t # kg/m^2/s
mm_m_0 = sum(y_i0 * mm) * 1 / 1000. # kg/mol
cp_m_0 = sum(y_i0 * cp_ig_durch_r(t0) * 8.3145) # J/mol/K
cp_g_0 = cp_m_0 / mm_m_0 # J/kg/K
# Anzahl an Übertragungseinheiten (NTU)
ntu = l_r * 1 / (g * cp_g_0) * 2 * u / (d_t / 2)
# m * m^2 s/kg * kg K /J * J/s/m^2/K *1/m = [dimensionslose Einheiten]
# Stöchiometrische Zahl
sn = (y_i0[namen.index('H2')] - y_i0[namen.index('CO2')]) / (
y_i0[namen.index('CO2')] + y_i0[namen.index('CO')]
)
ratio_h2_co2 = y_i0[namen.index('H2')] / y_i0[namen.index('CO2')]
verhaeltnis_co_co2 = y_i0[namen.index('CO')] / y_i0[namen.index('CO2')]
# Restliche Koeffizienten und stoffliche Parameter
# bleiben ungeändert, wie oben.
y_0 = np.empty([len(y_i0) + 1 + 1])
y_0[:-2] = y_i0
y_0[-2] = p0
y_0[-1] = t0
z_d_l_r = np.linspace(0, 1, 100)
dlr = 1 / (len(z_d_l_r) - 1) * l_r # m
def opt_dp():
# Partikeldurchmesser nach Parameter Delta_P=3bar optimisieren.
# Es wird direkt Fixpunkt-Iteration angewendet, nach der Form der Ergun Gl.
# D_p{n+1} = D_p{n} * int(1/D_p{n}*f(D_p{n}) dz) / -3bar
global d_p
for i in range(5):
soln_dp = odeint(df_dt, y_0, z_d_l_r)
dp_dz = (soln_dp[-1, -2] - soln_dp[0, -2]).item()
d_p = dp_dz * d_p / -3.0
# opt_dp() # momentan Partikeldurchmesser-Optimisierung
# überspringen
soln = odeint(df_dt, y_0, z_d_l_r)
y_i_soln = soln[:, :len(y_i0)]
p_soln = soln[:, -2]
t_soln = soln[:, -1]
n_i_soln = np.zeros_like(y_i_soln)
m_i_soln = np.zeros_like(y_i_soln)
n_soln = np.zeros_like(z_d_l_r)
mm_m_soln = np.zeros_like(z_d_l_r)
m_soln = np.zeros_like(z_d_l_r)
v_soln = np.zeros_like(z_d_l_r)
m_km_soln = np.zeros_like(z_d_l_r)
ums_soln = np.zeros_like(z_d_l_r)
for i in range(len(z_d_l_r)):
mm_m_soln[i] = sum(y_i_soln[i] * mm * 1 / 1000.) # kg/mol
n_soln[i] = g * n_t * 60**2 * (np.pi / 4 * d_t**2) / mm_m_soln[i]
# kg/s/m^2 * 60^2s/h * m^2 / kg*mol = mol/h
m_soln[i] = g * n_t * 60**2 * (np.pi / 4 * d_t**2) # kg/h
n_i_soln[i] = n_soln[i] * y_i_soln[i] # mol/h
m_i_soln[i] = n_soln[i] * y_i_soln[i] * (mm * 1 / 1000.)
# mol/h * g/mol * 1kg/1000g
v_soln[i] = n_soln[i] * 8.3145 * 1e-5 * t_soln[i] / p_soln[i]
ums_soln[i] = (n_i_soln[0][namen.index('CO')] -
n_i_soln[i][namen.index('CO')]
) / n_i_soln[0][namen.index('CO')]
m_km_soln[i] = u * (2 / (d_t / 2)) * (t_soln[i] - t_r) * (
np.pi / 4 * d_t**2) / delta_h_sat * n_t * 60**2 / 1000.
# J/s/K/m^2 * 1/m * K * m^2 * kg/kJ * 60^2s/h * 1kJ/(1000J) = kg/h/m
# Energie-Analyse
t_m_1 = 1 / 2 * (36696 / 8.3145 - np.sqrt(36696 /
8.3145 * (36696 / 8.3145 - 4 * t_r)))
t_m_2 = 1 / 2 * (-94765 / 8.3145 - np.sqrt(-94765 /
8.3145 * (-94765 / 8.3145 - 4 * t_r)))
vars_1 = [
[r'\rho_b', rho_b, r'\frac{kg_{Kat}}{m^3_{Schüttung}}'],
['\phi', phi, r'\frac{m^3_{Gas}}{m^3_{Schüttung}}'],
['D_p', d_p, 'm'],
['D_t', d_t, 'm'],
['L_R', l_r, 'm'],
]
vars_2 = [
['n_T', n_t, ''],
['U', u, r'\frac{W}{m^2\cdot K}'],
['\dot m', m_dot, 'kg/s'],
['C_{p_g}', cp_g_0 / 1000., r'\frac{kJ}{kg\cdot K}'],
['NTU', ntu, ''],
]
vars_3 = [
['T_0', t0 - 273.15, '°C'],
['P_0', p0, 'bar'],
['T_r', t_r - 273.15, '°C_{Kühlmittel}'],
['P_{Sät}', p_sat, 'bar_{Kühlmittel}'],
['\Delta H_{Sät}', delta_h_sat, r'\frac{kJ}{kg_{Kühlmittel}}'],
['SN', sn, '']
]
text_1 = '\n'.join(['$' + ' = '.join([line[0], '{:g}'.format(line[1]) +
' ' + line[2]]) + '$'
for line in vars_1])
text_2 = '\n'.join(['$' + ' = '.join([line[0], '{:g}'.format(line[1]) +
' ' + line[2]]) + '$'
for line in vars_2])
text_3 = '\n'.join(['$' + ' = '.join([line[0], '{:g}'.format(line[1]) +
' ' + line[2]]) + '$'
for line in vars_3])
fig = plt.figure(7, figsize=(8, 8))
fig.suptitle('Lösung der Zusammensetzung ' +
'{:.2f}'.format(round(ratio_h2_co2.item(), 2)) +
':1:' +
'{:.2f}'.format(round(verhaeltnis_co_co2.item(), 2)) +
'(H2:CO2:CO)')
fig.text(0.05, 0.945, text_1, va='top', fontsize=10)
fig.text(0.33, 0.935, text_2, va='top', fontsize=10)
fig.text(0.66, 0.930, text_3, va='top', fontsize=10)
fig.text(0.33, 0.725, 'System jetzt richtig dimensioniert',
fontdict=dict((('size', 7.5), ('family', 'monospace'),
('ha', 'left'), ('va', 'bottom')))
)
ax = plt.subplot2grid([2, 3], [0, 0])
ax.plot(z_d_l_r, v_soln, label='$\dot V$')
ax.set_ylabel(r'$\frac{\dot V}{m^3/h}$')
ax.set_xlabel('Reduzierte Position, $z/L_R$')
ax2 = plt.subplot2grid([2, 3], [1, 0])
ax2.plot(z_d_l_r, m_km_soln)
ax2.fill_between(z_d_l_r, 0, m_km_soln, color='orange')
ax2.text(0.3, 1 / 2. * (m_km_soln[0] + m_km_soln[-1]),
'{:g}'.format(sum(m_km_soln * dlr)) + 'kg/h')
ax2.set_ylabel(r'$\frac{\dot m_{Kuehlmittel}}{kg/h}$')
ax2.set_xlabel('Reduzierte Position, $z/L_R$')
ax3 = plt.subplot2grid([2, 3], [1, 1], colspan=2)
ax3.set_ylabel('Massenstrom / (kg/h)')
ax3.set_xlabel('Reduzierte Position, $z/L_R$')
for item in ['CO', 'H2O', 'MeOH', 'CO2']:
marker = np.random.choice(list(lines.lineMarkers.keys()))
index = namen.index(item)
ax3.plot(z_d_l_r, m_i_soln[:, index], label=item,
marker=marker)
ax3.legend(loc=1)
ax4 = plt.subplot2grid([2, 3], [0, 1])
ax4_1 = ax4.twinx()
ax4_1.set_ylabel('Umsatz (CO)')
ax4.set_ylabel('Temperatur / °C')
ax4.set_xlabel('Reduzierte Position, $z/L_R$')
ax4.plot(z_d_l_r, t_soln - 273.15, label='T / °C')
ax4_1.plot(z_d_l_r, ums_soln, label='Umsatz (CO)', ls='--', color='gray')
ax5 = plt.subplot2grid([2, 3], [0, 2], colspan=2)
ax5.set_ylabel('Druck / bar')
ax5.set_xlabel('Reduzierte Position, $z/L_R$')
ax5.plot(z_d_l_r, p_soln, label='p / bar')
plt.tight_layout(rect=[0, 0, 0.95, 0.75])
print('')
print('NTU= ' + '{:g}'.format(ntu))
print('\n'.join([
namen[i] + ': ' + '{:g}'.format(x) + ' kg/h'
for i, x in enumerate(m_i_soln[-1])
]))
print('T=' + str(t_soln[-1] - 273.15) + '°C')
print('P=' + str(p_soln[-1]) + 'bar')
print('V0=' + str(v_soln[0]) + 'm^3/h')
print('V=' + str(v_soln[-1]) + 'm^3/h')
print('Kühlmittel: Gesättigtes $H_2O(l)$' +
' bei ' + '{:g}'.format(p_sat) + ' bar' + '\n' + 'Verdampfungsenthalpie: ' +
'{:g}'.format(delta_h_sat) +
'kJ/kg' + '\n' + 'Kühlmittelmassenstrom: ' +
'{:g}'.format(sum(m_km_soln * dlr)) + 'kg/h')
print('Partikeldurchmesser für DeltaP=' +
'{:g}'.format(p_soln[0] - p_soln[-1]) + ' bar: ' +
'{:g}'.format(d_p) + ' m'
)
plt.tight_layout(rect=[0, 0, 0.95, 0.75])
In [6]:
from matplotlib.lines import Line2D
t_m_1 = 1 / 2 * (36696 / 8.3145 - np.sqrt(36696 /
8.3145 * (36696 / 8.3145 - 4 * t_r)))
t_m_2 = 1 / 2 * (-94765 / 8.3145 - np.sqrt(-94765 /
8.3145 * (-94765 / 8.3145 - 4 * t_r)))
fig = plt.figure(3, figsize=[10, 20])
fig.suptitle('Nicht adiabates System' +
'\nPhasenebenen' +
'\nCO/CO2 Verhältnis' +
'\nbei $p_{CO}+p_{CO_2}= $' +
'{:g}'.format(p0 * sum(
y_i0[[namen.index('CO'),
namen.index('CO2')]])) + ' bar',
x=0.8)
ax = plt.subplot(321)
ax.set_xlabel(r'T / K')
ax.set_ylabel('$p_{CO} / bar$')
ax2 = plt.subplot(322)
ax2.set_xlabel(r'T / K')
ax2.set_ylabel('$p_{CO_2} / bar$')
ax3 = plt.subplot(325)
ax3.set_xlabel(r'Reduzierte Position, $z/L_R$')
ax3.set_ylabel('$T / K$')
ax3.set_xlim([0, 1])
ax4 = plt.subplot(323)
ax4.set_xlabel(r'Reduzierte Position, $z/L_R$')
ax4.set_ylabel('$X(CO)$')
ax4.set_xlim([0, 1])
ax4.set_ylim([0, ax4.get_ylim()[1]])
ax5 = plt.subplot(324)
ax5.set_xlabel(r'Reduzierte Position, $z/L_R$')
ax5.set_ylabel('$log(Y(CH_3OH/CO))$')
ax5.set_xlim([0, 1])
ax5.set_ylim([0.5, 1000])
ax5.set_yscale('log')
ax6 = plt.subplot(326)
ax6.set_xlabel(r'Reduzierte Position, $z/L_R$')
ax6.set_ylabel(r'$\frac{\dot m(CH_3OH/CO)}{kg/h}$')
ax6.set_xlim([0, 1])
indexes_co_co2 = [namen.index('CO'), namen.index('CO2')]
indexes_not_co_co2 = [
i for i in range(len(y_i0))
if i not in [namen.index('CO'), namen.index('CO2')]
]
basis_ratio = y_i0[namen.index('CO')] / y_i0[namen.index('CO2')]
y_co_plus_co2 = sum(y_i0[[namen.index('CO'), namen.index('CO2')]])
lines1 = []
base_cm_1 = plt.get_cmap('viridis')
base_cm_2 = plt.get_cmap('inferno')
color_samples = np.array([
*base_cm_1(np.linspace(0, 1, 32))[:, :3],
*base_cm_2(np.linspace(0, 1, 32))[:, :3]
])
unfilled_markers = [m for m in Line2D.markers
if type(m) == str and
m not in ['None', 'nothing', None, ' ', '']]
for verhaeltnis_co_co2 in [1e-6, 0.15, 0.45, basis_ratio, 1.5, 7.5, 15]:
y_0[:-2][[namen.index('CO'), namen.index('CO2')]] = np.array(
[verhaeltnis_co_co2 / (1 + verhaeltnis_co_co2), 1 / (1 + verhaeltnis_co_co2)]) * y_co_plus_co2
p0_co2_co = sum(y_0[:-2][indexes_co_co2] * p0).item()
# Keep mass of all other components m' = m (1-y1-y2)/(1-y1'-y2')
g_new = (1 - sum(y_i0[indexes_co_co2])) / (
1 - sum(y_0[:-2][indexes_co_co2])
) * sum(y_0[:-2] * mm / 1000.) / sum(y_i0 * mm / 1000.) * g
soln = odeint(df_dt, y_0, z_d_l_r)
y_i_soln = soln[:, :len(y_i0)]
p_soln = soln[:, -2]
t_soln = soln[:, -1]
mm_m_soln = np.sum(y_i_soln * mm * 1 / 1000., axis=1) # kg/mol
n_soln = g_new * n_t * 60 ** 2 * (np.pi / 4 * d_t ** 2) / mm_m_soln
m_soln = g_new * n_t * 60 ** 2 * (np.pi / 4 * d_t ** 2) # kg/h
m_i_soln = np.array([
n_soln * y_i_soln[:, i] * (mm[i] * 1 / 1000.)
for i in range(y_i_soln.shape[1])]).T
ums_soln = (y_co_plus_co2 * verhaeltnis_co_co2 / (1 + verhaeltnis_co_co2) * n_soln[0] - (
y_i_soln[:, namen.index('CO')]) * n_soln) / (
y_co_plus_co2 * verhaeltnis_co_co2 / (1 + verhaeltnis_co_co2) * n_soln[0]
)
ausb_soln = (y_i_soln[1:, namen.index('MeOH')] * n_soln[1:] -
y_i_soln[0, namen.index('MeOH')] * n_soln[0]) / (
ums_soln[1:] * y_co_plus_co2 * verhaeltnis_co_co2 / (1 + verhaeltnis_co_co2) * n_soln[0])
# Anzahl an Übertragungseinheiten (NTU)
mm_m_0 = sum(y_0[:-2] * mm) * 1 / 1000. # kg/mol
cp_m_0 = sum(y_0[:-2] * cp_ig_durch_r(t0) * 8.3145) # J/mol/K
cp_g_0 = cp_m_0 / mm_m_0 # J/kg/K
ntu = l_r * 1 / (g * cp_g_0) * 2 * u / (d_t / 2)
sn = (y_0[:-2][namen.index('H2')] - y_0[:-2][namen.index('CO2')]) / (
y_0[:-2][namen.index('CO2')] + y_0[:-2][namen.index('CO')]
)
# m * m^2 s/kg * kg K /J * J/s/m^2/K *1/m = [dimensionslose Einheiten]
marker = np.random.choice(unfilled_markers)
color = color_samples[np.random.choice(len(color_samples))]
marker_color = color_samples[np.random.choice(len(color_samples))]
marker_fill = np.random.choice(list(lines.fillStyles))
line = ax.plot(t_soln, p_soln * y_i_soln[:, namen.index('CO')],
label='$p_{CO}/p_{CO_2}=' + str(verhaeltnis_co_co2) +
'\quad NTU= ' + '{:g}'.format(ntu) +
'\quad SN= ' + '{:g}'.format(sn) + '$',
marker=marker, markersize=3, color=color,
fillstyle=marker_fill, markerfacecolor=marker_color)[0]
lines1.append(line)
ax2.plot(t_soln, p_soln * y_i_soln[:, namen.index('CO2')],
label='$p_{CO}/p_{CO_2}=' + str(verhaeltnis_co_co2) + '$',
marker=marker, markersize=5, color=color,
fillstyle=marker_fill, markerfacecolor=marker_color)
ax3.plot(z_d_l_r, t_soln,
label='$p_{CO}/p_{CO_2}=' + str(verhaeltnis_co_co2) + '$',
marker=marker, markersize=5, color=color,
fillstyle=marker_fill, markerfacecolor=marker_color)
ax4.plot(z_d_l_r, ums_soln,
label='$p_{CO}/p_{CO_2}=' + str(verhaeltnis_co_co2) + '$',
marker=marker, markersize=5, color=color,
fillstyle=marker_fill, markerfacecolor=marker_color)
ax5.plot(z_d_l_r[1:], ausb_soln,
label='$p_{CO}/p_{CO_2}=' + str(verhaeltnis_co_co2) + '$',
marker=marker, markersize=5, color=color,
fillstyle=marker_fill, markerfacecolor=marker_color)
ax6.plot(z_d_l_r, m_i_soln[:, namen.index('MeOH')],
label='$p_{CO}/p_{CO_2}=' + str(verhaeltnis_co_co2) + '$',
marker=marker, markersize=5, color=color,
fillstyle=marker_fill, markerfacecolor=marker_color)
ax.plot([t_m_1, t_m_1], ax.get_ylim(), color='black', ls='--')
ax2.plot([t_m_1, t_m_1], ax2.get_ylim(), color='black', ls='--')
ax3.plot(ax2.get_xlim(), [t_m_1, t_m_1], color='black', ls='--')
fig.legend(lines1, [x.get_label() for x in lines1], 'upper left',
fontsize='large')
plt.tight_layout(rect=[0, 0, 1, 0.85])
plt.show()
Definition
$SN=\frac{H-CO_2}{CO_2+CO}$
Eine stöchiometrische Zahl leicht höher als 2 ist wünschenswert [8], und zwar wird 2,05 empfohlen [3].
Die Kennzahl besagt, dass genug Wasserstoff eingeführt wird, um das CO2 umzusetzen. Ein wesentlicher Überschuss an Wasserstoff ist auch von Nachteil, denn er führt bei niedrigerer Wärmekapazität konvektiv Wärme ab, und kann verursachen, dass der Temperaturprofil nicht erreicht werde.
In [46]:
fig = plt.figure(3, figsize=[10, 20])
fig.suptitle('Nicht adiabates System' +
'\nPhasenebenen' +
'\nStöchiometrische Zahl' +
'\nbei $p_{CO}/p_{CO_2}$ Verhältnis ' +
'{:g}'.format(y_i0[namen.index('CO')] / y_i0[namen.index('CO2')]) +
'\n$L_R = ' + '{:g}'.format(l_r) + 'm$',
horizontalalignment='left', x=0.8,
fontsize='xx-large')
ax = plt.subplot(321)
ax.set_xlabel(r'T / K')
ax.set_ylabel('$p_{CO} / bar$')
ax2 = plt.subplot(322)
ax2.set_xlabel(r'T / K')
ax2.set_ylabel('$p_{CO_2} / bar$')
ax3 = plt.subplot(325)
ax3.set_xlabel(r'Reduzierte Position, $z/L_R$')
ax3.set_ylabel('$T / K$')
ax3.set_xlim([0, 1])
ax4 = plt.subplot(323)
ax4.set_xlabel(r'Reduzierte Position, $z/L_R$')
ax4.set_ylabel('$X(CO)$')
ax4.set_xlim([0, 1])
ax4.set_ylim([0, 0.4])
ax5 = plt.subplot(324)
ax5.set_xlabel(r'Reduzierte Position, $z/L_R$')
ax5.set_ylabel('$log(Y(CH_3OH/CO))$')
ax5.set_xlim([0, 1])
ax5.set_yscale('symlog')
ax6 = plt.subplot(326)
ax6.set_xlabel(r'Reduzierte Position, $z/L_R$')
ax6.set_ylabel(r'$\frac{\dot m(CH_3OH/CO)}{kg/h}$')
ax6.set_xlim([0, 1])
verhaeltnis_co_co2 = y_i0[namen.index('CO')] / y_i0[namen.index('CO2')]
ratio_h2_co2 = y_i0[namen.index('H2')] / y_i0[namen.index('CO2')]
basis_y_co_plus_co2 = sum(y_i0[[namen.index('CO'), namen.index('CO2')]])
# Theoretisch pptimales H2-CO Verhältnis im Prinzip bei SN=2,05
# SN = 2,05 = (y_{H_2}-y_{CO_2})/(y_{CO}+y_{CO_2})
# y_CO
optimaler_faktor = ratio_h2_co2 / (
(1 - basis_y_co_plus_co2) * (2.05 * (1 + verhaeltnis_co_co2) + 1) +
ratio_h2_co2 * basis_y_co_plus_co2
)
factors_for_sn_0 = [3.835e-17, 0.375,
optimaler_faktor, 1.425,
1.754, 2.071, 3.289]
lines1 = [] # reinit lines
for factor in [3.835e-17, 0.375, 0.7,
1, optimaler_faktor, 1.425, 1.9, 2.071]:
y_0[:-2][indexes_co_co2] = np.array(
[verhaeltnis_co_co2 / (1 + verhaeltnis_co_co2), 1 / (1 + verhaeltnis_co_co2)]
) * basis_y_co_plus_co2 * factor
y_0[:-2][indexes_not_co_co2] = y_i0[indexes_not_co_co2] * (
1 - factor * basis_y_co_plus_co2
) / (
1 - basis_y_co_plus_co2
)
sn = (y_0[:-2][namen.index('H2')] - y_0[:-2][namen.index('CO2')]) / (
y_0[:-2][namen.index('CO2')] + y_0[:-2][namen.index('CO')]
)
p0_co2_co = sum(y_0[:-2][indexes_co_co2] * p0).item()
y_co_plus_co2 = sum(y_0[:-2][indexes_co_co2]).item()
# Keep mass of all other components m' = m (1-y1-y2)/(1-y1'-y2')
g_new = (1 - sum(y_i0[indexes_co_co2])) / (
1 - sum(y_0[:-2][indexes_co_co2])
) * sum(y_0[:-2] * mm / 1000.) / sum(y_i0 * mm / 1000.) * g
soln = odeint(df_dt, y_0, z_d_l_r)
y_i_soln = soln[:, :len(y_i0)]
p_soln = soln[:, -2]
t_soln = soln[:, -1]
mm_m_soln = np.sum(y_i_soln * mm * 1 / 1000., axis=1) # kg/mol
n_soln = g_new * n_t * 60 ** 2 * (np.pi / 4 * d_t ** 2) / mm_m_soln
m_soln = g_new * n_t * 60 ** 2 * (np.pi / 4 * d_t ** 2) # kg/h
m_i_soln = np.array([
n_soln * y_i_soln[:, i] * (mm[i] * 1 / 1000.)
for i in range(y_i_soln.shape[1])]).T
ums_soln = (y_co_plus_co2 * verhaeltnis_co_co2 / (1 + verhaeltnis_co_co2) * n_soln[0] - (
y_i_soln[:, namen.index('CO')]) * n_soln) / (
y_co_plus_co2 * verhaeltnis_co_co2 / (1 + verhaeltnis_co_co2) * n_soln[0]
)
ausb_soln = (y_i_soln[1:, namen.index('MeOH')] * n_soln[1:] -
y_i_soln[0, namen.index('MeOH')] * n_soln[0]) / (
ums_soln[1:] * y_co_plus_co2 * verhaeltnis_co_co2 / (1 + verhaeltnis_co_co2) * n_soln[0])
# Anzahl an Übertragungseinheiten (NTU)
mm_m_0 = sum(y_0[:-2] * mm) * 1 / 1000. # kg/mol
cp_m_0 = sum(y_0[:-2] * cp_ig_durch_r(t0) * 8.3145) # J/mol/K
cp_g_0 = cp_m_0 / mm_m_0 # J/kg/K
ntu = l_r * 1 / (g * cp_g_0) * 2 * u / (d_t / 2)
# m * m^2 s/kg * kg K /J * J/s/m^2/K *1/m = [dimensionslose Einheiten]
marker = np.random.choice(unfilled_markers)
color = color_samples[np.random.choice(len(color_samples))]
marker_color = color_samples[np.random.choice(len(color_samples))]
marker_fill = np.random.choice(list(lines.fillStyles))
line = ax.plot(t_soln, p_soln * y_i_soln[:, namen.index('CO')],
label='$p_{CO}+p_{CO_2}=' + '{:0.3g}'.format(p0_co2_co) +
'\quad SN= ' + '{:0.3g}'.format(sn) +
'\quad NTU= ' + '{:g}'.format(ntu) + '$',
marker=marker, markersize=7, color=color,
fillstyle=marker_fill, markerfacecolor=marker_color)[0]
lines1.append(line)
ax2.plot(t_soln, p_soln * y_i_soln[:, namen.index('CO2')],
label='$p_{CO}+p_{CO_2}=' + '{:0.3g}'.format(p0_co2_co) +
'\quad SN= ' + '{:0.3g}'.format(sn) +
'\quad NTU= ' + '{:g}'.format(ntu) + '$',
marker=marker, markersize=7, color=color,
fillstyle=marker_fill, markerfacecolor=marker_color)
ax3.plot(z_d_l_r, t_soln,
label='$p_{CO}+p_{CO_2}=' + '{:0.3g}'.format(p0_co2_co) +
'\quad SN= ' + '{:0.3g}'.format(sn) +
'\quad NTU= ' + '{:g}'.format(ntu) + '$',
marker=marker, markersize=7, color=color,
fillstyle=marker_fill, markerfacecolor=marker_color)
ax4.plot(z_d_l_r, ums_soln,
label='$p_{CO}+p_{CO_2}=' + '{:0.3g}'.format(p0_co2_co) +
'\quad SN= ' + '{:0.3g}'.format(sn) +
'\quad NTU= ' + '{:g}'.format(ntu) + '$',
marker=marker, markersize=7, color=color,
fillstyle=marker_fill, markerfacecolor=marker_color)
ax5.plot(z_d_l_r[1:], ausb_soln,
label='$p_{CO}+p_{CO_2}=' + '{:0.3g}'.format(p0_co2_co) +
'\quad SN= ' + '{:0.3g}'.format(sn) +
'\quad NTU= ' + '{:g}'.format(ntu) + '$',
marker=marker, markersize=7, color=color,
fillstyle=marker_fill, markerfacecolor=marker_color)
ax6.plot(z_d_l_r, m_i_soln[:, namen.index('MeOH')],
label='$p_{CO}+p_{CO_2}=' + '{:0.3g}'.format(p0_co2_co) +
'\quad SN= ' + '{:0.3g}'.format(sn) +
'\quad NTU= ' + '{:g}'.format(ntu) + '$',
marker=marker, markersize=7, color=color,
fillstyle=marker_fill, markerfacecolor=marker_color)
a = rho_b * sum(y_i0 * mm / 1000.) * \
y_i0[namen.index('H2')] * p0 * 1e5 # kgKat/mol kg/m^³ bar
b = -delta_h_r(t0)[0] * rho_b / (sum(y_i0 * cp_ig_durch_r(t0) * 8.3145) / sum(y_i0 * mm / 1000.)) * (
y_i0[namen.index('H2')] * p0) # kgKat/mol kg/m^3 K bar
c = 4 * u / (sum(y_i0 * cp_ig_durch_r(t0) * 8.3145) /
sum(y_i0 * mm / 1000.)) / d_t
t_m_reihe = np.linspace(t0, t_m_1, 50)
k_reihe = r_i(t_m_reihe, p0 * y_i0)[0] / (
p0**2 * (y_i0[namen.index('CO2')] * y_i0[namen.index('H2')]))
p_m_reihe = (t_m_reihe - t_r) / (b / c * (k_reihe))
ax2.plot(t_m_reihe, p_m_reihe, label='$T_m$', color='black', ls='--')
ax.plot([t_m_1, t_m_1], ax.get_ylim(), color='black', ls='--')
ax2.plot([t_m_1, t_m_1], ax2.get_ylim(), color='black', ls='--')
ax3.plot(ax2.get_xlim(), [t_m_1, t_m_1], color='black', ls='--')
fig.legend(lines1, [x.get_label() for x in lines1], 'upper left',
fontsize='xx-large')
plt.tight_layout(rect=[0, 0, 0.95, 0.85])
plt.show()
Daten aus der Grafik X(CO) gegen z für unterschiedliche SN exportieren:
In [47]:
import locale
import re
import os
# float-Zahl, ausschließlichfür Titeln
pattern=re.compile('.*SN=\s*(\-?[0-9]+\.[0-9]+(e\+?\-?[0-9]+)?)')
locale.setlocale(locale.LC_ALL, '')
[item.get_xdata() for item in ax4.get_lines()]
xdata = ax4.get_lines()[0].get_xdata()
header = [pattern.match(item.get_label()).groups()[0] for item in ax4.get_lines()]
with open('./logs/output_sn.txt', mode='w') as file:
file.write('z/L_R '+'X_{CO}(SN='+') X_{CO}(SN='.join(header)+')')
file.write('\n')
for i in range(len(ax4.get_lines()[0].get_xdata())):
xdata = ax4.get_lines()[0].get_xdata()[i]
temp_str = locale.format('%.20g', xdata)
for j in range(len(ax4.get_lines())):
ydata = ax4.get_lines()[j].get_ydata()[i]
temp_str=temp_str+' '+locale.format('%.20g', ydata)
file.write(temp_str+'\n')
file.close()
if os.path.exists(os.path.abspath('./logs/output_sn.txt')):
print('Datei expotiert: ' + './logs/output_sn.txt')
else:
print('Datei nicht exportiert.')
Das vollständige Verfahren besteht aus:
Um Den Auslegungsmassenstrom ($\dot m_1$) aus dem Makeup-Strom ($\dot m_0$) zu bestimmen ist daher eine Bilanz nötig, die die beiden mitsamt dem Reaktorprodukt verbindet:
Um den Mischer
$n_1 = n_0 + \frac{V}{F}n_2$
$m_{\text{Rückführung}}=\frac{V}{F}\cdot n_2\cdot\underbrace{\sum\limits_i^{\mathcal R}{y_{i, \text{Verf}}\cdot M_i}}_{M_V}$
Um den Reaktor
$m_1=m_2 \\ n_1\cdot M_1=n_2\cdot M_2$
Um den Mischer
$m_1=m_0+\frac{V}{F}\cdot n_2\cdot M_V$
Die Mengenstrombilanz zusammen mit der Massenstrombilanz werden bei jeder Iteration erfüllt.
Dazu liefern sie einen Ausdruck für die molare Masse des Eintritt-Stroms, und für den Massenstrom mit Rückführung:
$M_1=\frac{\frac{V}{F}\cdot n_2\cdot M_V+M_0}{\frac{V}{F}\frac{n_2}{n_0}+1}$
Massenstrom mit Rückführung:
$n_1\cdot M_1 = n_0\cdot M_0 + \frac{V}{F}\frac{m_1}{M_2}\cdot RV$
$\Rightarrow m_1 = m_0\cdot\frac{M_1\cdot M_2}{M_2 \cdot M_0-\frac{V}{F}\cdot RV\cdot M_1 \cdot M_0}$
Kennt man über die Betrachtung des Reaktors die Zusammensetzung des Produktes $M_2$; und über die Verflüssigung die Zusammensetzung der Rückführung $M_R$, so kann man exakt den Auslegungsmassenstrom $m_1$ bestimmen.
Bei jeder Rechnung weicht sich das Verflüssigungsverhältnis und die K-Werte von den unten angegebenen Werten kaum ab, wegen konstanter Ziel-Produktzusammensetzung. Es wird angenommen, dass $\frac{V}{F}=0,93-0,95$ bei einer Verflüssigungstemperatur von 60°C, und die Zusammensetzung der Dampf- bzw. Flüssugphase werden somit nach der Raschford-Rice Prozedur berechnet [11]:
$\begin{array}{ll} z_i &= \psi y_i + (1-\psi) x_i \\ &= \psi x_i K_i + (1-\psi) x_i \\ x_i &={{z_i}\over{1+\psi(K_i-1)}} \hspace{3cm} y_i &={{z_i K_i}\over{1+\psi(K_i-1)}}\\ \end{array}$
60°C, $\Psi=\frac{V}{F}=0,93 - 0,95$; K-Werte
$\begin{array}{ll}K_{CO}=&180.43464\\ K_{CO2}=&5.15534\\ K_{H2}=&244.60348\\ K_{H2O}=&0.01387\\ K_{MeOH}=&0.02512\\ K_{CH4}=&46.25476\\ K_{N2}=&196.12189\\ K_{EthOH}=&0.01978\\ K_{PrOH}=&0.01314\\ K_{METHF}=&0.39871\\ \end{array} $
In [48]:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import lines
from scipy.integrate import odeint
import locale
# Theorie und Aufstellung der Gleichungen:
# https://git.io/fdKBI
sn_param = 4.6 # 2.05 # dimensionslos
ntu_param = 3.09 # Parameter
verhaeltnis_co_co2 = 0.71 # Parameter
optimisieren_nach_param = ['n_t', 'l_r_n_t'][1]
namen = ['CO', 'CO2', 'H2', 'H2O', 'MeOH',
'CH4', 'N2', 'EthOH', 'PrOH', 'METHF']
# Katalysator
rho_b = 1190 # kg Kat/m^3 Feststoff
phi = 0.3 # m^3 Gas/m^3 Feststoff
m_kat = 1190 * (1 - 0.3) * np.pi / 4 * (
0.03)**2 * 7 # kg Kat (pro Rohr)
# Partikeldurchmesser wählen, damit
# Delta P gesamt=-3bar, denn dies ist der Parameter
d_p = 0.0054 / 0.0054 * 0.16232576224693065 # m Feststoff
# Reaktor
n_t = 1620 # Rohre
d_t = 0.04 # m Rohrdurchmesser
l_r = 7. # m Rohrlänge
# Betriebsbedingungen
t0 = 220 + 273.15 # K
p0 = 50 # bar
# Wärmetauschparameter
u = 118.44 # W/m^2/K
# Kühlmitteleigenschaften (VDI-WA)
t_r = (240 - 230) / (33.467 - 27.968) * (
29 - 33.467) + 240 + 273.15 # K
p_sat = 29 # bar
h_sat_l = (1037.5 - 990.21) / (33.467 - 27.968) * (
29 - 33.467) + 1037.5 # kJ/kg
h_sat_v = (2803.1 - 2803.0) / (33.467 - 27.968) * (
29 - 33.467) + 2803.1 # kJ/kg
delta_h_sat = (h_sat_v - h_sat_l)
# Zulaufbedingungen
n_i_0 = np.array([
1171.15, 1178.85, 3058.84,
100.65, 0, 0,
0, 0, 0,
0
]) / 60**2 * 1000 # mol/s
mm = np.array([
28.01, 44.01, 2.016,
18.015, 32.042, 16.043,
28.014, 46.069, 60.096,
60.053
], dtype=float) # g/mol
# Verflüssigung: K-Werte und V/F
v_f_flash = 0.9330904968264068
k_i_flash = np.array([
1.80434639e+02, 5.15534054e+00, 2.44603482e+02,
1.38677436e-02, 2.51200966e-02, 4.62547635e+01,
1.96121892e+02, 1.97755333e-02, 1.31435672e-02,
3.98712201e-01])
# Parameter der Wärmekapazität-Gleichung aus dem VDI-Wärmeatlas
cp_konstanten = np.zeros([len(namen), 7])
cp_konstanten_text = [
'407,9796 3,5028 2,8524 –2,3018 32,9055 –100,1815 106,1141',
'514,5073 3,4923 –0,9306 –6,0861 54,1586 –97,5157 70,9687',
'392,8422 2,4906 –3,6262 –1,9624 35,6197 –81,3691 62,6668',
'706,3032 5,1703 –6,0865 –6,6011 36,2723 –63,0965 46,2085',
'846,6321 5,7309 –4,8842 –12,8501 78,9997 –127,3725 82,7107',
'1530,8043 4,2038 –16,6150 –3,5668 43,0563 –86,5507 65,5986',
'432,2027 3,5160 2,8021 –4,1924 42,0153 –114,2500 111,1019',
'1165,8648 4,7021 9,7786 –1,1769 –135,7676 322,7596 –247,7349',
'506,0032 12,1539 0,0041 –36,1389 175,9466 –276,1927 171,3886',
'650,0705 5,0360 –0,4487 –13,5453 126,9502 –219,5695 151,4586'
]
for i in range(len(namen)):
cp_konstanten[i, :] = np.array(
cp_konstanten_text[i].replace(
'–', '-').replace(',', '.').split(' '),
dtype=float)
# Quelle: The properties of Gases and Liquids Poling Prausnitz
# Lennard-Jones Parameter Anhang B
# epsilon / kB
l_j_epsilon_d_k = np.array([
91.7, 195.2, 59.7,
809.1, 481.8, 148.6,
71.4, 362.6, 576.7,
469.8
]) # K
# sigma
l_j_sigma = np.array([
3.69, 3.941, 2.827,
2.641, 3.626, 3.758,
3.798, 4.53, 4.549,
4.936
]) # Angstrom
# T* = k T / epsilon
h_298 = np.array([
-110.53, -393.51, 0,
-241.81, -200.94, -74.52,
0, -234.95, -255.20,
-352.40
], dtype=float) * 1000. # J/mol
tc = np.array([
132.85, 304.12, 32.98,
647.14, 512.64, 190.56,
126.2, 513.92, 536.78,
487.2
]) # K
pc = np.array([
34.94, 73.74, 12.93,
220.64, 80.97, 45.99,
33.98, 61.48, 51.75,
60.0
]) # bar
omega_af = np.array([
0.045, 0.225, -0.217,
0.344, 0.565, 0.011,
0.037, 0.649, 0.629,
0.0
])
nuij = np.zeros([len(namen), 3])
# Hydrierung von CO2
nuij[[
namen.index('CO2'),
namen.index('H2'),
namen.index('MeOH'),
namen.index('H2O'),
namen.index('CO'),
namen.index('N2'),
], 0] = np.array([-1, -3, +1, +1, 0, 0], dtype=float)
# Hydrierung von CO
nuij[[
namen.index('CO2'),
namen.index('H2'),
namen.index('MeOH'),
namen.index('H2O'),
namen.index('CO'),
namen.index('N2'),
], 1] = np.array([0, -2, +1, 0, -1, 0], dtype=float)
# RWGS Reverse-Wassergasshiftreaktion (muss gleich sein als die Vorwärtsreaktion,
# wofür die Kinetik verfügbar ist)
nuij[[
namen.index('CO2'),
namen.index('H2'),
namen.index('MeOH'),
namen.index('H2O'),
namen.index('CO'),
namen.index('N2'),
], 2] = - np.array([+1, +1, 0, -1, -1, 0], dtype=float)
delta_h_r_298 = nuij.T.dot(h_298) # J/mol
def cp_ig_durch_r(t):
a, b, c, d, e, f, g = cp_konstanten.T
gamma_var = t / (a + t)
return b + (c - b) * gamma_var ** 2 * (
1 + (gamma_var - 1) * (
d + e * gamma_var + f * gamma_var ** 2 + g * gamma_var ** 3
)) # dimensionslos
def int_cp_durch_r_dt_minus_const(t):
a, b, c, d, e, f, g = cp_konstanten.T
return b * t + (c - b) * (
t - (d + e + f + g + 2) * a * np.log(a + t) +
-(2 * d + 3 * e + 4 * f + 5 * g + 1) * a**2 / (a + t) +
+(1 * d + 3 * e + 6 * f + 10 * g) * a**3 / 2 / (a + t)**2 +
-(1 * e + 4 * f + 10 * g) * a**4 / 3 / (a + t)**3 +
+(1 * f + 5 * g) * a**5 / 4 / (a + t)**4 +
- g * a**6 / 5 / (a + t)**5
)
def int_cp_durch_rt_dt_minus_const(t):
a, b, c, d, e, f, g = cp_konstanten.T
return b * np.log(t) + (c - b) * (
np.log(a + t) + (1 + d + e + f + g) * a / (a + t) +
-(d / 2 + e + 3 * f / 2 + 2 * g) * a**2 / (a + t)**2 +
+(e / 3 + f + 2 * g) * a**3 / (a + t)**3 +
-(f / 4 + g) * a**4 / (a + t)**4 +
+(g / 5) * a**5 / (a + t)**5
)
def mcph(x, t0_ref, t):
return sum(x * (
int_cp_durch_r_dt_minus_const(t) -
int_cp_durch_r_dt_minus_const(t0_ref))
) / sum(x) / (t - t0_ref)
def delta_h_r(t):
cp_m = (
int_cp_durch_r_dt_minus_const(t) -
int_cp_durch_r_dt_minus_const(298.15)
) / (t - 298.15) * 8.3145 # J/mol/K * K = J/mol
return delta_h_r_298 + nuij.T.dot(cp_m) * (t - 298.15)
def mu(t, y_i):
# T* = k T / epsilon
t_st = t / l_j_epsilon_d_k # J/K
# Stoßintegral (Bird Tabelle E.2)
omega_mu = 1.16145 / t_st**0.14874 + \
0.52487 / np.exp(0.77320 * t_st) + \
2.16178 / np.exp(2.43787 * t_st)
konst_1 = 5 / 16. * np.sqrt(
8.3145 * 1000 * 100 ** 2 / np.pi
) * 10**16 / 6.022e23 # g/cm/s
mu_i = konst_1 * np.sqrt(mm * t) / (
l_j_sigma**2 * omega_mu
) * 100 / 1000
# g/cm/s * 100cm/m * 1kg/1000g = kg/m/s = Pa s
phi_ab = 1 / np.sqrt(8) * (
1 + np.outer(mm, 1 / mm)) ** (-1 / 2.) * (
1 + np.outer(mu_i, 1 / mu_i) ** (1 / 2.) *
np.outer(1 / mm, mm) ** (1 / 4.)
) ** 2
mu_mix = np.sum(y_i * mu_i / phi_ab.dot(y_i))
return mu_mix # Pa s
# T-abhängige Parameter, dem Artikel nach
def k_t(t):
r = 8.3145 # Pa m^3/mol/K
# Geschwindigkeitskonstanten der 3 Reaktionen
k_1 = 10 ** (3066 / t - 10.592)
k_2 = 10 ** (5139 / t - 12.621)
k_3 = 10 ** (2073 / t - 2.029)
# Angepasste Parameter des kinetischen Modells
# A(i) exp(B(i)/RT)
a = np.array([
0.499, 6.62e-11, 3453.38, 1.07, 1.22e10
], dtype=float)
b = np.array([
17197, 124119, 0, 36696, -94765
], dtype=float)
k_h2 = a[0] ** 2 * np.exp(2 * b[0] / (r * t))
k_h2o = a[1] * np.exp(b[1] / (r * t))
k_h2o_d_k_8_k_9_k_h2 = a[2] * np.exp(b[2] / (r * t))
k5a_k_2_k_3_k_4_k_h2 = a[3] * np.exp(b[3] / (r * t))
k1_strich = a[4] * np.exp(b[4] / (r * t))
return np.array([
k_1, k_2, k_3,
k_h2, k_h2o, k_h2o_d_k_8_k_9_k_h2,
k5a_k_2_k_3_k_4_k_h2, k1_strich
])
def r_i(t, p_i):
p_co2 = p_i[namen.index('CO2')]
p_co = p_i[namen.index('CO')]
p_h2 = p_i[namen.index('H2')]
p_meoh = p_i[namen.index('MeOH')]
p_h2o = p_i[namen.index('H2O')]
[k_1, _, k_3,
k_h2, k_h2o, k_h2o_d_k_8_k_9_k_h2,
k5a_k_2_k_3_k_4_k_h2, k1_strich] = k_t(t)
r_meoh = k5a_k_2_k_3_k_4_k_h2 * p_co2 * p_h2 * (
1 - 1 / k_1 * p_h2o * p_meoh / (
p_h2 ** 3 * p_co2
)
) / (
1 + k_h2o_d_k_8_k_9_k_h2 * p_h2o / p_h2 +
np.sqrt(k_h2 * p_h2) + k_h2o * p_h2o
) ** 3
r_rwgs = k1_strich * p_co2 * (
1 - k_3 * p_h2o * p_co / (p_co2 * p_h2)
) / (
1 + k_h2o_d_k_8_k_9_k_h2 * p_h2o / p_h2 +
np.sqrt(k_h2 * p_h2) + k_h2o * p_h2o
) ** 1
return np.array([r_meoh, 0., r_rwgs])
def df_dt(y, _, g, d_p, l_r):
# FIXME: Normalize y_i on each iteration.
# Inerts should not change at all, but there is a 0.00056% increase already
# on the first time step. Test:
# g*60**2*n_t*(np.pi / 4 * d_t**2)/sum(y_i*mm/1000.)*y_i[6]*mm[6]/1000.
y_i = y[:-2]
p = y[-2]
t = y[-1]
mm_m = sum(y_i * mm) * 1 / 1000. # kg/mol
cp_m = sum(y_i * cp_ig_durch_r(t) * 8.3145) # J/mol/K
cp_g = cp_m / mm_m # J/kg/K
c_t = p / (8.3145 * 1e-5 * t) * 1
# bar * mol/bar/m^3/K*K = mol/m^3
p_i = y_i * p # bar
delta_h_r_t = delta_h_r(t)
mu_t_y = mu(t, y_i)
r_j = r_i(t, p_i) # mol/kg Kat/s
dyi_dz = l_r * rho_b * mm_m / g * (
nuij.dot(r_j) - y_i * sum(nuij.dot(r_j))
) # m * kgKat/m^3 * kg/mol * m^2 s/kg * mol/kgKat/s = dimlos
dp_dz = -l_r * 1 / 10**5 * g / (
c_t * mm_m * d_p
) * (1 - phi) / phi**3 * (
150 * (1 - phi) * mu_t_y / d_p + 1.75 * g
) # Pa * 1bar/(1e5 Pa) = bar
dt_dz = l_r / (
g * cp_g
) * (
2 * u / (d_t / 2) * (t_r - t) + rho_b * -delta_h_r_t.dot(r_j)
)
result = np.empty_like(y)
result[:-2] = dyi_dz
result[-2] = dp_dz
result[-1] = dt_dz
return result
def profile(n_i_1_ein, d_p_ein, optimisieren_nach='n_t'):
m_dot_i_ein = n_i_1_ein * mm / 1000. # kg/s
m_dot_ein = sum(m_dot_i_ein)
y_i1_ein = m_dot_i_ein / mm / sum(m_dot_i_ein / mm)
# Berechnung der Parameter
mm_m_1_ein = sum(y_i1_ein * mm) * 1 / 1000. # kg/mol
cp_m_1_ein = sum(y_i1_ein * cp_ig_durch_r(t1) * 8.3145) # J/mol/K
cp_g_1_ein = cp_m_1_ein / mm_m_1_ein # J/kg/K
if optimisieren_nach == 'n_t':
# Anzahl an Rohre anpassen
n_t_aus = ntu_param / (2 * np.pi * d_t / 2 * l_r * u) * (
m_dot_ein * cp_g_1_ein
)
l_r_aus = l_r
elif optimisieren_nach == 'l_r_n_t':
# Anzahl an Rohre und Länge anpassen
verhaeltnis = 1 / (l_r * n_t) * (
ntu_param * m_dot_ein * cp_g_1_ein / (
2 * np.pi * d_t / 2 * u)
)
l_r_aus = np.sqrt(verhaeltnis) * l_r
n_t_aus = np.sqrt(verhaeltnis) * n_t
g_neu = m_dot_ein / (np.pi / 4 * d_t ** 2) / n_t_aus # kg/m^2/s
y_0 = np.empty([len(namen) + 1 + 1])
y_0[:-2] = y_i1_ein
y_0[-2] = p1
y_0[-1] = t1
soln = odeint(lambda y, z0: df_dt(y, z0, g_neu, d_p_ein, l_r_aus),
y_0, z_d_l_r)
y_i_soln = soln[:, :len(y_i1_ein)]
p_soln = soln[:, -2]
t_soln = soln[:, -1]
mm_m_soln = np.sum(y_i_soln * mm * 1 / 1000., axis=1) # kg/mol
n_soln = g_neu * n_t_aus * (np.pi / 4 * d_t ** 2) / mm_m_soln
# kg/s/m^2 * m^2 / kg*mol = mol/s
n_i_soln = (y_i_soln.T * n_soln).T # mol/s
n_i_2 = n_i_soln[-1] # mol/s
n_0 = sum(n_i_0) # mol/s
n_2 = sum(n_i_2) # mol/s
t_2 = 60 + 273.15 # K
p_2 = p_soln[-1] # bar
z_i = y_i_soln[-1] # dimensionslos
#v_f_flash = v_f_flash
x_i_flash = z_i /(1 + v_f_flash*(k_i_flash-1))
y_i_flash = z_i * k_i_flash /(1 + v_f_flash*(k_i_flash-1))
y_co2_r = y_i_flash[namen.index('CO2')]
y_co_r = y_i_flash[namen.index('CO')]
y_h2_r = y_i_flash[namen.index('H2')]
y_h2_0 = n_i_0[namen.index('H2')] / n_0
y_co_0 = n_i_0[namen.index('CO')] / n_0
y_co2_0 = n_i_0[namen.index('CO2')] / n_0
n_i_1_aus = n_i_0 + v_f_flash * n_2 * y_i_flash
n_1 = sum(n_i_1_aus)
y_i_1 = n_i_1_aus / n_1
y_co_1 = y_i_1[namen.index('CO')]
y_co2_1 = y_i_1[namen.index('CO2')]
n_co_zus_aus = n_1 * (
verhaeltnis_co_co2 * y_co2_1 - y_co_1
)
n_h2_zus_aus = v_f_flash * n_2 * (
sn_param * (y_co2_r + y_co_r) +
y_co2_r - y_h2_r
) + n_0 * (
sn_param * (y_co2_0 + y_co_0) +
y_co2_0 - y_h2_0
) + sn_param * n_co_zus_aus
n_i_1_aus[namen.index('H2')] = \
n_i_1_aus[namen.index('H2')] + n_h2_zus_aus
n_i_1_aus[namen.index('CO')] = \
n_i_1_aus[namen.index('CO')] + n_co_zus_aus
# Partikeldurchmesser nach Parameter Delta_P=3bar optimisieren.
# Es wird direkt Fixpunkt-Iteration angewendet, nach der Form der Ergun Gl.
# D_p{n+1} = D_p{n} * int(1/D_p{n}*f(D_p{n}) dz) / -3bar
deltap_deltaz = (soln[-1][-2] - soln[0][-2]).item()
d_p_aus = deltap_deltaz * d_p_ein / -3.0
n_i_r_rueckf = v_f_flash * sum(n_i_2) * y_i_flash
n_i_0_vollst = n_i_0.copy()
n_i_0_vollst[namen.index('H2')] += n_h2_zus_aus
n_i_0_vollst[namen.index('CO')] += n_co_zus_aus
mm_0 = sum(n_i_0 / sum(n_i_0) * mm) / 1000. # kg/mol
mm_0_vollst = sum(n_i_0_vollst / sum(n_i_0_vollst) * mm) / 1000. # kg/mol
mm_2 = sum(n_i_2 / sum(n_i_2) * mm) / 1000. # kg/mol
mm_v = sum(n_i_r_rueckf / sum(n_i_r_rueckf) * mm) / 1000. # kg/mol
#mm_1 = (v_f_flash * sum(n_i_2) / sum(n_i_0_vollst) * mm_v + mm_0_vollst) / (
# v_f_flash * sum(n_i_2) / sum(n_i_0_vollst) + 1
#)
mm_1 = sum(n_i_1_aus / sum(n_i_1_aus) * mm) / 1000. # kg/mol
massenbilanz = mm_1 * sum(n_i_1_aus) - mm_v * sum(n_i_r_rueckf) \
- mm_0_vollst * sum(n_i_0_vollst)
mengenbilanz = sum(n_i_1_aus) - sum(n_i_r_rueckf) \
- sum(n_i_0_vollst)
aend = sum(n_i_1_aus - n_i_1_ein) / sum(n_i_1_ein) * 100
umsatz = 1 - n_i_2[namen.index('CO')] / n_i_1_aus[namen.index('CO')]
print('It. ' + '{:2d}'.format(i) + ' Änderung(n_1): ' +
'{:5.4f}'.format(np.sqrt(aend**2)) + '%\t' + ', Massenbilanz: ' +
'{:3.3e}'.format(np.sqrt(massenbilanz**2)) + '\t' + ', Mengenbilanz: ' +
'{:3.3e}'.format(np.sqrt(mengenbilanz**2)) + '\t' +
'Umsatz: ' +
'{:3.4f}'.format(umsatz) + '\t' + 'd_p= ' +
'{:3.4f}'.format(d_p_aus) + 'm')
return n_i_1_aus, y_i1_ein, n_t_aus, l_r_aus, d_p_aus, g_neu, \
cp_g_1_ein, m_dot_ein, t_soln, p_soln, \
n_soln, n_i_soln, y_i_soln, \
n_i_r_rueckf, n_h2_zus_aus, n_co_zus_aus
# Init.
m_dot_i_0 = n_i_0 * mm / 1000. # kg/s
m_dot_0 = sum(m_dot_i_0)
y_i0 = m_dot_i_0 / mm / sum(m_dot_i_0 / mm)
# Berechnung der Parameter
mm_m_0 = sum(y_i0 * mm) * 1 / 1000. # kg/mol
cp_m_0 = sum(y_i0 * cp_ig_durch_r(t0) * 8.3145) # J/mol/K
cp_g_0 = cp_m_0 / mm_m_0 # J/kg/K
n_i_1 = n_i_0 # mol/s
p1 = p0
t1 = t0
z_d_l_r = np.linspace(0, 1, 100)
dlr = 1 / (len(z_d_l_r) - 1) * l_r # m
for i in range(25):
n_i_1, y_i1, n_t, l_r, d_p, g, cp_g_1, \
m_dot, t_z, p_z, n_z, n_i_z, \
y_i_z, n_i_r, n_h2_zus, n_co_zus = profile(
n_i_1, d_p,
optimisieren_nach=optimisieren_nach_param
)
ntu = l_r * 1 / (g * cp_g_1) * 2 * u / (d_t / 2)
# m * m^2 s/kg * kg K /J * J/s/m^2/K *1/m = [dimensionslose Einheiten]
m_i_z = n_i_z * (mm * 1 / 1000.) * 60**2 # kg/h
# mol/h * g/mol * 1kg/1000g = 1/1000 kg/h
m_z = g * n_t * (np.pi / 4 * d_t ** 2) # kg/s
v_z = n_z * 8.3145 * 1e-5 * t_z / p_z
# mol/s * 8,3145Pa m^3/mol/K * 1e-5bar/Pa * K/bar = m^3/s
ums_z = 1 - n_i_z[:, namen.index('CO')] / \
n_i_z[0, namen.index('CO')]
m_km_z = u * (2 / (d_t / 2)) * (t_z - t_r) * (
np.pi / 4 * d_t ** 2) / delta_h_sat * n_t * 60 ** 2 / 1000.
# J/s/K/m^2 * 1/m * K * m^2 * kg/kJ * 60^2s/h * 1kJ/(1000J) = kg/h/m
# Tatsächliche stöchiometrische Zahl
sn = (y_i1[namen.index('H2')] - y_i1[namen.index('CO2')]) / (
y_i1[namen.index('CO2')] + y_i1[namen.index('CO')]
)
verhaeltnis_h2_co2 = y_i1[namen.index('H2')] / y_i1[namen.index('CO2')]
verhaeltnis_co_co2 = y_i1[namen.index('CO')] / y_i1[namen.index('CO2')]
# Energie-Analyse
t_m_1 = 1 / 2 * (36696 / 8.3145 - np.sqrt(36696 /
8.3145 * (36696 / 8.3145 - 4 * t_r)))
t_m_2 = 1 / 2 * (-94765 / 8.3145 - np.sqrt(-94765 /
8.3145 * (-94765 / 8.3145 - 4 * t_r)))
vars_1 = [
[r'\rho_b', rho_b, r'\frac{kg_{Kat}}{m^3_{Schüttung}}'],
['\phi', phi, r'\frac{m^3_{Gas}}{m^3_{Schüttung}}'],
['D_p', d_p, 'm'],
['D_t', d_t, 'm'],
['L_R', l_r, 'm'],
]
vars_2 = [
['n_T', n_t, ''],
['U', u, r'\frac{W}{m^2\cdot K}'],
['\dot m', m_dot, 'kg/s'],
['C_{p_g}', cp_g_1 / 1000., r'\frac{kJ}{kg\cdot K}'],
['NTU', ntu, ''],
]
vars_3 = [
['T_0', t0 - 273.15, '°C'],
['P_0', p0, 'bar'],
['T_r', t_r - 273.15, '°C_{Kühlmittel}'],
['P_{Sät}', p_sat, 'bar_{Kühlmittel}'],
['\Delta H_{Sät}', delta_h_sat, r'\frac{kJ}{kg_{Kühlmittel}}'],
['SN', sn, '']
]
text_1 = '\n'.join(['$' + ' = '.join([line[0], '{:g}'.format(line[1]) +
' ' + line[2]]) + '$'
for line in vars_1])
text_2 = '\n'.join(['$' + ' = '.join([line[0], '{:g}'.format(line[1]) +
' ' + line[2]]) + '$'
for line in vars_2])
text_3 = '\n'.join(['$' + ' = '.join([line[0], '{:g}'.format(line[1]) +
' ' + line[2]]) + '$'
for line in vars_3])
fig = plt.figure(1, figsize=(10,10))
font = {'size' : 14}
plt.rc('font', **font)
fig.suptitle('Lösung der Zusammensetzung ' +
'{:g}'.format(round(verhaeltnis_h2_co2.item(), 2)) +
':1:' +
'{:g}'.format(round(verhaeltnis_co_co2.item(), 2)) +
'(H2:CO2:CO)')
fig.text(0.05, 0.945, text_1, va='top', fontsize=14)
fig.text(0.33, 0.935, text_2, va='top', fontsize=14)
fig.text(0.66, 0.930, text_3, va='top', fontsize=14)
ax = plt.subplot2grid([2, 3], [0, 0])
ax.plot(z_d_l_r, v_z, label='$Idealgas$')
ax.set_ylabel(r'$\frac{\dot V}{m^3/s}$')
ax.set_xlabel('Reduzierte Position, $z/L_R$')
ax.legend(fontsize='xx-small')
ax2 = plt.subplot2grid([2, 3], [1, 0])
ax2.plot(z_d_l_r, m_km_z)
ax2.fill_between(z_d_l_r, 0, m_km_z, color='orange')
ax2.text(0.3, 1 / 2. * (m_km_z[0] + m_km_z[-1]),
'{:g}'.format(sum(m_km_z * dlr)) + 'kg/h')
ax2.set_ylabel(r'$\frac{\dot m_{Kuehlmittel}}{kg/h}$')
ax2.set_xlabel('Reduzierte Position, $z/L_R$')
ax3 = plt.subplot2grid([2, 3], [1, 1], colspan=2)
ax3.set_ylabel('Massenstrom / (kg/h)')
ax3.set_xlabel('Reduzierte Position, $z/L_R$')
for item in ['CO', 'H2O', 'MeOH', 'CO2']:
marker = np.random.choice(list(lines.lineMarkers.keys()))
index = namen.index(item)
ax3.plot(z_d_l_r, m_i_z[:, index], label=item,
marker=marker)
ax3.legend(loc=1)
ax4 = plt.subplot2grid([2, 3], [0, 1])
ax4_1 = ax4.twinx()
ax4_1.set_ylabel('Umsatz (CO)')
ax4.set_ylabel('Temperatur / °C')
ax4.set_xlabel('Reduzierte Position, $z/L_R$')
ax4.plot(z_d_l_r, t_z - 273.15, label='T / °C')
ax4_1.plot(z_d_l_r, ums_z, label='Umsatz (CO)', ls='--', color='gray')
ax5 = plt.subplot2grid([2, 3], [0, 2], colspan=2)
ax5.set_ylabel('Druck / bar')
ax5.set_xlabel('Reduzierte Position, $z/L_R$')
ax5.plot(z_d_l_r, p_z, label='p / bar')
plt.tight_layout(rect=[0, 0, 0.95, 0.75])
locale.setlocale(locale.LC_ALL, '')
print('')
print('=== GAS AUS REAKTOR ===')
print('SN= ' + str(sn))
print('NTU= ' + str(ntu))
print('y(CO)/y(CO2)= ' + str(verhaeltnis_co_co2))
print('\n'.join([
namen[i] + ': ' + locale.format('%.8g', x) + ' kg/h'
for i, x in enumerate(m_i_z[-1])
]))
print('T= ' + str(t_z[-1] - 273.15) + '°C')
print('P= ' + str(p_z[-1]) + 'bar')
print('V0= ' + str(v_z[0]) + 'm^3/s')
print('V= ' + str(v_z[-1]) + 'm^3/s')
print('Cpg= ' + str(cp_g_1) + 'J/kg/K')
print('Partikeldurchmesser für (DeltaP= ' +
'{:g}'.format(p_z[0] - p_z[-1]) + ' bar): ' +
'{:g}'.format(d_p) + ' m'
)
print('Kühlmittel: Gesättigtes H_2O(l) ' +
' bei ' + '{:g}'.format(p_sat) + ' bar' + '\n' +
'Verdampfungsenthalpie: ' + '{:g}'.format(delta_h_sat) +
'kJ/kg' + '\n' + 'Kühlmittelmassenstrom: ' +
'{:g}'.format(sum(m_km_z * dlr)) + 'kg/h')
print('')
print('======' * 3)
print('=== MAKEUP-STROM ===')
print('\n'.join([
namen[i] + ': ' + locale.format('%.8g', x) + ' kg/h'
for i, x in enumerate(n_i_0 * mm / 1000. * 60**2)
]))
print('')
print('======' * 3)
print('=== RÜCKLAUFSTROM ===')
print('\n'.join([
namen[i] + ': ' + locale.format('%.8g', x) + ' kg/h'
for i, x in enumerate(n_i_r * mm / 1000. * 60**2)
]))
print('')
print('======' * 3)
print('=== ERFORDERLICHE CO UND H2 STRÖME, UM SN UND CO/CO2 ANZUPASSEN ===')
print('H2: ' +
locale.format('%.8g', n_h2_zus.item() * mm[namen.index('H2')] /
1000. * 60 ** 2) + ' kg/h')
print('CO: ' +
locale.format('%.8g', n_co_zus.item() * mm[namen.index('CO')] /
1000. * 60**2) + ' kg/h')
print('auf kmol/h')
print('H2: ' +
locale.format('%.8g', n_h2_zus.item() / 1000. * 60 ** 2) + ' kmol/h')
print('CO: ' +
locale.format('%.8g', n_co_zus.item() / 1000. * 60 ** 2) + ' kmol/h')
print('')
print('======' * 3)
print('=== REAKTOR ZULAUFSTROM ===')
print('\n'.join([
namen[i] + ': ' + locale.format('%.8g', x) + ' kg/h'
for i, x in enumerate(n_i_1 * mm / 1000. * 60**2)
]))
plt.show()
Für den Rücklauf wird die isotherme Verflüssigungs anhand einer kubischen Zustandsgleichung berechnet.
Vorgehensweise bei der Rechnug: Siehe Anhang B unten.
Numerische Lösung: Link
Nach kinetischer Gastheorie (Bird 1.4)
$\require{cancel}$ $\begin{array}{ll} \overline c &=\sqrt{\frac{8 R T}{\pi M}} \quad \quad \text{Mittlere Geschwindigkeit}\\ \lambda &= \frac{k_B T}{\sqrt{2}\sigma p}=\frac{k_B T}{\sqrt{2}\cdot \pi d^2 \cdot p}\quad \quad \text{Mittlere freie Weglänge}\\ \mathcal N &= \frac{n\cdot N_A}{V} \quad\quad\text{...Teilchendichte}\\ m &= \frac{M}{N_A} \quad\quad\text{...Teilchenmasse}\\ \mu&=\frac{1}{3}\lambda\overline c m \mathcal N=\frac{\overline c m}{3\sqrt{2}\cdot\sigma N_A}\\ &=\frac{1}{3}\cdot\frac{k_B T}{\sqrt{2}\cdot \pi d^2\cdot (n R T/V))}\cdot\sqrt{\frac{8 R T}{\pi M}}\cdot m\cdot \frac{n\cdot N_A}{V}\\ &=\frac{2}{3}\cdot\frac{1}{\pi d^2}\cdot\sqrt{\frac{k_B N_A T}{\pi m N_A}}\cdot\cancelto{1}{\frac{k_B\cdot N_A}{R}}\cdot m\\ &=\frac{2}{3}\cdot\frac{\sqrt{m k_B T /\pi}}{\pi d^2}=\frac{2}{3}\cdot\frac{\sqrt{M R T/\pi}}{N_A\pi d^2} \quad \text{Maxwell 1860 (nicht so gut)}\\ \mu&=\frac{5}{16}\cdot\frac{\sqrt{\pi m k_B T}}{\pi \sigma^2} \quad \text{Chapman Enskog XIX Jh. (besser) Bird 1.4-9} \\ \end{array}$
Als Funktion der molaren Masse ($m=M/N_A$), und Trennung der Faktoren in dimensionslosem Teil und dimensionsbedingtem Teil:
$\begin{array}{ll} \mu &=\frac{5}{16}\cdot\frac{\sqrt{\pi M R T}}{N_A\pi \sigma^2}= \underbrace{\frac{5}{16}\cdot\frac{\sqrt{8,3145\frac{J}{mol K}\cdot\frac{1 kg m^2/s^2}{J}\cdot\frac{1000 g}{kg}\cdot\left(\frac{100cm}{m}\right)^2/\pi}}{\pi\cdot 6,0221\cdot10^{23} \frac{1}{mol}\cdot \mathring{A}^2\cdot\left(\frac{1m}{10^{10}\mathring A}\right)^2\cdot\left(\frac{100cm}{m}\right)^2}\cdot}_{\frac{5}{16}\sqrt{8.3145\cdot1000\cdot100^2/\pi}\cdot\frac{10^{16}}{6,0221\cdot 10^{23}}\cdot\frac{g}{cm s}} \frac{\sqrt{\left[\frac{M}{g/mol}\right]\cdot\left[\frac{T}{K}\right]}}{\left(\frac{\sigma}{\mathring{A}}\right)^2\Omega_\mu}\\ &=2,6696\cdot10^{-5}\frac{g}{cm s}\cdot\frac{\sqrt{\left[\frac{M}{g/mol}\right]\cdot\left[\frac{T}{K}\right]}}{\left(\frac{\sigma}{\mathring{A}}\right)^2\Omega_\mu}\\ \end{array}$
Funktion aus VDI-Wärmeatlas (Kapitel D3, Gleichung 10) liefert die spezifische Wärmekapazität idealer Gase bei T mit den Koeffizienten aus D3.1: $\frac{C_P^{id}(T)}{R}$. Stammfunktionen für die Berechnung von $\int\limits_{T_0}^T{\frac{C_P^{id}(T)}{R}dT}$ sind schon in [7] angegeben.
Definition des gemittelten Wertes für die Enthalpie-Berechnung, nach Smith-Van Ness-Ansatz:
$\left<\frac{C_P^{id}(T)}{R} \right>_{T_0}^T=\frac{\int\limits_{T_0}^T{\frac{C_P^{id}(T)}{R}dT}}{T-T_0}$
Definition des gemittelten Wertes für die Entropie-Berechnung, nach Smith-Van Ness-Ansatz:
$\left<\frac{C_P^{id}(T)}{R} \right>_{T_0}^T=\frac{\int\limits_{T_0}^T{\frac{C_P^{id}(T)}{R T}dT}}{ln(T/T_0)}$
Bei 298,15K setzt man die Bindungsenthalpien der Jeweiligen Komponenten (aus Gmehling-Kolbe) in der folgenden Gleichung ein:
$\Delta_R H^{\circ}_{j}(298,15K)=\sum\limits_i^{\mathcal K}\nu_{ji}\Delta_B H^{\circ}_{i}(298,15K)$
Bei T berechnet man die gemittelten $C_P$ Werte zwischen T und 298.15K, und führt die Enthalpieunterschiede der jeweiligen Komponenten auf die Temperatur T mittels der folgenden Gleichung:
$\Delta_R H^{\circ}_{j}(T)=\Delta_R H^{\circ}_{j}(298,15K)+\sum\limits_i^{\mathcal K}\nu_{ji}\left<\frac{C_P^{id}(T)}{R} \right>_{T_0}^T \cdot (T-T_0)$
Rachford - Rice [11]
Variablen: (N) Molbrüche in der Flüssigkeit $\{\dot x_i\}$, (N) Molbrüche im Gas $\{y_i\}$, (1) Dampf/Flüssigkeit Verhältnis
Gleichungen: (N) Stoffbilanzen, (N) Isotherme Gleichgewichtsbedingungen, (1) Zwangsbedingung
$\dot F z_i = \dot V y_i + \dot L x_i $
$(\hat{\phi_i}^V P) y_i^V=(\hat{\phi_i}^L P) x_i^L \hspace{2cm} K_i\equiv{{y_i}\over{x_i}}={{\hat{\phi_i}^L}\over{\hat{\phi_i}^V}}$
$0 = \sum\limits_i y_i - \sum\limits_i x_i$
System bei T und P zu einer Funktion einer Veränderlichen reduziert: $\psi\equiv(V/F)$
$\begin{array}{ll} z_i &= \psi y_i + (1-\psi) x_i \\ &= \psi x_i K_i + (1-\psi) x_i \\ x_i &={{z_i}\over{1+\psi(K_i-1)}} \hspace{3cm} y_i &={{z_i K_i}\over{1+\psi(K_i-1)}}\\ \Rightarrow 0 &= \sum\limits_i y_i - \sum\limits_i x_i = \sum\limits_i {{z_i(K_i-1)}\over{1+\psi(K_i-1)}}\\ \end{array}$
(Ojektivfunktion)
Lösung nach Rachford Rice-Algorithmus [11]