Metoda Elementelor Finite

158
Dan-Sorin Com¸ sa METODA ELEMENTELOR FINITE Curs introductiv Editura U.T. PRES Cluj-Napoca, 2007

Transcript of Metoda Elementelor Finite

Page 1: Metoda Elementelor Finite

Dan-Sorin Comsa

METODA ELEMENTELOR FINITE

Curs introductiv

Editura U.T. PRES

Cluj-Napoca, 2007

Page 2: Metoda Elementelor Finite
Page 3: Metoda Elementelor Finite

Prefata

In ultimele decenii, metoda elementelor finite (MEF) a devenitun instrument de calcul numeric tot mai des utilizat de catre in-gineri. Elaborata initial ca procedeu de rezolvare aproximativa aproblemelor de mecanica structurilor (ansambluri de bare, grinzisi/sau placi), MEF si-a extins cu repeziciune aria de aplicabilitate.Deja la mijlocul anilor 1970 era demonstrata capacitatea acesteimetode de a servi la abordarea numerica a oricarei probleme expri-mabile sub forma unui sistem de ecuatii diferentiale sau cu derivatepartiale. Sunt o multitudine de aplicatii ingineresti care apeleaza laun asemenea limbaj. Le vom aminti aici doar pe cele cateva care sedisting printr-o deosebita importanta:

• probleme de mecanica solidelor (deformabile elastic, elasto-plastic etc.);

• probleme de mecanica fluidelor;

• probleme de transfer termic;

• probleme de electricitate si magnetism.

Extinderea universului aplicativ al metodei a stimulat elaborareaa numeroase programe de analiza cu elemente finite. Sub aspectulfunctionalitatii, aceste programe se dovedesc totusi greu de stapanitatunci cand utilizatorul nu dispune de o pregatire minimala. Chiarmai grav, lipsa cunostintelor poate conduce la interpretari completeronate ale unui rezultat numeric.

iii

Page 4: Metoda Elementelor Finite

iv

Aspectele anterior amintite, dar si preocuparile de natura di-dactica au motivat elaborarea acestei carti. Prin continut, ea seadreseaza ın special studentilor de la sectiile ingineresti cu profilmecanic, a caror programa de ınvatamant prevede parcurgerea unuicurs introductiv de metoda elementelor finite. Orientandu-ne spreacest grup de cititori, am cautat sa mentinem expunerea la un nivelcat mai accesibil. Notiunile de matematica pe care le presupunemcunoscute se regasesc ın programa oricarei facultati de inginerie.De asemenea, exemplificarile din tot cuprinsul lucrarii apeleaza laun bagaj de informatie tehnica pe care studentii ıl poseda la fineleprimului an universitar. Desi a fost conceputa ca material de curs,cartea ofera numeroase indicatii aplicative. Fiecare exemplu esteurmat de comentarii si sfaturi destinate celor interesati de aspec-tele practice ale analizei cu elemente finite.

Dan-Sorin Comsa

Cluj-Napoca, ianuarie 2007

Page 5: Metoda Elementelor Finite

Cuprins

1 Metoda elementelor finite. Prezentare generala 1

1.1 Rezolvarea numerica a problemelor ingineresti . . . 1

1.2 Prezentarea notiunilor de baza ale metodei elemente-lor finite folosind ca exemplu modelul unei structuride tip cadru plan articulat . . . . . . . . . . . . . . 7

1.2.1 Formularea problemei-exemplu. Modelul cuelemente finite asociat . . . . . . . . . . . . 8

1.2.2 Ecuatiile care descriu comportarea unui ele-ment al retelei . . . . . . . . . . . . . . . . . 12

1.2.3 Sistemul de ecuatii care descrie raspunsul me-canic al retelei de elemente finite . . . . . . 20

1.2.4 Starea de eforturi proprie fiecarui element alretelei . . . . . . . . . . . . . . . . . . . . . 25

1.3 Concluzii . . . . . . . . . . . . . . . . . . . . . . . . 26

2 Caracteristici distinctive ale metodei elementelorfinite 29

2.1 Formularea problemei-exemplu. Modelul sau analitic 31

2.2 Modelul cu elemente finite atasat problemei-exemplu 35

2.3 Rezolvarea numerica a modelului cu elemente finite 41

2.4 Concluzii . . . . . . . . . . . . . . . . . . . . . . . . 50

3 Prezentarea catorva tipuri uzuale de elemente finite 53

3.1 Clasificare . . . . . . . . . . . . . . . . . . . . . . . 54

v

Page 6: Metoda Elementelor Finite

vi

3.1.1 Elemente finite unidimensionale . . . . . . . 553.1.2 Elemente finite bidimensionale . . . . . . . . 563.1.3 Elemente finite tridimensionale . . . . . . . 59

3.2 Constructia aproximantelor pentru cateva tipuri deelemente finite . . . . . . . . . . . . . . . . . . . . . 613.2.1 Elementul triunghiular liniar . . . . . . . . . 633.2.2 Elementul patrulater biliniar . . . . . . . . . 703.2.3 Elementul tetraedric liniar . . . . . . . . . . 753.2.4 Elementul hexaedric triliniar . . . . . . . . . 81

3.3 Aproximarea marimilor vectoriale . . . . . . . . . . 873.4 Expandarea aproximantelor elementale . . . . . . . 90

3.4.1 Expandarea aproximantelor unei marimi detip scalar . . . . . . . . . . . . . . . . . . . 91

3.4.2 Expandarea aproximantelor unei marimi detip vectorial . . . . . . . . . . . . . . . . . . 95

3.5 Concluzii . . . . . . . . . . . . . . . . . . . . . . . . 98

4 Modelul cu elemente finite al problemelor de elasti-citate liniara 994.1 Formularea problemelor de elasticitate liniara . . . 1004.2 Construirea modelului cu elemente finite al proble-

melor de elasticitate liniara utilizand teorema de mi-nim a energiei potentiale . . . . . . . . . . . . . . . 107

4.3 Probleme bidimensionale de elasticitate . . . . . . . 1234.3.1 Probleme de stare plana de tensiuni . . . . . 1244.3.2 Probleme de stare plana de deformatii . . . 1284.3.3 Probleme cu simetrie axiala . . . . . . . . . 131

4.4 Exemplu de utilizare a metodei elementelor finite larezolvarea unei probleme de elasticitate . . . . . . . 135

4.5 Concluzii . . . . . . . . . . . . . . . . . . . . . . . . 146

A Formule de derivare matriceala 147

Bibliografie 151

Page 7: Metoda Elementelor Finite

Capitolul 1

Metoda elementelor finite.Prezentare generala

1.1 Rezolvarea numerica a problemelor

ingineresti

Calculul ingineresc este o componenta esentiala a procesu-lui de proiectare. Orice echipament trebuie realizat ın asa felıncat functionarea lui sa fie sigura, ın conditiile ındeplinirii unorcerinte adeseori contradictorii ın ceea ce priveste costul, dimensiu-nile de gabarit, tehnologia de executie, fiabilitatea etc. Satisfacereaexigentelor anterior mentionate duce la impunerea unor restrictiipe care calculul ingineresc trebuie sa le ındeplineasca (de exemplu,nedepasirea unor tensiuni, deformatii sau deplasari maxime, ga-rantarea unei anumite rezistente la oboseala etc.). Calculele furni-zeaza informatii din care proiectantul poate deduce cat de aproapesunt componentele structurale de stari limita care ar putea periclitafunctionalitatea ıntregului ansamblu.

Marea diversitate a problemelor pe care inginerul trebuie sa lerezolve ın vederea proiectarii unui echipament se reflecta ın multi-tudinea metodelor de calcul utilizate la ora actuala. Aceste metode

1

Page 8: Metoda Elementelor Finite

2 Capitolul 1

s-au dezvoltat treptat, pe masura acumularii de cunostinte teore-tice si tehnologice si, ın ultima jumatate de secol, concomitent cuevolutia calculatoarelor electronice.

O prima clasificare a metodelor de calcul folosite ın practicainginereasca distinge doua mari categorii:

• metode exacte (sau analitice);

• metode aproximative.

Metodele exacte sunt aplicabile numai pentru rezolvarea unuinumar relativ redus de probleme simple. Restrictiile de utilizaresunt ın general introduse de geometria modelului studiat si de tipulconditiilor la limita. Metodele aproximative sunt folosite la rezol-varea problemelor mai complexe, pentru care nu se poate obtine osolutie analitica.

Exista o mare diversitate de metode aproximative pentru re-zolvarea problemelor ingineresti. Oricare ar fi metoda adoptata deproiectant, ea trebuie sa furnizeze o solutie suficient de precisa pen-tru problema practica analizata.

La ora actuala, ın ingineria mecanica sunt folosite preponderentmetodele aproximative de tip numeric. Din categoria acestora, mairaspandite ca arie de utilizare sunt:

• metoda diferentelor finite (MDF);

• metoda elementelor finite (MEF);

• metoda elementelor de frontiera (MEFr).

La prima vedere, cele trei metode par foarte diferite una de cealalta.La o analiza mai atenta, se constata ınsa ca ele sunt strans legateprin fundamentul matematic (toate apeleaza la notiuni din teoriaecuatiilor diferentiale, calculul variational sau metoda reziduurilorponderate). In continuare vom face o prezentare succinta a aces-tor metode numerice, pentru a le compara din punct de vedere alavantajelor si dezavantajelor specifice.

Page 9: Metoda Elementelor Finite

Rezolvarea numerica a problemelor ingineresti 3

Metoda diferentelor finite (MDF)

Prima metoda larg utilizata pentru rezolvarea unor problemede calcul ingineresc a fost metoda diferentelor finite. MDF este celmai simplu procedeu de rezolvare numerica a sistemelor de ecuatiidiferentiale. In demersul sau, MDF aproximeaza derivatele necunos-cutelor prin diferente finite [6]. Pentru aplicarea practica a MDF,domeniul spatial ocupat de sistemul fizic analizat este acoperit cuo retea rectangulara de noduri. Ecuatiile diferentiale care descriucomportarea sistemului sunt aproximate cu diferente finite calcu-late ın noduri. Se obtine astfel un set de ecuatii algebrice, pentrurezolvarea caruia se pot folosi procedee numerice consacrate (deexemplu, metoda eliminarii Gauss [3]). Solutia aproximativa estereprezentata de valorile necunoscutelor ın noduri. Calitatea solutieieste cu atat mai buna cu cat reteaua de noduri este mai densa.

Principalele avantaje ale MDF sunt simplitatea conceptuala siusurinta transpunerii ın programe de calcul. Totusi, ea prezinta siunele dezavantaje care ıi limiteaza serios aria de aplicabilitate:

• metoda nu furnizeaza decat valorile necunoscutelor ın noduri,neoferind informatii despre variatia acestora ıntre noduri;

• discretizarea unor corpuri de forma complexa folosind numairetele rectangulare conduce de multe ori la aproximari ne-satisfacatoare ın zonele de colt sau cu variatii importante alesectiunii;

• metoda creeaza dificultati ın introducerea unor conditii lalimita complexe.

Din cauza acestor dezavantaje, MDF se aplica astazi mai multpentru rezolvarea problemelor de transfer termic sau de mecanicafluidelor si mai putin ın elasticitate sau elastoplasticitate.

Page 10: Metoda Elementelor Finite

4 Capitolul 1

Metoda elementelor finite (MEF)

Metoda elementelor finite este, la ora actuala, cel mai raspanditprocedeu de rezolvare numerica a problemelor ingineresti. La apli-carea MEF, domeniul ocupat de sistemul fizic supus analizei estediscretizat ın subdomenii de dimensiune finita marginite de fron-tiere rectilinii sau curbilinii. Prin aceasta operatie, sistemul real esteınlocuit cu o retea de asa-numite elemente finite [4, 5, 8, 9, 10].Ecuatiile diferentiale care descriu comportarea sistemului vor fiverificate ın medie pe fiecare element. Constructia matematica aelementelor finite asigura un anumit grad de continuitate al aproxi-mantei la traversarea frontierei dintre elemente. Continuitatea esterealizata cu ajutorul unor puncte remarcabile asociate elementelor(cunoscute sub numele de noduri). De fapt, aproximanta solutieiexacte a problemei rezulta ca o functie de valorile necunoscutelorın respectivele noduri. Prin aplicarea MEF se obtine un sistem deecuatii care se rezolva numeric ın raport cu valorile necunoscutelornodale. Din cauza ca aproximantele solutiei exacte verifica ın medieecuatiile diferentiale si ın interiorul elementelor, ele pot fi folositepentru obtinerea unor informatii referitoare la variatia necunoscu-telor problemei ın tot domeniul analizat (nu numai ın noduri, asacum se ıntampla ın cazul utilizarii MDF).

Principalele avantaje ale MEF sunt urmatoarele:

• flexibilitatea (prin faptul ca permite discretizarea unor cor-puri de forma oricat de complexa si manipularea naturala aunor conditii la limita dintre cele mai diverse);

• posibilitatea de a modela corpuri neomogene din punct devedere al proprietatilor fizice;

• usurinta implementarii ın programe de calcul generale.

Cel mai important dezavantaj al MEF consta ın volumul mare dedate de intrare necesar pentru construirea si rezolvarea numerica

Page 11: Metoda Elementelor Finite

Rezolvarea numerica a problemelor ingineresti 5

a modelului. Cea mai mare pondere o detin datele de intrare refe-ritoare la configuratia retelei de elemente finite (coordonatele no-durilor, respectiv asocierea dintre elemente si noduri). Programelede analiza cu elemente finite existente pe piata la ora actuala de-greveaza utilizatorul de sarcina neplacuta a discretizarii manuale,transferand-o unor module specializate care automatizeaza aceastaoperatie.

In ultimele trei decenii au aparut multe asemenea programecomerciale. Cele mai multe dintre ele sunt interfatate cu programelede proiectare asistata de calculator, astfel ıncat utilizarea lor decatre ingineri a devenit relativ simpla.

Metoda elementelor de frontiera (MEFr)

Metoda elementelor de frontiera este un procedeu de rezolvarenumerica a problemelor de calcul ingineresc elaborat relativ recent.La baza MEFr se afla ideea de a atasa unui sistem de ecuatiidiferentiale definit pe un anumit domeniu o ecuatie integrala echi-valenta definita numai pe frontiera respectivului domeniu [2, 7]. For-mularea integrala pe frontiera aduce avantaje importante, fiindcanu mai impune discretizarea interiorului domeniului analizat. Sereduce astfel substantial numarul necunoscutelor problemei si vo-lumul datelor de intrare. Din punct de vedere al utilizarii practice,MEFr se aseamana mult cu metoda elementelor finite, ın sensul cafrontiera este discretizata folosind elemente construite dupa prin-cipii asemanatoare cu cele din MEF. Se obtine tot un sistem deecuatii care se rezolva numeric ın raport cu valorile necunoscutelordin nodurile elementelor de pe frontiera. Odata obtinute valorilenecunoscutelor nodale, se pot determina valorile acestora ın oricepunct din interiorul domeniului folosind un set de formule specifice.

In afara de avantajul anterior mentionat, MEFr mai are o seriede alte caracteristici pozitive:

• se poate aplica ın cazul unor domenii care se extind la infinitdupa una sau mai multe directii;

Page 12: Metoda Elementelor Finite

6 Capitolul 1

• permite obtinerea unor solutii precise ın zonele de colt saucu variatii rapide ale geometriei fara a necesita o densificareexcesiva a retelei (spre deosebire de MEF care impune demulte ori rafinarea pronuntata a retelei ın asemenea portiuniale domeniului analizat);

• poate fi integrata cu usurinta ın programele de proiectare asis-tata de calculator.

MEFr are si cateva dezavantaje care ıi reduc simtitor aria deaplicabilitate:

• dificultati serioase ın modelarea corpurilor neomogene dinpunct de vedere al proprietatilor fizice;

• scaderea preciziei rezultatelor numerice ın cazul corpurilorcare prezinta discrepante majore ale dimensiunilor pe diferitedirectii (de exemplu, bare, grinzi, placi sau ınvelitori).

Totusi, trebuie tinut cont de faptul ca MEFr este ınca la ınceputuriledezvoltarii sale, iar cercetarile viitoare vor ameliora multe dintredeficientele mentionate mai sus.

Din cele prezentate ın acest paragraf putem concluziona ca, laora actuala, metoda elementelor finite este cel mai adecvat procedeude rezolvare numerica a problemelor ingineresti. Raspandirea MEFın toate domeniile proiectarii a fost facilitata de generalitatea siflexibilitatea sa, precum si de numarul mare de programe comercialeinterfatate cu pachetele de proiectare asistata de calculator.

Metoda elementelor finite este departe de a-si fi atins limiteledezvoltarii. In prezent sunt ın desfasurare cercetari care vizeazaconstruirea unor noi tipuri de elemente care sa vina ın ıntampinareacelor mai diverse cerinte practice. De asemenea, apar continuu noiprograme comerciale care tind sa exploateze la maximum facilitatilede procesare ale calculatoarelor moderne (cum ar fi calculul paralel)si sa asigure o integrare cat mai naturala a analizei cu elementefinite ın procesul de proiectare asistata de calculator.

Page 13: Metoda Elementelor Finite

Notiuni de baza ale metodei elementelor finite 7

1.2 Prezentarea notiunilor de baza ale

metodei elementelor finite folosind

ca exemplu modelul unei structuri

de tip cadru plan articulat

Cadrele articulate sunt structuri familiare multor categorii deingineri. Din aceasta cauza, vom folosi problema determinarii de-plasarilor si a starii de eforturi dintr-un asemenea cadru pentruprezentarea notiunilor fundamentale cu care opereaza MEF. Simpli-tatea exemplului confera multa naturalete procesului de construirea modelului cu elemente finite. Vom vedea ca notiuni cum sunt celede element finit, respectiv nod au o semnificatie foarte clara: ele-mentele sunt pur si simplu barele care formeaza cadrul, iar nodurilesunt articulatiile prin care sunt legate barele ıntre ele.

Vom aborda ın continuare urmatoarele aspecte specifice aplicariimetodei elementelor finite:

• procesul de discretizare (definirea elementelor, a nodurilor sistabilirea unei numerotari a acestora);

• sistemul de coordonate global, respectiv sistemul de coordo-nate asociat elementelor (asa-numitul sistem de coordonatelocal);

• variabile globale, respectiv variabile locale ale modelului cuelemente finite;

• ecuatiile care descriu comportamentul unui element finit (maiconcret, comportamentul elastic al unei bare oarecare) expri-mate ın termenii variabilelor locale;

• trecerea de la variabile locale la variabile globale ın ecuatiilecare descriu comportarea unui element;

Page 14: Metoda Elementelor Finite

8 Capitolul 1

• asamblarea ecuatiilor care descriu comportarea elementelorindividuale ıntr-un sistem de ecuatii care descrie comportareastructurii analizate (sistemul de ecuatii global);

• introducerea conditiilor la limita (efectul reazemelor);

• rezolvarea sistemului de ecuatii global si prelucrarea rezulta-telor (ın vederea obtinerii starii de eforturi si deformatii dinfiecare element).

Abordarea problemei-exemplu va fi foarte intuitiva. Vom face apella cunostinte elementare de mecanica aplicata si rezistenta mate-rialelor. O prezentare mai riguroasa din punct de vedere matematicva face obiectul capitolelor urmatoare ale cursului.

1.2.1 Formularea problemei-exemplu. Modelulcu elemente finite asociat

Ne propunem sa determinam starea de eforturi si deformatiice apare ın cadrul din figura 1.1 datorita solicitarii lui de catreforta exterioara P = 1414 N, al carei suport este ınclinat la 45

fata de orizontala. Vom considera ca toate barele cadrului suntconfectionate dintr-un acelasi material metalic avand o compor-tare liniar elastica (descrisa prin modulul de elasticitate longitudi-nal E = 2 × 1011 Nm−2). Pentru a nu complica excesiv modelulcu elemente finite, consideram ca toate barele au aceeasi sectiunetransversala, a carei arie este A = 5 × 10−4 m2. In figura 1.1 seobserva ca lungimea fiecarei bare se obtine pornind de la un moduldimensional ℓ = 1 m. Vom mai admite ca articulatiile A, B, C siD ale cadrului sunt lipsite de frecare. Aceasta ipoteza este frecventutilizata ın calculele ingineresti. Ea are drept consecinta faptul cafiecare bara este solicitata numai de forte axiale.

Notiunea fundamentala cu care opereaza MEF este cea de ele-

ment. In exemplul de fata, fiecare bara a cadrului va reprezenta un

Page 15: Metoda Elementelor Finite

Notiuni de baza ale metodei elementelor finite 9

A B

CD

P

45

ℓ ℓℓ√

2

P = 1414 Nℓ = 1 mA = 5 × 10−4 m2

E = 2 × 1011 Nm−2

Figura 1.1: Cadru plan articulat

Articulatie Bara Articulatie

Nod Element Nod

Figura 1.2: Echivalenta dintre o bara si elementul corespunzatordin model

Page 16: Metoda Elementelor Finite

10 Capitolul 1

astfel de element. Multimea tuturor elementelor care formeaza mo-delul se numeste retea. O alta notiune importanta a MEF este cea denod. Din punct de vedere matematic, nodurile sunt puncte remarca-bile asociate elementelor si folosite pentru construirea functiilor careaproximeaza necunoscutele modelului. In cazul de fata, nodurile auo semnificatie fizica foarte simpla: ele sunt centrele articulatiilorprin care sunt interconectate barele cadrului. Figura 1.2 arataechivalenta dintre o bara reala si elementul folosit ın model.

Analizand figura 1.1, constatam ca avem cinci elemente si patrunoduri. Pentru a identifica ın mod unic aceste entitati vom utilizao numerotare separata a nodurilor, respectiv elementelor (fig. 1.3).Deocamdata nu ne punem problema folosirii unei strategii de nu-merotare anume. Vom vedea ın capitolele urmatoare ca, ın timp ceordinea numerotarii elementelor nu are importanta, ordinea nume-rotarii nodurilor influenteaza mult eficienta realizarii calculelor decatre programele MEF. Din fericire, programele existente la ora ac-tuala pe piata poseda module specializate care realizeaza automatdiscretizarea ın elemente finite si genereaza o numerotare optimalaa nodurilor.

Peste tot ın cadrul cursului vom utiliza urmatoarea conventiede notare (fig. 1.3):

• numerele de ordine ale elementelor vor fi ıncadrate de paran-teze rotunde;

• numerele de ordine ale nodurilor vor fi prezentate simplu (faraa fi deci ıncadrate de alte simboluri).

In felul acesta nu va exista pericolul confuziei ıntre cele doua nu-merotari.

Corespondenta dintre fiecare element si nodurile asociate estedescrisa de asa-numita matrice de conexiuni a retelei de elementefinite (vezi tabelul 1.1 si figura 1.3). Prin utilizarea matricii de co-nexiuni, programele de analiza cu elemente finite dobandesc o mareflexibilitate. Intr-adevar, odata generata aceasta matrice, pentru

Page 17: Metoda Elementelor Finite

Notiuni de baza ale metodei elementelor finite 11

(3)

(2)

(4)

(1) (5)

1

2

3

4

x

y

P

45

Y1

X1

Y3

u2, v2

U2, V2

u4, v4

U4, V4

u1, v1

U1, V1

u3, v3

U3, V3

Figura 1.3: Numerotarea elementelor si a nodurilor, respectiv siste-mul de coordonate global (x, y) ce va fi utilizat ın model

Tabelul 1.1: Matricea de conexiuni a retelei de elemente finite dinfigura 1.3

Element Nod Nod(e) i j

(1) 1 2(2) 1 3(3) 2 3(4) 2 4(5) 3 4

Page 18: Metoda Elementelor Finite

12 Capitolul 1

Tabelul 1.2: Coordonatele carteziene ale nodurilor din figura 1.3

Nod xk yk

k

1 0 02 0 ℓ3 ℓ 04 ℓ ℓ

precizarea completa a geometriei unei retele nu mai sunt necesaredecat coordonatele nodurilor.

Pentru a specifica respectivele coordonate, trebuie ales un sistemde referinta global. Acest sistem poate fi pozitionat oricum ın raportcu modelul. In general, utilizatorii MEF orienteaza axele sistemuluide referinta global ın maniera cea mai convenabila din punct devedere al geometriei modelului (urmarind ın principal exploatareasimetriilor acestuia si determinarea cat mai usoara a coordonatelornodale). In figura 1.3 este prezentat sistemul de referinta globalasociat problemei-exemplu. Se observa ca, prin aceasta orientare aaxelor, coordonatele nodurilor rezulta foarte simplu (tab. 1.2).

1.2.2 Ecuatiile care descriu comportarea unuielement al retelei

Elementele care formeaza modelul prezentat ın paragraful ante-rior sunt bare articulate realizate dintr-un material liniar elastic. Elesunt solicitate numai de forte axiale. Ecuatiile care descriu compor-tarea mecanica a elementelor vor fi deduse folosind legea lui Hookescrisa pentru cazul solicitarii uniaxiale.

Exprimarea acestor ecuatii direct ın sistemul de coordonateglobal este incomoda, fiindca elementele retelei din figura 1.3 auorientari diferite. Este mult mai convenabil sa se ataseze fiecarui ele-

Page 19: Metoda Elementelor Finite

Notiuni de baza ale metodei elementelor finite 13

ment un sistem de coordonate propriu (un asa-numit sistem de co-

ordonate local) orientat cat mai favorabil. Ecuatiile mecanice vor fiexprimate ın acest sistem, ıntr-o faza urmatoare ele fiind transferateın sistemul de coordonate global printr-o simpla rotatie. Trecereade la sistemul de coordonate local la sistemul global este esentiala,fiindca pana la urma ecuatiile de comportament ale elementelortrebuie asamblate ıntr-un sistem care sa descrie raspunsul ıntregiistructuri. Pentru asamblarea corecta a ecuatiilor elementale, varia-bilele asociate nodurilor trebuie exprimate ıntr-un acelasi sistem decoordonate global.

In figura 1.4 este prezentat sistemul de coordonate propriu unuielement oarecare de tip bara. Axele sistemului sunt notate ξ, res-pectiv η. Originea lui este plasata prin conventie ın nodul i, iar axaξ este paralela cu axa longitudinala a barei. Orientarea acestui sis-tem de coordonate local ın raport cu sistemul de coordonate globaleste definita prin unghiul θ(e) format de axa longitudinala a bareicu axa globala x. Pentru exprimarea lui θ(e) se adopta conventiauzuala din trigonometrie: unghiul θ(e) este masurat ıntotdeauna ınsens antiorar.

Variabilele care intervin ın ecuatiile elementale sunt deplasarilenodurilor si fortele care actioneaza asupra nodurilor din partea ele-mentelor la care sunt conectate. In continuare vom folosi urmatorulset de notatii (vezi figurile 1.4 si 1.5):

• ui, vi: deplasarile nodului i ın directiile axelor globale x, res-pectiv y;

• u(e)′

i , v(e)′

i : deplasarile nodului i ın directiile axelor locale ξ,respectiv η ale elementului (e);

• U(e)i , V

(e)i : fortele care actioneaza asupra nodului i din partea

elementului (e) ın directiile axelor globale x, respectiv y;

• U(e)′

i , V(e)′

i : fortele care actioneaza asupra nodului i din parteaelementului (e) ın directiile axelor locale ξ, respectiv η.

Page 20: Metoda Elementelor Finite

14 Capitolul 1

x

θ(e)

ξ

i

y

η jℓ(e)

(e)

Figura 1.4: Sistemul de coordonate local asociat unui element detip bara

ξ

η

i j(e)

ℓ(e)

u(e)′

i , v(e)′

i

U(e)′

i , V(e)′

i

u(e)′

j , v(e)′

j

U(e)′

j , V(e)′

j

Figura 1.5: Variabilele nodale exprimate ıntr-un sistem de coordo-nate local

Page 21: Metoda Elementelor Finite

Notiuni de baza ale metodei elementelor finite 15

Principiul actiunilor reciproce impune respectarea egalitatii

U(e)′

i + U(e)′

j = 0. (1.1)

De asemenea, ipoteza articulatiilor lipsite de frecare (vezi §1.2.1)are drept consecinta

V(e)′

i = 0 , V(e)′

j = 0. (1.2)

Daca deformatiile elementului (e) sunt mici, legea lui Hooke sepoate scrie sub forma

U(e)′

j

A= E

u(e)′

j − u(e)′

i

ℓ(e), (1.3)

unde ℓ(e) este lungimea elementului (fig. 1.4 - 1.5).Relatiile (1.1) - (1.3) permit scrierea urmatorului set de ecuatii

ce reflecta comportarea mecanica a unui element oarecare (e):

AE

ℓ(e)u

(e)′

i + 0 · v(e)′

i − AE

ℓ(e)u

(e)′

j + 0 · v(e)′

j = U(e)′

i ,

0 · u(e)′

i + 0 · v(e)′

i + 0 · u(e)′

j + 0 · v(e)′

j = V(e)′

i ,

− AE

ℓ(e)u

(e)′

i + 0 · v(e)′

i +AE

ℓ(e)u

(e)′

j + 0 · v(e)′

j = U(e)′

j ,

0 · u(e)′

i + 0 · v(e)′

i + 0 · u(e)′

j + 0 · v(e)′

j = V(e)′

j .

(1.4)

MEF utilizeaza de obicei o exprimare matriceala a ecuatiilorelementale. Aceasta din urma asigura mai multa conciziune si ge-neralitate relatiilor, fiind totodata si mai usor de implementat ınprograme. Drept consecinta, vom proceda si noi la rescrierea siste-mului de ecuatii (1.4) ın limbaj matriceal:

AE

ℓ(e)

26666664 1 0 | −1 00 0 | 0 0

−− −− −|− −− −−−1 0 | 1 00 0 | 0 0

37777775 266666664 u(e)′

i

v(e)′

i

−−u

(e)′

j

v(e)′

j

377777775 =

266666664 U(e)′

i

V(e)′

i

−−U

(e)′

j

V(e)′

j

377777775 . (1.5)

Page 22: Metoda Elementelor Finite

16 Capitolul 1

In continuare vor fi folosite notatiile

[M ] =

1 00 0

, (1.6)n

d(e)′

i

o=hu

(e)′

i , v(e)′

i

iT,

nd

(e)′

j

o=hu

(e)′

j , v(e)′

j

iT, (1.7)n

F(e)′

i

o=hU

(e)′

i , V(e)′

i

iT,

nF

(e)′

j

o=hU

(e)′

j , V(e)′

j

iT. (1.8)n

d(e)′

i

osind

(e)′

j

osunt vectori-coloana ın care sunt grupate de-

plasarile celor doua noduri exprimate ın sistemul de coordonate

local. In mod corespunzator,nF

(e)′

i

osinF

(e)′

j

osunt vectori-coloana

ın care sunt grupate componentele fortelor care actioneaza asupracelor doua noduri din partea elementului (e), exprimate de aseme-nea ın sistemul de coordonate local.

Cu ajutorul notatiilor (1.6) - (1.8), ecuatia matriceala (1.5) sepoate scrie astfel:

AE

ℓ(e)

24 [M ] − [M ]

− [M ] [M ]

35 264 nd(e)′

i

ond

(e)′

j

o 375 =

264 nF (e)′

i

onF

(e)′

j

o 375 . (1.9)

(1.9) descrie comportarea mecanica a unui element oarecare (e)ın sistemul de coordonate local. Spuneam la ınceputul prezentuluiparagraf ca este esentiala transferarea unei asemenea ecuatii ma-triceale ın sistemul de coordonate global. Pentru realizarea acesteitransformari se aplica o rotatie de unghi −θ(e) sistemului de co-ordonate local, ın asa fel ıncat el sa se suprapuna peste sistemul decoordonate global (fig. 1.4). Se pot scrie atunci urmatoarele relatiiıntre componentele deplasarilor si fortelor exprimate ın cele douasisteme de coordonate:

Page 23: Metoda Elementelor Finite

Notiuni de baza ale metodei elementelor finite 17

