The ode45 function only outputting constant value
10 views (last 30 days)
Show older comments
Braydon Lantz
on 25 Apr 2019
Commented: Braydon Lantz
on 26 Apr 2019
I am currently trying to create a program which accepts user inputs for parameters, passes those to a function, and then passes that function to the ode45 solver and plots the results. However, the ode solver seems to only give a constant output equal to the specified initial conditions. I was wondering whether there were some kind of error in my script that was leading to this result.
%% Units: g, mol, dm^3, K, Kpa, s
clear
clc
%Specifying operating conditions
R = 8.314472; %Ideal Gas Constant
m_a = 60.052; % Molar mass acetic acid
m_b = 74.1; % Molar mass Butanol
phi = 0.2; %Porosity of Catalyst
rho_c = 19300; %Density of Catalyst
P_0 = input('Please specify an operating pressure in Kpa.');
while P_0 <= 0
P_0 = input('Please specify an operating pressure greater than zero.');
end
T = input('Please specify an operating temperature in K.');
while T <= 0
T = input('Please specify an operating temperature greater than zero.');
end
Dp = input('Please specify a pipe diameter in dm.');
while Dp <= 0
Dp = input('Please specify a pipe diameter greater than zero.');
end
F_t_0 = input('Please specify incoming molar flow rate in mol/s.');
while F_t_0 <= 0
F_t_0 = input('Please specify an incoming molar flow rate greater than zero.');
end
x_b = input('Please specify mole fraction of Butanol in feed.');
while x_b <= 0 || x_b > 1
x_b = input('Please specify a mole fraction between zero and one.');
end
x_a = 1 - x_b; %Mole fraction of acetic acid
%Calculating parameters for pressure drop equation
Ap = pi*(Dp/2)^2; %Cross Sectional Area of Pipe
Cto = P_0/(R*T); %Initial concentration
rho_a = x_a*Cto*m_a; %Density of acetic acid
rho_b = x_b*Cto*m_b; %Density of butanol
rho = x_a*rho_a + x_b*rho_b; %Approximate total density
mu_a = 0.001155; %Viscosity of acetic acid
mu_b = 0.002593; %Viscosity of butanol
mu = x_a*mu_a + x_b*mu_b; %Approximate total viscosity
m_flow = x_a*F_t_0*m_a + x_b*F_t_0*m_b; %Incoming mass flow rate
G = m_flow/Ap; %Mass velocity
beta_1 = (G*(1-phi)/(rho*Dp*phi^3));
beta_2 = ((150*(1-phi)*mu)/Dp) + 1.75*G;
beta = beta_1*beta_2;
rho_bulk_c = rho_c*(1 - phi); %Bulk density of Catalyst
alpha = ((2*beta)/(Ap*rho_bulk_c*P_0));
%Calculating parameters for rate law
k_f = 0.0090953*exp(-45.49/(R*T));
k_r = 0.008643*exp(-23.90/(R*T));
%ODE solver
w_span = [0:1:100]; %In grams
F_a_0 = F_t_0*x_a;
F_b_0 = F_t_0*x_b;
[W,F] = ode45(@(w,y) f(w,y,Cto,k_f,k_r,F_t_0,alpha), w_span,[F_a_0; F_b_0; 0; 0; 1;]);
plot(W,F(:,1),W,F(:,2),W,F(:,3),W,F(:,4),W,F(:,5))
title('Component Flow Rates and Pressure Drop with Respect to Catalyst Mass');
xlabel('Mass W (g)');
legend('Fa','Fb','Fc','Fd','p')
function dydw = f(w,y,Cto,k_f,k_r,F_t_0,alpha)
dydw(1) = (Cto^2*y(5)^2/(y(1)+y(2)+y(3)+y(4))^2)*(k_f*y(1)*y(2) - k_r*y(3)*y(4)); %Mole Balance for Species A
dydw(2) = (Cto^2*y(5)^2/(y(1)+y(2)+y(3)+y(4))^2)*(k_f*y(1)*y(2) - k_r*y(3)*y(4)); %Mole Balance for Species B
dydw(3) = -1*(Cto^2*y(5)^2/(y(1)+y(2)+y(3)+y(4)^2)*(k_f*y(1)*y(2) - k_r*y(3)*y(4))); %Mole Balance for Species C
dydw(4) = -1*(Cto^2*y(5)^2/(y(1)+y(2)+y(3)+y(4))^2)*(k_f*y(1)*y(2) - k_r*y(3)*y(4)); %Mole Balance for Species D
dydw(5) = (-1*alpha/2*y(5))*((y(1)+y(2)+y(3)+y(4))/F_t_0); %Pressure Drop Equation
dydw = dydw(:);
end
2 Comments
Accepted Answer
Walter Roberson
on 26 Apr 2019
exp(-45.49/R*T) is exp(-1500-ish) which underflows to 0. Your input temperature would have to be below about 136 K in order for that expression to not underflow to 0.
Perhaps those should be exp(-45.49/(R*T)) ?
3 Comments
Walter Roberson
on 26 Apr 2019
What are some realistic inputs? With the random values I am inputting, I am seeing slight curves.
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!