0

I have 9 equations and 9 unknowns that I need to solve for, but when I set up the code and run it I get the following:

"Warning: Unable to find explicit solution. For options, see help.

In sym/solve (line 317) In HW2_PartA (line 65)

ans =

Empty sym: 0-by-1"

I have tried using both the solve functions and the vpasolve functions and I receive the same message each time.

The code is:

clear all
clc



syms T0 nO2 nH2 nH nOH nH2O nO lO lH
t = T0/1000;
R = 8.314;
nF = 1;
nOX = 11/32;
P0 = 5.3e6; % Pa
Pstp = 101.3e3;%Pa

hfH2 = 0; %j/mol
hfO2 = 0; %j/mol
hfH2O = -241.82e3;% j/mol
hfO = 249.19e3; %j/mol
hfH = 218e3;%J/mol
hfOH = 39.46e3;%j/mol

oxygengas = [20.9111 10.721071 -2.020498 0.046449 9.245722 5.337651 237.6185 0.00];
hydrogengas = [43.4136   -4.2931    1.2724   -0.0969  -20.5339  -38.5152  162.0814   0.00];
hydrogen = [20.7860    0.0000   -0.0000    0.0000    0.0000  211.8020  139.8711  217.9990];
OHgas = [28.7470    4.7145   -0.8147    0.0547   -2.7478   26.4144  214.1166   38.9871];
water = [41.9643    8.6221   -1.4998    0.0981  -11.1576 -272.1797  219.7809 -241.8264];

%Enthalpy
hO2 = 1000*(oxygengas*t + oxygengas(2)*(t^2)/2 + oxygengas(3)*(t^3)/3 + oxygengas(4)*(t^4)/4 - 
oxygengas(5)/t + oxygengas(6) - oxygengas(8));
hH2 = 1000*(hydrogengas(1)*t + hydrogengas(2)*(t^2)/2 + hydrogengas(3)*(t^3)/3 + hydrogengas(4)* 
(t^4)/4 - hydrogengas(5)/t + hydrogengas(6) - hydrogengas(8));
hH = 1000*(hydrogen(1)*t + hydrogen(2)*(t^2)/2 + hydrogen(3)*(t^3)/3 + hydrogen(4)*(t^4)/4 - 
hydrogen(5)/t + hydrogen(6) - hydrogen(8));
hOH = 1000*(OHgas(1)*t + OHgas(2)*(t^2)/2 + OHgas(3)*(t^3)/3 + OHgas(4)*(t^4)/4 - OHgas(5)/t + 
OHgas(6) - OHgas(8));
hO = hH;
hH2O = 1000*(water(1)*t + water(2)*(t^2)/2 + water(3)*(t^3)/3 + water(4)*(t^4)/4 - water(5)/t + 
water(6) - water(8));

%Entropy
sO2 = oxygengas(1)*log(t) + oxygengas(2)*t + oxygengas(3)*(t^2)/2 +oxygengas(4)*(t^3)/3 - 
oxygengas(5)/(2*t^2) +oxygengas(7);
sH2 = hydrogengas(1)*log(t) + hydrogengas(2)*t + hydrogengas(3)*(t^2)/2 +hydrogengas(4)*(t^3)/3 - 
hydrogengas(5)/(2*t^2) +hydrogengas(7);
sH = hydrogen(1)*log(t) + hydrogen(2)*t + hydrogen(3)*(t^2)/2 +hydrogen(4)*(t^3)/3 - 
hydrogen(5)/(2*t^2) +hydrogen(7);
sOH = OHgas(1)*log(t) + OHgas(2)*t + OHgas(3)*(t^2)/2 +OHgas(4)*(t^3)/3 - OHgas(5)/(2*t^2) +OHgas(7);
sO =  (5/2)*R*log(T0/298.15) + 161.003;
sH2O = water(1)*log(t) + water(2)*t + water(3)*(t^2)/2 +water(4)*(t^3)/3 - water(5)/(2*t^2) 
+water(7);

