TECHNISCHE UNIVERSITEIT EINDHOVEN FACULTEIT ELEKTROTECHNIEK
tlB
VAKGROEP Elektromagnetisme
COO pakket voor eendimensionale verstrooiing: integraalvergelijkingsmethode in het tijddomein door R.R.W. de Jager EM-14-94
Verslag van een stage, verricht in de vakgroep EM, onder leiding van prof.dr. A.G. Tijhuis, in de periode december 1993 - oktober 1994.
Eindhoven, 19 oktober 1994.
Samenvatting Om het eendimensionale verstrooiingsprobleem van een gepulste golf die invalt op een inhomogene dielektrische plak op te lossen, gaan we uit van de Maxwell vergelijkingen. Om deze vergelijkingen op te lossen schrijven we ze in contrastvorm, waarbij we gebruik maken van de magnetische en elektrische contraststroombronnen. M.b.v. deze contrast vergelijkingen is het mogelijk voor het verstrooide veld integraalrepresentaties te vinden. De integraalrepresentaties voor het elektrisch en het magnetisch veld worden bepaald door Greense funkties en gewichtsfunkties. Hier zijn de Greense funkties de oplossingen van een excitatie met eenheidsbronnen op het systeem en zijn de gewichtfunkties contraststroombronnen, waarmee we het systeem in feite exciteren. Deze contraststroombronnen wekken namelijk het verstrooide veld op. Voor punten in de plak kunnen de integraalrepresentaties worden opgevat als twee gekoppelde integraalvergelijkingen voor het elektrische en magnetische veld. Beide veldgrootheden worden uitgedrukt in veldwaarden in de hele plak op voorgaande tijdstippen. Dit betekent dat ze kunnen worden opgelost met behulp van de zogenaamde 'marching-on-intime' methode. In deze methode wordt zowel de ruimte als de tijd gediscretiseerd. In het geval van de plak komt dit neer op het invoeren van een ruimte-tijd rooster, het benaderen van een integraal over de dwarsdoorsnede van de plak met behulp van een samengestelde trapeziumregel en een tijdsafgeleide met een drie-punts achterwaartse interpolatieformule. Voor elk ruimte-tijd punt resulteert deze procedure in twee lineaire vergelijkingen. Voor punten in de plak blijven op een natuurlijke manier twee ontkoppelde vergelijkingen over, een lineaire vergelijking voor de elektrische veldsterkte en een voor de magnetische. Voor punten op de rand van de plak moet een twee bij twee systeem worden opgelost. In aIle gevallen volgen beide velden uit zowel elektrische als magnetische veldsterkten in de hele plak op eerdere tijdstippen. Met behulp van een recurrente betrekking is het mogelijk de rekentijd aanzienlijk te verkorten. Deze recurrente betrekking en de besproken 'marching-on-time' methode zijn tenslotte toegepast in een computerprogramma. Om het computerprogramma te testen zijn een aantal configuraties doorgerekend. Zo kunnen we aan de hand van een bepaalde waarde voor de materieparameters de doorlooptijd van de invallende golf controleren, daar deze bekend is. Aan de hand van een verliesvrije homogene configuratie is ook de nauwkeurigheid getest.
1
Inhoudsopgave 1 Inleiding
3
2 Formulering van het probleem
4
3
Afleiding van de integraalvergelijking
4
N umerieke implementatie
7
15
5 N umerieke resultaten
25
6
34
Conclusies
A Het computerprogramma verklaard
36
B Het computerprogramma
38
C Berekening van het elektrisch veld
44
2
Hoofdstuk 1 Inleiding Binnen de vakgroep EM loopt er bij Prof.dr.A.G.Tijhuis het projekt, 'scattering of waves in one dimension'. In dit projekt worden vier verschillende methoden toegepast om het veld te berekenen in een inhomogene, verliezende dielektrische laag in vacuum. Dit veld wordt veroorzaakt door een gepulste vlakke golf die loodrecht invalt. Het probleem is eendimensionaal en is een van de eenvoudigste voorbeelden van propagatie van elektromagnetische golven in inhomogene media. De methoden zijn numerieke oplostechnieken, waarbij er twee het verstrooide veld via het tijddomein en de andere twee het veld via het frequentiedomein bepalen. De technieken in het tijddomein zijn de volgende: • De 'finite-difference time domain' methode. Deze oplostechniek maakt gebruik van ruimte- en tijddifferenties. Het is een lokale methode, waarbij het veld in een bepaald punt in de laag wordt bepaald uit omliggende bekende veldwaarden op eerdere tijdstippen. • De 'marching-on-in-time' methode. Deze techniek maakt gebruik van integraalrepresentaties. Ret is een globale methode, waarbij het veld in een bepaald punt wordt bepaald uit voorgaande reeds bekende veldwaarden in de hele laag. De bedoeling van dit projekt is nu om de studenten ervaring op te laten doen met het gebruik van numerieke oplosmethoden. Verder kan de ontwikkelde programmatuur toegepast worden in een simulatie, waarbij de verstrooiing van de golven en het gedrag van de golf in de laag m.b.v. een computer zichtbaar kan worden gemaakt. Deze simulatie zou dan gebruikt kunnen worden in een demonstratiemodel, waarmee de studenten thuis op hun eigen PC het verschil in de oplostechnieken en het gedrag van de golven bij verschillende inhomogene lagen kunnen bestuderen. De oplostechniek die ik tijdens mijn stage heb uitgewerkt is de 'marching-on-in-time' methode.
3
Hoofdstuk 2 Formulering van het probleem
E(z,t)
f\
Figuur 2.1: Invallende golf op inhomogene plak We bekijken een gepulste golf van een eindige duur die invalt op een vlakke, inhomogene dielektrische plak. De plak met een dikte L, bevindt zich in vacuum (zie fig. (2.1)). Ret invallend elektrisch veld, £i en het bijbehorende magnetisch veld, 1{i zien er als voIgt uit:
{i = F
iii = Rier is F
(t -
~) Uy
-YoF (t
(t -
(2.1)
Co
- :,,) Ux
(2.2)
:J een naar rechts lopende gepulste golf in de z-richting,
Co
de snelheid in
u u
vacuum, Yo = y&/1-0 de karakteristieke admittantie van vacuum en x, y eenheidsvectoren in resp. de x- en y-richting. De funktie F(t) is integreerbaar en verschilt aIleen van nul in het interval 0 < t < T. In de figuur (2.1) aangegeven gebieden met VI, V 2 , V 3 zijn gebieden met verschillende materiaal eigenschappen zoals vermeld in tabel (2.1). Wanneer we op deze configuratie een golf laten invallen zijn we benieuwd naar het veld in de plak, het doorgelaten veld en het gere:fl.ecteerde veld. We zullen hiervoor integraalvergelijkingen a:fl.eiden, waarbij we uitgaan van de volgende Maxwell vergelijkingen: V' x { = -J-L(z)8t
ii
(2.3) 4
Domein VI V2 V3
z-coordinaten -00 < Z < 0 O
Permittiviteit
c(Z)=CO c(Z) = c2r(Z)cO c(Z) = Co
Permeabiliteit J-l(Z) = J-lo J-l(Z) = J-l2 r (Z) J-lo J-l(Z) = J-lo
Geleidbaarheid a(z) = 0 a(z) = a2(Z) a(z) = 0
Tabel 2.1: Verdeling van configuratie in domeinen (2.4) Rierin is V' = (8x , 8y,8z ), de gradient operator. Ret invallend veld is gedefinieerd als de oplossing van deze Maxwell vergelijkingen voor het speciale geval waar c(z) = co, a(z) = 0, J-l(z) = J-lo, d.w.z. in afwezigheid van de plak. Ret totale elektromagnetisch veld wordt nu:
t=
Uyey(z, t)
(2.5)
if
ux 1ix (z, t)
(2.6)
=
Voor het totale elektromagnetisch veld in de verschillende domeinen, weergegeven in tabel (2.1), kunnen we het volgende opschrijven:
VI
ey(z, t) = eY1 (z, t) = ei(z, t)
V2
ey(z, t)
= eY2 (z, t)
V3
ey(z, t)
= ey3 (z, t) = et(z, t)
+ er(z, t)
Soortgelijke vergelijkingen gelden voor 1ix (z, t). De superscripts r en t staan voor resp. het gereflecteerde (reflected) en het doorgelaten (transmitted) veld. Substitutie van vergelijking (2.5) in (2.3) levert:
Ux uy Uz 8x 8y 8z
o ey
0
:::} 8z ey = J-l(z)8t1i x
(2.7)
Substitutie van (2.6) in (2.4) geeft soortgelijk:
Ux uy Uz 8x 8y 8z 1ix 0 0
= 8z 1ix uy - 8y1ix u z = 8z 1ix uy = {a(z)
+ c(z)8deyuy (2.8)
5
Om de integraalvergelijkingen eenvoudig af te leiden passen we op de overgebleven scalaire vergelijkingen (2.7) en (2.8) een enkelzijdige Laplace transformatie naar de tijd toe. Voor het elektrisch veld is deze gedefinieerd als:
J 00
Ey(z, s) =
exp (-st)£y(z, t)dt
to(z)
Hierin is s E
at
azE y = M(Z)sHx
(2.9)
azHx = {o-(z) +c(z)s}Ey
(2.10)
Om deze vergelijkingen op te lossen schrijven we ze in contrastvorm: (2.11) (2.12) Hier is Kx(z, s) de secundaire magnetische stroomdichtheid en Jy(z, s) de secundaire elektrische stroomdichtheid. Beide stroomdichtheden representeren de invloed van het di(~lek trische medium in gebied '0 2 • M.b.v. de vergelijkingen (2.11) en (2.12) zuIlen we in het volgende hoofdstuk de integraalvergelijkingen afleiden.
6
Hoofdstuk 3 Afleiding van de integraalvergelijking Het invallend veld is het veld dat zou bestaan wanneer er geen plak aanwezig is. Voor het invallend veld gelden dus de veldvergelijkingen in vacuum: (3.1) (3.2) De oplossingen van deze vergelijkingen worden gegeven door de Laplace getransformeerden van (2.1) en (2.2). Deze getransformeerden zijn:
E~ ( Z, s) = exp (- ~: ) F (s ) H~(z, s)
= -Yoexp
(-~:)F(s)
Hierin is: T
F(s)
=
J
exp (-st):F(t)dt
o
Ais we deze getransformeerden in de veldvergelijkingen (3.1) en (3.2) invullen zien we direkt dat ze voldoen. Om het gerefiecteerde, het doorgelaten en het veld in de plak bepalen, vatten we het totale veld op als een superpositie van een invallend veld en een verstrooid veld. Dit verstrooide veld loopt zowel in gebied VI als in gebied V 3 van de plak weg (zie fig. (2.1)). Een representatie van het verstrooide veld krijgen we door van de vergelijkingen (2.11) en (2.12) de corresponderende veldvergelijkingen in vacuum af te trekken:
8z E; - /.LosH: = Kx(z, s)
(3.3)
8z H: - cosE;
(3.4)
=
Jy(z, s) 7
In deze vergelijkingen is E; = Ey - E~ en H~ = H x - H~; hierin staat de superscript s voor het Engelse woord 'scattered', wat in het Nederlands 'verstrooid' betekent. De subscripts x en y in de elektromagnetische velden zullen we in het vervolg weglaten. We kunnen nu m.b.v. vergelijking (3.3) en (3.4) het verstrooide veld bepalen, dat veroorzaakt wordt door de contrast stroombronnen K(z, s) en J(z, s). Deze bronnen staan als het ware in de vrije ruimte te stralen en wekken het verstrooide veld op. Hierdoor voldoet het veld aan de uitstralingsvoorwaarden. In de afieiding van de integraalrepresentatie beschouwen we de bronnen eerst als bekend. Vervolgens passen we het superpositiebeginsel toe. Veronderstel dat we een elektromagnetisch veld (Et, Ht) exciteren met een magnetische stroombron K 1 en een veld (E2, H2) met een magnetische stroombron K 2 • Dan voIgt uit de lineariteit van de vergelijkingen (3.3) en (3.4) dat een magnetische stroombron K 1 + K 2 het veld (Et + E 2, Ht + Hn genereert. Hetzelfde geldt voor combinaties van twee elektrische stroombronnen en van een elektrische en een magnetische stroombron. Op grand van deze conclusie is het voldoende om te kijken naar de excitatie van twee eenheidsbronnen, een magnetische en een elektrische. Ais de bijbehorende responsies C.q. velden, de zogenaamde Greense funkties, bekend zijn kan het complete verstrooide veld vervolgens geschreven worden als een lineaire combinatie van deze oplossingen. Hiertoe schrijven we eerst de contrast stroombronnen in een andere vorm:
J=K(z', s)8(z - z')dz' JK(z', s)8(z - z')dz' L
K(z, s)
J(z, s)
=
=
=
-=
0
=
L
J J(z', s)8(z - z')dz' = JJ(z', s)8(z - z')dz' -=
(3.5)
(3.6)
0
Uit de formules (2.11) en (2.12) voIgt dat de contrast stroombronnen alleen in het domein V 2 (zie tabel (2.1)) van nul verschillen. In vergelijking (3.5) en (3.6) zijn K(z, s) en J(z, s) nu als continue som 1 van resp. de gewichtsfunkties, K(z', s) en J(z', s) en de eenheidsbron, 8(z - z'), geschreven. Nu we de bronnen als een 'continue' som hebben geschreven, mogen we ook de velden volgens het superpositiebeginsel als 'continue' som schrijven. Vervolgens gaan we het systeem exciteren met de twee eenheidsbronnen en bepalen hieruit de Greense funkties. We exciteren nu alleen met de magnetische of elektrische eenheidsbron 8(z - z'), waarbij we de gewichtsfunkties K(z', s) en J(z', s) voor het gemak weglaten. Daar het systeem zich lineair gedraagt kunnen we hiermee altijd nog achteraf vermenigvuldigen. Uit de vergelijkingen van het verstrooide veld, namelijk (3.3) en (3.4), voIgt bij excitatie met een magnetische eenheidsbron:
(3.7) 1
We kunnen een integraal zien als een continue som.
8
De superscripts E en H in deze Greense funktie geven resp. het soort veld en de soort excitatie weer. In dit geval wordt geexciteerd met een magnetische eenheidsbron als elementaire bouwsteen van K(z, s) en dus van een magnetische contrastbron evenredig met H. CE,H is een elektrisch veld en CH,H een magnetisch veld. Bij excitatie met een elektrische eenheidsbron voIgt:
ozCE,E - SJ.LoCH,E = 0 { ozCH,E - SEoCE,E = b(z - z')
(3.8)
Als we de velden tengevolge van deze excitaties met eenheidsbronnen in de vorm van Greense funkties weten, kunnen we m.b.v. het superpositiebeginsel een integraalrepresentatie voor het elektrisch veld afieiden, waarin de Greense funkties voor het E-veld, CE,E en CE,H nog met de betreffende gewichtsfunkties, K(z', s) en J(z', s) moeten worden vermenigvuldigd (zie formule (3.5) en (3.6)). L
ES(z,s) =
L
JCE,E(Z, z'; s)J(z',s)dz' + JCE,H(Z, z'; s)K(z', s)dz' o
(3.9)
0
De integraalrepresentatie voor het magnetisch veld kunnen we op dezelfde manier afieiden. L
HS(z,s) =
L
JCH,E(Z,Z';s)J(z',s)dz' + JCH,H(Z,Z'; s)K(z', s)dz' o
(3.10)
0
Deze integraalrepresentaties zijn geldig voor -00 < z < 00. Rest ons aIleen de Greense funkties uit stelsel (3.7) en stelsel (3.8) te bepalen. Uit stelsel (3.7) voIgt:
O;CH,H - SEoOzCE,H H,H _ 02C z
s2 c~
= O;CH,H -
SE o{SJ.LoCH,H
CH,H = SE 0 b(z - z')
+ b(z -
z')}
=0 (3.11)
Voor z =j:. z' geldt in formule (3.11):
02C H,H _ £CH,H = 0 z
(3.12)
c~
De algemene oplossing van deze tweedegraads homogene differentiaalvergelijking luidt als voIgt:
CH,H = Aexp (-~:) + Bexp (~:) We zien dat de oplossing bestaat uit een naar rechts en een naar links lopende golf. De Greense funktie CH,H geeft de oplossing van een magnetische eenheidsbron (contrastbron),
9
die op het punt z = z' in de vrije ruimte staat te stralen. Vanwege de uitstralingsvoorwaarden 2 kunnen we nu voor de funktie GH,H het volgende sehrijven:
z > z'
GH,H
A exp ( - s(zc~z'»)
z < z'
GH,H
Bexp CCz:,:z'»)
(3.13)
We zien dus voor z < z' een naar links lopende golf, en voor z > z' een naar reehts lopende golf. Uit het ongerijmde kunnen we eoncluderen dat GH,H in formule (3.12) eontinu 3 moet zijn. Hieruit voIgt dat B = A. We kunnen nu (3.13) hersehrijven als: GH,H
= Aexp (_s1zc~z'l)
(3.14)
Om de eonstante A te berekenen gaan we vergelijking (3.11) integreren van z' - 6z naar z' + 6z en nemen we de limiet voor 6z 1 o. z'+L'l.z
J
lim
L'l.zl0
z'-L'l.z
lim
f)z L'l.zl0 {
z'+L'l.z
{riGH,H - s2 GH,H} dz = lim z C~ L'l.zl0
GH,H Z'+L'l.Z} - 0 = Iz'-L'l.z
J
Sc a
8(z - z')dz
z'-L'l.z
SCa
Wanneer we in deze vergelijking formule (3.14) substitueren, geeft dit: lim
L'l.zl0 { _2s Co
IZ'+L'l.Z} f)zA exp ( - ~) Co z'-L'l.z
A=
SCa ::::}
A=
-~ 2
=
=
SCa
-~ 2
Hieruit voIgt: GH,H = _ ~ exp (_ slz-z'l) 2
(3.15)
Co
Om nu de funktie GE,H te bepalen substitueren we funktie (3.15) in de tweede regel van stelsel (3.7). _l_f)zGH,H
=
_XsLf)z exp (_ slz-z'l)
SCo
2sco
- ..YsLf)z exp (- s(z-z'») 2sco Co
Co
voor z > z'
~ exp ( - s(zc~z'») _XsLf)z exp (sCz-z'») 2sco
Co
voor z
< z'
_1. exp (s(z-z'») 2
Co
2De oplossing bestaat uit een naar links lopende golf voor z
< z'
en een naar rechts lopende golf voor
z > z' en voldoet hiermee aan de uitstralingsvoorwaarden. 3De afgeleide van een discontinue funktie levert een 8-funktie. Daar het rechterlid van for mule (3.12) nul is kan alleen een continue funktie de juiste oplossing zijn.
10
We kunnen dus de Greense funktie CE,H
CE,H
verkort schrijven als:
= lsgn(z - z') exp (_SIZ-ZII) 2
(3.16)
~
Hierin is sgn(z) de zogenaamde tekenfunktie:
~
= {
sgn(z)
voor z > 0 voor z = 0 voor z < 0
-1
(3.17)
De Greense funkties bij excitatie met de eenheidsbron van de elektrische stroomdichtheid kunnen we op een zelfde manier afleiden. Wanneer we echter gebruik maken van het dualiteitsprincipe 4 kunnen we deze op een veel gemakkelijkere wijze bepalen. Met formule (3.15) en het dualiteitsprincipe voIgt: CE,E
=
_~ 2
exp (_ slz-Z11)
(3.18)
Co
Hierin is Zo = ~o de karakteristieke impedantie van vacuum. Met vergelijking (3.16) en het dualiteitsprincipe voIgt: CH,E
z') exp (_ s/z-zll)
= lsgn(z -
(3.19)
~
2
De Greense funkties zijn nu allemaal bekend. We kunnen nu de integraalrepresentaties (3.9) en (3.10) als voIgt schrijven: L
ES(z,s)
=
Jexp (_s1zc~zll)J(z',s)dz' +k Jsgn(z - z') exp (- s1zc~zll)K(z', s)dz' _~o
o
L
(3.20)
o L
kJsgn(z -
z') exp ( - s1zc~z'l) J(z', s)dz'
o L
-If Jexp (- s1zc~zll)K(z', s)dz'
(3.21)
o
Hierin substitueren we de contrast stroomdichtheden gedefinieerd in (2.11) en (2.12):
K(z, s) J(z, s)
= [{L(z) - {Lo]sH(z, s) = {Lo[{Lr(Z) - l]sH(z, s) =
{a(z)
+ Eo[varepsilonr(z) -
l)s]E(z, s)
4Bij het dualiteitsprincipe maken we gebruik van de eigenschap dat in de veldvergelijkingen van Maxwell de volgende grootheden mogen worden verwisseld: E - H, Zo - Yo, Co - /-Lo' Wanneer we stelsels (3.7) en (3.8) met elkaar vergelijken zien we dit dualiteitsprincipe naar voren komen.
11
De integraalrepresentaties worden dan: L
ES(z, s)
=
_co;o J c't') + [Er(Z') -
1]s} exp ( - s1zc~z'l )E(z', s)dz'
o L
+T J
sgn(z - Z')[ILr(Z') - 1]sexp
(_s1zc~z'I)H(z', s)dz'
(3.22)
o L
HS(z,s)
=
-Jl:°ro J[ILr(Z') -1]sexp (_s1zc~z'I)H(z',s)dz' o L
+T J
sgn(z - z')ct')
+ [Er(Z') -1]s}exp (_s1zc~z'I)E(z',s)dz'
(3.23)
o
Hierin is Xe(z) = Er(Z) - 1 de elektrische susceptibiliteit, en Xm(z) = ILr(Z) - 1 de magnetische susceptibiliteit. Om de integraalrepresentaties eenvoudiger te maken normeren we aIle afstanden naar de lengte van de plak, en aIle tijden naar de tijd die nodig is om de plak te doorlopen met de snelheid in vacuum. We krijgen dan de volgende dimensieloze grootheden :
z=
~
L'
t
= sJ. s = sL L ' Co
Voor de elektromagnetische velden geldt:
E(z, s)
= E(z, s), fI(z, s) = ZoH(z, s)
De eenheid van fI is hier gelijk aan de eenheid van E = ~. In termen van deze genormeerde grootheden kunnen de integraalrepresentaties worden geschreven als: 1
ES(z, s)
=
T Jsgn(z -
z')Xm(,Z')~ exp (_sco:~~-z'I)YofI(z', s)Ldz'
o 1
_co;o J{YoZ;:') + Xe(z')~} exp (_sco:~~-z'I)E(z', s)Ldz'
(3.24)
o 1
T Jsgn(z -
z'HYoZ}:')
+ Xe(z')~} exp (_Sco:~~-z'I)E(z', s)Ldz'
o 1
- Jl:°ro
JXm(z')~
exp ( - sco:~~-z'I)YofI(z', s)Ldz'
o
12
(3.25)
In het vervolg zuBen we de normeerstrepen in de formules weglaten om schrijfwerk te vermijden. Wanneer we integraalrepresentatie (3.24) en (3.25) verder uitwerken, krijgen we : 1
ES(z, s)
=
J +~ Jsgn(z - z')Xm(z')sexp(-slz - z'I)H(z',s)dz' -~ {a(z') + Xe(z')s} exp (-slz - z'I)E(z', s)dz' o
1
(3.26)
o 1
JXm(z')s exp (-slz - z'I)H(z', s)dz' +~ Jsgn(z - z'){a(z') + Xe(z')S} exp (-slz - z'I)E(z', s)dz'
HS(z, s) = -~
o
1
(3.27)
o
Deze uitdrukkingen zien er gelukkig al een stuk eenvoudiger uit. Het enige dat ons nog te doen staat is het terugtransformeren naar het tijddomein. De complexe hoekfrequentie s gaat bij terugtransformatie over in een differentiatie naar de tijd 8t . De e-macht in bovenstaande integraalrepresentaties levert bij terugtransformatie een tijdsvertraging op. De vertraagde tijd bedraagt:
t' = t - Iz - z'l Voor de integraalrepresentaties, waarbij we t:S(z, t) en 1i S(z, t) schrijven als resp. t:(z, t) t:i(z, t) en 1i(z, t) - 1ii (z, t) geldt nu: 1
t:(z, t)
=
t:i(z, t) - ~
J{a(z') + Xe(z')8dt:(z', t')dz' o
1
+~
Jsgn(z - z')xm(z')8 1i(z', t')dz' t
(3.28)
o 1
1i(z, t)
=
1ii (z, t) - ~
JXm(z')8 1i(z', t')dz' t
o 1
+~
Jsgn(z - z'){a(z') + Xe(z')8 }t:(z', t')dz' t
(3.29)
o
Voor z < 0 en z > 1 is het veld nu uitgedrukt in termen van het totale veld in de plak. Wanneer we dus het veld in de plak weten, kunnen we hieruit het verstrooide veld bepalen. Het gerefiecteerde veld 1iT (z, t) voor z < 0 komt overeen met het rechterlid van 13
de integraalrepresentatie (3.29) min het invallende veld 1{i(z, t), terwijl het doorgelaten veld 1{t(z, t) voor z > 1 overeen komt met totale veld 1{(z, t). Voor [r(z, t) en [t(z, t) geldt een soortgelijke afleiding uit integraalrepresentatie (3.28). M.b.v. het voorgaande en het gegeven dat t' = t -Iz - z'l kunnen we het volgende afleiden:
1{r(z, t) = 1{r(o, t + z) 1{t(z, t) = 1{t(1, t - z + 1)
voor z < 0 voor z > 1
Voor het verstrooide elektrisch veld gelden soortgelijke vergelijkingen. Het verstrooide veld op een bepaalde plaats, z buiten de inhomogene plak en tijdstip t is dus gelijk aan het verstrooide veld op de randen van de plak op het tijdstip waarop de bijbehorende golf de plak verlaten heeft. Voor 0 :::; z :::; 1 en 0 :::; t < 00 gaan de integraalrepresentaties over in integraalvergelijkingen. M.b.v. deze integraalvergelijkingen kunnen we het veld in de plak bepalen. In figuur (3.1) zien we de plaats-tijd punten die van belang zijn voor het bepalen van het veld op een bepaalde (z, t) in de plak.
(z,t)
Figuur 3.1: Relevante plaats-tijd punten We zien in deze figuur, dat het veld in de plak op een bepaalde plaats en tijd gevonden kan worden uit de velden op voorgaande plaats-tijd punten. Om nu het veld op een bepaalde (z, t) te berekenen, moeten we al deze veldwaarden integreren en differentii:~ren. Dit lijkt in eerste instantie een heel karwei. Echter de numerieke wiskunde kan ons hierbij uitkomst bieden. In het volgende hoofdstuk zullen we dan ook de integraalvergelijkingen (3.28) en (3.29) numeriek gaan oplossen.
14
Hoofdstuk 4 N umerieke implementatie De integraalvergelijkingen (3.28) en (3.29) zullen eerst verder uitgewerkt worden om de numerieke implementatie duidelijker te maken. In het bijzonder brengen we de tekenfunktie sgn(z - z') in rekening: 1
J{O"(z') + Xe(z')&t}£(z', t')dz' +~ JXm(z')&/H(z', t')dz' - ~ JXm(z')&/H(z', t')dz' o
£(z, t) = £i(Z, t) - ~
o
z
1
(4.1)
z
1
H(z, t) = Hi(z, t) - ~
JXm(Z')&tH(Z', t')dz' o 1
z
J
J
+~ {O"(Z') + Xe(z')&d£(z', t')dz' - ~ {O"(Z') + Xe(z')&d£(z', t')dz' (4.2) o
z
We hebben in deze integraalvergelijkingen de integralen met een tekenfunktie (zie (3.28) en (3.29)) gesplitst in een postitief gedeelte voor z' :S z en een negatief gedeelte voor z' 2:: z. We mogen in deze integralen het punt z' = z meenemen, daar de waarde van de integrand in een punt geen invloed heeft op de waarde van de integraal. Om nu het stelsel integraalvergelijkingen (4.1) en (4.2) numeriek op te lossen, discretiseren we ruimte en tijd. Hiertoe voeren we het rooster (mh, nl::::.t) in, waarbij m = 0,1,2, ... ,N en n = 0,1,2, ... ,00 zijn met stapgrootten h = l::::.t = ~ en N het aantal stappen in de plak. We kunnen nu de integraal in deze vergelijkingen discretiseren met een samengestelde trapeziumregel. De tijdsafgeleide in deze vergelijkingen zal gediscretiseerd worden met een drie-punts achterwaartse interpolatieformule. Eerst zullen de samengestelde trapeziumregel en de interpolatieformule afgeleid worden. De samengestelde trapeziumregel kunnen we voor een willekeurige twee maal differentieerbare
15
funktie f(z) op de volgende manier afleiden: N-I (n+l)h
I
Jf(z)dz
=
o
J f(z)dz
L n=O
nh
t
Hier is h = de stapgrootte, N het aantal stappen. De integraal over een enkel subinterval kunnen we nu benaderen met een trapeziumregel: N-I (n+l)h
L n=O
J f(z)dz
N-I
L
=
nh
~{f((n + 1)h) + f(nh)} + O(h 2 )
n=O
Als h 1 o. Dit wordt de samengestelde trapeziumregel genoemd. De fout die hierbij optreedt is gelijk aan O(h 2 ). De funktie moet dan weI twee maal differentieerbaar zijn. In ons geval is de funktie f(z) inderdaad twee maal differentieerbaar. Als we de trapeziumregel nog verder uitwerken, krijgen we: N-I
L
N-I
~{f((n
+ 1)h) + f(nh)}
= ~{L
n=O
N-I
f(nh)
+
n=O N-I
= ~{L
N
f(nh)
+
n=O
L
L
n=O
f((n
+ 1) h)}
~ m
N-I
f(mh)}
= ~{f(O)
+ 2 L f(nh) + f(!:!}J)}
m=1
n=1
I
N-I
= ~f(O)
+ h L f(nh) + ~f(1)
(4.3)
n=1
Om de tijdsafgeleide te discretiseren, gaan we uit van Taylorreeksen. In figuur (4.1) zijn de punten weergegeven, die we nodig hebben voor een Taylor reeksontwikkeling. In deze
f(t)
At t,
At t,
to
t
Figuur 4.1: Ligging punten reeksontwikkeling. figuur is f(t) een drie maal differentieerbare, continue funktie. Taylorontwikkeling van f(tl) rond to geeft:
f(t l ) = f(t o) + (-6.t)J'(t o) + (-~t? !"(to) + O((6.t)3)
16
Taylorontwikkeling van f(t 2) rand to:
f(t 2) = f(t o) + (-26t)f'(t o) + (-2~t)2 f"(t o) + O((6t)3) Om f'(t o) te bepalen kiezen we een zodanige lineaire combinatie van f(t 1 ) en f(t 2) dat daaruit de tweede afgeleide f"(t o ) wegvalt:
4f(t 1 )
-
f'(t o) =
f(t 2) = 3f(t o) - 26tf'(t o) + O((6t)3)
L Hf(t o) -
2f(t 1 )
+ ~f(t2)} + O((6t)2)
(4.4)
We hebben nu in feite een drie-punts achterwaartse interpolatieformule gevonden met een fout gelijk aan O((6t)2). Daar de integraalvergelijkingen (4.1) en (4.2) mooi differentieerbaar zijn en de derde afgeleide bestaat mogen we deze door gebruik te maken van de formules (4.3) en (4.4) discretiseren. Als we met een dubbele tijdstap werken, kunnen we bovendien de rekentijd met een factor twee reduceren. We leggen dan een beperking op aan het aantal plaats-tijd punten en alleen de punten zoals in figuur (4.2) weergegeven zijn nog van belang. Deze plaatstijd punten (mh, n6t) zijn punten met m + n even. Hierin zijn m = 0,1,2, ... , N en n = 0,1,2, ... ,00. We kunnen nu de tijddiscrete integraalvergelijking voor het elektrisch
Figuur 4.2: Plaats-tijd punten in het discrete geval veld afleiden, waarin de tijdsafgeleide door een drie-punts achterwaartse interpolatieformule met dubbele tijdstap 6t = 2h wordt benaderd: N
E[m, n]
= Ei[m, n] -
~
L
Wm,
o-[m']E[m', n']
m'=O
17
N
-~ L
Wm, Xe [m'] AH£[m', n'] - 2£[m', n' - 2] + ~£[m', n' - 4]}
m'=O m-l
+~
L
Wm, Xm[m'] AHil[m',n'] - 2il[m',n' - 2] + ~il[m',n' - 4]}
m'=O
+ 41hDm,N WN N
-~
xm[N]{~il[N, n] - 2il[N, n - 2]
+ ~il[N, n - 4]}
Wm, Xm[m'] AHil[m',n'] - 2il[m',n' - 2] + ~il[m',n' - 4]}
L m'=m+l
- 4\ Dm,O Wo Xm[O]{~il[O, n] -
2il[0, n - 2] + ~il[O, n - 4]}
(4.5)
De tijddiscrete integraalvergelijking voor het magnetisch veld wordt:
il[m, n]
=
ili[m, n] N
-~
L m'=O m-l
+~
L
Wm, Xm[m'] AHil[m', n'] - 2il[m', n' - 2] + ~il[m', n' - 4]} Wm, Xe [m'] A{~£[m', n'] - 2£[m', n' - 2] + ~£[m', n' - 4]}
m'=O
+ 41h Dm,N WN Xe[N]{ ~£[N, n] - 2£[N, n - 2] + ~£[N, n - 4]} N
-~
L
Wm, Xe[m'] AH£[m', n'] - 2£[m', n' - 2] + ~£[m', n' - 4]}
m'=m+l
- 4\ Dm,o Wo Xe[O]{ ~£[O, n] - 2£[0, n - 2] + ~£[O, n +~ L Wm, iT[m']£[m', n'] + ~Dm,N WN iT[N]£[N, n]
4]}
m-l
m'=O N
-~
L
Wm, iT[m']£[m', n'] - ~Dm,o Wo iT[O]£[O, n]
(4.6)
m'=m+l
In (4.5) en (4.6) is Dm,r het Kronecker symbool: I
Dm,r = { 0
alsm=r als m =1= r
(4.7)
n' = n - 1m - m'l, £i[m, n] = £i(mh, n6t), ili[m, n] = Hi(mh, n6t), Xe[m] = Xe(mh), Xm[m] = Xm(mh), iT[m] = CJ(mh), Wm = h voor 0 < m < N en Wo = WN = ~. In vergelijking (4.5) is een sommatie van m' = 0 naar m-l en een sommatie van m' = m+l naar N te zien. Uit deze twee sommaties is de waarde op het punt m' = m tegen elkaar weggevallen, waardoor de bovengrens van de ene sommatie gelijk aan m-l en de ondergrens van de andere sommatie gelijk aan m+ 1 is geworden. Uitzonderingen hierbij zijn de punten m = 0 en m = N. Op deze twee punten valt er geen waarde meer uit de sommaties weg, daar een van de twee sommaties wordt uitgesloten. Hierdoor zijn we genoodzaakt deze waarde via de Kronecker delta Dm,r toe te voegen. Een zelfde redenering geldt voor de tijddiscrete 'magnetische' integraalvergelijking (4.6). 18
We onderscheiden nu twee gevallen. In het eerste geval, voor 0 < m < N geldt dat de Kronecker symbolen 8m ,N en 8m ,o gelijk aan nul zijn. De tijddiscrete elektrische integraalvergelijking (4.5) bevat nu een onbekende, de E[m, n]. In het rechterlid van deze vergelijking heeft E[m, n] een gewichtsfactor gelijk aan :
De tijddiscrete 'elektrische' integraalvergelijking is nu eenvoudig oplosbaar, daar we een vergelijking hebben met een onbekende. Dit geldt ook voor de tijddiscrete 'magnetische' integraalvergelijking (4.6), waarbij de gewichtsfactor in het rechterlid gelijk is aan:
WH[m]
hWm Xm[m]
=
_g3
In het tweede geval, op de punten m = 0 en m = N komt er een tweede onbekende bij in de tijddiscrete integraalvergelijkingen. Dit is voor de tijddiscrete 'elektrische' integraalvergelijking, het magnetisch veld op het punt (m, n). Deze H[m, n] heeft in het rechterlid van vergelijking (4.5) een gewichtsfactor gelijk aan:
WH[O]
=
_g3h
Wo Xm[O]
voor m =
-WH[N] = +g3h WN Xm[N]
voor m
°
=N
Bij de tijddiscrete 'magnetische' integraalvergelijking krijgen we ook een tweede onbekende, het elektrisch veld op het punt (m, n). Deze waarde E[m, n] heeft in het rechterlid van vergelijking (4.6) een gewichtsfactor gelijk aan:
WE[O]
= -~wo
{o-[O] + 43hXe [O]}
-WE[N] = +~WN {o-[N] + 43hXe [N]}
voor m = voor m
°
=N
We hebben in dit geval te maken met twee vergelijkingen in twee onbekenden, waaruit het elektromagnetisch veld te berekenen is. Schrijven we nu in de vergelijkingen (4.5) en (4.6) de somtermen van m' = 0 naar m' = N als voIgt: N
L
m'=O
N
m
L+ L
m'=O m'=m+l
dan kunnen we uit de tijddiscrete 'elektromagnetische' vergelijkingen de volgende definities voor de somtermen afleiden: m
S;[m,n]
~
L
m'=O
Wm' o-[m']E[m', n -
m
+ m']
N
S;[m,n]
~
L
m'=m+l
Wm' o-[m']E[m', n + m
19
-
m']
m
S~n1[m, n]
AL
Wm, Xm[m']il[m', n -
m
+ m']
m'=O N
A L
Wm'Xm[m']il[m',n+m-m']
m'=m+l m
S~e[m,n]
L
4~
W
m, Xe[m']£[m', n -
m
+ m']
m'=O N
S~e[m,n]
L
l 4h
W
m, Xe[m']£[m', n + m
-
m']
m'=m+l
Bij deze afleiding is de definitie n' = n - 1m - m'l gebruikt. In figuur (4.2) zijn de met schuine lijnstukken verbonden plaats-tijd punten weergegeven, die van belang zijn bij het berekenen van het elektrisch of het magnetisch veld op een bepaald punt (m, n). Deze lijnstukken aangegeven met de letters a tim f komen overeen met de volgende somtermen:
a : S~(m -
1, n
-
1], S~Jm
-
1, n
-
1], S~Jm
-
1, n
-
1]
b : S~e [m, n - 2], S~n1 [m, n - 2] c : S~e [m, n - 4], S~n1 [m, n - 4]
d: S';[m, n], S~Jm, n], S~Jm, n] e : S~e [m, n - 2], S~n1 [m, n - 2]
f : S~Jm, n - 4], S~Jm, n - 4] M.b.v. deze definities en de gewichtsfactoren kunnen de tijddiscrete elektromagnetische integraalvergelijkingen als voIgt herschreven worden. Om de integraalvergelijkingen overzichtelijk te houden, onderscheiden we drie gevallen. De elektromagnetische vergelijkingen in de plak, voor 0 < m < N, en de vergelijkingen op de randen van de plak, voor m = 0 en m= N. 1. De elektromagnetische vergelijkingen in de plak, voor 0 < m < N:
£[m,n]
=
{1- WE[m]}-l {£i[m,n] -~(S~e
[m - 1, n - 1] + S~e [m, n] - S~n1 [m - 1, n - 1] + S~n1 [m, n]) +2(S~Jm, n - 2] + S~Jm, n - 2] - S~Jm - 1, n - 3] + S~Jm, n - 2]) -HS~e [m, n - 4] + S~Jm, n - 4] - S~Jm - 1, n - 5] + S~Jm, n - 4]) -S~[m - 1, n - 1] - S';[m, nJ} (4.8)
20
H[m,n]
= {1- WH[m]}-l {Hi[m,n] -HS~",[m -1,n -1] + S~",[m,n] +2(S~",[m, n - 2] + S~",[m, n - 2] -HS~",[m, n - 4] + S~",[m, n - 4] -
S~e[m -1,n -1]
+ S~e[m, n]) S~Jm - 1, n - 3] + S~Jm, n - 2]) S~Jm - 1, n - 5] + S~Jm, n - 4])
+S,;[m - 1, n - 1] - S;[m, n]}
(4.9)
= 0: = tiro, n] - HS~JO, n] + S~",[O, n]) +2(S~e [-1, n - 1] + S~", [-1, n _1(S2 [-1 n - 3] + S2x",, [-1 n 2 Xe'
2. De elektromagnetische vergelijkingen op de rand, voor m
{I - We[O]} £[0, n] - WH[O]H[O, n]
1]) 3])
-S;[O, n]
{I - WH[O]}H[O, n] - We[O]£[O, n]
=
(4.10)
Hi[O, n] - HS~",[O, n] + S~JO, n]) +2(S~",[-I, n - 1] + S~J-l, n - 1]) _1(S2 [-1 n - 3] + S2X[-1 n - 3]) 2 x",, e' -S;[O, n]
(4.11)
3. De elektromagnetische vergelijkingen op de rand, voor m = N:
{I - WE [N]}£[N, n]
+ WH[N]H[N, n]
=
£i[N, n] _.1(Sl [N - 1 n - 1] - SlX'" [N - 1' n - 1]) 2 Xe' +2(S~JN, n - 2] - S~",[N, n - 2])
_l(Sl [N ' n - 4] - SlX'" [N , n - 4]) 2 Xe -S';[N - 1, n - 1]
{1- WH[N]}H[N,n]
+ We[N]£[N,n]
=
(4.12)
Hi[N,n] _.1(Sl 2 x[N "-, 1 , n - 1] - SlXe [N - 1' n - 1]) +2(S~",[N, n - 2] - S~JN, n - 2])
_l(Sl [N n - 4] - Sl [N n - 4]) 2
X'"
'
+S,;[N - 1, n - 1]
Xe
'
(4.13)
In vergelijking (4.8) of (4.9) zien we voor 0 < m < N een onbekende. We kunnen dus uit vergelijking (4.8) het elektrisch veld op het punt (m, n) bepalen en uit vergelijking (4.9) het magnetisch veld. Voor m = 0 of m = N komen er echter twee onbekenden Elm, n] en H[m, n] in een vergelijking voor. Om nu het elektrisch en magnetisch veld te bepalen, zullen we de twee vergelijkingen voor m = 0 of m = N simultaan moeten oplossen. Om deze afieiding eenvoudiger te maken, definieren we de volgende somtermen. Uit het rechterlid van vergelijking (4.10) voIgt voor m = 0:
SEo = tiro, n] - ~( ... ) + ... - S;[O, n] 21
Uit vergelijking (4.11) volgt voor m = 0: SHo = j{i[O, n] - ~( ...) + ...
-
S;[O, n]
Uit vergelijking (4.12) volgt voor m = N: SEN
= Ei[N,n] -
~( ...) + ...
-
S;[N -1,n -1]
En tenslotte volgt uit vergelijking (4.13) voor m = N: -i
3
SHN = 1-l [N,n] - 2"("')
1 + ... + Sq[N -1,n -1]
Uit de vergelijkingen (4.10) en (4.11) volgt het volgende stelsel voor m = 0:
{1 - WE [O]}EJO, n] - WH[O]~[O, n] = SEo { {1 - WH[O]} 7-l[0, n] - WE[O]£[O, n] = SHo Uit stelsel (4.14) wordt nu het elektrisch veld voor m
(4.14)
= °bepaald.
Voor dit veld geldt:
{1- WH [0]}{1- WE[O]}E[O,n] - WH[O]WdO]E[O,n] = {1 - WH[O]}SEo
+ WH[O]SHo
Hieruit volgt: E[O n] _ {1- WH[O]}SEo + WH[O]SHo ,1 - WH[O] - WElD]
(4.15)
Voor het magnetisch veld geldt uit stelsel (4.14) een soortgelijke afleiding: j{[0 n] = {1 - WE[O]}SHo + WE [0] SEo , 1- WH[O] - WElD]
(4.16)
Voor m = N voIgt uit de vergelijkingen (4.12) en (4.13) ook een stelsel, waaruit het elektromagnetisch veld bepaald kan worden:
{1 - WE [N]}EJN, n] + WH[N]~[N, n] { {1- WH[N]}7-l[N,n] + WE[N]£[N,n]
= SEN = SHN
(4.17)
Uit stelsel (4.17) wordt nu het elektromagnetisch veld voor m = N bepaald. Voor het elektrisch veld geldt:
{1- WH [N]}{1- WE[N]}E[N,n] - WH[N]WE[N]E[N,n]
=
{l- WH[N]}SE N - WH[N]SHN Hieruit voIgt: E[N n] _ {1- WH[N]}SEN - WH[N]SHN , 1- WH[N] - WE[N]
22
(4.18)
Voor het magnetisch veld geldt:
il[N n]
{I - WE[N]}SHN - WE[N]SEN 1- WH[N] - WE[N]
=
,
(4.19)
Met deze uitdrukkingen (4.15), (4.16), (4.18) en (4.19) zijn we nu in staat om op de randen van de plak het elektromagnetisch veld te berekenen. Ret elektromagnetisch veld in de plak (voor 0 < m < N) voIgt uit de vergelijkingen (4.8) en (4.9). Ret elektromagnetisch veld is nu uitgedrukt in voorgaande elektrische en magnetische veldwaarden. Dit wordt de 'marching-on-in-time' methode genoemd. am de rekentijd nog verder te beperken, bewaren we de somtermen berekend in het punt (m - 1, n - 1). We kunnen dan bij de berekening van het veld in het punt (m, n) de somtermen weergegeven door de schuine lijnen a,b en c (zie fig. (4.2)) weer opnieuw gebruiken. Retzelfde geldt voor somtermen berekend in het punt (m + 1, n - 1). Deze somtermen zijn in figuur (4.2) weergegeven door de schuine lijnen d,e en f. De somtermen berekend in het punt (m - 1, n - 1) moeten we bij hergebruik in het punt (m, n) weI eerst aanpassen. We moeten namelijk de invloed van het elektromagnetisch veld op het punt (m - 1, n - 1) er eerst nog bijtellen. We leiden hiervoor de volgende recurrente betrekkingen af, die volgen uit de definities van de somtermen: m
S~[m, n] = ~
L
Wm,
o-[m']£[m', n -
m
+ m']
Wm ,
o-[m']£[m', n -
m
+ m'] + ~o-[m]£[m, n]
m'=O m-l
~
L
m'=O
= S~[m
- 1, n - 1] + ~o-[m]£[m, n]
N
~
S;[m,n]
L
Wm ,
o-[m']£[m', n + m
-
m']
Wm ,
o-[m']£[m', n + m
-
m'] + ~o-[m + l]£[m + 1, n - 1]
m'=m+l
N
~
L m'=m+2
=
S; [m + 1, n -
1]
+ ~o- [m + 1]£[m + 1, n -
1]
Vergelijkbare recurrente betrekkingen gelden voor S~e' S~e' S~m en S~m' Op de randen van de plak geldt het volgende: S~[-l, n]
S;[N, n]
= S~J-1, n] = S~J-1, n] = 0
= S~JN, n] = S~JN, n] = 0
Door gebruik te maken van deze recurrente betrekkingen kunnen we de rekentijd aanzienlijk verkorten. Stel we hebben te maken met N punten in de lengte van de plak en aN punten in de lengte van het tijdinterval, dan moeten we bij ieder van deze aN2 plaats-tijd punten een som uitrekenen van N punten als N --+ 00. Dit komt overeen met een rekentijd in de orde
23
van N 3 als N ----+ 00. Wanneer we echter gebruik maken van de recurrente betrekkingen, hoeven we bij ieder van deze plaats-tijd punten de som alleen met een punt uit te breiden. Dit zou een rekentijd in de orde van N 2 vergen. In het volgende hoofdstuk zijn de resultaten te zien van de 'marching-on-in-time' methode, toegepast in een computerprogramma. In dit programma zijn de hier besproken recurrente betrekkingen gebruikt om de rekentijd te verkorten.
24
Hoofdstuk 5 N umerieke resultaten Het computerprogramma, waarmee de numerieke resultaten behaald zijn is opgenomen in appendix B. In appendix A wordt de werking van het programma verklaard. Om het programma grondig te testen zijn we in eerste instantie uitgegaan van een homogene plak. Hierbij hebben we het voordeel, dat voor een homogene plak het doorgelaten en gereflecteerde veld analytisch berekenend kunnen worden (zie appendix C). We kunnen dus bij de homogene plak de numerieke resultaten toetsen en hieruit de nauwkeurigheid afschatten. Vervolgens zijn we gaan kijken naar verscheidene inhomogene configuraties. H
1 rt";-------::-------------;:==:::;-,
r,., !\'" t,I~i'" I E
,
I
oI \ I
1\
'
\, ,
I
1. -I,
"
IMlo'll\ H(O,I)
I
-0.5
I I
·0.5
.
I
1..,
.
I , \
I , I,
H\1,ll\
H(l,II"
-I '--------'--'- - - - - - - - - - - - - - - - ' 2 3 4 5 6 7 8 9 10 o
,
I \ -1 ' - - - - - - - - - - - - - - - - - - - - - ' 2 3 4 5 6 7 8 9 10 0
t-.
t -----+
Figuur 5.1: Invallende sinuskwadraat puis op homogene plak De homogene plak waar we in eerste instantie van uit gaan is een verliesvrij medium met de volgende materiaaleigenschappen: Xe = 1,25, Xm = 0, a = O. In figuur (5.1) is het elektromagnetisch veld op de randen (z = 0 en z = 1) van de plak te zien bij een invallende sinuskwadraat puIs, :F = sin2(~)PT(t). Hier is PT een rechthoekige puIs gedefinieerd door:
PT(t)
=
{~
voor t < 0 voor 0 < t < T voor t > T
(5.1)
25
We nemen hier een sinuskwadraat puIs, omdat dit signaal een mooie gerefl.ecteerde en doorgelaten golf opwekt. Toch treedt bij deze sinuskwadraat puIs al een lichte overshoot op. Deze is in figuur (5.1) te zien op het moment dat de doorgelaten golf (£t(l, t) of 1{t(l, t)) de plak verlaat. In figuur (5.1) is aIleen het invallend elektrisch veld weergegeven. Ret invallend magnetisch veld is weggelaten, daar het aIleen van teken verschilt van het elektrisch invallend veld (zie formule (2.2). T geeft de tijdsduur van de puIs weer. In figuur (5.1) is T = 1 en het aantal stappen in de plak N = 60. Dit geldt ook voor de rest van de behaalde numerieke resultaten. Ook is in deze figuur het verstrooide veld weergegeven van een invallende sinuskwadraat puIs op een plak met de volgende materieparameters: Xe = 0, Xm = 1,25, (j = O. Bij deze configuratie is de rol van Xe en Xm omgedraaid. Als we de figuur aandachtig bekijken dan zien we voor het doorgelaten veld, £t(l, t) en 1{t(l, t) geen verandering optreden. Ret gerefl.ecteerde veld, t:r(0, t) en 1{T(O, t) veranderd daarentegen weI van teken. De amplitudes van de golven blijven gelijk. In deze twee figuren geldt dat het ingaande vermogen gelijk is aan het uitgaande vermogen. We hebben immers te maken met een verliesvrije plak ((j = 0). We kunnen nu de nauwkeurigheid van de numerieke resultaten uit figuur (5.1) afschatten. Deze afschatting maken we aIleen voor de verliesvrije plak met Xe = 1,25 en Xm = O. We berekenen hierbij de 'Root-Mean-Square', oftewel de RMS-fout. Deze is voor het gerefl.ecteerde elektrisch veld als voIgt te bepalen:
RMS-fout
=
(5.2)
Rierin is £T(O, t) het numeriek benaderde veld en t:r(0, t) het exact berekende veld. Voor het gerefl.ecteerde en het doorgelaten elektrische veld zijn in tabel (5.1) voor verschillende waarden van N de RMS-fouten bepaald. De RMS-fouten van het magnetisch veld zijn niet bepaald, daar deze fouten even groot zullen blijken te zijn. Dit voIgt uit de symmetrie van de integraalrepresentaties en de bijbehorende gediscretiseerde versies. We kunnen dit ook conc1uderen aan de hand van figuur (5.1). Ret gerefl.ecteerde elektrisch en magnetisch veld zijn namelijk even groot en het doorgelaten veld verschilt aIleen van teken. In tabel (5.1) zijn tevens een aantal karakteristieke waarden weergegeven. Rierin is r de doorlooptijd, die een golf nodig heeft om van z = 0 naar z = 1 te komen. In de configuratie weergegeven in figuur (5.1) is r = 1,5. Deze r voIgt eenvoudig uit het funktievoorschrift van de invallende golf in de plak. Deze funktie is evenredig met:
T12 F
(t - c:)
Rierin is T12 de transmissiecoefficient van gebied 1 (een vacuum) naar gebied 2 (de plak) (zie tabel (2.1). Wanneer we deze funktie normeren, voIgt: j:
(t -
S!.~) L C2
=
j:
(t - S!.z) C2
26
i r (0,27 + f)
it (1,7 + f)
RMS-
fout (%) exact N= 10 N=20 N=40 N=60 N=80 N= 100
37,4 21,6 7,81 3,84 2,29 1,53
0,192 0,0955 0,1547 0,1869 0,1912 0,19158 0,19188
RMS-
fout (%)
Rekentijd Trek (ms)
38,5 18,1 5,93 2,90 1,73 1,16
33 111 408 932 1561 2420
0,96 0,4550 0,8615 0,9492 0,9564 0,95858 0,95925
Tabel 5.1: De RMS-fouten De doorlooptijd 7 is de tijd, die deze golf nodig heeft om van Deze 7 is gelijk aan: 7
=
~=
vcr/1r = V(Xe
z=
°naar z = 1 te komen.
+ l)(Xm + 1)
Hierin is 7 een dimensieloze grootheid. De doorlooptijd is afhankelijk van de materieparameters Cr en /1r. De karakteristieke waarden zijn de lokale maxima. Het is leuk om in tabel (5.1) te zien, hoe bij hogere N deze maxima steeds dichter bij de exaete waarden komen. Als we naar de RMS-fout kijken zien we een afname bij grotere N, maar de afwijking is echter groter dan de procentuele afwijking tussen de exacte maxima en de numeriek bepaalde. Dit heeft te maken met het feit dat bij deze maxima de numeriek bepaalde waarde dichter bij de exacte waarde ligt. Er waren ook waarden, waarbij de fout groter was. De grootste fout doet zich voor bij de 'staart' van de direkt propagerende puIs. De verklaring hiervoor is dat op deze tijdstappen de tweede afgeleide naar de tijd van £(z, t) en H(z, t) discontinu is. Dit betekent dat de afleiding die leidt tot formule (4.4),niet langer geldig is. Aangezien dit slechts om enkele 'lijnen' in het (z, t) domein gaat blijven de veldwaarden op zich echter weI betrouwbaar. Van de kolom rekentijd Trek in tabel (5.1) hebben we de tijden logaritmisch uitgezet tegen het aantal stappen in de plak N. De grafiek is te zien figuur (5.2). Elke keer dat we de discretisatie met een factor 2 verfijnen (dus N twee maal zo groot) blijkt de rekentijd ongeveer met een factor 4 toe te nemen. Dit is in figuur (5.2) duidelijk te zien daar ~ log(Trek ) tegen log(N) is uitgezet. De hellingshoek a van de lijn is dan ook tan(a) = 0,93 is ongeveer gelijk aan 45°. Dit bevestigt de conclusie, dat de rekentijd door gebruik te maken van de recurrente betrekkingen evenredig is met N 2 • De rekentijden weergegeven in de kolom rekentijd van tabel (5.1) zijn behaald met een PC (80486DX, 33 MHz). In figuur (5.3) is voor Xe = 3 en Xe = 15 het verstrooide veld bepaald. We kunnen voor deze twee gevallen de doorlooptijd, 7 als voIgt berekenen: Xe = 3 Xm =
° /1r Cr
= 4 } :::} = 1
7
=
.J4 =
2
27
1.75..,-----------------------."
1/2'0 log(T,ek)
f'''
1.35
.
1.15
0.95
0.75 - f ' - - - - - - - . - - - - - - - , - - - - . - - - - - - - - I 1.75 1.25 1.0 1.50 2.00
'Olog(N) ---+
Figuur 5.2: Rekentijd Trek logaritmisch uitgezet tegen N Xe = 15 Xm = 0
= 16 } J.Lr = 1
Cr
=} T
=
Vi6 =
4
In figuur (5.3) is de doorlooptijd met een verticale gestippelde lijn weergegeven. We zien dat op deze tijdspunten de doorgelaten golf inderdaad de plak verlaat. Bij Xe = 15 is zoals verwacht de langzaamste tijd waarneembaar. In deze figuur is te zien dat bij Xe = 15 de overshoot in het verstrooide veld (in de 'staarten' van de golven) toeneemt. Dit heeft te maken met het feit dat het signaal voor een hogere Xe meer vertraagd wordt, waardoor de golven in de plak als het ware meer worden samengeknepen. Hierdoor neemt de fout in het berekende verstrooide veld toe. Om dit te vermijden zouden we de discretisatie in de plak moeten verfijnen. In figuur (5.4) is een leuke situatie weergegeven. De materieparameters Xe en Xm zijn gelijk aan elkaar. In dit geval treden er geen refiecties op. Omdat we te maken hebben met een verliesvrije plak (cr = 0) is de amplitude van de doorgelaten golf gelijk aan de invallende golf. De T in deze configuratie is gelijk aan 2. We kunnen dit gedrag theoretisch als voIgt verklaren. Hierbij gaan we uit van een twee lagen medium, om de verklaring te vergemakkelijken. Het twee lagen medium bestaat uit een vacuum in domein 1 voor -00 < z < 0 en een verliesvrije halfruimte met materieparameters C2 en J.L2 in domein 2 voor 0 < z < 00. Voor de refiectiecoefficient geldt nu:
Hierin is ZOI = V~ de karakteristieke impedantie van domein 1, een vacuum en Z02 = V!iii go g2 de karakteristieke impedantie van domein 2. Willen er nu geen refiecties optreden dan moet R = 0 zijn. Dit is het geval wanneer Z02 = ZO!' Hieruit voIgt:
28
H
1
:',
1 I, 0.5 ;.
~
I ' ,
,:
I I I,
" 1"
0\ I.
:\
V,H'IO,I)
-0.5
E(O,I)
/\
!. ~ ~(O~t)
0.5
o
'-,
I
,
I
,
I
\
. . . • . • • . . .:. ,':'" • • . • • . . . • . . . . . .a.
,' ,
EII,I)
'
,:
I I I H'(l,l) '1' i .,.
I \ . . . . . . . . . . . . . . '1' I HIO,I)
1 , I
.
H
E
I
I
I
e'iO,I)
.
~~ .
1
'\ .
,\
I
5
, .-
"
I,
,,
~:;~
, 1
a =0
• ,
•... :. ,: . ',.-E\1'!l ........•..••••..........
E
I,
IXe=31 %,,;:0
I I ,,',
I I E\O,l)
\ I
I I
,
J . .. ~(~,I)• . . . • .
_ •• _ . • . _ •
\I
I
-1 1---+--+--+-+--1--+--1---+--+_---1 -11---+-+--+-+--+-+--I--+--1-----l 2 3 4 5 6 7 8 9 10 2 3 4 5 6 7 8 9 10 0 o t---+ t---+
Figuur 5.3: De doorlooptijd
7
afhankelijk van Xe
Dus in het geval van gelijke materieparameters Xe = Xm treden er geen reflecties op. In figuur (5.5) is de invloed van de materieparameter (J weergegeven. Rier is duidelijk zicht baar dat de doorlooptijd 7 niet afhankelijk is van deze materieparameter. Indien het doorgelaten veld nog zichtbaar is, is 7 in deze plaatjes gelijk aan 1. We zien verder dat bij een steeds grotere (J, de amplitude van de doorgelaten golf afneemt. De verliezen worden dus bij een hogere waarde voor (J groter. In het geval van (J = 100 is de doorgelaten golf zelfs helemaal niet meer waarneembaar. Ret invaIlend dringt nog weI een stukje in de plak door, maar wordt uiteindelijk geheel gereflecteerd. Rierbij moet worden opgemerkt dat we het hier hebben over de genormeerde 0-. De werkelijke (J is gelijk aan:
Rierin is L de lengte van de plak en Zo de karakteristieke impedantie van vacuum. In figuur (5.6) is het elektromagnetisch verstrooide veld zichtbaar van een invaIlende sinuskwadraat puIs op een inhomogene verliesvrije plak. In het ene geval hebben we te maken met een Xe = 1 + Z, Xm = 0 en in het andere geval hebben we de materieparameters weer verwisseld om de werking van het computerprogramma te testen. In aIle twee de gevaIlen hebben we te maken met een lineair profiel. We zien, zoals ook bij het homogene geval (zie fig. (5.1)), dat het doorgelaten veld niet en het gereflecteerde aIleen van teken verandert. In figuur (5.7) is nog een leuke situatie weergegeven. Ook hier hebben we te maken met een inhomogene verliesvrije plak, maar nu hebben Xe en Xm aIlebei een lineair verloop. Ret verloop van deze materieparameters is echter zodanig dat de een maximaal is waar de ander minimaal is, en omgekeerd. Wanneer we nu het gerefl.ecteerde veld, £T(O, t) en 1{T(O, t) in figuur (5.7) vergelijken met het gereflecteerde veld in figuur (5.1) dan kunnen we het verloop als voIgt verklaren. In figuur (5.7) is het gerefl.ecteerde veld in eerste instantie positief. Dit komt overeen met het gerefl.ecteerde veld voor Xm = 1,25 en Xe = 0 in 29
',
H
X
11
:
r
,
I
~ E(O,t)
IX.,,=1 e=11
: :: E(1,I)
C1 =0
I
,
E 0.5:···~··········:··~·········· I ,
, I
I
I
I ,
, I
,
I
.
01'------'---<----..;.,....--------------1
,
\ \
I
I I
I ,
1 1
, ,
I
I
······i··j··
-0.5
•
I
,
I
I
'
.
I
H(1,t)I , 1 I
-1
L-
o
---'
"
~
2
3
4
5
6
7
8
9
t
10
----+
Figuur 5.4: Geen gerefiecties bij gelijke materieparameters Xe = Xm figuur (5.1), daar de materieparameter Xm in figuur (5.7) voor z = 0 groter is. De tweede positieve 'bult' in het gerefl.ecteerde veld kunnen we op analoge wijze verklaren. De Xe is voor z = 1 maximaal en daarom is het gerefl.ecteerde veld nu weer positief. Vergelijk met figuur (5.1) voor Xe = 1,25 en Xm = O. Bij een plak met een inhomogene (7 is het soort profiel, als (7 niet al te groot is, af te lezen uit het gerefl.ecteerde veld. In figuur (5.8) zijn als voorbeeld de verstrooide elektromagnetische velden van een invallende sinuskwadraat puIs op drie verschillende profielen weergegeven. Bij ieder van deze profielen is nu een karakteristiek gerefl.ecteerd veld waarneembaar, waaruit het soort profiel kan worden afgelezen. Een zelfde redenering kunnen we maken voor een plak met een inhomogene Xe' In dit geval is het gerefl.ecteerde veld ongeveer evenredig met de afgeleide van Xe(z). De resultaten zijn te zien in figuur (5.9). Wanneer we deze resultaten vergelijken met figuur (5.8) dan kunnen we concluderen dat het gerefl.ecteerde veld in figuur (5.9) voor het kwadratisch profiel overeenstemt met het gerefl.ecteerde veld in figuur (5.8) voor het lineaire profiel. Ret soort profiel kan nu ook weer uit het gerefl.ecteerde veld worden afgelezen. Nu levert echter een kwadratisch profiel weergegeven door de materieparameter Xe een lineair verloop in het gerefl.ecteerde veld op en een lineair profiel een constant verloop (zie fig. (5.9)). Er zijn nog talloze verschillende soorten configuraties bedenkbaar. Ropelijk zijn de hier gegeven voorbeelden voldoende om een algemene indruk te krijgen over het functioneren van het gebruikte algoritme, de 'marching-on-in-time' methode. De numerieke resultaten, zoals in dit hoofdstuk weergegeven, hebben allemaal betrekking op het verstrooide veld. In hoofdstuk 2 hadden we echter gezegd dat we ook het veld in de plak wilden weten. Wat doen we nu met deze veldwaarden in de plak? Met deze waarden kunnen we bijvoorbeeld een real-time simulatie van de golf in het eendimensionale verstrooiingsprobleem maken, waarbij ook het veld in de plak zichtbaar kan worden gemaakt.
30
H
~=o x..=o H
1
I,\
t1
I E(O,') I \ .. El",) .. \ .. ,
lo.5 ,
-
-
.
I . \'
,
E
-
.,
=1
I
~
1 'Iit~g \
0'
E
,
\.
~
1
o
o .......
,
\
•• e'(l.Q
"
....
\~ E(O~O,')
.
-0.5
-1
L-
-1
--'
o
2
3
4
5
6
7
8
9
- -
L-
o
10
.
---'
2
3
4
8
6
H
E
I~~ I
,\
!-
0.5
_\~'(~,')__ .. ,
' ,
\ \
_
~.=.2~.
_
E\"I)
J\(.. ~ L-
o
.1
--'
2
3
.
. ~:\f:
E(O,I)
.1
10
t-.
t-.
1
9
4
6
7
8
9
L-
o
10
---'
2
3
4
6
t-.
7
8
9
10
t-.
Figuur 5.5: Invloed van
(J'
op het verstrooide veld
Zulke simulaties geven spectaculaire beelden over het verstrooien van een invallende golf op een plak.
31
H
1
I
0.5 ,'.
E
1X.=1+Z x..,: 0 I
IIII
~~'(9") . ; . :.E(l,~
'I
,,
I
o
·0.5
a =0
I
_ ..
1 '.1~E'lD'~) i.~ H
1~m=1+Z
I
E
0.5 ,;
1\
I
I~.=o
,\
.
,
I
I
I
Ell,!)
~ ~~ ..
-
i"{
.
Olll-~==-""'.);--'-,--;;,~""----------l
~D,I)\ HID,I) "
"
I
I
t
.~
,
--
- . - - - •...•.•.•....•••
\
I
I
,
I
,
,E(D,t) H(D,I)
I , I . f"
-0.5 .
,
,Hll,l)
I'
-
.
,H(1.1)
I, II
II
-1 ' - - - - - - - - - - - - - - - - - - - - - - ' 8 9 10 o 2 4 5 6 7 3
-1 ' - - -
o
t ----+
,
.....J
4
3
2
5
6
7
8
9
10
t ----+
Figuur 5.6: Invallende sinuskwadraat puls op inhomogene plak
"
I ,
X.=1 +z
I I
Xm =2-z
, , I
:
I
(J
=0
'IE(O,t)
0.5 ','" \',." ..... , I
,
I I I
I I ,
,~"
H(0~(0,1)
0F------'.1=====i;====t-=-==~L.----"'-----------j '.
-0.5
I
I
I
I
I
I
:
I
I
.•.....•.•........ l, ...•. , •...••••••••••••.••••••....•..•••••••••••••••
,
'H(1,1)
I
I
I,
'
I
,'
-1
L-
a
-',,_'- - - - - - - - - - - - - - - - - - - '
2
3
4
5
6
7
8
9
10
Figuur 5.7: Inhomogene verliesvrije plak met twee lineaire materieparameters
32
KWADRATISCH PROFIEL
L1NEAIR PROFIEL
1.-,-----------;===;1
IfJ I ', . E i .:"
!:,.,
HI
, I~Oll ::
E 0.5 ,.' I
E 0.5 "' I..•.' .' ........••.••••••••
,, . , ,: 0l~, . ~ I
I
.'
,
•
~
1
~I~
"
1 I
.•
,
I .:
:.
'I
I ; H
:EioA ;':
.
'.
, I
1!'""T"-----------.====O
Ilr.=, %,r0 2 0' =1'4(1·0,5)
'1
H
EIO,Q :: EII,Q
SINUSKWADRAAT PROREL
1,.-------------.::;=:'=9
1\
,
I
:
I
I
.
o
f\
.
H~Qe'~,11 :
,
I I
J~
.
.(l.5 •••
I,
J~
..
I, I,
I,
~I,Q'
~l.Q'
·1'--------------------'
1
.
'I .'
I ,
o
~1»
i . : .:
E 0.5 ,.' I
H~II 0.11, ,
.(l.5
r.=1
%,r0 0' =sinlz,zj
2
3
4
5
6
7
8
9
10
t--+
·1'--------------------' 012345678910
·1'--------------------'
012345678910
t--+
t--+
Figuur 5.8: Zichtbaar maken van het profiel in de plak met inhomogene
L1NEAIR PROFIEL
1,------..----------------.===
H
"
lO.5 f E
II
.'t·
'I I I
1X.=1+zl
'.
!
~~
.. '>\1;11.•••...•••••••••••••
f\
I
\
KWADRATISCH PROFIEL
H
-0.5
,I
1 E
,
I
I,Xe=2'4(Z'0.5f
Xm=o
.,
. .
1..::"_=_0_ _---'1
:":~11:11
I,
__ .. . ..
.
.
:~./\ 0"1\-r"=,,;...........,."""'--"--'-""-------~ I
l.f.YElo,l) M(O.Q
I
I
..
I I EtO.I)
I
ok"\----:/1=";'---:=/........'-----'------.:':...;,----------i HiVE(O.I) (O.I)
1r - - < - - - - - - - - - - r : ; " = = = j 1 J
0.5 : .:
••
(7
I
1 I
. . • . • . 1..'
-0.5
1. .1. . . . . . . . . . . . . . . . . • . . . • . . . . .
I
1 IH\1.I)
.
I IH\1.') \ 1
I I
,I
II
J
•
-1 ' - - - - - - - - - - - - - - - - - - - - - - - ' 2 3 4 5 6 7 8 9
-1'-----------------------'
t ------.
t ------.
o
10
o
2
3
4
5
6
7
8
9
Figuur 5.9: Zichtbaar maken van het profiel in de plak met inhomogene Xe
33
10
Hoofdstuk 6 Conclusies Aan de hand van de resultaten, die gepresenteerd zijn in hoofdstuk 5, is te concluderen dat het computerprogramma naar behoren werkt. Zo hebben we de doorlooptijd van de invallende elektromagnetische golven gecontroleerd door bepaalde waarden voor de materieparameters Xe en Xm te nemen. Deze doorlooptijd moet nu gelijk zijn aan JCrf.Lr en daaraan werd voldaan. Ook is gekeken naar de invloed van de materieparameter (j. We weten dat bij (j = 00 de plak ondoordringbaar wordt voor elektromagnetische golven. Ret bleek dat bij hoge waarden van (j de amplitude van de doorgelaten golf bijna nul is en de invallende golf bijna geheel gereflecteerd werd. Het is dus aannemelijk om te zeggen, dat wanneer (j --+ 00 de plak volledig ondoordringbaar zal worden. Ret leuke geval waarbij de materieparameters Xe en Xm gelijk aan elkaar zijn kan ook worden gebruikt om het computerprogramma te testen. Bij gelijke materieparameters treden er namelijk geen reflecties op. Om de nauwkeurigheid van het programma te testen hebben we het verstrooide elektrisch veld van een invallende golf op een homogene verliesvrije plak berekend (zie appendix C). Deze waarden zijn vervolgens vergeleken met de door het computerprogramma bepaalde numerieke waarden en hiermee hebben we de fout kunnen afschatten. Ret bleek dat de fout niet voor alle waarden hetzelfde was. De kleinste fout trad op bij de topwaarden van de golven. Hierdoor zijn de waarden in tabel (5.1) enigszins misleidend, daar hier de RMSfout, de berekende en de numeriek bepaalde topwaarden zijn weergegeven. De RMS-fout is dan ook groter, omdat deze een globale indruk van de opgetreden fout weergeeft. Toch kunnen we zeggen dat de fout bij grotere waarden voor N zodanig klein wordt dat de berekening nauwkeurig verloopt. De toegepaste 'marching-on-in-time' methode in het computerprogramma bij het bepalen van het elektromagnetisch veld levert goede resultaten op. De nauwkeurigheid is niet alleen goed, maar aan de hand van de verschillende testconfiguraties kunnen we ook concluderen dat de methode zeker betrouwbaar is. Verder is de berekening op een PC snel uit te voeren (min. 80386 of 80486). De rekentijd was zoals verwacht evenredig met N2. De numerieke resultaten kunnen verder nog gecontroleerd worden aan de hand van een reeds bestaand werk [1]. Dit werk heb ik tevens gebruikt als voorbeeld bij het behalen van de numerieke resultaten. Verder diende het als leidraad bij het schrijven van dit 34
stageverslag. M.b.v. het bijgeleverde computerprogramma (zie appendix B) zijn aIle veldwaarden voor de verschillende figuren in hoofdstuk 5 bepaald. Dit computerprogramma is echter een subroutine. Wil men de numerieke resultaten controleren, dan zal er een hoofdprogramma gebruikt moeten worden waarin de invallende elektromagnetische velden, de materieparameters en het aantal stappen in de plak worden klaargezet.
35
Bijlage A Het computerprogramma verklaard Het computerprogramma dat te vinden is in appendix B is geprogrammeerd in 'Fortran 77'. Dit programma maakt gebruik van de recurrente betrekkingen, zoals vermeld in hoofdstuk 4. De somtermen die in hoofdstuk 4 vermeld staan komen overeen met de somtermen vermeld in het computerprogramma. Zo is bijvoorbeeld S~e in de tekst gelijk aan SCHIEl in het programma. Bij de berekening van het elektromagnetisch veld lopen we schuin van links naar rechts door de plak. Wanneer we figuur (4.2) bekijken dan starten we in het plaats-tijd punt [0, 0] en eindigen in het punt [N, N]. Vervolgens beginnen we weer in het punt [0,2] en eindigen in het punt [N, N + 2], etc. Hierin is N het aantal stappen in de plak. Dit heeft als voordeel dat we het gedeelte van de plaats-tijd ruimte waar het elektromagnetisch veld nul is (in figuur (4.2) weergegeven met een driehoek) overslaan. Wanneer we horizontaal door de plak van links naar rechts het elektromagnetisch veld zouden berekenen dan moeten we een uitzondering maken voor de plaats-tijd punten met n oneven. Bij deze punten moeten we namelijk beginnen op het plaats-tijd punt met m = 1 en eindigen in het punt met m = N - 1 (zie fig. (4.2)). Wanneer we schuin door de plak lopen is deze uitzondering niet nodig. Het programma is onderverdeeld in een aantal blokken. Met een blok bedoelen we een gedeelte van het programma afgescheiden door asterix's. Boven zo'n blok staat tevens een beknopte beschrijving. Als voorbeeld het eerste blok, 'Trapeziumregel'. In dit gedeelte wordt de trapeziumregel zoals vermeld in fomule (4.3) met de halve waarde op de randen van de plak klaar gezet. De rest van de blokken zijn bijna allemaal in hoofdstuk 4 terug te vinden. Zo is de berekening van het elektromagnetisch veld in het blok, 'Voorkant plaat' terug te vinden in formule (4.15) en (4.16). De somtermen die in het programma worden opgeslagen in de array's onder de namen, SCHIEl, SCHIE2, SCHIM1 en SCHIM2. Deze array's bestaan ieder weer uit drie rijen waar aIle berekende veldwaarden, die voor de recurrente betrekkingen nodig zijn, in worden opgeslagen. In figuur (4.2) is te zien, welke veldwaarden allemaal bekend moeten zijn voor de berekening van het elektromagnetisch veld op een plaats-tijd punt. Wanneer we bijvoorbeeld het elektrisch veld op een plaats-tijd punt [m, n] willen berekenen, dan staan de betreffende somtermen weergegeven met schuine lijnen in de volgende array's: 36
lijnstuk a: SCHIEl(AO,M-l),SCHIMl(AO,M-l) lijnstuk b: SCHIEl(Al,M),SCHIMl(Al,M) lijnstuk c: SCHIEl (A2,M) ,SCHIMl (A2,M) Voor de lijnstukken d,e en f geldt hetzelfde, maar nu met de array's: SCHIE2 en SCHIM2. De variabelen AO,Al en A2 krijgen in het programma een zodanige waarde, dat de waarden in de array SCHIEl(AO,M), berekend bij de punten van bijv. [0,0] tot [N, N], bij de volgende schuine lijn met de punten van [0,2] tot [N, N + 2] in de array SCHIEl(Al,M) terecht komen. De waarden in de array's worden dus doorgeschoven. Hiermee hebben we het voordeel dat we de waarden in de array's niet hoeven te kopieren.
37
Bijlage B Het computerprogramma C
C C C C C
****************************************************** Dit programma is gesehreven door: Remeo de Jager. T.U. Eindhoven, Faeulteit Elektroteehniek, vakgroep Elektromagnetisme (EM), 5 oktober 1994. ******************************************************
C
SUBROUTINE GL(EH,HH,EI,HI,CHIE,CHIM,SIG,NM,NN,HALFDEL) IMPLICIT NONE DOUBLE PRECISION E(O:100,O:500),H(O:100,O:500), & EH(O:50,O:500),HH(O:50,O:500), & EI(O:500),HI(O:500), & SCHIE1(O:100,O:2),SCHIE2(O:100,O:2), & SCHIM1(O:100,O:2),SCHIM2(O:100,O:2), & SSIG1(O:100),SSIG2(O:100) , & WE(O:100),WH(O:100), & CHIE(O:100),CHIM(O:100),SIG(O:100), & SOM1,SOM2,SOM3,SOM4,HALFDEL, & SEO,SHO,SEN,SHN, & WEO,WHO,WEN,WHN, & CO,DO,CN,DN,SOMTO,SOMTN C
INTEGER M,N,NN,NM,HM,L,AO,A1,A2 C
C C C
********************* Trapeziumregel *********************
C
SIG(O)=5D-1*SIG(O) SIG(NM)=5D-1*SIG(NM) 38
CHIE(0)=5D-l*CHIE(0) CHIE(NM)=5D-l*CHIE(NM) CHIM(0)=5D-l*CHIM(0) CHIM(NM)=5D-l*CHIM(NM) C
C
*********************
C
Zelftermen
C
*********************
C
100
DO 100,M=0,NM WE(M)=lDO/(lDO+375D-3*CHIE(M)+HALFDEL*SIG(M)) WH(M)=lDO/(lDO+375D-3*CHIM(M)) CONTINUE
C
C
**********************
C C
Sommen nodig berekening rand plaat
C
**********************
C
WEO=-375D-3*CHIE(0)-HALFDEL*SIG(0) WHO=-375D-3*CHIM(0) CO=lDO-WHO DO=lDO-WEO SOMTO=lDO-WEO-WHO WEN=-375D-3*CHIE(NM)-HALFDEL*SIG(NM) WHN=-375D-3*CHIM(NM) CN=lDO-WHN DN=lDO-WEN SOMTN=lDO-WEN-WHN C
C
*********************
C
Nodig voor berekening
C
*********************
C
200
DO 200,M=0,NM CHIE(M)=CHIE(M)/4DO CHIM(M)=CHIM(M)/4DO CONTINUE
C
C
*********************
C
Sommen op nul zetten
C
*********************
C
39
300
400
DO 400,M=O,NM DO 300,L=O,2 SCHIM1(M,L)=ODO SCHIM2(M,L)=ODO SCHIE1(M,L)=ODO SCHIE2(M,L)=ODO CONTINUE SSIG1(M)=ODO SSIG2(M)=ODO CONTINUE
C
C C C
********************* Initialisatie ARRAYS *********************
C
AO=2 A1=1 A2=0 C
C C C
********************* Berekening *********************
C
DO 600,N=O,NN C
C C C
********************** Voorkant plaat **********************
C
SOM2=SCHIM2(1,AO)+SCHIE2(1,AO) SOM3=SCHIM2(O,AO)+SCHIE2(O,AO) SOM4=SCHIM2(O,A1)+SCHIE2(O,A1) C
&
SEO=EI(N)-15D-1*SOM2 +2DO*SOM3-5D-1*SOM4-SSIG2(1)
&
SHO=HI(N)-15D-1*SOM2 +2DO*SOM3-5D-1*SOM4-SSIG2(1)
C
C
E(O,N)=(CO*SEO+WHO*SHO)/SOMTO H(O,N)=(DO*SHO+WEO*SEO)/SOMTO C
SCHIE1(O,AO)=CHIE(0)*E(O,N)
40
SCHIE2(0,A2)=SCHIE2(1,AO)+CHIE(0)*E(0,N) SCHIM1(0,AO)=CHIM(0)*H(0,N) SCHIM2(0,A2)=SCHIM2(1,AO)+CHIM(0)*H(0,N) SSIG1(0)=HALFDEL*SIG(0)*E(0,N) SSIG2(0)=SSIG2(1)+HALFDEL*SIG(0)*E(0,N) C C C C C
******************* Plaat ******************* DO 500,M=1,NM-l SOM1=SCHIE1(M-l,AO)-SCHIM1(M-1,AO) SOM2=SCHIM2(M+l,AO)+SCHIE2(M+1,AO) SOM3=SCHIM2(M+l,Al)+SCHIE2(M+1,Al) SOM4=SCHIM2(M+l,A2)+SCHIE2(M+1,A2)
C
& & &
E(M,N)=(EI(N)-15D-l*(SOM1+S0M2) +2DO*(SOM3-SCHIM1(M-l,Al)+SCHIE1(M,Al)) -5D-l*(SCHIE1(M,A2)-SCHIM1(M-1,A2)+SOM4) -SSIG1(M-l)-SSIG2(M+l))*WE(M)
& & &
H(M,N)=(HI(N)-15D-l*(SOM2-S0Ml) +2DO*(SOM3+SCHIM1(M,Al)-SCHIE1(M-1,Al)) -5D-l*(SCHIM1(M,A2)-SCHIE1(M-1,A2)+SOM4) +SSIG1(M-l)-SSIG2(M+l))*WH(M)
C
C
500 C C C C C
SCHIE1(M,AO)=SCHIE1(M-1,AO)+CHIE(M)*E(M,N) SCHIE2(M,A2)=SCHIE2(M+1,AO)+CHIE(M)*E(M,N) SCHIM1(M,AO)=SCHIM1(M-1,AO)+CHIM(M)*H(M,N) SCHIM2(M,A2)=SCHIM2(M+1,AO)+CHIM(M)*H(M,N) SSIG1(M)=SSIG1(M-l)+HALFDEL*SIG(M)*E(M,N) SSIG2(M)=SSIG2(M+l)+HALFDEL*SIG(M)*E(M,N) CONTINUE ********************** Achterkant plaat ********************** SOM1=SCHIE1(NM-l,AO)-SCHIM1(NM-1,AO)
C
& &
SEN=EI(N)-15D-l*SOMl +2DO*(SCHIE1(NM,Al)-SCHIM1(NM,Al)) -5D-l*(SCHIE1(NM,A2)-SCHIM1(NM,A2))
41
&
-SSIG1(NM-l)
C
& & &
SHN=HI(N)+15D-l*SOMl +2DO*(SCHIM1(NM,Al)-SCHIE1(NM,Al)) -5D-l*(SCHIM1(NM,A2)-SCHIE1(NM,A2)) +SSIG1(NM-l)
C
E(NM,N)=(CN*SEN-WHN*SHN)/SOMTN H(NM,N)=(DN*SHN-WEN*SEN)/SOMTN C
SCHIE1(NM,AO)=SCHIE1(NM-l,AO)+CHIE(NM)*E(NM,N) SCHIE2(NM,A2)=CHIE(NM)*E(NM,N) SCHIM1(NM,AO)=SCHIM1(NM-l,AO)+CHIM(NM)*H(NM,N) SCHIM2(NM,A2)=CHIM(NM)*H(NM,N) SSIG1(NM)=SSIG1(NM-l)+HALFDEL*SIG(NM)*E(NM,N) SSIG2(NM)=HALFDEL*SIG(NM)*E(NM,N) C
C
**********************
C
Doorschuiven ARRAY's
C
**********************
C
AO=AO+l IF(AO.GT.2)THEN AO=O ENDIF Al=Al+l IF (Al. GT .2)THEN Al=O ENDIF A2=A2+1 IF(A2.GT.2)THEN A2=0 END IF C
600
CONTINUE
C
C
**********************
C C
Recht trekken van de plaats-tijd punten
C
**********************
C
L=O DO 800,M=0,NM,2 42
700 800
HM=M/2 DO 700,N=0,NN IF (N-L.LT.O) THEN EH(HM,N)=O HH(HM,N)=O ELSE EH(HM,N)=E(M,N-L) HH(HM,N)=H(M,N-L) END IF CONTINUE L=L+l CONTINUE RETURN END
43
Bijlage C Berekening van het elektrisch veld We gaan uit van een homogene, verliesvrije plak in vacuum. In dit geval kunnen we het verstrooide veld als voIgt berekenen. Voor het invallend veld geldt:
j} = F(w)uy exp( -~) Dit voIgt ook uit de Fourier-getransformeerde van formule (2.1). Voor het gereflecteerde veld geldt dan:
ET =
R l1 F(w)uy exp(~)
= ETuy
Hierin is R l1 de reflectiecoefficient [2]. Het gereflecteerde veld in het tijddomein kunnen we nu vinden door deze vergelijking terug te transformeren naar de tijd:
£T =
2~
JETexp(jwt)dw = 2~ J 00
00
-00
-00
R l1 F(w)exp(jw{t+ :))dw
00
2~
R-Rexp(-
J JOO{
1- R2 exp( _
2jwL )
~ F(w)exp(jw{t+: })dw
-00
R 271"
1 R2
-
0
F(w)
'li!:!l. exp(Jw{t + ~ }) .
Z
1- R2 exp( _ )
-00
_
)
c2
0
C2
F(w). zwL) exp (2 - C2
exp(jw{t + .L
_
Co
2L})} dw C2
Hierin is :
R is de reflectiecoefficient voor een golf die vanuit medium 1 loodrecht invalt op medium 2. Dit wordt ook weI genoteerd als R 1-+ 2 •
44
Wanneer we nu de reeksontwikkeling l~x
e"
=
:,;
gebruiken, krijgen we:
-I {f., R2nF(w)exp(_2j~;L)exp(jw{t+
-f.,
R 2n F(w) exp( -
~ R 2n 2~
-
= 2:: x n
+l {
-I
2~
-I
2j~2nL) exp(jw{ t + c: - ~~}) }dw
F(w) exp(jw{t + ,:,
F(w) exp(jw{t + ,:,
00
" R 2n +1 {F(t + ~ L...J n=O Co
,:,})
-
2~,L}) tM
- 2(n;;I)L})tM }
_ 2nL) _ F(t + ~ _ 2(n+1)L)} C2
Co
C2
00
00
RF(t + ~) Co
+ "L..J R 2n+1 F(t + ~ _ 2nL) _ "L..J R 2n+1 F(t + ~ _ 2(n+1)L)
RF(t + ~)
+ "L...J R 2n+1F(t + ~ _ 2nL) _ "L...J R 2n - 1F(t + ~ _ 2nL)
Co
Co
n=l
C2
Co
n=O
00
C2
00
Co
n=l
C2
Co
n=l
C2
00
RF(t + ~) Co
+ "L...J _R2n - 1(1 - R 2)F(t + ~ _ 2nL)
RF(t + c:)
+ L _R 2n - 1(1 - R)(l + R)F(t + ~ - 2;2L )
Co
n=l
C2
00
n=l
00
R 1->2 F (t + cZo )
+ "T1->2(R1->2)2n-1T2->lF(t +~ _ 2nL) L...J ~ C2 n=l
Hierin zijn:
R 1->2 = R R 2->1 =-R Tk->l = 1 + Rk->l Tk->l is de transmissiecoefficient van een golf die vanuit medium k loodrecht invalt op medium t. Voor de linkerrand van de plak, z
= a geldt nu:
00
£T(O, t) = R 1->2F (t)
+ LT1->2(R2->1)2n- 1T2->lF(t - 2;2L ) n=l
Normeren met
t = 7:": 00
tT(O, f) = R 1->2F (f)
+ L Tl->2(R2->d2n-lT2->lF(t - 2n..jE2 fJ-2.) r
n=l 45
Voor het doorgelaten veld geldt:
it
= T13F(w)uyexp(-~) = Etuy
Hierin is T 13 de transmissiecoefficient [3] van medium 1 naar 3. We kunnen nu net zoals bij het gerefiecteerde veld het doorgelaten veld berekenen. Voor het genormeerde doorgelaten veld op de rechterrand van de plak, Z = 1 voIgt: 00
tt(l, f) =
L
T1->2(R2->1?nT2->lF(f - (2n
+ 1)Jc2 JL2J r
n=O
Met de gevonden formules voor het doorgelaten en gerefiecteerde veld kunnen we het exacte verstrooide veld berekenen.
46
Bibliografie [1] Tijhuis, A.G. Elektromagnetic inverse profiling: theory and numerical implementation. VNU Science Press, Utrecht, Nederland, 1987, pp. 211-222. [2] Scharten, T. Collegedictaat: Elektromagnetisme 3 T.U. Eindhoven, Eindhoven, Nederland, 1990, p. 3.28. [3] Scharten, T. Collegedictaat: Elektromagnetisme 3 T.U. Eindhoven, Eindhoven, Nederland, 1990, p. 3.29.
47
