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

@topgun · Регистрация: 22.05.2013

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

Всем привет. Пытался перевести задачу с кода Pascal на Matlab, но выдает ошибки. Подскажите , пожалуйста в чем ошибка?
До этого обсуждалась эта тема. Различие только в добавлении нелинейности в коэффициенте lamda. Соответственно и прогоночные коэффициенты изменились.

Исходный код Pascal:

Pascal

uses crt;
const mf=500;
eps=1e-5;
type
vector=array[1..mf] of real;
var {раздел описания переменных, которые мы будем использовать в
программе}
i, j, N, s : integer;
T, Ts, Tn, alfa, beta : vector;
ai, bi, ci, fi, max1, max2 : real;
ro, c, h, tau : real;
Th, T0, Tc, L, t_end, time : real;
f, g, f1, f2 : text;
function lamda(x:real):real;
{функция вычисления коэффициента теплопроводности по формуле
(31)}
begin
lamda:=5500/(560+x)+0.942*(1e-10)*x*sqr(x);
end;
begin
clrscr;
{с клавиатуры вводим все необходимые входные параметры}
Writeln('Введите количество пространственных узлов, N');
Readln(N);
Writeln('Введите окончание по времени, t_end');
Readln(t_end);
Writeln('Введите толщину пластины, L');
Readln(L);
Writeln('Введите плотность материала пластины, ro');
Readln(ro);
Writeln('Введите теплоемкость материала пластины, c');
Readln(c);
Writeln('Введите начальную температуру в К, T0');
Readln(T0);
Writeln('Введите температуру в К на границе х=0, Th');
Readln(Th);
Writeln('Введите температуру в К на границе х=L, Tc');
Readln(Tc);
{определяем расчетный шаг сетки по пространственной координате}
h:=L/(N-1);
{определяем расчетный шаг сетки по времени}
tau:=t_end/100.0;
{определяем поле температуры в начальный момент времени}
for i:= 2 to N-1 do
T[i]:=T0;
{заводим файл содержащий количество итераций на каждом шаге по
времени}
Assign(f1,'iter.txt');
Rewrite(f1);
{проводим интегрирование нестационарного уравнения
теплопроводности}
time:=0;
while time<t_end do {используем цикл с предусловием}
begin
{увеличиваем переменную времени на шаг τ}
time:=time+tau;
{запоминаем поле температуры на предыдущем временном слое}
for i:= 1 to N do
Tn[i]:=T[i];
s:=0; {будем считать число итераций}
{цикл с постусловием, позволяющий итерационно вычислять поле
температуры}
repeat
inc(s);
{запоминаем поле температуры на предыдущей итерации}
for i:= 1 to N do
Ts[i]:=T[i];
{определяем начальные прогоночные коэффициенты на основе левого
граничного условия}
alfa[1]:=0.0;
beta[1]:=Th;
{цикл с параметром для определения прогоночных коэффициентов по
формуле (8)}
for i:= 2 to N-1 do
begin
{ai, bi, ci, fi – коэффициенты канонического представления системы
уравнений с трехдиагональной матрицей}
ai:=0.5*(lamda(T[i])+lamda(T[i+1]))/sqr(h);
ci:=0.5*(lamda(T[i-1])+lamda(T[i]))/sqr(h);
bi:=ai+ci+ro*c/tau;
fi:=-ro*c*Tn[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]:=Tc;
{используя соотношение (7) определяем неизвестное поле
температуры}
for i:= N-1 downto 1 do
T[i]:=alfa[i]*T[i+1]+beta[i];
{определяем максимум модуля разности температур на данной и
предыдущей итерации}
max1:=abs(T[1]-Ts[1]);
for i:= 2 to N do
if max1 < abs(T[i]-Ts[i]) then max1:=abs(T[i]-Ts[i]);
{определяем максимальное по модулю значение температуры на данной
итерации}
max2:=abs(T[1]);
for i:= 2 to N do
if max2 < abs(T[i]) then max2:=abs(T[i]);
until max1/max2<=eps; {цикл с постусловием окончен}
Writeln(f1,'В момент времени ',time:6:4,' проведено ',s,' итераций');
end; {цикл с предусловием окончен}
{выводим результат в файл}
Assign(f,'res.txt');
Rewrite(f);
Writeln(f,'Толщина пластины L = ',L:6:4);
Writeln(f,'Число узлов по координате N = ',N);
Writeln(f,'Плотность материала пластины ro = ',ro:6:4);
Writeln(f,'Теплоемкость материала пластины с = ',c:6:4);
Writeln(f,'Начальная температура T0 = ',T0:6:4);
Writeln(f,'Температура на границе x = 0, Th = ',Th:6:4);
Writeln(f,'Температура на границе x = L, Tc = ',Tc: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]-273:8:5);
close(g);
close(f1);
end.

