Teoria 12

Da CdM_unimore.
Versione del 14 lug 2014 alle 10:55 di 190020 (Discussione | contributi) (Sistemi di equazioni nonlineari: metodo di Newton-Raphson)

Gli argomenti trattati in questa lezione riguardano:

  • Sistemi di coordinate locali.
  • Algoritmo Newton Raphson N-D.
  • Fenomeno dello shear-locking iso4, su elemento a geometria semplificata rettangolare (vedi Cook).
  • Analisi dell'errore.


Sistemi di Coordinate Locali

Diapositiva5.JPG


Considerando una generica struttura, ad ogni nodo corrispondono spostamenti in x(u) e y(v) secondo il classico sistema di riferimento Cartesiano x-y. E' opportuno talvolta "distaccarsi" da tale forma di rappresentazione delle coordinate in quanto risulterebbe scomodo trattare casi particolari come quello, ad esempio, del seguente caso di carrello inclinato:


Diapositiva2.JPG


Il problema, in tale situazione, non è dato dall'applicazione dei carichi, dato che le forze possono comunque essere scomposte secondo le diverse componenti, quanto piuttosto dall'impossibilità di imporre che lo spostamento sia nullo in direzione obliqua. Risulta quindi necessario considerare un sistema di riferimento locale s-r ed imporre che sul nodo lo spostamento valga 0, (s=0).

In sostanza, considerato un generico nodo, è indifferente descriverne il moto in uno o nell'altro sistema di riferimento, può essere altresì comodo definire un sistema di riferimento locale. Un generico nodo J può utilizzare un sistema di riferimento diverso da quello di un generico nodo I e da quello globale, su cui viene definita la matrice di rigidezza della struttura.

Impostiamo, di seguito un metodo che consenta di imporre spostamento "personalizzato" nodo per nodo. Consideriamo un generico nodo, libero di avere un vettore spostamento incognito δ:

dove sono i versori dell'asse x e dell'asse y rispettivamente.

Diapositiva3.JPG

Se definisco un secondo sistema di assi, con versori e , ruotato rispetto al primo, lo stesso vettore δ deve essere riscritto come :

Il problema consiste nella trasformazione delle coordinate di un generico nodo, da coordinate globali a coordinate locali. Rappresentiamo e , ovvero i versori del sistema di riferimento locale nel sistema di riferimento gloale xy:

dove, , , e sono i coseni direttori dei versori nel sistema di coordinate x,y.

Sostituendo queste definizioni nell'espressione di δ, si ottiene:

Quest'ultima è un'equazione vettoriale in due componenti, quindi un sistema di due equazioni scalari. Perchè sia soddisfatta l'uguaglianza con dovrà risultare:

in forma matriciale:


Diapositiva4.JPG


Tale matrice indica la trasformazione delle coordinate in un piano. Si noti che apparirebbe più corretto aggiungere l'indice "i", (riferito al nodo i-esimo su cui si sta operando la trasformazione), all'espressione dei coseni direttori contenuti nella matrice, poichè possono essere diversi nodo per nodo. A questo punto è opportuno "stendere" una versione globale valida per tutti i nodi della struttura. Supponendo un cambiamento del sistema di riferimento solo per il nodo i-esimo (per tutti gli altri nodi non è necessaria una traformazione)si ottiene:

MatriceT.png

dove n è il numero di nodi della struttura. Quindi, scrivendo il tutto in forma compatta, l'espressione del vettore di spostamento δ secondo il sistema di riferimento globale, si può scrivere come:

δ = T * δ*

dove T è la matrice di trasformazione T e δ* è il vettore spostamento secondo il sistema di riferimento locale.

Quindi, ricordando che l'espressione "pre-vincolamento" era:

K * δ = F

sostituendo il δ ottenuto precedentemente si ottiene:

K *  T * δ* = F

Poichè la matrice K * T non è simmetrica si deve premoltiplicare ambo i membri per TT , ottenendo:

TT * K *  T * δ* = TT * F

da cui:

K* * δ* = F*

dove K* è simmetrica. Si noti che le normali procedure di vincolamento vengono applicate solo dopo che il sistema sia stato riscritto secondo l'ultima espressione.

Giunti a questo punto si introduce una forma di cambio di sistema di coordinate più complessa. Consideriamo un nodo j del sistema il cui spostamento δj in una qualsiasi direzione non sia nullo. Supponiamo inoltre di voler esprimere tale grado di libertà come combinazione lineare di altri gradi di libertà della struttura (es.: il nodo i_esimo deve muoversi in modo uguale e opposto al nodo j_esimo). In pratica si vuole ricondurre lo spostamento di un nodo in una data direzione ad una funzione. In formule:



