Beispiel 4.2: Berechnung von Gleichgewichtskonstanten¶

(Dieses Beispiel finden Sie im Lehrbuch auf Seite 42)

In diesem Berechnungsbeispiel sollen exemplarisch die thermodynamische Gleichgewichtskonstante K$°$ sowie die spezielle Gleichgewichtskonstante K$_{x}$ für die Ammoniaksynthese berechnet werden. Die erforderlichen Gleichungen (4.1) bis (4.15) sind auf den Seiten 37 bis 43 beschrieben.

Laden der benötigten Pakete¶

Innerhalb dieses Beispiels werden die Pakete numpy und matplotlib.pyplot benötigt. Während numpy wichtige Rechenoperationen und Funktionalitäten bereitstellt, werden mit Hilfe des Pakets matplotlib.pyplot die Ergebnisse grafisch dargestellt.

In [1]:

import numpy as np
import matplotlib.pyplot as plt

Definition der Arrays¶

Zuerst müssen die jeweiligen Arrays initialisiert werden. Die Komponenten Stickstoff (N₂), Wasserstoff (H₂) und Ammoniak (NH₃) werden mit den Indizes 0, 1 und 2 beschrieben.

In [2]:

comps = np.linspace(0,2,3) # 0: N_2, 1: H_2, 2: NH_3
nu = np.array((comps))
cp = np.array((comps))
H_f = np.array((comps))
S_g = np.array((comps))

Definition der Konstanten¶

Stöchiometrische Koeffizienten¶

Die stöchiometrischen Koeffizienten $\nu_i$ können in der Reaktionsgleichung (Tabelle 4.1, Seite 38) abgelesen werden:

N₂ + 3 H₂ ⇌ 2 NH₃

Während die stöchiometrischen Koeffizienten der Edukte ein negatives Vorzeichen erhalten, ist das Vorzeichen der Produkte positiv.

In [3]:

nu[0] = -1 # N_2
nu[1] = -3 # H_2
nu[2] = 2 # NH_3

Thermodynamische Daten¶

Außerdem werden für die Berechnung einige thermodynamische Daten benötigt, welche folgender Quelle entnommen wurden: https://webbook.nist.gov/chemistry (2017)

Wärmekapazität $c_\mathrm{p}$¶

In [4]:

cp[0] = 29.1 # J / (mol * K)
cp[1] = 28.8 # J / (mol * K)
cp[2] = 35.1 # J / (mol * K)

Bildungsenthalpie $H_\mathrm{f}$¶

In [5]:

H_f[0] = 0 # J / (mol)
H_f[1] = 0 # J / (mol)
H_f[2] = -45.9e3 # J / (mol)

Entropie $S$¶

In [6]:

S_g[0] = 191.61 # J / (mol * K)
S_g[1] = 130.68 # J / (mol * K)
S_g[2] =  192.77 # J / (mol * K)

Referenztemperatur T$°$, Referenzdruck p$°$ sowie universelle Gaskonstante R¶

In [7]:

T_S = 298.15 # K
p_S = 1e5 # Pa
R = 8.314 # J (mol * K)

Berechnung der Reaktionsdaten¶

Die Wärmekapazität $\Delta_\mathrm{R}c_\mathrm{p}$ (Gleichung 4.5), die Reaktionsenthalpie $\Delta_\mathrm{R}H$ (Gleichung 4.1) und die Reaktionsentropie $\Delta_\mathrm{R}S$ (Gleichung 4.11) unter Standardbedingungen sind das Skalarprodukt (numpy-Funktion np.dot()) aus den thermodynamischen Daten und den stöchiometrischen Koeffizienten.

In [8]:

Delta_R_cp_S = np.dot(nu,cp) # J / (mol * K)
Delta_R_H_S = np.dot(nu,H_f) # J / (mol * K)
Delta_R_S_S = np.dot(nu,S_g) # J / (mol * K)

Definition der temperaturabhängigen Funktionen¶

Im weiteren Verlauf müssen die temperaturabhängigen Funktionen für die Reaktionsenthalpie $\Delta_\mathrm{R}H(T)$ (Gleichung 4.4), die Reaktionsentropie $\Delta_\mathrm{R}S(T)$ (Gleichung 4.10) sowie für die freie Reaktionsenthalpie $\Delta_\mathrm{R}G(T)$ (Gleichung 4.9) definiert werden.

In [9]:

def Delta_R_H(T):
    return (Delta_R_H_S + Delta_R_cp_S * (T - T_S))

def Delta_R_S(T):
    return (Delta_R_S_S + Delta_R_cp_S * np.log(T/T_S))

def Delta_R_G(T):
    return( Delta_R_H(T) - T * Delta_R_S(T))

Definition der Gleichgewichtskonstanten¶

Die thermodynamische Gleichgewichtskonstante K$°$ berechnet sich nach Gleichung (4.12).

In [10]:

def K_S(T):
    return( np.exp( -Delta_R_G(T) / (R * T) ) )

Für die spezielle Gleichgewichtskonstante K$_{x}$ kann Gleichung (4.14) herangezogen werden.

In [11]:

def K_x(T,p):
    return (K_S(T) * (p/p_S)**2)

Darstellung der Ergebnisse¶

Für die Darstellung der Ergebnisse muss der zu betrachtende Temperaturbereich $T_\mathrm{range}$ angegeben werden. In diesem Fall ist dies der Bereich von 400 °C bis 500 °C in Schritten von 1 K.

In [12]:

T_range = np.linspace(273.15 + 400, 273.15 + 500, 101) # K

Anschließend werden die Funktionen über die Temperatur mit Hilfe des Pakets matplotlib.pyplot dargestellt.

In [13]:

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12,5))

ax1.plot(T_range-273.15, K_S(T_range)*1e4, color='black')
ax1.set_xlabel("T / °C")
ax1.set_ylabel("K$°$ $\cdot 10^{4}$ / 1")
ax1.set_xlim(400,500)
ax1.set_ylim(0,1.5)

ax2.plot(T_range-273.15, K_x(T_range, 100e5), color='black', linestyle='dotted', label="100 bar")
ax2.plot(T_range-273.15, K_x(T_range, 200e5), color='black', linestyle='dashed', label="200 bar")
ax2.plot(T_range-273.15, K_x(T_range, 300e5), color='black', linestyle='solid', label="300 bar")
ax2.set_xlabel("T / °C")
ax2.set_ylabel("K$_{x}$ / 1")
ax2.set_xlim(400,500)
ax2.set_ylim(0,14)
ax2.legend(loc='best')

plt.show()

Die Berechnungsergebnisse sind ebenfalls im Lehrbuch in Abbildung 4.1 auf Seite 43 dargestellt.