Processing math: 100%

Strumenti Utente

Strumenti Sito


wikipaom2015:lez11

ELEMENTI FINITI ISOPARAMETRICI

Gli elementi finiti isoparametrici si basano su una teoria che permette di ottenere elementi a lati curvi; nella pratica, sono particolarmente adatti a ricopiare bordi e superfici curve di componenti meccanici. Il termine isoparametrico significa “ad uguale parametro” e deriva dal fatto che due tipi di funzioni, impiegate in due punti distinti della teoria degli elementi finiti isoparametrici, sono tra loro uguali. Queste due funzioni riguardano il collegamento tra coordinate globali e locali ed il legame tra spostamenti globali e locali. Un ulteriore chiarimento sull’idea di elementi finiti, relativo al concetto di coordinate locali, è quello di considerare un elemento finito (ad esempio quadrilatero) che viene guardato da una lente attraverso la quale risulta come un quadrato con i lati di lunghezza 2 centrato nell’origine di un sistema di riferimento; in questo modo applicare qualsiasi teoria risulterebbe molto semplice. Una volta scritta la teoria è possibile riportare l’elemento finito alla sua forma originale attraverso una”controlente”. (Vedi Figura 1)

SPOSTAMENTI

Nella seguente trattazione del metodo ad elementi finiti si farà riferimento ad un elemento isoparametrico quadrilatero a 4 nodi ottenendo dei risultati che potranno essere utilizzati allo stesso modo per tutti gli altri elementi.

Definendo l'elemento finito del quadrilatero con coordinate parametriche si ottengono tre casi :

  • 1 caso

u(x,y)=a+bx+cy+dx2

Il vantaggio di questo caso è che gli spostamenti vengono descritti meglio lungo asse x che sull’asse y , ma il termine quadratico (dx2)introduce un'anisotropia del modello , tale che non rende l’equazione utilizzabile.

  • 2 caso

u(x,y)=a+bx+cy+dy2

Il fenomeno dell’anisotropia riguarda anche in questo punto ,per la presenza del termine quadratico (dy2)

  • 3 caso

u(x,y)=a+bx+cy+dyx

Il parametro dxy causa la deformazione dei lati in modo non rettilineo creando una discontinuità tra un elemento e l’altro.

Per semplificare la trattazione si costruisce l’elemento con un metodo alternativo.

Si consideri un DOPPIO piano :

- fisico, ( di coordinate x,y ): nel modello fisico è contenuto il modello trattato con coordinate nodali (1,2,3,4) nel senso antiorario.

- A coordinate locali (ξ,η): nel metodo di coordinate locali il quadrilatero è centrato nel punto (0,0) ed ha lato 2.

Per poter associare l’elemento dal piano di coordinate locali al piano fisico serve una funzione di mappatura BIUNIVOCA ED INVERTIBILE x(ξ,η),y(ξ,η) . Tale funzione non è fornita in modo esplicito, ma è costruibile.

Per definire la legge di mappatura si definiscono quattro funzioni di forma sul piano di coordinate locali:

Nn(ξ,η)=a+bξ+cη+dξη con n= numero di nodi.

Tali funzioni assumono valore 1 sul nodo indicato dal pedice e 0 in tutti gli altri. Sono di forma bilineare in ξ,η, non in x, y, ovvero sono lineari con ξoη costanti, mentre sono non lineari (ovvero quadratiche) in caso contrario. Nel caso in esame con 4 nodi ottengo le seguenti funzioni di forma:

N1(ξ,η)=(ξ1)(η1)4(1)

N2(ξ,η)=(ξ+1)(η1)4(2)

N3(ξ,η)=(ξ+1)(η+1)4(3)

N4(ξ,η)=(ξ1)(η+1)4(4)

ESEMPIO

Se si volesse disegnare, ad esempio, la funzione di forma associata al nodo 1 (N1) si otterrebbe la Figura 2.

La Funzione di forma vale 1 in corrispondenza al nodo 1 e 0 sugli altri tre nodi; inoltre si noti che, essendo la funzione bilineare, i 4 lati del dominio in coordinate locali sono lati in cui una delle due coordinate ha valori costanti. Infine avendo costruito la superficie da essa si possono ricavare i valori della funzione in qualunque punto del dominio, come mostrato in figura.

Quindi la funzione di mappatura per le coordinate x e y è quindi:

x(ξ,η)=x1N1(ξ,η)+x2N2(ξ,η)+x3N3(ξ,η)+x4N4(ξ,η)(5)

y(ξ,η)=y1N1(ξ,η)+y2N2(ξ,η)+y3N3(ξ,η)+y4N4(ξ,η)(6)

Espandendo, ad esempio, la (5) si ottiene:

x(ξ,η)=4i=1xiNi(ξ,η)(7)

Viene definito un vettore colonna x_ una matrice N_, ottenendo:

x(ξ,η)=ˉxTˉN(ξ,η)(8)

con:

ˉx=[x1x2x3x4]

