RADIO-AEROSZOLOK LÉGÚTI LOKÁLIS KIÜLEPEDÉSÉNEK VIZSGÁLATA NUMERIKUS ÁRAMLÁSTANI MÓDSZEREKKEL

RADIO-AEROSZOLOK LÉGÚTI LOKÁLIS KIÜLEPEDÉSÉNEK VIZSGÁLATA NUMERIKUS ÁRAMLÁSTANI MÓDSZEREKKEL

Farkas Árpád

EÖTVÖS LORÁND TUDOMÁNYEGYETEM TERMÉSZETTUDOMÁNYI KAR FIZIKA DOKTORI ISKOLA STATISZTIKUS FIZIKA, BIOLÓGIAI FIZIKA ÉS KVANTUMRENDSZEREK FIZIKÁJA PROGRAM

DOKTORI ISKOLA VEZETŐJE: Dr. Horváth Zalán PROGRAMVEZETŐ: Dr. Vicsek Tamás

TÉMAVEZETŐ Dr. Balásházy Imre Tudományos Főmunkatárs, MTA KFKI AEKI

A DOLGOZAT KÉSZÍTÉSÉNEK HELYE Magyar Tudományos Akadémia KFKI Atomenergia Kutatóintézet

Tartalomjegyzék 1. BEVEZETÉS, CÉLKITŰZÉS......................................................................................4 I. RÉSZ: A NUMERIKUS ÁRAMLÁSTANI ÉS RÉSZECSKE KIÜLEPEDÉSI MODELL FELÉPÍTÉSE ÉS VALIDÁLÁSA...............................................................11 2. A NUMERIKUS MODELLEK SZÜKSÉGESSÉGE...............................................11 3. LÉGÚTI GEOMETRIÁK MODELLEZÉSE...........................................................13 3.1. A légzőrendszer szerepe és felépítése......................................................................13 3.2. A tüdőfa matematikai modelljei...............................................................................14 3.3. Digitalizált háromdimenziós centrális légúti-geometria modellek..........................15 3.4. A felsőlégutak modellezése.....................................................................................23 4. A NUMERIKUS RÁCS……………………………………………………………...25 5. A LEVEGŐ ÁRAMLÁSI TERÉNEK NUMERIKUS ÁRAMLÁSTANI MODELLJE.....................................................................................................................29 5.1. Megmaradási törvények...........................................................................................29 5.1.1. Tömegmegmaradás (folytonosság) egyenlete.................................................29 5.1.2. Impulzusmegmaradás egyenlete (mozgásegyenlet)........................................29 5.2. Az áramlási tér meghatározása................................................................................30 5.2.1. A transzportegyenlet diszkretizálása...............................................................32 5.2.2. A konvekciós tag kiszámítása..........................................................................33 5.2.3. A diffúziós tag kiszámítása..............................................................................34 5.2.4. Peremfeltételek................................................................................................35 5.2.5. Lineáris algebrai egyenletek és megoldásuk..................................................36 5.2.6. Idődiszkretizálás instacioner áramlás esetén.................................................38 5.2.7. A turbulencia kérdése.....................................................................................39 6. RÉSZECSKE TRAJEKTÓRIA ÉS KIÜLEPEDÉSI MODELL…………………44 6.1. Kiülepedési mechanizmusok és matematikai leírásmódjuk....................................44 6.2. A mozgásegyenlet numerikus megoldása................................................................51 6.3. Részecskesorsolás és peremfeltétel.........................................................................52 7. A NUMERIKUS MODELL VALIDÁLÁSA………………………………………53 7.1. A Navier-Stokes egyenlet megoldásának validálása...............................................53 7.2. A trajektória modell validálása................................................................................57 II. RÉSZ: A NUMERIKUS MODELL ALKALMAZÁSA.........................................63 8. LÉGÚTI LEVEGŐÁRAMOK SZÁMÍTÁSA..........................................................63 9. RADIONUKLIDOK LÉGÚTI KIÜLEPEDÉSÉNEK LEÍRÁSA..........................70 9.1. Extrathoratikus aeroszol kiülepedés........................................................................70 9.2. A részecskeméret hatása a bronchiális kiülepedésre...............................................72 9.3. A térfogatáram hatása a bronchiális kiülepedésre...................................................73

2

9.4. A légúti morfológia hatása a bronchiális kiülepedésre............................................75 9.5. Részecskekiülepedés belégzéskor és kilégzéskor....................................................76 9.6. Regionális depozíció................................................................................................78 9.7. Légúti elágazásonkénti depozíció sűrűség...............................................................79 9.8. Lokális depozíció.....................................................................................................80 9.9. Időfüggő részecske kiülepedés................................................................................82 10. KIÜLEPEDETT RADONSZÁRMAZÉKOK LOKÁLIS TERHELÉSÉNEK MODELLEZÉSE URÁNBÁNYÁKNAK ÉS LAKÁSOKNAK MEGFELELŐ EXPOZÍCIÓS KÖRNYEZETBEN................................................................................85 11 KITEKINTÉS.............................................................................................................94 12. ÖSSZEFOGLALÁS………………………………………………………………...96 13. SUMMARY…………………………………………………………………………97 14. IRODALOMJEGYZÉK...........................................................................................98 15. KÖSZÖNETNYILVÁNÍTÁS.................................................................................118 16. FÜGGELÉK.............................................................................................................119 1. A dolgozatban használt rövidítések listája.............................................................119 2. A dolgozatban használt jelölések listája................................................................121

3

1. Bevezetés, célkitűzés Az emberiség radioaktív sugárzások „özönében” él. Bárhol is tartózkodunk, mindig ér minket több-kevesebb természetes, esetleg mesterséges forrásból származó ionizáló sugárzás. Az élő szervezetben lefékeződő ionizáló sugárzás először fizikai és kémiai jelenségeket okoz. Ilyenek az atomok ionizációja, a víz vagy más molekulák kémiai felbomlása, molekulaszerkezeti változások. Ezek a folyamatok gyakorlatilag azonnal lezajlanak. A fizikai és kémiai hatások indukálta biológiai hatások kialakulása azonban időben elnyújtva, órák, napok vagy akár évek múltán figyelhetők meg. Azokat a hatásokat, amelyek többnyire rövid időn belül és vitathatatlanul a kapott sugárterhelés miatt lépnek fel, determinisztikus hatásoknak nevezzük. A determinisztikus hatások egy küszöbdózis felett mindenkinél fellépnek. Ez valószínűleg azért van így, mert egy adott terhelés fölött a szervezet védekezési mechanizmusai kimerülnek. Ezzel ellentétben a sztochasztikus hatásoknál, amelyek a kiváltó sugárterhelés elszenvedése után jóval később lépnek fel, adott dózis esetén megmondható a sztochasztikus effektusok fellépésének valószínűsége, vagy gyakorisága egy nagyobb népességre, de soha nem mondható meg, hogy konkrétan kinél lépett fel a sugárzás miatt az adott hatás. Ezek a hatások ugyanis többlet sugárzásnak nem kitett populációban is gyakran előfordulnak. A sugárvédelem mai gyakorlatában úgy tekintjük, hogy már az átlagos természetes háttérsugárzás (~2,4 mSv, [UNSCEAR 2000]) nagyságánál kisebb többletterhelés is megnövelheti a sztochasztikus hatásokat, a daganatokat és az örökletes károsodásokat. Ugyanakkor a világ számos területén a háttérsugárzás értéke a világátlag tízszeresét sőt néhol akár százszorosát is eléri, de a megbetegedések gyakoriságának valódi növekedését ezidáig nem mutatták ki ezeken a területeken [Sohrabi 1990, Kesavan 1996]. A sztochasztikus hatások tehát tipikusan a kis dózisok tartományában érvényesülnek. Ugyanakkor úgy tűnik, a kis dózisok biológiai hatásának megértése és a lineáris küszöb nélküli dózis-hatás hipotézis (LNT) megtartása vagy elvetése jelenti a jelenkori sugárvédelem legnagyobb kihívását. Ezt látszik alátámasztani az a tény is, hogy az utóbbi évek sugárvédelmi és sugárbiológiai konferenciái többnyire az LNT jegyében zajlanak. A probléma megoldásának társadalmi jelentősége azért is nagy lehet, mert a széles tömegeket érő sugárterhelés ebbe a tartományba esik. A jelenlegi hivatalos

4

sugárvédelmi álláspont szerint a kockázat a kis dózisok tartományában is lineárisan változik a terheléssel, amit manapság sok kutató megkérdőjelez. Ime néhány pro és kontra vélemény, ami jól mutatja a szakmabeliek megosztottságát a kérdésben. Zbigniew Jaworowski, az Egyesült Nemzetek Szervezete Atomi Sugárzások Hatásait Vizsgáló Bizottságának (UNSCEAR) volt elnöke: „Az LNT tudományos képtelenség, a radiofóbia fő kiváltója.” Az Amerikai Egyesült Államok Sugárvédelmi Tanácsának (NCRP) 136-os számú jelentése: „A tudomány mai állása alapján az LNT-nek nincs alternatívája.” Myron Pollycove, az Egyesült Államok Nukleáris Szabályozó Bizottságának (NRC) tagja: „Az LNT nem más, mint egy sok milliárd dolláros önfenntartó politikai és gazdasági erő.” Lássuk, mi lehet az oka - a gazdasági és politikai érdekeken túlmenően - az ellentétes véleményeknek? Mint tudjuk, a sugárforrások feloszthatók természetes és mesterséges forrásokra. A természetes források lehetnek külső- (kozmikus és földkéregből jövő) és belső források (pl. 238U, 232Th leányelemei, 40K, de főleg belélegzett radon bomlástermékek). A legfontosabb mesterséges források az orvosi diagnosztika vagy kezelés során használt röntgensugarak, az atomrobbantások és jóval kisebb mértékben az atomerőművek. A sugárzások hatásának megértéséhez kézenfekvőnek tűnik az álagosnál nagyobb terhelésnek kitett csoportok vizsgálata. Ezt teszi az epidemiológia. A legnagyobb és talán legtöbbet tanulmányozott és vitatott adatbázist a Japán városokra (Hiroshima és Nagasaki) ledobott atombomba túlélőiről évtizedeken át lejegyzett epidemiológiai adatok (ún. élettartam tanulmányok) képezik. Ezek analízise azt mutatja, hogy a magas sugárterhelés következtében először a leukémia, majd másfél évtized után a rákos daganatok is megjelentek [Pierce 1996, Thompson 1994, Preston 1994]. Mivel azonban a túlélők általában nagyobb dózisokat kaptak, a kis dózisok tartományába eső terhelések hatásait inkább csak extrapolálni tudjuk. A kis dózisok esetében a bizonytalanság annyira nagy, hogy amíg ugyanazon nyers adatok analízisével egyesek az LNT-t igazolják [Sinclair 1996], addig mások egyenesen a hormézis elmélete

5

- vagyis, hogy a kis dózisok jótékonyan hatnának az egészségre - fényes bizonyítékának tartják [Cohen 2000]. A hétköznapi kis dózisokra történő extrapolálás azért is megkérdőjelezhető, mert az ionizáló sugárzások biológiai hatása nem csupán a dózistól, de a dózisteljesítménytől is függ [UNSCEAR 1994], ami több nagyságrenddel nagyobb volt az atombombák esetében, mint aminek általában ki vagyunk téve. A módszer alkalmazhatatlanságát jól mutatja, hogy míg a Japán epidemiológiai adatokból extrapolálva több mint 50000 halálesetet jósoltak az 1986-os Csernobili balesetnek köszönhetően az ezt követő 50 évre [Goldman 1987], addíg a nemrég (2005 szeptember 5) a Nemzetközi Atomenergia Ügynökség (IAEA), az Egészségügyi Világszervezet (WHO) és az Egyesült Nemzetek Szervezetének Fejlesztési Programja (UNDP) által közölt átfogó tanulmány szerint ezeknek száma 4000 körüli lesz, amiből az első 20 évben minden valószínűség szerint 60-an haltak meg [IAEA/WHO/UNDP 2005]. Egyéb nukleáris tesztrobbantások (az UNSCEAR több mint 2400-ról tud, főleg a hidegháború idejéből) kevés támpontot szolgáltatnak az alacsony dózisok biológiai hatásának vizsgálatához. Egyrészt mert túlnyomó többségük föld alatt, levegőben vagy lakatlan területen történt és az átlagolt többletterhelés a természetes háttérnek az 1%-át sem teszi ki [Jaworowski 1999], másrészt, ahol emiatt emberek jelentős dózist kaptak (pl. Marshall szigetek) sokáig titkosították az adatokat. Később (a 90-es években) számos ország (Kazahsztán, Marshall Szigetek Köztársaság, Oroszország stb.) kért segítséget a Nemzetközi Atomenergia Ügynökségtől a feltételezetten szennyezett területek radiológiai felméréséhez, de szinte sehol nem mértek feltűnően magas terhelési szinteket lakott területen [Gonzalez 1998]. Ennél sokkal több adat és epidemiológiai tanulmány áll rendelkezésre az egykori uránbányák munkásait ért sugárdózisokról és ezek következményeiről. Már a 16. század elején leírták, hogy szokatlanul sok tüdő-megbetegedéssel kapcsolatos haláleset fordul elő a Közép-Európai Érchegység földalatti bányáiban dolgozók között. Az akkoriban hegyibetegségnek nevezett kór valószínűleg magába foglalta a tüdőrákot, szilikózist és a tuberkulózisos megbetegedéseket egyaránt. Ezzel nagyjából egy időben a németországi Schneeberg környéki bányászok között is szokatlanul magas "tüdőbaj" miatti halálozási arányt figyeltek meg. Csupán a 19. század vége felé sikerült ezt a betegséget a tüdőrákkal azonosítani. A 20. század első évtizedeiben mérték először a radon nagy koncentrációjú

6

jelenlétét a csehországi Jachimov (Joachimstahl) bányáiban. 1932-ben írták le először, hogy a megnövekedett tüdőrákos esetek és a nagy radon koncentráció között összefüggés van [Pirchan és Sikl 1932]. Ennek epidemiológiai bizonyítékát szolgáltatják azok a vizsgálatok, melyeket a mai napig világszerte végeztek a földalatti bányászok különböző csoportjain. Ezek a 20. század utolsó évtizedeiben végzett széleskörű felmérések központi szerepet játszanak azokban a törekvésekben, amelyek a radon és bomlástermékeitől származó kockázatok meghatározására irányulnak. Az egyik legnagyobb ilyen tanulmány az Ionizáló Sugárzások Biológiai Hatásaival foglalkozó amerikai bizottság (BEIR) 1999-ben kiadott közleménye, amely a világ különböző országaiból származó 11, többségében uránbányában végzett részletes felmérés eredményeit kísérli meg egységes szemléletmóddal értelmezni és összefoglalni [BEIR VI 1999]. A tanulmány rámutat, hogy a legnagyobb átlagos sugárterhelést elszenvedett bányászcsoportok között fordult elő a legtöbb tüdőrákos eset. A Nemzetközi Sugárvédelmi Bizottság (ICRP) ezen témával foglalkozó kiadványa is leközölte egy 7 bányára kiterjedő csoportvizsgálat eredményeit. Ebből az derül ki, hogy a radonbomlástermékektől származó halmozott sugárterhelés növekedésével a tüdőrák többlet relatív kockázata monoton növekszik [ICRP 65 1993]. Ugyanakkor az uránbányász adatok kis dózisok tartományára történő extrapolálása megint csak nem járható út, mivel az expozíciós körülmények mások bányára mint például lakásra. Elég ha csak a különböző egyensúlyi tényezőkre, ki nem tapadt hányadra és légzési paraméterekre gondolunk [BEIR VI 1999, NRC 1991]. Fokozottabb sugárterhelésnek vannak kitéve a magas háttérsugárzású helyeken élők, valamint a magas aktivitású építőanyagból készült lakások lakói is. Elméletileg ilyen csoportok vizsgálata lenne ideális a kis dózis dilemma eldöntéséhez, azonban a mai napig végzett epidemiológiai vizsgálatok eredményei nagyon szórnak. Sok a megtévesztő tényező és nehéz elkülöníteni a többletterhelésnek tulajdonítható hatásokat a többitől. Egyéb karcinogének (Ni, Zn, Pb, azbeszt, nem bioszolubilis üveggyapot, ultraibolya sugárzás, dohányfüst, vírusok stb.) jelenléte és szinergikus hatása szintén zavaró. Így aztán nem meglepő, hogy míg több tanulmány is többletkockázatot mutatott ki magas rezidenciális radonkoncentráció esetén [Schoenberg 1990, Field 2000, Ruosteenoja

7

1996], addíg más tanulmányok nem találtak semmilyen összefüggést [Alavanja 1994, Kreienbrock 2001, Blot 1990]. Mint az a fentiekből is látszik, az epidemiológiai tanulmányok bár egy adott érték fölött ok-okozati összefüggést mutattak ki a magas aktivitás-koncentráció és a kóros egészségügyi elváltozások között, nem szolgáltatnak érdemi információt a kis dózisok biológiai hatásának megértéséhez és az LNT dózis-hatás hipotézis érvényességének eldöntéséhez. Erre a megállapításra jutnak az UNSCEAR 2000 jelentés összeállítói, valamint Roger Clarke, az ICRP elnöke is [Clarke 1999, Clarke 2001]. Egy nemrég megjelent, az elmúlt két évtized sugárbiológiai állatkísérleteit átfogó tanulmány rámutat, hogy az eddigi eredmények alapján semmilyen biológiailag plauzibilis dózis-hatás görbe nem zárható ki és egyértelmű dózis-hatás összefüggés megállapítása kis dózisokra gyakorlatilag lehetetlen [Duport 2003]. Ugyanakkor több mint húsz éve a világ számos laboratóriumában intenzív in vitro sejt-besugárzásos kisérletek folynak a biológiai válasz követésének érdekében [Pálfalvi 2003]. A kilencvenes években hat európai laboratóriumban C3H10T1/2 sejtvonalat tettek ki különböző dózisú röntgen sugaraknak. A csoport arra a következtetésre jutott, hogy 200 mGy dózis felett a transzformáció valószínűsége lineáris a dózissal [Mill 1998]. Ezen dózisérték alatti terhelés válaszát több csoport is megkísérelte megállapítani (pl. kromoszóma aberrációkra humán limfocitákban). E mérések analíziséből arra a következtetésre juthatunk, hogy a nagyon kis dózisok tartományában nem valószínű, hogy valaha is megmérhető lesz az igazi válasz [Chadwick 2002]. Mindezen túl az is kérdéses, hogy az in vitro eredmények mennyire tekinthetők relevánsnak a magasabb szintű szerveződésű in vivo környezetben. Az epidemiológia és a mérések mérsékelt sikere a kis dózisok tartományában azt sugallja, hogy csak a karcinogenezis sejtszinten lejátszódó folyamatainak megértése vezethet el a biológiai kockázatok pontos felméréséhez. Ezen a téren a mechanisztikus modellezés, vagyis a biofizikai mechanizmuson alapuló kockázati modellek jelenthetik a jövőt. Az ilyen modellek általában a rákhoz vezető folyamatokat és állapotokat írják le. A radioaktív sugárzás által indukált kóros elváltozások fő stádiumainak az inicializálást, a promóciót és a progressziót tekintik és olyan mechanizmusokat vesznek figyelembe mint például a specifikus és nem specifikus DNS kettős törés, transzlokáció, fixáció, sejt

8

inaktiválás, DNS károsodás kijavítás, apoptózis, mitózis stb [Klebanov 1993, CrawfordBrown 2002]. E modellek bemenő adatai olyan mikrodozimetriai mennyiségek, mint a sejtdózis vagy a találati valószínűség. Ezek ismeretében megmondható, hogy egy sejt milyen valószínűséggel transzformálódik, vagy lesz rákos [Schöllnberger 2002]. Az említett input paramétereket mikrodozimetriai modellekkel lehet kiszámítani. Mivel a természetes sugárterhelés több mint a felét a belélegzett radon bomlástermékeitől kapjuk, ami a tüdőre jelenti a legnagyobb kockázatot, a mikrodozimetriai és rákkeletkezési modellek alkotói is főleg a tüdőre koncentrálnak. A jelenlegi tüdő-depozíciós és dozimetriai modellek azonban a valóságot erősen leegyszerűsítő feltevésekkel élnek a légutak morfológiáját és a radio-aeroszolok légúti kiülepedését tekintve. A légutakat henger alakú egyenes csövekkel írják le és a légutak felületén egyenletes depozíciót feltételeznek, tehát átlag depozíció-sűrűség értékekkel számolnak [Nikezic 2002, Bohm 2003]. Ennek oka pedig az, hogy az analitikus depozíciós modellek bár hatékonyak a teljes és regionális kiülepedés jellemzésében, nem képesek a depozíciót lokálisan kvantálni [Koblinger 1990, Asgharian 2004]. Ugyanakkor, már a korai hisztológiai tanulmányok is kimutatták, hogy a malignáns tumorok eredete néhány sejtnyi lokális mikrokörnyezetre vezethető vissza [Auerbach 1961, Kotin 1959, Macklin 1956]. A légúti depozíció inhomogenitását kísérleti mérések is egyértelműen igazolták [Kinsara 1995, Kim 1999, Oldham 2000]. A fentiekből következik, hogy a rákkeletkezési modellek bemenő adatait képező mikrodozimetriai paraméterek pontos kiszámításához a mai modellek helyett szükség van egy a morfológiailag realisztikus légutak felületi lokális kiülepedését leíró új depozíciós és mikrodozimetriai modellre. Az aeroszolok légúti kiülepedésének inhomogenitását a numerikus áramlástan (CFD) alapú módszerekkel már jó évtizede sikerült kimutatni [Balásházy 1993a,b]. A számítógép kapacitások növekedése és az egyre precízebb numerikus algoritmusok elterjedése további lendületet adott e területnek. Mindezek ellenére a mikrodozimetriában mindmáig nem terjedt el a CFD alapú modellezés. Jelen munkámban egy ilyen modell kifejlesztését és alkalmazását tűztem ki célul. A következő fejezetben bemutatom a numerikus modell felépítése során felmerült problémákat és ezek megoldásait. Ezt követi a légúti levegőáramok és a radionuklid kiülepedés különböző fizikai és biológiai paraméterektől való függésének vizsgálata a modell segítségével. Ezután néhány konkrét

9

példa segítségével bemutatom hogyan alkalmazható az általam felállított modell a tüdőhámsejtek radioaktív terhelésének becslésére. Végül fontos következtetéseket vonok le a makroszkópikus és a sejtszintű terhelés kapcsolatáról, a dózis-hatás görbékről, valamint kitekintést adok a modell távlati szerepéről a sugárterheléshez kapcsolódó kockázat elemzésében. Megmutatom, hogy a modell hogyan kapcsolódhat a jövőben a jelenlegi rákkeletkezési és egyéb kockázati modellekhez, valamint rámutatok a modell továbbfejlesztési lehetőségeire is.

10

I. Rész: A numerikus áramlástani és részecske kiülepedési modell felépítése és validálása 2. A numerikus modellek szükségessége A levegőben található részecskék (aeroszolok) belégzés által a tüdőbe juthatnak. Innen rövid idő alatt kikerülhetnek a kilégzés során, vagy kitapadhatnak a légutak nyákkal borított falára. Az egészségre káros aeroszolok, mint például a gumipor, a fosszilis tüzelőanyag égetéséből származó pernye [Balásházy 2003d], a cigrettafüstben található karcinogén részecskék és sok más toxikus vagy radioaktív anyag a kiülepedés után károsíthatják a tüdőt, vagy a véráramba kerülve az egész emberi szervezetet. Mint azt egy korábbi cikkünkben is kimutattuk [Balásházy 2002a], a légutakba bekerült és ott kiülepedett rövid felezési idejű alfa-bomló radon leányelemek különös veszélyt jelenthetnek a tüdőhámsejtekre. Mivel alfa-bomlás révén kis távolságon nagy energiát adnak le, károsíthatják az epithel sejtmagok DNS-ét és rákot is okozhatnak [Trosko 1996, Farkas 2003c]. A belélegzett radio-aeroszolok biológiai hatásának tanulmányozásához szükség van a részecskék légúti transzportjának és kiülepedésének pontos modellezésére. Mivel a belélegzett részecskék sorsa a légutak geometriája mellett a levegővel való kölcsönhatásuktól is függ, egy olyan kétfázisú áramlási problémával állunk szemben, ahol az egyik fázist a levegő, a másikat pedig a belélegzett részecskék képezik. Analitikus módszerekkel még egyfázisú áramlást is csak a légutak geometriájára és az áramlási profilra vonatkozó erős egyszerűsítésekkel (elágazások helyett egyenes vagy hajlított csövek és egyenletes vagy parabolikus sebességprofil) lehet számítani [Kim 1989, Balásházy 1990]. Az analitikus modellek másik hátránya, hogy nem írják le a részecske kiülepedés légúti elágazásokon belüli inhomogenitását [Asgharian 2004, Cassee 1999, Koblinger 1990]. Márpedig az ionizáló sugárzások egészségügyi hatásának vizsgálata szempontjából a lokális kiülepedés ismerete rendkívüli fontossággal bír [Balásházy 2003a]. A megoldást a numerikus modellek jelenthetik, amelyek a számítógép kapacitások növekedésével egyre hatékonyabbá válnak. A kétfázisú áramlás numerikus leírásához az Euler-Lagrange és az Euler-Euler módszereket lehet használni

11

