GW Matlab

%model Greenwood and Williamson (GW)
%function [p_g]=greenwood
%====================================
%===Greenwood and Williamson Model===
%====================================
%Input data: vector data_z with the roughness points and vector
%data_x with the respective coordinates
%Output data: percentage of contact real area, deformation, plot with the
%original profile and deformed profile
x=data_x;
rug=data_z;
%load aplied
load=1.9894; %N/mm
%properties of the material
H=2785; %MPa
E1=205000; %MPa
E2=62750; %MPa (%Glass)
niu1=0.29;
niu2=0.2;
Ecom=1/((1-niu1^2)/E1+(1-niu2^2)/E2); %[MPa]
%rug is the vector with one profile of the roughness topography
%the profile will be approached by polynomial functions using the
%Aramki formulation
%determination of ACF length and the coefficient of ACF, ACF length is the
%length where autocorrelation function is 0.368 (=1/e)
[ACF,Lags,Bounds] = autocorr(rug,length(x)-1);
index_ACF_0368=1;
while ACF(index_ACF_0368)>0.368
index_ACF_0368=index_ACF_0368+1;
end
%plot(x,ACF); %plot whith the fuction of autorrelation
length_ACF=x(index_ACF_0368)-x(1);
alfa=1/length_ACF;
%standard deviation:
sigma=std(rug);
%definition of a vector L_peak (peaks), obtained considering the cross with the reference line:
n=length(x);
k=1;
for i=1:n-1
if ((rug(i)<=0)&(rug(i+1)>0));
j=i+1;
while ((rug(j)>=0)&(j+1<n))
if rug(j+1)<0
Lpeak(k,1)=x(i)-rug(i)*(x(i+1)-x(i))/(rug(i+1)-rug(i));
Lpeak(k,2)=x(j)-rug(j)*(x(j+1)-x(j))/(rug(j+1)-rug(j));
L_peak(k)=Lpeak(k,2)-Lpeak(k,1);
k=k+1;
end
j=j+1;
end
end
end
%definition of a vector L_valley (valleys), obtained considering the cross with the reference line:
k=1;
for i=1:n-1
if ((rug(i)>=0)&(rug(i+1)<0));
j=i+1;
while ((rug(j)<=0)&(j+1<n))
if rug(j+1)>0
Lvalley(k,1)=x(i)-rug(i)*(x(i+1)-x(i))/(rug(i+1)-rug(i));
Lvalley(k,2)=x(j)-rug(j)*(x(j+1)-x(j))/(rug(j+1)-rug(j));
L_valley(k)=Lvalley(k,2)-Lvalley(k,1);
k=k+1;
end
j=j+1;
end
end
end
%create one vector whith the x positions of the all crossings whith the
%reference line
for i=1:(length(L_peak))
Lp(i)=Lpeak(i,1);
Lp(i+length(L_peak))=Lpeak(i,2);
end
for i=1:(length(L_valley))
Lv(i)=Lvalley(i,1);
Lv(i+length(L_valley))=Lvalley(i,2);
end
%vector X that contain all x positions (positions of rough points and the
%crossings)
X=[];
X=[x';Lp';Lv'];
X=unique(X);
X=sort(X);
%create one new vector RUG with the same length that X
RUG=[];
for i=1:length(X)
for j=1:length(x)
if X(i)==x(j)
RUG(i)=rug(j); %the others positions RUG=0
end
end
end
%generation of the one profile approach by parabolas
csi_peak=L_peak*sqrt(2/pi)*alfa*sigma; %equation 8 Aramki part I
csi_valley=L_valley*sqrt(2/pi)*alfa*sigma;
mean_L_peak=(mean(L_peak));
mean_L_valley=(mean(L_valley));
mean_L=1/2*(mean(L_peak)+mean(L_valley));
K1_peak=8*(csi_peak)./(L_peak.^2); %equation 9-b Aramki part I
K1_valley=8*(csi_valley)./(L_valley.^2);
%generation of the vector with points that represent parabolas
%start the vectors with zeros and the same length that X
parabola=zeros(1,length(X));
parabola_peak=zeros(1,length(X));
parabola_valley=zeros(1,length(X));
for i=1:length(L_peak)
j=find(X==(Lpeak(i,1)));
while (X(j)>=Lpeak(i,1)&X(j)<=Lpeak(i,2))
parabola(j)=-(4*csi_peak(i)/(L_peak(i)^2))*(X(j)-Lpeak(i,1)-L_peak(i)/2)^2+csi_peak(i);
parabola_peak(j)=parabola(j);
j=j+1;
end
end
for i=1:length(L_valley)
j=find(X==(Lvalley(i,1)));
while (X(j)>=Lvalley(i,1)&X(j)<=Lvalley(i,2))
parabola(j)=(4*csi_valley(i)/(L_valley(i)^2))*(X(j)-Lvalley(i,1)-L_valley(i)/2)^2-csi_valley(i);
parabola_valley(j)=parabola(j);
j=j+1;
end
end
temp_rq=0;
for i=1:length(rug)
temp_rq=temp_rq+(rug(i))^2;
end
increment=0.0001; %increment of the displacement delta [micron]

%critic interference for each peak
delta_c(i)=(pi*k*H/(2*Ecom))^2*(1/K1_peak(i));%micron
end
y=max(parabola);
lt=zeros(length(L_peak),1); %vector that indicate if the deformation is elastic (0) or plastic (1)
ltemp=0;
f=zeros(length(L_peak),1); %load applied in each asperity
n_steps=0;
while sum(f)<load
y=y-increment;
n_steps=n_steps+1;
%elastic
if lt(i,1)==0 & (csi_peak(i)-y)>0
f(i)=(4/3)*Ecom*((1e-3/K1_peak(i))^(0.5))*((csi_peak(i)-y)*1e-3)^(3/2);
if (csi_peak(i)-y)>delta_c(i)
lt(i,1)=1;
end
end
%plastic
if lt(i,1)==1 & (csi_peak(i)-y)>0
f(i)=2*pi*(1e-3/K1_peak(i))*(csi_peak(i)-y)*1e-3*H;
end
end
end
%contact area
A_cont=0;
for i=2:length(RUG)
if parabola(i)>=y
A_cont=A_cont+(X(i)-X(i-1));
end
end
%percentage of contact area
A_cont_a=A_cont/X(length(X))
%plot deformed profile
for i=1:length(RUG)
if parabola(i)>y
parabola_de(i)=y;
else
parabola_de(i)=RUG(i);
end
end
plot(X,parabola,'k');
hold on;
plot(X,parabola_de,'LineWidth',1.4);
hold off;
%deformation
deformation=max(parabola)-abs(y)
Undefined function or variable 'data_x'.
Published with MATLAB® 7.0

GW Matlab

Uploaded by

Copyright:

Available Formats

GW Matlab

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

GW Matlab

Uploaded by

Copyright:

Available Formats

%model Greenwood and Williamson (GW)

increment=0.0001; %increment of the displacement delta [micron]

Undefined function or variable 'data_x'.

Published with MATLAB® 7.0

You might also like