Can someone helps me to solve the error in this code?
1 view (last 30 days)
Show older comments
il = 4;
R = 8.314;
F = 96485;
Alpha = 0.5;
n = 2;
Eexc = 53990.065;
sigma_ref = 0.020;
Epro = 18912.42;
delta = 7*2.54*10^-3;
I = 0:50:350;
for i = 1:length(I)
t=0;
for T = [313 333 353]
t=t+1;
% Anode
ia = (10^-7)(exp(-(Eexc/R)((1/T)-(1/298))));
Eta_Anode(i,t) = ((R*T)/(Alpha*n*F))*(log(I(i)/ia));
% Cathode
ic = (10^-3)(exp(-(Eexc/R)((1/T)-(1/353))));
Eta_Cathode(i,t) = -((R*T)/(Alpha*n*F))*(log(I(i)/ic));
%Ohmic
sigma(i,T) = sigma_ref*(exp(-(Epro/R)*((1/T)-(1/298))));
Rohm = delta/sigma(i,T);
Eta_Ohmic(i,t) = I(i)*Rohm;
%Concentration
Eta_Conc(i,t) = ((R*T)/(Alpha*n*F))*(log(il/(il-I(i))));
%Voltage
H_h2(i,T)=0+(28.84*(T-298));
H_o2(i,T)=0+(28.91*(T-298));
H_h2o(i,T)=-285830+(75.37*(T-298));
H(i,T)=H_h2(i,T)+(0.5*H_o2(i,T))-H_h2o(i,T);
S_h2(i,T)=130.68+(28.84*(log(T/298)));
S_o2(i,T)=205+(28.91*(log(T/298)));
S_h2o(i,T)=69.95+(75.37*(log(T/298)));
S(i,T)=S_h2(i,T)+(0.5*S_o2(i,T))-S_h2o(i,T);
G(i,T)=H(i,T)-(T*S(i,T));
E(i,T)=(G(i,T)/(2*F));
ET(i,t) = 1.5184-(1.5421*10^-3*T)+(9.523*10^-5*T*log(T))+(9.84*10^-8*T^2);
EP(i,t)= 1.23+(((R*T)/(n*F))*(log(0.5*0.5^0.5)));
VT(i,t) = ET(i,t) + Eta_Anode(i,t) - Eta_Cathode(i,t) + Eta_Ohmic(i,t) + Eta_Conc(i,t);
VP(i,t) = EP(i,t) + Eta_Anode(i,t) - Eta_Cathode(i,t) + Eta_Ohmic(i,t) + Eta_Conc(i,t);
V(i,t) = E(i,t) + Eta_Anode(i,t) - Eta_Cathode(i,t) + Eta_Ohmic(i,t) + Eta_Conc(i,t);
%Mass Flow Rate
A(i) = I(i)* 50;
M_H2(i) = (A(i)*0.99)/(2*F);
M_O2(i) = (A(i)*0.99)/(4*F);
M_H2O(i) = (A(i)*0.99)/(2*F);
disp (['At Current = ',num2str(I(i)),'A',', Mass Flow Rate of H2 = ', num2str(M_H2(i)),' mol/s',', Mass Flow Rate of O2 = ', num2str(M_O2(i)),' mol/s',', Mass Flow Rate of H2O = ', num2str(M_H2O(i)),' mol/s']);
end
end
figure(1);
plot (I,V(:,1),I,V(:,2),I,V(:,3),I,V(:,4),I,V(:,5),I,V(:,6))
legend('T=25','T=40','T=56','T=64','T=55','T=48')
xlabel('Current Density (A/cm2)');
ylabel('Cell Voltage (Volt)');
title('I-V Curve (General Equation of E reversible)');
figure(2);
plot (I,VP(:,1),I,VP(:,2),I,VP(:,3),I,VP(:,4),I,VP(:,5),I,VP(:,6))
legend('T=25','T=40','T=56','T=64','T=55','T=48')
xlabel('Current Density (A/cm2)');
ylabel('Cell Voltage (Volt)');
title('I-V Curve (Pressure & Temperature variation)');
figure(3);
plot (I,VT(:,1),I,VT(:,2),I,VT(:,3),I,VT(:,4),I,VT(:,5),I,VT(:,6))
legend('T=25','T=40','T=56','T=64','T=55','T=48')
xlabel('Current Density (A/cm2)');
ylabel('Cell Voltage (Volt)');
title('I-V Curve (Temperature variation)');
figure(4);
plot(I,M_H2(:),I,M_H2(:),I,M_H2(:),I,M_H2(:),I,M_H2(:),I,M_H2(:))
xlabel('Current Density (A/cm2)');
ylabel('Molar Flow Rate (mol/s)');
title('Current vs Molar Flow Rate of H2');
figure(5);
plot(I,Eta_Anode(:,1),I,Eta_Anode(:,2),I,Eta_Anode(:,3),I,Eta_Anode(:,4),I,Eta_Anode(:,5),I,Eta_Anode(:,6),I,Eta_Cathode(:,1),I,Eta_Cathode(:,2),I,Eta_Cathode(:,3),I,Eta_Cathode(:,4),I,Eta_Cathode(:,5),I,Eta_Cathode(:,6),I,Eta_Ohmic(:,1),I,Eta_Ohmic(:,2),I,Eta_Ohmic(:,3),I,Eta_Ohmic(:,4),I,Eta_Ohmic(:,5),I,Eta_Ohmic(:,6),I,Eta_Conc(:,1),I,Eta_Conc(:,2),I,Eta_Conc(:,3),I,Eta_Conc(:,4),I,Eta_Conc(:,5),I,Eta_Conc(:,6))
legend('T=25','T=40','T=56','T=64','T=55','T=48')
xlabel('Current Density (A/cm2)');
ylabel('Overpotential (Volt)');
title('Different Overpotential');
figure(6);
plot(I,Eta_Anode(:,1),I,Eta_Anode(:,2),I,Eta_Anode(:,3),I,Eta_Anode(:,4),I,Eta_Anode(:,5),I,Eta_Anode(:,6))
xlabel('Current Density (A/cm2)');
ylabel('Overpotential (Volt)');
title('Activation Anode Overpotential');
figure(7); plot(I,Eta_Cathode(:,1),I,Eta_Cathode(:,2),I,Eta_Cathode(:,3),I,Eta_Cathode(:,4),I,Eta_Cathode(:,5),I,Eta_Cathode(:,6))
xlabel('Current Density (A/cm2)');
ylabel('Overpotential (Volt)');
title('Activation Cathode Overpotential');
figure(8); plot(I,Eta_Ohmic(:,1),I,Eta_Ohmic(:,2),I,Eta_Ohmic(:,3),I,Eta_Ohmic(:,4),I,Eta_Ohmic(:,5),I,Eta_Ohmic(:,6))
xlabel('Current Density (A/cm2)');
ylabel('Overpotential (Volt)');
title('Ohmic Overpotential');
figure(9); plot(I,Eta_Conc(:,1),I,Eta_Conc(:,2),I,Eta_Conc(:,3),I,Eta_Conc(:,4),I,Eta_Conc(:,5),I,Eta_Conc(:,6))
xlabel('Current Density (A/cm2)');
ylabel('Overpotential (Volt)');
title('Concentration Overpotential');
0 Comments
Answers (2)
Torsten
on 26 Jan 2024
Edited: Torsten
on 26 Jan 2024
You mean
ia = (10^-7)*(exp(-(Eexc/R)*((1/T)-(1/298))));
ic = (10^-3)*(exp(-(Eexc/R)*((1/T)-(1/353))));
instead of
ia = (10^-7)(exp(-(Eexc/R)((1/T)-(1/298))));
ic = (10^-3)(exp(-(Eexc/R)((1/T)-(1/353))));
?
And the second dimension of V, VP, VT, Eta_anode,... is only 3, not 6, so that you will get a dimension error in the plots.
0 Comments
Dyuman Joshi
on 26 Jan 2024
There are missing operators at many places in your code, e.g. -
% vv vv
ia = (10^-7)(exp(-(Eexc/R)((1/T)-(1/298))));
% vv vv
ic = (10^-3)(exp(-(Eexc/R)((1/T)-(1/353))));
Use appropriate operators and this issue will be resolved.
2 Comments
Dyuman Joshi
on 26 Jan 2024
Additionally, V is 8x3 where as you trying to access the 4th, 5th and 6th column of it, which don't exist.
Same for VT and VP.
I suggest you go through your code thoroughly once and check each line.
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!