close all ; %hertzian contact calculator for a cylinder
clc;
format long
F= 5000 ;
v1=0.3;
v2=0.3;
E1=700000;
E2=700000;
d1=10 ;
d2=100 ;
L= 25;
b=(1/d1)+(1/d2);
m1=(1-v1^2)/E1;
m2=(1-v2^2)/E2;
a=0.173;
Pmax=(2*F)/(pi*a*L);
z= 0.786.*a;
ZA=(1+(z^2/a^2))^(1/2);
qx=(-2*v1*Pmax)*(ZA-abs(z/a));
qy=-Pmax*(((1+2*(z^2/a^2))/ZA)-2*abs(z/a));
qz=-Pmax/ZA;
tx=(qx-qz)/2;
ty=(qy-qz)/2;
%Plot creation
hu=[0:0.05:3]; %the used values
ZI=(1+hu.^2).^(1/2);
qx2=(-2*v1*Pmax)*(ZI-hu);
qy2=-Pmax*(((1+2.*hu.^2)./ZI)-2*hu);
qz2=-Pmax./ZI;
tmax1=(qx2-qz2)./2; %equation to be used when "hu" variable is lower than 0.436
tmax2=(qy2-qz2)./2; %equation to be used while "hu" variable is above 0.436
tmax = zeros(size(hu));
tmax(hu < 0.436) = tmax1(hu < 0.436)/Pmax;
tmax(hu >= 0.436) = tmax2(hu >= 0.436)/Pmax;
plot(hu,tmax)
legend('ty')
title('Stress distribution along cylindrical bodies')
