MATLAB: Help with using fsolve

fefsolveimplicit equation

I can not figure out why matlab can not execute my function:

function fx = CSTR(x)
%CSTR
%CO + 0.5 O2 --> CO2
T = x(1);
phi = x(2);
nu_a = -1;
nu_b = -0.5;
nu_c = 1;
%intial reactor conditions
Po = 10^5; %Pa
To = 700; %K
Qo = 0.005; %m^3/s
Rg = 8.314; %Pa.m^3/mol.K
V = 0.3; %m^3
%moles
Ftoto = Po*Qo/(Rg*To); %mol/s






Fao = 0.02*Ftoto; %mol/s
Fbo = 0.5*Fao; %mol/s
Fco = 0; %mol/s
Fdo = 3.5*Fao; %mol/s
Feo = 45*Fao; %mol/s
Fa = Fao*phi; %mol/s
%change in flowrate
Q = ((0.5*phi+49.5)/50)*Qo*(T/To);
%rate equation
k =(1.26*10^10)*exp(-20131/T);
ra = -(0.5^0.25)*(3.5^0.5)*k*((Fao/Q)^1.75)*(phi.^1.25); 
%standard enthalpy of reaction at 298K
fH_a = -110.6*10^3; %J/mol




fH_b = 0; %J/mol
fH_c = -393.8*10^3;%J/mol
rH_298 = nu_a*fH_a + nu_b*fH_b + nu_c*fH_c;%J/mol
%heat capacity data
a_a = 28.11; b_a = 0.1672*10^(-2); c_a = 0.5363*10^(-5); d_a = -2.218*10^(-9);
a_b = 25.44; b_b = 1.518*10^(-2); c_b = -0.7144*10^(-5); d_b = 1.310*10^(-9);
a_c = 22.22; b_c = 5.9711*10^(-2); c_c = -3.495*10^(-5); d_c = 7.457*10^(-9);
a_d = 32.19; b_d = 0.1920*10^(-2); c_d = 1.054*10^(-5); d_d = -3.589*10^(-9);
a_e = 28.85; b_e = -0.1569*10^(-2); c_e = 0.8067*10^(-5); d_e = -2.868*10^(-9);
da = nu_a*a_a + nu_b*a_b + nu_c*a_c;
db = nu_a*b_a + nu_b*b_b + nu_c*b_c;
dc = nu_a*c_a + nu_b*c_b + nu_c*c_c;
dd = nu_a*d_a + nu_b*d_b + nu_c*d_c;
%Let g_i(T) = Int_To to T of Cpi)
g_a = a_a(T-To) + (1/2)*b_a*(T^2-To^2) + (1/3)*c_a*(T^3-To^3) + (1/4)*d_a*(T^3-To^3);
g_b = a_b(T-To) + (1/2)*b_b*(T^2-To^2) + (1/3)*c_b*(T^3-To^3) + (1/4)*d_b*(T^3-To^3);
g_c = a_c(T-To) + (1/2)*b_c*(T^2-To^2) + (1/3)*c_c*(T^3-To^3) + (1/4)*d_c*(T^3-To^3);
g_d = a_d(T-To) + (1/2)*b_d*(T^2-To^2) + (1/3)*c_d*(T^3-To^3) + (1/4)*d_d*(T^3-To^3);
g_e = a_e(T-To) + (1/2)*b_e*(T^2-To^2) + (1/3)*c_e*(T^3-To^3) + (1/4)*d_e*(T^3-To^3);
Sum_g_i= Fao*g_a + Fbo*g_b + Fco*g_c + Fdo*g_d + Feo*g_e; %J/s
 
%enthalpy of reaction
rH = rH_298 + da*(T-298) + db*(1/2)*(T^2-298^2) + dc*(1/3)*(T^3-298^3) + dd*(1/4)*(T^2-298^2); %J/mol
%energy balance and mol balance
fx(1) = Sum_g_i - ra*V*rH;
fx(2) = Fao - Fa + ra*V;
end

I am calling this function using the following:

xguess = [860 0.1];
y = fsolve(@CSTR,xguess);
T = y(1)
phi = y(2)

