1. Limited time only! Sign up for a free 30min personal tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Concentration of CO in Los Angeles Basin (matlab)

  1. Sep 10, 2014 #1

    Maylis

    User Avatar
    Gold Member

    1. The problem statement, all variables and given/known data
    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*10^{12} \frac{ft^3}{hr}##. 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 ##2*10^{-10} \frac{lb mol}{ft^3}##.

    Given Information:

    [tex]F_{CO,A} + F_{CO,S} - Q*C_{CO} = V\frac{dC_{CO}}{dt}[/tex]
    [tex]F_{CO,A} = a + b sin(\pi*t/6)[/tex]
    Where ##a = 20254815 \frac {lb mol}{hr}## and ##b = 0.5*a \frac {lb mol}{hr}##
    [tex]V=4*10^{13} ft^3[/tex]
    ##Q=1.67*10^{13} \frac{ft^3}{hr}## for ##t \lt 72## and ##Q= 1.67*10^{12} \frac{ft^3}{hr}## for ##t \ge 72##
    [tex]F_{CO,S}=C_{CO,S}*Q[/tex]
    [tex]C_{CO,S} = 2*10^{-10} \frac{lb mol}{ft^3}[/tex]
    Initial condition: [tex]C_{CO,t=0} = 2*10^{-10} \frac{lb mol}{ft^3}[/tex]

    2. Relevant equations
    The crux of my problem is using ode45 to solve this differential equation
    ##F_{CO,A} + F_{CO,S} - Q*C_{CO} = V\frac{dC_{CO}}{dt}##

    3. The attempt at a solution
    Hello, I want to give a little context for this problem so that you can understand what is going on. ##F_{CO,S}## is the flow rate of carbon monoxide entering the basin with a concentration of carbon monoxide ##C_{CO,s}## from the Santa Ana winds. ##F_{CO,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\frac{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 \ge 72##, and I'm not sure how to remedy this in the ode45 built in function. It is that ##Q## is a piecewise function of ##t##, as well as ##C_{CO} = f(t)##. Anyways, here is what I wrote so far. I know the function handle is supposed to have @(t,x) for a derivative ##\frac {dx}{dt}##, so I used in my function handle ##@(t,C_{CO})## for ##\frac {dC_{CO}}{dt}##. I keep getting the error Undefined function or variable 't'.


    Code (Text):
    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');
     
    Last edited: Sep 10, 2014
  2. jcsd
  3. Sep 10, 2014 #2

    DrClaude

    User Avatar

    Staff: Mentor

    What you should do is define your function using "function", such that you can incorporate the code defining Q
    Code (Text):

    function dC_COdt = fHan (t, C_CO)

       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

       dC_COdt = (F_COa + F_COs - C_CO*Q)/Volume;

    end
     
     
  4. Sep 11, 2014 #3

    Maylis

    User Avatar
    Gold Member

    I'm confused how to make it work how you are describing

    Code (Text):
    function dC_COdt = fHan (t, C_CO)

       if t < 72
           Q = @(t) 1.67e13; % ft^3/hr
       else
           Q = @(t) 1.67e12; % ft^3/hr
       end

       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

       dC_COdt = @(t,C_CO)(F_COa(t) + F_COs(t) - C_CO*Q(t))/Volume;

    end
    [T, C_CO] = ode45(dC_COdt,tSpan,C_COs);
    plot(T,C_CO)
    xlabel('time (hours)');
    ylabel('CO Concentration (lb mol/ft^3)');
    title('CO Concentration vs. time');
    Will not work, I get this error

    Code (Text):
    fHan
    Error: File: fHan.m Line: 21 Column: 1
    This statement is not inside any function.
     (It follows the END that terminates the definition of the
     function "fHan".)
     
    Last edited: Sep 11, 2014
  5. Sep 11, 2014 #4

    DrClaude

    User Avatar

    Staff: Mentor

    The file fHan.m should contain only the function, starting from the "function" statement and terminating with the "end".

    This function will then be called by your main program. You still need to initialize the time and set the initial conditions. Note also that the function is called "fHan", not "dC_COdt". You simply need to remove some lines from your original posting, leaving:
    Code (Text):

    tSpan = [0 100]; % hr

    C_COs = 2e-10; % lb mol/ft^3

    [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');

     
     
  6. Sep 11, 2014 #5

    Maylis

    User Avatar
    Gold Member

    Something isn't right still, I just copy paste your code in #2

    Code (Text):
    function dC_COdt = fHan (t, C_CO)

       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

       dC_COdt = (F_COa + F_COs - C_CO*Q)/Volume;

    end
    Then make a separate script, and I put in arguments for fHan in my ode45 call,
    Code (Text):
    tSpan = [0 100]; % hr

    C_COs = 2e-10; % lb mol/ft^3

    [T, C_CO] = ode45(fHan(0,C_COs),tSpan,C_COs);
    plot(T,C_CO)
    xlabel('time (hours)');
    ylabel('CO Concentration (lb mol/ft^3)');
    title('CO Concentration vs. time');
    But I get the error
    Code (Text):
    Error using exist
    The first input to exist must be a string.

    Error in odearguments (line 60)
      if (exist(ode)==2)

    Error in ode45 (line 114)
    [neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs,
    odeFcn, ...
    the ode45 function doesn't seem to like the way I input the first argument
     
  7. Sep 11, 2014 #6

    DrClaude

    User Avatar

    Staff: Mentor

    Do not put arguments, it's ode45 which supplies the correct arguments.

    What was missing in my post is to indicate that fHan is a function:
    Code (Text):
    [T, C_CO] = ode45(@fHan,tSpan,C_COs);
     
  8. Sep 11, 2014 #7

    Maylis

    User Avatar
    Gold Member

    Thanks, ode45 arguments behave somewhat mysteriously to me. It just magically knows what to put it
     
  9. Sep 11, 2014 #8

    DrClaude

    User Avatar

    Staff: Mentor

    Of course, there is nothing magical about it :smile:

    It is basically a stepper method: starting from some initial conditions ##y_0## at ##x_0##, it uses a discrete form of the differential equation to find the solution ##y_1## at ##x_1##. For that, it needs to know the differential equations ##dy/dx##, which is what you supply. In the simplest case, you would have
    $$
    y_1 = y_0 + (x_1 - x_0) \left. \frac{dy}{dx} \right|_{x_0}
    $$
    which the algorithm can calculate by feeding ##x_0## and ##y_0## into your function.

    In reality, the algorithm used by ode45, the Runge-Kutta-Fehlberg 4-5 method, is a more sophisticated version of the equation above with a smaller error. In addition, it can make an estimate of the error on ##y_1##, and it reduces it by subdividing the interval between ##x_0## and ##x_1## into smaller intervals, until a desired accuracy is obtained. For that, it will need to call your function many times to get ##dy/dx## at different values of ##x##. In the end, ode45 returns the values of ##x## for all the intermediate steps it took, along with the corresponding value of ##y(x)##.

    A good introduction to the subject can be found in the book Numerical Recipes in C, see http://apps.nrbook.com/c/index.html, starting at page 707.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted



Similar Discussions: Concentration of CO in Los Angeles Basin (matlab)
  1. Stress concentration (Replies: 1)

  2. Stress concentrator (Replies: 3)

Loading...