[Gradon 1996]. Az Euler-Lagrange közelítés esetén a levegőt kontinuumnak tekintjük, benne diszkrét részecskékkel, amelyek pályáját egyénileg követjük. Ezzel szemben az Euler-Euler modell két folytonos keveredett közegként kezeli a két fázist. Természetesen mindkét módszernek vannak előnyei és hátrányai és a konkrét fizikai probléma dönti el, hogy épp melyiket érdemes alkalmaznunk. Mivel a klasszikus Euler-Euler módszerrel nem lehet impakciót (és számos egyéb hatást) modellezni, jelen munkában az EulerLagrange eljárást alkalmaztam. A módszer választását az is indokolta, hogy addig míg Euler-Euler közelítéssel csak egy számítási cella méretű véges felületen kiülepedett részecskék számát kaphatjuk meg, az Euler-Lagrange módszert alkalmazva lehetőség van a részecskék kiülepedési helyének pontos meghatározására is. Ez azért is fontos mert reális gépidők mellett ma még nem lehet néhány tíz sejt felületét kitevő cellánál kisebb tartományban számolni, holott a munka célja a sejtszintű terhelés modellezése. Végül megemlítem, hogy az általam használt FLUENT CFD kereskedelmi kód is az EulerLagrange módszert ajánlja a belélegzett részecskék alacsony térfogatfrakciója esetén. A fenti közelítés legfőbb lépései a háromdimenziós légúti geometriák számítógépes előállítása és berácsozása, a levegő sebességterének a kiszámítása és a belélegzett részecskék transzportjának és kiülepedésének a modellezése. A modellépítés lépéseit valamint a numerikus eljárás validálását a következő alfejezetekben mutatom be.

12

3. Légúti geometriák modellezése Mivel a radon egy belélegezhető gáz és leányelemei közül néhány rövid felezési idejű alfa-bomló izotóp, ezek a legnagyobb kockázatot a tüdőre jelentik [Tóth 1983]. Ezért a kockázatelemző modellek a radonszármazékok légúti transzportját és kiülepedését és a tüdőhámsejtekre gyakorolt hatását tanulmányozzák. A radon leányelemek légutakon belüli pályája és kiülepedése nagymértékben függhet a légutak geometriájától [Balásházy 1996]. Ezért ezeknek minél pontosabb ismerete és modellezése kiemelten fontos feladat.

3.1. A légzőrendszer szerepe és felépítése A légzőrendszer elsődleges élettani szerepe a vér gázcseréjének biztosítása. A gázcsere egy állandó, a szervezet pillanatnyi szükségleteihez alkalmazkodó bonyolult biológiai folyamat. Eredményeként az oxigén a vérbe kerül, hogy aztán a vörösvértestek tovább szállítsák a test minden részébe. Ugyanakkor a test sejtjeinek anyagcsere végtermékeként keletkező széndioxid gáz a légutakon keresztül kikerül a szervezetből [Magyar 1998]. Ezért a szállító rendszernek számos feltételnek kell eleget tennie. Ilyenek például a levegő szállítása során fellépő dinamikai veszteségek minimalizálása, vagy a levegőnek egy igen nagy (kb 60 négyzetméter) felületet befedő komplex hajszálérrendszerbe juttatása. Ezért a légutak szerveződése egy morfológiailag rendkívül bonyolult geometriát eredményez (I.1. ábra). A légutakhoz tartozó szervek az orrüreg, a szájüreg, a garat, a gége, a légcső és a két tüdő, jobbról három, balról két lebennyel [Horváth 1998]. Anatómiailag a légcső előtti légutak a felsőlégutak csoportjába tartoznak és szerepük többek között a levegő szűrése, melegítése és nedvesítése. A légcső két főhörgőre oszlik, majd azok sorra újabb és újabb hörgőkre (bronchus) oszlanak miközben méretük csökken [Kassai 1950]. A légcső és a hörgők rendszerét szokás még tracheobronchiális fának vagy tüdőfának is nevezni. A közel milliméter átmérőjű csövecskéket hörgőcskéknek (bronchiolus) nevezzük. Rajtuk adott generációszám után megjelennek a kiöblösödések (alveolus). Az ilyen csövecskéket bronchioli respiratoriinak nevezzük. Ezek szintén elágaznak ductus alveolisokká, amelyeken már kosárszerűen kiöblösödő kicsiny hólyagocskák találhatók. A levegő az alveolusok falán keresztül

13

diffundál az alveolus körüli kötőszövetet sűrűn beszövő, a kis vérkörhöz tartozó kapilláris rendszerbe, illetve a vérgáz onnan vissza az alveolusba.

garat (pharynx)

orr száj

légcső (trachea)

gége (larynx)

hörgők (bronchusok)

bronchiolus terminalis

bronchiolii respiratorii ductus alveolaris alveolusok

I.1. ábra A légzőrendszer felépítése.

3.2. A tüdőfa matematikai modelljei Mivel a szakirodalmi adatok szerint az uránbányászok esetében a tüdőrák a centrális légutakban keletkezett [Martin 1972, Schlesinger 1978, Churg 1996], elsősorban e régió geometriájának modellezését céloztam meg. A tracheobronchiális fa szerkezetének matematikai leírásával számos publikáció foglalkozik. Ezek egy része a mért adatokra (csőátmérők, hosszak, elágazási és gravitációs szögek) támaszkodik. Az egyik legrégebbi, de még ma is gyakran használatos modell a szimmetrikus, dichotómikus Weibel modell [Weibel 1963]. Később Yeh [Yeh 1980], valamint Phalen [Phalen 1985] is közölte a mért, felnőtt tüdőre vonatkozó geometriai adatokat. Mivel a három publikáció adatai eléggé eltértek egymástól, James [James 1988] átlagolta őket. A tüdőfa e modelljét tekintik referenciamodellnek a mai sugárvédelemben is [ICRP66 1994]. A modell egyik

14

előnye, hogy alkalmazásával a dozimetriai számítások viszonylag leegyszerűsödnek. Azonban a tüdő valós, aszimmetrikus morfológiáját nem írja le. Más szerzők ugyancsak mért adatok alapján próbálnak a tüdő szerkezeti felépítésére vonatkozó törvényszerűségeket levezetni [Koblinger 1985, Philips 1997, Kitaoka 1999]. Ezek már a mért adatok valamilyen szintű statisztikai analízisén alapuló morfometriai összefüggések. Újabban, az orvosi képalkotó módszerekkel nyert képek alapján, a morfológiai paraméterek mérése könnyebb és pontosabb [Sauret 2002]. Egy másik irányzat a Horsfield [Horsfield 1967] által indított, úgynevezett tüdőfa tervezés, ami azt jelenti, hogy biofizikai meggondolások alapján (áramlási rezisztencia minimalizálása, tüdő térfogatának optimalizálása stb.) építenek fel virtuális tüdőfákat, amelyeket aztán összehasonlítanak az igazival [Wechsatol 2004, Mauroy 2003, Kitaoka 1997, Reis 2004]. Ezen irányzatnak nagy lendületet adott a komplex rendszerek fraktál és multifraktál elméletének a tüdőfára adaptálása [Suki 2002, Schlesinger 1991, Sujeer 1997, Barabási 1996, Benett 2003].

3.3. Digitalizált háromdimenziós centrális légúti-geometria modellek A fent leírt egyre fejlettebb tüdőfa modellek nem alkalmasak numerikus áramlástani számítások elvégzésére, legfeljebb a digitalizált modellek elkészítéséhez nyújthatnak némi támpontot. Mivel a későbbiekben alkalmazott számítási modell numerikus fluidumés részecsketranszport számításokon alapszik, a munka egyik fontos eleme a minél realisztikusabb, reprodukálható digitális geometria elkészítése volt. Morfológiai adatokból kiindulva egyedi légúti elágazások generálhatók, amelyek egymáshoz csatolásával tüdőfa szegmensek építhetők. Ezért első lépésben matematikailag egzaktul definiált, reprodukálható és morfológiailag realisztikus elágazásokat hoztam létre. Munkám kezdetekor már léteztek a „narrow” és „wide” nevű, matematikailag is leírt elágazás geometriák (I.2. ábra), azonban ezek az elágazások csúcsainak nyereg alakját (karina régió) nem tükrözték. Egyes szerzők (pl. [Heistracher 1995]) lekerekített karinát alkalmaztak, azonban matematikailag nem írták le az áramlás és részecske-kiülepedés szempontjából kulcsfontosságú régiót.

15

I.2. ábra Háromdimenziós légúti elágazások. a) „narrow”, b) „wide” A Morfológiailag Realisztikus Elágazás (MRB) nevet viselő modellt matematikailag a [Hegedűs 2004] cikkünkben írtuk le részletesen. A geometria struktúráját az I.3. ábra szemlélteti. Az ábrán Rp az anyaág sugarát, Dda és Ddb a leányágak átmérőit, Lp’ az anyaág, L’da és L’db a leányágak hengeres részének hosszát jelöli. Hammersley megfigyelései alapján [Hammersley 1992], jelen modellben a hengeres rész a teljes hossz 80%-át teszi ki mind az anyaág, mind a leányágak esetében. Továbbá, Φa és Φb (ϕsa és ϕsb szaggitális szögek maximumai) a leányágak görbületi szögeit, K pedig a lekerekítési kör középpontját jelöli. Az anya és leányágak hengeres részei közötti átmenetet a BB’D’D és a BB’GG’ gyűrűs csövek képezik. B, B’, D, D’ és G,G’ pontokban a sima átmenetet (C1 kontinuitást: a függvény és az első rendű deriváltja is folytonos) a BD, B’D’, B’G és BG’ görbéket leíró szigmoid függvények biztosítják. A fenti paraméterek megfelelő megválasztásával tetszőleges tüdőelágazás szerkeszthető. Így modellezhetővé válnak az egyének és nemek közti eltérések, valamint figyelembe vehetők a kor szerinti és az egyénen belüli morfológiai különbségek is. A pontos matematikai leíráson túlmenően, szükség van a felületek numerikus létrehozására, hogy később levegő és részecske áramlást modellezhessünk e felületek által definiált háromdimenziós térrészben. Ezt a feladatot két módszerrel is megoldhatjuk, melyek közül az első numerikus pontgeneráláson és háromdimenziós függvényillesztésen alapszik. 16

I.3. ábra Morfológiailag realisztikus elágazás (MRB) fősíkja. A pontgeneráláshoz egy FORTRAN program íródott, melynek bemenő adatait 11 geometriai paraméter képezi. Ezek az anya és leányágak hosszai és átmérői, a leányágak görbületi sugarai és elágazási szögei valamint a karina lekerekítő kör sugara az elágazás fősíkjában. A program a fősíkra merőlegesen megforgat egy egyenest és numerikus gyökkereső eljárással megtalálja, majd fájlba írja, a forgatott egyenes és az elágazásgeometriát határoló felület metszéspontjainak térkoordinátáit. Ezeket az (x,y,z) koordinátákat beolvastam a GAMBIT nevű geometria szerkesztő és hálózó programba, amely a FLUENT CFD kereskedelmi szoftvercsomag része (lásd [GAMBIT Manual 2001]). GAMBIT-es környezetben a beolvasott pontokra háromdimenziós felületet illesztettem, majd létrehoztam a felület által határolt térfogatot, ami már rácsozással számítási tartománnyá alakítható.

17

A második módszer a geometriát leíró egyenleteknek a UNIGRAPHICS CAD/CAM/CAE nevű kereskedelmi szoftverbe való illesztésén alapszik. Mivel az adaptáláshoz az egyenleteket parametrikusan kell megadni, levezettem és egy PERL nyelven írt programmal automatizáltam az eredeti leírásmód és a UNIGRAPHICS által megkövetelt szintaxis közötti átmenetet. Az adaptáció után a UNIGRAPHICS automatikus felületgenerálási eljárással létrehozza a légúti elágazás numerikus felületét illetve térfogatát. Amint azt a I.4. ábra is mutatja, a két módszer ugyanahhoz a geometriához vezet.

I.4. ábra A 4-5. légúti generáció numerikus elágazás-modellje (ha a légcső sorszáma 1) két különböző nézetben. A pontsorok az első, a felület a második módszerrel készült. A két eljárás egyazon felülethez vezetett: a FORTRAN program által generált pontsorok a UNIGRAPHICS által generált felületen nyugszanak. A fenti geometriát validációs és levegőáramlási, valamint részecskekiülepedési számításoknál alkalmaztam, ezért a továbbiakban megadom a morfometriai paramétereit olyan formában, ahogy azok a fent leírt FORTRAN program bemeneti fájljában szerepelnek.

18

Anyaág hossza

(cm)

0,76

A oldali leányág hossza

(cm)

1,27

B oldali leányág hossza

(cm)

1,27

Anyaág átmérője

(cm)

0,56

A oldali leányág átmérője

(cm)

0,45

B oldali leányág átmérője

(cm)

0,45

A oldali elágazási szög

(°)

35

B oldali elágazási szög

(°)

35

A oldali görbületi sugár

(cm)

1,42

B oldali görbületi sugár

(cm)

1,42

Lekerekítési görbületi sugár

(cm)

0,5

1

9 6

8 7

3

2

5

11 10

4

1: légcső (trachea) 2: bal főhörgő 3: jobb főhörgő 4: bronchus intermedius, a jobb főhörgő folytatása, 1: légcső a jobb-közép és jobb-alsó lebenyhörgők anyaága 2: bal főhörgő 5: lebenyhörgő 3: jobb-felső jobb főhörgő 4:7: jobb-közép lebenyhörgő 6, szegmentális hörgők, a jobb-felső lebenyhörgő 5: jobb-felsőleányágai lebenyhörgő 6, 7: szegmentális hörgők 8, 9, 10, 11: szubszegmentális hörgők, a 8, 9, 10, 11: szubszegmentális szegmentális hörgők leányágai hörgők

I.5. ábra Numerikus ötgenerációs centrális légúti szegmens. A leírt két módszer bármelyikét alkalmazva bonyolultabb légúti szegmensek is előállíthatók. A szakirodalmi adatok alapján [Martin 1972, Schlesinger 1978, Churg

19

1996] az uránbányászok esetében a tüdőrák a legnagyobb gyakorisággal a jobb felső lebeny 3-5. generációiban keletkezett, ezért munkám során elsősorban ezt a légúti szegmenst tanulmányoztam. Az I.5. ábra a jobb felső lebenyhez vezető ötgenerációs fát mutatja, amely magába foglalja a légcsövet, a főhörgőket, a jobb felső lebenyhörgőt, a bronchus intermediust, valamint a jobb felső lebenyhörgőből ágazó szegmentális és szubszegmentális hörgőket. A geometria paramétereit a sztochasztikus tüdőmodell segítségével generáltam [Koblinger 1990]. A modell a tüdő morfológiai adatait tartalmazó

Lovelace

adatbázis

[Raabe

1976]

statisztikai

eloszlásaiból

sorsol

csőhosszakat, átmérőket, elágazási és gravitációs szögeket. Az I.5. ábrán szereplő tracheobronchiális légúti geometria adatait az I.1. táblázat tartalmazza. I.1. táblázat Az öt elágazásegység geometriai adatai. A csövek sorszámai megegyeznek az I.5. ábrán szemléltetettekkel, a generációs felosztást pedig az I.6. ábra mutatja. gen 1-2

gen 2-3

gen 3-4

gen 4-5.1

gen 4-5.2

Anyaág hossza (cm)

9,68

1,22

0,52

0,72

0,7

A oldali leányág hossza (cm)

1,22

1,24

0,72

0,59

0,61

B oldali leányág hossza (cm)

2,15

0,49

0,7

0,57

0,59

Anyaág átmérője (cm)

1,61

1,36

0,94

0,69

0,83

A oldali leányág átmérője (cm)

1,36

1,14

0,83

0,63

0,69

B oldali leányág átmérője (cm)

1,1

0,94

0,69

0,51

0,63

A oldali elágazási szög (°)

35

5

45

5

45

B oldali elágazási szög (°)

50

85

50

45

45

A oldali görbületi sugár (cm)

35,42

24

1,6

1,6

1,6

B oldali görbületi sugár (cm)

35,42

1,2

1,6

2

1,6

Karina görbületi sugara (cm)

0,01

0,1

0,1

0,1

0,3

Cső sorszáma (I.5. ábra) Gravitációs szög

1

2

3

4

5

6

7

8

9

10

11

0

50

35

30

120

108,7

110,7

113,3

67,3

133,1

79,5

20

gen 1-2

gen 4-5.2 gen 3-4

gen 4-5.1 gen 2-3 I.6. ábra Az I.5. ábrán látható geometria légúti generációs felosztása. A fenti két légúti geometrián kívül számos más centrális légúti elágazásban és elágazásrendszerben számoltam. Közülük néhányat a I.7. ábra szemléltet.

I.7. ábra Centrális légúti elágazás-geometriák.

21

A légúti geometriák modellezésében egy új irányzat van kibontakozóban. Ez a háromdimenziós geometria rekonstrukcióját jelenti orvosi képalkotó eszközök, komputer tomográf (CT), mágneses magrezonanciás készülék (MRI) stb. által készített kétdimenziós szeletekből [Szőke 2004a]. Mivel a képpontok árnyalata a szövet, csont, levegő stb. sűrűségétől függ, az intenzitások analízise révén lehetőség van a légutak szegmentálására [Perlz 1996]. Megfelelő programmal a légutakat határoló görbékből háromdimenziós felület vagy térfogat generálható (I.8. ábra) [Moody 1997]. Munkám során e célra a 3D-DOCTOR és a SURFDRIVER programokat alkalmaztam. A jelenlegi legjobb CT készülékek felbontásával is csak a centrális légutak első három-négy generációja rekonstruálható élethűen. A módszer leginkább a felsőlégutak esetében hatékony, ahol a méretek nagyobbak. A módszer kidolgozása folyamatban van (lásd [Szőke 2002d,f], a közeljövőben reményeink szerint már rekonstruált geometriákban is számolhatunk.

I.8. ábra Tracheobronchiális fa rekonstrukciója CT felvételekből. Fent egy kétdimenziós szelet (balról) és annak szegmentációja (kinagyítva jobbról) látható. A légutak kétdimenziós metszetei lent balról, míg a háromdimenziós rekonstruált felület lent jobbról látható. 22

3.4. A felsőlégutak modellezése Mint említettem, radon leányelemek dozimetriája szempontjából elsősorban a bronchiális régió modellezése fontos. Azonban ahhoz, hogy e régióban pontosan modellezzük a radonszármazékok kiülepedését, szükség van a belégzéskor a felsőlégutakból a légcsőbe jutó, illetve kilégzéskor a tüdő perifériájáról visszajövő részecskék számának ismeretére. Jelenleg ez például a sztochasztikus tüdőmodellel számítható ki, de e területen is a numerikus áramlástani modelleké a jövő. Felsőlégúti számításaimhoz a Cheng által leírt [Cheng 1999] száj-garat-gége modell leegyszerűsített változatát használtam, ami az I.9. ábrán látható. Ugyancsak az I.9. ábra szemlélteti a CT felvételekből rekonstruált háromdimenziós orrgarat-gége modellt. Ez utóbbiban számításokat még nem végeztem, ugyanis számítási tartománnyá alakítása még csak most van folyamatban.

I.9. ábra Numerikus száj-garat-gége (balról) és orr-garat-gége (jobbról) modellek. A bemutatott felső és centrális légúti geometriákon kívül kidolgoztam egy alveolus [Farkas 2004a,b,c], meg több beteg tüdő (tumoros, beszűkült, elzáródott) [Farkas 2001b, 2005a] modellt is. Ezek a geometriák fontos szerepet játszanak pl. a toxikus aeroszolok vagy aeroszol gyógyszerek [Farkas 2005b] kiülepedésének modellezésében, de kevésbé

23

jelentősek a radio-aeroszolok kiülepedése és hatása tanulmányozásának szempontjából. Ezért ezek bemutatásától itt eltekintek.

24

4. A numerikus rács A 3.3. fejezetben bemutatott légúti geometriák hálózással számítási tartománnyá alakíthatók. Ez gyakorlatilag azt jelenti, hogy a légtereket diszkrét tartományokra, cellákra osztjuk. A jó rács elkészítése elengedhetetlen előfeltétele a pontos számításoknak. Ezért, e feladatra különös figyelmet fordítottam. Rácsozáshoz a GAMBIT nevű CFD pre-processzor programot használtam, amely egyaránt alkalmas strukturált és nem strukturált matematikai hálók generálására. A strukturált és nem strukturált háló közötti különbség az, hogy addig míg a strukturált háló minden csomópontjához szükségszerűen egy i, j, k (két dimenzióban i, j) index tartozik, addíg a strukturálatlan hálónál nincsenek tengelyek, ezért topológiailag rugalmasabb. Ugyanakkor a strukturálatlan hálónál megengedett az ún. függő csomópontok jelenléte is (lásd I.10. ábra).

A (i,j)

F

(2,3 ) (2,2 (2,1) (1,3 (1,2 ) ) (1,1) ) F - függő csomópont I.10. ábra Struktúrált (bal oldal) és struktúrálatlan (jobb oldal) hálók. A FLUENT, mint egy véges térfogat módszert alkalmazó kód, nem követel meg semmilyen rácsstruktúrát vagy topológiát [Fluent Manual 2001]. Ezért célszerűnek láttam a nagyobb flexibilitást nyújtó strukturálatlan háló alkalmazását. Bár strukturálatlan háló

25

hasáb alakú cellákból is létrehozható, a geometria megfelelő altartományokra való osztásával, a leggyakrabban tetraéderes rácsot használtam. Cellatípuson belül a cellák alakja is nagymértékben befolyásolja a számítások pontosságát. A nagyon lapos vagy ferde cellák interpolációs hibákat hozhatnak be, ezért minden esetben ajánlatos ezek mellőzése. Általános szabályként arra törekedtem, hogy a cellák egyetlen éle se legyen a többi él egyötödénél rövidebb, illetve az élek által bezárt szög 60 fokhoz közeli legyen, de semmiképp se legyen nagyobb mint 90 fok. Emellett hangsúlyt fektettem arra is, hogy a szomszédos cellák térfogatai közötti különbség ne legyen nagy. A rácsgenerálás másik kulcskérdése a cellák száma. A levegőáramlás és részecsketranszport leírásának pontossága általában nő a térfelbontással, vagyis a numerikus cellák számával. Ugyanakkor a gépidő is a számítási kontrollpontok (csomópontok, vagy cella közepek) számával nő. A jelenlegi számítógép kapacitások lehetővé tették, hogy sűrűbb háló alkalmazásával az áramlási térnek olyan részletei is vizsgálhatók legyenek, amire a korábbi ritka háló nem adott lehetőséget. Az I.11. ábrán sebesség szintvonalas reprezentáció látható egy légúti elágazás fősíkjában homogén (bal oldal) illetve inhomogén (jobb oldal) háló esetén. A szintvonalas ábrázolás esetén folytonos vonallal kötjük össze az azonos sebességértékkel jellemzett pontokat. A szomszédos szintvonalaknak megfelelő sebességkülönbség állandó. Ezért a sűrű szintvonalakkal jellemzett régiókban a sebességgradiens nagy, míg a ritkább szintvonalú helyeken a sebességprofil laposabb. Látható, hogy a sebességtér mennyire érzékeny a rács finomságára.

I.11. ábra Rács sűrűségének hatása az áramlási képre. Bal oldal: ritka rács (30 000 cella), jobb oldal: sűrű rács (800 000 cella). 26

A számítási kapacitások végesek. A fenti cellasűrűséget használva egy öt-hat generációt magába foglaló légúti szegmens esetén a számítások akár évekig is eltartanának. Ezért meg kell találni egy optimális cellaszámot, ami már matematikailag és fizikailag elfogadható pontosságú eredményt produkál, de a szükséges gépidő nem irreálisan nagy. Ez hatékonyan úgy valósítható meg, hogy sűrűbb hálót készítünk azokon a helyeken, ahol a sebességgradiens nagy és durvább rácsot ott, ahol a sebesség közel állandó. Mivel ez feltételezi a sebességtér előzetes ismeretét, mérési és más számítási tapasztalatokra támaszkodhatunk, vagy az is lehetséges, hogy durva hálóval számolunk, majd visszatérünk a hálózáshoz és finomítjuk ahol szükség. Sebességgradiens szerinti automatizált hálóadaptációra a FLUENT-ben is lehetőség van, de mivel az eljárással szemben kifogásaim voltak, inkább saját hálófinomítást végeztem. Az automatikus adaptációnál a program a meglévő durva hálót úgy finomítja, hogy csak a számítási tartomány belsejébe tesz új csomópontokat és közben a durva rácsozás miatti eredeti geometriától való eltérést nem orvosolja (I.12. ábra, b panel). A „size function” technikát alkalmazva a GAMBIT-ben viszont már eleve olyan hálót generáltam, amely a kritikus helyeken sűrű és optimálisan tölti ki a légutak belsejét [Gambit Manual 2001]. Ugyancsak

célszerű

sűrű

rácsot

használni

azokon

a

helyeken,

ahol

a

részecskekiülepedés várhatóan intenzívebb lesz. Mivel mind a kisérletek [Oldham 2000, Miguel 2004, Kinsara 1995], mind pedig az eddigi modellszámítások [Gradon 1990, Zhang 2002, Broday 2004] azt mutatják, hogy az elágazások csúcsában, az ún. karina régióban a kiülepedési sűrűség nagy, lehetőség szerint a hálót itt is besűrítettem (I.13. ábra). Hálózáskor érdemes még a numerikus diffúziót is figyelembe venni. A kontinuumról diszkrét tartományokra való áttérés elkerülhetetlenül numerikus hibához vezet [Hirsch 1988]. Mivel ez úgy nyilvánul meg, mintha a valós diffúzió megnőne, szokás még hamis diffúziónak is nevezni. A konvekció dominált áramlásnál a szerepe megnő. Nagyságát csökkenteni a rács finomításával, illetve másodrendű sémák alkalmazásával lehet [Fluent Manual 2001]. Ezért munkám során kizárólag ilyen algoritmusokat használtam.

27

légút fala

a)

b)

c)

I.12. ábra Egy légúti elágazás rácsozott leányágának keresztmetszete: a) homogén rács, b) automatikusan adaptált rács, c) „size function” technikával készült rács. Fent kinagyítva a háló peremhez közeli része látható.

lokálisan (karina régió) sűrű rács

légút falának közelében sűrű rács

I.13. ábra Inhomogén matematikai rács a centrális légutak egy ötgenerációs szegmensén.

28

5. A levegő áramlási terének numerikus áramlástani modellje 5.1. Megmaradási törvények A radio-aeroszolokat szállító levegő áramlási terének meghatározása matematikailag a megmaradási törvényeket leíró egyenletek megoldását jelenti. Ezek elsősorban a tömeg és az impulzus megmaradás törvényei. Hőtranszport vagy összenyomható fluidum esetén szükség van még az energia megmaradás törvényére és turbulencia esetén plusz transzport egyenletekre is [Chen 1998, Davidson 2003]. 5.1.1. Tömegmegmaradás (folytonosság) egyenlete A tömegmegmaradás vagy kontinuitás egyenletének differenciális alakja a következőképpen írható:

∂ρ + ∇ ⋅ ( ρv) = S m ∂t

(5.1)

ahol ρ a sűrűséget, t az időt, v a sebességvektort, Sm pedig a forrástagot jelöli. A törvény érvényes összenyomható és összenyomhatatlan fluidumokra egyaránt. A forrástag a fluidumhoz hozzáadott, vagy elvett (pl. párolgás, lecsapódás által) tömegmennyiséget jelenti. Ha nincsenek források meg nyelők és a sűrűség időben és térben állandó, akkor a ∇⋅v = 0

(5.2)

lesz a folytonossági egyenlet. 5.1.2. Impulzusmegmaradás egyenlete (mozgásegyenlet) Az impulzus megmaradását egy inerciarendszerben a ∂ ( ρv) + ∇( ρv ⋅ v) = −∇p + ∇ ⋅ (τ ) + ρg ∂t 29

(5.3)

egyenlet írja le, ahol p a statikus nyomás, ρg pedig a gravitációs erő. Newtoni fluidumra (amilyen a levegő is) a τ feszültség tenzor komponenseit a

τ ij = µ (

∂u i ∂u j 2 ∂u + ) − µ k δ ij ∂x j ∂xi 3 ∂xk

(5.4)

képlet adja meg, ahol µ a viszkozitás, ui , uj , uk és xi , xj , xk a sebesség és koordináta komponensek, δij pedig a Kronecker szimbólum. A fenti képletben a nabla formalizmus helyett indexes írásmódot illetve az Einstein-konvenciót használtam. A továbbiakban mindkét írásmódot használom. Ha a kőzeg összenyomhatatlan, akkor a mozgásegyenlet a

µ ∂v 1 + ( v ⋅ ∇) v = g − ∇p + ∆v ∂t ρ ρ

(5.5)

alakot veszi fel, amit a szakirodalomban Navier-Stokes egyenletnek is neveznek [Landau 1971]. A µ/ρ arányt kinematikai viszkozitásnak nevezzük. Ha a közeg súrlódásmentes is akkor az (5.5) jobboldali utolsó tagja hiányzik és az Euler egyenlethez jutunk, ez utóbbit azonban a számítások során nem használtam, így be sem mutatom. Ezenkívül, ha az áramlás stacionárius, akkor a fenti egyenletekből az idő szerinti parciális deriváltak hiányozni fognak. A dolgozatban stacionárius és időfüggő számításokat egyaránt végeztem.

5.2. Az áramlási tér meghatározása Az impulzus egyenlet három skaláris egyenletet jelent és így a folytonossági egyenlettel négy egyenletből álló rendszert alkot. A sűrűséget időben és térben állandónak tekintve az áramlási teret jellemző másik négy mennyiség (három sebességkomponens és a nyomás) meghatározható [Lajos 2000]. A valóságban minden áramlási tér végtelen belső pontot tartalmaz, tehát a négy mennyiség megadása a tér minden pontjában lehetetlen. Elvileg erre nem is lenne szükség, ha minden esetben

30

találnánk analitikus megoldást, de sajnos ez olyan bonyolult geometriákra mint például a légúti elágazások nem lehetséges. Marad tehát az a megoldás, hogy a teret diszkretizálva minél több pontban megadjuk a sebesség és nyomás értékeket, a meg nem adott pontokban pedig, ha szükség van rájuk, az ismert értékekből kiindulva valamilyen módszerrel megbecsüljük. Erre a célra numerikus áramlástani (CFD) eljárások alkalmasak. A ma már klasszikusnak számító véges differencia, véges elem és véges térfogat módszereken kívül megjelentek alternatív numerikus számítási modellek is. Ilyen például a numerikus rácsot nem igénylő simított részecske modell (SPH) [Benz 1990, Monaghan 1992], a spektrális módszerek (SM) [Canuto 1988, Peyret 2002], a rács Boltzmann módszer (LBM) [Chen 1992] vagy a spektrális és véges elem módszereket ötvöző spektrális elem módszer (SEM) [Bernardi 1996]. E módszereket alkalmazták már az asztrofizikában, sőt napi mérnöki áramlástani problémák megoldásában is. Mivel a FLUENT nevű véges térfogat kód adott volt, számításaimhoz is ezt használtam. A három alapmódszer (véges differencia, véges elem és véges térfogat) történetét, matematikai leírását és összehasonlítását számos szakkönyv (pl. [Roache 1972, Baker 1983, Hirsch 1988]) tartalmazza, ezért itt csak az általam alkalmazott véges térfogat módszert ismertetem röviden. A véges térfogat módszer lényege abban áll, hogy miután a fizikai teret véges tartományokra (számítási cellákra) osztottuk, a megmaradási törvények egyenleteit minden cellára felírjuk, majd ezek integrál alakjait diszkretizáljuk és a kapott egyenleteket megoldjuk [Versteg 1995]. A módszer előnye a véges differencia módszerhez képest, hogy itt a celláknak tetszőleges alakja és elrendeződése lehet [Vinokur 1989]. Az áramlási tér jellemző mennyiségeit cella közepében vagy a csúcspontjaiban tárolhatjuk. Mivel az általam használt tetraéder alakú cellák száma általában sokkal nagyobb volt mint a csomópontoké, ugyanannyi számítási befektetéssel jobb felbontást értem el ha a cella-közép alapú tárolást használtam.

31

5.2.1. A transzportegyenlet diszkretizálása Mivel a megmaradási tételeket kifejező tömeg és impulzus transzport egyenletek egységes formába hozhatók, a továbbiakban csak egyetlen Φ skalár transzportját leíró egyenletet használok. A Φ megfelelő helyettesítésével visszakaphatjuk a tömegre és impulzusra vonatkozó egyenleteket. Az egyesített transzportegyenlet differenciális alakja ∂ ∂ ∂Φ ( ρui Φ ) = (Γ ) + SΦ ∂xi ∂xi ∂xi

(5.6)

lesz, ahol Γ diffúziós együttható, SΦ pedig Φ forrása az egységnyi térfogatban. A bal oldali kifejezést a transzportegyenlet konvektív, a jobb oldali első tagot pedig az egyenlet diffúziós tagjának nevezzük. A fenti egyenletet véges V térfogatra integrálva és diszkretizálva a

∑J f

f

Φ f = ∑ D f + SΦV

(5.7)

f

egyenlethez jutunk [Mathur 1997], ahol Jf a tömegáram, Φf a Φ mennyiség értéke a cella határán, Df pedig a diffúziós tag, szintén a cella határán. Az összegzések a cella lapjaira (két dimenzióban a cella éleire) vonatkoznak. Látható, hogy az egyenletben a cella határaira jellemző értékek szerepelnek. Mivel az áramlási mennyiségeknek csak a cellaközépre jellemző értékei ismertek, a határokon valamilyen módszerrel ki kell számítani őket. A fenti egyenletek minden véges térfogat módszer esetében érvényesek. Az említett mennyiségek cellahatár értékeinek a kiszámítására már sokféle eljárás létezik és általában ebben különböznek egymástól a véges térfogat módszerek. A FLUENT CFD kódba több eljárást is beépítettek, a felhasználóra bízva az adott feladat megoldásához legalkalmasabb módszer kiválasztását. Mivel néhány megoldás eléggé speciális és a legmegfelelőbb eljárás kiválasztása sem triviális, szükségesnek látom röviden bemutatni az általam alkalmazottakat. A klasszikus, vagy jól dokumentált eljárásoknál a részletektől természetesen eltekintek.

32

5.2.2. A konvekciós tag kiszámítása A Jf tömegáramot egy adott iterációs lépésnél ismertnek tekinthetjük. Ezért a konvektív fluxus meghatározása a Φf értékek meghatározására redukálódik. Ezt sok esetben a szélirányú cellaértékkel helyettesítik [Brooks 1982]. Szélirányú cellának azt a cellát nevezzük, amelyből a sebességvektor az aktuális cella felé mutat (lásd az I.14. ábrát).

a

f Φszi et es

dr

C1

v

C0

b

I.14. ábra Két szomszédos számítási cella; f – a két cella határlapja (két dimenzióban határél); C0, C1 – a két cella közepe; dr - szélirányú cella közepét a határlap közepével összekötő helyzetvektor; es , et – C0C1 és ba irányokba mutató egységvektorok; v – sebesség a cellák határán; Φszi – a Φ értéke a szélirányú cellában. Tehát ezen elsőrendű közelítésnél az I.14. ábra szerinti C0 és C1 középpontú cellák határán a

Φ f = Φ szi

(5.8)

közelítés érvényes. Ez a módszer viszont a legtöbb esetben jelentős numerikus hibát hoz be. Ezért számításaimban kizárólag magasabb rendű közelítést alkalmaztam. Ilyen eljárás például a rekonstruált gradiensen alapuló másodrendű szélirányú séma. Ebben az esetben a Φ értékét a cella határán a 33

Φ f = Φ szi + ∇Φ szi dr

(5.9)

képlet adja meg, ahol dr a szélirányú cella közepét a határlap közepével összekötő vektor, ∇Φ szi pedig a szélirányú cellára jellemző rekonstruált gradiens. A strukturált hálókkal ellentétben itt a gradiens kiszámítása nem egyértelmű. Értékét lehetőleg olyan módszerrel kell megbecsülni, ami független a cella alakjától. Ilyen eljárás a divergencia tételen alapuló gradiens rekonstrukció [Hyman 1992], aminek értelmében

∇Φ szi =

α0 V

∑ (Φ

f

A) ,

(5.10)

f

ahol az összegzés a cella összes lapja szerint van és a Φ f értékeit a szomszédos cellák Φ értékeinek átlagolásából kapjuk. Az α0 tényezőt úgy határozzuk meg, hogy a csomópontokra vetített Φ értékek ne hozzanak be új lokális szélsőértékeket [Venkatakrishnan 1993]. Tapasztalataim szerint, a másodrendű közelítés lassítja a konvergenciát. Ezért bizonyos esetekben elsőrendű közelítéssel indítottam, majd miután az bekonvergált, átkapcsoltam másodrendűre. 5.2.3. A diffúziós tag kiszámítása A diffúziós tag az f cellalapnál D f = Γ f ∇Φ ⋅ A ,

(5.11)

ahol A a lap területvektorát jelöli (kétdimenziós esetben az él hosszvektorát). Df -et a C0 és C1 cellákra jellemző Φ értékek segítségével fejezzük ki. Ennek érdekében a ∇Φ ⋅ A szorzatot a határlap középpontjából kiinduló és a cellák középpontjait összekötő

C0C1 valamint az f lap mentén húzott egységvektorokkal írjuk fel (I.14. ábra):

34

∇Φ ⋅ A =

Φ1 − Φ 0 A ⋅ A Φ a − Φ b A ⋅ A et ⋅ e s . + ds A ⋅ e s A A ⋅ es

(5.12)

A fenti kifejezés hasonló, mint amit a strukturált hálókra is használnak [Peric 1985], de a Φa és Φb meghatározása, a kétdimenziós esetet leszámítva [Jiang 1994], strukturálatlan hálókra már nem egyértelmű. Ezért a diffúziós tagot a

D f = Γf

A⋅A Φ1 − Φ 0 A ⋅ A + Γ f (∇Φ ⋅ A ) − ∇Φ ⋅ e s ) ds A ⋅ e s A ⋅ es

(5.13)

képlettel közelítjük, ahol ∇Φ értékét a szomszédos celláknak a bemutatott eljárás szerint rekonstruált ∇Φ értékeiből átlagoljuk. 5.2.4. Peremfeltételek Az áramlási tér meghatározásához szükség van a peremfeltételek ismeretére. Ezeket a felhasználó adja meg. A peremcellák esetében a cella fali határán is tároljuk Φ értékeit. A megfelelő diffúziós fluxus itt

Db = Γb

Φb − Φ0 A ⋅ A A⋅A + Γb (∇Φ 0 ⋅ A − ∇Φ 0 ⋅ e b ) db A ⋅ e b A ⋅ eb

(5.14)

lesz. Itt Φb a Φ fali értékét jelöli, db pedig a távolság a cella középpontjától a fali lap középpontjáig (I.15. ábra). A FLUENT CFD kódban számos peremfeltétel közül választhatunk. Jelen számításoknál a légutak falánál nulla levegő sebességértéket írtam elő. A bemeneten sebességprofilt adtam meg, amit általában parabolikusnak vettem. Stacionárius áramlásnál a kimeneti peremeket „pressure outlet”-nek vagy „outflow”-nak definiáltam. A „pressure outlet” tipusú peremnél a Gauge típusú nyomást (statikus túlnyomást) kell megadni. Mivel erre nincs mért adat, ezt nullának vesszük minden kimeneten. Ha figyelembe vesszük, hogy a légúti levegőáramokat a nyomáskülönbségek vezérlik és néhány generáción át a nyomásesés kicsi [Pedley 1977] akkor ez a feltétel

35

elfogadhatónak tűnik. Ha a geometria aszimmetrikus és több elágazásból áll akkor ez már egy erősen durva közelítés. Ilyenkor jobb „outflow” peremet venni, amelynél a kimenetre merőleges gradiensek nullák, kivéve a nyomásgradienst. Itt viszont meg kell adni a térfogatáramot minden egyes kimenetre, amit sokszor nem tudhatunk előre. Számításaim során azt vettem alapul, hogy a tüdő különböző részeinek ventilációja azonos kell, hogy legyen. Ennek fényében feltételezhető, hogy az egyes ágakban a térfogatáramok arányosak a mögöttük lévő ventilált térfogattal. A megfelelő térfogatértékeket a szakirodalomból vettem [Horsfield 1971]. Időfüggő áramlásnál a bemeneten időfüggő sebességprofilt, a kimeneten túlnyomást adtam meg.

db c0

eb

b

I.15. ábra Numerikus cella a számítási tartomány határán; C0 – a cella középpontja; b – fali lap (él) közepe; db – C0 és b közötti távolság 5.2.5. Lineáris algebrai egyenletek és megoldásuk A cellák középpontjaiban tárolt Φ változóra a diszkretizációs eljárás a következő lineáris egyenletrendszert adja: a p Φ p = ∑ anbΦ nb + S p

(5.15)

nb

Itt az összegzés a p cella összes nb szomszédja szerint történik. Az Sp forrástag magába foglalja Φ összes térfogati forrását, a konvektív fluxus másodrendű tagjait, valamint a másodrendű diffúziós fluxusokat. Az egyenletrendszer iteratív megoldásához a konvergencia gyorsítása érdekében a hagyományos Jacobi módszer [Demidovich 1981] helyett egy lokálisan simító Gauss-Seidel [Dringó 1992] és egy nagyobb skálán relaxáló 36

multigrid eljárást [Hutchinson 1986] alkalmaztam. A sebesség-nyomás csatolást egy fél-

Skálázott reziduálisok

implicit algoritmus (SIMPLE) biztosította [Patankar 1980]. — x sebesség — y sebesség — z sebesség — folytonosság

Iteráció száma

I.16. ábra Skálázott reziduálisok A konvergencia követésekor a skálázott reziduálisok viselkedését vettem alapul amelyeket valamely mennyiség két egymás utáni iterációs lépésben számított relatív különbségének összes cellára vett összegeként értelmezünk [Ralston 1969]. Általában egy futást akkor tekintettem konvergáltnak, ha a skálázott reziduálisok 10-5 - 10-6 alá estek. Egy tipikusan jól konvergáló számítás reziduálisait az I.16. ábra mutatja. A FLUENT tartalmaz szegregált és csatolt megoldó modult is. A szegregált eljárás esetén a transzport egyenleteket külön-külön (szegregálva) oldjuk meg, csatolt eljárásnál pedig együtt. A csatolt eljárás alkalmasabb összenyomható és nagy sebességű áramlások esetén. Összenyomhatatlan fluidumra a két eljárás egyformán jónak bizonyult. Mivel a csatolt megoldó memóriaigénye a szegregálténak közel a kétszerese, választásom ez utóbbira esett. A szegregált eljárás vázlatát az I.17. ábra szemlélteti. Mint látható, a módszer első lépésként frissíti az áramlási tér jellemzőit. Ha csak az első iterációs eljárásnál tartunk, akkor ezeket a megadott kezdeti értékekből számítja ki. Utána a három impulzus egyenletet oldja meg a nyomás és a tömegáram cellalapi aktuális értékeinek felhasználásával. Az így kapott sebességértékek azonban lokálisan nem mindig teljesítik a kontinuitás feltételt. Ezért, a folytonosság nyomás korrekciós egyenleteit oldjuk meg a sima kontinuitási egyenlet helyett [Rhie 1983]. A korrekciót alkalmazva új nyomás,

37

sebesség és tömegáram értékeket kapunk. Utolsó lépésként a konvergencia kritériumok teljesítését ellenőrizzük és ha ezek nem teljesülnek egy újabb iterációs ciklus következik. Áramlási tér mennyiségeinek frissítése

Impulzus egyenletek megoldása

Nyomás korrekciós, (folytonossági) egyenlet megoldása

Konvergált?

STOP

I.17. ábra A szegregált módszer áttekintése 5.2.6. Idődiszkretizálás instacioner áramlás esetén Ha az áramlás időfüggő, az ismeretleneknek az idő szerinti deriváltja is megjelenik

∂Φ = F (Φ ) . ∂t

(5.16)

Az egyenletben F(Φ) magába foglalhat bármilyen térdiszkretizációt. Elsőrendű „backward” differencia közelítést használva [Hirsch 1998] a fenti egyenlet Φ n +1 = Φ n + ∆tF (Φ n +1 ) alakba írható. Másodrendű közelítéssel pedig a

38

(5.17)

3Φ n +1 − 4Φ n + Φ n −1 = F (Φ n +1 ) 2 ∆t

(5.18)

sémát használjuk [Collatz 1960]. Ezek implicit közelítések, mert Φn+1 értéke egy cellában a szomszédos cellák Φn+1 értékeitől is függ. A munka során az elsőrendű közelítés is elégségesnek bizonyult és mivel ennek jóval kisebb a gépidő igénye, ezt a módszert alkalmaztam. Ez gyakorlatilag a Φ i = Φ n + ∆tF (Φ i )

(5.19)

képlet szerinti iterálást jelenti miután a Φi - t Φn – re inicializáltuk, Φn+1 pedig a már konvergált Φi értékét veszi fel. 5.2.7. A turbulencia kérdése

A légúti levegőáramok lehetnek laminárisak, tranziensek vagy turbulensek, annak függvényében, hogy milyen légzési paraméterekről és melyik légúti szegmensről van szó. A felsőlégutakban a viszonylag nagy sebességek és az összetett geometria miatt az áramlás szinte mindig turbulens [Dekker 1961, Katz 1999]. A bronchiólusokban, a levegő nagyon lassú, ezért az áramlás lamináris, míg a centrális légutak első néhány generációjában nehéz fizikai munkának megfelelő légzésre szintén turbulens, ülő vagy fekvő helyzetnek megfelelő légzés mellett átmeneti vagy lamináris [Finlay 1996]. Értékelhető mérések hiányában csak elméleti megfontolásokat követhetünk. A turbulencia kihalásának egyik oka lehet a légutaknak, mint konvergáló csöveknek (konfúzoroknak) a laminarizáló hatása [Braga 2004]. Másrészt, az összenyomhatatlan fluidumra felírt (5.2) egyenletből és a morfometriai adatokból egyenesen következik, hogy a leányágakra jellemző Reynolds számok kisebbek mint az anyaágé. A turbulenciának tehát valahol ki kell halnia, kérdés, hogy hol? Ha a csőátmérőket mint a generációszám függvénye vizsgáljuk, akkor a hatványfüggvényt mutató görbének valahol a 4. generáció környékén van egy törése (I.18. ábra). A görbéhez ebben a pontban húzott két érintő iránytényezője 0.29 és 0.47. Ugyanakkor a légutakban az aerodinamikai 39

veszteségeket minimalizáló elmélet lamináris áramlásra 0.33 [Horsfield 1967, Reis 2004], turbulensre pedig 0.43-at ad [Bejan 2001], amik eléggé közel vannak a fenti értékekhez. Ha tehát elfogadjuk, hogy a légutak geometriája olyan, hogy a szállítási veszteségek minél kisebbek legyenek akkor a turbulencia valahol a harmadik generáció után szűnhet meg. E munkában a felsőlégutakban turbulens, a centrálisban lamináris modellt alkalmaztam. Az első turbulencia modellek [Stapleton 2000, Matida 2004] gyenge, sőt kiábrándító eredményeket produkáltak a légutakban. Az alacsony Reynolds számú modellek megjelenése, úgy tűnik, áttörést hozott. Kleinstreuer nemrég sikerrel alkalmazott egy alacsony Reynolds számú k-ω modellt a nasopharyngeális légúti levegőáramok modellezésére [Kleinstreuer 2003]. Egy több modellt összevető tanulmány [Zhang 2003] szerint is az alacsony Reynolds számokra adaptált k-ω modellek a leginkább alkalmasak a tranziens légúti levegőáramok kiszámítására. Munkám során a k-

ε és a k-ω modellek eredményeit a kísérleti eredményekkel [Cheng 1999] hasonlítottam össze. Végül a a mérésekhez sokkal közelebb álló eredményeket adó k-ω modellt választottam.

log2Dd (cm)

1

0.5

0.25

0.125

0.0625 0

2

4

6

8

10

12

14

16

N

I.18. ábra Légúti csőátmérő (Dd), mint a generációszám (N) függvénye.

40

A k-ω turbulencia modellben a k turbulens kinetikus energia és ω fajlagos disszipációs ráta transzportját a ∂ ∂ ∂k ∂ ( ρk ) + ( ρkui ) = ( Γk ) + Gk − Yk ∂t ∂xi ∂x j ∂x j

(5.20)

∂ ∂ ∂ ∂ω ( ρω ) + ( ρωui ) = ( Γω ) + Gω −Y ω ∂t ∂xi ∂x j ∂x j

(5.21)

és

egyenletek írják le. A fenti egyenletekben ui a sebességkomponenseket, xi pedig a koordinátákat jelöli. Gk és Gω a k és az ω keletkezése, Yk és Yω pedig ezen mennyiségek disszipációja. Továbbá Γk és Γω jelöli a k és ω effektív diffúzivitását, amit a Γk = µ +

µt σk

(5.22)

Γω = µ +

µt σω

(5.23)

és

képletekkel adunk meg. A képletekben σκ és σω k-ra és ω-ra vonatkozó turbulens Prandtl számok. A µt turbulens viszkozitást a

µt = α *

ρk ω

(5.24)

összefüggéssel modellezzük. Az α* együttható nagy Reynolds számú k-ω modellek esetében egynek vehető [Wilcox 1993]. A kis Reynolds számú k-ω modellek α* -ot az ún. korrekciós függvényekkel közelítik [Peng 1997, Bredberg 2002]. A FLUENT-be épített korrekciót az

41

α* =

α ∞* (0.024 + Ret / 6) 1 + Ret / 6

(5.25)

egyenlet jelenti, ahol

Ret =

ρk . µω

(5.26)

A transzportegyenletekben szereplő Gk és Gω tagok modellezéséhez a Gk = µt S 2

(5.27)

és Gω = α

ω k

Gk

(5.28)

képleteket használtam, ahol S-et az Sij feszültségtenzor segítségével az

S = 2Sij Sij

(5.29)

egyenletből kaptam. Az α együtthatót a jelen kis Reynolds számú k-ω modellnél az

α=

α ∞ 0.11 + Ret / 2.95 ( ) α * 1 + Ret / 2.95

(5.30)

összefüggés adja meg. A turbulens disszipációs tagokat (Yk és Yω) a Yk = ρkω

(5.31)

Yω = 0.83ρω 2

(5.32)

és

42

kifejezések írják le. A falközeli réteg modellezésére az ún. „enhanced wall treatment” modellt használtam, amely egy Kader által javasolt eljárással [Kader 1993] sikeresen kombinálja a viszkózus alréteg lineáris és a turbulens régió logaritmikus faltörvényeit [Fluent Manual 2001].

43

6. Részecske trajektória és kiülepedési modell

6.1. Kiülepedési mechanizmusok és matematikai leírásmódjuk Az Euler-Lagrange eljárás lehetővé teszi, hogy a levegő sebességterének és a fő kiülepedési mechanizmusoknak az ismeretében, a belélegzett részecskék pályáját kövessük és kiülepedés esetén annak pontos helyét meghatározzuk. Az általam figyelembe vett depozíciós hatások az impakció, a gravitációs ülepedés és a Brown diffúzió voltak. Mint ismeretes [Ounis 1991, Crowe 1998], normál hőmérsékleti és nyomás viszonyok mellett a Brown mozgás (vagy Wiener folyamat [Wikipedia 2005]) csak mikron alatti részecskék esetében jelentős. Nagyobb részecskéknél az impakció dominál, de a gravitáció is fontos lehet. Mivel a levegő és a légutak hőmérsékletét azonosnak vettem, a fal melletti termoforetikus hatástól [Talbot 1980] e jelen közelítésben eltekintettem. Az egyéb hidrodinamikai erők nagyságrendi okok miatt hanyagolhatók el. A levegő és a részecskék sűrűsége közötti különbség miatt a felhajtó erő, a virtuális tömegerő és a nyomás erő igen kicsik [Ounis 1990]. A Saffmann erő szintén elhanyagolható az alacsony nyíróhatások miatt [Saffmann 1965, Li 1992). Így a részecske mozgásegyenlete (az x komponensekre felírva) du p dt

= β (u − u p ) + g x

(6.1)

lesz, ahol u és up a levegő illetve a részecske sebességének, gx pedig a gravitációs vektornak az x komponense. A jobb oldalon szereplő determinisztikus erők közül a β(uup) az egységnyi tömegre jutó súrlódási erőt jelenti, ahol a mikronnál nagyobb

részecskékre fennáll a

β=

18µ CD Rer ρ p d p 2 24

44

(6.2)

