Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem

AIM: -

To solve the 2D heat conduction problem for the steady and unsteady state using Jacobi, gauss seidel, and successive over-relaxation iterative techniques.

MAT LAB CODE: -

There is 2 MATLAB code written separately for steady and unsteady conditions. But for both conditions’ inputs are the same, while programming it is assumed that, the domain is the surface of area (1mX1m). Discretized into 10 nodes both horizontally and vertically. Constant values like alpha and time step 'data, by keeping the explicit condition in mind. Relaxation factor omega is also assumed while programming only.

Inputs: -

Top boundary temperature	600 K
Bottom Boundary temperature	900 K
Left Boundary temperature	400 K
Right Boundary temperature	800 K
error	1e-4

Steady-State 2D Heat conduction: -

 

<!-- [if !supportLists]-->·         <!--[endif]-->First inputs are defined into variables and the temperature array is created.

<!-- [if !supportLists]-->·         <!--[endif]-->Boundary conditions are introduced into the temperature array, and at corners, boundary values are averaged.

<!-- [if !supportLists]-->·         <!--[endif]-->Constants like alpha, omega, tolerance are defined to run the iteration.

<!-- [if !supportLists]-->·         <!--[endif]-->Input must be provided by the user, to choose the iterative method.

<!-- [if !supportLists]-->·         <!--[endif]-->The if-else loop is used to solve the equation using the selected (input from the user) iterative technique.

<!-- [if !supportLists]-->·         <!--[endif]-->while loop is used to run the iterative technique till the required accuracy is achieved.

<!-- [if !supportLists]-->·         <!--[endif]-->Two for loop is used to calculate the temperature values at every node at each iterative step.

<!-- [if !supportLists]-->·         <!--[endif]-->Finally, the temperature distribution graph is obtained using 'conturf' command.

MAT LAB Code

 

% Domain and discretization length

number_node = 10;

domain_length = 1;

 

%value of x and y at nodes in a discretization plane

x = linspace (0, domain_length,number_node);

y = linspace (0, domain_length,number_node);

dx = domain_length/(number_node-1);

dy = domain_length/(number_node-1);

 

 

% Let define a temperature value at all nodes with initial guess 300k

T = 300*ones(number_node,number_node);

 

% Initial boundary conditions at left and right

T (1: end,1) = 400;

T (1: end,end) = 800;

 

% Initial boundary conditions at bottom and top

T (end,1: end) = 900;

T (1,1: end) = 600;

 

%average values of T at corner boundaries

T (1,1) = 500;

T(end,1) = 650;

T (1, end) = 700;

T (end, end) = 850;

 

%copy of Temperature to compute to compute

T_old = T;

 

%constant

K = (2*(dx^2+dy^2))/(dx^2*dy^2);

omega = 1.2;

 

% to enter while loop initialize error value randomly

error = 10;

tol = 1e-4;

no_iteration = 0;

 

method = input ('Enter (1) for Jacobi, (2) for gauss seidel, (3) for sos');

if method == 1

    while(error>tol)

        for i = 2: number_node-1

            for j = 2: number_node-1

            term1 = (T_old (i+1, j) +T_old(i-1, j))/dx^2;

            term2 = (T_old (i, j+1) +T_old(i, j-1))/dy^2;

            T(i,j) = (term1+term2)/K;

            end

        end

        error = max(max(abs(T_old-T)));

       

        % to upadte the previous velcity with new velocity

        %got from the iteration

        T_old = T;

        no_iteration = no_iteration+1;

       

    end

elseif method == 2

    while(error>tol)

        for i = 2:number_node-1

            for j = 2:number_node-1

            term1 = (T_old(i+1,j)+T(i-1,j))/dx^2;

            term2 = (T_old(i,j+1)+T(i,j-1))/dy^2;

            T(i,j) = (term1+term2)/K;

            end

        end

        error = max(max(abs(T_old-T)));

       

        % to upadte the previous velcity with new velocity

        %got from the iteration

        T_old = T;

        no_iteration = no_iteration+1;

    end

   

else

    while(error>tol)

        for i = 2:number_node-1

            for j = 2:number_node-1

            term1 = (T_old(i+1,j)+T(i-1,j))/dx^2;

            term2 = (T_old(i,j+1)+T(i,j-1))/dy^2;

            T(i,j) = ((1-omegha)*T_old(i,j))+(omegha*((term1+term2)/K));

            end

        end

        error = max(max(abs(T_old-T)));

       

        % to update the previous velocity with new velocity

        %got from the iteration

        T_old = T;

        no_iteration = no_iteration+1;

       

    end

end

 

% to visualize the result

 

[X,Y] = meshgrid(x,y);

contourf(X,Y,T)

colorbar

set(gca,'ydir','reverse')

Unsteady-State 2D Heat conduction: -

