Solving equations of two dimensional heat conduction and Laplace equation with no internal heat generation on square plate using Matlab/Octave

Objective 

Simulating two dimensional Fourier biot equation with no internal heat generation on a square plate with the temperature of 400,800,900,600 degrees Celcius on left, right, bottom, top respectively. To solve this differential equation finite difference method is being used, and to solve linear equations gauss-seidel, Jacobi, a successive over-relaxation method is being used.

Steps to solve

Unsteady 2d heat equation or fourier biot 2d heat equation

euation

`(dT)/dt=alpha*((d^2T)/dx^2+(d^2T)/dy^2)`

`alpha=k/(cp*rho)`

`k= thermal conductivity`

cp= specific heat capacity

rho = dencity of that material.

•  discetize forier equation

implicite (n+1) term of time with forward in time central in space 

`T^(n+1)(i,j)=Tp^(n+1)`

`Tp(n+1)(i,j)=(1+2k1+2k2)^-1*((Tp^(n)(i,j))+[((T(i-1,j))+T(i+1,j))+((T(i,j-1))+T(i,j+1))]^(n+1))`

`k1=alpha*(dt/dx^2)`

`k2=alpha*(dt/dy^2)`

Jacobi iteration 

clear all
close all
clc

%dT/dt-alpha(d^2(T)/dx^2+d^2(T)/dy^2)=0 


% number of grid point in x direction
nx = 10;
% number of grid point in y direction
ny=nx;



x=linspace(0,1,nx); % 0 to 1 equal dirubuted values are 10 beacause nx is 10 size 1*10
y=linspace(0,1,ny);

% assign  bouundary condition
T(ny,nx)=zeros; % size arry 10*10 all value are 1



% all row 1 st column left 400
T(:,1)=400; 
% all row last last column right 800
T(:,nx)=800;
% first row all column top 600 temp
T(1,:)=600; 
% last row all column bottom 900 temp
T(nx,:)=900;

% corner points average
T(1,1)=(600+400)*0.5;
T(end,1)=(900+400)*0.5;
T(end,end)=(800+900)*0.5;
T(1,end)=(600+800)*0.5;
%copy of T
Told=T;
Tprev_dt=Told;

% thermal diffusivity k/rho*cp 
alpha=0.5;

% space between two point on x or y axis 
dx=x(2)-x(1);
dy=y(2)-y(1);

% time increment
dt=0.9;

%cfl
k1=alpha*(dt/dx^2);
k2=alpha*(dt/dy^2);




%  error
error=9e9;
%jacobi iterative solver
% till tolerance

tol=1e-4;

% time loop


