A Motorcycle Nonlinear Model for Virtual Rider design.

R.Lot and M.Massaro, Copyright (C) 2010.

> restart:

read MBSymba package for multibody modeling

(Ref. R. Lot, M. Da Lio, A Symbolic Approach for Automatic Generation of the Equations of Motion of Multibody Systems, Multibody System Dynamics n° 12, Kluwer Academic Publishers. pp. 147-172, 2004)

> libname:="C:\libs",libname:with(MBSymba_r4):

Warning, Non-Linear Modeling option has been choosen

use compact form for cos and sin

> alias(s=sin,c=cos);

Motorcycle Model

Degrees of Fredom (dof)

> u(t),v(t), longitudinal and lateral velocity (dof 1, 2)

> z(t), vertical motion of the chassis (dof 3)

> psi(t), phi(t),mu(t), yaw, roll and pitch angles of the chassis (dof 4,5,6)

> zsr(t), rear suspension travel (dof 7)

> delta(t), steering angle (dof 8)

> beta(t), twist deflection (dof 9)

> zsf(t); front suspension travel (10)

> PDEtools[declare](%,prime=t,quiet):

Inputs

> uvars:=[tau(t), steering torque

Mr(t), rear wheel torque

Mf(t) front wheel torque

];PDEtools[declare](%,prime=t,quiet):

the following variables are assumed "small"

> linear_modeling({ delta(t), psi4(t), phi4(t)=phi(t), beta(t) });

Warning, Linear Modeling option has been choosen for the following variables:

Frames and Bodies

Define a moving reference frame with longitudinal velocity u(t), lateral velocity v(t) and angular velocity psi[dot](t). Note that this is not an inertial frame.

> USE_master_frame(u(t),v(t),0, 0,0,psi[dot](t));

> T1 := master_frame:

Warning, The base frame <moving_frame> has been selected. The moving velocity matrix is:

define a proceure for compensating the inertia already included in the body "motorcycle"

> anti_body := proc(BB::BODY)

local anti_BB:

anti_BB := copy(BB):

anti_BB[frame] := map(eval,subs( seq(args[i],i=2..nargs), evalm(BB[frame]))):

anti_BB[mass] := -BB[mass]:

anti_BB[inertia] := evalm(-BB[inertia]):

copy(anti_BB):

end:

body "motorcycle"

> T2 := T1 * rotate('X',phi(t)) * translate(0,0,z(t)-H) * rotate('Y',mu(t)):

> G := origin(T2): #show(G);

> CYZ := 0: CXY:=0:because of the simmetry with respect XY plane

> motorcycle := make_BODY(G,M,IX,IY,IZ,CYZ,CXZ,CXY): show(motorcycle);

body "rear_vheel"

> TWR0 := T2 * translate(-B,0,H-Rr+zsr(t)) * rotate('Y',-mu(t)):

> TWR := TWR0 * rotate('Y',theta1(t)):

> rear_wheel := make_BODY(TWR,umr,0,iry,0): show(rear_wheel,'synthetic');

> anti_rear_wheel := anti_body(rear_wheel, zsr(t)=0, theta1(t)=0): #show(anti_rear_wheel);

body "front_chassis"

> T4 := T2 * rotate('Y',-mu(t))* rotate('Y',mu(t)+epsilon) * translate(a1,0,0) * rotate('X',beta(t)) * rotate('Z',delta(t)):

> T40:= T2 * rotate('Y',-mu(t))* rotate('Y',mu(t)+epsilon) * translate(a1,0,0) * rotate('X',beta(t)):

> Gf := make_POINT(T4,gfx,0,gfz): show(Gf,'syntetic');

> front_chassis := make_BODY(Gf,mf,Ifx,Ify,Ifz): show(front_chassis,'syntetic');

> anti_front_chassis := anti_body(front_chassis,delta(t)=0,beta(t)=0): #show(anti_front_chassis);

body "front_wheel"

> psi4(t):='psi4(t)':front wheel yaw

> phi4(t):='phi4(t)':front wheel roll

> T4p := rotate('Z',psi4(t)) * rotate('X',phi4(t)): #show(T4p);

> zero := project( make_VECTOR(T4,0,1,0)-make_VECTOR(T4p,0,1,0), T2 ):

front wheel yaw and roll as a function of other dof

> front_angles := combine (solve ( {Xcomp(zero),Zcomp(zero)} , {psi4(t),phi4(t)}));

> WF := make_POINT(T4,lx,0,lz+zsf(t)): WF := project(WF,T1):

