UNIVERSITATEA POLITEHNICA TIMIŞOARA

FACULTATEA DE ELECTRONICĂ ŞI TELECOMUNICAŢII

REFERAT DE DOCTORAT NR. 3

MODELARE ŞI SIMULARE

Coordonator ştiinţific: prof.dr.ing. Ignea Alimpie Drd. as.ing. Iftode Cora

2008

Cap. 1 Rezolvarea ecuaţiilor lui Maxwell prin metoda diferenţelor finite în timp (FDTD)

Introducere

Metoda FDTD este folosită tot mai mult în analiza interacţiunii câmpului electromagnetic cu materia (analiza circuitelor şi antenelor microstrip, a ghidurilor de undă, analiza efectelor biologice ale câmpului electromagnetic, etc) Faţă de alte metode numerice, tehnica diferenţelor finite se bucură de câteva avantaje, între care:

Ecuaţiile lui Maxwell pot fi discretizate foarte uşor. Ecuaţiile discretizate sunt rezolvate în manieră secvenţială, fiind uşor de implementat pe un

sistem numeric de calcul. Metoda se poate aplica problemelor ce implică structuri complexe, care ar fi foarte greu de

rezolvat pe cale analitică sau prin alte metode numerice. Fiind o metodă de analiză în timp, oferă în mod direct soluţiile de regim tranzitoriu, putând

apoi transforma uşor soluţia din domeniul timp în domeniul frecvenţă. În acest fel se poate obţine comportarea în frecvenţă a diverselor dispozitive.

Metoda poate fi aplicată pentru analiza mediilor neomogene, cu pierderi, anizotrope, dispersive şi cu parametri variabili în timp, ceea ce-i oferă un avantaj net faţă de alte metode numerice care nu au un spectru de utilizare atât de larg.

Pentru simulările în care regiunea analizată are întindere infinită, se utilizează condiţii de absorbţie la frontieră (engl. ABC - Absorbing Boundary Condition), impuse planelor care delimitează regiunea modelată, plane care apar din cauza trunchierii spaţiului.

Desigur, metoda prezintă şi dezavantaje, care însă pot fi minimizate. Principalele dezavantaje sunt:

• La capetele domeniului analizat pot apărea false reflexii ale câmpului (din cauza algoritmului)

• Dispersie numerică a soluţiei poate apărea ca urmare a faptului că discretizarea nu se face in paşi infinit mici.

Totuşi, dezavantajul principal al metodei îl constituie, în cazul problemelor complexe, timpul mare de calcul şi necesitatea de a dispune de spaţiu de memorare mare.

Tehnica FDTD poate fi adaptată pentru modelarea cât mai fidelă a mediilor reale. Una dintre modificările importante aduse algoritmului iniţial, dezvoltat de K. Yee [1], a fost implementată de către Luebbers ş.a. [2] pentru mediile în care parametrii dielectrici depind de frecvenţă, cum este cazul celor biologice. Metoda propusă de Luebbers, denumită (FD)2-TD, este o îmbunătăţire a algoritmului FDTD, în care se ţine cont de dependenţa spatială de frecvenţă. Se poate determina astfel, cu precizie bună, răspunsul unui mediu la semnale de bandă largă, fără a creşte semnificativ timpul de calcul. Metoda a fost folosită [3] pentru a determina curentul indus şi distribuţia SAR într-un model eterogen al corpului uman atunci când este expus unui impuls electromagnetic de bandă largă. Alte lucrări [4] au folosit diferenţele finite în timp pentru a modela fantome de corp uman care sunt expuse iradierii cu diferite câmpuri. Parametrii pentru fantoma de corp uman sunt recomandaţi de standardele internaţionale [5]. Modelul utilizat este

un cilindru de plastic, cu diametru 0.305 m şi înaltime de 1.7 m umplut cu soluţie salină având permitivitatea 77

rε = şi conductivitatea 0.7 /S mσ = . Rezultatele simulării sunt în concordanţă

cu valorile obţinute experimental. Pentru a modela cât mai fidel mediul real în cazul în care se impune trunchierea acestuia (limitarea simulării la un domeniu restrâns din spaţiu) J.P. Berenger [6] a modificat algoritmul Yee [1] ţinând cont de condiţiile de absorbţie la frontieră, introducând conceptul de strat perfect adaptat (engl. PML – Perfectly Matched Layer). Definirea, în acest fel, a condiţiilor de frontieră pentru domeniile cu întindere (teoretic) infinită a permis atât îmbunătăţirea rezultatelor cât şi a timpului de calcul la analiza interacţiunii între câmpul electromagnetic produs de telefoanele celulare şi modelele de cap uman [7]. Tot din cauza utilizării pe scară largă a telefoanelor mobile, problema interferenţei electromagnetice între diverse echipamente a devenit tot mai importantă. Fiind tot o problemă de câmp, se poate trata folosind metoda diferenţelor finite în timp. În [8] se analizează interferenţa produsă de telefoanele mobile asupra stimulatoarelor cardiace. Şi în acest caz rezultatele numerice sunt în concordanţă cu cele obţinute pe cale experimentală. În [9] s-au publicat rezultatele recente ale unui studiu privind efectul câmpului electromagnetic emis de un telefon mobil asupra unui fetus de 6 luni. Modelul geometric al mamei, fetusului şi interacţiunea acestora cu câmpul radiat de antena telefonului au fost analizate utilizând metoda FDTD.

1.1. Scurtă descriere a metodei FDTD

Metoda diferenţelor finite în domeniul timp (FDTD) permite modelarea numerică a interacţiunii dintre câmpul electromagnetic şi materie şi se bazează pe rezolvarea ecuaţiilor lui Maxwell în forma diferenţială. În esenţă, metoda constă în înlocuirea ecuaţiilor diferenţiale cu ecuaţii cu diferenţe finite. Cea mai cunoscută implementare a FDTD este algoritmul Yee [1], în care spaţiul este discretizat într-o retea cubică, în nodurile căreia se memorează proprietăţile dielectrice şi valoarea câmpului electric/magnetic. Derivatele spaţiale şi temporale sunt înlocuite cu diferenţe centrale. Pentru actualizarea câmpului electromagnetic în fiecare punct din spaţiu, se foloseşte un algoritm de tip leap-frog (întreţeserea în timp a componentelor E şi H). Se cunoaşte excitaţia aplicată la momentul iniţial în spaţiul 3D, iar apoi se calculează evoluţia în timp a câmpului electromagnetic.

Într-o analiză FDTD de câmp trebuie avute în vedere în special următoarele aspecte: • Precizia rezultatelor: spunem despre o metodă că este precisă, dacă soluţia numerică este

foarte apropiată de soluţia exactă. În general sunt 3 tipuri de erori care afectează acurateţea unui rezultat: erorile de modelare, erorile de discretizare şi erorile de rotunjire [10].

• Dispersia numerică: apare din cauza discretizării spaţiului. Discretizarea are ca efect apariţia de componente spectrale noi, care se propagă cu viteze diferite.

• Stabilitatea algoritmului: o soluţie numerică este stabilă dacă produce un rezultat mărginit pentru o variaţie mărginita a intrării. Dacă notăm cu nξ , eroarea la momentul n, eroarea la

momentul n+1 este: 1n ngξ ξ+ = , unde g se numeşte factor de amplificare (a erorii). Pentru ca

soluţia numerică să fie stabilă este necesar ca 1n nξ ξ+ ≤ , altfel spus 1g ≤ .

Precizia, dispersia şi stabilitatea numerică sunt influenţate de rezoluţiile de eşantionare spaţială şi temporală. Rezoluţia spaţială trebuie să fie mai mică decât 0.1 minλ ( minλ fiind

lungimea de undă minimă ce se propagă în domeniul analizat). Rezoluţia temporală se alege ulterior, astfel încat să fie asigurată stabilitatea numerică a algoritmului. Cerinţele de memorie pentru FDTD sunt proporţionale cu N - numărul de celule din domeniul de calcul.

1.2 Ecuaţiile lui Maxwell în 3 dimensiuni

Se consideră o regiune din spaţiu în care mediul material poate absorbi energie electrică sau magnetică. Ecuaţiile lui Maxwell în forma diferentială sunt:

Legea lui Faraday: t

∂= −∇ × −

∂

BE M (1)

Legea lui Ampere: t

∂= ∇× −

∂

DH J (2)

Legea lui Gauss pentru câmpul electric:v

ρ∇ ⋅ =D (3)

Legea lui Gauss pentru câmpul magnetic: 0∇⋅ =B (4)

Unde: E - intensitatea câmpului electric [V/m] H - intensitatea câmpului magnetic [A/m]

D - inducţia electrică [C/m2] B - inducţia magnetică [Wb/m2]

M - densitatea de curent magnetic echivalentă [V/m2] J - densitatea de curent electric [A/m2]

vρ -densitatea de sarcină [C/m3]

În mediile liniare, izotrope şi nedispersive sunt valabile următoarele relaţii de proporţionalitate: 0r

ε ε ε= =D E E (5)

0rµ µ µ= =B H H (6)

-120 8.85 x 10 F/mε = - permitivitatea aerului,

-60 1.256 x 10 H/mµ = - permeabilitatea aerului,

rε si r

µ - permitivitatea, respectiv permeabilitatea relativă.

J si M pot conţine componente care acţionează ca surse independente pentru câmpurile E şi H (se vor nota în continuare

sourceJ şi

sourceM ). Dacă există pierderi prin transformarea energiei

