BlogFor BusinessTalk to a counsellor
CoursesIIT ProgramsPlacementsLog in
BlogFor Business
Talk to a counsellor
Challenge

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

PPratyush Ranjan Meher28 Dec 2021

The equation given is:`(delT)/(delt)+alpha((del^2T)/(delx^2)+(del^2T)/(dely^2))=0`

For Steady state case, this will reduce to:`(del^2T)/(delx^2)+(del^2T)/(dely^2)=0`

The assumtions for solving the problem is:

nx=ny=10

Tempertaure for diffrent boundaries:

  • Left=400K
  • Top=600K
  • Right=800K
  • Bottom=900K

Tolerance=1e-4

Starting Error=9e9

`alpha`=1.88-5 `m^2/s`

The Steady state equation after applying central differencing method became:

`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+dy^2)/(dx^2*dy^2))`

Code for Solving the equations(steady state)

close all 
clc

nx=10;
ny=nx;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
[xx,yy]=meshgrid(x,y);

T=ones(nx,ny);
T(1,:)=900;
T(end,:)=600;
T(:,1)=400;
T(:,end)=800;
Told=T;
err_tot=1e-4;
strt_err=9e9;
k=2*((dx^2+dy^2)/(dx^2*dy^2));
T_jcob=jcob(strt_err,err_tot,xx,yy,nx,ny,dx,dy,T,k,Told);
T_gs=gs(strt_err,err_tot,x,y,nx,ny,dx,dy,T,k,Told);
T_sor=sor(strt_err,err_tot,x,y,nx,ny,dx,dy,T,k,Told);

Function code for SOR(steady)

function T=sor(strt_err,err_tol,x,y,nx,ny,dx,dy,T,k,Told)
sor_iter=1;
omega=1.7;
while(strt_err>err_tol)
    for i=2:nx-1
        for j=2:ny-1
             T_gs=(((T(i-1,j)+Told(i+1,j))/dx^2)+((T(i,j-1)+Told(i,j+1))/dy^2))/k;
             T(i,j)=omega*(T_gs)+(1-omega)*(Told(i,j));
        end
    end
    strt_err=max(max(abs(Told-T)));
    Told=T;
    sor_iter=sor_iter+1;
    figure(3)
    [c,d]=contourf(x,y,T);
    colormap(jet)
    colorbar
    title_txt=sprintf('SOR iteration number %d in staedy state',sor_iter);
    title(title_txt)
    xlabel('X-Axis')
    ylabel('Y-Axis')
    clabel(c,d)
    pause(0.03)
end
fprintf('No of iterations in SOR %d',sor_iter);
end

Function code for Jacobi(steady)

function T=jcob(strt_err,err_tol,xx,yy,nx,ny,dx,dy,T,k,Told)
jcobi_iter=1;
while(strt_err>err_tol)
    for i=2:nx-1
        for j=2:ny-1
            T(i,j)=(((Told(i-1,j)+Told(i+1,j))/dx^2)+((Told(i,j-1)+Told(i,j+1))/dy^2))/k;
        end
    end
    strt_err=max(max(abs(Told-T)));
    Told=T;
    jcobi_iter=jcobi_iter+1;
    figure(1)
    [c,d]=contourf(xx,yy,T);
    colormap(jet)
    colorbar
    title_txt=sprintf('Jacobi iteration number %d in staedy state',jcobi_iter);
    title(title_txt)
    xlabel('X-Axis')
    ylabel('Y-Axis')
    clabel(c,d)
    pause(0.03)
end
fprintf('No of iterations in Jacobi %d',jcobi_iter);
end

Function code for Gauss-Seidal(steady)

unction T=gs(strt_err,err_tol,x,y,nx,ny,dx,dy,T,k,Told)
gs_iter=1;
while(strt_err>err_tol)
    for i=2:nx-1
        for j=2:ny-1
            T(i,j)=(((T(i-1,j)+Told(i+1,j))/dx^2)+((T(i,j-1)+Told(i,j+1))/dy^2))/k;
        end
    end
    strt_err=max(max(abs(Told-T)));
    Told=T;
    gs_iter=gs_iter+1;
    figure(2)
    [c,d]=contourf(x,y,T);
    colormap(jet)
    colorbar
    title_txt=sprintf('Gauss Seidel iteration number %d in staedy state',gs_iter);
    title(title_txt)
    xlabel('X-Axis')
    ylabel('Y-Axis')
    clabel(c,d)
    pause(0.03)
