Simulation of 2D heat conduction in steady and unsteady forms

Simulation of a 2D Heat Conduction problem in steady and unsteady/transient forms using iterative methods.

Project Objectives:

1. Solving the 2 Dimensional Heat conduction equation in the generalized form using various iterative techniques:

i. Explicit Solver (for unsteady state)

ii. Implicit Solvers (for both steady and unsteady state)

(a). Gauss Jacobi method

(b). Gauss Siedel method

(c). Successive Over Relaxation method

2. Analysis of steady state solution and unsteady state solution and comparison of them:

3. Stability Analysis of the unsteady state solution: ______________________________________________________________________

[Part: 1/3]: SOLVING THE 2D HEAT CONDUCTION EQUATION:

Introduction: Thermal conduction is the transfer of heat by molecular collisions of particles and movement of electrons with high kinetic energy within a body. It is an important mode of heat transfer that occurs between two bodies in contact. Conduction requires an temperature gradient between two bodies to take place, the other important modes of heat transfer are convection and radiation.

Conduction is essentially a time-dependant phenomena, the transient nature arises from the governing energy equation in 3 dimensions and transient conduction eventually moves towards steady state. (i.e. when the temperature gradients that have been established within a body do not change with time)

Problem specifications:

i. Assume that the domain is an unit square

ii. Assume Nx = Ny (Number of gridpoints is same in both x and y directions or the mesh is a square mesh)

iii. Boundary Conditions: (Drichlet type)

1. Left wall = 400K

2. Top wall = 600K

3. Right wall = 800K

4. Bottom wall = 900K

The Plate is subjected to the boundary conditions at time t = 0. The Temperature of the plate is considered to be the ambient temperature i.e. T = 298K

TRANSIENT / UNSTEADY STATE SOLUTION:

Derivation of the Numerical Scheme:

From the generalized Energy equation in 3 Dimensions we make the following assumptions:

i. No convection occurs

ii. No internal heat generation or dissipation

iii. Constant properties of system i.e. `alpha != f(t,x,y)`

Hence for a 2 Dimensional condition we have:

`(del^2T)/(delx^2)+(del^2T)/(dely^2) = alpha(delT)/(delt)`

We Discretize the spatial terms by central differencing (`O(h^2)`) i.e.

`(del^2T)/(delx^2) = (T_(i-1,j)-2T_(i,j)+T_(i+1,j))/(Deltax^2)`

`(del^2T)/(dely^2) = (T_(i,j-1)-2T_(i,j)+T_(i,j+1))/(Deltay^2)`

We discretize the temporal terms by forward differencing (`O(h)`) i.e.

`(delT)/(delt) = ((T^(n+1))_(i,j)-(T^(n))_(i,j))/(Deltat)`

Hence we have:

`(T_(i-1,j)-2T_(i,j)+T_(i+1,j))^(n+1)/(Deltax^2) + (T_(i,j-1)-2T_(i,j)+T_(i,j+1))^(n+1)/(Deltay^2) = ((T^(n+1))_(i,j)-(T^(n))_(i,j))/(Deltat)`

We define two constants `K1` and `K2` such that:

`K1 = alpha(Deltat)/((Deltax)^2)`

`K2 = alpha(Deltat)/((Deltay)^2)`

Here `Deltax = Deltay`

Hence the equation:

`(T^(n+1))_(i,j)=(T^(n))_(i,j)+K1((T^(n+1))_(i-1,j)-2(T^(n+1))_(i,j)+(T^(n+1))_(i+1,j))+K2((T^(n+1))_(i,j-1)-2(T^(n+1))_(i,j)+(T^(n+1))_(i,j+1))`

`(T^(n+1))_(i,j)=((T^(n))_(i,j)+(K1(T_(i-1,j)+T_(i+1,j))+K2(T_(i,j-1)+T_(i,j+1)))^(n+1))/(1+2K1+2K2)`

Here `n` denotes the time index and `(i,j)` denotes the space index

initial_cond

Other conditions are specified for the problem:

Thermal Diffusivity `alpha` = `10^-4`

Time-step = 0.01 seconds

Number of grid-points = 51

dx = dy = 0.02

tolerance value = `10^-4`

1. Explicit Iteration Scheme:

(A) MATLAB PROGRAM:

%% MATLAB Program to generate temperature profile of 2D Transient State Heat Conduction Equation explicitly:
clear all;
close all;
clc;

%% Initializing Variables:
Nx = 51;
Ny = Nx;
x = linspace(0,1,Nx);
y = linspace(0,1,Ny);
dx = 1/(Nx-1);
dy = 1/(Ny-1);
alpha = 1e-4;
t_start = 0;
t_end = 1/alpha;
dt = 0.01;
error = 6e9;
tol = 1e-4;

