I tried to establish a dynamic model of the reactor using PDEs, but I couldn’t solve it. I don’t know why. The code is as follows：
In this problem, my boundary condition only has the left boundary, which is the u value and dudx when z=0. How should I write it?
function reactor_dynamic_model(obj)
global R miu comp_set A_for A_rev Ea_for Ea_rev Da lambda
Da = 1.8e-4;
lambda = 5e-4;
F_molar = obj.F_molar;
A = obj.cat_vol / obj.z0;
z0 = obj.z0;
y0 = obj.y0;
epsilon = obj.epsilon;
P = obj.P;
T0 = obj.T;
Fi0 = F_molar * y0;
Cpc = obj.Cpc;
rho = obj.rho;
epsilon = obj.epsilon;
nfe = obj.nfe;
ct0 = 0.001/pr(T0,P,y0);
Cpm0 = get_Cpm(T0,y0);
ci0=ct0.*y0;
uv0 = (epsilon*ct0.*Cpm0+Cpc*rho).*T0;
zspan = 0:z0/(nfe-1):z0;
tspan = 0:10:100;
sol = pdepe(0,@pdefun,@pdeic,@pdebc,zspan,tspan);
a=1;
function [c,f,s] = pdefun(z,t,u,dudx)
c1 = u(1);
c2 = u(2);
c3 = u(3);
c4 = u(4);
uv = u(5);
ct = c1+c2+c3+c4;
y1=c1/ct;
y2=c2/ct;
y3=c3/ct;
y4=c4/ct;
y = [y1,y2,y3,y4];
T = fzero(@solve_T_uv, T0,[],y,uv,ct);
function f0 = solve_T_uv(T,y,uv,ct)
f0 = (epsilon*ct*get_Cpm(T,y)+Cpc*rho)*T-uv;
end
dc1dz = dudx(1);
dc2dz = dudx(2);
dc3dz = dudx(3);
dc4dz = dudx(4);
dTdz = dudx(5);
x = (Fi0(1)-F_molar*y1) / (1-2*y1);
F1 = Fi0(1)-x;
F2 = Fi0(2) -3*x;
F3 = Fi0(3)+2*x;
F4 = Fi0(4);
Ft = F1+F2+F3+F4;
Cpm = get_Cpm(T, y);
h = Ft*Cpm*T;
X = (Fi0(1) – F1) / Fi0(1);
eta = get_eta(T, P, X);
r = get_r(T, P, y, eta, miu,A_for, A_rev, Ea_for, Ea_rev);
r1 = r(1);
r2 = r(2);
r3 = r(3);
r4 = r(4);
Hr = get_Hr(T, P);
c = [1;1;1;1;1];
f = [Da*dc1dz/epsilon;
Da*dc2dz/epsilon;
Da*dc3dz/epsilon;
Da*dc4dz/epsilon;
lambda*dTdz];
s = [r1/epsilon-F1/A/epsilon;
r2/epsilon-F2/A/epsilon;
r3/epsilon-F3/A/epsilon;
r4/epsilon-F4/A/epsilon;
r3*(-Hr)-h/A];
end
function u0 = pdeic(z)
class1 = PFR;
result = class1.reactor_steadystate_model(z);
T_init = result(1);
y_init = result(10:13);
ct_init = 0.001/pr(T_init,P,y_init);
ci_init = ct_init * y_init;
c1_init = ci_init(1);
c2_init = ci_init(2);
c3_init = ci_init(3);
c4_init = ci_init(4);
Cpm_init = get_Cpm(T_init,y_init);
uv_init = (epsilon*ct_init*Cpm_init+Cpc*rho)*T_init;
u0 = [c1_init;c2_init;c3_init;c4_init;uv_init];
end