> TWF0 := translate( Comps(WF) ) * linearize( evalm(T4p )):

> TWF := TWF0 * rotate('Y',theta4(t)):

> front_wheel := make_BODY(TWF,umf,0,ify,0): show(front_wheel,'syntethic');

> anti_front_wheel := anti_body(front_wheel, zsf(t)=0,theta4(t)=0):

Forces

> _gravity := make_VECTOR(master_frame,0,0,+g): show(_gravity);

rider steering torque + steering damper

> steering_torque := make_TORQUE(T4, STX,STY,tau(t)-Cdelta*diff(delta(t),t), front_chassis,motorcycle):

> show(steering_torque,'syntetic');

drag Force

> CA := make_POINT(T2,0,0,H-HCA): show(CA,'synthetic');

> drag_force := make_FORCE(T1,-kD*VX(t)^2,0,0,CA,motorcycle):

> show(drag_force);

suspensions

> v0 :=1/20:

> cr := (crc+cre)/2 +(cre-crc)*arctan(diff(zsr(t),t)/v0)/Pi:

> cf := (cfc+cfe)/2 +(cfe-cfc)*arctan(diff(zsf(t),t)/v0)/Pi:

> - ( RSZ0 + kr*zsr(t) + cr*diff(zsr(t),t) ): molla + ammortizzatore

> RSZ(t) := subs(diff(zsr(t),t)=zsr_dot(t), %):

> - ( FSZ0 + kf*zsf(t) + cf*diff(zsf(t),t) ): molla + ammortizzatore

> FSZ(t) := subs(diff(zsf(t),t)=zsf_dot(t), %):

> rear_suspension := make_FORCE(T2, RSX,RSY,RSZ(t), origin(rear_wheel [frame]), rear_wheel ): show(rear_suspension,synthetic);

> front_suspension := make_FORCE(T4, FSX,FSY,FSZ(t), origin(front_wheel[frame]), front_wheel): show(front_suspension,synthetic);

rear tyre

> Fr(t) := fr(t)*M*g: Sr(t) := -Mr(t)/Rr:

> rear_tire_geometry := [YT1 = rt1*sin(phi(t)), ZT1 = rho1 + rt1*cos(phi(t))]: carcass_geometry=Vector(%);

> rear_tire := make_TIRE(TWR0,theta1(t), YT1,ZT1, 0,0 ,Sr(t),Fr(t),-Nr(t), 0,

uR*Nr(t),rolling resistance

mrR * phi(t) *Nr(t)* (1+twR* phi(t)^2) +Mzsar twisting torque + self-aligning torque

,phi(t),0):

show(rear_tire,synthetic);

front tyre

> Ff(t):=ff(t)*M*g: Sf(t):=-Mf(t)*M*g:

front_tire_geometry := [YT4 = rt4*sin(phi4(t)), ZT4 = rho4 + rt4*cos(phi4(t)) ]:

front_tire_geometry := linearize( subs(front_angles,front_tire_geometry) , small_vars union {mu(t)},t):

'carcass_geometry'=Vector(%);

> front_tire := make_TIRE(TWF0,theta4(t), YT4,ZT4, 0,0 ,Sf(t),Ff(t),-Nf(t), 0,

uF*Nf(t),

mrF * subs(front_angles,phi4(t)) *Nf(t)* (1+twF* subs(front_angles,phi4(t))^2)+Mzsaf

,subs(front_angles,phi4(t)),0):

show(front_tire,synthetic);

twist torque

> flex_torque := make_TORQUE(T40, -kb*beta(t)-cb*diff(beta(t),t) ,0,0,front_chassis,motorcycle):show(flex_torque,'synthetic');

Newton-Euler Equations of Motion

assembly multibody model

> Nr:='Nr': Nf:='Nf':

> bodies := {motorcycle,front_chassis, anti_front_chassis, rear_wheel, anti_rear_wheel, front_wheel, anti_front_wheel};

> forces := {rear_tire, front_tire, drag_force, rear_wheel_torque,front_wheel_torque, flex_torque};

Newton Equations for the Whole Vehicle

> eqnsN := newton_equations(bodies union forces,verbose):

BODIES:

1 - anti_rear_wheel

2 - rear_tire_body

3 - motorcycle

4 - rear_wheel

5 - front_chassis

6 - anti_front_chassis

7 - anti_front_wheel

8 - front_wheel

9 - front_tire_body

EXTERNAL FORCES:

