Electrotehnica - INTRODUCERE LA METODA DIFERENTELOR FINITE

W(M) fiind o functie scalara de punct, iar - un camp de vectori ales astfel incat rotorul sau sa coincida cu campul densitatii de curent

Se observa ca, de fapt, exista numai doua formulari distincte: folosind un potential scalar sau unul vectorial.

Tratarea potentialului scalar va fi exemplificata in cazul campului electrostatic.

Potentialul vector, folosit in cazul campului magnetic stationar (produs de curenti), creeaza anumite dificultati la problemele tridimensionale si aici va fi exemplificat numai pentru camp in 2D.

Folosind potentialul electric scalar, se poate stabili ecuatia cu derivate partiale a potentialului electric

In cazul potentialului magnetic vector se obtine ecuatia

unde cu n = 1/µ s-a notat reluctivitatea, functie de H sau B si de punct. In mediu liniar si omogen (cel putin pe portiuni), cu conditia de 'etalonare' se obtine ecuatia lui Poisson in forma vectoriala

Atunci cand problema este in 2D (plan-paralela sau plan-radiala) exista un sistem de referinta in care campul are o singura componenta, care depinde numai de coordonatele din planul perpendicular pe acesta componenta: Atunci solutia se descrie intr-un plan perpendicular pe si ecuatia vectoriala se reduce la una scalara ( , are numai o componenta, paralela cu

Pentru ca problema de camp stationar sa fie formulata complet, adica sa aiba solutie unica, mai trebuie precizate conditiile la limita, care pot fi date sub formele cunoscute: Dirichlet, Neumann sau combinata (Robin).

2.2. Metoda diferentelor finite

Metoda diferentelor finite stabileste un sistem de ecuatii algebrice care aproximeaza solutia (de regula functia potential) printr-un numar finit de valori in puncte discrete ale domeniului de camp, reprezentand noduri ale unei retele uni-, bi- sau tri- dimensionale, dupa cum problema este in 1D, in 2D sau in 3D. Caracterul liniar sau neliniar al sistemului algebric obtinut este dat de relatia constitutiva.

Metoda diferentelor finite implica alegerea unei retele, in ale carei noduri vor fi definite valorile functiei potential cu care se descrie solutia. Reteaua multidimensionala se alege, de regula, ortogonala. Ea poate fi carteziana, polara, sferica, sau chiar curbilinie. Aici se va aborda numai cazul retelei carteziene, cea mai des folosita.

Reteaua pote avea pas neuniform si se recomanda ca ea sa fie astfel construita incat frontiera domeniului si liniile sau suprafetele de discontinuitate a proprietatilor de material sa fie formate din laturi si/sau diagonale in 2D, respectiv din fete si/sau patrulatere formate din muchii si diagonale in 3D. Reteaua cu pas constant dupa fiecare directie (o coordonata) asigura o mai buna precizie de aproximare prin diferente si o convergenta mai buna a procesului de calcul. Ecuatiile asociate metodei diferentelor finite se pot stabili pe doua cai: aproximand ecuatia diferentiala prin diferente (divizate) sau aproximand formele integrale ale ecuatiilor campului.

2.3. Aproximarea ecuatiei diferentiale

A) Se considera intai medii liniare si omogene (cel putin pe portiuni).

La problema unidimensionala, trebuie aproximata solutia ecuatiei

cu f(x) cunoscuta. Se aleg ca necunoscute valorile potentialului V in punctele discrete x_k, k = 1,2,n, si se noteaza V(x_k)=V_k. Intai se considera pasul de discretizare h constant

Daca solutia poate fi dezvoltata in serie Taylor in jurul fiecarui punct, se stabilesc relatiile

si atunci

Astfel se obtine ecuatia discreta corespunzatoare nodului k

si ecuatii similare pentru toate nodurile.

Relatia obtinuta este citata ca 'ecuatia in trei puncte' si este asociata 'punctului central' de indice k

Daca reteaua nu are pas constant, se noteaza h_k = x_k+1 - x_k si relatiile se modifica astfel

