Уравнение теплопроводности в matlabe

@leno4k@ · Регистрация: 04.11.2009

Author24 — интернет-сервис помощи студентам

Помогите разобраться. Не могу понять в чем ошибка. Переписывала код с Pascal на Matlab. Вроде результаты выдает, но не верные. Причем графически не получается вывести даже не правильные данные.
Уравнение теплопроводности самое просто pc(dT/dt)=lamda(d²T/dx²)(рс=плотность*теплоемкость,lamda- коэффициент теплопроводности) . Граничные условия тоже простые 1) t=0: T = T₀, 0<=x<=L; 2) x=0: T = T_л, t>0; 3)x=L: T = T_п, t>0. Уравнение решается простым методом прогонки.

Matlab M

clear, clc
 
N =10; % количество узлов в сетке;
t_end = 60; %окончание по времени;
L = 0.1; %длина толщены пластинки;
lamba = 46; %коэффициент теплопроводности;
ro = 7800; % плотность;
c = 460; %теплоемкость;
T0 = 20; %начальная температура;
T1 = 300; % температура на границе x=0);
Tr = 100; %температура на границе x=L ;
 
T = (1:1:500);
beta = (1:1:500);
alfa = (1:1:500);
 
h = L/(N-1); % определяем расчетный шаг сетки по пространственной координате
tau = t_end/100; %определяем расчетный шаг сетки по времени
 
for i=1:N
    T(i)= T0;
    time = 0;
    
    while time < t_end
        time = time + tau;
        alfa(1) = 0;
        beta(1) = T1;
        
        for j=2:1:N-1
            ai = lamba/(h^2);
            bi = 2*lamba/(h^2)+ro*c/tau;
            ci = lamba/(h^2);
            fi = -ro*c*T(j)/tau;
            alfa(j) = ai/(bi-ci*alfa(j-1)); % прогоночные коэффициенты
            beta(j) = (ci*beta(j-1)-fi)/(bi-ci*alfa(j-1));% прогоночные коэффициенты            
        end
    
        T(N) = Tr; %определяем значение температуры на правой границе
        
        for k = N-1:-1:1
            T(k)= alfa(k)*T(k+1)+ beta(k); 
        end
 
        subplot(111);
        plot(T(i), h*(i-1)); 
        grid on; 
    end
  
end;

Привожу код на паскале с чего переписывала

Pascal

uses crt;
const mf = 500;
type
    vector = array[1..mf] of real;
var
   i, j, N: integer;
   T, alfa, beta: vector;
   ai, bi, ci, fi: real;
   lamda, ro, c, h, tau: real;
   Tl, T0, Tr, L, t_end, time: real;
   f, g: text;
 
begin
     clrscr;
     Writeln('ВВЕДИТЕ КОЛИЧЕСтво пространственных узлов,N');
     Readln(N);
     Writeln('введите окончание по времени,t_end');
     Readln(t_end);
     Writeln('введите толщину пластины,L');
     Readln(L);
     Writeln('введитекоэффициент теплопроводности,lamda');
     Readln(lamda);
     Writeln('введите плотность ,ro');
     Readln(ro);
     Writeln('введите теплоемкость,c');
     Readln(c);
     Writeln('введите начальную температуру,T0');
     Readln(T0);
     Writeln('введите температуру на границе x=0,Tl');
     Readln(Tl);
     Writeln('введите температуру на границе x=L,Tr');
     Readln(Tr);
 
     h:= L/(N-1);
     tau:= t_end/100.0;
 
     for i:=1 to N do
        T[i]:=T0;
        time:=0;
 
        while time<t_end do {используем цикл с предусловием}
             begin
                  time:=time+tau;
                  alfa[1]:= 0.0;
                  beta[1]:= Tl;
 
                  for i:=2 to N-1 do
                  begin
                       ai:=lamda/sqr(h);
                       bi:= 2.0*lamda/sqr(h)+ro*c/tau;
                       ci:=lamda/sqr(h);
                       fi:= - ro*c*T[i]/tau;
                       alfa[i]:= ai/(bi-ci*alfa[i-1]);
                       beta[i]:= (ci*beta[i-1]-fi)/(bi-ci*alfa[i-1]);
                  end;
             T[N]:=Tr;{определяем значение температуры на правой границе}
 
             for i:=N-1 downto 1 do
                T[i]:= alfa[i]*T[i+1]+beta[i];
             end;{цикл с предусловием окончен}
 
     Assign(f, 'res.txt');
     Rewrite(f);
     Writeln(f, 'Толщина пластины L = ', L:6:4);
     Writeln(f, 'Число узлов по координате N = ', N);
     Writeln(f, 'коэффициент теплопроводности lamda = ', lamda:6:4);
     Writeln(f, 'плотность ro = ', ro:6:4);
     Writeln(f, 'теплоемкость пластины c = ', c:6:4);
     Writeln(f,'начальная температура,T0 =', T0:6:4);
     Writeln(f, 'температура на границе x=0,Tl', Tl:6:4);
     Writeln(f,'температура на границе x=L,Tr', Tr:6:4);
     Writeln(f,'результат получен с шагом по координате h=', h:6:4);
     Writeln(f,'результат получен с шагом по времени tau=', tau:6:4);
     Writeln(f,'температуратурное поле в момент времени t=', t_end:6:4);
     close(f);
 
     Assign(g, 'tempr.txt');
     Rewrite(g);
 
     for i:=1 to N do
       writeln (g,' ',h*(i-1):6:3,' ', T[i]:8:5);
     close(g);
