pulisco la memoria

kill(all);

assume(mm>0 , scala >0);

giochini con assume/forget

sqrt(mm^2);
assume(mm>0);
sqrt(mm^2);
forget(mm>0);
sqrt(mm^2);

inserisco un dimensionamento di riferimento, per controllare
numericamente i risultati

dim : [
    dm= 39 *mm,
    t = 2 *mm,
    E = 210000 *N/mm^2,
    nu = 0.3,
    eta = 2
];
dim : append( 
    dim , 
    ev( 
        [ 
            a=4/3*(dm+t/2)*scala,
            b=    (dm+t/2)*scala
        ]  ,  
        dim 
    )    
);
dim : append( 
    dim , 
    ev( 
        [ 
            l = sqrt(a^2+b^2),
            A = %pi*dm*t,
            J = integrate( dm^3/8 * sin(phi)^2 * t , phi,0,2*%pi),
            G = E / 2 / (1+nu)

        ]  ,  
        dim 
    )    
);

esempio uso append
append( [pippo, pluto] , [ 1 ,2 ,3 , paperino ]  );

inserisco i diagrammi di sollecitazione
in funzione di P, CO, CA, CB

definisco delle funzioni ausiliarie lineari

una prima funzione vale 1 sul primo estremo del tratto di
trave, e 0 sul secondo

variabile indipendente è un'ascissa curvilinea
adimensionale xi che vale 0 al primo estremo e 1 al secondo

N1 : 1 - xi ;

la seconda funzione vale 0 sul primo estremo 
e 1 sul secondo estremo. Tali funzioni
hanno andamento lineare

N2 : xi ;

ora posso definire una quantità linearmente variabile
da un f1 a un f2 in xi da 0 a 1 come

f1*N1 + f2*N2;

vado a definire le caratteristiche di sollecitazione

contributi
            |---P---|-----CB-----|-----CA-----|-----CO-----|

N_OA : +P      +CB/b        +0           -CO/b        $
N_AB : -P*l/a  -CB*a/b/l    -CA*b/a/l    +CO*l/a/b    $
N_BO : +P*b/a  +0           +CA/a        -CO/a        $

Mf_OA: +0      +0           + 0*N1+CA*N2 +CO*N1+ 0*N2 $
Mf_AB: +0      + 0*N1+CB*N2 +CA*N1+ 0*N2 +0           $
Mf_BO: +0      +CB*N1+0 *N2 +0           + 0*N1+CO*N2 $

T_OA : +0      +0           -CA/a        +CO/a        $
T_AB : +0      -CB/l        +CA/l        +0           $
T_BO : +0      +CB/b        +0           -CO/b        $

manca taglio

definisco l'energia potenziale elastica sulla struttura
in funzione di P, CO, CB, CA

U:
    integrate(
        (
            +  N_OA^2/2/E/A
            + Mf_OA^2/2/E/J
            +  T_OA^2/2/G/A * eta
        ) *a 
        ,xi, 0, 1
    )
+
    integrate(
        (
            +  N_AB^2/2/E/A
            + Mf_AB^2/2/E/J
            +  T_AB^2/2/G/A * eta
        ) *l 
        ,xi, 0, 1
    )
+
    integrate(
        (
            +  N_BO^2/2/E/A
            + Mf_BO^2/2/E/J
            +  T_BO^2/2/G/A * eta
        ) *b 
        ,xi, 0, 1
    )
$

U : fullratsimp(U)$

devo imporre che gli angoli tra le barre non varino, in
coerenza con il vincolo originario di incastro interno

posso usare castigliano, considerando che

diff(U, CB) = angolo di apertura al giunto in B, come variazione
              rispetto all'indeformata

dthetaB : diff(U,CB)$
dthetaO : diff(U,CO)$
dthetaA : diff(U,CA)$

impongo che tali variazioni dell'angolo di incidenza tra le
varie barre siano nulle, in coerenza con i giunti saldati

eqns : [
 dthetaB=0,
 dthetaO=0,
 dthetaA=0
]$

risolvo tali equazioni nelle incognite CB,CA,CO
e trovo quella definizione delle stesse nel caso iperstatico

def_iper : linsolve(eqns , [CB,CA,CO]) $

provo a semplificare con fullratsimp

def_iper : fullratsimp (def_iper)$

provo valutazione numerica, per controllo

ev(def_iper, dim, numer , infeval);

per confronto, definisco CB,CA, CO nel caso reso isostatico
tramite le cerniere

def_isos : [ CA = 0, CB=0, CO=0 ];

sostituisco tali definizioni nell'energia potenziale elastica;
U risulterà funzione solamente di P, e avrò il cedimento della struttura
- al punto di applicazione di P
- in direzione e verso di P
con Castigliano

U_iper : ev(U, def_iper ,infeval) ;

U_isos : ev(U, def_isos ,infeval) ;

posso quindi calcolare la freccia in corrispondenza di P 
nei due casi

delta_iper : diff(U_iper,P);

delta_isos : diff(U_isos,P);

confronto i due risultati mediante un rapporto

delta_ratio : delta_isos / delta_iper $

controllo numericamente, in un caso specifico

ev(delta_ratio , scala=10 ,dim, numer , infeval );

grafico la sovrastima del cedimento se calcolato
in forma semplificata isostatica

definisco la grandezza da plottare, ossia l'errore in
sovrastima del cedimento

plottami : ev(delta_ratio - 1  ,dim, numer , infeval )$

dovrebbe ridursi a valore numerico, una volta
assegnato un valore specifico per la variabile indipendente scala,
ascissa del grafico.

ev ( plottami , scala=10 , numer );

plot2d (
    plottami,
    [scala, 1, 100],
    [logx],[logy]
);