Skip to content

Commit

Permalink
version 1.0
Browse files Browse the repository at this point in the history
  • Loading branch information
polgia0 authored and gaborcsardi committed Dec 23, 2014
0 parents commit eeee759
Show file tree
Hide file tree
Showing 57 changed files with 3,137 additions and 0 deletions.
17 changes: 17 additions & 0 deletions DESCRIPTION
@@ -0,0 +1,17 @@
Package: HomoPolymer
Type: Package
Title: Theoretical Model to Simulate Radical Polymerization
Version: 1.0
Date: 2014-12-23
Author: Gianmarco Polotti
Maintainer: Gianmarco Polotti <gianmarco.polotti@gmail.com>
Description: A theoretical model to simulate radical polymerization. Material, energy and population balances are integrated for batch, semi batch and continuous process in a ideally mixed reactor. Limitations: single monomer (i.e.homo polymer), one phase (organic, aqueous). Datasets with chemical-physical data for the most common monomers is included. Some background in polymer science is suggested for its use. Graphical interface for a quick and friendly use is available.
Depends: R (>= 3.0.0), RGtk2, MenuCollection
Imports: deSolve
LazyData: true
LazyLoad: true
License: GPL-2
Packaged: 2015-01-11 19:21:42 UTC; PolGia0
NeedsCompilation: no
Repository: CRAN
Date/Publication: 2015-01-12 10:27:23
56 changes: 56 additions & 0 deletions MD5
@@ -0,0 +1,56 @@
e6d285fe519b02823fd5c5f318e269a7 *DESCRIPTION
92ebad46ec23b191b25255c90e3d3144 *NAMESPACE
c1f9a4f720a0d90331ce550151f49c63 *R/PlotBasic.R
b8fb117ecce67b5eed9018ea748b38b4 *R/PlotChoice.r
88b2ce113a8c17b34f0dcfc72ec4a17b *R/StopConditions.R
01363e3fa086be4df2411f7c96258ad1 *R/func_homo.R
3849e33a02ee68d8a96f0ff4435961f8 *R/gel.R
62e125e291fde88ce56e522afae96648 *R/inpbox_computational_options.R
2570e2ccaf9fef2d1eb6eb4db02f2841 *R/inpbox_flow_options.R
0c2c727f147fbc8c084460d13b71c2c3 *R/inpbox_reaction_formulation.R
1b459145f1f5e4c864ec2efd414266ac *R/integral.R
a8d8256ad9de7c26539dfada4538eabf *R/polymer_toolbar.r
606f5028a713d31fc5e46e88dab96432 *R/prepare.R
043fce87272549b9e83e81db1585cab6 *build/vignette.rds
326fb3abe08f21af0afd36533c56ac9b *data/DBk.rda
d3150d004dfb47e2198e6d2b69df69c5 *data/DBpp.rda
b71ab0cbae065ccab0ec0eeba0cd22b9 *inst/doc/Data_and_model_equations.Rnw
47e1aae01299a26524f5d476e672f261 *inst/doc/Data_and_model_equations.pdf
f71ac6bf27f738ff993984af5eecd58b *inst/doc/User_manual.Rnw
a5ec67fdd3045cbf3efea8e1420c1721 *inst/doc/User_manual.pdf
614e2627ba3a278927b45343e71da2fb *man/DBk.Rd
ad3de7a4a3ade104ae0036671bf463ed *man/DBpp.Rd
6b48bcb01e7bdebdc2bbe74f5a4c249a *man/PlotBasic.Rd
efb4e3fb852c51a64c6679c6a778265a *man/PlotChoice.Rd
1915c6b7578b61fdb4cb52a43dc0969f *man/StopConditions.Rd
9cde43c3138650a60bdb0b892379de36 *man/figures/flow.png
9cdaf63788ce5a277a9a512ed4e8fe98 *man/figures/process.png
4d76ab531454b085f5ba9d42d986b7cb *man/figures/recipe.png
7100825215483e5547a75540356c2833 *man/func_homo.Rd
35ae90114220d6b31b4e227e04afce04 *man/gel.Rd
490f1d5d73f4139b76603e0544b48cd2 *man/homopolymer-package.Rd
77158e0ac8eed3e248877d844fafd6c8 *man/inpbox_computational_options.Rd
b0d43d879e58c8cdaab21fb346f9e551 *man/inpbox_flow_options.Rd
993783d7781d5c514260822dba20dd45 *man/inpbox_reaction_formulation.Rd
31bd837c1a971b912a24425b1858a425 *man/integral.Rd
46709e6c1dc20d937bbf96f40f66f95e *man/polymer_toolbar.Rd
a29ceb63bf5bea4772e867c3ad38ce31 *man/prepare.Rd
b71ab0cbae065ccab0ec0eeba0cd22b9 *vignettes/Data_and_model_equations.Rnw
f71ac6bf27f738ff993984af5eecd58b *vignettes/User_manual.Rnw
f8e59ceb3b83d7b74902b5745de2ce4d *vignettes/images/AA.png
cf15f53682d14fffbc7deaf4a00c3281 *vignettes/images/AA1.png
9889a948ffdb26858a7435f96ae28ed5 *vignettes/images/AA1r.png
f8ea1046aecc0ca0ce6d9ccb86c1643d *vignettes/images/AAr.png
0f85949842e5f74757d07a2e08b7d95b *vignettes/images/EA.png
f992842ff4cf9aa967d3b6472ca81a73 *vignettes/images/EAr.png
54d16c5550db5b11a6df5926141d7529 *vignettes/images/MMA.png
38b15b2646597f9c1b706997bd1ec66f *vignettes/images/MMAr.png
1cf386627bf04e0ee6813a5f135fa71c *vignettes/images/fig1.png
3755c4313d5e0d97fada832a947a339b *vignettes/images/fig2.png
78626c61a7cc324fa07f6bce947bd062 *vignettes/images/fig3.png
2a14cca51b38e365422c7b6ccd5df0e5 *vignettes/images/fig4.png
2e872b9f91bdc6077bc2d256b3daf09a *vignettes/images/fig5.png
84bb73b12161f799cd561ebe7ee6dc22 *vignettes/images/fig6.png
9905ffa694b18d7620503674fffffc46 *vignettes/images/fig7.png
6938d59146f1dadd4a3197df87d32915 *vignettes/images/schema.png
85db6ff7fe5e58d5b14a8427443e5583 *vignettes/images/toolbar.png
16 changes: 16 additions & 0 deletions NAMESPACE
@@ -0,0 +1,16 @@
import(RGtk2)
import(MenuCollection)
import(deSolve)
## Exported functions:
export(StopConditions)
export(prepare)
export(polymer_toolbar)
export(PlotChoice)
export(PlotBasic)
export(integral)
export(gel)
export(func_homo)
export(inpbox_reaction_formulation)
export(inpbox_flow_options)
export(inpbox_computational_options)