összefüggés. Látható, hogy a súrlódási erő a levegő viszkozitásán (µ) túlmenően a részecske átmérőjétől (dp) és sűrűségétől (ρp), valamint a részecskére jellemző relatív Reynolds számtól is függ, amelyet a

Re r =

ρd p | up − u | µ

(6.3)

képlet ad meg. A CD súrlódási együtthatót a

CD = a1 +

a 2 a3 + Re Re 2

(6.4)

kifejezéssel közelíthetjük, ahol a1, a2 és a3 empirikus konstansokat jelölnek (lásd [Haider 1989, Morsi 1972]). Mikrométernél kisebb részecskékre a

β=

18µ 2 d p ρ p Cc

(6.5)

Stokes képlet érvényes, ahol a nevezőben szereplő Cc mennyiség a Cunningham csúszási faktor. A faktor nagyságát a részecskeátmérő függvényében a

Cc = 1 +

2λ − (1.1d p / 2 λ ) (1.257 + 0.4e ) dp

(6.6)

empirikus képlet adja meg. A modellezés során használt Cc értékeket (ami 69,1 nm-es közepes szabad úthossznak, azaz 106 N/m2 nyomásnak és 298.15 K hőmérsékletnek felel meg, [Lide 1998]) az I.2. táblázatban tüntettem fel

45

I.2. táblázat Cunningham csúszási faktorok különböző részecskeátmérőkre

10-3

dp (µm)

2 x 10-3

229.56 115.06

Cc

5 x 10-3 62.3

10-2

2 x 10-2

5 x 10-2

10-1

2 x 10-1

5 x 10-1

23.86

12.04

5.22

2.99

1.92

1.35

Mikrométernél kisebb átmérőjű részecskéknél a determinisztikus erőkön kívül egy sztochasztikus hatás, a Brown mozgás is érvényesül. A termikus energiával rendelkező levegőmolekulák rendezetlen mozgásuk során ütköznek a részecskével (I.19. ábra).

I.19. ábra A Brown mozgás szemléltetése.

Az ütközéssel együtt járó impulzuscsere kis részecskék esetében a részecske elmozdulását eredményezi, amit a részecske trajektóriák modellezésénél figyelembe kell venni.

A

levegőben

normál

nyomáson

és

hőmérsékleten

másodpercenként

nagyságrendileg 1010 ütközés történik, így nincs esély a részecske útvonalának részletes követésére, ami ráadásul soktest-problémára vezetődne vissza. Ilyen körülmények között csak a sztochasztikus megközelítés lehet célravezető. Ennek következtében a részecske mozgásegyenletében a determinisztikus tagok mellett megjelenik egy a részecske sebességtől független nulla átlagú A(t) fluktuációs tag is [Gupta 1985]. Ennek megfelelően a mozgásegyenlet a

46

du p 1 = (u − u p ) + A(t ) dt τ p

(6.7)

alakot ölti. A fenti Langevin típusú egyenletben τp a részecske relaxációs ideje, a β reciproka. A mozgásegyenlet megoldása egy W(up, t, up0) valószínűség eloszlás lesz, ami megadja az up érték előfordulásának valószínűségét a t időpillanatban, ha az a t = 0 pillanatban up0 volt. Az egyenlet megoldásakor feltételezzük, hogy a ∆t időlépés sokkal nagyobb a két egymás utáni ütközés közötti időtartamnál. Ugyanakkor azt is feltételezzük, hogy a W valószínűségi eloszlás a t → 0 határesetben a Dirac eloszláshoz, t → ∞ esetben pedig a Maxwell eloszláshoz tart [Chandrasekhar 1943]. A szabad

részecskére vonatkozó Langevin egyenlet formális megoldása t

u p −u p 0 e − β t = e − β t ∫ e βξ A(ξ )d ξ .

(6.8)

0

t

Ennek következtében az u p −u p 0 e− β t és a e − β t ∫ e βξ A(ξ )d ξ kifejezéseknek statisztikai 0

tulajdonságai azonosak kell, hogy legyenek. Ez azt jelenti, hogy a

lim e t →∞ 

− βt

mu 2p

m − 2 kBT = ξ ξ e ( ) e A d ∫0  2k BπT t

βξ

(6.9)

feltételnek is teljesülni kell. Itt T az abszolút hőmérsékletet, kB pedig a Boltzmann konstanst jelöli. Ha az időlépés annyira kicsi, hogy azalatt a determinisztikus hatásokból származó sebesség változás elhanyagolható akkor az e βξ tényező kihozható az integrál elé és a B ( ∆t ) = ∫

t + ∆t

t

A(ξ )dξ

(6.10)

jelölést bevezetve kimutatható, hogy a (6.9) feltétel akkor teljesül ha a B gyorsulás valószínűség-eloszlását a

47

1

w [ B( ∆t )] =

(4π q∆t )

1 2

e

−

B ( ∆t ) 4 q∆t

2

(6.11)

képlet adja meg [Chandrasekhar 1943]. Itt q a

q=

k BT mτ p

(6.12)

képlettel megadott mennyiség, m pedig a részecske tömege. A FLUENT CFD kódba épített Brown modell a B értékeit a fenti eloszlásból sorsolja, azt minden időlépésben a mozgásegyenlethez csatolja és kiintegrálja a mozgásegyenletet. Ugyanehhez az eljáráshoz juthatunk, ha a Brown mozgást Gauss-féle fehér zajként modellezzük melynek spektrális intenzitása 2kTB/τpπ m [Liu 2005]. Mint a fenti levezetésből is kiderült, a gyorsulások eloszlása olyan időlépésekre igaz, amelyek jóval nagyobbak a két ütközés között eltelt időnél, de nem nagyobbak a részecske relaxációs idejénél. Mint említettem, normál körülmények között a levegő molekulák ütközései között 10-10 s nagyságrendű idő telik el. Mivel a részecskeátmérője és ezáltal ütközési hatáskeresztmetszete a levegő molekuláiénál nagyobb, a részecske által elszenvedett ütközések közötti idő még ennél is rövidebb. A relaxációs időnek néhány jellegzetes részecskeméretre megadott értékét a I.3. táblázat tartalmazza I.3. táblázat A részecske relaxációs ideje, mint a részecskeméret függvénye dp (µm) τp (s)

1

0,1

0,01

3x10-6

8,97x10-8

7,16x10-9

0,001 6,89x10-10

Mint látható, a jelenlegi számítási kapacitások mellett még realisztikus gépidőt eredményező 10-6 másodperces időlépés általában jóval a fenti értékek fölött van, ezért a modell nem lehet alkalmas belélegzett részecskék Brown mozgásának szimulálására. A problémát egy új módszer kidolgozásával oldottam meg. A modell felépítése során olyan összefüggéseket használtam, amelyek ∆t >> τp időlépésre érvényesek és olyan

48

eljárást eredményeznek, amely beépíthető a FLUENT trajektória kódjába. Kézenfekvő lenne Brown erő helyett elmozdulást számítani és ezt minden időlépésben hozzáadni a determinisztikus erők által okozott elmozduláshoz. Ez gyakorlatilag azt jelentené, hogy a mozgásegyenletnek megfelelő diffúziós (vagy Fokker-Planck, vagy külső erőmentes Smoluchowski) ∂W ( x, t ) ∂ 2W ( x, t ) =D ∂t ∂x 2

(6.13)

egyenlet W (∆xB , t ) =

1 4π Dt

e

−|∆xB |2 4 Dt

(6.14)

megoldásából [Uhlenbeck 1930] sorsolunk elmozdulás értékeket, ahol a D diffúziós együtthatót a

D=

k BTCC 3πµ d p

(6.15)

fluktuáció - disszipáció egyenlet (vagy Einstein képlet, lásd [Einstein 1905]) adja meg. Ez az eljárás azért nem alkalmazható, mert nincs mód a FLUENT kódba elmozdulást bevinni. Ezért egyedüli megoldásnak az tűnt, hogy kidolgozzak egy olyan Brown erőt szimuláló eljárást, ami nagy időlépésekre is érvényes képleteken alapszik. A ∆t >> τp esetben a Brown erőre az

u p = ς (0,1) β

k BT m

(6.16)

kifejezés adódik, amit úgy kapunk, hogy a (6.9) egyenlet jobb oldalán lévő Maxwell eloszlásból sorsolunk sebességet, majd kiintegráljuk az egyenletet. A sebesség sorsolást úgy oldottam meg, hogy a nulla átlagú és egységnyi szórású Gauss eloszlásból ς (0,1) véletlen számokat sorsoltam, majd ezeket megszoroztam a Maxwell eloszlás szórásával.

49

Egy a fenti eljárással előállított Brown erő minta az I.20. ábrán látható.

A(t)

t

I.20. ábra Szimulált Brown erő minta. Az így előállított Brown erőt minden időlépésnél a (6.7) mozgásegyenletbe helyettesítve egy közönséges differenciál egyenletet kapunk, amit tovább a FLUENT trajektória kódja old meg. Érdemes megjegyezni, hogy a Brown sebesség eloszlását tetszőleges időre a 1

 2 − m W (u p , t , u p 0 ) =  e −2 β t   2π k BT (1 − e ) 

m ( u p − u p 0 e − β t )2 2 k BT (1− e−2 β t )

(6.17)

képlet adja meg [Balescu 1975, Kac 1987]. Belátható, hogy az eloszlásból sorsolt sebesség ∆t<<τp és ∆t>>τp határesetekben visszaadja, a FLUENT illetve az általam használt értékeket. A Brown modellt C++ nyelven programoztam be [Press 2002], mielőtt a FLUENT trajektória kódjához csatoltam. A véletlen számokat Visual C++ 6.0-val generáltam miután a generátort teszteltem. Mivel a FLUENT-nek van saját C++ fordítója, a csatolás egy speciális könyvtárstruktúra létrehozásával viszonylag könnyen megoldható [FLUENT Manual 2001]. A modell validálását a 7.2 fejezet tartalmazza.

50

6.2. A mozgásegyenlet numerikus megoldása A részecskepályát a dx p dt

= up

(6.18)

egyenlet megoldásából kaphatjuk meg numerikus integrálással, ahol up a du p dt

=

1 (u − u p ) lp

(6.19)

linearizált egyenlet megoldása (lp az arányossági tényező reciproka). A fenti egyenletet a

up

n +1

− up

∆t

n

=

1 * n +1 (u − u p ) lp

(6.20)

iterációs formula segítségével oldhatjuk meg, ahol

u* =

1 n +1 (u + u n ) 2

(6.21)

és n

u n +1 = u n + ∆tu p ∇u n .

51

(6.22)

6.3. Részecskesorsolás és peremfeltétel A fenti iterációs technika akkor alkalmazható, ha a kezdeti és peremfeltételek adottak. A modellezett légúti szegmensekbe jutó részecskék kezdeti helyét és sebességét egy Monte-Carlo eljárással sorsoltam. Ez azt jelenti, hogy a részecskék koordinátáit és sebességét a bemeneti sebesség profilnak megfelelően állítottam elő. Ennek eredményeként, a bejáraton ott lesznek sűrűbben a részecskék, ahol a levegő sebessége nagyobb. A részecskék sebességét minden esetben az adott pontra jellemző levegő sebességgel tettem egyenlővé. A modellezés során gömb alakú, különböző aerodinamikai átmérőjű és 1000 kg/m3 sűrűségű részecskéket feltételeztem. Statisztikai meggondolásból a részecskék száma sosem volt kisebb, mint 10 000, de a legtöbb esetben 50 000, 100 000 vagy éppen 107 részecske pályáját követtem. Peremfeltételként azt írtam elő, hogy a fallal ütköző részecske kitapadjon (teljesen rugalmatlan ütközés) és ott is maradjon. Eszerint a trajektória fallal való metszéspontja egyben a pálya végpontja is, aminek koordinátáit rögzíteni lehet. Figyelembe véve, hogy a légutak falát egy nagy viszkozitású vékony nyákréteg borítja [Asgharian 2001, Rubin 2002] és az adhéziós erők annyira nagyok, hogy nincs reszuszpenzió [Auvinen 2004, Grzybowski 2004], a feltevés realisztikusnak mondható.

52

7. A numerikus modell validálása

A későbbi modellszámítások eredményeinek hitelessége szempontjából az 5. és 6. fejezetben

bemutatott

számítási

modell

alkalmasságának

széleskörű

vizsgálata

kulcsfontosságú. Ezért munkám során különös hangsúlyt fektettem a numerikus modellem által számított eredményeknek a szakirodalmi mérési adatokkal és szimulációs eredményekkel történő összehasonlítására. A validáció során légúti elágazásokban mások által mért, vagy számított sebességprofilokat és kísérleti vagy szimulációs kiülepedési eloszlásokat, illetve a kiülepedett és a bemenő részecskeszám arányaként értelmezett kiülepedési hatásfokokat hasonlítottam össze az általam számítottakkal.

7.1. A Navier-Stokes egyenlet megoldásának validálása

Légúti elágazások áramlási terének jellemzéséhez szükség van a levegő sebesség és nyomás értékeire az elágazás belső pontjaiban. Ezen értékek mérése élő ember esetében rendkívül nehéz feladat. Modell elágazásokban is igen nehéz a nyomáseloszlás mérése, mivel a nyomásesés egy elágazásra nagyságrendileg nem nagyobb, mint a csőátmérő mentén fellépő nyomáskülönbségek [Pedley 1977]. Ellenben a sebességméréshez számos módszer áll rendelkezésre. Ilyenek a részecskekép alapú sebességmérés (PIV) a lézer Doppler anemometria (LDA), de újabban a MicroPIV módszerrel már bronchioláris méreteknek megfelelő elágazásokban is hatékonyan mérnek sebességprofilokat [Lee 2005]. Az áramlási modell validációja érdekében kézenfekvő tehát, hogy a mért sebességértékeket az azonos geometriájú és térfogatáramú esetekre számítottakkal összehasonlítsuk. Elsőként a 3.3. fejezetben leírt, a felnőtt tüdő Weibel modellje 3-4. generációjának megfelelő elágazás áramlásterét vizsgáltam. E szimmetrikus légúti elágazásban 10 illetve 60 l/min tracheális térfogatáram mellett az általam vizsgált numerikus modell az I.21. ábrán látható sebességteret eredményezte.

53

60 l/perc

10 l/perc

Belégzés

0.5 m/s

3.6 m/s

10 l/perc

60 l/perc

Kilégzés

I.21. ábra Sebesség szintvonalak a Weibel modell 3-4. generációs légúti elágazásának fősíkjában és radiális sebességvektorok a leányág (fent) vagy anyaág (lent) felénél található kör keresztmetszetekben. A tracheális térfogatáram 10 l/perc (bal oldal), illetve 60 l/perc (jobb oldal). A felső panel belégzést, az alsó panel kilégzést szemléltet.

Úgy a fenti (belégzés) mint a lenti (kilégzés) ábrarészen szintvonalas ábrázolás látható az elágazás szimmetria síkjában. Alatta pedig radiális sebességvektor mezőt ábrázoltam az egyik leányág (belégzés) vagy az anyaág (kilégzés) közepénél egy a 54

hengeres rész tengelyére merőleges síkban. Az I.21. ábra alapján a numerikus modell visszaadja a Pedley által kísérletileg tapasztalt valamennyi jelenséget [Pedley 1971]. Így belégzéskor, különösen a nagyobb térfogatáram esetében, megjelenik az alacsony sebességértékekkel jellemzett szeparációs zóna a centrális rész után a leányágak külső szélénél. Ugyanakkor az anyaágban még parabolikus sebességprofilok a leányágakban aszimmetrikussá válnak és a sebesség maximumok a leányágak belső pereme felé tolódnak el. Az eltolódás fokozódik a bemeneti sebesség növekedésével. Ez összhangban áll többek között Schroter, Zhao és Caro méréseivel [Schroter 1969, Zhao 1994, Caro 2001], valamint Liu, Balásházy és Comer [Liu 2002, Balásházy 1993a,b, Comer 2001a] szimulációs eredményeivel is. Továbbá, a radiális belégzési sebesség komponenseket mutató ábrákon fellelhető az Isabey által kimutatott Dean típusú szimmetrikus örvénypáros [Isabey 1982]. Ugyanez igaz a kilégzéskor megjelenő négyes örvényre. Zhao a fősíkra merőleges síkban M alakú profilt mért a leányág felénél. Amint azt az I.22. ábra szemlélteti, az I.21. ábra szerinti AB szakaszra merőleges egyenes mentén a számítások alapján is megjelenik az M alakú sebességprofil. Az AB szakasz mentén pedig a már említett ferde profil észlelhető.

I.22. ábra Axiális sebességprofilok az I.21. ábra szerinti elágazásában az egyik leányág közepén, az AB szakasz mentén (bal oldali ábra) és rá merőlegesen (jobb oldali ábra). A vektorok mindkét esetben a maximális axiális sebességre lettek normálva.

Miután minőségileg az áramlási tér kitűnő hasonlóságot mutat a mások által mért és számított terekkel, célszerű megvizsgálni a mennyiségi egyezést is. Ennek érdekében az általam számított sebességprofilokat összehasonlítottam a Raick által kísérletileg meghatározottakkal [Raick 2003]. Raick a méréseket egy háromgenerációs légúti geometriában végezte PIV módszerrel. A kísérleti eredmények és a számított profilok

55

összehasonlítását az I.23. ábra mutatja a geometria fősíkjában, a bal oldali dupla elágazáson jelölt helyeken. 1.0 0.8 0.6

Dimenziótlan sebesség, v/v(-) max (-) Dimensionless velocity

0.4

◊ PIV mérés [Raick 2003] PIV measurement + Jelen szimuláció simulation

0.2 0.0 0.0

0.2

0.4

0.6

0.8

1.0

1.0 0.8 0.6 0.4

◊ PIVPIVmérés [Raick 2003] measurement simulation + Jelen szimuláció

0.2 0.0 0.0

0.2

0.4

0.6

0.8

1.0

1.0 0.8 0.6 0.4

◊ PIV mérés [Raick 2003] PIV measurements + Jelensimulation szimuláció

0.2 0.0 0.0

0.2

0.4

0.6

0.8

1.0

Dimensionless width Dimenziótlan távolság, d/D(-) d (-)

I.23. ábra Mért és számított sebességprofilok egy háromgenerációs légúti szegmensben. Az ábrán relatív sebességértékek szerepelnek. A bemenetre számolt Reynolds szám 700, a bemeneti sebességprofil parabolikus; v – sebesség; vmax – maximális sebesség; d – távolság; Dd – csőátmérő.

Az I.23. ábrán relatív sebességértékeket láthatunk redukált távolság függvényében a második és harmadik generációs ágak közepén, a megjelölt egyenesek mentén. Látható, hogy a modellezett profilok jó egyezést mutatnak a mértekkel. Az eltérések nem nagyobbak, mint a PIV és LDA módszerekkel mért profilok közötti különbségek [Ramuzat 1998]. A validálás során a modellel számított eredmények rácsfüggetlenségét is vizsgáltam. A korai hálóknál ezt nem mindig sikerült elérni, de az ún. „size function” technika alkalmazásával ez a probléma megoldódott. 56

7.2. A trajektória modell validálása

Miután az áramlási modell jó hasonlóságot mutatott a mért és számított adatokkal, megvizsgáltam, hogy az áramlási térben számolt részecskepályák és a trajektóriák fallal való metszetéből származó kiülepedési helyek mennyire egyeznek az irodalmi értékekkel. Viszonyítási alapul ismét csak az elérhető kísérleti eredményeket és modellszámításokat használtam. Kísérletileg Johnston és később Kim is inhomogén kiülepedést észlelt a Weibel modell szerinti 3-4. generációs légúti elágazásban [Johnston 1977, Kim 1989]. Az elágazások csúcsában és a leányágak belső falán a kiülepedés intenzívebb volt az átlagosnál. Ugyanebben a geometriában 60 l/perc – es térfogatáramra kiszámoltam az 1, 10 és 100 nanométeres valamint 1, 5 és 10 mikrométer aerodinamikai átmérőjű részecskék kiülepedési helyét. A háromdimenziós kiülepedési képek megjelenítéséhez (lásd I.24. ábra) egy C++ nyelvű programot írtam, amit a FLUENT-hez csatoltam.

dp = 1 nm

dp = 1 µm

dp = 10 nm

dp = 5 µm

dp = 100 nm

dp = 10 µm

I.24. ábra Inhalációs kiülepedési eloszlások. Bemenő részecskék száma 50 000, a belépő profil parabolikus, a térfogatáram pedig 60 l/perc; dp – részecskeátmérő.

57

Az

I.24.

ábra

tanúsága

szerint,

a

kiülepedett

részecskék

eloszlása

minden

részecskeméretre inhomogén, de legfőképpen a mikro-részecskékre. A kísérleti megfigyelésekkel összhangban a kiülepedés az elágazások csúcsában a legintenzívebb. A számszerű egyezés vizsgálatához a kiülepedési hatásfokot ábrázoltam az anyaágra vonatkoztatott Stokes szám függvényében (I.25. ábra ). A Stokes szám az inerciális és a viszkózus erők viszonyszáma, az impakció mértéke. Értékének kiszámításához a

St =

ρ p d p2 u 18µDd

(7.1)

képletet használtam, ahol ρp és dp a részecske sűrűsége, illetve átmérője, µ és u a levegő viszkozitása, illetve átlagsebessége a bejáraton, Dd pedig az anyaág átmérője.

Deposition hatásfok efficiency (%) Kiülepedési (%)

60

Kísérlet Kisérlet [Kim 1989] Kim & Iglesias Kísérlet Johnston et[Johnston al. Kisérlet 1977] Present model Jelen modell Jelen modell

50 40 30 20 10 0 0.0

0.1

0.2

0.3

0.4

0.5

0.6

Stokes Stokesnumber szám(-) (-)

I.25. ábra A mért és modellezett kiülepedési hatásfokok összehasonlítása.

E mennyiségi összehasonlításból kiderül, hogy a kísérletileg meghatározott és a számított hatásfokok között jó az egyezés. Ennek ellenére az összehasonlítás a nagyon alacsony Stokes számokra (0.1 alatt) nem ad érdemi információt, mert a mérések ezen érték felett történtek. Viszont e tartományban mért kiülepedést Kim az áramlások validálásánál már 58

bemutatott dupla elágazásban [Kim 1999]. Az általa végzett kísérletek és saját modellszámításaim eredményeinek összehasonlítását az I.26. ábra szemlélteti.

Deposition efficiency (%) (%) Kiülepedési hatásfok

30 25

