PK FJBH mimetypetext/x-wxmathmlPK FJ)^9| 9| content.xml
Instabilità in strutture elastiche (sistema nonlineare conservativo)Implementazione elemento finito "barra deformabile elasticamente incernierata a terra"in termini di energia potenziale elasticaelemento 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 rigidezzeassume(a>0,b>0,k>0,k1>0,k2>0);definisco "non enormi" gli spostamentiassume(u<a,v<a,u<b,v<b,-u<a,-v<a,-u<b,-v<b);assume(l0>0);funzioni ausiliarie utilifunzione lunghezzalen(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) Ot(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));
controllot(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 esempioCaso di esempioesempio: Y A O --------- | | | O--> X | k | / | | a / |- | / | \ s | /____O___|____|_ ____O__________/////////////// ///////////////A y||+---> xGradi di libertà:u : spostamento dir. x del punto di applicazione del carico v : spostamento dir. y del punto di applicazione del caricodefinisco una configurazione di riferimento della struttura, relativa ad una pura rotazione di %pi/6ueq : a*sin(%pi/6)$%,numer;veq : a*(cos(%pi/6)-1)$%,numer;Analisi dell'intorno dell'indeformata, carichi costantiricavo componenti di azione esterna tali da equilibrare tale configurazioneUTOT : 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 equilibrioeqns:[diff(UTOT,u)=0,diff(UTOT,v)=0];provo a linearizzarle nell'intorno di una configurazione scarica/indeformatataylor(eqns,[u,v],[0,0],1);estraggo matrice di sistema e termine noto (cambiato di segno) della forma linearizzataaugcoefmatrix(%,[u,v]);le prime due colonne rappresentano la matrice di rigidezza in piccoli spostamenti, la terza il termine noto cambiato di segnoprovo ad estrarre l'hessiana della funzione potenziale totale UTOTla valuto nell'intorno della configurazione scaricahessian(UTOT,[u,v])$ev(%,u=0,v=0);trovo la stessa matrice di rigidezza.Analisi dell'intorno dell'indeformata, carichi variabili CONSERVATIVINOTA: i carichi non sono funzione della configurazione, da cui un contributo nullo all'hessiana.se fosse invece ad esempio X,Y da forza centrifuga, avreiX(u,v):=u*m*omega^2;Y(u,v):=(a+v)*m*omega^2;percorso di integrazione del potenziale da posizione deformata aconfigurazione di riferimento indeformata (pot. U0) su retta congiungenteUTOTbis : 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 riferimentoricerca delle forze equilibranti la condizione di riferimentopotenzialeUTOT : U(0,0,0,a,u,v,k,s) + X*(u0-u)+Y*(v0-v)+U0;equazioni di equilibrio, genericheeqns:[diff(UTOT,u)=0,diff(UTOT,v)=0];equazioni di equilibrio, alla configurazione di riferimentofullratsimp(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 equilibrioev(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 pienoAnalisi variematrice hessiana dell'energia potenziale elasticahessian(U,[u,v])$concide con la matrice di rigidezza tangenteforma nell'intorno di [u,v]=[0,0]ev(%,u=0,v=0);definisco matrice di rigidezza elementodefine(K(x0,y0,x1,y1,k,s),%);verifica su esempioK(0,0,0,a,k,s);suppongo ora l'elemento precaricato nella configurazione indeformata da uno sforzo normale N, positivo se trattivovado 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 formahessian(L,[u,v]);valuto nell'intorno di spostamenti nulliev(%,u=0,v=0);definisco la matrice di rigidezza geometrica = matrice di rigidezza da precaricodefine(G(x0,y0,x1,y1,N),%);FORMA ALTERNATIVA:in alternativa (e per controllo) posso definire tale lavoro come la differenza traa) l'energia potenziale elastica CON precarico iniziale eb) l'energia potenziale elastica SENZA precarico inizialeLbis : 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 nulliev(%,u=0,v=0);FORMA ALTERNATIVA:definisco la matrice di rigidezza geometrica = matrice di rigidezza da precaricodefine(Gbis(x0,y0,x1,y1,N),%);FORMA ALTERNATIVA:le due forme sono EQUIVALENTIfullratsimp(Gbis(x0,y0,x1,y1,N)-G(x0,y0,x1,y1,N));esempio: | P | V O --------- | | | k | | | a |- | | \ s |____O___|____|_/////////////// A y||+---> xGradi di libertà:u : spostamento dir. x del punto di applicazione del carico v : spostamento dir. y del punto di applicazione del caricomatrice di rigidezza geometrica associata a precarico compressivo unitarioG(0,0,0,a,-1);matrice di rigidezza elasticaK(0,0,0,a, k,s);Linearized Pre Buckling analysisK(0,0,0,a, k,s)+P*G(0,0,0,a,-1);determinant(%);solve(%,P);in forma di estrazione autovalori/autovettoririduco a problema agli autovalori standard utilizzando(K + P*G).v =0(inv(K).K + P*inv(K).G).v = inv(K).0con K invertibile se correttamente vincolato( I + P*inv(K).G).v = 0( inv(K).G + (1/P) I ).v =0formalmente analogo a ( A - lambda I ).v =0con lambda = -1/P autovalorev 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/Pvals[1];mults[1];carico critico-1/vals[1];autovettori associati, potenzialmente multiplivecs[1];modo instabile con puro spostamento u, mentre v=0il 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||+---> xGradi di libertà:u : spostamento dir. x del punto di applicazione del carico v : spostamento dir. y del punto di applicazione del caricomatrice di rigidezza elastica data dall'assemblaggio dei due elementiKel : K(0,0,0,b,k2,0) + K(0,b,a,b,k1,0);matrice di rigidezza geometrica associata a precarico compressivo unitarioKG : G(0,0,0,b,0) + P * G(0,b,a,b,-1);Linearized Pre Buckling analysisdeterminant(Kel + KG);solve(%,P);Caso | P | V O -------------- k + + k | + + | b____O____ ____O____ ____|_///////// ///////// | | | | a | a | |-----|-----| A y||+---> xGradi di libertà:u : spostamento dir. x del punto di applicazione del carico v = -delta : spostamento dir. y del punto di applicazione del caricodimensionamento : [b=1,a=2,k=1];matrice di rigidezza elastica data dall'assemblaggio dei due elementiKel : 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 unitarioKG : ( 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 analysisdeterminant(Kel + KG);fullratsimp(%);solve(%,P);es. dimensionamento specifico, vedi sotto.Pcrit_LPBA : ev(%,dimensionamento,numer);rimuovo un'eventuale precedente assegnazione di Pkill(P);analisi nonlinearesomma contributi- energia potenziale elastica elemento (-a,0)-(0,b)- energia potenziale elastica elemento (+a,0)-(0,b)- energia potenziale carico FUtot: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 potenzialeequilibrio:[diff(Utot,u)=0,diff(Utot,v)=0];definisco P come equilibrante del caso u=0ev(equilibrio,u=0);linsolve(%[2],P),globalsolve=True;plotto P(delta) = P(-v) per un dimensionamento specificoplotme: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=0hessian(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 secantedefinisco un Hessiano "ridotto" 1x1, con u=0 e unica variabile vHred: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 secantecondizione 0; delta = -v = 0.3v0:-0.3;P0:ev(P,dimensionamento,v=v0,numer);Hred0 : ev(Hred,v=v0,dimensionamento,numer);condizione 1; delta = -v = 0.4v1:-0.4;P1:ev(P,dimensionamento,v=v1,numer);Hred1 : ev(Hred,v=v1,dimensionamento,numer);ricerco condizione di Hessiano ridotto singolarelinsolve( Hred0 + zeta * (Hred1-Hred0), zeta),globalsolve=True,numer;carico a cui ottengo tale condizioneP0 + zeta*(P1-P0);non lontanissimo dal massimo ottenuto su grafico se i punti di campionamento sono prossimi.
PK FJgn=> =>
image1.pngPNG
IHDR X W sBITO IDATx`Ms#3X36JgVE2Bu}T+a([
ǏL/Q&vOk{y>r;yϙ=uHj@DD$!i!i!i!i!i!i!i!i!iD`4[WfO>DG6FDDvA
a^^óǏĉRmƍgϞ.]IHDDA
aBBSXXN
.fԨQkB}#:u{yy-Rm
Sn݂+#"5rBh6]\\fs5SPP :ԺKt:U Fz}nnnbNN/Qe8H!NMM-ZLNNd weɒ%%!!!...$$@LLLIDDT\?r"""
k>ĉ4h ??VZk3f̢EJXwqZIKf ]vݷo_9rkvCz-
wuqdu[mҝ:RzXHXHXHXHXHXHXHXHXHXHXHXHXHXHXHgzRIVt3.\+/\
NNh͛
h|}ѢDD#