Andycar/Matlab-Copolymerization-Modeling
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Add files via upload

Browse Source
This commit is contained in:
Georgi Tancev 2019-02-14 18:37:52 +01:00 committed by GitHub
commit 4a403e86af
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
3 changed files with 300 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

71
batch.m Normal file
View File

@ -0,0 +1,71 @@
function [ F ] = batch( t,c )
cI = c(1);
cA = c(2);
cA_R = c(3);
cB = c(4);
cB_R = c(5);
%%
V = 15; % L
cA_0 = 3.5; % mol/L
pA = 0.90; % kg/L
pB = 0.94; % kg/L
MA = 0.104; % kg/mol
MB = 0.1; % kg/mol
mA_0 = cA_0*V*MA; % kg
VA_0 = mA_0/pA; % L
mB_0 = (V-VA_0)*pB; % kg
cB_0 = mB_0/(MB*V); % mol/L
f = 0.5;
C_eta = 25;
C_RD = 180;
kd = 6.77e-6; % s^(-1)
kp_AA_0 = 4.1e2; % L/(mol s)
kp_AD_0 = 2e11; % L/(mol s)
kt_AA_0 = 2.4e7; % L/(mol s)
kt_AD_0 = 5e8; % L/(mol s)
kp_BB_0 = 9.3e2; % L/(mol s)
kp_BD_0 = 2e11; % L/(mol s)
kt_BB_0 = 9.2e6; % L/(mol s)
kt_BD_0 = 5e8; % L/(mol s)
rA = 0.52;
rB = 0.46;
%%
wp = (MA*V*(cA_0-c(2))+MB*V*(cB_0-c(4)))/(mA_0+mB_0);
kp_AA = 1/(1/kp_AA_0 + exp(C_eta*wp)/kp_AD_0);
kt_AA = 1/(1/kt_AA_0 + exp(C_eta*wp)/kt_AD_0) + C_RD*kp_AA*(1-wp);
kp_BB = 1/(1/kp_BB_0 + exp(C_eta*wp)/kp_BD_0);
kt_BB = 1/(1/kt_BB_0 + exp(C_eta*wp)/kt_BD_0) + C_RD*kp_BB*(1-wp);
kp_AB = kp_AA/rA;
kt_AB = sqrt(kt_AA*kt_BB);
kp_BA = kp_BB/rB;
kt_BA = kt_AB;
%%
dcI = - kd*cI;
dcA = - kp_AA*cA*cA_R - kp_BA*cA*cB_R;
dcA_R = + 2*f*kd*cI*(cA/(cA+cB)) - kt_AA*cA_R^2 - kt_AB*cA_R*cB_R - kp_AB*cA_R*cB + kp_BA*cB_R*cA;
dcB = - kp_BB*cB*cB_R - kp_AB*cB*cA_R;
dcB_R = + 2*f*kd*cI*(cB/(cA+cB)) - kt_BB*cB_R^2 - kt_BA*cA_R*cB_R - kp_BA*cB_R*cA + kp_AB*cA_R*cB;
%%
F = [dcI; dcA; dcA_R; dcB; dcB_R];
end

112
cumulative.m Normal file
View File

