guneysus
5/25/2013 - 5:34 PM

The problem is solving a heat equation for wing has circle geometry fixed to 200 C wall and isolated at x = L Thomas algoritm for 1D stea

The problem is solving a heat equation for wing has circle geometry fixed to 200 C wall and isolated at x = L

Thomas algoritm for 1D steady state heat equation with T0 = 200 C; and isolated wall at x=L dT/dx = 0;

and given heat equation in differential equation d2T/dx2 -hP/kA ( T - T_inf ) = 0

T_inf = 20 C

k = 200 W/mK D = 20 mm L = 150 mm h = 50 W/m^2K

Heat flow Q = -kA dT/dx

main.m
content_copyfile_download
function main()
clc;
k = 200 ; %W/mK
D = 20e-3 ; %m
L = 150e-3 ; %m
h = 50 ; %W/m^2K

T0 = 200;
Tinf = 20;

A = pi * D^2 / 4 ; %m^2
P = pi * D ;

C = -(h*P) / (k*A) ;

n = 100; % number of nodes
dx = L / n ;

b ( 1:n-1 ) = -2/dx^2 + C ;
b ( n ) = 1/dx ;

a ( 2:n-1 ) = 1/dx^2 ;
a ( n ) = -1/dx ;

c ( 1:n-1 ) = 1/dx^2 ;

d ( 1 ) = -T0/dx^2 + C*Tinf ;
d ( 2:n-1 ) = + C * Tinf ;
d ( n ) = 0 ;

T = [T0 TDMAsolver(a,b,c,d)];
x = 0:dx:L;

Q = k * A *  ( T(2) - T(1) ) / dx;
disp (['Heat flow from wing : ' num2str(Q) ' W']);
plot(x,T,'r--','LineWidth',2.0),grid on, title('1D Steady State Heat Flow from wing'),...
    xlabel('x (m)'),ylabel('T (K)');

    function x = TDMAsolver(a,b,c,d)
    %a, b, c are the column vectors for the compressed tridiagonal matrix, d is the right vector
    n = length(b); % n is the number of rows

    % Modify the first-row coefficients
    c(1) = c(1) / b(1);    % Division by zero risk.
    d(1) = d(1) / b(1);    % Division by zero would imply a singular matrix.

        for i = 2:n-1
            temp = b(i) - a(i) * c(i-1);
            c(i) = c(i) / temp;
            d(i) = (d(i) - a(i) * d(i-1))/temp;
        end

    d(n) = (d(n) - a(n) * d(n-1))/( b(n) - a(n) * c(n-1));

    % Now back substitute.
    x(n) = d(n);
        for i = n-1:-1:1
            x(i) = d(i) - c(i) * x(i + 1);
        end

    end

end