câmpurilor E şi H în căldură, rezultă:

sourceσ= +J J E ;

source σ ∗= +M M H (7)

unde σ -conductivitatea electrică [S/m] σ ∗ -pierderi magnetice echivalente [Ω/m]

Înlocuind (5) şi (6) în (1) şi (2), obţinem:

1 1( )

sourceE

tσ

µ µ∗∂

= − ∇× − +∂

HM H (8)

1 1

( )sourcet

σε ε

∂= − ∇ × − +

∂

EH J E (9)

Se rescriu ecuaţiile (8) şi (9), folosind proiecţiile pe axele de coordonate:

1[ ( )]

x

yx zsource x

EH EM H

t z yσ

µ∗

∂∂ ∂= − − +

∂ ∂ ∂ (10a)

1

[ ( )]y

y xzsource y

H EEM H

t x zσ

µ∗

∂ ∂∂= − − +

∂ ∂ ∂ (10b)

1

[ ( )]z

yxzsource z

EEHM H

t y xσ

µ∗

∂∂∂= − − +

∂ ∂ ∂ (10c)

1[ ( )]

x

yx zsource x

HE HJ E

t y zσ

ε

∂∂ ∂= − − +

∂ ∂ ∂ (11a)

1

[ ( )]y

y x zsource y

E H HJ E

t z xσ

ε

∂ ∂ ∂= − − +

∂ ∂ ∂ (11b)

1

[ ( )]z

y xzsource z

H HEJ E

t x yσ

ε

∂ ∂∂= − − +

∂ ∂ ∂ (11c)

Ecuaţiile (10) şi (11) stau la baza algoritmului numeric FDTD pentru descrierea interacţiunilor dintre undele electromagnetice şi mediile considerate.

Soluţiile ecuaţiilor lui Maxwell sunt adesea imposibil de exprimat analitic. K. Yee a dezvoltat o metodă numerică de determinare a soluţiilor ecuaţiilor, în cazul în care condiţiile de frontieră sunt cele pentru un conductor perfect [1].

1.3 Algoritmul lui Yee

Metoda numerică dezvoltată de Yee oferă posibilitatea determinării în timp şi spaţiu atât a câmpului electric cât şi a câmpului magnetic, prin rezolvarea celor două ecuaţii de rotor. Spaţiul în care se efectuează analiza de câmp este împărţit în celule cubice. După cum se observă în figura 1, componentele E şi H sunt astfel poziţionate în spaţiul tridimensional, încât fiecare componentă de câmp electric să fie înconjurată de 4 componente de câmp magnetic, şi fiecare componentă de câmp magnetic să fie înconjurată de 4 componente de câmp electric. În spaţiul propus de Yee (o reţea 3D), continuitatea componentelor tangenţiale ale lui E şi ale lui H este păstrată la interfaţa a două medii cu proprietăţi distincte, atâta timp cât interfaţa este paralelă cu una dintre axele de coordonate ale reţelei. Prin urmare, nu mai trebuie specificate condiţiile de frontieră la interfaţă, dar trebuie specificate permitivitatea şi permeabilitatea în fiecare punct în care se calculează câmpul. Deşi algoritmul Yee foloseşte doar ecuaţiile de rotor ale lui Maxwell, soluţiile pe care acesta le furnizează satisfac implicit celelalte două ecuaţii (relaţiile (3) şi (4) ) datorită modului în care sunt localizate componentele lui E şi H în reţeaua Yee şi a operaţiilor cu diferenţe centrale efectuate [11].

Fig 1. Plasarea vectorilor de câmp electric şi magnetic într-o celula Yee cubică. Componentele H sunt situate la mijlocul fiecărei laturi a cubului, iar componentele E, în centrul fiecărei feţe a cubului.

Componentele vectorilor E şi H sunt calculate în timp în maniera leap-frog, aşa cum se poate vedea în fig. 2. Astfel, valoarea E , într-un punct din spaţiu şi la un moment dat de timp, se calculează şi stochează în memorie utilizând valoarea lui E în acelaşi punct, dar la momentul precedent, precum şi valorile lui H la momentul precendent, în puncte adiacente celui de interes. Valoarea lui H se calculează apoi în funcţie de valoarea sa în acelaşi punct, dar la momentul precedent, şi de valorile lui E (anterior calculate) în punctele adiacente, iar rezultatul se stochează în memorie. Ciclul continuă cu recalcularea valorilor E pentru momentul de timp imediat următor, pe baza valorilor E şi H anterior determinate, până la momentul de timp dorit, şamd.

Fig. 2 Diagrama spaţio-temporală folosită pentru determinarea E şi H într-un spatiu 1-D. Se observă utilizarea diferenţelor centrale pentru derivatele spaţiale şi metoda leap-frog pentru derivatele temporale.

H

E E E E

E E E E

E E E E

H H

H H H

0x = 2

xx

∆= x x= ∆

3

2

xx

∆= 2x x= ∆

5

2

xx

∆= 3x x= ∆

0t =

2

tt

∆=

t t= ∆

3

2

tt

∆=

2t t= ∆

Condiţii la frontieră În cazul unui conductor perfect, condiţia de frontieră este anularea componentei tangente

a câmpului electric la suprafaţa conductoare. Această condiţie implică şi anularea componentei normale a câmpului magnetic la suprafaţa conductoare. Suprafaţa conductoare este aproximată printr-un set de feţe ale unor cuburi, paralele cu axele de coordonate. Feţele perpendiculare pe axa x se vor alege astfel încât să conţină puncte în care sunt definite

yE şi z

E (vezi fig.1). În

mod similar se aleg planele perpendiculare pe celelalte axe. Ecuaţiile cu diferenţe finite folosite pentru derivatele temporale sunt de tip “diferenţe

centrale” şi produc erori de ordinul 2, atât în timp cât şi în spaţiu (v. Anexa 1 – fdtd_1d.m). Algoritmul nu este disipativ temporal, adică undele numerice ce se propagă în reţea nu se atenuează, în timp, din cauza metodei numerice folosite.

1.4 Notaţii pentru diferenţe finite

Unui punct din spaţiul discretizat, de coordonate ( , , )i j k , i se va asocia cantitatea ( , , )i x j y k z∆ ∆ ∆ O funcţie dependentă de timp şi spaţiu se va nota astfel:

.

, ,( , , , )not

n

i j ku i x j y k z n t u∆ ∆ ∆ ∆ = (12)

unde: t∆ - incrementul temporal, presupus uniform în intervalul de observaţie, iar n –nr întreg

x∆ , y∆ , z∆ - incremenţii spaţiali pe cele trei direcţii ale sistemului de coordonate, iar i, j, k -numere întregi.

Astfel, dacă considerăm derivata lui u în raport cu x, evaluată la momentul n

t n t= ∆ ,

aceasta se va scrie astfel, folosind ecuaţiile cu diferenţe finite:

1/ 2, , 1/ 2, , 2 2 2 2( , , , ) [( ) , ( ) , ( ) , ( ) ]n n

i j k i j ku uui x j y k z n t O x y z t

x x

+ −−∂∆ ∆ ∆ ∆ = + ∆ ∆ ∆ ∆

∂ ∆ (13a)

Incrementul de 1/ 2 la indicele (pe axa ) denotă o diferenţă spaţială de 1/ 2i x x± ± ∆ . Yee a ales această notaţie pentru a intercala componentele E şi H în spaţiu la intervale de / 2x∆ . Similar se scriu derivatele parţiale ale lui u în funcţie de y, respectiv de z. Derivata lui u în raport cu timpul, evaluată într-un punct dat din spaţiu, (i, j, k), este:

1/ 2 1/ 2, , , , 2 2 2 2( , , , ) [( ) , ( ) , ( ) , ( ) ]n n

i j k i j ku uui x j y k z n t O x y z t

t t

+ −−∂∆ ∆ ∆ ∆ = + ∆ ∆ ∆ ∆

∂ ∆ (13b)

1.4.1 Ecuaţiile lui Maxwell în trei dimensiuni, scrise cu diferenţe finite

Ecuaţia (11a) se va scrie, folosind diferenţele finite, astfel:

1/2 1/2, 1/2, 1/2 , 1/2, 1/2 , 1, 1/2 , , 1/2 , 1/2, 1 , 1/2,

, 1/2, 1/2

, 1/2, 1/2 , 1/2, 1/2 , 1/2, 1/2

1(

)x

n n n n n n

x i j k x i j k z i j k z i j k y i j k y i j k

i j k

n n

source i j k i j k x i j k

E E H H H H

t y z

J E

ε

σ

+ −+ + + + + + + + + +

+ +

+ + + + + +

− − −= −

∆ ∆ ∆

− − ⋅

(14)

Prin ecuaţia (14) s-a dorit evaluarea componentei Ex a câmpului electric la momentul n,

în punctul de coordinate (i, j+1/2, k+1/2). Se observă că toate valorile componentelor câmpului din membrul drept al ecuaţiei sunt evaluate la momentul n, inclusiv valoarea câmpului electric Ex

ce apare datorită conductivităţii σ a materialului. Întrucât valorile lui Ex la momentul n nu sunt stocate în memorie, ci doar valorile lui Ex la momentul n-1/2, este necesară estimarea acestor valori.

O metodă de estimare a lui Ex la momentul n este aşa numita aproximaţie semi-implicită (medie aritmetică):

1/ 2 1/ 2, 1/ 2, 1/ 2 , 1/ 2, 1/ 2

, 1/ 2, 1/ 2 2

n n

x i j k x i j kn