Код на Matlab:

Matlab M

clc; close all; clear all;
 
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;  % шаг сетки по времени
p1      = 0.01;       % коэффициент в lamda
eps     = 1e-5;       % точность 
a0      = 0.5;          % коэффициент в lamda
 
function
    lamda   = a0*(1-p1*T(i));
end
% поле температуры в начальный момент времени
for i=2:N-1
   T(i)=T0; 
end
 
% интегрирование нестационарного уравнения теплопроводности
time = 0;
while time < t_end
    time = time+tau;
    
    %запоминаем поле температуры на предыдущем временном слое
    for i=1:N
        Tn(i)=T(i)
        s=0 %будем считать число итераций
        %repeat 
       % inc(s)
        for i=1:N
            Ts(i)=T(i)
    
    % определение начальных прогоночных коэффициентов на основе левого
    % граничного условия
 
    alfa(1)=0;
    beta(1)=Tl;
    % определение остальных прогоночных коэффициентов 
    for i=2:N-1
        % ai,bi,ci,fi - коэффициенти канонического представления СЛАУ с
        % трехдиагональной матрицей
        ai=0.5*(lamda(T(i))+lamda(T(i+1)))/sqr(h);
        ci=0.5*(lamda(T(i-1))+lamda(T(i)))/sqr(h);
        bi=ai+ci+r0*c/tau;
        
        fi=-ro*c*Tn(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);  
      
     % определяем максимальное по модулю значение температуры на данной
     % итерации
     max1=abs(T(1)-Ts(1));
     for i=2:N
         if max1<abs(T(i)-Ts(i)) then max1=abs(T(i)-Ts(i));
             max2=abs(T(1));
             for i=2:N
                 if max2<abs(T(i)) then max2=abs(T(i));
                     until max1/max2<=eps;
    
        end
    end
end
 
plot(h*(0:N-1),T);
grid on

@Centurio · 02.11.2016, 14:13

topgun, лучше бы вы формулы привели, чем код на Паскале.

@topgun · 04.11.2016, 10:01 **[ТС]**

$https://www.cyberforum.ru/cgi-bin/latex.cgi?\rho c\frac{dT}{dt}=\lambda \frac{{d}^{2}T}{{dx}^{2}}$

@Centurio · 07.11.2016, 17:44

topgun, вот один вариант, с использованием прогоночных коэффициентов:

Matlab M

 N       = 1000;       % количество пространственных узлов
t_end   = 60;         % окончание по времени
L       = 0.1;        % толщина пластины
ro      = 7800;       % плотность материала пластины
c       = 460;        % теплоемкость материала пластины
T0      = 20;         % начальная температура
Tl      = 300;        % температура на границе, x=0
Tr      = 100;        % температура на границе, x=L
h       = L/(N-1);    % шаг сетки по пространственной координате
tau     = t_end/100;  % шаг сетки по времени
p1      = 0.01;       % коэффициент в lamda
a0      = 0.5;          % коэффициент в lamda
t=t_end+1;
T=repmat(T0,t,N);
T(:,end)=Tr;
alfa=zeros(1,N);
bet=repmat(Tl,1,N);
k=ro*c/tau;
for I=2:t
    for J=2:N-1
        lamda=a0*(1-p1*T(I-1,J));
        A=lamda/h^2;
        B=2*A+k;
        C=A;
        F=-k*T(I-1,J);
        s=B-C*alfa(J-1);
        alfa(J)=A/s;
        bet(J)=(C*bet(J-1)-F)/s;
    end
    for J=N-1:-1:1
        T(I,J)=alfa(J)*T(I,J+1)+bet(J);
    end
end
x=linspace(0,0.1,N);
plot(x,T(end,:)),grid on

Вариант с использованием трёхдиагональной матрицы:

Matlab M

