The follow code was create to calculate the time to increase the temperature, but when I run it the result is infinite. Anyone can help me?

3 次查看（过去 30 天）
AOrtenzi
AOrtenzi 2024-6-17
评论： AOrtenzi 2024-6-18
Ran in:
Below are the tow functions used to calculate the time to increase the temp. according to the inputs: L1, L2, L3 and heat+power. The initial values were: 0.05m; 0.1 m; 0.05m and 500W respectively. But with any other values the result returns the same: inifinite seconds! I need the results in terms of plotting.
function time_to_increase_temperature = calculate_time_to_increase_temperature( ...
rho_concrete, rho_air, k_concrete, k_air, ~, cp_concrete, ...
cp_air, T_initial, T_surface, L1, L2, L3, volume_P1_P3, ~, ...
Nx, Nt, delta_T, heat_power)
% Calculate total thickness
total_thickness = L1 + L2 + L3;
% Thermal diffusivities
alpha_concrete = k_concrete / (rho_concrete * cp_concrete); % m²/s
alpha_air = k_air / (rho_air * cp_air); % m²/s
% Discretization parameters
dx = total_thickness / Nx; % Spatial step size
dt = 6 / Nt; % Time step size
% Preallocate temperature arrays
T = T_initial * ones(Nx, 1);
T_new = T;
% Define material properties at each point
len_P1 = round(Nx * L1 / total_thickness); % Length of concrete layer P1
len_P2 = round(Nx * L2 / total_thickness); % Length of air layer P2
len_P3 = Nx - len_P1 - len_P2; % Length of concrete layer P3
% Ensure len_P3 is correctly adjusted if rounding caused any issues
if len_P1 + len_P2 + len_P3 ~= Nx
len_P3 = Nx - len_P1 - len_P2;
end
% Preallocate material array
material = zeros(Nx, 1);
material(1:len_P1) = alpha_concrete;
material(len_P1+1:len_P1+len_P2) = alpha_air;
material(len_P1+len_P2+1:Nx) = alpha_concrete;
time_to_increase_temperature = inf; % Initialize to store time to increase temperature
% Calculate heat input per spatial point
heat_input_per_point = heat_power * dt / (rho_concrete * cp_concrete * volume_P1_P3 / Nx);
% Time integration using Finite Difference Method (FDM)
time = 0; % Initialize time
for n = 1:Nt
for i = 2:Nx-1
T_new(i) = T(i) + material(i) * dt / dx^2 * (T(i+1) - 2*T(i) + T(i-1));
end
% Boundary conditions
T_new(1) = T(1) + alpha_concrete * dt / dx^2 * (T(2) - T(1)) + heat_input_per_point * (T_surface - T(1));
T_new(Nx) = T(Nx-1);
% Update temperature
T = T_new;
time = time + dt;
% Check if the temperature at the right side of P3 reaches T_initial + delta_T
if T(Nx) >= T_initial + delta_T
time_to_increase_temperature = time;
break;
end
end
% Return the time to increase temperature by 1 K at the right side of P3
end
Function main_process:
function main_process()
% Parameters with fixed values
rho_concrete = 2500; % kg/m³
rho_air = 1293; % kg/m³
k_concrete = 1.35; % W/(m·K)
k_air = 0.024; % W/(m·K) assuming typical value for air
cp_concrete = 0.92 * 1000; % J/(kg·K)
cp_air = 1.005 * 1000; % J/(kg·K)
T_initial = 300; % Initial temperature in K
T_surface = 411; % Surface temperature in K
volume_P1_P3 = 0.00488; % Volume of each concrete plate in m³
Nx = 100; % Number of spatial points
Nt = 1000; % Number of time steps
delta_T = 1; % Temperature increase in K
% Prompt the user for specific input parameters
% L1 = input('Enter the thickness of the first concrete layer (m): ');
% L2 = input('Enter the thickness of the air layer (m): ');
% L3 = input('Enter the thickness of the second concrete layer (m): ');
% heat_power = input('Enter the heat power (W): ');
L1 = 1;
L2 = 1;
L3 = 1;
heat_power = 1;
% Call the function
time_to_increase_temperature = calculate_time_to_increase_temperature( ...
rho_concrete, rho_air, k_concrete, k_air, [], cp_concrete, ...
cp_air, T_initial, T_surface, L1, L2, L3, volume_P1_P3, [], ...
Nx, Nt, delta_T, heat_power);
% Display the result
fprintf('Time to increase temperature by 1 K at the right side of P3: %.2f seconds\n', time_to_increase_temperature);
% Plot the result
total_thickness = L1 + L2 + L3;
x = linspace(0, total_thickness, Nx); % Spatial points
y = linspace(0, time_to_increase_temperature, Nx); % Temporal points for plotting
plot(x, y, 'DisplayName', 'Time to Increase Temperature by 1 K');
xlabel('Position (m)');
ylabel('Time (s)');
legend('show');
title('Time to Increase Temperature by 1 K at the right side of P3');
end
main_process()
Time to increase temperature by 1 K at the right side of P3: Inf seconds
  1 个评论