x i j k

E EE

+ −+ + + +

+ +

+= (15)

Obs.: se presupun stocate în memorie valorile lui Ex la momentul n-1/2, nu şi la momentul n+1/2. Substituind (15) în (14) şi rearanjând termenii, rezultă:

, 1/ 2, 1/ 2

, 1/ 2, 1/ 21/ 2 1/ 2, 1/ 2, 1/ 2 , 1/ 2, 1/ 2

, 1/ 2, 1/ 2

, 1/ 2, 1/ 2

, 1, 1/ 2, 1/ 2, 1/ 2

, 1/ 2, 1/ 2

, 1/ 2, 1/ 2

12

12

(

12

i j k

i j kn n

x i j k x i j ki j k

i j k

n

z i j k zi j k

i j k

i j k

t

E Et

t

H H

t

σ

ε

σ

ε

ε

σ

ε

+ +

+ ++ −+ + + +

+ +

+ +

+ ++ +

+ +

+ +

∆ −

= +∆

+

∆ − ⋅

∆ +

, , 1/ 2 , 1/ 2, 1 , 1/ 2,

, 1/ 2, 1/ 2 )x

n n n

i j k y i j k y i j k

n

source i j k

H H

y z

J

+ + + +

+ +

−−

∆ ∆

−

(16a)

În mod similar se pot deduce ecuaţiile cu diferenţe finite pentru determinarea lui

yE :

1/2, 1, 1/ 2

1/ 2, 1, 1/21/2 1/21/2, 1, 1/2 1/2, 1, 1/2

1/2, 1, 1/2

1/2, 1, 1/2

1/2, 1, 1/2

1/2, 1, 1/2

1/2, 1, 1/ 2

12

12

12

i j k

i j kn n

y i j k y i j ki j k

i j k

i j k

i j k

i j k

t

E Et

t

t

σ

ε

σ

ε

ε

σ

ε

− + +

− + ++ −− + + − + +

− + +

− + +

− + +

− + +

− + +

∆ −

= +∆

+

∆ ⋅

∆ +

1/2, 1, 1 1/2, 1, , 1, 1/ 2 1, 1, 1/2

1/ 2, 1, 1/2

(

)y

n n n n

x i j k x i j k z i j k z i j k

n

source i j k

H H H H

z x

J

− + + − + + + − + +

− + +

− −−

∆ ∆

−

(16b)

si pentru :z

E

1/2, 1/2, 1

1/2, 1/2, 11/2 1/21/2, 1/2, 1 1/2, 1/2, 1

1/ 2, 1/2, 1

1/2, 1/2, 1

1/2, 1/ 2, 1

1/2, 1/ 2, 1

1/2, 1/2, 1

12

12

12

i j k

i j kn n

z i j k y i j ki j k

i j k

i j k

i j k

i j k

t

E Et

t

t

σ

ε

σ

ε

ε

σ

ε

− + +

− + ++ −− + + − + +

− + +

− + +

− + +

− + +

− + +

∆ −

=∆

+

∆ + ⋅

∆ +

, 1/2, 1 1, 1/ 2, 1 1/ 2, 1, 1 1/2, , 1

1/ 2, 1/2, 1

(

)y

n n n n

y i j k y i j k x i j k x i j k

n

source i j k

H H H H

x y

J

+ + − + + − + + − +

− + +

− −−

∆ ∆

−

(16c)

În mod analog, se pot deduce ecuaţiile cu diferenţe finite pentru

xH , yH si z

H . De exemplu, pentru componenta Hx situată în colţul din dreapta sus al celulei Yee cubice prezentate în fig.1, avem:

1/2, 1, 1

1/2, 1, 111/2, 1, 1 1/2, 1, 1

1/2, 1, 1

1/2, 1, 1

1/2, 1, 3/1/2, 1, 1

1/2, 1, 1

1/2, 1, 1

12

12

(

12

i j k

i j kn n

x i j k x i j k

i j k

i j k

y i j ki j k

i j k

i j k

t

H Ht

t

E

t

σ

µ

σ

µ

µ

σ

µ

∗− + +

− + ++− + + − + +∗

− + +

− + +

− + +− + +

∗− + +

− + +

∆−

= ∆ +

∆ + ⋅ ∆ +

1/2 1/2 1/2 1/22 1/2, 1, 1/2 1/2, 3/2, 1 1/2, 1/2, 1

1/21/2, 1, 1)x

n n n n

y i j k z i j k z i j k

n

source i j k

E E E

z y

M

+ + + +− + + − + + − + +

+− + +

− −−

∆ ∆

−

(17a)

Pentru componenta Hy situată în partea superioară-frontală a celulei Yee, rezultă:

, 1/2, 1

, 1/2, 11, 1/2, 1 , 1/2, 1

, 1/2, 1

, 1/2, 1

1/21/2, 1/2, 1 1/2, 1, 1/2, 1

, 1/2, 1

, 1/2, 1

12

12

(

12

i j k

i j kn n

y i j k y i j k

i j k

i j k

n

z i j k z i ji j k

i j k

i j k

t

H Ht

t

E E

t

σ

µ

σ

µ

µ

σ

µ

∗+ +

+ +++ + + +∗

+ +

+ +

++ + + − ++ +

∗+ +

+ +

∆−

= ∆ +

∆

− + ⋅ ∆ +

1/2 1/2 1/2/2, 1 , 1/2, 3/2 , 1/2, 1/2

1/2, 1/2, 1)y

n n n

k x i j k x i j k

n

source i j k

E E

x z

M

+ + ++ + + + +

++ +

−−

∆ ∆

−

(17b)

Pentru zH :

, 1, 1/2

, 1, 1/21, 1, 1/2 , 1, 1/2

, 1, 1/2

, 1, 1/2

1/2, 3/2, 1/2 , 1/2, 1, 1, 1/2

, 1, 1/2

, 1, 1/2

12

12

(

12

i j k

i j kn n

z i j k z i j k

i j k

i j k

n

x i j k x i j ki j k

i j k

i j k

t

H Ht

t

E E

t

σ

µ

σ

µ

µ

σ

µ

∗+ +

+ +++ + + +∗

+ +

+ +

++ + + ++ +

∗+ +

+ +

∆−

= ∆ +

∆

− + ⋅ ∆ +

1/2 1/2 1/2/2 1/2, 1, 1/2 1/2, 1, 1/2

1/2, 1, 1/2)

z

n n n

y i j k y i j k

n

source i j k

E E

y x

M

+ + ++ + + − + +

++ +

−−

∆ ∆

−

(17c)

Din ecuaţiile (16) şi (17) se observă (fig.3) că valoarea câmpului într-un punct din spaţiu,

la un moment de timp, depinde doar de valoarea sa la momentul precedent şi de valorile adiacente ale câmpurilor.

Fig. 3 Valoarea câmpului într-un punct, la un moment particular. Pentru implementarea FDTD într-o regiune ai cărei parametri dielectrici variază în mod

continuu cu poziţia spaţială, se definesc şi stochează următorii coeficienţi de actualizare ai câmpului electromagnetic, la momentul iniţial (t=0)

Coeficienţii de actualizare pentru câmpul electric din punctul (i,j,k):

, , , ,, ,

, , , ,

1 12 2

i j k i j k

a i j k

i j k i j k

t tC

σ σ

ε ε

∆ ∆= − +

(18a)

, ,1 , ,

, , 1 , ,

12

i j k

b i j k

i j k i j k

ttC

σ

ε ε

∆∆= + ∆

(18b)

H

E E E E

E E E E

E E E E

H H

H H H

0x = 2

xx

∆= x x= ∆

3

2

xx

∆= 2x x= ∆

5

2

xx

∆= 3x x= ∆

0t =

2

tt

∆=

t t= ∆

3

2

tt

∆=

2t t= ∆

, ,2 , ,

, , 2 , ,

12

i j k

b i j k

i j k i j k

ttC

σ

ε ε

∆∆= + ∆

(18c)

Coeficienţii de actualizare pentru câmpul magnetic din punctul (i,j,k):

, , , ,, ,

, , , ,

1 12 2

i j k i j k

a i j k

i j k i j k

t tD

σ σ

µ µ

∗ ∗ ∆ ∆= − +

(19a)

, ,1 , ,

, , 1 , ,

12

i j k

b i j k

i j k i j k

ttD

σ

µ µ

∗ ∆∆= + ∆

(19b)

, ,2 , ,

, , 2 , ,

12

i j k

b i j k

i j k i j k

ttD

σ

µ µ

∗ ∆∆= + ∆

(19c)

unde: ∆1 şi ∆2 reprezintă cei doi incremenţi spaţiali care se modifică atunci când se calculează câmpul într-un punct. De exemplu: în ecuaţia (16a), se modifică doar indicii după y şi după z, rezultă că ∆1=∆y şi ∆2=∆z. În ecuaţia (16b) se modifică doar indicii după z şi x. Rezultă că: ∆1=∆z şi ∆2=∆x, ş.a.m.d. Pentru o reţea cubică: ∆x=∆y=∆z=∆ şi, deci: ∆1=∆2=∆, iar

1 2b b bC C C= = şi 1 2b b b

D D D= = .

1.4.2. Mărimea celulei şi criteriul de stabilitate

Mărimea unei celule 3D în algoritmul FDTD este un parametru critic în simulare. Ea