for w=1:3

      % convergence loop
      
                  % initial jacobi iter 
                     jacobi_iter=1;
      while (error>tol)
                     % space loop    
           for i=2:nx-1  
        
               for j=2:ny-1
                   term1=((1+(2*k1))+(2*k2))^-1;
                   term2=k1*term1;
                   term3=k2*term1;
                   H=(Told(i-1,j))+Told(i+1,j);
                   V=(Told(i,j-1))+Told(i,j+1);
                   T(i,j)=(Tprev_dt(i,j)*term1)+(H*term2)+(V*term3);
               end
           end
        
           error=max(max(abs(Told-T)))
          
    
           % updating error term to next itervation
           Told=T;
           jacobi_iter=jacobi_iter+1;
           contourf(x,y,T)
           colorbar
           xlabel(\'x axis\')
           ylabel(\'yaxis\')
    
           title_text = sprintf(\'iteration number=%d\',jacobi_iter);
           title(title_text)
           drawnow
        end
       
       
    T=Tprev_dt;
end

output

The implicit method here is stable and gives iteration of 196 for Jacobi which more than Gauss and successive over relaxation.

k1 and k2 = 36.4500

gauss seidel itertion in transient heat conduction in implicite discetization

clear all
close all
clc

%dT/dt-alpha(d^2(T)/dx^2+d^2(T)/dy^2)=0 


% number of grid point in x direction
nx = 10;
% number of grid point in y direction
ny=nx;



x=linspace(0,1,nx); % 0 to 1 equal dirubuted values are 10 beacause nx is 10 size 1*10
y=linspace(0,1,ny);

% assign  bouundary condition
T(ny,nx)=zeros; % size arry 10*10 all value are 1



% all row 1 st column left 400
T(:,1)=400; 
% all row last last column right 800
T(:,nx)=800;
% first row all column top 600 temp
T(1,:)=600; 
% last row all column top 900 temp
T(nx,:)=900;

% corner points average
T(1,1)=(600+400)*0.5;
T(end,1)=(900+400)*0.5;
T(end,end)=(800+900)*0.5;
T(1,end)=(600+800)*0.5;
%copy of T
Told=T;
Tprev_dt=Told;

% thermal diffusivity k/rho*cp 
alpha=0.5;

% space between two point on x or y axis 
dx=x(2)-x(1);
dy=y(2)-y(1);

% time increment
dt=0.9;

%cfl
k1=alpha*(dt/dx^2);
k2=alpha*(dt/dy^2);




%  error
error=9e9;
%jacobi iterative solver
% till tolerance

tol=1e-4;

% time loop
% forward 3 steps

for w=1:3

      % convergence loop
      
                  % initial jacobi iter 
                     gs_implicite=1;
      while (error>tol)
                     % space loop    
           for i=2:nx-1  
        
               for j=2:ny-1
                   term1=((1+(2*k1))+(2*k2))^-1;
                   term2=k1*term1;
                   term3=k2*term1;
                   H=(T(i-1,j))+T(i+1,j);
                   V=(T(i,j-1))+T(i,j+1);
                   T(i,j)=(T(i,j)*term1)+(H*term2)+(V*term3);
               end
           end
        
           error=max(max(abs(Told-T)))
          
    
           % updating error term to next itervation
           Told=T;
           gs_implicite=gs_implicite+1;
           contourf(x,y,T)
           colorbar
    
           title_text = sprintf(\'iteration number=%d\',gs_implicite);
           title(title_text)
           drawnow
        end
       
       
    T=Tprev_dt;
end

plotting

here a number of the iteration is less because previously calculated values of T are used to find next iteration 

Sucessive over relaxation implicite on transient

`T^(n+1)(i,j) = (1-w)T(i,j)+w*[(1+2*k1+2*k2)^(-1)][T(i,j)+k1*T(i-1,j)+k1*(i+1,j)+k2*(i,j-1)+k2*(i,j+1)]`

clear all
close all
clc

%dT/dt-alpha(d^2(T)/dx^2+d^2(T)/dy^2)=0 


% number of grid point in x direction
nx = 10;
% number of grid point in y direction
ny=nx;



x=linspace(0,1,nx); % 0 to 1 equal dirubuted values are 10 beacause nx is 10 size 1*10
y=linspace(0,1,ny);

% assign  bouundary condition
T(ny,nx)=zeros; % size arry 10*10 all value are 1



% all row 1 st column left 400
T(:,1)=400; 
% all row last last column right 800
T(:,nx)=800;
% first row all column top 600 temp
T(1,:)=600; 
% last row all column top 900 temp
T(nx,:)=900;

% corner points average
T(1,1)=(600+400)*0.5;
T(end,1)=(900+400)*0.5;
T(end,end)=(800+900)*0.5;
T(1,end)=(600+800)*0.5;
%copy of T
Told=T;
Tprev_dt=Told;

% thermal diffusivity k/rho*cp 
alpha=0.5;

% space between two point on x or y axis 
dx=x(2)-x(1);
dy=y(2)-y(1);

% time increment
dt=0.9;

%cfl
k1=alpha*(dt/dx^2);
k2=alpha*(dt/dy^2);


%omega w>1
w=1.5;
%  error
error=9e9;
%jacobi iterative solver
% till tolerance

tol=1e-4;

% time loop


for Y=1:3

      % convergence loop
      
                  % initial jacobi iter 
                     gs_implicite=1;
      while (error>tol)
                     % space loop    
           for i=2:nx-1  
        
               for j=2:ny-1
                   term1=((1+(2*k1))+(2*k2))^-1;
                   term2=k1*term1;
                   term3=k2*term1;
                   H=(T(i-1,j))+T(i+1,j);
                   V=(T(i,j-1))+T(i,j+1);
                   J(i,j)=(T(i,j)*term1)+(H*term2)+(V*term3);
                   T(i,j)=((1-w)*T(i,j))+(w*J(i,j));
               end
           end
        
           error=max(max(abs(Told-T)))
          
    
           % updating error term to next itervation
           Told=T;
           gs_implicite=gs_implicite+1;
           contourf(x,y,T)
           colorbar
    
           title_text = sprintf(\'iteration number=%d\',gs_implicite);
           title(title_text)
           drawnow
        end
       
       
    T=Tprev_dt;
end

plotting

successive over-relaxation required less number of iteration than GS and Jacobi.

explicite (n+1) term of time with forward in time central in space 

`T^(n+1)(i,j)=Tp^(n+1)`

`Tp(n+1)(i,j)=(1-2k1-2k2)*[((Tp^(n)(i,j))]+[{((T(i-1,j))+T(i+1,j))}*k1+*k2{((T(i,j-1))+{T(i,j+1))}]^(n))`

k1=alpha*(dt/dx^2)

k2=alpha*(dt/dy^2)

explicite jacobi in transient state heat equation   grid points 10 

clear all
close all
clc

%dT/dt-alpha(d^2(T)/dx^2+d^2(T)/dy^2)=0 


% number of grid point in x direction
nx = 10;
% number of grid point in y direction
ny=nx;



x=linspace(0,1,nx); % 0 to 1 equal dirubuted values are 10 beacause nx is 10 size 1*10
y=linspace(0,1,ny);

% assign  bouundary condition
T(ny,nx)=zeros; % size arry 10*10 all value are 1



% all row 1 st column left 400
T(:,1)=400; 
% all row last last column right 800
T(:,nx)=800;
% first row all column top 600 temp
T(1,:)=600; 
% last row all column top 900 temp
T(nx,:)=900;

% corner points average
T(1,1)=(600+400)*0.5;
T(end,1)=(900+400)*0.5;
T(end,end)=(800+900)*0.5;
T(1,end)=(600+800)*0.5;
%copy of T
Told=T;
Tprev_dt=Told;

% thermal diffusivity k/rho*cp (aluminium)
alpha=0.5;

% space between two point on x or y axis 
dx=x(2)-x(1);
dy=y(2)-y(1);

% time increment
dt=0.0062;

%cfl
k1=alpha*(dt/dx^2);
k2=alpha*(dt/dy^2);




%  error
error=9e9;
%jacobi iterative solver
% till tolerance

tol=1e-4;

% time loop 30 cycle
tic;
for w=1:30

      % convergence loop
      
                  % initial jacobi iter 
                     jacob_iter_expli=1;
      while (error>tol)
                     % space loop    
           for i=2:nx-1  
        
               for j=2:ny-1
                   term1=((1-(2*k1))-(2*k2));
                   H=(Told(i-1,j))+Told(i+1,j);
                   V=(Told(i,j-1))+Told(i,j+1);
                   
                   T(i,j)=(Tprev_dt(i,j)*(term1))+k1*(H)+k2*(V);
               end
           end
        
           error=max(max(abs(Told-T)))
          
    
           % updating error term to next itervation
           Told=T;
           jacob_iter_expli=jacob_iter_expli+1;
           contourf(x,y,T)
           colorbar
            xlabel('x axis')
           ylabel('y axis')
           title_text = sprintf('iteration number=%d',jacob_iter_expli);
           title(title_text)
           drawnow
        end
       
     T=Tprev_dt;
    
end
t1_jacob_explicite_trans=toc;


figure(2)
bar(t1_jacob_explicite_trans)
ylabel('simulation time ')

title('jacobi iteration in unsteady simulation time(explicite)')

PLOTTING of explicite jacobi transient 2d heat 

here solution looks like it gets stable

Explicit solution using  gauss using grid point 10

clear all
close all
clc

%dT/dt-alpha(d^2(T)/dx^2+d^2(T)/dy^2)=0 


% number of grid point in x direction
nx = 10;
% number of grid point in y direction
ny=nx;



x=linspace(0,1,nx); % 0 to 1 equal dirubuted values are 10 beacause nx is 10 size 1*10
y=linspace(0,1,ny);

% assign  bouundary condition
T(ny,nx)=zeros; % size arry 10*10 all value are 1



% all row 1 st column left 400
T(:,1)=400; 
% all row last last column right 800
T(:,nx)=800;
% first row all column top 600 temp
T(1,:)=600; 
% last row all column top 900 temp
T(nx,:)=900;

% corner points average
T(1,1)=(600+400)*0.5;
T(end,1)=(900+400)*0.5;
T(end,end)=(800+900)*0.5;
T(1,end)=(600+800)*0.5;
%copy of T
Told=T;


% thermal diffusivity k/rho*cp (aluminium)
alpha=0.5;

% space between two point on x or y axis 
dx=x(2)-x(1);
dy=y(2)-y(1);

% time increment
dt=0.0062;

%cfl
k1=alpha*(dt/dx^2);
k2=alpha*(dt/dy^2);




%  error
error=9e9;
%jacobi iterative solver
% till tolerance

tol=1e-4;

% time loop

for w=1:3

      % convergence loop
      
                  % initial jacobi iter 
                     expli_gs_iter=1;
      while (error>tol)
                     % space loop    
           for i=2:nx-1  
        
               for j=2:ny-1
                   term1=((1-(2*k1))-(2*k2));
                   H=(T(i-1,j))+T(i+1,j);
                   V=(T(i,j-1))+T(i,j+1);
                   
                   T(i,j)=(T(i,j)*(term1))+k1*(H)+k2*(V);
               end
           end
        
           error=max(max(abs(Told-T)))
          
    
           % updating error term to next itervation
           Told=T;
           expli_gs_iter=expli_gs_iter+1;
           contourf(x,y,T)
           colorbar
           xlabel(\'x axis\')
           ylabel(\'y axis\')
           title_text = sprintf(\'iteration number=%d\',expli_gs_iter);
           title(title_text)
           drawnow
        end
       
       T=Told;
    
end

plotting (explicit with transient gauss iteration method)

here again, compare to implicit its stable less iteration than Jacobi used for explicit solver\\.

explicit solver and successive over-relaxation using gauss with w=1.5

`T(i,j)=(1-w)T(i,j)+(w)*[(T(i,j)(1-2*k1-2*k2)+k1*(T(i-1,j)-T(i+1,j))+k2*(T(i,j-1)-T(i,j+1))}`

clear all
close all
clc

%dT/dt-alpha(d^2(T)/dx^2+d^2(T)/dy^2)=0 


% number of grid point in x direction
nx = 10;
% number of grid point in y direction
ny=nx;



x=linspace(0,1,nx); % 0 to 1 equal dirubuted values are 10 beacause nx is 10 size 1*10
y=linspace(0,1,ny);

% assign  bouundary condition
T(ny,nx)=zeros; % size arry 10*10 all value are 1



% all row 1 st column left 400
T(:,1)=400; 
% all row last last column right 800
T(:,nx)=800;
% first row all column top 600 temp
T(1,:)=600; 
% last row all column top 900 temp
T(nx,:)=900;

% corner points average
T(1,1)=(600+400)*0.5;
T(end,1)=(900+400)*0.5;
T(end,end)=(800+900)*0.5;
T(1,end)=(600+800)*0.5;
%copy of T
Told=T;


% thermal diffusivity k/rho*cp (aluminium)
alpha=0.5;

% space between two point on x or y axis 
dx=x(2)-x(1);
dy=y(2)-y(1);

% time increment
dt=0.062;

%cfl
k1=alpha*(dt/dx^2);
k2=alpha*(dt/dy^2);

%sucessive over relaxation value w>1
w=1.5;


%  error
error=9e9;
%jacobi iterative solver
% till tolerance

tol=1e-4;

% time loop

for O=1:3

      % convergence loop
      
                  % initial jacobi iter 
                     expli_gs_iter=1;
      while (error>tol)
                     % space loop    
           for i=2:nx-1  
        
               for j=2:ny-1
                   term1=((1-(2*k1))-(2*k2));
                   H=(T(i-1,j))+T(i+1,j);
                   V=(T(i,j-1))+T(i,j+1);
                   J(i,j)=(T(i,j)*(term1))+k1*(H)+k2*(V);
                   T(i,j)=(1-w)*T(i,j)+(w*(J(i,j));
               end
           end
        
           error=max(max(abs(Told-T)))
          
    
           % updating error term to next itervation
           Told=T;
           expli_gs_iter=expli_gs_iter+1;
           contourf(x,y,T)
           colorbar
           xlabel(\'x axis\')
           ylabel(\'y axis\')
           title_text = sprintf(\'iteration number=%d\',expli_gs_iter);
           title(title_text)
           drawnow
        end
       
       T=Told;
    
end

Plotting(gauss on SOR in a transient equation)

above plot is also stable 

equation for steady state 

`T(i,j)=(1/k)*(((Told(i+1,j)+Told(i-1,j))/(dx^2)+(Told(i,j+1)+Told(i,j-1))/(dy)^2))`

`k=(2/(dx)^2)+(2/(dy)^2)`

using jacobi iteration

clear all
close all
clc

%(d^2(T)/dx^2+d^2(T)/dy^2)=0 


% number of grid point in x direction
nx = 10;
% number of grid point in y direction
ny=nx;



x=linspace(0,1,nx); % 0 to 1 equal dirubuted values are 10 beacause nx is 10 size 1*10
y=linspace(0,1,ny);

% assign  bouundary condition
T(ny,nx)=300*ones; % size arry 10*10 all value are 300kelivn



% all row 1 st column left 400
T(:,1)=400; 
% all row last last column right 800
T(:,nx)=800;
% first row all column top 600 temp
T(1,:)=600; 
% last row all column top 900 temp
T(nx,:)=900;

% corner points average
T(1,1)=(600+400)*0.5;
T(end,1)=(900+400)*0.5;
T(end,end)=(800+900)*0.5;
T(1,end)=(600+800)*0.5;
%copy of T
Told=T;




% space between two point on x or y axis 
dx=x(2)-x(1);
dy=y(2)-y(1);

k=(2/(dx)^2)+(2/(dy)^2);

%  error
error=9e9;
%jacobi iterative solver
% till tolerance

tol=1e-4;



      % convergence loop
      
                  % initial jacobi iter 
                     jacobi_iter=1;
      while (error>tol)
                     % space loop    
           for i=2:nx-1  
        
               for j=2:ny-1
                  T(i,j)=(1/k)*(((Told(i+1,j)+Told(i-1,j))/(dx^2)+(Told(i,j+1)+Told(i,j-1))/(dy)^2));
               end
           end
        
           error=max(max(abs(Told-T)));
          
    
           % updating error term to next itervation
           Told=T;
           jacobi_iter=jacobi_iter+1;
           contourf(x,y,T)
           colorbar
           xlabel(\'x axis\')
           ylabel(\'y axis\')
           title_text = sprintf(\'iteration number=%d\',jacobi_iter);
           title(title_text)
           drawnow
        end
       

plot

here iteration is 217 for Jacobi steady state

gauss seidel on steady-state 

clear all
close all
clc

%(d^2(T)/dx^2+d^2(T)/dy^2)=0 


% number of grid point in x direction
nx = 10;
% number of grid point in y direction
ny=nx;



x=linspace(0,1,nx); % 0 to 1 equal dirubuted values are 10 beacause nx is 10 size 1*10
y=linspace(0,1,ny);

% assign  bouundary condition
T(ny,nx)=300*ones; % size arry 10*10 all value are 300kelivn



% all row 1 st column left 400
T(:,1)=400; 
% all row last last column right 800
T(:,nx)=800;
% first row all column top 600 temp
T(1,:)=600; 
% last row all column top 900 temp
T(nx,:)=900;

% corner points average
T(1,1)=(600+400)*0.5;
T(end,1)=(900+400)*0.5;
T(end,end)=(800+900)*0.5;
T(1,end)=(600+800)*0.5;
%copy of T
Told=T;




% space between two point on x or y axis 
dx=x(2)-x(1);
dy=y(2)-y(1);

k=(2/(dx)^2)+(2/(dy)^2);

%  error
error=9e9;
%jacobi iterative solver
% till tolerance

tol=1e-4;



      % convergence loop
      
                  % initial jacobi iter 
                     jacobi_iter=1;
      while (error>tol)
                     % space loop    
           for i=2:nx-1  
        
               for j=2:ny-1
                  T(i,j)=(1/k)*(((T(i+1,j)+T(i-1,j))/(dx)^2)+(T(i,j+1)+T(i,j-1))/(dy)^2);
               end
           end
        
           error=max(max(abs(Told-T)));
          
    
           % updating error term to next itervation
           Told=T;
           jacobi_iter=jacobi_iter+1;
           contourf(x,y,T)
           colorbar
           xlabel(\'x axis\')
           ylabel(\'y axis\')
           title_text = sprintf(\'iteration number=%d\',jacobi_iter);
           title(title_text)
           drawnow
        end
       

Plotting

here gauss methods take less iteration than Jacob, for the steady-state 2d heat equation

steady-state solving with sor on gauss 

`T(i,j)=(1-w)T+Tgs(w)`

clear all
close all
clc

%(d^2(T)/dx^2+d^2(T)/dy^2)=0 


% number of grid point in x direction
nx = 10;
% number of grid point in y direction
ny=nx;



x=linspace(0,1,nx); % 0 to 1 equal dirubuted values are 10 beacause nx is 10 size 1*10
y=linspace(0,1,ny);

% assign  bouundary condition
T(ny,nx)=300*ones; % size arry 10*10 all value are 300kelivn



% all row 1 st column left 400
T(:,1)=400; 
% all row last last column right 800
T(:,nx)=800;
% first row all column top 600 temp
T(1,:)=600; 
% last row all column top 900 temp
T(nx,:)=900;

% corner points average
T(1,1)=(600+400)*0.5;
T(end,1)=(900+400)*0.5;
T(end,end)=(800+900)*0.5;
T(1,end)=(600+800)*0.5;
%copy of T
Told=T;


%omega;for sor 
w=1.4;

% space between two point on x or y axis 
dx=x(2)-x(1);
dy=y(2)-y(1);

k=(2/(dx)^2)+(2/(dy)^2);

%  error
error=9e9;
%jacobi iterative solver
% till tolerance

tol=1e-4;



      % convergence loop
      
                  % initial sor on gs 
                     sor_gs_iter=1;
      while (error>tol)
                     % space loop    
           for i=2:nx-1  
        
               for j=2:ny-1
                   
                  T(i,j)=(T(i,j)*(1-w))+(w*((1/k)*(((T(i+1,j)+T(i-1,j))/(dx)^2)+(T(i,j+1)+T(i,j-1))/(dy)^2)));
               end
           end
        
           error=max(max(abs(Told-T)));
          
    
           % updating error term to next itervation
           Told=T;
           sor_gs_iter=sor_gs_iter+1;
           contourf(x,y,T)
           colorbar
           xlabel(\'x axis\')
           ylabel(\'y axis\')
           title_text = sprintf(\'iteration number=%d\',sor_gs_iter);
           title(title_text)
           drawnow
        end
       

Plotting

here less number of iteration than all other methods.

conclusion

1. transient state = implicit solver is stable at alliterative method (less iteration means less computing time )
2. transient state = explicit solver is stable  at alliterative method because `dt<=(dx)^2/(4*alpha)`
3. Here learned that implicit solution anyway stable.

1. steady-state = Jacobi takes more iteration than all other methods and hence high computing time
2. SOR iterative solver takes less iteration to achieve final state and hence less computing time required for it.