ˉN(ξ,η)=[N1(ξ,η)N2(ξ,η)N3(ξ,η)N4(ξ,η)]

La formulazione (8) apre immediatamente all’estensione a funzioni di forma di ordine maggiore. Ad esempio si può definire 8 funzioni di forma opportune per ottenere un elemento isoparametrico ad 8 nodi, come mostrato in Figura 3.

FIXME la funzione di forma deve essere nulla in corrispondenza del nodo 6

Si definisce il CAMPO DEGLI SPOSTAMENTI in funzione del piano di coordinate locali:

u(ξ,η)=4i=1uiNi(ξ,η)(9)

v(ξ,η)=4i=1viNi(ξ,η)(10)

che in forma matriciale diventano :

u(ξ,η)=u_TN_(ξ,η)

v(ξ,η)=v_TN_(ξ,η)

con:

v_=[v1v2v3v4]

u_=[u1u2u3u4]

Il grande vantaggio della formulazione precedente è che gli spostamenti hanno un andamento lineare lungo i lati di un elemento per cui, se il lato di un elemento nasce rettilineo ed ai punti di quel lato applico uno spostamento lineare, un lato rettilineo rimane rettilineo.

Poichè queste considerazioni valgono per tutti gli elementi della mesh, si esclude la possibilità che due elementi che condividono un lato si spostino in maniera differenziata sui punti del lato stesso.

Gli spostamenti alla frontiera sono continui poichè valgono due ipotesi:

  • essendo x ed y funzioni lineari sul lato allora sono lineari anche le loro inverse, ovvero ξedη. Essendo gli spostamenti u(ξ(x,y),η(x,y))ev(ξ(x,y),η(x,y)) funzione di funzioni lineari allora sono anch’essi lineari.
  • Gli spostamenti nodali sono univoci ovvero indipendenti dall’elemento utilizzato come riferimento.

Combinando quanto detto non solo si ha la continuità degli spostamenti dei nodi, ma si ottiene anche la garanzia di poter affiancare un elemento isoparametrico 4 nodi ad uno triangolare 3 nodi senza avere discontinuità di spostamenti sull’interfaccia.

Si noti come il termine isoparametrico venga spiegato dal fatto che si usano le stesse funzioni di forma per la mappatura delle coordinate x,y,ξedη e per gli spostamenti u(ξ,η)ev(ξ,η)

DEFORMAZIONI

Si vogliono calcolare le deformazioni generalizzate εx,εyeγxy a partire dagli spostamenti globali u e v.

εx=ux

εy=vy

γxy=uy+vx

Essendo gli spostamenti funzioni composte verrebbe automatico scrivere

ux=uξξx+uηηx(11)

Il problema è che le derivate ξxeηx non sono facili da calcolare in quanto è disponibile solo x in funzione di ξedη, non il viceversa. Si decide quindi di utilizzare il seguente approccio solo sugli spostamenti u, quelli per v sono analoghi.

