Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
chemical equilibrium of a gasifier using fsolve

Subject: chemical equilibrium of a gasifier using fsolve

From: Gavin

Date: 15 Jun, 2010 11:08:05

Message: 1 of 11

I'm trying to model a biomass gasifier using chemical equilibrium. I have a system of equations and unknowns that I would like to solve with Matlab using the "fsolve" command. I have two m-files which I use. The first m-file is named "gasf_funcB.m" where I create my system of equations to be used by fsolve. The second m-file is named "gasf_solB.m" where I state my inputs (total moles of C, H, O, N into the system) and initial guess for the outputs (mole fractions) obtained by fsolve. I am trying to solve the system for the unknowns which are the compounds created from the reaction. I am getting complex answers from using fsolve. Please help and let me know what I'm doing wrong. All m-files are shown below. My output to the workspace from running gasf_solB.m is also shown below.

---------------------------------------- Here is my gasf_funcB.m file --------------------------------------
% gasifier functions

function F = gasf_funcB(x,nC,nH,nO,nN,T)

% Overall reaction is
% CHON + air -> CO + CO2 + H2O + H2 + H + N + N2 + NO + O2 + O + OH
yCO = x(1);
yCO2 = x(2);
yH2O = x(3);
yH2 = x(4);
yH = x(5);
yN = x(6);
yN2 = x(7);
yNO = x(8);
yO2 = x(9);
yO = x(10);
yOH = x(11);

% total number of moles
nt = yCO+yCO2+yH2O+yH2+yH+yN+yN2+yNO+yO2+yO+yOH;

% determine equilibrium constant at specified temperature
[K] = K_est(T);
K1 = K(1);
K2 = K(2);
K3 = K(3);
K4 = K(4);
K5 = K(5);
K6 = K(6);
K7 = K(7);
K8 = K(8);

% balance equations to solve
F(1) = yCO+yCO2-nC/nt; % C balance
F(2) = 2*yH2O+2*yH2+yH+yOH-nH/nt; % H balance
F(3) = yCO+2*yCO2+yH2O+2*yO2+yO+yNO+yOH-nO/nt; % O balance
F(4) = 2*yN2+yN+yNO-nN/nt; % N balance

% equilibrium reactions to solve where 1) refers to K1 etc...
% assume P = 1 atm therefore P/Pref = 1
% 1) CO2 <-> CO + 0.5O2
% 2) CO2 + H2 <-> CO + H2O
% 3) N2 <-> 2N
% 4) H2 <-> 2H
% 5) O2 <-> 2O
% 6) 0.5O2 + 0.5N2 <-> NO
% 7) H2O <-> H2 + 0.5O2
% 8) H2O <-> OH + 0.5H2
F(5) = (yCO*yO2^0.5)/yCO2*(1/nt)^0.5-K1;
F(6) = (yCO*yH2O)/(yCO2*yH2)-K2;
F(7) = yN^2/yN2*(1/nt)-K3;
F(8) = yH^2/yH2*(1/nt)-K4;
F(9) = yO^2/yO2*(1/nt)-K5;
F(10) = yNO/(yO2^0.5*yN2^0.5)-K6;
F(11) = (yO2^0.5*yH2)/yH2O*(1/nt)^0.5-K7;
F(12) = (yOH*yH2^0.5)/yH2O*(1/nt)^0.5-K8;

end

% estimate of equilibrium constant (K) at temperature T in Kelvin for
% various reactions

function [K] = K_est(T)

% database of temperatures in Kelvin
Tdb = [298 500 1000 1200 1400 1600 1700 1800 1900 2000 2100 2200 2300 ...
    2400 2500 2600 2700 2800 2900 3000];

% database values for log base 10 equilibrium
% these correspond to equilibrium reactions stated previously
log10K1 = [-45.066 -25.025 -10.221 -7.764 -6.014 -4.706 -4.169 -3.693...
          -3.267 -2.884 -2.539 -2.226 -1.940 -1.679 -1.440 -1.219 -1.015...
          -0.825 -0.649 -0.485];
log10K2 = [-5.018 -2.139 -0.159 0.135 0.333 0.474 0.530 0.577 0.619...
          0.656 0.688 0.716 0.742 0.764 0.784 0.802 0.818 0.833 0.846...
          0.858];
log10K3 = [-159.6 -92.672 -43.056 -34.754 -28.812 -24.350 -22.512...
          -20.874 -19.410 -18.092 -16.898 -15.810 -14.818 -13.908...
          -13.070 -12.298 -11.580 -10.914 -10.294 -9.716];