%% Initializing the Temperature grid:
T = 298*ones(Nx,Ny);
T(1,:) = 900; % Bottom Temperature
T(end,:) = 600; % Top Temperature
T(:,1) = 400; % Left Temperature
T(:,end) = 800; % Right Temperature
% Edge refining
T(1,1) = (900+400)/2;
T(1,end) = (900+800)/2;
T(end,1) = (400+600)/2;
T(end,end) = (800+600)/2;

%% CFL Number calculation and adaptive Time-step control:
CFL_Number = alpha*dt/(dx^2);

%% Calculation of Temperature Distribution explicitly using iterative solvers:

k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
[xx,yy] = meshgrid(x,y);
Told = T;
T_prev_dt = Told;
tic;
counter = 1;
for z = 1:20000
    %% Gauss Jacobi explicit Iteration Scheme
    for i = 2:Nx-1
        for j = 2:Ny-1
            term1 = Told(i,j);
            term2 = k1*(Told(i-1,j)-2*Told(i,j)+Told(i+1,j));
            term3 = k2*(Told(i,j-1)-2*Told(i,j)+Told(i,j+1));
            T(i,j) = term1+term2+term3;
        end
    end
    Told = T;
    counter = counter+1;
end
time_counter = toc;
figure(1)
[C,h] = contourf(xx,yy,T);
colorbar;
colormap(jet);
clabel(C,h);
title(sprintf('2D Unsteady state heat conduction (Explicit Method) n Number of iterations: %d (Time: %0.2fs) n Computation Time: %0.4f n CFL Number: %0.4f',counter-1,z*dt,time_counter,CFL_Number));
xlabel('xx');
ylabel('yy');

 (B) RESULTS:

explicit

 2. Implicit Iteration Scheme:

(A) MATLAB CODE:

%% MATLAB Program to solve the unsteady 2D Heat Conduction equation Implictly:
clear all;
close all;
clc;

%% Initializing Variables:
Nx = 51;
Ny = Nx;
x = linspace(0,1,Nx);
y = linspace(0,1,Ny);
dx = 1/(Nx-1);
dy = 1/(Ny-1);
alpha = 1e-4;
t_start = 0;
t_end = 1/alpha;
dt = 0.01;
Nt = [t_start:dt*50:t_end];
error = 6e9;
tol = 1e-4;

%% Initializing the Temperature grid:
T = 298*ones(Nx,Ny);
T(1,:) = 900; % Bottom Temperature
T(end,:) = 600; % Top Temperature
T(:,1) = 400; % Left Temperature
T(:,end) = 800; % Right Temperature
% Edge refining
T(1,1) = (900+400)/2;
T(1,end) = (900+800)/2;
T(end,1) = (400+600)/2;
T(end,end) = (800+600)/2;

%% Calculation of optimal value of over-relaxation parameter (omega):
omega = 2/(1+sin(pi*dx));

%% Calculation of Temperature Distribution explicitly using iterative solvers:

k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
[xx,yy] = meshgrid(x,y);
Told = T;
T_prev_dt = Told;
counter = 1;
solver = input('Enter the number for the dedicated solver type: n 1 for Gauss Jacobi n 2 for Gauss Siedel n 3 for SOR n');
tic;
for k = 1:20000
    %% Gauss Jacobi implicit Iteration Scheme:
    if solver == 1
        term1 = (1+2*k1+2*k2)^(-1);
        term2 = k1*term1;
        term3 = k2*term1;
        error = 1e9;
        while (error > tol)
            for i = 2:Nx-1
                for j = 2:Ny-1
                    H = (Told(i-1,j)+Told(i+1,j));
                    V = (Told(i,j-1)+Told(i,j+1));
                    T(i,j) = term1*T_prev_dt(i,j)+term2*H+term3*V;
                end
            end
            error = max(max(abs(Told-T)));
            Told = T;
            counter = counter + 1;
        end
    end
    if solver == 2
        term1 = (1+2*k1+2*k2)^(-1);
        term2 = k1*term1;
        term3 = k2*term1;
        error = 1e9;
        while (error > tol)
            for i = 2:Nx-1
                for j = 2:Ny-1
                    H = (T(i-1,j)+Told(i+1,j));
                    V = (T(i,j-1)+Told(i,j+1));
                    T(i,j) = term1*T_prev_dt(i,j)+term2*H+term3*V;
                end
            end
            error = max(max(abs(Told-T)));
            Told = T;
            counter = counter + 1;
        end
    end
    if solver == 3
        term1 = (1+2*k1+2*k2)^(-1);
        term2 = k1*term1;
        term3 = k2*term1;
        error = 1e9;
        while (error > tol)
            for i = 2:Nx-1
                for j = 2:Ny-1
                    %H = (T(i-1,j)+Told(i+1,j));
                    %V = (T(i,j-1)+Told(i,j+1));
                    %T(i,j) = (omega)*(term1*T_prev_dt(i,j)+term2*H+term3*V) + (1-omega)*T_prev_dt(i,j);
                    H=(omega*T(i-1,j)+(1-omega)*Told(i-1,j)+Told(i+1,j));
                    V=(omega*T(i,j-1)+(1-omega)*Told(i,j-1)+Told(i,j+1));
                    T(i,j)=T_prev_dt(i,j)*term1+(term2*H+term3*V);
                end
            end
            error = max(max(abs(Told-T)));
            Told = T;
            counter = counter + 1;
        end
    end
    T_prev_dt = T;