end
fprintf('No of iterations in Gauss Seidel %d',gs_iter);
end

Code for Solving the equations(transient state)

close all 
clc

nx=41;
ny=nx;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
dt=0.01;
alpha=1.88e-5;
[xx,yy]=meshgrid(x,y);

T=298*ones(nx,ny);
T(1,:)=900;
T(end,:)=600;
T(:,1)=400;
T(:,end)=800;
T(1,1)=(400+900)/2;
T(1,ny)=(900+800)/2;
T(nx,1)=(400+600)/2;
T(nx,ny)=(600+800)/2;
Told=T;
Tprev=Told;
err_tol=1e-4;
strt_err=9e+9;
k1=alpha*(dt/(dx^2));
k2=alpha*(dt/(dy^2));
omega=2/(1+sin(pi*dx));
CFL=alpha*(dt/(dx^2));
T_jcob=Trns_jcob(strt_err,err_tol,xx,yy,nx,ny,T,k1,k2,Told,Tprev);
T_gs=Trns_gs(strt_err,err_tol,xx,yy,nx,ny,T,k1,k2,Told,Tprev);
T_sor=Trns_sor(strt_err,err_tol,xx,yy,nx,ny,T,k1,k2,Told,Tprev,omega);

Function code for Jacobi(transient)

function T=Trns_jcob(strt_err,err_tol,xx,yy,nx,ny,T,k1,k2,Told,Tprev)
jcobi_iter=1;
    for k=1:120000
        while(strt_err>err_tol)
            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(i,j)*term1+term2*H+term3*V;
                end
            end
            strt_err=max(max(abs(Told-T)));
            Told=T;
            jcobi_iter=jcobi_iter+1;
        end
        Tprev=T;
    end
     
            figure(1)
            contourf(xx,yy,T,'ShowText','on');
            colormap(jet)
            colorbar
            xlabel('X-Axis')
            ylabel('Y-Axis')
            colorbar
            title_txt=sprintf('Solution for 2D Heat conduction Equation in transient state with number of iterations %d',jcobi_iter);
            title(title_txt)
            fprintf('No of iterations in Jacobi %d',jcobi_iter);
end

Function code for SOR(transient)

function T=Trns_sor(strt_err,err_tol,xx,yy,nx,ny,T,k1,k2,Told,Tprev,omega)
sor_iter=1;
    for k=1:120000
        while(strt_err>err_tol)
            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)=omega*(Tprev(i,j)*term1+term2*H+term3*V)+(1-omega)*Tprev(i,j);
                end
            end
            strt_err=max(max(abs(Told-T)));
            Told=T;
            sor_iter=sor_iter+1;
        end
        Tprev=T;
    end
            figure(3)
            contourf(xx,yy,T,'ShowText','on');
            colormap(jet)
            colorbar
            xlabel('X-Axis')
            ylabel('Y-Axis')
            colorbar
            title_txt=sprintf('SOR iteration number %d in transient state',sor_iter);
            title(title_txt)
            fprintf('No of iterations in SOR %d',sor_iter);
end

Function code for Gauss-Seidal(transient)

function T=Trns_gs(strt_err,err_tol,xx,yy,nx,ny,T,k1,k2,Told,Tprev)
gs_iter=1;
    for k=1:12000
        while(strt_err>err_tol)
            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)+Told(i+1,j);
                    V=T(i,j-1)+Told(i,j+1);
                    T(i,j)=Tprev(i,j)*term1+term2*H+term3*V;
                end
            end
            strt_err=max(max(abs(Told-T)));
            Told=T;
            gs_iter=gs_iter+1;
        end
        Tprev=T;
    end
            figure(2)
            contourf(xx,yy,T,'ShowText','on');
            colormap(jet)
            colorbar
            xlabel('X-Axis')
            ylabel('Y-Axis')
            colorbar
            title_txt=sprintf('Gaussian iteration number %d in transient state',gs_iter);
            title(title_txt)
            fprintf('No of iterations in Gaussian %d',gs_iter);
end

Sample Output for Steady state

1.SOR

2.Gauss

3.Jacobi

Conclusion

  • Steady State method yields results very quickly with the fastest being SOR method
  • In transient case the time and no of iterations required is more
Your portfolio, next

Ship work like this
in your own course.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.