trebuie să fie suficient de mică încât să permită obţinerea soluţiei cu precizie bună în cazul frecvenţei maxime ce prezintă inters. Pe de altă parte, discretizarea domeniului analizat în paşi foarte fini, ar duce la cerinţe de memorie şi timp de calcul excesiv de mari. Întrucât lungimea de undă se modifică în funcţie de parametrii electrici ai mediului prin care se propagă semnalul, dimensiunea celulei va depinde şi ea de aceşti parametri.Astfel, odată cu creşterea permitivităţii şi conductivităţii, lungimea de undă corespunzătoare unei anumite frecvenţe scade, prin urmare trebuie micşorată şi dimensiunea celulei.

Rezoluţia spaţială sau pasul de eşantionare se alege o fracţiune din lungimea de undă minimă ce se propagă în reţea. Se observă din figura 4 că, în cazul unei eşantionări mai grosiere (mai puţin de 2 eşantione pe perioadă spaţială), în reţea se pot propaga moduri de undă cu viteză supraluminală, însă acestea sunt rapid atenuate. Dacă se prelevează mai puţin de 3 eşantioane pe o perioadă spaţială (λmin) atenuarea undelor numerice ce se propagă în reţea creşte exponenţial [11].

1 2 3 4 5 6 7 8 90

1

2

3

4

5

6

Variatia vitezei de faza normate, vp/c, in functie de rezolutia de esantionare spatiala, lambda/dx

viteza d

e f

aza n

orm

ata

si ate

nuare

a in r

ete

a [

Np/c

elu

la]

rezolutia spatiala

viteza de fazaatenuarea in retea

Fig. 4 Dependenţa vitezei de fază şi a atenuării de perioada de eşantionare spaţiala. În concluzie, rezoluţia spaţială (distanţa între 2 puncte succesive) trebuie aleasă astfel

încât, de la un punct la altul, câmpul electromagnetic să nu se modifice semnificativ. Se alege iniţial x y z∆ = ∆ = ∆ . Pentru asigurarea stabilităţii numerice a algoritmului este necesară satisfacerea unei relaţii [12] între incrementul spaţial şi cel temporal, t∆ . Dacă ε şi µ sunt constante, stabilitatea numerică se obţine dacă:

( ) ( ) ( )

max

2 2 2

1

1 1 1v t

x y z

∆ ≤

+ +∆ ∆ ∆

(20)

unde maxv - viteza maximă de propagare în reţea.

Dacă ε şi µ variază, este mai dificil de obţinut un criteriu de stabilitate.

Condiţia (20) impune restricţii asupra incrementului temporal, pentru incremenţi spaţiali impuşi.

În cazul propagării doar pe directia x, se notează cu c t

Sx

∆=

∆ceea ce vom numi în

continuare factor de stabilitate (numit şi număr Courant )[11]. În mod similar se definesc şi factorii de stabilitate pentru celelalte direcţii. Se demonstrează că dacă 0 1S< ≤ , unda care se propagă în reţea nu se atenuează în timp

(în condiţiile în care mediul de propagare nu are pierderi), iar dacă numărul Courant devine supraunitar, amplitudinea undei creşte exponenţial cu fiecare increment temporal [11].

În principiu, incremenţii spaţiali trebuie astfel aleşi încât să redea cu acurateţe detaliile geometrice ale structurii analizate şi pentru a distinge între componentele spectrale ale lungimilor de undă ce se propagă în reţea. Odată ce s-au fixat incremenţii spaţiali, incrementul temporal este ales astfel încat să se asigure condiţia de stabilitate numerică.

Cazuri extreme:

• Dacă eşantionarea spaţială şi temporală este foarte fină ( 0, 0x t∆ → ∆ → ), soluţia numerică devine tot mai precisă, dar creşte considerabil volumul şi timpul de calcul. Se poate arăta că atât viteza de fază cât şi viteza de grup sunt egale cu c (viteza luminii în vid), indiferent de frecvenţă (propagare nedispersivă), prin urmare algoritmul numeric este nedispersiv [11]

• Eşantionarea cu “pasul magic” (magic time-step): c t x∆ = ∆ .În acest caz soluţia obtinută este exactă, indiferent de pasul de eşantionare în timp şi spaţiu.

• Propagare dispersivă. Aceasta este situaţia cea mai des întâlnită şi are loc dacă 1S < şi x∆ e comparabil cu minλ . S-a observat că unda “numerică” se propagă în reţea cu o

viteză mai mică decât unda “analogică” echivalentă. Pe masură ce x∆ devine mult mai mic decât

minλ , dispersia numerică devine tot mai mică.

Cap 2. Analiza propagării câmpului într-un spaţiu unidimensional

Se va analiza în continuare propagarea unui câmp electromagnetic într-un spaţiu

unidimensional. Presupunem că nu există variaţii după y nici ale câmpului, nici ale geometriei reţelei. Prin urmare, derivatele parţiale în raport cu y şi z sunt nule. Propagarea se face în acest caz de-a lungul axei x. De asemenea, presupunem că nu există surse ale câmpului ( 0

xsourceJ = ,

0xsourceM = )

Ecuaţiile lui Maxwell (ecuaţiile 1 şi 2) devin:

yzHE

t xε

∂∂=

∂ ∂, y z

H E

t xµ

∂ ∂=

∂ ∂, (21)

(mod TMz) respectiv,

yzEH

t xµ

∂∂− =

∂ ∂, yz

EH

x tε

∂∂− =

∂ ∂ (22)

(mod TEz)

Ecuaţiile (23) se rescriu astfel, utilizând diferenţe finite:

( )1/ 2 1/ 2 1

1/ 2 1/ 21/ 2 1/ 2i i i i

n n n n

z z y ya i b iE C E C H H

− − −

+ −

− −= + − (23)

( )1/ 2 1/ 2

1 1/ 2 1/ 2

i i i i

n n n n

y y z za i b iH D H D E E

+ −

+ + += + − (24)

În continuare se prezintă rezultatele simulării (program Matlab în Anexa 2) unor cazuri

de propagare într-un mediu unidimensional pentru diverse forme de undă, propagarea făcându-se după direcţia x.

Notaţii: Spaţiul discretizat este format din puncte de calcul, notate cu i. Între două puncte succesive distanţa este x∆ . Câmpul este evaluat la momente de timp notate cu t. Distanţa temporală între două momente succesive este t∆ . În continuare, reprezentările grafice sunt făcute în funcţie de i şi t, nu în funcţie de x∆ şi t∆ .

În figurile 5 şi 6 s-a reprezentat propagarea unui impuls gaussian într-un spaţiu 1D,

impunând în toate punctele spaţiului un factor Courant: S=1.0008. S-a reprezentat impulsul în reţea la momentul t=200, respectiv la momentul t=250. Se observă că, pe măsura propagării în

0 50 100 150 200 250 300 350 400 450 500-0.2

0

0.2

0.4

0.6

0.8

1

1.2

S=1.0008

54 56 58 60 62-0.3

-0.25

-0.2

-0.15

-0.1

-0.05

0

0.05

zoom in

Propagarea unui impuls gaussian intr-un spatiu 1D

space points

E[V

/m]

reţea, impulsul este urmat de oscilaţii a căror perioadă spaţială este 2 xλ = ∆% [11], aşa cum se poate vedea în detaliul evidenţiat. Aceasta oscilaţie, care se suprapune peste unda utilă, are o creştere exponenţială cu timpul de analiză [11] şi stă la originea instabilităţii numerice. Din condiţia (20) rezultă că, pentru a nu apărea instabilitatea numerică la propagarea unei unde într-un spaţiu unidimensional, este necesar ca factorul Courant să fie subunitar.

Fig 5. Apariţia instabilităţii numerice pentru factor Courant: S=1.0008

0 50 100 150 200 250 300 350 400 450 500-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

1.2

Propagarea unui impuls gaussian intr-un spatiu 1D

space points

Ez [

V/m

]

impulsul la t=200

impulsul la t=250

impulsul la t=280

Fig.6 Impuls gaussian în reţea, la 3 momente de timp distincte

Din figura 6 se observă că, pe măsura propagării în spaţiu, oscilaţiile ce urmează impulsului au amplitudine tot mai mare. În figura 7, criteriul de stabilitate este încălcat într-un singur punct al spaţiului (s-a impus

S=1.0008 la i=90), punct în care şi apar oscilaţiile cu frecvenţa 1

2 x∆. În timp, amplitudinea

acestor oscilaţii va creşte exponenţial, înecand practic semnalul util şi ducând la depăşiri ale limitelor de calcul ale calculatorului. (“run-time floating-point overflow”).

0 50 100 150 200 250 300 350 400 450 500-5

-4

-3

-2

-1

0

1

2

3

4Propagarea unui impuls gaussian in spatiul 1-D

E [

V/m

]

Space points

impulsul la t=80 impulsul la t=150 impulsul la t=240

Fig. 7 Propagarea impulsului gaussian în cazul în care S>1 doar într-un singur punct.

Se observă (fig.8) că, dacă factorul de stabilitate este subunitar, unda suferă o uşoară dispersie (tot mai accentuată pe măsură ce S scade). Pentru a putea face compararea, s-a urmărit situaţia impulsului în reţea la acelaşi moment de timp, considerând însă 2 cazuri particulare: S=1 şi S=0.5. Pentru S=1, atât forma cât şi întinderea spaţială a impulsului se păstrează. În acest caz discontinuităţile impulsului sunt perfect modelate. Dacă S scade (de exemplu S=0.5) discontinuităţile impulsului generează oscilaţii semnificative [11]. S-a observat că durata şi perioada acestor oscilaţii cresc pe masură ce S scade.