13 changes: 13 additions & 0 deletions R/PlotBasic.R
@@ -0,0 +1,13 @@
PlotBasic <-function(out,pars){
dev.new()
op<-par(mfrow=c(2,2),pty = "s",mar=c(4,3,1,1),cex.axis=0.8,tck=-0.01,cex.lab=0.8,font=2,mgp=c(2,1,0))
plot(out$time,out$X*100,type='l',xlab='Time(min)',ylab='X(%)',col='red');grid()
plot(out$time,out$M*pars['MWM']/1000*out$Vl,type='l',xlab='Time(min)',ylab='M(Kg)',col='red');grid()
plot(out$time,out$T-273.16,type='l',xlab='Time(min)',ylab='Temperature(C)',col='red');grid()
lgMn<-log10(out$Mnm)
lgMw<-log10(out$Mwm)
plot(out$time,lgMn,type='l',xlab='Time(min)',ylab='log Mn(red), log Mw(blue) in (g/mol)',col='red',
ylim=c(min(lgMn[-1],lgMw[-1]),max(lgMn[-1],lgMw[-1])));grid()
lines(out$time,lgMw,col='blue')
par(op)
}
40 changes: 40 additions & 0 deletions R/PlotChoice.r
@@ -0,0 +1,40 @@
PlotChoice<-function(out,otab){
plotgraph<-function(out,vx,vy,vz=NULL){
par(mar=c(5,4,4,5)+.1)
xlim<-c(min(out[,vx]),max(out[,vx]))
ylim<-c(min(out[,vy]),max(out[,vy]))
plot(out[,vx],out[,vy[1]],xlim=xlim,ylim=ylim,type='n',
xlab=otab$text[ix],ylab=otab$text[iy],col.lab='red')
grid()
for(i in 1:length(vy))points(out[,vx],out[,vy[i]],col='red')
for(i in 1:length(vy))lines(out[,vx],out[,vy[i]],col='red')
if(!is.null(vz)){
par(new=TRUE)
zlim<-c(min(out[,vz]),max(out[,vz]))
plot(out[,vx],out[,vz[1]],type='n',xaxt="n",yaxt="n",xlab=otab$text[ix],
ylab=otab$text[iy],,col.lab='red',xlim=xlim,ylim=zlim)
for(i in 1:length(vz))points(out[,vx],out[,vz[i]],col='blue')
for(i in 1:length(vz))lines(out[,vx],out[,vz[i]],col='blue')
axis(4)
zlab<-otab$text[iz]
mtext(zlab,side=4,line=3,col='blue')
}
}
var.list<-otab$text
ans<-inpboxck(c('Variable on x-axis','Secondary y-axis'),var.list,c('-1','FALSE'))
ix<-as.numeric(ans[[1]])
vx<-which(names(out)==otab$symbol[ix])
vz<-ans[[2]]
ans<-inpboxlist(var.list,'Variables on y main axis')
iy<-as.numeric(ans[[1]])
vy<-which(names(out)==otab$symbol[iy])
if(as.logical(vz)){
ans<-inpboxlist(var.list,'Variables on y secondary axis')
iz<-as.numeric(ans[[1]])
vz<-which(names(out)==otab$symbol[iz])
plotgraph(out,vx,vy,vz)
}else{
plotgraph(out,vx,vy)
}
}

