Západočeská univerzita v Plzni Fakulta aplikovaných věd Katedra mechaniky. Optimalizace akustického pole

Za´padoˇceska´ univerzita v Plzni Fakulta aplikovanyćh vˇed Katedra mechaniky

´ RSK ˇ ´ PRACE ´ BAKALA A Optimalizace akustick´ eho pole

Plzeˇ n 2013

Zdenˇek Novotn´ y

Prohl´ aˇ sen´ı Prohlaˇsuji, ˇze jsem bakaláˇrskou práci vypracoval samostatnˇe a v´ yhradnˇe s pouˇzit´ım citovan´ ych pramen˚ u. V Plzni dne 28. kvˇetna 2013 Zdenˇek Novotn´ y

Podˇ ekov´ an´ı Pˇri této pˇr´ıleˇzitosti bych rád podˇekoval vedouc´ımu bakaláˇrské práce Prof. Dr. Ing. Eduardu Rohanovi, DrSc. za motivaci ke studiu a veden´ı práce a Ing. Vladim´ıru Lukeˇsovi, Ph.D. za vˇsechny cenné rady a vˇenovan´ y ˇcas.

2

Abstrakt Hlavn´ım c´ılem této práce je proveden´ı tvarové optimalizace akustického pole popsaného omezenou 2D oblast´ı Ω. Nejdˇr´ıve se odvodily rovnice akustiky potˇrebné k z´ıskán´ı vlnové rovnice pro tekutiny. Z tˇechto rovnic se urˇcila vlnová rovnice. Hledán´ım jednofrekvenˇcn´ıho ˇreˇsen´ı se z n´ı odvodila Helmholtzova rovnice. Pak se zde ukázaly moˇzné pˇr´ıstupy k ˇreˇsen´ı vlnové a Helmholtzovy rovnice pro 1D kontinuum, kde se vyuˇzilo Fourierovy metody. Z u ´vah o ˇs´ıˇren´ı akustick´ ych vln se odvodily okrajové podm´ınky na hranici ∂Ω. Naformulovala se u ´loha tvarové optimalizace. Jako stavovou rovnici se vyuˇzila slabá formulace Helmholtzovy rovnice. Provedla se citlivostn´ı anal´ yza. Zparametrizoval se design oblasti Ω parametry α pomoc´ı spline-boxu. Optimalizaˇcn´ı u ´loze se pˇriˇradila Lagrangeova funkce a posléze se pˇreˇslo k u ´loze adjungované. Touto cestou se z´ıská celková derivace neboli citlivost u ´ˇcelové funkce Φ na zmˇenu parametr˚ u α. Dále jsou porovnány pro kontrolu v´ ysledky citlivostn´ı anal´ yzy a koneˇcn´ ych diferenc´ı. Pomoc´ı softwaru SfePy a Matlabu jsme provedli nˇekolik optimalizaˇcn´ıch v´ ypoˇct˚ u ve 2D. V´ ysledky se zobrazily programem ParaView a následnˇe vyhodnotily.

Abstract The main objective of this study is to implement shape optimization of an acoustic field described by 2D domain Ω. The very first step was the familiarization with derivation of acoustic equations necessary to obtain wave equation for fluids. From these equations, the wave equation was derived. By searching for a single-frequency solution, the Helmhotz equation was found. Possible approach in solving wave and Helmholtz equation for 1D continuum, where Fourier method was used, was shown. From our knowledge about acoustic propagation boundary conditions at boundary ∂Ω were derived. A problem of shape optimization was defined. As a state equation weak formulation of Helmhotz equation was applied. Then the sensitive analysis was performed. Design of Ω domain was parameterized by parameters α using spline-box. Lagrange equation was assigned to the optimization problem. Consecutively, there was a proceeding to solve an adjugate problem. This way allowed to obtain a total derivation of objective function that is a sensitivity of objective function to the change of parameters α. Furthermore, there were compared results from sensitive analysis and finite differences. Using software SfePy and Matlab we performed a few optimization calculations in 2D. At the end the results were displayed with software ParaView and also appraised.

Obsah ´ 1 Uvod 2 Odvozen´ı rovnic akustiky ´ 2.1 Uvodem . . . . . . . . . . . . . . . . . . . 2.2 Eulerova rovnice . . . . . . . . . . . . . . 2.3 Rovnice kontinuity . . . . . . . . . . . . . 2.4 Stavová rovnice . . . . . . . . . . . . . . . 2.5 Vlnová rovnice v kartézsk´ ych souˇradnic´ıch 2.6 Helmholtzova rovnice . . . . . . . . . . . .

1

. . . . . .

3 3 3 6 7 8 10

ˇ sen´ı u 3 Reˇ ´ loh pro 1D kontinuum 3.1 Vlnová rovnice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Helmholtzova rovnice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

11 11 17

4 Form. u ´ loh akust. pro danou frekv. ve 3D 4.1 Okrajová u ´loha . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Slabé ˇreˇsen´ı . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Numerické ˇreˇsen´ı MKP . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

21 21 22 23

5 Formulace u ´ loh tvarov´ e optimalizace ´ 5.1 Uvod . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Citlivostn´ı anal´ yza . . . . . . . . . . . . . . . . . . . . . . 5.2.1 Metoda adjungované promˇenné . . . . . . . . . . . 5.2.2 V´ ypoˇcet citlivosti pomoc´ı materiálové derivace . . . 5.3 Implementace . . . . . . . . . . . . . . . . . . . . . . . . . 5.3.1 Adjungovaná u ´loha . . . . . . . . . . . . . . . . . . 5.3.2 Pˇr´ıprava oblasti Ω . . . . . . . . . . . . . . . . . . 5.3.3 Spline-box . . . . . . . . . . . . . . . . . . . . . . . 5.3.4 Porovnán´ı koneˇcn´ ych diferenc´ı a citlivostn´ı anal´ yzy 5.3.5 Optimalizace . . . . . . . . . . . . . . . . . . . . . 5.4 Optimalizaˇcn´ı v´ ypoˇcty . . . . . . . . . . . . . . . . . . . . ´ 5.4.1 Uvodem . . . . . . . . . . . . . . . . . . . . . . . . 5.4.2 Vlastn´ı v´ ypoˇcty . . . . . . . . . . . . . . . . . . . .

27 27 28 28 30 31 31 32 32 34 35 35 35 36

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . . . . . . . . .

. . . . . .

. . . . . . . . . . . . .

. . . . . .

. . . . . . . . . . . . .

. . . . . .

. . . . . . . . . . . . .

. . . . . .

. . . . . . . . . . . . .

. . . . . .

. . . . . . . . . . . . .

. . . . . .

. . . . . . . . . . . . .

. . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

OBSAH

OBSAH

6 Z´ avˇ er

40

5

´ 1 Uvod ˇ sen´ı problém˚ Reˇ u akustiky, napˇr´ıklad tvarovou optimalizac´ı, má potencionáln´ı aplikace v mnoha oblastech mechaniky. Napˇr´ıklad pˇri optimalizaci ˇcást´ı stroj˚ u a mechanick´ ych zaˇr´ızen´ı, která sv´ ym chodem zp˚ usobuj´ı pˇr´ıliˇsnou hluˇcnost. Pak je nutné tyto hlukové emise sn´ıˇzit. Nicménˇe tato práce nen´ı vázána na konkrétn´ı u ´lohu nebo aplikaci, ale zab´ yvá se optimalizac´ı sp´ıˇse teoreticky. V druhé kapitole se seznám´ıme s odvozen´ım rovnic potˇrebn´ ych k z´ıskán´ı vlnové rovnice pro tekutiny. Tyto rovnice odvod´ıme s pˇredpokladem, ˇze prostˇred´ı je stlaˇcitelné, spojité, homogenn´ı, izotropn´ı a neviskózn´ı. Akustické pole budeme pˇredpokládat za nev´ırové. Odvod´ıme si Eulerovu pohybovou rovnici pro proudˇen´ı, která zanedbává tˇren´ı a vazkost. Rovnice kontinuity vyjadˇruje zákon zachován´ı hmoty. Posledn´ı stavová rovnice vycház´ı z popsán´ı termodynamického chován´ı plynu pˇri ˇs´ıˇren´ı zvukové vlny. Za dalˇs´ıch pˇredpoklad˚ u tˇechto rovnic urˇc´ıme tvar vlnové rovnice v kartézsk´ ych souˇradnic´ıch pro neznámou rychlostn´ı potenciál φ a akustick´ y tlak P . Pˇri hledán´ı jednofrekvenˇcn´ıho ˇreˇsen´ı vlnové rovnice pro frekvenci ω z n´ı odvod´ıme Helmholtzovu rovnici, která popisuje stojaté vlnˇen´ı. Zde ˇ ˇcerpáme zejména z [Skvor(2001)]. V tˇret´ı kapitole si ukáˇzeme nˇekteré pˇr´ıstupy k analytickému ˇreˇsen´ı vlnové a Helmholtzovy rovnice v jednorozmˇern´ ych pˇr´ıpadech, kde pouˇzijeme Fourierovu metodu ˇreˇsen´ı parciáln´ıch ˇ diferenciáln´ıch rovnic, viz [M´ıka(1983)] a [Skvor(2001)]. Budeme uvaˇzovat ˇs´ıˇren´ı zvukov´ ych vln mezi dvˇema rovnobˇeˇzn´ ymi dokonale odraziv´ ymi deskami. Z toho vyplynou okrajové a poˇca´teˇcn´ı podm´ınky. Pomoc´ı Matlabu si zobraz´ıme nˇekterá ˇreˇsen´ı. Ve ˇctvrté kapitole formulujeme u ´lohu akustiky pro danou frekvenci ω. Budeme tedy ˇreˇsit slabou formulaci, viz [M´ıka(2007)], Helmholtzovy rovnice. Zavedeme si oblast Ω. V softwaru MSC.Marc Mentat si vygenerujeme jej´ı MKP s´ıt’. Z u ´vah o ˇs´ıˇren´ı akustick´ ych vln si odvod´ıme okrajové podm´ınky na hranici ∂Ω = Γin ∪Γout ∪Γ0 , viz , viz [E. Bängtsson(2002)]. Pomoc´ı MKP provedeme ˇreˇsen´ı v systému SfePy. Na dané oblasti vyzkouˇs´ıme r˚ uzné frekvence ω incidenˇcn´ı vlny. V´ yslednou mapu ˇreˇsen´ı na oblasti si zobraz´ıme v programu ParaView. Ve páté kapitole formulujeme u ´lohu tvarové optimalizace pro akustické pole. Uvaˇzujeme oblast Ω a rozdˇel´ıme ji na dvˇe podoblasti, designovou ΩD a zbylou ΩC . MKP s´ıt’ oblasti vygenerujeme v MSC.Marc Mentat. Podél oblasti ΩD zavedeme hranici ΓD ⊂ Γ0 . ´ Uvaha byla taková, ˇze zmˇenou tvaru hranice ΓD , zprostˇredkovanou zmˇenou parametr˚ u α, se nám bude mˇenit ˇreˇsen´ı p v celé oblasti. Zmˇena se nám projev´ı i na u ´ˇcelové funkci Φ. Definujeme optimalizaˇcn´ı u ´lohu, kde jako stavovou rovnici vyuˇzijeme Helmholtzovu rovnici a jej´ı slabou formulaci. Provedeme citlivostn´ı anal´ yzu, viz [Rohan(2012)]. T´ım se rozum´ı v´ ypoˇcet citlivosti zmˇeny u ´ˇcelové funkce v závislosti na zmˇenˇe optimalizaˇcn´ıch parametr˚ u, kdyˇz na nich závis´ı nepˇr´ımo prostˇrednictv´ım stavové promˇenné. Zparametrizu1

