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
% N=1000;

%% 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

%%
%Initializing the vectors for interations containing the END-OF-STROKE 
%pysical and geometrical quantities
%I have to take into account that the GIWM of H. is on a disc that is
%moving. Which geometry should i take for the sliding distance?

s= zeros(N,1);              %initializing the incremental sliding distance vector
a= zeros(N,1);              %initializing the major contact axis vector
b= zeros(N,1);              %initializing the minor contact axis vector
hw= zeros(N,1);             %initializing the underfomed wear depth vector
he= zeros(N,1);             %initializing the elastic deformation depth vector
h_tot= zeros(N,1);         %initializing the wear depth vector + elastic deformation
p= zeros(N,1);             %initializing the underfomed wear depth vector
                            
%% The first stroke

Einv=(1-ni_s^2)/E_s+(1-ni_d^2)/E_d;
Ec=1/Einv;
Rp=1/r_s;

a_0=((3*F*Rp)/(4*Ec))^(1/3);
he_0=F/(2*Ec*a_0);

for j=1:N       %j -> stroke counter    
    
    if j==1 %First iteration (first stroke)
        hw_0=0;
        hw(j)=hw_0;
        s_0=0;
        s(j)=s_0;
        a_0=((3*F*Rp)/(4*Ec))^(1/3);
        a(j)=a_0;
        b(j)=a_0;
        he_0=F/(2*Ec*a_0);
        he(j)=he_0;
        p_0=F/(pi*a(j)^2);
        p(j)=p_0;
        
    else
        s(j)=s(j-1)+l;
        hw(j)=hw(j-1)+2*k(j-1)*p(j-1)*b(j-1);       %in the paper k=kd has some issues on the hardness of the material
        h_tot(j)=he(j-1)+hw(j);
        a(j)=(2*Rp*h_tot(j)-h_tot(j)^2)^(1/2);
        p(j)=F/(pi*a(j)*b(j-1));
        he(j)=F/(2*Ec*(a(j)*b(j-1))^(1/2));
        b(j)=2*((F*Rp)/(pi*a(j)*pi*Ec))^(1/2);
        
    end

end


%% Results


figure()
plot(s,(h_tot*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','best')
hold off