@ -0,0 +1,112 @@
function [ D ] = cumulative( t,c )
%% data
V = 15; % L
cA_0 = 3.5; % mol/L
pA = 0.90; % kg/L
pB = 0.94; % kg/L
MA = 0.104; % kg/mol
MB = 0.1; % kg/mol
mA_0 = cA_0*V*MA; % kg
VA_0 = mA_0/pA; % L
mB_0 = (V-VA_0)*pB; % kg
cB_0 = mB_0/(MB*V); % mol/L
f = 0.5;
C_eta = 25;
C_RD = 180;
kd = 6.77e-6; % s^(-1)
kp_AA_0 = 4.1e2; % L/(mol s)
kp_AD_0 = 2e11; % L/(mol s)
kt_AA_0 = 2.4e7; % L/(mol s)
kt_AD_0 = 5e8; % L/(mol s)
kp_BB_0 = 9.3e2; % L/(mol s)
kp_BD_0 = 2e11; % L/(mol s)
kt_BB_0 = 9.2e6; % L/(mol s)
kt_BD_0 = 5e8; % L/(mol s)
rA = 0.52;
rB = 0.46;
C_t = 1000;
%%
h = length(t);
for t1=1:1:h;
c_A = c(t1,2);
c_A_dot = c(t1,3);
c_B = c(t1,4);
c_B_dot = c(t1,5);
c_R_dot = c_A_dot+c_B_dot;
wp = (MA*V*(cA_0-c_A)+MB*V*(cB_0-c_B))/(mA_0+mB_0);
kp_AA = 1/(1/kp_AA_0 + exp(C_eta*wp)/kp_AD_0);
kt_AA = 1/(1/kt_AA_0 + exp(C_eta*wp)/kt_AD_0) + C_RD*kp_AA*(1-wp);
kp_BB = 1/(1/kp_BB_0 + exp(C_eta*wp)/kp_BD_0);
kt_BB = 1/(1/kt_BB_0 + exp(C_eta*wp)/kt_BD_0) + C_RD*kp_BB*(1-wp);
kp_AB = kp_AA/rA;
kt_AB = sqrt(kt_AA*kt_BB);
kp_BA = kp_BB/rB;
kt_BA = kt_AB;
p_A = c_A_dot/c_R_dot;
p_B = c_B_dot/c_R_dot;
kpM = (kp_AA*p_A+kp_BA*p_B)*c_A+(kp_BB*p_B+kp_AB*p_A)*c_B;
kt = kt_AA*p_A^2+2*kt_AB*p_A*p_B+kt_BB*p_B^2;
ktc = kt/(1+C_t);
ktd = kt-ktc;
beta = ktc*c_R_dot/(kpM);
gamma = ktd*c_R_dot/(kpM);
alph = beta + gamma;
R_p_A = (kp_AA*p_A+kp_BA*p_B)*c_A*c_R_dot;
R_p_B = (kp_BB*p_B+kp_AB*p_A)*c_B*c_R_dot;
R_p = R_p_A+R_p_B;
dPdt(t1) = R_p*(gamma+beta/2);
n_N(t1) = 1/(gamma+beta/2);
n_W(t1) = 2*(gamma+1.5*beta)/alph^2;
S1(t1) = n_W(t1)/n_N(t1);
mu1(t1) = n_N(t1)*dPdt(t1);
mu2(t1) = n_W(t1)*mu1(t1);
F_A_1(t1) = (((rA*c_A+c_B)*c_A)/((rA*c_A+c_B)*c_A+(rB*c_B+c_A)*c_B));
F_A(t1) = (((rA*c_A+c_B)*c_A)/((rA*c_A+c_B)*c_A+(rB*c_B+c_A)*c_B))*dPdt(t1);
end
tspan = 0:1:length(t);
for t1=2:1:length(t)-1;
P_t(t1) = trapz(tspan(2:t1+1),dPdt(2:t1+1));
F_A_c(t1) = trapz(tspan(2:t1+1),F_A(2:t1+1))/P_t(t1);
mu1c(t1) = trapz(tspan(2:t1+1),mu1(2:t1+1))/P_t(t1);
mu2c(t1) = trapz(tspan(2:t1+1),mu2(2:t1+1))/P_t(t1);
n_N_c(t1) = mu1c(t1);
n_W_c(t1) = mu2c(t1)/mu1c(t1);
S2(t1) = mu2c(t1)/(mu1c(t1)^2);
end
D = [F_A_1(1:end-1)' F_A_c' n_N_c' n_W_c' S2' n_N(1:end-1)' n_W(1:end-1)' S1(1:end-1)'];
end

117
main.m Normal file
View File