I get this error:

Error in A8 (line 3)

y =fsolve(@CSTR,xguess);

Caused by:

Failure in initial objective function evaluation. FSOLVE cannot continue.

I am unsire why the initial objective function can not be evaluted.

Would therer be a better way to evaulate these implicit and explicit functions?

thank you!

Best Answer

a_a = 28.11;

So a_a is a scalar.

g_a = a_a(T-To) + (1/2)*b_a*(T^2-To^2) + (1/3)*c_a*(T^3-To^3) + (1/4)*d_a*(T^3-To^3);

And there you are indexing the scalar at (T-To)

Similar problem on the lines that follow that.

Related Solutions

MATLAB: Error when using ode15

It is good practice to try and avoid globals, but rather pass them as auxilary variables into the function.

Fao = 1487.351; % kmol/h
Tr = 273; % K

z = 1;
a = (3 * 1.5) / ((3 * 1.5) + 1);
b1 = (0.5 + (0.5 / 1.5)) / 1;
b2 = (-0.5+(1.5*1.5))/1;
R = 8.314 / 1000; % kJ/mol*K
K30 = 0.0307;
E2 = 72100.82; % kJ/mol

E3 = -80989.98; % kJ/mol
omega = 1.564;
alpha = 0.640;
k20 = 2.12 * (10 ^ 13);  % should be 2.12e13 I suspect 
Bo = 0.254626; %atm/m
rho_c = 1600; %kg/m^3
Phi = 0.5;
D = 0.05; %Diameter of tube, m
Ac = pi * (D^2)/4; %m^2
X0 = 0;
T0 = 673; % K
P0 = 150; %atm
x0 = [X0 T0 P0];
Wspan = [0 5];
[W, x] = ode15s(@(t,x) fchem(t,x,Fao, Tr, z, a, b1, b2, R, K30, E2, E3, omega, alpha, k20, Bo, rho_c, Phi, Ac, T0, P0) ... 
           , Wspan, x0);
       
plot(W,x)

