Форум программистов, компьютерный форум, киберфорум
Matlab
Форум программистов и сисадминов Киберфорум > Научный форум > Математические программы > Matlab
Войти
Регистрация
Восстановить пароль
Карта форума Темы раздела Блоги Сообщество Поиск Заказать работу  
 
Рейтинг 4.60/55: Рейтинг темы: голосов - 55, средняя оценка - 4.60
177 / 143 / 50
Регистрация: 07.02.2014
Сообщений: 489
1

Метод прогонки для двумерного уравнения теплопроводности (задать температуру в центре пластины)

24.02.2017, 18:11. Показов 10074. Ответов 11
Метки нет (Все метки)
Author24 — интернет-сервис помощи студентам
Добрый день, форумчане!

Решаю уравнение двумерное теплопроводности. Для случая, когда задаю ГУ в виде температур и нулевых тепловых потоков -- все хорошо.

Но вопрос, как мне задать ячейку в центре пластины с определенной температурой (это для начала, потом в виде ГУ 3-го рода)?

Привожу свой код:
Кликните здесь для просмотра всего текста
Matlab M
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

clc, clear all, close all
tic
Lx = 10; dx = 0.10; Nx = Lx/dx;
Ly = 10; dy = 0.10; Ny = Ly/dy;
lamda = 384.0; rho = 8800; C = 381; at = lamda/(rho*C);
T0 = 5; T_left = 50; T_right = -10; T_up = 75; T_down = -10;
t_end = 3600*24; tau = 1800;
h1 = 50; T_ext1 = 60; h2 = 50; T_ext2 = 60; q = 1e07;
q1 = 5000;
 
T = ones(Nx,Ny)*T0;
 
a = zeros(1,Nx-1);
b = zeros(1,Nx-1);
c = zeros(1,Nx-1);
f = zeros(1,Nx-1);
alfa = zeros(1,Nx-1);
beta = zeros(1,Nx-1);
 
time = 0;
 
while time < t_end
    
    time = time + tau;
    
    time1 = tau/2;
    
    for k = 1:Ny
        alfa(1) = 0;
        beta(1) = T_right;
        for j = 2:Nx-1
            a(j) = lamda/dx^2;
            b(j) = 2.0*lamda/dx^2+rho*C/time1;
            c(j) = lamda/dx^2;
            f(j) = -rho*C*T(j,k)/time1;
            
            alfa(j) = a(j)/(b(j)-c(j)*alfa(j-1));
            beta(j) = (c(j)*beta(j-1)-f(j))/(b(j)-c(j)*alfa(j-1));
        end
        
        T(Nx,k) = T_left;
        
        for j = Nx-1:-1:1
            T(j,k) = alfa(j)*T(j+1,k)+beta(j);
        end
    end
    
    time2 = tau/2;
    
    for j = 1:Nx
        alfa(1) = 1;
        beta(1) = 0*dy/lamda;
        for k = 2:Ny-1
            a(k) = lamda/dy^2;
            b(k) = 2.0*lamda/dy^2+rho*C/time2;
            c(k) = lamda/dy^2;
            f(k) = -rho*C*T(j,k)/time2;
            
            alfa(k) = a(k)/(b(k)-c(k)*alfa(k-1));
            beta(k) = (c(k)*beta(k-1)-f(k))/(b(k)-c(k)*alfa(k-1));
        end
        T(j,Ny) = T(j,Ny-1);
        for k = Ny-1:-1:1
            T(j,k) = alfa(k)*T(j,k+1) + beta(k);
        end
    end
end
toc
x = linspace(0,Lx,Nx);
y = linspace(0,Ly,Ny);
[X,Y] = meshgrid(x,y);
surf(X,Y,T) 
axis equal tight
caxis([0 50])
set(gca,'ydir','normal')
colorbar, grid on % закрепляем сетку и цвет.шкалу
colormap('jet')
hold on
view(90,-90)


В данном куске реализован метод расщепления по пространственной координате и каждая из подсхем решается методом прогонки. При запуске кода, получается вполне ожидаемая картина:
Метод прогонки для двумерного уравнения теплопроводности (задать температуру в центре пластины)


Теперь я хочу задать температуру в центральном узле --
Matlab M
1

T(Nx/2,Ny/2) = 80
То получаю бред:
Метод прогонки для двумерного уравнения теплопроводности (задать температуру в центре пластины)


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

Matlab M
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

