(Array indices must be positive integers or logical values.) any help plz

1 Ansicht (letzte 30 Tage)
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
%Réparation de la pression.
P(i)=p0*[1-(x./a).^2].^(1/2);
% Ci-dessous la boucle qui calcule les matrices de contrainte sous la surface :
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

Akzeptierte Antwort

the cyclist
the cyclist am 25 Jun. 2020
The specific reason you get that error is that sx is complex number.
sx is complex because s is complex.
s is complex because when you take the square root, elements 66 and 67 are negative, so you are taking the square root of a negative number.
I did not try to figure out any more than that.
If you don't know some of the tools for debugging MATLAB, you might want to read this documentation.

Weitere Antworten (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by