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
define system inputs
> inputs:= [U = 25, delta(t)=0.05*sin(Pi*t)];
#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):
> plot_solz(solz);