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?

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

[AOrtenzi]

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.

Answer 1

[Animesh]

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.

[AOrtenzi]

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.

Answer 2

[surya venu]

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.

[AOrtenzi]

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.