kill(all);

contributi ai momenti flettenti e torcenti sui vari tratti
dovuti a F, VA, PsiA, ThetaA;
le altre reazioni vincolari in A sono nulle in quanto
reazioni simmetriche su caso di carico antisimmetrico.

Mf_AB : ( 0)*F + (    -s)*VA + ( 0)*PsiA + (-1)*ThetaA$
Mt_AB : ( 0)*F + (     0)*VA + (-1)*PsiA + ( 0)*ThetaA$
Mf_BC : (+s)*F + (    -s)*VA + ( 1)*PsiA + ( 0)*ThetaA$
Mt_BC : ( 0)*F + (    -b)*VA + ( 0)*PsiA + (-1)*ThetaA$
Mf_CD : (+s)*F + (+(b-s))*VA + ( 0)*PsiA + (+1)*ThetaA$
Mt_CD : (+a)*F + (    -a)*VA + (+1)*PsiA + ( 0)*ThetaA$

le caratteristiche di sollecitazione N, T1 e T2 sono omesse
in quanto ipotizzo trascurabili i relativi contributi energetici.

calcolo l'energia potenziale elastica sulla struttura

U : 
  integrate( Mf_AB^2/2/E/J + Mt_AB^2/2/G/Kt , s, 0, b)
+ integrate( Mf_BC^2/2/E/J + Mt_BC^2/2/G/Kt , s, 0, a)
+ integrate( Mf_CD^2/2/E/J + Mt_CD^2/2/G/Kt , s, 0, b);

ho ottenuto l'energia potenziale
 U(F, VA, ThetaA, PsiA )
in funzione del carico applicato e delle reazioni
vincolari iperstatiche;

utilizzo Castigliano per derivare spostamenti e rotazioni
del punto A; essendo il punto A inizialmente incastrato,
tali cedimenti devono essere imposti nulli.

eq_comp :[
    diff( U , VA     ) = 0,
    diff( U , ThetaA ) = 0,
    diff( U , PsiA   ) = 0 
];

risolvo il sistema di equazioni nelle incognite iperstatiche
e assegno alle stesse il loro valore soluzione

linsolve( eq_comp , [ VA, ThetaA, PsiA ] ) , globalsolve=true;

procedo calcolando il cedimento del punto di applicazione
della forza esterna F nella direzione e nel verso della forza

devo ottenere una forma dell'energia potenziale elastica
nella quale siano esplicitati TUTTI i contributi associati ad F

U(F, VA(F), ThetaA(F), PsiA(F)) = U(F)

U : ev( U ) $

procedo con il calcolo della freccia

delta_F : fullratsimp( diff ( U , F ) );

inserisco un dimensionamento di riferimento
per avere un dato numerico

dim:	[
            a=800,
            b=1000, 
            E=210000,
            G=210000/2/(1+3/10),
            J=(40^4-36^4)*%pi/64,
            Kt=(40^4-36^4)*%pi/32
        ];

ev( [ VA, ThetaA, PsiA, delta_F ] , dim);

valutazione numerica

ev(%,numer);