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+e) + 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];

sols : linsolve( eqlst , unlst_red ) ;

option A for assignment:

''unlst_red : ev( unlst_red , sols ) ;

massive assignment of the values that solve the equations to the 
unknowns. The ''unlst_red syntax forces an immediate evaluation of the
symbol, such that the actual assignment is

[UA,WA,PhiA,VD,WD] : [UA=0,WA=((e-c+b)*F+ThetaD)/(e+b),
                      PhiA=-(a*c*F-a*ThetaD)/(e+b),VD=0,
                      WD=(c*F-ThetaD)/(e+b)];

option B for assignment

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

N : 0 $

Mf_eta : VD * s $

Mf_xi : WD * s ;

Mt : ThetaD;

S_eta : WD;

S_xi  : VD$

we now evaluate the internal strain energy for the DB beam

U : integrate ( 
               ( 
                       J_etaeta* Mf_xi^2
                   +    J_xixi * Mf_eta^2
                   + 2*J_xieta * Mf_xi*Mf_eta
                )/( 
                    2*E*(J_xixi*J_etaeta-J_xieta^2) * bending_stiff_factor
                )
            +
                N^2/2/E/A
            +
                ( 
                      chi_xi * S_xi^2
                    + chi_eta * S_eta^2
                    + chi_etaxi * S_xi * S_eta
                )/(
                    2*G*A  * shear_stiff_factor
                )
            +
                Mt^2 / 2 / G/ ( Kt * torsion_stiff_factor )
 , 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 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 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 )
) ;

i copy and paste the section properties, the material properties and
the chassis dimensions

Channel section dimensions:
h : web heigth, measured at midsurface
b : wing width, mesured at midsurface
t : web and wing wall thickness
(e,f) : xi,eta position of the shear center, with respect to the center of gravity
(m,n) : xi,eta position of the center of gravity, with respect to the web midpoint
upsilon: b/h
epsilon: t/h
Chassis dimensions:
a   : half wheelbase
b   : half the distance between the profiles
c   : half the axle track

dim : 
[
    a=1200,
    b=400-m,
    c=800,
    h=80,
    upsilon=40/80,
    epsilon=6/80
];

mat : 
[
    G  = E / 2 / (1 +nu),
    E  = 210000,
    nu = 3/10    
];

sec:
[
    A        = 2*epsilon*h^2*upsilon+epsilon*h^2,
    J_xixi   =(epsilon*h^4+6*epsilon*h^4*upsilon)/12,
    J_etaeta =(2*epsilon*h^4*upsilon^3+epsilon*h^4*upsilon^4)/(6*upsilon+3),
    J_xieta  =0,
    Kt       =(epsilon^3*h^3*(h+2*h*upsilon))/3,
    chi_xi   =(15+78*upsilon+138*upsilon^2+96*upsilon^3+24*upsilon^4)/(20*upsilon^4+80*upsilon^3+80*upsilon^2),
    chi_eta  =(6+72*upsilon+300*upsilon^2+480*upsilon^3+240*upsilon^4)/(180*upsilon^2+60*upsilon+5),
    chi_xieta=0,
    e        = (4*h*upsilon^2+12*h*upsilon^3)/(12*upsilon^2+8*upsilon+1),
    f        = 0  ,
    m        = (h*upsilon^2)/(2*upsilon+1),
    n        = 0
];

some experiments

limit for null torsional stiffness

limit ( ts , torsion_stiff_factor, 0);

ev( 
    % , dim , sec , mat , 
    bending_stiff_factor=1,
    shear_stiff_factor=1,
    infeval 
);

  profile torsional stiffness lost , 
+ profile bending stiffness lost

limit (
    limit ( 
        ts , torsion_stiff_factor, 0
    ),
    bending_stiff_factor, 0
);

  profile torsional stiffness lost , 
+ profile shear stiffness lost

limit (
    limit ( 
        ts , torsion_stiff_factor, 0
    ),
    shear_stiff_factor, 0
);

bending stiffness lost

limit ( 
        ts , bending_stiff_factor, 0
    );

residual chassis stiffness is due to the torsional stiffness
of the profile

shear stiffness lost

limit ( 
        ts , shear_stiff_factor, 0
    );

infinite profile torsional stiffness

limit ( ts , torsion_stiff_factor , infinity );

infinite profile bending stiffness

limit ( ts , bending_stiff_factor , infinity );

infinite profile bending and shear stiffness

limit(
    limit ( 
        ts , bending_stiff_factor , infinity 
    ) , shear_stiff_factor , infinity
);