Annexe 1.
Code Matlab
Fonction générale : Bilan_total.m
% Calcul de la surface de capteur nécessaire à maintenir la température de
% la piscine constante. Cette fonction dépend de plusieurs paramètres :
% * L'orientation et l'inclinaison du capteur (or et i)
% * Le débit dans le capteur (D) en Litres par heure
% * La surface de piscine considérée et la température voulue (Spiscine et
% Ti)
% * Les conditions extérieures : vitesse du vent, taux d'humidité,
% température de l'air ambiant (vent, humidite, Text)
%
% Un dernier paramètre sera le jour considéré.
function []= Bilan_total (TIME, i, or, Spiscine, D, vent, humidite, Ti, Text)
Ti = Ti + 273.15; % Les températures, entrées en °C, sont ici exprimées en Kelvin
Text = Text + 273.15;
D=D/3600; % Débit en kg/s sachant que la masse volumique de l'eau est de 1 kg/L
% Constantes
cp=4186; % Capacité calorifique de l'eau en J/kg.K
sigma = 5.68 * 10^-8; % Constante de Stephan-Boltzmann exprimée en J/m².s.K^4
Eair = 0.8 ; % Emissivité de l’air
% Calcul du bilan de la piscine, dépendant de la fonction "bilan_piscine" et
% de la surface du capteur. La fonction exprime l'énergie totale d'un mètre
% carré d'eau en faisant la différence du flux arrivant sur la surface et
% des pertes par évaporation, convection et rayonnement.
Episcine = Spiscine * bilan_piscine (TIME, vent, humidite, Ti, Text);
if Episcine > 0 % Si le bilan de la piscine est positif, il n'est pas nécessaire d'avoir un capteur
Scapteur = 0;
disp ('Aucune surface de capteur n''est nécessaire pour chauffer la piscine')
else
% Pour maintenir la température de l'eau constante, il nous faut
% Ecapteur + Episcine = 0. De plus, Ecapteur = D*cp * (Ts-Ti).
Ts = Ti - Episcine / (D*cp);
% Connaissant la température de sortie de l'eau du capteur (Ts), nous
% pouvons calculer la température de l'absorbeur et les pertes par
% rayonnement associées.
Tabs = (Ts+Ti)/2;
PCapteur = sigma * (Tabs^4 - Eair * Text^4); % Pertes par rayonnement du capteur sans
% vitre
F = flux(TIME,i,or);
if (F - PCapteur) < 0 % Si les pertes du capteur sont
% supérieures aux apports par énergie solaire, le capteur ne pourra
% pas chauffer l'eau de la piscine.
fprintf ('Le flux arrivant sur le capteur est de %f W/m². \n', flux(TIME,i,or));
fprintf ('La température de sortie de l''eau sera de %d °C.\n',round(Ts-273.15));
fprintf ('La température de l''absorbeur sera de %d °C. \n', round(Tabs-273.15));
fprintf ('Les pertes de l''absorbeur seront de %f W/m².\n',PCapteur);
disp ('Les pertes dans le capteur sont trop importantes pour pouvoir chauffer l''eau de la
piscine')
else
% L'énergie du capteur = Scapteur / (Flux arrivant sur la surface du
% capteur - Pertes par rayonnement).
Scapteur = -Episcine/(F-PCapteur);
fprintf ('Le flux arrivant sur le capteur est de %f W/m². \n', F);
fprintf ('La température de sortie de l''eau sera de %d °C.\n',round(Ts-273.15));
fprintf ('La température de l''absorbeur sera de %d °C. \n', round(Tabs-273.15));
fprintf ('Les pertes de l''absorbeur seront de %f W/m².\n',PCapteur);
fprintf ('La surface de capteur nécessaire pour chauffer la piscine est alors de \n %f m².
\n', Scapteur);
end
end
Pertes de la piscine : bilan_piscine.m
% Calcul du bilan thermique d'un plan d'eau en fonction des conditions
% extérieures, du jour considéré et de la température de l'eau voulue.
function [total] = bilan_piscine(TIME, vent, humid_air, Teau, Tair);
% Constantes :
h = 5.5; % Coefficient d'échange thermique, enW/m²K.
sigma = 5.68 * 10^-8; % Constante de Stephan-Boltzmann J/m².s.K^4.
Eair = 0.8; % Emissivité de l'air
%%%%%%%%%%%
%% Evaporation %%
%%%%%%%%%%%
% Calcul de la pression de vapeur saturante à une température donnée en
% fonction d'une valeur de référence.
% ln (p_sat_25) + DHvap / R * (1/T_25 - 1/T)
p_sat_25 = 3.17 * 10^3 ; % Pression de vapeur saturante à 25°C (en Pascal).
DHvap = 43.9 * 10^3; % Chaleur latente à 25°C (J/mol).
R = 8.315 ; % Constante des gaz parfaits (J/mol K).
Psat = inline('exp (log (3.17 * 10^3) + 43.9 * 10^3 / 8.315 * (1/298 - 1/T))', 'T');
% Où log est le logarithme népérien.
% Les pertes par évaporation du plan d'eau sont donc : (en valeur absolue et en W/m²)
Qevap = 0.8 * (0.089 + 0.078 * vent) * (Psat(Teau) - humid_air * Psat(Tair));
%%%%%%%%%%
%% Convection %%
%%%%%%%%%%
Flux_chaleur = h * (Teau - Tair); % Pertes par convection exprimées en W/m².
%%%%%%%%%%%
%% Rayonnement %%
%%%%%%%%%%%
Qray = sigma * (Teau^4 - Eair * Tair^4); % W/m²
%%%%%%%
%% Total %%
%%%%%%%
% Le bilan total s'exprime par la différence entre l'énergie reçue et les
% pertes de la piscine :
total = flux(TIME,0) - (Qevap + Flux_chaleur + Qray); % W/m²
fprintf('\nLes pertes de la piscine sont de %f W/m².', - (Qevap + Flux_chaleur+Qray))
fprintf('\nLe bilan total de la piscine est de %f W/m². \n', total)
Flux d’énergie solaire à Bruxelles : flux.m
% Calcul du flux solaire à un moment donné dans une zone urbaine de faible altitude.
% Paramètre à introduire : inclinaison de la surface(par rapport à
% l'horizontale donc 90 pour une surface verticale), l'orientation par
% rapport au sud, compté positivement vers l'ouest (par défaut, plein sud)
% Retourne le flux total au moment précis (temps donné) et l'angle entre la
% normale au panneau et les rayons solaires.
function [G] = flux(TIME, i,or)
%Par défaut, orientation plein sud
if nargin == 2
or=0;
end
%Vérifie les valeurs introduites
if i>90 || i<0
error('L''angle par rapport à l''horizontale doit être compris entre 0 et 90 °')
end
if or>180 || or<0
error('L''orientation doit être compris entre 0 et 180 °')
end
% Initialisation des valeurs du lieu étudié
lat = 50 + 50/60; %Latitude de Bruxelles : 50°50'
lat = deg2rad(lat); %On passe tous les angles en radian pour les calculs avec les fonctions
trigonométriques.
%On définit la longueur d'une journée
day = [23 56 4];
% Pas nécessaire mais permet d'avoir l'heure d'hiver
TIME(1,4)= TIME(1,4)-1;
% On calule le nombre de jours depuis le premier janvier à l'aide de la
% fonction datenum
J = datenum(TIME(1,1),TIME(1,2),TIME(1,3))-datenum(TIME(1,1),1,1) + 1;
J = J + TIME(1,4)/24 + TIME(1,5)/(24*60); % On tient compte des heures, minutes pour le
"jour"
J = J*360/365.25;
J = deg2rad(J);
% On calcule la declinaison en fonction du jour
dec = 0.32281 - 22.984*cos(J) - 0.34990*cos(2*J)- 0.13980*cos(3*J) + 3.7872*sin(J) +
0.03205*sin(2*J) + 0.07187*sin(3*J);
dec = deg2rad(dec);
% Angle solaire AH. Approximation en accord avec le cahier des charges
AH = -180+ (heuretonum(TIME(1,[4 5 6])) / heuretonum(day) * 360);
AH = deg2rad(AH);
% Hauteur du soleil
h = asin (sin(lat)*sin(dec) + cos(lat)*cos(dec)*cos(AH));
% Azimut du soleil
a = asin( cos(dec)*sin(AH)/cos(h) );
if h>0
% Initialisation des valeurs propres au panneau
i = deg2rad(i);
or = deg2rad(or);
%%%%%%%%%%%%%%%%%
% Calcul du rayonnement direct %
%%%%%%%%%%%%%%%%%
% Calcul de l'angle entre la normale au panneau et les rayons
% solaires
beta = acos( cos(h)*sin(i)*cos(a - or) + sin(h)*cos(i) );
cos(beta);
% Calcul du flux arrivant à la surface de la terre sur une surface
% perpendiculaire au flux
I = 1260*exp(-1/(2.3*(h+3)*sin(h)));
% Calcul du flux (S) sur la surface inclinée
if cos(beta) <= 0
S = 0;
else
S = I*cos(beta);
end
%%%%%%%%%%%%%%%%%
% Calcul du rayonnement diffus %
%%%%%%%%%%%%%%%%%
% Initialisation des paramètres
alb = 0.2; % Albédo
% Calcul du rayonnement diffus (D)
D = 125*(sin(h)^0.4)*( (1+cos(i))/2 ) + alb * 1080* (sin(h)^1.22) * ( (1-cos(i))/2);
%%%%%%%%%%%%%%%%%
% Calcul du rayonnement global %
%%%%%%%%%%%%%%%%%
G = S + D;
else
G=0
end
heuretonum.m
% Cette fonction permet de convertir une heure sous format d'une matrice du
% type [H M S] en heure décimale
function a = heuretonum(G)
% Vérifie la taille de la matrice introduite
if size(G) ~= [1 3]
error('La taille de la matrice introduite n''est pas bonne')
end
% Fait la convertion et retourne le résultat
a = G(1) + G(2)/60 + G(3)/3600
Lever et coucher du soleil
% Script dans le cadre du projet BA2
% Calcul du flux solaire à un moment donné dans une zone urbaine de faible altitude.
% Initialisation des valeurs
lat = 50 + 50/60; % Latitude de Bruxelles : 50°50'
lat = degtorad(lat); % On passe tous les angles en radian pour les calculs avec les fonctions
trigonométriques.
% Initialisation des valeurs propres au panneau
i = 10; % Inclinaison du panneau par rapport à l'horizontale
i = degtorad(i);
or = 0; % Orientation du panneau (0° pour le sud, compté positivement vers l'ouest)
or = degtorad(or);
% On définit la longueur d'une journée
day = [23 56 4];
% On fixe le jour pour lequel on fait les calculs [Y M D H M S]
TIME = [2006 6 21 16 0 0];
% On calcule le nombre de jours depuis le premier janvier à l'aide de la
% fonction datenum
jour = datenum(TIME(1,1),TIME(1,2),TIME(1,3))-datenum(TIME(1,1),1,1) + 1; % On
calcule le jour de l'année
J = jour;
J = J*360/365.25;
J = degtorad(J);
% On initialise un compteur qui s'incrémente à chaque passage dans la boucle
k=0;
for heure=0:23
for min = [Link]
k = k+1;
TIME(1,[4 5 6]) = [heure min 0];
% On recalcule le jour en tenant compte des heures, minutes
J = jour + TIME(1,4)/24 + TIME(1,5)/(24*60);
J = J*360/365.25;
J = degtorad(J);
% On calcule la déclinaison en fonction du jour
dec = 0.32281 - 22.984*cos(J) - 0.34990*cos(2*J)- 0.13980*cos(3*J) + 3.7872*sin(J) +
0.03205*sin(2*J) + 0.07187*sin(3*J);
dec = degtorad(dec);
% Angle solaire AH. Approximation en accord avec le cahier des charges
AH = -180+(((heuretonum(TIME(1,[4 5 6])))/heuretonum(day))*360);
AH = degtorad(AH);
% Hauteur du soleil
h = asin (sin(lat)*sin(dec) + cos(lat)*cos(dec)*cos(AH));
hauteur(k)=radtodeg(h);
% Azimut du soleil
a = asin( cos(dec)*sin(AH)/cos(h) );
azimut(k)=radtodeg(a);
%%%%%%%%%%%%%%%%%
% Calcul du rayonnement direct %
%%%%%%%%%%%%%%%%%
% Calcul de l'angle entre la normale au panneau et les rayons solaires
beta = acos( cos(h)*sin(i)*cos(a - or) + sin(h)*cos(i) );
% Calcul du flux arrivant à la surface de la terre sur une surface
% perpendiculaire au flux
if radtodeg(h)>0
I = 1260*exp(-1./(2.3*(h+3)*sin(h)));
else I = 0;
end
% Calcul du flux (S) sur la surface inclinée
if cos(beta) <= 0
S = 0;
else
S = I*cos(beta);
end
fdirect(k) = S;
%%%%%%%%%%%%%%%%%
% Calcul du rayonnement diffus %
%%%%%%%%%%%%%%%%%
if h>0
% Initialisation des paramètres
alb = 0.2; %Albédo
% Calcul du rayonnement diffus (D)
D = 125*(sin(h)^0.4)*( (1+cos(i))/2 ) + alb * 1080* (sin(h)^1.22) * ( (1-cos(i))/2);
else D=0;
end
fdiffus(k) = D;
%%%%%%%%%%%%%%%%%
% Calcul du rayonnement global %
%%%%%%%%%%%%%%%%%
G = S + D;
flux(k) = G;
end
end
t = linspace(heuretonum([0 0 0]),heuretonum([23 50 0]),k);
plot(t,hauteur)
hold on
plot([0 23],[0 0],'r')
tmp = ephem(TIME);
plot(heuretonum(tmp(1,[4 5 6])),0,'ro')
plot(heuretonum(tmp(2,[4 5 6])),0,'ro')
title('Hauteur du soleil en fonction de l''heure')
xlabel('Heure')
ylabel('Hauteur en degré')
grid on
figure
plot(t,azimut)
title('Azimut du soleil en fonction de l''heure')
xlabel('Heure')
ylabel('Azimut en degré')
grid on
figure
subplot(2,2,[1 2])
plot(t,flux)
grid on
title('Valeur du rayonnement direct et diffus en fonction de l''heure')
xlabel('Heure')
ylabel('Rayonnement global en W/m²')
subplot(2,2,3)
plot(t,fdirect)
grid on
title('Valeur du rayonnement direct en fonction de l''heure')
xlabel('Heure')
ylabel('Rayonnement direct en W/m²')
subplot(2,2,4)
plot(t,fdiffus)
grid on
title('Valeur du rayonnement diffus en fonction de l''heure')
xlabel('Heure')
ylabel('Rayonnement diffus en W/m²')
