# Instabilità in strutture elastiche (sistema nonlineare conservativo)

## Implementazione elemento finito "barra deformabile elasticamente incernierata a terra"

in termini di energia potenziale elastica

elemento utile alla presente dissertazione in quanto
- già vincolato (non occorre impostare una procedura di vincolamento per ottenere un esempio)
- solo 2 gdl (semplice assemblaggio)

kill(all);

definisco >0 proprietà geometriche e rigidezze

assume(a>0,b>0,k>0,k1>0,k2>0);

definisco "non enormi" gli spostamenti

assume(u<a,v<a,u<b,v<b,-u<a,-v<a,-u<b,-v<b);

assume(l0>0);

funzioni ausiliarie utili

funzione lunghezza

len(x1,y1,x2,y2):=sqrt((x2-x1)^2+(y2-y1)^2);

componenti seno e coseno, definizione nonlineare
                      P (x2,y2)  -------
                    +                 ^
                  +                   | sin
                +                     |
      (x1,y1) O        -----------------
              |       |
              |  cos  |
              |------>|
attenti ai segni!!

s(x1,y1,x2,y2):=(y2-y1)/len(x1,y1,x2,y2)$
c(x1,y1,x2,y2):=(x2-x1)/len(x1,y1,x2,y2)$

funzione rotazione angolare rispetto a indeformato

           (x3,y3) Q 
                  /   P (x2,y2)
                 /t +
                / +
               /+
      (x1,y1) O

t(x1,y1,x2,y2,x3,y3):= atan2(

(x2*y3+x3*y1+x1*y2-x2*y1-x3*y2-x1*y3)
,
(x3-x1)*(x2-x1)+(y3-y1)*(y2-y1)
);

controllo

t(1,2,1+1,2+1,1+0,2+1);

t(1,2,1+1,2+1,1+1,2+0);

t(0,0,0,1,1,-1);

t(0,0,0,1,-1,-1);

energia potenziale elastica dell'elemento barra incernierata a telaio al punto (x0,y0),
connesso al punto mobile (x1,y1) con posizione deformata (x1+u,y1+v)

U : 1/2*k*(len(x0,y0,x1+u,y1+v)-len(x0,y0,x1,y1))^2 
  + 1/2*s*t(x0,y0,x1,y1,x1+u,y1+v)^2;

funzione potenziale, nel caso serva.

define(U(x0,y0,x1,y1,u,v,k,s) , U);

## Analisi della condizione di equilibrio in casi di esempio

### Caso di esempio

esempio:

                            Y A 
    O ---------               |
    |        |                O--> X
    | k      |               /
    |        | a            /
    |-       |             /
    |  \ s   |            /
____O___|____|_      ____O__________
///////////////      ///////////////
A y
|
|
+---> x
Gradi di libertà:
u : spostamento dir. x del punto di applicazione del carico 
v : spostamento dir. y del punto di applicazione del carico

definisco una configurazione di riferimento della struttura, relativa ad una pura rotazione di %pi/6

ueq : a*sin(%pi/6)$
%,numer;
veq : a*(cos(%pi/6)-1)$
%,numer;

### Analisi dell'intorno dell'indeformata, carichi costanti

ricavo componenti di azione esterna tali da equilibrare tale configurazione

UTOT : U(0,0,0,a,u,v,k,s) + X*(u0-u)+Y*(v0-v)+U0;

cerco punto stazionario dell'energia potenziale totale <=> punto di equilibrio

eqns:[diff(UTOT,u)=0,diff(UTOT,v)=0];

provo a linearizzarle nell'intorno di una configurazione scarica/indeformata

taylor(eqns,[u,v],[0,0],1);

estraggo matrice di sistema e termine noto (cambiato di segno) della forma linearizzata

augcoefmatrix(%,[u,v]);

le prime due colonne rappresentano la matrice di rigidezza in piccoli spostamenti, 
la terza il termine noto cambiato di segno

provo ad estrarre l'hessiana della funzione potenziale totale UTOT
la valuto nell'intorno della configurazione scarica

hessian(UTOT,[u,v])$
ev(%,u=0,v=0);

trovo la stessa matrice di rigidezza.

### Analisi dell'intorno dell'indeformata, carichi variabili CONSERVATIVI

NOTA: i carichi non sono funzione della configurazione, da cui un contributo nullo all'hessiana.
se fosse invece ad esempio X,Y da forza centrifuga, avrei

X(u,v):=u*m*omega^2;

Y(u,v):=(a+v)*m*omega^2;

percorso di integrazione del potenziale da posizione deformata a
configurazione di riferimento indeformata (pot. U0) su retta congiungente

UTOTbis : U(0,0,0,a,u,v,k,s) 
+ integrate(X(u*aux,v*aux)*u,aux,1,0)
+ integrate(Y(u*aux,v*aux)*v,aux,1,0)
+ U0
;

eqnsbis:[diff(UTOTbis,u)=0,diff(UTOTbis,v)=0];

taylor(eqnsbis,[u,v],[0,0],1);

augcoefmatrix(%,[u,v]);

