global PMMA_parameter ;
%%%%%%%%%% Kinetic parameters used in simulation %%%%%%%%%
theta = 50*60
T = 86+273.15 %K
R = 1.987 %Cal/mol*k
kd = (1.69*10^14)*exp(-30000/(R*T)) %(s^-1)
kp0 = (4.92*10^5)*exp(-4353/(R*T)) %(L/=gmol*s)
kt0 = (9.8*10^7)*exp(-701/(R*T)) %(L/=gmol*s)
f = 0.5 ; %(無單位)
kf = 0.043 ; %(L/=gmol*s)
kfs = 0.091 ;%(L/=gmol*s)
PMMA_parameter.theta = theta ;
PMMA_parameter.T = T ;
PMMA_parameter.R = R ;
PMMA_parameter.kd = kd ;
PMMA_parameter.kp0 = kp0 ;
PMMA_parameter.kt0 = kt0 ;
PMMA_parameter.f = f ;
PMMA_parameter.kf = kf ;
PMMA_parameter.kfs = kfs ;
%%%%%%%%%% gel effect %%%%%%%%%
%alphp = 0.00048 ; %(無單位)
%alphm = 0.001 ; %(無單位)
%alphs = 0.001 ; %(無單位)
%Tgp = 387 ; %K
%Tgm = 167 ; %K
%Tgs = 181 ; %K
Vfp = 0.025+0.00048*(T-387) ;
Vfm = 0.025+0.001*(T-167) ;
Vfs = 0.025+0.001*(T-181) ;
PMMA_parameter.Vfp = Vfp ;
PMMA_parameter.Vfm = Vfm ;
PMMA_parameter.Vfs = Vfs ;
rhoMMA = 0.9654-0.00109*T-(9.7*10^-7)*(T^2) ; %gm/ml 1984 schmidt
rhoI = 1.334; %g/L
PMMA_parameter.rhoMMA = rhoMMA;
PMMA_parameter.rhoI = rhoI;
%%%%%%%%%% initial concentration %%%%%%%%%
If = 0.0032 ; %mole/l
%Sf = 0.13166 ;
Ii = If/(1+theta*kd)
PMMA_parameter.If = If ;
PMMA_parameter. Ii = Ii ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rhoPMMA = rhoMMA/(0.754-(9*10^-4)*(T-70)); %gm/ml
Wmma = 100.121 ; %g/mol
PMMA_parameter.rhoPMMA = rhoPMMA ;
PMMA_parameter.Wmma = Wmma ;
[t,x] = ode45 ('pmma3',[0.1 10],[0.002499057 3.039996046719E-01 0.131657143 0 0])
以上是先建立一個參數檔
以下是方程式檔
function schmidt_pmma = pmma3(t,x)
global PMMA_parameter;
fim = (x(2)*100.121/rhoMMA)/100
fis = 0.3
fii = (x(1)*242.23/rhoI)/100
fip = 1-fim-fis-fii
Vf = Vfp*fip+Vfm*fim+Vfs*fis
if Vf >= (0.1856-2.965*10^(-4)*(T-273.2))
gt = 0.10575*exp(17.15*Vf-0.01715*(T-273.2))
else
gt = (2.3*10^(-6))*exp(75*Vf)
end
kt = gt*kt0
if Vf >= 0.05
gp = 1
else
gp = 7.1*10^(-5)*exp(171.53*Vf)
end
kp = gp*kp0
ktd = kt/(1+1/8.32)
ktc = ktd/8.32
P0 = ((2*f*kd*Ii)/kt)^1/2
%Sf = Si ;
Mf = If*95 ; %Mf/If=95
Mi = Mf/(1+theta*kp*P0)
M0 = Mi;
epsilonp = (M0*Wmma/1000)*((1/rhoPMMA)-(1/rhoMMA))
%%%%%%%%%% polymer mass concentration in the reactor %%%%%%%%%
Cpol = Wmma*x(5)
%%%%%%%%%% solvent-free mass fraction %%%%%%%%%
Xp = Cpol/(Cpol+Wmma*x(2))
P = (((2*f*kd*x(1))/kt)^1/2)
alph = kp*x(2)/(kp*x(2)+kfs*x(3)+kf*x(2)+kt*P)
schmidt_pmma = [
(1/theta)*(If-((1+epsilonp*Xp)*x(1)))-(kd*x(1)); %I
(1/theta)*(Mf-((1+epsilonp*Xp)*x(2)))-(kp*x(2)*P); %M
(1/theta)*(Sf-((1+epsilonp*Xp)*x(3))); %S
-(1/theta)*(1+epsilonp*Xp)*x(4)+((kfs*x(3)+kf*x(2)+ktd*P)*alph*P)+(1/theta)*(1/2*ktc*(P^2)); %Lambda0
-(1/theta)*(1+epsilonp*Xp)*x(5)+((P/(1-alph))*((kfs*x(3)+kf*x(2)+ktd*P)*alph*(2-alph)+ktc*P))]; %Lambda1
end
結果為
t =
0.1000
0.3475
0.5950
0.8425
1.0900
1.3375
1.5850
1.8325
2.0800
2.3275
2.5750
2.8225
3.0700
3.3175
3.5650
3.8125
4.0600
4.3075
4.5550
4.8025
5.0500
5.2975
5.5450
5.7925
6.0400
6.2875
6.5350
6.7825
7.0300
7.2775
7.5250
7.7725
8.0200
8.2675
8.5150
8.7625
9.0100
9.2575
9.5050
9.7525
10.0000
x =
0.0025 0.3040 0.1317 0 0
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
求大老解答
请把能运行的代码粘贴一下