Numerical Code: Prediction of Temperature distribution within Aluminium Alloy chip during orthogonal turning process

clear all
clc

% Eng. Ammar Aziz   Zohaib Hassan     Ahmad Bilal      Ali Hassan  
%-------------------------------------------------------------------------
% Modeling of heat generation and Temperature distribution in machining
%-------------------------------------------------------------------------
   
   
%--------------------------------------------------------------------------
    %--------------------------------------------------------------
    % Heat generated in the Primary and Secondary deformation Zones
    %--------------------------------------------------------------

% The tool is standard C2 WC tool
% The chip is Al2024-T351.
hehe = input('Press any numeric key and then enter to start working=');

% Rake face angle of tool is

    alpha=input('rake face angle of tool(0 Radians is recommended for turning)=');
    rake_face_angle = num2str(alpha);
    display(rake_face_angle);

% Coefficient of friction between tool and chp 
  u = 0.57271;
  coefficient_of_friction = num2str(u);
  display(coefficient_of_friction);

% Normal friction angle in Radians
% Depends upon the coefficient of friction b/w tool and chip
    beta = atan(u);
    normal_friction_angle = num2str(beta);
    display(normal_friction_angle);
   
% Thickness of the cutted chip in meters is (found experimentally....)
t2 = 0.000333;
thickness_of_cutted_chip = num2str(t2);
display(thickness_of_cutted_chip);
% Contact length
cntct = 1.5*t2;

% Feed per revolution(in case of turning), in meters is given by
h = input('feed per revolution(0.000165 m/rev is standard)=');
feed_per_revolution = num2str(h);
display(feed_per_revolution);

% Chip thickness ratio is
    rat = h/t2;
    chip_thickness_ratio = num2str(rat);
    display(chip_thickness_ratio);

% Shear angle in the shear plane
% Depends upon the rake face angle and chip thickness ratio
    phi = atan((rat*cos(alpha))/(1-rat*sin(alpha)));
    shear_angle = num2str(phi);
    display(shear_angle);

% Width of the cutted chip in meters
w = 0.00254;
width_of_cutted_chip = num2str(w);
display(width_of_cutted_chip);

% Entering the shear strength of the chip283e6
S = 276e8;
shear_strength_of_chip = num2str(S);
display(shear_strength_of_chip);

 

