45 views (last 30 days)

Oonitejas Sahoo
on 10 Apr 2018

- %% its the mathematical approach by which i sloved the thermal problems
- clear all
- close all
- %Specify grid size
- Nx = 10;
- Ny = 10;
- %Specify boundary conditions
- Tbottom = 50
- Ttop = 150
- Tleft = 50
- Tright = 50
- % initialize coefficient matrix and constant vector with zeros
- A = zeros(Nx*Ny);
- C = zeros(Nx*Ny,1);
- % initial 'guess' for temperature distribution
- T(1:Nx*Ny,1) = 100;
- % Build coefficient matrix and constant vector
- % inner nodes
- for n = 2:(Ny-1)
- for m = 2:(Nx-1)
- i = (n-1)*Nx + m;
- A(i,i+Nx) = 1;
- A(i,i-Nx) = 1;
- A(i,i+1) = 1;
- A(i,i-1) = 1;
- A(i,i) = -4;
- end
- end
- % Edge nodes
- % bottom
- for m = 2:(Nx-1)
- %n = 1
- i = m;
- A(i,i+Nx) = 1;
- A(i,i+1) = 1;
- A(i,i-1) = 1;
- A(i,i) = -4;
- C(i) = -Tbottom;
- end
- %top:
- for m = 2:(Nx-1)
- % n = Ny
- i = (Ny-1)*Nx + m;
- A(i,i-Nx) = 1;
- A(i,i+1) = 1;
- A(i,i-1) = 1;
- A(i,i) = -4;
- C(i) = -Ttop;
- end
- %left:
- for n=2:(Ny-1)
- %m = 1
- i = (n-1)*Nx + 1;
- A(i,i+Nx) = 1;
- A(i,i+1) = 1;
- A(i,i-Nx) = 1;
- A(i,i) = -4;
- end
- %right:
- for n=2:(Ny-1)
- %m = Nx
- i = (n-1)*Nx + Nx;
- A(i,i+Nx) = 1;
- A(i,i-1) = 1;
- A(i,i-Nx) = 1;
- A(i,i) = -4;
- C(i) = -Tright;
- end
- % Corners
- %bottom left (i=1):
- i=1;
- A(i,Nx+i) = 1;
- A(i,2) = 1;
- A(i,1) = -4;
- C(i) = -(Tbottom + Tleft);
- %bottom right:
- i = Nx;
- A(i,i+Nx) = 1;
- A(i,i-1) = 1;
- A(i,i) = -4;
- C(i) = -(Tbottom + Tright);
- %top left:
- i = (Ny-1)*Nx + 1;
- A(i,i+1) = 1;
- A(i,i) = -4;
- A(i,i-Nx) = 1;
- C(i) = -(Ttop + Tleft);
- %top right:
- i = Nx*Ny;
- A(i,i-1) = 1;
- A(i,i) = -4;
- A(i,i-Nx) = 1;
- C(i) = -(Tright + Ttop);
- %Solve using Gauss-Seidel
- residual = 100;
- iterations = 0;
- while (residual > 0.0001) % The residual criterion is 0.0001 in this example
- % You can test different values
- iterations = iterations+1;
- %Transfer the previously computed temperatures to an array Told
- Told = T;
- %Update estimate of the temperature distribution
- for n=1:Ny
- for m=1:Nx
- i = (n-1)*Nx + m;
- Told(i) = T(i);
- end
- end
- % iterate through all of the equations
- for n=1:Ny
- for m=1:Nx
- i = (n-1)*Nx + m;
- %sum the terms based on updated temperatures
- sum1 = 0;
- for j=1:i-1
- sum1 = sum1 + A(i,j)*T(j);
- end
- %sum the terms based on temperatures not yet updated
- sum2 = 0;
- for j=i+1:Nx*Ny
- sum2 = sum2 + A(i,j)*Told(j);
- end
- % update the temperature for the current node
- T(i) = (1/A(i,i)) * (C(i) - sum1 - sum2);
- end
- end
- residual = max(T(i) - Told(i));
- end
- %compute residual
- deltaT = abs(T - Told);
- residual = max(deltaT);
- iterations; % report the number of iterations that were executed
- %Now transform T into 2-D network so it can be plotted.
- delta_x = 0.03/(Nx+1);
- delta_y = 0.03/(Ny+1);
- for n=1:Ny
- for m=1:Nx
- i = (n-1)*Nx + m;
- T2d(m,n) = T(i);
- x(m) = m*delta_x;
- y(n) = n*delta_y;
- end
- end
- T2d;
- surf(x,y,T2d)
- figure
- contour(x,y,T2d)

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!
## 0 Comments

Sign in to comment.