In questo modo si potrebbe ad esempio semplificare il problema del tubo in pressione, imponendo ad un dato nodo di spostarsi della stessa quantità, ma in direzione opposta, al nodo ad esso diametralmente opposto, eliminando l'eventuale presenza di moti rigidi.


Immagine.png
u1= 1* u2

Alternativamente, è opportuno considerare una variante dello spostamento δj. Secondo tale variante lo spostamento viene visto come la combinazione lineare degli spostamenti degli altri nodi della struttura sommata ad un termine aggiuntivo Δδj. il cui significato è associabile ad una deviazione dalla legge lineare dello spostamento. In formule:


Supponiamo di avere il vettore δ degli spostamenti in x e y dei gradi di libertà della struttura. Allora:

Matrice2.png

Ottengo una ripartizione della matrice in 9 blocchi, in cui le I sono 'sottomatrici quadrate identità' e gli 0 sono 'sottomatrici quadrate nulle'. La riga centrale rappresenta la definizione dello spostamento del nodo j-esimo come combinazione lineare degli spostamenti degli altri nodi e del termine di deviazione. Nel caso ci fossero altri nodi il cui spostamento fosse così espresso, si aggiungerebbero altre righe nella matrice, come quella in figura. Chiamiamo tale matrice 'L'. Riscrivendo il tutto in forma compatta:

δ = L * δ*

Svolgendo i conti:

LT * K * L * δ* = LT * F

che equivale a:

K* * δ* = F*

Questa forma, con deviazione degli spostamenti, può essere utile ad esempio per l'impostazione dell'analisi di un accoppiamento albero-mozzo, in cui si voglia imporre che i nodi corrispondenti su albero e mozzo non possano nè staccarsi nè compenetrarsi. A questo proposito consideriamo un modello di accoppiamento tra albero e mozzo avente stesso numero di nodi per ciascuno degli elementi e un sistema di coordinate radiale e tangenziale.

Immagine2.png

La condizione matematica per esprimere il vincolo per cui i nodi sull'albero non possono staccarsi 'radialmente' dai nodi sul mozzo è:

ra,i = um,i

Se si volesse impedire anche lo scorrimento relativo dei due corpi si dovrebbe imporre che le componenti tangenziali dello spostamento sono uguali:

ta,i= tm,i

E' importante che i nodi di superficie dell'albero corrispondano ai nodi di superficie del mozzo.

In alternativa, se si vuole considerare l'interferenza si deve imporre che il termine di deviazione sia diverso da zero. Questo oggetto nel Mark Mentat si chiama "servolink" e lega un grado di libertà ad n altri. Di fatto, introducendo il termine di deviazione Δδj, vedo il problema come spostamento relativo tra albero e mozzo, anzichè come spostamento del mozzo rispetto all'albero.

Algoritmo Newton-Raphson per sistemi di equazioni non lineari : Iterazione base

Si consideri un sistema lineare di n equazioni

R(u) = F (u)

nelle n componenti incognite del vettore u, con

  • R : u → Rn , u ∈ C ⊆ Rn
  • F : u → Rn , u ∈ C ⊆ Rn

funzioni vettoriali di variabile vettoriale.

In cui, nel caso specifico della soluzione di sistemi di equazioni derivate dagli equilibri nodali di strutture discretizzate con metodo FEM, si definiscono le seguenti grandezze:

  • u: vettore contenente le componenti di spostamento/rotazione nodale dalla configurazione indeformata (incognite);
  • F(u): vettore contenente le componenti di forza/coppia nodale applicate dall'esterno sul sistema, supposte note per una data configurazione della struttura.
  • R(u): vettore contenente le componenti di azione nodale (uguali e contrarie alle reazioni elastiche della struttura costretta in stato deformato) necessarie a mantenere la struttura in equilibrio nello stato deformativo associato al vettore spostamenti nodali u.

Nel caso particolare di sistema elastico lineare si ha che:

R(u) = K * u, con K matrice di rigidezza.

Si nota che tale interpretazione dei termini dell’equazione R(u) = F(u) è appropriata nel caso le condizioni al contorno siano alle sole forze. Nel caso in cui siano definiti vincoli di spostamento nodale imposto, alcune coppie di termini coniugati Ri(u)-Fi(u) risulteranno modificate in quanto all'equazione di equilibrio i-esima si sostituisce l’identità cinematica tra spostamenti incognito ed imposto.

Una scrittura alternativa prevede la definizione e l’annullamento di un termine di residuo

r(u) = R(u) − F(u) = 0 

Tale scrittura permette di riassumere in un unico termine le variazioni in u di forze e reazioni elastiche, per cui risulta vantaggioso procedere con tale notazione.