´ Uvod jeme design pomoc´ı spline-boxu, viz [Rohan(2007)]. Zavedeme pole designov´ ych rychlost´ı. Optimalizaˇcn´ı u ´loze pˇriˇrad´ıme Lagrangeovu funkci a posléze pˇrejdeme k ˇreˇsen´ı adjungované u ´lohy. Touto cestou, poˇc´ıtanou v systému SfePy, z´ıskáme celkovou derivaci neboli citlivost u ´ˇcelové funkce Φ na zmˇenu parametr˚ u α. Následnˇe porovnáme pro kontrolu v´ ysledky citlivostn´ı anal´ yzy a koneˇcn´ ych diferenc´ı. Dalˇs´ı zdroje v této kapitole byly [E. Bängtsson(2002)],[B. Engquist(1977)],[Rohan(2007)],[D. Cioranescu(2008)], [J. Haslinger(1996)],[E. Rohan(2011b)],[E. Rohan(2006)],[E. Rohan(2011a)]. Na závˇer provedeme nˇekolik optimalizaˇcn´ıch v´ ypoˇct˚ u pro 2D oblast Ω. Optimalizován´ım tvaru hranice ΓD budeme hledat minimum dvou u ´ˇcelov´ ych funkc´ı ΦI,II . V´ ypoˇcet provedeme pomoc´ı Matlabu a SfePy. V´ ysledky zobraz´ıme v ParaView a zhodnot´ıme.

2

2 Odvozen´ı rovnic akustiky 2.1

´ Uvodem

V tekutém prostˇred´ı je zvuk vyvolán zmˇenami, které v ˇcase a prostoru zp˚ usob´ı zmˇeny tlaku, hustoty a kmitán´ı jednotliv´ ych ˇcástic prostˇred´ı s lokáln´ımi v´ ychylkami ξ~ = (x, y, z) a rychlostmi ~v = (vx , vy , vz ). O prostˇred´ı pˇredpokládáme, ˇze je stlaˇcitelné, spojité, homogenn´ı, izotropn´ı a neviskózn´ı. V´ ychylky ˇca´stic podle naˇseho pˇredpokladu jsou malé a vˇsechny sledované jevy budeme povaˇzovat za lineárn´ı. (Pˇri zpracován´ı této kapitoly se vycház´ı vesmˇes ˇ z knihy [Skvor(2001)] a z [Linhart(2009)].) V nepohybuj´ıc´ım se prostˇred´ı je akustické pole popsáno vektory akustick´ ych rychlost´ı jednotliv´ ych ˇca´stic prostˇred´ı. V prostˇred´ı s uvaˇzován´ım pohybu ˇca´stic je vyjádˇren´ı sloˇzitˇejˇs´ı o sloˇzku proudˇen´ı a pohyb popisujeme Lagrangeov´ ymi nebo Eulerov´ ymi promˇenn´ ymi. Zde nebudeme dále proudˇen´ı ˇcástic uvaˇzovat. Akustické pole pokládáme na nev´ırové, a tak se tedy pˇredpokládá, ˇze rot ~v = ~0.

(2.1)

V´ yj´ımku tvoˇr´ı pˇr´ıpady velk´ ych rychlost´ı u reáln´ ych plyn˚ u a kapalin, kdy se uplatˇ nuje vliv viskozity. Vektorové pole lze rozdˇelit na dvˇe sloˇzky, na sloˇzku nev´ırovou ~vn a v´ırovou (rotaˇcn´ı) ~vr . ~v = ~vn + ~vr (2.2) Pro nev´ırovou sloˇzku plat´ı rot~vn = 0 a div~vn 6= 0. Na základˇe toho zavedeme skalárn´ı fci φ [m2 s−1 ]. Nazveme ji rychlostn´ım potenciálem, jehoˇz gradient je roven rychlosti v~n [ms−1 ] grad φ = v~n .

(2.3)

Podle pˇredpokladu je pole rychlost´ı nev´ırové, a bude tedy platit ~v = v~n = grad φ div ~v = div v~n 6= 0 .

(2.4)

Pˇri odvozen´ı vlnové rovnice vyjdeme z pohybové (Eulerovy) rovnice , rovnice kontinuity (spojitosti) a stavové rovnice.

2.2

Eulerova rovnice

Eulerova rovnice je pohybová rovnice pro proudˇen´ı, která zanedbává tˇren´ı a vazkost. Odvozen´ı je provedeno pro elementárn´ı krychli dxdydz. Vyjdeme ze známého DÁlembertova 3

Odvozen´ı rovnic akustiky

Eulerova rovnice

principu, kter´ y ˇriká, ˇze v kaˇzdém okamˇziku je soustava vnitˇrn´ıch a vnˇejˇs´ıch sil v rovnováze X ~i, dm ~a = dF i

kde levá strana reprezentuje vnitˇrn´ı (setrvaˇcné) s´ıly a pravá vnˇejˇs´ı s´ıly p˚ usob´ıc´ı na element. −2 Pro zrychlen´ı ~a [ms ] plat´ı d~v dvx dvy dvz = , , , (2.5) ~a = ~i ax + ~j ay + ~k az = (ax , ay , az ) = dt dt dt dt kde ~i, ~j a ~k jsou jednotkové vektory orientovány v kladném smyslu ve smˇerech os x, y a z. Elementárn´ı krychle má hmotnost dm a tekutina v n´ı hustotu ρ [kg m−3 ] dm = ρdV = dxdydz . Setrvaˇcná s´ıla p˚ usob´ıc´ı na elementárn´ı krychli ve smˇerech os x, y a z je

Obrázek 2.1: S´ıly p˚ usob´ıc´ı na elementárn´ı krychli dvx , dF~x1 = ~i dm dt ~y1 = ~j dm dvy , dF dt dv ~z1 = ~k dm z . dF dt V´ ysledná setrvaˇcná s´ıla bude tedy rovna souˇctu jednotliv´ ych setrvaˇcn´ ych sil (2.6) ~ s = dF~x1 + dF ~y1 + dF ~z1 = dm ~i dvx + ~j dvy + ~k dvz = dm d~v . dF dt dt dt dt 4

(2.6)

(2.7)


Eulerova rovnice

V kartézsk´ ych souˇradnic´ıch je akustick´ y tlak p [P a] funkc´ı ˇcasu a prostorov´ ych promˇenn´ ych x, y a z. Takˇze p = p(t, x, y, z). Pˇrepokládejme podle obr.(2.2), ˇze na protilehl´ ych ploˇskách ∂p dy. Tlakové elementárn´ı krychle kolm´ ych na osu y p˚ usob´ı akustické tlaky p a proti p + ∂y

Obrázek 2.2: Tlak na elementárn´ı krychli s´ıly p˚ usob´ıc´ı na elemetárn´ı ploˇsku dxdz pˇredstavuj´ı u ´ˇcinek okoln´ı tekutiny nacházej´ıc´ı se vnˇe elementárn´ı krychle. V´ ysledná sloˇzka s´ıly ~j dFy1 ve smˇeru osy y bude rovna v´ yslednici povrchov´ ych sil p˚ usob´ıc´ıch na protilehl´ ych stˇenách h i ~j dFy2 = ~j p − p + ∂p dy dxdy ∂y ∂p ~y2 = −~j dydxdz . dF ∂y

(2.8)

Podobnˇe se odvod´ı ∂p dF~x2 = −~i dxdydz , ∂x

~z2 = −~k ∂p dzdydx . dF ∂z

(2.9)

~ v = −~i ∂p dxdydz − ~j ∂p dydxdz − ~k ∂p dzdydx . dF ∂x ∂y ∂z

(2.10)

~ v se rovná V´ ysledná s´ıla dF

~ v = dF ~ s . Dosad´ıme tedy z (2.10) a (2.7) Z DÁlembertova principu mus´ı platit dF ∂p ∂p ~ ∂p ∂~v ~ ~ − i +j +k =ρ , ∂x ∂y ∂z ∂t −grad p = ρ

∂~ v

. (2.11) ∂t Rovnice (2.11) je Eulerova rovnice dynamiky ideáln´ı tekutiny bez pˇrihl´ıˇzen´ı k jej´ımu proudˇen´ı a p˚ usoben´ı vnˇejˇs´ıch sil z jej´ıho okol´ı.

5


2.3

Rovnice kontinuity

Rovnice kontinuity

Dalˇs´ı z rovnic je rovnice kontinuity vyjadˇruj´ıc´ı zákon zachován´ı hmoty. Budeme uvaˇzovat podle obr.(4.4) elementárn´ı objem. Pˇredpokládejme, ˇze podél osy y v kladném smyslu vtéká a vytéká tekutina ploˇskami dx, dz. Hmotnost pˇritékaj´ıc´ı tekutiny za dt je ρ vy dx dz dt = ρ dwy dt .

(2.12)

kde wy [m3 s−1 ] je pr˚ utoková rychlost ve smˇeru osy y. Hmotnost vytékaj´ıc´ı tekutiny za dt je ∂(ρ dwy dt) dy . (2.13) ρ dwy dt + ∂y Je vidˇet, ˇze vyteˇce z elementárn´ı krychle o hmotnost

∂(ρ dwy dt) dy ∂y

v´ıce neˇz vteklo. Obdobnˇe

Obrázek 2.3: Pr˚ utok elementárn´ı krychl´ı tento rozd´ıl urˇc´ıme pro proudˇen´ı podél zbyl´ ych os x a z. Celkem tedy ∂(ρ dwx dt) dx , ∂x

∂(ρ dwy dt) dy , ∂y

∂(ρ dwz dt) dz . ∂z

(2.14)

Fakt, ˇze v´ıc hmotnosti odteˇce neˇz vteˇce, se mus´ı projevit poklesem hustoty v elementárn´ım ˇ dt a elementárn´ı objemy dxdydz se zkrát´ı objemu. Cas ∂(ρ dwx dt) ∂(ρ dwy dt) ∂(ρ dwz dt) ∂ρ dx + dy + dz = − dt dx dy dz ∂x ∂y ∂z ∂t ∂(ρ dvx ) ∂(ρ dvy ) ∂(ρ dvz ) ∂ρ + + = − ∂x ∂y ∂z ∂t ∂ρ div (ρ ~v ) = − ∂t 6

(2.15)


Stavová rovnice

V lineárn´ı akustice uvaˇzujeme v celém objemu tekutiny ve stejném ˇcase pˇribliˇznˇe stejnou hustotu ρ = ρ0 + ρ0 . Odchylky ρ0 od poˇca´teˇcn´ı hustoty ρ0 jsou velmi malé . Proto mohu celou rovnici vydˇelit ρ0 a z´ıskám t´ım rovnici kontinuity ve zjednoduˇseném tvaru div ~ v=−

2.4

1 ∂ρ ρ0 ∂t

.

(2.16)

Stavov´ a rovnice

Termodynamické chován´ı plynu je popsáno vzájemnou souvislost´ı tlaku p [P a], hustoty ρ [kg m−3 ] a teploty T [K]. Stavová rovnice ideáln´ıho plynu má tvar p=ρrT ,

(2.17)

kde r [Jkg −1 K −1 ] je mˇerná plynová konstanta. Jevy prob´ıhaj´ıc´ı v akustice jsou zpravidla tak rychlé, ˇze nedocház´ı ke sd´ılen´ı tepla, neboli dQ = 0, a dˇeje se stávaj´ı adiabatick´ ymi. Pˇri n´ızk´ ych kmitoˇctech jsou zmˇeny pomalé a termodynamické chován´ı plynu se bl´ıˇz´ı izotermickému dˇeji, tedy dˇeji, kter´ y prob´ıhá pˇri T = konst. V tˇechto pˇr´ıpadech je hustota plynu pouze funkc´ı tlaku ρ = ρ(p) a hovoˇr´ıme o baratropn´ım dˇeji. Pro izotermickou stavovu zmˇenu plat´ı podle Boyleova-Mariottova zákona p = konst . ρ

(2.18)

Pro adiabatickou stavovou zmˇenu plat´ı Poissonova rovnice pro tlak a objem p V κ = konst .

