flex motorcycle model Madrid.html

Motorcycle Model with Frame Compliance

Roberto Lot & Matteo Massaro 2005 5 17

REFERENCES:

The Significance of Frame Compliance and Rider Mobility on the Motorcycle Stability

Multibody Dynamics 2005 - International Conference on Advances in Computational Multibody Dynamics

> restart:

> mylib :="C://libs":
libname := mylib,libname: with(MBsymba):

> alias(c=cos,s=sin); 'sin' and 'cosine' functions

Motorcycle Model

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:

$small_vars = {y(t), psi(t), phi(t), delta(t), yp(t)...$

Semplificative assumptions

> 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

>

Definition of 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');

>

Wheels

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');

>

Kinematics

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( {Xcomp(O1),Zcomp(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( Ycomp(sr),Xcomp(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,'absolute'),TSAEr):show(VCr,'syntetic');

rear wheel sideslip angle

> lambda[R] := linearize( -arctan((Ycomp(VCr)-vyR(t))/Xcomp(VCr)) );

no longitudinal slipping constraint test for the rear wheel

> make_POINT(TWR,Xcomp(project(Cr,TWR)),Ycomp(project(Cr,TWR)),Zcomp(project(Cr,TWR))):
if Xcomp(velocity(%,'eulerian'))<>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( Ycomp(sf),Xcomp(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,'absolute'),TSAEf): show(VCf,'syntetic');

front wheel sideslip angle

> lambda[F] := linearize( -arctan((Ycomp(VCf)-vyF(t))/Xcomp(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 rear wheel

> make_POINT(TWF,Xcomp(project(Cf,TWF)),Ycomp(project(Cf,TWF)),Zcomp(project(Cf,TWF)));
if Xcomp(velocity(%,'eulerian'))<>0 then ERROR("the wheel is slipping"); end if;

>

Forces

define acceleration due to the gravity

> _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

> rear_tire_forces := make_FORCE (TSAEr, Sr,YrS,-Nr, 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)

> rear_tire_torque := make_TORQUE(TSAEr, 0, Mry, Mrz, rear_wheel, ground): show(rear_tire_torque,'synthetic');

front tire forces and torques

> front_tire_forces := make_FORCE (TSAEf, Sf,YfS,-Nf, Cf, front_wheel): show(front_tire_forces,'synthetic');
front_tire_torque := make_TORQUE(TSAEf, 0,Mfy,Mfz,front_wheel, ground): 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

> aerodynamic_force := make_FORCE(TRC, -FD, 0, FA, 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)

> steering_torque := make_TORQUE(T4, 0,0, tau(t), front_chassis, motorcycle): show(steering_torque,'syntetic');

motorcycle steering damper torque (between front and rear frame)

> steering_damper := make_TORQUE(T4, MX, MY, tau(t)-C[delta]*diff(delta(t),t), front_chassis, motorcycle):
show(steering_damper, 'syntetic');

rider steering damper torque (between rider and front frame)

> angular_velocity(T4) - angular_velocity(Tup):
omegaZP := Zcomp(project( % , T4));

> body_steering_torque := make_TORQUE(T4, 0, 0, -Cb[delta]*omegaZP, front_chassis, upper_rider):
show(body_steering_torque,'syntetic');

>

rider-saddle interaction

> SP := origin(TS):

> saddle_force := make_FORCE(TS, FX, -kyp*yp(t)-cyp*diff(yp(t),t), FZ, SP, rider, motorcycle): show(saddle_force,'syntetic');

> saddle_torque := make_TORQUE(TS, -kRxp*phi[P](t)-cRxp*diff(phi[P](t),t), MY, -kRzp*psi[P](t)-cRzp*diff(psi[P](t),t), rider,motorcycle): show(saddle_torque,'syntetic');

>

fork flexibility

> fork_flex_torque := make_TORQUE(T60, -kFFx*alpha[F](t)-cFFx*diff(alpha[F](t),t), MY, -kFFz*beta[F](t)-cFFz*diff(beta[F](t),t),
flex_fork,front_chassis): show(fork_flex_torque,'synthetic');

>

swingarm flexibility

> swingarm_flex_torque := make_TORQUE(T70, -kRFFx*beta[R](t)-cRFFx*diff(beta[R](t),t), MY, -kRFFz*alpha[R](t)-cRFFz*diff(alpha[R](t),t),flex_swingarm,motorcycle): show(swingarm_flex_torque,'synthetic');

>

Linear 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
};

$whole_vehicle := {flex_fork, motorcycle, aerodynami...$
$whole_vehicle := {flex_fork, motorcycle, aerodynami...$
$whole_vehicle := {flex_fork, motorcycle, aerodynami...$

translational equations for the whole vehicle

> eqns_T_motorcycle := newton_equations(whole_vehicle, TRC, verbose):

BODIES:

  1 - flex_fork

  2 - motorcycle

  3 - _front_chassis

  4 - front_wheel

  5 - upper_rider

  6 - _flex_fork

  7 - _rear_wheel

  8 - _front_wheel

  9 - front_chassis

 10 - rear_wheel

 11 - _flex_swingarm

 12 - rider

 13 - flex_swingarm

 14 - _upper_rider

EXTERNAL FORCES:

 15 - aerodynamic_force

 16 - rear_tire_forces

 17 - front_tire_forces

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

 18 - front_tire_torque

 19 - fork_flex_torque

 20 - body_steering_torque

 21 - swingarm_flex_torque

 22 - steering_torque

 23 - saddle_force

 24 - rear_tire_torque

 25 - steering_damper

 26 - saddle_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;

rotational equations f or the whole vehicle

> TEQ := TRC * translate(-B,0,H):
eqns_R_motorcycle := euler_equations(whole_vehicle, origin(TEQ), TEQ, 'verbose'):

BODIES:

  1 - flex_fork

  2 - motorcycle

  3 - _front_chassis

  4 - front_wheel

  5 - upper_rider

  6 - _flex_fork

  7 - _rear_wheel

  8 - _front_wheel

  9 - front_chassis

 10 - rear_wheel

 11 - _flex_swingarm

 12 - rider

 13 - flex_swingarm

 14 - _upper_rider

EXTERNAL FORCES:

 15 - aerodynamic_force

 16 - rear_tire_forces

 17 - front_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 - fork_flex_torque

 21 - body_steering_torque

 22 - swingarm_flex_torque

 23 - steering_torque

 24 - saddle_force

 25 - steering_damper

 26 - saddle_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;

>

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};

$front_assembly := {flex_fork, front_tire_torque, fo...$

rotational equations for the front frame

> eqns_R_front_assembly := euler_equations( front_assembly ,origin(T4),T4,'verbose'):

BODIES:

  1 - flex_fork

  2 - front_wheel

  3 - _flex_fork

  4 - front_chassis

EXTERNAL FORCES:

  5 - front_tire_forces

EXTERNAL TORQUES:

  6 - front_tire_torque

  7 - body_steering_torque

  8 - steering_torque

  9 - steering_damper

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;

>

rider eqns

> rider_assembly := {rider,upper_rider,_upper_rider, saddle_force, saddle_torque, steering_torque, body_steering_torque};

$rider_assembly := {body_steering_torque, upper_ride...$

translational equations for the suspended rider

> eqns_T_rider := newton_equations(rider_assembly,TS,verbose):
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 - upper_rider

  2 - rider

  3 - _upper_rider

EXTERNAL FORCES:

  4 - saddle_force

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

  5 - body_steering_torque

  6 - steering_torque

  7 - saddle_torque

rotational equations for the suspended rider

> eqns_R_rider := euler_equations(rider_assembly,SP,TS,verbose):
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 - upper_rider

  2 - rider

  3 - _upper_rider

EXTERNAL FORCES:

  4 - saddle_force

EXTERNAL TORQUES:

  5 - body_steering_torque

  6 - saddle_torque

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

  7 - steering_torque

>

fork eqns

rotational equations for the flexible fork

> fork_assembly:={fork_flex_torque, flex_fork, front_wheel, front_tire_forces, front_tire_torque};

$fork_assembly := {flex_fork, front_tire_torque, for...$

> eqns_R_flex_fork := euler_equations(fork_assembly,origin(T60),T60,verbose):

BODIES:

  1 - flex_fork

  2 - front_wheel

EXTERNAL FORCES:

  3 - front_tire_forces

EXTERNAL TORQUES:

  4 - front_tire_torque

  5 - fork_flex_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;

>

swingarm eqns

rotational equations for the flexible swingarm

> swingarm_assembly:={rear_tire_forces, rear_tire_torque, rear_wheel, flex_swingarm, swingarm_flex_torque};

$swingarm_assembly := {swingarm_flex_torque, rear_ti...$

> eqns_R_flex_swingarm := euler_equations(swingarm_assembly,origin(T70),T70,verbose):
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 - swingarm_flex_torque

  5 - rear_tire_torque

>

collection and organization of the equations of motion

> mot_eqns := [
whole motorcycle
Ycomp(eqns_T_motorcycle), y
Zcomp(eqns_R_motorcycle), psi
Xcomp(eqns_R_motorcycle), phi
front chassis
Zcomp(eqns_R_front_assembly), delta
rider
Ycomp(eqns_T_rider), yP
Xcomp(eqns_R_rider), phiP
fork
Xcomp(eqns_R_flex_fork), alphaF
Zcomp(eqns_R_flex_fork), betaF
swingarm
Zcomp(eqns_R_flex_swingarm), alphaR
Xcomp(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)), {Zcomp(eqns_T_motorcycle),Ycomp(eqns_R_motorcycle)})):

> tire_loads := simplify(subs(wheelbase,solve(eqns,{Nr,Nf})));

$tire_loads := {Nr = -(-ha*FD-Mry+bp*g*mp+B*g*M+ba*F...$

>

rear tire

force due to the carcass elasticity

> YrC := zeta[R](t)*kyr + Nr*'phi[R](t)';

force due to the contact 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 equations

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;