{uξ=uxxξ+uyyξuη=uxxη+uyyη(12)

In questo modo si esplicitano delle identità in cui i termini sono comodi da calcolare a partire dalla relazione seguente:

uξ=iuiNi(ξ,η)ξ(13)

Gli unici due termini non facili da calcolare nella relazione (12) sono ux,uy:

Per trovarli è possibile utilizzare un approccio matriciale definendo il vettore di 2 componenti che sono funzione lineare degli spostamenti nodali: {uξuη}=[xξyξxηyη]{uxuy}=J{uxuy}(14)

In cui la matrice viene chiamata Jacobiana (J) ed è moltiplicata per il vettore delle incognite.

NOTA BENE: E' opportuno notare che in programi quali Maxima e Mathematica la matrice Jacobiana scritta in questa forma viene indicata come Matrice Jacobiana trasposta.

Se la matrice J è invertibile le derivate degli spostamenti di u e v rispetto x ed y valgono: {uxuy}=J1{uξuη}(15)

Similmente {vxvy}=J1[vξvη](16)

Dove: J1=1|J|[J22J12J21J11](17)

con

J11=xξ

J12=yξ

J21=xη

J11=yη

Quindi passando per l'inversa della matrice Jacobiana posso ricavare, partendo dalla forma degli spostamenti, le derivate ux,uy,vxevy che servono per calcolare le deformazioni:

[uxuyvxvy]=[J1__00J1__][uξuηvξvη](18)

La matrice appena definita prende i nome di Jinv__(ξ,η) (J indica che è una matrice che deriva dallo Jacobiano, inv indica che è l'inversa e * indica che è una variazione dello Jacobiano).

Mentre il vettore di 4 componenti che è funzione lineare degli spostamenti nodali si calcola [uξuηvξvη]=[η1401η40η+140η140ξ140ξ140ξ+1401ξ40 0η1401η401+η401η4 0ξ140ξ140ξ+1401ξ4]Q__[u1v1u2v2u3v3u4v4](19)

Infine per ricavare il vettore delle deformazioni, bisogna introdurre una matrice H 3×4 (non più funzione di ξedη) tale per cui : [εxεyγxy]=[100000010110][uxuyvxvy](20)

In definitiva i vettore deformazione epsilon sarà definito come:

ε_=H__Jinv__(ξ,η)Q__(ξ,η)B__δ_(21)

con B__(ξ,η) la matrice di operatori differenziali che lega gli spostamenti nodali alle deformazioni

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

ε_=B__(ξ,η)δ_(22)

SFORZI

Si parte dal calcolo dell'energia potenziale elastica sull'elemento :

U=AREA12ε_Tσ_dA=12AREA[δ_TB__T(ξ,η)D__B__(ξ,η)δ_]dA(23)

Portando fuori le quantità costanti si ottiene:

U=12δ_TAREA[B__T(ξ,η)D__B__(ξ,η)]dAK__δ_(24)

Dove K__ è la matrice di rigidezza dell'elemento.

Il problema è che B e BT sono funzione di ξ ed η e non di x , y ; naturalmente verrebbe logico integrare in dx e dy ma queste coordinate complicano l'integrazione. Quindi si decide di declinare il dA in funzione di ξ ed η (dimostrazione in Appendice 1):

dA=|J(ξ,η)|dξdη(25)

Ottenendo:

U=12δ_Ta[B__T(ξ,η)D__B__(ξ,η)|J__|dξdη]δ_=12δ_T1111[B__T(ξ,η)D__B__(ξ,η)|J__|dξdη]δ_(25)

Essendo l'integrando non costante sull'elemento si decide di passare ad un'integrazione numerica gaussiana a 4 punti (di coordinate ξ,η=±13) ottenendo:

K__ele4j=1wjB__T(ξj,ηj)D__B__(ξj,ηj)|J__(ξj,ηj)|(26)

considerando:

wj peso di ogni punto

\underline{\underline{D}} che indica l'uniformità del materiale e che in genere si suppone costante, ma può anche non esserlo se l'elemento non è uniforme.

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

σ=D__ε(27)

Appendice 1

Se si è sul dominio ξ ed η prendo un elemento di area (a), presa nell'intorno di uno specifico punto (ξ;η), e trasformata sul piano x ed y in un oggetto di forma generica la cui area(a') è:

a=|J(ξ,η)|a

Il quadrilatero dξ,dη viene trasformato in un quadrilatero di coordinate:

1:x(ξ,η);y(ξ,η)

2:x(ξ+dξ,η);y(ξ+dξ,η)

3:x(ξ+dξ,η+dη);y(ξ+dξ,η+dη)

4:x(ξ,η+dη);y(ξ,η+dη)

Trascurando gli infinitesimi di ordine superiore:

2:≅x(ξ,η)+xξ|ξ,ηdξ;y(ξ,η)+yξ|ξ,ηdξ

3:≅x(ξ,η)+xξ|ξ,ηdξ+xη|ξ,ηdη;y(ξ,η)+yξ|ξ,ηdξ+yη|ξ,ηdη

4:≅x(ξ,η)+xη|ξ,ηdη;y(ξ,η)+yη|ξ,ηdη

Se si considera un troncamento al primo ordine della trasformazione precedente risulta che il quadrato viene trasformato in un parallelogramma costituito da due triangoli di aree uguali (vedi Figura 5 ):

L'area del quadrato evidenziata (a) vale:

a=12dξdη(1.1)

Poichè il parallelogramma è costituito da due triangoli di aree uguali si considera solo quello evidenziato e si calcola l'area con la regola del determinante: A=12|111x_x_+x_ξ|ξ,ηdξx_+x_η|ξ,ηdηy_y_+y_ξ|ξ,ηdξy_+y_η|ξ,ηdη|(1.2)

Dalle proprietà dei determinanti si sa che : “se si aggiunge ad una riga di una matrice una combinazione lineare delle rimanenti righe il determinante non cambia” ; quindi si sottrae alla seconda riga la prima moltiplicata per x_, in seguito si sottrae alla terza riga la prima moltiplicata per y segnato ottenendo:

a=12|1110x_ξ|ξ,ηdξx_η|ξ,ηdη0y_ξ|ξ,ηdξy_η|ξ,ηdη|(1.3)

Sapendo che “detta G' la matrice ottenuta moltiplicando una riga (o una colonna) per un elemento h di K, si ha : det G' = h det G” si possono portare fuori i termini dξ e dη dalla matrice ottenendo la matrice jacobiana:

a=12[x_ξ|ξ,ηx_η|ξ,ηy_ξ|ξ,ηy_η|ξ,η]J__(ξ,η)dξdη(1.4)

In conclusione si può quindi affermare che :

a=|J(ξ,η)|a

Bibliografia

1) “PROGETTAZIONE ASSISTITA DI STRUTTURE MECCANICHE” dalle lezioni del Prof.Antonio Strozzi

2) link a “PROPRIETA' DEI DETERMINANTI” : pdf proprietà determinante citate a lezione

wikipaom2015/lez11.txt · Ultima modifica: 2016/02/03 08:34 da ebertocchi