11 changes: 11 additions & 0 deletions R/StopConditions.R
@@ -0,0 +1,11 @@
StopConditions<-function(t,y,pars,...){
yroot<-rep(1,6)
yroot[1]<- y['M']
yroot[2]<-y['X']-pars['Xlim']
yroot[3]<- y['I']
if(pars['pH0']!=0){
yroot[4]<-y['HA']-1e-8
# yroot[5]<-y['M']-y['HA']
}
return(yroot)
}
220 changes: 220 additions & 0 deletions R/func_homo.R
@@ -0,0 +1,220 @@
func_homo <-function(t,y,pars,Kpf,Kdf,Kfz,Kfm,Kfs,Kfp,Kps,Kpss,KfCTA,KfCCTA,Ktf,Ktd,Ktc,roM,roS,roP){
with(as.list(c(y,pars)),{
#browser()
Kd1<-function(T){return(1.178855e16*exp(-28925.227/1.987/T))}
Kd2<-function(T){return(1.3e15*exp(-27756.9/1.987/T))}
Kc<-function(X){return(3*X^2-3.17*X+2.9)}
FWin<-FinM+FinI+FinZ+FinS+FinCTA+FinCCTA+FinNa
Fin<-FinM+FinI+FinZ #mol/min
#composition of in flow
if(FinS!=0)Fin<-Fin+FinS
if(FinCTA!=0)Fin<-Fin+FinCTA
if(FinCCTA!=0)Fin<-Fin+FinCCTA
if(FinNa!=0)Fin<-Fin+FinNa
uM<-0
uI<-0
uZ<-0
uS<-0
uCTA<-0
uCCTA<-0
uNa<-0
Cpin<-0
roin<-0
if(Fin!=0){
uM<-FinM/Fin #molar fraction
uI<-FinI/Fin #molar fraction
uZ<-FinZ/Fin #molar fraction
if(FinS!=0)uS<-FinS/Fin #molar fraction
if(FinCTA!=0)uCTA<-FinCTA/Fin #molar fraction
if(FinCCTA!=0)uCCTA<-FinCCTA/Fin #molar fraction
if(FinNa!=0)uNa<-FinNa/Fin #molar fraction
Cpin<-FinM/FWin*MWM/1000*CpM+FinS/FWin*MWS/1000*CpS #cal/mol/k
roin<-uM*roM(Tin)*1000/MWM #mol/L
if(uS!=0)roin<-roin+uS*roS(Tin)*1000/MWS #mol/L
}
Hin<-H0+Cpin*(Tin-T0) #cal/mol
#composition of out flow
Mtm<-M+I+S+Z+P+Na+CTA+CCTA #mol/L
zM<-M/Mtm #molar fraction
zI<-I/Mtm #molar fraction
zS<-S/Mtm #molar fraction
zZ<-Z/Mtm #molar fraction
zP<-P/Mtm #molar fraction
zNa<-Na/Mtm #molar fraction
zCTA<-CTA/Mtm #molar fraction
zCCTA<-CCTA/Mtm #molar fraction
MWm<-(zM+zP)*MWM+zI*MWI+zS*MWS+zZ*MWZ+zNa*MWNa+zCTA*MWCTA+zCCTA*MWCCTA #g/mol
Fout<-FWout*1000/MWm #mol/min
#gel effect parameter determination
Cgel<-gel(y,pars,roM,roS,roP)
Ktb<-Ktf(T)
Ktt<-Ktb
Ktr<-0
Kts<-0
if(DCT)Ktt<-Ktb*Cgel$Ft
if(ST)Kts<-Ktb*Cgel$Fs
if(RT){
Ktrmin<-Kpf(T)*M*Cgel$Ktrmin
Ktrmax<-Kpf(T)*M*Cgel$Ktrmax
Ktr<-Ktrmin+(Ktrmax-Ktrmin)*X
}
if((Ktr!=0)&(Kts!=0))Ktb<-1/(1/Ktt+1/Kts+1/Ktr)
if((Ktr==0)&(Kts!=0))Ktb<-1/(1/Ktt+1/Kts)
if((Ktr!=0)&(Kts==0))Ktb<-1/(1/Ktt+1/Ktr)
if((Ktr==0)&(Kts==0))Ktb<-Ktt
Ktd<-Ktd(T)/Ktf(T)*Ktb
Ktc<-Ktc(T)/Ktf(T)*Ktb
Kp<-Kpf(T)
if(DCP)Kp<-Kp*Cgel$Fp
Kd<-Kdf(T)
if(DCI)Kd<-Kd*Cgel$Fd
Kd1<-Kd1(T)
if(DCI)Kd1<-Kd1*Cgel$Fd
Kd2<-Kd2(T)
if(DCI)Kd2<-Kd2*Cgel$Fd
Cp<-zM*CpM*MWM/1000+zS*MWS/1000*CpS+zP*CpP*MWM/1000 #cal/mol/k
roout<-zM*roM(T)*1000/MWM+zP*roP(T)*1000/MWM #mol/L
if(MWS!=0)roout<-roout+zS*roS(T)*1000/MWS #mol/L
#(A3.14)
if(pH0!=0){
# case of water soluble monomers
AM<-M-HA
if(DPH){
pH<--log10(HA/AM)+pH0+log10(HA0/(M0-HA0))
if(pH>14)pH<-14
}else{
pH<-pH0
}
#(A3.8)
RI<-2*Kd*I+2*Kd1*HA+2*Kd2*AM #mol/L/min
#(A1.25) Radical Balance under SSA
R<-0.5*(sqrt((Kfz(T)*Z/Ktb)^2+4*RI/Ktb)-Kfz(T)*Z/Ktb) #mol/L
#(12.4)
alfa2<-1/(1+10^(KdisP-pH))
alfa3<-1/(1+10^(KdisP-pH))
HP<-P/(1+alfa2)
PM<-alfa2*P/(1+alfa2)
HR<-R/(1+alfa3)
RM<-alfa3*R/(1+alfa3)
#(A3.22)
B<-Kc(X)*(HP-Na-P)-1
C<-Kc(X)*Na*(P-HP)
Pc<-(-B-sqrt(B^2-4*Kc(X)*C))/(2*Kc(X))
#(A3.23)
Rc<-Kc(X)*(Na-Pc)*RM
#(A3.16)
Kpp<-Kp*(HA/M+Kratio*HR*AM/(R*M)+Kratio*Rc*AM/(R*M))
}else{
# case of oil soluble monomers
AM<-0
HA<-M
pH<-0
KdisM<-0
KdisP<-0
RI<-2*Kd*I #mol/L/min
R<-0.5*(sqrt((Kfz(T)*Z/Ktb)^2+4*RI/Ktb)-Kfz(T)*Z/Ktb) #mol/L
alfa2<-1
alfa3<-1
HP<-P
PM<-0
HR<-R
RM<-0
Pc<-0
Rc<-0
Kpp<-Kp
}
#(A1.29) Enthalpy
H<-H0+Cp*(T-T0) #cal/mol
#(A1.32)
Tau<-Ktd*(Kpp*M*R)/(Kpp*M)^2+Kfm(T)/Kpp+Kfs(T)*S/(Kpp*M)+KfCTA(T)*CTA/(Kpp*M)+KfCCTA(T)*CCTA/(Kpp*M)+Kfz(T)*Z/(Kpp*M)
#(A1.33)
Beta<-Ktc*(Kpp*M*R)/(Kpp*M)^2
#(A1.42)
Gam<-1+RI/(Kpp*M*R)+Kfm(T)/Kpp+Kfs(T)*S/(Kpp*M)+KfCTA(T)*CTA/(Kpp*M)+KfCCTA(T)*CCTA/(Kpp*M)+Kfz(T)*Z/(Kpp*M)
#(A1.43)
Lam<-Ktb*R/(Kpp*M)+Kfm(T)/Kpp+Kfs(T)*S/(Kpp*M)+KfCTA(T)*CTA/(Kpp*M)+KfCCTA(T)*CCTA/(Kpp*M)+Kfz(T)*Z/(Kpp*M)+Kfp(T)*Mu1/(Kpp*M)
#(A1.44)
Eta<-1+Kfm(T)/Kpp+Kfs(T)*S/(Kpp*M)+KfCTA(T)*CTA/(Kpp*M)+KfCCTA(T)*CCTA/(Kpp*M)+Kfz(T)*Z/(Kpp*M)+Kps(T)*Mu1/(Kpp*M)+RI/(Kpp*M*R)
+(Kfp(T)/Kpp+Kpss(T)/Kpp)*Mu2/M
if(LIN){
#(A1.30)
Mn<-MWM/(Tau+Beta/2)
#(A1.31)
Mw<-MWM*(2*Tau+3*Beta)/(Tau+Beta)^2
}else{
#(A1.45)
Mn<-MWM*Mu1/Mu0
if(is.na(Mn))Mn<-MWM
#(A1.46)
Mw<-MWM*Mu2/Mu1
if(is.na(Mw))Mw<-MWM
}
if(Mn<MWM)Mn<-MWM
if(Mw<MWM)Mw<-MWM
#Derivative output vector
yp<-rep(0,19)
#(A1.18) Monomer
dM<-Kpp*M*R #mol/L/min
yp[1]<-Fin*uM-dM*Vl-Fout*zM #mol/min
#molar conversion
yp[2]<--yp[1]/M0 #1/min
#(A3.12) undissociated monomer
if(pH0!=0){
yp[4]<--Kp*R*HA*Vl #mol/min
}else{
yp[4]<-yp[1]
}
#(A1.19) Initiator
yp[5]<-Fin*uI-2*Kd*I*Vl-Fout*zI #mol/min
#(A1.20) Solvent
if(S0!=0)yp[6]<-Fin*uS-Kfs(T)*S*R*Vl-Fout*zS #mol/min
#(A1.21) Polymer
yp[7]<-dM*Vl-Fout*zP #mol/min
#(A1.22) Inhibitor
if(Z0!=0)yp[8]<-Fin*uZ-Kfz(T)*Z*R*Vl-Fout*zZ #mol/min
#(A1.23) CTA
if(CTA0!=0)yp[9]<-Fin*uCTA-KfCTA(T)*CTA*R*Vl-Fout*zCTA #mol/min
#(A1.24) CCTA
if(CCTA0!=0)yp[10]<-Fin*uCCTA-KfCCTA(T)*CCTA*R*Vl-Fout*zCCTA #mol/min
# Counter ions balance
if(Na0!=0)y[11]<-Fin*uNa-Fout*zNa #mol/min
#(A1.27) Energy balance
if(!ISOT){
yp[12]<-Fin*Hin+dM*Vl*DHp-UA*(T-Tj)-Fout*H
yp[12]<-yp[12]/(Vl*roout*Cp) #k/min
}else{
yp[12]<-0
}
#total volume
if(Fin!=0)yp[3]<-yp[3]+Fin/roin
if(Fout!=0)yp[3]<-yp[3]-Fout/roout
yp[3]<-yp[3]+MWM/1000*(1/roM(T)-1/roP(T))*yp[1]
yp[3]<-yp[3]+MWM/1000*M*(roP2/roP(T)^2-roM2/roM(T)^2)*yp[12] #L/min
#(A1.45)
yp[13]<-dM*Vl/Mn #mol^2/min/g
#(A1.46)
yp[14]<-dM*Vl*Mw #g/min
if(LIN){
yp[15]<-0
yp[16]<-0
yp[17]<-0
yp[18]<-0
yp[19]<-0
}else{
#(A1.39) 0th moment
yp[15]<-(Tau+Beta/2-Kpss(T)*Mu1/(Kpp*M)-Kps(T)*Mu0/(Kpp*M))*Kpp*M*R*Vl-Mu0*Fout*Vl
#(A1.40) 1st moment
yp[16]<-(Tau-Ktd*R/(Kpp*M))*Kpp*M*R*Vl-Mu1*Fout*Vl
#(A1.41) 2nd moment
yp[17]<-(Gam+2*(1+(Kps(T)*Mu1+Kpss(T)*Mu2)/(Kpp*M))*Eta/Gam+Ktc*R/(Kpp*M)*(Eta/Lam)^2)*Kpp*M*R*Vl-Mu2*Fout*Vl
#(A1.47)
yp[18]<-(Kfp(T)*Mu1*R+Kps(T)*R*Mu0-Mu0*BN3*Fout)*Vl
#(A1.48)
yp[19]<-(Kpss(T)*Mu2*R-Mu0*BN4*Fout)*Vl
}
der<-c(AM,R,Kp,Ktb,Kd,Mn,Mw,Cp,roin,roout,pH,Kpp)
attributes(der)<-NULL
names(der)<-c("AM","R","Kp","Ktb","Kd","Mn","Mw","Cp","roin","roout","pH","Kpp")
return(list(yp,der))
})
}

0 comments on commit eeee759

Please sign in to comment.