Teoria 11

Da CdM_unimore.

Teoria Degli Elementi Isoparametrici

Gli elementi finiti isoparametrici, piani o tridimensionali, sono elementi finiti "distorti" che si adattano facilmente allo studio e alla riproduzione di geometrie relativamente complesse, essendo, come vedremo, molto semplici da integrare e derivare. Una spiegazione introduttiva sull'idea degli elementi finiti isoparametrici relativa al concetto di coordinate locali puo' essere questa: un elemento finito viene guardato con una lente che lo rende quadrato con i lati di lunghezza 2(dove la variabile scorre tra -1 ed 1). La teoria risulta semplificata dalla forma quadrata e dalla lunghezza normalizzata dei lati. Una volta scritta la teoria, una "controlente" riporta l'elemento finito alla sua forma originale.

Definizione degli spostamenti nodali

Un elemento finito isoparametrico può esser descritto tramite due sistemi di coordinate:

1) LOCALE: di coordinate ,

2) GLOBALE: di coordinate ,

Si consideri un elemento finito isoparametrico a quattro nodi.

Grafico2 8apr.PNG

Di sopra, è riportata la rappresentazione dell’elemento finito in coordinate globali e : come si può notare, esso assume una forma quadrangolare generica. Di seguito, invece, è riportato l’elemento finito a quattro nodi rispetto alle coordinate locali e : l'elemento isoparametrico, definito sul sistema locale, è un quadrato, i cui lati possiedono lunghezza pari a 2. Si precisa che i quattro nodi dell'elemento sono stati definiti mediante i numeri 1, 2, 3, 4, sia nella numerazione globale, sia in quella locale, coerentemente con la scelta adottata nella dispensa.

Grafico1 8apr.png

Come illustrato in figura, è possibile visualizzare nel piano , una versione trasformata degli assi ,, congiungendo i punti medi dei lati. Il passaggio dalla geometria dell’elemento distorto al quadrato, ovvero dalle coordinate globali alle coordinate locali, avviene tramite funzioni di forma che descrivono la coordinata globale, ad esempio , in termini delle coordinate globali nodali , tramite funzioni delle coordinate locali (,) del tipo . Analogo discorso per la coordinata globale . Per ogni nodo dell’elemento finito si avranno le seguenti funzioni di forma, associate ai quattro nodi dell'elemento:

Si può allora, per comodità, definire il vettore colonna:

La coordinata globale può esser allora espressa come semplice prodotto scalare tra il vettore contenente le funzioni di forma trasposto (in modo che risulti un vettore riga) e il vettore colonna che contiene le coordinate dei nodi. Infatti, si ha:

In pratica, questa scrittura è equivalente alla seguente sommatoria:

che rappresenta la mappatura della coordinata ; vale a dire: preso un punto qualunque sul dominio locale, si prendono le specifiche coordinate e e si trova l'associato punto nel dominio globale xy, calcolando e .

Per trovare la mappatura della coordinata si procede allo stesso modo:

L'elemento è definito isoparametrico perchè le stesse funzioni di forma utilizzate per definire la trasformazione di coordinate da globali a locali vengono usate anche per definire gli spostamenti incogniti. Quindi, ad esempio, lo spostamento in direzione x sarà definito in funzione di coordinate locali (,). In forma matriciale si ha:

Analogamente, si procede per il calcolo degli spostamenti in direzione .

Si sottolinea che si tratta di spostamenti definiti in funzione delle coordinate locali e , quindi sono delle forme bilineari in e

ma non in e , poichè la loro forma è più complessa, in quanto non deriva da una forma bilineare.

Calcolo della matrice di rigidezza

Vediamo ora come si costruisce la matrice di rigidezza di un elemento isoparametrico a quattro nodi. Per gli elementi triangolari si è trovata l'energia interna dell'elemento, si è ricavato il lavoro delle forze esterne sui tre nodi, e si è partiti dall'equivalenza di queste ultime. Si definisce l'energia interna in forma integrale dell'elemento considerato:

con =area dell’elemento. Quest'area è calcolata nel dominio xy, perciò non è pari a 4; e soprattutto, non trattandosi di un quadrilatero regolare, il calcolo dell'area risulta piuttosto complesso. Invece utilizzando il quadrato unitario in coordinate e l'integrazione risulta facilitata.

Negli elementi isoparametrici, a differenza degli elementi triangolari a 3 nodi dove le tensioni e le deformazioni restano costanti all'interno dell'integrale suddetto, risulta difficile calcolare esattamente tale integrale. Infatti poichè discende direttamente da tramite la matrice secondo la relazione:

resta difficoltoso il computo dell'integrale e della determinazione della . La questione integrale è stata storicamente risolta usando un metodo numerico, nello specifico un'integrazione approssimata con metodo di Gauss.

Si consideri una generica funzione (per esempio, polinomiale), di cui si vuole calcolare l'integrale da -1 a 1. Si applichi l'integrazione Gaussiana al dominio monodimensionale, rappresentato in figura.

Fig5.jpg

La funzione è definita sull'asse ; è perciò funzione di una sola variabile indipendente . Andando a cercare la migliore approssimazione dell'integrale per un polinomio, volendo campionare la funzione solo in due punti, tali due punti si trovano in corrispondenza di = e = .

In particolare, si ha:

in cui l'integrale, secondo la teoria di Gauss, è approssimabile mediante la sommatoria della funzione calcolata nei punti di Gauss e moltiplicata per un opportuno peso (nello specifico, unitario).

Facendo tesoro delle conclusioni effettuate per il caso monodimensionale, è possibile estendere tale metodologia a Domini bidimensionali di forma semplice, regolare, quadrilatera, e quindi anche al Dominio dell'elemento isoparametrico a quattro nodi in coordinate locali, in cui si ha un'area di integrazione da -1 a 1 in e da -1 a 1 in .

Si consideri una funzione , definita sul seguente dominio. Si definiscano i quattro punti di Gauss nelle coordinate ed , come illustrato nella figura sottostante:

Fig8.jpg

punto 1 =

punto 2 =

punto 3 =

punto 4 =

Si ottiene quindi l’approssimazione:

Pertanto, nel caso in cui si voglia integrare la funzione , su un dominio quadrato, nelle coordinate ed che vanno da -1 ad 1, è possibile utilizzare un'integrazione a quattro punti di Gauss, le cui coordinate sono prese dall'ottimizzazione fatta per i punti di Gauss nel caso monodimensionale. Il peso è unitario, in quanto, per la simmetria che regola il dominio, ogni nodo possiede lo stesso peso. Se si campiona una funzione unitaria, cioè se vale una costante unitaria, siccome l'integrale deve risultare 4 e la funzione campionata vale 1, il peso con cui queste unità devono entrare è unitario.

Si possono effettuare integrazioni di Gauss a più punti: ad esempio, si può considerare un unico punto posto nel centroide dell’elemento (cioè a coordinate ed nulle) ed associargli un peso pari a 4.

Fig6.jpg

Oppure è possibile considerare, oltre ai punti ai vertici, ulteriori punti appartenenti agli assi , fino ad arrivare a 9 punti di Gauss.

Fig7.jpg

Tutto ciò tenendo presente che, poiché l’area dell’elemento a 4 nodi è pari a 4, la somma dei pesi dovrà uguagliare la superficie dell’elemento.

Per l'elemento isoparametrico a 4 nodi si parlerà di full integration se i punti di integrazione sono 4; si parlerà invece di reduced integration se il numero di punti di integrazione è pari a 1. Per gli elementi a spostamenti quadratici (triangolo 6 nodi o quadrilatero 8 nodi con nodi di centro lato), i punti di integrazione sono di più (il quadrilatero 8 nodi è integrato in 9 punti nella versione full integration, mentre è integrato in soli 4 punti nella versione reduced integration). In genere si tende ad utilizzare, per l'elemento isoparametrico a 4 nodi, l'integrazione alla gauss a 4 punti, ovvero l'integrazione più povera che non dia risultati spuri.

Come detto il problema si riduce al calcolo delle (,), ovvero del tensore/vettore:

