Lucrarea nr. 5

1

Toolbox-ul PDE

Este un instrument destinat studiului ecuaţiilor diferenţiale cu derivate parţiale (PDE) bidimensionale. Toolbox-ul conţine un set de funcţii specializate precum si o interfaţă grafică (GUI) intuitivă care permit definirea si rezolvarea unei clase largi de probleme întâlnite în ingineria electrică. Aceste probleme pot fi formulate pe domenii cu geometrii complicate.

Cu toolboxul PDE se rezolvă ecuaţia eliptică (cu necunoscuta u):

fauuc pe un domeniu bidimensional . Condiţiile pe frontiera reprezintă valori ale soluţiei u

sau ale derivatelor sale. Alte ecuaţii care pot fi rezolvate sunt:

fauuctud - ecuaţie parabolică

fauuctud

2

2

- ecuaţie hiperbolică

)()()( ufuuauuc - ecuaţie eliptică neliniară

Coeficienţii c, a, f, d pot fi funcţii scalare sau complexe definite pe domeniul , care pot depinde de: coordonatele x şi y soluţia u

derivatele yu ,

xu ale soluţiei u

timp t (pt. ecuaţii parabolice şi hiperbolice)

Condiţiile pe frontieră pot fi :

Dirichlet: hu = r pe Neumann generalizată: gquucn

Pentru rezolvarea unei probleme se parcurg următorii paşi:

Desenarea domeniului Stabilirea condiţiilor pe frontiere Definirea ecuaţiilor Generarea reţelei de discretizare (mesh) Rezolvarea ecuaţiilor în nodurile reţelei de discretizare Reprezentarea grafică a rezultatelor

Lucrarea nr. 5

2

Există câteva aplicaţii predefinite: Electrostatics, Magnetostatics, AC Power Electromagnetics, Conductive Media DC, Heat Transfer, Diffusion.

Exemplul 1. Rezolvarea unei probleme de electrostatică: ρVε unde

r0 εεε ; ( F/m108.854ε 120

), pe un domeniu limitat la interior de o suprafaţă patratică cu latura de 0.25 m şi limitat la exterior de o suprafaţă pătratică având latura de 0.5m. Nu există sarcină electrică în domeniul specificat. Pe frontiera interioară potenţialul electric scalar este 1000 V, iar pe frontiera exterioară potenţialul este 0.

Se apelează Toolboxul PDE din fereastra de comenzi Matlab: >> pdetool Options-Application-Generic Scalar Se desenează două pătrate R1, R2 cu latura de 0.25 respectiv 0.5m. Pentru stabilirea domeniului de calcul se introduce formula R2-R1.

Fereastra pentru stabilirea dimensiunilor lui R1.

Lucrarea nr. 5

3

Se introduc condiţiile pe frontierele domeniului de calcul: Boundary – Specify Boundary

Conditions ...

Pe frontiera interioară se introduce condiţia Dirichlet V=1000, iar pe frontiera exterioară

condiţia Dirichlet V=0.

Lucrarea nr. 5

4

Ecuaţia de rezolvat: 0V Se introduc coeficienţii ecuaţiei ca în fereastra următoare PDE – PDE Specification...

Se discretizează domeniul de calcul Mesh – Initialize Mesh şi apoi Refine Mesh.

Lucrarea nr. 5

5

Se stabileşte modul de rezolvare adaptiv: Solve – Parameters...

Se rezolvă ecuaţia Solve – Solve PDE. Se reprezintă grafic soluţia obţinută.

Lucrarea nr. 5

6

Exemplul 2: Modelarea unui electromagnet – problemă neliniară Un electromagnet simplu este alcătuit dintr-un miez magnetic înconjurat de o bobină. Se

consideră un miez lung şi se neglijează efectele de capăt. Domeniul de calcul se reduce la o secţiune plană. Ecuaţia câmpului electromagnetic exprimată în potenţial magnet vector este:

01

JA

unde J este densitatea de curent - pemeabilitate magnetică (funcţie neliniară de A în cazul materialelor feromagnetice)

205.015000200

A

>> pdetool Options - Application-Magnetostatics

Options – Grid Spacing

Lucrarea nr. 5

7

Se desenează dreptunghiurile R1 (0.2 0.05), R2 (0.2 0.1) şi R3(1 1)

Se fixează condiţia pe frontieră A = 0

Lucrarea nr. 5

8

PDE – Show Subdomain Labels

Regiunea 1 – miezul feromagnetic Regiunile 2, 3 – părţi ale bobinei Regiunea 4 – aer Se introduc separat coeficienţii ecuaţiei pentru fiecare regiune. Regiunea 1:

Coeficientul mu are valoarea: 200+5000./(1+0.05*(ux.^2+uy.^2))

Lucrarea nr. 5

9

Regiunea 2:

Regiunea 3:

Regiunea 4:

Se discretizează domeniul de calcul Mesh – Initialize Mesh şi apoi Refine Mesh. Se stabileşte modul de rezolvare adaptiv: Solve – Parameters... Se rezolvă ecuaţia Solve – Solve PDE. Se reprezintă grafic soluţia obţinută.

Lucrarea nr. 5

10

Zoom

Lucrarea nr. 5

11

Dacă se inversează sensul lui J în regiunile 2 şi 3 se obţine:

În orice moment al rezolvării unei probleme se pot exporta rezultatele din GUI în spaţiul de

lucru unde se poate beneficia de întreaga flexibilitate a mediului MATLAB şi de unde se pot accesa funcţiile care formează Toolboxul PDE.

Lucrarea nr. 5

12

>>pdetool- Se deseneaza geometria- Se introduc conditiile pe frontiere

Boundary - Export Decomposed Geometry, Boundary Cond's...

Matricea geometrieig

Geometrie M-file

Conditii pe frontiereM-file

Matricea conditiilorpe frontiere

wbound

b

wgeom

[p,e,t]=initmesh(g)p,e,t - datele retelei de discretizare

u=assempde(b,p,e,t,c,a,f)-div(c*grad(u))+a*u=f

pdeplot(p,e,t,'Property',Value)

Matrici coeficienti CoeficientiM-file

Lucrarea nr. 5

13

Aplicaţii 1. Modelul numeric al unui traductor de nivel de tip capacitiv Sunt cazuri în care acurateţea procesului de măsurare sau de urmărire a evoluţiei în timp a

nivelului lichidului într-un recipient cu ajutorul unui traductor de tip capacitiv presupun cunoaşterea cu precizie suficient de ridicată din punct de vedere tehnic a dependenţei dintre mărimile de ieşire şi de intrare ale traductorului, în cazul de faţă dintre capacitate şi nivel. Acest lucru este posibil doar folosind un model numeric pentru determinarea câmpului electric, respectiv a capacităţii condensatorului, deoarece un astfel de model este în măsură să considere configuraţia geometrică concretă şi proprietăţile fizice reale.

Această aplicaţie prezintă un model numeric în element finit, dezvoltat în Matlab, al condensatorului unui traductor pentru măsurarea nivelului unui lichid cu proprietăţi electroizolante. Formularea problemei Să se construiască modelul numeric în element finit al unui traductor de nivel de tip capacitiv, având structura unui condensator plan cu dielectric mobil, în vederea evaluării caracteristicii statice a traductorului – capacitatea electrică între bornele (a) şi (b), Fig..1, în funcţie de nivelul lichidului, C(y) şi a sensibilităţii s = dC/dy.

x

z

a

y

g

O

A

l

Fig. 1 Modelul fizic al traductorului de nivel

Date geometrice Aplicaţiile aferente modelului numeric au următoarele date geometrice comune:

- distanţa dintre armăturile condensatorului, a = 40 mm;

Lucrarea nr. 5

14

- grosimea armăturilor, g = 5 mm; - înălţimea armăturilor, l = 300 mm; Poziţia de referinţă pentru construcţia geometriei este aceea în care nivelul lichidului se

află la 60 mm faţă de limita inferioară a armăturilor condensatorului.

Modelul de câmp Extinderea armăturilor condensatorului în direcţia Oz, Fig. 1, fiind mult mai mare decât

dimensiunile în planul yOz, modelul pentru calculul câmpului electrostatic atunci când există o diferenţă de potenţial între plăci este unul 2D plan-paralel. Aceasta înseamnă că în evaluarea mărimilor globale precum de exemplu energia câmpului electric nu se iau în considerare contribuţia regiunii din afara spaţiului paralelipipedic dintre plăci şi scăderea câmpului către marginile perpendiculare pe axa Oz în raport cu zona centrală.

Fig. 2. Domeniul de calcul al câmpului electrostatic

Modelul matematic al câmpului electrostatic exprimat în potenţial electric scalar este reprezentat de ecuaţia:

div(·gradV) = 0 ,

unde este permitivitatea electrică, iar V este mărimea de stare asociată câmpului electrostatic – potenţialul electric scalar. Modelul numeric prezentat în continuare presupune determinarea necunoscutei V(x,z) prin metoda elementelor finite atunci când se cunosc potenţialele celor două armături metalice,

Lucrarea nr. 5

15

evaluarea intensităţii câmpului electric E = - gradV, a energiei câmpului electrostatic şi apoi a capacităţii traductorului.

Definirea domeniului de calcul al câmpului electrostatic Domeniul de calcul din Fig. 2 exploatează simetria câmpului plan-paralel (acelaşi în orice plan paralel cu xOy, Fig. 1) în raport cu axa Oz. El conţine următoarele două regiuni: - AER, delimitată de conturul ABCDEFGH; - LICHID, delimitată de conturul ABKLEFMN Conturul BCDELK al armăturii metalice este o frontieră echipotenţială a domeniului de calcul, ca orice suprafaţă a unui corp metalic în regim electrostatic; interiorul acestui contur nu face parte din domeniul de calcul, ca orice secţiune într-un corp metalic, în care câmpul electric este nul în regim electrostatic. Domeniul de calcul este mărginit de linia MN , care face parte din axa de simetrie a câmpului plan-paralel şi de conturul HGMN, situat suficient de departe de zona de câmp electric intens dintre armături astfel încât aproximaţia de câmp nul să fie acceptată. Acest contur modelează de fapt frontiera de la infinit a domeniului de existenţă a câmpului electrostatic atunci când în întreg spaţiul se află numai condensatorul.

Construirea geometriei şi a reţelei de elemente finite Se reprezintă în Toolboxul PDE trei dreptunghiuri: - R1 – corespunzător conturului AFGH, Fig. 2; - R2 – corespunzător conturului AFMN; - R3 – corespunzător conturului CDLK.

Domeniul de calcul din Fig. 3 corespunde operaţiei (R1+R2)-R3. Matricea geometriei se exportă din PDE, iar pe structura ei se scrie geometria parametrizată. Reţeaua de discretizare a domeniului de calcul în elemente finite, Fig. 4, se realizează cu funcţia PDE initmesh care construieşte elemente finite triunghiulare pe baza algoritmului Delaunay.

Lucrarea nr. 5

16

Fig. 3. Geometria Matlab a domeniului de calcul

Fig. 4. Reţeaua de elemente finite a domeniului de calcul

Lucrarea nr. 5

17

Proprietăţi fizice, condiţii de frontieră Permitivitatea relativă a materialului asociat regiunii LICHID este r = 80, iar a regiunii

AER este 1. Referitor la condiţiile la limită, variabila de stare potenţial electric scalar are următoarele

valori: - V = 100 V, pe conturul BCDELKB; - V = 0 V, pe frontiera HAN.

Conturul HGFMN este o frontiera de tip Neuman omogen, unde derivata normală a potenţialului

electric este nulă, 0nV

; altfel spus, liniile câmpului electric sunt tangente la acest contur, sau

valoarea locală a fluxul electric este nulă. Modelul numeric al caracteristicii statice a traductorului Linii echipotenţiale şi harta în degrade de culori a potenţialului electrostatic sunt reprezentate în Fig. 5, iar harta modulului vectorului intensitate a câmpului electric în Fig. 6.

Fig. 5. Linii echipotenţiale şi harta potenţialului electric

Lucrarea nr. 5

18

Fig. 6. Harta intensităţii câmpului electric Capacitatea C a condensatorului se determină din egalitatea:

2CUdV

2εE 2

V

2

care exprimă energia câmpului electrostatic ca integrală pe volumul domeniului de calcul a densităţii de volum E2/2, respectiv în funcţie de capacitate şi tensiunea U = 200 V a tensiunii la bornele condensatorului. Rezultatul numeric este C = 1,62 nF. Pentru comparaţie, în aproximaţia nelijării şi a efectului de capăt în raport cu coordonata y, Fig. 1, sunt două condensatoare ideale (cu câmp uniform între armături) în paralel, cu valori diferite ale permitivităţii electrice, a căror capacitate echivalentă este Cid = 1,12 nF. Dependenţa capacităţii condensatorului de nivelul lichidului, C(y), Fig. 7, este rezultatul unui studiul parametric al modelul numeric. Se observă o lejeră neliniaritate a caracteristicii pentru valori mici ale variabilei y. În rest, respectiv în domeniul y (50, 300) mm, sensibilitatea traductorului are o valoare practic constantă s = dC/dy = 0,017 nF/mm.

Lucrarea nr. 5

19

Fig. 7. Dependenţa capacităţii de nivelul lichidului

Lucrarea nr. 5

20

2. Model numeric pentru studiul unui traductor de proximitate inductiv de tip bobină simplă

Cel mai simplu traductor de tip inductiv pentru evidenţierea vecinătăţii sau pentru măsurarea distanţei faţă de un corp metalic se bazează pe modificarea parametrilor echivalenţi ai unei bobine parcurse de curent alternativ, în principal a inductivităţii acestei bobine în funcţie de poziţia relativă traductor – corp metalic. Această modificare este rezultatul prezenţei curenţilor induşi în corpul din proximitatea bobinei a căror intensitatea depinde între altele de poziţia relativă dintre bobina inductor şi corpul indus, respectiv de cuplajul electromagnetic inductor – indus. În ipoteza că nici-un alt corp conductor şi/sau magnetic nu se află în imediata vecinătate a ansamblului bobină – corp proxim, parametrii principali ai cuplajului inductor - indus într-un traductor inductiv tip bobină simplă sunt următorii:

- configuraţia geometrică şi numărul de spire ale bobinei; - configuraţia geometrică a corpului proxim; - rezistivitatea şi permeabilitatea – mai precis, proprietăţile de magnetizare ale corpului

proxim; - poziţia relativă bobină – corp; - frecvenţa de alimentare a bobinei.

În aplicaţia de faţă se dezvoltă un model numeric destinat studiului modelui fizic simplu al traductorului de proximitate din Fig. 1, în care bobina are formă toroidală, are contur circular şi secţiune dreptunghiulară, corpul din proximitate este semispaţiul z < zA, iar planul transversal de simetrie al bobinei este paralel cu suprafaţa semispaţiului.

Fig. 1. Modelul fizic al unui traductor de proximitate tip bobină simplă

Chiar şi pentru configuraţia geometrică simplă din Fig.1, dar cu atât mai mult în cazul unor configuraţii geometrice neregulate sau al corpurilor magnetic sau electric neliniare, un model numeric este singurul care poate să ofere de exemplu caracteristica statică a traductorului, respectiv variaţia inductivităţii în funcţie de poziţia relativă bobină – corp, cu acurateţea impusă de precizia măsurării electrice a mărimilor neelectrice.

z

O

A

I . ~ r

.

Lucrarea nr. 5

21

Formularea problemei Să se construiască modelul numeric în element finit destinat studiului unui traductor de proximitate cu curenţi induşi tip bobină simplă de formă toroidală, paralelă cu suprafaţă unei plăci din oţel nemagnetic de grosime considerabilă. Să se determine dependenţa inductivităţii bobinei – inductor şi a sensibilităţii traductorului în funcţie de distanţa la care aceasta se află faţă de suprafaţa plană a plăcii - indus. Să se studieze modificarea caracteristicilor traductorului atunci când placa este din oţel magnetic.

Date geometrice şi electrice Bobina traductorului se caracterizează prin următoarele dimensiuni:

- diametrul interior, di = 50 mm; - diametrul exterior, de = 70 mm; - grosimea bobinei, g = 5 mm;

În poziţia de referinţă pentru construcţia geometriei, distanţa dintre bobină şi suprafaţă plană a plăcii, Fig. 1, este OA = 30 mm. Numărul de spire al bobinei este N = 200.

Modelul de câmp al traductorului Modelul fizic al traductorului de proximitate se caracterizează prin simetrie de rotaţie şi

ca urmare, modelul de câmp asociat este unul 2D axisimetric în coordonate cilindrice (r, z). Densitatea curentului inductor în bobină, J1 – sursa câmpului electromagnetic şi potenţialul magnetic vector A – mărimea de stare necunoscută a câmpului electromagnetic, au orientare azimutală, perpendiculară pe planul rOz, Fig. 1.

Modelul matematic al câmpului electromagnetic cvasistaţionar armonic, exprimat prin potenţialul magnetic vector complex A, este reprezentat de ecuaţia:

rot[(1/ rot A] + jA/ = J1

unde = 2f, f este frecvenţa, este rezistivitatea şi este permeabilitatea magnetică. Termenul al doilea din partea stângă a ecuaţiei este densitatea curenţilor induşi, J2, nenulă doar în regiunea plăcii indus (unde are valoare finită). Densitatea de curent J1 este diferită de zero doar în cele două regiuni (suprafeţe dreptunghiulare) care reprezintă secţiunea axială prin bobina inductor. Bobina este realizată din conductor filiform, deci aceste regiuni sunt de tipul fără curenţi induşi, motiv pentru care modelul de material care se asociază unor astfel de regiuni este neconductor (echivalent cu în ecuaţia de mai sus) şi nemagnetic.

Modelul numeric prezentat în continuare, dezvoltat în Matlab, presupune determinarea prin metoda elementelor finite a necunoscutei A(r, z).

Definirea domeniului de calcul al câmpului electromagnetic Domeniul de calcul al structurii 2D axisimetrice a câmpului electromagnetic din Fig. 2

conţine următoarele regiuni, cu proprietăţi fizice distincte: - BOBINA, interiorul conturului EFGH , regiune care modelează bobina din conductor

de tip filiform, fără curenţi induşi, în care se cunoaşte densitatea J1 a curentului inductor;

- PLACA, interiorul conturului ABCD, regiune conductoare şi nemagnetică cu curenţi induşi.

- AER, regiune neconductoare şi nemagnetică.

Lucrarea nr. 5

22

Conturul domeniului de calcul este format din: - linia KD, care aparţine axei de simetrie Oz, Fig. 1; - conturul KLCD, care modelează infinitul regiunii de existenţă a câmpului

electromagnetic. În orice punct al acestui contur câmpul electromagnetic are o valoare neglijabilă în raport cu zona de câmp intens (de exemplu, inducţia magnetică are o valoare neglijabilă în raport cu valoarea maximă pe conturul EFGH). Cât priveşte segmentul DC, condiţia de câmp neglijabil este asigurată dacă acesta se află la o disatanţă AD = BC 3, unde este adâncimea de pătrundere a câmpului în placă.

Fig. 2. Domeniul de calcul al câmpului electromagnetic

Construirea geometriei şi a reţelei de elemente finite În prima etapă se reprezintă în Toolboxul PDE trei dreptunghiuri: - R1 – corespunzător conturului EFGH, Fig. 2; - R2 – corespunzător conturului ABCD; - R3 – corespunzător conturului KLBA.

Domeniul de calcul corespunde reuniunii R1+R2+R3. Matricea geometriei se exportă din PDE, iar pe structura ei se scrie geometria sub formă matriceală, complet parametrizată. Reţeaua de discretizare se realizează cu funcţia PDE initmesh, care realizează divizarea domeniului în elemente finite triunghiulare, Fig.3, pe baza algoritmului Delaunay.

Lucrarea nr. 5

23

Fig. 3. Reţeaua de elemente finite

Proprietăţi fizice, condiţii de frontieră

Câmpul electromagnetic se determină pentru: - frecvenţa f = 20 kHz; - densitatea curentului în bobină J1 = 0,1 A/mm2

Rezistivitatea oţelului nemagnetic este = 10-6 m, iar permeabilitatea sa magnetică este aceeaşi cu valoarea corespunzătoare întregului domeniu de calcul, 0 = 410-7 H/m (r = 1). Pentru aplicaţiile în care placa este din oţel magnetic se consideră aceeaşi valoare a rezistivităţii şi model neliniar de magnetizare cu permeabilitatea magnetică relativă iniţială ri = 1000. Fiind o aplicaţie de semnal, câmpul magnetic este foarte slab şi adâncimea de pătrundere are valoarea corespunzătoare lui ri. Variabila de stare potenţial magnetic vector satisface o condiţia A = 0 pe întreg conturul KLCDK al domeniului de calcul. Semnificaţia fizică este aceea că axa de simetrie într-o problemă 2D axisimetrică este linie de câmp, iar pe restul conturului, care modelează infinitul, valoarea locală a fluxului magnetic este nulă (componenta normală a inducţiei magnetice este nulă). Modelul numeric al traductorului Liniile câmpului magnetic şi harta inducţiei magnetice pentru placa nemagnetică sunt reprezentate în Fig. 4 şi 5 pentru OA = 30 mm.

Lucrarea nr. 5

24

Fig. 4 Liniile câmpului magnetic

Fig. 5 Harta inducţiei magnetice

Lucrarea nr. 5

25

În Fig. 6 se prezintă harta densităţii J2 a curenţilor induşi în regiunea PLACA pentru cazurile oţel nemagnetic, respectiv oţel magnetic.

a)

b)

Fig. 6. Harta densităţii curenţilor induşi în regiunea PLACA a – placă nemagnetică b – placă magnetică

Lucrarea nr. 5

26

Inductivitatea L a bobinei se determină din egalitatea 2

dV2

2

V

2 LIB

, care exprimă

energia câmpului magnetic ca integrală pe volumul domeniului de calcul a densităţii de volum B2/2, respectiv în funcţie de inductivitate si curent prin bobină.

Modelul numeric permite efectuarea studiului parametric al dependenţei inductivităţii L a bobinei de distanţa OA dintre bobină şi placă, distanţă notată cu z.

Fig. 7. Variaţia inductivităţii cu distanţa bobină – placă

a – placă nemagnetică b – placă magnetică

Lucrarea nr. 5

27

Fig.8. Variaţia sensibilităţii cu distanţa bobină – placă