ui = u(e)′

i cos θ(e) − v(e)′

i sin θ(e),

vi = u(e)′

i sin θ(e) + v(e)′

i cos θ(e),

uj = u(e)′

j cos θ(e) − v(e)′

j sin θ(e),

vj = u(e)′

j sin θ(e) + v(e)′

j cos θ(e),

(1.10)

respectiv

Ui = U(e)′

i cos θ(e) − V(e)′

i sin θ(e),

Vi = U(e)′

i sin θ(e) + V(e)′

i cos θ(e),

Uj = U(e)′

j cos θ(e) − V(e)′

j sin θ(e),

Vj = U(e)′

j sin θ(e) + V(e)′

j cos θ(e).

(1.11)

(1.10) si (1.11) sunt echivalente cu ecuatiile matriceale

di =T (e)

nd

(e)′

i

o, dj =

T (e)

nd

(e)′

j

o, (1.12)n

F(e)i

o=T (e)

nF

(e)′

i

o,

nF

(e)j

o=T (e)

nF

(e)′

j

o, (1.13)

ın care am folosit notatiileT (e)

=

24 cos θ(e) − sin θ(e)

sin θ(e) cos θ(e)

35 , (1.14)

di = [ui, vi]T , dj = [uj, vj]

T , (1.15)nF

(e)i

o=hU

(e)i , V

(e)i

iT,

nF

(e)j

o=hU

(e)j , V

(e)j

iT. (1.16)

T (e)

este asa-numita matrice de rotatie locala. di si dj suntvectori-coloana care grupeaza deplasarile celor doua noduri ex-primate ın sistemul de coordonate global. In mod corespunzator,

Page 24: Metoda Elementelor Finite

18 Capitolul 1nF

(e)i

osinF

(e)j

osunt vectori-coloana care grupeaza componentele

fortelor ce se exercita asupra celor doua noduri din partea elemen-tului (e), exprimate de asemenea ın sistemul de coordonate global.

MatriceaT (e)

este ortogonala, adica are proprietatea

T (e)T

T (e)=T (e)

T (e)

T= [I2] , (1.17)

unde [I2] reprezinta matricea unitate de ordinul al doilea. Relatiile(1.17) pot fi verificate prin calcul direct pornind de la expresia(1.14) a lui

T (e)

.

Transformarile (1.12) - (1.13) urmeaza a fi ınlocuite ın ecuatia(1.9). Dupa cateva operatii matematice foarte simple, ın care se uti-lizeaza proprietatea de ortogonalitate (1.17), se obtine urmatoarearelatie matriceala care descrie comportarea mecanica a unui elementoarecare (e) ın termenii variabilelor nodale raportate la sistemul decoordonate global:

AE

ℓ(e)

264 T (e)

[M ]

T (e)

T −T (e)

[M ]

T (e)

T−T (e)

[M ]

T (e)

T T (e)

[M ]

T (e)

T 375××24 didj

35 =

264 nF (e)i

onF

(e)j

o 375 .

(1.18)

Se observa simetria si simplitatea acestei formule. Pentru utilizarea

ei este necesara evaluarea produsuluiT (e)

[M ]

T (e)

T. Daca se tine

cont de relatiile (1.6) si (1.14), rezulta imediatT (e)

[M ]

T (e)

T=

24 cos2 θ(e) sin θ(e) · cos θ(e)

sin θ(e) · cos θ(e) sin2 θ(e)

35 . (1.19)

Ecuatiile (1.18) - (1.19) sunt valabile pentru oricare din celecinci elemente ale retelei care modeleaza cadrul articulat. Daca setine cont de informatiile din figura 1.3, respectiv tabelele 1.1 - 1.2,se pot scrie urmatoarele particularizari ale lui (1.18):

Page 25: Metoda Elementelor Finite

Notiuni de baza ale metodei elementelor finite 19

• elementul (1): i = 1, j = 2, ℓ(1) = ℓ, θ(1) = π/2

AE

26664 0 0 0 00 1 0 −10 0 0 00 −1 0 1

37775 26664 u1

v1

u2

v2

37775 =

2666664 U(1)1

V(1)1

U(1)2

V(1)2

3777775 ; (1.20)

• elementul (2): i = 1, j = 3, ℓ(2) = ℓ, θ(2) = 0

AE

26664 1 0 −1 00 0 0 0

−1 0 1 00 0 0 0

37775 26664 u1

v1

u3

v3

37775 =

2666664 U(2)1

V(2)1

U(2)3

V(2)3

3777775 ; (1.21)

• elementul (3): i = 2, j = 3, ℓ(3) = ℓ√

2, θ(3) = 7π/4

AE

2√

2 ℓ

26664 1 −1 −1 1−1 1 1 −1−1 1 1 −1

1 −1 −1 1

37775 26664 u2

v2

u3

v3

37775 =

2666664 U(3)2

V(3)2

U(3)3

V(3)3

3777775 ; (1.22)

• elementul (4): i = 2, j = 4, ℓ(4) = ℓ, θ(4) = 0

AE

26664 1 0 −1 00 0 0 0

−1 0 1 00 0 0 0

37775 26664 u2

v2

u4

v4

37775 =

2666664 U(4)2

V(4)2

U(4)4

V(4)4

3777775 ; (1.23)

• elementul (5): i = 3, j = 4, ℓ(5) = ℓ, θ(5) = π/2

AE

26664 0 0 0 00 1 0 −10 0 0 00 −1 0 1

37775 26664 u3

v3

u4

v4

37775 =

2666664 U(5)3

V(5)3

U(5)4

V(5)4

3777775 . (1.24)

Page 26: Metoda Elementelor Finite

20 Capitolul 1

1.2.3 Sistemul de ecuatii care descrie raspunsulmecanic al retelei de elemente finite

Ecuatiile elementale se combina ıntr-un sistem care modeleazaraspunsul retelei ca ansamblu. Fiecare articulatie (nod) este supusaactiunii fortelor din barele (elementele) pe care le conecteaza, pre-cum si actiunii eventualelor forte exterioare. Pentru ca structura safie ın echilibru, este necesar ca rezultantele fortelor din articulatiisa fie zero. Aceasta observatie conduce la urmatoarea regula simplade asamblare a ecuatiilor elementale: ele trebuie combinate ın asafel ıncat membrii drepti ai ecuatiilor din sistemul care rezulta sa fiesuma fortelor interne asociate unui nod. Daca nodul nu este supusactiunii nici unei ıncarcari exterioare, suma fortelor interne trebuiesa fie evident nula. In schimb, daca asupra nodului actioneaza oıncarcare exterioara, conditia de echilibru impune ca suma fortelorinterne asociate sa fie egala cu acea forta exterioara.

Pentru aplicarea strategiei de asamblare prezentate mai sus sefoloseste informatia de conectivitate a elementelor din tabelul 1.1.Acest tabel ne arata, de exemplu, ca forta interna asociata nodului1 rezulta prin ınsumarea contributiilor elementelor (1) si (2). La fel,forta interna din nodul 3 este suma contributiilor elementelor (2),(3) si (4) (vezi si figura 1.3).

Exploatand ecuatiile elementale (1.20) - (1.24) si tabelul 1.1,obtinem urmatorul sistem:

AE

2666666666666666641 0 0 0 −1 0 0 0

0 1 0 −1 0 0 0 0

0 0 1 + 12√

2− 1

2√

2− 1

2√

21

2√

2−1 0

0 −1 − 12√

21 + 1

2√

21

2√

2− 1

2√

20 0

−1 0 − 12√

21

2√

21 + 1

2√

2− 1

2√

20 0

0 0 12√

2− 1

2√

2− 1

2√

21 + 1

2√

20 −1

0 0 −1 0 0 0 1 0

0 0 0 0 0 −1 0 1

377777777777777775×

Page 27: Metoda Elementelor Finite

Notiuni de baza ale metodei elementelor finite 21

×

266666666666666666664u1

v1

u2

v2

u3

v3

u4

v4

377777777777777777775 =

266666666666666666666664U

(1)1 + U

(2)1

V(1)1 + V

(2)1

U(1)2 + U

(3)2 + U

(4)2

V(1)2 + V

(3)2 + V

(4)2

U(2)3 + U

(3)3 + U

(5)3

V(2)3 + V

(3)3 + V

(5)3

U(4)4 + U

(5)4

V(4)4 + V

(5)4

377777777777777777777775. (1.25)

Din figura 1.3 rezulta de asemenea urmatoarele ecuatii de echilibrupentru noduri:

U(1)1 + U

(2)1 = X1,

V(1)1 + V

(2)1 = Y1,

U(1)2 + U

(3)2 + U

(4)2 = 0,

V(1)2 + V

(3)2 + V

(4)2 = 0,

U(2)3 + U

(3)3 + U

(5)3 = 0,

V(2)3 + V

(3)3 + V

(5)3 = Y3,

U(4)4 + U

(5)4 = − P√

2,

V(4)4 + V

(5)4 = − P√

2,

(1.26)

ın care X1 si Y1 sunt reactiunile din reazemul fix asociat nodului 1,iar Y3 este reactiunea din reazemul glisant asociat nodului 3.

(1.25) - (1.26) permit scrierea sistemului de ecuatii global (numitasa pentru ca modeleaza raspunsul ıntregii retele):

Page 28: Metoda Elementelor Finite

22 Capitolul 1

AE

2666666666666666641 0 0 0 −1 0 0 0

0 1 0 −1 0 0 0 0

0 0 1 + 12√

2− 1

2√

2− 1

2√

21

2√

2−1 0

0 −1 − 12√

21 + 1

2√

21

2√

2− 1

2√

20 0

−1 0 − 12√

21

2√

21 + 1

2√

2− 1

2√

20 0

0 0 12√

2− 1

2√

2− 1

2√

21 + 1

2√

20 −1

0 0 −1 0 0 0 1 0

0 0 0 0 0 −1 0 1

377777777777777775×

×

266666666666666666664u1

v1

u2

v2

u3

v3

u4

v4

377777777777777777775 =

266666666666666666664X1

Y1

0

0

0

Y3

− P√2

− P√2

377777777777777777775 . (1.27)

Trebuie observat faptul ca deplasarile u1 si v1 ale nodului 1,precum si deplasarea v3 a nodului 3 sunt cunoscute. Intr-adevar,reazemele reprezentate ın figura 1.3 sunt echivalente cu urmatoareleconditii:

u1 = 0, v1 = 0, v3 = 0. (1.28)

In schimb, reactiunile asociate acestor grade de libertate nodalesunt necunoscute (X1, Y1 si Y3). Deci (1.27) este un sistem de optecuatii liniare avand ca necunoscute marimile X1, Y1, u2, v2, u3,Y3, u4 si v4, repartizate ın ambii membri ai ecuatiilor (deplasarilese afla ın membrul stang, iar fortele ın membrul drept).

In practica se adopta urmatoarea strategie de rezolvare a siste-mului (1.27):

Page 29: Metoda Elementelor Finite

Notiuni de baza ale metodei elementelor finite 23

• se rezolva mai ıntai subsistemul format din ecuatiile carecontin numai necunoscutele de tip deplasare (el este formatdin a treia, a patra, a cincea, a saptea si a opta ecuatie asistemului (1.27));

• ın etapa urmatoare, cand deja se cunosc deplasarile tuturornodurilor, pot fi utilizate si restul ecuatiilor (adica prima,a doua si a sasea ecuatie din (1.27)) pentru determinareareactiunilor X1, Y1 si Y3.

In continuare vom vedea cum se aplica practic aceasta strategie. Din(1.27), daca se tine cont de conditiile (1.28), se obtine urmatorulsubsistem de ecuatii, ale carui necunoscute sunt deplasarile u2, v2,u3, u4 si v4:

AE

266666666641 + 1

2√

2− 1

2√

2− 1

2√

2−1 0

− 12√

21 + 1

2√

21

2√

20 0

− 12√

21

2√

21 + 1

2√

20 0

−1 0 0 1 0

0 0 0 0 1

3777777777526666666664

u2

v2

u3

u4

v4

37777777775 =

=

266666666640

0

0

− P√2

− P√2

37777777775 . (1.29)

Programele de analiza cu elemente finite rezolva numeric sis-temele de ecuatii de tipul (1.29). In acest scop, ele utilizeaza deregula metoda eliminarii Gauss [3, 4, 5, 8, 9, 10]. Am putea si noisa procedam ın felul acesta, calculand ın prealabil valorile numericeale coeficientilor si termenilor liberi din (1.29). Totusi, consideramca un asemenea calcul este lipsit de valoare didactica. Prezentam

Page 30: Metoda Elementelor Finite

24 Capitolul 1

deci mai jos expresiile simbolice ale solutiei sistemului (1.29), ıncare vom ınlocui valorile marimilor A, ℓ, P si E pentru a determinadeplasarile nodale necunoscute:

u2 = −2 +

√2 Pℓ

AE, v2 = − 1√

2

Pℓ

AE,

u3 = − 1√2

Pℓ

AE,

u4 = −

2 +3√2

Pℓ

AE, v4 = − 1√

2

Pℓ

AE.

(1.30)

Daca se tine cont de valorile numerice ale marimilor A, ℓ, P si E,se obtine Pℓ/(AE) =

√2 × 10−5 m. Din (1.30) rezulta atunci

u2 = −2 + 2

√2× 10−5 = −4,83 × 10−5 m,

v2 = −10−5 m,

u3 = −10−5 m,

u4 = −3 + 2

√2× 10−5 = −5,83 × 10−5 m,

v4 = −10−5 m.

(1.31)

Folosind acum (1.28) si (1.31), din prima, a doua si a saseaecuatie a sistemului (1.27) obtinem reactiunile (vezi si (1.30))

X1 =AE

ℓ(u1 − u3) =

P√2,

Y1 =AE

ℓ(v1 − v2) =

P√2,

Y3 =AE

2√

2 ℓ

u2 − v2 − u3 +

1 + 2

√2v3 − 2

√2v4

= 0.

(1.32)

Inlocuind valoarea lui P ın (1.32), deducem

X1 = 1000 N, Y1 = 1000 N, Y3 = 0. (1.33)

Page 31: Metoda Elementelor Finite

Notiuni de baza ale metodei elementelor finite 25

1.2.4 Starea de eforturi proprie fiecarui elemental retelei

Dupa rezolvarea sistemului de ecuatii global dispunem de toatedatele necesare pentru evaluarea tensiunilor axiale din fiecare ele-ment (bara). In acest scop, vom folosi ecuatiile matriceale (1.20) -(1.24).

De exemplu, din (1.22) rezulta componentele fortei pe care oexercita elementul (3) asupra nodurilor 2, respectiv 3 (componen-tele ei sunt exprimate ın sistemul de coordonate global). De fapt,

este suficient sa calculam U(3)3 si V

(3)3 , fiindca U

(3)2 = −U

(3)3 , iar

V(3)2 = −V

(3)3 . Din (1.22) obtinem urmatoarele expresii ale lui U

(3)3

si V(3)3 :

U(3)3 =

AE

2√

2 ℓ(−u2 + v2 + u3 − v3) ,

V(3)3 =

AE

2√

2 ℓ(u2 − v2 − u3 + v3) .

(1.34)

Prin ınlocuirea deplasarilor nodale date de (1.30) ın (1.34), rezultaın continuare

U(3)3 =

P√2

= 1000 N, V(3)3 = − P√

2= −1000 N. (1.35)

Stim ca barele cadrului sunt solicitate doar pe directie longitudi-nala. Pentru a putea evalua tensiunea axiala corespunzatoare avemnevoie de componentele fortei exprimate ın sistemul de coordonatelocal al elementului (3). Din (1.13) - (1.14) rezulta24 U

(3)′

3

V(3)′

3

35 =

24 cos θ(3) − sin θ(3)

sin θ(3) cos θ(3)

35T 24 U(3)3

V(3)3

35 , (1.36)

unde θ(3) = 7π/4. (1.36) este echivalenta cu

U(3)′

3 =1√2

U

(3)3 − V

(3)3

, V

(3)′

3 =1√2

U

(3)3 + V

(3)3

. (1.37)

Page 32: Metoda Elementelor Finite

26 Capitolul 1

Prin ınlocuirea valorilor lui U(3)3 si V

(3)3 date de (1.35) ın (1.37)

obtinem

U(3)′

3 =2000√

2= 1414 N, V

(3)′

3 = 0. (1.38)

Daca se tine cont de faptul ca axa locala ξ este orientata de la nodul2 spre nodul 3, din (1.38) rezulta ca elementul (3) este solicitat laıntindere. Tensiunea axiala corespunzatoare este

σ(3) =U

(3)′

3

A=

1414

5 × 10−4Nm−2 = 2,83 × 106 Nm−2. (1.39)

Metodologia prezentata mai sus poate fi aplicata pentru fiecaredin celelalte patru elemente ale retelei ın vederea determinarii ten-siunilor axiale.

1.3 Concluzii

Exemplul prezentat ın paragrafele anterioare evidentiaza unelecaracteristici importante ale MEF.

Trebuie sa remarcam ın primul rand faptul ca MEF construiestemodelul numeric prin discretizare. Discretizarea nu ınseamna nu-mai divizarea unui corp continuu ın subdomenii. Ea are, de fapt,implicatii matematice mult mai profunde ın ceea ce priveste rezol-varea numerica a problemei studiate. In esenta, discretizarea presu-pune, ın cazul MEF, ınlocuirea distributiei reale a necunoscutelorproblemei cu o aproximare care depinde de valorile acestora ıntr-unnumar finit de puncte (nodurile retelei de elemente finite). Astfel,problema initiala se reduce la determinarea valorilor necunoscutelorıntr-un numar finit de puncte.

Am vazut ın paragrafele anterioare ca valorile necunoscutelordin noduri sunt solutiile unui sistem de ecuatii. Acest sistem seobtine impunand ca aproximarile necunoscutelor sa verifice un setde legi specifice problemei studiate. In exemplul prezentat a fostvorba de legea lui Hooke, prin care era descrisa comportarea elastica

Page 33: Metoda Elementelor Finite

Concluzii 27

a fiecarei bare din componenta cadrului articulat. In alte problemepot aparea legi fizice diferite. De exemplu, ın problemele de trans-fer termic este folosita ecuatia caldurii. Totusi, oricare ar fi tipulproblemei analizate, strategia de aplicare a MEF ramane aceeasi.

Daca legile care guverneaza procesul studiat sunt liniare, siste-mul de ecuatii care se obtine este si el liniar (am vazut acest lucru ınparagrafele anterioare). Daca ınsa legile de pornire sunt neliniare,sistemul de ecuatii va fi ın general neliniar. Pentru rezolvarea luinumerica vor trebui adoptate atunci metodologii diferite fata decea prezentata ın problema-exemplu de mai sus. Nu vom intra ındetalii, fiindca modelele neliniare nu fac obiectul acestui curs.

Sistemul de ecuatii obtinut ın §1.2.3 are o serie de particularitatiinteresante. In primul rand, trebuie sa remarcam faptul ca matriceacoeficientilor este simetrica (vezi (1.27)). In plus, se mai observa caun numar mare de coeficienti sunt nuli (se spune ca matricea esterara). Daca se adopta o strategie buna de numerotare a nodurilor,toti coeficientii nenuli ai sistemului de ecuatii se strang ıntr-o zonadin vecinatatea diagonalei principale a matricii coeficientilor (ınacest caz, se spune ca matricea este de tip banda).

Proprietatile matricii coeficientilor au drept consecinta o seriede avantaje de natura numerica. In primul rand, programele MEFnu vor fi obligate sa stocheze ın memoria calculatorului ıntreagamatrice a coeficientilor. Va fi suficient sa memoreze doar coeficientiinenuli aflati pe diagonala principala si deasupra acesteia. Rezultaastfel economii de memorie semnificative si o crestere a vitezei derezolvare a sistemului de ecuatii (prin eliminarea operatiilor asupracoeficientilor nuli nestocati).

In final, mai trebuie sa mentionam faptul ca modelul cu ele-mente finite construit ın paragrafele anterioare nu este decat o va-rianta a metodei deplasarilor [1]. In forma aceasta, MEF se poateaplica numai pentru analiza unor structuri de tipul cadrelor. Vomvedea ın capitolele urmatoare ca se pot construi modele cu elementefinite mult mai generale (de exemplu, de tip variational sau rezi-dual [4, 5, 8, 9, 10]). Utilizarea acestora nu mai este restransa la

Page 34: Metoda Elementelor Finite

28 Capitolul 1

o singura clasa de probleme. De fapt, ın forma sa reziduala, MEFeste o metoda numerica aplicabila oricarei probleme care poate fiexprimata ın limbaj matematic sub forma unui sistem de ecuatiidiferentiale sau cu derivate partiale.

Page 35: Metoda Elementelor Finite

Capitolul 2

Caracteristici distinctive alemetodei elementelor finite

Din punct de vedere matematic, multe probleme de calcul in-gineresc se prezinta sub forma unor sisteme de ecuatii cu derivatepartiale (evident, ın cazurile cele mai simple, ele se pot reduce laecuatii diferentiale). Metoda elementelor finite (MEF) este unuldintre procedeele numerice pe care le putem folosi pentru rezolva-rea aproximativa a acestor probleme.

MEF construieste aproximarea solutiei exacte ın felul urmator:

• domeniul spatial pe care ıl ocupa sistemul supus analizei estedivizat ıntr-un numar finit de subdomenii reciproc disjunctenumite elemente;

• ın interiorul fiecarui element se defineste o aproximare de tip

polinomial a necunoscutei (necunoscutelor) problemei;

• fiecare asemenea polinom depinde de valorile necunoscutei(necunoscutelor) ıntr-un numar finit de puncte asociate ele-mentului si numite noduri ;

• aproximantele elementale sunt asamblate ımpreuna pentru aforma o aproximanta globala a necunoscutei (necunoscutelor).

29

Page 36: Metoda Elementelor Finite

30 Capitolul 2

Aproximanta globala trebuie ınteleasa ca o functie care depindede un numar finit de parametri deocamdata nedeterminati. Acestiparametri sunt valorile necunoscutei (necunoscutelor) ın nodurileretelei. Atribuind valori arbitrare parametrilor, obtinem o multimeinfinita de aproximante globale virtual acceptabile pentru problemaanalizata. Aceasta multime este cunoscuta ın literatura de specia-litate sub denumirea de familie a functiilor de test [4]. MEF cautaprintre functiile de test cea mai buna aproximare globala a solutieiexacte. Propriu-zis, solutia numerica obtinuta prin aplicarea MEFse apropie cat poate ea de bine de solutia exacta pe ansambluldomeniului spatial discretizat (nu numai ın anumite puncte).

In acest scop, MEF foloseste un criteriu global de minimizarea erorii solutiei numerice. Sistemul de ecuatii cu derivate partialenu poate fi utilizat direct pentru obtinerea unui astfel de crite-riu de optim global. Aceasta din cauza ca operatorii diferentialidescriu doar comportarea locala a marimilor necunoscute. Crite-riul de optimizare folosit de MEF este construit sub forma uneiintegrale definite pe ıntregul domeniu ocupat de sistemul supusanalizei. Este cunoscut faptul ca, spre deosebire de derivate, inte-grala descrie o comportare a uneia sau mai multor marimi mediatala scara unui domeniu. Ea este deci un instrument matematic adec-vat pentru obtinerea unui criteriu de minimizare global. Vom vedeaın paragrafele urmatoare ca, ın cazul problemelor de elasticitate,acest criteriu are o semnificatie fizica foarte clara, exprimata printeorema de minim a energiei potentiale.

Impunand conditia ca solutia numerica sa corespunda minimuluierorii, MEF construieste un sistem de ecuatii algebrice. Rezolvandacest sistem, obtinem valorile necunoscutei (necunoscutelor) ın no-durile retelei de elemente finite. Inlocuind mai departe valorile de-terminate ın functiile de aproximare elementala, rezulta aproximariale variatiei fiecarei necunoscute la nivelul elementelor, deci pana laurma gasim o aproximare a solutiei problemei pe ıntregul domeniu.

Vom ıntelege mult mai bine acest algoritm daca ıl vom aplicape un exemplu simplu. Vom considera problema sagetilor unei corzi

Page 37: Metoda Elementelor Finite

Formularea problemei-exemplu. Modelul sau analitic 31

flexibile pretensionate. In limbaj matematic, aceasta problema seexprima sub forma unei ecuatii diferentiale care se poate integraexact.

Dupa ce vom obtine solutia analitica, vom aplica MEF pentrurezolvarea numerica a problemei. Vom vedea pe rand toate etapelede construire si rezolvare a unui model cu elemente finite.

Am ales intentionat o problema pentru care putem obtine sisolutia analitica, pentru ca astfel vom avea posibilitatea sa evaluamsi precizia solutiei numerice. Observatiile de pe parcurs vor furnizacateva reguli practice ce se vor dovedi utile ın contextul rezolvariiunor probleme mai complexe.

2.1 Formularea problemei-exemplu.

Modelul sau analitic

Ne propunem sa determinam sagetile unei corzi flexibile pre-tensionate de o forta P = 1000 N (fig. 2.1) si supusa concomitentactiunii propriei greutati. Coarda are urmatoarele caracteristici:

• lungimea ℓ = 60 m;

• aria sectiunii transversale A = 10−5 m2;

• greutatea unitatii de lungime q = 0,8 Nm−1;

• materialul corzii are un comportament elastic.

Vom admite o flexibilitate perfecta a corzii. Cu alte cuvinte, con-sideram ca ın coarda apar numai solicitari de ıntindere (nu siıncovoieri sau forfecari). De asemenea, vom neglija micile variatiiale tensiunilor din coarda cauzate de greutatea proprie. Aceastaınseamna ca tensiunea va fi presupusa practic constanta pe toatalungimea corzii:

σ ≈ P

A= const. (2.1)

Page 38: Metoda Elementelor Finite

32 Capitolul 2

q

B

xP

ℓ = 60 m

A = 10−5 m2

q = 0,8 Nm−1

P = 1000 N

y

A

Figura 2.1: Coarda flexibila pretensionata si supusa actiunii proprieigreutati

Pentru a descrie deformatia corzii vom folosi un sistem de co-ordonate carteziene ale carui axe x si y sunt orientate asa cum searata ın figura 2.1. Vom admite ca deformatiile sunt mici. Se potaplica atunci urmatoarele ipoteze simplificatoare:

• sectiunile transversale plane ale corzii se deplaseaza practicpe verticala;

• distorsiunile sectiunilor transversale sunt neglijabile.

In aceste conditii, deplasarile verticale ale particulelor materialeaflate initial pe axa x pot fi considerate reprezentative pentru toateparticulele din sectiunile transversale corespondente. Altfel spus,deplasarea verticala v este functie numai de coordonata x (nu si decoordonata y):

v = v (x) , x ∈ [0, ℓ] . (2.2)

Dintre componentele tensorului deformatiei, singura care esteneneglijabila este deformatia longitudinala. Pentru a gasi expresiaei folosim informatia din figura 2.2.

In aceasta figura am reprezentat un element liniar infinitezimalde coarda ın configuratia initiala, respectiv ın stare deformata. Lun-gimea elementului de coarda deformata poate fi exprimata astfel:

Page 39: Metoda Elementelor Finite

Formularea problemei-exemplu. Modelul sau analitic 33

x

ds

σA ≈ P

σA ≈ P

v

v+

dv

q

x dx

y

Figura 2.2: Element de coarda deformata sub actiunea proprieigreutati

ds =

Ìv +

dv

dxdx − v

2

+ dx2 =

Ì1 +

dv

dx

2

dx. (2.3)

Deformatia longitudinala a corzii va fi atunci

ε =ds − dx

dx=

Ì1 +

dv

dx

2

− 1. (2.4)

Ipoteza micilor deformatii este echivalenta cu

dvdx

2 ≪ 1. In acesteconditii, se poate folosi urmatoarea aproximare a lui ε obtinuta printrunchierea dezvoltarii ın serie Taylor a membrului drept din (2.4):

ε ≈ 1

2

dv

dx

2

. (2.5)

Formularea matematica a problemei se obtine impunand con-ditia de echilibru mecanic pentru elementul liniar de coarda dinfigura 2.2. Mai exact, suma proiectiilor verticale ale fortelor careactioneaza asupra elementului de lungime initiala dx trebuie sa fienula:

Page 40: Metoda Elementelor Finite

34 Capitolul 2

P

dv

dx+

d2v

dx2dx

!− P

dv

dx+ q dx = 0. (2.6)

(2.6) este echivalenta cu urmatoarea ecuatie diferentiala:

Pd2v

dx2+ q = 0, x ∈ [0, ℓ] . (2.7)

(2.7) exprima deci echilibrul mecanic al corzii pe directia axei y.Acestei ecuatii diferentiale trebuie sa i se ataseze conditiile la limita

v (0) = 0, v (ℓ) = 0, (2.8)

care exprima efectul reazemelor din punctele A si B (fig. 2.1).Rezolvarea analitica a ecuatiei (2.7) cu conditiile (2.8) nu ridica

nici o dificultate. Prin doua integrari succesive ale lui (2.7) se obtineurmatoarea solutie generala:

v (x) = − q

2Px2 + c1x + c2, x ∈ [0, ℓ] , (2.9)

unde c1 si c2 sunt constante reale care se determina impunandrespectarea conditiilor (2.8). Propriu-zis, din (2.8) si (2.9) rezultaurmatorul sistem de ecuatii ın necunoscutele c1 si c2:

c2 = 0

− qℓ2

2P+ c1ℓ + c2 = 0.

(2.10)

Solutia lui este

c1 =qℓ

2P, c2 = 0. (2.11)

Din (2.9) si (2.11) se obtine atunci urmatoarea distributie a de-plasarilor verticale (sagetilor):

v (x) =q

2Px (ℓ − x) , x ∈ [0, ℓ] . (2.12)

(2.12) este ecuatia unei parabole. Ea reprezinta configuratia de echi-libru a corzii deformate sub actiunea propriei greutati.

Page 41: Metoda Elementelor Finite

Modelul cu elemente finite atasat problemei-exemplu 35

2.2 Modelul cu elemente finite atasat

problemei-exemplu

O prima etapa care trebuie parcursa la elaborarea unui modelcu elemente finite este identificarea clara a necunoscutei (necunos-cutelor) problemei analizate. In cazul concret al corzii din figura2.1 avem de-a face cu o singura necunoscuta – deplasarea verti-cala v (sageata). In problemele complexe pot aparea mai multenecunoscute care, ın plus, pot fi ierarhizate. Vom avea atunci ne-

cunoscute primare, respectiv necunoscute derivate. De exemplu, ıncazul deformarii elastice a unui corp tridimensional (batiul uneimasini-unelte, o biela etc.), necunoscutele primare sunt componen-tele u, v si w ale campului vectorial al deplasarilor particulelorfata de pozitia lor initiala. Necunoscutele derivate sunt legate decele primare prin relatii ın care intervin operatori diferentiali (deaici si denumirea lor). In cazul problemelor de elasticitate pentrucorpuri tridimensionale, necunoscutele derivate sunt componenteleεx, εy, εz, γxy, γyz si γzx ale tensorului deformatiei, respectiv com-ponentele σx, σy, σz, τxy, τyz si τzx ale tensorului tensiunii. Ierar-hia necunoscutelor este definita de ordinea determinarii lor. Astfel,ın cazul unei probleme de elasticitate tridimensionala, ierarhia ne-cunoscutelor este urmatoarea:

• componentele u, v si w ale campului deplasarilor;

• componentele εx, εy, εz, γxy, γyz si γzx ale tensorului de-formatiei (ele rezulta din u, v si w folosind relatiile diferentialedintre deformatii si deplasari);

• componentele σx, σy, σz, τxy, τyz si τzx ale tensorului tensiunii(ele sunt legate de componentele tensorului deformatiei prinlegea lui Hooke).

In orice situatie, pe locul ıntai al ierarhiei se vor afla necunoscuteleprimare.

Page 42: Metoda Elementelor Finite

36 Capitolul 2