10 - front_tire_force

11 - rear_tire_force

12 - drag_force

WARNING: the following objects will not appear in the equations of motion:

13 - rear_wheel_torque

14 - front_wheel_torque

15 - flex_torque

store equilibrium in the X direction

> eqn_X := Xcomp(eqnsN,master_frame):

store equilibrium in the Y direction

> eqn_Y := Ycomp(eqnsN,master_frame):

store equilibrium in the Z direction

> eqn_Z := Zcomp(eqnsN,master_frame):

Euler Equations for the Whole Vehicle

> eqnsE := euler_equations(bodies union forces, G, verbose): eqnsE := project(eqnsE, motorcycle[frame]):

BODIES:

1 - anti_rear_wheel

2 - motorcycle

3 - rear_wheel

4 - front_chassis

5 - rear_tire_bodyT

6 - front_tire_bodyF

7 - rear_tire_bodyF

8 - anti_front_chassis

9 - anti_front_wheel

10 - front_wheel

11 - front_tire_bodyT

EXTERNAL FORCES:

12 - front_tire_force

13 - rear_tire_force

14 - drag_force

EXTERNAL TORQUES:

15 - rear_tire_torque

16 - front_tire_torque

WARNING: the following objects will not appear in the equations of motion:

17 - rear_wheel_torque

18 - front_wheel_torque

19 - flex_torque

store roll equilibrium

> eqn_phi := Xcomp(eqnsE):

store pitch equilibrium

> eqn_mu := Ycomp(eqnsE):

store yaw equilibrium

> eqn_psi := Zcomp(eqnsE):

kinematic equations for 1st order formulation

collect equations

> eqns_NE := [ eqn_X,eqn_Y,eqn_Z, eqn_phi, eqn_mu, eqn_psi]: #Vector(%);

linear velocities of th chassis

> VG2 := velocity(origin(T2)): #show(VG2);

angular velocities of the chassis

> Omega := angular_velocity(motorcycle[frame]): #show(Omega, synthetic);

"kinematic equations" for the 1st order formulation

