Heat Transfer

clc; clear all; close all;

%%%%%% Heat transfer across a thin plate with a hole %%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% Material Properties %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Aluminum
a     = 9.7e-5; % m^2/s Thermal Diffusivity
cp    = 0.91;   % kJ/kg*K Heat Capacity
p     = 2700;   % kg/m^3 Density




%%%%%% Input for plate geometry %%%%%%

Lx        = 0.13;        % length in the x-direction in meters
Ly        = 0.167;        % length in the y-direction in meters
thickness = 0.15;       % thickness of the plate
delta     = 0.0025;       % distance between grid points in meters
R         = 0.05;      % radius of the hole in the plate in meters
C         = [Lx/2;Ly/2]; % center of the hole
Mass      = Lx*Ly*thickness*p; % Mass of Plate

%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% Temperatures %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%

% Gasoline
Tc    = 530.372;    % degrees Kelvin (Temp of Combustion)
% Ambient Temperature
Ta    = 288.706;    % degrees Kelvin (Temp Outside)
% Effective Temperature of Coolant
Te    = 274.817;    % degrees Kelvin (Temp of the coolant)



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% Create plate geometry %%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

x = [0:delta:Lx];   % vector of all x values
y = [0:delta:Ly];   % vector of all y values


M = length(x);      % Number of points in x direction
N = length(y);      % Number of points in y direction

for i=1:M           % Loop over all x values
     for j=1:N      % Loop over all y values

         X(i,j)= x(i); % Gridded data for all x values
         Y(i,j)= y(j); % Gridded data for all y values
         Conditional(i,j) = 1; % all points exist by default

         %%% Initial conditions %%%

         Q0(i,j) = Ta*Mass*cp; %Set heat equal to zero

         if j == N
             Q0(i,j) = Te*Mass*cp % Temperature of cooled side
         end






         %%% Create hole %%%
         D = sqrt( (C(1) - X(i,j))^2 + (C(2) - Y(i,j))^2 ); % distance from current point to center of the hole
         if D < R
             Conditional(i,j) = 0; % if point is in the hole, point does not exist
         end
    end
end
for i=2:M-1 %Loop through all x-points not on the edge of plate
    for j=2:N-1 %Loop through all y-points not on the edge of plate
    LeftN   = Conditional(i-1,j);
    RightN  = Conditional(i+1,j);
    TopN    = Conditional(i,j+1);
    BottomN = Conditional(i,j-1);

    %If any neighboring values DNE and the point exists, they are on the edge of the hole
        if LeftN ==0 || RightN == 0 || TopN == 0 || BottomN == 0
            if Conditional(i,j) == 1
             Q0(i,j) = Tc*Mass*cp; %Edge of hole is hot, like a piston would be
            end
        end
    end
end

%%%%% Extend Plate %%%%%

% Add length of plate in x direction plus delta to current x-values
%X1 = X + Lx + delta;

% Paste X and X1 together vertically where X1 is pasted just below X
%X = vertcat(X,X1);
% Paste two copies of Y vertically
%Y = vertcat(Y,Y);
% Paste two copies of Conditional vertically
%Conditional = vertcat(Conditional,Conditional);
% Paste two copies of Q0 vertically and set the initial heat of the joining
% point to zero
%Q0 = vertcat(Q0,Q0);
%Q0(M,:) = 0;
%Q0(M+1,:) = 0;

% Redefine M and N
%M = length(X(:,1)); % M equals the number of rows in X
%N = length(Y(1,:)); % N equals the number of columns in Y

%%%%%% Solve Heat Transfer Equation %%%%%%

Q = Q0;    % Let current Q equal initial Q to start with
t = 0;     % Time starts at zero seconds
dt = 8e-5; % Change in time in seconds

Tavg0 = 0;
dTavg = 10;
eps = 1e-12;
while dTavg > eps

while t < 10
   for i=2:M-1     % Loop over all x values not on boundary
       for j=2:N-1 % Loop over all y values not on boundary

           if Conditional(i,j) == 1 % If point exists, solve heat equation

                LeftN   = Conditional(i-1,j);
                RightN  = Conditional(i+1,j);
                TopN    = Conditional(i,j+1);
                BottomN = Conditional(i,j-1);

               if LeftN == 1 && RightN == 1 && TopN == 1 && BottomN == 1
               Q(i,j) = ( (Q0(i+1,j) - ...
                   2*Q0(i,j) + Q0(i-1,j))/(0.02^2) + ...
                   (Q0(i,j+1) - 2*Q0(i,j) + Q0(i,j-1))/(0.02^2))*dt + Q0(i,j);
               end
           end
       end
   end
t = t + dt; % Update time
Q0 = Q;     % Make previous heat values equal current heat values for next loop


figure(1)
cla
hold on
surf(X,Y,Q,'EdgeColor','none')
contour3(X,Y,Q,20,'k')
xlabel('X')
ylabel('Y')
shading interp
colormap(jet)
view(0,90)
pause(1e-16)

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%% Calculate Average Temperature Of Engine Head %%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

T = Q./(Mass*cp);

Tsum    = 0;

Npoints = 0;


for i=1:M

    for j=1:N

        if Conditional(i,j) == 1



            Tsum    = T(i,j) + Tsum;

            Npoints = Npoints + 1;



        end

    end

end

Tavg = Tsum/Npoints;
dTavg = Tavg - Tavg0;
Tavg0 = Tavg;

%tcooling = [200; 250; 274.81; 300; 350];

%taverage = [374.64; 384.81; 389.86; 394.98; 405.15];


end