MHP akustická analýza řešena s využitím programů ANSYS a Matlab BEM acoustics analysis solved using programs ANSYS and Matlab Jan SZWEDA 1)
Abstrakt Tento článek je zaměřen na popis jednoho ze způsobů numerického řešení šíření zvuku, které je vyvoláváno vibracemi povrchu strojních celků. Základem pro numerické akustické analýzy je metoda hraničních prvků. Smyslem práce bylo zpřístupnit uživateli programu ANSYS rychlé a nenáročné řešení akustické analýzy pro určování vyzářené akustické energie od mechanických vibrací tělesa. Cílem práce je posouzení míry nepřesnosti výsledků získaných užitím navrženého postupu, který je založen na jednoduchém typu hraničních prvků. Pro naplnění cíle je provedeno řešení akustické úlohy dvěma způsoby, každým pro jiný typy hraničních elementů, kdy jedno z řešení je provedeno programem Sysnoise v5.6. Průměrná relativní chyba v hodnotách hladiny akustického výkonu se pohybuje na 5% hranici, ojedinělá maximální hodnota této chyby činí 16 %. Klíčová slova: vibrace, MKP; Helmholcova rovnice, MHP, akustická analýza; ANSYS, Sysnoise Summary This paper deals with description of the way to numerical solution of sound waves propagation caused by vibrating surfaces of the bodies. The numerical acoustic analyses are based on the boundary element method. The purpose of this work is to provide an efficient method for solution of emitted acoustic energy, which could be used by programme ANSYS. The aim of this work is accuracy rate evaluation of results obtained by suggested technique that approximated boundary by constant boundary elements. Two ways of acoustic analysis solution, each with another type of boundary elements, are made to meet the defined goal. One solution is done by programme Sysnoise 5.6. The average relative error of acoustic energy level is approximately 5 %; single maximum value of this error reaches 16 %. Key words: vibration, FEM; Helmholc equation, BEM, acoustic analysis; ANSYS, Sysnoise
1)
VŠB-Technická univerzita Ostrava, Fakulta strojní, katedra mechaniky, 17. listopadu 15, 708 33 Ostrava-Poruba, Email:
[email protected]
13. ANSYS Users’ Meeting, 21. – 23. září 2005 Přerov -1-
1.
ÚVOD
Jedním ze zákonitých projevů strojů a strojních dílů jsou jejich mechanické vibrace. Důsledkem vibrací je jednak nepříznivé namáhání strojních částí na únavu, ale také vznik a šíření zvuku, který v případě jeho nepříjemné intenzity nazýváme hlukem. Tento článek je zaměřen na využití numerických metod k řešení šíření zvuku vyzařovaného vibracemi povrchu strojních celků. Základem pro numerické řešení šíření zvuku je metoda hraničních prvků (zkr. MHP). Ve článku jsou porovnávány výsledky získané užitím komerčního programu LMS Sysnoise a výsledky získané autorským užitím MHP, které bylo zpracováno v prostředí programu Matlab. Pro oba přístupy byly zdrojem zvuku rychlosti kmitání povrchu železničního kola, které bylo v radiálním směru buzeno sílou harmonického průběhu. Rychlosti vibrací vyzařovací plochy byly zjištěny pro různé frekvence budicí síly, k řešení byl použit program ANSYS. Cílem práce je zjištění míry nepřesnosti výsledků získaných užitím jednoduchého typu hraničních prvků a posouzení prezentovaného zjednodušeného modelování vyzařování hluku pro vzájemné porovnávání železničních kol z hlediska jejich akustických vlastností. Rozhodujícím kritériem jsou grafy a hodnoty hladin vyzařované akustické energie a velikosti hladin akustického tlaku v referenční rovině. 2.
AKUSTICKÉ ÚLOHY – MHP Je uvažováno akustické vlnové pole, ve kterém jsou akustické vlny vztaženy k lokálnímu pohybu tekutiny a nezávisí na jejím celkovém pohybu. Pro proměnné veličiny vlnového pole, tj. tlak p, rychlost v a hustota ρ, se uvažuje jejich fluktuace kolem stření hodnoty p0, v0 , ρ0. Pro nevizkózní tekutinu má pohybová rovnice, linearizovaná vzhledem k fluktuacím p, v, ρ, tvar tzv. Eulerovy pohybové rovnice ρ0 ·∂ v/∂ t = –∇p, (1) kde p a v značí fluktuace tlaku, resp. rychlosti (tzv. akustický tlak a rychlost), t značí čas a ρ0 značí střední hodnotu hustoty tekutiny. Pro určení konstitutivního vztahu lze šíření akustických vln považovat za adaibatický děj s konstantní entropií. Zanedbáním vyšších členů Taylorova rozvoje hustoty nabude konstitutivní vztah izotropní tekutiny tvar p = –ρ0 ·c2·∇v, (2) kde p a v značí fluktuace tlaku, resp. rychlosti, ρ0 střední hodnotu hustoty a c značí rychlost šíření akustických vln. Řešení problému šíření akustických vln může být provedeno zavedením funkce tzv. rychlostního potenciálu v = –∇Φ. Její užitím při úpravě vztahů (1) a (2) se úloha šíření akustických vln vyjádři v časové oblasti ve tvaru 1 ∂ 2Φ ∇ Φ− 2 2 =0. (3) c ∂t Za předpokladu, že rychlostní potenciál lze vyjádřit harmonickou funkcí Φ(r, t) = Φ(r)·e iω t, je možno rovnici (3) transformovat do frekvenční oblasti na tvar (∇2 + k2)·Φ(r) = 0, (4) kde k = ω /c značí vlnové číslo. V tomto vyjádření je Φ pouze funkcí polohy, tj. Φ(r), a akustický tlak lze určit podle vztahu p = –iω ρ0 Φ . (5) 2
Hledané řešení rovnice (4) musí splňovat okrajové podmínky: a) speciální podmínku v nekonečnu, tzv. Sommerfeldovu radiační podmínku ∂Φ
∂n
= ikΦ,
13. ANSYS Users’ Meeting, 21. – 23. září 2005 Přerov -2-
b) na povrchu tělesa může nastat některá z následujících podmínek: – Φ = ΦB: podmínka Dirichletova typu, – ∂Φ ∂n = q = qB: podmínka Neumannova typu, – lineární kombinace podmínek předchozích dvou typů. Užitím postupu vážených residuí lze řešení rovnice (4) převést do tvaru integrální rovnice hraničních elementů v podobě c j Φ j + ∫ Φ q*dΓ = ∫ q Φ *dΓ , Γ
(6)
Γ
kde j značí bod hranice Γ a Φ*, resp. q* je fundamentální řešení rovnice (4), které má pro 3D úlohu v izotropním médiu tvar 1⎞ e − ikr ⎛ e − ikr , resp. q * = (7) Φ* = ⎜ ik + ⎟ cos β n . 4π r ⎝ r⎠ 4π r Ve vztazích (7) a (8) značí r polohový vektor elementu dΓ vzhledem k bodu j a βn značí úhel mezi vektory r a n, což je vektor normály elementu dΓ. Jelikož fundamentální řešení je z oboru komplexních čísel jsou i ostatní veličiny rovnice (7) obecně komplexní. Lze odvodit, že imaginární část konstanty cj je rovna nule, kdežto hodnota reálné části této konstanty závisí na řádu hraničních prvků a na geometrií hranice v bodu j. Pro hraniční prvky 0. řádu, tzv. konstantní elementy, je Re(cj) = 1/2, zatímco pro hraniční prvky 1. řádu, tzv. lineární, je hodnota Re(cj) dána vztahem cj = θ 2 π , (8) kde θ značí úhel, který v bodě j svírají sousední elementy hranice. Diskretizací hranice Γ na hraniční elementy přejde integrální rovnice (6) na tvar sumační, který může být vyjádřen maticovým zápisem [H]·{Φ} = [G]·{q}. (9) Uplatněním okrajových podmínek na hranici tělesa a úpravou rovnice (9) obdržíme rovnici pro řešení neznámých veličin na hranici. Pro úlohu šíření vln vybuzených harmonickými vibracemi povrchu tělesa nabude rovnice (9) tvaru [H]·{Φ} = [G]·{q}B = {b}, (10) kde {q}B jsou známé hodnoty normálových rychlostí na vyzařovací hranici a {Φ} jsou hledané hodnoty rychlostního potenciálu na vyzařovací hranici, ze kterých lze akustický tlak určit zpětně podle (5). Tímto je získáno úplné řešení akustické úlohy na vyzařovací hranici a toto řešení může být dále použito pro zjištění hodnot akustických parametrů v požadovaných místech akustického prostoru, podrobně viz [2]. Hodnoty okrajových podmínek pro sestavení rovnice (9), tzn. velikosti normálových rychlostí vyzařovací hranice lze najít řešením odezvy na harmonické buzení vyzařovacího tělesa. Pomocí MKP lze pro buzení {f(t)} = {fa}·e iω t najít odezvu ve tvaru {u(t)} = {ua}·e iω t řešením rovnice − ω 2 ⋅ [M ]⋅ {u a }+ iω ⋅ [B]⋅ {u a }+ [K ]⋅ {u a } = {f a } , (11) kde fa jsou amplitudy budicích silových účinků kruhové frekvence ω ; M, B a K jsou matice setrvačných, tlumicích a tuhostních vlastností tělesa a ua jsou hledané amplitudy kmitání.
3.
UŽITÉ PROGRAMOVÉ ZPRACOVANÍ MHP AKUSTICKÉ ANALÝZY Stejně jako MKP je i MHP metodou diskretizační, tzn. že numerické řešení je prováděno na síti prvků, kterými je popsána předmětná oblast kontinua, a to způsobem odpovídajícím požadavkům vlastní metody. V případě MHP akustické analýzy je potřeba diskretizovat hranici vyšetřovaného akustického prostoru. U volného otevřeného akustického prostoru, ve 13. ANSYS Users’ Meeting, 21. – 23. září 2005 Přerov -3-
kterém se šíří vlny vznikající kmitáním povrchu pružného tělesa, je hranice totožná s vyzařovacím povrchem tělesa. Pro tento typ úloh jsou okrajovými podmínkami normálové složky rychlosti kmitání vyzařovacího povrchu a řešením akustické analýzy je získáno rozložení akustického tlaku na vyzařovacím povrchu. Dalšími výpočty s akustickým řešením na hranici lze mezi jinými zjistit množství vyzářené akustické energie, a to ve složkách její aktivní a reaktivní části. Pomocí akustického řešení na hranici je možné zjistit i hodnoty akustických veličin v definovaných bodech akustického prostoru, tzn. že pro tuto situaci se definuje rovněž síť bodů, resp. prvků v akustickém prostoru. V těchto bodech se nepředepisují žádné okrajové podmínky. Řešení akustické analýzy představuje po teoretické stránce sestavení a řešení rovnice (11) pro zjištění normálových rychlostí kmitání vyzařovacího povrchu a rovnice (10) pro zjištění akustických tlaků na vyzařovacím povrchu. Po praktické stránce je obvyklé, že každá z rovnic je řešena v jiném programu, přičemž geometrie hranice a okrajové podmínky jsou do programu pro akustickou analýzu exportovány z předchozí analýzy mechanického kmitání tělesa. Z komerčně dostupných programů bylo pro tento účel použito dvojice programu ANSYS 5.7 a Sysnoise 5.6. Tímto postupem bylo získáno referenční řešení problému, které bude dále použito pro posouzení přesnosti výsledků získaných autorským zpracováním akustické analýzy pomocí programu Matlab. MKP model tělesa
MKP model tělesa vč. definice vazby na akustickou analýzu
MKP model tělesa + info pro export vn
export geometrie akustické úlohy
Geometrie akustické úlohy pro MHP analýzu
Řešení úlohy harmonického kmitání tělesa export okrajových podmínek
MHP analýza – spustit Sestavení MHP modelu akustického problému Zpracování okrajových podmínek pro MHP model Řešení MHP úlohy a export výsledků
import výsledků akustické analýzy
Zpracování výsledků MHP výsledků a jejich vizualizace
Obr.1 Blokové schéma pro užití ANSYS k řešení akustické analýzy pomocí MHP
Způsob řešení akustické analýzy pomocí komerčního programu Sysnoise sice představuje relativně jednoduchou a zaručenou variantu řešení, s výhodami i nevýhodami plynoucími z užití komerčního programu. Jako nevýhodu lze uvést již kombinaci dvou programů různých výrobců, která má za následek téměř znemožnění automatizaci výpočtů pro různé tvarové varianty tělesa, nebo přímo optimalizaci tvaru tělesa s ohledem na množství vyzářené akus13. ANSYS Users’ Meeting, 21. – 23. září 2005 Přerov -4-
tické energie. Z tohoto důvodu bylo navrženo a odzkoušeno alternativní řešení MHP akustické analýzy, ve kterém je převážná část operací prováděna v programu ANSYS a pouze pro sestavení a samotné řešení rovnice akustické analýzy (10) je použito programu Matlab. Vzhledem k vlastnostem obou programů je možno v krajním případě využít této jejich kombinace i k řešení zmíněné optimalizace. Blokové schéma alternativního řešení je znázorněno na obr. 1. Na tomto diagramu jsou jednotlivé kroky uspořádány chronologicky ve svislém směru, zatím co vodorovné rozmístění ve 3 krocích znázorňuje postupné užití proramů ANSYS–Matlab–ANSYS. Mimo to jsou tučným orámováním zvýrazněny úkony, které musí provést uživatel a kurzívou jsou zvýrazněny procesy výměny dat, které jsou prováděny pomocí speciálních souborů.
4.
UKÁZKA APLIKACE MHP NA ŘEŠENÍ VYZAŘOVÁNÍ HLUKU
Aplikace popsaného postupu řešení úlohy šíření akustických vln bude předvedena na problému vyzařování hluku železničním kolem. Uvažovaný akustický problém byl řešen pro případ vibrací kola od harmonického buzení radiální silou. Úloha byla řešena pro vybrané frekvence buzení z úzkého pásma v blízkosti každé vlastní frekvence kola v rozsahu 200 ÷ 6500 Hz. Uvedenou volbou budicích frekvencí je sledován cíl, kterým je stanovení extrémních hodnot spektra akustické odezvy kola vyvolané radiálním buzením. Geometrie úlohy a okrajové podmínky jsou patrny z obr. 2. Pomocí MKP bylo v programu ANSYS provedeno řešení mechanických vibrací, jehož výsledky v podobě rychlosti kmitání povrchu byly použity jako okrajová podmínka v úloze šíření akustických vln. Tato úloha byla řešena dvěma nástroji: prvním byl program LMS Sysnoise, který diskretizuje hranici lineárními hraničními prvky a druhým nástrojem byl autorský program pro MHP zpracovaný v prostředí Matlab, kde byly použity konstantní hraniční prvky.
F(t) = Fa·eiω t
Obr.2 Model vyzařovací struktury kola
Obr.3 Rozložení Lp na povrchu kola– 960 Hz
Výsledky řešení akustické úlohy oběma nástroji jsou prezentovány na obr. 3 až 5. Pro názornost je na obr. 3 zobrazeno rozložení hladiny akustického tlaku na vyzařovacím povrchu kola. Z hlediska posouzení akustických vlastností kola je důležitějším výsledkem graf na obr. 4. Graf na obr. 5 dále kvantifikuje velikost odchylky ve vypočteném vyzářeném výkonu pro řešení pomocí lineárních a konstantních hraničních prvků.
13. ANSYS Users’ Meeting, 21. – 23. září 2005 Přerov -5-
ACTIVE POWER LEVEL Lw (dB Lin)
F R E Q U E N C Y (H z)
Obr.4 Spektrum hladin vyzářeného akustického výkonu – řešeno pro diskrétní frekvence
Obr.5 Odchylky ve výsledcích Lw mezi lineárními a konstantními hraničními prvky
5.
ZÁVĚR
Cílem práce je posouzení nepřesnosti výsledků akustické analýzy získaných konstantními hraničními prvky. Porovnáním spekter a hodnot hladiny akustického výkonu bylo zjištěno, že průměrná relativní chyba se pohybuje na 5% hranici. Ojedinělé maximum této chyby činí 16 %. Z ohledem na tyto výsledy lze konstatovat, že užití konstantních elementů MHP pro řešení akustických problémů dává výsledky dostatečně přesné pro účely inženýrské praxe.
PODĚKOVÁNÍ Práce vznikla na námět a za podpory společnosti Bonatrans, a.s., jak rovněž za podpory MŠMT ČR prostřednictvím projektu CEZ:J17/98:2724019.
LITERATURA [1] Ciskowski R. D., Brebbia C. A.: Boundary Element Methods in Acoustics. Computational Mechanics Publications, Southampton, 1991. ISBN 1-85312-104-5. [2] Mišun V.: Vibrace a hluk. Fakulta strojní VUT v Brně, Brno 1998. ISBN 80-214-1262-3. [3] Biskup M.: Numerické modelování šíření zvuku a inženýrské aplikace. VŠB-TU Ostrava, Ostrava, 2003. ISBN 80-248-0332-1.
13. ANSYS Users’ Meeting, 21. – 23. září 2005 Přerov -6-