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 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(xi) := 1-xi;

g(xi) :=   xi;

## Bending and torsional moments, with strain energy
accumulation

both in plane and out of plane

strain energy accumulation variable

U : 0 $

### 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 current beam segment

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

strain energy contribution due to the current beam segment

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 current beam segment

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

strain energy contribution due to the current beam segment

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 current beam segment

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

strain energy contribution due to the current beam segment

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 current beam segment

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

strain energy contribution due to the current beam segment

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 the strain energy along the whole beam

U:expand(U);

it's a quadratic form in the actions

P,F,XA,YA,PsiA , H,I,ZA,ThetaA,PhiA

whose matrix may be derived in maxima 
by observing that if I write the quadratic form as

1/2 * P' * M * P = U

the M matrix is consistent with the Hessian matrix of the second order derivatives

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

## Displacement calculation

Now I have the strain energy, I can proceed with the application of the 2nd Castigliano
theorem
in order to retrieve displacements and rotations

[P,F,XA,YA,PsiA , H,I,ZA,ThetaA,PhiA]

deflections at application points of the external (real or fictitious) forces

dP : diff( U, P )$
dF : diff( U, F );
dH : diff( U, H )$
dI : diff( U, I )$

from this point on, I can stick F and I to the null value, since they were introduced just to 
evaluate deflections through the Castigliano theorem, as in the command above

F:0$
I:0$

Redefine the quantities of interest based on the now declared null value of F and I

[dP,dF,dH,dI,U]:ev([dP,dF,dH,dI,U]);

deflections at the A initially clamped rollbar end

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

impose null deflections/rotations at the clamped end

linsolve(
    [ uA=0 , vA , wA=0  , thetaA=0, phiA  , psiA /*=0 implicitly, if not explicited*/ ],
    [ XA, YA,ZA, ThetaA, PhiA, PsiA ]
), globalsolve=true ;

derive reactions at other clamped end E by considering the equilibrium of the
complete rollbar

linsolve([
        /* add the other equations here*/
            /*eq. tx*/
            /*eq. ty*/
            /*eq. tz*/
            /*eq. rotEx*/
        PhiE+PhiA+0*I+H*b-ZA*b=0 /*eq. rotEy*/
            /*eq. rotEz*/
    ],
    [ /* ... */ PhiE /*...*/ ]
),globalsolve=true;

I redefine quantities of interest based on the now known value of the reactions at A

[dP,dF,dH,dI] : factor(ratsimp(ev([dP,dF,dH,dI])));

I can similarly proceed with bending/torsional moment values at nodes, which
may be also interpolated if needed

factor(ev([M1Cp,M2Cp,MtCp]));
factor(ev([M1Dm,M2Dm,MtDm]));

let's find a suitable normalization to further simplify the obtained expression

reference displacements for in plane and o.o.p. loadings

[dref_ip,dref_oop] : a^3*[P,H]/3/E/J;

normalize deflections accordingly

[ndP,ndF,ndH,ndI] : factor([dP/dref_ip, dF/dref_ip, dH/dref_oop, dI/dref_oop]);

let's try to further reduce the number of parameters by e.g. defining

b = r*a , with r=b/a aspect ratio of the rollbar

if r =1 the rollbar is square,
if r<<1 the rollbar is tall and/or narrow
if r>>1 the rollbar is short and/or wide

h = q*a, with q=h/a a normalize paramater defining where we sample 
the lateral displacement along the directly loaded upright member

if q=1 we sample at the C upper node
if q=0 we sample at the A clamp
if q=1/2 we sample halfway the vertical member

G*Kt = t*E*J, where t = (G*Kt)/(E*J)

t is a function 
    i) of the E/G ratio, which in turn is a function of the poisson ratio
in isotropic materials, and 
    ii) of the cross sectiona shape

if t <<1  the cross section is patologically more deformable in torsion with respect to its
compliance in bending, e.g. open thin walled cross section (I-beams, C-beams) or possibly
closed thin walled cross sections lacking ±45° reinforcements
if t >>1 the cross section in patologicall more deformable in bending with respect to torsion
(e.g. not enough 0° axial fibers)

for the thin walled hollow circular cross section under scrutiny

myt: ( E/2/(1+nu) * (douter^4-dinner^4)*%pi/32 ) / (E*(douter^4-dinner^4)*%pi/64);

myt,nu=0.3;

[ndP,ndF,ndH,ndI]$
ratsubst(r*a,b,%)$
ratsubst(q*a,h,%)$
ratsubst(t*E*J,G*Kt,%)$
factor(%);
[ndP,ndF,ndH,ndI] : % $

now, all the quantities of interest are reduced to functions of normalized
parameters, i.e. ratios between dimensionally consistent quantities

Let's analize the stiffness of the rollbar structure vs. the stiffness of
the reference - single upright member cantilever beam - structure

normalized stiffness for in plane loads
defined as the inverse of the normalized compliance ratio

ndP <- norm. displacement at the loaded point of the rollbar
---
 1  <- norm. displacement at the loaded point of the reference cantilever beam

ipkratio : (1)/(ndP);

let's plot this ratio

wxplot2d(
    ipkratio,
    [ r , 0.001 , 1000 ],
    [logx],
    [ylabel,"ipkratio"]
);

ev(ipkratio,r=1,numer);

let's evaluate the out of plane stiffness ratio

oopkratio: (1)/(ndH);

wxplot2d(
    [
        ev(oopkratio,t=myt*0.01), /* seriously lacking of +-45 deg. fibers*/
        ev(oopkratio,t=myt*0.1), /* lacking of +-45 deg. fibers*/
        ev(oopkratio,t=myt*1), /*isotropic material, hollow circular cross section*/
        ev(oopkratio,t=myt*10), /* Superalloy-Z +-45 deg. fiber*/
        ev(oopkratio,t=myt*100) 
    ],
    [r,0.001,1000],
    [logx],
    [legend,"t/tref = 0.01","t/tref = 0.1","t/tref = 1", "t/tref = 10","t/tref = 100"],
    [ylabel,"oopkratio"]
), nu=0.3;

ev(oopkratio,r=1);
makelist(ev(%,t=tcurrent) , tcurrent , myt*[0.01,0.1,1,10,100]),nu=0.3,eval;

wxplot2d(
    ev([ndH,ndI],t=myt),
    [r,0.001,1000],
    [logx],
    [legend,"ndH","ndI"]
),nu=0.3;