Steady-State Heat Conduction Model using FDM Solution to 2D Laplace Equation

Topic: Finite Difference Method of Laplace Equation in 2D;

Subject: Heat Transfer and Numerical Methods;

Tool: Scilab;

By: Gani Comia, Dec. 2024;

The steady-state heat conduction model can be represented mathematically by Laplace equation. Shown is the Laplace equation in 2D:

$$\frac{\partial^2 T}{\partial x^2}+\frac{\partial^2 T}{\partial y^2}\;=\;0$$

The solution, \(T(x,y)\), represents the temperature distribution. The Laplace equation can be discretized using the Finite Difference Method (FDM) to approximate the solution numerically. This article will present the Scilab script of the FDM solution to the Laplace equation.

The FDM solution for a uniform grid, \((\Delta x = \Delta y)\), with the boundary temperatures specified based on Dirichlet conditions, is shown:

$$T_{i,j}\;=\;\frac{T_{i+1,j}+T_{i-1,j}+T_{i,j+1}+T_{i,j-1}}{4}$$

The final value of  \(T(x,y)\) is the result of iteration based on the convergence criterion:

$$max\left\vert T_{i,j}^{new}\;-\;T_{i,j}^{old}\right\vert\;<\;\epsilon$$

Consider a sample 2D domain of dimensions 30x30 units with the following constant temperature conditions. 

Fig. 1. Sample 2D Domain with Constant Temperature.

Shown below is the temperature profile as represented by \(T(x,y)\) using the FDM solution to the Laplace equation.

Fig. 2. Temperature Profile \(T(x,y)\).

The Scilab script as a toolbox can be used to implement the FDM numerical solution for the steady-state heat conduction model represented mathematically using the Laplace equation.

Scilab Script

// Copyright (C) 2024 - Gani Comia
// Date of creation: 8 Dec 2024
// FDM Solution to 2D Laplace Eq, d2T/dx2 + d2T/dy2 = 0
clear;clc;
// domain and grid parameters
Lx = 30;                    // units, domain length at x-axis
Ly = 30;                    // units, domain length at y-axis
dx = 1;                     // grid spacing at x-axis
dy = 1;                     // grid spacing at y-axis
nx = Lx./dx;                // no. of grid points in x-direction
ny = Ly./dy;                // no. of grid points in y-direction
max_iter = 1000;            // max no. of iterations
tolerance = 1e-6;           // convergence tol.

// temperature field initial value
T = zeros(nx, ny);
// boundary value conditions
    T(:,1) = 0;             // T = 0 C @ y = 1
T(1,:) = 0;                 // T = 0 C @ x = 1
T(:,ny) = 0;                // T = 0 deg C @ y = 30
T(nx,:) = 100;              // T = 100 deg C @ x = 30

// finite difference method (Jacobi iteration)
for iter = 1:max_iter
    Tnew = T;
    for i = 2:nx-1
        for j = 2:ny-1
            Tnew(i,j)=0.25*(T(i+1,j)+T(i-1,j)+T(i,j+1)+T(i,j-1));
        end
    end
    // convergence check
    if norm(Tnew - T) < tolerance then
        break;
    end
    T = Tnew;
end
// calculation status 
mprintf("\n\tCalculation completed!\n")

// plot properties
clf;
f = gcf()
f.figure_size = [700,650]
ax = gca()
ax.tight_limits = ["on" "on" "off"];

// plotting results
surf(T);
title(["$Temperature\;Profile$","$T(x,y)$"],"fontsize",4);
xlabel("x-axis");
ylabel("y-axis");
zlabel("$T$");
colormap(jet);
colorbar(0,100);

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