Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
Aim :- Solving Steady State 2D heat conduction equation.
Driver script for program.
%Script for solving 2d steady state heat conduction equationclear allclose allclcnx= 50;ny = nx;x=linspace(0,1,nx);y=linspace(0,1,ny);dx= x(2)-x(1);dy= y(2)-y(1);u = 400*ones(nx,ny);u(2:end-1,1)=400; u(1,2:end-1)=600; u(2:end-1,ny)=800;u(nx,2:end-1)=900;u(1,1) = 500;u(end,1) = 650;u(1,end)= 700;u(end,end)=850;tol = 1e-4;error = 9e9;iterative_solver = 1;k=(2*(dx^2+dy^2))/(dx^2*dy^2);if iterative_solver==1;u = jacobi_iter(u,x,y,nx,ny,dx,dy,tol,error,k);elseif iterative_solver==2;u = gauss_iter(u,x,y,nx,ny,dx,dy,tol,error,k);elseif iterative_solver==3;u = sor_iter(u,x,y,nx,ny,dx,dy,tol,error,k);endFollowing are the functions for Jacobi, Gauss-Seidel, and SOR method.
function u = jacobi_iter(u,x,y,nx,ny,dx,dy,tol,error,k)[X, Y] = meshgrid(x,y);uold=u;tic jacobi_iterr = 1; while(error>tol) for i = 2:nx-1 for j = 2:ny-1 u_1 = (uold(i-1,j)+uold(i+1,j))/(k*(dx^2)); u_2 = (uold(i,j-1)+uold(i,j+1))/(k*(dy^2)); u(i,j) = u_1+u_2; end end error = max(max(abs(uold-u))); uold = u; jacobi_iterr = jacobi_iterr+1; end time_required=toc;figure(1)C = contourf(X,Y,u);set(gca,'YDIR','reverse');colormap(jet);clabel(C,'FontSize',14);title(sprintf('Plot for steady state heat conduction equation using Jacobi method \n No. of Iterations: %d',jacobi_iterr));xlabel('X')ylabel('Y')hold on;saveas(gca,'2D steady state heat conduction Jacobi method.png');fprintf('Number of iterations done are %d \n',jacobi_iterr);fprintf('Time taken for solving equation is %f \n',time_required);fprintf('The error obtained is %f\n',error);endfunction u = gauss_iter(u,x,y,nx,ny,dx,dy,tol,error,k) uold=u; [X, Y] = meshgrid(x,y); gs_iter = 1; while(error>tol) for i = 2:nx-1 for j = 2:ny-1 u_1 = (u(i-1,j)+u(i+1,j))/(k*(dx^2)); u_2 = (u(i,j-1)+u(i,j+1))/(k*(dy^2)); u(i,j) = u_1+u_2; end end error =max(max(abs(uold-u))); uold = u; gs_iter = gs_iter+1; end solving_time=toc;figure(1)C = contourf(X,Y,u);set(gca,'YDIR','reverse');colormap(jet);clabel(C,'FontSize',14);title(sprintf('Plot for steady state heat conduction equation using Gauss Sidel method. \n No. of Iterations: %d',gs_iter));xlabel('X')ylabel('Y')hold on;saveas(gca,'2D steady state heat conduction Gauss Sidel method.png');fprintf('Number of iterations done are %d \n',gs_iter);fprintf('Time taken for solving equation is %f \n',solving_time);fprintf('The error obtained is %f\n',error);endfunction u = sor_iter(u,x,y,nx,ny,dx,dy,tol,error,k) % relaxation factor calculationz=cos(pi/nx)+cos(pi/ny);omega = 4/(2+sqrt(4-z^2));uold = u;[X Y] = meshgrid(x,y); sor_iter = 1; while(error>tol) for i = 2:nx-1 for j = 2:ny-1 u1 = (u(i-1,j)+uold(i+1,j))/(k*(dx^2)); u2 = (u(i,j-1)+uold(i,j+1))/(k*(dy^2)); u(i,j) = omega*(u1+u2)-(omega-1)*uold(i,j); end end error = max(max(abs(uold-u))); uold = u; sor_iter = sor_iter+1; end solving_time=toc;figure(1)C = contourf(X,Y,u);set(gca,'YDIR','reverse');colormap(jet);clabel(C,'FontSize',14);title(sprintf('Plot for steady state heat conduction equation using SOR method. \n No. of Iterations: %d\n',sor_iter));xlabel('X')ylabel('Y')hold on;saveas(gca,'2D steady state heat conduction SOR method.png');fprintf('Relaxation factor considered is %d \n',omega);fprintf('Number of iterations done are %d \n',sor_iter);fprintf('Time taken for solving equation is %f \n',solving_time);fprintf('The error obtained is %f\n',error);Results Obtained are as follows.
Following observations were made.
| Method | No. of iterations | Time required (seconds) | Error value |
| Jacobi | 4437 | 1.238763 | 0.000100 |
| Gauss Seidel | 2397 | 0.586370 | 0.000100 |
| SOR | 142 | 0.043705 | 0.000097 |
Conclusion :-
So we can conclude that SOR method gives fastest solution with minimum number of iterations and less error value as compared to others.
Solving Transient state heat conduction equation :-
Driver script for the program is as follows.
%Script for solving 2d transient state heat conduction equation by explicit%and implicit method.clear allclose allclcnx= 50;ny = nx;x=linspace(0,1,nx);y=linspace(0,1,ny);dx=x(2)-x(1);dy= y(2)-y(1);t_start = 0;dt = 1e-3;t_end =3;nt =[t_start:dt:t_end];u = 400*ones(nx,ny);u(2:end-1,1)=400; u(1,2:end-1)=600; u(2:end-1,ny)=800;u(nx,2:end-1)=900;u(1,1) = 500;u(end,1) = 650;u(1,end)= 700;u(end,end)=850;tol = 1e-4;error = 9e9;alpha=0.05; %assumedgk1=(alpha*dt)/(dx^2);k2=(alpha*dt)/(dy^2);uold = u;u_prev_dt = uold;iterative_solver = 1; % use 1,2,3 for explicit methods and 4,5,6 for implicit methods Jacobi, Gauss Sidel and SOR respectively.k=(2*(dx^2+dy^2))/(dx^2*dy^2);if iterative_solver==1;u = jacobi_iter_explicit(u,x,y,nx,ny,tol,error,k1,k2,nt);elseif iterative_solver==2;u = gauss_iter_explicit(u,x,y,nx,ny,tol,error,k1,k2,nt);ielseif iterative_solver==3;u = sor_iter_explicit(u,x,y,nx,ny,tol,error,k1,k2,nt);elseif iterative_solver==4;u = jacobi_2d_transient_implicit(u,x,y,nx,ny,tol,error,k1,k2,nt);elseif iterative_solver==5;u = gauss_2d_transient_implicit(u,x,y,nx,ny,tol,error,k1,k2,nt);elseif iterative_solver==6;u = sor_2d_transient_implicit(u,x,y,nx,ny,tol,error,k1,k2,nt);endFollowing are the functions for Jacobi, Gauss Seidel and SOR method.
Solving explicitly
function u = jacobi_iter_explicit(u,x,y,nx,ny,tol,error,k1,k2,nt)%jacobi method[X, Y] = meshgrid(x,y);uold=u;ticfor z=length(nt); jacobi_iterr = 1; while(error>tol) for i = 2:nx-1 for j = 2:ny-1 term1=uold(i,j); term2=k1*(uold(i-1,j)-2*uold(i,j)+uold(i+1,j)); term3=k2*(uold(i,j-1)-2*uold(i,j)+uold(i,j+1)); u(i,j)=term1+term2+term3; end end error = max(max(abs(uold-u))); uold = u; jacobi_iterr = jacobi_iterr+1; endend %Saving the plot and printing the results time_required=toc;figure(1)C = contourf(X,Y,u);set(gca,'YDIR','reverse');colormap(jet);clabel(C,'FontSize',14);title(sprintf('Plot for solving Transient state heat conduction equation explicitly\n using Jacobi method \n No. of Iterations: %d\n',jacobi_iterr));xlabel('X');ylabel('Y');hold on;saveas(gca,'2D Transient state heat conduction Jacobi method EXPLICITLY.png');fprintf('Number of iterations done are %d \n',jacobi_iterr);fprintf('Time taken for solving equation is %f \n',time_required);fprintf('The error obtained is %f\n',error);end function u = gauss_iter_explicit(u,x,y,nx,ny,tol,error,k1,k2,nt)%Gauss Sidel method[X, Y] = meshgrid(x,y);uold=u;ticfor z=length(nt); gauss_iter = 1; while(error>tol) for i = 2:nx-1 for j = 2:ny-1 term1=uold(i,j); term2=k1*(u(i-1,j)-2*u(i,j)+u(i+1,j)); term3=k2*(u(i,j-1)-2*u(i,j)+u(i,j+1)); u(i,j)=term1+term2+term3; end end error = max(max(abs(uold-u))); uold = u; gauss_iter = gauss_iter+1; endend %Saving the plot and printing the results time_required=toc;figure(1)C = contourf(X,Y,u);set(gca,'YDIR','reverse');colormap(jet);clabel(C,'FontSize',14);title(sprintf('Plot for solving Transient state heat conduction equation explicitly\n using Gauss Sidel method. \n No. of Iterations: %d\n',gauss_iter));xlabel('X')ylabel('Y')hold on;saveas(gca,'2D Transient state heat conduction Gauss Sidel method EXPLICITLY.png');fprintf('Number of iterations done are %d \n',gauss_iter);fprintf('Time taken for solving equation is %f \n',time_required);fprintf('The error obtained is %f\n',error);endfunction u = sor_iter_explicit(u,x,y,nx,ny,tol,error,k1,k2,nt)%SOR method % relaxation factor calculationz=cos(pi/nx)+cos(pi/ny);omega = 4/(2+sqrt(4-z^2));uold = u;[X,Y] = meshgrid(x,y);ticfor z = length(nt); sor_iter = 1; while(error>tol) for i = 2:nx-1; for j = 2:ny-1; term1=uold(i,j); term2=k1*(u(i-1,j)-2*uold(i,j)+uold(i+1,j)); term3=k2*(u(i,j-1)-2*uold(i,j)+uold(i,j+1)); u(i,j)=omega*(term1+term2+term3)-(omega-1)*term1; end end error = max(max(abs(uold-u))); uold = u; sor_iter = sor_iter+1; end end %Saving the plot and printing the results. time_required=toc;figure(1)C = contourf(X,Y,u);set(gca,'YDIR','reverse');colormap(jet);clabel(C,'FontSize',14);title(sprintf('Plot for transient state heat conduction equation using SOR method. \n No. of Iterations: %d \n',sor_iter));xlabel('X')ylabel('Y')hold on;saveas(gca,'2D transient state heat conduction SOR method EXPLICITLY.png');fprintf('Relaxation factor considered is %d \n',omega);fprintf('Number of iterations done are %d \n',sor_iter);fprintf('Time taken for solving equation is %f \n',time_required);fprintf('The error obtained is %f\n',error);Solving Implicitly
function u = jacobi_2_transient_implicit(u,x,y,nx,ny,tol,error,k1,k2,nt)uold = u;u_prev_dt = uold;term1=(1+(2*k1)+(2*k2))^-1;term2=k1*term1;term3=k2*term1;ticfor z = 1:length(nt) %Time loop jacobi_iter = 1; error = 9e9; while (error>tol) % Convergence loop % Nodal loop for i = 2:nx-1 for j = 2:ny-1 H=(uold(i-1,j)+uold(i+1,j)); V=(uold(i,j-1)+uold(i,j+1)); u(i,j) =((u_prev_dt(i,j)*term1)+(H*term2)+(V*term3)); end end error=max(max(abs(uold-u))); uold = u; jacobi_iter = jacobi_iter+1; end u_prev_dt = uold; endtime_required=toc;figure(1)[X, Y] = meshgrid(x,y);C = contourf(X,Y,u);set(gca,'YDIR','reverse');colormap(jet);clabel(C,'FontSize',14);title(sprintf('Plot for solving Transient state heat conduction equation implicitly.\n using Jacobi method. \n No. of Iterations: %d\n',jacobi_iter));xlabel('X');ylabel('Y');hold on;saveas(gca,'2D Transient state heat conduction Jacobi method IMPLICITLY1.png');fprintf('Number of iterations done are %d \n',jacobi_iter);fprintf('Time taken for solving equation is %f \n',time_required);fprintf('The error obtained is %f\n',error); function u = gauss_2_transient_implicit(u,x,y,nx,ny,tol,error,k1,k2,nt)term1=(1+(2*k1)+(2*k2))^-1;term2=k1*term1;term3=k2*term1;uold=u;u_prev_dt = uold;ticfor z = 1:length(nt) %Time loop gs_iter = 1; error = 9e9; while (error>tol) % Convergence loop % Nodal loop for i = 2:nx-1 for j = 2:ny-1 H=(u(i-1,j)+u(i+1,j)); V=(u(i,j-1)+u(i,j+1)); u(i,j) =((u_prev_dt(i,j)*term1)+(H*term2)+(V*term3)); end end error=max(max(abs(uold-u))); uold = u; gs_iter = gs_iter+1; end u_prev_dt = uold; endtime_required=toc;figure(1)[X, Y] = meshgrid(x,y);C = contourf(X,Y,u);set(gca,'YDIR','reverse');colormap(jet);clabel(C,'FontSize',14);title(sprintf('Plot for solving Transient state heat conduction equation implicitly. By using gauss method. \n No. of Iterations: %d\n',gs_iter));xlabel('X');ylabel('Y');hold on;saveas(gca,'2D Transient state heat conduction Gauss sidel method IMPLICITLY.png');fprintf('Number of iterations done are %d \n',gs_iter);fprintf('Time taken for solving equation is %f \n',time_required);fprintf('The error obtained is %f\n',error); function u = sor_2_transient_implicit(u,x,y,nx,ny,tol,error,k1,k2,nt)uold = u;u_prev_dt = uold;z=cos(pi/nx)+cos(pi/ny);omega = 4/(2+sqrt(4-z^2));term1=(1+(2*k1)+(2*k2))^-1;term2=k1*term1;term3=k2*term1;ticfor z = 1:length(nt) %Time loop sor_iter = 1; error = 9e9; while (error>tol) % Convergence loop % Nodal loop for i = 2:nx-1 for j = 2:ny-1 H=(u(i-1,j)+u(i+1,j)); V=(u(i,j-1)+u(i,j+1)); u(i,j) =omega*((u_prev_dt(i,j)*term1)+(H*term2)+(V*term3))-(omega-1)*u_prev_dt(i,j); end end error=max(max(abs(uold-u))); uold = u; sor_iter = sor_iter+1; end u_prev_dt = uold; end time_required=toc; figure(1)[X, Y] = meshgrid(x,y);C = contourf(X,Y,u);set(gca,'YDIR','reverse');colormap(jet);clabel(C,'FontSize',14);title(sprintf('Plot for solving Transient state heat conduction implicitly using SOR method. \n No. of Iterations: %d\n',sor_iter));xlabel('X');ylabel('Y');hold on;saveas(gca,'2D Transient state heat conduction SOR method IMPLICITLY1.png');fprintf('Number of iterations done are %d \n',sor_iter);fprintf('Time taken for solving equation is %f \n',time_required);fprintf('The error obtained is %f\n',error); Following are the results obtained.
Following observations were made.
FOR EXPLICIT METHOD
| Method | No. of iterations | Time required (seconds) |
Error value |
| Jacobi | 8499 | 2.149332 | 0.000100 |
| Gauss-Seidel | 6681 | 1.689989 | 0.000100 |
| SOR | 2849 | 0.987259 | 0.000100 |
FOR IMPLICIT METHOD
| Method | No. of iterations | Time required (seconds) |
Error value |
| Jacobi | 7 | 6.523987 | 0.000055 |
| Gauss-Seidel | 6 | 4.775224 | 0.000026 |
| SOR | 2 | 6.207440 | 0.000083 |