clc, clear all, close all
tic
Lx = 10; dx = 0.10; Nx = Lx/dx;
Ly = 10; dy = 0.10; Ny = Ly/dy;
lamda = 384.0; rho = 8800; C = 381; at = lamda/(rho*C);
T0 = 5; T_left = 50; T_right = -10; T_up = 75; T_down = -10;
t_end = 3600*24; tau = 1800;
h1 = 50; T_ext1 = 60; h2 = 50; T_ext2 = 60; q = 1e07;
q1 = 5000;
 
T = ones(Nx,Ny)*T0;
 
a = zeros(1,Nx-1);
b = zeros(1,Nx-1);
c = zeros(1,Nx-1);
f = zeros(1,Nx-1);
alfa = zeros(1,Nx-1);
beta = zeros(1,Nx-1);
 
time = 0;
 
while time < t_end
    
    time = time + tau;
    
    time1 = tau/2;
    
    T(Nx/2,Ny/2) = 80;
    
    for k = 1:Ny
        alfa(1) = 0;
        beta(1) = T_right;
        for j = 2:Nx-1
            a(j) = lamda/dx^2;
            b(j) = 2.0*lamda/dx^2+rho*C/time1;
            c(j) = lamda/dx^2;
            f(j) = -rho*C*T(j,k)/time1;
            
            alfa(j) = a(j)/(b(j)-c(j)*alfa(j-1));
            beta(j) = (c(j)*beta(j-1)-f(j))/(b(j)-c(j)*alfa(j-1));
        end
        
        T(Nx,k) = T_left;
        
        for j = Nx-1:-1:1
            T(j,k) = alfa(j)*T(j+1,k)+beta(j);
        end
    end
    
    time2 = tau/2;
    
    for j = 1:Nx
        alfa(1) = 1;
        beta(1) = 0*dy/lamda;
        for k = 2:Ny-1
            a(k) = lamda/dy^2;
            b(k) = 2.0*lamda/dy^2+rho*C/time2;
            c(k) = lamda/dy^2;
            f(k) = -rho*C*T(j,k)/time2;
            
            alfa(k) = a(k)/(b(k)-c(k)*alfa(k-1));
            beta(k) = (c(k)*beta(k-1)-f(k))/(b(k)-c(k)*alfa(k-1));
        end
        T(j,Ny) = T(j,Ny-1);
        for k = Ny-1:-1:1
            T(j,k) = alfa(k)*T(j,k+1) + beta(k);
        end
    end
end
toc
x = linspace(0,Lx,Nx);
y = linspace(0,Ly,Ny);
[X,Y] = meshgrid(x,y);
surf(X,Y,T) 
axis equal tight
caxis([0 50])
set(gca,'ydir','normal')
colorbar, grid on % закрепляем сетку и цвет.шкалу
colormap('jet')
hold on
view(90,-90)


Как видно из картинки, что-то в узле T(Nx/2,Ny/2) меняется, но совсем не так, как должно

Поэтому у меня такой вопрос: Как все-таки правильно задать температуру в центре пластины? В каком месте? Может, ее как-то нужно запихнуть в коэффициенты прогонки?

Жду Ваших ответов! Очень нужна помощь
1
Лучшие ответы (1)
Сообщение: #10151181
Programming
Эксперт
94731 / 64177 / 26122
Регистрация: 12.04.2006
Сообщений: 116,782
 24.02.2017, 18:11
Ответы с готовыми решениями:

Решение уравнения теплопроводности по явной схеме (метод прогонки)
Возник вопрос при реализации метода прогонки при решении уравнения теплопроводности по неявной...

с++ метод прогонки, уравнение теплопроводности
∂T/∂t=a (∂^2 T)/(∂x^2 ) , где а - коэффициент температуропроводности (м^2/с) Необходимо решить...

Неявная схема решения уравнения теплопроводности методом прогонки
Всем доброго времени суток! Есть построение уравнения Лапласа методом прогонки. Есть явная схема...

Смешанная задача для уравнения теплопроводности метод разделения переменных
http://s2.ipicture.ru/uploads/20130617/h76YT5a0.jpg

11
6830 / 4890 / 2065
Регистрация: 02.02.2014
Сообщений: 13,048
24.02.2017, 18:56 2
Цитата Сообщение от Norwall Посмотреть сообщение
хочу задать температуру в центральном узле --
задайте ее после цикла, т.е. после сформированного массива
Изображения
 
