From 4a403e86af4acef2fd8e77ffd092c8f8df3f0e16 Mon Sep 17 00:00:00 2001 From: Georgi Tancev <30585013+gtancev@users.noreply.github.com> Date: Thu, 14 Feb 2019 18:37:52 +0100 Subject: [PATCH] Add files via upload --- batch.m | 71 +++++++++++++++++++++++++++++++ cumulative.m | 112 ++++++++++++++++++++++++++++++++++++++++++++++++ main.m | 117 +++++++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 300 insertions(+) create mode 100644 batch.m create mode 100644 cumulative.m create mode 100644 main.m diff --git a/batch.m b/batch.m new file mode 100644 index 0000000..acc101a --- /dev/null +++ b/batch.m @@ -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 + diff --git a/cumulative.m b/cumulative.m new file mode 100644 index 0000000..ef3b3a0 --- /dev/null +++ b/cumulative.m @@ -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 \ No newline at end of file diff --git a/main.m b/main.m new file mode 100644 index 0000000..22ae9ac --- /dev/null +++ b/main.m @@ -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'); + +%%