@ -0,0 +1,117 @@
%% Polymerization Reaction Model
% George Tancev
clear all; close all; clc;
%% data
V = 15; % L
cA_0 = 3.5; % mol/L
pA = 0.90; % kg/L
pB = 0.94; % kg/L
wI_0 = 0.01;
MA = 0.104; % kg/mol
MB = 0.1; % kg/mol
MI = 0.164; % kg/mol
mA_0 = cA_0*V*MA; % kg
VA_0 = mA_0/pA; % L
mB_0 = (V-VA_0)*pB; % kg
cB_0 = mB_0/(MB*V); % mol/L
X_A_0 = cA_0/(cA_0+cB_0);
mI_0 = wI_0*(mA_0+mB_0); % kg
cI_0 = mI_0/(MI*V); % mol/L
f = 0.5;
C_eta = 25;
C_RD = 180;
kd = 6.77e-6; % s^(-1)
kp_AA_0 = 4.1e2; % L/(mol s)
kp_AD_0 = 2e11; % L/(mol s)
kt_AA_0 = 2.4e7; % L/(mol s)
kt_AD_0 = 5e8; % L/(mol s)
kp_BB_0 = 9.3e2; % L/(mol s)
kp_BD_0 = 2e11; % L/(mol s)
kt_BB_0 = 9.2e6; % L/(mol s)
kt_BD_0 = 5e8; % L/(mol s)
rA = 0.52;
rB = 0.46;
kp_AB_0 = kp_AA_0/rA; % L/(mol s)
kp_BA_0 = kp_BB_0/rB; % L/(mol s)
%% Solving the system of ODEs
tspan = 1:1:7200; % s
[t,c] = ode15s(@(t,c)batch(t,c),tspan,[cI_0 cA_0 0 cB_0 0]);
t = t/3600;
%% a)
X_A = linspace(0,1,100);
F_A = ((rA-1).*X_A.^2+X_A)./((rA-2).*X_A.^2+2.*X_A+rB.*(1-X_A).^2);
DIAG = X_A;
x = 1-(c(:,2)+c(:,4))/(cA_0+cB_0);
F_A_c = cumulative( t,c );
figure(1);
subplot(3,2,1);
plot(X_A,F_A,X_A,DIAG,'--');
title('Mayo-Lewis diagram');
xlabel({'X_A'});
ylabel({'F_A'});
figure(1);
subplot(3,2,2);
plot(x(1:end-1),F_A_c(1:end,1),x(3:end),F_A_c(2:end,2));
axis([0 1 0 0.6])
title('cumulative composition distribution');
xlabel({'conversion','x'});
ylabel({'F_A, F_A^c'});
legend('F_A','F_A^c','Location','best');
%% b)
figure(1);
subplot(3,2,3);
plot(t,c(:,2),t,c(:,4));
axis([0 2 0 6]);
title('concentration profile');
xlabel({'time','h'});
ylabel({'concentration','mol / L'});
legend('styrene','methyl methacrylate','Location','best');
figure(1);
subplot(3,2,4);
plot(t,x);
title('conversion profile');
xlabel({'time','h'});
ylabel({'conversion','x'});
%% c)
figure(1);
subplot(3,2,5);
plot(x(3:end),F_A_c(2:end,3),x(3:end),F_A_c(2:end,4),x(3:end),F_A_c(2:end,6),x(3:end),F_A_c(2:end,7));
title({'cumulative and instantanenous', 'number/ weight average'});
xlabel({'conversion','x'});
ylabel({'n_N, n_W'});
legend('n^c_N','n^c_W','n_N','n_W','Location','best');
figure(1);
subplot(3,2,6);
plot(x(3:end),F_A_c(2:end,5),x(3:end),F_A_c(2:end,8));
title({'cumulative and instantaneous', 'polydispersity'});
xlabel({'conversion','x'});
ylabel({'\sigma'});
legend('\sigma^c','\sigma','Location','best');
%%