log10K4 = [-71.224 -40.316 -17.292 -13.414 -10.63 -8.532 -7.666 -6.896...
          -6.204 -5.58 -5.016 -4.502 -4.032 -3.6 -3.202 -2.836 -2.494...
          -2.178 -1.882 -1.606];
log10K5 = [-81.208 -45.88 -19.614 -15.208 -12.054 -9.684 -8.706 -7.836...
          -7.058 -6.356 -5.72 -5.142 -4.614 -4.13 -3.684 -3.272 -2.892...
          -2.536 -2.206 -1.898];
log10K6 = [-15.171 -8.783 -4.062 -3.275 -2.712 -2.29 -2.116 -1.962...
           -1.823 -1.699 -1.586 -1.484 -1.391 -1.305 -1.227 -1.154...
           -1.087 -1.025 -0.967 -0.913];
log10K7 = [-40.048 -22.886 -10.062 -7.899 -6.347 -5.18 -4.699 -4.27...
           -3.886 -3.54 -3.227 -2.942 -2.682 -2.443 -2.224 -2.021...
           -1.833 -1.658 -1.495 -1.343];
log10K8 = [-46.054 -26.130 -11.280 -8.811 -7.021 -5.677 -5.124 -4.613...
           -4.19 -3.776 -3.434 -3.091 -2.809 -2.52 -2.27 -2.038 -1.823...
           -1.624 -1.438 -1.265];
       
% interpolate equilibrium constant at T
K1_interp = interp1(Tdb,log10K1,T,'spline','extrap');
K1 = 10^K1_interp;
K2_interp = interp1(Tdb,log10K2,T,'spline','extrap');
K2 = 10^K2_interp;
K3_interp = interp1(Tdb,log10K3,T,'spline','extrap');
K3 = 10^K3_interp;
K4_interp = interp1(Tdb,log10K4,T,'spline','extrap');
K4 = 10^K4_interp;
K5_interp = interp1(Tdb,log10K5,T,'spline','extrap');
K5 = 10^K5_interp;
K6_interp = interp1(Tdb,log10K6,T,'spline','extrap');
K6 = 10^K6_interp;
K7_interp = interp1(Tdb,log10K7,T,'spline','extrap');
K7 = 10^K7_interp;
K8_interp = interp1(Tdb,log10K8,T,'spline','extrap');
K8 = 10^K8_interp;

K = [K1 K2 K3 K4 K5 K6 K7 K8];

end

----------------------------------------------------------------------------------------------------------------------

--------------------------------------------- Here is my gasf_solB.m file -----------------------------------
% gasifier equilibrium solution
clear
clc

% Inputs to system, total moles of C, H, O, N coming from CHON and air
nC = 2; % C total moles
nH = 2; % H total moles
nO = 4; % O total moles
nN = 1; % N total moles
T = 1000; % equilibrium temperature, K

% initial guess of outputs
x0 = [.1 .1 .1 .1 .1 .1 .1 .1 .1 .1 .1];

options = optimset('Display','iter','MaxIter',10000,'MaxFunEvals',10000);
[x,exitflag] = fsolve(@(x) gasf_funcB(x,nC,nH,nO,nN,T),x0,options);

yCO = x(1)
yCO2 = x(2)
yH2O = x(3)
yH2 = x(4)
yH = x(5)
yN = x(6)
yN2 = x(7)
yNO = x(8)
yO2 = x(9)
yO = x(10)
yOH = x(11)
nt = yCO+yCO2+yH2O+yH2+yH+yN+yN2+yNO+yO2+yO+yOH
T

----------------------------------------------------------------------------------------------------------------------

-------------------------------- Output from running gasf_solB.m file ---------------------------------
Warning: Trust-region-dogleg algorithm of FSOLVE cannot handle non-square systems; using
Levenberg-Marquardt algorithm instead.
> In fsolve at 324
  In gasf_solB at 16

                                        First-Order Norm of
 Iteration Func-count Residual optimality Lambda step
     0 12 13.2407 25.4 0.01
     1 24 0.924578 4.48 0.001 0.478901
     2 36 0.0103891 0.336 0.0001 0.288275
     3 48 0.00146625 0.142 1e-05 0.103138
     4 65 0.00115348 0.146 1 0.0106871
     5 77 0.000951684 0.143 0.1 0.00882184
     6 90 0.000838769 0.148 1 0.00729546
     7 102 0.000739851 0.145 0.1 0.00620988
     8 115 0.00067335 0.148 1 0.00564338
     9 127 0.000609141 0.146 0.1 0.00504862
    10 140 0.000564525 0.148 1 0.00464933
   143 1806 3.28161e-06 1.65 10000 4.43205e-06
   144 1819 2.42909e-06 1.18 100000 1.00569e-06
   145 1831 2.26005e-06 0.0609 10000 1.69937e-07
