Rayleigh-Bénard konvektion
Förord
Föreliggande rapport är ett examenskandidatarbete som ingår i civilingenjörsut-bildningen vid Kungl. Tekniska Högskolan (KTH) och ges av mekanikinstitutionen, våren 2012. Det utförs inom ämnet strömningsmekanik och behandlar Rayleigh-Bénard-konvektion.
Projektgruppen vill tacka Philipp Schlatter för vägledning och goda råd. Stock-holm i Maj, 2012.
Abstract
Consider a uid being heated from below. The heating leads to an upward con-vective force that is counteracted by the viscous forces of the uid. If the concon-vective force is large enough in comparison to the viscous forces the uid will be put in an unstable state. This means that a small disturbance will give rise to a ow driven by a temperature gradient. This ow is characterised by a pattern of convection cells. The phenomenon is called Rayleigh-Bénard convection. An example of this can be seen when heating a pot of oil from below. A part of the contribution to the formation of these cells is attributed to the variation of surface tension due to heating. This contribution is of less signi cance when the uid layer is thicker. In this report the studied ow eld lies between two plates where the convective force drives the motion. The in uence of surface tension is eliminated since the uid lacks a free surface in this problem. The boundary between stability and insta-bility is investigated both theoretically, using simpli ed Navier-Stokes equations, and by simulation using a DNS-code with the program Simson (Chevalier et al., 2007). The simulation also makes it possible to see the shape of the convection cells. The results is presented in stability diagrams that describe how the stability boundary is aected by the wavelength, related to the wave number K, of the applied disturbance and the dimensionless Rayleigh number, Ra. The critical va-lue for the two parameters is found to be Ra = 1708 when K = 3.12
Sammanfattning
Betrakta en uid som värms underifrån. Uppvärmningen leder till en uppåtrik-tad konvektiv kraft som motverkas av viskösa krafter i uiden. Om de konvektiva krafterna är tillräckligt stora i förhållande till de viskösa krafterna försätts ui-den i ett instabilt tillstånd vilket leder till att en liten störning ger upphov till en strömning. Denna strömning karaktäriseras av ett mönster av så kallade konvek-tionsceller. Fenomenet kallas Rayleigh-Bénard-konvektion. Exempel på detta kan ses i en kastrull med olja som värms underifrån. Dock bidrar även oljans ytspän-ning till bildytspän-ningen av celler i detta exempel men bidraget blir mindre viktigt när uidskiktets tjocklek ökas.
I föreliggande rapport betraktas ett strömningsfält mellan två plattor där den ter-miska konvektionen är den drivande kraften. Eftersom uiden då inte har någon fri yta elimineras ytspänningens inverkan. Gränsen mellan stabilitet och instabilitet undersöks både teoretiskt, med hjälp av förenklingar av Navier-Stokes ekvationer, och genom en simulering i en DNS-kod i programmet Simson (Chevalier et al., 2007). Simuleringen möjliggör också att utseendet på konvektionscellerna kan ses. Resultat presenteras i stabilitetsdiagram som beskriver hur stabilitetsgränsen på-verkas av en störnings våglängd, kopplat till vågtalet K, och det dimensionslösa Rayleigh talet, Ra. Den kritiska stabilitetsgränsen för dessa två parametrar visar sig vara Ra = 1708 då K = 3, 12.
Innehåll
1 Inledning 1 2 Grundläggande ekvationer 2 3 Boussinesq approximation 2 3.1 Kontinuitetsekvationen . . . 3 3.2 Rörelsemängdsekvationen . . . 4 3.3 Energiekvationen . . . 54 Icke-dimensionalisering och dynamisk similaritet 6 4.1 Dimensionslöst med Ri, Re och P e . . . 7
4.2 Dimensionslöst med Ra och P r . . . 8
5 Instabilitet 9 5.1 Teori . . . 10
5.1.1 Små störningars utveckling . . . 11
5.1.2 Egensvängningsmetoden . . . 14
5.1.3 Gränsen till instabilitet . . . 15
6 Simulering med Simson 16 7 Resultat 17 7.1 Resultat från teorin . . . 17 7.2 Resultat från simulering . . . 19 7.3 Jämförelse . . . 24 8 Diskussion 24 9 Slutsats 25 Referenser 26 A Bilagor 27 A.1 Kontinuitetsekvationen . . . 27 A.2 Rörelsemängdsekvationen . . . 28
A.2.1 Rörelsemängdens bevarande . . . 28
A.2.2 Konstitutiv ekvation för Newtonsk uid . . . 28
A.2.3 Navier-Stokes ekvation . . . 30
A.3 Energiekvationen . . . 31
1 Inledning
När en kopp varm choklad svalnar kan det bildas ett mönster av celler på ytan. Dessa celler uppkommer på grund av temperaturskillnaden mellan chokladen i koppens botten och den svalare chokladen vid ytan. En densitetsskillnad uppkom-mer vilket leder till en ytkraft som vill driva en rörelse i chokladen. Det är denna kraft som bildar de så kallade konvektionscellerna på ytan. Viskositet och vär-mediusion motverkar denna rörelse, men om temperaturgradienten är tillräckligt stor bildas dessa celler. Cellernas uppkomst beror till stor del på ett fenomen som kallas Rayleigh-Bénard-konvektion, med bidrag från ytspänning.
Ett annat exempel på Rayleigh-Bénard-konvektion är vid bildningen av cumulus-moln i atmosfären. Då den undre delen av atmosfären värms upp av värmen från jordytan och den övre delen kyls av då den avger mer infraröd strålning än den absorberar, bildas en konvektiv rörelse (Weidauer, Pauluis & Schumacher, 2010). I detta arbete betraktas enbart Rayleigh-Bénard-konvektion där bland annat yt-spänningens bidrag elimineras då systemet som betraktas saknar en fri yta. Pro-blemet består av en uid som be nner sig mellan två horisontella, oändligt långa plattor. Temperaturerna vid plattorna hålls konstanta, och den undre plattan hålls vid en högre temperatur än den övre plattan. Beroende på uidens egenskaper och hur stor temperaturgradienten är, kan uiden försättas i ett instabilt läge och konvektionsceller kan bildas (se gur 1). Stabilitetsgränsen för Rayleigh-Bénard-konvektion, det vill säga övergången från vila till rörelse och därmed bildning av konvektionsceller, undersöks teoretiskt och med simuleringsprogrammet Simson. Detta görs med hjälp av ett antal dimensionslösa tal där Prandtltalet och Rayle-ightalet visar sig vara extra betydelsefulla i detta fall.
2 Grundläggande ekvationer
De tre grundläggande ekvationerna som används i detta arbete är kontinuitet-sekvationen, Navier-Stokes ekvationer och en energiekvation. Dessa har härletts utifrån principerna om mass-, rörelsemängds- och energikonservering. Här nedan är dessa ekvationer givna.
Kontinuitetsekvationen: Dρ Dt + ρ∇ · uuu = 0. (1) Rörelsemängdsekvationen: ρDuuu Dt = −∇p + ρggg + µ∇ 2uuu. (2) Här antas att Dρ
Dt ej är extremt stor som i stötvågor. Newtonsk uid antas och
viskositeten µ ses som konstant då temperaturskillnaden i uiden antas vara liten. Energiekvationen:
ρDe
Dt = −∇ · qqq − p∇ · uuu + φ. (3)
I dessa ekvationer är ρ densitet, uuu hastighet, p tryck, ggg gravitation (i z-led), µ vis-kositet, e inre energi, qqq värme öde och φ viskös energidissipation. I denna rapport används D för att beteckna materiell derivata. För härledningar se appendix A.1 - A.3.
3 Boussinesq approximation
1903 föreslog Boussinesq att man i samtliga termer i ekvationerna (1) - (3), utom den som är beroende av gravitationen g, kan bortse från densitetskillnader (Spiegel, Veronis 1960. s.442). Senare visade det sig att antaget var korrekt om nedanstående faktorer stämde (Pijush, Kundu & Cohen, 2008);
Machtalet är litet (lägre än 0,3).
Spridning av ljud-, och chockvågor inte behöver tas hänsyn till. Strömningen har en liten vertikal längdskala.
De tre översta punkterna har direkt koppling till tryck och kompressibilitet. För uider som strömmar med ett Machtal på över 0,3 är strömningen inte längre in-kompressibel och därmed kan inte densitetsskillnader försummas (K. Andersson et al. 2007, s.106). Denna approximation är väl anpassad för vätskor då de sällan kommer upp i så höga hastigheter.
Olika former av vågor är även de beroende av kompressibilitet. Om densitetsskill-nader inte tas hänsyn till, kan vågor fortplanta sig och nå oändligt höga hastigheter. Den tredje faktorn som är direkt beroende av kompressibilitet är begränsningen av strömningens vertikala längdskala. Här är skillnader i hydrostatiskt tryck vä-sentliga då de ger upphov till densitetsskillnader vid stora vertikala längdskalor. Den vertikala längdskalan måste följa sambandet
L = c
2
g ,
där c är ljudhastigheten i uiden. Den fjärde faktorn betraktas senare. Boussinesq approximation är alltså en förenkling av de tre grundekvationerna från ovanstående kapitel.
3.1 Kontinuitetsekvationen
Den vänstra termen i högerledet i kontinuitetsekvationen (1) är liten i jämförelse med hastighetsgradienten och kan därmed försummas. Observera att termen antas vara liten, inte noll. Det antas alltså inte att uiden är inkompressibel. Detta kan visas genom att betrakta små densitetsändringar δρ och temperaturändringar δT . Då gäller det linjära sambandet
δρ
ρ = −αδT,
där α är den termiska expansionskoe cienten och de nieras som
α = 1 ρ ∂ρ ∂T p .
För de esta uider är α av storleksordningen 10−4 10−3 K−1 (Pijush, Kundu &
Cohen, 2008, s.126). Därmed fås
αδT 1
1 ρ Dρ Dt ∇ · uuu ∼ 1 ρ u∂ρ∂x ∂u ∂x ∼ U ρ δρ L U L = δρ ρ = αδT 1,
där U och δT motsvarar en hastighets- och en temperaturskala över en längdskala L. Alltså är
∇ · uuu 1 ρ
Dρ Dt, vilket ger den förenklade kontinuitetsekvationen
∇ · uuu = 0. (4)
3.2 Rörelsemängdsekvationen
Anta ett referenssystem där ρ0 är densiteten och p0 är hydrostatiska trycket så
att ∇p0 = ρ0g. Beskriv systemet i en superposition av ett grundtillstånd och en
störning p = p0+ p0 samt ρ = ρ0+ ρ0, där (0) denoterar störningen. Ekvation (2)
kan då skrivas som
(ρ0+ ρ0)
Duuu
Dt = −∇p
0+ ρ0ggg + µ∇2uuu.
Genom att dividera denna med ρ0 fås
(1 + ρ 0 ρ0 )Duuu Dt = − 1 ρ0 ∇p0+ ρ 0 ρ0 g g g + ν∇2uuu.
Här är ν = µ/ρ0. Då små densitetsändringar betraktas är det givet att ρ0 ρ0
vilket innebär att kvoten ρ0/ρ
0 kan försummas i vänsterledet. Detta kan dock inte
göras när denna är multiplicerad med g, eftersom denna term beskriver ytkraftens inverkan och inte kan försummas.
Duuu Dt = − 1 ρ0 ∇p0+ ρ 0 ρ0 g gg + ν∇2uuu.
I denna ekvation betraktas alltså konvektion i vänsterled, i högerled nns de fak-torer som påverkar konvektionen;
−1 ρ0∇p
0 - Denna term behandlar tryckets inverkan.
ρ0
ρ0ggg - Denna term beskriver ytkraftens inverkan.
Komponentvis blir ekvationen Du Dt = − 1 ρ0 ∂p ∂x + ν∇ 2u, Dv Dt = − 1 ρ0 ∂p ∂y + ν∇ 2 v, Dw Dt = − 1 ρ0 ∂p ∂z − ρ ρ0 g + ν∇2w.
Observera att ytkraftens inverkan enbart verkar i z-led. Detta på grund av gra-vitationens riktining.
3.3 Energiekvationen
Vid Boussinesq approximation av energiekvationen (3) kan inte den förenklade kontinuitetsekvationen (4) appliceras på volymsexpansionstermen p∇ · uuu. Detta eftersom det inte är inkompressibla uider som betraktas, volymsexpansionstermen är i samma storleksordning som de andra termerna i ekvation (3). Därför används istället den vanliga kontinuitetsekvationen vilket ger
−p∇ · uuu = p ρ Dρ Dt ≈ p ρ ∂ρ ∂T DT Dt = −pα DT Dt. Antar ideal gas; p = ρRT , Cp− Cv = R samt α = T1. Detta ger
−p∇ · uuu = −ρRT αDT
Dt = −ρ(Cp− Cv) DT
Dt. Med detta antagande fås att
ρCp
DT
Dt = −∇ · qqq + φ, (5)
ty e = CvT. Ansätt Fouriers lag om värmeledning,
qqq = −k∇T. (6)
Ekvation (6) i ekvation (5) tillsammans med de nitionen av värmediusivitet
Genom att jämföra termen innehållande den viskösa dissipationen φ med vänster-ledet i ekvationen (7) fås φ ρCpDTDt ∼ 2µeijeij ρCpuj(∂x∂T j) ∼ µU 2/L2 ρ0CpU δT /L = ν Cp U δT L.
Denna kvot är mycket liten, ∼ 10−7 (Pijush, Kundu & Cohen, 2008, s.128). Därför
kan φ försummas vilket ger
DT
Dt = κ∇
2
T.
Efter Boussinesq approximation fås alltså ett ekvationssystem som beskriver rö-relsen i uiden; ∇ · uuu = 0, (8) Du Dt = − 1 ρ0 ∂p ∂x + ν∇ 2u, Dv Dt = − 1 ρ0 ∂p ∂y + ν∇ 2v, Dw Dt = − 1 ρ0 ∂p ∂z − ρ ρ0 g + ν∇2w, (9) DT Dt = κ∇ 2T. (10)
4 Icke-dimensionalisering och dynamisk similaritet
Strömning av uider i olika fall kan se lika ut även om de inte har samma densi-tet, viskosidensi-tet, längdskala etc. Detta antyder att det nns vissa kombinationer av de ingående storheterna som leder till lösningar som ser lika ut. Det visar sig att kravet för att lösningarna ska vara lika är att vissa dimensionslösa tal måste vara lika. Vilka tal som betraktas och vilka kombinationer av storheter som dessa tal består av beror på situationen som betraktas och vilka antaganden som görs. Att ha de aktuella ekvationerna på dimensionslös form är mycket användbart eftersom det är enklare att utföra experiment på modeller än på det fullskaliga problemet. Detta är vanligt vid utveckling av ygplan.
Ra = Flytkrafter Viskösa krafter = gαΓd4 κν , P r = Rörelsemängdsdiusion Värmediusion = ν κ, Re = Tröghetskrafter Viskösa krafter = U h ν , Ri = Flytkrafter Tröghetskrafter = α∆T hg U2 , P e = Re · P r = U h κ .
Notera att P r-talet endast beror på uidens egenskaper. För härledning av Rayleigh-talet, Ra, se bilaga A.4.
Här nedan följer två omskrivningar av ekvationerna (8) - (10). Dessa är centrala för stabilitetsundersökningen av Bénards problem. Genom att sedan jämföra de resulterande versionerna av Navier-Stokes visas att
Re = 1 P r, vilket används i simuleringen.
4.1 Dimensionslöst med Ri, Re och P e
Introducerar följande dimensionslösa parametrar w∗ = w U, t ∗ = U ht, p ∗ = p U2ρ 0 , z∗ = z h, T ∗ = T − T0 ∆T .
Här är U en karakteristisk hastighet för uiden, h är en karaktäristisk längd, ∆T är temperaturskillnaden och p är trycket.
Här betraktas z-komponenten av (9) vid omskrivningen av Navier-Stokes. Om-skrivningarna i de andra två riktningarna är helt analoga men saknar gravitations termen eftersom den är parallell med z-riktningen. Med ovan nämmnda dimen-sionslösa parametrar insatta i (9) fås
D (U w∗) D Uht∗ = − 1 ρ0 ∂ (U ρ0p∗) ∂ (hz∗) + α∆T T ∗ g + ν∂ 2(U w∗) ∂ (hz∗)2 .
Genom division med U2/h och identi ering av dimensionslösa tal fås Dw∗ Dt∗ = − ∂p∗ ∂z∗ + Ri · T ∗ + 1 Re∇ 2 w∗. (11)
På samma sätt skrivs (10) om vilket ger U ∆T h DT∗ Dt∗ = κ ∆T h2 ∇ 2T∗ . Förkortning med ∆T och identi ering av Pecletalet ger
DT∗ Dt∗ = 1 P e∇ 2T∗ .
Med denna framställning blir alltså de aktuella ekvationerna skrivna i vektornota-tion ∇ · uuu, Duuu∗ Dt∗ = −∇ · ppp ∗ + Ri · T∗eeez+ 1 Re∇ 2uuu∗ , DT∗ Dt∗ = 1 P e∇ 2 T∗.
4.2 Dimensionslöst med Ra och P r
På samma sätt som med Ri, Re och P e introduceras här dimensionslösa paramet-rar w∗ = w Uκ = wh κ , t ∗ = Uκ h t = κ h2t, p ∗ = p U2 κρ0 = ph 2 κ2ρ 0 , z∗ = z h, T ∗ = T − T0 ∆T .
I det betraktade fallet nns ingen yttre referenshastighet då uiden är stillstående till en början. Referenshastigheten Uκ de nieras därför som Uκ = κ/h och kan
ses som den hastighet som orsakas av temperaturgradienten. Med ovanstående parametrar insatta i (9) fås D(κhw∗) D(hκ2t∗) = − 1 ρ0 ρ(κ2ρ0 h2 p ∗) ρ(hz∗) + αδT T ∗ g + ν∂ 2(κ hw ∗) ∂(hx∗ i)2 .
Genom multiplikation med h3/νκfås κ ν Dw∗ Dt∗ = − κ ν ρp∗ ρz∗ + h3αδT g κν T ∗ + ∇2w∗. Efter identi ering av de dimensionslösa talen Ra och P r fås
1 P r Dw∗ Dt∗ = − 1 P r ∂p∗ ∂z∗ + RaT ∗ + ∇2w∗. (12)
På samma sätt skrivs (10) om vilket ger κ∆T h2 ∂T∗ ∂t∗ + κ∆T h2 u ∗ i ∂T∗ ∂x∗i = κ∆T h2 ∂2T∗ ∂(x∗i)2.
Omskrivning med matriell derivata ger DT∗
Dt∗ = ∇ 2T∗
. De aktuella ekvationerna blir alltså
∇ · uuu, 1 P r Dw∗ Dt∗ = − 1 P r ∂p∗ ∂z∗ + RaT ∗ + ∇2w∗, DT∗ Dt∗ = ∇ 2T∗ .
Eftersom ekvation (11) och (12) bara är olika sätt att skriva samma ekvation ger jämförelse av dessa att
Re = 1 P r.
5 Instabilitet
in hans experiment.
Antag ett experiment, likt det som Bénard utförde, med två oändligt långa plattor med en uid emellan, där den undre plattan har en högre temperatur än den övre. Detta leder till att densiteten vid den undre plattan blir lägre än vid den övre. En ytkraft uppstår som vill driva en rörelse i uiden. Denna motverkas dock av viskositet och värmediusion i uiden.
Skillnaden mellan detta experiment och det utfört av Bénard är i synnerhet den fria övre ytan. Med en fri yta uppkommer två eekter som även de påverkar stabilite-ten. Marangonieekten som kommer från den temperaturberoende ytspänningen, samt Bioteekten som beskriver värme ödet över den fria ytan. Det har även se-nare bevisats att desto tunse-nare lager av uid som används, desto större inverkan har Marangonieekten. Vilket leder till, att i Bénards fall, var det denna som var den största anledningen till resultatet. För vidare läsning om detta rekommende-ras boken Convection in uids (R. Kh. Zeytounian 2009).
En teoretiskt förklaring av vad som inverkar formulerades av Rayleigh år 1916 (R. Kh. Zeytounian 2009, kap. 3) i form av det så kallade Rayleightalet,
Ra = gαΓd
4
κν ,
där Γ = −d ˜T /dz är den vertikala temperaturgradienten i grundstillståndet. Tem-peraturgradienten måste vara tillräckligt stor för att kvoten ovan ska uppnå ett speci kt värde som svarar mot instabilitet. Att detta tal är direkt kopplat till frå-geställningen är tydligt ty Rayleightalet beskriver förhållandet mellan ytkraften och viskösa krafterna (för bevis se bilaga A.4).
5.1 Teori
tankegång och struktur som Fluid Mechanics (Pijush, Kundu & Cohen, 2008, kap. 5).
5.1.1 Små störningars utveckling
I detta kapitel används indexnotation där (x1, x2, x3) = (x, y, z)är ortsvektorn och
(u1, u2, u3) = (u, v, w) är hastighetsvektorn. Genom att utgå från de Boussinesqa
ekvationerna (8) - (10) med omskrivningen ˜
ρ = ρ0[1 − α( ˜T − T0)],
där ρ0 representerar en referensdensitet vid en referenstemperatur T0 fås
∂ ˜ui ∂t + ˜uj ∂ ˜ui ∂xj = −1 ρ0 ∂ ˜p ∂xi − g[1 − α( ˜T − T0)]δi3+ ν∇2u˜i, (13) ∂ ˜T ∂t + ˜uj ∂ ˜T ∂xj = κ∇2T ,˜ (14) ∂ ˜ui ∂xi = 0.
I ekvationerna ovan denoterar (˜u) superpositionen av ett grundtillstånd och en störning. Dessa skrivs som
˜ ui = 0 + ui(xxx, t), ˜ T = ¯T (x3) + T0(xxx, t), ˜ p = P (x3) + p(xxx, t).
Här är alltså vänsterleden de totala tillstånden, versaler betyder grundtillstånd och gemener, samt T0, är störningarna. Här är U = 0 ty uiden är stillastående i
grundtillståndet. Insatt i ekvation (13) och ekvation (14) fås ∂ui ∂t + uj ∂ui ∂xj = −1 ρ0 ∂ ∂xi (P + p) − g[1 − α( ¯T + T0− T0)]δi3+ ν∇2ui, ∂T0 ∂t + uj ∂ ∂xj ( ¯T + T0) = κ∇2( ¯T + T0). (15)
Genom att subtrahera grundtillståndet (16) från det totala tillståndet (15) och försumma de icke-linjära termerna för störningarna störningarna fås
∂ui ∂t = − 1 ρ0 ∂p ∂xi + gαT0δi3+ ν∇2ui, ∂T0 ∂t + u3 d ¯T dx3 = κ∇2T0. (17)
I denna förenkling har även använts att ¯T = ¯T (x3). För att få uttrycket för ¯T så
integreras energiekvationen i (16) vilket ger ¯
T = C1x3+ C2. (18)
Integrationskonstanterna C1 och C2 bestäms av randvillkoren som fås med hjälp
av gur 2 nedan.
Figur 2
Eftersom plattorna hålls vid konstant temperatur blir randvillkoren för tempera-turen
¯
T (x3 = 0) = T0+ ∆T,
¯
T (x3 = d) = T0.
Detta ger konstanterna
C2 = T0+ ∆T,
C1 = −
∆T d , vilka insatta i ekvation (18) ger temperaturen som
Vid insättning av detta i energiekvationen (17) fås ∂T0
∂t − Γu3 = κ∇
2T0, (19)
där Γ = ∆T/d. Γ är alltså den linjära vertikala temperaturgradienten. För att skriva om ekvationerna i termer av u3 och T0 görs följande matematiska
manipu-lationer. Först utvärderas Laplacianen av x3-komponenten i den första ekvationen
i (17). ∂ ∂t ∇ 2 u3 = − 1 ρ0 ∇2 ∂p ∂x3 + gα∇2T0+ ν∇4u3 (20)
Ekvation (20) kan förenklas ytterliggare. Betrakta divergensen av första ekvationen i (17). ∂ ∂t ∂ui ∂xi = −1 ρ0 ∂2p ∂xi∂xi + gα∂T 0 ∂xi δi3+ ν ∂2 ∂xj∂xj ∂ui ∂xi
Genom att derivera denna med avseende på x3 och eftersom ∂ui/∂xi = 0 fås
0 = −1 ρ0 ∇2 ∂p ∂x3 + gα∂ 2T0 ∂x2 3 . (21)
Subtraheras (21) från (20) fås de ekvationer som beskriver störningarna utveckling som ∂ ∂t ∇ 2u 3 = gα ∂2 ∂x2 1 + ∂ 2 ∂x2 2 T0+ ν∇4u3, ∂T0 ∂t − Γu3 = κ∇ 2 T0. (22)
Här görs de oberoende variablerna x3 och t i (22) dimensionslösa på följande sätt.
t∗ = κ d2t, x ∗ i = xi d. Dessa insatta i (22) ger
κ d4 ∂ ∂t∗ ∇ 2u 3 = gα d2 ∂2 ∂(x∗ 1)2 + ∂ 2 ∂(x∗ 2)2 T0+ ν d4∇ 4u 3, κ d2 ∂T0 ∂t∗ − Γu3 = κ∇ 2T0 .
där P r = ν/κ. Randvillkoren formuleras som u3 = ∂u3 ∂x∗3 = T 0 = 0, (24) då x∗
3 = 0 och då x∗3 = 1 eftersom plattorna hålls vid konstant temperatur och
vidhäftningsvillkoret1 gäller överallt vid plattorna.
5.1.2 Egensvängningsmetoden
Störningarna i hastighet och temperatur, u3 och T0, skrivs på formen
u3 = ˆu3(x∗3)e ikx∗ 1+ilx∗2+σt∗, T0 = ˆT (x∗3)eikx∗1+ilx ∗ 2+σt ∗ , (25)
där ˆT och ˆu3 är komplexa. Då strömningsfältet som undersöks är obegränsat i x1
-och x2-led måste k och l vara reella tal, det vill säga störningarna är periodiska i
x1- och x2-led. Huruvida störningen är växande eller avtagande beror således på
σ, vilken består av en realdel och en imaginärdel; σ = σr+ iσi.
Huruvida strömningen är stabil eller inte kan ses genom följande: σr < 0 : stabilt tillstånd
σr > 0 : instabilt tillstånd
Utöver dessa två lägen tillkommer givetvis fallet då σr = 0, detta läge kallas
mar-ginaltillståndet. Marginaltillståndet är gränsen mellan stabilitet och instabilitet. Det nns två former av marginaltillstånd; då σi = 0 respektive då σi 6= 0.
Skill-naden mellan dessa två är att om σi 6= 0 fås en oscillerande störning med ökande
amplitud. Om istället σi = 0 fås ett stationärt öde. I detta fallet, ett cellulärt
mönster. I det betraktade fallet gäller att σi = 0. För bevis av detta se Fluid
Me-chanics (Pijush, Kundu & Cohen, 2008, s.475-477). Med störningarna de nierade enligt (25) motsvaras operatorerna i (23) av nya operatorer enligt
∂ ∂t∗ → σ, ∂2 ∂(x∗1)2 + ∂2 ∂(x∗2)2 → −K 2, ∇2 → d2 d(x∗3)2 − K 2
där K = k2+ l2. Detta ger att ekvationerna i (23) kan skrivas som σ − D2− K2 ˆ T = Γd 2 κ uˆ3, h σ P r − D 2− K2i D2− K2 ˆu3 = − gαd2K2 ν ˆ T , där D = d/dx∗
3. Med (Γd2/κ) ˆu3 = W samt gαΓd4/ (κν) = Ra fås att
σ − D2− K2 ˆ T = W, (26) h σ P r − D 2− K2i D2− K2 W = −RaK2ˆ T . (27)
Detta ger att randvillkoren (24) skrivs som
W = DW = ˆT = 0, (28)
då x∗
3 = 0 och då x∗3 = 1.
5.1.3 Gränsen till instabilitet
För att hitta marginaltillståndet i det betraktade fallet med två oöndligt stora plattor används ekvationerna (26) och (27). Vid marginaltillståndet gäller att σ = 0 och ˆT kan då elimineras genom att lösa ut ˆT ur (26) och sätta in i (27). Då fås
(D2 − K2)3W = −RaK2W, (29) med randvillkoren W = DW = (D2− K2)2W = 0, (30) då x∗ 3 = 0och då x ∗
3 = 1. För att lösningen till problemet ska vara nollskild krävs
ett speci kt förhållande mellan Ra och K. Det är tydligt att ekvation (29) är en sjätteordningens dierentialekvation och därmed genererar sex rötter.
Problemets utformning kan skilja sig. Antingen kan problemets vertikalhastighet vara symmetrisk kring x1x2-planet vid x3 = 12, eller också kan den vara assymetrisk
kring detta plan. Det förstnämnda kallas udda mod, det andra kallas jämn mod. Jämn mod ger en strömning med en rad av konvektionsceller, medans udda mod ger två rader (Pijush, Kundu & Cohen, 2008, s.477).
Då ekvation (29) inte är beroende av x∗
3 kan lösningen ansättas som
vilket ger (q2− K2)3 = RaK2. Rötterna till detta tredjegradspolynom är q2 = −K2 " Ra K4 1/3 − 1 # , q2 = K2 " 1 + 1 2 Ra K4 1/3 1 ± i√3 # .
Vilket ger rötterna ±iq0, ±q samt ±q0, där q0 är konjugatet till q. q0 är
q0 = K[(
Ra K4)
1/3− 1]1/2.
För det jämna läget blir därför lösningen till ekvation (29) W = A cosh q0x∗3+ B cosh qx
∗
3+ C cosh q 0
x∗3,
där A,B och C är konstanter. Genom att derivera ovanstånde lösning kan rand-villkoren (30) användas vilket leder till ett ekvationssystem på formen
cosq0 2 cosh q 2 cosh q0 2
−q0sinq20 q sinhq2 q0sinhq
0 2 (q02+ K2)2cosq0 2 (q 2− K2)2coshq 2 q 02− K22 coshq20 A B C = 0. (32) För icke-trivial lösning krävs att A, B och C inte alla får vara noll på samma gång, vilket innebär att determinanten måste vara noll. Med Matlab fås de K och Ra som uppfyller ovanstående ekvationssystem.
6 Simulering med Simson
För att kontrollera huruvida den förenklade teoretiska ansatsen stämmer används programmet Simson (Chevalier et al., 2007). Det är en DNS-kod som utvecklats för att lösa inkompressibla gränsskiktsströmningsfall. DNS (Direct Numerical Sol-ver) innebär en diskretisering av både tid- och rumsdomänen varvid Navier-Stokes ekvationer löses numeriskt (Moin, Mahesh, 1998).
En rad simuleringar för olika vågtal K och Ra-tal körs. Resultaten som sparas är huruvida strömningsfältet är stabilt eller instabilt. Resultaten presenteras sedan tillsammans i form av ett stabilitetdiagram.
Utöver att bestämma stabilitet/instabilitet möjliggör Simson att det faktiska strömningsfältet kan visualiseras i form av temperaturfördelning och hastighets-fält.
7 Resultat
7.1 Resultat från teorin
Ekvationssytemet (32) ger förhållandet mellan Ra och K för marginaltillståndent. Detta presenteras i ett stabilitetsdiagram ( gur 3).
Stabilt (σ < 0) Instabilt (σ > 0) ← Marginaltillstånd (σ = 0) ← (1708, 3.12) Ra K 1000 2000 3000 4000 5000 6000 7000 8000 1 2 3 4 5 6 7 8
Figur 3: Teoretiskt stabilitetsdiagram.
7.2 Resultat från simulering
Utifrån beräkningarna i Simson plottas ännu ett stabilitetsdiagram ( gur 4).
0 1000 2000 3000 4000 5000 6000 7000 8000 0 1 2 3 4 5 6 7 8 Stabilt (σ < 0) Instabilt (σ > 0) Ra K Instabila punkter Stabila punkter
Figur 4: Stabilitetsdiagram från simulering.
0 50 100 150 200 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 Tid Hastighet i y−led Ra = 1700, Stabilt Ra = 1730, Instabilt (a) K = 3.14. 0 50 100 150 200 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 Tid Hastighet i y−led K=6.51, Stabilt K=6.41, Instabilt (b) Ra = 4000.
Figur 6 visar strömningen i xy-planet. 0 0.5 1 1.5 2 2.5 3 3.5 4 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 x y
Temperaturfördelningen i samma plan ritas upp för ett instabilt och ett stabilt fall ( gurer 7 och 8). 0 0.5 1 1.5 2 2.5 3 3.5 4 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 x y
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 x y
7.3 Jämförelse
Det simulerade stabilitetsdiagrammet läggs över det teoretiska stabilitetsdiagram-met för att se hur väl det båda resultaten stämmer överens.
Stabilt (σ < 0) Instabilt (σ > 0) Ra K 0 1000 2000 3000 4000 5000 6000 7000 8000 0 1 2 3 4 5 6 7 8 Teoretiska marginaltillståndet Instabila punkter Stabila punkter
Figur 9: En jämförelse av det teoretiska stabilitetsdiagrammet och simuleringar.
8 Diskussion
I detta arbete har vi undersökt huruvida dessa lösningar med förenklingar och antaganden stämmer överens med verkligheten genom att simulera systemet i en DNS-kod med programmet Simson. Som tidigare nämnt löser denna kod Navier-Stokes ekvationer numeriskt genom diskretisering av tid- och rumsdomänen. Ko-den ger verklighetstrogna lösningar och har i Ko-denna rapport använts som substitut för experiment. Det bör dock tilläggas att, som visas i härledningen av Navier-Stokes ekvationer, det nns antaganden även här. Som till exempel att uiden är en Newtonsk uid. Dessa antaganden är viktiga att ta hänsyn till vid användandet av Navier-Stokes ekvationer. I det betraktade fallet används dessa antaganden bå-de i teorin och i simuleringen då bägge utgår från Navier-Stokes ekvationer. Även Boussinesq approximation används i de båda lösningarna.
Den stora skillnaden mellan simuleringen och teorin är att i teorin linjäriseras problemet. Trots att teorin är linjäriserad och simuleringen är icke-linjär överen-stämmer mätpunkterna i simuleringen väl med den teoretiska kurvan, se gur 9 i kapitel 7.3. Detta tyder på att de antaganden och förenklingar som gjorts har varit rimliga och att de metoder som används i teorin är väl överensstämmande med verkligheten för detta fall.
9 Slutsats
Följande slutsatser kan dras ur denna rapport:
Kritiska Rayleightalet fås till Ra=1708, med K= 3.12. Under kritiska Ray-leightalet är uiden stabil för alla vågtal K.
Referenser
Andersson K., Artman K., Astell M., Axberg S., Liwång H., Lundberg A., Norsell M. & Tornérhielm L. (2007) Lärobok i Militärteknik, Vol. 1: Grunder, Stockholm: Försvarshögskolan.
Mattias Chevalier, Philipp Schlatter, Anders Lundbladh, Dan S. Henningson. (2007) SIMSON A Pseudo-Spectral Solver for Incompressible Boundary Layer Flows. Stockholm: KTH Mechanics.
De nitions.net. 19 april 2012. http://www.de nitions.net/de nition/stability Parviz Moin, Krishnan Mahesh. (1998) DIRECT NUMERICAL SIMULA-TION: A Tool in Turbulence Research. Annual Review of Fluid Mechanics, Vol. 30, s.539-578.
Pijush K. Kundu, Ira M. Cohen. (2008) Fluid Mechanics. 4th ed. Amsterdam: Academic Press.
Spiegel E.A., Veronis G. (1960) On the Boussinesq Approximation for a Compressible Fluid. Astrophysical Journal, Vol. 131. Massachusetts.
A Bilagor
Vid följande härledningar används följande beteckningar. ρ =densitet, p =tryck, u uu = uidens hastighetsvektor, t =tid, V = x kontrollvolym, V =materiell kontrollvolym, n n n =utåtriktad enhetsnormal, S =volymens mantelarea.
A.1 Kontinuitetsekvationen
Kontinuitets ekvationen kan härledas genom att betrakta massbalansen för en x kontrollvolym placerad i en strömmande uid. En x kontrollvolym är en tänkt godtycklig volym som ej deformeras och är xerad i rummet. Massökningen i denna volym måste alltså vara lika med den inströmmande uidens massa,
d dt Z V ρdV = − Z S ρuuu · nnndS. (33) Vänsterledet i ekvation (33) kan skrivas om som
d dt Z V ρdV = Z V ∂ρ ∂tdV. (34)
Högerledet i ekvation (33) skrivs om med hjälp av Gauss sats, − Z S ρuuu · nnndS = − Z V ∇ · (ρuuu) dV. (35)
Med omskrivningarna (34) och (35) insatta i ekvation (33) och genom samling av
termer fås Z V ∂ρ ∂t + ∇ · (ρuuu) dV = 0.
Eftersom denna ska gälla för en godtycklig x kontrollvolym V så måste integran-den vara noll vilket ger följande ekvation som kallas kontinuitetsekvationen,
∂ρ
Ekvationen kan skrivas om med hjälp av materiell derivata och blir då Dρ
Dt + ρ∇ · uuu = 0.
A.2 Rörelsemängdsekvationen
A.2.1 Rörelsemängdens bevarande
För ett uidelement gäller att summan av yt- och volymskrafter är lika med pro-dukten av dess massa och acceleration. Ytkrafter i i-riktning är
∂τji
∂xj
dV.
Spänningstensorn τττ är symmetrisk ty uiden är isotrop, vilket innebär att τji = τij.
Newtons andra lag för ett uid element ger Cauchys rörelseekvation ρDui Dt = ρgi + δτij δxj , (36) där ρggg är volymskraft.
A.2.2 Konstitutiv ekvation för Newtonsk uid
På en uid i vila verkar enbart spänningar i dess ytors normalriktningar. Spän-ningen är oberoende av ytornas riktning. Spänningstensorn är därför isotrop och måste vara proportionell mot Kroneckers delta,
δδδ = 1 0 0 0 1 0 0 0 1 . Spänningstensorn kan därför skrivas som
τij = −pδij,
Hos en uid i rörelse nns även skjuvspänningar σij. Den dynamiska
spännings-tensorn blir
τij = −pδij + σij. (37)
σσσ är anisotrop och funktion av hastighetsgradienten δui
δxj som kan delas upp i en
symmetrisk och en antisymmetrisk del, ∂ui ∂xj = 1 2 ∂ui ∂xj + ∂uj ∂xi +1 2 ∂ui ∂xj − ∂uj ∂xi .
Den antisymmetriska delen representerar uidens rotation utan deformation och ger därför inte upphov till någon spänning. Spänningen måste därför helt och hållet komma från den symmetriska delen hos hastighetsgradienten, vilken kan kallas töjningstensorn, eij ≡ 1 2 ∂ui ∂xj +∂uj ∂xi . Spänningen antas bero linjärt av töjningen,
σij = Kijmnemn. (38)
Då uiden är isotrop och spänningstensorn isotrop måste även Kijmn vara isotrop
på formen
Kijmn = λδijδmn+ µδimδjn+ γδinδjm, (39)
där λ, µ och γ är konstanter som beror på uidens termodynamiska tillstånd. Då σij är symmetrisk i ij måste också Kijmn vara symmetrisk i ij vilket innebär att
µ = γ,
och med denna likhet insatt i (39) som sen sätts in i (38) fås σij = 2µeij+ λemmδij,
där emm = ∇ · uuu är töjning per volymsenhet.
Spänningstensorn (37) blir således
Diagonaltermerna i eij är inte nödvändigtvis lika. I sådant fall kan även τij ha
olika diagonaltermer. Därför tas medelvärdet av diagonaltermerna hos τττ som me-deltrycket ¯ p ≡ −1 3τii. (43) (43) insatt i (42) ger p − ¯p = 2 3µ + λ ∇ · uuu.
Skillnaden mellan termodynamiskt tryck p och medeltryck ¯p är proportionellt mot κ = λ +2
3µ.
Observera att detta κ kallas volymsviskositetskoe cienten och bör ej förväxlas med värmediusiviteten som betecknas med κ i övriga delar av rapporten. För att κ ska uppnå en mätbar storhet krävs att DρDt är extremt stor så som i stötvågor. Inom många tillämpningar, Rayleigh-Bènard konvektion inräknad, gäller Stokes antagande att
κ = λ +2
3µ = 0. (44)
Med Stokes antagande (44) reduceras (40) till τij = − p + 2 3µ∇ · uuu δij + 2µeij. (45)
A.2.3 Navier-Stokes ekvation
Rörelseekvationen för en Newtonsk uid fås då (45) sätts in i Cauchys ekvation (36), ρDui Dt = − ∂p ∂xi + ρgi+ ∂ ∂xj 2µeij − 2µ 3 (∇ · uuu) δij .
Viskositeten µ kan tas ut ur derivatatermen om temperaturskillnaden i uiden är tillräckligt liten, vilket ger
ρDui Dt = − ∂p ∂xi + ρgi+ 2µ ∂eij ∂xj −2µ 3 ∂ ∂xi (∇ · uuu) ⇔ ρDui Dt = − ∂p ∂xi + ρgi+ µ ∇2u i+ 1 3 ∂ ∂xi (∇ · uuu) , (46) där ∇2u i = ∂ 2u i ∂x2 j.
För en inkompressibel uid där ∇ · uuu = 0 reduceras (46) till ρDuuu
Dt = −∇p + ρggg + µ∇
2
A.3 Energiekvationen
Termodynamikens första lag applicerad på en materiell kontrollvolym, V, ger D Dt Z V ρ e +u 2 i 2 dV = Z V ρgiuidV + Z A τijuidAj− Z A qidAi.
där e är inre energi per volymsenhet, qi är värme öde per areaenhet.
Ändringen per tidsenhet av lagrad energi i den materiella kontrollvolymen är lika med ödet av arbete och värme genom dess ränder. Genom att skriva alla termer som volymsintegraler enligt Gauss sats fås
Z V ρD Dt e +u 2 i 2 dV = Z V ρgiui+ ∂ ∂xj (τijui) − ∂qi ∂xi dV ⇔ ρD Dt e + u 2 i 2 = ρgiui+ ∂ ∂xj (τijui) − ∂qi ∂xi . Genom att subtrahera mekaniska energins bidrag
ρD Dt u2 i 2 = ρggguuu + ∂ ∂xj (τijui) + p∇ · uuu − φ,
där φ är energidissipation från viskösa krafter, från ekvation (A.3) fås den termiska energiekvationen
ρDe
Dt = −∇ · qqq − p∇ · uuu + φ.
A.4 Rayleightalet
Rayleightalet beskriver kvoten av ytkraften med den viskösa kraften. Med andra ord hur dessa två krafters storlek förhåller sig till varandra. Detta kan inses genom att jämföra storleken av advektions- och diusionstermen i ekvation (19),
u3 ∼ κT0/d2 Γ ∼ κΓ/d Γ = κ d.