Pentru a reduce cat mai mult efortul de calcul implicat de re-zolvarea numerica a modelului cu elemente finite, acesta va trebuisa includa un set de necunoscute cat mai restrans. In acest scop, seva urmari exprimarea modelului exclusiv ın termenii necunoscutei(necunoscutelor) primare.

In cazul corzii din figura 2.1, aplicarea regulii de mai sus esteredundanta pentru ca problema nu are decat o necunoscuta (depla-sarea verticala v) care joaca automat rolul de necunoscuta primara.Pentru problemele mai complicate, aceasta regula poate da nasterela anumite dificultati ın etapa de construire a modelului cu elementefinite.

Fiindca necunoscuta v depinde de o singura variabila (coor-donata x), domeniul pe care va trebui sa ıl discretizam este uni-dimensional (intervalul [0, ℓ]). Dupa cum am mentionat ın parteaintroductiva a acestui capitol, elementele finite sunt portiuni reci-proc disjuncte ale domeniului de analiza. In cazul nostru, ele vor fisubintervale ale lui [0, ℓ], discretizarea neridicand nici o dificultate.Este important sa mentionam faptul ca ın practica avem de-a facede cele mai multe ori cu situatii mai complicate. De exemplu, ıncazul problemelor de elasticitate bidimensionala sau tridimensio-nala, domeniul de analiza va fi plan, respectiv spatial. Elementelefinite folosite pentru discretizarea lui vor fi regiuni ale planului sauale spatiului (vom vedea ın capitolul urmator cateva asemenea ele-mente). Operatia de discretizare nu va mai fi atunci chiar atat desimpla, ea facand de multe ori apel la experienta si simtul inginerescal analistului.

In interiorul fiecarui element se defineste o aproximare de tippolinomial a necunoscutei (necunoscutelor) din model. Pentru re-zolvarea numerica a problemei enuntate ın §2.1 vom folosi cel maisimplu element unidimensional si anume, elementul de ordinul ıntai(fig. 2.3). Denumirea lui este legata de gradul polinomului de apro-ximare utilizat:

v(e)(x) = a(e)0 + a

(e)1 x. (2.13)

Page 43: Metoda Elementelor Finite

Modelul cu elemente finite atasat problemei-exemplu 37

xj

v(e)(x) = a(e)0 + a

(e)1 x

v i

v j

(e)

xi ℓ(e)

xj

y

i

Figura 2.3: Elementul unidimensional de ordinul ıntai folosit pentrumodelarea corzii

Elementul din figura 2.3 este reprezentat ın sistemul de coordo-nate global. Evident, am fi putut folosi si un sistem de coordonatelocal (vezi problema rezolvata ın capitolul 1), ınsa, ın afara ge-neralitatii, nu am fi obtinut nici un avantaj, coarda flexibila fiindmodelata ca un singur interval de pe aceeasi axa. In cele ce urmeazavom renunta deci la utilizarea unui sistem de coordonate local.

Coeficientii a(e)0 si a

(e)1 din (2.13) se determina impunand

conditia ca polinomul v(e)(x) sa ia valorile necunoscutei v ın no-

durile elementului. Pentru calculul coeficientilor a(e)0 si a

(e)1 avem

nevoie de doua conditii, deci elementul va trebui sa aiba tot atateanoduri. In figura 2.3, nodurile au fost etichetate cu simbolurile i sij, iar coordonatele lor au fost notate xi, respectiv xj (lungimea ele-mentului este ℓ(e)). Valorile necunoscutei v ın noduri au fost notatecorespunzator cu vi, respectiv vj.

Page 44: Metoda Elementelor Finite

38 Capitolul 2

Impunand conditia mentionata mai sus, din (2.13) obtinem un

sistem de doua ecuatii ın necunoscutele a(e)0 si a

(e)1 :

a(e)0 + a

(e)1 xi = vi,

a(e)0 + a

(e)1 xj = vj.

(2.14)

Solutia lui este

a(e)0 =

xj

ℓ(e)vi −

xi

ℓ(e)vj , a

(e)1 = − 1

ℓ(e)vi +

1

ℓ(e)vj , (2.15)

unde (vezi figura 2.3)

ℓ(e) = xj − xi. (2.16)

Inlocuind (2.15) ın (2.13), obtinem expresia polinomului care apro-ximeaza sageata v la nivelul unui element finit oarecare (e):

v(e)(x) = viN(e)i (x) + vjN

(e)j (x), (2.17)

unde

N(e)i (x) =

xj − x

ℓ(e), N

(e)j (x) =

x − xi

ℓ(e)(2.18)

sunt functiile de forma asociate nodurilor elementale [4, 5, 8, 9, 10].Aceasta denumire se justifica prin faptul ca respectivele functii nudepind decat de caracteristicile geometrice ale elementului (ın cazulde fata, ele depind numai de coordonatele nodurilor si de lungimeaelementului).

Pentru discretizarea intervalului [0, ℓ] care modeleaza coarda dinfigura 2.1 vom folosi patru elemente unidimensionale de ordinulıntai. Am optat pentru un numar atat de redus de elemente finitenumai pentru ca vom rezolva manual problema. In practica se uti-lizeaza de cele mai multe ori un numar considerabil mai mare deelemente. Singurele restrictii asupra gabaritului retelei vor fi intro-duse atunci numai de performantele calculatorului si ale programu-lui MEF disponibil.

Page 45: Metoda Elementelor Finite

Modelul cu elemente finite atasat problemei-exemplu 39

x

y

1 2 3 4 5

(1) (2) (3) (4)v1 = 0 v5 = 0

v2 v3v4

v(1)(x) v(2)(x) v(3)(x) v(4)(x)

Figura 2.4: Reteaua de elemente finite care modeleaza coarda

Tabelul 2.1: Matricea de conexiuni a retelei din figura 2.4

Element Nod Nod(e) i j

(1) 1 2(2) 2 3(3) 3 4(4) 4 5

Tabelul 2.2: Coordonatele nodurilor din figura 2.4

Nod xk

k [m]

1 02 203 304 455 60

Page 46: Metoda Elementelor Finite

40 Capitolul 2

Reteaua de elemente finite pe care o vom folosi ın continuareeste prezentata ın figura 2.4. Notatiile respecta conventiile pe carele-am adoptat ın capitolul 1. Matricea de conexiuni a retelei dinfigura este prezentata ın tabelul 2.1, iar coordonatele celor cincinoduri sunt date ın tabelul 2.2.

Se observa ca elementele nu sunt egale ca lungime (nodurile dinfigura 2.4 nu sunt plasate echidistant). De fapt, din acest punct devedere, MEF este cat se poate de flexibila. Ea permite utilizatoru-lui sa stabileasca densitatea retelei si repartitia nodurilor ın functiede necesitati. Unul dintre parametrii prin care avem control asupracalitatii solutiei numerice este densitatea discretizarii. Propriu-zis,cu cat reteaua de elemente finite este mai densa, cu atat solutiaobtinuta este mai precisa. Din aceasta cauza, ın practica se utili-zeaza urmatoarea strategie de discretizare:

• ın zonele unde aproximantele polinomiale trebuie sa urma-reasca o variatie foarte accentuata a necunoscutei (necunos-cutelor), reteaua va fi mai densa;

• ın zonele unde variatia necunoscutei (necunoscutelor) este mailina, reteaua poate fi mai rara.

Folosirea unor retele cu densitate diferentiata permite obtinereaunor solutii numerice bune, ın conditiile unui efort de calcul minim.

In legatura cu repartitia nodurilor din figura 2.4, mai trebuiesa facem o observatie. Fiindca ne intereseaza sageata maxima acorzii si ne asteptam ca aceasta sa apara ın punctul de coordonatax = ℓ/2, am plasat acolo un nod. Avem astfel garantia ca valoareadeterminata va fi cea mai precisa posibil pentru tipul de elementfinit si reteaua din figura.

Informatia din tabelele 2.1 si 2.2 ne permite sa particularizamexpresia generala (2.17) a aproximantei v(e) pentru fiecare din celepatru elemente. Vom obtine astfel patru polinoame de gradul ıntai,v(1), v(2), v(3) si v(4), fiecare definit ın interiorul cate unui element.v(1) contine ca parametri nedeterminati valorile v1 si v2 ale sagetii.

Page 47: Metoda Elementelor Finite

Rezolvarea numerica a modelului cu elemente finite 41

In mod asemanator, v(2) contine v2 si v3, v(3) contine v3 si v4, iar v(4)

contine v4 si v5. Cele patru polinoame sunt componente ale unei sin-gure aproximante globale v care contine ca parametri nedeterminativalorile v1, v2, v3, v4 si v5 ale necunoscutei ın nodurile retelei. Dandvalori arbitrare acestor parametri, obtinem functii dintr-o aceeasimultime (asa-numita familie a functiilor de test despre care am dis-cutat la ınceputul capitolului). Vom vedea ın paragraful urmatorcum reuseste MEF sa aleaga acea functie de test care asigura oaproximare globala optima a solutiei.

2.3 Rezolvarea numerica a modelului

cu elemente finite

Spuneam ın partea introductiva a acestui capitol ca MEF cautaacea functie de test care se apropie cel mai bine de solutia exactape ıntregul domeniu de analiza. In acest scop, ea impune conditiaca respectiva functie sa minimizeze o cantitate scalara definita subforma unei integrale.

In cazul problemelor de mecanica structurilor, procedura de maisus ısi are justificarea ın teorema de minim a energiei potentialeelastice. Se poate demonstra [10] ca starea de tensiuni si deformatiicare verifica ecuatiile de echilibru ale unui corp elastic realizeaza mi-nimul asa-numitei energii potentiale Π. In general, aceasta marimescalara are urmatoarea expresie:

Π = Πdef − Wext, (2.19)

ın care Πdef este energia de deformatie, iar Wext este lucrul mecanic

al ıncarcarilor exterioare.In cazul particular al corzii din figura 2.1, cei doi termeni din

membrul drept al lui (2.19) se pot scrie astfel (vezi si (2.1)):

Πdef =Z ℓ

0Pε dx, (2.20)

Page 48: Metoda Elementelor Finite

42 Capitolul 2

Wext =Z ℓ

0qv dx. (2.21)

Cu ajutorul lui (2.5), relatia (2.20) poate fi adusa la forma

Πdef =Z ℓ

0

P

2

dv

dx

2

dx. (2.22)

Din (2.19) si (2.21) - (2.22) rezulta atunci urmatoarea expresie aenergiei potentiale elastice a corzii:

Π =Z ℓ

0

"P

2

dv

dx

2

− qv

#dx. (2.23)

Π trebuie ınteleasa ca o aplicatie de un tip mai special. Ea esteo functie cu valori reale definita pe o familie de functii v. Dacaprivim relatia (2.23), constatam ca, atunci cand ınlocuim ın expre-sia lui Π diferite functii v, obtinem diferite valori reale. Teoremaenuntata mai sus ne asigura ca functia v : [0, ℓ] → R care verificaecuatia diferentiala (2.7) si conditiile la limita (2.8) realizeaza mi-nimul energiei potentiale elastice Π definite de (2.23).

Aceasta proprietate ne ofera criteriul global de minimizare aerorii de care avem nevoie la aplicarea MEF. Ea ne spune ca putemgasi solutia ecuatiei diferentiale care descrie problema corzii flexibilecautand acea functie care minimizeaza cantitatea Π(v). Adoptandaceasta cale alternativa, vom ınlocui ın expresia (2.23) a lui Π(v)aproximanta v construita ın §2.2. Evident, vom obtine astfel o apro-ximare a lui Π(v) care nu mai depinde decat de nedeterminatele luiv, adica de valorile nodale v1, v2, v3, v4 si v5 ale sagetii:

Π(v) → Π (v1, v2, v3, v4, v5) . (2.24)

Acesteia ıi vom impune conditia de minim. De fapt, vom cautavalorile nodale v1, v2, v3, v4 si v5 care realizeaza minimul lui Π. Inacest scop, vom folosi teorema lui Lagrange care ne da o conditienecesara de extrem pentru functii de mai multe variabile:

∂ Π

∂ vk

= 0, k = 1, . . . , 5. (2.25)

Page 49: Metoda Elementelor Finite

Rezolvarea numerica a modelului cu elemente finite 43

(2.25) nu este altceva decat un sistem de ecuatii ın necunoscutelev1, v2, v3, v4 si v5. Prin rezolvarea lui numerica obtinem acele valorinodale ale sagetii care minimizeaza cantitatea Π. Rezulta astfel siaproximanta globala v care minimizeaza Π. Fiindca Π este o aproxi-mare a lui Π, v tocmai determinat poate fi considerat o aproximareglobala optima a solutiei exacte v.

Aceasta ar fi, ın cateva cuvinte, descrierea traseului pe care ılvom parcurge ın continuare.

Aplicand proprietatea de aditivitate a integralelor ın raportcu domeniul, energia potentiala elastica poate fi scrisa ca sumacontributiilor celor patru elemente:

Π(v) =4X

e=1

Π(e)(v), (2.26)

unde

Π(e)(v) =Z

ℓ(e)

"P

2

dv

dx

2

− qv

#dx. (2.27)

Daca introducem ın expresia (2.27) a lui Π(e)(v) aproximanta (2.17)a sagetii, obtinem o aproximare a lui Π valabila la nivelul ele-mentului (e):

Π(e)(v) → Π(e) (vi, vj) =Z

ℓ(e)

8<:P

2

"dv(e)

dx

#2

− qv(e)

9=;dx, (2.28)

ın care vi si vj sunt valorile necunoscute ale sagetii ın cele douanoduri asociate.

Pentru evaluarea lui Π(e) (vi, vj) avem nevoie de derivata lui v(e)

ın raport cu x. Din (2.17) si (2.18) rezulta printr-un calcul foartesimplu

dv(e)

dx=

dN(e)i

dxvi +

dN(e)j

dxvj, (2.29)

undedN

(e)i

dx= − 1

ℓ(e),

dN(e)j

dx=

1

ℓ(e)· (2.30)

Page 50: Metoda Elementelor Finite

44 Capitolul 2

Cu ajutorul lui (2.17) si (2.29) - (2.30), expresia (2.28) a lui Π(e)

poate fi adusa la forma

Π(e) (vi, vj) =Z xj

xi

8<:P

2

vi − vj

ℓ(e)

2−

− qhN

(e)i vi + N

(e)j vj

i9=;dx.

(2.31)

Descompunand integrala din membrul drept al lui (2.31) ın termeni,obtinem

Π(e) (vi, vj) =P

2

vi − vj

l(e)

2 Z xj

xi

dx−

− qZ xj

xi

N(e)i dx

vi − q

Z xj

xi

N(e)j dx

vj.

(2.32)

Daca tinem cont de (2.16) si (2.18), putem scrie urmatoarele expre-sii ale integralelor din (2.32):Z xj

xi

dx = ℓ(e),Z xj

xi

N(e)i dx =

Z xj

xi

N(e)j dx =

ℓ(e)

2.

(2.33)

Din (2.32) - (2.33) rezulta atunci

Π(e) (vi, vj) =P

2ℓ(e)(vi − vj)

2 − qℓ(e)

2(vi + vj) . (2.34)

Expresia (2.34) este valabila pentru un element oarecare (e). Folo-sind informatia din tabelul 2.1, putem obtine particularizari ale eipentru fiecare din cele patru elemente ale retelei:

Page 51: Metoda Elementelor Finite

Rezolvarea numerica a modelului cu elemente finite 45

Π(1) (v1, v2) =P

2ℓ(1)(v1 − v2)

2 − qℓ(1)

2(v1 + v2) ,

Π(2) (v2, v3) =P

2ℓ(2)(v2 − v3)

2 − qℓ(2)

2(v2 + v3) ,

Π(3) (v3, v4) =P

2ℓ(3)(v3 − v4)

2 − qℓ(3)

2(v3 + v4) ,

Π(4) (v4, v5) =P

2ℓ(4)(v4 − v5)

2 − qℓ(4)

2(v4 + v5) .

(2.35)

Aproximarea lui Π la nivelul ıntregului domeniu de analiza re-zulta ca suma a celor patru cantitati definite de relatiile (2.35):

Π (v1, v2, v3, v4, v5) = Π(1) (v1, v2) + Π(2) (v2, v3) +

+ Π(3) (v3, v4) + Π(4) (v4, v5) .(2.36)

Conditiile necesare de minim (2.25) vor avea atunci urmatoareleexpresii:

∂ Π

∂ v1

=∂ Π(1)

∂ v1

= 0,

∂ Π

∂ v2

=∂ Π(1)

∂ v2

+∂ Π(2)

∂ v2

= 0,

∂ Π

∂ v3

=∂ Π(2)

∂ v3

+∂ Π(3)

∂ v3

= 0,

∂ Π

∂ v4

=∂ Π(3)

∂ v4

+∂ Π(4)

∂ v4

= 0,

∂ Π

∂ v5

=∂ Π(4)

∂ v5

= 0.

(2.37)

Ecuatiile (2.37) pot fi explicitate prin derivarea relatiilor (2.35):

Page 52: Metoda Elementelor Finite

46 Capitolul 2

∂ Π(1)

∂ v1

=P

ℓ(1)v1 −

P

ℓ(1)v2 −

qℓ(1)

2,

∂ Π(1)

∂ v2

= − P

ℓ(1)v1 +

P

ℓ(1)v2 −

qℓ(1)

2,

∂ Π(2)

∂ v2

=P

ℓ(2)v2 −

P

ℓ(2)v3 −

qℓ(2)

2,

∂ Π(2)

∂ v3

= − P

ℓ(2)v2 +

P

ℓ(2)v3 −

qℓ(2)

2,

∂ Π(3)

∂ v3

=P

ℓ(3)v3 −

P

ℓ(3)v4 −

qℓ(3)

2,

∂ Π(3)

∂ v4

= − P

ℓ(3)v3 +

P

ℓ(3)v4 −

qℓ(3)

2,

∂ Π(4)

∂ v4

=P

ℓ(4)v4 −

P

ℓ(4)v5 −

qℓ(4)

2,

∂ Π(4)

∂ v5

= − P

ℓ(4)v4 +

P

ℓ(4)v5 −

qℓ(4)

2.

(2.38)

Putem rescrie atunci sistemul (2.37) sub forma

P

26666666641

ℓ(1)− 1

ℓ(1)0 0 0

− 1ℓ(1)

1ℓ(1)

+ 1ℓ(2)

− 1ℓ(2)

0 0

0 − 1ℓ(2)

1ℓ(2)

+ 1ℓ(3)

− 1ℓ(3)

0

0 0 − 1ℓ(3)

1ℓ(3)

+ 1ℓ(4)

− 1ℓ(4)

0 0 0 − 1ℓ(4)

1ℓ(4)

3777777775××

2666666664 v1

v2

v3

v4

v5

3777777775 =q

2

2666666664 ℓ(1)

ℓ(1) + ℓ(2)

ℓ(2) + ℓ(3)

ℓ(3) + ℓ(4)

ℓ(4)

3777777775 .

(2.39)

Page 53: Metoda Elementelor Finite

Rezolvarea numerica a modelului cu elemente finite 47

Inainte de a trece la rezolvarea lui (2.39) trebuie sa discutamdespre impunerea conditiilor la limita (2.8). Pana ın acest moment,respectivele conditii nu s-au reflectat nicaieri ın modelul cu elementefinite. De fapt, daca am calcula determinantul matricii coeficientilorsistemului (2.39), am constata ca el este nul. Aceasta arata ca mo-delul nostru sufera de o nedeterminare matematica. Lipsa deter-minarii provine tocmai din neintroducerea conditiilor la limita (ıninterpretare mecanica, ea exprima absenta reazemelor de la cape-tele corzii, deci posibilitatea modelului de a fi translatat ın oricepozitie a spatiului).

Vom vedea ın continuare cum se introduc conditiile (2.8) ınmodelul cu elemente finite si cum afecteaza ele structura sistemuluide ecuatii (2.39).

Daca folosim informatia din figurile 2.1 si 2.4, constatam ca,ın modelul nostru cu elemente finite, conditiile la limita (2.8) suntechivalente cu

v1 = 0, v5 = 0. (2.40)

Relatiile (2.40) ne spun ca prima si a cincea ecuatie din sistemul(2.39) sunt redundante. Introducerea conditiilor la limita (2.8) ınmodelul cu elemente finite este deci echivalenta cu eliminarea aces-tor doua ecuatii din (2.39). Evident, si ın cele trei ecuatii ramase vatrebui sa operam anumite modificari si anume, sa transferam ter-menii corespunzatori lui v1 si v5 ın membrul drept (de fapt, aceastasituatie este naturala, fiindca, datorita conditiilor (2.40), respectiviitermeni reprezinta cantitati cunoscute). Obtinem astfel un sistemredus avand ca necunoscute numai sagetile din nodurile 2, 3 si 4:

P

2664 1ℓ(1)

+ 1ℓ(2)

− 1ℓ(2)

0

− 1ℓ(2)

1ℓ(2)

+ 1ℓ(3)

− 1ℓ(3)

0 − 1ℓ(3)

1ℓ(3)

+ 1ℓ(4)

3775 2664 v2

v3

v4

3775 =

=q

2

2664 ℓ(1) + ℓ(2)

ℓ(2) + ℓ(3)

ℓ(3) + ℓ(4)

3775 .

(2.41)

Page 54: Metoda Elementelor Finite

48 Capitolul 2

(2.41) nu mai sufera de nici o nedeterminare (putem verifica acestlucru calculand determinantul matricii coeficientilor). Solutia luieste

v2 =q

2Pℓ(1)

ℓ(2) + ℓ(3) + ℓ(4)

,

v3 =q

2P

ℓ(1) + ℓ(2)

ℓ(3) + ℓ(4)

,

v4 =q

2P

ℓ(1) + ℓ(2) + ℓ(3)

ℓ(4).

(2.42)

Folosind valorile numerice ale parametrilor P si q (vezi §2.1),obtinem q/(2P ) = 4×10−4 m−1. De asemenea, din tabelele 2.1 - 2.2rezulta lungimile elementelor: ℓ(1) = 20 m, ℓ(2) = 10 m, ℓ(3) = 15 msi ℓ(4) = 15 m. Relatiile (2.42) ne dau atunci urmatoarele valori alesagetii ın nodurile 2, 3 si 4:

v2 = 0,32 m, v3 = 0,36 m, v4 = 0,27 m. (2.43)

La acestea trebuie adaugate conditiile (2.40) care exprima sagetileın nodurile de capat.

Avand determinate valorile nodale ale necunoscutei v, dinrelatiile (2.17) - (2.18), ın care ınlocuim datele din tabelele 2.1 -2.2, obtinem imediat polinoamele care aproximeaza variatia sagetiila nivelul fiecarui element:

v(1) =x2 − x

ℓ(1)v1 +

x − x1

ℓ(1)v2 =

x

200,32, x ∈ [0, 20],

v(2) =x3 − x

ℓ(2)v2 +

x − x2

ℓ(2)v3 =

30 − x

100,32 +

x − 20

100,36,

x ∈ [20, 30],

v(3) =x4 − x

ℓ(3)v3 +

x − x3

ℓ(3)v4 =

45 − x

150,36 +

x − 30

150,27,

x ∈ [30, 45],

v(4) =x5 − x

ℓ(4)v4 +

x − x4

ℓ(4)v5 =

60 − x

150,27, x ∈ [45, 60].

(2.44)

Page 55: Metoda Elementelor Finite

Rezolvarea numerica a modelului cu elemente finite 49

v(1), v(2), v(3) si v(4) sunt componentele aproximantei globale v. Seobserva ın (2.44) ca fiecare componenta este definita pe un sub-interval care corespunde unuia din elementele retelei.

Pentru a evalua calitatea solutiei numerice obtinute prin rezol-varea modelului cu elemente finite, vom reprezenta pe acelasi graficpolinoamele v(1), v(2), v(3), v(4) si solutia analitica (2.12). Folosindvalorile numerice ale parametrilor ℓ, q si P (vezi §2.1), obtinemurmatoarea expresie a lui (2.12):

v(x) = 4 · 10−4x(60 − x), x ∈ [0, 60]. (2.45)

Un prim aspect pe care ıl vom verifica este precizia valorilornodale ale sagetii. Din (2.45) rezulta

v(0) = 0,

v(20) = 0,32 m,

v(30) = 0,36 m,

v(45) = 0,27 m,

v(60) = 0.

(2.46)

Comparand (2.46) cu (2.40) si (2.43), constatam ca sagetile nodalese suprapun perfect peste valorile date de solutia analitica (2.45):

v1 = v(0),

v2 = v(20),

v3 = v(30),

v4 = v(45),

v5 = v(60).

(2.47)

Totusi, trebuie sa notam faptul ca solutia numerica nu se mai supra-pune peste solutia exacta ın interiorul elementelor. Intr-adevar,daca privim relatiile (2.44) - (2.45), observam ca aproximanteleelementale ale sagetii sunt polinoame de gradul ıntai, iar solutia

Page 56: Metoda Elementelor Finite

50 Capitolul 2

y

1 2 3 4 5 x

Calcul numeric (MEF)

Solutie analitica

(1) (2) (3) (4)

Figura 2.5: Comparatie ıntre solutia numerica si solutia analitica aproblemei corzii flexibile

exacta este o distributie parabolica. Deci aproximanta globala veste formata dintr-o succesiune de segmente de dreapta ale carorcapete se afla pe curba ce reprezinta solutia analitica a problemei(fig. 2.5). Diagrama din figura 2.5 ne conduce la urmatoarea con-cluzie practica: solutia furnizata de MEF poate fi ımbunatatita dinpunct de vedere al preciziei prin utilizarea unei retele mai dense.Intr-adevar, daca am fi folosit mai mult de patru elemente pentrumodelarea corzii, solutia aproximativa v ar fi fost formata din maimulte segmente de dreapta si s-ar fi apropiat mai mult de distributiaanalitica a sagetii. Dar chiar si solutia obtinuta cu numai patruelemente finite de ordinul ıntai este satisfacatoare, fiindca ne davaloarea exacta a sagetii maxime si se apropie suficient de bine desolutia analitica ın restul domeniului.

2.4 Concluzii

Acest capitol a fost consacrat unei prezentari mai riguroase aprocedurii de aplicare a MEF. Am vazut ın detaliu care sunt eta-pele construirii si rezolvarii unui model cu elemente finite pentru

Page 57: Metoda Elementelor Finite

Concluzii 51

o problema de mecanica structurilor. In sinteza, aceste etape ar fiurmatoarele:

• alegerea tipului de element finit adecvat pentru problemadata;

• discretizarea domeniului de analiza (stabilirea numarului deelemente finite care vor fi utilizate si a pozitiei nodurilor aso-ciate; rezulta astfel matricea de conexiuni a retelei si tabelacoordonatelor nodale);

• construirea aproximantelor elementale ale necunoscutei (ne-cunoscutelor) problemei;

• ınlocuirea aproximantelor elementale ın expresia energiei po-tentiale;

• impunerea conditiei de minim asupra aproximarii cu elementefinite a energiei potentiale si asamblarea sistemului de ecuatiiglobal care rezulta prin aplicarea acestei conditii;

• introducerea conditiilor la limita (prin reducerea corespun-zatoare a sistemului de ecuatii global);

• rezolvarea sistemului de ecuatii global ın raport cu valorilenodale ale necunoscutei (necunoscutelor) problemei;

• reconstituirea aproximantelor elementale si asamblarea aces-tora ıntr-o aproximanta globala;

• interpretarea rezultatelor numerice.

In practica, atunci cand se foloseste un program MEF comer-cial, cele mai multe dintre etapele enumerate mai sus sunt par-curse automat, fara a mai fi necesara interventia analistului. Deregula, utilizatorul construieste un model geometric al domeniuluide analiza (folosind un program CAD sau chiar facilitatile de de-senare puse la dispozitie de programul MEF) si apoi alege tipul

Page 58: Metoda Elementelor Finite

52 Capitolul 2

de element finit. In acest stadiu, programul MEF realizeaza auto-mat discretizarea (tinand cont de anumite recomandari ale analis-tului ın ceea ce priveste densitatea locala a retelei). In continuare,utilizatorul introduce ın modelul discretizat ıncarcarile mecanicesi conditiile la limita de tipul restrictiilor de miscare. ProgramulMEF dispune acum de toate informatiile necesare pentru construi-rea sistemului de ecuatii global. Asamblarea si rezolvarea numericaa acestuia sunt etape complet automatizate. Toate programele MEFcomerciale folosite la ora actuala ofera si facilitati de reprezentaregrafica a solutiilor numerice sub forma a diferite tipuri de diagrame.In felul acesta, analistul este asistat si ın etapa de interpretare arezultatelor.

Procedura de construire a modelului cu elemente finite bazatape minimizarea energiei potentiale este general aplicabila pentrurezolvarea problemelor de elasticitate. Din punct de vedere ma-tematic, ea apartine asa-numitei familii a metodelor variationale

[4, 5, 8, 9, 10].

Page 59: Metoda Elementelor Finite

Capitolul 3

Prezentarea catorva tipuriuzuale de elemente finite

Discretizarea domeniului de analiza este prima etapa care tre-buie parcursa la rezolvarea numerica a unei probleme prin metodaelementelor finite (MEF). Discretizarea presupune luarea unor de-cizii ın legatura cu tipul, numarul si marimea elementelor utilizate.Ca ingineri suntem confruntati cu problema delicata a echilibruluidintre calitatea solutiei numerice si efortul de calcul necesar pen-tru obtinerea acesteia. In general, cu cat reteaua de elemente finiteeste mai densa, cu atat solutia este mai precisa. Totusi, o retea preafina conduce la sisteme de ecuatii foarte mari, a caror rezolvare estelaborioasa. Analistul trebuie sa faca apel la cunostintele teoreticesi la simtul sau ingineresc pentru a rafina reteaua numai ın zoneleunde sunt de asteptat gradienti mari ai necunoscutelor. In zoneleın care variatia acestora este mai putin accentuata, reteaua poatefi mai rara, fiindca si un numar redus de elemente poate conducela un rezultat acceptabil ca precizie.

Calitatea solutiei numerice este influentata si de caracteristicilematematice ale elementului finit utilizat. Concret, analistul trebuiesa aleaga un element ale carui proprietati (ın special gradul poli-nomului de aproximare) sa fie adecvate tipului de problema care

53

Page 60: Metoda Elementelor Finite

54 Capitolul 3

urmeaza a fi rezolvata. Importanta alegerii elementului finit nu tre-buie niciodata subestimata, fiindca o decizie incorecta ın aceastaetapa conduce adeseori la rezultate nesatisfacatoare sau eronate.

In practica, se constata ca majoritatea greselilor de discretizareau la origine urmatoarele cauze:

• insuficienta cunostintelor teoretice ın domeniul problemeiabordate;

• lipsa unor informatii referitoare la proprietatile matematiceale elementelor finite utilizate.

Chiar daca avem la dispozitie unul dintre cele mai performanteprograme de analiza cu elemente finite, nu suntem deloc degrevatide obligatia de a avea asemenea cunostinte. Datorita specificuluiacestui curs, ın cele ce urmeaza ne vom limita la prezentarea carac-teristicilor matematice ale catorva tipuri de elemente finite folositefrecvent la rezolvarea problemelor de inginerie mecanica. Spatiul li-mitat nu ne permite nicidecum o descriere a tuturor categoriilor deelemente finite care pot sa apara ın aplicatii. Cititorii interesati aula dispozitie o bogata bibliografie de specialitate [4, 5, 8, 9, 10]. Deasemenea, este foarte importanta consultarea manualelor de utili-zare ale programelor MEF.

3.1 Clasificare

O clasificare foarte generala (ınsa, din pacate, destul de vaga)ımparte elementele finite ın trei mari categorii, ın functie de dimen-sionalitatea domeniului la discretizarea caruia sunt folosite. Vomdistinge astfel:

• elemente finite unidimensionale;

• elemente finite bidimensionale;

• elemente finite tridimensionale.

Page 61: Metoda Elementelor Finite

Clasificare 55

