In [1]:
import scikits.bvp_solver as bvp
import numpy as np
In [2]:
T10 = 130
T2Ahx = 70
Ahx = 5
U = 1.0
In [3]:
def function(a, T):
q = (T[0] - T[1]) * U # calculate the heat transfer from stream 1 to stream 2
return np.array([-q , # evaluate dT1/dA
q/-2.0]) # evaluate dT2/dA
In [4]:
def boundary_conditions(Ta,Tb):
return (np.array([Ta[0] - T10]), #evaluate the difference between the temperature of the hot
#stream on theleft and the required boundary condition
np.array([Tb[1] - T2Ahx])) #evaluate the difference between the temperature of the cold
#stream on the right and the required boundary condition
In [5]:
problem = bvp.ProblemDefinition(num_ODE=2,
num_parameters=0,
num_left_boundary_conditions=1,
boundary_points=(0, Ahx),
function=function,
boundary_conditions=boundary_conditions)
In [6]:
solution = bvp.solve(problem,solution_guess=((T10 + T2Ahx)/2.0,
(T10 + T2Ahx)/2.0))
In [7]:
A = np.linspace(0,Ahx, 45)
T = solution(A)
In [8]:
%matplotlib inline
import matplotlib.pyplot as plt
plt.plot(A, T[0,:],'-')
plt.plot(A, T[1,:],'-')
plt.show()
In [ ]:
%matplotlib inline
import scikits.bvp_solver as bvp
import numpy as np
import matplotlib.pyplot as plt
def rhs_bvp (x, y):
return np.array([y[1], (y[2]-100)*y[0] - y[0]**3])
def bc_bvp(y1, yr):
return np.array([y1[0], y1[1]-1, yr[0]])
problem = bvp.ProblemDefinition(num_ODE=2,
num_parameters=0,
num_left_boundary_conditions=1,
boundary_points=(0, 0, -1),
function=rhs_bvp,
boundary_conditions=bc_bvp)
solution = bvp.solve(problem,solution_guess=(0., 0.), trace=2)
In [9]:
%matplotlib inline
import scikits.bvp_solver
from numpy import zeros, array, linspace
import pylab as pl
# Definition of the problem parameters
PeL = 16.0
Ci0 = 1.0
kt = 2.0
# By converting the original problem into a set of 2 firts order
# ODE's. Reaction rate is then Ri = -k*Ci0*(1-vi) the problem is defined as:
dvdz = zeros(2)
def systemode(z, v):
dvdz[0] = v[1]
dvdz[1] = PeL*(v[1]-kt*(1.0-v[0]))
return dvdz
# Defining the BC on the left (z=0) and the right (z=1) side:
BCleft = zeros(1)
BCright = zeros(1)
def boundary_conditions(va, vb):
BCleft[0] = va[1] - PeL * va[0]
BCright[0] = vb[1]
return (BCleft), (BCright)
# Problem definition: a set of 2 first order ODE with a BC on the left side.
# Because the total length was scaled, the solution is reduced within an
# interval of (0.0 , 1.0)
problem = scikits.bvp_solver.ProblemDefinition(num_ODE = 2,
num_parameters = 0,
num_left_boundary_conditions = 1,
boundary_points = (0.0, 1.0),
function = systemode,
boundary_conditions = boundary_conditions)
# By using conversion instead of concentration the solution is narrowed to
# an interval where the minimum value is 0.0 and the maximum is 1.0.
# Thus (0.0 , 1.0) is a sensible guess to initialize the problem
guess = array ([0.0 , 1.0])
# Once the problem and the BC's are being defined the solution of the problem
# is done by calling
solution = scikits.bvp_solver.solve(problem, solution_guess = guess, max_subintervals = 300)
# Getting the solution from z =0 to z = 1.0
z = linspace (0.0, 1.0, 100, endpoint = True )
v = solution(z)
# Comparison of the results
print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
print 'The analytical solution of the problem at the outlet is: '
print 'Final Conversion= 0.83605373297390484'
print 'Final Concentration = 0.16394626702609519'
print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
print 'The numerical solution of the problem at the outlet is: '
print 'Final Conversion= ', v[0,99]
print 'Final Concentration = ', 1.0- v[0,99]
print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
# Plotting of the solution along z
pl.figure()
pl.subplot(1,2,1)
pl.plot(z, v[1,:], 'b-', z, v[0,:], 'g-')
pl.legend(('dv/dz[0]', 'dv/dz[1]'))
pl.xlabel('Normalized Reactor Length')
pl.ylabel('dv/dz')
pl.title('System of ODEs for the Axial Dispersion Model')
pl.grid('on')
pl.subplot(1,2,2)
pl.plot(z, v[0,:], 'g-')
pl.xlabel('Normalized Reactor Length')
pl.ylabel('Conversion of Reactant')
pl.title('Conversion produced by the Axial Dispersion Model')
pl.grid('on')
pl.show()