1. No category

239
9 Simulointitehtävän asettaminen
9.1 Simulointitehtävä
Simulointitehtävissä on tietty määrä samantyyppisiä vaiheita. Virtauslaskennan yhteydessä näissä vaiheissa on luonnollisesti jossain määrin eroja virtaustyypin mukaan. On hyödyksi identifioida simulointitehtävän vaiheet, vaikka ne yleensä menevätkin käytännön työssä sekaisin: ’Ideointivaihe’ ei suinkaan tule valmiiksi tehtävää aloitettaessa, vaan ensimmäiset laskelmat saattavat muuttaa tilannetta huomattavastikin ja koko suunnitelma menee uusiksi. On myös syytä pitää mielessä, mitä
tietokone oikein tekee: valtavan suuren määrän peruslaskentaoperaatioita ja lopputuloksen pitäisi olla tietyn tyyppisten osittaisdifferentiaaliyhtälöiden likimääräinen
ratkaisu. Simulointi ei sisällä mitään ’älyä’ sen enempää kuin tehtävän asettelija
on koneelle kertonut. Tämä seikka helposti unohtuu käytettäessä valmis-ohjelmia,
joista on kehittynyt tavaramerkkejä ja joista ainakin osa pyrkii markkinoinnissaan
painottamaan helppokäyttöisyyttä.
Virtauksen yhteydessä osittaisdifferentiaaliyhtälöitä yleensä joudutaan mallintamaan, joten tuloksessa on numeeristen approksimaatioiden lisäksi fysikaalisia approksimaatioita. Osittaisdifferentiaaliyhtälöiden ratkaisun asettamisessa on seuraavat vaiheet:
1. on määriteltävä tietty pistejoukko, jota käytetään approksimoitaessa numeerisesti osittaisdifferentiaaliyhtälöitä. Tämä pätee kaikkiin ratkaisumenetelmiin.
FLUENTissa ja useimmissa kaupallisissa koodeissa käytetään kontrollitilavuusmenetelmää, jossa pistejoukosta muodostetaan tasealueita (’kontrollitilavuuksia’).
2. on asetettava reunaehdot. Reunaehdot ovat osittain luonteeltaan fysikaalisia,
osittain numeerisia. Tämä aiheutuu siitä, että reunaehdot on asetettava fysikaalisesti oikein. Koska ehtoja tarvitaan numeerisista syistä kaikilla reunoilla yhtä paljon kuin ratkaistavien yhtälöiden määrä, osa ehdoista on luonteel-
9.1. SIMULOINTITEHTÄVÄ
240
taan numeerisia. Fysikaalisia ehtoja ei voida asettaa liikaa, koska tilanne tulisi
ylispesifioiduksi. Tasapainotilassa ei siis voida esimerkiksi asettaa toisistaan
riippumattomia arvoja sisään- ja ulostulomassavirralle.
3. on asetettava virtaavan aineen ominaisuudet kuten tiheys, viskositeetti jne.
Monimutkaisissa tilanteissa virtaukseen liittyy fysiikkaa, joka on mallinnettava.
4. lopuksi on asetettava simulointia säätelevät parametrit. Eri tilanteet vaativat
konvergoituakseen erilaisia parametriasetuksia ja viime kädessä vain laskemalla on mahdollista saada tuntuma sopivista parametrien arvoista. Manuaaleissa voidaan antaa vain suuruusluokka-arvioita ja suosituksia, joista lähdetään liikkeelle.
Kolmannen vaiheen jälkeen on saatu koottua epälineaarinen yhtälöryhmä, joka
on ratkaistava. Ratkaisuun liittyvät parametrit annetaan neljännessä vaiheessa. Ratkaisun tuloksena saadaan osittaisdifferentiaaliyhtälöille approksimatiivinen ratkaisu
määritellyissä pisteissä. Virtaussimulointi on siis oikeastaan osittaisdifferentiaaliyhtälöiden ratkaisemista ja koostuu yllä olevista tehtävistä. Tehtäväasettelua kuvattiin
jo aiemmin ensimmäisessä luvussa, mutta otetaan työvaiheet vielä lyhyesti esille
(Shaw’n mukaan):
• Alussa tehtävä ideointi. On ehdottoman tärkeää ensin miettiä tilannetta paperin ja kynän kanssa. Käytännön laskentatehtävät ovat yleensä työläitä ja pitkiä. Koska tehtävät ovat yleensä vaativia ei kannata syöksyä koneen pariin,
vaan aloittaa kartoittamalla tilanteeseen liittyvä fysiikka, geometrian vaatimukset jne. Tässä vaiheessa hahmotellaan geometriaa, yksinkertaistetaan sitä, pohditaan alustavasti reunaehtoja jne. Ja kuten aikaisemmin on todettu,
lasketaan tilanteeseen liittyvät tärkeimmät dimensiottomat luvut.
• Hilangenerointi on yleensä vielä nykyisin aikaa vievin vaihe. Tämänkin vaiheen perusteet mietitään etukäteen. On silti muistettava, että vasta ratkaisun
jälkeen on tiedossa ensimmäisen kopin y + -arvot, joten hyvin usein hilaa joudutaan korjaamaan.
• Virtaustilanteen määrittely. Tässä vaiheessa asetetaan täsmällisesti reuna- ja
alkuehdot sekä aineominaisuudet.
9.1. SIMULOINTITEHTÄVÄ
241
• Numeerinen ratkaisu saattaa viedä aikaa viikkoja. Perussyy pitkään laskentaaikaan on se, että tietokoneiden kapasiteetin kasvaessa halutaan laskea monimutkaisempia tilanteita paremmalla resoluutiolla. Hyvin usein tulee konvergenssivaikeuksia ja laskentavaihe venyy. Siksi on hyvä tehdä samanaikaisesti
myös muita simulointitehtävän vaiheita, kuten alustavaa tulosten analysointia.
• Tulosten analyysi on periaatteessa viimeinen vaihe, mutta kuten edellä todettiin, laskennan kuluessa tilanteen kehittymistä on syytä seurata päivittäin.
aloitus
ongelman
alkutarkastelu
virtaustilanteen
määrittely virhe?
hilan generointi
huono laskenta−
hila
virtauksen määrittely
(reunaehdot)
virtaustilanne ei
ole määritetty
oikein?
numeerisen ratkaisun
laskenta
laskentatulosten
analysointi
tulokset eivät ole
kunnossa
lopetus
Kuva 9.1: Virtaussimulointitehtävän kulku.
Ihanteellisessa tilanteessa voisimme kuvitella tehtävän suorittamisen etenevän
kronologisesti vaiheittain. Käytännössä joudutaan tekemään erilaisia tarkistuksia ja
palaamaan edellisiin vaiheisiin. Prosessia havainnollistetaan kuvassa 9.1. Etenemällä kuvan osoittamalla tavalla, tehdään tarkistuksia, jotka muodostavat eräänlaisen
laatujärjestelmän. Vaikka kuvassa ei olekaan tarkastuspisteitä hilangenerointivaiheen jälkeen, on sellainenkin syytä tehdä. Tätä tarkoitusta varten hilangenerointivaiheesta saadaan informaatiota eri muodoissa. Joskus hilan tarkistaminen tai ainakin osa siitä tapahtuu ratkaisuvaiheen initialisoinnissa. Mikäli hilaa koskevissa
9.2. HARJOITUSHÄVITTÄJÄN YLÖSVETOTILANNE
242
tunnusluvuissa näkyy jotain epäilyttävää, on palattava takaisinpäin jo ennen varsinaisen laskennan aloittamista.
Suurin osa ongelmista paljastuu vasta tuloksia analysoitaessa, jolloin on tehtävä
kuvan 9.1 tarkastusketju. Monet virheet näkyvät jo konvergoitumattomissa tuloksissa, minkä vuoksi laskennan tiivis seuraaminen säästää usein paljon turhaa laskentaaikaa.
Seuraavaksi käydään lyhyesti läpi muutama simulointitehtävä. Laskuja ei ole
tehty FLUENTilla eikä edellä esitettyjä työvaiheita selosteta tarkasti. Mahdollisuuksien mukaan tutkitaan, miten tehtäväasettelu tapahtuisi FLUENT-ympäristössä.
9.2 Harjoitushävittäjän ylösvetotilanne
Tehtävän tarkoituksena on selvittää BAe Hawk Mk 51 suihkuharjoituskoneen rakenteisiin kohdistuvat aerodynaamiset kuormat lentotiloissa, joissa kone on ns. ylösvetotilanteessa. Tyypillisesti tällaisia tilanteita esiintyy ilmataisteluharjoituksen alkaessa, kun ohjaaja havaitsee ’viholliskoneen’ ja kaartaa sen perään. Laskennallisesti on tarkoitus määrittää pintoihin vaikuttavat voimat, ts. paine- ja kitkavoimat.
( x 0 , y0 )
( x 1 , y1 )
R
α
q
8
V
α
( x cg , ycg )
Kuva 9.2: Ylösvetotilanteessa oleva lentokone voidaan analysoida laittamalla laskentahila
pyörimisliikkeeseen.
9.2. HARJOITUSHÄVITTÄJÄN YLÖSVETOTILANNE
243
Ylösvedossa koneen lentorata on kaareva (kts. kuva 9.2). Koneeseen kohdistuu
tällöin keskeiskiihtyvyys, jonka suuruus on
v2
(9.1)
R
missä R on kaarevuussäde. Yleensä aerodynamiikassa käytetään ns. tuulitunnelia=
koordinaatistoa, jossa kone, ts. laskentahila, on paikallaan ja reunaehtona annetaan
vapaan virtauksen nopeus. Ylösvetotilanteessa ei tätä laskentatapaa voida käyttää
kiihtyvän liikkeen vuoksi. Simulointi voidaan tehdä asettamalla laskentahila pyörimään. Koneen kohtauskulma otetaan huomioon siirtämällä kaarevuussäteen keskipistettä pystysuorasta hieman eteenpäin kuvan 9.2 osoittamalla tavalla. FLUENTissakin laskenta onnistuu tällä tavoin. Tapauksessa käytetään siis pyörivää hilaa ja
karteesisia nopeuskomponentteja.
Virtaustilanteessa Machin luku on hieman alle yhden ja Reynoldsin luku luokkaa Re ≈ O(107). Tilanteeseen sopii siis FLUENTin tiheyspohjainen ratkaisu. Koska kyseessä on rajakerrostyyppinen virtaus, käytetään mieluiten pienen Reynoldsin
luvun mallia. Myös FLUENTin kaksikerros -mallia voidaan periaatteessa käyttää.
Eräs ongelma, jota ei aina tiedosteta, on rajakerroksen transitio laminaarista turbulenttiin virtaukseen. Nyt Reynoldsin luku on niin suuri, että transitio tapahtuu
lähellä siiven johtoreunaa. Yleinen käytäntö on tällöin laskea siipi kokonaan turbulenttina.
Kuva 9.3: Lentokoneen pinta alkuperäisessä IGES-formaatissa (vasemmalla) ja laskentaa
varten tehty pintahila (oikealla).
Aerodynamiikassa tarkkuusvaatimus on hyvin suuri ja yleensä pyritään käyttämään rakenteellista hilaa rajakerrosten vuoksi. Toisena mahdollisuutena on raja-
9.2. HARJOITUSHÄVITTÄJÄN YLÖSVETOTILANNE
244
kerroksen kuvaaminen prismoilla. Vaikka lentokone on suhteellisen monimutkainen, sille on mahdollista laatia monilohkoinen rakenteellinen hila. FLUENTissa ei
ole lohkokäsitettä ja hila voidaan laatia tällöin hieman eri periaatteella, vaikka pääpyrkimyksenä olisi ainakin pintojen lähellä käyttää rakenteellista topologiaa. Ongelmallisissa tilanteissa voidaan käyttää rakenteetonta hilaa. Hilan laadinta alkaa
pinnoista, jotka näin monimutkaisen geometrian tapauksessa on saatava valmiina.
Tässä tilanteessa on geometria saatu IGES-formaatissa. Laskentahila on laadittu
IGG-ohjelmistolla (kts. kuva 9.3). Hilan laadinnassa on arvioitu ensimmäisen pintaa vasten olevan kopin korkeutta. Korkeus ja vasta ratkaisun jälkeen tiedossa oleva
y + -jakauma ovat kuvassa 9.4. Tässä tapauksessa y + -jakauma on varsin hyvin onnis-
Kuva 9.4: Lentokoneen hilan ensimmäisen kopin korkeus (vasemmalla) ja y + -jakauma (oikealla).
tunut. Kuvan skaala on välillä 0...3 ja y + -arvot ovat ainoastaan paikallisesti lähellä
ylärajaa.
Virtausyhtälöiden ratkaisussa kannattaa ensimmäisenä vaihtoehtona kokeilla implisiittistä ratkaisua. Mikäli sen kanssa syntyy konvergenssiongelmia, on mahdollista yrittää myös eksplisiittistä ratkaisua. Kuten aiemmin todettiin, mikäli eksplisiittisen ratkaisun yhteydessä voidaan käyttää useaa monihilatasoa, se saattaa olla jopa
implisiittistä ratkaisua tehokkaampi.
Diskretointi on tehtävä toisen kertaluvun menetelmällä. Tästä saattaa FLUENTin yhteydessä syntyä ongelmia, koska kaikki FLUENTin diskretointitavat ovat ra-
9.2. HARJOITUSHÄVITTÄJÄN YLÖSVETOTILANNE
245
Kuva 9.5: Hawk-harjoitushävittäjälle laskettuja virtaviivoja, pinnalle pakotettuja virtaviivoja sekä painejakauma.
joitettuja, ts. niihin lisätään hieman numeerista vaimennusta keinotekoisesti. Tämä
on haitallista erityisesti pintoja vastaan kohtisuorassa suunnassa. Kuvassa 9.5 olevat tulokset on laskettu FINFLO-ohjelmalla käyttäen pintaa vastaan kohtisuorassa
suunnassa muodollisesti kolmannen kertaluvun menetelmää, päävirtauksen suunnassa toisen kertaluvun rajoitettua diskretointia ja kolmannessa suunnassa toisen
kertaluvun rajoittamatonta menettelyä. Tällä tavoin saatiin rajakerros tarkasti lasketuksi ja numeriikasta aiheutuva keinotekoinen kitka mahdollisimman pieneksi.
Yhteenvetona voidaan todeta, että periaatteessa FLUENTilla voidaan laskea ylösvetotilanteessa oleva lentokone transoonisella nopeusalueella. Ongelmia melko todennäköisesti kuitenkin aiheutuisi siitä, että vaikka FLUENTissa periaatteessa on
mallit tilanteen laskemiselle, ohjelmaa ei ole kovin paljon sovellettu tällaisiin tarkoituksiin. Selvinä ongelmina on nähtävissä sopivien diskretointimenetelmien puute, turbulenssimallit ja se, ettei niitä ole käytetty transoonisella alueella, mahdolliset
konvergointiongelmat ja melko todennäköisesti myös tulosten jälkikäsittely.
9.3. SAVUPELTILAITTEEN PAINEHÄVIÖIDEN SELVITTÄMINEN
246
9.3 Savupeltilaitteen painehäviöiden selvittäminen
Kuva 9.6: Monilapainen savupelti.
Toisena laskentatapauksena tarkastellaan savupeltilaitteen (kts. kuva 9.6) painehäviöitä. Savukaasukanavan suunnittelussa tarvitaan peltilaitteen kertavastuskerroin. Periaatteessa tämä voidaan määrittää kokeellisesti, mutta erityisesti pellin ollessa auki-asennossa painehäviö on pieni ja vaikea mitata. Peltejä valmistetaan vaihtelevilla lapamäärillä ja valmistaja halusi tietää myös lapamäärän vaikutuksen vastukseen. Lapamäärä puolestaan vaikuttaa yksittäisten lapojen sivusuhteeseen ja Reynoldsin lukuun. Tilanteeseen vaikuttaa myös virtauskanavan hydrauliseen halkaisijaan referoitu Reynoldsin luku.
Etukäteen voitiin arvella lapamäärän ja sivusuhteen olevan tärkeitä parametreja. Lisäksi painehäviö on tietenkin lapojen aukeamiskulman funktio. Hyvin pian havaittiin, että tehtävästä tulee varsin mittava. Lapakulman vaikutuksen selvittämiseksi ajettiin kahdeksan eri lapamäärää neljällä kulmalla kaksiulotteisena. Ajomääräksi
tuli 32, mutta kaksidimensioisina ne olivat kevyitä. Reynoldsin luvun vaikutuksen
selvittämiseksi ajettiin yksi kaksiulotteinen tapaus viidellä eri Reynoldsin luvulla
tilanteessa, jossa pelti on täysin auki. Sivusuhteen vaikutusta tutkittiin laskemalla
9.3. SAVUPELTILAITTEEN PAINEHÄVIÖIDEN SELVITTÄMINEN
Tulppa−
virtaus
247
Symmetriataso
Lapa
Symmetriataso
50 m
Lapa
Kuva 9.7: 3D reunaehdot savupellille.
kolme eri sivusuhdetta kolmidimensioisena. Lisäksi käytössä oli ääretön sivusuhde
tehdyistä kaksidimensioisista simuloinneista. Kokonaisuudessaan tarvittavien tietokoneajojen määrä nousi yllättävän suureksi ja työn kuluessa ilmeni vielä tarvetta
tehdä lisälaskuja pelkälle virtauskanavalle.
Reunaehdoista tiedettiin, että keskimääräinen virtausnopeus on 10 m/s. Koska käytännön tilanteissa nopeusjakauma voi olla hyvinkin erilainen, käytettiin tässä täysin kehittyneen kanavavirtauksen jakaumia nopeuden ja turbulenssisuureiden
osalta. Kaksiulotteisissa laskuissa kanavavirtauksen profiili laskettiin erillisellä ohjelmalla ja asetettiin se reunaehdoksi. Koska vastaavaa kolmidimensioista kanavavirtauksen ratkaisutapaa ei aikoinaan ollut käytettävissä, jouduttiin kolmidimensioisissa laskuissa säätölaitteen eteen asettamaan 50 m pitkä tyhjä kanava, jonka toiseen
päähän asetettiin tasainen nopeusjakauma. Kolmidimensioisissa tapauksissa käytettiin symmetriaehtoja siten, että vain yksi neljäsosa kanavasta mallinnettiin. Kolmidimensioisen tilanteen reunaehdot on esitetty kuvassa 9.7. Symmetriaehtoa käytettiin myös kuvatasoa vasten kohtisuorassa olevassa suunnassa. FLUENTissa voitaisiin käyttää samoja reunaehtoja, sisäänvirtauksen osalta massavirtaehtoa ja ulosvirtauksessa painereunaehtoa. Nyt tehdyissä analyyseissä ei ulostuloreunalla tapahtunut takaisin virtausta, joten tässä tapauksessa olisi toiminut FLUENTin yksinkertaisempikin vaihtoehto, jossa kaikki suureet annetaan vain sisäänvirtausreunalla.
Laskentahila savupellin läheisyydessä on esitetty kuvassa 9.8. Tässäkin tapauksessa käytettiin pienen Reynoldsin luvun turbulenssimallia, joten hila on voimakkaasti tihennetty seinille. Pellin pituuteen referoitu Reynoldsin luku on luokkaa
O(105). Virtaus pellin rajakerroksessa on osittain laminaaria ja tällöin olisi arvioitava transitiokohta. Nyt esillä olevassa tapauksessa muut epävarmuustekijät ovat
9.3. SAVUPELTILAITTEEN PAINEHÄVIÖIDEN SELVITTÄMINEN
248
Kuva 9.8: Savupellin laskentahila, jossa lohkot on vedetty erilleen.
hyvin suuria ja laskuilla pyritään löytämään enemmän trendejä kuin laskemaan painehäviön eri osatekijät absoluuttisen tarkasti.
Suurilla lapakulmilla virtaus irtoaa ja on mitä ilmeisemmin voimakkaasti ajasta riippuva. Jotta tehtävä pysyisi laajuudeltaan jonkinlaisissa järkevissä rajoissa,
oli perusteltua laskea tämäkin tapaus kokonaan turbulenttina. Kuvasta 9.8 nähdään
myös, että virtauskanavaa ympäröi kapea lista, johon pelti nojaa kiinni-asennossa.
Tämä lista on tärkeää mallintaa tarkasti, koska pellin ollessa täysin auki -asennnossa,
painehäviöstä suuri osa aiheutuu listasta. Painehäviö on sinänsä hyvin pieni aukiasennossa.
Kaksidimensioisista simuloinneista on esimerkkinä kuvan 9.9 nopeusjakauma
ja virtaviivat. Lapojen yläpinnoilla virtaus irtoaa voimakkaasti kuvan lapakulmalla.
Virtaus irtoaa myös kanavan seinällä olevan listan jälkeen ja tätä ilmiötä näyttäisi voimistavan lähellä sijaitseva lapa. Kyseessä on tasapainotilan lasku, mutta silti
lavoista alavirtaan on irtoavien pyörteiden kaltaisia ilmiöitä. Tilanne pitäisi periaatteessa laskea ajan suhteen tarkasti. Koska nyt ollaan lähinnä kiinnostuneita voimakertoimista, oli mahdollista käyttää iteraatiokierrosten suhteen keskiarvotettuja
suureita.
9.3. SAVUPELTILAITTEEN PAINEHÄVIÖIDEN SELVITTÄMINEN
249
Kuva 9.9: Virtaviivat ja nopeus lapojen ympäristössä.
Savupellin laskenta on tilanne, jossa tulee kaikilla ohjelmilla konvergenssiongelmia, koska fysikaalinen tilanne oikeastaan edellyttää sitä. Laskentaa ei pidä ryhtyä väkisin saattamaan tasapainotilaan iteraatioparametreja muuttamalla tai vaihtamalla esimerkiksi ensimmäisen kertaluvun diskretointiin. Jos virtausjakaumien suhteen tarvitaan tarkempaa kuin kvalitatiivista tietoa, on tehtävä ajan suhteen tarkkoja laskuja. Usein tällöin käy niin, että tasapainotilan simuloinnissa voimakkaana
esiintynyt värähtely pienenee, joskus jopa häviää kokonaan. Myös turbulenssimalli
vaikuttaa tilanteeseen. Jollain turbulenssimallilla tilanne saattaa konvergoida tasapainotilaan, jollain toisella taas ei. Esillä olevissa laskuissa käytettiin Chienin pienen Reynoldsin luvun k − ǫ-mallia. Jotta eri simulointien tulokset olisivat keskenään vertailukelpoisia, on käytettävä systemaattisesti samaa mallia kaikissa ajoissa.
Turbulenssimallia ei siis pidä vaihtaa sen mukaan konvergoiko laskenta tasapainotilaan. Sen sijaan, jos ajo jatkuvasti kaatuu, on ehkä vaihdettava johonkin robustimpaan turbulenssimalliin ja käytettävä sen jälkeen systemaattisesti sitä. Turbulenssin
ja ajasta riippuvien ilmiöiden välillä on Reynolds-keskiarvotetuissa simuloinneissa
vaikea tehdä eroa ja tältä osin tehdään simulointimalleissa paljon kehitystyötä.
Koska tässä tilanteessa on useita parametreja ja lisäksi virtauskanava ja sitä kier-
9.4. RADIAALIKOMPRESSORIN AJASTA RIIPPUVA SIMULOINTI
250
tävä lista, osoittautui yksinkertaisen laskentakaavan laadinta kertavastuskertoimelle
vaikeaksi. Tehtävässä jouduttiin siten palamaan useasti ’alkuideointi’-vaiheeseen.
Kertavastuskerroin siipien lukumäärän ja lapakulman funktiona on kuvassa 9.9.
Vastuskertoimissa on hieman hajontaa, mikä osittain johtuu siitä, ettei yleensä ole
saavutettu tasapainotilaa, vaan käytetään keskiarvotettuja kertoimia. Siipien lukumäärällä havaitaan olevan jonkinasteinen, mutta heikko vaikutus vastuskertoimeen.
Tarkempi tutkimus osoitti, että siipimäärä ei suoraan vaikuta, vaan siipimäärän kasvaessa kanavaa reunustavan listan osuus vastukseen heikkenee. Lasketuista tuloksista oli lopulta mahdollista saada käyttökelpoinen menettely kertavastuskertoimen
määritykseen hajottamalla kertavastus eri tekijöistä aiheutuviin komponentteihin.
Kertavastus riippuu lapakulmasta, lavan sivusuhteesta, säätölaitteen kehyksen listan Reynoldsin luvusta ja säätölaitteen lavan jänteen Reynoldsin luvusta.
Kuva 9.10: Laskettu kertavastuskerroin lapamäärän funktiona.
9.4 Radiaalikompressorin ajasta riippuva simulointi
Aiemmin todettiin pyörimisliikkeen yhteydessä, että FLUENT on hyvin varustettu reunaehtojen ja muidenkin mallien osalta pyöriviä virtauksia varten. Seuraava-
9.4. RADIAALIKOMPRESSORIN AJASTA RIIPPUVA SIMULOINTI
251
o
boundary condition
measurement 1 m
before impeller
Pipe inlet
section height
45 (1)
o
R pipe
90 (2)
Paineen
mittaus
R inlet
Conical inlet
section height
modelled
inlet
Inlet height
o
J
180 (3)
Outlet as modelled
(not in scale)
Sh
ub
ud
H
ro
K
Splitter LE
o
(5) 360
Blade LE
Volute
Diffusor
o
270 (4)
R (impeller)
R (Stator)
Volute
o
o
ω
Kuva 9.11: Radiaalikompressorin geometria.
na esimerkkinä tarkastellaan radiaalikompressorin laskentaa ajasta riippuvassa tilanteessa. Kompressorin geometria on esitetty kuvassa 9.11. Yksinkertaisimmassa
simulointitavassa mallinnetaan yksi kompressorin siipisola ja jonkin matkaa tuloputkea ja diffuusoria reunaehtojen asettamista varten. Tällöin laskenta voidaan tehdä tasapainotilan oletuksella, koska reunaehdot ovat täysin symmetriset. Nyt esillä
olevassa esimerkissä tarkoituksena oli tutkia ajasta riippuvia ilmiöitä. Kompressorin
siivistöä seuraa diffuusori ja sen ulkopuolella spiraali (’volute’). Roottorin ulosvirtausreunalla tilanne ei ole todellisuudessa symmetrinen, vaan vaihtelee kehäkulman
funktiona. Ilmeisesti virtaus on roottorissakin sykkivää.
Ajasta riippuva tilanne voidaan laskea myös kvasistaattisena, kuten aiemmin on
jo esitetty. Reunaehtotyyppejä on silloin kaksi (kts. luku 8.3), toisessa suoritetaan
kehän suunnassa keskiarvottaminen, toisessa tilanne jäädytetään paikalleen. Nyrkkisääntönä on, että jos siivellinen staattori näkee roottoripuolen virtauksen nopeasti
sykkivänä, niin on syytä tehdä keskiarvottaminen kehän suunnassa. Esillä olevassa geometriassa staattori on siivetön, jolloin tarkemmalta kvasistaattiselta approksimaatiolta vaikuttaisi se, ettei keskiarvottamista tehdä. Tässä tapauksessa simulointi
aloitettiin näin tehdystä kvasistaattisesta laskusta. FLUENTissa voidaan tehdä täsmälleen samoin. Lasketaan jälleen pyörivässä koordinaatistossa (karteesisilla nopeuskomponenteilla!) kvasistaattisesti tasapainotila, josta käynnistetään ajan suh-
9.4. RADIAALIKOMPRESSORIN AJASTA RIIPPUVA SIMULOINTI
252
teen tarkka laskenta.
6
7
5
4
2
22
20 21
9
8
25
3
17 19
1
18
14
12
13
16
15
10
11
23
24
Sliding mesh boundary
Kuva 9.12: Radiaalikompressorin lohkojako.
Kompressorin geometria on vaikeahko, mutta monilohkoinen rakenteellinen hila on kuitenkin mahdollinen ja sellaista voitiin käyttää tässäkin. Suurimmat hankaluudet rakenteellisuuden vuoksi tulevat ns. kielen alueella. Kuvassa 9.12 kielen alue
näkyy numeron 25 oikealla puolella. Tässä kohden FLUENTissa kannattaisi ehkä
käyttää rakenteetonta hilaa. Koska jälleen kyseessä on rajakerrostyyppinen virtaus,
seiniä vastaan kohtisuorissa suunnissa käytetään silloin rakenteellista tyyppiä olevaa topologiaa. Kielen alueen ulkopuolella geometria on säännöllinen, joten hilangenerointi ei tuota ongelmia. Esimerkissä käytetty pintahila on nähtävissä kuvassa 9.13. Kuvaan, samoin kuin kuvaan 9.11 on merkitty myös paineen mittausreiät,
joista oli saatavilla mitattuja paineita validointiin.
Hilan suunnitteluvaiheessa kannattaa topologiasta tehdä periaatekaavio, jossa
tulevat myös reunaehdot esille (kts. kuva 9.14). Kyseisessä kompressorissa on seitsemän siipisolaa, jotka on puolesta välistä jaettu lyhyemmällä siivellä (’splitterillä’) kahteen osaan. Rakenteellisella koodilla tulee näin ollen lohkojen lukumääräksi
suuri, tässä tapauksessa 25. Roottorin mukana oleva osa hilasta pyörii (myös kvasistaattisessa analyysissä), loppuosan ollessa paikallaan. Ajan suhteen tarkassa laskennassa näiden välillä on liukuhila (sliding mesh, luku 8.4). Sisäänvirtausreunalla
käytettiin massavirtareunaehtoa ja ulosvirtausreunalla lohkossa 25 painereunaehtoa.
Sisäänvirtausreunalla voitaisiin käyttää myös paine-ehtoa, jolloin kompressorin yli
oleva massavirta laskettaisiin paine-eron ollessa kiinnitetty. Usein halutaan kiinnittää massavirta ja annetaan koodin laskea paine-ero. Tätä tapaa noudatettiin tässäkin
9.4. RADIAALIKOMPRESSORIN AJASTA RIIPPUVA SIMULOINTI
253
5
1
4
3
2
Kuva 9.13: Radiaalikompressorin pintahila ja paineen mittausreiät.
1
Rotating mesh
Inlet, Block 1
2
5
20
A
B
A
B
3
4
6
7
(23)
21 22
Sliding mesh boundary
Stationary
(24)
Outlet
}
}
1’
Impeller
2
Diffusor
3
Volute, Block 25
5
Kuva 9.14: Radiaalikompressorin lohkojaon topologia.
simuloinnissa.
Ajan suhteen tarkka laskenta on vielä nykyisilläkin tietokoneresursseilla aikaa
vievää. Tässä tapauksessa on käytettävä toisen kertaluvun implisiittistä aikaintegrointimenetelmää. Koska kompressorivirtaus on kokoonpuristuvaa, FLUENTissa käytettäisiin tiheyspohjaista ratkaisijaa. Aika-askeleen sisällä on suoritettava iterointeja implisiittisyyden vuoksi, mikä pidentää laskenta-aikaa. Aika-askeleen on oltava varsin lyhyt, jotta halutut ilmiöt saataisiin esille. Tässä tapauksessa käytettiin
∆t = 2µs, jolloin roottori pyörähtää 0,216◦ kulman aika-askelta kohden. Alkutilanteena käytettiin kvasistaattisen laskennan tulosta, jonka jälkeen syklinen virtaustilanne saatiin stabiloitumaan, kun oli laskettu noin 3,8 pyörähdystä. Tämä vaati 6395
9.4. RADIAALIKOMPRESSORIN AJASTA RIIPPUVA SIMULOINTI
254
aika-askelta.
Kuva 9.15: Hetkellinen staattisen paineen jakauma radiaalikompressorissa.
Ajasta riippuvien tilanteiden simulointi on hankalaa paitsi pitkän laskenta-ajan,
myös tulosten käsittelyn ja hallinnan osalta. Laskentatuloksia tulee niin paljon, ettei niitä kaikkia voida aina tallentaa jokaiselta aika-askeleelta. Tulosten havainnollistamisessa animaatiot ovat erittäin hyvä keino, mutta niidenkin yhteydessä pitää
etukäteen pystyä spesifioimaan pinnat ja suureet, joita halutaan jatkossa käsitellä. Monien suureiden osalta on myös laskettava keskiarvoja, joita voidaan verrata
kvasistaattisiin tuloksiin. Ajan suhteen tarkat analyysit olivat aiemmin enemmän
tutkimuksen kuin suunnittelun apuväline, mutta nyt niitä tehdään enenevässä määrin myös käytännön laskennassa. Niiden avulla saadaan mm. tietoa kvasistaattisten
laskujen tarkkuudesta ja ratkaisussa olevien ajan suhteen tapahtuvien värähtelyjen
amplitudeista. Nyt tehdystä laskennasta ilmeni, että siipisola kokee huomattavan
suuren paine- ja vastaavan massavirtaheilahtelun pyörimisliikkeen aikana. Suurin
syy tähän on paineen vaihtelu diffuusorissa kehän suuntaisen aseman funktiona.
Diffuusorissa paikallisesti tapahtuva paineen muutos siiven ohituksen tapahtuessa
on tätä ilmiötä huomattavasti vähäisempi. Esimerkkinä laskentatuloksista on hetkellinen liikemäärän jakautuma kuvassa 9.15. Paineen vertailussa mittaustuloksiin
saatiin kohtalaisen hyvä yhtäpitävyys.
9.5. KANAVAVIRTAUS ISOJEN PYÖRTEIDEN MENETELMÄLLÄ
255
9.5 Kanavavirtaus isojen pyörteiden menetelmällä
Viimeisenä esimerkkinä tarkastellaan yksinkertaisinta mahdollista isojen pyörteiden menetelmällä laskettavaa tilannetta, kanavavirtausta. Tilanteen tekee yksinkertaiseksi se, että reunaehdot voidaan käsitellä joko periodisina tai kiinteän pinnan
ehtoina. Reunaehtojen asettaminen on LESissä huomattavasti RANS-yhtälöitä monimutkaisempaa, koska reunoilla tapahtuva suureiden värähtely vaikuttaa koko laskenta-alueen ratkaisuun. Esimerkiksi kanavavirtauksella sisääntuloehdon on edustettava turbulenttia virtausta. Tällainen saadaan helposti kuvatuksi ottamalla nopeusjakauma ulosvirtauksesta, edellyttäen, että reunojen välimatka on riittävän pitkä. FLUENTissa on LES-mallinnus mahdollista, mutta tilanteen vaatimukset on
otettava huomioon ennen kuin laskentaan ryhdytään! Mikä tahansa simulointi ei ole
LESiä, vaan lähinnä jonkinlaista RANS-laskentaa. FLUENTissa on isojen pyörteiden menetelmää varten reunaehto, jossa tulovirtausta perturboidaan satunnaisesti.
Menettelyä voidaan käyttää, mutta sisäänvirtausreuna on tällöin asetettava hyvin
kauaksi varsinaisesta laskenta-alueesta, muuten virtaus ei ole täysin kehittynyt tar-
log(E(k))
kasteltavalla alueella.
inertial
subrange
large scale
eddies
dissipation
subrange
η
l
log(k)
Kuva 9.16: Turbulenssin kineettisen energian spektri.
Palautetaan mieliin suoran simuloinnin ja isojen pyörteiden menetelmän perusideat: Suorassa simuloinnissa laskennan on pystyttävä kuvaamaan koko turbulenssin spektri (kuva 9.16). LESissä osa spektriä mallinnetaan. Joskus käytetään hyväksi numeriikan vaimennusominaisuuksia, mutta tämä laskentatapa on hyvin kiistanalainen. Perinteisesti on pyritty siihen, että numeriikka ei saa vaimentaa eikä
9.5. KANAVAVIRTAUS ISOJEN PYÖRTEIDEN MENETELMÄLLÄ
256
vääristää ajan suhteen tärkeitä värähtelyjä. Tavanomaisilla numeerisilla menetelmillä on hyvin vaikea erottaa numeerista virhettä ja fysikaalista vaimennusta toisistaan. Tämä onkin johtanut siihen, että sekä suoraa simulointia että LESiä on harrastettu paljolti ns. spektraalimenetelmillä. Numeeriikan vaikutusta voidaan tukia
aaltoluvun avulla. Aaltoluku tavallaan vastaa turbulenssin pituusskaalan käänteisarvoa. Esimerkiksi funktion sin kx derivaatta on k cos kx ja se saadaan tarkasti
spektraalimenetelmällä. Sen sijaan tarkallakin differenssimenetelmällä, kuten neljännen kertaluvun keskeisdifferensiillä, derivaatan kertoimeen k tulee virhettä (kuva
9.17). Spektraalimenetelmässä tätä virhettä ei tapahdu, minkä vuoksi monet pitävät
Kuva 9.17: Kahden differenssimenetelmän ja spektraalimenetelmän vaimennukset aaltoluvun funktiona.
sitä ideaalisena teoreettisten turbulenssitutkimusten yhteydessä. Tällöin reunaehtojen on aina oltava periodisia, joten käytännön monimutkaisia tapauksia ei tietenkään voida laskea. Muilla menetelmillä esiintyy aina virhettä, mikä tekee LESin tai
DNS:n tyyppiset simuloinnit erittäin vaikeiksi ja tulokset tulkinnanvaraisiksi.
Tarkastellaan seuraavaksi kanavavirtauksen laskentaa. Kanavalla tarvitaan fysikaaliset ehdot vain ’katolle’ ja ’lattialle’, muissa suunnissa voidaan käyttää periodisuutta. Tämä edellyttää laskenta-alueelta riittävää kokoa, jotta sisäänvirtausreunan reunaehto ei näy enää turbulenssin rakenteessa lähestyttäessä ulosvirtausreunaa.
Virtauksen suunnassa tapahtuu painehäviö, joka on otettava huomioon. FLUENTissa on mahdollisuus juuri tämän tyyppiseen reunaehtoon, mutta siinä paine-ero kanavan päiden välillä on kiinnitettävä, mikä RANS-laskuissa on järkevää, jos massavirta annetaan. Sen sijaan ajasta riippuvassa tilanteessa paine-erokin voi väräh-
9.5. KANAVAVIRTAUS ISOJEN PYÖRTEIDEN MENETELMÄLLÄ
257
dellä. Nyt esillä olevassa laskennassa painehäviötä säädettiin PID-säätimellä, jossa
parametrinä oli massavirran suhteellinen virhe. Toteutunut painehäviö on kuvassa 9.18. Virtausta vastaan kohtisuorassa suunnassa käytettiin periodisuutta eli siinä
suunnassa virtaustilanne toistaa itseään äärettömyyteen asti. Kanavan leveyden on
oltava riittävän suuri, jotta turbulenssin rakenteilla ei ole vuorovaikutusta reunojen
välillä.
Kuva 9.18: Suorassa simuloinnissa laskettu painehäviö kanavavirtauksessa.
Kuva 9.19: Suoralla simuloinnilla saatuja nopeusjakaumia kanavavirtaukselle.
Laskentamenetelmänä käytettiin eksplisiittistä toisen kertaluvun aikaintegrointia (Adams-Bashford). Nopeuteen referoidut Courantin luvut olivat välillä 0,1...0,2.
Jotta virtaus etenisi alueen päästä päähän, tarvittiin noin 500 askelta. Simuloinnissa
9.5. KANAVAVIRTAUS ISOJEN PYÖRTEIDEN MENETELMÄLLÄ
258
tehtiin kymmeniä tuhansia aika-askelia, joilta tulokset keskiarvotettiin. Jotta virtaus
saataisiin turbulentiksi mahdollisimman nopeasti, alkuarvoiksi annettiin voimakkaasti häiritty tilanne.
Numeerisen vaimennuksen välttämiseksi konvektiotermeille käytettiin keskeisdifferenssimenetelmää. Aiemmin FLUENTissa voitiin käyttää vain rajoittimella varustettuja toisen kertaluvun menetelmiä ja on epäselvää vaimentavatko ne ratkaisua
liiaksi, ts. voidaanko niitä edes käyttää isojen pyörteiden menetelmän yhteydessä.
Tämän vuoksi jo FLUENT 6 ohjelmassa oli mukana keskeisdifferenssimenetelmä,
jota tulisi käyttää LES-laskennassa.
Simulointi tehtiin sekä alihilamallin kanssa että ilman. Viimeksi mainitussa tapauksessa kyseessä onkin DNS, jos hilatiheys on riittävä. Harvemmilla hilatiheyksillä olisi puhuttava jonkinlaisesta näennäis-DNS:stä tai eräänlaisesta LESistä. Laskentaa tehtiin kolmella eri hilatiheydellä, joista harvin johti epätarkkaan tulokseen
ilman alihilamallia. Sen sijaan jo tiheydellä 32 × 64 × 32 saatiin varsin hyvä tulos
virtausnopeuksille, mutta turbulenssisuureissa virheet olivat merkittävämmät. Tiheimmällä tasolla saatiin käytännössä sama tulos kuin kirjallisuudessa raportoidussa suorassa simuloinnissa (Kim, Moin ja Moser). Keskiarvotetut nopeusjakaumat
eri hilatiheyksillä on esitetty kuvassa 9.19. Nopeusjakautumien muodosta voidaan
havaita virtauksen jo olevan selvästi turbulentin. Tiheimmillä hilatasoilla saadut tulokset ovat hyvässä sopusoinnussa logaritmisen nopeusjakauman kanssa.
Tässä tapauksessa kanavan korkeuteen referoitu Reynoldsin luku oli varsin pieni Re = 2800. Laskettaessa tilannetta alihilamallin avulla, tulos ei parantunut, vaan
huononi. Tähän on ehkä syynä käytetty hyvin alhainen Reynoldsin luku. Toisaalta
tämä kuvaa osaltaan LESin yhteydessä syntyviä hankaluuksia. Jotta suora tai isojen
pyörteiden simulointi kuvaisivat turbulenttia virtausta, on niiden toteutettava tietyt
ehdot. Mikäli lähdetään kyseisillä menetelmillä esimerkiksi FLUENTia käyttäen
laskemaan monimutkaisia tapauksia, on ehdottomasti ensin tarkistettava selviääkö
koodi edes kanavavirtauksesta alhaisella Reynoldsin luvulla. Laskennan on myös
näytettävä oikealta. Kuvassa 9.20 on esitetty hetkellinen nopeusjakauma, joka vastaa visuaalisesti mielikuviamme turbulentista virtauksesta. Jotta varmistutaan siitä,
että turbulenssi on oikein kuvattu, on syytä laskea myös Reynoldsin jännitykset.
Reynoldsin jännityksiä voidaan verrata alihilamallin antamiin jännityksiin. Kokonaisjännitykset ovat näiden summa. Hyväksyttävässä LESissä alihilajännitysten tulisi olla selvästi pienemmällä tasolla kuin nopeusheilahteluista lasketut jännitykset.
Laskentatuloksen tulisi myös toteuttaa turbulentin virtauksen spektri (kuva 9.16).
9.6. KERTAUS
259
Sekä simuloinnin suorittamiseen että tulosten analyysiin kohdistuu isojen pyörteiden menetelmässä ja suorassa simuloinnissa suuret RANS-simuloinneista poikkeavat vaatimukset.
Kuva 9.20: Hetkellinen nopeusjakauma kanavassa.
9.6 Kertaus
Simulointitehtävän vaiheet ovat:
• ongelman alkutarkastelu
• hilan generointi (esikäsittely)
• virtaustilanteen määrittely
• numeerinen ratkaisu
• laskentatulosten analysointi (jälkikäsittely)
9.6. KERTAUS
260
• aerodynaamisissa simuloinneissa on yleensä tuulitunnelikoordinaatisto, mutta ylösvetotilanteessa käytetään pyörivää hilaa
• monimutkaisen kappaleen pinta on saatava valmiina, esimerkiksi IGES-formaatissa
• kiinteää pintaa vasten kohtisuora suunta on tarkkuudeltaan kriittinen pienen
Reynoldsin luvun turbulenssimallin yhteydessä
• kanavavirtauksen reunaehtojen asettaminen on pohdittava tarkasti, muuten
reuna vaikuttaa laskenta-alueeseen häiritsevästi
• symmetriaehdolla saadaan usein laskentahila pienenemään yhteen neljäsosaan
• jos virtaus vaikuttaa ajasta riippuvalta, sitä ei pidä keinotekoisesti yrittää saattaa tasapainotilaan
• verratessa fysikaalisten parametrien vaikutusta on turbulenssimallin oltava sama
• ajasta riippuva simulointi kannattaa aloittaa kvasistaattisesta tuloksesta
• RANS-laskuissa aikaintegrointina voidaan käyttää vain toisen kertaluvun implisiittistä menetelmää. Isojen pyörteiden menetelmässä joidenkin termien osalta eksplisiittinen menetelmä saattaa olla tehokkaampi.
• reunaehtojen asettaminen LESissä on erittäin vaikeaa
• mikä tahansa ajasta riippuva lasku ei ole LESiä, vaan eräänlaista RANSsimulointia omituisella turbulenssimallilla
• ajasta riippuvissa tilanteissa paine-erotkin muuttuvat ajan funktiona, mikä
osaltaan vaikeuttaa reunaehtojen asettamista
• suorassa simuloinnissa tai LESissä kannattaa alkutilaksi usein antaa voimakkaasti häiritty virtaustilanne
• Ylävirtapainotteiset diskretointimenetelmät ovat huonoja LESin yhteydessä
• LES-laskentaa on aina ensin testattava jollain mahdollisimman yksinkertaisella esimerkillä
Päivitetty 11.3.2014