Having problem to simulate polymerization model using ODE

Hello.I'm trying to simulate two different model for free radical polymerization. I've successfully simulated the first model but I need guidelines to solve the second model which is more difficult. Here is the matlab mfile coding for the first model.

1.function mfile:

function Y=CCSmodel(t,x) global T Io %variables X = x(1); V = x(2); M = x(3); I = x(4); S = x(5); Lo = x(6); L1 = x(7); L2 = x(8); uo = x(9); u1 = x(10); u2 = x(11);

%volume contraction : e=(dp-dm)/dp dp = 1084.8 -0.605*(T-273.15); dm = 924 - 0.918*(T-273.15); e = (dp-dm)/dp;

%initiator rate constant & efficiency kd = 1.58e15*exp((-15508)./T)*60; f = 89.81*exp((-1713)/T);

%volume fraction fs = 0.0; beta = fs/(1-fs);

%% R = 1.98722; Bi = 0.03; Tgp = 378; Zp = (X*(1-e))./(1-(e*X)+beta); A = 0.168-(8.21e-6*(((T-273.15)-(Tgp-273.15))^2)); C = exp((2.303*(1-Zp))/(A+(Bi*(1-Zp))));

kp0 = 1.051e7*exp((-3557)/T)*60; kt0 = 1.255e9*exp((-844)/T)*60;

%gel effect equation teta_T = ((1.1353e-22)/Io)*exp(17420/T); kt = (kt0*C)/(C+(teta_T*kt0*Lo)); ktd = 0; ktc = kt-ktd;

%glass effect equation teta_P = 5.4814e-16*exp(13982/T); kp = (kp0*C)/(C+(teta_P*kp0*Lo));

%chain transfer to monomer kf = 6.128e8*exp(-13450/(R*T)); %chain transfer to solvent ks = 3.148e9*exp(-16264/(R*T))*60;

%% rI = kd*I; rLo = 2*f*kd*I - kt*Lo.^2; rL1 = 2*f*kd*I - kt*Lo.*L1 + ((ks*S+kf*M).*(Lo-L1)) + (kp*M.*Lo); rL2 = 2*f*kd*I - kt*Lo.*L2 + (ks*S.*(Lo-L2)) + (kp*M.*(2*L1+Lo)) + (kf*M.*(Lo-L2)); ruo = ktd*Lo.^2 + 0.5*ktc*Lo.^2+kf*M.*Lo+ks*S.*Lo; ru1 = kt*(Lo.*L1)+kf*M.*L1+ks*S.*L1; ru2 = kt*(Lo.*L2) + ktc*(L1.^2)+kf*M.*L2+ks*S.*L2;

%%%%%%%%%list of ODE%%%%%%% %balances dXdt = (kp+kf)*(1-X)*Lo; dVdt = (-V*e*dXdt)/(1+beta); dMdt = -kp*M.*Lo - (M./V).*dVdt; dIdt = - rI - (I./V).*dVdt; dSdt = -(ks*S*Lo)-(S./V).*dVdt; % moment of live radical distribution dLodt = rLo - (Lo/V).*dVdt; dL1dt = rL1 - (L1/V).*dVdt; dL2dt = rL2 - (L2/V).*dVdt; % moment of dead polymer distribution duodt = ruo - (uo/V).*dVdt; du1dt = ru1 - (u1/V).*dVdt; du2dt = ru2 - (u2/V).*dVdt;

Y=[dXdt;dVdt;dMdt;dIdt;dSdt;dLodt;dL1dt;dL2dt;duodt;du1dt;du2dt];

2.ODE task mfile:

close all clear all clc

global T Io Tc = 50; %changeable T = Tc+273.15; Io = 0.0258; Vo = 1.2; Mo = 8.85; So = 3.60; %% xo = [0 Vo Mo Io So 0 0 0 0 0 0]; TSPAN = (0:0.1:1400); [t,x] = ode15s(@CCSmodel,TSPAN,xo);

%% average MW % Mn = MW(L1+u1)/(Lo+uo) & Mw = MW(L2+u2)/(L1+u1) MW = 104.14; Mn = ((x(:,7)+x(:,10))./(x(:,6)+x(:,9))).*MW; %number average molecular weight Mw = ((x(:,8)+x(:,11))./(x(:,7)+x(:,10))).*MW; %weight average molecular weight

For the second model (free volume theory),the difference is for the gel effect equation which is depends on:

i) Kcr=9.44exp(3833/(RT)) ii) K=((Mw)^m)exp(A/Vf) where Vf=Vfm+Vfp+Vfs

when K is less than Kcr;

kt=kt0(1+delta*cp*MW)

when K becomes equal to Kcr, Vf and Mw are taken as being their critical values and;

kt=kt0((Mwcr/Mw)^n)exp(-A((1/Vf)-(1/Vfcr)))

