te maken. Zo kan er bijvoorbeeld gekozen worden voor een hogere orde schema of een aangepast schema om de accuraatheid te vergroten, waarvan enkele in het boek van Cohen gepresenteerd worden.[4] Voor deze scriptie hebben we echter genoeg aan dit schema aangezien deze niet complex is, maar toch goed te gebruiken is in het alternatieve schema dat gepresenteerd zal worden.
4.2
k-Space
We hebben net gekeken naar een methode die een numerieke benadering geeft voor de golf- vergelijking. Er zal nu een alternatief numeriek schema gepresenteerd worden dat k-space genoemd wordt omdat de plaatsafgeleides zal worden benaderd met behulp van Fouriertransfor- maties, waar het frequentiedomein (Fourierdomein/k-domein) een rol speelt. De tijdsafgeleide zal wederom met behulp van tweede orde eindige differenties benaderd worden. Ondanks de fout van de eindige differenties blijkt deze methode een stuk nauwkeuriger te zijn dan de vorige om twee redenen. Behalve dat de plaatsafgeleide met Fouriertransformaties exact opgelost kan worden, zal er bovendien nog een term geïntroduceerd worden in het schema dat de fout van de eindige differenties verkleint. We zullen zien dat deze methode zelfs exact is wanneer de geluidssnelheid in het medium constant is. Onder de aanname dat de eigenschappen van het medium maximaal 5% afwijken van die van water, heeft deze methode een grote accuraatheid voor relatief weinig rasterpunten. Wanneer er naar media gekeken wordt waar de geluidssnelheid grote sprongen maakt, zal eerder voor hogere orde eindige differenties bij de tijdsafgeleide gekozen worden.[1]
Ook bij het stelsel van eerste orde vergelijkingen kan k-space toegepast worden; de plaatsaf- geleides kunnen exact bepaald worden met Fouriertransformaties en dezelfde soort correctieterm kan in het schema geïntroduceerd.[2]Het voordeel van numerieke schema’s gebaseerd op het stelsel eerste orde vergelijkingen is dat er gebruik gemaakt kan worden van staggered grids (zie Appendix), wat de nauwkeurigheid van het schema nog meer vergroot. Tevens is voor het eerste orde stelsel gemakkelijk een Perfectly Matched Layer (PML) te implementeren, die er voor zorgt dat de golven op de randen van het computationele domein geabsorbeerd worden (en er dus geen reflecties ontstaan), wat in de praktijk van belang zal zijn. Bij de tweede orde vergelijking k-space methode is een dergelijke PML niet gemakkelijk in het schema te implementeren. Een reden om de k-space methode voor de tweede orde vergelijking te bestuderen is dat het concept minder complex en computationeel efficiënter is; het akoestische veld wordt namelijk alleen gedefiniëerd door de akoestische druk, in plaats van de akoestische druk en de snelheid van de deeltjes.[2]
4.2.1 Afleiding k-space methode
In dit deel zal bekeken worden hoe de benaderingen uitgevoerd worden en waar de correctieterm vandaan komt. Op deze manier kunnen we het tijdsstapschema van k-space afleiden. We zullen de afleiding maken vanuit de golfvergelijking met een constante c0. De golfvergelijking is dan
4.2. k-SPACE van de vorm: (4.3) ∂ 2u ∂t2(x, t) = c 2 0 ∂2u ∂x2(x, t).
Deze vergelijking zullen we weer numeriek benaderen op een raster met dezelfde onderlinge afstanden∆t en ∆x. Om te beginnen zal de tijdsafgeleide met eindige differentie te benaderd worden. Zoals we net zagen kan de tweede afgeleide benaderd worden door
∂2u
∂t2(x, t) ≈
u(x, t +∆t)−2u(x, t)+ u(x, t −∆t) (∆t)2 .
De plaatsafgeleide zal op een totaal andere manier benaderd worden. Zoals in het hoofdstuk Fouriertransformaties getoond is, geldt F £∂2u
∂x2(x, t)¤ = −k2u(k, t) = −kˆ 2F [u(x, t)]. Hieruit volgt
dat dat het rechterlid van (4.3) te schrijven is als
∂2u
∂x2(x, t) = −F
−1£k2F [u(x, t)]¤.
De numerieke benadering op het raster voor deze Fouriertransformaties kan op zeer accuraat uitgevoerd worden door gebruik te maken van de (Inverse) Fast Fourier Transform.[1] Het is een techniek die in 1965 door Cooley en Tukey gepubliceerd is.[21]De fout van de FFT bij een zekere afstand ∆x is haast verwaarloosbaar vergeleken met bijvoorbeeld de fout van tweede eindige differenties waar dezelfde afstand gebruikt zou worden. Omdat voorwaartse en inverse transformatie ingebouwde functies inMATLABzijn, zal ik me hier niet verder aan wijden. Op deze manier komt al een verbeterde numerieke benadering voor de golfvergelijking tot stand:
u(x, t +∆t)−2u(x, t)+ u(x, t −∆t) = −c20F−1£k2F [u(x, t)]¤.
Nu is er alleen nog het probleem dat de benadering voor de tijdsafgeleide een relatief grote fout met zich meebrengt. We zagen net namelijk al dat de benadering met eindige differenties een fout ter grootte vanO¡(∆t)2¢ met zich meebrengt. Zoals al eerder is vermeld, is het mogelijk naar eigenschappen te kijken van oplossingen door slechts te kijken naar de oplossing van de vorm ei(kx+ωt). Wat gebeurt er met deze oplossing als we hierop de tweede orde centrale eindige differentie toepassen?
ei(kx+ω(t+∆t))− 2ei(kx+ωt)+ ei(kx+ω(t−∆t)) (∆t)2 = e i(kx+ωt)eiω∆t− 2 + e−iω∆t (∆t)2 = ei(kx+ωt)2 cos(ω∆t)−2 (∆t)2 = −ω2ei(kx+ωt)sinc2(ω∆t 2 ).
Hier wordt bij de laatste gelijkheid eerst gebruik gemaakt van 2 − 2cos(ω∆t) = 4sin2(ω∆t2 ) en vervolgens de term sinc(x) =sin(x)x geintroduceerd. We zien dat de afgeleide niet exact is, maar er een extra term bij de afgeleide komt: sinc2(ω∆t2 ). Als we deze wegdelen, krijgen we een dus
4.2. k-SPACE
gecorrigeerde eindige differenties waar de fout opgeheven is en dus een exacte waarde van de tweede afgeleide geeft.
Als we bovenstaande toepassen op de golfvergelijking, kan er een exacte benadering van de golfvergelijking gegeven worden in de volgende termen:
u(x, t +∆t)−2u(x, t)+ u(x, t −∆t) (∆t)2sinc2(ω∆t
2 )
= −c20F−1£k2F [u(x, t)]¤.
Met behulp van de dispersierelatieω= ck kan bovenstaande vergelijking omgeschreven worden. Hieruit volgt de volgende vergelijking, waar de gehele operatie die op u wordt toegepast aan de rechterkant de tweede orde k-space operator genoemd wordt:
(4.4) u(x, t +∆t)−2u(x, t)+ u(x, t −∆t) = −c2F−1£k2(∆t)2sinc2(crefk∆t
2 )F [u(x, t)]¤.
Voor de afleiding werd een constante c0gebruikt, maar die is nu vervangen door een c die variabel
kan zijn. Wat waarschijnlijk ook opvalt is dat er nu cref in de correctieterm staat. Hier moet
namelijk een juiste referentiesnelheid gekozen worden, die zo veel mogelijk in de buurt van c ligt, zodat de correctieterm daadwerkelijk de fout verkleint. Voor c0 constant wordt uiteraard
cref= c0 gekozen. Omdat crefconstant is kan deze niet gelijk aan c genomen worden als deze
laatste variabel is; dan moet er een passende waarde voor crefgekozen worden.
Bovenstaande vorm kan tot slot omgeschreven worden naar een tijdsstapschema, waarmee het drukveld op een volgend tijdstip bepaald kan worden aan de hand van het drukveld op voorgaande tijdstippen. Voor het (x, t)-raster zoals eerder gedefiniëerd kan namelijk de verzameling van drukwaardes um+1= {um1, u2m, . . . , umn}op het tijdstip tm+1bepaald worden door:
(4.5) um+1= 2um− um−1− c20F−1£k2(∆t)2sinc2(crefk∆t
2 )F [u
m]¤.
Let hier wel op het verschil van de drukwaardes tussen dit schema en het schema met eindige differenties (4.2). Bij eindige differenties kan de drukwaarde op het punt (xn, tm+1) bepaald
worden met behulp van waardes op enkele omliggende punten. Bij k-space wordt echter niet per punt, maar gelijk de gehele verzameling drukwaardes op het volgende tijdstip berekend vanuit de twee verzamelingen van drukwaardes op voorgaande tijdstippen.
tm−1 tm
tm+1
x1 x2 x3 . . . xn−1 xn
4.3. BEGINCONDITIES IMPLEMENTEREN
4.2.2 Vervolg k-space methode
Als we nu een oplossing invullen voor de golfvergelijking met constante snelheid, klopt ver- gelijking (4.4) dan? Om onszelf te overtuigen hiervan kunnen we dit direct nagaan door naar de oplossing u(x, t) = F(x − ct) + G(x + ct) te kijken, wat volgens (2.5) een algemene oplossing is. Om onderscheid te maken tussen Fouriertransformaties van tijd en plaats schrijven we respectievelijkFt en Fx. Tevens kunnen we voor het gemak de volgende notatie gebruiken: Fx[u(x, t)] = Fxu(k, t). Met deze schrijfwijze komt er bij het rechterlid het volgende:
F [u(x, t)] =p1 2π
Z ∞
−∞
u(x, t)e−ikx dx =p1 2π Z ∞ −∞F(x − c0 t)e−ikx dx +p1 2π Z ∞ −∞G(x + c0 t)e−ikx dx =p1 2π Z ∞ −∞ F(z)e−ik(z+c0t)d z + 1 p 2π Z ∞ −∞ G(z)e−ik(z−c0t) d z = e−ikc0t 1 p 2π Z ∞ −∞
F(z)e−ikz d z + eikc0t 1
p 2π Z ∞ −∞ G(z)e−ikz d z = e−ikc0tF xF(k) + eikc0tFxG(k).
Na multiplicatie met onder andere de correctieterm moet er vervolgens weer terug getrans- formeerd worden. Als de equivalentie gebruikt wordt dat
(4.6) k2(∆t)2sinc2(c0k∆t 2 ) = 4 c20sin 2(c0k∆t 2 ) = 1 c20(2 − 2cos(c0k∆t)) = 1 c20(2 − e ic0k∆t − e−ic0k∆t),
kan met cref= c0 de inverse transformatie als volgt uitgeschreven worden:
F−1hk2(∆t)2sinc2(c0k∆t 2 ) ³ e−ikc0tF xF(k) + eikc0tFxG(k) ´i =−1 c20 ³
u(x, t +∆t)−2u(x, t)+ u(x, t −∆t)´.
De uitgebreide uitwerking van deze laatste gelijkheid staat in A2. Duidelijk is nu dat het rechterlid en linkerlid van (4.4) overeenkomt; de term met c0 valt weg tegen zichzelf in het
rechterlid die we nog niet meegenomen hebben in de berekening.
4.3
Begincondities implementeren
Er zijn twee methodes geïntroduceerd waarmee de drukwaardes op een volgend tijdstip benaderd kunnen worden wanneer de drukwaardes op de voorgaande twee tijdstippen bekend zijn. Hoe zit het dan op het begintijdstip t0? Op dit tijdstip is er geen drukwaardes bekend op het voorgaande
tijdstip, maar er zijn wel begincondities gegeven. Voor de eerste stap van t0 naar t1zullen de formules aangepast worden zodat de formules uitgedrukt worden in de begincondities. Stel dat
4.3. BEGINCONDITIES IMPLEMENTEREN
de beginwaardes als volgt gegeven zijn:
u(x, 0) = f (x),
∂
∂tu(x, 0) = g(x).
Voor de eerste stap op t0is de waarde u0ngegeven door f te discretiseren. Dat wil zeggen dat we
f bekijken op alle rasterpunt {xn}, zodat volgt:
fn= f (xn) = u0n.
Net zo zal g gediscretiseerd worden door alleen naar de waardes op de rasterpunten {xn}te kijken.
Met behulp van eindige differenties voor eerste afgeleides kunnen we de waarde van g op elk rasterpunt (xn) schrijven als
gn= g(xn) = ∂
∂tu(xn, 0) ≈
u1n− u−1n 2∆t .
Op deze manier kunnen de fictieve waardes u−1n gesubstitueerd worden met behulp van de volgende vergelijking, die wegens de benadering wel een kleine fout met zich meebrengt:
u−1n = u1n− 2∆tgn.
Voor de eerste tijdsstap kunnen deze waardes dus in de schema’s geïmplementeerd worden. Het schema van eindige differentie zal dan voor elk punt (xn, t1) de drukwaarde kunnen berekenen
met de formule:
(4.7) u1n=s
2( fn−1+ fn+1) + (1 − s)fn+∆tgn.
Het schema van k-space zal met de volgende formule de hele verzameling drukwaardes berekenen op het tijdstip t1: (4.8) u1= f +∆tg −1 2c 2¡ F−1£k2(∆t)2sinc2¡ crefk∆t 2 ¢ F [f ]¤¢.
5
ANALYSE VAN DE METHODES
Er zijn twee numerieke tijdsstapschema’s geïntroduceerd waarmee golven gesimuleerd kunnen worden. Er is echter al opgemerkt dat eindige differenties in een homogeen medium al een fout geven. Er kan een Gaussiaanse beginconditie u = e−x2 en ∂t∂u = 0 bekeken worden om het verschil tussen eindige differenties en k-space globaal te kunnen aanschouwen. Deze kan met beide methodes voor een bepaald aantal tijdstappen gesimuleerd worden en vergeleken worden met d’Alemberts oplossing. Voor bijvoorbeeld∆x = 0,5 kunnen we uit deze figuren opmaken dat de methode van eindige differenties al zeer grote fouten oplevert, terwijl de k-space methode gelijk loopt met de exacte oplossing. Eindige differenties zal uiteraard nauwkeuriger worden
Eindige differenties k-space methode d’Alembert oplossing
Figuur 5.1: Voorbeeldsimulatie van een Gaussiaanse puls in het (u, x)-vlak na een bepaalde tijd op een raster met∆x = 0,5.
naarmate∆x kleiner wordt, maar daardoor wordt de rekentijd ook langer. Bovendien zien we al uit Figuur 5.1 dat k-space duidelijk nauwkeuriger is op hetzelfde raster. Om deze reden zal eindige differenties aan de kant geschoven worden en k-space beter bestudeerd worden.
Bij k-space moet er een referentiesnelheid gekozen worden. Wanneer deze gelijk aan de constante geluidssnelheid in het medium gekozen wordt, zal deze methode gelijk lopen met de exacte golf (voor een klein genoeg raster, zie volgend hoofdstuk). Maar wat gebeurt er wanneer de referentiesnelheid ongelijk aan de geluidssnelheid in het medium is? Als voorbeeld zijn de exacte golf en de k-space gesimuleerde golf na een bepaalde tijd in één figuur gezet voor∆x = 0,5 en is de golf die naar links propageert uitvergroot. In Figuur 5.2a zijn de twee golven niet van
5.1. HARMONISCHE GOLVEN
elkaar te onderscheiden. Hier is dan ook de referentiesnelheid gelijk aan de geluidssnelheid van het medium gekozen. In Figuur 5.2b is de golf wederom met k-space gesimuleerd, maar nu is de referentiesnelheid cref= 1000 m/s, terwijl de geluidssnelheid in het medium 1500 m/s is.
(a) Uitvergroting k-space met cref= c0 en
exacte oplossing. De golven lopen precies gelijk.
(b) Uitvergroting k-space met cref6= c0 en
exacte oplossing. Nu geeft k-space duidelijk een fout.
Figuur 5.2: Simulatie met k-space en de exacte oplossing in dezelfde figuren op een raster met ∆x = 0,5. Het is dus van belang om de referentiesnelheid goed te kiezen om (grote) fouten te voorkomen.
In de praktijk zal het medium waarin de golf zich voorplant niet homogeen zijn, omdat er naar biologisch weefsel gekeken zal worden. Er zullen dan delen zijn in het medium waar de geluidssnelheid niet overeenkomt met de gekozen referentiesnelheid. De correctieterm zal de fout dan niet geheel opheffen. We kunnen in kaart brengen hoeveel k-space zal verschillen van de werkelijke golf, door te kijken naar de fasefout (phase error). Bij dit type fout kijk je naar het faseverschil van de twee golven op een bepaald tijdstip. Om te bekijken hoe de nauwkeurigheid in termen van de fase van k-space verandert, kan er gekeken worden naar de exacte oplossing van een golf die met een vaste snelheid van 1500 m/s propageert en cref laten variëren. Een
handige manier om het verschil in fase te bekijken is door de golven op het frequentiedomein te bekijken. We zullen ons verdiepen in een resultaat uit het artikel van Treeby et al.[1]over de fout die zal ontstaan bij k-space wanneer crefongelijk aan c0is. In het volgende hoofdstuk zal de
rasterafstand∆x aan bod komen.
5.1
Harmonische golven
We zullen de analyse maken met behulp van de oplossing u(x, t) = ei(kx−ωt)+ ei(kx+ωt) die met een snelheid van c0= 1500 m/s propageren. Er kan dan gekeken worden naar het verschil in
5.1. HARMONISCHE GOLVEN
fase tussen de numerieke simulatie en de exacte oplossing na bijvoorbeeld 50 golflengtes. Als we golflengteλnemen, dan kunnen er enkele andere eigenschappen van de golf berekend worden:
1. De lengte van het interval waar exact 50 golflengtes in zitten: L = 50 ·λ. 2. De periode: T =1500λ .
3. Het frequentiegetal: kl=2Lπ· 50.
4. Het golfgetal:ωl= 1500 · kl
Voor de k-space methode is er tevens een geschikt raster nodig. Een verhouding tussen de ruimte- en tijdsafstand wordt gegeven door het zogenaamdeCFL-getal. Dit getal bepaalt op de volgende wijze de verhouding:
(5.1) CFL= c20∆t
∆x
Met de waarde CFL= 0, 3 wordt een geschikte verhouding tussen de afstanden verkregen; dit
geeft een goede balans tussen nauwkeurigheid en rekensnelheid in een zwak heterogeen medium. Verder zal er gebruik gemaakt worden van 4 rasterpunten per golflengte in de analyse van de fasefout.[1]
We kunnen de representatie van de harmonische golf op het frequentiedomein bekijken door de Fouriercoëfficiënten van de Fourierreeks te berekenen van u. Het deel van de oplossing met de harmonische golf die naar links beweegt wordt in het Fourierdomein gerepresenteerd wordt door:
cj(t) = 1 L Z L/2 −L/2 ei(klx+ωlt)e−i j2πx/L dx =1 Le iωlt Z L/2 −L/2 ei(kl− j2π/L)xdx.
De laatste integraal in de vergelijking is wegens orthogonaliteit gelijk aan 0 als kl6= j2Lπ en gelijk aan 1 als kl= j2Lπ; ditzelfde geldt voor de andere harmonische golf. Oftewel, de enige coëfficient
ˆ
ujdie geen niet-nul oplossing is, is voor j = 50 en wordt gegeven door
ˆ uj(t) = 1 Le iωlt +1 Le −iωlt.
We zullen voor het analyseren van het faseverschil tussen de originele oplossing en de numerieke oplossing kijken naar één richting van de golf, dus ei(kx+ωlt). Deze wordt in het
Fourierdomein gerepresenteerd door 1Leiωlt. De golf beweegt naar links en daarbij zal duidelijk
ook de oplossing in het Fourierdomein veranderd worden (want t verandert). Echter, na een bepaalde tijd zal de golf er weer hetzelfde uitzien omdat deze periodiek is; om precies te zijn na één periode T =2ωπ. De oplossing ziet er dan in het Fourierdomein ook weer hetzelfde uit, want
eiωl(t+T)
= ei(ωlt+2π)
5.2. NUMERIEKE OPLOSSING IN HOMOGEEN MEDIUM
Er is dus een één op één relatie tussen het aantal periodes dat afgelegd is in het ruimtelijke domein en hoeveel keer de eenheidscirkel op het complexe vlak is afgelegd door de representatie in het Fourierdomein. Er zal nu gekeken worden naar de originele golf en een door k-space numeriek gesimuleerde golf. Wat we zullen merken is dat de amplitude van de golven kunnen gaan verschillen (amplitudefout) en dat een golf verder is dan de andere (fasefout). In deze tweede type fout zijn we geïnteresseerd en kan bestudeerd worden aan de hand van de representaties in het complexe Fourierdomein.
5.2
Numerieke oplossing in homogeen medium
Als bij (4.4) aan beide kanten de Fouriertransformatie naar x wordt uitgevoerd, wordt de volgende differentievergelijking verkregen:
Fxu(k, t −∆t)−2Fxu(k,∆t)+Fxu(k, t +∆t) = −c2k2(∆t)2sinc2(
crefk∆t
2 )Fxu(k, t).
Voor k = klen cref= c = c0zal ˆujaan bovenstaande vergelijking voldoen. We kunnen daarentegen
wat voor oplossing er wordt verkregen als we de differentievergelijking oplossen. Deze oplossing zou voor k = klen cref= c = c0exact gelijk moeten zijn aan ˆuj; in het bijzonder willen we dat de
fase dan gelijk is.
Als crefgelijk gekozen wordt aan c = c0 zal hieruit de exacte oplossing moeten komen (want k-space lost exact op voor constante c). Voor het gemak maken we gebruik van de notatie op een (k, t)-raster; op het rasterpunt (kl, tm) nemen we ˆuml := ˆu(kl, tm), waar kl nog steeds het
golfgetal van u voorstelt. Met deze notatie en omschrijven van bovenstaande vergelijking, komt de volgende vergelijking tot stand:
(5.2) Fxum−1l − AFxuml + Fxum+1l = 0,
waar A = 2 − c2k2l(∆t)2sinc2¡crefkl∆t
2 ¢. Een bekende manier om dit op te lossen is door aan te
nemen dat de oplossing van de vormFxuml = a1· bm+ a2· bmis, waar m nu een macht is. Deze
vorm invullen in (5.2) geeft voor m = 1 de volgende gelijkheid:
1 − A · b + b2= 0.
Om te controleren dat deze oplossing van de differentievergelijking exact gelijk is aan ˆuj(t) als
cref= c constant is, zullen we het voor deze constante snelheid oplossen. Ten eerste valt met
behulp van de gelijkheden (4.6) A gemakkelijker te schrijven:
A = 2 − c20k 2 l(∆t) 2sinc2³c0kl∆t 2 ´ = 2 cos(c0kl∆t)
5.3. NUMERIEKE OPLOSSING IN HETEROGEEN MEDIUM
Als we met deze A de differentievergelijking (5.2) oplossen met behulp van de ABC-formule, volgen de twee oplossingen:
b =2 cos(c0kl∆t) 2 ± p 4 cos2(c 0kl∆t)−4 2 = cos(c0kl∆t)± q −4 sin2(c0kl∆t) 2 = cos(c0kl∆t)± i sin(c0kl∆t) = e±ic0kl∆t
De algemene oplossing voor de differentievergelijking wordt dus gegeven door:
a1· e(ic0kl∆t)m+ a2· e(−ic0kl∆t)m
Het nulpunt met het plus-teken komt exact overeen met de harmonische golf met plus-teken (en het nulpunt met het min-teken komt overeen met de harmonische oplossing met min-teken). Dit is in te zien door één van de twee nulpunten te kiezen, een bepaald aantal tijdstappen te laten doorlopen en te kijken of op dit eindpunt de numerieke oplossing dezelfde fase heeft als de overeenkomende harmonische golf. De tijdsduur D om voor de harmonische golf 50 golflengtes te overbruggen is gelijk aan 50 periodes. Dus de fase van de golf na deze tijd kan bekeken worden door te kijken naar eiωD = ei1500·50·2Lπ·50T = ei100π. Voor deze tijdsduur zijn er precies n = D
(∆t)
aantal tijdsstappen nodig. De fase die de oplossing uit de differentievergelijking heeft na deze tijdsduur is is af te leiden uit e(ic0k∆t)n= ei100π. De golven zijn na tijd D dus nog steeds in fase;
voor slechts vier golfpunten per golflengte is dat een indrukwekkend resultaat.
5.3
Numerieke oplossing in heterogeen medium
Wanneer de snelheid in het medium niet constant is, zal er echter wel een fasefout op gaan treden. We kunnen namelijk net als hiervoor via de differentievergelijking een oplossing vinden, maar nu nemen we dat cref niet per se gelijk hoeft te zijn aan c0. Om te beginnen zal de differentievergelijking (5.2) bekeken worden waar nu A gegeven wordt door:
A = 2 − c20k2l(∆t)
2sinc2³crefkl∆t
2 ´
De wortels die met deze A verkregen worden uit de differentievergelijking, worden in termen van A als volgt gegeven:
b = A 2 ±
p A2− 4
2
Omdat c06= crefkan deze wortel niet op een kortere manier geschreven worden. We kunnen deze
wortel echter wel analyseren door te kijken naar een bepaalde macht ervan. Als c0= 1500, kan bnbepaald worden voor verschillende waardes van cref. In het bijzonder kan het faseverschil
5.3. NUMERIEKE OPLOSSING IN HETEROGEEN MEDIUM
Figuur 5.3: Uitkomsten van bnvoor de waardes cref= 10, 20, ..., 2000.
Figuur 5.4: Het faseverschil in procenten tussen bnen de exacte oplossing.
tussen de uitkomsten en de originele oplossing bekeken worden. Merk hiervoor op de originele oplossing na 50 golflengtes op ei100π= 1 zit.
Uit deze figuren is op te maken dat er een faseverschil optreedt wanneer de referentiesnelheid ongelijk aan de geluidssnelheid in het medium is. Na propagatie over een afstand van 50 golflengtes kan het faseverschil al oplopen tot 50% als de referentiesnelheid erg verschilt van c0. In de praktijk zullen ultrasone golven met golflengtes van 1 − 20 MHz gebruikt worden.[15]
Om een indruk te geven van de hoeveelheid golflengtes die in een medisch proces vanaf het weefsel zal worden afgelegd tot de detectoren, kan er bijvoorbeeld gekeken worden naar golven van rond de 3 MHz. Tot een diepte van 15 cm kunnen deze golven opgevangen worden door de juiste detectoren.[1]Een golf van rond de 3 MHz zal in zacht weefsel over de afstand van 15 cm op zijn fundamentele frequentie al ongeveer 300 golflengtes hebben. Een fout bij één golflengte kan dramatische gevolgen hebben wanneer er naar 300 golflengtes gekeken wordt. Om deze reden is dan ook de aanname gemaakt de k-space in zwak heterogene media gebruikt wordt, waarvan de