Aceasta clasificare nu surprinde nici pe departe toate caracteristi-cile elementelor. De exemplu, ea nu face nici o referire la gradulpolinomului de aproximare. Din aceasta cauza, ın cadrul fiecareiadin categoriile mentionate mai sus pot fi realizate subclasificari ınfunctie de o serie de alte criterii. Daca se are ın vedere gradul poli-noamelor cu ajutorul carora sunt aproximate necunoscutele, pot fidistinse:

• elemente finite de ordinul ıntai (sau liniare);

• elemente finite de ordinul al doilea (sau patratice);

• elemente finite de ordinul al treilea (sau cubice);

• . . .

In practica sunt foarte rar utilizate elemente cu functii de aproxi-mare de grad superior lui trei.

Desi acest al doilea criteriu precizeaza o serie ıntreaga de as-pecte, nici el nu este exhaustiv. Propriu-zis, marea diversitate aelementelor ın uz la ora actuala face imposibila gasirea unui crite-riu de clasificare unic. Mai ales elementele finite de complexitatematematica ridicata sunt greu de inclus ıntr-o anumita categorie.

Pentru scopurile noastre, clasificarea de mai sus este totusisuficienta. Ne vom folosi de ea atunci cand vom prezenta unele din-tre elementele cel mai frecvent utilizate la rezolvarea problemelorde inginerie mecanica.

3.1.1 Elemente finite unidimensionale

Elementele finite unidimensionale sunt folosite atunci cand ma-rimea ce trebuie aproximata depinde de o singura variabila (veziproblema-exemplu rezolvata ın capitolul 2). Drept consecinta, ele-mentele unidimensionale sunt segmente de dreapta sau arce decurba de-a lungul carora ia valori variabila independenta a proble-mei (fig. 3.1). In capitolul 2 am folosit un astfel de element pentrudiscretizarea unei corzi perfect flexibile.

Page 62: Metoda Elementelor Finite

56 Capitolul 3

(a) (b) (c)

Figura 3.1: Elemente finite unidimensionale: (a) element de ordinulıntai (liniar); (b) elemente de ordinul al doilea (patratice); (c) ele-mente de ordinul al treilea (cubice)

Cel mai simplu element unidimensional este segmentul dedreapta prevazut cu doua noduri, cate unul la fiecare capat (fig.3.1.a). Un asemenea element are asociata o functie de aproximarede ordinul ıntai, motiv pentru care se si numeste liniar. De altfel,am si vazut cum se foloseste el ın exemplul din capitolul 2.

Mai complexe sunt elementele unidimensionale de ordinul al doi-lea (fig. 3.1.b), respectiv al treilea (fig. 3.1.c). In general, ele au nunumai noduri de capat, ci si noduri intermediare. Elementele deordinul al doilea (sau patratice) au trei noduri, iar functia lor deaproximare este de tip parabolic. Cele de ordinul al treilea au pa-tru noduri, iar functia lor de aproximare este o parabola cubica.Atat elemente unidimensionale de ordinul al doilea, cat si cele deordinul al treilea pot avea si variante curbilinii, asa cum se vede ınfigura 3.1.

3.1.2 Elemente finite bidimensionale

Elementele finite bidimensionale sunt utilizate atunci candmarimea ce trebuie aproximata depinde de doua variabile. De exem-plu, ele pot fi folosite pentru rezolvarea problemelor de elasticitatebidimensionale (probleme de stare plana de tensiuni, probleme destare plana de deformatii sau probleme cu simetrie axiala).

Page 63: Metoda Elementelor Finite

Clasificare 57

(a) (b) (c)

Figura 3.2: Elemente finite triunghiulare: (a) element de ordinulıntai (liniar); (b) elemente de ordinul al doilea (patratice); (c) ele-mente de ordinul al treilea (cubice)

Elementele finite bidimensionale pot fi subımpartite ın douamari clase, ın functie de forma lor:

• elemente triunghiulare (fig. 3.2);

• elemente patrulatere (fig. 3.3).

Elementele triunghiulare sunt, dupa cum arata chiar numele lor,regiuni de forma triunghiulara ale planului. Cel mai simplu estetriunghiul cu laturi rectilinii si trei noduri, cate unul ın fiecare varf(fig. 3.2.a). Functia de aproximare asociata acestui element este unpolinom de gradul ıntai ın doua variabile, motiv pentru care se sinumeste element triunghiular de ordinul ıntai (sau liniar).

Elemente triunghiulare mai complexe sunt cele de ordinul al doi-lea (fig. 3.2.b), respectiv cele de ordinul al treilea (fig. 3.2.c). Ele

Page 64: Metoda Elementelor Finite

58 Capitolul 3

(a) (b) (c)

Figura 3.3: Elemente finite patrulatere: (a) element biliniar; (b) ele-mente bipatratice; (c) elemente bicubice

au noduri nu numai ın varfuri, ci si pe laturi. Elementul de ordi-nul al treilea are chiar un nod interior. Elementele de ordinul aldoilea au ca functie de aproximare un polinom de gradul al doileacomplet, iar cele de ordinul al treilea au un polinom de gradul altreilea complet (evident, ın doua variabile). Se observa ın figura 3.2ca elementele triunghiulare de ordin superior lui unu au si variantecu laturi curbilinii. In practica, elementele triunghiulare de ordinulal treilea (sau mai mare) sunt rareori folosite, fiindca, din cauzanumarului prea mare de noduri, devin neeconomice din punct devedere al timpului de calcul. Aceasta remarca este valabila si pen-tru alte familii de elemente. In general, utilizatorii tind sa adopteelemente cu numar de noduri cat mai redus, dar care sa ofere ungrad de aproximare suficient pentru tipul de problema abordat.

Page 65: Metoda Elementelor Finite

Clasificare 59

Elementele patrulatere sunt regiuni cu patru laturi ale planului.Cel mai simplu este patrulaterul cu laturi rectilinii si patru noduri,cate unul ın fiecare varf (fig. 3.3.a). Functia de aproximare asociataeste de tip biliniar (vom reveni asupra acestui element ın §3.2.2, ast-fel ıncat semnificatia matematica a aproximantei va deveni clara).

In figura 3.3.b este prezentat un element patrulater cu opt no-duri, ımpreuna cu varianta sa curbilinie. Functia de aproximareasociata este de tip bipatratic. In sfarsit, ın figura 3.3.c apare unelement patrulater cu douasprezece noduri atat ın varianta recti-linie, cat si ın cea curbilinie. Functia lui de aproximare este de tipbicubic. Multi utilizatori considera ca deja elementul bipatratic areprea multe noduri. Din aceasta cauza, ın practica se constata opreferinta aproape exclusiva pentru elementul patrulater biliniar.

3.1.3 Elemente finite tridimensionale

Din marea clasa a elementelor tridimensionale, cel mai des fo-losite sunt urmatoarele doua familii:

• elementele finite tetraedrice (fig. 3.4);

• elementele finite hexaedrice (fig. 3.5).

In spatiul cu trei dimensiuni, numarul nodurilor elementelor tindesa creasca foarte repede odata cu cresterea gradului polinomului deaproximare. Din aceasta cauza, utilizatorii ısi limiteaza optiunea laelementele de ordin mai redus. Si noi ne vom rezuma la a descrieacele elemente tridimensionale care sunt mai des ıntalnite ın biblio-tecile programelor MEF comerciale.

In figura 3.4.a este prezentat cel mai simplu element tetraedric.El are fete plane si patru noduri, cate unul ın fiecare varf. Functia deaproximare asociata este un polinom de gradul ıntai ın trei variabile,motiv pentru care elementul se si numeste tetraedric de ordinul ıntai(sau liniar).

In figura 3.4.b este prezentat elementul tetraedric de ordinul aldoilea (sau patratic) atat ın varianta cu fete plane, cat si ın cea cu

Page 66: Metoda Elementelor Finite

60 Capitolul 3

(a) (b)

Figura 3.4: Elemente finite tetraedrice: (a) element de ordinul ıntai(liniar); (b) elemente de ordinul al doilea (patratice)

fete curbe. Se observa ca el are si noduri plasate pe muchii. Functiade aproximare asociata este un polinom complet de gradul al doileaın trei variabile.

Practic singurul element hexaedric utilizat ın aplicatii este celcu opt noduri si muchii rectilinii (fig. 3.5). Functia de aproximareasociata este de tip triliniar (vom reveni asupra acestui elementın §3.2.4, astfel ıncat semnificatia matematica a aproximantei vadeveni clara).

Exista si elemente hexaedrice de ordin superior (care au variantecurbilinii), ınsa utilizarea lor este foarte rara, fiindca au deja unnumar de noduri excesiv de mare.

Page 67: Metoda Elementelor Finite

Constructia aproximantelor elementale 61

Figura 3.5: Elementul finit hexaedric triliniar

3.2 Constructia aproximantelor pentru

cateva tipuri de elemente finite

Din cele discutate ın §3.1 ne dam seama ca un element finitnu este doar o regiune a spatiului, ci este o entitate mai complexareprezentata prin urmatoarele caracteristici:

• domeniul spatial ocupat;

• setul de noduri asociat;

• functia folosita pentru aproximarea necunoscutei (necunoscu-telor) problemei ın interiorul sau.

Spuneam ın capitolele precedente ca, de obicei, functia de aproxi-mare este de tip polinomial. Motivul pentru care se utilizeaza cu

Page 68: Metoda Elementelor Finite

62 Capitolul 3

preponderenta polinoame rezida ın usurinta cu care acestea pot fiderivate sau integrate.

Am vazut ca elementele finite nu trebuie sa aiba neaparat oforma regulata (de exemplu, elementele triunghiulare pot fi ne-echilatere, iar cele de ordin superior pot avea chiar laturi curbi-linii). Prin aceasta, MEF ısi dovedeste odata ın plus flexibilitatea,oferindu-ne mijloace eficiente pentru discretizarea unor domeniispatiale a caror forma poate fi oricat de complicata.

Totusi, construirea unor functii de aproximare de tip polinomialpe domenii neregulate ale spatiului nu este o problema chiar simpla.Aceste functii trebuie sa garanteze prin structura lor anumite pro-prietati de continuitate ale aproximarii necunoscutei (necunoscute-lor) problemei la traversarea frontierei dintre elemente. Construireaunor asemenea functii pe domenii de forma simpla nu ridica pro-bleme majore. In schimb, atunci cand se ıncearca definirea lor directpe domenii neregulate apar dificultati de natura matematica.

Vom vedea ın acest paragraf maniera foarte eleganta ın careMEF depaseste aceste dificultati. In esenta, definirea functiilor deaproximare elementale se realizeaza ın doua etape:

• se construieste mai ıntai o aproximanta polinomiala pe unelement-tip de forma foarte simpla (numit element-parinte);

• elementul-parinte este apoi”strambat” ın asa fel ıncat sa ca-

pete forma si pozitia elementului din spatiul real (spatiul ıncare se desfasoara operatia de discretizare propriu-zisa).

”Strambarea” se traduce ın limbaj matematic prin utilizarea unei

transformari de coordonate care leaga spatiul elementului-parintede spatiul real. Si aici MEF recurge la un artificiu ingenios. Pentrudefinirea transformarii de coordonate sunt folosite chiar functiile deforma ale elementului-parinte.

In capitolele precedente am vazut deja cum se construiesc sise utilizeaza elementele unidimensionale (mai precis, am ıntalnitelementul unidimensional liniar). Din aceasta cauza, nu vom relua

Page 69: Metoda Elementelor Finite

Constructia aproximantelor elementale 63

prezentarea respectivelor elemente. In schimb, vom insista asupraelementelor bi- si tridimensionale, cu atat mai mult cu cat ele aparcel mai frecvent ın problemele practice pe care trebuie sa le rezolveinginerii mecanici.

Spatiul limitat nu ne va permite sa analizam decat un numar re-lativ restrans de elemente: elementul triunghiular liniar, elementulpatrulater biliniar, elementul tetraedric liniar si elementul hexae-dric triliniar. Le-am ales pe acestea din trei motive: ele sunt foartedes utilizate ın practica, sunt cele mai simple din familiile lor si, dinpunct de vedere al tratarii matematice, sunt reprezentative pentrualte tipuri de elemente mai complexe (ın sensul ca ne dau o ideeclara asupra manierei ın care pot fi construite functiile de aproxi-mare si pentru elementele de ordin superior).

3.2.1 Elementul triunghiular liniar

Am vazut ca elementul triunghiular liniar este o regiune planade forma triunghiulara, cu laturi drepte si trei noduri plasate ınvarfuri (fig. 3.2.a). Pentru a construi functia de aproximare poli-nomiala asociata acestui element vom parcurge pe rand cele douaetape mentionate mai sus.

Vom ıncepe prin a construi o astfel de functie pentru un elementtriunghiular de forma foarte simpla. Acest element este prezentatın figura 3.6. Observam ca el are forma unui triunghi dreptunghicisoscel ale carui catete au lungimea egala cu unu. Elementul esteplasat cu unghiul drept ın originea unui sistem de coordonate nor-malizate (ξ, η). Nodurile sale au fost notate conventional cu 1, 2, res-pectiv 3. Nu am mai folosit notatia literala din capitolele precedentedin dorinta de a simplifica scrierea relatiilor matematice. Sa retinemınsa ca etichetele numerice din figura 3.6 sunt conventionale, deci nureprezinta neaparat numerele de ordine ale nodurilor unui anumitelement din spatiul real.

Vom nota cu ϕ marimea necunoscuta care trebuie aproximatala nivelul elementului. Valorile lui ϕ ın cele trei noduri vor fi de-

Page 70: Metoda Elementelor Finite

64 Capitolul 3

ξ13

2

η

(0, 1)

(0, 0) (1, 0)

Figura 3.6: Elementul-parinte folosit pentru construirea functiei deaproximare a unui element triunghiular liniar

semnate prin simbolurile ϕi(i = 1, 2, 3). Functia polinomiala care ılaproximeaza pe ϕ va fi notata ϕ.

In general, pentru determinarea coeficientilor lui ϕ se impuneconditia ca aceasta sa ia valorile ϕi(i = 1, 2, 3) ın nodurile cores-punzatoare ale elementului. Deci va exista o relatie bine definitaıntre numarul nodurilor elementului si gradul polinomului de apro-ximare asociat. In cazul elementului din figura 3.6, polinomul ϕtrebuie sa aiba trei coeficienti (exact numarul de noduri). Expresialui va fi atunci

ϕ(ξ, η) = a0 + a1ξ + a2η, ξ ∈ [0, 1], η ∈ [0, 1 − ξ]. (3.1)

Observam ca ϕ este un polinom de gradul ıntai ın variabilele ξ, η.

Page 71: Metoda Elementelor Finite

Constructia aproximantelor elementale 65

Pentru determinarea coeficientilor a0, a1 si a2 impunem con-ditiile (vezi figura 3.6)

ϕ(1, 0) = ϕ1,

ϕ(0, 1) = ϕ2,

ϕ(0, 0) = ϕ3.

(3.2)

Cu ajutorul lui (3.1) ele pot fi transcrise sub forma unui sistem deecuatii liniare ın necunoscutele a0, a1 si a2:

a0 +a1 = ϕ1,

a0 +a2 = ϕ2,

a0 = ϕ3.

(3.3)

Solutia lui este

a0 = ϕ3, a1 = ϕ1 − ϕ3, a2 = ϕ2 − ϕ3. (3.4)

Inlocuind (3.4) ın (3.1), obtinem urmatoarea expresie a aproximan-tei ϕ:

ϕ(ξ, η) =3X

i=1

ϕiNi(ξ, η), ξ ∈ [0, 1], η ∈ [0, 1 − ξ], (3.5)

ın care

N1(ξ, η) = ξ, N2(ξ, η) = η, N3(ξ, η) = 1 − ξ − η (3.6)

sunt functiile de forma asociate celor trei noduri ale elementului dinfigura 3.6. Am mai ıntalnit aceasta notiune ın capitolul 2, cand amprezentat elementul unidimensional liniar folosit pentru modelareacorzii flexibile.

Functia ϕ definita prin relatia (3.5) este rezultatul primei etapede construire a aproximantei elementale. Ea este capabila sa descrievariatia necunoscutei ϕ la nivelul elementului de forma particularadin figura 3.6.

Page 72: Metoda Elementelor Finite

66 Capitolul 3

ξ

η

1

2

3

y

1

3

2

x

x = x(ξ, η)y = y(ξ, η) (Ω)

Figura 3.7: Transpunerea elementului-parinte ın spatiul real al ele-mentului triunghiular liniar

Vom continua cu etapa”strambarii” elementului-parinte ın asa

fel ıncat el sa devina coincident cu un element oarecare din spatiulreal. Va trebui sa realizam operatia reprezentata ın figura 3.7.

In limbaj matematic, aceasta operatie ınseamna construireaunei transformari biunivoce care sa lege coordonatele reale (x, y)de coordonatele normalizate (ξ, η) si care sa garanteze obtinereadomeniului plan Ω ocupat de elementul real. O asemenea trans-formare de coordonate are expresia generala

x = x(ξ, η),

y = y(ξ, η), ξ ∈ [0, 1], η ∈ [0, 1 − ξ].(3.7)

Pentru ca ea sa fie biunivoca, este necesar ca jacobianul asociat safie strict pozitiv:

Page 73: Metoda Elementelor Finite

Constructia aproximantelor elementale 67

J =D(x, y)

D(ξ, η)=

∂ x∂ ξ

∂ x∂ η

∂ y∂ ξ

∂ y∂ η

> 0. (3.8)

In aceste conditii, transformarea (3.7) poate fi inversata, deci sepoate gasi corespondenta

ξ = ξ(x, y),

η = η(x, y), (x, y) ∈ Ω.(3.9)

Inlocuind dependentele (3.9) ın (3.5), se obtine expresia aproximan-tei ϕ ca functie de coordonatele reale x si y:

ϕ(x, y) =3X

i=1

ϕiNi [ξ(x, y), η(x, y)] . (3.10)

Aceasta ar fi calea de urmat. Problema-cheie aici este gasireaefectiva a transformarii de coordonate (3.7). Pentru rezolvarea ei,MEF propune folosirea unor functii identice cu aproximanta poli-nomiala (3.5). Cu alte cuvinte, se adopta urmatoarea transformarede coordonate:

x(ξ, η) =3X

i=1

xiNi(ξ, η),

y(ξ, η) =3X

i=1

yiNi(ξ, η), ξ ∈ [0, 1], η ∈ [0, 1 − ξ],

(3.11)

ın care Ni(ξ, η)(i = 1, 2, 3) sunt functiile de forma definite derelatiile (3.6), iar xi si yi (i = 1, 2, 3) sunt coordonatele celortrei noduri ale elementului real. Aceste coordonate sunt ın gene-ral cunoscute odata ce discretizarea domeniului de analiza a fostrealizata.

Observam ca, de fapt, transformarea (3.11) nu este altceva decato interpolare liniara a coordonatelor nodurilor asociate elementului.Ne dam astfel seama ca forma elementului real este complet definitaprin pozitiile nodurilor sale.

Page 74: Metoda Elementelor Finite

68 Capitolul 3

Sa vedem acum ın ce conditii transformarea (3.11) este bi-univoca. In acest scop, vom explicita jacobianul (3.8). Din (3.11) si(3.6) rezulta prin calcul direct urmatoarele expresii ale derivatelor∂ x/∂ ξ, ∂ x/∂ η, ∂ y/∂ ξ si ∂ y/∂ η:

∂ x

∂ ξ=

3Xi=1

xi

∂ Ni

∂ ξ= x1 − x3,

∂ x

∂ η=

3Xi=1

xi

∂ Ni

∂ η= x2 − x3,

∂ y

∂ ξ=

3Xi=1

yi

∂ Ni

∂ ξ= y1 − y3,

∂ y

∂ η=

3Xi=1

yi

∂ Ni

∂ η= y2 − y3.

(3.12)

Din (3.8) si (3.12) se obtine atunci

J =

x1 − x3 x2 − x3

y1 − y3 y2 − y3

= 1 1 1

x1 x2 x3

y1 y2 y3

. (3.13)

Pentru ca jacobianul (3.13) sa fie strict pozitiv este necesar sisuficient ca nodurile elementului real sa fie necoliniare si nume-rotate ın succesiune antiorara (vezi figura 3.7). In aceste conditii, Jare o interpretare geometrica foarte simpla. El este, de fapt, dublulariei elementului real:

J = 2A > 0. (3.14)

Elementul triunghiular liniar este un caz fericit, fiindca trans-formarea de coordonate (3.11) poate fi inversata analitic. Pentruelementele triunghiulare de ordin superior inversarea nu se maipoate face decat numeric si, din aceasta cauza, este lasata pe seamaprogramului MEF.

Page 75: Metoda Elementelor Finite

Constructia aproximantelor elementale 69

Vom vedea ın continuare cum se inverseaza (3.11) si care esteexpresia aproximantei ϕ ca dependenta explicita de coordonatelereale (x, y).

Cu ajutorul lui (3.6) se poate obtine urmatoarea expresie atransformarii (3.11):

(x1 − x3)ξ + (x2 − x3)η = x − x3,

(y1 − y3)ξ + (y2 − y3)η = y − y3.(3.15)

(3.15) poate fi privit ca un sistem de ecuatii liniare ın necunoscuteleξ si η. Observam ca determinantul sau este chiar jacobianul J (vezi(3.13)). Fiindca el este strict pozitiv, (3.15) are o solutie unica.Aceasta poate fi obtinuta cu usurinta fie prin eliminare, fie cu regulalui Cramer:

ξ(x, y) =1

2A[(x2y3 − x3y2) + (y2 − y3) x + (x3 − x2) y] ,

η(x, y) =1

2A[(x3y1 − x1y3) + (y3 − y1) x + (x1 − x3) y] ,

(x, y) ∈ Ω.

(3.16)

(3.16) este, de fapt, transformarea inversa (3.9) pe care trebuie sao ınlocuim ın expresia (3.5) a aproximantei polinomiale ϕ, pentru ao exprima ca dependenta de coordonatele reale (x, y). Vom obtineastfel o functie de aproximare

ϕ(x, y) =3X

i=1

ϕiNi(x, y), (x, y) ∈ Ω, (3.17)

ın careNi(x, y) = Ni [ξ(x, y), η(x, y)] , i = 1, 2, 3, (3.18)

sunt functiile de forma asociate celor trei noduri ale elementului,scrise de data aceasta ın termenii coordonatelor reale (x, y). Expre-siile lor rezulta prin ınlocuirea relatiilor (3.16) ın (3.6):

Ni(x, y) =1

2A(bi + cix + diy) , i = 1, 2, 3, (3.19)

unde bi, ci, di (i = 1, 2, 3) sunt coeficienti definiti astfel:

Page 76: Metoda Elementelor Finite

70 Capitolul 3

b1 = x2y3 − x3y2, c1 = y2 − y3, d1 = x3 − x2,

b2 = x3y1 − x1y3, c2 = y3 − y1, d2 = x1 − x3,

b3 = x1y2 − x2y1, c3 = y1 − y2, d3 = x2 − x1.

(3.20)

Daca analizam (3.19) - (3.20), constatam din nou ca denumireade functii de forma data cantitatilor Ni(x, y) (i = 1, 2, 3) este pedeplin justificata, fiindca ın expresiile lor intervin doar coordonatelenodale si aria elementului real.

3.2.2 Elementul patrulater biliniar

Elementul patrulater biliniar este o regiune plana de formapatrulatera, cu laturi rectilinii si patru noduri plasate ın varfuri(fig. 3.3.a).

La fel ca ın §3.2.1, vom ıncepe prin a construi functia de apro-ximare pentru un element patrulater biliniar foarte simplu. Acestaeste prezentat ın figura 3.8. Observam ca elementul are forma unuipatrat cu laturi de lungime egala cu doi, orientate paralalel cu axeleunui sistem de coordonate normalizate (ξ, η). Centrul elementuluicoincide cu originea sistemului de coordonate. Nodurile sale au fostnumerotate conventional cu 1, 2, 3, respectiv 4. In cele ce urmeazavom folosi notatii asemanatoare cu cele din §3.2.1.

Fiindca elementul din figura 3.8 are patru noduri, polinomulde aproximare ϕ va trebui sa aiba patru coeficienti. Putem adoptaatunci urmatoarea expresie a lui ϕ:

ϕ(ξ, η) = a0 + a1ξ + a2η + a3ξη, ξ, η ∈ [−1, 1]. (3.21)

Functia definita de (3.21) are o serie de proprietati interesante:

• ea este o dependenta de gradul al doilea ın raport cu ansam-blul variabilelor ξ si η, fara a fi ınsa un polinom de gradul aldoilea complet (ıi lipsesc termenii corespunzatori lui ξ2 si η2);

• daca ıl pastram constant pe η, ramanem cu o dependentaliniara de ξ;

Page 77: Metoda Elementelor Finite

Constructia aproximantelor elementale 71

ξ

(+1,−1)

(+1, +1)(−1, +1)

(−1,−1)

3

21

4

η

Figura 3.8: Elementul-parinte folosit pentru construirea functiei deaproximare a unui element patrulater biliniar

• daca ıl pastram constant pe ξ, ramanem cu o dependentaliniara de η.

Un polinom ın doua variabile care are proprietatile enumerate maisus este cunoscut sub numele de functie biliniara. De aici provinesi denumirea elementului patrulater asociat.

Pentru a determina coeficientii a0, a1, a2 si a3, impunem conditiaca polinomul ϕ definit de (3.21) sa ia valorile necunoscutei ϕ ın celepatru noduri ale elementului din figura 3.8:

Page 78: Metoda Elementelor Finite

72 Capitolul 3

ϕ(−1,−1) = ϕ1,

ϕ(+1,−1) = ϕ2,

ϕ(+1, +1) = ϕ3,

ϕ(−1, +1) = ϕ4.

(3.22)

Din (3.22) si (3.21) obtinem un sistem de ecuatii liniare ın ne-cunoscutele a0, a1, a2 si a3:

a0 − a1 − a2 + a3 = ϕ1,

a0 + a1 − a2 − a3 = ϕ2,

a0 + a1 + a2 + a3 = ϕ3,

a0 − a1 + a2 − a3 = ϕ4.

(3.23)

Solutia lui este

a0 =1

4(ϕ1 + ϕ2 + ϕ3 + ϕ4) ,

a1 =1

4(−ϕ1 + ϕ2 + ϕ3 − ϕ4) ,

a2 =1

4(−ϕ1 − ϕ2 + ϕ3 + ϕ4) ,

a3 =1

4(ϕ1 − ϕ2 + ϕ3 − ϕ4) .

(3.24)

Prin ınlocuirea lui (3.24) ın (3.21) rezulta urmatoarea expresie alui ϕ:

ϕ(ξ, η) =4X

i=1

ϕiNi(ξ, η), ξ, η ∈ [−1, 1], (3.25)

ın care

Page 79: Metoda Elementelor Finite

Constructia aproximantelor elementale 73

N1(ξ, η) =1

4(1 − ξ) (1 − η) ,

N2(ξ, η) =1

4(1 + ξ) (1 − η) ,

N3(ξ, η) =1

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

N4(ξ, η) =1

4(1 − ξ) (1 + η)

(3.26)

sunt functiile de forma asociate celor patru noduri ale elementuluidin figura 3.8.

Vom trece acum la”strambarea” elementului-parinte ın asa fel

ıncat el sa devina coincident cu un element oarecare din spatiul real(fig. 3.9). In acest scop, vom folosi o transformare de coordonatesimilara ca structura aproximantei definite de (3.25):

x(ξ, η) =4X

i=1

xiNi(ξ, η),

y(ξ, η) =4X

i=1

yiNi(ξ, η), ξ, η ∈ [−1, 1],

(3.27)

ın care Ni(ξ, η)(i = 1, 2, 3, 4) sunt functiile de forma definite derelatiile (3.26), iar xi si yi (i = 1, 2, 3, 4) sunt coordonatele celorpatru noduri ale elementului real din figura 3.9 (ın general, acestecoordonate sunt marimi cunoscute).

Fara a mai intra ın detalii, vom mentiona ca jacobianul asociatlui (3.27) este strict pozitiv daca nodurile elementului real suntnumerotate ın succesiune antiorara, sunt necoliniare toate patru,iar laturile succesive ale elementului nu formeaza unghiuri internemai mari decat 180. Ca exemplu, ın figura 3.10 prezentam unelement patrulater biliniar degenerat ca urmare a ıncalcarii ulti-mei conditii. Intr-o astfel de situatie, transformarea de coordonate(3.27) ısi pierde biunivocitatea (jacobianul ei devine negativ).

In ciuda caracterului sau neliniar, transformarea (3.27) poate fiinversata analitic. Totusi, expresia inversei este atat de complicata,

Page 80: Metoda Elementelor Finite

74 Capitolul 3

ξ

x

y

η

(Ω)

2

4

1

3

3

x = x(ξ, η)y = y(ξ, η) 1

4

2

Figura 3.9: Transpunerea elementului-parinte ın spatiul real al ele-mentului patrulater biliniar

Page 81: Metoda Elementelor Finite

Constructia aproximantelor elementale 75

α > 180

Figura 3.10: Exemplu de element patrulater biliniar degenerat

ıncat explicitarea ei devine cu totul lipsita de atractivitate, maiales daca stim ca operatia poate fi realizata numeric de programeleMEF. Pentru scopurile noastre este suficient sa retinem conditiilede inversabilitate ale lui (3.27) si sa cunoastem modul ın care estefolosita inversa mai departe.

3.2.3 Elementul tetraedric liniar

Elementul tetraedric liniar este o regiune spatiala marginita depatru fete plane triunghiulare. El are patru noduri, cate unul plasatın fiecare varf (fig. 3.4.a).

Procedand ın maniera deja uzuala, vom ıncepe prin a construifunctia de aproximare pentru cazul unui element tetraedric liniarfoarte simplu. Acest element este prezentat ın figura 3.11. Obser-vam ca el are trei muchii reciproc perpendiculare de lungime egalacu unu. Varful comun acestor muchii coincide cu originea unui sis-tem de coordonate normalizate (ξ, η, ζ). Nodurile elementului aufost notate conventional cu 1, 2, 3, respectiv 4.

Fiindca elementul din figura 3.11 are patru noduri, polinomulde aproximare ϕ poate primi urmatoarea expresie:

Page 82: Metoda Elementelor Finite

76 Capitolul 3

η

3

4

(0, 0, 1)

ξ

ζ

1

(1, 0, 0)

(0, 0, 0)

2

(0, 1, 0)

Figura 3.11: Elementul-parinte folosit pentru construirea functieide aproximare a unui element tetraedric liniar

ϕ(ξ, η, ζ) = a0 + a1ξ + a2η + a3ζ,

ξ ∈ [0, 1], η ∈ [0, 1 − ξ], ζ ∈ [0, 1 − ξ − η].(3.28)

Functia definita de (3.28) este un polinom de gradul ıntai ın vari-abilele ξ, η si ζ, fapt care justifica denumirea elementului tetraedricasociat.

Pentru determinarea coeficientilor a0, a1, a2 si a3 impunemconditia ca functia ϕ sa ia valorile necunoscutei ϕ ın cele patrunoduri din figura 3.11:

Page 83: Metoda Elementelor Finite

Constructia aproximantelor elementale 77

ϕ(1, 0, 0) = ϕ1,

ϕ(0, 1, 0) = ϕ2,

ϕ(0, 0, 1) = ϕ3,

ϕ(0, 0, 0) = ϕ4.

(3.29)

Din (3.29) si (3.28) obtinem un sistem de ecuatii ın necunoscutelea0, a1, a2 si a3:

a0 +a1 = ϕ1,

a0 +a2 = ϕ2,

a0 +a3 = ϕ3,

a0 = ϕ4.

(3.30)

Solutia lui este

a0 = ϕ4,

a1 = ϕ1 − ϕ4,

a2 = ϕ2 − ϕ4,

a3 = ϕ3 − ϕ4.

(3.31)

