Integrating equations in python

| categories: | tags: | View Comments

A common need in engineering calculations is to integrate an equation over some range to determine the total change. For example, say we know the volumetric flow changes with time according to $$d\nu/dt = \alpha t$$, where $$\alpha = 1$$ L/min and we want to know how much liquid flows into a tank over 10 minutes if the volumetric flowrate is $$\nu_0 = 5$$ L/min at $$t=0$$. The answer to that question is the value of this integral: $$V = \int_0^{10} \nu_0 + \alpha t dt$$.

import scipy

nu0 = 5     # L/min
alpha = 1.0 # L/min
def integrand(t):
return nu0 + alpha * t

t0 = 0.0
tfinal = 10.0
V, estimated_error = quad(integrand, t0, tfinal)
print('{0:1.2f} L flowed into the tank over 10 minutes'.format(V))

100.00 L flowed into the tank over 10 minutes


That is all there is too it!

org-mode source

Integrating a batch reactor design equation

| categories: integration | tags: reaction engineering | View Comments

For a constant volume batch reactor where $$A \rightarrow B$$ at a rate of $$-r_A = k C_A^2$$, we derive the following design equation for the length of time required to achieve a particular level of conversion :

$$t(X) = \frac{1}{k C_{A0}} \int_{X=0}^X \frac{dX}{(1-X)^2}$$

if $$k = 10^{-3}$$ L/mol/s and $$C_{A0}$$ = 1 mol/L, estimate the time to achieve 90% conversion.

We could analytically solve the integral and evaluate it, but instead we will numerically evaluate it using scipy.integrate.quad. This function returns two values: the evaluated integral, and an estimate of the absolute error in the answer.

from scipy.integrate import quad

def integrand(X):
k = 1.0e-3
Ca0 = 1.0  # mol/L
return 1./(k*Ca0)*(1./(1-X)**2)

sol, abserr = quad(integrand, 0, 0.9)
print 't = {0} seconds ({1} hours)'.format(sol, sol/3600)
print 'Estimated absolute error = {0}'.format(abserr)

t = 9000.0 seconds (2.5 hours)
Estimated absolute error = 2.12203274482e-07


You can see the estimate error is very small compared to the solution.