pulisco la memoria

kill(  all  )  ;

inserisco le equazioni di equilibrio 
del quarto di supporto

equilibrio alle traslazioni rispetto a z,y,x

eqtz : WA + WD = 0 ;

eqtx : UA = 0;

eqty : VD = 0;

equilibrio alle rotazioni rispetto agli assi 
- passanti per O
- orientati in x,y,z

eqrOx : + C/2      + UA*0        + WA*0
        + ThetaD          + VD*0 + WD*(b+e) 
        = 0;

eqrOy : PhiA - a * WA ;

eqrOz : 0 = 0;

definisco il sistema di equazioni

eqns : [ eqtx , eqty , eqtz , eqrOx , eqrOy , eqrOz ]  ;

unks : [ UA , WA , PhiA , VD, WD, ThetaD ] ;

risolvo il sistema di equazioni

linsolve ( eqns , unks ) ;

mi accorgo che il sistema è a rango non pieno
ho solo 5 equazioni indipendenti per le 6 incognite
decido di trattare ThetaD come parametro (e non come incognita)
ridefinendo la lista delle incognite

unks : [ UA , WA , PhiA , VD, WD] ;

adesso ho un sistema di 6 equazioni, di cui 5 indipendenti
in 5 incognite e 1 parametro 
(o 2 parametri se considero anche ThetaA)

sols : linsolve ( eqns, unks );

assegno alle incognite l'espressione delle stesse
che risolve il sistema di equazioni

opzione 1:

WA : ev ( WA , sols );
UA : ev ( UA , sols );
PhiA : ev ( PhiA , sols );
VD : ev ( VD , sols );
WD : ev ( WD , sols );

opzione 2

risolvo le equazioni con opzione aggiuntiva

globalsolve=True

globalsolve è una variabile che il codice del linsolve testa
durante l'esecuzione

se globalsolve=False (default), non assegna alle incognite
il valore soluzione del sistema; l'assegnazione è da eseguirsi
a posteriori come da opzione 1

altrimenti (globalsolve=True) assegno automaticamente alle 
incognite il valore soluzione per il sistema.

linsolve ( eqns, unks )  , globalsolve=True ;

controllo che l'assegnazione sia avvenuta con successo.

WA;

WD;

ok, pare avvenuta con successo.

Inserisco ora le caratteristiche di sollecitazione

Sforzo normale e tagli

N : 0;

T_eta : -WD ;

T_xi  : -VD ;

Momenti flettenti e torcente, in funzione di un'ascissa
"curvilinea" dimensionale che scorre lungo il tratto di trave
0 <= s <= a

Mt : -ThetaD;

Mf_xi : -WD*s;

Mf_eta : +VD*s;

Una volta ottenute le espressioni delle caratteristiche di
sollecitazione su tutti i tratti di trave, procedo alla
valutazione dell'energia potenziale elastica sulla struttura.

Tale energia è la sommatoria di tutti i contributi energetici
associati alle varie porzioni di struttura.

La sezione in oggetto ha assi xi, eta principali d'inerzia in quanto
simmetrica; il materiale è inoltre supposto omogeneo e isotropo

densità lineica di en. pot. elastica

dU_ds :             Mf_xi ^2 / 2 / E / J_xixi
        +           Mf_eta^2 / 2 / E / J_etaeta
        +           Mt    ^2 / 2 / G / Kt
        +           N     ^2 / 2 / E / area
        + chi_xi  * T_xi  ^2 / 2 / G / area
        + chi_eta * T_eta ^2 / 2 / G / area ;

integro la densità lineica di en. pot. el. e ottengo l'en. pot. 
elastica della struttura.

U : integrate ( dU_ds , s , 0, a );

nota l'espressione per l'energia pot. el. della struttura,
applico il teoreoma di Castigliano per ricavare la 
componente x della rotazione della sezione in D

esempio

d^2(f(a*t)*sin(a*t + b))
--------------------
      da db

diff( f(a*t)*sin(a*t + b), a, 1 , b, 1);

thetaD : diff(  U , ThetaD );

ricavo quel valore di ThetaD che annulla la rotazione thetaD
ripristinando la continuità del materiale al piano di 
antisimmetria yz tra la porzione di struttura modellata
e la sua immagine

eq_thetaD : thetaD = 0;

prima di usare linsolve, controllo che l'equazione
sia effettivamente lineare

diff(eq_thetaD , ThetaD, 2);

ok, è effettivamente un'equazione lineare.

linsolve( eq_thetaD , ThetaD ) , globalsolve=True;

forzo la sostituzione dei termini ThetaD presenti entro U
secondo la sua più recente definizione con la forma 

ThetaD = ThetaD(C)

appena trovata

U :
    fullratsimp(
        ev( U ) 
    );

l'energia pot. elastica appena trovata è relativa al
quarto di struttura.

L'energia pot. el. per la struttua completa è 4 volte quella del
quarto.

Ustruttcomp : 4 * U;

calcolo ora la rigidezza a torsione della struttura di supporto
completa, ottendo 

A) la rotazione relativa tra i due corpi rigidi se suppongo 
l'oggetto non vincolato nello spazio

B) se immagino uno dei corpi rigidi incastrato, la rotazione
assoluta dell'altro

theta_strutt : diff( Ustruttcomp , C );

definisco ora la rigidezza torsionale dell'oggetto come rapporto
tra coppia torsionale applicata e rotaz. torsionale indotta

Kt_strutt : fullratsimp( 
                C/theta_strutt
            );

non è strettamente necessario che i profili presi singolarmente
abbiano rigidezza torsionale elevata per garantire
una rigidezza torsionale elevata alla struttura, ho infatti

limit ( Kt_strutt, Kt , 0 );