function [Qnor Qtag] = FContacto(q,v)
global MBody
global Parametr
Param =ParametrosMotocicleta;
a = Param.a_n;
b = Param.b_n;
beta=Parametr(31);
Rt=Parametr(3);
Rd=Parametr(4);
r2x=Parametr(32);
r2z=Parametr(33);
r4x=Parametr(34);
r4z=Parametr(35);
r5x=Parametr(36);
r5z=Parametr(37);
r2=[r2x 0 r2z];
r4=[r4x 0 r4z];
r5=[r5x 0 r5z];
Qnor = zeros(MBody.ncoord,1); %en principio 0 y luego se añaden las =!0
Qtag = zeros(MBody.ncoord,1);
for i = 1:MBody.nce %hay 2 nce
% DETERMINA LA LOCALIZACION DE LOS PUNTOS DE CONTACTO Y LA IND%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if i==1
indent1=q(3) - 0.99 + 0.15*tan( q(6) );
indent=indent1;
MBody.indent(1,MBody.i)=indent;
else
indent2=q(3) - 0.99 - ( 0.65+ 0.45*sin(beta) )*tan( q(6) );
indent=indent2;
MBody.indent(2,MBody.i)=indent;
end
%CALCULO DE LA FUERZA NORMAL DE CONTACTO%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if indent<0 %%----->
Kn = MBody.ContactoElastico(i).kn;
fn = - Kn*( indent );
if i==1
Aw = A2(q);
sj=[0 -pi/2+q(6)+q(9)];
Rprima=-Aw(3,2)/( (Aw(3,1)*cos(sj(2))) + (Aw(3,3)*sin(sj(2))) ) ;
sj(1)=(a*Rprima)/( sqrt( ( Rprima*Rprima) +(b/a)^2) ); %parametro lineal contacto rueda
sj(1)=abs(sj(1));
if q(5)>0
sj(1)=-sj(1);
end
r = Rt-b+b*sqrt(1-(sj(1)/a)^2);
Hp=H2p(q,r*cos(sj(2)), sj(1), r*sin(sj(2)));
H=H2(q);
else
Aw = A5(q);
sj=[0 -pi/2+q(6)-beta+q(8)];
Rprima=-Aw(3,2)/( (Aw(3,1)*cos(sj(2))) + (Aw(3,3)*sin(sj(2))) ) ;
sj(1)=(a*Rprima)/( sqrt( ( Rprima*Rprima) +(b/a)^2) ); %parametro lineal contacto rueda
sj(1)=abs(sj(1));
if q(5)>0
sj(1)=-sj(1);
end
r = Rd-b+b*sqrt(1-(sj(1)/a)^2);
Hp=H5p(q, r*cos(sj(2)), sj(1), r*sin(sj(2)));
H=H5(q);
end
nk=[0 0 1]';
Fcon= nk*fn;
c=30000;
vptocontacto=Hp*v;
vindent=vptocontacto'*nk;
Fdisip=c*vindent*indent*nk;
Fnor= Fcon + Fdisip;
%Wcon = Fcon'*Vp = Fcon' *Hp*dp-->Qcon'=Fcon'*Hp
Qnor = Qnor + Hp' * Fnor;
%CALCULO DE LAs FUERZAs TANGS%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
kc = MBody.ContactoElastico(i).kc;
mux = MBody.ContactoElastico(i).mux;
alfac =MBody.ContactoElastico(i).alfac;
muy = MBody.ContactoElastico(i).muy;
u=Aw(:,2); %eje rotacion rueda es 2nd col
bb=cross(u,nk)/( norm(cross(u,nk)) );
vc=H*v; %veloc centro ruedas
vr= Hp*v ; %veloc pto cont
vslip=vr;
Lslip=-bb'*vslip/( (bb'*vc) );
%slip angle
vector=vc-(nk'*vc)*nk;
vector=cross( bb,vector/norm(vector));
slipANG=-asin( nk'*vector );
if Lslip<=kc %mientras sea menor que el critico
Fx= mux*fn*Lslip/kc;
else
Fx= mux*fn;
end
if slipANG<=alfac%mientras sea menor que el critico
Fy=muy*fn*slipANG/alfac;
else
Fy=muy*fn;
end
%elipse saturacion%
%%%%%%%%%%%%%%%%%%%
MBody.ang(MBody.i)=slipANG;
if (Fx/mux)^2+(Fy/muy)^2 >= fn^2 %si la supera se reemplazan por los sat!!
froz=(Fx/mux)^2+(Fy/muy)^2;
froz=sqrt(froz);
Fxsat=fn*Fx/froz;
Fysat=fn*Fy/froz;
Fx=Fxsat;
Fy=Fysat;
end
Ftag = Fx*bb + Fy*u;
Qtag= Qtag + Hp' * Ftag;
else
fn=0;
end %%%%%<-----
end
end
e-REdING. Biblioteca de la Escuela Superior de Ingenieros de Sevilla.
