1) Oberflächenreaktion A+2B==>>2P 2) Oberflächenreaktion A+3B==>> P 3) Freie Ausströmung

ttp_01.py

from numpy import array, linspace, exp, zeros, diff
from scipy.integrate import odeint, ode
from matplotlib import pyplot as plt
import ctypes  # Needed to set the app icon correctly

r1, r2, r3 = [0.33, 1, -0.33]
d21 = d12 = d23 = d32 = 1e-5 # m^2/s
d13 = d31 = 1/10 * d12
y0 = [0.3, 0.7, 0]

zeta = linspace(0,2.5e-6, 40)

def dy_dzeta(y, t):
    y1, y2, y3 = y
    return array([
        1/d12*(y1*r2 - y2*r1) + 1/d13*(y1*r3 - y3*r1),
        1/d21*(y2*r1 - y1*r2) + 1/d23*(y2*r3 - y3*r2),
        1/d31*(y3*r1 - y1*r3) + 1/d32*(y3*r2 - y2*r3),
    ])

soln = odeint(dy_dzeta, y0, zeta)
lines = plt.plot(zeta, soln)
lines[0].set_label(r'$y_{C_6H_6}$')
lines[1].set_label(r'$y_{H_2}$')
lines[2].set_label(r'$y_{C_{6}H_{12}}$')
plt.xlabel(r'$\zeta \quad [-]$')
plt.ylabel(r'$y \quad [-]$')
plt.title('Skript 16.1.1')
plt.legend()

plt.figure()

#20140310
r1, r2, r3 = [-1, 2, 0]
d21 = d12 = d23 = d32 = 1e-5 # m^2/s
d13 = d31 = 1/1 * d12
y0 = [0, 2/3, 1/3]

zeta = linspace(0,2.5e-6*2, 40)

def dy_dzeta(t, y):
    y1, y2, y3 = y
    return array([
        1/d12*(y1*r2 - y2*r1) + 1/d13*(y1*r3 - y3*r1),
        1/d21*(y2*r1 - y1*r2) + 1/d23*(y2*r3 - y3*r2),
        1/d31*(y3*r1 - y1*r3) + 1/d32*(y3*r2 - y2*r3),
    ])

soln = ode(dy_dzeta).set_initial_value(y0)
i=0
dzeta = diff(zeta)[0]
y = zeros([len(zeta), 3])
y[0] = y0
while soln.successful() and i < zeta.size-1 and (y>=0).all():
    i+=1
    soln.integrate(soln.t+dzeta)
    y[i] = soln.y

lines = plt.plot(zeta[:i+1], y[:i+1], 'o', fillstyle='none')
plt.plot(zeta[:i+1], r1-exp(zeta[:i+1]/d12)*(r1-y0[0]), '-', color=lines[0].get_color(), 
    label=r'$\widetilde{y}_{1,S}=\dot{r}_1-exp\left(\frac{\zeta}{\delta} \right)(\dot{r}_1-\widetilde{y}_{1,0})$')
plt.plot(zeta[:i+1], r2-exp(zeta[:i+1]/d12)*(r2-y0[1]), '-', color=lines[1].get_color(),
    label=r'$\widetilde{y}_{2,S}=\dot{r}_2-exp\left(\frac{\zeta}{\delta} \right)(\dot{r}_2-\widetilde{y}_{2,0})$')
plt.plot(zeta[:i+1], r3-exp(zeta[:i+1]/d12)*(r3-y0[2]), '-', color=lines[2].get_color(),
    label=r'$\widetilde{y}_{3,S}=\dot{r}_3-exp\left(\frac{\zeta}{\delta} \right)(\dot{r}_3-\widetilde{y}_{3,0})$')
lines[0].set_label(r'$y_{O_2}$')
lines[1].set_label(r'$y_{CO}$')
lines[2].set_label(r'$y_{N_2}$')
plt.xlabel(r'$\zeta \quad [-]$')
plt.ylabel(r'$y \quad [-]$')
plt.title(
    r'$\delta=10^{-5}\frac{m^2}{s} \quad'+
    r'\dot{r}_1=-1, \dot{r}_2=+2, \dot{r}_3=0 \quad'+
    r'\widetilde{y}_{O_2,0}=0, '+
    r'\widetilde{y}_{CO,0}=2/3, '+
    r'\widetilde{y}_{N_2,0}=1/3$'
    )
plt.suptitle('20140310')
plt.legend()

plt.figure()

t = linspace(0,4,30)
plt.plot(t, exp(-t), '--', label='exp(-t)')
plt.plot(t, (1-t/2)**2, '--', label=r'$(1-t/2)^2$')

soln1 = odeint(lambda m, t: -m, 1, t)
soln2 = ode(lambda t, m: -m**(1/2.))
soln2.set_initial_value(1, 0).set_integrator('vode')
dt = (t[-1]-t[-2])
y = zeros(t.size)
y[0] = 1
i = 0
while soln2.successful() and i < t.size-1:
    i+=1
    soln2.integrate(soln2.t+dt)
    y[i] = soln2.y

plt.plot(t, soln1, 'o', fillstyle='none', label=r'$\frac{dM}{dt}=-c_{lam} M$')
plt.plot(t, y, 'o', fillstyle='left', label=r'$\frac{dM}{dt}=-c_{turb} M^{1/2}$')
plt.xlabel(r'$t \quad [s] $')
plt.ylabel(r'$\frac{M}{M_0} \quad [-]$')
plt.legend()

# following 2 lines for setting app icon correctly
myappid = u'ssv.ttp.01.0'  # arbitrary string
ctypes.windll.shell32.SetCurrentProcessExplicitAppUserModelID(myappid)
# ended lines for setting app icon correctly
plt.show()