2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки




Начать новую тему Ответить на тему
 
 Явная схема трехмерной теплопроводности
Сообщение01.04.2017, 17:20 


30/05/16
3
Здравствуйте. Написал явную схему трехмерной теплопроводности самым простым методом сеток, с ГУ третьего рода.. Но решение находит неправильно, а именно: вывожу температурное поле в сечениях по х, у, z. И нигде нет совпадения температуры. Вот код в MATLAB:
Код:
clear all; clc
% Matetial
rho = 7800;
c = 500;
lambda = 46;

a = lambda / (rho * c);       % thermal diffusitivy

% Parameters
A = 1.5;           % i          %  height
B = 0.25;          % j           %  width 
C = 0.5 ;          % k          %  length
t = 6;             % p


%  Bordering
T0 = 20;

hx = 0.01;
hy = 0.01;
hz = 0.1;
tau = 1 / (a * 30 * (1 / (hx * hx) + 1 / (hy * hy)) + 1 / (hz * hz)) ;
%hy = (A / B) * hx;

Nx = round(A / (hx));
Ny = round(B / (hy));
Nz = round(C / (hz));
Nt = round(t / tau);

%Tempearture of environment
Tl = 1300; % left
Tr = 1300; % right
Tu = 1300; % up
Td = 1300; % down
Tf = 1300; % front
Tb = 1300; % back

% flow
alpha_u = 1900;
alpha_l = 850;
alpha_r = 850;
alpha_d = 1900/5;
alpha_f = 850;
alpha_b = 850;

%
tic
for i = 1:Nx+1
    for j = 1:Ny+1
        for k = 1:Nz+1
            T(i, j, k, 1) = T0;
        end
    end
end


for p = 1:Nt-1
   
    for i = 1:1:Nx
        for j = 1:1:Ny           
            for k =1:1:Nz
               
                % left side
                if i == 1
                    T(i, j, k, p+1) = T(i, j, k, p) + (alpha_l * (Tl - T(i, j, k, p))...
                        + (lambda * (T(i+1, j, k, p) - T(i, j, k, p)) / hx)) * ...
               (2 * tau / (c * rho * hx));
           
             
                 % right side
                else if i == A
                       T(i, j, k, p+1) = T(i, j, k, p) + (alpha_r * (Tr - T(i, j, k, p))...
                        + (lambda * (T(i-1, j, k, p) - T(i, j, k, p)) / hx)) * ...
               (2 * tau / (c * rho * hx));
           
                 
                 % down side
                    else if j == 1
                      T(i, j, k, p+1) = T(i, j, k, p) + (alpha_d * (Td - T(i, j, k, p))...
                        +(lambda * (T(i, j+1, k, p) - T(i, j, k, p)) / hx)) * ...
               (2 * tau / (c * rho * hy));
                       
           
                   % up side
                        else  if j == B
                    T(i, j, k, p+1) = T(i, j, k, p) + (alpha_u * (Tu - T(i, j, k, p))...
                        +(lambda * (T(i, j-1, k, p) - T(i, j, k, p)) / hx)) * ...
               (2 * tau / (c * rho * hy));
           
           
                   % front side
                            else  if  k == 1
                    T(i, j, k, p+1) = T(i, j, k, p) + (alpha_f * (Tf - T(i, j, k, p))...
                        +(lambda * (T(i, j, k+1, p) - T(i, j, k, p)) / hx)) * ...
               (2 * tau / (c * rho * hz));
           
           
                    % back side
                                else  if  k == C
                    T(i, j, k, p+1) = T(i, j, k, p) + (alpha_b * (Tb - T(i, j, k, p))...
                        +(lambda * (T(i, j, k-1, p) - T(i, j, k, p)) / hz)) * ...
               (2 * tau / (c * rho * hz));
           
                       
                    % implement points
                                    else   
                           T(i , j, k, p+1)  = T(i, j, k, p) +...
                           a * tau * (T(i+1, j, k, p) - 2 * T(i, j, k, p) + T(i-1, j, k, p))/hx^2 +...                               
                           a * tau * (T(i, j+1, k, p) - 2 * T(i, j, k, p) + T(i, j-1, k, p))/hy^2 +...                             
                           a * tau * (T(i, j, k+1, p) - 2 * T(i, j, k, p) + T(i, j, k-1, p))/hz^2;
                                       
                                    end
                               
                                        end
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
   
 

% section in a yz plane
x = round(A / 2);

% section in a xz plane
y = 125;

% section in a xy plane
z = 3;


for i = 1:Nx+1
     for j = 1:Ny
         for k = 1:Nz
     Txyend(i, j) = T(i, j, z, Nt);
     Tyzend(j, k) = T(x, j, k, Nt);
     Txzend(i, k) = T(i, 1, k, Nt);
         end
     end
end
   

    pcolor(Txyend)
  colorbar
  grid on
  s1 = ['section with z = ', z];
  title(s1)
  xlabel('x')
  ylabel('y')
 
figure;
    pcolor(Tyzend)
  colorbar
  grid on
  s2 = ['section with x = ', x];
  title(s2)
  xlabel('y')
  ylabel('z')
 
  figure;
    pcolor(Txzend)
  colorbar
  grid on
  s3 = ['section with y = ', y];
  title(s3)
  xlabel('x')
  ylabel('z')
   

  toc
 

Кто может заметил, или знает, что не так, подскажите, пожалуйста.

 Профиль  
                  
 
 Posted automatically
Сообщение01.04.2017, 17:24 
Модератор


19/10/15
1196
 i  Тема перемещена из форума «Помогите решить / разобраться (М)» в форум «Computer Science»

 Профиль  
                  
 
 Re: Явная схема трехмерной теплопроводности
Сообщение01.04.2017, 17:48 
Заслуженный участник


12/07/07
4438
1. trofimovep, запишите исходную задачу (формулы с описаниями).
2. Обычно задают число узлов сетки (вдоль направления), и уже потом по этому числу находят длину шага (вдоль направления).
3. i, j, k изменяются в цикле от 1 до Nx, Ny, Nz. Проверка if i == A вызывает удивление.
4. Я бы не вставлял в основной цикл вычисление значений в граничных точках, а выполнил бы это отдельными циклами. (Для скорости вычислений.)
5. Я бы задавал размер массива при помощи zeros . Это позволит избежать изменения размера массива в цикле.

 Профиль  
                  
Показать сообщения за:  Поле сортировки  
Начать новую тему Ответить на тему  [ Сообщений: 3 ] 

Модераторы: Karan, Toucan, PAV, maxal, Супермодераторы



Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей


Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете добавлять вложения

Найти:
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group