0 20 40 60 80 100 120 140 160 180 200-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

1.2

1.4

Propagarea unui impuls dreptunghiular intr-un spatiu 1-D

E[V

/m]

Space points

S=1S=0.5

Fig. 8 Propagarea unui impuls dreptunghiular în reţea unidimensională, când S<1, respectiv S=1 Dacă însă în reţea se transmite un impuls gaussian, după cum se observă din fig 9, acesta va fi afectat nesemnificativ de fenomenul de dispersie, atunci când S<1.

0 20 40 60 80 100 120 140 160 180 200-0.2

0

0.2

0.4

0.6

0.8

1

1.2Propagarea unui impuls gaussian intr-un spatiu 1-D

E[V

/m]

Space points

S=1S=0.5

Fig. 9. Propagarea unui impuls gaussian în reţea unidimensională pentru 1S ≤ În cazul impulsului dreptunghiular, având armonici de frecvenţă foarte înaltă (adică lungimi de undă tot mai mici ce se propagă în reţea), procesul de eşantionare spaţială nu mai poate distinge între lungimile de undă ale armonicilor de ordin superior. În cazul în care unda trece dintr-un mediu material în altul, o parte a undei va fi transmisă în noul mediu, iar o parte a undei se reflectă. De asemenea, viteza de propagare în noul mediu se modifică. Am verificat acest lucru, utilizând un semnal triunghiular, cu durata 188T t= ∆ . S-a determinat evoluţia în timp a câmpului într-un spaţiu discretizat (980 de puncte), cuprinzând două medii: aer ( 1ε = ) în primele 327 de puncte ale spaţiului şi un mediu oarecare (caracterizat de 9ε = ) în restul punctelor din spaţiu. În figurile 10-12 se prezintă evoluţia în timp a semnalului pe măsura propagării în spaţiul numeric. Întrucât al doilea mediu este caracterizat de

9ε = , la trecerea în acest mediu, viteza de propagare a undei trebuie sa fie de 3 ori mai mică decât cea corespunzătoare propagării prin aer. Pe de altă parte, datorită discontinuităţii de impedanţă, doar jumătate din energia semnalului va fi transmisă în al doilea mediu, cealaltă jumătate fiind reflectată (unda reflectată suferă un defazaj de 180°)

În figura 11 se prezintă variaţia în spaţiu a impulsului la momentul particular 515t = . Impulsul generat la 0t = ajunge după 327 t∆ în punctul i=327, adică în punctul de discontinuitate a impedanţei. La 515 t∆ , impulsul reflectat a ajuns în punctul i=327-(515-327) ≅ 139, iar cel transmis, întrucât se propagă cu viteza de 3 ori mai mică, a ajuns doar în punctul i=327+(515-327)/3 ≅ 390. Se observă, de asemenea, că cele două unde, reflectată şi transmisă, au amplitudinile de două ori mai mici decât ale semnalului transmis în reţea, generat la i=1 (fig.10)

0 500 1000 1500 2000 2500 30000

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1Semnalul in timp, pentru i=1

time samples

Ez [

V/m

]

Fig. 10 Generarea unui semnal triunghiular, la i=1.

0 100 200 300 400 500 600 700 800 900 1000-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

space samples

Ez [

V/m

]

Semnalul in retea la momentul t = 517

Fig. 11 Variaţia în spaţiu a semnalului, la momentul de timp t=515.

În figura 12 se reprezintă atât câmpul electric, cât şi cel magnetic, la momentul 900t = , presupunând că în tot spaţiul permitivitatea electrică este 9ε = , iar durata impulsului ce se propagă în reţea este de 936 t∆ . Datorită permitivităţii supraunitare, la momentul 900t = , impulsul va fi ajuns doar în punctul 300i = . Se ştie că raportul între câmpul electric şi cel magnetic într-un punct din spaţiu este egal cu impedanţa mediului în acel punct. Întrucât 9ε = iar

1µ = , impedanţa mediului este 1/3. Se observă în figura 12 că raportul între intensitatea câmpului electric şi-a celui magnetic este 1/3.

0 100 200 300 400 500 600 700 800 900 1000-3

-2.5

-2

-1.5

-1

-0.5

0

0.5

1

space points

Ez [

V/m

], H

y[A

/m]

Campul electric/magnetic in retea la momentul t = 900

câmpul electric, Ezcâmpul magnetic, Hy

Fig 12. Intensitatea câmpului electric, respectiv magnetic într-un mediu cu 9ε =

În simulările anterioare, s-a considerat conductivitatea mediului nulă. În cazul în care conductivitatea este nenulă, se observă (figura 13) că unda este atenuată pe măsura propagării în reţea şi este afectată de dispersie.

0 100 200 300 400 500 600 700 800 900 10000

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

space samples

Ez [

V/m

], H

y[A

/m]

Campul electric intr-un mediu cu conductivitatea unitara

t=200t=900t=1900

Fig. 13 Propagarea unui semnal triunghiular cu durata 900T t= ∆ într-un mediu cu conductivitatea unitară.

Cap 3. Analiza propagării câmpului într-un spaţiu bidimensional

Se analizează în continuare propagarea unui câmp într-un spaţiu cu două dimensiuni. Se presupune că, în orice punct din spaţiu, componentele câmpului nu depind de coordonata z a punctului. Dacă unda incidentă este uniformă în direcţia z, toate derivatele parţiale ale câmpului în raport cu z sunt nule. Pentru simplificare, ε şi µ se consideră constante, iar J=0 şi M=0 (nu există surse de câmp). Astfel, singura sursă este unda incidentă.

În cazul în care anumite componente de câmp sunt nule, rezultă două moduri de propagare, astfel:

• Dacă 0z

E = şi 0x yH H= = , modul de propagare se numeşte transversal electric.

• Dacă 0z

H = şi 0x yE E= = , modul de propagare se numeşte transversal magnetic.

Se prezintă în continuare ecuaţiile celor două moduri de propagare:

1) Propagare transversal electrică (TEz)

0x yH H= = , 0z

E = ,

( )1x

x zsource x

E HJ E

t yσ

ε

∂ ∂= − + ∂ ∂

(25a)

( )1y

y zsource y

E HJ E

t xσ

ε

∂ ∂ = − − + ∂ ∂

(25b)

( )1z

yxzsource z

EEHM H

t y xσ

µ∗

∂ ∂∂= − − +

∂ ∂ ∂ (25c)

2) Propagare transversal magnetică (TMz)

0x yE E= = , 0z

H = ,

( )1x

x zsource x

H EM H

t yσ

µ∗ ∂ ∂

= − − + ∂ ∂ (26a)

( )1y

y zsource y

H EM H

t xσ

µ∗

∂ ∂ = − − + ∂ ∂

(26b)

( )1z

y xzsource z

H HEJ E

t x yσ

ε

∂ ∂∂= − − +

∂ ∂ ∂ (26c)

Se poate observa că modurile TEz şi TMz nu conţin componente comune de câmp şi, prin

urmare, pot exista simultan fără a interacţiona.

Se consideră în continuare modul TMz. Ecuaţiile (26) se rescriu folosind diferenţele finite, în care considerăm că nu există surse de câmp (

sourceJ =0 şi

sourceM =0)

1/2 1/21/ 2, 1/2 1 1/2, 1/2 1/2, 1/2

, 1/2 1, 1/2 1/2, 1 1/2,

1/2, 1/2 ( )

n n

z i j a i j y i j

n n n n

y i j y i j x i j x i j

b i j

E C E

H H H HC

x y

+ −− + − + − +

+ − + − + −

− +

=

− −+ ⋅ −

∆ ∆

(27a)

1/2 1/21/2, 3/2 1/2, 1/21

1/2, 1 1/2, 1 1/2, 1 1/2, 1

n n

z i j z i jn n

x i j a i j x i j b i j

E EH D H D

y

+ +− + − ++

− + − + − + − +

−= − ⋅

∆ (27b)

1/2 1/21/2, 1/2, 1 1/2, 1/2, 11

, 1/2 , 1/2 , 1/2, , 1/2

n n

z i j k z i j kn n

y i j a i j y i j b i j

E EH D H D

x

+ +

+ + + − + ++

+ + + +

−= + ⋅

∆ (27c)

În conditiile în care rezoluţiile de eşantionare spaţială pe cele două direcţii sunt egale, criteriul de stabilitate (20) impune ca perioada maximă de eşantionare temporală să fie:

max2

tc

∆∆ =

unde x y∆ = ∆ = ∆ , reprezintă rezoluţia de eşantionare spaţială. Aceasta înseamnă că, între două momente succesive de timp (este vorba de timpul discretizat din

FDTD), unda parcurge o distanţă maximă de 2

∆

În continuare este ilustrată propagarea unui impuls gaussian(fig. 14) în spaţiul

bidimensional. S-a considerat: x y∆ = ∆ = ∆ şi 2

tc

∆∆ = . În vederea ilustrării propagării

câmpului în spaţiul 2D, s-a realizat programul Matlab prezentat în anexa 2 (fdtd_2d.m) Impulsul este generat în centrul reţelei

Fig. 14 Impuls gaussian

În figurile 15 (vizualizare 3D) şi 16 (vizualizare în plan) este redat câmpul Ez la momentul 70t = . Se observă că impulsul generat în centrul reţelei se propaga radial.

Fig. 15 Impulsul în reţea la 70t =

Fig. 16 Impulsul în reţea la 70t = , vedere în plan.

