clc
clear all
close all


%user defined variables
Tl = 20;            %°C
To = 20;            %°C
R_tank = 0.5;       %m
S_heater = 0.06;    %m2

%constants
ml = 1000;  %kg
cl = 4180;  %J/kg*K
mh = 5;     %kg
ch = 445;   %J/kg*K


%tank dimensions
V_tank = 1;                                          %m3
v_tank = V_tank/(pi*R_tank*R_tank);                  %m
S_tank = 2*pi*R_tank*R_tank + v_tank*2*pi*R_tank;    %m2

%temp simulation
time = 3600*16;    %s
P = 12000;         %W
Qsum = 0;          %J


%heat transfer coeficients
Kl = 11.3*S_tank;       %W/°C
Kh = 400*S_heater;      %W/°C

Temp_plot = 1:time;
Temp_plot2 = 1:time;
Tv = 100;
Td = 300;
Inflow = 6;
InflowT = 20;

t1 = 3600*5;
t2 = 3600*7;
t3 = 3600*10;

    
    
for t=1:Tv:time
    
    %heating - temp increase
    if t<t1
        a2 = (mh*ch*ml*cl)/Kh;
        a1 = mh*ch*((Kl/Kh)+1) + ml*cl;
        a0 = Kl;
        x = roots([a2 a1 a0]);
        a = -x(1);
        b = -x(2);
        A = P/a0;
        C = ((P/a2)*(b*a2 - a1))/(a0*(a-b));
        B = (-P/a0)-C;
        
        T_increment = A + B*exp(-a*t) + C*exp(-b*t);
        Tl = To+T_increment;
        Qsum = ml*cl*T_increment;         
               
    %no action (water cools through tank) - small temp decrease  
    elseif t>t1 && t<t2
        Qout = (Kl/S_tank)*Tv*(Tl-To);
        Qsum = Qsum - Qout;
        Tl = (Qsum/(ml*cl)) + To;
      
        
    
    %inflow of cooler water into tank - temp decrease
    elseif t>t2 && t < t3
        ml = ml + Inflow;
        Tl = ((Inflow*InflowT)+(ml*Tl))/(Inflow+ml);
        Tonew = Tl;
        
    %heating + inwlow of cooler water - small temp increase
    else
        ml = ml + Inflow;
        a2 = (mh*ch*ml*cl)/Kh;
        a1 = mh*ch*((Kl/Kh)+1) + ml*cl;
        a0 = Kl;
        x = roots([a2 a1 a0]);
        a = -x(1);
        b = -x(2);
        A = P/a0;
        C = ((P/a2)*(b*a2 - a1))/(a0*(a-b));
        B = (-P/a0)-C;
        
        T_increment = A + B*exp(-a*(t-t3)) + C*exp(-b*(t-t3));  %-t3 because heating process must start from t=0s
        Tl = Tonew+T_increment;
        Tl = ((Inflow*InflowT)+(ml*Tl))/(Inflow+ml);
    end
   
    Temp_plot(1,t) = Tl;
      
       
end

figure
for t=1:Tv:time
    plot(t,Temp_plot(1,t), 'b.')
    hold on
   
end
xlabel('time[s]')
ylabel('temperature[°C]')
axis([-100 time 20 70])
grid







 