N       = 1000;       % количество пространственных узлов
t_end   = 60;         % окончание по времени
L       = 0.1;        % толщина пластины
ro      = 7800;       % плотность материала пластины
c       = 460;        % теплоемкость материала пластины
T0      = 20;         % начальная температура
Tl      = 300;        % температура на границе, x=0
Tr      = 100;        % температура на границе, x=L
h       = L/(N-1);    % шаг сетки по пространственной координате
tau     = t_end/100;  % шаг сетки по времени
p1      = 0.01;       % коэффициент в lamda
a0      = 0.5;          % коэффициент в lamda
t=t_end+1;
k=ro*c/tau;
T=repmat(T0,N-2,t);
for I=2:t
lamda=a0*(1-p1*T(:,I-1));
  A=lamda/h^2;
  B=2*A+k;
  C=A;
  M1=diag(C(2:end),-1);
  M2=diag(-B);
  M3=diag(A(1:end-1),1);
  M=M1+M2+M3;
  F=-k*T(:,I-1);
  F(1)=F(1)-A(1)*Tl;
  F(end)=F(end)-C(end)*Tr;
  T(:,I)=M\F;
end
T=[repmat(Tl,1,t);T;repmat(Tr,1,t)];
x=linspace(0,0.1,N);
plot(x',T(:,end)),grid on

@topgun · 09.11.2016, 19:51 **[ТС]**

Спасибо. Но не смогли бы вы подсказать как изменить код. Это код для того же уравнения, отличие в начальном условие, когда вначале соединяются два тела (for i=1:N/2 T(i)=400 end
for i=N/2:N T(i)=300 end ) и добавляется слагаемое f, которое представляет собой кривую (fi=(0.1+i*h*2)*ro*c;
). Т.е уравнение имеет вид: $https://www.cyberforum.ru/cgi-bin/latex.cgi?\rho c\frac{dT}{dt}=\lambda \frac{{d}^{2}T}{d{x}^{2}}+f$
Сам код: Неявная схема

Matlab M

clc; close all; clear all; 
 
N = 1000; % количество пространственных узлов 
t_end = 60; % окончание по времени 
L = 0.5; % толщина пластины 
lamda = 60; % коэффициент теплопроводности материала пластины 
ro = 7800; % плотность материала пластины 
c = 460; % теплоемкость материала пластины 
Tl = 300; % температура на границе, x=0 
w = 0.1; 
h = L/(N-1); % шаг сетки по пространственной координате 
tau = t_end/100; % шаг сетки по времени 
% поле температуры в начальный момент времени 
 
for i=1:N/2 
T(i)=400 
 
end 
for i=N/2:N 
T(i)=300 
end 
% интегрирование нестационарного уравнения теплопроводности 
time = 0; 
while time < t_end 
time = time+tau; 
%Условие на границе
Tr = 400*(1+0.1*(cos(w*time))); 
% определение начальных прогоночных коэффициентов на основе левого 
% граничного условия 
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; 
%f в виде кривой
fi=(0.1+i*h*2)*ro*c; 
 
 
 
% alfa(i), beta(i) - прогоночные коэффициенты 
alfa(i)= ai/(bi-ci*alfa(i-1)); 
beta(i)= (((T(i)*ro*c/tau))+ci*beta(i-1)+fi)/(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 
 
plot(h*(0:N-1),T); 
grid on

Как изменить сейчас, если $https://www.cyberforum.ru/cgi-bin/latex.cgi?\lambda =a0*(1-a1*T)$ . Т.е. уравнение принимает нелинейный вид.

@Centurio · 09.11.2016, 20:36

topgun, вы даже не попытались разобраться в моём коде, а опять просто сделали кальку с кода в Паскале с кучей лишних команд и операций.

@topgun · 09.11.2016, 21:14 **[ТС]**

Centurio, я попытался изменить ваш код, но у меня не совпал график при задании lamda постоянным. Код, который я прикрепил, полностью рабочий. Как я и писал выше, изменил начальное условие и граничное, а также функцию f в виде кривой.

@Centurio · 10.11.2016, 06:48

Matlab M

N = 1000; % количество пространственных узлов 
t_end = 60; % окончание по времени 
L = 0.5; % толщина пластины 
lamda = 60; % коэффициент теплопроводности материала пластины 
ro = 7800; % плотность материала пластины 
c = 460; % теплоемкость материала пластины 
Tl = 300; % температура на границе, x=0 
w = 0.1; 
h = L/(N-1); % шаг сетки по пространственной координате 
tau = t_end/100; % шаг сетки по времени 
 
n=fix(N/2);
T=[repmat(400,1,n),repmat(300,1,N-n)];% вектор начальных температур
Tr=400*(1+0.1*(cos(w*(0:tau:t_end))));% вектор температур правой границы
x=0:h:L;% координаты узлов сетки
% вспомогательные коэффициенты:
k=ro*c/tau;
A=lamda/h^2;
B=2*A+k;
C=A;
f=(0.1+2*x(2:end))*ro*c;
%-----------------------
alfa=zeros(1,N);% подготовка массива для коэффициентов alfa
bet=repmat(Tl,1,N);% подготовка массива для коэффициентов beta
for I=2:length(Tr)
    % Вычисление коэффициентов alfa и beta:
    for J=2:N-1
        s=B-C*alfa(J-1);
        alfa(J)=A/s;
        bet(J)=(T(J)*k+C*bet(J-1)+f(J))/s;
    end
    % Определение неизвестного поля температуры:
    T(N)=Tr(I);
    for J=N-1:-1:1
        T(J)=alfa(J)*T(J+1)+bet(J);
    end
end
plot(x,T),grid on

@topgun · 10.11.2016, 19:49 **[ТС]**

Centurio, спасибо, я хотел уточнить, а что нужно изменить, чтобы вместо постоянного значения lamda указать зависимость: lamda=a0*(1-a1*Ti); Пытался вставить в коде, который вы скинули в последнем сообщение, он выдаёт ошибку. Вы до этого скидывали код, с учётом этой зависимости, но я не смог внести свои условия на температуры. А в последнем коде все работает, только не получилось в него добавить зависимость lamda. Вы бы не смогли подсказать, пожалуйста.

@Centurio · 10.11.2016, 20:42

Нужно вычисление lamda внести в цикл J в виде lambda=a0*(1-a1*T(J-1)). Также в цикл внести вычисление коэффициентов A, B, C.

Зосима · 11.11.2016, 15:08

Matlab M

lamda = @(T) a0*(1-p1*T); % создали функцию lambda(T)

@Centurio Модератор 1700 / 1552 / 520 Регистрация: 13.09.2015 Сообщений: 5,370
	02.11.2016, 14:13	2
	topgun, лучше бы вы формулы привели, чем код на Паскале. 0

@topgun 0 / 0 / 0 Регистрация: 22.05.2013 Сообщений: 40
	04.11.2016, 10:01 [ТС]	3
	$https://www.cyberforum.ru/cgi-bin/latex.cgi?\rho c\frac{dT}{dt}=\lambda \frac{{d}^{2}T}{{dx}^{2}}$ 0

@Centurio Модератор 1700 / 1552 / 520 Регистрация: 13.09.2015 Сообщений: 5,370
	09.11.2016, 20:36	6
	topgun, вы даже не попытались разобраться в моём коде, а опять просто сделали кальку с кода в Паскале с кучей лишних команд и операций. 0

@topgun 0 / 0 / 0 Регистрация: 22.05.2013 Сообщений: 40
	09.11.2016, 21:14 [ТС]	7
	Centurio, я попытался изменить ваш код, но у меня не совпал график при задании lamda постоянным. Код, который я прикрепил, полностью рабочий. Как я и писал выше, изменил начальное условие и граничное, а также функцию f в виде кривой. 0

@topgun 0 / 0 / 0 Регистрация: 22.05.2013 Сообщений: 40
	10.11.2016, 19:49 [ТС]	9
	Centurio, спасибо, я хотел уточнить, а что нужно изменить, чтобы вместо постоянного значения lamda указать зависимость: lamda=a0(1-a1Ti); Пытался вставить в коде, который вы скинули в последнем сообщение, он выдаёт ошибку. Вы до этого скидывали код, с учётом этой зависимости, но я не смог внести свои условия на температуры. А в последнем коде все работает, только не получилось в него добавить зависимость lamda. Вы бы не смогли подсказать, пожалуйста. 0

@Centurio Модератор 1700 / 1552 / 520 Регистрация: 13.09.2015 Сообщений: 5,370
	10.11.2016, 20:42	10
	Нужно вычисление lamda внести в цикл J в виде lambda=a0(1-a1T(J-1)). Также в цикл внести вычисление коэффициентов A, B, C. 0

Зосима 5242 / 3570 / 379 Регистрация: 02.04.2012 Сообщений: 6,473 Записей в блоге: 17
	11.11.2016, 15:08	11
	Matlab M lamda = @(T) a0(1-p1T); % создали функцию lambda(T) 0