Free Radical Polymerization for Copolymers AA and Styrene using pseudokinetic rate constants

9 次查看（过去 30 天）
Oreoluwa Ogunnaike
评论： Davide Masiello 2022-5-5
Hello, I have a MATLAB project for my polymer engineering class where I am supposed to use pseudokinetics to model the MW distribution and polymer concentration, etc of a styrene and acrylic acid copolymer. The first part of the project involved only acrylic acid and I will include the code below. Since this new case involves two polymers, I'm having a hard time and many sleepless nights making it into just one ode. I don't have much of a background in MATLAB (i have rarely used it prior to this semester, so ode45 is the only function I know for solving ODEs. If there's a better one I could use, I'd love to know, and I'd love any help with this project in general.
Part 1 code:
global kp kt kd f MWo
tf=15000;
tvect=[0 tf];
Mo=3.2; %mol/L changing monomer concentration +/- 25% changes polymer MW
Io=0.01; %mol/L
Mu=0;
MWo=72; %MWn=DPn*MWo, MWw=DPw*MWo
kp=950; %L/mol.s
Ep=6300; %cal/mol
T=348;%K
r=1.9872; %cal/K.mol
kt=110000000*exp(-2800/(r*T)); %L/mol.s
kd=140000000000000*exp(-30600/(r*T)); %s^-1
f=0.7;
%Ro=sqrt((2*f*kd*Io)/kt); %initial condition=Y1=Y2 ICs
%Y1= (kp*Mo/kt)+Yo; %initial condition
%Y2= ((2*kp*Mo*Y1)+(kp*Mo*Yo)+(kt*Yo^2))/(kt*Yo);
Ci0=[Mo Io Mu Mu Mu];
[t,Ci]=ode45(@monomer,tvect,Ci0);
M=Ci(:,1);I=Ci(:,2);muo=Ci(:,3);mu1=Ci(:,4);mu2=Ci(:,5);
%dmu1./dmu0=Y1./Yo;
%DPw=Y2./Y1;
%tau=kt.*Yo./(kp.*M);
Yo=sqrt((2*f*kd.*I)/kt);
tau=kt.*Yo./(kp.*M);
Y1= (kp.*M/kt)+Yo;
Y2=((2*kp.*M.*Y1)+(kp.*M.*Yo)+(kt.*Yo.^2))./(kt.*Yo);
MWn=(Y1./Yo)*MWo;
MWw=(Y2./Y1)*MWo;
DPnd=mu1./muo;
DPwd=mu2./mu1;
DPw=Y2./Y1;
DPn=Y1./Yo;
MWdn=DPnd*MWo;
MWdw=DPwd*MWo;
figure(1)
plot(t,Yo)
title('Yo vs t')
figure(6)
plot(t,Y1)
title('Y1 vs t')
figure(7)
plot(t,Y2)
title('Y2 vs t')
figure (10)
plot(t,muo)
title('Muo vs t')
figure (11)
plot(t,mu1)
title('Mu1 vs t')
figure (12)
plot(t,mu2)
title('Mu2 vs t')
figure(2)
plot(t,tau)
title('tau vs t')
figure(3)
plot(t,M)
title('[M] vs t')
figure(20)
plot(t,I)
title('[I] vs t')
figure(4)
plot(t,MWn,t,MWw)
legend('live MWn','live MWw')
title('live molecular weight distributions vs t')
figure(9)
plot(t,MWdn,t,MWdw)
legend('dead MWn','dead MWw')
title('dead molecular weight distributions vs time')
X=((Mo-M)./Mo)*100;
figure(15)
plot(t,X)
title('Monomer conversion vs time')
function dCdT=monomer(t,vars)
global kp kt kd f
M1=vars(1); I=vars(2);muo=vars(3); mu1=vars(4); mu2=vars(5);
dM=-kp*M*sqrt((2*f*kd.*I)/kt);
dI=-kd*I;
dmu0=2*f*kd*I;
dmu1=((kp*M/kt)+sqrt(2*f*kd.*I/kt))*kt.*sqrt((2*f*kd.*I)/kt);
dmu2=((2*kp.*M.*((kp.*M/kt)+sqrt(2*f*kd.*I/kt))+(kp.*M.*sqrt(2*f*kd.*I)/kt))+(kt.*(2*f*kd.*I/kt)))./(kt.*sqrt(2*f*kd.*I/kt))*kt.*sqrt(2*f*kd.*I/kt);
%muo=sqrt((2*f*kd.*I)/kt);
%mu1= (kp.*M/kt)+muo;
%mu2=((2*kp.*M.*mu1)+(kp.*M.*muo)+(kt.*muo.^2))./(kt.*muo);
%tau;
dCdT=[dM;dI;dmu0;dmu1;dmu2];
end
Part 2 (idk what I'm doing) code:
global kp11 kp12 kp21 kp22 ktc22 ktd11 ktd12 ktd22 f kd
tf=15000;
tvect=[0 tf];
M1o=1.6; %mol/L changing monomer concentration +/- 25% changes polymer MW
M2o=1.6; %styrene
Mo=M1o+M2o;
Io=0.01; %mol/L
Mu=0;
MW1=72; %MWn=DPn*MWo, MWw=DPw*MWo
MW2=104;
T = 348;
R = 1.987;
r1 = 0.27;
r2 = 0.72;
kp11 = 950;
kp12 = kp11/r1;
kp22 = (4.3e7)*exp(-7769/(R*T));
kp21 = kp22/r2;
ktc22 = (1.06e9)*exp(-1496/(R*T));
ktd11 = (1.1e8)*exp(-2800/(R*T));
ktd22 = (1.52e8)*exp(-1496/(R*T));
ktd12 = sqrt(ktd11*ktd22);
f = 0.7;
kd = (1.4e14)*exp(-30600/(R*T));
%Ro=sqrt((2*f*kd*Io)/kt); %initial condition=Y1=Y2 ICs
%Y1= (kp*Mo/kt)+Yo; %initial condition
%Y2= ((2*kp*Mo*Y1)+(kp*Mo*Yo)+(kt*Yo^2))/(kt*Yo);
Ci0=[Mo Io Mu Mu Mu];
[t,Ci]=ode45(@monomer,tvect,Ci0);
M=Ci(:,1);I=Ci(:,2);muo=Ci(:,3);mu1=Ci(:,4);mu2=Ci(:,5);
%dmu1./dmu0=Y1./Yo;
%DPw=Y2./Y1;
%tau=kt.*Yo./(kp.*M);
phi1 = (kp21.*M1)./((kp21.*M1)+(kp12.*f2.*M));
phi2 = 1-phi1;
f1 = M1./(M);
f2 = 1-f1;
M2=M-M1;
ktc = (ktc22.*phi2.^2);
ktd = (ktd11*phi1.^2)+(2*ktd12*phi1*phi2)+(ktd22*phi2.^2);
Yo = sqrt((2.*f.*kd.*I)./(ktc+ktd));
kp = (kp11*phi1*f1)+(kp12*phi1*f2)+(kp21*phi2*f1)+(kp22*phi2*f2);
Y1 = (kp.*M./(ktc+ktd))+Yo;
Y2 = (2.*kp.*M.*Y1./((ktc+ktd).*Yo))+(kp.*M/(ktc+ktd))+Yo;
tau=ktd.*Yo./(kp.*M);
beta=ktc*Yo./(kp.*M);
%M=M1+M2;
%MWn=(Y1./Yo)*MWo;
%MWw=(Y2./Y1)*MWo;
%DPnd=mu1./muo;
%DPwd=mu2./mu1;
%DPw=Y2./Y1;
%DPn=Y1./Yo;
%MWdn=DPnd*MWo;
%MWdw=DPwd*MWo;
%figure(1)
%plot(t,Yo)
%title('Yo vs t')
%figure(6)
%plot(t,Y1)
%title('Y1 vs t')
%figure(7)
%plot(t,Y2)
%title('Y2 vs t')
%figure (10)
%plot(t,muo)
%title('Muo vs t')
%figure (11)
%plot(t,mu1)
%title('Mu1 vs t')
%figure (12)
%plot(t,mu2)
%title('Mu2 vs t')
%figure(2)
%plot(t,tau)
%title('tau vs t')
%figure(3)
%plot(t,M)
%title('[M] vs t')
%figure(20)
%plot(t,I)
%title('[I] vs t')
%figure(4)
%plot(t,MWn,t,MWw)
%legend('live MWn','live MWw')
%title('live molecular weight distributions vs t')
%figure(9)
%plot(t,MWdn,t,MWdw)
%legend('dead MWn','dead MWw')
%title('dead molecular weight distributions vs time')
%X=((Mo-M)./Mo)*100;
%figure(15)
%plot(t,X)
%title('Monomer conversion vs time')
function dCdT=monomer(t,vars,kp11, kp12, kp21, kp22, ktc22, ktd11, ktd12, ktd22, f, kd)
M=vars(1);I=vars(2);%Y1=vars(8);Y2=vars(9);
dM=-((kp11*phi1*Yo*f1*M)+(kp12*phi1*Yo*f2*M)+(kp21*phi2*Yo*f1*M)+(kp22*phi2*Yo*f2*M));
dI=-kd*I;
dmu0=2*f*kd*I;
dmu1= kp*(M1+M2)*Y1*(tau+beta);
dmu2= kp*(M1+M2)*((tau*Y2)+(beta*(Y2+((Y1^2)/Yo))));
%muo=sqrt((2*f*kd.*I)/kt);
%mu1= (kp.*M/kt)+muo;
%mu2=((2*kp.*M.*mu1)+(kp.*M.*muo)+(kt.*muo.^2))./(kt.*muo);
%tau;
dCdT=[dM;dI;dmu0;dmu1;dmu2];
end
  7 个评论
显示 5更早的评论
Oreoluwa Ogunnaike
Oreoluwa Ogunnaike 2022-5-5
Hi Davide, I cannot find any other possible constraint, but it made me wonder why i couldnt have two separate ODEs for M1 and M2. I have adjusted the code for that. I still have the issue of figuring out how to use the equations in the ODE.
clear clc
%Constants
M1o=1.6; %mol/L changing monomer concentration +/- 25% changes polymer MW
M2o=1.6; %styrene
Io=0.01; %mol/L
Mu=0;
MW1=72; %MWn=DPn*MWo, MWw=DPw*MWo
MW2=104;
T = 348;
R = 1.987;
r1 = 0.27;
r2 = 0.72;
kp11 = 950;
kp12 = kp11/r1;
kp22 = (4.3e7)*exp(-7769/(R*T));
kp21 = kp22/r2;
ktc22 = (1.06e9)*exp(-1496/(R*T));
ktd11 = (1.1e8)*exp(-2800/(R*T));
ktd22 = (1.52e8)*exp(-1496/(R*T));
ktd12 = sqrt(ktd11*ktd22);
f = 0.7;
kd = (1.4e14)*exp(-30600/(R*T));
%Numerical solution
Ci0=[M1o M2o Io Mu Mu Mu];
tvect = [0 15000];
[t,Ci]=ode45(@(t,vars)monomer(t,vars,kd,kp11,kp12,kp21,kp22,ktc22,ktd11,ktd12,ktd22,f),tvect,Ci0);
% Post-processing
M1=Ci(:,1);
M2=Ci(:,2);
I=Ci(:,3);
muo=Ci(:,4);
mu1=Ci(:,5);
mu2=Ci(:,6);
%dmu1./dmu0=Y1./Yo;
%DPw=Y2./Y1;
%tau=kt.*Yo./(kp.*M);
phi1 = (kp21.*M1)./((kp21.*M1)+(kp12.*M2));
phi2 = 1-phi1;
f1 = M1./(M1+M2);
f2 = 1-f1;
ktc = (ktc22.*phi2.^2);
ktd = (ktd11*phi1.^2)+(2*ktd12*phi1*phi2)+(ktd22*phi2.^2);
Yo = sqrt((2.*f.*kd.*I)./(ktc+ktd));
kp = (kp11*phi1*f1)+(kp12*phi1*f2)+(kp21*phi2*f1)+(kp22*phi2*f2);
Y1 = (kp.*(M1+M2)./(ktc+ktd))+Yo;
Y2 = (2.*kp.*(M1+M2).*Y1./((ktc+ktd).*Yo))+(kp.*(M1+M2)/(ktc+ktd))+Yo;
M=M1+M2;
tau=ktd.*Yo./(kp.*M);
beta=ktc*Yo./(kp.*M);
%MWn=(Y1./Yo)*MWo;
%MWw=(Y2./Y1)*MWo;
%DPnd=mu1./muo;
%DPwd=mu2./mu1;
%DPw=Y2./Y1;
%DPn=Y1./Yo;
%MWdn=DPnd*MWo;
%MWdw=DPwd*MWo;
%figure(1)
%plot(t,Yo)
%title('Yo vs t')
%figure(6)
%plot(t,Y1)
%title('Y1 vs t')
%figure(7)
%plot(t,Y2)
%title('Y2 vs t')
%figure (10)
%plot(t,muo)
%title('Muo vs t')
%figure (11)
%plot(t,mu1)
%title('Mu1 vs t')
%figure (12)
%plot(t,mu2)
%title('Mu2 vs t')
%figure(2)
%plot(t,tau)
%title('tau vs t')
%figure(3)
%plot(t,M)
%title('[M] vs t')
%figure(20)
%plot(t,I)
%title('[I] vs t')
%figure(4)
%plot(t,MWn,t,MWw)
%legend('live MWn','live MWw')
%title('live molecular weight distributions vs t')
%figure(9)
%plot(t,MWdn,t,MWdw)
%legend('dead MWn','dead MWw')
%title('dead molecular weight distributions vs time')
%X=((Mo-M)./Mo)*100;
%figure(15)
%plot(t,X)
%title('Monomer conversion vs time')
%ODE system
function dCdT=monomer(t,vars,kd,kp11,kp12,kp21,kp22,ktc22,ktd11,ktd12,ktd22,f)
M1=vars(1);M2=vars(2);I=vars(3);mu0=vars(4);mu1=vars(5);mu2=vars(6);
dM1=-kp*M1*Yo;
dM2=-kp*M2*Yo;
dI=-kd*I;
dmu0=2*f*kd*I;
dmu1= kp*(M1+M2)*Y1*(tau+beta);
dmu2= kp*(M1+M2)*((tau*Y2)+(beta*(Y2+((Y1^2)/Yo))));
dCdT=[dM1;dM2;dI;dmu0;dmu1;dmu2];
end
Davide Masiello
Davide Masiello 2022-5-5
Well, then it's easy enough.
I have attached an example in my answer below, so that you can click the accept button if you think it solves your problem.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Davide Masiello
Davide Masiello 2022-5-5
You just need to calculate all the quantitites inside the ODE function as well.
If you need those quantities to be plotted, then you may recompute them after the integration as been performed (see the post-processing section in the example below).
clear, clc
% Constants
M1o = 1.6;
M2o = 1.6;
Mo = M1o+M2o;
Io = 0.01;
Mu = 0;
MW1 = 72;
MW2 = 104;
T = 348;
R = 1.987;
r1 = 0.27;
r2 = 0.72;
kp11 = 950;
kp12 = kp11/r1;
kp22 = (4.3e7)*exp(-7769/(R*T));
kp21 = kp22/r2;
ktc22 = (1.06e9)*exp(-1496/(R*T));
ktd11 = (1.1e8)*exp(-2800/(R*T));
ktd22 = (1.52e8)*exp(-1496/(R*T));
ktd12 = sqrt(ktd11*ktd22);
f = 0.7;
kd = (1.4e14)*exp(-30600/(R*T));
% Numerical solution
Ci0 = [Mo Mo Io Mu Mu Mu];
tvect = [0 15000];
[t,Ci] = ode45(@(t,vars)monomer(t,vars,kd,kp11,kp12,kp21,kp22,ktc22,ktd11,ktd12,ktd22,f),tvect,Ci0);
% Post-processing
M1 = Ci(:,1);
M2 = Ci(:,2);
I = Ci(:,3);
muo = Ci(:,4);
mu1 = Ci(:,5);
mu2 = Ci(:,6);
M = M1+M2;
f1 = M1./(M);
f2 = 1-f1;
M2 = M-M1;
phi1 = (kp21.*f1.*M)./((kp21.*f1.*M)+(kp12.*f2.*M));
phi2 = 1-phi1;
ktc = (ktc22.*phi2.^2);
ktd = (ktd11.*phi1.^2)+(2*ktd12.*phi1.*phi2)+(ktd22.*phi2.^2);
Yo = sqrt((2.*f.*kd.*I)./(ktc+ktd));
kp = (kp11.*phi1.*f1)+(kp12.*phi1.*f2)+(kp21.*phi2.*f1)+(kp22.*phi2.*f2);
Y1 = (kp.*M./(ktc+ktd))+Yo;
Y2 = (2.*kp.*M.*Y1./((ktc+ktd).*Yo))+(kp.*M./(ktc+ktd))+Yo;
tau = ktd.*Yo./(kp.*M);
beta = ktc.*Yo./(kp.*M);
%ODE system
function dCdT=monomer(t,vars,kd,kp11,kp12,kp21,kp22,ktc22,ktd11,ktd12,ktd22,f)
M1 = vars(1);
M2 = vars(2);
I = vars(3);
mu0 = vars(4);
mu1 = vars(5);
mu2 = vars(6);
M = M1+M2;
f1 = M1./(M);
f2 = 1-f1;
M2 = M-M1;
phi1 = (kp21.*f1.*M)./((kp21.*f1.*M)+(kp12.*f2.*M));
phi2 = 1-phi1;
ktc = (ktc22.*phi2.^2);
ktd = (ktd11*phi1.^2)+(2*ktd12*phi1*phi2)+(ktd22*phi2.^2);
Yo = sqrt((2.*f.*kd.*I)./(ktc+ktd));
kp = (kp11*phi1*f1)+(kp12*phi1*f2)+(kp21*phi2*f1)+(kp22*phi2*f2);
Y1 = (kp.*M./(ktc+ktd))+Yo;
Y2 = (2.*kp.*M.*Y1./((ktc+ktd).*Yo))+(kp.*M/(ktc+ktd))+Yo;
tau = ktd.*Yo./(kp.*M);
beta = ktc*Yo./(kp.*M);
dM1 = -kp*M1*Yo;
dM2 = -kp*M2*Yo;
dI = -kd*I;
dmu0 = 2*f*kd*I;
dmu1 = kp*M*Y1*(tau+beta);
dmu2 = kp*M*((tau*Y2)+(beta*(Y2+((Y1^2)/Yo))));
dCdT = [dM1;dM2;dI;dmu0;dmu1;dmu2];
end

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Numerical Integration and Differential Equations 的更多信息

标签

产品


版本

R2021a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by