
clc
clear all
close all

%% Geometric and material parameters
r_s= 0.005;         %Sphere radius [m]
d_s= 2*r_s;         %Sphere diameter [m]
E_s= 390*10^9;      %[Pa] Sphere elastic modulus
ni_s= 0.26;         %Sphere Poisson ration

E_d= 218*10^9;      %[Pa] Disc Elastic modulus
ni_d= 0.30;         %Disc Poisson ration

F= 10;              %[N] Applied vertical load
Er= 0.5*( (1-ni_s^2)/E_s + (1-ni_d^2)/E_d ); %Reduced elastic modulus of material pair

l= 0.016;           %[m] stroke length -> wear increment
L= 1000;            %[m] test duration - total sliding distance
N= round(L/l);      %overall number of duration spans -> number of cycles doubled

%% Experimental wear law
% Archard coefficients from excel analysis of wear curves
k1= 8*10^(-15); % [m3/(m*N)] dimensional wear rate k = V/(s*F)- wear-in
L1= 23; %[m] end of wear-in
k2= 3*10^(-14); % [m3/(m*N)] dimensional wear rate k = V/(s*F)- wear phase I
L2= 175; %[m] end of wear phase I
k3= 5*10^(-14); % [m3/(m*N)] dimensional wear rate k = V/(s*F)- wear phase II

k= zeros(N,1); %initializing the dimensional wear coefficient vector
k( 1:round(L1/l) )= k1;
k( (round(L1/l)+1) : round(L2/l) )= k2;
k( (round(L2/l)+1) : end )= k3;

% Experimental data import
opts = delimitedTextImportOptions("NumVariables", 2);
opts.DataLines = [1, Inf]; % Specify range and delimiter
opts.Delimiter = "\t";
opts.VariableNames = ["VarName1", "VarName2"]; % Specify column names and types
opts.VariableTypes = ["double", "double"];
opts.ExtraColumnsRule = "ignore"; % Specify file level properties
opts.EmptyLineRule = "read";
opts = setvaropts(opts, ["VarName1", "VarName2"], "DecimalSeparator", ",");% Specify variable properties

% Import the data of wear track cross-sectional area
WearTrackArea = readtable('Wear track area.txt', opts);     %File must be in the same folder as the script
% Convert to numeric matrix
WearTrackArea = table2array(WearTrackArea);

% Import the data of wear track depth
WearTrackDepth = readtable('Wear track depth.txt', opts);   %File must be in the same folder as the script
% Convert to numeric matrix
WearTrackDepth = table2array(WearTrackDepth);

% Clear temporary variables
clear opts

% Table's data import, first have to change load options
opts = delimitedTextImportOptions("NumVariables", 5);
opts.DataLines = [1, Inf]; % Specify range and delimiter
opts.Delimiter = "\t";
opts.VariableNames = ["VarName1", "VarName2", "VarName3", "VarName4", "VarName5"]; % Specify column names and types
opts.VariableTypes = ["double", "double", "double", "double", "double"];
opts.ExtraColumnsRule = "ignore"; % Specify file level properties
opts.EmptyLineRule = "read";
opts = setvaropts(opts, ["VarName1", "VarName2", "VarName3", "VarName4", "VarName5"], "DecimalSeparator", ",");% Specify variable properties

% Import the table
CosTauTable = readtable('Hertz table.txt', opts);   %File must be in the same folder as the script
% Convert to numeric matrix
CosTauTable = table2array(CosTauTable);

% Clear temporary variables
clear opts

%%
%Initializing the vectors for interations containing the END-OF-STROKE pysical and geometrical quantities

s= zeros(N,1);          %initializing the incremental sliding distance vector
a= zeros(N+1,1);          %initializing the major contact axis vector
b= zeros(N,1);          %initializing the minor contact axis vector
delta= zeros(N,1);      %initializing the contact body convergence vector
d= zeros(N+1,1);        %initializing the wear track width vector
Aw= zeros(N,1);         %initializing the cross-sectional area of wear track vector
r_t= zeros(N,1);        %initializing the wear track profile radius vector
theta= zeros(N,1);      %initializing the wear track angular amplitude vector
arc= zeros(N,1);        %initializing the wear track arc profile length vector
hw= zeros(N+1,1);         %initializing the underfomed wear depth vector
h_tot= zeros(N,1);     %initializing the wear depth vector + elastic deformation
corrections= zeros(N,2);    %initializing the correction flags matrix
%column 1 for r_t limited to r_s
%column 2 for d limited to linearized arc profile

%initializing the major contact axis vector, as it's done for the Hamrock
%and dawson simplified formulas, it's done for the handbook's formulas
% _t is the subscript to indicate the values came from the Table of
% Schaeffler handbook

a_t= zeros(N,1);
b_t= zeros(N,1);
delta_t= zeros(N,1);
hw_t= zeros(N,1);
hw_tot_t= zeros(N,1);
cos_tau=zeros(N,1);

%A_0= pi*(a_0)^2; Hertzian contact area sphere-plane
%pm_0= F/A_0; Hertzian mean contact pressure

p= zeros(N,1);              %initializing the pressure of contact of the pin
Var= zeros(N,1);            %initializing the Volume of wear by Archard
Vsph= zeros(N,1);           %initializing the Volume of wear of the pin
he= zeros(N+1,1);           %initializing the displacement elastic of the pin

k=k/13*0.1;                 %I don't know the value of k for the pin, as first tentative I set 0.1 k of the disc


%% Start of the computing

