%Control of floor heating process

%Simple model

function [x0,par]=getData

%% Simulation parameters

%Data

l_wall = 2.4; %[m]
l_floor = 5; %[m]
A_floor = 25; %[m^2]

x_wood = 0.015; %[m]  %Assumption: Wood = oak
x_rock_wool = 0.25; %[m]
x_concrete = 0.05; %[m]

%Densitiy
d_concrete = 2400; %[kg/m^3]
d_wood = 825; %[kg/m^3]
d_rock_wool = 192; %[kg/m^3]
d_air = 1.166; %[kg/m^3] average between 20 and 40C 

%Molcular weight
Mr = 28.97/1000; %[kg/mol] Air. Assume 78% N2, 21% O2 and 1% Ar

co2air=0.04/100; %[vol %]
N2=0.78;
O2=0.21;

wco2=44.01/1000; % [kg/mol]
wN2=28/1000; % [kg/mol]
wO2=32/1000; % [kg/mol]
co2in=co2air/(co2air+N2*wN2/wco2+O2*wO2/wco2) ;%[kg/kg];


co2max=0.001/(0.001+N2*wN2/wco2+O2*wO2/wco2);
breathrate=8e-3/60; %[m^3/s]-> [4-8] l/min at rest (per person)
co2breath=0.04;%[ %Vol of CO2/s]
co2wbreath=co2breath/(co2breath+N2*wN2/wco2+O2*wO2/wco2);%[mass fraction CO2 produced by breathing];
co2rate=co2wbreath*breathrate*d_air; %[Kg CO2/ s]


%Heat capacity
Cp_concrete = 0.63; %[kJ/kgK]  
Cpr = 1.005; % [kJ/kgK] Cp_air at -50 - 40C
Cpf = Cp_concrete;

%Heat transfer coefficients
hair = 15/10^3; %[kJ/(s x m^2 x K)] (Still air)
k_wood = 0.208/10^3; %[kJ/(s x m x K)] Oak
k_rock_wool = 0.0391/10^3; %[kJ/(s x m x K)]  At 37.8C
k_concrete = 0.762/10^3; %[kJ/(s x m x K)]

%Other
R = 8.314/(10^5); %[kJ / K x mol]
Po = 1.01; %[bar]
k = 100;
% ------------------------------------------------------------------------

%Mass = density x volume
%Floor
m_concrete = d_concrete*(x_concrete*A_floor);
%m_wood = d_wood*(x_wood*A_floor);
mf = m_concrete;
%Room
%m_wood1 = d_wood*(x_wood*l_wall*l_floor);
%m_rock_wool = d_rock_wool*(x_rock_wool*l_wall*l_floor);
mr = d_air*(l_wall*A_floor);
Vr=l_wall*A_floor; %m^3
% ------------------------------------------------------------------------

%Floor
Af = A_floor; % [m^2]
Ufr = 1/(1/hair+x_wood/k_wood); %[KJ/(s x m^2 x K)]

UAfr = Ufr*Af;
mCpf = mf*Cpf;

%Room
Ar = 4*l_wall*l_floor; % [m^2] Varme  fra fire vegger
Uro = 1/((1/hair)+(x_wood/k_wood)+(x_rock_wool/k_rock_wool)+...
 (x_wood/k_wood)+(1/hair)) ;%[J/(h x m^2 x K)] 
 
% Change by Vinicius- Increased by factor of 10
UAro = 3*Uro*Ar;
mCpr = mr*Cpr;

Vr = l_wall*A_floor ;%[m^3]

%% Inputs and initial conditions

%Inputs
qhr = 2; % No heat on: 0 %max: 2000*60*60; % [KJ/s]
qhf = 0; %No heat on: 0 %max: 575*60; % [J/h]
qin = 22/3600*1.166; % [Kg/s] typical: q=22m^3/h= 22/3600 m^3/s => 22/3600*1.166 Kg/s; 
To = 273; % [K]

%Initial conditions
Tf = 273+17; % [K] value at steady state with no heating on
Tr = 273+17; % [K] value at steady state with no heating on
mr = d_air*(l_wall*A_floor); % [kg]
init = [Tf;Tr;mr];
init2=[Tf;Tr;mr;co2in];
%Set point
Trset = 293;
Tmin = 283; %[K]
Tmax = 293; %[K]

%% Simulation

% Simulation parameters
par.UAfr=UAfr;
par.UAro=UAro;
par.Cpr=Cpr;
par.Cpf=Cpf;
par.Po=Po;
par.mf=mf;
par.Mr=Mr;
par.R=R;
par.k=k;
par.Vr=Vr;
par.co2in=co2in;
par.co2rate=co2rate;
par.co2max=co2max;

%states + initial cost of energy
x0=[init2 ;0];

% sim('simulink_advancedmdl')

%Outputs
% Tf = simout(:,1);
% Tr = simout(:,2);
% mr = simout(:,3);


%% Plotting
% 
% Subplot(3,1,1)
% plot(t,Tf,'b')
% ylabel('Temperature [K]')
% xlabel('Time [h]')
% title('Tf') 
% 
% Subplot(3,1,2)
% plot(t,Tr,'r')
% ylabel('Temperature [K]')
% xlabel('Time [h]')
% title('Tr') 
% 
% Subplot(3,1,3)
% plot(t,mr,'g')
% ylabel('Mass [kg]')
% xlabel('Time [h]')
% title('mr')
% 
% 