(2.19)

Objem je tedy pˇribliˇznˇe nepˇr´ımo u ´mˇern´ y hustotˇe V ∼ ρ−1 . Pak m˚ uˇzeme zapsat p ρ−κ = konst ,

(2.20)

kde konstanta κ je pomˇer mˇerné tepelné kapacity plynu pˇri stálém tlaku cp a stálém objemu cv cp κ= ,. (2.21) cv Podle adiabatického dˇeje popsaného v (2.20) jsou hustota ρ a tlak p v urˇcitém m´ıstˇe závislé na ˇcase, a tak lze stanovit závislost ˇcasové zmˇeny tlaku na zmˇenˇe hustoty ∂p −κ ∂ρ ρ − κ p ρ−κ−1 = 0 ∂t ∂t ∂p κ p ∂ρ = . ∂t ρ ∂t

7

(2.22)


Vlnová rovnice v kartézských souˇradnic´ıch

Kdyˇz budeme uvaˇzovat tlak p a hustotu ρ jako superpozici klidov´ ych hodnot p0 , ρ0 a zmˇen p0 , ρ 0 ρ = ρ 0 + ρ 0 , p = p0 + p0 . (2.23) A pˇredpokládáme p0 p0 , ρ0 ρ0 . Potom bude m´ıt rovnice (2.22) tvar ∂p ∂t

2.5

=

κ p0 ∂ρ ρ0 ∂t

.

(2.24)

Vlnov´ a rovnice v kart´ ezsk´ ych souˇ radnic´ıch

V pˇredchoz´ıch odstavc´ıch jsme odvodili tˇri rovnice, které vyuˇzijeme ke stanoven´ı vlnové rovnice. Jde o Eulerovo pohybovou rovnici (2.11), rovnici kontinuity(2.16) a stavovou rovnici (2.24) ∂~v 1 ∂ρ ∂p κ p0 ∂ρ − grad p = ρ , div ~v = − , = . ∂t ρ0 ∂t ∂t ρ0 ∂t Podle dˇr´ıvˇejˇs´ıho pˇredpokladu je rot~v = 0 a m´ısto rychlosti ~v zavedeme rychlostn´ı potenciál φ. Do rovnice (2.11) dosad´ıme rychlostn´ı potenciál a dostaneme ∂φ ∂ grad p = −ρ0 (grad φ) = −ρ0 grad (2.25) ∂t ∂t a dále po u ´pravˇe p = −ρ0

∂φ + konst . ∂t

(2.26)

V rovnici (2.16) také pˇrejdeme k rychlostn´ımu potenciálu a dostaneme div grad φ = −

1 ∂ρ . ρ0 ∂t

(2.27)

Zavedeme Laplaceo˚ uv operátor div grad φ = ∆ φ , pak rovnice dostane tvar ∆φ=−

1 ∂ρ . ρ0 ∂t

(2.28)

Rovnici (2.26) zderivujeme podle ˇcasu ∂p ∂ 2φ = −ρ0 2 ∂t ∂t

8

(2.29)


Vlnová rovnice v kartézských souˇradnic´ıch

a pak dáme v rovnost pravou stranu (2.29) s pravou stranou z (2.24) κ p0 ∂ρ ∂ 2φ = −ρ0 2 ρ0 ∂t ∂t ρ0 ∂ 2 φ 1 ∂ρ = − ρ0 ∂t κ p0 ∂t2

(2.30)

Do (2.30) dosad´ıme z (2.28) a koneˇcnˇe z´ıskáme vlnovou rovnici pro rychlostn´ı potenciál φ ∆φ=

ρ0 ∂ 2 φ κ p0 ∂t2

(2.31)

Dále vyuˇzijeme toho, ˇze známe vzorec pro v´ ypoˇcet rychlosti zvuku r κ p0 . c0 = ρ0 Rovnice (2.31) po dosazen´ı za rychlost zvuku bude vypadat takto ∆φ=

1 ∂ 2φ c20 ∂t2

,

(2.32)

nebo po rozepsán´ı Laplaceova operátoru do kartézsk´ ych souˇradnic na jednotlivé sloˇzky 1 ∂ 2φ ∂ 2φ ∂ 2φ ∂ 2φ + + = . ∂x2 ∂y 2 ∂z 2 c20 ∂t2

(2.33)

Uveden´ ym zp˚ usobem m˚ uˇzeme ze tˇr´ı základn´ıch rovnic stanovit téˇz vlnovou rovnici pro akustick´ y tlak. Vztah (2.11) zdivergujeme −div grad p = ρ0

∂ div ~v ∂t

a dosad´ıme z (2.16) div grad p =

∂ 2ρ ∂ 2ρ , ∆ p = . ∂t2 ∂t2

Z rovnice (2.24) urˇc´ıme ∂ 2p κ p0 ∂ 2 ρ = ∂t2 ρ0 ∂t2 ∂ 2 p ρ0 ∂ 2ρ = . ∂t2 κ p0 ∂t2 Po dosazen´ı z pˇredchoz´ı rovnice dostaneme vlnovou rovnici pro akustický tlak, kde jeˇstˇe pouˇzijeme opˇet vztahu pro rychlost zvuku ∆p=

1 ∂ 2p c20 ∂t2 9

.

(2.34)


2.6

Helmholtzova rovnice


Pˇri odvozen´ı Helmholtzovy rovnice (HR) uvaˇzujeme jednofrekvenˇcn´ı ˇcasovˇe harmonické ˇreˇsen´ı pro jednu frekvenci ω. Celkové ˇreˇsen´ı vlnové rovnice 2.32 pro tˇri prostorové promˇenné x, y, z a promˇennou ˇcas t budeme pˇredpokládat ve tvaru φ(x, y, z, t) = p(x, y, z) eiωt ,

(2.35)

To dosad´ıme do (2.20). Po zkrácen´ı eiωt dostaneme ω 2 p + c2 ∇2 p = 0 ∇2 p + κ2 p = 0 ,

(2.36)

kde κ = ωc . T´ım jsme z´ıskali Helmholtzovu rovnici, která popisuje ˇs´ıˇren´ı vlnˇen´ı ve tˇrech dimenz´ıch. Zápis pro HR s jednou prostorovou promˇennou x bude vypadat p00 (x) + κ2 p(x) = 0 .

10

(2.37)

ˇ sen´ı u 3 Reˇ ´ loh pro 1D kontinuum 3.1

Vlnov´ a rovnice

Vlnová rovnice (2.32) je parciáln´ı diferenciáln´ı rovnice hyperbolického typu, jej´ıˇz ˇreˇsen´ı popisuje ˇs´ıˇren´ı akustického potenciálu. Pro jednoduchost ji budeme ˇreˇsit v 1D pro jednu ˇ promˇennou x viz (3.1), viz [M´ıka(1983)] a [Skvor(2001)]. c

2

∂ 2φ ∂x2

−

∂ 2φ ∂t2

=0

(3.1)

Pˇri ˇreˇsen´ı rovnice (3.1) vyuˇzijeme Fourierovy metody neboli metody oddˇelen´ı (separace) promˇenn´ ych. Tato metoda nahrad´ı parciáln´ı derivace obyˇcejn´ ymi. U u ´loh s pevnˇe stanoven´ ymi okrajov´ ymi podm´ınkami, se hledá ˇreˇsen´ı ve tvaru souˇcinu v´ıce funkc´ı, kde kaˇzdá z nich je funkc´ı jedné nezávisle promˇenné. Tato metoda se pouˇz´ıvá hlavnˇe k ˇreˇsen´ı poˇca´teˇcnˇeokrajov´ ych u ´loh pro hyperbolické i parabolické rovnice. Zde budeme rychlostn´ı potenciál φ(x, t) pˇredpokládat jako souˇcin dvou funkc´ı, funkce X(x), která je funkc´ı jedné nezávisle ˇ sen´ı lze tedy promˇenné x, a funkce T (t), která je funkc´ı jedné nezávisle promˇenné t. Reˇ zapsat ve tvaru φ(x, t) = T (t) X(x) . (3.2) Po dosazen´ı do rovnice (3.1) z´ıskáme c2 ·

T¨ X 00 = . X T

(3.3)

Má-li platit rovnice (3.3) , tak mus´ı b´ yt obˇe jej´ı strany rovny jedné konstantˇe − λ = −κ2 , kde κ je vlnové ˇc´ıslo a uvaˇzujeme λ ≥ 0. Zavedli jsme zápornou konstantu, protoˇze kladná by vedla na hyperbolické funkce. Pˇredpokládáme, ˇze ani jedna z funkc´ı, nebo jej´ı druhá derivace, nen´ı v uvaˇzovaném oboru rovna nule. Z (3.3) z´ıskáme soustavu dvou obyˇcejn´ ych homogenn´ıch diferenciáln´ıch rovnic s konstantn´ımi koeficienty c2 X 00 + λX = 0 , T¨ + λT = 0 .

(3.4) (3.5)

Vyˇreˇsen´ım soustavy bychom mˇeli z´ıskat systém (φ1 (x, t), φ2 (x, t), ... , φk (x, t)) ˇreˇsen´ı p˚ uvodn´ı vlnové rovnice. Jako celkové ˇreˇsen´ı mus´ıme volit vhodnou kombinaci funkc´ı φ(x, t) =

∞ X

Dk φk (x, t) =

k=0

∞ X k=0

11

Dk Xk (x) Tk (t),

ˇ sen´ı u Reˇ ´loh pro 1D kontinuum

Vlnová rovnice

kde konstatn´ı koeficienty Dk se vol´ı tak, aby ˇreˇsen´ı φ(x, t) vyhovovalo okrajov´ ym i poˇcáteˇcn´ım podm´ınkám. Dále budeme ˇreˇsit soustavu diferenciáln´ıch rovnic jako vyˇsetˇrován´ı kmit˚ u mezi dvˇema dokonale odraziv´ ymi plochami. Jedna bude um´ıstˇena v poˇca´tku x = 0 a druhá ve vzdálenosti x = l. Vyjdeme z pˇredpokladu, ˇze akustická rychlost je na obou plochách nulová. Pro akustickou rychlost v plat´ı vztah (2.3). V 1D pro jedinou prostorovou promˇennou x postaˇcuje ∂ φ. Pouˇzijeme tyto okrajové podm´ınky v = ∂x ∂ φ(0, t) = X 0 (0) T (t) = 0 , ∂x ∂ φ(l, t) = X 0 (l) T (t) = 0 . ∂x

(3.6)

Jako poˇca´teˇcn´ı podm´ınky uvaˇzujeme ∂ φ(x, 0) = 0 , ∂t φ(x, 0) = g(x) ,

(3.7)

kde g(x) popisuje rozloˇzen´ı akustické rychlosti mezi dvˇema deskami v ˇcase t = 0. ˇ sen´ı povede k u Nyn´ı budeme ˇreˇsit (3.4) s okrajov´ ymi podm´ınkami. Reˇ ´loze na vlastn´ı ˇc´ısla a vlastn´ı funkce. c2 X 00 + λX = 0 , X 0 (0) = 0 , X 0 (l) = 0 .

(3.8)

X(x) = A eαx .

(3.9)

c2 α2 eαx = λeαx λ α2 = − 2 c 1√ α=± −λ c 1√ α = ±i λ. c

(3.10)

1√ 1√ = B · cos λx + C · sin λx . c c

(3.11)

Hledáme netriviáln´ı ˇreˇsen´ı ve tvaru

Po dosazen´ı (3.9) do (3.8) z´ıskáme

Po dosazen´ı α do (3.9) ±i 1c

X(x) = A · e

√

λx

12


Vlnová rovnice

Nyn´ı s pouˇzit´ım okrajov´ ych podm´ınek urˇc´ıme konstanty B a C. Najdeme prvn´ı derivaci X 0 (x) 1√ 1√ 1√ 1√ λx · λ + C · cos λx · λ. (3.12) X 0 (x) = −B · sin c c c c Ted’ dosad´ıme prvn´ı okrajovou podm´ınku z (3.6) 0 = C · cos(0) ·