Il metodo Newton Raphson è costruito a partire dallo sviluppo in serie di Taylor al primo ordine dell’equazione 'r(u) = 0' di un punto di iterato i-esimo u^i, ossia:

r (u∗) = r (u^i) + Jr (u^i) · (u^∗ - u^i) + o (u^∗ − u^i) = 0.

Sistemi di equazioni nonlineari: metodo di Newton-Raphson

Si considera il problema monodimensionale di un albero rotante con velocità Ω con una massa con eccentricità δ

Massa eccentrica.png


Tale problema è Non Lineare in quanto per piccoli spostamenti lavora a flessione, mentre, per spostamenti via via crescenti essa tende a deformarsi come riportato nella figura di seguito.

Grandi deformazioni.png

In questo caso l'eccessiva deformazione fa sì che si deformino anche i vincoli, e che quindi la trave lavori prevalentemente a trazione. La sua rigidezza, in tal caso, è molto maggiore nel caso di sollecitazione a flessione. Si parla in questo caso di Stiffening, ovvero di situazioni in cui la rigidezza cresce al crescere della deformazione.

un altro aspetto da considerare è la forza centrifuga, dovuta alla massa eccentrica mantenuta in rotazione, che inizialmente vale:

F = m Ω2 δ

e, dopo la deformazione:

F = m Ω2 (δ + u)

Avendo quindi:

  • F ≠ COST
  • EJ ≠ COST

il caso considerato è palesemente non lineare.

Problema generico

r(u*)=0 con u*=soluzione esatta

Il sistema non lineare sarà del tipo: r1(u*)=0 r2(u*)=0 . . . . rn(u*)=0

Partiamo da una soluzione ui di primo tentativo e svolgo una serie di Taylor di r(u*)=0 nell'intorno di ui r(u*)=r(ui)+Jr(ui)*(u*-ui)+o(u*-ui)=0

Se trascuro o(...) l'espressione non vale più con u* soluzione esatta ma "non so cos'è e lo cancello" cit. Bertocchi. r(ui)+Jr(ui)*(u(i+1)-ui)=0 con Jr definito come nella Wiki.

[Jr(ui)]l,m= [d(rl)/d(um)]u=ui estraggo riga l e colonna m

lo Jacobiano è funzione di r e del punto in cui vado a calcolarlo

Jr(ui)=JR(ui)-JF(ui) con JR(ui)= K e se le forze sono costanti JF(ui)=0

r(ui)+Jr(ui)*(u(i+1)-ui)=0 può essere riscritta come

Jr(ui)*(u(i+1)-ui)=-r(ui) con Jr(ui)= coefficiente, (u(i+1)-ui)= incognita, e -r(ui)= termine noto

In Mathlab la notazione corrispondente sarebbe: A x = b quindi x = A\b

quindi:

(u(i+1)-ui)=Jr(ui)\(-r(ui))

u(i+1)=ui-Jr(ui)\r(ui)

Condizione per la convergenza:

||ui-u(i+1)||<eps (agli spostamenti) ||-r(ui)||< eps (residui)

Per l'analisi dettagliata del metodo vedi dispensa provvisoria 2013, sezione 2

Fenomeno dello Shear Locking

slide di riferimento proiettata in classe:

Slide shear locking.png

Locking: E’ un problema tipico degli elementi basati sulla formulazione del principio dei lavori virtuali (metodo degli spostamenti). Gli elementi affetti da locking presentano risultati poco accurati e una scarsa convergenza. Il locking consiste in una sottostima degli spostamenti (responso strutturale troppo rigido).La struttura blocca (“locks”) se stessa contro le deformazioni Matematicamente l’origine del locking è dovuta al malcondizionamento del sistema di equazioni differenziali (alle derivate parziali); all’interno delle equazioni ci potrebbero essere parametri “molto piccoli” (parametro critico) che portano ad elevati coefficienti nel sistema discretizzato di equazioni. Esempio è lo Shear Locking nel caso dei gusci sottili. Nel caso in cui l’elemento finito è un elemento guscio (shell) o si devono gestire materiali incomprimibili (o quasi) come le gomme nel caso delle strutture o fluidi incomprimibili nel caso fluidodinamico, la formulazione basata sugli spostamenti risulta inadeguata e bisogna allora passare ad una FORMULAZIONE MISTA.


Consideriamo di avere uno stato di deformazione lineare ed un elemento isoparametrico a 4 nodi, rettangolare, avente larghezza 2a e altezza 2b:

Iso4.png


Lo Jacobiano di questo elemento è: Jacobiano.png

Questo elemento può essere modellato per rappresentare un concio di trave sottoposto a puro momento flettente. Il concio assumerà la configurazione riportata in figura:


Concio flessione.png


