Transient Heat Conduction using Explicit FDM Solution to Thermal Diffusion

Topic: Transient Heat Conduction & Thermal Diffusion

Subject: Heat Transfer

Tool: Scilab

By: Gani Comia, Dec. 2024

This article presents the solution to the transient heat conduction system represented by the thermal diffusion equation in one dimension shown in equation (1).

$$\frac{\partial u}{\partial t}\;=\;D\;\frac{\partial^2 u}{\partial x^2} \tag{1}$$

The diffusion equation can be solved with the explicit finite difference method (FDM) given by the equation (2).

$$u_{i}^{n+1}\;=\;u_{i}^{n}\;+\;\alpha\left(u_{i+1}^{n}\,-\,2 u_{i}^{n}\,+\,u_{i-1}^{n}\right) \tag{2}$$

In solving the explicit FDM, equation (3) which is defined as a stability factor, must be satisfied to maintain numerical stability:

$$\alpha\;\leq\frac{1}{2} \tag{3}$$

or:

$$\frac{D\,\Delta t}{\left(\Delta x\right)^2}\,\leq\,\frac{1}{2}\tag{4}$$

One may refer to the detailed discussion and derivation of numerical methods of its solution. Here we are going to present the application of solving a thermal diffusion problem in one dimension using FDM. The thermal diffusion equation governs transient heat conduction where the temperature within the material changes over time and space.

Consider an example material shown in Figure 1; the temperature at each finite or nodal point at a particular time as represented by \(T(x,t)\) can be solved by using the explicit FDM.

Figure 1. Sample Material Block with Initial and Boundary Temperature. 

Figure 2 is the solution for the illustrated problem. Shown is the temperature profile, \(T(x,t)\), for a particular period of time. For example, at t = 1 sec, the material had the initial temperature of 30 deg C, and as the time progressed, the temperature at each node point became equal to the boundary temperature. The thermal diffusion equation solves for the time the material temperature reaches the boundary condition from the initial state.

Figure 2. Graphical Plot of the Explicit FDM Solution.

Scilab is used to simulate the solution with the use of FDM. An explanation of each chunk of the program script is provided here for more clarity. First on the script are parameters that defined the physical properties of the problem.

// (1) parameters
L = 160;                        // mm, length of workpiece
Nx = 40;                        // unit, number of spatial points
Time = 1000;                    // sec, time consideration
D = 12;                         // mm^2/s thermal diffusivity of 4340

Derived parameters are secondary variables calculated for explicit FDM.

// (2) derived parameters
dx = L / (Nx - 1);              // mm, spatial step
alpha = 0.5;                    // stability factor
dt = alpha .* dx.^2 ./ D;       // sec, time step for stability factor

Initial and boundary temperature conditions are provided as row vectors of the size of spatial points.

// (3) initial and boundary conditions
u = zeros(Nx, 1) + 30;          // initial temperature @ 30 deg C
u(1) = 200;                     // left boundary temperature
u(Nx) = 200;                    // right boundary temperature

Temperature distributions at each time step are calculated using the for loops of the program. The plot functions are used for the visualization of the results.

// (4) temperature distribution for time steps
for n = 1:Time
    uNew = u;                   // current state of temperature
    for i = 2:Nx-1
        uNew(i) = u(i) + alpha * (u(i+1) - 2*u(i) + u(i-1));
    end
    u = uNew;                   // updated state of temperature
    
    // (5) plotting of space, time & temperature
    f = gcf()
    f.figure_size = [600,700]
    
     function plotFunc(n)
        plot(linspace(0, L, Nx), u, 'b-o');
        title("Temperature Profile @ Time = "+ string(n) + " sec")
        ylabel("Temp [ C ]")
        xgrid(1)
        a = gca()
        a.data_bounds = [0 0; L max(u)]
        a.children.children.mark_background = 9;
    endfunction
        
    select n
        case 1 then                     // @ time = 1 sec
            subplot(3,1,1)
            plotFunc(n)
        case round(Time/5) then         // @ time = 200 sec
            subplot(3,1,2)
            plotFunc(n)
        case Time then                  // @ time = 1000 sec
            subplot(3,1,3)
            plotFunc(n)
            xlabel("Length [ mm ]")
            note1 = ["https://gani-mech-toolbox.blogspot.com"]
            xstring(85,58,note1)
    end
end

Transient heat conduction is the process where the material temperature changes for a short period of time until it reaches the steady-state heat conduction, which depends on the temperature difference between boundaries.

Feel free to comment for inquiry, clarification, or suggestion for improvement. Drop your email to request the soft copy of the file.