Poiché l'integrale associato all'energia viene svolto numericamente usando i punti di Gauss, lo stato di deformazione dell'elemento deve essere campionato solo sui punti di integrazione. Quindi, non sarà necessario avere il valore della su tutto l'elemento, ma, poiché si integra campionando solo in quattro punti, sarà necessario avere la solo su quei quattro punti. I punti in integrazione sono sempre definiti sul sistema di coordinate locali e hanno un loro punto associato sul sistema e .

Il nostro solutore agli elementi finiti restituisce come risultato lo stato di deformazione e di tensione ai punti di Gauss. Questo ha un serie di ripercussioni: se, ad esempio, il materiale arriva a snervamento, ciò accade solo nei punti di Gauss, perché, campionando le tensioni solamente ai punti di Gauss, i fenomeni legati all'eccesso di stato tensionale, tipo la plasticizzazione, vengono testati solo ai punti in cui si ha il campionamento dello stato tensionale, quindi nei punti di Gauss. Ad esempio, un elemento quadrilatero a 4 nodi full integretion con 4 punti di Gauss può essere parzialmente plasticizzato in 4 stati differenti (1 punto, 2 punti, 3 punti o tutti e 4 i punti). Non esiste una parzializzazione ulteriore dello stato di plasticizzazione, infatti qualora l'elemento plasticizzasse lo farebbe come minimo per un quarto di elemento ovverosia al più in un punto. Tale fenomeno può costituire in alcuni casi un limite nell'analisi fem di un elemento.