% Cutting force in Newtons is given by (Merchant's circle approximation ..)
Fc = 573;
Cutting_Force = num2str(Fc);
display(Cutting_Force);

% Feed or Thrust force is Newtons is given by (Merchant's approximation..)
Ft = 329;
thrust_force = num2str(Ft);
display(thrust_force);

% Shear Force
Fs = Fc*cos(phi) - Ft*sin(phi);
shear_force = num2str(Fs);
display(shear_force);

% Shear area is given by
A_s = (h*w)/sin(phi);
shear_area = num2str(A_s);
display(shear_area);

% Shear stress within the primary zone is
tau = Fs/A_s;
shear_stress = num2str(tau);
display(shear_stress);

% The cutting velocity in meters per second is
Vw= 1.36;
cutting_velocity = num2str(Vw);
display(cutting_velocity);

% heat generation rate per unit depth of cut in primary zone :
Qs =(tau*h*Vw*cos(alpha))/(sin(phi)*cos(phi - alpha));
rate_of_heat_generation_in_primary_zone = num2str(Qs);
display(rate_of_heat_generation_in_primary_zone);
% heat generation rate for, length of cut as thickness of chip, in primary zone
Qso = Qs*t2;

% Frictional force acting between the rake face and chip is given by
Ff = Fc*sin(alpha) + Ft*cos(alpha);
frictional_force = num2str(Ff);
display(frictional_force);

% Rate of heat generation in the secondary zone is:
Qf =(tau*h*Vw*sin(beta))/(cos(phi + beta - alpha)*sin(phi - alpha));
rate_of_heat_generation_in_secondary_zone = num2str(Qf);
display(rate_of_heat_generation_in_secondary_zone);
% Rate of heat generation within contact length of secondary zone is:
Qfo = Qf*cntct;

%--------------------------------------------------------------------------

               

%--------------------------------------------------------------------------
                %------------------------------------
                %Average temperature rise calculation
                %------------------------------------
% Thermal conductivity of the workpiece in watts per meter per kelvin
  Kc=177;

% Mass density of chip in kilograms per cubic meter
  p=2785;

% specific heat capacity of chip in joules per kilogram per kelvin
  Cc=875;

% Thermal Diffusivity
Sa = Kc/(p*Cc);

% Thermal number
Rt = (h*Vw)/Sa;
   
% fraction of heat flowing into the tool
X = 0.5 - (0.35* (log10(Rt*tan(phi))));

% Average temperature rise of the (((chip))) per unit depth of cut due to shearing
delta_T = ((1-X)*Qs)/(p*Cc*h*Vw);
average_temp_rise = num2str(delta_T);
display(average_temp_rise);
%This average temperature rise on the shear plane is used as a boundary
%condition at the upper left corner of the chip

 

%-------------------------------------------------------------------------
                   
 
 
 
%--------------------------------------------------------------------------                   
                    %------------------
                    %Mesh generation
                    %------------------

% Chip velocity, for unit depth of cut, is given as
Vc = Qf/Ff;
chip_velocity = num2str(Vc);
display(chip_velocity);

 

% No of grids in x-direction
nbx = input('please enter no. of grids in x-direction (4 is default)=');
%No of grids in y-direction
nby = input('please enter no. of grids in y-direction (4 is default)=');

% No of nodes in x-direction
nbnx = nbx + 1;
% No of nodes in y-direction
nbny = nby + 1;

% proportion of frictional heat entering into the tool
Tio = 0.03;

% Grid spacing in x-direction
grsx = cntct/nbx;

% Grid spacing in y-direction
grsy = t2/nby;

% No. of nodes in the mesh
nod = nbnx*nbny;

% heat flow rate in chip's primary control zone due to primary zone heat source is
Qrs = ((1-Tio)*Qs*grsy)/t2;
heat_flow_in_chip_node_due_to_shear_source=num2str(Qrs);
display(heat_flow_in_chip_node_due_to_shear_source);

% heat generation rate for a  secondary node in chip due to frictional heat source
Qrf = ((1 - Tio)*Qf*grsx)/cntct;
heat_flow_in_chip_due_to_frictional_source=num2str(Qrf);
display(heat_flow_in_chip_due_to_frictional_source);

 

% Heat generation array
hga = [];
for i= 1:((sqrt(nod)));
    hga(i,sqrt(nod)) = Qrf;
end

for j= 1:((sqrt(nod))-1)
    for k = 1:(sqrt(nod)-1)
    hga(k,j) = 0;
    end
end

for m=1:(sqrt(nod)-1)
    hga(sqrt(nod),m) = Qrs;
end

hga(sqrt(nod),sqrt(nod)) = Qrf + Qrs;
display(hga);
   

%temperature distribution within the chip

Am0 = (p*Cc*Vc)/(Kc);
Am1 =  Am0/grsx;
Am2 = - (2/(grsx^2)) - (2/(grsy^2)) + Am1;
Am3 = (1/(grsx^2)) - (Am1);
Am4 = (grsx^2)*(grsy^2);
Am5 = Am4 * Am2;
Am6 = Am4 * Am3;
Am7 = grsy^2;
Am8 = grsx^2;
Am9 = Qrf/Kc;
A_1 = Am6/Am5;
A_2 = Am7/Am5;
A_3 = Am8/Am5;
A_4 = Am9/Am5;
A1 = -A_1;
A2 = -A_2;
A3 = -A_3;
A4 = -A_4;

for M = 1: (nbx - 1)
    xM = M*grsx;
end

for N = 1: (nby - 1)
    xN = N*grsy;
end

tempi = [];

% Boundary Condition
tempi(sqrt(nod),1) = delta_T;

 

% Data is imported from Solidworks 2012 by applying the boundary conditions
tempi = [344.1 345 346.1 347.3 350.3; 343 344.2 345.4 347.8 350.3; 339 336.9 343.9 348 351.2; 329.9 335.6 341.3 354.4 356.1; 234.5 343.1 352.7 362.4 370.6];
display(tempi);

% matrix of coefficients is given by
coeffik = (inv(tempi))*hga;
display(coeffik);

% Temperature distribution within the chip is given by multiplying inverse
% of coefficient matrix with heat generation array.

% plotting temperature
[x,y]=meshgrid(1:sqrt(nod),1:sqrt(nod));
surf(x,y,tempi)

 

%--------------------------------------------------------------------------
                    %THE END
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------