hessian(UTOTbis,[u,v])$
ev(%,u=0,v=0);

notate che il carico in questo caso entra nell'hessiana del potenziale, 
aggiungendo contributi alla matrice di rigidezza già nell'intorno di spostamento nullo.

### Analisi nell'intorno della configurazione di riferimento

ricerca delle forze equilibranti la condizione di riferimento

potenziale

UTOT : U(0,0,0,a,u,v,k,s) + X*(u0-u)+Y*(v0-v)+U0;

equazioni di equilibrio, generiche

eqns:[diff(UTOT,u)=0,diff(UTOT,v)=0];

equazioni di equilibrio, alla
 configurazione di riferimento

fullratsimp(ev(eqns,u=ueq,v=veq));

linsolve(%,[X,Y]);

[Xeq,Yeq] : ev([X,Y],%);

ora ho le forze da applicare per mantenere in equilibrio tale configurazione.

considero ora una perturbazione in forza e spostamento nell'intorno di tale configurazione di equilibrio

ev(UTOT,u=ueq+du,v=veq+dv,X=Xeq+dX,Y=Yeq+dY);

eqns_pert:[diff(%,du)=0,diff(%,dv)=0];

taylor(eqns_pert,[du,dv],[0,0],1);

augcoefmatrix(%,[du,dv]);

questo sistema ha soluzione univoca solo se la matrice di rigidezza TANGENTE ha rango pieno

## Analisi varie

matrice hessiana dell'energia potenziale elastica

hessian(U,[u,v])$

concide con la matrice di rigidezza tangente

forma nell'intorno di [u,v]=[0,0]

ev(%,u=0,v=0);

definisco matrice di rigidezza elemento

define(
K(x0,y0,x1,y1,k,s)
,
%
);

verifica su esempio

K(0,0,0,a,k,s);

suppongo ora l'elemento precaricato nella configurazione indeformata da uno sforzo normale N, positivo se trattivo

vado a considerare il lavoro svolto da tale forzi sugli spostamenti (u,v)

L :
  N * (  len(x0,y0,x1+u,y1+v)
        -len(x0,y0,x1  ,y1  )
      )
;

tale lavoro è da considerarsi un lavoro di forze interne, e si addiziona alla variazione di energia potenziale, 
da cui una correzione alla matrice di rigidezza nella forma

hessian(L,[u,v]);

valuto nell'intorno di spostamenti nulli

ev(%,u=0,v=0);

definisco la matrice di rigidezza geometrica = matrice di rigidezza da precarico

define(
G(x0,y0,x1,y1,N)
,
%
);

FORMA ALTERNATIVA:
in alternativa (e per controllo) posso definire tale lavoro come la differenza tra
a) l'energia potenziale elastica CON precarico iniziale e
b) l'energia potenziale elastica SENZA precarico iniziale

Lbis : 
fullratsimp(
    ( 
          1/2*k*(len(x0,y0,x1+u,y1+v)-len(x0,y0,x1,y1) + N/k )^2 
        + 1/2*s*t(x0,y0,x1,y1,x1+u,y1+v)^2
   )
    -
    (
          1/2*k*(len(x0,y0,x1+u,y1+v)-len(x0,y0,x1,y1) + 0/k )^2 
        + 1/2*s*t(x0,y0,x1,y1,x1+u,y1+v)^2
    )
);

hessian(Lbis,[u,v]);

FORMA ALTERNATIVA:
valuto nell'intorno di spostamenti nulli

ev(%,u=0,v=0);

FORMA ALTERNATIVA:
definisco la matrice di rigidezza geometrica = matrice di rigidezza da precarico

define(
Gbis(x0,y0,x1,y1,N)
,
%
);

FORMA ALTERNATIVA:
le due forme sono EQUIVALENTI

fullratsimp(Gbis(x0,y0,x1,y1,N)-G(x0,y0,x1,y1,N));

esempio:

    | P
    |
    V
    O ---------
    |        |
    | k      |
    |        | a
    |-       |
    |  \ s   |
____O___|____|_
/////////////// 


A y
|
|
+---> x

Gradi di libertà:
u : spostamento dir. x del punto di applicazione del carico 
v : spostamento dir. y del punto di applicazione del carico

matrice di rigidezza geometrica associata a precarico compressivo unitario

G(0,0,0,a,-1);

matrice di rigidezza elastica

K(0,0,0,a, k,s);

Linearized Pre Buckling analysis

K(0,0,0,a, k,s)+P*G(0,0,0,a,-1);

determinant(%);

solve(%,P);

in forma di estrazione autovalori/autovettori

riduco a problema agli autovalori standard utilizzando

(K + P*G).v =0

(inv(K).K + P*inv(K).G).v = inv(K).0

con K invertibile se correttamente vincolato

( I + P*inv(K).G).v = 0
( inv(K).G + (1/P) I ).v =0

formalmente analogo a 

( A - lambda I ).v =0

con  

lambda = -1/P autovalore
v autovettore

[[vals,mults],vecs] : eigenvectors(invert(K(0,0,0,a, k,s)).G(0,0,0,a,-1));

