/
integrated.py
243 lines (207 loc) · 7.08 KB
/
integrated.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
"""
This module collects a few analytic solutions of initial value problems (IVPs) in
chemical kinetics (i.e. integrated rate expressions in closed form).
The expressions are useful in e.g. regression or for comparison with numerical
solution of the corresponding ODE system.
"""
from __future__ import (absolute_import, division, print_function)
from .._util import get_backend
def dimerization_irrev(t, kf, initial_C, P0=1, t0=0):
return 1/(1/initial_C + 2*kf*(t-t0))
def pseudo_irrev(t, kf, prod, major, minor, backend=None):
""" Analytic product transient of a irreversible pseudo first order reaction.
Product concentration vs time from pseudo-first order irreversible kinetics.
Parameters
----------
t : float, Symbol or array_like
Time.
kf : number or Symbol
Forward (bimolecular) rate constant.
kb : number or Symbol
Backward (unimolecular) rate constant.
prod : number or Symbol
Initial concentration of the complex.
major : number or Symbol
Initial concentration of the more abundant reactant.
minor : number or Symbol
Initial concentration of the less abundant reactant.
backend : module or str
Default is 'numpy', can also be e.g. ``sympy``.
"""
be = get_backend(backend)
return prod + minor*(1 - be.exp(-major*kf*t))
pseudo_irrev.name = 'Pseudo first order irreversible'
def pseudo_rev(t, kf, kb, prod, major, minor, backend=None):
""" Analytic product transient of a reversible pseudo first order reaction.
Product concentration vs time from pseudo-first order reversible kinetics.
Parameters
----------
t : float, Symbol or array_like
Time.
kf : number or Symbol
Forward (bimolecular) rate constant.
kb : number or Symbol
Backward (unimolecular) rate constant.
prod : number or Symbol
Initial concentration of the complex.
major : number or Symbol
Initial concentration of the more abundant reactant.
minor : number or Symbol
Initial concentration of the less abundant reactant.
backend : module or str
Default is 'numpy', can also be e.g. ``sympy``.
"""
be = get_backend(backend)
return (
-kb*prod + kf*major*minor + (kb*prod - kf*major*minor) *
be.exp(-t*(kb + kf*major))
)/(kb + kf*major)
pseudo_rev.name = 'Pseudo first order reversible'
def binary_irrev(t, kf, prod, major, minor, backend=None):
""" Analytic product transient of a irreversible 2-to-1 reaction.
Product concentration vs time from second order irreversible kinetics.
Parameters
----------
t : float, Symbol or array_like
kf : number or Symbol
Forward (bimolecular) rate constant.
prod : number or Symbol
Initial concentration of the complex.
major : number or Symbol
Initial concentration of the more abundant reactant.
minor : number or Symbol
Initial concentration of the less abundant reactant.
backend : module or str
Default is 'numpy', can also be e.g. ``sympy``.
"""
be = get_backend(backend)
return prod + major*(1 - be.exp(-kf*(major-minor)*t))/(
major/minor - be.exp(-kf*t*(major-minor)))
binary_irrev.name = 'Second order irreversible'
def binary_rev(t, kf, kb, prod, major, minor, backend=None):
""" Analytic product transient of a reversible 2-to-1 reaction.
Product concentration vs time from second order reversible kinetics.
Parameters
----------
t : float, Symbol or array_like
Time.
kf : number or Symbol
Forward (bimolecular) rate constant.
kb : number or Symbol
Backward (unimolecular) rate constant.
prod : number or Symbol
Initial concentration of the complex.
major : number or Symbol
Initial concentration of the more abundant reactant.
minor : number or Symbol
Initial concentration of the less abundant reactant.
backend : module or str
Default is 'numpy', can also be e.g. ``sympy``.
"""
# see _integrated.ipynb for derivation
be = get_backend(backend)
X, Y, Z = prod, major, minor
x0 = Y*kf
x1 = Z*kf
x2 = 2*X*kf
x3 = -kb - x0 - x1
x4 = -x2 + x3
x5 = be.sqrt(-4*kf*(X**2*kf + X*x0 + X*x1 + Z*x0) + x4**2)
x6 = kb + x0 + x1 + x5
x7 = (x3 + x5)*be.exp(-t*x5)
x8 = x3 - x5
return (x4*x8 + x5*x8 + x7*(x2 + x6))/(2*kf*(x6 + x7))
binary_rev.name = 'Second order reversible'
def unary_irrev_cstr(t, k, r, p, fr, fp, fv, backend=None):
""" Analytic solution for ``A -> B`` in a CSTR.
Analytic solution for a first order process in a continuously
stirred tank reactor (CSTR).
Parameters
----------
t : array_like
k : float_like
Rate constant
r : float_like
Initial concentration of reactant.
p : float_like
Initial concentration of product.
fr : float_like
Concentration of reactant in feed.
fp : float_like
Concentration of product in feed.
fv : float_like
Feed rate / tank volume ratio.
backend : module or str
Default is 'numpy', can also be e.g. ``sympy``.
Returns
-------
length-2 tuple
concentrations of reactant and product
"""
# See _kinetics_cstr.ipynb
be = get_backend(backend)
x0 = fr*fv
x1 = fv + k
x2 = 1/x1
x3 = fv*r + k*r - x0
x4 = fr*k
x5 = be.exp(-fv*t)
return (
x0*x2 + x2*x3*be.exp(-t*x1),
-x2*x3*x5*(-1 + be.exp(-k*t)) + x2*x5*(-fp*fv - fp*k + fv*p + k*p - x4) + x2*(fp*x1 + x4)
)
def binary_irrev_cstr(t, k, r, p, fr, fp, fv, n=1, backend=None):
""" Analytic solution for ``2 A -> n B`` in a CSTR.
Parameters
----------
t : array_like
k : float_like
Rate constant
r : float_like
Initial concentration of reactant.
p : float_like
Initial concentration of product.
fr : float_like
Concentration of reactant in feed.
fp : float_like
Concentration of product in feed.
fv : float_like
Feed rate / tank volume ratio.
n : int
backend : module or str
Default is 'numpy', can also be e.g. ``sympy``.
Returns
-------
length-2 tuple
concentrations of reactant and product
"""
# Mathematica source:
# FortranForm[
# DSolve[{Derivative[1][A][t] == a0 f - f A[t] - 2 k A[t]^2,
# Derivative[1][B][t] == b0 f + n*k A[t]^2 - f B[t], A[0] == x,
# B[0] == y}, {A[t], B[t]}, {t}]]
# Post processed using sympy's cse function.
# (see _derive_analytic_cstr_bireac.ipynb)
be = get_backend(backend)
atanh = getattr(be, 'atanh', be.arctanh)
three = 3*be.cos(0)
x0 = 1/k
x1 = be.sqrt(fv)
x2 = 8*k
x3 = fr*x2
x4 = be.sqrt(fv + x3)
x5 = x1*x4
x6 = x1*x4/2
x7 = atanh((-fv**(three/2)*x4 - 4*k*r*x5)/(fv**2 + fv*x3))
x8 = fv*t
x9 = fp*x2
x10 = 4*k*n
x11 = fr*x10
x12 = be.exp(x8)
x13 = n*x12
return (
x0*(-fv + x5*be.tanh(t*x6 - x7))/4,
x0*(fv*x13 + 8*k*p + r*x10 - x1*x13*x4*be.tanh(
x6*(t - 2*x7/(x1*x4))
) + x11*x12 - x11 + x12*x9 - x9)*be.exp(-x8)/8
)