Experiment (Kim1999] & Fisher 1999) Kísérlet Kisérlet [Kim Simulation (Comer2001] et al 2001) Elmélet [Comer Present model Jelen modell

20 15 10 5 0

0.04

0.06

0.08

0.10

0.12

Stokes number Stokes szám(-)(-)

I.26. ábra Szimulációs kiülepedési hatásfokok összehasonlítása a mások által mért és számított értékekkel.

Az I.26. ábra szerint, a depozíciós számítási modellem kis Stokes számokra is jól adja vissza a mért kiülepedési hatásfok értékeket. A Comer által ugyanabban a légúti geometriában numerikusan számított hatásfokokkal szintén jó az egyezés [Comer 2001b]. A részecske depozíció inhomogenitását realisztikus légúti elágazásokban is ellenőriztem. Viszonyításul a patkánytüdő egy légúti elágazásának központi részéhez közeli metszeteiről röntgensugaras fluoreszcencia (XRF) technikával készült felvételeket használtam. A metszetek eozin aeroszol inhaláltatását követően készültek. A méréseket Dr. Török Szabina, Dr. Osán János és Dr. Alföldy Bálint végezték. A metszeteket szegmentáltam (a zöld színű szegmentációs kontúrokat lásd az I.27. ábra bal oldalán) és az így nyert kontúrvonalakra a 3D-DOCTOR szoftverrel háromdimenziós numerikus felületet fektettem. Ez a felület egy légúti elágazás központi részét jelenti, amihez anyaágat

és

leányágakat

varrtam,

hogy

részecskekiülepedést modellezhessek benne.

59

numerikus

levegőáramlást

és

-500

-400

-300

-200

-100

0

100

200

300

400

500

600

700

800

900

1000

2000 1800 1600 1400 1200 1000 800 600 400 200 0 -200 -400 -600 -800 -1000 2000

1800

1600

1400

1200

1000

800

600

400

200

0

-200

-400

-600

-800

-1000

-400 -200

0

200

400

600

800 1000

2000

1800

1600

1400

1200

1000

800

600

400

200

0

-200

-400

-600

-800

-1000

-500-400-300-200-100 0 1002003004005006007008009001000

I.27. ábra Szegmentált tüdőmetszetek (bal oldal), kísérletileg nyert (középső) és számított (jobb oldal) kiülepedési képek. 60

Az így készült geometria az I.28. ábrán látható. Az I.27. ábra középső oszlopában a kiülepedett részecskék láthatók a bal oldali oszlopnak megfelelő metszeteken. A középső panel XRF-es felvételeinek bal oldalán látható fehér foltok vérereket jelölnek. A kiülepedési foltokat a két leányág között körök jelölik. A jobb oldali oszlopban a számított kiülepedési képek láthatók. Az ábra középső és jobb oszlopait összehasonlítva látható, hogy a karina környékén a kiülepedési eloszlások hasonlóak. Látható, hogy a számítási eredmények jól modellezik a kiülepedést az első három metszeten. Míg a mérések szerint a negyedik-ötödik metszeten már szinte nincs kiülepedés, a szimulált metszeteken még látható egy enyhe depozíció. Ennek az lehet az oka, hogy adat hiányában a levegőáramot csak megbecsülni tudtam, valamint a bemeneti sebességprofilt ugyancsak konkrét mérés hiányában parabolikusnak vettem. Ennek ellenére az összehasonlítás minőségileg jó egyezést mutat, a karina régió fokozott kiülepedését a számítási modell jól írja le.

I.28. ábra A rekonstruált központi rész (bal oldal) és a teljes légúti elágazás (jobb oldal). A 6.1. fejezetben leírt eljárás validációjához 1000 darab 10 nanométeres részecske pályáját követtem, amelyek a nulla időpillanatban a termikus egyensúlyban lévő 300K hőmérsékletű levegő ugyanazon pontjában helyezkedtek el. E kezdeti ponttól vett 61

távolságok x komponenseinek négyzetátlagát az idő függvényében az I.29. ábrán tüntettem fel. Az egzakt megoldás az x=(2Dt)1/2 képletnek felel meg. A elmozdulásokat a 6.1. fejezetben leírt módszerrel sorsoltam.

-2 cm) cm) Sqrtsqrt(<x (<x22>) >)( x( 10 x -210

10

8

egzakt megoldás szimuláció

6

4

2

0

0

20

40

60

80

100

t (s)

Idő (s)

I.29. ábra Elméleti és szimulációs elmozdulás négyzetátlagok összehasonlítása. A szimulált részecskék száma 1000, a részecskék átmérője 0.01 µm, a diffúziós állandó 5x10-4 cm2/s.

Mint látható, a szimulált értékek jól egyeznek az egzakt megoldásnak megfelelő értékekkel. A fent leírt kísérleti és modellszámítások közötti egyezések alapján joggal remélhető, hogy a használt numerikus modell kellő pontossággal leírja a belélegzett részecskék pályáját és kiülepedését a légutakban.

62

II. Rész: A numerikus modell alkalmazása 8. Légúti levegőáramok számítása A rövid felezési idejű radon leányelemek által a centrális légutak érzékeny hámsejtjeire gyakorolt sugárterhelés pontos kiszámításához elengedhetetlen ezek lokális légúti kiülepedésének ismerete. Ez gyakorlatilag minden egyes radio-aeroszol részecske követését feltételezi a belépéstől a kiülepedésig, vagy amíg a részeske elhagyja a tanulmányozott légúti geometriát. A felsőlégutak szűrőként viselkednek és csak a belélegzett részecskék egy bizonyos hányadát engedik tovább. Ezért, a radionuklidok hörgők falára történő kiülepedésének és aktivitás-eloszlásának modellezéséhez szükség van a felsőlégúti kiülepedés ismeretére is [Farkas 2003d, 2004d]. Mivel a radon bomlástermékek az őket szállító levegővel együtt jutnak a felsőlégutakba és a levegő áramlása további útjukat is befolyásolja, ezért először a levegő áramlásterét tanulmányoztam. A felsőlégúti levegőáramok modellezéséhez az 5.2.7. fejezetben bemutatott k-ω turbulencia modellt használtam. A II.1. ábra a szájon át belélegzett levegő sebességvektorait szemlélteti a 3.3. fejezetben leírt légúti geometria fősíkjában. Az ábra tanúsága szerint a kialakuló sebességmező meglehetősen bonyolult. A falak mentén több helyen visszaáramlást tapasztalhatunk. Ilyen zónákat mutatnak a kinagyított részek is. A bemenet közelében a közel félkör alakú részben a Tsuda által is leírt [Tsuda 2002] örvény található. A második kinagyított részben látható sebességmező a szűkület és a könyök hatásaként jelenik meg. Szembetűnő továbbá, hogy a kezdeti egyenletes sebességeloszlás gyorsan megváltozik és lokálisan lassú és gyors áramlási zónák alakulnak ki. A nagy sebességekkel jellemzett helyeken előreláthatóan igen jelentős lesz a nagy részecskék impakciós kiülepedése. Ilyen térrész például az ábrán B-vel jelölt környezet, ahol a szűkület okozta gyorsuláson kívül még a kanyaros geometria is elősegíti az impakciós kiülepedést. Más helyeken (pl. A környezet) a nagyon kis részecskék diffúziós kiülepedése lesz számottevő.

63

A

m/s

B

II.1. ábra Felsőlégúti levegőáram sebességvektorai a geometria fősíkjában. A kinagyított részek örvényt és visszaáramlást szemléltetnek. A belépő sebességprofil egyenletes, a térfogatáram 18 l/perc. Légzés során a levegő be- és kiáramlása a légutak különböző pontjai között fellépő nyomáskülönbségek hatására történik. A nyomásesések a térfogati levegőárammal és a légúti

ellenállással

vannak

összefüggésben.

A

tanulmányozott

száj-garat-gége

geometriában különböző térfogatáram értékekre végeztem nyomásesés számításokat. Az eredményeket a II.2. ábra foglalja össze. Hyatt és Wilcox [Hyatt és Wilcox 1961] kísérletileg kimutatták, hogy a térfogatáram növekedésével a felső légutak ellenálása csökken. Numerikus eredményeim összhangban vannak a fenti megállapítással. A II.2 ábráról az is látszik, hogy az extrathorakális nyomásesés még a legnagyobb térfogatáramokra is kicsi a normál légköri nyomáshoz képest. A négy nagyságrendnyi eltérésből következik, hogy joggal vehetjük a levegőt összenyomhatatlannak a légúti levegőáram és részecsketranszport számításaink során. Ha figyelembe vesszük, hogy Wilcox és Hyatt szerint a felső légúti nyomásesés a teljes légúti nyomásesésnek durván a

64

negyven százalékát teszi ki, akkor a fenti nagyságrendi megállapítások nemcsak a felső

Térfogatáram (l/perc)

légutakra, de a teljes légúti rendszerre is igazak.

Nyomásesés (Pa) II.2. ábra A térfogatáram és a nyomásesés közötti összefüggés a tanulmányozott szájgarat-gége modellben. Addig míg a száj-garat-gége térfogati levegőárama kinondottan a légzési paraméterek (légzési frekvencia, légzési térfogat) függvénye, a centrális légutak esetében a levegőáram

az

előző

elágazásokon

belüli

hozam

eloszlásától,

a

bemeneti

sebességprofiltól, a leányágak átmérőitől és az elágazási szögektől is függ. A II.3. ábra térfogatáram eloszlást szemléltet a légcső utáni néhány elágazásban. Látható, hogy a szimmetrikus bemeneti sebességprofil ellenére, az egyes csövekben áramló levegő mennyisége csőről-csőre változik még az azonos generációszámú csöveknél is. Várható tehát, hogy az egyes légúti járatokba jutó radioaktív részecskék száma (és ezáltal az ott kiülepedett részecskék száma is) nagyon különböző lesz az egyes hörgőkben. Ahhoz viszont, hogy ezt pontosan kiszámíthassuk ismernünk kell a tracheobronchiális áramlási teret [Farkas 2001a,c,d].

65

Térfogatáram (l/perc)

18

18

16 14

1

12

10.57 10,57

10

9

8,38

8

7.4 7,43

8

11

6

10

4 2,19

2

0,8

1

2

3

4

5

6

2 7

6

3

5 4

1,39

7

0,21 0,59

8

9

0,51

10

0,88

11

A légúti szakasz száma II.3. ábra Levegőáram eloszlás a tanulmányozott aszimmetrikus légutakban. Az abszcisszán feltűntetett számok az ábra jobb oldalán látható számozásnak felelnek meg. A tracheális térfogatáram 18 l/perc. A bemeneti sebességprofil parabolikus.

A centrális légutakban lamináris áramlást tételeztem fel. A II.4. ábra bal oldalán a 4. fejezetben értelmezett sebesség szintvonalakat láthatjuk a fenti légúti geometria elágazásainak fősíkjaiban. A korábbi, egyszerű szimmetrikus elágazásokban végzett számítások eredményeivel ellentétben [Heistracher 1996] és más, aszimmetrikus geometriákat figyelembe vevő szakirodalmi publikációkkal [Horsfield 1971, Philips 1997, Mauroy 2003] összhangban a belélegzett levegő sebességtere messzemenően nem szimmetrikus. Mivel olyan komplex geometriákban, mint amilyen a II.3. ábrán is látható, még a szimmetrikus elágazás (gen 4-5.1) sebességtere is erősen aszimmetrikus, egyértelmű, hogy a radionuklidok pontos trajektóriáit csak realisztikus alakú geometriákban tudjuk valósághűen számolni. Érdekes jelenség továbbá, hogy az egyszerű 66

geometriák esetében az elágazás után a külső fal mellett kialakuló lassú határréteg [Pedley 1971] jelen esetben csak az 1-2, részben a 3-4 és 4-5.1 elágazásoknál figyelhető meg, de szinte teljesen hiányzik a 2-3 és 4-5.2 elágazásoknál. Az aszimmetrikus sebességprofilok (a II.4. ábra jobb oldala) arra engednek következtetni, hogy a részecske kiülepedés nem csak bronchusról-bronchusra lesz eltérő, de egy adott bronchus falán is inhomogén lesz. m/s

A3

gen 3-4

A5

B3

A1

B1

A2

B2

A3

B3

A4

B4

A5

B5

B5

gen 4-5.2

A1

A4

B4

gen 4-5.1

A2

B1

gen 1-2

B2

gen 2-3

II.4. ábra A belélegzett levegő sebességterének szintvonalas ábrázolása öt elágazás fősíkjában (bal oldal) és az öt légúti elágazás bemeneti sebességprofilja (jobb oldal). A légcsőre vonatkoztatott térfogati levegőáram 18 l/perc.

A centrális légúti elágazások morfológiájának minél realisztikusabb modellezését célzó, a 3.3. fejezetben leírt karina lekerekítő eljárás is újabb jelenség megfigyelésére

67

adott lehetőséget. A II.5. ábrán látható radiális sebességvektorok az elágazás csúcsának közelében felvett síkmetszetben a faltól kifelé mutatnak. Ugyanakkor ismert, hogy a fal közeli másodlagos áramlások nagy mértékben befolyásolják a részecske kiülepedést [Broday 2004]. Mindezt figyelembe véve valószínűsíthető, hogy az eddigi idealizált, éles karinájú geometriák felülbecsülték az elágazások csúcsa körüli lokális aeroszol kiülepedést.

II.5. ábra Radiális sebességvektorok egy légúti elágazás csúcsának közelében. A keretben a karina csúcsában megfigyelhető visszáram kinagyítva látható. A légcsőre vonatkoztatott térfogati levegőáram 18 l/perc, a légcső bemeneti sebességprofilja parabolikus.

A sebességvektorok iránya és nagysága időben is változik. Igen érdekes momentum ebből a szempontból a be és kilégzés határa, amikor a levegő egy pillanatra úgymond megáll (II.6. ábra). Az ábrán az ún. „start-up” (kezdeti) innerciális effektusok elkerülése céljából a második légzési ciklus eredményei vannak feltüntetve. Egyes szerzők (pl. [Moskal 2002]) kimutatták, hogy a részecskekiülepedés ilyenkor az átlagosnál intenzívebb. Ezért indokolt megvizsgálni, hogy az aeroszol depozíció szempontjából 68

mennyire fontos az időfüggő modellezés. Ugyan Gurman [Gurman 1994] és Kim és Garcia [Kim 1991] is bebizonyították, hogy időfüggő ciklikus áramlásnál az aeroszolok kiülepedési hatásfoka nagyobb, mint stacionárius esetben, a stacionárius modellezés is pontos kiülepedési hatásfokot adhat, ha megfelelő Reynolds szám korrekciót alkalmazunk [Zhang 2002]. Szinuszos légzési mintát feltételezve [Tobin 1983, Paek 1992, Hausermann 2000], ez az átlagos Reynolds szám 30%-os növelését jelenti. Ez a 30 % nem nagyobb az egyének vagy nemek közötti térfogatáram (és Reynolds szám) eltérésnél [ICRP66 1994]. A stacionárius modellezés mellett szól az is, hogy sajnos az időfüggő számítások nagyon gépidő igényesek. Egy-egy futás a számítások pontosságához szükséges rácsfelbontás mellett akár hónapokig eltarthat. Ezért áramlási számításaim túlnyomórészt stacionáriusak és átlagos felnőttre jellemző légúti levegőáramlást adnak meg. (m/s) 10.2

t=0.5 s

t=2.5

t=1 s

t=3 s

t=1.5

t=2 s

t=3.5

t=4 s

0

II.6. ábra A levegő sebességtere egy aszimmetrikus légúti elágazásban a légzési ciklus különböző időpontjaiban. A bemeneti sebességprofil egy olyan paraboloid, amelynek minden pontjában a sebesség az idővel szinuszosan változik. A térfogatáram 7.5 l/perc, a légzési frekvencia 15 perc-1.

69

9. Radionuklidok légúti kiülepedésének leírása

Miután az előző fejezetben rámutattam, hogy a belélegzett radioaktív részecskéket szállító levegő sebességtere erősen függ a légzési paraméterektől és a légúti geometriától, jelen fejezetben a részecskekiülepedést tanulmányozom. Először is megvizsgálom, hogy a felsőlégutak milyen mértékben engedik tovább a belélegzett részecskéket. Az extrathoratikus légutak által ki nem szűrt részecskék lokális kiülepedésének az említett légzési módon és légúti morfológián kívül a részecskék paramétereitől való függését is tanulmányozom. Továbbá, mivel a radonszármazékok egészségre gyakorolt hatása elsősorban a lokális terheléstől függ [Balásházy 2003a], mennyiségileg is jellemzem a helyi részecske kiülepedést.

9.1. Extrathoratikus aeroszol kiülepedés Pontos centrális légúti kiülepedés és dozimetria számítások csak a légcsőbe belépő részecskeszám és eloszlás ismeretében végezhetők, melynek egyik feltétele az aeroszolok felsőlégúti kiülepedésének ismerete. Ennek érdekében megvizsgáltam, hogy hogyan tapadnak ki a részecskék az I.9. ábrán is látható oropharingeális légúti geometriában. A gáz halmazállapotú radon bomlása során keletkező izotópok nagy része általában rátapad a levegőben lévő aeroszol részecskékre. A levegő porszemcséi adott helyen és időben egy nanométertől néhányszor tíz mikronig terjedő részecskeméretet magába foglaló eloszlást mutatnak [Reineking 1994]. Mivel a kiülepedési mechanizmusok részecskeméret függőek, várható, hogy a depozíciós hatásfok is függvénye a részecskeméretnek. Ezenkívül a részecskéknek a légutak falával való ütközését az áramlási sebesség, vagyis a térfogatáram nagysága is befolyásolja. Ezért a felsőlégúti depozíciót e két paraméter függvényében vizsgáltam. A számítások során aerodinamikai átmérővel jellemzett gömb alakú részecskéket feltételeztem. A II.7. ábra szájlégzés esetén a mikrométernél nagyobb részecskék extrathoratikus kiülepedési hatásfokát mutatja, mint a részecskeméret és a térfogatáram függvénye. Az eredményeket az irodalmi mérési adatokkal [Cheng 1999] is összehasonlítottam és - amint azt az ábra is tanúsítja - jó egyezést találtam. 70

Megállapítható, hogy a kiülepedési hatásfok minden térfogatáramra a részecskeméret növekedésével nő. Ez annak tulajdonítható, hogy az inerciális impakció (jelen esetben mint fő depozíciós mechanizmus) az átmérővel négyzetesen növekszik. Megfigyelhető továbbá egy telítődési jelenség is, ami már térfogatáram függő és nagy sebességeke már kisebb átmérőknél is megnyílvánul. Az ábra alapján a térfogatáram növekedése egyértelműen a kiülepedési hatásfok növekedéséhez vezet, ami szintén az impakciós kiülepedés dominanciájával magyarázható. Mivel a kisérleti mérésekhez csak mikrométernél nagyobb részecskéket használtak az ábrán nem tüntettem fel a szubmikronos tartományt. E tartományra vonatkozó eredményeim is jól egyeztek a mások által [Zhang 2005] hasonló geometriában számoltakkal.

100

Kiülepedési hatásfok (%) (%) Deposition efficiency

80

60

40

Experiment l/min60 l/perc Kisérlet [Cheng60 1999], Kísérlet Experiment l/min30 l/perc Kisérlet [Cheng30 1999], Kísérlet Experiment l/min15 l/perc Kísérlet Kisérlet [Cheng15 1999], Simulation 60 60l/perc l/min Szimuláció Simulation 30 30l/perc l/min Szimuláció Simulation 15 15l/perc l/min Szimuláció

20

0 0

5

10

15

20

25

30

35

40

45

Particle diameter ((µm) µm) Részecskeátmérő II.7. ábra Számított és mért oro-pharyngeális kiülepedési hatásfokok összehasonlítása különböző térfogatáramok esetében. A szimulált belépő sebességprofil egyenletes.

71

Az orrlégzés során történő kiülepedési hatásfokokat a sztochasztikus tüdőmodellel [Koblinger 1990] számítottam ki, mivel a nasopharingeális geometria sokkal bonyolultabb

és

orvosi

képalkotó

módszerekkel

nyert

metszetekből

történő

háromdimenziós rekonstrukciója csak most van folyamatban. Így a közeljövőben várhatóan itt is numerikus áramlástan alapú módszerekkel számított egzakt kiülepedési eredmények állnak majd rendelkezésre.

9.2. A részecskeméret hatása a bronchiális kiülepedésre A belélegzett részecskék centrális légúti kiülepedésének részecskeméret függését a jobb felső lebeny 4-5. generációs elágazásában (ahol a légcső jelenti az első generációt) vizsgáltam. Az áramlási tér sebesség értékeinek ismeretében követtem a belélegzett részecskék pályáit és kiülepedési helyét. A számításokat széles részecskeméret

Kiülepedési hatásfok (%)

tartományra (1 nm-25 µm) végeztem. Az eredményeket a II.8. ábra tünteti fel.

Részecskeátmérő (µm) II.8. ábra Részecskeméret függő kiülepedési hatásfokok a tüdő jobb felső lebenyének 45. generációs elágazásában. A bemeneti sebességprofil parabolikus, a térfogatáram 7.5 l/perc. 72

Látható, hogy a nanométeres részecskék kiülepedési hatásfoka csökken a részecskeméret növekedésével. Ennek oka az, hogy ebben a mérettartományban a diffúziós kiülepedés a legjelentősebb és a részecskék diffuzivitása csökken az átmérő növekedésével. A 100 nm és 1 µm közötti tartományba eső részecskék kiülepedése a legkevésbé intenzív és ezeknek a részecskéknek van a legnagyobb esélyük arra, hogy elérjék a tüdő mélyebb régióit. A mikron feletti részecskék kiülepedési hatásfoka a felsőlégúti kiülepedés tanulmányozásánál is leírt okok miatt nő az átmérő növekedésével.

9.3. A térfogatáram hatása a bronchiális kiülepedésre A fizikai tevékenység, emocionális állapot stb. függvényében változhat az egy légzési periódus alatt belélegzett levegő mennyisége (légzési térfogat) és az egységnyi időre jutó ki-belégzések száma (légzési frekvencia). E két légzési paraméter egyértelműen meghatározza a térfogatáramot a trachea szintjén. A fenti 4-5. generációs centrális légúti geometriában különböző részecskeméretekre számított kiülepedési hatásfokokat a térfogatáram függvényében a II.9. ábra mutatja. A felső panel a nano, az alsó panel pedig a mikrométeres tartomány kiülepedési hatásfokát szemlélteti. A részecskeméret tartomány mikrométernél történő szétválasztását az indokolja, hogy míg a nagy részecskék depozíciós hatásfoka nő a térfogatáram növekedésével, addíg a kis részecskéknél ennek az ellenkezője figyelhető meg. Az ok ismét a domináns depozíciós mechanizmusokra vezethető vissza. A 10 nm alatti átmérőjű részecskék esetében a Brown diffúzió dominál és a diffúziós kiülepedés annál erősebb minél több időt tartózkodnak a részecskék a légúti elágazásban. Ez utóbbi fordítottan arányos a térfogatárammal. A 10 nm és 1µm közötti tartományt a diffúzió és az impakció versenye jellemzi. A hatásfokgörbék is e tartományon belül metszik egymást. A nagy részecskék impakciós kiülepedése a térfogatáram növekedésével nő. A legnagyobb részecskék esetében a gravitációs kiülepedés is jelentős lehet.

73

Kiülepedési hatásfok (%) Deposition efficiency (%)

10

18 18 l/perc l/min 30 30l/perc l/min 60 60 l/perc l/min

1

1

10

100

1000

Kiülepedési hatásfok (%)

Particle diameter (nm) Részecskeátmérő (nm)

18 l/perc 30 l/perc 60 l/perc

Részecskeátmérő (µm) II.9. ábra Belélegzett részecskék kiülepedési hatásfoka különböző térfogatáramok esetén a tüdő jobb felső lebenyének 4-5. generációs elágazásában. A bemeneti sebességprofil parabolikus. A térfogatáramok a légcsőre vonatkoznak.

9.4. A légúti morfológia hatása a bronchiális kiülepedésre A légúti levegőáramok tanulmányozásakor kiderült, hogy az áramlási tér sebességértékei érzékenyek a geometriára. Várható, hogy a megfelelő kiülepedés szintén morfológiafüggő lesz. Ennek vizsgálatára összehasonlítottam a részecske kiülepedést az

74

előbbi 4-5. generációs elágazásban és egy ugyanakkora anyaág átmérőjű, de erősen aszimmetrikus geometriában (II.10. ábra).

dp = 1 nm η = 16,5 %

dp = 1 nm η = 15,6 %

dp = 10 µm η = 59,6 %

dp = 10 µm η = 66,6 %

II.10. ábra Részecske kiülepedés eloszlások szimmetrikus (bal oldal) és aszimmetrikus (jobb oldal) centrális légúti elágazásokban. A bejáraton sorsolt részecskék száma 50 000, a térfogatáram 7,5 l/perc. dp-részecskeátmérő; η - kiülepedési hatásfok. Mint látható, a kis és nagy részecskék kiülepedési hatásfoka egyaránt kissé különböző a két geometriára. A nanorészecskék esetében a különbség kisebb, mivel a diffúziós kiülepedés dominál, amely egyenletesebb kiülepedéshez vezet, mint az impakciós. Mivel az aszimmetrikus elágazás felülete kisebb, itt valamivel kevesebb részecske ülepszik ki. A nagy részecskéknél az aszimmetrikus geometria nagyobb elágazási szögének köszönhető az intenzívebb részecske kiülepedés. Az ábráról az is látható, hogy nem csak a kiülepedett részecskék száma, de a lokális kiülepedés is különböző. Mindez azt jelenti, 75

hogy a kiülepedés geometria függő és a pontos modellezéshez minél realisztikusabb geometriát kell alkalmaznunk. Ennek kapcsán felmerül, hogy fontos lenne a modellt továbbfejleszteni, hogy egyenes és merev falak helyett lokálisan nem sima és flexibilis falakkal számoljon. Ennek azonban egyelőre komoly akadályai vannak, ugyanis addig míg a lokális gyűrűsség, ráncosság stb. karakterisztikus mérete tizedmiliméter vagy az alatti a jelenlegi orvosi képalkotó eszközök rétegtávolsága miliméter körüli, tehát egyelőre ezek nem szolgáltatnak érdemi információt ilyen tekintetben. A légutak szisztematikus

mozgásának

modellezése

minden

egyes

elmozduló

rácspont

mozgásegyenletének figyelembe vételét jelentené. Ez a mai számítógép kapacitásokat figyelembe véve még két dimenziós modellnél is jelentősen megnövelné a gépidőt, de a processzorok és memóriák jelen fejlődési ütemét látva nem tűnik irrealisztikusnak, hogy öt-tíz év múlva a falak dinamikáját is modellezhetjük majd.

9.5. Részecskekiülepedés belégzéskor és kilégzéskor A be- és kilégzéshez tartozó részecskekiülepedés összehasonlításához a fenti megállapítások alapján olyan geometriát választottam, amely a radonszármazékok sugárterhelése és hatása szempontjából a legtöbb információt szolgáltatja. Az irodalom szerint [Veeze 1968] az uránbányászok esetében a tüdőrák a legnagyobb gyakorisággal a jobb felső lebeny 3-5. légúti generációiban keletkezett, ezért a kiülepedést egy olyan geometriában számoltam, amely ezen részeket is magába foglalja. Mivel a belépő levegősebesség és részecske eloszlás függ az előző légúti szegmenstől is, az említett légúti generációkhoz a légcsövet és a főbronchusokat is hozzácsatoltam. A különböző nagyságú részecskék be- és kilégzési kiülepedési hatásfokát a II.11. ábra mutatja. Látható, hogy a be- és kilégzési kiülepedési hatásfok értékek igazán csak 15 µm-nél nagyobb részecskék esetében térnek el jelentősen. További vizsgálataim azt mutatják, hogy nagy részecskékre a lokális kiülepedés merően más be- és kilégzésnél [Farkas 2005c]. Addig amíg belégzésnél az elágazások csúcsán kívül a leányágak belső felületén nagy a kiülepedési sűrűség, addíg a kilégzésnél az anyaág felső és alsó felülete a leginkább terhelt (lásd a II.12. ábrát).

76

Kiülepedési hatásfok (%) Deposition efficiency (%)

50

Belégzés Inhalation Exhalation Kilégzés

40

30

20

10

1E-3

0.01

0.1

1

10

Particle diameter (µm)

Részecske átmérő (µm)

II.11. ábra Be és kilégzési kiülepedési hatásfokok a részecskeátmérő függvényében. A légcsőre vonatkoztatott térfogatáram 18 l/perc, a bemeneti sebességprofil parabolikus. η=49%

g

η=33,1% dp = 25 µm

g

Kilégzés

Belégzés

II.12. ábra Nagy részecskék be (bal oldal) és kilégzési (jobb oldal) kiülepedés eloszlásai. A bemeneten (kilégzéskor a hat bemeneten) sorsolt részecskék száma 50 000, a térfogatáram a légcsőnél 18 l/perc. A nyilak a gravitáció irányát mutatják. dp részecskeátmérő; η - kiülepedési hatásfok. 77

9.6. Regionális depozíció A jobb felső lebeny első öt elágazását magába foglaló geometriában és a felsőlégutakban számolt depozíciós hatásfokokat összesítve lehetőség van az e rendszerre vonatkoztatott regionális és teljes depozíció jellemzésére. Ez egyrészt azért hasznos mert láthatóvá válik, hogy a különböző méretű részecskék a légutak mely részén ülepednek ki, másrészt az eredmények összehasonlíthatóak a korábbi analitikus modellek regionális depozíciós eredményeivel. Kiülepedési hatásfok helyett a regionális kiülepedést depozíció frakcióval szokták jellemezni, ami az adott tartományban kiülepedett és a belélegzett részecskék számának az aránya. A regionális és teljes kiülepedett frakció

Kiülepedett Deposited frakció

értékeket a II.13. ábrán tüntettem fel.

0.7

extrathoratikus

extrathoratikus

0.6

tracheobronchiális teljes (gen. 1-5)

tracheobronchiális teljes (gen 1-5)

0.5 0.4

15 l/perc

60 l/perc

0.3 0.2 0.1

1 nm

1 µm

10 µm

1 nm

1 µm

10 µm

II.13. ábra 1 nm, 1 µm és 10 µm-es részecskék felsőlégúti, centrális légúti és teljes kiülepedett frakciói 15 l/perc és 60 l/perc térfogatáramú szájlégzés esetén. A II.13. ábrán szemléltetetteket összesítve elmondható, hogy a néhány nanométeres és több mikronos részecskék diffúziós, illetve impakciós kiülepedés során kitapadnak a 78

felső és centrális légutak falára és csak kis hányaduk jut le az ötödik generációnál mélyebb régiókba. Ezzel ellentétben a közepes méretű részecskék nagy eséllyel az alveolusokat fogják terhelni. Ez az eredmény összhangban van több korábbi számítás eredményével, amit analitikus modellekkel számítottak (pl. [Koblinger 1990]). A CFD alapú modellek erőssége viszont a lokális kiülepedés jellemzésénél mutatkozik meg, amire az analitikus modellek nem képesek. A numerikus elágazás-szintű és annál kisebb skálán végzett jellemzést a továbbiakban mutatom be.

9.7. Légúti elágazásonkénti depozíció sűrűség Komplex, sokgenerációs légúti geometriák alkalmazásával pontosabban kiszámítható a kiülepedés, mintha csak egyetlen elágazást vennénk figyelembe. A kapott kiülepedési hatásfok viszont a teljes rendszert és nem az egyes elágazásokat jellemzi. Az elágazásonkénti hatásfok számítására egy PERL nyelven írt saját programot használtam. A program mindegyik kiülepedett részecskét hozzárendeli valamelyik elágazáshoz majd elágazás-szintű kiülepedési hatásfokokat számít. A program segítségével a II.11. ábrán látható globális be és kilégzési kiülepedési hatásfokokat elágazásonkénti depozíciós hatásfokokká alakítottam. Mivel a nagy felülettel rendelkező elágazásokon több részecske tapad meg, ezek kiülepedési hatásfoka is nagyobb. Az egészségügyi hatások szempontjából viszont több információt nyújthat a felületegységre vett kiülepedési hatásfok. A kiülepedési hatásfokot a megfelelő elágazás területével osztva megkapjuk, hogy az adott elágazás egységnyi felületén a bemenő részecskék hány százaléka ülepszik ki. Az eredmény a II.14. ábrán látható. Mint az már az áramlási tér ismeretében várható volt, mint belégzésnél mint kilégzésnél jelentős különbségeket tapasztalunk az egyes elágazások terheltségét illetően.

79

Kiülepedési hatásfok sűrűség (%/mm2)

0.006

gen 1-2

gen 4-5.1

gen 2-3

gen 4-5.2

gen 3-4 0.004

1 nm

1 nm

25 µm

25 µm

0.002

Belégzés

Kilégzés

Belégzés

Kilégzés

II.14. ábra Be és kilégzési elágazás-specifikus kiülepedési hatásfok sűrűség a II.12. ábrán is látható centrális légúti geometriában. A légcsőre számított térfogatáram 18 l/perc.

9.8. Lokális depozíció Mint azt a II.12. ábrán látható kiülepedési eloszlások is mutatják, a depozíció eloszlás egyetlen elágazáson belül sem egyenletes. A legerősebb expozíciót minden esetben az elágazások nyereg alakú csúcsában láthatjuk (lásd például [Farkas 2002a, Balásházy 2002b,c,e Balásházy 2004a vagy Farkas 2003a] munkáinkat). Ez megegyezik az uránbányászoknál megfigyelt, radon leányelemek által indukált malignáns tumorok keletkezési helyével [Churg 1996]. Nyilvánvaló tehát, hogy a kockázatbecslés alapja sokkal inkább a lokális mint az átlagos terhelés kell, hogy legyen. A lokális kiülepedés kvantifikálásához bevezettem egy ún. fokozott kiülepedési tényezőt (FKT), ami azt mutatja, hogy a lokális kiülepedési sűrűség hányszorosa az átlagosnak. Lokális kiülepedési sűrűség alatt egy kis felületre kiülepedett részecskék számának és az adott felület területének arányát értjük. Az átlagos kiülepedési sűrűséget úgy számítjuk ki,

80

hogy a tanulmányozott geometriában kiülepedett összes részecskék számát elosztjuk a teljes felület területével.

dp= 10 µm

212 0

FKT Maximum enhancement Maximális FKT factor

0

dp = 1 nm

129

FKT

2000 1800

Q60 = l/perc 60 l/min

1600 1400 1200 1000 800 600 400 200 0 1E-3

0.01

0.1

1

10

Particle diameter ((µm) µm) Részecskeátmérő Részecske átmérő

II.15. ábra 10 µm (bal felső panel) és 1 nm (jobb felső panel) aerodinamikai átmérőjű részecskék fokozott kiülepedési tényezői (FKT) és a maximális fokozott kiülepedési tényező részecskeméret függése (jobb alsó panel) 60 l/perc tracheális térfogatáram esetében, a felnőtt tüdő 4-5. generációs (ha a légcső az első) elágazásában. Az elemi felület oldalhossza 45 µm. dp - részecskeátmérő.

81

Értelemszerűen, a depozíció inhomogenitása miatt a fokozott kiülepedési tényező attól is függ, hogy mekkora lokális elemi felületet választunk ki. Ennek eldöntésekor a sejtszintű terhelést tartottam szem előtt. A szakirodalmi adatok alapján egy bronchiális epithel sejt nagysága durván 10 x 10 µm [Mercer 1991]. Tekintettel a sejtek közötti kommunikációra (lásd bystander hatás, [Hall 2003]) néhány tíz sejtből álló csoport területének megfelelő nagyságú elemi felület kiválasztása tűnik indokoltnak. Ezen elemi felülettel végigpásztáztam a teljes geometrián és egy saját program segítségével kiszámítottam a megfelelő fokozott kiülepedési tényezőket. A II.15. ábrán az elemi felületek húsz sejtnyi területű háromszögek. Látható, hogy addig míg a sejtcsoportok egy részére nem ülepszik ki egyetlen részecske sem (sötét kék színű háromszögek), addíg mások az átlagosnál két vagy akár három nagyságrenddel nagyobb terhelést kapnak. A legnagyobb fokozott kiülepedési tényezővel jellemzett helyek minden esetben a karina régióba esnek, pont ott ahol az uránbányászok esetében a pre-neoplasztikus és neoplasztikus károsodásokat észlelték. Szembetűnő továbbá, hogy a maximális fokozott kiülepedési tényező kisebb a nanorészecskék esetében. Ez minden bizonnyal a diffúziónak köszönhető, amely úgymond szétszórja a nagyon kis részecskéket. A maximális fokozott kiülepedési tényező valahol az 1 µm-es részecskék körül a legnagyobb, mivel e részecskékre az átlagos depozíció sűrűség igen kicsi.

9.9. Időfüggő részecske kiülepedés A stacioner számításokból kiderült, hogy a kiülepedési hatásfok részecskeméret függő. A nanométeres részecskék diffúziós és a mikrométernél nagyobb átmérőjű részecskék impakciós kiülepedése jóval hatékonyabb, mint a 0,1-1 µm tartományba eső részecskéké. Az időfüggő kiülepedés tanulmányozására e három tartományból választottam egy-egy részecskeméretet (1 nm és 1 illetve 10 µm). A megfelelő kiülepedési hatásfokok időfüggését a II.16. ábra mutatja.

82

Térfogatáram (l/perc)

100

50

0

-5 0

-1 0 0

0 .0

0 .5

1 .0

1 .5

2 .0

2 .5

3 .0

3 .5

4 .0

3 .5

4 .0

100

Kiülepedési hatásfok (%)

90

1 0 µm µm 10 11µµm m 11nnm m

80 70 60 50 40 30 20 10 0 0 .0

0 .5

1 .0

1 .5

2 .0

2 .5

3 .0

Idő (s)

II.16. ábra 1 nm-es, valamint 1 és 10 µm-es részecskék kiülepedési hatásfokának időfüggése egy 4-5. generációs centrális légúti elágazásban. A térfogatáram a bemeneten 7,5 l/perc, a légzési periódus 4 s.

Látható, hogy a kiülepedési hatásfok itt is részecskeméret függő. Továbbá, az inhalációs kiülepedési hatásfok minden részecskeméretre nagyobb az exhalációsnál, a különbség a mikrométeres részecskéknél a legkevésbé jelentős. E részecskék belégzési és kilégzési kiülepedési hatásfoka egyaránt alacsony. Ezért a becslés relatív hibája itt a legnagyobb. A nanorészecskék kiülepedése a félperiódusok közelében a legintenzívebb. Ennek oka valószínűleg az, hogy ilyenkor a levegő sebessége lecsökken és ezáltal a 83

részecskék diffúzivitása (és diffúziós kiülepedése) megnő. Mint látható, a kiülepedés maximuma a félperiódustól kicsit jobbra tolódik. Ez annak köszönhető, hogy a kis levegősebesség miatt a részecskék rezidens ideje nagy, ezért egy részük később ülepszik ki. A 10 µm aerodinamikai átmérőjű részecskék kiülepedése a legintenzívebb. A levegő sebességének növekedésével a kiülepedési hatásfok is nő és belégzésnél a negyedperiódus, kilégzésnél pedig a háromnegyed periódus környékén a legnagyobb. Ez az impakcióval magyarázható, amely domináns kiülepedési mechanizmus ennél a részecskeméretnél és amely arányos a sebességgel.

84

10. Kiülepedett radonszármazékok lokális terhelésének modellezése uránbányáknak és lakásoknak megfelelő expozíciós környezetben Amint a bevezetőben már említettem, a radoninhalációt követő kockázat becslésében fontos szerepet játszanak az uránbányákra és lakásokra vonatkozó epidemiológiai adatok. A kis dózisok tartományában viszont ezek statisztikája nagyon szór, az eredmények bizonytalansága nagy. A számítási modellek fő hibája, hogy a légutak falán egyenletes aktivitás eloszlásokat feltételeznek. Ezzel ellentétben az általam alkalmazott modell kimutatta, hogy a radionuklidok tracheobronchiális kiülepedése erősen inhomogén. Az inhomogeneitás következményeinek vizsgálatához e jelen fejezetben a szakirodalomból vett bányára illetve lakásra vonatkozó expozíciós adatokra végzek modellszámításokat. Az expozíciós körülmények lakásról-lakásra, illetve bányáról bányára különböznek. Ezenkívül az expozíciós paraméterek ugyanazon bánya, vagy lakás esetében is időrőlidőre változnak. Ezért nehéz tipikus lakásról vagy bányáról beszélni. Mindennek ellenére sok mérést figyelembe véve expozíciós átlagértékeket számíthatunk ki. A lakások esetében ezt könnyebb megtenni, mivel sokkal több adat áll rendelkezésre. Az uránbányák egykori légteréről viszonylag kevés adat áll rendelkezésre. Ezért amíg lakások esetében a világ több mint 2000 különböző helyén mért aktivitás koncentráció adatok átlagát használtam, addíg a bányáknál a viszonylag jól dokumentált Új-Mexikói uránbányát vettem alapul. A számítási modell a konkrét mérési adatok ismeretében bármilyen expozíciós környezetre alkalmazható. Az irodalomból vett expozíciós adatokat a II.1. táblázat tartalmazza. Az ICRP 66 [ICRP66 1994] publikáció szerint nyugalmi helyzetben, de még könnyű fizikai munkánál is orrlégzés a jellemző. A lakás esetében nyugalmi, míg bányászok esetében a könnyű fizikai munka légzési paramétereivel számoltam, mintegy kiátlagolva a pihenés és a nehéz fizikai munka paramétereit. Bár a tanulmányozott bányák közül az Új-Mexikóira találtam a legtöbb adatot, konkrét részecskeméret eloszlást itt sem adtak meg. Ezért az [NRC 1991] által is alapul vett bányákat jellemző részecskeméret eloszlás két csúcsának megfelelő aktivitás átmérőkkel dolgoztam [Cooper 1973]. Számításaim során kitapadt radionuklidként a radon három rövid felezési idejű bomlástermékét (218Po, 214Pb és 214Bi) vettem figyelembe

85

II.1. táblázat A modellszámítások során használt irodalmi expozíciós adatok.

Lakás

Uránbánya

Adat

Forrás

Légzés típusa: Térfogatáram: Légzési frekvencia: Részecskeátmérő: - kitapadt - ki nem tapadt Ki nem tapadt hányad: Izotóp aktivitás koncentráció arány (218Po/214Pb/214Bi) Potenciális alfa-energia koncentráció

orrlégzés 50 l/perc 20/perc 250 nm 1nm 1% 0,6/0,29/0,21 5,7 WL

[BEIR VI 1999]

Légzés típusa: Térfogatáram: Légzési frekvencia: Részecskeátmérő: - kitapadt - ki nem tapadt Ki nem tapadt hányad: Izotóp aktivitás koncentráció arány (218Po/214Pb/214Bi) Potenciális alfa-energia koncentráció

orrlégzés 18 l/perc 12/perc 200 nm 1,2 nm 6% 4/3/2

[ICRP 66 1994] [ICRP 66 1994] [ICRP 66 1994] [Haninger 1997] [Reineking 1988] [Haninger 1997] [Unscear 2000]

0,0072 WL

[ICRP 66 1994] [ICRP 66 1994] [ICRP 66 1994] [Cooper 1973] [NRC 1991] [Samet 1989]

[Unscear 2000]

A gáznemű radon légúti terhelését a viszonylag nagy felezési idő miatt elhanyagoltam. A 214Pb és 214Bi izotópok bár nem alfa-aktívak (lásd a II.17. ábrát), a tüdő falára kitapadva viszonylag rövid idő alatt 214Po-re bomlanak, ami viszont alfa-részecskét bocsájt ki. Ez utóbbi gyakorlatilag azonnal elbomlik, ezért légúti transzportja és kiülepedése modellezésének nincs értelme. Az eleme a 210Po, amely a 214Po-ből a

210

238

U bomlási sorának további alfa-bomló

Pb és 210Bi-en át jön létre amelyek együttes felezési

ideje közel 21 év, saját felezési ideje pedig 138,4 nap. A tüdő gyors tisztulási mechanizmusait figyelembe véve [Balásházy 2003b,e Szőke 2002a,b,c,e,f, 2003a] ezen izotóp terhelése is elhanyagolható. Modellszámításaimban a ki nem tapadt hányadot teljes mértékben a 218Po jelenti. Ennek a közelítésnek az alapját az a tény képezi, hogy a többi ki nem tapadt izotóp a ki nem tapadt hányad aktivitásának csupán a tizedét adja [UNSCEAR 2000]. Ez valószínűleg azért van így, mert a 218Po nagy része még elbomlása előtt kitapad a levegőben található aeroszolokra.

86

6 MeV 5.5 α

222

Rn

3.823 nap 3.823 d

5.49 MeV 5.49 α 218 γ

Po

66 MeV MeV α

Po

(stabil)

214

Pb

26.8 26.8perc

3.05 3.05perc

206

0.020,02 %

5.3 5.3MeV α γ

γ β

214

Bi

19.7 19.7perc

210

Tl

99.8 99.98% γ β

1.32 1.32perc

214

Po

β

164 µs 164 µs

γ

7.697.69 MeV

210

Po

β

138.4 138.4nap

210

Bi

5.01 nap 5.01

γ β

α

210

Pb 21 év

II.17. ábra Az 238U bomlási sora a 222Rn-től kezdődően. Statisztikai okokból indokolt minél nagyobb számú inhalált részecske pályáját követni. Ez jelen esetben tízmillió belélegzett részecskét jelentett úgy lakásra mint bányára, amelynek trajektória-számítási gépidő nagyságrendje a hónap. Ebből a légcsőbe jutó részecskék számát a felsőlégúti kiülepedést figyelembe véve kapjuk meg. A sztochasztikus tüdőmodellel [Koblinger 1990] végzett számítások a bányára jellemző légzési paraméterek mellett 3,3 és 82 % -os felsőlégúti kiülepedést adtak a kitapadt illetve ki nem tapadt frakcióra. Ugyanezek a felsőlégúti kiülepedési hatásfokok lakásra 4,37 és 89,8 %. A belélegzett és a légcsőbe jutó részecskék számát kitapadt és ki nem tapadt részecskeszámra bontva a II.2. táblázat adja meg. Ugyancsak a II.2. táblázat tartalmazza II.1 táblázatban feltüntetett aktivitás koncentráció arányoknak és potenciális alfa-energia koncentrációknak megfelelő izotóp szám és aktivitás koncentrációkat. A részecske trajektória számításokat a II.12. ábrán is látható centrális légúti geometriában végeztem. A kiülepedési képek és hatásfokok a II.18. ábrán láthatók. Az ábra szerint a radioaktív részecskék kiülepedése úgy bányára mint lakásra inhomogén. A ki nem tapadt hányad kiülepedése intenzívebb lakásra mint bányára. Ez egyrészt annak köszönhető, hogy a lakás légterében nagyobb arányban van jelen ez a hányad, másrészt a kisebb sebességek miatt a diffúziós kiülepedés lakásra erősebb.

87

II.2. táblázat Rövid felezési idejű radon leányelemek aktivitás és számkoncentrációja, valamint a belélegzett és a légcsőbe jutó izotópok száma a) bánya és b) lakás esetén. a) Uránbánya Izotóp