end.

Зосима · 24.05.2013, 11:55

Лапа, у тебя в цикле ставится по одной точке, да и та затирается следующей

Попробуй перенести plot в саамый-самый конец, после всех циклов:

Matlab M

clear, clc
 
N =10; % количество узлов в сетке;
t_end = 60; %окончание по времени;
L = 0.1; %длина толщены пластинки;
lamba = 46; %коэффициент теплопроводности;
ro = 7800; % плотность;
c = 460; %теплоемкость;
T0 = 20; %начальная температура;
T1 = 300; % температура на границе x=0);
Tr = 100; %температура на границе x=L ;
 
T = 1:500;
beta = 1:500;
alfa = 1:500;
 
h = L/(N-1); % определяем расчетный шаг сетки по пространственной координате
tau = t_end/100; %определяем расчетный шаг сетки по времени
 
for i=1:N
    T(i)= T0;
    time = 0;
    
    while time < t_end
        time = time + tau;
        alfa(1) = 0;
        beta(1) = T1;
        
        for j=2:N-1
            ai = lamba/(h^2);
            bi = 2*lamba/(h^2)+ro*c/tau;
            ci = lamba/(h^2);
            fi = -ro*c*T(j)/tau;
            alfa(j) = ai/(bi-ci*alfa(j-1)); % прогоночные коэффициенты
            beta(j) = (ci*beta(j-1)-fi)/(bi-ci*alfa(j-1));% прогоночные коэффициенты
        end
        
        T(N) = Tr; %определяем значение температуры на правой границе
        
        for k = N-1:-1:1
            T(k)= alfa(k)*T(k+1)+ beta(k);
        end
    end
end
 
plot(h*(0:499),T);
grid on

Получается такой график:

@leno4k@ · 24.05.2013, 13:49 **[ТС]**

Зосима, у меня такой же график выводил но он не правильный, он должен получится примерно таким.

У меня он не получается. Скорее всего где ошибка но найти ее я не могу.

Зосима · 24.05.2013, 14:20

а ты обрати внимание на пределы по Х

возможно если верхний предел Х будет 0.1 то там как раз такая горка получится

Поидее длинна массива T должна равняться N=10, а их - 500!

Ошибка в строке 13.
Чуток подправил:

Кликните здесь для просмотра всего текста

Matlab M

clear, clc
 
N =10; % количество узлов в сетке;
t_end = 60; %окончание по времени;
L = 0.1; %длина толщены пластинки;
lamba = 46; %коэффициент теплопроводности;
ro = 7800; % плотность;
c = 460; %теплоемкость;
T0 = 20; %начальная температура;
T1 = 300; % температура на границе x=0);
Tr = 100; %температура на границе x=L ;
 
T = zeros(1,N);
alfa = zeros(1,N);
beta = zeros(1,N);
 
h = L/(N-1); % определяем расчетный шаг сетки по пространственной координате
tau = t_end/100; %определяем расчетный шаг сетки по времени
 
for i=1:N
    T(i)= T0;
    time = 0;
    
    while time < t_end
        time = time + tau;
        alfa(1) = 0;
        beta(1) = T1;
        
        for j=2:N-1
            ai = lamba/(h^2);
            bi = 2*lamba/(h^2)+ro*c/tau;
            ci = lamba/(h^2);
            fi = -ro*c*T(j)/tau;
            alfa(j) = ai/(bi-ci*alfa(j-1)); % прогоночные коэффициенты
            beta(j) = (ci*beta(j-1)-fi)/(bi-ci*alfa(j-1));% прогоночные коэффициенты
        end
        
        T(N) = Tr; %определяем значение температуры на правой границе
        
        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

Что скажешь?

можно попробовать увеличить L.
Я ставил L=0.3 и получалось вот так

@leno4k@ · 24.05.2013, 14:31 **[ТС]**

СПАСИБО ОГРОМНОЕ!!!!!!!!!!!!!! Я уже второй день голову ломаю.