end
execution_time = toc;
figure(1);
str_array = ["Gauss-Jacobi","Gauss-Siedel","SOR"];
formatSpec = "Unsteady state 2D Heat Transfer by implicit: (%s Method) n Time: %0.3fs n Iteration Counter: %d n Computation time: %0.3fs";

[C,h] = contourf(xx,yy,T);
colorbar;
colormap(jet);
clabel(C,h);
xlabel('xx');
ylabel('yy');
title(sprintf(formatSpec,str_array(solver),k*dt,counter,execution_time));

 (B) RESULTS:

 i. Gauss Jacobi Method:

imp_jac

 ii. Gauss Siedel Method:

imp_gs

iii. Successive Over Relaxation Method:

imp_sor

STEADY STATE SOLUTION:

Derivation of the numerical scheme:

From the generalized Energy equation in 3 Dimensions we make the following assumptions:

i. No convection occurs

ii. No internal heat generation or dissipation

iii. `(delT)/(delt) = 0`

Hence the equation becomes:

`(del^2T)/(delx^2)+(del^2T)/(dely^2) = 0`

We Discretize the spatial terms by central differencing (`O(h^2)`) i.e.

`(del^2T)/(delx^2) = (T_(i-1,j)-2T_(i,j)+T_(i+1,j))/(Deltax^2)`

`(del^2T)/(dely^2) = (T_(i,j-1)-2T_(i,j)+T_(i,j+1))/(Deltay^2)`

hence:`(T_(i-1,j)-2T_(i,j)+T_(i+1,j))/(Deltax^2) + (T_(i,j-1)-2T_(i,j)+T_(i,j+1))/(Deltay^2) = 0`

Also `Deltax = Deltay`

`T_i = 1/4(T_(i-1,j)+T_(i+1,j)+T_(i,j-1)+T_(i,j+1))`

Number of gridpoints Nx = Ny = 81

other problem specifications are same as before.

SOLVING THE SYSTEM IMPLICTLY:

(A) MATLAB PROGRAM:

%% MATLAB Program to generate the 2D temperature distribution at steady state:
clear all;
close all;
clc;

%% Initial Conditions:
Nx = 81;
Ny = Nx;
x = linspace(0,1,Nx);
y = linspace(0,1,Ny);
dx = 1/(Nx-1);
dy = 1/(Ny-1);

k = 2*(1/(dx^2)+1/(dy^2));
error = 1e8;
tol = 1e-4;

%% Initializing the temperature matrix:

T = 298*ones(Nx,Ny);
T(1,:) = 900; % Bottom Temperature
T(end,:) = 600; % Top Temperature
T(:,1) = 400; % Left Temperature
T(:,end) = 800; % Right Temperature
% Edge refining
T(1,1) = (900+400)/2;
T(1,end) = (900+800)/2;
T(end,1) = (400+600)/2;
T(end,end) = (800+600)/2;

%% Calculation of over-relaxation parameter Omega:

omega = 2/(1+sin(pi*dx)); 
%% Solving the 2D steady state system iteratively:
[xx,yy] = meshgrid(x,y);
Told = T;

solver = input('Enter the desired number: ')

tic;
%% Gauss Jacobi Iteration Scheme:
if solver == 1
    iter_count_jac = 1;
    while (error > tol)
        for i = 2:Nx-1
            for j = 2:Ny-1
                term1 = (Told(i-1,j)+Told(i+1,j))/(k*(dx^2));
                term2 = (Told(i,j-1)+Told(i,j+1))/(k*(dy^2));
                T(i,j) = term1+term2;
            end
        end
        error = max(max(abs(Told-T)));
        Told = T;
        iter_count_jac = iter_count_jac + 1;
    end
end
%% Gauss Siedel Iteration Scheme:
if solver == 2
    iter_count_gs = 1;
    while (error > tol)
        for i = 2:Nx-1
            for j = 2:Ny-1
                term1 = (T(i-1,j)+T(i+1,j))/(k*(dx^2));
                term2 = (T(i,j-1)+T(i,j+1))/(k*(dy^2));
                T(i,j) = term1+term2;
            end
        end
        error = max(max(abs(Told-T)));
        Told = T;
        iter_count_gs = iter_count_gs + 1;
    end