218

Po Pb 214 Bi Összes 214

Aktivitás konc. (Bq/m3)

Szám konc. (l-1)

Belélegzett részecskék száma

43 456 21 004 15 207

11 474 48 735 25 941 86 150

1 332 000 5 657 000 3 011 000 10 000 000

Aktivitás konc.. (Bq/m3)

Szám konc. (l-1)

Belélegzett részecskék száma

40 30 20

10,36 69,61 34,13 114,1

908 000 6 101 000 2 991 000 10 000 000

Belélegzett kitapadt részecskék száma 1 232 000 5 657 000 3 011 000 9 900 000

Légcsőbe jutó részecskék száma 1 219 344 5 470 319 2 911 637 9 591 300

Légcsőbe jutó kitapadt részecskék száma 1 191 344 5 470 319 2 911 637 9 573 300

Belélegzett kitapadt részecskék száma 308 000 6 101 000 2 991 000 9 400 000

Légcsőbe jutó részecskék száma 355 740 5 834 386 2 860 293 9 050 419

Légcsőbe jutó kitapadt részecskék száma 294 540 5 834 386 2 860 293 8 989 219

b) Lakás Izotóp

218

Po Pb 214 Bi Összes 214

Az egészségügyi hatások szempontjából itt is fontos a lokális kiülepedés jellemzése [Farkas 2002b, Balásházy 2002d]. A helyi kiülepedést kvantáló fokozott kiülepedési tényezőket az előző fejezetben ismertetett módon számítottam ki. A II.19. ábra ezek eloszlását mutatja bányára és lakásra. Látható, hogy a lakásnak megfelelő expozíciós feltételek mellett a körülbelül 50 000 cella (a 75 000-ből) egyáltalán nem kap részecskét, ugyanakkor lesz az átlagoshoz képest több mint 300-szorosan terhelt elemi terület is. Uránbányára jellemző körülmények között több mint 60 000 elemi cellára egyetlen részecske sem ülepszik ki, de néhányra az átlagosnak több mint 400-szorosa.

88

a) ki nem tapadt

kitapadt

Q = 50 l/perc η = 27,8 %

Q = 50 l/perc η = 0,94 %

b) ki nem tapadt

kitapadt

Q = 18 l/perc η = 36,88 %

Q = 18 l/perc η = 0,73 %

II.18. ábra Kitapadt (218Po,

214

Pb és

214

Bi) és ki nem tapadt (218Po) radon leányelemek

centrális légúti kiülepedési eloszlásai a) uránbányában és b) lakásban uralkodó expozíciós körülmények mellett. A belélegzett részecskék száma mindkét esetben tízmillió. Q – térfogatáram; η − kiülepedési hatásfok. Mivel a kis és nagy részecskék szövetbehatoló képessége különböző, érdemes a maximális fokozott kiülepedési tényezőt a kitapadt és ki nem tapadt hányadra külön is kiszámítani. A számítások eredményeit a II.20. ábra mutatja. Mint az várható volt, a ki nem tapadt hányad maximális fokozott kiülepedési tényezői úgy bányára mint lakásra kisebbek a kitapadt hányadra jellemző értékeknél. Továbbá, míg a kitapadt hányad esetében a bányára és lakásra jellemző értékek közötti eltérés nem nagy, addig a ki nem tapadt hányad esetében a bányai a lakásra jellemző értéknek több mint kétszerese.

89

Gyakoriság

Bánya (Q = 50 l/perc)

FKT

Gyakoriság

Lakás (Q = 18 l/perc)

FKT II.19. ábra Fokozott kiülepedési tényezők (FKT) eloszlása lakás (fent) és uránbánya (lent) esetén. Az elemi felület területe 1,375 x 10-7 m2; Q – térfogatáram.

90

EF max max FKT

Bánya 487

500

Lakás 461,1

460

400

343,4 300 200 100

64,2 25,9

ki nem tapadt

II.20. ábra Kitapadt (218Po,

kitapadt

214

Pb és

214

összes

Bi) és ki nem tapadt (218Po) radon

bomlástermékek maximális fokozott kiülepedési tényezői (FKTmax) uránbányára és lakásra. A fentiek alapján egyértelmű, hogy az egyenletes kiülepedés eloszlás feltételezése a mikrodozimetriai paraméterek hibás becsléséhez vezethet. Ugyanakkor a sejtszintű terhelés nem csupán a kiülepedett radionuklidok számától, de az egyes nuklidok által kibocsájtott alfa-részecskék energiájától is függ. Míg a kiülepedett 13,69 MeV energiával fogja terhelni a hámsejteket, addíg a

214

218

Po potenciálisan

Pb és a

214

Bi izotópok

csupán 7,69 MeV-al. A különböző izotópok potenciális alfa-energia kibocsájtó képességét egy saját program segítségével vettem figyelembe, amellyel így sikerült kvantifikálnom a lokális potenciális alfa-energiát, vagyis a tüdőhámsejtek radon bomlástermékek okozta lokális terhelését Ezt az ún. fokozott potenciális alfa-energia tényezőkkel (FPAT) fejeztem ki, amelyek az egységnyi felületre jutó lokális és átlagos potenciális alfa-energia hányadát adják meg. Ezek maximális értékeire bányára 480,5,

91

lakásra pedig 325,82 adódott, vagyis a leginkább terhelt sejtcsoportokra az átlagosnál ennyiszer több energia jut bányában illetve lakásban. A fenti eredmények elemzésével néhány fontos következtetésre juthatunk. Ma már a legtöbb szakmabeli belátja, hogy nem helyes a lakásokra jellemző dózis–hatás összefüggést az uránbányász adatok direkt extrapolációjával meghatározni. Ennek egyik oka az, hogy a dózisteljesítmény a két környezetben nagyon különböző. Az általam vizsgált esetben például, míg a tízmillió radioaktív részecskét lakásban körülbelül egy hét alatt inhaláljuk, addig bányában valamivel több mint négy perc alatt (a kettő között egy 2094-es faktor van). Ezt bizonyos korrekciós faktorok alkalmazásával szokták kiküszöbölni [Cohen 2000]. Amint arra a jelen számítások rámutatnak, ezen túlmenően azt is figyelembe kell venni, hogy ugyanaz a makroszkópikus dózis, mást terhelést jelent bányában mint lakásban. Jelen esetben például a tízmillió inhalált radioaktív részecske összesen 70 243 080 Mev potenciális alfa-energiát jelent bánya és 66 982 520 MeV-ot lakás esetében. Vagyis ugyanakkora makroszkópikus dózisra a biológiai válasz várhatóan még egyenletes kiülepedést feltételezve is más lesz bányában mint lakásban. Erre rátevődik még az, hogy a leginkább terhelt sejtcsoportoknak az átlagos terheléshez viszonyított aránya (maximális fokozott potenciális alfa-energia tényezője) más bányára mint lakásra. Ez utóbbi szempontot nem veszi figyelembe semmilyen eddigi modell. Továbbá, jelen modellszámítások eredményeit figyelembe véve, legalábbis kis dózisokra, megkérdőjelezhető a dózis-hatás görbe unicitása. Ugyanis, ha figyelembe vesszük a radio-aeroszolok kiülepedésének inhomogenitását, akkor egyes sejtcsoportok az átlagosnál jóval nagyobb dózist kapnak és a legtöbb sejtcsoport az átlagosnál kevesebbet. A hisztológiai tanulmányok szerint a rák az elágazások csúcsában kezdődött. Ez a hely megegyezik a számítások szerint leginkább terhelt sejtcsoportok helyével. Ez azt sugallja, hogy a kockázat becslésnél elsősorban ezeket a klasztereket kell figyelembe venni. A forró területek terhelése más és más lehet ugyanakkora makroszkopikus dózisra (a fokozott potenciális alfa-energia tényező ugyanannyi belélegzett részecskére más és más lehet pl a részecskeméret eloszlás, légzési mód stb függvényében), vagyis a lokális terhelés a dózison kívül más paraméterektől is függ. A mikroszkópikus terhelés attól függ, hogy valaki milyen egyéb paraméterekkel jellemzett helyen és mi módon kapja az adott dózist. Ez azt jelenti, hogy a forró területek sejtcsoportjaira jutó alfa-energia (és az

92

ebből származó mikroszkópikus dózis) egy sokváltozós függvény, aminek egyik paramétere a makroszkópikus dózis. A makroszkópikus dózis-lokális dózis konverziós görbének a többi paraméter függvényében más és más alakja lehet, mert ez tulajdonképpen egy sokdimenziós felületnek a dózisok síkjával való metszete. Ez azt is előrevetíti, hogy még a lokális terheléssel lineárisan változó biológiai válasz mellett is adott makroszkópikus dózisra több válasz lehetséges, vagyis több dózis-hatás görbe.

