clear
%Propriétés des matériaux
R1=630; %Rayon de la roue en mm
R2=inf ; % largeur du rail en mm
E1=210000; %Module d’élasticité du matériau 1 en MPa
E2=210000; %Module d’élasticité du matériau 2 en MPa
v1=0.3; %Coefficient de Poisson du matériau 1
v2=0.3; %Coefficient de Poisson du matériau 2
u=0.5 %Coefficient de frottement
p= 104125; % Charge appliquée en N
L=135; %largeur entre roue/rail en mm
q=p/L ; %charge appliqué entre roue rail en N/mm
%Mode de la sortie
%Input choice=1 pour résultat adimensionné
%Choice=2 pour résultat normale (avec dimension)
choice=2;
R=1/(1/R1+1/R2); %Rayon composé en mm
E=1/((1-v1^2)/E1+(1-v2^2)/E2); %Module d’élasticité composé en MPa
%Rayon de contact
a=(4*q*R/(pi*E))^(1/2); %Rayon de contact en mm
x=[-2*a:.01*3*a:2*a]; %Discrétisation de l’axe x
z=[0:.005*3*a:2*a]; % Discrétisation de l’axe z
p0=(q*E/(pi*R))^(1/2); %Pression maximal pour contact linéique en MPa
% Ci-dessous la boucle qui calcule les matrices de contrainte sous la surface :
for i=1:length(z);
for j=1:length(x);
xx(i, j)=x(j);
zz(i, j)=z(i);
end
end
s=[1/2*[((a^2)-(z.^2)-((x.^2))).^2+4*(a^2)*(x.^2).^1/2-((a^2)-(x.^2)-((z.^2)))]].^(1/2);
sx=(u*p0*(2*x./a)).*(1-(s./((a^2+s.^2).^(1/2))))+((x.*z.^2.*s.*a)./((a^2+s.^2).^(1/2)).*(s.^4+z.^2.*a^2)) ;
sy=(2*v1*x./a).*(1-(s./((a^2+s.^2)))) ;
sz=(u*p0)+((x.*z.^2.*s.*a)./ ((a^2+s.^2).^(1/2)).*(s.^4+a^2*z.^2)) ;
sxz=((u*p0)*(z./a)).*(2-(s./((a^2+s.^2).^(1/2))- ((a^2+s.^2).^(1/2)./s).*(x.^2.*s.^3*a.^2)./(a^2+s.^2).^(3/2).*(s.^4+z.^2*a^2))) ;
%OUTPUT
if choice==1
figure('name',' Contours contraintes sous surface xz, Modèle de Hertz');
subplot(231)
contour(xx/a, -zz/a, sx/p0)
C = contour(xx/a, -zz/a, sx/p0);
clabel(C)
xlabel('x/a')
ylabel('z/a')
title('Contrainte Sx')
subplot(232)
contour(xx/a, -zz/a, sy/p0)
C = contour(xx/a, -zz/a, sy/p0);
clabel(C)
xlabel('x/a')
ylabel('z/a')
title('Contrainte Sy')
subplot(233)
contour(xx/a, -zz/a, sz/p0)
C = contour(xx/a, -zz/a, sz/p0);
clabel(C)
xlabel('x/a')
ylabel('z/a')
title('Contrainte Sz')
subplot(234)
contour(xx/a, -zz/a, sxz/p0)
C = contour(xx/a, -zz/a, sxz/p0);
clabel(C)
xlabel('x/a')
ylabel('z/a')
title('txz')
subplot(235)
contour(xx/a, -zz/a, sxz/p0)
C = contour(xx/a, -zz/a, sxz/p0);
clabel(C)
xlabel('x/a')
ylabel('z/a')
title('tmax')
subplot(236)
plot(x/a,P./p0)
title('Pression normale')
xlabel('x/a')
ylabel('P/p0')
end
if choice==2
figure('name',' Contours contraintes sous surface xz: Modèle de Hertz');
subplot(231)
contour(xx, -zz, sx)
C = contour(xx, -zz, sx);
clabel(C)
xlabel('x en mm')
ylabel('z en mm')
title('Sxx [MPa]')
subplot(232)
contour(xx, -zz, sy)
C = contour(xx, -zz, sy);
clabel(C)
xlabel('x en mm')
ylabel('z en mm')
title('Syy [MPa]')
subplot(233)
contour(xx, -zz, sz)
C = contour(xx, -zz, sz);
clabel(C)
xlabel('x en mm')
ylabel('z en mm')
title('Szz [MPa]')
subplot(234)
contour(xx, -zz, txz)
C = contour(xx, -zz, sxz);
clabel(C)
xlabel('x en mm')
ylabel('z en mm')
title('txz [MPa]')
subplot(235)
contour(xx, -zz, sxz)
C = contour(xx, -zz, sxz);
clabel(C)
xlabel('x en mm')
ylabel('z en mm')
title('tmax [MPa]')
subplot(236)
plot(x,P)
title('P[MPa]')
xlabel('x en mm')
ylabel('P en MPa')
end
Best Answer