i need to get my variable y to vary....as its a key variable for this entire system. i cant seem to find a way. would you help me in that regard...

   clc
clear
clf

%Initial values
x0 = [300; 700; 0; 200; 0; 20; 1; 423]; %Iterate initial values
%Initial values:
%nNO_0 = 300 mol
%nH2_0 = 700 mol
%nN2_0 = 0 mol
%nH2O_0 = 0 mol
%nN2O_0 = 0 mol
%nO2_0 = 20 mol
%P_0 = 101325 Pa
%T_0 = 423 K
%Integrate ODE's
[W, x] = ode45(@(W,x) PBR(W,x), [0 1000], x0);

%Plot results
%Create data and 2-by-1 tiled chart layout



plot(W,x(:,1),'LineWidth',2)
title('NO moles VS catalyst weight');
xlabel('W (g)');
ylabel('NO (mol)');
plot(W,x(:,2),'LineWidth',2);
title('H2 moles VS catalyst weight');
xlabel('W (g)');
ylabel('H2 (mol)');
shg
plot(W,x(:,3),'LineWidth',2);
title('N2 moles VS catalyst weight');
xlabel('W (g)');
ylabel('N2 (mol)');
shg
plot(W,x(:,4),'LineWidth',2);
title('H2O moles VS catalyst weight');
xlabel('W (g)');
ylabel('H2O (mol)');
shg
plot(W,x(:,5),'LineWidth',2);
title('N2O moles VS catalyst weight');
xlabel('W (g)');
ylabel('N2O (mol)');

shg
plot(W,x(:,6),'LineWidth',2);
title('O2 moles VS catalyst weight');
xlabel('W (g)');
ylabel('O2 (mol)');
shg
plot(W,x(:,7),'LineWidth',2);
title('y VS catalyst weight');
xlabel('W (g)');
ylabel('y ');
shg
plot(W,x(:,8),'LineWidth',2);
title('Temp VS catalyst weight');
xlabel('W (g)');
ylabel('T (K)');
shg
%Steady state values to use in Objective 1
% nNOss = x(end,1) %mol NO
% nH2ss = x(end,2) %mol H2
% nN2ss = x(end,3)  % mol N2
% nH2Oss = x(end, 4) %mol H2O
% nTss = x(end,5) %K
% nPss = x(end,6) %Pa