Now you function is (which I've renamed to make it more user-friendly)

 function func = fchem(t,x,Fao, Tr, z, a, b1, b2, R, K30, E2, E3, omega, alpha, k20, Bo, rho_c, Phi, Ac, T0, P0)
%% Identification of global variables
%global Fao Tr z a b1 b2 R K30 E2 E3 omega alpha k20 Bo rho_c Phi Ac T0 P0
%% Decleration of the non-linear equations
func = [((k20*(-E2/R*x(2))*(x(3)*(a/3)*(1-b2*z)*((10^(-5.519265*(10^(-5))*x(2)+...
    1.848863*(10^(-7))*(x(2)^2)+(2001.6/x(2))+2.6899)/(x(2)^(5.519265)))^2)-...
    (((x(3)*z)^2)/((x(3)*a*(1-b1*z))^3))))/((1+(K30*exp(-E3/R*x(2)))*(x(3)*z)/...
    ((x(3)*a*(1-b1*z))^omega))^(2*alpha)))/Fao; % dX/dW, x(1) = X
    
    ((2.12*(10^13)*(-E2/R*x(2))*(x(3)*(a/3)*(1-b2*z)*((10^(-5.519265*(10^(-5))*x(2)...
    + 1.848863*(10^(-7))*(x(2)^2)+(2001.6/x(2))+2.6899)/(x(2)^(5.519265)))^2)-...
    (((x(3)*z)^2)/((x(3)*a*(1-b1*z))^3))))/((1+(K30*exp(-E3/R*x(2)))*(x(3)*z)/...
    ((x(3)*a*(1-b1*z))^omega))^(2*alpha)))/(-(-91.63+(-12.96*(x(2)-Tr)+0.00917*((x(2)-Tr)^2)/2)))/...
    (Fao*(6.50+0.00100*x(2))+x(1)*2*(6.70+ 0.00630*x(2))-...
    3*(6.62+0.00081*x(2))-(6.50+0.00100*x(2))); % dT/dW, x(2) = T
    
    -(Bo/(Ac*(1-Phi)*rho_c))*(P0/x(3))*(x(2)/T0)*(1-0.5*x(1)) %dp/dW, x(3) = P
    ];
end

When I plot the solution, it looks pretty boring, but that's a different problem.

MATLAB: What could be the issue if a global variable gets an error “Undefined function or variable ‘F_FEED’ ” when called from a function. The global variable is defined in the file from which the function is called.

Maybe you mean

function main
global F_FEED T0 P0 WEIGHT_cat R Ua F_G0 F_H2O0 F_N20 F_A0 F_AC0 F_AA0 F_FA0 F_T0;
F_FEED = 100;           %TOTAL FLOWRATE OF THE FEED (mol/s)
T0 = 650;          %INLET TEMPERATURE OF THE REACTOR (K)
P0 = 303975;    %INLET PRESSURE OF THE REACTOR (Pa)
WEIGHT_cat = 100;       %CATALYST WEIGHT (kg)
R = ToBeSpecified;       
Ua = ToBeSpecified;
Feed = zeros(8,1);
Feed(1) = 0.014*F_FEED;           %INLET FLOWRATE OF G (mol/s)
Feed(2) = 0.634*F_FEED;           %INLET FLOWRATE OF H2O (mol/s)
Feed(3) = 0.352*F_FEED;           %INLET FLOWRATE OF N2 (mol/s)
Feed(8) = T0; 
F_G0 = Feed(1);                 
F_H2O0 = Feed(2);              
F_N20 = Feed(3);                
F_A0 = Feed(4);               
F_AC0 = Feed(5);                
F_AA0 = Feed(6);                
F_FA0= Feed(7);                
F_T0 = F_FEED;
[W,F] = ode45(@Reactors , [0 WEIGHT_cat], Feed);
end
function dFdW = Reactors(t,F)
global F_FEED T0 P0 WEIGHT_cat R Ua F_G0 F_H2O0 F_N20 F_A0 F_AC0 F_AA0 F_FA0 F_T0;
dFdW = zeros(8,1);
%dFdW = [r_G]
%       [r_H2O]
%       [r_N2]
%       [r_O2]
%       [r_A]
%       [r_AC]
%       [r_AA]
%       [r_FA]
%       [T]     Temperature
%reaction
Hf_G = -577.9;                    %Heat of formation: G (kJ/mol)
Hf_H2O = -241.818;                %Heat of formation: H2O (kJ/mol)
Hf_N2 = 0;                        %Heat of formation: N2 (kJ/mol)
Hf_A = -65;                       %Heat of formation: A (kJ/mol)
Hf_AC = -370.06;                  %Heat of formation: AC (kJ/mol)
Hf_AA = -166.1;                   %Heat of formation: AA (kJ/mol)
Hf_FA = -108.6;                   %Heat of formation: FA (kJ/mol)
Cp_G = -35.3267443 + 0.645106658*F(8) - 0.00061741*F(8)^2 + 0.000000239229*F(8)^3;               %provide integrated form
Cp_H2O = 32.95266198 - 0.00391977*F(8) + 0.0000222314*F(8)^2 - 0.0000000104953*F(8)^3;
Cp_N2 = 28.883 -0.00157*F(8)+ 0.00000808*F(8)^2 -0.000000002871*F(8)^3;
Cp_A = -14.9617445 + 0.373240705*F(8) - 0.00034559*F(8)^2 + 0.000000129748*F(8)^3;
Cp_AC = -11.7479865 + 0.402868987*F(8) - 0.00031782*F(8)^2 + 0.000000109142*F(8)^3;
Cp_AA = -5.82005966 + 0.23787576*F(8) - 0.00017362*F(8)^2 + 0.0000000543415*F(8)^3;
Cp_FA = 19.10519375 + 0.048680586*F(8) + 0.00000878405*F(8)^2 - 0.0000000150776*F(8)^3;
integraldeltaCp_rxn_1 = 2*(32.95266198*(F(8)-298) - 0.5*0.00391977*(F(8)^2-298^2) + (1/3)*0.0000222314*(F(8)^3-298^3) - (1/4)*0.0000000104953*(F(8)^4-298^4)) + (-14.9617445*(F(8)-298) + 0.5*0.373240705*(F(8)^2-298^2) - (1/3)*0.00034559*(F(8)^3-298^3) + (1/4)*0.000000129748*(F(8)^4-298^4)) - ((-35.3267443*(F(8)-298) + 0.5*0.645106658*(F(8)^2-298^2) - (1/3)*0.00061741*(F(8)^3-298^3) + (1/4)*0.000000239229*(F(8)^4-298^4))); 
integraldeltaCp_rxn_2 = (32.95266198*(F(8)-298) - 0.5*0.00391977*(F(8)^2-298^2) + (1/3)*0.0000222314*(F(8)^3-298^3) - (1/4)*0.0000000104953*(F(8)^4-298^4) + (-11.7479865*(F(8)-298) + 0.5*0.402868987*(F(8)^2-298^2) - (1/3)*0.00031782*(F(8)^3-298^3) + (1/4)*0.000000109142*(F(8)^4-298^4)) - (32.95266198*(F(8)-298) - 0.5*0.00391977*(F(8)^2-298^2) + (1/3)*0.0000222314*(F(8)^3-298^3) - (1/4)*0.0000000104953*(F(8)^4-298^4)));
integraldeltaCp_rxn_3 = (32.95266198*(F(8)-298) - 0.5*0.00391977*(F(8)^2-298^2) + (1/3)*0.0000222314*(F(8)^3-298^3) - (1/4)*0.0000000104953*(F(8)^4-298^4) + (19.10519375*(F(8)-298) + 0.5*0.048680586*(F(8)^2-298^2) + (1/3)*0.00000878405*(F(8)^3-298^3) - (1/4)*0.0000000150776*(F(8)^4-298^4)) + (-5.82005966*(F(8)-298) + 0.5*0.23787576*(F(8)^2-298^2) - (1/3)*0.00017362*(F(8)^3-298^3) + (1/4)*0.0000000543415*(F(8)^4-298^4)) - (-35.3267443*(F(8)-298) + 0.5*0.645106658*(F(8)^2-298^2) - (1/3)*0.00061741*(F(8)^3-298^3) + (1/4)*0.000000239229*(F(8)^4-298^4)));
integraldeltaCp_rxn_4 = (32.95266198*(F(8)-298) - 0.5*0.00391977*(F(8)^2-298^2) + (1/3)*0.0000222314*(F(8)^3-298^3) - (1/4)*0.0000000104953*(F(8)^4-298^4) + (19.10519375*(F(8)-298) + 0.5*0.048680586*(F(8)^2-298^2) + (1/3)*0.00000878405*(F(8)^3-298^3) - (1/4)*0.0000000150776*(F(8)^4-298^4)) + (-5.82005966*(F(8)-298) + 0.5*0.23787576*(F(8)^2-298^2) - (1/3)*0.00017362*(F(8)^3-298^3) + (1/4)*0.0000000543415*(F(8)^4-298^4)) - (-35.3267443*(F(8)-298) + 0.5*0.645106658*(F(8)^2-298^2) - (1/3)*0.00061741*(F(8)^3-298^3) + (1/4)*0.000000239229*(F(8)^4-298^4)));
Hrxn_1 =  2*Hf_H2O + Hf_A - Hf_G + (integraldeltaCp_rxn_1);               %Heat of reaction: reaction 1
Hrxn_2 =  2*Hf_H2O + Hf_AC - Hf_G + (integraldeltaCp_rxn_2);              %Heat of reaction: reaction 2
Hrxn_3 =  Hf_H2O + Hf_AA + Hf_FA - Hf_G + (integraldeltaCp_rxn_3);        %Heat of reaction: reaction 3
Hrxn_4 =  Hf_H2O + Hf_AA + Hf_FA - Hf_G + (integraldeltaCp_rxn_4);        %Heat of reaction: reaction 4
%INITIAL CONCENTRATIONS  
VOLUMETRICFLOWRATEin = F_FEED*R*T0;                         %Initial volumetric flowrate into the reactor (m^3/s) (assuming ideal gas in)   NEGLECTING P/P0 to start...
C_G0 =      F_G0/VOLUMETRICFLOWRATEin;                      %Inital concentration of G
C_H2O0 =    F_H2O0/VOLUMETRICFLOWRATEin;                    %Inital concentration of H2O
C_N20 =     F_N20/VOLUMETRICFLOWRATEin;                     %Inital concentration of N2
C_A0 =  0;                                                  %Outlet concentration of A

C_AC0 = 0;                                                  %Outlet concentration of AC

C_AA0 = 0;                                                  %Outlet concentration of AA

C_FA0 = 0;                                                  %Outlet concentration of FA

C_T0 =      F_FEED/VOLUMETRICFLOWRATEin;                    %Initial total concentration
%FINAL Concentrations
F_T = F(1) + F(2) + F(3) + F(4) + F(5) + F(6) + F(7);       
C_G =   C_T0*(F(1)/F_T)*(T0/F(8));    %Outlet concentration of G              
C_H2O = C_T0*(F(2)/F_T)*(T0/F(8));      %Outlet concentration of H2O
C_N2 =  C_T0*(F(3)/F_T)*(T0/F(8));      %Outlet concentration of N2              NEGLECTING P/P0 to start...
C_A =   C_T0*(F(4)/F_T)*(T0/F(8));    %Outlet concentration of A
C_AC =  C_T0*(F(5)/F_T)*(T0/F(8));      %Outlet concentration of AC
C_AA =  C_T0*(F(6)/F_T)*(T0/F(8));      %Outlet concentration of AA
C_FA =  C_T0*(F(7)/F_T)*(T0/F(8));      %Outlet concentration of FA
%k-VALUES 
k1 = 23.3244*exp(-45007.1/(R*F(8)));           %Reaction rate constant k1 (m^3/kgcat*s)
k2 = 2528.8*exp(-79244.89/(R*F(8)));           %Reaction rate constant k2 (m^3/kgcat*s)
k3 = 2150.595*exp(-73702.78/(R*F(8)));         %Reaction rate constant k3 (m^3/kgcat*s)
k4 = 0.146768*exp(-24596.97/(R*F(8)));         %Reaction rate constant k4 (m^3/kgcat*s)
%reaction rates
r_1 = -k1*C_G;                                     %Rate of reaction 1
r_2 = -k2*C_G;                                     %Rate of reaction 2
r_3 = -k3*C_G;                                     %Rate of reaction 3
r_4 = -k4*C_G;                                     %Rate of reaction 4
r_G = -(k1+k2+k3+k4)*C_G;                  %Rate of formation of G
r_N2 = 0;                                  %Rate of formation of N2
r_H2O = k1*(C_G)^2+(k3+k4)*G_G;            %Rate of formation of H2O
r_A = k1*C_G;                              %Rate of formation of A
r_AC = k2*C_G;                             %Rate of formation of AC
r_AA = (k3+k4)*C_G;                        %Rate of formation of AA
r_FA = (k3+k4)*C_G;                        %Rate of formation of FA
%Differential Equations to Solve
dFdW(1) = (r_G);
dFdW(2) = (r_H2O);
dFdW(3) = (r_N2);
dFdW(4) = (r_A);
dFdW(5) = (r_AC);
dFdW(6) = (r_AA);
dFdW(7) = (r_FA);
dFdW(8) = (Ua(Ta-T0)+(-r_1)*(-Hrxn_1)+(-r_2)*(-Hrxn_2)+(-r_3)*(-Hrxn_3)+(-r_4)*(-Hrxn_4))/(F(1)*Cp_G + F(2)*Cp_H2O + F(3)*Cp_N2 + F(4)*Cp_A + F(5)*Cp_AC + F(6)*Cp_AA + F(7)*Cp_FA); %dTdW from pg 545 in Fogler;   
end

Best Answer

Related Solutions

MATLAB: Error when using ode15

MATLAB: What could be the issue if a global variable gets an error “Undefined function or variable ‘F_FEED’ ” when called from a function. The global variable is defined in the file from which the function is called.

Related Question