1√ λ. c

(3.13)

Konstanta C mus´ı tedy b´ yt rovna nule. Dále vyuˇzijeme druhou okrajovou podm´ınku z (3.6) 1√ 1√ λl · λ. (3.14) 0 = −B · sin c c √ ˇ sen´ı pro B 6= 0 vyhovuje (3.14), jestliˇze sin( 1 λ l) = Z (3.14) vypl´ yvaj´ı dvˇe moˇznosti ˇreˇsen´ı. Reˇ c 0. Odtud plyne 1p λk l = kπ c (3.15) kπc 2 k∈N. λk = l Nyn´ı jiˇz známe systém vlastn´ıch fc´ı Xk (x) a vlastn´ıch ˇc´ısel λk πx kπc 2 Xk (x) = cos k , λk = ; k∈N. l l Vypoˇcteme si normu vlastn´ıch funkc´ı pro k > 0 a k = 0 Z l 2 cos kπ x dx = l l 2 0 w kπ w r l w w x w= , k>0 w cos l 2 Z l |cos(0)|2 dx = l 0 √ k cos(0)k = l , k = 0 Normovaná vlastn´ı funkce nabude tvaru r kπ 2 Xk (x) = cos x l l r 1 X0 (x) = k=0. l

(3.16)

(3.17)

k>0, (3.18)

Obecné ˇreˇsen´ı rovnice (3.4) má tvar X(x) =

∞ X

Bk Xk (x) ,

k=0

13

(3.19)


Vlnová rovnice

