Learning R Using a Chemical Reaction Engineering Book: Part 4

The links to previous parts are listed here. (Part 1, Part 2, Part 3). In this part, I tried to recreate the examples in sections A.3.1 of the computational appendix in the reaction engineering book (by Rawlings and Ekerdt).

Solving a system of ordinary differential equations

This example involves reaction (Benzene pyrolysis) in a plug flow reactor. The actual reactions happening are:

\begin{aligned}  (Rxn1) \; &2B = D + H \\  (Rxn2) \; &B + D = T+H    \end{aligned}

The rate of each reaction is given by:

r_1=k_1(c_B^2-\frac{c_Dc_H}{K_1})
r_2 = k_2(c_Bc_D-\frac{c_Tc_H}{K_2})

The feed to the reactor consists of 60kmol/hr of Benzene (B). The temperature of the reactor is T=1033K and the pressure is P=1atm.The rate constants and equilibrium constants for this example are:
k_1=7\times 10^5\; L/mol.hr,\; k_2=4\times 10^5 \; L/mol.hr,\; K_1=0.31,\; K_2=0.48


# Appendix A.3.1: Solution of Differential Equations

# Benzene pyrolysis example

# Parameters
# NBf - feed benzene rate - mol/h
# R - Universal gas constant
# T - Reactor temperature K
# P - Reactor pressure atm
# k1 - rxn1 forward rate constant L/mol.h
# k2 - rxn2 forward rate constant L/mol.h
# Keq1 - rxn1 equilibrium constant
# Keq2 - rxn2 equilibrium constant
pars=c(
NBf=60e3,
R=0.08205,
T=1033,
P=1,
k1=7e5,
k2=4e5,
Keq1=0.31,
Keq2=0.48
)

The governing equations for conversion versus volume in a plug flow reactor is based on extent of each of the reactions:
\frac{d\epsilon_1}{dV}=r_1,\; \frac{d\epsilon_2}{dV}=r_2
The initial conditions (corresponding to feed conditions N_B(0)=60kmol/h,\;N_D(0)=N_H(0)=N_T(0)=0) are that the extent of reaction is zero.
\epsilon_1(0)=0, \; \epsilon_2(0)=0

The flow rates of each component along the reactor volume can be calculated from reaction extent
N_B=N_B(0)-2\epsilon_1-\epsilon_2, \; N_D=\epsilon_1-\epsilon_2, \; N_H=\epsilon_1+\epsilon_2, \; N_T=\epsilon_2

These are setup in a function that can be passed to an ODE solver. In this case the ODE solver we use is lsode from the R package deSolve. The inputs to the function are:

  • Variable over which the integration is done (Volume in this case)
  • The state variables of the system (Extent of the two reactions)
  • Parameters that are needed for description of the system (Rate constants, Temperature, Pressure, etc.)

The output from this function is the rate of change as described by the equations previously.

# load library deSolve for solving ODEs
library(deSolve)

# function that will be passed to odesolver
# vol is the variable over which the system is integrated (equivalent of time in batch reactions)
# ext is the extent of reactions 1 and 2
# params are the parameters passed to the system
rxnrate=function(vol,ext,params) {
    with(as.list(c(ext,params)),{
      NB=NBf-2*ext1-ext2
      ND=ext1-ext2
      NH=ext1+ext2
      NT=ext2
      Q=NBf*R*T/P
      cB=NB/Q
      cD=ND/Q
      cT=NT/Q
      cH=NH/Q
      dext1=k1*(cB*cB-cD*cH/Keq1)
      dext2=k2*(cB*cD-cT*cH/Keq2)
      return(list(c(dext1=dext1,dext2=dext2)))
    })
}

Since the reaction start only after the feed enters the reactor, the extent of reaction is zero for both reactions at the beginning of the reactor (V=0L). The set of volumes where the concentration and reaction extent is computed is chosen in this case to be from 0L to 1600L at every 50L. The ODE solver lsode from deSolve package is used to solve the system of equations.


# initial extent of reaction (zero in this case for both reactions)
extinit=c(ext1=0,ext2=0)
# Volumes where the concentration is reported (in this case 0 to 1600L at every 50L)
vols=seq(0,1600,length=50)
# Solution of the set of differential equations using lsode solver in deSolve package
extout=lsode(times=vols,y=extinit,func=rxnrate,parms=pars)

extout contains the extent of reaction vs volume data. That is used to compute mole fraction and conversion at different volumes along the reactor.


# Calculation of mole fraction and conversion from extent of reaction at different volumes
extoutdf=data.frame(extout)
NBf=pars["NBf"]
extoutdf$conv=(extoutdf$ext1*2+extoutdf$ext2)/NBf
extoutdf$yB=(NBf-2*extoutdf$ext1-extoutdf$ext2)/NBf
extoutdf$yD=(extoutdf$ext1-extoutdf$ext2)/NBf
extoutdf$yT=(extoutdf$ext2)/NBf
extoutdf$yH=(extoutdf$ext1+extoutdf$ext2)/NBf

Next conversion and mole fraction is plotted as a function of reaction volume

# load library ggplot2 for plotting
library(ggplot2)
# load library reshape2 for data reshaping
library(reshape2)
# plot of conversion vs volume
ggplot(extoutdf,aes(x=time,y=conv))+geom_line()+
scale_x_continuous(breaks=seq(0,1600,by=200))+xlab('Volume (L)')+ylab('Conversion')+theme_bw(20)

# plot of mole fraction vs volume
tmp=melt(extoutdf[,c("time","yB","yD","yT","yH")],id.vars=c("time"),variable.name="moleFraction")
ggplot(tmp,aes(x=time,y=value,color=moleFraction))+geom_line()+
scale_x_continuous(breaks=seq(0,1600,by=200))+xlab('Volume (L)')+ylab('moleFraction')+theme_bw(20)

conversion

molefraction


sessionInfo()
R version 2.15.1 (2012-06-22)
Platform: i386-apple-darwin9.8.0/i386 (32-bit)

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] reshape2_1.2.1  ggplot2_0.9.2.1 deSolve_1.10-4

loaded via a namespace (and not attached):
[1] colorspace_1.1-1   dichromat_1.2-4    digest_0.5.2       grid_2.15.1
[5] gtable_0.1.1       labeling_0.1       MASS_7.3-18        memoise_0.1
[9] munsell_0.3        plyr_1.7.1         proto_0.3-9.2      RColorBrewer_1.0-5
[13] scales_0.2.2       stringr_0.6.1

About these ads
This entry was posted in R. Bookmark the permalink.

2 Responses to Learning R Using a Chemical Reaction Engineering Book: Part 4

  1. anspiess says:

    Nice example! Is it possible to do exactly the opposite, i.e. have some empirical curve data and the fit the differential equations to obtain the rate and parameter constants from the fit if they are not known?
    Cheers,
    Andrej

    • rdabbler says:

      Thanks for your comment. Yes, it is possible to fit the parameters from the data. That is in fact the next section in the computational appendix of the book and I plan to post on it in the next couple of weeks.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s