% Wear increment loop
for j=1:N                                       %j -> stroke counter
    s(j)= j*l;                                  %cumulative sliding distance
    %     deltaAw= k(j)*F*0.1;                  %At the moment I don't know the exact value of k for the pin so I use 0.1 as first tentative
    
    % It is assumed that the wear volume is inscribed into a circular profile if elastic deformation is recovered
    % Circular segment area formula: A = 0.5 * r_t^2 * (2*theta - sen(2*theta))
    % Wear track radius: r_t= a_0 * 1/sen(theta)
    % Wear track angular amplitude -> theta
    % Wear track widht --> d
    % r_t, theta and d are unknown -> system of 2 equations for 3 unknowns: one must add an hypotesys on wear track width d
    
    if j==1 %First iteration (first stroke)
        
        a_0= ( 3*F/(4/d_s)*Er )^(1/3);                  %Hertzian contact RADIUS sphere-plane at first contact (Shaeffler handbook formula)
        a(j)=a_0;
        d(j) =2*a_0;                                    %I know the track amplitude from the first iteration -> the width of the sphere-plane hertzian contact
        p(j)=F/(3.14*a(j)^2);                           %Pressure of contact of spherical pin
        hw(j+1)=k(j)*p(j)*l;                            %The height of the cutted part of the pin, instead of s there is a delta s --> l
        
        
        
        a(j+1)=(d_s*hw(j+1)-hw(j+1)^2)^(1/2);     %The radius of the cutted part of the pin
        
        Var(j)= k(j)*F*s(j);                            %Volume of wear for Archard
        Vsph(j)=3.14*hw(j)/6*(3*a(j)^2+hw(j)^2);        %Volume of wear computed by the geometry of the Sphere
        
        
    else %strokes else than the first
        
        
        p(j)=F/(3.14*a(j)^2);                           %Pressure of contact of spherical pin
        hw(j+1)=k(j)*p(j)*l + hw(j);                    %The height of the cutted part of the pin, instead of s there is a delta s --> l
        d(j)=2*a(j);                                     %The dimension of contact zone is now equal to 2*a
        
            
        a(j+1)=(d_s*hw(j+1)-hw(j+1)^2)^(1/2);     %The radius of the cutted part of the pin
        
        
        Var(j)= k(j)*F*s(j);                            %Volume of wear for Archard
        Vsph(j)=3.14*hw(j)/6*(3*a(j)^2+hw(j)^2);        %Volume of wear computed by the geometry of the Sphere
        
        
        
    end
    
end

%% Results
close all
% figure()
% plot(s,(hw_tot*10^6),'DisplayName','Total wear penetration depth');
% title('Total wear penetration depth [um]');
% xlabel('Test duration [m]');
%
% figure()
% plot(s,(Aw*10^12),'DisplayName','Simulated cross-sectional area')
% hold on
% plot(WearTrackArea(:,1),WearTrackArea(:,2),'xr','DisplayName','Exprimental values')
% title('Cross-sectional wear track area [um^2]')
% xlabel('Test duration [m]')
% legend('location','southeast')
% hold off
%
% figure()
% loglog(s,(r_t*10^3),'DisplayName','Wear track profile radius')
% title('Wear track profile radius [mm]')
% xlabel('Test duration [m]')

figure()
plot(s,(d(1:end-1)*10^3),'DisplayName','Wear track amplitude')
title('Wear track width [mm]')
xlabel('Test duration [m]')

figure()
plot(s,(hw(1:end-1)*10^6),'DisplayName','Simulated wear track depth')
hold on
plot(WearTrackDepth(:,1),WearTrackDepth(:,2),'xr','DisplayName','Exprimental values')
title('Wear track maximum depth [um]')
xlabel('Test duration [m]')
legend('location','southeast')
hold off

figure()
plot(s,Var,'DisplayName','Volume of wear Archard')
hold on
plot(s,Vsph,'DisplayName','Volume of wear Spherical part')
title('Volume of wear Comparison [um3]')
xlabel('Test duration [m]')
legend('location','southeast')
hold off


% figure()
% plot(x,z,'DisplayName','Wear track geometry')
% title('Wear track maximum depth [um]')
% xlabel('Test duration [m]')
% legend('location','southeast')
% hold off

% figure()
% plot(1:size_1, results(:,1)*10^6, 'DisplayName', 'Hamrock & Dawson')
% hold on
% plot(1:size_1, results(:,4)*10^6,'.','MarkerIndices',1:10:size_1, 'DisplayName', 'Hertz theory')
% title('Major contact axis [um]')
% xlabel('Stroke count')
% legend('location','northwest')
% hold off
%
% figure()
% plot(1:size_1, results(:,2)*10^6, 'DisplayName', 'Hamrock & Dawson')
% hold on
% plot(1:size_1, results(:,5)*10^6, '.','MarkerIndices',1:10:size_1, 'DisplayName', 'Hertz theory')
% title('Minor contact axis [um]')
% xlabel('Stroke count')
% legend('location','southwest')
% hold off

% figure()
% plot(1:size_1, results(:,3), 'DisplayName', 'Hamrock & Dawson')
% hold on
% plot(1:size_1, results(:,6), 'DisplayName', 'Hertz theory')
% title('contact body convergence')
% xlabel('Stroke count')
% legend('location','northeast')
% hold off

% figure()
% plot(1:size_1, hw_tot(1:size_1)*10^6, 'DisplayName', 'Hamrock & Dawson')
% hold on
% plot(1:size_1, hw_tot_t(1:size_1)*10^6,'r--', 'DisplayName', 'Hertz theory')
% title('Wear track maximum depth [um]')
% xlabel('Stroke count')
% legend('location','southwest')
% hold off

% function h = circle2(x,y,r)
% d = r_t(N)*2;
% px = -r;
% py = -h_w(N)-r_t(N);
% h = rectangle('Position',[px py d d],'Curvature',[1,1]);
% daspect([1,1,1])
% end