Pe măsură ce timpul de propagare creşte, amplitudinea impulsului scade, invers proporţional cu distanţa parcursă. În figura 17 este redat impulsul în reţea la momentul 130t = . Am considerat că viteza de propagare este cea din spaţiul liber,

0 0

1c

ε µ= . În simularile Matlab, pentru simplificare, am considerat viteza luminii normată,

adică 1c = . Impulsul s-a deplasat, faţă de momentul 70t = , cu 60 / 2 42≅ celule. Acest lucru se observă în figura 18.

Fig. 17 Impulsul în reţea la 130t =

Fig. 18 Impulsul în reţea la 130t = , vedere in plan.

În figura 19 este redată variaţia în spaţiu a impulsului, de-a lungul unei direcţii date (am impus ct.y = ) la momentul 130t = .

Fig. 19 Variaţia spaţială a impulsului, pe o direcţie impusă.

Odată cu creşterea timpului de analiză, semnalul ce se propagă în spaţiul bidimensional va ajunge la un moment dat la capetele domeniului analizat. În algoritmul propus de Yee, spaţiul în care se propagă semnalul este terminat cu pereţi conductori (σ = ∞ ). Prin urmare, când unda ajunge la marginile domeniului analizat, va fi reflectată în întregime. În cazul în care se doreşte analiza unui domeniu cu întindere infinită (exemplu: propagarea în aer) şi oricând se impune trunchierea spaţiului analizat, trebuie modificate condiţiile la frontieră, astfel încât să nu apară reflexii. S-au dezvoltat diverse metode de simulare a absorbţiei undelor, denumite generic ABCs (Absorbing Boundary Conditions), dintre care cea mai cunoscută a fost introdusă de J. P. Berenger, şi anume, metoda PML (Perfectly Matched Layer). În cazul condiţiilor PML, între domeniul analizat şi frontieră (care rămâne perfect conductoare) se inserează un strat absorbant, cu conductivitatea variabilă (sigma creşte, liniar sau după altă lege, către marginile domeniului). Astfel, dacă stratul absorbant are grosime suficient de mare, undele care ajung la capătul domeniului vor fi puternic atenuate, iar reflexiile neglijabile. Pentru implementarea metodei PML, se introduc, în plus, în ecuaţiile ce caracterizează mediul absorbant, conductivitatea electrică şi conductivitatea magnetică (mărime fictivă), cu scopul de a realiza adaptarea de impedanţe între medii. Se rescriu ecuaţiile de rotor pentru mediul absorbant:

2t

ε σ∂

+ ∇ ×∂

EE = H (29a)

*2 -

tµ σ

∂+ ∇ ×

∂

HH = E (29b)

unde *

2 2, , ,ε µ σ σ sunt parametrii dielectrici ai mediului absorbant

Impedanţa de undă pentru mediul absorbant, are expresia:

*2

2PML

jZ

j

µ σ ω

ε σ ω

+=

+ (30)

La incidenţa normală a unei unde pe acest mediu coeficientul de reflexie a câmpului electric este:

1

1z

PMLE

PML

Z Z

Z Z

−Γ =

+ (31)

unde 11

1

Zµ

ε= este impedanţa mediului din care vine unda.

Pentru a nu apărea reflexii, este suficient ca:

*

0 0

σ σ

µ ε= şi 2 1µ µ= , respectiv 2 1ε ε= (32)

rezultă că, dacă sunt îndeplinite condiţiile (32), 1 PML

Z Z= şi nu apar reflexii.

În cazul incidenţei oblice (cum se întamplă în mediile bi- şi tri-dimensionale), reflexiile sunt dificil de evitat. Pentru aceasta, cazul incidenţei oblice se reduce la cel al incidenţei normale prin descompunerea undei în două component ortogonale. Pentru un spaţiu bidimensional, se consideră în continuare modul de propagare TMz. În mediul absorbant, componenta de câmp electric, Ez, este împărţită în două subcomponente [6]:

z zx zyE E E= + (33)

Ecuaţiile (26) se rescriu pentru mediul absorbant:

*2

x zy x

H EH

t yµ σ

∂ ∂+ = −

∂ ∂ (34a)

*2

y zx y

H EH

t xµ σ

∂ ∂+ = −

∂ ∂ (34b)

2yzx

x zx

HEE

t xε σ

∂∂+ =

∂ ∂ (34c)

2zy x

y zy

E HE

t yε σ

∂ ∂+ = −

∂ ∂ (34d)

unde: 2µ , 2ε , σ şi *σ caracterizează mediul absorbant.

Mediul absorbant propus de Berenger este ilustrat în figura 20, pentru o reţea bidimensională. În stânga şi dreapta reţelei sunt adăugaţi pereţi absorbanţi (PML), a căror conductivitate

electrică/magnetică este x

σ / *xσ şi îndeplineşte condiţia de adaptare:

*

2 2

x xσ σ

µ ε= . Reţeaua este

delimitată şi superior, respectiv inferior, de către pereţi absorbanţi ai căror conductivităţi satisfac,

pentru adaptare, relaţia: *

2 2

y yσ σ

µ ε= . În colţurile reţelei, unde se suprapun două PML, sunt prezente

toate cele 4 tipuri de pierderi: * *, , ,x x y yσ σ σ σ , egale cu pierderile din pereţii adiacenţi.

Fig. 20. Reţea bidimensională cu PML

În cazul în care nu se folosesc pereţi absorbanţi, atunci când unda ajunge la marginea domeniului analizat, se va reflecta total. Domeniul considerat este o matrice 250 x 250, adică întinderea sa pe direcţia x este 250 x⋅∆ , iar pe direcţia y este 250 y⋅ ∆ ( x∆ şi y∆ fiind distanţele între 2 eşantioane spaţiale succesive). Sursa câmpului este plasata excentric, concentrată într-un punct şi are o variaţie sinusoidală. În continuare este ilustrată propagarea undei în situaţia în care se folosesc condiţii PML (fig 21), respectiv fară condiţii PML (fig. 22).

Fig. 21 Câmpul electric la 290t = , fără perete absorbant.

În cazul utilizării PML, pereţii absorbanţi au grosimea de 20 x⋅∆ , respectiv 20 y⋅ ∆ . În acest mediu, datorită conductivităţii nenule, unda este atenuată, astfel că în momentul în care ajunge la marginea domeniului (perete reflector) amplitudinea ei este neglijabilă şi, mai mult, reflexiile sunt, la rândul lor, atenuate, astfel că pot fi practic neglijate. Am considerat că, în peretele absorbant, conductivitatea variază proporţional cu grosimea peretelui, între 0.001 S şi 0.9 S. Se observa că, la acelaşi moment de timp, în cazul folosirii pereţilor absorbanţi, unda nu mai suferă reflexii.

Fig. 22 Câmpul electric la 290t = , cu perete absorbant.

În cazul în care unda care se propagă întâlneşte un obstacol, se observă şi un fenomen de difracţie. În figurile 23 şi 24 am considerat că unda care se propagă într-un mediu cu prmitivitatea unitară (aer), întalneşte un alt mediu (o suprafaţă circulară) a cărei permitivitate relativă este: 4ε = . Sursa câmpului electric este concentrată în centrul reţelei şi emite un impuls gaussian.

Fig. 23 Propagarea unui impuls gaussian în medii cu permitivităţi diferite. Vedere în plan la momentul t=120.

Fig. 24 Propagarea unui impuls gaussian în medii cu permitivităţi diferite. Vedere în plan la momentul t=190.

Cap. 4. Modelarea 3D a mediilor reale.

Simulatorul CST Microwave Studio

CST Microwave Studio –prescurtat CST MWS- este un pachet software dedicat analizei de câmp electromagnetic şi proiectării în domeniul frecvenţelor înalte. Modelarea 3D este facilitată de existenţa unui nucleu ACIS (software dedicat pentru modelarea 3D). După modelarea a obiectului de analizat, se generează, cu ajutorul unui sistem expert, o reţea tridimensională în care este plasat obiectul. CST MWS utilizează patru tipuri de analiză, adaptate diferitelor aplicaţii (analiză în domeniul timp, în domeniul frecvenţă, a modurilor de propagare în ghiduri de undă, a valorilor proprii).

Caracteristici ale CST MWS

• Se permite importul-exportul datelor 3D din/în alte medii de proiectare (ex. AutoCAD, STEP etc) şi parametrizarea datelor importate.

• Calcularea, într-o singură simulare, a parametrilor S de bandă largă, a distribuţiei de câmp în funcţie de timp.

• Se pot analiza materiale izotrope/anizotrope, precum şi materiale ale căror proprietăţi variază cu frecvenţa.

• Condiţii de radiaţie/absorbţie la frontieră.

• Calcularea diferiţilor parametri de câmp: câmp electric/magnetic, curenţi superficiali, transfer de putere, densitate de curent, densitate de energie electrică/magnetică, variaţia în timp şi frecvenţă a tensiunii.

• Calcularea parametrilor de câmp apropiat şi depărtat ai antenelor (cîştig, directivitate etc).

• Calcularea distribuţiei SAR

• Posibilitatea de a plasa elemente discrete în structura analizată.

• Reflectometrie în timp.

• Semnale de excitaţie predefinite şi configurabile de către utilizator.

• Modificarea adaptivă a reţelei (a densităţii de celule) [14] Metoda de simulare

CST MWS se bazează pe tehnica integrării finite (FIT). Aceasta presupune o schemă de discretizare valabilă în tot spectrul de frecvenţe şi potrivită atât analizei in domeniul timp, cât şi analizei în domeniul frecvenţă. FIT presupune discretizarea ecuaţiilor lui Maxwell, exprimate în formă integrală.

