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
1079
 i  Тема перемещена из форума «Помогите решить / разобраться (М)» в форум «Computer Science»

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


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

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

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



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

Сейчас этот форум просматривают: volchenok


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

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