Akoestische modellering van axisymmetrische openingen met Rayleigh elementen J.P. Starreveld 420726
Stage-Rapport 97.044 WFW-dynamica Eindhoven, Juli 1997
Begeleider: Afstudeerdocent:
ir. A.H.W.M. Kuijpers prof.dr.ir. D.H. van Campen
Technische Universiteit Eindhoven Faculteit Werktuigbouwkunde
Samenvatting Dit verslag behandelt de ontwikkeling van een Boundary Elementen Methode (BEM) formulering voor de akoestische modellering van een cilinder met aan de uiteinden vlakke oneindige platen (baffles). Met behulp van deze BEM formulering in combinatie met het BEM pakket BARD is het mogelijk om de geluidsproductie t e berekenen voor zowel binnen als buiten een gebaffelde cilinder veroorzaakt door niet-axisymmetrische trillingen. Dit model kan mogelijk uitkomst bieden voor het beschrijven van de akoestische afstraling van een MRI-systeem. De BEM formulering is afgeleid van de zogenaamde ‘Rayleigh-integraal’ en is gebaseerd op een theorie voor de beschrijving van de akoestische afstraling voor axisymmetrisch constructies veroorzaakt door niet-axisymmetrische trillingen. Voor het bepalen van de ‘Rayleighintegralen’ wordt gebruik gemaakt van drie verschillende integratie methoden, zijnde een hogere orde, Gauss-Legendre- en Log-Gauss-integratie. Log-Gauss-integratie levert de beste resultaten bij integratie over de ‘Rayleigh-elementen’ in de singuliere knooppunten. De koppeling is verkregen door compatibiliteit t e eisen op de contactvlakken van beide BEM sommen. Bij gebruik van de standaard instellingen kunnen goede tot zeer goede resultaten worden bereikt. Verhoging van het aantal integratie-punten en het kiezen van een kleinere toleîaritie leidt tot een verbeteringen van de resultaten. Daar staat tegenover dat d e rekentijd hierdoor toeneemt. Het toepassen van meshverfijning is noodzakelijk, wanneer er minder dan 6 knooppunten per golflengte zijn. Het is verstandig om de hier ontwikkelde BEM formulering op basis van ‘Rayleigh-elementen’ binnen MATLAB te vertalen naar C++. Dit zal leiden tot afname van de rekentijden, een betere efficiëntie en het kunnen gebruiken van slechts één invoer-file. Daarnaast moet worden onderzocht of het gebaffelde cilinder model gebruikt kan worden voor de modellering van de akoestische afstraling van een MRI-systeem.
Inhoudsopgave Samenvatting
1
1 Inleiding
1
.................................. Opdrachtomschrijving . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Gevolgde strategie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.1 Probleemstelling
1
i .2
1
1.3
2 Theoretische afleiding van de geluidsafstraling van een baffle
2.1 Rayleigh integraal
.................................
2.2 Rayleigh integraal voor axisymmetrische constructies met niet-axisymmetrische randvoorwaarden . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
...... Berekening van de Rayleigh integraal over het contactvlak . . . . . . . Berekenen van de integraal H , . . . . . . . . . . . . . . . . . . . . . .
2.3 Constructie van de benaderingsoplossing voor de Rayleigh integraal 2.3.1 2.3.2
3 Koppeling van Rayleigh integraal aan BEM
3.1 Theoretische afleiding ten behoeve van koppeling
1 3
3
5 7
7
8 11
................
3.2 Koppelen van axisymmetrische constructies met niet-axisymmetrische randvoorwaarden . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Numerieke implementatie
................................ Numerieke integratie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.1 Berekenen van de integraal H , . . . . . . . . . . . . . . . . . . . . . . 4.2.2 Berekening van de Rayleigh integraal over het contactvlak . . . . . . . Koppeling van Rayleigh integraal aan BEM . . . . . . . . . . . . . . . . . . .
12 12 15
4.1 Programma opbouw
15
4.2
17
4.3
17 18 20
iv
5 Verificatie met analytische oplossingen 5.1
................
23
... 5.1.2 Geluidsdruk in het verre veld voor een zeer locale wandtrilling . . . . 5.1.3 Geluidsdruk op de as van symmetrie . . . . . . . . . . . . . . . . . . . 5.1.4 De ontwikkeling van ‘druklobben’ bij toenemende axiale afstand . . . Geluidsafstraling van een pijp met een opening aan één zijde . . . . . . . . . Pulserende cilinderwand in een akoestisch baffle . . . . . . . . . . . . . . . . .
23
Oscillerende zuiger in een vlakke oneindige plaat 5.1.1
5.2 5.3 6
23
Geluidsdruk op het wandoppervlak voor laag frequente oscillaties
Conclusies en aanbevelingen 6.1 6.2
...................................... Aanbevelingen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Conclusies
25 26 26 28 29
31 31 31
A Opbouw van het BEM-pakket (MATLAB-procedures)
35
B Numerieke integratie
37
....................... Formule van Simpson . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . QUAD8 .standaard MATLAB-procedure . . . . . . . . . . . . . . . . . . . . GUASSHM .BARD procedure . . . . . . . . . . . . . . . . . . . . . . . . . .
B.l Formule van Euler (trapeziumregel)
37
B.2
37
B.3 B.4
C Waarden van f en
wint
ten behoeve van Gauss-Legendre-integratie
D Log-Gauss-transformatie D.1 G e v a l 1 D.2 Geval2 D.3 Geval3
....................................... ....................................... .......................................
E Waarden van 7 en
wint
ten behoeve van Log-Gauss-integratie
38 38
39 41
41 42
43 45
Hoofdstuk i
Inleiding 1.P
Probleemstelling
Het project ‘Het geluidsbewust ontwerpen van MRI-systemen’ heeft als doelstelling om de mechanica en akoestiek van ‘Magnetic Resonance Imaging’ (MR1)-scanners beter te begrijpen en te verbeteren. De in dit project opgedane kennis moet worden gebruikt om al in de ontwerpfase van deze scanner gericht te kunnen zoeken naar een zo stil mogelijk MRI-systeem. In het kader van het bovengenoemde project is in eigen beheer het Boundary Elementen Methode (BEM) pakket BARD geschreven, waarmee de geluidsproduktie kan worden berekend van axisymmetrische constructies veroorzaakt door (niet-axisymmetrische) trillingen. Voor de modellering van de akoestische afstraling van een MRI-systeem wordt in eerste instantie uitgegaan van een zogenaamd gebaffelde cilinder model - een model van de afstraling van een cilinder met aan de uiteinden vlakke oneindige platen (baffles). De beschrijving van die baffles is echter nog niet in het programma aanwezig. Een adequate beschrijving van de akoestische baffle aan de uiteinden van een cilinder moet daarom worden ontwikkeld.
1.2
Opdracht omschrijving
Het doel van de stageopdracht is de ontwikkeling, implementatie en het testen van een (numerieke) formulering waarmee de afstraling vanuit de uiteinden van een axisymmetrische constructie in een akoestische baffle kan worden beschreven. Deze formulering moet gekoppeld worden aan de BEM die reeds geïmplementeerd is voor de berekening van de geluidsafstraling in de cilinder.
1.3
Gevolgde strategie
In hoofdstuk 2 zal een theoretische afleiding volgen van de formulering waarmee de geluidsafstraling vanuit het uiteinde van een axisymmetrische constructie in een akoestische baffle kan worden beschreven. Daarbij is gebruik gemaakt van theorie uit [i] en [2]. Deze afleiding
2
Inleiding
zal leiden tot de zogenaamde ‘Rayleigh-integraal’. Vervolgens zal overgegaan worden t o t de constructie van de benaderingsoplossing ten behoeve van de numerieke implementatie. [i], [3] en [4] is de hierbij gebruikte theorie. Hoofdstuk 3 behandelt de koppeling van de ‘Rayleigh-integraal’ aan de BEM. In [3] en [5] wordt eeri methode beschreven om twee afzonderlijke BEM sommen op de contactvlakken t e koppelen. Deze Koppeling wordt verkiegeri door compatibliteit t e eisen op de cûntacivIakken. De implementatie van bovengenoemde formulering is uitgevoerd binnen MATLAB. In hoofdstuk 4 wordt de opbouw van het programma beschreven. Daarnaast zullen enkele bijzonderheden binnen het programma aan het licht worden gebracht, zoals de benodigde integratie methoden voor het bepalen van de ‘Rayleigh-integraal’. Ook zal worden ingegaan hoe de numerieke koppeling tussen de twee afzonderlijke BEM sommen tot stand komt. In hoofdstuk 5 worden enkele analytische oplossingen van de ‘Rayleigh-integraal’ gepresenteerd. Ook wordt de oplossing van een gekoppeld probleem gepresenteerd. Deze oplossingen worden vergeleken met de hier ontwikkelde methode. Hoofdstuk 6 beschrijft enkele conclusies en aanbevelingen.
Hoofdstuk 2
Theoretische afleiding van de geluidsafstraling van een baffle In dit ho’ofdstuk zal de formulering worden afgeleid waarmee de geluidsafstraling vanuit het uiteinde van een axisymmetrische constructie in een akoestisch baffle kan worden beschreven.
2. I
Rayleigh int egraal
De algemene definitie van geluid is: ‘De mechanische trillingen die zich als golven uitbreiden in elastische media in het gebied van de hoorbare frequenties (16 [Hz] - 16 [kHz])’. De voortplanting van deze harmonische akoestische golven kunnen worden beschreven met de Helmholtz-vergelijking:
waarin 4 de scalaire snelheidspotentiaal die evenredig is met de amplitude van de harmonische drukverstoring (8 = ikpcq5) en k het golfgetal welke gedefinieerd is als:
IC=-
W C f
waarin w en c respectievelijk de golffrequentie en de geluidsnelheid zijn. Met het theorema volgens Kirchhoff en Helmholtz [2], kan dit ook geschreven worden als:
Q(V2
+ k 2 ) 4 - q5(V2+ k2)P= v
(QVq5 - 4VP),
waarin P een functie van de positie is. Door integratie over het volume V en toepassen van het Gauss-theorema, en vanwege het feit dat (V2 k 2 ) 4 = O binnen V , volgt:
+
waar IR en n S respectievelijk de oppervlakte integraal over het buiten oppervlak en de eenheidsnormaal op dit oppervlak zijn. Door Q te definiëren als zijnde de Green’s functie
Theoretische afleiding van de g.eluidsafstraling: van een baffle
4
Q(P,Q) = e-ikR(p7Q)/R(P,Q) ( R = IP - QI) - een functie die de oplossing is voor een puntbron in een 3 dimensionaal oneindig systeem - wordt het linkerlid gelijk aan C(P)q5(P) als gegeven is dat P binnen V ligt. Daarnaast volgt met de Sommerfeld stralingsconditie - ik$)] = O) dat IR gelijk nul is: (limr+.m [r
(2
(2.1’I
Dit vormt de basis voor de BEM. Daarbij is C ( P ) = 47r wanneer P in het vrije veld ligt en C ( P )= 27r als P op het oppervlak ligt. De algemene definitie van C ( P ) luidt: C ( P )=
47r + sss
& (&)
WQ).
De modellering van de gebaffelde cilinder wordt als een axisymmetrisch probleem beschouwd en is schematisch weergegeven in figuur 2.1. Hierbij wordt een cilindrisch coördinaten stelsel
Z
Figuur 2.1: Schematische weergave van de gebaffelde cilinder gehanteerd. Ten gevoige van de in axiale richting volledig refiecterende wand kunnen we het zogenaamde spiegeleffect toepassen. Het spiegelpunt van P is in de figuur aangegeven met Pl. Het punt Q ligt op het contactvlak tussen cilinder en baffle. Door het spiegeleffect toe t e passen voor het punt Q kan de Green’s functie geschreven worden als:
met:
+ R(P’,Q ) = JP’- Q J= [rpr + R(P,Q ) = IP - &I = [Y;
T&
- 2rprQ COS@
T&
+ (Zp + ( Z p t - ZQ)’];, ZQ)2]i
- 2TptfQ C O S 8
waarbij : r p ; rp’; TQ
z p ; zpt; ZQ 8 = @Q - @p
radiale afstand van respectievelijk P , PI en Q tot de rotatie as 3 axiale afstand van respectievelijk P , PI en Q tot de wand en P 5 hoek tussen
2.2 Rayleigh integraal voor axisymmetrische constructies met
niet-axisymmetrische randvoorwaarden
5
Omdat voor een gebaffelde cilinder r p = r p i , z p = -zpi en functie geschreven worden als:
R(P,Q ) = [rg + r i
- 2rprQ cos
ZQ
= O geldt, kun de Green's
e + ,413.
De factor 2 in Q(P,&) kan fysisch worden verklaard door het wandoppervlak te beschouwen als een infinitesimale dunne permeabele plaat. De normaalsnelheid is dan aan beide azQ zijde van de plaat van gelijke grootte en richting (het medium links van de plaat beweegt naar links en het medium rechts van de plaat beweegt ook naar links). Vandaar d a t Q(P,Q ) voor vermenigvuldigd wordt met de waarde 2 ((2.2)). Vanwege het feit dat Q' # Q'(zQ) kan de afgeleide van de Green's functie geschreven worden als:
Beschouwen we weer de infinitesimale dunne permeabele plaat, dan kan deze formule samen met (2.2) worden verklaard. Wanneer ter plaatse P en P' een gelijke druk w $ ( P ) wordt opgewekt, dan wordt in het punt Q links en rechts van de plaat een gelijke druk waargenomen. De netto bijdrage van de oppervlakte druk in de Kirchhoff-Helmholtz integraal is daarom gelijk aan nul. Vandaar dat de integraal (2.1) nu overgaat in:
Deze vergelijking wordt ook wel de Rayleigh integraal genoemd. Omdat C ( P ) geen functie van ZQ is, geldt altijd dat C ( P )= 47c. Dit volgt uit de algemene definitie van C ( P )(pagina 4).
2.2
Rayleigh integraal voor axisymmetrische constructies met niet-axisymmetrische randvoorwaarden
Het BEM pakket BARD is gebaseerd op een theorie die de geluidsafstraling berekend voor axisymmetrische constructies veroorzaakt door niet-axisymmetrische trillingen. Vandaar dat een dergelijke theorie ook moet worden toegepast voor de afstraling vanuit d e uiteinden van een cilindische constructie in een akoestische baffle. Ten gevolge van de axisymmetrie kunnen de niet-axisymmetrisch variabelen 4 ( P ) ,e - z W
P,Q )
R(P,Q)
en 2.4( Q )worden uitgedrukt in een Fourier-reeks over een omwentelingshoek met een periode az, van 27c volgens:
6
Theoretische afleiding van de aeluidsafstraline: van een baffle
waarin: 1
K A = -sinmepH, K& =
1 T
cos mOpH,,
met: 2~ ,-ikR(P,Q)
R(P,Q)
cosm0 do.
Door de introductie van deze Fourier-reeksen kan (2.3) nu worden geschreven als (met dS(Q) = ' Q deQ dLQ):
Rangschikken van termen en met gebruikmaking van het feit d a t de integralen over de termen sin mOQ sin nOQ, sin mOQ cos nOQ en cos mOQ cos nOQ slechts ongelijk nul zijn wanneer m = n) leidt na integratie over de hoek OQ tot:
Rangschikken van de coëffiënten in deze vergelijking leidt tot de volgende 3 vergelijkingen:
Uit deze vergelijkingen dienen voor elk gewenst punt P , $o(P),4 i ( P ) en )'I(& bepaald, wanneer +O'(&), Pn'(Q)en &'(Q) bekend zijn, of juist andersom. Tot dusver is er nog geen sprake van een benadering.
te worden
2.3 Constructie van de benaderingsoplossing voor de Rayleigh integraal
2.3
7
Constructie van de benaderingsoplossing voor de Rayleigh integraal
Een eerste benadering van (2.4) is het meenemen van slechts M Fourier-termen in plaats van eeI; oneir,dig aantal.
' s de "-(P)'s op het contactvlak. In de BEM som zijn we vooral geïnteresseerd in de ~ ( P ) en 8ZQ Deze zijn nodig voor de koppeling van het gebaffelde cilinder model aan BARD (zie hoofdstuk 3). Daartoe verdelen we het contactvlak op in N elementen. Daarbij wordt gewoonlijk gebruik gemaakt van kwadratische lijn-elementen ([i]). Met dit 3-knoopselement zal (2.4) worden gediscretiseerd.
2.3.1
Berekening van de Rayleigh integraal over het contactvlak
In paragraaf 2.2 hebben we gezien dat de Rayleigh integraal voor een axisymmetrische constructie met niet-axisymmetrische randvoorwaarden wordt gegeven door (2.4) :
Zoals zojuist vermeldt, wordt het contactvlak verdeelt in N elementen. We concentreren ons nu even op één element. Binnen dit element introduceren we het locale coördinatensysteem zoals weergegeven in figuur 2.2. Binnen dit element kan voor respectievelijk straal r(6) en
4
2 I
&1 r=rl
15
r=r2
3
a
t=l r=r3
Figuur 2.2: Locale coördinatensysteem Jacobiaan J ,
(6) geschreven
worden:
8
Theoretische afleiding van de geluidsafstraling van een baffle
O p dezelfde wijze kunnen nu ook de overige variabelen uit (2.5) worden uitgeschreven (element m):
Door nu de gedefinieerde kolommen te substitueren in (2.5) levert dit na sommatie over alle elementen:
met:
Hierbij dienen de integralen numeriek t e worden bepaald. Zoals nader in hoofdstuk 4 zal worden besproken, dient hier een Gauss-Legendre- en Log-Gauss-integratie schema te worden gehanteerd.
2.3.2
Berekenen van de integraal H,
Zoals in paragraaf 2.2 werd beschreven, wordt H , gegeven door:
met:
2.3 Constructie van de benaderingsoplossing voor de Rayleigh integraal
9
Hierin kunnen weer de gedefinieerde kolommen uit de vorige paragraaf worden gesubstitueerd, waardoor deze gebruikt kan worden in combinatie met (2.6). Bij het oplossen van deze integraal kan bijvoorbeeld gebruik worden gemaakt van de formule van Euler (trapezium regel), de formule van Simpson of een hogere orde integratie schema. Over welke methode het beste kan worden toegepast, wordt in hoofdstuk 4 nader ingegaan.
Hoofdstuk 3
Koppeling van Rayleigh integraal aan BEM Met behulp van (2.1):
C(P>$(P> =
JJS (Sb(Q>VQ'(P,Q>- Q(P,Q)V4(Q))ns d S ( Q ) *
is het in principe mogelijk om de akoestische afstraling van een MRI-systeem t e modelleren. Tenminste in het geval dat het golfgetal k in domein 1 en 2 gelijk zijn (zie figuur 3.1). Het op
P N elementen
o N elementen
Z
Figuur 3.1: Wijze van modellering van het MRI-systeem deze manier modelleren heeft echter twee nadelen:
i. In elk punt P dient bovengenoemde integraal geïntegreerd te worden over het oppervlak Si en
SZ.
2. Men kan zich afvragen, hoe men de vlakke oneindige plaat kan modelleren.
Dit hoofdstuk zal de koppeling van de Rayleigh integraal aan het waardoor de zojuist genoemde nadelen komen te vervallen.
BEM pakket BARD behandelen,
Koppeling van Rayleigh integraal aan BEM
12
3.1
Theoretische afleiding ten behoeve van koppeling
In plaats van slechts één integraalvergelijking uit t e gaan kan het akoestisch probleem ook worden opgedeeld in tweeën. Eén voor de modellering van de geluidsafstraling binnen de cilinder en één voor de modellering van de geluidsafstraling naar de vrije ruimte. De integraalvergelijking voor de akoestische afstraling binnen de c i h d e r beschreven middels:
-
domein 1 - wordt
Deze geluidsproduktie kan worden berekend met behulp van het BEM pakket BARD. Voor de geluidsafstraling naar de vrije ruimte wordt uitgegaan van de in hoofdstuk 2 afgeleide Rayleigh integraal (2.3):
O p Si gelden d e voorgeschreven randvoorwaarden, terwijl op s compatibiliteit wordt geëist. Dat wil zeggen dat op het contantvlak s voor de continuïteit van de normaalsnelheden wordt geëist dat:
en voor de continuïteit van de druk wordt geëist dat:
pidi(&) = ~ 2 4 2 ( Q ) -
3.2
Koppelen van axisymrnetrische constructies met niet-axisymmetrische randvoorwaarden
De binnen het BEM pakket BARD aanwezige integraalvergelijkingen zijn zoals gezegd gebaseerd op een theorie die de geluidsafstraling berekend voor axisymmetrische constructies veroorzaakt door niet-axisymmetrische trillingen. Ter beschrijving van deze niet-axisymmetrische trillingen zijn net zoals in hoofdstuk 2 Fourier-reeksen ingevoerd. Daarbij worden weer slechts M Fourier-coëfficiënten meegenomen. Ook wordt het vlak L I en I binnen domein 1opgedeeld iri Ni N kwadratische lijn-elementen. Voor domein 1 komen de integraalvergelijkingen er dan als volgt uit de t e zien ([l]):
+
voor n = 1,.. ., M ,
(3.1)
3.2 Koppelen van axisymmetrische constructies met niet-axisymmetrische randvoorwaarden
13
met:
en voor de Rayleigh integraal gelden de vergelijkingen (2.6). Deze vergelijkingen samen met de compatibiliteitseisen maken vervolgens de koppeling mogelijk. Op de numerieke implimentatie zal in hoofdstuk 4 uitvoerig worden ingegaan.
HooÍästuk 4
Numerieke implementatie In dit hoofdstuk zal besproken worden hoe het boundary elementen methode programma is opgebouwd. Daarnaast wordt ingegaan op de benodigde numerieke integratie-methoden. Ook wordt de numerieke koppeling van beide BEM sommen beschreven.
4.1
Programma opbouw
Voor de manier van opbouw is gekeken naar andere programma's. Neem bijvoorbeeld de programmatuur beschreven in [3], de .mud file van het EEM pakket MARC, de invoerfile van ANSYS en de programmatuur behorende bij het vak: 'Applied Computational Mechanics'. Het hier beschreven BEM pakket is globaal als volgt opgebouwd: 1. Het opvragen of creëren van een mesh met randvoorwaarden en deze weergeven op het scherm. De mogelijke randvoorwaarden zijn: 40, $;, $i, @,',4: en/of 4;' voorschrijven voor n Fourier-coëfficiënten. 2. Het invoeren van het golfgetal k en de materiaaleigenschappen (binnen het hoofdprogramma (RAYLE1GH.M) kan het aantal integratie punten ten behoeve van GaussLegendre- en Log-Gauss-integratie worden opgegeven en kan de tolerantie van de integratie van H , worden opgegeven).
3. Het bepalen van de Rayleigh integralen voor de elementen en deze assembleren tot de matrices G en H (dus de integralen in de knooppunten Q ) .
4. Oplossen van het stelsel vergelijkingen:
&)6(&) = G(PopQ ,&)a(&) of het koppel probleem
H ( P o pQ ,
Dit levert de onbekende
40,
&,, 4,;
c', 4; en/of 4:
op de elementen op.
5. Het bepalen van de Rayleigh integralen voor de veldpunten en deze assembleren tot de matrices G en H (dus de integralen in de veldpunten P ) . 6. Oplossen van het stelsel vergelijkingen:
H ( P , &)&P) = G(P,Q)@(Q)
Numerieke implementatie
16
Dit levert de onbekende
40,
&, + g', i, 4:
en/of
4: in de veldpunten op.
In figuur 4.1 staat de opbouw nog wat meer in detail weergegeven. Hierin is duidelijk het
M voor de afstralvan geluidsgolven
A
/
/
J Vraag een bestaande mesh op
Maak een nieuwe mesh
Vraag een bArd-mesh op
r i ,
Plot axisymmetrische mesh & Plot randvoorwaarden
Plot contact-elementen
I
1
Invoeren van materiaal eigenschappen
Invoeren van materiaal eigenschappen
I
f
Bepalen van Rayleigh-integralen voor contact-elementen & Assembleren tot matrices
Bepalen van Rayleigh-integralen voor rand-elementen & Assembleren tot matrices
I
1
Bepalen van drukken en snelheden op de contact-
Bepalen van drukken en snelheden op de randelementen
I Bepa!er! van Ray!eigh-integralen I voor veld-punten
Bepalen van drukken en
Einde BEM
Figuur 4.1: Opbouw van het BEM pakket onderscheid tussen een zuiver Rayleigh probleem en een gekoppeld probleem zichtbaar. In
4.2 Numerieke integratie
17
bijlage A staat hetzelfde schema nogmaals gegeven, maar nu is hierin aangegeven welke procedures achtereenvolgens worden aangeroepen.
4.2
Numerieke integratie
In paragraaf 2.2 is reeds ter sprake gekomen, dat voor het oplossen van d e Rayleigh integraal diverse integratie-methoden moeten worden toegepast. In deze paragraaf zullen deze methoden worden besproken.
4.2.1
Berekenen van de integraal H,
Bij het oplossen van de integraal:
kan gebruik worden gemaakt van de formule van Euler (trapezium regel), de formule van Simpson of een hogere orde integratie methode, zoals QUAD8 - een standaard MATLABprocedure - en GAUSSHM - een procedure die onder andere gehanteerd wordt binnen BARD. Deze integratie methoden worden in bijlage B kort besproken. In figuur 4.2 zijn de resultaten van de absolute waarde van H , (n = O) weergegeven als functie van Y-Q( r p = 1).In de linker grafiek staan de resultaten weergegeven wanneer voor de Euleren Simpson-integratie 24 intervallen worden genomen. In de rechter grafiek staat de relatieve fout van IH,] voor de diverse methoden ten opzichte van GAUSSHM weergegeven. Zoals uit
Figuur 4.2: Integratie van H , met Euler, Simpson, QUAD8 en GAUSSHM de grafieken blijkt, treedt bij Y-Q = 1 een singulariteit op (R(P,&)= IQ - PI = O). Hierbij treedt de singulatiteit slechts op in het reële deel van H,, want: lim e-ihR(p,Q) = lim cos(-kR(P,Q)) R-+o
R(P,&)
R-+û
R(P,&)
+
lim isin(-kR(P, R-to
R(P,Q)
Q))
= 00
+ io.
Daarnaast is waar te nemen dat de resultaten met behulp van de formule van Euler en Simpson afwijken van de resultaten met de hogere orde methoden. Een toename van het
Numerieke implementatie
18
aantal intervallen blijkt echter wel te leiden tot een afname van de relatieve fout. Doordat de procedure QUAD8 de mogelijkheid biedt tot het opgeven van een tolerantie, is deze bij de verdere numerieke implementatie toegepast.
Zoals uit de vorige paragraaf blijkt treedt er een i-singulatiteit op. De singulariteit treedt alleen op in het element met de node waarvoor R(P,Q ) = IQ - PI = O. Met deze wetenschap kunnen we nu twee gevallen onderscheiden: 1. integratie over elementen zonder &singulariteit
2. integratie over elementen met ;-singulariteit Volgens [3] is het daarbij gebruikelijk om voor geval 1 Gauss-Legendre-integratie toe te passen. Voor geval 2 dient echter een nauwkeurigere integratie methode t e worden toegepast, zoals een hogere orde integratie, of een speciale formule (bijvoorbeeld Log-Gauss) . Hier is gekozen voor Log-Gauss, omdat het reële deel van Hn zich in de buurt van de singulariteit logaritmisch gedraagt. Opgemerkt dient t e worden dat in plaats van Log-Gauss-integratie ook gewoon Gauss-Legendreintegratie kan worden toegepast. Echter het grote voordeel van Log-Gauss is dat de rekentijd sterk wordt gereduceerd bij gelijk blijvende nauwkeurigheid (minder aantal integratie-punten nodig).
Gauss-Legendre-integratie Gauss-Legendre-integratie is niets anders dan het invullen van een waarde (hier [) in de te integreren vergelijking, gevolgt door een vermenigvuldiging met een weegfactor wint. De waarde van [ en wintzijn afiiankelijk van het aantal integratie-punten. De Rayleigh integraal (2.6) komt er dan als volgt uit t e zien:
Deze integratie is terug t e vinden in de procedure GEXT.M. In bijlage C is een tabel met waarden van [ en wint gegeven voor een aantal integratie punten van i tot en met 12. Dit staat ook in de procedure GAUSS.M.
4.2 Numerieke integratie
19
Log- Gauss-integratie
De te integreren termen in de Rayleigh integraal luiden voor een kwadratisch element:
si(<)= -LniT;(~)H,(J)fQ(J)Jrn(J)
met i = 1,273.
Eierbij treedt de &ngp!arit& ~p in -Vn, zo& we 21 eerder in n z r' r a'sor-a'a f Z"ab" s"e n Wanriper P Q dicht nadert kan het reële deel van Hn beschouwd worden als een functie van de vorm In ($) ([3]). De integratie grenzen van deze functie lopen van O tot 1, in plaats van -1 tot 1.
Log-Gauss-integratie is van de vorm:
Deze t e integreren functie lijkt sterk op si((). Vandaar dat we deze Log-Gauss-integratie methode willen transformeren t o t in de vorm van de Rayleigh integraal. Daarbij zijn de in figuur 4.3 aangegeven gevallen te onderscheiden. In bijlage D zijn de transformaties voor deze drie gevallen uitgevoerd. In tabel 4.1 staan de
--geval 1
geval 2
C=O
C=-l q=O
C=l q=l
.=knooppunt
C)-1 q =1
Coördinaten-transformatie
1 2
i7= + ( J + l ) 7' = -Q en = f
q=
5'+0
Ti
5-1 q=l
C=l q=1
q=o
5=0
X=singulariteit
geval
3
geval 3
i(-(+1)
Log-Gauss-transformatie ~ ' 1 g(t) Jol
g(F') f l
de =
-2wint,ip
+ Jtg ( F ) d~ =
s(t)de = E
2
9 2% -1)
d-vi
)+g(i7ip) ~ 2 - W i n 1t , i p i?,,> 1 -2Wint,i?,d-217ipfl)
weergeven, zien we d a t - ten gevolge van het gedrag van de vormfuncties - de Log-Gaussintegratie alleen mag worden toegepast op het reële deel van H, in het knooppunt waarin de singulariteit optreedt (zie figuur 4.4). Voor de overige knooppunten en het imaginaire deel van Hn dient gewoon Gauss-Legendre-integratie t e worden toegepast. Deze integraties zijn terug t e vinden in de procedure GL0C.M. Standaard worden daarbij 4 integratie punten gebruikt. Het aantal integratie punten voor Gauss-Legendre-integratie kiezen we hieraan gelijk.
5=1 q=O
Numerieke implementatie
20
4
i
i
i
m
m
i
i
i
i
2-
.. .
0 -.<.
.i_-
,..,
. I
-2 -
--
-4-
Singulariteitin knwppunt3
U..
,. ...__*
-g d
-2
-4-
(r
-6
,
,
,
,
,
,
,
O
0.2
0.4
0.6
0.8
.:.,.
L3 -10 -1
-0.8
-0.6
4.4
42
n
Figuur 4.4: Grafische weergave van g; voor de drie gevallen In bijlage E zijn de waarden van 7 en wint in een tabel weergegeven voor een aantal integratie punten van 1 t o t en met 10. Dit staat ook in de procedure L0GGAUSS.M.
4.3
Koppeling van Rayleigh integraal aan BEM
Met de nu geïntroduceerde integratie-methoden kunnen de integraties:
uit (2.6) en (3.1) worden uitgevoerd. De punten P worden nu zo gekozen dat deze samenvallen met de knooppunten Q uit de domeinen 1 en 2. Op deze manier kunnen de bovenstaande integralen voor elk element in elk punt P worden bepaald. We definiëren nu de vectoren (5'O ? @n ? &,& en Hierin staan respectievelijk de bekende en onbekende normaalsnelheden en drukken in de knooppunten Q. Door de definitie van deze vectoren kunnen de integralen {g} in de punten P uit domein 1worden opgeslagen in de matrix G,, enerzijds en de integralen {g} in de punten P uit domein 2 in de matrix G,, anderzijds. In de matrices Hnl en H,, worden de waarden van C ( P )en { h }opgeslagen. Op deze manier kan het stelsel vergelijkingen
2.
40,
4.3 Koppeling van Rayleigh integraal aan BEM
21
voor domein 1 en 2 geschreven worden als:
Na zorgvuldig rangschikken kunnen de matrices G,, en H,, worden opgesplitst in een gedeelte d a t deel uitmaakt van het oppervlak L i en een gedeelte dat deel uitmaakt van het contactvlak 1. Met behulp van de in hoofdstuk 3 geïntroduceerde compatibiliteitseisen kan voor iedere Fourier-coëfficiënt n het volgende stelsel vergelijkingen worden gevonden die moeten worden opgelost :
Met behulp van (2.6) respectievelijk (3.1) kunnen met behulp van de nu bepaalde normaalsnelheden en drukken de akoestische afstraling binnen de cilinder en naar de vrije ruimte worden bepaald. Tot nu toe hebben we ons geconcenteerd op één opening. Het MRI-systeem heeft echter zowel links als rechts een opening. Op dezelfde manier kan hiervoor worden afgeleid:
[ [?I
1.
22
Numerieke implement at ie
Deze theorie staat ook binnen de MATLAB-procedure C0UPLE.M. Binnen C0UPLE.M worden G,,, GE, en Gn, gegeven door B/(-ikpc) en worden H,,, H E en H t l gegeven door A . De matrices A en B zijn berekend middels het BEM pakket BARD. Het verschil tussen de matrix G en B wordt veroorzaakt door een verschil in definitie van de basisvergelijking voor de BEM binnen BARD en RAYLEIGH.
Hoofdstuk 5
Verificatie met analytische oplossingen In dit hoofdstuk zullen enkele analytische oplossingen voor de Rayleigh integraal worden vergeleken met oplossingen verkregen met de hier ontwikkelde BEM. Tevens worden twee voorbeelden gegeven van een gekoppeld probleem.
5.1
Oscillerende zuiger in een vlakke oneindige plaat
Een oscillerende zuiger in een vlakke oneindige plaat is een probleem dat uitstekend door de Rayleigh integraal beschreven kan worden. Voor een dergelijk probleem bestaan, door het doen van enige aannamen, analytische oplossingen. Deze analytische oplossingen worden gebruikt om het hier ontwikkelde BEM pakket t e controleren.
5.1 .I
Geluidsdruk op het wandoppervlak voor laag frequente oscillaties
In een vlakke oneindige plaat bevindt zich een zuiger met een straal a van 1 [m]. Deze zuiger oscilleert met een zeer lage frequentie ( k a = 0.01 < 1 [-3). Met de aannamen d a t ka < 1 E-], dat 8, (snelheidsamplitude) constant is over het zuigeroppervlak (Fourier-coëfficiënt O is gelijk aan 2 ) en d a t we slechts geïnteresseerd zijn in de geluidsdruk j? op het wandoppervlak ( z = O), kan de Rayleigh integraal worden vereenvoudigd tot ([2], pagina’s 218-219):
waarin:
De integralen E ( m ) en K ( m ) worden respectievelijk de elliptische integralen van de eerste en tweede soort genoemd en worden gegeven door:
Verificatie met analytische oplossingen
24
In figuur 5.1 staat links de mesh weergegeven. Deze mesh bestaat uit 5 3-knoopselementen. Rechts staat de dimensieloze drukamplitude als functie van voor zowel de analytische
p”
a= 1 [m]
_
-
- ‘Z
’O
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
rla
Figuur 5.1: Mesh en resultaten van het drukverloop op het wandoppervlak voor een laag frequent oscillerende zuiger in een vlakke oneindige plaat (ka = 0.01 [-I)
als de BEM oplossing gegeven. De BEM oplossing is verkregen met de standaard instellingen van RAYLE1GH.M (aantal integratie = 4; tolerantie = le-3). Zoals uit d e resultaten blijkt, komt de BEM oplossing uitstekend overeen met de analytische oplossing. Het is echter voldoende klein wordt genomen. belangrijk d a t de frequentie f = Uit de resultaten blijkt dat bij een laag frequente oscillatie van een zuiger in een baffle de geluidsdruk op de as van symmetrie het grootste is. In radiale richting neemt de druk aan de wand af. De verklaring hiervoor is dat de lucht voor de zuiger op de symmetrie-as zich uitsluitend in axiale richting verplaatst, terwijl bij toename van T- de lucht steeds meer in radiale richting gaat wegvloeien. In tabel 5.1 staat de invloed van meshverfijning, aantal integratie-punten en tolerantie op het verkregen resultaat weergegeven voor zowel Gauss-Legendre- als Log-Gauss-integratie. Zoals GaussLegendre
LogGauss
codes 11 11 11 11 6 11 11 11 11 6
nint
4 2 6 4 4 4 2 6 4 4
tdrrantie le-3 le-3 le-3 le-6 le-3 le-3 le-3 le-3 le-6 le-3
relatieve fout [%] 3.1877e- 1 1.0570e+0 1.5449e-1 3.1972e-1 3.187?e- 1 1.5165e-1 8.2722e-1 6.0766e-2 1.5070e-1 1.3174e-1
Tabel 5.1: Invloed van meshverfijning, aantal integratie-punten en tolerantie op het verkregen resultaat uit de resultaten blijkt, kan, zoals eerder reeds opgemerkt, op de elementen het beste LogGauss-integratie worden toegepast. Een toename van het aantal integratie-punten leidt voor
5.1 Oscillerende zuiger in een vlakke oneindige plaat
25
beide integratie-methoden tot een verbetering van het resultaat. Meshverfijning leidt niet tot een noemenswaardige verbetering.
5.1.2
Geluidsdruk in het verre veld voor een zeer locale wandtrilling
vve DeScnouwen een vlakke ûneindigc. and wa2ri2 zich e r n osdlerende s ~ i g e rhevindt. met een straal a van 1 [m]. Deze zuiger trilt met een constante frequentie ( k u = 8 [-I) en heeft een snelheidsamplitude fin = 1 [m/s] (Fourier-coëfficiënt O is gelijk aan 2). Wanneer we aannemen d a t de afstand van de zuiger tot de veldpunten vele malen groter is dan a of ka2 (bijv. tangentiale afstand = 800a) en de snelheidsamplitude fin constant is over het zuigeroppervlak, kan de Rayleigh integraal worden benaderd door ([2], pagina’s 225-227) : -77
1
1
met:
m 4 >= -8
.pcÛnka2 2Jl(kasin0) kasin0 *
Jl wordt wel de Bessel functie van de eerste orde genoemd en wordt gegeven door: J l ( 7 ) = ~ / TO c o s ( q c o s ~ ) s i n 2 qd4. 5 De zuiger wordt gemodelleerd met 5 kwadratische elementen en in het verre veld liggen 30 veldpunten (zie figuur 5.2). In de rechter figuur staat met een getrokken lijn de analytische x 104
I a=l[m] I
x x
I
! !
rn
I
x
x
x
x
x x
-
x _.-.i-
Figuur 5.2: Mesh en resultaten van het drukverloop op grootte tangetiale afstand van een oscillerende zuiger in een vlakke oneindige plaat ( k a = 8 tangetiale afstand= 800a)
[-I,
oplossing weergegeven en met o de oplossing met behulp van de BEM aangegeven. De resultaten komen zeer goed overeen. Wordt de tangentiale afstand echter niet ver genoeg gekozen, dan treden in de buurt van z = O afwijkingen op. Deze afwijkingen worden veroorzaakt door de aannamen die binnen de analytische oplossing zijn gedaan. Het blijkt d a t bij toename van ka het aantal ‘druklobben’ toeneemt. Binnen een ‘lob’ plant de geluidsdruk zich gemakkelijk voort, terwijl daar buiten zich geen geluidsdruk voortplant.
Verificatie met analytische oplossingen
26
5.1.3
Geluidsdruk op de as van symmetrie
Beschouwd wordt een vlakke oneindige wand waarin zich een oscillerende zuiger bevindt met een straal a van 1 [m]. Deze zuiger trilt met een constante frequentie (ka = 1 1 ~ en heeft een snelheidsamplitude U, = 1 [m/s] (Fourier-coëfficiënt 0 is gelijk aan 2). Wanneer we aannemen d a t de veldpunten zich uitsluitend op de symmetrie-as bevinden, kan de Rayleigh integraal worden vereenvoudigd tot ([2], pagina’s 232-233) :
[-I)
Links in figuur 5.3 staat de gebruikte mesh gegeven. In de rechter figuur staat de dimensieloze
Figuur 5.3: Mesh en resultaten van het drukverloop op de symmetrieas voor een oscillerende zuiger in een vlakke oneindige plaat ( k a = 1 1 [-I) ~
p”
drukamplitude als functie van de axiale afstand voor de analytische en de BEM oplossing weergegeven. De resultaten komen goed overeen. Alleen dicht bij het oppervlak zien we d a t er een kleine fout (relatieve fout is f 2 . 4 3 [%o]) optreedt. Door gebruik van meer integratie-punten (standaard 4)’ kan deze fout worden verkleind. Zo levert bijvoorbeeld 5 integratie-punten een relatieve fout van 60.085 [%o]. (Let op: Bij gebruik van kwadratische elementen kunnen alleen een even aantal integratie-punten worden gekozen.) Zoais waarneembaar is, blijkt dicht bij het zuiger oppervlak een sterke druk-interferentie op te treden. In het verre veld echter neemt de axiale druk monotoon af (asymptotisch volgens: $= 1 2 PcÛ,~ka).
5.1.4
De ontwikkeling van ‘druklobben’ bij toenemende axiale afstand
Er wordt uitgegaan van een vlakke oneindige wand waarin zich een oscillerende zuiger bevindt met een straal a van 1 [m]. Deze zuiger trilt met een constante frequentie van ka = 20 en heeft een snelheidsamplitude 8, = 1 [m/s] (Fourier-coëfficiënt 0 is gelijk aan 2). We zijn nu geïnteresseerd in de ontwikkeling van de ‘druklobben’ bij toenemende axiale afstand. Deze
[-I
5.1 Oscillerende zuiger in een vlakke oneindige plaat
27
kan beschreven worden met de volgende analytische oplossing ([2], pagina's 236-237) :
met: 1 voor r O voor r
H=
a
+
R, = [ ( a 21; R1 = [ ( a w)2 + 212
+
X = H wordt de Heaviside-eenheidstapfunctie genoemd en R, en RI zijn respectievelijk de kleinste en de grootste afstand van een veldpunt tot de omtrek van de zuiger. AD is de Diffusie integraal en kan met de Fresnel integralen: C ( X )=
1
1
X
X
cos (;t2)
dt
S ( X )=
sin ( i t 2 ) dt,
worden geschreven als:
Links in figuur 5.4 staat de mesh met veldpunten op een axiale afstand van z/a =0.5, 1.0, lP21 als 1.5 en 2.0 weergegeven. In de rechter figuur staat de dimensieloze drukamplitude ____ (Pc)21%l I I
I
I I
! ! I
! I I
! I I
x
Y
Y
x
Y
Y
x x Y Y Y
x Y
x
x x Y Y
x Y Y
x Y
x Y Y
x Y
x Y
Y
a=hi Y Y Y Y Y
x
Figuur 5.4: Mesh en resultaten van de visualisatie van de ontwikkeling van 'druklobben' bij toenemende axiale afstand voor een oscillerende zuiger in een vlakke oneindige plaat ( k a = 20
[-I ) functie van de radiale afstand f voor de analytische en de BEM oplossing weergegeven. De resultaten zijn uitstekend.
Verificatie met analytische oplossingen
28
De resultaten geven duidelijk de ontwikkeling van de ‘druklobben’ weer. Een toename van de axiale afstand gaat gepaard met een verschuiving van de ‘lobben’ in radiale richting.
5.2
Geluidsafstraling van een pijp m-et een openhg aan één zijde
Aan de ene zijde van een pijp bevindt zich een oscillerende zuiger en aan de andere zijde een opening. De binnenstraal van de pijp is gelijk aan de straal van de zuiger en heeft een straal a van 1 [m]. De buitenstraal van de pijp is 2a en de pijplengte I is 3a. De zuiger oscilleert met een constante dimensieloze frequentie ka van 0.85 [-l. De snelheidsamplitude 6, is 1 [m/s] (Fourier-coëfficiënt O is gelijk aan 2). In het artikel [5] zijn resultaten voor een dergelijk probleem gegeven. Deze zijn verkregen op basis van een BEM d a t het volledige probleem (‘enkel-domein’) kan modelleren. Met het BEM pakket BARD en de hier ontwikkelde BEM is een koppeling voor dit probleem uitgevoerd. De pijp is gemodelleerd binnen BARD wat de matrices A en B en een mesh oplevert. Aan de ene zijde van de pijp wordt de constante snelheidamplitude U, voorgeschreven, terwijl de andere zijde wordt gekoppeld aan een baffle. Het beschouwen van het uiteinde van de pijp als zijnde een baffle is een benadering. De resultaten van beide methoden staan in figuur 5.5 weergegeven. Links staat de dimen-
10.9
-
0.8 -
--
20.7-
e
-YaOB 0.5 -
-gekoppelde BEM oplming (balRe) 1 aplmring cilindrischepijp -1.5
-1
4.5
O
ua
0.5
X
1
1.5
0‘
oplossingcilirdiischa pijp 0.005
0.01
0.015
0.02 0.025 lpl/(~twclvnl)
0.03
0.035
0.04
15
Figuur 5.5: Resultaten van de geluidsdruk binnen de pijp en de afstraling naar de vrije ruimte ( k a = 0.85 [-I; tangentiale afstand =loa)
&
sieloze drukamplitude op de pijpwand als functie van de plaats weergegeven. Rechts is de drukamplitude op een tangentiale afstand van 10a weergegeven. Er treden behoorlijke verschillen op in de resultaten. Door het gebruik van ‘Rayleigh-elementen’ treedt er meer reflectie op, dan bij de pijp-opening met een buitenstraal van 2a. Hierdoor is de druk aan de buitenwand voor het gekoppelde probleem hoger. Deze extra reflectie leidt tevens tot een verandering van de drukamplitude binnen de pijp. Dit verklaart de verschillen in de resultaten.
5.3 Pulserende cilinderwand in een akoestisch baffle
5.3
29
Pulserende cilinderwand in een akoestisch baffle
Binnen BARD is een optie (‘BaffledDuct’) aanwezig om de geluidsafstraling binnen een gebaffelde cilinder t e modelleren met niet-axisymmetrische randvoorwaarden. Met behulp van deze optie zal het gebruik van niet-axisymmetrische randvoorwaarden binnen de koppeling tiissen het RER4 =&ket BARD en de hier ontwikkelde REM worden gecontroleerd. Daartoe wordt binnen BARD de cilinder gemodelleerd. Dit levert de vereiste matrices A en B en een mesh. Vervolgens worden beide zijden gekoppeld aan de ‘Rayleigh-elementen’ en wordt op de cilinderwand de gewenste mode - Fourier-coëfficiënt - voorgeschreven. Deze procedure is uitgevoerd voor mode O, 1 en 3 voor de frequenties ka = 0.85 [-I en ka = 10 [-l. Voor mode 1 en 3 worden alleen de cosinus-termen van de Fourier-coëfficiënt voorgeschreven. De cilinder heeft een lengte van 3 [m] en heeft een binnenstraal u van 1 [m]. De dichtheid p en de geluidsnelheid c zijn voor de eenvoud gelijk aan 1 gekozen. In figuur 5.6 staan de resultaten voor de verschillende modes en frequenties weergegeven. Per rij staan respectievelijk de modes O, 1 en 3 gegeven, terwijl in de afzonderlijke kolommen respectievelijk de frequenties ku = 0.85 [-I en ku = 10 [-] staan. Mode 1 is alleen voor een frequentie van ka = 0.85 [-I uitgevoerd. De resultaten komen voor lage frequenties zeer goed overeen. Voor hogere frequenties zien we dat bij een toenemend mode-nummer afwijkingen gaan optreden. Het verhogen van het aantal integratie-punten en/of de tolerantie leidt daarbij niet tot noemenswaardige verbeteringen. We zien ook dat bij toename van 5 en n de periodetijd van een trilling afneemt, waardoor ook de golflengte afneemt. Een afname van de golflengte bij constant blijvende elementlengte leidt tot het minder nauwkeurig kunnen beschrijven van een trilling. Zo geldt er de regel: ‘Er zijn minstens 6 knooppunten per golflengte noodzakelijk om een trilling goed t e kunnen beschrijven’. Meshverfijning leidt echter alleen tot een visuele verbetering en niet tot het kleiner worden van de afwijkingen. Bij de resultaten die verkregen zijn met de BARD-optie ‘BaffledDuct’, is een trend waarneembaar d a t bij toename van het aantal Gauss-integratie-punten de afwijking kleiner wordt. Het grote voordeel van de gekoppelde BEM is dat het de mogelijkheid biedt om zowel de geluidsafstraling binnen de cilinder als naar de vrije ruimte te beschrijven.
Verificatie met analvtische oDlossinnen
30
mads=O; k a 8 5
m&e=O; k=10
8
I
1
‘I
I -!.5
6
-1
I
o
-0.5
0.5
1
1.5
O
-1.5
-1
-0.5
O
0.5
1
1.5
0.5
1
1.5
z
2
made=l: ki0.85
p a v (ab$ gekBEM 1
05
O
z
mcde=3,k=O 85
0.8 -
0.6
mods=3; k=10
p 2. v (abs) gikBEM
-
0.4 ., . .... .. .. .
. .
l
-?.5
-1
-0.5
O
. . .
0.5
-
0.6
-
0.4
1.5
-1.5
. ..
1
-i
-0.5
O
Figuur 5.6: Resultaten van de koppeling tussen BARD en de ‘Rayleigh-elementen’ voor mode O, 1 en 3 en d e frequenties ka = 0.85 [-] en ka = 10 [-l.
Hoofdstuk 6
Conclusies en aanbevelingen 6. P
Condusks
Met het hier ontwikkelde Boundary Elementen Methoden (BEM) pakket binnen MATLAB is het mogelijk om de geluidsafstraling vanuit de uiteinden van een cilindrische constructie in een akoestische baffle te beschrijven. Daarnaast bestaat de mogelijkheid om deze formulering t e koppelen aan het BEM pakket BARD. Hiermee hebben we een gereedschap om de geluidsproduktie t e berekenen voor zowel binnen als buiten een gebaffelde cilinder veroorzaakt door niet-axisymmetrische trillingen. Dit model kan mogelijk uitkomst bieden voor het beschrijven van de akoestische afstraling van een MRI-systeem. Voor de modellering van het akoestisch baffle is uitgegaan van de zogenaamde ‘Rayleighintegraal’. Voor numerieke implementatie van deze integraal is over gegaan o p een BEMformulering. Omdat het BEM pakket BARD gebaseerd is op een theorie voor de beschrijving van de akoestische afstraling voor axisymmetrische constructies veroorzaakt door nietaxisymmetrische trillingen, is de hier ontwikkelde BEM-formulering hier ook op gebaseerd. De koppeling is verkregen door compatibiliteit te eisen op de contactvlakken van beide BEM sommen. Zoals uit de resultaten blijkt kan voor de bepaling van de ‘Rayleigh-integralen’ in de singuliere knooppunten o p de elementen het beste Log-Gauss-integratie worden toegepast. Daarnaast blijkt d a t bij een toename van het aantal integratie-punten de nauwkeurigheid van het eindresultaat wordt verbeterd. Ook het verkleinen van de tolerantie leidt tot een verbetering. Daar staat echter tegenover dat de rekentijd toeneemt. De standaard instelling, = 4 en tol = i e - 3, is echter een goed compromis. Meshverfijning is zinvol wanneer er minder dan 6 knooppunten per golflengte zijn. Er zijn namelijk minimaal 6 knooppunten nodig om een trilling t e kunnen beschrijven.
Conclusies en aanbevelingen
32
6.2
Aanbevelingen
Het hier ontwikkelde BEM pakket op basis van ‘Rayleigh-elementen’ is in principe een gereedschap, d a t zonder enige beperking met BARD kan worden gekoppeld. Toch tot slot nog enkele aanbevelingen: o
Tijdens d e ontwikkeling van het BEM pakket RAYLE1GH.M is gebleken dat MATLAB, vanwege de lange rekentijden, niet de ideale omgeving is om de modellering uit te voeren. Een vertaalslag naar een programmeertaal als C++ kan dit probleem verhelpen.
o
Een tweede reden om over t e gaan op de programmeertaal C++, is dat de koppeling direkt binnen BARD (geschreven in C++) kan worden uitgevoerd. Dit zal leiden tot een betere efficiëntie van de BEM en maakt het mogelijk om gebruik te maken van slechts één invoer-file.
o
Onderzocht moet worden of het gebaffelde cilinder model gebruikt kan worden voor de modellering van de akoestische afstraling van een MRI-systeem. Dit onderzoek vereist mogelijk de eerste prioriteit.
Bibliografie B. Soenarko. A boundary element formulation f o r radiation of acoustic waves from axisymmetric bodies with arbitrary boundary conditions. Journal of the Acoustical Society of America, Vol. 93, Nr. 2, pag. 631-639,1993.
A.D. Pierce. A C 0 USTICS An Introduction to Its Physical Principles and Applications. McGraw-Hill Book Company, New York, 1981, ISBN 0-07-049961-6. C.A. Brebbia en J. Dominguez. Boundary Elements An Introductory Course. Computational Mechanics Publications, Southampton, 1992, ISBN 1-85312-160-6. R.D. Ciskowski en C.A. Brebbia. Boundary Element Methods in Acoustics. Elsevier Applied Science, London, New York, 1991, ISBN 1085166-679-6. A.F. Seybert, C.Y.R. Cheng en T.W. Wu. The solution of coupled interior/exterior acoustic problems using the boundary element method. Journal of the Acoustical Society of America, Vol. 88, Nr. 3, pag. 1612-1618,1990.
Bijlage A
Opbouw van het BEM-pakket (MATLAB-procedures) De geschreven procedures kunnen grofweg in 7 groepen worden verdeeld: o
hoofdprogramma (Rayleigh function)
o
creëren van mesh-procedures (mesh function)
e
plot-procedures (plot function)
o
vormfunctie-procedures (shape function)
o
hulp-procedures (utility function)
o
koppel-procedures (utility function)
o
extra-procedures (ten behoeve van controle) (extra function)
In figuur A.l staat het schema uit hoofdstuk 4 nogmaals weergegeven. Nu wordt echter aangegeven welke procedures achtereenvolgens worden aangeroepen wanneer het schema wordt doorlopen. Deze procedures zijn op diskette bijgevoegd. Binnen MATLAB kan met het commando: ‘HELP procedure-name’, een korte verklaring van de functie gevraagd worden. Als voorwaarden voor het kunnen werken met RAYLEIGH.M, dient eerst het programma START.M t e worden opgestart. Dit programma zorgt voor het actief maken van de juiste werk-directory’s. Vervolgens kan RAYLE1GH.M worden gestart. Bij wijziging van directorystructuur dienen zowel in START.M als in RAYLE1GH.M de juiste directory’s te worden gezet.
Opbouw van het BEM-pakket (MATLAB-procedures)
36
t
geheqgen. hline3.m
sshm.m .m ass.m
Figuur A.l: Opbouw van het BEM pakket (MATLAB-procedures)
Numerieke integratie B .1 Formule van Eder (trapeziumregel) Interval (a,b ) is verdeelt in n gelijke intervallen, zoals in figuur B.l staat weergegeven.
Figuur B.l: Numeriek integrerei?
Bij een toename van het aantal intervallen neemt de nauwkeurigheid van de integratie toe.
B.2
Formule van Simpson
Interval ( a , b) is verdeelt in n gelijke intervallen (zie figuur B.l).
Numerieke integratie
38
Het aantal intervallen moet even zijn! Ook hier neemt de nauwkeurigheid van d e integratie toe bij een toenemend aantal intervallen.
B.3
QUAD8 - standaard MATLAB-procedure
De procedure QUAD8.M is eer, hogere orde numerieke integratie methode. QUAD8 biedt de mogelijkheid om een tolerantie (standaard = 1 zodoende rekening met het verloop van de functie.
B.4
GUASSHM
op te geven. Deze procedure houdt
- BARD procedure
Deze GAUSSHM.M routine berust op een Gauss-integratie (besproken in hoofdstuk 4). Ook hier geldt d a t bij het vergroten van het aantal intervallen de nauwkeurigheid toeneemt. De helling van d e functie H , is daarbij van invloed op de nauwkeurigheid. Een toename van de helling leidt tot een afname van de nauwkeurigheid. Zodoende wordt binnen deze routine bij een stijler wordende helling het aantal intervallen vergroot. Het vergroten van het aantal intervallen is ook noodzakelijk, wanneer het golfgetal k en de Fourier-coëfficiënt n toeneemt, omdat de periodetijd van H , dan kleiner wordt.
Bijlage C
Waarden van en wint ten behoeve van Gauss-Legendre-integratie nip -
.Wint,ip
.tip
1 0.00000 00000 00000 2 -0.57735 02691 89626 0.57735 02691 89626 3 -0.77459 66692 41483 0.00000 00000 00000 0.77459 66692 41483 4 -0.86113 63115 94053 -0.33998 10435 84856 0.33998 10435 84856 - 0.86113 63115 94053 5 -0.90617 98459 38664 -0.53846 93101 05683 c!.c!c!ooo 00000 00000 0.53846 93101 05683 0.90617 98459 38664 6 -0.93246 95142 03152 -0.66120 93864 66265 -0.23861 91860 83197 0.23861 91860 83197 0.66120 93864 66265 - 0.93246 95142 03152 7 -0.94910 79123 42759 -0.74153 11855 99394 -0.40584 51513 77397 0.00000 00000 00000 0.40584 51513 77397 0.74153 11855 99394 0.94910 79123 42759
~
~
2.00000 00000 00000 1.00000 00000 00000 1.00000 00000 00000 0.88888 88888 88888 0.55555 55555 55555 0.88888 88888 88888 0.34785 48451 37454 0.65214 51548 62546 0.65214 51548 62546 0.34785 48451 37454 0.23692 68850 56189 0.47862 86704 99366 0.56888 88888 88889 0.47862 86704 99366 0.23692 68850 56189 0.17132 44923 79170 0.36076 15730 48139 0.46791 39345 72691 0.46791 39345 72691 0.36076 15730 48139 0.17132 44923 79170 0.12948 49661 68870 0.27970 53914 89277 0.38183 00505 05119 0.41795 91836 73469 0.38183 00505 05119 0.27970 53914 89277 0.12948 49661 68870
40
Waarden van f en
wint
ten behoeve van Gauss-Legendre-integratie
.Wint ,ip
.tip
-0.96028 98564 97536 0.10122 85362 90376 -0.79666 64774 13627 0.22238 10344 53374 n CO KC^ onnnn i c m o 0.31379 66458 77887 -u.dLiddV i Y U . 7 . 7 -0.18343 46424 95650 0.36268 37833 78362 0.18343 46424 95650 0.36268 37833 78362 0.52553 24099 16329 0.31370 66458 77887 0.79666 64774 13627 0.22238 10344 53374 0.96028 98564 97536 0.10122 85362 90376 -0.96816 02395 07626 0.08127 43883 61574 -0.83603 11073 26636 0.18064 81606 94857 -0.61337 14327 00590 0.26061 06964 02935 -0.32425 34234 03809 0.31234 70770 40003 0.00000 00000 00000 0.33023 93550 01260 0.32425 34234 03809 0.31234 70770 40003 0.61337 14327 00590 0.26061 06964 02935 0.83603 11073 26636 0.18064 81606 94857 0.96816 02395 07626 0.08127 43883 61574 0.9739065285 17172 0.06667 13443 08688 -0.86506 33666 88985 0.14945 13491 50581 -0.67940 95682 99024 0.21908 63625 15982 -0.43339 53941 29247 0.26926 67193 09996 -0.14887 43389 81631 0.29552 42247 14753 0.14887 43389 81631 0.29552 42247 14753 0.43339 53941 29247 0.26926 67193 09996 0.67940 95682 99024 0.21908 63625 15982 0.86506 33666 88985 0.14945 13491 50581 0.97390 65285 17172 0.06667 13443 08688 -û.98156 U6342 46719 Û.Û4717 53363 865112 -0.90411 72563 70475 0.10693 93259 95318 -0.76990 26741 94305 0.16007 83285 43346 -0.58731 79542 86617 0.20316 74267 23066 -0.36783 14989 98180 0.23349 25365 38355 -0.12523 34085 i1469 0.24914 7Û458 13403 0.12523 34085 11469 0.24914 70458 13403 0.36783 14989 98180 0.23349 25365 38355 0.58731 79542 86617 0.20316 74267 23066 0.76990 26741 94305 0.16007 83285 43346 0.90411 72563 70475 0.10693 93259 95318 0.98156 06342 46719 0.04717 53363 86512 IUV.&.J
Tabel C.l: Gauss-Legendre-integratie
Bijlage D
Log-Gauss-transformatie Log-Gauss-integratie is van de vorm:
Deze integratie methode willen we transformeren tot in de vorm van de Rayleigh integraal.
D.l
Geval 1
In figuur D.l staan de beide locale coördinatensystemen van geval 1weergegeven. Met behulp
geval 1 P 5-1 q=O
e=knooppunt
5=0
4=1 q=1
X=singulariteit
Figuur D.l: Geval 1 hiervan kunnen de volgende transformatie-vergelijkingen worden afgeleid. Tevens is ook diens afgeleide bepaald: 1 7=s(E+1)
e <=2’1-1,
Log-Gauss-t ransformatie
42
Substitutie van de uitdrukkingen in de definitie vergelijking van Log-Gauss levert (f ( q i P )=
1.
d217ip-1) ln(l/vip)
D.2
Geval 2
In figuur D.2 staan de beide locale coördinatensystemen van geval 2 weergegeven. Hiermee
geval 2
I q' =1
.=knooppunt
5'+0
Ti =q=o
k=l q=l
Xsingulariteit
Figuur D.2: Geval 2 kunnen de voigende transformatie-vergelijkingen en diens afgeleide worden afgeleid:
Substitutie van de uitdrukkingen in de definitie vergelijking van Log-Gauss levert (f(@
rO
Pl
=
D.3 Geval 3
D.3
43
Geval 3
-
In figuur D.3 staan de beide locale coördinatensystemen van geval 3 weergegeven. Nu gelden
geval 3
5-1 q=1
.=knooppunt
5=1 q=O
5=0
X=singulariteit
Figuur D.3: Geval 3 de volgende transformatie-vergelij kingen en afgeleide: 1 2
7/=-(-[+1)
w
[=-27+1,
Substitutie van deze uitdrukkingen in de definitie vergelijking van Log-Gauss levert ( f ( q i p )=
Bijlage E
Waarden van 1;1 en wint ten behoeve van Log-Gauss-int egrat ie .Wint ,ip
I 0.60227691 3
4
5
6
7
0.63890792e-1 0.36899706 0.76688030 0.41448480e-1 0.24527491 O .556 16545 0.84898239 0.29134472e-1 O. 17397721 0.41170251 c!.57731417 0.89477136 0.216344005e-1 0.12958339 O .31402045 0.53865721 0.75691533 0.92266884 0.16719355e-1 0.10018568 0.24629424 0.43346349 0.63235098 0.81111862 0.94084816
0.71853931 0.28146068 O .5 1340455 0.39198004 0.94615406e-1 0.38346406 0.38687532 0.19043513 O .39225487e-l 0.29789346 0.34977622 0.23448829 0.98930460e-1 0.18911552e-1 0.23876366 0.30828657 O -24531742 0.14200875 O .55454622e- 1 0.10168958e-1 0.19616938 0.27030264 0.23968187 0.16577577 0.88943226e-1 0.33 194304e-1 0.59327869e-2
I
46
Waarden van 7 en wint ten behoeve van Log-Gauss-integratie
*Yip
0.13320243e-1 0.79750427e-1 0.19787102 0.35415398 0.52945857 0 .70 181452 0.84937932 0.95332645 0.10869338e-1 0.64983682e-1 0.16222943 0.29374996 0.44663195 0 .60548172 0.75411017 0.87726585 0.96225056 0.90425944e-2 0.53971054e-1 0.13531134 0.24705169 0.38021i71 0.52379159 0.66577472 0.79419019 0.89816102 0.96884798
*Wint,ip
0.16441660 0.23752560 0.22684198 0.17575408 0.11292402 0.57872212e-1 0.20979074e-1 0.36864071e-2 0.14006846 0.20977224 0.21142716 0.17715622 0.12779920 0.78478879e-1 0.39022490e-1 0.13867290e-1 0.24080402e-2 0.12095474 0.18636310 0.19566066 0.17357723 0.i3569597 0.93647084e-1 0.55787938e-1 0.27159893e-1 0.95151992e-2 0.1638i586e-2
Tabel E.l: Log-Gauss-integratie