Prin ınlocuirea lui (3.31) ın (3.28), rezulta urmatoarea expresie aaproximantei ϕ:

ϕ(ξ, η, ζ) =4X

i=1

ϕiNi(ξ, η, ζ),

ξ ∈ [0, 1], η ∈ [0, 1 − ξ], ζ ∈ [0, 1 − ξ − η],

(3.32)

ın care

N1(ξ, η, ζ) = ξ,

N2(ξ, η, ζ) = η,

N3(ξ, η, ζ) = ζ,

N4(ξ, η, ζ) = 1 − ξ − η − ζ

(3.33)

sunt functiile de forma asociate nodurilor elementului din figura3.11.

Page 84: Metoda Elementelor Finite

78 Capitolul 3

Vom trece acum la”strambarea” elementului-parinte ın asa fel

ıncat el sa devina coincident cu un element oarecare din spatiul real(fig. 3.12). In acest scop, vom folosi o transformare de coordonatesimilara ca structura cu (3.11) si (3.27). Evident ınsa, ea va implicaacum trei ecuatii:

x(ξ, η, ζ) =4X

i=1

xiNi(ξ, η, ζ),

y(ξ, η, ζ) =4X

i=1

yiNi(ξ, η, ζ),

z(ξ, η, ζ) =4X

i=1

ziNi(ξ, η, ζ),

ξ ∈ [0, 1], η ∈ [0, 1 − ξ], ζ ∈ [0, 1 − ξ − η],

(3.34)

ın care Ni(ξ, η, ζ) (i = 1, 2, 3, 4) sunt functiile de forma (3.33), iarxi, yi si zi (i = 1, 2, 3, 4) sunt coordonatele celor patru noduri aleelementului real din figura 3.12 (ın general, aceste coordonate suntmarimi cunoscute). Pentru ca transformarea (3.34) sa fie biunivocaeste necesar ca jacobianul asociat sa fie strict pozitiv:

J =D(x, y, z)

D(ξ, η, ζ)=

∂ x∂ ξ

∂ x∂ η

∂ x∂ ζ

∂ y∂ ξ

∂ y∂ η

∂ y∂ ζ

∂ z∂ ξ

∂ z∂ η

∂ z∂ ζ

> 0. (3.35)

Observam din nou ca transformarea (3.34) nu este altceva decato interpolare liniara a coordonatelor nodale. Ne dam astfel seamaca forma elementului real este complet determinata prin pozitiilenodurilor sale.

Sa vedem acum ın ce conditii transformarea (3.34) este bi-univoca. Pentru aceasta vom explicita jacobianul (3.35). Din (3.34)si (3.33) rezulta prin calcul direct urmatoarele expresii ale deri-vatelor ∂ x/∂ ξ, ∂ x/∂ η, ∂ x/∂ ζ, ∂ y/∂ ξ, ∂ y/∂ η, ∂ y/∂ ζ, ∂ z/∂ ξ,∂ z/∂ η si ∂ z/∂ ζ:

Page 85: Metoda Elementelor Finite

Constructia aproximantelor elementale 79

ξ

24

y

x

z

3

(Ω)

z = z(ξ, η, ζ)

ζ

3

1

η

4y = y(ξ, η, ζ)x = x(ξ, η, ζ)

1

2

Figura 3.12: Transpunerea elementului-parinte ın spatiul real al ele-mentului tetraedric liniar

Page 86: Metoda Elementelor Finite

80 Capitolul 3

∂ x

∂ ξ=

4Xi=1

xi

∂ Ni

∂ ξ= x1 − x4,

∂ x

∂ η=

4Xi=1

xi

∂ Ni

∂ η= x2 − x4,

∂ x

∂ ζ=

4Xi=1

xi

∂ Ni

∂ ζ= x3 − x4,

∂ y

∂ ξ=

4Xi=1

yi

∂ Ni

∂ ξ= y1 − y4,

∂ y

∂ η=

4Xi=1

yi

∂ Ni

∂ η= y2 − y4,

∂ y

∂ ζ=

4Xi=1

yi

∂ Ni

∂ ζ= y3 − y4,

∂ z

∂ ξ=

4Xi=1

zi

∂ Ni

∂ ξ= z1 − z4,

∂ z

∂ η=

4Xi=1

zi

∂ Ni

∂ η= z2 − z4,

∂ z

∂ ζ=

4Xi=1

zi

∂ Ni

∂ ζ= z3 − z4.

(3.36)

Din (3.35) - (3.36) se obtine atunci

J =

x1 − x4 x2 − x4 x3 − x4

y1 − y4 y2 − y4 y3 − y4

z1 − z4 z2 − z4 z3 − z4

= 1 1 1 1

x1 x2 x3 x4

y1 y2 y3 y4

z1 z2 z3 z4

. (3.37)

Pentru ca jacobianul (3.37) sa fie strict pozitiv este necesar sisuficient ca nodurile elementului real sa fie necoplanare si numero-tate ın succesiunea prezentata ın figura 3.12. In practica se adopta

Page 87: Metoda Elementelor Finite

Constructia aproximantelor elementale 81

urmatoarea regula de etichetare a nodurilor: se alege un nod oare-care al elementului tetraedric si i se atribuie indicele 1; ın conti-nuare, privind din acest nod ınspre fata opusa a tetraedrului, nodu-rile ramase primesc indicii 2, 3, respectiv 4 ın succesiune antiorara.

Daca este strict pozitiv, jacobianul J are o interpretare geo-metrica foarte simpla. El este, de fapt, de sase ori volumul elementu-lui tetraedric real:

J = 6V > 0. (3.38)

Si ın cazul elementului tetraedric liniar transformarea de co-ordonate (3.34) poate fi inversata analitic. Totusi, obtinerea inver-sei pe aceasta cale presupune un volum de calcule destul de mare,motiv pentru care ea este lasata pe seama programului MEF. Pen-tru elementele tetraedrice de ordin superior inversarea nu se maipoate face decat numeric.

3.2.4 Elementul hexaedric triliniar

Elementul hexaedric triliniar este o regiune spatiala cu sase fetepatrulatere si muchii drepte, avand opt noduri plasate ın varfuri(fig. 3.5).

Vom ıncepe prin a construi functia de aproximare pentru cazulunui element hexaedric triliniar foarte simplu. Acest element esteprezentat ın figura 3.13. Observam ca el are forma unui cub culaturile de lungime egala cu doi, orientate paralel cu axele unuisistem de coordonate normalizate (ξ, η, ζ). Centrul elementului estesuprapus peste originea sistemului de coordonate. Nodurile sale aufost notate conventional cu 1, 2, . . . , 8.

Fiindca elementul din figura 3.13 are opt noduri, polinomul deaproximare ϕ va trebui sa aiba opt coeficienti. Putem adopta atunciurmatoarea expresie a lui ϕ:

ϕ (ξ, η, ζ) = a0 + a1ξ + a2η + a3ζ+

+ a4ξη + a5ηζ + a6ζξ + a7ξηζ,

ξ, η, ζ ∈ [−1, 1].

(3.39)

Page 88: Metoda Elementelor Finite

82 Capitolul 3

1 → (−1,−1,−1)2 → (+1,−1,−1)3 → (+1, +1,−1)4 → (−1, +1,−1)5 → (−1,−1, +1)6 → (+1,−1, +1)7 → (+1, +1, +1)8 → (−1, +1, +1)

ξ

η

ζ

1

3

4

2

5

6 7

8

Figura 3.13: Elementul-parinte folosit pentru construirea functieide aproximare a unui element hexaedric triliniar

Page 89: Metoda Elementelor Finite

Constructia aproximantelor elementale 83

Functia definita de (3.39) are o serie de proprietati interesante:

• ea este o dependenta de gradul al treilea ın raport cu ansam-blul variabilelor ξ, η si ζ, fara a fi ınsa un polinom de gradulal treilea complet;

• daca pastram constante doua dintre coordonate, ramanem cuo dependenta liniara ın raport cu cea de a treia.

Un polinom ın trei variabile care are proprietatile enumerate maisus este cunoscut sub numele de functie triliniara. De aici provinesi denumirea elementului hexaedric asociat.

Pentru determinarea coeficientilor a0, a1, . . . , a7 impunemconditia ca polinomul ϕ sa ia valorile necunoscutei ϕ ın cele optnoduri din figura 3.13:

ϕ(−1,−1,−1) = ϕ1, ϕ(+1,−1,−1) = ϕ2,

ϕ(+1, +1,−1) = ϕ3, ϕ(−1, +1,−1) = ϕ4,

ϕ(−1,−1, +1) = ϕ5, ϕ(+1,−1, +1) = ϕ6,

ϕ(+1, +1, +1) = ϕ7, ϕ(−1, +1, +1) = ϕ8.

(3.40)

Din (3.40) si (3.39) obtinem un sistem de ecuatii liniare ın necu-noscutele a0, a1, . . . , a7:

a0 − a1 − a2 − a3 + a4 + a5 + a6 − a7 = ϕ1,

a0 + a1 − a2 − a3 − a4 + a5 − a6 + a7 = ϕ2,

a0 + a1 + a2 − a3 + a4 − a5 − a6 − a7 = ϕ3,

a0 − a1 + a2 − a3 − a4 − a5 + a6 + a7 = ϕ4,

a0 − a1 − a2 + a3 + a4 − a5 − a6 + a7 = ϕ5,

a0 + a1 − a2 + a3 − a4 − a5 + a6 − a7 = ϕ6,

a0 + a1 + a2 + a3 + a4 + a5 + a6 + a7 = ϕ7,

a0 − a1 + a2 + a3 − a4 + a5 − a6 − a7 = ϕ8.

(3.41)

Page 90: Metoda Elementelor Finite

84 Capitolul 3

Solutia lui este

a0 = (ϕ1 + ϕ2 + ϕ3 + ϕ4 + ϕ5 + ϕ6 + ϕ7 + ϕ8) / 8,

a1 = (−ϕ1 + ϕ2 + ϕ3 − ϕ4 − ϕ5 + ϕ6 + ϕ7 − ϕ8) / 8,

a2 = (−ϕ1 − ϕ2 + ϕ3 + ϕ4 − ϕ5 − ϕ6 + ϕ7 + ϕ8) / 8,

a3 = (−ϕ1 − ϕ2 − ϕ3 − ϕ4 + ϕ5 + ϕ6 + ϕ7 + ϕ8) / 8,

a4 = (ϕ1 − ϕ2 + ϕ3 − ϕ4 + ϕ5 − ϕ6 + ϕ7 − ϕ8) / 8,

a5 = (ϕ1 + ϕ2 − ϕ3 − ϕ4 − ϕ5 − ϕ6 + ϕ7 + ϕ8) / 8,

a6 = (ϕ1 − ϕ2 − ϕ3 + ϕ4 − ϕ5 + ϕ6 + ϕ7 − ϕ8) / 8,

a7 = (−ϕ1 + ϕ2 − ϕ3 + ϕ4 + ϕ5 − ϕ6 + ϕ7 − ϕ8) / 8.

(3.42)

Prin ınlocuirea lui (3.42) ın (3.39) rezulta urmatoarea expresie apolinomului aproximant ϕ:

ϕ (ξ, η, ζ) =8X

i=1

ϕiNi (ξ, η, ζ) , ξ, η, ζ ∈ [−1, 1], (3.43)

ın care

N1 (ξ, η, ζ) = (1 − ξ) (1 − η) (1 − ζ) / 8,

N2 (ξ, η, ζ) = (1 + ξ) (1 − η) (1 − ζ) / 8,

N3 (ξ, η, ζ) = (1 + ξ) (1 + η) (1 − ζ) / 8,

N4 (ξ, η, ζ) = (1 − ξ) (1 + η) (1 − ζ) / 8,

N5 (ξ, η, ζ) = (1 − ξ) (1 − η) (1 + ζ) / 8,

N6 (ξ, η, ζ) = (1 + ξ) (1 − η) (1 + ζ) / 8,

N7 (ξ, η, ζ) = (1 + ξ) (1 + η) (1 + ζ) / 8,

N8 (ξ, η, ζ) = (1 − ξ) (1 + η) (1 + ζ) / 8

(3.44)

sunt functiile de forma asociate celor opt noduri ale elementului dinfigura 3.13.

Page 91: Metoda Elementelor Finite

Constructia aproximantelor elementale 85

Vom trece acum la”strambarea” elementului-parinte ın asa fel

ıncat el sa devina coincident cu un element oarecare din spatiul real(fig. 3.14). In acest scop, vom folosi o transformare de coordonateasemanatoare ca structura aproximantei definite de (3.43):

x(ξ, η, ζ) =8X

i=1

xiNi(ξ, η, ζ),

y(ξ, η, ζ) =8X

i=1

yiNi(ξ, η, ζ),

z(ξ, η, ζ) =8X

i=1

ziNi(ξ, η, ζ),

ξ, η, ζ ∈ [−1, 1],

(3.45)

ın care Ni(ξ, η, ζ) (i = 1, 2, . . . , 8) sunt functiile de forma (3.44), iarxi, yi si zi (i = 1, 2, . . . , 8) sunt coordonatele celor opt noduri aleelementului real din figura 3.14 (ın general, aceste coordonate suntmarimi cunoscute).

Fara a mai intra ın detalii, vom mentiona ca jacobianul asociatlui (3.45) este strict pozitiv daca nodurile elementului real sunt nu-merotate ın succesiunea prezentata ın figura 3.14, sunt necoplanaretoate opt, iar unghiurile diedre interne formate de fetele adiacentenu sunt mai mari decat 180. In practica se utilizeaza urmatoarearegula de etichetare a nodurilor: plasandu-ne imaginar ın interio-rul elementului, selectam doua fete opuse ale acestuia; privind spreuna din aceste fete, alegem un nod oarecare si ıi atribuim indicele1, apoi asociem indicii 2, 3, respectiv 4 celorlalte noduri ale feteiparcurgandu-i conturul ın sens antiorar; ıntorcandu-ne spre cea-lalta fata selectata, atribuim indicele 5 nodului legat prin muchiede nodul 1, dupa care, parcurgand conturul fetei ın sens antiorar,atribuim celorlalte noduri indicii 6, 7, respectiv 8.

In general, inversarea transformarii de coordonate (3.45) esterealizata pe cale numerica de programul MEF atat pentru elementultriliniar, cat si pentru elementele hexaedrice de ordin superior.

Page 92: Metoda Elementelor Finite

86 Capitolul 3

ξ

η

ζ

1

2 3

4

5

6 7

8

yx

z

(Ω)

1

2

3

4

5

6

7

8

x = x(ξ, η, ζ)y = y(ξ, η, ζ)z = z(ξ, η, ζ)

Figura 3.14: Transpunerea elementului-parinte ın spatiul real al ele-mentului hexaedric triliniar

Page 93: Metoda Elementelor Finite

Aproximarea marimilor vectoriale 87

3.3 Aproximarea marimilor vectoriale

La construirea aproximantelor prezentate ın paragrafele ante-rioare am presupus tacit ca marimea necunoscuta ϕ este un scalar.Exista ınsa un mare numar de aplicatii ın care apar necunoscute detip vectorial. De exemplu, ın cazul problemelor de elasticitate, ne-cunoscuta primara este campul vectorilor deplasare ai particulelor.La abordarea unor asemenea probleme, va fi deci necesara aproxi-marea unor marimi de tip vectorial la nivelul fiecarui element.

In general, pentru construirea aproximantei unei functii vecto-riale se procedeaza ın felul urmator: respectiva marime se descom-pune ıntr-un sistem de coordonate convenabil ales, iar componenteleastfel obtinute se aproximeaza ca si marimile scalare obisnuite.

Peste tot ın acest curs vom nota cu ~d campul vectorilor depla-sare. Componentele lui ~d ıntr-un sistem de coordonate cartezieneortogonale vor fi simbolizate astfel:

u – deplasarea ın directia axei x;

v – deplasarea ın directia axei y;

w – deplasarea ın directia axei z.

Evident, ın cazul problemelor bidimensionale (cum sunt cele destare plana de tensiuni sau stare plana de deformatii), doar doua

componente ale lui ~d sunt nenule (noi vom considera ca, ın aseme-nea situatii, componentele nenule sunt u si v, deci w = 0).

Sa ne ocupam acum ın detaliu de problema aproximarii necunos-cutelor de tip vectorial. Vom ıncepe prin a aborda un caz bidimen-sional, dupa care vom generaliza procedura la cazul tridimensional.

Ne propunem sa construim aproximanta marimii vectoriale~d(u, v) la nivelul elementului triunghiular liniar din figura 3.15. Sa

observam pe aceasta figura faptul ca vectorul ~d asociat fiecarui nodal elementului a fost descompus ın raport cu axele unui sistem decoordonate plan. Componentele astfel obtinute au fost etichetatefolosind indicele nodului corespunzator.

Page 94: Metoda Elementelor Finite

88 Capitolul 3

x

y

1

2

3

u1

v1

u2

v2

u3

v3

(Ω)

Figura 3.15: Aproximarea unei marimi vectoriale ın interiorul unuielement triunghiular liniar (valorile nodale ale componentelor u siv sunt reprezentate prin sageti paralele cu axele x, respectiv y)

Dupa cum am mentionat mai sus, componentele u si v alemarimii necunoscute ~d sunt aproximate ca si marimile scalareobisnuite. De fapt, ele sunt propriu-zis scalari, fiindca au rezultatca proiectii ale unui vector pe axele de coordonate. Folosind atuncirelatia (3.17), putem scrie urmatoarele formule de aproximare pen-tru respectivele componente:

u(x, y) =3X

i=1

uiNi(x, y), v(x, y) =3X

i=1

viNi(x, y),

(x, y) ∈ Ω,

(3.46)

ın care apar functiile de forma definite de (3.19) si (3.20).

Page 95: Metoda Elementelor Finite

Aproximarea marimilor vectoriale 89

y

z

1

2

3

4

v1

w1

u2

u1

v2

w2

u3

v3

w3

u4

v4

w4

(Ω)

x

Figura 3.16: Aproximarea unei marimi vectoriale ın interiorul unuielement tetraedric liniar (valorile nodale ale componentelor u, v siw sunt reprezentate prin sageti paralele cu axele x, y, respectiv z)

Putem scrie ecuatii asemanatoare cu (3.46) pentru oricare alt tipde element bidimensional. Evident ınsa ca atunci numarul nodurilorasociate va fi altul, iar functiile de forma vor avea si ele alte expresii.

Sa consideram acum cazul elementului tetraedric liniar (fig.3.16), pentru care dorim sa construim aproximanta unei marimi~d(u, v, w). De aceasta data proiectam necunoscuta vectoriala afiecarui nod pe axele unui sistem cartezian triortogonal. Componen-tele astfel obtinute au ca eticheta indicele nodului corespunzator.

Page 96: Metoda Elementelor Finite

90 Capitolul 3

Pentru fiecare dintre cele trei componente ale marimii ~d vomputea scrie atunci relatii de aproximare de tipul (3.32):

u(x, y, z) =4X

i=1

uiNi(x, y, z), v(x, y, z) =4X

i=1

viNi(x, y, z),

w(x, y, z) =4X

i=1

wiNi(x, y, z), (x, y, z) ∈ Ω,

(3.47)

ın care Ni(x, y, z) (i = 1, 2, 3, 4) sunt functiile de forma asociate no-durilor elementului. Ele au fost deduse ın §3.2.3 ca dependente decoordonatele normalizate ξ, η si ζ (vezi (3.33)). Prin folosirea trans-formarii de coordonate (3.34), functiile de forma pot fi explicitateca dependente de coordonatele reale x, y si z. Am aratat ın §3.2.3ca aceasta operatie este destul de laborioasa si, ın consecinta, estelasata pe seama programului MEF. Sa retinem ınsa ca, atat timpcat transformarea (3.34) este inversabila, exprimarea functiilor deforma (3.33) ca dependente de coordonatele reale este posibila, deciecuatiile (3.47) sunt valide.

Putem scrie relatii asemanatoare cu (3.47) pentru oricare alt tipde element tridimensional. Evident ınsa ca atunci numarul noduri-lor va fi altul, iar functiile de forma vor avea si ele alte expresii.

3.4 Expandarea aproximantelor

elementale

In paragrafele anterioare au fost prezentate functiile de aproxi-mare asociate catorva tipuri de elemente finite uzuale. Expresiileacestor functii contin ca parametri valorile necunoscutei din nodu-rile unui singur element (vezi ecuatiile (3.5), (3.25), (3.32), (3.43),(3.46) si (3.47)). Manevrarea aproximantelor ıntr-o asemenea formada nastere la multe complicatii atunci cand se pune problema ela-borarii teoretice a unui model cu elemente finite. Dificultatile suntgenerate de faptul ca orice model opereaza asupra unor retele care

Page 97: Metoda Elementelor Finite

Expandarea aproximantelor elementale 91

sunt formate din mai multe elemente. Daca functiile de aproximareasociate elementelor sunt exprimate cu ajutorul unor seturi diferitede valori nodale ale necunoscutei, ele devin foarte greu combina-bile ıntr-un sistem de ecuatii global (prin acest termen ıntelegando aproximare extinsa la nivelul ıntregii retele). Sa ne gandim lafaptul ca modelele cu elemente finite sunt ın general rezolvate nu-meric, deci ele trebuie sa fie implementate ın programe. Pentruscrierea unor programe eficiente, este esentiala gasirea unei expre-sii unitare a acestor aproximante, care sa poata fi exploatata dupareguli simple si, mai ales, repetitive.

Rescrierea aproximantelor elementale ıntr-o asemenea manierareprezinta obiectivul prezentului paragraf. Vom vedea ca strategiaadoptata de MEF presupune extinderea formulelor de aproximarededuse de noi anterior ın asa fel ıncat ele sa contina ca parametrivalorile necunoscutei din toate nodurile retelei. Termenii suplimen-tari se obtin pur si simplu prin ınmultirea valorilor necunoscuteidin nodurile neasociate elementului cu zero. Acesti termeni numodifica deci sensul aproximarii elementale. Prin introducerea lorse obtine ınsa o exprimare unitara a tuturor aproximantelor dinmodel, fiindca ele vor opera acum cu toate valorile nodale ale ne-cunoscutei. Transformarile de natura matematica vor fi astfel multsimplificate si usor implementabile ın programe.

Pentru a ıntelege mai bine strategia utilizata de MEF vom re-curge din nou la exemple. Vom aborda mai ıntai problema ex-pandarii aproximantelor pentru cazul marimilor de tip scalar. Intr-oetapa urmatoare, vom generaliza procedura de expandare la cazulnecunoscutelor vectoriale.

3.4.1 Expandarea aproximantelor unei marimide tip scalar

Vom detalia procedura de expandare pe cazul unui domeniu planpentagonal discretizat ın cinci elemente triunghiulare liniare (fig.3.17). Matricea de conexiuni a retelei este prezentata ın tabelul 3.1.

Page 98: Metoda Elementelor Finite

92 Capitolul 3

x

y

1

2

4

5

63

(1)

(2)

(3)(4)

(5)

Figura 3.17: Domeniu bidimensional discretizat ın cinci elementetriunghiulare liniare

Tabelul 3.1: Matricea de conexiuni a retelei de elemente finite dinfigura 3.17

Element Nod Nod Nod(e) 1 2 3

(1) 2 3 1(2) 3 2 4(3) 5 3 4(4) 6 3 5(5) 3 6 1

Page 99: Metoda Elementelor Finite

Expandarea aproximantelor elementale 93

Asocierea dintre noduri si elemente a fost facuta respectand regulaetichetarii ın sens antiorar (vezi §3.2.1). Notam cu ϕ necunoscutascalara care trebuie aproximata la nivelul fiecarui element. Valorileei ın cele sase noduri ale retelei vor fi desemnate prin simbolurileϕk (k = 1, 2, . . . , 6).

Aproximantele elementale ale marimii ϕ rezulta prin particula-rizarea ecuatiei (3.17) cu ajutorul informatiei din tabelul 3.1:

ϕ(1) = ϕ2N(1)2 + ϕ3N

(1)3 + ϕ1N

(1)1 ,

ϕ(2) = ϕ3N(2)3 + ϕ2N

(2)2 + ϕ4N

(2)4 ,

ϕ(3) = ϕ5N(3)5 + ϕ3N

(3)3 + ϕ4N

(3)4 ,

ϕ(4) = ϕ6N(4)6 + ϕ3N

(4)3 + ϕ5N

(4)5 ,

ϕ(5) = ϕ3N(5)3 + ϕ6N

(5)6 + ϕ1N

(5)1 .

(3.48)

Pentru a nu aglomera excesiv notatia, la scrierea relatiilor (3.48)am omis precizarea variabilelor independente x si y. Fiecare dintreaproximantele (3.48) este valabila la nivelul unui element. Pentru ale putea distinge, le-am atasat cate un indice superior ıncadrat deparanteze rotunde.

Din relatiile (3.48) ne putem da seama ın ce consta pro-blema unificarii expresiei aproximantelor elementale. Privind acesteecuatii, observam imediat ca ın fiecare din ele apare alt set de valorinodale ale necunoscutei ϕ (de exemplu, ın prima intervin ϕ1, ϕ2 siϕ3, iar ın a patra intervin ϕ3, ϕ5 si ϕ6). Intr-o asemenea forma,aproximantele sunt greu manevrabile ın etapa de construire teo-retica a modelului cu elemente finite. Sa ne gandim numai cat dedificil de surprins ıntr-o formula este banala operatie de adunare afunctiilor (3.48).

Ideala ar fi deci unificarea expresiei aproximantelor (3.48), ınasa fel ıncat ele sa contina acelasi set de necunoscute nodale. Putemrealiza acest lucru adaugand la fiecare din relatiile (3.48) termeniobtinuti prin ınmultirea valorilor nodale absente cu zero. In felulacesta, aproximarea dobandeste un plus de generalitate. Aplicand

Page 100: Metoda Elementelor Finite

94 Capitolul 3

strategia de mai sus, obtinem urmatorul set de aproximante ele-mentale (termenii au fost ordonati crescator dupa numarul noduluiasociat):

ϕ(1) = ϕ1N(1)1 + ϕ2N

(1)2 + ϕ3N

(1)3 + ϕ4 · 0 + ϕ5 · 0 + ϕ6 · 0,

ϕ(2) = ϕ1 · 0 + ϕ2N(2)2 + ϕ3N

(2)3 + ϕ4N

(2)4 + ϕ5 · 0 + ϕ6 · 0,

ϕ(3) = ϕ1 · 0 + ϕ2 · 0 + ϕ3N(3)3 + ϕ4N

(3)4 + ϕ5N

(3)5 + ϕ6 · 0,

ϕ(4) = ϕ1 · 0 + ϕ2 · 0 + ϕ3N(4)3 + ϕ4 · 0 + ϕ5N

(4)5 + ϕ6N

(4)6 ,

ϕ(5) = ϕ1N(5)1 + ϕ2 · 0 + ϕ3N

(5)3 + ϕ4 · 0 + ϕ5 · 0 + ϕ6N

(5)6 .

(3.49)

Formulele (3.49) definesc asa-numita forma expandata a apro-

ximantelor elementale ale necunoscutei ϕ. Se observa generalitatealor. Intr-adevar, putem sa le reprezentam pe toate prin urmatoareaecuatie matriceala:

ϕ(e) =N (e)

ϕ , e = 1, 2, . . . , 5, (3.50)

ın careϕ = [ϕ1, ϕ2, ϕ3, ϕ4, ϕ5]

T (3.51)

este un vector-coloana ce grupeaza toate necunoscutele nodale, iarN (e)

este un vector-linie care contine, ın anumite pozitii, functiile

de forma asociate nodurilor elementului (e), ın rest el fiind plin cuzerouri. De fapt, pozitiile pe care se afla coeficientii nenuli ın

N (e)

se deduc usor din matricea de conexiuni (tab. 3.1). De exemplu,vectorul-linie

N (5)

are urmatoarea structura:

N (5)=hN

(5)1 , 0, N

(5)3 , 0, 0, N

(5)6

i. (3.52)

Manevrarea aproximantelor elementale sub forma expandata(3.50) este mult mai usoara, fiindca acum vectorul-coloana ϕactioneaza ca un factor comun pentru toate elementele din retea.Chiar si programele MEF opereaza conceptual tot cu forma expan-data (3.50), desi ele nu stocheaza vectorul-linie

N (e)

ın ıntregime.

Page 101: Metoda Elementelor Finite

Expandarea aproximantelor elementale 95

Ar fi neeconomic sa se pastreze ın memoria calculatorului o struc-tura de date plina aproape peste tot cu valori nule. In realitate, pro-gramele MEF retin doar coeficientii nenuli din

N (e)

si matricea de

conexiuni (ea furnizeaza pozitiile ın care se afla acesti coeficienti).

3.4.2 Expandarea aproximantelor unei marimide tip vectorial

Pentru a ıntelege mai usor modul ın care se generalizeaza strate-gia de expandare atunci avem de-a face cu necunoscute de tip vec-torial, vom considera ca exemplu tot domeniul plan din figura 3.17.Marimea care trebuie acum aproximata la nivelul fiecarui elementeste un camp vectorial ~d cu doua componente (notate u, respectivv). Vom desemna valorile nodale ale acestor componente prin sim-bolurile uk si vk (k = 1, 2, . . . , 6). Evident, matricea de conexiuni aretelei ramane cea din tabelul 3.1.

Aproximantele elementale ale componentelor lui ~d rezulta parti-cularizand formulele (3.46) cu ajutorul informatiei din tabelul 3.1:

u(1) = u2N(1)2 + u3N

(1)3 + u1N

(1)1 ,

v(1) = v2N(1)2 + v3N

(1)3 + v1N

(1)1 ,

u(2) = u3N(2)3 + u2N

(2)2 + u4N

(2)4 ,

v(2) = v3N(2)3 + v2N

(2)2 + v4N

(2)4 ,

u(3) = u5N(3)5 + u3N

(3)3 + u4N

(3)4 ,

v(3) = v5N(3)5 + v3N

(3)3 + v4N

(3)4 ,

u(4) = u6N(4)6 + u3N

(4)3 + u5N

(4)5 ,

v(4) = v6N(4)6 + v3N

(4)3 + v5N

(4)5 ,

u(5) = u3N(5)3 + u6N

(5)6 + u1N

(5)1 ,

v(5) = v3N(5)3 + v6N

(5)6 + v1N

(5)1 .

(3.53)

Page 102: Metoda Elementelor Finite

96 Capitolul 3

Din nou, pentru a putea face distinctie ıntre aproximante, le-amatasat un indice superior ıncadrat de paranteze rotunde.

Pentru unificarea relatiilor (3.53) procedam ıntr-o manieraasemanatoare cu cea pe care am utilizat-o ın paragraful anterior:

u(1) = u1N(1)1 + v1 · 0 + u2N

(1)2 + v2 · 0 + u3N

(1)3 + v3 · 0+

+ u4 · 0 + v4 · 0 + u5 · 0 + v5 · 0 + u6 · 0 + v6 · 0,v(1) = u1 · 0 + v1N

(1)1 + u2 · 0 + v2N

(1)2 + u3 · 0 + v3N

(1)3 +

+ u4 · 0 + v4 · 0 + u5 · 0 + v5 · 0 + u6 · 0 + v6 · 0,u(2) = u1 · 0 + v1 · 0 + u2N

(2)2 + v2 · 0 + u3N

(2)3 + v3 · 0+

+ u4N(2)4 + v4 · 0 + u5 · 0 + v5 · 0 + u6 · 0 + v6 · 0,

v(2) = u1 · 0 + v1 · 0 + u2 · 0 + v2N(2)2 + u3 · 0 + v3N

(2)3 +

+ u4 · 0 + v4N(2)4 + u5 · 0 + v5 · 0 + u6 · 0 + v6 · 0,

u(3) = u1 · 0 + v1 · 0 + u2 · 0 + v2 · 0 + u3N(3)3 + v3 · 0+

+ u4N(3)4 + v4 · 0 + u5N

(3)5 + v5 · 0 + u6 · 0 + v6 · 0,

