bonjour,
ci dessous du code MATLAB qui fonctionne parfaitement:
mais bon, je dois rendre ça sous Scilab, j'utilise don la fonction "import MATLAB file" et j'obtiens ce code:
Code : Sélectionner tout - Visualiser dans une fenêtre à part
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 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*pi()/43300;%secondes rammenés en radians //a: angle que fait le soleil avec la normale à la surface if a<(pi()/2)% faces est exposé au rayonnement t2=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 t5=278.13; d5=ones(1,50)*t5;%west/radiateur elseif t<tf 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=omega-k; end end contour(omegac);
et c'est là que ça deconne, Scilab trouve des bugs dans ce truc.
Code : Sélectionner tout - Visualiser dans une fenêtre à part
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);
merci
Partager