kill(all);

## Auxiliary definitions

I define two auxiliary functions for expressing the moments, where linear
on a branch

f( ascissanormalizzata ) := 1 - ascissanormalizzata ;

g( ascissanormalizzata ) := ascissanormalizzata ;

## Bending and torsional moments

### AB beam segment

Moment components at nodes

monospaced font
iiiiii
mmmmmm
Courier New

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));

Strain energy calculation

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

### 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)));

Strain energy calculation AND ACCUMULATION

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

### 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));

Strain energy calculation AND ACCUMULATION

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

### 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));

Strain energy calculation AND ACCUMULATION

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

## some further checks on the internal action components at nodes

double check all the nodal values!

Fig 2 graphs

l:a-h$

 [ M1A, M1B, M1C, M1D, M1E] :
 [  0 ,  0 ,  0 ,  0 ,  a ]*P
+[  0 ,  0 ,  0 , -b , -b ]*YA
+[  0 ,  h ,  a ,  a ,  0 ]*XA
+[ -1 , -1 , -1 , -1 , -1 ]*PsiA
+[  0 ,  0 , -l , -l , +h ]*F;

comparison with the previously defined M1

fullratsimp(
    [M1A,M1B,M1C,M1D]-[M1Ap,M1Bp,M1Cp,M1Dp]
);

fullratsimp(
    [M1B,M1C,M1D,M1E]-[M1Bm,M1Cm,M1Dm,M1Em]
);

a null difference, WOW!!! 

Figure 3 graphs, bending moment component 2

fullratsimp(
  - [M2Ap,M2Bm,M2Bp,M2Cm,M2Cp,M2Dm,M2Dp,M2Em]
  + [   0,   0,   0,   0,   0,  -b,   0,  -a]*H
  + [   0,  +h,  +h,  +a,   0,  +b,  -a,   0]*ZA
  + [  +1,  +1,  +1,  +1,   0,   0,  -1,  -1]*ThetaA
  + [   0,   0,   0,   0,  -1,  -1,   0,   0]*PhiA
  + [   0,   0,   0,   0,   0,   0,   0,  -a]*I
);

torsional moment

fullratsimp(
  - [MtAp,MtBm,MtBp,MtCm,MtCp,MtDm,MtDp,MtEm]
  + [   0,   0,   0,   0,   0,   0,  -b,  -b]*H
  + [   0,   0,   0,   0,  +a,  +a,  +b,  +b]*ZA
  + [   0,   0,   0,   0,  +1,  +1,   0,   0]*ThetaA
  + [  +1,  +1,  +1,  +1,   0,   0,  -1,  -1]*PhiA
  + [   0,   0,   0,   0,   0,   0,   0,   0]*I
);

## Castigliano theorem application

### spostamenti e rotazioni al punto A

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

### spostamenti richiesti come risultati dell'esercizio

utB : diff(U,F)$
utC : diff(U,P)$
wtC : diff(U,H)$
wtD : diff(U,I)$

## Impongo la compatibilità cinematica al punto A

x     y   z   assi
u     v   w   spostamenti
theta phi psi rotazioni
X     Y   Z   forze
U     V   W   forze
Theta Phi Psi coppie

def_ipRA: 
linsolve(
  [ uA = 0 , vA = 0, psiA = 0 ],
  [ XA ,     YA,     PsiA     ]  
);

def_oopRA:
linsolve(
  [ wA = 0 , thetaA = 0 , phiA = 0 ],
  [ ZA ,     ThetaA ,     PhiA     ]
);

sostituisco entro gli spostamenti di interesse i.p. le espressioni ricavate per le reazioni vincolare i.p.

    [ utB , utC ] 
: 
ev( [ utB , utC ] , def_ipRA )$

ripeto per spostamenti e rotazioni o.o.p.

    [ wtC , wtD ] 
: 
ev( [ wtC , wtD ] , def_oopRA )$

esplicito la natura fittizia delle azioni F e I

    [ utB , utC ] 
: 
ev( [ utB , utC ] , F=0 );

    [ wtC , wtD ] 
: 
ev( [ wtC , wtD ] , I=0 );

[utB,utC] : factor( [utB,utC] );

[wtC,wtD] : factor( [wtC,wtD] );

ingolosito dalla possibilità di semplificare oltre wtC e wtD
provo ad esprimere G*Kt come eta*E*J

ratsubst( 
    ESPR. CON CUI SOSTITUIRE , 
    ESPR. DA CERCARE ,  
    ESPR. SU CUI OPERARE
);

[wtC,wtD] : factor(
    ratsubst( eta*E*J , G*Kt , [wtC,wtD] )
);

descrivo la risposta del rollbar in funzione del rapporto
adimensionale r = b/a e p=h/a

ratsubst( r*a , b , espressione di interesse )
ratsubst( p*a , h , espressione di interesse )

[utB, utC , wtC, wtD] $
ratsubst( r*a , b , % ) $
ratsubst( p*a , h , % ) $
factor( % ) $ 
[utB, utC , wtC, wtD] : % ;

derivo degli spostamenti di riferimento per i casi i.p. e o.o.p. 
                        | P,H
\\|                     V
\\|====================== 
\\|
  |--------- a ---------|

refipd : P/3/E/J*a^3;

refipd : diff(
    integrate( (P*s)^2/2/E/J , s, 0, a),
    P
);

refoopd : diff(
    integrate( (H*s)^2/2/E/J , s, 0, a),
    H
);

normalized utB

nutB : utB/refipd;

nutC : utC/refipd;

controlli di base

limit( nutB , p , 0 );

limit( nutB/nutC , p , 1 );

risposta i.p. del portale in termini di cedevolezza al punto C

wxplot2d(
    [ nutC , ev(nutB,p=1/2) ],
    [ r , 0.0001 , 10000 ],
    logx,
    [ legend , "norm. utC" , "norm. utB, p=0.5" ]
);

plotto rigidezza relativa
[spostamento trave di riferimento]/[spostamento portale]

wxplot2d(
    [ 1/nutC  ],
    [ r , 0.0001 , 10000 ],
    logx,
    [ legend , "norm. utC" ]
);

1/nutC,r=1;
%,numer;

analizzo il rapporto tra le rigidezze di portale e solo montante
nel caso di carico fuori piano; monitoro lo spostamento in C

normalized wtC e wtD

nwtC : wtC/refoopd;