93

11. Kitekintés A modell jövőbeni továbbfejlesztése több irányban is lehetséges. Egyrészt a növekvő számítógép kapacitások egyre pontosabb számításokat engednek majd meg, egyre összetettebb geometriákon. Másrészt az alkalmazott numerikus eljárások tökéletesítése újabb minőségi javulást hozhat a levegő és részecskeáramlás modellezésében. Az orvosi képalkotó eszközök felbontóképességének javulása a jövőben lehetővé teszi majd a realisztikusabb légúti geometriák numerikus rekonstrukcióját. Sugárterhelés Makroszkópikus expozíciós adatok Numerikus áramlástan alapú részecsketranszport ás depozíciós modell

Lokális terhelés Alfanyomokat és sejttalálatokat számító modell (HPM) Mikrodozimetriai paraméterek Állapotvektor modell (SVM) Sejttranszformációs és rákkeletkezési valószínűségek Biológiai válasz

II.21. ábra A fejlesztésünk alatt álló biofizikai mechanizmusokon alapuló komplex kockázati modell vázlata.

94

Az általam számított lokális terhelés és az érzékeny hámsejtek eloszlásának ismeretében olyan mikrodozimetriai paraméterek számíthatók mint alfa-nyom hossza, találati valószínűség vagy sejtdózis [Farkas 2002c,2003b, Szőke 2003b,2004b,c, Balásházy 2003c, 2004b,c]. Ezekből a biológiai válasz a bevezetőben is említett mechanisztikus sejttranszformációs és rákkeletkezési modellekkel számítható ki. A pontos kockázatbecslés az általam kifejlesztett numerikus áramlástan alapú részecske kiülepedési modell és a most említett modellek jövőbeni integrálásával valósítható meg (II.21. ábra). Egy ilyen, biofizikai mechanizmusokon alapuló komplex kockázati modell felépítése (aktív részvételemmel) folyamatban van.

95

12. Összefoglalás A jelenlegi hivatalos sugárvédelmi álláspont szerint a biológiai kockázat lineáris a dózissal az ún. kis dózisok tartományában. A hipotézis helyességének elméleti vizsgálata magába foglalja a makroszkópikus és a sejtszintű terhelés közötti kapcsolatok feltárását, valamint a lokális terhelés biológiai válaszának vizsgálatát. E munka szerves részeként, jelen dolgozatban a belélegzett radio-aeroszolok légúti kiülepedésének numerikus modellezését és az expozíciós körülmények ismeretében a sejtcsoportok lokális terhelésének kvantifikálását tűztem ki célul. A kutatási téma felvetése és a kitűzött célok megfogalmazása után felépítettem egy fluidum és részecskedinamikai modellt a belélegzett levegő áramlásának és a légutakba jutó radioizotópok transzportjának, valamint kiülepedésének tanulmányozására. A modellfejlesztés főbb mozzanatai a numerikus légúti geometriák megszerkesztése és behálózása, az áramlási mező meghatározásához szükséges algoritmusok felépítése és a részecskék követését célzó numerikus eljárás kidolgozása voltak. A felépített modell kísérleti és szimulációs adatok segítségével történt validálása megmutatta, hogy a modell alkalmas a légúti levegőáramok és részecsketrajektóriák valósághű modellezésére. A kidolgozott és validált modellt alkalmazva jellemeztem a felső és centrális légutak áramlási terét. Megállapítottam, hogy a radon leányelemeket szállító levegő sebességtere érzékeny a légzési feltételekre és erősen függ a légutak morfológiájától. Az áramlási mező ismeretében megvizsgáltam az inhalált radionuklidok regionális, elágazás-specifikus és lokális légúti kiülepedését. Különböző felső és centrális légúti geometriákban tanulmányoztam a részecskeméret és a térfogatáram hatását a radioaeroszolok kiülepedésére. Azt találtam, hogy a részecskekiülepedés a légutakban inhomogén, a legintenzívebb az elágazások csúcsában. A kiülepedési hatásfok a nanométeres tartományban csökken, míg a mikrométeres tartományban nő a térfogatáram növekedésével. A depozíciós hatásfok csökken a részecskeméret növekedésével a nanométeres tartományban és nő annak növekedésével a mikrométeres mérettartományban. A modellt bányára és lakásra jellemző irodalmi expozíciós adatokra is alkalmaztam. A lokális kiülepedés kvantifikálása lehetőséget adott a két környezetnek megfelelő lokális terhelés összehasonlítására és a bányász adatoknak kisebb dózisokra történő extrapolálása során fellépő problémák feltárására. Rámutattam továbbá, hogy a lokális terhelés a dózison kívül több más paramétertől is függ, ami megkérdőjelezi a dózis-hatás görbe unicitását is.

96

13. Summary

The Linear-Nonthreshold (LNT) dose-effect theory presumes that the biological effects of radiation are proportional to the dose even in the low dose range. However, epidemiology and experiments do not provide unambiguous information regarding the health effects of low doses. Theoretical investigation of the validity of LNT hypothesis at low doses involves the clarification of relationships between the macroscopic and cellular burdens and the modeling of biological responses. As a part of these efforts, the primary objectives of this dissertation were to model the transport and deposition of the inhaled radio-aerosols within the respiratory tract and the quantification of local particle deposition and cellular radiation burdens. After a short introduction, a computational fluid and particle dynamics based model has been described for the characterization of airflow fields and particle transport within the airways. The model developments involved the region of the upper and central airways, the construction of numerical meshes and elaboration of flow field and particle trajectory describing numerical techniques. The model was validated by experimental and theoretical data of the published literature. By the application of the model, airflow fields in the extrathoracic and bronchial airways were characterized. The main statement of these investigations is that the airflow is very sensitive to the airway morphology and breathing conditions. Inspiratory and expiratory regional, generation specific and local deposition efficiencies have been computed as a function of flow rate and particle size. Highly inhomogeneous deposition patterns were found with primary hot spots at the carinal ridges of the bronchial airway bifurcations. The application of the model for exposure conditions characteristic to uranium mines and homes revealed that direct extrapolation of high dose health effects to low doses is not appropriate because of the different exposure conditions. Beside the dose, the local burden proved to be sensitive to several other parameters (particle size distribution, attached fraction, physical activity etc.), which suggests the existence of multiple doseeffect curves or rather multi-dimensional surfaces.

97

14 Irodalomjegyzék

Alavanja M C, Brownson R C, Lubin J H, Berger E, Chang J and Boice JD Jr 1994 Residential radon exposure and lung cancer among nonsmoking women J. Natl. Cancer Inst. 86 1829-37 Asgharian B, Hofmann W and Miller F J 2001 Mucociliary clearance of insoluble particles from the tracheobronchial airways of the human lung J. Aerosol Sci. 32 81732 Asgharian B, Menache M G and Miller F J 2004 Modeling age-related particle deposition in humans J. Aerosol Med. 17 213-24 Auerbach O, Stout A P, Hammond E C and Garfinkel L 1961 Changes in bronchial epithelium in relation to cigarette smoking and in relation to lung cancer N. Engl J. Med. 265 253-67 Auvinen A and Jokiniemi J K 2004 Modelling resuspension of particle layer in an internal tube flow J. Aerosol Sci. 1 S417-8 Baker A J 1983 Finite element computational fluid mechanics, Hemisphere, New-York Balásházy I, Hofmann W and Martonen T B 1990 Inertial impaction and gravitational deposition of aerosols in curved tubes and airway bifurcations Aerosol Sci. Technol.

13 308-21 Balásházy I and Hofmann W 1993a Particle deposition in airway bifurcations I. Inspiratory flow, J. Aerosol Sci. 24 745-72 Balásházy I and Hofmann W 1993b Particle deposition in airway bifurcations II. Expiratory flow, J. Aerosol Sci. 24 773-86 Balásházy I, Heistracher T and Hofmann W 1996 Airflow and particle deposition patterns in bronchial airway bifurcations. The effect of different CFD models and bifurcation geometries J. Aerosol Med. 9 287-301

98

Balásházy I, Hofmann W, Farkas Á and Szőke I 2002a Modelling carcinogenic effects of low doses of radon progenies J. Radiol. Prot. 22 89-93 Balásházy I, Farkas Á, Hofmann W and Kurunczi S 2002b Local deposition distributions of inhaled radionuclides in the human tracheobronchial tree Rad. Prot. Dos. 99 469-70 Balásházy I, Hofmann W and Farkas Á 2002c Numerical modelling of deposition of inhaled particles in central human airways Ann. Occup. Hyg. 46 353-7 Balásházy I, Hofmann W, Farkas Á, Pálfalvi J, Lohász M and Kristóf G 2002d Microdosimetric consequences of the inhomogeneity of the inhaled radionuclide deposition in the human lung. IRPA Regional Congress on Radiation Protection in Cenral Europe, Radiation Protection in Health, Dubrovnik, Croatia, Proceedings, 3p07, 1-6, ISBN 953-96133-3-7 Balásházy I, Hofmann W and Heistracher T 2003a Local particle deposition patterns may play a key role in the development of lung cancer J. Appl. Physiol. 94 1719-25 Balásházy I, Farkas Á, Szőke I, Hofmann W and Sturm R 2003b Simulation of deposition and clearance of inhaled particles in central human airways Rad. Prot. Dos. 105 129-32 Balásházy I, Farkas Á, Szőke I and Hofmann W 2003c Numerical simulation of alpha hit probability distributions in sensitive bronchial epithelial cells by inhaling radon progenies. 12th International Congress of Radiation Research, Brisbane, Queensland, Australia, Book of Abstract, 169 Balásházy I, Alföldy B, Osán J, Farkas Á, Szőke I and Török Sz 2003d Numerical simulation of the deposition of toxic elements originating from fossil burning in the human airway system, 12th International Conference on Fluid Flow Technologies, Conference on Modelling Fluid Flow, CMFF’03, Budapest, Hungary, Book of Proceedings 750-756, ISBN 963 420 777 4ö, ISBN 963 420 778 2 Balásházy I, Farkas Á, Szőke I and Hofmann W 2003e Simulation of airflow, aerosol deposition and clearance in central human airways. Particulate Matter and Health, 5th International Technion Symposium “Technology for Peace – Science for Mankind”.

99

Vienna, Austria, 2003. Proceedings, 57-65, ISBN 3-9501023-2-9. Editor: H. Gutmann. Published: Austrian Technion Society Balásházy I, Farkas Á, Szőke I and Hofmann W 2004a Computational fluid dynamics simulations of radioaerosol deposition and related health effects in central human airways J. Aerosol Sci. 2 S1205-6 Balásházy I, Farkas Á and Szõke I 2004b Simulation of deposition and alpha-hit distributions of inhaled radon progenies in bronchial airways. Conference of the Aerosol Society of the United Kingdom. Lung Modelling: Numerical and Experimental. University of Sheffield, United Kingdom. Book of Abstract Balásházy I, Szőke I and Farkas Á 2004c Numerical modelling of alpha-hit probability distributions in airways following the inhalation of radon progenies. General Meeting of the Roland Eötvös Physical Society, Szombathely, Hungary. Book of Abstracts 26 Balescu R 1975 Equilibrium and nonequilibrium statistical mechanics, John Wiley & Sons, Toronto Barabási A L, Buldyrev S V, Stanley H E and Suki B 1996 Avalanches in the lung: A statistical mechanical model Phys. Rev. Letters 76 2192-5 BEIR VI Report 1999 Health effects of exposure to radon, National Academy Press, Washington, D.C. Bejan A and Lorente S 2001 Thermodynamic optimization of flow geometry in mechanical and civil engineering J. Non-Equilib. Thermodyn. 26 305-54 Benett S H, Eldridge M W, Puente C E, Riedi R H, Nelson T R, Goetzman B W, Milstein J M, Singhal S S, Horsfield K and Woldenberg M J 2003 Origin of fractal branchig complexity in the lung, http://neonatology.ucdavis.edu/fractal Benz W 1990 Smoothed particle hydrodynamics: A review, Buchler Bernardi C and Maday Y 1996 Spectral element methods in Handbook of numerical analysis, North-Holland, Amsterdam Blot W J, Xu Z Y, Boice J D Jr, Zhao D Z, Stone B J, Sun J, Jing L B, Fraumeni J F Jr 1990 Indoor radon and lung cancer in China J. Natl. Cancer Inst. 7 1722-3

100

Bohm R, Nikodemova D and Holy K 2003 Use of various microdosimetric models for the prediction of radon induced damage in human lungs Radiat. Prot. Dosim. 104 127-37 Braga E J and Lemos M J S 2004 Numerical simulation of turbulent flow in small-angle diffusers and contractions using a new wall treatment and a linear high Reynolds k–ε model Numer. Heat Transfer 45 911-33 Bredberg J, Peng S H and Davidson L 2002 An improved k-ω turbulence model applied to recirculating flows Int. J. Heat and Fluid Flow 23 731-43 Broday D M 2004 Deposition of ultrafine particles at carinal ridges of the upper bronchial airways Aerosol Sci. Tech. 38 991-1000 Brooks A N and Hughes T J R 1982 Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with special emphasis on the incompressible NavierStokes equations Comput. Meth. Appl. Mech. Eng. 32 199-259 Canuto C, Hussaini M Y, Quarteroni A and Zhang T A 1988 Spectral methods in fluid dynamics, Springer, New-York Caro C G, Schroter R C, Watkins N and Sherwin S J 2001 Steady inspiratory flows in planar and non-planar models of human bronchial airways J. Physiol. 531 3P Cassee F R, Freijer J I, Subramaniam R, Asgharian B, Miller F J, Bree L and Rombout P J A 1999 Development of a model for human and rat airway particle deposition: implications for risk assessment. RIVM, Research for Man and Environment Report 650010018 Chadwick K H and Leenhouts H P 2002 What can we say about the dose-effect relationship at very low doses? J. Radiol. Prot. 22 A155-8 Chandrasekhar S 1943 Stochastic problems in physics and astronomy Rev. Modern Phys.

15 1-89 Chen H and Patel V 1998 Near wall turbulence models for complex flows including separation AIAA Journal 26 641-8

101

Chen S, Wang Z, Shan X, and Doolen G D 1992 Lattice Boltzmann computational fluid dynamics in three dimensions J. Stat. Phys. 68 379-400 Cheng Y S, Thou Y and Chen B T 1999 Particle deposition in a cast of human oral airways Aerosol Sci. Technol. 31 286-300 Churg A and Vedal S 1996 Carinal and tubular airway particle concentrations in the large airways of non-smokers in the general population: Evidence for high particle concentration at airway carinas Occup. Env. Med. 53 553-8 Clarke R 1999 Control of low level radiation exposure: time for a change? J. Radiol. Prot. 19 107-15 Clarke R 2001 Control of low level radiation exposure: what is the problem and how can it be solved? Health Physics 80 391-6 Cohen B L 2000 Validity of the linear no-threshold theory of radiation carcinogenesis at low doses Energy Environ. 11 149-66 Collatz L 1960 The numerical treatment of differential equations, Springer, Berlin Comer J K, Kleinstreuer C and Zhang Z 2001a Flow structures and particle deposition patterns in double-bifurcation airway models. Part 1. Air flow fields J. Fluid Mech.

435 25-54 Comer J K, Kleinstreuer C and Zhang Z 2001b Flow structures and particle deposition patterns in double-bifurcation airway models. Part 2: Aerosol transport and deposition J. Fluid Mech. 435 55-80 Cooper J A, Jackson P O, Langford J C, Petersen M R and Stuart B O 1973 Characteristics of attached radon-222 daughters under both laboratory and field conditions with particular emphasis upon underground mine environments. Report to the U.S. Bureau of Mines under contract H0220029, Richland, Crawford-Brown D J and Hofmann W 2002 Analysis of radon-induced lung cancer risk by a stochastic state-vector model of radiation carcinogenesis J. Radiol. Prot. 22 A61-5

102

Crowe C, Sommerfeld M and Tsuji Y 1998 Multiphase flows with droplets and particles, CRC Press, Boca Raton Davidson L 2003 An introduction to turbulence models Chalmers Publication 97/2: http://www.tfd.chalmers.se/~lada Dekker E 1961 Transition between laminar and turbulent flow in human trachea J. Appl. Physiol. 16 1060-4 Demidovich B P and Maron I A 1981 Computational mathematics, Mir Publishers, Moscow Dringó L 1992 Numerikus analízis I., Tankönyvkiadó, Budapest Duport P 2003 Low-dose radiation carcinogenesis in mammals: A review of a database of cancer induction by low-dose radiation in mammals: Initial observations and comparison with some human data. International symposium on biological effects of low

dose

radiation,

2002,

Japan,

Book

of

Proceedings.

http://www.ie.uottawa.ca/English/Reports/Lowdose%20carcinogenesis%20in%20ani mals-Rokkasho%202002.pdf Einstein A 1905 Molecular-kinetic theoretic aspects of the motion caused by heat of particles suspended in motionless fluids Ann. D. Physik 17 549-60 Farkas Á Hegedűs Cs and Pálfalvi J 2001a Realistic 3D flow calculations in the human tracheobronchial tree J. Aerosol Sci. 32, S803-4 Farkas Á, Balásházy I, Hegedűs Cs, Varga L and Barta L 2001b Numerical modelling of local deposition patterns of inhaled aerosol drugs in the tracheobronchial tree J. Aerosol Med. 14, 389 Farkas Á, Balásházy I, Hofmann W, Varga L and Pálfalvi J 2001c Flow structure computations in airways for the modeling of radionuclide inhalation. XXVI. Hungarian Health Physics Conference, Balatonkenese, Book of Abstracts 32 and 34 Farkas Á, Hegedűs Cs, Pálfalvi J, Kristóf G, Lohász M and Szabó J 2001d CFD simulations of air-flow in human airways for the modeling of inhalation of radionuclides. IRPA

103

Regional Congress on Radiation Protection in Cenral Europe, Radiation Protection in Health, Dubrovnik, Croatia. Proceedings, 3p-06, 1-6 Farkas Á, Balásházy I and Szőke I 2002a Numerical modeling of airflows and deposition patterns of radio-aerosol in central human airways J. Aerosol Sci. 34 651 Farkas Á, Balásházy I, Szőke I, Hofmann W and Golser R 2002b Simulation of deposition and activity distributions of radionuclides in human airways. European IRPA Congress, Towards Harmonization of Radiation Protection in Europe, Florence, Italy. Proceedings, 023-K, 1-9, ISBN 88-886-48-09-7 Farkas Á, Balásházy I and Szőke I 2002c Modeling the cellular distribution of inhaled dose in bronchial airways following radon exposure. Annual Hungarian Radiochemical Conference, Gyula, Hungary. Book of Abstract 61 Farkas Á, Balásházy I and Szőke I 2003a Simulation of flow fields and particle deposition patterns in human airways by computational fluid dynamics methods. 12th International Conference on Fluid Flow Technologies, Conference on Modelling Fluid Flow, CMFF’03, Budapest, Hungary, Book of Proceeding 743-749, ISBN 963 420 777 4ö, ISBN 963 420 778 2 Farkas Á, Balásházy I and Szőke I 2003b Numerical modelling of local deposition patterns, activity distributions and cellular hit probabilities of inhaled radon progenies in human airways. IRPA Regional Congress on Radiation Protection in Neighbouring Countries of Central Europe, Bratislava, Slovakia. Proceedings: VI.P1., p1.-4., ISBN 80-88806-43-7 Farkas Á, Szőke I and Balásházy I 2003c Simulation of activity and DNA alpha hit probability distributions deposited radon progenies in the epithelium of central human airways. EU-US Workshop on Molecular signatures of DNA damage induced stress responses. Italy, Book of Abstract 84-5 Farkas Á, Balásházy I and Szőke I 2003d Numerical modelling of deposition and activity distribution of inhaled radon progenies. Annual Radiochemical Conference Balatonföldvár, Book of Abstract 22 and 46 Farkas Á, Balásházy I and Szőke I 2004a Numerical modelling of airflow and aerosol deposition in a human alveolus J. Aerosol Sci. 2 S1113-4 104

Farkas Á 2004b Alveoláris levegőáramok és részecskekiülepedések numerikus modellezése FLUENT CFD kóddal, nem publikált dolgozat Farkas Á and Balásházy I 2004c CFD simulation of airflow and particle deposition patterns in diseased human airways. Conference of the Aerosol Society of the United Kingdom. Lung Modelling: Numerical and Experimental. University of Sheffield, United Kingdom Farkas Á, Balásházy I and Szőke I 2004d CFD simulation of activity distributions of deposited radon progenies in central human airways. European Radiation Research Conference, 33rd Annual Meeting of the European Society for Radiation Biology. Book of Abstracts 100 Farkas Á and Balásházy I 2005a Simulation of airflow field and particle deposition in the diseased central human airways J. Appl. Physiol (under preparation) Farkas Á, Balásházy I and Szőcs K 2005b Characterization of regional and local deposition of inhaled aerosol drugs in the respiratory system by computational fluid and particle dynamics methods J. Aerosol Med. (invited paper, in press) Farkas Á and Balásházy I 2005c Computational modelling of pulmonary airflow field and particle deposition Phys. Med. Biol. (submitted) Field R W, Lynch C F, Steck D J, Smith B J, Bruce C P, Neuberger J S, Woolson R F, Fisher E F, Platz C E and Robinson R A 2000 Residential gas exposure and lung cancer: The Iowa radon lung cancer study Am. J. Epidemiol. 151 1091-102 Finlay W H, Stapleton K W and Yokota J 1996 On the use of computational fluid dynamics for simulating flow and particle deposition in the human respiratory tract J. Aerosol Med. 9 329-41 FLUENT Manual 2001, Sheffield GAMBIT Manual 2001, Sheffield Goldman M, Catlin R J and Anspaugh L 1987 Health and Environment Impact of the Chernobyl Accident US Department of Energy research report, DOE/RR-0232

105

Gonzalez A J 1998 Radioactive residues of the cold war periode: A radiological legacy IAEA Bulletin 40/4/1998 Gradon L and Orlicki D 1990 Deposition of inhaled aerosol particles in a generation of the tracheobronchial tree J. Aerosol Sci. 21 3-19 Gradon L and Podgorski A 1996 Deposition of inhaled particles: Discussion of present modeling techniques J. Aerosol Med. 9 343-55 Grzybowski K and Gradon L Resuspension of particles from the powder structure J. Aerosol Sci. 1 S617-8 Gupta D and Peters M 1985 A Brownian dynamics simulation of aerosol deposition onto spherical collectors J. Colloid Interface Sci. 104 375-89 Gurman J L, Lippmann N and Schlesinger L B 1984 Particle deposition in replicate casts of the human upper tracheobronchial tree under constant and cyclic inspiratory flow. I. Experimental Aerosol Sci. Technol. 3 245-52 Haider A and Levenspiel O 1989 Drag coefficient and terminal velocity of spherical and nonspherical particles Pow. Tech. 58 63-70 Hall E J 2003 The bystander effect Health Physics 85 31-5 Hammersley J R and Olson D E 1992 Physical models of the smaller pulmonary airways J. Appl. Physiol. 72 2402-14 Haninger T. 1997 Size distributions of radon progeny and their influence on lung dose. Radon and Thoron in the Human Environment, Proceedings of the 7th Tohwa University International Symposium, Editors: Katase A. and Shima M., Publisher: World Scientific Singapore, New Jersey, London, Hong Kong Hausermann S, Bailey A G and Maul C 2000 A system to reproduce human breathing patterns: Its development and validation J. Aerosol Med. 13 199-204 Hegedűs Cs J, Balásházy I and Farkas Á 2004. Detailed mathematical description of the geometry of airway bifurcations Resp. Physiol. Neurobiology 141 99-114 Heistracher T and Hofmann W 1995 Physiologically realistic models of bronchial airway bifurcations J. Aerosol Sci. 26 497-509 106

Heistracher T 1996 Numerical simulation of airflow and particle deposition in bronchial airway bifurcation models, PhD Thesis, University of Salzburg, Salzburg, Austria Hirsch C 1988 Numerical computation of internal and external flows Vol. 1: Fundamentals of Numerical Discretization, John Wiley & Sons, Chichester Horsfield K and Cumming G 1967 Angles of branching and diameters of branches in the human bronchial tree Bul. Math. Biophysics 29 245-59 Horsfield K, Dart G, Olson D E, Filley G F and Cumming G 1971 Models of the human bronchial tree J. Appl. Physiol. 31 207-17 Horváth L 1998 Funkcionális anatómia, Nemzeti Tankönyvkiadó, Budapest Hutchinson B R and Raithby G D 1986 A multigrid method based on the additive correction strategy Numer. Heat Transfer 9 511-37 Hyatt R E and Wilcox R E 1961 Extrathoracic airway resistance in man J. Appl. Physiol.

16 326-30 Hyman J, Knapp R and Scovel J C 1992 High order finite volume approximations of differential operators on nonuniform grids Physica D 60 112-38 IAEA/WHO/UNDP report 2005 Chernobyl: The true scale of the accident, Wienna ICRP Publication 65 1993 Protection against Radon-222 at home and work. Annals of the ICRP 23 (2) ICRP Publication 66 1994 Human respiratory tract model for radiological protection, Annals of the ICRP 24 Isabey D and Chang HK 1982 A model study of flow dynamics in human central airways. Part II: Secondary flow velocities Resp. Physiol. 49 97-113 James A C 1988 Lung dosimetry. In: Radon and its decay products in indoor air, Wiley and Sons Inc., New York Jaworowski Z 1999 Radiation risk and ethics Phys. Today 52 24-9 Jiang Y and Przewkas A J 1994 Implicit, pressure-based incompressible Navier-Stokes equations solver for unstructured meshes AIAA-94-0305

107

