Clear the memory, or, kill all the variable definitions

kill( all ) ;

	equilibrium equations for the quarter
	(simplified) ladder frame chassis


eqtX : UA = 0 ;


eqtY : VD = 0 ;


eqtZ : WA + WD - F = 0;


eqrOX : - F * c + WD * b + ThetaD = 0;


eqrOY : PhiA - WA*a + F*a = 0;


eqrOZ : 0 = 0 ;


	define a list of equations and a list of unknowns


eqlst : [ eqtX  , eqtY  , eqtZ ,
          eqrOX , eqrOY , eqrOZ ] ;


unlst : [ UA , WA , PhiA , VD , WD, ThetaD ];


	solve the linear system of equations


linsolve( eqlst , unlst ) ;


	there is one dependent equation, so I'm unable to solve 
	the problem of obtaining all the reaction forces/moments by equilibrium alone


unlst_red : [ UA , WA , PhiA , VD , WD];


	Solution of the system of linear equation and assignemtn to the variables in one step.
	linsolve(eqn,[x_1,..x_n]), globalsolve = true;


linsolve( eqlst , unlst_red ) , globalsolve = true ;


	we now define the cross section resultants for the D->B beam, 
	as a function of the s = 0 ... a arclength coordinate
	directions x and y constitute an in-section reference system.


	traslation equilibrium along direction z 
	N:0;
	traslation equilibrium along direction x and y
	Sx and Sy are the shear forces acting at the shear center (Point C)
	Sx: -VD;
	Sy: -WD;
	N acts at the centroid of section (Point G)

N  :  0 $

Sx : -VD;

Sy : -WD;


	equilibrium with respect to rotation along the (G,x) axis


linsolve(Mfx + WD * s = 0, Mfx ), globalsolve = true ;

linsolve(Mfy - VD * s = 0, Mfy ), globalsolve = true ;

	torque must be evaluated by the equilibrium along the axis z at the shear center.

linsolve(Mt + ThetaD - WD * e = 0, Mt), globalsolve = true ;

	we now evaluate the internal strain energy for the DB beam
	Kt coincides with Jp for circular or hollowed circular cross-section. 

U : integrate ( 
               ( 
                       Jyy * Mfy^2
                   +   Jxx * Mfx^2
                   + 2*J12 * Mfx*Mfy
                )/( 
                    2*E*(Jxx*Jyy-Jxy^2)
                )
            +
                N^2/2/E/A
            +
                ( 
                      chix   * Sx^2
                    + chiy   * Sy^2
                    + chixy  * Sx * Sy
                )/(
                    2*G*A
                )
            +
                Mt^2 / 2 / G/ ( Kt )
 , s , 0 , a );


	The strain energy is now known as a quadratic function of F and ThetaD
	We recur to the Castigliano's second theorem to obtain the rotation at D, 
	component x

thetaD : diff( U , ThetaD )   ;
	
	we now force that rotation to be zero in order to reinstate the kinematic
	compatibility with respect to the skew-symmetry constraint acting in D.

linsolve ( thetaD = 0 , ThetaD ) , globalsolve = true;

	I now redefine the internal strain energy as a function of F alone,
	i.e. I now force a substitution of the expression above within U

U : ev( U )$

	We employ Castigliano theorem again to obtain the deflection of
	the quarter ladder frame chassis at the application point of the external
	force F


delta : diff( U , F )$

	we now evaluate from delta the torsional stiffness (ts) of the chassis;
	the deflection delta of the whole chassis, supported at three
	wheel centers and loaded at the fourth, is four times delta.


ts : fullratsimp(
 (F * (2*c)^2) / ( 4 * delta ) ) ;