Inputs are defined into variables and the temperature array is created. Boundary conditions are introduced into the temperature array, and at corners, boundary values are averaged. Constants like alpha, omega, tolerance are defined to run the iteration. Here a user must provide, one or two inputs. The first input is to whether solve the problem using an explicit approach or an implicit approach. If an implicit approach is selected, then the user must provide another input, to select the iterative technique. The if-else loop is used to solve the equation using the selected (input from the user) approach and iterative technique. For loop is used to run the iteration for a large number of times to reach steady-state conditions. While loop is used to run the iterative technique till the required accuracy is achieved. Two for loop is used to calculate the temperature values at every node at each iterative step. Finally, a temperature distribution graph is obtained using 'conturf' command.

MAT LAB Code

% Domain and discretization length

number_node = 10;

domain_length = 1;

 

%value of x and y at nodes in a discretization plane

x = linspace(0,domain_length,number_node);

y = linspace(0,domain_length,number_node);

dx = (domain_length/number_node-1);

dy = (domain_length/number_node-1);

dt = 1/10;

 

 

% Let define a temperature value at all nodes with initial guess 300k

T = 300*ones(number_node,number_node);

 

% Initia boundary conditions at left and right

T(1:end,1) = 400;

T(1:end,end) = 800;

% Initia boundary conditions at bottom and top

T(end,1:end) = 900;

T(1,1:end) = 600;

%average values of T at corner boundaries

T(1,1) = 500;

T(end,1) = 650;

T(1,end) = 700;

T(end,end) = 850;

 

%copy of Temperature to compute for previous space terms

T_old = T;

 

% copy of T for previous time_step

T_prev = T;

 

%constants

alpha = 1.5;

omegha = 1.15;

K1 = (alpha*dt)/(dx^2);

K2 = (alpha*dt)/(dy^2);

iteration = 0;

 

 

approach = input('Enter 1 for implicit or Enter 2 for explicit')

 

if approach == 1

    method = input('Enter 1 for jacobi or 2 for gauss siedel or 3 for sos')

    if method == 1

        for nt = 1:1500

            error = 10;

            tol = 1e-4;

            while(error>tol)

                for i = 2:number_node-1

                    for j = 2:number_node-1

                        term1 = 1/(1+(2*K1)+(2*K2));

                        term2 = K1*(T_old(i+1,j)+T_old(i-1,j));

                        term3 = K2 *(T_old(i,j+1)+T_old(i,j-1));

                        T(i,j) = term1*(T_prev(i,j)+term2+term3);

                    end

                end

                iteration = iteration+1;

                error = max(max((abs(T_old-T))));

                T_old = T;

            end

            T_prev = T;

        end

    elseif method ==2

        for nt = 1:1500

            error = 10;

            tol = 1e-4;

            while(error>tol)

                for i = 2:number_node-1

                    for j = 2:number_node-1

                        term1 = 1/(1+(2*K1)+(2*K2));

                        term2 = K1*(T_old(i+1,j)+T(i-1,j));

                        term3 = K2 *(T_old(i,j+1)+T(i,j-1));

                        T(i,j) = term1*(T_prev(i,j)+term2+term3);

                    end

                end

                iteration = iteration+1;

                error = max(max(abs(T_old-T)));

                T_old = T;

               

            end

            T_prev = T;

        end

    else

        for nt = 1:1500

            error = 10;

            tol = 1e-4;

            while(error>tol)

                for i = 2:number_node-1

                    for j = 2:number_node-1

                        term1 = 1/(1+(2*K1)+(2*K2));

                        term2 = K1*(T_old(i+1,j)+T(i-1,j));

                        term3 = K2 *(T_old(i,j+1)+T(i,j-1));

                        term4 = term1*(T_prev(i,j)+term2+term3);

                        T(i,j)=(1-omegha)*T_old(i,j)+(omegha*term4);

                    end

                end

                iteration = iteration+1;

                error = max(max(abs(T_old-T)));

                T_old = T;

            end

            T_prev = T;

        end

    end

else

     for nt = 1:1500

        for i = 2:number_node-1

            for j = 2:number_node-1

                term1 = K1*(T_prev(i+1,j)-2*T_prev(i,j)+T_prev(i-1,j));

                term2 = K2*(T_prev(i,j+1)-2*T_prev(i,j)+T_prev(i,j-1));

                T(i,j) = T_prev(i,j)+term1+term2;

            end

        end

        % to update the previous timestep value with new timestep value

        %got from the iteration

        T_prev = T;

        iteration = iteration+1;

    end

end

 

% to visualize the result

 

[X,Y] = meshgrid(x,y);

figure(1)

contourf(X,Y,T)

colorbar

set(gca,'ydir','reverse')

RESULT: -

Steady-State analysis:

<!-- [if !supportLists]-->·         <!--[endif]-->Jacobi

<!-- [if !supportLists]-->·         <!--[endif]-->Gauss-Siedel

<!-- [if !supportLists]-->·         <!--[endif]-->Successive over relaxation

Unsteady analysis:

Out of a number of iterations,1500 iterations is due to time step

<!-- [if !supportLists]-->·         <!--[endif]-->Jacobi

<!-- [if !supportLists]-->·         <!--[endif]-->Gauss-Siedel

<!-- [if !supportLists]-->·         <!--[endif]-->Successive over relaxation

 

Explicit

1500 iterations are due to time step

According to CFL stability criteria, K1+K2<0>