> equ_omega := [

VX(t) - Xcomp(VG2),

VY(t) - Ycomp(VG2),

VZ(t) - Zcomp(VG2),

omega_X(t) - Xcomp(Omega),

omega_Y(t) - Ycomp(Omega),

omega_Z(t) - Zcomp(Omega)

> Sspeed := expand(op(solve(equ_omega, [u(t),v(t),diff(z(t),t),diff(phi(t),t), diff(mu(t),t), psi[dot](t)]))):

> Sspeed := Sspeed[1..6]:Vector(%);

> simplify_speed := x -> expand(simplify(subs(Sspeed,expand(subs(Sspeed,x))))):

> eqns_NE0 := expand(subs(front_angles, eqns_NE)):

> eqns_NE0 := simplify(subs(expand(map(simplify_speed, eqns_NE0)))):

Euler Equations about the steering axis and twist axis

> eqnsE_steer := euler_equations( {front_chassis, front_wheel, anti_front_wheel, steering_torque,front_tire,

flex_torque}, origin(T4),verbose): eqnsE_steer := project(eqnsE_steer,T4):

BODIES:

1 - front_chassis

2 - front_tire_bodyF

3 - anti_front_wheel

4 - front_wheel

5 - front_tire_bodyT

EXTERNAL FORCES:

6 - front_tire_force

EXTERNAL TORQUES:

7 - steering_torque

8 - front_tire_torque

9 - flex_torque

equation related to twist dof

> eqn_beta:= Xcomp(eqnsE_steer,T4):

eqn_beta:= expand( subs(front_angles, eqn_beta)):

eqn_beta:= simplify(subs(expand(map(simplify_speed, eqn_beta)))):

equation related to steer dof

> eqn_delta := simplify(expand(Zcomp(eqnsE_steer,T4))):

> eqn_delta0 := expand( subs(front_angles, eqn_delta)):

> eqn_delta0 := simplify(subs(expand(map(simplify_speed, eqn_delta0)))):

Newton Equations for suspensions

rear

> eqnsN_r_wheel := newton_equations({rear_wheel,rear_suspension,rear_tire},verbose):

store equation

> eqn_zsr := Zcomp(eqnsN_r_wheel,T2):

BODIES:

1 - rear_tire_body

2 - rear_wheel

EXTERNAL FORCES:

3 - rear_tire_force

4 - rear_suspension

front

> eqnsN_f_wheel := newton_equations({front_wheel,front_suspension,front_tire},verbose):

store equation

> eqn_zsf := Zcomp(eqnsN_f_wheel,T4):

eqn_zsf := expand( subs(front_angles, eqn_zsf)):

BODIES:

1 - front_wheel

2 - front_tire_body

EXTERNAL FORCES:

3 - front_tire_force

4 - front_suspension

Tyre Equations

rear

kinematics

> (CC1, wC1,VS1,VN1,VR1) := tire_kinematics(rear_tire):

contact point

> CC1 := project(CC1,T1): show(CC1);

contact point velocity

> VS1 := VX(t);

> VN1 := subs( diff(zsr(t),t)=zsr_dot(t),simplify_speed(VN1));

vertical deflection and load

> xir := Zcomp(CC1,master_frame);

> Nr(t) := M*g*(1-B/p) + kpr * xir:

> nr(t) := (1-B/p) + kpr * xir /(M*g);

non-linear force

> mgcfr(t) := Dr*sin( Crl*arctan(Brl*'slipR'-Erl*(Brl*'slipR'-arctan(Brl*'slipR')))

+ Crp*arctan(Brp* phi(t) -Erp*(Brp* phi(t) -arctan(Brp* phi(t) ))) );

relaxation equation

> eqn_rtire :=(sigmaR0+dsigmaR*(nr(t)-nr0)*M*g)/VS1*diff(fr(t),t) + fr(t) - nr(t) * ( mgcfr(t)-iphi*mgcfr[phi](t) );

yaw torque

> Mzsar:= Nr(t)*Drsa*sin(Crsa*atan(Brsa*slipR-Ersa*(Brsa*slipR-atan(Brsa*slipR))));

front

> (CC4, wC4,VS4,VN4,VR4) := tire_kinematics(front_tire):

> CC4 := project(CC4,T1): #show(CC4,synthetic); #show(wC4);

> VS4 := VX(t):

> VN4 := expand(subs(front_angles,VN4)):

> xif := subs(front_angles,Zcomp(CC4,master_frame)):

> Nf(t) := M*g*B/p + kpf * xif :

> nf(t) := B/p + kpf * xif / (M*g):

> mgcff(t) := Df*sin( Cfl*arctan(Bfl*'slipF'-Efl*(Bfl*'slipF'-arctan(Bfl*'slipF')))

+ Cfp*arctan(Bfp*subs(front_angles,phi4(t))-Efp*(Bfp*subs(front_angles,phi4(t))-arctan(Bfp*subs(front_angles,phi4(t)) ))) ):

> eqn_ftire := (sigmaF0+dsigmaF*(nf(t)-nf0)*M*g)/VS4*diff(ff(t),t) + ff(t) - nf(t) * ( mgcff(t) -iphi*mgcff[phi](t)):

> Mzsaf:= Nf(t)*Dfsa*sin(Cfsa*atan(Bfsa*slipF-Efsa*(Bfsa*slipF-atan(Bfsa*slipF)))):

Assembly Equations and Variables

state variables

> tvars0 := [VX(t),z(t),VZ(t),mu(t),omega_Y(t),

VY(t),omega_Z(t), phi(t),omega_X(t), delta(t), delta_dot(t),

beta(t),beta_dot(t),

fr(t), ff(t),

zsr(t), zsr_dot(t), zsf(t), zsf_dot(t)];

> acc_vars := [delta(t), zsr(t), zsf(t),beta(t)]:

collect equations

> teqns0 :=

map(numer, [ lhs(Sspeed[3])-rhs(Sspeed[3]), lhs(Sspeed[4])-rhs(Sspeed[4]), lhs(Sspeed[5])-rhs(Sspeed[5])])kinematic equation

union down_order(eqns_NE0, acc_vars, t) whole vehicle equations

union [diff(delta(t),t)-delta[dot](t) , down_order(eqn_delta0,acc_vars, t)]steer equation

union [diff(beta(t),t)-beta[dot](t) , down_order(eqn_beta,acc_vars, t)]twist equation

union down_order([eqn_rtire, eqn_ftire],acc_vars,t)tyre equations

union [diff(zsr(t),t)-zsr_dot(t)] union down_order([eqn_zsr],acc_vars,t)rear suspension equation

union [diff(zsf(t),t)-zsf_dot(t)] union down_order([eqn_zsf],acc_vars,t):front suspension equation

nops(teqns0);

compute wheel spin acceleration

> tire_geometry := rear_tire_geometry union subs(delta(t)=0,front_tire_geometry): Vector(%);

> wheels_spin := subs(tire_geometry,[

diff(theta1(t),t,t) = -diff(VX(t),t)/ZT1, diff(theta1(t),t) = -VX(t)/ZT1,

diff(theta4(t),t,t) = -diff(VX(t),t)/ZT4, diff(theta4(t),t) = -VX(t)/ZT1

]): Vector(%);

update equations

> teqns0 := subs(tire_geometry,wheels_spin,teqns0):

rear_rolling_radius:=rho1+rt1*cos(phi(t));

front_rolling_radius:=rho4+rt4*cos(phi(t));

display equations

> for i from 1 to nops(teqns0) do i=tvars0[i]; teqns0[i]; end do;

M(t)*dx/dt = F(x,u)

> (MM,FF) := LinearAlgebra[GenerateMatrix](teqns0 , diff(tvars0,t)):

plots[sparsematrixplot](MM);display mass matrix

> Example of Vehicle Parameters and Vibration Modes

g = 9.8070E+00

p = 1.4150E+00

epsilon = 4.1840E-01

an = 9.1240E-02

M = 2.5560E+02

B = 7.0980E-01

H = 6.3980E-01

CXZ =-1.1300E+00

IX = 1.8650E+01

IY = 5.0460E+01

IZ = 3.7140E+01

mf = 3.4630E+01

bf = 1.3150E+00

hf = 4.8170E-01

Ifx = 2.1410E+00

Ify = 2.1980E+00

Ifz = 5.0840E-01

Rr = 3.1670E-01

rho1 = 2.2270E-01

rt1 = 9.4000E-02

KlambdaR = 1.0020E+01

KphiR = 1.2170E+00

umr = 2.5260E+01

mwr = 5.3341E+00

Rf = 2.9100E-01

rho4 = 2.3100E-01

rt4 = 6.0000E-02

umf = 2.5880E+01

mwf = 5.3849E+00

KlambdaF = 1.2880E+01

KphiF = 1.0390E+00

HCA = 8.3270E-01

kD = 3.0000E-01

kpr = 1.3500E+05

kpf = 1.4000E+05

RSZ0 =-1.0630E+03

kr = 2.1320E+04

crc = 9.2840E+02

cre = 2.4836E+03

FSZ0 =-9.1720E+02

kf = 1.8500E+04

cfc = 5.2530E+02

cfe = 9.7470E+02

sigmaR0 = 1.0430E-01

dsigmaR = 7.0000E-05

nr0 = 5.7367E-01

sigmaF0 = 1.0800E-01

dsigmaF = 1.0000E-04

nf0 = 4.2646E-01

Cdelta = 5.0000E+00

Brsa = 7.1500E-01

Crsa = 1.9960E+01

Drsa =-2.0000E-02

Ersa = 1.1000E-01

Bfsa = 4.6650E+00

Cfsa = 4.6510E+00

Dfsa =-1.7000E-02

Efsa = 5.4980E+00

Brl = 8.1628E-01

Crl = 9.4280E+00

Dr = 1.3020E+00

Erl =-1.9000E-01

Brp = 7.4302E-01

Crp = 1.2580E+00

Erp = 3.3400E-01

Bfl = 1.4340E+00

Cfl = 7.2550E+00

Df = 1.2380E+00

Efl = 7.5500E-01

Bfp = 9.1983E-02

Cfp = 9.1240E+00

Efp =-1.5400E+00

kb = 4.5100E+04

cb = 2.5000E+01

uR = 1.1000E-02

uF = 1.1000E-02

mrF = 2.1000E-02

twF = 5.5300E-01

mrR = 2.1000E-02

twR = 2.0430E+00

a1 = 4.7566E-01

lz = 6.0523E-01

gfz = 3.9035E-01

gfx = 1.3100E-02

lx = 2.6993E-02

Root-locus for speeds from 20 to 40 m/s.

-black crosses: lateral acceleration = longitudinal acceleration = 0

-blue circles: lateral acceleration = 5 m/s2

-red triangles: longitudinal acceleration = 5 m/s2

-green squares: longitudinal deceleration = 5 m/s2 (only front wheel braking torque used)

All of the most important vehicle vibration modes are visible.

For additional details on vibration modes see e.g.

V. Cossalter, Motorcycle Dynamics, Lulu.com, 2006, ISBN: 978-1-4303-0861-4, http://www.lulu.com/content/paperback-book/motorcycle-dynamics/221690