%Mass continuity
eqn1 = 2*nO2 + nOH + nO +nH2O - 2*nOX == 0;
eqn2 = 2*nH2 + 2*nH2O + nOH + nH - 2*nF == 0;

%Minimizing Gibbs Free Energy
%Oxygen Gas
eqn3 = hO2 + hfO2 - T0*sO2 + R*T0*log(P0/Pstp) + R*T0*log(nO2/(nO2+nH2+nH+nOH+nO+nH2O)) + 2*lO ==0 ;
%Hydrogen Gas
eqn4 = hH2 + hfH2 - T0*sH2 + R*T0*log(P0/Pstp) + R*T0*log(nH2/(nO2+nH2+nH+nOH+nO+nH2O)) + 2*lH == 0;
%Hydrogen
eqn5 = hH + hfH - T0*sH + R*T0*log(P0/Pstp) + R*T0*log(nH/(nO2+nH2+nH+nOH+nO+nH2O)) + lH == 0;
%OH Gas
eqn6 = hOH + hfOH - T0*sOH + R*T0*log(P0/Pstp) + R*T0*log(nOH/(nO2+nH2+nH+nOH+nO+nH2O)) + lO + lH == 
0;
%Oxygen
eqn7 = hO + hfO - T0*sO + R*T0*log(P0/Pstp) + R*T0*log(nO/(nO2+nH2+nH+nOH+nO+nH2O)) + lO == 0;
%Water
eqn8 = hH2O + hfH2O - T0*sO2 + R*T0*log(P0/Pstp) + R*T0*log(nH2O/(nO2+nH2+nH+nOH+nO+nH2O)) + 2*lH + 
lO == 0;

% Enthalpy Balance
eqn9 = nO2*hO2 + nH2*hH2 + nH*hH + nOH*hOH + nO*hO + nH2O*hH2O - nOX*hfO2 - nF*hfH2 + nO2*hfO2 + 
nH2*hfH2 + nH*hfH + nOH*hfOH + nO*hfOH + nH2O*hfH2O == 0;

eqns = [eqn1,eqn2,eqn3,eqn4,eqn5, eqn6, eqn7, eqn8, eqn9];
sol = solve([eqn1,eqn2,eqn3,eqn4,eqn5,eqn6,eqn7,eqn8,eqn9], [T0 nO2 nH2 nH nOH nH2O nO lO lH]);
sol.T0
  • If you want a numerical solution, why are you using the symbolic mechanisms of Matlab and not a numerical solver? – Lutz Lehmann Oct 24 '20 at 07:09
  • I have tried using the numerical solver vpasolve and have set it up as examples have showed using vpasolve but i get the following errors: Error using mupadengine/feval_internal More equations than variables is only supported for polynomial systems. Error in sym/vpasolve (line 172) sol = eng.feval_internal('symobj::vpasolve',eqns,vars,X0); Error in HW2_PartA (line 69) s = vpasolve([eqn1,eqn2,eqn3,eqn4,eqn5,eqn6,eqn7,eqn8,eqn9], [T0 nO2 nH2 nH nOH nH2O nO lO lH]) – Jacob Pearson Oct 24 '20 at 13:43
  • @JacobPearson: Shouldn't `hO2 = 1000*(oxygengas*t + ...` actually be `hO2 = 1000*(oxygengas(1)*t + ...`? – njuffa Nov 11 '20 at 03:35
  • @JacobPearson Could you review your equations? My attempts of solving as closely as possible (least-square error criterion) are not successful so far. Do we know for sure that a solution must exist? Is there anything in terms of domain knowledge that lets us put bounds on some or all variables? – njuffa Nov 11 '20 at 08:07
  • @JacobPearson Should `eqn8 = hH2O + hfH2O - T0*sO2 ...` actually be `eqn8 = hH2O + hfH2O - T0*sH2O ...` ? – njuffa Nov 12 '20 at 02:05

0 Answers0