si apoi

Se observa ca in cazul retelei cu pas neuniform eroarea de aproximare este mai mare, intrucat restul O(h_m) variaza liniar cu un pas mediu h_m si nu patratic, ca in cazul pasului constant.

Pentru probleme multidimensionale se procedeaza in acelsi mod dupa fiecare directie (dimensiune), folosind retele de discretizare independente. Astfel, de exemplu, pentru o retea carteziana bidimensionala, cu pasii constanti h_x, h_y si indicii k, i, se obtine aproximarea

cu care se stabileste imediat ecuatia discreta corespunzatoare unui nod (k,i). Expresia obtinuta este citata ca 'ecuatie in 5 puncte' si este asociata punctului central (k,i

Generalizarea pentru cazul tridimensional este evidenta, obtinand 'o ecuatie in 7 puncte', asociata unui punct central (k,i,j

Expresiile de mai sus se simplifica daca pasul retelelor are aceeasi valoare dupa toate coordonatele h_x = h_y = h_z = h, fara a afecta precizia de aproximare.

Pentru retele cu pas variabil relatiile se stabilesc in mod similar.

In 'punctele centrale' (k,i,j) situate pe suprafete de discontinuitate a proprietatilor de material nu se mai pot folosi relatiile de mai sus, ci trebuie scrise ecuatii speciale, care exprima relatiile de trecere intre medii diferite: conservarea componentei normale a inductiei electrice si a componentei tangentiale a intensitatii campului electric; ultima conditie este echivalenta cu continuitatea valorii potentialului independent de mediu. In legatura cu exprimarea componentei normale a inductiei apar anumite dificultati, care se rezolva folosind pentru aceste componente aproximarea prin diferente inainte, respectiv inapoi. Aceste dificultati dispar la metoda aproximarii formei integrale a ecuatiilor.

B) In medii neomogene sau neliniare, proprietatea mediului care intervine in relatia constitutiva variaza o data cu pozitia punctului in camp si trebuie aproximata prin diferente o ecuatie diferentiala de forma

Pentru ca operatia sa fie consecventa, derivatele partiale din interiorul parantezelor ar trebui evaluate in acele puncte in care se da si proprietatea e (sau µ, n la medii magnetice). Acest fapt determina complicatii suplimentare deosebite in aproximarea ecuatiei diferentiale. Cum se va arata mai departe, metoda aproximarii formei integrale satisface cu usurinta aceasta cerinta, ca si cele legate de simularea conditiilor la limita, fapt pentru care se recomanda folosirea ultimei metode.

2.4. Aproximarea formei integrale a ecuatiilor

Se pleaca de la forma ecuatiilor campului exprimata cu ajutorul unui potential, scalar sau vector; in ultimul caz aici se va aborda numai cazul campului in 2D.

2.4.1. Campul electric

Pentru campul electric, ecuatiile se prezinta sub forma integrala

Prima integrala se poate evalua cu o formula de cvadratura prin diferente. Folosind expresia (1-28) cu diferente centrale pentru o componenta E_x a intensitatii campului, se obtine

sau, mai concis

cu o aproximatie de ordinul , independent de natura, omogenitatea si liniaritatea mediului.

Fiind data o retea rectangulara si valorile potentialului V in nodurile retelei, in acest mod se pot determina componentele intensitatii campului electric de-a lungul tuturor laturilor, la mijloacele laturilor (v. fig.2-1).

Cunoscand valorile permitivitatii, care poate avea valori diferite pe ochiuri in 2D sau pe hexaedre in 3D, se pot determina inductiile electrice corespunzatoare.

Fie problema de camp in 2D, in coordonatele (x,y). Se considera conturul rectangular G_k,i care trece prin mijloacele laturilor concurente in nodul (k,i). In modul descris mai sus se determina Se noteaza cu e_k,i permitivitatea ochiului (k,k+1,i,i+1), cu r_k,i densitatea de suprafata a sarcinii pe acest ochi si cu q_k,i sarcina punctiforma din nodul (k,i). Aplicand conturului G_k,i legea fluxului electric (a doua integrala din 2-19) si estimand integrala pe fiecare latura cu primul termen al relatiei (1-26), se obtine ecuatia asociata nodului central (i,k

In cazul unei retele cu pas constant, respectiv in medii omogene, se obtin relatii mult mai simple.

In ecuatia (2-22) nu au fost explicitate valorile intensi­tatilor campului, dar ele se exprima prin relatii de forma (2-21) si deriva din 5 valori ale potentialului: din cel al nodului central si din cele ale nodurilor de la extermitatile laturilor concurente in nodul central.

Ecuatia obtinuta (2-22) este remarcabila prin faptul ca ea este valabila si pentru puncte centrale (k,i) situate pe suprafete de discontinuitate a proprietatilor de material, respectiv pentru medii in care proprietatile de material sunt functiune de punct (inclusiv in medii neliniare), deci nu mai trebuie folosite si alte ecuatii 'speciale' pentru aceste zone. Cum se va vedea mai departe, o parte a aceleiasi ecuatii va fi folosita si la frontierele cu conditii Neumann.

Pentru probleme de camp in 3D, in coordonatele (x,y,z) se scrie legea fluxului electric pe un hexaedru S_k,i,j centrat pe nodul (k,i,j), care trece prin mijloacele laturilor concurente in nodul central (k,i,j). La intersectia laturilor retelei cu fetele hexaedrului se calculeaza 6 componente ale intensitatii campului electric, perpendiculare pe fetele respective. In cele 8 hexaedre in care este impartit hexaedrul principal prin suprafetele de coordonate care trec prin nodul central, se definesc 8 conduc­tivitati, cu ajutorul carora se determina inductiile electrice (6*4=24 valori) si apoi se calculeaza fluxurile electrice ale fetelor. In sfarsit se pot determina si sarcinile punctuale, de volum si de suprafata care intervin in membrul drept.

Ecuatia astfel obtinuta are aceleasi proprietati ca cea in 2D, dar are 24 termeni in membrul stang si se exprima cu ajutorul potentialelor a 7 noduri: cel al nodului central si cele ale nodurilor de la extremitatile laturilor concurente in nodul central.

Metoda aproximarii formei integrale a ecuatiilor are deci doua etape, corespunzatoare celor doua ecuatii din (2-19):

- in prima se calculeaza la mijlocul fiecarei laturi a retelei componenta intensitatii campului electric din lungul acelei laturi, cu care se pot calcula componentele longitudinale ale inductiilor electrice din toate diedrele adiacente;

- in a doua etapa se calculeaza fluxul electric prin supra­fete care trec perpendicular prin mijloacele laturilor retelei in jurul cate unui nod 'central'. Pentru membrul drept se calculeaza sarcina electrica continuta de aceste suprafete.

2.4.2. Campul magnetic

Pentru campul magnetic, in 2D, potentialul vector are o singura componenta - fie ea dupa axa z - si atunci

Campul magnetic va fi cuprins in planul z = const. Se noteaza cu un versor in planul z = const, perpendicular pe elementul de arc , astfel incat produsul vectorial sa fie omoparalel cu . Atunci ecuatiile campului magnetic in forma integrala sunt

Prima integrala se calculeaza cu primul termen al formulei de cvadratura (1-28) si pentru componenta B_h a campului, perpendiculara pe segmentul x_k,k+1, da rezultatul

sau mai concis

cu o aproximatie de ordinul independent de natura, omogenitatea si liniaritatea mediului.

Fiind data o retea rectangulara si valorile potentialului A in nodurile retelei, in acest mod se pot determina componentele inductiei magnetice perpendiculare pe laturi, la mijloacele laturilor (v. fig.2-2).

Cunoscand valorile reluctivitatii, care poate avea valori diferite pe cele patru ochiuri adiacente, se pot determina inten­sitatile corespunzatoare ale campului magnetic.

Fie problema de camp in 2D, in coordonatele (x,y). Se consi­dera conturul rectangular G_k,i care trece prin mijloacele laturilor concurente in nodul (k,i). In modul descris mai sus se determina Se noteaza cu A_k,i reluctivitatea ochiului (k,k+1,i,i+1), cu J_k,i densitatea de curent pe suprafata acestui ochi si cu i_k,i curentul conductorului filiform din nodul (k,i). Aplicand conturului G_k,i teorema lui Ampčre (a doua integrala din 2-24) si estimand integrala pe fiecare latura cu primul termen al relatiei (1-26), se obtine ecuatia corespunzatoare nodului central (i,k

In cazul unei retele cu pas constant, respectiv in medii omogene, se obtin relatii mult mai simple.

In ecuatia (2-27) nu au fost explicitate valorile intensitatilor campului, dar ele se exprima prin relatii de forma (2-26) si deriva din 5 valori ale potentialului: din cel al nodului central si din cele ale nodurilor de la extermitatile laturilor concurente in nodul central.

Ecuatia obtinuta (2-27) este remarcabila prin faptul ca ea este valabila si pentru puncte centrale (k,i) situate pe suprafete de discontinuitate a proprietatilor de material, respectiv pentru medii in care proprietatile de material sunt functiune de punct (inclusiv in medii neliniare), deci nu mai trebuie folosite si alte ecuatii 'speciale' pentru aceste zone. Cum se va vedea mai departe, o parte a aceleiasi ecuatii va fi folosita si la frontierele cu conditii Neumann.

Metoda aproximarii formei integrale a ecuatiilor are deci doua etape, corespunzatoare celor doua ecuatii din (2-24):

- in prima se calculeaza la mijlocul fiecarei laturi a retelei componenta inductiei magnetice perpendiculara pe acea latura, cu care se pot calcula componentele transversale ale intensitatilor campului magnetic in cele doua diedre adiacente;

- in a doua etapa se calculeaza circulatia intensitatii campului magnetic (tensiunea magnetomotoare) pe conturul care trece prin mijloacele laturilor retelei in jurul cate unui nod 'central'. Pentru membrul drept se calculeaza curentul electric care trece prin suprafata sprijinita pe acest contur.

Nota. Pentru un camp magnetic in 3D s-ar putea proceda in mod asemanator, cu observatia ca in acest caz potentialul vector are trei componente, care trebuie sa indeplineasca o conditie de etalonare (de obicei divergenta nula). Atunci a doua ecuatie din (2-24) ar trebui scrisa pe trei contururi, dar pentru a nu obtine un sistem cu mai multe ecuatii decat necunoscute, conditia de etalonare trebuie folosita ca relatie de eliminare etc, ceea ce complica foarte mult aplicarea. Din aceasta cauza la problemele de camp magnetic in 3D se prefera metoda T-W sau similara, in care apare numai o functiune scalara (W) ca necunoscuta.

2.5. Simularea conditiilor la limita

Conditiile la limita sunt de trei feluri: Dirichlet, Neumann sau mixte. Ultimele apar mai rar la campurile electrice si magnetice stationare.

Daca conditia Dirichlet - cu potential dat - este pusa in puncte care coincid cu noduri ale retelei, nodurile respective au potentialul cunoscut, pentru aceste noduri nu se scrie ecuatia in care ar aparea ca nod central, iar termenii celorlalte ecuatii care contin potentialele cunoscute trec in membrul drept al sistemului (cu semn schimbat).

Daca conditia Dirichlet este data in puncte care nu coincid cu noduri ale retelei, se pot stabili relatii de legatura cu potentialele nodurilor vecine, folosind, de regula, interpolari liniare. Prin introducerea acestor noi relatii se strica simetria sistemului de ecuatii, ceea ce are influente importante asupra metodei de rezolare si asupra spatiului de stocare necesar. Din aceste motive se evita aceasta imprejurare, alegand o retea adaptata frontierei.

Nodurile aflate pe frontiera cu conditie Neumann au poten­tiale variabile (necunoscute a priori).

Daca conditia Neumann - adica cu derivata dupa normala data - este in noduri ale retelei, ea poate fi luata in consideratie cu usurinta la metoda aproximarii formei integrale a ecuatiilor, intrucat marimea impusa reprezinta fie intensitatea unui camp electric, fie intensitatea unui camp magnetic (cand se foloseste un potential magnetic scalar), fie o inductie magnetica, marimi care intervin direct in forma integrala a ecuatiilor si se iau in consideratie fara nici o dificultate. In fig. 2-3 s-a figurat o retea si trei zone (A, B, C) in care se dau conditii la limita tip Neumann pe segmentele de frontiera marcate cu sageti.

De exemplu, pentru un domeniu cu camp electric, daca in nodul 6 din zona A se da E V x, atunci in ecuatia asociata nodului central 6 se vor include numai termenii corespunzatori domeniului de camp (nodurile 1, 6, 7, 11, cu ariile si permitivitatile din ochiurile 1-2-7-6 si 6-7-12-11), iar la membrul drept se va aduna fluxul electric care intra prin portiunea de frontiera cuprinsa intre nodul central 6 si mijloacele laturilor 1-6 si 6-11, adica E e e )h_y/2, permiti­vitatile fiind numerotate dupa nodul stanga-jos al ochiului.

In zona B, ecuatia nodului central 19 va cuprinde termenii corespunzatori nodurilor 14, 18, 19, 20, cu ariile si permitivi­tatile din ochiurile 18-13-14-19 si 19-14-15-20, iar la membrul drept se aduna fluxul electric care intra prin portiunea de frontiera cuprinsa intre nodul central 19 si mijloacele laturilor 18-19 si 19-20, adica -E₁₉(e e )h_x/2, unde E₁₉ = V y in nodul 19.

In zona de colt C, ecuatia nodului central 16 va include numai termenii corespunzatori nodurilor 11, 12, 16, 17, cu aria si permitivitatea ochiului 11-12-17-16, iar la membrul drept se aduna fluxul electric care intra prin portiunea de frontiera cuprinsa intre nodul 16 si mijloacele laturilor 11-16 si 16-17.

Daca conditia tip Neumann este pusa unui domeniu de camp magnetic stationar, in care se foloseste potentialul magnetic vector, ecuatia nodului aflat pe frontiera corespunde unei integrale de contur, care se efectueaza numai pe portiunea de contur aflata in domeniul de camp si se completeaza cu tensiunea magnetica stabilita cu conditia la limita. Deci pentru zona A se aduna in membrul drept B A A )h_y, unde B₆ = A x in nodul 6, iar pentru zona B se aduna in membrul drept B₁₉(A A )h_x, unde B₁₉ = A y in nodul 19. Pentru zona C vor interveni doua inductii tangentiale etc.

Sistemul de ecuatii Si rezolvarea

Necunoscutele problemei de camp descrisa cu ajutorul dife­rentelor finite sunt potentialele nodurilor interioare si ale nodurilor de pe portiunile de frontiera cu conditii Neumann. Pentru fiecare nod cu potential necunoscut se scrie cate o ecuatie, corespunzatoare cazului in care acel nod este considerat central. Se obtine, astfel, un sistem algebric liniar cu atatea ecuatii cate necunoscute. In cazul cand intreaga frontiera este numai cu conditii Neumann, trebuie fixat potentialul unui nod oarecare (conditie Dirichlet), altfel solutia nu este unica si pot aparea dificultati de calcul numeric.

Daca toate conditiile la limita sunt precizate in noduri ale retelei, atunci sistemul de ecuatii algebrice obtinut este simetric si pozitiv definit. Pentru rezolvarea sa se pot folosi metode directe (eliminare Gauss, factorizare Choleschi etc) sau iterative (SOR = Suprarexare succesiva, gradient conjugat etc). Cea mai simpla din punct de vedere al programarii este SOR.

Fie sistemul de ecuatii prezentat sub forma

Se demonstreza ca pentru sisteme de ecuatii stabilite cu metoda diferentelor finite matricea N este pozitiv definita si are raza spectrala subunitara, adica permite un calcul iterativ convergent. Relatia de calcul iterativ in forma Jacobi este

Efectuarea unei iteratii comporta multiplicarea vectorului solutie V^(k) de la iteratia k cu matricea N si adunarea cu termenul liber b', adica calcule simple. Totusi metoda necesita rezervarea de spatiu pentru doi vectori solutie, V^(k) si V^(k+1), iar apoi transferul celui de al doilea vector in primul, pentru reluarea calculelor la pasul iterativ urmator.

Componentele noului vector solutie se calculeaza succesiv. Componenta de ordin p a solutiei este calculata facand produsul scalar al solutiei existente cu linia p a matricei N

Pentru eliminarea vectorului solutie suplimentar V^(k+1) noua componenta se poate stoca in spatiul rezervat vechiului vector. Astfel se obtine schema de calcul Gauss-Seidel, care se poate prezenta sub o forma matriceala.

Fie matricea M descompusa in forma

in care L si U sunt matrici strict inferior si superior triunghiulare. Atunci procesul iterativ Gauss-Seidel se descrie prin relatia

Intrucat componentele ale vectorului solutie sunt 'noi' prima suma trece in membrul stang si relatia devine

Si in acest caz iteratiile sunt convergente, iar matricea de iterare (- (L+D)^-1 U) are o raza spectrala mai mica decat matricea de iterare Jacobi, adica convergenta este mai rapida.

Relatia de iterare se mai pune sub forma

Norma vectorului variatie de la o iteratie la alta DV^(k+1) este o masura a convergentei. Cand norma este nula s-a atins solutia 'exacta'. Daca convergenta este in progresie geometrica atunci 'eroarea' fata de solutie se estimeaza prin raportul

In relatia (2-34) se poate folosi orice norma. Este simpla si cu semnificatie bine precizata norma max(abs), care determina maximul valorii absolute a componentelor vectorului variatie. Convergenta se imbunatateste prin suprarelaxare succesiva (SOR = successive over-relaxation), cu relatia

unde w este factorul de relaxare, cu valori intre 0 si 2; pentru w I (1,2) este suprarelaxare (iteratiile sunt divergente la w ³ 2), iar pentru w I (0,1) este subrelaxare. Ultima se foloseste la stabilizarea unor procese iterative neliniare, altfel instabile (divergente).

Pentru suprarelaxare in cazul sistemelor rezultate din metoda diferentelor finite se lucreaza cu w = 1,71,95 cu valori mai mari la sisteme mari. Exista o teorie pentru determinarea factorului de suprarelaxare optim w_o, care asigura cea mai rapida convergenta.

Rata de convergenta asimptotica (la numar mare de iteratii) este in relatie directa cu raza spectrala r a matricii de iterare. Intr-adevar, fiecare iteratie multiplica vectorul eroare cu matricea de iterare. Prezentand matricea de iterare sub forma sa modala, rezulta ca la iteratia de ordinul p valorile proprii se prezinta ridicate la puterea p. Astfel dupa un numar mare de iteratii va ramane predominanta numai cea mai mare valoare proprie - raza spectrala. Mai departe procesul iterativ continua dupa o progresie geometrica, cu ratia r

Fie dw w w_o abaterea factorului de relaxare de la valoarea optima. La apropierea de optim prin valori mai mici, dw < 0 sau w < w_o, rata de convergenta variaza mult cu dw, pe cand la depasirea valorii optime dw > 0 se obtine o variatie mult mai lenta si aproape proportionala cu dw a ratei de convergenta, iar solutia iterativa incepe sa oscileze in jurul solutiei exacte.

Solutia se considera aproximata cu eroarea e data de (2-34). Atunci cand eroarea scade sub o valoare tolerata, se considera atinsa solutia 'satisfacatoare', adica este rezolvat sistemul.

Document Info

Accesari: 4572
Apreciat:

Comenteaza documentul:

Nu esti inregistrat
Trebuie sa fii utilizator inregistrat pentru a putea comenta

Creaza cont nou