function [pl,ql,pr,qr] = pdebc(zl,ul,zr,ur,t)
pl = [ul(1)-ci0(1);ul(2)-ci0(2);ul(3)-ci0(3);ul(4)-ci0(4);ul(5)-uv0];
ql = [0;0;0;0;0];
pr = [0;0;0;0;0];
qr = [0;0;0;0;0];
end
endI tried to establish a dynamic model of the reactor using PDEs, but I couldn’t solve it. I don’t know why. The code is as follows：
In this problem, my boundary condition only has the left boundary, which is the u value and dudx when z=0. How should I write it?
function reactor_dynamic_model(obj)
global R miu comp_set A_for A_rev Ea_for Ea_rev Da lambda
Da = 1.8e-4;
lambda = 5e-4;
F_molar = obj.F_molar;
A = obj.cat_vol / obj.z0;
z0 = obj.z0;
y0 = obj.y0;
epsilon = obj.epsilon;
P = obj.P;
T0 = obj.T;
Fi0 = F_molar * y0;
Cpc = obj.Cpc;
rho = obj.rho;
epsilon = obj.epsilon;
nfe = obj.nfe;
ct0 = 0.001/pr(T0,P,y0);
Cpm0 = get_Cpm(T0,y0);
ci0=ct0.*y0;
uv0 = (epsilon*ct0.*Cpm0+Cpc*rho).*T0;
zspan = 0:z0/(nfe-1):z0;
tspan = 0:10:100;
sol = pdepe(0,@pdefun,@pdeic,@pdebc,zspan,tspan);
a=1;
function [c,f,s] = pdefun(z,t,u,dudx)
c1 = u(1);
c2 = u(2);
c3 = u(3);
c4 = u(4);
uv = u(5);
ct = c1+c2+c3+c4;
y1=c1/ct;
y2=c2/ct;
y3=c3/ct;
y4=c4/ct;
y = [y1,y2,y3,y4];
T = fzero(@solve_T_uv, T0,[],y,uv,ct);
function f0 = solve_T_uv(T,y,uv,ct)
f0 = (epsilon*ct*get_Cpm(T,y)+Cpc*rho)*T-uv;
end
dc1dz = dudx(1);
dc2dz = dudx(2);
dc3dz = dudx(3);
dc4dz = dudx(4);
dTdz = dudx(5);
x = (Fi0(1)-F_molar*y1) / (1-2*y1);
F1 = Fi0(1)-x;
F2 = Fi0(2) -3*x;
F3 = Fi0(3)+2*x;
F4 = Fi0(4);
Ft = F1+F2+F3+F4;
Cpm = get_Cpm(T, y);
h = Ft*Cpm*T;
X = (Fi0(1) – F1) / Fi0(1);
eta = get_eta(T, P, X);
r = get_r(T, P, y, eta, miu,A_for, A_rev, Ea_for, Ea_rev);
r1 = r(1);
r2 = r(2);
r3 = r(3);
r4 = r(4);
Hr = get_Hr(T, P);
c = [1;1;1;1;1];
f = [Da*dc1dz/epsilon;
Da*dc2dz/epsilon;
Da*dc3dz/epsilon;
Da*dc4dz/epsilon;
lambda*dTdz];
s = [r1/epsilon-F1/A/epsilon;
r2/epsilon-F2/A/epsilon;
r3/epsilon-F3/A/epsilon;
r4/epsilon-F4/A/epsilon;
r3*(-Hr)-h/A];
end
function u0 = pdeic(z)
class1 = PFR;
result = class1.reactor_steadystate_model(z);
T_init = result(1);
y_init = result(10:13);
ct_init = 0.001/pr(T_init,P,y_init);
ci_init = ct_init * y_init;
c1_init = ci_init(1);
c2_init = ci_init(2);
c3_init = ci_init(3);
c4_init = ci_init(4);
Cpm_init = get_Cpm(T_init,y_init);
uv_init = (epsilon*ct_init*Cpm_init+Cpc*rho)*T_init;
u0 = [c1_init;c2_init;c3_init;c4_init;uv_init];
end