0
177 / 143 / 50
Регистрация: 07.02.2014
Сообщений: 489
25.02.2017, 13:14  [ТС] 3
Krasme, спасибо за ответ. Я так тоже пробовал, но, как видно, температура ячейки не оказывает никакого влияния на температуры соседних ячеек не зависимо от времени расчета.
0
6830 / 4890 / 2065
Регистрация: 02.02.2014
Сообщений: 13,048
25.02.2017, 13:30 4
мне не непонятно, что вы хотите получить, а пока пальцем в небо...
Название: 2017-02-25_132930.png Просмотров: 270 Размер: 1.3 Кб
0
177 / 143 / 50
Регистрация: 07.02.2014
Сообщений: 489
25.02.2017, 13:45  [ТС] 5
Krasme, хочу получить следующее:
Метод прогонки для двумерного уравнения теплопроводности (задать температуру в центре пластины)


Как видно из рисунка, центральная ячейка контактирует с соседними и в них происходит изменение температуры посредствам теплопередачи. Данная картинка получена путем решения уравнения теплопроводности по явной схеме.
0
6830 / 4890 / 2065
Регистрация: 02.02.2014
Сообщений: 13,048
25.02.2017, 13:53 6
Лучший ответ Сообщение было отмечено Norwall как решение

Решение

может, одна точка слишком мала...
Цитата Сообщение от Norwall Посмотреть сообщение
немного измененный код:
заменить 28-ую строку на
Matlab M
1

T(Nx/2-2:Nx/2+2,Ny/2-2:Ny/2+2) = 80;
у меня получилась эта картинка
Метод прогонки для двумерного уравнения теплопроводности (задать температуру в центре пластины)
1
6830 / 4890 / 2065
Регистрация: 02.02.2014
Сообщений: 13,048
25.02.2017, 13:53 7
вот и картинки стали похожи
0
177 / 143 / 50
Регистрация: 07.02.2014
Сообщений: 489
25.02.2017, 14:16  [ТС] 8
Krasme, мде... забавно) Спасибо за помощь)

Добавлено через 10 минут
Насколько я понял, нельзя обращаться к дробной ячейке. Если записать так
Matlab M
1

T(round(Nx/2),round(Ny/2)) = 80;
то все нормально считает:
0
6830 / 4890 / 2065
Регистрация: 02.02.2014
Сообщений: 13,048
25.02.2017, 14:27 9
Nx=100
Nx/2=50
откуда дробь?

Добавлено через 1 минуту
и потом, если дробь бы оказалась в индексации массива, то матлаб заругался бы..
0
177 / 143 / 50
Регистрация: 07.02.2014
Сообщений: 489
25.02.2017, 15:15  [ТС] 10
Krasme, блин, это я уже совсем попутал)))

[удалено]
0
0 / 0 / 0
Регистрация: 05.12.2017
Сообщений: 9
24.12.2019, 21:34 11
Norwall, Алексей,не могли бы вы помочь с одним вопросом. У вас очень хороший код,но не подскажите как добавить в него разрывность на пластину? У меня есть код разрыва для одномерного случая+стержня,код я приведу ниже,но вот не могу объединить их.
Matlab M
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

clear all;
clc;
 
a=3; 
 
 h=0.1; %шаг по пространству
x=0:h:a; N=length(x); 
Nx = N;% по x количество шагов
t_end = 10; %время
Lx = 1; %размер по оси x
lamda = 100.0;%теплопроводность
rho = 7200;%плотность на участке
C = 460;%теплоемкость
T0 = 220;%температура пластины в момент времени T0
Tl = 0; %темп верхнего края (граничные условия)
Tr = 0;%(граничные условия) низ
at = lamda/(rho*C);%коэф уравнения 
hx = Lx/(Nx-1); %шаг по x
 
tau = 0.25*(hx^2)/at;%шаг по времени, вычисляемый для устойчивости уравнения
 
tau_t(Nx)=0;
 for i=1:Nx
       if i< ceil(Nx/4)  rho(i)=7200; 
       elseif i> ceil(2*Nx/3) rho(i)=40;
       else rho(i)=400;
       end;
 end;
 
for i = 1:Nx
          at(i) = lamda/(rho(i)*C);
         tau_t(i) = 0.25*(hx^2)/at(i);
end;
 
tau=min(tau_t); 
 
T=ones(Nx,1)*T0; %если начальная температура стержня - равномерная
%for i=1:Nx        %Строим начальное неравномерное распределение температуры
 %   if x(i)<=0.5*a 
 %       T(i)=((2*T0)/a)*x(i); T(i)=T0;
 %   end 
  %  if x(i)>0.5*a 
  %      T(i)=((2*T0)/a)*(a-x(i)); T(i)=T0;
  %  end 
