Errors in transient heat transfer node

48 次查看（过去 30 天）
준 최
준 최 2024-12-2，10:13
% Simulation Parameters
dt =10; % Time step [s]
t_max = 2000; % Maximum simulation time [s]
time_steps = t_max / dt;
Nx = 34; Ny = 40; % Grid dimensions
dx = 1; dy = 1; % Uniform grid spacing
h = ones(Ny, Nx) * 15;
h_air=15;
h_water=300;
h_in=2;
% Initialize material map and properties
material_map = 10 * ones(Ny, Nx); % Initialize material map
% ABS
material_map(1:13, 6:8) = 1;
material_map(1:13, 28:30) = 1;
material_map(1:4, 9:27) = 1;
% Handle
material_map(29:30, 6:9) = 2;
material_map(20:28, 6:7) = 2;
% Vacuum
material_map(10:31,11) = 3; % Row 11, Columns 10 to 31
material_map(10, 11:12) = 3; % Row 10, Columns 11 to 12
material_map(6:10, 12) = 3; % Rows 6 to 10, Column 12
material_map(6, 12:24) = 3; % Row 6, Columns 12 to 24
material_map(6:10, 24) = 3; % Rows 6 to 10, Column 24
material_map(10, 24:25) = 3; % Row 10, Columns 24 to 25
material_map(10:31,25) = 3;
% Water
material_map(12:30,13) = 4; % Row 13, Columns 12 to 30
material_map(8:30, 14) = 4; % Rows 8 to 30, Column 14
material_map(8:26, 15:21) = 4; % Rows 8 to 26, Columns 15 to 21
material_map(8:30, 22) = 4; % Rows 8 to 30, Column 22
material_map(12:30, 23) = 4; % Rows 12 to 30, Column 23
% Ice
material_map(27:30,15:21 ) = 5; % Rows 15 to 21, Columns 27 to 30
% Internal Air
material_map(31:34, 15:21) = 6; % Rows 31 to 34, Columns 15 to 21
material_map(31:32, 14) = 6; % Rows 31 to 32, Column 14
material_map(31, 13) = 6; % Row 31, Column 13
material_map(31, 23) = 6; % Row 31, Column 23
material_map(31:32, 22) = 6; % Rows 31 to 32, Column 22
% Lid
material_map(35:36, 14:22) = 7; % Rows 35 to 36, Columns 14 to 22
% Stainless Steel (Inside)
material_map(33, 13:14) = 8; % Row 33, Columns 13 to 14
material_map(33, 22:23) = 8; % Row 33, Columns 22 to 23
material_map(32, 23:24) = 8; % Row 32, Columns 23 to 24
material_map(32, 12:13) = 8; % Row 32, Columns 12 to 13
material_map(11:31, 12) = 8; % Rows 11 to 31, Column 12
material_map(11:31, 24) = 8; % Rows 11 to 31, Column 24
material_map(11, 23:24) = 8; % Row 11, Columns 23 to 24
material_map(11, 12:13) = 8; % Row 11, Columns 12 to 13
material_map(7:10, 13) = 8; % Rows 7 to 10, Column 13
material_map(7:10, 23) = 8; % Rows 7 to 10, Column 23
material_map(7, 13:23) = 8; % Row 7, Columns 13 to 23
% Stainless Steel (Outside) 위치 정의
material_map(34, 12:14) = 9; % Row 34, Columns 12 to 14
material_map(34, 22:24) = 9; % Row 34, Columns 22 to 24
material_map(33, 11:12) = 9; % Row 33, Columns 11 to 12
material_map(33, 24:25) = 9; % Row 33, Columns 24 to 25
material_map(32, 25:26) = 9; % Row 32, Columns 25 to 26
material_map(32, 10:11) = 9; % Row 32, Columns 10 to 11
material_map(9:32, 26) = 9; % Rows 9 to 32, Column 26
material_map(31:32, 11) = 9; % Rows 31 to 32, Column 11
material_map(9:32, 10) = 9; % Rows 9 to 32, Column 10
material_map(5:9, 11) = 9; % Rows 5 to 9, Column 11
material_map(5:9, 25) = 9; % Rows 5 to 9, Column 25
material_map(5, 11:25) = 9; % Row 5, Columns 11 to 25
k(material_map == 1) = 0.23; % ABS
k(material_map == 2) = 0.3; % Handle
k(material_map == 3) = 0.0003; % Vacuum
k(material_map == 4) = 0.575; % Water
k(material_map == 5) = 2.2; % Ice
k(material_map == 6) = 0.025; % Internal Air
k(material_map == 7) = 0.2; % Lid
k(material_map == 8) = 16.2; % Stainless Steel (Inside)
k(material_map == 9) = 15.2; % Stainless Steel (Outside)
k(material_map == 10) = 0.025; % External Air
rho(material_map == 1) = 1040;
rho(material_map == 2) = 950;
rho(material_map == 3) = 0; % Vacuum
rho(material_map == 4) = 1000;
rho(material_map == 5) = 917;
rho(material_map == 6) = 1.2;
rho(material_map == 7) = 950;
rho(material_map == 8) = 8000;
rho(material_map == 9) = 8000;
rho(material_map == 10) = 1.2;
c(material_map == 1) = 1.5;
c(material_map == 2) = 2.0;
c(material_map == 3) = 0; % Vacuum
c(material_map == 4) = 4.186;
c(material_map == 5) = 2.06;
c(material_map == 6) = 1.005;
c(material_map == 7) = 1.75;
c(material_map == 8) = 0.5;
c(material_map == 9) = 0.5;
c(material_map == 10) = 1.005;
T(material_map == 1) = 21;
T(material_map == 2) = 19;
T(material_map == 3) = 14; % Vacuum
T(material_map == 4) = 4;
T(material_map == 5) = -15;
T(material_map == 6) = 8;
T(material_map == 7) = 15;
T(material_map == 8) = 5;
T(material_map == 9) = 19;
T(material_map == 10) = 20; % External
% Initialize material properties
[k, rho, c, T] = initialize_material_properties(material_map, Nx, Ny);
% Constants
T_infinite = 20; % Ambient temperature
sigma = 5.67e-8; % Stefan-Boltzmann constant
epsilon = 0.9; % Emissivity
video = VideoWriter('Temperature_Simulation.avi'); % 동영상 파일 이름
video.FrameRate = 10; % 프레임 속도 설정
open(video); % 동영상 파일 열기
% Simulation loop
all_temperatures = zeros(Ny, Nx, t_max / dt); % Ny x Nx 크기에 TimeSteps 추가
for t = 1:t_max
% 새로운 온도 저장 배열
T_new = T;
% 내부 노드 및 경계 처리
for i = 2:Ny-1
for j = 2:Nx-1
if rho(i, j) > 0 && c(i, j) > 0
% 물질 내부의 전도 (Fourier's Law)
Txx = (T(i+1, j) - 2*T(i, j) + T(i-1, j)) / dx^2;
Tyy = (T(i, j+1) - 2*T(i, j) + T(i, j-1)) / dy^2;
alpha = k(i, j) / (rho(i, j) * c(i, j)); % 열확산계수
T_cond_internal = alpha * (Txx + Tyy);
% 물질 간 경계에서의 열 전달
T_boundary = 0;
T_new(i, j) = T(i, j) + dt * T_cond_internal;
end
% 1. ABS - External Air (대류)
if material_map(i, j) == 1 && is_material_boundary(i, j, material_map, 10)
T_conv = h(i, j) * (T_infinite - T(i, j)) / (rho(i, j) * c(i, j));
T_new(i, j) = T(i, j) + dt * T_conv;
end
% 2. ABS - External Stainless Steel (전도)
if material_map(i, j) == 1 && is_material_boundary(i, j, material_map, 9)
k_eff = 2 * k(i, j) * k(i+1, j) / (k(i, j) + k(i+1, j)); % Effective conductivity
T_cond = k_eff * (T(i+1, j) - T(i, j)) / dx^2;
T_new(i, j) = T(i, j) + dt * T_cond / (rho(i, j) * c(i, j));
end
% 3. External Stainless Steel - External Air (대류)
if material_map(i, j) == 9 && is_material_boundary(i, j, material_map, 10)
T_conv = h(i, j) * (T_infinite - T(i, j)) / (rho(i, j) * c(i, j));
T_new(i, j) = T(i, j) + dt * T_conv;
end
% 4. External Stainless Steel - Vacuum (전도)
if material_map(i, j) == 9 && is_material_boundary(i, j, material_map, 3)
k_eff = k(i+1, j); % Vacuum's conductivity
T_cond = k_eff * (T(i+1, j) - T(i, j)) / dx^2;
T_new(i, j) = T(i, j) + dt * T_cond / (rho(i, j) * c(i, j));
end
% 5. Internal Stainless Steel - Water (대류)
if material_map(i, j) == 8 && is_material_boundary(i, j, material_map, 4)
T_conv = h_water * (T(i+1, j) - T(i, j)) / (rho(i, j) * c(i, j));
T_new(i, j) = T(i, j) + dt * T_conv;
end
% 6. Internal Stainless Steel - Vacuum (전도)
if material_map(i, j) == 8 && is_material_boundary(i, j, material_map, 3)
k_eff = k(i+1, j); % Vacuum's conductivity
T_cond = k_eff * (T(i+1, j) - T(i, j)) / dx^2;
T_new(i, j) = T(i, j) + dt * T_cond / (rho(i, j) * c(i, j));
end
% 7. Water - Internal Air (대류)
if material_map(i, j) == 4 && is_material_boundary(i, j, material_map, 6)
T_conv = 50 * (T(i+1, j) - T(i, j)) / (rho(i, j) * c(i, j));
T_new(i, j) = T(i, j) + dt * T_conv;
end
% 8. Water - Ice (전도)
if material_map(i, j) == 4 && is_material_boundary(i, j, material_map, 5)
k_eff = 2 * k(i, j) * k(i+1, j) / (k(i, j) + k(i+1, j));
T_cond = k_eff * (T(i+1, j) - T(i, j)) / dx^2;
T_new(i, j) = T(i, j) + dt * T_cond / (rho(i, j) * c(i, j));
end
% 9. Ice - Internal Air (대류)
if material_map(i, j) == 5 && is_material_boundary(i, j, material_map, 6)
T_conv = h_in * (T(i+1, j) - T(i, j)) / (rho(i, j) * c(i, j));
T_new(i, j) = T(i, j) + dt * T_conv;
end
% 10. Internal Air - Lid (대류)
if material_map(i, j) == 6 && is_material_boundary(i, j, material_map, 7)
T_conv = h_in* (T(i+1, j) - T(i, j)) / (rho(i, j) * c(i, j));
T_new(i, j) = T(i, j) + dt * T_conv;
end
% 11. Lid - External Stainless Steel (전도)
if material_map(i, j) == 7 && is_material_boundary(i, j, material_map, 9)
k_eff = 2 * k(i, j) * k(i+1, j) / (k(i, j) + k(i+1, j));
T_cond = k_eff * (T(i+1, j) - T(i, j)) / dx^2;
T_new(i, j) = T(i, j) + dt * T_cond / (rho(i, j) * c(i, j));
end
% 12. Lid - External Air (대류)
if material_map(i, j) == 7 && is_material_boundary(i, j, material_map, 10)
T_conv = h(i, j) * (T_infinite - T(i, j)) / (rho(i, j) * c(i, j));
T_new(i, j) = T(i, j) + dt * T_conv;
end
% 13. Handle - External Stainless Steel (전도)
if material_map(i, j) == 2 && is_material_boundary(i, j, material_map, 9)
k_eff = 2 * k(i, j) * k(i+1, j) / (k(i, j) + k(i+1, j));
T_cond = k_eff * (T(i+1, j) - T(i, j)) / dx^2;
T_new(i, j) = T(i, j) + dt * T_cond / (rho(i, j) * c(i, j));
end
% 14. Handle - External Air (대류)
if material_map(i, j) == 2 && is_material_boundary(i, j, material_map, 10)
T_conv = h(i, j) * (T_infinite - T(i, j)) / (rho(i, j) * c(i, j));
T_new(i, j) = T(i, j) + dt * T_conv;
end
if material_map(i, j) == 4 % 물
% 상하좌우 인접 노드와의 대류 효과
T_up = T(i-1, j);
T_down = T(i+1, j);
T_left = T(i, j-1);
T_right = T(i, j+1);
if is_material_boundary(i, j, material_map, 6)
if material_map(i, j) == 4 % 물과 공기
T_conv = h_water * (T(i+1, j) - T(i, j)) / (rho(i, j) * c(i, j));
else % 기본 전도
T_cond = k_eff * (T(i+1, j) - T(i, j)) / dx^2;
end
T_new(i, j) = T(i, j) + dt * (T_cond + T_conv);
end
% 대류에 의한 열 이동 추가
T_conv_internal = h_water/ (rho(i, j) * c(i, j)) * ...
((T_up - T(i, j)) + (T_down - T(i, j)) + ...
(T_left - T(i, j)) + (T_right - T(i, j)));
T_new(i, j) = T(i, j) + dt * T_conv_internal;
end
% 내부 전도와 경계에서의 열 전달을 모두 반영하여 온도 계산
T_new(i, j) = T(i, j) + dt * (T_cond_internal + T_boundary);
end
end
% 추가 경계 조건 적용
% 하단 경계 (대류 경계, T_infinite = 20°C)
for j = 1:Nx
T_conv_bottom = h_air * (T_infinite - T(Ny, j)) / (rho(Ny, j) * c(Ny, j));
T_new(Ny, j) = T(Ny, j) + dt * T_conv_bottom;
end
% 왼쪽 경계 (대류 경계, T_infinite = 20°C)
for i = 1:Ny
T_conv_left = h_air * (T_infinite - T(i, 1)) / (rho(i, 1) * c(i, 1));
T_new(i, 1) = T(i, 1) + dt * T_conv_left;
end
% 상단 경계 (단열 경계)
T_new(1, :) = T_new(2, :);
% 우측 경계 (고정 온도)
T_new(:, Nx) = 400;
% 온도 업데이트
T = T_new;
if mod(t, 10) == 0 % High-resolution colormap
imagesc(flipud(T), [-20, 40]); % Show the initial temperature field
colorbar;
title(['Time = ', num2str(t*dt), ' seconds']);
xlabel('X-Axis (Nodes)');
ylabel('Y-Axis (Nodes)');
drawnow; % Ensure the figure updates properly
frame = getframe(gcf); % Capture the current figure
writeVideo(video, frame); % Save the frame to the video
if mod(t, 500) == 0
writematrix(T, ['Node_Temperatures_Time_', num2str(t * dt), '.csv']);
disp(['Saved node temperatures at Time = ', num2str(t * dt), ' seconds']);
end
end
if t * dt >= t_max
break; % 종료
end
end
writematrix(T, 'Final_Node_Temperatures.csv');
disp('Final node temperatures saved to Final_Node_Temperatures.csv');
close(video);
% Local Function: is_material_boundary
function is_boundary = is_material_boundary(i, j, material_map, target_material)
% 기본적으로 경계 아님
is_boundary = false;
% 인접 노드를 검사 (상하좌우)
if i > 1 && material_map(i-1, j) == target_material % 위쪽
is_boundary = true;
elseif i < size(material_map, 1) && material_map(i+1, j) == target_material % 아래쪽
is_boundary = true;
elseif j > 1 && material_map(i, j-1) == target_material % 왼쪽
is_boundary = true;
elseif j < size(material_map, 2) && material_map(i, j+1) == target_material % 오른쪽
is_boundary = true;
end
end
% Local Function: initialize_material_properties
function [k, rho, c, T] = initialize_material_properties(material_map, Nx, Ny)
% Initialize properties
k = zeros(Ny, Nx); rho = zeros(Ny, Nx); c = zeros(Ny, Nx);
T = zeros(Ny, Nx);
% ABS
k(material_map == 1) = 0.23;
rho(material_map == 1) = 1040;
c(material_map == 1) = 1.5;
T(material_map == 1) = 21;
% Handle
k(material_map == 2) = 0.3;
rho(material_map == 2) = 950;
c(material_map == 2) = 2.0;
T(material_map == 2) = 19;
% Vacuum
k(material_map == 3) = 0.0003;
rho(material_map == 3) = 0;
c(material_map == 3) = 0;
T(material_map == 3) = 14;
% Water
k(material_map == 4) = 0.575;
rho(material_map == 4) = 1000;
c(material_map == 4) = 4.186;
T(material_map == 4) = 4;
% Ice
k(material_map == 5) = 2.2;
rho(material_map == 5) = 917;
c(material_map == 5) = 2.06;
T(material_map == 5) = -15;
% Internal Air
k(material_map == 6) = 0.025;
rho(material_map == 6) = 1.2;
c(material_map == 6) = 1.005;
T(material_map == 6) = 8;
% Lid
k(material_map == 7) = 0.2;
rho(material_map == 7) = 950;
c(material_map == 7) = 1.75;
T(material_map == 7) = 15;
% Stainless Steel (Inside)
k(material_map == 8) = 16.2;
rho(material_map == 8) = 8000;
c(material_map == 8) = 0.5;
T(material_map == 8) = 5;
% Stainless Steel (Outside)
k(material_map == 9) = 15.2;
rho(material_map == 9) = 8000;
c(material_map == 9) = 0.5;
T(material_map == 9) = 19;
% External Air
k(material_map == 10) = 0.025;
rho(material_map == 10) = 1.2;
c(material_map == 10) = 1.005;
T(material_map == 10) = 20;
end
I have two questions. firstone is that if i change the dx,and dy to 0.5, the figure flashes with check patterns.
also i wonder why the water near inside stainless steel goes below zero. please help

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Introduction to Installation and Licensing 的更多信息

标签

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by