Week 5 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
Objectives of the project:-
1) Obtain temperature distribution for a thin square plate by solving 2D conduction equation.
2) Solve for both- steady state as well as the transient state of the 2D heat conduction equation and compare the results.
3) Solve the transient state equation by explicit and implicit methods.
4) The implicit method employs iterative solvers like Jacobi, Gauss-Seidel and Successive over-relaxation. Use it in your code and compare these methods.
Theory:-
The Conduction equation is given as,
where alpha is the thermal diffusivity of material (m^2/sec). We assume that thermal diffusivity is not affected by temperature, thereby the resulting differential equation is linear in nature.
For a thin plate, the conduction equation is written as,
as the dimension in the z-direction is very small as the plate is thin. This form is the Transient form of the 2D conduction equation.
In the Steady-state form of the conduction equation, temperature does not change with time. Hence, the time derivative term is zero. The Steady-state conduction equation is given as,
(1) Explicit:
The Forward Euler procedure gives the explicit form of the governing equation.
By definition, in an explicit approach, each difference equation contains only one unknown and can be solved explicitly for this one unknown in a straightforward manner.
The first step to solve a differential equation by computational means is to discretize the equation. The time derivative is discretized by using the first-order forward difference approximation, whereas the spatial terms are discretized by using the second-order central difference approximation.
The value of T at time level 'n+1' is thus calculated at all grid points by changing x-coordinate and y-coordinate in the grid. Thus once all values of T are known at time level 'n+1', we then calculate T at next time-step i.e. at the time 'n+2'. This is carried out until the end of simulation time or until the solution converges.
This method is called as Time marching solution, as we increment the time step one at a time.
The code is as follows:-
% Transient state Heat conduction clear allclose allclc%Domain conditiosnL = 1;X_L = L; %length in metersY_L = L; %width in meters%time stepT = 500;dt = 1e-2;alpha = 1e-3;%Grid Generationn = 11;Grid = linspace(0,L,n);x = Grid;dx = x(2)-x(1);y = Grid;dy = y(2)-y(1);error = 1e7;tol = 1e-4;%inputs Boundry conditionsu = 300*ones(n,n);u(:,1) = 400;u(1,:) = 600;u(:,end) = 800;u(end,:) = 900;uold = u;k1 = ((alpha*dt)/dx^2);k2 = ((alpha*dt)/dy^2);ct = 1;%time loop and node loopfor i = 1:T while (error > tol) %% Jacobi Iterative method for j = 2:length(x)-1 for k = 2:length(y)-1 term_1 = (uold((j+1),k) - (2*uold(j,k)) + uold((j-1),k))/dx^2; term_2 = (uold(j,(k+1)) - (2*uold(j,k)) + uold(j,(k-1)))/dy^2; u(j,k) = uold(j,k) + ((k1*term_1) + (k2*term_2)); end end uold = u; [C,h] =contourf(x,y,u); set(gca,'xTick',(0:dx:1)) % Setting x tick intervals as step size in i direction set(gca,'yTick',(0:dy:1)) % Setting y tick intervals as step size in j direcion clabel(C, h,'FontSize',8,'FontWeight','bold','color','k') colorbar; colormap(jet); title(sprintf('Temperature profile@ %d timestep is',ct)) pause(0.03) Q(ct) = getframe(gcf); ct = ct+1;endendResults:-
(2) Implicit:
The Backward Euler method or the Crank-Nicolson method gives the implicit form of the equation.
In Implicit approach, each difference equation contains more than one unknown. Hence it cannot be solved directly. It results in a system of algebraic equations for which the unknowns can be solved simultaneously.
Crank-Nicolson form is used mostly to solve the parabolic governing equations in CFD. In this project, backward Euler method is used for discretization.
The time derivative is discretized by using the first-order forward difference approximation, whereas the spatial terms are discretized by using the second-order central difference approximation.
Iterative Methods:
The equation is written in the given form,
where
= ,,,,
;
(a) Jacobi:
In this method, the unknowns on the right-hand side are assumed to have an initial value, i.e. T at 'n+1' at all nodes is assumed. Using the equation, then the LHS is calculated. On the next value of i,j similar method is applied to calculate LHS. Thus in the first iteration, value of T for time 'n+1' is calculated at all nodes. The difference between T at 'n+1' which is obtained and T at 'n+1' which was assumed previously is taken. If it is greater than a specified tolerance value, iterations are carried out to calculate T at 'n+1'. But in the second iteration, T at 'n+1' is assumed as the value of T which was obtained in the previous iteration.
Once the value of T is found out iteratively whose error is less than the tolerance value, the value of T is calculated at the next time step 'n+2' using the same iterative method.
(b) Gauss-Seidel:
It is similar to the Jacobi method, but the only difference is- For the first iteration when the value of T at a particular node is calculated, this updated value at that node is used to calculate the value at the next node. But, the error term is calculated similar to Jacobi method, i.e. T at 'n+1' calculated minus the T at 'n+1' previously assumed.
Once the value of T is found out iteratively whose error is less than the tolerance value, the value of T is calculated at the next time step 'n+2' using the same iterative method.
(c) Successive over-relaxation:
This method can be applied to both Jacobi and Gauss-Seidel. Here a weighting factor 'w' is considered, due to which the increment in values of T is more rigorous, thereby the number of iterations for convergence of solution is reduced.
The equation for SOR is modified as,
The code for the Transient Implicit Computation is as follows:-
%Implicit Unsteady solver - 2D heat conduction problems% MATLAB Program to solve the unsteady 2D Heat Conduction equation Implictly:clear all;close all;clc;%% Initializing Variables:Nx = 11;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 TemperatureT(end,:) = 600; % Top TemperatureT(:,1) = 400; % Left TemperatureT(:,end) = 800; % Right Temperature% Edge refiningT(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:2000 %% 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;endexecution_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));Results:-
For the Jacobi Method:
For the Gauss-Seidel Method :-
For the SOR Method:-
Solving Steady-state equation:
Only the spatial terms are present in steady-state conduction equation. The time derivative is zero. Practically, it is the condition where the temperature distribution will no longer be affected by time. The steady-state equation is given as,
The equation is discretized and solved by the code as follows:-
% ========================== 2D Steady State Heat Conduction Equation - Explicit Schemes ====================================== %clear allclose allclc I = input('Enter the desired number: ') % For loop to use all the point iterative techniques % Starting the timer to calculate the computational time %% Problem Setup Lx = 1; % Lenghth of x grid Ly = Lx; % Length of y grid nx = 11; % Number of nodes along x grid ny = nx; % Number of nodes along y grid dx = Lx/(nx-1); % Step size for spatial x co-ordinate dy = Ly/(ny-1); % Step size for spatial y - co-ordinate k = 2*(1/(dx^2)+1/(dy^2)); tol = 1e-4; % Error criteria error = 1e4; % Initialising error l = 0; % Initialising iteration w = 1.7; %% Creating 1D and 2D grids x = linspace(0, Lx, nx); % Creating 1D x grid y = linspace(0, Ly, ny); % Creating 1D y grid [X Y] = meshgrid(x, y); % Creating 2D x-y grid %% Initialising temprature array and assigning boundary conditions T = ones(nx,ny); % Initialising temprature array %Assigning BCs T(1,:) = 900; % T at bottom phase T(end,:) = 600; % T at top phase T(:,1) = 400; % T at left phase T(:,end) = 800; % T at right phase % Creating a copy of Temperature array after assigning BCs Told = T; while error > tol % While loop runs untill error > tolerance value for i = 2:nx-1 % For loop runs for x values from 2 to last node - 1 for j = 2:ny-1 % For loop runs for y values from 2 to last node - 1 %% Iterative solvers if I == 1 % Jacobi figure(1) T(i,j) = (1/(k*(dx^2)))*(Told(i-1,j) + Told(i+1,j) + Told(i,j-1) + Told(i,j+1)); elseif I == 2 % Gauss Seidel figure(2) T(i,j) = (1/(k*(dx^2)))*(T(i-1,j) + T(i+1,j) + T(i,j-1) + T(i,j+1)); elseif I ==3 % Successive over relaxation figure(3) T(i,j) = w*((1/(k*(dx^2)))*(T(i-1,j) + T(i+1,j) + T(i,j-1) + T(i,j+1))) + (1-w)*Told(i,j); end end % Ending for loop for y values end % Ending for loop for x values error_T = max(abs(Told - T)); % Calculating error to decide wether while loop should stop or not error = max(error_T); % Calculating maximum error between absolute value of old array of T and new array of T Told = T; % Updating Temprature array for next iteration l = l + 1; % Incrementing iterations by 1. %Generating the contour map plotfilled with colours representing distribution of Temperature in the 1 x 1 domain [C, h] = contourf(X,Y,T); % Contourf function used set(gca,'XTick',(0:dx:1)) % Setting x tick intervals as step size in i direction set(gca,'YTick',(0:dy:1)) % Setting y tick intervals as step size in j direcion clabel(C, h,'FontSize',8,'FontWeight','bold','color','k') % Labeling the lines on the plot xlabel('Spatial X Co-ordinate','FontSize',12) % Labeling x axis ylabel('Spatial Y Co-ordinate','FontSize',12) % Labeling y axis colorbar title(sprintf('Temprature distribution on grid using GS Method at iteration no %d',l),'FontSize',15); % Giving title to the plot that changes with iteration number pause(0.003) % Pausing the plot for 0.003 seconds to observe the changes end % Ending while loop Results:-
CONCLUSIONS:
1. Transient state equation requires more computational time than the steady-state equation. Also when the time required for attaining steady state is large.
2. The time for convergence for an explicit method is very less as compared to the implicit method. But the explicit method seems to be unstable for CFL greater than 0.25 in this case. Whereas, the implicit method is conditionally unstable. If CFL number is greater than 0.25 then implicit method is preferred as we know that higher courant numbers can be used to get stable results for implicit method.
3. For Steady-state equation- Jacobi solver takes less time for convergence than the Gauss-Seidel solver. But, the number of iterations for Gauss-Seidel is less than Jacobi solver. But the most efficient solver amongst these three is the SOR solver. It has the least number of iterations(50 for my case) and the minimum convergence time.
4. For Transient state equation- Explicit method has the least convergence time. But it is only preferred when CFL<0.25, otherwise Implicit method is used. Jacobi has the least time but the most number of iterations. The SOR method is preferred over Gauss-seidel, as it takes less time for convergence and also has far lesser iterations than Gauss-Seidel.
N.B: Due to the lack of computational resources available to me, I have only used Nx=Ny=11..ie. a 10x10 element mesh.