by Roberto Lot and Matteo Massaro

Vehicle modelling

> restart: libname := libname, "C:/MBSymba": with(MBSymba_r6):LA:=LinearAlgebra:

> alias(S=sin,C=cos);
PDEtools[declare]({x(t),y(t),u(t),v(t),psi(t),omega(t),delta(t)},prime=t,quiet);

Definition of the Moving Frame

> TV := translate(x(t),y(t),0) . rotate('Z',psi(t)):

> track_XY := use_moving_frame(u(t),v(t),0,0,0,omega(t),TV):

Warning, The base frame <moving_frame> has been selected.

Trajectory equations

> Vector(track_XY);

1) vehicle

define main frame

> TV := moving_frame:

define vehicle body

> vehicle := make_BODY(TV,M,IX,IY,IZ): show(vehicle,'synthetic');

position, velocity and acceleration of the vehicle CoM

> G := CoM(vehicle): show(G);

> VG := velocity(G): show(VG);

> AG := velocity(VG): show(AG);

2) Active Forces

rear tire

> R := make_POINT(TV,-b,0,0): show(R);

force

> make_VECTOR(TV, Sr, Fr ,0):
rear_tire := make_FORCE(%, R, vehicle): show(rear_tire);

front tire

> TF := TV * Translate(a,0,0) * Rotate('Z',delta(t));

contact point

> F := origin(TF): show(F,'synthetic');

force

> make_VECTOR(TF, 0, Ff ,0):
front_tire := make_FORCE(%, F, vehicle): show(front_tire);

3) Newton-Euler equations

Newton

> eqnsN := newton_equations({vehicle,rear_tire,front_tire}): show(eqnsN);

Euler

> eqnsE := euler_equations({vehicle,rear_tire,front_tire},CoM(vehicle)):
show(eqnsE);

4) ODE formulation

Dynamics

> dyna_vars := [u(t),v(t),omega(t)];

implicit equations

> dyna_eqns := [comp_X(eqnsN),comp_Y(eqnsN),comp_Z(eqnsE)]: Vector(%);

> #dyna_eqns; diff(dyna_vars,t); op(solve( %% , % ));

explicit equations: ODE formulation

> dyna_eqns := op(solve(dyna_eqns, diff(dyna_vars,t))): Vector(%);

Kinematics

> kin_vars := [x(t),y(t), psi(t)]:
kin_eqns := solve(track_XY,diff(kin_vars,t)):

> Vector(%);

ODE formulation

> #<Vector(dyna_eqns) , Vector(kin_eqns)>;

5) Tire model: lateral force vs sideslip

simplified Pacejka formulation:

> 'sin(arctan(lambda))'; %;

> FT := Dy*sin(arctan(KL/Dy*lambda)) * N;

> FT := unapply(FT,(lambda,N,KL,Dy));

> FT(lambda, 1,KL,1);

> plot( subs(Dy=1.2,KL=12,[FT(lambda,1,KL,Dy),KL*lambda,Dy*signum(lambda)]),lambda=-0.3..0.3,
view=[-0.3..0.3,-1.3..1.3],color=[red,black,blue],thickness=[2,1,1],labels=["lambda[rad]","F/N"]
,title="Dy=1.2 - KL=12");

> stiffness=subs(lambda=0,diff(FT(lambda,1,KL,Dy),lambda));

rear tire kinematics

> VR:=velocity(R):show(VR);

> lambda[R] := -arctan( comp_Y(VR,TV)/comp_X(VR,TV) );

> #Fr := FT('lambda[R]',Nr, KLr,Dyr);

front tire kinematics

> VF := velocity(F):show(VF);

> project(VF,TF): 'VF' = show(%);

> lambda[F] := -arctan( comp_Y(VF,TF)/comp_X(VF,TF) );
#Ff := FT('lambda[F]',Nf, KLf,Dyf);

tire forces

> tire_forces := [Fr = FT(lambda[R],Nr, KLr,Dyr), Ff = FT(lambda[F],Nf, KLf,Dyf)]:

#tire_forces := [Fr = KLr*lambda[R], Ff = KLf*lambda[F]]:
Vector(%);

Save the model

> save (dyna_vars,dyna_eqns,tire_forces,lambda, kin_vars,kin_eqns, "Single-Track-UV");

Numerical Simulations

> restart: libname := libname, "C:/MBSymba": with(MBSymba_r6):LA:=LinearAlgebra:
read("Single-Track-UV"):

define vehicle characteristics

> vehicle_data :=[
a=1.3, from Cr to G
b=1.3, from G to Cf
M=1000,mass
Nr=5000,rear load
Nf=5000,front
IZ=1200,yaw inertia
KLr=24,rear tire cornering stiffness / D
KLf=20, front
Dyr=1.2,rear tire Pacejeka D coefficient
Dyf=1.2 front
];

> #Vector(dyna_eqns);

constant speed

> u(t):=U;

out of plane dynamics

> ODEeqns := dyna_eqns[2..3]: Vector(%);

> vars := [v(t),omega(t)]:

numeric eqns

> assign(tire_forces): #Vector(ODEeqns);
ODEeqnsN:= subs(vehicle_data,ODEeqns):

define simulation Initial conditions

> IC := [v(0)=0, omega(0)=0 ]: 'IC'=Vector(IC);

procedure for plotting solution

> plot_solz := proc(solz)
local p_v,p_delta,p_psidot, p_latacc,p_force,p_slip,th,npts:
th:=1: npts:=800:

yaw angle and steer angle

> p_v := plots[odeplot](solz,[ [ t, v(t), color=blue,thickness=th ]] , numpoints=npts, labels=["t","lateral speed [m/s]"]):
p_delta := plot(subs(inputs,delta(t)),t=0..tf ,color=black,thickness=th, numpoints=npts, labels=["t","steering angle [rad]"] ):
p_psidot := plots[odeplot](solz,[ [ t, omega(t),color=blue,thickness=th] ],numpoints=npts):
p_latacc := plots[odeplot](solz,[ [ t, subs(inputs,U*omega(t)),color=blue,thickness=th] ],numpoints=npts, labels=["t","lateral acceration (approximate)"]):

sideslip angles and lateral forces
under-steering behaviour: front slip > rear

> p_slip:= plots[odeplot](solz,[
[t, subs(vehicle_data,inputs,lambda[R]),color=red , thickness=th],
[t, subs(vehicle_data,inputs,lambda[F]),color=blue ,thickness=th]
],title="sidelslip [rad]"):# ,legend=["rear","front"]):

> p_force:= plots[odeplot](solz,[
[t, subs(vehicle_data,inputs,Fr), color=red , thickness=th],
[t, subs(vehicle_data,inputs,Ff), color=blue ,thickness=th]
],title="lateral force [N]"):# ,legend=["rear","front"]):

> plots[display]( array(1..3,1..2,[ [p_v,p_delta],[p_psidot,p_latacc],[p_slip,p_force]]));
end:

Lane Change

define system inputs

> delta0 := 0.05*5:
inputs:= [
U = 25,
#delta(t)=0.10*Heaviside(t-1)
delta(t)= piecewise(t<1,0,t<2,delta0*(t-1),delta0)*0.2
];
#plot( subs(inputs,delta(t)) , t=0..5,labels=["t[s]","delta[rad]"]);

integratedifferential equations of motion

> tf:=10.0;simulation final time

> solz:=dsolve( subs(inputs,ODEeqnsN union IC), vars, type=numeric, range=0..tf):

trajectory

> plot_solz(solz);

Slalom