Зосима · 24.05.2013, 14:37

leno4k@, тут еще одну странность увидел:
хотел увеличить число узлов N чтоб график был глаже, но ложбинка уехала!

пришлось увеличивать L до 1.3 чтоб ее увидеть...
Так же вроде не должно быть?

почему так - сказать не могу, т.к. я очень слабо понимаю как оно считает.

@Sabbat · 25.05.2013, 17:44

Смотрел когда-то такое. Первый for, для инициализации начального условия, а интегрирование идёт уже в while.

Matlab M

N       = 1000;       % количество пространственных узлов
t_end   = 60;         % окончание по времени
L       = 0.1;        % толщина пластины
lamda   = 46;         % коэффициент теплопроводности материала пластины
ro      = 7800;       % плотность материала пластины
c       = 460;        % теплоемкость материала пластины
T0      = 20;         % начальная температура
Tl      = 300;        % температура на границе, x=0
Tr      = 100;        % температура на границе, x=L
h       = L/(N-1);    % шаг сетки по пространственной координате
tau     = t_end/100;  % шаг сетки по времени
% поле температуры в начальный момент времени
for i=1:N
   T(i)=T0; 
end
% интегрирование нестационарного уравнения теплопроводности
time = 0;
while time < t_end
    time = time+tau;
    % определение начальных прогоночных коэффициентов на основе левого
    % граничного условия
    alfa(1)=0;
    beta(1)=Tl;
    % определение остальных прогоночных коэффициентов 
    for i=2:N-1
        % ai,bi,ci,fi - коэффициенти канонического представления СЛАУ с
        % трехдиагональной матрицей
        ai=lamda/h^2;
        bi=2*lamda/h^2+ro*c/tau;
        ci=lamda/h^2;
        fi=-ro*c*T(i)/tau;
        % alfa(i), beta(i) - прогоночные коэффициенты
        alfa(i)=               ai/(bi-ci*alfa(i-1));
        beta(i)=(ci*beta(i-1)-fi)/(bi-ci*alfa(i-1));
    end
    % определение значения температуры на правой границе
      T(N)=Tr;
    % определение неизвестного поля температуры
    for i=N-1:-1:1
      T(i)=alfa(i)*T(i+1)+beta(i);  
    end    
end

Зосима · 25.05.2013, 17:58

Sabbat, слух, может ты знаешь, почему при разных N расположение ложбинки меняется?

@Sabbat · 25.05.2013, 18:19

Сообщение от Зосима

Sabbat, слух, может ты знаешь, почему при разных N расположение ложбинки меняется?

Зосима, ничего не меняется....по крайней мере в моем исправленном варианте.

Сообщение от Зосима

Sabbat, слух, может ты знаешь, почему при разных N расположение ложбинки меняется?

В вашем последнем коде, перенесите end со строки 44 на 23...и всё, и ничего там не меняется. (итегрирует - while.)
Что бы понятней было, это просто один из способов решения диф. уравнения представленного неявной разностной схемой.

Зосима · 25.05.2013, 20:16

Огромнейшее мерси!

@leno4k@ · 26.05.2013, 16:56 **[ТС]**

Спасибо всем огромное кто откликнулся на мою задачу. Очень выручаете!!!!!!!!!

11.04.2014, 08:20

Зосима, leno4k@, Sabbat, ребята спасибо, за то что Вы обсуждайте такой тему тема актуальная

Еще у меня вопрос если Вам не трудна можете ответит. Например тут решено уравнения теплопроводности без химии можно ли еще добавит химии и получит график что, то у меня не получается.

@Sabbat · 11.04.2014, 23:19

Maikl2014, ну смотря что Вы представляете под понятием "химия", и чтобы решить что-то по химии должна быть какая то интерпретация того что нужно решить. Если есть наработки можете их показать.

@Maikl2020 · 13.04.2014, 18:55

Sabbat, Извиняюсь я тут новинки по этому не знаю как тут набирать математические формулы

Ну вот такой уравнения dT/dt=d^2T/dx-dT/dx+k0*tz*exp(-E/(RT)) а вот код

@topgun · 01.11.2016, 07:23

Здравствуйте. Подскажите, пожалуйста. Как дополнить начальное условие массивом из двух тел с различными температурами. Например, если разделить область на 100, то на расстоянии 50 м они соединились. У одного тела была температура 100 и другого 200. И как добавить в уравнение f(x,t): pc(dT/dt)=lamda(d2T/dx2)+f(x,t). f(x,t) имеет вид кривой.

Добавлено через 23 часа 26 минут
Разобрался с температурой, а как изменить значение F источника на прямую линию?

@apmfiit · 14.02.2019, 15:57