v(3) = u1 · 0 + v1 · 0 + u2 · 0 + v2 · 0 + u3 · 0 + v3N(3)3 +

+ u4 · 0 + v4N(3)4 + u5 · 0 + v5N

(3)5 + u6 · 0 + v6 · 0,

u(4) = u1 · 0 + v1 · 0 + u2 · 0 + v2 · 0 + u3N(4)3 + v3 · 0+

+ u4 · 0 + v4 · 0 + u5N(4)5 + v5 · 0 + u6N

(4)6 + v6 · 0,

v(4) = u1 · 0 + v1 · 0 + u2 · 0 + v2 · 0 + u3 · 0 + v3N(4)3 +

+ u4 · 0 + v4 · 0 + u5 · 0 + v5N(4)5 + u6 · 0 + v6N

(4)6 ,

u(5) = u1N(5)1 + v1 · 0 + u2 · 0 + v2 · 0 + u3N

(5)3 + v3 · 0+

+ u4 · 0 + v4 · 0 + u5 · 0 + v5 · 0 + u6N(5)6 + v6 · 0,

v(5) = u1 · 0 + v1N(5)1 + u2 · 0 + v2 · 0 + u3 · 0 + v3N

(5)3 +

+ u4 · 0 + v4 · 0 + u5 · 0 + v5 · 0 + u6 · 0 + v6N(5)6 .

(3.54)

Relatiile (3.54) definesc forma expandata a aproximantelor elemen-

tale ale marimii ~d. Ele au expresia matriceala

Page 103: Metoda Elementelor Finite

Expandarea aproximantelor elementale 97nd(e)

o=N (e)

d , e = 1, 2, . . . , 5, (3.55)

ın care nd(e)

o=u(e), v(e)

T, e = 1, 2, . . . , 5, (3.56)

este un vector-coloana ce grupeaza aproximantele componentelor usi v la nivelul elementului (e),

d = [u1, v1, u2, v2, u3, v3, u4, v4, u5, v5, u6, v6]T (3.57)

este un vector-coloana ın care sunt stocate ın ordine valorile nodaleale componentelor u si v, iar

N (e)

este o matrice cu doua linii

si un numar de coloane egal cu dublul numarului total de noduri(2×6 noduri = 12 coloane, pentru exemplul din figura 3.17).

N (e)

contine, ın anumite pozitii, functiile de forma asociate nodurilorelementului (e), peste tot ın rest avand zerouri. Din nou, pozitiileın care se afla coeficientii nenuli pot fi deduse din matricea de co-nexiuni (tab. 3.1). De exemplu,

N (5)

are urmatoarea structura:

N (5)2×12

=

=h

N(5)1

i2×2

, [0]2×2 ,hN

(5)3

i2×2

, [0]2×2 , [0]2×2 ,hN

(5)6

i2×2

,

(3.58)

unde [0]2×2 este matricea nula de ordinul doi, iarhN

(5)1

i2×2

,hN

(5)3

i2×2

sihN

(5)6

i2×2

sunt la randul lor matrici diagonale definite

astfel: hN

(5)1

i2×2

=

"N

(5)1 0

0 N(5)1

#,h

N(5)3

i2×2

=

"N

(5)3 0

0 N(5)3

#,h

N(5)6

i2×2

=

"N

(5)6 0

0 N(5)6

#.

(3.59)

Procedura de expandare pe care am aplicat-o mai sus ın cazulunei marimi vectoriale cu numai doua componente poate fi genera-lizata cu usurinta pentru marimile vectoriale cu trei componente.

Page 104: Metoda Elementelor Finite

98 Capitolul 3

3.5 Concluzii

Am prezentat ın acest capitol cateva familii de elemente adec-vate rezolvarii problemelor de inginerie mecanica. In afara discutieinoastre au ramas o serie de elemente a caror structura este maicomplexa (elementele de tip grinda, membrana, placa, respectivınvelitoare). Totusi, informatia prezentata este suficienta pentru afurniza o imagine generala asupra constructiei elementelor finite.

Am vazut astfel ca un element este o entitate complexa re-prezentata prin trei caracteristici: domeniul spatial ocupat, numarulsi pozitia nodurilor asociate, respectiv polinomul de aproximareutilizat. Numarul nodurilor si gradul polinomului aproximant suntıntr-o relatie precisa, fiindca procedura de evaluare a coeficientilorpolinomiali face apel la valorile nodale ale necunoscutei.

In general, calitatea aproximarii este cu atat mai buna cu catgradul polinomului este mai mare. Am aratat ınsa ca numarul no-durilor creste repede daca marim gradul polinomului. Va trebui decisa cautam un compromis ıntre precizia solutiei si efortul de calculnecesar pentru obtinerea ei. Regula tocmai enuntata ne conduce laurmatoarele concluzii:

• pentru rezolvarea problemelor bidimensionale este recoman-dabila folosirea elementului patrulater biliniar ın locul celuitriunghiular liniar (primul, avand un polinom aproximant degrad superior, va asigura obtinerea unei solutii mai precise,fara sa mareasca exagerat numarul nodurilor retelei);

• ın mod asemanator, ın cazul problemelor tridimensionale, ele-mentul hexaedric triliniar este mai precis decat elementultetraedric liniar (aceasta ımbunatatire a performantelor im-plica ınsa un numar dublu de noduri la nivelul elementului).

Vom retine deci ca este bine sa evitam pe cat posibil folosirea unorelemente cu polinoame de aproximare de gradul ıntai. Totusi, atuncicand ne decidem asupra tipului de element trebuie sa estimamvolumul de calcule necesar pentru obtinerea solutiei numerice.

Page 105: Metoda Elementelor Finite

Capitolul 4

Modelul cu elementefinite al problemelorde elasticitate liniara

In acest capitol vom vedea cum se aplica metoda elementelorfinite (MEF) la rezolvarea numerica a problemelor de elasticitateliniara. Intr-o prima etapa vom prezenta sistemul de ecuatii cu de-rivate partiale care descrie deformarea unui corp elastic. Vom aratacare sunt necunoscutele problemei si vom insista asupra ierarhizariilor. De asemenea, vom analiza motivele care ımpiedica rezolvareaanalitica a sistemului de ecuatii cu derivate partiale ın marea ma-joritate a aplicatiilor practice. In continuare, vom arata cum seconstruieste modelul cu elemente finite al problemelor de elastici-tate liniara pornind de la teorema de minim a energiei potentiale.Pe parcursul expunerii vom face referiri si la procedeele numericefolosite pentru construirea si rezolvarea sistemului de ecuatii glo-bal. De asemenea, vom arata cum se prelucreaza solutia sistemuluide ecuatii global ın vederea determinarii deformatiilor si tensiuni-lor. Un paragraf distinct va fi consacrat prezentarii unor problemeparticulare: starea plana de tensiuni, starea plana de deformatii,respectiv probleme cu simetrie axiala.

99

Page 106: Metoda Elementelor Finite

100 Capitolul 4

4.1 Formularea problemelor de

elasticitate liniara

Ne propunem spre rezolvare urmatoarea problema generala:

Sa se determine starea de tensiuni si deformatii a unuicorp liniar elastic rezultata ca urmare a solicitarii salede catre un sistem de forte exterioare.

Inainte de a trece la exprimarea problemei ın limbaj matematic,vom prezenta cateva notatii si ipoteze pe care le vom folosi frecventın continuare.

Vom nota cu Ω domeniul spatial ocupat de corp ın configuratianedeformata (fig. 4.1). Vom presupune ca frontiera Σ a domeniuluiΩ este o suprafata ınchisa si neteda pe portiuni. Pentru exprimareapozitiilor particulelor corpului vom folosi un sistem de coordonatecarteziene ortogonale (x, y, z).

Vom admite ca asupra corpului actioneaza numai ıncarcari denatura mecanica. Vom ımparti aceste ıncarcari ın doua categorii:

• forte volumice (numite astfel pentru ca actioneaza asupra tu-turor particulelor corpului);

• forte de suprafata (numite astfel pentru ca actioneaza asupraparticulelor situate ın anumite portiuni ale frontierei Σ).

Greutatea proprie este un exemplu tipic de forta volumica. Inmulte situatii, mai ales cand este vorba de structuri cu gabaritredus, deformatiile si tensiunile datorate propriei greutati pot fineglijate. Sunt ınsa si cazuri (de exemplu, constructii foarte ma-sive de tipul podurilor sau al turnurilor) cand efectele mecaniceale greutatii nu mai pot fi catusi de putin trecute cu vederea, iarintroducerea lor ın model devine absolut necesara.

Fortele de suprafata sunt asociate de regula interactiunilor decontact cu alte corpuri (aflate ın stare solida, lichida sau gazoasa).

Page 107: Metoda Elementelor Finite

Formularea problemelor de elasticitate liniara 101

(Σ)

(Ω)

P

P′

~d

~s

~s

~n

y

z

x

Figura 4.1: Notatii utilizate ın modelul unui corp liniar elastic

Un exemplu de asemenea forte sunt cele exercitate de un mediulichid asupra peretelui interior al unui cilindru hidraulic.

Din punct de vedere matematic, fortele volumice sunt descriseprin densitatea lor spatiala ~b, iar fortele de suprafata – prin densi-tatea lor superficiala ~s. Vom nota cu Σs portiunea de pe frontieraΣ pe care actioneaza forte de suprafata (echivalentul acestor forteeste reprezentat de ansamblul conditiilor la limita de tip ıncarcare

superficiala).

Page 108: Metoda Elementelor Finite

102 Capitolul 4

Corpul studiat poate fi supus si unor restrictii de miscare (deexemplu, prin prezenta unor reazeme de diferite tipuri). In gene-ral, restrictiile de miscare se exprima sub forma unor conditii im-puse deplasarilor particulelor de pe anumite portiuni ale frontiereiΣ. Vom nota cu Σd zona de pe frontiera pe care sunt specificaterestrictii de miscare (acestea se mai numesc conditii la limita de tip

deplasare).

Σd si Σs trebuie sa se afle ın urmatoarea relatie:

Σd ∪ Σs = Σ, Σd ∩ Σs = ∅. (4.1)

Ea arata ca ın fiecare punct de pe frontiera corpului trebuie pre-cizata fie o ıncarcare exterioara (chiar daca aceasta este nula), fieo restrictie de miscare, ınsa ın nici un caz ambele deodata. Dacane gandim putin, ne dam seama ca regula de mai sus este foartenaturala. Intr-adevar, acolo unde sunt impuse deplasarile particu-lelor (din cauza unei rezemari, de exemplu) ne intereseaza calcu-lul reactiunilor asociate. Invers, acolo unde sunt impuse anumiteıncarcari ne intereseaza determinarea deformatiilor rezultante. Nici-odata nu este corect sa introducem simultan restrictii de miscare siıncarcari ıntr-un acelasi punct, fiindca putem veni ın contradictiecu ecuatiile care descriu comportarea mecanica a corpului.

Sub actiunea fortelor exterioare, particulele corpului sufera de-plasari. Aceste deplasari sunt descrise de un camp vectorial ~d careleaga pozitia initiala a fiecarei particule de pozitia ei finala (fig.

4.1). Vom nota cu u, v si w componentele campului vectorial ~dın raport cu axele x, y, respectiv z. In tot ceea ce urmeaza vomadmite ca deplasarile si deformatiile corpului sunt foarte mici. Inaceste conditii, vom putea scrie ecuatiile de echilibru pe domeniulΩ si frontiera Σ, fara a comite erori semnificative. Daca ar fi saoperam cu toata rigoarea, ar trebui sa exprimam respectivele ecuatiiın configuratia finala a corpului, ınsa o asemenea abordare ar danastere la complicatii matematice serioase si nu s-ar justifica ıncazul unor deformatii foarte mici.

Page 109: Metoda Elementelor Finite

Formularea problemelor de elasticitate liniara 103

Pentru descrierea distorsiunii corpului vom folosi o marime ten-soriala numita deformatie specifica (sau, pur si simplu, deformatie).Ea se exprima prin intermediul a sase componente distincte [1]:

εx =∂ u

∂ x,

εy =∂ v

∂ y,

εz =∂ w

∂ z,

γxy =∂ u

∂ y+

∂ v

∂ x,

γyz =∂ v

∂ z+

∂ w

∂ y,

γzx =∂ w

∂ x+

∂ u

∂ z·

(4.2)

εx, εy si εz se numesc deformatii liniare, iar γxy, γyz si γzx se numescdeformatii unghiulare.

Observam ca relatiile (4.2) exprima deformatiile ca derivate

spatiale ale componentelor campului vectorial ~d. Sub aspect fizic,ele descriu discrepantele dintre deplasarile particulelor ınvecinate.Mai concret, daca doua particule foarte apropiate au tendinta dea se misca la fel, deformatia din zona respectiva a corpului va finula. In cazul ın care ele tind sa se miste diferit, corpul va suferidistorsiuni, iar deformatiile date de relatiile (4.2) vor fi nenule.

In general, ın masa corpului apar solicitari interne ca urmarea deformarii. Pentru descrierea acestora vom folosi o alta marimetensoriala numita tensiune. Vom nota componentele distincte aletensiunii cu σx, σy, σz, τxy, τyz si τzx. Marimile σx, σy si σz senumesc tensiuni normale, iar τxy, τyz si τzx – tensiuni tangentiale.

Se poate demonstra ca, ın orice punct al corpului, componenteleenumerate mai sus trebuie sa verifice urmatorul sistem de ecuatiicu derivate partiale [1]:

Page 110: Metoda Elementelor Finite

104 Capitolul 4

∂ σx

∂ x+

∂ τxy

∂ y+

∂ τzx

∂ z+ bx = 0,

∂ τxy

∂ x+

∂ σy

∂ y+

∂ τyz

∂ z+ by = 0,

∂ τzx

∂ x+

∂ τyz

∂ y+

∂ σz

∂ z+ bz = 0, (x, y, z) ∈ Ω,

(4.3)

ın care bx, by si bz sunt componentele densitatii spatiale ~b a fortelorvolumice. Daca neglijam efectul mecanic al ıncarcarilor volumice,termenii bx, by si bz se anuleaza.

Sistemul (4.3) exprima echilibrul fortelor care actioneaza asupravolumelor elementare ale corpului. El este cunoscut ın literatura despecialitate sub denumirea de sistemul ecuatiilor de echilibru ale lui

Cauchy.Tensiunile care apar ın masa corpului sunt legate de distorsiunile

acestuia prin asa-numitele ecuatii constitutive. Vom admite ca ma-terialul din care este realizat corpul are o comportare liniar elasticasi, ın plus, este izotrop. In aceste conditii, ecuatiile constitutive potfi scrise sub forma data de legea lui Hooke [1]:

σx =E(1 − ν)

(1 + ν)(1 − 2ν)εx +

(1 + ν)(1 − 2ν)(εy + εz) ,

σy =E(1 − ν)

(1 + ν)(1 − 2ν)εy +

(1 + ν)(1 − 2ν)(εz + εx) ,

σz =E(1 − ν)

(1 + ν)(1 − 2ν)εz +

(1 + ν)(1 − 2ν)(εx + εy) ,

τxy =E

2(1 + ν)γxy,

τyz =E

2(1 + ν)γyz,

τzx =E

2(1 + ν)γzx.

(4.4)

Page 111: Metoda Elementelor Finite

Formularea problemelor de elasticitate liniara 105

E si ν din relatiile (4.4) sunt modulul lui Young (sau modulul de

elasticitate longitudinala), respectiv coeficientul lui Poisson (saucoeficientul de contractie transversala). Aceste doua marimi suntconstante de material.

Deformarea unui corp elastic este descrisa de ansamblul ecu-atiilor (4.2), (4.3) si (4.4). Ele formeaza un sistem de cincisprezeceecuatii cu derivate partiale definite pe domeniul Ω. Necunoscutelelui sunt:

• componentele u, v si w ale campului vectorilor deplasare;

• deformatiile εx, εy, εz, γxy, γyz si γzx;

• tensiunile σx, σy, σz, τxy, τyz si τzx.

Acestui sistem de ecuatii i se ataseaza doua tipuri de conditii lalimita:

• conditii de tip deplasare

u(x, y, z) = u(x, y, z),

v(x, y, z) = v(x, y, z),

w(x, y, z) = w(x, y, z), (x, y, z) ∈ Σd,

(4.5)

ın care u, v si w sunt componentele unui camp de vectorideplasare cunoscut pe portiunea Σd a frontierei;

• conditii de tip ıncarcare superficiala

σxnx + τxyny + τzxnz = sx(x, y, z),

τxynx + σyny + τyznz = sy(x, y, z),

τzxnx + τyzny + σznz = sz(x, y, z), (x, y, z) ∈ Σs,

(4.6)

ın care sx, sy si sz sunt componentele campului ıncarcarilorsuperficiale care actioneaza pe portiunea Σs a frontierei, iarnx, ny si nz sunt componentele versorului normalei exterioare

la frontiera (fig. 4.1). In general, sx, sy si sz sunt marimi cu-noscute. Daca unele dintre ele (sau toate) sunt nule, membruldrept corespunzator din (4.6) devine pur si simplu zero.

Page 112: Metoda Elementelor Finite

106 Capitolul 4

Sistemul de ecuatii cu derivate partiale prezentat mai sus poatefi redus ın asa fel ıncat sa nu mai contina decat trei necunoscute.Intr-adevar, daca ınlocuim expresiile (4.2) ale deformatiilor ın legealui Hooke (4.4), obtinem relatii directe ıntre tensiuni si componen-tele campului vectorilor deplasare. Inlocuind mai departe acesterelatii ın ecuatiile de echilibru (4.3), ajungem la un sistem de treiecuatii cu derivate partiale de ordinul al doilea avand ca necu-noscute numai deplasarile u, v si w ale particulelor. Conditiile lalimita (4.5) pot fi pastrate asa cum sunt, ıntrucat implica tocmainecunoscutele ramase (componente ale campului vectorilor depla-sare). In schimb, asupra conditiilor (4.6) va trebui sa operam otransformare similara celei la care a fost supus sistemul de ecuatii,pentru a le exprima numai ın termenii deplasarilor.

S-ar parea ca problema s-a simplificat mult. In loc sa fim obligatila a determina cincisprezece functii necunoscute definite pe dome-niul Ω, am ramas cu numai trei si anume, componentele campuluivectorilor deplasare. Evident, odata ce vom fi obtinut functiile u,v si w, vom putea sa le ınlocuim ın relatiile (4.2) pentru a cal-cula deformatiile εx, εy, εz, γxy, γyz si γzx. Mai departe, ınlocuinddeformatiile ın legea lui Hooke (4.4), putem calcula tensiunile σx,σy, σz, τxy, τyz si τzx. Aceasta ordine defineste ierarhia necunos-cutelor unei probleme de elasticitate. Pe prima pozitie ın ierarhiesunt deplasarile (ca necunoscute primare). Pe nivelul urmator seafla deformatiile, iar pe ultima pozitie sunt tensiunile (deformatiilesi tensiunile sunt incluse ın categoria asa-numitelor necunoscutederivate).

Chiar si ın forma sa redusa, sistemul de ecuatii care descriedeformarea corpurilor elastice poate fi rezolvat exact numai pen-tru cateva cazuri relativ simple (mult mai simple decat majori-tatea problemelor ıntalnite ın practica inginereasca). Principalelemotive care ımpiedica obtinerea unei solutii analitice a probleme-lor de elasticitate sunt complexitatea formei domeniului spatial Ω,respectiv complexitatea conditiilor la limita de tipul deplasarilorsi/sau ıncarcarilor superficiale.

Page 113: Metoda Elementelor Finite

Modelul cu elemente finite al problemelor de elasticitate 107

In situatiile ın care nu pot determina o solutie exacta, ingineriiapeleaza de regula la un procedeu de rezolvare aproximativa a pro-blemei cu care se confrunta. Vom vedea ın continuare cum poatefi utilizata metoda elementelor finite pentru obtinerea unei solutiinumerice a problemelor de elasticitate.

4.2 Construirea modelului cu elemente

finite al problemelor de elasticitate

liniara utilizand teorema de minim

a energiei potentiale

Am aratat ın capitolul 2 ca MEF cauta solutia aproximativaa unei probleme sub forma unei functii care se apropie cat maibine de solutia exacta pe ıntregul domeniu de analiza. In acestscop, este folosit un criteriu de minimizare globala a erorii construitca o marime de tip integral. Am mentionat ın respectivul capitolca, ın cazul deformarii unui corp elastic, criteriul poate fi obtinutpornind de la teorema de minim a energiei potentiale. Am aratatchiar maniera practica ın care se construieste un model cu elementefinite bazat pe aceasta teorema, folosind ca exemplu o problemafoarte simpla (deformarea unei corzi perfect flexibile sub actiuneapropriei greutati).

Nu vom relua aici enuntul teoremei de minim a energieipotentiale elastice (el poate fi gasit ın §2.3). Vom reaminti doar fap-tul ca aceasta marime este diferenta dintre energia de deformatieΠdef si lucrul mecanic al ıncarcarilor exterioare Wext (vezi relatia

(2.19)). In cazul general al unui corp liniar elastic, energia de defor-matie are expresia [4, 5, 8, 9, 10]

Πdef =ZΩ

1

2(σxεx + σyεy + σzεz+

+τxyγxy + τyzγyz + τzxγzx) dΩ.

(4.7)

Page 114: Metoda Elementelor Finite

108 Capitolul 4

In ceea ce priveste lucrul mecanic al ıncarcarilor exterioare, el poatefi scris ca o suma a lucrurilor mecanice efectuate de diferitele cate-gorii de forte care actioneaza asupra corpului:

Wext = Wb + Ws + Wp. (4.8)

Marimile din relatia de mai sus au urmatoarea semnificatie: Wb

este lucrul mecanic al fortelor volumice, Ws este lucrul mecanic alfortelor de suprafata, iar Wp este lucrul mecanic al fortelor punc-tuale (sau concentrate). Primele doua categorii de ıncarcari le-amıntalnit deja ın §4.1. In schimb, fortele punctuale nu au fost incluseın sistemul de ecuatii care descrie deformarea corpurilor elastice.Vom vedea ceva mai jos de unde provine incapacitatea modeluluidiferential de a manevra astfel de ıncarcari. Insa ınainte de aceastane vom ocupa de expresiile lui Wb si Ws.

Cu ajutorul componentelor bx, by si bz ale densitatii spatiale aıncarcarilor volumice, ıl putem scrie pe Wb sub forma

Wb =ZΩ

(bxu + byv + bzw) dΩ. (4.9)

In mod asemanator, pentru Ws putem scrie relatia

Ws =ZΣs

(sxu + syv + szw) dΣ, (4.10)

ın care sx, sy si sz sunt componentele densitatii superficiale aıncarcarilor care actioneaza pe portiunea Σs a frontierei corpului.

Sa vedem acum de ce lipsesc fortele punctuale din sistemulde ecuatii prezentat ın §4.1. Motivul absentei lor rezida ın faptulca operatorii diferentiali nu pot surprinde efectul unor ıncarcari acaror zona de aplicare se reduce la un punct (ın realitate, forteleactioneaza ıntotdeauna pe suprafete sau volume finite, ınsa uneoriacestea sunt atat de restranse ıncat pot fi considerate punctiforme).Includerea unor forte concentrate ın sistemul de ecuatii cu deri-vate partiale conduce la aparitia unor singularitati care distrug pro-prietatile de continuitate si derivabilitate ale functiilor din model

Page 115: Metoda Elementelor Finite

Modelul cu elemente finite al problemelor de elasticitate 109

(propriu-zis, sistemul de ecuatii ar deveni neviabil, fiindca este in-corect sa aplicam un operator de derivare unei functii care nu estediferentiabila ın toate punctele domeniului Ω).

Ne putem ıntreba atunci pe buna dreptate: daca ıncarcarilepunctuale nu pot fi manevrate prin sistemul de ecuatii cu deri-vate partiale, cum de pot aparea totusi ıntr-un model cu elementefinite? Obtinem raspunsul la aceasta ıntrebare daca sesizam o carac-teristica esentiala a modelelor construite pornind de la teorema deminim a energiei potentiale elastice. Relatiile (4.7), (4.9) si (4.10),precum si aspectele discutate ın capitolul 2 ne arata ca energiapotentiala elastica este o marime de tip integral. In consecinta,modelul cu elemente finite asociat va fi unul global, ın sensul caopereaza la nivelul unui domeniu si nu ın maniera localizata, asacum se ıntampla cu modelul diferential. Singularitatile introduse defortele punctuale vor fi atenuate prin efectul de mediere al integra-lelor. Concluzionam ca modelul cu elemente finite este mai tolerantdecat cel diferential si, din acest punct de vedere, mai flexibil. Defapt, putem afirma ca este chiar mai general, fiindca accepta o gamamai larga de ıncarcari exterioare.

In cele ce urmeaza vom presupune ca asupra corpului actioneazaun sistem de np forte punctuale. Vom nota aceste forte cu ~Pk (k =1, 2, . . . , np). Lucrul mecanic efectuat de ele va avea atunci expresia

Wp =npX

k=1

(Px,kuk + Py,kvk + Pz,kwk) , (4.11)

ın care Px,k, Py,k si Pz,k sunt componentele fortei ~Pk (k =1, 2, . . . , np), iar uk, vk si wk sunt deplasarile punctului ei de aplicatieın lungul axelor de coordonate.

Relatiile (2.19) si (4.8) ne permit sa rescriem energia potentialaelastica sub urmatoarea forma generala:

Π = Πdef − Wb − Ws − Wp. (4.12)

Termenii din membrul drept al lui (4.12) sunt definiti de (4.7), (4.9),(4.10) si (4.11).

Page 116: Metoda Elementelor Finite

110 Capitolul 4

Acestea fiind spuse, suntem pregatiti sa trecem la construireamodelului cu elemente finite al problemelor de elasticitate.

Presupunem ca domeniul Ω a fost discretizat ın nelem elementefinite spatiale de un tip oarecare. Notam cu nnod numarul totalal nodurilor din retea. Admitem ca discretizarea s-a facut ın asafel ıncat ın punctul de aplicare al fiecarei forte concentrate ~Pk

(k = 1, 2, . . . , np) exista un nod. Vom vedea ca respectarea acesteiconditii are drept consecinta simplificarea expresiei lucrului meca-nic Wp.

O etapa esentiala ın aplicarea MEF pentru rezolvarea numericaa unei probleme este construirea aproximantelor elementale ale ne-cunoscutelor. In cazul problemelor de elasticitate, marimile caretrebuie aproximate sunt:

• componentele u, v si w ale campului vectorilor deplasare;

• deformatiile εx, εy, εz, γxy, γyz si γzx;

• tensiunile σx, σy, σz, τxy, τyz si τzx.

Incepem prin a ne ocupa de necunoscutele primare u, v si w.Ele fiind componentele unei marimi de tip vectorial, le vom apro-xima prin relatii de tipul (3.47). Putem scrie respectivele ecuatii suburmatoarea forma matriceala expandata (vezi discutia din §3.4):n

d(e)o

=N (e)

d , e = 1, 2, . . . , nelem, (4.13)

ın care nd(e)

o=u(e), v(e), w(e)

T(4.14)

este un vector-coloana ce grupeaza aproximantele marimilor u, v siw la nivelul unui element oarecare (e),

d = [u1, v1, w1, u2, v2, w2, . . . , unnod, vnnod

, wnnod]T (4.15)

este un vector-coloana ın care sunt stocate ın ordine valorile no-dale ale lui u, v si w, iar

N (e)

este o matrice cu trei linii si un

Page 117: Metoda Elementelor Finite

Modelul cu elemente finite al problemelor de elasticitate 111

numar de coloane egal cu triplul numarului total de noduri.N (e)

contine, ın anumite pozitii, functiile de forma asociate nodurilorelementului (e), peste tot ın rest avand zerouri. Pozitiile ın care seafla coeficientii nenuli rezulta din matricea de conexiuni a retelei.

Avand construite aproximantele u(e), v(e) si w(e) ale componen-telor campului vectorilor deplasare, le putem introduce ın relatiile(4.2). Obtinem astfel aproximari elementale ale deformatiilor εx, εy,εz, γxy, γyz si γzx:

ε(e)x =

∂ u(e)

∂ x,

ε(e)y =

∂ v(e)

∂ y,

ε(e)z =

∂ w(e)

∂ z,

γ(e)xy =

∂ u(e)

∂ y+

∂ v(e)

∂ x,

γ(e)yz =

∂ v(e)

∂ z+

∂ w(e)

∂ y,

γ(e)zx =

∂ w(e)

∂ x+

∂ u(e)

∂ z·

(4.16)

Daca tinem cont de (4.14), putem rescrie ecuatiile (4.16) suburmatoarea forma matriceala:

ε(e)©

= [g]nd(e)

o, e = 1, 2, . . . , nelem, (4.17)

ın care ε(e)

©=ε(e)

x , ε(e)y , ε(e)

z , γ(e)xy , γ(e)

yz , γ(e)zx

T(4.18)

este un vector-coloana ce grupeaza aproximantele elementale aledeformatiilor, iar

Page 118: Metoda Elementelor Finite

112 Capitolul 4

[g] =

266666666666664∂∂ x

0 0

0 ∂∂ y

0

0 0 ∂∂ z

∂∂ y

∂∂ x

0

0 ∂∂ z

∂∂ y

∂∂ z

0 ∂∂ x

377777777777775 (4.19)

este asa-numita matrice a gradientilor. Liniile ei contin operatori dederivare care actioneaza asupra componentelor vectorului-coloanand(e)

o. (4.13) si (4.17) permit explicitarea aproximantelor elemen-

tale ale deformatiilor ca dependente de deplasarile nodale:ε(e)

©=B(e)

d , e = 1, 2, . . . , nelem, (4.20)

unde B(e)

= [g]

N (e)

(4.21)

este o matrice cu sase linii si un numar de coloane egal cu triplulnumarului total de noduri. Ea contine, ın anumite pozitii, derivatespatiale ale functiilor de forma ale elementului (e), peste tot ın restavand zerouri. Pozitiile ın care se afla coeficientii nenuli rezulta dinmatricea de conexiuni a retelei.

In sfarsit, avand determinate aproximantele deformatiilor, leputem introduce ın legea lui Hooke (4.4). Obtinem astfel aproximarielementale ale tensiunilor σx, σy, σz, τxy, τyz si τzx:

σ(e)x =

E(e)

[1 + ν(e)] [1 − 2ν(e)]

1 − ν(e)

ε(e)

x + ν(e)ε(e)

y + ε(e)z

©,

σ(e)y =

E(e)

[1 + ν(e)] [1 − 2ν(e)]

1 − ν(e)

ε(e)

y + ν(e)ε(e)

z + ε(e)x

©,

σ(e)z =

E(e)

[1 + ν(e)] [1 − 2ν(e)]

1 − ν(e)

ε(e)

z + ν(e)ε(e)

x + ε(e)y

©,

Page 119: Metoda Elementelor Finite

Modelul cu elemente finite al problemelor de elasticitate 113

τ (e)xy =

E(e)

2 [1 + ν(e)]γ(e)

xy ,

τ (e)yz =

E(e)

2 [1 + ν(e)]γ(e)

yz ,

τ (e)zx =

E(e)

2 [1 + ν(e)]γ(e)

zx .

(4.22)

In relatiile de mai sus, E(e) si ν(e) sunt modulul lui Young, respec-tiv coeficientul lui Poisson corespunzatoare elementului (e). Acestemarimi nu trebuie sa fie neaparat identice pentru toate elementeleretelei. MEF ne permite asadar sa modelam cu usurinta compor-tarea elastica a unor corpuri neomogene.

Daca tinem cont de (4.18), putem rescrie ecuatiile (4.22) suburmatoarea forma matriceala:

σ(e)©

=C(e)

ε(e)

©, e = 1, 2, . . . , nelem, (4.23)

ın care σ(e)

©=σ(e)

x , σ(e)y , σ(e)

z , τ (e)xy , τ (e)

