kill(all);

equilibrium equations

eqtrx : UA = 0;
eqtry : VD = 0;
eqtrz : -F+WA+WD = 0;

eqrOx : ThetaD +0*VD +b*WD -c*F ;

eqrOy : PhiA + a*F -a*WA + 0*UA ;

eqrOz : 0 = 0 ;

blindly solve the equations

linsolve(
[ eqtrx , eqtry, eqtrz , eqrOx , eqrOy, eqrOz ],
[UA,PhiA, WD,VD,ThetaD]
) , globalsolve=True;

the WA unknown reaction force has been removed from
the unknowns list, thus obtaining a parametric
representation of the reaction forces and moments, 
with WA as a parameter (like F)

to obtain the actual value of WA, the deflection in the 
z direction of the A point has to be evaluated and
zeroed. Castigliano's Theorem fits.

N  : 0;
Sx : -VD;
Sy : -WD;

equilibrium with respect to rotation along the
local Cz axis 

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

equilibrium with respect to rotation along the
local Gx and Gy axes

linsolve( Mx +WD*s =0, Mx), globalsolve=True;

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

Internal energy for the DB beam

U : integrate(
(Jyy*Mx^2+Jxx*My^2+2*Jxy*Mx*My)
 /
(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);

Castigliano's Theorem is now applied to obtain the 
z deflection at A

linsolve(
    diff( U , WA ) = 0,
    WA
), globalsolve=true;

U : ev(U);

I check that all the WA terms hidden within U are now
expressed as functions of F

diff( U , 'WA) ;

the z deflection at the L point is evaluated for the
quarter ladder frame structure, constrained according 
to two skew-symm planes

d : diff(U , F);

dw : diff(4*U , F);