Sabbat, как plot написать для твоего графика?

Зосима · 15.02.2019, 14:50

apmfiit, все также как было ранее: plot(h*(0:N-1),T);
где h*(0:N-1) - значения узлов длинны х. Мы получаем распределение температуры по длинне в последний момент времени.

@leno4k@ 0 / 0 / 0 Регистрация: 04.11.2009 Сообщений: 10
	24.05.2013, 13:49 [ТС]	3
	Зосима, у меня такой же график выводил но он не правильный, он должен получится примерно таким. У меня он не получается. Скорее всего где ошибка но найти ее я не могу. 0

@leno4k@ 0 / 0 / 0 Регистрация: 04.11.2009 Сообщений: 10
	24.05.2013, 14:31 [ТС]	5
	СПАСИБО ОГРОМНОЕ!!!!!!!!!!!!!! Я уже второй день голову ломаю. 0

Зосима 5242 / 3570 / 379 Регистрация: 02.04.2012 Сообщений: 6,473 Записей в блоге: 17
	24.05.2013, 14:37	6
	leno4k@, тут еще одну странность увидел: хотел увеличить число узлов N чтоб график был глаже, но ложбинка уехала! пришлось увеличивать L до 1.3 чтоб ее увидеть... Так же вроде не должно быть? почему так - сказать не могу, т.к. я очень слабо понимаю как оно считает. 0

Зосима 5242 / 3570 / 379 Регистрация: 02.04.2012 Сообщений: 6,473 Записей в блоге: 17
	25.05.2013, 17:58	8
	Sabbat, слух, может ты знаешь, почему при разных N расположение ложбинки меняется? 0

@Sabbat 135 / 22 / 1 Регистрация: 19.10.2012 Сообщений: 42
	25.05.2013, 18:19	9
	Сообщение от Зосима Sabbat, слух, может ты знаешь, почему при разных N расположение ложбинки меняется? Зосима, ничего не меняется....по крайней мере в моем исправленном варианте. Сообщение от Зосима Sabbat, слух, может ты знаешь, почему при разных N расположение ложбинки меняется? В вашем последнем коде, перенесите end со строки 44 на 23...и всё, и ничего там не меняется. (итегрирует - while.) Что бы понятней было, это просто один из способов решения диф. уравнения представленного неявной разностной схемой. 1

Зосима 5242 / 3570 / 379 Регистрация: 02.04.2012 Сообщений: 6,473 Записей в блоге: 17
	25.05.2013, 20:16	10
	Огромнейшее мерси! 1

@leno4k@ 0 / 0 / 0 Регистрация: 04.11.2009 Сообщений: 10
	26.05.2013, 16:56 [ТС]	11
	Спасибо всем огромное кто откликнулся на мою задачу. Очень выручаете!!!!!!!!! 0

Maikl2014
	11.04.2014, 08:20	12
	Зосима, leno4k@, Sabbat, ребята спасибо, за то что Вы обсуждайте такой тему тема актуальная Еще у меня вопрос если Вам не трудна можете ответит. Например тут решено уравнения теплопроводности без химии можно ли еще добавит химии и получит график что, то у меня не получается.

@Sabbat 135 / 22 / 1 Регистрация: 19.10.2012 Сообщений: 42
	11.04.2014, 23:19	13
	Maikl2014, ну смотря что Вы представляете под понятием "химия", и чтобы решить что-то по химии должна быть какая то интерпретация того что нужно решить. Если есть наработки можете их показать. 0

@Maikl2020 Всегда онлайн 49 / 49 / 10 Регистрация: 13.04.2014 Сообщений: 1,428
	13.04.2014, 18:55	14
	Sabbat, Извиняюсь я тут новинки по этому не знаю как тут набирать математические формулы Ну вот такой уравнения dT/dt=d^2T/dx-dT/dx+k0tzexp(-E/(RT)) а вот код 0

	15.02.2019, 14:50

@topgun 0 / 0 / 0 Регистрация: 22.05.2013 Сообщений: 40
	01.11.2016, 07:23	15
	Здравствуйте. Подскажите, пожалуйста. Как дополнить начальное условие массивом из двух тел с различными температурами. Например, если разделить область на 100, то на расстоянии 50 м они соединились. У одного тела была температура 100 и другого 200. И как добавить в уравнение f(x,t): pc(dT/dt)=lamda(d2T/dx2)+f(x,t). f(x,t) имеет вид кривой. Добавлено через 23 часа 26 минут Разобрался с температурой, а как изменить значение F источника на прямую линию? 0

@apmfiit 0 / 0 / 0 Регистрация: 28.11.2018 Сообщений: 3
	14.02.2019, 15:57	16
	Sabbat, как plot написать для твоего графика? 0