function [pl,ql,pr,qr] = pdebc(zl,ul,zr,ur,t)
pl = [ul(1)-ci0(1);ul(2)-ci0(2);ul(3)-ci0(3);ul(4)-ci0(4);ul(5)-uv0];
ql = [0;0;0;0;0];
pr = [0;0;0;0;0];
qr = [0;0;0;0;0];
end
end I tried to establish a dynamic model of the reactor using PDEs, but I couldn’t solve it. I don’t know why. The code is as follows：
In this problem, my boundary condition only has the left boundary, which is the u value and dudx when z=0. How should I write it?
function reactor_dynamic_model(obj)
global R miu comp_set A_for A_rev Ea_for Ea_rev Da lambda
Da = 1.8e-4;
lambda = 5e-4;
F_molar = obj.F_molar;
A = obj.cat_vol / obj.z0;
z0 = obj.z0;
y0 = obj.y0;
epsilon = obj.epsilon;
P = obj.P;
T0 = obj.T;
Fi0 = F_molar * y0;
Cpc = obj.Cpc;
rho = obj.rho;
epsilon = obj.epsilon;
nfe = obj.nfe;
ct0 = 0.001/pr(T0,P,y0);
Cpm0 = get_Cpm(T0,y0);
ci0=ct0.*y0;
uv0 = (epsilon*ct0.*Cpm0+Cpc*rho).*T0;
zspan = 0:z0/(nfe-1):z0;
tspan = 0:10:100;
sol = pdepe(0,@pdefun,@pdeic,@pdebc,zspan,tspan);
a=1;
function [c,f,s] = pdefun(z,t,u,dudx)
c1 = u(1);
c2 = u(2);
c3 = u(3);
c4 = u(4);
uv = u(5);
ct = c1+c2+c3+c4;
y1=c1/ct;
y2=c2/ct;
y3=c3/ct;
y4=c4/ct;
y = [y1,y2,y3,y4];
T = fzero(@solve_T_uv, T0,[],y,uv,ct);
function f0 = solve_T_uv(T,y,uv,ct)
f0 = (epsilon*ct*get_Cpm(T,y)+Cpc*rho)*T-uv;
end
dc1dz = dudx(1);
dc2dz = dudx(2);
dc3dz = dudx(3);
dc4dz = dudx(4);
dTdz = dudx(5);
x = (Fi0(1)-F_molar*y1) / (1-2*y1);
F1 = Fi0(1)-x;
F2 = Fi0(2) -3*x;
F3 = Fi0(3)+2*x;
F4 = Fi0(4);
Ft = F1+F2+F3+F4;
Cpm = get_Cpm(T, y);
h = Ft*Cpm*T;
X = (Fi0(1) – F1) / Fi0(1);
eta = get_eta(T, P, X);
r = get_r(T, P, y, eta, miu,A_for, A_rev, Ea_for, Ea_rev);
r1 = r(1);
r2 = r(2);
r3 = r(3);
r4 = r(4);
Hr = get_Hr(T, P);
c = [1;1;1;1;1];
f = [Da*dc1dz/epsilon;
Da*dc2dz/epsilon;
Da*dc3dz/epsilon;
Da*dc4dz/epsilon;
lambda*dTdz];
s = [r1/epsilon-F1/A/epsilon;
r2/epsilon-F2/A/epsilon;
r3/epsilon-F3/A/epsilon;
r4/epsilon-F4/A/epsilon;
r3*(-Hr)-h/A];
end
function u0 = pdeic(z)
class1 = PFR;
result = class1.reactor_steadystate_model(z);
T_init = result(1);
y_init = result(10:13);
ct_init = 0.001/pr(T_init,P,y_init);
ci_init = ct_init * y_init;
c1_init = ci_init(1);
c2_init = ci_init(2);
c3_init = ci_init(3);
c4_init = ci_init(4);
Cpm_init = get_Cpm(T_init,y_init);
uv_init = (epsilon*ct_init*Cpm_init+Cpc*rho)*T_init;
u0 = [c1_init;c2_init;c3_init;c4_init;uv_init];
end

function [pl,ql,pr,qr] = pdebc(zl,ul,zr,ur,t)
pl = [ul(1)-ci0(1);ul(2)-ci0(2);ul(3)-ci0(3);ul(4)-ci0(4);ul(5)-uv0];
ql = [0;0;0;0;0];
pr = [0;0;0;0;0];
qr = [0;0;0;0;0];
end
end pdepe, boundary condition, spatial discretization, elliptic, pdes, initial condition, pfr, dynamic model, ammonia synthesis MATLAB Answers — New Questions

​