1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
|
// Display mode
mode(0);
// Display warning for floating point exception
ieee(1);
max = 1000;
//nbr max d''iterations
l = 6;
//languer du mur
h = 3;
//hauteur du mur
e = 0.3;
//epaisseur du mur
n1 = 1/0.29;
//n=1/(r-1)//r=30
n2 = 1/2.99;
//n=1/(r-1)//r=300
n3 = 1/2.49;
//n=1/(r-1)//r=300-50=250
lambdab = 1.8;
//lambde du beton
rho = 2300;
// rho du beton
E1 = 925;
//cte solaire
emis = 0.6;
c = 800;
k = ones(30,300)*273.13;
//matrice qui sert à rammener les K en °C
sigma = 5.67*(10^(-8));
//cte de staphane-boltzmann
t0 = 278.13;
//T° ambiante de 5°C
tc = 308.13;
//chauffage à 35°C
p0 = (t0^4)*sigma;
s = l*h;
td = 14400;
//debut du chauffage à 10H
tf = 28800;
//fin du chauffage à 14H
t = 21600;
//instant auquel on veut connaitre la T° // à exprimer en secondes à partir de 6H du matin
omega = ones(30,300)*278.13;
//initialisation du domaine
//d2: face est:
a = (t*mtlb_double(%pi))/43300;
//secondes rammenés en radians //a: angle que fait le soleil avec la normale à la surface
if mtlb_logic(a,"<",mtlb_double(%pi)/2) then
// faces est exposé au rayonnement
t2 = mtlb_a(278.13,((((emis*E1)*sin(a))*cos(a))/(sigma*s))^0.25);
d2 = ones(1,300)*t2;
//est
else
t2 = 278.13;
d2 = ones(1,300)*t2;
//est
end;
//d5: face west en contact avec le radiateur:
if t<td then
t5 = 278.13;
d5 = ones(1,50)*t5;
//west/radiateur
elseif t<tf then
t5 = tc;
d5 = ones(1,50)*t5;
//west/radiateur
else
t5 = 278.13;
d5 = ones(1,50)*t5;
//west/radiateur
end;
//d1: base:
d1 = ones(30,1)*278.13;
//base
d1(1,1) = t5;
//base
d1(30,1) = t2;
//base
for m = 1:1:max
for i = 2:1:29
d1(i,1) = d1(i,1)+(-2*d1(i,1)+d1(i-1,1)+d1(i+1,1))/(2+((n1*rho)*c)/lambdab); // base
end;
end;
//d4: face west hors radiateur:
d4 = ones(1,250)*278.13;
//west/hors radiateur
d4(1,1) = t5;
//west/hors radiateur
d4(1,250) = 278.13;
//west/hors radiateur
for m = 1:1:max
for j = 2:1:249
d4(1,j) = d4(1,j)+(-2*d4(1,j)+d4(1,j-1)+d4(1,j+1))/(2+((n3*rho)*c)/lambdab); // west/hors radiateur
end;
end;
//d3: face superieure
d3 = ones(30,1)*278.13;
//haut
d3(1,1) = 278.13;
//haut
d3(30,1) = t2;
//haut
for m = 1:1:max
for i = 2:1:29
d3(i,1) = d1(i,1)+(-2*d1(i,1)+d1(i-1,1)+d1(i+1,1))/(2+((n1*rho)*c)/lambdab); //haut
end;
end;
d6 = ones(1,300)*0;
//d6 rassemble d4 et d5
d6(1,1:50) = d5;
d6(1,51:300) = d4;
omega(30,:) = d6;
omega(:,300) = d1;
omega(1,:) = d2;
omega(:,1) = d3;
for j = 2:1:299
for i = 2:1:29
omega(i,j) = omega(i,j)+(-4*omega(i,j)+omega(i-1,j)+omega(i,j-1)+omega(i+1,j)+omega(i,j+1))/(4+(((n1*n2)*rho)*c)/lambdab); //eq de la chaleur à deux dimensions
omegac = mtlb_s(omega,k);
end;
end;
// !! L.108: Matlab function contour not yet converted, original calling sequence used
contour(omegac); |
