> restart: libname := libname,"C:/MBSymba": with(MBSymba_r6):
> alias(c=cos,s=sin); 'sin' and 'cosine' functions
Semplificative assumptions
Reference
V. Cossalter, R. Lot, M. Massaro, The influence of Frame Compliance and Rider Mobility on the Scooter Stability. Vehicle Syst. Dyn. 45 (2007), pp. 315-326, link: dx.doi.org/10.1080/00423110600976100
small vars of the model
> linear_modeling({
y(t),lateral displacement
psi(t),yaw angle
phi(t),roll angle
delta(t),steer angle
yp(t),phi[P](t),lateral displacement (whole rider) and lean angle of the rider (upper body)
zeta[F](t),zeta[R](t),front and rear carcass deflections (lateral)
alpha[F](t),beta[F](t),bending and torsion angles of the fork
alpha[R](t),beta[R](t),bending and torsion angles of the swingarm
vyR(t),vyF(t)shaker actions (INPUT)
});
Warning, Linear Modeling option has been choosen for the following variables:
> PDEtools[declare](small_vars,prime=t,quiet):
> x(t) := V*t;straight-running motion
> theta1(t) := -(V*t)/Rr;no slip constraint for the rear wheel
theta4(t) := -(V*t)/Rf;no slip constraint for the front wheel
Bodies
Whole motorcycle chassis
frame TRC: fixed to the rear motorcycle chassis and with the origin in the motorcycle center of mass
> TRC := translate(x(t),y(t),0) * rotate('Z',psi(t)) * rotate('X',phi(t)) * translate(B,0,-H);
> GM := origin(TRC): show(GM,'syntetic');whole motorcycle center of mass
> M,IX,IY,IZ,0,CXZ,inertia properties of the whole motorcycle (rear and front chassis, wheels)
CYZ=0, CXY=0: because of the simmetry with respect to XZ plane
> motorcycle := make_BODY(GM,M,IX,IY,IZ,0,CXZ,0): show(motorcycle);
Front chassis
> TRC * translate(-B,0,-Rr+H) * translate(lf*cos(mu),0,-lf*sin(mu)):swingar pinion
> T40 := % * rotate('Y',epsilon) * translate(l23,0,0):steering axis
> T4 := T40 * rotate('Z',delta(t));frame T4, fixed to the front chassis
> Gf := make_POINT(T4,e,0,f): show(Gf,'syntetic');front chassis center of mass
> mf,Ifx,Ify,Ifz,Cfxz,inertia properties of the front chassis (no wheel)
Cfxy=0, Cfyz:because of the simmetry with respect to the XZ plane
front_chassis := make_BODY(Gf,mf,Ifx,Ify,Ifz,0,Cfxz,0): show(front_chassis,'syntetic');
in order to avoid of computing twice the front chassis inertia (already included in whole_motorcycle), a dummy body with negative inertia is created
> Gf0 := make_POINT(T40,e,0,f):
> _front_chassis := make_BODY(Gf0,-mf,-Ifx,-Ify,-Ifz,0,-Cfxz,0): show(_front_chassis,'syntetic');
Front flexible fork
> T60 := T4 * translate(lx6,0,lz6):lumped stiffness position (front)
T6 := T60 * rotate('X',alpha[F](t)) * rotate('Z',beta[F](t)):T6: fixed to the flexible fork and with the origin in the front lumped stiffness
G6 := make_POINT(T6,gx6,0,gz6): show(G6,'synthetic');flexible fork center of mass
> m6,I6x,I6y,I6z,0,C6xz,inertia properties of the flexible fork (no wheel)
C6xy=0, C6yz=0:because of the simmetry with respect to the XZ plane
flex_fork := make_BODY(G6,m6,I6x,I6y,I6z,0,C6xz,0): show(flex_fork,'synthetic');
in order to avoid of computing twice the front flexible fork inertia (already included in front_chassis), a dummy body with negative inertia is created
> G60 := make_POINT(T60,gx6,0,gz6):
> _flex_fork := make_BODY(G60,-m6,-I6x,-I6y,-I6z,0,-C6xz,0):show(_flex_fork,'synthetic');
>
Flexible swingarm
> T70 := TRC * translate(-B,0,-Rr+H) * rotate('Y',mu) * translate(lx7, 0, lz7):lumped stiffness position (rear)
T7 := T70 * rotate('X',beta[R](t)) * rotate('Z',alpha[R](t)):T7: fixed to flexible swingarm and with the origin in the lumped stiffness
G7 := make_POINT(T7,gx7,0,gz7): show(G7,'synthetic');flexible fork center of mass
> G7,m7,I7x,I7y,I7z,0,C7xz,inertia properties of the flexible fork (no wheel)
C7xy=0, C7yz=0:because of the simmetry with respect to the XZ plane
flex_swingarm := make_BODY(G7,m7,I7x,I7y,I7z,0,C7xz,0):show(flex_swingarm,'synthetic');
in order to avoid of computing twice the flexible swingarm inertia (already included in whole_motorcycle), a dummy body with negative inertia is created
> G70 := make_POINT(T70,gx7,0,gz7):
_flex_swingarm := make_BODY(G70,-m7,-I7x,-I7y,-I7z,0,-C7xz,0):show(_flex_swingarm,'synthetic');
REAR WHEEL
> TR := T7 * translate(-lx7, 0, -lz7) * rotate('Y',-mu):rear wheel reference frame with the origin in the wheel center
> TWR := TR * rotate('Y',theta1(t)):rear wheel rolling frame
mrw,Ird,Ira:inertia properties of the rear wheel
rear_wheel := make_BODY(TWR ,mrw,Ird,Ira,Ird): show(rear_wheel,'syntetic');
REAR DUMMY WHEEL
in order to avoid of computing twice the rear wheel inertia (already included in whole_motorcycle), a dummy body with negative inertia is created
> Ira := mry*Rr^2; spin inertia
TWR0 := TRC * translate(-B,0,-Rr+H):TWR0: fixed to a wheel which don't rotate and with the assumption of rigid swingarm
> _rear_wheel := make_BODY(TWR0,-mrw,-Ird,-Ira,-Ird): show(_rear_wheel,'syntetic');
FRONT WHEEL
> TF := T6 * translate(-lx6+lx , 0 , -lz6+lz) * rotate('Y',-epsilon):front wheel reference frame with the origin in the wheel center
TWF:= TF * rotate('Y',theta4(t)):front wheel rolling frame
front_wheel := make_BODY(TWF,mfw,Ifd,Ifa,Ifd): show(front_wheel,'syntethic');
FRONT DUMMY WHEEL
in order to avoid of computing twice the front wheel inertia (already included in whole_motorcycle), a dummy body with negative inertia is created
> Ifa := mfy*Rf^2; spin inertia
mu:swingarm angle with respect to the horizontal axis X
> TWF0 := T40 * translate(lx,0,lz):
_front_wheel := make_BODY(TWF0,-mfw,-Ifd,-Ifa,-Ifd): show(_front_wheel,'syntethic');
Rider
> bs,hs:lumped stiffness position
> TS := TRC * translate(-B+bs,yp(t),H-hs): introduce lateral mobility
> GP := make_POINT(TS,bp-bs,0,-hp+hs): show(GP,'syntetic');whole rider center of mass
WHOLE BODY
> mp,Ipx,Ipy,Ipz,Cpxz,inertia properties of the whole rider
Cpyz=0,Cpxy=0:because of the simmetry with respect to XZ plane
rider := make_BODY(GP,mp,Ipx,Ipy,Ipz,0,Cpxz,0): show(rider,'syntetic');
UPPER BODY
> Tup:= TS * rotate('X',phi[P](t)):Tup: fixed to the upper rider body, which rotate with respect to his hip
GPup:=make_POINT(Tup,bpUP-bs,0,-hpUP+hs):
mpUP,IpxUP,IpyUP,IpzUP,0,CpxzUP:inertia properties of the upper rider body
upper_rider:=make_BODY(GPup,mpUP,IpxUP,IpyUP,IpzUP,0,CpxzUP,0): show(upper_rider,'syntetic');
DUMMY UPPER RIDER
in order to avoid of computing twice the upper rider inertia (already included in rider body), a dummy body with negative inertia is created
> GPup0:=make_POINT(TS,bpUP-bs,0,-hpUP+hs): #show(GPup0);
_upper_rider := make_BODY(GPup0,-mpUP,-IpxUP,-IpyUP,-IpzUP,0,-CpxzUP,0): show(_upper_rider,'syntetic');
Constraints
Dependent geometrical parameters
trail constraint
> an: normal trail
> lx: offset
> normal_trail := lx=Rf*sin(epsilon)-an;
wheelbase constraint
> origin( translate(lf*cos(mu),0,-lf*sin(mu)-Rr) * rotate('Y',epsilon) * translate(l23+lx,0,lz) * rotate('Y',-epsilon) * translate(-p,0,Rf) ):
O1 := project(%,ground): #show(OO);
> wheelbase := solve( {comp_X(O1),comp_Z(O1),lx = Rf*s(epsilon)-an},{lx,lz,l23});
Rear tire kinematics
> TWR,rear wheel rolling frame
TR: rear wheel no-rolling frame
rear wheel camber angle
> phi[R](t) := linearize( arcsin(TWR[3,2]) );
> rear_spin_axis:= project(make_VECTOR(TWR,0,1,0),ground): #show(rear_spin_axis);wheel spin axis (unit vector)
> cz := make_VECTOR(ground,0,0,1): #show(cz);direction of travel of the rear wheel (unit vector)
> sr := cross_prod(rear_spin_axis,cz): show(sr);
> if simplify(sqrt(dot_prod(sr,sr))) <> 1 then ERROR("the vector sr has an incorrect magnitude"); end if;check magnitude
rear wheel yaw angle
> psi[R](t) := linearize( arctan( comp_Y(sr),comp_X(sr) ));
SAE reference frame (x axis in agreement with sr)
> TSAEr := rotate('Z','psi[R](t)');
rear wheel contact point
> rho[R]:toroid radius
POSITION
Cr := make_POINT (TR, 0 , rho[R]*'phi[R](t)' , Rr)geometry
+ make_VECTOR(TR, 0, zeta[R](t), 0):carcass deflection
> show(Cr,'synthetic');
SPEED
VCr := project(velocity(Cr),TSAEr):show(VCr,'syntetic');
rear wheel sideslip angle
> lambda[R] := linearize( -arctan((comp_Y(VCr)-vyR(t))/comp_X(VCr)) );
no longitudinal slipping constraint test for the rear wheel
> make_POINT(TWR,comp_X(project(Cr,TWR)),comp_Y(project(Cr,TWR)),comp_Z(project(Cr,TWR))):
if simplify(comp_X(eulerian_velocity(%)))<>0 then ERROR("the wheel is slipping"); end if;
Front tire kinematics
front wheel camber angle
> phi[F](t) := linearize( arcsin(TWF[3,2]) );
> front_spin_axis := project(make_VECTOR(TWF,0,1,0),ground): #show(front_spin_axis);
> sf := cross_prod(front_spin_axis,cz): show(sf);
> if simplify(sqrt(dot_prod(sf,sf))) <> 1 then ERROR("the vector sf has an incorrect magnitude"); end if:
front yaw angle
> psi[F](t) := linearize( arctan( comp_Y(sf),comp_X(sf) ));
SAE reference frame for the front tire (x axis in agreement with sf)
> TSAEf := rotate('Z','psi[F](t)');
front wheel contatc point
> rho[F]:toroid radius
POSITION
Cf := make_POINT(TF, 0, rho[F]*'phi[F](t)' , Rf ) + make_VECTOR(TF, 0, zeta[F](t),0):
show(Cf,'synthetic');
SPEED
VCf := project(velocity(Cf),TSAEf): show(VCf,'syntetic');
front wheel sideslip angle
> lambda[F] := linearize( -arctan((comp_Y(VCf)-vyF(t))/comp_X(VCf)) ):
'lambda[F]' = collect(simplify(subs(wheelbase,lambda[F])),diff(small_vars,t,t) union diff(small_vars,t) union small_vars );
no longitudinal slipping constraint test for the front wheel
> make_POINT(TWF,comp_X(project(Cf,TWF)),comp_Y(project(Cf,TWF)),comp_Z(project(Cf,TWF))):
if simplify(comp_X(eulerian_velocity(%)),trig)<>0 then ERROR("the wheel is slipping"); end if;
Forces
gravity acceleration
> _gravity := make_VECTOR(ground,0,0,g): show(%,'syntetic');
tire forces
rear tire forces and torques
tire forces are applied on the center of the contact patch
the position of contact patch is calculated according to 1) the tire geometry and 2) the tire elastic deformation
> make_VECTOR(TSAEr, Sr,YrS,-Nr):
rear_tire_forces := make_FORCE (%, Cr, rear_wheel): show(rear_tire_forces,'synthetic');
the normal pressure and shear stress distribution move the actual acting point of tire forces from the center of the contact patch.
This effect is taken into account by introducing rolling resistance torque (for normal pressure) and yaw torque (for lateral shear stresses)
> make_VECTOR(TSAEr, 0, Mry, Mrz):
rear_tire_torque := make_TORQUE(%, rear_wheel): show(rear_tire_torque,'synthetic');
front tire forces and torques
> make_VECTOR(TSAEf, Sf,YfS,-Nf):
front_tire_forces := make_FORCE (%, Cf, front_wheel): show(front_tire_forces,'synthetic');
make_VECTOR(TSAEf, 0,Mfy,Mfz):
front_tire_torque := make_TORQUE(%,front_wheel): show(front_tire_torque,'synthetic');
aerodynamics forces
> CA := make_POINT(TRC, -B+ba,0,H-ha): show(CA,'syntetic');aerodynamic center
> FD, FL: drag and lift force
> make_VECTOR(TRC, -FD, 0, FA):
aerodynamic_force := make_FORCE(%, CA, motorcycle): show(aerodynamic_force,'syntetic');
rider steering torque and steering damper
> tau(t), rider active torque
> C[delta], motorcycle steering damping coefficient
> Cb[delta],rider steering damping coefficient
> MX,MY; reactive torques
active steering torque(betwwen front and rear frame)
> make_VECTOR(T4, 0,0, tau(t)):
steering_torque := make_TORQUE(%, front_chassis, motorcycle): show(steering_torque,'syntetic');
motorcycle steering damper torque(between front and rear frame)
> make_VECTOR(T4, MX, MY, tau(t)-C[delta]*diff(delta(t),t)):
steering_damper := make_TORQUE(%, front_chassis, motorcycle):
show(steering_damper, 'syntetic');
rider steering damper torque(between rider and front frame)
> angular_velocity(T4) - angular_velocity(Tup):
omegaZP := comp_Z(project( % , T4));
> make_VECTOR(T4, 0, 0, -Cb[delta]*omegaZP):
body_steering_torque := make_TORQUE(%, front_chassis, upper_rider):
show(body_steering_torque,'syntetic');
rider-saddle interaction
> SP := origin(TS):
> make_VECTOR(TS, FX, -kyp*yp(t)-cyp*diff(yp(t),t), FZ):
saddle_force := make_FORCE(%, SP, rider, motorcycle): show(saddle_force,'syntetic');
> make_VECTOR(TS, -kRxp*phi[P](t)-cRxp*diff(phi[P](t),t), MY, -kRzp*psi[P](t)-cRzp*diff(psi[P](t),t)):
saddle_torque := make_TORQUE(%, rider,motorcycle): show(saddle_torque,'syntetic');
fork flexibility
> make_VECTOR(T60, -kFFx*alpha[F](t)-cFFx*diff(alpha[F](t),t), MY, -kFFz*beta[F](t)-cFFz*diff(beta[F](t),t)):
fork_flex_torque := make_TORQUE(%,flex_fork,front_chassis): show(fork_flex_torque,'synthetic');
swingarm flexibility
> make_VECTOR(T70, -kRFFx*beta[R](t)-cRFFx*diff(beta[R](t),t), MY, -kRFFz*alpha[R](t)-cRFFz*diff(alpha[R](t),t)):
swingarm_flex_torque := make_TORQUE(%,flex_swingarm,motorcycle): show(swingarm_flex_torque,'synthetic');
Newton-Euler's equations
whole vehicle's (rider included) eqns
> whole_vehicle := {
motorcycle,
rear_wheel,_rear_wheel,
front_wheel,_front_wheel,
front_chassis,_front_chassis,
flex_fork,_flex_fork,
flex_swingarm,_flex_swingarm,
rider, upper_rider,_upper_rider,
rear_tire_forces, rear_tire_torque, front_tire_forces, front_tire_torque,
steering_torque, steering_damper, body_steering_torque,
fork_flex_torque, swingarm_flex_torque,
saddle_force, saddle_torque, aerodynamic_force
};
translational equationsfor the whole vehicle
> newton_equations(whole_vehicle,true):
eqns_T_motorcycle := project(%, TRC):
BODIES:
1 - rider
2 - _flex_fork
3 - _rear_wheel
4 - flex_fork
5 - front_wheel
6 - motorcycle
7 - rear_wheel
8 - upper_rider
9 - _flex_swingarm
10 - _front_chassis
11 - _front_wheel
12 - _upper_rider
13 - flex_swingarm
14 - front_chassis
EXTERNAL FORCES:
15 - aerodynamic_force
16 - front_tire_forces
17 - rear_tire_forces
WARNING: the following objects will not appear in the equations of motion:
18 - saddle_force
19 - saddle_torque
20 - steering_damper
21 - steering_torque
22 - fork_flex_torque
23 - front_tire_torque
24 - rear_tire_torque
25 - body_steering_torque
26 - swingarm_flex_torque
> eqns_T_motorcycle[comps] := map(collect, %[comps], diff(small_vars,t,t) union diff(small_vars,t) union small_vars ):
> direction := [X,Y,Z]:
for i from 1 to 3 do printf(direction[i]); eqns_T_motorcycle[comps][i,1]=0; end;
X
Y
Z
rotational equationsfor the whole vehicle
> TEQ := TRC * translate(-B,0,H):
euler_equations(whole_vehicle, origin(TEQ), true):
eqns_R_motorcycle := project(%, TEQ):
BODIES:
1 - rider
2 - _flex_fork
3 - _rear_wheel
4 - flex_fork
5 - front_wheel
6 - motorcycle
7 - rear_wheel
8 - upper_rider
9 - _flex_swingarm
10 - _front_chassis
11 - _front_wheel
12 - _upper_rider
13 - flex_swingarm
14 - front_chassis
EXTERNAL FORCES:
15 - aerodynamic_force
16 - front_tire_forces
17 - rear_tire_forces
EXTERNAL TORQUES:
18 - front_tire_torque
19 - rear_tire_torque
WARNING: the following objects will not appear in the equations of motion:
20 - saddle_force
21 - saddle_torque
22 - steering_damper
23 - steering_torque
24 - fork_flex_torque
25 - body_steering_torque
26 - swingarm_flex_torque
> eqns_R_motorcycle[comps] := map(collect, simplify(expand(subs(normal_trail,eqns_R_motorcycle[comps]))), diff(small_vars,t,t) union diff(small_vars,t) union small_vars ):
> for i from 1 to 3 do printf(direction[i]); eqns_R_motorcycle[comps][i,1]=0; end;
X
Y
Z
front frame eqns
> front_assembly := {steering_torque, steering_damper, body_steering_torque, front_chassis, flex_fork, _flex_fork, fork_flex_torque, front_wheel, front_tire_forces,front_tire_torque};
rotational equationsfor the front frame
> euler_equations( front_assembly ,origin(T4),true):
eqns_R_front_assembly := project(%,T4):
BODIES:
1 - _flex_fork
2 - flex_fork
3 - front_wheel
4 - front_chassis
EXTERNAL FORCES:
5 - front_tire_forces
EXTERNAL TORQUES:
6 - steering_damper
7 - steering_torque
8 - front_tire_torque
9 - body_steering_torque
WARNING: the following objects will not appear in the equations of motion:
10 - fork_flex_torque
> eqns_R_front_assembly[comps] := map(collect, simplify(expand(subs(normal_trail,eqns_R_front_assembly[comps]))), diff(small_vars,t,t) union diff(small_vars,t) union small_vars ):
> for i from 1 to 3 do printf(direction[i]); eqns_R_front_assembly[comps][i,1]=0; end;
X
Y
Z
rider eqns
> rider_assembly := {rider,upper_rider,_upper_rider, saddle_force, saddle_torque, steering_torque, body_steering_torque};
translational equationsfor the suspended rider
> newton_equations(rider_assembly,true):
eqns_T_rider := project(%,TS):
eqns_T_rider[comps]:= map(collect,%[comps],diff(small_vars,t,t) union diff(small_vars,t) union small_vars ):
> for i from 1 to 3 do printf(direction[i]); eqns_T_rider[comps][i,1]=0; end;
BODIES:
1 - rider
2 - upper_rider
3 - _upper_rider
EXTERNAL FORCES:
4 - saddle_force
WARNING: the following objects will not appear in the equations of motion:
5 - saddle_torque
6 - steering_torque
7 - body_steering_torque
X
Y
Z
rotationalequations for the suspended rider
> euler_equations(rider_assembly,SP,true):
eqns_R_rider := project(%,TS):
eqns_R_rider[comps]:= map(collect,%[comps], diff(small_vars,t,t) union diff(small_vars,t) union small_vars ):
for i from 1 to 3 do printf(direction[i]); eqns_R_rider[comps][i,1]=0; end;
BODIES:
1 - rider
2 - upper_rider
3 - _upper_rider
EXTERNAL FORCES:
4 - saddle_force
EXTERNAL TORQUES:
5 - saddle_torque
6 - body_steering_torque
WARNING: the following objects will not appear in the equations of motion:
7 - steering_torque
X
Y
Z
fork eqns
rotational equations for the flexible fork
> fork_assembly:={fork_flex_torque, flex_fork, front_wheel, front_tire_forces, front_tire_torque};
> euler_equations(fork_assembly,origin(T60),true):
eqns_R_flex_fork := project(%,T60):
BODIES:
1 - flex_fork
2 - front_wheel
EXTERNAL FORCES:
3 - front_tire_forces
EXTERNAL TORQUES:
4 - fork_flex_torque
5 - front_tire_torque
> eqns_R_flex_fork[comps]:= map(collect,simplify(map(expand,(eqns_R_flex_fork[comps]))), diff(small_vars,t,t) union diff(small_vars,t) union small_vars ):
for i from 1 to 3 do printf(direction[i]); eqns_R_flex_fork[comps][i,1]=0; end;
X
Y
Z
swingarm eqns
rotational equations for the flexible swingarm
> swingarm_assembly:={rear_tire_forces, rear_tire_torque, rear_wheel, flex_swingarm, swingarm_flex_torque};
> euler_equations(swingarm_assembly,origin(T70),true):
eqns_R_flex_swingarm := project(%,T70):
eqns_R_flex_swingarm[comps]:= map(collect,%[comps], diff(small_vars,t,t) union diff(small_vars,t) union small_vars ):
for i from 1 to 3 do printf(direction[i]); eqns_R_flex_swingarm[comps][i,1]=0; end;
BODIES:
1 - rear_wheel
2 - flex_swingarm
EXTERNAL FORCES:
3 - rear_tire_forces
EXTERNAL TORQUES:
4 - rear_tire_torque
5 - swingarm_flex_torque
X
Y
Z
collection and organization of the equations of motion
> mot_eqns := [
whole motorcycle
comp_Y(eqns_T_motorcycle),y
comp_Z(eqns_R_motorcycle),psi
comp_X(eqns_R_motorcycle),phi
front chassis
comp_Z(eqns_R_front_assembly),delta
rider
comp_Y(eqns_T_rider),yP
comp_X(eqns_R_rider),phiP
fork
comp_X(eqns_R_flex_fork),alphaF
comp_Z(eqns_R_flex_fork),betaF
swingarm
comp_Z(eqns_R_flex_swingarm),alphaR
comp_X(eqns_R_flex_swingarm) betaR
]:
variables
> qq := [y(t),psi(t),phi(t),delta(t),yp(t),phi[P](t),alpha[F](t),beta[F](t),alpha[R](t),beta[R](t),zeta[R](t),zeta[F](t)];
Tires equations
vertical loads
vertical loads on the tires, with speed <>0
> eqns := map(evalm,subs(seq(small_vars[i]=0,i=1..nops(small_vars)), {comp_Z(eqns_T_motorcycle),comp_Y(eqns_R_motorcycle)})):
> tire_loads := simplify(subs(wheelbase,solve(eqns,{Nr,Nf})));
>
rear tire
force due to the carcass elasticity
> YrC := zeta[R](t)*kyr + Nr*'phi[R](t)';
force due to the contacty patch sliding
> Cr1,Cr2;non-dimensionl cornering and sideslip tire stiffnesses
YrS := (Cr1*'lambda[R]' + Cr2*'phi[R](t)') * Nr;
equilibrium equation
> 'h1' = V*simplify('YrC'-'YrS');
h1 := collect(collect(rhs(%),diff(small_vars,t,t) union diff(small_vars,t) union small_vars ),Nr);
tire torques
> Mry := Nr*ur;rolling friction
> Mrz := (CRsa*'lambda[R]' + mrr*'phi[R](t)') * Nr;
front tire
force due to the carcass elasticity
> YfC := zeta[F](t)*kyf + Nf*'phi[F](t)';
force due to the contact patch sliding
> Cf1,Cf2;non-dimensionl cornering and sideslip tire stiffnesses
YfS := (Cf1*'lambda[F]' + Cf2*'phi[F](t)') * Nf;
equilibrium equation
> 'h2' = V*simplify('YfC'-'YfS');
h2 := collect(subs(normal_trail, collect(rhs(%),diff(small_vars,t,t) union diff(small_vars,t) union small_vars )),Nf);
tire torques
> Mfy := Nf*uf;rolling friction
> Mfz := (CFsa*'lambda[F]' + mff*'phi[F](t)') * Nf;
First order formulation
assembly equations and variables
> all_vars := [y(t), psi(t), phi(t), delta(t), yp(t), phi[P](t), alpha[F](t), beta[F](t), alpha[R](t), beta[R](t), zeta[R](t), zeta[F](t)];
all_eqns := [ op(mot_eqns), h1, h2 ]:
reduction to 1st order
> (xx1,state_eqns) := first_order(all_eqns,all_vars,t,flag):
check imposed and calculated variables
> xx := [y[dot](t),psi(t),phi(t),delta(t),yp(t),phi[P](t),alpha[F](t),beta[F](t),alpha[R](t),beta[R](t),zeta[R](t),zeta[F](t), psi[dot](t),phi[dot](t),delta[dot](t),yp[dot](t),phi[P][dot](t),alpha[F][dot](t),beta[F][dot](t),alpha[R][dot](t),beta[R][dot](t)];
if convert(xx,set) minus convert(xx1,set) <> {} then ERROR("WRONG VARIABLES SET") end if;
> nops(state_eqns);
for i from 1 to nops(state_eqns) do i,xx[i];
collect (state_eqns[i],xx union diff(xx,t));
od;
state space formulation
matrix form A*Xdot = B*(X-X0) + D*(U-U0)
> A0:='A0':AA:='AA':BB:='BB':
uu := [tau(t),vyR(t),vyF(t)]; driving vector: steering torque, shaker lateral speed at rear and front wheel
> (A0,AA,BB) := implicit_state_space(state_eqns,t,xx,uu);
> 'A0'=evalm(A0); 'AA'=evalm(AA); 'BB'=evalm(BB);
Dependent data
> Ifa := 'Ifa': Ira := 'Ira':
> tire_loads:=solve(tire_loads,{Nf,Nr}):
dependent_params := tire_loads union wheelbase union {mfy='Ifa'/Rf^2,mry='Ira'/Rr^2};
> tire_loads_static:=subs(ur=0,uf=0,tire_loads);
saving
> #save(all_eqns,h1,h2,A0,AA,BB,state_eqns,xx,uu,dependent_params,tire_loads,"flex-motorcycle-model.m"):
Stability Analysis
geometry, mass, inertia and tyre parameters
> motorcycle_data := {Sr = 68.9, Sf = -68.9, ba = .5, ha = 1, FA = 0, FD = 0, M = 149.30, B = .560, H = .469, IX = 11.624, IZ = 39.568, C[delta] = 0, Cb[delta] = 0, mf = 17.15, Ifx = 1.859, Ifz = .3540, Cfxz = -.780e-2, e = -.8048e-3, f = -.2707, lx6 = 0, lz6 = 0, m6 = 17.15, I6x = 1.859, I6z = .3540, C6xz = .780e-2, gx6 = -.8048e-3, gz6 = -.2707, g = 9.81, CXZ = -1.025, p = 1.381, mu = -.1404277777, epsilon = .4954222221, an = .810e-1, lf = .5096, lx7 = .5096, lz7 = 0, m7 = 48.85, I7x = 1, I7z = 2.45, C7xz = 0, gx7 = -.2966, gz7 = -.46e-1, mrw = 10.15, Ird = .2734, Ira = .4597, Rr = .291, rho[R] = .8e-1, kyr = 57930, Cr1 = 15.5274, Cr2 = 1.29759, ur = .22e-1, mrr = .4e-1, CRsa = -.4458828e-1, mfw = 8.65, Ifd = .1897, Ifa = .3126, Rf = .271, rho[F] = .5e-1, kyf = 57930, Cf1 = 15.5274, Cf2 = 1.29759, uf = .22e-1, mff = .4e-1, CFsa = -.4458828e-1, mp = 75, Ipx = 8.56, Ipz = 3.03, Cpxz = 0, mpUP = 45, IpxUP = 4.75, IpzUP = 1.50, CpxzUP = 0, bs = .55, hs = .80, bp = .64, hp = 1.06, bpUP = .64, hpUP = 1.20}:
numerical state space, with flex parameters to be defined
> A0N := convert(evalf(subs(dependent_params,motorcycle_data,evalm(A0))),Matrix):
> AAN := convert(evalf(subs(dependent_params,motorcycle_data,evalm(AA))),Matrix):
flex parameters
> rigid_params:=[
> kFFx=1e8,cFFx=0,
kFFz=1e8,cFFz=0,
> kRFFx=1e8,cRFFx=0,
> kRFFz=1e8,cRFFz=0,
> kyp=1e8,cyp=0,
kRxp=1e8,cRxp=0]:
numerical state space, when flexibilities are not considered
> A0Nr := subs(rigid_params,A0N):
> AANr := subs(rigid_params,AAN):
speed vector
> speed := [ seq(1.*i,i=1..60)]:
> np := nops(speed):
nn := nops(xx):
rigid vs flex fork (bending)
> flex_fork_params:=[
> kFFx=23e3,cFFx=0,
kFFz=1e8,cFFz=0,
> kRFFx=1e8,cRFFx=0,
> kRFFz=1e8,cRFFz=0,
> kyp=1e8,cyp=0,
kRxp=1e8,cRxp=0]:
numerical state space
> A0Nf := subs(flex_fork_params,A0N):
> AANf := subs(flex_fork_params,AAN):
eigenvalues
> eigen_vals := 'eigen_vals':
for i from 1 to np do
(eigen_vals_r[i],eigen_vect_r[i]) := LinearAlgebra[Eigenvectors](subs(V=speed[i],eval(AANr)),subs(V=speed[i],eval(A0Nr)));
(eigen_vals_f[i],eigen_vect_f[i]) := LinearAlgebra[Eigenvectors](subs(V=speed[i],eval(AANf)),subs(V=speed[i],eval(A0Nf)));
od:
root-locus
> Xrange := -25..5: Yrange := 0..70:
> pr:= plot( [seq( [[Re(eigen_vals_r[n][i]),Im(eigen_vals_r[n][i])] $n=1..np ] ,i=1..nn)], real = Xrange, imag = Yrange, style=point, symbol=circle, color=blue):
pf:= plot( [seq( [[Re(eigen_vals_f[n][i]),Im(eigen_vals_f[n][i])] $n=1..np ] ,i=1..nn)], real = Xrange, imag = Yrange, style=point, symbol=circle, color=red):
plots[display](pr,pf);
rigid vs flex swing (torsion)
> flex_swing_params:=[
> kFFx=1e8,cFFx=0,
kFFz=1e8,cFFz=0,
> kRFFx=10e3,cRFFx=0,
> kRFFz=1e8,cRFFz=0,
> kyp=1e8,cyp=0,
kRxp=1e8,cRxp=0]:
numerical state space
> A0Nf := subs(flex_swing_params,A0N):
> AANf := subs(flex_swing_params,AAN):
eigenvalues
> eigen_vals := 'eigen_vals':
for i from 1 to np do
#(eigen_vals_r[i],eigen_vect_r[i]) := LinearAlgebra[Eigenvectors](subs(V=speed[i],eval(AANr)),subs(V=speed[i],eval(A0Nr)));
(eigen_vals_f[i],eigen_vect_f[i]) := LinearAlgebra[Eigenvectors](subs(V=speed[i],eval(AANf)),subs(V=speed[i],eval(A0Nf)));
od:
root-locus
> Xrange := -25..5: Yrange := 0..70:
> pr:= plot( [seq( [[Re(eigen_vals_r[n][i]),Im(eigen_vals_r[n][i])] $n=1..np ] ,i=1..nn)], real = Xrange, imag = Yrange, style=point, symbol=circle, color=blue):
pf:= plot( [seq( [[Re(eigen_vals_f[n][i]),Im(eigen_vals_f[n][i])] $n=1..np ] ,i=1..nn)], real = Xrange, imag = Yrange, style=point, symbol=circle, color=red):
plots[display](pr,pf);
rigid vs flex swing (bending)
> flex_swing_params:=[
> kFFx=1e8,cFFx=0,
kFFz=1e8,cFFz=0,
> kRFFx=1e8,cRFFx=0,
> kRFFz=20e3,cRFFz=0,
> kyp=1e8,cyp=0,
kRxp=1e8,cRxp=0]:
numerical state space
> A0Nf := subs(flex_swing_params,A0N):
> AANf := subs(flex_swing_params,AAN):
eigenvalues
> eigen_vals := 'eigen_vals':
for i from 1 to np do
#(eigen_vals_r[i],eigen_vect_r[i]) := LinearAlgebra[Eigenvectors](subs(V=speed[i],eval(AANr)),subs(V=speed[i],eval(A0Nr)));
(eigen_vals_f[i],eigen_vect_f[i]) := LinearAlgebra[Eigenvectors](subs(V=speed[i],eval(AANf)),subs(V=speed[i],eval(A0Nf)));
od:
root-locus
> Xrange := -25..5: Yrange := 0..70:
> pr:= plot( [seq( [[Re(eigen_vals_r[n][i]),Im(eigen_vals_r[n][i])] $n=1..np ] ,i=1..nn)], real = Xrange, imag = Yrange, style=point, symbol=circle, color=blue):
pf:= plot( [seq( [[Re(eigen_vals_f[n][i]),Im(eigen_vals_f[n][i])] $n=1..np ] ,i=1..nn)], real = Xrange, imag = Yrange, style=point, symbol=circle, color=red):
plots[display](pr,pf);
rigid vs vibrating rider
> flex_rider_params:=[
> kFFx=1e8,cFFx=0,
kFFz=1e8,cFFz=0,
> kRFFx=1e8,cRFFx=0,
> kRFFz=1e8,cRFFz=0,
> kyp=1e8,cyp=0,
kRxp=800,cRxp=80]:
numerical state space
> A0Nf := subs(flex_rider_params,A0N):
> AANf := subs(flex_rider_params,AAN):
eigenvalues
> eigen_vals := 'eigen_vals':
for i from 1 to np do
#(eigen_vals_r[i],eigen_vect_r[i]) := LinearAlgebra[Eigenvectors](subs(V=speed[i],eval(AANr)),subs(V=speed[i],eval(A0Nr)));
(eigen_vals_f[i],eigen_vect_f[i]) := LinearAlgebra[Eigenvectors](subs(V=speed[i],eval(AANf)),subs(V=speed[i],eval(A0Nf)));
od:
root-locus
> Xrange := -25..5: Yrange := 0..70:
> pr:= plot( [seq( [[Re(eigen_vals_r[n][i]),Im(eigen_vals_r[n][i])] $n=1..np ] ,i=1..nn)], real = Xrange, imag = Yrange, style=point, symbol=circle, color=blue):
pf:= plot( [seq( [[Re(eigen_vals_f[n][i]),Im(eigen_vals_f[n][i])] $n=1..np ] ,i=1..nn)], real = Xrange, imag = Yrange, style=point, symbol=circle, color=red):
plots[display](pr,pf);
>