Optimization terminated: the magnitude of the search direction is less than options.TolX
yCO =
      0.32134 - 1.4945e-06i
yCO2 =
      0.74179 + 1.223e-06i
yH2O =
      0.32147 + 1.2992e-06i
yH2 =
      0.20083 - 5.1165e-07i
yH =
     0.018548 + 2.2209e-07i
yN =
     0.022879 + 1.2783e-07i
yN2 =
      0.25434 - 8.5025e-07i
yNO =
  -5.2241e-09 + 3.6983e-08i
yO2 =
  -5.1624e-07 + 1.7704e-08i
yO =
  -1.6125e-05 - 4.0885e-07i
yOH =
  -2.9233e-06 - 1.6344e-07i
nt =
       1.8812 - 5.0192e-07i
T =
        1000
----------------------------------------------------------------------------------------------------------------------

Subject: chemical equilibrium of a gasifier using fsolve

From: Torsten Hennig

Date: 15 Jun, 2010 11:33:41

Message: 2 of 11

> I'm trying to model a biomass gasifier using chemical
> equilibrium. I have a system of equations and
> unknowns that I would like to solve with Matlab using
> the "fsolve" command. I have two m-files which I
> use. The first m-file is named "gasf_funcB.m" where
> I create my system of equations to be used by fsolve.
> The second m-file is named "gasf_solB.m" where I
> I state my inputs (total moles of C, H, O, N into the
> system) and initial guess for the outputs (mole
> fractions) obtained by fsolve. I am trying to solve
> the system for the unknowns which are the compounds
> created from the reaction. I am getting complex
> answers from using fsolve. Please help and let me
> know what I'm doing wrong. All m-files are shown
> below. My output to the workspace from running
> gasf_solB.m is also shown below.
>
> ---------------------------------------- Here is my
> gasf_funcB.m file
> --------------------------------------
> % gasifier functions
>
> function F = gasf_funcB(x,nC,nH,nO,nN,T)
>
> % Overall reaction is
> % CHON + air -> CO + CO2 + H2O + H2 + H + N + N2 + NO
> + O2 + O + OH
> yCO = x(1);
> yCO2 = x(2);
> yH2O = x(3);
> yH2 = x(4);
> yH = x(5);
> yN = x(6);
> yN2 = x(7);
> yNO = x(8);
> yO2 = x(9);
> yO = x(10);
> yOH = x(11);
>
> % total number of moles
> nt = yCO+yCO2+yH2O+yH2+yH+yN+yN2+yNO+yO2+yO+yOH;
>
> % determine equilibrium constant at specified
> temperature
> [K] = K_est(T);
> K1 = K(1);
> K2 = K(2);
> K3 = K(3);
> K4 = K(4);
> K5 = K(5);
> K6 = K(6);
> K7 = K(7);
> K8 = K(8);
>
> % balance equations to solve
> F(1) = yCO+yCO2-nC/nt;
> % C
> % C
> % C balance
> F(2) = 2*yH2O+2*yH2+yH+yOH-nH/nt;
> % H balance
> F(3) = yCO+2*yCO2+yH2O+2*yO2+yO+yNO+yOH-nO/nt; % O
> balance
> F(4) = 2*yN2+yN+yNO-nN/nt;
> % N balance
>
> % equilibrium reactions to solve where 1) refers to
> K1 etc...
> % assume P = 1 atm therefore P/Pref = 1
> % 1) CO2 <-> CO + 0.5O2
> % 2) CO2 + H2 <-> CO + H2O
> % 3) N2 <-> 2N
> % 4) H2 <-> 2H
> % 5) O2 <-> 2O
> % 6) 0.5O2 + 0.5N2 <-> NO
> % 7) H2O <-> H2 + 0.5O2
> % 8) H2O <-> OH + 0.5H2
> F(5) = (yCO*yO2^0.5)/yCO2*(1/nt)^0.5-K1;
> F(6) = (yCO*yH2O)/(yCO2*yH2)-K2;
> F(7) = yN^2/yN2*(1/nt)-K3;
> F(8) = yH^2/yH2*(1/nt)-K4;
> F(9) = yO^2/yO2*(1/nt)-K5;
> F(10) = yNO/(yO2^0.5*yN2^0.5)-K6;
> F(11) = (yO2^0.5*yH2)/yH2O*(1/nt)^0.5-K7;
> F(12) = (yOH*yH2^0.5)/yH2O*(1/nt)^0.5-K8;
>
> end
>
> % estimate of equilibrium constant (K) at temperature
> T in Kelvin for
> % various reactions
>
> function [K] = K_est(T)
>
> % database of temperatures in Kelvin
> Tdb = [298 500 1000 1200 1400 1600 1700 1800 1900
> 2000 2100 2200 2300 ...
> 2400 2500 2600 2700 2800 2900 3000];
>
> % database values for log base 10 equilibrium
> % these correspond to equilibrium reactions stated
> previously
> log10K1 = [-45.066 -25.025 -10.221 -7.764 -6.014
> -4.706 -4.169 -3.693...
> -3.267 -2.884 -2.539 -2.226 -1.940 -1.679
> 940 -1.679 -1.440 -1.219 -1.015...
> -0.825 -0.649 -0.485];
> log10K2 = [-5.018 -2.139 -0.159 0.135 0.333 0.474
> 0.530 0.577 0.619...
> 0.656 0.688 0.716 0.742 0.764 0.784 0.802
> .784 0.802 0.818 0.833 0.846...
> 0.858];
> log10K3 = [-159.6 -92.672 -43.056 -34.754 -28.812
> -24.350 -22.512...
> -20.874 -19.410 -18.092 -16.898 -15.810
> 98 -15.810 -14.818 -13.908...
> -13.070 -12.298 -11.580 -10.914 -10.294
> 14 -10.294 -9.716];
> log10K4 = [-71.224 -40.316 -17.292 -13.414 -10.63
> -8.532 -7.666 -6.896...
> -6.204 -5.58 -5.016 -4.502 -4.032 -3.6
> 4.032 -3.6 -3.202 -2.836 -2.494...
> -2.178 -1.882 -1.606];
> log10K5 = [-81.208 -45.88 -19.614 -15.208 -12.054
> -9.684 -8.706 -7.836...
> -7.058 -6.356 -5.72 -5.142 -4.614 -4.13
> .614 -4.13 -3.684 -3.272 -2.892...
> -2.536 -2.206 -1.898];
> log10K6 = [-15.171 -8.783 -4.062 -3.275 -2.712 -2.29
> -2.116 -1.962...
> -1.823 -1.699 -1.586 -1.484 -1.391 -1.305
> .391 -1.305 -1.227 -1.154...
> -1.087 -1.025 -0.967 -0.913];
> log10K7 = [-40.048 -22.886 -10.062 -7.899 -6.347
> -5.18 -4.699 -4.27...
> -3.886 -3.54 -3.227 -2.942 -2.682 -2.443
> .682 -2.443 -2.224 -2.021...
> -1.833 -1.658 -1.495 -1.343];
> log10K8 = [-46.054 -26.130 -11.280 -8.811 -7.021
> -5.677 -5.124 -4.613...
> -4.19 -3.776 -3.434 -3.091 -2.809 -2.52
> 2.809 -2.52 -2.27 -2.038 -1.823...
> -1.624 -1.438 -1.265];
>
> % interpolate equilibrium constant at T
> K1_interp = interp1(Tdb,log10K1,T,'spline','extrap');
> K1 = 10^K1_interp;
> K2_interp = interp1(Tdb,log10K2,T,'spline','extrap');
> K2 = 10^K2_interp;
> K3_interp = interp1(Tdb,log10K3,T,'spline','extrap');
> K3 = 10^K3_interp;
> K4_interp = interp1(Tdb,log10K4,T,'spline','extrap');
> K4 = 10^K4_interp;
> K5_interp = interp1(Tdb,log10K5,T,'spline','extrap');
> K5 = 10^K5_interp;
> K6_interp = interp1(Tdb,log10K6,T,'spline','extrap');
> K6 = 10^K6_interp;
> K7_interp = interp1(Tdb,log10K7,T,'spline','extrap');
> K7 = 10^K7_interp;
> K8_interp = interp1(Tdb,log10K8,T,'spline','extrap');
> K8 = 10^K8_interp;
>
> K = [K1 K2 K3 K4 K5 K6 K7 K8];
>
> end
>
> ------------------------------------------------------
> ------------------------------------------------------
> ----------
>
> --------------------------------------------- Here is
> my gasf_solB.m file
> -----------------------------------
> % gasifier equilibrium solution
> clear
> clc
>
> % Inputs to system, total moles of C, H, O, N coming
> from CHON and air
> nC = 2; % C total moles
> nH = 2; % H total moles
> nO = 4; % O total moles
> nN = 1; % N total moles
> T = 1000; % equilibrium temperature, K
>
> % initial guess of outputs
> x0 = [.1 .1 .1 .1 .1 .1 .1 .1 .1 .1 .1];
>
> options =
> optimset('Display','iter','MaxIter',10000,'MaxFunEvals
> ',10000);
> [x,exitflag] = fsolve(@(x)
> gasf_funcB(x,nC,nH,nO,nN,T),x0,options);
>
> yCO = x(1)
> yCO2 = x(2)
> yH2O = x(3)
> yH2 = x(4)
> yH = x(5)
> yN = x(6)
> yN2 = x(7)
> yNO = x(8)
> yO2 = x(9)
> yO = x(10)
> yOH = x(11)
> nt = yCO+yCO2+yH2O+yH2+yH+yN+yN2+yNO+yO2+yO+yOH
> T
>
> ------------------------------------------------------
> ------------------------------------------------------
> ----------
>
> -------------------------------- Output from running
> gasf_solB.m file ---------------------------------
> Warning: Trust-region-dogleg algorithm of FSOLVE
> cannot handle non-square systems; using
> Levenberg-Marquardt algorithm instead.
> > In fsolve at 324
> In gasf_solB at 16
>
> First-Order
> First-Order
> First-Order Norm
> Norm of
> Iteration Func-count Residual optimality
> y Lambda step
> 0 12 13.2407 25.4
> 25.4 0.01
> 1 24 0.924578 4.48
> 4.48 0.001 0.478901
> 2 36 0.0103891 0.336
> 0.336 0.0001 0.288275
> 3 48 0.00146625 0.142
> 0.142 1e-05 0.103138
> 4 65 0.00115348 0.146
> 0.146 1 0.0106871
> 5 77 0.000951684 0.143
> 0.143 0.1 0.00882184
> 6 90 0.000838769 0.148
> 0.148 1 0.00729546
> 7 102 0.000739851 0.145
> 0.145 0.1 0.00620988
> 8 115 0.00067335 0.148
> 0.148 1 0.00564338
> 9 127 0.000609141 0.146
> 0.146 0.1 0.00504862
> 10 140 0.000564525 0.148
> .148 1 0.00464933
> 143 1806 3.28161e-06 1.65
> .65 10000 4.43205e-06
> 144 1819 2.42909e-06 1.18
> .18 100000 1.00569e-06
> 145 1831 2.26005e-06 0.0609
> 609 10000 1.69937e-07
> Optimization terminated: the magnitude of the search
> direction is less than options.TolX
> yCO =
> 0.32134 - 1.4945e-06i
> yCO2 =
> 0.74179 + 1.223e-06i
> yH2O =
> 0.32147 + 1.2992e-06i
> yH2 =
> 0.20083 - 5.1165e-07i
> yH =
> 0.018548 + 2.2209e-07i
> yN =
> 0.022879 + 1.2783e-07i
> yN2 =
> 0.25434 - 8.5025e-07i
> yNO =
> -5.2241e-09 + 3.6983e-08i
> yO2 =
> -5.1624e-07 + 1.7704e-08i
> yO =
> -1.6125e-05 - 4.0885e-07i
> yOH =
> -2.9233e-06 - 1.6344e-07i
> nt =
> 1.8812 - 5.0192e-07i
> T =
> 1000
> ------------------------------------------------------
> ------------------------------------------------------
> ----------

Hi Gavin,

during iteration with fsolve some of your mole fractions
seem to become negative - and taking the square root
gives complex values.
Whereever there is a square root of a mole fraction y,
you should substitute y := y~^2 which makes the
square root disappear.

Best wishes
Torsten.

Subject: chemical equilibrium of a gasifier using fsolve

From: Gavin

Date: 15 Jun, 2010 13:45:25

Message: 3 of 11

Torsten Hennig <Torsten.Hennig@umsicht.fhg.de> wrote in message <963770834.348103.1276601652274.JavaMail.root@gallium.mathforum.org>...
> >
> > ------------------------------------------------------
> > ----------
>
> Hi Gavin,
>
> during iteration with fsolve some of your mole fractions
> seem to become negative - and taking the square root
> gives complex values.
> Whereever there is a square root of a mole fraction y,
> you should substitute y := y~^2 which makes the
> square root disappear.
>
> Best wishes
> Torsten.

Torsten,

I have changed the following equations to the following to prevent complex numbers:

F(5) = (yCO*(yO2^2)^0.5)/yCO2*(1/nt)^0.5-K1;
F(6) = (yCO*yH2O)/(yCO2*yH2)-K2;
F(7) = yN^2/yN2*(1/nt)-K3;
F(8) = yH^2/yH2*(1/nt)-K4;
F(9) = yO^2/yO2*(1/nt)-K5;
F(10) = yNO/((yO2^2)^0.5*yN2)^0.5-K6;
F(11) = ((yO2^2)^0.5*yH2)/yH2O*(1/nt)^0.5-K7;
F(12) = (yOH*(yH2^2)^0.5)/yH2O*(1/nt)^0.5-K8;

The complex answers are gone but I'm getting negative values for some of my answers which is not correct:
yO2 =
  -2.2808e-10
yO =
  -9.1383e-11

Is there a way to tell fsolve that a term cannot be negative when evaluating a system of equations?

Subject: chemical equilibrium of a gasifier using fsolve

From: Torsten Hennig

Date: 15 Jun, 2010 14:17:21

Message: 4 of 11

> Torsten Hennig <Torsten.Hennig@umsicht.fhg.de> wrote
> in message
> <963770834.348103.1276601652274.JavaMail.root@gallium.
> mathforum.org>...
> > >
> > >
> ------------------------------------------------------
> > > ----------
> >
> > Hi Gavin,
> >
> > during iteration with fsolve some of your mole
> fractions
> > seem to become negative - and taking the square
> root
> > gives complex values.
> > Whereever there is a square root of a mole fraction
> y,
> > you should substitute y := y~^2 which makes the
> > square root disappear.
> >
> > Best wishes
> > Torsten.
>
> Torsten,
>
> I have changed the following equations to the
> following to prevent complex numbers:
>
> F(5) = (yCO*(yO2^2)^0.5)/yCO2*(1/nt)^0.5-K1;
> F(6) = (yCO*yH2O)/(yCO2*yH2)-K2;
> F(7) = yN^2/yN2*(1/nt)-K3;
> F(8) = yH^2/yH2*(1/nt)-K4;
> F(9) = yO^2/yO2*(1/nt)-K5;
> F(10) = yNO/((yO2^2)^0.5*yN2)^0.5-K6;
> F(11) = ((yO2^2)^0.5*yH2)/yH2O*(1/nt)^0.5-K7;
> F(12) = (yOH*(yH2^2)^0.5)/yH2O*(1/nt)^0.5-K8;
>
> The complex answers are gone but I'm getting negative
> values for some of my answers which is not correct:
> yO2 =
> -2.2808e-10
> yO =
> -9.1383e-11
>
> Is there a way to tell fsolve that a term cannot be
> negative when evaluating a system of equations?

Hi Gavin,

when returning from fsolve, did you retransform your
solution, i.e. squared the mole fraction for
which you inserted the square in the algebraic
equations ?

Best wishes
Torsten.

Subject: chemical equilibrium of a gasifier using fsolve

From: Gavin

Date: 15 Jun, 2010 15:45:19

Message: 5 of 11

Torsten Hennig <Torsten.Hennig@umsicht.fhg.de> wrote in message <800314764.348942.1276611471963.JavaMail.root@gallium.mathforum.org>...
> >
> > Is there a way to tell fsolve that a term cannot be
> > negative when evaluating a system of equations?
>
> Hi Gavin,
>
> when returning from fsolve, did you retransform your
> solution, i.e. squared the mole fraction for
> which you inserted the square in the algebraic
> equations ?
>
> Best wishes
> Torsten.

I'm not sure what you mean by "retransform your solution when returning to fsolve". Also, another way of dealing with the complex numbers was to make sure the mole fraction could not be less than 0. I added the following to the gasf_funcB.m file:

if yCO<=0; yCO=1e-20; end
if yCO2<=0; yCO2=1e-20; end
if yH2O<=0; yH2O=1e-20; end
if yH2<=0; yH2=1e-20; end
if yH<=0; yH=1e-20; end
if yN<=0; yN=1e-20; end
if yN2<=0; yN2=1e-20; end
if yNO<=0; yNO=1e-20; end
if yO2<=0; yO2=1e-20; end
if yO<=0; yO=1e-20; end
if yOH<=0; yOH=1e-20; end

Subject: chemical equilibrium of a gasifier using fsolve

From: Gavin

Date: 15 Jun, 2010 22:46:04

Message: 6 of 11

"Gavin " <wigging@gmail.com> wrote in message <hv876f$qgn$1@fred.mathworks.com>...
> Torsten Hennig <Torsten.Hennig@umsicht.fhg.de> wrote in message <800314764.348942.1276611471963.JavaMail.root@gallium.mathforum.org>...
>
> if yCO<=0; yCO=1e-20; end
> if yCO2<=0; yCO2=1e-20; end
> if yH2O<=0; yH2O=1e-20; end
> if yH2<=0; yH2=1e-20; end
> if yH<=0; yH=1e-20; end
> if yN<=0; yN=1e-20; end
> if yN2<=0; yN2=1e-20; end
> if yNO<=0; yNO=1e-20; end
> if yO2<=0; yO2=1e-20; end
> if yO<=0; yO=1e-20; end
> if yOH<=0; yOH=1e-20; end

Even when using the "if" statements I still get negative values at low temperature ranges. This is still incorrect.

Subject: chemical equilibrium of a gasifier using fsolve

From: Torsten Hennig

Date: 16 Jun, 2010 06:45:09

Message: 7 of 11

> "Gavin " <wigging@gmail.com> wrote in message
> <hv876f$qgn$1@fred.mathworks.com>...
> > Torsten Hennig <Torsten.Hennig@umsicht.fhg.de>
> wrote in message
> <800314764.348942.1276611471963.JavaMail.root@gallium.
> mathforum.org>...
> >
> > if yCO<=0; yCO=1e-20; end
> > if yCO2<=0; yCO2=1e-20; end
> > if yH2O<=0; yH2O=1e-20; end
> > if yH2<=0; yH2=1e-20; end
> > if yH<=0; yH=1e-20; end
> > if yN<=0; yN=1e-20; end
> > if yN2<=0; yN2=1e-20; end
> > if yNO<=0; yNO=1e-20; end
> > if yO2<=0; yO2=1e-20; end
> > if yO<=0; yO=1e-20; end
> > if yOH<=0; yOH=1e-20; end
>
> Even when using the "if" statements I still get
> negative values at low temperature ranges. This is
> still incorrect.

If by 'negative' you mean values for the mole fractions
in the order of -1e-6, you can interprete them as 0.
If you don't want to get 'negative' values, just
substitute _all_ mole fraction in _all_ your equations
by their squares.
When you then get the solution vector [y1,...,y_n]
from fsolve, your true solution is [(y1)^2,...,(y_n)^2].
That's what I meant by: don't forget to retransform
the variables of which you took the square in your
problem formulation.

Best wishes
Torsten.

Subject: chemical equilibrium of a gasifier using fsolve

From: Gavin

Date: 17 Jun, 2010 03:20:21

Message: 8 of 11

Torsten Hennig <Torsten.Hennig@umsicht.fhg.de> wrote in message <324982408.352424.1276670740011.JavaMail.root@gallium.mathforum.org>...
> >
> > Even when using the "if" statements I still get
> > negative values at low temperature ranges. This is
> > still incorrect.
>
> If by 'negative' you mean values for the mole fractions
> in the order of -1e-6, you can interprete them as 0.
> If you don't want to get 'negative' values, just
> substitute _all_ mole fraction in _all_ your equations
> by their squares.
> When you then get the solution vector [y1,...,y_n]
> from fsolve, your true solution is [(y1)^2,...,(y_n)^2].
> That's what I meant by: don't forget to retransform
> the variables of which you took the square in your
> problem formulation.
>
> Best wishes
> Torsten.

Torsten,

It seems safe to assume that the solution value is zero if it's negative. I have added this assumption to the solve file. I plan to make more changes to the system and will reply if I have any other questions or problems. Thank you for all the help you have given and please monitor this thread for updates on this.

Subject: chemical equilibrium of a gasifier using fsolve

From: Steve Amphlett

Date: 17 Jun, 2010 09:58:05

Message: 9 of 11

"Gavin " <wigging@gmail.com> wrote in message <hvc49l$6e0$1@fred.mathworks.com>...
> Torsten Hennig <Torsten.Hennig@umsicht.fhg.de> wrote in message <324982408.352424.1276670740011.JavaMail.root@gallium.mathforum.org>...
> > >
> > > Even when using the "if" statements I still get
> > > negative values at low temperature ranges. This is
> > > still incorrect.
> >
> > If by 'negative' you mean values for the mole fractions
> > in the order of -1e-6, you can interprete them as 0.
> > If you don't want to get 'negative' values, just
> > substitute _all_ mole fraction in _all_ your equations
> > by their squares.
> > When you then get the solution vector [y1,...,y_n]
> > from fsolve, your true solution is [(y1)^2,...,(y_n)^2].
> > That's what I meant by: don't forget to retransform
> > the variables of which you took the square in your
> > problem formulation.
> >
> > Best wishes
> > Torsten.
>
> Torsten,
>
> It seems safe to assume that the solution value is zero if it's negative. I have added this assumption to the solve file. I plan to make more changes to the system and will reply if I have any other questions or problems. Thank you for all the help you have given and please monitor this thread for updates on this.

This system of equations is very difficult to solve using traditional methods. When your concentrations approach zero, some of the equations are either no longer valid or just go crazy (division by near zero). I'd be suprised if simply packaging them up and passing them to a canned solver is going to work without a lot of fiddling. I would search the literature for custom solvers for this well-known problem and implement them directly in Matlab.

Subject: chemical equilibrium of a gasifier using fsolve

From: Gavin

Date: 17 Jun, 2010 13:22:06

Message: 10 of 11

"Steve Amphlett" <Firstname.Lastname@Where-I-Work.com> wrote in message <hvcrjd$oug$1@fred.mathworks.com>...
>
> This system of equations is very difficult to solve using traditional methods. When your concentrations approach zero, some of the equations are either no longer valid or just go crazy (division by near zero). I'd be suprised if simply packaging them up and passing them to a canned solver is going to work without a lot of fiddling. I would search the literature for custom solvers for this well-known problem and implement them directly in Matlab.

Steve,
Would you happen to know of any such solvers that I could apply to this problem? Thanks.

Subject: chemical equilibrium of a gasifier using fsolve

From: Steve Amphlett

Date: 17 Jun, 2010 14:47:05

Message: 11 of 11

"Gavin " <wigging@gmail.com> wrote in message <hvd7hu$202$1@fred.mathworks.com>...
> "Steve Amphlett" <Firstname.Lastname@Where-I-Work.com> wrote in message <hvcrjd$oug$1@fred.mathworks.com>...
> >
> > This system of equations is very difficult to solve using traditional methods. When your concentrations approach zero, some of the equations are either no longer valid or just go crazy (division by near zero). I'd be suprised if simply packaging them up and passing them to a canned solver is going to work without a lot of fiddling. I would search the literature for custom solvers for this well-known problem and implement them directly in Matlab.
>
> Steve,
> Would you happen to know of any such solvers that I could apply to this problem? Thanks.

I think most engine simulation (1D, 3D) programs have some form of solver for these equations, but none would be prepared to hand it out. When I worked on ours, I was heavily guided by Heywood's "Internal Combustion Engine Fundamentals", probably not quite in your field, but the same problems to solve nevertheless. He says: "Standardized computer methodsfor the calculationof complex chemical equilibrium have been developed. A generally available and well documented example is the NASA program of this type." He then goes on to spend a few pages describing how it works and shows some plots of outputs. The NASA ref is Svelhla, R. A. and McBride, B, J.

http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19730007227_1973007227.pdf

It's old (1973) and it's in FORTRAN IV.

Others exist, but they generally are for sale, not free. Googling "Thermochemical Equilibrium" provides some interesting links that have all the right words. E.g. (possible wrapping problems here): http://www.sciencedirect.com/science?_ob=ArticleURL&_udi=B6V2P-4KCHCY6-1&_user=10&_coverDate=01%2F31%2F2007&_rdoc=1&_fmt=high&_orig=search&_sort=d&_docanchor=&view=c&_searchStrId=1372914681&_rerunOrigin=google&_acct=C000050221&_version=1&_urlVersion=0&_userid=10&md5=7d7e09d97b2f5234f0a07d3ed958145b

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us