MATLAB Answers

ODE45 problem with concentration of CO

1 view (last 30 days)
Rick
Rick on 10 Sep 2014
Commented: Star Strider on 11 Sep 2014
Smog begins to build up again immediately after a Santa Ana wind passes through the basin. The volumetric flow rate through the basin has dropped to 1.67 x 10^12 ft3/hr (1.5 mph). Plot the concentration of carbon monoxide in the basin as a function of time for up to 72 hours and just after a Santa Ana wind. The initial concentration of CO is 0.08 ppm (2 x 10^-10 lbmol/ft3).
Given F_CO,A + F_CO,S - Q*C_CO = V(dC_CO/dt)
F_CO,A = a + b sin(π*t/6)
Where a = 35,000 lb mol/h and b = 30,000 lb mol/h
V=4x10^13 ft3
Q =1.67x10^13 ft3/hr for t <72, Q= 1.67x10^12 ft3/hr for t ≥72
F_CO,S= C_CO,S*Q
C_CO,S = 2x10^-10 lb mol/ft3
Initial condition: C_CO(t=0)= 2x10^-10 lb mol/ft3
Hello, I want to give a little context for this problem so that you can understand what is going on. FCO,S is the flow rate of carbon monoxide entering the basin with a concentration of carbon monoxide (CO) CCO,s. FCO,A is the flow rate of carbon monoxide created by the cars driving inside the basin. Q*C_CO is the flow rate of carbon monoxide exiting the basin, and C_CO is changing with time because this is an unsteady state process. My problem is essentially solving the differential equation given in the problem
F_CO,A + F_CO,S - Q*C_CO = V(dC_CO/dt)
There are a couple of problems for me. First, there is the fact that Q changes depending on if t < 72 or if t >= 72, and I'm not sure how to remedy this in the ode45 built in function. Anyways, here is what I wrote so far. Keep in mind my a and b are different from the problem statement, but it should have no effect on where my issue lies here.
tSpan = [0 100]; % hr
if t < 72
Q = 1.67e13; % ft^3/hr
else
Q = 1.67e12; % ft^3/hr
end
C_COs = 2e-10; % lb mol/ft^3
F_COs = C_COs*Q; % lb mol/hr
a = 20254815; % lb mol/hr
b = 0.5*a; % lb mol/hr
F_COa = a + b*sin(pi*t/6); % lb mol/hr
Volume = 4e13; % ft^3
fHan = @(t,C_CO) (F_COa + F_COs - C_CO*Q)/Volume;
[T, C_CO] = ode45(fHan,tSpan,C_COs);
plot(T,C_CO)
xlabel('time (hours)');
ylabel('CO Concentration (lb mol/ft^3)');
title('CO Concentration vs. time');
and I get an error Undefined function or variable 't'. Well, I remember with ode45 you need to have your function handle be @(t,x). Also, how do I make it so that it knows with my if/else conditions how to handle the changing volumetric flow rate, Q?

Accepted Answer

Star Strider
Star Strider on 10 Sep 2014
What you want to do is a bit challenging with an anonymous function for your ODE function, but it can be done.
You need to create three functions of t:
Q = @(t) [1.67e13*(t<72) + 1.67e12*(t>=72)];
F_COs = @(t) C_COs*Q(t); % lb mol/hr
F_COa = @(t) a + b*sin(pi*t/6); % lb mol/hr
The Q function uses a combination of numerical and logical statements to force it to take on the appropriate values at various values of t. If (t < 72), that logical value equates to 1, so Q takes equates to 1.67e13. Since (t >= 72) is by definition false at the same value of t, that term evaluates to zero. The reverse occurs when (t >= 72). It’s a one-line if-then-else with only two conditions.
The other two lines simply create your previous statements as functions of t so your ODE function can call them with the values of t it gets from ode45.
Your ODE then becomes:
fHan = @(t,C_CO) (F_COa(t) + F_COs(t) - C_CO*Q(t))/Volume;
Your resulting code (with some statements commented-out) is:
tSpan = [0 100]; % hr
% if t < 72
% Q = 1.67e13; % ft^3/hr
% else
% Q = 1.67e12; % ft^3/hr
% end
Q = @(t) [1.67e13*(t<72) + 1.67e12*(t>=72)];
C_COs = 2e-10; % lb mol/ft^3
F_COs = @(t) C_COs*Q(t); % lb mol/hr
a = 20254815; % lb mol/hr
b = 0.5*a; % lb mol/hr
F_COa = @(t) a + b*sin(pi*t/6); % lb mol/hr
Volume = 4e13; % ft^3
% fHan = @(t,C_CO) (F_COa + F_COs - C_CO*Q)/Volume;
fHan = @(t,C_CO) (F_COa(t) + F_COs(t) - C_CO*Q(t))/Volume;
[T, C_CO] = ode45(fHan,tSpan,C_COs);
plot(T,C_CO)
xlabel('time (hours)');
ylabel('CO Concentration (lb mol/ft ^3)');
title('CO Concentration vs. time');
It runs, and seems to give the correct result, with the appropriate jump starting at (t = 72). Only you know if it does.
  3 Comments
Star Strider
Star Strider on 11 Sep 2014
I didn’t run your newest code, but with function files, you need to refer to them as function handles in the code that calls them. If they aren’t already declared as functions (with function handles) in your code, you have to declare them as such with a leading ‘@’ sign.
This should work:
[T, C_CO] = ode45(@fHan,tSpan,C_COs);
It’s not necessary to include the argument list since you aren’t passing other arguments to your ODE function ‘fHan’, and in the example you posted, t would always be 0 anyway. That’s not what you want.
If you were going to pass other arguments (in this illustration a vector of parameters called ‘parms’) to your ODE function using an anonymous function construction, you would do it this way:
[t,y] = ode45(@(t,y) odefun(t,y,parms),tspan,y0);

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!