(35a)

(35b)

(35c)

(35d)

Spre deosebire de FDTD, FIT elimină eroarea de aproximare prin trepte, prin integrarea unei tehnici suplimentare, denumite PBA (Perfect Boundary Approximation). În figura 25 se prezintă două dintre tehnicile de discretizare numerică cele mai folosite şi neajunsurile lor (segmentarea, în cazul FEM şi aproximarea prin trepte, în cazul FDTD), precum si rezultatul implementării FIT+PBA în cadrul CST MWS.

Fig. 25 Ilustrare FEM, FIT+PBA, FDTD

În figura 26 se prezintă distribuţia câmpului electric într-un model de cap uman, produs de un telefon mobil la frecvenţa de 1.8 GHz, realizat în CST MWS.

Fig. 26 Distribuţia câmpului electric într-un model de cap uman.

În continuare, se va utiliza programul CST MWS pentru modelarea şi simularea interacţiunii între câmpul electromagnetic de înalta frecvenţă şi ţesutul biologic, urmând a se compara rezultatele cu cele ce se vor obţine pe cale experimentală.

Anexa 1

Erori de ordinul 2

Fie o funcţie de forma ( , )u u x t= , dependentă de timp, t, şi spaţiu, x, după o lege oarecare. Într-o simulare numerică, soluţia unei ecuaţii nu poate fi decât o aproximare a soluţiei analitice, din cauza discretizării inerente oricărui calcul efectuat cu un sistem numeric. Prin urmare, dacă ( , )u x t este soluţia unei ecuaţii, aceasta va putea fi determinată doar în puncte

discrete, de forma ( , )i n

x t , unde i-indice spaţial, n-indice temporal.

Dacă se face dezvoltare în serie Taylor a funcţiei ( , )u x t , în jurul punctului i

x , pentru un

moment de timp fixat, n

t , se obţine: 2 2 3 3

2 3,, , ,

( ) ( )( ) ...

2 6n i n

i n i n

i t x txi tn x t x t

u x u x uu x x u x

x x x

∂ ∆ ∂ ∆ ∂+ ∆ = + ∆ + + +

∂ ∂ ∂

2 2 3 3

2 3,, , ,

( ) ( )( ) ...

2 6n i n

i n i n

i t x txi tn x t x t

u x u x uu x x u x

x x x

∂ ∆ ∂ ∆ ∂+ ∆ = − ∆ + − +

∂ ∂ ∂

Însumând cele două ecuaţii şi împărţind la 2( )x∆ , se obţine o aproximare pentru 2

2

u

x

∂

∂,

astfel:

2

, 22 2

( ) 2 ( )[( ) ]

( )n i n ni t x t i tu x x u u x xu

O xx x

+ ∆ − + − ∆∂= + ∆

∂ ∆

unde 2[( ) ]O x∆ -eroare de ordinul 2, din cauza neglijării termenilor de ordin superior lui 2, din seria Taylor.

Anexa 2

FDTD_1D.m

%program propagare impuls gaussian într-un spaţiu unidimensional, folosind FDTD. Conditiile la frontiera sunt cele pentru un conductor perfect. clear all close all %definire constante dx=1; %increment spatial dt=1.0008; %increment temporal c=1; %viteza luminii in vid, normata. %S=c*dt/dx - factorul Courant M=500; N=700; i=1:M; %puncte in spatiu pentru Ez n=1:N; %momente de timp p=1:M-1; %initializare campuri ez(n,i)=0; hy(n,p)=0; %definire impuls medt=N/(5.25); sig=N/30; ez(:,1)=exp(-(n(:)-medt).*(n(:)-medt)/(2*sig*sig)); figure(1) plot(n, ez(:,1)) title('impulsul in timp, ptr i=1') xlabel('time samples') ylabel('Ez') grid minor hold on %definire constante eps(1:M)=1; %permitivitatea relativa a aerului miu=1; %permeabilitatea relativa a aerului sig=0; %conductivitate electrica sig_m=0; %conductivitate magnetica d1=dx; d2=dx; %initializarea si calculul coeficientilor de actualizare ai campului Ca(1:M)=1; Cb1(1:M)=1; Cb2(1:M)=1;

Ca(:)=(1-sig*dt./(2*eps(:)))./(1+sig*dt./(2*eps(:))); Cb1(:)=(dt./(eps(:)*d1))./(1+sig*dt./(2*eps(:))); Cb2(:)=(dt./(eps(:)*d2))./(1+sig*dt./(2*eps(:))); Da=(1-sig_m*dt/(2*miu))/(1+sig_m*dt/(2*miu)); Db1=(dt/(miu*d1))/(1+sig_m*dt/(2*miu)); Db2=(dt/(miu*d2))/(1+sig_m*dt/(2*miu)); %calculul campului for tt=2:N; k=2:M-1; p=1:M-1; ez(tt,k)=Ca(k).*ez(tt-1,k)+Cb1(k).*(hy(tt-1,k)-hy(tt-1,k-1)); hy(tt,p)=Da*hy(tt-1,p)+Db1*(ez(tt,p+1)-ez(tt,p)); end figure(2) plot(ez(:,60),'r') xlabel('time samples') ylabel('Ez') title('variatia campului in punctul i = 60') grid on figure(3) plot(i,ez(200,:),'r',i,ez(250,:),'b') xlabel('space samples') ylabel('Ez') title('impulsul in grid la momentul t = 200, respectiv t=250' ) grid on

FDTD_2D.m clear all close all %definire constante Mx=250; My=250; g=20; %grosimea peretelui absorbant Mx1=g; Mx2=Mx-(g-1); My1=g; My2=My-(g-1); N=320; %momente de timp x=1:Mx; y=1:My; xx=1:Mx-1; yy=1:My-1; t=1:N; lam0=10; %cm c=1; %viteza normata; dx=lam0/20; dy=dx; dt=dx/(sqrt(2)*c); d1=dx; d2=dy; %generare impuls ezx(x,y,t)=0; ezy(x,y,t)=0; ez(x,y,t)=0; hx(xx,yy,t)=0; hy(xx,yy,t)=0; medt=30; sig=5; c1=round(2*Mx/3); c2=round(2*My/3); f=1/(47*dt); ezx(c1,c2,:)=sin(2*pi*f*t(:)); ezx(c1,c2,:)=sin(2*pi*f*t(:)); %ezx(c1,c2,:)=exp(-(t(:)-medt).*(t(:)-medt)/(2*sig*sig)); %ezy(c1,c2,:)=exp(-(t(:)-medt).*(t(:)-medt)/(2*sig*sig)); ezx2(c1,c2,:)=ezx(c1,c2,:); ezy2(c1,c2,:)=ezy(c1,c2,:); s=zeros(1,N); s(:)=ezx(c1,c2,:); figure(1) plot(t,s)