yz , τ (e)zx

T(4.24)

este un vector-coloana ce grupeaza aproximantele elementale aletensiunilor, iar

C(e)

=

E(e)

1 + ν(e)

266666666666666641−ν(e)

1−2ν(e)ν(e)

1−2ν(e)ν(e)

1−2ν(e) 0 0 0

ν(e)

1−2ν(e)1−ν(e)

1−2ν(e)ν(e)

1−2ν(e) 0 0 0

ν(e)

1−2ν(e)ν(e)

1−2ν(e)1−ν(e)

1−2ν(e) 0 0 0

0 0 0 12

0 0

0 0 0 0 12

0

0 0 0 0 0 12

37777777777777775 (4.25)

este asa-numita matrice constitutiva elastica a elementului (e).Dispunem acum de aproximantele elementale ale tuturor ne-

cunoscutelor problemei. Putem merge mai departe, ınlocuindu-leın expresia (4.12) a energiei potentiale elastice. Vom obtine astfelo aproximare a lui Π care depinde numai de deplasarile nodale:

Page 120: Metoda Elementelor Finite

114 Capitolul 4

Π → Π (u1, v1, w1, u2, v2, w2, . . . , unnod, vnnod

, wnnod) . (4.26)

Pentru a vedea care este expresia lui Π, vom lua separat fiecaretermen din membrul drept al lui (4.12) si ıi vom construi aproxi-marea de tip element finit.

Incepem cu energia de deformatie Πdef data de relatia (4.7).Aplicand proprietatea de aditivitate a integralei ın raport cu dome-niul, obtinem urmatoarea aproximare a lui Πdef :

Πdef =nelemXe=1

Π(e)def , (4.27)

ın care

Π(e)def =

ZΩ(e)

1

2

σ(e)

x ε(e)x + σ(e)

y ε(e)y + σ(e)

z ε(e)z +

+τ (e)xy γ(e)

xy + τ (e)yz γ(e)

yz + τ (e)zx γ(e)

zx

(4.28)

este contributia unui element oarecare (e). Am notat cu Ω(e) dome-niul spatial ocupat de elementul (e).

Cu ajutorul lui (4.18) si (4.24), putem rescrie (4.28) sub formamatriceala:

Π(e)def =

ZΩ(e)

1

2

ε(e)

©T σ(e)

©dΩ. (4.29)

Inlocuind ın (4.29) aproximantele (4.23) si (4.20), obtinem urma-toarea expresie a energiei de deformatie pentru elementul (e):

Π(e)def = dT

ZΩ(e)

1

2

B(e)

T C(e)

B(e)

dΩd . (4.30)

Fiindca deplasarile nodale sunt marimi asociate unei multimi dis-crete de puncte (nodurile retelei), vectorul-coloana d poate fi pla-sat ın afara expresiei de integrare.

Revenind cu (4.30) ın (4.27), gasim aproximarea de tip elementfinit a energiei de deformatie Πdef :

Πdef = dT

nelemXe=1

ZΩ(e)

1

2

B(e)

T C(e)

B(e)

!d . (4.31)

Page 121: Metoda Elementelor Finite

Modelul cu elemente finite al problemelor de elasticitate 115

Lucrul mecanic al fortelor volumice definit de relatia (4.9)fiind tot o marime de tip integral, poate fi aproximat ın manieraasemanatoare:

Wb =nelemXe=1

W(e)b , (4.32)

unde

W(e)b =

ZΩ(e)

hb(e)x u(e) + b(e)

y v(e) + b(e)z w(e)

idΩ (4.33)

este contributia unui element oarecare (e). In (4.33) am notat cub(e)x , b(e)

y si b(e)z aproximantele componentelor fortei volumice aso-

ciate elementului (e). Daca aceasta forta este constanta ın domeniulΩ(e), marimile b(e)

x , b(e)y si b(e)

z sunt chiar componentele ei. Se poateıntampla ınsa ca forta volumica sa se modifice de la un punct laaltul. In astfel de situatii, b(e)

x , b(e)y si b(e)

z se construiesc ca aproxi-mante de tip element finit, folosind valorile fortei volumice ın noduri(vezi [9, 10] pentru detalii suplimentare).

Cu ajutorul lui (4.14), putem rescrie (4.33) sub forma matri-ceala:

W(e)b =

ZΩ(e)

nd(e)

oT nb(e)

odΩ, (4.34)

unde nb(e)

o=hb(e)x , b(e)

y , b(e)z

iT(4.35)

este un vector-coloana ce grupeaza aproximantele elementale alecomponentelor fortei volumice. Inlocuind (4.13) ın (4.34), obtinemurmatoarea expresie a lucrului mecanic al fortelor volumice asociateelementului (e):

W(e)b = dT

ZΩ(e)

N (e)

T nb(e)

odΩ. (4.36)

Revenind cu (4.36) ın (4.32), gasim aproximarea de tip elementfinit a lui Wb:

Wb = dT

nelemXe=1

ZΩ(e)

N (e)

T nb(e)

odΩ

!. (4.37)

Page 122: Metoda Elementelor Finite

116 Capitolul 4

Lucrul mecanic al fortelor de suprafata definit de (4.10) poatefi aproximat astfel:

Ws =nelemXe=1

W (e)s , (4.38)

unde

W (e)s =

(e)s

s(e)

x u(e) + s(e)y v(e) + s(e)

z w(e)dΣ (4.39)

este contributia unui element oarecare (e). Am notat cu Σ(e)s

portiunea din Σs care revine elementului (e). Evident, daca acestanu are nici o regiune comuna cu Σs, suprafata Σ(e)

s este o multimevida, iar W (e)

s se anuleaza. Marimile s(e)x , s(e)

y si s(e)z sunt aproximante

ale componentelor fortei de suprafata care actioneaza asupra ele-mentului (e). Observatiile pe care le-am facut ın legatura cu b(e)

x ,b(e)y si b(e)

z raman valabile si pentru s(e)x , s(e)

y si s(e)z .

Cu ajutorul lui (4.14), putem rescrie (4.39) sub forma matri-ceala:

W (e)s =

(e)s

nd(e)

oT s(e)

©dΣ, (4.40)

unde s(e)

©=s(e)

x , s(e)y , s(e)

z

T(4.41)

este un vector-coloana ce grupeaza aproximantele elementale alecomponentelor fortei de suprafata. Inlocuind (4.13) ın (4.40),obtinem urmatoarea expresie a lucrului mecanic al fortelor desuprafata asociate elementului (e):

W (e)s = dT

(e)s

N (e)

T s(e)

©dΣ. (4.42)

Revenind cu (4.42) ın (4.38), gasim aproximarea de tip elementfinit a lui Ws:

Ws = dT

nelemXe=1

(e)s

N (e)

T s(e)

©dΣ

!. (4.43)

Page 123: Metoda Elementelor Finite

Modelul cu elemente finite al problemelor de elasticitate 117

Ultimul termen din membrul drept al lui (4.12) este lucrul me-canic al fortelor punctuale definit de (4.11). Observam ca el nu maieste o marime de tip integral, ci o suma discreta. Pentru a-l in-troduce ın modelul cu elemente finite, ne vom folosi de faptul caın punctul de aplicatie al fiecarei forte concentrate se afla un nod.Vom construi atunci un vector-coloana P, ın care vom asociacate trei pozitii succesive fiecarui nod al retelei. Daca ın respecti-vul nod actioneaza o forta punctuala, cele trei pozitii vor gazduichiar componentele fortei. Daca nodul nu are asociata nici o ase-menea forta, pozitiile corespunzatoare vor contine zerouri. Asadar,vectorul-coloana P are structura

P = [Px,1, Py,1, Pz,1, Px,2, Py,2, Pz,2, . . . ,

Px,nnod, Py,nnod

, Pz,nnod]T .

(4.44)

Cu ajutorul lui si al vectorului-coloana d definit de (4.15), lucrulmecanic al fortelor punctuale poate fi exprimat sub forma matri-ceala:

Wp = dT P . (4.45)

Am determinat astfel aproximarile tuturor termenilor care in-tervin ın expresia energiei potentiale elastice. Propriu-zis, prinınlocuirea lor ın (4.12), putem scrie

Π = Πdef − Wb − Ws − Wp (4.46)

sau, ıntr-o formulare explicita,

Π = dT

nelemXe=1

ZΩ(e)

1

2

B(e)

T C(e)

B(e)

!d−

− dT

nelemXe=1

ZΩ(e)

N (e)

T nb(e)

odΩ

!−

− dT

nelemXe=1

(e)s

N (e)

T s(e)

©dΣ

!− dT P .

(4.47)

Page 124: Metoda Elementelor Finite

118 Capitolul 4

(4.47) ne arata clar faptul ca marimea Π nu mai contine ca para-metri nedeterminati decat deplasarile nodale grupate ın vectorul-coloana d.

Noi stim ınsa ca solutia sistemului de ecuatii cu derivate partialeprezentat ın §4.1 realizeaza minimul lui Π. Vom cauta atunci oaproximare a acestei solutii impunand conditia ca deplasarile nodalesa minimizeze marimea Π definita de (4.47).

In acest scop, vom folosi teorema lui Lagrange, care ne daurmatoarea conditie necesara de extrem pentru Π (interpretata cao functie de mai multe variabile):

∂ Π

∂ u1

= 0,∂ Π

∂ v1

= 0,∂ Π

∂ w1

= 0,

∂ Π

∂ u2

= 0,∂ Π

∂ v2

= 0,∂ Π

∂ w2

= 0,

. . .

∂ Π

∂ unnod

= 0,∂ Π

∂ vnnod

= 0,∂ Π

∂ wnnod

= 0.

(4.48)

De fapt, (4.48) este un sistem de ecuatii avand ca necunoscutedeplasarile nodale. Prin rezolvarea lui numerica putem determinacomponentele vectorului-coloana d care minimizeaza marimea Π.Realizarea acestui deziderat face necesara deducerea unei forme ex-plicite a sistemului (4.48).

Pentru a pastra concizia relatiilor matematice, vom folosi unoperator matriceal de derivare definit ın felul urmator:

∂ d =

∂ u1

,∂

∂ v1

,∂

∂ w1

,∂

∂ u2

,∂

∂ v2

,∂

∂ w2

, . . . ,

∂ unnod

,∂

∂ vnnod

,∂

∂ wnnod

T

.

(4.49)

Dupa cum se observa, coeficientii sai sunt derivate partiale ın raportcu deplasarile nodale. Ei actioneaza asupra marimii scalare careia

Page 125: Metoda Elementelor Finite

Modelul cu elemente finite al problemelor de elasticitate 119

i se aplica respectivul operator. De exemplu, (4.49) ne permite sarescriem sistemul de ecuatii (4.48) sub urmatoarea forma concisa:

∂ Π

∂ d = 0 , (4.50)

ın care 0 este un vector-coloana nul cu un numar de componenteegal cu triplul numarului total de noduri din retea.

Operatorul matriceal de derivare ∂/∂ d are proprietati asema-natoare cu ale derivatei partiale clasice. Din (4.46) rezulta atunci

∂ Π

∂ d =∂ Πdef

∂ d − ∂ Wb

∂ d − ∂ Ws

∂ d − ∂ Wp

∂ d· (4.51)

Va trebui deci sa explicitam fiecare termen din membrul drept allui (4.51) aplicand operatorul de derivare ∂/∂ d asupra relatiilor(4.31), (4.37), (4.43) si (4.45).

Folosind formulele de derivare matriceala deduse ın anexa A,putem scrie imediat urmatoarele ecuatii pentru evaluarea marimilormentionate mai sus:

∂ Πdef

∂ d =

nelemXe=1

ZΩ(e)

B(e)

T C(e)

B(e)

!d , (4.52)

∂ Wb

∂ d =nelemXe=1

ZΩ(e)

N (e)

T nb(e)

odΩ, (4.53)

∂ Ws

∂ d =nelemXe=1

(e)s

N (e)

T s(e)

©dΣ, (4.54)

∂ Wp

∂ d = P . (4.55)

(4.51) - (4.55) permit rescrierea sistemului (4.50) sub forma

[K] d = F , (4.56)

unde

Page 126: Metoda Elementelor Finite

120 Capitolul 4

[K] =nelemXe=1

K(e)

, (4.57)

F = P +nelemXe=1

F (e)

©, (4.58)

K(e)=ZΩ(e)

B(e)

T C(e)

B(e)

dΩ, (4.59)

F (e)©

=ZΩ(e)

N (e)

T nb(e)

odΩ +

(e)s

N (e)

T s(e)

©dΣ. (4.60)

(4.56) se numeste sistem de ecuatii global (sau sistem de ecuatii

nodale). Numele sau se justifica prin faptul ca descrie raspunsulıntregii retele de elemente finite. Matricile si vectorii-coloanadin structura lui (4.56) au de asemenea denumiri consacrate:[K] – matricea globala de rigiditate;F – vectorul global al fortelor ;d – vectorul gradelor de libertate nodale;K(e)

– matricea de rigiditate a elementului (e);

F (e)©

– vectorul fortelor asociate elementului (e).Teminologia de mai sus ısi are explicatia ın analogia dintre sis-

temul (4.56) si formula care descrie comportarea mecanica a unuiresort liniar. Din perspectiva acestei analogii, matricea [K] joacarolul constantei elastice a resortului (cunoscuta si sub numele derigiditate), vectorul-coloana F este echivalentul fortei exterioarecare solicita resortul, iar d reprezinta elongatia.

Este important sa subliniem unele proprietati deosebit de avan-tajoase ale sistemului de ecuatii global.

In primul rand, sa remarcam maniera repetitiva ın care potfi construite matricea globala de rigiditate si vectorul global alfortelor. Intr-adevar, din (4.57) se constata ca matricea [K] re-zulta prin ınsumarea unor contributii elementale. La fel, (4.58) nearata ca vectorul-coloana F poate fi asamblat la randul sau dintermeni asociati elementelor individuale. Face exceptie vectorul-coloana P, care trebuie construit printr-o procedura de asam-blare la nivel de noduri si apoi adunat la F (totusi, aceasta estela fel de repetitiva ca si asamblarea la nivel de elemente).

Page 127: Metoda Elementelor Finite

Modelul cu elemente finite al problemelor de elasticitate 121

O alta caracteristica pozitiva a sistemului de ecuatii (4.56) estelegata de procedura folosita pentru evaluarea contributiilor elemen-tale

K(e)

siF (e)

©. Relatiile (4.59) si (4.60) ne arata ca acestea

sunt integrale definite pe domeniile ocupate de elementele indivi-duale sau pe frontierele lor. In capitolul 3 am vazut ca orice elementpoate fi construit pornind de la un

”parinte” de forma foarte simpla.

Devine atunci posibil transferul operatiei de integrare din spatiulreal ın spatiul coordonatelor normalizate ale

”parintelui” (ın care

forma domeniului de integrare este regulata si, lucrul cel mai impor-tant, ıntotdeauna aceeasi). Programele MEF folosesc fara exceptieaceasta tehnica, fiindca ea le permite sa trateze ın maniera unitaraelementele retelei, indiferent de forma lor reala. Pentru calculul inte-gralelor din relatiile (4.59) si (4.60) sunt utilizate procedee numerice(formule de cuadratura de tip Gauss, Newton-Cotes, Simpson etc.[3, 4, 5, 8, 9, 10]).

In sfarsit, sa mai observam faptul ca matricea globala de rigi-ditate este simetrica (vezi relatiile (4.57) si (4.59)). Am evidentiataceasta proprietate a sistemului de ecuatii global si ın exemplul tra-tat ın capitolul 1. Acum ne dam seama ca ea este o caracteristicagenerala a modelelor cu elemente finite pentru probleme de elasti-citate. Am mentionat cu prilejul rezolvarii numerice a problemei-exemplu din capitolul 1 ca matricea globala de rigiditate mai este side tip banda (ın sensul ca are toti coeficientii nenuli grupati ıntr-ofasie situata de o parte si de alta a diagonalei). Programele MEFvor putea atunci sa stocheze ın memorie doar partea utila a lui[K] (ın general, este pastrata semibanda superioara). Prin aceasta,ele economisesc resursele calculatorului si castiga ın viteza, fiindcaevita efectuarea operatiilor ın care intervin coeficienti nuli ai ma-tricii globale de rigiditate.

O problema foarte importanta pe care trebuie sa o discutameste aceea a introducerii conditiilor la limita de tip deplasare ınmodelul cu elemente finite. Intr-adevar, daca analizam etapele par-curse pana acum ın vederea construirii sistemului de ecuatii glo-bal, constatam ca restrictiile de miscare (4.5) nu se reflecta nicaieri

Page 128: Metoda Elementelor Finite

122 Capitolul 4

ın demersul nostru. Din aceasta cauza, modelul cu elemente finiteobtinut sufera de o nedeterminare matematica. Propriu-zis, matri-cea [K] este singulara, deci sistemul (4.56) are o infinitate de solutii.Aceasta situatie evidentiaza o realitate fizica: prin neimpunerea niciunei restrictii de miscare, reteaua de elemente finite poate fi depla-sata ca un rigid ın orice pozitie a spatiului (cu alte cuvinte, putemobtine o infinitate de solutii valide ale lui (4.56) prin simpla adu-nare a componentelor unui vector de translatie arbitrar la o solutieparticulara).

In metoda elementelor finite, procedura de impunere a con-ditiilor la limita de tip deplasare opereaza asupra nodurilor carese afla pe portiunea Σs a frontierei. Cade ın sarcina analistu-lui specificarea restrictiilor de miscare pentru respectivele noduri.Totusi, trebuie sa remarcam faptul ca aproximantele (4.13) ale com-ponentelor campului vectorilor deplasare extind efectul restrictiilornodale asupra unor ıntregi portiuni de pe frontiera. Constatam ast-fel ca, desi opereaza la nivelul nodurilor, strategia adoptata de MEFreflecta caracterul de suprafata al conditiilor (4.5).

Exemplele din capitolele 1 si 2 ne-au permis sa prezentam ma-niera practica ın care se introduc restrictiile de miscare ıntr-unmodel cu elemente finite. In esenta, procedura se reduce la a su-prima din sistemul (4.56) ecuatiile asociate deplasarilor nodalerestrictionate. Rezulta astfel un sistem redus care nu mai suferade nici o nedeterminare matematica. Pentru rezolvarea lui nume-rica, programele MEF folosesc de regula variante ale algoritmuluide eliminare Gauss [3, 4, 5, 8, 9, 10]).

Solutia lui (4.56) este reprezentata de vectorul-coloana d careminimizeaza cantitatea Π. Cu ajutorul lui si al formulelor de apro-ximare (4.13), (4.20), respectiv (4.23), se pot obtine ın continuarevalorile necunoscutelor problemei ın orice punct al domeniului deanaliza. Intr-o prima etapa, se ınlocuieste d ın (4.13). Rezulta

astfel aproximantele elementalend(e)

oale componentelor campului

vectorilor deplasare. De asemenea, prin ınlocuirea lui d ın (4.20)se obtin aproximantele elementale

ε(e)

©ale deformatiilor. In etapa

Page 129: Metoda Elementelor Finite

Probleme bidimensionale de elasticitate 123

finala, se introduc vectorii-coloanaε(e)

©ın (4.23). Rezulta astfel

aproximantele elementaleσ(e)

©ale tensiunilor. In general, progra-

mele MEF comerciale stocheaza rezultatele numerice ale analizei ınfisiere pe disc, dar le reprezinta si sub forma unor diagrame (carepot fi vizualizate pe ecran sau tiparite la imprimanta). Folosireaunei interfete grafice simplifica mult interpretarea datelor de iesire.

Discutia de mai sus a urmarit evidentierea principalelor carac-teristici ale procesului de rezolvare numerica a problemelor de elas-ticitate prin metoda elementelor finite. Din dorinta de a pastrageneralitatea expunerii, am considerat cazul deformarii corpurilortridimensionale. In practica, inginerii mecanici ıntalnesc totusi des-tul de des probleme de elasticitate mai particulare, a caror tratarespatiala se dovedeste greoaie sau ineficienta. In paragraful urmatorvom descrie trei tipuri de asemenea probleme care se pot reduce lamodele cu elemente finite bidimensionale.

4.3 Probleme bidimensionale de

elasticitate

Se ıntampla frecvent sa ıntalnim probleme de elasticitate alecaror necunoscute depind de numai doua coordonate spatiale. Deexemplu, sunt situatii ın care vectorii deplasare, deformatiile si ten-siunile din corp depind numai de coordonatele x si y, dar nu depindde z. Domeniul de analiza Ω asociat unor astfel de probleme va aveanumai doua dimensiuni. Drept consecinta, si modelul cu elementefinite va trebui sa utilizeze elemente triunghiulare sau patrulatere(de felul celor prezentate ın capitolul 3).

Este foarte important ca analistii sa identifice corect problemelepractice pe care le pot asimila unor cazuri bidimensionale. Din acestpunct de vedere, lucrurile nu sunt ıntotdeauna foarte clare. Ambi-guitatea ısi are originea ın principal ın faptul ca multe probleme cucare se confrunta inginerii nu pot fi considerate ca riguros bidimen-sionale. In astfel de situatii, analistii cu oarecare experienta cauta

Page 130: Metoda Elementelor Finite

124 Capitolul 4

sa evalueze erorile pe care le-ar introduce utilizarea unui model bi-dimensional. Ei vor adopta un astfel de model numai daca abaterileestimate se ıncadreaza ıntr-un interval acceptabil. Pentru a fi ca-pabili sa aprecieze ın mod realist aceste erori, analistii trebuie saposede suficiente informatii despre problemele bidimensionale stan-dard din teoria elasticitatii.

Iata de ce am considerat oportuna includerea unui paragraf con-sacrat prezentarii caracteristicilor esentiale ale urmatoarelor tipuride probleme bidimensionale ce apar mai frecvent ın aplicatii:

• probleme de stare plana de tensiuni;

• probleme de stare plana de deformatii;

• probleme cu simetrie axiala.

Spatiul limitat nu ne permite sa descriem ın toate detaliile modelelecu elemente finite asociate. Vom discuta numai particularitatile careindividualizeaza respectivele probleme bidimensionale si, mai ales,vom indica unele cazuri practice pe care le putem asimila lor.

4.3.1 Probleme de stare plana de tensiuni

Consideram un corp subtire de tip placa plana (fig. 4.2) supusunor ıncarcari orientate paralel cu fetele sale. Presupunem ca res-pectivele ıncarcari sunt repartizate uniform pe grosime.

Vom numi plan median locul geometric al punctelor situate laegala distanta de fetele placii. Pentru descrierea matematica a de-formarii elastice a unui corp de felul celui din figura 4.2 folosim unsistem de coordonate carteziene ortogonale ale carui axe x si y seafla ın planul median.

Particularitatile geometrice si mecanice mentionate mai sus, neconduc la concluzia ca starea de tensiuni si deformatii este aceeasiın toate punctele amplasate de-a lungul unei normale la planul me-dian. Cu alte cuvinte, punctele aflate ın planul median sunt repre-zentative din acest punct de vedere pentru toate particulele aflate

Page 131: Metoda Elementelor Finite

Probleme bidimensionale de elasticitate 125

x

y

z

Plan median

~s

==

h

Figura 4.2: Placa plana subtire aflata ın stare plana de tensiuni

pe normalele corespondente. Deformarea elastica a placii poate fideci modelata ın maniera bidimensionala, folosind numai sectiuneamediana.

Sub aspect matematic, observatiile de mai sus sunt echivalentecu faptul ca deformatiile si tensiunile nu depind de coordonata z(ele sunt functii numai de x si y). In plus, caracteristicile sistemu-lui de forte exterioare au drept consecinta anularea unora dintredeformatii si tensiuni ın toate punctele placii:

γyz = 0, γzx = 0, (4.61)

σz = 0, τyz = 0, τzx = 0. (4.62)

Relatiile (4.61) - (4.62) definesc asa-numita stare plana de tensiuni.Modelul cu elemente finite al problemelor de stare plana de

tensiuni este mai simplu decat cel prezentat ın §4.2. Aceasta dato-rita faptului ca numarul necunoscutelor este mult mai mic (o serieıntreaga de deformatii si tensiuni fiind nule, calculul lor nu-si maiare sensul).

De exemplu, vectorul-coloanaσ(e)

©are acum urmatoarea

structura:

Page 132: Metoda Elementelor Finite

126 Capitolul 4σ(e)

©=σ(e)

x , σ(e)y , τ (e)

xy

T. (4.63)

Observam ca ın el nu mai sunt stocate decat tensiunile nenule (vezirelatiile (4.62)). Vectorul-coloana

ε(e)

©se prezinta sub o forma

corespunzator redusa:ε(e)

©=ε(e)

x , ε(e)y , γ(e)

xy

T. (4.64)

Sa remarcam faptul ca ın structura lui nu intra deformatia ε(e)z ,

desi ea este ın general diferita de zero (vezi relatiile (4.61)). Acestlucru se explica prin faptul ca ıngrosarea sau subtierea placii nuare efect energetic, ca urmare a anularii tensiunii σ(e)

z . In concluzie,introducerea lui ε(e)

z ın vectorul-coloanaε(e)

©ar fi un balast pen-

tru modelul cu elemente finite. Oricum, dupa evaluarea tensiunilornenule, deformatia ε(e)

z poate fi determinata cu urmatoarea relatie(obtinuta prin inversarea legii lui Hooke generale):

ε(e)z = − ν(e)

E(e)

σ(e)

x + σ(e)y

. (4.65)

Matricea constitutiva elastica asociata problemelor de stareplana de tensiuni are structura

C(e)=

E(e)

1 + ν(e)

26664 11−ν(e)

ν(e)

1−ν(e) 0

ν(e)

1−ν(e)1

1−ν(e) 0

0 0 12

37775 . (4.66)

Sistemul de ecuatii global care se obtine prin aplicarea MEFare tot o forma de tipul (4.56). Trebuie sa retinem ınsa faptul caın vectorul-coloana d nu mai sunt stocate decat deplasarile u siv ale nodurilor:

d = [u1, v1, u2, v2, . . . , unnod, vnnod

]T . (4.67)

Pentru a ıntelege de ce nu este necesara utilizarea deplasarilor no-dale pe directia axei z, este suficient sa ne reamintim ca discreti-zarea se realizeaza exclusiv ın planul median al placii (folosind ele-mente finite triunghiulare sau patrulatere), iar particulele situate

Page 133: Metoda Elementelor Finite

Probleme bidimensionale de elasticitate 127

ın aceasta sectiune au deplasari w nule. In plus, pentru calcululdeformatiilor ε(e)

x , ε(e)y si γ(e)

xy avem nevoie numai de aproximanteledeplasarilor u si v (vezi relatiile (4.16)).

Integralele de volum sau de suprafata care intra ın structuramatricii globale de rigiditate si a vectorului global al fortelor (vezi(4.57) - (4.60)) se reduc acum la integrale de arie, respectiv decontur ınmultite cu grosimea h a placii (fig. 4.2).

Problema starii plane de tensiuni, asa cum este ea formulata ınteoria elasticitatii, se dovedeste restrictiva. Se ıntampla destul derar ca o aplicatie practica sa aiba toate caracteristicile geometricesi mecanice mentionate la ınceputul acestui paragraf.

De exemplu, sunt situatii ın care grosimea placii nu este perfectuniforma. Din aceasta cauza, chiar daca toate celelalte conditii aleproblemei sunt respectate, relatiile (4.61) - (4.62) nu mai sunt ri-guros valabile. Cade ın sarcina analistului sa evalueze importantavariatiilor de grosime si gradul ın care ele perturba starea plana detensiuni. In functie de nivelul estimat al erorilor si de obiectivelesale, utilizatorul poate adopta un model bidimensional aproximativsau un model tridimensional mai riguros. De fapt, majoritatea pro-gramelor MEF comerciale admit construirea unor modele de stareplana de tensiuni ın care grosimea elementelor nu este uniforma.

O alta situatie delicata este aceea ın care fortele care solicitaplaca nu sunt riguros paralele cu planul ei median. Din nou, ana-listul trebuie sa evalueze cu atentie importanta tensiunilor normaleinduse de componentele neplanare ale ıncarcarilor (ın principal prinraportare la tensiunile σx, σy si τxy). Un model bidimensional se do-vedeste acceptabil numai daca perturbatiile constatate sunt mici.

In sfarsit, multi utilizatori de programe MEF ısi pun ıntrebarea:cand putem considera ca o placa este cu adevarat subtire? Altfelspus, exista vreo limita dincolo de care o placa trebuie modelataca un corp tridimensional? Majoritatea analistilor experimentatiadmit ca ipotezele starii plane de tensiuni sunt bine aproximate atattimp cat raportul dintre grosimea placii si cea mai mica dimensiuneplanara a acesteia este inferior lui 1/10. La valori mai mari ale

Page 134: Metoda Elementelor Finite

128 Capitolul 4

respectivului raport, starea de tensiuni si deformatii din corpul realse ındeparteaza din ce ın ce mai mult de conditiile (4.61) - (4.62),astfel ıncat devine necesara adoptarea unui model mai precis (deexemplu, pot fi utilizate elemente finite de tip placa groasa [11] sauchiar elemente finite tridimensionale).

4.3.2 Probleme de stare plana de deformatii

Un corp se afla ın stare plana de deformatii daca toate parti-culele sale au deplasari paralele cu un acelasi plan numit planul

problemei si, ın plus, deplasarile particulelor situate de-a lungulunei normale oarecare la acest plan sunt egale. Pentru a descriedeformarea elastica a unui asemenea corp, vom folosi un sistemde coordonate carteziene ortogonale ale carui axe x si y se afla ınplanul problemei.

Starea plana de deformatii este imposibil de obtinut ın prac-tica. Se poate ıntampla cel mult sa ıntalnim situatii care se apropiesuficient de bine de conditiile enumerate mai sus. Un astfel de exem-plu, este cel al unui corp cilindric foarte lung, cu sectiune transver-sala oarecare, supus unor ıncarcari uniform distribuite de-a lungulgeneratoarei si orientate perpendicular fata de aceasta (fig. 4.3). Lao distanta suficient de mare de capetele cilindrului, se poate con-sidera ca ipotezele starii plane de deformatii sunt bine aproximate.In general, sunt asimilate unor stari plane de deformatii proble-mele ın care sunt implicate corpuri de lungime mare si sectiunetransversala practic constanta, supuse actiunii unor ıncarcari uni-form distribuite ın lungul generatoarei si orientate perpendicular peaceasta (de exemplu, barajele solicitate de presiunea apei din laculde acumulare, rolele cilindrice ale laminoarelor etc.).

Datorita faptului ca deplasarile particulelor situate pe aceeasinormala la planul coordonatelor x, y sunt identice si perpendicularepe axa z, lungimea corpului nu influenteaza solutia problemei. Defapt, se constata ca sectiunile transversale ale corpului sunt perfectechivalente din punct de vedere al starii de tensiuni si deformatii.

Page 135: Metoda Elementelor Finite

Probleme bidimensionale de elasticitate 129

Planul problemei

x

y

z

~s

Figura 4.3: Cilindru lung aflat ın stare plana de deformatii

Page 136: Metoda Elementelor Finite

130 Capitolul 4

Din aceasta cauza, programele MEF ignora pur si simplu lungi-mea corpurilor aflate ın stare plana de deformatii. Utilizatorul tre-buie asadar sa discretizeze doar o sectiune transversala, folosindelemente finite bidimensionale.

Ca urmare a particularitatilor mentionate mai sus, componen-tele u si v ale campului vectorilor deplasare nu depind de coordo-nata z, iar componenta w este identic nula pentru toate particulelecorpului:

u = u(x, y), v = v(x, y), w = 0. (4.68)

Din (4.68) si (4.2) constatam imediat ca o parte dintre deformatiisunt nule:

εz = 0, γyz = 0, γzx = 0. (4.69)

Mai departe, din (4.69) si legea lui Hooke (4.4) rezulta ca si douatensiuni tangentiale sunt egale cu zero:

τyz = 0, τzx = 0. (4.70)