All the parameters given above have their own specific values. Mwcr and Vfcr are the critical values for Mw and Vf.

Mw used in this model is the weight average molecular weight where based on the first model, I calculated it after getting ODE's results. I'm aware that we can use command if to solve equation with conditions but I dont know how to make Mw as a part of function mfile and set the critical value. I tried using global but doesnt work.

Any help for me to solve this would be greatly appreciated. Thank you

[aza]

I dont know why I can't directly paste the coding, but here is the correct one

1. function mfie

function Y=CCSmodel(t,x)

global T Io

%variables

X = x(1);

V = x(2);

M = x(3);

I = x(4);

S = x(5);

Lo = x(6);

L1 = x(7);

L2 = x(8);

uo = x(9);

u1 = x(10);

u2 = x(11);

%volume contraction : e=(dp-dm)/dp

dp = 1084.8 -0.605*(T-273.15);

dm = 924 - 0.918*(T-273.15);

e = (dp-dm)/dp;

%initiator rate constant & efficiency

kd = 1.58e15*exp((-15508)./T)*60;

f = 89.81*exp((-1713)/T);

%volume fraction

fs = 0.0;

beta = fs/(1-fs);

%%

R = 1.98722;

Bi = 0.03;

Tgp = 378;

Zp = (X*(1-e))./(1-(e*X)+beta);

A = 0.168-(8.21e-6*(((T-273.15)-(Tgp-273.15))^2));

C = exp((2.303*(1-Zp))/(A+(Bi*(1-Zp))));

kp0 = 1.051e7*exp((-3557)/T)*60;

kt0 = 1.255e9*exp((-844)/T)*60;

%gel effect equation

teta_T = ((1.1353e-22)/Io)*exp(17420/T);

kt = (kt0*C)/(C+(teta_T*kt0*Lo));

ktd = 0;

ktc = kt-ktd;

%glass effect equation

teta_P = 5.4814e-16*exp(13982/T);

kp = (kp0*C)/(C+(teta_P*kp0*Lo));

%chain transfer to monomer

kf = 6.128e8*exp(-13450/(R*T));

%chain transfer to solvent

ks = 3.148e9*exp(-16264/(R*T))*60;

%%

rI = kd*I;

rLo = 2*f*kd*I - kt*Lo.^2;

rL1 = 2*f*kd*I - kt*Lo.*L1 + ((ks*S+kf*M).*(Lo-L1)) + (kp*M.*Lo);

rL2 = 2*f*kd*I - kt*Lo.*L2 + (ks*S.*(Lo-L2)) + (kp*M.*(2*L1+Lo)) + (kf*M.*(Lo-L2));

ruo = ktd*Lo.^2 + 0.5*ktc*Lo.^2+kf*M.*Lo+ks*S.*Lo;

ru1 = kt*(Lo.*L1)+kf*M.*L1+ks*S.*L1;

ru2 = kt*(Lo.*L2) + ktc*(L1.^2)+kf*M.*L2+ks*S.*L2;

%%%%%%%%%list of ODE%%%%%%%

%balances

dXdt = (kp+kf)*(1-X)*Lo;

dVdt = (-V*e*dXdt)/(1+beta);

dMdt = -kp*M.*Lo - (M./V).*dVdt;

dIdt = - rI - (I./V).*dVdt;

dSdt = -(ks*S*Lo)-(S./V).*dVdt;

% moment of live radical distribution

dLodt = rLo - (Lo/V).*dVdt;

dL1dt = rL1 - (L1/V).*dVdt;

dL2dt = rL2 - (L2/V).*dVdt;

% moment of dead polymer distribution

duodt = ruo - (uo/V).*dVdt;

du1dt = ru1 - (u1/V).*dVdt;

du2dt = ru2 - (u2/V).*dVdt;

Y=[dXdt;dVdt;dMdt;dIdt;dSdt;dLodt;dL1dt;dL2dt;duodt;du1dt;du2dt];

2. ODE task mfile:

close all

clear all

clc

global T Io

Tc = 50; %changeable

T = Tc+273.15;

Io = 0.0258;

Vo = 1.2;

Mo = 8.85;

So = 3.60;

%%

xo = [0 Vo Mo Io So 0 0 0 0 0 0];

TSPAN = (0:0.1:1400);

[t,x] = ode15s(@CCSmodel,TSPAN,xo);

%% average MW

% Mn = MW(L1+u1)/(Lo+uo) & Mw = MW(L2+u2)/(L1+u1)

MW = 104.14;

Mn = ((x(:,7)+x(:,10))./(x(:,6)+x(:,9))).*MW; %number average molecular weight

Mw = ((x(:,8)+x(:,11))./(x(:,7)+x(:,10))).*MW; %weight average molecular weight