progMatlab.ml

tic

clear; clf; clc; hold off;

%*****************************************

% Outside air characteristics *

%*****************************************

Ta=305.8; Cpa=1006; roa=1.293;hE=10;

%*****************************************

% Characteristics of the PVC casing *

%*****************************************

Ct1=0.2; ro1=1380; Cp1=1046; T01=Ta;

%*****************************************

% Characteristics of the insulating PSE *

%*****************************************

Ct2=0.047;Cppe=1300;rope=12;T02=Ta;

%*****************************************

% Specifications indoor air *

%*****************************************

CpI=Cpa; roI=roa;hI=200;T03=Ta;

I=0.3;% 

Ctal=204;roal=2707; Cpal=896; T04=Ta;% 

Ctce=26.35; Cpce=900; roce=3350; T05=Ta;% 

L=1.5; alpha=200*1e-6; Ct=1.5; sig=10;r=1.82;R=3.9; T0=Ta;


Ct6=Ctce; Cp6=900; ro6=roce; T06=Ta;


Ct7=Ctal;ro7=roal; Cpa7=Cpal; T07=300.6;

%*****************************************

%        Cooling water      *

%*****************************************

qmee=0.000556;Ve=0.027; roe=1000; Cpe=4*Cpa; he=991.4;Tee=300.6;

%*****************************************

% Characteristics of the tank water PVC *          

%*****************************************

Ct9=Ct1; ro9=ro1; Cp9=1046; T010=Tee;

%*****************************************

dx=5; N=10;

%******************************************************************

%Expression of resistance thermal and capacity thermal *          

%******************************************************************

R01=1/hE;R12=dx/Ct2; R23=1/hI; R34=dx/Ctal;R45=dx/Ctce;R56=dx/Ct;

R67=R45; R78=R34; R89=1/he;R910=R12; R10a=R01;

R01;

R23;

R34;

R45;

R56;

R67;

R78;

R89;

R910;

R10a;

%******************************************************************

C1=ro1*dx*Cp1; C2=C1; C3=roal*dx*Cpal; C4=roce*dx*Cpce; C5=C4; C6=C5;

C7=C3; C8=C7; C9=C1; C10=C9;

C1;

C2;

C3;

C4;

C5;

C6;

C7;

C8;

C9;

C10;

%********************************

% Construction of matrix B *

%********************************

a1=1/C1; a11=-((1/R01)+(1/R12))*a1;

a21=(1/R12)*a1; a12=(1/R12)*a1;

a22=-((1/R12)+(1/R23))*a1;

a2=1/C3;a32=(1/R23)*a2;

a23=(1/R23)*a1;

a33=-((1/R23)+(1/R34))*a2;

a3=1/C4;a43=(1/R34)*a3;

a34=(1/R34)*a2;

a44=-((1/R34)+(1/R45))*a3;

a54=(1/R45)*a3;

a45=(1/R45)*a3;

a55=((1/R)+(1/R45)+alpha*I)*a3;

a65=(1/R)*a3;

a56=-(1/R)*a3;

a66=(alpha*I-(1/R)-(1/R67))*a3;

a76=(1/R67)*a2;

a67=(1/R67)*a3;

a77=-((1/R67)+(1/R78))*a2;

a87=(1/R78)*a2;

a78=(1/R78)*a2;

a88=-((1/R78)+(1/R89))*a2;

a98=(1/R89)*a1;

a89=(1/R89)*a2;

a99=-((1/R89)+(1/R910))*a1;

a109=(1/R910)*a1;

a910=(1/R910)*a1;

a1010=-((1/R910)+(1/R10a))*a1;

%******************************************************************

T(1)=Ta; T(2:N-3)=Ta;T(8:N)=260;

M=diag([a11, a22, a33, a44, a55, a66, a77, a88, a99, a1010]);

Mup=diag([a12, a23, a34, a45, a56, a67, a78, a89, a910],1);

Mdown=diag([a21, a32, a43, a54, a65, a76, a87,a98, a109],-1);

B=M+Mup+Mdown;

%******************************************************************

% Definition of the unicolumn matrix D

%******************************************************************

D=[(1/R01)*a1*Ta, 0, 0, 0, -(1/2)*(r/C5)*I^2,(1/2)*(r/C6)*I^2 , 0, 0, 0, (1/R01)*a1*Ta]';

%*******************************************

% Resolution using Runge Kutta 4

%********************************************

M=B;

M;

S=D;

S;

T=T';

T;

n=0; t=0; h=0.01; m=0;

axis([0 N 260 Ta]);

j=[1,2:length(T)+2];

T_p=[Ta,T',Ta+10];

T_p;

plot(j,T_p);

text(j(2),T_p(2),['t=',int2str(t),'s']);

for k=1:10        %Number of iterations

   for m=1:10*N

       n=n+1;

       k1=h*(M*T+S);

       k2=h*(M*T+S+k1/2);

       k3=h*(M*T+S+k2/2);

       k4=h*(M*T+S+k3);

       T=T+(k1+2*k2+2*k3+k4)/6;

       t=h*n;

   end

   hold on;

   j=[0,2:length(T)+2];

   T_p=[Ta,T',Ta+10];

   plot(j,T_p);

   grid on;

   xlabel('time : per s');

   ylabel('T (in °K)');

   pause(0.1)

end

T_p;

t_mis=toc