L'integrale dell'energia non viene svolto sul dominio e , ossia non è integrazione doppia in e , per due motivi: è difficile trovare gli estremi di integrazione ed è più comodo definire queste relazioni in ed . Invece di integrare in e , si fa un cambio di coordinate corretto con il determinante dello Jacobiano (definisce il fattore di variazione dell'area tra il quadratino e l'equivalente su e , è il rapporto tra le due aree, è il fattore di scala dell'area).


Cambio coordinate.png


E' opportuno prestare particolare attenzione alla definizione di matrice Jacobiana così come ci viene proposta sulla dispensa di "Progettazione Assistita di Strutture Meccaniche".Dalla lettura della stessa, infatti, se ne desume che tale matrice NON È QUELLA UTILIZZATA SUI LIBRI DI TESTO, bensì ne è LA TRASPOSTA!!


Di seguito è indicata la matrice jacobiana così come compare sulla dispensa del Professor Strozzi, mentre la sua trasposta rappresenta la definizione più diffusa:

A questo punto il problema si riduce a calcolare che, per comodità, è un tensore e un vettore di tre componenti:

Nel seguito, si vedrà come si può calcolare questa .

Si ha che: , ,

In realtà, essendo gli spostamenti definiti in ed , risulta molto più comodo calcolare: , , ,

In particolare, si noti che il calcolo è abbastanza semplice: è sufficiente derivare la formula 1.3 pag 163 della Dispensa, trovando una formula lineare piuttosto semplice. Per scriverla in una forma più generica, è possibile utilizzare un approccio matriciale. Si defnisce un vettore con 4 componenti che sono funzione lineare (quindi rappresentabile come prodotto tra matrice-vettore) degli spostamenti nodali.

=

La matrice ha 4 righe e 8 colonne. Il numero di righe è sempre uguale al numero degli elementi del vettore prodotto e il numero di colonne è sempre uguale al numero di elementi del vettore parte della moltiplicazione (fattore). La prima colonna verrà moltiplicata per il termine . La seconda colonna sarà associata a , e così via. Per quanto riguarda , si prenda la formula e la si derivi rispetto a , mentre per quanto riguarda , si prenda la formula e la si derivi rispetto a . Si applichi lo stesso procedimento per e , andando così a riempire al matrice 4x8.

Si noti che non vi è nessun termine misto tra e , quindi tutti i termini di correlazione tra spostamenti nodali di tipo con sono nulli.

Si nomini tale matrice 4x8 in funzione di e .

La forma compatta risulta essere:

=

Si è esplicitato il legame tra le derivate e gli spostamenti nodali.

Il passaggio da derivate di e in e a derivate di e in e (sono quelle necessarie per ottenere le componenti di deformazione) si ottiene passando attraverso la matrice jacobiana, in particolare per la sua inversa. Si ricordi che:

=

Analogamente, per lo spostamento verticale, si ha che:

=

Risulta comodo scrivere le due relazioni matriciali nella seguente formula, accoppiando i vari termini in un unico vettore:

=

La matrice appena definita prende il nome di (J indica che è una matrice che deriva dallo Jacobiano, inv indica che è l'inverso e * indica che è una variazione dello Jacobiano).

Per ricavare il vettore delle deformazioni, bisogna introdurre una matrice 3x4 (non più funzione di ed ), tale per cui:

=

In definitiva, il vettore deformazione sarà definito come:

=

Rifacendosi alla scrittura già vista per elemento triangolare, si ha:

=

con

=

Il vantaggio di tali passaggi è che il legame tra il vettore delle deformazioni e il vettore degli spostamenti nodali è dato dalla matrice .

La matrice di rigidezza dell'elemento sarà:

formula presente anche nella dispensa di PA a pagina 169.

L'integrale viene svolto per via numerica, effettuando la sommatoria delle quattro componenti dell'integrando, valutate di volta in volta nei 4 punti di campionamento, pesati con peso unitario. Sommando i quattro contributi, si ottiene la matrice di rigidezza dell'elemento isoparametrico.

La differenza principale tra l'elemento triangolare a tre nodi e l'elemento isoparametrico a quattro nodi consiste nelle diverse caratteristiche della matrice : nel caso dell'elemento triangolare, la matrice è costante sull'area dell'elemento, ovvero si può scegliere qualsiasi punto dell'area come punto di campionamento. Pertanto, non ci si pone il problema di integrare con i punti di Gauss. Ciò è equivalente a pensare che il triangolo a tre nodi venga integrato numericamente con un punto di Gauss al suo centro o con un punto di Gauss in qualunque punto; essendo la funzione costante, è possibile campionarla in qualunque punto.

Ricavata la matrice di rigidezza e quindi le deformazioni, è possibile risalire alle tensioni attraverso la formula:

Quindi, i campi tensionali sulla struttura vengono risolti elemento per elemento all'interno del modello.


Riduzione ai nodi dei carichi distribuiti

Si consideri un generico elemento (indifferentemente di tipo triangolare a tre nodi o isoparametrico a quattro nodi. Per semplicità, si supponga che tale elemento abbia un lato verticale, lungo . Su tale lato, si applichi una pressione distribuita generica, non costante. Si definisca un asse , tale per cui il carico distribuito possa essere descritto come .

Wiki2.jpg

Tuttavia, siccome i carichi devono essere necessariamente nodali, occorre ricondurre il carico distribuito a carico nodale. Trascurando la parte destra dell'elemento (che potrebbe avere 1 o 2 nodi, a seconda che sia un triangolo a tre nodi o un quadrilatero a quattro nodi), si opera sugli unici 2 nodi disponibili. L'unico vincolo è che lo spostamento sia lineare lungo il bordo. Nel seguito, si presenta una procedura valida nel caso generale e non per uno specifico tipo di elemento. (Si precisa che, nel caso di triangolo a tre nodi, sarebbero sufficienti le equivalenze in termini di risultante e di momento risultante, al fine di ottenere i valori delle forze ai nodi dell'elemento).

La definizione utilizzata nel seguito è a pari lavoro delle forze esterne. Si consideri il lato verticale. La forza agisce perpendicolarmente al lato, quindi è orizzontale. Si trascurino pertanto le componenti di spostamento verticale. Il primo nodo ha uno spostamento orizzontale . Il secondo nodo ha uno spostamento orizzontale . Poiché le funzioni di forma sono di tipo lineare, o, comunque, sono lineari sui bordi (caso dell'elemento isoparametrico a quattro nodi), tutti i punti su cui agisce la pressione distribuita si muovono secondo la seguente legge:

Tale funzione spostamento è una funzione spostamento del bordo, secondo i vincoli imposti dalle funzioni di forma dell'elemento stesso. Si ricorda che ogni elemento ha un modo diverso di ripartire i carichi nodali, in quanto la ripartizione dei carichi nodali passa attraverso le funzioni di forma. Nel caso in esame, la funzione di forma è lineare, quindi il moto dei punti sul bordo caricato segue una legge lineare. Cambiando la funzione di forma dell'elemento, si cambia il legame tra gli spostamenti nodali e gli spostamenti di bordo, quindi cambia la ripartizione delle forze esterne sui nodi.

Si calcoli il lavoro svolto dalle forze esterne per compiere uno spostamento di entità del primo nodo e uno spostamento di entità del secondo nodo:

La quantità rappresenta la forza infinitesima ; è possibile riscrivere la formula in maniera più semplice:

Il lavoro delle forze concentrate nodali vale, supponendo le forze costanti:

I valori calcolati per il lavoro devono essere necessariamente uguali . Di conseguenza, ne deriva che:

;

Nel caso in cui le forze su un nodo derivano da lati di elementi adiacenti, la forza risultante sul nodo è semplicemente la somma vettoriale delle forze nodali calcolate per il singolo elemento.

5-copia.PNG

Calcolo delle reazioni vincolari

Le reazioni vincolari si ricavano, una volta calcolata la matrice di rigidezza utilizzando la formula:

=

tenendo presente che questa relazione è valida in condizioni di pre-vincolamento. Nello specifico, la matrice di rigidezza è la matrice di rigidezza prima delle modifiche associate alla procedura di vincolamento. Analogamente, il vettore dei termini noti è il vettore delle forze nodali prima che intervengano le modifiche associate alla procedura di vincolamento. Lo stesso vale per il vettore incognito . Tale sistema, nel suo complesso, è un sistema composto da equazioni di equilibrio nodale.

Una volta eseguito il vincolamento, si determina una , cioè una versione di K vincolata, tale che:

=

Tale sistema non è più un sistema di equazioni di equilibrio. Più precisamente, non lo è più completamente, dal momento che alcune equazioni di equilibrio sono state soppresse; in particolare, sono state soppresse quelle associate ai gradi di libertà vincolati, per far posto a dei vincoli cinematici.

La soluzione del sistema lineare è data da:

=

Tale simbologia appartiene alla tipica notazione Matlab per definire la soluzione di un sistema lineare di equazioni, in cui si esplicita il non-calcolo della matrice inversa. Si nota il NON utilizzo della seguente forma alternativa:

=

in quanto tale notazione tenderebbe a far pensare che si debba operare un'inversione della matrice , cosa che non si fa realmente mai. Il solutore andrebbe a calcolare l'inversa della matrice, utilizzando un algoritmo che l'esperienza dice essere più gravoso dal punto di vista computazionale.

Una volta ricavato , si procede, ricavando le reazioni vincolari con la seguente formula:

= -

meglio espresse nella forma:

= +

dove la matrice e il vettore delle forze nodali sono entrambi pre-vincolamento. Nel seguito, si illustra la definizione delle reazioni vincolari sopra esplicitata. Se si prende il vettore , lo si premoltiplica per la matrice e si sottrae il vettore , rimane un residuo, in quanto non è più soluzione del sistema = , ma è soluzione del sistema = . Ciò fa sì che - non sia pari a zero, ma sia pari ad un residuo . Si definiscono, pertanto, le reazioni vincolari come tale residuo. In tal modo, le reazioni vincolari appaiono analoghe alle forze esterne. Invece di essere applicate in forma di forza, sono passate dai vincoli e sono forze non note a priori, che agiscono dall'esterno sulla struttura. Tuttavia, la natura delle reazioni vincolari è analoga alle forze esterne, in quanto sono entrambe forze nodali.


Note sugli elementi isoparametrici

  • la definizione di matrice Jacobiana a p.165 della dispensa non coincide con l'usuale definizione della stessa, vedi voce wikipedia o le reference guide dei vari codici di calcolo. In particolare la forma presentata in dispensa risulta, secondo la definizione standard, una Jacobiana trasposta.
  • integrazione gaussiana applicata agli isoparametrici 4 nodi
  • integrazione gaussiana ridotta