kill(all);

The rollbar is loaded at the C node by the P and the H forces.

Rollbar intrusion has to be monitored at the point B, 
as long as the reaction forces at the rollbar supports A and E.
Lateral and longitudinal displacements are requested at C, 
and the longitudinal displacement at D.

s is a curvilinear abscissa that spans the beam axis from point to point, 
e.g. from A to B, from B to C, etc.

/home/ebertocchi/_shared/ffcd24/002_rollbar/rollbar_ip_oop_dimensions_and_loads.png

linear auxiliary (shape) functions in the [0,1] interval;

f(0)=1, f(1)=0
g(0)=0, g(1)=1

such that the linear interpolation along the xi∈[0,1] domain may be cast as

a*f(xi) + b*g(xi)

f(pippo) := 1 - pippo;

g(clarabella) := clarabella;

Accumulation variable for the strain energy

U:0;

## Bending and torsional moments 

both in plane and out of plane

### AB beam segment

Moment components at nodes

M1Ap :  fullratsimp( 0*P + 0*YA + 0*XA     - 1*PsiA + 0*F );
M2Ap :  fullratsimp( 0*H + 0*ZA + 1*ThetaA + 0*PhiA + 0*I );
MtAp :  fullratsimp( 0*H + 0*ZA + 0*ThetaA + 1*PhiA + 0*I );

M1Bm :  fullratsimp( 0*P + 0*YA + h*XA     - 1*PsiA + 0*F );
M2Bm :  fullratsimp( 0*H + h*ZA + 1*ThetaA + 0*PhiA + 0*I );
MtBm :  fullratsimp( 0*H + 0*ZA + 0*ThetaA + 1*PhiA + 0*I );

Moment components along the beam segment

M1 : fullratsimp( M1Ap*f(s/h) + M1Bm*g(s/h));
M2 : fullratsimp( M2Ap*f(s/h) + M2Bm*g(s/h));
Mt : fullratsimp( MtAp*f(s/h) + MtBm*g(s/h));

Accumulate on "U" the strain energy due to the AB segment deformation

U : U + integrate(    
     M1^2/2/E/J
    +M2^2/2/E/J
    +Mt^2/2/G/Kt,
    s,0,h
)$

### BC beam segment

Moment components at nodes

M1Bp : M1Bm;
M2Bp : fullratsimp( 0*H + h*ZA + 1*ThetaA + 0*PhiA + 0*I);
MtBp : fullratsimp( 0*H + 0*ZA + 0*ThetaA + 1*PhiA + 0*I);

M1Cm : fullratsimp( 0*P + 0*YA + a*XA -PsiA - (a-h)*F);
M2Cm : fullratsimp( 0*H + a*ZA + 1*ThetaA + 0*PhiA + 0*I);
MtCm : fullratsimp( 0*H + 0*ZA + 0*ThetaA + 1*PhiA + 0*I);

Moment components along the beam segment

M1 : fullratsimp( M1Bp*f(s/(a-h)) + M1Cm*g(s/(a-h)));
M2 : fullratsimp( M2Bp*f(s/(a-h)) + M2Cm*g(s/(a-h)));
Mt : fullratsimp( MtBp*f(s/(a-h)) + MtCm*g(s/(a-h)));

U : U + integrate(    
     M1^2/2/E/J
    +M2^2/2/E/J
    +Mt^2/2/G/Kt,
    s,0,a-h
)$

### CD beam segment

Moment components at nodes

M1Cp : M1Cm ;
M2Cp : fullratsimp( 0*H + 0*ZA + 0*ThetaA - 1*PhiA + 0*I);
MtCp : fullratsimp( 0*H + a*ZA + 1*ThetaA + 0*PhiA + 0*I);

M1Dm : fullratsimp( 0*P - b*YA + a*XA -PsiA - (a-h)*F);
M2Dm : fullratsimp(-b*H + b*ZA + 0*ThetaA - 1*PhiA + 0*I);
MtDm : fullratsimp( 0*H + a*ZA + 1*ThetaA + 0*PhiA + 0*I);

Moment components along the beam segment

M1 : fullratsimp( M1Cp*f(s/b) + M1Dm*g(s/b));
M2 : fullratsimp( M2Cp*f(s/b) + M2Dm*g(s/b));
Mt : fullratsimp( MtCp*f(s/b) + MtDm*g(s/b));

U : U + integrate(    
     M1^2/2/E/J
    +M2^2/2/E/J
    +Mt^2/2/G/Kt,
    s,0,b
)$

### DE beam segment

Moment components at nodes

M1Dp : M1Dm;
M2Dp : fullratsimp( 0*H - a*ZA - 1*ThetaA + 0*PhiA + 0*I);
MtDp : fullratsimp(-b*H + b*ZA + 0*ThetaA - 1*PhiA + 0*I);