l'inverso degli autovalori coincide con l'opposto del carico critico -1/P

vals[1];

mults[1];

carico critico

-1/vals[1];

autovettori associati, potenzialmente multipli

vecs[1];

modo instabile con puro spostamento u, mentre v=0

il secondo autovalore 

vals[2];

è nullo ed è quindi associato a carico critico infinito.

geometria di riferimento:
      |----- a ----|
     n2          n1|/
P ==> o-----k1----o|/ ----    
      |            |/   |
      |                 |
      k2                |
      |                 | b
      |                 |
  n3  o             -----
   -------
   ///////
nodo 2 spostamenti "u" in direzione x+ e "v" in direzione y+

A y
|
|
+---> x

Gradi di libertà:
u : spostamento dir. x del punto di applicazione del carico 
v : spostamento dir. y del punto di applicazione del carico

matrice di rigidezza elastica data dall'assemblaggio dei due elementi

Kel : K(0,0,0,b,k2,0) + K(0,b,a,b,k1,0);

matrice di rigidezza geometrica associata a precarico compressivo unitario

KG : G(0,0,0,b,0) + P * G(0,b,a,b,-1);

Linearized Pre Buckling analysis

determinant(Kel + KG);

solve(%,P);

Caso

          | P
          |
          V
          O   --------------
     k  +   +  k          |
      +       +           | b
____O____   ____O____ ____|_
/////////   /////////
    |     |     |
    |  a  |  a  |
    |-----|-----|        


A y
|
|
+---> x

Gradi di libertà:
u : spostamento dir. x del punto di applicazione del carico 
v = -delta : spostamento dir. y del punto di applicazione del carico

dimensionamento : [b=1,a=2,k=1];

matrice di rigidezza elastica data dall'assemblaggio dei due elementi

Kel : K(-a,0,0,b,k,0) + K(a,0,0,b,k,0);

Kel:fullratsimp(Kel);

matrice di rigidezza geometrica associata a precarico compressivo unitario

KG : ( G(-a,0,0,b,-1) + G(a,0,0,b,-1) ) * (P*sqrt(a^2+b^2)/2/b);

KG :fullratsimp(KG);

Linearized Pre Buckling analysis

determinant(Kel + KG);

fullratsimp(%);

solve(%,P);

es. dimensionamento specifico, vedi sotto.

Pcrit_LPBA : ev(%,dimensionamento,numer);

rimuovo un'eventuale precedente assegnazione di P

kill(P);

analisi nonlineare

somma contributi
- energia potenziale elastica elemento (-a,0)-(0,b)
- energia potenziale elastica elemento (+a,0)-(0,b)
- energia potenziale carico F

Utot:U(-a,0,0,b,u,v,k,0)+U(a,0,0,b,u,v,k,0)+P*v;

definisco condizione di equilibrio come punto stazionario energia potenziale

equilibrio:
[diff(Utot,u)=0,diff(Utot,v)=0];

definisco P come equilibrante del caso u=0

ev(equilibrio,u=0);

linsolve(%[2],P),globalsolve=True;

plotto P(delta) = P(-v) per un dimensionamento specifico

plotme:ev(P,dimensionamento,v=-delta);

wxplot2d(
plotme
,
[delta,0,2]
);

osservo un punto stazionario di P(-v)=P(delta)

deltacrit : find_root(diff(plotme,delta),delta,0.3,0.5);

Pcrit: ev(plotme, delta=deltacrit),numer;

noto che questo valore è notevolmente inferiore a quelli predetti dalla procedura 
Linearized pre-buckling analysis.

Pcrit_LPBA;

calcolo la matrice Hessiana del potenziale nel dimensionamento specifico e in v=-deltacrit,u=0

hessian(Utot,[u,v]);
ev(%,u=0,v=-deltacrit,dimensionamento),numer;

si nota che l'hessiana della funzione potenziale (=matrice di rigidezza tangente) 
risulta singolare nell'intorno della condizione critica P=Pcrit, v=-deltacrit,u=0.

Potrei quindi aspettarmi che il metodo proposto predica correttamente tale condizione.

provo il metodo secante

definisco un Hessiano "ridotto" 1x1, con u=0 e unica variabile v

Hred:
diff(
    U(-a,0,0,b,0,v,k,0)+U(a,0,0,b,0,v,k,0),
    v,
    2
);

definisco le due condizioni di riferimento per la secante

condizione 0; delta = -v = 0.3

v0:-0.3;

P0:ev(P,dimensionamento,v=v0,numer);

Hred0 : ev(Hred,v=v0,dimensionamento,numer);

condizione 1; delta = -v = 0.4

v1:-0.4;

P1:ev(P,dimensionamento,v=v1,numer);

Hred1 : ev(Hred,v=v1,dimensionamento,numer);

ricerco condizione di Hessiano ridotto singolare

linsolve(
    Hred0 + zeta * (Hred1-Hred0),
    zeta
),globalsolve=True,numer;

carico a cui ottengo tale condizione

P0 + zeta*(P1-P0);

non lontanissimo dal massimo ottenuto su grafico se i punti di campionamento sono prossimi.