%end 
 
T(1) = 0;   % и Т1 на одной грани, 
T(end) = 0; % Т2 - на противоположной грани
 
time = 0;
 for i=1:Nx
   k_x(i)= (lamda/(rho(i)*C))*tau/hx^2; 
 end;
 
%рисуем начальный профиль температуры 
%толстой красной линией 
plot(x,T,'Color','red','LineWidth',3); 
hold on 
 
 
 while time < t_end
    
    time = time + tau;
 
   TT=T;
  
for i = 2:Nx-1
        %T(i,1) = TT(i,1) + (lamda/(rho*C))*tau/hx^2*(TT(i+1,1)-2*TT(i,1)+TT(i-1,1)) + ...
         %   (lamda/(rho*C))*tau/hy^2*(TT(i,2)-TT(i,1));
       % T(i,1) = TT(i,1) + k_x*(TT(i+1,1)-2*TT(i,1)+TT(i-1,1)) + ...
          %  k_y*(TT(i,2)-TT(i,1));
        T(i) = TT(i) + k_x(i+1)*TT(i+1)-k_x(i)*2*TT(i)+k_x(i-1)*TT(i-1);
       % T(Nx) = k_x(Nx)*(-TT(Nx))+k_x(Nx)*TT(Nx);
    end
   
    %рисуем вычисленный профиль температуры 
    plot(x,T); title(num2str(time));
    hold on
    pause(0.2);
   
end 
 
%рисуем конечный профиль температуры 
%толстой зеленой линией 
plot(x,T,'Color','green','LineWidth',3);
0
177 / 143 / 50
Регистрация: 07.02.2014
Сообщений: 489
25.12.2019, 06:16  [ТС] 12
m3nka, добрый день!

Под разрывностью вы имеете разность теплофизических свойств? Или фазовые переходы?

Если первое, то задайте свойства в виде двумерного массива, тогда вы сможете задать локально любые свойства. Главное не забыть, что в цикле по времени нужно совершать операции с матрицами свойств. Ниже пример цикла нахождения энтальпии в 3D случае:
Matlab M
1
2
3
4
5
6
7
8
9
10
11
12

for i = 2:Nx-1
    for j = 2:Ny-1
        for k = 2:Nz-1
            H(i,j,k) = H(i,j,k) + tau/hx*((-1)*(T(i,j,k)-T(i-1,j,k))/(hx/2*(1/lamda(i-1,j,k)+1/lamda(i,j,k)))-...
                (-1)*(T(i+1,j,k)-T(i,j,k))/(hx/2*(1/lamda(i,j,k)+1/lamda(i+1,j,k)))) + ...
                tau/hy*((-1)*(T(i,j,k)-T(i,j-1,k))/(hy/2*(1/lamda(i,j-1,k)+1/lamda(i,j,k)))-...
                (-1)*(T(i,j+1,k)-T(i,j,k))/(hy/2*(1/lamda(i,j,k)+1/lamda(i,j+1,k)))) + ...
                tau/hz*((-1)*(T(i,j,k)-T(i,j,k-1))/(hz/2*(1/lamda(i,j,k-1)+1/lamda(i,j,k)))-...
                (-1)*(T(i,j,k+1)-T(i,j,k))/(hz/2*(1/lamda(i,j,k)+1/lamda(i,j,k+1))));
        end
    end
end
1
25.12.2019, 06:16
IT_Exp
Эксперт
87844 / 49110 / 22898
Регистрация: 17.06.2006
Сообщений: 92,604
 25.12.2019, 06:16
Помогаю со студенческими работами здесь

Метод конечных элементов для уравнения теплопроводности. Как найти производную (градиент)?
Здравствуйте. Помогите, пожалуйста, разобраться. Решаю нестационарную двумерную задачу...

Решение уравнения теплопроводности. Вывод двумерного массива
Проблема следующая. Написал алгоритм решения уравнения теплопроводности в виде явной схемы:...

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

Разностное решение уравнения теплопроводности.Явный метод
Здравствуйте. Помогите пожалуйста. Нужно немного доделать и ,если не сложно, подписать что в какой...


Искать еще темы с ответами

Или воспользуйтесь поиском по форуму:
12
Ответ Создать тему
КиберФорум - форум программистов, компьютерный форум, программирование
Поддержать - Реклама - Обратная связь
Powered by vBulletin
Copyright ©2000 - 2024, CyberForum.ru