how to I solve 2 eqns with 2 unknowns and plot 5 plots according to these solutions?

1 view (last 30 days)
clc;
clear;
deltaGrxn_std_rxn1 = -24800; %(J/mol)
deltaHrxn_std_rxn1 = -90200; %(J/mol)
deltaSrxn_std_rxn1 = -219;
deltaGrxn_std_rxn2 = -28600;
deltaHrxn_std_rxn2 = -41200;
deltaSrxn_std_rxn2 = -42;
R = 8.314;
T = 298:10:800;
for j = 1:length(T)
deltaGrxn1_funcT = (deltaHrxn_std_rxn1) - (T(j) .* deltaSrxn_std_rxn1);
Ka_rxn1(j) = exp((-deltaGrxn1_funcT)./(R.*T(j)));
deltaGrxn2_funcT = (deltaHrxn_std_rxn2) - (T(j) .* deltaSrxn_std_rxn2);
Ka_rxn2(j) = exp((-deltaGrxn2_funcT)./(R.*T(j)));
syms ksi1
syms ksi2
y_CO = (1-ksi1-ksi2)./(3-2.*ksi1);
y_H2 = (2-2.*ksi1+ksi2)./(3-2.*ksi1);
y_CH3OH = (ksi1)./(3-2.*ksi1);
y_H2O = (1-ksi2)./(3-2.*ksi1);
y_CO2 = (ksi2)./(3-2.*ksi1);
eqn1 = ((y_CH3OH)./((y_H2O.^2).*(y_CO))) - Ka_rxn1(j) == 0
eqn2 = ((y_CO2.*y_H2)./(y_CO.*y_H2O)) - Ka_rxn2(j) ==0
E = [eqn1, eqn2];
S = vpasolve(E, ksi1, ksi2, [0 1]);
end
plot(T,y_CO)
hold on
plot(T,y_H2)
hold on
plot(T,y_CH3OH)
hold on
plot(T,y_H2O)
hold on
plot(T,y_CO2)

Answers (1)

Alan Stevens
Alan Stevens on 20 Jun 2021
Here's a possible way. Only you will know if the results make sense
deltaGrxn_std_rxn1 = -24800; %(J/mol)
deltaHrxn_std_rxn1 = -90200; %(J/mol)
deltaSrxn_std_rxn1 = -219;
deltaGrxn_std_rxn2 = -28600;
deltaHrxn_std_rxn2 = -41200;
deltaSrxn_std_rxn2 = -42;
R = 8.314;
T = 298:10:800;
deltaGrxn1_funcT = (deltaHrxn_std_rxn1) - (T.* deltaSrxn_std_rxn1);
Ka_rxn1 = exp((-deltaGrxn1_funcT)./(R.*T));
deltaGrxn2_funcT = (deltaHrxn_std_rxn2) - (T.* deltaSrxn_std_rxn2);
Ka_rxn2 = exp((-deltaGrxn2_funcT)./(R.*T));
y_CO = @(ksi1,ksi2) (1-ksi1-ksi2)./(3-2.*ksi1);
y_H2 = @(ksi1,ksi2) (2-2.*ksi1+ksi2)./(3-2.*ksi1);
y_CH3OH = @(ksi1,ksi2) (ksi1)./(3-2.*ksi1);
y_H2O = @(ksi1,ksi2) (1-ksi2)./(3-2.*ksi1);
y_CO2 = @(ksi1,ksi2) ksi2./(3-2.*ksi1);
yco = zeros(size(T));
yh2 = zeros(size(T));
ych3oh = zeros(size(T));
yh2o = zeros(size(T));
yco2 = zeros(size(T));
K = [0 1]; % Initial guesses
for j = 1:length(T)
K0 = K; % Update initial guesses
F = @(K) norm(((y_CH3OH(K(1),K(2)))./((y_H2O(K(1),K(2)).^2).*(y_CO(K(1),K(2))))) - Ka_rxn1(j))...
+ norm(((y_CO2(K(1),K(2)).*y_H2(K(1),K(2)))./(y_CO(K(1),K(2)).*y_H2O(K(1),K(2)))) - Ka_rxn2(j));
K = fminsearch(F,K0);
yco(j) = y_CO(K(1),K(2));
yh2(j) = y_H2(K(1),K(2));
ych3oh(j) = y_CH3OH(K(1),K(2));
yh2o(j) = y_H2O(K(1),K(2));
yco2(j) = y_CO2(K(1),K(2));
end
figure(1)
plot(T,yh2,T,yco2),grid
xlabel('T'), ylabel('y')
legend('H2','CO2')
figure(2)
plot(T,yco,T,yh2o),grid
xlabel('T'),ylabel('y')
legend('CO','H2O')
figure(3)
plot(T,ych3oh),grid
xlabel('T'),ylabel('y_CH3OH')
CO and H2O are numerically too close to separate graphically.

Categories

Find more on General Applications in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!