AOrtenzi
AOrtenzi 2024-6-18
I have changed the code, with your suggestions, but didn't work. I will try other approach, since the temperature did not change with time, or the time remais unchanged.
Thanks.

请先登录，再进行评论。

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

回答（2 个）

Animesh
Animesh 2024-6-17
Hi @Altibano Ortenzi,
Based on my understanding, you are attempting to calculate the time required to increase the temperature by 1 K. In reviewing your code, I noticed that the variable time_to_increase_temperature is intended for this purpose. It is initially set to inf (infinity) and is later updated within an if-statement as follows:
% Check if the temperature at the right side of P3 reaches T_initial + delta_T
if T(Nx) >= T_initial + delta_T
time_to_increase_temperature = time;
break;
end
Given the inputs: L1(0.05), L2(0.1), L3(0.05), and for heat power(500), it appears that T(Nx) never exceeds the sum of T_initial(300) and delta_T(1).
Consequently, time_to_increase_temperature remains unchanged and retains its initial value of infinity. I suggest reviewing the code to ensure everything is configured correctly and to verify that the conditions for changing time_to_increase_temperature are met as expected.
  1 个评论
AOrtenzi
AOrtenzi 2024-6-17
Thanks Animesh. Actually, the temperature over the P3 surface canot be >= T(Nx), since it must be increased by the heat power source. So, I should rewrite the code to avoid this comparison.

请先登录，再进行评论。

surya venu
surya venu 2024-6-17
Hi,
To address the infinite time result in your temperature increase calculation, consider these points for correction:
  1. Heat Input Calculation: The formula for "heat_input_per_point" might not be correctly distributing the heat power across the spatial points. Simplify this calculation.
  2. Boundary Conditions: The way heat input is added in the boundary condition might not accurately model the physical process. Review and adjust this to ensure it reflects how heat is actually being applied to the system.
Here are the things you can try:
  • Revise Heat Input: Directly distribute the heat power across the spatial points without the complex division. Ensure it's added correctly at the boundary or throughout the volume as intended.
  • Adjust Boundary Condition: Ensure the boundary condition accurately models the heat addition. This could mean directly increasing the temperature based on the heat input or correctly modeling the heat flux.
Here's a simplified approach for the heat input calculation:
% Simplified heat input calculation
heat_input_per_point = heat_power / (rho_concrete * cp_concrete * volume_P1_P3);
And for the boundary condition, directly add the heat input:
T_new(1) = T(1) + (alpha_concrete * dt / dx^2 * (T(2) - T(1))) + heat_input_per_point;
Ensure these changes reflect the actual physical setup and re-evaluate the model.
Hope it helps.
  1 个评论
AOrtenzi
AOrtenzi 2024-6-17
Thanks @surya venu ! I'll try to modify the code with the lines you have to modify the calculi, but now I I'm seeing that the point described by @Animesh may be one of the code issues too. Anyway, I'll try your both suggestions.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

产品


版本

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by