Johnston J R, Isles K D and Muir D C F 1977 Inertial deposition of particles in human branching airways In: Inhaled Particles IV, Pergamon press, Oxford Kac M and Logan J 1987 Fluctuations in Fluctuation Phenomena, Elsevier Publishers, Amsterdam Kader B 1993 Temperature and concentration profiles in fully turbulent boundary layers Int. J. Heat Mass Transfer 24 1541-4 Kassai D 1950 A tüdő szegmentumai, Akadémiai Kiadó, Budapest Katz I M and Martonen T B 1999 A numerical study of particle motion within human larynx and trachea J. Aerosol Sci. 30 173-83 Kesavan P C 1996 in High levels of natural radiation, Wei L, Sugahara T and Tao Z eds. Elsevier, Amsterdam Kim C S and Garcia L 1991 Particle deposition in cyclic bifurcating tube flow Aerosol Sci. Technol. 14 302-15 Kitaoka H and Suki B 1997 Branching design of the bronchial tree based on diameterflow relationship J. Appl. Physiol. 82 968-76 Kitaoka H and Takari R 1999 Fractal analysis of the human fetal lung development, Forma 14 205-12 Klebanov L, Rachev S and Yakovlev A 1993 A stochastic model of radiation carcinogenesis: Latent distributions and their properties Math. Biosci. 113 51-75 Kleinstreuer C and Zhang Z 2003. Laminar-to-turbulent fluid-particle flows in a human airway model. Int. J. Multiphase Flow. 29 271-89 Kim C S and Iglesias A J 1989 Deposition of inhaled particles in bifurcating airway models J. Aerosol Med. 2 1-14 Kim C S and Garcia L 1991 Particle deposition in cyclic bifurcating tube flow Aerosol Sci. Technol. 14 302-15 Kim C S and Fisher D M 1999 Deposition characteristics of aerosol particles in sequentially bifurcating airway models Aerosol Sci. Technol. 31 198-220

108

Kinsara A A, Loyalka S K, Tompson R V, Miller W H and Holub R F 1995 Deposition patterns of molecular phase radon progeny 218Po in Lung Bifurcations Health Phys.

68 371-81 Koblinger L and Hofmann W 1985 Analysis of human morphometric data for stochastic aerosol deposition calculations Phys. Med. Biol. 30 541-56 Koblinger L and Hofmann W 1990 Monte Carlo modeling of aerosol deposition in human lungs. Part I: Simulation of particle transport in a stochastic lung structure J. Aerosol Sci. 21 661-74 Kotin P and Falk H L 1959 The role and activation of environmental agents in the pathogenesis of lung cancer. I. Air pollution Cancer 12 147-63 Kreienbrock L, Kreuzer M, Gerken M, Dingerkus G, Wellmann J, Keller G and Wichmann H E 2001 Case-control study on lung cancer and residential radon in western Germany Am. J. Epidemiol. 153 42-52 Lajos T 2000 Az áramlástan alapjai, Műegyetemi Kiadó, Budapest Landau L and Lifchitz E 1971 Mecanique des fluides, Editions MIR, Moscou Lee W J, Kawahashi M and Hirahara H 2005 Experimental investigation of oscillatory flow in a micro-channel of bronchiole model. 6th World Conference on Experimental Heat Transfer, Fluid Mechanics and Thermodynamics, Matsushima, Miyagi, Japan Li A and Ahmadi G 1992 Dispersion and deposition of spherical particles from point sources in a turbulent channel flow Aerosol Sci. Tech. 16 209-26 Lide D R 1998-1999 Handbook of Chemistry and Physics, 79th ed., CRC Press Liu C and Ahmadi G 2005 Transport and deposition of particles near a building model Build. Environ. (in press) Liu Y, So R M C and Zhang C H 2002 Modeling the bifurcating flow in a human lung airway J. Biomechanics 35 465-73 Macklin C C 1956 Induction of bronchial cancer by local massing of carcinogen in outdrifting mucus J. Thorac. Surg. 31 238-44

109

Magyar P, Hutás I és Vastag E 1998 Pulmonológia, Medicina Kiadó, Budapest Martin D and Jacobi W 1972 Diffusion deposition of small-sized particles in the bronchial tree Health Phys. 23 23-9 Mathur S R and Murthy J Y 1997 A pressure-based method for unstructured meshes Numer. Heat Transfer B 31 195-215 Matida E A, Finlay W H, Lange C F and Grgic B 2004 Improved numerical simulation of aerosol deposition in an idealized mouth-throat J. Aerosol Sci. 35 1-19 Mauroy B, Filoche M, Adrade Jr. J S and Sapoval B 2003 Interplay between geometry and flow distribution in an airway tree Phys. Rev. Letters 90, 1481-61 Mercer R M, Michael L R and James D C 1991 Radon dosimetry based on the depth distribution of nuclei in human and rat lungs Health Phys. 61 117-30 Miguel A F, Reis A H and Aydin M 2004 Aerosol particle deposition and distribution in bifurcating ventilation ducts J. Hazard. Mater. 116 249-55 Mill A J, Frankenberg D, Bettega D, Hieber L, Saran A, Allen L A, Calzolari P, Frankenberg-Schwager M, Lehane M M, Morgan G R, Pariset L, Pazzaglia S, Roberts C J and Tallone L 1998 Transformation of C3H10T1/2 cells by low doses of ionizing radiation: a collaborative study of six European laboratories strongly supporting the linear dose-response relationship J. Radiol. Prot. 18 79-100 Monaghan J J 1992 Smoothed particle hydrodynamics Ann. Rev. Astron. Astrophys. 30 543-74 Moody D. and Lozanoff S 1997 A practical computer program for generating threedimensional models of anatomical structures Proc. of 14th Annual Meeting of the American

Association

of

Clinical

Anatomists,

Honolulu,

Hawai’i

http://www.surfdriver.com/docs/reconstruction.html Morsi S A and Alexander A J 1972 An investigation of particle trajectories in two-phase flow systems J. Fluid Mech. 55 193-208 Moskal A and Gradon L 2002 Temporary and spatial deposition of aerosol particles in the upper human airways during a breathing cycle J. Aerosol Sci. 33 1525-39

110

Nikezic D, Haque A K M M and Yu K N 2002 Absorbed dose delivered by alpha particles calculated in cylindrical geometry J. Environ. Radioactivity 60 293-305 NRC 1991 Comparative dosimetry of radon in homes and mines, National Academy Press, Washington, D.C. Oldham M J, Phalen R F and Heistracher T 2000 Computational fluid dynamic predictions and experimental results for particle deposition in an airway model Aerosol Sci. Technol. 32 61-71 Ounis H and Ahmadi G 1990 Analysis of dispersion of small spherical particles in a random velocity field J. Fluids Engineering 112 114-20 Ounis H, Ahmadi G and McLaughlin J B 1991 Brownian diffusion of submicrometer particles in the viscous sublayer J. Colloid Interface. Sci. 143 266-77 Paek D and McCool D 1992 Breathing patterns during varied activities J. Appl. Physiol.

73 887-93 Pálfalvi J, Dám A M, Bogdándi E M, Polónyi I, Szabó J, Balásházy I and Farkas Á 2003 Etched track detectors and the low dose problem. Rad. Prot. Dos. 103 229-34 Patankar S V 1980 Numerical heat transfer and fluid flow, McGraw-Hill, New-York Pedley T J, Schroter R C and Sudlow 1971 Flow and pressure drop in systems of repeatedly branching tubes J. Fluid Mech. 46 365-83 Pedley T J 1977 Pulmonary fluid dynamics, Ann. Rev. Fluid Mech. 9, 229-74 Peng S H, Davidson L and Holmberg S 1997 A modified low Reynolds number k-ω model for recirculating flows ASME J. Fluid Eng. 119 867-75 Peric M 1985 A finite volume method for the prediction of three-dimensional fluid flow in complex ducts, Ph.D. thesis, University of London, UK Perzl M A, Schulz H, Paretzke H G, Englmeier K H and Heyder J 1996 Reconstruction of the lung geometry for the simulation of aerosol transport J. Aerosol Med. 9 409-18 Peyret R 2002 Spectral method for incompressible viscous flows, Springer, New-York

111

Phalen R F, Oldham M J, Beaucage C B, Crocker T T and Mortensen J D 1985 Postnatal enlargement of human tracheobronchial airways. Anat. Rec. 212 368-80 Philips C G and Kaye S R 1997 On the assymetry of the bifurcations in the bronchial tree Resp. Physiol. 107 85-98 Pierce D A, Shimizu Y, Preston D L, Vaeth M and Mabuchi K 1996 Studies of the mortality of atomic bomb survivers Report 12, Part 1 Canc. Rad. Res. 146 1-27 Pirchan A and Sikl H 1932. Cancer of the lung in miners of Jachymov (Joachimsthal) Report of cases observed in 1929-1930 Am. J. Cancer 16 681-722 Press W H, Vetterling W T, Teukolsky S A and Flannery B P 2002 Numerical recipes in C++: the art of scientific computing, William H. Press, 2nd edition Preston D L, Kusumi S, Tomonaga M, Izumi S, Ron E, Kuramoto A, Kamada N, Dohy H, Matsui T, Nonaka H, Thompson D H, Soda M and Mabuchi K 1994 Cancer incidence in atomic bomb survivers, part III: leukemia, lymphoma and multiple myeloma, 1950-1987 Radiat. Res. 137 68-97 Raabe O G, Yeh H C, Schum G M and Phalen R F 1976 Tracheobronchial Geometry: Human, Dog, Rat, Hamster, LF-53. Lovelace Foundation Report, Albuquerque, New Mexico, USA Raick C, Ramuzat A, Corieri P and Riethmuller M L 2003 Numerical and experimental study of spatial periodicity of steady air flows in pulmonary bifurcations Proc. FLUCOME Ralston A 1969 Bevezetés a numerikus analízisbe, Műszaki Könyvkiadó, Budapest Ramuzat A, Day S and Riethmuller M L 1998 Steady and unsteady LDV and PIV investigations of flow within a 2D lung bifurcations model Proc. ISALTFM Reineking A, Becker K H, Porstendörfer J 1988 Measurements of activity size distributions of the short-lived radon daughters in the indoor and outdoor environment. Rad. Prot. Dosim. 24 245-50 Reineking A, Knutson E O, George A C, Solomon S B, Kesten J, Butterweck G and Pörstendörfer J 1994 Size distribution of unattached and aerosol-attached short-lived

112

radon decay products: some results of intercomparison measurements Rad. Prot. Dosim. 56 113-8 Reis A H, Miguel A F and Aydin M 2004 Constructal theory of flow architecture of the lungs Med. Phys. 31 1135-40 Rhie C M and Chow W L 1983 Numerical study of the turbulent flow past an airfoil with trailing edge separation AIAA Journal 21 1525-32 Roache P J 1972 Computational fluid dynamics, Hermosa, Albuquerque, New Mexico Rubin B K 2002 Physiology of airway mucus clearance Respir. Care 47 761-8 Ruosteenoja E, Makelainen I, Rytomaa T, Hakulinen T and Hakama M 1996 Radon and lung cancer in Finland Health Phys. 71 185-9 Saffman P G 1965 The lift on a small sphere in a slow shear flow J. Fluid Mech. 22 385400 Samet J M, pathak D R, Morgan M V, Marbury M C, Key C R and Valdivia A A 1989 Radon progeny exposure and lung cancer risk in New Mexico U miners: a casecontrol study Health Phys. 56 415-21 Sauret V, Halson P M, Brown I W, Fleming J S and Bailey A G 2002 Study of the threedimensional geometry of the central conducting airways in man using computed tomographic (CT) images J. Anatomy 200 123-34 Schlesinger R B and Lippmann M 1978 Selective particle deposition and bronchogenic carcinoma Environ. Res. 15 424-31 Schlesinger M F and West B J 1991 Complex fractal dimensions of the bronchial tree, Phys. Rev. Letters 67 2106-8 Schoenberg J B, Klotz J B, Wilcox H B, Nicholls G P, Gil M T, Stemhagen A and Mason T J 1990 A case-control study of residential radon and lung cancer among New Jersey women. Cancer Res. 50 6520-24 Schöllnberger H, Mitchell R E J, Azzam E I, Crawford-Brown D J and Hofmann W 2002 Explanation of protective effects of low doses of γ-radiation with a mechanistic radiobiological model Int. J. Radiat. Biol. 78 1159-73 113

Schroter R C and Sudlow M F 1969 Flow patterns in models of the human bronchial airways Respir. Physiol. 7 341-55 Sinclair W K 1996 Linear-nonthreshold model vindicated by new data Nucl. Week 6 1-3 Sohrabi M 1990 in High levels of natural radiation, Sohrabi M, Durrani S A eds. International Atomic Energy Authority, Vienna, Austria Stapleton K W, Guentsch E, Hoskinson M K and Finlay W H 2000 On the suitability of k-ε turbulence modeling for aerosol deposition in the mouth and throat: A comparison with experiment J. Aerosol Sci. 31 739-49 Sujeer M K, Buldyrev S V, Zapperi S, Andrade J S, Stanley H E and Suki B 1997 Volume distributions of avalanches in lung inflation: A statistical mechanical approach Phys. Rev. E 56 3385-94 Suki B 2002 Fluctuations and Power Laws in Pulmonary Physiology Am. J. Resp.Crit.Care Med. 166 133-7 Szőke I, Balásházy I and Farkas Á 2002a Numerical modeling of mucociliary clearance of inhaled radionuclides. 27th Annual Meeting on Radiation protection. Mátrafüred, Hungary. Book of Abstract, 33 and 37 Szőke I, Balásházy I, Farkas Á, Hofmann W and Sturm R 2002b Numerical modeling of mucociliary clearance of inhaled particles. Extended abstracts of the Sixth International Aerosol Conference, 1139-40, Editor: Chiu-Sen Wang, ISBN 986-80544-1-9 Szőke I, Balásházy I and Farkas Á 2002c Effect of mucociliary clearance on the distribution of bronchial aerosol deposition J. Aerosol Sci. 34 664-5 Szőke I, Balásházy I, Farkas Á, Patonay L, Hrabák K and Kerényi T 2002d Deposition of inhaled radionuclides in computer tomographically reconstructed human airways. European IRPA Congress, Towards Harmonization of Radiation Protection in Europe, Florence, Italy. Proceedings, 119-K, 1-7, ISBN 88-886-48-09-7 Szőke I, Balásházy I and Farkas Á 2002e Numerical simulation of clearance in tracheobronchial airway bifurcations. Annual Hungarian Radiochemical Conference, Gyula, Hungary. Book of Abstract 63

114

Szőke I, Balásházy I, Sárkány Z and Farkas Á 2002f Modeling the deposition of inhaled radioaerosols in airway geometries reconstructed from CT images. Annual Hungarian Radiochemical Conference, Gyula, Hungary. Book of Abstract 62 Szőke I, Balásházy I, Farkas Á and Patonay L 2003a Deposition and mucociliary clearance patterns in airway geometries reconstructed by medical imaging techniques J. Aerosol Sci. S417-8. ISBN 0021-8502 Szőke I, Balásházy I and Farkas Á 2003b Alpha hit probability distributions of deposited radon progenies in the cell nuclei of the central airway epithelium. Annual Radiochemical Conference, Balatonföldvár. Book of Abstract 20 and 45 Szőke I, Balásházy I, Szabó J, Karlinger K, Patonay L, Petneházy Ö and Kerényi T 2004a Human airway models constructed by medical imaging techniques J. Aerosol Sci. S1135-36 Szőke I, Balásházy I and Farkas Á 2004b Alpha-hit probability distributions of deposited radon progenies in cell nuclei and cell surroundings of the central airway epithelium J. Aerosol Sci. 2 S1133-4 Szőke I, Balásházy I and Farkas Á 2004c Alpha-hit probability distributions of deposited radon progenies in cell nuclei, cells and cell surroundings of the central airway epithelium. European Radiation Research Conference, 33rd Annual Meeting of the European Society for Radiation Biology. Book of Abstracts 274 Talbot L, Cheng R K, Schefer R W and Willis D R 1980 Thermophoresis of particles in a heated boundary layer J. Fluid Mech. 101 737-58 Thompson D E, Mabuchi K, Ron E, Soda M, Tokunaga M, Ochicubo S, Sugimoto S, Ikeda T, Terasaki M, Izumi S and Preston D L 1994 Cancer incidence in atomic bomb survivers, part II: solid tumors, 1958-1987 Radiat. Res. 137 17-67 Tobin M J, Chadha T S, Jenouri G, Birch S J, Gazeroglu H B and Sackner M A Breathing patterns. 1. Normal subjects Chest 84 202-5 Tóth Á 1983 A lakosság természetes sugárterhelése, Akadémiai Kiadó, Budapest

115

Trosko J E 1996 Role of low level ionising radiation in multi-step carcinogenic process Health Phys. 70 812-22 Tsuda A, Rogers R A, Hydon P E and Butler J P 2002 Chaotic mixing deep in the lung Proc. Natl. Acad. Sci. 99 10173–8 Uhlenbeck G E and Ornstein L S 1930 On the theory of the Brownian motion Phys. Rev.

36 823-41 UNSCEAR 1994 Report Sources and effects of ionizing radiation Annex A. Epidemiological studies of radiation carcinogenesis, New York UNSCEAR 2000 Sources and effects of ionizing radiation Report to the General Assembly, with scientific annexes, New-York Veeze P 1968 Rationale and methods of early detection in lung cancer. Assen, The Netherlands, Van Gorcum Venkatakrishnan V 1993 On the accuracy of limiters and convergence to steady state solutions AIAA 93-0880 Versteg H and Malalasekera W 1995 An introduction to computational fluid dynamics: The finite volume method, Addison-Wesley Vinokur M 1989 An analysis of finite-difference and finite-volume formulations of conservation laws J. Comput. Physics 81 1-52 Wechsatol W, Lorente S and Bejan A 2004 Tree-shaped flow structures: are both thermal-resistance and flow-resistance minimisations necessary? Int. J. Energy 1 2-17 Weibel E R 1963 Morphometry of the human lung. Springer, Academic Press, Berlin, New-York Wikipedia 2005 http://en.wikipedia.org/wiki/Wiener_process Wilcox D C 1993 Turbulence modeling for CFD, DCW Industries Inc., LA Canada, CA Yeh H C and Schum G M 1980. Models of human lung airways and their application to inhaled particle deposition Bull. Math. Biol. 42 461-80

116

Zhang Z, Kleinstreuer C and Kim C S 2002 Cyclic micron-size particle inhalation and deposition in a triple bifurcation lung airway model J. Aerosol Sci. 33 257-81 Zhang Z and Kleinstreuer C 2003 Modeling of low Reynolds number turbulent flows in locally constricted conduits: A comparison study AIAA Journal 41 831-40 Zhang Z, Kleinstreuer C, Donohue D F and Kim C S 2005 Comparison of micro- and nano-size particle deposition in a human upper airway model. J. Aerosol Sci. 36 21133 Zhao Y and Lieber B B 1994 Steady inspiratory flow in a model symmetric bifurcation ASME J. Biomech. Eng. 116 488-96

117

15. Köszönetnyilvánítás

Köszönettel és hálával tartozom mindenekelőtt témavezetőmnek, Dr. Balásházy Imrének az elmúlt években nyújtott segítségéért, hasznos megjegyzéseiért, tanácsaiért és bátorításáért. Köszönet Dr. Gadó Jánosnak, a Magyar Tudományos Akadémia KFKI Atomenergia Kutatóintézet igazgatójának, aki biztosította számomra a kutatómunkához szükséges tárgyi és anyagi feltételeket. Hasonló köszönettel tartozom Dr. Török Szabinának, a Sugárvédelmi és Környezetfizikai Laboratórium vezetőjének és Dr. Zombori Péternek, a Sugárvédelmi Csoport vezetőjének. Köszönet a Sugárvédelmi Laboratórium valamennyi munkatársának az évek során nyújtott segítségért. Hálával tartozom továbbá családomnak, a nyugodt körülményekért, az erkölcsi támogatásért és a sok türelemért.

Farkas Árpád

118

16. Függelék 1. A dolgozatban használt rövidítések listája BEIR

-

Biological Effects of Ionizing Radiation – project of the US National Academies (Egyesült Államok Nemzeti Akadémiáinak az Ionizáló Sugárzások Biológiai Hatásait Vizsgáló projektje)

CAD

-

Computer Aided Design (Számítógép-Vezérelt Tervezés)

CAE

-

Computer Aided Engineering (Számítógép-Vezérelt Technika)

CAM

-

Computer Aided Manufacturing (Számítógép-Vezérelt Gyártás)

CFD

-

Computational Fluid Dynamics (Numerikus Áramlástan)

CT

-

Computer Tomography (Komputer Tomográfia)

FKT

-

Fokozott Kiülepedési Tényező

FPAT

-

Fokozott Potenciális Alfa-energia Tényező

HPM

-

Hit Probability Model (Találati Valószínűség Modell)

LBM

-

Lattice Boltzmann Method (Rács-Boltzmann Módszer)

LDA

-

Laser Doppler Anemometer (Lézer Doppler Anemométer)

DNS

-

Dezoxiribo-Nuklein-Sav

IAEA

-

International Atomic Energy Agency (Nemzetközi Atomenergia Ügynökség)

ICRP

-

International Comission on Radiological Protection (Nemzetközi Sugárvédelmi Bizottság)

LNT

-

Linear-Nonthreshold (Lineáris Küszöb Nélküli)

MRB

-

Morphologically Realistic Bifurcation (Morfológiailag Realisztikus Elágazás)

MRI

-

Magnetic Resonance Imaging (Mágneses Magrezonanciás Képalkotás)

NCRP

-

National Council on Radiation Protection and Measurements (Egyesült Államok Sugárvédelmi Nemzeti Tanácsa)

119

NRC

-

Nuclear Regulatory Comission (Egyesült Államok Nukleáris Szabályozó Bizottsága)

PIV

-

Particle Image Velocimetry (Részecskekép Alapú Sebességmérés)

SEM

-

Spectral Element Method (Spektrális Elem Módszer)

SIMPLE

-

Semi-Implicit Methods for Pressure-Linked Equations (NyomásCsatolt Egyenletek Fél-Implicit Módszerei)

SPH

-

Smoothed Particle Hydrodynamics (Simított Részecske Hidrodinamika)

SVM

-

State Vector Model (Állapotvektor Modell)

UNDP

-

United Nations Development Programme (Egyesült Nemzetek Szervezetének Fejlesztési Programja)

UNSCEAR

-

United Nations Scientific Committee on the Effects of Atomic Radiation (Egyesült Nemzetek Szervezetének Atomi Sugárzások Hatásait Vizsgáló Tudományos Bizottsága)

WHO

-

World Health Organization (Egészségügyi Világszervezet)

XRF

-

X Ray Fluorescence (Rötntgen Fluoreszcencia)

120

2. A dolgozatban használt jelölések listája Latin szimbólumok a1, a2, a3

empirikus konstansok a súrlódási együttható képletében

ap, anb

a p cellára vonatkozó linearizált algebrai együtthatói

d

távolság

dp

részecskeátmérő

dr

helyzetvektor

es, et

egységvektorok

g

gravitációs gyorsulásvektor

gx

gravitációs gyorsulásvektor x tengelyre vett vetülete

k

turbulens kinetikus energia

kB

Boltzmann állandó

lp

linearizációs együttható reciproka

m

tömeg

p

nyomás

q

Chandrasekhar által bevezetett mennyiség, az egységnyi tömegre jutó termikus energia és a relaxációs idő aránya [Chandrasekhar 1943]

t

idő

u n

részecskesebesség x komponense n+1

u ,u

a levegősebesség x komponense az n illetve (n+1)-edik iteráció után

upn, upn+1

a részecskesebesség x komponense az n illetve (n+1)-edik iteráció után

u*

a levegősebesség x komponensének két egymásutáni iterációs lépésre vett számtani átlaga

up

részecskesebesség x komponense

ui, uj. uk

sebességkomponensek

v

sebességvektor

xi, xj. xk

koordináták

xp

részecske x koordinátája

121

A

területvektor

A, B

sztochasztikus gyorsulások

CC

Cunningham csúszási faktor

CD

súrlódási együttható

D

diffúziós együttható

Dda

légúti elágazás A oldali leányágának átmérője

Ddb

légúti elágazás B oldali leányágának átmérője

Dd

csőátmérő

Df

diffúziós tag a cella határán

Gk

turbulens kinetikus energia keletkezése

Gω

fajlagos disszipációs ráta keletkezése

Jf

tömegáram a cella határán

Lda’

légúti elágazás A oldali leányágának hossza

Ldb’

légúti elágazás B oldali leányágának hossza

Lp’

légúti elágazás anyaágának hossza

Q

térfogatáram

Rer

részecske Reynolds szám

Rp

légúti elágazás anyaágának sugara

S

feszültségtenzor

Sij

feszültségtenzor komponensek

SΦ

Φ mennyiség forrása

Sm

forrástag (tömeg)

Sp

p cellán linearizált egyenlet forrástagja

St

Stokes szám

T

abszolút hőmérséklet

V

térfogat

Yk

turbulens kinetikus energia disszipációja

Yω

fajlagos disszipációs ráta disszipációja

122

Görög szimbólumok α, α∞

fajlagos disszipációs ráta keletkezésének becslésénél használt állandók

α0

a divergenciatétel arányossági tényezője

α*, α*∞

turbulens viszkozitás becslésénél használt állandók

β

a relaxációs idő reciproka

ε

turbulens kinetikus energia disszipációja

ζ(0,1)

nulla várható értékű, egységnyi szórású Gauss eloszlás

η

kiülepedési hatásfok

λ

közepes szabad úthossz

µ

viszkozitás

µt

turbulens viszkozitás

ρ

sűrűség

ρp

részecskesűrűség

σk , σ ω

turbulens Prandtl számok

τ

feszültségtenzorok

τij

feszültségtenzor komponensek

τp

relaxációs idő

ϕa

légúti elágazás A oldali szaggitális szöge

ϕb

légúti elágazás B oldali szaggitális szöge

ω

turbulens kinetikus energia fajlagos disszipációs rátája

Γ

tetszőleges diffúziós együttható

Γk

turbulens kinetikus energia diffúziója

Γω

fajlagos disszipációs ráta diffúziója

Φ

tetszőleges skalármennyiség

Φa

légúti elágazás A oldali szaggitális szögének maximuma

Φb

légúti elágazás B oldali szaggitális szögének maximuma

Φf

a Φ mennyiség értéke a cella határán

Φszi

Φ értéke a szélirányú cellában

123