function dxdW = PBR(~,x)
dxdW= zeros(8,1);
%Rewrite variables into vector form
nNO = x(1);
nH2 = x(2);
nN2 = x(3);
nH2O = x(4);
nN2O = x(5);
nO2 = x(6);
y = x(7);
T = x(8);
%Parameters
Ta = 20 +273;          % K, Heating jacket temperature
U = 110;            % W/m^2.K, Overall heat transfer coefficient
a = 80;              %m^2/g, heat transfer area per mass catalyst: Calculate
Nt = 4;             % -, Number of tubes
L = 8;                 % m, Tube length
D = 0.0266;           % m, Tube diameter
A = Nt*pi()*D*L;            %Heat transfer area
Ac = 5;                 %m2, cross-sectional area
Q = U*A*(Ta-T);
ff = 0.1;             % -, Tube friction factor: Get value
R = 8.314;            % Pa*m^3/mol.K, Universal gas constant
Dp = 1.25*10^-3;        %m, particle diameter
MNO = 0.03001;        % kg/mol, Molar mass of NO
MH2 = 0.001007;         % kg/mol, Molar mass of H2
MN2 = 0.016;         % kg/mol, Molar mass of N2
MN2O = 0.044013;
MO2 = 0.015999;
MH2O = 0.0180158;         % kg/mol, Molar mass of H2O
%Heat of reactions using enthalpies
CpNO = 29.93;           % J/mol.K, Specific heat capacity, NIST at 400 K
CpH2 = 29.18;          % J/mol.K, Specific heat capacity of H2
CpN2 = 29.25;          % J/mol.K, Specific heat capacity of N2
CpH2O = 79.59;         % J/mol.K, Specific heat capacity of H2O
CpN2O = 42.69;
CpO2 = 30.1;
deltaCp1 = 2*CpH2O + CpN2 - 2*CpH2 - 2*CpNO;
deltaCp2 = CpN2O + CpH2O - 2*CpNO - CpH2;
deltaCp3 = CpH2O - CpH2 - 0.5*CpO2;
deltaHrx1_ref = -241800-91300; %J/mol
deltaHrx2_ref = 0.5*-241800 +0.5*81600 -91300; %J/mol
deltaHrx3_ref = 2*-241800; %J/mol
deltaHrx1 = deltaHrx1_ref + deltaCp1*(T-298);        % J/mol.K, Heat of reaction for reaction 1: FUnctions of cp and T
deltaHrx2 = deltaHrx2_ref + deltaCp2*(T-298);        % J/mol.K, Heat of reaction for reaction 2
deltaHrx3 = deltaHrx3_ref + deltaCp3*(T-298);        % J/mol.K, Heat of reaction for reaction 3
%Initial mole compositions
yNO_0 = 0.7;            % Initial mole composition of NO
yO2_0 = 0.012;
yH2_0 = 1-yNO_0 - yO2_0; % Initial mole composition of H2
yN2_0 = 0;              % Initial mole composition of N2
yH2O_0 = 0;             % Initial mole composition of H2O
yN2O_0 = 0;
% External variables at inlet
T0 = 146+273;          % K, Inlet temperature
P0 = 162000;           % Pa, Inlet pressure
v0 = 150/(10^6*60);
CT_0 = P0/(R*T0);
P =y*P0;
%Relations
nT = nNO +nH2 + nN2 + nH2O + nN2O + nO2;
nT0 = 800000;
v = v0*(nT/nT0)*P0/P*T/T0;
yNO = nNO/nT;
yH2 = nH2/nT;
yN2 = nN2/nT;
yH2O = nH2O/nT;
yN2O = nN2O/nT;
yO2 = nO2/nT;
%Concentrations
cNO = CT_0*(nNO/nT)*(y)*(T0/T);
cH2 = CT_0*(nH2/nT)*(y)*(T0/T);
cN2 = CT_0*(nN2/nT)*(y)*(T0/T);
cH2O = CT_0*(nN2/nT)*(y)*(T0/T);
cO2 = CT_0*(nO2/nT)*(y)*(T0/T);
cN2O = CT_0*(nN2O/nT)*(y)*(T0/T);
%Partial pressures
pNO_0 = yNO_0*P0;
pH2_0 = yH2_0*P0;
pN2_0 = yN2_0*P0;
pH2O_0 = yH2O_0*P0;
pO2_0 = yO2_0*P0;
pN2O_0 = yN2O_0*P0;
pNO = cNO*R*T;
pH2 = cH2*R*T;
pN2 = cN2*R*T;
pH2O = cH2O*R*T;
pO2 = cO2*R*T;
pN2O = cN2O*R*T;
%Rate laws
p0 = 0.01; %((-4.95*10^(-4)*T) +0.244)*exp(4.59*pNO+(8.25*10^(-2)*pO2));
p0prime = 0.5*p0;
%Reaction rate constants (mol/g^-1s^-1)
kN2 = 8157*exp((-76.9*10^3)/(8.314*T));
kN2O = 0.716*exp((-45*10^3)/(8.314*T));
kH2O = 4.92*10^(11)*((-45*10^3)/(8.314*T));
%Equilibrium rate constants of adsorption (1/Pa)
KH2 = 2.12*10^(-4)*exp(58600/(8.314*T));
KNO = 2.83*10^(-6)*exp(77200/(8.314*T));
KO2 = 1.19*10^(-10)*exp(9650/(8.314*T));
%Langmuir-Hinshelwood model
RN2 = kN2*((KNO*pNO*KH2*(pH2-p0))/((1+KNO*pNO+KH2*(pH2-p0)+KO2*pO2)^2)) ;
RN2O = kN2O*((KNO*pNO*KH2*(pH2-p0))/((1+KNO*pNO+KH2*(pH2-p0)+KO2*pO2)^2)) ;
RH2O = kH2O*((KO2*pO2*KH2*(pH2-p0prime))/((1+KNO*pNO+KH2*(pH2-p0prime)+KO2*pO2)^2));
%Pressure drop (Help)
y = P/P0;
pc = 4000; %kg/m3, density of solid catalyst particles
phi = 0.7; %fraction of liquid
pb = pc*(1-phi); %mass of catalyst/volume of reactor bed, bulk density of catalyst
p_0 = yNO_0*pNO_0 + yH2_0*pH2_0 + yN2_0*pN2_0 + yH2O_0*pH2O_0 + yO2_0*pO2_0 + yN2O_0*pN2O_0;
u = v/Ac;
p = yNO*850 + yH2*600 + yN2*600 + yH2O*600 + yN2O*600 + nO2*600; %kg/m3, gas density: get as a function of temp
G = p*u;
gc = 32.174; %conversion factor, get correct units p. 178
visc = yNO*0.00001 + yH2*0.00001 + yN2*0.00001 + yH2O*0.00001 + yN2O*0.00001 + nO2*0.00001; %viscosity of gass passing through bed
B0 = (G*(1-phi))/(p_0*gc*Dp*phi^3)*(((150*(1-phi)*visc)/Dp)+ 1.75*G);
alpha = 2*B0/(Ac*pc*(1-phi)*P0);
%Mole balance R1
dnNO_dW = -2*RN2 -2*RN2O;
dnH2_dW = -2*RN2 -RN2O - RH2O;
dnN2_dW = RN2;
dnH2O_dW = 2*RN2 + RN2O + RH2O;
dnN2O_dW = RN2O;
dnO2_dW = -0.5*RH2O;
dy_dW = (-alpha./(2*y)).*(T/T0).*(nT/nT0);
dT_dW = ((U*a*(Ta-T))/pb + (-RN2)*(-deltaHrx1) + (-RN2O)*(-deltaHrx2) + (-RH2O)*(deltaHrx3))/(nNO*CpNO + nH2*CpH2 + nN2*CpN2 + nH2O*CpH2O + nN2O*CpN2O + nO2*CpO2);
dxdW = [dnNO_dW; dnH2_dW; dnN2_dW; dnH2O_dW; dnN2O_dW; dnO2_dW; dy_dW; dT_dW];
end

[Torsten]

Your code produces NaN for the vector dxdW (see above).

Check which operation in PBR is responsible for this behaviour.

[Thabo Ngwenya]

thank you very much................i have used the dxdW(isnan(dxdW(:,1)))= 100 operation by assigning 100to every Nan value , this is just a compromise for , apparently 7 out of 8 function operations are Nan and my guess is the first 7 are

[Benjamin Thompson]

Write a new question that is focused on your new issue, and attach your code and do not show so much MATLAB output. Include a description of what you are now trying to plot.

[Thabo Ngwenya]

my bad!.... thank you... i have rephrased my issue and put up a proper code

Answer 1

[Torsten]

In x0, set y to the value you want.

In PBR, set

dxdW = [dnNO_dW; dnH2_dW; dnN2_dW; dnH2O_dW; dnN2O_dW; dnO2_dW; 0.0; dT_dW];

[Thabo Ngwenya]

Thank very much ... i have done so and am getting no change. i will try optimising the y0 value , i must have chosen too big ,a big value