Using ODE45 to solve 2 dependent variable : dQ (solar radiation at specific time of the day(hour)) and dT/dt (temperature change of a particle ) of a particle during solar drying

Hi,

I am trying to use differential equation to solve how temperature changes with respect to time however the temperature is dependent on the solar radiation per unit area Q

so the formula for Q is given as:

    Time2=(1:t:12); %hours
    Q=960*(sin(2*pi*Time2/24)).^2; %W/m2

The equation for the temperature change is :

    dT/dt = ((Q*A)-(mw*lw))/(m*cp));

where

    T0 = 373.15;   %Initial temperature in Kelvin
    mw = 0.706;   
    m = 15; % Mass of product to dry (Kg)
    lw = 2260000; % Latent heat of vaporisation of water (J/Kg)
    A = 1;      % Surface Area of Collector (m^2)
    cp= 3746; % Specific heat capacity of product (J/Kg/K)

so at every time step (which is in hours ), I want the equation dT/dt to use the value of Q at that time to solve dT/dt

I tried this

    deltaT =@(t,T)((Q*A)-(mw*lw))/(m*cp);
    [t,T] = ode45('deltaT',tspan,T0);
where 

    tspan = [1 12] %hours

however I am getting so many errors

Any help will be appreciated, I am fairly new to MATLAB

Answer 1

[Torsten]

function main
  T0 = 373.15;    %Initial temperature in Kelvin
  tspan = [1 12]; %hours 
  dt = 0.1;
  Time2 = (1:dt:12); %hours
  Q = 960*(sin(2*pi*Time2/24)).^2; %W/m2
  [t,T] = ode45(@(t,T)deltaT(t,T,Time2,Q),tspan,T0);  
  plot(t,T)
end

function dT = deltaT(t,T,Time2,Q)
  Q_actual = interp1(Time2,Q,t);
  mw = 0.706;   
  m = 15; % Mass of product to dry (Kg)
  lw = 2260000; % Latent heat of vaporisation of water (J/Kg)
  A = 1;      % Surface Area of Collector (m^2)
  cp= 3746; % Specific heat capacity of product (J/Kg/K)
  dT = ((Q_actual*A)-(mw*lw))/(m*cp);
end

[Evans123]

Hi,

Thank you very much for the quick reply.

I expected the temperature to be proportional to the solar radiation however this is not the case since temperature is decreasing. for example I expect the temperature to peak at the 6th hour since solar radiation is highest at that time.

I may be thinking of this wrong but I will be glad if you could explain the reason behind your code.

[Torsten]

dT/dt is always negative - thus temperature will decrease.

This is your code - so you will have to ask yourself why the results you obtain are not physical.

[Evans123]

Hi Torsten,

This might be a long shot but I will really appreciate your help.

I'm currently working on a project to build a solar food dryer and I need to model on Matlab how temperature of the product will change with respect to change in solar radiation, Q .

Q is given by ;

   Q =960*(sin(2*pi*Time2/24)).^2; %W/m2
where 
    Time2=(1:t:12); %hours

The heat flow equation is given by

    Q(t)A = mcp*(T2-T1) + (mw*lw)

where :

    mw = 0.706; % Mass of water loss from product in hr (Kg/h)

    m = 15; % Mass of product to dry (Kg)

    lw = 2260000; % Latent heat of vaporisation of water (J/Kg)

    A = 1; % Surface Area of collector (m^2)

    cp= 3746; % Specific heat capacity of product (J/Kg/C)

    T1 = temperature at t

    T2 = temperature at t + dt

Manipulating the heat flow give T2 as ;

    T2= (((Q*A*3600) -(mw*lw))/(m*cp)) + T1; 
    % 3600 is there to convert j/s to J/h

however implementing this on Matlab is proving a challenge for me that's why i went the ode45 route but that was clearly the wrong way to go about it.

This is what I have so far ;

    close all

    clear;

    mw = 0.706; % Mass of water loss from product in hr (Kg/h)

    m = 15; % Mass of product to dry (Kg)

    lw = 2260000; % Latent heat of vaporisation of water (J/Kg)

    A = 1; % Surface Area of Collector (m^2)

    cp= 3746; % Specific heat capacity of product (J/Kg/C)

    t = 1; % Time step

    T = 24; % Initial Temperature (°C)

    Time2=(1:t:12); hours

    Q=960*(sin(2*pi*Time2/24)).^2; % Solar irradiation (W/m2)

    for j = 1:12

     T(j+1)= ((((QQ2(j)*A*j*3600))-(mw*lw))/(m*cp))+ T(1);

    end

    figure(2)

    plot(T)

    title('Temperature')

    xlabel('time (hours)')

    ylabel('Temperature (°C)')

This seems wrong since the mass, m should decrease by mw after each hr and the temperature profile should follow the profile of the solar radiation. i.e peak at the same time

I have been spending days to get my head around this but i'm pretty bad at Matlab so I haven't made any meaningful progress. I would be grateful for your help.

[Torsten]

Your equation reads (with the values you gave inserted):

(15*3746 + m_Water(t)*cp_Water)*dT/dt = -0.706/3600*2260000 + 960*sin^2(2*pi*t/(24*3600))*A(t)

with t in seconds.

Solving this will give you more realistic temperature profiles.

So the error was primary related to the unit conversion.

Nevertheless, you should take into account the time-dependent mass of water you have to heat up in the m_Water(t)*cp_Water - term and the decreasing surface area over time in the A(t)-term.

Best wishes

Torsten.