Conditiile (4.68) - (4.70) caracterizeaza sub aspect matematicstarea plana de deformatii.

Datorita anularii unora dintre necunoscute, modelul cu elementefinite asociat problemelor de stare plana de deformatii este mai sim-plu decat cel prezentat ın §4.2. Vectorii-coloana

σ(e)

©siε(e)

©sunt

formal identici cu cei corespunzatori starii plane de tensiuni (vezi(4.63) - (4.64)). Este important sa retinem ınsa ca asemanarea estepur matematica. De fapt, se observa imediat ca vectorul-coloanaσ(e)

©nu include tensiunea σ(e)

z , desi ea este ın general diferita dezero (vezi (4.70)). Acest lucru se explica prin faptul ca ea nu areefect energetic, ca urmare a anularii deformatiei corespunzatoareε(e)

z (vezi (4.69)). Introducerea lui σ(e)z ın vectorul-coloana

σ(e)

©ar

fi astfel un balast pentru modelul cu elemente finite. Oricum, dupaevaluarea deformatiilor nenule, tensiunea σ(e)

z poate fi determinatacu urmatoarea relatie (obtinuta direct din legea lui Hooke (4.4)):

σ(e)z =

E(e)ν(e)

[1 + ν(e)] [1 − 2ν(e)]

ε(e)

x + ε(e)y

. (4.71)

Page 137: Metoda Elementelor Finite

Probleme bidimensionale de elasticitate 131

Matricea constitutiva elastica asociata problemelor de stareplana de deformatii este o forma redusa a lui (4.25):

C(e)=

E(e)

1 + ν(e)

26664 1−ν(e)

1−2ν(e)ν(e)

1−2ν(e) 0

ν(e)

1−2ν(e)1−ν(e)

1−2ν(e) 0

0 0 12

37775 . (4.72)

Sistemul de ecuatii global care se obtine prin aplicarea MEFeste tot de tipul (4.56). La fel ca ın cazul starii plane de tensiuni,vectorul-coloana d contine numai deplasarile u si v ale noduri-lor (vezi (4.67)). Deplasarea w fiind oricum nula, includerea ei ınmodelul cu elemente finite ar fi redundanta.

Integralele de volum sau de suprafata care intra ın structuramatricii globale de rigiditate si a vectorului global al fortelor (vezi(4.57) - (4.60)) se reduc acum la integrale de arie, respectiv decontur ınmultite cu o grosime fictiva egala cu unu.

4.3.3 Probleme cu simetrie axiala

Problemele din aceasta categorie sunt caracterizate de faptul caatat geometria corpului, cat si rezemarile sale (conditiile la limita detip deplasare), respectiv ıncarcarile exterioare au simetrie cilindrica(fig. 4.4). In consecinta, starea de tensiuni si deformatii rezultataare si ea simetrie axiala. O astfel de simetrie exista numai dacafiecare particula a corpului ramane ın acelasi plan axial pe toatadurata procesului de deformare elastica.

In aceste conditii, este mai convenabila ınlocuirea sistemului decoordonate carteziene cu un sistem cilindric (r, θ, z) a carui axa zcoincide cu axa de simetrie a corpului (fig. 4.4). Fiindca sectiunileaxiale sunt perfect echivalente din punct de vedere al starii de ten-siuni si deformatii, toate necunoscutele problemei trebuie sa fie in-dependente de coordonata unghiulara θ (ele vor fi deci functii numaide r si z).

Page 138: Metoda Elementelor Finite

132 Capitolul 4

x

y

directieaxiala

directiecircumferentiala

Particulele aflate initial peaxa de simetrie a corpuluitrebuie sa aiba deplasariradiale nule

r

zz

directie radiala

θ

r

~s

Figura 4.4: Problema cu simetrie axiala si tratarea ei bidimensionalaıntr-un sistem de coordonate cilindric

Page 139: Metoda Elementelor Finite

Probleme bidimensionale de elasticitate 133

Problemele din aceasta categorie vor putea fi atunci tratate ınmaniera bidimensionala prin considerarea unei sectiuni axiale oare-care pozitionate ıntr-un plan de coordonate (r, z). Utilizatorul MEFva discretiza o astfel de sectiune cu ajutorul unor elemente tri-unghiulare sau patrulatere de tipul celor prezentate ın capitolul 3.

Notam cu ur, vθ si wz componentele radiala, circumferentiala,respectiv axiala ale campului vectorilor deplasare (le-am prevazutcu indici pentru a preciza explicit faptul ca sunt exprimate ıntr-unsistem de coordonate cilindric si nu ıntr-unul cartezian). Fiindcaparticulele corpului raman ın acelasi plan axial, este necesar ca vθ

sa fie o cantitate identic nula. In ceea ce priveste componentele ur siwz, ele trebuie sa fie independente de coordonata θ. Aceste conditiise exprima ın limbaj matematic sub forma

ur = ur(r, z), vθ = 0, wz = wz(r, z). (4.73)

In cazul problemelor cu simetrie axiala vor fi nenule numai urma-toarele deformatii:

εr – deformatia liniara ın directie radiala;

εθ – deformatia liniara ın directie circumferentiala;

εz – deformatia liniara ın directie axiala;

γzr – deformatia unghiulara ın plan axial.

Aceste marimi sunt definite prin relatiile [5, 9, 10]

εr =∂ ur

∂ r, εθ =

ur

r, εz =

∂ wz

∂ z, γzr =

∂ wz

∂ r+

∂ ur

∂ z· (4.74)

In mod corespunzator, vor fi nenule numai urmatoarele tensiuni:

σr – tensiunea normala ın directie radiala;

σθ – tensiunea normala ın directie circumferentiala;

Page 140: Metoda Elementelor Finite

134 Capitolul 4

σz – tensiunea normala ın directie axiala;

τzr – tensiunea tangentiala ın plan axial.

Ele sunt legate de deformatii prin legea lui Hooke [5, 9, 10]:

σr =E(1 − ν)

(1 + ν)(1 − 2ν)εr +

(1 + ν)(1 − 2ν)(εθ + εz) ,

σθ =E(1 − ν)

(1 + ν)(1 − 2ν)εθ +

(1 + ν)(1 − 2ν)(εz + εr) ,

σz =E(1 − ν)

(1 + ν)(1 − 2ν)εz +

(1 + ν)(1 − 2ν)(εr + εθ) ,

τzr =E

2(1 + ν)γzr.

(4.75)

Datorita numarului de necunoscute mai mic, modelul cu ele-mente finite asociat problemelor cu simetrie axiala este mai simpludecat cel prezentat ın §4.2. De exemplu, vectorul-coloana

ε(e)

©are

acum structura ε(e)

©=hε(e)

r , ε(e)θ , ε(e)

z , γ(e)zr

iT. (4.76)

Vectorul-coloanaσ(e)

©se prezinta si el sub o forma corespunzator

redusa: σ(e)

©=hσ(e)

r , σ(e)θ , σ(e)

z , τ (e)zr

iT. (4.77)

Matricea constitutiva elastica asociata problemelor cu simetrieaxiala are atunci structura (vezi (4.75))

C(e)=

E(e)

1 + ν(e)

266666664 1−ν(e)

1−2ν(e)ν(e)

1−2ν(e)ν(e)

1−2ν(e) 0

ν(e)

1−2ν(e)1−ν(e)

1−2ν(e)ν(e)

1−2ν(e) 0

ν(e)

1−2ν(e)ν(e)

1−2ν(e)1−ν(e)

1−2ν(e) 0

0 0 0 12

377777775 . (4.78)

Page 141: Metoda Elementelor Finite

Rezolvarea unei probleme de elasticitate 135

Sistemul de ecuatii care se obtine prin aplicarea MEF are tot oforma de tipul (4.56). Trebuie sa retinem ınsa faptul ca ın vectorul-coloana d nu mai sunt stocate decat deplasarile ur si wz alenodurilor:

d = [ur,1, wz,1, ur,2, wz,2, . . . , ur,nnod, wz,nnod

]T . (4.79)

Fiindca deplasarea vθ este identic nula (vezi (4.73)), includerea eiın modelul cu elemente finite devine inutila.

Sa retinem asadar ca problemele cu simetrie axiala pot fi re-zolvate eficient discretizand o sectiune axiala oarecare a corpuluistudiat. Totusi, utilizatorul programului MEF trebuie sa aiba ınvedere simetria axiala a domeniului de analiza atunci cand impuneconditiile la limita de tip deplasare. Luand ca exemplu corpul repre-zentat ın figura 4.4, constatam ca particulele aflate initial pe axa ztrebuie sa ramana pe aceasta pe toata durata deformarii (deci elevor avea deplasarea ur nula).

In practica, inginerii mecanici ıntalnesc foarte des problemecu simetrie axiala. Aproape ın orice echipament exista piese derevolutie supuse unor constrangeri de miscare si ıncarcari simetricrepartizate ın jurul unei axe. Un exemplu ın acest sens ıl reprezintacilindrii hidraulici solicitati de presiunea lichidului de lucru.

4.4 Exemplu de utilizare a metodei ele-

mentelor finite la rezolvarea unei

probleme de elasticitate

Consideram cazul unei placi dreptunghiulare de grosime con-stanta, ıncastrata pe una din laturi si supusa actiunii unei fortepunctuale P = 5000 N (fig. 4.5). Placa este confectionata dintr-unaliaj pe baza de aluminiu avand urmatoarele constante elastice:E = 0,72 × 1011 Nm−2, ν = 0,33 ≈ 1/3. In ceea ce urmeaza vomfolosi metoda elementelor finite pentru a determina tensiunile sideformatiile induse de forta P .

Page 142: Metoda Elementelor Finite

136 Capitolul 4

2ℓ

PGros. h

P = 5000 NE = 0,72 × 1011 Nm−2

ν = 0,33 ≈ 1/3ℓ = 1 mh = 0,1 m

Figura 4.5: Placa dreptunghiulara ıncastrata pe una din laturi sisupusa unei ıncarcari punctuale

Prin caracteristicile sale (geometria si dimensiunile corpului, na-tura si repartitia ıncarcarilor, modul de rezemare etc.), problemaenuntata anterior este asimilabila unei stari plane de tensiuni (vezi§4.3.1). Putem deci adopta un model ın a carui structura intervinmarimi dependente de numai doua coordonate. Acestea din urmacorespund axelor x si y din planul median al placii. In figura 4.6 estereprezentata grafic discretizarea domeniului de analiza. Se observafaptul ca amplasarea axelor x si y a fost aleasa de asa maniera ıncatexprimarea coordonatelor nodale sa fie cat mai convenabila. Dato-rita faptului ca urmarim elaborarea unui model pur ilustrativ sinicidecum obtinerea unei solutii numerice foarte precise, am folosito retea formata din numai doua elemente triunghiulare de ordinul

Page 143: Metoda Elementelor Finite

Rezolvarea unei probleme de elasticitate 137

P

x

y

31

2 4

(1)

(2)

Figura 4.6: Modelul cu elemente finite al placii dreptunghiulare(stare plana de tensiuni)

Tabelul 4.1: Matricea de conexiuni a retelei din figura 4.6

Element Nod Nod Nod(e) 1 2 3

(1) 1 3 2(2) 2 3 4

Tabelul 4.2: Coordonatele nodurilor din figura 4.6

Nod xk yk

k [m] [m]

1 0 02 0 ℓ = 13 2ℓ = 2 04 2ℓ = 2 ℓ = 1

Page 144: Metoda Elementelor Finite

138 Capitolul 4

ıntai (vezi §3.1.2 si §3.2.1). Va rezulta astfel un sistem de ecuatii cuun gabarit suficient de redus pentru a fi rezolvabil manual. Tabelele4.1 - 4.2 definesc matricea de conexiuni a retelei din figura 4.6 sicoordonatele celor patru noduri asociate.

Sistemul de ecuatii care descrie raspunsul elastic al modeluluirezulta prin particularizarea expresiei generale (4.56). In structuraacestui sistem, vectorul-coloana d stocheaza cate doua compo-nente de tip deplasare (uk, vk) pentru fiecare nod k = 1, . . . , 4 (vezi§3.4.2, mai ales exemplul (3.57)):

d = [u1, v1, u2, v2, u3, v3, u4, v4]T . (4.80)

Fiindca setul de ıncarcari exterioare este limitat la o forta con-centrata P (fig. 4.5 - 4.6), vectorul-coloana F are forma (vezi(4.58) si (4.60))

F = P = [0, 0, 0, 0, 0, 0, 0, −P = −5000 N]T . (4.81)

Remarcam faptul ca F poseda o singura componenta nenula (co-respunzatoare fortei ce actioneaza pe verticala ın jos asupra nodului4 – fig. 4.6).

Matricea globala de rigiditate [K] ınsumeaza contributia fie-caruia din cele doua elemente ale retelei (vezi (4.57) si (4.59)):

[K] =2X

e=1

K(e)

, (4.82)

K(e)= h

ZS(e)

B(e)

T C(e)

B(e)

dx dy, e = 1, 2. (4.83)

In (4.83), notatia S(e) simbolizeaza subdomeniul plan asociat ele-mentului (e), e = 1, 2, iar h reprezinta grosimea placii (fig. 4.5).

Observam ca formula (4.83) impune evaluarea matricilorB(e)

,

e = 1, 2. Acestea din urma se obtin aplicand operatorul [g] –vezi (4.19) – asupra matricilor corespunzatoare

N (e)

. In cazul

starii plane de tensiuni, [g] degenereaza la o forma considerabil mairestransa:

Page 145: Metoda Elementelor Finite

Rezolvarea unei probleme de elasticitate 139

[g] =

26664 ∂∂ x

0

0 ∂∂ y

∂∂ y

∂∂ x

37775 . (4.84)

(4.84) retine exclusiv derivatele ∂ /∂ x si ∂ /∂ y pentru ca toate ne-cunoscutele problemei depind numai de coordonatele (x, y).

Structura matricilor de tipN (e)

a fost precizata ın §3.4.2 (vezi

mai ales exemplul (3.58)). Pentru cele doua elemente triunghiularedin figura 4.6, respectivele matrici se exprima ın felul urmator (vezisi tabelul 4.1):

N (1)2×8

=h

N(1)1

i2×2

,hN

(1)2

i2×2

,hN

(1)3

i2×2

, [0]2×2

,

N (2)2×8

=[0]2×2 ,

hN

(2)2

i2×2

,hN

(2)3

i2×2

,hN

(2)4

i2×2

,

(4.85)

unde (vezi exemplele (3.59))hN

(e)k

i2×2

=

24 N(e)k 0

0 N(e)k

35 (4.86)

sunt submatrici diagonale care grupeaza functii de forma N(e)k

definite ın maniera (3.19) - (3.20).Din (4.21) si (4.84) - (4.86) rezulta urmatoarele expresii ale ma-

tricilorB(e)

, e = 1, 2:

B(1)3×8

=h

B(1)1

i3×2

,hB

(1)2

i3×2

,hB

(1)3

i3×2

, [0]3×2

,

B(2)3×8

=[0]3×2 ,

hB

(2)2

i3×2

,hB

(2)3

i3×2

,hB

(2)4

i3×2

,

(4.87)

unde hB

(e)k

i=

2666664 ∂N(e)k

∂ x0

0∂N

(e)k

∂ y

∂N(e)k

∂ y

∂N(e)k

∂ x

3777775 . (4.88)

Page 146: Metoda Elementelor Finite

140 Capitolul 4

In cazul elementelor triunghiulare de ordinul ıntai, functiile deforma N

(e)k sunt dependente liniare de x si y (vezi (3.19) - (3.20)).

Drept consecinta, atat derivatele ∂ N(e)k /∂ x si ∂ N

(e)k /∂ y, cat si

matricileB(e)

sunt constante (vezi tabelele 4.3 - 4.4). In aceste

conditii, formula (4.83) devineK(e)

= h

B(e)

T C(e)

B(e)

ZS(e)

dx dy =

= hA(e)B(e)

T C(e)

B(e)

, e = 1, 2.

(4.89)

Folosind (4.89) si (4.66) ımpreuna cu parametrii dimensionali siconstantele elastice din figura 4.5, obtinem prin calcul direct ma-tricile elementale

K(e)

, e = 1, 2 (tab. 4.5). Acestea din urma,

prin ınlocuire ın formula (4.82), furnizeaza matricea de rigiditate aretelei:

[K] = 6,75 × 108×

×

2666666666666666666647 4 −4 −2 −3 −2 0 0

4 13 −2 −12 −2 −1 0 0

−4 −2 7 0 0 4 −3 −2

−2 −12 0 13 4 0 −2 −1

−3 −2 0 4 7 0 −4 −2

−2 −1 4 0 0 13 −2 −12

0 0 −3 −2 −4 −2 7 4

0 0 −2 −1 −2 −12 4 13

377777777777777777775 .(4.90)

Explicitarile (4.90) si (4.81) definesc complet sistemul de ecuatii(4.56). Dintre cele opt ecuatii ale acestuia, primele patru sunt re-dundante. Ele corespund nodurilor 1 si 2, ale caror deplasari suntobiectul restrictiilor (vezi figura 4.6)

u1 = 0, v1 = 0,

u2 = 0, v2 = 0.(4.91)

Page 147: Metoda Elementelor Finite

Rezolvarea unei probleme de elasticitate 141

Tabelul 4.3: Gradientii functiilor de forma asociate elementelor dinfigura 4.6 (vezi tabelele 4.1 - 4.2, precum si relatiile (3.19) - (3.20))

Element Gradientii functiilor de forma asociate(e)

∂N(1)1

∂ x= y3−y2

2A(1) = − 12ℓ

= -0,5 m−1

∂N(1)1

∂ y= x2−x3

2A(1) = −1ℓ

= -1 m−1

(1)∂N

(1)2

∂ x= y1−y3

2A(1) = 0

∂N(1)2

∂ y= x3−x1

2A(1) = 1ℓ

= 1 m−1

∂N(1)3

∂ x= y2−y1

2A(1) = 12ℓ

= 0,5 m−1

∂N(1)3

∂ y= x1−x2

2A(1) = 0

∂N(2)2

∂ x= y3−y4

2A(2) = − 12ℓ

= -0,5 m−1

∂N(2)2

∂ y= x4−x3

2A(2) = 0

(2)∂N

(2)3

∂ x= y4−y2

2A(2) = 0

∂N(2)3

∂ y= x2−x4

2A(2) = −1ℓ

= -1 m−1

∂N(2)4

∂ x= y2−y3

2A(2) = 12ℓ

= 0,5 m−1

∂N(2)4

∂ y= x3−x2

2A(2) = 1ℓ

= 1 m−1

Page 148: Metoda Elementelor Finite

142 Capitolul 4

Tabelul 4.4: MatricileB(e)

asociate elementelor din figura 4.6 (vezi

relatiile (4.87) - (4.88) si tabelul 4.3)

Element Matricea asociata(e) [B(e)]

(1)

26664 -0,5 0 0 0 0,5 0 0 0

0 -1 0 1 0 0 0 0

-1 -0,5 1 0 0 0,5 0 0

37775(2)

26664 0 0 -0,5 0 0 0 0,5 0

0 0 0 0 0 -1 0 1

0 0 0 -0,5 -1 0 1 0,5

37775Suprimand ecuatiile asociate constrangerilor (4.91), obtinem urma-torul sistem de ecuatii redus (vezi (4.56), (4.90) si (4.80) - (4.81)):

6,75 × 108

2666664 7 0 −4 −2

0 13 −2 −12

−4 −2 7 4

−2 −12 4 13

3777775 2666664 u3

v3

u4

v4

3777775 =

2666664 0

0

0

−5000

3777775 . (4.92)

Solutia acestuia poate fi determinata numeric, folosind metoda eli-minarii Gauss [3, 4, 5, 8, 9, 10]. Rezulta astfel deplasarile nodurilornerestrictionate:

u3 = −1,023 × 10−6 m, v3 = −5,875 × 10−6 m,

u4 = 1,518 × 10−6 m, v4 = −6,618 × 10−6 m.(4.93)

Page 149: Metoda Elementelor Finite

Rezolvarea unei probleme de elasticitate 143

Tabelul 4.5: MatricileK(e)

asociate elementelor din figura 4.6

(vezi relatiile (4.89) si (4.66), tabelul 4.4, respectiv constantele elas-tice si parametrii dimensionali din figura 4.5)

Element Matricea asociata(e) [K(e)]

(1) 6,75 × 108

2666666666666666647 4 −4 −2 −3 −2 0 0

4 13 −2 −12 −2 −1 0 0

−4 −2 4 0 0 2 0 0

−2 −12 0 12 2 0 0 0

−3 −2 0 2 3 0 0 0

−2 −1 2 0 0 1 0 0

0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0

377777777777777775(2) 6,75 × 108

2666666666666666640 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0

0 0 3 0 0 2 −3 −2

0 0 0 1 2 0 −2 −1

0 0 0 2 4 0 −4 −2

0 0 2 0 0 12 −2 −12

0 0 −3 −2 −4 −2 7 4

0 0 −2 −1 −2 −12 4 13

377777777777777775

Page 150: Metoda Elementelor Finite

144 Capitolul 4

Tabelul 4.6: Starea de tensiuni si deformatii asociata elementelordin figura 4.6 (vezi relatiile (4.63) - (4.65), (4.20), (4.23) si (4.66),tabelul 4.4, respectiv constantele elastice din figura 4.5)

Element Starea de tensiuni si deformatii asociata(e)

(1)

ε(1)x = −0,5115 × 10−6

ε(1)y = 0

γ(1)xy = −2,9375 × 10−6

ε(1)z = 0,2558 × 10−6

σ(1)x = −41431,5 Nm−2

σ(1)y = −13810,5 Nm−2

τ (1)xy = −79312,5 Nm−2

(2)

ε(2)x = 0,7590 × 10−6

ε(2)y = −0,7430 × 10−6

γ(2)xy = −0,7680 × 10−6

ε(2)z = −0,0800 × 10−6

σ(2)x = 41418,0 Nm−2

σ(2)y = −39690,0 Nm−2

τ (2)xy = −20736,0 Nm−2

Page 151: Metoda Elementelor Finite

Rezolvarea unei probleme de elasticitate 145

In actuala etapa a calculelor, d este complet determinat.Propriu-zis, ınlocuind valorile (4.91) si (4.93) ın (4.80), obtinem

d = 10−6×[0; 0; 0; 0; −1,023; −5,875; 1,518; −6,618]T . (4.94)

Vectorul-coloana d sub forma (4.94) si matricileB(e)

, e = 1, 2,

deja calculate (vezi tabelul 4.4) urmeaza a fi introduse ın relatia(4.20). Rezulta astfel deformatiile elementale planare stocate ınvectori-coloana

ε(e)

©, e = 1, 2, avand structura de tip (4.64). In

stadiul imediat urmator, legea lui Hooke (4.23) ne ajuta sa evaluamsi tensiunile elementale planare. Acestea din urma sunt grupate ınvectori-coloana

σ(e)

©, e = 1, 2, avand structura de tip (4.63). Ime-

diat ce dispunem de tensiunile planare, suntem ın masura sa cal-culam si deformatia pe grosime a celor doua elemente (vezi formula(4.65)). Tabelul 4.6 centralizeaza rezultatele acestei suite de calcule.

Datorita faptului caB(e)

, e = 1, 2, sunt matrici constante

(vezi tabelul 4.4), aproximarea elementala a deformatiilor si ten-siunilor va avea aceeasi particularitate. O asemenea caracteristicaeste proprie tuturor elementelor care descriu campul deplasarilorprin polinoame de gradul ıntai. In general, ordinul de aproximareal necunoscutelor derivate este mai mic decat cel atasat necunoscu-telor primare. Aceasta regula se aplica si elementelor de ordin supe-rior lui unu. Atunci cand polinomul care descrie deplasarile continetermeni neliniari, aproximanta deformatiilor si a tensiunilor, chiardiminuata ca ordin, ramane totusi variabila. Iata unul dintre mo-tivele, pentru care elementele de tip patrulater si hexadric pot ficonsiderate calitativ superioare elementelor liniare de tip triunghi,respectiv tetraedru. Se impune ınsa o mentiune. Constructia mate-matica a metodei elementelor finite nu garanteaza continuitatea latraversarea frontierelor interelementale decat pentru functiile careaproximeaza necunoscutele primare (deplasarile, ın cazul nostru).Altfel spus, distributia necunoscutelor derivate (deformatii si ten-siuni, de exemplu) ar fi ramas discontinua si daca am fi utilizatpolinoame de grad superior. Acest neajuns este remediabil prin uti-lizarea unor algoritmi de netezire a solutiei numerice [9].

Page 152: Metoda Elementelor Finite

146 Capitolul 4

Rezultatele din tabelul 4.6, asa putine cum sunt, dau locunei observatii interesante si deosebit de utile sub aspect practic.Remarcam diferenta semnificativa dintre starile de tensiuni sideformatii asociate celor doua elemente. Acest lucru este un in-diciu clar al obtinerii unei solutii imprecise. In general, atunci canddistributia necunoscutelor derivate se distinge prin salturi majorela traversarea unei frontiere interelementale, putem fi siguri ca, lanivelul acelei regiuni, rezultatele numerice sunt afectate de erori se-rioase. La ora actuala, exista programe de analiza cu elemente finitecare pun la dispozitia utilizatorului tehnica discretizarii adaptive.Un asemenea program porneste de la o solutie determinata pe oretea mai grosolana, ale carei erori sunt estimate folosind criteriuldiscrepantelor la traversarea frontierelor interelementale. Pe bazaacestei evaluari, se procedeaza la o rafinare locala a retelei. Rezulta-tele obtinute pe noul model vor servi la reevaluarea erorii numerice.Procedura este parcursa iterativ pana la determinarea unei solutiiacceptabile din punct de vedere al preciziei.

4.5 Concluzii

In acest capitol am prezentat structura generala a modeluluicu elemente finite asociat problemelor de elasticitate liniara. Aminsistat asupra unor aspecte foarte importante pentru utilizatoriiprogramelor MEF, cum ar fi identificarea clara a necunoscutelorsi specificarea corecta a conditiilor la limita de diferite tipuri. Deasemenea, am analizat cateva cazuri particulare de probleme bi-dimensionale de elasticitate care sunt frecvent ıntalnite ın practica.

Spatiul limitat nu ne-a permis sa intram ın prea multe detalii.Totusi, informatiile prezentate precum si exemplul numeric suntsuficiente pentru a furniza o imagine de ansamblu asupra moduluiın care metoda elementelor finite poate servi la rezolvarea unei clasede probleme care, prin natura lor, au cea mai mare importanta ındomeniul ingineriei mecanice.

Page 153: Metoda Elementelor Finite

Anexa A

Formule de derivarematriceala

In cadrul discutiei din §4.2 (mai precis, atunci cand am stu-diat procedura de minimizare a marimii Π), am ıntalnit problemaderivarii unor produse matriceale de tipul dT V , respectivdT [M ] d ın raport cu d. In expresiile anterioare, d si V sunt vectori-coloana, iar [M ] este o matrice simetrica.

Ne propunem sa deducem formule matriceale pentru calcululderivatelor mentionate mai sus.

Vom nota cu n dimensiunea vectorilor-coloana d si V , res-pectiv ordinul matricii [M ]. Facand legatura cu discutia din §4.2,n este egal cu triplul numarului total de noduri din retea (totusi,sa retinem faptul ca relatiile pe care le vom obtine mai jos au va-labilitate generala, indiferent de valoarea lui n). In aceste conditii,marimile de tip matriceal care apar ın produsele ce trebuie derivateau urmatoarea structura:

d = [d1, d2 , . . . , dn]T , (A.1)

V = [V1, V2, . . . , Vn]T , (A.2)

147

Page 154: Metoda Elementelor Finite

148 Anexa A

[M ] =

2666664 M11 M12 . . . M1n

M21 M22 . . . M2n

......

. . ....

Mn1 Mn2 . . . Mnn

3777775 . (A.3)

Conditia de simetrie a matricii [M ] se scrie astfel:

Mij = Mji, (∀)i, j = 1, 2, . . . , n. (A.4)

Incepem prin a ne ocupa de produsul dT V . Acesta poatefi explicitat algebric sub forma

dT V = d1V1 + d2V2 + . . . + dnVn. (A.5)

Din (A.5) rezulta printr-un calcul simplu

∂di

dT V

= Vi, (∀)i = 1, 2, . . . , n, (A.6)

astfel ıncat ajungem imediat la prima formula de derivare matri-ceala cautata:

∂ ddT V

= V . (A.7)

In mod asemanator, produsul matriceal dT [M ] d este ex-plicitabil sub urmatoarea forma algebrica:

dT [M ] d = d1 (M11d1 + M12d2 + . . . + M1ndn) +

+ d2 (M21d1 + M22d2 + . . . + M2ndn) +

...

+ dn (Mn1d1 + Mn2d2 + . . . + Mnndn) .

(A.8)

Page 155: Metoda Elementelor Finite

Formule de derivare matriceala 149

Din (A.8) rezulta

∂ di

dT [M ] d

= (Mi1d1 + Mi2d2 + . . . + Mindn) +

+ d1∂

∂ di

(M11d1 + M12d2 + . . . + M1ndn) +

+ d2∂

∂ di

(M21d1 + M22d2 + . . . + M2ndn) +

...

+ dn

∂ di

(Mn1d1 + Mn2d2 + . . . + Mnndn) ,

(∀)i = 1, 2, . . . , n,

(A.9)

sau, dupa derivarea tuturor termenilor din membrul drept al ega-litatii de mai sus,

∂ di

dT [M ] d

= (Mi1d1 + Mi2d2 + . . . + Mindn) +

+ (d1M1i + d2M2i + . . . + dnMni) ,

(∀)i = 1, 2, . . . , n.

(A.10)

Daca tinem cont de proprietatea de simetrie (A.4), putem re-organiza relatia (A.10) sub forma

∂ di

dT [M ] d

= 2 (Mi1d1 + Mi2d2 + . . . + Mindn) ,

(∀)i = 1, 2, . . . , n.

(A.11)

Ajungem ın felul acesta la a doua formula de derivare matricealacautata:

∂ ddT [M ] d

= 2 [M ] d . (A.12)

Page 156: Metoda Elementelor Finite

150 Anexa A

(A.7) si (A.12) ne-au servit ın §4.2 la calculul termenilor dinmembrul drept al egalitatii (4.51) – vezi setul de relatii (4.52) -(4.55).

Page 157: Metoda Elementelor Finite

Bibliografie

[1] Bia, C., Ille, V., Soare, M. Rezistenta materialelor si teoria

elasticitatii. Bucuresti: Editura Didactica si Pedagogica, 1983.

[2] Brebbia, C.A. The Boundary Element Method for Engineers.Londra: Pentech Press, 1978.

[3] Demidovitch, B.P., Maron, I.A. Elements de calcul numerique.

Moscova: Mir, 1979.

[4] Henwood, D., Bonet, J. Finite Elements. A Gentle Introduc-

tion. Londra: MacMillan, 1996.

[5] Hutton, D.V. Fundamentals of Finite Element Analysis. NewYork: McGraw-Hill, 2004.

[6] Iacob, C., Homentcovschi, D., Marcov, N., Nicolau, A. Mate-

matici clasice si moderne, vol. IV. Bucuresti: Editura Tehnica,1983.

[7] Lazar, I. Metoda elementelor de frontiera ın inginerie. Cluj-Napoca: Presa Universitara Clujeana, 1997.

[8] Pascariu, I. Elemente finite. Concepte si aplicatii. Bucuresti:Editura Militara, 1985.

[9] Segerlind, L.J. Applied Finite Element Analysis. New York:John Wiley, 1976.

151

Page 158: Metoda Elementelor Finite

152 Bibliografie

[10] Zienkiewicz, O.C., Taylor, R.L. The Finite Element Method,

vol. I. Londra: McGraw-Hill, 1989.

[11] Zienkiewicz, O.C., Taylor, R.L. The Finite Element Method,

vol. II. Londra: McGraw-Hill, 1991.