xlabel ('time samples') ylabel('\itEz [V/m]') grid on miu=1; epsilon(x,y)=1; sigx(x,y)=0.0; sig_mx(xx,yy)=0; sigy(x,y)=0; sig_my(xx,yy)=0; %generare sigma variabil smax=0.9; smin=0.1; p=(smax-smin)/g; %panta s1=p*(1:g); K=1; s2_sjx=repmat(s1,My1,1); s2_sjy=repmat(s1',1,Mx1); s2_sy=repmat(s1,My2-My1-1,1); s2_sx=repmat(s1',1,Mx2-Mx1-1); %stanga jos sigx(1:Mx1,1:My1)=fliplr(s2_sjx)'; sigy(1:Mx1,1:My1)=fliplr(s2_sjy'); sig_mx(1:Mx1-1,1:My1)=fliplr(s2_sjx(1:g-1,1:g))'/K; sig_my(1:Mx1-1,1:My1)=fliplr(s2_sjx(1:g-1,1:g))/K; %stanga sigx(1:Mx1,My1+1:My2-1)=fliplr(s2_sy)'; sigy(1:Mx1,My1+1:My2-1)=0; sig_mx(1:Mx1-1,My1+1:My2-1)=fliplr(s2_sy)'/K; sig_my(1:Mx1-1,My1+1:My2-1)=0; %stanga sus sigx(1:Mx1,My2:My)=fliplr(s2_sjx)'; sigy(1:Mx1,My2:My)=(s2_sjy'); sig_mx(1:Mx1,My2:My-1)=fliplr(s2_sjx)'/K; sig_my(1:Mx1,My2:My-1)=s2_sjx/K; %dreapta jos sigx(Mx2:Mx,1:My1)=s2_sjx'; sigy(Mx2:Mx,1:My1)=fliplr(s2_sjx); sig_mx(Mx2:Mx-1,1:My1-1)=s2_sjx'/K; sig_my(Mx2:Mx-1,1:My1-1)=fliplr(s2_sjx)/K; %dreapta sigx(Mx2:Mx,My1+1:My2-1)=s2_sy'; sigy(Mx2:Mx,My1+1:My2-1)=0; sig_mx(Mx2:Mx-1,My1:My2-1)=s2_sy'/K; sig_my(Mx2:Mx-1,My1:My2-1)=0; %dreapta sus sigx(Mx2:Mx,My2:My)=s2_sjx'; sigy(Mx2:Mx,My2:My)=s2_sjy'; sig_mx(Mx2:Mx-1,My2:My-1)=s2_sjx'/K;

sig_my(Mx2:Mx-1,My2:My-1)=s2_sjx/K; %jos sigx(Mx1+1:Mx2-1,1:My1)=0; sigy(Mx1+1:Mx2-1,1:My1)=fliplr(s2_sy); sig_mx(Mx1:Mx2-1,1:My1-1)=0; sig_my(Mx1:Mx2-1,1:My1-1)=fliplr(s2_sy)/K; %sus sigx(Mx1+1:Mx2-1,My2:My)=0; sigy(Mx1+1:Mx2-1,My2:My)=s2_sx'; sig_mx(Mx1:Mx2-1,My2:My-1)=0; sig_my(Mx1:Mx2-1,My2:My-1)=s2_sx'/K; Caqx(1:Mx,1:My)=1; Cb1qx(1:Mx,1:My)=1; Cb2qx(1:Mx,1:My)=1; Daqx(1:Mx-1,1:My-1)=1; Db1qx(1:Mx-1,1:My-1)=1; Db2qx(1:Mx-1,1:My-1)=1; Caqy(1:Mx,1:My)=1; Cb1qy(1:Mx,1:My)=1; Cb2qy(1:Mx,1:My)=1; Daqy(1:Mx-1,1:My-1)=1; Db1qy(1:Mx-1,1:My-1)=1; Db2qy(1:Mx-1,1:My-1)=1; Caqx(:,:)=(1-sigx(:,:)*dt./(2*epsilon(:,:)))./(1+sigx(:,:)*dt./(2*epsilon(:,:))); Cb1qx(:,:)=(dt./(epsilon(:,:)*d1))./(1+sigx(:,:)*dt./(2*epsilon(:,:))); Cb2qx(:,:)=(dt./(epsilon(:,:)*d2))./(1+sigx(:,:)*dt./(2*epsilon(:,:))); Daqx=(1-sig_mx(:,:)*dt/(2*miu))./(1+sig_mx(:,:)*dt/(2*miu)); Db1qx=(dt/(miu*d1))./(1+sig_mx(:,:)*dt/(2*miu)); Db2qx=(dt/(miu*d2))./(1+sig_mx(:,:)*dt/(2*miu)); Caqy(:,:)=(1-sigy(:,:)*dt./(2*epsilon(:,:)))./(1+sigy(:,:)*dt./(2*epsilon(:,:))); Cb1qy(:,:)=(dt./(epsilon(:,:)*d1))./(1+sigy(:,:)*dt./(2*epsilon(:,:))); Cb2qy(:,:)=(dt./(epsilon(:,:)*d2))./(1+sigy(:,:)*dt./(2*epsilon(:,:))); Daqy=(1-sig_my(:,:)*dt/(2*miu))./(1+sig_my(:,:)*dt/(2*miu)); Db1qy=(dt/(miu*d1))./(1+sig_my(:,:)*dt/(2*miu)); Db2qy=(dt/(miu*d2))./(1+sig_my(:,:)*dt/(2*miu)); %calcul ptr ez, hx si hy for tt=2:N; kkx=[2:Mx-1]; kky=[2:My-1]; ppx=[1:Mx-1]; ppy=[1:My-1];

ezx(kkx,kky,tt)=Caqx(kkx,kky).*ezx(kkx,kky,tt-1)+Cb1qx(kkx,kky).*(hy(kkx,kky,tt-1)-hy(kkx-1,kky,tt-1)); ezx(c1,c2,tt)=ezx2(c1,c2,tt); ezy(kkx,kky,tt)=Caqy(kkx,kky).*ezy(kkx,kky,tt-1)-Cb2qy(kkx,kky).*(hx(kkx,kky-1,tt-1)-hx(kkx,kky,tt-1)); ezy(c1,c2,tt)=ezy2(c1,c2,tt); ez(kkx,kky,tt)=ezx(kkx,kky,tt)+ezy(kkx,kky,tt); hx(ppx,ppy,tt)=Daqx(ppx,ppy).*hx(ppx,ppy,tt-1)-Db2qx(ppx,ppy).*(ez(ppx,ppy,tt)-ez(ppx,ppy+1,tt)); hy(ppx,ppy,tt)=Daqy(ppx,ppy).*hy(ppx,ppy,tt-1)+Db1qy(ppx,ppy).*(ez(ppx+1,ppy,tt)-ez(ppx,ppy,tt)); end t1=280; figure(2) %hy [x1,y1]=meshgrid(yy,xx); hl1=surf(x1,y1,hy(:,:,t1),'FaceColor','interp',... 'EdgeColor','interp',... 'AlphaDataMapping','none'); alpha(0.1) figure(3) %ez [x1,y1]=meshgrid(y,x); contour(x1,y1,ez(:,:,t1)); %alpha(0.2); figure(4) %vedere-n lungul unei axe a=size(ez(:,round(My/2),60)); b=a(1); s=zeros(N,b); s(t1,:)=ez(:,round(My/2),t1); plot(s(t1,:)) grid minor xlabel('space samples in \itx direction') ylabel('\itEz [V/m]') grid on %reprezentare sigx, sigy -doar ptr test figure(5) mesh(sigy) hold on mesh(sigx) hold off axis equal

Bibliografie [1] K. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media”, IEEE Transactions on Antennas and Propagation, vol. AP-14, nr. 3, pp 302-307. Mai 1966. [2] R. J. Luebbers, F. P. Hunsberger, S. Kunz, „FDTD formulation for frequency dependent permittivity”, Antennas and Propagation Society International Symposium, vol. 1, pp. 50-53, Jun. 1989. [3] C. M. Furse, Jin-Yuan Chen, Om P. Gandhi, „The use of the frequency-dependent Finite-Difference Time-Domain method for induced current and SAR calculations for a heterogeneous model of the human body”, IEEE Transactions On Electromagnetic Compatibility, vol. 36, nr. 2,

pp. 128 – 133, May 1994.

[4] C. Chen, T.M. Babij, „FD-TD analysis of scattering of electromagnetic fields close to phantom humans”, Antennas and Propagation Society International Symposium, vol. 3, pp. 1674 – 1677, Jul. 1996. [5] Paging Systems; European Radio Message System (ERMES) Part 5: Receiver Conformance Specification,” ETS 300 133-5, ETSI, Valbonne, France, July 1992, Amended (Al), Jan. 1994. [6] J. Berenger, “Perfectly matched layer for the FDTD solution of wave-structure interaction problemes”, IEEE Transactions on Antennas and Propagation, vol. 44, nr. 1, pp. 110-117, Jan. 1996 [7] G. Lazzi, Om P. Gandhi, D. Sullivan, „Use of PML absorbing layers for the truncation of the head model in cellular telephone simulations“, IEEE Transactions on Microwave Theory and

Techniques, vol. 48, nr. 11, pp. 2033-2039, Nov. 2000 [8] J. Wang, O. Fujiwara, T. Nojima, „A model for predicting elctromagnetic interference of implanted cardiac pacemakers by mobile telephones“, IEEE Transactions on Microwave Theory

and Techniques, vol. 48, nr. 11, Nov. 2000 [9] T. Togashi, T. Nagaoka, S. Kikuchi s.a., “FDTD Calculations of Specific Absorption Rate in Fetus Caused by Electromagnetic Waves From Mobile Radio Terminal Using Pregnant Woman Model”, IEEE Transactions on Antennas and Propagation, vol. 44, nr. 1, pp. 554-559, Ian. 1996 [10] M. Sadiku, “Numerical techniques in electromagnetics, 2nd Ed.”, CRC Press, 2001. [11] A. Taflove, S. Hagness, “Computational electrodynamics: Finite Difference Time Domain method, 2nd Ed.”, Artech House, 2000.

[12] S. Kunz, R. Luebbers, “Finite Difference Time Domain method for electromagnetics”, CRC Press, 1993. [13] D. Sullivan, “Electronic simulation using FDTD”, IEEE Press, 2000. [14] Y. Liu, C. Sarris, “Fast and efficient modeling of microwave integrated circuit geometries via an adaptive mesh refinement (AMR) - FDTD Technique”, IEEE Transactions on microwave

theory and techniques, Vol. 54, nr. 2, pp. 689-703, Feb. 2006 [15] CST Microwave Studio, “Getting Started” - www.cst.com

CUPRINS

Cap. 1 Rezolvarea ecuaţiilor lui Maxwell prin metoda diferenţelor finite în timp (FDTD) ...................................................................................................................... 1

Introducere .................................................................................................................................. 2 1.1. Scurtă descriere a metodei FDTD ........................................................................................ 3 1.2. Ecuaţiile lui Maxwell în 3 dimensiuni ................................................................................. 4 1.3. Algoritmul lui Yee ............................................................................................................... 6 1.4. Notaţii pentru diferenţe finite ............................................................................................... 8

1.4.1. Ecuaţiile lui Maxwell în trei dimensiuni, scrise cu diferenţe finite .............................. 9 1.4.2. Mărimea celulei şi criteriul de stabilitate .................................................................... 13

Cap 2. Analiza propagării câmpului într-un spaţiu unidimensional ........................ 16 Cap 3. Analiza propagării câmpului într-un spaţiu bidimensional .......................... 24 Cap 4. Modelarea 3D a mediilor reale. Simulatorul CST Microwave Studio ........ 36 Anexa 1 .................................................................................................................... 39 Erori de ordinul 2 ..................................................................................................... 39 Anexa 2 .................................................................................................................... 40 FDTD_1D.m ............................................................................................................ 40 FDTD_2D.m ............................................................................................................ 42 Bibliografie .............................................................................................................. 46