end
%% Successive Over Relaxation (SOR) Iteration Scheme:
if solver == 3
    iter_count_sor = 1;
    while (error > tol)
        for i = 2:Nx-1
            for j = 2:Ny-1
                term1 = (T(i-1,j)+Told(i+1,j))/(k*(dx^2));
                term2 = (T(i,j-1)+Told(i,j+1))/(k*(dy^2));
                 T(i,j) = omega*(term1+term2)-(omega-1)*Told(i,j);
            end
        end
        error = max(max(abs(Told-T)));
        Told = T;
        iter_count_sor = iter_count_sor + 1;
    end
end
time_required = toc;
figure(1)
[C,h] = contourf(xx,yy,T);
colorbar;
colormap(jet);
clabel(C,h);
if solver == 1
    title(sprintf('Temperature profile of steady state heat conduction using Gauss Jacobi method n No. of Iterations: %d n Computation Time: %0.3f s',iter_count_jac,time_required));
else if solver == 2
        title(sprintf('Temperature profile of steady state heat conduction using Gauss Siedel method n No. of Iterations: %d n Computation Time: %0.3f s',iter_count_gs,time_required));
    else if solver == 3
            title(sprintf('Temperature profile of steady state heat conduction using SOR method n No. of Iterations: %d n Computation Time: %0.3f s n Relaxation Parameter: %f',iter_count_sor,time_required,real(omega)));
        end
    end
end
xlabel('xx');
ylabel('yy');

(B) RESULTS:

i. Gauss Jacobi Method:

steady_jac

ii. Gauss Siedel Method:

steady_gs

iii. Successive Over Relaxation Method:

steady_sor

The Optimal Value of `omega` is given by the formula:

`omega_(opt)=2/(1+sin(pidx)`


Projects by Priyotosh Bairagya

LAMINAR INCOMPRESSIBLE FLOW SIMULATION THROUGH A PIPE IN OPENFOAM: [PART 2/2]: SYMMETRY BOUNDARY CONDITION: OBJECTIVES:1. Creating the Mesh-Script (blockMeshDict file) with symmetry boundary condition.2. Simulation Results and post processing the velocity profile at Read more

LAMINAR INCOMPRESSIBLE FLOW SIMULATION THROUGH A PIPE IN OPENFOAM: [PART 1/2]: WEDGE BOUNDARY CONDITION: OBJECTIVES:1. Calculation of Pre-eliminary quantinites related to the flow.2. Creating the Mesh script (blockMeshDict file) for specifying the geometry and bounda Read more

I. Interpolation Schemes in Finite Volume Method: The approximation of surface and volume integrals may require values of the variable at locations other than the computational nodes of the CV. Values at these locations are obtained using interpolation formulae. Some o Read more

ANALYSIS OF NUMERICAL STABILITY OF VARIOUS ITERATIVE SOLVERS FOR TRANSIENT 2D HEAT CONDUCTION: [Part: 3/3] INTRODUCTION: The criterion of stability of a numerical scheme is determined by the way the errors propagate while the solution moves from one time-step to th Read more

MESH GENERATION AND ANALYSIS USING BLOCKMESH FOR FLOW OVER A BACKWARD FACING STEP: The purpose of the following project is to generate the geometry for a variation of the incompressible cavity flow problem in OpenFOAM. For this purpose we have modified the lid-driven c Read more

ANALYSIS OF VARIOUS ITERATIVE SCHEMES FOR THE SOLUTION OF A 2D HEAT CONDUCTION PROBLEM: [PART: 2/3] In the previous part we had explained the problem statement and the MATLAB Program in detail. In this Part we are going to explain the outputs from the 2D Heat Conduct Read more

NUMERICAL SOLUTION OF 1D SUPERSONIC NOZZLE FLOW SIMULATION BY MAC-CORMACK METHOD PROJECT OBJECTIVES: i. Numerical solution of the governing equations in both conservative and non-conservative forms. ii. Creating user defined functions for calculating the flow quan Read more

UNDERSTANDING LINEAR SYSTEMS(ANALYSIS OF VARIOUS ITERATIVE SCHEMES TO SOLVE A SYSTEM OF LINEAR EQUATIONS TO FIND THE EIGEN VALUES AND SPECTRAL RADIUS) (A) PROBLEM STATEMENT: Given coefficient matrix: `A = [[5,1,2],[-3,9,4],[1,2,-7]]` Given Solution Matrix: `X = [[x Read more

WEEK 4: (Effect of Grid-Size on output for the solution of 1D linear wave equation) 1. Problem Setup: Given Partial Differential Equation: `(delu)/(delt) + c(delu)/(delx) = 0` Numerical Discretization:  `u_(i,n+1) = u_(i,n) + (cDeltat)/(Deltax)*(u_(i-1,n+1) Read more


Loading...

The End