kde Bk jsou neznámé koeficienty. Pro známé koeficienty λk lze nyn´ı ˇreˇsit rovnici (3.5) s poˇca´teˇcn´ımi podm´ınkami (3.7) T¨ + λk T = 0 , ˙ 0) = 0 , φ(x,

(3.20)

φ(x, 0) = g(x) . Uvaˇzujeme jej´ı ˇreˇsen´ı ve tvaru Tk (t) = C eiωk t .

(3.21)

Dosad´ıme (3.21) do (3.20), odkud z´ıskáme ωk2 = λk p kπc ωk = ± λk = ± . l

(3.22)

Nyn´ı m˚ uˇzeme vyjádˇrit ˇreˇsen´ı pro k − tou funkci Tk Tk (t) = Dk cos(ωk t) + Ek sin(ωk t) , kde Dk a Ek jsou neznáme konstanty. Potom celkové ˇreˇsen´ı φ(x, t) =

∞ X

Bk Xk (x) Tk (t)

k=0

φ(x, t) =

∞ X

Bk [Dk cos(ωk t) + Ek sin(ωk t)] Xk (x) .

k=0

Dále pˇreznaˇc´ıme konstanty a budeme uvaˇzovat Dk = Bk · Dk a Ek = Bk · Ek . Vyjádˇr´ıme prvn´ı derivaci ˙ t) = φ(x,

∞ X

[−Dk ωk sin(ωk t) + Ek ωk cos(ωk t)] Xk (x) .

k=0

ˇ sen´ı mus´ı splˇ Reˇ novat poˇca´teˇcn´ı podm´ınky z (3.20). Pouˇzijeme prvn´ı z nich ˙ 0) = φ(x,

∞ X (0 + Ek ωk ) Xk (x) = 0 . k=0

Z toho vypl´ yvá, ˇze Ek = 0, ∀k ∈ N, a tedy plat´ı φ(x, t) =

∞ X

Dk Tk (t) Xk (x) = 0 ,

k=0

14

(3.23)


Vlnová rovnice

kde zb´ yvá urˇcit koeficienty Dk . U druhé podm´ınky z (3.20) zvol´ıme fci g(x) jako nekoneˇcnou ˇradu vlastn´ıch fc´ı, kde gk jsou konstantn´ı koeficienty vlastn´ıch funkc´ı Xk (x) g(x) =

∞ X

gk Xk (x) .

k=0

Celou rovnost vynásob´ıme j-tou vlastni funkc´ı a zintegrujeme v intervalu [0, l]. Z lX Z l ∞ g(x) Xj (x) dx = gk Xk (x) Xj (x) dx 0

0 k=0

=

∞ X

Z

Xk (x) Xj (x) dx .

gk

k=0

(3.24)

l

0

V´ıme, ˇze systém vlastn´ıch funkc´ı je ortonormáln´ı. Pro skalárn´ı souˇcin j-té a k-té vlastn´ı funkce potom plat´ı, kdyˇz j 6= k, ˇze Z l Xk (x) Xj (x) dx = 0 0

a kdyˇz j = k l

Z

Xk (x) Xk (x) dx = 1 . 0

Na pravé stranˇe (3.24) zbude jedin´ y ˇclen ˇrady, kdy k = j, a dostaneme pˇredpis pro koeficinty gj ˇrady zastupuj´ıc´ı funkci g(x) Z l g(x) Xj (x) dx . gj = 0

Nyn´ı m˚ uˇzeme g(x) zapsat g(x) =

∞ Z X k=0

l

g(x) Xk (x)dx Xk (x) .

(3.25)

0

Aplikujeme druhou poˇca´teˇcn´ı podm´ınku z (3.20) na (3.23) ∞ X k=0

gk Xk (x) =

∞ X

Dk cos(0) Xk (x)

k=0

gk Xk (x) = Dk Xk (x) r Z l Z l 2 D k = gk = g(x) Xk (x)dx = g(x) cos (κx) dx , l 0 0 r Z l Z l 1 D0 = g(x) X0 (x)dx = g(x) dx , l 0 0 15

(3.26)


kde

q

l 2

Vlnová rovnice

je norma vlastn´ı fce. Celkové ˇreˇsen´ı rovnice (3.1) bude tedy m´ıt tento tvar # 2 φ(x, t) = Dk cos(ωk t) · cos (κx) l k=0 ∞ X 2 kπc kπ φ(x, t) = gk cos t cos x . l l l k=0 ∞ X

"

r

(3.27)

Pokud chceme nyn´ı volit funkci g(x), tak mus´ı splˇ novat okrajové podm´ınky (3.6). Kdyˇz 3π si zvol´ıme tˇreba g(x) = cos( l x), kde v´ıme, ˇze jde vlastn´ı funkci pro k = 3. Z kolmosti vlastn´ıch funkc´ı Xk (x) a z (3.26) vypl´ yvá, ˇze jen pro k = 3 bude Dk nenulové 2 Z lr 2 3π D3 = cos x dx (3.28) l l 0 r r 2 l 2 = . D3 = l 2 l Celkové ˇreˇsen´ı tedy nabude tvaru φ(x, t) = D3 cos

kπc t l

cos

kπ x l

.

(3.29)

ˇ sen´ı si zobraz´ıme pomoc´ı Matlabu viz obr.(3.1) pro následuj´ıc´ı hodnoty parametr˚ Reˇ u l = 10 [m] , c = 100 [m/s] . (3.30)

ˇ sen´ı vlnové rovnice pro 1D pˇr´ıpad Obrázek 3.1: Reˇ 16



Z obr.(3.1) je vidˇet, jak zvolená funkce g(x) splˇ nuje z okrajov´ ych podm´ınek (3.6) nulové derivace.

3.2


Uvaˇzujeme jednofrekvenˇcn´ı ˇcasovˇe harmonické ˇreˇsen´ı pro u ´hlovou frekvenci ω. Jiˇz dˇr´ıve odvozenou HR (2.36) budeme ˇreˇsit pro jednu prostorovou neznámou x. HR nabude tvaru ω 2 p(x) + c2 p00 (x) = 0 p00 (x) + κ2 p(x) = 0 .

(3.31)

ˇ s´ıme opˇet 1D u Reˇ ´lohu, kterou si lze pˇredstavit jako vlnˇen´ı mezi dvˇema deskami. Prvn´ı deska v x = 0 bude kmitat frekvenc´ı ω s amplitudou pA . P˚ ujde tedy o pˇr´ıpad vibroakustiky. Druhá deska bude nehybnou dokonale odrazivou pˇrekáˇzkou, takˇze akustická rychlost zde mus´ı b´ yt nulová. p0 (l) = 0 , p0 (0) = iωpA .

(3.32)

Dále si p(x) zvol´ıme jako souˇcet dvou funkc´ı p(x) = q(x) + q¯(x) .

(3.33)

Pro funkce q(x) a q¯(x) vol´ıme okrajové podm´ınky ve tvaru q¯0 (0) = iωpA , q¯0 (l) = 0 , q 0 (0) = 0 , q 0 (l) = 0 . Funkci q¯(x) zvol´ıme tak, aby splˇ novala v´ yˇse uvedené okrajové podm´ınky π x . q¯(x) = iωpA sin 2l

(3.34)

(3.35)

A funkci q(x) budume uvaˇzovat jako lineárn´ı kombinaci k vlastn´ıch funkc´ı urˇcen´ ych Fourierovou metodou (3.18) s konstatn´ımi koeficienty Hk q(x) =

∞ X

Hk Xk (x) .

(3.36)

k=0

Kdyˇz dosad´ıme za p(x) fuknkce q(x) a q¯(x), tak rovnice (3.31) dostane tvar ω 2 q(x) + ω 2 q¯(x) + c2 q 00 (x) + c2 q¯00 (x) = 0 ω 2 q(x) + c2 q(x)00 = −(ω 2 q¯(x) + c2 q¯00 (x)) . 17

(3.37)



Pravou stranu (3.31) oznaˇc´ıme f (x) a zvol´ıme zase jako nekoneˇcnou ˇradu vlastn´ıch funkc´ı s konstantami fk ∞ X 2 2 00 f (x) = −(ω q¯(x) + c q¯ (x)) = fk Xk (x) . (3.38) k=0

Koeficienty fk budeme uvaˇzovat ve tvaru Z l fk = f (x) Xk (x) dx .

(3.39)

0

Po dosazen´ı za pravou stranu m˚ uˇzeme celou rovnici vynásobit j-tou vlastn´ı fc´ı a integrovat pˇres interval [0, l] ω 2 q(x) + c2 q 00 (x) =

∞ X

fk Xk (x)

(3.40)

fk Xk (x) Xj (x) dx .

(3.41)

k=0

ω2

l

Z

q(x) Xj (x) dx + c2

0

Z

l

q 00 (x) Xj (x) dx =

∞ Z X

0

k=0

l

0

V´ıme, ˇze vlastn´ı fce jsou ortonormáln´ı, takˇze integrál ze souˇcinu dvou vlastn´ıch fc´ı, k-té a j-té, je roven nule, pokud j 6= k. Pokud plat´ı j = k, skalárn´ı souˇcin je roven jedné. Na pravé stranˇe zbude tedy fj . Na levé dosad´ıme za q(x) a q 00 (x). ω

2

Z lX ∞

Hk Xk (x) Xj (x) dx + c

2

Z lX ∞

00

Hk Xk (x) Xj (x) dx = fj

(3.42)

0 k=0

0 k=0

V´ıme, ˇze

00

c2 Xk (x) = (−λk ) Xk (x) ,

(3.43)

takˇze po dosazen´ı (3.42) do (3.43) ω

2

Z lX ∞

Hk Xk (x) Xj (x) dx +

Z lX ∞

Hk (−λk ) Xk (x) Xj (x) dx = fj

0 k=0

0 k=0

(3.44)

Hj (ω 2 − λj ) = fj . Vyjádˇr´ıme konstanty Hj (dále uˇz Hk ) fk = Hk = 2 ω − λk

Rl

f (x) Xk dx 0 = ω 2 − λk

Rl 0

(−ω 2 q¯(x) − c2 q¯00 (x)) Xk dx . ω 2 − λk

(3.45)

N´ıˇze urˇc´ıme druhou derivaci funkce q¯(x) π 2l q¯(x) = iωpA sin x π 2l π π q¯00 (x) = −iωpA sin x . 2l 2l 18

(3.46)



Dopoˇcteme koeficienty fk Z l 2ω 2 l c2 π π kπ r 2 (−iωpA )· − · sin x · cos x · dx fk = π 2l 2l l l 0 r Z l r π kπ 2 2 sin fk = −µ · x · cos x dx = −µ ·I l 0 2l l l (3.47) r 2l 2 fk = −µ k>0 l π· (1 − 4k 2 ) r 1 2l f0 = −µ k=0, l π kde µ je konstanta 2 2ω l c2 π − . (3.48) iωpA π 2l V´ ypoˇcet integrálu I jsme provedli metodou Per partes Z l π kπ sin I = x · cos x dx 2l l 0 π kπ π kπ il h l l sin x · sin x + cos x · cos x + = πk 2l l 2πk 2 2l l 0 Z l π kπ 1 + 2 sin x · cos x dx 4k 0 2l l 1 −l I · 1− 2 = 4k 2πk 2 2l I = k ∈N. (3.49) π(1 − 4k 2 ) Konstatn´ı koeficienty Hk maj´ı tedy tvar Hk =

ω2

fk k>0 − ( kπc )2 l

f0 H0 = 2 k = 0 . ω Nyn´ı jiˇz známe vˇse a m˚ uˇzeme vyjádˇrit hledané ˇreˇsen´ı HR (3.31)

(3.50)

p(x) = q(x) + q¯(x) ∞ π X p(x) = Hk Xk (x) + iωpA sin x . 2l k=0 Na obr.(3.2) je zobrazeno ˇreˇsen´ı. Pouˇzili jsme tyto hodnoty parametr˚ u l = 10 [m] , ω = 750 [rad/s] , c = 342 [m/s] , pA = 20 [m/rad] . 19

(3.51)



ˇ sen´ı Hemholtzovy rovnice pro 1D pˇr´ıpad Obrázek 3.2: Reˇ

Je vidˇet, ˇze vˇsechny reálné ˇca´sti vycház´ı nulové. Je to zp˚ usobeno t´ım, ˇze jsme okrajovou podm´ınku z (3.32)2 , neboli buzen´ı, zvolili imaginárn´ı.

20

4 Formulace u´loh akustiky pro danou frekvenci ve 3D Pokud chceme ˇreˇsit vlnovou rovnici pro akustick´ y tlak (2.34) pro jednu danou frekvenci ω (jednofrekvenˇcn´ı ˇreˇsen´ı), mus´ıme vlastnˇe ˇreˇsit Helmholtzovu rovnici (2.36) s vhodn´ ymi okrajov´ ymi podm´ınkami jako okrajovou u ´lohu.

4.1

Okrajov´ au ´ loha

K odvozen´ı okrajov´ ych podm´ınek vyuˇzijeme oblast Ω (obr.(4.1)). Hranici oznaˇc´ıme ∂Ω. Plat´ı ∂Ω = Γin ∪ Γout ∪ Γ0 , kde Γin je hranice vstupu incidenˇcn´ı vlny. V´ ystup z Ω pˇredstavuje Γout . Zb´ yvaj´ıc´ı hranici oznaˇc´ıme Γ0 . P˚ ujde o pevnou dokonale odrazivou pˇrekáˇzkou. ~n je vnˇejˇs´ı jednotkov´ y normálov´ y vektor hranice ∂Ω. (Zde vyjdeme hlavnˇe z [E. Bängtsson(2002)].)

Obrázek 4.1: Obecná oblast Ω se vstupn´ı hranic´ı Γin , v´ ystupn´ı Γout a dokonale odrazivou Γ0

V praktick´ ych v´ ypoˇctech je ˇcasto dobré pouˇz´ıt umˇelé okrajové podm´ınky, které omez´ı oblast v´ ypoˇctu. Na Γin , kde vniká incidenˇcn´ı vlna do oblasti Ω, by mˇela b´ yt urˇcena jej´ı amplituda, zat´ımco odchoz´ı, odraˇzená vlna, by mˇela z˚ ustat nepoznamenána. Necht’ ~x popisuje prostorové souˇradnice. Pˇredpokládáme, ˇze jednofrekvenˇcn´ı rovinné vlny, procházej´ıc´ı Γin , m˚ uˇzeme zapsat ve tvaru P (x, t) = A ei(κ ~x·~n+ωt) + B ei(−κ ~x·~n+ωt)

, A, B ∈ C ,

(4.1)

kde prvn´ı ˇclen odpov´ıdá incidenˇcn´ı vlnˇe a druh´ y vlnˇe odraˇzené. Skalárn´ı souˇcin vektor˚ u

21

Form. u ´loh akust. pro danou frekv. ve 3D

Slabé ˇreˇsen´ı

~x · ~n = n je vlastnˇe pr˚ umˇet ~x do ~n. Derivován´ım rovnice (4.1) dostaneme ∂P ∂t ∂P ∂~n

= Aiω ei(κ ~x·~n+ωt) + Biω ei(−κ ~x·~n+ωt) ,

(4.2)

= Aiκ ei(κ ~x·~n+ωt) − Biκ ei(−κ ~x·~n+ωt) ,

(4.3)

kde ∂P = ~n · gradP popisuje derivaci ve smˇeru vnˇejˇs´ı normály. Seˇcten´ım (4.2) a (4.3) ∂n vynásobenou c a vyuˇzit´ım ω = κ c z´ıskáme okrajovou podm´ınku na Γin ∂P ∂P +c = 2iωA ei(κ ~x·~n+ωt) . ∂t ∂n

(4.4)

Podm´ınka (4.4) je splnˇena pro kaˇzdou vlnu typu (4.1) a m˚ uˇze tedy b´ yt pouˇzita k nastaven´ı amplitudy A incidenˇcn´ı vlny bez vlivu na amplitudu B vlny odraˇzené. Vhodná absorbˇcn´ı okrajová podm´ınka na Γout je radiaˇcn´ı podm´ınka prvn´ıho stupnˇe [B. Engquist(1977)]. Nebo m˚ uˇzeme uvaˇzovat, ˇze na v´ ystupu nen´ı ˇzádná incidenˇcn´ı vlna a tedy A = 0 ∂P ∂P +c =0. ∂t ∂n

(4.5)

Tato podm´ınka nezp˚ usob´ı odraz vln na hranici Γout jdouc´ıch pˇr´ımo pˇresnˇe ve smˇeru normály. Na zbyl´ ych hranic´ıch Γ0 pˇredpokládáme akustick´ y tvrd´ y materiál, takˇze budeme uvaˇzovat podm´ınku ∂P =0. ∂n

(4.6)

Pokud chceme ˇreˇsit u ´lohu pro jednu frekvenci, uvaˇzujeme ˇreˇsen´ı (2.35) jako v kapitole, kde jsme Helmhotzovu rovnici odvodili. Dosad´ıme do okrajov´ ych podm´ınek za P . Po zkrácen´ı iωt e a za pˇredpokladu, ˇze na Γin je v ~x = 0 z´ıskáme okrajové podm´ınky ∂p + iωp = 0 na Γout , ∂n ∂p c + iωp = 2iωA na Γin , ∂n ∂p = 0 na Γ0 . ∂n

c

4.2

(4.7)

Slab´ eˇ reˇ sen´ı

Pˇri odvozen´ı slabého ˇreˇsen´ı, viz [M´ıka(2007)], vyjdeme z HR (2.36) pro akustick´ y tlak p = p(x). Kde κ je vlnové ˇc´ıslo, ω je u ´hlová frekvence a c je rychlost ˇs´ıˇren´ı vlny prostˇred´ım (v tomhle pˇr´ıpadˇe jde o rychlost ve vzduchu), pro které plat´ı κ = ωc . HR celou vynásob´ıme 22


Numerické ˇreˇsen´ı MKP

testovac´ım tlakem q, kde q ∈ Q a Qje mnoˇzina vˇsech pˇr´ıpustn´ ych tlak˚ u, které vyhovuj´ı okrajov´ ym podm´ınkám (4.7) O2 p q + κ2 p q = 0

∀q ∈ Q . (4.8)

Celou rovnici integrujeme pˇres oblast Ω Z Z 2 O p q dΩ + κ2 p q dΩ = 0 Ω

∀q ∈ Q .

(4.9)

Ω

Vyuˇzijeme následuj´ıc´ı Greenovy vˇety Z Z Z ∂p 2 O p q dΩ = q d∂Ω − OqOp dΩ . ∂n Ω ∂Ω Ω Dosad´ıme (4.10) do (4.9) a dostaneme Z Z Z ∂p 2 q κ p q dΩ + d∂Ω − OqOp dΩ = 0 ∂n ∂Ω Ω Ω

∀q ∈ Q .

(4.10)

(4.11)

Po dosazen´ı okrajov´ ych podm´ınek (4.7) do (4.11) dostaneme tento tvar slabého ˇreˇsen´ı Z Z Z Z 2 iκpq dΓout = 0 ∀q ∈ Q(4.12) . (2iκ˜ u − iκp)q dΓin + OqOp dΩ − κ pq dΩ − Ω

4.3

Ω

Γout

Γin

Numerick´ eˇ reˇ sen´ı MKP s r˚ uzn´ ymi okrajov´ ymi podm´ınkami

V´ ypoˇcet budeme provádˇet v oblasti Ω (obr.(4.2)), kde ∂Ω bude jej´ı hranice. Plat´ı ∂Ω = Γin ∪ Γout ∪ Γ0 , kde Γin je hranice vstupu (ˇcervená linie) , kde vzniká incidenˇcn´ı vlna. Γout pˇredstavuje v´ ystup z Ω (modrá linie) a Γ0 je dokonalé odrazivá stˇena.

23



Obrázek 4.2: Popsaná pouˇzitá oblast Ω

Obrázek 4.3: Detail MKP s´ıtˇe oblasti Ω

Na obr.(4.3) vid´ıme pouˇzitou s´ıt’ MKP. Jako elementy jsme pouˇzili trojuheln´ıky. Celkem s´ıt’ obsahuje 18724 uzl˚ u a 36556 elemet˚ u, a tud´ıˇz má 56172 stupˇ n˚ u volnosti. Bud´ıc´ı funkci oznaˇc´ım jako u˜. Pouˇzijeme následuj´ıc´ı okrajovou podm´ınku na Γin . Budeme uvaˇzovat buzen´ı konstantou u˜ = uA .

(4.13)

Buzen´ı funkc´ı (4.13) interpretujeme jako kmitán´ı desky vˇe smˇeru osy x s amplitudou uA a u ´hlovou frekvenc´ı ω. Parametry u ´lohy jsme si navolili c = 342 [m/s] , ω = 500, 1000 [rad/s] , ω κ = [rad/m] , c uA = 50 000 [−] , l = 10 [m] , a = 4 [m] , b = 2 [m] ,

(4.14)

kde l je délka hrany oblasti Ω, b ˇs´ıˇrka zvukovod˚ u, a délka zvukovod˚ u, c rychlost ˇs´ıˇren´ı zvuku ve vzduchu, ω u ´hlová frekvence a κ vlnové ˇc´ıslo. Pro u ´hlovou frekvenci ω jsme volili dvˇe hodnoty. V´ ypoˇcet jsme provádˇeli v systému SfePy a v´ ysledky zobrazili v programu ParaView. S´ıt’ MKP (obr.(4.3)) jsme vygenerovali v MSC.Marc/Mentat. Pro buzen´ı jsme provedli v´ ypoˇcty pro dvˇe r˚ uzné hodnoty ω viz (4.14).

24


Obrázek 4.4: Imaginárn´ı ω = 500 [rad/s], u˜ = uA

ˇca´st


p; Obrázek 4.5: Reálná ω = 500 [rad/s], u˜ = uA

ˇca´st

p;

Obrázek 4.6: Graf znázorˇ nuj´ıc´ı hodnoty podél pˇr´ımky z A do B. ω = 500 [rad/s], u˜ = uA

Vyjdeme-li pro kontrolu ze vzorc˚ u T = λ/c a T = 2π/ω, kde λ je vlnová délka a T perioda, z´ıskáme vzorec pro λ: λ = 2πc/ω. Po dosazen´ı ω = 500 [rad/s] a c = 342 [m/s] zjist´ıme, ˇze λ = 4, 298 [m], coˇz odpov´ıdá obr.(4.4).

25


Obrázek 4.7: Imaginárn´ı ω = 1000 [rad/s], u˜ = uA

ˇca´st


p; Obrázek 4.8: Reálná ω = 1000 [rad/s], u˜ = uA

ˇca´st

p;

Obrázek 4.9: Graf znázorˇ nuj´ıc´ı hodnoty podél pˇr´ımky z bodu A do bodu B. ω = 1000 [rad/s], u˜ = uA

Na obr.(4.9) je uˇz lépe vidˇet, jak se vlivem kratˇs´ı vlnové délky mˇen´ı hladkost zobrazen´ı gradient˚ u. Nedokonalost je zp˚ usobena nedostateˇcnou hustotou s´ıtˇe oblasti.

26

5 Formulace u´loh tvarové optimalizace 5.1

´ Uvod

V této kapitole formulujeme u ´lohu tvarové optimalizace pro akustické pole. Oblast Ω (obr.(5.1)) je rozdˇelena na dvˇe disjunktn´ı podoblasti: plat´ı Ω = ΩD ∪ ΩC , kde ΩD pˇredstavuje takzvanou designovou oblast, a pro ΩC plat´ı: ΩC = Ω \ ΩD . Hranice ∂Ω je rozdˇelena na d´ılˇc´ı hranice ∂Ω = Γ0 ∪ Γin ∪ Γout . Na pravé hranici je designová hranice ΓD ⊂ Γ0 ⊂ ∂Ω. Tato designová hranice bude optimalizac´ı mˇenit tvar a to následnˇe ovlivn´ı s´ıt’ MKP v ΩD , ale uˇz ne v ΩC . Zmˇenou tvaru hranice ΓD , zprostˇredkovanou zmˇenou parametr˚ u α, se bude mˇenit ˇreˇsen´ı p v celé oblasti. Zmˇena se projev´ı také na u ´ˇcelové funkci Φ, která se vyhodnocuje v ΩC .

Obrázek 5.1: Popsaná oblast Ω

Jako stavovou rovnici budeme uvaˇzovat Hemholtzovu rovnici (2.36) pro akustick´ y tlak ´ p. Ulohu budeme ˇreˇsit na Ω a pouˇzijeme na jej´ı hranici ∂Ω jiˇz odvozené okrajové podm´ınky (4.7). Dále zavedeme slabou formulaci u ´lohy (viz (4.12)) Pomoc´ı bilineárn´ı formy zap´ıˇseme naˇsi stavovou u ´lohu jako (∇p, ∇q)Ω − κ2 (p, q)Ω + hiκp, qiΓin ∪Γout − h2iκ˜ u, qi = 0 , | {z } | {z Γin} Ψ(p,q)

(5.1)

f (q)

kterou m˚ uˇzeme zapsat v obstraktn´ı formˇe Ψ(p, q) = f (q) 27

p, ∀q ∈ Q .

(5.2)

Formulace u ´loh tvarové optimalizace

Citlivostn´ı analýza

Uvaˇzujeme následuj´ıc´ı optimalizaˇcn´ı u ´lohu min Φ(p) , α

p ∈ Q : Ψ(p, q) = f (q) ∀q ∈ Q , α∈U , γ : α −→ ΓD ,

(5.3)

kde p splˇ nuje stavovou u ´lohu (5.2), Φ je u ´ˇcelová funkce a U je mnoˇzina pˇr´ıpustn´ ych parametr˚ u α. Zvol´ıme u ´ˇcelovou funkci Z |p|2 = kpk2Ωc [P a2 m2 ] , ΦI (p) = Ωc ! (5.4) kpkΩA [dB] . ΦII (p) = 20 log kpkΩB

5.2

Citlivostn´ı anal´ yza

Proveden´ım citlivostn´ı anal´ yzy rozum´ıme v´ ypoˇcet citlivosti zmˇeny u ´ˇcelové funkce Φ v závislosti na zmˇenˇe optimalizaˇcn´ıch parametr˚ u α, kdyˇz na nich závis´ı u ´ˇcelová funkce nepˇr´ımo prostˇrednictv´ım stavové promˇenné p. Potˇrebujeme tedy vypoˇc´ıtat gradient funkce Φ(α, p(α)) podle jednotliv´ ych optimalizaˇcn´ıch parametr˚ u α. Máme zde tedy implicitn´ı funkci p = p(α). V oblasti Ω definujeme vektorové pole V, s pomoc´ı Spline-boxu [Rohan(2007)].

5.2.1

Metoda adjungovan´ e promˇ enn´ e

Zab´ yvejme se nyn´ı u ´lohou, kdy tvar hranice závis´ı na jediném parametru τ ∈ R Γ(τ ) = ΓD + τ {V(x)}x∈ΓD ,

(5.5)

kde x je prostorová souˇradnice. Mus´ı platit, ˇze V(x) = 0 pro x ∈ ΩC a pro body na linii zv´ yraznˇené ˇcernou pˇreruˇsovanou ˇca´rou (viz obr.(5.1)). K problému (5.3)2 , v nˇemˇz je α nahrazeno parametrem τ , m˚ uˇzeme pˇriˇradit Lagrangeovu funkci L Lτ (p, q) = Φ(p) + Ψτ (p, q) − f (q) + Ψ∗τ (p∗ , q ∗ ) − f ∗ (q ∗ ) ,

(5.6)

kde q je Lagrangeov˚ uv multiplikátor vazby splnˇen´ı stavové u ´lohy. Poznamenejme, ˇze po∗ ∗ kud plat´ı (5.3)2 pro Ψ(p, q) = f (q), tak i pro Ψ (p , q) = f ∗ (q) , kde * znaˇc´ı komplexnˇe 28



sdruˇzenou hodnotu. Nyn´ı formulujeme u ´lohu (5.3) jako u ´lohu sedlového bodu max q∈Q

min L(α, p, q) .

(5.7)

α∈U,p∈Q

Rovnici (5.6) zderivujeme podle ˇreˇsen´ı p. Gâteauxovu derivaci, viz[Drábek(1994)], podle p oznaˇc´ıme δp (derivace ve smˇeru q) δp Lτ (p, q; δp) = δΦ(p; δp) + Ψτ (δp, q) + Ψ∗τ (δp∗ , q ∗ ) . | {z }

(5.8)

2Re{Ψτ (δp,q)}

Pro splnˇen´ı podm´ınky optimality se mus´ı (5.8) rovnat nule. Dostaneme se k ˇreˇsen´ı adjungované u ´lohy pro adjungovanou promˇennou q ∈ Q 2Re {Ψτ (δp, q)} δp ∈ Q = −δΦ(p; δp) ∀ .

(5.9)

ˇ sen´ı Adjungovanou u ´lohu m˚ uˇzeme ˇreˇsit za pˇredpokladu, ˇze Φ ∈ R, pro p ∈ H 1 (Ω; C). Reˇ ´ celovou funkci adjungované u ´lohy je efektivn´ı ˇreˇsen´ı nalezen´ı gradientu u ´ˇcelové funkce. Uˇ ´ celovou budeme vyhodnocovat na ΩC , tedy oblast´ı, která nen´ı ovlivnˇena zmˇenou ΓD . Uˇ funkci zderivujeme podle ˇreˇsen´ı p Z Z ∗ ∗ ∗ (p δp) . [p(δp) + p δp] = 2Re δΦ(p; δp) = ΩC

ΩC

Proto m˚ uˇzeme (5.9) pˇrepsat ve tvaru (∇p, ∇q)Ω − κ2 (p, q)Ω + hiκp, qiΓin ∪Γout + + ∇p∗ , ∇q ∗ )Ω − κ2 (p∗ , q ∗ )Ω − hiκp∗ , q ∗ iΓin ∪Γout = 2Re {Ψ(p, q)} . Zap´ıˇseme adjungovanou u ´lohu takto pro u ´ˇcelovou funkci Φ definovanou dle ΦI , viz (5.4), Z Z 1 ∗ p q =− Re {Ψ(ˆ p, q)} = −Re (pq ∗ + p∗ q) ∀q ∈ Q , (5.10) 2 ΩC ΩC kde pˆ oznaˇcuje odjungovanou promˇennou. Vypoˇcteme-li adjungovanou promˇennou, vyuˇzijeme ji dále v (5.6) jako Lagrange˚ uv multiplikátor a provedeme totáln´ı drivaci rovnice (5.6) podle τ p∗ )) .(5.11) δτ L(p, pˆ) + δp L(p, pˆ; δp) = δτ Φ(p) + δτ (Ψτ (p, pˆ) − f (ˆ p)) + δτ (Ψ∗τ (p∗ , pˆ∗ ) − f ∗ (ˆ | {z } =0

Je-li splnˇena adjungovaná rovnice (5.9), tak druh´ y ˇclen na levé stranˇe (5.11) je roven nule a plat´ı δΦ = δΨ.

29


5.2.2


V´ ypoˇ cet citlivosti pomoc´ı materi´ alov´ e derivace

V tomto odstavci vyjádˇr´ıme ˇcleny v citlivostn´ım vztahu (5.11) s vyuˇzit´ım materiálové derivace známé z mechaniky kontinua, viz [Rohan(2012)],[E. Rohan(2012)],[E. Rohan(2010a)], [E. Rohan(2010b)]. Zavedeme zobrazen´ı F parametrizované pomoc´ı τ F(τ, ·) : Ω 7→ Ω(τ ), z(τ ) = F(τ, x) ,

(5.12)

z = x + τ V(x), x ∈ Ω , z ∈ Ω(τ ) , Ω(τ ) = Ω + τ {V(x)}x∈Ω , J = det(∂x F) ,

(5.13)

kde

kde ∂x F je Jacobiova matice zobrazen´ı F. Plat´ı, viz [Rohan(2012)] Z (∇z p, ∇z q)Ω(τ ) = ∂iz xk ∂kx p ∂jz xl ∂lx qδij J(F) , ZΩ (p, q)Ω(τ ) = pqJ((F ))

(5.14)

Ω

a odtud plyne Z d (∇z p, ∇z q)Ω(τ ) = (div(∇V)δij − ∂k Vj δki − δej ∂e Vi )∂j p∂i q ∗ , dτ |τ =0 ZΩ d (p, q)Ω(τ ) = pq div(∇V) . dτ |τ =0 Ω Nyn´ı oznaˇc´ıme Z a(p, q) = (div(∇V)δij − ∂k Vj δki − δej ∂e Vi )∂j p∂i q ∗ , ZΩ m(p, q) = pq div(∇V) .

(5.15)

(5.16)

Ω

Vyuˇzili jsme tˇechto skuteˇcnost´ı · ∂zi (∂j Fi ) = = (∂j Vi ) = (∇V)ij ∂xj ∂xk ∂(zk − τ Vk ) ∂Vk = = δki − τ , ∂zi zi ∂zi ∂Fi Fij = , ∂xj F F −1 = I F˙ F −1 + F (F −1 )· = 0 ·

−1

F˙ F −1

) = −F ) |τ =0 = −F˙ |τ =0 = −∇V ,