M1Em : fullratsimp( a*P - b*YA + 0*XA -PsiA + h*F);
M2Em : fullratsimp(-a*H + 0*ZA - 1*ThetaA + 0*PhiA - a*I);
MtEm : fullratsimp(-b*H + b*ZA + 0*ThetaA - 1*PhiA + 0*I);

Moment components along the beam segment

M1 : fullratsimp( M1Dp*f(s/a) + M1Em*g(s/a));
M2 : fullratsimp( M2Dp*f(s/a) + M2Em*g(s/a));
Mt : fullratsimp( MtDp*f(s/a) + MtEm*g(s/a));

U : U + integrate(    
     M1^2/2/E/J
    +M2^2/2/E/J
    +Mt^2/2/G/Kt,
    s,0,a
)$

Now I should have collected all the strain energy contributions

U : expand(U)$

I observe the uncoupled nature of symmetric and skew-symmetric loads
in the strain energy expression.
The matrix of the quadratic form is 2x2 block defined, and the out-of-diagonal
blocks are in fact  zero.

hessian(U , [ P,F,XA,YA,PsiA , H,I,ZA,ThetaA,PhiA])$

with the aid of the second Castigliano theorem, I evaluate notable deflections
and rotations

dP : diff(U,P)$

dF : diff(U,F)$

dH : diff(U,H)$

dI : diff(U,I)$

uA : diff(U,XA)$
vA : diff(U,YA)$
wA : diff(U,ZA)$
thetaA : diff(U,ThetaA)$
  phiA : diff(U,  PhiA)$
  psiA : diff(U,  PsiA)$

deflections at A are monitored in order to tune the clamp reactions such that 
they vanish

I define a list of equations, which are just expressions 
with an "=" sign within; an expression missing the "=" character, when 
treated as an equation, is set to zero, i.e. a "=0" is implicitly 
appended.

equazioni : [ uA=0 , vA=0, wA  , thetaA=0, phiA  , psiA=0 ]$

incognite : [ XA,YA,ZA,ThetaA,PhiA,PsiA];

"solve" is the general [system of] equations solver, "linsolve" is a
specialized implementation dedicated to linear [system of] equations.

sols : linsolve( equazioni, incognite )$

[dP,dF,dH,dI,uA,vA,wA,thetaA,phiA,psiA] :
factor(
ratsimp(
ev(
    [dP,dF,dH,dI,uA,vA,wA,thetaA,phiA,psiA],
    sols,
    F=0, /* I declare the fictitious nature of F and I*/
    I=0,
    eval
)
)
)$

I proceed with some normalizations; I need reference displacements for
in plane and out of plane deflections, I take them from the 
single upright member  substructur (I fictitiously assume that
the directly loaded vertical member is left alone in bearing the load)

M1 : P*s$
U:integrate(M1^2/2/E/J,s,0,a)$
refdP :diff(U,P);
refdH:ratsubst(H,P,refdP);

let's proceed with normalization, and a chain of manipulations

r=b/a aspect ratio of the frame/rollbar 
r<<1 is tall and narrow
r=1 is square
r>>1 is short and wide

q=h/a normalized location of the lateral displacement evaluation
q=0 I sample at A
q=1 I sample at C
q=1/2 I sample at the AC midpoint

I define the "t" adimensional parameter as
t= (G*Kt)/(E*J)

for a solid/hollow circular cross section, 
and isotropic material
G*Kt = E/2/(1+nu) * %pi*(Do^4-Di^4)/32
E*J  = E          * %pi*(Do^4-Di^4)/32/2
t=1/(nu+1)

[dP,dF,dH,dI]$
% / [refdP,refdP,refdH,refdH]$
ratsubst( r*a  , b   ,%)$
ratsubst( q*a  , h   ,%)$
ratsubst( t*E*J, G*Kt,% )$
ratsimp(%)$
[ndP,ndF,ndH,ndI]:%;

I plot the inverse of the deflection ratio between the full frame
and the single upright member, with deflections taken at the loaded point
for in plane loading

kframe / kref = (P/dP) / (P/dPref) = dPref/dP = 1/ndP

ip_kratio:1/ndP$

wxplot2d( 
    ip_kratio , 
    [r,0.001,1000], 
    logx , 
    [legend, "i.p. kframe/kref"],
    [ylabel, "i.p. kframe/kref"]
);

some reference values

[   
    ev(ip_kratio,r=1/5),
    ev(ip_kratio,r=1  ),
    ev(ip_kratio,r=5  )
],numer;

check the behavior for B->C

limit for q approaching 1 of the ndF/ndP ratio

limit( ndF/ndP, q, 1 );

ndF/dnP graph for varying q, and a few r values

ndF/ndP;
[  
    ev(%,r=1/5),
    ev(%,r=1  ),
    ev(%,r=5  )
];
wxplot2d( % , [q,0,1] , [legend, "r=1/5", "r=1","r=5" ] , [ylabel,"dF/dP"]);