N=10;
M=60;
alfa1=0.5;
alfa2=0.5;
lamda1=0.084;
lamda2=4;
cp=1000;
c2=660;
ro1=0.3;
v0=1;
ro2=1500;
Tgor=1480;
T0=315;
Q=15130000;
k0=500000000;
E=126000;
R=8.34;
eta0=0.077;
kap=(alfa1*lamda1+alfa2*lamda2)/(ro1*cp+ro2*c2);
L=0.1;
tz=(L*L)/kap;
h=L/(N-1);
a1=(tz*kap)/(L*L);
a2=(tz*v0)/L;
alfa = zeros(1,N);
beta = zeros(1,N);
h = L/(N-1); % определяем расчетный шаг сетки по пространственной координате
tau =M/100; %определяем расчетный шаг сетки по времени
for i=1:N
T(i)= T0;
time = 0;
while time < M
time = time + tau;
alfa(1) = 0;
beta(1) = T0;
for j=2:N-1
ai = a2/h-a1/h/h;
bi = (-2*a1)/h/h-1/tau+a2/h;
ci = -a1/h;
fi = -T(j)/tau-tz*k0*(1-T(j))*exp(-(E)/(R*(T0+T(j)*(Tgor-T0))));
alfa(j) = ai/(bi-ci*alfa(j-1)); % прогоночные коэффициенты
beta(j) = (ci*beta(j-1)-fi)/(bi-ci*alfa(j-1));% прогоночные коэффициенты
end
T(N) = T0; %определяем значение температуры на правой границе
for k = N-1:-1:1
T(k)= alfa(k)*T(k+1)+ beta(k);
end
end
end
plot(h*(0:N-1),T);
grid on