Secondo la teoria dell'elasticità per un concio di trave sottoposto a tale stato tensionale la soluzione esatta è:

Formule 1.jpg

Poichè lo stato di flessione è ottenuto per puro momento, e non applicazione di forza normale, le τ taglianti sono nulle e, quindi, anche ϒ_xy.


Applicando invece tale stato tensionale all'elemento isoparametrico 4 nodi (e non al tratto di trave) si ottiene, per la conservazione della funzione di forma, la seguente configurazione deformata:


Iso flessione.png


In particolare, in formule, si ha:

Formula2.jpg

che rappresenta lo stato deformativo ottenuto sull'elemento, avendo imposto la seguente funzione di spostamento:

Formula3.jpg

in cui:

  • u: spostamento in direzione orizzontale (asse x)del nodo considerato
  • v: spostamento in direzione verticale (asse y)del nodo considerato


Come si può notare dalle formulazioni sopra espresse, nell' elemento isoparametrico a 4 nodi si genera un taglio che NON è previsto dalla teoria. Questa τ deriva dal fatto che, per ottenere l'annullamento delle concavità, in modo da rispettare le funzioni di forma, compare il termine di taglio spurio da cui consegue una deformazione parabolica sull'elemento. In particolare, osservando l'elemento deformato, si nota che agli estremi dello stesso si ha una deformazione a parallelepipedo causata dall'azione tagliante stessa, cioè gli angoli non si mantengono retti.

La deformazione tagliante spuria ϒ_xy è, quindi, causata dall'azione tagliante lineare lungo il lato dell'elemento. Il problema fondamentale è che queste ϒ, associate alle τ, assorbono un quantitativo di energia. Questo comporta che si venga a creare una sostanziale discrepanza tra la previsione della deformazione dell' elemento data da teoria classica e teoria degli elementi finiti. In particolare, per deformare l' elemento isoparametrico, ruotandone ad esempio due lati di un angolo pari a ϴ/2, è necessaria più energia di quella prevista dalla teoria classica, sottostimando la freccia. Si dice quindi che nasce una forza parassita di deformazione.

Si analizza, quindi, di seguito l'errore che può derivare da questo fenomeno.

Analisi dell'errore

Si consideri il caso precedente. A causa del fenomeno dello "Shear Locking" è necessario spendere più energia per avere sull'elemento isoparametrico una deformazione pari a quella prevista dalla teoria classica. A parità di deformazione il rapporto tra la coppia prevista dal modello agli elementi finiti e dalla teoria classica è maggiore di uno e pari a:

Errore.png

in cui:

  • C_el: coppia prevista dalla teoria degli elementi finiti
  • C_b: coppia prevista dalla teoria classica
  • a: lunghezza dell'elemento
  • b: altezza dell'elemento


Anche da questa relazione si evince che, spesso, si sottostima la freccia.


Si noti che il rapporto C_el/C_b varia in funzione della geometria dell'elemento che scelgo per modellare lo spessore della trave:

  • Scegliendo a/b grandi, quindi elementi molto "tozzi", eccito il fenomeno dello "Shear Locking"
  • Scegliendo a/b piccoli, quindi elementi molto "snelli", minimizzo l'errore introdotto da questo fenomeno

Inoltre, variando il rapporto a/b dell'elemento scelto varierà, di conseguenza, il numero di elementi necessari per modellare lo spessore della trave. Quindi, l'errore è indirettamente dipendente dal numero di elementi che utilizzo per modellare la stessa geometria.


N.B. Bisogna fare attenzione alla direzione rispetto a cui definisco l'elemento. Se definisco un elemento avente a/b piccolo ma sbaglio la direzione rispetto a cui lo definisco, eccito il fenomeno dello Shear Locking.

Tornando al caso in analisi di trave sottoposta a momento flettente, si possono fare delle considerazioni derivanti dal tipo e dal numero di elementi usati per la modellazione dello spessore della trave stessa. Ad esempio, gli elementi triangolari a 3 nodi non sono adatti alla modellazione del problema poichè, per avere una soluzione con errori accettabili, sono necessari circa 8 strati di elementi triangolari nello spessore. Al contrario, gli elementi triangolari a 6 nodi ( aventi nodo a centro lato) risultano particolarmente indicati per la modellazione del problema dato che, con un solo elemento nello spessore, si arriva ad una soluzione avente errore accettabile.

Bibliografia

  • Appunti per il corso di Progettazione Assistita, Enrico Bertocchi, 2013.
  • Introduzione ai codici di calcolo agli elementi finiti, - Strutture Aerospaziali -, Sapienza Università di Roma, Luca Lampani.
  • Costruzione di Macchine 2, Lez. 8, Politecnico di Milano, Stefano Beretta.