(F (F

−1 ·

−1 ·

30

(5.17)

(5.18)


Implementace

∂Vj (5.19) . ∂xk S vyuˇzit´ım vztah˚ u (5.16) jiˇz m˚ uˇzeme dosadit do (5.11). Dostaneme totáln´ı derivaci δΦ vyjadˇruj´ıc´ı derivaci u ´ˇcelové funkce ve smˇeru vektorového pole V ∂k Vj =

δΦ = a(p∗ , pˆ∗ ) − κ2 m(p∗ , pˆ∗ ) + a(p, pˆ) − κ2 m(p, pˆ) ,

(5.20)

kde ∗

∗

∗

∗

Z

a(p , pˆ ) = ZΩ

m(p , pˆ ) = ZΩ a(p, pˆ) = ZΩ m(p, pˆ) =

(div(∇V)δij − ∂k Vj δki − δej ∂e Vi )∂j p∗ ∂i pˆ , p∗ pˆ div(∇V) , (5.21) (div(∇V)δij − ∂k Vj δki − δej ∂e Vi )∂j p∂i pˆ∗ , pˆ p∗ div(∇V) .

Ω

Touto cestou z´ıskáme u ´plnou derivaci u ´ˇcelové funkce Φ v˚ uˇci α, tj. jej´ı citlivost. V (5.20) derivace δΦ ve skuteˇcnosti znamená dτd Φ|τ =0 , kde Φ je implicitn´ı funkce hranice Γ(τ ) parametrizované vztahem (5.5) a závis´ı tedy na vektorovém poli V.

5.3 5.3.1

Implementace Adjungovan´ au ´ loha

V´ ypoˇcet citlivosti zmˇeny u ´ˇcelové funkce Φ jsme provedli ve SfePy. ˇ sen´ı p a zkuˇsebn´ı tlak q jsou komplexn´ı hodnoty. Hvˇezdiˇcka znamená komplexnˇe sdruReˇ ˇzenou hodnotu q = qr + i qi , q ∗ = qr − i qi , p = pr + i pi , p∗ = pr − i pi .

(5.22)

Pozdˇeji vyuˇzijeme tˇechto rovnost´ı p∗ q ∗ = (pq)∗ = ((pr qr − qi pi ) + i(pr qi + pi qr ))∗ , (∇p∗ , ∇q ∗ ) = (∇p, ∇q)∗ .

(5.23)

Reálné ˇcásti v adjungované u ´loze (5.10) bude vhodné pro implementaci zapsat takto Re {(∇(pr + i pi ), ∇(qr − i qi ))Ω } = (∇pr , ∇qr )Ω + (∇pi , ∇qi )Ω , Re κ2 (p, q)Ω = κ2 [(pr , qr )Ω + (pi , qi )Ω ] , Re (i hκp, qiΓin ∪Γout ) = −κ hpr , −qi iΓin ∪Γout + hpi , qr iΓin ∪Γout . 31

(5.24)


Implementace

Zu ´ˇcelové funkce ΦI také vybereme reálnou ˇcást Z Z ∗ 2Re p q =2 (pr qr + pi qi ) . ΩC

(5.25)

ΩC

Vyuˇzijeme (5.24), (5.25) a rozdˇel´ıme adjungovanou u ´lohu (5.10) na ˇreˇsen´ı zvláˇst’ reálné a imaginárn´ı ˇcásti. Zap´ıˇseme maticovˇe T (∇·, ∇·)Ω − κ2 (·, ·)Ω ( · , pr )Ωc −κ h·, ·iΓin ∪Γout pˆr qr qi , = − qr qi ( · , pi )Ωc κ h·, ·iΓin ∪Γout (∇·, ∇·)Ω − κ2 (·, ·)Ω pî (5.26) 2 coˇz plat´ı pro ∀(qi , qr ) ∈ Q .

5.3.2

Pˇ r´ıprava oblasti Ω

V MSC.Marc/Mentat jsme vygenerovali MKP s´ıt’ oblasti Ω se ˇctvercov´ ymi elementy

Obrázek 5.2: S´ıt’ MKP oblasti Ω; 12301 Obrázek 5.3: Oblast Ω s podoblastmi ΩC , uzl˚ u; 12000 element˚ u; 36903 stupˇ n˚ u volΩD , ΩA a ΩB nosti

V oblasti Ω budeme uvaˇzovat, jak jiˇz bylo zm´ınˇeno dˇr´ıve, podoblasti ΩC a ΩD . V ΩC zavedeme jeˇstˇe ΩA ⊂ ΩC a ΩB ⊂ ΩC obˇe o ˇs´ıˇrce jen dva elementy. ´ celové funkce ΦI a ΦII , viz (5.4), budeme vyhodnocovat pˇres ΩA a ΩB . Uˇ

5.3.3

Spline-box

Spline-box pomáhá vytváˇret geometrickou parametrizaci 2D i 3D oblast´ı. Pouˇzili jsme spline-box implementovan´ y v systému SfePy. Vyuˇzijeme ho k vytvoˇren´ı parametrizace tvaru hranice ΓD pˇri hledán´ı jej´ıho optimáln´ıho designu. Plat´ı ΩD ⊂ Ω, viz obr.(5.1), která je 32


Implementace

tvoˇrena 2D MKP s´ıt´ı, viz obr.(5.2). Jako zdroj k tomuto tématu jsem vyuˇzil [Rohan(2007)]. Uvaˇzujeme oblast ΩD , viz obr.(5.1), a parametrizaci jej´ı hranice ΓD , viz (5.5). Spline 2D reprezentaci pro nezmˇenˇenou Ω0D zap´ıˇseme x~0 (τ, s) =

4 X 4 X i

Ni (τ )Nj (s)~b0ij ,

(5.27)

j

kde N (τ ) je B-spline báze svisl´ ych splin˚ u a N (s) vodorovn´ ych. Poˇca´teˇcn´ı polohu ˇr´ıd´ıc´ıch 0 bod˚ u urˇcuje ~bij . Polohu uzl˚ u nezdeformované s´ıtˇe popisuje ~x. Vytvoˇrili jsem tedy nad Ω0D 2D Spline-box s 16 ˇr´ıd´ıc´ımi body, viz obr.(5.4)

Obrázek 5.4: Um´ıstˇen´ı ˇr´ıd´ıc´ıch bod˚ u po vytvoˇren´ı spline-boxu nad ΩD

Zmˇenu polohy ˇr´ıd´ıc´ıch bod˚ u pop´ıˇseme ~bij = ~b0 + d~ij αij . ij

(5.28)

Smˇer zmˇeny um´ıstˇen´ı (pohybu) ˇr´ıd´ıc´ıch bod˚ u urˇcuj´ı jednotkové vektory d~ij a velikost této zmˇeny parametry αij . Pro nenulové αij plat´ı ~x(τ, s) =

4 X 4 X i

Ni (τ )Nj (s)~bij .

(5.29)

j

Kdyˇz zderivujeme (5.29) podle αij , dostaneme ∂~x = Ni (τ )Nj (s)d~ij . ∂αij

(5.30)

Z´ıskáme tedy v ΩD vektorové pole V(τ, s)ij od parametru αij V(τ, s)ij = Ni (τ )Nj (s)d~ij . 33

(5.31)


Implementace

V naˇsem pˇr´ıpadˇe budeme mˇenit tvar jen hranice ΩD . P˚ ujde tedy o 1D u ´lohu a budeme uvˇzovat jen 4 parametry αi , i = 1, 2, 3, 4. Ztotoˇzn´ıme parametry αi s τ . Z toho vypl´ yvá, ˇze u ´plná derivace u ´ˇcelové funkce Φ je vlastnˇe derivac´ı ve smˇeru vektorového pole V. δαi Φ = δτ Φ = δΦ ◦ V i .

(5.32)

S ohledem na obr.(5.4) budeme uvaˇzovat posun ˇr´ıd´ıc´ıch bod˚ u ve vodorovném smˇeru. Paˇ ıd´ıc´ı body rametr αi proporcionálnˇe ˇrekne, o kolik se posune i-tá ˇrada ˇr´ıd´ıc´ıch bod˚ u. R´ oznaˇc´ım Bi,j . Kdyˇz bod B4,4 posuneme o α4 ve smˇeru d~44 = (−1, 0), tak B3,4 se posune o 2/3α4 , B2,4 o 1/3α4 a B1,4 o 0, viz obr.(5.5).

Obrázek 5.5: Posun ˇrady ˇr´ıd´ıc´ıch bod˚ u pro α4

5.3.4

Porovn´ an´ı koneˇ cn´ ych diferenc´ı a citlivostn´ı anal´ yzy

Z´ıskán´ı citlivosti u ´ˇcelové funkce Φ na parametrech α znamená z´ıskán´ı totáln´ı derivace δτ Φ. To jde nejen pomoc´ı citlivostn´ı anal´ yzy (CA), ale také numericky pomoc´ı koneˇcn´ ych diferenc´ı (KD). Druhá moˇznost v´ ypoˇctu je vˇsak v optimalizaci v´ ypoˇcetnˇe nároˇcnˇejˇs´ı a ménˇe pˇresná. Porovnáme KD s CA a ovˇeˇr´ıme si t´ım správnost v´ ypoˇctu citlivosti u ´ˇcelové funkce pomoc´ı CA. Pˇri pouˇzit´ı koneˇcn´ ych diferenc´ı budeme uvaˇzovat velmi malou zmˇenu dτ na hranici ΓD . Posuneme ˇr´ıd´ıc´ı body, viz obr.(5.4), na ΓD o dτ nejdˇr´ıve ve smˇeru d~+ = (1, 0) a poté d~− (−1, 0). Ω+ D = ΩD + dτ {V(x)}x∈ΩD , Ω− D = ΩD − dτ {V(x)}x∈ΩD .

(5.33)

Postupnˇe vyhodnot´ıme u ´ˇcelovou funkci (5.4)1 v nov´ ych oblastech Ω+ a Ω− . Z´ıskáme Φ+ a − Φ . Provedeme centráln´ı pomˇernou diferenci, jej´ıˇz v´ ysledkem bude hledaná citlivost δτ Φ =

Φ+ − Φ− . 2dτ

(5.34)

Pˇresnost závis´ı samoˇzˇrejmˇe na volbˇe velikosti dτ . Porovnali jsme v´ ysledky cilivostn´ı anal´ yzy a koneˇcn´ ych diferenc´ı pro nˇekolik frekvenc´ı ω na oblasti obr.(5.3). V´ ysledné citlivosti a jejich rozd´ıly jsou v následuj´ıc´ı tabulce.

34


Optimalizaˇcn´ı výpoˇcty

dτ

ω[rad s− 1]

Citlivostn´ı anal´ yza

Koneˇ cn´ e diference

Rozd´ıl(CA-KD)

1e-6

200 400 600 800 1000

-1813.56452583+0i 16562.4516291-8.73e-11i 20675.4802695-5.82e-11i 19822.959634-5.82e-11i 44579.6204201-1.75e-11i

-1813.55618406+0i 16561.4799553+0i 20675.2863232+0i 19822.8649024+0i 44578.836579+0i

-0.00834177+0i 0.97167380-8.73e-11i 0.19394625-5.82e-11i 0.09473160-5.82e-11i 0.78384415-1.75e-11i

5.3.5

Optimalizace

V´ ypoˇcet provedeme v softwaru SfePy a Matlab. SfePy bude vypoˇc´ıtávat citlivost u ´ˇcelové funkce na zmˇenu parametr˚ u α. Citlivost se pˇredá Matlabu, kde se vyuˇzije knihovn´ı optimalizaˇcn´ı funkce fmincon, viz [MathWorks(2013)], která pouˇzije citlivost k optimalizaci. Matlab bude vracet nové hodnoty optimalizovan´ ych parametr˚ u α do SfePy, kde zaˇcne dalˇs´ı iterace.

5.4 5.4.1

Optimalizaˇ cn´ı v´ ypoˇ cty ´ Uvodem

V´ ypoˇcty budeme provádˇet na oblasti Ω (obr.(5.1)), rozdˇelené na podoblasti (obr.(5.3)). Na hranici Γ = Γin ∪ Γout ∪ Γ0 pouˇzijeme okrajové podm´ınky (4.7). Zavedeme spline-box nad ΩD , viz obr.(5.4). Optimáln´ı tvar hranice ΓD budeme hledat pro minimum u ´ˇcelov´ ych funkc´ı ΦI a ΦII , viz (5.4). Plat´ı Z |p|2 = kpk2ΩA ∪ΩB [P a2 m2 ] , ΦI (p) = ΩA ∪ΩB ! (5.35) kpkΩA [dB] , ΦII (p) = TAB (p) = 20 log kpkΩB kde ΦII se rovná funkci transmisn´ıch ztrát TAB , která ˇriká, o kolik decibel˚ u [dB] se utlum´ı ´ celovou funkci ΦI (TAB > 0) nebo zes´ıl´ı (TAB < 0) zvuková vlna jdouc´ı z ΩA do ΩB . Uˇ budeme poˇc´ıtat pˇres ΩA ∪ΩB a ΦII dle (5.35)2 . V´ ypoˇcet provedeme pro tˇri u ´hlové frekvence −1 ω = 200, 600, 1000 [rads ].

35



Parametry oblasti Ω navol´ıme následovnˇe c = 342 [m/s] , ω κ = [rad/m] , c uA = 50 [−] , l = 15 [m] , v = 10 [m] , a = 5 [m] , b = 2 [m] .

5.4.2

(5.36)

Vlastn´ı v´ ypoˇ cty

Nejprve vypoˇcteme ˇreˇsen´ı |p| v neoptimalizované oblasti Ω a z´ıskáme optimalizac´ı neoN vlivnˇené hodnoty u ´ˇcelov´ ych funkc´ı ΦN I a ΦII . Pak provedene optimalizace pro ΦI a ΦII . O Z optimalizované oblasti Ω z´ıskáme optimalizované hodnoty u ´ˇcelov´ ych funkc´ı ΦO I a ΦII . Zobraz´ıme závislost velikosti Φ na iteraˇcn´ım kroku. Toto provedeme pro vˇsechny tˇri ω.

Obrázek 5.6: Zobrazeno |p| v neoptimalizované oblasti Ω; ω = 200 [rads−1 ]

Obrázek 5.7: Zobrazeno |p| v optimalizo- Obrázek 5.8: Zobrazeno |p| v optimalizované oblasti Ω pro ΦI ; ω = 200 [rads−1 ] vané oblasti Ω pro ΦII ; ω = 200 [rads−1 ]

36



Obrázek 5.9: Graf závislosti hodnoty ΦI na Obrázek 5.10: Graf závislosti hodnoty ΦII iteraˇcn´ım kroku optimalizace na iteraˇcn´ım kroku optimalizace



37



Obrázek 5.14: Graf závislosti hodnoty ΦI Obrázek 5.15: Graf závislosti hodnoty ΦII na iteraˇcn´ım kroku optimalizace na iteraˇcn´ım kroku optimalizace



38



Obrázek 5.19: Graf závislosti hodnoty ΦI Obrázek 5.20: Graf závislosti hodnoty ΦII na iteraˇcn´ım kroku optimalizace na iteraˇcn´ım kroku optimalizace

N´ıˇze vid´ıme tabulku, která ukazuje hodnoty u ´ˇcelov´ ych funkc´ı (5.35) pˇred (ΦN ) a po (Φ ) optimalizaci a jejich procentuáln´ı srovnán´ı, které ˇriká, o kolik procent se zvˇetˇsila ˇci zmenˇsila hodnota Φ. O

2 2 2 2 O ω [rad s− 1] ΦN Pa m Φ Pa m I I

±%

ΦN II [dB]

ΦO II [dB]

±%

200 600 1000

-11.92 -6.49 -10.55

0.921 9.833 3.8986

-11.7089 0.9845 -4.8617

-1371.32 -89.99 -224.70

5.1383e+4 5.1869e+4 5.2944e+4

4.5426e+4 4.85002e+4 4.7361e+4

Podle tabulky vid´ıme, ˇze optimalizace, podle oˇcekáván´ı, sn´ıˇzila hodnoty obou u ´ˇcelov´ ych funkc´ı pro vˇsechny u ´hlové frekvence ω. Postup minimalizován´ı ΦI,II je vidˇet z graf˚ u, viz obrázky (5.9), (5.10), (5.14), (5.15), (5.19) a (5.20). Na obrázc´ıch (5.7), (5.8), (5.12), (5.13), (5.17) a (5.18) je viditelná zmˇena geometrie designové hranice ΓD , coˇz zapˇr´ıˇcinilo zˇrejmou zmˇenu barevné mapy ˇreˇsen´ı. Nˇekteré v´ ysledné optimalizované hodnoty ΦII jsou záporné, protoˇze optimalizac´ı se dosáhlo toho, ˇze akustick´ y tlak v ΩA je menˇs´ı neˇz v ΩB neboli ˇ kpkΩA < kpkΩB . C´ımˇz nab´ yvá logaritmus (5.35)2 záporné hodnoty.

39

6 Závˇer V druhé kapitole jsme se seznámili s odvozen´ım rovnic potˇrebn´ ych k z´ıskán´ı vlnové rovnice pro tekutiny: Eulerovo pohybovou rovnic´ı, rovnic´ı kontinuity, stavovou rovnic´ı. Z nich jsme odvodili vlnovou rovnici a následnˇe i Helmholtzovu rovnici pro stojaté vlnˇen´ı. V tˇret´ı kapitole jsme analyticky vyˇreˇsili vlnovou a Helmholtzovu rovnici v 1D pomoc´ı Fourierovu metody a zobrazili nˇekterá ˇreˇsen´ı v Matlabu. Ve ˇctvrté kapitole jsme ˇreˇsili slabé formulace Helmholtzovy rovnice na oblasti Ω. Z u ´vah o ˇs´ıˇren´ı akustick´ ych vln jsme si odvodili okrajové podm´ınky na hranici ∂Ω, pomoc´ı MKP provedli ˇreˇsen´ı v systému SfePy. Na dané oblasti jsme vyzkouˇseli r˚ uzné frekvence ω inciˇ sen´ı jsme si zobrazili v programu ParaView. denˇcn´ı vlny. Reˇ V páté kapitole jsme naformulovali u ´lohu tvarové optimalizace pro akustické pole, zparametrizovali design pomoc´ı spline-boxu. Zavedli pole designov´ ych rychlost´ı. Optimalizaˇcn´ı u ´loze jsme pˇriˇradili Lagrangeovu funkci a pˇreˇsli k ˇreˇsen´ı adjungované u ´lohy. Porovnali jsme pro kontrolu v´ ysledky citlivostn´ı anal´ yzy a koneˇcn´ ych diferenc´ı. Oba pˇr´ıstupy jsme implementovali ve SfePy. Ujistili jsme se, ˇze v´ ysledky citlivostn´ı anal´ yzy a koneˇcn´ ych diferenc´ı si celkem odpov´ıdaj´ı. V závˇeru kapitoly jsme provedli nˇekolik optimalizaˇcn´ıch v´ ypoˇct˚ u ve 2D oblasti Ω v Matlabu a SfePy pro nˇekolik frekvenc´ı ω. Vyuˇzili jsme dvˇe u ´ˇcelové funkce ΦI,II . Naˇsli jsme nˇekolik optimáln´ıch tvar˚ u designové hranice ΩD . V´ ysledky jsme zobrazili programem ParaView. Vykreslili jsme závislosti hodnot ΦI,II na iteraci optimalizace. Z graf˚ u je vidˇet, ˇze hodnoty správnˇe klesly u vˇsech v´ ypoˇct˚ u. Vypoˇcetli jsme, o kolik procent klesly optimalizac´ı hodnoty ΦI,II . Zprovoznili jsme tedy u ´spˇeˇsnˇe tvarovou optimalizaci akustického pole. Práce nemá zat´ım urˇcitou konkrétn´ı aplikaci. Jde v n´ı sp´ıˇse jen o zprovoznˇen´ı a vyzkouˇsen´ı samotné optimalizace. Má urˇcitˇe mnoho prostoru pro vylepˇsen´ı a rozˇs´ıˇren´ı. Tˇreba pˇreveden´ı na v´ıcerozmˇernou u ´lohu se sloˇzitˇejˇs´ı geometri´ı, kde se objev´ı v´ıce optimalizaˇcn´ıch parametr˚ u. Poté by se naˇsla urˇcitˇe nˇejaká potenciáln´ı aplikace pro tuto práci. Jednou z moˇznost´ı je optimalizace tvaru skˇr´ınˇe nˇejakého stroje. Ve skˇr´ıni bude zdroj vibrac´ı (nˇejaká nevyváˇzená rotuj´ıc´ı hmotnost), které se budou konstrukc´ı pˇrenáˇset na skˇr´ıˇ n, t´ım by vznikala emise hluku, kterou bychom chtˇeli regulovat. V tomto smˇeru by mohla pokraˇcovat má dalˇs´ı práce.

40

Literatura [B. Engquist(1977)] B. ENGQUIST, A. M. Absorbing boundary conditions fot the numerical simulation of waves. 1977. Dostupné z: http://www.math.mcgill.ca/~gantumur/ docs/down/Engquist77.pdf. [D. Cioranescu(2008)] D. CIORANESCU, G. G. A. D. The periodic unfolding method in homogenization. 2008. ´ [Drábek(1994)] DRáBEK, P. Uvod do funkcionáln´ı analýzy. Plzeˇ n : Západoˇceská univerzita v Pzni, 1994. ISBN 80-7082-124-8. [E. Bängtsson(2002)] E. BäNGTSSON, M. B. D. N. Shape optimization of an acoustic horn. 2002. Dostupné z: http://www.sciencedirect.com/science/article/pii/ S0045782502006564. [E. Rohan(2006)] E. ROHAN, B. M. Homogenization and shape sensitivity of microstructures for design of piezoelectric bio-materials. 2006. [E. Rohan(2012)] E. ROHAN, V. L. Sensitivity analysis for optimal design of perforated plates in vibro-acoustics: homogenization approach. 2012. [E. Rohan(2010a)] E. ROHAN, V. L. Homogenization of the acoustic transmission through perforated layer. 2010a. [E. Rohan(2010b)] E. ROHAN, V. L. Sensitivity analysis for acoustic waves propagating through homogenized thin perforated layer. 2010b. [E. Rohan(2011a)] E. ROHAN, V. L. Homogenized perforated interface in acoustic wave propagation – modeling and optimization. 2011a. [E. Rohan(2011b)] E. ROHAN, V. L. Homogenization of the vibro–acoustic transmission on perforated Reissner-Mindlin plate. 2011b. [J. Haslinger(1996)] J. HASLINGER, P. N. Finite Element Approximation for Optimal Shape, Material and Topology Design. 1996.

41

LITERATURA

LITERATURA

ˇ ˇ [Skvor(2001)] SKVOR, Z. AKUSTIKA A ELEKTRO-AKUSTIKA. Praha : Academia, 2001. ISBN 80-200-0461-0. [Linhart(2009)] LINHART, J. Mechanika tekutin I. Plzeˇ n : Západoˇceská univerzita v Pzni, 2009. ISBN 978-80-7043-766-7. [MathWorks(2013)] MATHWORKS. Constrained optimization - fmincon, 2013. [M´ıka(1983)] M´ıKA, S. Parciáln´ı diferenciáln´ı rovnice I. Praha : SNTL, 1983. [M´ıka(2007)] M´ıKA, S. Numerické metody ˇreˇsen´ı okrajových u ´loh pro ODR. Plzeˇ n : Západoˇceská univerzita v Pzni, 2007. [Rohan(2007)] ROHAN, E. SPBOX, 2007. [Rohan(2012)] ROHAN, E. Citlivostn´ı analýza pro optimalizaci v mechanice kontinua, 2012.

42

Západočeská univerzita v Plzni Fakulta aplikovaných věd Katedra mechaniky. Optimalizace akustického pole

Recommend Documents